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