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
\*---------------------------------------------------------------------------*/
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.
Implementation of a parametrized full order unsteady Boussinesq problem and preparation of the the ...
{
public:
:
{}
autoPtr< volScalarField > _p_rgh
Shifted Pressure field.
UnsteadyBB()
Null constructor.
autoPtr< volVectorField > _U
Velocity field.
autoPtr< volScalarField > _p
Pressure field.
tutorial10(int argc, char *argv[])
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)
{
List<scalar> mu_now(1);
{
}
PtrList< volScalarField > Tfield
List of pointers used to form the temperature snapshots matrix.
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.
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++)
{
}
Eigen::MatrixXd mu
Row matrix of parameters.
Eigen::MatrixXi inletIndexT
void assignBC(volVectorField &s, label BC_ind, Vector< double > &value)
Assign Boundary Condition to a volVectorField.
}
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);
for (label j = 0; j < par_BC.cols(); j++)
{
}
void createSymLink(word folder)
Creates symbolic links to 0, system and constant.
}
}
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)
{
{
}
PtrList< volVectorField > Ufield_on
List of pointers used to form the temperature snapshots matrix.
PtrList< volScalarField > Tfield_on
List of pointers used to form the temperature snapshots matrix.
else
{
}
}
The liftingfunctions for temperature in this problem are determined by solving a steady state laplacian function
{
{
dimensionedScalar&
Pr =
_Pr();
volScalarField Tlift(
"Tlift" + name(k),
T);
instantList Times =
runTime.times();
Info << "Solving a lifting Problem" << endl;
scalar t1 = 1;
scalar t0 = 0;
for (label j = 0; j <
T.boundaryField().size(); j++)
{
if (j == BCind)
{
}
else if (
T.boundaryField()[BCind].type() ==
"fixedValue")
{
}
else
{
}
}
while (
simple.correctNonOrthogonal())
{
alphat.correctBoundaryConditions();
(
);
Info <<
"ExecutionTime = " <<
runTime.elapsedCpuTime() <<
" s"
<<
" ClockTime = " <<
runTime.elapsedClockTime() <<
" s"
<< nl << endl;
}
Tlift.write();
volScalarField alphaEff("alphaEff", turbulence->nu()/Pr+alphat)
autoPtr< fvMesh > _mesh
Mesh.
PtrList< volScalarField > liftfieldT
List of pointers used to form the list of lifting functions.
void assignIF(T &s, G &value)
Assign internal field condition.
autoPtr< Time > _runTime
Time.
simpleControl simple(mesh)
}
}
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:
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");
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");
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.
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;
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:
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:
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.
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;
}
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
for (
int i = 0;
i < List_of_modes.rows();
i++)
{
example.Umodes,
List_of_modes(
i, 0) + example.liftfield.size() + NmodesSUPproj);
List_of_modes(
i, 0) + example.liftfieldT.size());
example.Umodes, coeffU, List_of_modes(
i, 0));
TLmodes, coeffT, List_of_modes(
i, 0) + example.liftfieldT.size());
rec_fieldU);
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
"./ITHACAoutput/l2error");
"./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:
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
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;
HFonline2.onlineSolveFull(par_on_BC, 1,
"./ITHACAoutput/HFonline2");
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;
HFonline3.onlineSolveFull(par_on_BC, 2,
"./ITHACAoutput/HFonline3");
example.onlineSolveRead("./ITHACAoutput/Offline/");
example.onlineSolveRead("./ITHACAoutput/HFonline2/");
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
#include <chrono>
#include <math.h>
#include <iomanip>
{
public:
:
{}
{
List<scalar> mu_now(1);
{
}
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++)
{
}
}
}
}
}
void onlineSolveFull(Eigen::MatrixXd par_BC, label para_set_BC, fileName folder)
{
{
}
else
{
mkDir(folder);
for (label j = 0; j < par_BC.cols(); j++)
{
}
}
}
{
{
}
else
{
}
}
{
{
dimensionedScalar&
Pr =
_Pr();
volScalarField Tlift(
"Tlift" + name(k),
T);
instantList Times =
runTime.times();
Info << "Solving a lifting Problem" << endl;
scalar t1 = 1;
scalar t0 = 0;
for (label j = 0; j <
T.boundaryField().size(); j++)
{
if (j == BCind)
{
}
else if (
T.boundaryField()[BCind].type() ==
"fixedValue")
{
}
else
{
}
}
while (
simple.correctNonOrthogonal())
{
alphat.correctBoundaryConditions();
(
);
Info <<
"ExecutionTime = " <<
runTime.elapsedCpuTime() <<
" s"
<<
" ClockTime = " <<
runTime.elapsedClockTime() <<
" s"
<< nl << endl;
}
Tlift.write();
}
}
};
int main(
int argc,
char* argv[])
{
word par_offline_BC("./par_offline_BC");
word par_online_BC("./par_online_BC");
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;
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;
example.offlineSolve(par_off_BC);
example.liftSolveT();
example.computeLiftT(example.Tfield, example.liftfieldT, example.Tomfield);
example.podex, 0, 0,
NmodesOut);
example.podex, 0, 0,
NmodesOut);
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;
}
"./ITHACAoutput/l2error");
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());
}
Eigen::MatrixXd L2errorProjMatrixU(example.Ufield.size(), List_of_modes.rows());
Eigen::MatrixXd L2errorProjMatrixT(example.Tfield.size(), List_of_modes.rows());
for (
int i = 0;
i < List_of_modes.rows();
i++)
{
example.Umodes,
List_of_modes(
i, 0) + example.liftfield.size() + NmodesSUPproj);
List_of_modes(
i, 0) + example.liftfieldT.size());
example.Umodes, coeffU, List_of_modes(
i, 0));
TLmodes, coeffT, List_of_modes(
i, 0) + example.liftfieldT.size());
rec_fieldU);
rec_fieldT);
L2errorProjMatrixU.col(
i) = L2errorProjU;
L2errorProjMatrixT.col(
i) = L2errorProjT;
}
"./ITHACAoutput/l2error");
"./ITHACAoutput/l2error");
example.projectSUP("./Matrices", NmodesUproj, NmodesPproj, NmodesTproj,
NmodesSUPproj);
example.Tmodes.resize(NmodesTproj);
example.Umodes.resize(NmodesUproj);
reduced.nu = 0.00001;
reduced.Pr = 0.71;
reduced.tstart = 0.0;
reduced.finalTime = 10;
reduced.dt = 0.005;
Eigen::MatrixXd vel_now_BC(0, 0);
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);
}
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;
HFonline2.onlineSolveFull(par_on_BC, 1,
"./ITHACAoutput/HFonline2");
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;
HFonline3.onlineSolveFull(par_on_BC, 2,
"./ITHACAoutput/HFonline3");
example.onlineSolveRead("./ITHACAoutput/Offline/");
example.onlineSolveRead("./ITHACAoutput/HFonline2/");
example.onlineSolveRead("./ITHACAoutput/HFonline3/");
example.Ufield_on, reduced.UREC);
example.Tfield_on, reduced.TREC);
"./ITHACAoutput/l2error");
"./ITHACAoutput/l2error");
exit(0);
}
autoPtr< volScalarField > _T
Temperature field.
autoPtr< dimensionedScalar > _Prt
dimensionedScalar Prt;
autoPtr< dimensionedScalar > _Pr
dimensionedScalar Pr;
autoPtr< volScalarField > _alphat
dimensionedScalar alphat;
autoPtr< incompressible::turbulenceModel > turbulence
Turbulence model.