Tempus Version of the Day
Time Integration
Loading...
Searching...
No Matches
Tempus_BDF2Test.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 "Tempus_config.hpp"
15#include "Tempus_IntegratorBasic.hpp"
16#include "Tempus_StepperBDF2.hpp"
17
18#include "../TestModels/SinCosModel.hpp"
19#include "../TestModels/CDR_Model.hpp"
20#include "../TestModels/VanDerPolModel.hpp"
22
23#include "Stratimikos_DefaultLinearSolverBuilder.hpp"
24#include "Thyra_LinearOpWithSolveFactoryHelpers.hpp"
25
26#ifdef Tempus_ENABLE_MPI
27#include "Epetra_MpiComm.h"
28#else
29#include "Epetra_SerialComm.h"
30#endif
31
32#include <fstream>
33#include <limits>
34#include <sstream>
35#include <vector>
36
37namespace Tempus_Test {
38
39using Teuchos::RCP;
40using Teuchos::rcp;
41using Teuchos::rcp_const_cast;
42using Teuchos::ParameterList;
43using Teuchos::sublist;
44using Teuchos::getParametersFromXmlFile;
45
49
50
51// ************************************************************
52// ************************************************************
53TEUCHOS_UNIT_TEST(BDF2, ParameterList)
54{
55 // Read params from .xml file
56 RCP<ParameterList> pList =
57 getParametersFromXmlFile("Tempus_BDF2_SinCos.xml");
58
59 // Setup the SinCosModel
60 RCP<ParameterList> scm_pl = sublist(pList, "SinCosModel", true);
61 auto model = rcp(new SinCosModel<double> (scm_pl));
62
63 RCP<ParameterList> tempusPL = sublist(pList, "Tempus", true);
64
65 // Test constructor IntegratorBasic(tempusPL, model)
66 {
67 RCP<Tempus::IntegratorBasic<double> > integrator =
68 Tempus::createIntegratorBasic<double>(tempusPL, model);
69
70 RCP<ParameterList> stepperPL = sublist(tempusPL, "Default Stepper", true);
71 RCP<const ParameterList> defaultPL =
72 integrator->getStepper()->getValidParameters();
73 bool pass = haveSameValuesSorted(*stepperPL, *defaultPL, true);
74 if (!pass) {
75 out << std::endl;
76 out << "stepperPL -------------- \n" << *stepperPL << std::endl;
77 out << "defaultPL -------------- \n" << *defaultPL << std::endl;
78 }
79 TEST_ASSERT(pass)
80 }
81
82 // Test constructor IntegratorBasic(model, stepperType)
83 {
84 RCP<Tempus::IntegratorBasic<double> > integrator =
85 Tempus::createIntegratorBasic<double>(model, std::string("BDF2"));
86
87 RCP<ParameterList> stepperPL = sublist(tempusPL, "Default Stepper", true);
88 RCP<const ParameterList> defaultPL =
89 integrator->getStepper()->getValidParameters();
90
91 bool pass = haveSameValuesSorted(*stepperPL, *defaultPL, true);
92 if (!pass) {
93 out << std::endl;
94 out << "stepperPL -------------- \n" << *stepperPL << std::endl;
95 out << "defaultPL -------------- \n" << *defaultPL << std::endl;
96 }
97 TEST_ASSERT(pass)
98 }
99}
100
101
102// ************************************************************
103// ************************************************************
104TEUCHOS_UNIT_TEST(BDF2, ConstructingFromDefaults)
105{
106 double dt = 0.1;
107 std::vector<std::string> options;
108 options.push_back("Default Parameters");
109 options.push_back("ICConsistency and Check");
110
111 for(const auto& option: options) {
112
113 // Read params from .xml file
114 RCP<ParameterList> pList =
115 getParametersFromXmlFile("Tempus_BDF2_SinCos.xml");
116 RCP<ParameterList> pl = sublist(pList, "Tempus", true);
117
118 // Setup the SinCosModel
119 RCP<ParameterList> scm_pl = sublist(pList, "SinCosModel", true);
120 //RCP<SinCosModel<double> > model = sineCosineModel(scm_pl);
121 auto model = rcp(new SinCosModel<double>(scm_pl));
122
123 // Setup Stepper for field solve ----------------------------
124 auto stepper = rcp(new Tempus::StepperBDF2<double>());
125 stepper->setModel(model);
126 if ( option == "ICConsistency and Check") {
127 stepper->setICConsistency("Consistent");
128 stepper->setICConsistencyCheck(true);
129 }
130 stepper->initialize();
131
132 // Setup TimeStepControl ------------------------------------
133 auto timeStepControl = rcp(new Tempus::TimeStepControl<double>());
134 ParameterList tscPL = pl->sublist("Default Integrator")
135 .sublist("Time Step Control");
136 timeStepControl->setInitIndex(tscPL.get<int> ("Initial Time Index"));
137 timeStepControl->setInitTime (tscPL.get<double>("Initial Time"));
138 timeStepControl->setFinalTime(tscPL.get<double>("Final Time"));
139 timeStepControl->setInitTimeStep(dt);
140 timeStepControl->initialize();
141
142 // Setup initial condition SolutionState --------------------
143 auto inArgsIC = model->getNominalValues();
144 auto icSolution = rcp_const_cast<Thyra::VectorBase<double> > (inArgsIC.get_x());
145 auto icState = Tempus::createSolutionStateX(icSolution);
146 icState->setTime (timeStepControl->getInitTime());
147 icState->setIndex (timeStepControl->getInitIndex());
148 icState->setTimeStep(0.0);
149 icState->setOrder (stepper->getOrder());
150 icState->setSolutionStatus(Tempus::Status::PASSED); // ICs are passing.
151
152 // Setup SolutionHistory ------------------------------------
153 auto solutionHistory = rcp(new Tempus::SolutionHistory<double>());
154 solutionHistory->setName("Forward States");
155 solutionHistory->setStorageType(Tempus::STORAGE_TYPE_STATIC);
156 solutionHistory->setStorageLimit(3);
157 solutionHistory->addState(icState);
158
159 // Ensure ICs are consistent and stepper memory is set (e.g., xDot is set).
160 stepper->setInitialConditions(solutionHistory);
161
162 // Setup Integrator -----------------------------------------
163 RCP<Tempus::IntegratorBasic<double> > integrator =
164 Tempus::createIntegratorBasic<double>();
165 integrator->setStepper(stepper);
166 integrator->setTimeStepControl(timeStepControl);
167 integrator->setSolutionHistory(solutionHistory);
168 //integrator->setObserver(...);
169 integrator->initialize();
170
171
172 // Integrate to timeMax
173 bool integratorStatus = integrator->advanceTime();
174 TEST_ASSERT(integratorStatus)
175
176
177 // Test if at 'Final Time'
178 double time = integrator->getTime();
179 double timeFinal =pl->sublist("Default Integrator")
180 .sublist("Time Step Control").get<double>("Final Time");
181 TEST_FLOATING_EQUALITY(time, timeFinal, 1.0e-14);
182
183 // Time-integrated solution and the exact solution
184 RCP<Thyra::VectorBase<double> > x = integrator->getX();
185 RCP<const Thyra::VectorBase<double> > x_exact =
186 model->getExactSolution(time).get_x();
187
188 // Calculate the error
189 RCP<Thyra::VectorBase<double> > xdiff = x->clone_v();
190 Thyra::V_StVpStV(xdiff.ptr(), 1.0, *x_exact, -1.0, *(x));
191
192 // Check the order and intercept
193 out << " Stepper = " << stepper->description()
194 << "\n with " << option << std::endl;
195 out << " =========================" << std::endl;
196 out << " Exact solution : " << get_ele(*(x_exact), 0) << " "
197 << get_ele(*(x_exact), 1) << std::endl;
198 out << " Computed solution: " << get_ele(*(x ), 0) << " "
199 << get_ele(*(x ), 1) << std::endl;
200 out << " Difference : " << get_ele(*(xdiff ), 0) << " "
201 << get_ele(*(xdiff ), 1) << std::endl;
202 out << " =========================" << std::endl;
203 TEST_FLOATING_EQUALITY(get_ele(*(x), 0), 0.839732, 1.0e-4 );
204 TEST_FLOATING_EQUALITY(get_ele(*(x), 1), 0.542663, 1.0e-4 );
205 }
206}
207
208
209// ************************************************************
210// ************************************************************
211TEUCHOS_UNIT_TEST(BDF2, SinCos)
212{
213 RCP<Tempus::IntegratorBasic<double> > integrator;
214 std::vector<RCP<Thyra::VectorBase<double>>> solutions;
215 std::vector<RCP<Thyra::VectorBase<double>>> solutionsDot;
216 std::vector<double> StepSize;
217
218 // Read params from .xml file
219 RCP<ParameterList> pList = getParametersFromXmlFile("Tempus_BDF2_SinCos.xml");
220 //Set initial time step = 2*dt specified in input file (for convergence study)
221 double dt = pList->sublist("Tempus")
222 .sublist("Default Integrator")
223 .sublist("Time Step Control").get<double>("Initial Time Step");
224 dt *= 2.0;
225
226 // Setup the SinCosModel
227 RCP<ParameterList> scm_pl = sublist(pList, "SinCosModel", true);
228 const int nTimeStepSizes = scm_pl->get<int>("Number of Time Step Sizes", 7);
229 std::string output_file_string =
230 scm_pl->get<std::string>("Output File Name", "Tempus_BDF2_SinCos");
231 std::string output_file_name = output_file_string + ".dat";
232 std::string ref_out_file_name = output_file_string + "-Ref.dat";
233 std::string err_out_file_name = output_file_string + "-Error.dat";
234 double time = 0.0;
235 for (int n=0; n<nTimeStepSizes; n++) {
236
237 auto model = rcp(new SinCosModel<double>(scm_pl));
238
239 dt /= 2;
240
241 // Setup the Integrator and reset initial time step
242 RCP<ParameterList> tempusPL =
243 getParametersFromXmlFile("Tempus_BDF2_SinCos.xml");
244 RCP<ParameterList> pl = sublist(tempusPL, "Tempus", true);
245 pl->sublist("Default Integrator")
246 .sublist("Time Step Control").set("Initial Time Step", dt);
247 integrator = Tempus::createIntegratorBasic<double>(pl, model);
248
249 // Initial Conditions
250 // During the Integrator construction, the initial SolutionState
251 // is set by default to model->getNominalVales().get_x(). However,
252 // the application can set it also by integrator->initializeSolutionHistory.
253 RCP<Thyra::VectorBase<double> > x0 =
254 model->getNominalValues().get_x()->clone_v();
255 integrator->initializeSolutionHistory(0.0, x0);
256
257 // Integrate to timeMax
258 bool integratorStatus = integrator->advanceTime();
259 TEST_ASSERT(integratorStatus)
260
261 // Test if at 'Final Time'
262 time = integrator->getTime();
263 double timeFinal = pl->sublist("Default Integrator")
264 .sublist("Time Step Control").get<double>("Final Time");
265 TEST_FLOATING_EQUALITY(time, timeFinal, 1.0e-14);
266
267 // Plot sample solution and exact solution
268 if (n == 0) {
269 RCP<const SolutionHistory<double> > solutionHistory =
270 integrator->getSolutionHistory();
271 writeSolution(output_file_name, solutionHistory);
272
273 auto solnHistExact = rcp(new Tempus::SolutionHistory<double>());
274 for (int i=0; i<solutionHistory->getNumStates(); i++) {
275 double time_i = (*solutionHistory)[i]->getTime();
276 auto state = Tempus::createSolutionStateX(
277 rcp_const_cast<Thyra::VectorBase<double> > (
278 model->getExactSolution(time_i).get_x()),
279 rcp_const_cast<Thyra::VectorBase<double> > (
280 model->getExactSolution(time_i).get_x_dot()));
281 state->setTime((*solutionHistory)[i]->getTime());
282 solnHistExact->addState(state);
283 }
284 writeSolution(ref_out_file_name, solnHistExact);
285 }
286
287 // Store off the final solution and step size
288 StepSize.push_back(dt);
289 auto solution = Thyra::createMember(model->get_x_space());
290 Thyra::copy(*(integrator->getX()),solution.ptr());
291 solutions.push_back(solution);
292 auto solutionDot = Thyra::createMember(model->get_x_space());
293 Thyra::copy(*(integrator->getXDot()),solutionDot.ptr());
294 solutionsDot.push_back(solutionDot);
295 if (n == nTimeStepSizes-1) { // Add exact solution last in vector.
296 StepSize.push_back(0.0);
297 auto solutionExact = Thyra::createMember(model->get_x_space());
298 Thyra::copy(*(model->getExactSolution(time).get_x()),solutionExact.ptr());
299 solutions.push_back(solutionExact);
300 auto solutionDotExact = Thyra::createMember(model->get_x_space());
301 Thyra::copy(*(model->getExactSolution(time).get_x_dot()),
302 solutionDotExact.ptr());
303 solutionsDot.push_back(solutionDotExact);
304 }
305 }
306
307 // Check the order and intercept
308 if (nTimeStepSizes > 1) {
309 double xSlope = 0.0;
310 double xDotSlope = 0.0;
311 std::vector<double> xErrorNorm;
312 std::vector<double> xDotErrorNorm;
313 RCP<Tempus::Stepper<double> > stepper = integrator->getStepper();
314 double order = stepper->getOrder();
315 writeOrderError(err_out_file_name,
316 stepper, StepSize,
317 solutions, xErrorNorm, xSlope,
318 solutionsDot, xDotErrorNorm, xDotSlope);
319
320 TEST_FLOATING_EQUALITY( xSlope, order, 0.01 );
321 TEST_FLOATING_EQUALITY( xDotSlope, order, 0.01 );
322 TEST_FLOATING_EQUALITY( xErrorNorm[0], 5.13425e-05, 1.0e-4 );
323 TEST_FLOATING_EQUALITY( xDotErrorNorm[0], 5.13425e-05, 1.0e-4 );
324 }
325
326 Teuchos::TimeMonitor::summarize();
327}
328
329
330// ************************************************************
331// ************************************************************
332TEUCHOS_UNIT_TEST(BDF2, SinCosAdapt)
333{
334 RCP<Tempus::IntegratorBasic<double> > integrator;
335 std::vector<RCP<Thyra::VectorBase<double>>> solutions;
336 std::vector<RCP<Thyra::VectorBase<double>>> solutionsDot;
337 std::vector<double> StepSize;
338
339 // Read params from .xml file
340 RCP<ParameterList> pList =
341 getParametersFromXmlFile("Tempus_BDF2_SinCos_AdaptDt.xml");
342 //Set initial time step = 2*dt specified in input file (for convergence study)
343 double dt = pList->sublist("Tempus")
344 .sublist("Default Integrator")
345 .sublist("Time Step Control").get<double>("Initial Time Step");
346 dt *= 2.0;
347
348 // Setup the SinCosModel
349 RCP<ParameterList> scm_pl = sublist(pList, "SinCosModel", true);
350 const int nTimeStepSizes = scm_pl->get<int>("Number of Time Step Sizes", 7);
351 std::string output_file_string =
352 scm_pl->get<std::string>("Output File Name", "Tempus_BDF2_SinCos");
353 std::string output_file_name = output_file_string + ".dat";
354 std::string err_out_file_name = output_file_string + "-Error.dat";
355 double time = 0.0;
356 for (int n=0; n<nTimeStepSizes; n++) {
357
358 auto model = rcp(new SinCosModel<double>(scm_pl));
359
360 dt /= 2;
361
362 RCP<ParameterList> tempusPL =
363 getParametersFromXmlFile("Tempus_BDF2_SinCos_AdaptDt.xml");
364 RCP<ParameterList> pl = sublist(tempusPL, "Tempus", true);
365
366 // Setup the Integrator and reset initial time step
367 pl->sublist("Default Integrator")
368 .sublist("Time Step Control").set("Initial Time Step", dt/4.0);
369 // Ensure time step does not get larger than the initial time step size,
370 // as that would mess up the convergence rates.
371 pl->sublist("Default Integrator")
372 .sublist("Time Step Control").set("Maximum Time Step", dt);
373 // Ensure time step does not get too small and therefore too many steps.
374 pl->sublist("Default Integrator")
375 .sublist("Time Step Control").set("Minimum Time Step", dt/4.0);
376 // For the SinCos problem eta is directly related to dt
377 pl->sublist("Default Integrator")
378 .sublist("Time Step Control")
379 .sublist("Time Step Control Strategy")
380 .set("Minimum Value Monitoring Function", dt*0.99);
381 integrator = Tempus::createIntegratorBasic<double>(pl, model);
382
383 // Initial Conditions
384 // During the Integrator construction, the initial SolutionState
385 // is set by default to model->getNominalVales().get_x(). However,
386 // the application can set it also by integrator->initializeSolutionHistory.
387 RCP<Thyra::VectorBase<double> > x0 =
388 model->getNominalValues().get_x()->clone_v();
389 integrator->initializeSolutionHistory(0.0, x0);
390
391 // Integrate to timeMax
392 bool integratorStatus = integrator->advanceTime();
393 TEST_ASSERT(integratorStatus)
394
395 // Test if at 'Final Time'
396 time = integrator->getTime();
397 double timeFinal = pl->sublist("Default Integrator")
398 .sublist("Time Step Control").get<double>("Final Time");
399 TEST_FLOATING_EQUALITY(time, timeFinal, 1.0e-14);
400
401 // Time-integrated solution and the exact solution
402 RCP<Thyra::VectorBase<double> > x = integrator->getX();
403 RCP<const Thyra::VectorBase<double> > x_exact =
404 model->getExactSolution(time).get_x();
405
406 // Plot sample solution and exact solution
407 if (n == 0) {
408 std::ofstream ftmp(output_file_name);
409 //Warning: the following assumes serial run
410 FILE *gold_file = fopen("Tempus_BDF2_SinCos_AdaptDt_gold.dat", "r");
411 RCP<const SolutionHistory<double> > solutionHistory =
412 integrator->getSolutionHistory();
413 RCP<const Thyra::VectorBase<double> > x_exact_plot;
414 for (int i=0; i<solutionHistory->getNumStates(); i++) {
415 char time_gold_char[100];
416 fgets(time_gold_char, 100, gold_file);
417 double time_gold;
418 sscanf(time_gold_char, "%lf", &time_gold);
419 RCP<const SolutionState<double> > solutionState = (*solutionHistory)[i];
420 double time_i = solutionState->getTime();
421 //Throw error if time does not match time in gold file to specified tolerance
422 TEST_FLOATING_EQUALITY( time_i, time_gold, 1.0e-5 );
423 RCP<const Thyra::VectorBase<double> > x_plot = solutionState->getX();
424 x_exact_plot = model->getExactSolution(time_i).get_x();
425 ftmp << time_i << " "
426 << get_ele(*(x_plot), 0) << " "
427 << get_ele(*(x_plot), 1) << " "
428 << get_ele(*(x_exact_plot), 0) << " "
429 << get_ele(*(x_exact_plot), 1) << std::endl;
430 }
431 ftmp.close();
432 }
433
434 // Store off the final solution and step size
435 StepSize.push_back(dt);
436 auto solution = Thyra::createMember(model->get_x_space());
437 Thyra::copy(*(integrator->getX()),solution.ptr());
438 solutions.push_back(solution);
439 auto solutionDot = Thyra::createMember(model->get_x_space());
440 Thyra::copy(*(integrator->getXDot()),solutionDot.ptr());
441 solutionsDot.push_back(solutionDot);
442 if (n == nTimeStepSizes-1) { // Add exact solution last in vector.
443 StepSize.push_back(0.0);
444 auto solutionExact = Thyra::createMember(model->get_x_space());
445 Thyra::copy(*(model->getExactSolution(time).get_x()),solutionExact.ptr());
446 solutions.push_back(solutionExact);
447 auto solutionDotExact = Thyra::createMember(model->get_x_space());
448 Thyra::copy(*(model->getExactSolution(time).get_x_dot()),
449 solutionDotExact.ptr());
450 solutionsDot.push_back(solutionDotExact);
451 }
452 }
453
454 // Check the order and intercept
455 if (nTimeStepSizes > 1) {
456 double xSlope = 0.0;
457 double xDotSlope = 0.0;
458 std::vector<double> xErrorNorm;
459 std::vector<double> xDotErrorNorm;
460 RCP<Tempus::Stepper<double> > stepper = integrator->getStepper();
461 //double order = stepper->getOrder();
462 writeOrderError("Tempus_BDF2_SinCos-Error.dat",
463 stepper, StepSize,
464 solutions, xErrorNorm, xSlope,
465 solutionsDot, xDotErrorNorm, xDotSlope);
466
467 TEST_FLOATING_EQUALITY( xSlope, 1.932, 0.01 );
468 TEST_FLOATING_EQUALITY( xDotSlope, 1.932, 0.01 );
469 TEST_FLOATING_EQUALITY( xErrorNorm[0], 0.000192591, 1.0e-4 );
470 TEST_FLOATING_EQUALITY( xDotErrorNorm[0], 0.000192591, 1.0e-4 );
471 }
472
473 Teuchos::TimeMonitor::summarize();
474}
475
476
477// ************************************************************
478// ************************************************************
480{
481 // Create a communicator for Epetra objects
482 RCP<Epetra_Comm> comm;
483#ifdef Tempus_ENABLE_MPI
484 comm = rcp(new Epetra_MpiComm(MPI_COMM_WORLD));
485#else
486 comm = rcp(new Epetra_SerialComm);
487#endif
488
489 RCP<Tempus::IntegratorBasic<double> > integrator;
490 std::vector<RCP<Thyra::VectorBase<double>>> solutions;
491 std::vector<RCP<Thyra::VectorBase<double>>> solutionsDot;
492 std::vector<double> StepSize;
493
494 // Read params from .xml file
495 RCP<ParameterList> pList =
496 getParametersFromXmlFile("Tempus_BDF2_CDR.xml");
497 //Set initial time step = 2*dt specified in input file (for convergence study)
498 //
499 RCP<ParameterList> pl = sublist(pList, "Tempus", true);
500 double dt = pl->sublist("Demo Integrator")
501 .sublist("Time Step Control").get<double>("Initial Time Step");
502 dt *= 2.0;
503 RCP<ParameterList> model_pl = sublist(pList, "CDR Model", true);
504
505 const int nTimeStepSizes = model_pl->get<int>("Number of Time Step Sizes", 5);
506
507 for (int n=0; n<nTimeStepSizes; n++) {
508
509 // Create CDR Model
510 const int num_elements = model_pl->get<int>("num elements");
511 const double left_end = model_pl->get<double>("left end");
512 const double right_end = model_pl->get<double>("right end");
513 const double a_convection = model_pl->get<double>("a (convection)");
514 const double k_source = model_pl->get<double>("k (source)");
515
516 auto model = rcp(new Tempus_Test::CDR_Model<double>(comm,
517 num_elements,
518 left_end,
519 right_end,
520 a_convection,
521 k_source));
522
523 // Set the factory
524 ::Stratimikos::DefaultLinearSolverBuilder builder;
525
526 auto p = rcp(new ParameterList);
527 p->set("Linear Solver Type", "Belos");
528 p->set("Preconditioner Type", "None");
529 builder.setParameterList(p);
530
531 RCP< ::Thyra::LinearOpWithSolveFactoryBase<double> >
532 lowsFactory = builder.createLinearSolveStrategy("");
533
534 model->set_W_factory(lowsFactory);
535
536 // Set the step size
537 dt /= 2;
538
539 // Setup the Integrator and reset initial time step
540 pl->sublist("Demo Integrator")
541 .sublist("Time Step Control").set("Initial Time Step", dt);
542 integrator = Tempus::createIntegratorBasic<double>(pl, model);
543
544 // Integrate to timeMax
545 bool integratorStatus = integrator->advanceTime();
546 TEST_ASSERT(integratorStatus)
547
548 // Test if at 'Final Time'
549 double time = integrator->getTime();
550 double timeFinal =pl->sublist("Demo Integrator")
551 .sublist("Time Step Control").get<double>("Final Time");
552 double tol = 100.0 * std::numeric_limits<double>::epsilon();
553 TEST_FLOATING_EQUALITY(time, timeFinal, tol);
554
555 // Store off the final solution and step size
556 StepSize.push_back(dt);
557 auto solution = Thyra::createMember(model->get_x_space());
558 Thyra::copy(*(integrator->getX()),solution.ptr());
559 solutions.push_back(solution);
560 auto solutionDot = Thyra::createMember(model->get_x_space());
561 Thyra::copy(*(integrator->getXDot()),solutionDot.ptr());
562 solutionsDot.push_back(solutionDot);
563
564 // Output finest temporal solution for plotting
565 // This only works for ONE MPI process
566 if ((n == nTimeStepSizes-1) && (comm->NumProc() == 1)) {
567 std::ofstream ftmp("Tempus_BDF2_CDR.dat");
568 ftmp << "TITLE=\"BDF2 Solution to CDR\"\n"
569 << "VARIABLES=\"z\",\"T\"\n";
570 const double dx = std::fabs(left_end-right_end) /
571 static_cast<double>(num_elements);
572 RCP<const SolutionHistory<double> > solutionHistory =
573 integrator->getSolutionHistory();
574 int nStates = solutionHistory->getNumStates();
575 for (int i=0; i<nStates; i++) {
576 RCP<const SolutionState<double> > solutionState = (*solutionHistory)[i];
577 RCP<const Thyra::VectorBase<double> > x = solutionState->getX();
578 double ttime = solutionState->getTime();
579 ftmp << "ZONE T=\"Time="<<ttime<<"\", I="
580 <<num_elements+1<<", F=BLOCK\n";
581 for (int j = 0; j < num_elements+1; j++) {
582 const double x_coord = left_end + static_cast<double>(j) * dx;
583 ftmp << x_coord << " ";
584 }
585 ftmp << std::endl;
586 for (int j=0; j<num_elements+1; j++) ftmp << get_ele(*x, j) << " ";
587 ftmp << std::endl;
588 }
589 ftmp.close();
590 }
591 }
592
593 // Check the order and intercept
594 if (nTimeStepSizes > 2) {
595 double xSlope = 0.0;
596 double xDotSlope = 0.0;
597 std::vector<double> xErrorNorm;
598 std::vector<double> xDotErrorNorm;
599 RCP<Tempus::Stepper<double> > stepper = integrator->getStepper();
600 double order = stepper->getOrder();
601 writeOrderError("Tempus_BDF2_CDR-Error.dat",
602 stepper, StepSize,
603 solutions, xErrorNorm, xSlope,
604 solutionsDot, xDotErrorNorm, xDotSlope);
605 TEST_FLOATING_EQUALITY( xSlope, order, 0.35 );
606 TEST_COMPARE(xSlope, >, 0.95);
607 TEST_FLOATING_EQUALITY( xDotSlope, order, 0.35 );
608 TEST_COMPARE(xDotSlope, >, 0.95);
609
610 TEST_FLOATING_EQUALITY( xErrorNorm[0], 0.0145747, 1.0e-4 );
611 TEST_FLOATING_EQUALITY( xDotErrorNorm[0], 0.0563621, 1.0e-4 );
612 }
613
614 // Write fine mesh solution at final time
615 // This only works for ONE MPI process
616 if (comm->NumProc() == 1) {
617 RCP<ParameterList> pListCDR =
618 getParametersFromXmlFile("Tempus_BDF2_CDR.xml");
619 RCP<ParameterList> model_pl_CDR = sublist(pListCDR, "CDR Model", true);
620 const int num_elements = model_pl_CDR->get<int>("num elements");
621 const double left_end = model_pl_CDR->get<double>("left end");
622 const double right_end = model_pl_CDR->get<double>("right end");
623
624 const Thyra::VectorBase<double>& x = *(solutions[solutions.size()-1]);
625
626 std::ofstream ftmp("Tempus_BDF2_CDR-Solution.dat");
627 for (int n = 0; n < num_elements+1; n++) {
628 const double dx = std::fabs(left_end-right_end) /
629 static_cast<double>(num_elements);
630 const double x_coord = left_end + static_cast<double>(n) * dx;
631 ftmp << x_coord << " " << Thyra::get_ele(x,n) << std::endl;
632 }
633 ftmp.close();
634 }
635
636 Teuchos::TimeMonitor::summarize();
637}
638
639
640// ************************************************************
641// ************************************************************
642TEUCHOS_UNIT_TEST(BDF2, VanDerPol)
643{
644 std::vector<RCP<Thyra::VectorBase<double>>> solutions;
645 std::vector<double> StepSize;
646 std::vector<double> ErrorNorm;
647
648 // Read params from .xml file
649 RCP<ParameterList> pList =
650 getParametersFromXmlFile("Tempus_BDF2_VanDerPol.xml");
651 //Set initial time step = 2*dt specified in input file (for convergence study)
652 //
653 RCP<ParameterList> pl = sublist(pList, "Tempus", true);
654 double dt = pl->sublist("Demo Integrator")
655 .sublist("Time Step Control").get<double>("Initial Time Step");
656 dt *= 2.0;
657
658 RCP<ParameterList> vdpm_pl = sublist(pList, "VanDerPolModel", true);
659 const int nTimeStepSizes = vdpm_pl->get<int>("Number of Time Step Sizes", 3);
660 //const int nTimeStepSizes = 5;
661 double order = 0.0;
662
663 for (int n=0; n<nTimeStepSizes; n++) {
664
665 // Setup the VanDerPolModel
666 auto model = rcp(new VanDerPolModel<double>(vdpm_pl));
667
668 // Set the step size
669 dt /= 2;
670 if (n == nTimeStepSizes-1) dt /= 10.0;
671
672 // Setup the Integrator and reset initial time step
673 pl->sublist("Demo Integrator")
674 .sublist("Time Step Control").set("Initial Time Step", dt);
675 RCP<Tempus::IntegratorBasic<double> > integrator =
676 Tempus::createIntegratorBasic<double>(pl, model);
677 order = integrator->getStepper()->getOrder();
678
679 // Integrate to timeMax
680 bool integratorStatus = integrator->advanceTime();
681 TEST_ASSERT(integratorStatus)
682
683 // Test if at 'Final Time'
684 double time = integrator->getTime();
685 double timeFinal =pl->sublist("Demo Integrator")
686 .sublist("Time Step Control").get<double>("Final Time");
687 double tol = 100.0 * std::numeric_limits<double>::epsilon();
688 TEST_FLOATING_EQUALITY(time, timeFinal, tol);
689
690 // Store off the final solution and step size
691 auto solution = Thyra::createMember(model->get_x_space());
692 Thyra::copy(*(integrator->getX()),solution.ptr());
693 solutions.push_back(solution);
694 StepSize.push_back(dt);
695
696 // Output finest temporal solution for plotting
697 // This only works for ONE MPI process
698 if ((n == 0) || (n == nTimeStepSizes-1)) {
699 std::string fname = "Tempus_BDF2_VanDerPol-Ref.dat";
700 if (n == 0) fname = "Tempus_BDF2_VanDerPol.dat";
701 std::ofstream ftmp(fname);
702 RCP<const SolutionHistory<double> > solutionHistory =
703 integrator->getSolutionHistory();
704 int nStates = solutionHistory->getNumStates();
705 for (int i=0; i<nStates; i++) {
706 RCP<const SolutionState<double> > solutionState = (*solutionHistory)[i];
707 RCP<const Thyra::VectorBase<double> > x = solutionState->getX();
708 double ttime = solutionState->getTime();
709 ftmp << ttime << " " << get_ele(*x, 0) << " " << get_ele(*x, 1)
710 << std::endl;
711 }
712 ftmp.close();
713 }
714 }
715
716 // Calculate the error - use the most temporally refined mesh for
717 // the reference solution.
718 auto ref_solution = solutions[solutions.size()-1];
719 std::vector<double> StepSizeCheck;
720 for (std::size_t i=0; i < (solutions.size()-1); ++i) {
721 auto tmp = solutions[i];
722 Thyra::Vp_StV(tmp.ptr(), -1.0, *ref_solution);
723 const double L2norm = Thyra::norm_2(*tmp);
724 StepSizeCheck.push_back(StepSize[i]);
725 ErrorNorm.push_back(L2norm);
726 }
727
728 if (nTimeStepSizes > 2) {
729 // Check the order and intercept
730 double slope = computeLinearRegressionLogLog<double>(StepSizeCheck,ErrorNorm);
731 out << " Stepper = BDF2" << std::endl;
732 out << " =========================" << std::endl;
733 out << " Expected order: " << order << std::endl;
734 out << " Observed order: " << slope << std::endl;
735 out << " =========================" << std::endl;
736 TEST_FLOATING_EQUALITY( slope, order, 0.10 );
737 out << "\n\n ** Slope on BDF2 Method = " << slope
738 << "\n" << std::endl;
739 }
740
741 // Write error data
742 {
743 std::ofstream ftmp("Tempus_BDF2_VanDerPol-Error.dat");
744 double error0 = 0.8*ErrorNorm[0];
745 for (std::size_t n = 0; n < StepSizeCheck.size(); n++) {
746 ftmp << StepSizeCheck[n] << " " << ErrorNorm[n] << " "
747 << error0*(pow(StepSize[n]/StepSize[0],order)) << std::endl;
748 }
749 ftmp.close();
750 }
751
752 Teuchos::TimeMonitor::summarize();
753}
754
755} // 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...
BDF2 (Backward-Difference-Formula-2) time stepper.
TimeStepControl manages the time step size. There several mechanisms that effect the time step size a...
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.