Loading...
Searching...
No Matches
SteadyNSTurb.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
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 viscCoeff = ITHACAdict->lookupOrDefault<word>("viscCoeff", "RBF");
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 for (label k = 0; k < cSize; ++k)
215 {
216 lapRow[k] = fvc::laplacian(nutModes[j], L_U_SUPmodes[k]);
217 }
218
219 for (label i = 0; i < cSize; ++i)
220 {
221 for (label k = 0; k < cSize; ++k)
222 {
223 const volVectorField& lapField = lapRow[k]();
224 ct1Tensor(i, j, k) = fvc::domainIntegrate(L_U_SUPmodes[i] & lapField).value();
225 }
226 }
227 }
228
229 // Export the tensor
230 if (Pstream::master())
231 {
232 ITHACAstream::SaveDenseTensor(ct1Tensor, "./ITHACAoutput/Matrices/",
233 "ct1_" + name(liftfield.size()) + "_" + name(NUmodes) + "_" + name(
234 NSUPmodes) + "_" + name(nNutModes) + "_t");
235 }
236 return ct1Tensor;
237}
238
239List < Eigen::MatrixXd > SteadyNSTurb::turbulenceTerm2(label NUmodes,
240 label NSUPmodes, label nNutModes)
241{
242 label cSize = NUmodes + NSUPmodes + liftfield.size();
243 List < Eigen::MatrixXd > ct2Matrix;
244 ct2Matrix.setSize(cSize);
245
246 for (label j = 0; j < cSize; j++)
247 {
248 ct2Matrix[j].resize(nNutModes, cSize);
249 ct2Matrix[j] = ct2Matrix[j] * 0;
250 }
251
252 for (label i = 0; i < cSize; i++)
253 {
254 Info << "Filling layer number " << i + 1 << " in the matrix ct2Matrix" << endl;
255
256 for (label j = 0; j < nNutModes; j++)
257 {
258 for (label k = 0; k < cSize; k++)
259 {
260 ct2Matrix[i](j, k) = fvc::domainIntegrate(L_U_SUPmodes[i] & (fvc::div(
261 nutModes[j] * dev((fvc::grad(L_U_SUPmodes[k]))().T())))).value();
262 }
263 }
264 }
265
266 // Export the matrix
267 ITHACAstream::exportMatrix(ct2Matrix, "ct2Matrix", "python",
268 "./ITHACAoutput/Matrices/");
269 ITHACAstream::exportMatrix(ct2Matrix, "ct2Matrix", "matlab",
270 "./ITHACAoutput/Matrices/");
271 ITHACAstream::exportMatrix(ct2Matrix, "ct2Matrix", "eigen",
272 "./ITHACAoutput/Matrices/ct2");
273 return ct2Matrix;
274}
275
276Eigen::Tensor<double, 3> SteadyNSTurb::turbulenceTensor2(label NUmodes,
277 label NSUPmodes, label nNutModes)
278{
279 label cSize = NUmodes + NSUPmodes + liftfield.size();
280 Eigen::Tensor<double, 3> ct2Tensor;
281 ct2Tensor.resize(cSize, nNutModes, cSize);
282
283 for (label i = 0; i < cSize; i++)
284 {
285 for (label j = 0; j < nNutModes; j++)
286 {
287 for (label k = 0; k < cSize; k++)
288 {
289 ct2Tensor(i, j, k) = fvc::domainIntegrate(L_U_SUPmodes[i] & (fvc::div(
290 nutModes[j] * dev((fvc::grad(L_U_SUPmodes[k]))().T())))).value();
291 }
292 }
293 }
294
295 // Export the tensor
296 ITHACAstream::SaveDenseTensor(ct2Tensor, "./ITHACAoutput/Matrices/",
297 "ct2_" + name(liftfield.size()) + "_" + name(NUmodes) + "_" + name(
298 NSUPmodes) + "_" + name(nNutModes) + "_t");
299 return ct2Tensor;
300}
301
302Eigen::Tensor<double, 3> SteadyNSTurb::turbulenceTensor2_cache(label NUmodes,
303 label NSUPmodes, label nNutModes)
304{
305 label cSize = NUmodes + NSUPmodes + liftfield.size();
306 Eigen::Tensor<double, 3> ct2Tensor(cSize, nNutModes, cSize);
307
308 for (label j = 0; j < nNutModes; ++j)
309 {
310 // Cache only one j-row
311 List<tmp<volVectorField>> divRow(cSize);
312 for (label k = 0; k < cSize; ++k)
313 {
314 divRow[k] = fvc::div(nutModes[j] * dev((fvc::grad(L_U_SUPmodes[k]))().T()));
315 }
316
317 for (label i = 0; i < cSize; ++i)
318 {
319 for (label k = 0; k < cSize; ++k)
320 {
321 const volVectorField& divRowField = divRow[k]();
322 ct2Tensor(i, j, k) = fvc::domainIntegrate(L_U_SUPmodes[i] & divRowField).value();
323 }
324 }
325 }
326
327 // Export the tensor
328 if (Pstream::master())
329 {
330 ITHACAstream::SaveDenseTensor(ct2Tensor, "./ITHACAoutput/Matrices/",
331 "ct2_" + name(liftfield.size()) + "_" + name(NUmodes) + "_" + name(
332 NSUPmodes) + "_" + name(nNutModes) + "_t");
333 }
334 return ct2Tensor;
335}
336
337Eigen::Tensor<double, 3> SteadyNSTurb::turbulencePPETensor1(label NUmodes,
338 label NSUPmodes, label NPmodes, label nNutModes)
339{
340 label cSize = NUmodes + NSUPmodes + liftfield.size();
341 Eigen::Tensor<double, 3> ct1PPETensor;
342 ct1PPETensor.resize(NPmodes, nNutModes, cSize);
343
344 for (label i = 0; i < NPmodes; i++)
345 {
346 for (label j = 0; j < nNutModes; j++)
347 {
348 for (label k = 0; k < cSize; k++)
349 {
350 // ct1PPETensor(i, j, k) = fvc::domainIntegrate(2 * Pmodes[i] * (fvc::laplacian(
351 // L_U_SUPmodes[k]) & fvc::grad(nutModes[j]))).value();
352 // ct1PPETensor(i, j, k) = fvc::domainIntegrate(Pmodes[i] * (fvc::div(
353 // fvc::laplacian(
354 // nutModes[j], L_U_SUPmodes[k])))).value();
355 ct1PPETensor(i, j, k) = fvc::domainIntegrate(fvc::grad(Pmodes[i]) & (
356 fvc::laplacian(
357 nutModes[j], L_U_SUPmodes[k]))).value();
358 }
359 }
360 }
361
362 // Export the tensor
363 ITHACAstream::SaveDenseTensor(ct1PPETensor, "./ITHACAoutput/Matrices/",
364 "ct1PPE_" + name(liftfield.size()) + "_" + name(NUmodes) + "_" + name(
365 NSUPmodes) + "_" + name(NPmodes) + "_" + name(nNutModes) + "_t");
366 return ct1PPETensor;
367}
368
369Eigen::Tensor<double, 3> SteadyNSTurb::turbulencePPETensor1_cache(label NUmodes,
370 label NSUPmodes, label NPmodes, label nNutModes)
371{
372 label cSize = NUmodes + NSUPmodes + liftfield.size();
373 Eigen::Tensor<double, 3> ct1PPETensor(NPmodes, nNutModes, cSize);
374
375 List<tmp<volVectorField>> PmodesGrad(NPmodes);
376 for (label i = 0; i < NPmodes; ++i)
377 {
378 PmodesGrad[i] = fvc::grad(Pmodes[i]);
379 }
380
381 for (label j = 0; j < nNutModes; ++j)
382 {
383 // Cache one row (j fixed)
384 List<tmp<volVectorField>> lapRow(cSize);
385 for (label k = 0; k < cSize; ++k)
386 {
387 lapRow[k] = fvc::laplacian(nutModes[j], L_U_SUPmodes[k]);
388 }
389
390 for (label i = 0; i < NPmodes; ++i)
391 {
392 const volVectorField& gradPi = PmodesGrad[i]();
393 for (label k = 0; k < cSize; ++k)
394 {
395 const volVectorField& lapRowField = lapRow[k]();
396 ct1PPETensor(i, j, k) = fvc::domainIntegrate(gradPi & lapRowField).value();
397 }
398 }
399 }
400
401 // Export the tensor
402 if (Pstream::master())
403 {
404 ITHACAstream::SaveDenseTensor(ct1PPETensor, "./ITHACAoutput/Matrices/",
405 "ct1PPE_" + name(liftfield.size()) + "_" + name(NUmodes) + "_" + name(
406 NSUPmodes) + "_" + name(NPmodes) + "_" + name(nNutModes) + "_t");
407 }
408 return ct1PPETensor;
409}
410
411Eigen::Tensor<double, 3> SteadyNSTurb::turbulencePPETensor2(label NUmodes,
412 label NSUPmodes, label NPmodes, label nNutModes)
413{
414 label cSize = NUmodes + NSUPmodes + liftfield.size();
415 Eigen::Tensor<double, 3> ct2PPETensor;
416 ct2PPETensor.resize(NPmodes, nNutModes, cSize);
417
418 for (label i = 0; i < NPmodes; i++)
419 {
420 for (label j = 0; j < nNutModes; j++)
421 {
422 for (label k = 0; k < cSize; k++)
423 {
424 // ct2PPETensor(i, j, k) = fvc::domainIntegrate(Pmodes[i] * (fvc::grad(fvc::grad(
425 // nutModes[j])) && (dev2((fvc::grad(L_U_SUPmodes[k]))() + fvc::grad(
426 // L_U_SUPmodes[k]))().T()))).value();
427 // ct2PPETensor(i, j, k) = fvc::domainIntegrate(Pmodes[i] * ((fvc::div(fvc::div(
428 // nutModes[j] * dev2((fvc::grad(L_U_SUPmodes[k]))().T())))))).value();
429 ct2PPETensor(i, j, k) = fvc::domainIntegrate(fvc::grad(Pmodes[i]) & ((fvc::div(
430 nutModes[j] * dev2((fvc::grad(L_U_SUPmodes[k]))().T()))))).value();
431 }
432 }
433 }
434
435 // Export the tensor
436 ITHACAstream::SaveDenseTensor(ct2PPETensor, "./ITHACAoutput/Matrices/",
437 "ct2PPE_" + name(liftfield.size()) + "_" + name(NUmodes) + "_" + name(
438 NSUPmodes) + "_" + name(NPmodes) + "_" + name(nNutModes) + "_t");
439 return ct2PPETensor;
440}
441
442Eigen::Tensor<double, 3> SteadyNSTurb::turbulencePPETensor2_cache(label NUmodes,
443 label NSUPmodes, label NPmodes, label nNutModes)
444{
445 label cSize = NUmodes + NSUPmodes + liftfield.size();
446 Eigen::Tensor<double, 3> ct2PPETensor(NPmodes, nNutModes, cSize);
447
448 List<tmp<volVectorField>> PmodesGrad(NPmodes);
449 for (label i = 0; i < NPmodes; ++i)
450 {
451 PmodesGrad[i] = fvc::grad(Pmodes[i]);
452 }
453
454 for (label j = 0; j < nNutModes; ++j)
455 {
456 // Cache one row of div fields for this nutMode[j]
457 List<tmp<volVectorField>> divLow(cSize);
458 for (label k = 0; k < cSize; ++k)
459 {
460 divLow[k] = fvc::div(nutModes[j] * dev2((fvc::grad(L_U_SUPmodes[k]))().T()));
461 }
462
463 for (label i = 0; i < NPmodes; ++i)
464 {
465 const volVectorField& gradPi = PmodesGrad[i]();
466 for (label k = 0; k < cSize; ++k)
467 {
468 const volVectorField& divLowField = divLow[k]();
469 ct2PPETensor(i, j, k) = fvc::domainIntegrate(gradPi & divLowField).value();
470 }
471 }
472 }
473
474 // Export the tensor
475 if (Pstream::master())
476 {
477 ITHACAstream::SaveDenseTensor(ct2PPETensor, "./ITHACAoutput/Matrices/",
478 "ct2PPE_" + name(liftfield.size()) + "_" + name(NUmodes) + "_" + name(
479 NSUPmodes) + "_" + name(NPmodes) + "_" + name(nNutModes) + "_t");
480 }
481 return ct2PPETensor;
482}
483
484Eigen::MatrixXd SteadyNSTurb::btTurbulence(label NUmodes, label NSUPmodes)
485{
486 label btSize = NUmodes + NSUPmodes + liftfield.size();
487 Eigen::MatrixXd btMatrix(btSize, btSize);
488 btMatrix = btMatrix * 0;
489
490 // Project everything
491 for (label i = 0; i < btSize; i++)
492 {
493 for (label j = 0; j < btSize; j++)
494 {
495 btMatrix(i, j) = fvc::domainIntegrate(L_U_SUPmodes[i] & (fvc::div(dev((T(
496 fvc::grad(
497 L_U_SUPmodes[j]))))))).value();
498 }
499 }
500
501 // Export the matrix
502 ITHACAstream::SaveDenseMatrix(btMatrix, "./ITHACAoutput/Matrices/",
503 "bt_" + name(liftfield.size()) + "_" + name(NUmodes) + "_" + name(NSUPmodes));
504 return btMatrix;
505}
506
507void SteadyNSTurb::projectPPE(fileName folder, label NU, label NP, label NSUP,
508 label Nnut)
509{
510 NUmodes = NU;
511 NPmodes = NP;
512 NSUPmodes = 0;
513 nNutModes = Nnut;
514 L_U_SUPmodes.resize(0);
515
516 if (liftfield.size() != 0)
517 {
518 for (label k = 0; k < liftfield.size(); k++)
519 {
520 L_U_SUPmodes.append(liftfield[k].clone());
521 }
522 }
523
524 if (NUmodes != 0)
525 {
526 for (label k = 0; k < NUmodes; k++)
527 {
528 L_U_SUPmodes.append(Umodes[k].clone());
529 }
530 }
531
532 if (ITHACAutilities::check_folder("./ITHACAoutput/Matrices/"))
533 {
534 word bStr = "B_" + name(liftfield.size()) + "_" + name(NUmodes) + "_" + name(
535 NSUPmodes);
536
537 if (ITHACAutilities::check_file("./ITHACAoutput/Matrices/" + bStr))
538 {
539 ITHACAstream::ReadDenseMatrix(B_matrix, "./ITHACAoutput/Matrices/", bStr);
540 }
541 else
542 {
543 B_matrix = diffusive_term(NUmodes, NPmodes, NSUPmodes);
544 }
545
546 word btStr = "bt_" + name(liftfield.size()) + "_" + name(NUmodes) + "_" + name(
547 NSUPmodes);
548
549 if (ITHACAutilities::check_file("./ITHACAoutput/Matrices/" + btStr))
550 {
551 ITHACAstream::ReadDenseMatrix(btMatrix, "./ITHACAoutput/Matrices/", btStr);
552 }
553 else
554 {
555 btMatrix = btTurbulence(NUmodes, NSUPmodes);
556 }
557
558 word kStr = "K_" + name(liftfield.size()) + "_" + name(NUmodes) + "_" + name(
559 NSUPmodes) + "_" + name(NPmodes);
560
561 if (ITHACAutilities::check_file("./ITHACAoutput/Matrices/" + kStr))
562 {
563 ITHACAstream::ReadDenseMatrix(K_matrix, "./ITHACAoutput/Matrices/", kStr);
564 }
565 else
566 {
567 K_matrix = pressure_gradient_term(NUmodes, NPmodes, NSUPmodes);
568 }
569
570 word pStr = "P_" + name(liftfield.size()) + "_" + name(NUmodes) + "_" + name(
571 NSUPmodes) + "_" + name(NPmodes);
572
573 if (ITHACAutilities::check_file("./ITHACAoutput/Matrices/" + pStr))
574 {
575 ITHACAstream::ReadDenseMatrix(P_matrix, "./ITHACAoutput/Matrices/", pStr);
576 }
577 else
578 {
579 P_matrix = divergence_term(NUmodes, NPmodes, NSUPmodes);
580 }
581
582 word mStr = "M_" + name(liftfield.size()) + "_" + name(NUmodes) + "_" + name(
583 NSUPmodes);
584
585 if (ITHACAutilities::check_file("./ITHACAoutput/Matrices/" + mStr))
586 {
587 ITHACAstream::ReadDenseMatrix(M_matrix, "./ITHACAoutput/Matrices/", mStr);
588 }
589 else
590 {
591 M_matrix = mass_term(NUmodes, NPmodes, NSUPmodes);
592 }
593
594 word D_str = "D_" + name(NPmodes);
595
596 if (ITHACAutilities::check_file("./ITHACAoutput/Matrices/" + D_str))
597 {
598 ITHACAstream::ReadDenseMatrix(D_matrix, "./ITHACAoutput/Matrices/", D_str);
599 }
600 else
601 {
602 D_matrix = laplacian_pressure(NPmodes);
603 }
604
605 word bc3_str = "BC3_" + name(liftfield.size()) + "_" + name(
606 NUmodes) + "_" + name(NPmodes);
607
608 if (ITHACAutilities::check_file("./ITHACAoutput/Matrices/" + bc3_str))
609 {
610 ITHACAstream::ReadDenseMatrix(BC3_matrix, "./ITHACAoutput/Matrices/", bc3_str);
611 }
612 else
613 {
614 BC3_matrix = pressure_BC3(NUmodes, NPmodes);
615 }
616
617 word bc4_str = "BC4_" + name(liftfield.size()) + "_" + name(
618 NUmodes) + "_" + name(NPmodes);
619
620 if (ITHACAutilities::check_file("./ITHACAoutput/Matrices/" + bc4_str))
621 {
622 ITHACAstream::ReadDenseMatrix(BC4_matrix, "./ITHACAoutput/Matrices/", bc4_str);
623 }
624 else
625 {
626 BC4_matrix = pressure_BC4(NUmodes, NPmodes);
627 }
628
629 word C_str = "C_" + name(liftfield.size()) + "_" + name(NUmodes) + "_" + name(
630 NSUPmodes) + "_t";
631
632 if (ITHACAutilities::check_file("./ITHACAoutput/Matrices/" + C_str))
633 {
634 ITHACAstream::ReadDenseTensor(C_tensor, "./ITHACAoutput/Matrices/", C_str);
635 }
636 else
637 {
638 C_tensor = convective_term_tens(NUmodes, NPmodes, NSUPmodes);
639 }
640
641 word ct1Str = "ct1_" + name(liftfield.size()) + "_" + name(
642 NUmodes) + "_" + name(
643 NSUPmodes) + "_" + name(nNutModes) + "_t";
644
645 if (ITHACAutilities::check_file("./ITHACAoutput/Matrices/" + ct1Str))
646 {
647 ITHACAstream::ReadDenseTensor(ct1Tensor, "./ITHACAoutput/Matrices/", ct1Str);
648 }
649 else
650 {
651 ct1Tensor = turbulenceTensor1(NUmodes, NSUPmodes, nNutModes);
652 }
653
654 word ct2Str = "ct2_" + name(liftfield.size()) + "_" + name(
655 NUmodes) + "_" + name(
656 NSUPmodes) + "_" + name(nNutModes) + "_t";
657
658 if (ITHACAutilities::check_file("./ITHACAoutput/Matrices/" + ct2Str))
659 {
660 ITHACAstream::ReadDenseTensor(ct2Tensor, "./ITHACAoutput/Matrices/", ct2Str);
661 }
662 else
663 {
664 ct2Tensor = turbulenceTensor2(NUmodes, NSUPmodes, nNutModes);
665 }
666
667 word G_str = "G_" + name(liftfield.size()) + "_" + name(NUmodes) + "_" + name(
668 NSUPmodes) + "_" + name(NPmodes) + "_t";
669
670 if (ITHACAutilities::check_file("./ITHACAoutput/Matrices/" + G_str))
671 {
672 ITHACAstream::ReadDenseTensor(gTensor, "./ITHACAoutput/Matrices/", G_str);
673 }
674 else
675 {
676 gTensor = divMomentum(NUmodes, NPmodes);
677 }
678
679 if (bcMethod == "penalty")
680 {
681 bcVelVec = bcVelocityVec(NUmodes, NSUPmodes);
682 bcVelMat = bcVelocityMat(NUmodes, NSUPmodes);
683 }
684 }
685 else
686 {
687 L_U_SUPmodes.resize(0);
688
689 if (liftfield.size() != 0)
690 {
691 for (label k = 0; k < liftfield.size(); k++)
692 {
693 L_U_SUPmodes.append(liftfield[k].clone());
694 }
695 }
696
697 if (NUmodes != 0)
698 {
699 for (label k = 0; k < NUmodes; k++)
700 {
701 L_U_SUPmodes.append(Umodes[k].clone());
702 }
703 }
704
705 if (NSUPmodes != 0)
706 {
707 for (label k = 0; k < NSUPmodes; k++)
708 {
709 L_U_SUPmodes.append(supmodes[k].clone());
710 }
711 }
712
713 B_matrix = diffusive_term(NUmodes, NPmodes, NSUPmodes);
714 C_tensor = convective_term_tens(NUmodes, NPmodes, NSUPmodes);
715 K_matrix = pressure_gradient_term(NUmodes, NPmodes, NSUPmodes);
716 M_matrix = mass_term(NUmodes, NPmodes, NSUPmodes);
717 D_matrix = laplacian_pressure(NPmodes);
718 gTensor = divMomentum(NUmodes, NPmodes);
719 btMatrix = btTurbulence(NUmodes, NSUPmodes);
720 BC3_matrix = pressure_BC3(NUmodes, NPmodes);
721 // BC4_matrix = pressure_BC4(NUmodes, NPmodes);
722 ct1Tensor = turbulenceTensor1(NUmodes, NSUPmodes, nNutModes);
723 ct2Tensor = turbulenceTensor2(NUmodes, NSUPmodes, nNutModes);
724 ct1PPETensor = turbulencePPETensor1(NUmodes, NSUPmodes, NPmodes, nNutModes);
725 ct2PPETensor = turbulencePPETensor2(NUmodes, NSUPmodes, NPmodes, nNutModes);
726
727 if (bcMethod == "penalty")
728 {
729 bcVelVec = bcVelocityVec(NUmodes, NSUPmodes);
730 bcVelMat = bcVelocityMat(NUmodes, NSUPmodes);
731 }
732 }
733
734 // Export the matrices
735 if (para->exportPython)
736 {
737 ITHACAstream::exportMatrix(B_matrix, "B", "python", "./ITHACAoutput/Matrices/");
738 ITHACAstream::exportMatrix(K_matrix, "K", "python", "./ITHACAoutput/Matrices/");
739 ITHACAstream::exportMatrix(D_matrix, "D", "python", "./ITHACAoutput/Matrices/");
740 ITHACAstream::exportMatrix(M_matrix, "M", "python", "./ITHACAoutput/Matrices/");
741 ITHACAstream::exportMatrix(BC3_matrix, "BC3", "python",
742 "./ITHACAoutput/Matrices/");
743 ITHACAstream::exportMatrix(BC4_matrix, "BC4", "python",
744 "./ITHACAoutput/Matrices/");
745 ITHACAstream::exportMatrix(btMatrix, "bt", "python",
746 "./ITHACAoutput/Matrices/");
747 ITHACAstream::exportTensor(C_tensor, "C", "python", "./ITHACAoutput/Matrices/");
748 ITHACAstream::exportTensor(gTensor, "G", "python", "./ITHACAoutput/Matrices/");
749 ITHACAstream::exportTensor(ct1Tensor, "ct1", "python",
750 "./ITHACAoutput/Matrices/");
751 ITHACAstream::exportTensor(ct2Tensor, "ct2", "python",
752 "./ITHACAoutput/Matrices/");
753 ITHACAstream::exportTensor(ct1PPETensor, "ct1PPE", "python",
754 "./ITHACAoutput/Matrices/");
755 ITHACAstream::exportTensor(ct2PPETensor, "ct2PPE", "python",
756 "./ITHACAoutput/Matrices/");
757 }
758
759 if (para->exportMatlab)
760 {
761 ITHACAstream::exportMatrix(B_matrix, "B", "matlab", "./ITHACAoutput/Matrices/");
762 ITHACAstream::exportMatrix(K_matrix, "K", "matlab", "./ITHACAoutput/Matrices/");
763 ITHACAstream::exportMatrix(D_matrix, "D", "matlab", "./ITHACAoutput/Matrices/");
764 ITHACAstream::exportMatrix(M_matrix, "M", "matlab", "./ITHACAoutput/Matrices/");
765 ITHACAstream::exportMatrix(BC3_matrix, "BC3", "matlab",
766 "./ITHACAoutput/Matrices/");
767 ITHACAstream::exportMatrix(BC4_matrix, "BC4", "matlab",
768 "./ITHACAoutput/Matrices/");
769 ITHACAstream::exportMatrix(btMatrix, "bt", "matlab",
770 "./ITHACAoutput/Matrices/");
771 ITHACAstream::exportTensor(gTensor, "G", "matlab", "./ITHACAoutput/Matrices/");
772 ITHACAstream::exportTensor(ct1Tensor, "ct1", "matlab",
773 "./ITHACAoutput/Matrices/");
774 ITHACAstream::exportTensor(ct2Tensor, "ct2", "matlab",
775 "./ITHACAoutput/Matrices/");
776 ITHACAstream::exportTensor(ct1PPETensor, "ct1PPE", "matlab",
777 "./ITHACAoutput/Matrices/");
778 ITHACAstream::exportTensor(ct2PPETensor, "ct2PPE", "matlab",
779 "./ITHACAoutput/Matrices/");
780 }
781
782 if (para->exportTxt)
783 {
784 ITHACAstream::exportMatrix(B_matrix, "B", "eigen", "./ITHACAoutput/Matrices/");
785 ITHACAstream::exportMatrix(K_matrix, "K", "eigen", "./ITHACAoutput/Matrices/");
786 ITHACAstream::exportMatrix(D_matrix, "D", "eigen", "./ITHACAoutput/Matrices/");
787 ITHACAstream::exportMatrix(M_matrix, "M", "eigen", "./ITHACAoutput/Matrices/");
788 ITHACAstream::exportMatrix(BC3_matrix, "BC3", "eigen",
789 "./ITHACAoutput/Matrices/");
790 ITHACAstream::exportMatrix(BC4_matrix, "BC4", "eigen",
791 "./ITHACAoutput/Matrices/");
792 ITHACAstream::exportMatrix(btMatrix, "bt", "eigen", "./ITHACAoutput/Matrices/");
793 ITHACAstream::exportTensor(gTensor, "G", "eigen", "./ITHACAoutput/Matrices/G");
794 ITHACAstream::exportTensor(C_tensor, "C", "eigen", "./ITHACAoutput/Matrices/C");
795 ITHACAstream::exportTensor(ct1Tensor, "ct1_", "eigen",
796 "./ITHACAoutput/Matrices/ct1");
797 ITHACAstream::exportTensor(ct2Tensor, "ct2_", "eigen",
798 "./ITHACAoutput/Matrices/ct2");
799 ITHACAstream::exportTensor(ct1PPETensor, "ct1PPE_", "eigen",
800 "./ITHACAoutput/Matrices/ct1PPE");
801 ITHACAstream::exportTensor(ct2PPETensor, "ct2PPE_", "eigen",
802 "./ITHACAoutput/Matrices/ct2PPE");
803 }
804
805 bTotalMatrix = B_matrix + btMatrix;
806 label cSize = NUmodes + NSUPmodes + liftfield.size();
807 cTotalTensor.resize(cSize, nNutModes, cSize);
808 cTotalTensor = ct1Tensor + ct2Tensor;
809 cTotalPPETensor.resize(NPmodes, nNutModes, cSize);
810 cTotalPPETensor = ct1PPETensor + ct2PPETensor;
811 // Get the coeffs for interpolation (the orthonormal one is used because basis are orthogonal)
812 coeffL2 = ITHACAutilities::getCoeffs(nutFields,
813 nutModes, nNutModes);
814 ITHACAstream::exportMatrix(coeffL2, "coeffL2", "python",
815 "./ITHACAoutput/Matrices/");
816 ITHACAstream::exportMatrix(coeffL2, "coeffL2", "matlab",
817 "./ITHACAoutput/Matrices/");
818 ITHACAstream::exportMatrix(coeffL2, "coeffL2", "eigen",
819 "./ITHACAoutput/Matrices/");
820 // Export the matrix
821 ITHACAstream::SaveDenseMatrix(coeffL2, "./ITHACAoutput/Matrices/",
822 "coeffL2_nut_" + name(nNutModes));
823 samples.resize(nNutModes);
824 rbfSplines.resize(nNutModes);
825 Eigen::MatrixXd weights;
826
827 for (label i = 0; i < nNutModes; i++)
828 {
829 word weightName = "wRBF_N" + name(i + 1) + "_" + name(liftfield.size()) + "_"
830 + name(NUmodes) + "_" + name(NSUPmodes) ;
831
832 if (ITHACAutilities::check_file("./ITHACAoutput/weightsPPE/" + weightName))
833 {
834 samples[i] = new SPLINTER::DataTable(1, 1);
835
836 for (label j = 0; j < coeffL2.cols(); j++)
837 {
838 samples[i]->addSample(mu.row(j), coeffL2(i, j));
839 }
840
841 ITHACAstream::ReadDenseMatrix(weights, "./ITHACAoutput/weightsPPE/",
842 weightName);
843 rbfSplines[i] = new SPLINTER::RBFSpline(* samples[i],
844 SPLINTER::RadialBasisFunctionType::GAUSSIAN, weights);
845 std::cout << "Constructing RadialBasisFunction for mode " << i + 1 << std::endl;
846 }
847 else
848 {
849 samples[i] = new SPLINTER::DataTable(1, 1);
850
851 for (label j = 0; j < coeffL2.cols(); j++)
852 {
853 samples[i]->addSample(mu.row(j), coeffL2(i, j));
854 }
855
856 rbfSplines[i] = new SPLINTER::RBFSpline(* samples[i],
857 SPLINTER::RadialBasisFunctionType::GAUSSIAN);
858 ITHACAstream::SaveDenseMatrix(rbfSplines[i]->weights,
859 "./ITHACAoutput/weightsPPE/", weightName);
860 std::cout << "Constructing RadialBasisFunction for mode " << i + 1 << std::endl;
861 }
862 }
863}
864
865void SteadyNSTurb::projectSUP(fileName folder, label NU, label NP, label NSUP,
866 label Nnut)
867{
868 NUmodes = NU;
869 NPmodes = NP;
870 NSUPmodes = NSUP;
871 nNutModes = Nnut;
872 L_U_SUPmodes.resize(0);
873
874 if (liftfield.size() != 0)
875 {
876 for (label k = 0; k < liftfield.size(); k++)
877 {
878 L_U_SUPmodes.append(liftfield[k].clone());
879 }
880 }
881
882 if (NUmodes != 0)
883 {
884 for (label k = 0; k < NUmodes; k++)
885 {
886 L_U_SUPmodes.append(Umodes[k].clone());
887 }
888 }
889
890 if (NSUPmodes != 0)
891 {
892 for (label k = 0; k < NSUPmodes; k++)
893 {
894 L_U_SUPmodes.append(supmodes[k].clone());
895 }
896 }
897
898 if (ITHACAutilities::check_folder("./ITHACAoutput/Matrices/"))
899 {
900 word bStr = "B_" + name(liftfield.size()) + "_" + name(NUmodes) + "_" + name(
901 NSUPmodes);
902
903 if (ITHACAutilities::check_file("./ITHACAoutput/Matrices/" + bStr))
904 {
905 ITHACAstream::ReadDenseMatrix(B_matrix, "./ITHACAoutput/Matrices/", bStr);
906 }
907 else
908 {
909 B_matrix = diffusive_term(NUmodes, NPmodes, NSUPmodes);
910 }
911
912 word btStr = "bt_" + name(liftfield.size()) + "_" + name(NUmodes) + "_" + name(
913 NSUPmodes);
914
915 if (ITHACAutilities::check_file("./ITHACAoutput/Matrices/" + btStr))
916 {
917 ITHACAstream::ReadDenseMatrix(btMatrix, "./ITHACAoutput/Matrices/", btStr);
918 }
919 else
920 {
921 btMatrix = btTurbulence(NUmodes, NSUPmodes);
922 }
923
924 word kStr = "K_" + name(liftfield.size()) + "_" + name(NUmodes) + "_" + name(
925 NSUPmodes) + "_" + name(NPmodes);
926
927 if (ITHACAutilities::check_file("./ITHACAoutput/Matrices/" + kStr))
928 {
929 ITHACAstream::ReadDenseMatrix(K_matrix, "./ITHACAoutput/Matrices/", kStr);
930 }
931 else
932 {
933 K_matrix = pressure_gradient_term(NUmodes, NPmodes, NSUPmodes);
934 }
935
936 word pStr = "P_" + name(liftfield.size()) + "_" + name(NUmodes) + "_" + name(
937 NSUPmodes) + "_" + name(NPmodes);
938
939 if (ITHACAutilities::check_file("./ITHACAoutput/Matrices/" + pStr))
940 {
941 ITHACAstream::ReadDenseMatrix(P_matrix, "./ITHACAoutput/Matrices/", pStr);
942 }
943 else
944 {
945 P_matrix = divergence_term(NUmodes, NPmodes, NSUPmodes);
946 }
947
948 word mStr = "M_" + name(liftfield.size()) + "_" + name(NUmodes) + "_" + name(
949 NSUPmodes);
950
951 if (ITHACAutilities::check_file("./ITHACAoutput/Matrices/" + mStr))
952 {
953 ITHACAstream::ReadDenseMatrix(M_matrix, "./ITHACAoutput/Matrices/", mStr);
954 }
955 else
956 {
957 M_matrix = mass_term(NUmodes, NPmodes, NSUPmodes);
958 }
959
960 word C_str = "C_" + name(liftfield.size()) + "_" + name(NUmodes) + "_" + name(
961 NSUPmodes) + "_t";
962
963 if (ITHACAutilities::check_file("./ITHACAoutput/Matrices/" + C_str))
964 {
965 ITHACAstream::ReadDenseTensor(C_tensor, "./ITHACAoutput/Matrices/", C_str);
966 }
967 else
968 {
969 C_tensor = convective_term_tens(NUmodes, NPmodes, NSUPmodes);
970 }
971
972 word ct1Str = "ct1_" + name(liftfield.size()) + "_" + name(
973 NUmodes) + "_" + name(
974 NSUPmodes) + "_" + name(nNutModes) + "_t";
975
976 if (ITHACAutilities::check_file("./ITHACAoutput/Matrices/" + ct1Str))
977 {
978 ITHACAstream::ReadDenseTensor(ct1Tensor, "./ITHACAoutput/Matrices/", ct1Str);
979 }
980 else
981 {
982 ct1Tensor = turbulenceTensor1(NUmodes, NSUPmodes, nNutModes);
983 }
984
985 word ct2Str = "ct2_" + name(liftfield.size()) + "_" + name(
986 NUmodes) + "_" + name(
987 NSUPmodes) + "_" + name(nNutModes) + "_t";
988
989 if (ITHACAutilities::check_file("./ITHACAoutput/Matrices/" + ct2Str))
990 {
991 ITHACAstream::ReadDenseTensor(ct2Tensor, "./ITHACAoutput/Matrices/", ct2Str);
992 }
993 else
994 {
995 ct2Tensor = turbulenceTensor2(NUmodes, NSUPmodes, nNutModes);
996 }
997
998 if (bcMethod == "penalty")
999 {
1000 bcVelVec = bcVelocityVec(NUmodes, NSUPmodes);
1001 bcVelMat = bcVelocityMat(NUmodes, NSUPmodes);
1002 }
1003 }
1004 else
1005 {
1006 L_U_SUPmodes.resize(0);
1007
1008 if (liftfield.size() != 0)
1009 {
1010 for (label k = 0; k < liftfield.size(); k++)
1011 {
1012 L_U_SUPmodes.append(liftfield[k].clone());
1013 }
1014 }
1015
1016 if (NUmodes != 0)
1017 {
1018 for (label k = 0; k < NUmodes; k++)
1019 {
1020 L_U_SUPmodes.append(Umodes[k].clone());
1021 }
1022 }
1023
1024 if (NSUPmodes != 0)
1025 {
1026 for (label k = 0; k < NSUPmodes; k++)
1027 {
1028 L_U_SUPmodes.append(supmodes[k].clone());
1029 }
1030 }
1031
1032 B_matrix = diffusive_term(NUmodes, NPmodes, NSUPmodes);
1033 C_tensor = convective_term_tens(NUmodes, NPmodes, NSUPmodes);
1034 K_matrix = pressure_gradient_term(NUmodes, NPmodes, NSUPmodes);
1035 P_matrix = divergence_term(NUmodes, NPmodes, NSUPmodes);
1036 M_matrix = mass_term(NUmodes, NPmodes, NSUPmodes);
1037 btMatrix = btTurbulence(NUmodes, NSUPmodes);
1038 ct1Tensor = turbulenceTensor1(NUmodes, NSUPmodes, nNutModes);
1039 ct2Tensor = turbulenceTensor2(NUmodes, NSUPmodes, nNutModes);
1040
1041 if (bcMethod == "penalty")
1042 {
1043 bcVelVec = bcVelocityVec(NUmodes, NSUPmodes);
1044 bcVelMat = bcVelocityMat(NUmodes, NSUPmodes);
1045 }
1046 }
1047
1048 // Export the matrices
1049 if (para->exportPython)
1050 {
1051 ITHACAstream::exportMatrix(B_matrix, "B", "python", "./ITHACAoutput/Matrices/");
1052 ITHACAstream::exportMatrix(K_matrix, "K", "python", "./ITHACAoutput/Matrices/");
1053 ITHACAstream::exportMatrix(P_matrix, "P", "python", "./ITHACAoutput/Matrices/");
1054 ITHACAstream::exportMatrix(M_matrix, "M", "python", "./ITHACAoutput/Matrices/");
1055 ITHACAstream::exportMatrix(btMatrix, "bt", "python",
1056 "./ITHACAoutput/Matrices/");
1057 ITHACAstream::exportTensor(C_tensor, "C", "python", "./ITHACAoutput/Matrices/");
1058 ITHACAstream::exportTensor(ct1Tensor, "ct1", "python",
1059 "./ITHACAoutput/Matrices/");
1060 ITHACAstream::exportTensor(ct2Tensor, "ct2", "python",
1061 "./ITHACAoutput/Matrices/");
1062 }
1063
1064 if (para->exportMatlab)
1065 {
1066 ITHACAstream::exportMatrix(B_matrix, "B", "matlab", "./ITHACAoutput/Matrices/");
1067 ITHACAstream::exportMatrix(K_matrix, "K", "matlab", "./ITHACAoutput/Matrices/");
1068 ITHACAstream::exportMatrix(P_matrix, "P", "matlab", "./ITHACAoutput/Matrices/");
1069 ITHACAstream::exportMatrix(M_matrix, "M", "matlab", "./ITHACAoutput/Matrices/");
1070 ITHACAstream::exportMatrix(btMatrix, "bt", "matlab",
1071 "./ITHACAoutput/Matrices/");
1072 ITHACAstream::exportTensor(C_tensor, "C", "matlab", "./ITHACAoutput/Matrices/");
1073 ITHACAstream::exportTensor(ct1Tensor, "ct1", "matlab",
1074 "./ITHACAoutput/Matrices/");
1075 ITHACAstream::exportTensor(ct2Tensor, "ct2", "matlab",
1076 "./ITHACAoutput/Matrices/");
1077 }
1078
1079 if (para->exportTxt)
1080 {
1081 ITHACAstream::exportMatrix(B_matrix, "B", "eigen", "./ITHACAoutput/Matrices/");
1082 ITHACAstream::exportMatrix(K_matrix, "K", "eigen", "./ITHACAoutput/Matrices/");
1083 ITHACAstream::exportMatrix(P_matrix, "P", "eigen", "./ITHACAoutput/Matrices/");
1084 ITHACAstream::exportMatrix(M_matrix, "M", "eigen", "./ITHACAoutput/Matrices/");
1085 ITHACAstream::exportMatrix(btMatrix, "bt", "eigen", "./ITHACAoutput/Matrices/");
1086 ITHACAstream::exportTensor(C_tensor, "C", "eigen", "./ITHACAoutput/Matrices/C");
1087 ITHACAstream::exportTensor(ct1Tensor, "ct1_", "eigen",
1088 "./ITHACAoutput/Matrices/ct1");
1089 ITHACAstream::exportTensor(ct2Tensor, "ct2_", "eigen",
1090 "./ITHACAoutput/Matrices/ct2");
1091 }
1092
1093 bTotalMatrix = B_matrix + btMatrix;
1094 label cSize = NUmodes + NSUPmodes + liftfield.size();
1095 cTotalTensor.resize(cSize, nNutModes, cSize);
1096 cTotalTensor = ct1Tensor + ct2Tensor;
1097 // Get the coeffs for interpolation (the orthonormal one is used because basis are orthogonal)
1098 coeffL2 = ITHACAutilities::getCoeffs(nutFields,
1099 nutModes, nNutModes);
1100 ITHACAstream::exportMatrix(coeffL2, "coeffL2", "python",
1101 "./ITHACAoutput/Matrices/");
1102 ITHACAstream::exportMatrix(coeffL2, "coeffL2", "matlab",
1103 "./ITHACAoutput/Matrices/");
1104 ITHACAstream::exportMatrix(coeffL2, "coeffL2", "eigen",
1105 "./ITHACAoutput/Matrices/");
1106 // Export the matrix
1107 ITHACAstream::SaveDenseMatrix(coeffL2, "./ITHACAoutput/Matrices/",
1108 "coeffL2_nut_" + name(nNutModes));
1109 samples.resize(nNutModes);
1110 rbfSplines.resize(nNutModes);
1111 Eigen::MatrixXd weights;
1112
1113 for (label i = 0; i < nNutModes; i++)
1114 {
1115 word weightName = "wRBF_N" + name(i + 1) + "_" + name(liftfield.size()) + "_"
1116 + name(NUmodes) + "_" + name(NSUPmodes) ;
1117
1118 if (ITHACAutilities::check_file("./ITHACAoutput/weightsSUP/" + weightName))
1119 {
1120 samples[i] = new SPLINTER::DataTable(1, 1);
1121
1122 for (label j = 0; j < coeffL2.cols(); j++)
1123 {
1124 samples[i]->addSample(mu.row(j), coeffL2(i, j));
1125 }
1126
1127 ITHACAstream::ReadDenseMatrix(weights, "./ITHACAoutput/weightsSUP/",
1128 weightName);
1129 rbfSplines[i] = new SPLINTER::RBFSpline(* samples[i],
1130 SPLINTER::RadialBasisFunctionType::GAUSSIAN, weights);
1131 std::cout << "Constructing RadialBasisFunction for mode " << i + 1 << std::endl;
1132 }
1133 else
1134 {
1135 samples[i] = new SPLINTER::DataTable(1, 1);
1136
1137 for (label j = 0; j < coeffL2.cols(); j++)
1138 {
1139 samples[i]->addSample(mu.row(j), coeffL2(i, j));
1140 }
1141
1142 rbfSplines[i] = new SPLINTER::RBFSpline(* samples[i],
1143 SPLINTER::RadialBasisFunctionType::GAUSSIAN);
1144 ITHACAstream::SaveDenseMatrix(rbfSplines[i]->weights,
1145 "./ITHACAoutput/weightsSUP/", weightName);
1146 ITHACAstream::exportMatrix(rbfSplines[i]->weights,
1147 "wRBF_" + name(i), "eigen",
1148 "./ITHACAoutput/weightsSUP/"
1149 );
1150 std::cout << "Constructing RadialBasisFunction for mode " << i + 1 << std::endl;
1151 }
1152 }
1153}
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
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:116
SteadyNSTurb::btMatrix
Eigen::MatrixXd btMatrix
Turbulent viscosity matrix.
Definition SteadyNSTurb.H:89
SteadyNSTurb::cTotalPPETensor
Eigen::Tensor< double, 3 > cTotalPPETensor
Turbulent total viscosity tensor in the PPE equation.
Definition SteadyNSTurb.H:110
SteadyNSTurb::nNutModes
label nNutModes
Number of viscoisty modes used for the projection.
Definition SteadyNSTurb.H:125
SteadyNSTurb::bTotalMatrix
Eigen::MatrixXd bTotalMatrix
Total B Matrix.
Definition SteadyNSTurb.H:113
SteadyNSTurb::rbfSplines
std::vector< SPLINTER::RBFSpline * > rbfSplines
Create a RBF splines for interpolation.
Definition SteadyNSTurb.H:86
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:302
SteadyNSTurb::btTurbulence
Eigen::MatrixXd btTurbulence(label NUmodes, label NSUPmodes)
bt added matrix for the turbulence treatement
Definition SteadyNSTurb.C:484
SteadyNSTurb::_nut
autoPtr< volScalarField > _nut
Eddy viscosity field.
Definition SteadyNSTurb.H:128
SteadyNSTurb::samples
std::vector< SPLINTER::DataTable * > samples
Create a samples for interpolation.
Definition SteadyNSTurb.H:83
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:77
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:92
SteadyNSTurb::ct1PPETensor
Eigen::Tensor< double, 3 > ct1PPETensor
Turbulent viscosity tensor in the PPE equation.
Definition SteadyNSTurb.H:100
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:442
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:369
SteadyNSTurb::ct1Tensor
Eigen::Tensor< double, 3 > ct1Tensor
Definition SteadyNSTurb.H:97
SteadyNSTurb::ct1Matrix
List< Eigen::MatrixXd > ct1Matrix
Turbulent viscosity tensor.
Definition SteadyNSTurb.H:96
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:239
SteadyNSTurb::projectSUP
void projectSUP(fileName folder, label NUmodes, label NPmodes, label NSUPmodes, label nNutModes)
Project using a supremizer approach.
Definition SteadyNSTurb.C:865
SteadyNSTurb::ct2PPETensor
Eigen::Tensor< double, 3 > ct2PPETensor
Turbulent viscosity tensor in the PPE equation.
Definition SteadyNSTurb.H:103
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:337
SteadyNSTurb::cTotalTensor
Eigen::Tensor< double, 3 > cTotalTensor
Definition SteadyNSTurb.H:107
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:411
SteadyNSTurb::turbulenceTensor1
Eigen::Tensor< double, 3 > turbulenceTensor1(label NUmodes, label NSUPmodes, label nNutModes)
ct1 tensor for the turbulence treatement
Definition SteadyNSTurb.C:178
SteadyNSTurb::ct2Tensor
Eigen::Tensor< double, 3 > ct2Tensor
Definition SteadyNSTurb.H:93
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:276
SteadyNSTurb::projectPPE
void projectPPE(fileName folder, label NUmodes, label NPmodes, label NSUPmodes, label nNutModes)
Project using a PPE approach.
Definition SteadyNSTurb.C:507
SteadyNSTurb::viscCoeff
word viscCoeff
The way to compute the eddy viscosity snapshots.
Definition SteadyNSTurb.H:122
SteadyNSTurb::nutModes
volScalarModes nutModes
List of POD modes for eddy viscosity.
Definition SteadyNSTurb.H:80
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::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:1487
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:1452
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: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::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::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: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::bcVelocityMat
List< Eigen::MatrixXd > bcVelocityMat(label NUmodes, label NSUPmodes)
Boundary integral modes on boundary used by the penaly method.
Definition steadyNS.C:1519
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
T
volScalarField & T
Definition createT.H:46
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::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
simple
simpleControl simple(mesh)
phi
surfaceScalarField & phi
Definition createFields.H:68
U
volVectorField & U
Definition createFields.H:42
p
volScalarField & p
Definition createFields.H:42
NLsolveSteadyNSTurb.H
createFields.H
createFvOptions.H
laminarTransport
singlePhaseTransportModel & laminarTransport
Definition createFields.H:86
i
label i
Definition pEqn.H:46
steadyNS.H
Header file of the steadyNS class.
Generated by   doxygen 1.13.2