13 #include "Thyra_VectorStdOps.hpp"
15 #include "Tempus_IntegratorBasic.hpp"
16 #include "Tempus_WrapperModelEvaluatorPairIMEX_Basic.hpp"
17 #include "Tempus_StepperIMEX_RK.hpp"
19 #include "../TestModels/VanDerPol_IMEX_ExplicitModel.hpp"
20 #include "../TestModels/VanDerPol_IMEX_ImplicitModel.hpp"
21 #include "../TestUtils/Tempus_ConvergenceTestUtils.hpp"
26 namespace Tempus_Test {
28 using Teuchos::getParametersFromXmlFile;
32 using Teuchos::rcp_const_cast;
33 using Teuchos::sublist;
47 getParametersFromXmlFile(
"Tempus_IMEX_RK_VanDerPol.xml");
59 explicitModel, implicitModel));
63 stepper->setModel(model);
64 stepper->initialize();
70 timeStepControl->setInitIndex(tscPL.
get<
int>(
"Initial Time Index"));
71 timeStepControl->setInitTime(tscPL.
get<
double>(
"Initial Time"));
72 timeStepControl->setFinalTime(tscPL.
get<
double>(
"Final Time"));
73 timeStepControl->setInitTimeStep(dt);
74 timeStepControl->initialize();
77 auto inArgsIC = model->getNominalValues();
80 icState->setTime(timeStepControl->getInitTime());
81 icState->setIndex(timeStepControl->getInitIndex());
82 icState->setTimeStep(0.0);
83 icState->setOrder(stepper->getOrder());
88 solutionHistory->setName(
"Forward States");
90 solutionHistory->setStorageLimit(2);
91 solutionHistory->addState(icState);
95 Tempus::createIntegratorBasic<double>();
96 integrator->setStepper(stepper);
97 integrator->setTimeStepControl(timeStepControl);
98 integrator->setSolutionHistory(solutionHistory);
99 integrator->initialize();
102 bool integratorStatus = integrator->advanceTime();
106 double time = integrator->getTime();
107 double timeFinal = pl->
sublist(
"Default Integrator")
109 .
get<
double>(
"Final Time");
116 out <<
" Stepper = " << stepper->description() << std::endl;
117 out <<
" =========================" << std::endl;
118 out <<
" Computed solution: " << get_ele(*(x), 0) <<
" "
119 << get_ele(*(x), 1) << std::endl;
120 out <<
" =========================" << std::endl;
121 TEST_FLOATING_EQUALITY(get_ele(*(x), 0), 1.810210, 1.0e-4);
122 TEST_FLOATING_EQUALITY(get_ele(*(x), 1), -0.754602, 1.0e-4);
129 std::vector<std::string> stepperTypes;
130 stepperTypes.push_back(
"IMEX RK 1st order");
131 stepperTypes.push_back(
"SSP1_111");
132 stepperTypes.push_back(
"IMEX RK SSP2");
133 stepperTypes.push_back(
"SSP2_222");
134 stepperTypes.push_back(
"IMEX RK ARS 233");
135 stepperTypes.push_back(
"General IMEX RK");
136 stepperTypes.push_back(
"IMEX RK SSP3");
138 std::vector<double> stepperOrders;
139 stepperOrders.push_back(1.07964);
140 stepperOrders.push_back(1.07964);
141 stepperOrders.push_back(2.00408);
142 stepperOrders.push_back(2.76941);
143 stepperOrders.push_back(2.70655);
144 stepperOrders.push_back(2.00211);
145 stepperOrders.push_back(2.00211);
147 std::vector<double> stepperErrors;
148 stepperErrors.push_back(0.0046423);
149 stepperErrors.push_back(0.103569);
150 stepperErrors.push_back(0.0154534);
151 stepperErrors.push_back(0.000533759);
152 stepperErrors.push_back(0.000298908);
153 stepperErrors.push_back(0.0071546);
154 stepperErrors.push_back(0.0151202);
156 std::vector<double> stepperInitDt;
157 stepperInitDt.push_back(0.0125);
158 stepperInitDt.push_back(0.0125);
159 stepperInitDt.push_back(0.05);
160 stepperInitDt.push_back(0.05);
161 stepperInitDt.push_back(0.05);
162 stepperInitDt.push_back(0.05);
163 stepperInitDt.push_back(0.05);
169 std::vector<std::string>::size_type m;
170 for (m = 0; m != stepperTypes.size(); m++) {
171 std::string stepperType = stepperTypes[m];
172 std::string stepperName = stepperTypes[m];
173 std::replace(stepperName.begin(), stepperName.end(),
' ',
'_');
174 std::replace(stepperName.begin(), stepperName.end(),
'/',
'.');
177 std::vector<RCP<Thyra::VectorBase<double>>> solutions;
178 std::vector<RCP<Thyra::VectorBase<double>>> solutionsDot;
179 std::vector<double> StepSize;
180 std::vector<double> xErrorNorm;
181 std::vector<double> xDotErrorNorm;
183 const int nTimeStepSizes = 3;
184 double dt = stepperInitDt[m];
186 for (
int n = 0; n < nTimeStepSizes; n++) {
189 getParametersFromXmlFile(
"Tempus_IMEX_RK_VanDerPol.xml");
202 explicitModel, implicitModel));
206 if (stepperType ==
"General IMEX RK") {
208 pl->
sublist(
"Default Integrator")
209 .
set(
"Stepper Name",
"General IMEX RK");
212 pl->
sublist(
"Default Stepper").
set(
"Stepper Type", stepperType);
216 if (n == nTimeStepSizes - 1)
222 pl->
sublist(
"Default Integrator")
224 .
set(
"Initial Time Step", dt);
225 integrator = Tempus::createIntegratorBasic<double>(pl, model);
228 bool integratorStatus = integrator->advanceTime();
232 time = integrator->getTime();
233 double timeFinal = pl->sublist(
"Default Integrator")
234 .sublist(
"Time Step Control")
235 .get<
double>(
"Final Time");
236 double tol = 100.0 * std::numeric_limits<double>::epsilon();
240 StepSize.push_back(dt);
241 auto solution = Thyra::createMember(model->get_x_space());
242 Thyra::copy(*(integrator->getX()), solution.ptr());
243 solutions.push_back(solution);
244 auto solutionDot = Thyra::createMember(model->get_x_space());
245 Thyra::copy(*(integrator->getXDot()), solutionDot.ptr());
246 solutionsDot.push_back(solutionDot);
250 if ((n == 0) || (n == nTimeStepSizes - 1)) {
251 std::string fname =
"Tempus_" + stepperName +
"_VanDerPol-Ref.dat";
252 if (n == 0) fname =
"Tempus_" + stepperName +
"_VanDerPol.dat";
254 integrator->getSolutionHistory();
261 double xDotSlope = 0.0;
266 solutionsDot.clear();
268 writeOrderError(
"Tempus_" + stepperName +
"_VanDerPol-Error.dat", stepper,
269 StepSize, solutions, xErrorNorm, xSlope, solutionsDot,
270 xDotErrorNorm, xDotSlope, out);
Teuchos::RCP< SolutionState< Scalar > > createSolutionStateX(const Teuchos::RCP< Thyra::VectorBase< Scalar > > &x, const Teuchos::RCP< Thyra::VectorBase< Scalar > > &xdot=Teuchos::null, const Teuchos::RCP< Thyra::VectorBase< Scalar > > &xdotdot=Teuchos::null)
Nonmember constructor from non-const solution vectors, x.
van der Pol model formulated for IMEX-RK.
T & get(const std::string &name, T def_value)
ParameterList & set(std::string const &name, T const &value, std::string const &docString="", RCP< const ParameterEntryValidator > const &validator=null)
ModelEvaluator pair for implicit and explicit (IMEX) evaulations.
Implicit-Explicit Runge-Kutta (IMEX-RK) time stepper.
void writeOrderError(const std::string filename, Teuchos::RCP< Tempus::Stepper< Scalar >> stepper, std::vector< Scalar > &StepSize, std::vector< Teuchos::RCP< Thyra::VectorBase< Scalar >>> &solutions, std::vector< Scalar > &xErrorNorm, Scalar &xSlope, std::vector< Teuchos::RCP< Thyra::VectorBase< Scalar >>> &solutionsDot, std::vector< Scalar > &xDotErrorNorm, Scalar &xDotSlope, std::vector< Teuchos::RCP< Thyra::VectorBase< Scalar >>> &solutionsDotDot, std::vector< Scalar > &xDotDotErrorNorm, Scalar &xDotDotSlope, Teuchos::FancyOStream &out)
#define TEST_FLOATING_EQUALITY(v1, v2, tol)
void writeSolution(const std::string filename, Teuchos::RCP< const Tempus::SolutionHistory< Scalar >> solutionHistory)
TEUCHOS_UNIT_TEST(BackwardEuler, SinCos_ASA)
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
TimeStepControl manages the time step size. There several mechanisms that effect the time step size a...
SolutionHistory is basically a container of SolutionStates. SolutionHistory maintains a collection of...
Keep a fix number of states.
ParameterList & sublist(const std::string &name, bool mustAlreadyExist=false, const std::string &docString="")
van der Pol model formulated for IMEX.
#define TEUCHOS_ASSERT(assertion_test)
Solution state for integrators and steppers.