9#include "Teuchos_UnitTestHarness.hpp"
10#include "Teuchos_XMLParameterListHelpers.hpp"
11#include "Teuchos_TimeMonitor.hpp"
12#include "Teuchos_DefaultComm.hpp"
14#include "Tempus_config.hpp"
15#include "Tempus_IntegratorBasic.hpp"
16#include "Tempus_StepperTrapezoidal.hpp"
18#include "../TestModels/SinCosModel.hpp"
19#include "../TestModels/VanDerPolModel.hpp"
20#include "../TestUtils/Tempus_ConvergenceTestUtils.hpp"
22#include "Stratimikos_DefaultLinearSolverBuilder.hpp"
23#include "Thyra_LinearOpWithSolveFactoryHelpers.hpp"
25#ifdef Tempus_ENABLE_MPI
26#include "Epetra_MpiComm.h"
28#include "Epetra_SerialComm.h"
40using Teuchos::rcp_const_cast;
41using Teuchos::ParameterList;
42using Teuchos::sublist;
43using Teuchos::getParametersFromXmlFile;
55 RCP<ParameterList> pList =
56 getParametersFromXmlFile(
"Tempus_Trapezoidal_SinCos.xml");
59 RCP<ParameterList> scm_pl = sublist(pList,
"SinCosModel",
true);
62 RCP<ParameterList> tempusPL = sublist(pList,
"Tempus",
true);
66 RCP<Tempus::IntegratorBasic<double> > integrator =
67 Tempus::createIntegratorBasic<double>(tempusPL, model);
69 RCP<ParameterList> stepperPL = sublist(tempusPL,
"Default Stepper",
true);
71 RCP<const ParameterList> defaultPL =
72 integrator->getStepper()->getValidParameters();
73 bool pass = haveSameValuesSorted(*stepperPL, *defaultPL,
true);
76 out <<
"stepperPL -------------- \n" << *stepperPL << std::endl;
77 out <<
"defaultPL -------------- \n" << *defaultPL << std::endl;
84 RCP<Tempus::IntegratorBasic<double> > integrator =
85 Tempus::createIntegratorBasic<double>(model, std::string(
"Trapezoidal Method"));
87 RCP<ParameterList> stepperPL = sublist(tempusPL,
"Default Stepper",
true);
88 RCP<const ParameterList> defaultPL =
89 integrator->getStepper()->getValidParameters();
91 bool pass = haveSameValuesSorted(*stepperPL, *defaultPL,
true);
94 out <<
"stepperPL -------------- \n" << *stepperPL << std::endl;
95 out <<
"defaultPL -------------- \n" << *defaultPL << std::endl;
107 std::vector<std::string> options;
108 options.push_back(
"Default Parameters");
109 options.push_back(
"ICConsistency and Check");
111 for(
const auto& option: options) {
114 RCP<ParameterList> pList =
115 getParametersFromXmlFile(
"Tempus_Trapezoidal_SinCos.xml");
116 RCP<ParameterList> pl = sublist(pList,
"Tempus",
true);
119 RCP<ParameterList> scm_pl = sublist(pList,
"SinCosModel",
true);
138 stepper->setModel(model);
139 if ( option ==
"ICConsistency and Check") {
140 stepper->setICConsistency(
"Consistent");
141 stepper->setICConsistencyCheck(
true);
143 stepper->initialize();
147 ParameterList tscPL = pl->sublist(
"Default Integrator")
148 .sublist(
"Time Step Control");
149 timeStepControl->setInitIndex(tscPL.get<
int> (
"Initial Time Index"));
150 timeStepControl->setInitTime (tscPL.get<
double>(
"Initial Time"));
151 timeStepControl->setFinalTime(tscPL.get<
double>(
"Final Time"));
152 timeStepControl->setInitTimeStep(dt);
153 timeStepControl->initialize();
156 auto inArgsIC = model->getNominalValues();
157 auto icSoln = rcp_const_cast<Thyra::VectorBase<double> > (inArgsIC.get_x());
159 rcp_const_cast<Thyra::VectorBase<double> > (inArgsIC.get_x_dot());
161 icState->setTime (timeStepControl->getInitTime());
162 icState->setIndex (timeStepControl->getInitIndex());
163 icState->setTimeStep(0.0);
164 icState->setOrder (stepper->getOrder());
169 solutionHistory->setName(
"Forward States");
171 solutionHistory->setStorageLimit(2);
172 solutionHistory->addState(icState);
175 RCP<Tempus::IntegratorBasic<double> > integrator =
176 Tempus::createIntegratorBasic<double>();
177 integrator->setStepper(stepper);
178 integrator->setTimeStepControl(timeStepControl);
179 integrator->setSolutionHistory(solutionHistory);
181 integrator->initialize();
185 bool integratorStatus = integrator->advanceTime();
186 TEST_ASSERT(integratorStatus)
190 double time = integrator->getTime();
191 double timeFinal =pl->sublist(
"Default Integrator")
192 .sublist(
"Time Step Control").get<
double>(
"Final Time");
193 TEST_FLOATING_EQUALITY(time, timeFinal, 1.0e-14);
196 RCP<Thyra::VectorBase<double> > x = integrator->getX();
197 RCP<const Thyra::VectorBase<double> > x_exact =
198 model->getExactSolution(time).get_x();
201 RCP<Thyra::VectorBase<double> > xdiff = x->clone_v();
202 Thyra::V_StVpStV(xdiff.ptr(), 1.0, *x_exact, -1.0, *(x));
205 out <<
" Stepper = " << stepper->description()
206 <<
"\n with " << option << std::endl;
207 out <<
" =========================" << std::endl;
208 out <<
" Exact solution : " << get_ele(*(x_exact), 0) <<
" "
209 << get_ele(*(x_exact), 1) << std::endl;
210 out <<
" Computed solution: " << get_ele(*(x ), 0) <<
" "
211 << get_ele(*(x ), 1) << std::endl;
212 out <<
" Difference : " << get_ele(*(xdiff ), 0) <<
" "
213 << get_ele(*(xdiff ), 1) << std::endl;
214 out <<
" =========================" << std::endl;
215 TEST_FLOATING_EQUALITY(get_ele(*(x), 0), 0.841021, 1.0e-4 );
216 TEST_FLOATING_EQUALITY(get_ele(*(x), 1), 0.541002, 1.0e-4 );
225 RCP<Tempus::IntegratorBasic<double> > integrator;
226 std::vector<RCP<Thyra::VectorBase<double>>> solutions;
227 std::vector<RCP<Thyra::VectorBase<double>>> solutionsDot;
228 std::vector<double> StepSize;
229 std::vector<double> xErrorNorm;
230 std::vector<double> xDotErrorNorm;
231 const int nTimeStepSizes = 7;
234 for (
int n=0; n<nTimeStepSizes; n++) {
237 RCP<ParameterList> pList =
238 getParametersFromXmlFile(
"Tempus_Trapezoidal_SinCos.xml");
245 RCP<ParameterList> scm_pl = sublist(pList,
"SinCosModel",
true);
252 RCP<ParameterList> pl = sublist(pList,
"Tempus",
true);
253 pl->sublist(
"Default Integrator")
254 .sublist(
"Time Step Control").set(
"Initial Time Step", dt);
255 integrator = Tempus::createIntegratorBasic<double>(pl, model);
262 RCP<Thyra::VectorBase<double> > x0 =
263 model->getNominalValues().get_x()->clone_v();
264 RCP<Thyra::VectorBase<double> > xdot0 =
265 model->getNominalValues().get_x_dot()->clone_v();
266 integrator->initializeSolutionHistory(0.0, x0, xdot0);
270 bool integratorStatus = integrator->advanceTime();
271 TEST_ASSERT(integratorStatus)
274 time = integrator->getTime();
275 double timeFinal =pl->sublist(
"Default Integrator")
276 .sublist(
"Time Step Control").get<
double>(
"Final Time");
277 TEST_FLOATING_EQUALITY(time, timeFinal, 1.0e-14);
281 RCP<const SolutionHistory<double> > solutionHistory =
282 integrator->getSolutionHistory();
283 writeSolution(
"Tempus_Trapezoidal_SinCos.dat", solutionHistory);
286 for (
int i=0; i<solutionHistory->getNumStates(); i++) {
287 double time_i = (*solutionHistory)[i]->getTime();
290 model->getExactSolution(time_i).get_x()),
292 model->getExactSolution(time_i).get_x_dot()));
293 state->setTime((*solutionHistory)[i]->getTime());
294 solnHistExact->addState(state);
296 writeSolution(
"Tempus_Trapezoidal_SinCos-Ref.dat", solnHistExact);
300 StepSize.push_back(dt);
301 auto solution = Thyra::createMember(model->get_x_space());
302 Thyra::copy(*(integrator->getX()),solution.ptr());
303 solutions.push_back(solution);
304 auto solutionDot = Thyra::createMember(model->get_x_space());
305 Thyra::copy(*(integrator->getXDot()),solutionDot.ptr());
306 solutionsDot.push_back(solutionDot);
307 if (n == nTimeStepSizes-1) {
308 StepSize.push_back(0.0);
309 auto solutionExact = Thyra::createMember(model->get_x_space());
310 Thyra::copy(*(model->getExactSolution(time).get_x()),solutionExact.ptr());
311 solutions.push_back(solutionExact);
312 auto solutionDotExact = Thyra::createMember(model->get_x_space());
313 Thyra::copy(*(model->getExactSolution(time).get_x_dot()),
314 solutionDotExact.ptr());
315 solutionsDot.push_back(solutionDotExact);
320 double xDotSlope = 0.0;
321 RCP<Tempus::Stepper<double> > stepper = integrator->getStepper();
322 double order = stepper->getOrder();
325 solutions, xErrorNorm, xSlope,
326 solutionsDot, xDotErrorNorm, xDotSlope);
328 TEST_FLOATING_EQUALITY( xSlope, order, 0.01 );
329 TEST_FLOATING_EQUALITY( xErrorNorm[0], 0.000832086, 1.0e-4 );
330 TEST_FLOATING_EQUALITY( xDotSlope, order, 0.01 );
331 TEST_FLOATING_EQUALITY( xDotErrorNorm[0], 0.000832086, 1.0e-4 );
333 Teuchos::TimeMonitor::summarize();
341 RCP<Tempus::IntegratorBasic<double> > integrator;
342 std::vector<RCP<Thyra::VectorBase<double>>> solutions;
343 std::vector<RCP<Thyra::VectorBase<double>>> solutionsDot;
344 std::vector<double> StepSize;
345 std::vector<double> xErrorNorm;
346 std::vector<double> xDotErrorNorm;
347 const int nTimeStepSizes = 4;
350 for (
int n=0; n<nTimeStepSizes; n++) {
353 RCP<ParameterList> pList =
354 getParametersFromXmlFile(
"Tempus_Trapezoidal_VanDerPol.xml");
357 RCP<ParameterList> vdpm_pl = sublist(pList,
"VanDerPolModel",
true);
362 if (n == nTimeStepSizes-1) dt /= 10.0;
365 RCP<ParameterList> pl = sublist(pList,
"Tempus",
true);
366 pl->sublist(
"Demo Integrator")
367 .sublist(
"Time Step Control").set(
"Initial Time Step", dt);
368 integrator = Tempus::createIntegratorBasic<double>(pl, model);
371 bool integratorStatus = integrator->advanceTime();
372 TEST_ASSERT(integratorStatus)
375 time = integrator->getTime();
376 double timeFinal =pl->sublist(
"Demo Integrator")
377 .sublist(
"Time Step Control").get<
double>(
"Final Time");
378 double tol = 100.0 * std::numeric_limits<double>::epsilon();
379 TEST_FLOATING_EQUALITY(time, timeFinal, tol);
382 StepSize.push_back(dt);
383 auto solution = Thyra::createMember(model->get_x_space());
384 Thyra::copy(*(integrator->getX()),solution.ptr());
385 solutions.push_back(solution);
386 auto solutionDot = Thyra::createMember(model->get_x_space());
387 Thyra::copy(*(integrator->getXDot()),solutionDot.ptr());
388 solutionsDot.push_back(solutionDot);
392 if ((n == 0) || (n == nTimeStepSizes-1)) {
393 std::string fname =
"Tempus_Trapezoidal_VanDerPol-Ref.dat";
394 if (n == 0) fname =
"Tempus_Trapezoidal_VanDerPol.dat";
395 RCP<const SolutionHistory<double> > solutionHistory =
396 integrator->getSolutionHistory();
402 double xDotSlope = 0.0;
403 RCP<Tempus::Stepper<double> > stepper = integrator->getStepper();
404 double order = stepper->getOrder();
407 solutions, xErrorNorm, xSlope,
408 solutionsDot, xDotErrorNorm, xDotSlope);
410 TEST_FLOATING_EQUALITY( xSlope, order, 0.10 );
411 TEST_FLOATING_EQUALITY( xDotSlope, order, 0.10 );
412 TEST_FLOATING_EQUALITY( xErrorNorm[0], 0.00085855, 1.0e-4 );
413 TEST_FLOATING_EQUALITY( xDotErrorNorm[0], 0.00120695, 1.0e-4 );
415 Teuchos::TimeMonitor::summarize();
SolutionHistory is basically a container of SolutionStates. SolutionHistory maintains a collection of...
Solution state for integrators and steppers. SolutionState contains the metadata for solutions and th...
Trapezoidal method time stepper.
TimeStepControl manages the time step size. There several mechanisms that effect the time step size a...
Sine-Cosine model problem from Rythmos. This is a canonical Sine-Cosine differential equation.
van der Pol model problem for nonlinear electrical circuit.
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)
void writeSolution(const std::string filename, Teuchos::RCP< const Tempus::SolutionHistory< Scalar > > solutionHistory)
TEUCHOS_UNIT_TEST(BackwardEuler, SinCos_ASA)
@ STORAGE_TYPE_STATIC
Keep a fix number of states.
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.