Loading...
Searching...
No Matches
02compBump.C
Go to the documentation of this file.
1/*---------------------------------------------------------------------------*\
2 ██╗████████╗██╗ ██╗ █████╗ ██████╗ █████╗ ███████╗██╗ ██╗
3 ██║╚══██╔══╝██║ ██║██╔══██╗██╔════╝██╔══██╗ ██╔════╝██║ ██║
4 ██║ ██║ ███████║███████║██║ ███████║█████╗█████╗ ██║ ██║
5 ██║ ██║ ██╔══██║██╔══██║██║ ██╔══██║╚════╝██╔══╝ ╚██╗ ██╔╝
6 ██║ ██║ ██║ ██║██║ ██║╚██████╗██║ ██║ ██║ ╚████╔╝
7 ╚═╝ ╚═╝ ╚═╝ ╚═╝╚═╝ ╚═╝ ╚═════╝╚═╝ ╚═╝ ╚═╝ ╚═══╝
8
9 * In real Time Highly Advanced Computational Applications for Finite Volumes
10 * Copyright (C) 2017 by the ITHACA-FV authors
11-------------------------------------------------------------------------------
12License
13 This file is part of ITHACA-FV
14 ITHACA-FV is free software: you can redistribute it and/or modify
15 it under the terms of the GNU Lesser General Public License as published by
16 the Free Software Foundation, either version 3 of the License, or
17 (at your option) any later version.
18 ITHACA-FV is distributed in the hope that it will be useful,
19 but WITHOUT ANY WARRANTY; without even the implied warranty of
20 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
21 GNU Lesser General Public License for more details.
22 You should have received a copy of the GNU Lesser General Public License
23 along with ITHACA-FV. If not, see <http://www.gnu.org/licenses/>.
24Description
25 Example of steady NS Reduction Problem
26SourceFiles
27 03steadyNS.C
28\*---------------------------------------------------------------------------*/
29
30#include <torch/script.h>
31#include <torch/torch.h>
32#include "torch2Eigen.H"
35#include "RBFMotionSolver.H"
36#include "ITHACAstream.H"
37#include "ITHACAPOD.H"
38#include "forces.H"
39#include "IOmanip.H"
40#include "Foam2Eigen.H"
41
42using namespace ITHACAtorch::torch2Eigen;
43
45{
46 public:
48 CompressibleSteadyNN(int argc, char* argv[])
49 :
50 CompressibleSteadyNS(argc, argv)
51 {
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));
62 };
63
64 label NUmodes;
65 label NNutModes;
66
68 // Eddy viscosity initialization //
70
71 Eigen::MatrixXd bias_inp;
72 Eigen::MatrixXd scale_inp;
73 Eigen::MatrixXd bias_out;
74 Eigen::MatrixXd scale_out;
75
76 Eigen::MatrixXd coeffL2Nut;
77 Eigen::MatrixXd coeffL2U;
78
79 torch::Tensor coeffL2U_tensor;
80 torch::Tensor coeffL2Nut_tensor;
81
82 torch::nn::Sequential Net;
83 torch::optim::Optimizer* optimizer;
84 torch::jit::script::Module netTorchscript;
85
86 void loadNet(word filename)
87 {
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";
90 M_Assert(ITHACAutilities::check_file(filename), Msg.c_str());
91 netTorchscript = torch::jit::load(filename);
92 cnpy::load(bias_inp, "ITHACAoutput/NN/minAnglesInp_" + name(
93 NUmodes) + "_" + name(NNutModes) + ".npy");
94 cnpy::load(scale_inp, "ITHACAoutput/NN/scaleAnglesInp_" + name(
95 NUmodes) + "_" + name(NNutModes) + ".npy");
96 cnpy::load(bias_out, "ITHACAoutput/NN/minOut_" + name(NUmodes) + "_" + name(
97 NNutModes) + ".npy");
98 cnpy::load(scale_out, "ITHACAoutput/NN/scaleOut_" + name(NUmodes) + "_" + name(
99 NNutModes) + ".npy");
100 netTorchscript.eval();
101 }
102
103 // This function computes the coefficients which are later used for training
105 {
106 if (!ITHACAutilities::check_folder("ITHACAoutput/NN/coeffs"))
107 {
108 mkDir("ITHACAoutput/NN/coeffs");
109 // Read Fields for Train
110 PtrList<volVectorField> UfieldTrain;
111 PtrList<volScalarField> PfieldTrain;
112 PtrList<volScalarField> nutFieldsTrain;
113 ITHACAstream::readMiddleFields(UfieldTrain, _U(),
114 "./ITHACAoutput/Offline/");
115 ITHACAstream::readMiddleFields(PfieldTrain, _p(),
116 "./ITHACAoutput/Offline/");
117 auto nut_train = _mesh().lookupObject<volScalarField>("nut");
118 ITHACAstream::readMiddleFields(nutFieldsTrain, nut_train,
119 "./ITHACAoutput/Offline/");
121 std::cout << "Computing the coefficients for U train" << std::endl;
122 Eigen::MatrixXd coeffL2U_train = ITHACAutilities::getCoeffs(UfieldTrain,
123 Umodes,
124 0, true);
125 std::cout << "Computing the coefficients for p train" << std::endl;
126 Eigen::MatrixXd coeffL2P_train = ITHACAutilities::getCoeffs(PfieldTrain,
127 Pmodes,
128 0, true);
129 std::cout << "Computing the coefficients for nuT train" << std::endl;
130 Eigen::MatrixXd coeffL2Nut_train = ITHACAutilities::getCoeffs(nutFieldsTrain,
131 nutModes,
132 0, true);
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");
139 // Read Fields for Test
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");
149 ITHACAstream::readMiddleFields(nutFieldsTest, nut_test,
150 "./ITHACAoutput/checkOff/");
151 // Compute the coefficients for test
152 Eigen::MatrixXd coeffL2U_test = ITHACAutilities::getCoeffs(UfieldTest,
153 Umodes,
154 0, true);
155 Eigen::MatrixXd coeffL2P_test = ITHACAutilities::getCoeffs(PfieldTest,
156 Pmodes,
157 0, true);
158 Eigen::MatrixXd coeffL2Nut_test = ITHACAutilities::getCoeffs(nutFieldsTest,
159 nutModes,
160 0, true);
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");
167 }
168 }
169
170 Eigen::MatrixXd evalNet(Eigen::MatrixXd a, Eigen::MatrixXd mu_now)
171 {
172 Eigen::MatrixXd xpred(a.rows() + mu_now.rows(), 1);
173
174 if (xpred.cols() != 1)
175 {
176 throw std::runtime_error("Predictions for more than one sample not supported yet.");
177 }
178
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);
184 torch::Tensor out;
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();
191 return g;
192 }
193
194 // Function to eval the NN once the input is provided
195 Eigen::MatrixXd evalNet(Eigen::MatrixXd a)
196 {
197 netTorchscript.eval();
198 a.transposeInPlace();
199 a = (a.array() - bias_inp.array()) / scale_inp.array();
200 torch::Tensor aTensor = eigenMatrix2torchTensor(a);
201 torch::Tensor out;// = Net->forward(aTensor);
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();
208 return g;
209 }
210};
211
213{
214 public:
217 :
218 problem(&FOMproblem),
220 {}
221
224
225 void projectReducedOperators(int NmodesUproj, int NmodesPproj, int NmodesEproj)
226 {
227 PtrList<volVectorField> gradModP;
228
229 for (label i = 0; i < NmodesPproj; i++)
230 {
231 gradModP.append(fvc::grad(problem->Pmodes[i]));
232 }
233
234 projGradModP = problem->Umodes.project(gradModP,
235 NmodesUproj); // Modes without lifting
236 }
237
238 void solveOnlineCompressible(int NmodesUproj, int NmodesPproj, int NmodesEproj,
239 int NmodesNutProj, Eigen::MatrixXd mu_now,
240 word Folder = "./ITHACAoutput/Online/")
241 {
242 counter++;
243 // Residuals initialization
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(),
250 NmodesUproj));
251 Eigen::VectorXd eResidual(Eigen::Map<Eigen::VectorXd>(eResidualOld.data(),
252 NmodesEproj));
253 Eigen::VectorXd pResidual(Eigen::Map<Eigen::VectorXd>(pResidualOld.data(),
254 NmodesPproj));
255 // Parameters definition
257 float residualJumpLim =
258 para->ITHACAdict->lookupOrDefault<float>("residualJumpLim", 1e-5);
259 float normalizedResidualLim =
260 para->ITHACAdict->lookupOrDefault<float>("normalizedResidualLim", 1e-5);
261 int maxIter =
262 para->ITHACAdict->lookupOrDefault<float>("maxIter", 2000);
263 bool closedVolume = false;
264 label csolve = 0;
265 // Full variables initialization
266 fluidThermo& thermo = problem->pThermo();
267 volVectorField& U = problem->_U();
268 volScalarField& P = problem->_p();
269 volScalarField& E = problem->pThermo->he();
270 volScalarField& nut = const_cast<volScalarField&>
271 (problem->_mesh().lookupObject<volScalarField>("nut"));
272 volScalarField& rho = problem->_rho();
273 volScalarField& psi = problem->_psi();
274 surfaceScalarField& phi = problem->_phi();
275 Time& runTime = problem->_runTime();
276 fvMesh& mesh = problem->_mesh();
277 fv::options& fvOptions = problem->_fvOptions();
278 scalar cumulativeContErr = problem->cumulativeContErr;
279 // Reduced variables initialization
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);
283 Eigen::MatrixXd nutCoeff = ITHACAutilities::getCoeffs(nut, problem->nutModes,
284 NmodesNutProj, true);
285 //vector Uinlet(170,0,0); // Vector for the inlet boundary condition
286 label idInl =
287 problem->_mesh().boundaryMesh().findPatchID("inlet"); // ID of the inlet patch
288 vector Uinlet(problem->_U().boundaryFieldRef()[idInl][0][0], 0, 0);
289 P.storePrevIter();
290
291 while ((residualJump > residualJumpLim
292 || residualNorm > normalizedResidualLim) && csolve < maxIter)
293 {
294 csolve++;
295 Info << "csolve:" << csolve << endl;
296#if OFVER == 6
297 problem->_simple().loop(runTime);
298#else
299 problem->_simple().loop();
300#endif
301 uResidualOld = uResidual;
302 eResidualOld = eResidual;
303 pResidualOld = pResidual;
304 //Momentum equation phase
305 List<Eigen::MatrixXd> RedLinSysU;
306 ITHACAutilities::assignBC(U, idInl, Uinlet);
307 fvVectorMatrix UEqnR
308 (
309 fvm::div(phi, U)
310 - fvc::div((rho * problem->turbulence->nuEff()) * dev2(T(fvc::grad(U))))
311 - fvm::laplacian(rho * problem->turbulence->nuEff(), U)
312 ==
313 fvOptions(rho, U)
314 );
315 UEqnR.relax();
316 fvOptions.constrain(UEqnR);
317 RedLinSysU = problem->Umodes.project(UEqnR,
318 NmodesUproj); // Modes without lifting
319 Eigen::MatrixXd projGradP = projGradModP * p;
320 RedLinSysU[1] = RedLinSysU[1] - projGradP;
321 u = reducedProblem::solveLinearSys(RedLinSysU, u,
322 uResidual);//, "fullPivLu");//"bdcSvd");
323 problem->Umodes.reconstruct(U, u, "U");
324 ITHACAutilities::assignBC(U, idInl, Uinlet);
325 //solve(UEqnR == -problem->getGradP(P)); //For debug purposes only, second part only useful when using uEqn_global==-getGradP
326 fvOptions.correct(U);
327 //Energy equation phase
328 fvScalarMatrix EEqnR
329 (
330 fvm::div(phi, E)
331 + fvc::div(phi, volScalarField("Ekp", 0.5 * magSqr(U) + P / rho))
332 - fvm::laplacian(problem->turbulence->alphaEff(), E)
333 ==
334 fvOptions(rho, E)
335 );
336 EEqnR.relax();
337 fvOptions.constrain(EEqnR);
338 List<Eigen::MatrixXd> RedLinSysE = problem->Emodes.project(EEqnR, NmodesEproj);
339 e = reducedProblem::solveLinearSys(RedLinSysE, e, eResidual);
340 problem->Emodes.reconstruct(E, e, "e");
341 //EEqnR.solve(); //For debug purposes only
342 fvOptions.correct(E);
343 thermo.correct(); // Here are calculated both temperature and density based on P,U and he.
344 // Pressure equation phase
345 constrainPressure(P, rho, U, problem->getPhiHbyA(UEqnR, U, P),
346 problem->getRhorAUf(
347 UEqnR));// Update the pressure BCs to ensure flux consistency
348 surfaceScalarField phiHbyACalculated = problem->getPhiHbyA(UEqnR, U, P);
349 closedVolume = adjustPhi(phiHbyACalculated, U, P);
350 List<Eigen::MatrixXd> RedLinSysP;
351
352 while (problem->_simple().correctNonOrthogonal())
353 {
354 volScalarField rAU(1.0 /
355 UEqnR.A()); // Inverse of the diagonal part of the U equation matrix
356 volVectorField HbyA(constrainHbyA(rAU * UEqnR.H(), U,
357 P)); // H is the extra diagonal part summed to the r.h.s. of the U equation
358 surfaceScalarField phiHbyA("phiHbyA", fvc::interpolate(rho)*fvc::flux(HbyA));
359 surfaceScalarField rhorAUf("rhorAUf", fvc::interpolate(rho * rAU));
360 fvScalarMatrix PEqnR
361 (
362 fvc::div(phiHbyA)
363 - fvm::laplacian(rhorAUf, P)
364 ==
365 fvOptions(psi, P, rho.name())
366 );
367 PEqnR.setReference
368 (
369 problem->_pressureControl().refCell(),
370 problem->_pressureControl().refValue()
371 );
372 RedLinSysP = problem->Pmodes.project(PEqnR, NmodesPproj);
373 p = reducedProblem::solveLinearSys(RedLinSysP, p, pResidual);
374 problem->Pmodes.reconstruct(P, p, "p");
375
376 if (problem->_simple().finalNonOrthogonalIter())
377 {
378 phi = problem->getPhiHbyA(UEqnR, U, P) + PEqnR.flux();
379 }
380 }
381
382#include "continuityErrs.H"
383 P.relax();// Explicitly relax pressure for momentum corrector
384 U = problem->HbyA() - (1.0 / UEqnR.A()) * problem->getGradP(P);
385 U.correctBoundaryConditions();
386 fvOptions.correct(U);
387 bool pLimited = problem->_pressureControl().limit(P);
388
389 // For closed-volume cases adjust the pressure and density levels to obey overall mass continuity
390 if (closedVolume)
391 {
392 P += (problem->_initialMass() - fvc::domainIntegrate(psi * P))
393 / fvc::domainIntegrate(psi);
394 }
395
396 if (pLimited || closedVolume)
397 {
398 P.correctBoundaryConditions();
399 }
400
401 rho = thermo.rho(); // Here rho is calculated as p*psi = p/(R*T)
402 rho.relax();
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());
419 //problem->turbulence->correct(); // Resolution of the full turbulence (debug purposes only)
420 }
421
422 nutCoeff = problem->evalNet(u, mu_now);
423 problem->nutModes.reconstruct(nut, nutCoeff, "nut");
424 label k = 1;
425 // U.rename("Ur");
426 // P.rename("Pr");
427 // E.rename("Er");
428 // nut.rename("nutR");
429 ITHACAstream::exportSolution(U, name(counter), Folder);
430 ITHACAstream::exportSolution(P, name(counter), Folder);
431 ITHACAstream::exportSolution(E, name(counter), Folder);
433 }
434
435};
436
437class tutorial02 : public CompressibleSteadyNN
438{
439 public:
441 explicit tutorial02(int argc, char* argv[])
442 :
443 CompressibleSteadyNN(argc, argv)
444 {
445 dyndict = new IOdictionary
446 (
447 IOobject
448 (
449 "dynamicMeshDictRBF",
450 "./constant",
451 _mesh(),
452 IOobject::MUST_READ,
453 IOobject::NO_WRITE
454 )
455 );
458 // std::cout << _mesh().points().size() << std::endl;
459 ms = new RBFMotionSolver(_mesh(), *dyndict);
460 vectorField motion(ms->movingPoints().size(), vector::zero);
461 movingIDs = ms->movingIDs();
462 x0 = ms->movingPoints();
463 curX = x0;
464 point0 = ms->curPoints();
466 middleExport = para->ITHACAdict->lookupOrDefault<bool>("middleExport", true);
467 }
468
469 List<vector> top0;
470 List<vector> bot0;
471 labelList top0_ind;
472 labelList bot0_ind;
473 IOdictionary* dyndict;
474 RBFMotionSolver* ms;
475 labelList movingIDs;
476 List<vector> x0;
477 List<vector> curX;
478 vectorField point0;
479 vectorField point;
481
482
483 double f1(double chord, double x)
484 {
485 double res = chord * (std::pow((x) / chord,
486 0.5) * (1 - (x) / chord)) / (std::exp(15 * (x) / chord));
487 return res;
488 }
489
490 List<vector> moveBasis(const List<vector>& originalPoints, double par)
491 {
492 List<vector> movedPoints(originalPoints);
493
494 for (int i = 0; i < originalPoints.size(); i++)
495 {
496 movedPoints[i][2] += par * f1(1, movedPoints[i][0]);
497 }
498
499 return movedPoints;
500 }
501
502 void updateMesh(double parTop = 0, double parBot = 0)
503 {
504 _mesh().movePoints(point0);
505
506 if (parTop != 0 || parBot != 0)
507 {
508 // std::cout << parTop << std::endl;
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);
516 }
517 }
518
520 void offlineSolve(word folder = "./ITHACAoutput/Offline/")
521 {
523 volVectorField& U = _U();
525 volScalarField& p = _p();
527 volScalarField& E = _E();
528
529 // If the offline solution is already performed but POD modes are not present, then read the fields
530 if (offline && !ITHACAutilities::check_folder("./ITHACAoutput/POD/1"))
531 {
536 auto nut = _mesh().lookupObject<volScalarField>("nut");
538 mu_samples = ITHACAstream::readMatrix("./parsOff_mat.txt");
539 }
540 // If offline stage ha snot been performed, then perform it
541 else if (!offline)
542 {
543 double UIFinit = para->ITHACAdict->lookupOrDefault<double>("UIFinit", 170);
544 Vector<double> Uinl(UIFinit, 0, 0);
545
546 for (label i = 0; i < mu.rows(); i++)
547 {
548 updateMesh(mu(i, 0), mu(i, 1));
549 ITHACAstream::writePoints(_mesh().points(), folder, name(i + 1) + "/polyMesh/");
550 assignIF(_U(), Uinl);
551 truthSolve(folder);
552 label j = 1;
553 word polyMesh2beLinked = folder + name(i + 1) + "/" + "polyMesh/";
554
555 while (ITHACAutilities::check_folder(folder + name(i + 1) + "/" + name(j)))
556 {
557 word folderContLink = folder + name(i + 1) + "/" + name(j) + "/";
558 system("ln -s $(readlink -f " + polyMesh2beLinked + ") " + folderContLink +
559 " >/dev/null 2>&1");
560 j++;
561 }
562 }
563 }
564 }
565};
566
567int main(int argc, char* argv[])
568{
569 // Construct the tutorial object
570 tutorial02 example(argc, argv);
572 std::ifstream exFileOff("./parsOff_mat.txt");
573
574 if (exFileOff)
575 {
576 example.mu = ITHACAstream::readMatrix("./parsOff_mat.txt");
577 }
578 else
579 {
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);
583 Eigen::MatrixXd parTop = ITHACAutilities::rand(example.mu.rows(), 1, 0,
584 BumpAmp);
585 Eigen::MatrixXd parBot = ITHACAutilities::rand(example.mu.rows(), 1, -BumpAmp,
586 0);
587 example.mu.leftCols(1) = parTop;
588 example.mu.rightCols(1) = parBot;
589 ITHACAstream::exportMatrix(example.mu, "parsOff", "eigen", "./");
590 }
591
592 Eigen::MatrixXd parsOn;
593 std::ifstream exFileOn("./parsOn_mat.txt");
594
595 if (exFileOn)
596 {
597 parsOn = ITHACAstream::readMatrix("./parsOn_mat.txt");
598 }
599 else
600 {
601 int OnNum = para->ITHACAdict->lookupOrDefault<int>("OnNum", 20);
602 double BumpAmp = para->ITHACAdict->lookupOrDefault<double>("BumpAmp", 0.1);
603 parsOn.resize(OnNum, 2);
604 Eigen::MatrixXd parTopOn = ITHACAutilities::rand(OnNum, 1, 0, BumpAmp);
605 Eigen::MatrixXd parBotOn = ITHACAutilities::rand(OnNum, 1, -BumpAmp, 0);
606 parsOn.leftCols(1) = parTopOn;
607 parsOn.rightCols(1) = parBotOn;
608 ITHACAstream::exportMatrix(parsOn, "parsOn", "eigen", "./");
609 }
610
611 // Read some parameters from file
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();
621 //Perform the offline solve
622 example.offlineSolve();
623 // Move the mesh to the original geometry to get the modes into a mid mesh
624 example.updateMesh();
625 // Perform POD on velocity and pressure
626 ITHACAPOD::getModes(example.Ufield, example.Umodes, example._U().name(),
627 example.podex, 0, 0,
628 NmodesUout);
629 ITHACAPOD::getModes(example.Pfield, example.Pmodes, example._p().name(),
630 example.podex, 0, 0,
631 NmodesPout);
632 ITHACAPOD::getModes(example.Efield, example.Emodes, example._E().name(),
633 example.podex, 0, 0,
634 NmodesEout);
635 ITHACAPOD::getModes(example.nutFields, example.nutModes, "nut", example.podex,
636 0, 0,
637 NmodesNutOut);
638
639 if (!ITHACAutilities::check_folder("./ITHACAoutput/checkOff"))
640 {
641 example.updateMesh();
642 tutorial02 checkOff(argc, argv);
643 checkOff.mu = ITHACAstream::readMatrix("./parsOn_mat.txt");
644 //Perform the offline solve
645 checkOff.offline = false;
646 checkOff.restart();
647 checkOff.offlineSolve("./ITHACAoutput/checkOff/");
648 checkOff.offline = true;
649 // Eigen::MatrixXd snapsCheck = ITHACAstream::readMatrix("./ITHACAoutput/checkOff/snaps");
650 // label fieldNum = 0;
651 // for(label k=0; k<snapsCheck.rows(); k++)
652 // {
653 // fieldNum = fieldNum+snapsCheck(k,0);
654 // ITHACAstream::exportSolution(checkOff.Ufield[fieldNum-1], name(k+1), "./ITHACAoutput/checkOffSingle/");
655 // ITHACAstream::exportSolution(checkOff.Pfield[fieldNum-1], name(k+1), "./ITHACAoutput/checkOffSingle/");
656 // ITHACAstream::exportSolution(checkOff.Efield[fieldNum-1], name(k+1), "./ITHACAoutput/checkOffSingle/");
657 // ITHACAstream::exportSolution(checkOff.nutFields[fieldNum-1], name(k+1), "./ITHACAoutput/checkOffSingle/");
658 // ITHACAutilities::createSymLink("./ITHACAoutput/checkOff/"+name(k+1)+"/polyMesh", "./ITHACAoutput/checkOffSingle/"+name(k+1)+"/");
659 // }
660 // ITHACAutilities::createSymLink("./0", "./ITHACAoutput/checkOffSingle/");
661 // ITHACAutilities::createSymLink("./system", "./ITHACAoutput/checkOffSingle/");
662 // ITHACAutilities::createSymLink("./constant", "./ITHACAoutput/checkOffSingle/");
663 }
664
665 // Create the coefficients to train the net
666 example.getTurbNN();
667 //Before loading the net, it has to be created through the python script
668 example.loadNet("ITHACAoutput/NN/Net_" + name(example.NUmodes) + "_" + name(
669 example.NNutModes) + ".pt");
670 // Create the reduced object
671 ReducedCompressibleSteadyNN reduced(example);
672
673 //Perform the online solutions
674 for (label k = 0; k < parsOn.rows(); k++)
675 {
676 example.updateMesh();
677 example.updateMesh(parsOn(k, 0), parsOn(k, 1));
678 ITHACAstream::writePoints(example._mesh().points(),
679 "./ITHACAoutput/Online_" + name(NmodesUproj) + "_" + name(NmodesNutProj) + "/",
680 name(k + 1) + "/polyMesh/");
681 //std::cout << example.inletIndex.rows() << std::endl;
682 //reduced.setOnlineVelocity(vel);
683 reduced.projectReducedOperators(NmodesUproj, NmodesPproj, NmodesEproj);
684 example.restart();
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) + "/");
691 }
692
693 if (ITHACAutilities::check_folder("./ITHACAoutput/checkOff"))
694 {
695 PtrList<volVectorField> UfieldCheck;
696 PtrList<volScalarField> PfieldCheck;
697 PtrList<volScalarField> EfieldCheck;
698 PtrList<volScalarField> nutFieldsCheck;
699 ITHACAstream::readMiddleFields(UfieldCheck, example._U(),
700 "./ITHACAoutput/checkOff/");
701 ITHACAstream::readMiddleFields(PfieldCheck, example._p(),
702 "./ITHACAoutput/checkOff/");
703 ITHACAstream::readMiddleFields(EfieldCheck, example._E(),
704 "./ITHACAoutput/checkOff/");
705 auto nutCheck = example._mesh().lookupObject<volScalarField>("nut");
706 ITHACAstream::readMiddleFields(nutFieldsCheck, nutCheck,
707 "./ITHACAoutput/checkOff/");
708 // Info << "UfieldCheck.size" << UfieldCheck.size() << endl;
709 // Info << "PfieldCheck.size" << PfieldCheck.size() << endl;
710 // Info << "EfieldCheck.size" << EfieldCheck.size() << endl;
711 // Info << "nutFieldsCheck.size" << nutFieldsCheck.size() << endl;
712 // exit(0);
713 Eigen::MatrixXd snapsCheck =
714 ITHACAstream::readMatrix("./ITHACAoutput/checkOff/snaps");
715 label fieldNum = 0;
716
717 for (label k = 0; k < snapsCheck.rows(); k++)
718 {
719 fieldNum = fieldNum + snapsCheck(k, 0);
720 ITHACAstream::exportSolution(UfieldCheck[fieldNum - 1], name(k + 1),
721 "./ITHACAoutput/checkOffSingle/");
722 ITHACAstream::exportSolution(PfieldCheck[fieldNum - 1], name(k + 1),
723 "./ITHACAoutput/checkOffSingle/");
724 ITHACAstream::exportSolution(EfieldCheck[fieldNum - 1], name(k + 1),
725 "./ITHACAoutput/checkOffSingle/");
726 ITHACAstream::exportSolution(nutFieldsCheck[fieldNum - 1], name(k + 1),
727 "./ITHACAoutput/checkOffSingle/");
728 ITHACAutilities::createSymLink("./ITHACAoutput/checkOff/" + name(
729 k + 1) + "/polyMesh", "./ITHACAoutput/checkOffSingle/" + name(k + 1) + "/");
730 }
731
732 ITHACAutilities::createSymLink("./0", "./ITHACAoutput/checkOffSingle/");
733 ITHACAutilities::createSymLink("./system", "./ITHACAoutput/checkOffSingle/");
734 ITHACAutilities::createSymLink("./constant", "./ITHACAoutput/checkOffSingle/");
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;
743 ITHACAstream::read_fields(onlineU, example._U(),
744 "./ITHACAoutput/Online_" + name(NmodesUproj) + "_" + name(NmodesNutProj) + "/");
745 ITHACAstream::read_fields(onlineP, example._p(),
746 "./ITHACAoutput/Online_" + name(NmodesUproj) + "_" + name(NmodesNutProj) + "/");
747 ITHACAstream::read_fields(onlineE, example._E(),
748 "./ITHACAoutput/Online_" + name(NmodesUproj) + "_" + name(NmodesNutProj) + "/");
749 auto nut = example._mesh().lookupObject<volScalarField>("nut");
751 "./ITHACAoutput/Online_" + name(NmodesUproj) + "_" + name(NmodesNutProj) + "/");
752 ITHACAstream::read_fields(offlineU, example._U(),
753 "./ITHACAoutput/checkOffSingle/");
754 ITHACAstream::read_fields(offlineP, example._p(),
755 "./ITHACAoutput/checkOffSingle/");
756 ITHACAstream::read_fields(offlineE, example._E(),
757 "./ITHACAoutput/checkOffSingle/");
758 ITHACAstream::read_fields(offlineNut, nut, "./ITHACAoutput/checkOffSingle/");
759 Eigen::MatrixXd errorU = ITHACAutilities::errorL2Rel(offlineU,
760 onlineU);
761 Eigen::MatrixXd errorP = ITHACAutilities::errorL2Rel(offlineP,
762 onlineP);
763 Eigen::MatrixXd errorE = ITHACAutilities::errorL2Rel(offlineE,
764 onlineE);
765 Eigen::MatrixXd errorNut = ITHACAutilities::errorL2Rel(offlineNut,
766 onlineNut);
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) + "/");
784
785 for (label j = 0; j < parsOn.rows(); j++)
786 {
787 volVectorField Ue = offlineU[j] - onlineU[j];
788 //auto offU2Eigen = Foam2Eigen::field2Eigen(offlineU[j]);
789 auto u = ITHACAutilities::L2Norm(Ue);
790 Ue /= u;
792 volScalarField Pe = offlineP[j] - onlineP[j];
793 auto p = ITHACAutilities::L2Norm(Pe);
794 Pe /= p;
796 volScalarField Ee = offlineE[j] - onlineE[j];
797 auto e = ITHACAutilities::L2Norm(Ee);
798 Ee /= e;
799 volScalarField Nute = offlineNut[j] - onlineNut[j];
800 auto n = ITHACAutilities::L2Norm(Nute);
801 Nute /= n;
802 Ue.rename("Ue");
803 Pe.rename("Pe");
804 Ee.rename("Ee");
805 Nute.rename("Nute");
806 ITHACAstream::exportSolution(Ue, name(j + 1),
807 "./ITHACAoutput/ErrorFields_" + name(NmodesUproj) + "_" + name(
808 NmodesNutProj) + "/");
809 ITHACAstream::exportSolution(Pe, name(j + 1),
810 "./ITHACAoutput/ErrorFields_" + name(NmodesUproj) + "_" + name(
811 NmodesNutProj) + "/");
812 ITHACAstream::exportSolution(Ee, name(j + 1),
813 "./ITHACAoutput/ErrorFields_" + name(NmodesUproj) + "_" + name(
814 NmodesNutProj) + "/");
815 ITHACAstream::exportSolution(Nute, name(j + 1),
816 "./ITHACAoutput/ErrorFields_" + name(NmodesUproj) + "_" + name(
817 NmodesNutProj) + "/");
818 // ITHACAstream::exportSolution(Ue, name(j+1), "./ITHACAoutput/checkOffSingle/");
819 // ITHACAstream::exportSolution(Pe, name(j+1), "./ITHACAoutput/checkOffSingle/");
820 // ITHACAstream::exportSolution(Ee, name(j+1), "./ITHACAoutput/checkOffSingle/");
821 // ITHACAstream::exportSolution(Nute, name(j+1), "./ITHACAoutput/checkOffSingle/");
822 }
823 }
824 else
825 {
826 std::cerr << "checkOff folder is missing, error analysis cannot be performed."
827 << std::endl;
828 }
829
830 exit(0);
831}
int main(int argc, char *argv[])
Definition 02compBump.C:567
Header file of the steadyNS class.
Header file of the Foam2Eigen class.
Header file of the ITHACAPOD class.
Foam::fvMesh & mesh
Definition createMesh.H:47
Foam::Time & runTime
Definition createTime.H:33
#define M_Assert(Expr, Msg)
Header file of the ITHACAstream class, it contains the implementation of several methods for input ou...
fv::options & fvOptions
Definition NLsolve.H:25
label csolve
Definition NLsolve.H:42
maxIter
Definition NLsolve.H:75
bool closedVolume
Definition NLsolve.H:33
volScalarField & psi
Definition NLsolve.H:29
ITHACAparameters * para
Definition NLsolve.H:40
auto nut
Definition NLsolve.H:195
Header file of the reducedSteadyNS class.
CompressibleSteadyNN(int argc, char *argv[])
Constructor.
Definition 02compBump.C:48
Eigen::MatrixXd scale_inp
Definition 02compBump.C:72
Eigen::MatrixXd bias_out
Definition 02compBump.C:73
Eigen::MatrixXd coeffL2Nut
Definition 02compBump.C:76
Eigen::MatrixXd coeffL2U
Definition 02compBump.C:77
Eigen::MatrixXd bias_inp
Definition 02compBump.C:71
torch::optim::Optimizer * optimizer
Definition 02compBump.C:83
Eigen::MatrixXd evalNet(Eigen::MatrixXd a, Eigen::MatrixXd mu_now)
Definition 02compBump.C:170
torch::nn::Sequential Net
Definition 02compBump.C:82
torch::Tensor coeffL2U_tensor
Definition 02compBump.C:79
torch::Tensor coeffL2Nut_tensor
Definition 02compBump.C:80
Eigen::MatrixXd scale_out
Definition 02compBump.C:74
torch::jit::script::Module netTorchscript
Definition 02compBump.C:84
void loadNet(word filename)
Definition 02compBump.C:86
Eigen::MatrixXd evalNet(Eigen::MatrixXd a)
Definition 02compBump.C:195
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)
Definition 02compBump.C:225
CompressibleSteadyNN * problem
Full problem.
Definition 02compBump.C:223
ReducedCompressibleSteadyNN(CompressibleSteadyNN &FOMproblem)
Constructor.
Definition 02compBump.C:216
void solveOnlineCompressible(int NmodesUproj, int NmodesPproj, int NmodesEproj, int NmodesNutProj, Eigen::MatrixXd mu_now, word Folder="./ITHACAoutput/Online/")
Definition 02compBump.C:238
Eigen::MatrixXd projGradModP
Projected gradient of the pressure modes.
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.
Definition steadyNS.H:86
volVectorModes Umodes
List of pointers used to form the velocity modes.
Definition steadyNS.H:101
PtrList< volVectorField > Ufield
List of pointers used to form the velocity snapshots matrix.
Definition steadyNS.H:89
autoPtr< fvMesh > _mesh
Mesh.
Definition steadyNS.H:290
ITHACAparameters * para
Definition steadyNS.H:82
autoPtr< volVectorField > Uinl
Initial dummy field with all Dirichlet boundary conditions.
Definition steadyNS.H:281
autoPtr< volVectorField > _U
Velocity field.
Definition steadyNS.H:272
volScalarModes Pmodes
List of pointers used to form the pressure modes.
Definition steadyNS.H:98
Class where the tutorial number 2 is implemented.
void offlineSolve(word folder="./ITHACAoutput/Offline/")
Perform an Offline solve.
Definition 02compBump.C:520
void updateMesh(double parTop=0, double parBot=0)
Definition 02compBump.C:502
vectorField point0
Definition 02compBump.C:478
tutorial02(int argc, char *argv[])
Constructor.
Definition 02compBump.C:441
double f1(double chord, double x)
Definition 02compBump.C:483
List< vector > moveBasis(const List< vector > &originalPoints, double par)
Definition 02compBump.C:490
labelList bot0_ind
Definition 02compBump.C:472
List< vector > bot0
Definition 02compBump.C:470
List< vector > top0
Definition 02compBump.C:469
List< vector > curX
Definition 02compBump.C:477
List< vector > x0
Definition 02compBump.C:476
RBFMotionSolver * ms
Definition 02compBump.C:474
vectorField point
Definition 02compBump.C:479
IOdictionary * dyndict
Definition 02compBump.C:473
labelList movingIDs
Definition 02compBump.C:475
labelList top0_ind
Definition 02compBump.C:471
ITHACAparameters * para
Definition 02compBump.C:480
volScalarField & T
Definition createT.H:46
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.
Definition ITHACAPOD.C:90
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.
Definition torch2Eigen.C:39
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)
Definition ITHACAerror.C:41
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)
phiHbyA
Definition pcEqn.H:61
bool pLimited
Definition pcEqn.H:98
HbyA
Definition pcEqn.H:62
surfaceScalarField & phi
volVectorField & U
cumulativeContErr
volScalarField & p
volScalarField & rho
fluidThermo & thermo
surfaceScalarField rhorAUf("rhorAUf", fvc::interpolate(rho *rAU))
dimensionedVector & g
adjustPhi(phiHbyA, U, p)
label i
Definition pEqn.H:46
Header file of the torch2Eigen namespace.