Loading...
Searching...
No Matches
SteadyNSTurb.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
13License
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 "steadyNS.H"
32#include "SteadyNSTurb.H"
33#include "viscosityModel.H"
34
35// * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * * * //
36
37// Constructor
38SteadyNSTurb::SteadyNSTurb() {}
39
40SteadyNSTurb::SteadyNSTurb(int argc, char* argv[])
41{
42 _args = autoPtr<argList>
43 (
44 new argList(argc, argv)
45 );
46
47 if (!_args->checkRootCase())
48 {
49 Foam::FatalError.exit();
50 }
51
52 argList& args = _args();
53#include "createTime.H"
54#include "createMesh.H"
55 _simple = autoPtr<simpleControl>
56 (
57 new simpleControl
58 (
59 mesh
60 )
61 );
62 simpleControl& simple = _simple();
63#include "createFields.H"
64#include "createFvOptions.H"
65 //
66#pragma GCC diagnostic push
67#pragma GCC diagnostic ignored "-Wunused-variable"
68#include "initContinuityErrs.H"
69#pragma GCC diagnostic pop
70 //
71 turbulence->validate();
72 ITHACAdict = new IOdictionary
73 (
74 IOobject
75 (
76 "ITHACAdict",
77 runTime.system(),
78 mesh,
79 IOobject::MUST_READ,
80 IOobject::NO_WRITE
81 )
82 );
83 tolerance = ITHACAdict->lookupOrDefault<scalar>("tolerance", 1e-5);
84 maxIter = ITHACAdict->lookupOrDefault<scalar>("maxIter", 1000);
85 bcMethod = ITHACAdict->lookupOrDefault<word>("bcMethod", "lift");
86 M_Assert(bcMethod == "lift" || bcMethod == "penalty",
87 "The BC method must be set to lift or penalty in ITHACAdict");
88 viscDict = ITHACAdict->subDict("viscDict");
89 para = ITHACAparameters::getInstance(mesh, runTime);
90 offline = ITHACAutilities::check_off();
91 podex = ITHACAutilities::check_pod();
92 supex = ITHACAutilities::check_sup();
93}
94
95
96// * * * * * * * * * * * * * * Full Order Methods * * * * * * * * * * * * * * //
97
98// Method to performa a truthSolve
99void SteadyNSTurb::truthSolve(List<scalar> mu_now)
100{
101 Time& runTime = _runTime();
102 fvMesh& mesh = _mesh();
103 volScalarField& p = _p();
104 volVectorField& U = _U();
105 surfaceScalarField& phi = _phi();
106 fv::options& fvOptions = _fvOptions();
107 simpleControl& simple = _simple();
108 IOMRFZoneList& MRF = _MRF();
109 singlePhaseTransportModel& laminarTransport = _laminarTransport();
110#include "NLsolveSteadyNSTurb.H"
111 ITHACAstream::exportSolution(U, name(counter), "./ITHACAoutput/Offline/");
112 ITHACAstream::exportSolution(p, name(counter), "./ITHACAoutput/Offline/");
113 volScalarField _nut(turbulence->nut());
114 ITHACAstream::exportSolution(_nut, name(counter), "./ITHACAoutput/Offline/");
115 Ufield.append(U.clone());
116 Pfield.append(p.clone());
117 nutFields.append(_nut.clone());
118 counter++;
119 writeMu(mu_now);
120 // --- Fill in the mu_samples with parameters (mu) to be used for the PODI sample points
121 mu_samples.conservativeResize(mu_samples.rows() + 1, mu_now.size());
122
123 for (label i = 0; i < mu_now.size(); i++)
124 {
125 mu_samples(mu_samples.rows() - 1, i) = mu_now[i];
126 }
127
128 // Resize to Unitary if not initialized by user (i.e. non-parametric problem)
129 if (mu.cols() == 0)
130 {
131 mu.resize(1, 1);
132 }
133
134 if (mu_samples.rows() == mu.cols())
135 {
136 ITHACAstream::exportMatrix(mu_samples, "mu_samples", "eigen",
137 "./ITHACAoutput/Offline");
138 }
139}
140
141List < Eigen::MatrixXd > SteadyNSTurb::turbulenceTerm1(label NUmodes,
142 label NSUPmodes, label nNutModes)
143{
144 label cSize = NUmodes + NSUPmodes + liftfield.size();
145 List < Eigen::MatrixXd > ct1Matrix;
146 ct1Matrix.setSize(cSize);
147
148 for (label j = 0; j < cSize; j++)
149 {
150 ct1Matrix[j].resize(nNutModes, cSize);
151 ct1Matrix[j] = ct1Matrix[j] * 0;
152 }
153
154 for (label i = 0; i < cSize; i++)
155 {
156 Info << "Filling layer number " << i + 1 << " in the matrix ct1Matrix" << endl;
157
158 for (label j = 0; j < nNutModes; j++)
159 {
160 for (label k = 0; k < cSize; k++)
161 {
162 ct1Matrix[i](j, k) = fvc::domainIntegrate(L_U_SUPmodes[i] & fvc::laplacian(
163 nutModes[j], L_U_SUPmodes[k])).value();
164 }
165 }
166 }
167
168 // Export the matrix
169 ITHACAstream::exportMatrix(ct1Matrix, "ct1Matrix", "python",
170 "./ITHACAoutput/Matrices/");
171 ITHACAstream::exportMatrix(ct1Matrix, "ct1Matrix", "matlab",
172 "./ITHACAoutput/Matrices/");
173 ITHACAstream::exportMatrix(ct1Matrix, "ct1Matrix", "eigen",
174 "./ITHACAoutput/Matrices/ct1");
175 return ct1Matrix;
176}
177
178Eigen::Tensor<double, 3> SteadyNSTurb::turbulenceTensor1(label NUmodes,
179 label NSUPmodes, label nNutModes)
180{
181 label cSize = NUmodes + NSUPmodes + liftfield.size();
182 Eigen::Tensor<double, 3> ct1Tensor;
183 ct1Tensor.resize(cSize, nNutModes, cSize);
184
185 for (label i = 0; i < cSize; i++)
186 {
187 for (label j = 0; j < nNutModes; j++)
188 {
189 for (label k = 0; k < cSize; k++)
190 {
191 ct1Tensor(i, j, k) = fvc::domainIntegrate(L_U_SUPmodes[i] & fvc::laplacian(
192 nutModes[j], L_U_SUPmodes[k])).value();
193 }
194 }
195 }
196
197 // Export the tensor
198 ITHACAstream::SaveDenseTensor(ct1Tensor, "./ITHACAoutput/Matrices/",
199 "ct1_" + name(liftfield.size()) + "_" + name(NUmodes) + "_" + name(
200 NSUPmodes) + "_" + name(nNutModes) + "_t");
201 return ct1Tensor;
202}
203
204Eigen::Tensor<double, 3> SteadyNSTurb::turbulenceTensor1_cache(label NUmodes,
205 label NSUPmodes, label nNutModes)
206{
207 label cSize = NUmodes + NSUPmodes + liftfield.size();
208 Eigen::Tensor<double, 3> ct1Tensor(cSize, nNutModes, cSize);
209
210 for (label j = 0; j < nNutModes; ++j)
211 {
212 // Cache only one row (j fixed)
213 List<tmp<volVectorField>> lapRow(cSize);
214
215 for (label k = 0; k < cSize; ++k)
216 {
217 lapRow[k] = fvc::laplacian(nutModes[j], L_U_SUPmodes[k]);
218 }
219
220 for (label i = 0; i < cSize; ++i)
221 {
222 for (label k = 0; k < cSize; ++k)
223 {
224 const volVectorField& lapField = lapRow[k]();
225 ct1Tensor(i, j, k) = fvc::domainIntegrate(L_U_SUPmodes[i] & lapField).value();
226 }
227 }
228 }
229
230 // Export the tensor
231 if (Pstream::master())
232 {
233 ITHACAstream::SaveDenseTensor(ct1Tensor, "./ITHACAoutput/Matrices/",
234 "ct1_" + name(liftfield.size()) + "_" + name(NUmodes) + "_" + name(
235 NSUPmodes) + "_" + name(nNutModes) + "_t");
236 }
237
238 return ct1Tensor;
239}
240
241List < Eigen::MatrixXd > SteadyNSTurb::turbulenceTerm2(label NUmodes,
242 label NSUPmodes, label nNutModes)
243{
244 label cSize = NUmodes + NSUPmodes + liftfield.size();
245 List < Eigen::MatrixXd > ct2Matrix;
246 ct2Matrix.setSize(cSize);
247
248 for (label j = 0; j < cSize; j++)
249 {
250 ct2Matrix[j].resize(nNutModes, cSize);
251 ct2Matrix[j] = ct2Matrix[j] * 0;
252 }
253
254 for (label i = 0; i < cSize; i++)
255 {
256 Info << "Filling layer number " << i + 1 << " in the matrix ct2Matrix" << endl;
257
258 for (label j = 0; j < nNutModes; j++)
259 {
260 for (label k = 0; k < cSize; k++)
261 {
262 ct2Matrix[i](j, k) = fvc::domainIntegrate(L_U_SUPmodes[i] & (fvc::div(
263 nutModes[j] * dev((fvc::grad(L_U_SUPmodes[k]))().T())))).value();
264 }
265 }
266 }
267
268 // Export the matrix
269 ITHACAstream::exportMatrix(ct2Matrix, "ct2Matrix", "python",
270 "./ITHACAoutput/Matrices/");
271 ITHACAstream::exportMatrix(ct2Matrix, "ct2Matrix", "matlab",
272 "./ITHACAoutput/Matrices/");
273 ITHACAstream::exportMatrix(ct2Matrix, "ct2Matrix", "eigen",
274 "./ITHACAoutput/Matrices/ct2");
275 return ct2Matrix;
276}
277
278Eigen::Tensor<double, 3> SteadyNSTurb::turbulenceTensor2(label NUmodes,
279 label NSUPmodes, label nNutModes)
280{
281 label cSize = NUmodes + NSUPmodes + liftfield.size();
282 Eigen::Tensor<double, 3> ct2Tensor;
283 ct2Tensor.resize(cSize, nNutModes, cSize);
284
285 for (label i = 0; i < cSize; i++)
286 {
287 for (label j = 0; j < nNutModes; j++)
288 {
289 for (label k = 0; k < cSize; k++)
290 {
291 ct2Tensor(i, j, k) = fvc::domainIntegrate(L_U_SUPmodes[i] & (fvc::div(
292 nutModes[j] * dev((fvc::grad(L_U_SUPmodes[k]))().T())))).value();
293 }
294 }
295 }
296
297 // Export the tensor
298 ITHACAstream::SaveDenseTensor(ct2Tensor, "./ITHACAoutput/Matrices/",
299 "ct2_" + name(liftfield.size()) + "_" + name(NUmodes) + "_" + name(
300 NSUPmodes) + "_" + name(nNutModes) + "_t");
301 return ct2Tensor;
302}
303
304Eigen::Tensor<double, 3> SteadyNSTurb::turbulenceTensor2_cache(label NUmodes,
305 label NSUPmodes, label nNutModes)
306{
307 label cSize = NUmodes + NSUPmodes + liftfield.size();
308 Eigen::Tensor<double, 3> ct2Tensor(cSize, nNutModes, cSize);
309
310 for (label j = 0; j < nNutModes; ++j)
311 {
312 // Cache only one j-row
313 List<tmp<volVectorField>> divRow(cSize);
314
315 for (label k = 0; k < cSize; ++k)
316 {
317 divRow[k] = fvc::div(nutModes[j] * dev((fvc::grad(L_U_SUPmodes[k]))().T()));
318 }
319
320 for (label i = 0; i < cSize; ++i)
321 {
322 for (label k = 0; k < cSize; ++k)
323 {
324 const volVectorField& divRowField = divRow[k]();
325 ct2Tensor(i, j, k) = fvc::domainIntegrate(L_U_SUPmodes[i] &
326 divRowField).value();
327 }
328 }
329 }
330
331 // Export the tensor
332 if (Pstream::master())
333 {
334 ITHACAstream::SaveDenseTensor(ct2Tensor, "./ITHACAoutput/Matrices/",
335 "ct2_" + name(liftfield.size()) + "_" + name(NUmodes) + "_" + name(
336 NSUPmodes) + "_" + name(nNutModes) + "_t");
337 }
338
339 return ct2Tensor;
340}
341
342Eigen::Tensor<double, 3> SteadyNSTurb::turbulencePPETensor1(label NUmodes,
343 label NSUPmodes, label NPmodes, label nNutModes)
344{
345 label cSize = NUmodes + NSUPmodes + liftfield.size();
346 Eigen::Tensor<double, 3> ct1PPETensor;
347 ct1PPETensor.resize(NPmodes, nNutModes, cSize);
348
349 for (label i = 0; i < NPmodes; i++)
350 {
351 for (label j = 0; j < nNutModes; j++)
352 {
353 for (label k = 0; k < cSize; k++)
354 {
355 // ct1PPETensor(i, j, k) = fvc::domainIntegrate(2 * Pmodes[i] * (fvc::laplacian(
356 // L_U_SUPmodes[k]) & fvc::grad(nutModes[j]))).value();
357 // ct1PPETensor(i, j, k) = fvc::domainIntegrate(Pmodes[i] * (fvc::div(
358 // fvc::laplacian(
359 // nutModes[j], L_U_SUPmodes[k])))).value();
360 ct1PPETensor(i, j, k) = fvc::domainIntegrate(fvc::grad(Pmodes[i]) & (
361 fvc::laplacian(
362 nutModes[j], L_U_SUPmodes[k]))).value();
363 }
364 }
365 }
366
367 // Export the tensor
368 ITHACAstream::SaveDenseTensor(ct1PPETensor, "./ITHACAoutput/Matrices/",
369 "ct1PPE_" + name(liftfield.size()) + "_" + name(NUmodes) + "_" + name(
370 NSUPmodes) + "_" + name(NPmodes) + "_" + name(nNutModes) + "_t");
371 return ct1PPETensor;
372}
373
374Eigen::Tensor<double, 3> SteadyNSTurb::turbulencePPETensor1_cache(label NUmodes,
375 label NSUPmodes, label NPmodes, label nNutModes)
376{
377 label cSize = NUmodes + NSUPmodes + liftfield.size();
378 Eigen::Tensor<double, 3> ct1PPETensor(NPmodes, nNutModes, cSize);
379 List<tmp<volVectorField>> PmodesGrad(NPmodes);
380
381 for (label i = 0; i < NPmodes; ++i)
382 {
383 PmodesGrad[i] = fvc::grad(Pmodes[i]);
384 }
385
386 for (label j = 0; j < nNutModes; ++j)
387 {
388 // Cache one row (j fixed)
389 List<tmp<volVectorField>> lapRow(cSize);
390
391 for (label k = 0; k < cSize; ++k)
392 {
393 lapRow[k] = fvc::laplacian(nutModes[j], L_U_SUPmodes[k]);
394 }
395
396 for (label i = 0; i < NPmodes; ++i)
397 {
398 const volVectorField& gradPi = PmodesGrad[i]();
399
400 for (label k = 0; k < cSize; ++k)
401 {
402 const volVectorField& lapRowField = lapRow[k]();
403 ct1PPETensor(i, j, k) = fvc::domainIntegrate(gradPi & lapRowField).value();
404 }
405 }
406 }
407
408 // Export the tensor
409 if (Pstream::master())
410 {
411 ITHACAstream::SaveDenseTensor(ct1PPETensor, "./ITHACAoutput/Matrices/",
412 "ct1PPE_" + name(liftfield.size()) + "_" + name(NUmodes) + "_" + name(
413 NSUPmodes) + "_" + name(NPmodes) + "_" + name(nNutModes) + "_t");
414 }
415
416 return ct1PPETensor;
417}
418
419Eigen::Tensor<double, 3> SteadyNSTurb::turbulencePPETensor2(label NUmodes,
420 label NSUPmodes, label NPmodes, label nNutModes)
421{
422 label cSize = NUmodes + NSUPmodes + liftfield.size();
423 Eigen::Tensor<double, 3> ct2PPETensor;
424 ct2PPETensor.resize(NPmodes, nNutModes, cSize);
425
426 for (label i = 0; i < NPmodes; i++)
427 {
428 for (label j = 0; j < nNutModes; j++)
429 {
430 for (label k = 0; k < cSize; k++)
431 {
432 // ct2PPETensor(i, j, k) = fvc::domainIntegrate(Pmodes[i] * (fvc::grad(fvc::grad(
433 // nutModes[j])) && (dev2((fvc::grad(L_U_SUPmodes[k]))() + fvc::grad(
434 // L_U_SUPmodes[k]))().T()))).value();
435 // ct2PPETensor(i, j, k) = fvc::domainIntegrate(Pmodes[i] * ((fvc::div(fvc::div(
436 // nutModes[j] * dev2((fvc::grad(L_U_SUPmodes[k]))().T())))))).value();
437 ct2PPETensor(i, j, k) = fvc::domainIntegrate(fvc::grad(Pmodes[i]) & ((fvc::div(
438 nutModes[j] * dev2((fvc::grad(L_U_SUPmodes[k]))().T()))))).value();
439 }
440 }
441 }
442
443 // Export the tensor
444 ITHACAstream::SaveDenseTensor(ct2PPETensor, "./ITHACAoutput/Matrices/",
445 "ct2PPE_" + name(liftfield.size()) + "_" + name(NUmodes) + "_" + name(
446 NSUPmodes) + "_" + name(NPmodes) + "_" + name(nNutModes) + "_t");
447 return ct2PPETensor;
448}
449
450Eigen::Tensor<double, 3> SteadyNSTurb::turbulencePPETensor2_cache(label NUmodes,
451 label NSUPmodes, label NPmodes, label nNutModes)
452{
453 label cSize = NUmodes + NSUPmodes + liftfield.size();
454 Eigen::Tensor<double, 3> ct2PPETensor(NPmodes, nNutModes, cSize);
455 List<tmp<volVectorField>> PmodesGrad(NPmodes);
456
457 for (label i = 0; i < NPmodes; ++i)
458 {
459 PmodesGrad[i] = fvc::grad(Pmodes[i]);
460 }
461
462 for (label j = 0; j < nNutModes; ++j)
463 {
464 // Cache one row of div fields for this nutMode[j]
465 List<tmp<volVectorField>> divLow(cSize);
466
467 for (label k = 0; k < cSize; ++k)
468 {
469 divLow[k] = fvc::div(nutModes[j] * dev2((fvc::grad(L_U_SUPmodes[k]))().T()));
470 }
471
472 for (label i = 0; i < NPmodes; ++i)
473 {
474 const volVectorField& gradPi = PmodesGrad[i]();
475
476 for (label k = 0; k < cSize; ++k)
477 {
478 const volVectorField& divLowField = divLow[k]();
479 ct2PPETensor(i, j, k) = fvc::domainIntegrate(gradPi & divLowField).value();
480 }
481 }
482 }
483
484 // Export the tensor
485 if (Pstream::master())
486 {
487 ITHACAstream::SaveDenseTensor(ct2PPETensor, "./ITHACAoutput/Matrices/",
488 "ct2PPE_" + name(liftfield.size()) + "_" + name(NUmodes) + "_" + name(
489 NSUPmodes) + "_" + name(NPmodes) + "_" + name(nNutModes) + "_t");
490 }
491
492 return ct2PPETensor;
493}
494
495Eigen::MatrixXd SteadyNSTurb::btTurbulence(label NUmodes, label NSUPmodes)
496{
497 label btSize = NUmodes + NSUPmodes + liftfield.size();
498 Eigen::MatrixXd btMatrix(btSize, btSize);
499 btMatrix = btMatrix * 0;
500
501 // Project everything
502 for (label i = 0; i < btSize; i++)
503 {
504 for (label j = 0; j < btSize; j++)
505 {
506 btMatrix(i, j) = fvc::domainIntegrate(L_U_SUPmodes[i] & (fvc::div(dev((T(
507 fvc::grad(
508 L_U_SUPmodes[j]))))))).value();
509 }
510 }
511
512 // Export the matrix
513 ITHACAstream::SaveDenseMatrix(btMatrix, "./ITHACAoutput/Matrices/",
514 "bt_" + name(liftfield.size()) + "_" + name(NUmodes) + "_" + name(NSUPmodes));
515 return btMatrix;
516}
517
518void SteadyNSTurb::projectPPE(fileName folder, label NU, label NP, label NSUP,
519 label Nnut)
520{
521 NUmodes = NU;
522 NPmodes = NP;
523 NSUPmodes = 0;
524 nNutModes = Nnut;
525 L_U_SUPmodes.resize(0);
526
527 if (liftfield.size() != 0)
528 {
529 for (label k = 0; k < liftfield.size(); k++)
530 {
531 L_U_SUPmodes.append(liftfield[k].clone());
532 }
533 }
534
535 if (NUmodes != 0)
536 {
537 for (label k = 0; k < NUmodes; k++)
538 {
539 L_U_SUPmodes.append(Umodes[k].clone());
540 }
541 }
542
543 if (ITHACAutilities::check_folder("./ITHACAoutput/Matrices/"))
544 {
545 word bStr = "B_" + name(liftfield.size()) + "_" + name(NUmodes) + "_" + name(
546 NSUPmodes);
547
548 if (ITHACAutilities::check_file("./ITHACAoutput/Matrices/" + bStr))
549 {
550 ITHACAstream::ReadDenseMatrix(B_matrix, "./ITHACAoutput/Matrices/", bStr);
551 }
552 else
553 {
554 B_matrix = diffusive_term(NUmodes, NPmodes, NSUPmodes);
555 }
556
557 word btStr = "bt_" + name(liftfield.size()) + "_" + name(NUmodes) + "_" + name(
558 NSUPmodes);
559
560 if (ITHACAutilities::check_file("./ITHACAoutput/Matrices/" + btStr))
561 {
562 ITHACAstream::ReadDenseMatrix(btMatrix, "./ITHACAoutput/Matrices/", btStr);
563 }
564 else
565 {
566 btMatrix = btTurbulence(NUmodes, NSUPmodes);
567 }
568
569 word kStr = "K_" + name(liftfield.size()) + "_" + name(NUmodes) + "_" + name(
570 NSUPmodes) + "_" + name(NPmodes);
571
572 if (ITHACAutilities::check_file("./ITHACAoutput/Matrices/" + kStr))
573 {
574 ITHACAstream::ReadDenseMatrix(K_matrix, "./ITHACAoutput/Matrices/", kStr);
575 }
576 else
577 {
578 K_matrix = pressure_gradient_term(NUmodes, NPmodes, NSUPmodes);
579 }
580
581 word pStr = "P_" + name(liftfield.size()) + "_" + name(NUmodes) + "_" + name(
582 NSUPmodes) + "_" + name(NPmodes);
583
584 if (ITHACAutilities::check_file("./ITHACAoutput/Matrices/" + pStr))
585 {
586 ITHACAstream::ReadDenseMatrix(P_matrix, "./ITHACAoutput/Matrices/", pStr);
587 }
588 else
589 {
590 P_matrix = divergence_term(NUmodes, NPmodes, NSUPmodes);
591 }
592
593 word mStr = "M_" + name(liftfield.size()) + "_" + name(NUmodes) + "_" + name(
594 NSUPmodes);
595
596 if (ITHACAutilities::check_file("./ITHACAoutput/Matrices/" + mStr))
597 {
598 ITHACAstream::ReadDenseMatrix(M_matrix, "./ITHACAoutput/Matrices/", mStr);
599 }
600 else
601 {
602 M_matrix = mass_term(NUmodes, NPmodes, NSUPmodes);
603 }
604
605 word D_str = "D_" + name(NPmodes);
606
607 if (ITHACAutilities::check_file("./ITHACAoutput/Matrices/" + D_str))
608 {
609 ITHACAstream::ReadDenseMatrix(D_matrix, "./ITHACAoutput/Matrices/", D_str);
610 }
611 else
612 {
613 D_matrix = laplacian_pressure(NPmodes);
614 }
615
616 word bc3_str = "BC3_" + name(liftfield.size()) + "_" + name(
617 NUmodes) + "_" + name(NPmodes);
618
619 if (ITHACAutilities::check_file("./ITHACAoutput/Matrices/" + bc3_str))
620 {
621 ITHACAstream::ReadDenseMatrix(BC3_matrix, "./ITHACAoutput/Matrices/", bc3_str);
622 }
623 else
624 {
625 BC3_matrix = pressure_BC3(NUmodes, NPmodes);
626 }
627
628 word bc4_str = "BC4_" + name(liftfield.size()) + "_" + name(
629 NUmodes) + "_" + name(NPmodes);
630
631 if (ITHACAutilities::check_file("./ITHACAoutput/Matrices/" + bc4_str))
632 {
633 ITHACAstream::ReadDenseMatrix(BC4_matrix, "./ITHACAoutput/Matrices/", bc4_str);
634 }
635 else
636 {
637 BC4_matrix = pressure_BC4(NUmodes, NPmodes);
638 }
639
640 word C_str = "C_" + name(liftfield.size()) + "_" + name(NUmodes) + "_" + name(
641 NSUPmodes) + "_t";
642
643 if (ITHACAutilities::check_file("./ITHACAoutput/Matrices/" + C_str))
644 {
645 ITHACAstream::ReadDenseTensor(C_tensor, "./ITHACAoutput/Matrices/", C_str);
646 }
647 else
648 {
649 C_tensor = convective_term_tens(NUmodes, NPmodes, NSUPmodes);
650 }
651
652 word ct1Str = "ct1_" + name(liftfield.size()) + "_" + name(
653 NUmodes) + "_" + name(
654 NSUPmodes) + "_" + name(nNutModes) + "_t";
655
656 if (ITHACAutilities::check_file("./ITHACAoutput/Matrices/" + ct1Str))
657 {
658 ITHACAstream::ReadDenseTensor(ct1Tensor, "./ITHACAoutput/Matrices/", ct1Str);
659 }
660 else
661 {
662 ct1Tensor = turbulenceTensor1(NUmodes, NSUPmodes, nNutModes);
663 }
664
665 word ct2Str = "ct2_" + name(liftfield.size()) + "_" + name(
666 NUmodes) + "_" + name(
667 NSUPmodes) + "_" + name(nNutModes) + "_t";
668
669 if (ITHACAutilities::check_file("./ITHACAoutput/Matrices/" + ct2Str))
670 {
671 ITHACAstream::ReadDenseTensor(ct2Tensor, "./ITHACAoutput/Matrices/", ct2Str);
672 }
673 else
674 {
675 ct2Tensor = turbulenceTensor2(NUmodes, NSUPmodes, nNutModes);
676 }
677
678 word G_str = "G_" + name(liftfield.size()) + "_" + name(NUmodes) + "_" + name(
679 NSUPmodes) + "_" + name(NPmodes) + "_t";
680
681 if (ITHACAutilities::check_file("./ITHACAoutput/Matrices/" + G_str))
682 {
683 ITHACAstream::ReadDenseTensor(gTensor, "./ITHACAoutput/Matrices/", G_str);
684 }
685 else
686 {
687 gTensor = divMomentum(NUmodes, NPmodes);
688 }
689
690 if (bcMethod == "penalty")
691 {
692 bcVelVec = bcVelocityVec(NUmodes, NSUPmodes);
693 bcVelMat = bcVelocityMat(NUmodes, NSUPmodes);
694 }
695 }
696 else
697 {
698 L_U_SUPmodes.resize(0);
699
700 if (liftfield.size() != 0)
701 {
702 for (label k = 0; k < liftfield.size(); k++)
703 {
704 L_U_SUPmodes.append(liftfield[k].clone());
705 }
706 }
707
708 if (NUmodes != 0)
709 {
710 for (label k = 0; k < NUmodes; k++)
711 {
712 L_U_SUPmodes.append(Umodes[k].clone());
713 }
714 }
715
716 if (NSUPmodes != 0)
717 {
718 for (label k = 0; k < NSUPmodes; k++)
719 {
720 L_U_SUPmodes.append(supmodes[k].clone());
721 }
722 }
723
724 B_matrix = diffusive_term(NUmodes, NPmodes, NSUPmodes);
725 C_tensor = convective_term_tens(NUmodes, NPmodes, NSUPmodes);
726 K_matrix = pressure_gradient_term(NUmodes, NPmodes, NSUPmodes);
727 M_matrix = mass_term(NUmodes, NPmodes, NSUPmodes);
728 D_matrix = laplacian_pressure(NPmodes);
729 gTensor = divMomentum(NUmodes, NPmodes);
730 btMatrix = btTurbulence(NUmodes, NSUPmodes);
731 BC3_matrix = pressure_BC3(NUmodes, NPmodes);
732 // BC4_matrix = pressure_BC4(NUmodes, NPmodes);
733 ct1Tensor = turbulenceTensor1(NUmodes, NSUPmodes, nNutModes);
734 ct2Tensor = turbulenceTensor2(NUmodes, NSUPmodes, nNutModes);
735 ct1PPETensor = turbulencePPETensor1(NUmodes, NSUPmodes, NPmodes, nNutModes);
736 ct2PPETensor = turbulencePPETensor2(NUmodes, NSUPmodes, NPmodes, nNutModes);
737
738 if (bcMethod == "penalty")
739 {
740 bcVelVec = bcVelocityVec(NUmodes, NSUPmodes);
741 bcVelMat = bcVelocityMat(NUmodes, NSUPmodes);
742 }
743 }
744
745 // Export the matrices
746 if (para->exportPython)
747 {
748 ITHACAstream::exportMatrix(B_matrix, "B", "python", "./ITHACAoutput/Matrices/");
749 ITHACAstream::exportMatrix(K_matrix, "K", "python", "./ITHACAoutput/Matrices/");
750 ITHACAstream::exportMatrix(D_matrix, "D", "python", "./ITHACAoutput/Matrices/");
751 ITHACAstream::exportMatrix(M_matrix, "M", "python", "./ITHACAoutput/Matrices/");
752 ITHACAstream::exportMatrix(BC3_matrix, "BC3", "python",
753 "./ITHACAoutput/Matrices/");
754 ITHACAstream::exportMatrix(BC4_matrix, "BC4", "python",
755 "./ITHACAoutput/Matrices/");
756 ITHACAstream::exportMatrix(btMatrix, "bt", "python",
757 "./ITHACAoutput/Matrices/");
758 ITHACAstream::exportTensor(C_tensor, "C", "python", "./ITHACAoutput/Matrices/");
759 ITHACAstream::exportTensor(gTensor, "G", "python", "./ITHACAoutput/Matrices/");
760 ITHACAstream::exportTensor(ct1Tensor, "ct1", "python",
761 "./ITHACAoutput/Matrices/");
762 ITHACAstream::exportTensor(ct2Tensor, "ct2", "python",
763 "./ITHACAoutput/Matrices/");
764 ITHACAstream::exportTensor(ct1PPETensor, "ct1PPE", "python",
765 "./ITHACAoutput/Matrices/");
766 ITHACAstream::exportTensor(ct2PPETensor, "ct2PPE", "python",
767 "./ITHACAoutput/Matrices/");
768 }
769
770 if (para->exportMatlab)
771 {
772 ITHACAstream::exportMatrix(B_matrix, "B", "matlab", "./ITHACAoutput/Matrices/");
773 ITHACAstream::exportMatrix(K_matrix, "K", "matlab", "./ITHACAoutput/Matrices/");
774 ITHACAstream::exportMatrix(D_matrix, "D", "matlab", "./ITHACAoutput/Matrices/");
775 ITHACAstream::exportMatrix(M_matrix, "M", "matlab", "./ITHACAoutput/Matrices/");
776 ITHACAstream::exportMatrix(BC3_matrix, "BC3", "matlab",
777 "./ITHACAoutput/Matrices/");
778 ITHACAstream::exportMatrix(BC4_matrix, "BC4", "matlab",
779 "./ITHACAoutput/Matrices/");
780 ITHACAstream::exportMatrix(btMatrix, "bt", "matlab",
781 "./ITHACAoutput/Matrices/");
782 ITHACAstream::exportTensor(gTensor, "G", "matlab", "./ITHACAoutput/Matrices/");
783 ITHACAstream::exportTensor(ct1Tensor, "ct1", "matlab",
784 "./ITHACAoutput/Matrices/");
785 ITHACAstream::exportTensor(ct2Tensor, "ct2", "matlab",
786 "./ITHACAoutput/Matrices/");
787 ITHACAstream::exportTensor(ct1PPETensor, "ct1PPE", "matlab",
788 "./ITHACAoutput/Matrices/");
789 ITHACAstream::exportTensor(ct2PPETensor, "ct2PPE", "matlab",
790 "./ITHACAoutput/Matrices/");
791 }
792
793 if (para->exportTxt)
794 {
795 ITHACAstream::exportMatrix(B_matrix, "B", "eigen", "./ITHACAoutput/Matrices/");
796 ITHACAstream::exportMatrix(K_matrix, "K", "eigen", "./ITHACAoutput/Matrices/");
797 ITHACAstream::exportMatrix(D_matrix, "D", "eigen", "./ITHACAoutput/Matrices/");
798 ITHACAstream::exportMatrix(M_matrix, "M", "eigen", "./ITHACAoutput/Matrices/");
799 ITHACAstream::exportMatrix(BC3_matrix, "BC3", "eigen",
800 "./ITHACAoutput/Matrices/");
801 ITHACAstream::exportMatrix(BC4_matrix, "BC4", "eigen",
802 "./ITHACAoutput/Matrices/");
803 ITHACAstream::exportMatrix(btMatrix, "bt", "eigen", "./ITHACAoutput/Matrices/");
804 ITHACAstream::exportTensor(gTensor, "G", "eigen", "./ITHACAoutput/Matrices/G");
805 ITHACAstream::exportTensor(C_tensor, "C", "eigen", "./ITHACAoutput/Matrices/C");
806 ITHACAstream::exportTensor(ct1Tensor, "ct1_", "eigen",
807 "./ITHACAoutput/Matrices/ct1");
808 ITHACAstream::exportTensor(ct2Tensor, "ct2_", "eigen",
809 "./ITHACAoutput/Matrices/ct2");
810 ITHACAstream::exportTensor(ct1PPETensor, "ct1PPE_", "eigen",
811 "./ITHACAoutput/Matrices/ct1PPE");
812 ITHACAstream::exportTensor(ct2PPETensor, "ct2PPE_", "eigen",
813 "./ITHACAoutput/Matrices/ct2PPE");
814 }
815
816 bTotalMatrix = B_matrix + btMatrix;
817 label cSize = NUmodes + NSUPmodes + liftfield.size();
818 cTotalTensor.resize(cSize, nNutModes, cSize);
819 cTotalTensor = ct1Tensor + ct2Tensor;
820 cTotalPPETensor.resize(NPmodes, nNutModes, cSize);
821 cTotalPPETensor = ct1PPETensor + ct2PPETensor;
822 // Get the coeffs for interpolation (the orthonormal one is used because basis are orthogonal)
823 coeffL2 = ITHACAutilities::getCoeffs(nutFields,
824 nutModes, nNutModes);
825 ITHACAstream::exportMatrix(coeffL2, "coeffL2", "python",
826 "./ITHACAoutput/Matrices/");
827 ITHACAstream::exportMatrix(coeffL2, "coeffL2", "matlab",
828 "./ITHACAoutput/Matrices/");
829 ITHACAstream::exportMatrix(coeffL2, "coeffL2", "eigen",
830 "./ITHACAoutput/Matrices/");
831 // Export the matrix
832 ITHACAstream::SaveDenseMatrix(coeffL2, "./ITHACAoutput/Matrices/",
833 "coeffL2_nut_" + name(nNutModes));
834
835 // Create RBF interpolators for nut coefficient interpolation
836 rbfSplines.resize(nNutModes);
837 for (label i = 0; i < nNutModes; i++)
838 {
839 // Create ithacaInterpolator instance
840 rbfSplines[i] = std::make_shared<ithacaInterpolator>(viscDict);
841
842 // Prepare training data: X is parameter matrix (transposed), y is coefficient vector
843 Eigen::MatrixXd X = mu.transpose(); // Now each row is a parameter sample
844 Eigen::VectorXd y = coeffL2.row(i).transpose(); // Coefficient vector for this mode
845
846 rbfSplines[i]->fit(X, y);
847
848 Info << "Constructing ithacaInterpolator for mode " << i + 1 << endl;
849 }
850
851 Info<< "Info on interpolators for nut coefficients: "<< endl;
852 rbfSplines[0]->printInfo();
853}
854
855void SteadyNSTurb::projectSUP(fileName folder, label NU, label NP, label NSUP,
856 label Nnut)
857{
858 NUmodes = NU;
859 NPmodes = NP;
860 NSUPmodes = NSUP;
861 nNutModes = Nnut;
862 L_U_SUPmodes.resize(0);
863
864 if (liftfield.size() != 0)
865 {
866 for (label k = 0; k < liftfield.size(); k++)
867 {
868 L_U_SUPmodes.append(liftfield[k].clone());
869 }
870 }
871
872 if (NUmodes != 0)
873 {
874 for (label k = 0; k < NUmodes; k++)
875 {
876 L_U_SUPmodes.append(Umodes[k].clone());
877 }
878 }
879
880 if (NSUPmodes != 0)
881 {
882 for (label k = 0; k < NSUPmodes; k++)
883 {
884 L_U_SUPmodes.append(supmodes[k].clone());
885 }
886 }
887
888 if (ITHACAutilities::check_folder("./ITHACAoutput/Matrices/"))
889 {
890 word bStr = "B_" + name(liftfield.size()) + "_" + name(NUmodes) + "_" + name(
891 NSUPmodes);
892
893 if (ITHACAutilities::check_file("./ITHACAoutput/Matrices/" + bStr))
894 {
895 ITHACAstream::ReadDenseMatrix(B_matrix, "./ITHACAoutput/Matrices/", bStr);
896 }
897 else
898 {
899 B_matrix = diffusive_term(NUmodes, NPmodes, NSUPmodes);
900 }
901
902 word btStr = "bt_" + name(liftfield.size()) + "_" + name(NUmodes) + "_" + name(
903 NSUPmodes);
904
905 if (ITHACAutilities::check_file("./ITHACAoutput/Matrices/" + btStr))
906 {
907 ITHACAstream::ReadDenseMatrix(btMatrix, "./ITHACAoutput/Matrices/", btStr);
908 }
909 else
910 {
911 btMatrix = btTurbulence(NUmodes, NSUPmodes);
912 }
913
914 word kStr = "K_" + name(liftfield.size()) + "_" + name(NUmodes) + "_" + name(
915 NSUPmodes) + "_" + name(NPmodes);
916
917 if (ITHACAutilities::check_file("./ITHACAoutput/Matrices/" + kStr))
918 {
919 ITHACAstream::ReadDenseMatrix(K_matrix, "./ITHACAoutput/Matrices/", kStr);
920 }
921 else
922 {
923 K_matrix = pressure_gradient_term(NUmodes, NPmodes, NSUPmodes);
924 }
925
926 word pStr = "P_" + name(liftfield.size()) + "_" + name(NUmodes) + "_" + name(
927 NSUPmodes) + "_" + name(NPmodes);
928
929 if (ITHACAutilities::check_file("./ITHACAoutput/Matrices/" + pStr))
930 {
931 ITHACAstream::ReadDenseMatrix(P_matrix, "./ITHACAoutput/Matrices/", pStr);
932 }
933 else
934 {
935 P_matrix = divergence_term(NUmodes, NPmodes, NSUPmodes);
936 }
937
938 word mStr = "M_" + name(liftfield.size()) + "_" + name(NUmodes) + "_" + name(
939 NSUPmodes);
940
941 if (ITHACAutilities::check_file("./ITHACAoutput/Matrices/" + mStr))
942 {
943 ITHACAstream::ReadDenseMatrix(M_matrix, "./ITHACAoutput/Matrices/", mStr);
944 }
945 else
946 {
947 M_matrix = mass_term(NUmodes, NPmodes, NSUPmodes);
948 }
949
950 word C_str = "C_" + name(liftfield.size()) + "_" + name(NUmodes) + "_" + name(
951 NSUPmodes) + "_t";
952
953 if (ITHACAutilities::check_file("./ITHACAoutput/Matrices/" + C_str))
954 {
955 ITHACAstream::ReadDenseTensor(C_tensor, "./ITHACAoutput/Matrices/", C_str);
956 }
957 else
958 {
959 C_tensor = convective_term_tens(NUmodes, NPmodes, NSUPmodes);
960 }
961
962 word ct1Str = "ct1_" + name(liftfield.size()) + "_" + name(
963 NUmodes) + "_" + name(
964 NSUPmodes) + "_" + name(nNutModes) + "_t";
965
966 if (ITHACAutilities::check_file("./ITHACAoutput/Matrices/" + ct1Str))
967 {
968 ITHACAstream::ReadDenseTensor(ct1Tensor, "./ITHACAoutput/Matrices/", ct1Str);
969 }
970 else
971 {
972 ct1Tensor = turbulenceTensor1(NUmodes, NSUPmodes, nNutModes);
973 }
974
975 word ct2Str = "ct2_" + name(liftfield.size()) + "_" + name(
976 NUmodes) + "_" + name(
977 NSUPmodes) + "_" + name(nNutModes) + "_t";
978
979 if (ITHACAutilities::check_file("./ITHACAoutput/Matrices/" + ct2Str))
980 {
981 ITHACAstream::ReadDenseTensor(ct2Tensor, "./ITHACAoutput/Matrices/", ct2Str);
982 }
983 else
984 {
985 ct2Tensor = turbulenceTensor2(NUmodes, NSUPmodes, nNutModes);
986 }
987
988 if (bcMethod == "penalty")
989 {
990 bcVelVec = bcVelocityVec(NUmodes, NSUPmodes);
991 bcVelMat = bcVelocityMat(NUmodes, NSUPmodes);
992 }
993 }
994 else
995 {
996 L_U_SUPmodes.resize(0);
997
998 if (liftfield.size() != 0)
999 {
1000 for (label k = 0; k < liftfield.size(); k++)
1001 {
1002 L_U_SUPmodes.append(liftfield[k].clone());
1003 }
1004 }
1005
1006 if (NUmodes != 0)
1007 {
1008 for (label k = 0; k < NUmodes; k++)
1009 {
1010 L_U_SUPmodes.append(Umodes[k].clone());
1011 }
1012 }
1013
1014 if (NSUPmodes != 0)
1015 {
1016 for (label k = 0; k < NSUPmodes; k++)
1017 {
1018 L_U_SUPmodes.append(supmodes[k].clone());
1019 }
1020 }
1021
1022 B_matrix = diffusive_term(NUmodes, NPmodes, NSUPmodes);
1023 C_tensor = convective_term_tens(NUmodes, NPmodes, NSUPmodes);
1024 K_matrix = pressure_gradient_term(NUmodes, NPmodes, NSUPmodes);
1025 P_matrix = divergence_term(NUmodes, NPmodes, NSUPmodes);
1026 M_matrix = mass_term(NUmodes, NPmodes, NSUPmodes);
1027 btMatrix = btTurbulence(NUmodes, NSUPmodes);
1028 ct1Tensor = turbulenceTensor1(NUmodes, NSUPmodes, nNutModes);
1029 ct2Tensor = turbulenceTensor2(NUmodes, NSUPmodes, nNutModes);
1030
1031 if (bcMethod == "penalty")
1032 {
1033 bcVelVec = bcVelocityVec(NUmodes, NSUPmodes);
1034 bcVelMat = bcVelocityMat(NUmodes, NSUPmodes);
1035 }
1036 }
1037
1038 // Export the matrices
1039 if (para->exportPython)
1040 {
1041 ITHACAstream::exportMatrix(B_matrix, "B", "python", "./ITHACAoutput/Matrices/");
1042 ITHACAstream::exportMatrix(K_matrix, "K", "python", "./ITHACAoutput/Matrices/");
1043 ITHACAstream::exportMatrix(P_matrix, "P", "python", "./ITHACAoutput/Matrices/");
1044 ITHACAstream::exportMatrix(M_matrix, "M", "python", "./ITHACAoutput/Matrices/");
1045 ITHACAstream::exportMatrix(btMatrix, "bt", "python",
1046 "./ITHACAoutput/Matrices/");
1047 ITHACAstream::exportTensor(C_tensor, "C", "python", "./ITHACAoutput/Matrices/");
1048 ITHACAstream::exportTensor(ct1Tensor, "ct1", "python",
1049 "./ITHACAoutput/Matrices/");
1050 ITHACAstream::exportTensor(ct2Tensor, "ct2", "python",
1051 "./ITHACAoutput/Matrices/");
1052 }
1053
1054 if (para->exportMatlab)
1055 {
1056 ITHACAstream::exportMatrix(B_matrix, "B", "matlab", "./ITHACAoutput/Matrices/");
1057 ITHACAstream::exportMatrix(K_matrix, "K", "matlab", "./ITHACAoutput/Matrices/");
1058 ITHACAstream::exportMatrix(P_matrix, "P", "matlab", "./ITHACAoutput/Matrices/");
1059 ITHACAstream::exportMatrix(M_matrix, "M", "matlab", "./ITHACAoutput/Matrices/");
1060 ITHACAstream::exportMatrix(btMatrix, "bt", "matlab",
1061 "./ITHACAoutput/Matrices/");
1062 ITHACAstream::exportTensor(C_tensor, "C", "matlab", "./ITHACAoutput/Matrices/");
1063 ITHACAstream::exportTensor(ct1Tensor, "ct1", "matlab",
1064 "./ITHACAoutput/Matrices/");
1065 ITHACAstream::exportTensor(ct2Tensor, "ct2", "matlab",
1066 "./ITHACAoutput/Matrices/");
1067 }
1068
1069 if (para->exportTxt)
1070 {
1071 ITHACAstream::exportMatrix(B_matrix, "B", "eigen", "./ITHACAoutput/Matrices/");
1072 ITHACAstream::exportMatrix(K_matrix, "K", "eigen", "./ITHACAoutput/Matrices/");
1073 ITHACAstream::exportMatrix(P_matrix, "P", "eigen", "./ITHACAoutput/Matrices/");
1074 ITHACAstream::exportMatrix(M_matrix, "M", "eigen", "./ITHACAoutput/Matrices/");
1075 ITHACAstream::exportMatrix(btMatrix, "bt", "eigen", "./ITHACAoutput/Matrices/");
1076 ITHACAstream::exportTensor(C_tensor, "C", "eigen", "./ITHACAoutput/Matrices/C");
1077 ITHACAstream::exportTensor(ct1Tensor, "ct1_", "eigen",
1078 "./ITHACAoutput/Matrices/ct1");
1079 ITHACAstream::exportTensor(ct2Tensor, "ct2_", "eigen",
1080 "./ITHACAoutput/Matrices/ct2");
1081 }
1082
1083 bTotalMatrix = B_matrix + btMatrix;
1084 label cSize = NUmodes + NSUPmodes + liftfield.size();
1085 cTotalTensor.resize(cSize, nNutModes, cSize);
1086 cTotalTensor = ct1Tensor + ct2Tensor;
1087 // Get the coeffs for interpolation (the orthonormal one is used because basis are orthogonal)
1088 coeffL2 = ITHACAutilities::getCoeffs(nutFields,
1089 nutModes, nNutModes);
1090 ITHACAstream::exportMatrix(coeffL2, "coeffL2", "python",
1091 "./ITHACAoutput/Matrices/");
1092 ITHACAstream::exportMatrix(coeffL2, "coeffL2", "matlab",
1093 "./ITHACAoutput/Matrices/");
1094 ITHACAstream::exportMatrix(coeffL2, "coeffL2", "eigen",
1095 "./ITHACAoutput/Matrices/");
1096 // Export the matrix
1097 ITHACAstream::SaveDenseMatrix(coeffL2, "./ITHACAoutput/Matrices/",
1098 "coeffL2_nut_" + name(nNutModes));
1099
1100 // Create RBF interpolators for nut coefficient interpolation
1101 rbfSplines.resize(nNutModes);
1102 for (label i = 0; i < nNutModes; i++)
1103 {
1104 // Create ithacaInterpolator instance
1105 rbfSplines[i] = std::make_shared<ithacaInterpolator>(viscDict);
1106
1107 // Prepare training data: X is parameter matrix (transposed), y is coefficient vector
1108 Eigen::MatrixXd X = mu.transpose(); // Now each row is a parameter sample
1109 Eigen::VectorXd y = coeffL2.row(i).transpose(); // Coefficient vector for this mode
1110
1111 rbfSplines[i]->fit(X, y);
1112
1113 Info << "Constructing ithacaInterpolator for mode " << i + 1 << endl;
1114 }
1115 Info<< "Info on interpolators for nut coefficients: "<< endl;
1116 rbfSplines[0]->printInfo();
1117}
SteadyNSTurb.H
Header file of the SteadyNSTurb class.
ITHACAparameters::getInstance
static ITHACAparameters * getInstance()
Gets an instance of ITHACAparameters, to be used if the instance is already existing.
Definition ITHACAparameters.C:56
SteadyNSTurb::coeffL2
Eigen::MatrixXd coeffL2
The matrix of L2 projection coefficients for the eddy viscosity.
Definition SteadyNSTurb.H:114
SteadyNSTurb::btMatrix
Eigen::MatrixXd btMatrix
Turbulent viscosity matrix.
Definition SteadyNSTurb.H:87
SteadyNSTurb::cTotalPPETensor
Eigen::Tensor< double, 3 > cTotalPPETensor
Turbulent total viscosity tensor in the PPE equation.
Definition SteadyNSTurb.H:108
SteadyNSTurb::nNutModes
label nNutModes
Number of viscoisty modes used for the projection.
Definition SteadyNSTurb.H:123
SteadyNSTurb::bTotalMatrix
Eigen::MatrixXd bTotalMatrix
Total B Matrix.
Definition SteadyNSTurb.H:111
SteadyNSTurb::turbulenceTensor2_cache
Eigen::Tensor< double, 3 > turbulenceTensor2_cache(label NUmodes, label NSUPmodes, label nNutModes)
Method to compute one of the turbulence eddy viscosity tensors using the cached procedure.
Definition SteadyNSTurb.C:304
SteadyNSTurb::btTurbulence
Eigen::MatrixXd btTurbulence(label NUmodes, label NSUPmodes)
bt added matrix for the turbulence treatement
Definition SteadyNSTurb.C:495
SteadyNSTurb::_nut
autoPtr< volScalarField > _nut
Eddy viscosity field.
Definition SteadyNSTurb.H:126
SteadyNSTurb::turbulenceTerm1
List< Eigen::MatrixXd > turbulenceTerm1(label NUmodes, label NSUPmodes, label nNutModes)
ct1 added matrix for the turbulence treatement
Definition SteadyNSTurb.C:141
SteadyNSTurb::nutFields
PtrList< volScalarField > nutFields
List of snapshots for the solution for eddy viscosity.
Definition SteadyNSTurb.H:78
SteadyNSTurb::turbulenceTensor1_cache
Eigen::Tensor< double, 3 > turbulenceTensor1_cache(label NUmodes, label NSUPmodes, label nNutModes)
ct1 tensor for the turbulence treatement using the cached procedure
Definition SteadyNSTurb.C:204
SteadyNSTurb::ct2Matrix
List< Eigen::MatrixXd > ct2Matrix
Turbulent viscosity tensor.
Definition SteadyNSTurb.H:90
SteadyNSTurb::ct1PPETensor
Eigen::Tensor< double, 3 > ct1PPETensor
Turbulent viscosity tensor in the PPE equation.
Definition SteadyNSTurb.H:98
SteadyNSTurb::turbulencePPETensor2_cache
Eigen::Tensor< double, 3 > turbulencePPETensor2_cache(label NUmodes, label NSUPmodes, label NPmodes, label nNutModes)
ct2PPE added tensor for the turbulence treatement in the PPE method using the cached procedure
Definition SteadyNSTurb.C:450
SteadyNSTurb::turbulencePPETensor1_cache
Eigen::Tensor< double, 3 > turbulencePPETensor1_cache(label NUmodes, label NSUPmodes, label NPmodes, label nNutModes)
ct1PPE added tensor for the turbulence treatement in the PPE method using the cached procedure
Definition SteadyNSTurb.C:374
SteadyNSTurb::ct1Matrix
List< Eigen::MatrixXd > ct1Matrix
Turbulent viscosity tensor.
Definition SteadyNSTurb.H:94
SteadyNSTurb::turbulenceTerm2
List< Eigen::MatrixXd > turbulenceTerm2(label NUmodes, label NSUPmodes, label nNutModes)
Method to compute one of the turbulence eddy viscosity tensors.
Definition SteadyNSTurb.C:241
SteadyNSTurb::projectSUP
void projectSUP(fileName folder, label NUmodes, label NPmodes, label NSUPmodes, label nNutModes)
Project using a supremizer approach.
Definition SteadyNSTurb.C:855
SteadyNSTurb::ct2PPETensor
Eigen::Tensor< double, 3 > ct2PPETensor
Turbulent viscosity tensor in the PPE equation.
Definition SteadyNSTurb.H:101
SteadyNSTurb::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 SteadyNSTurb.C:342
SteadyNSTurb::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 SteadyNSTurb.C:419
SteadyNSTurb::turbulenceTensor1
Eigen::Tensor< double, 3 > turbulenceTensor1(label NUmodes, label NSUPmodes, label nNutModes)
ct1 tensor for the turbulence treatement
Definition SteadyNSTurb.C:178
SteadyNSTurb::turbulenceTensor2
Eigen::Tensor< double, 3 > turbulenceTensor2(label NUmodes, label NSUPmodes, label nNutModes)
Method to compute one of the turbulence eddy viscosity tensors.
Definition SteadyNSTurb.C:278
SteadyNSTurb::rbfSplines
std::vector< std::shared_ptr< ithacaInterpolator > > rbfSplines
Interpolators for nut coefficient interpolation.
Definition SteadyNSTurb.H:84
SteadyNSTurb::projectPPE
void projectPPE(fileName folder, label NUmodes, label NPmodes, label NSUPmodes, label nNutModes)
Project using a PPE approach.
Definition SteadyNSTurb.C:518
SteadyNSTurb::nutModes
volScalarModes nutModes
List of POD modes for eddy viscosity.
Definition SteadyNSTurb.H:81
SteadyNSTurb::viscDict
dictionary viscDict
The way to compute the eddy viscosity snapshots.
Definition SteadyNSTurb.H:120
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::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::BC4_matrix
Eigen::MatrixXd BC4_matrix
PPE BC4.
Definition steadyNS.H:203
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::_simple
autoPtr< simpleControl > _simple
simpleControl
Definition steadyNS.H:293
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::pressure_BC4
Eigen::MatrixXd pressure_BC4(label NPmodes, label NUmodes)
Term N° 4 given by the additional boundary condition using a PPE approach for time-dependent BCs.
Definition steadyNS.C:1477
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::divergence_term
Eigen::MatrixXd divergence_term(label NUmodes, label NPmodes, label NSUPmodes)
Divergence Term (supremizer approach).
Definition steadyNS.C:1112
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::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::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::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::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::turbulence
autoPtr< incompressible::turbulenceModel > turbulence
Turbulence model.
Definition steadyNS.H:308
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::P_matrix
Eigen::MatrixXd P_matrix
Div of velocity.
Definition steadyNS.H:173
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::bcVelocityMat
List< Eigen::MatrixXd > bcVelocityMat(label NUmodes, label NSUPmodes)
Boundary integral modes on boundary used by the penaly method.
Definition steadyNS.C:1547
steadyNS::_U
autoPtr< volVectorField > _U
Velocity field.
Definition steadyNS.H:272
steadyNS::Pmodes
volScalarModes Pmodes
List of pointers used to form the pressure modes.
Definition steadyNS.H:98
steadyNS::bcMethod
word bcMethod
Boundary Method.
Definition steadyNS.H:326
steadyNS::_p
autoPtr< volScalarField > _p
Pressure field.
Definition steadyNS.H:269
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::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::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
steadyNS.H
Header file of the steadyNS class.
Generated by   doxygen 1.16.1