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  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 // ************************************************************
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  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 // ************************************************************
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  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 // ************************************************************
479 TEUCHOS_UNIT_TEST(BDF2, CDR)
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 // ************************************************************
642 TEUCHOS_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
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:355
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.
Definition: Tempus_SolutionState_decl.hpp:137
Tempus_Test::CDR_Model
1D CGFEM model for convection/diffusion/reaction
Definition: CDR_Model_decl.hpp:45
Generated on Fri Sep 29 2023 07:58:46 for Tempus by   doxygen 1.8.5