Loading...
Searching...
No Matches
01simpleTurbGeomClosed.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 turbulent steady NS Reduction Problem solved by the use of the SIMPLE algorithm
26SourceFiles
27 01simpleTurbGeomClosed.C
28\*---------------------------------------------------------------------------*/
29
30#include <torch/script.h>
31#include <torch/torch.h>
32#include "torch2Eigen.H"
33#include "SteadyNSSimple.H"
34#include "ITHACAstream.H"
35#include "ITHACAPOD.H"
36#include "ReducedSimpleSteadyNS.H"
37#include "forces.H"
38#include "IOmanip.H"
39#include "RBFMotionSolver.H"
40
41
42using namespace ITHACAtorch::torch2Eigen;
43
44class SteadyNSSimpleNN : public SteadyNSSimple
45{
46 public:
48 SteadyNSSimpleNN(int argc, char* argv[])
49 :
50 SteadyNSSimple(argc, argv)
51 {
52 Net->push_back(torch::nn::Linear(NUmodes, 128));
53 Net->push_back(torch::nn::ReLU());
54 Net->push_back(torch::nn::Linear(128, 64));
55 Net->push_back(torch::nn::ReLU());
56 Net->push_back(torch::nn::Linear(64, NNutModes));
57 optimizer = new torch::optim::Adam(Net->parameters(),
58 torch::optim::AdamOptions(2e-2));
59 };
60
61 Eigen::MatrixXd bias_inp;
62 Eigen::MatrixXd scale_inp;
63 Eigen::MatrixXd bias_out;
64 Eigen::MatrixXd scale_out;
65
66 Eigen::MatrixXd coeffL2Nut;
67 Eigen::MatrixXd coeffL2U;
68
69 torch::Tensor coeffL2U_tensor;
70 torch::Tensor coeffL2Nut_tensor;
71
72 torch::nn::Sequential Net;
73 torch::optim::Optimizer* optimizer;
74 torch::jit::script::Module netTorchscript;
75
76 void loadNet(word filename)
77 {
78 std::string Msg = filename +
79 " is not existing, please run the training stage of the net with the correct number of modes for U and Nut";
80 M_Assert(ITHACAutilities::check_file(filename), Msg.c_str());
81 netTorchscript = torch::jit::load(filename);
82 cnpy::load(bias_inp, "ITHACAoutput/NN/minAnglesInp_" + name(
83 NUmodes) + "_" + name(NNutModes) + ".npy");
84 cnpy::load(scale_inp, "ITHACAoutput/NN/scaleAnglesInp_" + name(
85 NUmodes) + "_" + name(NNutModes) + ".npy");
86 cnpy::load(bias_out, "ITHACAoutput/NN/minOut_" + name(NUmodes) + "_" + name(
87 NNutModes) + ".npy");
88 cnpy::load(scale_out, "ITHACAoutput/NN/scaleOut_" + name(NUmodes) + "_" + name(
89 NNutModes) + ".npy");
90 netTorchscript.eval();
91 }
92
93 // This function computes the coefficients which are later used for training
94 void getTurbNN()
95 {
96 if (!ITHACAutilities::check_folder("ITHACAoutput/NN/coeffs"))
97 {
98 mkDir("ITHACAoutput/NN/coeffs");
99 // Read Fields for Train
100 PtrList<volVectorField> UfieldTrain;
101 PtrList<volScalarField> PfieldTrain;
102 PtrList<volScalarField> nutFieldsTrain;
103 ITHACAstream::readMiddleFields(UfieldTrain, _U(),
104 "./ITHACAoutput/Offline/");
105 ITHACAstream::readMiddleFields(PfieldTrain, _p(),
106 "./ITHACAoutput/Offline/");
107 auto nut_train = _mesh().lookupObject<volScalarField>("nut");
108 ITHACAstream::readMiddleFields(nutFieldsTrain, nut_train,
109 "./ITHACAoutput/Offline/");
111 std::cout << "Computing the coefficients for U train" << std::endl;
112 Eigen::MatrixXd coeffL2U_train = ITHACAutilities::getCoeffs(UfieldTrain,
113 Umodes,
114 0, true);
115 std::cout << "Computing the coefficients for p train" << std::endl;
116 Eigen::MatrixXd coeffL2P_train = ITHACAutilities::getCoeffs(PfieldTrain,
117 Pmodes,
118 0, true);
119 std::cout << "Computing the coefficients for nuT train" << std::endl;
120 Eigen::MatrixXd coeffL2Nut_train = ITHACAutilities::getCoeffs(nutFieldsTrain,
121 nutModes,
122 0, true);
123 coeffL2U_train.transposeInPlace();
124 coeffL2P_train.transposeInPlace();
125 coeffL2Nut_train.transposeInPlace();
126 cnpy::save(coeffL2U_train, "ITHACAoutput/NN/coeffs/coeffL2UTrain.npy");
127 cnpy::save(coeffL2P_train, "ITHACAoutput/NN/coeffs/coeffL2PTrain.npy");
128 cnpy::save(coeffL2Nut_train, "ITHACAoutput/NN/coeffs/coeffL2NutTrain.npy");
129 // Read Fields for Test
130 PtrList<volVectorField> UfieldTest;
131 PtrList<volScalarField> PfieldTest;
132 PtrList<volScalarField> nutFieldsTest;
134 ITHACAstream::readMiddleFields(UfieldTest, _U(),
135 "./ITHACAoutput/checkOff/");
136 ITHACAstream::readMiddleFields(PfieldTest, _p(),
137 "./ITHACAoutput/checkOff/");
138 auto nut_test = _mesh().lookupObject<volScalarField>("nut");
139 ITHACAstream::readMiddleFields(nutFieldsTest, nut_test,
140 "./ITHACAoutput/checkOff/");
141 // Compute the coefficients for test
142 Eigen::MatrixXd coeffL2U_test = ITHACAutilities::getCoeffs(UfieldTest,
143 Umodes,
144 0, true);
145 Eigen::MatrixXd coeffL2P_test = ITHACAutilities::getCoeffs(PfieldTest,
146 Pmodes,
147 0, true);
148 Eigen::MatrixXd coeffL2Nut_test = ITHACAutilities::getCoeffs(nutFieldsTest,
149 nutModes,
150 0, true);
151 coeffL2U_test.transposeInPlace();
152 coeffL2P_test.transposeInPlace();
153 coeffL2Nut_test.transposeInPlace();
154 cnpy::save(coeffL2U_test, "ITHACAoutput/NN/coeffs/coeffL2UTest.npy");
155 cnpy::save(coeffL2P_test, "ITHACAoutput/NN/coeffs/coeffL2PTest.npy");
156 cnpy::save(coeffL2Nut_test, "ITHACAoutput/NN/coeffs/coeffL2NutTest.npy");
157 }
158 }
159
160 Eigen::MatrixXd evalNet(Eigen::MatrixXd a, double mu_now)
161 {
162 Eigen::MatrixXd xpred(a.rows() + 1, 1);
163
164 if (xpred.cols() != 1)
165 {
166 throw std::runtime_error("Predictions for more than one sample not supported yet.");
167 }
168
169 xpred.bottomRows(a.rows()) = a;
170 xpred(0, 0) = mu_now;
171 xpred = xpred.array() * scale_inp.array() + bias_inp.array() ;
172 xpred.transposeInPlace();
173 torch::Tensor xTensor = eigenMatrix2torchTensor(xpred);
174 torch::Tensor out;
175 std::vector<torch::jit::IValue> XTensorInp;
176 XTensorInp.push_back(xTensor);
177 out = netTorchscript.forward(XTensorInp).toTensor();
178 Eigen::MatrixXd g = torchTensor2eigenMatrix<double>(out);
179 g.transposeInPlace();
180 g = (g.array() - bias_out.array()) / scale_out.array();
181 return g;
182 }
183
184 // Function to eval the NN once the input is provided
185 Eigen::MatrixXd evalNet(Eigen::MatrixXd a)
186 {
187 netTorchscript.eval();
188 a.transposeInPlace();
189 a = (a.array() - bias_inp.array()) / scale_inp.array();
190 torch::Tensor aTensor = eigenMatrix2torchTensor(a);
191 torch::Tensor out;// = Net->forward(aTensor);
192 std::vector<torch::jit::IValue> XTensorInp;
193 XTensorInp.push_back(aTensor);
194 out = netTorchscript.forward(XTensorInp).toTensor();
195 Eigen::MatrixXd g = torchTensor2eigenMatrix<double>(out);
196 g = g.array() * scale_out.array() + bias_out.array();
197 g.transposeInPlace();
198 return g;
199 }
200};
201
202class reducedSimpleSteadyNN : public reducedSimpleSteadyNS
203{
204 public:
206 explicit reducedSimpleSteadyNN(SteadyNSSimpleNN& FOMproblem)
207 :
208 problem(&FOMproblem)
209 {}
210
212 SteadyNSSimpleNN* problem;
213
214 // Function to perform the online phase
215 void solveOnline_Simple(scalar mu_now,
216 int NmodesUproj, int NmodesPproj, int NmodesNut = 0,
217 word Folder = "./ITHACAoutput/Reconstruct/")
218 {
219 counter++;
220
221 // For all the variables, in case one wants to use all the available modes it is just necessary to set
222 // the requested number into the ITHACAdict to zero
223 if (NmodesUproj == 0)
224 {
225 NmodesUproj = problem->Umodes.size();
226 }
227
228 if (NmodesPproj == 0)
229 {
230 NmodesPproj = problem->Pmodes.size();
231 }
232
233 if (NmodesNut == 0)
234 {
235 NmodesNut = problem->nutModes.size();
236 }
237
238 // Initializations
239 Eigen::VectorXd uresidualOld = Eigen::VectorXd::Zero(NmodesUproj);
240 Eigen::VectorXd presidualOld = Eigen::VectorXd::Zero(NmodesPproj);
241 Eigen::VectorXd uresidual = Eigen::VectorXd::Zero(NmodesUproj);
242 Eigen::VectorXd presidual = Eigen::VectorXd::Zero(NmodesPproj);
243 scalar U_norm_res(1);
244 scalar P_norm_res(1);
245 Eigen::MatrixXd a = Eigen::VectorXd::Zero(NmodesUproj);
246 Eigen::MatrixXd a0 = Eigen::VectorXd::Zero(NmodesUproj);
247 Eigen::MatrixXd b = Eigen::VectorXd::Zero(NmodesPproj);
248 Eigen::MatrixXd b0 = Eigen::VectorXd::Zero(NmodesPproj);
249 Eigen::MatrixXd bOld = b;
250 Eigen::MatrixXd nutCoeff = Eigen::VectorXd::Zero(NmodesNut);
251 Eigen::MatrixXd nutCoeffOld = Eigen::VectorXd::Zero(NmodesNut);
252 float residualJumpLim =
253 problem->para->ITHACAdict->lookupOrDefault<float>("residualJumpLim", 1e-5);
254 float normalizedResidualLim =
255 problem->para->ITHACAdict->lookupOrDefault<float>("normalizedResidualLim",
256 1e-5);
257 scalar residual_jump(1 + residualJumpLim);
258 volScalarField& P = problem->_p();
259 volVectorField& U = problem->_U();
260 volScalarField& nut = const_cast<volScalarField&>
261 (problem->_mesh().lookupObject<volScalarField>("nut"));
262 volVectorField u2 = U;
263 a0 = ITHACAutilities::getCoeffs(U, problem->Umodes, NmodesUproj, true);
264 b = ITHACAutilities::getCoeffs(P, problem->Pmodes, NmodesPproj, true);
265 nutCoeff = ITHACAutilities::getCoeffs(nut, problem->nutModes, NmodesNut, true);
266 a(0) = a0(0); // Do not remove: it is not working without this condition
267 b = b0;
268 fvMesh& mesh = problem->_mesh();
269 Time& runTime = problem->_runTime();
270 P.rename("p");
271 surfaceScalarField& phi(problem->_phi());
272 problem->Umodes.reconstruct(U, a, "U");
273 problem->Pmodes.reconstruct(P, b, "p");
274 vector v(1, 0, 0);
275 ITHACAutilities::assignBC(U, 0, v);
276 //problem->nutModes.reconstruct(nut, nutCoeff, "nut");
277 phi = fvc::flux(U);
278 int iter = 0;
279 simpleControl& simple = problem->_simple();
280 std::ofstream res_os_U;
281 std::ofstream res_os_P;
282 res_os_U.open(Folder + name(counter) + "/residualsU", std::ios_base::app);
283 res_os_P.open(Folder + name(counter) + "/residualsP", std::ios_base::app);
284
285 // SIMPLE algorithm starts here
286 while ((residual_jump > residualJumpLim
287 || std::max(U_norm_res, P_norm_res) > normalizedResidualLim) &&
288 iter < maxIterOn)
289 {
290 iter++;
291 std::cout << iter << std::endl;
292#if defined(OFVER) && (OFVER == 6)
293 simple.loop(runTime);
294#else
295 problem->_simple().loop();
296#endif
297
298 // If the case is turbulent, then the network is evaluated
299 if (ITHACAutilities::isTurbulent())
300 {
301 nutCoeff = problem->evalNet(a, mu_now);
302 //nutCoeff = nutCoeffOld + 0.7*(nutCoeff - nutCoeffOld);
303 volScalarField& nut = const_cast<volScalarField&>
304 (problem->_mesh().lookupObject<volScalarField>("nut"));
305 problem->nutModes.reconstruct(nut, nutCoeff, "nut");
306 //ITHACAstream::exportSolution(nut, name(counter), Folder);
307 }
308
309 volScalarField nueff = problem->turbulence->nuEff();
310 vector v(1, 0, 0);
311 ITHACAutilities::assignBC(U, 0, v);
312 // Momentum equation
313 fvVectorMatrix UEqn
314 (
315 fvm::div(phi, U)
316 - fvm::laplacian(nueff, U)
317 - fvc::div(nueff * dev2(T(fvc::grad(U))))
318 );
319 UEqn.relax();
320 //solve(UEqn == - fvc::grad(P));
321 List<Eigen::MatrixXd> RedLinSysU = problem->Umodes.project(UEqn, NmodesUproj);
322 volVectorField gradPfull = -fvc::grad(P);
323 Eigen::MatrixXd projGrad = problem->Umodes.project(gradPfull, NmodesUproj);
324 RedLinSysU[1] = RedLinSysU[1] + projGrad;
325 a = reducedProblem::solveLinearSys(RedLinSysU, a, uresidual);
326 problem->Umodes.reconstruct(U, a, "U");
327 ITHACAutilities::assignBC(U, 0, v);
328 volScalarField rAU(1.0 / UEqn.A());
329 volVectorField HbyA(constrainHbyA(1.0 / UEqn.A() * UEqn.H(), U, P));
330 surfaceScalarField phiHbyA("phiHbyA", fvc::flux(HbyA));
331 adjustPhi(phiHbyA, U, P);
332 tmp<volScalarField> rAtU(rAU);
333
334 if (simple.consistent())
335 {
336 rAtU = 1.0 / (1.0 / rAU - UEqn.H1());
337 phiHbyA +=
338 fvc::interpolate(rAtU() - rAU) * fvc::snGrad(P) * mesh.magSf();
339 HbyA -= (rAU - rAtU()) * fvc::grad(P);
340 }
341
342 List<Eigen::MatrixXd> RedLinSysP;
343 bOld = b;
344
345 while (simple.correctNonOrthogonal())
346 {
347 // Continuity equation
348 fvScalarMatrix pEqn
349 (
350 fvm::laplacian(rAtU(), P) == fvc::div(phiHbyA)
351 );
352 RedLinSysP = problem->Pmodes.project(pEqn, NmodesPproj);
353 //pEqn.solve();
354 b = reducedProblem::solveLinearSys(RedLinSysP, b, presidual);
355 problem->Pmodes.reconstruct(P, b, "p");
356
357 if (simple.finalNonOrthogonalIter())
358 {
359 phi = phiHbyA - pEqn.flux();
360 }
361 }
362
363 b = bOld + mesh.fieldRelaxationFactor("p") * (b - bOld);
364 problem->Pmodes.reconstruct(P, b, "p");
365 nutCoeffOld = nutCoeff;
366 // P.relax();
367 U = HbyA - rAtU() * fvc::grad(P);
368 U.correctBoundaryConditions();
369 uresidualOld = uresidualOld - uresidual;
370 presidualOld = presidualOld - presidual;
371 uresidualOld = uresidualOld.cwiseAbs();
372 presidualOld = presidualOld.cwiseAbs();
373 residual_jump = std::max(uresidualOld.sum(), presidualOld.sum());
374 uresidualOld = uresidual;
375 presidualOld = presidual;
376 uresidual = uresidual.cwiseAbs();
377 presidual = presidual.cwiseAbs();
378 U_norm_res = uresidual.sum() / (RedLinSysU[1].cwiseAbs()).sum();
379 P_norm_res = presidual.sum() / (RedLinSysP[1].cwiseAbs()).sum();
380
381 if (problem->para->debug)
382 {
383 std::cout << "Residual jump = " << residual_jump << std::endl;
384 std::cout << "Normalized residual = " << std::max(U_norm_res,
385 P_norm_res) << std::endl;
386 std::cout << "Final normalized residual for velocity: " << U_norm_res <<
387 std::endl;
388 std::cout << "Final normalized residual for pressure: " << P_norm_res <<
389 std::endl;
390 }
391
392 res_os_U << U_norm_res << std::endl;
393 res_os_P << P_norm_res << std::endl;
394 }
395
396 res_os_U.close();
397 res_os_P.close();
398 std::cout << "Solution " << counter << " converged in " << iter <<
399 " iterations." << std::endl;
400 std::cout << "Final normalized residual for velocity: " << U_norm_res <<
401 std::endl;
402 std::cout << "Final normalized residual for pressure: " << P_norm_res <<
403 std::endl;
404 problem->Umodes.reconstruct(U, a, "Uaux");
405 problem->Pmodes.reconstruct(P, b, "Paux");
406
407 if (ITHACAutilities::isTurbulent())
408 {
409 volScalarField& nut = const_cast<volScalarField&>
410 (problem->_mesh().lookupObject<volScalarField>("nut"));
411 nut.rename("nutAux");
412 ITHACAstream::exportSolution(nut, name(counter), Folder);
413 }
414
415 ITHACAstream::exportSolution(U, name(counter), Folder);
416 ITHACAstream::exportSolution(P, name(counter), Folder);
417 runTime.setTime(runTime.startTime(), 0);
418 }
419};
420
421class tutorial01cl : public SteadyNSSimpleNN
422{
423 public:
425 explicit tutorial01cl(int argc, char* argv[])
426 :
427 SteadyNSSimpleNN(argc, argv)
428 {
429 curX = _mesh().points();
430 point0 = _mesh().points();
431 }
432
434 vectorField point0;
436 List<vector> curX;
437
439 void offlineSolve(Eigen::MatrixXd Box, List<label> patches,
440 word folder = "./ITHACAoutput/Offline/")
441 {
442 Vector<double> inl(0, 0, 0);
443 List<scalar> mu_now(1);
444 volVectorField& U = _U();
445 volScalarField& p = _p();
446
447 // if the offline solution is already performed read the fields
448 if (offline && ITHACAutilities::isTurbulent() &&
449 !ITHACAutilities::check_folder("./ITHACAoutput/POD/1"))
450 {
451 ITHACAstream::readMiddleFields(Ufield, U, folder);
452 ITHACAstream::readMiddleFields(Pfield, p, folder);
453 auto nut = _mesh().lookupObject<volScalarField>("nut");
454 ITHACAstream::readMiddleFields(nutFields, nut, folder);
455 }
456 else if (offline && !ITHACAutilities::check_folder("./ITHACAoutput/POD/1"))
457 {
458 ITHACAstream::readMiddleFields(Ufield, U, folder);
459 ITHACAstream::readMiddleFields(Pfield, p, folder);
460 }
461 else if (!offline)
462 {
463 //Vector<double> Uinl(1, 0, 0);
464 for (label i = 0; i < mu.rows(); i++)
465 {
466 _mesh().movePoints(point0);
467 List<vector> points2Move;
468 labelList boxIndices = ITHACAutilities::getIndicesFromBox(_mesh(), patches, Box,
469 points2Move);
470 mu_now[0] = mu(i, 0);
471 linearMovePts(mu_now[0], points2Move);
472
473 for (int j = 0; j < boxIndices.size(); j++)
474 {
475 curX[boxIndices[j]] = points2Move[j];
476 }
477
478 Field<vector> curXV(curX);
479 _mesh().movePoints(curXV);
480 ITHACAstream::writePoints(_mesh().points(), folder,
481 name(i + 1) + "/polyMesh/");
482 //assignIF(U, Uinl);
483 truthSolve2(mu_now, folder);
484 int middleF = 1;
485
486 while (ITHACAutilities::check_folder(folder + name(
487 i + 1) + "/" + name(middleF)))
488 {
489 word command("ln -s $(readlink -f " + folder + name(
490 i + 1) + "/polyMesh/) " + folder + name(i + 1) + "/" + name(
491 middleF) + "/" + " >/dev/null 2>&1");
492 std::cout.setstate(std::ios_base::failbit);
493 system(command);
494 std::cout.clear();
495 middleF++;
496 }
497
498 restart();
499 }
500 }
501 }
502
503 void linearMovePts(double angle, List<vector>& points2Move)
504 {
505 double sMax;
506 scalarList x;
507 scalarList y;
508
509 for (label i = 0; i < points2Move.size(); i++)
510 {
511 x.append(points2Move[i].component(0));
512 y.append(points2Move[i].component(1));
513 }
514
515 double xMin = min(x);
516 double xMax = max(x);
517 double yMin = min(y);
518 double yMax = max(y);
519 sMax = (yMax - yMin) * std::tan(M_PI * angle / 180);
520
521 for (label i = 0; i < points2Move.size(); i++)
522 {
523 points2Move[i].component(0) = points2Move[i].component(0) +
524 (yMax - points2Move[i].component(1)) / (yMax - yMin) * (xMax -
525 points2Move[i].component(0)) / (xMax - xMin) * sMax;
526 }
527 }
528};
529
530int main(int argc, char* argv[])
531{
532 // Construct the tutorial object
533 tutorial01cl example(argc, argv);
534 // Read some parameters from file
535 ITHACAparameters* para = ITHACAparameters::getInstance(example._mesh(),
536 example._runTime());
537 // Read the files where the parameters are stored
538 std::ifstream exFileOff("./angOff_mat.txt");
539
540 if (exFileOff)
541 {
542 example.mu = ITHACAstream::readMatrix("./angOff_mat.txt");
543 }
544 else
545 {
546 example.mu = Eigen::VectorXd::LinSpaced(50, 0, 75);
547 ITHACAstream::exportMatrix(example.mu, "angOff", "eigen", "./");
548 }
549
550 Eigen::MatrixXd angOn;
551 std::ifstream exFileOn("./angOn_mat.txt");
552
553 if (exFileOn)
554 {
555 angOn = ITHACAstream::readMatrix("./angOn_mat.txt");
556 }
557 else
558 {
559 //angOn = Eigen::VectorXd::LinSpaced(20, 5, 70);
560 angOn = ITHACAutilities::rand(10, 1, 0, 70);
561 ITHACAstream::exportMatrix(angOn, "angOn", "eigen", "./");
562 }
563
564 //Set the box including the step
565 Eigen::MatrixXd Box(2, 3);
566 Box << 1.98, 0.01, 0.11,
567 7.02, -0.666669, -0.01;
568 //Select the patches to be moved
569 List<label> movPat;
570 movPat.append(3);
571 movPat.append(5);
572 // Set the maximum iterations number for the offline phase
573 example.maxIter = para->ITHACAdict->lookupOrDefault<int>("maxIter", 2000);
574 // Perform the offline solve
575 example.offlineSolve(Box, movPat);
576 List<vector> points2Move;
577 labelList boxIndices = ITHACAutilities::getIndicesFromBox(example._mesh(),
578 movPat, Box,
579 points2Move);
580 example.linearMovePts((example.mu.maxCoeff() + example.mu.minCoeff()) / 2,
581 points2Move);
582
583 for (int j = 0; j < boxIndices.size(); j++)
584 {
585 example.curX[boxIndices[j]] = points2Move[j];
586 }
587
588 Field<vector> curXV(example.curX);
589 example._mesh().movePoints(curXV);
590 // Perform POD on velocity and pressure and store the first 10 modes
591 ITHACAPOD::getModes(example.Ufield, example.Umodes, example._U().name(),
592 example.podex, 0, 0,
593 example.NUmodesOut);
594 ITHACAPOD::getModes(example.Pfield, example.Pmodes, example._p().name(),
595 example.podex, 0, 0,
596 example.NPmodesOut);
597 // Error analysis
598 tutorial01cl checkOff(argc, argv);
599 std::clock_t startOff;
600 double durationOff;
601 startOff = std::clock();
602
603 if (!ITHACAutilities::check_folder("./ITHACAoutput/checkOff"))
604 {
605 checkOff.restart();
606 ITHACAparameters* para = ITHACAparameters::getInstance(checkOff._mesh(),
607 checkOff._runTime());
608 checkOff.offline = false;
609 checkOff.mu = angOn;
610 checkOff.offlineSolve(Box, movPat, "./ITHACAoutput/checkOff/");
611 checkOff.offline = true;
612 }
613
614 durationOff = (std::clock() - startOff);
615
616 if (ITHACAutilities::isTurbulent())
617 {
618 ITHACAPOD::getModes(example.nutFields, example.nutModes, "nut",
619 example.podex, 0, 0, example.NNutModesOut);
620 // Create the RBF for turbulence
621 example.getTurbNN();
622 }
623
624 std::string wrng =
625 "The network for this problem has not been calculated yet. Perform the Python training (see train.py).";
626 M_Assert(ITHACAutilities::check_file("ITHACAoutput/NN/Net_" + name(
627 example.NUmodes) + "_" + name(
628 example.NNutModes) + ".pt"), wrng.c_str());
629 example.loadNet("ITHACAoutput/NN/Net_" + name(example.NUmodes) + "_" + name(
630 example.NNutModes) + ".pt");
631 example._mesh().movePoints(example.point0);
632 // Create the reduced object
633 reducedSimpleSteadyNN reduced(example);
634 PtrList<volVectorField> U_rec_list;
635 PtrList<volScalarField> P_rec_list;
636 // Reads inlet volocities boundary conditions.
637 word vel_file(para->ITHACAdict->lookup("online_velocities"));
638 Eigen::MatrixXd vel = ITHACAstream::readMatrix(vel_file);
639 // Set the maximum iterations number for the online phase
640 reduced.maxIterOn = para->ITHACAdict->lookupOrDefault<int>("maxIterOn", 2000);
641 std::cout << "Total amout of parameters to be solved online: " + name(
642 angOn.rows()) << std::endl;
643 //Perform the online solutions
644 std::clock_t startOn;
645 double durationOn;
646 startOn = std::clock();
647
648 for (label k = 0; k < angOn.rows(); k++)
649 {
650 std::cout << "Solving online for parameter number " + name(k + 1) << std::endl;
651 scalar mu_now = angOn(k, 0);
652 example.restart();
653 reduced.vel_now = vel;
654
655 if (ITHACAutilities::isTurbulent())
656 {
657 example._mesh().movePoints(example.point0);
658 List<vector> points2Move;
659 labelList boxIndices = ITHACAutilities::getIndicesFromBox(example._mesh(),
660 movPat, Box, points2Move);
661 example.linearMovePts(mu_now, points2Move);
662
663 for (int j = 0; j < boxIndices.size(); j++)
664 {
665 example.curX[boxIndices[j]] = points2Move[j];
666 }
667
668 Field<vector> curXV(example.curX);
669 example._mesh().movePoints(curXV);
670 ITHACAstream::writePoints(example._mesh().points(),
671 "./ITHACAoutput/Reconstruct_" + name(example.NUmodes) + "_" + name(
672 example.NPmodes) + "/", name(k + 1) + "/polyMesh/");
673 reduced.solveOnline_Simple(mu_now, example.NUmodes, example.NPmodes,
674 example.NNutModes, "./ITHACAoutput/Reconstruct_" + name(
675 example.NUmodes) + "_" + name(example.NPmodes) + "/");
676 }
677 else
678 {
679 example._mesh().movePoints(example.point0);
680 List<vector> points2Move;
681 labelList boxIndices = ITHACAutilities::getIndicesFromBox(example._mesh(),
682 movPat, Box, points2Move);
683 example.linearMovePts(mu_now, points2Move);
684
685 for (int j = 0; j < boxIndices.size(); j++)
686 {
687 example.curX[boxIndices[j]] = points2Move[j];
688 }
689
690 Field<vector> curXV(example.curX);
691 example._mesh().movePoints(curXV);
692 ITHACAstream::writePoints(example._mesh().points(),
693 "./ITHACAoutput/Reconstruct_" + name(example.NUmodes) + "_" + name(
694 example.NPmodes) + "/", name(k + 1) + "/polyMesh/");
695 reduced.solveOnline_Simple(mu_now, example.NUmodes, example.NPmodes,
696 example.NNutModes, "./ITHACAoutput/Reconstruct_" + name(
697 example.NUmodes) + "_" + name(example.NPmodes) + "/");
698 }
699 }
700
701 durationOn = (std::clock() - startOn);
702 std::cout << "ONLINE PHASE COMPLETED" << std::endl;
703 PtrList<volVectorField> Ufull;
704 PtrList<volScalarField> Pfull;
705 PtrList<volVectorField> Ured;
706 PtrList<volScalarField> Pred;
707 volVectorField Uaux("Uaux", checkOff._U());
708 volScalarField Paux("Paux", checkOff._p());
709 ITHACAstream::readConvergedFields(Ufull, checkOff._U(),
710 "./ITHACAoutput/checkOff/");
711 ITHACAstream::readConvergedFields(Pfull, checkOff._p(),
712 "./ITHACAoutput/checkOff/");
713 ITHACAstream::read_fields(Ured, Uaux,
714 "./ITHACAoutput/Reconstruct_" + name(example.NUmodes) + "_" + name(
715 example.NPmodes) + "/");
716 ITHACAstream::read_fields(Pred, Paux,
717 "./ITHACAoutput/Reconstruct_" + name(example.NUmodes) + "_" + name(
718 example.NPmodes) + "/");
719 Eigen::MatrixXd relErrorU(Ufull.size(), 1);
720 Eigen::MatrixXd relErrorP(Pfull.size(), 1);
721 dimensionedVector U_fs("U_fs", dimVelocity, vector(1, 0, 0));
722
723 for (label k = 0; k < Ufull.size(); k++)
724 {
725 volVectorField errorU = (Ufull[k] - Ured[k]).ref();
726 volVectorField devU = (Ufull[k] - U_fs).ref();
727 volScalarField errorP = (Pfull[k] - Pred[k]).ref();
728 relErrorU(k, 0) = ITHACAutilities::frobNorm(errorU) /
729 ITHACAutilities::frobNorm(devU);
730 relErrorP(k, 0) = ITHACAutilities::frobNorm(errorP) /
731 ITHACAutilities::frobNorm(Pfull[k]);
732 }
733
734 ITHACAstream::exportMatrix(relErrorU,
735 "errorU_" + name(example.NUmodes) + "_" + name(example.NPmodes) + "_" + name(
736 example.NNutModes), "python", ".");
737 ITHACAstream::exportMatrix(relErrorP,
738 "errorP_" + name(example.NUmodes) + "_" + name(example.NPmodes) + "_" + name(
739 example.NNutModes), "python", ".");
740 std::cout << "The online phase duration is equal to " << durationOn <<
741 std::endl;
742 std::cout << "The offline phase duration is equal to " << durationOff <<
743 std::endl;
744 exit(0);
745}
main
int main(int argc, char *argv[])
Definition 01simpleTurbGeomClosed.C:530
pEqn
fvScalarMatrix pEqn
Definition CFM.H:31
ITHACAPOD.H
Header file of the ITHACAPOD class.
mesh
Foam::fvMesh & mesh
Definition createMesh.H:47
runTime
Foam::Time & runTime
Definition createTime.H:33
M_Assert
#define M_Assert(Expr, Msg)
Definition ITHACAassert.H:46
ITHACAstream.H
Header file of the ITHACAstream class, it contains the implementation of several methods for input ou...
uresidual
scalar uresidual
Definition NLmsrProblem.H:27
presidual
scalar presidual
Definition NLmsrProblem.H:29
ReducedSimpleSteadyNS.H
Header file of the reducedSteadyNS class.
SteadyNSSimple.H
Header file of the steadyNS class.
ITHACAparameters
Class for the definition of some general parameters, the parameters must be defined from the file ITH...
Definition ITHACAparameters.H:41
ITHACAparameters::ITHACAdict
IOdictionary * ITHACAdict
Dictionary for input objects from file.
Definition ITHACAparameters.H:64
ITHACAparameters::getInstance
static ITHACAparameters * getInstance()
Gets an instance of ITHACAparameters, to be used if the instance is already existing.
Definition ITHACAparameters.C:55
Modes::reconstruct
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:275
Modes::project
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
SteadyNSSimpleNN::scale_inp
Eigen::MatrixXd scale_inp
Definition 01simpleTurbGeomClosed.C:62
SteadyNSSimpleNN::Net
torch::nn::Sequential Net
Definition 01simpleTurbGeomClosed.C:72
SteadyNSSimpleNN::coeffL2U_tensor
torch::Tensor coeffL2U_tensor
Definition 01simpleTurbGeomClosed.C:69
SteadyNSSimpleNN::coeffL2Nut
Eigen::MatrixXd coeffL2Nut
Definition 01simpleTurbGeomClosed.C:66
SteadyNSSimpleNN::SteadyNSSimpleNN
SteadyNSSimpleNN(int argc, char *argv[])
Constructor.
Definition 01simpleTurbGeomClosed.C:48
SteadyNSSimpleNN::optimizer
torch::optim::Optimizer * optimizer
Definition 01simpleTurbGeomClosed.C:73
SteadyNSSimpleNN::scale_out
Eigen::MatrixXd scale_out
Definition 01simpleTurbGeomClosed.C:64
SteadyNSSimpleNN::coeffL2Nut_tensor
torch::Tensor coeffL2Nut_tensor
Definition 01simpleTurbGeomClosed.C:70
SteadyNSSimpleNN::netTorchscript
torch::jit::script::Module netTorchscript
Definition 01simpleTurbGeomClosed.C:74
SteadyNSSimpleNN::loadNet
void loadNet(word filename)
Definition 01simpleTurbGeomClosed.C:76
SteadyNSSimpleNN::evalNet
Eigen::MatrixXd evalNet(Eigen::MatrixXd a, double mu_now)
Definition 01simpleTurbGeomClosed.C:160
SteadyNSSimpleNN::evalNet
Eigen::MatrixXd evalNet(Eigen::MatrixXd a)
Definition 01simpleTurbGeomClosed.C:185
SteadyNSSimple
Implementation of a parametrized full order steady NS problem and preparation of the the reduced ma...
Definition SteadyNSSimple.H:70
SteadyNSSimple::truthSolve2
void truthSolve2(List< scalar > mu_now, word Folder="./ITHACAoutput/Offline/")
Offline solver for the whole Navier-Stokes problem.
Definition SteadyNSSimple.C:114
SteadyNSSimple::nutModes
volScalarModes nutModes
List of POD modes for eddy viscosity.
Definition SteadyNSSimple.H:85
SteadyNSSimple::nutFields
PtrList< volScalarField > nutFields
List of snapshots for the solution for eddy viscosity.
Definition SteadyNSSimple.H:82
reducedProblem::solveLinearSys
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.
Definition ReducedProblem.C:64
reducedSimpleSteadyNN::reducedSimpleSteadyNN
reducedSimpleSteadyNN(SteadyNSSimpleNN &FOMproblem)
Constructor.
Definition 01simpleTurbGeomClosed.C:206
reducedSimpleSteadyNN::problem
SteadyNSSimpleNN * problem
Full problem.
Definition 01simpleTurbGeomClosed.C:212
reducedSimpleSteadyNN::solveOnline_Simple
void solveOnline_Simple(scalar mu_now, int NmodesUproj, int NmodesPproj, int NmodesNut=0, word Folder="./ITHACAoutput/Reconstruct/")
Definition 01simpleTurbGeomClosed.C:215
reducedSimpleSteadyNS
Class where it is implemented a reduced problem for the steady Navier-stokes problem.
Definition ReducedSimpleSteadyNS.H:60
reducedSimpleSteadyNS::maxIterOn
int maxIterOn
Maximum iterations number for the online step.
Definition ReducedSimpleSteadyNS.H:118
reducedSimpleSteadyNS::vel_now
Eigen::MatrixXd vel_now
Imposed boundary conditions.
Definition ReducedSimpleSteadyNS.H:115
reductionProblem::offline
bool offline
Boolean variable, it is 1 if the Offline phase has already been computed, else 0.
Definition reductionProblem.H:99
reductionProblem::mu
Eigen::MatrixXd mu
Row matrix of parameters.
Definition reductionProblem.H:89
reductionProblem::podex
bool podex
Boolean variable, it is 1 if the POD has already been computed, else 0.
Definition reductionProblem.H:97
steadyNS::maxIter
scalar maxIter
Number of maximum iterations to be done for the computation of the truth solution.
Definition steadyNS.H:125
steadyNS::restart
void restart()
set U and P back to the values into the 0 folder
Definition steadyNS.C:2253
steadyNS::NPmodes
label NPmodes
Number of pressure modes used for the projection.
Definition steadyNS.H:143
steadyNS::_phi
autoPtr< surfaceScalarField > _phi
Flux.
Definition steadyNS.H:293
steadyNS::_simple
autoPtr< simpleControl > _simple
simpleControl
Definition steadyNS.H:284
steadyNS::Pfield
PtrList< volScalarField > Pfield
List of pointers used to form the pressure snapshots matrix.
Definition steadyNS.H:86
steadyNS::_runTime
autoPtr< Time > _runTime
Time.
Definition steadyNS.H:290
steadyNS::Umodes
volVectorModes Umodes
List of pointers used to form the velocity modes.
Definition steadyNS.H:101
steadyNS::Ufield
PtrList< volVectorField > Ufield
List of pointers used to form the velocity snapshots matrix.
Definition steadyNS.H:89
steadyNS::_mesh
autoPtr< fvMesh > _mesh
Mesh.
Definition steadyNS.H:281
steadyNS::NUmodes
label NUmodes
Number of velocity modes used for the projection.
Definition steadyNS.H:140
steadyNS::para
ITHACAparameters * para
Definition steadyNS.H:82
steadyNS::NNutModesOut
label NNutModesOut
Number of nut modes to be calculated.
Definition steadyNS.H:137
steadyNS::NNutModes
label NNutModes
Number of nut modes used for the projection.
Definition steadyNS.H:149
steadyNS::turbulence
autoPtr< incompressible::turbulenceModel > turbulence
Turbulence model.
Definition steadyNS.H:299
steadyNS::NUmodesOut
label NUmodesOut
Number of velocity modes to be calculated.
Definition steadyNS.H:128
steadyNS::NPmodesOut
label NPmodesOut
Number of pressure modes to be calculated.
Definition steadyNS.H:131
steadyNS::_U
autoPtr< volVectorField > _U
Velocity field.
Definition steadyNS.H:263
steadyNS::Pmodes
volScalarModes Pmodes
List of pointers used to form the pressure modes.
Definition steadyNS.H:98
steadyNS::_p
autoPtr< volScalarField > _p
Pressure field.
Definition steadyNS.H:260
tutorial01cl::curX
List< vector > curX
List to store the moved coordinates.
Definition 01simpleTurbGeomClosed.C:436
tutorial01cl::linearMovePts
void linearMovePts(double angle, List< vector > &points2Move)
Definition 01simpleTurbGeomClosed.C:503
tutorial01cl::point0
vectorField point0
Initial coordinates of the grid points.
Definition 01simpleTurbGeomClosed.C:434
tutorial01cl::tutorial01cl
tutorial01cl(int argc, char *argv[])
Constructor.
Definition 01simpleTurbGeomClosed.C:425
tutorial01cl::offlineSolve
void offlineSolve(Eigen::MatrixXd Box, List< label > patches, word folder="./ITHACAoutput/Offline/")
Perform an Offline solve.
Definition 01simpleTurbGeomClosed.C:439
T
volScalarField & T
Definition createT.H:46
v
volScalarField & v
Definition createConstants.H:69
ITHACAPOD::getModes
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
ITHACAstream::writePoints
void writePoints(pointField points, fileName folder, fileName subfolder)
Write points of a mesh to a file.
Definition ITHACAstream.C:942
ITHACAstream::readConvergedFields
void readConvergedFields(PtrList< GeometricField< Type, PatchField, GeoMesh > > &Lfield, GeometricField< Type, PatchField, GeoMesh > &field, fileName casename)
Function to read a list of volVectorField from name of the field including only converged snapshots.
Definition ITHACAstream.C:659
ITHACAstream::exportSolution
void exportSolution(GeometricField< Type, PatchField, GeoMesh > &s, fileName subfolder, fileName folder, word fieldName)
Export a field to file in a certain folder and subfolder.
Definition ITHACAstream.C:849
ITHACAstream::readMiddleFields
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...
Definition ITHACAstream.C:643
ITHACAstream::exportMatrix
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 ...
Definition ITHACAstream.C:55
ITHACAstream::readMatrix
List< Eigen::MatrixXd > readMatrix(word folder, word mat_name)
Read a three dimensional matrix from a txt file in Eigen format.
Definition ITHACAstream.C:372
ITHACAstream::read_fields
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.
Definition ITHACAstream.C:454
ITHACAtorch::torch2Eigen::eigenMatrix2torchTensor
torch::Tensor eigenMatrix2torchTensor(Eigen::Matrix< type, Eigen::Dynamic, Eigen::Dynamic > _eigenMatrix)
Convert an eigen Matrix to a torch tensor.
Definition torch2Eigen.C:39
ITHACAtorch::torch2Eigen::torchTensor2eigenMatrix< double >
template Eigen::Matrix< double, Eigen::Dynamic, Eigen::Dynamic > torchTensor2eigenMatrix< double >(torch::Tensor &torchTensor)
ITHACAutilities::getCoeffs
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.
Definition ITHACAcoeffsMass.C:317
ITHACAutilities::frobNorm
double frobNorm(GeometricField< Type, PatchField, GeoMesh > &field)
Evaluate the Frobenius norm of a field.
Definition ITHACAerror.C:482
ITHACAutilities::rand
Eigen::MatrixXd rand(label rows, label cols, double min, double max)
Generates random matrix with random values in an interval.
Definition ITHACAutilities.C:44
ITHACAutilities::isTurbulent
bool isTurbulent()
This function checks if the case is turbulent.
Definition ITHACAutilities.C:90
ITHACAutilities::assignBC
void assignBC(GeometricField< scalar, fvPatchField, volMesh > &s, label BC_ind, double value)
Assign uniform Boundary Condition to a volScalarField.
Definition ITHACAassign.C:143
ITHACAutilities::check_folder
bool check_folder(word folder)
Checks if a folder exists.
Definition ITHACAsystem.C:89
ITHACAutilities::check_file
bool check_file(std::string fileName)
Function that returns true if a file exists.
Definition ITHACAsystem.C:136
ITHACAutilities::getIndicesFromBox
labelList getIndicesFromBox(const fvMesh &mesh, List< label > indices, Eigen::MatrixXd Box, List< vector > &points2Move)
Gives the indices conteined into a defined box.
Definition ITHACAgeometry.C:31
cnpy::save
void save(const Eigen::Matrix< T, -1, dim > &mat, const std::string fname)
cnpy::load
Eigen::Matrix< T, -1, dim > load(Eigen::Matrix< T, -1, dim > &mat, const std::string fname, std::string order="rowMajor")
simple
simpleControl simple(mesh)
phi
surfaceScalarField & phi
Definition createFields.H:68
U
volVectorField & U
Definition createFields.H:42
nut
volScalarField & nut
Definition createFields.H:72
g
dimensionedVector & g
Definition createFields.H:150
UEqn
fvVectorMatrix & UEqn
Definition UEqn.H:37
p
volScalarField & p
Definition createFields.H:41
adjustPhi
adjustPhi(phiHbyA, U, p)
rAtU
tmp< volScalarField > rAtU(rAU)
phiHbyA
surfaceScalarField phiHbyA("phiHbyA", fvc::flux(HbyA))
HbyA
volVectorField HbyA(constrainHbyA(rAU *UEqn.H(), U, p))
i
label i
Definition pEqn.H:46
rAU
volScalarField rAU(1.0/UEqn.A())
torch2Eigen.H
Header file of the torch2Eigen namespace.
Generated by   doxygen 1.11.0