Tempus Version of the Day
Time Integration
Loading...
Searching...
No Matches
Tempus_ExplicitRKTest.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#include "Teuchos_DefaultComm.hpp"
14
15#include "Thyra_VectorStdOps.hpp"
16
17#include "Tempus_IntegratorBasic.hpp"
18#include "Tempus_StepperFactory.hpp"
19#include "Tempus_StepperExplicitRK.hpp"
20
21#include "../TestModels/SinCosModel.hpp"
22#include "../TestModels/VanDerPolModel.hpp"
23#include "../TestUtils/Tempus_ConvergenceTestUtils.hpp"
24
25#include <fstream>
26#include <vector>
27
28namespace Tempus_Test {
29
30using Teuchos::getParametersFromXmlFile;
31using Teuchos::ParameterList;
32using Teuchos::RCP;
33using Teuchos::rcp;
34using Teuchos::rcp_const_cast;
35using Teuchos::sublist;
36
37using Tempus::IntegratorBasic;
38using Tempus::SolutionHistory;
39using Tempus::SolutionState;
40
41// ************************************************************
42// ************************************************************
43TEUCHOS_UNIT_TEST(ExplicitRK, ParameterList)
44{
45 std::vector<std::string> RKMethods;
46 RKMethods.push_back("General ERK");
47 RKMethods.push_back("RK Forward Euler");
48 RKMethods.push_back("RK Explicit 4 Stage");
49 RKMethods.push_back("RK Explicit 3/8 Rule");
50 RKMethods.push_back("RK Explicit 4 Stage 3rd order by Runge");
51 RKMethods.push_back("RK Explicit 5 Stage 3rd order by Kinnmark and Gray");
52 RKMethods.push_back("RK Explicit 3 Stage 3rd order");
53 RKMethods.push_back("RK Explicit 3 Stage 3rd order TVD");
54 RKMethods.push_back("RK Explicit 3 Stage 3rd order by Heun");
55 RKMethods.push_back("RK Explicit Midpoint");
56 RKMethods.push_back("RK Explicit Trapezoidal");
57 RKMethods.push_back("Heuns Method");
58 RKMethods.push_back("Bogacki-Shampine 3(2) Pair");
59 RKMethods.push_back("Merson 4(5) Pair");
60
61 for (std::vector<std::string>::size_type m = 0; m != RKMethods.size(); m++) {
62 std::string RKMethod = RKMethods[m];
63 std::replace(RKMethod.begin(), RKMethod.end(), ' ', '_');
64 std::replace(RKMethod.begin(), RKMethod.end(), '/', '.');
65
66 // Read params from .xml file
67 RCP<ParameterList> pList =
68 getParametersFromXmlFile("Tempus_ExplicitRK_SinCos.xml");
69
70 // Setup the SinCosModel
71 RCP<ParameterList> scm_pl = sublist(pList, "SinCosModel", true);
72 auto model = rcp(new SinCosModel<double>(scm_pl));
73
74 // Set the Stepper
75 RCP<ParameterList> tempusPL = sublist(pList, "Tempus", true);
76 if (RKMethods[m] == "General ERK") {
77 tempusPL->sublist("Demo Integrator")
78 .set("Stepper Name", "Demo Stepper 2");
79 }
80 else {
81 tempusPL->sublist("Demo Stepper").set("Stepper Type", RKMethods[m]);
82 }
83
84 // Test constructor IntegratorBasic(tempusPL, model, runInitialize)
85 {
86 RCP<Tempus::IntegratorBasic<double>> integrator =
87 Tempus::createIntegratorBasic<double>(tempusPL, model, false);
88
89 // check initialization
90 {
91 // TSC should not be initialized
92 TEST_ASSERT(!integrator->getNonConstTimeStepControl()->isInitialized());
93 // now initialize everything (TSC should also be initilized)
94 integrator->initialize();
95 // TSC should be initialized
96 TEST_ASSERT(integrator->getNonConstTimeStepControl()->isInitialized());
97 }
98
99 RCP<ParameterList> stepperPL = sublist(tempusPL, "Demo Stepper", true);
100 if (RKMethods[m] == "General ERK")
101 stepperPL = sublist(tempusPL, "Demo Stepper 2", true);
102 RCP<ParameterList> defaultPL =
103 Teuchos::rcp_const_cast<Teuchos::ParameterList>(
104 integrator->getStepper()->getValidParameters());
105
106 bool pass = haveSameValuesSorted(*stepperPL, *defaultPL, true);
107 if (!pass) {
108 out << std::endl;
109 out << "stepperPL -------------- \n"
110 << *stepperPL << std::endl;
111 out << "defaultPL -------------- \n"
112 << *defaultPL << std::endl;
113 }
114 TEST_ASSERT(pass)
115 }
116
117 // Adjust tempusPL to default values.
118 if (RKMethods[m] == "General ERK") {
119 tempusPL->sublist("Demo Stepper 2")
120 .set("Initial Condition Consistency", "None");
121 }
122 if (RKMethods[m] == "RK Forward Euler" ||
123 RKMethods[m] == "Bogacki-Shampine 3(2) Pair") {
124 tempusPL->sublist("Demo Stepper")
125 .set("Initial Condition Consistency", "Consistent");
126 tempusPL->sublist("Demo Stepper").set("Use FSAL", true);
127 }
128 else {
129 tempusPL->sublist("Demo Stepper")
130 .set("Initial Condition Consistency", "None");
131 }
132
133 // Test constructor IntegratorBasic(model, stepperType)
134 {
135 RCP<Tempus::IntegratorBasic<double>> integrator =
136 Tempus::createIntegratorBasic<double>(model, RKMethods[m]);
137
138 RCP<ParameterList> stepperPL = sublist(tempusPL, "Demo Stepper", true);
139 if (RKMethods[m] == "General ERK")
140 stepperPL = sublist(tempusPL, "Demo Stepper 2", true);
141 RCP<ParameterList> defaultPL =
142 Teuchos::rcp_const_cast<Teuchos::ParameterList>(
143 integrator->getStepper()->getValidParameters());
144
145 bool pass = haveSameValuesSorted(*stepperPL, *defaultPL, true);
146 if (!pass) {
147 out << std::endl;
148 out << "stepperPL -------------- \n"
149 << *stepperPL << std::endl;
150 out << "defaultPL -------------- \n"
151 << *defaultPL << std::endl;
152 }
153 TEST_ASSERT(pass)
154 }
155 }
156}
157
158// ************************************************************
159// ************************************************************
160TEUCHOS_UNIT_TEST(ExplicitRK, ConstructingFromDefaults)
161{
162 double dt = 0.1;
163 std::vector<std::string> options;
164 options.push_back("Default Parameters");
165 options.push_back("ICConsistency and Check");
166
167 for (const auto& option : options) {
168 // Read params from .xml file
169 RCP<ParameterList> pList =
170 getParametersFromXmlFile("Tempus_ExplicitRK_SinCos.xml");
171 RCP<ParameterList> pl = sublist(pList, "Tempus", true);
172
173 // Setup the SinCosModel
174 RCP<ParameterList> scm_pl = sublist(pList, "SinCosModel", true);
175 // RCP<SinCosModel<double> > model = sineCosineModel(scm_pl);
176 auto model = rcp(new SinCosModel<double>(scm_pl));
177
178 // Setup Stepper for field solve ----------------------------
179 RCP<Tempus::StepperFactory<double>> sf =
180 Teuchos::rcp(new Tempus::StepperFactory<double>());
181 RCP<Tempus::Stepper<double>> stepper =
182 sf->createStepper("RK Explicit 4 Stage");
183 stepper->setModel(model);
184 if (option == "ICConsistency and Check") {
185 stepper->setICConsistency("Consistent");
186 stepper->setICConsistencyCheck(true);
187 }
188 stepper->initialize();
189
190 // Setup TimeStepControl ------------------------------------
191 auto timeStepControl = rcp(new Tempus::TimeStepControl<double>());
192 ParameterList tscPL =
193 pl->sublist("Demo Integrator").sublist("Time Step Control");
194 timeStepControl->setInitIndex(tscPL.get<int>("Initial Time Index"));
195 timeStepControl->setInitTime(tscPL.get<double>("Initial Time"));
196 timeStepControl->setFinalTime(tscPL.get<double>("Final Time"));
197 timeStepControl->setInitTimeStep(dt);
198 timeStepControl->initialize();
199
200 // Setup initial condition SolutionState --------------------
201 auto inArgsIC = model->getNominalValues();
202 auto icSolution =
203 rcp_const_cast<Thyra::VectorBase<double>>(inArgsIC.get_x());
204 auto icState = Tempus::createSolutionStateX(icSolution);
205 icState->setTime(timeStepControl->getInitTime());
206 icState->setIndex(timeStepControl->getInitIndex());
207 icState->setTimeStep(0.0);
208 icState->setOrder(stepper->getOrder());
209 icState->setSolutionStatus(Tempus::Status::PASSED); // ICs are passing.
210
211 // Setup SolutionHistory ------------------------------------
212 auto solutionHistory = rcp(new Tempus::SolutionHistory<double>());
213 solutionHistory->setName("Forward States");
214 solutionHistory->setStorageType(Tempus::STORAGE_TYPE_STATIC);
215 solutionHistory->setStorageLimit(2);
216 solutionHistory->addState(icState);
217
218 // Setup Integrator -----------------------------------------
219 RCP<Tempus::IntegratorBasic<double>> integrator =
220 Tempus::createIntegratorBasic<double>();
221 integrator->setStepper(stepper);
222 integrator->setTimeStepControl(timeStepControl);
223 integrator->setSolutionHistory(solutionHistory);
224 // integrator->setObserver(...);
225 integrator->initialize();
226
227 // Integrate to timeMax
228 bool integratorStatus = integrator->advanceTime();
229 TEST_ASSERT(integratorStatus)
230
231 // Test if at 'Final Time'
232 double time = integrator->getTime();
233 double timeFinal = pl->sublist("Demo Integrator")
234 .sublist("Time Step Control")
235 .get<double>("Final Time");
236 TEST_FLOATING_EQUALITY(time, timeFinal, 1.0e-14);
237
238 // Time-integrated solution and the exact solution
239 RCP<Thyra::VectorBase<double>> x = integrator->getX();
240 RCP<const Thyra::VectorBase<double>> x_exact =
241 model->getExactSolution(time).get_x();
242
243 // Calculate the error
244 RCP<Thyra::VectorBase<double>> xdiff = x->clone_v();
245 Thyra::V_StVpStV(xdiff.ptr(), 1.0, *x_exact, -1.0, *(x));
246
247 // Check the order and intercept
248 out << " Stepper = " << stepper->description() << "\n with "
249 << option << std::endl;
250 out << " =========================" << std::endl;
251 out << " Exact solution : " << get_ele(*(x_exact), 0) << " "
252 << get_ele(*(x_exact), 1) << std::endl;
253 out << " Computed solution: " << get_ele(*(x), 0) << " "
254 << get_ele(*(x), 1) << std::endl;
255 out << " Difference : " << get_ele(*(xdiff), 0) << " "
256 << get_ele(*(xdiff), 1) << std::endl;
257 out << " =========================" << std::endl;
258 TEST_FLOATING_EQUALITY(get_ele(*(x), 0), 0.841470, 1.0e-4);
259 TEST_FLOATING_EQUALITY(get_ele(*(x), 1), 0.540303, 1.0e-4);
260 }
261}
262
263// ************************************************************
264// ************************************************************
265TEUCHOS_UNIT_TEST(ExplicitRK, useFSAL_false)
266{
267 double dt = 0.1;
268 std::vector<std::string> RKMethods;
269 RKMethods.push_back("RK Forward Euler");
270 RKMethods.push_back("Bogacki-Shampine 3(2) Pair");
271
272 for (std::vector<std::string>::size_type m = 0; m != RKMethods.size(); m++) {
273 // Setup the SinCosModel ------------------------------------
274 auto model = rcp(new SinCosModel<double>());
275
276 // Setup Stepper for field solve ----------------------------
277 auto sf = Teuchos::rcp(new Tempus::StepperFactory<double>());
278 auto stepper = sf->createStepper(RKMethods[m]);
279 stepper->setModel(model);
280 stepper->setUseFSAL(false);
281 stepper->initialize();
282
283 // Setup TimeStepControl ------------------------------------
284 auto timeStepControl = rcp(new Tempus::TimeStepControl<double>());
285 timeStepControl->setInitTime(0.0);
286 timeStepControl->setFinalTime(1.0);
287 timeStepControl->setInitTimeStep(dt);
288 timeStepControl->initialize();
289
290 // Setup initial condition SolutionState --------------------
291 auto inArgsIC = model->getNominalValues();
292 auto icSolution =
293 rcp_const_cast<Thyra::VectorBase<double>>(inArgsIC.get_x());
294 auto icState = Tempus::createSolutionStateX(icSolution);
295 icState->setTime(timeStepControl->getInitTime());
296 icState->setIndex(timeStepControl->getInitIndex());
297 icState->setTimeStep(0.0);
298 icState->setOrder(stepper->getOrder());
299 icState->setSolutionStatus(Tempus::Status::PASSED); // ICs are passing.
300
301 // Setup SolutionHistory ------------------------------------
302 auto solutionHistory = rcp(new Tempus::SolutionHistory<double>());
303 solutionHistory->setName("Forward States");
304 solutionHistory->setStorageType(Tempus::STORAGE_TYPE_STATIC);
305 solutionHistory->setStorageLimit(2);
306 solutionHistory->addState(icState);
307
308 // Setup Integrator -----------------------------------------
309 auto integrator = Tempus::createIntegratorBasic<double>();
310 integrator->setStepper(stepper);
311 integrator->setTimeStepControl(timeStepControl);
312 integrator->setSolutionHistory(solutionHistory);
313 integrator->initialize();
314
315 // Integrate to timeMax
316 bool integratorStatus = integrator->advanceTime();
317 TEST_ASSERT(integratorStatus)
318
319 // Test if at 'Final Time'
320 double time = integrator->getTime();
321 TEST_FLOATING_EQUALITY(time, 1.0, 1.0e-14);
322
323 // Time-integrated solution and the exact solution
324 auto x = integrator->getX();
325 auto x_exact = model->getExactSolution(time).get_x();
326
327 // Calculate the error
328 RCP<Thyra::VectorBase<double>> xdiff = x->clone_v();
329 Thyra::V_StVpStV(xdiff.ptr(), 1.0, *x_exact, -1.0, *(x));
330
331 // Check the order and intercept
332 out << " Stepper = " << stepper->description() << "\n with "
333 << "useFSAL=false" << std::endl;
334 out << " =========================" << std::endl;
335 out << " Exact solution : " << get_ele(*(x_exact), 0) << " "
336 << get_ele(*(x_exact), 1) << std::endl;
337 out << " Computed solution: " << get_ele(*(x), 0) << " "
338 << get_ele(*(x), 1) << std::endl;
339 out << " Difference : " << get_ele(*(xdiff), 0) << " "
340 << get_ele(*(xdiff), 1) << std::endl;
341 out << " =========================" << std::endl;
342 if (RKMethods[m] == "RK Forward Euler") {
343 TEST_FLOATING_EQUALITY(get_ele(*(x), 0), 0.882508, 1.0e-4);
344 TEST_FLOATING_EQUALITY(get_ele(*(x), 1), 0.570790, 1.0e-4);
345 }
346 else if (RKMethods[m] == "Bogacki-Shampine 3(2) Pair") {
347 TEST_FLOATING_EQUALITY(get_ele(*(x), 0), 0.841438, 1.0e-4);
348 TEST_FLOATING_EQUALITY(get_ele(*(x), 1), 0.540277, 1.0e-4);
349 }
350 }
351}
352
353// ************************************************************
354// ************************************************************
355TEUCHOS_UNIT_TEST(ExplicitRK, SinCos)
356{
357 std::vector<std::string> RKMethods;
358 RKMethods.push_back("General ERK");
359 RKMethods.push_back("RK Forward Euler");
360 RKMethods.push_back("RK Explicit 4 Stage");
361 RKMethods.push_back("RK Explicit 3/8 Rule");
362 RKMethods.push_back("RK Explicit 4 Stage 3rd order by Runge");
363 RKMethods.push_back("RK Explicit 5 Stage 3rd order by Kinnmark and Gray");
364 RKMethods.push_back("RK Explicit 3 Stage 3rd order");
365 RKMethods.push_back("RK Explicit 3 Stage 3rd order TVD");
366 RKMethods.push_back("RK Explicit 3 Stage 3rd order by Heun");
367 RKMethods.push_back("RK Explicit Midpoint");
368 RKMethods.push_back("RK Explicit Trapezoidal");
369 RKMethods.push_back("Heuns Method");
370 RKMethods.push_back("Bogacki-Shampine 3(2) Pair");
371 RKMethods.push_back("Merson 4(5) Pair"); // slope = 3.87816
372 RKMethods.push_back("General ERK Embedded");
373 RKMethods.push_back("SSPERK22");
374 RKMethods.push_back("SSPERK33");
375 RKMethods.push_back("SSPERK54"); // slope = 3.94129
376 RKMethods.push_back("RK2");
377
378 std::vector<double> RKMethodErrors;
379 RKMethodErrors.push_back(8.33251e-07);
380 RKMethodErrors.push_back(0.051123);
381 RKMethodErrors.push_back(8.33251e-07);
382 RKMethodErrors.push_back(8.33251e-07);
383 RKMethodErrors.push_back(4.16897e-05);
384 RKMethodErrors.push_back(8.32108e-06);
385 RKMethodErrors.push_back(4.16603e-05);
386 RKMethodErrors.push_back(4.16603e-05);
387 RKMethodErrors.push_back(4.16603e-05);
388 RKMethodErrors.push_back(0.00166645);
389 RKMethodErrors.push_back(0.00166645);
390 RKMethodErrors.push_back(0.00166645);
391 RKMethodErrors.push_back(4.16603e-05);
392 RKMethodErrors.push_back(1.39383e-07);
393 RKMethodErrors.push_back(4.16603e-05);
394 RKMethodErrors.push_back(0.00166645);
395 RKMethodErrors.push_back(4.16603e-05);
396 RKMethodErrors.push_back(3.85613e-07); // SSPERK54
397 RKMethodErrors.push_back(0.00166644);
398
399 TEST_ASSERT(RKMethods.size() == RKMethodErrors.size());
400
401 for (std::vector<std::string>::size_type m = 0; m != RKMethods.size(); m++) {
402 std::string RKMethod = RKMethods[m];
403 std::replace(RKMethod.begin(), RKMethod.end(), ' ', '_');
404 std::replace(RKMethod.begin(), RKMethod.end(), '/', '.');
405
406 RCP<Tempus::IntegratorBasic<double>> integrator;
407 std::vector<RCP<Thyra::VectorBase<double>>> solutions;
408 std::vector<RCP<Thyra::VectorBase<double>>> solutionsDot;
409 std::vector<double> StepSize;
410 std::vector<double> xErrorNorm;
411 std::vector<double> xDotErrorNorm;
412
413 const int nTimeStepSizes = 7;
414 double dt = 0.2;
415 double time = 0.0;
416 for (int n = 0; n < nTimeStepSizes; n++) {
417 // Read params from .xml file
418 RCP<ParameterList> pList =
419 getParametersFromXmlFile("Tempus_ExplicitRK_SinCos.xml");
420
421 // Setup the SinCosModel
422 RCP<ParameterList> scm_pl = sublist(pList, "SinCosModel", true);
423 // RCP<SinCosModel<double> > model = sineCosineModel(scm_pl);
424 auto model = rcp(new SinCosModel<double>(scm_pl));
425
426 // Set the Stepper
427 RCP<ParameterList> pl = sublist(pList, "Tempus", true);
428 if (RKMethods[m] == "General ERK") {
429 pl->sublist("Demo Integrator").set("Stepper Name", "Demo Stepper 2");
430 }
431 else if (RKMethods[m] == "General ERK Embedded") {
432 pl->sublist("Demo Integrator")
433 .set("Stepper Name", "General ERK Embedded Stepper");
434 }
435 else {
436 pl->sublist("Demo Stepper").set("Stepper Type", RKMethods[m]);
437 }
438
439 dt /= 2;
440
441 // Setup the Integrator and reset initial time step
442 pl->sublist("Demo Integrator")
443 .sublist("Time Step Control")
444 .set("Initial Time Step", dt);
445 integrator = Tempus::createIntegratorBasic<double>(pl, model);
446
447 // Initial Conditions
448 // During the Integrator construction, the initial SolutionState
449 // is set by default to model->getNominalVales().get_x(). However,
450 // the application can set it also by
451 // integrator->initializeSolutionHistory.
452 RCP<Thyra::VectorBase<double>> x0 =
453 model->getNominalValues().get_x()->clone_v();
454 integrator->initializeSolutionHistory(0.0, x0);
455
456 // Integrate to timeMax
457 bool integratorStatus = integrator->advanceTime();
458 TEST_ASSERT(integratorStatus)
459
460 // Test if at 'Final Time'
461 time = integrator->getTime();
462 double timeFinal = pl->sublist("Demo Integrator")
463 .sublist("Time Step Control")
464 .get<double>("Final Time");
465 TEST_FLOATING_EQUALITY(time, timeFinal, 1.0e-14);
466
467 // Time-integrated solution and the exact solution
468 RCP<Thyra::VectorBase<double>> x = integrator->getX();
469 RCP<const Thyra::VectorBase<double>> x_exact =
470 model->getExactSolution(time).get_x();
471
472 // Plot sample solution and exact solution
473 if (n == 0) {
474 RCP<const SolutionHistory<double>> solutionHistory =
475 integrator->getSolutionHistory();
476 writeSolution("Tempus_" + RKMethod + "_SinCos.dat", solutionHistory);
477
478 auto solnHistExact = rcp(new Tempus::SolutionHistory<double>());
479 for (int i = 0; i < solutionHistory->getNumStates(); i++) {
480 double time_i = (*solutionHistory)[i]->getTime();
481 auto state = Tempus::createSolutionStateX(
482 rcp_const_cast<Thyra::VectorBase<double>>(
483 model->getExactSolution(time_i).get_x()),
484 rcp_const_cast<Thyra::VectorBase<double>>(
485 model->getExactSolution(time_i).get_x_dot()));
486 state->setTime((*solutionHistory)[i]->getTime());
487 solnHistExact->addState(state);
488 }
489 writeSolution("Tempus_" + RKMethod + "_SinCos-Ref.dat", solnHistExact);
490 }
491
492 // Store off the final solution and step size
493 StepSize.push_back(dt);
494 auto solution = Thyra::createMember(model->get_x_space());
495 Thyra::copy(*(integrator->getX()), solution.ptr());
496 solutions.push_back(solution);
497 auto solutionDot = Thyra::createMember(model->get_x_space());
498 Thyra::copy(*(integrator->getXDot()), solutionDot.ptr());
499 solutionsDot.push_back(solutionDot);
500 if (n == nTimeStepSizes - 1) { // Add exact solution last in vector.
501 StepSize.push_back(0.0);
502 auto solutionExact = Thyra::createMember(model->get_x_space());
503 Thyra::copy(*(model->getExactSolution(time).get_x()),
504 solutionExact.ptr());
505 solutions.push_back(solutionExact);
506 auto solutionDotExact = Thyra::createMember(model->get_x_space());
507 Thyra::copy(*(model->getExactSolution(time).get_x_dot()),
508 solutionDotExact.ptr());
509 solutionsDot.push_back(solutionDotExact);
510 }
511 }
512
513 // Check the order and intercept
514 double xSlope = 0.0;
515 double xDotSlope = 0.0;
516 RCP<Tempus::Stepper<double>> stepper = integrator->getStepper();
517 double order = stepper->getOrder();
518 writeOrderError("Tempus_" + RKMethod + "_SinCos-Error.dat", stepper,
519 StepSize, solutions, xErrorNorm, xSlope, solutionsDot,
520 xDotErrorNorm, xDotSlope, out);
521
522 double order_tol = 0.01;
523 if (RKMethods[m] == "Merson 4(5) Pair") order_tol = 0.04;
524 if (RKMethods[m] == "SSPERK54") order_tol = 0.06;
525
526 TEST_FLOATING_EQUALITY(xSlope, order, order_tol);
527 TEST_FLOATING_EQUALITY(xErrorNorm[0], RKMethodErrors[m], 1.0e-4);
528 // xDot not yet available for ExplicitRK methods.
529 // TEST_FLOATING_EQUALITY( xDotSlope, order, 0.01 );
530 // TEST_FLOATING_EQUALITY( xDotErrorNorm[0], 0.0486418, 1.0e-4 );
531 }
532 // Teuchos::TimeMonitor::summarize();
533}
534
535// ************************************************************
536// ************************************************************
537TEUCHOS_UNIT_TEST(ExplicitRK, EmbeddedVanDerPol)
538{
539 std::vector<std::string> IntegratorList;
540 IntegratorList.push_back("Embedded_Integrator_PID");
541 IntegratorList.push_back("Demo_Integrator");
542 IntegratorList.push_back("Embedded_Integrator");
543 IntegratorList.push_back("General_Embedded_Integrator");
544 IntegratorList.push_back("Embedded_Integrator_PID_General");
545
546 // the embedded solution will test the following:
547 // using the starting stepsize routine, this has now decreased
548 const int refIstep = 36;
549
550 for (auto integratorChoice : IntegratorList) {
551 out << "Using Integrator: " << integratorChoice << " !!!" << std::endl;
552
553 // Read params from .xml file
554 RCP<ParameterList> pList =
555 getParametersFromXmlFile("Tempus_ExplicitRK_VanDerPol.xml");
556
557 // Setup the VanDerPolModel
558 RCP<ParameterList> vdpm_pl = sublist(pList, "VanDerPolModel", true);
559 auto model = rcp(new VanDerPolModel<double>(vdpm_pl));
560
561 // Set the Integrator and Stepper
562 RCP<ParameterList> pl = sublist(pList, "Tempus", true);
563 pl->set("Integrator Name", integratorChoice);
564
565 // Setup the Integrator
566 RCP<Tempus::IntegratorBasic<double>> integrator =
567 Tempus::createIntegratorBasic<double>(pl, model);
568
569 const std::string RKMethod =
570 pl->sublist(integratorChoice).get<std::string>("Stepper Name");
571
572 // Integrate to timeMax
573 bool integratorStatus = integrator->advanceTime();
574 TEST_ASSERT(integratorStatus);
575
576 // Test if at 'Final Time'
577 double time = integrator->getTime();
578 double timeFinal = pl->sublist(integratorChoice)
579 .sublist("Time Step Control")
580 .get<double>("Final Time");
581 TEST_FLOATING_EQUALITY(time, timeFinal, 1.0e-14);
582
583 // Numerical reference solution at timeFinal (for \epsilon = 0.1)
584 RCP<Thyra::VectorBase<double>> x = integrator->getX();
585 RCP<Thyra::VectorBase<double>> xref = x->clone_v();
586 Thyra::set_ele(0, -1.5484458614405929, xref.ptr());
587 Thyra::set_ele(1, 1.0181127316101317, xref.ptr());
588
589 // Calculate the error
590 RCP<Thyra::VectorBase<double>> xdiff = x->clone_v();
591 Thyra::V_StVpStV(xdiff.ptr(), 1.0, *xref, -1.0, *(x));
592 const double L2norm = Thyra::norm_2(*xdiff);
593
594 // Test number of steps, failures, and accuracy
595 if ((integratorChoice == "Embedded_Integrator_PID") ||
596 (integratorChoice == "Embedded_Integrator_PID_General")) {
597 const double absTol = pl->sublist(integratorChoice)
598 .sublist("Time Step Control")
599 .get<double>("Maximum Absolute Error");
600 const double relTol = pl->sublist(integratorChoice)
601 .sublist("Time Step Control")
602 .get<double>("Maximum Relative Error");
603
604 // Should be close to the prescribed tolerance (magnitude)
605 TEST_COMPARE(std::log10(L2norm), <, std::log10(absTol));
606 TEST_COMPARE(std::log10(L2norm), <, std::log10(relTol));
607
608 // get the number of time steps and number of step failure
609 // const int nFailure_c = integrator->getSolutionHistory()->
610 // getCurrentState()->getMetaData()->getNFailures();
611 const int iStep =
612 integrator->getSolutionHistory()->getCurrentState()->getIndex();
613 const int nFail = integrator->getSolutionHistory()
614 ->getCurrentState()
615 ->getMetaData()
616 ->getNRunningFailures();
617
618 // test for number of steps
619 TEST_EQUALITY(iStep, refIstep);
620 out << "Tolerance = " << absTol << " L2norm = " << L2norm
621 << " iStep = " << iStep << " nFail = " << nFail << std::endl;
622 }
623
624 // Plot sample solution and exact solution
625 std::ofstream ftmp("Tempus_" + integratorChoice + RKMethod +
626 "_VDP_Example.dat");
627 RCP<const SolutionHistory<double>> solutionHistory =
628 integrator->getSolutionHistory();
629 int nStates = solutionHistory->getNumStates();
630 // RCP<const Thyra::VectorBase<double> > x_exact_plot;
631 for (int i = 0; i < nStates; i++) {
632 RCP<const SolutionState<double>> solutionState = (*solutionHistory)[i];
633 double time_i = solutionState->getTime();
634 RCP<const Thyra::VectorBase<double>> x_plot = solutionState->getX();
635 // x_exact_plot = model->getExactSolution(time_i).get_x();
636 ftmp << time_i << " " << Thyra::get_ele(*(x_plot), 0) << " "
637 << Thyra::get_ele(*(x_plot), 1) << " " << std::endl;
638 }
639 ftmp.close();
640 }
641
642 Teuchos::TimeMonitor::summarize();
643}
644
645// ************************************************************
646// ************************************************************
647TEUCHOS_UNIT_TEST(ExplicitRK, stage_number)
648{
649 double dt = 0.1;
650
651 // Read params from .xml file
652 RCP<ParameterList> pList =
653 getParametersFromXmlFile("Tempus_ExplicitRK_SinCos.xml");
654 RCP<ParameterList> pl = sublist(pList, "Tempus", true);
655
656 // Setup the SinCosModel
657 RCP<ParameterList> scm_pl = sublist(pList, "SinCosModel", true);
658 // RCP<SinCosModel<double> > model = sineCosineModel(scm_pl);
659 auto model = rcp(new SinCosModel<double>(scm_pl));
660
661 // Setup Stepper for field solve ----------------------------
662 RCP<Tempus::StepperFactory<double>> sf =
663 Teuchos::rcp(new Tempus::StepperFactory<double>());
664 RCP<Tempus::Stepper<double>> stepper =
665 sf->createStepper("RK Explicit 4 Stage");
666 stepper->setModel(model);
667 stepper->initialize();
668
669 // Setup TimeStepControl ------------------------------------
670 auto timeStepControl = rcp(new Tempus::TimeStepControl<double>());
671 ParameterList tscPL =
672 pl->sublist("Demo Integrator").sublist("Time Step Control");
673 timeStepControl->setInitIndex(tscPL.get<int>("Initial Time Index"));
674 timeStepControl->setInitTime(tscPL.get<double>("Initial Time"));
675 timeStepControl->setFinalTime(tscPL.get<double>("Final Time"));
676 timeStepControl->setInitTimeStep(dt);
677 timeStepControl->initialize();
678
679 // Setup initial condition SolutionState --------------------
680 auto inArgsIC = model->getNominalValues();
681 auto icSolution = rcp_const_cast<Thyra::VectorBase<double>>(inArgsIC.get_x());
682 auto icState = Tempus::createSolutionStateX(icSolution);
683 icState->setTime(timeStepControl->getInitTime());
684 icState->setIndex(timeStepControl->getInitIndex());
685 icState->setTimeStep(0.0);
686 icState->setOrder(stepper->getOrder());
687 icState->setSolutionStatus(Tempus::Status::PASSED); // ICs are passing.
688
689 // Setup SolutionHistory ------------------------------------
690 auto solutionHistory = rcp(new Tempus::SolutionHistory<double>());
691 solutionHistory->setName("Forward States");
692 solutionHistory->setStorageType(Tempus::STORAGE_TYPE_STATIC);
693 solutionHistory->setStorageLimit(2);
694 solutionHistory->addState(icState);
695
696 // Setup Integrator -----------------------------------------
697 RCP<Tempus::IntegratorBasic<double>> integrator =
698 Tempus::createIntegratorBasic<double>();
699 integrator->setStepper(stepper);
700 integrator->setTimeStepControl(timeStepControl);
701 integrator->setSolutionHistory(solutionHistory);
702 // integrator->setObserver(...);
703 integrator->initialize();
704
705 // get the RK stepper
706 auto erk_stepper =
707 Teuchos::rcp_dynamic_cast<Tempus::StepperExplicitRK<double>>(stepper,
708 true);
709
710 TEST_EQUALITY(-1, erk_stepper->getStageNumber());
711 const std::vector<int> ref_stageNumber = {1, 4, 8, 10, 11, -1, 5};
712 for (auto stage_number : ref_stageNumber) {
713 // set the stage number
714 erk_stepper->setStageNumber(stage_number);
715 // make sure we are getting the correct stage number
716 TEST_EQUALITY(stage_number, erk_stepper->getStageNumber());
717 }
718}
719
720} // 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::SinCosModel
Sine-Cosine model problem from Rythmos. This is a canonical Sine-Cosine differential equation.
Definition SinCosModel_decl.hpp:110
Tempus_Test::VanDerPolModel
van der Pol model problem for nonlinear electrical circuit.
Definition VanDerPolModel_decl.hpp:113
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