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