Tempus  Version of the Day
Time Integration
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
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"
21 #include "../TestUtils/Tempus_ConvergenceTestUtils.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 
37 namespace Tempus_Test {
38 
39 using Teuchos::RCP;
40 using Teuchos::rcp;
41 using Teuchos::rcp_const_cast;
42 using Teuchos::ParameterList;
43 using Teuchos::sublist;
44 using Teuchos::getParametersFromXmlFile;
45 
46 using Tempus::IntegratorBasic;
47 using Tempus::SolutionHistory;
48 using Tempus::SolutionState;
49 
50 
51 // ************************************************************
52 // ************************************************************
53 TEUCHOS_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  std::cout << std::endl;
76  std::cout << "stepperPL -------------- \n" << *stepperPL << std::endl;
77  std::cout << "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, "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  std::cout << std::endl;
94  std::cout << "stepperPL -------------- \n" << *stepperPL << std::endl;
95  std::cout << "defaultPL -------------- \n" << *defaultPL << std::endl;
96  }
97  TEST_ASSERT(pass)
98  }
99 }
100 
101 
102 // ************************************************************
103 // ************************************************************
104 TEUCHOS_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  std::cout << " Stepper = " << stepper->description()
194  << "\n with " << option << std::endl;
195  std::cout << " =========================" << std::endl;
196  std::cout << " Exact solution : " << get_ele(*(x_exact), 0) << " "
197  << get_ele(*(x_exact), 1) << std::endl;
198  std::cout << " Computed solution: " << get_ele(*(x ), 0) << " "
199  << get_ele(*(x ), 1) << std::endl;
200  std::cout << " Difference : " << get_ele(*(xdiff ), 0) << " "
201  << get_ele(*(xdiff ), 1) << std::endl;
202  std::cout << " =========================" << 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 // ************************************************************
211 TEUCHOS_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 // ************************************************************
332 TEUCHOS_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  RCP<ParameterList> pl = sublist(pList, "Tempus", true);
344  double dt = pList->sublist("Tempus")
345  .sublist("Default Integrator")
346  .sublist("Time Step Control").get<double>("Initial Time Step");
347  dt *= 2.0;
348 
349  // Setup the SinCosModel
350  RCP<ParameterList> scm_pl = sublist(pList, "SinCosModel", true);
351  const int nTimeStepSizes = scm_pl->get<int>("Number of Time Step Sizes", 7);
352  std::string output_file_string =
353  scm_pl->get<std::string>("Output File Name", "Tempus_BDF2_SinCos");
354  std::string output_file_name = output_file_string + ".dat";
355  std::string err_out_file_name = output_file_string + "-Error.dat";
356  double time = 0.0;
357  for (int n=0; n<nTimeStepSizes; n++) {
358 
359  auto model = rcp(new SinCosModel<double>(scm_pl));
360 
361  dt /= 2;
362 
363  RCP<ParameterList> tempusPL =
364  getParametersFromXmlFile("Tempus_BDF2_SinCos_AdaptDt.xml");
365  RCP<ParameterList> pl = sublist(tempusPL, "Tempus", true);
366 
367  // Setup the Integrator and reset initial time step
368  pl->sublist("Default Integrator")
369  .sublist("Time Step Control").set("Initial Time Step", dt/4.0);
370  // Ensure time step does not get larger than the initial time step size,
371  // as that would mess up the convergence rates.
372  pl->sublist("Default Integrator")
373  .sublist("Time Step Control").set("Maximum Time Step", dt);
374  // Ensure time step does not get too small and therefore too many steps.
375  pl->sublist("Default Integrator")
376  .sublist("Time Step Control").set("Minimum Time Step", dt/4.0);
377  // For the SinCos problem eta is directly related to dt
378  pl->sublist("Default Integrator")
379  .sublist("Time Step Control")
380  .sublist("Time Step Control Strategy")
381  .set("Minimum Value Monitoring Function", dt*0.99);
382  integrator = Tempus::createIntegratorBasic<double>(pl, model);
383 
384  // Initial Conditions
385  // During the Integrator construction, the initial SolutionState
386  // is set by default to model->getNominalVales().get_x(). However,
387  // the application can set it also by integrator->initializeSolutionHistory.
388  RCP<Thyra::VectorBase<double> > x0 =
389  model->getNominalValues().get_x()->clone_v();
390  integrator->initializeSolutionHistory(0.0, x0);
391 
392  // Integrate to timeMax
393  bool integratorStatus = integrator->advanceTime();
394  TEST_ASSERT(integratorStatus)
395 
396  // Test if at 'Final Time'
397  time = integrator->getTime();
398  double timeFinal =pl->sublist("Default Integrator")
399  .sublist("Time Step Control").get<double>("Final Time");
400  TEST_FLOATING_EQUALITY(time, timeFinal, 1.0e-14);
401 
402  // Time-integrated solution and the exact solution
403  RCP<Thyra::VectorBase<double> > x = integrator->getX();
404  RCP<const Thyra::VectorBase<double> > x_exact =
405  model->getExactSolution(time).get_x();
406 
407  // Plot sample solution and exact solution
408  if (n == 0) {
409  std::ofstream ftmp(output_file_name);
410  //Warning: the following assumes serial run
411  FILE *gold_file = fopen("Tempus_BDF2_SinCos_AdaptDt_gold.dat", "r");
412  RCP<const SolutionHistory<double> > solutionHistory =
413  integrator->getSolutionHistory();
414  RCP<const Thyra::VectorBase<double> > x_exact_plot;
415  for (int i=0; i<solutionHistory->getNumStates(); i++) {
416  char time_gold_char[100];
417  fgets(time_gold_char, 100, gold_file);
418  double time_gold;
419  sscanf(time_gold_char, "%lf", &time_gold);
420  RCP<const SolutionState<double> > solutionState = (*solutionHistory)[i];
421  double time_i = solutionState->getTime();
422  //Throw error if time does not match time in gold file to specified tolerance
423  TEST_FLOATING_EQUALITY( time_i, time_gold, 1.0e-5 );
424  RCP<const Thyra::VectorBase<double> > x_plot = solutionState->getX();
425  x_exact_plot = model->getExactSolution(time_i).get_x();
426  ftmp << time_i << " "
427  << get_ele(*(x_plot), 0) << " "
428  << get_ele(*(x_plot), 1) << " "
429  << get_ele(*(x_exact_plot), 0) << " "
430  << get_ele(*(x_exact_plot), 1) << std::endl;
431  }
432  ftmp.close();
433  }
434 
435  // Store off the final solution and step size
436  StepSize.push_back(dt);
437  auto solution = Thyra::createMember(model->get_x_space());
438  Thyra::copy(*(integrator->getX()),solution.ptr());
439  solutions.push_back(solution);
440  auto solutionDot = Thyra::createMember(model->get_x_space());
441  Thyra::copy(*(integrator->getXDot()),solutionDot.ptr());
442  solutionsDot.push_back(solutionDot);
443  if (n == nTimeStepSizes-1) { // Add exact solution last in vector.
444  StepSize.push_back(0.0);
445  auto solutionExact = Thyra::createMember(model->get_x_space());
446  Thyra::copy(*(model->getExactSolution(time).get_x()),solutionExact.ptr());
447  solutions.push_back(solutionExact);
448  auto solutionDotExact = Thyra::createMember(model->get_x_space());
449  Thyra::copy(*(model->getExactSolution(time).get_x_dot()),
450  solutionDotExact.ptr());
451  solutionsDot.push_back(solutionDotExact);
452  }
453  }
454 
455  // Check the order and intercept
456  if (nTimeStepSizes > 1) {
457  double xSlope = 0.0;
458  double xDotSlope = 0.0;
459  std::vector<double> xErrorNorm;
460  std::vector<double> xDotErrorNorm;
461  RCP<Tempus::Stepper<double> > stepper = integrator->getStepper();
462  //double order = stepper->getOrder();
463  writeOrderError("Tempus_BDF2_SinCos-Error.dat",
464  stepper, StepSize,
465  solutions, xErrorNorm, xSlope,
466  solutionsDot, xDotErrorNorm, xDotSlope);
467 
468  TEST_FLOATING_EQUALITY( xSlope, 1.932, 0.01 );
469  TEST_FLOATING_EQUALITY( xDotSlope, 1.932, 0.01 );
470  TEST_FLOATING_EQUALITY( xErrorNorm[0], 0.000192591, 1.0e-4 );
471  TEST_FLOATING_EQUALITY( xDotErrorNorm[0], 0.000192591, 1.0e-4 );
472  }
473 
474  Teuchos::TimeMonitor::summarize();
475 }
476 
477 
478 // ************************************************************
479 // ************************************************************
480 TEUCHOS_UNIT_TEST(BDF2, CDR)
481 {
482  // Create a communicator for Epetra objects
483  RCP<Epetra_Comm> comm;
484 #ifdef Tempus_ENABLE_MPI
485  comm = rcp(new Epetra_MpiComm(MPI_COMM_WORLD));
486 #else
487  comm = rcp(new Epetra_SerialComm);
488 #endif
489 
490  RCP<Tempus::IntegratorBasic<double> > integrator;
491  std::vector<RCP<Thyra::VectorBase<double>>> solutions;
492  std::vector<RCP<Thyra::VectorBase<double>>> solutionsDot;
493  std::vector<double> StepSize;
494 
495  // Read params from .xml file
496  RCP<ParameterList> pList =
497  getParametersFromXmlFile("Tempus_BDF2_CDR.xml");
498  //Set initial time step = 2*dt specified in input file (for convergence study)
499  //
500  RCP<ParameterList> pl = sublist(pList, "Tempus", true);
501  double dt = pl->sublist("Demo Integrator")
502  .sublist("Time Step Control").get<double>("Initial Time Step");
503  dt *= 2.0;
504  RCP<ParameterList> model_pl = sublist(pList, "CDR Model", true);
505 
506  const int nTimeStepSizes = model_pl->get<int>("Number of Time Step Sizes", 5);
507 
508  for (int n=0; n<nTimeStepSizes; n++) {
509 
510  // Create CDR Model
511  const int num_elements = model_pl->get<int>("num elements");
512  const double left_end = model_pl->get<double>("left end");
513  const double right_end = model_pl->get<double>("right end");
514  const double a_convection = model_pl->get<double>("a (convection)");
515  const double k_source = model_pl->get<double>("k (source)");
516 
517  auto model = rcp(new Tempus_Test::CDR_Model<double>(comm,
518  num_elements,
519  left_end,
520  right_end,
521  a_convection,
522  k_source));
523 
524  // Set the factory
525  ::Stratimikos::DefaultLinearSolverBuilder builder;
526 
527  auto p = rcp(new ParameterList);
528  p->set("Linear Solver Type", "Belos");
529  p->set("Preconditioner Type", "None");
530  builder.setParameterList(p);
531 
532  RCP< ::Thyra::LinearOpWithSolveFactoryBase<double> >
533  lowsFactory = builder.createLinearSolveStrategy("");
534 
535  model->set_W_factory(lowsFactory);
536 
537  // Set the step size
538  dt /= 2;
539 
540  // Setup the Integrator and reset initial time step
541  pl->sublist("Demo Integrator")
542  .sublist("Time Step Control").set("Initial Time Step", dt);
543  integrator = Tempus::createIntegratorBasic<double>(pl, model);
544 
545  // Integrate to timeMax
546  bool integratorStatus = integrator->advanceTime();
547  TEST_ASSERT(integratorStatus)
548 
549  // Test if at 'Final Time'
550  double time = integrator->getTime();
551  double timeFinal =pl->sublist("Demo Integrator")
552  .sublist("Time Step Control").get<double>("Final Time");
553  double tol = 100.0 * std::numeric_limits<double>::epsilon();
554  TEST_FLOATING_EQUALITY(time, timeFinal, tol);
555 
556  // Store off the final solution and step size
557  StepSize.push_back(dt);
558  auto solution = Thyra::createMember(model->get_x_space());
559  Thyra::copy(*(integrator->getX()),solution.ptr());
560  solutions.push_back(solution);
561  auto solutionDot = Thyra::createMember(model->get_x_space());
562  Thyra::copy(*(integrator->getXDot()),solutionDot.ptr());
563  solutionsDot.push_back(solutionDot);
564 
565  // Output finest temporal solution for plotting
566  // This only works for ONE MPI process
567  if ((n == nTimeStepSizes-1) && (comm->NumProc() == 1)) {
568  std::ofstream ftmp("Tempus_BDF2_CDR.dat");
569  ftmp << "TITLE=\"BDF2 Solution to CDR\"\n"
570  << "VARIABLES=\"z\",\"T\"\n";
571  const double dx = std::fabs(left_end-right_end) /
572  static_cast<double>(num_elements);
573  RCP<const SolutionHistory<double> > solutionHistory =
574  integrator->getSolutionHistory();
575  int nStates = solutionHistory->getNumStates();
576  for (int i=0; i<nStates; i++) {
577  RCP<const SolutionState<double> > solutionState = (*solutionHistory)[i];
578  RCP<const Thyra::VectorBase<double> > x = solutionState->getX();
579  double ttime = solutionState->getTime();
580  ftmp << "ZONE T=\"Time="<<ttime<<"\", I="
581  <<num_elements+1<<", F=BLOCK\n";
582  for (int j = 0; j < num_elements+1; j++) {
583  const double x_coord = left_end + static_cast<double>(j) * dx;
584  ftmp << x_coord << " ";
585  }
586  ftmp << std::endl;
587  for (int j=0; j<num_elements+1; j++) ftmp << get_ele(*x, j) << " ";
588  ftmp << std::endl;
589  }
590  ftmp.close();
591  }
592  }
593 
594  // Check the order and intercept
595  if (nTimeStepSizes > 2) {
596  double xSlope = 0.0;
597  double xDotSlope = 0.0;
598  std::vector<double> xErrorNorm;
599  std::vector<double> xDotErrorNorm;
600  RCP<Tempus::Stepper<double> > stepper = integrator->getStepper();
601  double order = stepper->getOrder();
602  writeOrderError("Tempus_BDF2_CDR-Error.dat",
603  stepper, StepSize,
604  solutions, xErrorNorm, xSlope,
605  solutionsDot, xDotErrorNorm, xDotSlope);
606  TEST_FLOATING_EQUALITY( xSlope, order, 0.35 );
607  TEST_COMPARE(xSlope, >, 0.95);
608  TEST_FLOATING_EQUALITY( xDotSlope, order, 0.35 );
609  TEST_COMPARE(xDotSlope, >, 0.95);
610 
611  TEST_FLOATING_EQUALITY( xErrorNorm[0], 0.0145747, 1.0e-4 );
612  TEST_FLOATING_EQUALITY( xDotErrorNorm[0], 0.0563621, 1.0e-4 );
613  }
614 
615  // Write fine mesh solution at final time
616  // This only works for ONE MPI process
617  if (comm->NumProc() == 1) {
618  RCP<ParameterList> pListCDR =
619  getParametersFromXmlFile("Tempus_BDF2_CDR.xml");
620  RCP<ParameterList> model_pl_CDR = sublist(pListCDR, "CDR Model", true);
621  const int num_elements = model_pl_CDR->get<int>("num elements");
622  const double left_end = model_pl_CDR->get<double>("left end");
623  const double right_end = model_pl_CDR->get<double>("right end");
624 
625  const Thyra::VectorBase<double>& x = *(solutions[solutions.size()-1]);
626 
627  std::ofstream ftmp("Tempus_BDF2_CDR-Solution.dat");
628  for (int n = 0; n < num_elements+1; n++) {
629  const double dx = std::fabs(left_end-right_end) /
630  static_cast<double>(num_elements);
631  const double x_coord = left_end + static_cast<double>(n) * dx;
632  ftmp << x_coord << " " << Thyra::get_ele(x,n) << std::endl;
633  }
634  ftmp.close();
635  }
636 
637  Teuchos::TimeMonitor::summarize();
638 }
639 
640 
641 // ************************************************************
642 // ************************************************************
643 TEUCHOS_UNIT_TEST(BDF2, VanDerPol)
644 {
645  std::vector<RCP<Thyra::VectorBase<double>>> solutions;
646  std::vector<double> StepSize;
647  std::vector<double> ErrorNorm;
648 
649  // Read params from .xml file
650  RCP<ParameterList> pList =
651  getParametersFromXmlFile("Tempus_BDF2_VanDerPol.xml");
652  //Set initial time step = 2*dt specified in input file (for convergence study)
653  //
654  RCP<ParameterList> pl = sublist(pList, "Tempus", true);
655  double dt = pl->sublist("Demo Integrator")
656  .sublist("Time Step Control").get<double>("Initial Time Step");
657  dt *= 2.0;
658 
659  RCP<ParameterList> vdpm_pl = sublist(pList, "VanDerPolModel", true);
660  const int nTimeStepSizes = vdpm_pl->get<int>("Number of Time Step Sizes", 3);
661  //const int nTimeStepSizes = 5;
662  double order = 0.0;
663 
664  for (int n=0; n<nTimeStepSizes; n++) {
665 
666  // Setup the VanDerPolModel
667  auto model = rcp(new VanDerPolModel<double>(vdpm_pl));
668 
669  // Set the step size
670  dt /= 2;
671  if (n == nTimeStepSizes-1) dt /= 10.0;
672 
673  // Setup the Integrator and reset initial time step
674  pl->sublist("Demo Integrator")
675  .sublist("Time Step Control").set("Initial Time Step", dt);
676  RCP<Tempus::IntegratorBasic<double> > integrator =
677  Tempus::createIntegratorBasic<double>(pl, model);
678  order = integrator->getStepper()->getOrder();
679 
680  // Integrate to timeMax
681  bool integratorStatus = integrator->advanceTime();
682  TEST_ASSERT(integratorStatus)
683 
684  // Test if at 'Final Time'
685  double time = integrator->getTime();
686  double timeFinal =pl->sublist("Demo Integrator")
687  .sublist("Time Step Control").get<double>("Final Time");
688  double tol = 100.0 * std::numeric_limits<double>::epsilon();
689  TEST_FLOATING_EQUALITY(time, timeFinal, tol);
690 
691  // Store off the final solution and step size
692  auto solution = Thyra::createMember(model->get_x_space());
693  Thyra::copy(*(integrator->getX()),solution.ptr());
694  solutions.push_back(solution);
695  StepSize.push_back(dt);
696 
697  // Output finest temporal solution for plotting
698  // This only works for ONE MPI process
699  if ((n == 0) || (n == nTimeStepSizes-1)) {
700  std::string fname = "Tempus_BDF2_VanDerPol-Ref.dat";
701  if (n == 0) fname = "Tempus_BDF2_VanDerPol.dat";
702  std::ofstream ftmp(fname);
703  RCP<const SolutionHistory<double> > solutionHistory =
704  integrator->getSolutionHistory();
705  int nStates = solutionHistory->getNumStates();
706  for (int i=0; i<nStates; i++) {
707  RCP<const SolutionState<double> > solutionState = (*solutionHistory)[i];
708  RCP<const Thyra::VectorBase<double> > x = solutionState->getX();
709  double ttime = solutionState->getTime();
710  ftmp << ttime << " " << get_ele(*x, 0) << " " << get_ele(*x, 1)
711  << std::endl;
712  }
713  ftmp.close();
714  }
715  }
716 
717  // Calculate the error - use the most temporally refined mesh for
718  // the reference solution.
719  auto ref_solution = solutions[solutions.size()-1];
720  std::vector<double> StepSizeCheck;
721  for (std::size_t i=0; i < (solutions.size()-1); ++i) {
722  auto tmp = solutions[i];
723  Thyra::Vp_StV(tmp.ptr(), -1.0, *ref_solution);
724  const double L2norm = Thyra::norm_2(*tmp);
725  StepSizeCheck.push_back(StepSize[i]);
726  ErrorNorm.push_back(L2norm);
727  }
728 
729  if (nTimeStepSizes > 2) {
730  // Check the order and intercept
731  double slope = computeLinearRegressionLogLog<double>(StepSizeCheck,ErrorNorm);
732  std::cout << " Stepper = BDF2" << std::endl;
733  std::cout << " =========================" << std::endl;
734  std::cout << " Expected order: " << order << std::endl;
735  std::cout << " Observed order: " << slope << std::endl;
736  std::cout << " =========================" << std::endl;
737  TEST_FLOATING_EQUALITY( slope, order, 0.10 );
738  out << "\n\n ** Slope on BDF2 Method = " << slope
739  << "\n" << std::endl;
740  }
741 
742  // Write error data
743  {
744  std::ofstream ftmp("Tempus_BDF2_VanDerPol-Error.dat");
745  double error0 = 0.8*ErrorNorm[0];
746  for (std::size_t n = 0; n < StepSizeCheck.size(); n++) {
747  ftmp << StepSizeCheck[n] << " " << ErrorNorm[n] << " "
748  << error0*(pow(StepSize[n]/StepSize[0],order)) << std::endl;
749  }
750  ftmp.close();
751  }
752 
753  Teuchos::TimeMonitor::summarize();
754 }
755 
756 } // namespace Tempus_Test
Tempus::StepperBDF2
BDF2 (Backward-Difference-Formula-2) time stepper.
Definition: Tempus_StepperBDF2_decl.hpp:99
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:479
Teuchos::ParameterList::get
T & get(const std::string &name, T def_value)
Teuchos::ParameterList::set
ParameterList & set(std::string const &name, T const &value, std::string const &docString="", RCP< const ParameterEntryValidator > const &validator=null)
TEST_COMPARE
#define TEST_COMPARE(v1, comp, v2)
TEST_FLOATING_EQUALITY
#define TEST_FLOATING_EQUALITY(v1, v2, tol)
TEST_ASSERT
#define TEST_ASSERT(v1)
Teuchos::RCP::get
T * get() const
Tempus_Test::SinCosModel
Sine-Cosine model problem from Rythmos. This is a canonical Sine-Cosine differential equation with a...
Definition: SinCosModel_decl.hpp:91
Tempus_Test::writeSolution
void writeSolution(const std::string filename, Teuchos::RCP< const Tempus::SolutionHistory< Scalar > > solutionHistory)
Definition: Tempus_ConvergenceTestUtils.hpp:311
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::TEUCHOS_UNIT_TEST
TEUCHOS_UNIT_TEST(BackwardEuler, SinCos_ASA)
Definition: Tempus_BackwardEuler_ASA.cpp:47
Teuchos::rcp
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
Teuchos::TimeMonitor::summarize
static void summarize(Ptr< const Comm< int > > comm, std::ostream &out=std::cout, const bool alwaysWriteLocal=false, const bool writeGlobalStats=true, const bool writeZeroTimers=true, const ECounterSetOp setOp=Intersection, const std::string &filter="", const bool ignoreZeroTimers=false)
Teuchos::RCP::ptr
Ptr< T > ptr() const
Teuchos_UnitTestHarness.hpp
Tempus::TimeStepControl
TimeStepControl manages the time step size. There several mechanisms that effect the time step size a...
Definition: Tempus_Integrator.hpp:26
Tempus::SolutionHistory
SolutionHistory is basically a container of SolutionStates. SolutionHistory maintains a collection of...
Definition: Tempus_Integrator.hpp:25
Teuchos_XMLParameterListHelpers.hpp
Tempus::STORAGE_TYPE_STATIC
Keep a fix number of states.
Definition: Tempus_SolutionHistory_decl.hpp:23
Tempus_Test::VanDerPolModel
van der Pol model problem for nonlinear electrical circuit.
Definition: VanDerPolModel_decl.hpp:111
Teuchos::ParameterList::sublist
ParameterList & sublist(const std::string &name, bool mustAlreadyExist=false, const std::string &docString="")
Teuchos_TimeMonitor.hpp
Tempus::SolutionState
Solution state for integrators and steppers. SolutionState contains the metadata for solutions and th...
Definition: Tempus_SolutionState_decl.hpp:56
Tempus_Test::CDR_Model
1D CGFEM model for convection/diffusion/reaction
Definition: CDR_Model_decl.hpp:45
Generated by   doxygen 1.8.5