Tempus  Version of the Day
Time Integration
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::ParameterList;
41 using Teuchos::sublist;
42 using Teuchos::getParametersFromXmlFile;
43 
44 using Tempus::IntegratorBasic;
45 using Tempus::SolutionHistory;
46 using Tempus::SolutionState;
47 
48 // Comment out any of the following tests to exclude from build/run.
49 #define TEST_PARAMETERLIST
50 #define TEST_CONSTRUCTING_FROM_DEFAULTS
51 #define TEST_SINCOS
52 #define TEST_CDR
53 #define TEST_VANDERPOL
54 #define TEST_OPT_INTERFACE
55 
56 
57 #ifdef TEST_PARAMETERLIST
58 // ************************************************************
59 // ************************************************************
60 TEUCHOS_UNIT_TEST(BackwardEuler, ParameterList)
61 {
62  // Read params from .xml file
63  RCP<ParameterList> pList =
64  getParametersFromXmlFile("Tempus_BackwardEuler_SinCos.xml");
65 
66  // Setup the SinCosModel
67  RCP<ParameterList> scm_pl = sublist(pList, "SinCosModel", true);
68  RCP<SinCosModel<double> > model =
69  Teuchos::rcp(new SinCosModel<double> (scm_pl));
70 
71  RCP<ParameterList> tempusPL = sublist(pList, "Tempus", true);
72 
73  // Test constructor IntegratorBasic(tempusPL, model)
74  {
75  RCP<Tempus::IntegratorBasic<double> > integrator =
76  Tempus::integratorBasic<double>(tempusPL, model);
77 
78  RCP<ParameterList> stepperPL = sublist(tempusPL, "Default Stepper", true);
79  // Remove Predictor for comparison
80  stepperPL->remove("Predictor Name");
81  stepperPL->remove("Default Predictor");
82  RCP<ParameterList> defaultPL =
83  integrator->getStepper()->getDefaultParameters();
84  TEST_ASSERT(haveSameValues(*stepperPL, *defaultPL, true))
85  }
86 
87  // Test constructor IntegratorBasic(model, stepperType)
88  {
89  RCP<Tempus::IntegratorBasic<double> > integrator =
90  Tempus::integratorBasic<double>(model, "Backward Euler");
91 
92  RCP<ParameterList> stepperPL = sublist(tempusPL, "Default Stepper", true);
93  RCP<ParameterList> defaultPL =
94  integrator->getStepper()->getDefaultParameters();
95 
96  TEST_ASSERT(haveSameValues(*stepperPL, *defaultPL, true))
97  }
98 }
99 #endif // TEST_PARAMETERLIST
100 
101 
102 #ifdef TEST_CONSTRUCTING_FROM_DEFAULTS
103 // ************************************************************
104 // ************************************************************
105 TEUCHOS_UNIT_TEST(BackwardEuler, ConstructingFromDefaults)
106 {
107  double dt = 0.1;
108 
109  // Read params from .xml file
110  RCP<ParameterList> pList =
111  getParametersFromXmlFile("Tempus_BackwardEuler_SinCos.xml");
112  RCP<ParameterList> pl = sublist(pList, "Tempus", true);
113 
114  // Setup the SinCosModel
115  RCP<ParameterList> scm_pl = sublist(pList, "SinCosModel", true);
116  //RCP<SinCosModel<double> > model = sineCosineModel(scm_pl);
117  RCP<SinCosModel<double> > model =
118  Teuchos::rcp(new SinCosModel<double>(scm_pl));
119 
120  // Setup Stepper for field solve ----------------------------
121  RCP<Tempus::StepperBackwardEuler<double> > stepper =
122  Teuchos::rcp(new Tempus::StepperBackwardEuler<double>(model));
123  //{
124  // // Setup a linear NOX solve
125  // RCP<ParameterList> sPL = stepper->getNonconstParameterList();
126  // std::string solverName = sPL->get<std::string>("Solver Name");
127  // RCP<ParameterList> solverPL = Teuchos::sublist(sPL, solverName, true);
128  // stepper->setSolver(solverPL);
129  // stepper->initialize();
130  //}
131 
132  // Setup TimeStepControl ------------------------------------
133  RCP<Tempus::TimeStepControl<double> > timeStepControl =
134  Teuchos::rcp(new Tempus::TimeStepControl<double>());
135  ParameterList tscPL = pl->sublist("Default Integrator")
136  .sublist("Time Step Control");
137  timeStepControl->setStepType (tscPL.get<std::string>("Integrator Step Type"));
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  Thyra::ModelEvaluatorBase::InArgs<double> inArgsIC =
146  stepper->getModel()->getNominalValues();
147  RCP<Thyra::VectorBase<double> > icSolution =
148  Teuchos::rcp_const_cast<Thyra::VectorBase<double> > (inArgsIC.get_x());
149  RCP<Tempus::SolutionState<double> > icState =
150  Teuchos::rcp(new Tempus::SolutionState<double>(icSolution));
151  icState->setTime (timeStepControl->getInitTime());
152  icState->setIndex (timeStepControl->getInitIndex());
153  icState->setTimeStep(0.0);
154  icState->setOrder (stepper->getOrder());
155  icState->setSolutionStatus(Tempus::Status::PASSED); // ICs are passing.
156 
157  // Setup SolutionHistory ------------------------------------
158  RCP<Tempus::SolutionHistory<double> > solutionHistory =
159  Teuchos::rcp(new Tempus::SolutionHistory<double>());
160  solutionHistory->setName("Forward States");
161  solutionHistory->setStorageType(Tempus::STORAGE_TYPE_STATIC);
162  solutionHistory->setStorageLimit(2);
163  solutionHistory->addState(icState);
164 
165  // Setup Integrator -----------------------------------------
166  RCP<Tempus::IntegratorBasic<double> > integrator =
167  Tempus::integratorBasic<double>();
168  integrator->setStepperWStepper(stepper);
169  integrator->setTimeStepControl(timeStepControl);
170  integrator->setSolutionHistory(solutionHistory);
171  //integrator->setObserver(...);
172  integrator->initialize();
173 
174 
175  // Integrate to timeMax
176  bool integratorStatus = integrator->advanceTime();
177  TEST_ASSERT(integratorStatus)
178 
179 
180  // Test if at 'Final Time'
181  double time = integrator->getTime();
182  double timeFinal =pl->sublist("Default Integrator")
183  .sublist("Time Step Control").get<double>("Final Time");
184  TEST_FLOATING_EQUALITY(time, timeFinal, 1.0e-14);
185 
186  // Time-integrated solution and the exact solution
187  RCP<Thyra::VectorBase<double> > x = integrator->getX();
188  RCP<const Thyra::VectorBase<double> > x_exact =
189  model->getExactSolution(time).get_x();
190 
191  // Calculate the error
192  RCP<Thyra::VectorBase<double> > xdiff = x->clone_v();
193  Thyra::V_StVpStV(xdiff.ptr(), 1.0, *x_exact, -1.0, *(x));
194 
195  // Check the order and intercept
196  std::cout << " Stepper = BackwardEuler" << std::endl;
197  std::cout << " =========================" << std::endl;
198  std::cout << " Exact solution : " << get_ele(*(x_exact), 0) << " "
199  << get_ele(*(x_exact), 1) << std::endl;
200  std::cout << " Computed solution: " << get_ele(*(x ), 0) << " "
201  << get_ele(*(x ), 1) << std::endl;
202  std::cout << " Difference : " << get_ele(*(xdiff ), 0) << " "
203  << get_ele(*(xdiff ), 1) << std::endl;
204  std::cout << " =========================" << std::endl;
205  TEST_FLOATING_EQUALITY(get_ele(*(x), 0), 0.798923, 1.0e-4 );
206  TEST_FLOATING_EQUALITY(get_ele(*(x), 1), 0.516729, 1.0e-4 );
207 }
208 #endif // TEST_CONSTRUCTING_FROM_DEFAULTS
209 
210 
211 #ifdef TEST_SINCOS
212 // ************************************************************
213 // ************************************************************
214 TEUCHOS_UNIT_TEST(BackwardEuler, SinCos)
215 {
216  RCP<Tempus::IntegratorBasic<double> > integrator;
217  std::vector<RCP<Thyra::VectorBase<double>>> solutions;
218  std::vector<RCP<Thyra::VectorBase<double>>> solutionsDot;
219  std::vector<double> StepSize;
220  std::vector<double> xErrorNorm;
221  std::vector<double> xDotErrorNorm;
222  const int nTimeStepSizes = 7;
223  double dt = 0.2;
224  double time = 0.0;
225  for (int n=0; n<nTimeStepSizes; n++) {
226 
227  // Read params from .xml file
228  RCP<ParameterList> pList =
229  getParametersFromXmlFile("Tempus_BackwardEuler_SinCos.xml");
230 
231  //std::ofstream ftmp("PL.txt");
232  //pList->print(ftmp);
233  //ftmp.close();
234 
235  // Setup the SinCosModel
236  RCP<ParameterList> scm_pl = sublist(pList, "SinCosModel", true);
237  //RCP<SinCosModel<double> > model = sineCosineModel(scm_pl);
238  RCP<SinCosModel<double> > model =
239  Teuchos::rcp(new SinCosModel<double>(scm_pl));
240 
241  dt /= 2;
242 
243  // Setup the Integrator and reset initial time step
244  RCP<ParameterList> pl = sublist(pList, "Tempus", true);
245  pl->sublist("Default Integrator")
246  .sublist("Time Step Control").set("Initial Time Step", dt);
247  integrator = Tempus::integratorBasic<double>(pl, model);
248 
249  // Initial Conditions
250  // During the Integrator construction, the initial SolutionState
251  // is set by default to model->getNominalVales().get_x(). However,
252  // the application can set it also by integrator->setInitialState.
253  RCP<Thyra::VectorBase<double> > x0 =
254  model->getNominalValues().get_x()->clone_v();
255  integrator->setInitialState(0.0, x0);
256 
257  // Integrate to timeMax
258  bool integratorStatus = integrator->advanceTime();
259  TEST_ASSERT(integratorStatus)
260 
261  // Test if at 'Final Time'
262  time = integrator->getTime();
263  double timeFinal =pl->sublist("Default Integrator")
264  .sublist("Time Step Control").get<double>("Final Time");
265  TEST_FLOATING_EQUALITY(time, timeFinal, 1.0e-14);
266 
267  // Plot sample solution and exact solution
268  if (n == 0) {
269  RCP<const SolutionHistory<double> > solutionHistory =
270  integrator->getSolutionHistory();
271  writeSolution("Tempus_BackwardEuler_SinCos.dat", solutionHistory);
272 
273  RCP<Tempus::SolutionHistory<double> > solnHistExact =
274  Teuchos::rcp(new Tempus::SolutionHistory<double>());
275  for (int i=0; i<solutionHistory->getNumStates(); i++) {
276  double time = (*solutionHistory)[i]->getTime();
277  RCP<Tempus::SolutionState<double> > state =
278  Teuchos::rcp(new Tempus::SolutionState<double>(
279  model->getExactSolution(time).get_x(),
280  model->getExactSolution(time).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 solution = Thyra::createMember(model->get_x_space());
298  Thyra::copy(*(model->getExactSolution(time).get_x()),solution.ptr());
299  solutions.push_back(solution);
300  auto solutionDot = Thyra::createMember(model->get_x_space());
301  Thyra::copy(*(model->getExactSolution(time).get_x_dot()),
302  solutionDot.ptr());
303  solutionsDot.push_back(solutionDot);
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 = Teuchos::rcp(new Epetra_MpiComm(MPI_COMM_WORLD));
336 #else
337  comm = Teuchos::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  RCP<Tempus_Test::CDR_Model<double>> model =
363  Teuchos::rcp(new Tempus_Test::CDR_Model<double>(comm,
364  num_elements,
365  left_end,
366  right_end,
367  a_convection,
368  k_source));
369 
370  // Set the factory
371  ::Stratimikos::DefaultLinearSolverBuilder builder;
372 
373  Teuchos::RCP<Teuchos::ParameterList> p =
374  Teuchos::rcp(new Teuchos::ParameterList);
375  p->set("Linear Solver Type", "Belos");
376  p->set("Preconditioner Type", "None");
377  builder.setParameterList(p);
378 
379  Teuchos::RCP< ::Thyra::LinearOpWithSolveFactoryBase<double> >
380  lowsFactory = builder.createLinearSolveStrategy("");
381 
382  model->set_W_factory(lowsFactory);
383 
384  // Set the step size
385  dt /= 2;
386 
387  // Setup the Integrator and reset initial time step
388  RCP<ParameterList> pl = sublist(pList, "Tempus", true);
389  pl->sublist("Demo Integrator")
390  .sublist("Time Step Control").set("Initial Time Step", dt);
391  integrator = Tempus::integratorBasic<double>(pl, model);
392 
393  // Integrate to timeMax
394  bool integratorStatus = integrator->advanceTime();
395  TEST_ASSERT(integratorStatus)
396 
397  // Test if at 'Final Time'
398  double time = integrator->getTime();
399  double timeFinal =pl->sublist("Demo Integrator")
400  .sublist("Time Step Control").get<double>("Final Time");
401  double tol = 100.0 * std::numeric_limits<double>::epsilon();
402  TEST_FLOATING_EQUALITY(time, timeFinal, tol);
403 
404  // Store off the final solution and step size
405  StepSize.push_back(dt);
406  auto solution = Thyra::createMember(model->get_x_space());
407  Thyra::copy(*(integrator->getX()),solution.ptr());
408  solutions.push_back(solution);
409  auto solutionDot = Thyra::createMember(model->get_x_space());
410  Thyra::copy(*(integrator->getXdot()),solutionDot.ptr());
411  solutionsDot.push_back(solutionDot);
412 
413  // Output finest temporal solution for plotting
414  // This only works for ONE MPI process
415  if ((n == nTimeStepSizes-1) && (comm->NumProc() == 1)) {
416  std::ofstream ftmp("Tempus_BackwardEuler_CDR.dat");
417  ftmp << "TITLE=\"Backward Euler Solution to CDR\"\n"
418  << "VARIABLES=\"z\",\"T\"\n";
419  const double dx = std::fabs(left_end-right_end) /
420  static_cast<double>(num_elements);
421  RCP<const SolutionHistory<double> > solutionHistory =
422  integrator->getSolutionHistory();
423  int nStates = solutionHistory->getNumStates();
424  for (int i=0; i<nStates; i++) {
425  RCP<const SolutionState<double> > solutionState = (*solutionHistory)[i];
426  RCP<const Thyra::VectorBase<double> > x = solutionState->getX();
427  double ttime = solutionState->getTime();
428  ftmp << "ZONE T=\"Time="<<ttime<<"\", I="
429  <<num_elements+1<<", F=BLOCK\n";
430  for (int j = 0; j < num_elements+1; j++) {
431  const double x_coord = left_end + static_cast<double>(j) * dx;
432  ftmp << x_coord << " ";
433  }
434  ftmp << std::endl;
435  for (int j=0; j<num_elements+1; j++) ftmp << get_ele(*x, j) << " ";
436  ftmp << std::endl;
437  }
438  ftmp.close();
439  }
440  }
441 
442  // Check the order and intercept
443  double xSlope = 0.0;
444  double xDotSlope = 0.0;
445  RCP<Tempus::Stepper<double> > stepper = integrator->getStepper();
446  writeOrderError("Tempus_BackwardEuler_CDR-Error.dat",
447  stepper, StepSize,
448  solutions, xErrorNorm, xSlope,
449  solutionsDot, xDotErrorNorm, xDotSlope);
450 
451  TEST_FLOATING_EQUALITY( xSlope, 1.32213, 0.01 );
452  TEST_FLOATING_EQUALITY( xErrorNorm[0], 0.116919, 1.0e-4 );
453  TEST_FLOATING_EQUALITY( xDotSlope, 1.32052, 0.01 );
454  TEST_FLOATING_EQUALITY( xDotErrorNorm[0], 0.449888, 1.0e-4 );
455  // At small dt, slopes should be equal to order.
456  //double order = stepper->getOrder();
457  //TEST_FLOATING_EQUALITY( xSlope, order, 0.01 );
458  //TEST_FLOATING_EQUALITY( xDotSlope, order, 0.01 );
459 
460  // Write fine mesh solution at final time
461  // This only works for ONE MPI process
462  if (comm->NumProc() == 1) {
463  RCP<ParameterList> pList =
464  getParametersFromXmlFile("Tempus_BackwardEuler_CDR.xml");
465  RCP<ParameterList> model_pl = sublist(pList, "CDR Model", true);
466  const int num_elements = model_pl->get<int>("num elements");
467  const double left_end = model_pl->get<double>("left end");
468  const double right_end = model_pl->get<double>("right end");
469 
470  const Thyra::VectorBase<double>& x = *(solutions[solutions.size()-1]);
471 
472  std::ofstream ftmp("Tempus_BackwardEuler_CDR-Solution.dat");
473  for (int n = 0; n < num_elements+1; n++) {
474  const double dx = std::fabs(left_end-right_end) /
475  static_cast<double>(num_elements);
476  const double x_coord = left_end + static_cast<double>(n) * dx;
477  ftmp << x_coord << " " << Thyra::get_ele(x,n) << std::endl;
478  }
479  ftmp.close();
480  }
481 
482  Teuchos::TimeMonitor::summarize();
483 }
484 #endif // TEST_CDR
485 
486 
487 #ifdef TEST_VANDERPOL
488 // ************************************************************
489 // ************************************************************
490 TEUCHOS_UNIT_TEST(BackwardEuler, VanDerPol)
491 {
492  RCP<Tempus::IntegratorBasic<double> > integrator;
493  std::vector<RCP<Thyra::VectorBase<double>>> solutions;
494  std::vector<RCP<Thyra::VectorBase<double>>> solutionsDot;
495  std::vector<double> StepSize;
496  std::vector<double> xErrorNorm;
497  std::vector<double> xDotErrorNorm;
498  const int nTimeStepSizes = 4;
499  double dt = 0.05;
500  for (int n=0; n<nTimeStepSizes; n++) {
501 
502  // Read params from .xml file
503  RCP<ParameterList> pList =
504  getParametersFromXmlFile("Tempus_BackwardEuler_VanDerPol.xml");
505 
506  // Setup the VanDerPolModel
507  RCP<ParameterList> vdpm_pl = sublist(pList, "VanDerPolModel", true);
508  RCP<VanDerPolModel<double> > model =
509  Teuchos::rcp(new VanDerPolModel<double>(vdpm_pl));
510 
511  // Set the step size
512  dt /= 2;
513  if (n == nTimeStepSizes-1) dt /= 10.0;
514 
515  // Setup the Integrator and reset initial time step
516  RCP<ParameterList> pl = sublist(pList, "Tempus", true);
517  pl->sublist("Demo Integrator")
518  .sublist("Time Step Control").set("Initial Time Step", dt);
519  integrator = Tempus::integratorBasic<double>(pl, model);
520 
521  // Integrate to timeMax
522  bool integratorStatus = integrator->advanceTime();
523  TEST_ASSERT(integratorStatus)
524 
525  // Test if at 'Final Time'
526  double time = integrator->getTime();
527  double timeFinal =pl->sublist("Demo Integrator")
528  .sublist("Time Step Control").get<double>("Final Time");
529  double tol = 100.0 * std::numeric_limits<double>::epsilon();
530  TEST_FLOATING_EQUALITY(time, timeFinal, tol);
531 
532  // Store off the final solution and step size
533  StepSize.push_back(dt);
534  auto solution = Thyra::createMember(model->get_x_space());
535  Thyra::copy(*(integrator->getX()),solution.ptr());
536  solutions.push_back(solution);
537  auto solutionDot = Thyra::createMember(model->get_x_space());
538  Thyra::copy(*(integrator->getXdot()),solutionDot.ptr());
539  solutionsDot.push_back(solutionDot);
540 
541  // Output finest temporal solution for plotting
542  // This only works for ONE MPI process
543  if ((n == 0) or (n == nTimeStepSizes-1)) {
544  std::string fname = "Tempus_BackwardEuler_VanDerPol-Ref.dat";
545  if (n == 0) fname = "Tempus_BackwardEuler_VanDerPol.dat";
546  RCP<const SolutionHistory<double> > solutionHistory =
547  integrator->getSolutionHistory();
548  writeSolution(fname, solutionHistory);
549  }
550  }
551 
552  // Check the order and intercept
553  double xSlope = 0.0;
554  double xDotSlope = 0.0;
555  RCP<Tempus::Stepper<double> > stepper = integrator->getStepper();
556  double order = stepper->getOrder();
557  writeOrderError("Tempus_BackwardEuler_VanDerPol-Error.dat",
558  stepper, StepSize,
559  solutions, xErrorNorm, xSlope,
560  solutionsDot, xDotErrorNorm, xDotSlope);
561 
562  TEST_FLOATING_EQUALITY( xSlope, order, 0.10 );
563  TEST_FLOATING_EQUALITY( xErrorNorm[0], 0.571031, 1.0e-4 );
564  TEST_FLOATING_EQUALITY( xDotSlope, 1.74898, 0.10 );
565  TEST_FLOATING_EQUALITY( xDotErrorNorm[0], 1.0038, 1.0e-4 );
566  // At small dt, slopes should be equal to order.
567  //TEST_FLOATING_EQUALITY( xDotSlope, order, 0.01 );
568 
569  Teuchos::TimeMonitor::summarize();
570 }
571 #endif // TEST_VANDERPOL
572 
573 #ifdef TEST_OPT_INTERFACE
574 // ************************************************************
575 // ************************************************************
576 TEUCHOS_UNIT_TEST(BackwardEuler, OptInterface)
577 {
578  // Read params from .xml file
579  RCP<ParameterList> pList =
580  getParametersFromXmlFile("Tempus_BackwardEuler_SinCos.xml");
581 
582  // Setup the SinCosModel
583  RCP<ParameterList> scm_pl = sublist(pList, "SinCosModel", true);
584  RCP<SinCosModel<double> > model =
585  Teuchos::rcp(new SinCosModel<double>(scm_pl));
586 
587  // Setup the Integrator
588  RCP<ParameterList> pl = sublist(pList, "Tempus", true);
589  RCP<Tempus::IntegratorBasic<double> >integrator =
590  Tempus::integratorBasic<double>(pl, model);
591 
592  // Integrate to timeMax
593  bool integratorStatus = integrator->advanceTime();
594  TEST_ASSERT(integratorStatus);
595 
596  // Get solution history
597  RCP<const SolutionHistory<double> > solutionHistory =
598  integrator->getSolutionHistory();
599 
600  // Get the stepper and cast to optimization interface
601  RCP<Tempus::Stepper<double> > stepper = integrator->getStepper();
602  RCP<Tempus::StepperOptimizationInterface<double> > opt_stepper =
603  Teuchos::rcp_dynamic_cast< Tempus::StepperOptimizationInterface<double> >(
604  stepper, true);
605 
606  // Check stencil length
607  TEST_EQUALITY( opt_stepper->stencilLength(), 2);
608 
609  // Create needed vectors/multivectors
610  Teuchos::Array< RCP<const Thyra::VectorBase<double> > > x(2);
611  Teuchos::Array<double> t(2);
612  RCP< const Thyra::VectorBase<double> > p =
613  model->getNominalValues().get_p(0);
614  RCP< Thyra::VectorBase<double> > x_dot =
615  Thyra::createMember(model->get_x_space());
616  RCP< Thyra::VectorBase<double> > f =
617  Thyra::createMember(model->get_f_space());
618  RCP< Thyra::VectorBase<double> > f2 =
619  Thyra::createMember(model->get_f_space());
620  RCP< Thyra::LinearOpBase<double> > dfdx =
621  model->create_W_op();
622  RCP< Thyra::LinearOpBase<double> > dfdx2 =
623  model->create_W_op();
624  RCP< Thyra::MultiVectorBase<double> > dfdx_mv =
625  Teuchos::rcp_dynamic_cast< Thyra::MultiVectorBase<double> >(dfdx,true);
626  RCP< Thyra::MultiVectorBase<double> > dfdx_mv2 =
627  Teuchos::rcp_dynamic_cast< Thyra::MultiVectorBase<double> >(dfdx2,true);
628  const int num_p = p->range()->dim();
629  RCP< Thyra::MultiVectorBase<double> > dfdp =
630  Thyra::createMembers(model->get_f_space(), num_p);
631  RCP< Thyra::MultiVectorBase<double> > dfdp2 =
632  Thyra::createMembers(model->get_f_space(), num_p);
633  RCP< Thyra::LinearOpWithSolveBase<double> > W =
634  model->create_W();
635  RCP< Thyra::LinearOpWithSolveBase<double> > W2 =
636  model->create_W();
637  RCP< Thyra::MultiVectorBase<double> > tmp =
638  Thyra::createMembers(model->get_x_space(), num_p);
639  RCP< Thyra::MultiVectorBase<double> > tmp2 =
640  Thyra::createMembers(model->get_x_space(), num_p);
641  std::vector<double> nrms(num_p);
642  double err;
643 
644  // Loop over states, checking residuals and derivatives
645  const int n = solutionHistory->getNumStates();
646  for (int i=1; i<n; ++i) {
647  RCP<const SolutionState<double> > state = (*solutionHistory)[i];
648  RCP<const SolutionState<double> > prev_state = (*solutionHistory)[i-1];
649 
650  // Fill x, t stencils
651  x[0] = state->getX();
652  x[1] = prev_state->getX();
653  t[0] = state->getTime();
654  t[1] = prev_state->getTime();
655 
656  // Compute x_dot
657  const double dt = t[0]-t[1];
658  Thyra::V_StVpStV(x_dot.ptr(), 1.0/dt, *(x[0]), -1.0/dt, *(x[1]));
659 
660  // Create model inargs
661  typedef Thyra::ModelEvaluatorBase MEB;
662  MEB::InArgs<double> in_args = model->createInArgs();
663  MEB::OutArgs<double> out_args = model->createOutArgs();
664  in_args.set_x(x[0]);
665  in_args.set_x_dot(x_dot);
666  in_args.set_t(t[0]);
667  in_args.set_p(0,p);
668 
669  const double tol = 1.0e-14;
670 
671  // Check residual
672  opt_stepper->computeStepResidual(*f, x, t, *p, 0);
673  out_args.set_f(f2);
674  model->evalModel(in_args, out_args);
675  out_args.set_f(Teuchos::null);
676  Thyra::V_VmV(f.ptr(), *f, *f2);
677  err = Thyra::norm(*f);
678  TEST_FLOATING_EQUALITY(err, 0.0, tol);
679 
680  // Check df/dx_n
681  // df/dx_n = df/dx_dot * dx_dot/dx_n + df/dx_n = 1/dt*df/dx_dot + df/dx_n
682  opt_stepper->computeStepJacobian(*dfdx, x, t, *p, 0, 0);
683  out_args.set_W_op(dfdx2);
684  in_args.set_alpha(1.0/dt);
685  in_args.set_beta(1.0);
686  model->evalModel(in_args, out_args);
687  out_args.set_W_op(Teuchos::null);
688  Thyra::V_VmV(dfdx_mv.ptr(), *dfdx_mv, *dfdx_mv2);
689  Thyra::norms(*dfdx_mv, Teuchos::arrayViewFromVector(nrms));
690  err = 0.0;
691  for (auto nrm : nrms) err += nrm;
692  TEST_FLOATING_EQUALITY(err, 0.0, tol);
693 
694  // Check df/dx_{n-1}
695  // df/dx_{n-1} = df/dx_dot * dx_dot/dx_{n-1} = -1/dt*df/dx_dot
696  opt_stepper->computeStepJacobian(*dfdx, x, t, *p, 0, 1);
697  out_args.set_W_op(dfdx2);
698  in_args.set_alpha(-1.0/dt);
699  in_args.set_beta(0.0);
700  model->evalModel(in_args, out_args);
701  out_args.set_W_op(Teuchos::null);
702  Thyra::V_VmV(dfdx_mv.ptr(), *dfdx_mv, *dfdx_mv2);
703  Thyra::norms(*dfdx_mv, Teuchos::arrayViewFromVector(nrms));
704  err = 0.0;
705  for (auto nrm : nrms) err += nrm;
706  TEST_FLOATING_EQUALITY(err, 0.0, tol);
707 
708  // Check df/dp
709  opt_stepper->computeStepParamDeriv(*dfdp, x, t, *p, 0);
710  out_args.set_DfDp(
711  0, MEB::Derivative<double>(dfdp2, MEB::DERIV_MV_JACOBIAN_FORM));
712  model->evalModel(in_args, out_args);
713  out_args.set_DfDp(0, MEB::Derivative<double>());
714  Thyra::V_VmV(dfdp.ptr(), *dfdp, *dfdp2);
715  Thyra::norms(*dfdp, Teuchos::arrayViewFromVector(nrms));
716  err = 0.0;
717  for (auto nrm : nrms) err += nrm;
718  TEST_FLOATING_EQUALITY(err, 0.0, tol);
719 
720  // Check W
721  opt_stepper->computeStepSolver(*W, x, t, *p, 0);
722  out_args.set_W(W2);
723  in_args.set_alpha(1.0/dt);
724  in_args.set_beta(1.0);
725  model->evalModel(in_args, out_args);
726  out_args.set_W(Teuchos::null);
727  // note: dfdp overwritten above so dfdp != dfdp2
728  Thyra::solve(*W, Thyra::NOTRANS, *dfdp2, tmp.ptr());
729  Thyra::solve(*W2, Thyra::NOTRANS, *dfdp2, tmp2.ptr());
730  Thyra::V_VmV(tmp.ptr(), *tmp, *tmp2);
731  Thyra::norms(*tmp, Teuchos::arrayViewFromVector(nrms));
732  err = 0.0;
733  for (auto nrm : nrms) err += nrm;
734  TEST_FLOATING_EQUALITY(err, 0.0, tol);
735  }
736 
737  Teuchos::TimeMonitor::summarize();
738 }
739 #endif // TEST_OPT_INTERFACE
740 
741 
742 } // 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::StepperOptimizationInterface
Stepper interface to support full-space optimization.
Definition: Tempus_StepperOptimizationInterface.hpp:39
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::StepperBackwardEuler
Backward Euler time stepper.
Definition: Tempus_StepperBackwardEuler_decl.hpp:33
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