9 #ifndef Tempus_StepperNewmarkImplicitAForm_impl_hpp
10 #define Tempus_StepperNewmarkImplicitAForm_impl_hpp
12 #include "Tempus_config.hpp"
14 #include "Teuchos_VerboseObjectParameterListHelpers.hpp"
15 #include "NOX_Thyra.H"
23 template<
class Scalar>
class StepperFactory;
26 template<
class Scalar>
29 const Thyra::VectorBase<Scalar>& v,
30 const Thyra::VectorBase<Scalar>& a,
31 const Scalar dt)
const
33 #ifdef VERBOSE_DEBUG_OUTPUT
34 *out_ <<
"DEBUG: " << __PRETTY_FUNCTION__ <<
"\n";
37 Thyra::V_StVpStV(Teuchos::ptrFromRef(vPred), 1.0, v, dt*(1.0-gamma_), a);
40 template<
class Scalar>
43 const Thyra::VectorBase<Scalar>& d,
44 const Thyra::VectorBase<Scalar>& v,
45 const Thyra::VectorBase<Scalar>& a,
46 const Scalar dt)
const
48 #ifdef VERBOSE_DEBUG_OUTPUT
49 *out_ <<
"DEBUG: " << __PRETTY_FUNCTION__ <<
"\n";
51 Teuchos::RCP<const Thyra::VectorBase<Scalar> > tmp =
52 Thyra::createMember<Scalar>(dPred.space());
54 Scalar aConst = dt*dt/2.0*(1.0-2.0*beta_);
55 Thyra::V_StVpStV(Teuchos::ptrFromRef(dPred), dt, v, aConst, a);
57 Thyra::Vp_V(Teuchos::ptrFromRef(dPred), d, 1.0);
60 template<
class Scalar>
63 const Thyra::VectorBase<Scalar>& vPred,
64 const Thyra::VectorBase<Scalar>& a,
65 const Scalar dt)
const
67 #ifdef VERBOSE_DEBUG_OUTPUT
68 *out_ <<
"DEBUG: " << __PRETTY_FUNCTION__ <<
"\n";
71 Thyra::V_StVpStV(Teuchos::ptrFromRef(v), 1.0, vPred, dt*gamma_, a);
74 template<
class Scalar>
77 const Thyra::VectorBase<Scalar>& dPred,
78 const Thyra::VectorBase<Scalar>& a,
79 const Scalar dt)
const
81 #ifdef VERBOSE_DEBUG_OUTPUT
82 *out_ <<
"DEBUG: " << __PRETTY_FUNCTION__ <<
"\n";
85 Thyra::V_StVpStV(Teuchos::ptrFromRef(d), 1.0, dPred, beta_*dt*dt, a);
89 template<
class Scalar>
92 if (schemeName_ !=
"User Defined") {
93 *out_ <<
"\nWARNING: schemeName != 'User Defined' (=" <<schemeName_<<
").\n"
94 <<
" Not setting beta, and leaving as beta = " << beta_ <<
"!\n";
101 *out_ <<
"\nWARNING: Running (implicit implementation of) Newmark "
102 <<
"Implicit a-Form Stepper with Beta = 0.0, which \n"
103 <<
"specifies an explicit scheme. Mass lumping is not possible, "
104 <<
"so this will be slow! To run explicit \n"
105 <<
"implementation of Newmark Implicit a-Form Stepper, please "
106 <<
"re-run with 'Stepper Type' = 'Newmark Explicit a-Form'.\n"
107 <<
"This stepper allows for mass lumping when called through "
108 <<
"Piro::TempusSolver.\n";
111 TEUCHOS_TEST_FOR_EXCEPTION( (beta_ > 1.0) || (beta_ < 0.0),
113 "\nError in 'Newmark Implicit a-Form' stepper: invalid value of Beta = "
114 << beta_ <<
". Please select Beta >= 0 and <= 1. \n");
118 template<
class Scalar>
121 if (schemeName_ !=
"User Defined") {
122 *out_ <<
"\nWARNING: schemeName != 'User Defined' (=" <<schemeName_<<
").\n"
123 <<
" Not setting gamma, and leaving as gamma = " << gamma_ <<
"!\n";
129 TEUCHOS_TEST_FOR_EXCEPTION( (gamma_ > 1.0) || (gamma_ < 0.0),
131 "\nError in 'Newmark Implicit a-Form' stepper: invalid value of Gamma ="
132 <<gamma_ <<
". Please select Gamma >= 0 and <= 1. \n");
136 template<
class Scalar>
138 std::string schemeName)
140 schemeName_ = schemeName;
142 if (schemeName_ ==
"Average Acceleration") {
143 beta_= 0.25; gamma_ = 0.5;
145 else if (schemeName_ ==
"Linear Acceleration") {
146 beta_= 0.25; gamma_ = 1.0/6.0;
148 else if (schemeName_ ==
"Central Difference") {
149 beta_= 0.0; gamma_ = 0.5;
151 else if (schemeName_ ==
"User Defined") {
152 beta_= 0.25; gamma_ = 0.5;
155 TEUCHOS_TEST_FOR_EXCEPTION(
true,
157 "\nError in Tempus::StepperNewmarkImplicitAForm! "
158 <<
"Invalid Scheme Name = " << schemeName_ <<
". \n"
159 <<
"Valid Scheme Names are: 'Average Acceleration', "
160 <<
"'Linear Acceleration', \n"
161 <<
"'Central Difference' and 'User Defined'.\n");
164 this->isInitialized_ =
false;
168 template<
class Scalar>
170 out_(
Teuchos::VerboseObjectBase::getDefaultOStream())
172 #ifdef VERBOSE_DEBUG_OUTPUT
173 *
out_ <<
"DEBUG: " << __PRETTY_FUNCTION__ <<
"\n";
188 template<
class Scalar>
190 const Teuchos::RCP<
const Thyra::ModelEvaluator<Scalar> >& appModel,
192 const Teuchos::RCP<Thyra::NonlinearSolverBase<Scalar> >& solver,
194 std::string ICConsistency,
195 bool ICConsistencyCheck,
196 bool zeroInitialGuess,
197 std::string schemeName,
200 : out_(
Teuchos::VerboseObjectBase::getDefaultOStream())
214 if (appModel != Teuchos::null) {
222 template<
class Scalar>
224 const Teuchos::RCP<
const Thyra::ModelEvaluator<Scalar> >& appModel)
226 #ifdef VERBOSE_DEBUG_OUTPUT
227 *out_ <<
"DEBUG: " << __PRETTY_FUNCTION__ <<
"\n";
232 "Newmark Implicit a-Form"));
233 this->wrapperModel_ = wrapperModel;
235 this->isInitialized_ =
false;
239 template<
class Scalar>
245 int numStates = solutionHistory->getNumStates();
247 TEUCHOS_TEST_FOR_EXCEPTION(numStates < 1, std::logic_error,
248 "Error - setInitialConditions() needs at least one SolutionState\n"
249 " to set the initial condition. Number of States = " << numStates);
252 RCP<Teuchos::FancyOStream> out = this->getOStream();
253 Teuchos::OSTab ostab(out,1,
"StepperNewmarkImplicitAForm::setInitialConditions()");
254 *out <<
"Warning -- SolutionHistory has more than one state!\n"
255 <<
"Setting the initial conditions on the currentState.\n"<<std::endl;
258 RCP<SolutionState<Scalar> > initialState = solutionHistory->getCurrentState();
259 RCP<Thyra::VectorBase<Scalar> > x = initialState->getX();
260 RCP<Thyra::VectorBase<Scalar> > xDot = initialState->getXDot();
262 auto inArgs = this->wrapperModel_->getNominalValues();
263 TEUCHOS_TEST_FOR_EXCEPTION(
264 !((x != Teuchos::null && xDot != Teuchos::null) ||
265 (inArgs.get_x() != Teuchos::null &&
266 inArgs.get_x_dot() != Teuchos::null)), std::logic_error,
267 "Error - We need to set the initial conditions for x and xDot from\n"
268 " either initialState or appModel_->getNominalValues::InArgs\n"
269 " (but not from a mixture of the two).\n");
272 if ( x == Teuchos::null || xDot == Teuchos::null ) {
273 using Teuchos::rcp_const_cast;
274 TEUCHOS_TEST_FOR_EXCEPTION( (inArgs.get_x() == Teuchos::null) ||
275 (inArgs.get_x_dot() == Teuchos::null), std::logic_error,
276 "Error - setInitialConditions() needs the ICs from the initialState\n"
277 " or getNominalValues()!\n");
278 x = rcp_const_cast<Thyra::VectorBase<Scalar> >(inArgs.get_x());
279 initialState->setX(x);
280 xDot = rcp_const_cast<Thyra::VectorBase<Scalar> >(inArgs.get_x_dot());
281 initialState->setXDot(xDot);
285 if (initialState->getXDotDot() == Teuchos::null)
286 initialState->setXDotDot(initialState->getX()->clone_v());
289 std::string icConsistency = this->getICConsistency();
290 if (icConsistency ==
"None") {
291 if (initialState->getXDotDot() == Teuchos::null) {
292 RCP<Teuchos::FancyOStream> out = this->getOStream();
293 Teuchos::OSTab ostab(out,1,
294 "StepperNewmarkImplicitAForm::setInitialConditions()");
295 *out <<
"Warning -- Requested IC consistency of 'None' but\n"
296 <<
" initialState does not have an xDot.\n"
297 <<
" Setting a 'Zero' xDot!\n" << std::endl;
299 Thyra::assign(this->getStepperXDotDot(initialState).ptr(), Scalar(0.0));
302 else if (icConsistency ==
"Zero")
303 Thyra::assign(this->getStepperXDotDot(initialState).ptr(), Scalar(0.0));
304 else if (icConsistency ==
"App") {
305 auto xDotDot = Teuchos::rcp_const_cast<Thyra::VectorBase<Scalar> >(
306 inArgs.get_x_dot_dot());
307 TEUCHOS_TEST_FOR_EXCEPTION(xDotDot == Teuchos::null, std::logic_error,
308 "Error - setInitialConditions() requested 'App' for IC consistency,\n"
309 " but 'App' returned a null pointer for xDotDot!\n");
310 Thyra::assign(this->getStepperXDotDot(initialState).ptr(), *xDotDot);
312 else if (icConsistency ==
"Consistent") {
314 const Scalar time = initialState->getTime();
315 auto xDotDot = this->getStepperXDotDot(initialState);
319 if (this->initialGuess_ != Teuchos::null) {
320 TEUCHOS_TEST_FOR_EXCEPTION(
321 !((xDotDot->space())->isCompatible(*this->initialGuess_->space())),
323 "Error - User-provided initial guess for Newton is not compatible\n"
324 " with solution vector!\n");
325 Thyra::copy(*this->initialGuess_, xDotDot.ptr());
328 Thyra::put_scalar(0.0, xDotDot.ptr());
332 Teuchos::rcp_dynamic_cast<WrapperModelEvaluatorSecondOrder<Scalar> >(
333 this->wrapperModel_);
335 wrapperModel->initializeNewmark(xDot, x, 0.0, time, beta_, gamma_);
336 const Thyra::SolveStatus<Scalar> sStatus = this->solveImplicitODE(xDotDot);
338 TEUCHOS_TEST_FOR_EXCEPTION(
339 sStatus.solveStatus != Thyra::SOLVE_STATUS_CONVERGED, std::logic_error,
340 "Error - Solver failed while determining the initial conditions.\n"
344 TEUCHOS_TEST_FOR_EXCEPTION(
true, std::logic_error,
345 "Error - setInitialConditions() invalid IC consistency, "
346 << icConsistency <<
".\n");
351 initialState->setIsSynced(
true);
354 if (this->getICConsistencyCheck()) {
355 auto f = initialState->getX()->clone_v();
356 auto xDotDot = this->getStepperXDotDot(initialState);
358 typedef Thyra::ModelEvaluatorBase MEB;
359 MEB::InArgs<Scalar> appInArgs =
360 this->wrapperModel_->getAppModel()->createInArgs();
361 MEB::OutArgs<Scalar> appOutArgs =
362 this->wrapperModel_->getAppModel()->createOutArgs();
364 appInArgs.set_x (x );
365 appInArgs.set_x_dot (xDot );
366 appInArgs.set_x_dot_dot(xDotDot);
368 appOutArgs.set_f(appOutArgs.get_f());
370 appInArgs.set_W_x_dot_dot_coeff(Scalar(0.0));
371 appInArgs.set_alpha (Scalar(0.0));
372 appInArgs.set_beta (Scalar(0.0));
374 appInArgs.set_t (initialState->getTime() );
376 this->wrapperModel_->getAppModel()->evalModel(appInArgs, appOutArgs);
378 Scalar reldiff = Thyra::norm(*f);
379 Scalar normx = Thyra::norm(*x);
380 Scalar eps = Scalar(100.0)*std::abs(Teuchos::ScalarTraits<Scalar>::eps());
381 if (normx > eps*reldiff) reldiff /= normx;
384 RCP<Teuchos::FancyOStream> out = this->getOStream();
385 Teuchos::OSTab ostab(out,1,
386 "StepperNewmarkImplicitAForm::setInitialConditions()");
387 *out <<
"Warning -- Failed consistency check but continuing!\n"
388 <<
" ||f(x,xDot,xDotDot,t)||/||x|| > eps" << std::endl
389 <<
" ||f(x,xDot,xDotDot,t)|| = " << Thyra::norm(*f)<< std::endl
390 <<
" ||x|| = " << Thyra::norm(*x)<< std::endl
391 <<
" ||f(x,xDot,xDotDot,t)||/||x|| = " << reldiff << std::endl
392 <<
" eps = " << eps << std::endl;
396 if (!(this->getUseFSAL())) {
397 RCP<Teuchos::FancyOStream> out = this->getOStream();
398 Teuchos::OSTab ostab(out,1,
399 "StepperNewmarkImplicitAForm::setInitialConditions()");
400 *out <<
"\nWarning -- The First-Step-As-Last (FSAL) principle is "
401 <<
"part of the Newmark Implicit A-Form. The default is to "
402 <<
"set useFSAL=true, and useFSAL=false will be ignored." << std::endl;
407 template<
class Scalar>
411 #ifdef VERBOSE_DEBUG_OUTPUT
412 *out_ <<
"DEBUG: " << __PRETTY_FUNCTION__ <<
"\n";
414 this->checkInitialized();
418 TEMPUS_FUNC_TIME_MONITOR(
"Tempus::StepperNewmarkImplicitAForm::takeStep()");
420 TEUCHOS_TEST_FOR_EXCEPTION(solutionHistory->getNumStates() < 2,
422 "Error - StepperNewmarkImplicitAForm<Scalar>::takeStep(...)\n"
423 "Need at least two SolutionStates for NewmarkImplicitAForm.\n"
424 " Number of States = " << solutionHistory->getNumStates() <<
"\n"
425 "Try setting in \"Solution History\" \"Storage Type\" = \"Undo\"\n"
426 " or \"Storage Type\" = \"Static\" and \"Storage Limit\" = \"2\"\n");
428 RCP<SolutionState<Scalar> > workingState=solutionHistory->getWorkingState();
429 RCP<SolutionState<Scalar> > currentState=solutionHistory->getCurrentState();
432 Teuchos::rcp_dynamic_cast<WrapperModelEvaluatorSecondOrder<Scalar> >(
433 this->wrapperModel_);
436 RCP<const Thyra::VectorBase<Scalar> > d_old = currentState->getX();
437 RCP<const Thyra::VectorBase<Scalar> > v_old = currentState->getXDot();
438 RCP< Thyra::VectorBase<Scalar> > a_old = currentState->getXDotDot();
441 RCP<Thyra::VectorBase<Scalar> > d_new = workingState->getX();
442 RCP<Thyra::VectorBase<Scalar> > v_new = workingState->getXDot();
443 RCP<Thyra::VectorBase<Scalar> > a_new = workingState->getXDotDot();
446 const Scalar time = currentState->getTime();
447 const Scalar dt = workingState->getTimeStep();
452 if (!(this->getUseFSAL()) && workingState->getNConsecutiveFailures() == 0) {
453 wrapperModel->initializeNewmark(v_old, d_old, dt, time,
454 Scalar(0.0), Scalar(0.0));
455 const Thyra::SolveStatus<Scalar> sStatus = this->solveImplicitODE(a_old);
457 workingState->setSolutionStatus(sStatus);
461 predictDisplacement(*d_new, *d_old, *v_old, *a_old, dt);
462 predictVelocity(*v_new, *v_old, *a_old, dt);
465 wrapperModel->initializeNewmark(v_new,d_new,dt,t,beta_,gamma_);
468 const Thyra::SolveStatus<Scalar> sStatus = this->solveImplicitODE(a_new);
471 correctVelocity(*v_new, *v_new, *a_new, dt);
472 correctDisplacement(*d_new, *d_new, *a_new, dt);
474 workingState->setSolutionStatus(sStatus);
475 workingState->setOrder(this->getOrder());
476 workingState->computeNorms(currentState);
489 template<
class Scalar>
490 Teuchos::RCP<Tempus::StepperState<Scalar> >
494 #ifdef VERBOSE_DEBUG_OUTPUT
495 *out_ <<
"DEBUG: " << __PRETTY_FUNCTION__ <<
"\n";
497 Teuchos::RCP<Tempus::StepperState<Scalar> > stepperState =
503 template<
class Scalar>
505 Teuchos::FancyOStream &out,
506 const Teuchos::EVerbosityLevel verbLevel)
const
508 #ifdef VERBOSE_DEBUG_OUTPUT
509 *out_ <<
"DEBUG: " << __PRETTY_FUNCTION__ <<
"\n";
516 out <<
"--- StepperNewmarkImplicitAForm ---\n";
517 out <<
" schemeName_ = " << schemeName_ << std::endl;
518 out <<
" beta_ = " << beta_ << std::endl;
519 out <<
" gamma_ = " << gamma_ << std::endl;
520 out <<
"-----------------------------------" << std::endl;
524 template<
class Scalar>
527 bool isValidSetup =
true;
532 if (this->wrapperModel_->getAppModel() == Teuchos::null) {
533 isValidSetup =
false;
534 out <<
"The application ModelEvaluator is not set!\n";
537 if (this->wrapperModel_ == Teuchos::null) {
538 isValidSetup =
false;
539 out <<
"The wrapper ModelEvaluator is not set!\n";
542 if (this->solver_ == Teuchos::null) {
543 isValidSetup =
false;
544 out <<
"The solver is not set!\n";
551 template<
class Scalar>
552 Teuchos::RCP<const Teuchos::ParameterList>
555 #ifdef VERBOSE_DEBUG_OUTPUT
556 *out_ <<
"DEBUG: " << __PRETTY_FUNCTION__ <<
"\n";
558 Teuchos::RCP<Teuchos::ParameterList> pl = Teuchos::parameterList();
560 pl->set<std::string>(
"Scheme Name",
"Average Acceleration");
561 pl->set<
double> (
"Beta" , 0.25);
562 pl->set<
double> (
"Gamma", 0.5 );
563 pl->set<
bool> (
"Use FSAL", this->getUseFSALDefault());
564 pl->set<std::string>(
"Initial Condition Consistency",
565 this->getICConsistencyDefault());
566 pl->set<std::string>(
"Solver Name",
"Default Solver");
567 pl->set<
bool> (
"Zero Initial Guess",
false);
569 pl->set(
"Default Solver", *solverPL);
SolutionHistory is basically a container of SolutionStates. SolutionHistory maintains a collection of...
virtual void describe(Teuchos::FancyOStream &out, const Teuchos::EVerbosityLevel verbLevel) const
virtual void setDefaultSolver()
virtual void setZeroInitialGuess(bool zIG)
Set parameter so that the initial guess is set to zero (=True) or use last timestep (=False).
virtual void setSolver(Teuchos::RCP< Thyra::NonlinearSolverBase< Scalar > > solver)
Set solver.
StepperObserver class for Stepper class.
StepperState is a simple class to hold state information about the stepper.
Thyra Base interface for time steppers.
virtual bool getICConsistencyCheckDefault() const
void setICConsistencyCheck(bool c)
virtual void initialize()
Initialize after construction and changing input parameters.
void setStepperType(std::string s)
void setICConsistency(std::string s)
virtual void describe(Teuchos::FancyOStream &out, const Teuchos::EVerbosityLevel verbLevel) const
A ModelEvaluator for residual evaluations given a state. This ModelEvaluator takes a state,...
const std::string toString(const Status status)
Convert Status to string.
void getValidParametersBasic(Teuchos::RCP< Teuchos::ParameterList > pl, std::string stepperType)
Provide basic parameters to Steppers.
void validSecondOrderODE_DAE(const Teuchos::RCP< const Thyra::ModelEvaluator< Scalar > > &model)
Validate ME supports 2nd order implicit ODE/DAE evaluation, f(xdotdot,xdot,x,t) [= 0].
Teuchos::RCP< Teuchos::ParameterList > defaultSolverParameters()
Returns the default solver ParameterList for implicit Steppers.