30#include <torch/script.h>
31#include <torch/torch.h>
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());
61 torch::optim::AdamOptions(2e-2));
82 torch::nn::Sequential
Net;
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";
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 std::cout <<
"Computing the coefficients for U train" << std::endl;
125 std::cout <<
"Computing the coefficients for p train" << std::endl;
129 std::cout <<
"Computing the coefficients for nuT train" << std::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;
182 xpred.transposeInPlace();
185 std::vector<torch::jit::IValue> XTensorInp;
186 XTensorInp.push_back(xTensor);
189 g.transposeInPlace();
198 a.transposeInPlace();
202 std::vector<torch::jit::IValue> XTensorInp;
203 XTensorInp.push_back(aTensor);
207 g.transposeInPlace();
227 PtrList<volVectorField> gradModP;
229 for (label
i = 0;
i < NmodesPproj;
i++)
231 gradModP.append(fvc::grad(
problem->Pmodes[
i]));
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);
268 volScalarField& P =
problem->_p();
269 volScalarField& E =
problem->pThermo->he();
270 volScalarField&
nut =
const_cast<volScalarField&
>
271 (
problem->_mesh().lookupObject<volScalarField>(
"nut"));
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);
287 problem->_mesh().boundaryMesh().findPatchID(
"inlet");
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)
317 RedLinSysU =
problem->Umodes.project(UEqnR,
320 RedLinSysU[1] = RedLinSysU[1] - projGradP;
323 problem->Umodes.reconstruct(
U, u,
"U");
331 + fvc::div(
phi, volScalarField(
"Ekp", 0.5 * magSqr(
U) + P /
rho))
332 - fvm::laplacian(
problem->turbulence->alphaEff(), E)
338 List<Eigen::MatrixXd> RedLinSysE =
problem->Emodes.project(EEqnR, NmodesEproj);
340 problem->Emodes.reconstruct(E, e,
"e");
348 surfaceScalarField phiHbyACalculated =
problem->getPhiHbyA(UEqnR,
U, P);
350 List<Eigen::MatrixXd> RedLinSysP;
352 while (
problem->_simple().correctNonOrthogonal())
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));
369 problem->_pressureControl().refCell(),
370 problem->_pressureControl().refValue()
372 RedLinSysP =
problem->Pmodes.project(PEqnR, NmodesPproj);
374 problem->Pmodes.reconstruct(P,
p,
"p");
376 if (
problem->_simple().finalNonOrthogonalIter())
378 phi =
problem->getPhiHbyA(UEqnR,
U, P) + PEqnR.flux();
385 U.correctBoundaryConditions();
392 P += (
problem->_initialMass() - fvc::domainIntegrate(
psi * P))
393 / fvc::domainIntegrate(
psi);
398 P.correctBoundaryConditions();
403 std::cout <<
"Ures = " << (uResidual.cwiseAbs()).sum() /
404 (RedLinSysU[1].cwiseAbs()).sum() << std::endl;
405 std::cout <<
"Eres = " << (eResidual.cwiseAbs()).sum() /
406 (RedLinSysE[1].cwiseAbs()).sum() << std::endl;
407 std::cout <<
"Pres = " << (pResidual.cwiseAbs()).sum() /
408 (RedLinSysP[1].cwiseAbs()).sum() << std::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");
449 "dynamicMeshDictRBF",
460 vectorField motion(
ms->movingPoints().size(), vector::zero);
462 x0 =
ms->movingPoints();
466 middleExport =
para->ITHACAdict->lookupOrDefault<
bool>(
"middleExport",
true);
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]);
506 if (parTop != 0 || parBot != 0)
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++)
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);
668 example.
loadNet(
"ITHACAoutput/NN/Net_" + name(example.
NUmodes) +
"_" + name(
674 for (label k = 0; k < parsOn.rows(); k++)
677 example.
updateMesh(parsOn(k, 0), parsOn(k, 1));
679 "./ITHACAoutput/Online_" + name(NmodesUproj) +
"_" + name(NmodesNutProj) +
"/",
680 name(k + 1) +
"/polyMesh/");
686 Eigen::MatrixXd mu_now = parsOn.row(k);
687 mu_now.transposeInPlace();
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."
int main(int argc, char *argv[])
Header file of the steadyNS class.
Header file of the Foam2Eigen class.
Header file of the ITHACAPOD class.
#define M_Assert(Expr, Msg)
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.
Eigen::MatrixXd scale_inp
Eigen::MatrixXd coeffL2Nut
torch::optim::Optimizer * optimizer
Eigen::MatrixXd evalNet(Eigen::MatrixXd a, Eigen::MatrixXd mu_now)
torch::nn::Sequential Net
torch::Tensor coeffL2U_tensor
torch::Tensor coeffL2Nut_tensor
Eigen::MatrixXd scale_out
torch::jit::script::Module netTorchscript
void loadNet(word filename)
Eigen::MatrixXd evalNet(Eigen::MatrixXd a)
PtrList< volScalarField > Efield
List of pointers used to store the energy solutions.
CompressibleSteadyNS()
Null constructor.
autoPtr< volScalarField > _p
volScalarModes Emodes
List of pointers used to form the energy modes.
bool middleExport
Export also intermediate fields.
autoPtr< volScalarField > _E
autoPtr< compressible::turbulenceModel > turbulence
void restart()
set all variables back to the values into the 0 folder
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.
void projectReducedOperators(int NmodesUproj, int NmodesPproj, int NmodesEproj)
CompressibleSteadyNN * problem
Full problem.
ReducedCompressibleSteadyNN(CompressibleSteadyNN &FOMproblem)
Constructor.
void solveOnlineCompressible(int NmodesUproj, int NmodesPproj, int NmodesEproj, int NmodesNutProj, Eigen::MatrixXd mu_now, word Folder="./ITHACAoutput/Online/")
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.
ITHACAparameters * para
parameters to be read from the ITHACAdict file
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.
bool podex
Boolean variable, it is 1 if the POD has already been computed, else 0.
void truthSolve()
Perform a TruthSolve.
PtrList< volScalarField > Pfield
List of pointers used to form the pressure snapshots matrix.
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.
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.
Class where the tutorial number 2 is implemented.
void offlineSolve(word folder="./ITHACAoutput/Offline/")
Perform an Offline solve.
void updateMesh(double parTop=0, double parBot=0)
tutorial02(int argc, char *argv[])
Constructor.
double f1(double chord, double x)
List< vector > moveBasis(const List< vector > &originalPoints, double par)
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.
torch::Tensor eigenMatrix2torchTensor(Eigen::Matrix< type, Eigen::Dynamic, Eigen::Dynamic > _eigenMatrix)
Convert an eigen Matrix to a torch tensor.
template Eigen::Matrix< double, Eigen::Dynamic, Eigen::Dynamic > torchTensor2eigenMatrix< double >(torch::Tensor &torchTensor)
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 L2Norm(GeometricField< scalar, fvPatchField, volMesh > &field)
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.
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.
void save(const Eigen::Matrix< T, -1, dim > &mat, const std::string fname)
Eigen::Matrix< T, -1, dim > load(Eigen::Matrix< T, -1, dim > &mat, const std::string fname)
volScalarField rAU(1.0/Ueqn.A())
constrainPressure(p, rho, U, phiHbyA, rhorAtU, MRF)
surfaceScalarField rhorAUf("rhorAUf", fvc::interpolate(rho *rAU))
Header file of the torch2Eigen namespace.