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"
27 namespace Tempus_Test {
29 using Teuchos::getParametersFromXmlFile;
33 using Teuchos::rcp_const_cast;
34 using Teuchos::sublist;
45 std::vector<std::string> options;
46 options.push_back(
"Default Parameters");
47 options.push_back(
"ICConsistency and Check");
49 for (
const auto& option : options) {
52 getParametersFromXmlFile(
"Tempus_Leapfrog_SinCos.xml");
61 stepper->setModel(model);
62 if (option ==
"ICConsistency and Check") {
63 stepper->setICConsistency(
"Consistent");
64 stepper->setICConsistencyCheck(
true);
66 stepper->initialize();
72 timeStepControl->setInitIndex(tscPL.
get<
int>(
"Initial Time Index"));
73 timeStepControl->setInitTime(tscPL.
get<
double>(
"Initial Time"));
74 timeStepControl->setFinalTime(tscPL.
get<
double>(
"Final Time"));
75 timeStepControl->setInitTimeStep(dt);
76 timeStepControl->initialize();
79 auto inArgsIC = model->getNominalValues();
85 auto icState = Tempus::createSolutionStateX<double>(icX, icXDot, icXDotDot);
86 icState->setTime(timeStepControl->getInitTime());
87 icState->setIndex(timeStepControl->getInitIndex());
88 icState->setTimeStep(0.0);
89 icState->setOrder(stepper->getOrder());
94 solutionHistory->setName(
"Forward States");
96 solutionHistory->setStorageLimit(2);
97 solutionHistory->addState(icState);
100 stepper->setInitialConditions(solutionHistory);
104 Tempus::createIntegratorBasic<double>();
105 integrator->setStepper(stepper);
106 integrator->setTimeStepControl(timeStepControl);
107 integrator->setSolutionHistory(solutionHistory);
109 integrator->initialize();
112 bool integratorStatus = integrator->advanceTime();
116 double time = integrator->getTime();
117 double timeFinal = pl->
sublist(
"Default Integrator")
119 .
get<
double>(
"Final Time");
125 model->getExactSolution(time).get_x();
129 Thyra::V_StVpStV(xdiff.
ptr(), 1.0, *x_exact, -1.0, *(x));
132 out <<
" Stepper = " << stepper->description() <<
"\n with "
133 << option << std::endl;
134 out <<
" =========================" << std::endl;
135 out <<
" Exact solution : " << get_ele(*(x_exact), 0) << std::endl;
136 out <<
" Computed solution: " << get_ele(*(x), 0) << std::endl;
137 out <<
" Difference : " << get_ele(*(xdiff), 0) << std::endl;
138 out <<
" =========================" << std::endl;
139 TEST_FLOATING_EQUALITY(get_ele(*(x), 0), 0.167158, 1.0e-4);
148 std::vector<RCP<Thyra::VectorBase<double>>> solutions;
149 std::vector<RCP<Thyra::VectorBase<double>>> solutionsDot;
150 std::vector<double> StepSize;
151 std::vector<double> xErrorNorm;
152 std::vector<double> xDotErrorNorm;
153 const int nTimeStepSizes = 9;
158 getParametersFromXmlFile(
"Tempus_Leapfrog_SinCos.xml");
169 double dt = pl->
sublist(
"Default Integrator")
171 .
get<
double>(
"Initial Time Step");
174 for (
int n = 0; n < nTimeStepSizes; n++) {
177 out <<
"\n \n time step #" << n <<
" (out of " << nTimeStepSizes - 1
178 <<
"), dt = " << dt <<
"\n";
179 pl->
sublist(
"Default Integrator")
181 .
set(
"Initial Time Step", dt);
182 integrator = Tempus::createIntegratorBasic<double>(pl, model);
185 bool integratorStatus = integrator->advanceTime();
189 time = integrator->getTime();
190 double timeFinal = pl->sublist(
"Default Integrator")
191 .sublist(
"Time Step Control")
192 .
get<
double>(
"Final Time");
196 if (n == nTimeStepSizes - 1) {
198 integrator->getSolutionHistory();
199 writeSolution(
"Tempus_Leapfrog_SinCos.dat", solutionHistory);
202 for (
int i = 0; i < solutionHistory->getNumStates(); i++) {
203 double time_i = (*solutionHistory)[i]->getTime();
206 model->getExactSolution(time_i).get_x()),
208 model->getExactSolution(time_i).get_x_dot()));
209 state->setTime((*solutionHistory)[i]->getTime());
210 solnHistExact->addState(state);
212 writeSolution(
"Tempus_Leapfrog_SinCos-Ref.dat", solnHistExact);
217 StepSize.push_back(dt);
218 auto solution = Thyra::createMember(model->get_x_space());
219 Thyra::copy(*(integrator->getX()), solution.ptr());
220 solutions.push_back(solution);
221 auto solutionDot = Thyra::createMember(model->get_x_space());
222 Thyra::copy(*(integrator->getXDot()), solutionDot.ptr());
223 solutionsDot.push_back(solutionDot);
224 if (n == nTimeStepSizes - 1) {
225 StepSize.push_back(0.0);
226 auto solutionExact = Thyra::createMember(model->get_x_space());
227 Thyra::copy(*(model->getExactSolution(time).get_x()),
228 solutionExact.ptr());
229 solutions.push_back(solutionExact);
230 auto solutionDotExact = Thyra::createMember(model->get_x_space());
231 Thyra::copy(*(model->getExactSolution(time).get_x_dot()),
232 solutionDotExact.ptr());
233 solutionsDot.push_back(solutionDotExact);
239 double xDotSlope = 0.0;
241 double order = stepper->getOrder();
243 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.
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)
Consider the ODE: where is a constant, is a constant damping parameter, is a constant forcing par...
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.
ParameterList & sublist(const std::string &name, bool mustAlreadyExist=false, const std::string &docString="")
Solution state for integrators and steppers.