Tempus Version of the Day
Time Integration
Loading...
Searching...
No Matches
Tempus_DIRKTest.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 "Thyra_VectorStdOps.hpp"
15
16#include "Tempus_IntegratorBasic.hpp"
17#include "Tempus_StepperFactory.hpp"
18
19#include "../TestModels/SinCosModel.hpp"
20#include "../TestModels/VanDerPolModel.hpp"
21#include "../TestUtils/Tempus_ConvergenceTestUtils.hpp"
22
23#include <fstream>
24#include <vector>
25
26namespace Tempus_Test {
27
28using Teuchos::getParametersFromXmlFile;
29using Teuchos::ParameterList;
30using Teuchos::parameterList;
31using Teuchos::RCP;
32using Teuchos::rcp;
33using Teuchos::rcp_const_cast;
34using Teuchos::rcp_dynamic_cast;
35using Teuchos::sublist;
36
37using Tempus::IntegratorBasic;
38using Tempus::SolutionHistory;
39using Tempus::SolutionState;
40
41// ************************************************************
42// ************************************************************
43TEUCHOS_UNIT_TEST(DIRK, ParameterList)
44{
45 std::vector<std::string> RKMethods;
46 RKMethods.push_back("General DIRK");
47 RKMethods.push_back("RK Backward Euler");
48 RKMethods.push_back("DIRK 1 Stage Theta Method");
49 RKMethods.push_back("RK Implicit 1 Stage 1st order Radau IA");
50 RKMethods.push_back("RK Implicit Midpoint");
51 RKMethods.push_back("SDIRK 2 Stage 2nd order");
52 RKMethods.push_back("RK Implicit 2 Stage 2nd order Lobatto IIIB");
53 RKMethods.push_back("SDIRK 2 Stage 3rd order");
54 RKMethods.push_back("EDIRK 2 Stage 3rd order");
55 RKMethods.push_back("EDIRK 2 Stage Theta Method");
56 RKMethods.push_back("SDIRK 3 Stage 4th order");
57 RKMethods.push_back("SDIRK 5 Stage 4th order");
58 RKMethods.push_back("SDIRK 5 Stage 5th order");
59 RKMethods.push_back("SDIRK 2(1) Pair");
60 RKMethods.push_back("RK Trapezoidal Rule");
61 RKMethods.push_back("RK Crank-Nicolson");
62
63 for (std::vector<std::string>::size_type m = 0; m != RKMethods.size(); m++) {
64 std::string RKMethod = RKMethods[m];
65
66 // Read params from .xml file
67 RCP<ParameterList> pList =
68 getParametersFromXmlFile("Tempus_DIRK_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 RCP<ParameterList> tempusPL = sublist(pList, "Tempus", true);
75 tempusPL->sublist("Default Stepper").set("Stepper Type", RKMethods[m]);
76
77 if (RKMethods[m] == "DIRK 1 Stage Theta Method" ||
78 RKMethods[m] == "EDIRK 2 Stage Theta Method") {
79 // Construct in the same order as default.
80 RCP<ParameterList> stepperPL = sublist(tempusPL, "Default Stepper", true);
81 RCP<ParameterList> solverPL = parameterList();
82 *solverPL = *(sublist(stepperPL, "Default Solver", true));
83 if (RKMethods[m] == "EDIRK 2 Stage Theta Method")
84 tempusPL->sublist("Default Stepper").set<bool>("Use FSAL", 1);
85 tempusPL->sublist("Default Stepper").remove("Zero Initial Guess");
86 tempusPL->sublist("Default Stepper").remove("Default Solver");
87 tempusPL->sublist("Default Stepper").set<bool>("Zero Initial Guess", 0);
88 tempusPL->sublist("Default Stepper").remove("Reset Initial Guess");
89 tempusPL->sublist("Default Stepper").set<bool>("Reset Initial Guess", 1);
90 tempusPL->sublist("Default Stepper").set("Default Solver", *solverPL);
91 tempusPL->sublist("Default Stepper").set<double>("theta", 0.5);
92 }
93 else if (RKMethods[m] == "SDIRK 2 Stage 2nd order") {
94 // Construct in the same order as default.
95 RCP<ParameterList> stepperPL = sublist(tempusPL, "Default Stepper", true);
96 RCP<ParameterList> solverPL = parameterList();
97 *solverPL = *(sublist(stepperPL, "Default Solver", true));
98 tempusPL->sublist("Default Stepper").remove("Default Solver");
99 tempusPL->sublist("Default Stepper").remove("Zero Initial Guess");
100 tempusPL->sublist("Default Stepper").set<bool>("Zero Initial Guess", 0);
101 tempusPL->sublist("Default Stepper").remove("Reset Initial Guess");
102 tempusPL->sublist("Default Stepper").set<bool>("Reset Initial Guess", 1);
103 tempusPL->sublist("Default Stepper").set("Default Solver", *solverPL);
104 tempusPL->sublist("Default Stepper")
105 .set<double>("gamma", 0.2928932188134524);
106 }
107 else if (RKMethods[m] == "SDIRK 2 Stage 3rd order") {
108 // Construct in the same order as default.
109 RCP<ParameterList> stepperPL = sublist(tempusPL, "Default Stepper", true);
110 RCP<ParameterList> solverPL = parameterList();
111 *solverPL = *(sublist(stepperPL, "Default Solver", true));
112 tempusPL->sublist("Default Stepper").remove("Zero Initial Guess");
113 tempusPL->sublist("Default Stepper").set<bool>("Zero Initial Guess", 0);
114 tempusPL->sublist("Default Stepper").remove("Default Solver");
115 tempusPL->sublist("Default Stepper").remove("Reset Initial Guess");
116 tempusPL->sublist("Default Stepper").set<bool>("Reset Initial Guess", 1);
117 tempusPL->sublist("Default Stepper").set("Default Solver", *solverPL);
118 tempusPL->sublist("Default Stepper")
119 .set<std::string>("Gamma Type", "3rd Order A-stable");
120 tempusPL->sublist("Default Stepper")
121 .set<double>("gamma", 0.7886751345948128);
122 }
123 else if (RKMethods[m] == "RK Trapezoidal Rule") {
124 tempusPL->sublist("Default Stepper").set<bool>("Use FSAL", 1);
125 }
126 else if (RKMethods[m] == "RK Crank-Nicolson") {
127 tempusPL->sublist("Default Stepper").set<bool>("Use FSAL", 1);
128 // Match default Stepper Type
129 tempusPL->sublist("Default Stepper")
130 .set("Stepper Type", "RK Trapezoidal Rule");
131 }
132 else if (RKMethods[m] == "General DIRK") {
133 // Add the default tableau.
134 Teuchos::RCP<Teuchos::ParameterList> tableauPL = Teuchos::parameterList();
135 tableauPL->set<std::string>(
136 "A", "0.292893218813452 0; 0.707106781186548 0.292893218813452");
137 tableauPL->set<std::string>("b", "0.707106781186548 0.292893218813452");
138 tableauPL->set<std::string>("c", "0.292893218813452 1");
139 tableauPL->set<int>("order", 2);
140 tableauPL->set<std::string>("bstar", "");
141 tempusPL->sublist("Default Stepper").set("Tableau", *tableauPL);
142 }
143
144 // Test constructor IntegratorBasic(tempusPL, model)
145 {
146 RCP<Tempus::IntegratorBasic<double>> integrator =
147 Tempus::createIntegratorBasic<double>(tempusPL, model);
148
149 RCP<ParameterList> stepperPL = sublist(tempusPL, "Default Stepper", true);
150 RCP<ParameterList> defaultPL =
151 Teuchos::rcp_const_cast<Teuchos::ParameterList>(
152 integrator->getStepper()->getValidParameters());
153
154 // Do not worry about "Description" as it is documentation.
155 defaultPL->remove("Description");
156
157 bool pass = haveSameValuesSorted(*stepperPL, *defaultPL, true);
158 if (!pass) {
159 out << std::endl;
160 out << "stepperPL -------------- \n"
161 << *stepperPL << std::endl;
162 out << "defaultPL -------------- \n"
163 << *defaultPL << std::endl;
164 }
165 TEST_ASSERT(pass)
166 }
167
168 // Test constructor IntegratorBasic(model, stepperType)
169 {
170 RCP<Tempus::IntegratorBasic<double>> integrator =
171 Tempus::createIntegratorBasic<double>(model, RKMethods[m]);
172
173 RCP<ParameterList> stepperPL = sublist(tempusPL, "Default Stepper", true);
174 RCP<ParameterList> defaultPL =
175 Teuchos::rcp_const_cast<Teuchos::ParameterList>(
176 integrator->getStepper()->getValidParameters());
177
178 // Do not worry about "Description" as it is documentation.
179 defaultPL->remove("Description");
180
181 // These Steppers have 'Initial Condition Consistency = Consistent'
182 // set as the default, so the ParameterList has been modified by
183 // NOX (filled in default parameters). Thus solver comparison
184 // will be removed.
185 if (RKMethods[m] == "EDIRK 2 Stage Theta Method" ||
186 RKMethods[m] == "RK Trapezoidal Rule" ||
187 RKMethods[m] == "RK Crank-Nicolson") {
188 stepperPL->set<std::string>("Initial Condition Consistency",
189 "Consistent");
190 stepperPL->remove("Default Solver");
191 defaultPL->remove("Default Solver");
192 }
193
194 bool pass = haveSameValuesSorted(*stepperPL, *defaultPL, true);
195 if (!pass) {
196 out << std::endl;
197 out << "stepperPL -------------- \n"
198 << *stepperPL << std::endl;
199 out << "defaultPL -------------- \n"
200 << *defaultPL << std::endl;
201 }
202 TEST_ASSERT(pass)
203 }
204 }
205}
206
207// ************************************************************
208// ************************************************************
209TEUCHOS_UNIT_TEST(DIRK, ConstructingFromDefaults)
210{
211 double dt = 0.025;
212 std::vector<std::string> options;
213 options.push_back("Default Parameters");
214 options.push_back("ICConsistency and Check");
215
216 for (const auto& option : options) {
217 // Read params from .xml file
218 RCP<ParameterList> pList =
219 getParametersFromXmlFile("Tempus_DIRK_SinCos.xml");
220 RCP<ParameterList> pl = sublist(pList, "Tempus", true);
221
222 // Setup the SinCosModel
223 RCP<ParameterList> scm_pl = sublist(pList, "SinCosModel", true);
224 // RCP<SinCosModel<double> > model = sineCosineModel(scm_pl);
225 auto model = rcp(new SinCosModel<double>(scm_pl));
226
227 // Setup Stepper for field solve ----------------------------
228 RCP<Tempus::StepperFactory<double>> sf =
229 Teuchos::rcp(new Tempus::StepperFactory<double>());
230 RCP<Tempus::Stepper<double>> stepper =
231 sf->createStepper("SDIRK 2 Stage 2nd order");
232 stepper->setModel(model);
233 if (option == "ICConsistency and Check") {
234 stepper->setICConsistency("Consistent");
235 stepper->setICConsistencyCheck(true);
236 }
237 stepper->initialize();
238
239 // Setup TimeStepControl ------------------------------------
240 auto timeStepControl = rcp(new Tempus::TimeStepControl<double>());
241 ParameterList tscPL =
242 pl->sublist("Default Integrator").sublist("Time Step Control");
243 timeStepControl->setInitIndex(tscPL.get<int>("Initial Time Index"));
244 timeStepControl->setInitTime(tscPL.get<double>("Initial Time"));
245 timeStepControl->setFinalTime(tscPL.get<double>("Final Time"));
246 timeStepControl->setInitTimeStep(dt);
247 timeStepControl->initialize();
248
249 // Setup initial condition SolutionState --------------------
250 auto inArgsIC = model->getNominalValues();
251 auto icSolution =
252 rcp_const_cast<Thyra::VectorBase<double>>(inArgsIC.get_x());
253 auto icState = Tempus::createSolutionStateX(icSolution);
254 icState->setTime(timeStepControl->getInitTime());
255 icState->setIndex(timeStepControl->getInitIndex());
256 icState->setTimeStep(0.0);
257 icState->setOrder(stepper->getOrder());
258 icState->setSolutionStatus(Tempus::Status::PASSED); // ICs are passing.
259
260 // Setup SolutionHistory ------------------------------------
261 auto solutionHistory = rcp(new Tempus::SolutionHistory<double>());
262 solutionHistory->setName("Forward States");
263 solutionHistory->setStorageType(Tempus::STORAGE_TYPE_STATIC);
264 solutionHistory->setStorageLimit(2);
265 solutionHistory->addState(icState);
266
267 // Setup Integrator -----------------------------------------
268 RCP<Tempus::IntegratorBasic<double>> integrator =
269 Tempus::createIntegratorBasic<double>();
270 integrator->setStepper(stepper);
271 integrator->setTimeStepControl(timeStepControl);
272 integrator->setSolutionHistory(solutionHistory);
273 // integrator->setObserver(...);
274 integrator->initialize();
275
276 // Integrate to timeMax
277 bool integratorStatus = integrator->advanceTime();
278 TEST_ASSERT(integratorStatus)
279
280 // Test if at 'Final Time'
281 double time = integrator->getTime();
282 double timeFinal = pl->sublist("Default Integrator")
283 .sublist("Time Step Control")
284 .get<double>("Final Time");
285 TEST_FLOATING_EQUALITY(time, timeFinal, 1.0e-14);
286
287 // Time-integrated solution and the exact solution
288 RCP<Thyra::VectorBase<double>> x = integrator->getX();
289 RCP<const Thyra::VectorBase<double>> x_exact =
290 model->getExactSolution(time).get_x();
291
292 // Calculate the error
293 RCP<Thyra::VectorBase<double>> xdiff = x->clone_v();
294 Thyra::V_StVpStV(xdiff.ptr(), 1.0, *x_exact, -1.0, *(x));
295
296 // Check the order and intercept
297 out << " Stepper = " << stepper->description() << " with " << option
298 << std::endl;
299 out << " =========================" << std::endl;
300 out << " Exact solution : " << get_ele(*(x_exact), 0) << " "
301 << get_ele(*(x_exact), 1) << std::endl;
302 out << " Computed solution: " << get_ele(*(x), 0) << " "
303 << get_ele(*(x), 1) << std::endl;
304 out << " Difference : " << get_ele(*(xdiff), 0) << " "
305 << get_ele(*(xdiff), 1) << std::endl;
306 out << " =========================" << std::endl;
307 TEST_FLOATING_EQUALITY(get_ele(*(x), 0), 0.841470, 1.0e-4);
308 TEST_FLOATING_EQUALITY(get_ele(*(x), 1), 0.540304, 1.0e-4);
309 }
310}
311
312// ************************************************************
313// ************************************************************
314TEUCHOS_UNIT_TEST(DIRK, useFSAL_false)
315{
316 double dt = 0.1;
317 std::vector<std::string> RKMethods;
318 RKMethods.push_back("EDIRK 2 Stage Theta Method");
319 RKMethods.push_back("RK Trapezoidal Rule");
320
321 for (std::vector<std::string>::size_type m = 0; m != RKMethods.size(); m++) {
322 // Setup the SinCosModel ------------------------------------
323 auto model = rcp(new SinCosModel<double>());
324
325 // Setup Stepper for field solve ----------------------------
326 auto sf = Teuchos::rcp(new Tempus::StepperFactory<double>());
327 auto stepper = sf->createStepper(RKMethods[m]);
328 stepper->setModel(model);
329 stepper->setUseFSAL(false);
330 stepper->initialize();
331
332 // Setup TimeStepControl ------------------------------------
333 auto timeStepControl = rcp(new Tempus::TimeStepControl<double>());
334 timeStepControl->setInitTime(0.0);
335 timeStepControl->setFinalTime(1.0);
336 timeStepControl->setInitTimeStep(dt);
337 timeStepControl->initialize();
338
339 // Setup initial condition SolutionState --------------------
340 auto inArgsIC = model->getNominalValues();
341 auto icSolution =
342 rcp_const_cast<Thyra::VectorBase<double>>(inArgsIC.get_x());
343 auto icState = Tempus::createSolutionStateX(icSolution);
344 icState->setTime(timeStepControl->getInitTime());
345 icState->setIndex(timeStepControl->getInitIndex());
346 icState->setTimeStep(0.0);
347 icState->setOrder(stepper->getOrder());
348 icState->setSolutionStatus(Tempus::Status::PASSED); // ICs are passing.
349
350 // Setup SolutionHistory ------------------------------------
351 auto solutionHistory = rcp(new Tempus::SolutionHistory<double>());
352 solutionHistory->setName("Forward States");
353 solutionHistory->setStorageType(Tempus::STORAGE_TYPE_STATIC);
354 solutionHistory->setStorageLimit(2);
355 solutionHistory->addState(icState);
356
357 // Setup Integrator -----------------------------------------
358 auto integrator = Tempus::createIntegratorBasic<double>();
359 integrator->setStepper(stepper);
360 integrator->setTimeStepControl(timeStepControl);
361 integrator->setSolutionHistory(solutionHistory);
362 integrator->initialize();
363
364 // Integrate to timeMax
365 bool integratorStatus = integrator->advanceTime();
366 TEST_ASSERT(integratorStatus)
367
368 // Test if at 'Final Time'
369 double time = integrator->getTime();
370 TEST_FLOATING_EQUALITY(time, 1.0, 1.0e-14);
371
372 // Time-integrated solution and the exact solution
373 auto x = integrator->getX();
374 auto x_exact = model->getExactSolution(time).get_x();
375
376 // Calculate the error
377 RCP<Thyra::VectorBase<double>> xdiff = x->clone_v();
378 Thyra::V_StVpStV(xdiff.ptr(), 1.0, *x_exact, -1.0, *(x));
379
380 // Check the order and intercept
381 out << " Stepper = " << stepper->description() << "\n with "
382 << "useFSAL=false" << std::endl;
383 out << " =========================" << std::endl;
384 out << " Exact solution : " << get_ele(*(x_exact), 0) << " "
385 << get_ele(*(x_exact), 1) << std::endl;
386 out << " Computed solution: " << get_ele(*(x), 0) << " "
387 << get_ele(*(x), 1) << std::endl;
388 out << " Difference : " << get_ele(*(xdiff), 0) << " "
389 << get_ele(*(xdiff), 1) << std::endl;
390 out << " =========================" << std::endl;
391 if (RKMethods[m] == "EDIRK 2 Stage Theta Method") {
392 TEST_FLOATING_EQUALITY(get_ele(*(x), 0), 0.841021, 1.0e-4);
393 TEST_FLOATING_EQUALITY(get_ele(*(x), 1), 0.541002, 1.0e-4);
394 }
395 else if (RKMethods[m] == "RK Trapezoidal Rule") {
396 TEST_FLOATING_EQUALITY(get_ele(*(x), 0), 0.841021, 1.0e-4);
397 TEST_FLOATING_EQUALITY(get_ele(*(x), 1), 0.541002, 1.0e-4);
398 }
399 }
400}
401
402// ************************************************************
403// ************************************************************
404TEUCHOS_UNIT_TEST(DIRK, SinCos)
405{
406 std::vector<std::string> RKMethods;
407 RKMethods.push_back("General DIRK");
408 RKMethods.push_back("RK Backward Euler");
409 RKMethods.push_back("DIRK 1 Stage Theta Method");
410 RKMethods.push_back("RK Implicit 1 Stage 1st order Radau IA");
411 RKMethods.push_back("RK Implicit Midpoint");
412 RKMethods.push_back("SDIRK 2 Stage 2nd order");
413 RKMethods.push_back("RK Implicit 2 Stage 2nd order Lobatto IIIB");
414 RKMethods.push_back("SDIRK 2 Stage 3rd order");
415 RKMethods.push_back("EDIRK 2 Stage 3rd order");
416 RKMethods.push_back("EDIRK 2 Stage Theta Method");
417 RKMethods.push_back("SDIRK 3 Stage 4th order");
418 RKMethods.push_back("SDIRK 5 Stage 4th order");
419 RKMethods.push_back("SDIRK 5 Stage 5th order");
420 RKMethods.push_back("SDIRK 2(1) Pair");
421 RKMethods.push_back("RK Trapezoidal Rule");
422 RKMethods.push_back("RK Crank-Nicolson");
423 RKMethods.push_back("SSPDIRK22");
424 RKMethods.push_back("SSPDIRK32");
425 RKMethods.push_back("SSPDIRK23");
426 RKMethods.push_back("SSPDIRK33");
427 RKMethods.push_back("SDIRK 3 Stage 2nd order");
428
429 std::vector<double> RKMethodErrors;
430 RKMethodErrors.push_back(2.52738e-05);
431 RKMethodErrors.push_back(0.0124201);
432 RKMethodErrors.push_back(5.20785e-05);
433 RKMethodErrors.push_back(0.0124201);
434 RKMethodErrors.push_back(5.20785e-05);
435 RKMethodErrors.push_back(2.52738e-05);
436 RKMethodErrors.push_back(5.20785e-05);
437 RKMethodErrors.push_back(1.40223e-06);
438 RKMethodErrors.push_back(2.17004e-07);
439 RKMethodErrors.push_back(5.20785e-05);
440 RKMethodErrors.push_back(6.41463e-08);
441 RKMethodErrors.push_back(3.30631e-10);
442 RKMethodErrors.push_back(1.35728e-11);
443 RKMethodErrors.push_back(0.0001041);
444 RKMethodErrors.push_back(5.20785e-05);
445 RKMethodErrors.push_back(5.20785e-05);
446 RKMethodErrors.push_back(1.30205e-05);
447 RKMethodErrors.push_back(5.7869767e-06);
448 RKMethodErrors.push_back(1.00713e-07);
449 RKMethodErrors.push_back(3.94916e-08);
450 RKMethodErrors.push_back(2.52738e-05);
451
452 TEUCHOS_ASSERT(RKMethods.size() == RKMethodErrors.size());
453
454 for (std::vector<std::string>::size_type m = 0; m != RKMethods.size(); m++) {
455 std::string RKMethod = RKMethods[m];
456 std::replace(RKMethod.begin(), RKMethod.end(), ' ', '_');
457 std::replace(RKMethod.begin(), RKMethod.end(), '/', '.');
458
459 RCP<Tempus::IntegratorBasic<double>> integrator;
460 std::vector<RCP<Thyra::VectorBase<double>>> solutions;
461 std::vector<RCP<Thyra::VectorBase<double>>> solutionsDot;
462 std::vector<double> StepSize;
463 std::vector<double> xErrorNorm;
464 std::vector<double> xDotErrorNorm;
465
466 const int nTimeStepSizes = 2; // 7 for error plots
467 double dt = 0.05;
468 double time = 0.0;
469 for (int n = 0; n < nTimeStepSizes; n++) {
470 // Read params from .xml file
471 RCP<ParameterList> pList =
472 getParametersFromXmlFile("Tempus_DIRK_SinCos.xml");
473
474 // Setup the SinCosModel
475 RCP<ParameterList> scm_pl = sublist(pList, "SinCosModel", true);
476 auto model = rcp(new SinCosModel<double>(scm_pl));
477
478 // Set the Stepper
479 RCP<ParameterList> pl = sublist(pList, "Tempus", true);
480 pl->sublist("Default Stepper").set("Stepper Type", RKMethods[m]);
481 if (RKMethods[m] == "DIRK 1 Stage Theta Method" ||
482 RKMethods[m] == "EDIRK 2 Stage Theta Method") {
483 pl->sublist("Default Stepper").set<double>("theta", 0.5);
484 }
485 else if (RKMethods[m] == "SDIRK 2 Stage 2nd order") {
486 pl->sublist("Default Stepper").set("gamma", 0.2928932188134524);
487 }
488 else if (RKMethods[m] == "SDIRK 2 Stage 3rd order") {
489 pl->sublist("Default Stepper")
490 .set<std::string>("Gamma Type", "3rd Order A-stable");
491 }
492
493 dt /= 2;
494
495 // Setup the Integrator and reset initial time step
496 pl->sublist("Default Integrator")
497 .sublist("Time Step Control")
498 .set("Initial Time Step", dt);
499 integrator = Tempus::createIntegratorBasic<double>(pl, model);
500
501 // Initial Conditions
502 // During the Integrator construction, the initial SolutionState
503 // is set by default to model->getNominalVales().get_x(). However,
504 // the application can set it also by
505 // integrator->initializeSolutionHistory.
506 RCP<Thyra::VectorBase<double>> x0 =
507 model->getNominalValues().get_x()->clone_v();
508 integrator->initializeSolutionHistory(0.0, x0);
509
510 // Integrate to timeMax
511 bool integratorStatus = integrator->advanceTime();
512 TEST_ASSERT(integratorStatus)
513
514 // Test if at 'Final Time'
515 time = integrator->getTime();
516 double timeFinal = pl->sublist("Default Integrator")
517 .sublist("Time Step Control")
518 .get<double>("Final Time");
519 double tol = 100.0 * std::numeric_limits<double>::epsilon();
520 TEST_FLOATING_EQUALITY(time, timeFinal, tol);
521
522 // Plot sample solution and exact solution
523 if (n == 0) {
524 RCP<const SolutionHistory<double>> solutionHistory =
525 integrator->getSolutionHistory();
526 writeSolution("Tempus_" + RKMethod + "_SinCos.dat", solutionHistory);
527
528 auto solnHistExact = rcp(new Tempus::SolutionHistory<double>());
529 for (int i = 0; i < solutionHistory->getNumStates(); i++) {
530 double time_i = (*solutionHistory)[i]->getTime();
531 auto state = Tempus::createSolutionStateX(
532 rcp_const_cast<Thyra::VectorBase<double>>(
533 model->getExactSolution(time_i).get_x()),
534 rcp_const_cast<Thyra::VectorBase<double>>(
535 model->getExactSolution(time_i).get_x_dot()));
536 state->setTime((*solutionHistory)[i]->getTime());
537 solnHistExact->addState(state);
538 }
539 writeSolution("Tempus_" + RKMethod + "_SinCos-Ref.dat", solnHistExact);
540 }
541
542 // Store off the final solution and step size
543 StepSize.push_back(dt);
544 auto solution = Thyra::createMember(model->get_x_space());
545 Thyra::copy(*(integrator->getX()), solution.ptr());
546 solutions.push_back(solution);
547 auto solutionDot = Thyra::createMember(model->get_x_space());
548 Thyra::copy(*(integrator->getXDot()), solutionDot.ptr());
549 solutionsDot.push_back(solutionDot);
550 if (n == nTimeStepSizes - 1) { // Add exact solution last in vector.
551 StepSize.push_back(0.0);
552 auto solutionExact = Thyra::createMember(model->get_x_space());
553 Thyra::copy(*(model->getExactSolution(time).get_x()),
554 solutionExact.ptr());
555 solutions.push_back(solutionExact);
556 auto solutionDotExact = Thyra::createMember(model->get_x_space());
557 Thyra::copy(*(model->getExactSolution(time).get_x_dot()),
558 solutionDotExact.ptr());
559 solutionsDot.push_back(solutionDotExact);
560 }
561 }
562
563 // Check the order and intercept
564 double xSlope = 0.0;
565 double xDotSlope = 0.0;
566 RCP<Tempus::Stepper<double>> stepper = integrator->getStepper();
567 double order = stepper->getOrder();
568 writeOrderError("Tempus_" + RKMethod + "_SinCos-Error.dat", stepper,
569 StepSize, solutions, xErrorNorm, xSlope, solutionsDot,
570 xDotErrorNorm, xDotSlope, out);
571
572 TEST_FLOATING_EQUALITY(xSlope, order, 0.01);
573 TEST_FLOATING_EQUALITY(xErrorNorm[0], RKMethodErrors[m], 5.0e-4);
574 // xDot not yet available for DIRK methods.
575 // TEST_FLOATING_EQUALITY( xDotSlope, order, 0.01 );
576 // TEST_FLOATING_EQUALITY( xDotErrorNorm[0], 0.0486418, 1.0e-4 );
577 }
578}
579
580// ************************************************************
581// ************************************************************
582TEUCHOS_UNIT_TEST(DIRK, VanDerPol)
583{
584 std::vector<std::string> RKMethods;
585 RKMethods.push_back("SDIRK 2 Stage 2nd order");
586
587 std::string RKMethod = RKMethods[0];
588 std::replace(RKMethod.begin(), RKMethod.end(), ' ', '_');
589 std::replace(RKMethod.begin(), RKMethod.end(), '/', '.');
590
591 RCP<Tempus::IntegratorBasic<double>> integrator;
592 std::vector<RCP<Thyra::VectorBase<double>>> solutions;
593 std::vector<RCP<Thyra::VectorBase<double>>> solutionsDot;
594 std::vector<double> StepSize;
595 std::vector<double> xErrorNorm;
596 std::vector<double> xDotErrorNorm;
597
598 const int nTimeStepSizes = 3; // 8 for error plot
599 double dt = 0.05;
600 double time = 0.0;
601 for (int n = 0; n < nTimeStepSizes; n++) {
602 // Read params from .xml file
603 RCP<ParameterList> pList =
604 getParametersFromXmlFile("Tempus_DIRK_VanDerPol.xml");
605
606 // Setup the VanDerPolModel
607 RCP<ParameterList> vdpm_pl = sublist(pList, "VanDerPolModel", true);
608 auto model = rcp(new VanDerPolModel<double>(vdpm_pl));
609
610 // Set the Stepper
611 RCP<ParameterList> pl = sublist(pList, "Tempus", true);
612 pl->sublist("Default Stepper").set("Stepper Type", RKMethods[0]);
613 pl->sublist("Default Stepper").set("gamma", 0.2928932188134524);
614
615 // Set the step size
616 dt /= 2;
617 if (n == nTimeStepSizes - 1) dt /= 10.0;
618
619 // Setup the Integrator and reset initial time step
620 pl->sublist("Default Integrator")
621 .sublist("Time Step Control")
622 .set("Initial Time Step", dt);
623 integrator = Tempus::createIntegratorBasic<double>(pl, model);
624
625 // Integrate to timeMax
626 bool integratorStatus = integrator->advanceTime();
627 TEST_ASSERT(integratorStatus)
628
629 // Test if at 'Final Time'
630 time = integrator->getTime();
631 double timeFinal = pl->sublist("Default Integrator")
632 .sublist("Time Step Control")
633 .get<double>("Final Time");
634 double tol = 100.0 * std::numeric_limits<double>::epsilon();
635 TEST_FLOATING_EQUALITY(time, timeFinal, tol);
636
637 // Store off the final solution and step size
638 StepSize.push_back(dt);
639 auto solution = Thyra::createMember(model->get_x_space());
640 Thyra::copy(*(integrator->getX()), solution.ptr());
641 solutions.push_back(solution);
642 auto solutionDot = Thyra::createMember(model->get_x_space());
643 Thyra::copy(*(integrator->getXDot()), solutionDot.ptr());
644 solutionsDot.push_back(solutionDot);
645
646 // Output finest temporal solution for plotting
647 // This only works for ONE MPI process
648 if ((n == 0) || (n == nTimeStepSizes - 1)) {
649 std::string fname = "Tempus_" + RKMethod + "_VanDerPol-Ref.dat";
650 if (n == 0) fname = "Tempus_" + RKMethod + "_VanDerPol.dat";
651 RCP<const SolutionHistory<double>> solutionHistory =
652 integrator->getSolutionHistory();
653 writeSolution(fname, solutionHistory);
654 }
655 }
656
657 // Check the order and intercept
658 double xSlope = 0.0;
659 double xDotSlope = 0.0;
660 RCP<Tempus::Stepper<double>> stepper = integrator->getStepper();
661 double order = stepper->getOrder();
662
663 // xDot not yet available for DIRK methods, e.g., are not calc. and zero.
664 solutionsDot.clear();
665
666 writeOrderError("Tempus_" + RKMethod + "_VanDerPol-Error.dat", stepper,
667 StepSize, solutions, xErrorNorm, xSlope, solutionsDot,
668 xDotErrorNorm, xDotSlope, out);
669
670 TEST_FLOATING_EQUALITY(xSlope, order, 0.06);
671 TEST_FLOATING_EQUALITY(xErrorNorm[0], 1.07525e-05, 1.0e-4);
672 // TEST_FLOATING_EQUALITY( xDotSlope, 1.74898, 0.10 );
673 // TEST_FLOATING_EQUALITY( xDotErrorNorm[0], 1.0038, 1.0e-4 );
674
675 Teuchos::TimeMonitor::summarize();
676}
677
678// ************************************************************
679// ************************************************************
680TEUCHOS_UNIT_TEST(DIRK, EmbeddedVanDerPol)
681{
682 std::vector<std::string> IntegratorList;
683 IntegratorList.push_back("Embedded_Integrator_PID");
684 IntegratorList.push_back("Embedded_Integrator");
685
686 // the embedded solution will test the following:
687 const int refIstep = 217;
688
689 for (auto integratorChoice : IntegratorList) {
690 out << "Using Integrator: " << integratorChoice << " !!!" << std::endl;
691
692 // Read params from .xml file
693 RCP<ParameterList> pList =
694 getParametersFromXmlFile("Tempus_DIRK_VanDerPol.xml");
695
696 // Setup the VanDerPolModel
697 RCP<ParameterList> vdpm_pl = sublist(pList, "VanDerPolModel", true);
698 auto model = rcp(new VanDerPolModel<double>(vdpm_pl));
699
700 // Set the Integrator and Stepper
701 RCP<ParameterList> pl = sublist(pList, "Tempus", true);
702 pl->set("Integrator Name", integratorChoice);
703
704 // Setup the Integrator
705 RCP<Tempus::IntegratorBasic<double>> integrator =
706 Tempus::createIntegratorBasic<double>(pl, model);
707
708 const std::string RKMethod_ =
709 pl->sublist(integratorChoice).get<std::string>("Stepper Name");
710
711 // Integrate to timeMax
712 bool integratorStatus = integrator->advanceTime();
713 TEST_ASSERT(integratorStatus);
714
715 // Test if at 'Final Time'
716 double time = integrator->getTime();
717 double timeFinal = pl->sublist(integratorChoice)
718 .sublist("Time Step Control")
719 .get<double>("Final Time");
720 TEST_FLOATING_EQUALITY(time, timeFinal, 1.0e-14);
721
722 // Numerical reference solution at timeFinal (for \epsilon = 0.1)
723 RCP<Thyra::VectorBase<double>> x = integrator->getX();
724 RCP<Thyra::VectorBase<double>> xref = x->clone_v();
725 Thyra::set_ele(0, -1.5484458614405929, xref.ptr());
726 Thyra::set_ele(1, 1.0181127316101317, xref.ptr());
727
728 // Calculate the error
729 RCP<Thyra::VectorBase<double>> xdiff = x->clone_v();
730 Thyra::V_StVpStV(xdiff.ptr(), 1.0, *xref, -1.0, *(x));
731 const double L2norm = Thyra::norm_2(*xdiff);
732
733 // Test number of steps, failures, and accuracy
734 if (integratorChoice == "Embedded_Integrator_PID") {
735 const double absTol = pl->sublist(integratorChoice)
736 .sublist("Time Step Control")
737 .get<double>("Maximum Absolute Error");
738 const double relTol = pl->sublist(integratorChoice)
739 .sublist("Time Step Control")
740 .get<double>("Maximum Relative Error");
741
742 // get the number of time steps and number of step failure
743 // const int nFailure_c = integrator->getSolutionHistory()->
744 // getCurrentState()->getMetaData()->getNFailures();
745 const int iStep =
746 integrator->getSolutionHistory()->getCurrentState()->getIndex();
747 // const int nFail = integrator->getSolutionHistory()->
748 // getCurrentState()->getMetaData()->getNRunningFailures();
749
750 // Should be close to the prescribed tolerance
751 TEST_FLOATING_EQUALITY(std::log10(L2norm), std::log10(absTol), 0.3);
752 TEST_FLOATING_EQUALITY(std::log10(L2norm), std::log10(relTol), 0.3);
753 // test for number of steps
754 TEST_COMPARE(iStep, <=, refIstep);
755 }
756
757 // Plot sample solution and exact solution
758 std::ofstream ftmp("Tempus_" + integratorChoice + RKMethod_ +
759 "_VDP_Example.dat");
760 RCP<const SolutionHistory<double>> solutionHistory =
761 integrator->getSolutionHistory();
762 int nStates = solutionHistory->getNumStates();
763 // RCP<const Thyra::VectorBase<double> > x_exact_plot;
764 for (int i = 0; i < nStates; i++) {
765 RCP<const SolutionState<double>> solutionState = (*solutionHistory)[i];
766 double time_i = solutionState->getTime();
767 RCP<const Thyra::VectorBase<double>> x_plot = solutionState->getX();
768 // x_exact_plot = model->getExactSolution(time_i).get_x();
769 ftmp << time_i << " " << Thyra::get_ele(*(x_plot), 0) << " "
770 << Thyra::get_ele(*(x_plot), 1) << " " << std::endl;
771 }
772 ftmp.close();
773 }
774
775 Teuchos::TimeMonitor::summarize();
776}
777
778} // 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