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"
18#include "../TestUtils/Tempus_ConvergenceTestUtils.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
47using Tempus::IntegratorBasic;
48using Tempus::SolutionHistory;
49using Tempus::SolutionState;
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}
Tempus_ConvergenceTestUtils.hpp
Tempus::SolutionHistory
SolutionHistory is basically a container of SolutionStates. SolutionHistory maintains a collection of...
Definition Tempus_SolutionHistory_decl.hpp:121
Tempus::TimeStepControl
TimeStepControl manages the time step size. There several mechanisms that effect the time step size a...
Definition Tempus_TimeStepControl_decl.hpp:53
Tempus_Test::HarmonicOscillatorModel
Consider the ODE:
Definition HarmonicOscillatorModel_decl.hpp:55
Tempus_Test
Definition Tempus_BackwardEuler_ASA.cpp:34
Tempus_Test::writeOrderError
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)
Definition Tempus_ConvergenceTestUtils.hpp:189
Tempus_Test::writeSolution
void writeSolution(const std::string filename, Teuchos::RCP< const Tempus::SolutionHistory< Scalar > > solutionHistory)
Definition Tempus_ConvergenceTestUtils.hpp:311
Tempus_Test::TEUCHOS_UNIT_TEST
TEUCHOS_UNIT_TEST(BackwardEuler, SinCos_ASA)
Definition Tempus_BackwardEuler_ASA.cpp:47
Tempus::STORAGE_TYPE_STATIC
@ STORAGE_TYPE_STATIC
Keep a fix number of states.
Definition Tempus_SolutionHistory_decl.hpp:23
Tempus::createSolutionStateX
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.
Definition Tempus_SolutionState_impl.hpp:355
Tempus::createIntegratorBasic
Teuchos::RCP< IntegratorBasic< Scalar > > createIntegratorBasic()
Nonmember constructor.
Definition Tempus_IntegratorBasic_impl.hpp:730
Generated by doxygen 1.17.0