Skip to content

Commit d788f82

Browse files
giovastabilegiovastabile
committed
REFACTOR: moved coeff and mass matrix calculation functions from ITHACAutilities files to ITHACAcoeffsMass files.
1 parent 5733e82 commit d788f82

File tree

24 files changed

+556
-722
lines changed

24 files changed

+556
-722
lines changed

‎code_formatter.sh

+4
Original file line numberDiff line numberDiff line change
@@ -38,6 +38,10 @@ done
3838
echo $tutorial_files
3939
echo $src_files
4040

41+
for f in $(find . -name "*.H" -o -name "*.C"); do
42+
sed -i -e '$a\' $f
43+
done
44+
4145
# Here the important part: astyle formats the src files.
4246
astyle --style=bsd\
4347
--indent=spaces=4\

‎replace.sh

+6
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,6 @@
1+
#!/bin/bash
2+
3+
## template of a loop with sed to replace string within the library
4+
# for i in $(grep -Ril "stringToBeReplaced" folder/); do
5+
# sed -i 's/stringToBeReplaced/newString/g' $i
6+
# done

‎src/ITHACA_CORE/Containers/Modes.C

+3-3
Original file line numberDiff line numberDiff line change
@@ -96,7 +96,7 @@ Eigen::MatrixXd Modes<T>::project(GeometricField<T, fvPatchField, volMesh>&
9696
field, int numberOfModes)
9797
{
9898
Eigen::MatrixXd fieldEig = Foam2Eigen::field2Eigen(field);
99-
auto vol = ITHACAutilities::get_mass_matrix_FV(field);
99+
auto vol = ITHACAutilities::getMassMatrixFV(field);
100100
Eigen::MatrixXd projField;
101101

102102
if (EigenModes.size() == 0)
@@ -126,7 +126,7 @@ Eigen::MatrixXd Modes<T>::project(
126126
int numberOfModes)
127127
{
128128
Eigen::MatrixXd fieldEig = Foam2Eigen::PtrList2Eigen(fields);
129-
auto vol = ITHACAutilities::get_mass_matrix_FV(fields[0]);
129+
auto vol = ITHACAutilities::getMassMatrixFV(fields[0]);
130130
Eigen::MatrixXd projField;
131131

132132
if (EigenModes.size() == 0)
@@ -318,7 +318,7 @@ PtrList<GeometricField<T, fvPatchField, volMesh>>
318318

319319
if (innerProduct == "L2")
320320
{
321-
M_vol = ITHACAutilities::get_mass_matrix_FV(snapshots[i]);
321+
M_vol = ITHACAutilities::getMassMatrixFV(snapshots[i]);
322322
}
323323
else if (innerProduct == "Frobenius")
324324
{

‎src/ITHACA_CORE/ITHACAPOD/ITHACAPOD.C

+4-4
Original file line numberDiff line numberDiff line change
@@ -119,7 +119,7 @@ void ITHACAPOD::getModes(
119119
Eigen::MatrixXd SnapMatrix = Foam2Eigen::PtrList2Eigen(snapshots);
120120
List<Eigen::MatrixXd> SnapMatrixBC = Foam2Eigen::PtrList2EigenBC(snapshots);
121121
int NBC = snapshots[0].boundaryField().size();
122-
auto VM = ITHACAutilities::get_mass_matrix_FV(snapshots[0]);
122+
auto VM = ITHACAutilities::getMassMatrixFV(snapshots[0]);
123123
Eigen::MatrixXd _corMatrix = SnapMatrix.transpose() * VM.asDiagonal() *
124124
SnapMatrix;
125125

@@ -265,7 +265,7 @@ void ITHACAPOD::getWeightedModes(
265265
Eigen::MatrixXd SnapMatrix = Foam2Eigen::PtrList2Eigen(snapshots);
266266
List<Eigen::MatrixXd> SnapMatrixBC = Foam2Eigen::PtrList2EigenBC(snapshots);
267267
int NBC = snapshots[0].boundaryField().size();
268-
auto VM = ITHACAutilities::get_mass_matrix_FV(snapshots[0]);
268+
auto VM = ITHACAutilities::getMassMatrixFV(snapshots[0]);
269269
Eigen::MatrixXd _corMatrix = SnapMatrix.transpose() * VM.asDiagonal() *
270270
SnapMatrix;
271271
Eigen::VectorXd eigenValueseig;
@@ -409,7 +409,7 @@ void ITHACAPOD::getModesSVD(
409409
Info << "####### Performing POD using Singular Value Decomposition for " <<
410410
snapshots[0].name() << " #######" << endl;
411411
Eigen::MatrixXd SnapMatrix = Foam2Eigen::PtrList2Eigen(snapshots);
412-
Eigen::VectorXd V = ITHACAutilities::get_mass_matrix_FV(snapshots[0]);
412+
Eigen::VectorXd V = ITHACAutilities::getMassMatrixFV(snapshots[0]);
413413
Eigen::VectorXd V3dSqrt = V.array().sqrt();
414414
Eigen::VectorXd V3dInv = V3dSqrt.array().cwiseInverse();
415415
auto VMsqr = V3dSqrt.asDiagonal();
@@ -1588,7 +1588,7 @@ void ITHACAPOD::getModes(PtrList<Field_type>& snapshots,
15881588
Eigen::MatrixXd SnapMatrix = Foam2Eigen::PtrList2Eigen(fields2);
15891589
List<Eigen::MatrixXd> SnapMatrixBC = Foam2Eigen::PtrList2EigenBC(fields2);
15901590
int NBC = fields2[0].boundaryField().size();
1591-
auto VM = ITHACAutilities::get_mass_matrix_FV(fields2[0]);
1591+
auto VM = ITHACAutilities::getMassMatrixFV(fields2[0]);
15921592
Eigen::MatrixXd _corMatrix = SnapMatrix.transpose() * VM.asDiagonal() *
15931593
SnapMatrix;
15941594

‎src/ITHACA_CORE/ITHACAutilities/ITHACAcoeffsMass.C

+240
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,240 @@
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+
License
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+
27+
namespace ITHACAutilities
28+
{
29+
30+
template<typename TypeField>
31+
PtrList<TypeField> reconstructFromCoeff(
32+
PtrList<TypeField>& modes, Eigen::MatrixXd& coeff_matrix, label Nmodes)
33+
{
34+
PtrList<TypeField> rec_field;
35+
rec_field.resize(0);
36+
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+
else
45+
{
46+
rec_field[k] += modes[i] * coeff_matrix(i, k);
47+
}
48+
}
49+
50+
return rec_field;
51+
}
52+
53+
template PtrList<volScalarField> reconstructFromCoeff(
54+
PtrList<volScalarField>& modes, Eigen::MatrixXd& coeff_matrix, label Nmodes);
55+
template PtrList<volVectorField> reconstructFromCoeff(
56+
PtrList<volVectorField>& modes, Eigen::MatrixXd& coeff_matrix, label Nmodes);
57+
58+
template<typename T>
59+
Eigen::MatrixXd getMassMatrix(
60+
PtrList<GeometricField<T, fvPatchField, volMesh>>& modes, int Nmodes,
61+
bool consider_volumes)
62+
{
63+
label Msize;
64+
65+
if (Nmodes == 0)
66+
{
67+
Msize = modes.size();
68+
}
69+
else
70+
{
71+
Msize = Nmodes;
72+
}
73+
74+
M_Assert(modes.size() >= Msize,
75+
"The Number of requested modes is larger then the available quantity.");
76+
Eigen::MatrixXd F = Foam2Eigen::PtrList2Eigen(modes);
77+
Eigen::MatrixXd M;
78+
79+
if (consider_volumes)
80+
{
81+
Eigen::VectorXd V = getMassMatrixFV(modes[0]);
82+
M = F.transpose().topRows(Msize) * V.asDiagonal() * F.leftCols(Msize);
83+
}
84+
else
85+
{
86+
M = F.transpose().topRows(Msize) * F.leftCols(Msize);
87+
}
88+
89+
return M;
90+
}
91+
92+
template Eigen::MatrixXd getMassMatrix(
93+
PtrList<GeometricField<scalar, fvPatchField, volMesh>>& modes, int Nmodes,
94+
bool consider_volumes);
95+
96+
template Eigen::MatrixXd getMassMatrix(
97+
PtrList<GeometricField<vector, fvPatchField, volMesh>>& modes, int Nmodes,
98+
bool consider_volumes);
99+
100+
template<class TypeField>
101+
Eigen::VectorXd getMassMatrixFV(
102+
GeometricField<TypeField, fvPatchField, volMesh>& snapshot)
103+
{
104+
Eigen::MatrixXd snapEigen = Foam2Eigen::field2Eigen(snapshot);
105+
int dim = std::nearbyint(snapEigen.rows() / (snapshot.mesh().V()).size());
106+
Eigen::VectorXd volumes = Foam2Eigen::field2Eigen(snapshot.mesh().V());
107+
Eigen::VectorXd vol3 = volumes.replicate(dim, 1);
108+
return vol3;
109+
}
110+
111+
template Eigen::VectorXd getMassMatrixFV(
112+
GeometricField<scalar, fvPatchField, volMesh>& snapshot);
113+
template Eigen::VectorXd getMassMatrixFV(
114+
GeometricField<vector, fvPatchField, volMesh>& snapshot);
115+
116+
template<typename T>
117+
Eigen::VectorXd getCoeffs(GeometricField<T, fvPatchField, volMesh>&
118+
snapshot,
119+
PtrList<GeometricField<T, fvPatchField, volMesh>>& modes, int Nmodes,
120+
bool consider_volumes)
121+
{
122+
label Msize;
123+
124+
if (Nmodes == 0)
125+
{
126+
Msize = modes.size();
127+
}
128+
else
129+
{
130+
Msize = Nmodes;
131+
}
132+
133+
M_Assert(modes.size() >= Msize,
134+
"The Number of requested modes is larger then the available quantity.");
135+
Eigen::MatrixXd F = Foam2Eigen::PtrList2Eigen(modes);
136+
Eigen::MatrixXd M_matrix = getMassMatrix(modes, Nmodes, consider_volumes);
137+
Eigen::MatrixXd snapEigen = Foam2Eigen::field2Eigen(snapshot);
138+
Eigen::VectorXd a(Msize);
139+
Eigen::VectorXd b(Msize);
140+
141+
if (consider_volumes)
142+
{
143+
Eigen::VectorXd V = getMassMatrixFV(modes[0]);
144+
b = F.transpose().topRows(Msize) * V.asDiagonal() * snapEigen;
145+
}
146+
else
147+
{
148+
b = F.transpose().topRows(Msize) * snapEigen;
149+
}
150+
151+
a = M_matrix.fullPivLu().solve(b);
152+
return a;
153+
}
154+
155+
template Eigen::VectorXd getCoeffs(
156+
GeometricField<scalar, fvPatchField, volMesh>&
157+
snapshot, PtrList<GeometricField<scalar, fvPatchField, volMesh>>& modes,
158+
int Nmodes,
159+
bool consider_volumes);
160+
161+
template Eigen::VectorXd getCoeffs(
162+
GeometricField<vector, fvPatchField, volMesh>&
163+
snapshot, PtrList<GeometricField<vector, fvPatchField, volMesh>>& modes,
164+
int Nmodes,
165+
bool consider_volumes);
166+
167+
template<typename T>
168+
Eigen::MatrixXd getCoeffs(PtrList<GeometricField<T, fvPatchField, volMesh>>&
169+
snapshots,
170+
PtrList<GeometricField<T, fvPatchField, volMesh>>& modes, int Nmodes,
171+
bool consider_volumes)
172+
{
173+
label Msize;
174+
175+
if (Nmodes == 0)
176+
{
177+
Msize = modes.size();
178+
}
179+
else
180+
{
181+
Msize = Nmodes;
182+
}
183+
184+
M_Assert(modes.size() >= Msize,
185+
"The Number of requested modes is larger then the available quantity.");
186+
Eigen::MatrixXd coeff(Msize, snapshots.size());
187+
188+
for (auto i = 0; i < snapshots.size(); i++)
189+
{
190+
coeff.col(i) = getCoeffs(snapshots[i], modes, Nmodes, consider_volumes);
191+
}
192+
193+
return coeff;
194+
}
195+
196+
template Eigen::MatrixXd getCoeffs(
197+
PtrList<GeometricField<scalar, fvPatchField, volMesh>>&
198+
snapshot, PtrList<GeometricField<scalar, fvPatchField, volMesh>>& modes,
199+
int Nmodes,
200+
bool consider_volumes);
201+
202+
template Eigen::MatrixXd getCoeffs(
203+
PtrList<GeometricField<vector, fvPatchField, volMesh>>&
204+
snapshot, PtrList<GeometricField<vector, fvPatchField, volMesh>>& modes,
205+
int Nmodes,
206+
bool consider_volumes);
207+
208+
Eigen::MatrixXd parTimeCombMat(List<Eigen::VectorXd>
209+
acquiredSnapshotsTimes,
210+
Eigen::MatrixXd parameters)
211+
{
212+
int parsNum = parameters.cols();
213+
int parsSamplesNum = parameters.rows();
214+
M_Assert(parsSamplesNum == acquiredSnapshotsTimes.size(),
215+
"The list of time instants does not have the same number of vectors as the number of parameters samples");
216+
Eigen::MatrixXd comb;
217+
int totalSnapshotsNum = 0;
218+
219+
for (label k = 0; k < acquiredSnapshotsTimes.size(); k++)
220+
{
221+
totalSnapshotsNum += acquiredSnapshotsTimes[k].size();
222+
}
223+
224+
comb.resize(totalSnapshotsNum, parsNum + 1);
225+
label i = 0;
226+
227+
for (label j = 0; j < acquiredSnapshotsTimes.size(); j++)
228+
{
229+
for (label k = 0; k < acquiredSnapshotsTimes[j].size(); k++)
230+
{
231+
comb(i, parsNum) = (acquiredSnapshotsTimes[j])(k, 0);
232+
comb.block(i, 0, 1, parsNum) = parameters.row(j);
233+
i = i + 1;
234+
}
235+
}
236+
237+
return comb;
238+
}
239+
240+
}

0 commit comments

Comments
 (0)