9 #ifndef Tempus_StepperDIRK_impl_hpp
10 #define Tempus_StepperDIRK_impl_hpp
12 #include "Thyra_VectorStdOps.hpp"
14 #include "Tempus_WrapperModelEvaluatorBasic.hpp"
20 template<
class Scalar>
23 this->setUseEmbedded(
false);
24 this->setZeroInitialGuess(
false);
25 this->setStageNumber(-1);
27 this->setAppAction(Teuchos::null);
28 this->setDefaultSolver();
32 template<
class Scalar>
37 std::string ICConsistency,
38 bool ICConsistencyCheck,
40 bool zeroInitialGuess,
43 this->setUseFSAL( useFSAL);
44 this->setICConsistency( ICConsistency);
45 this->setICConsistencyCheck( ICConsistencyCheck);
46 this->setUseEmbedded( useEmbedded);
47 this->setZeroInitialGuess( zeroInitialGuess);
49 this->setStageNumber(-1);
51 this->setAppAction(stepperRKAppAction);
52 this->setSolver(solver);
54 if (appModel != Teuchos::null) {
55 this->setModel(appModel);
61 template<
class Scalar>
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());
77 template<
class Scalar>
81 this->tableau_ == Teuchos::null, std::logic_error,
82 "Error - Need to set the tableau, before calling "
83 "StepperDIRK::initialize()\n");
86 this->wrapperModel_==Teuchos::null, std::logic_error,
87 "Error - Need to set the model, setModel(), before calling "
88 "StepperDIRK::initialize()\n");
94 template<
class Scalar>
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());
107 xTilde_ = Thyra::createMember(this->wrapperModel_->get_x_space());
110 this->setEmbeddedMemory();
111 this->setErrorNorm();
113 this->isInitialized_ =
false;
117 template<
class Scalar>
120 if (this->getModel() == Teuchos::null)
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());
129 this->ee_ = Teuchos::null;
130 this->abs_u0 = Teuchos::null;
131 this->abs_u = Teuchos::null;
132 this->sc = Teuchos::null;
137 template<
class Scalar>
141 this->setStepperXDot(stageXDot_.back());
146 template<
class Scalar>
150 this->checkInitialized();
154 TEMPUS_FUNC_TIME_MONITOR(
"Tempus::StepperDIRK::takeStep()");
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");
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();
169 const int numStages = this->tableau_->numStages();
175 if ( this->getResetInitialGuess() && (!this->getZeroInitialGuess()) )
176 Thyra::assign(workingState->getX().ptr(), *(currentState->getX()));
178 RCP<StepperDIRK<Scalar> > thisStepper = Teuchos::rcpFromRef(*
this);
179 this->stepperRKAppAction_->execute(solutionHistory, thisStepper,
185 for (
int i=0; i < numStages; ++i) {
186 this->setStageNumber(i);
188 Thyra::assign(xTilde_.ptr(), *(currentState->getX()));
189 for (
int j=0; j < i; ++j) {
191 Thyra::Vp_StV(xTilde_.ptr(), dt*A(i,j), *(stageXDot_[j]));
194 this->setStepperXDot(stageXDot_[i]);
196 this->stepperRKAppAction_->execute(solutionHistory, thisStepper,
199 Scalar ts = time + c(i)*dt;
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)
208 if (isNeeded ==
false) {
212 if (i == 0 && this->getUseFSAL() &&
213 workingState->getNConsecutiveFailures() == 0) {
215 RCP<Thyra::VectorBase<Scalar> > tmp = stageXDot_[0];
216 stageXDot_[0] = stageXDot_.back();
217 stageXDot_.back() = tmp;
218 this->setStepperXDot(stageXDot_[0]);
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]);
230 this->wrapperModel_->getAppModel()->evalModel(inArgs,outArgs);
234 const Scalar alpha = 1.0/(dt*A(i,i));
235 const Scalar beta = 1.0;
240 alpha,xTilde_.getConst()));
245 this->stepperRKAppAction_->execute(solutionHistory, thisStepper,
248 sStatus = this->solveImplicitODE(workingState->getX(), stageXDot_[i], ts, p);
252 this->stepperRKAppAction_->execute(solutionHistory, thisStepper,
255 timeDer->compute(workingState->getX(), stageXDot_[i]);
257 this->stepperRKAppAction_->execute(solutionHistory, thisStepper,
259 this->stepperRKAppAction_->execute(solutionHistory, thisStepper,
263 this->setStageNumber(-1);
266 Thyra::assign((workingState->getX()).ptr(), *(currentState->getX()));
267 for (
int i=0; i < numStages; ++i) {
269 Thyra::Vp_StV((workingState->getX()).ptr(), dt*b(i), *(stageXDot_[i]));
273 if (this->tableau_->isEmbedded() && this->getUseEmbedded()) {
274 const Scalar tolRel = workingState->getTolRel();
275 const Scalar tolAbs = workingState->getTolAbs();
278 this->stepperErrorNormCalculator_->setRelativeTolerance(tolRel);
279 this->stepperErrorNormCalculator_->setAbsoluteTolerance(tolAbs);
284 errWght -= this->tableau_->bstar();
289 for (
int i=0; i < numStages; ++i) {
291 Thyra::Vp_StV(this->ee_.ptr(), dt*errWght(i), *(stageXDot_[i]));
295 Scalar err = this->stepperErrorNormCalculator_->computeWRMSNorm(currentState->getX(), workingState->getX(), this->ee_);
296 workingState->setErrorRel(err);
299 if (std::isinf(err) || std::isnan(err) || err > Teuchos::as<Scalar>(1.0))
306 workingState->setOrder(this->getOrder());
307 workingState->computeNorms(currentState);
308 this->stepperRKAppAction_->execute(solutionHistory, thisStepper,
320 template<
class Scalar>
331 template<
class Scalar>
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_ = "
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;
360 template<
class Scalar>
364 bool isValidSetup =
true;
369 if (this->tableau_ == Teuchos::null) {
370 isValidSetup =
false;
371 out <<
"The tableau is not set!\n";
374 if (this->stepperRKAppAction_ == Teuchos::null) {
375 isValidSetup =
false;
376 out <<
"The AppAction is not set!\n";
383 template<
class Scalar>
387 return this->getValidParametersBasicDIRK();
392 #endif // Tempus_StepperDIRK_impl_hpp
Teuchos::RCP< Teuchos::ParameterList > getValidParametersBasicDIRK() const
virtual void setModel(const Teuchos::RCP< const Thyra::ModelEvaluator< Scalar > > &appModel) override
Set the model.
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
virtual void initialize()
Initialize after construction and changing input parameters.
virtual void describe(Teuchos::FancyOStream &out, const Teuchos::EVerbosityLevel verbLevel) const override
virtual void describe(Teuchos::FancyOStream &out, const Teuchos::EVerbosityLevel verbLevel) const override
Thyra Base interface for time steppers.
virtual void setModel(const Teuchos::RCP< const Thyra::ModelEvaluator< Scalar > > &appModel) override
Set the model.
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
Application Action for StepperRKBase.
Time-derivative interface for DIRK.
virtual void setupDefault()
Default setup for constructor.
virtual Teuchos::RCP< const Teuchos::ParameterList > getValidParameters() const override
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 override
virtual void initialize() override
Initialize after construction and changing input parameters.
virtual void setInitialConditions(const Teuchos::RCP< SolutionHistory< Scalar > > &solutionHistory) override
Set the initial conditions and make them consistent.
virtual void setEmbeddedMemory() override
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 Teuchos::RCP< Tempus::StepperState< Scalar > > getDefaultStepperState() override
Get a default (initial) StepperState.
Solve for x and determine xDot from x.
virtual void takeStep(const Teuchos::RCP< SolutionHistory< Scalar > > &solutionHistory) override
Take the specified timestep, dt, and return true if successful.
std::string toString(const T &t)