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);
67 auto model = rcp(
new HarmonicOscillatorModel<double>(hom_pl));
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);
172 auto model = rcp(
new HarmonicOscillatorModel<double>(hom_pl));
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...
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)