Tempus  Version of the Day
Time Integration
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
Tempus_StepperBDF2_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_StepperBDF2_impl_hpp
10 #define Tempus_StepperBDF2_impl_hpp
11 
12 #include "Tempus_WrapperModelEvaluatorBasic.hpp"
14 #include "Tempus_StepperFactory.hpp"
15 
16 
17 namespace Tempus {
18 
19 
20 template<class Scalar>
22 {
23  this->setStepperName( "BDF2");
24  this->setStepperType( "BDF2");
25  this->setUseFSAL( false);
26  this->setICConsistency( "None");
27  this->setICConsistencyCheck( false);
28  this->setZeroInitialGuess( false);
29 
30  this->setAppAction(Teuchos::null);
31  this->setDefaultSolver();
32  this->setStartUpStepper("DIRK 1 Stage Theta Method");
33 
34 }
35 
36 template<class Scalar>
38  const Teuchos::RCP<const Thyra::ModelEvaluator<Scalar> >& appModel,
40  const Teuchos::RCP<Stepper<Scalar> >& startUpStepper,
41  bool useFSAL,
42  std::string ICConsistency,
43  bool ICConsistencyCheck,
44  bool zeroInitialGuess,
45  const Teuchos::RCP<StepperBDF2AppAction<Scalar> >& stepperBDF2AppAction)
46 {
47  this->setStepperName( "BDF2");
48  this->setStepperType( "BDF2");
49  this->setUseFSAL( useFSAL);
50  this->setICConsistency( ICConsistency);
51  this->setICConsistencyCheck( ICConsistencyCheck);
52  this->setZeroInitialGuess( zeroInitialGuess);
53 
54  this->setAppAction(stepperBDF2AppAction);
55  this->setSolver(solver);
56  this->setStartUpStepper(startUpStepper);
57 
58  if (appModel != Teuchos::null) {
59  this->setModel(appModel);
60  this->initialize();
61  }
62 }
63 
64 template<class Scalar>
66  const Teuchos::RCP<const Thyra::ModelEvaluator<Scalar> >& appModel)
67 {
69  // If the startUpStepper's model is not set, set it to the stepper model.
70  if (startUpStepper_->getModel() == Teuchos::null) {
71  startUpStepper_->setModel(appModel);
72  startUpStepper_->initialize();
73  }
74 
75  this->isInitialized_ = false;
76 }
77 
78 
79 template<class Scalar>
80 void StepperBDF2<Scalar>::setStartUpStepper(std::string startupStepperType)
81 {
82  auto sf = Teuchos::rcp(new StepperFactory<Scalar>());
83  if (this->wrapperModel_ != Teuchos::null &&
84  this->wrapperModel_->getAppModel() != Teuchos::null) {
85  startUpStepper_ = sf->createStepper(startupStepperType,
86  this->wrapperModel_->getAppModel());
87  } else {
88  startUpStepper_ = sf->createStepper(startupStepperType);
89  }
90 
91  this->isInitialized_ = false;
92 }
93 
94 
95 template<class Scalar>
97  Teuchos::RCP<Stepper<Scalar> > startUpStepper)
98 {
99  startUpStepper_ = startUpStepper;
100 
101  if (this->wrapperModel_ != Teuchos::null) {
103  this->wrapperModel_->getAppModel() == Teuchos::null, std::logic_error,
104  "Error - Can not set the startUpStepper to Teuchos::null.\n");
105 
106  if (startUpStepper->getModel() == Teuchos::null &&
107  this->wrapperModel_->getAppModel() != Teuchos::null) {
108  startUpStepper_->setModel(this->wrapperModel_->getAppModel());
109  startUpStepper_->initialize();
110  }
111  }
112 
113  this->isInitialized_ = false;
114 }
115 
116 template<class Scalar>
119 {
120  if (appAction == Teuchos::null) {
121  // Create default appAction
122  stepperBDF2AppAction_ =
124  }
125  else {
126  stepperBDF2AppAction_ = appAction;
127  }
128 }
129 
130 
131 
132 template<class Scalar>
134 {
136 }
137 
138 
139 template<class Scalar>
141  const Teuchos::RCP<SolutionHistory<Scalar> >& solutionHistory)
142 {
143  using Teuchos::RCP;
144 
145  RCP<SolutionState<Scalar> > initialState = solutionHistory->getCurrentState();
146 
147  // Check if we need Stepper storage for xDot
148  if (initialState->getXDot() == Teuchos::null)
149  this->setStepperXDot(initialState->getX()->clone_v());
150  else
151  this->setStepperXDot(initialState->getXDot());
152 
154 }
155 
156 
157 template<class Scalar>
159  const Teuchos::RCP<SolutionHistory<Scalar> >& solutionHistory)
160 {
161  this->checkInitialized();
162 
163  using Teuchos::RCP;
164 
165  TEMPUS_FUNC_TIME_MONITOR("Tempus::StepperBDF2::takeStep()");
166  {
167  int numStates = solutionHistory->getNumStates();
168 
169  RCP<Thyra::VectorBase<Scalar> > xOld;
170  RCP<Thyra::VectorBase<Scalar> > xOldOld;
171 
172  // If there are less than 3 states (e.g., first time step), call
173  // startup stepper and return.
174  if (numStates < 3) {
175  computeStartUp(solutionHistory);
176  return;
177  }
178  TEUCHOS_TEST_FOR_EXCEPTION( (numStates < 3), std::logic_error,
179  "Error in Tempus::StepperBDF2::takeStep(): numStates after \n"
180  << "startup stepper must be at least 3, whereas numStates = "
181  << numStates <<"!\n" << "If running with Storage Type = Static, "
182  << "make sure Storage Limit > 2.\n");
183 
184  //IKT, FIXME: add error checking regarding states being consecutive and
185  //whether interpolated states are OK to use.
186 
187  RCP<StepperBDF2<Scalar> > thisStepper = Teuchos::rcpFromRef(*this);
188  stepperBDF2AppAction_->execute(solutionHistory, thisStepper,
190 
191 
192  RCP<SolutionState<Scalar> > workingState=solutionHistory->getWorkingState();
193  RCP<SolutionState<Scalar> > currentState=solutionHistory->getCurrentState();
194 
195  RCP<Thyra::VectorBase<Scalar> > x = workingState->getX();
196  if (workingState->getXDot() != Teuchos::null)
197  this->setStepperXDot(workingState->getXDot());
198  RCP<Thyra::VectorBase<Scalar> > xDot = this->getStepperXDot();
199 
200  //get time, dt and dtOld
201  const Scalar time = workingState->getTime();
202  const Scalar dt = workingState->getTimeStep();
203  const Scalar dtOld = currentState->getTimeStep();
204 
205  xOld = solutionHistory->getStateTimeIndexNM1()->getX();
206  xOldOld = solutionHistory->getStateTimeIndexNM2()->getX();
207  order_ = Scalar(2.0);
208 
209  // Setup TimeDerivative
211  Teuchos::rcp(new StepperBDF2TimeDerivative<Scalar>(dt, dtOld, xOld, xOldOld));
212 
213  const Scalar alpha = getAlpha(dt, dtOld);
214  const Scalar beta = getBeta (dt);
215 
217  timeDer, dt, alpha, beta));
218  stepperBDF2AppAction_->execute(solutionHistory, thisStepper,
220 
221  const Thyra::SolveStatus<Scalar> sStatus =
222  this->solveImplicitODE(x, xDot, time, p);
223 
224  stepperBDF2AppAction_->execute(solutionHistory, thisStepper,
226 
227 
228  if (workingState->getXDot() != Teuchos::null)
229  timeDer->compute(x, xDot);
230 
231  workingState->setSolutionStatus(sStatus); // Converged --> pass.
232  workingState->setOrder(getOrder());
233  workingState->computeNorms(currentState);
234  stepperBDF2AppAction_->execute(solutionHistory, thisStepper,
236  }
237  return;
238 }
239 
240 template<class Scalar>
242  const Teuchos::RCP<SolutionHistory<Scalar> >& solutionHistory)
243 {
244  auto out = Teuchos::fancyOStream( this->getOStream() );
245  out->setOutputToRootOnly(0);
246  Teuchos::OSTab ostab(out,1,"StepperBDF2::computeStartUp()");
247  *out << "Warning -- Taking a startup step for BDF2 using '"
248  << startUpStepper_->getStepperType()<<"'!" << std::endl;
249 
250  //Take one step using startUpStepper_
251  startUpStepper_->takeStep(solutionHistory);
252 
253  order_ = startUpStepper_->getOrder();
254 }
255 
262 template<class Scalar>
266 {
268  rcp(new StepperState<Scalar>(this->getStepperType()));
269  return stepperState;
270 }
271 
272 
273 template<class Scalar>
276  const Teuchos::EVerbosityLevel verbLevel ) const
277 {
278  auto l_out = Teuchos::fancyOStream( out.getOStream() );
279  Teuchos::OSTab ostab(*l_out, 2, this->description());
280  l_out->setOutputToRootOnly(0);
281  *l_out << std::endl;
282  Stepper<Scalar>::describe(out, verbLevel);
283  StepperImplicit<Scalar>::describe(out, verbLevel);
284 
285  *l_out << "--- StepperBDF2 ---\n";
286  if (startUpStepper_ != Teuchos::null) {
287  *l_out << " startup stepper type = "
288  << startUpStepper_->description() << std::endl;
289  }
290  *l_out << " startUpStepper_ = "
291  << startUpStepper_ << std::endl;
292  *l_out << " startUpStepper_->isInitialized() = "
293  << Teuchos::toString(startUpStepper_->isInitialized()) << std::endl;
294  *l_out << " stepperBDF2AppAction_ = "
295  << stepperBDF2AppAction_ << std::endl;
296  *l_out << "----------------------------" << std::endl;
297  *l_out << " order_ = " << order_ << std::endl;
298  *l_out << "-------------------" << std::endl;
299 }
300 
301 
302 template<class Scalar>
304 {
305  bool isValidSetup = true;
306  auto l_out = Teuchos::fancyOStream( out.getOStream() );
307  l_out->setOutputToRootOnly(0);
308 
309  if ( !Stepper<Scalar>::isValidSetup(out) ) isValidSetup = false;
310  if ( !StepperImplicit<Scalar>::isValidSetup(out) ) isValidSetup = false;
311 
312  if ( !this->startUpStepper_->isInitialized() ) {
313  isValidSetup = false;
314  *l_out << "The startup stepper is not initialized!\n";
315  }
316  if (stepperBDF2AppAction_ == Teuchos::null) {
317  isValidSetup = false;
318  *l_out << "The BDF2 AppAction is not set!\n";
319  }
320  return isValidSetup;
321 }
322 
323 
324 template<class Scalar>
327 {
328  auto pl = this->getValidParametersBasicImplicit();
329  pl->set("Start Up Stepper Type", startUpStepper_->getStepperType());
330  return pl;
331 }
332 
333 
334 // Nonmember constructor - ModelEvaluator and ParameterList
335 // ------------------------------------------------------------------------
336 template<class Scalar>
339  const Teuchos::RCP<const Thyra::ModelEvaluator<Scalar> >& model,
341 {
342  auto stepper = Teuchos::rcp(new StepperBDF2<Scalar>());
343  stepper->setStepperImplicitValues(pl);
344 
345  std::string startUpStepperName = "DIRK 1 Stage Theta Method";
346  if (pl != Teuchos::null) startUpStepperName =
347  pl->get<std::string>("Start Up Stepper Type", startUpStepperName);
348  stepper->setStartUpStepper(startUpStepperName);
349 
350  if (model != Teuchos::null) {
351  stepper->setModel(model);
352  stepper->initialize();
353  }
354 
355  return stepper;
356 }
357 
358 
359 } // namespace Tempus
360 #endif // Tempus_StepperBDF2_impl_hpp
BDF2 (Backward-Difference-Formula-2) time stepper.
virtual void setInitialConditions(const Teuchos::RCP< SolutionHistory< Scalar > > &solutionHistory)
Set the initial conditions and make them consistent.
Time-derivative interface for BDF2.
virtual Teuchos::RCP< Tempus::StepperState< Scalar > > getDefaultStepperState()
Get a default (initial) StepperState.
virtual void setModel(const Teuchos::RCP< const Thyra::ModelEvaluator< Scalar > > &appModel)
Set the model.
T & get(const std::string &name, T def_value)
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
virtual void describe(Teuchos::FancyOStream &out, const Teuchos::EVerbosityLevel verbLevel) const
virtual void initialize()
Initialize after construction and changing input parameters.
virtual void describe(Teuchos::FancyOStream &out, const Teuchos::EVerbosityLevel verbLevel) const override
virtual void computeStartUp(const Teuchos::RCP< SolutionHistory< Scalar > > &solutionHistory)
Compute the first time step given the supplied startup stepper.
virtual void setModel(const Teuchos::RCP< const Thyra::ModelEvaluator< Scalar > > &appModel) override
Set the model.
Thyra Base interface for time steppers.
virtual void setInitialConditions(const Teuchos::RCP< SolutionHistory< Scalar > > &solutionHistory) override
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
virtual void setAppAction(Teuchos::RCP< StepperBDF2AppAction< Scalar > > appAction)
SolutionHistory is basically a container of SolutionStates. SolutionHistory maintains a collection of...
virtual bool isValidSetup(Teuchos::FancyOStream &out) const
RCP< std::basic_ostream< char_type, traits_type > > getOStream()
Teuchos::RCP< StepperBDF2< Scalar > > createStepperBDF2(const Teuchos::RCP< const Thyra::ModelEvaluator< Scalar > > &model, Teuchos::RCP< Teuchos::ParameterList > pl)
Nonmember constructor - ModelEvaluator and ParameterList.
Application Action for StepperBDF2.
Teuchos::RCP< const Teuchos::ParameterList > getValidParameters() const
virtual void takeStep(const Teuchos::RCP< SolutionHistory< Scalar > > &solutionHistory)
Take the specified timestep, dt, and return true if successful.
void setStartUpStepper(std::string startupStepperType)
Set the stepper to use in first step.
virtual void initialize()
Initialize during construction and after changing input parameters.
StepperBDF2()
Default constructor.
std::string toString(const T &t)