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"
28 namespace Tempus_Test {
32 using Teuchos::rcp_const_cast;
34 using Teuchos::sublist;
35 using Teuchos::getParametersFromXmlFile;
48 std::vector<std::string> options;
49 options.push_back(
"Default Parameters");
50 options.push_back(
"ICConsistency and Check");
52 for(
const auto& option: options) {
56 getParametersFromXmlFile(
"Tempus_Leapfrog_SinCos.xml");
65 stepper->setModel(model);
66 if ( option ==
"ICConsistency and Check") {
67 stepper->setICConsistency(
"Consistent");
68 stepper->setICConsistencyCheck(
true);
70 stepper->initialize();
76 timeStepControl->setInitIndex(tscPL.
get<
int> (
"Initial Time Index"));
77 timeStepControl->setInitTime (tscPL.
get<
double>(
"Initial Time"));
78 timeStepControl->setFinalTime(tscPL.
get<
double>(
"Final Time"));
79 timeStepControl->setInitTimeStep(dt);
80 timeStepControl->initialize();
83 auto inArgsIC = model->getNominalValues();
87 auto icState = Tempus::createSolutionStateX<double>(icX, icXDot, icXDotDot);
88 icState->setTime (timeStepControl->getInitTime());
89 icState->setIndex (timeStepControl->getInitIndex());
90 icState->setTimeStep(0.0);
91 icState->setOrder (stepper->getOrder());
96 solutionHistory->setName(
"Forward States");
98 solutionHistory->setStorageLimit(2);
99 solutionHistory->addState(icState);
102 stepper->setInitialConditions(solutionHistory);
106 Tempus::createIntegratorBasic<double>();
107 integrator->setStepper(stepper);
108 integrator->setTimeStepControl(timeStepControl);
109 integrator->setSolutionHistory(solutionHistory);
111 integrator->initialize();
115 bool integratorStatus = integrator->advanceTime();
120 double time = integrator->getTime();
121 double timeFinal =pl->
sublist(
"Default Integrator")
122 .
sublist(
"Time Step Control").
get<
double>(
"Final Time");
128 model->getExactSolution(time).get_x();
132 Thyra::V_StVpStV(xdiff.
ptr(), 1.0, *x_exact, -1.0, *(x));
135 out <<
" Stepper = " << stepper->description()
136 <<
"\n with " << option << std::endl;
137 out <<
" =========================" << std::endl;
138 out <<
" Exact solution : " << get_ele(*(x_exact), 0) << std::endl;
139 out <<
" Computed solution: " << get_ele(*(x ), 0) << std::endl;
140 out <<
" Difference : " << get_ele(*(xdiff ), 0) << std::endl;
141 out <<
" =========================" << std::endl;
142 TEST_FLOATING_EQUALITY(get_ele(*(x), 0), 0.167158, 1.0e-4 );
152 std::vector<RCP<Thyra::VectorBase<double>>> solutions;
153 std::vector<RCP<Thyra::VectorBase<double>>> solutionsDot;
154 std::vector<double> StepSize;
155 std::vector<double> xErrorNorm;
156 std::vector<double> xDotErrorNorm;
157 const int nTimeStepSizes = 9;
162 getParametersFromXmlFile(
"Tempus_Leapfrog_SinCos.xml");
173 double dt =pl->
sublist(
"Default Integrator")
174 .
sublist(
"Time Step Control").
get<
double>(
"Initial Time Step");
177 for (
int n=0; n<nTimeStepSizes; n++) {
181 out <<
"\n \n time step #" << n
182 <<
" (out of " << nTimeStepSizes-1 <<
"), dt = " << dt <<
"\n";
183 pl->
sublist(
"Default Integrator")
184 .
sublist(
"Time Step Control").
set(
"Initial Time Step", dt);
185 integrator = Tempus::createIntegratorBasic<double>(pl, model);
188 bool integratorStatus = integrator->advanceTime();
192 time = integrator->getTime();
193 double timeFinal =pl->sublist(
"Default Integrator")
194 .sublist(
"Time Step Control").
get<
double>(
"Final Time");
198 if (n == nTimeStepSizes-1) {
200 integrator->getSolutionHistory();
201 writeSolution(
"Tempus_Leapfrog_SinCos.dat", solutionHistory);
204 for (
int i=0; i<solutionHistory->getNumStates(); i++) {
205 double time_i = (*solutionHistory)[i]->getTime();
208 model->getExactSolution(time_i).get_x()),
210 model->getExactSolution(time_i).get_x_dot()));
211 state->setTime((*solutionHistory)[i]->getTime());
212 solnHistExact->addState(state);
214 writeSolution(
"Tempus_Leapfrog_SinCos-Ref.dat", solnHistExact);
220 StepSize.push_back(dt);
221 auto solution = Thyra::createMember(model->get_x_space());
222 Thyra::copy(*(integrator->getX()),solution.ptr());
223 solutions.push_back(solution);
224 auto solutionDot = Thyra::createMember(model->get_x_space());
225 Thyra::copy(*(integrator->getXDot()),solutionDot.ptr());
226 solutionsDot.push_back(solutionDot);
227 if (n == nTimeStepSizes-1) {
228 StepSize.push_back(0.0);
229 auto solutionExact = Thyra::createMember(model->get_x_space());
230 Thyra::copy(*(model->getExactSolution(time).get_x()),solutionExact.ptr());
231 solutions.push_back(solutionExact);
232 auto solutionDotExact = Thyra::createMember(model->get_x_space());
233 Thyra::copy(*(model->getExactSolution(time).get_x_dot()),
234 solutionDotExact.ptr());
235 solutionsDot.push_back(solutionDotExact);
241 double xDotSlope = 0.0;
243 double order = stepper->getOrder();
246 solutions, xErrorNorm, xSlope,
247 solutionsDot, xDotErrorNorm, xDotSlope);
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)
#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...
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_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.