Tempus Version of the Day
Time Integration
Loading...
Searching...
No Matches
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
17namespace Tempus {
18
19
20template<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
36template<class Scalar>
38 const Teuchos::RCP<const Thyra::ModelEvaluator<Scalar> >& appModel,
39 const Teuchos::RCP<Thyra::NonlinearSolverBase<Scalar> >& solver,
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
64template<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
79template<class Scalar>
80void 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
95template<class Scalar>
97 Teuchos::RCP<Stepper<Scalar> > startUpStepper)
98{
99 startUpStepper_ = startUpStepper;
100
101 if (this->wrapperModel_ != Teuchos::null) {
102 TEUCHOS_TEST_FOR_EXCEPTION(
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
116template<class Scalar>
118 Teuchos::RCP<StepperBDF2AppAction<Scalar> > appAction)
119{
120 if (appAction == Teuchos::null) {
121 // Create default appAction
122 stepperBDF2AppAction_ =
123 Teuchos::rcp(new StepperBDF2ModifierDefault<Scalar>());
124 }
125 else {
126 stepperBDF2AppAction_ = appAction;
127 }
128}
129
130
131
132template<class Scalar>
137
138
139template<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
157template<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
210 Teuchos::RCP<TimeDerivative<Scalar> > timeDer =
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
216 auto p = Teuchos::rcp(new ImplicitODEParameters<Scalar>(
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
240template<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
262template<class Scalar>
263Teuchos::RCP<Tempus::StepperState<Scalar> >
266{
267 Teuchos::RCP<Tempus::StepperState<Scalar> > stepperState =
268 rcp(new StepperState<Scalar>(this->getStepperType()));
269 return stepperState;
270}
271
272
273template<class Scalar>
275 Teuchos::FancyOStream &out,
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
302template<class Scalar>
303bool StepperBDF2<Scalar>::isValidSetup(Teuchos::FancyOStream & out) const
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
324template<class Scalar>
325Teuchos::RCP<const Teuchos::ParameterList>
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// ------------------------------------------------------------------------
336template<class Scalar>
337Teuchos::RCP<StepperBDF2<Scalar> >
339 const Teuchos::RCP<const Thyra::ModelEvaluator<Scalar> >& model,
340 Teuchos::RCP<Teuchos::ParameterList> pl)
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
Application Action for StepperBDF2.
void setStartUpStepper(std::string startupStepperType)
Set the stepper to use in first step.
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.
virtual bool isValidSetup(Teuchos::FancyOStream &out) const
virtual void setAppAction(Teuchos::RCP< StepperBDF2AppAction< Scalar > > appAction)
virtual void setModel(const Teuchos::RCP< const Thyra::ModelEvaluator< Scalar > > &appModel)
Set the model.
Teuchos::RCP< const Teuchos::ParameterList > getValidParameters() const
virtual void setInitialConditions(const Teuchos::RCP< SolutionHistory< Scalar > > &solutionHistory)
Set the initial conditions and make them consistent.
StepperBDF2()
Default constructor.
virtual Teuchos::RCP< Tempus::StepperState< Scalar > > getDefaultStepperState()
Get a default (initial) StepperState.
virtual void computeStartUp(const Teuchos::RCP< SolutionHistory< Scalar > > &solutionHistory)
Compute the first time step given the supplied startup stepper.
virtual void initialize()
Initialize during construction and after changing input parameters.
Thyra Base interface for implicit time steppers.
virtual void setInitialConditions(const Teuchos::RCP< SolutionHistory< Scalar > > &solutionHistory) override
Set the initial conditions and make them consistent.
virtual void describe(Teuchos::FancyOStream &out, const Teuchos::EVerbosityLevel verbLevel) const override
virtual void setModel(const Teuchos::RCP< const Thyra::ModelEvaluator< Scalar > > &appModel) override
Set the model.
Thyra Base interface for time steppers.
virtual void initialize()
Initialize after construction and changing input parameters.
virtual void describe(Teuchos::FancyOStream &out, const Teuchos::EVerbosityLevel verbLevel) const
Teuchos::RCP< StepperBDF2< Scalar > > createStepperBDF2(const Teuchos::RCP< const Thyra::ModelEvaluator< Scalar > > &model, Teuchos::RCP< Teuchos::ParameterList > pl)
Nonmember constructor - ModelEvaluator and ParameterList.