Loading...
Searching...
No Matches
UnsteadyNSTurb.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-------------------------------------------------------------------------------
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
33#include <cmath>
34#include <cstdlib>
35#include <fstream>
36#include <iostream>
37
38// * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * * * //
39
40UnsteadyNSTurb::UnsteadyNSTurb() {}
41
42UnsteadyNSTurb::UnsteadyNSTurb(int argc, char* argv[])
43{
44 _args = autoPtr<argList>(new argList(argc, argv));
45
46 if (!_args->checkRootCase())
47 {
48 Foam::FatalError.exit();
49 }
50
51 argList& args = _args();
52#include "createTime.H"
53#include "createMesh.H"
54 _pimple = autoPtr<pimpleControl>(new pimpleControl(mesh));
55#include "createFields.H"
56#include "createUfIfPresent.H"
57#include "createFvOptions.H"
58 ITHACAdict = new IOdictionary
59 (
60 IOobject
61 (
62 "ITHACAdict",
63 runTime.system(),
64 mesh,
65 IOobject::MUST_READ,
66 IOobject::NO_WRITE
67 )
68 );
69 tolerance = ITHACAdict->lookupOrDefault<scalar>("tolerance", 1e-5);
70 maxIter = ITHACAdict->lookupOrDefault<scalar>("maxIter", 1000);
71 bcMethod = ITHACAdict->lookupOrDefault<word>("bcMethod", "lift");
72 timeDerivativeSchemeOrder =
73 ITHACAdict->lookupOrDefault<word>("timeDerivativeSchemeOrder", "second");
74 M_Assert
75 (
76 bcMethod == "lift" || bcMethod == "penalty",
77 "The BC method must be set to lift or penalty in ITHACAdict"
78 );
79 M_Assert
80 (
81 timeDerivativeSchemeOrder == "first" || timeDerivativeSchemeOrder == "second",
82 "The time derivative approximation must be set to either first or second order scheme in ITHACAdict"
83 );
84 para = ITHACAparameters::getInstance(mesh, runTime);
85 offline = ITHACAutilities::check_off();
86 podex = ITHACAutilities::check_pod();
87 supex = ITHACAutilities::check_sup();
88 NUmodesOut = para->ITHACAdict->lookupOrDefault<label>("NmodesUout", 15);
89 NPmodesOut = para->ITHACAdict->lookupOrDefault<label>("NmodesPout", 15);
90 NNutModesOut = para->ITHACAdict->lookupOrDefault<label>("NmodesNutOut", 15);
91 NUmodes = para->ITHACAdict->lookupOrDefault<label>("NmodesUproj", 10);
92 NSUPmodes = para->ITHACAdict->lookupOrDefault<label>("NmodesSUPproj", 10);
93 NPmodes = para->ITHACAdict->lookupOrDefault<label>("NmodesPproj", 10);
94 NNutModes = para->ITHACAdict->lookupOrDefault<label>("NmodesNutProj", 0);
95}
96
97// Small construct to access member functions linked to Smagorinsky diffusion
98UnsteadyNSTurb::UnsteadyNSTurb(const Parameters* myParameters):
99 m_parameters(static_cast<const StoredParameters*>(myParameters))
100{
101}
102
103// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
104
105void UnsteadyNSTurb::truthSolve(List<scalar> mu_now, std::string& offlinepath)
106{
107 Time& runTime = _runTime();
108 surfaceScalarField& phi = _phi();
109 fvMesh& mesh = _mesh();
110#include "initContinuityErrs.H"
111 fv::options& fvOptions = _fvOptions();
112 pimpleControl& pimple = _pimple();
113 volScalarField& p = _p();
114 volVectorField& U = _U();
115 volScalarField& nut = _nut();
116 IOMRFZoneList& MRF = _MRF();
117 singlePhaseTransportModel& laminarTransport = _laminarTransport();
118 instantList Times = runTime.times();
119 label& pRefCell = _pRefCell;
120 scalar& pRefValue = _pRefValue;
121 mesh.setFluxRequired(p.name());
122 runTime.setEndTime(finalTime);
123 runTime.setTime(Times[1], 1);
124 runTime.setDeltaT(timeStep);
125 nextWrite = startTime + writeEvery;
126 label nsnapshots = 0;
127
128 while (runTime.run())
129 {
130#include "readTimeControls.H"
131#include "CourantNo.H"
132#include "setDeltaT.H"
133 ++runTime;
134 Info << "Time = " << runTime.timeName() << nl << endl;
135
136 while (pimple.loop())
137 {
138#include "UEqn.H"
139
140 while (pimple.correct())
141 {
142#include "pEqn.H"
143 }
144
145 if (pimple.turbCorr())
146 {
147 laminarTransport.correct();
148 turbulence->correct();
149 }
150 }
151
152 Info << "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
153 << " ClockTime = " << runTime.elapsedClockTime() << " s"
154 << nl << endl;
155
156 if (checkWrite(runTime))
157 {
158 ++nsnapshots;
159 // Produces error when uncommented
160 // volScalarField nut = turbulence->nut().ref();
161 nut = turbulence->nut();
162 ITHACAstream::exportSolution(U, name(counter), offlinepath);
163 ITHACAstream::exportSolution(p, name(counter), offlinepath);
164 ITHACAstream::exportSolution(nut, name(counter), offlinepath);
165 std::ofstream of(offlinepath + name(counter) + "/" + runTime.timeName());
166 Ufield.append(tmp<volVectorField>(U));
167 Pfield.append(tmp<volScalarField>(p));
168 nutFields.append(tmp<volScalarField>(nut));
169 ++counter;
170 nextWrite += writeEvery;
171 writeMu(mu_now);
172 mu_samples.conservativeResize(mu_samples.rows() + 1, mu_now.size() + 1);
173 mu_samples(mu_samples.rows() - 1, 0) = atof(runTime.timeName().c_str());
174
175 for (label i = 0; i < mu_now.size(); ++i)
176 {
177 mu_samples(mu_samples.rows() - 1, i + 1) = mu_now[i];
178 }
179 }
180 }
181
182 if (mu.cols() == 0)
183 {
184 mu.resize(1, 1);
185 }
186
187 if (mu_samples.rows() == nsnapshots * mu.cols())
188 {
189 ITHACAstream::exportMatrix(mu_samples, "mu_samples", "eigen", offlinepath);
190 }
191}
192
193// ====== SUP Full Tensor 1 ======
194Eigen::Tensor<double, 3>
195UnsteadyNSTurb::turbulenceTensor1(label NUmodes, label NSUPmodes,
196 label nNutModes)
197{
198 const label cSize = NUmodes + NSUPmodes + liftfield.size();
199 Eigen::Tensor<double, 3> ct1Tensor(cSize, nNutModes, cSize);
200
201 for (label i = 0; i < cSize; ++i)
202 {
203 for (label j = 0; j < nNutModes; ++j)
204 {
205 for (label k = 0; k < cSize; ++k)
206 {
207 ct1Tensor(i, j, k) =
208 fvc::domainIntegrate
209 (
210 L_U_SUPmodes[i]
211 & fvc::laplacian(nutModes[j], L_U_SUPmodes[k])
212 ).value();
213 }
214 }
215 }
216
217 ITHACAstream::SaveDenseTensor
218 (
219 ct1Tensor,
220 "./ITHACAoutput/Matrices/",
221 "ct1_" + name(liftfield.size()) + "_" + name(NUmodes) + "_" + name(
222 NSUPmodes) + "_" + name(nNutModes) + "_t"
223 );
224 return ct1Tensor;
225}
226
227// ====== SUP Average Tensor 1 ======
228Eigen::Tensor<double, 3>
229UnsteadyNSTurb::turbulenceAveTensor1(label NUmodes, label NSUPmodes)
230{
231 const label cSize = NUmodes + NSUPmodes + liftfield.size();
232 const label nAvg = nutAve.size();
233 Eigen::Tensor<double, 3> ct1AveTensor(cSize, nAvg, cSize);
234
235 for (label i = 0; i < cSize; ++i)
236 {
237 for (label j = 0; j < nAvg; ++j)
238 {
239 for (label k = 0; k < cSize; ++k)
240 {
241 ct1AveTensor(i, j, k) =
242 fvc::domainIntegrate
243 (
244 L_U_SUPmodes[i]
245 & fvc::laplacian(nutAve[j], L_U_SUPmodes[k])
246 ).value();
247 }
248 }
249 }
250
251 ITHACAstream::SaveDenseTensor
252 (
253 ct1AveTensor,
254 "./ITHACAoutput/Matrices/",
255 "ct1Ave_" + name(liftfield.size()) + "_" + name(NUmodes) + "_" + name(
256 NSUPmodes) + "_t"
257 );
258 return ct1AveTensor;
259}
260
261// ====== SUP Fluctuation Tensor 1 ======
262Eigen::Tensor<double, 3>
263UnsteadyNSTurb::turbulenceFluctTensor1(label NUmodes, label NSUPmodes)
264{
265 const label cSize = NUmodes + NSUPmodes + liftfield.size();
266 const label nFluct = nutFluctModes.size();
267 Eigen::Tensor<double, 3> ct1FluctTensor(cSize, nFluct, cSize);
268
269 for (label i = 0; i < cSize; ++i)
270 {
271 for (label j = 0; j < nFluct; ++j)
272 {
273 for (label k = 0; k < cSize; ++k)
274 {
275 ct1FluctTensor(i, j, k) =
276 fvc::domainIntegrate
277 (
278 L_U_SUPmodes[i]
279 & fvc::laplacian(nutFluctModes[j], L_U_SUPmodes[k])
280 ).value();
281 }
282 }
283 }
284
285 ITHACAstream::SaveDenseTensor
286 (
287 ct1FluctTensor,
288 "./ITHACAoutput/Matrices/",
289 "ct1Fluct_" + name(liftfield.size()) + "_" + name(NUmodes) + "_" + name(
290 NSUPmodes) + "_t"
291 );
292 return ct1FluctTensor;
293}
294
295// ====== PPE Full Tensor 1 ======
296Eigen::Tensor<double, 3>
297UnsteadyNSTurb::turbulencePPETensor1(label NUmodes, label NSUPmodes,
298 label NPmodes, label nNutModes)
299{
300 const label cSize = NUmodes + NSUPmodes + liftfield.size();
301 Eigen::Tensor<double, 3> ct1PPETensor(NPmodes, nNutModes, cSize);
302
303 for (label i = 0; i < NPmodes; ++i)
304 {
305 for (label j = 0; j < nNutModes; ++j)
306 {
307 for (label k = 0; k < cSize; ++k)
308 {
309 ct1PPETensor(i, j, k) =
310 fvc::domainIntegrate
311 (
312 fvc::grad(Pmodes[i])
313 & fvc::laplacian(nutModes[j], L_U_SUPmodes[k])
314 ).value();
315 }
316 }
317 }
318
319 ITHACAstream::SaveDenseTensor
320 (
321 ct1PPETensor,
322 "./ITHACAoutput/Matrices/",
323 "ct1PPE_" + name(liftfield.size()) + "_" + name(NUmodes) + "_" + name(
324 NSUPmodes) + "_"
325 + name(NPmodes) + "_" + name(nNutModes) + "_t"
326 );
327 return ct1PPETensor;
328}
329
330// ====== PPE Average Tensor 1 ======
331Eigen::Tensor<double, 3>
332UnsteadyNSTurb::turbulencePPEAveTensor1(label NUmodes, label NSUPmodes,
333 label NPmodes)
334{
335 const label cSize = NUmodes + NSUPmodes + liftfield.size();
336 const label nAvg = nutAve.size();
337 Eigen::Tensor<double, 3> ct1PPEAveTensor(NPmodes, nAvg, cSize);
338
339 for (label i = 0; i < NPmodes; ++i)
340 {
341 for (label j = 0; j < nAvg; ++j)
342 {
343 for (label k = 0; k < cSize; ++k)
344 {
345 ct1PPEAveTensor(i, j, k) =
346 fvc::domainIntegrate
347 (
348 fvc::grad(Pmodes[i])
349 & fvc::laplacian(nutAve[j], L_U_SUPmodes[k])
350 ).value();
351 }
352 }
353 }
354
355 ITHACAstream::SaveDenseTensor
356 (
357 ct1PPEAveTensor,
358 "./ITHACAoutput/Matrices/",
359 "ct1PPEAve_" + name(liftfield.size()) + "_" + name(NUmodes) + "_" + name(
360 NSUPmodes) + "_"
361 + name(NPmodes) + "_t"
362 );
363 return ct1PPEAveTensor;
364}
365
366// ====== PPE Fluctuation Tensor 1 ======
367Eigen::Tensor<double, 3>
368UnsteadyNSTurb::turbulencePPEFluctTensor1(label NUmodes, label NSUPmodes,
369 label NPmodes)
370{
371 const label cSize = NUmodes + NSUPmodes + liftfield.size();
372 const label nFluct = nutFluctModes.size();
373 Eigen::Tensor<double, 3> ct1PPEFluctTensor(NPmodes, nFluct, cSize);
374
375 for (label i = 0; i < NPmodes; ++i)
376 {
377 for (label j = 0; j < nFluct; ++j)
378 {
379 for (label k = 0; k < cSize; ++k)
380 {
381 ct1PPEFluctTensor(i, j, k) =
382 fvc::domainIntegrate
383 (
384 fvc::grad(Pmodes[i])
385 & fvc::laplacian(nutFluctModes[j], L_U_SUPmodes[k])
386 ).value();
387 }
388 }
389 }
390
391 ITHACAstream::SaveDenseTensor
392 (
393 ct1PPEFluctTensor,
394 "./ITHACAoutput/Matrices/",
395 "ct1PPEFluct_" + name(liftfield.size()) + "_" + name(NUmodes) + "_" + name(
396 NSUPmodes) + "_"
397 + name(NPmodes) + "_t"
398 );
399 return ct1PPEFluctTensor;
400}
401
402// ====== SUP Full Tensor 2 ======
403Eigen::Tensor<double, 3>
404UnsteadyNSTurb::turbulenceTensor2(label NUmodes, label NSUPmodes,
405 label nNutModes)
406{
407 const label cSize = NUmodes + NSUPmodes + liftfield.size();
408 Eigen::Tensor<double, 3> ct2Tensor(cSize, nNutModes, cSize);
409
410 for (label i = 0; i < cSize; ++i)
411 {
412 for (label j = 0; j < nNutModes; ++j)
413 {
414 for (label k = 0; k < cSize; ++k)
415 {
416 ct2Tensor(i, j, k) =
417 fvc::domainIntegrate
418 (
419 L_U_SUPmodes[i]
420 & fvc::div(nutModes[j] * dev2((fvc::grad(L_U_SUPmodes[k]))().T()))
421 ).value();
422 }
423 }
424 }
425
426 ITHACAstream::SaveDenseTensor
427 (
428 ct2Tensor,
429 "./ITHACAoutput/Matrices/",
430 "ct2_" + name(liftfield.size()) + "_" + name(NUmodes) + "_" + name(
431 NSUPmodes) + "_" + name(nNutModes) + "_t"
432 );
433 return ct2Tensor;
434}
435
436// ====== SUP Average Tensor 2 ======
437Eigen::Tensor<double, 3>
438UnsteadyNSTurb::turbulenceAveTensor2(label NUmodes, label NSUPmodes)
439{
440 const label cSize = NUmodes + NSUPmodes + liftfield.size();
441 const label nAvg = nutAve.size();
442 Eigen::Tensor<double, 3> ct2AveTensor(cSize, nAvg, cSize);
443
444 for (label i = 0; i < cSize; ++i)
445 {
446 for (label j = 0; j < nAvg; ++j)
447 {
448 for (label k = 0; k < cSize; ++k)
449 {
450 ct2AveTensor(i, j, k) =
451 fvc::domainIntegrate
452 (
453 L_U_SUPmodes[i]
454 & fvc::div(nutAve[j] * dev2((fvc::grad(L_U_SUPmodes[k]))().T()))
455 ).value();
456 }
457 }
458 }
459
460 ITHACAstream::SaveDenseTensor
461 (
462 ct2AveTensor,
463 "./ITHACAoutput/Matrices/",
464 "ct2Ave_" + name(liftfield.size()) + "_" + name(NUmodes) + "_" + name(
465 NSUPmodes) + "_t"
466 );
467 return ct2AveTensor;
468}
469
470// ====== SUP Fluctuation Tensor 2 ======
471Eigen::Tensor<double, 3>
472UnsteadyNSTurb::turbulenceFluctTensor2(label NUmodes, label NSUPmodes)
473{
474 const label cSize = NUmodes + NSUPmodes + liftfield.size();
475 const label nFluct = nutFluctModes.size();
476 Eigen::Tensor<double, 3> ct2FluctTensor(cSize, nFluct, cSize);
477
478 for (label i = 0; i < cSize; ++i)
479 {
480 for (label j = 0; j < nFluct; ++j)
481 {
482 for (label k = 0; k < cSize; ++k)
483 {
484 ct2FluctTensor(i, j, k) =
485 fvc::domainIntegrate
486 (
487 L_U_SUPmodes[i]
488 & fvc::div(nutFluctModes[j] * dev2((fvc::grad(L_U_SUPmodes[k]))().T()))
489 ).value();
490 }
491 }
492 }
493
494 ITHACAstream::SaveDenseTensor
495 (
496 ct2FluctTensor,
497 "./ITHACAoutput/Matrices/",
498 "ct2Fluct_" + name(liftfield.size()) + "_" + name(NUmodes) + "_" + name(
499 NSUPmodes) + "_t"
500 );
501 return ct2FluctTensor;
502}
503
504// ====== PPE Full Tensor 2 ======
505Eigen::Tensor<double, 3>
506UnsteadyNSTurb::turbulencePPETensor2(label NUmodes, label NSUPmodes,
507 label NPmodes, label nNutModes)
508{
509 const label cSize = NUmodes + NSUPmodes + liftfield.size();
510 Eigen::Tensor<double, 3> ct2PPETensor(NPmodes, nNutModes, cSize);
511
512 for (label i = 0; i < NPmodes; ++i)
513 {
514 for (label j = 0; j < nNutModes; ++j)
515 {
516 for (label k = 0; k < cSize; ++k)
517 {
518 ct2PPETensor(i, j, k) =
519 fvc::domainIntegrate
520 (
521 fvc::grad(Pmodes[i])
522 & fvc::div(nutModes[j] * dev2((fvc::grad(L_U_SUPmodes[k]))().T()))
523 ).value();
524 }
525 }
526 }
527
528 ITHACAstream::SaveDenseTensor
529 (
530 ct2PPETensor,
531 "./ITHACAoutput/Matrices/",
532 "ct2PPE_" + name(liftfield.size()) + "_" + name(NUmodes) + "_" + name(
533 NSUPmodes) + "_"
534 + name(NPmodes) + "_" + name(nNutModes) + "_t"
535 );
536 return ct2PPETensor;
537}
538
539// ====== PPE Average Tensor 2 ======
540Eigen::Tensor<double, 3>
541UnsteadyNSTurb::turbulencePPEAveTensor2(label NUmodes, label NSUPmodes,
542 label NPmodes)
543{
544 const label cSize = NUmodes + NSUPmodes + liftfield.size();
545 const label nAvg = nutAve.size();
546 Eigen::Tensor<double, 3> ct2PPEAveTensor(NPmodes, nAvg, cSize);
547
548 for (label i = 0; i < NPmodes; ++i)
549 {
550 for (label j = 0; j < nAvg; ++j)
551 {
552 for (label k = 0; k < cSize; ++k)
553 {
554 ct2PPEAveTensor(i, j, k) =
555 fvc::domainIntegrate
556 (
557 fvc::grad(Pmodes[i])
558 & fvc::div(nutAve[j] * dev2((fvc::grad(L_U_SUPmodes[k]))().T()))
559 ).value();
560 }
561 }
562 }
563
564 ITHACAstream::SaveDenseTensor
565 (
566 ct2PPEAveTensor,
567 "./ITHACAoutput/Matrices/",
568 "ct2PPEAve_" + name(liftfield.size()) + "_" + name(NUmodes) + "_" + name(
569 NSUPmodes) + "_"
570 + name(NPmodes) + "_t"
571 );
572 return ct2PPEAveTensor;
573}
574
575// ====== PPE Fluctuation Tensor 2 ======
576Eigen::Tensor<double, 3>
577UnsteadyNSTurb::turbulencePPEFluctTensor2(label NUmodes, label NSUPmodes,
578 label NPmodes)
579{
580 const label cSize = NUmodes + NSUPmodes + liftfield.size();
581 const label nFluct = nutFluctModes.size();
582 Eigen::Tensor<double, 3> ct2PPEFluctTensor(NPmodes, nFluct, cSize);
583
584 for (label i = 0; i < NPmodes; ++i)
585 {
586 for (label j = 0; j < nFluct; ++j)
587 {
588 for (label k = 0; k < cSize; ++k)
589 {
590 ct2PPEFluctTensor(i, j, k) =
591 fvc::domainIntegrate
592 (
593 fvc::grad(Pmodes[i])
594 & fvc::div(nutFluctModes[j] * dev2((fvc::grad(L_U_SUPmodes[k]))().T()))
595 ).value();
596 }
597 }
598 }
599
600 ITHACAstream::SaveDenseTensor
601 (
602 ct2PPEFluctTensor,
603 "./ITHACAoutput/Matrices/",
604 "ct2PPEFluct_" + name(liftfield.size()) + "_" + name(NUmodes) + "_" + name(
605 NSUPmodes) + "_"
606 + name(NPmodes) + "_t"
607 );
608 return ct2PPEFluctTensor;
609}
610
611// ====== BT Turbulence Matrix (SUP) ======
612Eigen::MatrixXd UnsteadyNSTurb::btTurbulence(label NUmodes, label NSUPmodes)
613{
614 const label btSize = NUmodes + NSUPmodes + liftfield.size();
615 Eigen::MatrixXd btMatrix(btSize, btSize);
616 btMatrix = btMatrix * 0;
617
618 for (label i = 0; i < btSize; ++i)
619 {
620 for (label j = 0; j < btSize; ++j)
621 {
622 btMatrix(i, j) =
623 fvc::domainIntegrate
624 (
625 L_U_SUPmodes[i]
626 & fvc::div(dev2((T(fvc::grad(L_U_SUPmodes[j])))))
627 ).value();
628 }
629 }
630
631 ITHACAstream::SaveDenseMatrix
632 (
633 btMatrix,
634 "./ITHACAoutput/Matrices/",
635 "bt_" + name(liftfield.size()) + "_" + name(NUmodes) + "_" + name(NSUPmodes)
636 );
637 return btMatrix;
638}
639
640Eigen::MatrixXd UnsteadyNSTurb::continuity_matrix(label NUmodes,
641 label NSUPmodes, label NPmodes)
642{
643 const label cSize = NUmodes + NSUPmodes + liftfield.size();
644 M_Assert(cSize > 0, "continuity_matrix: cSize=0");
645 M_Assert(NPmodes > 0, "continuity_matrix: NPmodes=0");
646 Eigen::MatrixXd Pmat(NPmodes, cSize);
647 Pmat.setZero();
648
649 for (label i = 0; i < NPmodes; ++i)
650 {
651 for (label j = 0; j < cSize; ++j)
652 {
653 Pmat(i, j) =
654 fvc::domainIntegrate(Pmodes[i] * fvc::div(L_U_SUPmodes[j])).value();
655 }
656 }
657
658 ITHACAstream::SaveDenseMatrix
659 (
660 Pmat,
661 "./ITHACAoutput/Matrices/",
662 "P_" + name(liftfield.size()) + "_" + name(NUmodes) + "_" + name(
663 NSUPmodes) + "_" + name(NPmodes)
664 );
665 return Pmat;
666}
667
668Eigen::MatrixXd UnsteadyNSTurb::pressurePPE_L(label NPmodes)
669{
670 Eigen::MatrixXd Lvec(NPmodes, 1);
671 Lvec.setZero();
672 const word muPath = "./ITHACAoutput/Offline/mu_samples_mat.txt";
673 M_Assert
674 (
675 ITHACAutilities::check_file(muPath),
676 "[PPE-L] mu_samples_mat.txt not found; cannot compute R_t."
677 );
678 Eigen::MatrixXd muMat = ITHACAstream::readMatrix(muPath);
679 M_Assert
680 (
681 muMat.rows() >= 2,
682 "[PPE-L] mu/time log has <2 rows; need two snapshots for R_t."
683 );
684 const double dt0 = muMat(1, 0) - muMat(0, 0);
685 M_Assert(std::abs(dt0) > SMALL, "[PPE-L] dt0 ~ 0; check mu_samples_mat.txt.");
686 M_Assert(Uomfield.size() >= 2,
687 "[PPE-L] Need at least 2 velocity snapshots to build R_t.");
688 const volVectorField& U0 = Uomfield[0];
689 const volVectorField& U1 = Uomfield[1];
690 tmp<volVectorField> tUdot = (U1 - U0) / dt0;
691 const volVectorField& Udot = tUdot();
692 const fvMesh& mesh = Pmodes[0]().mesh();
693
694 for (label i = 0; i < NPmodes; ++i)
695 {
696 scalar Li = 0.0;
697 forAll(mesh.boundary(), patchi)
698 {
699 const scalarField& chiF = Pmodes[i].boundaryField()[patchi];
700 const vectorField& UdotF = Udot.boundaryField()[patchi];
701 const vectorField& Sf = mesh.Sf().boundaryField()[patchi];
702 Li += gSum(chiF * (UdotF & Sf));
703 }
704 Lvec(i, 0) = Li;
705 }
706
707 ITHACAstream::SaveDenseMatrix(Lvec, "./ITHACAoutput/Matrices/",
708 "L_" + name(NPmodes));
709 return Lvec;
710}
711
712void UnsteadyNSTurb::projectSUP(fileName folder,
713 label NU,
714 label NP,
715 label NSUP,
716 label Nnut,
717 bool rbfInterp)
718{
719 Info << "[DEBUG] Entered projectSUP()" << endl;
720 NUmodes = NU;
721 NPmodes = NP;
722 NSUPmodes = NSUP;
723 nNutModes = Nnut;
724 L_U_SUPmodes.resize(0);
725
726 if (liftfield.size() != 0)
727 {
728 for (label k = 0; k < liftfield.size(); ++k)
729 {
730 L_U_SUPmodes.append(tmp<volVectorField>(liftfield[k]));
731 }
732 }
733
734 if (NUmodes != 0)
735 {
736 for (label k = 0; k < NUmodes; ++k)
737 {
738 L_U_SUPmodes.append(tmp<volVectorField>(Umodes[k]));
739 }
740 }
741
742 if (NSUPmodes != 0)
743 {
744 for (label k = 0; k < NSUPmodes; ++k)
745 {
746 L_U_SUPmodes.append(tmp<volVectorField>(supmodes[k]));
747 }
748 }
749
750 const label cSize = liftfield.size() + NUmodes + NSUPmodes;
751
752 if (ITHACAutilities::check_folder("./ITHACAoutput/Matrices/"))
753 {
754 {
755 const word M_str =
756 "M_" + name(liftfield.size()) + "_" + name(NUmodes) + "_" + name(NSUPmodes);
757
758 if (ITHACAutilities::check_file("./ITHACAoutput/Matrices/" + M_str))
759 {
760 ITHACAstream::ReadDenseMatrix(M_matrix, "./ITHACAoutput/Matrices/", M_str);
761 }
762 else
763 {
764 M_matrix = mass_term(NUmodes, NPmodes, NSUPmodes);
765 }
766 }
767 {
768 const word B_str =
769 "B_" + name(liftfield.size()) + "_" + name(NUmodes) + "_" + name(NSUPmodes);
770
771 if (ITHACAutilities::check_file("./ITHACAoutput/Matrices/" + B_str))
772 {
773 ITHACAstream::ReadDenseMatrix(B_matrix, "./ITHACAoutput/Matrices/", B_str);
774 }
775 else
776 {
777 B_matrix = diffusive_term(NUmodes, NPmodes, NSUPmodes);
778 }
779 }
780 {
781 const word btStr =
782 "bt_" + name(liftfield.size()) + "_" + name(NUmodes) + "_" + name(NSUPmodes);
783
784 if (ITHACAutilities::check_file("./ITHACAoutput/Matrices/" + btStr))
785 {
786 ITHACAstream::ReadDenseMatrix(btMatrix, "./ITHACAoutput/Matrices/", btStr);
787 }
788 else
789 {
790 btMatrix = btTurbulence(NUmodes, NSUPmodes);
791 }
792 }
793 {
794 const word K_str =
795 "K_" + name(liftfield.size()) + "_" + name(NUmodes) + "_" + name(
796 NSUPmodes) + "_" + name(NPmodes);
797
798 if (ITHACAutilities::check_file("./ITHACAoutput/Matrices/" + K_str))
799 {
800 ITHACAstream::ReadDenseMatrix(K_matrix, "./ITHACAoutput/Matrices/", K_str);
801 }
802 else
803 {
804 K_matrix = pressure_gradient_term(NUmodes, NPmodes, NSUPmodes);
805 }
806 }
807 {
808 const word P_str =
809 "P_" + name(liftfield.size()) + "_" + name(NUmodes) + "_" + name(
810 NSUPmodes) + "_" + name(NPmodes);
811
812 if (ITHACAutilities::check_file("./ITHACAoutput/Matrices/" + P_str))
813 {
814 ITHACAstream::ReadDenseMatrix(P_matrix, "./ITHACAoutput/Matrices/", P_str);
815 }
816 else
817 {
818 P_matrix = continuity_matrix(NUmodes, NSUPmodes, NPmodes);
819 }
820 }
821 {
822 const word C_str =
823 "C_" + name(liftfield.size()) + "_" + name(NUmodes) + "_" + name(
824 NSUPmodes) + "_t";
825
826 if (ITHACAutilities::check_file("./ITHACAoutput/Matrices/" + C_str))
827 {
828 ITHACAstream::ReadDenseTensor(C_tensor, "./ITHACAoutput/Matrices/", C_str);
829 }
830 else
831 {
832 C_tensor = convective_term_tens(NUmodes, NPmodes, NSUPmodes);
833 }
834 }
835 {
836 const word ct1Str =
837 "ct1_" + name(liftfield.size()) + "_" + name(NUmodes) + "_" + name(
838 NSUPmodes) + "_"
839 + name(nNutModes) + "_t";
840
841 if (ITHACAutilities::check_file("./ITHACAoutput/Matrices/" + ct1Str))
842 {
843 ITHACAstream::ReadDenseTensor(ct1Tensor, "./ITHACAoutput/Matrices/", ct1Str);
844 }
845 else
846 {
847 ct1Tensor = turbulenceTensor1(NUmodes, NSUPmodes, nNutModes);
848 }
849
850 const word ct2Str =
851 "ct2_" + name(liftfield.size()) + "_" + name(NUmodes) + "_" + name(
852 NSUPmodes) + "_"
853 + name(nNutModes) + "_t";
854
855 if (ITHACAutilities::check_file("./ITHACAoutput/Matrices/" + ct2Str))
856 {
857 ITHACAstream::ReadDenseTensor(ct2Tensor, "./ITHACAoutput/Matrices/", ct2Str);
858 }
859 else
860 {
861 ct2Tensor = turbulenceTensor2(NUmodes, NSUPmodes, nNutModes);
862 }
863 }
864
865 if (nutAve.size() != 0)
866 {
867 const word ct1AveStr =
868 "ct1Ave_" + name(liftfield.size()) + "_" + name(NUmodes) + "_" + name(
869 NSUPmodes) + "_t";
870
871 if (ITHACAutilities::check_file("./ITHACAoutput/Matrices/" + ct1AveStr))
872 {
873 ITHACAstream::ReadDenseTensor(ct1AveTensor, "./ITHACAoutput/Matrices/",
874 ct1AveStr);
875 }
876 else
877 {
878 ct1AveTensor = turbulenceAveTensor1(NUmodes, NSUPmodes);
879 }
880
881 const word ct2AveStr =
882 "ct2Ave_" + name(liftfield.size()) + "_" + name(NUmodes) + "_" + name(
883 NSUPmodes) + "_t";
884
885 if (ITHACAutilities::check_file("./ITHACAoutput/Matrices/" + ct2AveStr))
886 {
887 ITHACAstream::ReadDenseTensor(ct2AveTensor, "./ITHACAoutput/Matrices/",
888 ct2AveStr);
889 }
890 else
891 {
892 ct2AveTensor = turbulenceAveTensor2(NUmodes, NSUPmodes);
893 }
894 }
895
896 if (nutFluctModes.size() != 0)
897 {
898 const word ct1FluctStr =
899 "ct1Fluct_" + name(liftfield.size()) + "_" + name(NUmodes) + "_" + name(
900 NSUPmodes) + "_t";
901
902 if (ITHACAutilities::check_file("./ITHACAoutput/Matrices/" + ct1FluctStr))
903 {
904 ITHACAstream::ReadDenseTensor(ct1FluctTensor, "./ITHACAoutput/Matrices/",
905 ct1FluctStr);
906 }
907 else
908 {
909 ct1FluctTensor = turbulenceFluctTensor1(NUmodes, NSUPmodes);
910 }
911
912 const word ct2FluctStr =
913 "ct2Fluct_" + name(liftfield.size()) + "_" + name(NUmodes) + "_" + name(
914 NSUPmodes) + "_t";
915
916 if (ITHACAutilities::check_file("./ITHACAoutput/Matrices/" + ct2FluctStr))
917 {
918 ITHACAstream::ReadDenseTensor(ct2FluctTensor, "./ITHACAoutput/Matrices/",
919 ct2FluctStr);
920 }
921 else
922 {
923 ct2FluctTensor = turbulenceFluctTensor2(NUmodes, NSUPmodes);
924 }
925 }
926
927 if (bcMethod == "penalty")
928 {
929 bcVelVec = bcVelocityVec(NUmodes, NSUPmodes);
930 bcVelMat = bcVelocityMat(NUmodes, NSUPmodes);
931 }
932 }
933 else
934 {
935 M_matrix = mass_term(NUmodes, NPmodes, NSUPmodes);
936 B_matrix = diffusive_term(NUmodes, NPmodes, NSUPmodes);
937 btMatrix = btTurbulence(NUmodes, NSUPmodes);
938 K_matrix = pressure_gradient_term(NUmodes, NPmodes, NSUPmodes);
939 P_matrix = continuity_matrix(NUmodes, NSUPmodes, NPmodes);
940 C_tensor = convective_term_tens(NUmodes, NPmodes, NSUPmodes);
941 ct1Tensor = turbulenceTensor1(NUmodes, NSUPmodes, nNutModes);
942 ct2Tensor = turbulenceTensor2(NUmodes, NSUPmodes, nNutModes);
943
944 if (nutAve.size() != 0)
945 {
946 ct1AveTensor = turbulenceAveTensor1(NUmodes, NSUPmodes);
947 ct2AveTensor = turbulenceAveTensor2(NUmodes, NSUPmodes);
948 }
949
950 if (nutFluctModes.size() != 0)
951 {
952 ct1FluctTensor = turbulenceFluctTensor1(NUmodes, NSUPmodes);
953 ct2FluctTensor = turbulenceFluctTensor2(NUmodes, NSUPmodes);
954 }
955
956 if (bcMethod == "penalty")
957 {
958 bcVelVec = bcVelocityVec(NUmodes, NSUPmodes);
959 bcVelMat = bcVelocityMat(NUmodes, NSUPmodes);
960 }
961 }
962
963 bTotalMatrix = B_matrix + btMatrix;
964 cTotalTensor.resize(cSize, nNutModes, cSize);
965 cTotalTensor = ct1Tensor + ct2Tensor;
966
967 if (nutAve.size() != 0)
968 {
969 cTotalAveTensor.resize(cSize, nutAve.size(), cSize);
970 cTotalAveTensor = ct1AveTensor + ct2AveTensor;
971 }
972
973 if (nutFluctModes.size() != 0)
974 {
975 cTotalFluctTensor.resize(cSize, nutFluctModes.size(), cSize);
976 cTotalFluctTensor = ct1FluctTensor + ct2FluctTensor;
977 }
978
979 if (para->exportPython)
980 {
981 ITHACAstream::exportMatrix(M_matrix, "M", "python",
982 "./ITHACAoutput/Matrices/");
983 ITHACAstream::exportMatrix(B_matrix, "B", "python",
984 "./ITHACAoutput/Matrices/");
985 ITHACAstream::exportMatrix(btMatrix, "bt", "python",
986 "./ITHACAoutput/Matrices/");
987 ITHACAstream::exportMatrix(bTotalMatrix, "bTot", "python",
988 "./ITHACAoutput/Matrices/");
989 ITHACAstream::exportMatrix(K_matrix, "K", "python",
990 "./ITHACAoutput/Matrices/");
991 ITHACAstream::exportMatrix(P_matrix, "P", "python",
992 "./ITHACAoutput/Matrices/");
993 ITHACAstream::exportTensor(C_tensor, "C", "python",
994 "./ITHACAoutput/Matrices/");
995 ITHACAstream::exportTensor(ct1Tensor, "ct1", "python",
996 "./ITHACAoutput/Matrices/");
997 ITHACAstream::exportTensor(ct2Tensor, "ct2", "python",
998 "./ITHACAoutput/Matrices/");
999
1000 if (nutAve.size() != 0)
1001 {
1002 ITHACAstream::exportTensor(ct1AveTensor, "ct1Ave", "python",
1003 "./ITHACAoutput/Matrices/");
1004 ITHACAstream::exportTensor(ct2AveTensor, "ct2Ave", "python",
1005 "./ITHACAoutput/Matrices/");
1006 }
1007
1008 if (nutFluctModes.size() != 0)
1009 {
1010 ITHACAstream::exportTensor(ct1FluctTensor, "ct1Fluct", "python",
1011 "./ITHACAoutput/Matrices/");
1012 ITHACAstream::exportTensor(ct2FluctTensor, "ct2Fluct", "python",
1013 "./ITHACAoutput/Matrices/");
1014 }
1015 }
1016
1017 if (para->exportMatlab)
1018 {
1019 ITHACAstream::exportMatrix(M_matrix, "M", "matlab",
1020 "./ITHACAoutput/Matrices/");
1021 ITHACAstream::exportMatrix(B_matrix, "B", "matlab",
1022 "./ITHACAoutput/Matrices/");
1023 ITHACAstream::exportMatrix(btMatrix, "bt", "matlab",
1024 "./ITHACAoutput/Matrices/");
1025 ITHACAstream::exportMatrix(bTotalMatrix, "bTot", "matlab",
1026 "./ITHACAoutput/Matrices/");
1027 ITHACAstream::exportMatrix(K_matrix, "K", "matlab",
1028 "./ITHACAoutput/Matrices/");
1029 ITHACAstream::exportMatrix(P_matrix, "P", "matlab",
1030 "./ITHACAoutput/Matrices/");
1031 ITHACAstream::exportTensor(C_tensor, "C", "matlab",
1032 "./ITHACAoutput/Matrices/");
1033 ITHACAstream::exportTensor(ct1Tensor, "ct1", "matlab",
1034 "./ITHACAoutput/Matrices/");
1035 ITHACAstream::exportTensor(ct2Tensor, "ct2", "matlab",
1036 "./ITHACAoutput/Matrices/");
1037
1038 if (nutAve.size() != 0)
1039 {
1040 ITHACAstream::exportTensor(ct1AveTensor, "ct1Ave", "matlab",
1041 "./ITHACAoutput/Matrices/");
1042 ITHACAstream::exportTensor(ct2AveTensor, "ct2Ave", "matlab",
1043 "./ITHACAoutput/Matrices/");
1044 }
1045
1046 if (nutFluctModes.size() != 0)
1047 {
1048 ITHACAstream::exportTensor(ct1FluctTensor, "ct1Fluct", "matlab",
1049 "./ITHACAoutput/Matrices/");
1050 ITHACAstream::exportTensor(ct2FluctTensor, "ct2Fluct", "matlab",
1051 "./ITHACAoutput/Matrices/");
1052 }
1053 }
1054
1055 if (para->exportTxt)
1056 {
1057 ITHACAstream::exportMatrix(M_matrix, "M", "eigen",
1058 "./ITHACAoutput/Matrices/");
1059 ITHACAstream::exportMatrix(B_matrix, "B", "eigen",
1060 "./ITHACAoutput/Matrices/");
1061 ITHACAstream::exportMatrix(btMatrix, "bt", "eigen",
1062 "./ITHACAoutput/Matrices/");
1063 ITHACAstream::exportMatrix(bTotalMatrix, "bTot", "eigen",
1064 "./ITHACAoutput/Matrices/");
1065 ITHACAstream::exportMatrix(K_matrix, "K", "eigen",
1066 "./ITHACAoutput/Matrices/");
1067 ITHACAstream::exportMatrix(P_matrix, "P", "eigen",
1068 "./ITHACAoutput/Matrices/");
1069 ITHACAstream::exportTensor(C_tensor, "C", "eigen",
1070 "./ITHACAoutput/Matrices/C");
1071 ITHACAstream::exportTensor(ct1Tensor, "ct1_", "eigen",
1072 "./ITHACAoutput/Matrices/ct1");
1073 ITHACAstream::exportTensor(ct2Tensor, "ct2_", "eigen",
1074 "./ITHACAoutput/Matrices/ct2");
1075
1076 if (nutAve.size() != 0)
1077 {
1078 ITHACAstream::exportTensor(ct1AveTensor, "ct1Ave_", "eigen",
1079 "./ITHACAoutput/Matrices/ct1Ave");
1080 ITHACAstream::exportTensor(ct2AveTensor, "ct2Ave_", "eigen",
1081 "./ITHACAoutput/Matrices/ct2Ave");
1082 }
1083
1084 if (nutFluctModes.size() != 0)
1085 {
1086 ITHACAstream::exportTensor(ct1FluctTensor, "ct1Fluct", "eigen",
1087 "./ITHACAoutput/Matrices/");
1088 ITHACAstream::exportTensor(ct2FluctTensor, "ct2Fluct", "eigen",
1089 "./ITHACAoutput/Matrices/");
1090 }
1091 }
1092
1093 Info << "[DEBUG] SUP sizes: cSize=" << cSize
1094 << ", NP=" << NPmodes << ", nNut=" << nNutModes << '\n';
1095 Info << "[DEBUG] Shapes: M(" << M_matrix.rows() << "x" << M_matrix.cols()
1096 << ") B(" << B_matrix.rows() << "x" << B_matrix.cols()
1097 << ") bt(" << btMatrix.rows() << "x" << btMatrix.cols()
1098 << ") K(" << K_matrix.rows() << "x" << K_matrix.cols()
1099 << ") P(" << P_matrix.rows() << "x" << P_matrix.cols() << ")\n";
1100 Info << "[DEBUG] projectSUP() done. SUP matrices/tensors exported.\n";
1101
1102 if (rbfInterp && (!Pstream::parRun()))
1103 {
1104 Info
1105 << "\n==== [RBF DEBUG SUP] ENTERING OFFLINE RBF CONSTRUCTION (AVG linear-μ + FLUCT RBF) ====\n";
1106 const word coeffDir = "./ITHACAoutput/Coefficients/";
1107 const word muPath = "./ITHACAoutput/Offline/mu_samples_mat.txt";
1108 Eigen::MatrixXd muMat = ITHACAstream::readMatrix(muPath);
1109 Info << "[RBF DEBUG SUP] muMat shape: " << muMat.rows() << " x " <<
1110 muMat.cols() << endl;
1111 const int nPar = 12;
1112 const int nSnapshotsPerPar = 200;
1113 Eigen::VectorXd timeVec = muMat.col(0);
1114 Eigen::VectorXd muVec = muMat.col(1);
1115 Eigen::MatrixXd coeffNutAvg = ITHACAstream::readMatrix(
1116 coeffDir + "Nut_avg_coeffs_mat.txt");
1117 Eigen::MatrixXd coeffNutFluct = ITHACAstream::readMatrix(
1118 coeffDir + "Nut_fluct_coeffs_mat.txt");
1119 Info << "[RBF DEBUG SUP] Loaded coeffNutAvg shape: " << coeffNutAvg.rows()
1120 << " x " << coeffNutAvg.cols() << endl;
1121 Info << "[RBF DEBUG SUP] Loaded coeffNutFluct shape: " <<
1122 coeffNutFluct.rows() << " x " << coeffNutFluct.cols() << endl;
1123 const int nNutAvgModes = coeffNutAvg.rows();
1124 const int nNutFluctModes = coeffNutFluct.rows();
1125 const int nUniqueMu = nPar;
1126 Eigen::MatrixXd a = ITHACAutilities::getCoeffs(Uomfield, Umodes);
1127 a.transposeInPlace();
1128 Eigen::VectorXd initSnapInd(nPar);
1129 Eigen::VectorXd timeSnap(nPar);
1130
1131 for (int i = 0; i < nPar; ++i)
1132 {
1133 const int start = i * nSnapshotsPerPar;
1134 initSnapInd(i) = start;
1135 timeSnap(i) = timeVec(start + 1) - timeVec(start);
1136 Info << "[RBF DEBUG SUP] i=" << i << ", start=" << start << ", dt=" <<
1137 timeSnap(i) << endl;
1138 }
1139
1140 Eigen::VectorXd muVecUnique(nUniqueMu);
1141
1142 for (int i = 0; i < nUniqueMu; ++i)
1143 {
1144 muVecUnique(i) = muVec(static_cast<int>(initSnapInd(i)));
1145 }
1146
1147 Info << "[RBF DEBUG SUP] muVecUnique: [" << muVecUnique(0) << " ... "
1148 << muVecUnique(muVecUnique.size() - 1) << "] (M=" << muVecUnique.size() <<
1149 ")\n";
1150 Info <<
1151 "[RBF DEBUG SUP] Calling velDerivativeCoeff() for fluctuation part...\n";
1152 Eigen::MatrixXd Gfluct = coeffNutFluct.transpose();
1153 List<Eigen::MatrixXd> interpDataFluct = velDerivativeCoeff(a, Gfluct,
1154 initSnapInd, timeSnap);
1155 Eigen::MatrixXd velRBF_fluct = interpDataFluct[0];
1156 Eigen::MatrixXd coeffs_fluct = interpDataFluct[1];
1157 Info << "[RBF DEBUG SUP] velRBF_fluct shape: " << velRBF_fluct.rows() <<
1158 " x " << velRBF_fluct.cols() << endl;
1159 Info << "[RBF DEBUG SUP] coeffs_fluct shape: " << coeffs_fluct.rows() <<
1160 " x " << coeffs_fluct.cols() << endl;
1161 const int nSnapshots = velRBF_fluct.rows();
1162 const int nModes = velRBF_fluct.cols() / 2;
1163 Eigen::MatrixXd adot_only(nSnapshots, nModes);
1164
1165 for (int i = 0; i < nSnapshots; ++i)
1166 {
1167 adot_only.row(i) = velRBF_fluct.row(i).tail(nModes);
1168 }
1169
1170 ITHACAstream::exportMatrix(adot_only, "adot_coeffs", "python", coeffDir);
1171 ITHACAstream::exportMatrix(velRBF_fluct, "a_adot_concat", "python", coeffDir);
1172 ITHACAstream::exportMatrix(velRBF_fluct, "a_adot_concat", "eigen", coeffDir);
1173 const double eAvg = 3;
1174 const double eFluct = 0.05;
1175 Eigen::MatrixXd radiiAvg;
1176 Eigen::MatrixXd radiiFluct;
1177
1178 if (ITHACAutilities::check_file("./radii_avg.txt"))
1179 {
1180 radiiAvg = ITHACAstream::readMatrix("./radii_avg.txt");
1181 M_Assert(radiiAvg.size() == nNutAvgModes, "radiiAvg size mismatch");
1182 }
1183 else
1184 {
1185 radiiAvg = Eigen::MatrixXd::Ones(nNutAvgModes, 1) * eAvg;
1186 }
1187
1188 if (ITHACAutilities::check_file("./radii_fluct.txt"))
1189 {
1190 radiiFluct = ITHACAstream::readMatrix("./radii_fluct.txt");
1191 M_Assert(radiiFluct.size() == nNutFluctModes, "radiiFluct size mismatch");
1192 }
1193 else
1194 {
1195 radiiFluct = Eigen::MatrixXd::Ones(nNutFluctModes, 1) * eFluct;
1196 }
1197
1198 List<SPLINTER::DataTable*> samplesNutAvg;
1199 List<SPLINTER::RBFSpline*> rbfSplinesNutAvg;
1200 samplesNutAvg.resize(0);
1201 rbfSplinesNutAvg.resize(0);
1202 ITHACAstream::exportMatrix(muVecUnique, "NutAvg_mu_unique", "eigen",
1203 coeffDir);
1204 ITHACAstream::exportMatrix(coeffNutAvg, "NutAvg_coeffs_by_mu", "eigen",
1205 coeffDir);
1206 List<SPLINTER::DataTable*> samplesNutFluct;
1207 List<SPLINTER::RBFSpline*> rbfSplinesNutFluct;
1208 samplesNutFluct.resize(nNutFluctModes);
1209 rbfSplinesNutFluct.resize(nNutFluctModes);
1210 Info << ">>> [SUP] Building nut_fluct RBF splines...\n";
1211
1212 for (label i = 0; i < nNutFluctModes; ++i)
1213 {
1214 const word weightName = "wRBF_NUTFLUCT_" + name(i + 1);
1215 samplesNutFluct[i] = new SPLINTER::DataTable(velRBF_fluct.cols(), 1);
1216
1217 for (label j = 0; j < velRBF_fluct.rows(); ++j)
1218 {
1219 samplesNutFluct[i]->addSample(velRBF_fluct.row(j), coeffs_fluct(j, i));
1220 }
1221
1222 Eigen::MatrixXd weights;
1223
1224 if (ITHACAutilities::check_file("./ITHACAoutput/weightsSUP/" + weightName))
1225 {
1226 ITHACAstream::ReadDenseMatrix(weights, "./ITHACAoutput/weightsSUP/",
1227 weightName);
1228 rbfSplinesNutFluct[i] =
1229 new SPLINTER::RBFSpline
1230 (
1231 *samplesNutFluct[i],
1232 SPLINTER::RadialBasisFunctionType::GAUSSIAN,
1233 weights,
1234 radiiFluct(i)
1235 );
1236 Info << " [SUP] nut_fluct RBF " << i + 1 << "/" << nNutFluctModes
1237 << " loaded from ./ITHACAoutput/weightsSUP/" << weightName << "\n";
1238 }
1239 else
1240 {
1241 rbfSplinesNutFluct[i] =
1242 new SPLINTER::RBFSpline
1243 (
1244 *samplesNutFluct[i],
1245 SPLINTER::RadialBasisFunctionType::GAUSSIAN,
1246 false,
1247 radiiFluct(i)
1248 );
1249 ITHACAstream::SaveDenseMatrix
1250 (
1251 rbfSplinesNutFluct[i]->weights,
1252 "./ITHACAoutput/weightsSUP/",
1253 weightName
1254 );
1255 Info << " [SUP] nut_fluct RBF " << i + 1 << "/" << nNutFluctModes
1256 << " fitted & saved to ./ITHACAoutput/weightsSUP/" << weightName << "\n";
1257 }
1258 }
1259
1260 this->rbfSplinesNutAvg = rbfSplinesNutAvg;
1261 this->rbfSplinesNutFluct = rbfSplinesNutFluct;
1262 this->samplesNutAvg = samplesNutAvg;
1263 this->samplesNutFluct = samplesNutFluct;
1264 Info <<
1265 ">>> [SUP] Finished AVG linear-μ table export + FLUCT RBF build.\n";
1266 }
1267}
1268
1269void UnsteadyNSTurb::projectPPE(fileName folder,
1270 label NU,
1271 label NP,
1272 label NSUP,
1273 label Nnut,
1274 bool rbfInterp)
1275{
1276 NUmodes = NU;
1277 NPmodes = NP;
1278 NSUPmodes = 0;
1279 nNutModes = Nnut;
1280 L_U_SUPmodes.resize(0);
1281
1282 if (liftfield.size() != 0)
1283 {
1284 for (label k = 0; k < liftfield.size(); ++k)
1285 {
1286 L_U_SUPmodes.append(tmp<volVectorField>(liftfield[k]));
1287 }
1288 }
1289
1290 if (NUmodes != 0)
1291 {
1292 for (label k = 0; k < NUmodes; ++k)
1293 {
1294 L_U_SUPmodes.append(tmp<volVectorField>(Umodes[k]));
1295 }
1296 }
1297
1298 if (ITHACAutilities::check_folder("./ITHACAoutput/Matrices/"))
1299 {
1300 {
1301 const word B_str =
1302 "B_" + name(liftfield.size()) + "_" + name(NUmodes) + "_" + name(NSUPmodes);
1303
1304 if (ITHACAutilities::check_file("./ITHACAoutput/Matrices/" + B_str))
1305 {
1306 ITHACAstream::ReadDenseMatrix(B_matrix, "./ITHACAoutput/Matrices/", B_str);
1307 }
1308 else
1309 {
1310 B_matrix = diffusive_term(NUmodes, NPmodes, NSUPmodes);
1311 }
1312 }
1313 {
1314 const word btStr =
1315 "bt_" + name(liftfield.size()) + "_" + name(NUmodes) + "_" + name(NSUPmodes);
1316
1317 if (ITHACAutilities::check_file("./ITHACAoutput/Matrices/" + btStr))
1318 {
1319 ITHACAstream::ReadDenseMatrix(btMatrix, "./ITHACAoutput/Matrices/", btStr);
1320 }
1321 else
1322 {
1323 btMatrix = btTurbulence(NUmodes, NSUPmodes);
1324 }
1325 }
1326 {
1327 const word K_str =
1328 "K_" + name(liftfield.size()) + "_" + name(NUmodes) + "_" + name(
1329 NSUPmodes) + "_" + name(NPmodes);
1330
1331 if (ITHACAutilities::check_file("./ITHACAoutput/Matrices/" + K_str))
1332 {
1333 ITHACAstream::ReadDenseMatrix(K_matrix, "./ITHACAoutput/Matrices/", K_str);
1334 }
1335 else
1336 {
1337 K_matrix = pressure_gradient_term(NUmodes, NPmodes, NSUPmodes);
1338 }
1339 }
1340 {
1341 const word M_str =
1342 "M_" + name(liftfield.size()) + "_" + name(NUmodes) + "_" + name(NSUPmodes);
1343
1344 if (ITHACAutilities::check_file("./ITHACAoutput/Matrices/" + M_str))
1345 {
1346 ITHACAstream::ReadDenseMatrix(M_matrix, "./ITHACAoutput/Matrices/", M_str);
1347 }
1348 else
1349 {
1350 M_matrix = mass_term(NUmodes, NPmodes, NSUPmodes);
1351 }
1352 }
1353 {
1354 const word D_str = "D_" + name(NPmodes);
1355
1356 if (ITHACAutilities::check_file("./ITHACAoutput/Matrices/" + D_str))
1357 {
1358 ITHACAstream::ReadDenseMatrix(D_matrix, "./ITHACAoutput/Matrices/", D_str);
1359 }
1360 else
1361 {
1362 D_matrix = laplacian_pressure(NPmodes);
1363 }
1364 }
1365 {
1366 const word bc1_str =
1367 "BC1_" + name(liftfield.size()) + "_" + name(NUmodes) + "_" + name(NPmodes);
1368
1369 if (ITHACAutilities::check_file("./ITHACAoutput/Matrices/" + bc1_str))
1370 {
1371 ITHACAstream::ReadDenseMatrix(BC1_matrix, "./ITHACAoutput/Matrices/", bc1_str);
1372 }
1373 else
1374 {
1375 BC1_matrix = pressure_BC1(NUmodes, NPmodes);
1376 }
1377 }
1378 {
1379 const word bc2_str =
1380 "BC2_" + name(liftfield.size()) + "_" + name(NUmodes) + "_" + name(
1381 NSUPmodes) + "_" + name(NPmodes) + "_t";
1382
1383 if (ITHACAutilities::check_file("./ITHACAoutput/Matrices/" + bc2_str))
1384 {
1385 ITHACAstream::ReadDenseTensor(bc2Tensor, "./ITHACAoutput/Matrices/", bc2_str);
1386 }
1387 else
1388 {
1389 bc2Tensor = pressureBC2(NUmodes, NPmodes);
1390 }
1391 }
1392 {
1393 const word bc3_str =
1394 "BC3_" + name(liftfield.size()) + "_" + name(NUmodes) + "_" + name(NPmodes);
1395
1396 if (ITHACAutilities::check_file("./ITHACAoutput/Matrices/" + bc3_str))
1397 {
1398 ITHACAstream::ReadDenseMatrix(BC3_matrix, "./ITHACAoutput/Matrices/", bc3_str);
1399 }
1400 else
1401 {
1402 BC3_matrix = pressure_BC3(NUmodes, NPmodes);
1403 }
1404 }
1405 {
1406 const word C_str =
1407 "C_" + name(liftfield.size()) + "_" + name(NUmodes) + "_" + name(
1408 NSUPmodes) + "_t";
1409
1410 if (ITHACAutilities::check_file("./ITHACAoutput/Matrices/" + C_str))
1411 {
1412 ITHACAstream::ReadDenseTensor(C_tensor, "./ITHACAoutput/Matrices/", C_str);
1413 }
1414 else
1415 {
1416 C_tensor = convective_term_tens(NUmodes, NPmodes, NSUPmodes);
1417 }
1418 }
1419 {
1420 const word ct1Str =
1421 "ct1_" + name(liftfield.size()) + "_" + name(NUmodes) + "_" + name(
1422 NSUPmodes) + "_"
1423 + name(nNutModes) + "_t";
1424
1425 if (ITHACAutilities::check_file("./ITHACAoutput/Matrices/" + ct1Str))
1426 {
1427 ITHACAstream::ReadDenseTensor(ct1Tensor, "./ITHACAoutput/Matrices/", ct1Str);
1428 }
1429 else
1430 {
1431 ct1Tensor = turbulenceTensor1(NUmodes, NSUPmodes, nNutModes);
1432 }
1433
1434 const word ct2Str =
1435 "ct2_" + name(liftfield.size()) + "_" + name(NUmodes) + "_" + name(
1436 NSUPmodes) + "_"
1437 + name(nNutModes) + "_t";
1438
1439 if (ITHACAutilities::check_file("./ITHACAoutput/Matrices/" + ct2Str))
1440 {
1441 ITHACAstream::ReadDenseTensor(ct2Tensor, "./ITHACAoutput/Matrices/", ct2Str);
1442 }
1443 else
1444 {
1445 ct2Tensor = turbulenceTensor2(NUmodes, NSUPmodes, nNutModes);
1446 }
1447 }
1448 {
1449 const word G_str =
1450 "G_" + name(liftfield.size()) + "_" + name(NUmodes) + "_" + name(
1451 NSUPmodes) + "_" + name(NPmodes) + "_t";
1452
1453 if (ITHACAutilities::check_file("./ITHACAoutput/Matrices/" + G_str))
1454 {
1455 ITHACAstream::ReadDenseTensor(gTensor, "./ITHACAoutput/Matrices/", G_str);
1456 }
1457 else
1458 {
1459 gTensor = divMomentum(NUmodes, NPmodes);
1460 }
1461 }
1462
1463 if (nutAve.size() != 0)
1464 {
1465 const word ct1AveStr =
1466 "ct1Ave_" + name(liftfield.size()) + "_" + name(NUmodes) + "_" + name(
1467 NSUPmodes) + "_t";
1468
1469 if (ITHACAutilities::check_file("./ITHACAoutput/Matrices/" + ct1AveStr))
1470 {
1471 ITHACAstream::ReadDenseTensor(ct1AveTensor, "./ITHACAoutput/Matrices/",
1472 ct1AveStr);
1473 }
1474 else
1475 {
1476 ct1AveTensor = turbulenceAveTensor1(NUmodes, NSUPmodes);
1477 }
1478
1479 const word ct2AveStr =
1480 "ct2Ave_" + name(liftfield.size()) + "_" + name(NUmodes) + "_" + name(
1481 NSUPmodes) + "_t";
1482
1483 if (ITHACAutilities::check_file("./ITHACAoutput/Matrices/" + ct2AveStr))
1484 {
1485 ITHACAstream::ReadDenseTensor(ct2AveTensor, "./ITHACAoutput/Matrices/",
1486 ct2AveStr);
1487 }
1488 else
1489 {
1490 ct2AveTensor = turbulenceAveTensor2(NUmodes, NSUPmodes);
1491 }
1492 }
1493
1494 if (nutFluctModes.size() != 0)
1495 {
1496 const word ct1FluctStr =
1497 "ct1Fluct_" + name(liftfield.size()) + "_" + name(NUmodes) + "_" + name(
1498 NSUPmodes) + "_t";
1499
1500 if (ITHACAutilities::check_file("./ITHACAoutput/Matrices/" + ct1FluctStr))
1501 {
1502 ITHACAstream::ReadDenseTensor(ct1FluctTensor, "./ITHACAoutput/Matrices/",
1503 ct1FluctStr);
1504 }
1505 else
1506 {
1507 ct1FluctTensor = turbulenceFluctTensor1(NUmodes, NSUPmodes);
1508 }
1509
1510 const word ct2FluctStr =
1511 "ct2Fluct_" + name(liftfield.size()) + "_" + name(NUmodes) + "_" + name(
1512 NSUPmodes) + "_t";
1513
1514 if (ITHACAutilities::check_file("./ITHACAoutput/Matrices/" + ct2FluctStr))
1515 {
1516 ITHACAstream::ReadDenseTensor(ct2FluctTensor, "./ITHACAoutput/Matrices/",
1517 ct2FluctStr);
1518 }
1519 else
1520 {
1521 ct2FluctTensor = turbulenceFluctTensor2(NUmodes, NSUPmodes);
1522 }
1523 }
1524
1525 {
1526 const word ct1PPEStr =
1527 "ct1PPE_" + name(liftfield.size()) + "_" + name(NUmodes) + "_" + name(
1528 NSUPmodes) + "_"
1529 + name(NPmodes) + "_" + name(nNutModes) + "_t";
1530
1531 if (ITHACAutilities::check_file("./ITHACAoutput/Matrices/" + ct1PPEStr))
1532 {
1533 ITHACAstream::ReadDenseTensor(ct1PPETensor, "./ITHACAoutput/Matrices/",
1534 ct1PPEStr);
1535 }
1536 else
1537 {
1538 ct1PPETensor = turbulencePPETensor1(NUmodes, NSUPmodes, NPmodes, nNutModes);
1539 }
1540
1541 const word ct2PPEStr =
1542 "ct2PPE_" + name(liftfield.size()) + "_" + name(NUmodes) + "_" + name(
1543 NSUPmodes) + "_"
1544 + name(NPmodes) + "_" + name(nNutModes) + "_t";
1545
1546 if (ITHACAutilities::check_file("./ITHACAoutput/Matrices/" + ct2PPEStr))
1547 {
1548 ITHACAstream::ReadDenseTensor(ct2PPETensor, "./ITHACAoutput/Matrices/",
1549 ct2PPEStr);
1550 }
1551 else
1552 {
1553 ct2PPETensor = turbulencePPETensor2(NUmodes, NSUPmodes, NPmodes, nNutModes);
1554 }
1555 }
1556
1557 if (nutAve.size() != 0)
1558 {
1559 const word ct1PPEAveStr =
1560 "ct1PPEAve_" + name(liftfield.size()) + "_" + name(NUmodes) + "_" + name(
1561 NSUPmodes) + "_"
1562 + name(NPmodes) + "_t";
1563
1564 if (ITHACAutilities::check_file("./ITHACAoutput/Matrices/" + ct1PPEAveStr))
1565 {
1566 ITHACAstream::ReadDenseTensor(ct1PPEAveTensor, "./ITHACAoutput/Matrices/",
1567 ct1PPEAveStr);
1568 }
1569 else
1570 {
1571 ct1PPEAveTensor = turbulencePPEAveTensor1(NUmodes, NSUPmodes, NPmodes);
1572 }
1573
1574 const word ct2PPEAveStr =
1575 "ct2PPEAve_" + name(liftfield.size()) + "_" + name(NUmodes) + "_" + name(
1576 NSUPmodes) + "_"
1577 + name(NPmodes) + "_t";
1578
1579 if (ITHACAutilities::check_file("./ITHACAoutput/Matrices/" + ct2PPEAveStr))
1580 {
1581 ITHACAstream::ReadDenseTensor(ct2PPEAveTensor, "./ITHACAoutput/Matrices/",
1582 ct2PPEAveStr);
1583 }
1584 else
1585 {
1586 ct2PPEAveTensor = turbulencePPEAveTensor2(NUmodes, NSUPmodes, NPmodes);
1587 }
1588 }
1589
1590 if (nutFluctModes.size() != 0)
1591 {
1592 const word ct1PPEFluctStr =
1593 "ct1PPEFluct_" + name(liftfield.size()) + "_" + name(NUmodes) + "_" + name(
1594 NSUPmodes) + "_"
1595 + name(NPmodes) + "_t";
1596
1597 if (ITHACAutilities::check_file("./ITHACAoutput/Matrices/" + ct1PPEFluctStr))
1598 {
1599 ITHACAstream::ReadDenseTensor(ct1PPEFluctTensor, "./ITHACAoutput/Matrices/",
1600 ct1PPEFluctStr);
1601 }
1602 else
1603 {
1604 ct1PPEFluctTensor = turbulencePPEFluctTensor1(NUmodes, NSUPmodes, NPmodes);
1605 }
1606
1607 const word ct2PPEFluctStr =
1608 "ct2PPEFluct_" + name(liftfield.size()) + "_" + name(NUmodes) + "_" + name(
1609 NSUPmodes) + "_"
1610 + name(NPmodes) + "_t";
1611
1612 if (ITHACAutilities::check_file("./ITHACAoutput/Matrices/" + ct2PPEFluctStr))
1613 {
1614 ITHACAstream::ReadDenseTensor(ct2PPEFluctTensor, "./ITHACAoutput/Matrices/",
1615 ct2PPEFluctStr);
1616 }
1617 else
1618 {
1619 ct2PPEFluctTensor = turbulencePPEFluctTensor2(NUmodes, NSUPmodes, NPmodes);
1620 }
1621 }
1622
1623 if (NPmodes > 0)
1624 {
1625 const word L_str = "L_" + name(NPmodes);
1626
1627 if (ITHACAutilities::check_file("./ITHACAoutput/Matrices/" + L_str))
1628 {
1629 ITHACAstream::ReadDenseMatrix(L_vector, "./ITHACAoutput/Matrices/", L_str);
1630 }
1631 else
1632 {
1633 L_vector = pressurePPE_L(NPmodes);
1634 }
1635 }
1636
1637 if (bcMethod == "penalty")
1638 {
1639 bcVelVec = bcVelocityVec(NUmodes, NSUPmodes);
1640 bcVelMat = bcVelocityMat(NUmodes, NSUPmodes);
1641 }
1642 }
1643 else
1644 {
1645 B_matrix = diffusive_term(NUmodes, NPmodes, NSUPmodes);
1646 C_tensor = convective_term_tens(NUmodes, NPmodes, NSUPmodes);
1647 M_matrix = mass_term(NUmodes, NPmodes, NSUPmodes);
1648 K_matrix = pressure_gradient_term(NUmodes, NPmodes, NSUPmodes);
1649 D_matrix = laplacian_pressure(NPmodes);
1650 gTensor = divMomentum(NUmodes, NPmodes);
1651 BC1_matrix = pressure_BC1(NUmodes, NPmodes);
1652 bc2Tensor = pressureBC2(NUmodes, NPmodes);
1653 BC3_matrix = pressure_BC3(NUmodes, NPmodes);
1654 btMatrix = btTurbulence(NUmodes, NSUPmodes);
1655 ct1Tensor = turbulenceTensor1(NUmodes, NSUPmodes, nNutModes);
1656 ct2Tensor = turbulenceTensor2(NUmodes, NSUPmodes, nNutModes);
1657 ct1PPETensor = turbulencePPETensor1(NUmodes, NSUPmodes, NPmodes, nNutModes);
1658 ct2PPETensor = turbulencePPETensor2(NUmodes, NSUPmodes, NPmodes, nNutModes);
1659
1660 if (nutAve.size() != 0)
1661 {
1662 ct1AveTensor = turbulenceAveTensor1(NUmodes, NSUPmodes);
1663 ct2AveTensor = turbulenceAveTensor2(NUmodes, NSUPmodes);
1664 ct1PPEAveTensor = turbulencePPEAveTensor1(NUmodes, NSUPmodes, NPmodes);
1665 ct2PPEAveTensor = turbulencePPEAveTensor2(NUmodes, NSUPmodes, NPmodes);
1666 }
1667
1668 if (nutFluctModes.size() != 0)
1669 {
1670 ct1FluctTensor = turbulenceFluctTensor1(NUmodes, NSUPmodes);
1671 ct2FluctTensor = turbulenceFluctTensor2(NUmodes, NSUPmodes);
1672 ct1PPEFluctTensor = turbulencePPEFluctTensor1(NUmodes, NSUPmodes, NPmodes);
1673 ct2PPEFluctTensor = turbulencePPEFluctTensor2(NUmodes, NSUPmodes, NPmodes);
1674 }
1675
1676 if (NPmodes > 0)
1677 {
1678 L_vector = pressurePPE_L(NPmodes);
1679 }
1680
1681 if (bcMethod == "penalty")
1682 {
1683 bcVelVec = bcVelocityVec(NUmodes, NSUPmodes);
1684 bcVelMat = bcVelocityMat(NUmodes, NSUPmodes);
1685 }
1686 }
1687
1688 if (para->exportPython)
1689 {
1690 ITHACAstream::exportMatrix(B_matrix, "B", "python",
1691 "./ITHACAoutput/Matrices/");
1692 ITHACAstream::exportMatrix(K_matrix, "K", "python",
1693 "./ITHACAoutput/Matrices/");
1694 ITHACAstream::exportMatrix(D_matrix, "D", "python",
1695 "./ITHACAoutput/Matrices/");
1696 ITHACAstream::exportMatrix(M_matrix, "M", "python",
1697 "./ITHACAoutput/Matrices/");
1698 ITHACAstream::exportMatrix(BC1_matrix, "BC1", "python",
1699 "./ITHACAoutput/Matrices/");
1700 ITHACAstream::exportMatrix(BC3_matrix, "BC3", "python",
1701 "./ITHACAoutput/Matrices/");
1702 ITHACAstream::exportTensor(C_tensor, "C", "python",
1703 "./ITHACAoutput/Matrices/");
1704 ITHACAstream::exportTensor(gTensor, "G", "python",
1705 "./ITHACAoutput/Matrices/");
1706 ITHACAstream::exportTensor(bc2Tensor, "BC2", "python",
1707 "./ITHACAoutput/Matrices/");
1708 ITHACAstream::exportTensor(ct1Tensor, "ct1", "python",
1709 "./ITHACAoutput/Matrices/");
1710 ITHACAstream::exportTensor(ct2Tensor, "ct2", "python",
1711 "./ITHACAoutput/Matrices/");
1712 ITHACAstream::exportTensor(ct1PPETensor, "ct1PPE", "python",
1713 "./ITHACAoutput/Matrices/");
1714 ITHACAstream::exportTensor(ct2PPETensor, "ct2PPE", "python",
1715 "./ITHACAoutput/Matrices/");
1716
1717 if (nutAve.size() != 0)
1718 {
1719 ITHACAstream::exportTensor(ct1AveTensor, "ct1Ave", "python",
1720 "./ITHACAoutput/Matrices/");
1721 ITHACAstream::exportTensor(ct2AveTensor, "ct2Ave", "python",
1722 "./ITHACAoutput/Matrices/");
1723 ITHACAstream::exportTensor(ct1PPEAveTensor, "ct1PPEAve", "python",
1724 "./ITHACAoutput/Matrices/");
1725 ITHACAstream::exportTensor(ct2PPEAveTensor, "ct2PPEAve", "python",
1726 "./ITHACAoutput/Matrices/");
1727 }
1728
1729 if (nutFluctModes.size() != 0)
1730 {
1731 ITHACAstream::exportTensor(ct1FluctTensor, "ct1Fluct", "python",
1732 "./ITHACAoutput/Matrices/");
1733 ITHACAstream::exportTensor(ct2FluctTensor, "ct2Fluct", "python",
1734 "./ITHACAoutput/Matrices/");
1735 ITHACAstream::exportTensor(ct1PPEFluctTensor, "ct1PPEFluct", "python",
1736 "./ITHACAoutput/Matrices/");
1737 ITHACAstream::exportTensor(ct2PPEFluctTensor, "ct2PPEFluct", "python",
1738 "./ITHACAoutput/Matrices/");
1739 }
1740
1741 if (NPmodes > 0)
1742 {
1743 ITHACAstream::exportMatrix(L_vector, "L", "python", "./ITHACAoutput/Matrices/");
1744 }
1745 }
1746
1747 if (para->exportMatlab)
1748 {
1749 ITHACAstream::exportMatrix(B_matrix, "B", "matlab",
1750 "./ITHACAoutput/Matrices/");
1751 ITHACAstream::exportMatrix(K_matrix, "K", "matlab",
1752 "./ITHACAoutput/Matrices/");
1753 ITHACAstream::exportMatrix(D_matrix, "D", "matlab",
1754 "./ITHACAoutput/Matrices/");
1755 ITHACAstream::exportMatrix(M_matrix, "M", "matlab",
1756 "./ITHACAoutput/Matrices/");
1757 ITHACAstream::exportMatrix(BC1_matrix, "BC1", "matlab",
1758 "./ITHACAoutput/Matrices/");
1759 ITHACAstream::exportMatrix(BC3_matrix, "BC3", "matlab",
1760 "./ITHACAoutput/Matrices/");
1761 ITHACAstream::exportTensor(C_tensor, "C", "matlab",
1762 "./ITHACAoutput/Matrices/");
1763 ITHACAstream::exportTensor(gTensor, "G", "matlab",
1764 "./ITHACAoutput/Matrices/");
1765 ITHACAstream::exportTensor(bc2Tensor, "BC2", "matlab",
1766 "./ITHACAoutput/Matrices/");
1767 ITHACAstream::exportTensor(ct1Tensor, "ct1", "matlab",
1768 "./ITHACAoutput/Matrices/");
1769 ITHACAstream::exportTensor(ct2Tensor, "ct2", "matlab",
1770 "./ITHACAoutput/Matrices/");
1771 ITHACAstream::exportTensor(ct1PPETensor, "ct1PPE", "matlab",
1772 "./ITHACAoutput/Matrices/");
1773 ITHACAstream::exportTensor(ct2PPETensor, "ct2PPE", "matlab",
1774 "./ITHACAoutput/Matrices/");
1775
1776 if (nutAve.size() != 0)
1777 {
1778 ITHACAstream::exportTensor(ct1AveTensor, "ct1Ave", "matlab",
1779 "./ITHACAoutput/Matrices/");
1780 ITHACAstream::exportTensor(ct2AveTensor, "ct2Ave", "matlab",
1781 "./ITHACAoutput/Matrices/");
1782 ITHACAstream::exportTensor(ct1PPEAveTensor, "ct1PPEAve", "matlab",
1783 "./ITHACAoutput/Matrices/");
1784 ITHACAstream::exportTensor(ct2PPEAveTensor, "ct2PPEAve", "matlab",
1785 "./ITHACAoutput/Matrices/");
1786 }
1787
1788 if (nutFluctModes.size() != 0)
1789 {
1790 ITHACAstream::exportTensor(ct1FluctTensor, "ct1Fluct", "matlab",
1791 "./ITHACAoutput/Matrices/");
1792 ITHACAstream::exportTensor(ct2FluctTensor, "ct2Fluct", "matlab",
1793 "./ITHACAoutput/Matrices/");
1794 ITHACAstream::exportTensor(ct1PPEFluctTensor, "ct1PPEFluct", "matlab",
1795 "./ITHACAoutput/Matrices/");
1796 ITHACAstream::exportTensor(ct2PPEFluctTensor, "ct2PPEFluct", "matlab",
1797 "./ITHACAoutput/Matrices/");
1798 }
1799
1800 if (NPmodes > 0)
1801 {
1802 ITHACAstream::exportMatrix(L_vector, "L", "matlab", "./ITHACAoutput/Matrices/");
1803 }
1804 }
1805
1806 if (para->exportTxt)
1807 {
1808 ITHACAstream::exportMatrix(B_matrix, "B", "eigen",
1809 "./ITHACAoutput/Matrices/");
1810 ITHACAstream::exportMatrix(K_matrix, "K", "eigen",
1811 "./ITHACAoutput/Matrices/");
1812 ITHACAstream::exportMatrix(D_matrix, "D", "eigen",
1813 "./ITHACAoutput/Matrices/");
1814 ITHACAstream::exportMatrix(M_matrix, "M", "eigen",
1815 "./ITHACAoutput/Matrices/");
1816 ITHACAstream::exportMatrix(BC1_matrix, "BC1", "eigen",
1817 "./ITHACAoutput/Matrices/");
1818 ITHACAstream::exportMatrix(BC3_matrix, "BC3", "eigen",
1819 "./ITHACAoutput/Matrices/");
1820 ITHACAstream::exportTensor(C_tensor, "C", "eigen",
1821 "./ITHACAoutput/Matrices/C");
1822 ITHACAstream::exportTensor(gTensor, "G", "eigen",
1823 "./ITHACAoutput/Matrices/G");
1824 ITHACAstream::exportTensor(bc2Tensor, "BC2_", "eigen",
1825 "./ITHACAoutput/Matrices/BC2");
1826 ITHACAstream::exportMatrix(btMatrix, "bt", "eigen",
1827 "./ITHACAoutput/Matrices/");
1828 ITHACAstream::exportTensor(ct1Tensor, "ct1_", "eigen",
1829 "./ITHACAoutput/Matrices/ct1");
1830 ITHACAstream::exportTensor(ct2Tensor, "ct2_", "eigen",
1831 "./ITHACAoutput/Matrices/ct2");
1832 ITHACAstream::exportTensor(ct1PPETensor, "ct1PPE_", "eigen",
1833 "./ITHACAoutput/Matrices/ct1PPE");
1834 ITHACAstream::exportTensor(ct2PPETensor, "ct2PPE_", "eigen",
1835 "./ITHACAoutput/Matrices/ct2PPE");
1836
1837 if (nutAve.size() != 0)
1838 {
1839 ITHACAstream::exportTensor(ct1AveTensor, "ct1Ave_", "eigen",
1840 "./ITHACAoutput/Matrices/ct1Ave");
1841 ITHACAstream::exportTensor(ct2AveTensor, "ct2Ave_", "eigen",
1842 "./ITHACAoutput/Matrices/ct2Ave");
1843 ITHACAstream::exportTensor(ct1PPEAveTensor, "ct1PPEAve_", "eigen",
1844 "./ITHACAoutput/Matrices/ct1PPEAve");
1845 ITHACAstream::exportTensor(ct2PPEAveTensor, "ct2PPEAve_", "eigen",
1846 "./ITHACAoutput/Matrices/ct2PPEAve");
1847 }
1848
1849 if (nutFluctModes.size() != 0)
1850 {
1851 ITHACAstream::exportTensor(ct1FluctTensor, "ct1Fluct", "eigen",
1852 "./ITHACAoutput/Matrices/");
1853 ITHACAstream::exportTensor(ct2FluctTensor, "ct2Fluct", "eigen",
1854 "./ITHACAoutput/Matrices/");
1855 ITHACAstream::exportTensor(ct1PPEFluctTensor, "ct1PPEFluct", "eigen",
1856 "./ITHACAoutput/Matrices/");
1857 ITHACAstream::exportTensor(ct2PPEFluctTensor, "ct2PPEFluct", "eigen",
1858 "./ITHACAoutput/Matrices/");
1859 }
1860
1861 if (NPmodes > 0)
1862 {
1863 ITHACAstream::exportMatrix(L_vector, "L", "eigen", "./ITHACAoutput/Matrices/");
1864 }
1865 }
1866
1867 bTotalMatrix = B_matrix + btMatrix;
1868 const label cSize = NUmodes + NSUPmodes + liftfield.size();
1869 cTotalTensor.resize(cSize, nNutModes, cSize);
1870 cTotalTensor = ct1Tensor + ct2Tensor;
1871 cTotalPPETensor.resize(NPmodes, nNutModes, cSize);
1872 cTotalPPETensor = ct1PPETensor + ct2PPETensor;
1873
1874 if (nutAve.size() != 0)
1875 {
1876 cTotalAveTensor.resize(cSize, nutAve.size(), cSize);
1877 cTotalAveTensor = ct1AveTensor + ct2AveTensor;
1878 cTotalPPEAveTensor.resize(NPmodes, nutAve.size(), cSize);
1879 cTotalPPEAveTensor = ct1PPEAveTensor + ct2PPEAveTensor;
1880 }
1881
1882 if (nutFluctModes.size() != 0)
1883 {
1884 cTotalFluctTensor.resize(cSize, nutFluctModes.size(), cSize);
1885 cTotalFluctTensor = ct1FluctTensor + ct2FluctTensor;
1886 cTotalPPEFluctTensor.resize(NPmodes, nutFluctModes.size(), cSize);
1887 cTotalPPEFluctTensor = ct1PPEFluctTensor + ct2PPEFluctTensor;
1888 }
1889
1890 if (rbfInterp && (!Pstream::parRun()))
1891 {
1892 Info
1893 << "\n==== [RBF DEBUG] ENTERING OFFLINE RBF CONSTRUCTION (AVG linear-μ + FLUCT RBF) ====\n";
1894 const word coeffDir = "./ITHACAoutput/Coefficients/";
1895 const word muPath = "./ITHACAoutput/Offline/mu_samples_mat.txt";
1896 Eigen::MatrixXd muMat = ITHACAstream::readMatrix(muPath);
1897 Info << "[RBF DEBUG] muMat shape: " << muMat.rows() << " x " <<
1898 muMat.cols() << endl;
1899 const int nPar = 12;
1900 const int nSnapshotsPerPar = 200;
1901 Eigen::VectorXd timeVec = muMat.col(0);
1902 Eigen::VectorXd muVec = muMat.col(1);
1903 Eigen::MatrixXd coeffNutAvg = ITHACAstream::readMatrix(
1904 coeffDir + "Nut_avg_coeffs_mat.txt");
1905 Eigen::MatrixXd coeffNutFluct = ITHACAstream::readMatrix(
1906 coeffDir + "Nut_fluct_coeffs_mat.txt");
1907 Info << "[RBF DEBUG] Loaded coeffNutAvg shape: " << coeffNutAvg.rows() <<
1908 " x " << coeffNutAvg.cols() << endl;
1909 Info << "[RBF DEBUG] Loaded coeffNutFluct shape: " << coeffNutFluct.rows()
1910 << " x " << coeffNutFluct.cols() << endl;
1911 const int nNutAvgModes = coeffNutAvg.rows();
1912 const int nNutFluctModes = coeffNutFluct.rows();
1913 const int nUniqueMu = nPar;
1914 Eigen::MatrixXd a = ITHACAutilities::getCoeffs(Uomfield, Umodes);
1915 a.transposeInPlace();
1916 Eigen::VectorXd initSnapInd(nPar);
1917 Eigen::VectorXd timeSnap(nPar);
1918
1919 for (int i = 0; i < nPar; ++i)
1920 {
1921 const int start = i * nSnapshotsPerPar;
1922 initSnapInd(i) = start;
1923 timeSnap(i) = timeVec(start + 1) - timeVec(start);
1924 Info << "[RBF DEBUG] i=" << i << ", start=" << start << ", dt=" <<
1925 timeSnap(i) << endl;
1926 }
1927
1928 Eigen::VectorXd muVecUnique(nUniqueMu);
1929
1930 for (int i = 0; i < nUniqueMu; ++i)
1931 {
1932 muVecUnique(i) = muVec(static_cast<int>(initSnapInd(i)));
1933 }
1934
1935 Info << "[RBF DEBUG] muVecUnique: [" << muVecUnique(0) << " ... "
1936 << muVecUnique(muVecUnique.size() - 1) << "] (M=" << muVecUnique.size() <<
1937 ")\n";
1938 Info <<
1939 "[RBF DEBUG] Calling velDerivativeCoeff() for fluctuation part...\n";
1940 Eigen::MatrixXd Gfluct = coeffNutFluct.transpose();
1941 List<Eigen::MatrixXd> interpDataFluct = velDerivativeCoeff(a, Gfluct,
1942 initSnapInd, timeSnap);
1943 Eigen::MatrixXd velRBF_fluct = interpDataFluct[0];
1944 Eigen::MatrixXd coeffs_fluct = interpDataFluct[1];
1945 Info << "[RBF DEBUG] velRBF_fluct shape: " << velRBF_fluct.rows() << " x "
1946 << velRBF_fluct.cols() << endl;
1947 Info << "[RBF DEBUG] coeffs_fluct shape: " << coeffs_fluct.rows() << " x "
1948 << coeffs_fluct.cols() << endl;
1949 const int nSnapshots = velRBF_fluct.rows();
1950 const int nModes = velRBF_fluct.cols() / 2;
1951 Eigen::MatrixXd adot_only(nSnapshots, nModes);
1952
1953 for (int i = 0; i < nSnapshots; ++i)
1954 {
1955 adot_only.row(i) = velRBF_fluct.row(i).tail(nModes);
1956 }
1957
1958 ITHACAstream::exportMatrix(adot_only, "adot_coeffs", "python", coeffDir);
1959 ITHACAstream::exportMatrix(velRBF_fluct, "a_adot_concat", "python", coeffDir);
1960 ITHACAstream::exportMatrix(velRBF_fluct, "a_adot_concat", "eigen", coeffDir);
1961 Info << "\n[RBF DEBUG] First row of a (velocity coeffs): ";
1962
1963 for (int k = 0; k < a.cols(); ++k)
1964 {
1965 Info << a(0, k) << " ";
1966 }
1967
1968 Info << endl;
1969 Info << "[RBF DEBUG] First row of velRBF_fluct ([a, aDot]): ";
1970
1971 for (int k = 0; k < velRBF_fluct.cols(); ++k)
1972 {
1973 Info << velRBF_fluct(0, k) << " ";
1974 }
1975
1976 Info << endl;
1977 Info << "[RBF DEBUG] First row of adot_only (aDot): ";
1978
1979 for (int k = 0; k < adot_only.cols(); ++k)
1980 {
1981 Info << adot_only(0, k) << " ";
1982 }
1983
1984 Info << endl;
1985 const double eAvg = 3;
1986 const double eFluct = 0.05;
1987 Eigen::MatrixXd radiiAvg;
1988 Eigen::MatrixXd radiiFluct;
1989
1990 if (ITHACAutilities::check_file("./radii_avg.txt"))
1991 {
1992 radiiAvg = ITHACAstream::readMatrix("./radii_avg.txt");
1993 Info << "[RBF DEBUG] Loaded radii_avg.txt, shape: " << radiiAvg.rows() <<
1994 " x " << radiiAvg.cols() << endl;
1995 M_Assert(radiiAvg.size() == nNutAvgModes, "radiiAvg size mismatch");
1996 }
1997 else
1998 {
1999 radiiAvg = Eigen::MatrixXd::Ones(nNutAvgModes, 1) * eAvg;
2000 Info <<
2001 "[RBF DEBUG] Set default radii for nutAvg (unused with linear μ), e=" << eAvg
2002 << endl;
2003 }
2004
2005 if (ITHACAutilities::check_file("./radii_fluct.txt"))
2006 {
2007 radiiFluct = ITHACAstream::readMatrix("./radii_fluct.txt");
2008 Info << "[RBF DEBUG] Loaded radii_fluct.txt, shape: " << radiiFluct.rows()
2009 << " x " << radiiFluct.cols() << endl;
2010 M_Assert(radiiFluct.size() == nNutFluctModes, "radiiFluct size mismatch");
2011 }
2012 else
2013 {
2014 radiiFluct = Eigen::MatrixXd::Ones(nNutFluctModes, 1) * eFluct;
2015 Info << "[RBF DEBUG] Set default radii for nutFluct, e=" << eFluct <<
2016 endl;
2017 }
2018
2019 List<SPLINTER::DataTable*> samplesNutAvg;
2020 List<SPLINTER::RBFSpline*> rbfSplinesNutAvg;
2021 samplesNutAvg.resize(0);
2022 rbfSplinesNutAvg.resize(0);
2023 ITHACAstream::exportMatrix(muVecUnique, "NutAvg_mu_unique", "eigen",
2024 coeffDir);
2025 ITHACAstream::exportMatrix(coeffNutAvg, "NutAvg_coeffs_by_mu", "eigen",
2026 coeffDir);
2027 Info << ">>> Persisted νt_avg tables for linear μ-interpolation: "
2028 << "mu size=" << muVecUnique.size()
2029 << ", coeff table=" << coeffNutAvg.rows() << "x" << coeffNutAvg.cols() << "\n";
2030 List<SPLINTER::DataTable*> samplesNutFluct;
2031 List<SPLINTER::RBFSpline*> rbfSplinesNutFluct;
2032 samplesNutFluct.resize(nNutFluctModes);
2033 rbfSplinesNutFluct.resize(nNutFluctModes);
2034 Info << ">>> Building nut_fluct RBF splines...\n";
2035
2036 for (label i = 0; i < nNutFluctModes; ++i)
2037 {
2038 const word weightName = "wRBF_NUTFLUCT_" + name(i + 1);
2039 samplesNutFluct[i] = new SPLINTER::DataTable(velRBF_fluct.cols(), 1);
2040
2041 for (label j = 0; j < velRBF_fluct.rows(); ++j)
2042 {
2043 samplesNutFluct[i]->addSample(velRBF_fluct.row(j), coeffs_fluct(j, i));
2044 }
2045
2046 Eigen::MatrixXd weights;
2047
2048 if (ITHACAutilities::check_file("./ITHACAoutput/weightsPPE/" + weightName))
2049 {
2050 ITHACAstream::ReadDenseMatrix(weights, "./ITHACAoutput/weightsPPE/",
2051 weightName);
2052 rbfSplinesNutFluct[i] =
2053 new SPLINTER::RBFSpline
2054 (
2055 *samplesNutFluct[i],
2056 SPLINTER::RadialBasisFunctionType::GAUSSIAN,
2057 weights,
2058 radiiFluct(i)
2059 );
2060 Info << " nut_fluct RBF " << i + 1 << "/" << nNutFluctModes
2061 << " loaded from weights.\n";
2062 }
2063 else
2064 {
2065 rbfSplinesNutFluct[i] =
2066 new SPLINTER::RBFSpline
2067 (
2068 *samplesNutFluct[i],
2069 SPLINTER::RadialBasisFunctionType::GAUSSIAN,
2070 false,
2071 radiiFluct(i)
2072 );
2073 ITHACAstream::SaveDenseMatrix
2074 (
2075 rbfSplinesNutFluct[i]->weights,
2076 "./ITHACAoutput/weightsPPE/",
2077 weightName
2078 );
2079 Info << " nut_fluct RBF " << i + 1 << "/" << nNutFluctModes
2080 << " fitted & saved.\n";
2081 }
2082 }
2083
2084 Info << "[OFFLINE] Built nutAvgSplines: " << rbfSplinesNutAvg.size()
2085 << " (expected 0 for linear μ)\n";
2086 Info << "[OFFLINE] Built nutFluctSplines:" << rbfSplinesNutFluct.size() <<
2087 endl;
2088 Info << ">>> Finished AVG linear-μ table export + FLUCT RBF build.\n";
2089 this->rbfSplinesNutAvg = rbfSplinesNutAvg;
2090 this->rbfSplinesNutFluct = rbfSplinesNutFluct;
2091 this->samplesNutAvg = samplesNutAvg;
2092 this->samplesNutFluct = samplesNutFluct;
2093 Info << "[OFFLINE] Built nutAvgSplines: " << rbfSplinesNutAvg.size() <<
2094 endl;
2095 Info << "[OFFLINE] Built nutFluctSplines: " << rbfSplinesNutFluct.size() <<
2096 endl;
2097 }
2098}
2099
2100List<Eigen::MatrixXd> UnsteadyNSTurb::velDerivativeCoeff(Eigen::MatrixXd A,
2101 Eigen::MatrixXd G,
2102 Eigen::VectorXd initSnapInd,
2103 Eigen::VectorXd timeSnap)
2104{
2105 List<Eigen::MatrixXd> newCoeffs;
2106 newCoeffs.setSize(2);
2107 const label velCoeffsNum = A.cols();
2108 const label snapshotsNum = A.rows();
2109 const label parsSamplesNum = initSnapInd.size();
2110 const label timeSnapshotsPerSample = snapshotsNum / parsSamplesNum;
2111 const label newColsNum = 2 * velCoeffsNum;
2112 const label newRowsNum = snapshotsNum - parsSamplesNum;
2113 newCoeffs[0].resize(newRowsNum, newColsNum);
2114 newCoeffs[1].resize(newRowsNum, G.cols());
2115 int rowCount = 0;
2116
2117 for (label j = 0; j < parsSamplesNum; ++j)
2118 {
2119 const int i0 = j * timeSnapshotsPerSample;
2120 const int N = timeSnapshotsPerSample;
2121
2122 for (int n = 1; n < N; ++n, ++rowCount)
2123 {
2124 const Eigen::RowVectorXd a_now = A.row(i0 + n);
2125 const Eigen::RowVectorXd a_prev = A.row(i0 + n - 1);
2126 const Eigen::RowVectorXd adot = (a_now - a_prev) / timeSnap(j);
2127 newCoeffs[0].row(rowCount) << a_now, adot;
2128 newCoeffs[1].row(rowCount) = G.row(i0 + n);
2129 }
2130 }
2131
2132 interChoice = 3;
2133 return newCoeffs;
2134}
2135
2136List<Eigen::MatrixXd> UnsteadyNSTurb::velParCoeff(Eigen::MatrixXd A,
2137 Eigen::MatrixXd G)
2138{
2139 List<Eigen::MatrixXd> newCoeffs;
2140 newCoeffs.setSize(2);
2141 Eigen::MatrixXd pars = z.leftCols(z.cols() - 1);
2142 newCoeffs[0].resize(A.rows(), A.cols() + z.cols() - 1);
2143 newCoeffs[1].resize(G.rows(), G.cols());
2144 newCoeffs[0] << pars, A;
2145 newCoeffs[1] = G;
2146 interChoice = 2;
2147 return newCoeffs;
2148}
2149
2150List<Eigen::MatrixXd> UnsteadyNSTurb::velParDerivativeCoeff(Eigen::MatrixXd A,
2151 Eigen::MatrixXd G,
2152 Eigen::VectorXd initSnapInd,
2153 Eigen::VectorXd timeSnap)
2154{
2155 Info << "[velParDerivativeCoeff] ENTER" << endl;
2156 List<Eigen::MatrixXd> newCoeffs;
2157 newCoeffs.setSize(2);
2158 const label velCoeffsNum = A.cols();
2159 const label snapshotsNum = A.rows();
2160 const label parsSamplesNum = initSnapInd.size();
2161 const label timeSnapshotsPerSample = snapshotsNum / parsSamplesNum;
2162 Info << "[velParDerivativeCoeff] A shape: " << snapshotsNum << " x " <<
2163 velCoeffsNum << endl;
2164 Info << "[velParDerivativeCoeff] G shape: " << G.rows() << " x " <<
2165 G.cols() << endl;
2166 Info << "[velParDerivativeCoeff] nPars: " << parsSamplesNum
2167 << ", timeSnapshotsPerSample: " << timeSnapshotsPerSample << endl;
2168 Eigen::MatrixXd pars(snapshotsNum, 1);
2169
2170 for (label j = 0; j < parsSamplesNum; ++j)
2171 {
2172 pars.block(j * timeSnapshotsPerSample, 0, timeSnapshotsPerSample, 1) =
2173 Eigen::VectorXd::Constant(timeSnapshotsPerSample, mu(j));
2174 }
2175
2176 const label newColsNum = 2 * velCoeffsNum;
2177 const label newRowsNum = snapshotsNum - parsSamplesNum;
2178 Info << "[velParDerivativeCoeff] newCoeffs[0] size: " << newRowsNum <<
2179 " x " << (newColsNum + pars.cols()) << endl;
2180 Info << "[velParDerivativeCoeff] newCoeffs[1] size: " << newRowsNum <<
2181 " x " << G.cols() << endl;
2182 newCoeffs[0].resize(newRowsNum, newColsNum + pars.cols());
2183 newCoeffs[1].resize(newRowsNum, G.cols());
2184 int totalRowsWritten = 0;
2185
2186 for (label j = 0; j < parsSamplesNum; ++j)
2187 {
2188 const int start = j * timeSnapshotsPerSample;
2189 const int rowOffset = j * (timeSnapshotsPerSample - 1);
2190 Info << "[velParDerivativeCoeff] Group " << j
2191 << ": start=" << start << ", rowOffset=" << rowOffset << endl;
2192 Info << "[velParDerivativeCoeff] b0: rows " << start << " to "
2193 << (start + timeSnapshotsPerSample - 2) << endl;
2194
2195 if (start + timeSnapshotsPerSample - 1 > snapshotsNum)
2196 {
2197 std::cerr << "[velParDerivativeCoeff][ERROR] b0 access out of bounds! (start="
2198 << start << ", timeSnapshotsPerSample-1=" << (timeSnapshotsPerSample - 1)
2199 << ", snapshotsNum=" << snapshotsNum << ")" << endl;
2200 abort();
2201 }
2202
2203 const Eigen::MatrixXd b0 = A.middleRows(start, timeSnapshotsPerSample - 1);
2204 Info << "[velParDerivativeCoeff] b1: rows " << (start + 1) << " to "
2205 << (start + timeSnapshotsPerSample - 1) << endl;
2206
2207 if (start + 1 + timeSnapshotsPerSample - 2 >= snapshotsNum)
2208 {
2209 std::cerr << "[velParDerivativeCoeff][ERROR] b1 access out of bounds! (start+1="
2210 << (start + 1) << ", timeSnapshotsPerSample-1=" << (timeSnapshotsPerSample - 1)
2211 << ", snapshotsNum=" << snapshotsNum << ")" << endl;
2212 abort();
2213 }
2214
2215 const Eigen::MatrixXd b1 = A.middleRows(start + 1, timeSnapshotsPerSample - 1);
2216 Eigen::MatrixXd bNew(b0.rows(), b0.cols() + b1.cols());
2217 bNew << b1, (b1 - b0) / timeSnap(j);
2218 Info << "[velParDerivativeCoeff] pars block for input: rows " <<
2219 (start + 1) << " to "
2220 << (start + timeSnapshotsPerSample - 1) << endl;
2221
2222 if (start + 1 + timeSnapshotsPerSample - 2 >= snapshotsNum)
2223 {
2224 std::cerr << "[velParDerivativeCoeff][ERROR] pars access out of bounds!" <<
2225 endl;
2226 abort();
2227 }
2228
2229 newCoeffs[0].block(rowOffset, 0, timeSnapshotsPerSample - 1, pars.cols()) =
2230 pars.middleRows(start + 1, timeSnapshotsPerSample - 1);
2231 newCoeffs[0].block(rowOffset, pars.cols(), timeSnapshotsPerSample - 1,
2232 newColsNum) = bNew;
2233 Info << "[velParDerivativeCoeff] G block for output: rows " <<
2234 (start + 1) << " to "
2235 << (start + timeSnapshotsPerSample - 1) << endl;
2236
2237 if (start + 1 + timeSnapshotsPerSample - 2 >= G.rows())
2238 {
2239 std::cerr << "[velParDerivativeCoeff][ERROR] G access out of bounds!" <<
2240 endl;
2241 abort();
2242 }
2243
2244 newCoeffs[1].middleRows(rowOffset, timeSnapshotsPerSample - 1) =
2245 G.middleRows(start + 1, timeSnapshotsPerSample - 1);
2246 totalRowsWritten += timeSnapshotsPerSample - 1;
2247 Info << "[velParDerivativeCoeff] Finished group " << j
2248 << ", totalRowsWritten so far: " << totalRowsWritten << endl;
2249 }
2250
2251 Info << "[velParDerivativeCoeff] FINISHED. Total rows written: "
2252 << totalRowsWritten << "/" << newRowsNum << endl;
2253 interChoice = 4;
2254 return newCoeffs;
2255}
2256
2257Eigen::MatrixXd UnsteadyNSTurb::velParDerivativeCoeff(Eigen::MatrixXd A,
2258 Eigen::VectorXd par,
2259 double timeSnap)
2260{
2261 Eigen::MatrixXd newCoeffs;
2262 const label velCoeffsNum = A.cols();
2263 const label snapshotsNum = A.rows();
2264 const label parsSamplesNum = par.size();
2265 const label newColsNum = 2 * velCoeffsNum + parsSamplesNum;
2266 const label newRowsNum = snapshotsNum - 1;
2267 newCoeffs.resize(newRowsNum, newColsNum);
2268 const Eigen::MatrixXd b0 = A.topRows(A.rows() - 1);
2269 const Eigen::MatrixXd b1 = A.bottomRows(A.rows() - 1);
2270 Eigen::MatrixXd bNew(b0.rows(), b0.cols() + b1.cols());
2271 bNew << b1, ((b1 - b0) / timeSnap);
2272 newCoeffs.leftCols(parsSamplesNum) =
2273 Eigen::MatrixXd::Ones(newRowsNum, parsSamplesNum) * par;
2274 newCoeffs.rightCols(newColsNum - parsSamplesNum) = bNew;
2275 return newCoeffs;
2276}
2277
2278
2279
2280// =========================================================================
2281// ============================= Smagorinsky ===============================
2282// =========================================================================
2283
2284void UnsteadyNSTurb::computeNonLinearSnapshot_at_time(const word& snap_time,
2285 volVectorField& Smag, std::optional<PtrList<volVectorField>> modesU)
2286{
2287 if (m_parameters->get_DEIMInterpolatedField() == "fullStressFunction"
2288 || ITHACAutilities::containsSubstring(m_parameters->get_DEIMInterpolatedField(), "reducedFullStressFunction"))
2289 {
2290 Smag = computeSmagTerm_at_time(snap_time, modesU) ;
2291 }
2292 else
2293 {
2294 Info << "Error : DEIMInterpolatedField not valid : "
2295 << m_parameters->get_DEIMInterpolatedField() << endl;
2296 Info << "DEIM is available for fullStressFunction and nut only." << endl;
2297 abort();
2298 }
2299}
2300
2301// Init Smagorinsky term
2302volVectorField UnsteadyNSTurb::initSmagFunction()
2303{
2304 return volVectorField("fullStressFunction", computeSmagTerm_at_time("0") );
2305}
2306
2307// Init Smagorinsky term
2308volScalarField UnsteadyNSTurb::initSmagPhiFunction(const volScalarField template_field_phi)
2309{
2310 return volScalarField(m_parameters->get_DEIMInterpolatedField(), computeSmagTermPhi_at_time("0",template_field_phi));
2311}
2312
2313volVectorField UnsteadyNSTurb::computeSmagTerm_at_time(const word& snap_time,
2314 std::optional<PtrList<volVectorField>> modesU)
2315{
2316 // Read the j-th field
2317 volVectorField snapshotj = m_parameters->get_template_field_U();
2318 ITHACAstream::read_snapshot(snapshotj, snap_time, m_parameters->get_casenameData());
2319
2320 if (modesU)
2321 {
2322 volVectorField proj_snapshotj = ITHACAutilities::project_to_POD_basis(snapshotj, modesU.value(), m_parameters->get_meanU());
2323 snapshotj = proj_snapshotj;
2324 }
2325
2326 return computeSmagTerm_fromU(snapshotj);
2327}
2328
2329volScalarField UnsteadyNSTurb::computeSmagTermPhi_at_time(const word& snap_time, const volScalarField template_field_phi)
2330{
2331 // Read the j-th field
2332 volScalarField phij = template_field_phi;
2333 ITHACAstream::read_snapshot(phij, snap_time, m_parameters->get_casenameData());
2334 volVectorField snapshotj = m_parameters->get_template_field_U();
2335 ITHACAstream::read_snapshot(snapshotj, snap_time, m_parameters->get_casenameData());
2336 return computeSmagTermPhi_fromUPhi(snapshotj,phij);
2337}
2338
2339volVectorField UnsteadyNSTurb::computeSmagTerm_fromU(const volVectorField& snapshotj)
2340{
2341 volTensorField S=computeS_fromU(snapshotj);
2342 volScalarField nut = computeNut_fromS(S);
2343 return (fvc::div(2*nut*dev(S)));
2344}
2345
2346volScalarField UnsteadyNSTurb::computeSmagTermPhi_fromUPhi(const volVectorField& snapshotj, const volScalarField& phij)
2347{
2348 volTensorField S=computeS_fromU(snapshotj);
2349 volScalarField nut = computeNut_fromS(S);
2350 return (fvc::div(nut*fvc::grad(phij)));
2351}
2352
2353template<typename T>
2354volScalarField UnsteadyNSTurb::diffusion(const T& coefDiff, const volScalarField& phi)
2355{
2356 return ( fvc::laplacian(coefDiff , phi) );
2357}
2358template volScalarField UnsteadyNSTurb::diffusion(const volScalarField& coefDiff, const volScalarField& u);
2359template volScalarField UnsteadyNSTurb::diffusion(const volTensorField& coefDiff, const volScalarField& u);
2360
2361template<typename T>
2362volVectorField UnsteadyNSTurb::diffusion(const T& coefDiff, const volVectorField& u)
2363{
2364 volTensorField S_u=computeS_fromU(u);
2365 return ( fvc::div( 2 * ITHACAutilities::tensorFieldProduct( coefDiff , dev(S_u) ) ) );
2366}
2367template volVectorField UnsteadyNSTurb::diffusion(const volScalarField& coefDiff, const volVectorField& u);
2368template volVectorField UnsteadyNSTurb::diffusion(const volTensorField& coefDiff, const volVectorField& u);
2369
2370volVectorField UnsteadyNSTurb::computeSmagTermPhi_fromUPhi(const volVectorField& snapshotj, const volVectorField& phij)
2371{
2372 return computeSmagTerm_fromU(snapshotj);
2373}
2374
2375
2376// =========================================================================
2377// ======================= projFullStressFunction ==========================
2378// =========================================================================
2379
2380void UnsteadyNSTurb::computeNonLinearSnapshot_at_time(const word& snap_time, volScalarField& phi, volVectorField& modeU)
2381{
2382 if (ITHACAutilities::containsSubstring(m_parameters->get_HRSnapshotsField(), "projFullStressFunction")
2383 || ITHACAutilities::containsSubstring(m_parameters->get_HRSnapshotsField(), "projReducedFullStressFunction"))
2384 {
2385 phi = computeProjSmagTerm_at_time_fromMode(snap_time, modeU);
2386 }
2387}
2388
2389// Init projected Smagorinsky term
2390volScalarField UnsteadyNSTurb::initProjSmagFunction()
2391{
2392 return volScalarField("projFullStressFunction",
2393 m_parameters->get_template_field_fullStressFunction() & m_parameters->get_template_field_U());
2394}
2395
2396volScalarField UnsteadyNSTurb::computeProjSmagTerm_at_time_fromMode(const word& snap_time, const volVectorField& mode)
2397{
2398 // Read the j-th field
2399 volVectorField Smagj = m_parameters->get_template_field_fullStressFunction();
2400 if (!ITHACAutilities::containsSubstring(m_parameters->get_HRSnapshotsField(), "projReducedFullStressFunction"))
2401 {
2402 ITHACAstream::read_snapshot(Smagj, snap_time, m_parameters->get_casenameData());
2403 }
2404 else
2405 {
2406 word path = "./ITHACAoutput/Hyperreduction/reducedFullStressFunction/";
2407 ITHACAstream::read_snapshot(Smagj, snap_time, path);
2408 }
2409
2410 return Smagj & mode;
2411}
2412
2413void UnsteadyNSTurb::computeProjSmagTerm_fromSmag_fromMode(volScalarField& phi, const volVectorField& Smag,
2414 const volVectorField& mode)
2415{
2416 phi = Smag & mode;
2417}
2418
2419
2420// =========================================================================
2421// ========================== projSmagFromNut ==============================
2422// =========================================================================
2423
2424void UnsteadyNSTurb::computeNonLinearSnapshot_at_time(const word& snap_time, volScalarField& phi,
2425 volVectorField& modeU_proj , volVectorField& modeU_grad)
2426{
2427 if (ITHACAutilities::containsSubstring(m_parameters->get_HRSnapshotsField(), "projSmagFromNut")
2428 || ITHACAutilities::containsSubstring(m_parameters->get_HRSnapshotsField(), "projSmagFromReducedNut"))
2429 {
2430 phi = computeProjSmagFromNut_at_time_fromModes(snap_time, modeU_proj, modeU_grad);
2431 }
2432}
2433
2434// Init projected Smagorinsky term from nut
2435volScalarField UnsteadyNSTurb::initProjSmagFromNutFunction()
2436{
2437 return volScalarField("projSmagFromNut", initSmagFunction() & m_parameters->get_template_field_U());
2438}
2439
2440volScalarField UnsteadyNSTurb::computeProjSmagFromNut_at_time_fromModes
2441 (const word& snap_time, const volVectorField& modeU_proj, const volVectorField& modeU_grad)
2442{
2443 // Read the j-th field
2444 volScalarField Nutj(m_parameters->get_DEIMInterpolatedField(), m_parameters->get_template_field_nut());
2445 if (ITHACAutilities::containsSubstring(m_parameters->get_HRSnapshotsField(), "projSmagFromNut") )
2446 {
2447 ITHACAstream::read_snapshot(Nutj, snap_time, m_parameters->get_casenameData());
2448 }
2449 else if (ITHACAutilities::containsSubstring(m_parameters->get_HRSnapshotsField(), "projSmagFromReducedNut"))
2450 {
2451 word path = "./ITHACAoutput/Hyperreduction/reducedNut/";
2452 ITHACAstream::read_snapshot(Nutj, snap_time, path);
2453 }
2454 return projDiffusionIBP(Nutj, modeU_grad, modeU_proj);
2455}
2456
2457void UnsteadyNSTurb::computeProjSmagFromNut_fromNut_fromModes(volScalarField& phi,
2458 const volScalarField& Nut, const volVectorField& modeU_proj, const volVectorField& modeU_grad)
2459{
2460 phi = projDiffusionIBP(Nut, modeU_grad, modeU_proj);
2461}
2462
2463volScalarField UnsteadyNSTurb::projDiffusionIBP(const volScalarField& coefDiff,
2464 const volVectorField& u, const volVectorField& v)
2465{
2466 volTensorField symGradU = computeS_fromU(u);
2467 volTensorField symGradV = computeS_fromU(v);
2468 return 2*coefDiff*(dev(symGradU) && dev(symGradV));
2469}
2470
2471
2472// =========================================================================
2473// ============================= NUT =======================================
2474// =========================================================================
2475
2476void UnsteadyNSTurb::computeNonLinearSnapshot_at_time(const word& snap_time, volScalarField& phi,
2477 std::optional<PtrList<volVectorField>> modesU)
2478{
2479 if (m_parameters->get_DEIMInterpolatedField() == "nut"
2480 || ITHACAutilities::containsSubstring(m_parameters->get_DEIMInterpolatedField(), "reducedNut"))
2481 {
2482 phi = computeNut_at_time(snap_time, modesU) ;
2483 }
2484 else if (m_parameters->get_DEIMInterpolatedField() == "fullStressFunction_K")
2485 {
2486 phi = computeSmagTermPhi_at_time(snap_time, m_parameters->get_template_field_k());
2487 }
2488 else if (m_parameters->get_DEIMInterpolatedField() == "fullStressFunction_Omega")
2489 {
2490 phi = computeSmagTermPhi_at_time(snap_time, m_parameters->get_template_field_omega());
2491 }
2492 else
2493 {
2494 Info << "Error : DEIMInterpolatedField not valid : "
2495 << m_parameters->get_DEIMInterpolatedField() << endl;
2496 Info << "DEIM is available for fullStressFunction and nut only." << endl;
2497 abort();
2498 }
2499}
2500
2501volScalarField UnsteadyNSTurb::initNutFunction()
2502{
2503 return volScalarField(m_parameters->get_template_field_nut());
2504}
2505
2506volScalarField UnsteadyNSTurb::computeNut_at_time(const word& snap_time, std::optional<PtrList<volVectorField>> modesU)
2507{
2508 // Read the j-th field
2509 volVectorField snapshotj = m_parameters->get_template_field_U();
2510 ITHACAstream::read_snapshot(snapshotj, snap_time, m_parameters->get_casenameData());
2511
2512 if (modesU)
2513 {
2514 volVectorField proj_snapshotj = ITHACAutilities::project_to_POD_basis(snapshotj, modesU.value(), m_parameters->get_meanU());
2515 snapshotj = proj_snapshotj;
2516 }
2517
2518 return computeNut_fromU(snapshotj);
2519}
2520
2521volScalarField UnsteadyNSTurb::computeNut_fromU(const volVectorField& snapshotj)
2522{
2523 volTensorField S=computeS_fromU(snapshotj);
2524 return (computeNut_fromS(S));
2525}
2526
2527volScalarField UnsteadyNSTurb::computeNut_fromS(const volTensorField& S)
2528{
2529
2530 volScalarField delta = m_parameters->get_delta();
2531 float Ck = m_parameters->get_Ck();
2532 float Ce = m_parameters->get_Ce();
2533
2534 // // Incompressible flows:
2535 // float Cs = std::pow(Ck*std::pow(Ck/Ce,0.5),0.5);
2536 // return (pow(Cs*delta,2)*sqrt(2*S&&S));
2537
2538 // OpenFOAM like version
2539 // Piece of code strongly inspired by Smagorinsky OpenFOAM Methods
2540 // see https://develop.openfoam.com/Development/openfoam/blob/OpenFOAM-v2012/src/TurbulenceModels/turbulenceModels/LES/Smagorinsky/Smagorinsky.C
2541
2542 volScalarField a(Ce/delta);
2543 volScalarField b((2.0/3.0)*tr(S));
2544 volScalarField c(2*Ck*delta*(dev(S) && S));
2545
2546 volScalarField k(sqr((-b + sqrt(sqr(b) + 4*a*c))/(2*a)));
2547
2548 volScalarField nut = m_parameters->get_template_field_nut();
2549
2550 // Loop over inner values to conserve BC informations
2551 for (int i=0; i<nut.size(); i++)
2552 {
2553 nut[i] = Ck*delta[i]*std::pow(k[i],0.5);
2554 }
2555
2556 // TO DO : make it works
2557 nut.correctBoundaryConditions();
2558 return nut;
2559 }
2560
2561
2562 // void UnsteadyNSTurb::initTurbModel()
2563 // {
2564 // const volVectorField& U(m_parameters->get_template_field_U());
2565 //
2566 // Foam::Time runTime(Foam::Time::controlDictName, ".", m_parameters->get_casenameData());
2567 // wordList boundaryTypeNut = m_parameters->get_template_field_nut().boundaryField().types();
2568 //
2569 // const surfaceScalarField phi
2570 // (
2571 // IOobject
2572 // (
2573 // "phi",
2574 // runTime.timeName(),
2575 // m_parameters->get_mesh(),
2576 // IOobject::READ_IF_PRESENT,
2577 // IOobject::AUTO_WRITE
2578 // ),
2579 // m_parameters->get_mesh(),
2580 // dimensionedScalar("phi", dimensionSet(0,0,0,0,0,0,0), 1.0),
2581 // boundaryTypeNut
2582 // );
2583 //
2584 // const Foam::singlePhaseTransportModel transportModel_(U, phi);
2585 // turbModel = autoPtr<incompressible::turbulenceModel>
2586 // (
2587 // incompressible::turbulenceModel::New(U, phi, transportModel_)
2588 // );
2589 //
2590 // // turbModel->correct();
2591 // // volScalarField nut_test( turbModel->nut());
2592 //
2593 //
2594 // }
2595
2596
2597
2598 volTensorField UnsteadyNSTurb::computeS_fromU(const volVectorField& snapshotj)
2599 {
2600 volTensorField gradV=fvc::grad(snapshotj);
2601 return 0.5*(gradV+gradV.T());
2602 }
2603
2604 volVectorField UnsteadyNSTurb::computeS_fromU(const volScalarField& phij)
2605 {
2606 volVectorField gradV=fvc::grad(phij);
2607 return gradV;
2608 }
ITHACAparameters::getInstance
static ITHACAparameters * getInstance()
Gets an instance of ITHACAparameters, to be used if the instance is already existing.
Definition ITHACAparameters.C:56
StoredParameters
Class that contains all parameters of the stochastic POD.
Definition StoredParameters.H:37
UnsteadyNSTurb::computeSmagTerm_fromU
volVectorField computeSmagTerm_fromU(const volVectorField &snapshotj)
Compute Smagorinsky Snapshot for a given velocity field U.
Definition UnsteadyNSTurb.C:2339
UnsteadyNSTurb::computeProjSmagFromNut_fromNut_fromModes
void computeProjSmagFromNut_fromNut_fromModes(volScalarField &phi, const volScalarField &Nut, const volVectorField &modeU_proj, const volVectorField &modeU_grad)
Compute projSmagFromNut given Nut and 2 velocity modes.
Definition UnsteadyNSTurb.C:2457
UnsteadyNSTurb::ct1PPEFluctTensor
Eigen::Tensor< double, 3 > ct1PPEFluctTensor
Turbulent fluctuation viscosity tensor for the splitting approach in the PPE equation.
Definition UnsteadyNSTurb.H:161
UnsteadyNSTurb::cTotalTensor
Eigen::Tensor< double, 3 > cTotalTensor
Turbulent total viscosity tensor.
Definition UnsteadyNSTurb.H:170
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:438
UnsteadyNSTurb::ct1Tensor
Eigen::Tensor< double, 3 > ct1Tensor
Turbulent viscosity tensor.
Definition UnsteadyNSTurb.H:134
UnsteadyNSTurb::cTotalPPETensor
Eigen::Tensor< double, 3 > cTotalPPETensor
Turbulent total viscosity tensor in the PPE equation.
Definition UnsteadyNSTurb.H:182
UnsteadyNSTurb::computeNut_fromU
volScalarField computeNut_fromU(const volVectorField &snapshotj)
Compute turbulent viscosity Snapshot for a given velocity field U.
Definition UnsteadyNSTurb.C:2521
UnsteadyNSTurb::nutAve
PtrList< volScalarField > nutAve
List of for eddy viscosity time-averaged fields.
Definition UnsteadyNSTurb.H:96
UnsteadyNSTurb::_pRefValue
scalar _pRefValue
Pressure reference value.
Definition UnsteadyNSTurb.H:221
UnsteadyNSTurb::btTurbulence
Eigen::MatrixXd btTurbulence(label NUmodes, label NSUPmodes)
bt added matrix for the turbulence treatement
Definition UnsteadyNSTurb.C:612
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:2150
UnsteadyNSTurb::ct2Tensor
Eigen::Tensor< double, 3 > ct2Tensor
Turbulent viscosity tensor.
Definition UnsteadyNSTurb.H:146
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:297
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:229
UnsteadyNSTurb::ct2PPEAveTensor
Eigen::Tensor< double, 3 > ct2PPEAveTensor
Turbulent average viscosity tensor for the splitting approach in the PPE equation.
Definition UnsteadyNSTurb.H:155
UnsteadyNSTurb::initSmagFunction
volVectorField initSmagFunction()
Initialize Smagorinsky term with values at time 0.
Definition UnsteadyNSTurb.C:2302
UnsteadyNSTurb::ct1PPETensor
Eigen::Tensor< double, 3 > ct1PPETensor
Turbulent viscosity tensor in the PPE equation.
Definition UnsteadyNSTurb.H:140
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:712
UnsteadyNSTurb::computeNut_fromS
volScalarField computeNut_fromS(const volTensorField &S)
Compute turbulent viscosity Snapshot for a given tensor Field S.
Definition UnsteadyNSTurb.C:2527
UnsteadyNSTurb::computeS_fromU
static volTensorField computeS_fromU(const volVectorField &snapshotj)
Compute S deformation tensor for a given vector field U.
Definition UnsteadyNSTurb.C:2598
UnsteadyNSTurb::cTotalAveTensor
Eigen::Tensor< double, 3 > cTotalAveTensor
Turbulent total average viscosity tensor for the splitting approach.
Definition UnsteadyNSTurb.H:173
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:332
UnsteadyNSTurb::turbulenceFluctTensor1
Eigen::Tensor< double, 3 > turbulenceFluctTensor1(label NUmodes, label NSUPmodes)
ct1Fluct added tensor for approximation of the fluctuating part of the eddy viscosity
Definition UnsteadyNSTurb.C:263
UnsteadyNSTurb::computeSmagTermPhi_at_time
volScalarField computeSmagTermPhi_at_time(const word &snap_time, const volScalarField template_field_phi)
Compute Smagorinsky Snapshot of a tracer phi at a given time.
Definition UnsteadyNSTurb.C:2329
UnsteadyNSTurb::computeSmagTermPhi_fromUPhi
volScalarField computeSmagTermPhi_fromUPhi(const volVectorField &snapshotj, const volScalarField &phij)
Compute Smagorinsky Snapshot for given velocity field U and tracer phi.
Definition UnsteadyNSTurb.C:2346
UnsteadyNSTurb::initNutFunction
volScalarField initNutFunction()
Initialize Nut term with values at time 0.
Definition UnsteadyNSTurb.C:2501
UnsteadyNSTurb::ct1PPEAveTensor
Eigen::Tensor< double, 3 > ct1PPEAveTensor
Turbulent average viscosity tensor for the splitting approach in the PPE equation.
Definition UnsteadyNSTurb.H:143
UnsteadyNSTurb::ct1FluctTensor
Eigen::Tensor< double, 3 > ct1FluctTensor
Turbulent fluctuation viscosity tensor for the splitting approach.
Definition UnsteadyNSTurb.H:158
UnsteadyNSTurb::turbulenceTensor2
Eigen::Tensor< double, 3 > turbulenceTensor2(label NUmodes, label NSUPmodes, label nNutModes)
ct2 added tensor for the turbulence treatement
Definition UnsteadyNSTurb.C:404
UnsteadyNSTurb::turbulencePPEFluctTensor2
Eigen::Tensor< double, 3 > turbulencePPEFluctTensor2(label NUmodes, label NSUPmodes, label NPmodes)
ct2PPEFluct added tensor for approximation of the fluctuating part of the eddy viscosity with the usa...
Definition UnsteadyNSTurb.C:577
UnsteadyNSTurb::_nut
autoPtr< volScalarField > _nut
Eddy viscosity field.
Definition UnsteadyNSTurb.H:212
UnsteadyNSTurb::interChoice
label interChoice
Interpolation independent variable choice.
Definition UnsteadyNSTurb.H:209
UnsteadyNSTurb::computeNut_at_time
volScalarField computeNut_at_time(const word &snap_time, std::optional< PtrList< volVectorField > > modesU=std::nullopt)
Compute non linear Snapshot at a given time.
Definition UnsteadyNSTurb.C:2506
UnsteadyNSTurb::btMatrix
Eigen::MatrixXd btMatrix
Turbulent viscosity term.
Definition UnsteadyNSTurb.H:122
UnsteadyNSTurb::nutModes
volScalarModes nutModes
List of POD modes for eddy viscosity.
Definition UnsteadyNSTurb.H:93
UnsteadyNSTurb::turbulenceFluctTensor2
Eigen::Tensor< double, 3 > turbulenceFluctTensor2(label NUmodes, label NSUPmodes)
ct2Fluct added tensor for approximation of the fluctuating part of the eddy viscosity
Definition UnsteadyNSTurb.C:472
UnsteadyNSTurb::initSmagPhiFunction
volScalarField initSmagPhiFunction(const volScalarField template_field_phi)
Initialize Smagorinsky term of a tracer phi with values at time 0.
Definition UnsteadyNSTurb.C:2308
UnsteadyNSTurb::projDiffusionIBP
static volScalarField projDiffusionIBP(const volScalarField &coefDiff, const volVectorField &u, const volVectorField &v)
Compute projected diffusion term with the integration by parts (IBP) equation for a given vector fiel...
Definition UnsteadyNSTurb.C:2463
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:541
UnsteadyNSTurb::ct2FluctTensor
Eigen::Tensor< double, 3 > ct2FluctTensor
Turbulent fluctuation viscosity tensor for the splitting approach.
Definition UnsteadyNSTurb.H:164
UnsteadyNSTurb::e
double e
RBF functions radius.
Definition UnsteadyNSTurb.H:203
UnsteadyNSTurb::ct1AveTensor
Eigen::Tensor< double, 3 > ct1AveTensor
Turbulent average viscosity tensor for the splitting approach.
Definition UnsteadyNSTurb.H:137
UnsteadyNSTurb::ct2PPEFluctTensor
Eigen::Tensor< double, 3 > ct2PPEFluctTensor
Turbulent fluctuation viscosity tensor for the splitting approach in the PPE equation.
Definition UnsteadyNSTurb.H:167
UnsteadyNSTurb::nutFields
PtrList< volScalarField > nutFields
List of snapshots for the solution for eddy viscosity.
Definition UnsteadyNSTurb.H:90
UnsteadyNSTurb::initProjSmagFromNutFunction
volScalarField initProjSmagFromNutFunction()
Initialize projSmagFromNut with values at time 0.
Definition UnsteadyNSTurb.C:2435
UnsteadyNSTurb::diffusion
static volScalarField diffusion(const T &coefDiff, const volScalarField &phi)
Compute diffusion term for a given scalar field phi.
Definition UnsteadyNSTurb.C:2354
UnsteadyNSTurb::_pRefCell
label _pRefCell
Pressure reference cell.
Definition UnsteadyNSTurb.H:218
UnsteadyNSTurb::computeSmagTerm_at_time
volVectorField computeSmagTerm_at_time(const word &snap_time, std::optional< PtrList< volVectorField > > modesU=std::nullopt)
Compute Smagorinsky Snapshot at a given time.
Definition UnsteadyNSTurb.C:2313
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:2136
UnsteadyNSTurb::initProjSmagFunction
volScalarField initProjSmagFunction()
Initialize projFullStressFunction with values at time 0.
Definition UnsteadyNSTurb.C:2390
UnsteadyNSTurb::bTotalMatrix
Eigen::MatrixXd bTotalMatrix
Total B Matrix.
Definition UnsteadyNSTurb.H:191
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:2100
UnsteadyNSTurb::computeProjSmagFromNut_at_time_fromModes
volScalarField computeProjSmagFromNut_at_time_fromModes(const word &snap_time, const volVectorField &modeU_proj, const volVectorField &modeU_grad)
Compute projSmagFromNut snapshot at a given time from 2 velocity modes.
Definition UnsteadyNSTurb.C:2441
UnsteadyNSTurb::nNutModes
label nNutModes
Number of viscoisty modes used for the projection.
Definition UnsteadyNSTurb.H:206
UnsteadyNSTurb::ct2AveTensor
Eigen::Tensor< double, 3 > ct2AveTensor
Turbulent average viscosity tensor for the splitting approach.
Definition UnsteadyNSTurb.H:149
UnsteadyNSTurb::computeProjSmagTerm_fromSmag_fromMode
void computeProjSmagTerm_fromSmag_fromMode(volScalarField &phi, const volVectorField &Smag, const volVectorField &mode)
Compute projFullStressFunction given the Smagorinsky term and the velocity mode.
Definition UnsteadyNSTurb.C:2413
UnsteadyNSTurb::turbulencePPEFluctTensor1
Eigen::Tensor< double, 3 > turbulencePPEFluctTensor1(label NUmodes, label NSUPmodes, label NPmodes)
ct1PPEFluct added tensor for approximation of the fluctuating part of the eddy viscosity with the usa...
Definition UnsteadyNSTurb.C:368
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:1269
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:506
UnsteadyNSTurb::z
Eigen::MatrixXd z
Time-parameter combined matrix.
Definition UnsteadyNSTurb.H:131
UnsteadyNSTurb::turbulenceTensor1
Eigen::Tensor< double, 3 > turbulenceTensor1(label NUmodes, label NSUPmodes, label nNutModes)
ct1 added tensor for the turbulence treatement
Definition UnsteadyNSTurb.C:195
UnsteadyNSTurb::ct2PPETensor
Eigen::Tensor< double, 3 > ct2PPETensor
Turbulent viscosity tensor in the PPE equation.
Definition UnsteadyNSTurb.H:152
UnsteadyNSTurb::computeProjSmagTerm_at_time_fromMode
volScalarField computeProjSmagTerm_at_time_fromMode(const word &snap_time, const volVectorField &mode)
Compute projFullStressFunction snapshot at a given time projected on the given mode.
Definition UnsteadyNSTurb.C:2396
UnsteadyNSTurb::cTotalPPEAveTensor
Eigen::Tensor< double, 3 > cTotalPPEAveTensor
Turbulent total average viscosity tensor for the splitting approach in the PPE equation.
Definition UnsteadyNSTurb.H:185
UnsteadyNSTurb::computeNonLinearSnapshot_at_time
void computeNonLinearSnapshot_at_time(const word &snap_time, volVectorField &Smag, std::optional< PtrList< volVectorField > > modesU=std::nullopt)
Compute non linear Snapshot at a given time.
Definition UnsteadyNSTurb.C:2284
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:227
steadyNS::maxIter
scalar maxIter
Number of maximum iterations to be done for the computation of the truth solution.
Definition steadyNS.H:128
steadyNS::NPmodes
label NPmodes
Number of pressure modes used for the projection.
Definition steadyNS.H:146
steadyNS::supex
bool supex
Boolean variable to check the existence of the supremizer modes.
Definition steadyNS.H:265
steadyNS::BC1_matrix
Eigen::MatrixXd BC1_matrix
PPE BC1.
Definition steadyNS.H:191
steadyNS::diffusive_term
Eigen::MatrixXd diffusive_term(label NUmodes, label NPmodes, label NSUPmodes)
Diffusive Term.
Definition steadyNS.C:883
steadyNS::_phi
autoPtr< surfaceScalarField > _phi
Flux.
Definition steadyNS.H:302
steadyNS::BC3_matrix
Eigen::MatrixXd BC3_matrix
PPE BC3.
Definition steadyNS.H:200
steadyNS::C_tensor
Eigen::Tensor< double, 3 > C_tensor
Diffusion term.
Definition steadyNS.H:170
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:1515
steadyNS::_fvOptions
autoPtr< fv::options > _fvOptions
fvOptions
Definition steadyNS.H:296
steadyNS::_runTime
autoPtr< Time > _runTime
Time.
Definition steadyNS.H:299
steadyNS::Umodes
volVectorModes Umodes
List of pointers used to form the velocity modes.
Definition steadyNS.H:101
steadyNS::Ufield
PtrList< volVectorField > Ufield
List of pointers used to form the velocity snapshots matrix.
Definition steadyNS.H:89
steadyNS::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:1319
steadyNS::tolerance
scalar tolerance
Tolerance for the residual of the stationary problems, there is the same tolerance for velocity and p...
Definition steadyNS.H:125
steadyNS::mass_term
Eigen::MatrixXd mass_term(label NUmodes, label NPmodes, label NSUPmodes)
Mass Term.
Definition steadyNS.C:1085
steadyNS::_mesh
autoPtr< fvMesh > _mesh
Mesh.
Definition steadyNS.H:290
steadyNS::_laminarTransport
autoPtr< singlePhaseTransportModel > _laminarTransport
Laminar transport (used by turbulence model).
Definition steadyNS.H:311
steadyNS::NUmodes
label NUmodes
Number of velocity modes used for the projection.
Definition steadyNS.H:143
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:320
steadyNS::pRefCell
label pRefCell
Reference pressure cell.
Definition steadyNS.H:317
steadyNS::L_U_SUPmodes
volVectorModes L_U_SUPmodes
List of pointers containing the total number of lift, supremizer and velocity modes.
Definition steadyNS.H:119
steadyNS::Uomfield
PtrList< volVectorField > Uomfield
List of pointers used to form the homogeneous velocity snapshots.
Definition steadyNS.H:116
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:1439
steadyNS::NNutModesOut
label NNutModesOut
Number of nut modes to be calculated.
Definition steadyNS.H:140
steadyNS::B_matrix
Eigen::MatrixXd B_matrix
Diffusion term.
Definition steadyNS.H:160
steadyNS::D_matrix
Eigen::MatrixXd D_matrix
Laplacian term PPE.
Definition steadyNS.H:182
steadyNS::_MRF
autoPtr< IOMRFZoneList > _MRF
MRF variable.
Definition steadyNS.H:314
steadyNS::NSUPmodes
label NSUPmodes
Number of supremizer modes used for the projection.
Definition steadyNS.H:149
steadyNS::gTensor
Eigen::Tensor< double, 3 > gTensor
Divergence of momentum PPE.
Definition steadyNS.H:188
steadyNS::NNutModes
label NNutModes
Number of nut modes used for the projection.
Definition steadyNS.H:152
steadyNS::K_matrix
Eigen::MatrixXd K_matrix
Gradient of pressure matrix.
Definition steadyNS.H:166
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:997
steadyNS::NUmodesOut
label NUmodesOut
Number of velocity modes to be calculated.
Definition steadyNS.H:131
steadyNS::divMomentum
Eigen::Tensor< double, 3 > divMomentum(label NUmodes, label NPmodes)
Divergence of convective term (PPE approach).
Definition steadyNS.C:1176
steadyNS::pressure_gradient_term
Eigen::MatrixXd pressure_gradient_term(label NUmodes, label NPmodes, label NSUPmodes)
Gradient of pressure.
Definition steadyNS.C:935
steadyNS::bcVelMat
List< Eigen::MatrixXd > bcVelMat
Boundary term for penalty method - matrix.
Definition steadyNS.H:230
steadyNS::M_matrix
Eigen::MatrixXd M_matrix
Mass Matrix.
Definition steadyNS.H:163
steadyNS::laplacian_pressure
Eigen::MatrixXd laplacian_pressure(label NPmodes)
Laplacian of pressure term (PPE approach).
Definition steadyNS.C:1295
steadyNS::NPmodesOut
label NPmodesOut
Number of pressure modes to be calculated.
Definition steadyNS.H:134
steadyNS::bcVelocityMat
List< Eigen::MatrixXd > bcVelocityMat(label NUmodes, label NSUPmodes)
Boundary integral modes on boundary used by the penaly method.
Definition steadyNS.C:1547
steadyNS::bc2Tensor
Eigen::Tensor< double, 3 > bc2Tensor
PPE BC2.
Definition steadyNS.H:197
steadyNS::_U
autoPtr< volVectorField > _U
Velocity field.
Definition steadyNS.H:272
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:1396
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:326
steadyNS::_p
autoPtr< volScalarField > _p
Pressure field.
Definition steadyNS.H:269
unsteadyNS::turbulence
autoPtr< incompressible::turbulenceModel > turbulence
Turbulence model.
Definition unsteadyNS.H:76
unsteadyNS::_pimple
autoPtr< pimpleControl > _pimple
pimpleControl
Definition unsteadyNS.H:73
ITHACAstream::ReadDenseMatrix
void ReadDenseMatrix(MatrixType &Matrix, word folder, word MatrixName)
Read a dense matrix from a binary format file.
Definition ITHACAstream.H:740
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::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:723
ITHACAstream::ReadDenseTensor
void ReadDenseTensor(TensorType &Tensor, word folder, word MatrixName)
Read a dense tensor from file.
Definition ITHACAstream.H:787
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:766
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::tensorFieldProduct
volTensorField tensorFieldProduct(const volScalarField &coef, const volTensorField &S)
Tensor field product between a volScalarField and a volTensorField.
Definition ITHACAfieldsOperations.C:152
ITHACAutilities::containsSubstring
bool containsSubstring(std::string contain, std::string contained)
Returns true if contained is a substring of contain, false otherwise.
Definition ITHACAstring.C:81
ITHACAutilities::project_to_POD_basis
T project_to_POD_basis(T &field, PtrList< T > &modes)
Compute the projection of a field in a POD basis.
Definition ITHACAcoeffsMass.C:471
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
Generated by   doxygen 1.16.1