9 #ifndef Tempus_StepperNewmarkExplicitAForm_impl_hpp
10 #define Tempus_StepperNewmarkExplicitAForm_impl_hpp
12 #include "Thyra_VectorStdOps.hpp"
20 template <
class Scalar>
26 Thyra::V_StVpStV(Teuchos::ptrFromRef(vPred), 1.0, v, dt * (1.0 - gamma_), a);
29 template <
class Scalar>
33 const Scalar dt)
const
36 Thyra::createMember<Scalar>(dPred.
space());
38 Scalar aConst = dt * dt / 2.0;
39 Thyra::V_StVpStV(Teuchos::ptrFromRef(dPred), dt, v, aConst, a);
41 Thyra::Vp_V(Teuchos::ptrFromRef(dPred), d, 1.0);
44 template <
class Scalar>
50 Thyra::V_StVpStV(Teuchos::ptrFromRef(v), 1.0, vPred, dt * gamma_, a);
53 template <
class Scalar>
55 : gammaDefault_(Scalar(0.5)), gamma_(Scalar(0.5))
65 template <
class Scalar>
68 bool useFSAL, std::string ICConsistency,
bool ICConsistencyCheck,
72 : gammaDefault_(Scalar(0.5)), gamma_(Scalar(0.5))
83 if (appModel != Teuchos::null) {
89 template <
class Scalar>
95 int numStates = solutionHistory->getNumStates();
98 numStates < 1, std::logic_error,
99 "Error - setInitialConditions() needs at least one SolutionState\n"
100 " to set the initial condition. Number of States = "
104 RCP<Teuchos::FancyOStream> out = this->getOStream();
106 "StepperNewmarkExplicitAForm::setInitialConditions()");
107 *out <<
"Warning -- SolutionHistory has more than one state!\n"
108 <<
"Setting the initial conditions on the currentState.\n"
112 RCP<SolutionState<Scalar> > initialState = solutionHistory->getCurrentState();
113 RCP<Thyra::VectorBase<Scalar> > x = initialState->getX();
114 RCP<Thyra::VectorBase<Scalar> > xDot = initialState->getXDot();
119 !((initialState->getX() != Teuchos::null &&
120 initialState->getXDot() != Teuchos::null) ||
121 (this->inArgs_.get_x() != Teuchos::null &&
122 this->inArgs_.get_x_dot() != Teuchos::null)),
124 "Error - We need to set the initial conditions for x and xDot from\n"
125 " either initialState or appModel_->getNominalValues::InArgs\n"
126 " (but not from a mixture of the two).\n");
128 this->inArgs_ = this->appModel_->getNominalValues();
129 using Teuchos::rcp_const_cast;
131 if (initialState->getX() == Teuchos::null ||
132 initialState->getXDot() == Teuchos::null) {
134 (this->inArgs_.get_x_dot() == Teuchos::null),
136 "Error - setInitialConditions() needs the ICs "
137 "from the SolutionHistory\n"
138 " or getNominalValues()!\n");
140 initialState->setX(x);
143 initialState->setXDot(xDot);
146 this->inArgs_.set_x(x);
147 this->inArgs_.set_x_dot(xDot);
151 if (initialState->getXDotDot() == Teuchos::null)
152 initialState->setXDotDot(initialState->getX()->clone_v());
154 this->setStepperXDotDot(initialState->getXDotDot());
157 std::string icConsistency = this->getICConsistency();
158 if (icConsistency ==
"None") {
159 if (initialState->getXDotDot() == Teuchos::null) {
160 RCP<Teuchos::FancyOStream> out = this->getOStream();
162 "StepperForwardEuler::setInitialConditions()");
163 *out <<
"Warning -- Requested IC consistency of 'None' but\n"
164 <<
" initialState does not have an xDotDot.\n"
165 <<
" Setting a 'Consistent' xDotDot!\n"
168 this->evaluateExplicitODE(initialState->getXDotDot(), x, xDot,
169 initialState->getTime(), p);
171 initialState->setIsSynced(
true);
174 else if (icConsistency ==
"Zero")
175 Thyra::assign(initialState->getXDotDot().ptr(), Scalar(0.0));
176 else if (icConsistency ==
"App") {
178 this->inArgs_.get_x_dot_dot());
180 xDotDot == Teuchos::null, std::logic_error,
181 "Error - setInitialConditions() requested 'App' for IC consistency,\n"
182 " but 'App' returned a null pointer for xDotDot!\n");
183 Thyra::assign(initialState->getXDotDot().ptr(), *xDotDot);
185 else if (icConsistency ==
"Consistent") {
188 this->evaluateExplicitODE(initialState->getXDotDot(), x, xDot,
189 initialState->getTime(), p);
193 initialState->setIsSynced(
true);
197 true, std::logic_error,
198 "Error - setInitialConditions() invalid IC consistency, "
199 << icConsistency <<
".\n");
203 if (this->getICConsistencyCheck()) {
204 auto xDotDot = initialState->getXDotDot();
205 auto f = initialState->getX()->clone_v();
207 this->evaluateExplicitODE(f, x, xDot, initialState->getTime(), p);
208 Thyra::Vp_StV(f.ptr(), Scalar(-1.0), *(xDotDot));
209 Scalar reldiff = Thyra::norm(*f);
210 Scalar normxDotDot = Thyra::norm(*xDotDot);
213 if (normxDotDot > eps * reldiff) reldiff /= normxDotDot;
215 RCP<Teuchos::FancyOStream> out = this->getOStream();
217 "StepperNewmarkExplicitAForm::setInitialConditions()");
219 *out <<
"\n---------------------------------------------------\n"
220 <<
"Info -- Stepper = " << this->getStepperType() <<
"\n"
221 <<
" Initial condition PASSED consistency check!\n"
222 <<
" (||xDotDot-f(x,xDot,t)||/||x|| = " << reldiff <<
") > "
223 <<
"(eps = " << eps <<
")" << std::endl
224 <<
"---------------------------------------------------\n"
228 *out <<
"\n---------------------------------------------------\n"
229 <<
"Info -- Stepper = " << this->getStepperType() <<
"\n"
230 <<
"Initial condition FAILED consistency check but continuing!\n"
231 <<
" (||xDotDot-f(x,xDot,t)||/||x|| = " << reldiff <<
") > "
232 <<
"(eps = " << eps <<
")" << std::endl
233 <<
" ||xDotDot-f(x,xDot,t)|| = " << Thyra::norm(*f) << std::endl
234 <<
" ||x|| = " << Thyra::norm(*x) << std::endl
235 <<
"---------------------------------------------------\n"
241 template <
class Scalar>
245 this->checkInitialized();
249 TEMPUS_FUNC_TIME_MONITOR(
"Tempus::StepperNewmarkExplicitAForm::takeStep()");
252 solutionHistory->getNumStates() < 2, std::logic_error,
253 "Error - StepperNewmarkExplicitAForm<Scalar>::takeStep(...)\n"
254 <<
"Need at least two SolutionStates for NewmarkExplicitAForm.\n"
255 <<
" Number of States = " << solutionHistory->getNumStates()
256 <<
"\nTry setting in \"Solution History\" \"Storage Type\" = "
258 <<
" or \"Storage Type\" = \"Static\" and \"Storage Limit\" = "
261 auto thisStepper = Teuchos::rcpFromRef(*
this);
262 stepperNewmarkExpAppAction_->execute(
263 solutionHistory, thisStepper,
265 Scalar>::ACTION_LOCATION::BEGIN_STEP);
267 RCP<SolutionState<Scalar> > currentState =
268 solutionHistory->getCurrentState();
269 RCP<SolutionState<Scalar> > workingState =
270 solutionHistory->getWorkingState();
273 RCP<const Thyra::VectorBase<Scalar> > d_old = currentState->getX();
274 RCP<const Thyra::VectorBase<Scalar> > v_old = currentState->getXDot();
275 RCP<Thyra::VectorBase<Scalar> > a_old = currentState->getXDotDot();
278 const Scalar dt = workingState->getTimeStep();
279 const Scalar time_old = currentState->getTime();
282 if (!(this->getUseFSAL()) || workingState->getNConsecutiveFailures() != 0) {
284 this->evaluateExplicitODE(a_old, d_old, v_old, time_old, p);
288 currentState->setIsSynced(
true);
292 RCP<Thyra::VectorBase<Scalar> > d_new = workingState->getX();
293 RCP<Thyra::VectorBase<Scalar> > v_new = workingState->getXDot();
294 RCP<Thyra::VectorBase<Scalar> > a_new = workingState->getXDotDot();
297 predictDisplacement(*d_new, *d_old, *v_old, *a_old, dt);
298 predictVelocity(*v_new, *v_old, *a_old, dt);
300 stepperNewmarkExpAppAction_->execute(
301 solutionHistory, thisStepper,
303 Scalar>::ACTION_LOCATION::BEFORE_EXPLICIT_EVAL);
306 this->evaluateExplicitODE(a_new, d_new, v_new, time_old, p);
308 stepperNewmarkExpAppAction_->execute(
309 solutionHistory, thisStepper,
311 Scalar>::ACTION_LOCATION::AFTER_EXPLICIT_EVAL);
314 correctVelocity(*v_new, *v_new, *a_new, dt);
316 if (this->getUseFSAL()) {
318 const Scalar time_new = workingState->getTime();
319 this->evaluateExplicitODE(a_new, d_new, v_new, time_new, p);
323 workingState->setIsSynced(
true);
327 workingState->setIsSynced(
false);
331 workingState->setOrder(this->getOrder());
332 workingState->computeNorms(currentState);
334 stepperNewmarkExpAppAction_->execute(
335 solutionHistory, thisStepper,
337 Scalar>::ACTION_LOCATION::END_STEP);
348 template <
class Scalar>
357 template <
class Scalar>
366 out <<
"--- StepperNewmarkExplicitAForm ---\n";
367 out <<
" gamma_ = " << gamma_ << std::endl;
368 out <<
"-----------------------------------" << std::endl;
371 template <
class Scalar>
376 bool isValidSetup =
true;
384 template <
class Scalar>
388 auto pl = this->getValidParametersBasic();
389 pl->sublist(
"Newmark Explicit Parameters",
false,
"");
390 pl->sublist(
"Newmark Explicit Parameters",
false,
"")
391 .set(
"Gamma", gamma_,
"Newmark Explicit parameter");
395 template <
class Scalar>
399 if (appAction == Teuchos::null) {
401 stepperNewmarkExpAppAction_ =
405 stepperNewmarkExpAppAction_ = appAction;
408 this->isInitialized_ =
false;
413 template <
class Scalar>
420 stepper->setStepperExplicitValues(pl);
422 if (pl != Teuchos::null) {
423 Scalar gamma = pl->
sublist(
"Newmark Explicit Parameters")
424 .template get<double>(
"Gamma", 0.5);
425 stepper->setGamma(gamma);
428 if (model != Teuchos::null) {
429 stepper->setModel(model);
430 stepper->initialize();
437 #endif // Tempus_StepperNewmarkExplicitAForm_impl_hpp
virtual RCP< const VectorSpaceBase< Scalar > > space() const =0
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
void setStepperName(std::string s)
Set the stepper name.
virtual void initialize()
Initialize after construction and changing input parameters.
Thyra Base interface for time steppers.
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)
virtual void describe(Teuchos::FancyOStream &out, const Teuchos::EVerbosityLevel verbLevel) const
void setICConsistencyCheck(bool c)
virtual void describe(Teuchos::FancyOStream &out, const Teuchos::EVerbosityLevel verbLevel) const
SolutionHistory is basically a container of SolutionStates. SolutionHistory maintains a collection of...
virtual void setModel(const Teuchos::RCP< const Thyra::ModelEvaluator< Scalar > > &appModel)
Set model.
basic_FancyOStream & setOutputToRootOnly(const int rootRank)
Teuchos::RCP< StepperNewmarkExplicitAForm< Scalar > > createStepperNewmarkExplicitAForm(const Teuchos::RCP< const Thyra::ModelEvaluator< Scalar > > &model, Teuchos::RCP< Teuchos::ParameterList > pl)
Nonmember constructor - ModelEvaluator and ParameterList.
ParameterList & sublist(const std::string &name, bool mustAlreadyExist=false, const std::string &docString="")
Thyra Base interface for implicit time steppers.
void setStepperType(std::string s)
Set the stepper type.
void setICConsistency(std::string s)