Loading...
Searching...
No Matches
03compSteadyNS.C
1/*---------------------------------------------------------------------------*\
2 ██╗████████╗██╗ ██╗ █████╗ ██████╗ █████╗ ███████╗██╗ ██╗
3 ██║╚══██╔══╝██║ ██║██╔══██╗██╔════╝██╔══██╗ ██╔════╝██║ ██║
4 ██║ ██║ ███████║███████║██║ ███████║█████╗█████╗ ██║ ██║
5 ██║ ██║ ██╔══██║██╔══██║██║ ██╔══██║╚════╝██╔══╝ ╚██╗ ██╔╝
6 ██║ ██║ ██║ ██║██║ ██║╚██████╗██║ ██║ ██║ ╚████╔╝
7 ╚═╝ ╚═╝ ╚═╝ ╚═╝╚═╝ ╚═╝ ╚═════╝╚═╝ ╚═╝ ╚═╝ ╚═══╝
8
9 * In real Time Highly Advanced Computational Applications for Finite Volumes
10 * Copyright (C) 2017 by the ITHACA-FV authors
11-------------------------------------------------------------------------------
12License
13 This file is part of ITHACA-FV
14 ITHACA-FV is free software: you can redistribute it and/or modify
15 it under the terms of the GNU Lesser General Public License as published by
16 the Free Software Foundation, either version 3 of the License, or
17 (at your option) any later version.
18 ITHACA-FV is distributed in the hope that it will be useful,
19 but WITHOUT ANY WARRANTY; without even the implied warranty of
20 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
21 GNU Lesser General Public License for more details.
22 You should have received a copy of the GNU Lesser General Public License
23 along with ITHACA-FV. If not, see <http://www.gnu.org/licenses/>.
24Description
25 Example of steady NS Reduction Problem
26SourceFiles
27 03steadyNS.C
28\*---------------------------------------------------------------------------*/
29
30#include <torch/script.h>
31#include <torch/torch.h>
32#include "torch2Eigen.H"
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 Info << "Computing the coefficients for U train" << endl;
121 Eigen::MatrixXd coeffL2U_train = ITHACAutilities::getCoeffs(UfieldTrain,
122 Umodes,
123 0, true);
124 Info << "Computing the coefficients for p train" << endl;
125 Eigen::MatrixXd coeffL2P_train = ITHACAutilities::getCoeffs(PfieldTrain,
126 Pmodes,
127 0, true);
128 Info << "Computing the coefficients for nuT train" << 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 // Info << "Computing the coefficients for U train" << endl;
188 // Eigen::MatrixXd coeffL2U_train = ITHACAutilities::getCoeffs(UfieldTrain,
189 // Umodes,
190 // 0, true);
191 // Info << "Computing the coefficients for p train" << endl;
192 // Eigen::MatrixXd coeffL2P_train = ITHACAutilities::getCoeffs(PfieldTrain,
193 // Pmodes,
194 // 0, true);
195 // Info << "Computing the coefficients for nuT train" << 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 Info << "Ures = " << (uResidual.cwiseAbs()).sum() /
476 (RedLinSysU[1].cwiseAbs()).sum() << endl;
477 Info << "Eres = " << (eResidual.cwiseAbs()).sum() /
478 (RedLinSysE[1].cwiseAbs()).sum() << endl;
479 Info << "Pres = " << (pResidual.cwiseAbs()).sum() /
480 (RedLinSysP[1].cwiseAbs()).sum() << 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 Info << "Current mu = " << mu(i, 0) << endl;
560 changeViscosity(mu(i, 0));
561 assignIF(_U(), Uinl);
562 truthSolve(folder);
563 //durationOff = (std::clock() - startOff);
564 //cpuTimes << durationOff << 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" << 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" << 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" << endl;
629 example.offlineSolve();
630 std::cerr << "debug point 4" << 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 //Info << "############################" << endl;
738 example.restart();
739 example.turbulence->validate();
740 //Info << "##############################################################" << 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 << 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 //Info << "Ue = " << ITHACAutilities::L2Norm(Ue)/u << endl;
787 volScalarField Pe = offlineP[j] - onlineP[j];
788 auto p = ITHACAutilities::L2Norm(offlineP[j]);
789 Pe /= p;
790 //Info << "Pe = " << ITHACAutilities::L2Norm(Pe)/p << endl;
792 volScalarField Ee = offlineE[j] - onlineE[j];
793 auto e = ITHACAutilities::L2Norm(offlineE[j]);
794 Ee /= e;
795 //Info << "Ee = " << ITHACAutilities::L2Norm(Ee)/e << endl;
796 volScalarField Nute = offlineNut[j] - onlineNut[j];
797 auto n = ITHACAutilities::L2Norm(offlineNut[j]);
798 Nute /= n;
799 //Info << "Nute = " << ITHACAutilities::L2Norm(Nute)/n << 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 // //Info << "############################" << endl;
893 // example.restart();
894 // example.turbulence->validate();
895 // Info << "##############################################################" << 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 << 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." << endl;
953 // }
954 exit(0);
955}
CompressibleSteadyNS.H
Header file of the steadyNS class.
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...
ReducedCompressibleSteadyNS.H
Header file of the reducedSteadyNS class.
CompressibleSteadyNN::CompressibleSteadyNN
CompressibleSteadyNN(int argc, char *argv[])
Constructor.
Definition 03compSteadyNS.C:47
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::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
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::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:78
SteadyNSTurb::nutModes
volScalarModes nutModes
List of POD modes for eddy viscosity.
Definition SteadyNSTurb.H:81
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::truthSolve
void truthSolve()
Perform a TruthSolve.
Definition reductionProblem.C:97
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::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::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::U
volVectorField & U
Velocity field.
Definition 03steadyNS.C:45
tutorial03::tutorial03
tutorial03(int argc, char *argv[])
Constructor.
Definition 03compSteadyNS.C:514
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::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::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:337
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::L2Norm
double L2Norm(const T v)
Compute the L2 norm of v.
Definition ITHACAnorm.C:300
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
torch2Eigen.H
Header file of the torch2Eigen namespace.
Generated by   doxygen 1.16.1