Tempus Version of the Day
Time Integration
Loading...
Searching...
No Matches
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"
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
37namespace Tempus_Test {
38
39using Teuchos::RCP;
40using Teuchos::rcp;
41using Teuchos::rcp_const_cast;
42using Teuchos::ParameterList;
43using Teuchos::sublist;
44using Teuchos::getParametersFromXmlFile;
45
49
50
51// ************************************************************
52// ************************************************************
53TEUCHOS_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 out << std::endl;
77 out << "stepperPL -------------- \n" << *stepperPL << std::endl;
78 out << "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, std::string("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 out << std::endl;
97 out << "stepperPL -------------- \n" << *stepperPL << std::endl;
98 out << "defaultPL -------------- \n" << *defaultPL << std::endl;
99 }
100 TEST_ASSERT(pass)
101 }
102}
103
104
105// ************************************************************
106// ************************************************************
107TEUCHOS_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 out << " Stepper = " << stepper->description()
198 << " with " << option << std::endl;
199 out << " =========================" << std::endl;
200 out << " Exact solution : " << get_ele(*(x_exact), 0) << " "
201 << get_ele(*(x_exact), 1) << std::endl;
202 out << " Computed solution: " << get_ele(*(x ), 0) << " "
203 << get_ele(*(x ), 1) << std::endl;
204 out << " Difference : " << get_ele(*(xdiff ), 0) << " "
205 << get_ele(*(xdiff ), 1) << std::endl;
206 out << " =========================" << 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// ************************************************************
215TEUCHOS_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// ************************************************************
328TEUCHOS_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// ************************************************************
484TEUCHOS_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// ************************************************************
568TEUCHOS_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
SolutionHistory is basically a container of SolutionStates. SolutionHistory maintains a collection of...
Solution state for integrators and steppers. SolutionState contains the metadata for solutions and th...
TimeStepControl manages the time step size. There several mechanisms that effect the time step size a...
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)
void writeSolution(const std::string filename, Teuchos::RCP< const Tempus::SolutionHistory< Scalar > > solutionHistory)
TEUCHOS_UNIT_TEST(BackwardEuler, SinCos_ASA)
@ STORAGE_TYPE_STATIC
Keep a fix number of states.
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.