Loading...
Searching...
No Matches
03compSteadyNS.C
Go to the documentation of this file.
1/*---------------------------------------------------------------------------*\
2 ██╗████████╗██╗ ██╗ █████╗ ██████╗ █████╗ ███████╗██╗ ██╗
3 ██║╚══██╔══╝██║ ██║██╔══██╗██╔════╝██╔══██╗ ██╔════╝██║ ██║
4 ██║ ██║ ███████║███████║██║ ███████║█████╗█████╗ ██║ ██║
5 ██║ ██║ ██╔══██║██╔══██║██║ ██╔══██║╚════╝██╔══╝ ╚██╗ ██╔╝
6 ██║ ██║ ██║ ██║██║ ██║╚██████╗██║ ██║ ██║ ╚████╔╝
7 ╚═╝ ╚═╝ ╚═╝ ╚═╝╚═╝ ╚═╝ ╚═════╝╚═╝ ╚═╝ ╚═╝ ╚═══╝
8
9 * In real Time Highly Advanced Computational Applications for Finite Volumes
10 * Copyright (C) 2017 by the ITHACA-FV authors
11-------------------------------------------------------------------------------
12License
13 This file is part of ITHACA-FV
14 ITHACA-FV is free software: you can redistribute it and/or modify
15 it under the terms of the GNU Lesser General Public License as published by
16 the Free Software Foundation, either version 3 of the License, or
17 (at your option) any later version.
18 ITHACA-FV is distributed in the hope that it will be useful,
19 but WITHOUT ANY WARRANTY; without even the implied warranty of
20 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
21 GNU Lesser General Public License for more details.
22 You should have received a copy of the GNU Lesser General Public License
23 along with ITHACA-FV. If not, see <http://www.gnu.org/licenses/>.
24Description
25 Example of steady NS Reduction Problem
26SourceFiles
27 03steadyNS.C
28\*---------------------------------------------------------------------------*/
29
30#include <torch/script.h>
31#include <torch/torch.h>
32#include "torch2Eigen.H"
33#include "CompressibleSteadyNS.H"
34#include "ReducedCompressibleSteadyNS.H"
35#include "RBFMotionSolver.H"
36#include "ITHACAstream.H"
37#include "ITHACAPOD.H"
38#include "forces.H"
39#include "IOmanip.H"
40
41using namespace ITHACAtorch::torch2Eigen;
42
43class CompressibleSteadyNN : public CompressibleSteadyNS
44{
45 public:
47 CompressibleSteadyNN(int argc, char* argv[])
48 :
49 CompressibleSteadyNS(argc, argv)
50 {
51 ITHACAparameters* para = ITHACAparameters::getInstance();
52 NUmodes = para->ITHACAdict->lookupOrDefault<label>("NmodesUproj", 10);
53 NNutModes = para->ITHACAdict->lookupOrDefault<label>("NmodesNutProj", 10);
54 // Net->push_back(torch::nn::Linear(NUmodes, 128));
55 // Net->push_back(torch::nn::ReLU());
56 // Net->push_back(torch::nn::Linear(128, 64));
57 // Net->push_back(torch::nn::ReLU());
58 // Net->push_back(torch::nn::Linear(64, NNutModes));
59 // optimizer = new torch::optim::Adam(Net->parameters(),
60 // torch::optim::AdamOptions(2e-2));
61 };
62
63 label NUmodes;
64 label NNutModes;
65
67 // Eddy viscosity initialization //
69
70 Eigen::MatrixXd bias_inp;
71 Eigen::MatrixXd scale_inp;
72 Eigen::MatrixXd bias_out;
73 Eigen::MatrixXd scale_out;
74
75 Eigen::MatrixXd coeffL2Nut;
76 Eigen::MatrixXd coeffL2U;
77
78 torch::Tensor coeffL2U_tensor;
79 torch::Tensor coeffL2Nut_tensor;
80
81 torch::nn::Sequential Net;
82 torch::optim::Optimizer* optimizer;
83 torch::jit::script::Module netTorchscript;
84
85 void loadNet(word filename)
86 {
87 std::string Msg = filename +
88 " is not existing, please run the training stage of the net with the correct number of modes for U and Nut";
89 M_Assert(ITHACAutilities::check_file(filename), Msg.c_str());
90 netTorchscript = torch::jit::load(filename);
91 cnpy::load(bias_inp, "ITHACAoutput/NN/minAnglesInp_" + name(
92 NUmodes) + "_" + name(NNutModes) + ".npy");
93 cnpy::load(scale_inp, "ITHACAoutput/NN/scaleAnglesInp_" + name(
94 NUmodes) + "_" + name(NNutModes) + ".npy");
95 cnpy::load(bias_out, "ITHACAoutput/NN/minOut_" + name(NUmodes) + "_" + name(
96 NNutModes) + ".npy");
97 cnpy::load(scale_out, "ITHACAoutput/NN/scaleOut_" + name(NUmodes) + "_" + name(
98 NNutModes) + ".npy");
99 netTorchscript.eval();
100 }
101
102 // This function computes the coefficients which are later used for training
103 void getTurbNN()
104 {
105 if (!ITHACAutilities::check_folder("ITHACAoutput/NN/coeffs"))
106 {
107 mkDir("ITHACAoutput/NN/coeffs");
108 // Read Fields for Train
109 PtrList<volVectorField> UfieldTrain;
110 PtrList<volScalarField> PfieldTrain;
111 PtrList<volScalarField> nutFieldsTrain;
112 ITHACAstream::readMiddleFields(UfieldTrain, _U(),
113 "./ITHACAoutput/Offline/");
114 ITHACAstream::readMiddleFields(PfieldTrain, _p(),
115 "./ITHACAoutput/Offline/");
116 auto nut_train = _mesh().lookupObject<volScalarField>("nut");
117 ITHACAstream::readMiddleFields(nutFieldsTrain, nut_train,
118 "./ITHACAoutput/Offline/");
120 std::cout << "Computing the coefficients for U train" << std::endl;
121 Eigen::MatrixXd coeffL2U_train = ITHACAutilities::getCoeffs(UfieldTrain,
122 Umodes,
123 0, true);
124 std::cout << "Computing the coefficients for p train" << std::endl;
125 Eigen::MatrixXd coeffL2P_train = ITHACAutilities::getCoeffs(PfieldTrain,
126 Pmodes,
127 0, true);
128 std::cout << "Computing the coefficients for nuT train" << std::endl;
129 Eigen::MatrixXd coeffL2Nut_train = ITHACAutilities::getCoeffs(nutFieldsTrain,
130 nutModes,
131 0, true);
132 coeffL2U_train.transposeInPlace();
133 coeffL2P_train.transposeInPlace();
134 coeffL2Nut_train.transposeInPlace();
135 cnpy::save(coeffL2U_train, "ITHACAoutput/NN/coeffs/coeffL2UTrain.npy");
136 cnpy::save(coeffL2P_train, "ITHACAoutput/NN/coeffs/coeffL2PTrain.npy");
137 cnpy::save(coeffL2Nut_train, "ITHACAoutput/NN/coeffs/coeffL2NutTrain.npy");
138 // Read Fields for Test
139 PtrList<volVectorField> UfieldTest;
140 PtrList<volScalarField> PfieldTest;
141 PtrList<volScalarField> nutFieldsTest;
143 ITHACAstream::readMiddleFields(UfieldTest, _U(),
144 "./ITHACAoutput/checkOff/");
145 ITHACAstream::readMiddleFields(PfieldTest, _p(),
146 "./ITHACAoutput/checkOff/");
147 auto nut_test = _mesh().lookupObject<volScalarField>("nut");
148 ITHACAstream::readMiddleFields(nutFieldsTest, nut_test,
149 "./ITHACAoutput/checkOff/");
150 // Compute the coefficients for test
151 Eigen::MatrixXd coeffL2U_test = ITHACAutilities::getCoeffs(UfieldTest,
152 Umodes,
153 0, true);
154 Eigen::MatrixXd coeffL2P_test = ITHACAutilities::getCoeffs(PfieldTest,
155 Pmodes,
156 0, true);
157 Eigen::MatrixXd coeffL2Nut_test = ITHACAutilities::getCoeffs(nutFieldsTest,
158 nutModes,
159 0, true);
160 coeffL2U_test.transposeInPlace();
161 coeffL2P_test.transposeInPlace();
162 coeffL2Nut_test.transposeInPlace();
163 cnpy::save(coeffL2U_test, "ITHACAoutput/NN/coeffs/coeffL2UTest.npy");
164 cnpy::save(coeffL2P_test, "ITHACAoutput/NN/coeffs/coeffL2PTest.npy");
165 cnpy::save(coeffL2Nut_test, "ITHACAoutput/NN/coeffs/coeffL2NutTest.npy");
166 }
167 }
168
169 // // This function computes the coefficients which are later used for training
170 // void getTurbNN()
171 // {
172 // if (!ITHACAutilities::check_folder("ITHACAoutput/NN/coeffs"))
173 // {
174 // mkDir("ITHACAoutput/NN/coeffs");
175 // // Read Fields for Train
176 // PtrList<volVectorField> UfieldTrain;
177 // PtrList<volScalarField> PfieldTrain;
178 // PtrList<volScalarField> nutFieldsTrain;
179 // ITHACAstream::readMiddleFields(UfieldTrain, _U(),
180 // "./ITHACAoutput/Offline/");
181 // ITHACAstream::readMiddleFields(PfieldTrain, _p(),
182 // "./ITHACAoutput/Offline/");
183 // auto nut_train = _mesh().lookupObject<volScalarField>("nut");
184 // ITHACAstream::readMiddleFields(nutFieldsTrain, nut_train,
185 // "./ITHACAoutput/Offline/");
186 // /// Compute the coefficients for train
187 // std::cout << "Computing the coefficients for U train" << std::endl;
188 // Eigen::MatrixXd coeffL2U_train = ITHACAutilities::getCoeffs(UfieldTrain,
189 // Umodes,
190 // 0, true);
191 // std::cout << "Computing the coefficients for p train" << std::endl;
192 // Eigen::MatrixXd coeffL2P_train = ITHACAutilities::getCoeffs(PfieldTrain,
193 // Pmodes,
194 // 0, true);
195 // std::cout << "Computing the coefficients for nuT train" << std::endl;
196 // Eigen::MatrixXd coeffL2Nut_train = ITHACAutilities::getCoeffs(nutFieldsTrain,
197 // nutModes,
198 // 0, true);
199 // coeffL2U_train.transposeInPlace();
200 // coeffL2P_train.transposeInPlace();
201 // coeffL2Nut_train.transposeInPlace();
202 // cnpy::save(coeffL2U_train, "ITHACAoutput/NN/coeffs/coeffL2UTrain.npy");
203 // cnpy::save(coeffL2P_train, "ITHACAoutput/NN/coeffs/coeffL2PTrain.npy");
204 // cnpy::save(coeffL2Nut_train, "ITHACAoutput/NN/coeffs/coeffL2NutTrain.npy");
205 // // Read Fields for Test
206 // PtrList<volVectorField> UfieldTest;
207 // PtrList<volScalarField> PfieldTest;
208 // PtrList<volScalarField> nutFieldsTest;
209 // /// Compute the coefficients for test
210 // ITHACAstream::readMiddleFields(UfieldTest, _U(),
211 // "./ITHACAoutput/checkOff/");
212 // ITHACAstream::readMiddleFields(PfieldTest, _p(),
213 // "./ITHACAoutput/checkOff/");
214 // auto nut_test = _mesh().lookupObject<volScalarField>("nut");
215 // ITHACAstream::readMiddleFields(nutFieldsTest, nut_test,
216 // "./ITHACAoutput/checkOff/");
217
218 // Info << "UfieldTest.size" << UfieldTest.size() << endl;
219 // //Info << "PfieldCheck.size" << PfieldCheck.size() << endl;
220 // //Info << "EfieldCheck.size" << EfieldCheck.size() << endl;
221 // //Info << "nutFieldsCheck.size" << nutFieldsCheck.size() << endl;
222 // //exit(0);
223
224 // // Compute the coefficients for test
225 // Eigen::MatrixXd coeffL2U_test = ITHACAutilities::getCoeffs(UfieldTest,
226 // Umodes,
227 // 0, true);
228 // Eigen::MatrixXd coeffL2P_test = ITHACAutilities::getCoeffs(PfieldTest,
229 // Pmodes,
230 // 0, true);
231 // Eigen::MatrixXd coeffL2Nut_test = ITHACAutilities::getCoeffs(nutFieldsTest,
232 // nutModes,
233 // 0, true);
234 // coeffL2U_test.transposeInPlace();
235 // coeffL2P_test.transposeInPlace();
236 // coeffL2Nut_test.transposeInPlace();
237 // cnpy::save(coeffL2U_test, "ITHACAoutput/NN/coeffs/coeffL2UTest.npy");
238 // cnpy::save(coeffL2P_test, "ITHACAoutput/NN/coeffs/coeffL2PTest.npy");
239 // cnpy::save(coeffL2Nut_test, "ITHACAoutput/NN/coeffs/coeffL2NutTest.npy");
240 // }
241 // }
242
243 Eigen::MatrixXd evalNet(Eigen::MatrixXd a, Eigen::MatrixXd mu_now)
244 {
245 Eigen::MatrixXd xpred(a.rows() + mu_now.rows(), 1);
246
247 if (xpred.cols() != 1)
248 {
249 throw std::runtime_error("Predictions for more than one sample not supported yet.");
250 }
251
252 xpred.bottomRows(a.rows()) = a;
253 xpred.topRows(mu_now.rows()) = mu_now;
254 xpred = xpred.array() * scale_inp.array() + bias_inp.array() ;
255 xpred.transposeInPlace();
256 torch::Tensor xTensor = eigenMatrix2torchTensor(xpred);
257 torch::Tensor out;
258 std::vector<torch::jit::IValue> XTensorInp;
259 XTensorInp.push_back(xTensor);
260 out = netTorchscript.forward(XTensorInp).toTensor();
261 Eigen::MatrixXd g = torchTensor2eigenMatrix<double>(out);
262 g.transposeInPlace();
263 g = (g.array() - bias_out.array()) / scale_out.array();
264 return g;
265 }
266};
267
268class ReducedCompressibleSteadyNN : public ReducedCompressibleSteadyNS
269{
270 public:
272 explicit ReducedCompressibleSteadyNN(CompressibleSteadyNN& FOMproblem)
273 :
274 problem(&FOMproblem),
275 ReducedCompressibleSteadyNS(FOMproblem)
276 {}
277
279 CompressibleSteadyNN* problem;
280
281 void projectReducedOperators(int NmodesUproj, int NmodesPproj, int NmodesEproj)
282 {
283 PtrList<volVectorField> gradModP;
284
285 for (label i = 0; i < NmodesPproj; i++)
286 {
287 gradModP.append(fvc::grad(problem->Pmodes[i]));
288 }
289
290 projGradModP = problem->Umodes.project(gradModP,
291 NmodesUproj); // Modes without lifting
292 }
293
294 void solveOnlineCompressible(int NmodesUproj, int NmodesPproj, int NmodesEproj,
295 int NmodesNutProj, Eigen::MatrixXd mu_now,
296 word Folder = "./ITHACAoutput/Online/")
297 {
298 counter++;
299 // Residuals initialization
300 scalar residualNorm(1);
301 scalar residualJump(1);
302 Eigen::MatrixXd uResidualOld = Eigen::MatrixXd::Zero(1, NmodesUproj);
303 Eigen::MatrixXd eResidualOld = Eigen::MatrixXd::Zero(1, NmodesEproj);
304 Eigen::MatrixXd pResidualOld = Eigen::MatrixXd::Zero(1, NmodesPproj);
305 Eigen::VectorXd uResidual(Eigen::Map<Eigen::VectorXd>(uResidualOld.data(),
306 NmodesUproj));
307 Eigen::VectorXd eResidual(Eigen::Map<Eigen::VectorXd>(eResidualOld.data(),
308 NmodesEproj));
309 Eigen::VectorXd pResidual(Eigen::Map<Eigen::VectorXd>(pResidualOld.data(),
310 NmodesPproj));
311 // Parameters definition
312 ITHACAparameters* para = ITHACAparameters::getInstance();
313 float residualJumpLim =
314 para->ITHACAdict->lookupOrDefault<float>("residualJumpLim", 1e-5);
315 float normalizedResidualLim =
316 para->ITHACAdict->lookupOrDefault<float>("normalizedResidualLim", 1e-5);
317 int maxIter =
318 para->ITHACAdict->lookupOrDefault<float>("maxIter", 2000);
319 bool closedVolume = false;
320 label csolve = 0;
321 // Full variables initialization
322 fluidThermo& thermo = problem->pThermo();
323 volVectorField& U = problem->_U();
324 //volScalarField& P = problem->pThermo().p();
325 volScalarField& P = problem->_p();
326 volScalarField& E = problem->pThermo->he();
328 volScalarField nueff = problem->turbulence->nuEff();
329 volScalarField& nut = const_cast<volScalarField&>
330 (problem->_mesh().lookupObject<volScalarField>("nut"));
331 volScalarField& rho = problem->_rho();
332 volScalarField& psi = problem->_psi();
333 surfaceScalarField& phi = problem->_phi();
334 Time& runTime = problem->_runTime();
335 fvMesh& mesh = problem->_mesh();
336 fv::options& fvOptions = problem->_fvOptions();
337 scalar cumulativeContErr = problem->cumulativeContErr;
338 // Reduced variables initialization
339 //Eigen::MatrixXd u = Eigen::MatrixXd::Zero(NmodesUproj, 1);
340 //Eigen::MatrixXd e = Eigen::MatrixXd::Zero(NmodesEproj, 1);
341 Eigen::MatrixXd e = ITHACAutilities::getCoeffs(E, problem->Emodes, NmodesEproj,
342 true);
343 Eigen::MatrixXd u = ITHACAutilities::getCoeffs(U, problem->Umodes, NmodesUproj,
344 true);
345 //Eigen::MatrixXd p = Eigen::MatrixXd::Zero(NmodesPproj, 1);
346 Eigen::MatrixXd p = ITHACAutilities::getCoeffs(P, problem->Pmodes, NmodesPproj,
347 true);
348 //Eigen::MatrixXd nutCoeff = ITHACAutilities::getCoeffs(nut, problem->nutModes, NmodesNutProj, true);
349 //problem->nutModes.reconstruct(nut, nutCoeff, "nut");
350 Eigen::MatrixXd nutCoeff = problem->evalNet(u, mu_now);
351 problem->nutModes.reconstruct(nut, nutCoeff, "nut");
352 //vector Uinlet(170,0,0); // Vector for the inlet boundary condition
353 label idInl =
354 problem->_mesh().boundaryMesh().findPatchID("inlet"); // ID of the inlet patch
355 vector Uinlet(problem->_U().boundaryFieldRef()[idInl][0][0], 0, 0);
356 P.storePrevIter();
357
358 while ((residualJump > residualJumpLim
359 || residualNorm > normalizedResidualLim) && csolve < maxIter)
360 {
361 csolve++;
362 Info << "csolve:" << csolve << endl;
363#if OFVER == 6
364 problem->_simple().loop(runTime);
365#else
366 problem->_simple().loop();
367#endif
368 uResidualOld = uResidual;
369 eResidualOld = eResidual;
370 pResidualOld = pResidual;
371 //Momentum equation phase
372 List<Eigen::MatrixXd> RedLinSysU;
373 ITHACAutilities::assignBC(U, idInl, Uinlet);
375 nueff = nut + problem->turbulence->nu();
376 fvVectorMatrix UEqnR
377 (
378 fvm::div(phi, U)
379 - fvc::div((rho * nueff) * dev2(T(fvc::grad(U))))
380 - fvm::laplacian(rho * nueff, U)
381 ==
382 fvOptions(rho, U)
383 );
384 UEqnR.relax();
385 fvOptions.constrain(UEqnR);
386 RedLinSysU = problem->Umodes.project(UEqnR,
387 NmodesUproj); // Modes without lifting
388 //volVectorField gradpfull = -fvc::grad(P);
389 //Eigen::MatrixXd projGrad = problem->Umodes.project(gradpfull, NmodesUproj);
390 Eigen::MatrixXd projGradP = projGradModP * p;
391 RedLinSysU[1] = RedLinSysU[1] - projGradP; // projGradP;
392 u = reducedProblem::solveLinearSys(RedLinSysU, u,
393 uResidual);//, "fullPivLu");//"bdcSvd");
394 problem->Umodes.reconstruct(U, u, "U");
395 ITHACAutilities::assignBC(U, idInl, Uinlet);
396 //solve(UEqnR == -problem->getGradP(P)); //For debug purposes only, second part only useful when using uEqn_global==-getGradP
397 fvOptions.correct(U);
398 //Energy equation phase
399 fvScalarMatrix EEqnR
400 (
401 fvm::div(phi, E)
402 + fvc::div(phi, volScalarField("Ekp", 0.5 * magSqr(U) + P / rho))
403 - fvm::laplacian(problem->turbulence->alphaEff(), E)
404 ==
405 fvOptions(rho, E)
406 );
407 EEqnR.relax();
408 fvOptions.constrain(EEqnR);
409 List<Eigen::MatrixXd> RedLinSysE = problem->Emodes.project(EEqnR, NmodesEproj);
410 e = reducedProblem::solveLinearSys(RedLinSysE, e, eResidual);
411 problem->Emodes.reconstruct(E, e, "e");
412 //EEqnR.solve(); //For debug purposes only
413 fvOptions.correct(E);
414 thermo.correct(); // Here are calculated both temperature and density based on P,U and he.
415 // Pressure equation phase
416 constrainPressure(P, rho, U, problem->getPhiHbyA(UEqnR, U, P),
417 problem->getRhorAUf(
418 UEqnR));// Update the pressure BCs to ensure flux consistency
419 surfaceScalarField phiHbyACalculated = problem->getPhiHbyA(UEqnR, U, P);
420 closedVolume = adjustPhi(phiHbyACalculated, U, P);
421 List<Eigen::MatrixXd> RedLinSysP;
422
423 while (problem->_simple().correctNonOrthogonal())
424 {
425 volScalarField rAU(1.0 /
426 UEqnR.A()); // Inverse of the diagonal part of the U equation matrix
427 volVectorField HbyA(constrainHbyA(rAU * UEqnR.H(), U,
428 P)); // H is the extra diagonal part summed to the r.h.s. of the U equation
429 surfaceScalarField phiHbyA("phiHbyA", fvc::interpolate(rho)*fvc::flux(HbyA));
430 surfaceScalarField rhorAUf("rhorAUf", fvc::interpolate(rho * rAU));
431 fvScalarMatrix PEqnR
432 (
433 fvc::div(phiHbyA)
434 - fvm::laplacian(rhorAUf, P)
435 ==
436 fvOptions(psi, P, rho.name())
437 );
438 PEqnR.setReference
439 (
440 problem->_pressureControl().refCell(),
441 problem->_pressureControl().refValue()
442 );
443 RedLinSysP = problem->Pmodes.project(PEqnR, NmodesPproj);
444 p = reducedProblem::solveLinearSys(RedLinSysP, p, pResidual);
445 problem->Pmodes.reconstruct(P, p, "p");
446
447 if (problem->_simple().finalNonOrthogonalIter())
448 {
449 phi = problem->getPhiHbyA(UEqnR, U, P) + PEqnR.flux();
450 }
451 }
452
453 //#include "continuityErrs.H"
454#include "incompressible/continuityErrs.H"
455 P.relax();// Explicitly relax pressure for momentum corrector
456 U = problem->HbyA() - (1.0 / UEqnR.A()) * problem->getGradP(P);
457 U.correctBoundaryConditions();
458 fvOptions.correct(U);
459 bool pLimited = problem->_pressureControl().limit(P);
460
461 // For closed-volume cases adjust the pressure and density levels to obey overall mass continuity
462 if (closedVolume)
463 {
464 P += (problem->_initialMass() - fvc::domainIntegrate(psi * P))
465 / fvc::domainIntegrate(psi);
466 }
467
468 if (pLimited || closedVolume)
469 {
470 P.correctBoundaryConditions();
471 }
472
473 rho = thermo.rho(); // Here rho is calculated as p*psi = p/(R*T)
474 rho.relax();
475 std::cout << "Ures = " << (uResidual.cwiseAbs()).sum() /
476 (RedLinSysU[1].cwiseAbs()).sum() << std::endl;
477 std::cout << "Eres = " << (eResidual.cwiseAbs()).sum() /
478 (RedLinSysE[1].cwiseAbs()).sum() << std::endl;
479 std::cout << "Pres = " << (pResidual.cwiseAbs()).sum() /
480 (RedLinSysP[1].cwiseAbs()).sum() << std::endl;
481 residualNorm = max(max((uResidual.cwiseAbs()).sum() /
482 (RedLinSysU[1].cwiseAbs()).sum(),
483 (pResidual.cwiseAbs()).sum() / (RedLinSysP[1].cwiseAbs()).sum()),
484 (eResidual.cwiseAbs()).sum() / (RedLinSysE[1].cwiseAbs()).sum());
485 residualJump = max(max(((uResidual - uResidualOld).cwiseAbs()).sum() /
486 (RedLinSysU[1].cwiseAbs()).sum(),
487 ((pResidual - pResidualOld).cwiseAbs()).sum() /
488 (RedLinSysP[1].cwiseAbs()).sum()),
489 ((eResidual - eResidualOld).cwiseAbs()).sum() /
490 (RedLinSysE[1].cwiseAbs()).sum());
491 //problem->turbulence->correct(); // Resolution of the full turbulence (debug purposes only)
492 nutCoeff = problem->evalNet(u, mu_now);
493 problem->nutModes.reconstruct(nut, nutCoeff, "nut");
494 }
495
496 //label k = 1;
497 // U.rename("Ur");
498 // P.rename("Pr");
499 // E.rename("Er");
500 // nut.rename("nutR");
501 ITHACAstream::exportSolution(U, name(counter), Folder);
502 ITHACAstream::exportSolution(P, name(counter), Folder);
503 ITHACAstream::exportSolution(E, name(counter), Folder);
504 ITHACAstream::exportSolution(nut, name(counter), Folder);
505 }
506
507};
508
509
510class tutorial03 : public CompressibleSteadyNN
511{
512 public:
514 explicit tutorial03(int argc, char* argv[])
515 :
516 CompressibleSteadyNN(argc, argv)
517 {
519 middleExport = para->ITHACAdict->lookupOrDefault<bool>("middleExport", true);
520 }
521
522 ITHACAparameters* para = ITHACAparameters::getInstance();
523
525 void offlineSolve(word folder = "./ITHACAoutput/Offline/")
526 {
527 //std::ofstream cpuTimes;
528 //double durationOff;
529 //cpuTimes.open(folder + "/cpuTimes", std::ios_base::app);
531 volVectorField& U = _U();
533 volScalarField& p = _p();
535 volScalarField& E = _E();
536
537 // if the offline solution is already performed read the fields
538 if (offline && !ITHACAutilities::check_folder("./ITHACAoutput/POD/1"))
539 {
540 ITHACAstream::readMiddleFields(Ufield, U, folder);
541 ITHACAstream::readMiddleFields(Efield, E, folder);
542 ITHACAstream::readMiddleFields(Pfield, p, folder);
544 auto nut = _mesh().lookupObject<volScalarField>("nut");
545 ITHACAstream::readMiddleFields(nutFields, nut, folder);
546 mu_samples = ITHACAstream::readMatrix("./parsOff_mat.txt");
547 }
548 // else perform offline stage
549 else if (!offline)
550 {
551 double UIFinit = para->ITHACAdict->lookupOrDefault<double>("UIFinit", 250);
552 Vector<double> Uinl(UIFinit, 0, 0);
553 //Vector<double> Uinl(250, 0, 0);
554
555 for (label i = 0; i < mu.rows(); i++)
556 {
557 //std::clock_t startOff;
558 //startOff = std::clock();
559 std::cout << "Current mu = " << mu(i, 0) << std::endl;
560 changeViscosity(mu(i, 0));
561 assignIF(_U(), Uinl);
562 truthSolve(folder);
563 //durationOff = (std::clock() - startOff);
564 //cpuTimes << durationOff << std::endl;
565 }
566 }
567
568 //cpuTimes.close();
569 }
570
571};
572
573int main(int argc, char* argv[])
574{
575 // Construct the tutorial object
576 tutorial03 example(argc, argv);
577 //Info << example.pThermo().p() << endl;
578 //Info << example._p() << endl;
579 // exit(0);
580 //std::clock_t startOff, startOn;
581 //double durationOn, durationOff;
582 std::cerr << "debug point 1" << std::endl;
583 ITHACAparameters* para = ITHACAparameters::getInstance();
584 //Eigen::MatrixXd parOff;
585 std::ifstream exFileOff("./parsOff_mat.txt");
586
587 if (exFileOff)
588 {
589 example.mu = ITHACAstream::readMatrix("./parsOff_mat.txt");
590 }
591 else
592 {
593 //example.mu = ITHACAutilities::rand(20, 1, 1.00e-05, 1.00e-2);
594 label OffNum = para->ITHACAdict->lookupOrDefault<label>("OffNum", 25);
595 example.mu = Eigen::VectorXd::LinSpaced(OffNum, 1.00e-05, 1.00e-02);
596 ITHACAstream::exportMatrix(example.mu, "parsOff", "eigen", "./");
597 }
598
599 Eigen::MatrixXd parsOn;
600 std::ifstream exFileOn("./parsOn_mat.txt");
601
602 if (exFileOn)
603 {
604 parsOn = ITHACAstream::readMatrix("./parsOn_mat.txt");
605 }
606 else
607 {
608 label OnNum = para->ITHACAdict->lookupOrDefault<label>("OnNum", 20);
609 parsOn = ITHACAutilities::rand(OnNum, 1, 1.00e-05, 1.00e-02);
610 ITHACAstream::exportMatrix(parsOn, "parsOn", "eigen", "./");
611 }
612
613 std::cerr << "debug point 2" << std::endl;
614 // Read some parameters from file
615 int NmodesUout = para->ITHACAdict->lookupOrDefault<int>("NmodesUout", 15);
616 int NmodesPout = para->ITHACAdict->lookupOrDefault<int>("NmodesPout", 15);
617 int NmodesEout = para->ITHACAdict->lookupOrDefault<int>("NmodesEout", 15);
618 int NmodesNutOut = para->ITHACAdict->lookupOrDefault<int>("NmodesNutOut", 15);
619 int NmodesUproj = para->ITHACAdict->lookupOrDefault<int>("NmodesUproj", 10);
620 int NmodesPproj = para->ITHACAdict->lookupOrDefault<int>("NmodesPproj", 10);
621 int NmodesEproj = para->ITHACAdict->lookupOrDefault<int>("NmodesEproj", 10);
622 int NmodesNutProj = para->ITHACAdict->lookupOrDefault<int>("NmodesNutProj", 10);
623 //Set the inlet boundaries patch 0 directions x and y
624 // example.inletIndex.resize(1, 2);
625 // example.inletIndex(0, 0) = 1;
626 // example.inletIndex(0, 1) = 0;
627 //Perform the offline solve
628 std::cerr << "debug point 3" << std::endl;
629 example.offlineSolve();
630 std::cerr << "debug point 4" << std::endl;
631 //Read the lift field
632 // ITHACAstream::read_fields(example.liftfield, example._U(), "./lift/");
633 // ITHACAutilities::normalizeFields(example.liftfield);
634 // Homogenize the snapshots
635 // example.computeLift(example.Ufield, example.liftfield, example.Uomfield);
636 // Perform POD on velocity and pressure and store the first 10 modes
637 // ITHACAPOD::getModes(example.Uomfield, example.Umodes, example._U().name(), example.podex, 0, 0,
638 // NmodesUout);
639 ITHACAPOD::getModes(example.Ufield, example.Umodes, example._U().name(),
640 example.podex, 0, 0,
641 NmodesUout);
642 ITHACAPOD::getModes(example.Pfield, example.Pmodes, example._p().name(),
643 example.podex, 0, 0,
644 NmodesPout);
645 ITHACAPOD::getModes(example.Efield, example.Emodes, example._E().name(),
646 example.podex, 0, 0,
647 NmodesEout);
648 ITHACAPOD::getModes(example.nutFields, example.nutModes, "nut", example.podex,
649 0, 0,
650 NmodesNutOut);
651 // Info << Foam::mag(example.Umodes[10] ) << endl;
652 // Info << Foam::max(Foam::mag(example.Umodes[10] )) << endl;
653 // exit(0);
654
655 if (!ITHACAutilities::check_folder("./ITHACAoutput/checkOff"))
656 {
657 tutorial03 checkOff(argc, argv);
658 checkOff.mu = ITHACAstream::readMatrix("./parsOn_mat.txt");
659 //Perform the offline solve
660 //checkOff.middleExport = para->ITHACAdict->lookupOrDefault<bool>("middleExport", true);
661 checkOff.offline = false;
662 checkOff.restart();
663 checkOff.offlineSolve("./ITHACAoutput/checkOff/");
664 //durationOff = (std::clock() - startOff);
665 checkOff.offline = true;
666 }
667 else //(ITHACAutilities::check_folder("./ITHACAoutput/checkOff"))
668 {
669 PtrList<volVectorField> UfieldCheck;
670 PtrList<volScalarField> PfieldCheck;
671 PtrList<volScalarField> EfieldCheck;
672 PtrList<volScalarField> nutFieldsCheck;
673 ITHACAstream::readMiddleFields(UfieldCheck, example._U(),
674 "./ITHACAoutput/checkOff/");
675 ITHACAstream::readMiddleFields(PfieldCheck, example._p(),
676 "./ITHACAoutput/checkOff/");
677 ITHACAstream::readMiddleFields(EfieldCheck, example._E(),
678 "./ITHACAoutput/checkOff/");
679 auto nutCheck = example._mesh().lookupObject<volScalarField>("nut");
680 ITHACAstream::readMiddleFields(nutFieldsCheck, nutCheck,
681 "./ITHACAoutput/checkOff/");
682 // Info << "UfieldCheck.size" << UfieldCheck.size() << endl;
683 // Info << "PfieldCheck.size" << PfieldCheck.size() << endl;
684 // Info << "EfieldCheck.size" << EfieldCheck.size() << endl;
685 // Info << "nutFieldsCheck.size" << nutFieldsCheck.size() << endl;
686 // exit(0);
688 Eigen::MatrixXd snapsCheck =
689 ITHACAstream::readMatrix("./ITHACAoutput/checkOff/snaps");
690 label fieldNum = 0;
691
692 for (label k = 0; k < snapsCheck.rows(); k++)
693 {
694 //Info << "snapsCheck(" << k << ",0)=" << snapsCheck(k,0)<< endl;
695 fieldNum = fieldNum + snapsCheck(k, 0);
696 //Info << "fieldNum" << fieldNum << endl;
697 ITHACAstream::exportSolution(UfieldCheck[fieldNum - 1], name(k + 1),
698 "./ITHACAoutput/checkOffSingle/");
699 ITHACAstream::exportSolution(PfieldCheck[fieldNum - 1], name(k + 1),
700 "./ITHACAoutput/checkOffSingle/");
701 ITHACAstream::exportSolution(EfieldCheck[fieldNum - 1], name(k + 1),
702 "./ITHACAoutput/checkOffSingle/");
703 ITHACAstream::exportSolution(nutFieldsCheck[fieldNum - 1], name(k + 1),
704 "./ITHACAoutput/checkOffSingle/");
705 //ITHACAutilities::createSymLink("./ITHACAoutput/checkOff/"+name(k+1)+"/polyMesh", "./ITHACAoutput/checkOffSingle/"+name(k+1)+"/");
706 }
707 }
708
709 // Create the coefficients to train the net
710 example.getTurbNN();
711 //Before loading the net, it has to be created through the python script
712 example.loadNet("ITHACAoutput/NN/Net_" + name(NmodesUproj) + "_" + name(
713 NmodesNutProj) + ".pt");
714 // Create the reduced object
715 ReducedCompressibleSteadyNN reduced(example);
716
717 if (!ITHACAutilities::check_folder("./ITHACAoutput/Online_" + name(
718 NmodesUproj) + "_" + name(NmodesNutProj) + "/" ) )
719 {
720 //Perform the online solutions
721 std::ofstream cpuTimes;
722 word OnlineFolder = "./ITHACAoutput/Online_" + name(NmodesUproj) + "_" + name(
723 NmodesNutProj) + "/";
724 cpuTimes.open(OnlineFolder + "/cpuTimes", std::ios_base::app);
725
726 for (label k = 0; k < parsOn.rows(); k++)
727 {
728 //scalar mu_now = parOn(k, 0);
729 std::clock_t startOn;
730 startOn = std::clock();
731 double durationOn;
732 Eigen::MatrixXd mu_now = parsOn.row(k);
733 mu_now.transposeInPlace();
734 example.changeViscosity(mu_now(0, 0));
735 // reduced.setOnlineVelocity(vel);
736 reduced.projectReducedOperators(NmodesUproj, NmodesPproj, NmodesEproj);
737 //std::cout << "############################" << std::endl;
738 example.restart();
739 example.turbulence->validate();
740 //std::cout << "##############################################################" << std::endl;
741 // reduced.solveOnlineCompressible(mu_now, NmodesUproj, NmodesPproj, NmodesEproj);
742 reduced.solveOnlineCompressible(NmodesUproj, NmodesPproj, NmodesEproj,
743 NmodesNutProj, mu_now, "./ITHACAoutput/Online_" + name(NmodesUproj) + "_" +
744 name(NmodesNutProj) + "/");
745 durationOn = std::clock() - startOn;
746 cpuTimes << durationOn << std::endl;
747 }
748 }
750 else if (ITHACAutilities::check_folder("./ITHACAoutput/Online_" + name(
751 NmodesUproj) + "_" + name(NmodesNutProj) + "/" ) )
752 {
753 PtrList<volVectorField> offlineU, onlineU;
754 PtrList<volScalarField> offlineP, onlineP;
755 PtrList<volScalarField> offlineE, onlineE;
756 PtrList<volScalarField> offlineNut, onlineNut;
758 ITHACAstream::read_fields(onlineU, example._U(),
759 "./ITHACAoutput/Online_" + name(NmodesUproj) + "_" + name(NmodesNutProj) + "/");
760 ITHACAstream::read_fields(onlineP, example._p(),
761 "./ITHACAoutput/Online_" + name(NmodesUproj) + "_" + name(NmodesNutProj) + "/");
762 ITHACAstream::read_fields(onlineE, example._E(),
763 "./ITHACAoutput/Online_" + name(NmodesUproj) + "_" + name(NmodesNutProj) + "/");
764 auto nut = example._mesh().lookupObject<volScalarField>("nut");
766 ITHACAstream::read_fields(onlineNut, nut,
767 "./ITHACAoutput/Online_" + name(NmodesUproj) + "_" + name(NmodesNutProj) + "/");
769 ITHACAstream::read_fields(offlineU, example._U(),
770 "./ITHACAoutput/checkOffSingle/");
771 ITHACAstream::read_fields(offlineP, example._p(),
772 "./ITHACAoutput/checkOffSingle/");
773 ITHACAstream::read_fields(offlineE, example._E(),
774 "./ITHACAoutput/checkOffSingle/");
775 ITHACAstream::read_fields(offlineNut, nut, "./ITHACAoutput/checkOffSingle/");
776 //PtrList<volVectorField> Uerrflds;
777 //PtrList<volScalarField> Perrflds,Eerrflds,NuTerrflds;
778
779 for (label j = 0; j < parsOn.rows(); j++)
780 {
781 volVectorField Ue = offlineU[j] - onlineU[j];
782 //auto err = foam2Eigen::foam2Eigen
783 auto u = ITHACAutilities::L2Norm(offlineU[j]);
784 Ue /= u;
785 //std::cout << "Ue = " << ITHACAutilities::L2Norm(Ue)/u << std::endl;
787 volScalarField Pe = offlineP[j] - onlineP[j];
788 auto p = ITHACAutilities::L2Norm(offlineP[j]);
789 Pe /= p;
790 //std::cout << "Pe = " << ITHACAutilities::L2Norm(Pe)/p << std::endl;
792 volScalarField Ee = offlineE[j] - onlineE[j];
793 auto e = ITHACAutilities::L2Norm(offlineE[j]);
794 Ee /= e;
795 //std::cout << "Ee = " << ITHACAutilities::L2Norm(Ee)/e << std::endl;
796 volScalarField Nute = offlineNut[j] - onlineNut[j];
797 auto n = ITHACAutilities::L2Norm(offlineNut[j]);
798 Nute /= n;
799 //std::cout << "Nute = " << ITHACAutilities::L2Norm(Nute)/n << std::endl;
800 Ue.rename("Ue");
801 Pe.rename("Pe");
802 Ee.rename("Ee");
803 Nute.rename("Nute");
804 ITHACAstream::exportSolution(Ue, name(j + 1),
805 "./ITHACAoutput/ErrorFields_" + name(NmodesUproj) + "_" + name(
806 NmodesNutProj) + "/");
807 ITHACAstream::exportSolution(Pe, name(j + 1),
808 "./ITHACAoutput/ErrorFields_" + name(NmodesUproj) + "_" + name(
809 NmodesNutProj) + "/");
810 ITHACAstream::exportSolution(Ee, name(j + 1),
811 "./ITHACAoutput/ErrorFields_" + name(NmodesUproj) + "_" + name(
812 NmodesNutProj) + "/");
813 ITHACAstream::exportSolution(Nute, name(j + 1),
814 "./ITHACAoutput/ErrorFields_" + name(NmodesUproj) + "_" + name(
815 NmodesNutProj) + "/");
816 }
817
818 Eigen::MatrixXd errorU = ITHACAutilities::errorL2Rel(offlineU, onlineU);
819 Eigen::MatrixXd errorP = ITHACAutilities::errorL2Rel(offlineP, onlineP);
820 Eigen::MatrixXd errorE = ITHACAutilities::errorL2Rel(offlineE, onlineE);
821 Eigen::MatrixXd errorNut = ITHACAutilities::errorL2Rel(offlineNut, onlineNut);
823 ITHACAstream::exportMatrix(errorU,
824 "errorU" + name(NmodesUproj) + "_" + name(NmodesNutProj), "python",
825 "./ITHACAoutput/ErrorFields_" + name(NmodesUproj) + "_" + name(
826 NmodesNutProj) + "/");
827 ITHACAstream::exportMatrix(errorP,
828 "errorP" + name(NmodesUproj) + "_" + name(NmodesNutProj), "python",
829 "./ITHACAoutput/ErrorFields_" + name(NmodesUproj) + "_" + name(
830 NmodesNutProj) + "/");
831 ITHACAstream::exportMatrix(errorE,
832 "errorE" + name(NmodesUproj) + "_" + name(NmodesNutProj), "python",
833 "./ITHACAoutput/ErrorFields_" + name(NmodesUproj) + "_" + name(
834 NmodesNutProj) + "/");
835 ITHACAstream::exportMatrix(errorNut,
836 "errorNut" + name(NmodesUproj) + "_" + name(NmodesNutProj), "python",
837 "./ITHACAoutput/ErrorFields_" + name(NmodesUproj) + "_" + name(
838 NmodesNutProj) + "/");
839 }
840
841 // if(ITHACAutilities::check_folder("./ITHACAoutput/checkOff"))
842 // {
843 // PtrList<volVectorField> UfieldCheck;
844 // PtrList<volScalarField> PfieldCheck;
845 // PtrList<volScalarField> EfieldCheck;
846 // PtrList<volScalarField> nutFieldsCheck;
847 // ITHACAstream::readMiddleFields(UfieldCheck, example._U(),"./ITHACAoutput/checkOff/");
848 // ITHACAstream::readMiddleFields(PfieldCheck, example._p(),"./ITHACAoutput/checkOff/");
849 // ITHACAstream::readMiddleFields(EfieldCheck, example._E(),"./ITHACAoutput/checkOff/");
850 // auto nutCheck = example._mesh().lookupObject<volScalarField>("nut");
851 // ITHACAstream::readMiddleFields(nutFieldsCheck, nutCheck, "./ITHACAoutput/checkOff/");
852 // // Info << "UfieldCheck.size" << UfieldCheck.size() << endl;
853 // // Info << "PfieldCheck.size" << PfieldCheck.size() << endl;
854 // // Info << "EfieldCheck.size" << EfieldCheck.size() << endl;
855 // // Info << "nutFieldsCheck.size" << nutFieldsCheck.size() << endl;
856 // // exit(0);
857 // //////////
858 // Eigen::MatrixXd snapsCheck = ITHACAstream::readMatrix("./ITHACAoutput/checkOff/snaps");
859 // label fieldNum = 0;
860 // for(label k=0; k<snapsCheck.rows(); k++)
861 // {
862 // Info << "snapsCheck(" << k << ",0)=" << snapsCheck(k,0)<< endl;
863 // fieldNum = fieldNum + snapsCheck(k,0);
864 // Info << "fieldNum" << fieldNum << endl;
865 // ITHACAstream::exportSolution(UfieldCheck[fieldNum-1], name(k+1), "./ITHACAoutput/checkOffSingle/");
866 // ITHACAstream::exportSolution(PfieldCheck[fieldNum-1], name(k+1), "./ITHACAoutput/checkOffSingle/");
867 // ITHACAstream::exportSolution(EfieldCheck[fieldNum-1], name(k+1), "./ITHACAoutput/checkOffSingle/");
868 // ITHACAstream::exportSolution(nutFieldsCheck[fieldNum-1], name(k+1), "./ITHACAoutput/checkOffSingle/");
869 // //ITHACAutilities::createSymLink("./ITHACAoutput/checkOff/"+name(k+1)+"/polyMesh", "./ITHACAoutput/checkOffSingle/"+name(k+1)+"/");
870 // }
871 // exit(0);
872 // ITHACAutilities::createSymLink("./0", "./ITHACAoutput/checkOffSingle/");
873 // ITHACAutilities::createSymLink("./system", "./ITHACAoutput/checkOffSingle/");
874 // ITHACAutilities::createSymLink("./constant", "./ITHACAoutput/checkOffSingle/");
875 // if(!ITHACAutilities::check_folder("./ITHACAoutput/Online_"+name(NmodesUproj)+"_"+name(NmodesNutProj)+"/" ) )
876 // {
877 // //Perform the online solutions
878 // std::ofstream cpuTimes;
879 // word OnlineFolder = "./ITHACAoutput/Online_" + name(NmodesUproj) + "_" + name(NmodesNutProj) + "/";
880 // cpuTimes.open(OnlineFolder + "/cpuTimes", std::ios_base::app);
881 // for (label k = 0; k < parsOn.rows(); k++)
882 // {
883 // //scalar mu_now = parOn(k, 0);
884 // std::clock_t startOn;
885 // startOn = std::clock();
886 // double durationOn;
887 // Eigen::MatrixXd mu_now = parsOn.row(k);
888 // mu_now.transposeInPlace();
889 // example.changeViscosity(mu_now(0,0));
890 // // reduced.setOnlineVelocity(vel);
891 // reduced.projectReducedOperators(NmodesUproj, NmodesPproj, NmodesEproj);
892 // //std::cout << "############################" << std::endl;
893 // example.restart();
894 // example.turbulence->validate();
895 // std::cout << "##############################################################" << std::endl;
896 // // reduced.solveOnlineCompressible(mu_now, NmodesUproj, NmodesPproj, NmodesEproj);
897 // reduced.solveOnlineCompressible(mu_now, NmodesUproj, NmodesPproj, NmodesEproj, NmodesNutProj, "./ITHACAoutput/Online_" + name(NmodesUproj) + "_"+name(NmodesNutProj) + "/");
898 // durationOn = std::clock() - startOn;
899 // cpuTimes << durationOn << std::endl;
900 // }
901 // }
902 // //// Read the files
903 // else if (ITHACAutilities::check_folder("./ITHACAoutput/Online_" + name(NmodesUproj) + "_" + name(NmodesNutProj) + "/" ) )
904 // {
905 // PtrList<volVectorField> offlineU, onlineU;
906 // PtrList<volScalarField> offlineP, onlineP;
907 // PtrList<volScalarField> offlineE, onlineE;
908 // PtrList<volScalarField>offlineNut, onlineNut;
909 // //////
910 // ITHACAstream::read_fields(onlineU,example._U(),"./ITHACAoutput/Online_" + name(NmodesUproj) + "_" + name(NmodesNutProj) + "/");
911 // ITHACAstream::read_fields(onlineP,example._p(),"./ITHACAoutput/Online_" + name(NmodesUproj) + "_" + name(NmodesNutProj) + "/");
912 // ITHACAstream::read_fields(onlineE,example._E(),"./ITHACAoutput/Online_" + name(NmodesUproj) + "_" + name(NmodesNutProj) + "/");
913 // auto nut = example._mesh().lookupObject<volScalarField>("nut");
914 // /// Save the online fields
915 // ITHACAstream::read_fields(onlineNut,nut,"./ITHACAoutput/Online_" + name(NmodesUproj) + "_" + name(NmodesNutProj) + "/");
916 // /// Save the errors
917 // ITHACAstream::read_fields(offlineU,example._U(),"./ITHACAoutput/checkOffSingle/");
918 // ITHACAstream::read_fields(offlineP,example._p(),"./ITHACAoutput/checkOffSingle/");
919 // ITHACAstream::read_fields(offlineE,example._E(),"./ITHACAoutput/checkOffSingle/");
920 // ITHACAstream::read_fields(offlineNut,nut, "./ITHACAoutput/checkOffSingle/");
921 // //PtrList<volVectorField> Uerrflds;
922 // //PtrList<volScalarField> Perrflds,Eerrflds,NuTerrflds;
923 // for(label j=0; j<parsOn.rows(); j++)
924 // {
925 // volVectorField Ue = offlineU[j] - onlineU[j];
926 // //auto err = foam2Eigen::foam2Eigen
927 // volScalarField Pe = offlineP[j] - onlineP[j];
928 // volScalarField Ee = offlineE[j] - onlineE[j];
929 // volScalarField Nute = offlineNut[j] - onlineNut[j];
930 // Ue.rename("Ue");
931 // Pe.rename("Pe");
932 // Ee.rename("Ee");
933 // Nute.rename("Nute");
934 // ITHACAstream::exportSolution(Ue, name(j+1), "./ITHACAoutput/ErrorFields_" + name(NmodesUproj) + "_" + name(NmodesNutProj) + "/");
935 // ITHACAstream::exportSolution(Pe, name(j+1), "./ITHACAoutput/ErrorFields_" + name(NmodesUproj) + "_" + name(NmodesNutProj) + "/");
936 // ITHACAstream::exportSolution(Ee, name(j+1), "./ITHACAoutput/ErrorFields_" + name(NmodesUproj) + "_" + name(NmodesNutProj) + "/");
937 // ITHACAstream::exportSolution(Nute, name(j+1), "./ITHACAoutput/ErrorFields_" + name(NmodesUproj) + "_" + name(NmodesNutProj) + "/");
938 // }
939 // Eigen::MatrixXd errorU = ITHACAutilities::errorL2Rel(offlineU, onlineU);
940 // Eigen::MatrixXd errorP = ITHACAutilities::errorL2Rel(offlineP, onlineP);
941 // Eigen::MatrixXd errorE = ITHACAutilities::errorL2Rel(offlineE, onlineE);
942 // Eigen::MatrixXd errorNut = ITHACAutilities::errorL2Rel(offlineNut, onlineNut);
943 // ///
944 // ITHACAstream::exportMatrix(errorU,"errorU" + name(NmodesUproj) + "_" + name(NmodesNutProj), "python", "./ITHACAoutput/ErrorFields_" + name(NmodesUproj) + "_" + name(NmodesNutProj) + "/");
945 // ITHACAstream::exportMatrix(errorP,"errorP" + name(NmodesUproj) + "_" + name(NmodesNutProj), "python", "./ITHACAoutput/ErrorFields_" + name(NmodesUproj) + "_" + name(NmodesNutProj) + "/");
946 // ITHACAstream::exportMatrix(errorE,"errorE" + name(NmodesUproj) + "_" + name(NmodesNutProj), "python", "./ITHACAoutput/ErrorFields_" + name(NmodesUproj) + "_" + name(NmodesNutProj) + "/");
947 // ITHACAstream::exportMatrix(errorNut,"errorNut" + name(NmodesUproj) + "_" + name(NmodesNutProj), "python", "./ITHACAoutput/ErrorFields_" + name(NmodesUproj) + "_" + name(NmodesNutProj) + "/");
948 // }
949 // }
950 // else
951 // {
952 // std::cerr << "CheckOff folder is missing, error analysis cannot be performed." << std::endl;
953 // }
954 exit(0);
955}
main
int main(int argc, char *argv[])
Definition 03compSteadyNS.C:573
CompressibleSteadyNS.H
Header file of the steadyNS class.
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...
fvOptions
fv::options & fvOptions
Definition NLsolve.H:25
csolve
label csolve
Definition NLsolve.H:42
maxIter
maxIter
Definition NLsolve.H:75
cpuTimes
std::ofstream cpuTimes
Definition NLsolve.H:58
closedVolume
bool closedVolume
Definition NLsolve.H:33
psi
volScalarField & psi
Definition NLsolve.H:29
para
ITHACAparameters * para
Definition NLsolve.H:40
nut
auto nut
Definition NLsolve.H:195
ReducedCompressibleSteadyNS.H
Header file of the reducedSteadyNS class.
CompressibleSteadyNN::CompressibleSteadyNN
CompressibleSteadyNN(int argc, char *argv[])
Constructor.
Definition 03compSteadyNS.C:47
CompressibleSteadyNN::scale_inp
Eigen::MatrixXd scale_inp
Definition 02compBump.C:72
CompressibleSteadyNN::bias_out
Eigen::MatrixXd bias_out
Definition 02compBump.C:73
CompressibleSteadyNN::coeffL2Nut
Eigen::MatrixXd coeffL2Nut
Definition 02compBump.C:76
CompressibleSteadyNN::coeffL2U
Eigen::MatrixXd coeffL2U
Definition 02compBump.C:77
CompressibleSteadyNN::bias_inp
Eigen::MatrixXd bias_inp
Definition 02compBump.C:71
CompressibleSteadyNN::optimizer
torch::optim::Optimizer * optimizer
Definition 02compBump.C:83
CompressibleSteadyNN::evalNet
Eigen::MatrixXd evalNet(Eigen::MatrixXd a, Eigen::MatrixXd mu_now)
Definition 03compSteadyNS.C:243
CompressibleSteadyNN::Net
torch::nn::Sequential Net
Definition 02compBump.C:82
CompressibleSteadyNN::coeffL2U_tensor
torch::Tensor coeffL2U_tensor
Definition 02compBump.C:79
CompressibleSteadyNN::coeffL2Nut_tensor
torch::Tensor coeffL2Nut_tensor
Definition 02compBump.C:80
CompressibleSteadyNN::scale_out
Eigen::MatrixXd scale_out
Definition 02compBump.C:74
CompressibleSteadyNN::netTorchscript
torch::jit::script::Module netTorchscript
Definition 02compBump.C:84
CompressibleSteadyNN::loadNet
void loadNet(word filename)
Definition 03compSteadyNS.C:85
CompressibleSteadyNS
Implementation of a parametrized full order steady NS problem and preparation of the the reduced ma...
Definition CompressibleSteadyNS.H:67
CompressibleSteadyNS::Efield
PtrList< volScalarField > Efield
List of pointers used to store the energy solutions.
Definition CompressibleSteadyNS.H:158
CompressibleSteadyNS::CompressibleSteadyNS
CompressibleSteadyNS()
Null constructor.
Definition CompressibleSteadyNS.C:42
CompressibleSteadyNS::_p
autoPtr< volScalarField > _p
Definition CompressibleSteadyNS.H:127
CompressibleSteadyNS::Emodes
volScalarModes Emodes
List of pointers used to form the energy modes.
Definition CompressibleSteadyNS.H:161
CompressibleSteadyNS::changeViscosity
void changeViscosity(double mu_new)
Function to change the viscosity.
Definition CompressibleSteadyNS.C:116
CompressibleSteadyNS::middleExport
bool middleExport
Export also intermediate fields.
Definition CompressibleSteadyNS.H:186
CompressibleSteadyNS::_E
autoPtr< volScalarField > _E
Definition CompressibleSteadyNS.H:137
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
ReducedCompressibleSteadyNN::projectReducedOperators
void projectReducedOperators(int NmodesUproj, int NmodesPproj, int NmodesEproj)
Definition 03compSteadyNS.C:281
ReducedCompressibleSteadyNN::problem
CompressibleSteadyNN * problem
Full problem.
Definition 02compBump.C:223
ReducedCompressibleSteadyNN::ReducedCompressibleSteadyNN
ReducedCompressibleSteadyNN(CompressibleSteadyNN &FOMproblem)
Constructor.
Definition 03compSteadyNS.C:272
ReducedCompressibleSteadyNN::solveOnlineCompressible
void solveOnlineCompressible(int NmodesUproj, int NmodesPproj, int NmodesEproj, int NmodesNutProj, Eigen::MatrixXd mu_now, word Folder="./ITHACAoutput/Online/")
Definition 03compSteadyNS.C:294
ReducedCompressibleSteadyNS
Class where it is implemented a reduced problem for the steady Navier-stokes problem.
Definition ReducedCompressibleSteadyNS.H:63
ReducedCompressibleSteadyNS::projGradModP
Eigen::MatrixXd projGradModP
Projected gradient of the pressure modes.
Definition ReducedCompressibleSteadyNS.H:132
SteadyNSTurb::nutFields
PtrList< volScalarField > nutFields
List of snapshots for the solution for eddy viscosity.
Definition SteadyNSTurb.H:77
SteadyNSTurb::nutModes
volScalarModes nutModes
List of POD modes for eddy viscosity.
Definition SteadyNSTurb.H:80
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
reducedSteadyNS::para
ITHACAparameters * para
parameters to be read from the ITHACAdict file
Definition ReducedSteadyNS.H:104
reductionProblem::mu_samples
Eigen::MatrixXd mu_samples
Matrix of parameters to be used for PODI, where each row corresponds to a sample point....
Definition reductionProblem.H:93
reductionProblem::assignIF
void assignIF(T &s, G &value)
Assign internal field condition.
Definition reductionProblemTemplates.C:34
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
reductionProblem::truthSolve
void truthSolve()
Perform a TruthSolve.
Definition reductionProblem.C:97
steadyNS::restart
void restart()
set U and P back to the values into the 0 folder
Definition steadyNS.C:2277
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::para
ITHACAparameters * para
Definition steadyNS.H:82
steadyNS::turbulence
autoPtr< incompressible::turbulenceModel > turbulence
Turbulence model.
Definition steadyNS.H:308
steadyNS::Uinl
autoPtr< volVectorField > Uinl
Initial dummy field with all Dirichlet boundary conditions.
Definition steadyNS.H:281
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
tutorial03::offlineSolve
void offlineSolve()
Perform an Offline solve.
Definition 03steadyNS.C:52
tutorial03::p
volScalarField & p
Pressure field.
Definition 03steadyNS.C:47
tutorial03::offlineSolve
void offlineSolve(word folder="./ITHACAoutput/Offline/")
Perform an Offline solve.
Definition 03compSteadyNS.C:525
tutorial03::para
ITHACAparameters * para
Definition 03compSteadyNS.C:522
tutorial03::U
volVectorField & U
Velocity field.
Definition 03steadyNS.C:45
tutorial03::tutorial03
tutorial03(int argc, char *argv[])
Constructor.
Definition 03compSteadyNS.C:514
T
volScalarField & T
Definition createT.H:46
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::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::L2Norm
double L2Norm(GeometricField< scalar, fvPatchField, volMesh > &field)
Definition ITHACAerror.C:41
ITHACAutilities::errorL2Rel
double errorL2Rel(GeometricField< T, fvPatchField, volMesh > &field1, GeometricField< T, fvPatchField, volMesh > &field2, List< label > *labels)
Computes the relative error between two geometric Fields in L2 norm.
Definition ITHACAerror.C:369
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::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
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)
rAU
volScalarField rAU(1.0/Ueqn.A())
constrainPressure
constrainPressure(p, rho, U, phiHbyA, rhorAtU, MRF)
phiHbyA
phiHbyA
Definition pcEqn.H:61
pLimited
bool pLimited
Definition pcEqn.H:98
HbyA
HbyA
Definition pcEqn.H:62
phi
surfaceScalarField & phi
Definition createFields.H:68
U
volVectorField & U
Definition createFields.H:42
cumulativeContErr
cumulativeContErr
Definition continuityErrs.H:37
p
volScalarField & p
Definition createFields.H:42
rho
volScalarField & rho
Definition createFields.H:67
thermo
fluidThermo & thermo
Definition createFields.H:33
rhorAUf
surfaceScalarField rhorAUf("rhorAUf", fvc::interpolate(rho *rAU))
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