Tempus Version of the Day
Time Integration
Loading...
Searching...
No Matches
Tempus_UnitTest_NewmarkImplicitDForm.cpp
Go to the documentation of this file.
1//@HEADER
2// *****************************************************************************
3// Tempus: Time Integration and Sensitivity Analysis Package
4//
5// Copyright 2017 NTESS and the Tempus contributors.
6// SPDX-License-Identifier: BSD-3-Clause
7// *****************************************************************************
8//@HEADER
9
11
12#include "Teuchos_XMLParameterListHelpers.hpp"
13
14#include "Tempus_TimeStepControl.hpp"
15
16#include "Tempus_StepperNewmarkImplicitDForm.hpp"
21
22#include "../TestModels/HarmonicOscillatorModel.hpp"
23
24namespace Tempus_Unit_Test {
25
26using Teuchos::ParameterList;
27using Teuchos::RCP;
28using Teuchos::rcp;
29using Teuchos::rcp_const_cast;
30using Teuchos::rcp_dynamic_cast;
31using Teuchos::sublist;
32
34
35// ************************************************************
36// ************************************************************
37TEUCHOS_UNIT_TEST(NewmarkImplicitDForm, Constructing_From_Defaults)
38{
39 double dt = 1.0;
40
41 // Read params from .xml file
42 RCP<ParameterList> pList = Teuchos::getParametersFromXmlFile(
43 "Tempus_NewmarkImplicitDForm_HarmonicOscillator_Damped_SecondOrder.xml");
44 RCP<ParameterList> pl = sublist(pList, "Tempus", true);
45
46 // Setup the HarmonicOscillatorModel
47 RCP<ParameterList> hom_pl = sublist(pList, "HarmonicOscillatorModel", true);
48 auto model = Teuchos::rcp(
50 auto modelME = rcp_dynamic_cast<const Thyra::ModelEvaluator<double>>(model);
51
52 // Setup Stepper for field solve ----------------------------
53 auto stepper =
54 Tempus::createStepperNewmarkImplicitDForm(modelME, Teuchos::null);
55
56 // Setup TimeStepControl ------------------------------------
57 RCP<Tempus::TimeStepControl<double>> timeStepControl =
58 Teuchos::rcp(new Tempus::TimeStepControl<double>());
59 ParameterList tscPL =
60 pl->sublist("Default Integrator").sublist("Time Step Control");
61 timeStepControl->setInitIndex(tscPL.get<int>("Initial Time Index"));
62 timeStepControl->setInitTime(tscPL.get<double>("Initial Time"));
63 timeStepControl->setFinalTime(tscPL.get<double>("Final Time"));
64 timeStepControl->setInitTimeStep(dt);
65 timeStepControl->initialize();
66
67 // Setup initial condition SolutionState --------------------
68 using Teuchos::rcp_const_cast;
69 auto inArgsIC = model->getNominalValues();
70 RCP<Thyra::VectorBase<double>> icX =
71 rcp_const_cast<Thyra::VectorBase<double>>(inArgsIC.get_x());
72 RCP<Thyra::VectorBase<double>> icXDot =
73 rcp_const_cast<Thyra::VectorBase<double>>(inArgsIC.get_x_dot());
74 RCP<Thyra::VectorBase<double>> icXDotDot =
75 rcp_const_cast<Thyra::VectorBase<double>>(inArgsIC.get_x_dot_dot());
76 RCP<Tempus::SolutionState<double>> icState =
77 Tempus::createSolutionStateX(icX, icXDot, icXDotDot);
78 icState->setTime(timeStepControl->getInitTime());
79 icState->setIndex(timeStepControl->getInitIndex());
80 icState->setTimeStep(0.0);
81 icState->setOrder(stepper->getOrder());
82 icState->setSolutionStatus(Tempus::Status::PASSED); // ICs are passing.
83
84 // Setup SolutionHistory ------------------------------------
85 RCP<Tempus::SolutionHistory<double>> solutionHistory =
86 Teuchos::rcp(new Tempus::SolutionHistory<double>());
87 solutionHistory->setName("Forward States");
88 solutionHistory->setStorageType(Tempus::STORAGE_TYPE_STATIC);
89 solutionHistory->setStorageLimit(2);
90 solutionHistory->addState(icState);
91
92 // Setup Integrator -----------------------------------------
93 RCP<Tempus::IntegratorBasic<double>> integrator =
94 Tempus::createIntegratorBasic<double>();
95 integrator->setStepper(stepper);
96 integrator->setTimeStepControl(timeStepControl);
97 integrator->setSolutionHistory(solutionHistory);
98 // integrator->setObserver(...);
99 integrator->initialize();
100
101 // Integrate to timeMax
102 bool integratorStatus = integrator->advanceTime();
103 TEST_ASSERT(integratorStatus)
104
105 // Test if at 'Final Time'
106 double time = integrator->getTime();
107 double timeFinal = pl->sublist("Default Integrator")
108 .sublist("Time Step Control")
109 .get<double>("Final Time");
110 TEST_FLOATING_EQUALITY(time, timeFinal, 1.0e-14);
111
112 // Time-integrated solution and the exact solution
113 RCP<Thyra::VectorBase<double>> x = integrator->getX();
114 RCP<const Thyra::VectorBase<double>> x_exact =
115 model->getExactSolution(time).get_x();
116
117 // Calculate the error
118 RCP<Thyra::VectorBase<double>> xdiff = x->clone_v();
119 Thyra::V_StVpStV(xdiff.ptr(), 1.0, *x_exact, -1.0, *(x));
120
121 // Check the order and intercept
122 out << " Stepper = " << stepper->description() << std::endl;
123 out << " =========================" << std::endl;
124 out << " Exact solution : " << get_ele(*(x_exact), 0) << std::endl;
125 out << " Computed solution: " << get_ele(*(x), 0) << std::endl;
126 out << " Difference : " << get_ele(*(xdiff), 0) << std::endl;
127 out << " =========================" << std::endl;
128 TEST_FLOATING_EQUALITY(get_ele(*(x), 0), -0.222222, 1.0e-4);
129}
130
131// ************************************************************
132class StepperNewmarkImplicitDFormModifierTest
133 : virtual public Tempus::StepperNewmarkImplicitDFormModifierBase<double> {
134 public:
136 StepperNewmarkImplicitDFormModifierTest()
137 : testBEGIN_STEP(false),
138 testBEFORE_SOLVE(false),
139 testAFTER_SOLVE(false),
140 testEND_STEP(false),
141 testCurrentValue(-0.99),
142 testDt(-1.5),
143 testName("")
144 {
145 }
146
148 virtual ~StepperNewmarkImplicitDFormModifierTest() {}
149
151 virtual void modify(
152 Teuchos::RCP<Tempus::SolutionHistory<double>> sh,
155 double>::ACTION_LOCATION actLoc)
156 {
157 switch (actLoc) {
158 case StepperNewmarkImplicitDFormAppAction<double>::BEGIN_STEP: {
159 testBEGIN_STEP = true;
160 break;
161 }
162 case StepperNewmarkImplicitDFormAppAction<double>::BEFORE_SOLVE: {
163 testBEFORE_SOLVE = true;
164 testDt = sh->getWorkingState()->getTimeStep();
165 // sh->getWorkingState()->setTimeStep(testDt);
166 break;
167 }
168 case StepperNewmarkImplicitDFormAppAction<double>::AFTER_SOLVE: {
169 testAFTER_SOLVE = true;
170 testName = "Newmark Implicit A Form - Modifier";
171 stepper->setStepperName(testName);
172 break;
173 }
174 case StepperNewmarkImplicitDFormAppAction<double>::END_STEP: {
175 testEND_STEP = true;
176 auto x = sh->getWorkingState()->getX();
177 testCurrentValue = get_ele(*(x), 0);
178 break;
179 }
180 default:
181 TEUCHOS_TEST_FOR_EXCEPTION(true, std::logic_error,
182 "Error - unknown action location.\n");
183 }
184 }
185
186 bool testBEGIN_STEP;
187 bool testBEFORE_SOLVE;
188 bool testAFTER_SOLVE;
189 bool testEND_STEP;
190 double testCurrentValue;
191 double testDt;
192 std::string testName;
193};
194
195// ************************************************************
196// ************************************************************
197class StepperNewmarkImplicitDFormModifierXTest
198 : virtual public Tempus::StepperNewmarkImplicitDFormModifierXBase<double> {
199 public:
201 StepperNewmarkImplicitDFormModifierXTest()
202 : testX_BEGIN_STEP(false),
203 testX_BEFORE_SOLVE(false),
204 testX_AFTER_SOLVE(false),
205 testX_END_STEP(false),
206 testX(-0.99),
207 testXDot(-0.99),
208 testDt(-1.5),
209 testTime(-1.5)
210 {
211 }
212
214 virtual ~StepperNewmarkImplicitDFormModifierXTest() {}
215
217 virtual void modify(
218 Teuchos::RCP<Thyra::VectorBase<double>> x, const double time,
219 const double dt,
221 double>::MODIFIER_TYPE modType)
222 {
223 switch (modType) {
224 case StepperNewmarkImplicitDFormModifierXBase<double>::X_BEGIN_STEP: {
225 testX_BEGIN_STEP = true;
226 testX = get_ele(*(x), 0);
227 break;
228 }
229 case StepperNewmarkImplicitDFormModifierXBase<double>::X_BEFORE_SOLVE: {
230 testX_BEFORE_SOLVE = true;
231 testDt = dt;
232 break;
233 }
234 case StepperNewmarkImplicitDFormModifierXBase<double>::X_AFTER_SOLVE: {
235 testX_AFTER_SOLVE = true;
236 testTime = time;
237 testX = get_ele(*(x), 0);
238 break;
239 }
240 case StepperNewmarkImplicitDFormModifierXBase<double>::X_END_STEP: {
241 testX_END_STEP = true;
242 testX = get_ele(*(x), 0);
243 break;
244 }
245 default:
246 TEUCHOS_TEST_FOR_EXCEPTION(true, std::logic_error,
247 "Error - unknown action location.\n");
248 }
249 }
250
251 bool testX_BEGIN_STEP;
252 bool testX_BEFORE_SOLVE;
253 bool testX_AFTER_SOLVE;
254 bool testX_END_STEP;
255 double testX;
256 double testXDot;
257 double testDt;
258 double testTime;
259};
260
261// ************************************************************
262// ************************************************************
263TEUCHOS_UNIT_TEST(NewmarkImplicitDForm, Default_Construction)
264{
265 auto model = rcp(new Tempus_Test::HarmonicOscillatorModel<double>());
266
267 // Default construction.
268 auto stepper = rcp(new Tempus::StepperNewmarkImplicitDForm<double>());
269 stepper->setModel(model);
270 stepper->initialize();
271 TEUCHOS_TEST_FOR_EXCEPT(!stepper->isInitialized());
272
273 auto modifier =
275 auto modifierX =
277
278 // Default values for construction.
279 auto solver = rcp(new Thyra::NOXNonlinearSolver());
280 solver->setParameterList(Tempus::defaultSolverParameters());
281
282 bool useFSAL = stepper->getUseFSAL();
283 std::string ICConsistency = stepper->getICConsistency();
284 bool ICConsistencyCheck = stepper->getICConsistencyCheck();
285 bool zeroInitialGuess = stepper->getZeroInitialGuess();
286 std::string schemeName = "Average Acceleration";
287 double beta = 0.25;
288 double gamma = 0.5;
289
290 // Test the set functions.
291 stepper->setAppAction(modifier);
292 stepper->initialize();
293 TEUCHOS_TEST_FOR_EXCEPT(!stepper->isInitialized());
294 stepper->setAppAction(modifierX);
295 stepper->initialize();
296 TEUCHOS_TEST_FOR_EXCEPT(!stepper->isInitialized());
297 stepper->setSolver(solver);
298 stepper->initialize();
299 TEUCHOS_TEST_FOR_EXCEPT(!stepper->isInitialized());
300 stepper->setUseFSAL(useFSAL);
301 stepper->initialize();
302 TEUCHOS_TEST_FOR_EXCEPT(!stepper->isInitialized());
303 stepper->setICConsistency(ICConsistency);
304 stepper->initialize();
305 TEUCHOS_TEST_FOR_EXCEPT(!stepper->isInitialized());
306 stepper->setICConsistencyCheck(ICConsistencyCheck);
307 stepper->initialize();
308 TEUCHOS_TEST_FOR_EXCEPT(!stepper->isInitialized());
309 stepper->setZeroInitialGuess(zeroInitialGuess);
310 stepper->initialize();
311 TEUCHOS_TEST_FOR_EXCEPT(!stepper->isInitialized());
312
313 stepper->setSchemeName(schemeName);
314 stepper->initialize();
315 TEUCHOS_TEST_FOR_EXCEPT(!stepper->isInitialized());
316 stepper->setBeta(beta);
317 stepper->initialize();
318 TEUCHOS_TEST_FOR_EXCEPT(!stepper->isInitialized());
319 stepper->setGamma(gamma);
320 stepper->initialize();
321 TEUCHOS_TEST_FOR_EXCEPT(!stepper->isInitialized());
322
323 // Full argument list construction.
325 model, solver, useFSAL, ICConsistency, ICConsistencyCheck,
326 zeroInitialGuess, schemeName, beta, gamma, modifier));
327 TEUCHOS_TEST_FOR_EXCEPT(!stepper->isInitialized());
328
329 // Test stepper properties.
330 TEUCHOS_ASSERT(stepper->getOrder() == 2);
331}
332
333// ************************************************************
334// ************************************************************
335TEUCHOS_UNIT_TEST(NewmarkImplicitDForm, StepperFactory_Construction)
336{
337 auto model = rcp(new Tempus_Test::HarmonicOscillatorModel<double>());
338 testFactoryConstruction("Newmark Implicit d-Form", model);
339}
340
341TEUCHOS_UNIT_TEST(NewmarkImplicitDForm, AppAction_Modifier)
342{
343 using Teuchos::ParameterList;
344 using Teuchos::RCP;
345 using Teuchos::sublist;
346
347 double dt = 1.0;
348
349 // Read params from .xml file
350 RCP<ParameterList> pList = Teuchos::getParametersFromXmlFile(
351 "Tempus_NewmarkImplicitAForm_HarmonicOscillator_Damped_SecondOrder.xml");
352 RCP<ParameterList> pl = sublist(pList, "Tempus", true);
353
354 // Setup the HarmonicOscillatorModel
355 RCP<ParameterList> hom_pl = sublist(pList, "HarmonicOscillatorModel", true);
356 RCP<const Thyra::ModelEvaluator<double>> model =
357 Teuchos::rcp(new Tempus_Test::HarmonicOscillatorModel<double>(hom_pl));
358
359 // Setup Stepper for field solve ----------------------------
360 RCP<Tempus::StepperNewmarkImplicitDForm<double>> stepper =
361 Tempus::createStepperNewmarkImplicitDForm(model, Teuchos::null);
362
363 auto modifier = rcp(new StepperNewmarkImplicitDFormModifierTest());
364 stepper->setAppAction(modifier);
365 stepper->initialize();
366
367 // Setup TimeStepControl ------------------------------------
368 RCP<Tempus::TimeStepControl<double>> timeStepControl =
369 Teuchos::rcp(new Tempus::TimeStepControl<double>());
370 ParameterList tscPL =
371 pl->sublist("Default Integrator").sublist("Time Step Control");
372 timeStepControl->setInitIndex(tscPL.get<int>("Initial Time Index"));
373 timeStepControl->setInitTime(tscPL.get<double>("Initial Time"));
374 timeStepControl->setFinalTime(dt);
375 timeStepControl->setInitTimeStep(dt);
376 timeStepControl->initialize();
377
378 // Setup initial condition SolutionState --------------------
379 using Teuchos::rcp_const_cast;
380 auto inArgsIC = model->getNominalValues();
381 RCP<Thyra::VectorBase<double>> icX =
382 rcp_const_cast<Thyra::VectorBase<double>>(inArgsIC.get_x());
383 RCP<Thyra::VectorBase<double>> icXDot =
384 rcp_const_cast<Thyra::VectorBase<double>>(inArgsIC.get_x_dot());
385 RCP<Thyra::VectorBase<double>> icXDotDot =
386 rcp_const_cast<Thyra::VectorBase<double>>(inArgsIC.get_x_dot_dot());
387 RCP<Tempus::SolutionState<double>> icState =
388 Tempus::createSolutionStateX(icX, icXDot, icXDotDot);
389 icState->setTime(timeStepControl->getInitTime());
390 icState->setIndex(timeStepControl->getInitIndex());
391 icState->setTimeStep(0.0);
392 icState->setOrder(stepper->getOrder());
393 icState->setSolutionStatus(Tempus::Status::PASSED); // ICs are passing.
394
395 // Setup SolutionHistory ------------------------------------
396 RCP<Tempus::SolutionHistory<double>> solutionHistory =
397 Teuchos::rcp(new Tempus::SolutionHistory<double>());
398 solutionHistory->setName("Forward States");
399 solutionHistory->setStorageType(Tempus::STORAGE_TYPE_STATIC);
400 solutionHistory->setStorageLimit(2);
401 solutionHistory->addState(icState);
402
403 // Setup Integrator -----------------------------------------
404 RCP<Tempus::IntegratorBasic<double>> integrator =
405 Tempus::createIntegratorBasic<double>();
406 integrator->setStepper(stepper);
407 integrator->setTimeStepControl(timeStepControl);
408 integrator->setSolutionHistory(solutionHistory);
409 integrator->initialize();
410
411 // Integrate to timeMax
412 bool integratorStatus = integrator->advanceTime();
413 TEST_ASSERT(integratorStatus)
414
415 // Testing that each ACTION_LOCATION has been called.
416 TEST_COMPARE(modifier->testBEGIN_STEP, ==, true);
417 TEST_COMPARE(modifier->testBEFORE_SOLVE, ==, true);
418 TEST_COMPARE(modifier->testAFTER_SOLVE, ==, true);
419 TEST_COMPARE(modifier->testEND_STEP, ==, true);
420
421 // Testing that values can be set through the Modifier.
422 auto x = integrator->getX();
423 auto Dt = integrator->getTime();
424 TEST_FLOATING_EQUALITY(modifier->testDt, Dt, 1.0e-14);
425 TEST_FLOATING_EQUALITY(modifier->testCurrentValue, get_ele(*(x), 0), 1.0e-14);
426 TEST_COMPARE(modifier->testName, ==, stepper->getStepperName());
427}
428
429TEUCHOS_UNIT_TEST(NewmarkImplicitDForm, AppAction_ModifierX)
430{
431 using Teuchos::ParameterList;
432 using Teuchos::RCP;
433 using Teuchos::sublist;
434
435 double dt = 1.0;
436
437 // Read params from .xml file
438 RCP<ParameterList> pList = Teuchos::getParametersFromXmlFile(
439 "Tempus_NewmarkImplicitAForm_HarmonicOscillator_Damped_SecondOrder.xml");
440 RCP<ParameterList> pl = sublist(pList, "Tempus", true);
441
442 // Setup the HarmonicOscillatorModel
443 RCP<ParameterList> hom_pl = sublist(pList, "HarmonicOscillatorModel", true);
444 RCP<const Thyra::ModelEvaluator<double>> model =
445 Teuchos::rcp(new Tempus_Test::HarmonicOscillatorModel<double>(hom_pl));
446
447 // Setup Stepper for field solve ----------------------------
448 RCP<Tempus::StepperNewmarkImplicitDForm<double>> stepper =
449 Tempus::createStepperNewmarkImplicitDForm(model, Teuchos::null);
450
451 auto modifierX = rcp(new StepperNewmarkImplicitDFormModifierXTest());
452 stepper->setAppAction(modifierX);
453 stepper->initialize();
454
455 // Setup TimeStepControl ------------------------------------
456 RCP<Tempus::TimeStepControl<double>> timeStepControl =
457 Teuchos::rcp(new Tempus::TimeStepControl<double>());
458 ParameterList tscPL =
459 pl->sublist("Default Integrator").sublist("Time Step Control");
460 timeStepControl->setInitIndex(tscPL.get<int>("Initial Time Index"));
461 timeStepControl->setInitTime(tscPL.get<double>("Initial Time"));
462 timeStepControl->setFinalTime(dt);
463 timeStepControl->setInitTimeStep(dt);
464 timeStepControl->initialize();
465
466 // Setup initial condition SolutionState --------------------
467 using Teuchos::rcp_const_cast;
468 auto inArgsIC = model->getNominalValues();
469 RCP<Thyra::VectorBase<double>> icX =
470 rcp_const_cast<Thyra::VectorBase<double>>(inArgsIC.get_x());
471 RCP<Thyra::VectorBase<double>> icXDot =
472 rcp_const_cast<Thyra::VectorBase<double>>(inArgsIC.get_x_dot());
473 RCP<Thyra::VectorBase<double>> icXDotDot =
474 rcp_const_cast<Thyra::VectorBase<double>>(inArgsIC.get_x_dot_dot());
475 RCP<Tempus::SolutionState<double>> icState =
476 Tempus::createSolutionStateX(icX, icXDot, icXDotDot);
477 icState->setTime(timeStepControl->getInitTime());
478 icState->setIndex(timeStepControl->getInitIndex());
479 icState->setTimeStep(0.0);
480 icState->setOrder(stepper->getOrder());
481 icState->setSolutionStatus(Tempus::Status::PASSED); // ICs are passing.
482
483 // Setup SolutionHistory ------------------------------------
484 RCP<Tempus::SolutionHistory<double>> solutionHistory =
485 Teuchos::rcp(new Tempus::SolutionHistory<double>());
486 solutionHistory->setName("Forward States");
487 solutionHistory->setStorageType(Tempus::STORAGE_TYPE_STATIC);
488 solutionHistory->setStorageLimit(2);
489 solutionHistory->addState(icState);
490
491 // Setup Integrator -----------------------------------------
492 RCP<Tempus::IntegratorBasic<double>> integrator =
493 Tempus::createIntegratorBasic<double>();
494 integrator->setStepper(stepper);
495 integrator->setTimeStepControl(timeStepControl);
496 integrator->setSolutionHistory(solutionHistory);
497 integrator->initialize();
498
499 // Integrate to timeMax
500 bool integratorStatus = integrator->advanceTime();
501 TEST_ASSERT(integratorStatus)
502
503 // Testing that each ACTION_LOCATION has been called.
504 TEST_COMPARE(modifierX->testX_BEGIN_STEP, ==, true);
505 TEST_COMPARE(modifierX->testX_BEFORE_SOLVE, ==, true);
506 TEST_COMPARE(modifierX->testX_AFTER_SOLVE, ==, true);
507 TEST_COMPARE(modifierX->testX_END_STEP, ==, true);
508
509 // Testing that values can be set through the Modifier.
510 auto Dt = integrator->getTime();
511 TEST_FLOATING_EQUALITY(modifierX->testDt, Dt, 1.0e-14);
512
513 const auto x = integrator->getX();
514 TEST_FLOATING_EQUALITY(modifierX->testX, get_ele(*(x), 0), 1.0e-14);
515}
516
517} // namespace Tempus_Unit_Test
SolutionHistory is basically a container of SolutionStates. SolutionHistory maintains a collection of...
ACTION_LOCATION
Indicates the location of application action (see algorithm).
MODIFIER_TYPE
Indicates the location of application action (see algorithm).
virtual void modify(Teuchos::RCP< Thyra::VectorBase< double > >, const double, const double, const MODIFIER_TYPE modType)=0
Modify solution based on the MODIFIER_TYPE.
TimeStepControl manages the time step size. There several mechanisms that effect the time step size a...
TEUCHOS_UNIT_TEST(BackwardEuler, Default_Construction)
void testFactoryConstruction(std::string stepperType, const Teuchos::RCP< const Thyra::ModelEvaluator< double > > &model)
Unit test utility for Stepper construction through StepperFactory.
@ 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.
Teuchos::RCP< StepperNewmarkImplicitDForm< Scalar > > createStepperNewmarkImplicitDForm(const Teuchos::RCP< const Thyra::ModelEvaluator< Scalar > > &model, Teuchos::RCP< Teuchos::ParameterList > pl)
Nonmember constructor - ModelEvaluator and ParameterList.
Teuchos::RCP< Teuchos::ParameterList > defaultSolverParameters()
Returns the default solver ParameterList for implicit Steppers.