Loading...
Searching...
No Matches
24HyperReduction.C
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(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));
123 }
124
125 return ret;
126 }
127
128 Eigen::VectorXd onlineCoeffs(volVectorField& S, Eigen::MatrixXd mu)
129 {
130 Eigen::VectorXd f = evaluate_expression(S, mu, localNodePoints);
131 return pinvPU * f;
132 }
133
134};
135
136class HyperReduction_vectorScalarFunction : public
137 HyperReduction<PtrList<volVectorField> &, PtrList<volScalarField> & >
138{
139 public:
140 using HyperReduction::HyperReduction;
141 static std::pair<volVectorField, volScalarField> evaluate_expression(
142 volVectorField& V, volScalarField& S, Eigen::MatrixXd mu)
143 {
144 volScalarField yPos = S.mesh().C().component(vector::Y).ref();
145 volScalarField xPos = S.mesh().C().component(vector::X).ref();
146
147 for (auto i = 0; i < S.size(); i++)
148 {
149 V[i][0] = std::exp(- 2 * std::pow(xPos[i] - mu(0) - 1,
150 2) - 2 * std::pow(yPos[i] - mu(1) - 0.5, 2));
151 V[i][1] = std::exp(- 2 * std::pow(xPos[i] - mu(0) - 0.5,
152 2) - 2 * std::pow(yPos[i] - mu(1) - 0.5, 2));
153 V[i][2] = std::exp(- 2 * std::pow(xPos[i] - mu(0) - 1.5,
154 2) - 2 * std::pow(yPos[i] - mu(1) - 0., 2));
155 S[i] = std::exp(- 2 * std::pow(xPos[i] - mu(0) - 1,
156 2) - 2 * std::pow(yPos[i] - mu(1) - 0.5, 2));
157 }
158
159 return std::make_pair(V, S);
160 }
161
162 static Eigen::VectorXd evaluate_expression(volScalarField& S,
163 Eigen::MatrixXd mu, List<label> nodesList)
164 {
165 Eigen::VectorXd ret(nodesList.size() * 4);
166 volScalarField yPos = S.mesh().C().component(vector::Y).ref();
167 volScalarField xPos = S.mesh().C().component(vector::X).ref();
168
169 for (auto i = 0; i < nodesList.size(); i++)
170 {
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));
179 }
180
181 return ret;
182 }
183
184 Eigen::VectorXd onlineCoeffs(volScalarField& S, Eigen::MatrixXd mu)
185 {
186 Eigen::VectorXd f = evaluate_expression(S, mu, localNodePoints);
187 return pinvPU * f;
188 }
189
190};
191
192void test_scalar(ITHACAparameters* para, Foam::fvMesh& mesh,
193 Foam::Time& runTime)
194{
195 int n_modes = para->ITHACAdict->lookupOrDefault<int>("Modes", 15);
196 int n_nodes = para->ITHACAdict->lookupOrDefault<int>("Nodes", 15);
197 simpleControl simple(mesh);
198 word methodName = para->ITHACAdict->lookupOrDefault<word>("HyperReduction",
199 "GappyDEIM");
200#include "createFields.H"
201 // List of volScalarField where the snapshots are stored
202 PtrList<volScalarField> Sp;
203 // Non linear function
204 volScalarField S
205 (
206 IOobject
207 (
208 "S",
209 runTime.timeName(),
210 mesh,
211 IOobject::NO_READ,
212 IOobject::AUTO_WRITE
213 ),
214 T.mesh(),
215 dimensionedScalar("zero", dimensionSet(0, 0, -1, 1, 0, 0, 0), 0)
216 );
217 // Parameters used to train the non-linear function
218 Eigen::MatrixXd pars = ITHACAutilities::rand(100, 2, -0.5, 0.5);
219 //cnpy::load(pars, "trainingPars.npy");
220
221 // Perform the offline phase
222 for (int i = 0; i < 100; i++)
223 {
224 HyperReduction_function::evaluate_expression(S, pars.row(i));
225 Sp.append((S).clone());
226 ITHACAstream::exportSolution(S, name(i + 1), "./ITHACAoutput/Offline/");
227 }
228
229 // Define new online parameters
230 Eigen::MatrixXd parTest = ITHACAutilities::rand(100, 2, -0.5, 0.5);
231 //cnpy::load(parTest, "testingPars.npy");
232 // Create HyperReduction object with given number of basis functions
233 Eigen::VectorXi initSeeds;
234 // initSeeds = cnpy::load(initSeeds, "./mp.npy");//load mp to test initialSeeds
235
236 if (methodName == "GappyDEIM")
237 {
238 HyperReduction_function c(n_modes, n_nodes, initSeeds, "Gaussian_function", Sp);
239 Eigen::MatrixXd snapshotsModes;
240 Eigen::VectorXd normalizingWeights;
241 c.getModesSVD(c.snapshotsListTuple, snapshotsModes, normalizingWeights);
242 c.offlineGappyDEIM(snapshotsModes, normalizingWeights);
243 // Generate the submeshes with the depth of the layer
244 unsigned int layers{2};
245 c.generateSubmesh(layers, mesh);
246 auto sfield = c.interpolateField<volScalarField>(Sp[0]);
247 PtrList<volScalarField> testTrueFields;
248 PtrList<volScalarField> testRecFields;
249
250 for (unsigned int idTest = 0; idTest < parTest.rows(); idTest++)
251 {
252 // Online evaluation of the non linear function
253 Eigen::VectorXd aprfield = c.renormalizedBasisMatrix * c.onlineCoeffs(sfield(),
254 parTest.row(idTest));
255 // Transform to an OpenFOAM field and export
256 volScalarField S2("S_online", Foam2Eigen::Eigen2field(S, aprfield));
257 // Evaluate the full order function and export it
258 HyperReduction_function::evaluate_expression(S, parTest.row(idTest));
259 testRecFields.append(S2.clone());
260 testTrueFields.append(S.clone());
261 ITHACAstream::exportSolution(S2, name(idTest), "./ITHACAoutput/Online/");
262 ITHACAstream::exportSolution(S, name(idTest), "./ITHACAoutput/Online/");
263 }
264
265 // Compute the L2 error and print it
266 auto err = ITHACAutilities::errorL2Rel(testRecFields, testTrueFields);
267 Info << "GappyDEIM max relative error: " << err.maxCoeff() << endl;
268 Info << "GappyDEIM mean relative error: " << err.array().mean() << endl;
269 }
270 else if (methodName == "ECP")
271 {
272 HyperReduction_function c(n_modes, n_nodes, initSeeds, "Gaussian_function", Sp);
273 bool saveOFModes{true};
274 Eigen::MatrixXd snapshotsModes;
275 Eigen::VectorXd normalizingWeights;
276 c.getModesSVD(c.snapshotsListTuple, snapshotsModes, normalizingWeights,
277 saveOFModes);
278 c.offlineECP(snapshotsModes, normalizingWeights);
279 // Generate the submeshes with the depth of the layer
280 unsigned int layers{2};
281 c.generateSubmesh(layers, mesh);
282 auto sfield = c.interpolateField<volScalarField>(Sp[0]);
283 Eigen::VectorXd results(parTest.rows());
284
285 for (unsigned int idTest = 0; idTest < parTest.rows(); idTest++)
286 {
287 // Online evaluation of the non linear function
288 Eigen::VectorXd f = c.evaluate_expression(sfield(), parTest.row(idTest),
289 c.localNodePoints);
290 Eigen::VectorXd ff = Foam2Eigen::field2Eigen(c.evaluate_expression(Sp[0],
291 parTest.row(idTest)));
292 double trueIntegral = (normalizingWeights.cwiseInverse().asDiagonal() *
293 ff).array().sum();
294 double testIntegral = (c.wPU * f).array().sum();
295 // Info << "Integral: " << trueIntegral << endl;
296 // Info << "Reconstructed Integral: " << testIntegral << endl;
297 results(idTest) = trueIntegral - testIntegral;
298 }
299
300 Info << "ECP max error: " << results.cwiseAbs().maxCoeff() << endl;
301 Info << "ECP mean error: " << results.cwiseAbs().mean() << endl;
302 }
303}
304
305void test_vector(ITHACAparameters* para, Foam::fvMesh& mesh,
306 Foam::Time& runTime)
307{
308 int n_modes = para->ITHACAdict->lookupOrDefault<int>("Modes", 15);
309 int n_nodes = para->ITHACAdict->lookupOrDefault<int>("Nodes", 15);
310 simpleControl simple(mesh);
311 word methodName = para->ITHACAdict->lookupOrDefault<word>("HyperReduction",
312 "GappyDEIM");
313#include "createFields.H"
314 // List of volVectorField where the snapshots are stored
315 PtrList<volVectorField> Sp;
316 // Non linear function
317 volVectorField S
318 (
319 IOobject
320 (
321 "S",
322 runTime.timeName(),
323 mesh,
324 IOobject::NO_READ,
325 IOobject::AUTO_WRITE
326 ),
327 T.mesh(),
328 dimensionedVector("zero", dimensionSet(0, 0, -1, 1, 0, 0, 0), vector(0, 0, 0))
329 );
330 // Parameters used to train the non-linear function
331 Eigen::MatrixXd pars = ITHACAutilities::rand(100, 2, -0.5, 0.5);
332 // cnpy::load(pars, "trainingPars.npy");
333
334 // Perform the offline phase
335 for (int i = 0; i < 100; i++)
336 {
337 HyperReduction_vectorFunction::evaluate_expression(S, pars.row(i));
338 Sp.append((S).clone());
339 ITHACAstream::exportSolution(S, name(i + 1), "./ITHACAoutput/Offline/");
340 }
341
342 // Define new online parameters
343 Eigen::MatrixXd parTest = ITHACAutilities::rand(100, 2, -0.5, 0.5);
344 // cnpy::load(parTest, "testingPars.npy");
345 // Create HyperReduction object with given number of basis functions
346 Eigen::VectorXi initSeeds;
347 // initSeeds = cnpy::load(initSeeds, "./mp.npy");//load mp to test initialSeeds
348
349 if (methodName == "GappyDEIM")
350 {
351 HyperReduction_vectorFunction c(n_modes, n_nodes, initSeeds,
352 "Gaussian_function", Sp);
353 Eigen::MatrixXd snapshotsModes;
354 Eigen::VectorXd normalizingWeights;
355 c.getModesSVD(c.snapshotsListTuple, snapshotsModes, normalizingWeights);
356 c.offlineGappyDEIM(snapshotsModes, normalizingWeights);
357 // Generate the submeshes with the depth of the layer
358 unsigned int layers{2};
359 c.generateSubmesh(layers, mesh);
360 auto sfield = c.interpolateField<volVectorField>(Sp[0]);
361 PtrList<volVectorField> testTrueFields;
362 PtrList<volVectorField> testRecFields;
363
364 for (unsigned int idTest = 0; idTest < parTest.rows(); idTest++)
365 {
366 // Online evaluation of the non linear function
367 Eigen::VectorXd aprfield = c.renormalizedBasisMatrix * c.onlineCoeffs(sfield(),
368 parTest.row(idTest));
369 // Transform to an OpenFOAM field and export
370 volVectorField S2("S_online", Foam2Eigen::Eigen2field(S, aprfield));
371 // Evaluate the full order function and export it
372 HyperReduction_vectorFunction::evaluate_expression(S, parTest.row(idTest));
373 testRecFields.append(S2.clone());
374 testTrueFields.append(S.clone());
375 ITHACAstream::exportSolution(S2, name(idTest), "./ITHACAoutput/Online/");
376 ITHACAstream::exportSolution(S, name(idTest), "./ITHACAoutput/Online/");
377 }
378
379 // Compute the L2 error and print it
380 auto err = ITHACAutilities::errorL2Rel(testRecFields, testTrueFields);
381 Info << "GappyDEIM max relative error: " << err.maxCoeff() << endl;
382 Info << "GappyDEIM mean relative error: " << err.array().mean() << endl;
383 }
384 else if (methodName == "ECP")
385 {
386 HyperReduction_vectorFunction c(n_modes, n_nodes, initSeeds,
387 "Gaussian_function", Sp);
388 bool saveOFModes{true};
389 Eigen::MatrixXd snapshotsModes;
390 Eigen::VectorXd normalizingWeights;
391 c.getModesSVD(c.snapshotsListTuple, snapshotsModes, normalizingWeights,
392 saveOFModes);
393 c.offlineECP(snapshotsModes, normalizingWeights);
394 // Generate the submeshes with the depth of the layer
395 unsigned int layers{2};
396 c.generateSubmesh(layers, mesh);
397 auto sfield = c.interpolateField<volVectorField>(Sp[0]);
398 Eigen::VectorXd results(parTest.rows());
399
400 for (unsigned int idTest = 0; idTest < parTest.rows(); idTest++)
401 {
402 // Online evaluation of the non linear function
403 Eigen::VectorXd f = c.evaluate_expression(sfield(), parTest.row(idTest),
404 c.localNodePoints);
405 auto wholeField = c.evaluate_expression(Sp[0], parTest.row(idTest));
406 Eigen::VectorXd ff = Foam2Eigen::field2Eigen(wholeField);
407 double trueIntegral = (normalizingWeights.cwiseInverse().asDiagonal() *
408 ff).array().sum();
409 double testIntegral = (c.wPU * f).array().sum();
410 // Info << "Integral: " << trueIntegral << endl;
411 // Info << "Reconstructed Integral: " << testIntegral << endl;
412 results(idTest) = trueIntegral - testIntegral;
413 }
414
415 Info << "ECP max error: " << results.cwiseAbs().maxCoeff() << endl;
416 Info << "ECP mean error: " << results.cwiseAbs().mean() << endl;
417 }
418}
419
420void test_vector_scalar(ITHACAparameters* para, Foam::fvMesh& mesh,
421 Foam::Time& runTime)
422{
423 int n_modes = para->ITHACAdict->lookupOrDefault<int>("Modes", 15);
424 int n_nodes = para->ITHACAdict->lookupOrDefault<int>("Nodes", 15);
425 simpleControl simple(mesh);
426 word methodName = para->ITHACAdict->lookupOrDefault<word>("HyperReduction",
427 "GappyDEIM");
428#include "createFields.H"
429 // List of volVectorField where the snapshots are stored
430 PtrList<volVectorField> Vp;
431 PtrList<volScalarField> Sp;
432 // Non linear function
433 volVectorField V
434 (
435 IOobject
436 (
437 "V",
438 runTime.timeName(),
439 mesh,
440 IOobject::NO_READ,
441 IOobject::AUTO_WRITE
442 ),
443 T.mesh(),
444 dimensionedVector("zero", dimensionSet(0, 0, -1, 1, 0, 0, 0), vector(0, 0, 0))
445 );
446 // Non linear function
447 volScalarField S
448 (
449 IOobject
450 (
451 "S",
452 runTime.timeName(),
453 mesh,
454 IOobject::NO_READ,
455 IOobject::AUTO_WRITE
456 ),
457 T.mesh(),
458 dimensionedScalar("zero", dimensionSet(0, 0, -1, 1, 0, 0, 0), 0)
459 );
460 // Parameters used to train the non-linear function
461 Eigen::MatrixXd pars = ITHACAutilities::rand(100, 2, -0.5, 0.5);
462 // cnpy::load(pars, "trainingPars.npy");
463
464 // Perform the offline phase
465 for (int i = 0; i < 100; i++)
466 {
467 HyperReduction_vectorScalarFunction::evaluate_expression(V, S, pars.row(i));
468 Vp.append((V).clone());
469 Sp.append((S).clone());
470 ITHACAstream::exportSolution(V, name(i + 1), "./ITHACAoutput/Offline/");
471 ITHACAstream::exportSolution(S, name(i + 1), "./ITHACAoutput/Offline/");
472 }
473
474 // Define new online parameters
475 Eigen::MatrixXd parTest = ITHACAutilities::rand(100, 2, -0.5, 0.5);
476 // cnpy::load(parTest, "testingPars.npy");
477 // Create HyperReduction object with given number of basis functions
478 Eigen::VectorXi initSeeds;
479 // initSeeds = cnpy::load(initSeeds, "./mp.npy");//load mp to test initialSeeds
480
481 if (methodName == "GappyDEIM")
482 {
483 HyperReduction_vectorScalarFunction c(n_modes, n_nodes, initSeeds,
484 "Gaussian_function", Vp, Sp);
485 Eigen::MatrixXd snapshotsModes;
486 Eigen::VectorXd normalizingWeights;
487 c.getModesSVD(c.snapshotsListTuple, snapshotsModes, normalizingWeights);
488 c.offlineGappyDEIM(snapshotsModes, normalizingWeights);
489 // Generate the submeshes with the depth of the layer
490 unsigned int layers{2};
491 c.generateSubmesh(layers, mesh);
492 auto vfield = c.interpolateField<volVectorField>(Vp[0]);
493 auto sfield = c.interpolateField<volScalarField>(Sp[0]);
494 PtrList<volVectorField> testTrueFieldsV;
495 PtrList<volVectorField> testRecFieldsV;
496 PtrList<volScalarField> testTrueFieldsS;
497 PtrList<volScalarField> testRecFieldsS;
498
499 for (unsigned int idTest = 0; idTest < parTest.rows(); idTest++)
500 {
501 // Online evaluation of the non linear function
502 auto theta = c.onlineCoeffs(sfield(), parTest.row(idTest));
503 Eigen::VectorXd aprfield = c.renormalizedBasisMatrix * theta;
504 Eigen::VectorXd recvec = aprfield.head(3 * S.size());
505 Eigen::VectorXd recsca = aprfield.tail(S.size());
506 // Transform to an OpenFOAM field and export
507 volVectorField V2("V_online", Foam2Eigen::Eigen2field(V, recvec));
508 volScalarField S2("S_online", Foam2Eigen::Eigen2field(S, recsca));
509 // Evaluate the full order function and export it
510 HyperReduction_vectorScalarFunction::evaluate_expression(V, S,
511 parTest.row(idTest));
512 testRecFieldsV.append(V2.clone());
513 testTrueFieldsV.append(V.clone());
514 ITHACAstream::exportSolution(V2, name(idTest), "./ITHACAoutput/Online/");
515 ITHACAstream::exportSolution(V, name(idTest), "./ITHACAoutput/Online/");
516 testRecFieldsS.append(S2.clone());
517 testTrueFieldsS.append(S.clone());
518 ITHACAstream::exportSolution(S2, name(idTest), "./ITHACAoutput/Online/");
519 ITHACAstream::exportSolution(S, name(idTest), "./ITHACAoutput/Online/");
520 }
521
522 // Compute the L2 error and print it
523 auto errV = ITHACAutilities::errorL2Rel(testRecFieldsV, testTrueFieldsV);
524 Info << "GappyDEIM vector field max relative error: " << errV.maxCoeff() <<
525 endl;
526 Info << "GappyDEIM vector field mean relative error: " << errV.array().mean() <<
527 endl;
528 auto errS = ITHACAutilities::errorL2Rel(testRecFieldsS, testTrueFieldsS);
529 Info << "GappyDEIM scalar max relative error: " << errS.maxCoeff() << endl;
530 Info << "GappyDEIM scalar mean relative error: " << errS.array().mean() << endl;
531 }
532 else if (methodName == "ECP")
533 {
534 HyperReduction_vectorScalarFunction c(n_modes, n_nodes, initSeeds,
535 "Gaussian_function", Vp, Sp);
536 bool saveOFModes{true};
537 Eigen::MatrixXd snapshotsModes;
538 Eigen::VectorXd normalizingWeights;
539 c.getModesSVD(c.snapshotsListTuple, snapshotsModes, normalizingWeights,
540 saveOFModes);
541 c.offlineECP(snapshotsModes, normalizingWeights);
542 // Generate the submeshes with the depth of the layer
543 unsigned int layers{2};
544 c.generateSubmesh(layers, mesh);
545 auto vfield = c.interpolateField<volVectorField>(Vp[0]);
546 auto sfield = c.interpolateField<volScalarField>(Sp[0]);
547 Eigen::VectorXd results(parTest.rows());
548
549 for (unsigned int idTest = 0; idTest < parTest.rows(); idTest++)
550 {
551 // Online evaluation of the non linear function
552 Eigen::VectorXd f = c.evaluate_expression(sfield(), parTest.row(idTest),
553 c.localNodePoints);
554 HyperReduction_vectorScalarFunction::evaluate_expression(V, S,
555 parTest.row(idTest));
556 Eigen::VectorXd ffV = Foam2Eigen::field2Eigen(V);
557 Eigen::VectorXd ffS = Foam2Eigen::field2Eigen(S);
558 Eigen::VectorXd ff(ffV.rows() + ffS.rows());
559 ff.head(ffV.rows()) = ffV;
560 ff.tail(ffS.rows()) = ffS;
561 double trueIntegral = (normalizingWeights.cwiseInverse().asDiagonal() *
562 ff).array().sum();
563 double testIntegral = (c.wPU * f).array().sum();
564 // Info << "Integral: " << trueIntegral << endl;
565 // Info << "Reconstructed Integral: " << testIntegral << endl;
566 results(idTest) = trueIntegral - testIntegral;
567 }
568
569 Info << "ECP max error: " << results.cwiseAbs().maxCoeff() << endl;
570 Info << "ECP mean error: " << results.cwiseAbs().mean() << endl;
571 }
572}
573
574int main(int argc, char* argv[])
575{
576 // load stage to perform
577 argList::addOption("test", "scalar", "Perform scalar test");
578 argList::addOption("test", "vector", "Perform vector test");
579 argList::addOption("test", "vs", "Perform (vector, scalar) test");
580 // add options for parallel run
581 HashTable<string> validParOptions;
582 validParOptions.set
583 (
584 "test",
585 "scalar"
586 );
587 validParOptions.set
588 (
589 "test",
590 "vector"
591 );
592 validParOptions.set
593 (
594 "test",
595 "vs"
596 );
597 // Read parameters from ITHACAdict file
598#include "setRootCase.H"
599 Foam::Time runTime(Foam::Time::controlDictName, args);
600 Foam::fvMesh mesh
601 (
602 Foam::IOobject
603 (
604 Foam::fvMesh::defaultRegion,
605 runTime.timeName(),
606 runTime,
607 Foam::IOobject::MUST_READ
608 )
609 );
610 ITHACAparameters* para = ITHACAparameters::getInstance(mesh, runTime);
611
612 if (args.get("test").match("scalar"))
613 {
614 Info << "Init scalar testcase\n";
615 test_scalar(para, mesh, runTime);
616 }
617 else if (args.get("test").match("vector"))
618 {
619 Info << "Init vector testcase\n";
620 test_vector(para, mesh, runTime);
621 }
622 else if (args.get("test").match("vs"))
623 {
624 Info << "Init vector scalar testcase\n";
625 test_vector_scalar(para, mesh, runTime);
626 }
627 else
628 {
629 Info << "Pass '-test scalar', '-test vector', '-test vs'" << endl;
630 }
631
632 return 0;
633}
Foam2Eigen.H
Header file of the Foam2Eigen class.
ITHACAPOD.H
Header file of the ITHACAPOD class.
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.
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.
HyperReduction_function
Definition 24HyperReduction.C:44
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_vectorFunction
Definition 24HyperReduction.C:87
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_vectorScalarFunction
Definition 24HyperReduction.C:138
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::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
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
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::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:337
ITHACAutilities::rand
Eigen::MatrixXd rand(label rows, label cols, double min, double max)
Generates random matrix with random values in an interval.
Definition ITHACAutilities.C:44
Generated by   doxygen 1.16.1