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
36using Tempus::IntegratorBasic;
37using Tempus::SolutionHistory;
38using Tempus::SolutionState;
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
Tempus_ConvergenceTestUtils.hpp
Tempus::SolutionHistory
SolutionHistory is basically a container of SolutionStates. SolutionHistory maintains a collection of...
Definition Tempus_SolutionHistory_decl.hpp:121
Tempus::TimeStepControl
TimeStepControl manages the time step size. There several mechanisms that effect the time step size a...
Definition Tempus_TimeStepControl_decl.hpp:53
Tempus_Test::SinCosModel
Sine-Cosine model problem from Rythmos. This is a canonical Sine-Cosine differential equation.
Definition SinCosModel_decl.hpp:94
Tempus_Test::VanDerPolModel
van der Pol model problem for nonlinear electrical circuit.
Definition VanDerPolModel_decl.hpp:114
Tempus_Test
Definition Tempus_BackwardEuler_ASA.cpp:34
Tempus_Test::writeOrderError
void writeOrderError(const std::string filename, Teuchos::RCP< Tempus::Stepper< Scalar > > stepper, std::vector< Scalar > &StepSize, std::vector< Teuchos::RCP< Thyra::VectorBase< Scalar > > > &solutions, std::vector< Scalar > &xErrorNorm, Scalar &xSlope, std::vector< Teuchos::RCP< Thyra::VectorBase< Scalar > > > &solutionsDot, std::vector< Scalar > &xDotErrorNorm, Scalar &xDotSlope, std::vector< Teuchos::RCP< Thyra::VectorBase< Scalar > > > &solutionsDotDot, std::vector< Scalar > &xDotDotErrorNorm, Scalar &xDotDotSlope)
Definition Tempus_ConvergenceTestUtils.hpp:189
Tempus_Test::writeSolution
void writeSolution(const std::string filename, Teuchos::RCP< const Tempus::SolutionHistory< Scalar > > solutionHistory)
Definition Tempus_ConvergenceTestUtils.hpp:311
Tempus_Test::TEUCHOS_UNIT_TEST
TEUCHOS_UNIT_TEST(BackwardEuler, SinCos_ASA)
Definition Tempus_BackwardEuler_ASA.cpp:47
Tempus::STORAGE_TYPE_STATIC
@ STORAGE_TYPE_STATIC
Keep a fix number of states.
Definition Tempus_SolutionHistory_decl.hpp:23
Tempus::createSolutionStateX
Teuchos::RCP< SolutionState< Scalar > > createSolutionStateX(const Teuchos::RCP< Thyra::VectorBase< Scalar > > &x, const Teuchos::RCP< Thyra::VectorBase< Scalar > > &xdot=Teuchos::null, const Teuchos::RCP< Thyra::VectorBase< Scalar > > &xdotdot=Teuchos::null)
Nonmember constructor from non-const solution vectors, x.
Definition Tempus_SolutionState_impl.hpp:355
Tempus::createIntegratorBasic
Teuchos::RCP< IntegratorBasic< Scalar > > createIntegratorBasic()
Nonmember constructor.
Definition Tempus_IntegratorBasic_impl.hpp:730
Generated by doxygen 1.13.2