Loading...
Searching...
No Matches
24HyperReduction.C
Go to the documentation of this file.
1/*---------------------------------------------------------------------------*\
2 ██╗████████╗██╗ ██╗ █████╗ ██████╗ █████╗ ███████╗██╗ ██╗
3 ██║╚══██╔══╝██║ ██║██╔══██╗██╔════╝██╔══██╗ ██╔════╝██║ ██║
4 ██║ ██║ ███████║███████║██║ ███████║█████╗█████╗ ██║ ██║
5 ██║ ██║ ██╔══██║██╔══██║██║ ██╔══██║╚════╝██╔══╝ ╚██╗ ██╔╝
6 ██║ ██║ ██║ ██║██║ ██║╚██████╗██║ ██║ ██║ ╚████╔╝
7 ╚═╝ ╚═╝ ╚═╝ ╚═╝╚═╝ ╚═╝ ╚═════╝╚═╝ ╚═╝ ╚═╝ ╚═══╝
8
9 * In real Time Highly Advanced Computational Applications for Finite Volumes
10 * Copyright (C) 2017 by the ITHACA-FV authors
11-------------------------------------------------------------------------------
12License
13 This file is part of ITHACA-FV
14 ITHACA-FV is free software: you can redistribute it and/or modify
15 it under the terms of the GNU Lesser General Public License as published by
16 the Free Software Foundation, either version 3 of the License, or
17 (at your option) any later version.
18 ITHACA-FV is distributed in the hope that it will be useful,
19 but WITHOUT ANY WARRANTY; without even the implied warranty of
20 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
21 GNU Lesser General Public License for more details.
22 You should have received a copy of the GNU Lesser General Public License
23 along with ITHACA-FV. If not, see <http://www.gnu.org/licenses/>.
24Description
25 Example of the reconstruction of a non-linear function using the HyperReduction
26SourceFiles
27 24HyperReduction.C
28\*---------------------------------------------------------------------------*/
29
30#include "fvCFD.H"
31#include "fvOptions.H"
32#include "simpleControl.H"
33#include "simpleControl.H"
34#include "fvMeshSubset.H"
35#include "ITHACAutilities.H"
36#include "Foam2Eigen.H"
37#include "ITHACAstream.H"
38#include "ITHACAPOD.H"
39#include "hyperReduction.templates.H"
40#include <chrono>
41
42class HyperReduction_function : public
43 HyperReduction<PtrList<volScalarField> & >
44{
45 public:
46 using HyperReduction::HyperReduction;
47 static volScalarField evaluate_expression(volScalarField& S, Eigen::MatrixXd mu)
48 {
49 volScalarField yPos = S.mesh().C().component(vector::Y).ref();
50 volScalarField xPos = S.mesh().C().component(vector::X).ref();
51
52 for (auto i = 0; i < S.size(); i++)
53 {
54 S[i] = std::exp(- 2 * std::pow(xPos[i] - mu(0) - 1,
55 2) - 2 * std::pow(yPos[i] - mu(1) - 0.5, 2));
56 }
57
58 return S;
59 }
60
61 static Eigen::VectorXd evaluate_expression(volScalarField& S,
62 Eigen::MatrixXd mu, List<label> nodesList)
63 {
64 Eigen::VectorXd ret(nodesList.size());
65 volScalarField yPos = S.mesh().C().component(vector::Y).ref();
66 volScalarField xPos = S.mesh().C().component(vector::X).ref();
67
68 for (auto i = 0; i < nodesList.size(); i++)
69 {
70 ret[i] = std::exp(- 2 * std::pow(xPos[nodesList[i]] - mu(0) - 1,
71 2) - 2 * std::pow(yPos[nodesList[i]] - mu(1) - 0.5, 2));
72 }
73
74 return ret;
75 }
76
77 Eigen::VectorXd onlineCoeffs(volScalarField& S, Eigen::MatrixXd mu)
78 {
79 Eigen::VectorXd f = evaluate_expression(S, mu, localNodePoints);
80 return pinvPU * f;
81 }
82
83};
84
85class HyperReduction_vectorFunction : public
86 HyperReduction<PtrList<volVectorField> & >
87{
88 public:
89 using HyperReduction::HyperReduction;
90 static volVectorField evaluate_expression(volVectorField& S, Eigen::MatrixXd mu)
91 {
92 volScalarField yPos = S.mesh().C().component(vector::Y).ref();
93 volScalarField xPos = S.mesh().C().component(vector::X).ref();
94
95 for (auto i = 0; i < S.size(); i++)
96 {
97 S[i][0] = std::exp(- 2 * std::pow(xPos[i] - mu(0) - 1,
98 2) - 2 * std::pow(yPos[i] - mu(1) - 0.5, 2));
99 S[i][1] = std::exp(- 2 * std::pow(xPos[i] - mu(0) - 0.5,
100 2) - 2 * std::pow(yPos[i] - mu(1) - 0.5, 2));
101 S[i][2] = std::exp(- 2 * std::pow(xPos[i] - mu(0) - 1.5,
102 2) - 2 * std::pow(yPos[i] - mu(1) - 0., 2));
103 }
104
105 return S;
106 }
107
108 static Eigen::VectorXd evaluate_expression(volVectorField& S,
109 Eigen::MatrixXd mu, List<label> nodesList)
110 {
111 Eigen::VectorXd ret(nodesList.size() * 3);
112 volScalarField yPos = S.mesh().C().component(vector::Y).ref();
113 volScalarField xPos = S.mesh().C().component(vector::X).ref();
114
115 for (auto i = 0; i < nodesList.size(); i++)
116 {
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));
125 }
126
127 return ret;
128 }
129
130 Eigen::VectorXd onlineCoeffs(volVectorField& S, Eigen::MatrixXd mu)
131 {
132 Eigen::VectorXd f = evaluate_expression(S, mu, localNodePoints);
133 return pinvPU * f;
134 }
135
136};
137
138class HyperReduction_vectorScalarFunction : public
139 HyperReduction<PtrList<volVectorField> &, PtrList<volScalarField> & >
140{
141 public:
142 using HyperReduction::HyperReduction;
143 static std::pair<volVectorField, volScalarField> evaluate_expression(
144 volVectorField& V, volScalarField& S, Eigen::MatrixXd mu)
145 {
146 volScalarField yPos = S.mesh().C().component(vector::Y).ref();
147 volScalarField xPos = S.mesh().C().component(vector::X).ref();
148
149 for (auto i = 0; i < S.size(); i++)
150 {
151 V[i][0] = std::exp(- 2 * std::pow(xPos[i] - mu(0) - 1,
152 2) - 2 * std::pow(yPos[i] - mu(1) - 0.5, 2));
153 V[i][1] = std::exp(- 2 * std::pow(xPos[i] - mu(0) - 0.5,
154 2) - 2 * std::pow(yPos[i] - mu(1) - 0.5, 2));
155 V[i][2] = std::exp(- 2 * std::pow(xPos[i] - mu(0) - 1.5,
156 2) - 2 * std::pow(yPos[i] - mu(1) - 0., 2));
157 S[i] = std::exp(- 2 * std::pow(xPos[i] - mu(0) - 1,
158 2) - 2 * std::pow(yPos[i] - mu(1) - 0.5, 2));
159 }
160
161 return std::make_pair(V, S);
162 }
163
164 static Eigen::VectorXd evaluate_expression(volScalarField& S,
165 Eigen::MatrixXd mu, List<label> nodesList)
166 {
167 Eigen::VectorXd ret(nodesList.size() * 4);
168 volScalarField yPos = S.mesh().C().component(vector::Y).ref();
169 volScalarField xPos = S.mesh().C().component(vector::X).ref();
170
171 for (auto i = 0; i < nodesList.size(); i++)
172 {
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));
184 }
185
186 return ret;
187 }
188
189 Eigen::VectorXd onlineCoeffs(volScalarField& S, Eigen::MatrixXd mu)
190 {
191 Eigen::VectorXd f = evaluate_expression(S, mu, localNodePoints);
192 return pinvPU * f;
193 }
194
195};
196
197void test_scalar(ITHACAparameters* para, Foam::fvMesh& mesh,
198 Foam::Time& runTime)
199{
200 int n_modes = para->ITHACAdict->lookupOrDefault<int>("Modes", 15);
201 int n_nodes = para->ITHACAdict->lookupOrDefault<int>("Nodes", 15);
202 simpleControl simple(mesh);
203 word methodName = para->ITHACAdict->lookupOrDefault<word>("HyperReduction",
204 "GappyDEIM");
205#include "createFields.H"
206 // List of volScalarField where the snapshots are stored
207 PtrList<volScalarField> Sp;
208 // Non linear function
209 volScalarField S
210 (
211 IOobject
212 (
213 "S",
214 runTime.timeName(),
215 mesh,
216 IOobject::NO_READ,
217 IOobject::AUTO_WRITE
218 ),
219 T.mesh(),
220 dimensionedScalar("zero", dimensionSet(0, 0, -1, 1, 0, 0, 0), 0)
221 );
222 // Parameters used to train the non-linear function
223 Eigen::MatrixXd pars;
224 cnpy::load(pars, "trainingPars.npy");
225
226 // Perform the offline phase
227 for (int i = 0; i < 100; i++)
228 {
229 HyperReduction_function::evaluate_expression(S, pars.row(i));
230 Sp.append((S).clone());
231 ITHACAstream::exportSolution(S, name(i + 1), "./ITHACAoutput/Offline/");
232 }
233
234 // Define new online parameters
235 Eigen::MatrixXd parTest;// = ITHACAutilities::rand(100, 2, -0.5, 0.5);
236 cnpy::load(parTest, "testingPars.npy");
237 // Create HyperReduction object with given number of basis functions
238 Eigen::VectorXi initSeeds;
239 // initSeeds = cnpy::load(initSeeds, "./mp.npy");//load mp to test initialSeeds
240
241 if (methodName == "GappyDEIM")
242 {
243 HyperReduction_function c(n_modes, n_nodes, initSeeds, "Gaussian_function", Sp);
244 Eigen::MatrixXd snapshotsModes;
245 Eigen::VectorXd normalizingWeights;
246 c.getModesSVD(c.snapshotsListTuple, snapshotsModes, normalizingWeights);
247 c.offlineGappyDEIM(snapshotsModes, normalizingWeights);
248 // Generate the submeshes with the depth of the layer
249 unsigned int layers{2};
250 c.generateSubmesh(layers, mesh);
251 auto sfield = c.interpolateField<volScalarField>(Sp[0]);
252 PtrList<volScalarField> testTrueFields;
253 PtrList<volScalarField> testRecFields;
254
255 for (unsigned int idTest = 0; idTest < parTest.rows(); idTest++)
256 {
257 // Online evaluation of the non linear function
258 Eigen::VectorXd aprfield = c.renormalizedBasisMatrix * c.onlineCoeffs(sfield(),
259 parTest.row(idTest));
260 // Transform to an OpenFOAM field and export
261 volScalarField S2("S_online", Foam2Eigen::Eigen2field(S, aprfield));
262 // Evaluate the full order function and export it
263 HyperReduction_function::evaluate_expression(S, parTest.row(idTest));
264 testRecFields.append(S2.clone());
265 testTrueFields.append(S.clone());
266 ITHACAstream::exportSolution(S2, name(idTest), "./ITHACAoutput/Online/");
267 ITHACAstream::exportSolution(S, name(idTest), "./ITHACAoutput/Online/");
268 }
269
270 // Compute the L2 error and print it
271 auto err = ITHACAutilities::errorL2Rel(testRecFields, testTrueFields);
272 Info << "GappyDEIM max relative error: " << err.maxCoeff() << endl;
273 Info << "GappyDEIM mean relative error: " << err.array().mean() << endl;
274 }
275 else if (methodName == "ECP")
276 {
277 HyperReduction_function c(n_modes, n_nodes, initSeeds, "Gaussian_function", Sp);
278 bool saveOFModes{true};
279 Eigen::MatrixXd snapshotsModes;
280 Eigen::VectorXd normalizingWeights;
281 c.getModesSVD(c.snapshotsListTuple, snapshotsModes, normalizingWeights,
282 saveOFModes);
283 c.offlineECP(snapshotsModes, normalizingWeights);
284 // Generate the submeshes with the depth of the layer
285 unsigned int layers{2};
286 c.generateSubmesh(layers, mesh);
287 auto sfield = c.interpolateField<volScalarField>(Sp[0]);
288 Eigen::VectorXd results(parTest.rows());
289
290 for (unsigned int idTest = 0; idTest < parTest.rows(); idTest++)
291 {
292 // Online evaluation of the non linear function
293 Eigen::VectorXd f = c.evaluate_expression(sfield(), parTest.row(idTest),
294 c.localNodePoints);
295 Eigen::VectorXd ff = Foam2Eigen::field2Eigen(c.evaluate_expression(Sp[0],
296 parTest.row(idTest)));
297 double trueIntegral = (normalizingWeights.cwiseInverse().asDiagonal() *
298 ff).array().sum();
299 double testIntegral = (c.wPU * f).array().sum();
300 // Info << "Integral: " << trueIntegral << endl;
301 // Info << "Reconstructed Integral: " << testIntegral << endl;
302 results(idTest) = trueIntegral - testIntegral;
303 }
304
305 Info << "ECP max error: " << results.cwiseAbs().maxCoeff() << endl;
306 Info << "ECP mean error: " << results.cwiseAbs().mean() << endl;
307 }
308}
309
310void test_vector(ITHACAparameters* para, Foam::fvMesh& mesh,
311 Foam::Time& runTime)
312{
313 int n_modes = para->ITHACAdict->lookupOrDefault<int>("Modes", 15);
314 int n_nodes = para->ITHACAdict->lookupOrDefault<int>("Nodes", 15);
315 simpleControl simple(mesh);
316 word methodName = para->ITHACAdict->lookupOrDefault<word>("HyperReduction",
317 "GappyDEIM");
318#include "createFields.H"
319 // List of volVectorField where the snapshots are stored
320 PtrList<volVectorField> Sp;
321 // Non linear function
322 volVectorField S
323 (
324 IOobject
325 (
326 "S",
327 runTime.timeName(),
328 mesh,
329 IOobject::NO_READ,
330 IOobject::AUTO_WRITE
331 ),
332 T.mesh(),
333 dimensionedVector("zero", dimensionSet(0, 0, -1, 1, 0, 0, 0), vector(0, 0, 0))
334 );
335 // Parameters used to train the non-linear function
336 Eigen::MatrixXd pars;
337 cnpy::load(pars, "trainingPars.npy");
338
339 // Perform the offline phase
340 for (int i = 0; i < 100; i++)
341 {
342 HyperReduction_vectorFunction::evaluate_expression(S, pars.row(i));
343 Sp.append((S).clone());
344 ITHACAstream::exportSolution(S, name(i + 1), "./ITHACAoutput/Offline/");
345 }
346
347 // Define new online parameters
348 Eigen::MatrixXd parTest;// = ITHACAutilities::rand(100, 2, -0.5, 0.5);
349 cnpy::load(parTest, "testingPars.npy");
350 // Create HyperReduction object with given number of basis functions
351 Eigen::VectorXi initSeeds;
352 // initSeeds = cnpy::load(initSeeds, "./mp.npy");//load mp to test initialSeeds
353
354 if (methodName == "GappyDEIM")
355 {
356 HyperReduction_vectorFunction c(n_modes, n_nodes, initSeeds,
357 "Gaussian_function", Sp);
358 Eigen::MatrixXd snapshotsModes;
359 Eigen::VectorXd normalizingWeights;
360 c.getModesSVD(c.snapshotsListTuple, snapshotsModes, normalizingWeights);
361 c.offlineGappyDEIM(snapshotsModes, normalizingWeights);
362 // Generate the submeshes with the depth of the layer
363 unsigned int layers{2};
364 c.generateSubmesh(layers, mesh);
365 auto sfield = c.interpolateField<volVectorField>(Sp[0]);
366 PtrList<volVectorField> testTrueFields;
367 PtrList<volVectorField> testRecFields;
368
369 for (unsigned int idTest = 0; idTest < parTest.rows(); idTest++)
370 {
371 // Online evaluation of the non linear function
372 Eigen::VectorXd aprfield = c.renormalizedBasisMatrix * c.onlineCoeffs(sfield(),
373 parTest.row(idTest));
374 // Transform to an OpenFOAM field and export
375 volVectorField S2("S_online", Foam2Eigen::Eigen2field(S, aprfield));
376 // Evaluate the full order function and export it
377 HyperReduction_vectorFunction::evaluate_expression(S, parTest.row(idTest));
378 testRecFields.append(S2.clone());
379 testTrueFields.append(S.clone());
380 ITHACAstream::exportSolution(S2, name(idTest), "./ITHACAoutput/Online/");
381 ITHACAstream::exportSolution(S, name(idTest), "./ITHACAoutput/Online/");
382 }
383
384 // Compute the L2 error and print it
385 auto err = ITHACAutilities::errorL2Rel(testRecFields, testTrueFields);
386 Info << "GappyDEIM max relative error: " << err.maxCoeff() << endl;
387 Info << "GappyDEIM mean relative error: " << err.array().mean() << endl;
388 }
389 else if (methodName == "ECP")
390 {
391 HyperReduction_vectorFunction c(n_modes, n_nodes, initSeeds,
392 "Gaussian_function", Sp);
393 bool saveOFModes{true};
394 Eigen::MatrixXd snapshotsModes;
395 Eigen::VectorXd normalizingWeights;
396 c.getModesSVD(c.snapshotsListTuple, snapshotsModes, normalizingWeights,
397 saveOFModes);
398 c.offlineECP(snapshotsModes, normalizingWeights);
399 // Generate the submeshes with the depth of the layer
400 unsigned int layers{2};
401 c.generateSubmesh(layers, mesh);
402 auto sfield = c.interpolateField<volVectorField>(Sp[0]);
403 Eigen::VectorXd results(parTest.rows());
404
405 for (unsigned int idTest = 0; idTest < parTest.rows(); idTest++)
406 {
407 // Online evaluation of the non linear function
408 Eigen::VectorXd f = c.evaluate_expression(sfield(), parTest.row(idTest),
409 c.localNodePoints);
410 auto wholeField = c.evaluate_expression(Sp[0], parTest.row(idTest));
411 Eigen::VectorXd ff = Foam2Eigen::field2Eigen(wholeField);
412 double trueIntegral = (normalizingWeights.cwiseInverse().asDiagonal() *
413 ff).array().sum();
414 double testIntegral = (c.wPU * f).array().sum();
415 // Info << "Integral: " << trueIntegral << endl;
416 // Info << "Reconstructed Integral: " << testIntegral << endl;
417 results(idTest) = trueIntegral - testIntegral;
418 }
419
420 Info << "ECP max error: " << results.cwiseAbs().maxCoeff() << endl;
421 Info << "ECP mean error: " << results.cwiseAbs().mean() << endl;
422 }
423}
424
425void test_vector_scalar(ITHACAparameters* para, Foam::fvMesh& mesh,
426 Foam::Time& runTime)
427{
428 int n_modes = para->ITHACAdict->lookupOrDefault<int>("Modes", 15);
429 int n_nodes = para->ITHACAdict->lookupOrDefault<int>("Nodes", 15);
430 simpleControl simple(mesh);
431 word methodName = para->ITHACAdict->lookupOrDefault<word>("HyperReduction",
432 "GappyDEIM");
433#include "createFields.H"
434 // List of volVectorField where the snapshots are stored
435 PtrList<volVectorField> Vp;
436 PtrList<volScalarField> Sp;
437 // Non linear function
438 volVectorField V
439 (
440 IOobject
441 (
442 "V",
443 runTime.timeName(),
444 mesh,
445 IOobject::NO_READ,
446 IOobject::AUTO_WRITE
447 ),
448 T.mesh(),
449 dimensionedVector("zero", dimensionSet(0, 0, -1, 1, 0, 0, 0), vector(0, 0, 0))
450 );
451 // Non linear function
452 volScalarField S
453 (
454 IOobject
455 (
456 "S",
457 runTime.timeName(),
458 mesh,
459 IOobject::NO_READ,
460 IOobject::AUTO_WRITE
461 ),
462 T.mesh(),
463 dimensionedScalar("zero", dimensionSet(0, 0, -1, 1, 0, 0, 0), 0)
464 );
465 // Parameters used to train the non-linear function
466 Eigen::MatrixXd pars;
467 cnpy::load(pars, "trainingPars.npy");
468
469 // Perform the offline phase
470 for (int i = 0; i < 100; i++)
471 {
472 HyperReduction_vectorScalarFunction::evaluate_expression(V, S, pars.row(i));
473 Vp.append((V).clone());
474 Sp.append((S).clone());
475 ITHACAstream::exportSolution(V, name(i + 1), "./ITHACAoutput/Offline/");
476 ITHACAstream::exportSolution(S, name(i + 1), "./ITHACAoutput/Offline/");
477 }
478
479 // Define new online parameters
480 Eigen::MatrixXd parTest;// = ITHACAutilities::rand(100, 2, -0.5, 0.5);
481 cnpy::load(parTest, "testingPars.npy");
482 // Create HyperReduction object with given number of basis functions
483 Eigen::VectorXi initSeeds;
484 // initSeeds = cnpy::load(initSeeds, "./mp.npy");//load mp to test initialSeeds
485
486 if (methodName == "GappyDEIM")
487 {
488 HyperReduction_vectorScalarFunction c(n_modes, n_nodes, initSeeds,
489 "Gaussian_function", Vp, Sp);
490 Eigen::MatrixXd snapshotsModes;
491 Eigen::VectorXd normalizingWeights;
492 c.getModesSVD(c.snapshotsListTuple, snapshotsModes, normalizingWeights);
493 c.offlineGappyDEIM(snapshotsModes, normalizingWeights);
494 // Generate the submeshes with the depth of the layer
495 unsigned int layers{2};
496 c.generateSubmesh(layers, mesh);
497 auto vfield = c.interpolateField<volVectorField>(Vp[0]);
498 auto sfield = c.interpolateField<volScalarField>(Sp[0]);
499 PtrList<volVectorField> testTrueFieldsV;
500 PtrList<volVectorField> testRecFieldsV;
501 PtrList<volScalarField> testTrueFieldsS;
502 PtrList<volScalarField> testRecFieldsS;
503
504 for (unsigned int idTest = 0; idTest < parTest.rows(); idTest++)
505 {
506 // Online evaluation of the non linear function
507 auto theta = c.onlineCoeffs(sfield(), parTest.row(idTest));
508 Eigen::VectorXd aprfield = c.renormalizedBasisMatrix * theta;
509 Eigen::VectorXd recvec = aprfield.head(3 * S.size());
510 Eigen::VectorXd recsca = aprfield.tail(S.size());
511 // Transform to an OpenFOAM field and export
512 volVectorField V2("V_online", Foam2Eigen::Eigen2field(V, recvec));
513 volScalarField S2("S_online", Foam2Eigen::Eigen2field(S, recsca));
514 // Evaluate the full order function and export it
515 HyperReduction_vectorScalarFunction::evaluate_expression(V, S,
516 parTest.row(idTest));
517 testRecFieldsV.append(V2.clone());
518 testTrueFieldsV.append(V.clone());
519 ITHACAstream::exportSolution(V2, name(idTest), "./ITHACAoutput/Online/");
520 ITHACAstream::exportSolution(V, name(idTest), "./ITHACAoutput/Online/");
521 testRecFieldsS.append(S2.clone());
522 testTrueFieldsS.append(S.clone());
523 ITHACAstream::exportSolution(S2, name(idTest), "./ITHACAoutput/Online/");
524 ITHACAstream::exportSolution(S, name(idTest), "./ITHACAoutput/Online/");
525 }
526
527 // Compute the L2 error and print it
528 auto errV = ITHACAutilities::errorL2Rel(testRecFieldsV, testTrueFieldsV);
529 Info << "GappyDEIM vector field max relative error: " << errV.maxCoeff() <<
530 endl;
531 Info << "GappyDEIM vector field mean relative error: " << errV.array().mean() <<
532 endl;
533 auto errS = ITHACAutilities::errorL2Rel(testRecFieldsS, testTrueFieldsS);
534 Info << "GappyDEIM scalar max relative error: " << errS.maxCoeff() << endl;
535 Info << "GappyDEIM scalar mean relative error: " << errS.array().mean() << endl;
536 }
537 else if (methodName == "ECP")
538 {
539 HyperReduction_vectorScalarFunction c(n_modes, n_nodes, initSeeds,
540 "Gaussian_function", Vp, Sp);
541 bool saveOFModes{true};
542 Eigen::MatrixXd snapshotsModes;
543 Eigen::VectorXd normalizingWeights;
544 c.getModesSVD(c.snapshotsListTuple, snapshotsModes, normalizingWeights,
545 saveOFModes);
546 c.offlineECP(snapshotsModes, normalizingWeights);
547 // Generate the submeshes with the depth of the layer
548 unsigned int layers{2};
549 c.generateSubmesh(layers, mesh);
550 auto vfield = c.interpolateField<volVectorField>(Vp[0]);
551 auto sfield = c.interpolateField<volScalarField>(Sp[0]);
552 Eigen::VectorXd results(parTest.rows());
553
554 for (unsigned int idTest = 0; idTest < parTest.rows(); idTest++)
555 {
556 // Online evaluation of the non linear function
557 Eigen::VectorXd f = c.evaluate_expression(sfield(), parTest.row(idTest),
558 c.localNodePoints);
559 HyperReduction_vectorScalarFunction::evaluate_expression(V, S,
560 parTest.row(idTest));
561 Eigen::VectorXd ffV = Foam2Eigen::field2Eigen(V);
562 Eigen::VectorXd ffS = Foam2Eigen::field2Eigen(S);
563 Eigen::VectorXd ff(ffV.rows() + ffS.rows());
564 ff.head(ffV.rows()) = ffV;
565 ff.tail(ffS.rows()) = ffS;
566 double trueIntegral = (normalizingWeights.cwiseInverse().asDiagonal() *
567 ff).array().sum();
568 double testIntegral = (c.wPU * f).array().sum();
569 // Info << "Integral: " << trueIntegral << endl;
570 // Info << "Reconstructed Integral: " << testIntegral << endl;
571 results(idTest) = trueIntegral - testIntegral;
572 }
573
574 Info << "ECP max error: " << results.cwiseAbs().maxCoeff() << endl;
575 Info << "ECP mean error: " << results.cwiseAbs().mean() << endl;
576 }
577}
578
579int main(int argc, char* argv[])
580{
581 // load stage to perform
582 argList::addOption("test", "scalar", "Perform scalar test");
583 argList::addOption("test", "vector", "Perform vector test");
584 argList::addOption("test", "vs", "Perform (vector, scalar) test");
585 // add options for parallel run
586 HashTable<string> validParOptions;
587 validParOptions.set
588 (
589 "test",
590 "scalar"
591 );
592 validParOptions.set
593 (
594 "test",
595 "vector"
596 );
597 validParOptions.set
598 (
599 "test",
600 "vs"
601 );
602 // Read parameters from ITHACAdict file
603#include "setRootCase.H"
604 Foam::Time runTime(Foam::Time::controlDictName, args);
605 Foam::fvMesh mesh
606 (
607 Foam::IOobject
608 (
609 Foam::fvMesh::defaultRegion,
610 runTime.timeName(),
611 runTime,
612 Foam::IOobject::MUST_READ
613 )
614 );
615 ITHACAparameters* para = ITHACAparameters::getInstance(mesh, runTime);
616
617 if (args.get("test").match("scalar"))
618 {
619 Info << "Init scalar testcase\n";
620 test_scalar(para, mesh, runTime);
621 }
622 else if (args.get("test").match("vector"))
623 {
624 Info << "Init vector testcase\n";
625 test_vector(para, mesh, runTime);
626 }
627 else if (args.get("test").match("vs"))
628 {
629 Info << "Init vector scalar testcase\n";
630 test_vector_scalar(para, mesh, runTime);
631 }
632 else
633 {
634 Info << "Pass '-test scalar', '-test vector', '-test vs'" << endl;
635 }
636
637 return 0;
638}
main
int main(int argc, char *argv[])
Definition 24HyperReduction.C:579
test_vector
void test_vector(ITHACAparameters *para, Foam::fvMesh &mesh, Foam::Time &runTime)
Definition 24HyperReduction.C:310
test_vector_scalar
void test_vector_scalar(ITHACAparameters *para, Foam::fvMesh &mesh, Foam::Time &runTime)
Definition 24HyperReduction.C:425
test_scalar
void test_scalar(ITHACAparameters *para, Foam::fvMesh &mesh, Foam::Time &runTime)
Definition 24HyperReduction.C:197
Foam2Eigen.H
Header file of the Foam2Eigen class.
ITHACAPOD.H
Header file of the ITHACAPOD class.
mesh
Foam::fvMesh & mesh
Definition createMesh.H:47
runTime
Foam::Time & runTime
Definition createTime.H:33
ITHACAstream.H
Header file of the ITHACAstream class, it contains the implementation of several methods for input ou...
ITHACAutilities.H
Header file of the ITHACAutilities namespace.
para
ITHACAparameters * para
Definition NLsolve.H:40
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
HyperReduction_function
Definition 24HyperReduction.C:44
HyperReduction_function::onlineCoeffs
Eigen::VectorXd onlineCoeffs(volScalarField &S, Eigen::MatrixXd mu)
Definition 24HyperReduction.C:77
HyperReduction_function::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_function::evaluate_expression
static volScalarField evaluate_expression(volScalarField &S, Eigen::MatrixXd mu)
Definition 24HyperReduction.C:47
HyperReduction_function::evaluate_expression
static Eigen::VectorXd evaluate_expression(volScalarField &S, Eigen::MatrixXd mu, List< label > nodesList)
Definition 24HyperReduction.C:61
HyperReduction_vectorFunction
Definition 24HyperReduction.C:87
HyperReduction_vectorFunction::evaluate_expression
static volVectorField evaluate_expression(volVectorField &S, Eigen::MatrixXd mu)
Definition 24HyperReduction.C:90
HyperReduction_vectorFunction::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_vectorFunction::onlineCoeffs
Eigen::VectorXd onlineCoeffs(volVectorField &S, Eigen::MatrixXd mu)
Definition 24HyperReduction.C:130
HyperReduction_vectorFunction::evaluate_expression
static Eigen::VectorXd evaluate_expression(volVectorField &S, Eigen::MatrixXd mu, List< label > nodesList)
Definition 24HyperReduction.C:108
HyperReduction_vectorScalarFunction
Definition 24HyperReduction.C:140
HyperReduction_vectorScalarFunction::onlineCoeffs
Eigen::VectorXd onlineCoeffs(volScalarField &S, Eigen::MatrixXd mu)
Definition 24HyperReduction.C:189
HyperReduction_vectorScalarFunction::evaluate_expression
static Eigen::VectorXd evaluate_expression(volScalarField &S, Eigen::MatrixXd mu, List< label > nodesList)
Definition 24HyperReduction.C:164
HyperReduction_vectorScalarFunction::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_vectorScalarFunction::evaluate_expression
static std::pair< volVectorField, volScalarField > evaluate_expression(volVectorField &V, volScalarField &S, Eigen::MatrixXd mu)
Definition 24HyperReduction.C:143
HyperReduction::snapshotsListTuple
SnapshotsListTuple snapshotsListTuple
The snapshots matrix containing the nonlinear function or operator.
Definition hyperReduction.H:132
HyperReduction::generateSubmesh
void generateSubmesh(label layers, const fvMesh &mesh)
Compute the submesh common to all fields in SnapshotsLists.
Definition hyperReduction.templates.H:983
HyperReduction::interpolateField
autoPtr< FieldType > interpolateField(const FieldType &field)
TODO.
Definition hyperReduction.H:397
HyperReduction::renormalizedBasisMatrix
Eigen::MatrixXd renormalizedBasisMatrix
Renormalized basis of HR.
Definition hyperReduction.H:172
HyperReduction::offlineGappyDEIM
void offlineGappyDEIM(Eigen::MatrixXd &snapshotsModes, Eigen::VectorXd &normalizingWeights)
Methods implemented: 'GappyDEIM' from "DEIM, Chaturantabut, Saifon, and Danny C. Sorensen....
Definition hyperReduction.H:210
HyperReduction::offlineECP
void offlineECP(Eigen::MatrixXd &snapshotsModes, Eigen::VectorXd &normalizingWeights)
Methods implemented: 'ECP' from "ECP, Hernandez, Joaquin Alberto, Manuel Alejandro Caicedo,...
Definition hyperReduction.H:224
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::getModesSVD
void getModesSVD(SnapshotsListTuple &SnapshotsListTuple, Eigen::MatrixXd &modesSVD, Eigen::VectorXd &fieldWeights, bool saveModesFlag=false)
TODO.
Definition hyperReduction.templates.H:133
HyperReduction::wPU
Eigen::MatrixXd wPU
Definition hyperReduction.H:178
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
T
volScalarField & T
Definition createT.H:46
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:904
ITHACAutilities::errorL2Rel
double errorL2Rel(GeometricField< T, fvPatchField, volMesh > &field1, GeometricField< T, fvPatchField, volMesh > &field2, List< label > *labels)
Computes the relative error between two geometric Fields in L2 norm.
Definition ITHACAerror.C:369
cnpy::load
Eigen::Matrix< T, -1, dim > load(Eigen::Matrix< T, -1, dim > &mat, const std::string fname)
simple
simpleControl simple(mesh)
i
label i
Definition pEqn.H:46
xPos
volScalarField xPos
Definition createFields.H:38
yPos
volScalarField yPos
Definition createFields.H:37
S
volScalarField S(IOobject("S", runTime.timeName(), mesh, IOobject::NO_READ, IOobject::AUTO_WRITE), T.mesh(), dimensionedScalar("zero", dimensionSet(0, 0, -1, 1, 0, 0, 0), 0))
Generated by   doxygen 1.13.2