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