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