30#include <torch/script.h>
31#include <torch/torch.h>
35#include "RBFMotionSolver.H"
42using namespace ITHACAtorch::torch2Eigen;
53 NUmodes = para->ITHACAdict->lookupOrDefault<label>(
"NmodesUproj", 10);
54 NNutModes = para->ITHACAdict->lookupOrDefault<label>(
"NmodesNutProj", 10);
55 Net->push_back(torch::nn::Linear(NUmodes, 128));
56 Net->push_back(torch::nn::ReLU());
57 Net->push_back(torch::nn::Linear(128, 64));
58 Net->push_back(torch::nn::ReLU());
59 Net->push_back(torch::nn::Linear(64, NNutModes));
60 optimizer =
new torch::optim::Adam(Net->parameters(),
61 torch::optim::AdamOptions(2e-2));
71 Eigen::MatrixXd bias_inp;
72 Eigen::MatrixXd scale_inp;
73 Eigen::MatrixXd bias_out;
74 Eigen::MatrixXd scale_out;
76 Eigen::MatrixXd coeffL2Nut;
77 Eigen::MatrixXd coeffL2U;
79 torch::Tensor coeffL2U_tensor;
80 torch::Tensor coeffL2Nut_tensor;
82 torch::nn::Sequential Net;
83 torch::optim::Optimizer* optimizer;
84 torch::jit::script::Module netTorchscript;
86 void loadNet(word filename)
88 std::string Msg = filename +
89 " is not existing, please run the training stage of the net with the correct number of modes for U and Nut";
91 netTorchscript = torch::jit::load(filename);
92 cnpy::load(bias_inp,
"ITHACAoutput/NN/minAnglesInp_" + name(
94 cnpy::load(scale_inp,
"ITHACAoutput/NN/scaleAnglesInp_" + name(
96 cnpy::load(bias_out,
"ITHACAoutput/NN/minOut_" + name(
NUmodes) +
"_" + name(
98 cnpy::load(scale_out,
"ITHACAoutput/NN/scaleOut_" + name(
NUmodes) +
"_" + name(
100 netTorchscript.eval();
108 mkDir(
"ITHACAoutput/NN/coeffs");
110 PtrList<volVectorField> UfieldTrain;
111 PtrList<volScalarField> PfieldTrain;
112 PtrList<volScalarField> nutFieldsTrain;
114 "./ITHACAoutput/Offline/");
116 "./ITHACAoutput/Offline/");
117 auto nut_train =
_mesh().lookupObject<volScalarField>(
"nut");
119 "./ITHACAoutput/Offline/");
121 Info <<
"Computing the coefficients for U train" << endl;
125 Info <<
"Computing the coefficients for p train" << endl;
129 Info <<
"Computing the coefficients for nuT train" << endl;
133 coeffL2U_train.transposeInPlace();
134 coeffL2P_train.transposeInPlace();
135 coeffL2Nut_train.transposeInPlace();
136 cnpy::save(coeffL2U_train,
"ITHACAoutput/NN/coeffs/coeffL2UTrain.npy");
137 cnpy::save(coeffL2P_train,
"ITHACAoutput/NN/coeffs/coeffL2PTrain.npy");
138 cnpy::save(coeffL2Nut_train,
"ITHACAoutput/NN/coeffs/coeffL2NutTrain.npy");
140 PtrList<volVectorField> UfieldTest;
141 PtrList<volScalarField> PfieldTest;
142 PtrList<volScalarField> nutFieldsTest;
145 "./ITHACAoutput/checkOff/");
147 "./ITHACAoutput/checkOff/");
148 auto nut_test =
_mesh().lookupObject<volScalarField>(
"nut");
150 "./ITHACAoutput/checkOff/");
161 coeffL2U_test.transposeInPlace();
162 coeffL2P_test.transposeInPlace();
163 coeffL2Nut_test.transposeInPlace();
164 cnpy::save(coeffL2U_test,
"ITHACAoutput/NN/coeffs/coeffL2UTest.npy");
165 cnpy::save(coeffL2P_test,
"ITHACAoutput/NN/coeffs/coeffL2PTest.npy");
166 cnpy::save(coeffL2Nut_test,
"ITHACAoutput/NN/coeffs/coeffL2NutTest.npy");
170 Eigen::MatrixXd evalNet(Eigen::MatrixXd a, Eigen::MatrixXd mu_now)
172 Eigen::MatrixXd xpred(a.rows() + mu_now.rows(), 1);
174 if (xpred.cols() != 1)
176 throw std::runtime_error(
"Predictions for more than one sample not supported yet.");
179 xpred.bottomRows(a.rows()) = a;
180 xpred.topRows(mu_now.rows()) = mu_now;
181 xpred = xpred.array() * scale_inp.array() + bias_inp.array() ;
182 xpred.transposeInPlace();
183 torch::Tensor xTensor = eigenMatrix2torchTensor(xpred);
185 std::vector<torch::jit::IValue> XTensorInp;
186 XTensorInp.push_back(xTensor);
187 out = netTorchscript.forward(XTensorInp).toTensor();
188 Eigen::MatrixXd g = torchTensor2eigenMatrix<double>(out);
189 g.transposeInPlace();
190 g = (g.array() - bias_out.array()) / scale_out.array();
195 Eigen::MatrixXd evalNet(Eigen::MatrixXd a)
197 netTorchscript.eval();
198 a.transposeInPlace();
199 a = (a.array() - bias_inp.array()) / scale_inp.array();
200 torch::Tensor aTensor = eigenMatrix2torchTensor(a);
202 std::vector<torch::jit::IValue> XTensorInp;
203 XTensorInp.push_back(aTensor);
204 out = netTorchscript.forward(XTensorInp).toTensor();
205 Eigen::MatrixXd g = torchTensor2eigenMatrix<double>(out);
206 g = g.array() * scale_out.array() + bias_out.array();
207 g.transposeInPlace();
225 void projectReducedOperators(
int NmodesUproj,
int NmodesPproj,
int NmodesEproj)
227 PtrList<volVectorField> gradModP;
229 for (label i = 0; i < NmodesPproj; i++)
231 gradModP.append(fvc::grad(
problem->Pmodes[i]));
238 void solveOnlineCompressible(
int NmodesUproj,
int NmodesPproj,
int NmodesEproj,
239 int NmodesNutProj, Eigen::MatrixXd mu_now,
240 word Folder =
"./ITHACAoutput/Online/")
244 scalar residualNorm(1);
245 scalar residualJump(1);
246 Eigen::MatrixXd uResidualOld = Eigen::MatrixXd::Zero(1, NmodesUproj);
247 Eigen::MatrixXd eResidualOld = Eigen::MatrixXd::Zero(1, NmodesEproj);
248 Eigen::MatrixXd pResidualOld = Eigen::MatrixXd::Zero(1, NmodesPproj);
249 Eigen::VectorXd uResidual(Eigen::Map<Eigen::VectorXd>(uResidualOld.data(),
251 Eigen::VectorXd eResidual(Eigen::Map<Eigen::VectorXd>(eResidualOld.data(),
253 Eigen::VectorXd pResidual(Eigen::Map<Eigen::VectorXd>(pResidualOld.data(),
257 float residualJumpLim =
258 para->ITHACAdict->lookupOrDefault<
float>(
"residualJumpLim", 1e-5);
259 float normalizedResidualLim =
260 para->ITHACAdict->lookupOrDefault<
float>(
"normalizedResidualLim", 1e-5);
262 para->ITHACAdict->lookupOrDefault<
float>(
"maxIter", 2000);
263 bool closedVolume =
false;
266 fluidThermo& thermo =
problem->pThermo();
268 volScalarField& P =
problem->_p();
269 volScalarField& E =
problem->pThermo->he();
270 volScalarField& nut =
const_cast<volScalarField&
>
272 volScalarField& rho =
problem->_rho();
273 volScalarField& psi =
problem->_psi();
280 Eigen::MatrixXd u = Eigen::MatrixXd::Zero(NmodesUproj, 1);
281 Eigen::MatrixXd e = Eigen::MatrixXd::Zero(NmodesEproj, 1);
282 Eigen::MatrixXd p = Eigen::MatrixXd::Zero(NmodesPproj, 1);
284 NmodesNutProj,
true);
288 vector Uinlet(
problem->
_U().boundaryFieldRef()[idInl][0][0], 0, 0);
291 while ((residualJump > residualJumpLim
292 || residualNorm > normalizedResidualLim) && csolve < maxIter)
295 Info <<
"csolve:" << csolve << endl;
301 uResidualOld = uResidual;
302 eResidualOld = eResidual;
303 pResidualOld = pResidual;
305 List<Eigen::MatrixXd> RedLinSysU;
310 - fvc::div((rho *
problem->turbulence->nuEff()) * dev2(T(fvc::grad(U))))
311 - fvm::laplacian(rho *
problem->turbulence->nuEff(), U)
316 fvOptions.constrain(UEqnR);
320 RedLinSysU[1] = RedLinSysU[1] - projGradP;
326 fvOptions.correct(U);
331 + fvc::div(phi, volScalarField(
"Ekp", 0.5 * magSqr(U) + P / rho))
332 - fvm::laplacian(
problem->turbulence->alphaEff(), E)
337 fvOptions.constrain(EEqnR);
342 fvOptions.correct(E);
345 constrainPressure(P, rho, U,
problem->getPhiHbyA(UEqnR, U, P),
348 surfaceScalarField phiHbyACalculated =
problem->getPhiHbyA(UEqnR, U, P);
349 closedVolume = adjustPhi(phiHbyACalculated, U, P);
350 List<Eigen::MatrixXd> RedLinSysP;
354 volScalarField rAU(1.0 /
356 volVectorField HbyA(constrainHbyA(rAU * UEqnR.H(), U,
358 surfaceScalarField phiHbyA(
"phiHbyA", fvc::interpolate(rho)*fvc::flux(HbyA));
359 surfaceScalarField rhorAUf(
"rhorAUf", fvc::interpolate(rho * rAU));
363 - fvm::laplacian(rhorAUf, P)
365 fvOptions(psi, P, rho.name())
369 problem->_pressureControl().refCell(),
370 problem->_pressureControl().refValue()
378 phi =
problem->getPhiHbyA(UEqnR, U, P) + PEqnR.flux();
382#include "continuityErrs.H"
385 U.correctBoundaryConditions();
386 fvOptions.correct(U);
387 bool pLimited =
problem->_pressureControl().limit(P);
392 P += (
problem->_initialMass() - fvc::domainIntegrate(psi * P))
393 / fvc::domainIntegrate(psi);
396 if (pLimited || closedVolume)
398 P.correctBoundaryConditions();
403 Info <<
"Ures = " << (uResidual.cwiseAbs()).sum() /
404 (RedLinSysU[1].cwiseAbs()).sum() << endl;
405 Info <<
"Eres = " << (eResidual.cwiseAbs()).sum() /
406 (RedLinSysE[1].cwiseAbs()).sum() << endl;
407 Info <<
"Pres = " << (pResidual.cwiseAbs()).sum() /
408 (RedLinSysP[1].cwiseAbs()).sum() << endl;
409 residualNorm =
max(
max((uResidual.cwiseAbs()).sum() /
410 (RedLinSysU[1].cwiseAbs()).sum(),
411 (pResidual.cwiseAbs()).sum() / (RedLinSysP[1].cwiseAbs()).sum()),
412 (eResidual.cwiseAbs()).sum() / (RedLinSysE[1].cwiseAbs()).sum());
413 residualJump =
max(
max(((uResidual - uResidualOld).cwiseAbs()).sum() /
414 (RedLinSysU[1].cwiseAbs()).sum(),
415 ((pResidual - pResidualOld).cwiseAbs()).sum() /
416 (RedLinSysP[1].cwiseAbs()).sum()),
417 ((eResidual - eResidualOld).cwiseAbs()).sum() /
418 (RedLinSysE[1].cwiseAbs()).sum());
422 nutCoeff =
problem->evalNet(u, mu_now);
423 problem->nutModes.reconstruct(nut, nutCoeff,
"nut");
445 dyndict =
new IOdictionary
449 "dynamicMeshDictRBF",
459 ms =
new RBFMotionSolver(
_mesh(), *dyndict);
460 vectorField motion(ms->movingPoints().size(), vector::zero);
461 movingIDs = ms->movingIDs();
462 x0 = ms->movingPoints();
464 point0 = ms->curPoints();
466 middleExport = para->ITHACAdict->lookupOrDefault<
bool>(
"middleExport",
true);
473 IOdictionary* dyndict;
483 double f1(
double chord,
double x)
485 double res = chord * (std::pow((x) / chord,
486 0.5) * (1 - (x) / chord)) / (std::exp(15 * (x) / chord));
490 List<vector> moveBasis(
const List<vector>& originalPoints,
double par)
492 List<vector> movedPoints(originalPoints);
494 for (
int i = 0; i < originalPoints.size(); i++)
496 movedPoints[i][2] += par * f1(1, movedPoints[i][0]);
502 void updateMesh(
double parTop = 0,
double parBot = 0)
504 _mesh().movePoints(point0);
506 if (parTop != 0 || parBot != 0)
509 List<vector> top0_cur = moveBasis(top0, parTop);
510 List<vector> bot0_cur = moveBasis(bot0, parBot);
513 ms->setMotion(curX - x0);
514 point = ms->curPoints();
515 _mesh().movePoints(point);
523 volVectorField& U =
_U();
525 volScalarField& p =
_p();
527 volScalarField& E = _E();
536 auto nut =
_mesh().lookupObject<volScalarField>(
"nut");
543 double UIFinit = para->ITHACAdict->lookupOrDefault<
double>(
"UIFinit", 170);
544 Vector<double>
Uinl(UIFinit, 0, 0);
546 for (label i = 0; i <
mu.rows(); i++)
548 updateMesh(
mu(i, 0),
mu(i, 1));
553 word polyMesh2beLinked = folder + name(i + 1) +
"/" +
"polyMesh/";
557 word folderContLink = folder + name(i + 1) +
"/" + name(j) +
"/";
558 system(
"ln -s $(readlink -f " + polyMesh2beLinked +
") " + folderContLink +
567int main(
int argc,
char* argv[])
572 std::ifstream exFileOff(
"./parsOff_mat.txt");
580 int OffNum = para->ITHACAdict->lookupOrDefault<
int>(
"OffNum", 100);
581 double BumpAmp = para->ITHACAdict->lookupOrDefault<
double>(
"BumpAmp", 0.1);
582 example.mu.resize(OffNum, 2);
587 example.mu.leftCols(1) = parTop;
588 example.mu.rightCols(1) = parBot;
592 Eigen::MatrixXd parsOn;
593 std::ifstream exFileOn(
"./parsOn_mat.txt");
601 int OnNum = para->ITHACAdict->lookupOrDefault<
int>(
"OnNum", 20);
602 double BumpAmp = para->ITHACAdict->lookupOrDefault<
double>(
"BumpAmp", 0.1);
603 parsOn.resize(OnNum, 2);
606 parsOn.leftCols(1) = parTopOn;
607 parsOn.rightCols(1) = parBotOn;
612 int NmodesUout = para->ITHACAdict->lookupOrDefault<
int>(
"NmodesUout", 15);
613 int NmodesPout = para->ITHACAdict->lookupOrDefault<
int>(
"NmodesPout", 15);
614 int NmodesEout = para->ITHACAdict->lookupOrDefault<
int>(
"NmodesEout", 15);
615 int NmodesNutOut = para->ITHACAdict->lookupOrDefault<
int>(
"NmodesNutOut", 15);
616 int NmodesUproj = para->ITHACAdict->lookupOrDefault<
int>(
"NmodesUproj", 10);
617 int NmodesPproj = para->ITHACAdict->lookupOrDefault<
int>(
"NmodesPproj", 10);
618 int NmodesEproj = para->ITHACAdict->lookupOrDefault<
int>(
"NmodesEproj", 10);
619 int NmodesNutProj = para->ITHACAdict->lookupOrDefault<
int>(
"NmodesNutProj", 10);
620 example.updateMesh();
622 example.offlineSolve();
624 example.updateMesh();
641 example.updateMesh();
645 checkOff.offline =
false;
647 checkOff.offlineSolve(
"./ITHACAoutput/checkOff/");
648 checkOff.offline =
true;
668 example.loadNet(
"ITHACAoutput/NN/Net_" + name(example.NUmodes) +
"_" + name(
669 example.NNutModes) +
".pt");
674 for (label k = 0; k < parsOn.rows(); k++)
676 example.updateMesh();
677 example.updateMesh(parsOn(k, 0), parsOn(k, 1));
679 "./ITHACAoutput/Online_" + name(NmodesUproj) +
"_" + name(NmodesNutProj) +
"/",
680 name(k + 1) +
"/polyMesh/");
683 reduced.projectReducedOperators(NmodesUproj, NmodesPproj, NmodesEproj);
685 example.turbulence->validate();
686 Eigen::MatrixXd mu_now = parsOn.row(k);
687 mu_now.transposeInPlace();
688 reduced.solveOnlineCompressible(NmodesUproj, NmodesPproj, NmodesEproj,
689 NmodesNutProj, mu_now,
"./ITHACAoutput/Online_" + name(NmodesUproj) +
"_" +
690 name(NmodesNutProj) +
"/");
695 PtrList<volVectorField> UfieldCheck;
696 PtrList<volScalarField> PfieldCheck;
697 PtrList<volScalarField> EfieldCheck;
698 PtrList<volScalarField> nutFieldsCheck;
700 "./ITHACAoutput/checkOff/");
702 "./ITHACAoutput/checkOff/");
704 "./ITHACAoutput/checkOff/");
705 auto nutCheck = example._mesh().lookupObject<volScalarField>(
"nut");
707 "./ITHACAoutput/checkOff/");
713 Eigen::MatrixXd snapsCheck =
717 for (label k = 0; k < snapsCheck.rows(); k++)
719 fieldNum = fieldNum + snapsCheck(k, 0);
721 "./ITHACAoutput/checkOffSingle/");
723 "./ITHACAoutput/checkOffSingle/");
725 "./ITHACAoutput/checkOffSingle/");
727 "./ITHACAoutput/checkOffSingle/");
729 k + 1) +
"/polyMesh",
"./ITHACAoutput/checkOffSingle/" + name(k + 1) +
"/");
735 PtrList<volVectorField> onlineU;
736 PtrList<volScalarField> onlineP;
737 PtrList<volScalarField> onlineE;
738 PtrList<volScalarField> onlineNut;
739 PtrList<volVectorField> offlineU;
740 PtrList<volScalarField> offlineP;
741 PtrList<volScalarField> offlineE;
742 PtrList<volScalarField> offlineNut;
744 "./ITHACAoutput/Online_" + name(NmodesUproj) +
"_" + name(NmodesNutProj) +
"/");
746 "./ITHACAoutput/Online_" + name(NmodesUproj) +
"_" + name(NmodesNutProj) +
"/");
748 "./ITHACAoutput/Online_" + name(NmodesUproj) +
"_" + name(NmodesNutProj) +
"/");
749 auto nut = example._mesh().lookupObject<volScalarField>(
"nut");
751 "./ITHACAoutput/Online_" + name(NmodesUproj) +
"_" + name(NmodesNutProj) +
"/");
753 "./ITHACAoutput/checkOffSingle/");
755 "./ITHACAoutput/checkOffSingle/");
757 "./ITHACAoutput/checkOffSingle/");
769 "errorU" + name(NmodesUproj) +
"_" + name(NmodesNutProj),
"python",
770 "./ITHACAoutput/ErrorFields_" + name(NmodesUproj) +
"_" + name(
771 NmodesNutProj) +
"/");
773 "errorP" + name(NmodesUproj) +
"_" + name(NmodesNutProj),
"python",
774 "./ITHACAoutput/ErrorFields_" + name(NmodesUproj) +
"_" + name(
775 NmodesNutProj) +
"/");
777 "errorE" + name(NmodesUproj) +
"_" + name(NmodesNutProj),
"python",
778 "./ITHACAoutput/ErrorFields_" + name(NmodesUproj) +
"_" + name(
779 NmodesNutProj) +
"/");
781 "errorNut" + name(NmodesUproj) +
"_" + name(NmodesNutProj),
"python",
782 "./ITHACAoutput/ErrorFields_" + name(NmodesUproj) +
"_" + name(
783 NmodesNutProj) +
"/");
785 for (label j = 0; j < parsOn.rows(); j++)
787 volVectorField Ue = offlineU[j] - onlineU[j];
792 volScalarField Pe = offlineP[j] - onlineP[j];
796 volScalarField Ee = offlineE[j] - onlineE[j];
799 volScalarField Nute = offlineNut[j] - onlineNut[j];
807 "./ITHACAoutput/ErrorFields_" + name(NmodesUproj) +
"_" + name(
808 NmodesNutProj) +
"/");
810 "./ITHACAoutput/ErrorFields_" + name(NmodesUproj) +
"_" + name(
811 NmodesNutProj) +
"/");
813 "./ITHACAoutput/ErrorFields_" + name(NmodesUproj) +
"_" + name(
814 NmodesNutProj) +
"/");
816 "./ITHACAoutput/ErrorFields_" + name(NmodesUproj) +
"_" + name(
817 NmodesNutProj) +
"/");
826 std::cerr <<
"checkOff folder is missing, error analysis cannot be performed."
Header file of the steadyNS class.
Header file of the Foam2Eigen class.
Header file of the ITHACAPOD class.
Header file of the ITHACAstream class, it contains the implementation of several methods for input ou...
Header file of the reducedSteadyNS class.
CompressibleSteadyNN(int argc, char *argv[])
Constructor.
PtrList< volScalarField > Efield
List of pointers used to store the energy solutions.
CompressibleSteadyNS()
Null constructor.
volScalarModes Emodes
List of pointers used to form the energy modes.
bool middleExport
Export also intermediate fields.
Class for the definition of some general parameters, the parameters must be defined from the file ITH...
static ITHACAparameters * getInstance()
Gets an instance of ITHACAparameters, to be used if the instance is already existing.
GeometricField< Type, PatchField, GeoMesh > reconstruct(GeometricField< Type, PatchField, GeoMesh > &inputField, Eigen::MatrixXd Coeff, word Name)
Function to reconstruct the solution starting from the coefficients, in this case the field is passed...
List< Eigen::MatrixXd > project(fvMatrix< Type > &Af, label numberOfModes=0, word projType="G")
A function that project an FvMatrix (OpenFoam linear System) on the modes.
CompressibleSteadyNN * problem
Full problem.
ReducedCompressibleSteadyNN(CompressibleSteadyNN &FOMproblem)
Constructor.
Eigen::MatrixXd projGradModP
Projected gradient of the pressure modes.
ReducedCompressibleSteadyNS()
Construct Null.
PtrList< volScalarField > nutFields
List of snapshots for the solution for eddy viscosity.
volScalarModes nutModes
List of POD modes for eddy viscosity.
autoPtr< fvMesh > _mesh
Mesh.
static Eigen::MatrixXd solveLinearSys(List< Eigen::MatrixXd > LinSys, Eigen::MatrixXd x, Eigen::VectorXd &residual, const Eigen::MatrixXd &bc=Eigen::MatrixXd::Zero(0, 0), const std::string solverType="fullPivLu")
Linear system solver for the online problem.
Eigen::MatrixXd mu_samples
Matrix of parameters to be used for PODI, where each row corresponds to a sample point....
void assignIF(T &s, G &value)
Assign internal field condition.
bool offline
Boolean variable, it is 1 if the Offline phase has already been computed, else 0.
Eigen::MatrixXd mu
Row matrix of parameters.
void truthSolve()
Perform a TruthSolve.
autoPtr< surfaceScalarField > _phi
Flux.
autoPtr< simpleControl > _simple
simpleControl
PtrList< volScalarField > Pfield
List of pointers used to form the pressure snapshots matrix.
autoPtr< fv::options > _fvOptions
fvOptions
autoPtr< Time > _runTime
Time.
volVectorModes Umodes
List of pointers used to form the velocity modes.
PtrList< volVectorField > Ufield
List of pointers used to form the velocity snapshots matrix.
autoPtr< fvMesh > _mesh
Mesh.
label NUmodes
Number of velocity modes used for the projection.
label NNutModes
Number of nut modes used for the projection.
autoPtr< volVectorField > Uinl
Initial dummy field with all Dirichlet boundary conditions.
autoPtr< volVectorField > _U
Velocity field.
volScalarModes Pmodes
List of pointers used to form the pressure modes.
scalar cumulativeContErr
continuity error
autoPtr< volScalarField > _p
Pressure field.
Class where the tutorial number 2 is implemented.
void offlineSolve(word folder="./ITHACAoutput/Offline/")
Perform an Offline solve.
tutorial02(int argc, char *argv[])
Constructor.
T max(Eigen::SparseMatrix< T > &mat, label &ind_row, label &ind_col)
Find the maximum of a sparse Matrix (Useful for DEIM).
void getModes(PtrList< GeometricField< Type, PatchField, GeoMesh > > &snapshots, PtrList< GeometricField< Type, PatchField, GeoMesh > > &modes, word fieldName, bool podex, bool supex, bool sup, label nmodes, bool correctBC)
Computes the bases or reads them for a field.
void writePoints(pointField points, fileName folder, fileName subfolder)
Write points of a mesh to a file.
void exportSolution(GeometricField< Type, PatchField, GeoMesh > &s, fileName subfolder, fileName folder, word fieldName)
Export a field to file in a certain folder and subfolder.
void readMiddleFields(PtrList< GeometricField< Type, PatchField, GeoMesh > > &Lfield, GeometricField< Type, PatchField, GeoMesh > &field, fileName casename)
Funtion to read a list of volVectorField from name of the field including all the intermediate snapsh...
void exportMatrix(Eigen::Matrix< T, -1, dim > &matrix, word Name, word type, word folder)
Export the reduced matrices in numpy (type=python), matlab (type=matlab) and txt (type=eigen) format ...
List< Eigen::MatrixXd > readMatrix(word folder, word mat_name)
Read a three dimensional matrix from a txt file in Eigen format.
void read_fields(PtrList< GeometricField< Type, PatchField, GeoMesh > > &Lfield, word Name, fileName casename, int first_snap, int n_snap)
Function to read a list of fields from the name of the field and casename.
Eigen::VectorXd getCoeffs(GeometricField< Type, PatchField, GeoMesh > &snapshot, PtrList< GeometricField< Type, PatchField, GeoMesh > > &modes, label Nmodes, bool consider_volumes)
Projects a snapshot on a basis function and gets the coefficients of the projection.
double errorL2Rel(GeometricField< T, fvPatchField, volMesh > &field1, GeometricField< T, fvPatchField, volMesh > &field2, List< label > *labels)
Computes the relative error between two geometric Fields in L2 norm.
void setIndices2Value(labelList &ind2set, List< Type > &value2set, labelList &movingIDS, List< Type > &originalList)
Sets some given Indices of a list of objects to given values.
Eigen::MatrixXd rand(label rows, label cols, double min, double max)
Generates random matrix with random values in an interval.
double L2Norm(const T v)
Compute the L2 norm of v.
void assignBC(GeometricField< scalar, fvPatchField, volMesh > &s, label BC_ind, double value)
Assign uniform Boundary Condition to a volScalarField.
bool check_folder(word folder)
Checks if a folder exists.
bool check_file(std::string fileName)
Function that returns true if a file exists.
void createSymLink(word folder)
Creates symbolic links to 0, system and constant.
void getPointsFromPatch(fvMesh &mesh, label ind, List< vector > &points, labelList &indices)
Get the polabel coordinates and indices from patch.
Header file of the torch2Eigen namespace.