Loading...
Searching...
No Matches
01simpleTurbGeomClosed.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 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 Info << "Computing the coefficients for U train" << endl;
112 Eigen::MatrixXd coeffL2U_train = ITHACAutilities::getCoeffs(UfieldTrain,
113 Umodes,
114 0, true);
115 Info << "Computing the coefficients for p train" << endl;
116 Eigen::MatrixXd coeffL2P_train = ITHACAutilities::getCoeffs(PfieldTrain,
117 Pmodes,
118 0, true);
119 Info << "Computing the coefficients for nuT train" << 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 Info << iter << 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 Info << "Residual jump = " << residual_jump << endl;
384 Info << "Normalized residual = " << std::max(U_norm_res,
385 P_norm_res) << endl;
386 Info << "Final normalized residual for velocity: " << U_norm_res <<
387 endl;
388 Info << "Final normalized residual for pressure: " << P_norm_res <<
389 endl;
390 }
391
392 res_os_U << U_norm_res << endl;
393 res_os_P << P_norm_res << endl;
394 }
395
396 res_os_U.close();
397 res_os_P.close();
398 Info << "Solution " << counter << " converged in " << iter <<
399 " iterations." << endl;
400 Info << "Final normalized residual for velocity: " << U_norm_res <<
401 endl;
402 Info << "Final normalized residual for pressure: " << P_norm_res <<
403 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 Info << "Total amout of parameters to be solved online: " + name(
642 angOn.rows()) << 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 Info << "Solving online for parameter number " + name(k + 1) << 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 Info << "ONLINE PHASE COMPLETED" << 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 Info << "The online phase duration is equal to " << durationOn <<
741 endl;
742 Info << "The offline phase duration is equal to " << durationOff <<
743 endl;
744 return 0;
745}
ITHACAPOD.H
Header file of the ITHACAPOD class.
ITHACAstream.H
Header file of the ITHACAstream class, it contains the implementation of several methods for input ou...
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::getInstance
static ITHACAparameters * getInstance()
Gets an instance of ITHACAparameters, to be used if the instance is already existing.
Definition ITHACAparameters.C:56
SteadyNSSimpleNN::SteadyNSSimpleNN
SteadyNSSimpleNN(int argc, char *argv[])
Constructor.
Definition 01simpleTurbGeomClosed.C:48
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
reducedSimpleSteadyNS::maxIterOn
int maxIterOn
Maximum iterations number for the online step.
Definition ReducedSimpleSteadyNS.H:118
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
steadyNS::restart
void restart()
set U and P back to the values into the 0 folder
Definition steadyNS.C:2305
steadyNS::Pfield
PtrList< volScalarField > Pfield
List of pointers used to form the pressure snapshots matrix.
Definition steadyNS.H:86
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::NNutModes
label NNutModes
Number of nut modes used for the projection.
Definition steadyNS.H:152
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::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
EigenFunctions::max
T max(Eigen::SparseMatrix< T > &mat, label &ind_row, label &ind_col)
Find the maximum of a sparse Matrix (Useful for DEIM).
Definition EigenFunctions.H:300
EigenFunctions::min
T min(Eigen::SparseMatrix< T > &mat, label &ind_row, label &ind_col)
Find the minimum of a sparse Matrix (Useful for DEIM).
Definition EigenFunctions.H:324
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:1092
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:801
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:992
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:785
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:517
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:349
ITHACAutilities::frobNorm
double frobNorm(GeometricField< Type, PatchField, GeoMesh > &field)
Evaluate the Frobenius norm of a field.
Definition ITHACAnorm.C:355
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:80
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
torch2Eigen.H
Header file of the torch2Eigen namespace.
Generated by   doxygen 1.16.1