Loading...
Searching...
No Matches
ScalarTransport.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 a heat transfer Reduction Problem
26SourceFiles
27 02thermalBlock.C
28\*---------------------------------------------------------------------------*/
29
30//#include "laplacianProblem.H"
31#include "ScalarTransport.H"
32
33//#include "fvMeshSubsetter.H" // Not fvMeshSubset (need two-step subsetting)
37
38
39// * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * * * //
40// Constructor
41ScalarTransport::ScalarTransport() {}
42
43ScalarTransport::ScalarTransport(int argc, char* argv[])
44 :
45 UnsteadyProblem()
46{
47 _args = autoPtr<argList>
48 (
49 new argList(argc, argv)
50 );
51
52 if (!_args->checkRootCase())
53 {
54 Foam::FatalError.exit();
55 }
56
57 argList& args = _args();
58#include "createTime.H"
59#include "createMesh.H"
60 _simple = autoPtr<simpleControl>
61 (
62 new simpleControl
63 (
64 mesh
65 )
66 );
67 simpleControl& simple = _simple();
68#include "createFields.H"
69#include "createFvOptions.H"
70 ITHACAdict = new IOdictionary
71 (
72 IOobject
73 (
74 "ITHACAdict",
75 runTime.system(),
76 mesh,
77 IOobject::MUST_READ,
78 IOobject::NO_WRITE
79 )
80 );
81 para = ITHACAparameters::getInstance(mesh, runTime);
82 offline = ITHACAutilities::check_off();
83 podex = ITHACAutilities::check_pod();
84 setTimes(runTime);
85}
86
87// void offlineSolve(word folder = "./ITHACAoutput/Offline/")
88// {
89// if (offline)
90// {
91// ITHACAstream::read_fields(Tfield, T, "./ITHACAoutput/Offline/");
92// //ITHACAstream::readMiddleFields(Ufield, U, "./ITHACAoutput/Offline/");
93
94// }
95// else
96// {
97// truthSolve(folder);
98// }
99// }
100
101// void ELI(int NmodesT)
102// {
103// fvMesh& mesh = const_cast<fvMesh&>(T.mesh());
104// EliObject = autoPtr<EmpLagranInter> (new EmpLagranInter(Tfield, NmodesT,"Tsol", T.name()));
105// EliObject->generateSubmesh(1, T.mesh(), T);
106// // EliObject->generateSubmesh(2, U.mesh(), U);
107// }
108void ScalarTransport::truthSolve(word folder)
109{
110 Time& runTime = _runTime();
111 fvMesh& mesh = _mesh();
112 volScalarField& T = _T();
113 surfaceScalarField& phi = _phi();
114 fv::options& fvOptions = _fvOptions();
115 simpleControl& simple = _simple();
116 dimensionedScalar& nu = _nu();
117 counter = 1;
118 //ITHACAstream::exportSolution(U, name(counter), folder);
119 //counter++;
120 nextWrite = startTime;
121 nextWrite += writeEvery;
122
123 while (_simple().loop())
124 {
125 Info << "Time = " << _runTime().timeName() << nl << endl;
126
127 while (simple.correctNonOrthogonal())
128 {
129 fvScalarMatrix TEqn
130 (
131 fvm::ddt(T)
132 + fvm::div(phi, T)
133 - fvm::laplacian(nu, T)
134 );
135 TEqn.relax();
136 fvOptions.constrain(TEqn);
137 TEqn.solve();
138 fvOptions.correct(T);
139 }
140
141 //phi = linearInterpolate(U) & mesh.Sf();
142
143 if (checkWrite(runTime))
144 {
145 ITHACAstream::exportSolution(T, name(counter), folder );
146 counter++;
147 Tfield.append(T.clone());
148 nextWrite += writeEvery;
149 }
150 }
151}
152
153void ScalarTransport::restart()
154{
155 _T.clear();
156 _U.clear();
157 _phi.clear();
158 _fvOptions.clear();
159 _nu.clear();
160 _transportProperties.clear();
161 argList& args = _args();
162 Time& runTime = _runTime();
163 runTime.setTime(0, 1);
164 Foam::fvMesh& mesh = _mesh();
165#include "createFields.H"
166#include "createFvOptions.H"
167}
168
169// class EliBurgers: public reductionProblem
170// {
171// public:
172
173// //// Data members
174// tutorial23* problem;
175// scalar nextWrite;
176// scalar startTime;
177// scalar writeEvery;
178// PtrList<volVectorField> Ufield;
179// //volVectorModes Lagrangian;
180// List<Eigen::VectorXd> OnlineCoeffs;
181// EliBurgers(){}
182// explicit EliBurgers(tutorial23& FoamProblem) : problem(&FoamProblem)
183// {
184// nextWrite = 0;
185// startTime = problem->startTime;
186// writeEvery = problem->writeEvery;
187
188// Info << "######## reduced Constructor calling ##########"<< endl;
189// }
190
191// ~EliBurgers(){}
192
193
194// void OnlineBurgers(int NmodesUproj,word folder = "./ITHACAoutput/Online/")
195// {
196// // Time& runTime = problem->_runTime();
197// // //fvMesh& mesh = problem->_mesh();
198// //volVectorField& U = problem->_U();
199// // surfaceScalarField& phi = problem->_phi();
200// // //simpleControl& simple = problem->_simple();
201// //dimensionedScalar& nu = problem->_nu();
202// /// simpleControl
203// autoPtr<simpleControl> _simple;
204// // /// fvOptions
205// autoPtr<fv::options> _fvOptions;
206
207// auto submesh = problem->EliObject->submesh;
208// submesh->setCellSubset(problem->EliObject->uniqueMagicPoints() );
209// submesh->subMesh().fvSchemes::readOpt() = problem->_mesh().fvSchemes::readOpt();
210// submesh->subMesh().fvSolution::readOpt() = problem->_mesh().fvSolution::readOpt();
211
212// submesh->subMesh().fvSchemes::read();
213// submesh->subMesh().fvSolution::read();
214// //submesh().subMesh().write();
215// _simple = autoPtr<simpleControl>
216// (
217// new simpleControl
218// (
219// submesh().subMesh()
220// )
221// );
222
223// _fvOptions = autoPtr<fv::options>(new fv::options(submesh().subMesh()));
224// simpleControl& simple = _simple();
225// fv::options& fvOptions = _fvOptions();
226
227// autoPtr<volVectorField> _S;
228// _S = autoPtr<volVectorField>
229// (
230// new volVectorField
231// (
232// IOobject
233// (
234// "S",
235// problem->_runTime().timeName(),
236// submesh().subMesh(),
237// IOobject::NO_READ,
238// IOobject::NO_WRITE
239// ),
240// submesh().subMesh()
241// )
242// );
243// volVectorField& S = _S();
244// //Info << "############################################" << endl;
245// autoPtr<surfaceScalarField> _sphi;
246// _sphi = autoPtr<surfaceScalarField>
247// (
248// new surfaceScalarField
249// (
250// IOobject
251// (
252// "sphi",
253// problem->_runTime().timeName(),
254// submesh().subMesh(),
255// IOobject::READ_IF_PRESENT,
256// IOobject::NO_WRITE
257// ),
258// linearInterpolate(S) & submesh().subMesh().Sf()
259// )
260// );
261// surfaceScalarField& sphi = _sphi();
262// //Info << submesh->subMesh().schemesDict().size() << endl;
263// // Info << submesh->subMesh().solutionDict() << endl;
264// // exit(0);
265// counter = 1;
266// // int dim = 3;
267// // label nCells = problem->_mesh().nCells();
268// // //exit(0);
269// // /// Mask field to submesh
270// // Eigen::SparseMatrix<double> field2submesh;
271// // label M = problem->EliObject->uniqueMagicPoints().size();
272// // field2submesh.resize(nCells*dim, M*dim);
273
274// // for (unsigned int ith_subCell{0} ; ith_subCell < M; ith_subCell++)
275// // {
276// // for (unsigned int i= 0; i < dim; i++)
277// // {
278// // field2submesh.insert(problem->EliObject->uniqueMagicPoints()[ith_subCell] + nCells*i, ith_subCell + i*M) = 1;
279// // }
280// // }
281// // field2submesh.makeCompressed();
282// // /// Important matrices
283// // Eigen::MatrixXd X = problem->EliObject->U;
284// // Eigen::MatrixXd V = (X.transpose()*field2submesh*field2submesh.transpose()*X).inverse()*X.transpose()*field2submesh; // Ok
285// // /// Define the full matrix to the full field
286// // Eigen::MatrixXd Y = X*V;
287// ////
288// nextWrite = startTime;
289// nextWrite += writeEvery;
290
291// Eigen::SparseMatrix<double> Se;
292
293// Eigen::VectorXd se;
294// Eigen::SimplicialLLT<Eigen::SparseMatrix<double> > solver; //0k
295// // performs a Cholesky factorization of A
296// //Eigen::SimplicialCholesky<Eigen::SparseMatrix<double>> chol(Se);
297// //Eigen::SparseQR<Eigen::SparseMatrix<double> > solver;
298// //Eigen::SimplicialLDLT <Eigen::SparseMatrix<double> > solver;
299// //Eigen::SparseLU<Eigen::SparseMatrix<double>> solver;
300// // Eigen::VectorXd a;
301// //int r = problem->EliObject->MatrixOnline.cols();
302// //a.resize(r,1);
303// while (simple.loop() )
304// {
305// Info << "Time = " << problem->_runTime().timeName() << endl;
306
307// while (simple.correctNonOrthogonal())
308// {
309
310// fvVectorMatrix SEqn
311// (
312// fvm::ddt(S)
313// + fvm::div(sphi,S)
314// - fvm::laplacian(problem->_nu(),S)
315// );
316// //SEqn.relax();
317// Foam2Eigen::fvMatrix2Eigen(SEqn, Se, se);
318// /// Hyper-reduced system
319// /// TODO: make sure to include the values of the volumes.
320// Eigen::VectorXd a = solver.compute(Se).solve(se);
321// //Info << "a.rows() = " << a.rows() << endl;
322// //Info << "a.cols() = " << a.cols() << endl;
323// //Eigen::VectorXd z = Y*a;
324// //exit(0);
325// /// Solution in eigen format
326// S = Foam2Eigen::Eigen2field(_S(), a);
327// }
328// sphi = linearInterpolate(S) & submesh->subMesh().Sf();
329
330// if (checkWrite(problem->_runTime()))
331// {
332// //Info << "########## line 190" << endl;
333// //ITHACAstream::exportSolution(S, name(counter), folder);
334// counter++;
335// Ufield.append(S.clone());
336// //OnlineCoeffs.append(a);
337// nextWrite += writeEvery;
338// }
339// }
340
341// }
342
343// bool checkWrite(Time& timeObject)
344// {
345// scalar diffnow = mag(nextWrite - atof(timeObject.timeName().c_str()));
346// scalar diffnext = mag(nextWrite - atof(timeObject.timeName().c_str()) -
347// timeObject.deltaTValue());
348
349// if ( diffnow < diffnext)
350// {
351// return true;
352// }
353// else
354// {
355// return false;
356// }
357// }
358
359// };
360
361
362
363
364// int main(int argc, char* argv[])
365// {
366// // Create the train object of the tutorial02 type
367// tutorial24 train(argc, argv);
368// //tutorial23 test(argc, argv);
369// std::clock_t startOff;
370// double durationOff;
371
372// // Read some parameters from file
373// ITHACAparameters* para = ITHACAparameters::getInstance(train._mesh(),train._runTime());
374// int NmodesUout = readInt(para->ITHACAdict->lookup("NmodesUout"));
375// int NmodesUproj = readInt(para->ITHACAdict->lookup("NmodesUproj"));
376
377// startOff= std::clock();
378// train.offlineSolve();
379// durationOff = (std::clock() - startOff);
380// Info << "The Offline phase duration is equal to" << durationOff << endl;
381// // ITHACAPOD::getModes(train.Ufield, train.Umodes, train._U().name(),
382// // train.podex, 0, 0,NmodesUout);
383
384// // /// Construct the DEIM indices.
385// // train.ELI(NmodesUproj);
386
387// // int dim = 3;
388// // label nCells = train._mesh().nCells();
389// // //exit(0);
390// // /// Mask field to submesh
391// // Eigen::SparseMatrix<double> field2submesh;
392// // label M = train.EliObject->uniqueMagicPoints().size();
393// // field2submesh.resize(nCells*dim, M*dim);
394
395// // for (unsigned int ith_subCell{0} ; ith_subCell < M; ith_subCell++)
396// // {
397// // for (unsigned int i= 0; i < dim; i++)
398// // {
399// // field2submesh.insert(train.EliObject->uniqueMagicPoints()[ith_subCell] + nCells*i, ith_subCell + i*M) = 1;
400// // }
401// // }
402// // field2submesh.makeCompressed();
403
404// // /// Important matrices
405// // Eigen::MatrixXd X = train.EliObject->U;
406// // Eigen::MatrixXd V = (X.transpose()*field2submesh*field2submesh.transpose()*X).inverse()*X.transpose()*field2submesh; // Ok
407// // /// Define the full matrix to the full field
408// // Eigen::MatrixXd Y = X*V;
409// // //Info << "Y.rows() = " << Y.rows() << endl;
410// // //Info << "Y.cols() = " << Y.cols() << endl;
411// // train.EliObject->MatrixOnline = Y;
412
413// // // PtrList<volVectorField> subfields;
414
415// // // for (int i = 0; i < train.Ufield.size(); ++i)
416// // // {
417
418// // // auto subfld = train.EliObject->submesh().interpolate(train.Ufield[i] ).ref();
419// // // subfields.append(subfld);
420
421// // // }
422// // // ITHACAstream::exportFields(subfields, "./ITHACAoutput/SubFields", subfields[0].name() );
423
424// // //exit(0);
425// // std::clock_t startOn;
426// // double durationOn;
427// // startOn = std::clock();
428
429// // train.restart();
430// // EliBurgers reduced(train);
431
432// // reduced.OnlineBurgers(NmodesUproj);
433
434// // durationOn = (std::clock() - startOn);
435// // Info << "The Online phase duration is equal to " << durationOn << endl;
436// // //ITHACAstream::exportFields(reduced.Ufield, "./ITHACAoutput/Online", reduced.Ufield[0].name() );
437// // //PtrList<volVectorField> errflds;
438// // //exit(0);
439// // /// Mapping to the whole fields
440// // for (int i = 0; i < reduced.Ufield.size(); ++i)
441// // {
442// // Eigen::VectorXd c = Foam2Eigen::field2Eigen(reduced.Ufield[i]);
443// // Eigen::VectorXd z = Y*c; //reduced.OnlineCoeffs[i];
444// // auto fld = Foam2Eigen::Eigen2field(train._U(), z);
445// // fld.correctBoundaryConditions();
446// // ITHACAstream::exportSolution(fld, name(i+1), "./ITHACAoutput/Online");
447// // }
448// //ITHACAstream::exportFields(reduced.Ufield, "./ITHACAoutput/Online", "S");
449// exit(0);
450// }
451
452// // ########################################################
453// label n = train.EliObject->uniqueMagicPoints().size();
454// // int size = train.EliObject->magicPoints().size();
455// Info << train.EliObject->magicPoints() << endl;
456// Info << train.EliObject->uniqueMagicPoints() << endl;
457// // //Info << train.EliObject->uniqueMagicPoints()[5] << endl;
458
459// // for (int i = 0; i <n ; ++i)
460// // {
461// // //Info << "mesh.cells() = " << train._mesh().cells()[train.EliObject->uniqueMagicPoints()[5]] << endl;
462
463// // auto cell = train._mesh().cells()[train.EliObject->uniqueMagicPoints()[i]];
464// // Info << cell << endl;
465// // forAll (cell, faceI)
466// // {
467// // Info << cell[faceI] << endl;
468// // }
469
470// // }
471// // // ##############################################################
472
473
474
475//for(label i = 0; i < M; i++)
476// {
477
478// CoeffMat[i][i] = UEqn.diag()[problem->EliObject->magicPoints()[i]];
479// CoeffMat.source()[i] = UEqn.source()[problem->EliObject->magicPoints()[i]];
480// }
481// Info << "Time2 = " << problem->_runTime().timeName() << nl << endl;
482// // Assigning off-diagonal coefficients
483
484// for (label k = 0; k < M; k++)
485// {
486// //Info << "mesh.cells() = " << train._mesh().cells()[train.EliObject->uniqueMagicPoints()[5]] << endl;
487// auto cell = mesh.cells()[problem->EliObject->magicPoints()[k]];
488// //Info << cell << endl;
489// forAll (cell, faceI)
490// {
491// label l = UEqn.lduAddr().lowerAddr()[faceI];
492// label u = UEqn.lduAddr().upperAddr()[faceI];
493// CoeffMat[l][u] = UEqn.upper()[faceI];
494// CoeffMat[u][l] = UEqn.lower()[faceI];
495// //Info << cell[faceI] << endl;
496// }
497
498// }
499
500// auto x = CoeffMat.LUsolve();
501// Info << x << endl;
502// //volVectorField x(xx);
503
504// Info << problem->EliObject->MatrixOnline.cols() << endl;
505// Info << problem->EliObject->MatrixOnline.rows() << endl;
506
507
508
509// Eigen::VectorXd y = Foam2Eigen::field2Eigen(x);
510
511// Info << y.cols() << endl;
512// Info << y.rows() << endl;
513
514// exit(0);
515
516
517// Eigen::VectorXd z = problem->EliObject->MatrixOnline*y;
518
519//U = Foam2Eigen::Eigen2field(problem->_U(),z, true);
520//Info<< CoeffMat.solve() << endl;
521//Info<< x << endl;
522
523// Assigning contribution from BC
524// forAll(U.boundaryField(), patchI)
525// {
526// const fvPatch &pp = U.boundaryField()[patchI].patch();
527// forAll(pp, faceI)
528// {
529// label cellI = pp.faceCells()[faceI];
530// //CoeffMat[cellI][cellI] += UEqn.internalCoeffs()[patchI][faceI];
531// //CoeffMat.source()[cellI] +=UEqn.boundaryCoeffs()[patchI][faceI];
532// }
533// }
534// Info << CoeffMat << endl;
535// Foam2Eigen::fvMatrix2Eigen(UEqn, Ae, be);
536// auto Ar = P.transpose() * Ae * P;
537// auto br = P.transpose() * be;
538// auto Ap = Q.transpose()* Ar *Q;
539// auto bp = Q.transpose() * br;
540
541//Info << "P = " << P << endl;
542
543//Info << "Ar = " << Q.transpose()* Ar *Q << endl;
544//Info << "br = " << Q.transpose() * br << endl;
545//auto a = Ap.colPivHouseholderQr().solve(bp);
ScalarTransport.H
Header file of the Burgers class.
ITHACAparameters::getInstance
static ITHACAparameters * getInstance()
Gets an instance of ITHACAparameters, to be used if the instance is already existing.
Definition ITHACAparameters.C:56
ScalarTransport::_phi
autoPtr< surfaceScalarField > _phi
Flux.
Definition ScalarTransport.H:94
ScalarTransport::_runTime
autoPtr< Time > _runTime
Time.
Definition ScalarTransport.H:91
ScalarTransport::_U
autoPtr< volVectorField > _U
Velocity field.
Definition ScalarTransport.H:82
ScalarTransport::_mesh
autoPtr< fvMesh > _mesh
Mesh.
Definition ScalarTransport.H:88
ScalarTransport::Tfield
PtrList< volScalarField > Tfield
List of pointers used to form the T snapshots matrix.
Definition ScalarTransport.H:74
ScalarTransport::restart
void restart()
Function to restart the fields of the Burgers problem.
Definition ScalarTransport.C:153
ScalarTransport::ScalarTransport
ScalarTransport()
Null constructor.
Definition ScalarTransport.C:41
ScalarTransport::_fvOptions
autoPtr< fv::options > _fvOptions
fvOptions
Definition ScalarTransport.H:106
ScalarTransport::_simple
autoPtr< simpleControl > _simple
simpleControl
Definition ScalarTransport.H:109
ScalarTransport::_T
autoPtr< volScalarField > _T
Temperature field.
Definition ScalarTransport.H:80
UnsteadyProblem::checkWrite
bool checkWrite(Time &timeObject)
Function to check if the solution must be exported.
Definition UnsteadyProblem.C:37
UnsteadyProblem::startTime
scalar startTime
Start Time (initial time to start storing the snapshots).
Definition UnsteadyProblem.H:43
UnsteadyProblem::writeEvery
scalar writeEvery
Time step of the writing procedure.
Definition UnsteadyProblem.H:52
UnsteadyProblem::nextWrite
scalar nextWrite
Auxiliary variable to store the next writing instant.
Definition UnsteadyProblem.H:55
reductionProblem::counter
label counter
Counter used for the output of the full order solutions.
Definition reductionProblem.H:110
reductionProblem::offline
bool offline
Boolean variable, it is 1 if the Offline phase has already been computed, else 0.
Definition reductionProblem.H:99
reductionProblem::ITHACAdict
IOdictionary * ITHACAdict
dictionary to store input output infos
Definition reductionProblem.H:101
reductionProblem::_args
autoPtr< argList > _args
argList
Definition reductionProblem.H:103
reductionProblem::podex
bool podex
Boolean variable, it is 1 if the POD has already been computed, else 0.
Definition reductionProblem.H:97
reductionProblem::truthSolve
void truthSolve()
Perform a TruthSolve.
Definition reductionProblem.C:97
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::check_pod
bool check_pod()
Check if the POD data folder "./ITHACAoutput/POD" exists.
Definition ITHACAsystem.C:143
ITHACAutilities::check_off
bool check_off()
Check if the offline data folder "./ITHACAoutput/Offline" exists.
Definition ITHACAsystem.C:107
Generated by   doxygen 1.16.1