Tempus Version of the Day
Time Integration
Loading...
Searching...
No Matches
Tempus_StepperNewmarkExplicitAForm_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_StepperNewmarkExplicitAForm_impl_hpp
10#define Tempus_StepperNewmarkExplicitAForm_impl_hpp
11
12#include "Thyra_VectorStdOps.hpp"
13
15
16
17//#define DEBUG_OUTPUT
18
19namespace Tempus {
20
21
22template<class Scalar>
27 const Scalar dt) const
28{
29 //vPred = v + dt*(1.0-gamma_)*a
30 Thyra::V_StVpStV(Teuchos::ptrFromRef(vPred), 1.0, v, dt*(1.0-gamma_), a);
31}
32
33template<class Scalar>
39 const Scalar dt) const
40{
41 Teuchos::RCP<const Thyra::VectorBase<Scalar> > tmp =
42 Thyra::createMember<Scalar>(dPred.space());
43 //dPred = dt*v + dt*dt/2.0*a
44 Scalar aConst = dt*dt/2.0;
45 Thyra::V_StVpStV(Teuchos::ptrFromRef(dPred), dt, v, aConst, a);
46 //dPred += d;
47 Thyra::Vp_V(Teuchos::ptrFromRef(dPred), d, 1.0);
48}
49
50template<class Scalar>
53 const Thyra::VectorBase<Scalar>& vPred,
55 const Scalar dt) const
56{
57 //v = vPred + dt*gamma_*a
58 Thyra::V_StVpStV(Teuchos::ptrFromRef(v), 1.0, vPred, dt*gamma_, a);
59}
60
61template<class Scalar>
63 : gammaDefault_(Scalar(0.5)), gamma_(Scalar(0.5))
64{
65 this->setStepperName( "Newmark Explicit a-Form");
66 this->setStepperType( "Newmark Explicit a-Form");
67 this->setUseFSAL( true);
68 this->setICConsistency( "Consistent");
69 this->setICConsistencyCheck( false);
70 this->setAppAction(Teuchos::null);
71}
72
73template<class Scalar>
75 const Teuchos::RCP<const Thyra::ModelEvaluator<Scalar> >& appModel,
76 bool useFSAL,
77 std::string ICConsistency,
78 bool ICConsistencyCheck,
79 Scalar gamma,
80 const Teuchos::RCP<StepperNewmarkExplicitAFormAppAction<Scalar> >& stepperAppAction)
81 : gammaDefault_(Scalar(0.5)), gamma_(Scalar(0.5))
82{
83 this->setStepperName( "Newmark Explicit a-Form");
84 this->setStepperType( "Newmark Explicit a-Form");
85 this->setUseFSAL( useFSAL);
86 this->setICConsistency( ICConsistency);
87 this->setICConsistencyCheck( ICConsistencyCheck);
88
89 this->setAppAction(stepperAppAction);
90 setGamma(gamma);
91
92 if (appModel != Teuchos::null) {
93
94 this->setModel(appModel);
95 this->initialize();
96 }
97}
98
99template<class Scalar>
101 const Teuchos::RCP<SolutionHistory<Scalar> >& solutionHistory)
102{
103 using Teuchos::RCP;
104
105 int numStates = solutionHistory->getNumStates();
106
107 TEUCHOS_TEST_FOR_EXCEPTION(numStates < 1, std::logic_error,
108 "Error - setInitialConditions() needs at least one SolutionState\n"
109 " to set the initial condition. Number of States = " << numStates);
110
111 if (numStates > 1) {
112 RCP<Teuchos::FancyOStream> out = this->getOStream();
113 Teuchos::OSTab ostab(out,1,
114 "StepperNewmarkExplicitAForm::setInitialConditions()");
115 *out << "Warning -- SolutionHistory has more than one state!\n"
116 << "Setting the initial conditions on the currentState.\n"<<std::endl;
117 }
118
119 RCP<SolutionState<Scalar> > initialState = solutionHistory->getCurrentState();
120 RCP<Thyra::VectorBase<Scalar> > x = initialState->getX();
121 RCP<Thyra::VectorBase<Scalar> > xDot = initialState->getXDot();
122
123 // If initialState has x and xDot set, treat them as the initial conditions.
124 // Otherwise use the x and xDot from getNominalValues() as the ICs.
125 TEUCHOS_TEST_FOR_EXCEPTION(
126 !((initialState->getX() != Teuchos::null &&
127 initialState->getXDot() != Teuchos::null) ||
128 (this->inArgs_.get_x() != Teuchos::null &&
129 this->inArgs_.get_x_dot() != Teuchos::null)), std::logic_error,
130 "Error - We need to set the initial conditions for x and xDot from\n"
131 " either initialState or appModel_->getNominalValues::InArgs\n"
132 " (but not from a mixture of the two).\n");
133
134 this->inArgs_ = this->appModel_->getNominalValues();
135 using Teuchos::rcp_const_cast;
136 // Use the x and xDot from getNominalValues() as the ICs.
137 if ( initialState->getX() == Teuchos::null ||
138 initialState->getXDot() == Teuchos::null ) {
139 TEUCHOS_TEST_FOR_EXCEPTION( (this->inArgs_.get_x() == Teuchos::null) ||
140 (this->inArgs_.get_x_dot() == Teuchos::null), std::logic_error,
141 "Error - setInitialConditions() needs the ICs from the SolutionHistory\n"
142 " or getNominalValues()!\n");
143 x =rcp_const_cast<Thyra::VectorBase<Scalar> >(this->inArgs_.get_x());
144 initialState->setX(x);
145 xDot =rcp_const_cast<Thyra::VectorBase<Scalar> >(this->inArgs_.get_x_dot());
146 initialState->setXDot(xDot);
147 } else {
148 this->inArgs_.set_x(x);
149 this->inArgs_.set_x_dot(xDot);
150 }
151
152 // Check if we need Stepper storage for xDotDot
153 if (initialState->getXDotDot() == Teuchos::null)
154 initialState->setXDotDot(initialState->getX()->clone_v());
155 else
156 this->setStepperXDotDot(initialState->getXDotDot());
157
158 // Perform IC Consistency
159 std::string icConsistency = this->getICConsistency();
160 if (icConsistency == "None") {
161 if (initialState->getXDotDot() == Teuchos::null) {
162 RCP<Teuchos::FancyOStream> out = this->getOStream();
163 Teuchos::OSTab ostab(out,1,"StepperForwardEuler::setInitialConditions()");
164 *out << "Warning -- Requested IC consistency of 'None' but\n"
165 << " initialState does not have an xDotDot.\n"
166 << " Setting a 'Consistent' xDotDot!\n" << std::endl;
167 auto p = Teuchos::rcp(new ExplicitODEParameters<Scalar>(0.0));
168 this->evaluateExplicitODE(initialState->getXDotDot(), x, xDot,
169 initialState->getTime(), p);
170
171 initialState->setIsSynced(true);
172 }
173 }
174 else if (icConsistency == "Zero")
175 Thyra::assign(initialState->getXDotDot().ptr(), Scalar(0.0));
176 else if (icConsistency == "App") {
177 auto xDotDot = Teuchos::rcp_const_cast<Thyra::VectorBase<Scalar> >(
178 this->inArgs_.get_x_dot_dot());
179 TEUCHOS_TEST_FOR_EXCEPTION(xDotDot == Teuchos::null, std::logic_error,
180 "Error - setInitialConditions() requested 'App' for IC consistency,\n"
181 " but 'App' returned a null pointer for xDotDot!\n");
182 Thyra::assign(initialState->getXDotDot().ptr(), *xDotDot);
183 }
184 else if (icConsistency == "Consistent") {
185 // Evaluate xDotDot = f(x,t).
186 auto p = Teuchos::rcp(new ExplicitODEParameters<Scalar>(0.0));
187 this->evaluateExplicitODE(initialState->getXDotDot(), x, xDot,
188 initialState->getTime(), p);
189
190 // At this point, x, xDot and xDotDot are sync'ed or consistent
191 // at the same time level for the initialState.
192 initialState->setIsSynced(true);
193 }
194 else {
195 TEUCHOS_TEST_FOR_EXCEPTION(true, std::logic_error,
196 "Error - setInitialConditions() invalid IC consistency, "
197 << icConsistency << ".\n");
198 }
199
200 // Test for consistency.
201 if (this->getICConsistencyCheck()) {
202 auto xDotDot = initialState->getXDotDot();
203 auto f = initialState->getX()->clone_v();
204 auto p = Teuchos::rcp(new ExplicitODEParameters<Scalar>(0.0));
205 this->evaluateExplicitODE(f, x, xDot, initialState->getTime(), p);
206 Thyra::Vp_StV(f.ptr(), Scalar(-1.0), *(xDotDot));
207 Scalar reldiff = Thyra::norm(*f);
208 Scalar normxDotDot = Thyra::norm(*xDotDot);
209 //The following logic is to prevent FPEs
210 Scalar eps = Scalar(100.0)*std::abs(Teuchos::ScalarTraits<Scalar>::eps());
211 if (normxDotDot > eps*reldiff) reldiff /= normxDotDot;
212
213 RCP<Teuchos::FancyOStream> out = this->getOStream();
214 Teuchos::OSTab ostab(out,1,"StepperNewmarkExplicitAForm::setInitialConditions()");
215 if (reldiff > eps) {
216 *out << "\n---------------------------------------------------\n"
217 << "Info -- Stepper = " << this->getStepperType() << "\n"
218 << " Initial condition PASSED consistency check!\n"
219 << " (||xDotDot-f(x,xDot,t)||/||x|| = " << reldiff << ") > "
220 << "(eps = " << eps << ")" << std::endl
221 << "---------------------------------------------------\n"<<std::endl;
222 } else {
223 *out << "\n---------------------------------------------------\n"
224 << "Info -- Stepper = " << this->getStepperType() << "\n"
225 << "Initial condition FAILED consistency check but continuing!\n"
226 << " (||xDotDot-f(x,xDot,t)||/||x|| = " << reldiff << ") > "
227 << "(eps = " << eps << ")" << std::endl
228 << " ||xDotDot-f(x,xDot,t)|| = " << Thyra::norm(*f) << std::endl
229 << " ||x|| = " << Thyra::norm(*x) << std::endl
230 << "---------------------------------------------------\n"<<std::endl;
231 }
232 }
233}
234
235
236template<class Scalar>
238 const Teuchos::RCP<SolutionHistory<Scalar> >& solutionHistory)
239{
240 this->checkInitialized();
241
242 using Teuchos::RCP;
243
244 TEMPUS_FUNC_TIME_MONITOR("Tempus::StepperNewmarkExplicitAForm::takeStep()");
245 {
246 TEUCHOS_TEST_FOR_EXCEPTION(solutionHistory->getNumStates() < 2,
247 std::logic_error,
248 "Error - StepperNewmarkExplicitAForm<Scalar>::takeStep(...)\n"
249 "Need at least two SolutionStates for NewmarkExplicitAForm.\n"
250 " Number of States = " << solutionHistory->getNumStates() << "\n"
251 "Try setting in \"Solution History\" \"Storage Type\" = \"Undo\"\n"
252 " or \"Storage Type\" = \"Static\" and \"Storage Limit\" = \"2\"\n");
253
254 auto thisStepper = Teuchos::rcpFromRef(*this);
255 stepperNewmarkExpAppAction_->execute(solutionHistory, thisStepper,
257
258 RCP<SolutionState<Scalar> > currentState=solutionHistory->getCurrentState();
259 RCP<SolutionState<Scalar> > workingState=solutionHistory->getWorkingState();
260
261 //Get values of d, v and a from previous step
262 RCP<const Thyra::VectorBase<Scalar> > d_old = currentState->getX();
263 RCP<const Thyra::VectorBase<Scalar> > v_old = currentState->getXDot();
264 RCP< Thyra::VectorBase<Scalar> > a_old = currentState->getXDotDot();
265
266 //Get dt and time
267 const Scalar dt = workingState->getTimeStep();
268 const Scalar time_old = currentState->getTime();
269
270 auto p = Teuchos::rcp(new ExplicitODEParameters<Scalar>(dt));
271 if (!(this->getUseFSAL()) || workingState->getNConsecutiveFailures() != 0) {
272 // Evaluate xDotDot = f(x, xDot, t).
273 this->evaluateExplicitODE(a_old, d_old, v_old, time_old, p);
274
275 // For UseFSAL=false, x and xDot sync'ed or consistent
276 // at the same time level for the currentState.
277 currentState->setIsSynced(true);
278 }
279
280 //New d, v and a to be computed here
281 RCP<Thyra::VectorBase<Scalar> > d_new = workingState->getX();
282 RCP<Thyra::VectorBase<Scalar> > v_new = workingState->getXDot();
283 RCP<Thyra::VectorBase<Scalar> > a_new = workingState->getXDotDot();
284
285 //compute displacement and velocity predictors
286 predictDisplacement(*d_new, *d_old, *v_old, *a_old, dt);
287 predictVelocity(*v_new, *v_old, *a_old, dt);
288
289 stepperNewmarkExpAppAction_->execute(solutionHistory, thisStepper,
291
292 // Evaluate xDotDot = f(x, xDot, t).
293 this->evaluateExplicitODE(a_new, d_new, v_new, time_old, p);
294
295 stepperNewmarkExpAppAction_->execute(solutionHistory, thisStepper,
297
298 // Set xDot in workingState to velocity corrector
299 correctVelocity(*v_new, *v_new, *a_new, dt);
300
301 if ( this->getUseFSAL() ) {
302 // Evaluate xDotDot = f(x, xDot, t).
303 const Scalar time_new = workingState->getTime();
304 this->evaluateExplicitODE(a_new, d_new, v_new, time_new, p);
305
306 // For UseFSAL=true, x, xDot and xDotxDot are now sync'ed or consistent
307 // for the workingState.
308 workingState->setIsSynced(true);
309 } else {
310 assign(a_new.ptr(), Teuchos::ScalarTraits<Scalar>::zero());
311 workingState->setIsSynced(false);
312 }
313
314 workingState->setSolutionStatus(Status::PASSED);
315 workingState->setOrder(this->getOrder());
316 workingState->computeNorms(currentState);
317
318 stepperNewmarkExpAppAction_->execute(solutionHistory, thisStepper,
320 }
321 return;
322}
323
324
331template<class Scalar>
332Teuchos::RCP<Tempus::StepperState<Scalar> > StepperNewmarkExplicitAForm<Scalar>::
334{
335 Teuchos::RCP<Tempus::StepperState<Scalar> > stepperState =
336 rcp(new StepperState<Scalar>(this->getStepperType()));
337 return stepperState;
338}
339
340
341template<class Scalar>
343 Teuchos::FancyOStream &out,
344 const Teuchos::EVerbosityLevel verbLevel) const
345{
346 out.setOutputToRootOnly(0);
347 out << std::endl;
348 Stepper<Scalar>::describe(out, verbLevel);
349 StepperExplicit<Scalar>::describe(out, verbLevel);
350
351 out << "--- StepperNewmarkExplicitAForm ---\n";
352 out << " gamma_ = " << gamma_ << std::endl;
353 out << "-----------------------------------" << std::endl;
354}
355
356
357template<class Scalar>
358bool StepperNewmarkExplicitAForm<Scalar>::isValidSetup(Teuchos::FancyOStream & out) const
359{
360 out.setOutputToRootOnly(0);
361 bool isValidSetup = true;
362
363 if ( !Stepper<Scalar>::isValidSetup(out) ) isValidSetup = false;
364 if ( !StepperExplicit<Scalar>::isValidSetup(out) ) isValidSetup = false;
365
366 return isValidSetup;
367}
368
369
370template<class Scalar>
371Teuchos::RCP<const Teuchos::ParameterList>
373{
374 auto pl = this->getValidParametersBasic();
375 pl->sublist("Newmark Explicit Parameters", false, "");
376 pl->sublist("Newmark Explicit Parameters", false, "").set("Gamma",
377 gamma_, "Newmark Explicit parameter");
378 return pl;
379}
380
381template<class Scalar>
383 Teuchos::RCP<StepperNewmarkExplicitAFormAppAction<Scalar> > appAction)
384{
385
386 if (appAction == Teuchos::null) {
387 // Create default appAction
388 stepperNewmarkExpAppAction_ =
390 } else {
391 stepperNewmarkExpAppAction_ = appAction;
392 }
393
394 this->isInitialized_ = false;
395}
396
397
398// Nonmember constructor - ModelEvaluator and ParameterList
399// ------------------------------------------------------------------------
400template<class Scalar>
401Teuchos::RCP<StepperNewmarkExplicitAForm<Scalar> >
403 const Teuchos::RCP<const Thyra::ModelEvaluator<Scalar> >& model,
404 Teuchos::RCP<Teuchos::ParameterList> pl)
405{
406 auto stepper = Teuchos::rcp(new StepperNewmarkExplicitAForm<Scalar>());
407 stepper->setStepperExplicitValues(pl);
408
409 if (pl != Teuchos::null) {
410 Scalar gamma = pl->sublist("Newmark Explicit Parameters")
411 .template get<double>("Gamma", 0.5);
412 stepper->setGamma(gamma);
413 }
414
415 if (model != Teuchos::null) {
416 stepper->setModel(model);
417 stepper->initialize();
418 }
419
420 return stepper;
421}
422
423
424} // namespace Tempus
425#endif // Tempus_StepperNewmarkExplicitAForm_impl_hpp
SolutionHistory is basically a container of SolutionStates. SolutionHistory maintains a collection of...
Thyra Base interface for implicit time steppers.
virtual void setModel(const Teuchos::RCP< const Thyra::ModelEvaluator< Scalar > > &appModel)
Set model.
virtual void describe(Teuchos::FancyOStream &out, const Teuchos::EVerbosityLevel verbLevel) const
Application Action for StepperNewmarkExplicitAForm.
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.
void correctVelocity(Thyra::VectorBase< Scalar > &v, const Thyra::VectorBase< Scalar > &vPred, const Thyra::VectorBase< Scalar > &a, const Scalar dt) const
virtual bool isValidSetup(Teuchos::FancyOStream &out) const
void predictVelocity(Thyra::VectorBase< Scalar > &vPred, const Thyra::VectorBase< Scalar > &v, const Thyra::VectorBase< Scalar > &a, const Scalar dt) const
virtual Teuchos::RCP< Tempus::StepperState< Scalar > > getDefaultStepperState()
Get a default (initial) StepperState.
virtual void setAppAction(Teuchos::RCP< StepperNewmarkExplicitAFormAppAction< Scalar > > appAction)
Teuchos::RCP< const Teuchos::ParameterList > getValidParameters() const
virtual void setInitialConditions(const Teuchos::RCP< SolutionHistory< Scalar > > &solutionHistory)
Set the initial conditions and make them consistent.
void predictDisplacement(Thyra::VectorBase< Scalar > &dPred, const Thyra::VectorBase< Scalar > &d, const Thyra::VectorBase< Scalar > &v, const Thyra::VectorBase< Scalar > &a, const Scalar dt) const
StepperState is a simple class to hold state information about the stepper.
Thyra Base interface for time steppers.
void setICConsistencyCheck(bool c)
void setStepperName(std::string s)
Set the stepper name.
virtual void initialize()
Initialize after construction and changing input parameters.
void setStepperType(std::string s)
Set the stepper type.
void setICConsistency(std::string s)
virtual void describe(Teuchos::FancyOStream &out, const Teuchos::EVerbosityLevel verbLevel) const
Teuchos::RCP< StepperNewmarkExplicitAForm< Scalar > > createStepperNewmarkExplicitAForm(const Teuchos::RCP< const Thyra::ModelEvaluator< Scalar > > &model, Teuchos::RCP< Teuchos::ParameterList > pl)
Nonmember constructor - ModelEvaluator and ParameterList.