Loading...
Searching...
No Matches
10UnsteadyBBEnclosed.C

Introduction to tutorial 10

In this tutorial an unsteady Buoyant Boussinesq (BB) 2D problem with paramerized temperature boundary conditions is implemented. The physical problem represents a differentially heated cavity. A uniform temperature is set to the left (hot) and the right (cold) sides of the cavity while the other sides are set to adiabatic. The cavity aspect ratio is 1.0. The flow is laminar and the working fluid is air with Pr = 0.7. The ambient temperature is 300 K. The hot wall, Th, has a temperature of 301.5 K, while the cold wall, Tc, is set to 298.5 K. The initial condition for the velocity is (0.0001, 0, 0) m/s.

The following image illustrates the geometry with the boundary conditions.

A detailed look into the code

In this section we explain the main steps necessary to construct the tutorial N°10

The necessary header files

First of all let's have a look at the header files that need to be included and what they are responsible for.

The header files of ITHACA-FV necessary for this tutorial are: <UnsteadyBB.H> for the full order unsteady BB problem, <ITHACAPOD.H> for the POD decomposition, <ReducedUnsteadyBB.H> for the construction of the reduced order problem, and finally <ITHACAstream.H> for some ITHACA input-output operations.

10UnsteadyBBEnclosed.C
\*---------------------------------------------------------------------------*/
#include "UnsteadyBB.H"
#include "ITHACAPOD.H"
#include "ITHACAstream.H"
Header file of the ITHACAPOD class.
Header file of the ITHACAstream class, it contains the implementation of several methods for input ou...
Header file of the ReducedUnsteadyBB class.
Header file of the UnsteadyBB class.

Definition of the tutorial10 class

We define the tutorial10 class as a child of the UnsteadyBB class. The constructor is defined with members that are the fields need to be manipulated during the resolution of the full order problem using the BuoyancyBoussinesqPimpleFoam solver. Such fields are also initialized with the same initial conditions in the solver.

class tutorial10: public UnsteadyBB
Implementation of a parametrized full order unsteady Boussinesq problem and preparation of the the ...
Definition UnsteadyBB.H:60
{
public:
explicit tutorial10(int argc, char* argv[])
:
UnsteadyBB(argc, argv),
U(_U()),
p(_p()),
T(_T())
{}
autoPtr< volScalarField > _p_rgh
Shifted Pressure field.
Definition UnsteadyBB.H:131
UnsteadyBB()
Null constructor.
Definition UnsteadyBB.C:41
autoPtr< volVectorField > _U
Velocity field.
Definition steadyNS.H:263
autoPtr< volScalarField > _p
Pressure field.
Definition steadyNS.H:260
tutorial10(int argc, char *argv[])
_T
Definition createT.H:30
volScalarField & T
Definition createT.H:46
volVectorField & U
volScalarField & p_rgh
volScalarField & p

Inside the tutorial10 class we define the offlineSolve method according to the specific parametrized problem that needs to be solved. If the offline solve has been previously performed then the method simply reads in the existing snapshots from the Offline directory. If not, the offlineSolve changes the values of the temperature boundary conditions before performing the offline solve.

void offlineSolve(Eigen::MatrixXd par_BC)
void offlineSolve(Eigen::MatrixXd par_BC)
{
List<scalar> mu_now(1);
if (offline)
{
ITHACAstream::read_fields(Ufield, U, "./ITHACAoutput/Offline/");
ITHACAstream::read_fields(Tfield, T, "./ITHACAoutput/Offline/");
}
PtrList< volScalarField > Tfield
List of pointers used to form the temperature snapshots matrix.
Definition UnsteadyBB.H:82
bool offline
Boolean variable, it is 1 if the Offline phase has already been computed, else 0.
PtrList< volVectorField > Ufield
List of pointers used to form the velocity snapshots matrix.
Definition steadyNS.H:89
void read_fields(PtrList< GeometricField< Type, PatchField, GeoMesh > > &Lfield, word Name, fileName casename, int first_snap, int n_snap)
Function to read a list of fields from the name of the field and casename.
else
{
for (label k = 0; k < par_BC.rows(); k++)
{
for (label j = 0; j < par_BC.cols(); j++)
{
for (label i = 0; i < mu.cols(); i++)
{
mu_now[0] = mu(0, i);
}
Eigen::MatrixXd mu
Row matrix of parameters.
label i
Definition pEqn.H:46
assignBC(T, inletIndexT(j, 0), par_BC(k, j));
Eigen::MatrixXi inletIndexT
void assignBC(volVectorField &s, label BC_ind, Vector< double > &value)
Assign Boundary Condition to a volVectorField.
}
truthSolve(mu_now);
void truthSolve()
Perform a TruthSolve.
}
}
}

In order to calculate the L2 error between the online solved reduced order model and the corresponding full order solution at the end of this tutorial, the full order solve can be performed with the onlineSolveFull for a parameter set for which the reduced order model has been solved online.

void onlineSolveFull(Eigen::MatrixXd par_BC, label para_set_BC, fileName folder)
void onlineSolveFull(Eigen::MatrixXd par_BC, label para_set_BC, fileName folder)
{
{
}
bool check_folder(word folder)
Checks if a folder exists.
else
{
mkDir(folder);
label i = para_set_BC;
for (label j = 0; j < par_BC.cols(); j++)
{
assignBC(T, inletIndexT(j, 0), par_BC(i, j));
}
void createSymLink(word folder)
Creates symbolic links to 0, system and constant.
truthSolve(folder);
}
}

If the full order solve has been performed previously then the existing snapshots can be read in from the specified directory.

void onlineSolveRead(fileName folder)
void onlineSolveRead(fileName folder)
{
{
}
PtrList< volVectorField > Ufield_on
List of pointers used to form the temperature snapshots matrix.
Definition UnsteadyBB.H:91
PtrList< volScalarField > Tfield_on
List of pointers used to form the temperature snapshots matrix.
Definition UnsteadyBB.H:85
else
{
}
}

The liftingfunctions for temperature in this problem are determined by solving a steady state laplacian function

void liftSolveT()
{
for (label k = 0; k < inletIndexT.rows(); k++)
{
Time& runTime = _runTime();
fvMesh& mesh = _mesh();
volScalarField T = _T();
simpleControl simple(mesh);
volScalarField& alphat = _alphat();
dimensionedScalar& Pr = _Pr();
dimensionedScalar& Prt = _Prt();
label BCind = inletIndexT(k, 0);
volScalarField Tlift("Tlift" + name(k), T);
instantList Times = runTime.times();
runTime.setTime(Times[1], 1);
Info << "Solving a lifting Problem" << endl;
scalar t1 = 1;
scalar t0 = 0;
for (label j = 0; j < T.boundaryField().size(); j++)
{
assignIF(Tlift, t0);
if (j == BCind)
{
assignBC(Tlift, j, t1);
}
else if (T.boundaryField()[BCind].type() == "fixedValue")
{
assignBC(Tlift, j, t0);
}
else
{
}
}
while (simple.correctNonOrthogonal())
{
alphat = turbulence->nut() / Prt;
alphat.correctBoundaryConditions();
volScalarField alphaEff("alphaEff", turbulence->nu() / Pr + alphat);
fvScalarMatrix TEqn
(
fvm::laplacian(alphaEff, Tlift)
);
TEqn.solve();
Info << "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
<< " ClockTime = " << runTime.elapsedClockTime() << " s"
<< nl << endl;
}
Tlift.write();
liftfieldT.append((Tlift).clone());
Foam::fvMesh & mesh
Definition createMesh.H:47
Foam::Time & runTime
Definition createTime.H:33
volScalarField alphaEff("alphaEff", turbulence->nu()/Pr+alphat)
autoPtr< fvMesh > _mesh
Mesh.
Definition UnsteadyBB.H:128
PtrList< volScalarField > liftfieldT
List of pointers used to form the list of lifting functions.
Definition UnsteadyBB.H:94
void assignIF(T &s, G &value)
Assign internal field condition.
autoPtr< Time > _runTime
Time.
Definition steadyNS.H:290
fvScalarMatrix & TEqn
Definition TEqn.H:15
volScalarField & alphat
_Pr
dimensionedScalar & Prt
dimensionedScalar & Pr
_Prt
simpleControl simple(mesh)
turbulence
}
}

Definition of the main function

In this section we show the definition of the main function. First we construct the object "example" of type tutorial10:

tutorial10 example(argc, argv);

The parameter sets for the temperature BC are read in from file, both for the offline and online phase. Each row corresponds to a parameter set and each column to a specific boundary.

word par_offline_BC("./par_offline_BC");
Eigen::MatrixXd par_off_BC = ITHACAstream::readMatrix(par_offline_BC);
List< Eigen::MatrixXd > readMatrix(word folder, word mat_name)
Read a three dimensional matrix from a txt file in Eigen format.
word par_online_BC("./par_online_BC");
Eigen::MatrixXd par_on_BC = ITHACAstream::readMatrix(par_online_BC);

Then we parse the ITHACAdict file to determine the number of modes to be written out and also the ones to be used for projection of the velocity, pressure, temperature and the supremizer:

Class for the definition of some general parameters, the parameters must be defined from the file ITH...
static ITHACAparameters * getInstance()
Gets an instance of ITHACAparameters, to be used if the instance is already existing.
ITHACAparameters * para
Definition steadyNS.H:82
example._runTime());
int NmodesUproj = para->ITHACAdict->lookupOrDefault<int>("NmodesUproj", 5);
int NmodesPproj = para->ITHACAdict->lookupOrDefault<int>("NmodesPproj", 5);
int NmodesTproj = para->ITHACAdict->lookupOrDefault<int>("NmodesTproj", 5);
int NmodesSUPproj = para->ITHACAdict->lookupOrDefault<int>("NmodesSUPproj", 5);
int NmodesOut = para->ITHACAdict->lookupOrDefault<int>("NmodesOut", 15);
IOdictionary * ITHACAdict
Dictionary for input objects from file.

We note that a default value can be assigned in case the parser did not find the corresponding string in the ITHACAdict file.

Now the kinematic viscosity is set to 0.00001 m^2/s. It is possible to parametrize the viscosity. For more info take a look at tutorial04.

example.Pnumber = 1;
example.Tnumber = 1;
example.setParameters();
example.mu_range(0, 0) = 0.00001;
example.mu_range(0, 1) = 0.00001;
// Generate equispaced samples inside the parameter range
example.genEquiPar();

After that we set the boundaries where we have a non homogeneous BC for temperature. Patch 1 corresponds to the hot boundary and Patch 2 to the cold boundary.

example.inletIndexT.resize(2, 1);
example.inletIndexT << 1, 2;

Furthermore, we set the parameters for the time integration. In this case the simulation time is 10 seconds, with a step size = 0.005 seconds, and the data is written every 0.01 seconds, i.e.

example.startTime = 0.0;
example.finalTime = 10.0;
example.timeStep = 0.005;
example.writeEvery = 0.01;

Now we are ready to perform the offline stage:

example.offlineSolve(par_off_BC);

No lifting function has to be computed for velocity as the velocity is zero everywhere on the boundary. For the temperature, the lifting functions are determined by:

example.liftSolveT();

Then we create homogenuous basis functions for the temperature:

example.computeLiftT(example.Tfield, example.liftfieldT, example.Tomfield);

And after that, we obtain the modes for velocity and temperature:

ITHACAPOD::getModes(example.Ufield, example.Umodes, example._U().name(),
void getModes(PtrList< GeometricField< Type, PatchField, GeoMesh > > &snapshots, PtrList< GeometricField< Type, PatchField, GeoMesh > > &modes, word fieldName, bool podex, bool supex, bool sup, label nmodes, bool correctBC)
Computes the bases or reads them for a field.
Definition ITHACAPOD.C:93
ITHACAPOD::getModes(example.Tomfield, example.Tmodes, example._T().name(),

Before continuiting, the projection error is calculated for certain number of modes for both temperature and velocity.

Eigen::MatrixXd List_of_modes(NmodesOut - 5, 1);
for (int i = 0; i < List_of_modes.rows(); i++)
{
List_of_modes(i, 0) = i + 1;
}
ITHACAstream::exportMatrix(List_of_modes, "List_of_modes", "eigen",
void exportMatrix(Eigen::Matrix< T, -1, dim > &matrix, word Name, word type, word folder)
Export the reduced matrices in numpy (type=python), matlab (type=matlab) and txt (type=eigen) format ...

The temperature modes are created locally

PtrList<volScalarField> TLmodes;
for (label k = 0; k < example.liftfieldT.size(); k++)
{
TLmodes.append((example.liftfieldT[k]).clone());
}
for (label k = 0; k < List_of_modes.size(); k++)
{
TLmodes.append((example.Tmodes[k]).clone());
}

and the projection onto the POD modes is performed with:

Eigen::MatrixXd L2errorProjMatrixU(example.Ufield.size(), List_of_modes.rows());
Eigen::MatrixXd L2errorProjMatrixT(example.Tfield.size(), List_of_modes.rows());

Finally the L2 error between full order solution and projection of the basis are calculated

// Calculate the coefficients and L2 error and store the error in a matrix for each number of modes
for (int i = 0; i < List_of_modes.rows(); i++)
{
Eigen::MatrixXd coeffU = ITHACAutilities::getCoeffs(example.Ufield,
example.Umodes,
List_of_modes(i, 0) + example.liftfield.size() + NmodesSUPproj);
Eigen::MatrixXd coeffT = ITHACAutilities::getCoeffs(example.Tfield, TLmodes,
List_of_modes(i, 0) + example.liftfieldT.size());
PtrList<volVectorField> rec_fieldU = ITHACAutilities::reconstructFromCoeff(
example.Umodes, coeffU, List_of_modes(i, 0));
PtrList<volScalarField> rec_fieldT = ITHACAutilities::reconstructFromCoeff(
TLmodes, coeffT, List_of_modes(i, 0) + example.liftfieldT.size());
Eigen::MatrixXd L2errorProjU = ITHACAutilities::errorL2Rel(example.Ufield,
rec_fieldU);
Eigen::MatrixXd L2errorProjT = ITHACAutilities::errorL2Rel(example.Tfield,
rec_fieldT);
L2errorProjMatrixU.col(i) = L2errorProjU;
L2errorProjMatrixT.col(i) = L2errorProjT;
}
Eigen::VectorXd getCoeffs(GeometricField< Type, PatchField, GeoMesh > &snapshot, PtrList< GeometricField< Type, PatchField, GeoMesh > > &modes, label Nmodes, bool consider_volumes)
Projects a snapshot on a basis function and gets the coefficients of the projection.
PtrList< GeometricField< Type, PatchField, GeoMesh > > reconstructFromCoeff(PtrList< GeometricField< Type, PatchField, GeoMesh > > &modes, Eigen::MatrixXd &coeff_matrix, label Nmodes)
Exact reconstruction using a certain number of modes for a list of fields and the projection coeffici...
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.

and exported

ITHACAstream::exportMatrix(L2errorProjMatrixU, "L2errorProjMatrixU", "eigen",
"./ITHACAoutput/l2error");
ITHACAstream::exportMatrix(L2errorProjMatrixT, "L2errorProjMatrixT", "eigen",
"./ITHACAoutput/l2error");

Then the projection onto the POD modes is performed to get the reduced matrices

example.projectSUP("./Matrices", NmodesUproj, NmodesPproj, NmodesTproj,
NmodesSUPproj);

and the modes are resized to the number for which the projection has been performed

example.Tmodes.resize(NmodesTproj);
example.Umodes.resize(NmodesUproj);

Now that we obtained all the necessary information from the POD decomposition and the reduced matrices, we are ready to construct the dynamical system for the reduced order model (ROM). We proceed by constructing the object "reduced" of type ReducedUnsteadyBB:

ReducedUnsteadyBB reduced(example);
Class where it is implemented a reduced problem for the unsteady Navier-stokes problem.

We can use the new constructed ROM to perform the online procedure, from which we can simulate the problem at new set of parameters. For instance, we solve the problem for 10 seconds of physical time:

reduced.nu = 0.00001;
reduced.Pr = 0.71;
reduced.tstart = 0.0;
reduced.finalTime = 10;
reduced.dt = 0.005;

And then we can use the new constructed ROM to perform the online procedure, from which we can simulate the problem at new set of parameters. As the velocity is homogenous on the boundary, it is not parametrized and therefore an empty matrix is created:

Eigen::MatrixXd vel_now_BC(0, 0);

Lastly the online procedure is performed for all temperature boundary sets as defined in the par_online_BC file. And he ROM solution is reconstructed and exported:

for (label k = 0; k < (par_on_BC.rows()); k++)
{
Eigen::MatrixXd temp_now_BC(2, 1);
temp_now_BC(0, 0) = par_on_BC(k, 0);
temp_now_BC(1, 0) = par_on_BC(k, 1);
reduced.solveOnline_sup(temp_now_BC, vel_now_BC, k, par_on_BC.rows());
reduced.reconstruct_sup("./ITHACAoutput/ReconstructionSUP", 2);
}

For the second and third parameter set the full order solution is calculated

// Performing full order simulation for second parameter set - temp_BC
tutorial10 HFonline2(argc, argv);
HFonline2.Pnumber = 1;
HFonline2.Tnumber = 1;
HFonline2.setParameters();
HFonline2.mu_range(0, 0) = 0.00001;
HFonline2.mu_range(0, 1) = 0.00001;
HFonline2.genEquiPar();
HFonline2.inletIndexT.resize(2, 1);
HFonline2.inletIndexT << 1, 2;
HFonline2.startTime = 0.0;
HFonline2.finalTime = 10;
HFonline2.timeStep = 0.005;
HFonline2.writeEvery = 0.01;
// Reconstruct the online solution
HFonline2.onlineSolveFull(par_on_BC, 1,
"./ITHACAoutput/HFonline2");
// Performing full order simulation for third parameter set - temp_BC
tutorial10 HFonline3(argc, argv);
HFonline3.Pnumber = 1;
HFonline3.Tnumber = 1;
HFonline3.setParameters();
HFonline3.mu_range(0, 0) = 0.00001;
HFonline3.mu_range(0, 1) = 0.00001;
HFonline3.genEquiPar();
HFonline3.inletIndexT.resize(2, 1);
HFonline3.inletIndexT << 1, 2;
HFonline3.startTime = 0.0;
HFonline3.finalTime = 10.0;
HFonline3.timeStep = 0.005;
HFonline3.writeEvery = 0.01;
// Reconstruct the online solution
HFonline3.onlineSolveFull(par_on_BC, 2,
"./ITHACAoutput/HFonline3");
// Reading in the high-fidelity solutions for the parameter set
// for which the offline solve has been performed
example.onlineSolveRead("./ITHACAoutput/Offline/");
// Reading in the high-fidelity solutions for the second parameter set
example.onlineSolveRead("./ITHACAoutput/HFonline2/");
// Reading in the high-fidelity solutions for the second parameter set
example.onlineSolveRead("./ITHACAoutput/HFonline3/");

and the solutions are read in:

Finally the L2 error between full and reduced order solutions is calculated

The plain program

Here's the plain code

/*---------------------------------------------------------------------------*\
Copyright (C) 2017 by the ITHACA-FV authors
License
This file is part of ITHACA-FV
ITHACA-FV is free software: you can redistribute it and/or modify
it under the terms of the GNU Lesser General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
ITHACA-FV is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU Lesser General Public License for more details.
You should have received a copy of the GNU Lesser General Public License
along with ITHACA-FV. If not, see <http://www.gnu.org/licenses/>.
Description
Example of Boussinesq approximation for two way coupling NS-momentum equation
and heat transport equation for enclosed flows.
SourceFiles
10UnsteadyBBEnclosed.C
\*---------------------------------------------------------------------------*/
#include "UnsteadyBB.H"
#include "ITHACAPOD.H"
#include "ITHACAstream.H"
#include <chrono>
#include <math.h>
#include <iomanip>
class tutorial10: public UnsteadyBB
{
public:
explicit tutorial10(int argc, char* argv[])
:
UnsteadyBB(argc, argv),
U(_U()),
p(_p()),
T(_T())
{}
// Fields To Perform
volVectorField& U;
volScalarField& p;
volScalarField& p_rgh;
volScalarField& T;
void offlineSolve(Eigen::MatrixXd par_BC)
{
List<scalar> mu_now(1);
if (offline)
{
ITHACAstream::read_fields(Ufield, U, "./ITHACAoutput/Offline/");
ITHACAstream::read_fields(Tfield, T, "./ITHACAoutput/Offline/");
}
else
{
for (label k = 0; k < par_BC.rows(); k++)
{
for (label j = 0; j < par_BC.cols(); j++)
{
for (label i = 0; i < mu.cols(); i++)
{
mu_now[0] = mu(0, i);
}
assignBC(T, inletIndexT(j, 0), par_BC(k, j));
}
truthSolve(mu_now);
}
}
}
void onlineSolveFull(Eigen::MatrixXd par_BC, label para_set_BC, fileName folder)
{
{
}
else
{
mkDir(folder);
label i = para_set_BC;
for (label j = 0; j < par_BC.cols(); j++)
{
assignBC(T, inletIndexT(j, 0), par_BC(i, j));
}
truthSolve(folder);
}
}
void onlineSolveRead(fileName folder)
{
{
}
else
{
}
}
// Method to compute the lifting function for temperature
void liftSolveT()
{
for (label k = 0; k < inletIndexT.rows(); k++)
{
Time& runTime = _runTime();
fvMesh& mesh = _mesh();
volScalarField T = _T();
simpleControl simple(mesh);
volScalarField& alphat = _alphat();
dimensionedScalar& Pr = _Pr();
dimensionedScalar& Prt = _Prt();
label BCind = inletIndexT(k, 0);
volScalarField Tlift("Tlift" + name(k), T);
instantList Times = runTime.times();
runTime.setTime(Times[1], 1);
Info << "Solving a lifting Problem" << endl;
scalar t1 = 1;
scalar t0 = 0;
for (label j = 0; j < T.boundaryField().size(); j++)
{
assignIF(Tlift, t0);
if (j == BCind)
{
assignBC(Tlift, j, t1);
}
else if (T.boundaryField()[BCind].type() == "fixedValue")
{
assignBC(Tlift, j, t0);
}
else
{
}
}
while (simple.correctNonOrthogonal())
{
alphat = turbulence->nut() / Prt;
alphat.correctBoundaryConditions();
volScalarField alphaEff("alphaEff", turbulence->nu() / Pr + alphat);
fvScalarMatrix TEqn
(
fvm::laplacian(alphaEff, Tlift)
);
TEqn.solve();
Info << "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
<< " ClockTime = " << runTime.elapsedClockTime() << " s"
<< nl << endl;
}
Tlift.write();
liftfieldT.append((Tlift).clone());
}
}
};
/*---------------------------------------------------------------------------*\
Starting the MAIN
\*---------------------------------------------------------------------------*/
int main(int argc, char* argv[])
{
// Construct the tutorial object
tutorial10 example(argc, argv);
// the offline samples for the boundary conditions
word par_offline_BC("./par_offline_BC");
Eigen::MatrixXd par_off_BC = ITHACAstream::readMatrix(par_offline_BC);
// the samples which will be used for setting the boundary condition in the online stage
word par_online_BC("./par_online_BC");
Eigen::MatrixXd par_on_BC = ITHACAstream::readMatrix(par_online_BC);
// Read some parameters from file
example._runTime());
int NmodesUproj = para->ITHACAdict->lookupOrDefault<int>("NmodesUproj", 5);
int NmodesPproj = para->ITHACAdict->lookupOrDefault<int>("NmodesPproj", 5);
int NmodesTproj = para->ITHACAdict->lookupOrDefault<int>("NmodesTproj", 5);
int NmodesSUPproj = para->ITHACAdict->lookupOrDefault<int>("NmodesSUPproj", 5);
int NmodesOut = para->ITHACAdict->lookupOrDefault<int>("NmodesOut", 15);
example.Pnumber = 1;
example.Tnumber = 1;
example.setParameters();
example.mu_range(0, 0) = 0.00001;
example.mu_range(0, 1) = 0.00001;
// Generate equispaced samples inside the parameter range
example.genEquiPar();
example.inletIndexT.resize(2, 1);
example.inletIndexT << 1, 2;
example.startTime = 0.0;
example.finalTime = 10.0;
example.timeStep = 0.005;
example.writeEvery = 0.01;
// Perform the Offline Solve;
example.offlineSolve(par_off_BC);
// Search the lift function for the temperature
example.liftSolveT();
// Create homogeneous basis functions for temperature
example.computeLiftT(example.Tfield, example.liftfieldT, example.Tomfield);
// Perform a POD decomposition for velocity temperature and pressure fields
ITHACAPOD::getModes(example.Ufield, example.Umodes, example._U().name(),
example.podex, 0, 0,
NmodesOut);
ITHACAPOD::getModes(example.Tomfield, example.Tmodes, example._T().name(),
example.podex, 0, 0,
NmodesOut);
// Create a list with number of modes for which the projection needs to be performed
Eigen::MatrixXd List_of_modes(NmodesOut - 5, 1);
for (int i = 0; i < List_of_modes.rows(); i++)
{
List_of_modes(i, 0) = i + 1;
}
// Export with number of modes for which the projection needs to be performed
ITHACAstream::exportMatrix(List_of_modes, "List_of_modes", "eigen",
"./ITHACAoutput/l2error");
// Create locally the temperature modes
PtrList<volScalarField> TLmodes;
for (label k = 0; k < example.liftfieldT.size(); k++)
{
TLmodes.append((example.liftfieldT[k]).clone());
}
for (label k = 0; k < List_of_modes.size(); k++)
{
TLmodes.append((example.Tmodes[k]).clone());
}
// Perform the projection for all number of modes in list List_of_modes
Eigen::MatrixXd L2errorProjMatrixU(example.Ufield.size(), List_of_modes.rows());
Eigen::MatrixXd L2errorProjMatrixT(example.Tfield.size(), List_of_modes.rows());
// Calculate the coefficients and L2 error and store the error in a matrix for each number of modes
for (int i = 0; i < List_of_modes.rows(); i++)
{
Eigen::MatrixXd coeffU = ITHACAutilities::getCoeffs(example.Ufield,
example.Umodes,
List_of_modes(i, 0) + example.liftfield.size() + NmodesSUPproj);
Eigen::MatrixXd coeffT = ITHACAutilities::getCoeffs(example.Tfield, TLmodes,
List_of_modes(i, 0) + example.liftfieldT.size());
PtrList<volVectorField> rec_fieldU = ITHACAutilities::reconstructFromCoeff(
example.Umodes, coeffU, List_of_modes(i, 0));
PtrList<volScalarField> rec_fieldT = ITHACAutilities::reconstructFromCoeff(
TLmodes, coeffT, List_of_modes(i, 0) + example.liftfieldT.size());
Eigen::MatrixXd L2errorProjU = ITHACAutilities::errorL2Rel(example.Ufield,
rec_fieldU);
Eigen::MatrixXd L2errorProjT = ITHACAutilities::errorL2Rel(example.Tfield,
rec_fieldT);
L2errorProjMatrixU.col(i) = L2errorProjU;
L2errorProjMatrixT.col(i) = L2errorProjT;
}
// Export the matrix containing the error
ITHACAstream::exportMatrix(L2errorProjMatrixU, "L2errorProjMatrixU", "eigen",
"./ITHACAoutput/l2error");
ITHACAstream::exportMatrix(L2errorProjMatrixT, "L2errorProjMatrixT", "eigen",
"./ITHACAoutput/l2error");
// Get reduced matrices
example.projectSUP("./Matrices", NmodesUproj, NmodesPproj, NmodesTproj,
NmodesSUPproj);
// Resize the modes for projection
example.Tmodes.resize(NmodesTproj);
example.Umodes.resize(NmodesUproj);
// Online part
ReducedUnsteadyBB reduced(example);
// Set values of the online solve
reduced.nu = 0.00001;
reduced.Pr = 0.71;
reduced.tstart = 0.0;
reduced.finalTime = 10;
reduced.dt = 0.005;
// No parametrization of velocity on boundary
Eigen::MatrixXd vel_now_BC(0, 0);
// Set the online temperature BC and solve reduced model
for (label k = 0; k < (par_on_BC.rows()); k++)
{
Eigen::MatrixXd temp_now_BC(2, 1);
temp_now_BC(0, 0) = par_on_BC(k, 0);
temp_now_BC(1, 0) = par_on_BC(k, 1);
reduced.solveOnline_sup(temp_now_BC, vel_now_BC, k, par_on_BC.rows());
reduced.reconstruct_sup("./ITHACAoutput/ReconstructionSUP", 2);
}
// Performing full order simulation for second parameter set - temp_BC
tutorial10 HFonline2(argc, argv);
HFonline2.Pnumber = 1;
HFonline2.Tnumber = 1;
HFonline2.setParameters();
HFonline2.mu_range(0, 0) = 0.00001;
HFonline2.mu_range(0, 1) = 0.00001;
HFonline2.genEquiPar();
HFonline2.inletIndexT.resize(2, 1);
HFonline2.inletIndexT << 1, 2;
HFonline2.startTime = 0.0;
HFonline2.finalTime = 10;
HFonline2.timeStep = 0.005;
HFonline2.writeEvery = 0.01;
// Reconstruct the online solution
HFonline2.onlineSolveFull(par_on_BC, 1,
"./ITHACAoutput/HFonline2");
// Performing full order simulation for third parameter set - temp_BC
tutorial10 HFonline3(argc, argv);
HFonline3.Pnumber = 1;
HFonline3.Tnumber = 1;
HFonline3.setParameters();
HFonline3.mu_range(0, 0) = 0.00001;
HFonline3.mu_range(0, 1) = 0.00001;
HFonline3.genEquiPar();
HFonline3.inletIndexT.resize(2, 1);
HFonline3.inletIndexT << 1, 2;
HFonline3.startTime = 0.0;
HFonline3.finalTime = 10.0;
HFonline3.timeStep = 0.005;
HFonline3.writeEvery = 0.01;
// Reconstruct the online solution
HFonline3.onlineSolveFull(par_on_BC, 2,
"./ITHACAoutput/HFonline3");
// Reading in the high-fidelity solutions for the parameter set
// for which the offline solve has been performed
example.onlineSolveRead("./ITHACAoutput/Offline/");
// Reading in the high-fidelity solutions for the second parameter set
example.onlineSolveRead("./ITHACAoutput/HFonline2/");
// Reading in the high-fidelity solutions for the second parameter set
example.onlineSolveRead("./ITHACAoutput/HFonline3/");
// Calculate error between online- and corresponding full order solution
Eigen::MatrixXd L2errorMatrixU = ITHACAutilities::errorL2Rel(
example.Ufield_on, reduced.UREC);
Eigen::MatrixXd L2errorMatrixT = ITHACAutilities::errorL2Rel(
example.Tfield_on, reduced.TREC);
//Export the matrix containing the error
ITHACAstream::exportMatrix(L2errorMatrixU, "L2errorMatrixU", "eigen",
"./ITHACAoutput/l2error");
ITHACAstream::exportMatrix(L2errorMatrixT, "L2errorMatrixT", "eigen",
"./ITHACAoutput/l2error");
exit(0);
}
autoPtr< volScalarField > _T
Temperature field.
Definition UnsteadyBB.H:134
autoPtr< dimensionedScalar > _Prt
dimensionedScalar Prt;
Definition UnsteadyBB.H:146
autoPtr< dimensionedScalar > _Pr
dimensionedScalar Pr;
Definition UnsteadyBB.H:143
autoPtr< volScalarField > _alphat
dimensionedScalar alphat;
Definition UnsteadyBB.H:149
volScalarField & p
volScalarField & T
volScalarField & p_rgh
volVectorField & U
autoPtr< incompressible::turbulenceModel > turbulence
Turbulence model.
Definition unsteadyNS.H:76
int main(int argc, char *argv[])