Loading...
Searching...
No Matches
ITHACAcoeffsMass.C
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-------------------------------------------------------------------------------
12License
13 This file is part of ITHACA-FV
14 ITHACA-FV is free software: you can redistribute it and/or modify
15 it under the terms of the GNU Lesser General Public License as published by
16 the Free Software Foundation, either version 3 of the License, or
17 (at your option) any later version.
18 ITHACA-FV is distributed in the hope that it will be useful,
19 but WITHOUT ANY WARRANTY; without even the implied warranty of
20 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
21 GNU Lesser General Public License for more details.
22 You should have received a copy of the GNU Lesser General Public License
23 along with ITHACA-FV. If not, see <http://www.gnu.org/licenses/>.
24\*---------------------------------------------------------------------------*/
25#include "ITHACAcoeffsMass.H"
26
27namespace ITHACAutilities
28{
29
30template<class Type, template<class> class PatchField, class GeoMesh>
31PtrList<GeometricField<Type, PatchField, GeoMesh >> reconstructFromCoeff(
32 PtrList<GeometricField<Type, PatchField, GeoMesh >>& modes,
33 Eigen::MatrixXd& coeff_matrix, label Nmodes)
34{
35 PtrList<GeometricField<Type, PatchField, GeoMesh >> rec_field;
36 rec_field.resize(0);
37
38 for (label k = 0; k < coeff_matrix.cols(); k++)
39 {
40 for (label i = 0; i < Nmodes; i++)
41 if (i == 0)
42 {
43 rec_field.append(modes[i] * coeff_matrix(i, k));
44 }
45 else
46 {
47 rec_field[k] += modes[i] * coeff_matrix(i, k);
48 }
49 }
50
51 return rec_field;
52}
53
54template PtrList<GeometricField<scalar, fvPatchField, volMesh >>
56 PtrList<GeometricField<scalar, fvPatchField, volMesh >>& modes,
57 Eigen::MatrixXd& coeff_matrix, label Nmodes);
58template PtrList<GeometricField<vector, fvPatchField, volMesh >>
60 PtrList<GeometricField<vector, fvPatchField, volMesh >>& modes,
61 Eigen::MatrixXd& coeff_matrix, label Nmodes);
62template PtrList<GeometricField<tensor, fvPatchField, volMesh >>
64 PtrList<GeometricField<tensor, fvPatchField, volMesh >>& modes,
65 Eigen::MatrixXd& coeff_matrix, label Nmodes);
66
67template<class Type, template<class> class PatchField, class GeoMesh >
68Eigen::MatrixXd getMassMatrix(
69 PtrList<GeometricField<Type, PatchField, GeoMesh >>& modes, label Nmodes,
70 bool consider_volumes)
71{
72 label Msize;
73
74 if (Nmodes == 0)
75 {
76 Msize = modes.size();
77 }
78 else
79 {
80 Msize = Nmodes;
81 }
82
83 M_Assert(modes.size() >= Msize,
84 "The Number of requested modes is larger then the available quantity.");
85 Eigen::MatrixXd F = Foam2Eigen::PtrList2Eigen(modes);
86 Eigen::MatrixXd M = Eigen::MatrixXd::Zero(Msize, Msize);
87 constexpr bool check_vol = std::is_same<volMesh, GeoMesh>::value
88 || std::is_same<surfaceMesh, GeoMesh>::value;
89
90 if constexpr(check_vol)
91 {
92 Eigen::VectorXd V = getMassMatrixFV(modes[0]);
93
94 // If V is too big then we need a for loop to avoid RAM overload
95 if (V.size() > 1e6)
96 {
97 for (int i = 0; i < V.size(); i++)
98 {
99 M += V(i) * F.transpose().topRows(Msize).col(i) * F.leftCols(Msize).row(i);
100 }
101 }
102 else
103 {
104 M = F.transpose().topRows(Msize) * V.asDiagonal() * F.leftCols(Msize);
105 }
106 }
107 else if constexpr(std::is_same<pointMesh, GeoMesh>::value)
108 {
109 M = F.transpose().topRows(Msize) * F.leftCols(Msize);
110 }
111
112 if (Pstream::parRun())
113 {
114 reduce(M, sumOp<Eigen::MatrixXd>());
115 }
116
117 return M;
118}
119
120template Eigen::MatrixXd getMassMatrix(
121 PtrList<GeometricField<scalar, fvPatchField, volMesh >>& modes, label Nmodes,
122 bool consider_volumes);
123
124template Eigen::MatrixXd getMassMatrix(
125 PtrList<GeometricField<scalar, fvsPatchField, surfaceMesh >>& modes,
126 label Nmodes,
127 bool consider_volumes = false);
128
129template Eigen::MatrixXd getMassMatrix(
130 PtrList<GeometricField<vector, fvPatchField, volMesh >>& modes, label Nmodes,
131 bool consider_volumes);
132
133template Eigen::MatrixXd getMassMatrix(
134 PtrList<GeometricField<tensor, fvPatchField, volMesh >>& modes, label Nmodes,
135 bool consider_volumes);
136
137template<class Type, template<class> class PatchField, class GeoMesh>
138Eigen::MatrixXd getMassMatrix(
139 PtrList<GeometricField<Type, PatchField, GeoMesh >>& modes,
140 PtrList<GeometricField<Type, PatchField, GeoMesh >>& modes2, label Nmodes,
141 bool consider_volumes)
142{
143 label Msize, Msize2;
144
145 if (Nmodes == 0)
146 {
147 Msize = modes.size();
148 Msize2 = modes2.size();
149 }
150 else
151 {
152 Msize = Nmodes;
153 Msize2 = Nmodes;
154 }
155
156 M_Assert(modes.size() >= Msize,
157 "The Number of requested modes is larger then the available quantity.");
158 M_Assert(modes2.size() >= Msize2,
159 "The Number of requested modes is larger then the available quantity.");
160 Eigen::MatrixXd F = Foam2Eigen::PtrList2Eigen(modes);
161 Eigen::MatrixXd F2 = Foam2Eigen::PtrList2Eigen(modes2);
162 Eigen::MatrixXd M = Eigen::MatrixXd::Zero(Msize, Msize2);
163
164 if (consider_volumes)
165 {
166 Eigen::VectorXd V = getMassMatrixFV(modes[0]);
167
168 // If V is too big then we need a for loop to avoid RAM overload
169 if (V.size() > 1e6)
170 {
171 for (int i = 0; i < V.size(); i++)
172 {
173 M += V(i) * F.transpose().topRows(Msize).col(i) * F2.leftCols(Msize2).row(i);
174 }
175 }
176 else
177 {
178 // Classical M computation
179 M = F.transpose().topRows(Msize) * V.asDiagonal() * F2.leftCols(Msize2);
180 }
181 }
182 else
183 {
184 M = F.transpose().topRows(Msize) * F2.leftCols(Msize2);
185 }
186
187 if (Pstream::parRun())
188 {
189 reduce(M, sumOp<Eigen::MatrixXd>());
190 }
191
192 return M;
193}
194
195template Eigen::MatrixXd getMassMatrix(
196 PtrList<GeometricField<scalar, fvPatchField, volMesh >>& modes,
197 PtrList<GeometricField<scalar, fvPatchField, volMesh >>& modes2, label Nmodes,
198 bool consider_volumes);
199
200template Eigen::MatrixXd getMassMatrix(
201 PtrList<GeometricField<scalar, fvsPatchField, surfaceMesh >>& modes,
202 PtrList<GeometricField<scalar, fvsPatchField, surfaceMesh >>& modes2,
203 label Nmodes,
204 bool consider_volumes = false);
205
206template Eigen::MatrixXd getMassMatrix(
207 PtrList<GeometricField<vector, fvPatchField, volMesh >>& modes,
208 PtrList<GeometricField<vector, fvPatchField, volMesh >>& modes2, label Nmodes,
209 bool consider_volumes);
210
211template Eigen::MatrixXd getMassMatrix(
212 PtrList<GeometricField<tensor, fvPatchField, volMesh >>& modes,
213 PtrList<GeometricField<tensor, fvPatchField, volMesh >>& modes2, label Nmodes,
214 bool consider_volumes);
215
216template Eigen::MatrixXd getMassMatrix(
217 PtrList<GeometricField<vector, pointPatchField, pointMesh>>& modes,
218 label Nmodes,
219 bool consider_volumes);
220
221template<class Type, template<class> class PatchField, class GeoMesh>
222Eigen::MatrixXd getMassMatrix(
223 PtrList<GeometricField<Type, PatchField, GeoMesh >>
224 & modes,
225 PtrList<GeometricField<Type, PatchField, GeoMesh >>& modes2,
226 Eigen::VectorXd weights,
227 label Nmodes,
228 bool consider_volumes)
229{
230 label Msize, Msize2;
231
232 if (Nmodes == 0)
233 {
234 Msize = modes.size();
235 Msize2 = modes2.size();
236 }
237 else
238 {
239 Msize = Nmodes;
240 Msize2 = Nmodes;
241 }
242
243 M_Assert(modes.size() >= Msize,
244 "The Number of requested modes is larger then the available quantity.");
245 M_Assert(modes2.size() >= Msize2,
246 "The Number of requested modes is larger then the available quantity.");
247 Eigen::MatrixXd F = Foam2Eigen::PtrList2Eigen(modes);
248 Eigen::MatrixXd F2 = Foam2Eigen::PtrList2Eigen(modes2);
249 Eigen::MatrixXd M = Eigen::MatrixXd::Zero(Msize, Msize2);
250
251 if (consider_volumes)
252 {
253 Eigen::VectorXd V = getMassMatrixFV(modes[0]);
254
255 // If V is too big then we need a for loop to avoid RAM overload
256 if (V.size() > 1e6)
257 {
258 for (int i = 0; i < V.size(); i++)
259 {
260 M += (V(i) * weights(i)) * F.transpose().topRows(Msize).col(i) * F2.leftCols(
261 Msize2).row(i);
262 }
263 }
264 else
265 {
266 // Classical M computation
267 M = F.transpose().topRows(Msize) * V.asDiagonal() * ( weights.asDiagonal()
268 * F2.leftCols(Msize2) );
269 }
270 }
271 else
272 {
273 M = F.transpose().topRows(Msize) * weights.asDiagonal() * F2.leftCols(Msize2);
274 }
275
276 if (Pstream::parRun())
277 {
278 reduce(M, sumOp<Eigen::MatrixXd>());
279 }
280
281 return M;
282}
283
284template Eigen::MatrixXd getMassMatrix(
285 PtrList<GeometricField<scalar, fvPatchField, volMesh >>& modes,
286 PtrList<GeometricField<scalar, fvPatchField, volMesh >>& modes2,
287 Eigen::VectorXd weights,
288 label Nmodes = 0,
289 bool consider_volumes);
290
291template Eigen::MatrixXd getMassMatrix(
292 PtrList<GeometricField<scalar, fvsPatchField, surfaceMesh >>& modes,
293 PtrList<GeometricField<scalar, fvsPatchField, surfaceMesh >>& modes2,
294 Eigen::VectorXd weights,
295 label Nmodes = 0,
296 bool consider_volumes);
297
298template Eigen::MatrixXd getMassMatrix(
299 PtrList<GeometricField<vector, fvPatchField, volMesh >>& modes,
300 PtrList<GeometricField<vector, fvPatchField, volMesh >>& modes2,
301 Eigen::VectorXd weights,
302 label Nmodes = 0,
303 bool consider_volumes);
304
305template Eigen::MatrixXd getMassMatrix(
306 PtrList<GeometricField<tensor, fvPatchField, volMesh >>& modes,
307 PtrList<GeometricField<tensor, fvPatchField, volMesh >>& modes2,
308 Eigen::VectorXd weights,
309 label Nmodes = 0,
310 bool consider_volumes);
311
312template<class Type, template<class> class PatchField, class GeoMesh>
313Eigen::VectorXd getMassMatrixFV(
314 GeometricField<Type, PatchField, GeoMesh>& snapshot)
315{
316 constexpr bool check_vol = std::is_same<volMesh, GeoMesh>::value
317 || std::is_same<surfaceMesh, GeoMesh>::value;
318
319 if constexpr(check_vol)
320 {
321 Eigen::MatrixXd snapEigen = Foam2Eigen::field2Eigen(snapshot);
322 label dim = std::nearbyint(snapEigen.rows() / (snapshot.mesh().V()).size());
323 Eigen::VectorXd volumes = Foam2Eigen::field2Eigen(snapshot.mesh().V());
324 Eigen::VectorXd vol3 = EigenFunctions::repeatElements(volumes, dim);
325 return vol3;
326 }
327 else if constexpr(std::is_same<pointMesh, GeoMesh>::value)
328 {
329 Eigen::MatrixXd snapEigen = Foam2Eigen::field2Eigen(snapshot);
330 label dim = std::nearbyint(snapEigen.rows() / (snapshot.mesh()
331 ().points()).size());
332 Eigen::VectorXd pointsdata = Foam2Eigen::field2Eigen(snapshot);
333 Eigen::VectorXd points3 = EigenFunctions::repeatElements(pointsdata, dim);
334 return points3;
335 }
336}
337
338template Eigen::VectorXd getMassMatrixFV(
339 GeometricField<scalar, fvPatchField, volMesh>& snapshot);
340template Eigen::VectorXd getMassMatrixFV(
341 GeometricField<vector, fvPatchField, volMesh>& snapshot);
342template Eigen::VectorXd getMassMatrixFV(
343 GeometricField<vector, pointPatchField, pointMesh>& snapshot);
344
345template Eigen::VectorXd getMassMatrixFV(
346 GeometricField<tensor, fvPatchField, volMesh>& snapshot);
347
348template<class Type, template<class> class PatchField, class GeoMesh>
349Eigen::VectorXd getCoeffs(GeometricField<Type, PatchField, GeoMesh>&
350 snapshot,
351 PtrList<GeometricField<Type, PatchField, GeoMesh >>& modes, label Nmodes,
352 bool consider_volumes)
353{
354 label Msize;
355
356 if (Nmodes == 0)
357 {
358 Msize = modes.size();
359 }
360 else
361 {
362 Msize = Nmodes;
363 }
364
365 M_Assert(modes.size() >= Msize,
366 "The Number of requested modes is larger then the available quantity.");
367 Eigen::MatrixXd F = Foam2Eigen::PtrList2Eigen(modes);
368 Eigen::MatrixXd M_matrix = getMassMatrix(modes, Nmodes, consider_volumes);
369 Eigen::MatrixXd snapEigen = Foam2Eigen::field2Eigen(snapshot);
370 Eigen::VectorXd a(Msize);
371 Eigen::VectorXd b(Msize);
372 constexpr bool check_vol = std::is_same<volMesh, GeoMesh>::value
373 || std::is_same<surfaceMesh, GeoMesh>::value;
374
375 if constexpr(check_vol)
376 {
377 Eigen::VectorXd V = getMassMatrixFV(modes[0]);
378 b = F.transpose().topRows(Msize) * V.asDiagonal() * snapEigen;
379 }
380 else if constexpr(std::is_same<pointMesh, GeoMesh>::value)
381 {
382 b = F.transpose().topRows(Msize) * snapEigen;
383 }
384
385 if (Pstream::parRun())
386 {
387 reduce(b, sumOp<Eigen::VectorXd>());
388 }
389
390 a = M_matrix.fullPivLu().solve(b);
391 return a;
392}
393
394template Eigen::VectorXd getCoeffs(
395 GeometricField<scalar, fvPatchField, volMesh>&
396 snapshot, PtrList<GeometricField<scalar, fvPatchField, volMesh >>& modes,
397 label Nmodes,
398 bool consider_volumes);
399
400template Eigen::VectorXd getCoeffs(
401 GeometricField<scalar, fvsPatchField, surfaceMesh>&
402 snapshot, PtrList<GeometricField<scalar, fvsPatchField, surfaceMesh >>& modes,
403 label Nmodes,
404 bool consider_volumes = false);
405
406template Eigen::VectorXd getCoeffs(
407 GeometricField<vector, fvPatchField, volMesh>&
408 snapshot, PtrList<GeometricField<vector, fvPatchField, volMesh >>& modes,
409 label Nmodes,
410 bool consider_volumes);
411template Eigen::VectorXd getCoeffs(
412 GeometricField<vector, pointPatchField, pointMesh>&
413 snapshot, PtrList<GeometricField<vector, pointPatchField, pointMesh>>& modes,
414 label Nmodes,
415 bool consider_volumes);
416
417template<class Type, template<class> class PatchField, class GeoMesh >
418Eigen::MatrixXd getCoeffs(PtrList<GeometricField<Type, PatchField, GeoMesh >>&
419 snapshots,
420 PtrList<GeometricField<Type, PatchField, GeoMesh >>& modes, label Nmodes,
421 bool consider_volumes)
422{
423 label Msize;
424
425 if (Nmodes == 0)
426 {
427 Msize = modes.size();
428 }
429 else
430 {
431 Msize = Nmodes;
432 }
433
434 M_Assert(modes.size() >= Msize,
435 "The Number of requested modes is larger then the available quantity.");
436 Eigen::MatrixXd coeff(Msize, snapshots.size());
437
438 for (auto i = 0; i < snapshots.size(); i++)
439 {
440 coeff.col(i) = getCoeffs(snapshots[i], modes, Nmodes, consider_volumes);
441 }
442
443 return coeff;
444}
445
446template Eigen::MatrixXd getCoeffs(
447 PtrList<GeometricField<scalar, fvPatchField, volMesh >>&
448 snapshot, PtrList<GeometricField<scalar, fvPatchField, volMesh >>& modes,
449 label Nmodes,
450 bool consider_volumes);
451
452template Eigen::MatrixXd getCoeffs(
453 PtrList<GeometricField<vector, fvPatchField, volMesh >>&
454 snapshot, PtrList<GeometricField<vector, fvPatchField, volMesh >>& modes,
455 label Nmodes,
456 bool consider_volumes);
457
458template Eigen::MatrixXd getCoeffs(
459 PtrList<GeometricField<scalar, fvsPatchField, surfaceMesh >>&
460 snapshot, PtrList<GeometricField<scalar, fvsPatchField, surfaceMesh >>& modes,
461 label Nmodes,
462 bool consider_volumes);
463template Eigen::MatrixXd getCoeffs(
464 PtrList<GeometricField<vector, pointPatchField, pointMesh>>&
465 snapshot, PtrList<GeometricField<vector, pointPatchField, pointMesh>>& modes,
466 label Nmodes,
467 bool consider_volumes);
468
469
470template<typename T>
471T project_to_POD_basis(T& field, PtrList<T>& modes)
472{
473 bool consider_volumes = true;
474 Eigen::MatrixXd coeffsField = getCoeffs(field, modes, modes.size(), consider_volumes);
475 PtrList<T> reducedField = reconstructFromCoeff(modes, coeffsField, modes.size());
476
477 T projField = field;
478 setToZero(projField);
479 addFields(projField, reducedField[0]);
480
481 return projField;
482}
483template volVectorField project_to_POD_basis(volVectorField& field, PtrList<volVectorField>& modes);
484template volScalarField project_to_POD_basis(volScalarField& field, PtrList<volScalarField>& modes);
485
486template<typename T>
487T project_to_POD_basis(T& field, PtrList<T>& modes, const T& meanField)
488{
489 T fieldCentered = field;
490 subtractFields(fieldCentered, meanField);
491
492 bool consider_volumes = true;
493 Eigen::MatrixXd coeffsField = getCoeffs(fieldCentered, modes, modes.size(), consider_volumes);
494 PtrList<T> reducedField = reconstructFromCoeff(modes, coeffsField, modes.size());
495
496 T projField = meanField;
497 addFields(projField, reducedField[0]);
498
499 return projField;
500}
501template volVectorField project_to_POD_basis(volVectorField& field, PtrList<volVectorField>& modes, const volVectorField& meanField);
502template volScalarField project_to_POD_basis(volScalarField& field, PtrList<volScalarField>& modes, const volScalarField& meanField);
503
504
505Eigen::MatrixXd parTimeCombMat(List<Eigen::VectorXd>
506 acquiredSnapshotsTimes,
507 Eigen::MatrixXd parameters)
508{
509 label parsNum = parameters.cols();
510 label parsSamplesNum = parameters.rows();
511 M_Assert(parsSamplesNum == acquiredSnapshotsTimes.size(),
512 "The list of time instants does not have the same number of vectors as the number of parameters samples");
513 Eigen::MatrixXd comb;
514 label totalSnapshotsNum = 0;
515
516 for (label k = 0; k < acquiredSnapshotsTimes.size(); k++)
517 {
518 totalSnapshotsNum += acquiredSnapshotsTimes[k].size();
519 }
520
521 comb.resize(totalSnapshotsNum, parsNum + 1);
522 label i = 0;
523
524 for (label j = 0; j < acquiredSnapshotsTimes.size(); j++)
525 {
526 for (label k = 0; k < acquiredSnapshotsTimes[j].size(); k++)
527 {
528 comb(i, parsNum) = (acquiredSnapshotsTimes[j])(k, 0);
529 comb.block(i, 0, 1, parsNum) = parameters.row(j);
530 i = i + 1;
531 }
532 }
533
534 return comb;
535}
536
537}
Header file of the ITHACAcoeffsMass file.
static Eigen::MatrixXd field2Eigen(GeometricField< Type, PatchField, GeoMesh > &field)
Convert a vector OpenFOAM field into an Eigen Vector.
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
Eigen::VectorXd repeatElements(Eigen::VectorXd input, int n)
Repeat each element n times after themselves.
Namespace to implement some useful assign operation of OF fields.
Eigen::VectorXd getCoeffs(GeometricField< Type, PatchField, GeoMesh > &snapshot, PtrList< GeometricField< Type, PatchField, GeoMesh > > &modes, label Nmodes, bool consider_volumes)
Projects a snapshot on a basis function and gets the coefficients of the projection.
PtrList< GeometricField< Type, PatchField, GeoMesh > > reconstructFromCoeff(PtrList< GeometricField< Type, PatchField, GeoMesh > > &modes, Eigen::MatrixXd &coeff_matrix, label Nmodes)
Exact reconstruction using a certain number of modes for a list of fields and the projection coeffici...
void setToZero(T &f1)
Set a field of type vol[Scalar|Vector|Tensor]Field to 0.
T project_to_POD_basis(T &field, PtrList< T > &modes)
Compute the projection of a field in a POD basis.
Eigen::VectorXd getMassMatrixFV(GeometricField< Type, PatchField, GeoMesh > &snapshot)
Gets a vector containing the volumes of each cell of the mesh.
Eigen::MatrixXd parTimeCombMat(List< Eigen::VectorXd > acquiredSnapshotsTimes, Eigen::MatrixXd parameters)
A method to compute the time-parameter combined matrix whose any single element corresponds to a uniq...
void addFields(T &f1, const T &f2c, double alpha)
Perform the following operation f1 + f2 * alpha with f1 and f2 being of type vol[Scalar|Vector|Tensor...
Eigen::MatrixXd getMassMatrix(PtrList< GeometricField< Type, PatchField, GeoMesh > > &modes, label Nmodes, bool consider_volumes)
Gets the mass matrix using the eigen routine.
void subtractFields(T &f1, const T &f2)
Perform the substraction (f1 - f2) between two fields of type vol[Scalar|Vector|Tensor]Field and alph...