9#include "Teuchos_UnitTestHarness.hpp"
10#include "Teuchos_XMLParameterListHelpers.hpp"
11#include "Teuchos_TimeMonitor.hpp"
13#include "Thyra_VectorStdOps.hpp"
15#include "Tempus_config.hpp"
16#include "Tempus_IntegratorBasic.hpp"
17#include "Tempus_StepperLeapfrog.hpp"
19#include "../TestModels/HarmonicOscillatorModel.hpp"
20#include "../TestUtils/Tempus_ConvergenceTestUtils.hpp"
23#ifdef Tempus_ENABLE_MPI
24#include "Epetra_MpiComm.h"
26#include "Epetra_SerialComm.h"
38using Teuchos::rcp_const_cast;
39using Teuchos::ParameterList;
40using Teuchos::sublist;
41using Teuchos::getParametersFromXmlFile;
54 std::vector<std::string> options;
55 options.push_back(
"Default Parameters");
56 options.push_back(
"ICConsistency and Check");
58 for(
const auto& option: options) {
61 RCP<ParameterList> pList =
62 getParametersFromXmlFile(
"Tempus_Leapfrog_SinCos.xml");
63 RCP<ParameterList> pl = sublist(pList,
"Tempus",
true);
66 RCP<ParameterList> hom_pl = sublist(pList,
"HarmonicOscillatorModel",
true);
71 stepper->setModel(model);
72 if ( option ==
"ICConsistency and Check") {
73 stepper->setICConsistency(
"Consistent");
74 stepper->setICConsistencyCheck(
true);
76 stepper->initialize();
80 ParameterList tscPL = pl->sublist(
"Default Integrator")
81 .sublist(
"Time Step Control");
82 timeStepControl->setInitIndex(tscPL.get<
int> (
"Initial Time Index"));
83 timeStepControl->setInitTime (tscPL.get<
double>(
"Initial Time"));
84 timeStepControl->setFinalTime(tscPL.get<
double>(
"Final Time"));
85 timeStepControl->setInitTimeStep(dt);
86 timeStepControl->initialize();
89 auto inArgsIC = model->getNominalValues();
90 auto icX = rcp_const_cast<Thyra::VectorBase<double> > (inArgsIC.get_x());
91 auto icXDot = rcp_const_cast<Thyra::VectorBase<double> > (inArgsIC.get_x_dot());
92 auto icXDotDot = rcp_const_cast<Thyra::VectorBase<double> > (inArgsIC.get_x_dot_dot());
93 auto icState = Tempus::createSolutionStateX<double>(icX, icXDot, icXDotDot);
94 icState->setTime (timeStepControl->getInitTime());
95 icState->setIndex (timeStepControl->getInitIndex());
96 icState->setTimeStep(0.0);
97 icState->setOrder (stepper->getOrder());
102 solutionHistory->setName(
"Forward States");
104 solutionHistory->setStorageLimit(2);
105 solutionHistory->addState(icState);
108 stepper->setInitialConditions(solutionHistory);
111 RCP<Tempus::IntegratorBasic<double> > integrator =
112 Tempus::createIntegratorBasic<double>();
113 integrator->setStepper(stepper);
114 integrator->setTimeStepControl(timeStepControl);
115 integrator->setSolutionHistory(solutionHistory);
117 integrator->initialize();
121 bool integratorStatus = integrator->advanceTime();
122 TEST_ASSERT(integratorStatus)
126 double time = integrator->getTime();
127 double timeFinal =pl->sublist(
"Default Integrator")
128 .sublist(
"Time Step Control").get<
double>(
"Final Time");
129 TEST_FLOATING_EQUALITY(time, timeFinal, 1.0e-14);
132 RCP<Thyra::VectorBase<double> > x = integrator->getX();
133 RCP<const Thyra::VectorBase<double> > x_exact =
134 model->getExactSolution(time).get_x();
137 RCP<Thyra::VectorBase<double> > xdiff = x->clone_v();
138 Thyra::V_StVpStV(xdiff.ptr(), 1.0, *x_exact, -1.0, *(x));
141 out <<
" Stepper = " << stepper->description()
142 <<
"\n with " << option << std::endl;
143 out <<
" =========================" << std::endl;
144 out <<
" Exact solution : " << get_ele(*(x_exact), 0) << std::endl;
145 out <<
" Computed solution: " << get_ele(*(x ), 0) << std::endl;
146 out <<
" Difference : " << get_ele(*(xdiff ), 0) << std::endl;
147 out <<
" =========================" << std::endl;
148 TEST_FLOATING_EQUALITY(get_ele(*(x), 0), 0.167158, 1.0e-4 );
157 RCP<Tempus::IntegratorBasic<double> > integrator;
158 std::vector<RCP<Thyra::VectorBase<double>>> solutions;
159 std::vector<RCP<Thyra::VectorBase<double>>> solutionsDot;
160 std::vector<double> StepSize;
161 std::vector<double> xErrorNorm;
162 std::vector<double> xDotErrorNorm;
163 const int nTimeStepSizes = 9;
167 RCP<ParameterList> pList =
168 getParametersFromXmlFile(
"Tempus_Leapfrog_SinCos.xml");
171 RCP<ParameterList> hom_pl = sublist(pList,
"HarmonicOscillatorModel",
true);
176 RCP<ParameterList> pl = sublist(pList,
"Tempus",
true);
179 double dt =pl->sublist(
"Default Integrator")
180 .sublist(
"Time Step Control").get<
double>(
"Initial Time Step");
183 for (
int n=0; n<nTimeStepSizes; n++) {
187 out <<
"\n \n time step #" << n
188 <<
" (out of " << nTimeStepSizes-1 <<
"), dt = " << dt <<
"\n";
189 pl->sublist(
"Default Integrator")
190 .sublist(
"Time Step Control").set(
"Initial Time Step", dt);
191 integrator = Tempus::createIntegratorBasic<double>(pl, model);
194 bool integratorStatus = integrator->advanceTime();
195 TEST_ASSERT(integratorStatus)
198 time = integrator->getTime();
199 double timeFinal =pl->sublist(
"Default Integrator")
200 .sublist(
"Time Step Control").get<
double>(
"Final Time");
201 TEST_FLOATING_EQUALITY(time, timeFinal, 1.0e-14);
204 if (n == nTimeStepSizes-1) {
205 RCP<const SolutionHistory<double> > solutionHistory =
206 integrator->getSolutionHistory();
207 writeSolution(
"Tempus_Leapfrog_SinCos.dat", solutionHistory);
210 for (
int i=0; i<solutionHistory->getNumStates(); i++) {
211 double time_i = (*solutionHistory)[i]->getTime();
214 model->getExactSolution(time_i).get_x()),
216 model->getExactSolution(time_i).get_x_dot()));
217 state->setTime((*solutionHistory)[i]->getTime());
218 solnHistExact->addState(state);
220 writeSolution(
"Tempus_Leapfrog_SinCos-Ref.dat", solnHistExact);
226 StepSize.push_back(dt);
227 auto solution = Thyra::createMember(model->get_x_space());
228 Thyra::copy(*(integrator->getX()),solution.ptr());
229 solutions.push_back(solution);
230 auto solutionDot = Thyra::createMember(model->get_x_space());
231 Thyra::copy(*(integrator->getXDot()),solutionDot.ptr());
232 solutionsDot.push_back(solutionDot);
233 if (n == nTimeStepSizes-1) {
234 StepSize.push_back(0.0);
235 auto solutionExact = Thyra::createMember(model->get_x_space());
236 Thyra::copy(*(model->getExactSolution(time).get_x()),solutionExact.ptr());
237 solutions.push_back(solutionExact);
238 auto solutionDotExact = Thyra::createMember(model->get_x_space());
239 Thyra::copy(*(model->getExactSolution(time).get_x_dot()),
240 solutionDotExact.ptr());
241 solutionsDot.push_back(solutionDotExact);
247 double xDotSlope = 0.0;
248 RCP<Tempus::Stepper<double> > stepper = integrator->getStepper();
249 double order = stepper->getOrder();
252 solutions, xErrorNorm, xSlope,
253 solutionsDot, xDotErrorNorm, xDotSlope);
255 TEST_FLOATING_EQUALITY( xSlope, order, 0.02 );
256 TEST_FLOATING_EQUALITY( xErrorNorm[0], 0.0157928, 1.0e-4 );
257 TEST_FLOATING_EQUALITY( xDotSlope, 1.09387, 0.01 );
258 TEST_FLOATING_EQUALITY( xDotErrorNorm[0], 0.563002, 1.0e-4 );
260 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...
TimeStepControl manages the time step size. There several mechanisms that effect the time step size a...
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.