Tempus  Version of the Day
Time Integration
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
Tempus_StepperTrapezoidal_impl.hpp
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 #ifndef Tempus_StepperTrapezoidal_impl_hpp
10 #define Tempus_StepperTrapezoidal_impl_hpp
11 
13 #include "Tempus_WrapperModelEvaluatorBasic.hpp"
14 
15 namespace Tempus {
16 
17 template <class Scalar>
19 {
20  this->setStepperName("Trapezoidal Method");
21  this->setStepperType("Trapezoidal Method");
22  this->setUseFSAL(true);
23  this->setICConsistency("Consistent");
24  this->setICConsistencyCheck(false);
25  this->setZeroInitialGuess(false);
26 
27  this->setAppAction(Teuchos::null);
28  this->setDefaultSolver();
29 }
30 
31 template <class Scalar>
33  const Teuchos::RCP<const Thyra::ModelEvaluator<Scalar> >& appModel,
35  bool useFSAL, std::string ICConsistency, bool ICConsistencyCheck,
36  bool zeroInitialGuess,
38  stepperTrapAppAction)
39 {
40  this->setStepperName("Trapezoidal Method");
41  this->setStepperType("Trapezoidal Method");
42  this->setUseFSAL(useFSAL);
43  this->setICConsistency(ICConsistency);
44  this->setICConsistencyCheck(ICConsistencyCheck);
45  this->setZeroInitialGuess(zeroInitialGuess);
46 
47  this->setAppAction(stepperTrapAppAction);
48  this->setSolver(solver);
49 
50  if (appModel != Teuchos::null) {
51  this->setModel(appModel);
52  this->initialize();
53  }
54 }
55 
56 template <class Scalar>
59 {
60  if (appAction == Teuchos::null) {
61  // Create default appAction
62  stepperTrapAppAction_ =
64  }
65  else {
66  stepperTrapAppAction_ = appAction;
67  }
68  this->isInitialized_ = false;
69 }
70 
71 template <class Scalar>
73  const Teuchos::RCP<SolutionHistory<Scalar> >& solutionHistory)
74 {
75  using Teuchos::RCP;
76 
77  RCP<SolutionState<Scalar> > initialState = solutionHistory->getCurrentState();
78 
79  // Check if we need Stepper storage for xDot
80  if (initialState->getXDot() == Teuchos::null)
81  this->setStepperXDot(initialState->getX()->clone_v());
82  else
83  this->setStepperXDot(initialState->getXDot());
84 
86 
88  !(this->getUseFSAL()), std::logic_error,
89  "Error - The First-Same-As-Last (FSAL) principle is required\n"
90  " for the Trapezoidal Stepper (i.e., useFSAL=true)!\n");
91  // There are at least two ways around this, but are not implemented.
92  // - Do a solve for xDotOld, xDot_{n-1}, at each time step as for the
93  // initial conditions. This is expensive since you would be doing
94  // two solves every time step.
95  // - Use evaluateExplicitODE to get xDot_{n-1} if the application
96  // provides it. Explicit evaluations are cheaper but requires the
97  // application to implement xDot = f(x,t).
98 }
99 
100 template <class Scalar>
102  const Teuchos::RCP<SolutionHistory<Scalar> >& solutionHistory)
103 {
104  this->checkInitialized();
105 
106  using Teuchos::RCP;
107 
108  TEMPUS_FUNC_TIME_MONITOR("Tempus::StepperTrapezoidal::takeStep()");
109  {
111  solutionHistory->getNumStates() < 2, std::logic_error,
112  "Error - StepperTrapezoidal<Scalar>::takeStep(...)\n"
113  << "Need at least two SolutionStates for Trapezoidal.\n"
114  << " Number of States = " << solutionHistory->getNumStates()
115  << "\nTry setting in \"Solution History\" \"Storage Type\" = "
116  << "\"Undo\"\n"
117  << " or \"Storage Type\" = \"Static\" and \"Storage Limit\" = "
118  << "\"2\"\n");
119  RCP<StepperTrapezoidal<Scalar> > thisStepper = Teuchos::rcpFromRef(*this);
120  stepperTrapAppAction_->execute(
121  solutionHistory, thisStepper,
123 
124  RCP<SolutionState<Scalar> > workingState =
125  solutionHistory->getWorkingState();
126  RCP<SolutionState<Scalar> > currentState =
127  solutionHistory->getCurrentState();
128 
129  RCP<const Thyra::VectorBase<Scalar> > xOld = currentState->getX();
130  RCP<const Thyra::VectorBase<Scalar> > xDotOld = currentState->getXDot();
131  RCP<Thyra::VectorBase<Scalar> > x = workingState->getX();
132  if (workingState->getXDot() != Teuchos::null)
133  this->setStepperXDot(workingState->getXDot());
134  RCP<Thyra::VectorBase<Scalar> > xDot = this->getStepperXDot();
135 
136  const Scalar time = workingState->getTime();
137  const Scalar dt = workingState->getTimeStep();
138  const Scalar alpha = getAlpha(dt);
139  const Scalar beta = getBeta(dt);
140 
141  // Setup TimeDerivative
143  new StepperTrapezoidalTimeDerivative<Scalar>(alpha, xOld, xDotOld));
144 
145  auto p = Teuchos::rcp(
146  new ImplicitODEParameters<Scalar>(timeDer, dt, alpha, beta));
147  stepperTrapAppAction_->execute(
148  solutionHistory, thisStepper,
150 
151  const Thyra::SolveStatus<Scalar> sStatus =
152  this->solveImplicitODE(x, xDot, time, p);
153  stepperTrapAppAction_->execute(
154  solutionHistory, thisStepper,
156 
157  if (workingState->getXDot() != Teuchos::null) timeDer->compute(x, xDot);
158 
159  workingState->setSolutionStatus(sStatus); // Converged --> pass.
160  workingState->setOrder(this->getOrder());
161  workingState->computeNorms(currentState);
162  stepperTrapAppAction_->execute(
163  solutionHistory, thisStepper,
165  }
166  return;
167 }
168 
175 template <class Scalar>
178 {
180  rcp(new StepperState<Scalar>(this->getStepperType()));
181  return stepperState;
182 }
183 
184 template <class Scalar>
186  Teuchos::FancyOStream& out, const Teuchos::EVerbosityLevel verbLevel) const
187 {
188  out.setOutputToRootOnly(0);
189  out << std::endl;
190  Stepper<Scalar>::describe(out, verbLevel);
191  StepperImplicit<Scalar>::describe(out, verbLevel);
192 
193  out << "--- StepperTrapezoidal ---\n";
194  out << " stepperTrapAppAction_ = " << stepperTrapAppAction_ << std::endl;
195  out << "--------------------------" << std::endl;
196 }
197 
198 template <class Scalar>
200 {
201  out.setOutputToRootOnly(0);
202  bool isValidSetup = true;
203 
204  if (!Stepper<Scalar>::isValidSetup(out)) isValidSetup = false;
205  if (!StepperImplicit<Scalar>::isValidSetup(out)) isValidSetup = false;
206 
207  if (stepperTrapAppAction_ == Teuchos::null) {
208  isValidSetup = false;
209  out << "The Trapezoidal AppAction is not set!\n";
210  }
211 
212  return isValidSetup;
213 }
214 
215 // Nonmember constructor - ModelEvaluator and ParameterList
216 // ------------------------------------------------------------------------
217 template <class Scalar>
219  const Teuchos::RCP<const Thyra::ModelEvaluator<Scalar> >& model,
221 {
222  auto stepper = Teuchos::rcp(new StepperTrapezoidal<Scalar>());
223  stepper->setStepperImplicitValues(pl);
224 
225  if (model != Teuchos::null) {
226  stepper->setModel(model);
227  stepper->initialize();
228  }
229 
230  return stepper;
231 }
232 
233 } // namespace Tempus
234 #endif // Tempus_StepperTrapezoidal_impl_hpp
virtual Teuchos::RCP< Tempus::StepperState< Scalar > > getDefaultStepperState()
Get a default (initial) StepperState.
Trapezoidal method time stepper.
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
virtual void describe(Teuchos::FancyOStream &out, const Teuchos::EVerbosityLevel verbLevel) const override
Thyra Base interface for time steppers.
virtual void setInitialConditions(const Teuchos::RCP< SolutionHistory< Scalar > > &solutionHistory) override
Set the initial conditions and make them consistent.
Teuchos::RCP< StepperTrapezoidal< Scalar > > createStepperTrapezoidal(const Teuchos::RCP< const Thyra::ModelEvaluator< Scalar > > &model, Teuchos::RCP< Teuchos::ParameterList > pl)
Nonmember constructor - ModelEvaluator and ParameterList.
virtual void setInitialConditions(const Teuchos::RCP< SolutionHistory< Scalar > > &solutionHistory)
Set the initial conditions and make them consistent.
StepperState is a simple class to hold state information about the stepper.
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
Thyra Base interface for implicit time steppers.
virtual void describe(Teuchos::FancyOStream &out, const Teuchos::EVerbosityLevel verbLevel) const
Application Action for StepperTrapezoidal.
virtual void setAppAction(Teuchos::RCP< StepperTrapezoidalAppAction< Scalar > > appAction)
Time-derivative interface for Trapezoidal method.
SolutionHistory is basically a container of SolutionStates. SolutionHistory maintains a collection of...
basic_FancyOStream & setOutputToRootOnly(const int rootRank)
virtual bool isValidSetup(Teuchos::FancyOStream &out) const
virtual void describe(Teuchos::FancyOStream &out, const Teuchos::EVerbosityLevel verbLevel) const
virtual void takeStep(const Teuchos::RCP< SolutionHistory< Scalar > > &solutionHistory)
Take the specified timestep, dt, and return true if successful.