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