Tempus Version of the Day
Time Integration
Loading...
Searching...
No Matches
Tempus_UnitTest_EDIRK_TrapezoidalRule.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
10
12
13#include "../TestModels/DahlquistTestModel.hpp"
14
15
16namespace Tempus_Unit_Test {
17
18using Teuchos::RCP;
19using Teuchos::rcp;
20using Teuchos::rcp_const_cast;
21using Teuchos::rcp_dynamic_cast;
22using Teuchos::ParameterList;
23using Teuchos::sublist;
24
25
26// ************************************************************
27// ************************************************************
28TEUCHOS_UNIT_TEST(EDIRK_TrapezoidalRule, Default_Construction)
29{
30 auto stepper = rcp(new Tempus::StepperEDIRK_TrapezoidalRule<double>());
32
33 // Test stepper properties.
34 TEUCHOS_ASSERT(stepper->getOrder() == 2);
35}
36
37
38// ************************************************************
39// ************************************************************
40TEUCHOS_UNIT_TEST(EDIRK_TrapezoidalRule, StepperFactory_Construction)
41{
42 auto model = rcp(new Tempus_Test::SinCosModel<double>());
43 testFactoryConstruction("RK Trapezoidal Rule", model);
44}
45
46
47// ************************************************************
48// ************************************************************
49TEUCHOS_UNIT_TEST(EDIRK_TrapezoidalRule, AppAction)
50{
51 auto stepper = rcp(new Tempus::StepperEDIRK_TrapezoidalRule<double>());
52 auto model = rcp(new Tempus_Test::SinCosModel<double>());
53 testRKAppAction(stepper, model, out, success);
54}
55
56
57// ************************************************************
58// ************************************************************
59
60class StepperRKModifierEDIRK_TrapezoidaTest
61 : virtual public Tempus::StepperRKModifierBase<double>
62{
63public:
64
66 StepperRKModifierEDIRK_TrapezoidaTest(Teuchos::FancyOStream &Out, bool &Success)
67 : out(Out), success(Success)
68 {}
69
70 // FSAL
72 virtual ~StepperRKModifierEDIRK_TrapezoidaTest(){}
73
75 virtual void modify(
76 Teuchos::RCP<Tempus::SolutionHistory<double> > sh,
77 Teuchos::RCP<Tempus::StepperRKBase<double> > stepper,
79 {
80 const double relTol = 1.0e-14;
81 auto stageNumber = stepper->getStageNumber();
82 Teuchos::SerialDenseVector<int,double> c = stepper->getTableau()->c();
83
84 auto currentState = sh->getCurrentState();
85 auto workingState = sh->getWorkingState();
86 const double dt = workingState->getTimeStep();
87 double time = currentState->getTime();
88 if (stageNumber >= 0) time += c(stageNumber)*dt;
89
90 auto x = workingState->getX();
91 auto xDot = workingState->getXDot();
92 if (xDot == Teuchos::null) xDot = stepper->getStepperXDot();
93
94 switch (actLoc) {
95 case StepperRKAppAction<double>::BEGIN_STEP:
96 {
97 {
98 auto DME = Teuchos::rcp_dynamic_cast<
99 const Tempus_Test::DahlquistTestModel<double> > (stepper->getModel());
100 TEST_FLOATING_EQUALITY(DME->getLambda(), -1.0, relTol);
101 }
102 TEST_FLOATING_EQUALITY(dt, 1.0, relTol);
103
104 const double x_0 = get_ele(*(x), 0);
105 const double xDot_0 = get_ele(*(xDot), 0);
106 TEST_FLOATING_EQUALITY(x_0, 1.0, relTol); // Should be x_0
107 TEST_FLOATING_EQUALITY(xDot_0, -1.0, relTol); // Should be xDot_0
108 TEST_ASSERT(std::abs(time) < relTol);
109 TEST_FLOATING_EQUALITY(dt, 1.0, relTol);
110 TEST_COMPARE(stageNumber, ==, -1);
111 break;
112 }
113 case StepperRKAppAction<double>::BEGIN_STAGE:
114 case StepperRKAppAction<double>::BEFORE_SOLVE:
115 {
116 const double X_i = get_ele(*(x), 0);
117 const double f_i = get_ele(*(xDot), 0);
118
119 if (stageNumber == 0) {
120 TEST_FLOATING_EQUALITY(X_i, 1.0, relTol);
121 TEST_ASSERT(std::abs(f_i) < relTol);
122 TEST_ASSERT(std::abs(time) < relTol);
123 } else if (stageNumber == 1) {
124 TEST_FLOATING_EQUALITY(X_i, 1.0, relTol);
125 TEST_FLOATING_EQUALITY(f_i, -1.0, relTol);
126 TEST_FLOATING_EQUALITY(time, 1.0, relTol);
127 } else {
128 TEUCHOS_TEST_FOR_EXCEPT( !(-1 < stageNumber && stageNumber < 2));
129 }
130
131 break;
132 }
133 case StepperRKAppAction<double>::AFTER_SOLVE:
134 case StepperRKAppAction<double>::BEFORE_EXPLICIT_EVAL:
135 case StepperRKAppAction<double>::END_STAGE:
136 {
137 const double X_i = get_ele(*(x), 0);
138 const double f_i = get_ele(*(xDot), 0);
139
140 if (stageNumber == 0) {
141 // X_i = 1, f_1 = 0
142 TEST_FLOATING_EQUALITY(X_i, 1.0, relTol);
143 TEST_FLOATING_EQUALITY(f_i, -1.0, relTol);
144 TEST_ASSERT(std::abs(time) < relTol);
145 } else if (stageNumber == 1) {
146 // X_i = , f_i =
147 TEST_FLOATING_EQUALITY(X_i, 1.0/3.0, relTol);
148 TEST_FLOATING_EQUALITY(f_i, -1.0/3.0, relTol);
149 TEST_FLOATING_EQUALITY(time, 1.0, relTol);
150 } else {
151 TEUCHOS_TEST_FOR_EXCEPT( !(-1 < stageNumber && stageNumber < 2));
152 }
153
154 break;
155 }
156 case StepperRKAppAction<double>::END_STEP:
157 {
158 const double x_1 = get_ele(*(x), 0);
159 time = workingState->getTime();
160 TEST_FLOATING_EQUALITY(x_1, 1.0/3.0, relTol); // Should be x_1
161 TEST_FLOATING_EQUALITY(time, 1.0, relTol);
162 TEST_FLOATING_EQUALITY(dt, 1.0, relTol);
163 TEST_COMPARE(stageNumber, ==, -1);
164
165 if (stepper->getUseEmbedded() == true) {
166 TEST_FLOATING_EQUALITY(workingState->getTolRel(), 1.0, relTol);
167 TEST_ASSERT(std::abs(workingState->getTolAbs()) < relTol);
168 // e = 0 from doxygen above.
169 TEST_ASSERT(std::abs(workingState->getErrorRel()) < relTol);
170 }
171
172 }
173
174 }
175
176
177 }
178
179private:
180
181 Teuchos::FancyOStream & out;
182 bool & success;
183};
184
185
186// ************************************************************
187// ************************************************************
188TEUCHOS_UNIT_TEST(DIRK, EDIRK_Trapezoida_Modifier)
189{
190 auto stepper = rcp(new Tempus::StepperEDIRK_TrapezoidalRule<double>());
191 Teuchos::RCP<const Thyra::ModelEvaluator<double> >
193
194 auto modifier = rcp(new StepperRKModifierEDIRK_TrapezoidaTest(out, success));
195
196 stepper->setModel(model);
197 stepper->setAppAction(modifier);
198 stepper->setICConsistency("Consistent");
199 stepper->setUseFSAL(false);
200 stepper->initialize();
201
202 // Create a SolutionHistory.
203 auto solutionHistory = Tempus::createSolutionHistoryME(model);
204
205 // Take one time step.
206 stepper->setInitialConditions(solutionHistory);
207 solutionHistory->initWorkingState();
208 double dt = 1.0;
209 solutionHistory->getWorkingState()->setTimeStep(dt);
210 solutionHistory->getWorkingState()->setTime(dt);
211 stepper->takeStep(solutionHistory); // Primary testing occurs in
212 // modifierX during takeStep().
213 // Test stepper properties.
214 TEUCHOS_ASSERT(stepper->getOrder() == 2);
215}
216} // namespace Tempus_Test
SolutionHistory is basically a container of SolutionStates. SolutionHistory maintains a collection of...
RK Trapezoidal Rule (A.K.A. RK Crank-Nicolson)
ACTION_LOCATION
Indicates the location of application action (see algorithm).
Base class for Runge-Kutta methods, ExplicitRK, DIRK and IMEX.
The classic Dahlquist Test Problem.
Sine-Cosine model problem from Rythmos. This is a canonical Sine-Cosine differential equation.
void testDIRKAccessorsFullConstruction(const RCP< Tempus::StepperDIRK< double > > &stepper)
Unit test utility for ExplicitRK Stepper construction and accessors.
void testRKAppAction(const Teuchos::RCP< Tempus::StepperRKBase< double > > &stepper, const Teuchos::RCP< const Thyra::ModelEvaluator< double > > &model, Teuchos::FancyOStream &out, bool &success)
Unit test utility for Stepper RK AppAction.
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.
Teuchos::RCP< SolutionHistory< Scalar > > createSolutionHistoryME(const Teuchos::RCP< const Thyra::ModelEvaluator< Scalar > > &model)
Nonmember contructor from a Thyra ModelEvaluator.