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"
34 namespace Tempus_Test {
36 using Teuchos::getParametersFromXmlFile;
40 using Teuchos::rcp_const_cast;
41 using Teuchos::rcp_dynamic_cast;
42 using Teuchos::sublist;
54 getParametersFromXmlFile(
"Tempus_Subcycling_SinCos.xml");
69 Tempus::createIntegratorBasic<double>(tempusPL, model);
73 integrator->getStepper()->getValidParameters();
75 bool pass = haveSameValuesSorted(*stepperPL, *defaultPL,
true);
78 out <<
"stepperPL -------------- \n"
79 << *stepperPL << std::endl;
80 out <<
"defaultPL -------------- \n"
81 << *defaultPL << std::endl;
89 Tempus::createIntegratorBasic<double>(model,
90 std::string(
"Forward Euler"));
94 integrator->getStepper()->getValidParameters();
96 bool pass = haveSameValuesSorted(*stepperPL, *defaultPL,
true);
99 out <<
"stepperPL -------------- \n"
100 << *stepperPL << std::endl;
101 out <<
"defaultPL -------------- \n"
102 << *defaultPL << std::endl;
121 stepper->setSubcyclingStepper(stepperFE);
123 stepper->setSubcyclingMinTimeStep(0.1);
124 stepper->setSubcyclingInitTimeStep(0.1);
125 stepper->setSubcyclingMaxTimeStep(0.1);
126 stepper->setSubcyclingMaxFailures(10);
127 stepper->setSubcyclingMaxConsecFailures(5);
128 stepper->setSubcyclingScreenOutputIndexInterval(1);
129 stepper->setSubcyclingIntegratorObserver(
131 stepper->setSubcyclingPrintDtChanges(
true);
136 stepper->setSubcyclingTimeStepControlStrategy(subStrategy);
140 timeStepControl->setInitIndex(0);
141 timeStepControl->setFinalIndex(10);
142 timeStepControl->setInitTime(0.0);
143 timeStepControl->setFinalTime(1.0);
144 timeStepControl->setInitTimeStep(dt);
148 strategy->initialize();
149 timeStepControl->setTimeStepControlStrategy(strategy);
151 timeStepControl->initialize();
154 auto inArgsIC = stepper->getModel()->getNominalValues();
157 icState->setTime(timeStepControl->getInitTime());
158 icState->setIndex(timeStepControl->getInitIndex());
159 icState->setTimeStep(0.0);
164 solutionHistory->setName(
"Forward States");
166 solutionHistory->setStorageLimit(2);
167 solutionHistory->addState(icState);
170 stepper->setInitialConditions(solutionHistory);
171 stepper->initialize();
175 Tempus::createIntegratorBasic<double>();
176 integrator->setStepper(stepper);
177 integrator->setTimeStepControl(timeStepControl);
178 integrator->setSolutionHistory(solutionHistory);
179 integrator->setScreenOutputIndexInterval(1);
181 integrator->initialize();
184 bool integratorStatus = integrator->advanceTime();
188 double time = integrator->getTime();
189 double timeFinal = 1.0;
195 model->getExactSolution(time).get_x();
199 Thyra::V_StVpStV(xdiff.
ptr(), 1.0, *x_exact, -1.0, *(x));
202 out <<
" Stepper = " << stepper->description() << std::endl;
203 out <<
" =========================" << std::endl;
204 out <<
" Exact solution : " << get_ele(*(x_exact), 0) <<
" "
205 << get_ele(*(x_exact), 1) << std::endl;
206 out <<
" Computed solution: " << get_ele(*(x), 0) <<
" "
207 << get_ele(*(x), 1) << std::endl;
208 out <<
" Difference : " << get_ele(*(xdiff), 0) <<
" "
209 << get_ele(*(xdiff), 1) << std::endl;
210 out <<
" =========================" << std::endl;
220 std::vector<RCP<Thyra::VectorBase<double>>> solutions;
221 std::vector<RCP<Thyra::VectorBase<double>>> solutionsDot;
222 std::vector<double> StepSize;
227 const int nTimeStepSizes = 2;
228 std::string output_file_string =
"Tempus_Subcycling_SinCos";
229 std::string output_file_name = output_file_string +
".dat";
230 std::string err_out_file_name = output_file_string +
"-Error.dat";
232 for (
int n = 0; n < nTimeStepSizes; n++) {
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();
276 icState->setTime(timeStepControl->getInitTime());
277 icState->setIndex(timeStepControl->getInitIndex());
278 icState->setTimeStep(0.0);
283 solutionHistory->setName(
"Forward States");
285 solutionHistory->setStorageLimit(2);
286 solutionHistory->addState(icState);
289 stepper->setInitialConditions(solutionHistory);
290 stepper->initialize();
293 integrator = Tempus::createIntegratorBasic<double>();
294 integrator->setStepper(stepper);
295 integrator->setTimeStepControl(timeStepControl);
296 integrator->setSolutionHistory(solutionHistory);
297 integrator->setScreenOutputIndexInterval(10);
299 integrator->initialize();
302 bool integratorStatus = integrator->advanceTime();
306 time = integrator->getTime();
307 double timeFinal = 1.0;
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()),
352 solutionExact.ptr());
353 solutions.push_back(solutionExact);
358 if (nTimeStepSizes > 1) {
360 double xDotSlope = 0.0;
361 std::vector<double> xErrorNorm;
362 std::vector<double> xDotErrorNorm;
366 solutions, xErrorNorm, xSlope, solutionsDot, xDotErrorNorm,
383 std::vector<RCP<Thyra::VectorBase<double>>> solutions;
384 std::vector<RCP<Thyra::VectorBase<double>>> solutionsDot;
385 std::vector<double> StepSize;
386 std::vector<double> xErrorNorm;
387 std::vector<double> xDotErrorNorm;
388 const int nTimeStepSizes = 4;
391 for (
int n = 0; n < nTimeStepSizes; n++) {
394 if (n == nTimeStepSizes - 1) dt /= 10.0;
399 tmpModel->getValidParameters());
400 pl->
set(
"Coeff epsilon", 0.1);
412 stepperFE->setUseFSAL(
false);
413 stepperFE->initialize();
414 stepperSC->setSubcyclingStepper(stepperFE);
416 stepperSC->setSubcyclingMinTimeStep(0.00001);
417 stepperSC->setSubcyclingInitTimeStep(dt / 10.0);
418 stepperSC->setSubcyclingMaxTimeStep(dt / 10.0);
419 stepperSC->setSubcyclingMaxFailures(10);
420 stepperSC->setSubcyclingMaxConsecFailures(5);
421 stepperSC->setSubcyclingScreenOutputIndexInterval(1);
427 strategySC->setMinEta(0.000001);
428 strategySC->setMaxEta(0.01);
429 strategySC->initialize();
430 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);
461 timeStepControl->initialize();
464 auto inArgsIC = stepper->getModel()->getNominalValues();
469 icState->setTime(timeStepControl->getInitTime());
470 icState->setIndex(timeStepControl->getInitIndex());
471 icState->setTimeStep(0.0);
472 icState->setOrder(stepper->getOrder());
477 solutionHistory->setName(
"Forward States");
480 solutionHistory->setStorageLimit(3);
481 solutionHistory->addState(icState);
484 stepperSC->setInitialConditions(solutionHistory);
485 stepper->initialize();
488 integrator = Tempus::createIntegratorBasic<double>();
489 integrator->setStepper(stepper);
490 integrator->setTimeStepControl(timeStepControl);
491 integrator->setSolutionHistory(solutionHistory);
492 integrator->setScreenOutputIndexInterval(10);
494 integrator->initialize();
497 bool integratorStatus = integrator->advanceTime();
501 time = integrator->getTime();
502 double timeFinal = 2.0;
503 double tol = 100.0 * std::numeric_limits<double>::epsilon();
507 StepSize.push_back(dt);
508 auto solution = Thyra::createMember(implicitModel->get_x_space());
509 Thyra::copy(*(integrator->getX()), solution.ptr());
510 solutions.push_back(solution);
511 auto solutionDot = Thyra::createMember(implicitModel->get_x_space());
512 Thyra::copy(*(integrator->getXDot()), solutionDot.ptr());
513 solutionsDot.push_back(solutionDot);
517 if ((n == 0) || (n == nTimeStepSizes - 1)) {
518 std::string fname =
"Tempus_Subcycling_VanDerPol-Ref.dat";
519 if (n == 0) fname =
"Tempus_Subcycling_VanDerPol.dat";
527 double xDotSlope = 0.0;
530 writeOrderError(
"Tempus_Subcycling_VanDerPol-Error.dat", stepper, StepSize,
531 solutions, xErrorNorm, xSlope, solutionsDot, xDotErrorNorm,
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.
StepControlStrategy class for TimeStepControl.
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)
Sine-Cosine model problem from Rythmos. This is a canonical Sine-Cosine differential equation with a...
void writeSolution(const std::string filename, Teuchos::RCP< const Tempus::SolutionHistory< Scalar >> solutionHistory)
OperatorSplit stepper loops through the Stepper list.
Teuchos::RCP< StepperForwardEuler< Scalar > > createStepperForwardEuler(const Teuchos::RCP< const Thyra::ModelEvaluator< Scalar > > &model, Teuchos::RCP< Teuchos::ParameterList > pl)
Nonmember constructor - ModelEvaluator and ParameterList.
TEUCHOS_UNIT_TEST(BackwardEuler, SinCos_ASA)
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
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.
StepControlStrategy class for TimeStepControl.
Teuchos::RCP< StepperBackwardEuler< Scalar > > createStepperBackwardEuler(const Teuchos::RCP< const Thyra::ModelEvaluator< Scalar > > &model, Teuchos::RCP< Teuchos::ParameterList > pl)
Nonmember constructor - ModelEvaluator and ParameterList.
van der Pol model formulated for IMEX.
IntegratorObserverSubcycling class for time integrators. This basic class has simple no-op functions...
Grow the history as needed.
Solution state for integrators and steppers.