13 #include "Thyra_VectorStdOps.hpp"
15 #include "Tempus_IntegratorBasic.hpp"
16 #include "Tempus_WrapperModelEvaluatorPairPartIMEX_Basic.hpp"
17 #include "Tempus_StepperIMEX_RK_Partition.hpp"
19 #include "../TestModels/VanDerPol_IMEX_ExplicitModel.hpp"
20 #include "../TestModels/VanDerPol_IMEXPart_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");
52 const bool useProductVector =
true;
61 const int numExplicitBlocks = 1;
62 const int parameterIndex = 4;
64 explicitModel, implicitModel, numExplicitBlocks, parameterIndex));
68 stepper->setModel(model);
69 stepper->initialize();
75 timeStepControl->setInitIndex(tscPL.
get<
int>(
"Initial Time Index"));
76 timeStepControl->setInitTime(tscPL.
get<
double>(
"Initial Time"));
77 timeStepControl->setFinalTime(tscPL.
get<
double>(
"Final Time"));
78 timeStepControl->setInitTimeStep(dt);
79 timeStepControl->initialize();
82 auto inArgsIC = model->getNominalValues();
85 icState->setTime(timeStepControl->getInitTime());
86 icState->setIndex(timeStepControl->getInitIndex());
87 icState->setTimeStep(0.0);
88 icState->setOrder(stepper->getOrder());
93 solutionHistory->setName(
"Forward States");
95 solutionHistory->setStorageLimit(2);
96 solutionHistory->addState(icState);
100 Tempus::createIntegratorBasic<double>();
101 integrator->setStepper(stepper);
102 integrator->setTimeStepControl(timeStepControl);
103 integrator->setSolutionHistory(solutionHistory);
105 integrator->initialize();
108 bool integratorStatus = integrator->advanceTime();
112 double time = integrator->getTime();
113 double timeFinal = pl->
sublist(
"Default Integrator")
115 .
get<
double>(
"Final Time");
122 out <<
" Stepper = " << stepper->description() << std::endl;
123 out <<
" =========================" << std::endl;
124 out <<
" Computed solution: " << get_ele(*(x), 0) <<
" "
125 << get_ele(*(x), 1) << std::endl;
126 out <<
" =========================" << std::endl;
127 TEST_FLOATING_EQUALITY(get_ele(*(x), 0), 1.810210, 1.0e-4);
128 TEST_FLOATING_EQUALITY(get_ele(*(x), 1), -0.754602, 1.0e-4);
135 std::vector<std::string> stepperTypes;
136 stepperTypes.push_back(
"Partitioned IMEX RK 1st order");
137 stepperTypes.push_back(
"Partitioned IMEX RK SSP2");
138 stepperTypes.push_back(
"Partitioned IMEX RK ARS 233");
139 stepperTypes.push_back(
"General Partitioned IMEX RK");
141 std::vector<double> stepperOrders;
142 stepperOrders.push_back(1.07964);
143 stepperOrders.push_back(2.00408);
144 stepperOrders.push_back(2.70655);
145 stepperOrders.push_back(2.00211);
147 std::vector<double> stepperErrors;
148 stepperErrors.push_back(0.0046423);
149 stepperErrors.push_back(0.0154534);
150 stepperErrors.push_back(0.000298908);
151 stepperErrors.push_back(0.0071546);
153 std::vector<double> stepperInitDt;
154 stepperInitDt.push_back(0.0125);
155 stepperInitDt.push_back(0.05);
156 stepperInitDt.push_back(0.05);
157 stepperInitDt.push_back(0.05);
159 std::vector<std::string>::size_type m;
160 for (m = 0; m != stepperTypes.size(); m++) {
161 std::string stepperType = stepperTypes[m];
162 std::string stepperName = stepperTypes[m];
163 std::replace(stepperName.begin(), stepperName.end(),
' ',
'_');
164 std::replace(stepperName.begin(), stepperName.end(),
'/',
'.');
167 std::vector<RCP<Thyra::VectorBase<double>>> solutions;
168 std::vector<RCP<Thyra::VectorBase<double>>> solutionsDot;
169 std::vector<double> StepSize;
170 std::vector<double> xErrorNorm;
171 std::vector<double> xDotErrorNorm;
173 const int nTimeStepSizes = 3;
174 double dt = stepperInitDt[m];
176 for (
int n = 0; n < nTimeStepSizes; n++) {
179 getParametersFromXmlFile(
"Tempus_IMEX_RK_VanDerPol.xml");
183 const bool useProductVector =
true;
184 auto explicitModel =
rcp(
192 const int numExplicitBlocks = 1;
193 const int parameterIndex = 4;
196 explicitModel, implicitModel, numExplicitBlocks, parameterIndex));
201 if (stepperType ==
"General Partitioned IMEX RK") {
203 pl->
sublist(
"Default Integrator")
204 .
set(
"Stepper Name",
"General IMEX RK");
207 pl->
sublist(
"Default Stepper").
set(
"Stepper Type", stepperType);
211 if (n == nTimeStepSizes - 1)
217 pl->
sublist(
"Default Integrator")
219 .
set(
"Initial Time Step", dt);
220 integrator = Tempus::createIntegratorBasic<double>(pl, model);
223 bool integratorStatus = integrator->advanceTime();
227 time = integrator->getTime();
228 double timeFinal = pl->sublist(
"Default Integrator")
229 .sublist(
"Time Step Control")
230 .get<
double>(
"Final Time");
231 double tol = 100.0 * std::numeric_limits<double>::epsilon();
235 StepSize.push_back(dt);
236 auto solution = Thyra::createMember(model->get_x_space());
237 Thyra::copy(*(integrator->getX()), solution.ptr());
238 solutions.push_back(solution);
239 auto solutionDot = Thyra::createMember(model->get_x_space());
240 Thyra::copy(*(integrator->getXDot()), solutionDot.ptr());
241 solutionsDot.push_back(solutionDot);
245 if ((n == 0) || (n == nTimeStepSizes - 1)) {
246 std::string fname =
"Tempus_" + stepperName +
"_VanDerPol-Ref.dat";
247 if (n == 0) fname =
"Tempus_" + stepperName +
"_VanDerPol.dat";
249 integrator->getSolutionHistory();
256 double xDotSlope = 0.0;
261 solutionsDot.clear();
263 writeOrderError(
"Tempus_" + stepperName +
"_VanDerPol-Error.dat", stepper,
264 StepSize, solutions, xErrorNorm, xSlope, solutionsDot,
265 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 the partitioned 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)
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)
Partitioned Implicit-Explicit Runge-Kutta (IMEX-RK) time stepper.
TEUCHOS_UNIT_TEST(BackwardEuler, SinCos_ASA)
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
ModelEvaluator pair for implicit and explicit (IMEX) evaulations.
static void summarize(Ptr< const Comm< int > > comm, std::ostream &out=std::cout, const bool alwaysWriteLocal=false, const bool writeGlobalStats=true, const bool writeZeroTimers=true, const ECounterSetOp setOp=Intersection, const std::string &filter="", const bool ignoreZeroTimers=false)
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.
Solution state for integrators and steppers.