Loading...
Searching...
No Matches
UnsteadyNSTurb.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-------------------------------------------------------------------------------
12
13 License
14 This file is part of ITHACA-FV
15
16 ITHACA-FV is free software: you can redistribute it and/or modify
17 it under the terms of the GNU Lesser General Public License as published by
18 the Free Software Foundation, either version 3 of the License, or
19 (at your option) any later version.
20
21 ITHACA-FV is distributed in the hope that it will be useful,
22 but WITHOUT ANY WARRANTY; without even the implied warranty of
23 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
24 GNU Lesser General Public License for more details.
25
26 You should have received a copy of the GNU Lesser General Public License
27 along with ITHACA-FV. If not, see <http://www.gnu.org/licenses/>.
28
29\*---------------------------------------------------------------------------*/
30
31#include "UnsteadyNSTurb.H"
32
35
36// * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * * * //
37
38// Construct Null
39UnsteadyNSTurb::UnsteadyNSTurb() {}
40
41// Construct from zero
42UnsteadyNSTurb::UnsteadyNSTurb(int argc, char* argv[])
43{
44 _args = autoPtr<argList>
45 (
46 new argList(argc, argv)
47 );
48
49 if (!_args->checkRootCase())
50 {
51 Foam::FatalError.exit();
52 }
53
54 argList& args = _args();
55#include "createTime.H"
56#include "createMesh.H"
57 _pimple = autoPtr<pimpleControl>
58 (
59 new pimpleControl
60 (
61 mesh
62 )
63 );
64#include "createFields.H"
65#include "createUfIfPresent.H"
66#include "createFvOptions.H"
67 ITHACAdict = new IOdictionary
68 (
69 IOobject
70 (
71 "ITHACAdict",
72 runTime.system(),
73 mesh,
74 IOobject::MUST_READ,
75 IOobject::NO_WRITE
76 )
77 );
78 tolerance = ITHACAdict->lookupOrDefault<scalar>("tolerance", 1e-5);
79 maxIter = ITHACAdict->lookupOrDefault<scalar>("maxIter", 1000);
80 bcMethod = ITHACAdict->lookupOrDefault<word>("bcMethod", "lift");
81 timeDerivativeSchemeOrder =
82 ITHACAdict->lookupOrDefault<word>("timeDerivativeSchemeOrder", "second");
83 M_Assert(bcMethod == "lift" || bcMethod == "penalty",
84 "The BC method must be set to lift or penalty in ITHACAdict");
85 M_Assert(timeDerivativeSchemeOrder == "first"
86 || timeDerivativeSchemeOrder == "second",
87 "The time derivative approximation must be set to either first or second order scheme in ITHACAdict");
88 para = ITHACAparameters::getInstance(mesh, runTime);
89 offline = ITHACAutilities::check_off();
90 podex = ITHACAutilities::check_pod();
91 supex = ITHACAutilities::check_sup();
93 NUmodesOut = para->ITHACAdict->lookupOrDefault<label>("NmodesUout", 15);
95 NPmodesOut = para->ITHACAdict->lookupOrDefault<label>("NmodesPout", 15);
97 NNutModesOut = para->ITHACAdict->lookupOrDefault<label>("NmodesNutOut", 15);
99 NUmodes = para->ITHACAdict->lookupOrDefault<label>("NmodesUproj", 10);
101 NSUPmodes = para->ITHACAdict->lookupOrDefault<label>("NmodesSUPproj", 10);
103 NPmodes = para->ITHACAdict->lookupOrDefault<label>("NmodesPproj", 10);
105 NNutModes = para->ITHACAdict->lookupOrDefault<label>("NmodesNutProj", 0);
106}
107
108// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
109
110void UnsteadyNSTurb::truthSolve(List<scalar> mu_now, std::string& offlinepath)
111{
112 Time& runTime = _runTime();
113 surfaceScalarField& phi = _phi();
114 fvMesh& mesh = _mesh();
115#include "initContinuityErrs.H"
116 fv::options& fvOptions = _fvOptions();
117 pimpleControl& pimple = _pimple();
118 volScalarField& p = _p();
119 volVectorField& U = _U();
120 volScalarField& nut = _nut();
121 IOMRFZoneList& MRF = _MRF();
122 singlePhaseTransportModel& laminarTransport = _laminarTransport();
123 instantList Times = runTime.times();
124 label& pRefCell = _pRefCell;
125 scalar& pRefValue = _pRefValue;
126
127 mesh.setFluxRequired(p.name());
128
129 runTime.setEndTime(finalTime);
130 runTime.setTime(Times[1], 1);
131 runTime.setDeltaT(timeStep);
132 nextWrite = startTime + writeEvery;
133
134 label nsnapshots = 0;
135
136 // Start the time loop
137 while (runTime.run())
138 {
139 #include "readTimeControls.H"
140 #include "CourantNo.H"
141 #include "setDeltaT.H"
142
143 ++runTime;
144
145 Info << "Time = " << runTime.timeName() << nl << endl;
146
147 // --- Pressure-velocity PIMPLE corrector loop
148 while (pimple.loop())
149 {
150 #include "UEqn.H"
151 // --- Pressure corrector loop
152 while (pimple.correct())
153 {
154 #include "pEqn.H"
155 }
156
157 if (pimple.turbCorr())
158 {
159 laminarTransport.correct();
160 turbulence->correct();
161 }
162
163 }
164
165 Info << "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
166 << " ClockTime = " << runTime.elapsedClockTime() << " s"
167 << nl << endl;
168
169 if (checkWrite(runTime))
170 {
171 nsnapshots += 1;
172 // Produces error when uncommented
173 // volScalarField nut = turbulence->nut().ref();
174 nut = turbulence->nut();
175
176 ITHACAstream::exportSolution(U, name(counter), offlinepath);
177 ITHACAstream::exportSolution(p, name(counter), offlinepath);
178 ITHACAstream::exportSolution(nut, name(counter), offlinepath);
179 std::ofstream of(offlinepath + name(counter) + "/" +
180 runTime.timeName());
181 Ufield.append(tmp<volVectorField>(U));
182 Pfield.append(tmp<volScalarField>(p));
183 nutFields.append(tmp<volScalarField>(nut));
184 counter++;
185 nextWrite += writeEvery;
186 writeMu(mu_now);
187 // --- Fill in the mu_samples with parameters (time, mu) to be used for the PODI sample points
188 mu_samples.conservativeResize(mu_samples.rows() + 1, mu_now.size() + 1);
189 mu_samples(mu_samples.rows() - 1, 0) = atof(runTime.timeName().c_str());
190
191 for (label i = 0; i < mu_now.size(); i++)
192 {
193 mu_samples(mu_samples.rows() - 1, i + 1) = mu_now[i];
194 }
195 }
196 }
197
198 // Resize to Unitary if not initialized by user (i.e. non-parametric problem)
199 if (mu.cols() == 0)
200 {
201 mu.resize(1, 1);
202 }
203
204 if (mu_samples.rows() == nsnapshots * mu.cols())
205 {
206 ITHACAstream::exportMatrix(mu_samples, "mu_samples", "eigen",
207 offlinepath);
208 }
209}
210
211Eigen::Tensor<double, 3> UnsteadyNSTurb::turbulenceTensor1(label NUmodes,
212 label NSUPmodes, label nNutModes)
213{
214 label cSize = NUmodes + NSUPmodes + liftfield.size();
215 Eigen::Tensor<double, 3> ct1Tensor;
216 ct1Tensor.resize(cSize, nNutModes, cSize);
217
218 for (label i = 0; i < cSize; i++)
219 {
220 for (label j = 0; j < nNutModes; j++)
221 {
222 for (label k = 0; k < cSize; k++)
223 {
224 ct1Tensor(i, j, k) = fvc::domainIntegrate(L_U_SUPmodes[i] & fvc::laplacian(
225 nutModes[j], L_U_SUPmodes[k])).value();
226 }
227 }
228 }
229
230 // Export the tensor
231 ITHACAstream::SaveDenseTensor(ct1Tensor, "./ITHACAoutput/Matrices/",
232 "ct1_" + name(liftfield.size()) + "_" + name(NUmodes) + "_" + name(
233 NSUPmodes) + "_" + name(nNutModes) + "_t");
234 return ct1Tensor;
235}
236
237
238Eigen::Tensor<double, 3> UnsteadyNSTurb::turbulenceAveTensor1(label NUmodes,
239 label NSUPmodes)
240{
241 label cSize = NUmodes + NSUPmodes + liftfield.size();
242 Eigen::Tensor<double, 3> ct1AveTensor;
243 label samplesNumber = nutAve.size();
244 ct1AveTensor.resize(cSize, samplesNumber, cSize);
245
246 for (label i = 0; i < cSize; i++)
247 {
248 for (label j = 0; j < samplesNumber; j++)
249 {
250 for (label k = 0; k < cSize; k++)
251 {
252 ct1AveTensor(i, j, k) = fvc::domainIntegrate(L_U_SUPmodes[i] & fvc::laplacian(
253 nutAve[j], L_U_SUPmodes[k])).value();
254 }
255 }
256 }
257
258 // Export the tensor
259 ITHACAstream::SaveDenseTensor(ct1AveTensor, "./ITHACAoutput/Matrices/",
260 "ct1Ave_" + name(liftfield.size()) + "_" + name(NUmodes) + "_" + name(
261 NSUPmodes) + "_t");
262 return ct1AveTensor;
263}
264
265Eigen::Tensor<double, 3> UnsteadyNSTurb::turbulencePPETensor1(label NUmodes,
266 label NSUPmodes, label NPmodes, label nNutModes)
267{
268 label cSize = NUmodes + NSUPmodes + liftfield.size();
269 Eigen::Tensor<double, 3> ct1PPETensor;
270 ct1PPETensor.resize(NPmodes, nNutModes, cSize);
271
272 for (label i = 0; i < NPmodes; i++)
273 {
274 for (label j = 0; j < nNutModes; j++)
275 {
276 for (label k = 0; k < cSize; k++)
277 {
278 // ct1PPETensor(i, j, k) = fvc::domainIntegrate(2 * Pmodes[i] * (fvc::laplacian(
279 // L_U_SUPmodes[k]) & fvc::grad(nutModes[j]))).value();
280 // ct1PPETensor(i, j, k) = fvc::domainIntegrate(Pmodes[i] * (fvc::div(
281 // fvc::laplacian(
282 // nutModes[j], L_U_SUPmodes[k])))).value();
283 ct1PPETensor(i, j, k) = fvc::domainIntegrate(fvc::grad(Pmodes[i]) & (
284 fvc::laplacian(
285 nutModes[j], L_U_SUPmodes[k]))).value();
286 }
287 }
288 }
289
290 // Export the tensor
291 ITHACAstream::SaveDenseTensor(ct1PPETensor, "./ITHACAoutput/Matrices/",
292 "ct1PPE_" + name(liftfield.size()) + "_" + name(NUmodes) + "_" + name(
293 NSUPmodes) + "_" + name(NPmodes) + "_" + name(nNutModes) + "_t");
294 return ct1PPETensor;
295}
296
297Eigen::Tensor<double, 3> UnsteadyNSTurb::turbulencePPEAveTensor1(label NUmodes,
298 label NSUPmodes, label NPmodes)
299{
300 label cSize = NUmodes + NSUPmodes + liftfield.size();
301 Eigen::Tensor<double, 3> ct1PPEAveTensor;
302 label samplesNumber = nutAve.size();
303 ct1PPEAveTensor.resize(NPmodes, samplesNumber, cSize);
304
305 for (label i = 0; i < NPmodes; i++)
306 {
307 for (label j = 0; j < samplesNumber; j++)
308 {
309 for (label k = 0; k < cSize; k++)
310 {
311 // ct1PPEAveTensor(i, j, k) = fvc::domainIntegrate(2 * Pmodes[i] * (fvc::laplacian(
312 // L_U_SUPmodes[k]) & fvc::grad(nutAve[j]))).value();
313 // ct1PPEAveTensor(i, j, k) = fvc::domainIntegrate(Pmodes[i] * (fvc::div(
314 // fvc::laplacian(
315 // nutAve[j], L_U_SUPmodes[k])))).value();
316 ct1PPEAveTensor(i, j, k) = fvc::domainIntegrate(fvc::grad(Pmodes[i]) & (
317 fvc::laplacian(
318 nutAve[j], L_U_SUPmodes[k]))).value();
319 }
320 }
321 }
322
323 // Export the tensor
324 ITHACAstream::SaveDenseTensor(ct1PPEAveTensor, "./ITHACAoutput/Matrices/",
325 "ct1PPEAve_" + name(liftfield.size()) + "_" + name(NUmodes) + "_" + name(
326 NSUPmodes) + "_" + name(NPmodes) + "_t");
327 return ct1PPEAveTensor;
328}
329
330Eigen::Tensor<double, 3> UnsteadyNSTurb::turbulenceTensor2(label NUmodes,
331 label NSUPmodes, label nNutModes)
332{
333 label cSize = NUmodes + NSUPmodes + liftfield.size();
334 Eigen::Tensor<double, 3> ct2Tensor;
335 ct2Tensor.resize(cSize, nNutModes, cSize);
336
337 for (label i = 0; i < cSize; i++)
338 {
339 for (label j = 0; j < nNutModes; j++)
340 {
341 for (label k = 0; k < cSize; k++)
342 {
343 ct2Tensor(i, j, k) = fvc::domainIntegrate(L_U_SUPmodes[i] & (fvc::div(
344 nutModes[j] * dev2((fvc::grad(L_U_SUPmodes[k]))().T())))).value();
345 }
346 }
347 }
348
349 // Export the tensor
350 ITHACAstream::SaveDenseTensor(ct2Tensor, "./ITHACAoutput/Matrices/",
351 "ct2_" + name(liftfield.size()) + "_" + name(NUmodes) + "_" + name(
352 NSUPmodes) + "_" + name(nNutModes) + "_t");
353 return ct2Tensor;
354}
355
356Eigen::Tensor<double, 3> UnsteadyNSTurb::turbulenceAveTensor2(label NUmodes,
357 label NSUPmodes)
358{
359 label cSize = NUmodes + NSUPmodes + liftfield.size();
360 Eigen::Tensor<double, 3> ct2AveTensor;
361 label samplesNumber = nutAve.size();
362 ct2AveTensor.resize(cSize, samplesNumber, cSize);
363
364 for (label i = 0; i < cSize; i++)
365 {
366 for (label j = 0; j < samplesNumber; j++)
367 {
368 for (label k = 0; k < cSize; k++)
369 {
370 ct2AveTensor(i, j, k) = fvc::domainIntegrate(L_U_SUPmodes[i] & (fvc::div(
371 nutAve[j] * dev2((fvc::grad(L_U_SUPmodes[k]))().T())))).value();
372 }
373 }
374 }
375
376 // Export the tensor
377 ITHACAstream::SaveDenseTensor(ct2AveTensor, "./ITHACAoutput/Matrices/",
378 "ct2Ave_" + name(liftfield.size()) + "_" + name(NUmodes) + "_" + name(
379 NSUPmodes) + "_t");
380 return ct2AveTensor;
381}
382
383Eigen::Tensor<double, 3> UnsteadyNSTurb::turbulencePPETensor2(label NUmodes,
384 label NSUPmodes, label NPmodes, label nNutModes)
385{
386 label cSize = NUmodes + NSUPmodes + liftfield.size();
387 Eigen::Tensor<double, 3> ct2PPETensor;
388 ct2PPETensor.resize(NPmodes, nNutModes, cSize);
389
390 for (label i = 0; i < NPmodes; i++)
391 {
392 for (label j = 0; j < nNutModes; j++)
393 {
394 for (label k = 0; k < cSize; k++)
395 {
396 // ct2PPETensor(i, j, k) = fvc::domainIntegrate(Pmodes[i] * (fvc::grad(fvc::grad(
397 // nutModes[j])) && (dev2((fvc::grad(L_U_SUPmodes[k]))() + fvc::grad(
398 // L_U_SUPmodes[k]))().T()))).value();
399 // ct2PPETensor(i, j, k) = fvc::domainIntegrate(Pmodes[i] * ((fvc::div(fvc::div(
400 // nutModes[j] * dev2((fvc::grad(L_U_SUPmodes[k]))().T())))))).value();
401 ct2PPETensor(i, j, k) = fvc::domainIntegrate(fvc::grad(Pmodes[i]) & ((fvc::div(
402 nutModes[j] * dev2((fvc::grad(L_U_SUPmodes[k]))().T()))))).value();
403 }
404 }
405 }
406
407 // Export the tensor
408 ITHACAstream::SaveDenseTensor(ct2PPETensor, "./ITHACAoutput/Matrices/",
409 "ct2PPE_" + name(liftfield.size()) + "_" + name(NUmodes) + "_" + name(
410 NSUPmodes) + "_" + name(NPmodes) + "_" + name(nNutModes) + "_t");
411 return ct2PPETensor;
412}
413
414Eigen::Tensor<double, 3> UnsteadyNSTurb::turbulencePPEAveTensor2(label NUmodes,
415 label NSUPmodes, label NPmodes)
416{
417 label cSize = NUmodes + NSUPmodes + liftfield.size();
418 Eigen::Tensor<double, 3> ct2PPEAveTensor;
419 label samplesNumber = nutAve.size();
420 ct2PPEAveTensor.resize(NPmodes, samplesNumber, cSize);
421
422 for (label i = 0; i < NPmodes; i++)
423 {
424 for (label j = 0; j < samplesNumber; j++)
425 {
426 for (label k = 0; k < cSize; k++)
427 {
428 // ct2PPEAveTensor(i, j, k) = fvc::domainIntegrate(Pmodes[i] * (fvc::grad(
429 // fvc::grad(
430 // nutAve[j])) && (dev2((fvc::grad(L_U_SUPmodes[k]))() + fvc::grad(
431 // L_U_SUPmodes[k]))().T()))).value();
432 // ct2PPEAveTensor(i, j, k) = fvc::domainIntegrate(Pmodes[i] * ((fvc::div(fvc::div(
433 // nutAve[j] * dev2((fvc::grad(L_U_SUPmodes[k]))().T())))))).value();
434 ct2PPEAveTensor(i, j, k) = fvc::domainIntegrate(fvc::grad(Pmodes[i]) & ((
435 fvc::div(
436 nutAve[j] * dev2((fvc::grad(L_U_SUPmodes[k]))().T()))))).value();
437 }
438 }
439 }
440
441 // Export the tensor
442 ITHACAstream::SaveDenseTensor(ct2PPEAveTensor, "./ITHACAoutput/Matrices/",
443 "ct2PPEAve_" + name(liftfield.size()) + "_" + name(NUmodes) + "_" + name(
444 NSUPmodes) + "_" + name(NPmodes) + "_t");
445 return ct2PPEAveTensor;
446}
447
448Eigen::MatrixXd UnsteadyNSTurb::btTurbulence(label NUmodes, label NSUPmodes)
449{
450 label btSize = NUmodes + NSUPmodes + liftfield.size();
451 Eigen::MatrixXd btMatrix(btSize, btSize);
452 btMatrix = btMatrix * 0;
453
454 // Project everything
455 for (label i = 0; i < btSize; i++)
456 {
457 for (label j = 0; j < btSize; j++)
458 {
459 btMatrix(i, j) = fvc::domainIntegrate(L_U_SUPmodes[i] & (fvc::div(dev2((T(
460 fvc::grad(
461 L_U_SUPmodes[j]))))))).value();
462 }
463 }
464
465 // Export the matrix
466 ITHACAstream::SaveDenseMatrix(btMatrix, "./ITHACAoutput/Matrices/",
467 "bt_" + name(liftfield.size()) + "_" + name(NUmodes) + "_" + name(NSUPmodes));
468 return btMatrix;
469}
470
471void UnsteadyNSTurb::projectSUP(fileName folder, label NU, label NP, label NSUP,
472 label Nnut, bool rbfInterp)
473{
474 NUmodes = NU;
475 NPmodes = NP;
476 NSUPmodes = NSUP;
477 nNutModes = Nnut;
478 L_U_SUPmodes.resize(0);
479
480 if (liftfield.size() != 0)
481 {
482 for (label k = 0; k < liftfield.size(); k++)
483 {
484 L_U_SUPmodes.append(tmp<volVectorField>(liftfield[k]));
485 }
486 }
487
488 if (NUmodes != 0)
489 {
490 for (label k = 0; k < NUmodes; k++)
491 {
492 L_U_SUPmodes.append(tmp<volVectorField>(Umodes[k]));
493 }
494 }
495
496 if (NSUPmodes != 0)
497 {
498 for (label k = 0; k < NSUPmodes; k++)
499 {
500 L_U_SUPmodes.append(tmp<volVectorField>(supmodes[k]));
501 }
502 }
503
504 if (ITHACAutilities::check_folder("./ITHACAoutput/Matrices/"))
505 {
506 word bStr = "B_" + name(liftfield.size()) + "_" + name(NUmodes) + "_" + name(
507 NSUPmodes);
508
509 if (ITHACAutilities::check_file("./ITHACAoutput/Matrices/" + bStr))
510 {
511 ITHACAstream::ReadDenseMatrix(B_matrix, "./ITHACAoutput/Matrices/", bStr);
512 }
513 else
514 {
515 B_matrix = diffusive_term(NUmodes, NPmodes, NSUPmodes);
516 }
517
518 word btStr = "bt_" + name(liftfield.size()) + "_" + name(NUmodes) + "_" + name(
519 NSUPmodes);
520
521 if (ITHACAutilities::check_file("./ITHACAoutput/Matrices/" + btStr))
522 {
523 ITHACAstream::ReadDenseMatrix(btMatrix, "./ITHACAoutput/Matrices/", btStr);
524 }
525 else
526 {
527 btMatrix = btTurbulence(NUmodes, NSUPmodes);
528 }
529
530 word kStr = "K_" + name(liftfield.size()) + "_" + name(NUmodes) + "_" + name(
531 NSUPmodes) + "_" + name(NPmodes);
532
533 if (ITHACAutilities::check_file("./ITHACAoutput/Matrices/" + kStr))
534 {
535 ITHACAstream::ReadDenseMatrix(K_matrix, "./ITHACAoutput/Matrices/", kStr);
536 }
537 else
538 {
539 K_matrix = pressure_gradient_term(NUmodes, NPmodes, NSUPmodes);
540 }
541
542 word pStr = "P_" + name(liftfield.size()) + "_" + name(NUmodes) + "_" + name(
543 NSUPmodes) + "_" + name(NPmodes);
544
545 if (ITHACAutilities::check_file("./ITHACAoutput/Matrices/" + pStr))
546 {
547 ITHACAstream::ReadDenseMatrix(P_matrix, "./ITHACAoutput/Matrices/", pStr);
548 }
549 else
550 {
551 P_matrix = divergence_term(NUmodes, NPmodes, NSUPmodes);
552 }
553
554 word mStr = "M_" + name(liftfield.size()) + "_" + name(NUmodes) + "_" + name(
555 NSUPmodes);
556
557 if (ITHACAutilities::check_file("./ITHACAoutput/Matrices/" + mStr))
558 {
559 ITHACAstream::ReadDenseMatrix(M_matrix, "./ITHACAoutput/Matrices/", mStr);
560 }
561 else
562 {
563 M_matrix = mass_term(NUmodes, NPmodes, NSUPmodes);
564 }
565
566 word C_str = "C_" + name(liftfield.size()) + "_" + name(NUmodes) + "_" + name(
567 NSUPmodes) + "_t";
568
569 if (ITHACAutilities::check_file("./ITHACAoutput/Matrices/" + C_str))
570 {
571 ITHACAstream::ReadDenseTensor(C_tensor, "./ITHACAoutput/Matrices/", C_str);
572 }
573 else
574 {
575 C_tensor = convective_term_tens(NUmodes, NPmodes, NSUPmodes);
576 }
577
578 word ct1Str = "ct1_" + name(liftfield.size()) + "_" + name(
579 NUmodes) + "_" + name(
580 NSUPmodes) + "_" + name(nNutModes) + "_t";
581
582 if (ITHACAutilities::check_file("./ITHACAoutput/Matrices/" + ct1Str))
583 {
584 ITHACAstream::ReadDenseTensor(ct1Tensor, "./ITHACAoutput/Matrices/", ct1Str);
585 }
586 else
587 {
588 ct1Tensor = turbulenceTensor1(NUmodes, NSUPmodes, nNutModes);
589 }
590
591 word ct2Str = "ct2_" + name(liftfield.size()) + "_" + name(
592 NUmodes) + "_" + name(
593 NSUPmodes) + "_" + name(nNutModes) + "_t";
594
595 if (ITHACAutilities::check_file("./ITHACAoutput/Matrices/" + ct2Str))
596 {
597 ITHACAstream::ReadDenseTensor(ct2Tensor, "./ITHACAoutput/Matrices/", ct2Str);
598 }
599 else
600 {
601 ct2Tensor = turbulenceTensor2(NUmodes, NSUPmodes, nNutModes);
602 }
603
604 if (nutAve.size() != 0)
605 {
606 word ct1AveStr = "ct1Ave_" + name(liftfield.size()) + "_" + name(
607 NUmodes) + "_" + name(
608 NSUPmodes) + "_t";
609
610 if (ITHACAutilities::check_file("./ITHACAoutput/Matrices/" + ct1AveStr))
611 {
612 ITHACAstream::ReadDenseTensor(ct1AveTensor, "./ITHACAoutput/Matrices/",
613 ct1AveStr);
614 }
615 else
616 {
617 ct1AveTensor = turbulenceAveTensor1(NUmodes, NSUPmodes);
618 }
619
620 word ct2AveStr = "ct2Ave_" + name(liftfield.size()) + "_" + name(
621 NUmodes) + "_" + name(
622 NSUPmodes) + "_t";
623
624 if (ITHACAutilities::check_file("./ITHACAoutput/Matrices/" + ct2AveStr))
625 {
626 ITHACAstream::ReadDenseTensor(ct2AveTensor, "./ITHACAoutput/Matrices/",
627 ct2AveStr);
628 }
629 else
630 {
631 ct2AveTensor = turbulenceAveTensor2(NUmodes, NSUPmodes);
632 }
633 }
634
635 if (bcMethod == "penalty")
636 {
637 bcVelVec = bcVelocityVec(NUmodes, NSUPmodes);
638 bcVelMat = bcVelocityMat(NUmodes, NSUPmodes);
639 }
640 }
641 else
642 {
643 B_matrix = diffusive_term(NUmodes, NPmodes, NSUPmodes);
644 C_tensor = convective_term_tens(NUmodes, NPmodes, NSUPmodes);
645 K_matrix = pressure_gradient_term(NUmodes, NPmodes, NSUPmodes);
646 P_matrix = divergence_term(NUmodes, NPmodes, NSUPmodes);
647 M_matrix = mass_term(NUmodes, NPmodes, NSUPmodes);
648 btMatrix = btTurbulence(NUmodes, NSUPmodes);
649 ct1Tensor = turbulenceTensor1(NUmodes, NSUPmodes, nNutModes);
650 ct2Tensor = turbulenceTensor2(NUmodes, NSUPmodes, nNutModes);
651
652 if (nutAve.size() != 0)
653 {
654 ct1AveTensor = turbulenceAveTensor1(NUmodes, NSUPmodes);
655 ct2AveTensor = turbulenceAveTensor2(NUmodes, NSUPmodes);
656 }
657
658 if (bcMethod == "penalty")
659 {
660 bcVelVec = bcVelocityVec(NUmodes, NSUPmodes);
661 bcVelMat = bcVelocityMat(NUmodes, NSUPmodes);
662 }
663 }
664
665 // Export the matrices
666 if (para->exportPython)
667 {
668 ITHACAstream::exportMatrix(B_matrix, "B", "python", "./ITHACAoutput/Matrices/");
669 ITHACAstream::exportMatrix(K_matrix, "K", "python", "./ITHACAoutput/Matrices/");
670 ITHACAstream::exportMatrix(P_matrix, "P", "python", "./ITHACAoutput/Matrices/");
671 ITHACAstream::exportMatrix(M_matrix, "M", "python", "./ITHACAoutput/Matrices/");
672 ITHACAstream::exportMatrix(btMatrix, "bt", "python",
673 "./ITHACAoutput/Matrices/");
674 ITHACAstream::exportTensor(C_tensor, "C", "python", "./ITHACAoutput/Matrices/");
675 ITHACAstream::exportTensor(ct1Tensor, "ct1", "python",
676 "./ITHACAoutput/Matrices/");
677 ITHACAstream::exportTensor(ct2Tensor, "ct2", "python",
678 "./ITHACAoutput/Matrices/");
679
680 if (nutAve.size() != 0)
681 {
682 ITHACAstream::exportTensor(ct1AveTensor, "ct1Ave", "python",
683 "./ITHACAoutput/Matrices/");
684 ITHACAstream::exportTensor(ct2AveTensor, "ct2Ave", "python",
685 "./ITHACAoutput/Matrices/");
686 }
687 }
688
689 if (para->exportMatlab)
690 {
691 ITHACAstream::exportMatrix(B_matrix, "B", "matlab", "./ITHACAoutput/Matrices/");
692 ITHACAstream::exportMatrix(K_matrix, "K", "matlab", "./ITHACAoutput/Matrices/");
693 ITHACAstream::exportMatrix(P_matrix, "P", "matlab", "./ITHACAoutput/Matrices/");
694 ITHACAstream::exportMatrix(M_matrix, "M", "matlab", "./ITHACAoutput/Matrices/");
695 ITHACAstream::exportMatrix(btMatrix, "bt", "matlab",
696 "./ITHACAoutput/Matrices/");
697 ITHACAstream::exportTensor(C_tensor, "C", "matlab", "./ITHACAoutput/Matrices/");
698 ITHACAstream::exportTensor(ct1Tensor, "ct1", "matlab",
699 "./ITHACAoutput/Matrices/");
700 ITHACAstream::exportTensor(ct2Tensor, "ct2", "matlab",
701 "./ITHACAoutput/Matrices/");
702
703 if (nutAve.size() != 0)
704 {
705 ITHACAstream::exportTensor(ct1AveTensor, "ct1Ave", "matlab",
706 "./ITHACAoutput/Matrices/");
707 ITHACAstream::exportTensor(ct2AveTensor, "ct2Ave", "matlab",
708 "./ITHACAoutput/Matrices/");
709 }
710 }
711
712 if (para->exportTxt)
713 {
714 ITHACAstream::exportMatrix(B_matrix, "B", "eigen", "./ITHACAoutput/Matrices/");
715 ITHACAstream::exportMatrix(K_matrix, "K", "eigen", "./ITHACAoutput/Matrices/");
716 ITHACAstream::exportMatrix(P_matrix, "P", "eigen", "./ITHACAoutput/Matrices/");
717 ITHACAstream::exportMatrix(M_matrix, "M", "eigen", "./ITHACAoutput/Matrices/");
718 ITHACAstream::exportMatrix(btMatrix, "bt", "eigen", "./ITHACAoutput/Matrices/");
719 ITHACAstream::exportTensor(C_tensor, "C", "eigen", "./ITHACAoutput/Matrices/C");
720 ITHACAstream::exportTensor(ct1Tensor, "ct1_", "eigen",
721 "./ITHACAoutput/Matrices/ct1");
722 ITHACAstream::exportTensor(ct2Tensor, "ct2_", "eigen",
723 "./ITHACAoutput/Matrices/ct2");
724
725 if (nutAve.size() != 0)
726 {
727 ITHACAstream::exportTensor(ct1AveTensor, "ct1Ave_", "eigen",
728 "./ITHACAoutput/Matrices/ct1Ave");
729 ITHACAstream::exportTensor(ct2AveTensor, "ct2Ave_", "eigen",
730 "./ITHACAoutput/Matrices/ct2Ave");
731 }
732 }
733
734 bTotalMatrix = B_matrix + btMatrix;
735 label cSize = NUmodes + NSUPmodes + liftfield.size();
736 cTotalTensor.resize(cSize, nNutModes, cSize);
737 cTotalTensor = ct1Tensor + ct2Tensor;
738 // Define coeffL2
739 coeffL2 = ITHACAutilities::getCoeffs(nutFields,
740 nutModes, nNutModes);
741 ITHACAstream::exportMatrix(coeffL2, "coeffL2", "python",
742 "./ITHACAoutput/Matrices/");
743 ITHACAstream::exportMatrix(coeffL2, "coeffL2", "matlab",
744 "./ITHACAoutput/Matrices/");
745
746 if (nutAve.size() != 0)
747 {
748 cTotalAveTensor.resize(cSize, nutAve.size(), cSize);
749 cTotalAveTensor = ct1AveTensor + ct2AveTensor;
750 }
751
752 if (rbfInterp == true && (!Pstream::parRun()))
753 {
754 if (ITHACAutilities::check_file("./radii.txt"))
755 {
756 radii = ITHACAstream::readMatrix("./radii.txt");
757 M_Assert(radii.size() == nNutModes,
758 "Thes size of the shape parameters vector must be equal to the number of eddy viscosity modes nNutModes");
759 }
760 else
761 {
762 radii = Eigen::MatrixXd::Ones(nNutModes,
763 1) * e;
764 }
765
766 samples.resize(nNutModes);
767 rbfSplines.resize(nNutModes);
768 Eigen::MatrixXd weights;
769
770 for (label i = 0; i < nNutModes; i++)
771 {
772 word weightName = "wRBF_N" + name(i + 1) + "_" + name(liftfield.size()) + "_"
773 + name(NUmodes) + "_" + name(NSUPmodes) ;
774
775 if (ITHACAutilities::check_file("./ITHACAoutput/weightsSUP/" + weightName))
776 {
777 samples[i] = new SPLINTER::DataTable(1, 1);
778
779 for (label j = 0; j < coeffL2.cols(); j++)
780 {
781 samples[i]->addSample(velRBF.row(j), coeffL2(i, j));
782 }
783
784 ITHACAstream::ReadDenseMatrix(weights, "./ITHACAoutput/weightsSUP/",
785 weightName);
786 rbfSplines[i] = new SPLINTER::RBFSpline(*samples[i],
787 SPLINTER::RadialBasisFunctionType::GAUSSIAN, weights, radii(i));
788 std::cout << "Constructing RadialBasisFunction for mode " << i + 1 << std::endl;
789 }
790 else
791 {
792 samples[i] = new SPLINTER::DataTable(1, 1);
793
794 for (label j = 0; j < coeffL2.cols(); j++)
795 {
796 samples[i]->addSample(velRBF.row(j), coeffL2(i, j));
797 }
798
799 rbfSplines[i] = new SPLINTER::RBFSpline(*samples[i],
800 SPLINTER::RadialBasisFunctionType::GAUSSIAN, false, radii(i));
801 ITHACAstream::SaveDenseMatrix(rbfSplines[i]->weights,
802 "./ITHACAoutput/weightsSUP/", weightName);
803 std::cout << "Constructing RadialBasisFunction for mode " << i + 1 << std::endl;
804 }
805 }
806 }
807}
808
809
810
811void UnsteadyNSTurb::projectPPE(fileName folder, label NU, label NP, label NSUP,
812 label Nnut, bool rbfInterp)
813{
814 NUmodes = NU;
815 NPmodes = NP;
816 NSUPmodes = 0;
817 nNutModes = Nnut;
818 L_U_SUPmodes.resize(0);
819
820 if (liftfield.size() != 0)
821 {
822 for (label k = 0; k < liftfield.size(); k++)
823 {
824 L_U_SUPmodes.append(tmp<volVectorField>(liftfield[k]));
825 }
826 }
827
828 if (NUmodes != 0)
829 {
830 for (label k = 0; k < NUmodes; k++)
831 {
832 L_U_SUPmodes.append(tmp<volVectorField>(Umodes[k]));
833 }
834 }
835
836 if (ITHACAutilities::check_folder("./ITHACAoutput/Matrices/"))
837 {
838 word B_str = "B_" + name(liftfield.size()) + "_" + name(NUmodes) + "_" + name(
839 NSUPmodes);
840
841 if (ITHACAutilities::check_file("./ITHACAoutput/Matrices/" + B_str))
842 {
843 ITHACAstream::ReadDenseMatrix(B_matrix, "./ITHACAoutput/Matrices/", B_str);
844 }
845 else
846 {
847 B_matrix = diffusive_term(NUmodes, NPmodes, NSUPmodes);
848 }
849
850 word btStr = "bt_" + name(liftfield.size()) + "_" + name(NUmodes) + "_" + name(
851 NSUPmodes);
852
853 if (ITHACAutilities::check_file("./ITHACAoutput/Matrices/" + btStr))
854 {
855 ITHACAstream::ReadDenseMatrix(btMatrix, "./ITHACAoutput/Matrices/", btStr);
856 }
857 else
858 {
859 btMatrix = btTurbulence(NUmodes, NSUPmodes);
860 }
861
862 word K_str = "K_" + name(liftfield.size()) + "_" + name(NUmodes) + "_" + name(
863 NSUPmodes) + "_" + name(NPmodes);
864
865 if (ITHACAutilities::check_file("./ITHACAoutput/Matrices/" + K_str))
866 {
867 ITHACAstream::ReadDenseMatrix(K_matrix, "./ITHACAoutput/Matrices/", K_str);
868 }
869 else
870 {
871 K_matrix = pressure_gradient_term(NUmodes, NPmodes, NSUPmodes);
872 }
873
874 word M_str = "M_" + name(liftfield.size()) + "_" + name(NUmodes) + "_" + name(
875 NSUPmodes);
876
877 if (ITHACAutilities::check_file("./ITHACAoutput/Matrices/" + M_str))
878 {
879 ITHACAstream::ReadDenseMatrix(M_matrix, "./ITHACAoutput/Matrices/", M_str);
880 }
881 else
882 {
883 M_matrix = mass_term(NUmodes, NPmodes, NSUPmodes);
884 }
885
886 word D_str = "D_" + name(NPmodes);
887
888 if (ITHACAutilities::check_file("./ITHACAoutput/Matrices/" + D_str))
889 {
890 ITHACAstream::ReadDenseMatrix(D_matrix, "./ITHACAoutput/Matrices/", D_str);
891 }
892 else
893 {
894 D_matrix = laplacian_pressure(NPmodes);
895 }
896
897 word bc1_str = "BC1_" + name(liftfield.size()) + "_" + name(
898 NUmodes) + "_" + name(NPmodes);
899
900 if (ITHACAutilities::check_file("./ITHACAoutput/Matrices/" + bc1_str))
901 {
902 ITHACAstream::ReadDenseMatrix(BC1_matrix, "./ITHACAoutput/Matrices/", bc1_str);
903 }
904 else
905 {
906 BC1_matrix = pressure_BC1(NUmodes, NPmodes);
907 }
908
909 word bc2_str = "BC2_" + name(liftfield.size()) + "_" + name(
910 NUmodes) + "_" + name(
911 NSUPmodes) + "_" + name(NPmodes) + "_t";
912
913 if (ITHACAutilities::check_file("./ITHACAoutput/Matrices/" + bc2_str))
914 {
915 ITHACAstream::ReadDenseTensor(bc2Tensor, "./ITHACAoutput/Matrices/", bc2_str);
916 }
917 else
918 {
919 bc2Tensor = pressureBC2(NUmodes, NPmodes);
920 }
921
922 word bc3_str = "BC3_" + name(liftfield.size()) + "_" + name(
923 NUmodes) + "_" + name(NPmodes);
924
925 if (ITHACAutilities::check_file("./ITHACAoutput/Matrices/" + bc3_str))
926 {
927 ITHACAstream::ReadDenseMatrix(BC3_matrix, "./ITHACAoutput/Matrices/", bc3_str);
928 }
929 else
930 {
931 BC3_matrix = pressure_BC3(NUmodes, NPmodes);
932 }
933
934 word C_str = "C_" + name(liftfield.size()) + "_" + name(NUmodes) + "_" + name(
935 NSUPmodes) + "_t";
936
937 if (ITHACAutilities::check_file("./ITHACAoutput/Matrices/" + C_str))
938 {
939 ITHACAstream::ReadDenseTensor(C_tensor, "./ITHACAoutput/Matrices/", C_str);
940 }
941 else
942 {
943 C_tensor = convective_term_tens(NUmodes, NPmodes, NSUPmodes);
944 }
945
946 word ct1Str = "ct1_" + name(liftfield.size()) + "_" + name(
947 NUmodes) + "_" + name(
948 NSUPmodes) + "_" + name(nNutModes) + "_t";
949
950 if (ITHACAutilities::check_file("./ITHACAoutput/Matrices/" + ct1Str))
951 {
952 ITHACAstream::ReadDenseTensor(ct1Tensor, "./ITHACAoutput/Matrices/", ct1Str);
953 }
954 else
955 {
956 ct1Tensor = turbulenceTensor1(NUmodes, NSUPmodes, nNutModes);
957 }
958
959 word ct2Str = "ct2_" + name(liftfield.size()) + "_" + name(
960 NUmodes) + "_" + name(
961 NSUPmodes) + "_" + name(nNutModes) + "_t";
962
963 if (ITHACAutilities::check_file("./ITHACAoutput/Matrices/" + ct2Str))
964 {
965 ITHACAstream::ReadDenseTensor(ct2Tensor, "./ITHACAoutput/Matrices/", ct2Str);
966 }
967 else
968 {
969 ct2Tensor = turbulenceTensor2(NUmodes, NSUPmodes, nNutModes);
970 }
971
972 word G_str = "G_" + name(liftfield.size()) + "_" + name(NUmodes) + "_" + name(
973 NSUPmodes) + "_" + name(NPmodes) + "_t";
974
975 if (ITHACAutilities::check_file("./ITHACAoutput/Matrices/" + G_str))
976 {
977 ITHACAstream::ReadDenseTensor(gTensor, "./ITHACAoutput/Matrices/", G_str);
978 }
979 else
980 {
981 gTensor = divMomentum(NUmodes, NPmodes);
982 }
983
984 if (nutAve.size() != 0)
985 {
986 word ct1AveStr = "ct1Ave_" + name(liftfield.size()) + "_" + name(
987 NUmodes) + "_" + name(
988 NSUPmodes) + "_t";
989
990 if (ITHACAutilities::check_file("./ITHACAoutput/Matrices/" + ct1AveStr))
991 {
992 ITHACAstream::ReadDenseTensor(ct1AveTensor, "./ITHACAoutput/Matrices/",
993 ct1AveStr);
994 }
995 else
996 {
997 ct1AveTensor = turbulenceAveTensor1(NUmodes, NSUPmodes);
998 }
999
1000 word ct2AveStr = "ct2Ave_" + name(liftfield.size()) + "_" + name(
1001 NUmodes) + "_" + name(
1002 NSUPmodes) + "_t";
1003
1004 if (ITHACAutilities::check_file("./ITHACAoutput/Matrices/" + ct2AveStr))
1005 {
1006 ITHACAstream::ReadDenseTensor(ct2AveTensor, "./ITHACAoutput/Matrices/",
1007 ct2AveStr);
1008 }
1009 else
1010 {
1011 ct2AveTensor = turbulenceAveTensor2(NUmodes, NSUPmodes);
1012 }
1013 }
1014
1015 word ct1PPEStr = "ct1PPE_" + name(liftfield.size()) + "_" + name(
1016 NUmodes) + "_" + name(
1017 NSUPmodes) + "_" + name(NPmodes) + "_" + name(nNutModes) + "_t";
1018
1019 if (ITHACAutilities::check_file("./ITHACAoutput/Matrices/" + ct1PPEStr))
1020 {
1021 ITHACAstream::ReadDenseTensor(ct1PPETensor, "./ITHACAoutput/Matrices/",
1022 ct1PPEStr);
1023 }
1024 else
1025 {
1026 ct1PPETensor = turbulencePPETensor1(NUmodes, NSUPmodes, NPmodes, nNutModes);
1027 }
1028
1029 word ct2PPEStr = "ct2PPE_" + name(liftfield.size()) + "_" + name(
1030 NUmodes) + "_" + name(
1031 NSUPmodes) + "_" + name(NPmodes) + "_" + name(nNutModes) + "_t";
1032
1033 if (ITHACAutilities::check_file("./ITHACAoutput/Matrices/" + ct2PPEStr))
1034 {
1035 ITHACAstream::ReadDenseTensor(ct2PPETensor, "./ITHACAoutput/Matrices/",
1036 ct2PPEStr);
1037 }
1038 else
1039 {
1040 ct2PPETensor = turbulencePPETensor2(NUmodes, NSUPmodes, NPmodes, nNutModes);
1041 }
1042
1043 if (nutAve.size() != 0)
1044 {
1045 word ct1PPEAveStr = "ct1PPEAve_" + name(liftfield.size()) + "_" + name(
1046 NUmodes) + "_" + name(
1047 NSUPmodes) + "_" + name(NPmodes) + "_t";
1048
1049 if (ITHACAutilities::check_file("./ITHACAoutput/Matrices/" + ct1PPEAveStr))
1050 {
1051 ITHACAstream::ReadDenseTensor(ct1PPEAveTensor, "./ITHACAoutput/Matrices/",
1052 ct1PPEAveStr);
1053 }
1054 else
1055 {
1056 ct1PPEAveTensor = turbulencePPEAveTensor1(NUmodes, NSUPmodes, NPmodes);
1057 }
1058
1059 word ct2PPEAveStr = "ct2PPEAve_" + name(liftfield.size()) + "_" + name(
1060 NUmodes) + "_" + name(
1061 NSUPmodes) + "_" + name(NPmodes) + "_t";
1062
1063 if (ITHACAutilities::check_file("./ITHACAoutput/Matrices/" + ct2PPEAveStr))
1064 {
1065 ITHACAstream::ReadDenseTensor(ct2PPEAveTensor, "./ITHACAoutput/Matrices/",
1066 ct2PPEAveStr);
1067 }
1068 else
1069 {
1070 ct2PPEAveTensor = turbulencePPEAveTensor2(NUmodes, NSUPmodes, NPmodes);
1071 }
1072 }
1073
1074 if (bcMethod == "penalty")
1075 {
1076 bcVelVec = bcVelocityVec(NUmodes, NSUPmodes);
1077 bcVelMat = bcVelocityMat(NUmodes, NSUPmodes);
1078 }
1079 }
1080 else
1081 {
1082 B_matrix = diffusive_term(NUmodes, NPmodes, NSUPmodes);
1083 C_tensor = convective_term_tens(NUmodes, NPmodes, NSUPmodes);
1084 M_matrix = mass_term(NUmodes, NPmodes, NSUPmodes);
1085 K_matrix = pressure_gradient_term(NUmodes, NPmodes, NSUPmodes);
1086 D_matrix = laplacian_pressure(NPmodes);
1087 gTensor = divMomentum(NUmodes, NPmodes);
1088 BC1_matrix = pressure_BC1(NUmodes, NPmodes);
1089 bc2Tensor = pressureBC2(NUmodes, NPmodes);
1090 BC3_matrix = pressure_BC3(NUmodes, NPmodes);
1091 btMatrix = btTurbulence(NUmodes, NSUPmodes);
1092 ct1Tensor = turbulenceTensor1(NUmodes, NSUPmodes, nNutModes);
1093 ct2Tensor = turbulenceTensor2(NUmodes, NSUPmodes, nNutModes);
1094 ct1PPETensor = turbulencePPETensor1(NUmodes, NSUPmodes, NPmodes, nNutModes);
1095 ct2PPETensor = turbulencePPETensor2(NUmodes, NSUPmodes, NPmodes, nNutModes);
1096
1097 if (nutAve.size() != 0)
1098 {
1099 ct1AveTensor = turbulenceAveTensor1(NUmodes, NSUPmodes);
1100 ct2AveTensor = turbulenceAveTensor2(NUmodes, NSUPmodes);
1101 ct1PPEAveTensor = turbulencePPEAveTensor1(NUmodes, NSUPmodes, NPmodes);
1102 ct2PPEAveTensor = turbulencePPEAveTensor2(NUmodes, NSUPmodes, NPmodes);
1103 }
1104
1105 if (bcMethod == "penalty")
1106 {
1107 bcVelVec = bcVelocityVec(NUmodes, NSUPmodes);
1108 bcVelMat = bcVelocityMat(NUmodes, NSUPmodes);
1109 }
1110 }
1111
1112 // Export the matrices
1113 if (para->exportPython)
1114 {
1115 ITHACAstream::exportMatrix(B_matrix, "B", "python", "./ITHACAoutput/Matrices/");
1116 ITHACAstream::exportMatrix(K_matrix, "K", "python", "./ITHACAoutput/Matrices/");
1117 ITHACAstream::exportMatrix(D_matrix, "D", "python", "./ITHACAoutput/Matrices/");
1118 ITHACAstream::exportMatrix(M_matrix, "M", "python", "./ITHACAoutput/Matrices/");
1119 ITHACAstream::exportMatrix(BC1_matrix, "BC1", "python",
1120 "./ITHACAoutput/Matrices/");
1121 ITHACAstream::exportMatrix(BC3_matrix, "BC3", "python",
1122 "./ITHACAoutput/Matrices/");
1123 ITHACAstream::exportTensor(C_tensor, "C", "python", "./ITHACAoutput/Matrices/");
1124 ITHACAstream::exportTensor(gTensor, "G", "python", "./ITHACAoutput/Matrices/");
1125 ITHACAstream::exportTensor(bc2Tensor, "BC2", "python",
1126 "./ITHACAoutput/Matrices/");
1127 ITHACAstream::exportTensor(ct1Tensor, "ct1", "python",
1128 "./ITHACAoutput/Matrices/");
1129 ITHACAstream::exportTensor(ct2Tensor, "ct2", "python",
1130 "./ITHACAoutput/Matrices/");
1131 ITHACAstream::exportTensor(ct1PPETensor, "ct1PPE", "python",
1132 "./ITHACAoutput/Matrices/");
1133 ITHACAstream::exportTensor(ct2PPETensor, "ct2PPE", "python",
1134 "./ITHACAoutput/Matrices/");
1135
1136 if (nutAve.size() != 0)
1137 {
1138 ITHACAstream::exportTensor(ct1AveTensor, "ct1Ave", "python",
1139 "./ITHACAoutput/Matrices/");
1140 ITHACAstream::exportTensor(ct2AveTensor, "ct2Ave", "python",
1141 "./ITHACAoutput/Matrices/");
1142 ITHACAstream::exportTensor(ct1PPEAveTensor, "ct1PPEAve", "python",
1143 "./ITHACAoutput/Matrices/");
1144 ITHACAstream::exportTensor(ct2PPEAveTensor, "ct2PPEAve", "python",
1145 "./ITHACAoutput/Matrices/");
1146 }
1147 }
1148
1149 if (para->exportMatlab)
1150 {
1151 ITHACAstream::exportMatrix(B_matrix, "B", "matlab", "./ITHACAoutput/Matrices/");
1152 ITHACAstream::exportMatrix(K_matrix, "K", "matlab", "./ITHACAoutput/Matrices/");
1153 ITHACAstream::exportMatrix(D_matrix, "D", "matlab", "./ITHACAoutput/Matrices/");
1154 ITHACAstream::exportMatrix(M_matrix, "M", "matlab", "./ITHACAoutput/Matrices/");
1155 ITHACAstream::exportMatrix(BC1_matrix, "BC1", "matlab",
1156 "./ITHACAoutput/Matrices/");
1157 ITHACAstream::exportMatrix(BC3_matrix, "BC3", "matlab",
1158 "./ITHACAoutput/Matrices/");
1159 ITHACAstream::exportTensor(C_tensor, "C", "matlab", "./ITHACAoutput/Matrices/");
1160 ITHACAstream::exportTensor(gTensor, "G", "matlab", "./ITHACAoutput/Matrices/");
1161 ITHACAstream::exportTensor(bc2Tensor, "BC2", "matlab",
1162 "./ITHACAoutput/Matrices/");
1163 ITHACAstream::exportTensor(ct1Tensor, "ct1", "matlab",
1164 "./ITHACAoutput/Matrices/");
1165 ITHACAstream::exportTensor(ct2Tensor, "ct2", "matlab",
1166 "./ITHACAoutput/Matrices/");
1167 ITHACAstream::exportTensor(ct1PPETensor, "ct1PPE", "matlab",
1168 "./ITHACAoutput/Matrices/");
1169 ITHACAstream::exportTensor(ct2PPETensor, "ct2PPE", "matlab",
1170 "./ITHACAoutput/Matrices/");
1171
1172 if (nutAve.size() != 0)
1173 {
1174 ITHACAstream::exportTensor(ct1AveTensor, "ct1Ave", "matlab",
1175 "./ITHACAoutput/Matrices/");
1176 ITHACAstream::exportTensor(ct2AveTensor, "ct2Ave", "matlab",
1177 "./ITHACAoutput/Matrices/");
1178 ITHACAstream::exportTensor(ct1PPEAveTensor, "ct1PPEAve", "matlab",
1179 "./ITHACAoutput/Matrices/");
1180 ITHACAstream::exportTensor(ct2PPEAveTensor, "ct2PPEAve", "matlab",
1181 "./ITHACAoutput/Matrices/");
1182 }
1183 }
1184
1185 if (para->exportTxt)
1186 {
1187 ITHACAstream::exportMatrix(B_matrix, "B", "eigen", "./ITHACAoutput/Matrices/");
1188 ITHACAstream::exportMatrix(K_matrix, "K", "eigen", "./ITHACAoutput/Matrices/");
1189 ITHACAstream::exportMatrix(D_matrix, "D", "eigen", "./ITHACAoutput/Matrices/");
1190 ITHACAstream::exportMatrix(M_matrix, "M", "eigen", "./ITHACAoutput/Matrices/");
1191 ITHACAstream::exportMatrix(BC1_matrix, "BC1", "eigen",
1192 "./ITHACAoutput/Matrices/");
1193 ITHACAstream::exportMatrix(BC3_matrix, "BC3", "eigen",
1194 "./ITHACAoutput/Matrices/");
1195 ITHACAstream::exportTensor(C_tensor, "C", "eigen",
1196 "./ITHACAoutput/Matrices/C");
1197 ITHACAstream::exportTensor(gTensor, "G", "eigen",
1198 "./ITHACAoutput/Matrices/G");
1199 ITHACAstream::exportTensor(bc2Tensor, "BC2_", "eigen",
1200 "./ITHACAoutput/Matrices/BC2");
1201 ITHACAstream::exportMatrix(btMatrix, "bt", "eigen", "./ITHACAoutput/Matrices/");
1202 ITHACAstream::exportTensor(ct1Tensor, "ct1_", "eigen",
1203 "./ITHACAoutput/Matrices/ct1");
1204 ITHACAstream::exportTensor(ct2Tensor, "ct2_", "eigen",
1205 "./ITHACAoutput/Matrices/ct2");
1206 ITHACAstream::exportTensor(ct1PPETensor, "ct1PPE_", "eigen",
1207 "./ITHACAoutput/Matrices/ct1PPE");
1208 ITHACAstream::exportTensor(ct2PPETensor, "ct2PPE_", "eigen",
1209 "./ITHACAoutput/Matrices/ct2PPE");
1210
1211 if (nutAve.size() != 0)
1212 {
1213 ITHACAstream::exportTensor(ct1AveTensor, "ct1Ave_", "eigen",
1214 "./ITHACAoutput/Matrices/ct1Ave");
1215 ITHACAstream::exportTensor(ct2AveTensor, "ct2Ave_", "eigen",
1216 "./ITHACAoutput/Matrices/ct2Ave");
1217 ITHACAstream::exportTensor(ct1PPEAveTensor, "ct1PPEAve_", "eigen",
1218 "./ITHACAoutput/Matrices/ct1PPEAve");
1219 ITHACAstream::exportTensor(ct2PPEAveTensor, "ct2PPEAve_", "eigen",
1220 "./ITHACAoutput/Matrices/ct2PPEAve");
1221 }
1222 }
1223
1224 bTotalMatrix = B_matrix + btMatrix;
1225 label cSize = NUmodes + NSUPmodes + liftfield.size();
1226 cTotalTensor.resize(cSize, nNutModes, cSize);
1227 cTotalTensor = ct1Tensor + ct2Tensor;
1228 cTotalPPETensor.resize(NPmodes, nNutModes, cSize);
1229 cTotalPPETensor = ct1PPETensor + ct2PPETensor;
1230
1231 if (nutAve.size() != 0)
1232 {
1233 cTotalAveTensor.resize(cSize, nutAve.size(), cSize);
1234 cTotalAveTensor = ct1AveTensor + ct2AveTensor;
1235 cTotalPPEAveTensor.resize(NPmodes, nutAve.size(), cSize);
1236 cTotalPPEAveTensor = ct1PPEAveTensor + ct2PPEAveTensor;
1237 }
1238
1239 if (rbfInterp == true && (!Pstream::parRun()))
1240 {
1241 if (ITHACAutilities::check_file("./radii.txt"))
1242 {
1243 radii = ITHACAstream::readMatrix("./radii.txt");
1244 M_Assert(radii.size() == nNutModes,
1245 "Thes size of the shape parameters vector must be equal to the number of eddy viscosity modes nNutModes");
1246 }
1247 else
1248 {
1249 radii = Eigen::MatrixXd::Ones(nNutModes,
1250 1) * e;
1251 }
1252
1253 samples.resize(nNutModes);
1254 rbfSplines.resize(nNutModes);
1255 Eigen::MatrixXd weights;
1256
1257 for (label i = 0; i < nNutModes; i++)
1258 {
1259 word weightName = "wRBF_N" + name(i + 1) + "_" + name(liftfield.size()) + "_"
1260 + name(NUmodes) + "_" + name(NSUPmodes) ;
1261
1262 if (ITHACAutilities::check_file("./ITHACAoutput/weightsPPE/" + weightName))
1263 {
1264 samples[i] = new SPLINTER::DataTable(1, 1);
1265
1266 for (label j = 0; j < coeffL2.cols(); j++)
1267 {
1268 samples[i]->addSample(velRBF.row(j), coeffL2(i, j));
1269 }
1270
1271 ITHACAstream::ReadDenseMatrix(weights, "./ITHACAoutput/weightsPPE/",
1272 weightName);
1273 rbfSplines[i] = new SPLINTER::RBFSpline(*samples[i],
1274 SPLINTER::RadialBasisFunctionType::GAUSSIAN, weights, radii(i));
1275 std::cout << "Constructing RadialBasisFunction for mode " << i + 1 << std::endl;
1276 }
1277 else
1278 {
1279 samples[i] = new SPLINTER::DataTable(1, 1);
1280
1281 for (label j = 0; j < coeffL2.cols(); j++)
1282 {
1283 samples[i]->addSample(velRBF.row(j), coeffL2(i, j));
1284 }
1285
1286 rbfSplines[i] = new SPLINTER::RBFSpline(*samples[i],
1287 SPLINTER::RadialBasisFunctionType::GAUSSIAN, false, radii(i));
1288 ITHACAstream::SaveDenseMatrix(rbfSplines[i]->weights,
1289 "./ITHACAoutput/weightsPPE/", weightName);
1290 std::cout << "Constructing RadialBasisFunction for mode " << i + 1 << std::endl;
1291 }
1292 }
1293 }
1294}
1295
1296List < Eigen::MatrixXd > UnsteadyNSTurb::velDerivativeCoeff(Eigen::MatrixXd A,
1297 Eigen::MatrixXd G,
1298 Eigen::VectorXd initSnapInd, Eigen::VectorXd timeSnap)
1299{
1300 List < Eigen::MatrixXd > newCoeffs;
1301 newCoeffs.setSize(2);
1302 label velCoeffsNum = A.cols();
1303 label snapshotsNum = A.rows();
1304 Eigen::MatrixXd pars;
1305 label parsSamplesNum = initSnapInd.size();
1306 label timeSnapshotsPerSample = snapshotsNum / parsSamplesNum;
1307 label newColsNum = 2 * velCoeffsNum;
1308 label newRowsNum = snapshotsNum - parsSamplesNum;
1309 newCoeffs[0].resize(newRowsNum, newColsNum);
1310 newCoeffs[1].resize(newRowsNum, G.cols());
1311
1312 for (label j = 0; j < parsSamplesNum; j++)
1313 {
1314 Eigen::MatrixXd b0 = A.middleRows(j * timeSnapshotsPerSample,
1315 timeSnapshotsPerSample - 1);
1316 Eigen::MatrixXd b2 = A.middleRows(j * timeSnapshotsPerSample + 1,
1317 timeSnapshotsPerSample - 1);
1318 Eigen::MatrixXd bNew(b0.rows(), b0.cols() + b2.cols());
1319 bNew << b2, (b2 - b0) / (timeSnap(j, 0));
1320 newCoeffs[0].block(j * timeSnapshotsPerSample - j, 0,
1321 timeSnapshotsPerSample - 1, newColsNum) = bNew;
1322 newCoeffs[1].middleRows(j * timeSnapshotsPerSample - j,
1323 timeSnapshotsPerSample - 1) = G.middleRows(j * timeSnapshotsPerSample + 1,
1324 timeSnapshotsPerSample - 1);
1325 }
1326
1327 interChoice = 3;
1328 return newCoeffs;
1329}
1330
1331List < Eigen::MatrixXd > UnsteadyNSTurb::velParCoeff(Eigen::MatrixXd A,
1332 Eigen::MatrixXd G)
1333{
1334 List < Eigen::MatrixXd > newCoeffs;
1335 newCoeffs.setSize(2);
1336 Eigen::MatrixXd pars;
1337 pars = z.leftCols(z.cols() - 1);
1338 newCoeffs[0].resize(A.rows(), A.cols() + z.cols() - 1);
1339 newCoeffs[1].resize(G.rows(), G.cols());
1340 newCoeffs[0] << pars, A;
1341 newCoeffs[1] = G;
1342 interChoice = 2;
1343 return newCoeffs;
1344}
1345
1346List < Eigen::MatrixXd > UnsteadyNSTurb::velParDerivativeCoeff(
1347 Eigen::MatrixXd A, Eigen::MatrixXd G,
1348 Eigen::VectorXd initSnapInd, Eigen::VectorXd timeSnap)
1349{
1350 List < Eigen::MatrixXd > newCoeffs;
1351 newCoeffs.setSize(2);
1352 label velCoeffsNum = A.cols();
1353 label snapshotsNum = A.rows();
1354 Eigen::MatrixXd pars;
1355 pars = z.leftCols(z.cols() - 1);
1356 label parsSamplesNum = initSnapInd.size();
1357 label timeSnapshotsPerSample = snapshotsNum / parsSamplesNum;
1358 label newColsNum = 2 * velCoeffsNum;
1359 label newRowsNum = snapshotsNum - parsSamplesNum;
1360 newCoeffs[0].resize(newRowsNum, newColsNum + z.cols() - 1);
1361 newCoeffs[1].resize(newRowsNum, G.cols());
1362
1363 for (label j = 0; j < parsSamplesNum; j++)
1364 {
1365 Eigen::MatrixXd b0 = A.middleRows(j * timeSnapshotsPerSample,
1366 timeSnapshotsPerSample - 1);
1367 Eigen::MatrixXd b2 = A.middleRows(j * timeSnapshotsPerSample + 1,
1368 timeSnapshotsPerSample - 1);
1369 Eigen::MatrixXd bNew(b0.rows(), b0.cols() + b2.cols());
1370 bNew << b2, (b2 - b0) / (timeSnap(j, 0));
1371 newCoeffs[0].block(j * timeSnapshotsPerSample - j, 0,
1372 timeSnapshotsPerSample - 1, z.cols() - 1) = pars.middleRows(
1373 j * timeSnapshotsPerSample + 1,
1374 timeSnapshotsPerSample - 1);
1375 newCoeffs[0].block(j * timeSnapshotsPerSample - j, z.cols() - 1,
1376 timeSnapshotsPerSample - 1, newColsNum) = bNew;
1377 newCoeffs[1].middleRows(j * timeSnapshotsPerSample - j,
1378 timeSnapshotsPerSample - 1) = G.middleRows(j * timeSnapshotsPerSample + 1,
1379 timeSnapshotsPerSample - 1);
1380 }
1381
1382 interChoice = 4;
1383 return newCoeffs;
1384}
1385
1386Eigen::MatrixXd UnsteadyNSTurb::velParDerivativeCoeff(Eigen::MatrixXd A,
1387 Eigen::VectorXd par, double timeSnap)
1388{
1389 Eigen::MatrixXd newCoeffs;
1390 label velCoeffsNum = A.cols();
1391 label snapshotsNum = A.rows();
1392 label parsSamplesNum = par.size();
1393 label newColsNum = 2 * velCoeffsNum + parsSamplesNum;
1394 label newRowsNum = snapshotsNum - 1;
1395 newCoeffs.resize(newRowsNum, newColsNum);
1396 Eigen::MatrixXd b0 = A.topRows(A.rows() - 1);
1397 Eigen::MatrixXd b1 = A.bottomRows(A.rows() - 1);
1398 Eigen::MatrixXd bNew(b0.rows(), b0.cols() + b1.cols());
1399 bNew << b1, ((b1 - b0) / (timeSnap));
1400 newCoeffs.leftCols(parsSamplesNum) = Eigen::MatrixXd::Ones(newRowsNum,
1401 parsSamplesNum) * par;
1402 newCoeffs.rightCols(newColsNum - parsSamplesNum) = bNew;
1403 return newCoeffs;
1404}
mesh
Foam::fvMesh & mesh
Definition createMesh.H:47
runTime
Foam::Time & runTime
Definition createTime.H:33
createMesh.H
createTime.H
M_Assert
#define M_Assert(Expr, Msg)
Definition ITHACAassert.H:46
UnsteadyNSTurb.H
Header file of the UnsteadyNSTurb class.
ITHACAparameters::ITHACAdict
IOdictionary * ITHACAdict
Dictionary for input objects from file.
Definition ITHACAparameters.H:64
ITHACAparameters::getInstance
static ITHACAparameters * getInstance()
Gets an instance of ITHACAparameters, to be used if the instance is already existing.
Definition ITHACAparameters.C:55
UnsteadyNSTurb::ct1AveTensor
Eigen::Tensor< double, 3 > ct1AveTensor
Turbulent average viscosity tensor for the splitting approach.
Definition UnsteadyNSTurb.H:104
UnsteadyNSTurb::turbulenceAveTensor2
Eigen::Tensor< double, 3 > turbulenceAveTensor2(label NUmodes, label NSUPmodes)
ct2Ave added tensor for approximation of the averaged part of the eddy viscosity
Definition UnsteadyNSTurb.C:356
UnsteadyNSTurb::cTotalTensor
Eigen::Tensor< double, 3 > cTotalTensor
Turbulent total viscosity tensor.
Definition UnsteadyNSTurb.H:125
UnsteadyNSTurb::nutAve
PtrList< volScalarField > nutAve
List of for eddy viscosity time-averaged fields.
Definition UnsteadyNSTurb.H:86
UnsteadyNSTurb::_pRefValue
scalar _pRefValue
Pressure reference value.
Definition UnsteadyNSTurb.H:167
UnsteadyNSTurb::btTurbulence
Eigen::MatrixXd btTurbulence(label NUmodes, label NSUPmodes)
bt added matrix for the turbulence treatement
Definition UnsteadyNSTurb.C:448
UnsteadyNSTurb::velParDerivativeCoeff
List< Eigen::MatrixXd > velParDerivativeCoeff(Eigen::MatrixXd A, Eigen::MatrixXd G, Eigen::VectorXd initSnapInd, Eigen::VectorXd timeSnap)
A method to compute the two matrices needed for the RBF interpolation by combining the parameter valu...
Definition UnsteadyNSTurb.C:1346
UnsteadyNSTurb::turbulencePPETensor1
Eigen::Tensor< double, 3 > turbulencePPETensor1(label NUmodes, label NSUPmodes, label NPmodes, label nNutModes)
ct1PPE added tensor for the turbulence treatement in the PPE method
Definition UnsteadyNSTurb.C:265
UnsteadyNSTurb::rbfSplines
std::vector< SPLINTER::RBFSpline * > rbfSplines
Create a samples for interpolation.
Definition UnsteadyNSTurb.H:92
UnsteadyNSTurb::turbulenceAveTensor1
Eigen::Tensor< double, 3 > turbulenceAveTensor1(label NUmodes, label NSUPmodes)
ct1Ave added tensor for approximation of the averaged part of the eddy viscosity
Definition UnsteadyNSTurb.C:238
UnsteadyNSTurb::projectSUP
void projectSUP(fileName folder, label NUmodes, label NPmodes, label NSUPmodes, label nNutModes, bool rbfInterp=true)
Project using a supremizer approach.
Definition UnsteadyNSTurb.C:471
UnsteadyNSTurb::turbulencePPEAveTensor1
Eigen::Tensor< double, 3 > turbulencePPEAveTensor1(label NUmodes, label NSUPmodes, label NPmodes)
ct1PPEAve added tensor for approximation of the averaged part of the eddy viscosity with the usage of...
Definition UnsteadyNSTurb.C:297
UnsteadyNSTurb::UnsteadyNSTurb
UnsteadyNSTurb()
Construct Null.
Definition UnsteadyNSTurb.C:39
UnsteadyNSTurb::turbulenceTensor2
Eigen::Tensor< double, 3 > turbulenceTensor2(label NUmodes, label NSUPmodes, label nNutModes)
ct2 added tensor for the turbulence treatement
Definition UnsteadyNSTurb.C:330
UnsteadyNSTurb::_nut
autoPtr< volScalarField > _nut
Eddy viscosity field.
Definition UnsteadyNSTurb.H:158
UnsteadyNSTurb::ct2PPEAveTensor
Eigen::Tensor< double, 3 > ct2PPEAveTensor
Turbulent average viscosity tensor for the splitting approach in the PPE equation.
Definition UnsteadyNSTurb.H:122
UnsteadyNSTurb::interChoice
label interChoice
Interpolation independent variable choice.
Definition UnsteadyNSTurb.H:155
UnsteadyNSTurb::coeffL2
Eigen::MatrixXd coeffL2
The matrix of L2 projection coefficients for the eddy viscosity.
Definition UnsteadyNSTurb.H:140
UnsteadyNSTurb::btMatrix
Eigen::MatrixXd btMatrix
Turbulent viscosity term.
Definition UnsteadyNSTurb.H:95
UnsteadyNSTurb::nutModes
volScalarModes nutModes
List of POD modes for eddy viscosity.
Definition UnsteadyNSTurb.H:83
UnsteadyNSTurb::turbulencePPEAveTensor2
Eigen::Tensor< double, 3 > turbulencePPEAveTensor2(label NUmodes, label NSUPmodes, label NPmodes)
ct2PPEAve added tensor for approximation of the averaged part of the eddy viscosity with the usage of...
Definition UnsteadyNSTurb.C:414
UnsteadyNSTurb::cTotalPPEAveTensor
Eigen::Tensor< double, 3 > cTotalPPEAveTensor
Turbulent total average viscosity tensor for the splitting approach in the PPE equation.
Definition UnsteadyNSTurb.H:134
UnsteadyNSTurb::e
double e
RBF functions radius.
Definition UnsteadyNSTurb.H:149
UnsteadyNSTurb::nutFields
PtrList< volScalarField > nutFields
List of snapshots for the solution for eddy viscosity.
Definition UnsteadyNSTurb.H:80
UnsteadyNSTurb::ct1PPEAveTensor
Eigen::Tensor< double, 3 > ct1PPEAveTensor
Turbulent average viscosity tensor for the splitting approach in the PPE equation.
Definition UnsteadyNSTurb.H:110
UnsteadyNSTurb::ct2PPETensor
Eigen::Tensor< double, 3 > ct2PPETensor
Turbulent viscosity tensor in the PPE equation.
Definition UnsteadyNSTurb.H:119
UnsteadyNSTurb::ct1PPETensor
Eigen::Tensor< double, 3 > ct1PPETensor
Turbulent viscosity tensor in the PPE equation.
Definition UnsteadyNSTurb.H:107
UnsteadyNSTurb::ct2Tensor
Eigen::Tensor< double, 3 > ct2Tensor
Turbulent viscosity tensor.
Definition UnsteadyNSTurb.H:113
UnsteadyNSTurb::ct2AveTensor
Eigen::Tensor< double, 3 > ct2AveTensor
Turbulent average viscosity tensor for the splitting approach.
Definition UnsteadyNSTurb.H:116
UnsteadyNSTurb::_pRefCell
label _pRefCell
Pressure reference cell.
Definition UnsteadyNSTurb.H:164
UnsteadyNSTurb::ct1Tensor
Eigen::Tensor< double, 3 > ct1Tensor
Turbulent viscosity tensor.
Definition UnsteadyNSTurb.H:101
UnsteadyNSTurb::cTotalPPETensor
Eigen::Tensor< double, 3 > cTotalPPETensor
Turbulent total viscosity tensor in the PPE equation.
Definition UnsteadyNSTurb.H:131
UnsteadyNSTurb::velParCoeff
List< Eigen::MatrixXd > velParCoeff(Eigen::MatrixXd A, Eigen::MatrixXd G)
A method to compute the two matrices needed for the RBF interpolation by combining the parameter samp...
Definition UnsteadyNSTurb.C:1331
UnsteadyNSTurb::bTotalMatrix
Eigen::MatrixXd bTotalMatrix
Total B Matrix.
Definition UnsteadyNSTurb.H:137
UnsteadyNSTurb::velDerivativeCoeff
List< Eigen::MatrixXd > velDerivativeCoeff(Eigen::MatrixXd A, Eigen::MatrixXd G, Eigen::VectorXd initSnapInd, Eigen::VectorXd timeSnap)
A method to compute the two matrices needed for the RBF interpolation by combining the velocity L2 pr...
Definition UnsteadyNSTurb.C:1296
UnsteadyNSTurb::nNutModes
label nNutModes
Number of viscoisty modes used for the projection.
Definition UnsteadyNSTurb.H:152
UnsteadyNSTurb::cTotalAveTensor
Eigen::Tensor< double, 3 > cTotalAveTensor
Turbulent total average viscosity tensor for the splitting approach.
Definition UnsteadyNSTurb.H:128
UnsteadyNSTurb::projectPPE
void projectPPE(fileName folder, label NUmodes, label NPmodes, label NSUPmodes, label nNutModes, bool rbfInterp=true)
Project using the Poisson Equation for pressure.
Definition UnsteadyNSTurb.C:811
UnsteadyNSTurb::radii
Eigen::VectorXd radii
RBF shape parameters vector.
Definition UnsteadyNSTurb.H:146
UnsteadyNSTurb::turbulencePPETensor2
Eigen::Tensor< double, 3 > turbulencePPETensor2(label NUmodes, label NSUPmodes, label NPmodes, label nNutModes)
ct2PPE added tensor for the turbulence treatement in the PPE method
Definition UnsteadyNSTurb.C:383
UnsteadyNSTurb::z
Eigen::MatrixXd z
Time-parameter combined matrix.
Definition UnsteadyNSTurb.H:98
UnsteadyNSTurb::samples
std::vector< SPLINTER::DataTable * > samples
Create a Rbf splines for interpolation.
Definition UnsteadyNSTurb.H:89
UnsteadyNSTurb::velRBF
Eigen::MatrixXd velRBF
Velocity coefficients for RBF interpolation.
Definition UnsteadyNSTurb.H:143
UnsteadyNSTurb::turbulenceTensor1
Eigen::Tensor< double, 3 > turbulenceTensor1(label NUmodes, label NSUPmodes, label nNutModes)
ct1 added tensor for the turbulence treatement
Definition UnsteadyNSTurb.C:211
UnsteadyProblem::checkWrite
bool checkWrite(Time &timeObject)
Function to check if the solution must be exported.
Definition UnsteadyProblem.C:37
UnsteadyProblem::startTime
scalar startTime
Start Time (initial time to start storing the snapshots)
Definition UnsteadyProblem.H:43
UnsteadyProblem::writeEvery
scalar writeEvery
Time step of the writing procedure.
Definition UnsteadyProblem.H:52
UnsteadyProblem::timeStep
scalar timeStep
Time step of the simulation.
Definition UnsteadyProblem.H:49
UnsteadyProblem::nextWrite
scalar nextWrite
Auxiliary variable to store the next writing instant.
Definition UnsteadyProblem.H:55
UnsteadyProblem::finalTime
scalar finalTime
Final time (final time of the simulation and consequently of the acquisition of the snapshots)
Definition UnsteadyProblem.H:46
reductionProblem::writeMu
void writeMu(List< scalar > mu_now)
Write out a list of scalar corresponding to the parameters used in the truthSolve.
Definition reductionProblem.C:242
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::counter
label counter
Counter used for the output of the full order solutions.
Definition reductionProblem.H:110
reductionProblem::offline
bool offline
Boolean variable, it is 1 if the Offline phase has already been computed, else 0.
Definition reductionProblem.H:99
reductionProblem::ITHACAdict
IOdictionary * ITHACAdict
dictionary to store input output infos
Definition reductionProblem.H:101
reductionProblem::mu
Eigen::MatrixXd mu
Row matrix of parameters.
Definition reductionProblem.H:89
reductionProblem::_args
autoPtr< argList > _args
argList
Definition reductionProblem.H:103
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::bcVelVec
List< Eigen::MatrixXd > bcVelVec
Boundary term for penalty method - vector.
Definition steadyNS.H:218
steadyNS::maxIter
scalar maxIter
Number of maximum iterations to be done for the computation of the truth solution.
Definition steadyNS.H:125
steadyNS::NPmodes
label NPmodes
Number of pressure modes used for the projection.
Definition steadyNS.H:143
steadyNS::supex
bool supex
Boolean variable to check the existence of the supremizer modes.
Definition steadyNS.H:256
steadyNS::BC1_matrix
Eigen::MatrixXd BC1_matrix
PPE BC1.
Definition steadyNS.H:182
steadyNS::diffusive_term
Eigen::MatrixXd diffusive_term(label NUmodes, label NPmodes, label NSUPmodes)
Diffusive Term.
Definition steadyNS.C:861
steadyNS::_phi
autoPtr< surfaceScalarField > _phi
Flux.
Definition steadyNS.H:293
steadyNS::BC3_matrix
Eigen::MatrixXd BC3_matrix
PPE BC3.
Definition steadyNS.H:191
steadyNS::C_tensor
Eigen::Tensor< double, 3 > C_tensor
Diffusion term.
Definition steadyNS.H:167
steadyNS::Pfield
PtrList< volScalarField > Pfield
List of pointers used to form the pressure snapshots matrix.
Definition steadyNS.H:86
steadyNS::supmodes
volVectorModes supmodes
List of pointers used to form the supremizer modes.
Definition steadyNS.H:104
steadyNS::bcVelocityVec
List< Eigen::MatrixXd > bcVelocityVec(label NUmodes, label NSUPmodes)
Boundary integral modes on boundary used by the penaly method.
Definition steadyNS.C:1418
steadyNS::_fvOptions
autoPtr< fv::options > _fvOptions
fvOptions
Definition steadyNS.H:287
steadyNS::_runTime
autoPtr< Time > _runTime
Time.
Definition steadyNS.H:290
steadyNS::Umodes
volVectorModes Umodes
List of pointers used to form the velocity modes.
Definition steadyNS.H:101
steadyNS::Ufield
PtrList< volVectorField > Ufield
List of pointers used to form the velocity snapshots matrix.
Definition steadyNS.H:89
steadyNS::pressure_BC1
Eigen::MatrixXd pressure_BC1(label NPmodes, label NUmodes)
Term N° 1 given by the additional boundary condition using a PPE approach.
Definition steadyNS.C:1212
steadyNS::divergence_term
Eigen::MatrixXd divergence_term(label NUmodes, label NPmodes, label NSUPmodes)
Divergence Term (supremizer approach)
Definition steadyNS.C:1073
steadyNS::tolerance
scalar tolerance
Tolerance for the residual of the stationary problems, there is the same tolerance for velocity and p...
Definition steadyNS.H:122
steadyNS::mass_term
Eigen::MatrixXd mass_term(label NUmodes, label NPmodes, label NSUPmodes)
Mass Term.
Definition steadyNS.C:1041
steadyNS::_mesh
autoPtr< fvMesh > _mesh
Mesh.
Definition steadyNS.H:281
steadyNS::_laminarTransport
autoPtr< singlePhaseTransportModel > _laminarTransport
Laminar transport (used by turbulence model)
Definition steadyNS.H:302
steadyNS::NUmodes
label NUmodes
Number of velocity modes used for the projection.
Definition steadyNS.H:140
steadyNS::liftfield
PtrList< volVectorField > liftfield
List of pointers used to form the list of lifting functions.
Definition steadyNS.H:110
steadyNS::pRefValue
scalar pRefValue
Reference pressure value.
Definition steadyNS.H:311
steadyNS::pRefCell
label pRefCell
Reference pressure cell.
Definition steadyNS.H:308
steadyNS::L_U_SUPmodes
volVectorModes L_U_SUPmodes
List of pointers containing the total number of lift, supremizer and velocity modes.
Definition steadyNS.H:116
steadyNS::para
ITHACAparameters * para
Definition steadyNS.H:82
steadyNS::pressure_BC3
Eigen::MatrixXd pressure_BC3(label NPmodes, label NUmodes)
Term N° 3 given by the additional boundary condition using a PPE approach.
Definition steadyNS.C:1338
steadyNS::NNutModesOut
label NNutModesOut
Number of nut modes to be calculated.
Definition steadyNS.H:137
steadyNS::B_matrix
Eigen::MatrixXd B_matrix
Diffusion term.
Definition steadyNS.H:157
steadyNS::D_matrix
Eigen::MatrixXd D_matrix
Laplacian term PPE.
Definition steadyNS.H:173
steadyNS::_MRF
autoPtr< IOMRFZoneList > _MRF
MRF variable.
Definition steadyNS.H:305
steadyNS::NSUPmodes
label NSUPmodes
Number of supremizer modes used for the projection.
Definition steadyNS.H:146
steadyNS::gTensor
Eigen::Tensor< double, 3 > gTensor
Divergence of momentum PPE.
Definition steadyNS.H:179
steadyNS::NNutModes
label NNutModes
Number of nut modes used for the projection.
Definition steadyNS.H:149
steadyNS::K_matrix
Eigen::MatrixXd K_matrix
Gradient of pressure matrix.
Definition steadyNS.H:163
steadyNS::convective_term_tens
Eigen::Tensor< double, 3 > convective_term_tens(label NUmodes, label NPmodes, label NSUPmodes)
Export convective term as a tensor.
Definition steadyNS.C:995
steadyNS::NUmodesOut
label NUmodesOut
Number of velocity modes to be calculated.
Definition steadyNS.H:128
steadyNS::divMomentum
Eigen::Tensor< double, 3 > divMomentum(label NUmodes, label NPmodes)
Divergence of convective term (PPE approach)
Definition steadyNS.C:1147
steadyNS::pressure_gradient_term
Eigen::MatrixXd pressure_gradient_term(label NUmodes, label NPmodes, label NSUPmodes)
Gradient of pressure.
Definition steadyNS.C:923
steadyNS::bcVelMat
List< Eigen::MatrixXd > bcVelMat
Boundary term for penalty method - matrix.
Definition steadyNS.H:221
steadyNS::P_matrix
Eigen::MatrixXd P_matrix
Div of velocity.
Definition steadyNS.H:170
steadyNS::M_matrix
Eigen::MatrixXd M_matrix
Mass Matrix.
Definition steadyNS.H:160
steadyNS::laplacian_pressure
Eigen::MatrixXd laplacian_pressure(label NPmodes)
Laplacian of pressure term (PPE approach)
Definition steadyNS.C:1183
steadyNS::NPmodesOut
label NPmodesOut
Number of pressure modes to be calculated.
Definition steadyNS.H:131
steadyNS::bcVelocityMat
List< Eigen::MatrixXd > bcVelocityMat(label NUmodes, label NSUPmodes)
Boundary integral modes on boundary used by the penaly method.
Definition steadyNS.C:1455
steadyNS::bc2Tensor
Eigen::Tensor< double, 3 > bc2Tensor
PPE BC2.
Definition steadyNS.H:188
steadyNS::_U
autoPtr< volVectorField > _U
Velocity field.
Definition steadyNS.H:263
steadyNS::pressureBC2
Eigen::Tensor< double, 3 > pressureBC2(label NPmodes, label NUmodes)
Term N° 2 given by the additional boundary condition using a PPE approach.
Definition steadyNS.C:1293
steadyNS::Pmodes
volScalarModes Pmodes
List of pointers used to form the pressure modes.
Definition steadyNS.H:98
steadyNS::bcMethod
word bcMethod
Boundary Method.
Definition steadyNS.H:317
steadyNS::_p
autoPtr< volScalarField > _p
Pressure field.
Definition steadyNS.H:260
unsteadyNS::turbulence
autoPtr< incompressible::turbulenceModel > turbulence
Turbulence model.
Definition unsteadyNS.H:76
unsteadyNS::timeDerivativeSchemeOrder
word timeDerivativeSchemeOrder
Definition unsteadyNS.H:99
unsteadyNS::_pimple
autoPtr< pimpleControl > _pimple
pimpleControl
Definition unsteadyNS.H:73
T
volScalarField & T
Definition createT.H:46
A
volScalarField & A
Definition createConstants.H:144
ITHACAstream::ReadDenseMatrix
void ReadDenseMatrix(MatrixType &Matrix, word folder, word MatrixName)
Read a dense matrix from a binary format file.
Definition ITHACAstream.H:670
ITHACAstream::exportSolution
void exportSolution(GeometricField< Type, PatchField, GeoMesh > &s, fileName subfolder, fileName folder, word fieldName)
Export a field to file in a certain folder and subfolder.
Definition ITHACAstream.C:849
ITHACAstream::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::SaveDenseMatrix
void SaveDenseMatrix(MatrixType &Matrix, word folder, word MatrixName)
Save a dense matrix to a binary format file.
Definition ITHACAstream.H:655
ITHACAstream::ReadDenseTensor
void ReadDenseTensor(TensorType &Tensor, word folder, word MatrixName)
Read a dense tensor from file.
Definition ITHACAstream.H:717
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::exportTensor
void exportTensor(Eigen::Tensor< T, 3 > tensor, word Name, word type, word folder)
Export the reduced tensor in numpy (tipo=python), matlab (tipo=matlab) and txt (tipo=eigen) format.
Definition ITHACAstream.C:264
ITHACAstream::SaveDenseTensor
void SaveDenseTensor(TensorType &Tensor, word folder, word MatrixName)
Save a dense tensor to file.
Definition ITHACAstream.H:696
ITHACAutilities::getCoeffs
Eigen::VectorXd getCoeffs(GeometricField< Type, PatchField, GeoMesh > &snapshot, PtrList< GeometricField< Type, PatchField, GeoMesh > > &modes, label Nmodes, bool consider_volumes)
Projects a snapshot on a basis function and gets the coefficients of the projection.
Definition ITHACAcoeffsMass.C:317
ITHACAutilities::check_pod
bool check_pod()
Check if the POD data folder "./ITHACAoutput/POD" exists.
Definition ITHACAsystem.C:143
ITHACAutilities::check_off
bool check_off()
Check if the offline data folder "./ITHACAoutput/Offline" exists.
Definition ITHACAsystem.C:107
ITHACAutilities::check_folder
bool check_folder(word folder)
Checks if a folder exists.
Definition ITHACAsystem.C:89
ITHACAutilities::check_file
bool check_file(std::string fileName)
Function that returns true if a file exists.
Definition ITHACAsystem.C:136
ITHACAutilities::check_sup
bool check_sup()
Check if the supremizer folder exists.
Definition ITHACAsystem.C:173
phi
surfaceScalarField & phi
Definition createFields.H:68
U
volVectorField & U
Definition createFields.H:42
pimple
pimpleControl & pimple
Definition createFields.H:334
nut
volScalarField & nut
Definition createFields.H:72
UEqn.H
createFields.H
createFvOptions.H
pEqn.H
p
volScalarField & p
Definition createFields.H:41
laminarTransport
singlePhaseTransportModel & laminarTransport
Definition createFields.H:86
i
label i
Definition pEqn.H:46
Generated by   doxygen 1.11.0