Loading...
Searching...
No Matches
ITHACAPOD.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
33
34#include "ITHACAPOD.H"
35#include "EigenFunctions.H"
36
37namespace ITHACAPOD
38{
39
40template<class Type, template<class> class PatchField, class GeoMesh>
41void getNestedSnapshotMatrix(
42 PtrList<GeometricField<Type, PatchField, GeoMesh>>& snapshots,
43 PtrList<GeometricField<Type, PatchField, GeoMesh>>& ModesGlobal,
44 word fieldName,
45 label Npar, label NnestedOut)
46{
47 ITHACAparameters* para(ITHACAparameters::getInstance());
48 List<PtrList<GeometricField<Type, PatchField, GeoMesh>>>
49 SnapMatrixNested;
50 label Nt = snapshots.size() / Npar;
51 SnapMatrixNested.setSize(Nt);
52
53 for (label i = 0; i < Npar; i++)
54 {
55 SnapMatrixNested[i].resize(Nt);
56
57 for (label j = 0; j < Nt; j++)
58 {
59 SnapMatrixNested[i].set(j, snapshots[j + Nt * i].clone());
60 }
61 }
62
63 List<PtrList<GeometricField<Type, PatchField, GeoMesh>>> UModesNested;
64 PtrList<GeometricField<Type, PatchField, GeoMesh>> y;
65 UModesNested.setSize(Nt);
66
67 for (label i = 0; i < Npar; i++)
68 {
69 getWeightedModes(SnapMatrixNested[i], UModesNested[i], 0, 0, 0,
70 NnestedOut);
71 }
72
73 for (label i = 0; i < Npar; i++)
74 {
75 y = UModesNested[i];
76
77 for (label j = 0; j < y.size(); j++)
78 {
79 ModesGlobal.append(y[j].clone());
80 }
81 }
82}
83
84template void getNestedSnapshotMatrix(
85 PtrList<volScalarField>& snapshots, PtrList<volScalarField>& ModesGlobal,
86 word fieldName, label Npar, label NnestedOut);
87
88template void getNestedSnapshotMatrix(
89 PtrList<volVectorField>& snapshots, PtrList<volVectorField>& ModesGlobal,
90 word fieldName, label Npar, label NnestedOut);
91
92template<class Type, template<class> class PatchField, class GeoMesh>
93void getModes(
94 PtrList<GeometricField<Type, PatchField, GeoMesh>>& snapshots,
95 PtrList<GeometricField<Type, PatchField, GeoMesh>>& modes,
96 word fieldName, bool podex, bool supex, bool sup, label nmodes,
97 bool correctBC)
98{
99 ITHACAparameters* para(ITHACAparameters::getInstance());
100 word PODkey = "POD_" + fieldName;
101 word PODnorm = para->ITHACAdict->lookupOrDefault<word>(PODkey, "L2");
102 M_Assert(PODnorm == "L2" ||
103 PODnorm == "Frobenius", "The PODnorm can be only L2 or Frobenius");
104 Info << "Performing POD for " << fieldName << " using the " << PODnorm <<
105 " norm" << endl;
106
107 if ((podex == 0 && sup == 0) || (supex == 0 && sup == 1))
108 {
109 if (para->eigensolver == "spectra" )
110 {
111 if (nmodes == 0)
112 {
113 nmodes = snapshots.size() - 2;
114 }
115
116 M_Assert(nmodes <= snapshots.size() - 2,
117 "The number of requested modes cannot be bigger than the number of Snapshots - 2");
118 }
119 else
120 {
121 if (nmodes == 0)
122 {
123 nmodes = snapshots.size();
124 }
125
126 M_Assert(nmodes <= snapshots.size(),
127 "The number of requested modes cannot be bigger than the number of Snapshots");
128 }
129
130 Eigen::MatrixXd SnapMatrix = Foam2Eigen::PtrList2Eigen(snapshots);
131 List<Eigen::MatrixXd> SnapMatrixBC = Foam2Eigen::PtrList2EigenBC(snapshots);
132 label NBC = snapshots[0].boundaryField().size();
133 Eigen::MatrixXd _corMatrix;
134
135 if (PODnorm == "L2")
136 {
137 _corMatrix = ITHACAutilities::getMassMatrix(snapshots);
138 }
139 else if (PODnorm == "Frobenius")
140 {
141 _corMatrix = ITHACAutilities::getMassMatrix(snapshots, 0, false);
142 }
143
144 if (Pstream::parRun())
145 {
146 reduce(_corMatrix, sumOp<Eigen::MatrixXd>());
147 }
148
149 Eigen::VectorXd eigenValueseig;
150 Eigen::MatrixXd eigenVectoreig;
151 modes.resize(nmodes);
152 Info << "####### Performing the POD using EigenDecomposition " <<
153 fieldName << " #######" << endl;
154 label ncv = snapshots.size();
155 Spectra::DenseSymMatProd<double> op(_corMatrix);
156 Eigen::SelfAdjointEigenSolver<Eigen::MatrixXd> esEg;
157
158 if (para->eigensolver == "spectra")
159 {
160 Spectra::SymEigsSolver<double, Spectra::LARGEST_ALGE, Spectra::DenseSymMatProd<double>>
161 es(&op, nmodes, ncv);
162 std::cout << "Using Spectra EigenSolver " << std::endl;
163 es.init();
164 es.compute(1000, 1e-10, Spectra::LARGEST_ALGE);
165 M_Assert(es.info() == Spectra::SUCCESSFUL,
166 "The Eigenvalue Decomposition did not succeed");
167 eigenVectoreig = es.eigenvectors().real();
168 eigenValueseig = es.eigenvalues().real();
169 }
170 else if (para->eigensolver == "eigen")
171 {
172 std::cout << "Using Eigen EigenSolver " << std::endl;
173 esEg.compute(_corMatrix);
174 M_Assert(esEg.info() == Eigen::Success,
175 "The Eigenvalue Decomposition did not succeed");
176 eigenVectoreig = esEg.eigenvectors().real().rowwise().reverse().leftCols(
177 nmodes);
178 eigenValueseig = esEg.eigenvalues().real().array().reverse();
179 }
180
181 if (eigenValueseig.array().minCoeff() < 0)
182 {
183 eigenValueseig = eigenValueseig.array() + 2 * abs(
184 eigenValueseig.array().minCoeff());
185 }
186
187 Info << "####### End of the POD for " << snapshots[0].name() << " #######" <<
188 endl;
189 //Eigen::VectorXd eigenValueseigLam =
190 // eigenValueseig.real().array().abs().cwiseInverse().sqrt() ;
191 //Eigen::MatrixXd modesEig = (SnapMatrix * eigenVectoreig) *
192 // eigenValueseigLam.head(nmodes).asDiagonal();
193 Eigen::MatrixXd modesEig = (SnapMatrix * eigenVectoreig);
194 // Computing Normalization factors of the POD Modes
195 Eigen::VectorXd V = ITHACAutilities::getMassMatrixFV(snapshots[0]);
196 Eigen::MatrixXd normFact(nmodes, 1);
197
198 for (label i = 0; i < nmodes; i++)
199 {
200 if (PODnorm == "L2")
201 {
202 normFact(i, 0) = std::sqrt((modesEig.col(i).transpose() * V.asDiagonal() *
203 modesEig.col(i))(0, 0));
204
205 if (Pstream::parRun())
206 {
207 normFact(i, 0) = (modesEig.col(i).transpose() * V.asDiagonal() *
208 modesEig.col(i))(0, 0);
209 }
210 }
211 else if (PODnorm == "Frobenius")
212 {
213 normFact(i, 0) = std::sqrt((modesEig.col(i).transpose() * modesEig.col(i))(0,
214 0));
215
216 if (Pstream::parRun())
217 {
218 normFact(i, 0) = (modesEig.col(i).transpose() * modesEig.col(i))(0, 0);
219 }
220 }
221 }
222
223 if (Pstream::parRun())
224 {
225 reduce(normFact, sumOp<Eigen::MatrixXd>());
226 }
227
228 if (Pstream::parRun())
229 {
230 normFact = normFact.cwiseSqrt();
231 }
232
233 List<Eigen::MatrixXd> modesEigBC;
234 modesEigBC.resize(NBC);
235
236 for (label i = 0; i < NBC; i++)
237 {
238 modesEigBC[i] = (SnapMatrixBC[i] * eigenVectoreig);
239 }
240
241 std::cout << normFact << std::endl;
242
243 for (label i = 0; i < nmodes; i++)
244 {
245 modesEig.col(i) = modesEig.col(i).array() / normFact(i, 0);
246
247 for (label j = 0; j < NBC; j++)
248 {
249 modesEigBC[j].col(i) = modesEigBC[j].col(i).array() / normFact(i, 0);
250 }
251 }
252
253 for (label i = 0; i < modes.size(); i++)
254 {
255 GeometricField<Type, PatchField, GeoMesh> tmp2(snapshots[0].name(),
256 snapshots[0]);
257 Eigen::VectorXd vec = modesEig.col(i);
258 tmp2 = Foam2Eigen::Eigen2field(tmp2, vec, correctBC);
259
260 for (label k = 0; k < NBC; k++)
261 {
262 ITHACAutilities::assignBC(tmp2, k, modesEigBC[k].col(i));
263 }
264
265 modes.set(i, tmp2.clone());
266 }
267
268 eigenValueseig = eigenValueseig / eigenValueseig.sum();
269 Eigen::VectorXd cumEigenValues(eigenValueseig);
270
271 for (label j = 1; j < cumEigenValues.size(); ++j)
272 {
273 cumEigenValues(j) += cumEigenValues(j - 1);
274 }
275
276 Info << "####### Saving the POD bases for " << snapshots[0].name() <<
277 " #######" << endl;
278
279 if (sup)
280 {
281 ITHACAstream::exportFields(modes, "./ITHACAoutput/supremizer/",
282 snapshots[0].name());
283 }
284 else
285 {
286 ITHACAstream::exportFields(modes, "./ITHACAoutput/POD/", snapshots[0].name());
287 }
288
289 Eigen::saveMarketVector(eigenValueseig,
290 "./ITHACAoutput/POD/Eigenvalues_" + snapshots[0].name(), para->precision,
291 para->outytpe);
292 Eigen::saveMarketVector(cumEigenValues,
293 "./ITHACAoutput/POD/CumEigenvalues_" + snapshots[0].name(), para->precision,
294 para->outytpe);
295 }
296 else
297 {
298 Info << "Reading the existing modes" << endl;
299
300 if (sup == 1)
301 {
302 ITHACAstream::read_fields(modes, fieldName + "sup",
303 "./ITHACAoutput/supremizer/");
304 }
305 else
306 {
307 ITHACAstream::read_fields(modes, fieldName, "./ITHACAoutput/POD/");
308 }
309 }
310}
311
312template void getModes(
313 PtrList<volVectorField>& snapshots, PtrList<volVectorField>& modes,
314 word fieldName, bool podex, bool supex, bool sup, label nmodes,
315 bool correctBC);
316
317template void getModes(
318 PtrList<volScalarField>& snapshots, PtrList<volScalarField>& modes,
319 word fieldName, bool podex, bool supex, bool sup, label nmodes,
320 bool correctBC);
321
322template void getModes(
323 PtrList<surfaceScalarField>& snapshots, PtrList<surfaceScalarField>& modes,
324 word fieldName, bool podex, bool supex, bool sup, label nmodes,
325 bool correctBC);
326
327template<class Type, template<class> class PatchField, class GeoMesh>
328void getWeightedModes(
329 PtrList<GeometricField<Type, PatchField, GeoMesh>>& snapshots,
330 PtrList<GeometricField<Type, PatchField, GeoMesh>>& modes,
331 word fieldName, bool podex, bool supex, bool sup, label nmodes,
332 bool correctBC)
333{
334 ITHACAparameters* para(ITHACAparameters::getInstance());
335
336 if (nmodes == 0)
337 {
338 nmodes = snapshots.size() - 2;
339 }
340
341 M_Assert(nmodes <= snapshots.size() - 2,
342 "The number of requested modes cannot be bigger than the number of Snapshots - 2");
343
344 if ((podex == 0 && sup == 0) || (supex == 0 && sup == 1))
345 {
346 Eigen::MatrixXd SnapMatrix = Foam2Eigen::PtrList2Eigen(snapshots);
347 List<Eigen::MatrixXd> SnapMatrixBC = Foam2Eigen::PtrList2EigenBC(snapshots);
348 label NBC = snapshots[0].boundaryField().size();
349 Eigen::MatrixXd _corMatrix = ITHACAutilities::getMassMatrix(snapshots);
350 Eigen::VectorXd eigenValueseig;
351 Eigen::MatrixXd eigenVectoreig;
352 modes.resize(nmodes);
353 Info << "####### Performing the POD using EigenDecomposition " <<
354 snapshots[0].name() << " #######" << endl;
355 label ncv = snapshots.size();
356 Spectra::DenseSymMatProd<double> op(_corMatrix);
357 Eigen::SelfAdjointEigenSolver<Eigen::MatrixXd> esEg;
358 Spectra::SymEigsSolver<double, Spectra::LARGEST_ALGE, Spectra::DenseSymMatProd<double>>
359 es(&op, nmodes, ncv);
360
361 if (para->eigensolver == "spectra")
362 {
363 std::cout << "Using Spectra EigenSolver " << std::endl;
364 es.init();
365 es.compute(1000, 1e-10, Spectra::LARGEST_ALGE);
366 M_Assert(es.info() == Spectra::SUCCESSFUL,
367 "The Eigenvalue Decomposition did not succeed");
368 eigenVectoreig = es.eigenvectors().real();
369 eigenValueseig = es.eigenvalues().real();
370 }
371 else if (para->eigensolver == "eigen")
372 {
373 std::cout << "Using Eigen EigenSolver " << std::endl;
374 esEg.compute(_corMatrix);
375 M_Assert(esEg.info() == Eigen::Success,
376 "The Eigenvalue Decomposition did not succeed");
377 eigenVectoreig = esEg.eigenvectors().real().rowwise().reverse().leftCols(
378 nmodes);
379 eigenValueseig = esEg.eigenvalues().real().reverse();
380 }
381
382 Info << "####### End of the POD for " << snapshots[0].name() << " #######" <<
383 endl;
384 Eigen::VectorXd eigenValueseigLam =
385 eigenValueseig.real().array().cwiseInverse().sqrt() ;
386 Eigen::VectorXd eigenValueseigWeigted = eigenValueseig.head(
387 nmodes).real().array() ;
388 Eigen::MatrixXd modesEig = (SnapMatrix * eigenVectoreig) *
389 eigenValueseigLam.head(nmodes).asDiagonal() *
390 eigenValueseigWeigted.asDiagonal();
391 List<Eigen::MatrixXd> modesEigBC;
392 modesEigBC.resize(NBC);
393
394 for (label i = 0; i < NBC; i++)
395 {
396 modesEigBC[i] = (SnapMatrixBC[i] * eigenVectoreig) *
397 eigenValueseigLam.asDiagonal() * eigenValueseigWeigted.asDiagonal();
398 }
399
400 for (label i = 0; i < modes.size(); i++)
401 {
402 GeometricField<Type, PatchField, GeoMesh> tmp2(snapshots[0].name(),
403 snapshots[0]);
404 Eigen::VectorXd vec = modesEig.col(i);
405 tmp2 = Foam2Eigen::Eigen2field(tmp2, vec, correctBC);
406
407 for (label k = 0; k < tmp2.boundaryField().size(); k++)
408 {
409 ITHACAutilities::assignBC(tmp2, k, modesEigBC[k].col(i));
410 }
411
412 modes.set(i, tmp2.clone());
413 }
414
415 eigenValueseig = eigenValueseig / eigenValueseig.sum();
416 Eigen::VectorXd cumEigenValues(eigenValueseig);
417
418 for (label j = 1; j < cumEigenValues.size(); ++j)
419 {
420 cumEigenValues(j) += cumEigenValues(j - 1);
421 }
422
423 Info << "####### Saving the POD bases for " << snapshots[0].name() <<
424 " #######" << endl;
425
426 //exportBases(modes, snapshots, sup);
427 if (sup)
428 {
429 ITHACAstream::exportFields(modes, "./ITHACAoutput/supremizer/",
430 snapshots[0].name());
431 }
432 else
433 {
434 ITHACAstream::exportFields(modes, "./ITHACAoutput/POD/", snapshots[0].name());
435 }
436
437 Eigen::saveMarketVector(eigenValueseig,
438 "./ITHACAoutput/POD/Eigenvalues_" + snapshots[0].name(), para->precision,
439 para->outytpe);
440 Eigen::saveMarketVector(cumEigenValues,
441 "./ITHACAoutput/POD/CumEigenvalues_" + snapshots[0].name(), para->precision,
442 para->outytpe);
443 }
444 else
445 {
446 Info << "Reading the existing modes" << endl;
447
448 if (sup == 1)
449 {
450 ITHACAstream::read_fields(modes, fieldName + "sup",
451 "./ITHACAoutput/supremizer/");
452 }
453 else
454 {
455 ITHACAstream::read_fields(modes, fieldName, "./ITHACAoutput/POD/");
456 }
457 }
458}
459
460template void getWeightedModes(
461 PtrList<volScalarField>& snapshots, PtrList<volScalarField>& modes,
462 word fieldName, bool podex, bool supex, bool sup, label nmodes,
463 bool correctBC);
464
465template void getWeightedModes(
466 PtrList<volVectorField>& snapshots, PtrList<volVectorField>& modes,
467 word fieldName, bool podex, bool supex, bool sup, label nmodes,
468 bool correctBC);
469
470template<class Type, template<class> class PatchField, class GeoMesh>
471void getModesSVD(
472 PtrList<GeometricField<Type, PatchField, GeoMesh>>& snapshots,
473 PtrList<GeometricField<Type, PatchField, GeoMesh>>& modes,
474 word fieldName, bool podex, bool supex, bool sup, label nmodes,
475 bool correctBC)
476{
477 ITHACAparameters* para(ITHACAparameters::getInstance());
478
479 if ((podex == 0 && sup == 0) || (supex == 0 && sup == 1))
480 {
481 PtrList<volVectorField> Bases;
482 modes.resize(nmodes);
483 Info << "####### Performing POD using Singular Value Decomposition for " <<
484 snapshots[0].name() << " #######" << endl;
485 Eigen::MatrixXd SnapMatrix = Foam2Eigen::PtrList2Eigen(snapshots);
486 Eigen::VectorXd V = ITHACAutilities::getMassMatrixFV(snapshots[0]);
487 Eigen::VectorXd V3dSqrt = V.array().sqrt();
488 Eigen::VectorXd V3dInv = V3dSqrt.array().cwiseInverse();
489 auto VMsqr = V3dSqrt.asDiagonal();
490 auto VMsqrInv = V3dInv.asDiagonal();
491 Eigen::MatrixXd SnapMatrix2 = VMsqr * SnapMatrix;
492 Eigen::JacobiSVD<Eigen::MatrixXd> svd(SnapMatrix2,
493 Eigen::ComputeThinU | Eigen::ComputeThinV);
494 Info << "####### End of the POD for " << snapshots[0].name() << " #######" <<
495 endl;
496 Eigen::VectorXd eigenValueseig;
497 Eigen::MatrixXd eigenVectoreig;
498 eigenValueseig = svd.singularValues().real();
499 eigenVectoreig = svd.matrixU().real();
500 Eigen::MatrixXd modesEig = VMsqrInv * eigenVectoreig;
501 GeometricField<Type, PatchField, GeoMesh> tmb_bu(snapshots[0].name(),
502 snapshots[0] * 0);
503
504 for (label i = 0; i < nmodes; i++)
505 {
506 Eigen::VectorXd vec = modesEig.col(i);
507 tmb_bu = Foam2Eigen::Eigen2field(tmb_bu, vec, correctBC);
508 modes.set(i, tmb_bu.clone());
509 }
510
511 eigenValueseig = eigenValueseig / eigenValueseig.sum();
512 Eigen::VectorXd cumEigenValues(eigenValueseig);
513
514 for (label j = 1; j < cumEigenValues.size(); ++j)
515 {
516 cumEigenValues(j) += cumEigenValues(j - 1);
517 }
518
519 Info << "####### Saving the POD bases for " << snapshots[0].name() <<
520 " #######" << endl;
521
522 if (sup)
523 {
524 ITHACAstream::exportFields(modes, "./ITHACAoutput/supremizer/",
525 snapshots[0].name());
526 }
527 else
528 {
529 ITHACAstream::exportFields(modes, "./ITHACAoutput/POD/", snapshots[0].name());
530 }
531
532 //exportBases(modes, snapshots, sup);
533 Eigen::saveMarketVector(eigenValueseig,
534 "./ITHACAoutput/POD/Eigenvalues_" + snapshots[0].name(), para->precision,
535 para->outytpe);
536 Eigen::saveMarketVector(cumEigenValues,
537 "./ITHACAoutput/POD/CumEigenvalues_" + snapshots[0].name(), para->precision,
538 para->outytpe);
539 }
540 else
541 {
542 Info << "Reading the existing modes" << endl;
543
544 if (sup == 1)
545 {
546 ITHACAstream::read_fields(modes, fieldName + "sup",
547 "./ITHACAoutput/supremizer/");
548 }
549 else
550 {
551 ITHACAstream::read_fields(modes, fieldName, "./ITHACAoutput/POD/");
552 }
553 }
554}
555
556template void getModesSVD(
557 PtrList<volScalarField>& snapshots, PtrList<volScalarField>& modes,
558 word fieldName, bool podex, bool supex, bool sup, label nmodes,
559 bool correctBC);
560
561template void getModesSVD(
562 PtrList<volVectorField>& snapshots, PtrList<volVectorField>& modes,
563 word fieldName, bool podex, bool supex, bool sup, label nmodes,
564 bool correctBC);
565
567template<>
568Eigen::MatrixXd corMatrix(PtrList<volScalarField>& snapshots)
569{
570 Info << "########## Filling the correlation matrix for " << snapshots[0].name()
571 << "##########" << endl;
572 Eigen::MatrixXd matrix( snapshots.size(), snapshots.size());
573
574 for (label i = 0; i < snapshots.size(); i++)
575 {
576 for (label j = 0; j <= i; j++)
577 {
578 matrix(i, j) = fvc::domainIntegrate(snapshots[i] * snapshots[j]).value();
579 }
580 }
581
582 for (label i = 1; i < snapshots.size(); i++)
583 {
584 for (label j = 0; j < i; j++)
585 {
586 matrix(j, i) = matrix(i, j);
587 }
588 }
589
590 return matrix;
591}
592
593
595template<>
596Eigen::MatrixXd corMatrix(PtrList<volVectorField>& snapshots)
597{
598 Info << "########## Filling the correlation matrix for " << snapshots[0].name()
599 << "##########" << endl;
600 Eigen::MatrixXd matrix( snapshots.size(), snapshots.size());
601
602 for (label i = 0; i < snapshots.size(); i++)
603 {
604 for (label j = 0; j <= i; j++)
605 {
606 matrix(i, j) = fvc::domainIntegrate(snapshots[i] & snapshots[j]).value();
607 }
608 }
609
610 for (label i = 1; i < snapshots.size(); i++)
611 {
612 for (label j = 0; j < i; j++)
613 {
614 matrix(j, i) = matrix(i, j);
615 }
616 }
617
618 return matrix;
619}
620
622template<>
623Eigen::MatrixXd corMatrix(List<Eigen::SparseMatrix<double>>&
624 snapshots)
625{
626 Info << "########## Filling the correlation matrix for the matrix list ##########"
627 << endl;
628 Eigen::MatrixXd matrix( snapshots.size(), snapshots.size());
629
630 for (label i = 0; i < snapshots.size(); i++)
631 {
632 for (label j = 0; j <= i; j++)
633 {
634 double res = 0;
635
636 for (label k = 0; k < snapshots[i].cols(); k++)
637 {
638 res += snapshots[i].col(k).dot(snapshots[j].col(k));
639 }
640
641 matrix(i, j) = res;
642 }
643 }
644
645 for (label i = 1; i < snapshots.size(); i++)
646 {
647 for (label j = 0; j < i; j++)
648 {
649 matrix(j, i) = matrix(i, j);
650 }
651 }
652
653 return matrix;
654}
655
657template<>
658Eigen::MatrixXd corMatrix(List<Eigen::VectorXd>& snapshots)
659{
660 Info << "########## Filling the correlation matrix for the matrix list ##########"
661 << endl;
662 Eigen::MatrixXd matrix( snapshots.size(), snapshots.size());
663
664 for (label i = 0; i < snapshots.size(); i++)
665 {
666 for (label j = 0; j <= i; j++)
667 {
668 matrix(i, j) = (snapshots[i].transpose() * snapshots[j]).trace();
669 }
670 }
671
672 for (label i = 1; i < snapshots.size(); i++)
673 {
674 for (label j = 0; j < i; j++)
675 {
676 matrix(j, i) = matrix(i, j);
677 }
678 }
679
680 return matrix;
681}
682
683
684
686template<class Type, template<class> class PatchField, class GeoMesh>
687void exportBases(PtrList<GeometricField<Type, PatchField, GeoMesh>>& s,
688 PtrList<GeometricField<Type, PatchField, GeoMesh>>& bases,
689 word fieldName, bool sup)
690{
691 if (sup)
692 {
693 fileName fieldname;
694
695 for (label i = 0; i < s.size(); i++)
696 {
697 mkDir("./ITHACAoutput/supremizer/" + name(i + 1));
698 fieldname = "./ITHACAoutput/supremizer/" + name(i + 1) + "/" +
699 bases[i].name();
700 OFstream os(fieldname);
701 bases[i].writeHeader(os);
702 os << s[i] << endl;
703 }
704 }
705 else
706 {
707 fileName fieldname;
708
709 for (label i = 0; i < s.size(); i++)
710 {
711 mkDir("./ITHACAoutput/POD/" + name(i + 1));
712 fieldname = "./ITHACAoutput/POD/" + name(i + 1) + "/" + bases[i].name();
713 OFstream os(fieldname);
714 bases[i].writeHeader(os);
715 os << s[i] << endl;
716 }
717 }
718}
719template void exportBases(PtrList<volVectorField>& s,
720 PtrList<volVectorField>& bases, word fieldName, bool sup);
721template void exportBases(PtrList<volScalarField>& s,
722 PtrList<volScalarField>& bases, word fieldName, bool sup);
723
724void exportEigenvalues(scalarField Eigenvalues, fileName name,
725 bool sup)
726{
727 if (sup)
728 {
729 fileName fieldname;
730 mkDir("./ITHACAoutput/supremizer/");
731 fieldname = "./ITHACAoutput/supremizer/Eigenvalues_" + name;
732 OFstream os(fieldname);
733
734 for (label i = 0; i < Eigenvalues.size(); i++)
735 {
736 os << Eigenvalues[i] << endl;
737 }
738 }
739 else
740 {
741 fileName fieldname;
742 mkDir("./ITHACAoutput/POD/");
743 fieldname = "./ITHACAoutput/POD/Eigenvalues_" + name;
744 OFstream os(fieldname);
745
746 for (label i = 0; i < Eigenvalues.size(); i++)
747 {
748 os << Eigenvalues[i] << endl;
749 }
750 }
751}
752
753void exportcumEigenvalues(scalarField cumEigenvalues, fileName name,
754 bool sup)
755{
756 if (sup)
757 {
758 fileName fieldname;
759 mkDir("./ITHACAoutput/supremizer/");
760 fieldname = "./ITHACAoutput/supremizer/cumEigenvalues_" + name;
761 OFstream os(fieldname);
762
763 for (label i = 0; i < cumEigenvalues.size(); i++)
764 {
765 os << cumEigenvalues[i] << endl;
766 }
767 }
768 else
769 {
770 fileName fieldname;
771 mkDir("./ITHACAoutput/POD/");
772 fieldname = "./ITHACAoutput/POD/cumEigenvalues_" + name;
773 OFstream os(fieldname);
774
775 for (label i = 0; i < cumEigenvalues.size(); i++)
776 {
777 os << cumEigenvalues[i] << endl;
778 }
779 }
780}
781
782
783std::tuple<List<Eigen::SparseMatrix<double>>, List<Eigen::VectorXd>>
784 DEIMmodes(List<Eigen::SparseMatrix<double>>& A,
785 List<Eigen::VectorXd>& b, label nmodesA, label nmodesB, word MatrixName)
786{
787 ITHACAparameters* para(ITHACAparameters::getInstance());
788 List<Eigen::SparseMatrix<double>> ModesA(nmodesA);
789 List<Eigen::VectorXd> ModesB(nmodesB);
790
791 if (!ITHACAutilities::check_folder("./ITHACAoutput/DEIM/" + MatrixName))
792 {
793 if (nmodesA > A.size() - 2 || nmodesB > A.size() - 2 )
794 {
795 std::cout <<
796 "The number of requested modes cannot be bigger than the number of Snapshots - 2"
797 << std::endl;
798 exit(0);
799 }
800
801 scalarField eigenValuesA(nmodesA);
802 scalarField cumEigenValuesA(nmodesA);
803 scalarField eigenValuesB(nmodesB);
804 scalarField cumEigenValuesB(nmodesB);
805 List<scalarField> eigenVectorA(nmodesA);
806 List<scalarField> eigenVectorB(nmodesB);
807
808 for (label i = 0; i < nmodesA; i++)
809 {
810 eigenVectorA[i].setSize(A.size());
811 }
812
813 for (label i = 0; i < nmodesB; i++)
814 {
815 eigenVectorB[i].setSize(A.size());
816 }
817
818 Eigen::MatrixXd corMatrixA = corMatrix(A);
819 Eigen::MatrixXd corMatrixB = corMatrix(b);
820 Info << "####### Performing the POD for the Matrix List #######" << endl;
821 Spectra::DenseSymMatProd<double> opA(corMatrixA);
822 Spectra::DenseSymMatProd<double> opB(corMatrixB);
823 Spectra::SymEigsSolver<double, Spectra::LARGEST_ALGE, Spectra:: DenseSymMatProd<double>>
824 esA(&opA, nmodesA, nmodesA + 10);
825 Spectra::SymEigsSolver<double, Spectra::LARGEST_ALGE, Spectra:: DenseSymMatProd<double>>
826 esB(&opB, nmodesB, nmodesB + 10);
827 esA.init();
828 esB.init();
829 esA.compute();
830
831 if (nmodesB != 1)
832 {
833 esB.compute();
834 }
835
836 Info << "####### End of the POD for the Matrix List #######" << endl;
837 Eigen::VectorXd eigenValueseigA;
838 Eigen::MatrixXd eigenVectorseigA;
839 Eigen::VectorXd eigenValueseigB;
840 Eigen::MatrixXd eigenVectorseigB;
841
842 if (esA.info() == Spectra::SUCCESSFUL)
843 {
844 eigenValueseigA = esA.eigenvalues().real();
845 eigenVectorseigA = esA.eigenvectors().real();
846
847 if (esB.info() == Spectra::SUCCESSFUL && nmodesB != 1)
848 {
849 eigenValueseigB = esB.eigenvalues().real();
850 eigenVectorseigB = esB.eigenvectors().real();
851 }
852 else if (nmodesB == 1)
853 {
854 eigenValueseigB.resize(1);
855 eigenVectorseigB.resize(A.size(), nmodesB);
856 eigenValueseigB(0) = 1;
857 eigenVectorseigB = eigenVectorseigB * 0;
858 eigenVectorseigB(0, 0) = 1;
859 }
860 else
861 {
862 Info << "The Eigenvalue solver in ITHACAPOD.H did not converge, exiting the code"
863 << endl;
864 exit(0);
865 }
866 }
867 else
868 {
869 Info << "The Eigenvalue solver in ITHACAPOD.H did not converge, exiting the code"
870 << endl;
871 exit(0);
872 }
873
874 for (label i = 0; i < nmodesA; i++)
875 {
876 eigenValuesA[i] = eigenValueseigA(i) / eigenValueseigA.sum();
877 }
878
879 for (label i = 0; i < nmodesB; i++)
880 {
881 eigenValuesB[i] = eigenValueseigB(i) / eigenValueseigB.sum();
882 }
883
884 for (label i = 0; i < nmodesA; i++)
885 {
886 for (label k = 0; k < A.size(); k++)
887 {
888 eigenVectorA[i][k] = eigenVectorseigA(k, i);
889 }
890 }
891
892 for (label i = 0; i < nmodesB; i++)
893 {
894 for (label k = 0; k < A.size(); k++)
895 {
896 eigenVectorB[i][k] = eigenVectorseigB(k, i);
897 }
898 }
899
900 cumEigenValuesA[0] = eigenValuesA[0];
901 cumEigenValuesB[0] = eigenValuesB[0];
902
903 for (label i = 1; i < nmodesA; i++)
904 {
905 cumEigenValuesA[i] = cumEigenValuesA[i - 1] + eigenValuesA[i];
906 }
907
908 for (label i = 1; i < nmodesB; i++)
909 {
910 cumEigenValuesB[i] = cumEigenValuesB[i - 1] + eigenValuesB[i];
911 }
912
913 Eigen::SparseMatrix<double> tmp_A;
914 Eigen::VectorXd tmp_B;
915
916 for (label i = 0; i < nmodesA; i++)
917 {
918 tmp_A = eigenVectorA[i][0] * A[0];
919
920 for (label k = 1; k < A.size(); k++)
921 {
922 tmp_A += eigenVectorA[i][k] * A[k];
923 }
924
925 ModesA[i] = tmp_A;
926 }
927
928 for (label i = 0; i < nmodesB; i++)
929 {
930 tmp_B = eigenVectorB[i][0] * b[0];
931
932 for (label k = 1; k < A.size(); k++)
933 {
934 tmp_B += eigenVectorB[i][k] * b[k];
935 }
936
937 ModesB[i] = tmp_B;
938 }
939
940 ITHACAstream::exportList(eigenValuesA,
941 "./ITHACAoutput/DEIM/" + MatrixName + "/", "eigenValuesA_" + MatrixName);
942 ITHACAstream::exportList(cumEigenValuesA,
943 "./ITHACAoutput/DEIM/" + MatrixName + "/", "cumEigenValuesA_" + MatrixName);
944 ITHACAstream::exportList(eigenValuesB,
945 "./ITHACAoutput/DEIM/" + MatrixName + "/", "eigenValuesB_" + MatrixName);
946 ITHACAstream::exportList(cumEigenValuesB,
947 "./ITHACAoutput/DEIM/" + MatrixName + "/", "cumEigenValuesB_" + MatrixName);
948
949 for (label i = 0; i < ModesA.size(); i++)
950 {
951 ITHACAstream::SaveSparseMatrix(ModesA[i],
952 "./ITHACAoutput/DEIM/" + MatrixName + "/", "A_" + MatrixName + name(i));
953 }
954
955 for (label i = 0; i < ModesB.size(); i++)
956 {
957 ITHACAstream::SaveDenseMatrix(ModesB[i],
958 "./ITHACAoutput/DEIM/" + MatrixName + "/", "B_" + MatrixName + name(i));
959 }
960 }
961 else
962 {
963 for (label i = 0; i < nmodesA; i++)
964 {
965 ITHACAstream::ReadSparseMatrix(ModesA[i],
966 "./ITHACAoutput/DEIM/" + MatrixName + "/", "A_" + MatrixName + name(i));
967 }
968
969 for (label i = 0; i < nmodesB; i++)
970 {
971 ITHACAstream::ReadDenseMatrix(ModesB[i],
972 "./ITHACAoutput/DEIM/" + MatrixName + "/", "B_" + MatrixName + name(i));
973 }
974 }
975
976 std::tuple <List<Eigen::SparseMatrix<double>>, List<Eigen::VectorXd>> tupla;
977 tupla = std::make_tuple(ModesA, ModesB);
978 return tupla;
979}
980
981void GrammSchmidt(Eigen::MatrixXd& Matrix)
982{
983 Eigen::MatrixXd Ortho = Matrix;
984 Ortho = Matrix;
985
986 for (label i = 0; i < Matrix.cols(); i++)
987 {
988 for (label k = 0; k < i; k++)
989 {
990 double num = Ortho.col(k).transpose() * Matrix.col(i);
991 double den = (Ortho.col(k).transpose() * Ortho.col(k));
992 double fact = num / den;
993 Ortho.col(i) -= fact * Ortho.col(k) ;
994 }
995
996 Ortho.col(i).normalize();
997 }
998
999 Matrix = Ortho;
1000}
1001template<class Type, template<class> class PatchField, class GeoMesh>
1002void getModes(
1003 PtrList<GeometricField<Type, PatchField, GeoMesh>>& snapshots,
1004 PtrList<GeometricField<Type, PatchField, GeoMesh>>& modes,
1005 PtrList<volScalarField>& Volumes, word fieldName, bool podex, bool supex,
1006 bool sup, label nmodes, bool correctBC)
1007{
1008 ITHACAparameters* para(ITHACAparameters::getInstance());
1009
1010 if (nmodes == 0 && para->eigensolver == "spectra")
1011 {
1012 nmodes = snapshots.size() - 2;
1013 }
1014
1015 if (nmodes == 0 && para->eigensolver == "eigen")
1016 {
1017 nmodes = snapshots.size();
1018 }
1019
1020 if (para->eigensolver == "spectra")
1021 {
1022 M_Assert(nmodes <= snapshots.size() - 2,
1023 "The number of requested modes cannot be bigger than the number of Snapshots - 2");
1024 }
1025
1026 if ((podex == 0 && sup == 0) || (supex == 0 && sup == 1))
1027 {
1028 Eigen::MatrixXd SnapMatrix = Foam2Eigen::PtrList2Eigen(snapshots);
1029 List<Eigen::MatrixXd> SnapMatrixBC = Foam2Eigen::PtrList2EigenBC(snapshots);
1030 label NBC = snapshots[0].boundaryField().size();
1031 Eigen::MatrixXd V = Foam2Eigen::PtrList2Eigen(Volumes);
1032 Eigen::MatrixXd _corMatrix(snapshots.size(), snapshots.size());
1033 Info << "Filling the correlation matrix for field " << snapshots[0].name() <<
1034 endl;
1035
1036 for (label i = 0; i < snapshots.size(); i++)
1037 {
1038 for (label j = 0; j <= i; j++)
1039 {
1040 Eigen::VectorXd Mij = (V.col(i).array() * V.col(j).array());
1041 Mij = Mij.array().abs().sqrt();
1042 _corMatrix(i, j) = SnapMatrix.col(i).transpose() * Mij.asDiagonal() *
1043 SnapMatrix.col(j);
1044 }
1045 }
1046
1047 std::cout << std::endl;
1048
1049 for (label i = 1; i < snapshots.size(); i++)
1050 {
1051 for (label j = 0; j < i; j++)
1052 {
1053 _corMatrix(j, i) = _corMatrix(i, j);
1054 }
1055 }
1056
1057 Eigen::VectorXd eigenValueseig;
1058 Eigen::MatrixXd eigenVectoreig;
1059 modes.resize(nmodes);
1060 Info << "####### Performing the POD using EigenDecomposition for " <<
1061 snapshots[0].name() << " #######" << endl;
1062 label ncv = snapshots.size();
1063 Spectra::DenseSymMatProd<double> op(_corMatrix);
1064 Eigen::SelfAdjointEigenSolver<Eigen::MatrixXd> esEg;
1065
1066 if (para->eigensolver == "spectra")
1067 {
1068 Spectra::SymEigsSolver<double, Spectra::LARGEST_ALGE, Spectra::DenseSymMatProd<double>>
1069 es(&op, nmodes, ncv);
1070 std::cout << "Using Spectra EigenSolver " << std::endl;
1071 es.init();
1072 es.compute(1000, 1e-10, Spectra::LARGEST_ALGE);
1073 M_Assert(es.info() == Spectra::SUCCESSFUL,
1074 "The Eigenvalue Decomposition did not succeed");
1075 eigenVectoreig = es.eigenvectors().real();
1076 eigenValueseig = es.eigenvalues().real();
1077 }
1078 else if (para->eigensolver == "eigen")
1079 {
1080 std::cout << "Using Eigen EigenSolver " << std::endl;
1081 esEg.compute(_corMatrix);
1082 M_Assert(esEg.info() == Eigen::Success,
1083 "The Eigenvalue Decomposition did not succeed");
1084 eigenVectoreig = esEg.eigenvectors().real().rowwise().reverse().leftCols(
1085 nmodes);
1086 eigenValueseig = esEg.eigenvalues().real().reverse().head(nmodes);
1087 }
1088
1089 Info << "####### End of the POD for " << snapshots[0].name() << " #######" <<
1090 endl;
1091 Eigen::VectorXd eigenValueseigLam =
1092 eigenValueseig.real().array().cwiseInverse().sqrt() ;
1093 Eigen::MatrixXd modesEig = (SnapMatrix * eigenVectoreig) *
1094 eigenValueseigLam.asDiagonal();
1095 List<Eigen::MatrixXd> modesEigBC;
1096 modesEigBC.resize(NBC);
1097
1098 for (label i = 0; i < NBC; i++)
1099 {
1100 modesEigBC[i] = (SnapMatrixBC[i] * eigenVectoreig) *
1101 eigenValueseigLam.asDiagonal();
1102 }
1103
1104 for (label i = 0; i < modes.size(); i++)
1105 {
1106 GeometricField<Type, PatchField, GeoMesh> tmp2(snapshots[0].name(),
1107 snapshots[0] * 0);
1108 Eigen::VectorXd vec = modesEig.col(i);
1109 tmp2 = Foam2Eigen::Eigen2field(tmp2, vec, correctBC);
1110
1111 // Adjusting boundary conditions
1112 for (label k = 0; k < tmp2.boundaryField().size(); k++)
1113 {
1114 ITHACAutilities::assignBC(tmp2, k, modesEigBC[k].col(i));
1115 }
1116
1117 modes.set(i, tmp2.clone());
1118 }
1119
1120 eigenValueseig = eigenValueseig / eigenValueseig.sum();
1121 Eigen::VectorXd cumEigenValues(eigenValueseig);
1122
1123 for (label j = 1; j < cumEigenValues.size(); ++j)
1124 {
1125 cumEigenValues(j) += cumEigenValues(j - 1);
1126 }
1127
1128 Info << "####### Saving the POD bases for " << snapshots[0].name() << " #######"
1129 << endl;
1130 exportBases(modes, snapshots, fieldName, sup);
1131 Eigen::saveMarketVector(eigenValueseig,
1132 "./ITHACAoutput/POD/Eigenvalues_" + snapshots[0].name(), para->precision,
1133 para->outytpe);
1134 Eigen::saveMarketVector(cumEigenValues,
1135 "./ITHACAoutput/POD/CumEigenvalues_" + snapshots[0].name(), para->precision,
1136 para->outytpe);
1137 }
1138 else
1139 {
1140 Info << "Reading the existing modes" << endl;
1141
1142 if (sup == 1)
1143 {
1144 ITHACAstream::read_fields(modes, "Usup", "./ITHACAoutput/supremizer/");
1145 }
1146 else
1147 {
1148 ITHACAstream::read_fields(modes, fieldName, "./ITHACAoutput/POD/");
1149 }
1150 }
1151}
1152template void getModes(
1153 PtrList<volScalarField>& snapshots, PtrList<volScalarField>& modes,
1154 PtrList<volScalarField>& Volumes, word fieldName, bool podex, bool supex,
1155 bool sup, label nmodes, bool correctBC);
1156
1157template void getModes(
1158 PtrList<volVectorField>& snapshots, PtrList<volVectorField>& modes,
1159 PtrList<volScalarField>& Volumes, word fieldName, bool podex, bool supex,
1160 bool sup, label nmodes, bool correctBC);
1161
1162template<typename type_matrix>
1163std::tuple<List<Eigen::SparseMatrix<double>>, List<Eigen::VectorXd>>
1164 DEIMmodes(PtrList<type_matrix>& MatrixList, label nmodesA, label nmodesB,
1165 word MatrixName)
1166{
1167 ITHACAparameters* para(ITHACAparameters::getInstance());
1168 List<Eigen::SparseMatrix<double>> ModesA(nmodesA);
1169 List<Eigen::VectorXd> ModesB(nmodesB);
1170
1171 if (!ITHACAutilities::check_folder("./ITHACAoutput/DEIM/" + MatrixName))
1172 {
1173 M_Assert(nmodesA <= MatrixList.size() - 2
1174 && nmodesB <= MatrixList.size() - 2,
1175 "The number of requested modes cannot be bigger than the number of Snapshots - 2");
1176 std::tuple<List<Eigen::SparseMatrix<double>>, List<Eigen::VectorXd>> snapshots =
1177 Foam2Eigen::LFvMatrix2LSM(MatrixList);
1178 Eigen::MatrixXd corMatrixA = corMatrix(std::get<0>(snapshots));
1179 Eigen::MatrixXd corMatrixB = corMatrix(std::get<1>(snapshots));
1180 Eigen::VectorXd eigenValueseigA;
1181 Eigen::MatrixXd eigenVectorseigA;
1182 Eigen::VectorXd eigenValueseigB;
1183 Eigen::MatrixXd eigenVectorseigB;
1184
1185 if (para->eigensolver == "spectra")
1186 {
1187 Info << "####### Performing the POD decomposition for the Matrix List using Spectra #######"
1188 << endl;
1189 Spectra::DenseSymMatProd<double> opA(corMatrixA);
1190 Spectra::DenseSymMatProd<double> opB(corMatrixB);
1191 label ncvA = MatrixList.size();
1192 label ncvB = MatrixList.size();
1193 Spectra::SymEigsSolver<double, Spectra::LARGEST_ALGE, Spectra:: DenseSymMatProd<double>>
1194 esA(&opA, nmodesA, ncvA);
1195 Spectra::SymEigsSolver<double, Spectra::LARGEST_ALGE, Spectra:: DenseSymMatProd<double>>
1196 esB(&opB, nmodesB, ncvB);
1197 esA.init();
1198 esB.init();
1199 esA.compute(1000, 1e-10, Spectra::LARGEST_ALGE);
1200 M_Assert(esA.info() == Spectra::SUCCESSFUL,
1201 "The Eigenvalue Decomposition did not succeed");
1202
1203 if (nmodesB != 1)
1204 {
1205 esB.compute(1000, 1e-10, Spectra::LARGEST_ALGE);
1206 M_Assert(esB.info() == Spectra::SUCCESSFUL,
1207 "The Eigenvalue Decomposition did not succeed");
1208 }
1209
1210 Info << "####### End of the POD decomposition for the Matrix List #######" <<
1211 endl;
1212 eigenValueseigA = esA.eigenvalues().real();
1213 eigenVectorseigA = esA.eigenvectors().real();
1214
1215 if (nmodesB != 1)
1216 {
1217 eigenValueseigB = esB.eigenvalues().real();
1218 eigenVectorseigB = esB.eigenvectors().real();
1219 }
1220 else
1221 {
1222 eigenValueseigB.resize(1);
1223 eigenVectorseigB.resize(MatrixList.size(), nmodesB);
1224 eigenValueseigB(0) = 1;
1225 eigenVectorseigB = eigenVectorseigB * 0;
1226 eigenVectorseigB(0, 0) = 1;
1227 }
1228 }
1229 else if (para->eigensolver == "eigen")
1230 {
1231 Info << "####### Performing the POD decomposition for the Matrix List using Eigen #######"
1232 << endl;
1233 Eigen::SelfAdjointEigenSolver<Eigen::MatrixXd> esEgA;
1234 esEgA.compute(corMatrixA);
1235 Eigen::SelfAdjointEigenSolver<Eigen::MatrixXd> esEgB;
1236 esEgB.compute(corMatrixB);
1237 M_Assert(esEgA.info() == Eigen::Success,
1238 "The Eigenvalue Decomposition did not succeed");
1239 M_Assert(esEgB.info() == Eigen::Success,
1240 "The Eigenvalue Decomposition did not succeed");
1241 eigenVectorseigA = esEgA.eigenvectors().real().rowwise().reverse().leftCols(
1242 nmodesA);
1243 eigenValueseigA = esEgA.eigenvalues().real().reverse().head(nmodesA);
1244 eigenVectorseigB = esEgB.eigenvectors().real().rowwise().reverse().leftCols(
1245 nmodesB);
1246 eigenValueseigB = esEgB.eigenvalues().real().reverse().head(nmodesB);
1247 }
1248
1249 Eigen::SparseMatrix<double> tmp_A;
1250 Eigen::VectorXd tmp_B;
1251
1252 for (label i = 0; i < nmodesA; i++)
1253 {
1254 tmp_A = eigenVectorseigA(0, i) * std::get<0>(snapshots)[0];
1255
1256 for (label k = 1; k < MatrixList.size(); k++)
1257 {
1258 tmp_A += eigenVectorseigA(k, i) * std::get<0>(snapshots)[k];
1259 }
1260
1261 ModesA[i] = tmp_A;
1262 }
1263
1264 for (label i = 0; i < nmodesB; i++)
1265 {
1266 tmp_B = eigenVectorseigB(0, i) * std::get<1>(snapshots)[0];
1267
1268 for (label k = 1; k < MatrixList.size(); k++)
1269 {
1270 tmp_B += eigenVectorseigB(k, i) * std::get<1>(snapshots)[k];
1271 }
1272
1273 ModesB[i] = tmp_B;
1274 }
1275
1276 eigenValueseigA = eigenValueseigA / eigenValueseigA.sum();
1277 eigenValueseigB = eigenValueseigB / eigenValueseigB.sum();
1278 Eigen::VectorXd cumEigenValuesA(eigenValueseigA);
1279 Eigen::VectorXd cumEigenValuesB(eigenValueseigB);
1280
1281 for (label j = 1; j < cumEigenValuesA.size(); ++j)
1282 {
1283 cumEigenValuesA(j) += cumEigenValuesA(j - 1);
1284 }
1285
1286 for (label j = 1; j < cumEigenValuesB.size(); ++j)
1287 {
1288 cumEigenValuesB(j) += cumEigenValuesB(j - 1);
1289 }
1290
1291 for (label i = 0; i < ModesA.size(); i++)
1292 {
1293 ITHACAstream::SaveSparseMatrix(ModesA[i],
1294 "./ITHACAoutput/DEIM/" + MatrixName + "/", "A_" + MatrixName + name(i));
1295 }
1296
1297 for (label i = 0; i < ModesB.size(); i++)
1298 {
1299 ITHACAstream::SaveDenseMatrix(ModesB[i],
1300 "./ITHACAoutput/DEIM/" + MatrixName + "/", "B_" + MatrixName + name(i));
1301 }
1302
1303 Eigen::saveMarketVector(eigenValueseigA,
1304 "./ITHACAoutput/DEIM/" + MatrixName + "/eigenValuesA", para->precision,
1305 para->outytpe);
1306 Eigen::saveMarketVector(eigenValueseigB,
1307 "./ITHACAoutput/DEIM/" + MatrixName + "/eigenValuesB", para->precision,
1308 para->outytpe);
1309 Eigen::saveMarketVector(cumEigenValuesA,
1310 "./ITHACAoutput/DEIM/" + MatrixName + "/cumEigenValuesA", para->precision,
1311 para->outytpe);
1312 Eigen::saveMarketVector(cumEigenValuesB,
1313 "./ITHACAoutput/DEIM/" + MatrixName + "/cumEigenValuesB", para->precision,
1314 para->outytpe);
1315 }
1316 else
1317 {
1318 for (label i = 0; i < nmodesA; i++)
1319 {
1320 ITHACAstream::ReadSparseMatrix(ModesA[i],
1321 "./ITHACAoutput/DEIM/" + MatrixName + "/", "A_" + MatrixName + name(i));
1322 }
1323
1324 for (label i = 0; i < nmodesB; i++)
1325 {
1326 ITHACAstream::ReadDenseMatrix(ModesB[i],
1327 "./ITHACAoutput/DEIM/" + MatrixName + "/", "B_" + MatrixName + name(i));
1328 }
1329 }
1330
1331 std::tuple <List<Eigen::SparseMatrix<double>>, List<Eigen::VectorXd>> tupla;
1332 tupla = std::make_tuple(ModesA, ModesB);
1333 return tupla;
1334}
1335
1336template std::tuple<List<Eigen::SparseMatrix<double>>, List<Eigen::VectorXd>>
1337DEIMmodes(PtrList<fvScalarMatrix>& MatrixList, label nmodesA,
1338 label nmodesB,
1339 word MatrixName);
1340
1341template std::tuple<List<Eigen::SparseMatrix<double>>, List<Eigen::VectorXd>>
1342DEIMmodes(PtrList<fvVectorMatrix>& MatrixList, label nmodesA,
1343 label nmodesB,
1344 word MatrixName);
1345
1346template<class Type, template<class> class PatchField, class GeoMesh>
1347PtrList<GeometricField<Type, PatchField, GeoMesh>>DEIMmodes(
1348 PtrList<GeometricField<Type, PatchField, GeoMesh>>& snapshots, label nmodes,
1349 word FunctionName, word fieldName)
1350{
1351 ITHACAparameters* para(ITHACAparameters::getInstance());
1352 word PODkey = "POD_" + fieldName;
1353 word PODnorm = para->ITHACAdict->lookupOrDefault<word>(PODkey, "L2");
1354 M_Assert(PODnorm == "L2" ||
1355 PODnorm == "Frobenius", "The PODnorm can be only L2 or Frobenius");
1356 Info << "Performing POD for " << fieldName << " using the " << PODnorm <<
1357 " norm" << endl;
1358 PtrList<GeometricField<Type, fvPatchField, volMesh>> modes;
1359 bool correctBC = true;
1360
1361 if (nmodes == 0 && para->eigensolver == "spectra")
1362 {
1363 nmodes = snapshots.size() - 2;
1364 }
1365
1366 if (nmodes == 0 && para->eigensolver == "eigen")
1367 {
1368 nmodes = snapshots.size();
1369 }
1370
1371 if (para->eigensolver == "spectra")
1372 {
1373 M_Assert(nmodes <= snapshots.size() - 2,
1374 "The number of requested modes cannot be bigger than the number of Snapshots - 2");
1375 }
1376
1377 if (!ITHACAutilities::check_folder("./ITHACAoutput/DEIM/" + FunctionName))
1378 {
1379 Eigen::MatrixXd SnapMatrix = Foam2Eigen::PtrList2Eigen(snapshots);
1380 List<Eigen::MatrixXd> SnapMatrixBC = Foam2Eigen::PtrList2EigenBC(snapshots);
1381 label NBC = snapshots[0].boundaryField().size();
1382 Eigen::MatrixXd _corMatrix;
1383
1384 if (PODnorm == "L2")
1385 {
1386 _corMatrix = ITHACAutilities::getMassMatrix(snapshots);
1387 }
1388 else if (PODnorm == "Frobenius")
1389 {
1390 _corMatrix = ITHACAutilities::getMassMatrix(snapshots, 0, false);
1391 }
1392
1393 if (Pstream::parRun())
1394 {
1395 reduce(_corMatrix, sumOp<Eigen::MatrixXd>());
1396 }
1397
1398 Eigen::VectorXd eigenValueseig;
1399 Eigen::MatrixXd eigenVectoreig;
1400 modes.resize(nmodes);
1401 Info << "####### Performing the POD using EigenDecomposition " <<
1402 fieldName << " #######" << endl;
1403 label ncv = snapshots.size();
1404 Spectra::DenseSymMatProd<double> op(_corMatrix);
1405 Eigen::SelfAdjointEigenSolver<Eigen::MatrixXd> esEg;
1406
1407 if (para->eigensolver == "spectra")
1408 {
1409 Spectra::SymEigsSolver<double, Spectra::LARGEST_ALGE, Spectra::DenseSymMatProd<double>>
1410 es(&op, nmodes, ncv);
1411 std::cout << "Using Spectra EigenSolver " << std::endl;
1412 es.init();
1413 es.compute(1000, 1e-10, Spectra::LARGEST_ALGE);
1414 M_Assert(es.info() == Spectra::SUCCESSFUL,
1415 "The Eigenvalue Decomposition did not succeed");
1416 eigenVectoreig = es.eigenvectors().real();
1417 eigenValueseig = es.eigenvalues().real();
1418 }
1419 else if (para->eigensolver == "eigen")
1420 {
1421 std::cout << "Using Eigen EigenSolver " << std::endl;
1422 esEg.compute(_corMatrix);
1423 M_Assert(esEg.info() == Eigen::Success,
1424 "The Eigenvalue Decomposition did not succeed");
1425 eigenVectoreig = esEg.eigenvectors().real().rowwise().reverse().leftCols(
1426 nmodes);
1427 eigenValueseig = esEg.eigenvalues().real().array().reverse();
1428 }
1429
1430 if (eigenValueseig.array().minCoeff() < 0)
1431 {
1432 eigenValueseig = eigenValueseig.array() + 2 * abs(
1433 eigenValueseig.array().minCoeff());
1434 }
1435
1436 Info << "####### End of the POD for " << snapshots[0].name() << " #######" <<
1437 endl;
1438 Eigen::MatrixXd modesEig = (SnapMatrix * eigenVectoreig);
1439 Eigen::VectorXd V = ITHACAutilities::getMassMatrixFV(snapshots[0]);
1440 Eigen::MatrixXd normFact(nmodes, 1);
1441
1442 for (label i = 0; i < nmodes; i++)
1443 {
1444 if (PODnorm == "L2")
1445 {
1446 normFact(i, 0) = std::sqrt((modesEig.col(i).transpose() * V.asDiagonal() *
1447 modesEig.col(i))(0, 0));
1448
1449 if (Pstream::parRun())
1450 {
1451 normFact(i, 0) = (modesEig.col(i).transpose() * V.asDiagonal() *
1452 modesEig.col(i))(0, 0);
1453 }
1454 }
1455 else if (PODnorm == "Frobenius")
1456 {
1457 normFact(i, 0) = std::sqrt((modesEig.col(i).transpose() * modesEig.col(i))(0,
1458 0));
1459
1460 if (Pstream::parRun())
1461 {
1462 normFact(i, 0) = (modesEig.col(i).transpose() * modesEig.col(i))(0, 0);
1463 }
1464 }
1465 }
1466
1467 if (Pstream::parRun())
1468 {
1469 reduce(normFact, sumOp<Eigen::MatrixXd>());
1470 }
1471
1472 if (Pstream::parRun())
1473 {
1474 normFact = normFact.cwiseSqrt();
1475 }
1476
1477 List<Eigen::MatrixXd> modesEigBC;
1478 modesEigBC.resize(NBC);
1479
1480 for (label i = 0; i < NBC; i++)
1481 {
1482 modesEigBC[i] = (SnapMatrixBC[i] * eigenVectoreig);
1483 }
1484
1485 std::cout << normFact << std::endl;
1486
1487 for (label i = 0; i < nmodes; i++)
1488 {
1489 modesEig.col(i) = modesEig.col(i).array() / normFact(i, 0);
1490
1491 for (label j = 0; j < NBC; j++)
1492 {
1493 modesEigBC[j].col(i) = modesEigBC[j].col(i).array() / normFact(i, 0);
1494 }
1495 }
1496
1497 for (label i = 0; i < modes.size(); i++)
1498 {
1499 GeometricField<Type, PatchField, GeoMesh> tmp2(snapshots[0].name(),
1500 snapshots[0]);
1501 Eigen::VectorXd vec = modesEig.col(i);
1502 tmp2 = Foam2Eigen::Eigen2field(tmp2, vec, correctBC);
1503
1504 for (label k = 0; k < NBC; k++)
1505 {
1506 ITHACAutilities::assignBC(tmp2, k, modesEigBC[k].col(i));
1507 }
1508
1509 modes.set(i, tmp2.clone());
1510 }
1511
1512 eigenValueseig = eigenValueseig / eigenValueseig.sum();
1513 Eigen::VectorXd cumEigenValues(eigenValueseig);
1514
1515 for (label j = 1; j < cumEigenValues.size(); ++j)
1516 {
1517 cumEigenValues(j) += cumEigenValues(j - 1);
1518 }
1519
1520 Info << "####### Saving the POD bases for " << snapshots[0].name() <<
1521 " #######" << endl;
1522 ITHACAutilities::createSymLink("./ITHACAoutput/DEIM");
1523
1524 for (label i = 0; i < modes.size(); i++)
1525 {
1526 ITHACAstream::exportSolution(modes[i], name(i + 1), "./ITHACAoutput/DEIM",
1527 fieldName);
1528 }
1529
1530 Eigen::saveMarketVector(eigenValueseig,
1531 "./ITHACAoutput/DEIM/eigenValues_" + fieldName, para->precision,
1532 para->outytpe);
1533 Eigen::saveMarketVector(cumEigenValues,
1534 "./ITHACAoutput/DEIM/cumEigenValues_" + fieldName, para->precision,
1535 para->outytpe);
1536 }
1537 else
1538 {
1539 Info << "Reading the existing modes" << endl;
1540 ITHACAstream::read_fields(modes, fieldName, "./ITHACAoutput/DEIM/");
1541 }
1542
1543 return modes;
1544}
1545
1546template<class Field_type, class Field_type_2>
1547void getModes(
1548 PtrList<Field_type>& snapshots, PtrList<Field_type>& modes,
1549 PtrList<Field_type_2>& fields2, word fieldName, bool podex, bool supex,
1550 bool sup, label nmodes, bool correctBC)
1551{
1552 ITHACAparameters* para(ITHACAparameters::getInstance());
1553
1554 if ((podex == 0 && sup == 0) || (supex == 0 && sup == 1))
1555 {
1556 if (para->eigensolver == "spectra" )
1557 {
1558 if (nmodes == 0)
1559 {
1560 nmodes = snapshots.size() - 2;
1561 }
1562
1563 M_Assert(nmodes <= snapshots.size() - 2,
1564 "The number of requested modes cannot be bigger than the number of Snapshots - 2");
1565 }
1566 else
1567 {
1568 if (nmodes == 0)
1569 {
1570 nmodes = snapshots.size();
1571 }
1572
1573 M_Assert(nmodes <= snapshots.size(),
1574 "The number of requested modes cannot be bigger than the number of Snapshots");
1575 }
1576
1577 Eigen::MatrixXd SnapMatrix = Foam2Eigen::PtrList2Eigen(fields2);
1578 List<Eigen::MatrixXd> SnapMatrixBC = Foam2Eigen::PtrList2EigenBC(fields2);
1579 label NBC = fields2[0].boundaryField().size();
1580 auto VM = ITHACAutilities::getMassMatrixFV(fields2[0]);
1581 Eigen::MatrixXd _corMatrix = SnapMatrix.transpose() * VM.asDiagonal() *
1582 SnapMatrix;
1583
1584 if (Pstream::parRun())
1585 {
1586 reduce(_corMatrix, sumOp<Eigen::MatrixXd>());
1587 }
1588
1589 Eigen::VectorXd eigenValueseig;
1590 Eigen::MatrixXd eigenVectoreig;
1591 modes.resize(nmodes);
1592 Info << "####### Performing the POD using EigenDecomposition " <<
1593 fields2[0].name() << " #######" << endl;
1594 label ncv = snapshots.size();
1595 Spectra::DenseSymMatProd<double> op(_corMatrix);
1596 Eigen::SelfAdjointEigenSolver<Eigen::MatrixXd> esEg;
1597
1598 if (para->eigensolver == "spectra")
1599 {
1600 Spectra::SymEigsSolver<double, Spectra::LARGEST_ALGE, Spectra::DenseSymMatProd<double>>
1601 es(&op, nmodes, ncv);
1602 std::cout << "Using Spectra EigenSolver " << std::endl;
1603 es.init();
1604 es.compute(1000, 1e-10, Spectra::LARGEST_ALGE);
1605 M_Assert(es.info() == Spectra::SUCCESSFUL,
1606 "The Eigenvalue Decomposition did not succeed");
1607 eigenVectoreig = es.eigenvectors().real();
1608 eigenValueseig = es.eigenvalues().real();
1609 }
1610 else if (para->eigensolver == "eigen")
1611 {
1612 std::cout << "Using Eigen EigenSolver " << std::endl;
1613 esEg.compute(_corMatrix);
1614 M_Assert(esEg.info() == Eigen::Success,
1615 "The Eigenvalue Decomposition did not succeed");
1616 eigenVectoreig = esEg.eigenvectors().real().rowwise().reverse().leftCols(
1617 nmodes);
1618 eigenValueseig = esEg.eigenvalues().real().reverse().head(nmodes);
1619 }
1620
1621 Info << "####### End of the POD for " << snapshots[0].name() << " #######" <<
1622 endl;
1623 Eigen::VectorXd eigenValueseigLam =
1624 eigenValueseig.real().array().cwiseInverse().abs().sqrt() ;
1625 Eigen::MatrixXd modesEig = (SnapMatrix * eigenVectoreig) *
1626 eigenValueseigLam.asDiagonal();
1627 List<Eigen::MatrixXd> modesEigBC;
1628 modesEigBC.resize(NBC);
1629
1630 for (label i = 0; i < modes.size(); i++)
1631 {
1632 Field_type tmp2(snapshots[0].name(), snapshots[0] * 0);
1633
1634 for (label j = 0; j < snapshots.size(); j++)
1635 {
1636 tmp2 += eigenValueseigLam(i) * eigenVectoreig.col(i)(j) * snapshots[j];
1637 }
1638
1639 modes.set(i, tmp2.clone());
1640 }
1641
1642 eigenValueseig = eigenValueseig / eigenValueseig.sum();
1643 Eigen::VectorXd cumEigenValues(eigenValueseig);
1644
1645 for (label j = 1; j < cumEigenValues.size(); ++j)
1646 {
1647 cumEigenValues(j) += cumEigenValues(j - 1);
1648 }
1649
1650 Info << "####### Saving the POD bases for " << snapshots[0].name() <<
1651 " #######" << endl;
1652
1653 if (sup)
1654 {
1655 ITHACAstream::exportFields(modes, "./ITHACAoutput/supremizer/",
1656 snapshots[0].name());
1657 }
1658 else
1659 {
1660 ITHACAstream::exportFields(modes, "./ITHACAoutput/POD/", snapshots[0].name());
1661 }
1662
1663 Eigen::saveMarketVector(eigenValueseig,
1664 "./ITHACAoutput/POD/Eigenvalues_" + snapshots[0].name(), para->precision,
1665 para->outytpe);
1666 Eigen::saveMarketVector(cumEigenValues,
1667 "./ITHACAoutput/POD/CumEigenvalues_" + snapshots[0].name(), para->precision,
1668 para->outytpe);
1669 }
1670 else
1671 {
1672 Info << "Reading the existing modes" << endl;
1673
1674 if (sup == 1)
1675 {
1676 ITHACAstream::read_fields(modes, fieldName, "./ITHACAoutput/supremizer/");
1677 }
1678 else
1679 {
1680 ITHACAstream::read_fields(modes, fieldName, "./ITHACAoutput/POD/");
1681 }
1682 }
1683}
1684
1685template void getModes(
1686 PtrList<surfaceScalarField>& snapshots, PtrList<surfaceScalarField>& modes,
1687 PtrList<volVectorField>& fields2, word fieldName, bool podex, bool supex,
1688 bool sup, label nmodes, bool correctBC);
1689
1690template void getModes(
1691 PtrList<volScalarField>& snapshots, PtrList<volScalarField>& modes,
1692 PtrList<volVectorField>& fields2, word fieldName, bool podex, bool supex,
1693 bool sup, label nmodes, bool correctBC);
1694
1695template PtrList<volScalarField>
1696DEIMmodes(
1697 PtrList<volScalarField>& SnapShotsMatrix,
1698 label nmodes,
1699 word FunctionName, word FieldName);
1700
1701template PtrList<volVectorField>
1702DEIMmodes(
1703 PtrList<volVectorField>& SnapShotsMatrix,
1704 label nmodes,
1705 word FunctionName, word FieldName);
1706
1707}
EigenFunctions.H
Header file of the EigenFunctions class.
ITHACAPOD.H
Header file of the ITHACAPOD class.
M_Assert
#define M_Assert(Expr, Msg)
Definition ITHACAassert.H:46
Foam2Eigen::Eigen2field
static GeometricField< scalar, PatchField, GeoMesh > Eigen2field(GeometricField< scalar, PatchField, GeoMesh > &field, Eigen::VectorXd &eigen_vector, bool correctBC=true)
Convert a vector in Eigen format into an OpenFOAM scalar GeometricField.
Definition Foam2Eigen.C:506
Foam2Eigen::LFvMatrix2LSM
static std::tuple< List< Eigen::SparseMatrix< double > >, List< Eigen::VectorXd > > LFvMatrix2LSM(PtrList< fvMatrix< Type > > &MatrixList)
Convert a PtrList of OpenFOAM fvMatrix into a tuple of lists of Eigen Sparse Matrices and source vect...
Definition Foam2Eigen.C:1184
Foam2Eigen::PtrList2EigenBC
static List< Eigen::MatrixXd > PtrList2EigenBC(PtrList< GeometricField< scalar, PatchField, GeoMesh > > &fields, label Nfields=-1)
Convert an OpenFOAM scalar field to a List of Eigen Vectors, one for each boundary.
Definition Foam2Eigen.C:257
Foam2Eigen::PtrList2Eigen
static Eigen::MatrixXd PtrList2Eigen(PtrList< GeometricField< Type, PatchField, GeoMesh > > &fields, label Nfields=-1)
Convert a PtrList of snapshots to Eigen matrix (only internal field)
Definition Foam2Eigen.C:638
ITHACAparameters
Class for the definition of some general parameters, the parameters must be defined from the file ITH...
Definition ITHACAparameters.H:41
ITHACAparameters::eigensolver
word eigensolver
type of eigensolver used in the eigenvalue decomposition can be either be eigen or spectra
Definition ITHACAparameters.H:47
ITHACAparameters::precision
label precision
precision of the output Market Matrix objects (i.e. reduced matrices, eigenvalues,...
Definition ITHACAparameters.H:50
ITHACAparameters::ITHACAdict
IOdictionary * ITHACAdict
Dictionary for input objects from file.
Definition ITHACAparameters.H:63
ITHACAparameters::getInstance
static ITHACAparameters * getInstance()
Gets an instance of ITHACAparameters, to be used if the instance is already existing.
Definition ITHACAparameters.C:54
ITHACAparameters::outytpe
std::_Ios_Fmtflags outytpe
type of output format can be fixed or scientific
Definition ITHACAparameters.H:60
A
volScalarField & A
Definition createConstants.H:144
Eigen::saveMarketVector
bool saveMarketVector(const VectorType &vec, const std::string &filename, label prec, std::_Ios_Fmtflags outytpe=std::ios_base::scientific)
Definition EigenFunctions.H:441
ITHACAPOD
namespace for the computation of the POD, it exploits bot the method of snapshots and SVD.
Definition incrementalPOD.H:183
ITHACAPOD::exportcumEigenvalues
void exportcumEigenvalues(scalarField cumEigenvalues, fileName name, bool sup)
Export the cumulative eigenvalues as a txt file.
Definition ITHACAPOD.C:753
ITHACAPOD::exportEigenvalues
void exportEigenvalues(scalarField Eigenvalues, fileName name, bool sup)
Exports the eigenvalues as a txt file.
Definition ITHACAPOD.C:724
ITHACAPOD::getModes
void getModes(PtrList< GeometricField< Type, PatchField, GeoMesh > > &snapshots, PtrList< GeometricField< Type, PatchField, GeoMesh > > &modes, word fieldName, bool podex, bool supex, bool sup, label nmodes, bool correctBC)
Computes the bases or reads them for a field.
Definition ITHACAPOD.C:93
ITHACAPOD::DEIMmodes
std::tuple< List< Eigen::SparseMatrix< double > >, List< Eigen::VectorXd > > DEIMmodes(List< Eigen::SparseMatrix< double > > &A, List< Eigen::VectorXd > &b, label nmodesA, label nmodesB, word MatrixName)
Get the DEIM modes for a generic non a parametrized matrix coming from a differential operator functi...
Definition ITHACAPOD.C:784
ITHACAPOD::getNestedSnapshotMatrix
void getNestedSnapshotMatrix(PtrList< GeometricField< Type, PatchField, GeoMesh > > &snapshots, PtrList< GeometricField< Type, PatchField, GeoMesh > > &ModesGlobal, word fieldName, label Npar, label NnestedOut)
Nested-POD approach.
Definition ITHACAPOD.C:41
ITHACAPOD::exportBases
void exportBases(PtrList< GeometricField< Type, PatchField, GeoMesh > > &s, PtrList< GeometricField< Type, PatchField, GeoMesh > > &bases, word fieldName, bool sup)
Export the Bases.
Definition ITHACAPOD.C:687
ITHACAPOD::getModesSVD
void getModesSVD(PtrList< GeometricField< Type, PatchField, GeoMesh > > &snapshots, PtrList< GeometricField< Type, PatchField, GeoMesh > > &modes, word fieldName, bool podex, bool supex, bool sup, label nmodes, bool correctBC)
Gets the bases for a scalar field using SVD instead of the method of snapshots.
Definition ITHACAPOD.C:471
ITHACAPOD::corMatrix
Eigen::MatrixXd corMatrix(PtrList< volScalarField > &snapshots)
Construct the Correlation Matrix for Scalar Field.
Definition ITHACAPOD.C:568
ITHACAPOD::GrammSchmidt
void GrammSchmidt(Eigen::MatrixXd &Matrix)
Performs GrammSchmidt orthonormalization on an Eigen Matrix.
Definition ITHACAPOD.C:981
ITHACAPOD::getWeightedModes
void getWeightedModes(PtrList< GeometricField< Type, PatchField, GeoMesh > > &snapshots, PtrList< GeometricField< Type, PatchField, GeoMesh > > &modes, word fieldName, bool podex, bool supex, bool sup, label nmodes, bool correctBC)
Computes the weighted bases (using the nested-pod approach) or read them for a vector field.
Definition ITHACAPOD.C:328
ITHACAstream::SaveSparseMatrix
void SaveSparseMatrix(Eigen::SparseMatrix< T, Nrows, IND > &m, word folder, word MatrixName)
Export an Eigen sparse matrix into bynary format file.
ITHACAstream::ReadDenseMatrix
void ReadDenseMatrix(MatrixType &Matrix, word folder, word MatrixName)
Read a dense matrix from a binary format file.
Definition ITHACAstream.H:632
ITHACAstream::exportFields
void exportFields(PtrList< GeometricField< Type, PatchField, GeoMesh > > &field, word folder, word fieldname)
Function to export a scalar of vector field.
Definition ITHACAstream.C:715
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:746
ITHACAstream::exportList
void exportList(T &list, word folder, word filename)
Export a list to file.
Definition ITHACAstream.C:956
ITHACAstream::SaveDenseMatrix
void SaveDenseMatrix(MatrixType &Matrix, word folder, word MatrixName)
Save a dense matrix to a binary format file.
Definition ITHACAstream.H:617
ITHACAstream::read_fields
void read_fields(PtrList< GeometricField< Type, PatchField, GeoMesh > > &Lfield, word Name, fileName casename, int first_snap, int n_snap)
Function to read a list of fields from the name of the field and casename.
Definition ITHACAstream.C:454
ITHACAstream::ReadSparseMatrix
void ReadSparseMatrix(Eigen::SparseMatrix< T, Nrows, IND > &m, word folder, word MatrixName)
Read an Eigen sparse matrix from a bynary format file.
ITHACAutilities::getMassMatrixFV
Eigen::VectorXd getMassMatrixFV(GeometricField< Type, PatchField, GeoMesh > &snapshot)
Gets a vector containing the volumes of each cell of the mesh.
Definition ITHACAcoeffsMass.C:120
ITHACAutilities::getMassMatrix
Eigen::MatrixXd getMassMatrix(PtrList< GeometricField< Type, PatchField, GeoMesh > > &modes, label Nmodes, bool consider_volumes)
Gets the mass matrix using the eigen routine.
Definition ITHACAcoeffsMass.C:68
ITHACAutilities::assignBC
void assignBC(GeometricField< scalar, fvPatchField, volMesh > &s, label BC_ind, double value)
Assign uniform Boundary Condition to a volScalarField.
Definition ITHACAassign.C:143
ITHACAutilities::check_folder
bool check_folder(word folder)
Checks if a folder exists.
Definition ITHACAsystem.C:89
ITHACAutilities::createSymLink
void createSymLink(word folder)
Creates symbolic links to 0, system and constant.
Definition ITHACAsystem.C:30
col
singVal col(i)
svd
Info<< endl;Info<< "*********************************************************"<< endl;Info<< "Performing test for the parameterized BC inverse solver"<< endl;Info<< endl;word outputFolder="./ITHACAoutput/parameterizedBCtest_RBFparameter/";volScalarField gTrueField=example_paramBC.list2Field(example_paramBC.gTrue);ITHACAstream::exportSolution(gTrueField, "1", outputFolder, "gTrue");Eigen::VectorXd rbfWidth=EigenFunctions::ExpSpaced(0.001, 100, rbfWidthTest_size);ITHACAstream::exportMatrix(rbfWidth, "rbfShapeParameters", "eigen", outputFolder);Eigen::VectorXd residualNorms;residualNorms.resize(rbfWidthTest_size);Eigen::VectorXd heatFluxL2norm(rbfWidthTest_size);Eigen::VectorXd heatFluxLinfNorm=heatFluxL2norm;Eigen::VectorXd condNumber=heatFluxL2norm;Eigen::MatrixXd singVal;for(int i=0;i< rbfWidthTest_size;i++){ Info<< "*********************************************************"<< endl;Info<< "RBF parameter "<< rbfWidth(i)<< endl;Info<< "Test "<< i+1<< endl;Info<< endl;example_paramBC.set_gParametrized("rbf", rbfWidth(i));example_paramBC.parameterizedBCoffline(1);example_paramBC.parameterizedBC("fullPivLU");volScalarField gParametrizedField=example_paramBC.list2Field(example_paramBC.g);ITHACAstream::exportSolution(gParametrizedField, std::to_string(i+1), outputFolder, "gParametrized");volScalarField &T(example_paramBC._T());ITHACAstream::exportSolution(T, std::to_string(i+1), outputFolder, "T");residualNorms(i)=Foam::sqrt(example_paramBC.residual.squaredNorm());Eigen::MatrixXd A=example_paramBC.Theta.transpose() *example_paramBC.Theta;Eigen::JacobiSVD< Eigen::MatrixXd > svd(A, Eigen::ComputeThinU|Eigen::ComputeThinV)
i
label i
Definition pEqn.H:46
Generated by   doxygen 1.10.0