Tempus  Version of the Day
Time Integration
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
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 
16 namespace Tempus_Unit_Test {
17 
18 using Teuchos::RCP;
19 using Teuchos::rcp;
20 using Teuchos::rcp_const_cast;
21 using Teuchos::rcp_dynamic_cast;
23 using Teuchos::sublist;
24 
25 
26 // ************************************************************
27 // ************************************************************
28 TEUCHOS_UNIT_TEST(EDIRK_TrapezoidalRule, Default_Construction)
29 {
32 
33  // Test stepper properties.
34  TEUCHOS_ASSERT(stepper->getOrder() == 2);
35 }
36 
37 
38 // ************************************************************
39 // ************************************************************
40 TEUCHOS_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 // ************************************************************
49 TEUCHOS_UNIT_TEST(EDIRK_TrapezoidalRule, AppAction)
50 {
52  auto model = rcp(new Tempus_Test::SinCosModel<double>());
53  testRKAppAction(stepper, model, out, success);
54 }
55 
56 
57 // ************************************************************
58 // ************************************************************
59 
60 class StepperRKModifierEDIRK_TrapezoidaTest
61  : virtual public Tempus::StepperRKModifierBase<double>
62 {
63 public:
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(
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) {
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  }
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  }
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  }
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 
179 private:
180 
181  Teuchos::FancyOStream & out;
182  bool & success;
183 };
184 
185 
186 // ************************************************************
187 // ************************************************************
188 TEUCHOS_UNIT_TEST(DIRK, EDIRK_Trapezoida_Modifier)
189 {
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_Unit_Test
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.
void testFactoryConstruction(std::string stepperType, const Teuchos::RCP< const Thyra::ModelEvaluator< double > > &model)
Unit test utility for Stepper construction through StepperFactory.
The classic Dahlquist Test Problem.
#define TEST_COMPARE(v1, comp, v2)
#define TEST_FLOATING_EQUALITY(v1, v2, tol)
#define TEST_ASSERT(v1)
Sine-Cosine model problem from Rythmos. This is a canonical Sine-Cosine differential equation with a...
RK Trapezoidal Rule (A.K.A. RK Crank-Nicolson)
Teuchos::RCP< SolutionHistory< Scalar > > createSolutionHistoryME(const Teuchos::RCP< const Thyra::ModelEvaluator< Scalar > > &model)
Nonmember contructor from a Thyra ModelEvaluator.
Base class for Runge-Kutta methods, ExplicitRK, DIRK and IMEX.
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
TEUCHOS_UNIT_TEST(BackwardEuler, Default_Construction)
SolutionHistory is basically a container of SolutionStates. SolutionHistory maintains a collection of...
ACTION_LOCATION
Indicates the location of application action (see algorithm).
void testDIRKAccessorsFullConstruction(const RCP< Tempus::StepperDIRK< double > > &stepper)
Unit test utility for ExplicitRK Stepper construction and accessors.
#define TEUCHOS_ASSERT(assertion_test)
#define TEUCHOS_TEST_FOR_EXCEPT(throw_exception_test)