39#include <boost/format.hpp>
50 "set-start-time",
"",
"Set the starting time of the simulation.");
54 "set-start-chknumber",
"",
55 "Set the starting number of the checkpoint file.");
92 m_session->MatchSolverInfo(
"DIFFUSIONADVANCEMENT",
"Explicit",
94 m_session->MatchSolverInfo(
"ADVECTIONADVANCEMENT",
"Explicit",
96 m_session->MatchSolverInfo(
"REACTIONADVANCEMENT",
"Explicit",
107 if (
m_session->DefinesSolverInfo(
"TimeIntegrationMethod") ||
138 "Final condition not unique: "
139 "fintime > 0.0 and Nsteps > 0");
142 "Timestep not unique: timestep > 0.0 & CFL > 0.0");
154 "m_CFLEnd >= m_cflSafetyFactor required");
169 x.first,
m_session, shared_from_this(), x.second)));
219 for (i = 0; i < nfields; ++i)
223 nvariables = nfields;
233 for (i = 0; i < nfields; ++i)
247 for (i = 0; i < nvariables; ++i)
263 bool doCheckTime =
false;
275 string abortFile =
"abort";
276 if (
m_session->DefinesSolverInfo(
"CheckAbortFile"))
278 abortFile =
m_session->GetSolverInfo(
"CheckAbortFile");
287 tmp_cflSafetyFactor =
292 if (!
m_comm->IsParallelInTime())
350 m_session->GetComm()->GetSpaceComm()->GetRank() == 0 &&
353 cpuPrevious = cpuTime;
358 for (i = 0; i < nvariables; ++i)
383 if (
m_comm->GetRank() == 0)
385 cout <<
"Reached Steady State to tolerance "
396 for (i = 0; i < nvariables; ++i)
407 if (
m_session->GetComm()->GetSpaceComm()->GetRank() == 0)
409 if (fs::exists(abortFile))
411 fs::remove(abortFile);
416 m_session->GetComm()->GetSpaceComm()->AllReduce(
419 ASSERTL0(!abortFlags[0],
"NaN found during time integration.");
429 totFilterTime += elapsed;
434 m_session->DefinesCmdLineArgument(
"verbose"))
437 s0 << x.first <<
":";
439 s1 << elapsed <<
"s";
441 s2 << elapsed / cpuPrevious * 100 <<
"%";
442 cout <<
"CPU time for filter " << setw(25) << left << s0.str()
443 << setw(12) << left << s1.str() << endl
444 <<
"\t Percentage of time integration: " << setw(10)
445 << left << s2.str() << endl;
454 ss << totFilterTime <<
"s";
455 cout <<
"Total filters CPU Time:\t\t\t " << setw(10) << left
466 vector<bool> transformed(nfields,
false);
467 for (i = 0; i < nfields; i++)
474 transformed[i] =
true;
481 for (i = 0; i < nfields; i++)
512 for (i = 0; i < nfields; i++)
524 for (i = 0; i < nvariables; ++i)
549 if (
m_comm->IsParallelInTime())
551 cout <<
"RANK " <<
m_session->GetComm()->GetTimeComm()->GetRank()
552 <<
" Steps: " << setw(8) << left << step + 1 <<
" "
553 <<
"Time: " << setw(12) << left <<
m_time;
557 cout <<
"Steps: " << setw(8) << left << step + 1 <<
" "
558 <<
"Time: " << setw(12) << left <<
m_time;
563 cout <<
" Time-step: " << setw(12) << left <<
m_timestep;
567 ss << cpuTime <<
"s";
568 cout <<
" CPU Time: " << setw(8) << left << ss.str() << endl;
574 if (
m_session->GetComm()->GetRank() == 0)
582 if (
m_session->GetSolverInfo(
"Driver") !=
"SteadyState" &&
583 m_session->GetSolverInfo(
"Driver") !=
"Parareal" &&
584 m_session->GetSolverInfo(
"Driver") !=
"PFASST")
586 cout <<
"Time-integration : " << intTime <<
"s" << endl;
615 if (
m_session->GetSolverInfo(
"EQTYPE") ==
616 "SteadyAdvectionDiffusionReaction")
645 outfile.open(
"solution1D.txt");
648 outfile << scientific << setw(17) << setprecision(16) << x[i] <<
" "
651 outfile << endl << endl;
660 if (
m_session->DefinesFunction(
"InitialConditions"))
662 for (
int i = 0; i <
m_fields.size(); ++i)
666 vType =
m_session->GetFunctionType(
"InitialConditions",
671 std::string filename =
m_session->GetFunctionFilename(
672 "InitialConditions",
m_session->GetVariable(i));
674 fs::path pfilename(filename);
677 if (fs::is_directory(pfilename))
679 fs::path metafile(
"Info.xml");
680 fs::path fullpath = pfilename / metafile;
693 time = boost::lexical_cast<NekDouble>(iter->second);
699 nchk = boost::lexical_cast<NekDouble>(iter->second);
707 if (
m_session->DefinesCmdLineArgument(
"set-start-time"))
709 time = boost::lexical_cast<NekDouble>(
710 m_session->GetCmdLineArgument<std::string>(
"set-start-time")
713 if (
m_session->DefinesCmdLineArgument(
"set-start-chknumber"))
715 nchk = boost::lexical_cast<int>(
716 m_session->GetCmdLineArgument<std::string>(
"set-start-chknumber"));
719 "Starting time and checkpoint number should be >= 0");
774 int phystot =
m_fields[0]->GetTotPoints();
775 int nvel = vel.size();
780 Vmath::Vmul(phystot, vel[0], 1, vel[0], 1, varcoeff, 1);
781 for (
int n = 1; n < nvel; ++n)
783 Vmath::Vvtvp(phystot, vel[n], 1, vel[n], 1, varcoeff, 1, varcoeff, 1);
787 for (
int i = 0; i <
m_fields[0]->GetNumElmts(); ++i)
789 int offset =
m_fields[0]->GetPhys_Offset(i);
790 int nq =
m_fields[0]->GetExp(i)->GetTotPoints();
795 for (
int n = 0; n <
m_fields[0]->GetExp(i)->GetNumBases(); ++n)
797 nmodes = max(nmodes,
m_fields[0]->GetExp(i)->GetBasisNumModes(n));
803 Vmath::Smul(nq, h, varcoeff + offset, 1, tmp = varcoeff + offset, 1);
817 const int nPoints =
m_fields[0]->GetTotPoints();
821 for (
int i = 0; i <
m_fields.size(); ++i)
828 if (
m_comm->GetRank() == 0)
831 m_session->GetSessionName() + std::string(
".resd");
833 m_errFile << setw(26) << left <<
"# Time";
835 m_errFile << setw(26) << left <<
"CPU_Time";
839 for (
int i = 0; i <
m_fields.size(); ++i)
856 const int nFields =
m_fields.size();
871 int stepp = step + 1;
876 for (
int i = 0; i < nFields; ++i)
890 cout <<
"-- Maximum L^2 residual: " << maxL2 << endl;
898 for (
int i = 0; i <
m_fields.size(); ++i)
914 const int nFields =
m_fields.size();
921 for (
int i = 0; i < nFields; ++i)
939 for (
int i = 0; i < nFields; ++i)
941 reference[i] = (reference[i] == 0) ? 1 : reference[i];
942 L2[i] =
sqrt(residual[i] / reference[i]);
955 if (&inarray != &outarray)
957 for (
int i = 0; i < inarray.size(); ++i)
#define ASSERTL0(condition, msg)
#define NEKERROR(type, msg)
Assert Level 0 – Fundamental assert which is used whether in FULLDEBUG, DEBUG or OPT compilation mode...
static std::shared_ptr< FieldIO > CreateForFile(const LibUtilities::SessionReaderSharedPtr session, const std::string &filename)
Construct a FieldIO object for a given input filename.
tBaseSharedPtr CreateInstance(tKey idKey, tParam... args)
Create an instance of the class referred to by idKey.
static std::string RegisterCmdLineArgument(const std::string &pName, const std::string &pShortName, const std::string &pDescription)
Registers a command-line argument with the session reader.
Binds a set of functions for use by time integration schemes.
NekDouble TimePerTest(unsigned int n)
Returns amount of seconds per iteration in a test with n iterations.
A base class for describing how to solve specific equations.
int m_spacedim
Spatial dimension (>= expansion dim).
LibUtilities::CommSharedPtr m_comm
Communicator.
NekDouble m_timestep
Time step size.
int m_infosteps
Number of time steps between outputting status information.
NekDouble m_time
Current time of simulation.
Array< OneD, MultiRegions::ExpListSharedPtr > m_fields
Array holding all dependent variables.
NekDouble m_fintime
Finish time of the simulation.
NekDouble m_lastCheckTime
SOLVER_UTILS_EXPORT void Checkpoint_Output(const int n)
Write checkpoint file of m_fields.
NekDouble m_checktime
Time between checkpoints.
virtual SOLVER_UTILS_EXPORT void v_InitObject(bool DeclareFeld=true)
Initialisation object for EquationSystem.
LibUtilities::SessionReaderSharedPtr m_session
The session reader.
int m_initialStep
Number of the step where the simulation should begin.
enum HomogeneousType m_HomogeneousType
SOLVER_UTILS_EXPORT int GetNpoints()
int m_nchk
Number of checkpoints written so far.
NekDouble m_TimeIncrementFactor
int m_steps
Number of steps to take.
LibUtilities::FieldMetaDataMap m_fieldMetaDataMap
Map to identify relevant solver info to dump in output fields.
SOLVER_UTILS_EXPORT void SetBoundaryConditions(NekDouble time)
Evaluates the boundary conditions at the given time.
virtual SOLVER_UTILS_EXPORT void v_GenerateSummary(SummaryList &l)
Virtual function for generating summary information.
int m_checksteps
Number of steps between checkpoints.
SOLVER_UTILS_EXPORT void SetInitialConditions(NekDouble initialtime=0.0, bool dumpInitialConditions=true, const int domain=0)
Initialise the data in the dependent fields.
SOLVER_UTILS_EXPORT int GetTotPoints()
LibUtilities::TimeIntegrationSchemeOperators m_ode
The time integration scheme operators to use.
NekDouble m_CFLGrowth
CFL growth rate.
bool m_explicitReaction
Indicates if explicit or implicit treatment of reaction is used.
Array< OneD, Array< OneD, NekDouble > > m_previousSolution
Storage for previous solution for steady-state check.
SOLVER_UTILS_EXPORT void v_DoInitialise(bool dumpInitialConditions=true) override
Sets up initial conditions.
virtual SOLVER_UTILS_EXPORT bool v_UpdateTimeStepCheck()
std::vector< std::pair< std::string, FilterSharedPtr > > m_filters
static std::string cmdSetStartTime
NekDouble m_cflSafetyFactor
CFL safety factor (comprise between 0 to 1).
LibUtilities::TimeIntegrationSchemeSharedPtr m_intScheme
Wrapper to the time integration scheme.
SOLVER_UTILS_EXPORT void DoDummyProjection(const Array< OneD, const Array< OneD, NekDouble > > &inarray, Array< OneD, Array< OneD, NekDouble > > &outarray, const NekDouble time)
Perform dummy projection.
NekDouble m_steadyStateTol
Tolerance to which steady state should be evaluated at.
int m_abortSteps
Number of steps between checks for abort conditions.
bool m_explicitAdvection
Indicates if explicit or implicit treatment of advection is used.
SOLVER_UTILS_EXPORT NekDouble MaxTimeStepEstimator()
Get the maximum timestep estimator for cfl control.
bool m_explicitDiffusion
Indicates if explicit or implicit treatment of diffusion is used.
virtual SOLVER_UTILS_EXPORT bool v_PostIntegrate(int step)
NekDouble m_CFLEnd
Maximun cfl in cfl growth.
SOLVER_UTILS_EXPORT ~UnsteadySystem() override
Destructor.
virtual SOLVER_UTILS_EXPORT void v_PrintSummaryStatistics(const NekDouble intTime)
Print Summary Statistics.
virtual SOLVER_UTILS_EXPORT void v_PrintStatusInformation(const int step, const NekDouble cpuTime)
Print Status Information.
bool CheckSteadyState(int step, const NekDouble &totCPUTime=0.0)
Calculate whether the system has reached a steady state by observing residuals to a user-defined tole...
virtual SOLVER_UTILS_EXPORT NekDouble v_GetTimeStep(const Array< OneD, const Array< OneD, NekDouble > > &inarray)
Return the timestep to be used for the next step in the time-marching loop.
void InitializeSteadyState()
SOLVER_UTILS_EXPORT void v_GenerateSummary(SummaryList &s) override
Print a summary of time stepping parameters.
virtual SOLVER_UTILS_EXPORT void v_SteadyStateResidual(int step, Array< OneD, NekDouble > &L2)
SOLVER_UTILS_EXPORT void AppendOutput1D(void)
Print the solution at each solution point in a txt file.
virtual SOLVER_UTILS_EXPORT bool v_PreIntegrate(int step)
SOLVER_UTILS_EXPORT void CheckForRestartTime(NekDouble &time, int &nchk)
SOLVER_UTILS_EXPORT UnsteadySystem(const LibUtilities::SessionReaderSharedPtr &pSession, const SpatialDomains::MeshGraphSharedPtr &pGraph)
Initialises UnsteadySystem class members.
virtual SOLVER_UTILS_EXPORT bool v_RequireFwdTrans()
SOLVER_UTILS_EXPORT void SteadyStateResidual(int step, Array< OneD, NekDouble > &L2)
static std::string cmdSetStartChkNum
std::vector< int > m_intVariables
SOLVER_UTILS_EXPORT void v_InitObject(bool DeclareField=true) override
Init object for UnsteadySystem class.
SOLVER_UTILS_EXPORT LibUtilities::TimeIntegrationSchemeOperators & GetTimeIntegrationSchemeOperators()
Returns the time integration scheme operators.
int m_steadyStateSteps
Check for steady state at step interval.
SOLVER_UTILS_EXPORT LibUtilities::TimeIntegrationSchemeSharedPtr & GetTimeIntegrationScheme()
Returns the time integration scheme.
int m_filtersInfosteps
Number of time steps between outputting filters information.
SOLVER_UTILS_EXPORT void SVVVarDiffCoeff(const Array< OneD, Array< OneD, NekDouble > > vel, StdRegions::VarCoeffMap &varCoeffMap)
Evaluate the SVV diffusion coefficient according to Moura's paper where it should proportional to h t...
bool m_homoInitialFwd
Flag to determine if simulation should start in homogeneous forward transformed state.
SOLVER_UTILS_EXPORT NekDouble GetTimeStep()
SOLVER_UTILS_EXPORT void v_DoSolve() override
Solves an unsteady problem.
TimeIntegrationSchemeFactory & GetTimeIntegrationSchemeFactory()
std::shared_ptr< FieldIO > FieldIOSharedPtr
static std::string PortablePath(const fs::path &path)
create portable path on different platforms for std::filesystem path.
std::shared_ptr< SessionReader > SessionReaderSharedPtr
static FieldMetaDataMap NullFieldMetaDataMap
std::shared_ptr< TimeIntegrationScheme > TimeIntegrationSchemeSharedPtr
static const NekDouble kNekZeroTol
std::vector< std::pair< std::string, std::string > > SummaryList
FilterFactory & GetFilterFactory()
void AddSummaryItem(SummaryList &l, const std::string &name, const std::string &value)
Adds a summary item to the summary info list.
std::shared_ptr< MeshGraph > MeshGraphSharedPtr
std::map< StdRegions::VarCoeffType, VarCoeffEntry > VarCoeffMap
std::vector< double > z(NPUPPER)
void Vsqrt(int n, const T *x, const int incx, T *y, const int incy)
sqrt y = sqrt(x)
void Vmul(int n, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
Multiply vector z = x*y.
void Vvtvp(int n, const T *w, const int incw, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
vvtvp (vector times vector plus vector): z = w*x + y
T Vsum(int n, const T *x, const int incx)
Subtract return sum(x)
void Smul(int n, const T alpha, const T *x, const int incx, T *y, const int incy)
Scalar multiply y = alpha*x.
int Nnan(int n, const T *x, const int incx)
Return number of NaN elements of x.
T Vmax(int n, const T *x, const int incx)
Return the maximum element in x – called vmax to avoid conflict with max.
void Vcopy(int n, const T *x, const int incx, T *y, const int incy)
void Vsub(int n, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
Subtract vector z = x-y.
scalarT< T > sqrt(scalarT< T > in)
std::vector< NekDouble > freeParams