Loading...
Searching...
No Matches
02compBump.C
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 Info << "Computing the coefficients for U train" << endl;
122 Eigen::MatrixXd coeffL2U_train = ITHACAutilities::getCoeffs(UfieldTrain,
123 Umodes,
124 0, true);
125 Info << "Computing the coefficients for p train" << endl;
126 Eigen::MatrixXd coeffL2P_train = ITHACAutilities::getCoeffs(PfieldTrain,
127 Pmodes,
128 0, true);
129 Info << "Computing the coefficients for nuT train" << 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 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());
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);
432 ITHACAstream::exportSolution(nut, 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 );
456 ITHACAutilities::getPointsFromPatch(_mesh(), 0, top0, top0_ind);
457 ITHACAutilities::getPointsFromPatch(_mesh(), 1, bot0, bot0_ind);
458 // Info << _mesh().points().size() << 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 // Info << parTop << endl;
509 List<vector> top0_cur = moveBasis(top0, parTop);
510 List<vector> bot0_cur = moveBasis(bot0, parBot);
511 ITHACAutilities::setIndices2Value(top0_ind, top0_cur, movingIDs, curX);
512 ITHACAutilities::setIndices2Value(bot0_ind, bot0_cur, movingIDs, curX);
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 //Info << example.inletIndex.rows() << 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");
750 ITHACAstream::read_fields(onlineNut, 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 << endl;
828 }
829
830 return 0;
831}
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.
Definition 02compBump.C:48
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...
Definition Modes.C:276
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.
Definition Modes.C:63
CompressibleSteadyNN * problem
Full problem.
Definition 02compBump.C:223
ReducedCompressibleSteadyNN(CompressibleSteadyNN &FOMproblem)
Constructor.
Definition 02compBump.C:216
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.
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.
Definition steadyNS.H:302
autoPtr< simpleControl > _simple
simpleControl
Definition steadyNS.H:293
PtrList< volScalarField > Pfield
List of pointers used to form the pressure snapshots matrix.
Definition steadyNS.H:86
autoPtr< fv::options > _fvOptions
fvOptions
Definition steadyNS.H:296
autoPtr< Time > _runTime
Time.
Definition steadyNS.H:299
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
label NUmodes
Number of velocity modes used for the projection.
Definition steadyNS.H:143
label NNutModes
Number of nut modes used for the projection.
Definition steadyNS.H:152
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
scalar cumulativeContErr
continuity error
Definition steadyNS.H:323
autoPtr< volScalarField > _p
Pressure field.
Definition steadyNS.H:269
Class where the tutorial number 2 is implemented.
void offlineSolve(word folder="./ITHACAoutput/Offline/")
Perform an Offline solve.
Definition 02compBump.C:520
tutorial02(int argc, char *argv[])
Constructor.
Definition 02compBump.C:441
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.
Definition ITHACAPOD.C:93
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.
Definition ITHACAnorm.C:300
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.