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