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