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
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
startOff
std::clock_t startOff
Definition NLsolve.H:47
para
ITHACAparameters * para
Definition NLsolve.H:40
nut
auto nut
Definition NLsolve.H:195
durationOff
double durationOff
Definition NLsolve.H:49
presidual
scalar presidual
Definition NLsolve.H:38
ReducedSimpleSteadyNS.H
Header file of the reducedSteadyNS class.
SteadyNSSimple.H
Header file of the steadyNS class.
pEqn
fvScalarMatrix pEqn
Definition CFM.H:31
ITHACAparameters
Class for the definition of some general parameters, the parameters must be defined from the file ITH...
Definition ITHACAparameters.H:41
ITHACAparameters::getInstance
static ITHACAparameters * getInstance()
Gets an instance of ITHACAparameters, to be used if the instance is already existing.
Definition ITHACAparameters.C:56
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::SteadyNSSimple
SteadyNSSimple()
Null constructor.
Definition SteadyNSSimple.C:41
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::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
reducedSimpleSteadyNS::reducedSimpleSteadyNS
reducedSimpleSteadyNS()
Construct Null.
Definition ReducedSimpleSteadyNS.C:39
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:128
steadyNS::restart
void restart()
set U and P back to the values into the 0 folder
Definition steadyNS.C:2277
steadyNS::NPmodes
label NPmodes
Number of pressure modes used for the projection.
Definition steadyNS.H:146
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:299
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:290
steadyNS::NUmodes
label NUmodes
Number of velocity modes used for the projection.
Definition steadyNS.H:143
steadyNS::NNutModesOut
label NNutModesOut
Number of nut modes to be calculated.
Definition steadyNS.H:140
steadyNS::NNutModes
label NNutModes
Number of nut modes used for the projection.
Definition steadyNS.H:152
steadyNS::NUmodesOut
label NUmodesOut
Number of velocity modes to be calculated.
Definition steadyNS.H:131
steadyNS::NPmodesOut
label NPmodesOut
Number of pressure modes to be calculated.
Definition steadyNS.H:134
steadyNS::_U
autoPtr< volVectorField > _U
Velocity field.
Definition steadyNS.H:272
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:269
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:90
ITHACAstream::writePoints
void writePoints(pointField points, fileName folder, fileName subfolder)
Write points of a mesh to a file.
Definition ITHACAstream.C:997
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:717
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:904
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:701
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:516
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:325
ITHACAutilities::frobNorm
double frobNorm(GeometricField< Type, PatchField, GeoMesh > &field)
Evaluate the Frobenius norm of a field.
Definition ITHACAerror.C:477
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)
rAtU
volScalarField rAtU(1.0/(1.0/rAU - Ueqn.H1()))
rAU
volScalarField rAU(1.0/Ueqn.A())
phiHbyA
phiHbyA
Definition pcEqn.H:61
HbyA
HbyA
Definition pcEqn.H:62
simple
simpleControl simple(mesh)
phi
surfaceScalarField & phi
Definition createFields.H:68
U
volVectorField & U
Definition createFields.H:42
UEqn
fvVectorMatrix & UEqn
Definition UEqn.H:13
p
volScalarField & p
Definition createFields.H:42
g
dimensionedVector & g
Definition createFields.H:150
adjustPhi
adjustPhi(phiHbyA, U, p)
i
label i
Definition pEqn.H:46
torch2Eigen.H
Header file of the torch2Eigen namespace.
Generated by   doxygen 1.13.2