Skip to content

Commit a2c7e82

Browse files
RomainTiphaigneRomainTiphaigne
committed
Fixing hyperreduction (source code and tutorials) with the new Foam2Eigen convention
1 parent 923bbae commit a2c7e82

10 files changed

Lines changed: 111 additions & 79 deletions

File tree

‎src/ITHACA_CORE/EigenFunctions/EigenFunctions.C‎

Lines changed: 24 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -142,4 +142,28 @@ Eigen::VectorXd repeatElements(Eigen::VectorXd input, int n)
142142
return output;
143143
}
144144
}
145+
146+
Eigen::VectorXd reorderVectorFromDim(Eigen::VectorXd input, int dim)
147+
{
148+
if ((input.size() % dim) != 0)
149+
{
150+
Info << "EigenFunctions::reorderVectorFromDim has incompatible vector size " <<
151+
"and reordering dimension. Aborting" << endl;
152+
abort();
153+
}
154+
else
155+
{
156+
Eigen::VectorXd output(input.size());
157+
int n = int(input.size() / dim);
158+
for (int i = 0; i < n; i++)
159+
{
160+
for (int d = 0; d < dim; d++)
161+
{
162+
output(i * dim + d) = input(i + d * dim);
163+
}
164+
}
165+
return output;
166+
}
167+
}
168+
145169
}

‎src/ITHACA_CORE/EigenFunctions/EigenFunctions.H‎

Lines changed: 12 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -260,6 +260,18 @@ Eigen::Matrix<T, Eigen::Dynamic, Eigen::Dynamic> vectorTensorProduct(
260260
/// @param n int
261261
Eigen::VectorXd repeatElements(Eigen::VectorXd input, int n);
262262

263+
//------------------------------------------------------------------------------
264+
/// @brief Changes the order of the vector's elements given a dimension
265+
///
266+
/// @details Example: input = (a_1, b_1, a_2, b_2, a_3, b_3); dim = 3
267+
/// output = (a_1, a_2, a_3, b_1, b_2, b_3)
268+
/// Typically switch from a component then coordinate order to the opposite
269+
/// as did the change in the Foam2Eigen mapping
270+
///
271+
/// @param input Eigen::VectorXd
272+
/// @param dim int
273+
Eigen::VectorXd reorderVectorFromDim(Eigen::VectorXd input, int dim);
274+
263275
};
264276

265277
template <typename T>

‎src/ITHACA_CORE/ITHACAutilities/ITHACAutilities.C‎

Lines changed: 0 additions & 19 deletions
Original file line numberDiff line numberDiff line change
@@ -232,24 +232,5 @@ Eigen::MatrixXd invertMatrix(Eigen::MatrixXd& matrixToInvert,
232232
}
233233
}
234234

235-
Eigen::SparseMatrix<double> sparseBlockDiagonalRepeat(Eigen::SparseMatrix<double>& A, int n)
236-
{
237-
std::vector<Eigen::Triplet<double>> tripletList;
238-
tripletList.reserve(n*A.nonZeros());
239-
240-
for (int i = 0; i < A.outerSize(); i++)
241-
{
242-
for (Eigen::SparseMatrix<double>::InnerIterator j(A,i); j; ++j)
243-
{
244-
for (int k = 0; k < n; k++)
245-
{
246-
tripletList.push_back(Eigen::Triplet<double>(j.row() + k*A.rows(), j.col() + k*A.cols(), j.value()));
247-
}
248-
}
249-
}
250-
Eigen::SparseMatrix<double> out(n*A.rows(), n*A.cols());
251-
out.setFromTriplets(tripletList.begin(), tripletList.end());
252-
return out;
253-
}
254235

255236
}

‎src/ITHACA_CORE/ITHACAutilities/ITHACAutilities.H‎

Lines changed: 0 additions & 6 deletions
Original file line numberDiff line numberDiff line change
@@ -148,12 +148,6 @@ Eigen::MatrixXd pinv_eigen_based(Eigen::MatrixXd& origin, const float er = 0);
148148
Eigen::MatrixXd invertMatrix(Eigen::MatrixXd& matrixToInvert,
149149
const word inversionMethod);
150150

151-
//------------------------------------------------------------------------------
152-
/// @brief Builds a block diagonal sparse matrix from n repetitions of the sparse matrix A
153-
///
154-
/// @param A Eigen::SparseMatrix
155-
/// @param n int
156-
Eigen::SparseMatrix<double> sparseBlockDiagonalRepeat(Eigen::SparseMatrix<double>& A, int n);
157151

158152
};
159153

‎src/ITHACA_FOMPROBLEMS/UnsteadyNSTurb/TurbDiffusionHyperreduction.C‎

Lines changed: 4 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -544,15 +544,15 @@ void TurbDiffusionHyperreduction<volF,T,S>::computeTurbDiffusionHyperreduction()
544544
{
545545
Eigen::VectorXd vectorWeights = ITHACAutilities::getMassMatrixFV(f_spatialModesU[0]).array().sqrt();
546546
Eigen::MatrixXd spatialModesU = vectorWeights.asDiagonal() * Foam2Eigen::PtrList2Eigen(f_spatialModesU);
547-
Eigen::SparseMatrix<double> P = ITHACAutilities::sparseBlockDiagonalRepeat(ithacaHyperreduction->P, 3);
547+
Eigen::SparseMatrix<double> P = ithacaHyperreduction->maskToOtherDim(3);
548548
Eigen::MatrixXd K_DEIM = spatialModesU.transpose() * P;
549549
for (label m = 0; m < ithacaHyperreduction->n_nodes; m++)
550550
{
551551
for (label c = 0; c < f_spatialModesU.size(); c++)
552552
{
553-
K_DEIM(c,m) *= quadWeights[c][m];
554-
K_DEIM(c,ithacaHyperreduction->n_nodes + m) *= quadWeights[c][m];
555-
K_DEIM(c,2*ithacaHyperreduction->n_nodes + m) *= quadWeights[c][m];
553+
K_DEIM(c, 3*m) *= quadWeights[c][m];
554+
K_DEIM(c, 3*m + 1) *= quadWeights[c][m];
555+
K_DEIM(c, 3*m + 2) *= quadWeights[c][m];
556556
}
557557
}
558558
cnpy::save(K_DEIM, folder_HR + "/K_DEIM.npy");

‎src/ITHACA_HR/hyperReduction.H‎

Lines changed: 7 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -526,6 +526,13 @@ class HyperReduction
526526
void getSnapMatrix(Eigen::MatrixXd& snapMatrix, Eigen::VectorXd& fieldWeights,
527527
List<Eigen::MatrixXd>& snapMatrixBoundary,
528528
List<Eigen::VectorXd>& fieldWeightsBoundary);
529+
530+
/// @brief Compute the sparse mask matrix P with a vectorial dimension different
531+
/// from the snapshots (enlarged or reduced)
532+
///
533+
/// @param[in] newDim int, the vectorial dim for the output
534+
///
535+
Eigen::SparseMatrix<double> maskToOtherDim(int newDim);
529536
};
530537

531538
#endif

‎src/ITHACA_HR/hyperReduction.templates.H‎

Lines changed: 41 additions & 16 deletions
Original file line numberDiff line numberDiff line change
@@ -554,8 +554,8 @@ void HyperReduction<SnapshotsLists...>::offlineGappyDEIM(
554554
adding.rowwise() += norma.transpose();
555555
Eigen::VectorXd denom = (adding.array().pow(1. /
556556
V.cols()).matrix()).rowwise().prod();
557-
max = ((mp_not_mask.head(n_cells).asDiagonal() * num).array() /
558-
denom.array()).maxCoeff(& ind_max, & c1);
557+
max = ((mp_not_mask(Eigen::seq(0, Eigen::indexing::last, vectorial_dim)).asDiagonal() * num)
558+
.array() / denom.array()).maxCoeff(& ind_max, & c1);
559559
}
560560
else
561561
{
@@ -732,8 +732,8 @@ void HyperReduction<SnapshotsLists...>::offlineECP(List<Eigen::MatrixXd>& listSn
732732

733733
for (unsigned int ith_field = 0; ith_field < vectorial_dim; ith_field++)
734734
{
735-
Eigen::MatrixXd block = listSnapshotsModes[ith_subspace].block(ith_field* n_cells, 0, n_cells,
736-
n_modesCurrent).transpose();
735+
Eigen::MatrixXd block = listSnapshotsModes[ith_subspace](Eigen::seq(ith_field,
736+
Eigen::indexing::last, vectorial_dim), Eigen::indexing::all).transpose();
737737

738738
if (predictVolume)
739739
{
@@ -815,8 +815,8 @@ void HyperReduction<SnapshotsLists...>::offlineECP(List<Eigen::MatrixXd>& listSn
815815

816816
while (nodePoints->size() < n_nodes)
817817
{
818-
max = ((b.transpose() * A * mp_not_mask.head(
819-
n_cells).asDiagonal()).colwise().sum()).maxCoeff(& r1, & ind_max);
818+
max = ((b.transpose() * A * mp_not_mask(Eigen::seq(0, Eigen::indexing::last, vectorial_dim))
819+
.asDiagonal()).colwise().sum()).maxCoeff(& r1, & ind_max);
820820
updateNodes(P, ind_max, mp_not_mask);
821821
n_modesPrevious = 0;
822822
for (int ith_subspace = 0; ith_subspace < n_subspaces; ith_subspace++)
@@ -873,7 +873,7 @@ void HyperReduction<SnapshotsLists...>::offlineECP(List<Eigen::MatrixXd>& listSn
873873
if (n_subspaces == 1)
874874
{
875875
basisMatrix = listSnapshotsModes[0];
876-
quadratureWeights = quadratureWeightsSubspaces[0];
876+
quadratureWeights = EigenFunctions::reorderVectorFromDim(quadratureWeightsSubspaces[0], vectorial_dim);
877877
evaluateWPU(P, basisMatrix, normalizingWeights, quadratureWeights);
878878
cnpy::save(basisMatrix, folderMethod + "/basisMatrix.npy");
879879
cnpy::save(quadratureWeights, folderMethod + "/quadratureWeights.npy");
@@ -884,6 +884,8 @@ void HyperReduction<SnapshotsLists...>::offlineECP(List<Eigen::MatrixXd>& listSn
884884
Eigen::MatrixXd quadratureWeightsMerged(n_nodes, n_subspaces);
885885
for (int ith_subspace = 0; ith_subspace < n_subspaces; ith_subspace++)
886886
{
887+
quadratureWeightsSubspaces[ith_subspace] = EigenFunctions::reorderVectorFromDim(
888+
quadratureWeightsSubspaces[ith_subspace], vectorial_dim);
887889
cnpy::save(quadratureWeightsSubspaces[ith_subspace], folderMethod + "/quadratureWeights_subspace" +
888890
std::to_string(ith_subspace+1) + ".npy");
889891
quadratureWeightsMerged.col(ith_subspace) = quadratureWeightsSubspaces[ith_subspace];
@@ -942,9 +944,9 @@ void HyperReduction<SnapshotsLists...>::initSeeds(Eigen::VectorXd mp_not_mask,
942944

943945
for (unsigned int ith_field = 0; ith_field < vectorial_dim; ith_field++)
944946
{
945-
P.insert(index + ith_field * n_cells,
946-
vectorial_dim * trueSeeds + ith_field) = 1;
947-
mp_not_mask(index + ith_field * n_cells) = 0;
947+
P.insert(index * vectorial_dim + ith_field,
948+
trueSeeds * vectorial_dim + ith_field) = 1;
949+
mp_not_mask(index * vectorial_dim + ith_field) = 0;
948950
}
949951

950952
trueSeeds++;
@@ -982,14 +984,14 @@ void HyperReduction<SnapshotsLists...>::updateNodes(Eigen::SparseMatrix<double>
982984
{
983985
for (unsigned int ith_field = 0; ith_field < vectorial_dim; ith_field++)
984986
{
985-
P.insert(nodes[ith_node] + ith_field * n_cells,
986-
ith_node + ith_field * (step + 1)) = 1;
987+
P.insert(nodes[ith_node] * vectorial_dim + ith_field,
988+
ith_node * vectorial_dim + ith_field) = 1;
987989
}
988990
}
989991

990992
for (unsigned int ith_field = 0; ith_field < vectorial_dim; ith_field++)
991993
{
992-
mp_not_mask(ind + ith_field * n_cells) = 0;
994+
mp_not_mask(ind * vectorial_dim + ith_field) = 0;
993995
}
994996
}
995997

@@ -1021,8 +1023,8 @@ void HyperReduction<SnapshotsLists...>::updateNodesRemove(Eigen::SparseMatrix<do
10211023
{
10221024
for (unsigned int ith_field = 0; ith_field < vectorial_dim; ith_field++)
10231025
{
1024-
P.insert(nodes[ith_node] + ith_field * n_cells,
1025-
ith_node + ith_field*(step + 1)) = 1;
1026+
P.insert(nodes[ith_node] * vectorial_dim + ith_field,
1027+
ith_node * vectorial_dim + ith_field) = 1;
10261028
}
10271029
}
10281030
// Not updating mp_not_mask : we are not choosing again a node that we removed
@@ -1062,7 +1064,7 @@ void HyperReduction<SnapshotsLists...>::initReshapeMat(
10621064
{
10631065
for (unsigned int ith_field = 0; ith_field < vectorial_dim; ith_field++)
10641066
{
1065-
reshapeMat.insert(i, i + n_cells * ith_field) = 1;
1067+
reshapeMat.insert(i, i * vectorial_dim + ith_field) = 1;
10661068
}
10671069
}
10681070

@@ -1307,4 +1309,27 @@ List<label> HyperReduction<SnapshotsLists...>::global2local(
13071309
return localPoints;
13081310
}
13091311

1312+
template<typename... SnapshotsLists>
1313+
Eigen::SparseMatrix<double> HyperReduction<SnapshotsLists...>::maskToOtherDim(int newDim)
1314+
{
1315+
if (!(newDim > 0))
1316+
{
1317+
Info << "The integer for the HyperReduction::maskToOtherDim method must be positive. Aborting" << endl;
1318+
abort();
1319+
}
1320+
else
1321+
{
1322+
Eigen::SparseMatrix<double> output(n_cells * newDim, nodes.size() * newDim);
1323+
for (int ith_node = 0; ith_node < nodes.size(); ith_node++)
1324+
{
1325+
for (unsigned int ith_field = 0; ith_field < newDim; ith_field++)
1326+
{
1327+
output.insert(nodes[ith_node] * newDim + ith_field,
1328+
ith_node * newDim + ith_field) = 1;
1329+
}
1330+
}
1331+
return output;
1332+
}
1333+
}
1334+
13101335
#endif

‎tutorials/CFD/24HyperReduction/24HyperReduction.C‎

Lines changed: 14 additions & 19 deletions
Original file line numberDiff line numberDiff line change
@@ -114,14 +114,12 @@ class HyperReduction_vectorFunction : public
114114

115115
for (auto i = 0; i < nodesList.size(); i++)
116116
{
117-
ret(i) = std::exp(- 2 * std::pow(xPos[nodesList[i]] - mu(0) - 1,
118-
2) - 2 * std::pow(yPos[nodesList[i]] - mu(1) - 0.5, 2));
119-
ret(i + nodesList.size()) = std::exp(- 2 * std::pow(xPos[nodesList[i]] - mu(
120-
0) - 0.5,
121-
2) - 2 * std::pow(yPos[nodesList[i]] - mu(1) - 0.5, 2));
122-
ret(i + 2 * nodesList.size()) = std::exp(- 2 * std::pow(
123-
xPos[nodesList[i]] - mu(0) - 1.5,
124-
2) - 2 * std::pow(yPos[nodesList[i]] - mu(1) - 0., 2));
117+
ret(3 * i) = std::exp(- 2 * std::pow(xPos[nodesList[i]] - mu(0) - 1, 2)
118+
- 2 * std::pow(yPos[nodesList[i]] - mu(1) - 0.5, 2));
119+
ret(3 * i + 1) = std::exp(- 2 * std::pow(xPos[nodesList[i]] - mu(0) - 0.5, 2)
120+
- 2 * std::pow(yPos[nodesList[i]] - mu(1) - 0.5, 2));
121+
ret(3 * i + 2) = std::exp(- 2 * std::pow(xPos[nodesList[i]] - mu(0) - 1.5, 2)
122+
- 2 * std::pow(yPos[nodesList[i]] - mu(1) - 0., 2));
125123
}
126124

127125
return ret;
@@ -170,17 +168,14 @@ class HyperReduction_vectorScalarFunction : public
170168

171169
for (auto i = 0; i < nodesList.size(); i++)
172170
{
173-
ret(i) = std::exp(- 2 * std::pow(xPos[nodesList[i]] - mu(0) - 1,
174-
2) - 2 * std::pow(yPos[nodesList[i]] - mu(1) - 0.5, 2));
175-
ret(i + nodesList.size()) = std::exp(- 2 * std::pow(xPos[nodesList[i]] - mu(
176-
0) - 0.5,
177-
2) - 2 * std::pow(yPos[nodesList[i]] - mu(1) - 0.5, 2));
178-
ret(i + 2 * nodesList.size()) = std::exp(- 2 * std::pow(
179-
xPos[nodesList[i]] - mu(0) - 1.5,
180-
2) - 2 * std::pow(yPos[nodesList[i]] - mu(1) - 0., 2));
181-
ret(i + 3 * nodesList.size()) = std::exp(- 2 * std::pow(
182-
xPos[nodesList[i]] - mu(0) - 1,
183-
2) - 2 * std::pow(yPos[nodesList[i]] - mu(1) - 0.5, 2));
171+
ret(4 * i) = std::exp(- 2 * std::pow(xPos[nodesList[i]] - mu(0) - 1, 2)
172+
- 2 * std::pow(yPos[nodesList[i]] - mu(1) - 0.5, 2));
173+
ret(4 * i + 1) = std::exp(- 2 * std::pow(xPos[nodesList[i]] - mu(0) - 0.5, 2)
174+
- 2 * std::pow(yPos[nodesList[i]] - mu(1) - 0.5, 2));
175+
ret(4 * i + 2) = std::exp(- 2 * std::pow(xPos[nodesList[i]] - mu(0) - 1.5, 2)
176+
- 2 * std::pow(yPos[nodesList[i]] - mu(1) - 0., 2));
177+
ret(4 * i + 3) = std::exp(- 2 * std::pow(xPos[nodesList[i]] - mu(0) - 1, 2)
178+
- 2 * std::pow(yPos[nodesList[i]] - mu(1) - 0.5, 2));
184179
}
185180

186181
return ret;

‎tutorials/CFD/27SmagorinskyHyperreduction/27Offline.C‎

Lines changed: 5 additions & 11 deletions
Original file line numberDiff line numberDiff line change
@@ -139,16 +139,10 @@ void tutorial27_offline::project()
139139
K_proj_line_i = K_DEIM.row(i);
140140
}
141141

142-
for (label j_p = 0 ; j_p < d ; j_p++)
143-
{
144-
for (label j = 0 ; j < nMagicPoints ; j++)
145-
{
146-
K_proj(i, j + j_p * (nMagicPoints + 1)) = K_proj_line_i(j + j_p * nMagicPoints);
147-
}
148-
}
149-
K_proj(i, nMagicPoints) = ITHACAutilities::dot_product_L2(spatialModesU[i], meanSmag);
150-
K_proj(i, nMagicPoints + nMagicPoints + 1) = 0.0;
151-
K_proj(i, nMagicPoints + 2 * (nMagicPoints + 1)) = 0.0;
142+
K_proj.row(i).head(d*nMagicPoints) = K_proj_line_i;
143+
K_proj(i, d*nMagicPoints) = ITHACAutilities::dot_product_L2(spatialModesU[i], meanSmag);
144+
K_proj(i, d*nMagicPoints + 1) = 0.0;
145+
K_proj(i, d*nMagicPoints + 2) = 0.0;
152146
}
153147

154148
cnpy::save(K_proj, file_projK_DEIM);
@@ -172,7 +166,7 @@ void tutorial27_offline::project()
172166

173167
if (m_parameters->get_HRMethod() == "DEIM")
174168
{
175-
K_DEIM.array().colwise() /= ITHACAutilities::getMassMatrixFV(spatialModesU[0]).head(K_DEIM.rows()).array().sqrt();
169+
K_DEIM.array().colwise() /= Foam2Eigen::field2Eigen(m_parameters->get_volume()).array().sqrt();
176170
for (label i = 0; i < nModesU; i++)
177171
{
178172
for (label q = 0; q < nModesU+1; q++)

‎tutorials/CFD/27SmagorinskyHyperreduction/27Online.C‎

Lines changed: 4 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -117,11 +117,11 @@ Eigen::VectorXd tutorial27_online::computeApproxSmagMg(const Eigen::VectorXd& re
117117
{
118118
vector stressVector = weightAtMg(i) * (stressField)[localMagicPoints[i]];
119119

120-
output(i) = stressVector.x();
121-
output(i + nMagicPoints + 1) = stressVector.y();
122-
output(i + 2 * (nMagicPoints + 1)) = stressVector.z();
120+
output(i * d) = stressVector.x();
121+
output(i * d + 1) = stressVector.y();
122+
output(i * d + 2) = stressVector.z();
123123
}
124-
output(nMagicPoints) = 1.0;
124+
output(d * nMagicPoints) = 1.0;
125125
return output;
126126
}
127127

0 commit comments

Comments
 (0)