9#include "Teuchos_UnitTestHarness.hpp"
10#include "Teuchos_XMLParameterListHelpers.hpp"
11#include "Teuchos_TimeMonitor.hpp"
13#include "Thyra_VectorStdOps.hpp"
15#include "Tempus_IntegratorBasic.hpp"
16#include "Tempus_IntegratorObserverSubcycling.hpp"
18#include "Tempus_StepperFactory.hpp"
19#include "Tempus_StepperForwardEuler.hpp"
20#include "Tempus_StepperBackwardEuler.hpp"
21#include "Tempus_StepperSubcycling.hpp"
22#include "Tempus_StepperOperatorSplit.hpp"
26#include "../TestModels/SinCosModel.hpp"
27#include "../TestModels/VanDerPol_IMEX_ExplicitModel.hpp"
28#include "../TestModels/VanDerPol_IMEX_ImplicitModel.hpp"
29#include "../TestUtils/Tempus_ConvergenceTestUtils.hpp"
38using Teuchos::rcp_const_cast;
39using Teuchos::rcp_dynamic_cast;
40using Teuchos::ParameterList;
41using Teuchos::sublist;
42using Teuchos::getParametersFromXmlFile;
54 RCP<ParameterList> pList =
55 getParametersFromXmlFile(
"Tempus_Subcycling_SinCos.xml");
62 RCP<ParameterList> scm_pl = sublist(pList,
"SinCosModel",
true);
65 RCP<ParameterList> tempusPL = sublist(pList,
"Tempus",
true);
69 RCP<Tempus::IntegratorBasic<double> > integrator =
70 Tempus::createIntegratorBasic<double>(tempusPL, model);
72 RCP<ParameterList> stepperPL = sublist(tempusPL,
"Demo Stepper",
true);
73 RCP<const ParameterList> defaultPL =
74 integrator->getStepper()->getValidParameters();
76 bool pass = haveSameValuesSorted(*stepperPL, *defaultPL,
true);
79 out <<
"stepperPL -------------- \n" << *stepperPL << std::endl;
80 out <<
"defaultPL -------------- \n" << *defaultPL << std::endl;
87 RCP<Tempus::IntegratorBasic<double> > integrator =
88 Tempus::createIntegratorBasic<double>(model, std::string(
"Forward Euler"));
90 RCP<ParameterList> stepperPL = sublist(tempusPL,
"Demo Stepper",
true);
91 RCP<const ParameterList> defaultPL =
92 integrator->getStepper()->getValidParameters();
94 bool pass = haveSameValuesSorted(*stepperPL, *defaultPL,
true);
96 std::cout << std::endl;
97 std::cout <<
"stepperPL -------------- \n" << *stepperPL << std::endl;
98 std::cout <<
"defaultPL -------------- \n" << *defaultPL << std::endl;
113 auto modelME = rcp_dynamic_cast<const Thyra::ModelEvaluator<double> > (model);
118 stepper->setSubcyclingStepper(stepperFE);
120 stepper->setSubcyclingMinTimeStep (0.1);
121 stepper->setSubcyclingInitTimeStep (0.1);
122 stepper->setSubcyclingMaxTimeStep (0.1);
123 stepper->setSubcyclingMaxFailures (10);
124 stepper->setSubcyclingMaxConsecFailures(5);
125 stepper->setSubcyclingScreenOutputIndexInterval(1);
126 stepper->setSubcyclingIntegratorObserver(
128 stepper->setSubcyclingPrintDtChanges (
true);
132 stepper->setSubcyclingTimeStepControlStrategy(subStrategy);
136 timeStepControl->setInitIndex(0);
137 timeStepControl->setFinalIndex(10);
138 timeStepControl->setInitTime (0.0);
139 timeStepControl->setFinalTime(1.0);
140 timeStepControl->setInitTimeStep(dt);
144 strategy->initialize();
145 timeStepControl->setTimeStepControlStrategy(strategy);
147 timeStepControl->initialize();
150 auto inArgsIC = stepper->getModel()->getNominalValues();
151 auto icSolution = rcp_const_cast<Thyra::VectorBase<double> > (inArgsIC.get_x());
153 icState->setTime (timeStepControl->getInitTime());
154 icState->setIndex (timeStepControl->getInitIndex());
155 icState->setTimeStep(0.0);
160 solutionHistory->setName(
"Forward States");
162 solutionHistory->setStorageLimit(2);
163 solutionHistory->addState(icState);
166 stepper->setInitialConditions(solutionHistory);
167 stepper->initialize();
170 RCP<Tempus::IntegratorBasic<double> > integrator =
171 Tempus::createIntegratorBasic<double>();
172 integrator->setStepper(stepper);
173 integrator->setTimeStepControl(timeStepControl);
174 integrator->setSolutionHistory(solutionHistory);
175 integrator->setScreenOutputIndexInterval(1);
177 integrator->initialize();
181 bool integratorStatus = integrator->advanceTime();
182 TEST_ASSERT(integratorStatus)
186 double time = integrator->getTime();
187 double timeFinal = 1.0;
188 TEST_FLOATING_EQUALITY(time, timeFinal, 1.0e-14);
191 RCP<Thyra::VectorBase<double> > x = integrator->getX();
192 RCP<const Thyra::VectorBase<double> > x_exact =
193 model->getExactSolution(time).get_x();
196 RCP<Thyra::VectorBase<double> > xdiff = x->clone_v();
197 Thyra::V_StVpStV(xdiff.ptr(), 1.0, *x_exact, -1.0, *(x));
200 std::cout <<
" Stepper = " << stepper->description() << std::endl;
201 std::cout <<
" =========================" << std::endl;
202 std::cout <<
" Exact solution : " << get_ele(*(x_exact), 0) <<
" "
203 << get_ele(*(x_exact), 1) << std::endl;
204 std::cout <<
" Computed solution: " << get_ele(*(x ), 0) <<
" "
205 << get_ele(*(x ), 1) << std::endl;
206 std::cout <<
" Difference : " << get_ele(*(xdiff ), 0) <<
" "
207 << get_ele(*(xdiff ), 1) << std::endl;
208 std::cout <<
" =========================" << std::endl;
209 TEST_FLOATING_EQUALITY(get_ele(*(x), 0), 0.882508, 1.0e-4 );
210 TEST_FLOATING_EQUALITY(get_ele(*(x), 1), 0.570790, 1.0e-4 );
218 RCP<Tempus::IntegratorBasic<double> > integrator;
219 std::vector<RCP<Thyra::VectorBase<double>>> solutions;
220 std::vector<RCP<Thyra::VectorBase<double>>> solutionsDot;
221 std::vector<double> StepSize;
226 const int nTimeStepSizes = 2;
227 std::string output_file_string =
"Tempus_Subcycling_SinCos";
228 std::string output_file_name = output_file_string +
".dat";
229 std::string err_out_file_name = output_file_string +
"-Error.dat";
231 for (
int n=0; n<nTimeStepSizes; n++) {
237 auto modelME=rcp_dynamic_cast<const Thyra::ModelEvaluator<double> > (model);
242 stepper->setSubcyclingStepper(stepperFE);
244 stepper->setSubcyclingMinTimeStep (dt/10.0);
245 stepper->setSubcyclingInitTimeStep (dt/10.0);
246 stepper->setSubcyclingMaxTimeStep (dt);
247 stepper->setSubcyclingMaxFailures (10);
248 stepper->setSubcyclingMaxConsecFailures(5);
249 stepper->setSubcyclingScreenOutputIndexInterval(1);
256 strategy->setMinEta(0.02);
257 strategy->setMaxEta(0.04);
258 strategy->initialize();
259 stepper->setSubcyclingTimeStepControlStrategy(strategy);
263 timeStepControl->setInitIndex(0);
264 timeStepControl->setInitTime (0.0);
265 timeStepControl->setFinalTime(1.0);
266 timeStepControl->setMinTimeStep (dt);
267 timeStepControl->setInitTimeStep(dt);
268 timeStepControl->setMaxTimeStep (dt);
269 timeStepControl->initialize();
272 auto inArgsIC = stepper->getModel()->getNominalValues();
273 auto icSolution = rcp_const_cast<Thyra::VectorBase<double> > (inArgsIC.get_x());
275 icState->setTime (timeStepControl->getInitTime());
276 icState->setIndex (timeStepControl->getInitIndex());
277 icState->setTimeStep(0.0);
282 solutionHistory->setName(
"Forward States");
284 solutionHistory->setStorageLimit(2);
285 solutionHistory->addState(icState);
288 stepper->setInitialConditions(solutionHistory);
289 stepper->initialize();
292 integrator = Tempus::createIntegratorBasic<double>();
293 integrator->setStepper(stepper);
294 integrator->setTimeStepControl(timeStepControl);
295 integrator->setSolutionHistory(solutionHistory);
296 integrator->setScreenOutputIndexInterval(10);
298 integrator->initialize();
302 bool integratorStatus = integrator->advanceTime();
303 TEST_ASSERT(integratorStatus)
306 time = integrator->getTime();
307 double timeFinal = 1.0;
308 TEST_FLOATING_EQUALITY(time, timeFinal, 1.0e-14);
311 RCP<Thyra::VectorBase<double> > x = integrator->getX();
312 RCP<const Thyra::VectorBase<double> > x_exact =
313 model->getExactSolution(time).get_x();
344 StepSize.push_back(dt);
345 auto solution = Thyra::createMember(model->get_x_space());
346 Thyra::copy(*(integrator->getX()),solution.ptr());
347 solutions.push_back(solution);
348 if (n == nTimeStepSizes-1) {
349 StepSize.push_back(0.0);
350 auto solutionExact = Thyra::createMember(model->get_x_space());
351 Thyra::copy(*(model->getExactSolution(time).get_x()),solutionExact.ptr());
352 solutions.push_back(solutionExact);
357 if (nTimeStepSizes > 1) {
359 double xDotSlope = 0.0;
360 std::vector<double> xErrorNorm;
361 std::vector<double> xDotErrorNorm;
362 RCP<Tempus::Stepper<double> > stepper = integrator->getStepper();
366 solutions, xErrorNorm, xSlope,
367 solutionsDot, xDotErrorNorm, xDotSlope);
369 TEST_FLOATING_EQUALITY( xSlope, 1.00137, 0.01 );
371 TEST_FLOATING_EQUALITY( xErrorNorm[0], 0.00387948, 1.0e-4 );
375 Teuchos::TimeMonitor::summarize();
383 RCP<Tempus::IntegratorBasic<double> > integrator;
384 std::vector<RCP<Thyra::VectorBase<double>>> solutions;
385 std::vector<RCP<Thyra::VectorBase<double>>> solutionsDot;
386 std::vector<double> StepSize;
387 std::vector<double> xErrorNorm;
388 std::vector<double> xDotErrorNorm;
389 const int nTimeStepSizes = 4;
392 for (
int n=0; n<nTimeStepSizes; n++) {
396 if (n == nTimeStepSizes-1) dt /= 10.0;
400 auto pl = Teuchos::rcp_const_cast<Teuchos::ParameterList> (
401 tmpModel->getValidParameters());
402 pl->set(
"Coeff epsilon", 0.1);
403 RCP<const Thyra::ModelEvaluator<double> > explicitModel =
405 RCP<const Thyra::ModelEvaluator<double> > implicitModel =
413 stepperFE->setUseFSAL(
false);
414 stepperFE->initialize();
415 stepperSC->setSubcyclingStepper(stepperFE);
417 stepperSC->setSubcyclingMinTimeStep (0.00001);
418 stepperSC->setSubcyclingInitTimeStep (dt/10.0);
419 stepperSC->setSubcyclingMaxTimeStep (dt/10.0);
420 stepperSC->setSubcyclingMaxFailures (10);
421 stepperSC->setSubcyclingMaxConsecFailures(5);
422 stepperSC->setSubcyclingScreenOutputIndexInterval(1);
428 strategySC->setMinEta(0.000001);
429 strategySC->setMaxEta(0.01);
430 strategySC->initialize();
431 stepperSC->setSubcyclingTimeStepControlStrategy(strategySC);
438 stepper->addStepper(stepperSC);
439 stepper->addStepper(stepperBE);
443 timeStepControl->setInitIndex(0);
444 timeStepControl->setInitTime (0.0);
446 timeStepControl->setFinalTime(2.0);
447 timeStepControl->setMinTimeStep (0.000001);
448 timeStepControl->setInitTimeStep(dt);
449 timeStepControl->setMaxTimeStep (dt);
460 timeStepControl->initialize();
463 auto inArgsIC = stepper->getModel()->getNominalValues();
464 auto icX = rcp_const_cast<Thyra::VectorBase<double> > (inArgsIC.get_x());
465 auto icXDot = rcp_const_cast<Thyra::VectorBase<double> > (inArgsIC.get_x_dot());
467 icState->setTime (timeStepControl->getInitTime());
468 icState->setIndex (timeStepControl->getInitIndex());
469 icState->setTimeStep(0.0);
470 icState->setOrder (stepper->getOrder());
475 solutionHistory->setName(
"Forward States");
478 solutionHistory->setStorageLimit(3);
479 solutionHistory->addState(icState);
482 stepperSC->setInitialConditions(solutionHistory);
483 stepper->initialize();
486 integrator = Tempus::createIntegratorBasic<double>();
487 integrator->setStepper(stepper);
488 integrator->setTimeStepControl(timeStepControl);
489 integrator->setSolutionHistory(solutionHistory);
490 integrator->setScreenOutputIndexInterval(10);
492 integrator->initialize();
496 bool integratorStatus = integrator->advanceTime();
497 TEST_ASSERT(integratorStatus)
500 time = integrator->getTime();
501 double timeFinal = 2.0;
502 double tol = 100.0 * std::numeric_limits<double>::epsilon();
503 TEST_FLOATING_EQUALITY(time, timeFinal, tol);
506 StepSize.push_back(dt);
507 auto solution = Thyra::createMember(implicitModel->get_x_space());
508 Thyra::copy(*(integrator->getX()),solution.ptr());
509 solutions.push_back(solution);
510 auto solutionDot = Thyra::createMember(implicitModel->get_x_space());
511 Thyra::copy(*(integrator->getXDot()),solutionDot.ptr());
512 solutionsDot.push_back(solutionDot);
516 if ((n == 0) || (n == nTimeStepSizes-1)) {
517 std::string fname =
"Tempus_Subcycling_VanDerPol-Ref.dat";
518 if (n == 0) fname =
"Tempus_Subcycling_VanDerPol.dat";
526 double xDotSlope = 0.0;
527 RCP<Tempus::Stepper<double> > stepper = integrator->getStepper();
531 solutions, xErrorNorm, xSlope,
532 solutionsDot, xDotErrorNorm, xDotSlope);
534 TEST_FLOATING_EQUALITY( xSlope, 1.25708, 0.05 );
535 TEST_FLOATING_EQUALITY( xDotSlope, 1.20230, 0.05 );
536 TEST_FLOATING_EQUALITY( xErrorNorm[0], 0.37156, 1.0e-4 );
537 TEST_FLOATING_EQUALITY( xDotErrorNorm[0], 3.11651, 1.0e-4 );
539 Teuchos::TimeMonitor::summarize();
IntegratorObserverSubcycling class for time integrators. This basic class has simple no-op functions,...
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...
OperatorSplit stepper loops through the Stepper list.
StepControlStrategy class for TimeStepControl.
StepControlStrategy class for TimeStepControl.
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 formulated for IMEX.
van der Pol model formulated for IMEX-RK.
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.
@ STORAGE_TYPE_UNLIMITED
Grow the history as needed.
Teuchos::RCP< StepperBackwardEuler< Scalar > > createStepperBackwardEuler(const Teuchos::RCP< const Thyra::ModelEvaluator< Scalar > > &model, Teuchos::RCP< Teuchos::ParameterList > pl)
Nonmember constructor - ModelEvaluator and ParameterList.
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.
Teuchos::RCP< StepperForwardEuler< Scalar > > createStepperForwardEuler(const Teuchos::RCP< const Thyra::ModelEvaluator< Scalar > > &model, Teuchos::RCP< Teuchos::ParameterList > pl)
Nonmember constructor - ModelEvaluator and ParameterList.