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