Tempus  Version of the Day
Time Integration
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Macros Pages
Tempus_BackwardEulerTest.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_StepperBackwardEuler.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 <vector>
33 #include <fstream>
34 #include <sstream>
35 #include <limits>
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 
49 
50 // Comment out any of the following tests to exclude from build/run.
51 #define TEST_PARAMETERLIST
52 #define TEST_CONSTRUCTING_FROM_DEFAULTS
53 #define TEST_SINCOS
54 #define TEST_CDR
55 #define TEST_VANDERPOL
56 #define TEST_OPT_INTERFACE
57 
58 
59 #ifdef TEST_PARAMETERLIST
60 // ************************************************************
61 // ************************************************************
62 TEUCHOS_UNIT_TEST(BackwardEuler, ParameterList)
63 {
64  // Read params from .xml file
65  RCP<ParameterList> pList =
66  getParametersFromXmlFile("Tempus_BackwardEuler_SinCos.xml");
67 
68  // Setup the SinCosModel
69  RCP<ParameterList> scm_pl = sublist(pList, "SinCosModel", true);
70  auto model = 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 Predictor for comparison
81  stepperPL->remove("Predictor Name");
82  stepperPL->remove("Default Predictor");
83  RCP<ParameterList> defaultPL =
84  integrator->getStepper()->getDefaultParameters();
85 
86  bool pass = haveSameValues(*stepperPL, *defaultPL, true);
87  if (!pass) {
88  std::cout << std::endl;
89  std::cout << "stepperPL -------------- \n" << *stepperPL << std::endl;
90  std::cout << "defaultPL -------------- \n" << *defaultPL << std::endl;
91  }
92  TEST_ASSERT(pass)
93  }
94 
95  // Test constructor IntegratorBasic(model, stepperType)
96  {
97  RCP<Tempus::IntegratorBasic<double> > integrator =
98  Tempus::integratorBasic<double>(model, "Backward Euler");
99 
100  RCP<ParameterList> stepperPL = sublist(tempusPL, "Default Stepper", true);
101  RCP<ParameterList> defaultPL =
102  integrator->getStepper()->getDefaultParameters();
103 
104  bool pass = haveSameValues(*stepperPL, *defaultPL, true);
105  if (!pass) {
106  std::cout << std::endl;
107  std::cout << "stepperPL -------------- \n" << *stepperPL << std::endl;
108  std::cout << "defaultPL -------------- \n" << *defaultPL << std::endl;
109  }
110  TEST_ASSERT(pass)
111  }
112 }
113 #endif // TEST_PARAMETERLIST
114 
115 
116 #ifdef TEST_CONSTRUCTING_FROM_DEFAULTS
117 // ************************************************************
118 // ************************************************************
119 TEUCHOS_UNIT_TEST(BackwardEuler, ConstructingFromDefaults)
120 {
121  double dt = 0.1;
122 
123  // Read params from .xml file
124  RCP<ParameterList> pList =
125  getParametersFromXmlFile("Tempus_BackwardEuler_SinCos.xml");
126  RCP<ParameterList> pl = sublist(pList, "Tempus", true);
127 
128  // Setup the SinCosModel
129  RCP<ParameterList> scm_pl = sublist(pList, "SinCosModel", true);
130  //RCP<SinCosModel<double> > model = sineCosineModel(scm_pl);
131  auto model = rcp(new SinCosModel<double>(scm_pl));
132 
133  // Setup Stepper for field solve ----------------------------
134  auto stepper = rcp(new Tempus::StepperBackwardEuler<double>());
135  stepper->setModel(model);
136  stepper->initialize();
137 
138  // Setup TimeStepControl ------------------------------------
139  auto timeStepControl = rcp(new Tempus::TimeStepControl<double>());
140  ParameterList tscPL = pl->sublist("Default Integrator")
141  .sublist("Time Step Control");
142  timeStepControl->setStepType (tscPL.get<std::string>("Integrator Step Type"));
143  timeStepControl->setInitIndex(tscPL.get<int> ("Initial Time Index"));
144  timeStepControl->setInitTime (tscPL.get<double>("Initial Time"));
145  timeStepControl->setFinalTime(tscPL.get<double>("Final Time"));
146  timeStepControl->setInitTimeStep(dt);
147  timeStepControl->initialize();
148 
149  // Setup initial condition SolutionState --------------------
150  Thyra::ModelEvaluatorBase::InArgs<double> inArgsIC =
151  stepper->getModel()->getNominalValues();
152  auto icSolution =
153  rcp_const_cast<Thyra::VectorBase<double> > (inArgsIC.get_x());
154  auto icState = 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 ------------------------------------
163  solutionHistory->setName("Forward States");
165  solutionHistory->setStorageLimit(2);
166  solutionHistory->addState(icState);
167 
168  // Setup Integrator -----------------------------------------
169  RCP<Tempus::IntegratorBasic<double> > integrator =
170  Tempus::integratorBasic<double>();
171  integrator->setStepperWStepper(stepper);
172  integrator->setTimeStepControl(timeStepControl);
173  integrator->setSolutionHistory(solutionHistory);
174  //integrator->setObserver(...);
175  integrator->initialize();
176 
177 
178  // Integrate to timeMax
179  bool integratorStatus = integrator->advanceTime();
180  TEST_ASSERT(integratorStatus)
181 
182 
183  // Test if at 'Final Time'
184  double time = integrator->getTime();
185  double timeFinal =pl->sublist("Default Integrator")
186  .sublist("Time Step Control").get<double>("Final Time");
187  TEST_FLOATING_EQUALITY(time, timeFinal, 1.0e-14);
188 
189  // Time-integrated solution and the exact solution
190  RCP<Thyra::VectorBase<double> > x = integrator->getX();
191  RCP<const Thyra::VectorBase<double> > x_exact =
192  model->getExactSolution(time).get_x();
193 
194  // Calculate the error
195  RCP<Thyra::VectorBase<double> > xdiff = x->clone_v();
196  Thyra::V_StVpStV(xdiff.ptr(), 1.0, *x_exact, -1.0, *(x));
197 
198  // Check the order and intercept
199  std::cout << " Stepper = BackwardEuler" << std::endl;
200  std::cout << " =========================" << std::endl;
201  std::cout << " Exact solution : " << get_ele(*(x_exact), 0) << " "
202  << get_ele(*(x_exact), 1) << std::endl;
203  std::cout << " Computed solution: " << get_ele(*(x ), 0) << " "
204  << get_ele(*(x ), 1) << std::endl;
205  std::cout << " Difference : " << get_ele(*(xdiff ), 0) << " "
206  << get_ele(*(xdiff ), 1) << std::endl;
207  std::cout << " =========================" << std::endl;
208  TEST_FLOATING_EQUALITY(get_ele(*(x), 0), 0.798923, 1.0e-4 );
209  TEST_FLOATING_EQUALITY(get_ele(*(x), 1), 0.516729, 1.0e-4 );
210 }
211 #endif // TEST_CONSTRUCTING_FROM_DEFAULTS
212 
213 
214 #ifdef TEST_SINCOS
215 // ************************************************************
216 // ************************************************************
217 TEUCHOS_UNIT_TEST(BackwardEuler, SinCos)
218 {
219  RCP<Tempus::IntegratorBasic<double> > integrator;
220  std::vector<RCP<Thyra::VectorBase<double>>> solutions;
221  std::vector<RCP<Thyra::VectorBase<double>>> solutionsDot;
222  std::vector<double> StepSize;
223  std::vector<double> xErrorNorm;
224  std::vector<double> xDotErrorNorm;
225  const int nTimeStepSizes = 7;
226  double dt = 0.2;
227  double time = 0.0;
228  for (int n=0; n<nTimeStepSizes; n++) {
229 
230  // Read params from .xml file
231  RCP<ParameterList> pList =
232  getParametersFromXmlFile("Tempus_BackwardEuler_SinCos.xml");
233 
234  //std::ofstream ftmp("PL.txt");
235  //pList->print(ftmp);
236  //ftmp.close();
237 
238  // Setup the SinCosModel
239  RCP<ParameterList> scm_pl = sublist(pList, "SinCosModel", true);
240  //RCP<SinCosModel<double> > model = sineCosineModel(scm_pl);
241  auto model = rcp(new SinCosModel<double>(scm_pl));
242 
243  dt /= 2;
244 
245  // Setup the Integrator and reset initial time step
246  RCP<ParameterList> pl = sublist(pList, "Tempus", true);
247  pl->sublist("Default Integrator")
248  .sublist("Time Step Control").set("Initial Time Step", dt);
249  integrator = Tempus::integratorBasic<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("Tempus_BackwardEuler_SinCos.dat", 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 = rcp(new Tempus::SolutionState<double>(
279  model->getExactSolution(time_i).get_x(),
280  model->getExactSolution(time_i).get_x_dot()));
281  state->setTime((*solutionHistory)[i]->getTime());
282  solnHistExact->addState(state);
283  }
284  writeSolution("Tempus_BackwardEuler_SinCos-Ref.dat", 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  double xSlope = 0.0;
309  double xDotSlope = 0.0;
310  RCP<Tempus::Stepper<double> > stepper = integrator->getStepper();
311  double order = stepper->getOrder();
312  writeOrderError("Tempus_BackwardEuler_SinCos-Error.dat",
313  stepper, StepSize,
314  solutions, xErrorNorm, xSlope,
315  solutionsDot, xDotErrorNorm, xDotSlope);
316 
317  TEST_FLOATING_EQUALITY( xSlope, order, 0.01 );
318  TEST_FLOATING_EQUALITY( xErrorNorm[0], 0.0486418, 1.0e-4 );
319  TEST_FLOATING_EQUALITY( xDotSlope, order, 0.01 );
320  TEST_FLOATING_EQUALITY( xDotErrorNorm[0], 0.0486418, 1.0e-4 );
321 
322  Teuchos::TimeMonitor::summarize();
323 }
324 #endif // TEST_SINCOS
325 
326 
327 #ifdef TEST_CDR
328 // ************************************************************
329 // ************************************************************
330 TEUCHOS_UNIT_TEST(BackwardEuler, CDR)
331 {
332  // Create a communicator for Epetra objects
333  RCP<Epetra_Comm> comm;
334 #ifdef Tempus_ENABLE_MPI
335  comm = rcp(new Epetra_MpiComm(MPI_COMM_WORLD));
336 #else
337  comm = rcp(new Epetra_SerialComm);
338 #endif
339 
340  RCP<Tempus::IntegratorBasic<double> > integrator;
341  std::vector<RCP<Thyra::VectorBase<double>>> solutions;
342  std::vector<RCP<Thyra::VectorBase<double>>> solutionsDot;
343  std::vector<double> StepSize;
344  std::vector<double> xErrorNorm;
345  std::vector<double> xDotErrorNorm;
346  const int nTimeStepSizes = 5;
347  double dt = 0.2;
348  for (int n=0; n<nTimeStepSizes; n++) {
349 
350  // Read params from .xml file
351  RCP<ParameterList> pList =
352  getParametersFromXmlFile("Tempus_BackwardEuler_CDR.xml");
353 
354  // Create CDR Model
355  RCP<ParameterList> model_pl = sublist(pList, "CDR Model", true);
356  const int num_elements = model_pl->get<int>("num elements");
357  const double left_end = model_pl->get<double>("left end");
358  const double right_end = model_pl->get<double>("right end");
359  const double a_convection = model_pl->get<double>("a (convection)");
360  const double k_source = model_pl->get<double>("k (source)");
361 
362  auto model = rcp(new Tempus_Test::CDR_Model<double>(comm,
363  num_elements,
364  left_end,
365  right_end,
366  a_convection,
367  k_source));
368 
369  // Set the factory
370  ::Stratimikos::DefaultLinearSolverBuilder builder;
371 
372  auto p = rcp(new ParameterList);
373  p->set("Linear Solver Type", "Belos");
374  p->set("Preconditioner Type", "None");
375  builder.setParameterList(p);
376 
377  RCP< ::Thyra::LinearOpWithSolveFactoryBase<double> >
378  lowsFactory = builder.createLinearSolveStrategy("");
379 
380  model->set_W_factory(lowsFactory);
381 
382  // Set the step size
383  dt /= 2;
384 
385  // Setup the Integrator and reset initial time step
386  RCP<ParameterList> pl = sublist(pList, "Tempus", true);
387  pl->sublist("Demo Integrator")
388  .sublist("Time Step Control").set("Initial Time Step", dt);
389  integrator = Tempus::integratorBasic<double>(pl, model);
390 
391  // Integrate to timeMax
392  bool integratorStatus = integrator->advanceTime();
393  TEST_ASSERT(integratorStatus)
394 
395  // Test if at 'Final Time'
396  double time = integrator->getTime();
397  double timeFinal =pl->sublist("Demo Integrator")
398  .sublist("Time Step Control").get<double>("Final Time");
399  double tol = 100.0 * std::numeric_limits<double>::epsilon();
400  TEST_FLOATING_EQUALITY(time, timeFinal, tol);
401 
402  // Store off the final solution and step size
403  StepSize.push_back(dt);
404  auto solution = Thyra::createMember(model->get_x_space());
405  Thyra::copy(*(integrator->getX()),solution.ptr());
406  solutions.push_back(solution);
407  auto solutionDot = Thyra::createMember(model->get_x_space());
408  Thyra::copy(*(integrator->getXdot()),solutionDot.ptr());
409  solutionsDot.push_back(solutionDot);
410 
411  // Output finest temporal solution for plotting
412  // This only works for ONE MPI process
413  if ((n == nTimeStepSizes-1) && (comm->NumProc() == 1)) {
414  std::ofstream ftmp("Tempus_BackwardEuler_CDR.dat");
415  ftmp << "TITLE=\"Backward Euler Solution to CDR\"\n"
416  << "VARIABLES=\"z\",\"T\"\n";
417  const double dx = std::fabs(left_end-right_end) /
418  static_cast<double>(num_elements);
419  RCP<const SolutionHistory<double> > solutionHistory =
420  integrator->getSolutionHistory();
421  int nStates = solutionHistory->getNumStates();
422  for (int i=0; i<nStates; i++) {
423  RCP<const SolutionState<double> > solutionState = (*solutionHistory)[i];
424  RCP<const Thyra::VectorBase<double> > x = solutionState->getX();
425  double ttime = solutionState->getTime();
426  ftmp << "ZONE T=\"Time="<<ttime<<"\", I="
427  <<num_elements+1<<", F=BLOCK\n";
428  for (int j = 0; j < num_elements+1; j++) {
429  const double x_coord = left_end + static_cast<double>(j) * dx;
430  ftmp << x_coord << " ";
431  }
432  ftmp << std::endl;
433  for (int j=0; j<num_elements+1; j++) ftmp << get_ele(*x, j) << " ";
434  ftmp << std::endl;
435  }
436  ftmp.close();
437  }
438  }
439 
440  // Check the order and intercept
441  double xSlope = 0.0;
442  double xDotSlope = 0.0;
443  RCP<Tempus::Stepper<double> > stepper = integrator->getStepper();
444  writeOrderError("Tempus_BackwardEuler_CDR-Error.dat",
445  stepper, StepSize,
446  solutions, xErrorNorm, xSlope,
447  solutionsDot, xDotErrorNorm, xDotSlope);
448 
449  TEST_FLOATING_EQUALITY( xSlope, 1.32213, 0.01 );
450  TEST_FLOATING_EQUALITY( xErrorNorm[0], 0.116919, 1.0e-4 );
451  TEST_FLOATING_EQUALITY( xDotSlope, 1.32052, 0.01 );
452  TEST_FLOATING_EQUALITY( xDotErrorNorm[0], 0.449888, 1.0e-4 );
453  // At small dt, slopes should be equal to order.
454  //double order = stepper->getOrder();
455  //TEST_FLOATING_EQUALITY( xSlope, order, 0.01 );
456  //TEST_FLOATING_EQUALITY( xDotSlope, order, 0.01 );
457 
458  // Write fine mesh solution at final time
459  // This only works for ONE MPI process
460  if (comm->NumProc() == 1) {
461  RCP<ParameterList> pList =
462  getParametersFromXmlFile("Tempus_BackwardEuler_CDR.xml");
463  RCP<ParameterList> model_pl = sublist(pList, "CDR Model", true);
464  const int num_elements = model_pl->get<int>("num elements");
465  const double left_end = model_pl->get<double>("left end");
466  const double right_end = model_pl->get<double>("right end");
467 
468  const Thyra::VectorBase<double>& x = *(solutions[solutions.size()-1]);
469 
470  std::ofstream ftmp("Tempus_BackwardEuler_CDR-Solution.dat");
471  for (int n = 0; n < num_elements+1; n++) {
472  const double dx = std::fabs(left_end-right_end) /
473  static_cast<double>(num_elements);
474  const double x_coord = left_end + static_cast<double>(n) * dx;
475  ftmp << x_coord << " " << Thyra::get_ele(x,n) << std::endl;
476  }
477  ftmp.close();
478  }
479 
480  Teuchos::TimeMonitor::summarize();
481 }
482 #endif // TEST_CDR
483 
484 
485 #ifdef TEST_VANDERPOL
486 // ************************************************************
487 // ************************************************************
488 TEUCHOS_UNIT_TEST(BackwardEuler, VanDerPol)
489 {
490  RCP<Tempus::IntegratorBasic<double> > integrator;
491  std::vector<RCP<Thyra::VectorBase<double>>> solutions;
492  std::vector<RCP<Thyra::VectorBase<double>>> solutionsDot;
493  std::vector<double> StepSize;
494  std::vector<double> xErrorNorm;
495  std::vector<double> xDotErrorNorm;
496  const int nTimeStepSizes = 4;
497  double dt = 0.05;
498  for (int n=0; n<nTimeStepSizes; n++) {
499 
500  // Read params from .xml file
501  RCP<ParameterList> pList =
502  getParametersFromXmlFile("Tempus_BackwardEuler_VanDerPol.xml");
503 
504  // Setup the VanDerPolModel
505  RCP<ParameterList> vdpm_pl = sublist(pList, "VanDerPolModel", true);
506  auto model = rcp(new VanDerPolModel<double>(vdpm_pl));
507 
508  // Set the step size
509  dt /= 2;
510  if (n == nTimeStepSizes-1) dt /= 10.0;
511 
512  // Setup the Integrator and reset initial time step
513  RCP<ParameterList> pl = sublist(pList, "Tempus", true);
514  pl->sublist("Demo Integrator")
515  .sublist("Time Step Control").set("Initial Time Step", dt);
516  integrator = Tempus::integratorBasic<double>(pl, model);
517 
518  // Integrate to timeMax
519  bool integratorStatus = integrator->advanceTime();
520  TEST_ASSERT(integratorStatus)
521 
522  // Test if at 'Final Time'
523  double time = integrator->getTime();
524  double timeFinal =pl->sublist("Demo Integrator")
525  .sublist("Time Step Control").get<double>("Final Time");
526  double tol = 100.0 * std::numeric_limits<double>::epsilon();
527  TEST_FLOATING_EQUALITY(time, timeFinal, tol);
528 
529  // Store off the final solution and step size
530  StepSize.push_back(dt);
531  auto solution = Thyra::createMember(model->get_x_space());
532  Thyra::copy(*(integrator->getX()),solution.ptr());
533  solutions.push_back(solution);
534  auto solutionDot = Thyra::createMember(model->get_x_space());
535  Thyra::copy(*(integrator->getXdot()),solutionDot.ptr());
536  solutionsDot.push_back(solutionDot);
537 
538  // Output finest temporal solution for plotting
539  // This only works for ONE MPI process
540  if ((n == 0) or (n == nTimeStepSizes-1)) {
541  std::string fname = "Tempus_BackwardEuler_VanDerPol-Ref.dat";
542  if (n == 0) fname = "Tempus_BackwardEuler_VanDerPol.dat";
543  RCP<const SolutionHistory<double> > solutionHistory =
544  integrator->getSolutionHistory();
545  writeSolution(fname, solutionHistory);
546  }
547  }
548 
549  // Check the order and intercept
550  double xSlope = 0.0;
551  double xDotSlope = 0.0;
552  RCP<Tempus::Stepper<double> > stepper = integrator->getStepper();
553  double order = stepper->getOrder();
554  writeOrderError("Tempus_BackwardEuler_VanDerPol-Error.dat",
555  stepper, StepSize,
556  solutions, xErrorNorm, xSlope,
557  solutionsDot, xDotErrorNorm, xDotSlope);
558 
559  TEST_FLOATING_EQUALITY( xSlope, order, 0.10 );
560  TEST_FLOATING_EQUALITY( xErrorNorm[0], 0.571031, 1.0e-4 );
561  TEST_FLOATING_EQUALITY( xDotSlope, 1.74898, 0.10 );
562  TEST_FLOATING_EQUALITY( xDotErrorNorm[0], 1.0038, 1.0e-4 );
563  // At small dt, slopes should be equal to order.
564  //TEST_FLOATING_EQUALITY( xDotSlope, order, 0.01 );
565 
566  Teuchos::TimeMonitor::summarize();
567 }
568 #endif // TEST_VANDERPOL
569 
570 #ifdef TEST_OPT_INTERFACE
571 // ************************************************************
572 // ************************************************************
573 TEUCHOS_UNIT_TEST(BackwardEuler, OptInterface)
574 {
575  // Read params from .xml file
576  RCP<ParameterList> pList =
577  getParametersFromXmlFile("Tempus_BackwardEuler_SinCos.xml");
578 
579  // Setup the SinCosModel
580  RCP<ParameterList> scm_pl = sublist(pList, "SinCosModel", true);
581  auto model = rcp(new SinCosModel<double>(scm_pl));
582 
583  // Setup the Integrator
584  RCP<ParameterList> pl = sublist(pList, "Tempus", true);
585  RCP<Tempus::IntegratorBasic<double> >integrator =
586  Tempus::integratorBasic<double>(pl, model);
587 
588  // Integrate to timeMax
589  bool integratorStatus = integrator->advanceTime();
590  TEST_ASSERT(integratorStatus);
591 
592  // Get solution history
593  RCP<const SolutionHistory<double> > solutionHistory =
594  integrator->getSolutionHistory();
595 
596  // Get the stepper and cast to optimization interface
597  RCP<Tempus::Stepper<double> > stepper = integrator->getStepper();
598  RCP<Tempus::StepperOptimizationInterface<double> > opt_stepper =
599  Teuchos::rcp_dynamic_cast< Tempus::StepperOptimizationInterface<double> >(
600  stepper, true);
601 
602  // Check stencil length
603  TEST_EQUALITY( opt_stepper->stencilLength(), 2);
604 
605  // Create needed vectors/multivectors
606  Teuchos::Array< RCP<const Thyra::VectorBase<double> > > x(2);
607  Teuchos::Array<double> t(2);
608  RCP< const Thyra::VectorBase<double> > p =
609  model->getNominalValues().get_p(0);
610  RCP< Thyra::VectorBase<double> > x_dot =
611  Thyra::createMember(model->get_x_space());
612  RCP< Thyra::VectorBase<double> > f =
613  Thyra::createMember(model->get_f_space());
614  RCP< Thyra::VectorBase<double> > f2 =
615  Thyra::createMember(model->get_f_space());
616  RCP< Thyra::LinearOpBase<double> > dfdx =
617  model->create_W_op();
618  RCP< Thyra::LinearOpBase<double> > dfdx2 =
619  model->create_W_op();
620  RCP< Thyra::MultiVectorBase<double> > dfdx_mv =
621  Teuchos::rcp_dynamic_cast< Thyra::MultiVectorBase<double> >(dfdx,true);
622  RCP< Thyra::MultiVectorBase<double> > dfdx_mv2 =
623  Teuchos::rcp_dynamic_cast< Thyra::MultiVectorBase<double> >(dfdx2,true);
624  const int num_p = p->range()->dim();
625  RCP< Thyra::MultiVectorBase<double> > dfdp =
626  Thyra::createMembers(model->get_f_space(), num_p);
627  RCP< Thyra::MultiVectorBase<double> > dfdp2 =
628  Thyra::createMembers(model->get_f_space(), num_p);
629  RCP< Thyra::LinearOpWithSolveBase<double> > W =
630  model->create_W();
631  RCP< Thyra::LinearOpWithSolveBase<double> > W2 =
632  model->create_W();
633  RCP< Thyra::MultiVectorBase<double> > tmp =
634  Thyra::createMembers(model->get_x_space(), num_p);
635  RCP< Thyra::MultiVectorBase<double> > tmp2 =
636  Thyra::createMembers(model->get_x_space(), num_p);
637  std::vector<double> nrms(num_p);
638  double err;
639 
640  // Loop over states, checking residuals and derivatives
641  const int n = solutionHistory->getNumStates();
642  for (int i=1; i<n; ++i) {
643  RCP<const SolutionState<double> > state = (*solutionHistory)[i];
644  RCP<const SolutionState<double> > prev_state = (*solutionHistory)[i-1];
645 
646  // Fill x, t stencils
647  x[0] = state->getX();
648  x[1] = prev_state->getX();
649  t[0] = state->getTime();
650  t[1] = prev_state->getTime();
651 
652  // Compute x_dot
653  const double dt = t[0]-t[1];
654  Thyra::V_StVpStV(x_dot.ptr(), 1.0/dt, *(x[0]), -1.0/dt, *(x[1]));
655 
656  // Create model inargs
657  typedef Thyra::ModelEvaluatorBase MEB;
658  MEB::InArgs<double> in_args = model->createInArgs();
659  MEB::OutArgs<double> out_args = model->createOutArgs();
660  in_args.set_x(x[0]);
661  in_args.set_x_dot(x_dot);
662  in_args.set_t(t[0]);
663  in_args.set_p(0,p);
664 
665  const double tol = 1.0e-14;
666 
667  // Check residual
668  opt_stepper->computeStepResidual(*f, x, t, *p, 0);
669  out_args.set_f(f2);
670  model->evalModel(in_args, out_args);
671  out_args.set_f(Teuchos::null);
672  Thyra::V_VmV(f.ptr(), *f, *f2);
673  err = Thyra::norm(*f);
674  TEST_FLOATING_EQUALITY(err, 0.0, tol);
675 
676  // Check df/dx_n
677  // df/dx_n = df/dx_dot * dx_dot/dx_n + df/dx_n = 1/dt*df/dx_dot + df/dx_n
678  opt_stepper->computeStepJacobian(*dfdx, x, t, *p, 0, 0);
679  out_args.set_W_op(dfdx2);
680  in_args.set_alpha(1.0/dt);
681  in_args.set_beta(1.0);
682  model->evalModel(in_args, out_args);
683  out_args.set_W_op(Teuchos::null);
684  Thyra::V_VmV(dfdx_mv.ptr(), *dfdx_mv, *dfdx_mv2);
685  Thyra::norms(*dfdx_mv, Teuchos::arrayViewFromVector(nrms));
686  err = 0.0;
687  for (auto nrm : nrms) err += nrm;
688  TEST_FLOATING_EQUALITY(err, 0.0, tol);
689 
690  // Check df/dx_{n-1}
691  // df/dx_{n-1} = df/dx_dot * dx_dot/dx_{n-1} = -1/dt*df/dx_dot
692  opt_stepper->computeStepJacobian(*dfdx, x, t, *p, 0, 1);
693  out_args.set_W_op(dfdx2);
694  in_args.set_alpha(-1.0/dt);
695  in_args.set_beta(0.0);
696  model->evalModel(in_args, out_args);
697  out_args.set_W_op(Teuchos::null);
698  Thyra::V_VmV(dfdx_mv.ptr(), *dfdx_mv, *dfdx_mv2);
699  Thyra::norms(*dfdx_mv, Teuchos::arrayViewFromVector(nrms));
700  err = 0.0;
701  for (auto nrm : nrms) err += nrm;
702  TEST_FLOATING_EQUALITY(err, 0.0, tol);
703 
704  // Check df/dp
705  opt_stepper->computeStepParamDeriv(*dfdp, x, t, *p, 0);
706  out_args.set_DfDp(
707  0, MEB::Derivative<double>(dfdp2, MEB::DERIV_MV_JACOBIAN_FORM));
708  model->evalModel(in_args, out_args);
709  out_args.set_DfDp(0, MEB::Derivative<double>());
710  Thyra::V_VmV(dfdp.ptr(), *dfdp, *dfdp2);
711  Thyra::norms(*dfdp, Teuchos::arrayViewFromVector(nrms));
712  err = 0.0;
713  for (auto nrm : nrms) err += nrm;
714  TEST_FLOATING_EQUALITY(err, 0.0, tol);
715 
716  // Check W
717  opt_stepper->computeStepSolver(*W, x, t, *p, 0);
718  out_args.set_W(W2);
719  in_args.set_alpha(1.0/dt);
720  in_args.set_beta(1.0);
721  model->evalModel(in_args, out_args);
722  out_args.set_W(Teuchos::null);
723  // note: dfdp overwritten above so dfdp != dfdp2
724  Thyra::solve(*W, Thyra::NOTRANS, *dfdp2, tmp.ptr());
725  Thyra::solve(*W2, Thyra::NOTRANS, *dfdp2, tmp2.ptr());
726  Thyra::V_VmV(tmp.ptr(), *tmp, *tmp2);
727  Thyra::norms(*tmp, Teuchos::arrayViewFromVector(nrms));
728  err = 0.0;
729  for (auto nrm : nrms) err += nrm;
730  TEST_FLOATING_EQUALITY(err, 0.0, tol);
731  }
732 
733  Teuchos::TimeMonitor::summarize();
734 }
735 #endif // TEST_OPT_INTERFACE
736 
737 
738 } // namespace Tempus_Test
Sine-Cosine model problem from Rythmos. This is a canonical Sine-Cosine differential equation with a...
void writeSolution(const std::string filename, Teuchos::RCP< const Tempus::SolutionHistory< Scalar > > solutionHistory)
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)
TEUCHOS_UNIT_TEST(BackwardEuler, SinCos_ASA)
TimeStepControl manages the time step size. There several mechanicisms that effect the time step size...
Teuchos::RCP< SolutionHistory< Scalar > > solutionHistory(Teuchos::RCP< Teuchos::ParameterList > pList=Teuchos::null)
Nonmember constructor.
SolutionHistory is basically a container of SolutionStates. SolutionHistory maintains a collection of...
Keep a fix number of states.
van der Pol model problem for nonlinear electrical circuit.
Stepper interface to support full-space optimization.
Solution state for integrators and steppers. SolutionState contains the metadata for solutions and th...
1D CGFEM model for convection/diffusion/reaction