33#include "viscosityModel.H"
41 _args = autoPtr<argList>
43 new argList(argc, argv)
46 if (!
_args->checkRootCase())
48 Foam::FatalError.exit();
51 argList& args =
_args();
54 _simple = autoPtr<simpleControl>
65#pragma GCC diagnostic push
66#pragma GCC diagnostic ignored "-Wunused-variable"
67#include "initContinuityErrs.H"
68#pragma GCC diagnostic pop
86 "The BC method must be set to lift or penalty in ITHACAdict");
102 volScalarField&
p =
_p();
103 volVectorField&
U =
_U();
104 surfaceScalarField&
phi =
_phi();
107 IOMRFZoneList& MRF =
_MRF();
122 for (label
i = 0;
i < mu_now.size();
i++)
136 "./ITHACAoutput/Offline");
141 label NSUPmodes, label nNutModes)
147 for (label j = 0; j < cSize; j++)
153 for (label
i = 0;
i < cSize;
i++)
155 Info <<
"Filling layer number " <<
i + 1 <<
" in the matrix ct1Matrix" << endl;
159 for (label k = 0; k < cSize; k++)
169 "./ITHACAoutput/Matrices/");
171 "./ITHACAoutput/Matrices/");
173 "./ITHACAoutput/Matrices/ct1");
178 label NSUPmodes, label nNutModes)
184 for (label
i = 0;
i < cSize;
i++)
188 for (label k = 0; k < cSize; k++)
208 label NSUPmodes, label nNutModes)
214 for (label j = 0; j < cSize; j++)
220 for (label
i = 0;
i < cSize;
i++)
222 Info <<
"Filling layer number " <<
i + 1 <<
" in the matrix ct2Matrix" << endl;
226 for (label k = 0; k < cSize; k++)
236 "./ITHACAoutput/Matrices/");
238 "./ITHACAoutput/Matrices/");
240 "./ITHACAoutput/Matrices/ct2");
245 label NSUPmodes, label nNutModes)
251 for (label
i = 0;
i < cSize;
i++)
255 for (label k = 0; k < cSize; k++)
271 label NSUPmodes, label NPmodes, label nNutModes)
281 for (label k = 0; k < cSize; k++)
303 label NSUPmodes, label NPmodes, label nNutModes)
313 for (label k = 0; k < cSize; k++)
336 Eigen::MatrixXd
btMatrix(btSize, btSize);
340 for (label
i = 0;
i < btSize;
i++)
342 for (label j = 0; j < btSize; j++)
367 for (label k = 0; k <
liftfield.size(); k++)
375 for (label k = 0; k <
NUmodes; k++)
383 word bStr =
"B_" + name(
liftfield.size()) +
"_" + name(
NUmodes) +
"_" + name(
395 word btStr =
"bt_" + name(
liftfield.size()) +
"_" + name(
NUmodes) +
"_" + name(
407 word kStr =
"K_" + name(
liftfield.size()) +
"_" + name(
NUmodes) +
"_" + name(
419 word pStr =
"P_" + name(
liftfield.size()) +
"_" + name(
NUmodes) +
"_" + name(
431 word mStr =
"M_" + name(
liftfield.size()) +
"_" + name(
NUmodes) +
"_" + name(
443 word D_str =
"D_" + name(
NPmodes);
454 word bc3_str =
"BC3_" + name(
liftfield.size()) +
"_" + name(
466 word bc4_str =
"BC4_" + name(
liftfield.size()) +
"_" + name(
478 word C_str =
"C_" + name(
liftfield.size()) +
"_" + name(
NUmodes) +
"_" + name(
490 word ct1Str =
"ct1_" + name(
liftfield.size()) +
"_" + name(
503 word ct2Str =
"ct2_" + name(
liftfield.size()) +
"_" + name(
516 word G_str =
"G_" + name(
liftfield.size()) +
"_" + name(
NUmodes) +
"_" + name(
540 for (label k = 0; k <
liftfield.size(); k++)
548 for (label k = 0; k <
NUmodes; k++)
591 "./ITHACAoutput/Matrices/");
593 "./ITHACAoutput/Matrices/");
595 "./ITHACAoutput/Matrices/");
599 "./ITHACAoutput/Matrices/");
601 "./ITHACAoutput/Matrices/");
603 "./ITHACAoutput/Matrices/");
605 "./ITHACAoutput/Matrices/");
615 "./ITHACAoutput/Matrices/");
617 "./ITHACAoutput/Matrices/");
619 "./ITHACAoutput/Matrices/");
622 "./ITHACAoutput/Matrices/");
624 "./ITHACAoutput/Matrices/");
626 "./ITHACAoutput/Matrices/");
628 "./ITHACAoutput/Matrices/");
638 "./ITHACAoutput/Matrices/");
640 "./ITHACAoutput/Matrices/");
645 "./ITHACAoutput/Matrices/ct1");
647 "./ITHACAoutput/Matrices/ct2");
649 "./ITHACAoutput/Matrices/ct1PPE");
651 "./ITHACAoutput/Matrices/ct2PPE");
664 "./ITHACAoutput/Matrices/");
666 "./ITHACAoutput/Matrices/");
668 "./ITHACAoutput/Matrices/");
674 Eigen::MatrixXd weights;
678 word weightName =
"wRBF_N" + name(
i + 1) +
"_" + name(
liftfield.size()) +
"_"
683 samples[
i] =
new SPLINTER::DataTable(1, 1);
685 for (label j = 0; j <
coeffL2.cols(); j++)
693 SPLINTER::RadialBasisFunctionType::GAUSSIAN, weights);
694 std::cout <<
"Constructing RadialBasisFunction for mode " <<
i + 1 << std::endl;
698 samples[
i] =
new SPLINTER::DataTable(1, 1);
700 for (label j = 0; j <
coeffL2.cols(); j++)
706 SPLINTER::RadialBasisFunctionType::GAUSSIAN);
708 "./ITHACAoutput/weightsPPE/", weightName);
709 std::cout <<
"Constructing RadialBasisFunction for mode " <<
i + 1 << std::endl;
725 for (label k = 0; k <
liftfield.size(); k++)
733 for (label k = 0; k <
NUmodes; k++)
749 word bStr =
"B_" + name(
liftfield.size()) +
"_" + name(
NUmodes) +
"_" + name(
761 word btStr =
"bt_" + name(
liftfield.size()) +
"_" + name(
NUmodes) +
"_" + name(
773 word kStr =
"K_" + name(
liftfield.size()) +
"_" + name(
NUmodes) +
"_" + name(
785 word pStr =
"P_" + name(
liftfield.size()) +
"_" + name(
NUmodes) +
"_" + name(
797 word mStr =
"M_" + name(
liftfield.size()) +
"_" + name(
NUmodes) +
"_" + name(
809 word C_str =
"C_" + name(
liftfield.size()) +
"_" + name(
NUmodes) +
"_" + name(
821 word ct1Str =
"ct1_" + name(
liftfield.size()) +
"_" + name(
834 word ct2Str =
"ct2_" + name(
liftfield.size()) +
"_" + name(
859 for (label k = 0; k <
liftfield.size(); k++)
867 for (label k = 0; k <
NUmodes; k++)
905 "./ITHACAoutput/Matrices/");
908 "./ITHACAoutput/Matrices/");
910 "./ITHACAoutput/Matrices/");
920 "./ITHACAoutput/Matrices/");
923 "./ITHACAoutput/Matrices/");
925 "./ITHACAoutput/Matrices/");
937 "./ITHACAoutput/Matrices/ct1");
939 "./ITHACAoutput/Matrices/ct2");
950 "./ITHACAoutput/Matrices/");
952 "./ITHACAoutput/Matrices/");
954 "./ITHACAoutput/Matrices/");
960 Eigen::MatrixXd weights;
964 word weightName =
"wRBF_N" + name(
i + 1) +
"_" + name(
liftfield.size()) +
"_"
969 samples[
i] =
new SPLINTER::DataTable(1, 1);
971 for (label j = 0; j <
coeffL2.cols(); j++)
979 SPLINTER::RadialBasisFunctionType::GAUSSIAN, weights);
980 std::cout <<
"Constructing RadialBasisFunction for mode " <<
i + 1 << std::endl;
984 samples[
i] =
new SPLINTER::DataTable(1, 1);
986 for (label j = 0; j <
coeffL2.cols(); j++)
992 SPLINTER::RadialBasisFunctionType::GAUSSIAN);
994 "./ITHACAoutput/weightsSUP/", weightName);
996 "wRBF_" + name(
i),
"eigen",
997 "./ITHACAoutput/weightsSUP/"
999 std::cout <<
"Constructing RadialBasisFunction for mode " <<
i + 1 << std::endl;
#define M_Assert(Expr, Msg)
Header file of the SteadyNSTurb class.
static ITHACAparameters * getInstance()
Gets an instance of ITHACAparameters, to be used if the instance is already existing.
Eigen::MatrixXd coeffL2
The matrix of L2 projection coefficients for the eddy viscosity.
Eigen::MatrixXd btMatrix
Turbulent viscosity matrix.
Eigen::Tensor< double, 3 > cTotalPPETensor
Turbulent total viscosity tensor in the PPE equation.
label nNutModes
Number of viscoisty modes used for the projection.
Eigen::MatrixXd bTotalMatrix
Total B Matrix.
std::vector< SPLINTER::RBFSpline * > rbfSplines
Create a RBF splines for interpolation.
Eigen::MatrixXd btTurbulence(label NUmodes, label NSUPmodes)
bt added matrix for the turbulence treatement
autoPtr< volScalarField > _nut
Eddy viscosity field.
std::vector< SPLINTER::DataTable * > samples
Create a samples for interpolation.
List< Eigen::MatrixXd > turbulenceTerm1(label NUmodes, label NSUPmodes, label nNutModes)
ct1 added matrix for the turbulence treatement
PtrList< volScalarField > nutFields
List of snapshots for the solution for eddy viscosity.
List< Eigen::MatrixXd > ct2Matrix
Turbulent viscosity tensor.
Eigen::Tensor< double, 3 > ct1PPETensor
Turbulent viscosity tensor in the PPE equation.
Eigen::Tensor< double, 3 > ct1Tensor
List< Eigen::MatrixXd > ct1Matrix
Turbulent viscosity tensor.
List< Eigen::MatrixXd > turbulenceTerm2(label NUmodes, label NSUPmodes, label nNutModes)
Method to compute one of the turbulence eddy viscosity tensors.
void projectSUP(fileName folder, label NUmodes, label NPmodes, label NSUPmodes, label nNutModes)
Project using a supremizer approach.
Eigen::Tensor< double, 3 > ct2PPETensor
Turbulent viscosity tensor in the PPE equation.
Eigen::Tensor< double, 3 > turbulencePPETensor1(label NUmodes, label NSUPmodes, label NPmodes, label nNutModes)
ct1PPE added tensor for the turbulence treatement in the PPE method
Eigen::Tensor< double, 3 > cTotalTensor
Eigen::Tensor< double, 3 > turbulencePPETensor2(label NUmodes, label NSUPmodes, label NPmodes, label nNutModes)
ct2PPE added tensor for the turbulence treatement in the PPE method
Eigen::Tensor< double, 3 > turbulenceTensor1(label NUmodes, label NSUPmodes, label nNutModes)
ct1 tensor for the turbulence treatement
Eigen::Tensor< double, 3 > ct2Tensor
Eigen::Tensor< double, 3 > turbulenceTensor2(label NUmodes, label NSUPmodes, label nNutModes)
Method to compute one of the turbulence eddy viscosity tensors.
void projectPPE(fileName folder, label NUmodes, label NPmodes, label NSUPmodes, label nNutModes)
Project using a PPE approach.
word viscCoeff
The way to compute the eddy viscosity snapshots.
volScalarModes nutModes
List of POD modes for eddy viscosity.
void writeMu(List< scalar > mu_now)
Write out a list of scalar corresponding to the parameters used in the truthSolve.
Eigen::MatrixXd mu_samples
Matrix of parameters to be used for PODI, where each row corresponds to a sample point....
label counter
Counter used for the output of the full order solutions.
bool offline
Boolean variable, it is 1 if the Offline phase has already been computed, else 0.
IOdictionary * ITHACAdict
dictionary to store input output infos
Eigen::MatrixXd mu
Row matrix of parameters.
autoPtr< argList > _args
argList
bool podex
Boolean variable, it is 1 if the POD has already been computed, else 0.
void truthSolve()
Perform a TruthSolve.
List< Eigen::MatrixXd > bcVelVec
Boundary term for penalty method - vector.
scalar maxIter
Number of maximum iterations to be done for the computation of the truth solution.
label NPmodes
Number of pressure modes used for the projection.
bool supex
Boolean variable to check the existence of the supremizer modes.
Eigen::MatrixXd diffusive_term(label NUmodes, label NPmodes, label NSUPmodes)
Diffusive Term.
autoPtr< surfaceScalarField > _phi
Flux.
Eigen::MatrixXd BC4_matrix
PPE BC4.
Eigen::MatrixXd BC3_matrix
PPE BC3.
Eigen::Tensor< double, 3 > C_tensor
Diffusion term.
autoPtr< simpleControl > _simple
simpleControl
PtrList< volScalarField > Pfield
List of pointers used to form the pressure snapshots matrix.
volVectorModes supmodes
List of pointers used to form the supremizer modes.
List< Eigen::MatrixXd > bcVelocityVec(label NUmodes, label NSUPmodes)
Boundary integral modes on boundary used by the penaly method.
autoPtr< fv::options > _fvOptions
fvOptions
Eigen::MatrixXd pressure_BC4(label NPmodes, label NUmodes)
Term N° 4 given by the additional boundary condition using a PPE approach for time-dependent BCs.
autoPtr< Time > _runTime
Time.
volVectorModes Umodes
List of pointers used to form the velocity modes.
PtrList< volVectorField > Ufield
List of pointers used to form the velocity snapshots matrix.
Eigen::MatrixXd divergence_term(label NUmodes, label NPmodes, label NSUPmodes)
Divergence Term (supremizer approach)
scalar tolerance
Tolerance for the residual of the stationary problems, there is the same tolerance for velocity and p...
Eigen::MatrixXd mass_term(label NUmodes, label NPmodes, label NSUPmodes)
Mass Term.
autoPtr< fvMesh > _mesh
Mesh.
autoPtr< singlePhaseTransportModel > _laminarTransport
Laminar transport (used by turbulence model)
label NUmodes
Number of velocity modes used for the projection.
PtrList< volVectorField > liftfield
List of pointers used to form the list of lifting functions.
volVectorModes L_U_SUPmodes
List of pointers containing the total number of lift, supremizer and velocity modes.
Eigen::MatrixXd pressure_BC3(label NPmodes, label NUmodes)
Term N° 3 given by the additional boundary condition using a PPE approach.
Eigen::MatrixXd B_matrix
Diffusion term.
Eigen::MatrixXd D_matrix
Laplacian term PPE.
autoPtr< IOMRFZoneList > _MRF
MRF variable.
label NSUPmodes
Number of supremizer modes used for the projection.
Eigen::Tensor< double, 3 > gTensor
Divergence of momentum PPE.
Eigen::MatrixXd K_matrix
Gradient of pressure matrix.
Eigen::Tensor< double, 3 > convective_term_tens(label NUmodes, label NPmodes, label NSUPmodes)
Export convective term as a tensor.
autoPtr< incompressible::turbulenceModel > turbulence
Turbulence model.
Eigen::Tensor< double, 3 > divMomentum(label NUmodes, label NPmodes)
Divergence of convective term (PPE approach)
Eigen::MatrixXd pressure_gradient_term(label NUmodes, label NPmodes, label NSUPmodes)
Gradient of pressure.
List< Eigen::MatrixXd > bcVelMat
Boundary term for penalty method - matrix.
Eigen::MatrixXd P_matrix
Div of velocity.
Eigen::MatrixXd M_matrix
Mass Matrix.
Eigen::MatrixXd laplacian_pressure(label NPmodes)
Laplacian of pressure term (PPE approach)
List< Eigen::MatrixXd > bcVelocityMat(label NUmodes, label NSUPmodes)
Boundary integral modes on boundary used by the penaly method.
autoPtr< volVectorField > _U
Velocity field.
volScalarModes Pmodes
List of pointers used to form the pressure modes.
word bcMethod
Boundary Method.
autoPtr< volScalarField > _p
Pressure field.
void ReadDenseMatrix(MatrixType &Matrix, word folder, word MatrixName)
Read a dense matrix from a binary format file.
void exportSolution(GeometricField< Type, PatchField, GeoMesh > &s, fileName subfolder, fileName folder, word fieldName)
Export a field to file in a certain folder and subfolder.
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 ...
void SaveDenseMatrix(MatrixType &Matrix, word folder, word MatrixName)
Save a dense matrix to a binary format file.
void ReadDenseTensor(TensorType &Tensor, word folder, word MatrixName)
Read a dense tensor from file.
void exportTensor(Eigen::Tensor< T, 3 > tensor, word Name, word type, word folder)
Export the reduced tensor in numpy (tipo=python), matlab (tipo=matlab) and txt (tipo=eigen) format.
void SaveDenseTensor(TensorType &Tensor, word folder, word MatrixName)
Save a dense tensor to file.
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.
bool check_pod()
Check if the POD data folder "./ITHACAoutput/POD" exists.
bool check_off()
Check if the offline data folder "./ITHACAoutput/Offline" exists.
bool check_folder(word folder)
Checks if a folder exists.
bool check_file(std::string fileName)
Function that returns true if a file exists.
bool check_sup()
Check if the supremizer folder exists.
simpleControl simple(mesh)
singlePhaseTransportModel & laminarTransport
Header file of the steadyNS class.