Loading...
Searching...
No Matches
ITHACAcoeffsMass.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-------------------------------------------------------------------------------
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 for (label k = 0; k < coeff_matrix.cols(); k++)
38 {
39 for (label i = 0; i < Nmodes; i++)
40 if (i == 0)
41 {
42 rec_field.append(modes[i] * coeff_matrix(i, k));
43 }
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 >>
55reconstructFromCoeff(
56 PtrList<GeometricField<scalar, fvPatchField, volMesh >> & modes,
57 Eigen::MatrixXd& coeff_matrix, label Nmodes);
58template PtrList<GeometricField<vector, fvPatchField, volMesh >>
59reconstructFromCoeff(
60 PtrList<GeometricField<vector, fvPatchField, volMesh >> & modes,
61 Eigen::MatrixXd& coeff_matrix, label Nmodes);
62template PtrList<GeometricField<tensor, fvPatchField, volMesh >>
63reconstructFromCoeff(
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 if (Nmodes == 0)
74 {
75 Msize = modes.size();
76 }
77
78 else
79 {
80 Msize = Nmodes;
81 }
82 M_Assert(modes.size() >= Msize,
83 "The Number of requested modes is larger then the available quantity.");
84 Eigen::MatrixXd F = Foam2Eigen::PtrList2Eigen(modes);
85 Eigen::MatrixXd M = Eigen::MatrixXd::Zero(Msize, Msize);
86
87 if (consider_volumes)
88 {
89 Eigen::VectorXd V = getMassMatrixFV(modes[0]);
90
91 // If V is too big then we need a for loop to avoid RAM overload
92 if (V.size() > 1e6)
93 {
94 for (int i = 0; i < V.size(); i++)
95 {
96 M += V(i) * F.transpose().topRows(Msize).col(i) * F.leftCols(Msize).row(i);
97 }
98 }
99 else
100 {
101 M = F.transpose().topRows(Msize) * V.asDiagonal() * F.leftCols(Msize);
102 }
103 }
104 else
105 {
106 M = F.transpose().topRows(Msize) * F.leftCols(Msize);
107 }
108
109 if (Pstream::parRun())
110 {
111 reduce(M, sumOp<Eigen::MatrixXd>());
112 }
113
114 return M;
115}
116
117template Eigen::MatrixXd getMassMatrix(
118 PtrList<GeometricField<scalar, fvPatchField, volMesh >> & modes, label Nmodes,
119 bool consider_volumes);
120
121template Eigen::MatrixXd getMassMatrix(
122 PtrList<GeometricField<scalar, fvsPatchField, surfaceMesh >>& modes,
123 label Nmodes,
124 bool consider_volumes = false);
125
126template Eigen::MatrixXd getMassMatrix(
127 PtrList<GeometricField<vector, fvPatchField, volMesh >> & modes, label Nmodes,
128 bool consider_volumes);
129
130template Eigen::MatrixXd getMassMatrix(
131 PtrList<GeometricField<tensor, fvPatchField, volMesh >>& modes, label Nmodes,
132 bool consider_volumes);
133
134template<class Type, template<class> class PatchField, class GeoMesh>
135Eigen::MatrixXd getMassMatrix(
136 PtrList<GeometricField<Type, PatchField, GeoMesh >> & modes,
137 PtrList<GeometricField<Type, PatchField, GeoMesh >>& modes2, label Nmodes,
138 bool consider_volumes)
139{
140 label Msize, Msize2;
141
142 if (Nmodes == 0)
143 {
144 Msize = modes.size();
145 Msize2 = modes2.size();
146 }
147 else
148 {
149 Msize = Nmodes;
150 Msize2 = Nmodes;
151 }
152
153 M_Assert(modes.size() >= Msize,
154 "The Number of requested modes is larger then the available quantity.");
155 M_Assert(modes2.size() >= Msize2,
156 "The Number of requested modes is larger then the available quantity.");
157 Eigen::MatrixXd F = Foam2Eigen::PtrList2Eigen(modes);
158 Eigen::MatrixXd F2 = Foam2Eigen::PtrList2Eigen(modes2);
159 Eigen::MatrixXd M = Eigen::MatrixXd::Zero(Msize, Msize2);
160
161 if (consider_volumes)
162 {
163 Eigen::VectorXd V = getMassMatrixFV(modes[0]);
164
165 // If V is too big then we need a for loop to avoid RAM overload
166 if (V.size() > 1e6)
167 {
168 for (int i = 0; i < V.size(); i++)
169 {
170 M += V(i) * F.transpose().topRows(Msize).col(i) * F2.leftCols(Msize2).row(i);
171 }
172 }
173 else
174 {
175 // Classical M computation
176 M = F.transpose().topRows(Msize) * V.asDiagonal() * F2.leftCols(Msize2);
177 }
178 }
179 else
180 {
181 M = F.transpose().topRows(Msize) * F2.leftCols(Msize2);
182 }
183
184 if (Pstream::parRun())
185 {
186 reduce(M, sumOp<Eigen::MatrixXd>());
187 }
188
189 return M;
190}
191
192template Eigen::MatrixXd getMassMatrix(
193 PtrList<GeometricField<scalar, fvPatchField, volMesh >> & modes,
194 PtrList<GeometricField<scalar, fvPatchField, volMesh >>& modes2, label Nmodes,
195 bool consider_volumes);
196
197template Eigen::MatrixXd getMassMatrix(
198 PtrList<GeometricField<scalar, fvsPatchField, surfaceMesh >> & modes,
199 PtrList<GeometricField<scalar, fvsPatchField, surfaceMesh >>& modes2,
200 label Nmodes,
201 bool consider_volumes = false);
202
203template Eigen::MatrixXd getMassMatrix(
204 PtrList<GeometricField<vector, fvPatchField, volMesh >> & modes,
205 PtrList<GeometricField<vector, fvPatchField, volMesh >>& modes2, label Nmodes,
206 bool consider_volumes);
207
208template Eigen::MatrixXd getMassMatrix(
209 PtrList<GeometricField<tensor, fvPatchField, volMesh >> & modes,
210 PtrList<GeometricField<tensor, fvPatchField, volMesh >>& modes2, label Nmodes,
211 bool consider_volumes);
212
213
214template<class Type, template<class> class PatchField, class GeoMesh>
215Eigen::MatrixXd getMassMatrix(
216 PtrList<GeometricField<Type, PatchField, GeoMesh >>
217 & modes,
218 PtrList<GeometricField<Type, PatchField, GeoMesh >>& modes2,
219 Eigen::VectorXd weights,
220 label Nmodes,
221 bool consider_volumes)
222{
223 label Msize, Msize2;
224
225 if (Nmodes == 0)
226 {
227 Msize = modes.size();
228 Msize2 = modes2.size();
229 }
230 else
231 {
232 Msize = Nmodes;
233 Msize2 = Nmodes;
234 }
235
236 M_Assert(modes.size() >= Msize,
237 "The Number of requested modes is larger then the available quantity.");
238 M_Assert(modes2.size() >= Msize2,
239 "The Number of requested modes is larger then the available quantity.");
240 Eigen::MatrixXd F = Foam2Eigen::PtrList2Eigen(modes);
241 Eigen::MatrixXd F2 = Foam2Eigen::PtrList2Eigen(modes2);
242 Eigen::MatrixXd M = Eigen::MatrixXd::Zero(Msize, Msize2);
243
244 if (consider_volumes)
245 {
246 Eigen::VectorXd V = getMassMatrixFV(modes[0]);
247
248 // If V is too big then we need a for loop to avoid RAM overload
249 if (V.size() > 1e6)
250 {
251 for (int i = 0; i < V.size(); i++)
252 {
253 M += (V(i) * weights(i)) * F.transpose().topRows(Msize).col(i) * F2.leftCols(
254 Msize2).row(i);
255 }
256 }
257 else
258 {
259 // Classical M computation
260 M = F.transpose().topRows(Msize) * V.asDiagonal() * ( weights.asDiagonal()
261 * F2.leftCols(Msize2) );
262 }
263 }
264 else
265 {
266 M = F.transpose().topRows(Msize) * weights.asDiagonal() * F2.leftCols(Msize2);
267 }
268
269 if (Pstream::parRun())
270 {
271 reduce(M, sumOp<Eigen::MatrixXd>());
272 }
273
274 return M;
275}
276
277template Eigen::MatrixXd getMassMatrix(
278 PtrList<GeometricField<scalar, fvPatchField, volMesh >> & modes,
279 PtrList<GeometricField<scalar, fvPatchField, volMesh >>& modes2,
280 Eigen::VectorXd weights,
281 label Nmodes = 0,
282 bool consider_volumes);
283
284template Eigen::MatrixXd getMassMatrix(
285 PtrList<GeometricField<scalar, fvsPatchField, surfaceMesh >> & modes,
286 PtrList<GeometricField<scalar, fvsPatchField, surfaceMesh >>& modes2,
287 Eigen::VectorXd weights,
288 label Nmodes = 0,
289 bool consider_volumes);
290
291template Eigen::MatrixXd getMassMatrix(
292 PtrList<GeometricField<vector, fvPatchField, volMesh >> & modes,
293 PtrList<GeometricField<vector, fvPatchField, volMesh >>& modes2,
294 Eigen::VectorXd weights,
295 label Nmodes = 0,
296 bool consider_volumes);
297
298template Eigen::MatrixXd getMassMatrix(
299 PtrList<GeometricField<tensor, fvPatchField, volMesh >> & modes,
300 PtrList<GeometricField<tensor, fvPatchField, volMesh >>& modes2,
301 Eigen::VectorXd weights,
302 label Nmodes = 0,
303 bool consider_volumes);
304
305
306template<class Type, template<class> class PatchField, class GeoMesh>
307Eigen::VectorXd getMassMatrixFV(
308 GeometricField<Type, PatchField, GeoMesh>& snapshot)
309{
310 Eigen::MatrixXd snapEigen = Foam2Eigen::field2Eigen(snapshot);
311 label dim = std::nearbyint(snapEigen.rows() / (snapshot.mesh().V()).size());
312 Eigen::VectorXd volumes = Foam2Eigen::field2Eigen(snapshot.mesh().V());
313 Eigen::VectorXd vol3 = volumes.replicate(dim, 1);
314 return vol3;
315}
316
317template Eigen::VectorXd getMassMatrixFV(
318 GeometricField<scalar, fvPatchField, volMesh>& snapshot);
319template Eigen::VectorXd getMassMatrixFV(
320 GeometricField<vector, fvPatchField, volMesh>& snapshot);
321template Eigen::VectorXd getMassMatrixFV(
322 GeometricField<tensor, fvPatchField, volMesh>& snapshot);
323
324template<class Type, template<class> class PatchField, class GeoMesh>
325Eigen::VectorXd getCoeffs(GeometricField<Type, PatchField, GeoMesh>&
326 snapshot,
327 PtrList<GeometricField<Type, PatchField, GeoMesh >> & modes, label Nmodes,
328 bool consider_volumes)
329{
330 label Msize;
331 if (Nmodes == 0)
332 {
333 Msize = modes.size();
334 }
335
336 else
337 {
338 Msize = Nmodes;
339 }
340 M_Assert(modes.size() >= Msize,
341 "The Number of requested modes is larger then the available quantity.");
342 Eigen::MatrixXd F = Foam2Eigen::PtrList2Eigen(modes);
343 Eigen::MatrixXd M_matrix = getMassMatrix(modes, Nmodes, consider_volumes);
344 Eigen::MatrixXd snapEigen = Foam2Eigen::field2Eigen(snapshot);
345 Eigen::VectorXd a(Msize);
346 Eigen::VectorXd b(Msize);
347
348 if (consider_volumes)
349 {
350 Eigen::VectorXd V = getMassMatrixFV(modes[0]);
351 b = F.transpose().topRows(Msize) * V.asDiagonal() * snapEigen;
352 }
353 else
354 {
355 b = F.transpose().topRows(Msize) * snapEigen;
356 }
357
358 if (Pstream::parRun())
359 {
360 reduce(b, sumOp<Eigen::VectorXd>());
361 }
362
363 a = M_matrix.fullPivLu().solve(b);
364 return a;
365}
366
367template Eigen::VectorXd getCoeffs(
368 GeometricField<scalar, fvPatchField, volMesh>&
369 snapshot, PtrList<GeometricField<scalar, fvPatchField, volMesh >> & modes,
370 label Nmodes,
371 bool consider_volumes);
372
373template Eigen::VectorXd getCoeffs(
374 GeometricField<scalar, fvsPatchField, surfaceMesh>&
375 snapshot, PtrList<GeometricField<scalar, fvsPatchField, surfaceMesh >> & modes,
376 label Nmodes,
377 bool consider_volumes = false);
378
379template Eigen::VectorXd getCoeffs(
380 GeometricField<vector, fvPatchField, volMesh>&
381 snapshot, PtrList<GeometricField<vector, fvPatchField, volMesh >> & modes,
382 label Nmodes,
383 bool consider_volumes);
384
385template<class Type, template<class> class PatchField, class GeoMesh >
386Eigen::MatrixXd getCoeffs(PtrList<GeometricField<Type, PatchField, GeoMesh >> &
387 snapshots,
388 PtrList<GeometricField<Type, PatchField, GeoMesh >>& modes, label Nmodes,
389 bool consider_volumes)
390{
391 label Msize;
392
393 if (Nmodes == 0)
394 {
395 Msize = modes.size();
396 }
397 else
398 {
399 Msize = Nmodes;
400 }
401
402 M_Assert(modes.size() >= Msize,
403 "The Number of requested modes is larger then the available quantity.");
404 Eigen::MatrixXd coeff(Msize, snapshots.size());
405
406 for (auto i = 0; i < snapshots.size(); i++)
407 {
408 coeff.col(i) = getCoeffs(snapshots[i], modes, Nmodes, consider_volumes);
409 }
410
411 return coeff;
412}
413
414template Eigen::MatrixXd getCoeffs(
415 PtrList<GeometricField<scalar, fvPatchField, volMesh >> &
416 snapshot, PtrList<GeometricField<scalar, fvPatchField, volMesh >>& modes,
417 label Nmodes,
418 bool consider_volumes);
419
420template Eigen::MatrixXd getCoeffs(
421 PtrList<GeometricField<vector, fvPatchField, volMesh >> &
422 snapshot, PtrList<GeometricField<vector, fvPatchField, volMesh >>& modes,
423 label Nmodes,
424 bool consider_volumes);
425
426template Eigen::MatrixXd getCoeffs(
427 PtrList<GeometricField<scalar, fvsPatchField, surfaceMesh >> &
428 snapshot, PtrList<GeometricField<scalar, fvsPatchField, surfaceMesh >>& modes,
429 label Nmodes,
430 bool consider_volumes);
431
432Eigen::MatrixXd parTimeCombMat(List<Eigen::VectorXd>
433 acquiredSnapshotsTimes,
434 Eigen::MatrixXd parameters)
435{
436 label parsNum = parameters.cols();
437 label parsSamplesNum = parameters.rows();
438 M_Assert(parsSamplesNum == acquiredSnapshotsTimes.size(),
439 "The list of time instants does not have the same number of vectors as the number of parameters samples");
440 Eigen::MatrixXd comb;
441 label totalSnapshotsNum = 0;
442
443 for (label k = 0; k < acquiredSnapshotsTimes.size(); k++)
444 {
445 totalSnapshotsNum += acquiredSnapshotsTimes[k].size();
446 }
447
448 comb.resize(totalSnapshotsNum, parsNum + 1);
449 label i = 0;
450
451 for (label j = 0; j < acquiredSnapshotsTimes.size(); j++)
452 {
453 for (label k = 0; k < acquiredSnapshotsTimes[j].size(); k++)
454 {
455 comb(i, parsNum) = (acquiredSnapshotsTimes[j])(k, 0);
456 comb.block(i, 0, 1, parsNum) = parameters.row(j);
457 i = i + 1;
458 }
459 }
460
461 return comb;
462}
463
464}
M_Assert
#define M_Assert(Expr, Msg)
Definition ITHACAassert.H:46
ITHACAcoeffsMass.H
Header file of the ITHACAcoeffsMass file.
Foam2Eigen::field2Eigen
static Eigen::VectorXd field2Eigen(GeometricField< tensor, PatchField, GeoMesh > &field)
Convert a vector OpenFOAM field into an Eigen Vector.
Definition Foam2Eigen.C:40
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
ITHACAutilities
Namespace to implement some useful assign operation of OF fields.
Definition ITHACAassign.C:28
ITHACAutilities::getCoeffs
Eigen::VectorXd getCoeffs(GeometricField< Type, PatchField, GeoMesh > &snapshot, PtrList< GeometricField< Type, PatchField, GeoMesh > > &modes, label Nmodes, bool consider_volumes)
Projects a snapshot on a basis function and gets the coefficients of the projection.
Definition ITHACAcoeffsMass.C:325
ITHACAutilities::reconstructFromCoeff
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...
Definition ITHACAcoeffsMass.C:31
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::parTimeCombMat
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...
Definition ITHACAcoeffsMass.C:432
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
i
label i
Definition pEqn.H:46
Generated by   doxygen 1.13.2