Tempus Version of the Day
Time Integration
Loading...
Searching...
No Matches
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
16namespace Tempus {
17
18
19template<class Scalar>
21{
22 this->setStepperName( "Trapezoidal Method");
23 this->setStepperType( "Trapezoidal Method");
24 this->setUseFSAL( true);
25 this->setICConsistency( "Consistent");
26 this->setICConsistencyCheck( false);
27 this->setZeroInitialGuess( false);
28
29 this->setAppAction(Teuchos::null);
30 this->setDefaultSolver();
31}
32
33
34template<class Scalar>
36 const Teuchos::RCP<const Thyra::ModelEvaluator<Scalar> >& appModel,
37 const Teuchos::RCP<Thyra::NonlinearSolverBase<Scalar> >& solver,
38 bool useFSAL,
39 std::string ICConsistency,
40 bool ICConsistencyCheck,
41 bool zeroInitialGuess,
42 const Teuchos::RCP<StepperTrapezoidalAppAction<Scalar> >& stepperTrapAppAction)
43{
44 this->setStepperName( "Trapezoidal Method");
45 this->setStepperType( "Trapezoidal Method");
46 this->setUseFSAL( useFSAL);
47 this->setICConsistency( ICConsistency);
48 this->setICConsistencyCheck( ICConsistencyCheck);
49 this->setZeroInitialGuess( zeroInitialGuess);
50
51 this->setAppAction(stepperTrapAppAction);
52 this->setSolver(solver);
53
54 if (appModel != Teuchos::null) {
55 this->setModel(appModel);
56 this->initialize();
57 }
58}
59
60
61template<class Scalar>
63 Teuchos::RCP<StepperTrapezoidalAppAction<Scalar> > appAction)
64{
65 if (appAction == Teuchos::null) {
66 // Create default appAction
67 stepperTrapAppAction_ =
68 Teuchos::rcp(new StepperTrapezoidalModifierDefault<Scalar>());
69 } else {
70 stepperTrapAppAction_ = appAction;
71 }
72 this->isInitialized_ = false;
73}
74
75
76template<class Scalar>
78 const Teuchos::RCP<SolutionHistory<Scalar> >& solutionHistory)
79{
80 using Teuchos::RCP;
81
82 RCP<SolutionState<Scalar> > initialState = solutionHistory->getCurrentState();
83
84 // Check if we need Stepper storage for xDot
85 if (initialState->getXDot() == Teuchos::null)
86 this->setStepperXDot(initialState->getX()->clone_v());
87 else
88 this->setStepperXDot(initialState->getXDot());
89
91
92 TEUCHOS_TEST_FOR_EXCEPTION( !(this->getUseFSAL()), std::logic_error,
93 "Error - The First-Same-As-Last (FSAL) principle is required\n"
94 " for the Trapezoidal Stepper (i.e., useFSAL=true)!\n");
95// There are at least two ways around this, but are not implemented.
96// - Do a solve for xDotOld, xDot_{n-1}, at each time step as for the
97// initial conditions. This is expensive since you would be doing
98// two solves every time step.
99// - Use evaluateExplicitODE to get xDot_{n-1} if the application
100// provides it. Explicit evaluations are cheaper but requires the
101// application to implement xDot = f(x,t).
102}
103
104
105template<class Scalar>
107 const Teuchos::RCP<SolutionHistory<Scalar> >& solutionHistory)
108{
109 this->checkInitialized();
110
111 using Teuchos::RCP;
112
113 TEMPUS_FUNC_TIME_MONITOR("Tempus::StepperTrapezoidal::takeStep()");
114 {
115 TEUCHOS_TEST_FOR_EXCEPTION(solutionHistory->getNumStates() < 2,
116 std::logic_error,
117 "Error - StepperTrapezoidal<Scalar>::takeStep(...)\n"
118 "Need at least two SolutionStates for Trapezoidal.\n"
119 " Number of States = " << solutionHistory->getNumStates() << "\n"
120 "Try setting in \"Solution History\" \"Storage Type\" = \"Undo\"\n"
121 " or \"Storage Type\" = \"Static\" and \"Storage Limit\" = \"2\"\n");
122 RCP<StepperTrapezoidal<Scalar> > thisStepper = Teuchos::rcpFromRef(*this);
123 stepperTrapAppAction_->execute(solutionHistory, thisStepper,
125
126 RCP<SolutionState<Scalar> > workingState=solutionHistory->getWorkingState();
127 RCP<SolutionState<Scalar> > currentState=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
142 Teuchos::RCP<TimeDerivative<Scalar> > timeDer =
143 Teuchos::rcp(new StepperTrapezoidalTimeDerivative<Scalar>(
144 alpha, xOld, xDotOld));
145
146 auto p = Teuchos::rcp(new ImplicitODEParameters<Scalar>(
147 timeDer, dt, alpha, beta));
148 stepperTrapAppAction_->execute(solutionHistory, thisStepper,
150
151 const Thyra::SolveStatus<Scalar> sStatus =
152 this->solveImplicitODE(x, xDot, time, p);
153 stepperTrapAppAction_->execute(solutionHistory, thisStepper,
155
156 if (workingState->getXDot() != Teuchos::null)
157 timeDer->compute(x, xDot);
158
159 workingState->setSolutionStatus(sStatus); // Converged --> pass.
160 workingState->setOrder(this->getOrder());
161 workingState->computeNorms(currentState);
162 stepperTrapAppAction_->execute(solutionHistory, thisStepper,
164 }
165 return;
166}
167
168
175template<class Scalar>
176Teuchos::RCP<Tempus::StepperState<Scalar> >
179{
180 Teuchos::RCP<Tempus::StepperState<Scalar> > stepperState =
181 rcp(new StepperState<Scalar>(this->getStepperType()));
182 return stepperState;
183}
184
185
186template<class Scalar>
188 Teuchos::FancyOStream &out,
189 const Teuchos::EVerbosityLevel verbLevel) const
190{
191 out.setOutputToRootOnly(0);
192 out << std::endl;
193 Stepper<Scalar>::describe(out, verbLevel);
194 StepperImplicit<Scalar>::describe(out, verbLevel);
195
196 out << "--- StepperTrapezoidal ---\n";
197 out << " stepperTrapAppAction_ = " << stepperTrapAppAction_ << std::endl;
198 out << "--------------------------" << std::endl;
199}
200
201
202template<class Scalar>
203bool StepperTrapezoidal<Scalar>::isValidSetup(Teuchos::FancyOStream & out) const
204{
205 out.setOutputToRootOnly(0);
206 bool isValidSetup = true;
207
208 if ( !Stepper<Scalar>::isValidSetup(out) ) isValidSetup = false;
209 if ( !StepperImplicit<Scalar>::isValidSetup(out) ) isValidSetup = false;
210
211 if (stepperTrapAppAction_ == Teuchos::null) {
212 isValidSetup = false;
213 out << "The Trapezoidal AppAction is not set!\n";
214 }
215
216 return isValidSetup;
217}
218
219
220// Nonmember constructor - ModelEvaluator and ParameterList
221// ------------------------------------------------------------------------
222template<class Scalar>
223Teuchos::RCP<StepperTrapezoidal<Scalar> >
225 const Teuchos::RCP<const Thyra::ModelEvaluator<Scalar> >& model,
226 Teuchos::RCP<Teuchos::ParameterList> pl)
227{
228 auto stepper = Teuchos::rcp(new StepperTrapezoidal<Scalar>());
229 stepper->setStepperImplicitValues(pl);
230
231 if (model != Teuchos::null) {
232 stepper->setModel(model);
233 stepper->initialize();
234 }
235
236 return stepper;
237}
238
239
240} // namespace Tempus
241#endif // Tempus_StepperTrapezoidal_impl_hpp
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
Application Action for StepperTrapezoidal.
virtual void takeStep(const Teuchos::RCP< SolutionHistory< Scalar > > &solutionHistory)
Take the specified timestep, dt, and return true if successful.
virtual void setInitialConditions(const Teuchos::RCP< SolutionHistory< Scalar > > &solutionHistory)
Set the initial conditions and make them consistent.
virtual bool isValidSetup(Teuchos::FancyOStream &out) const
virtual void setAppAction(Teuchos::RCP< StepperTrapezoidalAppAction< Scalar > > appAction)
virtual void describe(Teuchos::FancyOStream &out, const Teuchos::EVerbosityLevel verbLevel) const
virtual Teuchos::RCP< Tempus::StepperState< Scalar > > getDefaultStepperState()
Get a default (initial) StepperState.
Thyra Base interface for time steppers.
virtual void describe(Teuchos::FancyOStream &out, const Teuchos::EVerbosityLevel verbLevel) const
Teuchos::RCP< StepperTrapezoidal< Scalar > > createStepperTrapezoidal(const Teuchos::RCP< const Thyra::ModelEvaluator< Scalar > > &model, Teuchos::RCP< Teuchos::ParameterList > pl)
Nonmember constructor - ModelEvaluator and ParameterList.