Loading...
Searching...
No Matches
hyperReduction.templates.H
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 GN3U 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#ifndef hyperReduction_templates_H
32#define hyperReduction_templates_H
33
34#include "hyperReduction.H"
35
36// TODO increase modularity to avoid heavy compiling costs
37
38// Template function constructor
39template <typename... SnapshotsLists>
40HyperReduction<SnapshotsLists...>::HyperReduction(label n_modes,
41 label n_nodes,
42 Eigen::VectorXi initialSeeds,
43 word problemName,
44 SnapshotsLists &&...snapshotsLists)
45 : vectorial_dim{compute_vectorial_dim(snapshotsLists...)},
46 n_modes{n_modes},
47 n_nodes{n_nodes},
48 initialSeeds{initialSeeds},
49 problemName{problemName},
50 snapshotsListTuple{std::forward_as_tuple(snapshotsLists...)}
51{
52 Info << "Init HyperReduction class with vectorial dim: " << vectorial_dim <<
53 endl;
54 para = ITHACAparameters::getInstance();
55 folderProblem = "ITHACAoutput/" + problemName;
56 mkDir(folderProblem);
57 // TODO check snapshotsLists is not empty
58 // TODO check that the first snapshotsList is not empty (includes above)
59 // TODO initialize vectorial_dim only from SnapshotsLists type
60 // vectorial_dim = std::apply(compute_vectorial_dim<SnapshotsLists&&...>, snapshotsListTuple);
61 // Get fields' names
62 std::apply([this](auto & ...snapList)
63 {
64 (..., stackNames(snapList));
65 }, snapshotsListTuple);
66 // Get fields' dimensions
67 std::apply([this](auto & ...snapList)
68 {
69 (..., stackDimensions(snapList));
70 }, snapshotsListTuple);
71 sumFieldsDim = 0;
72 std::apply([this](auto & ...snapList)
73 {
74 (..., sumDimensions(this->sumFieldsDim, snapList));
75 }, snapshotsListTuple);
76 Info << "Fields names: ";
77
78 for (unsigned int ith_field = 0; ith_field < fieldNames.size(); ith_field++)
79 {
80 Info << fieldNames[ith_field] << " (dim=" << fieldDims[ith_field] << "); ";
81 }
82
83 Info << endl;
84 n_snapshots = std::get<0>(snapshotsListTuple).size();
85 Info << "The number of snapshots is: " << n_snapshots << endl;
86 n_cells = std::get<0>(snapshotsListTuple)[0].size();
87 Info << "The number of cells is: " << n_cells << endl;
88 Info << "Initial seeds length: " << initialSeeds.rows() << endl;
89 // get boundaries info
90 n_boundary_patches = std::get<0>(snapshotsListTuple)[0].boundaryField().size();
91 n_boundary_cells_list.resize(n_boundary_patches);
92 n_boundary_cells = 0;
93
94 for (unsigned int ith_boundary_patch = 0;
95 ith_boundary_patch < n_boundary_patches; ith_boundary_patch++)
96 {
97 n_boundary_cells_list[ith_boundary_patch] = std::get<0>
98 (snapshotsListTuple)[0].boundaryField()[ith_boundary_patch].size();
99 n_boundary_cells = n_boundary_cells + n_boundary_cells_list[ith_boundary_patch];
100 }
101
102 Info << "Number of boundary cells: " << n_boundary_cells << endl;
103
104 volumes = ITHACAutilities::getMassMatrixFV(std::get<0>(snapshotsListTuple)[0]).head(n_cells);
105}
106
107
108// Template function constructor
109template <typename... SnapshotsLists>
110HyperReduction<SnapshotsLists...>::HyperReduction(label n_modes,
111 label n_nodes,
112 unsigned int vectorialDim,
113 label n_cells,
114 Eigen::VectorXi initialSeeds,
115 word problemName,
116 Eigen::VectorXd volumes)
117 : n_modes{n_modes},
118 n_nodes{n_nodes},
119 vectorial_dim{vectorialDim},
120 n_cells{n_cells},
121 initialSeeds{initialSeeds},
122 problemName{problemName},
123 volumes{volumes}
124{
125 Info << "Init HyperReduction class with vectorial dim: " << vectorial_dim <<
126 endl;
127 para = ITHACAparameters::getInstance();
128 folderProblem = problemName;
129 mkDir(folderProblem);
130 Info << "Initial seeds length: " << initialSeeds.rows() << endl;
131}
132
133template <typename... SnapshotsLists>
134void HyperReduction<SnapshotsLists...>::getModesSVD(
135 SnapshotsListTuple& snapshotsListTuple, Eigen::MatrixXd& modesSVD,
136 Eigen::VectorXd& fieldWeights, bool saveModesFlag)
137{
138 // TODO move inside ITHACAPOD
139 ITHACAparameters* para(ITHACAparameters::getInstance());
140 word folderSVD = "ITHACAoutput/" + problemName + "/ModesSVD/";
141
142 if (!ITHACAutilities::check_folder(folderSVD))
143 {
144 mkDir(folderSVD);
145 // Initialize stacked snapshotsMatrix and normalizing weights
146 Eigen::MatrixXd snapMatrix;
147 getSnapMatrix(snapMatrix, fieldWeights);
148 word dimensionReductionMethod =
149 para->ITHACAdict->lookupOrDefault<word>("DimensionReductionMethod", "SVD");
150 Eigen::MatrixXd eigenVectoreig;
151
152 if (dimensionReductionMethod == "RSVD")
153 {
154 Info << "####### Performing RSVD for " << problemName << " #######" << endl;
155 unsigned int r_modeSVD =
156 para->ITHACAdict->lookupOrDefault<unsigned int>("RSVDdim", snapMatrix.cols());
157 RedSVD::RedSVD svd(fieldWeights.array().cwiseInverse().matrix().asDiagonal() *
158 snapMatrix, r_modeSVD);
159 Info << "####### End of RSVD for " << problemName << " #######" <<
160 endl;
161 eigenValueseig = svd.singularValues().real();
162 eigenVectoreig = svd.matrixU().real();
163 }
164 else
165 {
166 Info << "####### Performing SVD for " << problemName << " #######" << endl;
167 Eigen::JacobiSVD<Eigen::MatrixXd> svd(
168 fieldWeights.array().cwiseInverse().matrix().asDiagonal() * snapMatrix,
169 Eigen::ComputeThinU | Eigen::ComputeThinV);
170 Info << "####### End of SVD for " << problemName << " #######" <<
171 endl;
172 eigenValueseig = svd.singularValues().real();
173 eigenVectoreig = svd.matrixU().real();
174 }
175
176 modesSVD = eigenVectoreig;
177 // TODO correct boundary conditions
178 cnpy::save(eigenValueseig,
179 "ITHACAoutput/" + problemName + "/ModesSVD/eigenvalues.npy");
180 eigenValueseig = eigenValueseig / eigenValueseig.sum();
181 Eigen::VectorXd cumEigenValues(eigenValueseig);
182
183 for (label j = 1; j < cumEigenValues.size(); ++j)
184 {
185 cumEigenValues(j) += cumEigenValues(j - 1);
186 }
187
188 Eigen::saveMarketVector(eigenValueseig,
189 "ITHACAoutput/" + problemName + "/ModesSVD/Eigenvalues_" + problemName,
190 para->precision,
191 para->outytpe);
192 Eigen::saveMarketVector(cumEigenValues,
193 "ITHACAoutput/" + problemName + "/ModesSVD/CumEigenvalues_" + problemName,
194 para->precision,
195 para->outytpe);
196 cnpy::save(modesSVD, "ITHACAoutput/" + problemName + "/ModesSVD/modes.npy");
197 cnpy::save(fieldWeights,
198 "ITHACAoutput/" + problemName + "/ModesSVD/normalizingWeights.npy");
199
200 if (saveModesFlag)
201 {
202 unsigned int rowIndex{0};
203
204 for (unsigned int modeIndex = 0; modeIndex < modesSVD.cols(); modeIndex++)
205 {
206 std::apply([this, & modesSVD, & rowIndex, & modeIndex,
207 & folderSVD](auto & ...snapList)
208 {
209 (..., saveModes(snapList, modesSVD, rowIndex, modeIndex, folderSVD));
210 }, snapshotsListTuple);
211 rowIndex = 0;
212 }
213 }
214 }
215 else
216 {
217 Info << "Reading the existing modes" << endl;
218 cnpy::load(modesSVD, folderSVD + "/modes.npy");
219 cnpy::load(fieldWeights, folderSVD + "/normalizingWeights.npy");
220 }
221}
222
223template <typename... SnapshotsLists>
224void HyperReduction<SnapshotsLists...>::getModesSVD(
225 SnapshotsListTuple& snapshotsListTuple, Eigen::MatrixXd& modesSVD,
226 Eigen::VectorXd& fieldWeights, Eigen::MatrixXd& modesSVDBoundary,
227 Eigen::VectorXd& fieldWeightsBoundary, bool saveModesFlag)
228{
229 // TODO move inside ITHACAPOD
230 ITHACAparameters* para(ITHACAparameters::getInstance());
231 word folderSVD = "ITHACAoutput/" + problemName + "/ModesSVD/";
232
233 if (!ITHACAutilities::check_folder(folderSVD))
234 {
235 mkDir(folderSVD);
236 // Initialize stacked snapshotsMatrix and normalizing weights
237 Eigen::MatrixXd snapMatrix;
238 List<Eigen::MatrixXd> snapMatrixBoundary;
239 List<Eigen::VectorXd> fieldWeightsBoundaryList;
240 getSnapMatrix(snapMatrix, fieldWeights, snapMatrixBoundary,
241 fieldWeightsBoundaryList);
242
243 for (int id = 0; id < n_boundary_patches; id++)
244 {
245 snapMatrix.conservativeResize(snapMatrix.rows() + snapMatrixBoundary[id].rows(),
246 snapMatrix.cols());
247 snapMatrix.bottomRows(snapMatrixBoundary[id].rows()) = snapMatrixBoundary[id];
248 fieldWeightsBoundary.conservativeResize(fieldWeightsBoundary.rows() +
249 snapMatrixBoundary[id].rows());
250 fieldWeightsBoundary.tail(fieldWeightsBoundaryList[id].rows()) =
251 fieldWeightsBoundaryList[id];
252 }
253
254 cnpy::save(fieldWeights,
255 "ITHACAoutput/" + problemName + "/ModesSVD/normalizingWeights.npy");
256 cnpy::save(fieldWeightsBoundary,
257 "ITHACAoutput/" + problemName + "/ModesSVD/normalizingWeightsBoundary.npy");
258 fieldWeights.conservativeResize(fieldWeights.rows() +
259 fieldWeightsBoundary.rows());
260 fieldWeights.tail(fieldWeightsBoundary.rows()) = fieldWeightsBoundary;
261 word dimensionReductionMethod =
262 para->ITHACAdict->lookupOrDefault<word>("DimensionReductionMethod", "SVD");
263 Eigen::MatrixXd eigenVectoreig;
264
265 if (dimensionReductionMethod == "RSVD")
266 {
267 Info << "####### Performing RSVD for " << problemName << " #######" << endl;
268 unsigned int r_modeSVD =
269 para->ITHACAdict->lookupOrDefault<unsigned int>("RSVDdim", snapMatrix.cols());
270 RedSVD::RedSVD svd(fieldWeights.array().cwiseInverse().matrix().asDiagonal() *
271 snapMatrix, r_modeSVD);
272 Info << "####### End of RSVD for " << problemName << " #######" << endl;
273 eigenValueseig = svd.singularValues().real();
274 eigenVectoreig = svd.matrixU().real();
275 }
276 else
277 {
278 Info << "####### Performing SVD for " << problemName << " #######" << endl;
279 Eigen::JacobiSVD<Eigen::MatrixXd> svd(
280 fieldWeights.array().cwiseInverse().matrix().asDiagonal() * snapMatrix,
281 Eigen::ComputeThinU | Eigen::ComputeThinV);
282 Info << "####### End of SVD for " << problemName << " #######" <<
283 endl;
284 eigenValueseig = svd.singularValues().real();
285 eigenVectoreig = svd.matrixU().real();
286 }
287
288 modesSVD = eigenVectoreig.bottomRows(vectorial_dim * n_cells);
289 modesSVDBoundary = eigenVectoreig.topRows(n_boundary_cells);
290 fieldWeights.conservativeResize(fieldWeights.rows() -
291 fieldWeightsBoundary.rows());
292 cnpy::save(eigenValueseig,
293 "ITHACAoutput/" + problemName + "/ModesSVD/eigenvalues.npy");
294 eigenValueseig = eigenValueseig / eigenValueseig.sum();
295 Eigen::VectorXd cumEigenValues(eigenValueseig);
296
297 for (label j = 1; j < cumEigenValues.size(); ++j)
298 {
299 cumEigenValues(j) += cumEigenValues(j - 1);
300 }
301
302 Eigen::saveMarketVector(eigenValueseig,
303 "ITHACAoutput/" + problemName + "/ModesSVD/Eigenvalues_" + problemName,
304 para->precision,
305 para->outytpe);
306 Eigen::saveMarketVector(cumEigenValues,
307 "ITHACAoutput/" + problemName + "/ModesSVD/CumEigenvalues_" + problemName,
308 para->precision,
309 para->outytpe);
310 cnpy::save(modesSVD, "ITHACAoutput/" + problemName + "/ModesSVD/modes.npy");
311 cnpy::save(modesSVDBoundary,
312 "ITHACAoutput/" + problemName + "/ModesSVD/modesBoundary.npy");
313
314 if (saveModesFlag)
315 {
316 unsigned int rowIndex{0};
317 unsigned int rowIndexBoundary{0};
318
319 for (unsigned int modeIndex = 0; modeIndex < modesSVD.cols(); modeIndex++)
320 {
321 std::apply([this, & modesSVD, & modesSVDBoundary, & rowIndex, & modeIndex,
322 & rowIndexBoundary, & folderSVD](auto & ...snapList)
323 {
324 (..., saveModes(snapList, modesSVD, modesSVDBoundary, rowIndex,
325 rowIndexBoundary, modeIndex, folderSVD));
326 }, snapshotsListTuple);
327 rowIndex = 0;
328 }
329 }
330 }
331 else
332 {
333 Info << "Reading the existing modes" << endl;
334 cnpy::load(modesSVD, folderSVD + "/modes.npy");
335 cnpy::load(modesSVDBoundary,
336 "ITHACAoutput/" + problemName + "/ModesSVD/modesBoundary.npy");
337 cnpy::load(fieldWeights, folderSVD + "/normalizingWeights.npy");
338 cnpy::load(fieldWeightsBoundary,
339 "ITHACAoutput/" + problemName + "/ModesSVD/normalizingWeightsBoundary.npy");
340 }
341}
342
343template <typename... SnapshotsLists>
344void HyperReduction<SnapshotsLists...>::getSnapMatrix(Eigen::MatrixXd&
345 snapMatrix, Eigen::VectorXd& fieldWeights)
346{
347 std::apply([this, & fieldWeights, & snapMatrix](auto & ...snapList)
348 {
349 (..., stackSnapshots(snapList, snapMatrix, fieldWeights));
350 }, snapshotsListTuple);
351}
352
353template <typename... SnapshotsLists>
354void HyperReduction<SnapshotsLists...>::getSnapMatrix(Eigen::MatrixXd&
355 snapMatrix, Eigen::VectorXd& fieldWeights,
356 List<Eigen::MatrixXd>& snapMatrixBoundary,
357 List<Eigen::VectorXd>& fieldWeightsBoundary)
358{
359 fieldWeightsBoundary.resize(n_boundary_patches);
360 snapMatrixBoundary.resize(n_boundary_patches);
361 std::apply([this, & fieldWeights, & snapMatrix, & fieldWeightsBoundary,
362 & snapMatrixBoundary](auto & ...snapList)
363 {
364 (..., stackSnapshots(snapList, snapMatrix, fieldWeights));
365 (..., stackSnapshotsBoundary(snapList, snapMatrixBoundary,
366 fieldWeightsBoundary));
367 },
368 snapshotsListTuple);
369}
370
371template <typename... SnapshotsLists>
372template <typename SnapshotsList>
373void HyperReduction<SnapshotsLists...>::saveModes(SnapshotsList sList,
374 Eigen::MatrixXd& snapshotsMatrix, unsigned int& rowIndex,
375 unsigned int& modeIndex, word folder)
376{
377 unsigned int field_dim = get_field_dim<typename SnapshotsList::value_type>();
378 auto fieldName = sList[0].name();
379 unsigned int fieldSize = field_dim * sList[0].size();
380 Eigen::VectorXd fieldBlock = snapshotsMatrix.block(rowIndex, modeIndex,
381 fieldSize, 1);
382 auto fieldOut = Foam2Eigen::Eigen2field(sList[0], fieldBlock, true);
383 rowIndex += fieldSize;
384 ITHACAstream::exportSolution(fieldOut, name(modeIndex), folder);
385}
386
387template <typename... SnapshotsLists>
388template <typename SnapshotsList>
389void HyperReduction<SnapshotsLists...>::saveModes(SnapshotsList sList,
390 Eigen::MatrixXd& snapshotsMatrix, Eigen::MatrixXd& snapshotsMatrixBoundary,
391 unsigned int& rowIndex, unsigned int& rowIndexBoundary,
392 unsigned int& modeIndex, word folder)
393{
394 unsigned int field_dim = get_field_dim<typename SnapshotsList::value_type>();
395 auto fieldName = sList[0].name();
396 unsigned int fieldSize = field_dim * sList[0].size();
397 Eigen::VectorXd fieldBlock = snapshotsMatrix.block(rowIndex, modeIndex,
398 fieldSize, 1);
399 List<Eigen::VectorXd> fieldBlockBoundary;
400 fieldBlockBoundary.resize(n_boundary_patches);
401
402 for (unsigned int id = 0; id < n_boundary_patches; id++)
403 {
404 unsigned int bfieldDim = field_dim * int(n_boundary_cells_list[id]);
405 fieldBlockBoundary[id] = snapshotsMatrixBoundary.block(rowIndexBoundary,
406 modeIndex, bfieldDim, 1);
407 rowIndexBoundary += bfieldDim;
408 }
409
410 auto fieldOut = Foam2Eigen::Eigen2field(sList[0], fieldBlock,
411 fieldBlockBoundary);
412 rowIndex += fieldSize;
413 ITHACAstream::exportSolution(fieldOut, name(modeIndex), folder);
414}
415
416template <typename... SnapshotsLists>
417void HyperReduction<SnapshotsLists...>::offlineGappyDEIM(
418 Eigen::MatrixXd& snapshotsModes, Eigen::VectorXd& weights,
419 word folderMethodName)
420{
421 folderMethod = folderMethodName;
422 Info << "FolderMethod : " << folderMethod << endl;
423 mkDir(folderMethod);
424 normalizingWeights = weights;
425 nodePoints = autoPtr<IOList<label >> (new IOList<label>(
426 IOobject("nodePoints", para->runTime.time().constant(), "../" + folderMethod,
427 para->mesh, IOobject::READ_IF_PRESENT, IOobject::NO_WRITE)));
428 auto greedyMetric = para->ITHACAdict->lookupOrDefault<word>("GreedyMetric",
429 "L2");
430 bool offlineStage = !nodePoints().headerOk();
431
432 if (offlineStage)
433 {
434 assert(n_modes > 0);
435 assert(n_nodes >= n_modes);
436 Eigen::VectorXd mp_not_mask = Eigen::VectorXd::Constant(n_cells * vectorial_dim,
437 1);
438 std::set<label> nodePointsSet;
439
440 // set initialSeeds
441 if (initialSeeds.rows() > 0)
442 {
443 initSeeds(mp_not_mask, nodePointsSet);
444 }
445
446 int na = n_nodes - initialSeeds.rows();
447
448 if (na > 0)
449 {
450 Eigen::SparseMatrix<double> reshapeMat;
451 initReshapeMat(reshapeMat);
452
453 if (greedyMetric == "L2")
454 {
455 Info << "####### Begin Greedy-L2 GappyDEIM #######" << endl;
456 Info << "####### Modes=" << n_modes
457 << ", nodePoints=" << n_nodes << " #######"
458 << endl;
459 int nb = 0;
460 int nit = std::min(n_modes, na);
461 int ncimin = std::floor(n_modes / nit);
462 int naimin = std::floor(na / n_modes);
463 Eigen::MatrixXd A;
464 Eigen::VectorXd b;
465 Eigen::VectorXd c;
466 Eigen::VectorXd r;
467 label ind_max, c1;
468 double max;
469
470 for (int i = 1; i <= nit; i++)
471 {
472 int nci = ncimin;
473
474 // add the remaining modes in the quotient n_modes / nite
475 if (i <= n_modes % nit)
476 {
477 nci = ncimin + 1;
478 }
479
480 int nai = naimin;
481
482 // add the remaining nodePoints in the quotient na / n_modes
483 if (i <= na % n_modes)
484 {
485 nai = naimin + 1;
486 }
487
488 Eigen::MatrixXd V;
489
490 // select basis
491 if (i == 1)
492 {
493 V = snapshotsModes.leftCols(nci);
494 }
495 else
496 {
497 for (int q = 1; q <= nci; q++)
498 {
499 A = P.transpose() * basisMatrix;
500 b = P.transpose() * snapshotsModes.col(nb + q - 1);
501 c = A.fullPivLu().solve(b);
502 r = snapshotsModes.col(nb + q - 1) - basisMatrix * c;
503 V.conservativeResize(snapshotsModes.rows(), q);
504 V.col(q - 1) = r;
505 }
506 }
507
508 for (int j = 1; j <= nai; j++)
509 {
510 if (P.cols() > 0)
511 {
512 max = (reshapeMat * (mp_not_mask.asDiagonal() * V)
513 .rowwise()
514 .lpNorm<2>().array().square().matrix())
515 .maxCoeff(& ind_max, & c1);
516 }
517 else
518 {
519 max = (reshapeMat * V.rowwise().lpNorm<2>().array().square().matrix()).maxCoeff(
520 & ind_max, & c1);
521 }
522
523 updateNodes(P, ind_max, mp_not_mask);
524 }
525
526 nb += nci;
527 basisMatrix = snapshotsModes.leftCols(nb);
528 }
529
530 Info << "####### End of greedy GappyDEIM #######\n";
531 }
532 else if (greedyMetric == "SOPT")
533 {
534 Info << "####### Begin SOPT GappyDEIM #######" << endl;
535 Info << "####### Modes=" << n_modes
536 << ", nodePoints=" << n_nodes << " #######"
537 << endl;
538 label ind_max, c1;
539 double max;
540 Eigen::MatrixXd V = snapshotsModes.leftCols(n_modes);
541
542 for (int ith_node = 0; ith_node < na; ith_node++)
543 {
544 if (P.cols() > 0)
545 {
546 // TODO make efficient
547 Eigen::MatrixXd tmp = P.transpose() * V;
548 Eigen::MatrixXd R = mp_not_mask.asDiagonal() * V;
549 Eigen::MatrixXd inv = (tmp.transpose() * tmp).fullPivLu().inverse();
550 Eigen::VectorXd num = 1 + (reshapeMat * ((R * inv).array() *
551 R.array()).rowwise().sum().matrix()).array();
552 Eigen::VectorXd norma = (tmp.array().square()).colwise().sum();
553 Eigen::MatrixXd adding = reshapeMat * R.array().square().matrix();
554 adding.rowwise() += norma.transpose();
555 Eigen::VectorXd denom = (adding.array().pow(1. /
556 V.cols()).matrix()).rowwise().prod();
557 max = ((mp_not_mask(Eigen::seq(0, Eigen::indexing::last, vectorial_dim)).asDiagonal() * num)
558 .array() / denom.array()).maxCoeff(& ind_max, & c1);
559 }
560 else
561 {
562 max = (reshapeMat * V.rowwise().lpNorm<2>().array().square().matrix()).maxCoeff(
563 & ind_max, & c1);
564 }
565
566 updateNodes(P, ind_max, mp_not_mask);
567 }
568
569 basisMatrix = V;
570 Info << "####### End of SOPT GappyDEIM #######\n";
571 }
572 else if (greedyMetric == "SOPTE")
573 {
574 Info << "####### Begin SOPT-exact GappyDEIM #######" << endl;
575 Info << "####### Modes=" << n_modes
576 << ", nodePoints=" << n_nodes << " #######"
577 << endl;
578 label ind_max, c1;
579 double max;
580 Eigen::MatrixXd A;
581 Eigen::VectorXd b;
582 Eigen::VectorXd c;
583 Eigen::VectorXd r;
584 Eigen::MatrixXd V = snapshotsModes.col(0);
585 max = (reshapeMat * (mp_not_mask.asDiagonal() * snapshotsModes.leftCols(
586 n_modes)).rowwise().lpNorm<2>().array().square().matrix()).maxCoeff(& ind_max,
587 & c1);
588 updateNodes(P, ind_max, mp_not_mask);
589
590 for (int ith_node = 1; ith_node < na; ith_node++)
591 {
592 if (ith_node < n_modes)
593 {
594 A = P.transpose() * V;
595 b = P.transpose() * snapshotsModes.col(ith_node);
596 c = A.fullPivLu().solve(b);
597 r = snapshotsModes.col(ith_node) - V * c;
598 V.conservativeResize(snapshotsModes.rows(), ith_node + 1);
599 V.col(ith_node) = snapshotsModes.col(ith_node);
600 max = (reshapeMat * r.rowwise().lpNorm<2>().array().square().matrix()).maxCoeff(
601 & ind_max, & c1);
602 }
603 else
604 {
605 // TODO optimize
606 Eigen::MatrixXd tmp = P.transpose() * snapshotsModes.leftCols(n_modes);
607 Info << "SOPT: " << s_optimality(tmp) << endl;
608 tmp.conservativeResize(tmp.rows() + vectorial_dim, tmp.cols());
609 Eigen::MatrixXd masked = mp_not_mask.asDiagonal() * snapshotsModes.leftCols(
610 n_modes);
611 Eigen::VectorXd results(n_cells);
612
613 for (int ith_cell = 0; ith_cell < n_cells; ith_cell++)
614 {
615 for (unsigned int ith_field = 0; ith_field < vectorial_dim; ith_field++)
616 {
617 tmp.row(tmp.rows() - vectorial_dim + ith_field) = masked.row(
618 ith_cell + ith_field * n_cells);
619 }
620
621 results(ith_cell) = s_optimality(tmp);
622 }
623
624 max = results.maxCoeff(& ind_max, & c1);
625 }
626
627 updateNodes(P, ind_max, mp_not_mask);
628 }
629
630 basisMatrix = snapshotsModes.leftCols(n_modes);
631 Info << "####### End of SOPTE GappyDEIM #######\n";
632 }
633 }
634 else
635 {
636 basisMatrix = snapshotsModes.leftCols(n_modes);
637 Info << "####### InitialSeeds equal to n_nodes #######\n";
638 }
639
640 evaluatePinv(P, basisMatrix, normalizingWeights);
641 renormalizedBasisMatrix = normalizingWeights.asDiagonal() * basisMatrix;
642 MatrixOnline = renormalizedBasisMatrix * pinvPU;
643 P.makeCompressed();
644 cnpy::save(basisMatrix, folderMethod + "/basisMatrix.npy");
645 cnpy::save(P, folderMethod + "/projectionMatrix.npz");
646 cnpy::save(normalizingWeights, folderMethod + "/normalizingWeights.npy");
647 cnpy::save(nodes, folderMethod + "/mp.npy");
648 cnpy::save(pinvPU, folderMethod + "/pinvPU.npy");
649 cnpy::save(MatrixOnline, folderMethod + "/MatrixOnline.npy");
650 nodePoints().write();
651 Info << "Projection Matrix shape: " << P.rows() << " " << P.cols() << endl;
652 Info << "Basis Matrix shape: " << basisMatrix.rows() << " " <<
653 basisMatrix.cols() << endl;
654 Info << "Pseudo Inverse Matrix shape: " << pinvPU.rows() << " " << pinvPU.cols()
655 << endl;
656 }
657 else
658 {
659 Info << "Read GappyDEIM\n";
660 cnpy::load(basisMatrix, folderMethod + "/basisMatrix.npy");
661 cnpy::load(MatrixOnline, folderMethod + "/MatrixOnline.npy");
662 cnpy::load(normalizingWeights, folderMethod + "/normalizingWeights.npy");
663 cnpy::load(P, folderMethod + "/projectionMatrix.npz");
664 cnpy::load(nodes, folderMethod + "/mp.npy");
665 cnpy::load(pinvPU, folderMethod + "/pinvPU.npy");
666 renormalizedBasisMatrix = normalizingWeights.asDiagonal() * basisMatrix;
667 Info << "Projection Matrix shape: " << P.rows() << " " << P.cols() << endl;
668 Info << "Basis Matrix shape: " << basisMatrix.rows() << " " <<
669 basisMatrix.cols() << endl;
670 Info << "Pseudo Inverse Matrix shape: " << pinvPU.rows() << " " << pinvPU.cols()
671 << endl;
672 }
673}
674
675template <typename... SnapshotsLists>
676void HyperReduction<SnapshotsLists...>::offlineECP(List<Eigen::MatrixXd>& listSnapshotsModes,
677 Eigen::VectorXd& weights, word folderMethodName)
678{
679 folderMethod = folderMethodName;
680 Info << "FolderMethod : " << folderMethod << endl;
681 mkDir(folderMethod);
682 // n_nodes = n_modes + 1;
683 normalizingWeights = weights;
684
685 label n_subspaces = listSnapshotsModes.size();
686 List<label> n_subspaceModes(0);
687 label n_totalModes = 0;
688 for (int ith_subspace = 0; ith_subspace < n_subspaces; ith_subspace++)
689 {
690 n_subspaceModes.append(listSnapshotsModes[ith_subspace].cols());
691 n_totalModes += listSnapshotsModes[ith_subspace].cols();
692 }
693 quadratureWeightsSubspaces.resize(n_subspaces);
694
695 Info << "####### Begin ECP #######" << endl;
696 Info << "####### Modes=" << n_subspaceModes
697 << ", nodePoints=" << n_nodes << " #######"
698 << endl;
699 nodePoints = autoPtr<IOList<label >> (new IOList<label>(
700 IOobject("nodePoints", para->runTime.time().constant(), "../" + folderMethod,
701 para->mesh, IOobject::READ_IF_PRESENT, IOobject::NO_WRITE)));
702 bool offlineStage = !nodePoints().headerOk();
703
704 if (offlineStage)
705 {
706
707 bool normalizedArgmaxForCubaturePointsSelection = 1;
708 bool cellVolumesInCubatureWeights = 0; // Gives an equivalent system, but changes the condition number.
709 bool predictVolume = 0; // Adding the prediction of the volume in the system
710 bool NNLS = 1; // Weights are all positive
711
712 assert(n_totalModes > 0);
713 // assert(n_nodes >= n_modes);
714
715 // matrices for quadratureWeights evaluation
716 Eigen::MatrixXd J;
717 List<Eigen::MatrixXd> Jwhole(n_subspaces);
718 List<Eigen::VectorXd> q(n_subspaces);
719 // matrices for greedy selection of the nodes
720 Eigen::MatrixXd A(vectorial_dim * n_totalModes, n_cells);
721 Eigen::VectorXd b = Eigen::VectorXd::Constant(vectorial_dim * n_totalModes, 1);
722 Eigen::VectorXd cellVolumes = volumes;
723 double volume = cellVolumes.array().sum();
724 label n_modesPrevious = 0;
725
726 for (int ith_subspace = 0; ith_subspace < n_subspaces; ith_subspace++)
727 {
728 label n_modesCurrent = n_subspaceModes[ith_subspace];
729 label n_modesCurrent_modif = n_modesCurrent + predictVolume;
730 q[ith_subspace] = Eigen::VectorXd::Zero(vectorial_dim * n_modesCurrent_modif);
731 Jwhole[ith_subspace] = Eigen::MatrixXd::Zero(vectorial_dim * n_modesCurrent_modif, n_cells);
732
733 for (unsigned int ith_field = 0; ith_field < vectorial_dim; ith_field++)
734 {
735 Eigen::MatrixXd block = listSnapshotsModes[ith_subspace](Eigen::seq(ith_field,
736 Eigen::indexing::last, vectorial_dim), Eigen::indexing::all).transpose();
737
738 if (predictVolume)
739 {
740 Eigen::MatrixXd volumeAwareMean = (block.rowwise().sum()/volume) * cellVolumes.transpose();
741 block -= volumeAwareMean;
742
743 q[ith_subspace](ith_field * n_modesCurrent_modif + n_modesCurrent) = volume;
744 Jwhole[ith_subspace].row(ith_field * n_modesCurrent_modif + n_modesCurrent) = cellVolumes;
745 }
746
747 q[ith_subspace].segment(ith_field * n_modesCurrent_modif, n_modesCurrent) = block.rowwise().sum();
748 Jwhole[ith_subspace].middleRows(ith_field * n_modesCurrent_modif, n_modesCurrent) = block;
749
750 if (cellVolumesInCubatureWeights)
751 {
752 Jwhole[ith_subspace].middleRows(ith_field * n_modesCurrent_modif, n_modesCurrent_modif)
753 .array().rowwise() /= cellVolumes.transpose().array().sqrt();
754 }
755
756 if (normalizedArgmaxForCubaturePointsSelection)
757 {
758 Eigen::VectorXd mean = block.rowwise().mean();
759 block.colwise() -= mean;
760 Eigen::VectorXd Anorm = block.colwise().lpNorm<2>();
761
762 if (Anorm.minCoeff() > 0)
763 {
764 block.array().rowwise() /= Anorm.transpose().array();
765 }
766
767 A.middleRows(n_modesPrevious + ith_field* n_modesCurrent, n_modesCurrent) = block;
768 }
769 else
770 {
771 A.middleRows(n_modesPrevious + ith_field * n_modesCurrent, n_modesCurrent) =
772 Jwhole[ith_subspace].middleRows(ith_field * n_modesCurrent_modif, n_modesCurrent);
773 }
774 }
775 n_modesPrevious += vectorial_dim * n_modesCurrent;
776 }
777
778 Eigen::VectorXd mp_not_mask = Eigen::VectorXd::Constant(n_cells * vectorial_dim,
779 1);
780 std::set<label> nodePointsSet;
781
782 // set initialSeeds
783 if (initialSeeds.rows() > 0)
784 {
785 initSeeds(mp_not_mask, nodePointsSet);
786 n_modesPrevious = 0;
787 for (int ith_subspace = 0; ith_subspace < n_subspaces; ith_subspace++)
788 {
789 label n_modesCurrent = n_subspaceModes[ith_subspace];
790 label n_modesCurrent_modif = n_modesCurrent + predictVolume;
791 quadratureWeightsSubspaces[ith_subspace].resize(vectorial_dim * nodePoints->size());
792
793 for (unsigned int ith_field = 0; ith_field < vectorial_dim; ith_field++)
794 {
795 Eigen::VectorXd b_tmp(n_modesCurrent);
796 Eigen::MatrixXd J_partial = Jwhole[ith_subspace].middleRows(ith_field * n_modesCurrent_modif, n_modesCurrent_modif);
797 Eigen::VectorXd q_partial = q[ith_subspace].segment(ith_field * n_modesCurrent_modif, n_modesCurrent_modif);
798 computeLS(J, J_partial, b_tmp, q_partial, NNLS);
799 b.segment(n_modesPrevious,n_modesCurrent) = b_tmp;
800 quadratureWeightsSubspaces[ith_subspace].segment(ith_field * nodePoints->size(),
801 nodePoints->size()) = quadratureWeights;
802 n_modesPrevious += n_modesCurrent;
803 }
804 }
805 }
806
807 int na = n_nodes - initialSeeds.rows();
808 Eigen::SparseMatrix<double> reshapeMat;
809 initReshapeMat(reshapeMat);
810
811 if (na > 0)
812 {
813 label ind_max, r1;
814 double max;
815
816 while (nodePoints->size() < n_nodes)
817 {
818 max = ((b.transpose() * A * mp_not_mask(Eigen::seq(0, Eigen::indexing::last, vectorial_dim))
819 .asDiagonal()).colwise().sum()).maxCoeff(& r1, & ind_max);
820 updateNodes(P, ind_max, mp_not_mask);
821 n_modesPrevious = 0;
822 for (int ith_subspace = 0; ith_subspace < n_subspaces; ith_subspace++)
823 {
824 label n_modesCurrent = n_subspaceModes[ith_subspace];
825 label n_modesCurrent_modif = n_modesCurrent + predictVolume;
826 quadratureWeightsSubspaces[ith_subspace].resize(vectorial_dim * nodePoints->size());
827
828 for (unsigned int ith_field = 0; ith_field < vectorial_dim; ith_field++)
829 {
830 Eigen::VectorXd b_tmp(n_modesCurrent);
831 Eigen::MatrixXd J_partial = Jwhole[ith_subspace].middleRows(ith_field * n_modesCurrent_modif, n_modesCurrent_modif);
832 Eigen::VectorXd q_partial = q[ith_subspace].segment(ith_field * n_modesCurrent_modif, n_modesCurrent_modif);
833
834 computeLS(J, J_partial, b_tmp, q_partial, NNLS);
835 b.segment(n_modesPrevious,n_modesCurrent) = b_tmp;
836 quadratureWeightsSubspaces[ith_subspace].segment(ith_field * nodePoints->size(),
837 nodePoints->size()) = quadratureWeights;
838 n_modesPrevious += n_modesCurrent;
839 }
840 }
841 if (NNLS) // Removes nodes with weight 0 from set of nodes
842 {
843 Eigen::VectorXd maxQuadWeights = Eigen::VectorXd::Zero(nodePoints->size());
844 for (int ith_subspace = 0; ith_subspace < n_subspaces; ith_subspace++)
845 {
846 for (unsigned int ith_field = 0; ith_field < vectorial_dim; ith_field++)
847 {
848 maxQuadWeights = maxQuadWeights.cwiseMax(quadratureWeightsSubspaces[ith_subspace]
849 .segment(ith_field * nodePoints->size(),nodePoints->size()));
850 }
851 }
852 updateNodesRemove(P, maxQuadWeights);
853 }
854 }
855 }
856
857 if (cellVolumesInCubatureWeights)
858 {
859 for (label i=0; i<nodePoints->size(); i++)
860 {
861 for (int ith_subspace = 0; ith_subspace < n_subspaces; ith_subspace++)
862 {
863 for (unsigned int ith_field = 0; ith_field < vectorial_dim; ith_field++)
864 {
865 quadratureWeightsSubspaces[ith_subspace](ith_field * nodePoints->size() + i) /=
866 std::sqrt(volumes(nodePoints()[i]));
867 }
868 }
869 }
870 }
871
872 Info << "####### End ECP #######" << endl;
873 if (n_subspaces == 1)
874 {
875 basisMatrix = listSnapshotsModes[0];
876 quadratureWeights = EigenFunctions::reorderVectorFromDim(quadratureWeightsSubspaces[0], vectorial_dim);
877 evaluateWPU(P, basisMatrix, normalizingWeights, quadratureWeights);
878 cnpy::save(basisMatrix, folderMethod + "/basisMatrix.npy");
879 cnpy::save(quadratureWeights, folderMethod + "/quadratureWeights.npy");
880 cnpy::save(wPU, folderMethod + "/wPU.npy");
881 }
882 else
883 {
884 Eigen::MatrixXd quadratureWeightsMerged(n_nodes, n_subspaces);
885 for (int ith_subspace = 0; ith_subspace < n_subspaces; ith_subspace++)
886 {
887 quadratureWeightsSubspaces[ith_subspace] = EigenFunctions::reorderVectorFromDim(
888 quadratureWeightsSubspaces[ith_subspace], vectorial_dim);
889 cnpy::save(quadratureWeightsSubspaces[ith_subspace], folderMethod + "/quadratureWeights_subspace" +
890 std::to_string(ith_subspace+1) + ".npy");
891 quadratureWeightsMerged.col(ith_subspace) = quadratureWeightsSubspaces[ith_subspace];
892 }
893 cnpy::save(quadratureWeightsMerged, folderMethod + "/quadratureWeights.npy");
894 }
895 cnpy::save(normalizingWeights, folderMethod + "/normalizingWeights.npy");
896 cnpy::save(nodes, folderMethod + "/mp.npy");
897 P.makeCompressed();
898 cnpy::save(P, folderMethod + "/projectionMatrix.npz");
899 nodePoints().write();
900 }
901 else
902 {
903 Info << "Read ECP\n";
904 if (n_subspaces == 1)
905 {
906 cnpy::load(basisMatrix, folderMethod + "/basisMatrix.npy");
907 cnpy::load(quadratureWeights, folderMethod + "/quadratureWeights.npy");
908 cnpy::load(wPU, folderMethod + "/wPU.npy");
909 }
910 else
911 {
912 for (int ith_subspace = 0; ith_subspace < n_subspaces; ith_subspace++)
913 {
914 cnpy::load(quadratureWeightsSubspaces[ith_subspace], folderMethod + "/quadratureWeights_subspace" +
915 std::to_string(ith_subspace+1) + ".npy");
916 }
917 }
918 cnpy::load(normalizingWeights, folderMethod + "/normalizingWeights.npy");
919 cnpy::load(nodes, folderMethod + "/mp.npy");
920 cnpy::load(P, folderMethod + "/projectionMatrix.npz");
921 }
922}
923
924template<typename... SnapshotsLists>
925void HyperReduction<SnapshotsLists...>::initSeeds(Eigen::VectorXd mp_not_mask,
926 std::set<label> nodePointsSet)
927{
928 P.resize(n_cells * vectorial_dim, initialSeeds.rows() * vectorial_dim);
929 P.reserve(Eigen::VectorXi::Constant(initialSeeds.rows() *
930 vectorial_dim /*n_cols*/, 1 /*n_non_zero_elements*/));
931 unsigned int trueSeeds{0};
932
933 for (int i = 0; i < initialSeeds.rows(); i++)
934 {
935 unsigned int index = initialSeeds(i) % n_cells;
936
937 // check that there are not repeated nodes in initialSeeds
938 if (nodePointsSet.find(index) == nodePointsSet.end())
939 {
940 nodePointsSet.insert(index);
941 nodePoints->append(index);
942 nodes.conservativeResize(nodes.rows() + 1);
943 nodes(nodes.rows() - 1) = index;
944
945 for (unsigned int ith_field = 0; ith_field < vectorial_dim; ith_field++)
946 {
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;
950 }
951
952 trueSeeds++;
953 }
954 }
955
956 P.conservativeResize(P.rows(), nodePoints->size() * vectorial_dim);
957 P.makeCompressed();
958 M_Assert(nodePoints->size() <= n_nodes,
959 "Size of 'initialSeeds' is greater than 'n_nodes'");
960}
961
962template <typename... SnapshotsLists>
963void HyperReduction<SnapshotsLists...>::updateNodes(Eigen::SparseMatrix<double>
964 & P, label& ind, Eigen::VectorXd& mp_not_mask)
965{
966 nodePoints->append(ind);
967 nodes.conservativeResize(nodes.rows() + 1);
968 nodes(nodes.rows() - 1) = ind;
969 unsigned int last_col{0};
970
971 if (P.rows() == 0)
972 {
973 P.resize(n_cells * vectorial_dim, vectorial_dim);
974 }
975 else
976 {
977 last_col = P.cols();
978 P.resize(P.rows(), P.cols() + vectorial_dim);
979 }
980
981 int step = P.cols() / vectorial_dim - 1;
982
983 for (int ith_node = 0; ith_node <= step; ith_node++)
984 {
985 for (unsigned int ith_field = 0; ith_field < vectorial_dim; ith_field++)
986 {
987 P.insert(nodes[ith_node] * vectorial_dim + ith_field,
988 ith_node * vectorial_dim + ith_field) = 1;
989 }
990 }
991
992 for (unsigned int ith_field = 0; ith_field < vectorial_dim; ith_field++)
993 {
994 mp_not_mask(ind * vectorial_dim + ith_field) = 0;
995 }
996}
997
998template <typename... SnapshotsLists>
999void HyperReduction<SnapshotsLists...>::updateNodesRemove(Eigen::SparseMatrix<double>& P, Eigen::VectorXd weights)
1000{
1001 IOList<label> nodePointsCopy(IOobject("nodePoints", para->runTime.time().constant(), "../" + folderMethod,
1002 para->mesh, IOobject::READ_IF_PRESENT, IOobject::NO_WRITE));
1003 label n = 0;
1004 for (int ith_node = 0; ith_node < weights.size(); ith_node++)
1005 {
1006 if (weights(ith_node)>0)
1007 {
1008 nodePointsCopy.append(nodePoints()[ith_node]);
1009 nodes(n) = nodePoints()[ith_node];
1010 n += 1;
1011 }
1012 }
1013 nodePoints->transfer(nodePointsCopy);
1014 nodes.conservativeResize(n);
1015
1016 P.resize(n_cells * vectorial_dim, nodes.size() * vectorial_dim);
1017 P.reserve(nodes.size() * vectorial_dim);
1018 P.setZero();
1019
1020 int step = P.cols() / vectorial_dim - 1;
1021
1022 for (int ith_node = 0; ith_node <= step; ith_node++)
1023 {
1024 for (unsigned int ith_field = 0; ith_field < vectorial_dim; ith_field++)
1025 {
1026 P.insert(nodes[ith_node] * vectorial_dim + ith_field,
1027 ith_node * vectorial_dim + ith_field) = 1;
1028 }
1029 }
1030 // Not updating mp_not_mask : we are not choosing again a node that we removed
1031}
1032
1033template<typename... SnapshotsLists>
1034void HyperReduction<SnapshotsLists...>::computeLS(Eigen::MatrixXd& J,
1035 Eigen::MatrixXd& Jwhole, Eigen::VectorXd& b, Eigen::VectorXd& q, bool NNLS)
1036{
1037 J.resize(vectorial_dim * Jwhole.rows(), nodePoints->size());
1038 J = Jwhole(Eigen::indexing::all, nodes);
1039
1040 if (NNLS)
1041 {
1042 Eigen::NNLS<Eigen::MatrixXd> NNLS_solver(J);
1043 Eigen::VectorXd NNLS_solution = NNLS_solver.solve(q);
1044 quadratureWeights = NNLS_solution;
1045 }
1046
1047 else
1048 {
1049 quadratureWeights = J.colPivHouseholderQr().solve(q);
1050 }
1051 b = q.segment(0,b.size()) - J.middleRows(0,b.size()) * quadratureWeights;
1052 b = b / b.lpNorm<2>();
1053}
1054
1055template<typename... SnapshotsLists>
1056void HyperReduction<SnapshotsLists...>::initReshapeMat(
1057 Eigen::SparseMatrix<double>& reshapeMat)
1058{
1059 reshapeMat.resize(n_cells, n_cells * vectorial_dim);
1060 reshapeMat.reserve(Eigen::VectorXi::Constant(n_cells * vectorial_dim /*n_cols*/,
1061 1 /*n_non_zero_elements*/));
1062
1063 for (int i = 0; i < n_cells; i++)
1064 {
1065 for (unsigned int ith_field = 0; ith_field < vectorial_dim; ith_field++)
1066 {
1067 reshapeMat.insert(i, i * vectorial_dim + ith_field) = 1;
1068 }
1069 }
1070
1071 reshapeMat.makeCompressed();
1072}
1073
1074template <typename... SnapshotsLists>
1075template <typename SnapshotsList>
1076void HyperReduction<SnapshotsLists...>::stackSnapshots(SnapshotsList sList,
1077 Eigen::MatrixXd& snapshotsMatrix, Eigen::VectorXd& fieldWeights)
1078{
1079 unsigned int field_dim = get_field_dim<typename SnapshotsList::value_type>();
1080 Eigen::MatrixXd tmpSnapshots = Foam2Eigen::PtrList2Eigen(sList);
1081 // get volumes
1082 Eigen::VectorXd V = ITHACAutilities::getMassMatrixFV(sList[0]);
1083 double maxVal = std::sqrt(tmpSnapshots.colwise().lpNorm<2>().maxCoeff());
1084 fieldWeights.conservativeResize(fieldWeights.rows() + field_dim * n_cells);
1085 fieldWeights.tail(n_cells * field_dim) = V.array().sqrt().cwiseInverse() *
1086 maxVal;
1087 snapshotsMatrix.conservativeResize(snapshotsMatrix.rows() + n_cells * field_dim,
1088 n_snapshots);
1089 snapshotsMatrix.bottomRows(n_cells * field_dim) = tmpSnapshots;
1090}
1091
1092template <typename... SnapshotsLists>
1093template <typename SnapshotsList>
1094void HyperReduction<SnapshotsLists...>::stackSnapshotsBoundary(
1095 SnapshotsList sList, List<Eigen::MatrixXd>& snapshotsMatrixBoundary,
1096 List<Eigen::VectorXd>& fieldWeightsBoundary)
1097{
1098 unsigned int field_dim = get_field_dim<typename SnapshotsList::value_type>();
1099 List<Eigen::MatrixXd> tmpBoundarySnapshots = Foam2Eigen::PtrList2EigenBC(sList);
1100 List<double> maxVal;
1101 maxVal.resize(n_boundary_patches);
1102
1103 for (label id = 0; id < n_boundary_patches; id++)
1104 {
1105 maxVal[id] = std::sqrt(
1106 tmpBoundarySnapshots[id].colwise().lpNorm<2>().maxCoeff());
1107 // get surface areas
1108 Eigen::VectorXd S = Foam2Eigen::field2Eigen(
1109 sList[0].mesh().magSf().boundaryField()[id]);
1110 S = S.replicate(field_dim, 1);
1111 unsigned int bSize = field_dim * int(n_boundary_cells_list[id]);
1112 fieldWeightsBoundary[id].conservativeResize(fieldWeightsBoundary[id].rows() +
1113 bSize);
1114 fieldWeightsBoundary[id].tail(bSize) = S.array().pow(3 / 4) * maxVal[id];
1115 snapshotsMatrixBoundary[id].conservativeResize(
1116 snapshotsMatrixBoundary[id].rows() + bSize, n_snapshots);
1117 snapshotsMatrixBoundary[id].bottomRows(bSize) = tmpBoundarySnapshots[id];
1118 }
1119}
1120
1121template<typename... SnapshotsLists>
1122void HyperReduction<SnapshotsLists...>::evaluatePinv(Eigen::SparseMatrix<double>
1123 & Projector, Eigen::MatrixXd& Modes, Eigen::VectorXd& fieldWeights)
1124{
1125 assert(Projector.cols() > 0);
1126 Eigen::MatrixXd restricted = Projector.transpose() * Modes;
1127 Info << "S-Optimalty: " << s_optimality(restricted) << endl;
1128 pinvPU = ITHACAutilities::invertMatrix(restricted,
1129 para->ITHACAdict->lookupOrDefault<word>("InversionMethod",
1130 "completeOrthogonalDecomposition"));
1131 pinvPU = pinvPU * (Projector.transpose() *
1132 fieldWeights.array().cwiseInverse().matrix()).asDiagonal();
1133}
1134
1135template<typename... SnapshotsLists>
1136void HyperReduction<SnapshotsLists...>::evaluateWPU(Eigen::SparseMatrix<double>
1137 & Projector, Eigen::MatrixXd& Modes, Eigen::VectorXd& fieldWeights,
1138 Eigen::VectorXd& quadratureWeights)
1139{
1140 assert(Projector.cols() > 0);
1141 Eigen::VectorXd quadratureWeightsOrderedAsProjector(quadratureWeights.rows());
1142 unsigned int n_weightsPerField = quadratureWeights.rows() / vectorial_dim;
1143 unsigned int reorderIndex{0};
1144
1145 for (unsigned int ith_field = 0; ith_field < vectorial_dim; ith_field++)
1146 {
1147 for (unsigned int ith_weight = 0; ith_weight < n_weightsPerField; ith_weight++)
1148 {
1149 quadratureWeightsOrderedAsProjector(reorderIndex) = quadratureWeights(
1150 ith_weight + ith_field * n_weightsPerField);
1151 reorderIndex++;
1152 }
1153 }
1154
1155 wPU = quadratureWeightsOrderedAsProjector.transpose() *
1156 (Projector.transpose() *
1157 fieldWeights.array().cwiseInverse().matrix()).asDiagonal();
1158}
1159
1160template<typename... SnapshotsLists>
1161void HyperReduction<SnapshotsLists...>::generateSubmesh(label layers,
1162 const fvMesh& mesh)
1163{
1164 Info << "####### Extract submesh #######\n";
1165 ITHACAparameters* para(ITHACAparameters::getInstance());
1166 totalNodePoints = autoPtr<IOList<labelList >> (new IOList<labelList>(IOobject(
1167 "totalNodePoints", para->runTime.time().constant(), "../" + folderMethod,
1168 para->mesh, IOobject::READ_IF_PRESENT, IOobject::NO_WRITE)));
1169 uniqueNodePoints = autoPtr<IOList<label >> (new IOList<label>(IOobject(
1170 "uniqueNodePoints", para->runTime.time().constant(), "../" + folderMethod,
1171 para->mesh, IOobject::READ_IF_PRESENT, IOobject::NO_WRITE)));
1172 volScalarField Indici
1173 (
1174 IOobject
1175 (
1176 problemName + "_indices",
1177 mesh.time().timeName(),
1178 mesh,
1179 IOobject::NO_READ,
1180 IOobject::NO_WRITE
1181 ),
1182 mesh,
1183 dimensionedScalar(problemName + "_indices",
1184 dimensionSet(0, 0, 0, 0, 0, 0, 0), 0.0));
1185 submesh = autoPtr<fvMeshSubset>(new fvMeshSubset(mesh));
1186 bool offlineStage = !totalNodePoints().headerOk();
1187
1188 if (offlineStage)
1189 {
1190 List<label> indices;
1191
1192 for (label i = 0; i < nodePoints().size(); i++)
1193 {
1194 indices = ITHACAutilities::getIndices(mesh, nodePoints()[i], layers);
1195 totalNodePoints().append(indices);
1196 }
1197
1198 labelList a = ListListOps::combine<labelList>(totalNodePoints(),
1199 accessOp<labelList>());
1200 inplaceUniqueSort(a);
1201 uniqueNodePoints() = a;
1202 scalar zerodot25 = 0.25;
1203 ITHACAutilities::assignIF(Indici, zerodot25,
1204 uniqueNodePoints().List<label>::clone()());
1205 ITHACAutilities::assignONE(Indici, nodePoints());
1206 totalNodePoints().write();
1207 uniqueNodePoints().write();
1208 ITHACAstream::exportSolution(Indici, "1", folderMethod);
1209 }
1210
1211 submesh->setCellSubset(uniqueNodePoints());
1212 submesh->subMesh().fvSchemes::readOpt() = mesh.fvSchemes::readOpt();
1213 submesh->subMesh().fvSolution::readOpt() = mesh.fvSolution::readOpt();
1214 submesh->subMesh().fvSchemes::read();
1215 submesh->subMesh().fvSolution::read();
1216 std::cout.clear();
1217 localNodePoints = global2local(nodePoints(), submesh());
1218 n_cellsSubfields = submesh().cellMap().size();
1219 Info << "####### End extract submesh size = " << n_cellsSubfields <<
1220 " #######\n";
1221 createMasks(offlineStage);
1222 Info << "####### End create masks #######\n";
1223}
1224
1225template<typename... SnapshotsLists>
1226void HyperReduction<SnapshotsLists...>::createMasks(bool offlineStage)
1227{
1228 if (offlineStage)
1229 {
1230 field2submesh.resize(submesh().cellMap().size() * vectorial_dim,
1231 n_cells * vectorial_dim);
1232 field2submesh.reserve(Eigen::VectorXi::Constant(n_cells * vectorial_dim, 1));
1233
1234 for (int ith_subCell{0} ; ith_subCell < submesh().cellMap().size();
1235 ith_subCell++)
1236 {
1237 for (unsigned int ith_field = 0; ith_field < vectorial_dim; ith_field++)
1238 {
1239 field2submesh.insert(ith_subCell + ith_field * submesh().cellMap().size(),
1240 submesh().cellMap()[ith_subCell] + n_cells * ith_field) = 1;
1241 }
1242 }
1243
1244 field2submesh.makeCompressed();
1245 submesh2nodes.resize(nodePoints().size() * vectorial_dim,
1246 submesh().cellMap().size() * vectorial_dim);
1247 submesh2nodes.reserve(Eigen::VectorXi::Constant(submesh().cellMap().size() *
1248 vectorial_dim, 1));
1249 submesh2nodesMask.resize(nodePoints().size() * vectorial_dim);
1250 int index_col = 0;
1251
1252 for (int ith_node{0} ; ith_node < nodePoints().size(); ith_node++)
1253 {
1254 for (int ith_subCell{0} ; ith_subCell < submesh().cellMap().size();
1255 ith_subCell++)
1256 {
1257 if (nodePoints()[ith_node] == submesh().cellMap()[ith_subCell])
1258 {
1259 for (unsigned int ith_field = 0; ith_field < vectorial_dim; ith_field++)
1260 {
1261 submesh2nodes.insert(ith_node + nodePoints().size() * ith_field,
1262 index_col + ith_field * submesh().cellMap().size()) = 1;
1263 submesh2nodesMask(ith_node + nodePoints().size() * ith_field) = index_col +
1264 ith_field * submesh().cellMap().size();
1265 }
1266
1267 break;
1268 }
1269 else
1270 {
1271 index_col++;
1272 }
1273 }
1274
1275 index_col = 0;
1276 }
1277
1278 submesh2nodes.makeCompressed();
1279 cnpy::save(field2submesh, folderMethod + "/field2submesh.npz");
1280 cnpy::save(submesh2nodes, folderMethod + "/submesh2nodes.npz");
1281 cnpy::save(submesh2nodesMask, folderMethod + "/submesh2nodesMask.npy");
1282 }
1283 else
1284 {
1285 cnpy::load(field2submesh, folderMethod + "/field2submesh.npz");
1286 cnpy::load(submesh2nodes, folderMethod + "/submesh2nodes.npz");
1287 cnpy::load(submesh2nodesMask, folderMethod + "/submesh2nodesMask.npy");
1288 }
1289}
1290
1291template<typename... SnapshotsLists>
1292List<label> HyperReduction<SnapshotsLists...>::global2local(
1293 List<label>& points, fvMeshSubset& submesh)
1294{
1295 List<label> localPoints;
1296
1297 for (label i = 0; i < points.size(); i++)
1298 {
1299 for (label j = 0; j < submesh.cellMap().size(); j++)
1300 {
1301 if (submesh.cellMap()[j] == points[i])
1302 {
1303 localPoints.append(j);
1304 break;
1305 }
1306 }
1307 }
1308
1309 return localPoints;
1310}
1311
1312template<typename... SnapshotsLists>
1313Eigen::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
1335#endif
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:533
Foam2Eigen::field2Eigen
static Eigen::MatrixXd field2Eigen(GeometricField< Type, PatchField, GeoMesh > &field)
Convert a vector OpenFOAM field into an Eigen Vector.
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:273
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:665
HyperReduction::stackSnapshotsBoundary
void stackSnapshotsBoundary(SnapshotsList sList, List< Eigen::MatrixXd > &snapshotsMatrixBoundary, List< Eigen::VectorXd > &fieldWeightsBoundary)
TODO.
Definition hyperReduction.templates.H:1094
HyperReduction::snapshotsListTuple
SnapshotsListTuple snapshotsListTuple
The snapshots matrix containing the nonlinear function or operator.
Definition hyperReduction.H:135
HyperReduction::n_modes
label n_modes
The maximum number of modes to be considered.
Definition hyperReduction.H:117
HyperReduction::evaluateWPU
void evaluateWPU(Eigen::SparseMatrix< double > &Projector, Eigen::MatrixXd &Modes, Eigen::VectorXd &fieldWeights, Eigen::VectorXd &quadratureWeights)
Compute the pseudo-inverse of the matrix M restricted with the projector P.
Definition hyperReduction.templates.H:1136
HyperReduction::n_nodes
label n_nodes
The maximum number of modes to be considered.
Definition hyperReduction.H:120
HyperReduction::stackNames
void stackNames(SnapshotsList sList)
TODO.
Definition hyperReduction.H:354
HyperReduction::fieldNames
List< word > fieldNames
Names of the fields.
Definition hyperReduction.H:138
HyperReduction::generateSubmesh
void generateSubmesh(label layers, const fvMesh &mesh)
Compute the submesh common to all fields in SnapshotsLists.
Definition hyperReduction.templates.H:1161
HyperReduction::submesh
autoPtr< fvMeshSubset > submesh
Submesh of the HyperReduction method.
Definition hyperReduction.H:206
HyperReduction::offlineGappyDEIM
void offlineGappyDEIM(Eigen::MatrixXd &snapshotsModes, Eigen::VectorXd &normalizingWeights)
Methods implemented: 'GappyDEIM' from "DEIM, Chaturantabut, Saifon, and Danny C. Sorensen....
Definition hyperReduction.H:217
HyperReduction< PtrList< volScalarField > & >::updateNodes
void updateNodes(Eigen::SparseMatrix< double > &P, label &ind_max, Eigen::VectorXd &mp_not_mask)
Definition hyperReduction.templates.H:963
HyperReduction::evaluatePinv
void evaluatePinv(Eigen::SparseMatrix< double > &Projector, Eigen::MatrixXd &Modes, Eigen::VectorXd &fieldWeights)
Compute the pseudo-inverse of the matrix M restricted with the projector P.
Definition hyperReduction.templates.H:1122
HyperReduction::quadratureWeights
Eigen::VectorXd quadratureWeights
Quadrature weights. Ordered in the same order of matrix P.
Definition hyperReduction.H:202
HyperReduction::totalNodePoints
autoPtr< IOList< labelList > > totalNodePoints
List of label lists of the nodes and corresponding surrounding nodes.
Definition hyperReduction.H:156
HyperReduction::sumDimensions
void sumDimensions(double sum, SnapshotsList sList)
TODO.
Definition hyperReduction.H:378
HyperReduction::global2local
List< label > global2local(List< label > &points, fvMeshSubset &submesh)
Get local indices in the submesh from indices in the global ones.
Definition hyperReduction.templates.H:1292
HyperReduction::stackDimensions
void stackDimensions(SnapshotsList sList)
TODO.
Definition hyperReduction.H:366
HyperReduction::compute_vectorial_dim
constexpr unsigned int compute_vectorial_dim(LastList x)
TODO.
Definition hyperReduction.H:449
HyperReduction::getSnapMatrix
void getSnapMatrix(Eigen::MatrixXd &snapMatrix, Eigen::VectorXd &fieldWeights)
TODO.
Definition hyperReduction.templates.H:344
HyperReduction::localNodePoints
List< label > localNodePoints
Indices of the local node points in the subMesh.
Definition hyperReduction.H:168
HyperReduction::get_field_dim
constexpr unsigned int get_field_dim()
TODO.
Definition hyperReduction.H:438
HyperReduction::offlineECP
void offlineECP(Eigen::MatrixXd &snapshotsModes, Eigen::VectorXd &normalizingWeights)
Methods implemented: 'ECP' from "ECP, Hernandez, Joaquin Alberto, Manuel Alejandro Caicedo,...
Definition hyperReduction.H:231
HyperReduction::problemName
word problemName
The name of the non-linear function e.g. HR_method/residual.
Definition hyperReduction.H:132
HyperReduction::HyperReduction
HyperReduction(label n_modes, label n_nodes, Eigen::VectorXi initialSeeds, word problemName, SnapshotsLists &&...snapshotsLists)
Construct HyperReduction class, interpolation-based.
Definition hyperReduction.templates.H:40
HyperReduction::updateNodesRemove
void updateNodesRemove(Eigen::SparseMatrix< double > &P, Eigen::VectorXd weights)
Remove from the nodes the weights equal to 0 (used with NNLS).
Definition hyperReduction.templates.H:999
HyperReduction::n_cells
label n_cells
Int Number of Cells;.
Definition hyperReduction.H:144
HyperReduction< PtrList< volScalarField > & >::initReshapeMat
void initReshapeMat(Eigen::SparseMatrix< double > &reshapeMat)
Definition hyperReduction.templates.H:1056
HyperReduction< PtrList< volScalarField > & >::initSeeds
void initSeeds(Eigen::VectorXd mp_not_mask, std::set< label > nodePointsSet)
Definition hyperReduction.templates.H:925
HyperReduction::volumes
Eigen::VectorXd volumes
Volume of each cell.
Definition hyperReduction.H:147
HyperReduction::maskToOtherDim
Eigen::SparseMatrix< double > maskToOtherDim(int newDim)
Compute the sparse mask matrix P with a vectorial dimension different from the snapshots (enlarged or...
Definition hyperReduction.templates.H:1313
HyperReduction::s_optimality
static double s_optimality(Eigen::MatrixXd &A)
TODO.
Definition hyperReduction.H:472
HyperReduction::getModesSVD
void getModesSVD(SnapshotsListTuple &SnapshotsListTuple, Eigen::MatrixXd &modesSVD, Eigen::VectorXd &fieldWeights, bool saveModesFlag=false)
TODO.
Definition hyperReduction.templates.H:134
HyperReduction::n_snapshots
label n_snapshots
The length of the snapshots lists.
Definition hyperReduction.H:123
HyperReduction::uniqueNodePoints
autoPtr< IOList< label > > uniqueNodePoints
List of the unique indices of the nodes that define the submesh.
Definition hyperReduction.H:159
HyperReduction::saveModes
void saveModes(SnapshotsList sList, Eigen::MatrixXd &snapshotsMatrix, unsigned int &rowIndex, unsigned int &modeIndex, word folder)
TODO.
Definition hyperReduction.templates.H:373
HyperReduction::folderProblem
word folderProblem
Folder for the HR problem.
Definition hyperReduction.H:162
HyperReduction::stackSnapshots
void stackSnapshots(SnapshotsList sList, Eigen::MatrixXd &snapshotsMatrix, Eigen::VectorXd &fieldWeights)
TODO.
Definition hyperReduction.templates.H:1076
HyperReduction::computeLS
void computeLS(Eigen::MatrixXd &J, Eigen::MatrixXd &JWhole, Eigen::VectorXd &b, Eigen::VectorXd &q, bool NNLS)
TODO.
Definition hyperReduction.templates.H:1034
HyperReduction::n_cellsSubfields
label n_cellsSubfields
Int Number of Cells in submeshes;.
Definition hyperReduction.H:150
HyperReduction::createMasks
void createMasks(bool offlineStage=true)
TODO.
Definition hyperReduction.templates.H:1226
ITHACAparameters
Class for the definition of some general parameters, the parameters must be defined from the file ITH...
Definition ITHACAparameters.H:41
ITHACAparameters::getInstance
static ITHACAparameters * getInstance()
Gets an instance of ITHACAparameters, to be used if the instance is already existing.
Definition ITHACAparameters.C:56
ITHACAparameters::getInstance
static ITHACAparameters * getInstance(const fvMesh &mesh, Time &localTime)
Gets an instance of ITHACAparameters, to be used if the instance is not existing.
Definition ITHACAparameters.C:45
Modes
Implementation of a container class derived from PtrList.
Definition Modes.H:69
EigenFunctions::reorderVectorFromDim
Eigen::VectorXd reorderVectorFromDim(Eigen::VectorXd input, int dim)
Changes the order of the vector's elements given a dimension.
Definition EigenFunctions.C:146
EigenFunctions::max
T max(Eigen::SparseMatrix< T > &mat, label &ind_row, label &ind_col)
Find the maximum of a sparse Matrix (Useful for DEIM).
Definition EigenFunctions.H:300
ITHACAstream::exportSolution
void exportSolution(GeometricField< Type, PatchField, GeoMesh > &s, fileName subfolder, fileName folder, word fieldName)
Export a field to file in a certain folder and subfolder.
Definition ITHACAstream.C:992
ITHACAutilities::assignIF
void assignIF(GeometricField< Type, fvPatchField, volMesh > &s, Type value)
Assign internal field.
Definition ITHACAassign.C:31
ITHACAutilities::getIndices
List< label > getIndices(const fvMesh &mesh, int index, int layers)
Gets the indices of the cells around a certain cell.
Definition ITHACAgeometry.C:62
ITHACAutilities::invertMatrix
Eigen::MatrixXd invertMatrix(Eigen::MatrixXd &matrixToInvert, const word inversionMethod)
Invert a matrix given the method name in the ITHACAdict.
Definition ITHACAutilities.C:167
ITHACAutilities::assignONE
void assignONE(volScalarField &s, List< label > &L)
Assign one to volScalarField.
Definition ITHACAassign.C:72
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:313
ITHACAutilities::check_folder
bool check_folder(word folder)
Checks if a folder exists.
Definition ITHACAsystem.C:89
Generated by   doxygen 1.16.1