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 constexpr bool check_vol = std::is_same<volMesh, GeoMesh>::value
101 || std::is_same<surfaceMesh, GeoMesh>::value;
102 word PODkey = "POD_" + fieldName;
103 word PODnorm = para->ITHACAdict->lookupOrDefault<word>(PODkey, "L2");
104 M_Assert(PODnorm == "L2" ||
105 PODnorm == "Frobenius", "The PODnorm can be only L2 or Frobenius");
106 Info << "Performing POD for " << fieldName << " using the " << PODnorm <<
107 " norm" << endl;
108
109 if ((podex == 0 && sup == 0) || (supex == 0 && sup == 1))
110 {
111 if (para->eigensolver == "spectra" )
112 {
113 if (nmodes == 0)
114 {
115 nmodes = snapshots.size() - 2;
116 }
117
118 M_Assert(nmodes <= snapshots.size() - 2,
119 "The number of requested modes cannot be bigger than the number of Snapshots - 2");
120 }
121 else
122 {
123 if (nmodes == 0)
124 {
125 nmodes = snapshots.size();
126 }
127
128 M_Assert(nmodes <= snapshots.size(),
129 "The number of requested modes cannot be bigger than the number of Snapshots");
130 }
131
132 Eigen::MatrixXd SnapMatrix = Foam2Eigen::PtrList2Eigen(snapshots);
133 List<Eigen::MatrixXd> SnapMatrixBC = Foam2Eigen::PtrList2EigenBC(snapshots);
134 label NBC = snapshots[0].boundaryField().size();
135 Eigen::MatrixXd _corMatrix;
136
137 if (PODnorm == "L2")
138 {
139 _corMatrix = ITHACAutilities::getMassMatrix(snapshots);
140 }
141 else if (PODnorm == "Frobenius")
142 {
143 _corMatrix = ITHACAutilities::getMassMatrix(snapshots, 0, false);
144 }
145
146 if (Pstream::parRun())
147 {
148 reduce(_corMatrix, sumOp<Eigen::MatrixXd>());
149 }
150
151 Eigen::VectorXd eigenValueseig;
152 Eigen::MatrixXd eigenVectoreig;
153 modes.resize(nmodes);
154 Info << "####### Performing the POD using EigenDecomposition " <<
155 fieldName << " #######" << endl;
156 label ncv = snapshots.size();
157 Spectra::DenseSymMatProd<double> op(_corMatrix);
158 Eigen::SelfAdjointEigenSolver<Eigen::MatrixXd> esEg;
159
160 if (para->eigensolver == "spectra")
161 {
162 Spectra::SymEigsSolver<Spectra::DenseSymMatProd<double>> es(op, nmodes, ncv);
163 Info << "Using Spectra EigenSolver " << endl;
164 es.init();
165 es.compute(Spectra::SortRule::LargestAlge);
166 M_Assert(es.info() == Spectra::CompInfo::Successful,
167 "The Eigenvalue Decomposition did not succeed");
168 eigenVectoreig = es.eigenvectors().real();
169 eigenValueseig = es.eigenvalues().real();
170 }
171 else if (para->eigensolver == "eigen")
172 {
173 Info << "Using Eigen EigenSolver " << endl;
174 esEg.compute(_corMatrix);
175 M_Assert(esEg.info() == Eigen::Success,
176 "The Eigenvalue Decomposition did not succeed");
177 eigenVectoreig = esEg.eigenvectors().real().rowwise().reverse().leftCols(
178 nmodes);
179 eigenValueseig = esEg.eigenvalues().real().array().reverse();
180 }
181
182 if (eigenValueseig.array().minCoeff() < 0)
183 {
184 eigenValueseig = eigenValueseig.array() + 2 * abs(
185 eigenValueseig.array().minCoeff());
186 }
187
188 Info << "####### End of the POD for " << snapshots[0].name() << " #######" <<
189 endl;
190 //Eigen::VectorXd eigenValueseigLam =
191 // eigenValueseig.real().array().abs().cwiseInverse().sqrt() ;
192 //Eigen::MatrixXd modesEig = (SnapMatrix * eigenVectoreig) *
193 // eigenValueseigLam.head(nmodes).asDiagonal();
194 Eigen::MatrixXd modesEig = (SnapMatrix * eigenVectoreig);
195 // Computing Normalization factors of the POD Modes
196 Eigen::VectorXd V = ITHACAutilities::getMassMatrixFV(snapshots[0]);
197 Eigen::MatrixXd normFact(nmodes, 1);
198
199 for (label i = 0; i < nmodes; i++)
200 {
201 if (PODnorm == "L2")
202 {
203 normFact(i, 0) = std::sqrt((modesEig.col(i).transpose() * V.asDiagonal() *
204 modesEig.col(i))(0, 0));
205
206 if constexpr(check_vol)
207 {
208 normFact(i, 0) = std::sqrt((modesEig.col(i).transpose() * V.asDiagonal() *
209 modesEig.col(i))(0, 0));
210 }
211 else if constexpr(std::is_same<pointMesh, GeoMesh>::value)
212 {
213 normFact(i, 0) = std::sqrt((modesEig.col(i).transpose() * modesEig.col(i))(0,
214 0));
215 }
216
217 if (Pstream::parRun())
218 {
219 normFact(i, 0) = (modesEig.col(i).transpose() * V.asDiagonal() *
220 modesEig.col(i))(0, 0);
221 }
222 }
223 else if (PODnorm == "Frobenius")
224 {
225 normFact(i, 0) = std::sqrt((modesEig.col(i).transpose() * modesEig.col(i))(0,
226 0));
227
228 if (Pstream::parRun())
229 {
230 normFact(i, 0) = (modesEig.col(i).transpose() * modesEig.col(i))(0, 0);
231 }
232 }
233 }
234
235 if (Pstream::parRun())
236 {
237 reduce(normFact, sumOp<Eigen::MatrixXd>());
238 }
239
240 if (Pstream::parRun())
241 {
242 normFact = normFact.cwiseSqrt();
243 }
244
245 List<Eigen::MatrixXd> modesEigBC;
246 modesEigBC.resize(NBC);
247
248 for (label i = 0; i < NBC; i++)
249 {
250 modesEigBC[i] = (SnapMatrixBC[i] * eigenVectoreig);
251 }
252
253 Info << endl << "####### Normalized Eigenvalues of " << snapshots[0].name()
254 << " #######" << endl << normFact << endl << endl;
255
256 for (label i = 0; i < nmodes; i++)
257 {
258 modesEig.col(i) = modesEig.col(i).array() / normFact(i, 0);
259
260 for (label j = 0; j < NBC; j++)
261 {
262 modesEigBC[j].col(i) = modesEigBC[j].col(i).array() / normFact(i, 0);
263 }
264 }
265
266 for (label i = 0; i < modes.size(); i++)
267 {
268 GeometricField<Type, PatchField, GeoMesh> tmp2(snapshots[0].name(),
269 snapshots[0]);
270 Eigen::VectorXd vec = modesEig.col(i);
271 tmp2 = Foam2Eigen::Eigen2field(tmp2, vec, correctBC);
272
273 for (label k = 0; k < NBC; k++)
274 {
275 ITHACAutilities::assignBC(tmp2, k, modesEigBC[k].col(i));
276 }
277
278 modes.set(i, tmp2.clone());
279 }
280
281 eigenValueseig = eigenValueseig / eigenValueseig.sum();
282 Eigen::VectorXd cumEigenValues(eigenValueseig);
283
284 for (label j = 1; j < cumEigenValues.size(); ++j)
285 {
286 cumEigenValues(j) += cumEigenValues(j - 1);
287 }
288
289 Info << "####### Saving the POD bases for " << snapshots[0].name() <<
290 " #######" << endl;
291
292 if (sup)
293 {
294 ITHACAstream::exportFields(modes, "./ITHACAoutput/supremizer/",
295 snapshots[0].name());
296 }
297 else
298 {
299 ITHACAstream::exportFields(modes, "./ITHACAoutput/POD/", snapshots[0].name());
300 }
301
302 Eigen::saveMarketVector(eigenValueseig,
303 "./ITHACAoutput/POD/Eigenvalues_" + snapshots[0].name(), para->precision,
304 para->outytpe);
305 Eigen::saveMarketVector(cumEigenValues,
306 "./ITHACAoutput/POD/CumEigenvalues_" + snapshots[0].name(), para->precision,
307 para->outytpe);
308 }
309 else
310 {
311 Info << "Reading the existing modes" << endl;
312
313 if (sup == 1)
314 {
315 ITHACAstream::read_fields(modes, fieldName + "sup",
316 "./ITHACAoutput/supremizer/");
317 }
318 else
319 {
320 ITHACAstream::read_fields(modes, fieldName, "./ITHACAoutput/POD/");
321 }
322 }
323}
324
325template void getModes(
326 PtrList<volVectorField>& snapshots, PtrList<volVectorField>& modes,
327 word fieldName, bool podex, bool supex, bool sup, label nmodes,
328 bool correctBC);
329
330template void getModes(
331 PtrList<volScalarField>& snapshots, PtrList<volScalarField>& modes,
332 word fieldName, bool podex, bool supex, bool sup, label nmodes,
333 bool correctBC);
334
335template void getModes(
336 PtrList<surfaceScalarField>& snapshots, PtrList<surfaceScalarField>& modes,
337 word fieldName, bool podex, bool supex, bool sup, label nmodes,
338 bool correctBC);
339template void getModes(
340 PtrList<pointVectorField>& snapshots, PtrList<pointVectorField>& modes,
341 word fieldName, bool podex, bool supex, bool sup, label nmodes,
342 bool correctBC);
343
344template<class Type, template<class> class PatchField, class GeoMesh>
345void getModesMemoryEfficient(
346 GeometricField<Type, PatchField, GeoMesh>& templateField,
347 word snapshotsPath,
348 PtrList<GeometricField<Type, PatchField, GeoMesh >>& modes,
349 word fieldName,
350 bool podex,
351 bool supex,
352 bool sup,
353 label nmodes,
354 bool correctBC,
355 autoPtr<GeometricField<Type, PatchField, GeoMesh >> meanField)
356{
357 // Get parameters instance for POD settings
358 ITHACAparameters* para(ITHACAparameters::getInstance());
359 word PODkey = "POD_" + fieldName;
360 word PODnorm = para->ITHACAdict->lookupOrDefault<word>(PODkey, "L2");
361 // Verify valid norm selection
362 M_Assert(PODnorm == "L2" || PODnorm == "Frobenius",
363 "The PODnorm can be only L2 or Frobenius");
364 Info << "Performing memory efficient POD for " << fieldName
365 << " using the " << PODnorm << " norm" << endl;
366
367 if ((podex == 0 && sup == 0) || (supex == 0 && sup == 1))
368 {
369 // Count number of snapshots in directory (excluding 0/ and constant/)
370 fileName rootPath(".");
371 Foam::Time runTime2(Foam::Time::controlDictName, rootPath, snapshotsPath);
372 label nSnaps = runTime2.times().size() - 2;
373 Info << "Found " << nSnaps << " time directories" << endl;
374
375 // Verify we have at least one snapshot
376 if (nSnaps < 1)
377 {
378 FatalError
379 << "Error: No time directories found in " << snapshotsPath
380 << exit(FatalError);
381 }
382
383 // Set number of modes based on eigensolver type
384 if (para->eigensolver == "spectra")
385 {
386 if (nmodes == 0)
387 {
388 nmodes = nSnaps - 2;
389 }
390
391 M_Assert(nmodes <= nSnaps - 2,
392 "The number of requested modes cannot be bigger than the number of snapshots - 2");
393 }
394 else
395 {
396 if (nmodes == 0)
397 {
398 nmodes = nSnaps;
399 }
400
401 M_Assert(nmodes <= nSnaps,
402 "The number of requested modes cannot be bigger than the number of snapshots");
403 }
404
405 // Initialize correlation matrix and boundary data structures
406 Eigen::MatrixXd _corMatrix(nSnaps, nSnaps);
407 _corMatrix.setZero();
408 List<Eigen::MatrixXd> SnapMatrixBC;
409 label NBC = templateField.boundaryField().size();
410 SnapMatrixBC.resize(NBC);
411
412 // Initialize matrices for boundary conditions
413 for (label i = 0; i < NBC; i++)
414 {
415 SnapMatrixBC[i].resize(templateField.boundaryField()[i].size(), nSnaps);
416 }
417
418 // Build correlation matrix by processing snapshots sequentially
419 for (label i = 0; i < nSnaps; i++)
420 {
421 // Read snapshot i
422 GeometricField<Type, PatchField, GeoMesh> snapI =
423 ITHACAstream::readFieldByIndex(templateField, snapshotsPath, i);
424
425 // Substract mean field if provided
426 if (meanField)
427 {
428 snapI -= *meanField;
429 }
430
431 // Store boundary field data for snapshot i
432 List<Eigen::VectorXd> snapIBC = Foam2Eigen::field2EigenBC(snapI);
433
434 for (label k = 0; k < NBC; k++)
435 {
436 SnapMatrixBC[k].col(i) = snapIBC[k];
437 }
438
439 // Compute correlations with all subsequent snapshots
440 for (label j = i; j < nSnaps; j++)
441 {
442 GeometricField<Type, PatchField, GeoMesh> snapJ =
443 ITHACAstream::readFieldByIndex(templateField, snapshotsPath, j);
444
445 // Subtract mean field if provided
446 if (meanField)
447 {
448 snapJ -= *meanField;
449 }
450
451 // Calculate correlation using specified norm
452 if (PODnorm == "L2")
453 {
454 _corMatrix(i, j) = computeInnerProduct(snapI, snapJ);
455 }
456 else // Frobenius norm
457 {
458 _corMatrix(i, j) = computeFrobeniusInnerProduct(snapI, snapJ);
459 }
460
461 // Matrix is symmetric - copy value to lower triangle
462 if (i != j)
463 {
464 _corMatrix(j, i) = _corMatrix(i, j);
465 }
466 }
467
468 Info << "Processed snapshot " << i + 1 << " of " << nSnaps << endl;
469 }
470
471 // Sum up correlation matrix across processors if running in parallel
472 if (Pstream::parRun())
473 {
474 reduce(_corMatrix, sumOp<Eigen::MatrixXd>());
475 }
476
477 // Solve eigenvalue problem using selected solver
478 Eigen::VectorXd eigenValues;
479 Eigen::MatrixXd eigenVectors;
480 Info << "####### Performing the POD using EigenDecomposition " <<
481 fieldName << " #######" << endl;
482
483 if (para->eigensolver == "spectra")
484 {
485 // Use Spectra solver for large eigenvalue problems
486 Info << "Using Spectra EigenSolver " << endl;
487 Spectra::DenseSymMatProd<double> op(_corMatrix);
488 Spectra::SymEigsSolver<Spectra::DenseSymMatProd<double >>
489 solver(op, nmodes, nSnaps);
490 solver.init();
491 solver.compute(Spectra::SortRule::LargestAlge);
492 M_Assert(solver.info() == Spectra::CompInfo::Successful,
493 "Eigenvalue decomposition failed");
494 eigenVectors = solver.eigenvectors().real();
495 eigenValues = solver.eigenvalues().real();
496 }
497 else if (para->eigensolver == "eigen")
498 {
499 // Use Eigen solver for smaller problems
500 Info << "Using Eigen EigenSolver " << endl;
501 Eigen::SelfAdjointEigenSolver<Eigen::MatrixXd> solver(_corMatrix);
502 M_Assert(solver.info() == Eigen::Success,
503 "Eigenvalue decomposition failed");
504 eigenVectors = solver.eigenvectors().real().rowwise().reverse().leftCols(
505 nmodes);
506 eigenValues = solver.eigenvalues().real().array().reverse();
507 }
508
509 // Handle negative eigenvalues if they occur
510 if (eigenValues.array().minCoeff() < 0)
511 {
512 eigenValues = eigenValues.array() + 2 * abs(
513 eigenValues.array().minCoeff());
514 }
515
516 Info << "####### End of the POD for " << fieldName << " #######" << endl;
517 // Construct POD modes
518 modes.resize(nmodes);
519 // Read first snapshot to get boundary conditions
520 GeometricField<Type, PatchField, GeoMesh> firstSnap =
521 ITHACAstream::readFieldByIndex(templateField, snapshotsPath, 0);
522
523 for (label i = 0; i < nmodes; i++)
524 {
525 // Initialize mode with proper dimensions and boundary conditions
526 GeometricField<Type, PatchField, GeoMesh> modeI
527 (
528 IOobject
529 (
530 templateField.name(),
531 templateField.time().timeName(),
532 templateField.mesh(),
533 IOobject::NO_READ,
534 IOobject::NO_WRITE
535 ),
536 templateField.mesh(),
537 dimensioned<Type>("zero", templateField.dimensions(), Zero),
538 firstSnap.boundaryField().types()
539 );
540
541 // Construct mode as linear combination of snapshots
542 for (label j = 0; j < nSnaps; j++)
543 {
544 GeometricField<Type, PatchField, GeoMesh> snapJ =
545 ITHACAstream::readFieldByIndex(templateField, snapshotsPath, j);
546 modeI += snapJ * eigenVectors(j, i);
547 }
548
549 // Calculate normalization factor based on selected norm
550 scalar normFactor;
551
552 if (PODnorm == "L2")
553 {
554 normFactor = computeInnerProduct(modeI, modeI);
555 }
556 else // Frobenius norm
557 {
558 normFactor = computeFrobeniusInnerProduct(modeI, modeI);
559 }
560
561 if (Pstream::parRun())
562 {
563 reduce(normFactor, sumOp<scalar>());
564 }
565
566 normFactor = Foam::sqrt(normFactor);
567 // Normalize the mode
568 modeI *= dimensionedScalar("normFactor", dimless, 1.0 / normFactor);
569
570 // Apply boundary conditions
571 for (label k = 0; k < NBC; k++)
572 {
573 Eigen::VectorXd bcValues = SnapMatrixBC[k] * eigenVectors.col(i);
574 bcValues = bcValues / normFactor;
575 ITHACAutilities::assignBC(modeI, k, bcValues);
576 }
577
578 if (correctBC)
579 {
580 modeI.correctBoundaryConditions();
581 }
582
583 modes.set(i, modeI.clone());
584 Info << "Constructed mode " << i + 1 << " of " << nmodes << endl;
585 }
586
587 // Save modes to appropriate directory
588 if (sup)
589 {
590 ITHACAstream::exportFields(modes, "./ITHACAoutput/supremizer/", fieldName);
591 }
592 else
593 {
594 ITHACAstream::exportFields(modes, "./ITHACAoutput/POD/", fieldName);
595 }
596
597 // Calculate and save eigenvalue data
598 eigenValues = eigenValues / eigenValues.sum();
599 Eigen::VectorXd cumEigenValues = eigenValues;
600
601 for (label j = 1; j < cumEigenValues.size(); ++j)
602 {
603 cumEigenValues(j) += cumEigenValues(j - 1);
604 }
605
606 // Export eigenvalues
607 Eigen::saveMarketVector(eigenValues,
608 "./ITHACAoutput/POD/Eigenvalues_" + fieldName, para->precision, para->outytpe);
609 Eigen::saveMarketVector(cumEigenValues,
610 "./ITHACAoutput/POD/CumEigenvalues_" + fieldName, para->precision,
611 para->outytpe);
612 }
613 else
614 {
615 // Read existing modes instead of computing new ones
616 Info << "Reading existing modes" << endl;
617
618 if (sup)
619 {
620 ITHACAstream::read_fields(modes, fieldName + "sup",
621 "./ITHACAoutput/supremizer/");
622 }
623 else
624 {
625 ITHACAstream::read_fields(modes, fieldName, "./ITHACAoutput/POD/");
626 }
627 }
628}
629
630template void getModesMemoryEfficient
631(
632 GeometricField<scalar, fvPatchField, volMesh>&,
633 word,
634 PtrList<GeometricField<scalar, fvPatchField, volMesh >>&,
635 word,
636 bool,
637 bool,
638 bool,
639 label,
640 bool,
641 autoPtr<GeometricField<scalar, fvPatchField, volMesh >>
642 );
643
644
645template void getModesMemoryEfficient
646(
647 GeometricField<vector, fvPatchField, volMesh>&,
648 word,
649 PtrList<GeometricField<vector, fvPatchField, volMesh >>&,
650 word,
651 bool,
652 bool,
653 bool,
654 label,
655 bool,
656 autoPtr<GeometricField<vector, fvPatchField, volMesh >>
657 );
658
659template<class Type, template<class> class PatchField, class GeoMesh>
660void getMeanMemoryEfficient(
661 GeometricField<Type, PatchField, GeoMesh>& templateField,
662 word snapshotsPath,
663 autoPtr<GeometricField<Type, PatchField, GeoMesh >>& meanField,
664 bool meanex)
665{
666 // Count number of snapshots in directory (excluding 0/ and constant/)
667 fileName rootPath(".");
668 Foam::Time runTime2(Foam::Time::controlDictName, rootPath, snapshotsPath);
669 label nSnaps = runTime2.times().size() - 2;
670 Info << "Found " << nSnaps << " time directories" << endl;
671
672 // Compute mean field
673 if (!meanex)
674 {
675 Info << "Computing the mean of snapshots" << endl;
676 // Initialize mean field to zero
677 *meanField = templateField * 0.;
678
679 for (int i = 0; i < nSnaps; i++)
680 {
681 // Read snapshot i
682 GeometricField<Type, PatchField, GeoMesh> snapI =
683 ITHACAstream::readFieldByIndex(templateField, snapshotsPath, i);
684 // Sum the snapshots
685 *meanField += snapI;
686 }
687
688 *meanField *= 1. / nSnaps;
689 ITHACAstream::exportSolution(*meanField, "ITHACAoutput", "mean");
690 }
691 else
692 {
693 Info << "Reading the mean of snapshots" << endl;
694 ITHACAstream::readFieldByIndex(templateField, "ITHACAoutput/mean/", 0);
695 }
696}
697
698template void getMeanMemoryEfficient
699(
700 GeometricField<scalar, fvPatchField, volMesh>&,
701 word snapshotsPath,
702 autoPtr<GeometricField<scalar, fvPatchField, volMesh >>&,
703 bool
704 );
705
706template void getMeanMemoryEfficient
707(
708 GeometricField<vector, fvPatchField, volMesh>&,
709 word snapshotsPath,
710 autoPtr<GeometricField<vector, fvPatchField, volMesh >>&,
711 bool
712 );
713
714template<class Type, template<class> class PatchField, class GeoMesh >
715void getWeightedModes(
716 PtrList<GeometricField<Type, PatchField, GeoMesh >>& snapshots,
717 PtrList<GeometricField<Type, PatchField, GeoMesh >>& modes,
718 word fieldName, bool podex, bool supex, bool sup, label nmodes,
719 bool correctBC)
720{
721 ITHACAparameters* para(ITHACAparameters::getInstance());
722
723 if (nmodes == 0)
724 {
725 nmodes = snapshots.size() - 2;
726 }
727
728 M_Assert(nmodes <= snapshots.size() - 2,
729 "The number of requested modes cannot be bigger than the number of Snapshots - 2");
730
731 if ((podex == 0 && sup == 0) || (supex == 0 && sup == 1))
732 {
733 Eigen::MatrixXd SnapMatrix = Foam2Eigen::PtrList2Eigen(snapshots);
734 List<Eigen::MatrixXd> SnapMatrixBC = Foam2Eigen::PtrList2EigenBC(snapshots);
735 label NBC = snapshots[0].boundaryField().size();
736 Eigen::MatrixXd _corMatrix = ITHACAutilities::getMassMatrix(snapshots);
737 Eigen::VectorXd eigenValueseig;
738 Eigen::MatrixXd eigenVectoreig;
739 modes.resize(nmodes);
740 Info << "####### Performing the POD using EigenDecomposition " <<
741 snapshots[0].name() << " #######" << endl;
742 label ncv = snapshots.size();
743 Spectra::DenseSymMatProd<double> op(_corMatrix);
744 Eigen::SelfAdjointEigenSolver<Eigen::MatrixXd> esEg;
745 Spectra::SymEigsSolver<Spectra::DenseSymMatProd<double >>
746 es(op, nmodes, ncv);
747
748 if (para->eigensolver == "spectra")
749 {
750 Info << "Using Spectra EigenSolver " << endl;
751 es.init();
752 es.compute(Spectra::SortRule::LargestAlge);
753 M_Assert(es.info() == Spectra::CompInfo::Successful,
754 "The Eigenvalue Decomposition did not succeed");
755 eigenVectoreig = es.eigenvectors().real();
756 eigenValueseig = es.eigenvalues().real();
757 }
758 else if (para->eigensolver == "eigen")
759 {
760 Info << "Using Eigen EigenSolver " << endl;
761 esEg.compute(_corMatrix);
762 M_Assert(esEg.info() == Eigen::Success,
763 "The Eigenvalue Decomposition did not succeed");
764 eigenVectoreig = esEg.eigenvectors().real().rowwise().reverse().leftCols(
765 nmodes);
766 eigenValueseig = esEg.eigenvalues().real().reverse();
767 }
768
769 Info << "####### End of the POD for " << snapshots[0].name() << " #######" <<
770 endl;
771 Eigen::VectorXd eigenValueseigLam =
772 eigenValueseig.real().array().cwiseInverse().sqrt() ;
773 Eigen::VectorXd eigenValueseigWeigted = eigenValueseig.head(
774 nmodes).real().array() ;
775 Eigen::MatrixXd modesEig = (SnapMatrix * eigenVectoreig) *
776 eigenValueseigLam.head(nmodes).asDiagonal() *
777 eigenValueseigWeigted.asDiagonal();
778 List<Eigen::MatrixXd> modesEigBC;
779 modesEigBC.resize(NBC);
780
781 for (label i = 0; i < NBC; i++)
782 {
783 modesEigBC[i] = (SnapMatrixBC[i] * eigenVectoreig) *
784 eigenValueseigLam.asDiagonal() * eigenValueseigWeigted.asDiagonal();
785 }
786
787 for (label i = 0; i < modes.size(); i++)
788 {
789 GeometricField<Type, PatchField, GeoMesh> tmp2(snapshots[0].name(),
790 snapshots[0]);
791 Eigen::VectorXd vec = modesEig.col(i);
792 tmp2 = Foam2Eigen::Eigen2field(tmp2, vec, correctBC);
793
794 for (label k = 0; k < tmp2.boundaryField().size(); k++)
795 {
796 ITHACAutilities::assignBC(tmp2, k, modesEigBC[k].col(i));
797 }
798
799 modes.set(i, tmp2.clone());
800 }
801
802 eigenValueseig = eigenValueseig / eigenValueseig.sum();
803 Eigen::VectorXd cumEigenValues(eigenValueseig);
804
805 for (label j = 1; j < cumEigenValues.size(); ++j)
806 {
807 cumEigenValues(j) += cumEigenValues(j - 1);
808 }
809
810 Info << "####### Saving the POD bases for " << snapshots[0].name() <<
811 " #######" << endl;
812
813 //exportBases(modes, snapshots, sup);
814 if (sup)
815 {
816 ITHACAstream::exportFields(modes, "./ITHACAoutput/supremizer/",
817 snapshots[0].name());
818 }
819 else
820 {
821 ITHACAstream::exportFields(modes, "./ITHACAoutput/POD/", snapshots[0].name());
822 }
823
824 Eigen::saveMarketVector(eigenValueseig,
825 "./ITHACAoutput/POD/Eigenvalues_" + snapshots[0].name(), para->precision,
826 para->outytpe);
827 Eigen::saveMarketVector(cumEigenValues,
828 "./ITHACAoutput/POD/CumEigenvalues_" + snapshots[0].name(), para->precision,
829 para->outytpe);
830 }
831 else
832 {
833 Info << "Reading the existing modes" << endl;
834
835 if (sup == 1)
836 {
837 ITHACAstream::read_fields(modes, fieldName + "sup",
838 "./ITHACAoutput/supremizer/");
839 }
840 else
841 {
842 ITHACAstream::read_fields(modes, fieldName, "./ITHACAoutput/POD/");
843 }
844 }
845}
846
847template void getWeightedModes(
848 PtrList<volScalarField>& snapshots, PtrList<volScalarField>& modes,
849 word fieldName, bool podex, bool supex, bool sup, label nmodes,
850 bool correctBC);
851
852template void getWeightedModes(
853 PtrList<volVectorField>& snapshots, PtrList<volVectorField>& modes,
854 word fieldName, bool podex, bool supex, bool sup, label nmodes,
855 bool correctBC);
856
857template<class Type, template<class> class PatchField, class GeoMesh>
858void getModesSVD(
859 PtrList<GeometricField<Type, PatchField, GeoMesh >>& snapshots,
860 PtrList<GeometricField<Type, PatchField, GeoMesh >>& modes,
861 word fieldName, bool podex, bool supex, bool sup, label nmodes,
862 bool correctBC)
863{
864 ITHACAparameters* para(ITHACAparameters::getInstance());
865
866 if ((podex == 0 && sup == 0) || (supex == 0 && sup == 1))
867 {
868 PtrList<volVectorField> Bases;
869 modes.resize(nmodes);
870 Info << "####### Performing POD using Singular Value Decomposition for " <<
871 snapshots[0].name() << " #######" << endl;
872 Eigen::MatrixXd SnapMatrix = Foam2Eigen::PtrList2Eigen(snapshots);
873 Eigen::VectorXd V = ITHACAutilities::getMassMatrixFV(snapshots[0]);
874 Eigen::VectorXd V3dSqrt = V.array().sqrt();
875 Eigen::VectorXd V3dInv = V3dSqrt.array().cwiseInverse();
876 auto VMsqr = V3dSqrt.asDiagonal();
877 auto VMsqrInv = V3dInv.asDiagonal();
878 Eigen::MatrixXd SnapMatrix2 = VMsqr * SnapMatrix;
879 Eigen::JacobiSVD<Eigen::MatrixXd> svd(SnapMatrix2,
880 Eigen::ComputeThinU | Eigen::ComputeThinV);
881 Info << "####### End of the POD for " << snapshots[0].name() << " #######" <<
882 endl;
883 Eigen::VectorXd eigenValueseig;
884 Eigen::MatrixXd eigenVectoreig;
885 eigenValueseig = svd.singularValues().real();
886 eigenVectoreig = svd.matrixU().real();
887 Eigen::MatrixXd modesEig = VMsqrInv * eigenVectoreig;
888 GeometricField<Type, PatchField, GeoMesh> tmb_bu(snapshots[0].name(),
889 snapshots[0] * 0);
890
891 for (label i = 0; i < nmodes; i++)
892 {
893 Eigen::VectorXd vec = modesEig.col(i);
894 tmb_bu = Foam2Eigen::Eigen2field(tmb_bu, vec, correctBC);
895 modes.set(i, tmb_bu.clone());
896 }
897
898 eigenValueseig = eigenValueseig / eigenValueseig.sum();
899 Eigen::VectorXd cumEigenValues(eigenValueseig);
900
901 for (label j = 1; j < cumEigenValues.size(); ++j)
902 {
903 cumEigenValues(j) += cumEigenValues(j - 1);
904 }
905
906 Info << "####### Saving the POD bases for " << snapshots[0].name() <<
907 " #######" << endl;
908
909 if (sup)
910 {
911 ITHACAstream::exportFields(modes, "./ITHACAoutput/supremizer/",
912 snapshots[0].name());
913 }
914 else
915 {
916 ITHACAstream::exportFields(modes, "./ITHACAoutput/POD/", snapshots[0].name());
917 }
918
919 //exportBases(modes, snapshots, sup);
920 Eigen::saveMarketVector(eigenValueseig,
921 "./ITHACAoutput/POD/Eigenvalues_" + snapshots[0].name(), para->precision,
922 para->outytpe);
923 Eigen::saveMarketVector(cumEigenValues,
924 "./ITHACAoutput/POD/CumEigenvalues_" + snapshots[0].name(), para->precision,
925 para->outytpe);
926 }
927 else
928 {
929 Info << "Reading the existing modes" << endl;
930
931 if (sup == 1)
932 {
933 ITHACAstream::read_fields(modes, fieldName + "sup",
934 "./ITHACAoutput/supremizer/");
935 }
936 else
937 {
938 ITHACAstream::read_fields(modes, fieldName, "./ITHACAoutput/POD/");
939 }
940 }
941}
942
943template void getModesSVD(
944 PtrList<volScalarField>& snapshots, PtrList<volScalarField>& modes,
945 word fieldName, bool podex, bool supex, bool sup, label nmodes,
946 bool correctBC);
947
948template void getModesSVD(
949 PtrList<volVectorField>& snapshots, PtrList<volVectorField>& modes,
950 word fieldName, bool podex, bool supex, bool sup, label nmodes,
951 bool correctBC);
952
954template<>
955Eigen::MatrixXd corMatrix(PtrList<volScalarField>& snapshots)
956{
957 Info << "########## Filling the correlation matrix for " << snapshots[0].name()
958 << "##########" << endl;
959 Eigen::MatrixXd matrix( snapshots.size(), snapshots.size());
960
961 for (label i = 0; i < snapshots.size(); i++)
962 {
963 for (label j = 0; j <= i; j++)
964 {
965 matrix(i, j) = fvc::domainIntegrate(snapshots[i] * snapshots[j]).value();
966 }
967 }
968
969 for (label i = 1; i < snapshots.size(); i++)
970 {
971 for (label j = 0; j < i; j++)
972 {
973 matrix(j, i) = matrix(i, j);
974 }
975 }
976
977 return matrix;
978}
979
980
982template<>
983Eigen::MatrixXd corMatrix(PtrList<volVectorField>& snapshots)
984{
985 Info << "########## Filling the correlation matrix for " << snapshots[0].name()
986 << "##########" << endl;
987 Eigen::MatrixXd matrix( snapshots.size(), snapshots.size());
988
989 for (label i = 0; i < snapshots.size(); i++)
990 {
991 for (label j = 0; j <= i; j++)
992 {
993 matrix(i, j) = fvc::domainIntegrate(snapshots[i] & snapshots[j]).value();
994 }
995 }
996
997 for (label i = 1; i < snapshots.size(); i++)
998 {
999 for (label j = 0; j < i; j++)
1000 {
1001 matrix(j, i) = matrix(i, j);
1002 }
1003 }
1004
1005 return matrix;
1006}
1007
1009template<>
1010Eigen::MatrixXd corMatrix(List<Eigen::SparseMatrix<double >>&
1011 snapshots)
1012{
1013 Info << "########## Filling the correlation matrix for the matrix list ##########"
1014 << endl;
1015 Eigen::MatrixXd matrix(snapshots.size(), snapshots.size());
1016
1017 for (label i = 0; i < snapshots.size(); i++)
1018 {
1019 for (label j = 0; j <= i; j++)
1020 {
1021 double res = 0;
1022
1023 for (label k = 0; k < snapshots[i].cols(); k++)
1024 {
1025 res += snapshots[i].col(k).dot(snapshots[j].col(k));
1026 }
1027
1028 matrix(i, j) = res;
1029 }
1030 }
1031
1032 for (label i = 1; i < snapshots.size(); i++)
1033 {
1034 for (label j = 0; j < i; j++)
1035 {
1036 matrix(j, i) = matrix(i, j);
1037 }
1038 }
1039
1040 return matrix;
1041}
1042
1044template<>
1045Eigen::MatrixXd corMatrix(List<Eigen::VectorXd>& snapshots)
1046{
1047 Info << "########## Filling the correlation matrix for the matrix list ##########"
1048 << endl;
1049 Eigen::MatrixXd matrix( snapshots.size(), snapshots.size());
1050
1051 for (label i = 0; i < snapshots.size(); i++)
1052 {
1053 for (label j = 0; j <= i; j++)
1054 {
1055 matrix(i, j) = (snapshots[i].transpose() * snapshots[j]).trace();
1056 }
1057 }
1058
1059 for (label i = 1; i < snapshots.size(); i++)
1060 {
1061 for (label j = 0; j < i; j++)
1062 {
1063 matrix(j, i) = matrix(i, j);
1064 }
1065 }
1066
1067 return matrix;
1068}
1069
1070
1071
1073template<class Type, template<class> class PatchField, class GeoMesh>
1074void exportBases(PtrList<GeometricField<Type, PatchField, GeoMesh >>& s,
1075 PtrList<GeometricField<Type, PatchField, GeoMesh >>& bases,
1076 word fieldName, bool sup)
1077{
1078 if (sup)
1079 {
1080 fileName fieldname;
1081
1082 for (label i = 0; i < s.size(); i++)
1083 {
1084 mkDir("./ITHACAoutput/supremizer/" + name(i + 1));
1085 fieldname = "./ITHACAoutput/supremizer/" + name(i + 1) + "/" +
1086 bases[i].name();
1087 OFstream os(fieldname);
1088 bases[i].writeHeader(os);
1089 os << s[i] << endl;
1090 }
1091 }
1092 else
1093 {
1094 fileName fieldname;
1095
1096 for (label i = 0; i < s.size(); i++)
1097 {
1098 mkDir("./ITHACAoutput/POD/" + name(i + 1));
1099 fieldname = "./ITHACAoutput/POD/" + name(i + 1) + "/" + bases[i].name();
1100 OFstream os(fieldname);
1101 bases[i].writeHeader(os);
1102 os << s[i] << endl;
1103 }
1104 }
1105}
1106
1107template void exportBases(PtrList<volVectorField>& s,
1108 PtrList<volVectorField>& bases, word fieldName, bool sup);
1109template void exportBases(PtrList<volScalarField>& s,
1110 PtrList<volScalarField>& bases, word fieldName, bool sup);
1111
1112void exportEigenvalues(scalarField Eigenvalues, fileName name,
1113 bool sup)
1114{
1115 if (sup)
1116 {
1117 fileName fieldname;
1118 mkDir("./ITHACAoutput/supremizer/");
1119 fieldname = "./ITHACAoutput/supremizer/Eigenvalues_" + name;
1120 OFstream os(fieldname);
1121
1122 for (label i = 0; i < Eigenvalues.size(); i++)
1123 {
1124 os << Eigenvalues[i] << endl;
1125 }
1126 }
1127 else
1128 {
1129 fileName fieldname;
1130 mkDir("./ITHACAoutput/POD/");
1131 fieldname = "./ITHACAoutput/POD/Eigenvalues_" + name;
1132 OFstream os(fieldname);
1133
1134 for (label i = 0; i < Eigenvalues.size(); i++)
1135 {
1136 os << Eigenvalues[i] << endl;
1137 }
1138 }
1139}
1140
1141void exportcumEigenvalues(scalarField cumEigenvalues, fileName name,
1142 bool sup)
1143{
1144 if (sup)
1145 {
1146 fileName fieldname;
1147 mkDir("./ITHACAoutput/supremizer/");
1148 fieldname = "./ITHACAoutput/supremizer/cumEigenvalues_" + name;
1149 OFstream os(fieldname);
1150
1151 for (label i = 0; i < cumEigenvalues.size(); i++)
1152 {
1153 os << cumEigenvalues[i] << endl;
1154 }
1155 }
1156 else
1157 {
1158 fileName fieldname;
1159 mkDir("./ITHACAoutput/POD/");
1160 fieldname = "./ITHACAoutput/POD/cumEigenvalues_" + name;
1161 OFstream os(fieldname);
1162
1163 for (label i = 0; i < cumEigenvalues.size(); i++)
1164 {
1165 os << cumEigenvalues[i] << endl;
1166 }
1167 }
1168}
1169
1170
1171std::tuple<List<Eigen::SparseMatrix<double >>, List<Eigen::VectorXd >>
1172DEIMmodes(List<Eigen::SparseMatrix<double >>& A,
1173 List<Eigen::VectorXd>& b, label nmodesA, label nmodesB, word MatrixName)
1174{
1175 ITHACAparameters* para(ITHACAparameters::getInstance());
1176 List<Eigen::SparseMatrix<double >> ModesA(nmodesA);
1177 List<Eigen::VectorXd> ModesB(nmodesB);
1178
1179 if (!ITHACAutilities::check_folder("./ITHACAoutput/DEIM/" + MatrixName))
1180 {
1181 if (nmodesA > A.size() - 2 || nmodesB > A.size() - 2 )
1182 {
1183 Info <<
1184 "The number of requested modes cannot be bigger than the number of Snapshots - 2"
1185 << endl;
1186 exit(0);
1187 }
1188
1189 scalarField eigenValuesA(nmodesA);
1190 scalarField cumEigenValuesA(nmodesA);
1191 scalarField eigenValuesB(nmodesB);
1192 scalarField cumEigenValuesB(nmodesB);
1193 List<scalarField> eigenVectorA(nmodesA);
1194 List<scalarField> eigenVectorB(nmodesB);
1195
1196 for (label i = 0; i < nmodesA; i++)
1197 {
1198 eigenVectorA[i].setSize(A.size());
1199 }
1200
1201 for (label i = 0; i < nmodesB; i++)
1202 {
1203 eigenVectorB[i].setSize(A.size());
1204 }
1205
1206 Eigen::MatrixXd corMatrixA = corMatrix(A);
1207 Eigen::MatrixXd corMatrixB = corMatrix(b);
1208 Info << "####### Performing the POD for the Matrix List #######" << endl;
1209 Spectra::DenseSymMatProd<double> opA(corMatrixA);
1210 Spectra::DenseSymMatProd<double> opB(corMatrixB);
1211 Spectra::SymEigsSolver<Spectra:: DenseSymMatProd<double >>
1212 esA(opA, nmodesA, nmodesA + 10);
1213 Spectra::SymEigsSolver<Spectra:: DenseSymMatProd<double
1214 >>
1215 esB(opB, nmodesB, nmodesB + 10);
1216 esA.init();
1217 esB.init();
1218 esA.compute(Spectra::SortRule::LargestAlge);
1219
1220 if (nmodesB != 1)
1221 {
1222 esB.compute();
1223 }
1224
1225 Info << "####### End of the POD for the Matrix List #######" << endl;
1226 Eigen::VectorXd eigenValueseigA;
1227 Eigen::MatrixXd eigenVectorseigA;
1228 Eigen::VectorXd eigenValueseigB;
1229 Eigen::MatrixXd eigenVectorseigB;
1230
1231 if (esA.info() == Spectra::CompInfo::Successful)
1232 {
1233 eigenValueseigA = esA.eigenvalues().real();
1234 eigenVectorseigA = esA.eigenvectors().real();
1235
1236 if (esB.info() == Spectra::CompInfo::Successful && nmodesB != 1)
1237 {
1238 eigenValueseigB = esB.eigenvalues().real();
1239 eigenVectorseigB = esB.eigenvectors().real();
1240 }
1241 else if (nmodesB == 1)
1242 {
1243 eigenValueseigB.resize(1);
1244 eigenVectorseigB.resize(A.size(), nmodesB);
1245 eigenValueseigB(0) = 1;
1246 eigenVectorseigB = eigenVectorseigB * 0;
1247 eigenVectorseigB(0, 0) = 1;
1248 }
1249 else
1250 {
1251 Info << "The Eigenvalue solver in ITHACAPOD.H did not converge, exiting the code"
1252 << endl;
1253 exit(0);
1254 }
1255 }
1256 else
1257 {
1258 Info << "The Eigenvalue solver in ITHACAPOD.H did not converge, exiting the code"
1259 << endl;
1260 exit(0);
1261 }
1262
1263 for (label i = 0; i < nmodesA; i++)
1264 {
1265 eigenValuesA[i] = eigenValueseigA(i) / eigenValueseigA.sum();
1266 }
1267
1268 for (label i = 0; i < nmodesB; i++)
1269 {
1270 eigenValuesB[i] = eigenValueseigB(i) / eigenValueseigB.sum();
1271 }
1272
1273 for (label i = 0; i < nmodesA; i++)
1274 {
1275 for (label k = 0; k < A.size(); k++)
1276 {
1277 eigenVectorA[i][k] = eigenVectorseigA(k, i);
1278 }
1279 }
1280
1281 for (label i = 0; i < nmodesB; i++)
1282 {
1283 for (label k = 0; k < A.size(); k++)
1284 {
1285 eigenVectorB[i][k] = eigenVectorseigB(k, i);
1286 }
1287 }
1288
1289 cumEigenValuesA[0] = eigenValuesA[0];
1290 cumEigenValuesB[0] = eigenValuesB[0];
1291
1292 for (label i = 1; i < nmodesA; i++)
1293 {
1294 cumEigenValuesA[i] = cumEigenValuesA[i - 1] + eigenValuesA[i];
1295 }
1296
1297 for (label i = 1; i < nmodesB; i++)
1298 {
1299 cumEigenValuesB[i] = cumEigenValuesB[i - 1] + eigenValuesB[i];
1300 }
1301
1302 Eigen::SparseMatrix<double> tmp_A;
1303 Eigen::VectorXd tmp_B;
1304
1305 for (label i = 0; i < nmodesA; i++)
1306 {
1307 tmp_A = eigenVectorA[i][0] * A[0];
1308
1309 for (label k = 1; k < A.size(); k++)
1310 {
1311 tmp_A += eigenVectorA[i][k] * A[k];
1312 }
1313
1314 ModesA[i] = tmp_A;
1315 }
1316
1317 for (label i = 0; i < nmodesB; i++)
1318 {
1319 tmp_B = eigenVectorB[i][0] * b[0];
1320
1321 for (label k = 1; k < A.size(); k++)
1322 {
1323 tmp_B += eigenVectorB[i][k] * b[k];
1324 }
1325
1326 ModesB[i] = tmp_B;
1327 }
1328
1329 ITHACAstream::exportList(eigenValuesA,
1330 "./ITHACAoutput/DEIM/" + MatrixName + "/", "eigenValuesA_" + MatrixName);
1331 ITHACAstream::exportList(cumEigenValuesA,
1332 "./ITHACAoutput/DEIM/" + MatrixName + "/", "cumEigenValuesA_" + MatrixName);
1333 ITHACAstream::exportList(eigenValuesB,
1334 "./ITHACAoutput/DEIM/" + MatrixName + "/", "eigenValuesB_" + MatrixName);
1335 ITHACAstream::exportList(cumEigenValuesB,
1336 "./ITHACAoutput/DEIM/" + MatrixName + "/", "cumEigenValuesB_" + MatrixName);
1337
1338 for (label i = 0; i < ModesA.size(); i++)
1339 {
1340 ITHACAstream::SaveSparseMatrix(ModesA[i],
1341 "./ITHACAoutput/DEIM/" + MatrixName + "/", "A_" + MatrixName + name(i));
1342 }
1343
1344 for (label i = 0; i < ModesB.size(); i++)
1345 {
1346 ITHACAstream::SaveDenseMatrix(ModesB[i],
1347 "./ITHACAoutput/DEIM/" + MatrixName + "/", "B_" + MatrixName + name(i));
1348 }
1349 }
1350 else
1351 {
1352 for (label i = 0; i < nmodesA; i++)
1353 {
1354 ITHACAstream::ReadSparseMatrix(ModesA[i],
1355 "./ITHACAoutput/DEIM/" + MatrixName + "/", "A_" + MatrixName + name(i));
1356 }
1357
1358 for (label i = 0; i < nmodesB; i++)
1359 {
1360 ITHACAstream::ReadDenseMatrix(ModesB[i],
1361 "./ITHACAoutput/DEIM/" + MatrixName + "/", "B_" + MatrixName + name(i));
1362 }
1363 }
1364
1365 std::tuple <List<Eigen::SparseMatrix<double >>, List<Eigen::VectorXd >> tupla;
1366 tupla = std::make_tuple(ModesA, ModesB);
1367 return tupla;
1368}
1369
1370void GrammSchmidt(Eigen::MatrixXd& Matrix)
1371{
1372 Eigen::MatrixXd Ortho = Matrix;
1373 Ortho = Matrix;
1374
1375 for (label i = 0; i < Matrix.cols(); i++)
1376 {
1377 for (label k = 0; k < i; k++)
1378 {
1379 double num = Ortho.col(k).transpose() * Matrix.col(i);
1380 double den = (Ortho.col(k).transpose() * Ortho.col(k));
1381 double fact = num / den;
1382 Ortho.col(i) -= fact * Ortho.col(k) ;
1383 }
1384
1385 Ortho.col(i).normalize();
1386 }
1387
1388 Matrix = Ortho;
1389}
1390
1391template<class Type, template<class> class PatchField, class GeoMesh>
1392void getModes(
1393 PtrList<GeometricField<Type, PatchField, GeoMesh >>& snapshots,
1394 PtrList<GeometricField<Type, PatchField, GeoMesh >>& modes,
1395 PtrList<volScalarField>& Volumes, word fieldName, bool podex, bool supex,
1396 bool sup, label nmodes, bool correctBC)
1397{
1398 ITHACAparameters* para(ITHACAparameters::getInstance());
1399
1400 if (nmodes == 0 && para->eigensolver == "spectra")
1401 {
1402 nmodes = snapshots.size() - 2;
1403 }
1404
1405 if (nmodes == 0 && para->eigensolver == "eigen")
1406 {
1407 nmodes = snapshots.size();
1408 }
1409
1410 if (para->eigensolver == "spectra")
1411 {
1412 M_Assert(nmodes <= snapshots.size() - 2,
1413 "The number of requested modes cannot be bigger than the number of Snapshots - 2");
1414 }
1415
1416 if ((podex == 0 && sup == 0) || (supex == 0 && sup == 1))
1417 {
1418 Eigen::MatrixXd SnapMatrix = Foam2Eigen::PtrList2Eigen(snapshots);
1419 List<Eigen::MatrixXd> SnapMatrixBC = Foam2Eigen::PtrList2EigenBC(snapshots);
1420 label NBC = snapshots[0].boundaryField().size();
1421 Eigen::MatrixXd V = Foam2Eigen::PtrList2Eigen(Volumes);
1422 Eigen::MatrixXd _corMatrix(snapshots.size(), snapshots.size());
1423 Info << "Filling the correlation matrix for field " << snapshots[0].name() <<
1424 endl;
1425
1426 for (label i = 0; i < snapshots.size(); i++)
1427 {
1428 for (label j = 0; j <= i; j++)
1429 {
1430 Eigen::VectorXd Mij = (V.col(i).array() * V.col(j).array());
1431 Mij = Mij.array().abs().sqrt();
1432 _corMatrix(i, j) = SnapMatrix.col(i).transpose() * Mij.asDiagonal() *
1433 SnapMatrix.col(j);
1434 }
1435 }
1436
1437 Info << endl;
1438
1439 for (label i = 1; i < snapshots.size(); i++)
1440 {
1441 for (label j = 0; j < i; j++)
1442 {
1443 _corMatrix(j, i) = _corMatrix(i, j);
1444 }
1445 }
1446
1447 Eigen::VectorXd eigenValueseig;
1448 Eigen::MatrixXd eigenVectoreig;
1449 modes.resize(nmodes);
1450 Info << "####### Performing the POD using EigenDecomposition for " <<
1451 snapshots[0].name() << " #######" << endl;
1452 label ncv = snapshots.size();
1453 Spectra::DenseSymMatProd<double> op(_corMatrix);
1454 Eigen::SelfAdjointEigenSolver<Eigen::MatrixXd> esEg;
1455
1456 if (para->eigensolver == "spectra")
1457 {
1458 Spectra::SymEigsSolver< Spectra::DenseSymMatProd<double >>
1459 es(op, nmodes, ncv);
1460 Info << "Using Spectra EigenSolver " << endl;
1461 es.init();
1462 es.compute(Spectra::SortRule::LargestAlge);
1463 M_Assert(es.info() == Spectra::CompInfo::Successful,
1464 "The Eigenvalue Decomposition did not succeed");
1465 eigenVectoreig = es.eigenvectors().real();
1466 eigenValueseig = es.eigenvalues().real();
1467 }
1468 else if (para->eigensolver == "eigen")
1469 {
1470 Info << "Using Eigen EigenSolver " << endl;
1471 esEg.compute(_corMatrix);
1472 M_Assert(esEg.info() == Eigen::Success,
1473 "The Eigenvalue Decomposition did not succeed");
1474 eigenVectoreig = esEg.eigenvectors().real().rowwise().reverse().leftCols(
1475 nmodes);
1476 eigenValueseig = esEg.eigenvalues().real().reverse().head(nmodes);
1477 }
1478
1479 Info << "####### End of the POD for " << snapshots[0].name() << " #######" <<
1480 endl;
1481 Eigen::VectorXd eigenValueseigLam =
1482 eigenValueseig.real().array().cwiseInverse().sqrt() ;
1483 Eigen::MatrixXd modesEig = (SnapMatrix * eigenVectoreig) *
1484 eigenValueseigLam.asDiagonal();
1485 List<Eigen::MatrixXd> modesEigBC;
1486 modesEigBC.resize(NBC);
1487
1488 for (label i = 0; i < NBC; i++)
1489 {
1490 modesEigBC[i] = (SnapMatrixBC[i] * eigenVectoreig) *
1491 eigenValueseigLam.asDiagonal();
1492 }
1493
1494 for (label i = 0; i < modes.size(); i++)
1495 {
1496 GeometricField<Type, PatchField, GeoMesh> tmp2(snapshots[0].name(),
1497 snapshots[0] * 0);
1498 Eigen::VectorXd vec = modesEig.col(i);
1499 tmp2 = Foam2Eigen::Eigen2field(tmp2, vec, correctBC);
1500
1501 // Adjusting boundary conditions
1502 for (label k = 0; k < tmp2.boundaryField().size(); k++)
1503 {
1504 ITHACAutilities::assignBC(tmp2, k, modesEigBC[k].col(i));
1505 }
1506
1507 modes.set(i, tmp2.clone());
1508 }
1509
1510 eigenValueseig = eigenValueseig / eigenValueseig.sum();
1511 Eigen::VectorXd cumEigenValues(eigenValueseig);
1512
1513 for (label j = 1; j < cumEigenValues.size(); ++j)
1514 {
1515 cumEigenValues(j) += cumEigenValues(j - 1);
1516 }
1517
1518 Info << "####### Saving the POD bases for " << snapshots[0].name() << " #######"
1519 << endl;
1520 exportBases(modes, snapshots, fieldName, sup);
1521 Eigen::saveMarketVector(eigenValueseig,
1522 "./ITHACAoutput/POD/Eigenvalues_" + snapshots[0].name(), para->precision,
1523 para->outytpe);
1524 Eigen::saveMarketVector(cumEigenValues,
1525 "./ITHACAoutput/POD/CumEigenvalues_" + snapshots[0].name(), para->precision,
1526 para->outytpe);
1527 }
1528 else
1529 {
1530 Info << "Reading the existing modes" << endl;
1531
1532 if (sup == 1)
1533 {
1534 ITHACAstream::read_fields(modes, "Usup", "./ITHACAoutput/supremizer/");
1535 }
1536 else
1537 {
1538 ITHACAstream::read_fields(modes, fieldName, "./ITHACAoutput/POD/");
1539 }
1540 }
1541}
1542
1543template void getModes(
1544 PtrList<volScalarField>& snapshots, PtrList<volScalarField>& modes,
1545 PtrList<volScalarField>& Volumes, word fieldName, bool podex, bool supex,
1546 bool sup, label nmodes, bool correctBC);
1547
1548template void getModes(
1549 PtrList<volVectorField>& snapshots, PtrList<volVectorField>& modes,
1550 PtrList<volScalarField>& Volumes, word fieldName, bool podex, bool supex,
1551 bool sup, label nmodes, bool correctBC);
1552
1553template<typename type_matrix>
1554std::tuple<List<Eigen::SparseMatrix<double >>, List<Eigen::VectorXd >>
1555DEIMmodes(PtrList<type_matrix>& MatrixList, label nmodesA, label nmodesB,
1556 word MatrixName)
1557{
1558 ITHACAparameters* para(ITHACAparameters::getInstance());
1559 List<Eigen::SparseMatrix<double >> ModesA(nmodesA);
1560 List<Eigen::VectorXd> ModesB(nmodesB);
1561
1562 if (!ITHACAutilities::check_folder("./ITHACAoutput/DEIM/" + MatrixName))
1563 {
1564 M_Assert(nmodesA <= MatrixList.size() - 2
1565 && nmodesB <= MatrixList.size() - 2,
1566 "The number of requested modes cannot be bigger than the number of Snapshots - 2");
1567 std::tuple<List<Eigen::SparseMatrix<double >>, List<Eigen::VectorXd >> snapshots
1568 =
1569 Foam2Eigen::LFvMatrix2LSM(MatrixList);
1570 Eigen::MatrixXd corMatrixA = corMatrix(std::get<0>(snapshots));
1571 Eigen::MatrixXd corMatrixB = corMatrix(std::get<1>(snapshots));
1572 Eigen::VectorXd eigenValueseigA;
1573 Eigen::MatrixXd eigenVectorseigA;
1574 Eigen::VectorXd eigenValueseigB;
1575 Eigen::MatrixXd eigenVectorseigB;
1576
1577 if (para->eigensolver == "spectra")
1578 {
1579 Info << "####### Performing the POD decomposition for the Matrix List using Spectra #######"
1580 << endl;
1581 Spectra::DenseSymMatProd<double> opA(corMatrixA);
1582 Spectra::DenseSymMatProd<double> opB(corMatrixB);
1583 label ncvA = MatrixList.size();
1584 label ncvB = MatrixList.size();
1585 Spectra::SymEigsSolver<Spectra:: DenseSymMatProd<double >>
1586 esA(opA, nmodesA, ncvA);
1587 Spectra::SymEigsSolver<Spectra:: DenseSymMatProd<double
1588 >>
1589 esB(opB, nmodesB, ncvB);
1590 esA.init();
1591 esB.init();
1592 esA.compute(Spectra::SortRule::LargestAlge);
1593 M_Assert(esA.info() == Spectra::CompInfo::Successful,
1594 "The Eigenvalue Decomposition did not succeed");
1595
1596 if (nmodesB != 1)
1597 {
1598 esB.compute(Spectra::SortRule::LargestAlge);
1599 M_Assert(esB.info() == Spectra::CompInfo::Successful,
1600 "The Eigenvalue Decomposition did not succeed");
1601 }
1602
1603 Info << "####### End of the POD decomposition for the Matrix List #######" <<
1604 endl;
1605 eigenValueseigA = esA.eigenvalues().real();
1606 eigenVectorseigA = esA.eigenvectors().real();
1607
1608 if (nmodesB != 1)
1609 {
1610 eigenValueseigB = esB.eigenvalues().real();
1611 eigenVectorseigB = esB.eigenvectors().real();
1612 }
1613 else
1614 {
1615 eigenValueseigB.resize(1);
1616 eigenVectorseigB.resize(MatrixList.size(), nmodesB);
1617 eigenValueseigB(0) = 1;
1618 eigenVectorseigB = eigenVectorseigB * 0;
1619 eigenVectorseigB(0, 0) = 1;
1620 }
1621 }
1622 else if (para->eigensolver == "eigen")
1623 {
1624 Info << "####### Performing the POD decomposition for the Matrix List using Eigen #######"
1625 << endl;
1626 Eigen::SelfAdjointEigenSolver<Eigen::MatrixXd> esEgA;
1627 esEgA.compute(corMatrixA);
1628 Eigen::SelfAdjointEigenSolver<Eigen::MatrixXd> esEgB;
1629 esEgB.compute(corMatrixB);
1630 M_Assert(esEgA.info() == Eigen::Success,
1631 "The Eigenvalue Decomposition did not succeed");
1632 M_Assert(esEgB.info() == Eigen::Success,
1633 "The Eigenvalue Decomposition did not succeed");
1634 eigenVectorseigA = esEgA.eigenvectors().real().rowwise().reverse().leftCols(
1635 nmodesA);
1636 eigenValueseigA = esEgA.eigenvalues().real().reverse().head(nmodesA);
1637 eigenVectorseigB = esEgB.eigenvectors().real().rowwise().reverse().leftCols(
1638 nmodesB);
1639 eigenValueseigB = esEgB.eigenvalues().real().reverse().head(nmodesB);
1640 }
1641
1642 Eigen::SparseMatrix<double> tmp_A;
1643 Eigen::VectorXd tmp_B;
1644
1645 for (label i = 0; i < nmodesA; i++)
1646 {
1647 tmp_A = eigenVectorseigA(0, i) * std::get<0>(snapshots)[0];
1648
1649 for (label k = 1; k < MatrixList.size(); k++)
1650 {
1651 tmp_A += eigenVectorseigA(k, i) * std::get<0>(snapshots)[k];
1652 }
1653
1654 ModesA[i] = tmp_A;
1655 }
1656
1657 for (label i = 0; i < nmodesB; i++)
1658 {
1659 tmp_B = eigenVectorseigB(0, i) * std::get<1>(snapshots)[0];
1660
1661 for (label k = 1; k < MatrixList.size(); k++)
1662 {
1663 tmp_B += eigenVectorseigB(k, i) * std::get<1>(snapshots)[k];
1664 }
1665
1666 ModesB[i] = tmp_B;
1667 }
1668
1669 eigenValueseigA = eigenValueseigA / eigenValueseigA.sum();
1670 eigenValueseigB = eigenValueseigB / eigenValueseigB.sum();
1671 Eigen::VectorXd cumEigenValuesA(eigenValueseigA);
1672 Eigen::VectorXd cumEigenValuesB(eigenValueseigB);
1673
1674 for (label j = 1; j < cumEigenValuesA.size(); ++j)
1675 {
1676 cumEigenValuesA(j) += cumEigenValuesA(j - 1);
1677 }
1678
1679 for (label j = 1; j < cumEigenValuesB.size(); ++j)
1680 {
1681 cumEigenValuesB(j) += cumEigenValuesB(j - 1);
1682 }
1683
1684 for (label i = 0; i < ModesA.size(); i++)
1685 {
1686 ITHACAstream::SaveSparseMatrix(ModesA[i],
1687 "./ITHACAoutput/DEIM/" + MatrixName + "/", "A_" + MatrixName + name(i));
1688 }
1689
1690 for (label i = 0; i < ModesB.size(); i++)
1691 {
1692 ITHACAstream::SaveDenseMatrix(ModesB[i],
1693 "./ITHACAoutput/DEIM/" + MatrixName + "/", "B_" + MatrixName + name(i));
1694 }
1695
1696 Eigen::saveMarketVector(eigenValueseigA,
1697 "./ITHACAoutput/DEIM/" + MatrixName + "/eigenValuesA", para->precision,
1698 para->outytpe);
1699 Eigen::saveMarketVector(eigenValueseigB,
1700 "./ITHACAoutput/DEIM/" + MatrixName + "/eigenValuesB", para->precision,
1701 para->outytpe);
1702 Eigen::saveMarketVector(cumEigenValuesA,
1703 "./ITHACAoutput/DEIM/" + MatrixName + "/cumEigenValuesA", para->precision,
1704 para->outytpe);
1705 Eigen::saveMarketVector(cumEigenValuesB,
1706 "./ITHACAoutput/DEIM/" + MatrixName + "/cumEigenValuesB", para->precision,
1707 para->outytpe);
1708 }
1709 else
1710 {
1711 for (label i = 0; i < nmodesA; i++)
1712 {
1713 ITHACAstream::ReadSparseMatrix(ModesA[i],
1714 "./ITHACAoutput/DEIM/" + MatrixName + "/", "A_" + MatrixName + name(i));
1715 }
1716
1717 for (label i = 0; i < nmodesB; i++)
1718 {
1719 ITHACAstream::ReadDenseMatrix(ModesB[i],
1720 "./ITHACAoutput/DEIM/" + MatrixName + "/", "B_" + MatrixName + name(i));
1721 }
1722 }
1723
1724 std::tuple <List<Eigen::SparseMatrix<double >>, List<Eigen::VectorXd >> tupla;
1725 tupla = std::make_tuple(ModesA, ModesB);
1726 return tupla;
1727}
1728
1729template std::tuple<List<Eigen::SparseMatrix<double >>, List<Eigen::VectorXd >>
1730DEIMmodes(PtrList<fvScalarMatrix>& MatrixList, label nmodesA,
1731 label nmodesB,
1732 word MatrixName);
1733
1734template std::tuple<List<Eigen::SparseMatrix<double >>, List<Eigen::VectorXd >>
1735DEIMmodes(PtrList<fvVectorMatrix>& MatrixList, label nmodesA,
1736 label nmodesB,
1737 word MatrixName);
1738
1739template<class Type, template<class> class PatchField, class GeoMesh >
1740PtrList<GeometricField<Type, PatchField, GeoMesh >> DEIMmodes(
1741 PtrList<GeometricField<Type, PatchField, GeoMesh >>& snapshots, label nmodes,
1742 word FunctionName, word fieldName)
1743{
1744 ITHACAparameters* para(ITHACAparameters::getInstance());
1745 word PODkey = "POD_" + fieldName;
1746 word PODnorm = para->ITHACAdict->lookupOrDefault<word>(PODkey, "L2");
1747 M_Assert(PODnorm == "L2" ||
1748 PODnorm == "Frobenius", "The PODnorm can be only L2 or Frobenius");
1749 Info << "Performing POD for " << fieldName << " using the " << PODnorm <<
1750 " norm" << endl;
1751 PtrList<GeometricField<Type, fvPatchField, volMesh >> modes;
1752 bool correctBC = true;
1753
1754 if (nmodes == 0 && para->eigensolver == "spectra")
1755 {
1756 nmodes = snapshots.size() - 2;
1757 }
1758
1759 if (nmodes == 0 && para->eigensolver == "eigen")
1760 {
1761 nmodes = snapshots.size();
1762 }
1763
1764 if (para->eigensolver == "spectra")
1765 {
1766 M_Assert(nmodes <= snapshots.size() - 2,
1767 "The number of requested modes cannot be bigger than the number of Snapshots - 2");
1768 }
1769
1770 if (!ITHACAutilities::check_folder("./ITHACAoutput/DEIM/" + FunctionName))
1771 {
1772 Eigen::MatrixXd SnapMatrix = Foam2Eigen::PtrList2Eigen(snapshots);
1773 List<Eigen::MatrixXd> SnapMatrixBC = Foam2Eigen::PtrList2EigenBC(snapshots);
1774 label NBC = snapshots[0].boundaryField().size();
1775 Eigen::MatrixXd _corMatrix;
1776
1777 if (PODnorm == "L2")
1778 {
1779 _corMatrix = ITHACAutilities::getMassMatrix(snapshots);
1780 }
1781 else if (PODnorm == "Frobenius")
1782 {
1783 _corMatrix = ITHACAutilities::getMassMatrix(snapshots, 0, false);
1784 }
1785
1786 if (Pstream::parRun())
1787 {
1788 reduce(_corMatrix, sumOp<Eigen::MatrixXd>());
1789 }
1790
1791 Eigen::VectorXd eigenValueseig;
1792 Eigen::MatrixXd eigenVectoreig;
1793 modes.resize(nmodes);
1794 Info << "####### Performing the POD using EigenDecomposition " <<
1795 fieldName << " #######" << endl;
1796 label ncv = snapshots.size();
1797 Spectra::DenseSymMatProd<double> op(_corMatrix);
1798 Eigen::SelfAdjointEigenSolver<Eigen::MatrixXd> esEg;
1799
1800 if (para->eigensolver == "spectra")
1801 {
1802 Spectra::SymEigsSolver< Spectra::DenseSymMatProd<double >>
1803 es(op, nmodes, ncv);
1804 Info << "Using Spectra EigenSolver " << endl;
1805 es.init();
1806 es.compute(Spectra::SortRule::LargestAlge);
1807 M_Assert(es.info() == Spectra::CompInfo::Successful,
1808 "The Eigenvalue Decomposition did not succeed");
1809 eigenVectoreig = es.eigenvectors().real();
1810 eigenValueseig = es.eigenvalues().real();
1811 }
1812 else if (para->eigensolver == "eigen")
1813 {
1814 Info << "Using Eigen EigenSolver " << endl;
1815 esEg.compute(_corMatrix);
1816 M_Assert(esEg.info() == Eigen::Success,
1817 "The Eigenvalue Decomposition did not succeed");
1818 eigenVectoreig = esEg.eigenvectors().real().rowwise().reverse().leftCols(
1819 nmodes);
1820 eigenValueseig = esEg.eigenvalues().real().array().reverse();
1821 }
1822
1823 if (eigenValueseig.array().minCoeff() < 0)
1824 {
1825 eigenValueseig = eigenValueseig.array() + 2 * abs(
1826 eigenValueseig.array().minCoeff());
1827 }
1828
1829 Info << "####### End of the POD for " << snapshots[0].name() << " #######" <<
1830 endl;
1831 Eigen::MatrixXd modesEig = (SnapMatrix * eigenVectoreig);
1832 Eigen::VectorXd V = ITHACAutilities::getMassMatrixFV(snapshots[0]);
1833 Eigen::MatrixXd normFact(nmodes, 1);
1834
1835 for (label i = 0; i < nmodes; i++)
1836 {
1837 if (PODnorm == "L2")
1838 {
1839 normFact(i, 0) = std::sqrt((modesEig.col(i).transpose() * V.asDiagonal() *
1840 modesEig.col(i))(0, 0));
1841
1842 if (Pstream::parRun())
1843 {
1844 normFact(i, 0) = (modesEig.col(i).transpose() * V.asDiagonal() *
1845 modesEig.col(i))(0, 0);
1846 }
1847 }
1848 else if (PODnorm == "Frobenius")
1849 {
1850 normFact(i, 0) = std::sqrt((modesEig.col(i).transpose() * modesEig.col(i))(0,
1851 0));
1852
1853 if (Pstream::parRun())
1854 {
1855 normFact(i, 0) = (modesEig.col(i).transpose() * modesEig.col(i))(0, 0);
1856 }
1857 }
1858 }
1859
1860 if (Pstream::parRun())
1861 {
1862 reduce(normFact, sumOp<Eigen::MatrixXd>());
1863 }
1864
1865 if (Pstream::parRun())
1866 {
1867 normFact = normFact.cwiseSqrt();
1868 }
1869
1870 List<Eigen::MatrixXd> modesEigBC;
1871 modesEigBC.resize(NBC);
1872
1873 for (label i = 0; i < NBC; i++)
1874 {
1875 modesEigBC[i] = (SnapMatrixBC[i] * eigenVectoreig);
1876 }
1877
1878 Info << normFact << endl;
1879
1880 for (label i = 0; i < nmodes; i++)
1881 {
1882 modesEig.col(i) = modesEig.col(i).array() / normFact(i, 0);
1883
1884 for (label j = 0; j < NBC; j++)
1885 {
1886 modesEigBC[j].col(i) = modesEigBC[j].col(i).array() / normFact(i, 0);
1887 }
1888 }
1889
1890 for (label i = 0; i < modes.size(); i++)
1891 {
1892 GeometricField<Type, PatchField, GeoMesh> tmp2(snapshots[0].name(),
1893 snapshots[0]);
1894 Eigen::VectorXd vec = modesEig.col(i);
1895 tmp2 = Foam2Eigen::Eigen2field(tmp2, vec, correctBC);
1896
1897 for (label k = 0; k < NBC; k++)
1898 {
1899 ITHACAutilities::assignBC(tmp2, k, modesEigBC[k].col(i));
1900 }
1901
1902 modes.set(i, tmp2.clone());
1903 }
1904
1905 eigenValueseig = eigenValueseig / eigenValueseig.sum();
1906 Eigen::VectorXd cumEigenValues(eigenValueseig);
1907
1908 for (label j = 1; j < cumEigenValues.size(); ++j)
1909 {
1910 cumEigenValues(j) += cumEigenValues(j - 1);
1911 }
1912
1913 Info << "####### Saving the POD bases for " << snapshots[0].name() <<
1914 " #######" << endl;
1915 ITHACAutilities::createSymLink("./ITHACAoutput/DEIM");
1916
1917 for (label i = 0; i < modes.size(); i++)
1918 {
1919 ITHACAstream::exportSolution(modes[i], name(i + 1), "./ITHACAoutput/DEIM",
1920 fieldName);
1921 }
1922
1923 Eigen::saveMarketVector(eigenValueseig,
1924 "./ITHACAoutput/DEIM/eigenValues_" + fieldName, para->precision,
1925 para->outytpe);
1926 Eigen::saveMarketVector(cumEigenValues,
1927 "./ITHACAoutput/DEIM/cumEigenValues_" + fieldName, para->precision,
1928 para->outytpe);
1929 }
1930 else
1931 {
1932 Info << "Reading the existing modes" << endl;
1933 ITHACAstream::read_fields(modes, fieldName, "./ITHACAoutput/DEIM/");
1934 }
1935
1936 return modes;
1937}
1938
1939template<class Field_type, class Field_type_2>
1940void getModes(
1941 PtrList<Field_type>& snapshots, PtrList<Field_type>& modes,
1942 PtrList<Field_type_2>& fields2, word fieldName, bool podex, bool supex,
1943 bool sup, label nmodes, bool correctBC)
1944{
1945 ITHACAparameters* para(ITHACAparameters::getInstance());
1946
1947 if ((podex == 0 && sup == 0) || (supex == 0 && sup == 1))
1948 {
1949 if (para->eigensolver == "spectra" )
1950 {
1951 if (nmodes == 0)
1952 {
1953 nmodes = snapshots.size() - 2;
1954 }
1955
1956 M_Assert(nmodes <= snapshots.size() - 2,
1957 "The number of requested modes cannot be bigger than the number of Snapshots - 2");
1958 }
1959 else
1960 {
1961 if (nmodes == 0)
1962 {
1963 nmodes = snapshots.size();
1964 }
1965
1966 M_Assert(nmodes <= snapshots.size(),
1967 "The number of requested modes cannot be bigger than the number of Snapshots");
1968 }
1969
1970 Eigen::MatrixXd SnapMatrix = Foam2Eigen::PtrList2Eigen(fields2);
1971 List<Eigen::MatrixXd> SnapMatrixBC = Foam2Eigen::PtrList2EigenBC(fields2);
1972 label NBC = fields2[0].boundaryField().size();
1973 auto VM = ITHACAutilities::getMassMatrixFV(fields2[0]);
1974 Eigen::MatrixXd _corMatrix = SnapMatrix.transpose() * VM.asDiagonal() *
1975 SnapMatrix;
1976
1977 if (Pstream::parRun())
1978 {
1979 reduce(_corMatrix, sumOp<Eigen::MatrixXd>());
1980 }
1981
1982 Eigen::VectorXd eigenValueseig;
1983 Eigen::MatrixXd eigenVectoreig;
1984 modes.resize(nmodes);
1985 Info << "####### Performing the POD using EigenDecomposition " <<
1986 fields2[0].name() << " #######" << endl;
1987 label ncv = snapshots.size();
1988 Spectra::DenseSymMatProd<double> op(_corMatrix);
1989 Eigen::SelfAdjointEigenSolver<Eigen::MatrixXd> esEg;
1990
1991 if (para->eigensolver == "spectra")
1992 {
1993 Spectra::SymEigsSolver<Spectra::DenseSymMatProd<double >>
1994 es(op, nmodes, ncv);
1995 Info << "Using Spectra EigenSolver " << endl;
1996 es.init();
1997 es.compute(Spectra::SortRule::LargestAlge);
1998 M_Assert(es.info() == Spectra::CompInfo::Successful,
1999 "The Eigenvalue Decomposition did not succeed");
2000 eigenVectoreig = es.eigenvectors().real();
2001 eigenValueseig = es.eigenvalues().real();
2002 }
2003 else if (para->eigensolver == "eigen")
2004 {
2005 Info << "Using Eigen EigenSolver " << endl;
2006 esEg.compute(_corMatrix);
2007 M_Assert(esEg.info() == Eigen::Success,
2008 "The Eigenvalue Decomposition did not succeed");
2009 eigenVectoreig = esEg.eigenvectors().real().rowwise().reverse().leftCols(
2010 nmodes);
2011 eigenValueseig = esEg.eigenvalues().real().reverse().head(nmodes);
2012 }
2013
2014 Info << "####### End of the POD for " << snapshots[0].name() << " #######" <<
2015 endl;
2016 Eigen::VectorXd eigenValueseigLam =
2017 eigenValueseig.real().array().cwiseInverse().abs().sqrt() ;
2018 Eigen::MatrixXd modesEig = (SnapMatrix * eigenVectoreig) *
2019 eigenValueseigLam.asDiagonal();
2020 List<Eigen::MatrixXd> modesEigBC;
2021 modesEigBC.resize(NBC);
2022
2023 for (label i = 0; i < modes.size(); i++)
2024 {
2025 Field_type tmp2(snapshots[0].name(), snapshots[0] * 0);
2026
2027 for (label j = 0; j < snapshots.size(); j++)
2028 {
2029 tmp2 += eigenValueseigLam(i) * eigenVectoreig.col(i)(j) * snapshots[j];
2030 }
2031
2032 modes.set(i, tmp2.clone());
2033 }
2034
2035 eigenValueseig = eigenValueseig / eigenValueseig.sum();
2036 Eigen::VectorXd cumEigenValues(eigenValueseig);
2037
2038 for (label j = 1; j < cumEigenValues.size(); ++j)
2039 {
2040 cumEigenValues(j) += cumEigenValues(j - 1);
2041 }
2042
2043 Info << "####### Saving the POD bases for " << snapshots[0].name() <<
2044 " #######" << endl;
2045
2046 if (sup)
2047 {
2048 ITHACAstream::exportFields(modes, "./ITHACAoutput/supremizer/",
2049 snapshots[0].name());
2050 }
2051 else
2052 {
2053 ITHACAstream::exportFields(modes, "./ITHACAoutput/POD/", snapshots[0].name());
2054 }
2055
2056 Eigen::saveMarketVector(eigenValueseig,
2057 "./ITHACAoutput/POD/Eigenvalues_" + snapshots[0].name(), para->precision,
2058 para->outytpe);
2059 Eigen::saveMarketVector(cumEigenValues,
2060 "./ITHACAoutput/POD/CumEigenvalues_" + snapshots[0].name(), para->precision,
2061 para->outytpe);
2062 }
2063 else
2064 {
2065 Info << "Reading the existing modes" << endl;
2066
2067 if (sup == 1)
2068 {
2069 ITHACAstream::read_fields(modes, fieldName, "./ITHACAoutput/supremizer/");
2070 }
2071 else
2072 {
2073 ITHACAstream::read_fields(modes, fieldName, "./ITHACAoutput/POD/");
2074 }
2075 }
2076}
2077
2078template void getModes(
2079 PtrList<surfaceScalarField>& snapshots, PtrList<surfaceScalarField>& modes,
2080 PtrList<volVectorField>& fields2, word fieldName, bool podex, bool supex,
2081 bool sup, label nmodes, bool correctBC);
2082
2083template void getModes(
2084 PtrList<volScalarField>& snapshots, PtrList<volScalarField>& modes,
2085 PtrList<volVectorField>& fields2, word fieldName, bool podex, bool supex,
2086 bool sup, label nmodes, bool correctBC);
2087
2088template PtrList<volScalarField>
2089DEIMmodes(
2090 PtrList<volScalarField>& SnapShotsMatrix,
2091 label nmodes,
2092 word FunctionName, word FieldName);
2093
2094template PtrList<volVectorField>
2095DEIMmodes(
2096 PtrList<volVectorField>& SnapShotsMatrix,
2097 label nmodes,
2098 word FunctionName, word FieldName);
2099
2100template<>
2101scalar computeInnerProduct(
2102 const GeometricField<scalar, fvPatchField, volMesh>& field1,
2103 const GeometricField<scalar, fvPatchField, volMesh>& field2)
2104{
2105 return fvc::domainIntegrate(field1 * field2).value();
2106}
2107
2108template<>
2109scalar computeInnerProduct(
2110 const GeometricField<vector, fvPatchField, volMesh>& field1,
2111 const GeometricField<vector, fvPatchField, volMesh>& field2)
2112{
2113 return fvc::domainIntegrate(field1 & field2).value();
2114}
2115
2116template<>
2117scalar computeFrobeniusInnerProduct(
2118 const GeometricField<scalar, fvPatchField, volMesh>& field1,
2119 const GeometricField<scalar, fvPatchField, volMesh>& field2)
2120{
2121 return sum(field1.primitiveField() * field2.primitiveField());
2122}
2123
2124template<>
2125scalar computeFrobeniusInnerProduct(
2126 const GeometricField<vector, fvPatchField, volMesh>& field1,
2127 const GeometricField<vector, fvPatchField, volMesh>& field2)
2128{
2129 return sum(field1.primitiveField()&field2.primitiveField());
2130}
2131
2132}
EigenFunctions.H
Header file of the EigenFunctions class.
ITHACAPOD.H
Header file of the ITHACAPOD class.
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:533
Foam2Eigen::field2EigenBC
static List< Eigen::VectorXd > field2EigenBC(GeometricField< scalar, PatchField, GeoMesh > &field)
Convert an OpenFOAM scalar field to a List of Eigen Vectors, one for each boundary.
Definition Foam2Eigen.C:253
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:1354
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:273
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:665
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::outytpe
std::_Ios_Fmtflags outytpe
type of output format can be fixed or scientific
Definition ITHACAparameters.H:62
ITHACAparameters::getInstance
static ITHACAparameters * getInstance(const fvMesh &mesh, Time &localTime)
Gets an instance of ITHACAparameters, to be used if the instance is not existing.
Definition ITHACAparameters.C:45
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:1141
ITHACAPOD::getModesMemoryEfficient
void getModesMemoryEfficient(GeometricField< Type, PatchField, GeoMesh > &templateField, word snapshotsPath, PtrList< GeometricField< Type, PatchField, GeoMesh > > &modes, word fieldName, bool podex, bool supex, bool sup, label nmodes, bool correctBC, autoPtr< GeometricField< Type, PatchField, GeoMesh > > meanField)
Gets the modes in a memory-efficient manner.
Definition ITHACAPOD.C:345
ITHACAPOD::exportEigenvalues
void exportEigenvalues(scalarField Eigenvalues, fileName name, bool sup)
Exports the eigenvalues as a txt file.
Definition ITHACAPOD.C:1112
ITHACAPOD::getMeanMemoryEfficient
void getMeanMemoryEfficient(GeometricField< Type, PatchField, GeoMesh > &templateField, word snapshotsPath, autoPtr< GeometricField< Type, PatchField, GeoMesh > > &meanField, bool meanex)
Gets the mean field in a memory-efficient manner.
Definition ITHACAPOD.C:660
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:1172
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:1074
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:858
ITHACAPOD::corMatrix
Eigen::MatrixXd corMatrix(PtrList< volScalarField > &snapshots)
Construct the Correlation Matrix for Scalar Field.
Definition ITHACAPOD.C:955
ITHACAPOD::GrammSchmidt
void GrammSchmidt(Eigen::MatrixXd &Matrix)
Performs GrammSchmidt orthonormalization on an Eigen Matrix.
Definition ITHACAPOD.C:1370
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:715
ITHACAstream::SaveSparseMatrix
void SaveSparseMatrix(Eigen::SparseMatrix< T, Nrows, IND > &m, word folder, word MatrixName)
Export an Eigen sparse matrix into bynary format file.
ITHACAstream::readFieldByIndex
GeometricField< Type, PatchField, GeoMesh > readFieldByIndex(const GeometricField< Type, PatchField, GeoMesh > &field, fileName casename, label index)
Function to read a single field by index from a folder.
Definition ITHACAstream.C:454
ITHACAstream::ReadDenseMatrix
void ReadDenseMatrix(MatrixType &Matrix, word folder, word MatrixName)
Read a dense matrix from a binary format file.
Definition ITHACAstream.H:740
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:958
ITHACAstream::exportSolution
void exportSolution(GeometricField< Type, PatchField, GeoMesh > &s, fileName subfolder, fileName folder, word fieldName)
Export a field to file in a certain folder and subfolder.
Definition ITHACAstream.C:992
ITHACAstream::exportList
void exportList(T &list, word folder, word filename)
Export a list to file.
Definition ITHACAstream.C:1240
ITHACAstream::SaveDenseMatrix
void SaveDenseMatrix(MatrixType &Matrix, word folder, word MatrixName)
Save a dense matrix to a binary format file.
Definition ITHACAstream.H:723
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:517
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:313
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:80
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
Generated by   doxygen 1.16.1