Tempus Version of the Day
Time Integration
Loading...
Searching...
No Matches
Tempus_HHTAlphaTest.cpp
Go to the documentation of this file.
1// @HEADER
2// ****************************************************************************
3// Tempus: Copyright (2017) Sandia Corporation
4//
5// Distributed under BSD 3-clause license (See accompanying file Copyright.txt)
6// ****************************************************************************
7// @HEADER
8
9#include "Teuchos_UnitTestHarness.hpp"
10#include "Teuchos_XMLParameterListHelpers.hpp"
11#include "Teuchos_TimeMonitor.hpp"
12
13#include "Tempus_config.hpp"
14#include "Tempus_IntegratorBasic.hpp"
15#include "Tempus_StepperHHTAlpha.hpp"
16
17#include "../TestModels/HarmonicOscillatorModel.hpp"
19
20#include "Stratimikos_DefaultLinearSolverBuilder.hpp"
21#include "Thyra_LinearOpWithSolveFactoryHelpers.hpp"
22#include "Thyra_DetachedVectorView.hpp"
23#include "Thyra_DetachedMultiVectorView.hpp"
24#include "NOX_Thyra.H"
25
26
27#ifdef Tempus_ENABLE_MPI
28#include "Epetra_MpiComm.h"
29#else
30#include "Epetra_SerialComm.h"
31#endif
32
33#include <fstream>
34#include <limits>
35#include <sstream>
36#include <vector>
37
38namespace Tempus_Test {
39
40using Teuchos::RCP;
41using Teuchos::rcp;
42using Teuchos::rcp_const_cast;
43using Teuchos::ParameterList;
44using Teuchos::sublist;
45using Teuchos::getParametersFromXmlFile;
46
50
51
52// ************************************************************
53// ************************************************************
54TEUCHOS_UNIT_TEST(HHTAlpha, BallParabolic)
55{
56 //Tolerance to check if test passed
57 double tolerance = 1.0e-14;
58 // Read params from .xml file
59 RCP<ParameterList> pList =
60 getParametersFromXmlFile("Tempus_HHTAlpha_BallParabolic.xml");
61
62 // Setup the HarmonicOscillatorModel
63 RCP<ParameterList> hom_pl = sublist(pList, "HarmonicOscillatorModel", true);
64 auto model = rcp(new HarmonicOscillatorModel<double>(hom_pl));
65
66 // Setup the Integrator and reset initial time step
67 RCP<ParameterList> pl = sublist(pList, "Tempus", true);
68
69 RCP<Tempus::IntegratorBasic<double> > integrator =
70 Tempus::createIntegratorBasic<double>(pl, model);
71
72 // Integrate to timeMax
73 bool integratorStatus = integrator->advanceTime();
74 TEST_ASSERT(integratorStatus)
75
76 // Test if at 'Final Time'
77 double time = integrator->getTime();
78 double timeFinal =pl->sublist("Default Integrator")
79 .sublist("Time Step Control").get<double>("Final Time");
80 TEST_FLOATING_EQUALITY(time, timeFinal, 1.0e-14);
81
82 // Time-integrated solution and the exact solution
83 RCP<Thyra::VectorBase<double> > x = integrator->getX();
84 RCP<const Thyra::VectorBase<double> > x_exact =
85 model->getExactSolution(time).get_x();
86
87 // Plot sample solution and exact solution
88 std::ofstream ftmp("Tempus_HHTAlpha_BallParabolic.dat");
89 ftmp.precision(16);
90 RCP<const SolutionHistory<double> > solutionHistory =
91 integrator->getSolutionHistory();
92 bool passed = true;
93 double err = 0.0;
94 RCP<const Thyra::VectorBase<double> > x_exact_plot;
95 for (int i=0; i<solutionHistory->getNumStates(); i++) {
96 RCP<const SolutionState<double> > solutionState = (*solutionHistory)[i];
97 double time_i = solutionState->getTime();
98 RCP<const Thyra::VectorBase<double> > x_plot = solutionState->getX();
99 x_exact_plot = model->getExactSolution(time_i).get_x();
100 ftmp << time_i << " "
101 << get_ele(*(x_plot), 0) << " "
102 << get_ele(*(x_exact_plot), 0) << std::endl;
103 if (abs(get_ele(*(x_plot),0) - get_ele(*(x_exact_plot), 0)) > err)
104 err = abs(get_ele(*(x_plot),0) - get_ele(*(x_exact_plot), 0));
105 }
106 ftmp.close();
107 out << "Max error = " << err << "\n \n";
108 if (err > tolerance)
109 passed = false;
110
111 TEUCHOS_TEST_FOR_EXCEPTION(!passed, std::logic_error,
112 "\n Test failed! Max error = " << err << " > tolerance = " << tolerance << "\n!");
113
114}
115
116
117// ************************************************************
118// ************************************************************
119TEUCHOS_UNIT_TEST(HHTAlpha, ConstructingFromDefaults)
120{
121 double dt = 0.05;
122
123 // Read params from .xml file
124 RCP<ParameterList> pList =
125 getParametersFromXmlFile("Tempus_HHTAlpha_SinCos_SecondOrder.xml");
126 RCP<ParameterList> pl = sublist(pList, "Tempus", true);
127
128 // Setup the HarmonicOscillatorModel
129 RCP<ParameterList> hom_pl = sublist(pList, "HarmonicOscillatorModel", true);
130 auto model = rcp(new HarmonicOscillatorModel<double>(hom_pl));
131
132 // Setup Stepper for field solve ----------------------------
133 auto stepper = rcp(new Tempus::StepperHHTAlpha<double>());
134 stepper->setModel(model);
135 stepper->initialize();
136
137 // Setup TimeStepControl ------------------------------------
138 auto timeStepControl = rcp(new Tempus::TimeStepControl<double>());
139 ParameterList tscPL = pl->sublist("Default Integrator")
140 .sublist("Time Step Control");
141 timeStepControl->setInitIndex(tscPL.get<int> ("Initial Time Index"));
142 timeStepControl->setInitTime (tscPL.get<double>("Initial Time"));
143 timeStepControl->setFinalTime(tscPL.get<double>("Final Time"));
144 timeStepControl->setInitTimeStep(dt);
145 timeStepControl->initialize();
146
147 // Setup initial condition SolutionState --------------------
148 auto inArgsIC = model->getNominalValues();
149 auto icX = rcp_const_cast<Thyra::VectorBase<double> > (inArgsIC.get_x());
150 auto icXDot = rcp_const_cast<Thyra::VectorBase<double> > (inArgsIC.get_x_dot());
151 auto icXDotDot = rcp_const_cast<Thyra::VectorBase<double> > (inArgsIC.get_x_dot_dot());
152 auto icState = Tempus::createSolutionStateX(icX, icXDot, icXDotDot);
153 icState->setTime (timeStepControl->getInitTime());
154 icState->setIndex (timeStepControl->getInitIndex());
155 icState->setTimeStep(0.0);
156 icState->setOrder (stepper->getOrder());
157 icState->setSolutionStatus(Tempus::Status::PASSED); // ICs are passing.
158
159 // Setup SolutionHistory ------------------------------------
160 auto solutionHistory = rcp(new Tempus::SolutionHistory<double>());
161 solutionHistory->setName("Forward States");
162 solutionHistory->setStorageType(Tempus::STORAGE_TYPE_STATIC);
163 solutionHistory->setStorageLimit(2);
164 solutionHistory->addState(icState);
165
166 // Setup Integrator -----------------------------------------
167 RCP<Tempus::IntegratorBasic<double> > integrator =
168 Tempus::createIntegratorBasic<double>();
169 integrator->setStepper(stepper);
170 integrator->setTimeStepControl(timeStepControl);
171 integrator->setSolutionHistory(solutionHistory);
172 //integrator->setObserver(...);
173 integrator->initialize();
174
175
176 // Integrate to timeMax
177 bool integratorStatus = integrator->advanceTime();
178 TEST_ASSERT(integratorStatus)
179
180
181 // Test if at 'Final Time'
182 double time = integrator->getTime();
183 double timeFinal =pl->sublist("Default Integrator")
184 .sublist("Time Step Control").get<double>("Final Time");
185 TEST_FLOATING_EQUALITY(time, timeFinal, 1.0e-14);
186
187 // Time-integrated solution and the exact solution
188 RCP<Thyra::VectorBase<double> > x = integrator->getX();
189 RCP<const Thyra::VectorBase<double> > x_exact =
190 model->getExactSolution(time).get_x();
191
192 // Calculate the error
193 RCP<Thyra::VectorBase<double> > xdiff = x->clone_v();
194 Thyra::V_StVpStV(xdiff.ptr(), 1.0, *x_exact, -1.0, *(x));
195
196 // Check the order and intercept
197 out << " Stepper = " << stepper->description() << std::endl;
198 out << " =========================" << std::endl;
199 out << " Exact solution : " << get_ele(*(x_exact), 0) << std::endl;
200 out << " Computed solution: " << get_ele(*(x ), 0) << std::endl;
201 out << " Difference : " << get_ele(*(xdiff ), 0) << std::endl;
202 out << " =========================" << std::endl;
203 TEST_FLOATING_EQUALITY(get_ele(*(x), 0), 0.144918, 1.0e-4 );
204}
205
206
207// ************************************************************
208TEUCHOS_UNIT_TEST(HHTAlpha, SinCos_SecondOrder)
209{
210 RCP<Tempus::IntegratorBasic<double> > integrator;
211 std::vector<RCP<Thyra::VectorBase<double>>> solutions;
212 std::vector<RCP<Thyra::VectorBase<double>>> solutionsDot;
213 std::vector<double> StepSize;
214 std::vector<double> xErrorNorm;
215 std::vector<double> xDotErrorNorm;
216 const int nTimeStepSizes = 7;
217 double time = 0.0;
218
219 // Read params from .xml file
220 RCP<ParameterList> pList =
221 getParametersFromXmlFile("Tempus_HHTAlpha_SinCos_SecondOrder.xml");
222
223 // Setup the HarmonicOscillatorModel
224 RCP<ParameterList> hom_pl = sublist(pList, "HarmonicOscillatorModel", true);
225 auto model = rcp(new HarmonicOscillatorModel<double>(hom_pl));
226
227 //Get k and m coefficients from model, needed for computing energy
228 double k = hom_pl->get<double>("x coeff k");
229 double m = hom_pl->get<double>("Mass coeff m");
230
231 // Setup the Integrator and reset initial time step
232 RCP<ParameterList> pl = sublist(pList, "Tempus", true);
233
234 //Set initial time step = 2*dt specified in input file (for convergence study)
235 //
236 double dt =pl->sublist("Default Integrator")
237 .sublist("Time Step Control").get<double>("Initial Time Step");
238 dt *= 2.0;
239
240 for (int n=0; n<nTimeStepSizes; n++) {
241
242 //Perform time-step refinement
243 dt /= 2;
244 out << "\n \n time step #" << n << " (out of "
245 << nTimeStepSizes-1 << "), dt = " << dt << "\n";
246 pl->sublist("Default Integrator")
247 .sublist("Time Step Control").set("Initial Time Step", dt);
248 integrator = Tempus::createIntegratorBasic<double>(pl, model);
249
250 // Integrate to timeMax
251 bool integratorStatus = integrator->advanceTime();
252 TEST_ASSERT(integratorStatus)
253
254 // Test if at 'Final Time'
255 time = integrator->getTime();
256 double timeFinal =pl->sublist("Default Integrator")
257 .sublist("Time Step Control").get<double>("Final Time");
258 TEST_FLOATING_EQUALITY(time, timeFinal, 1.0e-14);
259
260 // Time-integrated solution and the exact solution
261 RCP<Thyra::VectorBase<double> > x = integrator->getX();
262 RCP<const Thyra::VectorBase<double> > x_exact =
263 model->getExactSolution(time).get_x();
264
265 // Plot sample solution and exact solution at most-refined resolution
266 if (n == nTimeStepSizes-1) {
267 RCP<const SolutionHistory<double> > solutionHistory =
268 integrator->getSolutionHistory();
269 writeSolution("Tempus_HHTAlpha_SinCos_SecondOrder.dat", solutionHistory);
270
271 auto solnHistExact = rcp(new Tempus::SolutionHistory<double>());
272 for (int i=0; i<solutionHistory->getNumStates(); i++) {
273 double time_i = (*solutionHistory)[i]->getTime();
274 auto state = Tempus::createSolutionStateX(
275 rcp_const_cast<Thyra::VectorBase<double> > (
276 model->getExactSolution(time_i).get_x()),
277 rcp_const_cast<Thyra::VectorBase<double> > (
278 model->getExactSolution(time_i).get_x_dot()));
279 state->setTime((*solutionHistory)[i]->getTime());
280 solnHistExact->addState(state);
281 }
282 writeSolution("Tempus_HHTAlpha_SinCos_SecondOrder-Ref.dat", solnHistExact);
283
284 // Problem specific output
285 {
286 std::ofstream ftmp("Tempus_HHTAlpha_SinCos_SecondOrder-Energy.dat");
287 ftmp.precision(16);
288 RCP<const Thyra::VectorBase<double> > x_exact_plot;
289 for (int i=0; i<solutionHistory->getNumStates(); i++) {
290 RCP<const SolutionState<double> > solutionState = (*solutionHistory)[i];
291 double time_i = solutionState->getTime();
292 RCP<const Thyra::VectorBase<double> > x_plot = solutionState->getX();
293 RCP<const Thyra::VectorBase<double> > x_dot_plot = solutionState->getXDot();
294 x_exact_plot = model->getExactSolution(time_i).get_x();
295 //kinetic energy = 0.5*m*xdot*xdot
296 double ke = Thyra::dot(*x_dot_plot, *x_dot_plot);
297 ke *= 0.5*m;
298 //potential energy = 0.5*k*x*x
299 double pe = Thyra::dot(*x_plot, *x_plot);
300 pe *= 0.5*k;
301 double te = ke + pe;
302 //Output to file the following:
303 //[time, x computed, x exact, xdot computed, ke, pe, te]
304 ftmp << time_i << " "
305 << get_ele(*(x_plot), 0) << " "
306 << get_ele(*(x_exact_plot), 0) << " "
307 << get_ele(*(x_dot_plot), 0) << " "
308 << ke << " " << pe << " " << te << std::endl;
309 }
310 ftmp.close();
311 }
312 }
313
314 // Store off the final solution and step size
315 StepSize.push_back(dt);
316 auto solution = Thyra::createMember(model->get_x_space());
317 Thyra::copy(*(integrator->getX()),solution.ptr());
318 solutions.push_back(solution);
319 auto solutionDot = Thyra::createMember(model->get_x_space());
320 Thyra::copy(*(integrator->getXDot()),solutionDot.ptr());
321 solutionsDot.push_back(solutionDot);
322 if (n == nTimeStepSizes-1) { // Add exact solution last in vector.
323 StepSize.push_back(0.0);
324 auto solutionExact = Thyra::createMember(model->get_x_space());
325 Thyra::copy(*(model->getExactSolution(time).get_x()),solutionExact.ptr());
326 solutions.push_back(solutionExact);
327 auto solutionDotExact = Thyra::createMember(model->get_x_space());
328 Thyra::copy(*(model->getExactSolution(time).get_x_dot()),
329 solutionDotExact.ptr());
330 solutionsDot.push_back(solutionDotExact);
331 }
332 }
333
334 // Check the order and intercept
335 double xSlope = 0.0;
336 double xDotSlope = 0.0;
337 RCP<Tempus::Stepper<double> > stepper = integrator->getStepper();
338 double order = stepper->getOrder();
339 writeOrderError("Tempus_HHTAlpha_SinCos_SecondOrder-Error.dat",
340 stepper, StepSize,
341 solutions, xErrorNorm, xSlope,
342 solutionsDot, xDotErrorNorm, xDotSlope);
343
344 TEST_FLOATING_EQUALITY( xSlope, order, 0.02 );
345 TEST_FLOATING_EQUALITY( xErrorNorm[0], 0.00644755, 1.0e-4 );
346 TEST_FLOATING_EQUALITY( xDotSlope, order, 0.01 );
347 TEST_FLOATING_EQUALITY( xDotErrorNorm[0], 0.104392, 1.0e-4 );
348}
349
350
351// ************************************************************
352// ************************************************************
353TEUCHOS_UNIT_TEST(HHTAlpha, SinCos_FirstOrder)
354{
355 RCP<Tempus::IntegratorBasic<double> > integrator;
356 std::vector<RCP<Thyra::VectorBase<double>>> solutions;
357 std::vector<RCP<Thyra::VectorBase<double>>> solutionsDot;
358 std::vector<double> StepSize;
359 std::vector<double> xErrorNorm;
360 std::vector<double> xDotErrorNorm;
361 const int nTimeStepSizes = 7;
362 double time = 0.0;
363
364 // Read params from .xml file
365 RCP<ParameterList> pList =
366 getParametersFromXmlFile("Tempus_HHTAlpha_SinCos_FirstOrder.xml");
367
368 // Setup the HarmonicOscillatorModel
369 RCP<ParameterList> hom_pl = sublist(pList, "HarmonicOscillatorModel", true);
370 auto model = rcp(new HarmonicOscillatorModel<double>(hom_pl));
371
372 //Get k and m coefficients from model, needed for computing energy
373 double k = hom_pl->get<double>("x coeff k");
374 double m = hom_pl->get<double>("Mass coeff m");
375
376 // Setup the Integrator and reset initial time step
377 RCP<ParameterList> pl = sublist(pList, "Tempus", true);
378
379 //Set initial time step = 2*dt specified in input file (for convergence study)
380 //
381 double dt =pl->sublist("Default Integrator")
382 .sublist("Time Step Control").get<double>("Initial Time Step");
383 dt *= 2.0;
384
385 for (int n=0; n<nTimeStepSizes; n++) {
386
387 //Perform time-step refinement
388 dt /= 2;
389 out << "\n \n time step #" << n << " (out of "
390 << nTimeStepSizes-1 << "), dt = " << dt << "\n";
391 pl->sublist("Default Integrator")
392 .sublist("Time Step Control").set("Initial Time Step", dt);
393 integrator = Tempus::createIntegratorBasic<double>(pl, model);
394
395 // Integrate to timeMax
396 bool integratorStatus = integrator->advanceTime();
397 TEST_ASSERT(integratorStatus)
398
399 // Test if at 'Final Time'
400 time = integrator->getTime();
401 double timeFinal =pl->sublist("Default Integrator")
402 .sublist("Time Step Control").get<double>("Final Time");
403 TEST_FLOATING_EQUALITY(time, timeFinal, 1.0e-14);
404
405 // Time-integrated solution and the exact solution
406 RCP<Thyra::VectorBase<double> > x = integrator->getX();
407 RCP<const Thyra::VectorBase<double> > x_exact =
408 model->getExactSolution(time).get_x();
409
410 // Plot sample solution and exact solution at most-refined resolution
411 if (n == nTimeStepSizes-1) {
412 RCP<const SolutionHistory<double> > solutionHistory =
413 integrator->getSolutionHistory();
414 writeSolution("Tempus_HHTAlpha_SinCos_FirstOrder.dat", solutionHistory);
415
416 auto solnHistExact = rcp(new Tempus::SolutionHistory<double>());
417 for (int i=0; i<solutionHistory->getNumStates(); i++) {
418 double time_i = (*solutionHistory)[i]->getTime();
419 auto state = Tempus::createSolutionStateX(
420 rcp_const_cast<Thyra::VectorBase<double> > (
421 model->getExactSolution(time_i).get_x()),
422 rcp_const_cast<Thyra::VectorBase<double> > (
423 model->getExactSolution(time_i).get_x_dot()));
424 state->setTime((*solutionHistory)[i]->getTime());
425 solnHistExact->addState(state);
426 }
427 writeSolution("Tempus_HHTAlpha_SinCos_FirstOrder-Ref.dat", solnHistExact);
428
429 // Problem specific output
430 {
431 std::ofstream ftmp("Tempus_HHTAlpha_SinCos_FirstOrder-Energy.dat");
432 ftmp.precision(16);
433 RCP<const Thyra::VectorBase<double> > x_exact_plot;
434 for (int i=0; i<solutionHistory->getNumStates(); i++) {
435 RCP<const SolutionState<double> > solutionState = (*solutionHistory)[i];
436 double time_i = solutionState->getTime();
437 RCP<const Thyra::VectorBase<double> > x_plot = solutionState->getX();
438 RCP<const Thyra::VectorBase<double> > x_dot_plot = solutionState->getXDot();
439 x_exact_plot = model->getExactSolution(time_i).get_x();
440 //kinetic energy = 0.5*m*xdot*xdot
441 double ke = Thyra::dot(*x_dot_plot, *x_dot_plot);
442 ke *= 0.5*m;
443 //potential energy = 0.5*k*x*x
444 double pe = Thyra::dot(*x_plot, *x_plot);
445 pe *= 0.5*k;
446 double te = ke + pe;
447 //Output to file the following:
448 //[time, x computed, x exact, xdot computed, ke, pe, te]
449 ftmp << time_i << " "
450 << get_ele(*(x_plot), 0) << " "
451 << get_ele(*(x_exact_plot), 0) << " "
452 << get_ele(*(x_dot_plot), 0) << " "
453 << ke << " " << pe << " " << te << std::endl;
454 }
455 ftmp.close();
456 }
457 }
458
459 // Store off the final solution and step size
460 StepSize.push_back(dt);
461 auto solution = Thyra::createMember(model->get_x_space());
462 Thyra::copy(*(integrator->getX()),solution.ptr());
463 solutions.push_back(solution);
464 auto solutionDot = Thyra::createMember(model->get_x_space());
465 Thyra::copy(*(integrator->getXDot()),solutionDot.ptr());
466 solutionsDot.push_back(solutionDot);
467 if (n == nTimeStepSizes-1) { // Add exact solution last in vector.
468 StepSize.push_back(0.0);
469 auto solutionExact = Thyra::createMember(model->get_x_space());
470 Thyra::copy(*(model->getExactSolution(time).get_x()),solutionExact.ptr());
471 solutions.push_back(solutionExact);
472 auto solutionDotExact = Thyra::createMember(model->get_x_space());
473 Thyra::copy(*(model->getExactSolution(time).get_x_dot()),
474 solutionDotExact.ptr());
475 solutionsDot.push_back(solutionDotExact);
476 }
477 }
478
479 // Check the order and intercept
480 double xSlope = 0.0;
481 double xDotSlope = 0.0;
482 RCP<Tempus::Stepper<double> > stepper = integrator->getStepper();
483 //double order = stepper->getOrder();
484 writeOrderError("Tempus_HHTAlpha_SinCos_FirstOrder-Error.dat",
485 stepper, StepSize,
486 solutions, xErrorNorm, xSlope,
487 solutionsDot, xDotErrorNorm, xDotSlope);
488
489 TEST_FLOATING_EQUALITY( xSlope, 0.977568, 0.02 );
490 TEST_FLOATING_EQUALITY( xErrorNorm[0], 0.048932, 1.0e-4 );
491 TEST_FLOATING_EQUALITY( xDotSlope, 1.2263, 0.01 );
492 TEST_FLOATING_EQUALITY( xDotErrorNorm[0], 0.393504, 1.0e-4 );
493
494}
495
496
497// ************************************************************
498// ************************************************************
499TEUCHOS_UNIT_TEST(HHTAlpha, SinCos_CD)
500{
501 RCP<Tempus::IntegratorBasic<double> > integrator;
502 std::vector<RCP<Thyra::VectorBase<double>>> solutions;
503 std::vector<RCP<Thyra::VectorBase<double>>> solutionsDot;
504 std::vector<double> StepSize;
505 std::vector<double> xErrorNorm;
506 std::vector<double> xDotErrorNorm;
507 const int nTimeStepSizes = 7;
508 double time = 0.0;
509
510 // Read params from .xml file
511 RCP<ParameterList> pList =
512 getParametersFromXmlFile("Tempus_HHTAlpha_SinCos_ExplicitCD.xml");
513
514 // Setup the HarmonicOscillatorModel
515 RCP<ParameterList> hom_pl = sublist(pList, "HarmonicOscillatorModel", true);
516 auto model = rcp(new HarmonicOscillatorModel<double>(hom_pl));
517
518 //Get k and m coefficients from model, needed for computing energy
519 double k = hom_pl->get<double>("x coeff k");
520 double m = hom_pl->get<double>("Mass coeff m");
521
522 // Setup the Integrator and reset initial time step
523 RCP<ParameterList> pl = sublist(pList, "Tempus", true);
524
525 //Set initial time step = 2*dt specified in input file (for convergence study)
526 //
527 double dt =pl->sublist("Default Integrator")
528 .sublist("Time Step Control").get<double>("Initial Time Step");
529 dt *= 2.0;
530
531 for (int n=0; n<nTimeStepSizes; n++) {
532
533 //Perform time-step refinement
534 dt /= 2;
535 out << "\n \n time step #" << n << " (out of "
536 << nTimeStepSizes-1 << "), dt = " << dt << "\n";
537 pl->sublist("Default Integrator")
538 .sublist("Time Step Control").set("Initial Time Step", dt);
539 integrator = Tempus::createIntegratorBasic<double>(pl, model);
540
541 // Integrate to timeMax
542 bool integratorStatus = integrator->advanceTime();
543 TEST_ASSERT(integratorStatus)
544
545 // Test if at 'Final Time'
546 time = integrator->getTime();
547 double timeFinal =pl->sublist("Default Integrator")
548 .sublist("Time Step Control").get<double>("Final Time");
549 TEST_FLOATING_EQUALITY(time, timeFinal, 1.0e-14);
550
551 // Time-integrated solution and the exact solution
552 RCP<Thyra::VectorBase<double> > x = integrator->getX();
553 RCP<const Thyra::VectorBase<double> > x_exact =
554 model->getExactSolution(time).get_x();
555
556 // Plot sample solution and exact solution at most-refined resolution
557 if (n == nTimeStepSizes-1) {
558 RCP<const SolutionHistory<double> > solutionHistory =
559 integrator->getSolutionHistory();
560 writeSolution("Tempus_HHTAlpha_SinCos_ExplicitCD.dat", solutionHistory);
561
562 auto solnHistExact = rcp(new Tempus::SolutionHistory<double>());
563 for (int i=0; i<solutionHistory->getNumStates(); i++) {
564 double time_i = (*solutionHistory)[i]->getTime();
565 auto state = Tempus::createSolutionStateX(
566 rcp_const_cast<Thyra::VectorBase<double> > (
567 model->getExactSolution(time_i).get_x()),
568 rcp_const_cast<Thyra::VectorBase<double> > (
569 model->getExactSolution(time_i).get_x_dot()));
570 state->setTime((*solutionHistory)[i]->getTime());
571 solnHistExact->addState(state);
572 }
573 writeSolution("Tempus_HHTAlpha_SinCos_ExplicitCD-Ref.dat", solnHistExact);
574
575 // Problem specific output
576 {
577 std::ofstream ftmp("Tempus_HHTAlpha_SinCos_ExplicitCD-Energy.dat");
578 ftmp.precision(16);
579 RCP<const Thyra::VectorBase<double> > x_exact_plot;
580 for (int i=0; i<solutionHistory->getNumStates(); i++) {
581 RCP<const SolutionState<double> > solutionState = (*solutionHistory)[i];
582 double time_i = solutionState->getTime();
583 RCP<const Thyra::VectorBase<double> > x_plot = solutionState->getX();
584 RCP<const Thyra::VectorBase<double> > x_dot_plot = solutionState->getXDot();
585 x_exact_plot = model->getExactSolution(time_i).get_x();
586 //kinetic energy = 0.5*m*xdot*xdot
587 double ke = Thyra::dot(*x_dot_plot, *x_dot_plot);
588 ke *= 0.5*m;
589 //potential energy = 0.5*k*x*x
590 double pe = Thyra::dot(*x_plot, *x_plot);
591 pe *= 0.5*k;
592 double te = ke + pe;
593 //Output to file the following:
594 //[time, x computed, x exact, xdot computed, ke, pe, te]
595 ftmp << time_i << " "
596 << get_ele(*(x_plot), 0) << " "
597 << get_ele(*(x_exact_plot), 0) << " "
598 << get_ele(*(x_dot_plot), 0) << " "
599 << ke << " " << pe << " " << te << std::endl;
600 }
601 ftmp.close();
602 }
603 }
604
605 // Store off the final solution and step size
606 StepSize.push_back(dt);
607 auto solution = Thyra::createMember(model->get_x_space());
608 Thyra::copy(*(integrator->getX()),solution.ptr());
609 solutions.push_back(solution);
610 auto solutionDot = Thyra::createMember(model->get_x_space());
611 Thyra::copy(*(integrator->getXDot()),solutionDot.ptr());
612 solutionsDot.push_back(solutionDot);
613 if (n == nTimeStepSizes-1) { // Add exact solution last in vector.
614 StepSize.push_back(0.0);
615 auto solutionExact = Thyra::createMember(model->get_x_space());
616 Thyra::copy(*(model->getExactSolution(time).get_x()),solutionExact.ptr());
617 solutions.push_back(solutionExact);
618 auto solutionDotExact = Thyra::createMember(model->get_x_space());
619 Thyra::copy(*(model->getExactSolution(time).get_x_dot()),
620 solutionDotExact.ptr());
621 solutionsDot.push_back(solutionDotExact);
622 }
623 }
624
625 // Check the order and intercept
626 double xSlope = 0.0;
627 double xDotSlope = 0.0;
628 RCP<Tempus::Stepper<double> > stepper = integrator->getStepper();
629 double order = stepper->getOrder();
630 writeOrderError("Tempus_HHTAlpha_SinCos_ExplicitCD-Error.dat",
631 stepper, StepSize,
632 solutions, xErrorNorm, xSlope,
633 solutionsDot, xDotErrorNorm, xDotSlope);
634
635 TEST_FLOATING_EQUALITY( xSlope, order, 0.02 );
636 TEST_FLOATING_EQUALITY( xErrorNorm[0], 0.00451069, 1.0e-4 );
637 TEST_FLOATING_EQUALITY( xDotSlope, order, 0.01 );
638 TEST_FLOATING_EQUALITY( xDotErrorNorm[0], 0.0551522, 1.0e-4 );
639}
640
641
642}
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.