Loading...
Searching...
No Matches
Modes.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
13License
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
31
34
35#include "Modes.H"
36
37template<class Type, template<class> class PatchField, class GeoMesh>
38List<Eigen::MatrixXd> Modes<Type, PatchField, GeoMesh>::toEigen()
39{
40 NBC = 0;
41
42 for (label i = 0; i < (this->first()).boundaryFieldRef().size(); i++)
43 {
44 if ((this->first()).boundaryFieldRef()[i].type() != "processor")
45 {
46 NBC++;
47 }
48 }
49
50 EigenModes.resize(NBC + 1);
51 EigenModes[0] = Foam2Eigen::PtrList2Eigen(this->toPtrList());
52 List<Eigen::MatrixXd> BC = Foam2Eigen::PtrList2EigenBC(this->toPtrList());
53
54 for (label i = 0; i < NBC; i++)
55 {
56 EigenModes[i + 1] = BC[i];
57 }
58
59 return EigenModes;
60}
61
62template<class Type, template<class> class PatchField, class GeoMesh>
63List<Eigen::MatrixXd> Modes<Type, PatchField, GeoMesh>::project(
64 fvMatrix<Type>& Af, label numberOfModes,
65 word projType)
66{
67 M_Assert(projType == "G" || projType == "PG",
68 "Projection type can be G for Galerkin or PG for Petrov-Galerkin");
69 List<Eigen::MatrixXd> LinSys;
70 LinSys.resize(2);
71
72 if (EigenModes.size() == 0)
73 {
74 toEigen();
75 }
76
77 Eigen::SparseMatrix<double> Ae;
78 Eigen::VectorXd be;
79 Foam2Eigen::fvMatrix2Eigen(Af, Ae, be);
80
81 if (numberOfModes == 0)
82 {
83 if (projType == "G")
84 {
85 LinSys[0] = EigenModes[0].transpose() * Ae * EigenModes[0];
86 LinSys[1] = EigenModes[0].transpose() * be;
87 }
88
89 if (projType == "PG")
90 {
91 LinSys[0] = (Ae * EigenModes[0]).transpose() * Ae * EigenModes[0];
92 LinSys[1] = (Ae * EigenModes[0]).transpose() * be;
93 }
94 }
95 else
96 {
97 M_Assert(numberOfModes <= EigenModes[0].cols(),
98 "Number of required modes for projection is higher then the number of available ones");
99
100 if (projType == "G")
101 {
102 LinSys[0] = ((EigenModes[0]).leftCols(numberOfModes)).transpose() * Ae *
103 (EigenModes[0]).leftCols(numberOfModes);
104 LinSys[1] = ((EigenModes[0]).leftCols(numberOfModes)).transpose() * be;
105 }
106
107 if (projType == "PG")
108 {
109 LinSys[0] = (Ae * ((EigenModes[0]).leftCols(numberOfModes))).transpose() * Ae *
110 (EigenModes[0]).leftCols(numberOfModes);
111 LinSys[1] = (Ae * ((EigenModes[0]).leftCols(numberOfModes))).transpose() * be;
112 }
113 }
114
115 return LinSys;
116}
117
118template<class Type, template<class> class PatchField, class GeoMesh>
119Eigen::MatrixXd Modes<Type, PatchField, GeoMesh>::project(
120 GeometricField<Type, PatchField, GeoMesh>&
121 field, label numberOfModes, word projType, fvMatrix<Type>* Af)
122{
123 M_Assert(projType == "F" || projType == "G" || projType == "PG",
124 "Projection type can be F for Frobenius, G for Galerkin or PG for Petrov-Galerkin");
125 Eigen::MatrixXd fieldEig = Foam2Eigen::field2Eigen(field);
126 auto vol = ITHACAutilities::getMassMatrixFV(field);
127 Eigen::MatrixXd projField;
128 Eigen::MatrixXd M;
129 Eigen::MatrixXd b;
130
131 if (EigenModes.size() == 0)
132 {
133 toEigen();
134 }
135
136 if (numberOfModes == 0)
137 {
138 if (projType == "F")
139 {
140 vol = Eigen::VectorXd::Ones(vol.size());
141 b = EigenModes[0].transpose() * vol.asDiagonal() * fieldEig;
142 M = EigenModes[0].transpose() * EigenModes[0].transpose();
143 projField = M.fullPivLu().solve(b);
144 }
145 else if (projType == "G")
146 {
147 projField = EigenModes[0].transpose() * vol.asDiagonal() * fieldEig;
148 }
149 else if (projType == "PG")
150 {
151 M_Assert(Af != NULL,
152 "Using a Petrov-Galerkin projection you have to provide also the system matrix");
153 Eigen::SparseMatrix<double> Ae;
154 Eigen::VectorXd be;
155 Foam2Eigen::fvMatrix2Eigen(* Af, Ae, be);
156 projField = (Ae * EigenModes[0]).transpose() * vol.asDiagonal() * fieldEig;
157 }
158 }
159 else
160 {
161 M_Assert(numberOfModes <= EigenModes[0].cols(),
162 "Number of required modes for projection is higher then the number of available ones");
163
164 if (projType == "F")
165 {
166 vol = Eigen::VectorXd::Ones(vol.size());
167 M = EigenModes[0].leftCols(numberOfModes).transpose() * EigenModes[0].leftCols(
168 numberOfModes);
169 b = ((EigenModes[0]).leftCols(numberOfModes)).transpose() *
170 vol.asDiagonal() * fieldEig;
171 projField = M.fullPivLu().solve(b);
172 }
173 else if (projType == "G")
174 {
175 projField = ((EigenModes[0]).leftCols(numberOfModes)).transpose() *
176 vol.asDiagonal() * fieldEig;
177 }
178 else if (projType == "PG")
179 {
180 M_Assert(Af != NULL,
181 "Using a Petrov-Galerkin projection you have to provide also the system matrix");
182 Eigen::SparseMatrix<double> Ae;
183 Eigen::VectorXd be;
184 Foam2Eigen::fvMatrix2Eigen(* Af, Ae, be);
185 projField = (Ae * ((EigenModes[0]).leftCols(numberOfModes))).transpose() *
186 vol.asDiagonal() * fieldEig;
187 }
188 }
189
190 return projField;
191}
192
193template<class Type, template<class> class PatchField, class GeoMesh>
194GeometricField<Type, PatchField, GeoMesh>
195Modes<Type, PatchField, GeoMesh>::projectSnapshot(
196 GeometricField<Type, PatchField, GeoMesh>&
197 field, label numberOfModes, word projType, fvMatrix<Type>* Af)
198{
199 Eigen::MatrixXd proj = project(field, numberOfModes, projType, Af);
200 GeometricField<Type, PatchField, GeoMesh> projSnap = field;
201 reconstruct(projSnap, proj, projSnap.name());
202 return projSnap;
203}
204
205template<class Type, template<class> class PatchField, class GeoMesh>
206Eigen::MatrixXd Modes<Type, PatchField, GeoMesh>::project(
207 PtrList<GeometricField<Type, PatchField, GeoMesh >> &
208 fields,
209 label numberOfModes, word projType, fvMatrix<Type>* Af)
210{
211 M_Assert(projType == "G" || projType == "PG",
212 "Projection type can be G for Galerking or PG for Petrov-Galerkin");
213 Eigen::MatrixXd fieldEig = Foam2Eigen::PtrList2Eigen(fields);
214 auto vol = ITHACAutilities::getMassMatrixFV(fields[0]);
215 Eigen::MatrixXd projField;
216
217 if (EigenModes.size() == 0)
218 {
219 toEigen();
220 }
221
222 if (numberOfModes == 0)
223 {
224 if (projType == "F")
225 {
226 vol = Eigen::VectorXd::Ones(vol.size());
227 projField = EigenModes[0].transpose() * vol.asDiagonal() * fieldEig;
228 }
229 else if (projType == "G")
230 {
231 projField = EigenModes[0].transpose() * vol.asDiagonal() * fieldEig;
232 }
233 else if (projType == "PG")
234 {
235 M_Assert(Af != NULL,
236 "Using a Petrov-Galerkin projection you have to provide also the system matrix");
237 Eigen::SparseMatrix<double> Ae;
238 Eigen::VectorXd be;
239 Foam2Eigen::fvMatrix2Eigen(* Af, Ae, be);
240 projField = (Ae * EigenModes[0]).transpose() * vol.asDiagonal() * fieldEig;
241 }
242 }
243 else
244 {
245 M_Assert(numberOfModes <= EigenModes[0].cols(),
246 "Number of required modes for projection is higher then the number of available ones");
247
248 if (projType == "F")
249 {
250 vol = Eigen::VectorXd::Ones(vol.size());
251 projField = ((EigenModes[0]).leftCols(numberOfModes)).transpose() *
252 vol.asDiagonal() * fieldEig;
253 }
254 else if (projType == "G")
255 {
256 projField = ((EigenModes[0]).leftCols(numberOfModes)).transpose() *
257 vol.asDiagonal() * fieldEig;
258 }
259 else if (projType == "PG")
260 {
261 M_Assert(Af != NULL,
262 "Using a Petrov-Galerkin projection you have to provide also the system matrix");
263 Eigen::SparseMatrix<double> Ae;
264 Eigen::VectorXd be;
265 Foam2Eigen::fvMatrix2Eigen(* Af, Ae, be);
266 projField = (Ae * ((EigenModes[0]).leftCols(numberOfModes))).transpose() *
267 vol.asDiagonal() * fieldEig;
268 }
269 }
270
271 return projField;
272}
273
274template<class Type, template<class> class PatchField, class GeoMesh>
275GeometricField<Type, PatchField, GeoMesh>
276Modes<Type, PatchField, GeoMesh>::reconstruct(
277 GeometricField<Type, PatchField, GeoMesh>& inputField,
278 Eigen::MatrixXd Coeff,
279 word Name)
280{
281 if (EigenModes.size() == 0)
282 {
283 toEigen();
284 }
285
286 label Nmodes = Coeff.rows();
287 Eigen::VectorXd InField = EigenModes[0].leftCols(Nmodes) * Coeff;
288
289 if (inputField.name() == "nut")
290 {
291 InField = (InField.array() < 0).select(0, InField);
292 }
293
294 inputField = Foam2Eigen::Eigen2field(inputField, InField);
295 inputField.rename(Name);
296
297 for (label i = 0; i < NBC; i++)
298 {
299 Eigen::VectorXd BF = EigenModes[i + 1].leftCols(Nmodes) * Coeff;
300
301 if (inputField.name() == "nut")
302 {
303 BF = (BF.array() < 0).select(0, BF);
304 }
305
306 ITHACAutilities::assignBC(inputField, i, BF);
307 }
308
309 return inputField;
310}
311
312template<class Type, template<class> class PatchField, class GeoMesh>
313PtrList<GeometricField<Type, PatchField, GeoMesh >>
314Modes<Type, PatchField, GeoMesh>::reconstruct(
315 GeometricField<Type, PatchField, GeoMesh>& inputField,
316 List < Eigen::MatrixXd> Coeff,
317 word Name)
318{
319 PtrList<GeometricField<Type, PatchField, GeoMesh >> inputFields;
320 inputFields.resize(0);
321 for (label i = 0; i < Coeff.size(); i++)
322 {
323 inputField = reconstruct(inputField, Coeff[i], Name);
324 inputFields.append(inputField.clone());
325 }
326
327 return inputFields;
328}
329
330
331template<class Type, template<class> class PatchField, class GeoMesh >
332void Modes<Type, PatchField, GeoMesh>::projectSnapshots(
333 PtrList<GeometricField<Type, PatchField, GeoMesh >> snapshots,
334 PtrList<GeometricField<Type, PatchField, GeoMesh >>& projSnapshots,
335 PtrList<volScalarField> Volumes,
336 label numberOfModes,
337 word innerProduct)
338{
339 if (EigenModes.size() == 0)
340 {
341 toEigen();
342 }
343
344 M_Assert(snapshots.size() == Volumes.size(),
345 "The number of snapshots and the number of volumes vectors must be equal");
346 M_Assert(numberOfModes <= this->size(),
347 "The number of Modes used for the projection cannot be bigger than the number of available modes");
348 M_Assert(innerProduct == "L2" || innerProduct == "Frobenius",
349 "The chosen inner product is not implemented");
350 projSnapshots.resize(snapshots.size());
351 label dim = std::nearbyint(EigenModes[0].rows() /
352 Volumes[0].size()); //Checking if volumes and modes have the same size that means check if the problem is vector or scalar
353 Eigen::MatrixXd totVolumes(Volumes[0].size() * dim, Volumes.size());
354
355 for (label i = 0; i < Volumes.size(); i++)
356 {
357 totVolumes.col(i) = Foam2Eigen::field2Eigen(Volumes[i]);
358 }
359
360 totVolumes.replicate(dim, 1);
361 Eigen::MatrixXd Modes;
362
363 if (numberOfModes == 0)
364 {
365 Modes = EigenModes[0];
366 }
367 else
368 {
369 Modes = EigenModes[0].leftCols(numberOfModes);
370 }
371
372 Eigen::MatrixXd M;
373 Eigen::MatrixXd projSnapI;
374 Eigen::MatrixXd projSnapCoeff;
375
376 for (label i = 0; i < snapshots.size(); i++)
377 {
378 GeometricField<Type, PatchField, GeoMesh> Fr = snapshots[0];
379 Eigen::MatrixXd F_eigen = Foam2Eigen::field2Eigen(snapshots[i]);
380
381 if (innerProduct == "L2")
382 {
383 M = Modes.transpose() * (totVolumes.col(i)).asDiagonal() * Modes;
384 projSnapI = Modes.transpose() * (totVolumes.col(i)).asDiagonal() * F_eigen;
385 }
386 else //Frobenius
387 {
388 M = Modes.transpose() * Modes;
389 projSnapI = Modes.transpose() * F_eigen;
390 }
391
392 projSnapCoeff = M.fullPivLu().solve(projSnapI);
393 reconstruct(Fr, projSnapCoeff, "projSnap");
394 projSnapshots.set(i, Fr.clone());
395 }
396}
397
398template<class Type, template<class> class PatchField, class GeoMesh>
399void Modes<Type, PatchField, GeoMesh>::projectSnapshots(
400 PtrList<GeometricField<Type, PatchField, GeoMesh >> snapshots,
401 PtrList<GeometricField<Type, PatchField, GeoMesh >>& projSnapshots,
402 PtrList<volScalarField> Volumes, word innerProduct)
403{
404 label numberOfModes = 0;
405 projectSnapshots(snapshots, projSnapshots, Volumes, numberOfModes,
406 innerProduct);
407}
408
409template<class Type, template<class> class PatchField, class GeoMesh>
410void Modes<Type, PatchField, GeoMesh>::projectSnapshots(
411 PtrList<GeometricField<Type, PatchField, GeoMesh >> snapshots,
412 PtrList<GeometricField<Type, PatchField, GeoMesh >>& projSnapshots,
413 label numberOfModes,
414 word innerProduct)
415{
416 if (EigenModes.size() == 0)
417 {
418 toEigen();
419 }
420
421 M_Assert(numberOfModes <= this->size(),
422 "The number of Modes used for the projection cannot be bigger than the number of available modes");
423 M_Assert(innerProduct == "L2" || innerProduct == "Frobenius",
424 "The chosen inner product is not implemented");
425 projSnapshots.resize(snapshots.size());
426 Eigen::MatrixXd Modes;
427
428 if (numberOfModes == 0)
429 {
430 Modes = EigenModes[0];
431 }
432 else
433 {
434 Modes = EigenModes[0].leftCols(numberOfModes);
435 }
436
437 Eigen::MatrixXd M_vol;
438 Eigen::MatrixXd M;
439 Eigen::MatrixXd projSnapI;
440 Eigen::MatrixXd projSnapCoeff;
441
442 for (label i = 0; i < snapshots.size(); i++)
443 {
444 GeometricField<Type, PatchField, GeoMesh> Fr = snapshots[0];
445 Eigen::MatrixXd F_eigen = Foam2Eigen::field2Eigen(snapshots[i]);
446
447 if (innerProduct == "L2")
448 {
449 M_vol = ITHACAutilities::getMassMatrixFV(snapshots[i]);
450 }
451 else if (innerProduct == "Frobenius")
452 {
453 M_vol = Eigen::VectorXd::Identity(F_eigen.rows(), 1);
454 }
455 else
456 {
457 std::cout << "Inner product not defined" << endl;
458 exit(0);
459 }
460
461 M = Modes.transpose() * M_vol.asDiagonal() * Modes;
462 projSnapI = Modes.transpose() * M_vol.asDiagonal() * F_eigen;
463 projSnapCoeff = M.fullPivLu().solve(projSnapI);
464 reconstruct(Fr, projSnapCoeff, "projSnap");
465 projSnapshots.set(i, Fr.clone());
466 }
467}
468
469template<class Type, template<class> class PatchField, class GeoMesh>
470void Modes<Type, PatchField, GeoMesh>::projectSnapshots(
471 PtrList<GeometricField<Type, PatchField, GeoMesh >> snapshots,
472 PtrList<GeometricField<Type, PatchField, GeoMesh >>& projSnapshots,
473 word innerProduct)
474{
475 label numberOfModes = 0;
476 projectSnapshots(snapshots, projSnapshots, numberOfModes, innerProduct);
477}
478
479template<class Type, template<class> class PatchField, class GeoMesh>
480void Modes<Type, PatchField, GeoMesh>::operator=(const
481 PtrList<GeometricField<Type, PatchField, GeoMesh >> & modes)
482{
483 this->resize(modes.size());
484 for (label i = 0; i < modes.size(); i++)
485 {
486 (* this).set(i, modes[i].clone());
487 }
488}
489
490
491template class Modes<scalar, fvPatchField, volMesh>;
492template class Modes<vector, fvPatchField, volMesh>;
493template class Modes<scalar, fvsPatchField, surfaceMesh>;
M_Assert
#define M_Assert(Expr, Msg)
Definition ITHACAassert.H:46
Modes.H
Header file of the Modes class.
Foam2Eigen::Eigen2field
static GeometricField< scalar, PatchField, GeoMesh > Eigen2field(GeometricField< scalar, PatchField, GeoMesh > &field, Eigen::VectorXd &eigen_vector, bool correctBC=true)
Convert a vector in Eigen format into an OpenFOAM scalar GeometricField.
Definition Foam2Eigen.C:514
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::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
Foam2Eigen::fvMatrix2Eigen
static void fvMatrix2Eigen(fvMatrix< type_foam_matrix > foam_matrix, type_A &A, type_B &b)
Convert a FvMatrix OpenFOAM matrix (Linear System) into a Eigen Matrix A and a source vector b.
Modes
Implementation of a container class derived from PtrList.
Definition Modes.H:69
Modes< scalar, fvPatchField, volMesh >::reconstruct
GeometricField< scalar, fvPatchField, volMesh > reconstruct(GeometricField< scalar, fvPatchField, volMesh > &inputField, Eigen::MatrixXd Coeff, word Name)
Modes::EigenModes
List< Eigen::MatrixXd > EigenModes
List of Matrices that contains the internalField and the additional matrices for the boundary patches...
Definition Modes.H:73
Modes::toEigen
List< Eigen::MatrixXd > toEigen()
Method that convert a PtrList of modes into Eigen matrices filling the EigenModes object.
Definition Modes.C:38
Modes::NBC
label NBC
Number of patches.
Definition Modes.H:76
Modes::project
List< Eigen::MatrixXd > project(fvMatrix< Type > &Af, label numberOfModes=0, word projType="G")
A function that project an FvMatrix (OpenFoam linear System) on the modes.
Definition Modes.C:63
Modes::projectSnapshots
void projectSnapshots(PtrList< GeometricField< Type, PatchField, GeoMesh > > snapshots, PtrList< GeometricField< Type, PatchField, GeoMesh > > &projSnapshots, label numberOfModes=0, word innerProduct="L2")
Function to project a list of fields into the modes manifold.
Definition Modes.C:410
Modes::toPtrList
PtrList< GeometricField< Type, PatchField, GeoMesh > > & toPtrList()
Function that returns the Modes object as a standard PtrList.
Definition Modes.H:86
Modes::projectSnapshot
GeometricField< Type, PatchField, GeoMesh > projectSnapshot(GeometricField< Type, PatchField, GeoMesh > &field, label numberOfModes=0, word projType="G", fvMatrix< Type > *Af=NULL)
A function that project and reconstruct a snapshot on the modes.
Definition Modes.C:195
Modes::operator=
void operator=(const PtrList< GeometricField< Type, PatchField, GeoMesh > > &modes)
Definition Modes.C:480
resize
Tf resize(coldSideSize)
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::assignBC
void assignBC(GeometricField< scalar, fvPatchField, volMesh > &s, label BC_ind, double value)
Assign uniform Boundary Condition to a volScalarField.
Definition ITHACAassign.C:143
i
label i
Definition pEqn.H:46
Generated by   doxygen 1.13.2