Tempus Version of the Day
Time Integration
Loading...
Searching...
No Matches
Tempus_StepperDIRK_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_StepperDIRK_impl_hpp
10#define Tempus_StepperDIRK_impl_hpp
11
12#include "Thyra_VectorStdOps.hpp"
13
14#include "Tempus_WrapperModelEvaluatorBasic.hpp"
15
16
17namespace Tempus {
18
19
20template<class Scalar>
22{
23 this->setUseEmbedded( false);
24 this->setZeroInitialGuess( false);
25 this->setStageNumber(-1);
26
27 this->setAppAction(Teuchos::null);
28 this->setDefaultSolver();
29}
30
31
32template<class Scalar>
34 const Teuchos::RCP<const Thyra::ModelEvaluator<Scalar> >& appModel,
35 const Teuchos::RCP<Thyra::NonlinearSolverBase<Scalar> >& solver,
36 bool useFSAL,
37 std::string ICConsistency,
38 bool ICConsistencyCheck,
39 bool useEmbedded,
40 bool zeroInitialGuess,
41 const Teuchos::RCP<StepperRKAppAction<Scalar> >& stepperRKAppAction)
42{
43 this->setUseFSAL( useFSAL);
44 this->setICConsistency( ICConsistency);
45 this->setICConsistencyCheck( ICConsistencyCheck);
46 this->setUseEmbedded( useEmbedded);
47 this->setZeroInitialGuess( zeroInitialGuess);
48
49 this->setStageNumber(-1);
50 this->setErrorNorm();
51 this->setAppAction(stepperRKAppAction);
52 this->setSolver(solver);
53
54 if (appModel != Teuchos::null) {
55 this->setModel(appModel);
56 this->initialize();
57 }
58}
59
60
61template<class Scalar>
62Teuchos::RCP<Teuchos::ParameterList>
64{
65 auto pl = this->getValidParametersBasicImplicit();
66 pl->template set<bool>("Use Embedded", this->getUseEmbedded(),
67 "'Whether to use Embedded Stepper (if available) or not\n"
68 " 'true' - Stepper will compute embedded solution and is adaptive.\n"
69 " 'false' - Stepper is not embedded(adaptive).\n");
70 pl->template set<std::string>("Description", this->getDescription());
71 pl->template set<bool>("Reset Initial Guess", this->getResetInitialGuess());
72
73 return pl;
74}
75
76
77template<class Scalar>
79{
80 TEUCHOS_TEST_FOR_EXCEPTION(
81 this->tableau_ == Teuchos::null, std::logic_error,
82 "Error - Need to set the tableau, before calling "
83 "StepperDIRK::initialize()\n");
84
85 TEUCHOS_TEST_FOR_EXCEPTION(
86 this->wrapperModel_==Teuchos::null, std::logic_error,
87 "Error - Need to set the model, setModel(), before calling "
88 "StepperDIRK::initialize()\n");
89
91}
92
93
94template<class Scalar>
96 const Teuchos::RCP<const Thyra::ModelEvaluator<Scalar> >& appModel)
97{
99
100 // Set the stage vectors
101 const int numStages = this->tableau_->numStages();
102 stageXDot_.resize(numStages);
103 for (int i=0; i<numStages; ++i) {
104 stageXDot_[i] = Thyra::createMember(this->wrapperModel_->get_f_space());
105 assign(stageXDot_[i].ptr(), Teuchos::ScalarTraits<Scalar>::zero());
106 }
107 xTilde_ = Thyra::createMember(this->wrapperModel_->get_x_space());
108 assign(xTilde_.ptr(), Teuchos::ScalarTraits<Scalar>::zero());
109
110 this->setEmbeddedMemory();
111 this->setErrorNorm();
112
113 this->isInitialized_ = false;
114}
115
116
117template<class Scalar>
119{
120 if (this->getModel() == Teuchos::null)
121 return; // Embedded memory will be set when setModel() is called.
122
123 if (this->tableau_->isEmbedded() && this->getUseEmbedded()) {
124 this->ee_ = Thyra::createMember(this->wrapperModel_->get_f_space());
125 this->abs_u0 = Thyra::createMember(this->wrapperModel_->get_f_space());
126 this->abs_u = Thyra::createMember(this->wrapperModel_->get_f_space());
127 this->sc = Thyra::createMember(this->wrapperModel_->get_f_space());
128 } else {
129 this->ee_ = Teuchos::null;
130 this->abs_u0 = Teuchos::null;
131 this->abs_u = Teuchos::null;
132 this->sc = Teuchos::null;
133 }
134}
135
136
137template<class Scalar>
139 const Teuchos::RCP<SolutionHistory<Scalar> >& solutionHistory)
140{
141 this->setStepperXDot(stageXDot_.back());
143}
144
145
146template<class Scalar>
148 const Teuchos::RCP<SolutionHistory<Scalar> >& solutionHistory)
149{
150 this->checkInitialized();
151
152 using Teuchos::RCP;
153
154 TEMPUS_FUNC_TIME_MONITOR("Tempus::StepperDIRK::takeStep()");
155 {
156 TEUCHOS_TEST_FOR_EXCEPTION(solutionHistory->getNumStates() < 2,
157 std::logic_error,
158 "Error - StepperDIRK<Scalar>::takeStep(...)\n"
159 "Need at least two SolutionStates for DIRK.\n"
160 " Number of States = " << solutionHistory->getNumStates() << "\n"
161 "Try setting in \"Solution History\" \"Storage Type\" = \"Undo\"\n"
162 " or \"Storage Type\" = \"Static\" and \"Storage Limit\" = \"2\"\n");
163
164 RCP<SolutionState<Scalar> > currentState=solutionHistory->getCurrentState();
165 RCP<SolutionState<Scalar> > workingState=solutionHistory->getWorkingState();
166 const Scalar dt = workingState->getTimeStep();
167 const Scalar time = currentState->getTime();
168
169 const int numStages = this->tableau_->numStages();
170 Teuchos::SerialDenseMatrix<int,Scalar> A = this->tableau_->A();
171 Teuchos::SerialDenseVector<int,Scalar> b = this->tableau_->b();
172 Teuchos::SerialDenseVector<int,Scalar> c = this->tableau_->c();
173
174 // Reset non-zero initial guess.
175 if ( this->getResetInitialGuess() && (!this->getZeroInitialGuess()) )
176 Thyra::assign(workingState->getX().ptr(), *(currentState->getX()));
177
178 RCP<StepperDIRK<Scalar> > thisStepper = Teuchos::rcpFromRef(*this);
179 this->stepperRKAppAction_->execute(solutionHistory, thisStepper,
181
182 // Compute stage solutions
183 bool pass = true;
184 Thyra::SolveStatus<Scalar> sStatus;
185 for (int i=0; i < numStages; ++i) {
186 this->setStageNumber(i);
187
188 Thyra::assign(xTilde_.ptr(), *(currentState->getX()));
189 for (int j=0; j < i; ++j) {
190 if (A(i,j) != Teuchos::ScalarTraits<Scalar>::zero()) {
191 Thyra::Vp_StV(xTilde_.ptr(), dt*A(i,j), *(stageXDot_[j]));
192 }
193 }
194 this->setStepperXDot(stageXDot_[i]);
195
196 this->stepperRKAppAction_->execute(solutionHistory, thisStepper,
198
199 Scalar ts = time + c(i)*dt;
200
201 // Check if stageXDot_[i] is needed.
202 bool isNeeded = false;
203 for (int k=i+1; k<numStages; ++k) if (A(k,i) != 0.0) isNeeded = true;
204 if (b(i) != 0.0) isNeeded = true;
205 if (this->tableau_->isEmbedded() && this->getUseEmbedded() &&
206 this->tableau_->bstar()(i) != 0.0)
207 isNeeded = true;
208 if (isNeeded == false) {
209 assign(stageXDot_[i].ptr(), Teuchos::ScalarTraits<Scalar>::zero());
210 } else if (A(i,i) == Teuchos::ScalarTraits<Scalar>::zero()) {
211 // Explicit stage for the ImplicitODE_DAE
212 if (i == 0 && this->getUseFSAL() &&
213 workingState->getNConsecutiveFailures() == 0) {
214 // Reuse last evaluation for first step
215 RCP<Thyra::VectorBase<Scalar> > tmp = stageXDot_[0];
216 stageXDot_[0] = stageXDot_.back();
217 stageXDot_.back() = tmp;
218 this->setStepperXDot(stageXDot_[0]);
219 } else {
220 // Calculate explicit stage
221 typedef Thyra::ModelEvaluatorBase MEB;
222 MEB::InArgs<Scalar> inArgs = this->wrapperModel_->getInArgs();
223 MEB::OutArgs<Scalar> outArgs = this->wrapperModel_->getOutArgs();
224 inArgs.set_x(xTilde_);
225 if (inArgs.supports(MEB::IN_ARG_t)) inArgs.set_t(ts);
226 if (inArgs.supports(MEB::IN_ARG_x_dot))
227 inArgs.set_x_dot(Teuchos::null);
228 outArgs.set_f(stageXDot_[i]);
229
230 this->wrapperModel_->getAppModel()->evalModel(inArgs,outArgs);
231 }
232 } else {
233 // Implicit stage for the ImplicitODE_DAE
234 const Scalar alpha = 1.0/(dt*A(i,i));
235 const Scalar beta = 1.0;
236
237 // Setup TimeDerivative
238 Teuchos::RCP<TimeDerivative<Scalar> > timeDer =
239 Teuchos::rcp(new StepperDIRKTimeDerivative<Scalar>(
240 alpha,xTilde_.getConst()));
241
242 auto p = Teuchos::rcp(new ImplicitODEParameters<Scalar>(
243 timeDer, dt, alpha, beta, SOLVE_FOR_X, i));
244
245 this->stepperRKAppAction_->execute(solutionHistory, thisStepper,
247
248 sStatus = this->solveImplicitODE(workingState->getX(), stageXDot_[i], ts, p);
249
250 if (sStatus.solveStatus != Thyra::SOLVE_STATUS_CONVERGED) pass=false;
251
252 this->stepperRKAppAction_->execute(solutionHistory, thisStepper,
254
255 timeDer->compute(workingState->getX(), stageXDot_[i]);
256 }
257 this->stepperRKAppAction_->execute(solutionHistory, thisStepper,
259 this->stepperRKAppAction_->execute(solutionHistory, thisStepper,
261 }
262 // reset the stage number
263 this->setStageNumber(-1);
264
265 // Sum for solution: x_n = x_n-1 + Sum{ dt*b(i) * f(i) }
266 Thyra::assign((workingState->getX()).ptr(), *(currentState->getX()));
267 for (int i=0; i < numStages; ++i) {
268 if (b(i) != Teuchos::ScalarTraits<Scalar>::zero()) {
269 Thyra::Vp_StV((workingState->getX()).ptr(), dt*b(i), *(stageXDot_[i]));
270 }
271 }
272
273 if (this->tableau_->isEmbedded() && this->getUseEmbedded()) {
274 const Scalar tolRel = workingState->getTolRel();
275 const Scalar tolAbs = workingState->getTolAbs();
276
277 // update the tolerance
278 this->stepperErrorNormCalculator_->setRelativeTolerance(tolRel);
279 this->stepperErrorNormCalculator_->setAbsoluteTolerance(tolAbs);
280
281 // just compute the error weight vector
282 // (all that is needed is the error, and not the embedded solution)
283 Teuchos::SerialDenseVector<int,Scalar> errWght = b ;
284 errWght -= this->tableau_->bstar();
285
286 // compute local truncation error estimate: | u^{n+1} - \hat{u}^{n+1} |
287 // Sum for solution: ee_n = Sum{ (b(i) - bstar(i)) * dt*f(i) }
288 assign(this->ee_.ptr(), Teuchos::ScalarTraits<Scalar>::zero());
289 for (int i=0; i < numStages; ++i) {
290 if (errWght(i) != Teuchos::ScalarTraits<Scalar>::zero()) {
291 Thyra::Vp_StV(this->ee_.ptr(), dt*errWght(i), *(stageXDot_[i]));
292 }
293 }
294
295 Scalar err = this->stepperErrorNormCalculator_->computeWRMSNorm(currentState->getX(), workingState->getX(), this->ee_);
296 workingState->setErrorRel(err);
297
298 // test if step should be rejected
299 if (std::isinf(err) || std::isnan(err) || err > Teuchos::as<Scalar>(1.0))
300 pass = false;
301 }
302
303 if (pass) workingState->setSolutionStatus(Status::PASSED);
304 else workingState->setSolutionStatus(Status::FAILED);
305
306 workingState->setOrder(this->getOrder());
307 workingState->computeNorms(currentState);
308 this->stepperRKAppAction_->execute(solutionHistory, thisStepper,
310 }
311 return;
312}
313
320template<class Scalar>
321Teuchos::RCP<Tempus::StepperState<Scalar> >
324{
325 Teuchos::RCP<Tempus::StepperState<Scalar> > stepperState =
326 rcp(new StepperState<Scalar>(this->getStepperType()));
327 return stepperState;
328}
329
330
331template<class Scalar>
333 Teuchos::FancyOStream &out,
334 const Teuchos::EVerbosityLevel verbLevel) const
335{
336 out.setOutputToRootOnly(0);
337 out << std::endl;
338 Stepper<Scalar>::describe(out, verbLevel);
339 StepperImplicit<Scalar>::describe(out, verbLevel);
340
341 out << "--- StepperDIRK ---\n";
342 out << " tableau_ = " << this->tableau_ << std::endl;
343 if (this->tableau_ != Teuchos::null) this->tableau_->describe(out, verbLevel);
344 out << " stepperRKAppAction_ = " << this->stepperRKAppAction_ << std::endl;
345 out << " xTilde_ = " << xTilde_ << std::endl;
346 out << " stageXDot_.size() = " << stageXDot_.size() << std::endl;
347 const int numStages = stageXDot_.size();
348 for (int i=0; i<numStages; ++i)
349 out << " stageXDot_["<<i<<"] = " << stageXDot_[i] << std::endl;
350 out << " useEmbedded_ = "
351 << Teuchos::toString(this->useEmbedded_) << std::endl;
352 out << " ee_ = " << this->ee_ << std::endl;
353 out << " abs_u0 = " << this->abs_u0 << std::endl;
354 out << " abs_u = " << this->abs_u << std::endl;
355 out << " sc = " << this->sc << std::endl;
356 out << "-------------------" << std::endl;
357}
358
359
360template<class Scalar>
361bool StepperDIRK<Scalar>::isValidSetup(Teuchos::FancyOStream & out) const
362{
363 out.setOutputToRootOnly(0);
364 bool isValidSetup = true;
365
366 if ( !Stepper<Scalar>::isValidSetup(out) ) isValidSetup = false;
367 if ( !StepperImplicit<Scalar>::isValidSetup(out) ) isValidSetup = false;
368
369 if (this->tableau_ == Teuchos::null) {
370 isValidSetup = false;
371 out << "The tableau is not set!\n";
372 }
373
374 if (this->stepperRKAppAction_ == Teuchos::null) {
375 isValidSetup = false;
376 out << "The AppAction is not set!\n";
377 }
378
379 return isValidSetup;
380}
381
382
383template<class Scalar>
384Teuchos::RCP<const Teuchos::ParameterList>
386{
387 return this->getValidParametersBasicDIRK();
388}
389
390
391} // namespace Tempus
392#endif // Tempus_StepperDIRK_impl_hpp
SolutionHistory is basically a container of SolutionStates. SolutionHistory maintains a collection of...
Time-derivative interface for DIRK.
virtual void initialize() override
Initialize after construction and changing input parameters.
virtual void setupDefault()
Default setup for constructor.
virtual bool isValidSetup(Teuchos::FancyOStream &out) const override
virtual void takeStep(const Teuchos::RCP< SolutionHistory< Scalar > > &solutionHistory) override
Take the specified timestep, dt, and return true if successful.
Teuchos::RCP< Teuchos::ParameterList > getValidParametersBasicDIRK() const
virtual void describe(Teuchos::FancyOStream &out, const Teuchos::EVerbosityLevel verbLevel) const override
virtual void setInitialConditions(const Teuchos::RCP< SolutionHistory< Scalar > > &solutionHistory) override
Set the initial conditions and make them consistent.
virtual Teuchos::RCP< const Teuchos::ParameterList > getValidParameters() const override
virtual void setModel(const Teuchos::RCP< const Thyra::ModelEvaluator< Scalar > > &appModel) override
Set the model.
virtual Teuchos::RCP< Tempus::StepperState< Scalar > > getDefaultStepperState() override
Get a default (initial) StepperState.
virtual void setup(const Teuchos::RCP< const Thyra::ModelEvaluator< Scalar > > &wrapperModel, const Teuchos::RCP< Thyra::NonlinearSolverBase< Scalar > > &solver, bool useFSAL, std::string ICConsistency, bool ICConsistencyCheck, bool useEmbedded, bool zeroInitialGuess, const Teuchos::RCP< StepperRKAppAction< Scalar > > &stepperRKAppAction)
Setup for constructor.
virtual void setEmbeddedMemory() override
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.
Application Action for StepperRKBase.
StepperState is a simple class to hold state information about the stepper.
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
@ SOLVE_FOR_X
Solve for x and determine xDot from x.