35 #include <boost/core/ignore_unused.hpp>
55 #include <boost/format.hpp>
67 std::string EquationSystem::equationSystemTypeLookupIds[2] = {
68 LibUtilities::SessionReader::RegisterEnumValue(
"DEALIASING",
"True", 0),
69 LibUtilities::SessionReader::RegisterEnumValue(
"DEALIASING",
"False", 1)};
96 EquationSystem::EquationSystem(
99 : m_comm(pSession->GetComm()), m_session(pSession), m_graph(pGraph),
103 const vector<std::string> filenames =
m_session->GetFilenames();
105 for (
int i = 0; i < filenames.size(); ++i)
107 string sessionname =
"SessionName";
108 sessionname += boost::lexical_cast<std::string>(i);
144 if (0 ==
m_comm->GetRank())
149 if (
m_session->DefinesSolverInfo(
"HOMOGENEOUS"))
151 std::string HomoStr =
m_session->GetSolverInfo(
"HOMOGENEOUS");
154 if ((HomoStr ==
"HOMOGENEOUS1D") || (HomoStr ==
"Homogeneous1D") ||
155 (HomoStr ==
"1D") || (HomoStr ==
"Homo1D"))
161 if (
m_session->DefinesSolverInfo(
"ModeType"))
163 m_session->MatchSolverInfo(
"ModeType",
"SingleMode",
167 m_session->MatchSolverInfo(
"ModeType",
"MultipleModes",
172 if (
m_session->DefinesSolverInfo(
"ModeType"))
188 ASSERTL0(
false,
"SolverInfo ModeType not valid");
197 if ((HomoStr ==
"HOMOGENEOUS2D") || (HomoStr ==
"Homogeneous2D") ||
198 (HomoStr ==
"2D") || (HomoStr ==
"Homo2D"))
208 if ((HomoStr ==
"HOMOGENEOUS3D") || (HomoStr ==
"Homogeneous3D") ||
209 (HomoStr ==
"3D") || (HomoStr ==
"Homo3D"))
231 m_session->MatchSolverInfo(
"SPECTRALHPDEALIASING",
"True",
235 m_session->MatchSolverInfo(
"SPECTRALHPDEALIASING",
"On",
241 if (
m_session->DefinesSolverInfo(
"PROJECTION"))
243 std::string ProjectStr =
m_session->GetSolverInfo(
"PROJECTION");
245 if ((ProjectStr ==
"Continuous") || (ProjectStr ==
"Galerkin") ||
246 (ProjectStr ==
"CONTINUOUS") || (ProjectStr ==
"GALERKIN"))
250 else if ((ProjectStr ==
"MixedCGDG") ||
251 (ProjectStr ==
"Mixed_CG_Discontinuous"))
255 else if (ProjectStr ==
"DisContinuous")
261 ASSERTL0(
false,
"PROJECTION value not recognised");
266 cerr <<
"Projection type not specified in SOLVERINFO,"
267 "defaulting to continuous Galerkin"
276 int nvariables =
m_session->GetVariables().size();
277 bool DeclareCoeffPhysArrays =
true;
305 for (i = 0; i <
m_fields.size(); i++)
308 MultiRegions ::ContField3DHomogeneous2D>::
309 AllocateSharedPtr(
m_session, BkeyY, BkeyZ,
317 for (i = 0; i <
m_fields.size(); i++)
343 for (i = 0; i <
m_fields.size(); i++)
346 MultiRegions ::ContField3DHomogeneous1D>::
369 for (i = 0; i <
m_fields.size(); i++)
371 if (
m_session->GetVariable(i).compare(
"w") == 0)
375 ContField3DHomogeneous1D>::
386 ContField3DHomogeneous1D>::
403 for (i = 0; i <
m_fields.size(); i++)
406 MultiRegions ::ContField3DHomogeneous1D>::
422 DeclareCoeffPhysArrays,
425 for (i = 1; i <
m_fields.size(); i++)
427 if (
m_graph->SameExpansionInfo(
436 DeclareCoeffPhysArrays,
446 DeclareCoeffPhysArrays,
477 DeclareCoeffPhysArrays,
481 for (i = 1; i <
m_fields.size(); i++)
483 if (
m_graph->SameExpansionInfo(
492 DeclareCoeffPhysArrays,
502 DeclareCoeffPhysArrays,
521 for (i = 1; i <
m_fields.size(); ++i)
529 ASSERTL0(
false,
"Expansion dimension not recognised");
552 for (i = 0; i <
m_fields.size(); i++)
555 MultiRegions ::DisContField3DHomogeneous2D>::
556 AllocateSharedPtr(
m_session, BkeyY, BkeyZ,
564 for (i = 0; i <
m_fields.size(); i++)
585 for (i = 0; i <
m_fields.size(); i++)
588 MultiRegions ::DisContField3DHomogeneous1D>::
597 for (i = 0; i <
m_fields.size(); i++)
615 "3D fully periodic problems not implemented yet");
619 for (i = 0; i <
m_fields.size(); i++)
631 ASSERTL0(
false,
"Expansion dimension not recognised");
663 "Only one of IO_CheckTime and IO_CheckSteps "
719 for (
int i = 0; i < nvariables; ++i)
722 m_fields[i]->EvaluateBoundaryConditions(time, varName);
741 if (
m_fields[field]->GetPhysState() ==
false)
747 if (exactsoln.size())
752 else if (
m_session->DefinesFunction(
"ExactSolution"))
767 if (Normalised ==
true)
772 L2error =
sqrt(L2error * L2error / Vol);
797 if (
m_fields[field]->GetPhysState() ==
false)
803 if (exactsoln.size())
808 else if (
m_session->DefinesFunction(
"ExactSolution"))
827 Linferror = L2INF[1];
847 LibUtilities::eGaussRadauMAlpha1Beta0);
870 Tkeys.push_back(BkeyT1);
871 Tkeys.push_back(BkeyT2);
874 Qkeys.push_back(BkeyQ1);
875 Qkeys.push_back(BkeyQ2);
883 int ErrorCoordim = ErrorExp->GetCoordim(0);
884 int ErrorNq = ErrorExp->GetTotPoints();
890 switch (ErrorCoordim)
893 ErrorExp->GetCoords(ErrorXc0);
896 ErrorExp->GetCoords(ErrorXc0, ErrorXc1);
899 ErrorExp->GetCoords(ErrorXc0, ErrorXc1, ErrorXc2);
903 m_session->GetFunction(
"ExactSolution", field);
908 exSol->Evaluate(ErrorXc0, ErrorXc1, ErrorXc2,
m_time, ErrorSol);
912 ErrorExp->BwdTrans(
m_fields[field]->GetCoeffs(), ErrorExp->UpdatePhys());
914 L2INF[0] = ErrorExp->L2(ErrorExp->GetPhys(), ErrorSol);
915 L2INF[1] = ErrorExp->Linf(ErrorExp->GetPhys(), ErrorSol);
927 bool dumpInitialConditions,
930 boost::ignore_unused(initialtime);
932 if (
m_session->GetComm()->GetRank() == 0)
934 cout <<
"Initial Conditions:" << endl;
937 if (
m_session->DefinesFunction(
"InitialConditions"))
945 for (
int i = 0; i <
m_fields.size(); ++i)
954 if (
m_session->GetComm()->GetRank() == 0)
957 for (
int i = 0; i <
m_fields.size(); ++i)
959 std::string varName =
m_session->GetVariable(i);
960 cout <<
" - Field " << varName <<
": "
962 ->Describe(varName, domain)
970 for (
int i = 0; i <
m_fields.size(); i++)
976 if (
m_session->GetComm()->GetRank() == 0)
978 cout <<
" - Field " <<
m_session->GetVariable(i)
979 <<
": 0 (default)" << endl;
996 "ExactSolution array size mismatch.");
998 if (
m_session->DefinesFunction(
"ExactSolution"))
1001 ->Evaluate(
m_session->GetVariable(field), outfield, time);
1063 for (
int i = 0; i <
m_fields.size(); i++)
1074 for (
int i = 0; i <
m_fields.size(); i++)
1088 std::string outname =
1100 std::vector<std::string> &variables)
1102 std::string outname =
1104 WriteFld(outname, field, fieldcoeffs, variables);
1113 std::string outname =
1114 m_sessionName +
"_BaseFlow_" + boost::lexical_cast<std::string>(n);
1125 std::vector<Array<OneD, NekDouble>> fieldcoeffs(
m_fields.size());
1126 std::vector<std::string> variables(
m_fields.size());
1128 for (
int i = 0; i <
m_fields.size(); ++i)
1132 fieldcoeffs[i] =
m_fields[i]->UpdateCoeffs();
1137 m_fields[0]->ExtractCoeffsToCoeffs(
1158 std::vector<std::string> &variables)
1160 std::vector<LibUtilities::FieldDefinitionsSharedPtr> FieldDef =
1161 field->GetFieldDefinitions();
1162 std::vector<std::vector<NekDouble>> FieldData(FieldDef.size());
1165 for (
int j = 0; j < fieldcoeffs.size(); ++j)
1167 for (
int i = 0; i < FieldDef.size(); ++i)
1170 FieldDef[i]->m_fields.push_back(variables[j]);
1171 field->AppendFieldData(FieldDef[i], FieldData[i], fieldcoeffs[j]);
1185 boost::lexical_cast<std::string>(
m_nchk);
1195 mapping->Output(fieldMetaDataMap, outname);
1201 std::vector<std::string> vSuffix = {
"_x",
"_y",
"_z"};
1202 for (
int i = 0; i < 3; ++i)
1204 std::string sTheta =
"Theta" + vSuffix[i];
1210 #ifdef NEKTAR_DISABLE_BACKUPS
1211 bool backup =
false;
1216 m_fld->Write(outname, FieldDef, FieldData, fieldMetaDataMap, backup);
1226 const std::string &infile,
1229 std::vector<LibUtilities::FieldDefinitionsSharedPtr> FieldDef;
1230 std::vector<std::vector<NekDouble>> FieldData;
1233 field_fld->Import(infile, FieldDef, FieldData);
1236 for (
int j = 0; j < pFields.size(); ++j)
1240 for (
int i = 0; i < FieldDef.size(); ++i)
1243 std::string(
"Order of ") + infile +
1244 std::string(
" data and that defined in "
1245 "m_boundaryconditions differs"));
1247 pFields[j]->ExtractDataToCoeffs(FieldDef[i], FieldData[i],
1249 pFields[j]->UpdateCoeffs());
1251 pFields[j]->BwdTrans(pFields[j]->GetCoeffs(), pFields[j]->UpdatePhys());
1264 const std::string &infile,
1267 std::vector<LibUtilities::FieldDefinitionsSharedPtr> FieldDef;
1268 std::vector<std::vector<NekDouble>> FieldData;
1275 ndomains * nvariables == pFields.size(),
1276 "Number of fields does not match the number of variables and domains");
1279 for (
int j = 0; j < ndomains; ++j)
1281 for (
int i = 0; i < nvariables; ++i)
1284 pFields[j * nvariables + i]->UpdateCoeffs(), 1);
1286 for (
int n = 0; n < FieldDef.size(); ++n)
1289 std::string(
"Order of ") + infile +
1290 std::string(
" data and that defined in "
1291 "m_boundaryconditions differs"));
1293 pFields[j * nvariables + i]->ExtractDataToCoeffs(
1294 FieldDef[n], FieldData[n], FieldDef[n]->
m_fields[i],
1295 pFields[j * nvariables + i]->UpdateCoeffs());
1297 pFields[j * nvariables + i]->BwdTrans(
1298 pFields[j * nvariables + i]->GetCoeffs(),
1299 pFields[j * nvariables + i]->UpdatePhys());
1311 std::string &pFieldName)
1313 std::vector<LibUtilities::FieldDefinitionsSharedPtr> FieldDef;
1314 std::vector<std::vector<NekDouble>> FieldData;
1318 field_fld->Import(infile, FieldDef, FieldData);
1321 Vmath::Zero(pField->GetNcoeffs(), pField->UpdateCoeffs(), 1);
1323 for (
int i = 0; i < FieldDef.size(); ++i)
1326 for (
int j = 0; j < FieldData.size(); ++j)
1328 if (FieldDef[i]->
m_fields[j] == pFieldName)
1333 ASSERTL1(idx >= 0,
"Field " + pFieldName +
" not found.");
1335 pField->ExtractDataToCoeffs(FieldDef[i], FieldData[i],
1337 pField->UpdateCoeffs());
1339 pField->BwdTrans(pField->GetCoeffs(), pField->UpdatePhys());
1350 std::vector<std::string> &fieldStr,
1354 ASSERTL0(fieldStr.size() <= coeffs.size(),
1355 "length of fieldstr should be the same as pFields");
1357 std::vector<LibUtilities::FieldDefinitionsSharedPtr> FieldDef;
1358 std::vector<std::vector<NekDouble>> FieldData;
1362 field_fld->Import(infile, FieldDef, FieldData);
1365 for (
int j = 0; j < fieldStr.size(); ++j)
1368 for (
int i = 0; i < FieldDef.size(); ++i)
1370 m_fields[0]->ExtractDataToCoeffs(FieldDef[i], FieldData[i],
1371 fieldStr[j], coeffs[j]);
1386 m_fields[0]->EvalBasisNumModesMax());
1388 if (
m_session->GetComm()->GetSize() > 1)
1428 if (
m_session->DefinesSolverInfo(
"UpwindType"))
1431 m_session->GetSolverInfo(
"UpwindType"));
1434 if (
m_session->DefinesSolverInfo(
"AdvectionType"))
1436 std::string AdvectionType;
1437 AdvectionType =
m_session->GetSolverInfo(
"AdvectionType");
1439 s,
"Advection Type",
1454 "Mixed Continuous Galerkin and Discontinuous");
1457 if (
m_session->DefinesSolverInfo(
"DiffusionType"))
1459 std::string DiffusionType;
1460 DiffusionType =
m_session->GetSolverInfo(
"DiffusionType");
1462 s,
"Diffusion Type",
1474 ASSERTL0(
false,
"This function is not valid for the Base class");
1481 std::vector<std::string> &variables)
1483 boost::ignore_unused(fieldcoeffs, variables);
#define ASSERTL0(condition, msg)
#define ASSERTL1(condition, msg)
Assert Level 1 – Debugging which is used whether in FULLDEBUG or DEBUG compilation mode....
static GLOBAL_MAPPING_EXPORT MappingSharedPtr Load(const LibUtilities::SessionReaderSharedPtr &pSession, const Array< OneD, MultiRegions::ExpListSharedPtr > &pFields)
Return a pointer to the mapping, creating it on first call.
Describes the specification for a Basis.
static std::shared_ptr< FieldIO > CreateForFile(const LibUtilities::SessionReaderSharedPtr session, const std::string &filename)
Construct a FieldIO object for a given input filename.
static std::shared_ptr< FieldIO > CreateDefault(const LibUtilities::SessionReaderSharedPtr session)
Returns an object for the default FieldIO method.
Provides a generic Factory class.
Defines a specification for a set of points.
General purpose memory allocation routines with the ability to allocate from thread specific memory p...
static std::shared_ptr< DataType > AllocateSharedPtr(const Args &...args)
Allocate a shared pointer from the memory pool.
virtual SOLVER_UTILS_EXPORT void v_DoInitialise()
Virtual function for initialisation implementation.
int m_spacedim
Spatial dimension (>= expansion dim).
SOLVER_UTILS_EXPORT int GetTraceNpoints()
bool m_useFFT
Flag to determine if FFT is used for homogeneous transform.
SOLVER_UTILS_EXPORT int GetNvariables()
int m_expdim
Expansion dimension.
LibUtilities::FieldIOSharedPtr m_fld
Field input/output.
SOLVER_UTILS_EXPORT void SessionSummary(SummaryList &vSummary)
Write out a session summary.
virtual SOLVER_UTILS_EXPORT void v_SetInitialConditions(NekDouble initialtime=0.0, bool dumpInitialConditions=true, const int domain=0)
virtual SOLVER_UTILS_EXPORT ~EquationSystem()
Destructor.
SpatialDomains::MeshGraphSharedPtr m_graph
Pointer to graph defining mesh.
LibUtilities::CommSharedPtr m_comm
Communicator.
NekDouble m_timestep
Time step size.
SOLVER_UTILS_EXPORT void ImportFld(const std::string &infile, Array< OneD, MultiRegions::ExpListSharedPtr > &pFields)
Input field data from the given file.
NekDouble m_time
Current time of simulation.
bool m_multipleModes
Flag to determine if use multiple homogenenous modes are used.
Array< OneD, MultiRegions::ExpListSharedPtr > m_fields
Array holding all dependent variables.
virtual SOLVER_UTILS_EXPORT void v_TransPhysToCoeff()
Virtual function for transformation to coefficient space.
NekDouble m_fintime
Finish time of the simulation.
virtual SOLVER_UTILS_EXPORT NekDouble v_L2Error(unsigned int field, const Array< OneD, NekDouble > &exactsoln=NullNekDouble1DArray, bool Normalised=false)
Virtual function for the L_2 error computation between fields and a given exact solution.
virtual SOLVER_UTILS_EXPORT MultiRegions::ExpListSharedPtr v_GetPressure(void)
int m_npointsZ
number of points in Z direction (if homogeneous)
SOLVER_UTILS_EXPORT void Checkpoint_Output(const int n)
Write checkpoint file of m_fields.
std::map< std::string, SolverUtils::SessionFunctionSharedPtr > m_sessionFunctions
Map of known SessionFunctions.
NekDouble m_checktime
Time between checkpoints.
virtual SOLVER_UTILS_EXPORT void v_InitObject(bool DeclareFeld=true)
Initialisation object for EquationSystem.
virtual SOLVER_UTILS_EXPORT NekDouble v_LinfError(unsigned int field, const Array< OneD, NekDouble > &exactsoln=NullNekDouble1DArray)
Virtual function for the L_inf error computation between fields and a given exact solution.
virtual SOLVER_UTILS_EXPORT void v_EvaluateExactSolution(unsigned int field, Array< OneD, NekDouble > &outfield, const NekDouble time)
SOLVER_UTILS_EXPORT void WriteFld(const std::string &outname)
Write field data to the given filename.
SOLVER_UTILS_EXPORT void ImportFldToMultiDomains(const std::string &infile, Array< OneD, MultiRegions::ExpListSharedPtr > &pFields, const int ndomains)
Input field data from the given file to multiple domains.
NekDouble m_LhomY
physical length in Y direction (if homogeneous)
int m_npointsY
number of points in Y direction (if homogeneous)
std::string m_sessionName
Name of the session.
SOLVER_UTILS_EXPORT void Checkpoint_BaseFlow(const int n)
Write base flow file of m_fields.
bool m_specHP_dealiasing
Flag to determine if dealisising is usde for the Spectral/hp element discretisation.
LibUtilities::SessionReaderSharedPtr m_session
The session reader.
bool m_singleMode
Flag to determine if single homogeneous mode is used.
virtual SOLVER_UTILS_EXPORT bool v_NegatedOp()
Virtual function to identify if operator is negated in DoSolve.
SOLVER_UTILS_EXPORT void ZeroPhysFields()
int m_HomoDirec
number of homogenous directions
virtual SOLVER_UTILS_EXPORT void v_ExtraFldOutput(std::vector< Array< OneD, NekDouble >> &fieldcoeffs, std::vector< std::string > &variables)
Array< OneD, Array< OneD, NekDouble > > m_traceNormals
Array holding trace normals for DG simulations in the forwards direction.
enum HomogeneousType m_HomogeneousType
SOLVER_UTILS_EXPORT Array< OneD, NekDouble > ErrorExtraPoints(unsigned int field)
Compute error (L2 and L_inf) over an larger set of quadrature points return [L2 Linf].
SOLVER_UTILS_EXPORT int GetNpoints()
Array< OneD, bool > m_checkIfSystemSingular
Flag to indicate if the fields should be checked for singularity.
bool m_homogen_dealiasing
Flag to determine if dealiasing is used for homogeneous simulations.
SOLVER_UTILS_EXPORT int GetNcoeffs()
int m_nchk
Number of checkpoints written so far.
enum MultiRegions::ProjectionType m_projectionType
Type of projection; e.g continuous or discontinuous.
NekDouble m_TimeIncrementFactor
virtual SOLVER_UTILS_EXPORT void v_DoSolve()
Virtual function for solve implementation.
virtual SOLVER_UTILS_EXPORT Array< OneD, bool > v_GetSystemSingularChecks()
int m_steps
Number of steps to take.
int m_NumQuadPointsError
Number of Quadrature points used to work out the error.
Array< OneD, NekDouble > m_movingFrameTheta
Moving frame of reference angles with respect to the.
LibUtilities::FieldMetaDataMap m_fieldMetaDataMap
Map to identify relevant solver info to dump in output fields.
SOLVER_UTILS_EXPORT int GetNumExpModes()
SpatialDomains::BoundaryConditionsSharedPtr m_boundaryConditions
Pointer to boundary conditions object.
SOLVER_UTILS_EXPORT void SetBoundaryConditions(NekDouble time)
Evaluates the boundary conditions at the given time.
NekDouble m_LhomZ
physical length in Z direction (if homogeneous)
virtual SOLVER_UTILS_EXPORT void v_GenerateSummary(SummaryList &l)
Virtual function for generating summary information.
SOLVER_UTILS_EXPORT void ExtraFldOutput(std::vector< Array< OneD, NekDouble >> &fieldcoeffs, std::vector< std::string > &variables)
bool m_halfMode
Flag to determine if half homogeneous mode is used.
SOLVER_UTILS_EXPORT void FwdTransFields()
int m_checksteps
Number of steps between checkpoints.
SOLVER_UTILS_EXPORT SessionFunctionSharedPtr GetFunction(std::string name, const MultiRegions::ExpListSharedPtr &field=MultiRegions::NullExpListSharedPtr, bool cache=false)
Get a SessionFunction by name.
virtual SOLVER_UTILS_EXPORT void v_Output(void)
virtual SOLVER_UTILS_EXPORT void v_TransCoeffToPhys()
Virtual function for transformation to physical space.
GLOBAL_MAPPING_EXPORT typedef std::shared_ptr< Mapping > MappingSharedPtr
A shared pointer to a Mapping object.
std::vector< BasisKey > BasisKeyVector
Name for a vector of BasisKeys.
std::shared_ptr< FieldIO > FieldIOSharedPtr
void Import(const std::string &infilename, std::vector< FieldDefinitionsSharedPtr > &fielddefs, std::vector< std::vector< NekDouble >> &fielddata, FieldMetaDataMap &fieldinfomap, const Array< OneD, int > &ElementIDs)
This function allows for data to be imported from an FLD file when a session and/or communicator is n...
std::map< std::string, std::string > FieldMetaDataMap
std::shared_ptr< SessionReader > SessionReaderSharedPtr
std::shared_ptr< Equation > EquationSharedPtr
static FieldMetaDataMap NullFieldMetaDataMap
@ eFourierEvenlySpaced
1D Evenly-spaced points using Fourier Fit
@ eGaussLobattoLegendre
1D Gauss-Lobatto-Legendre quadrature points
@ eFourierSingleModeSpaced
1D Non Evenly-spaced points for Single Mode analysis
@ eModified_B
Principle Modified Functions .
@ eFourierSingleMode
Fourier ModifiedExpansion with just the first mode .
@ eModified_A
Principle Modified Functions .
@ eFourierHalfModeIm
Fourier Modified expansions with just the imaginary part of the first mode .
@ eFourierHalfModeRe
Fourier Modified expansions with just the real part of the first mode .
@ eFourier
Fourier Expansion .
@ eMixed_CG_Discontinuous
std::shared_ptr< ExpList > ExpListSharedPtr
Shared pointer to an ExpList object.
std::shared_ptr< ContField > ContFieldSharedPtr
AdvectionFactory & GetAdvectionFactory()
Gets the factory for initialising advection objects.
std::vector< std::pair< std::string, std::string > > SummaryList
DiffusionFactory & GetDiffusionFactory()
EquationSystemFactory & GetEquationSystemFactory()
void AddSummaryItem(SummaryList &l, const std::string &name, const std::string &value)
Adds a summary item to the summary info list.
std::shared_ptr< SessionFunction > SessionFunctionSharedPtr
std::shared_ptr< ExpansionInfoMap > ExpansionInfoMapShPtr
std::shared_ptr< MeshGraph > MeshGraphSharedPtr
std::map< int, ExpansionInfoShPtr > ExpansionInfoMap
The above copyright notice and this permission notice shall be included.
NekMatrix< NekMatrix< NekMatrix< NekDouble, StandardMatrixTag >, ScaledMatrixTag >, BlockMatrixTag > DNekScalBlkMat
NekMatrix< NekMatrix< NekDouble, StandardMatrixTag >, ScaledMatrixTag > DNekScalMat
void Zero(int n, T *x, const int incx)
Zero vector.
scalarT< T > sqrt(scalarT< T > in)
Used to lookup the create function in NekManager.