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