9 #ifndef Tempus_StepperNewmarkImplicitDForm_impl_hpp
10 #define Tempus_StepperNewmarkImplicitDForm_impl_hpp
19 template <
class Scalar>
24 #ifdef VERBOSE_DEBUG_OUTPUT
25 *out_ <<
"DEBUG: " << __PRETTY_FUNCTION__ <<
"\n";
28 Thyra::V_StVpStV(Teuchos::ptrFromRef(vPred), 1.0, v, dt * (1.0 - gamma_), a);
31 template <
class Scalar>
35 const Scalar dt)
const
37 #ifdef VERBOSE_DEBUG_OUTPUT
38 *out_ <<
"DEBUG: " << __PRETTY_FUNCTION__ <<
"\n";
41 Thyra::createMember<Scalar>(dPred.
space());
43 Scalar aConst = dt * dt / 2.0 * (1.0 - 2.0 * beta_);
44 Thyra::V_StVpStV(Teuchos::ptrFromRef(dPred), dt, v, aConst, a);
46 Thyra::Vp_V(Teuchos::ptrFromRef(dPred), d, 1.0);
49 template <
class Scalar>
54 #ifdef VERBOSE_DEBUG_OUTPUT
55 *out_ <<
"DEBUG: " << __PRETTY_FUNCTION__ <<
"\n";
58 Thyra::V_StVpStV(Teuchos::ptrFromRef(v), 1.0, vPred, dt * gamma_, a);
61 template <
class Scalar>
66 #ifdef VERBOSE_DEBUG_OUTPUT
67 *out_ <<
"DEBUG: " << __PRETTY_FUNCTION__ <<
"\n";
70 Thyra::V_StVpStV(Teuchos::ptrFromRef(d), 1.0, dPred, beta_ * dt * dt, a);
73 template <
class Scalar>
78 #ifdef VERBOSE_DEBUG_OUTPUT
79 *out_ <<
"DEBUG: " << __PRETTY_FUNCTION__ <<
"\n";
82 Scalar
const c = 1.0 / beta_ / dt / dt;
83 Thyra::V_StVpStV(Teuchos::ptrFromRef(a), c, d, -c, dPred);
86 template <
class Scalar>
89 if (schemeName_ !=
"User Defined") {
90 *out_ <<
"\nWARNING: schemeName != 'User Defined' (=" << schemeName_
92 <<
" Not setting beta, and leaving as beta = " << beta_ <<
"!\n";
99 *out_ <<
"\nWARNING: Running (implicit implementation of) Newmark "
100 <<
"Implicit a-Form Stepper with Beta = 0.0, which \n"
101 <<
"specifies an explicit scheme. Mass lumping is not possible, "
102 <<
"so this will be slow! To run explicit \n"
103 <<
"implementation of Newmark Implicit a-Form Stepper, please "
104 <<
"re-run with 'Stepper Type' = 'Newmark Explicit a-Form'.\n"
105 <<
"This stepper allows for mass lumping when called through "
106 <<
"Piro::TempusSolver.\n";
110 (beta_ > 1.0) || (beta_ < 0.0), std::logic_error,
111 "\nError in 'Newmark Implicit a-Form' stepper: invalid value of Beta = "
112 << beta_ <<
". Please select Beta >= 0 and <= 1. \n");
114 this->isInitialized_ =
false;
117 template <
class Scalar>
120 if (schemeName_ !=
"User Defined") {
121 *out_ <<
"\nWARNING: schemeName != 'User Defined' (=" << schemeName_
123 <<
" Not setting gamma, and leaving as gamma = " << gamma_ <<
"!\n";
130 (gamma_ > 1.0) || (gamma_ < 0.0), std::logic_error,
131 "\nError in 'Newmark Implicit a-Form' stepper: invalid value of Gamma ="
132 << gamma_ <<
". Please select Gamma >= 0 and <= 1. \n");
134 this->isInitialized_ =
false;
137 template <
class Scalar>
140 schemeName_ = schemeName;
142 if (schemeName_ ==
"Average Acceleration") {
146 else if (schemeName_ ==
"Linear Acceleration") {
150 else if (schemeName_ ==
"Central Difference") {
154 else if (schemeName_ ==
"User Defined") {
160 true, std::logic_error,
161 "\nError in Tempus::StepperNewmarkImplicitDForm! "
162 <<
"Invalid Scheme Name = " << schemeName_ <<
". \n"
163 <<
"Valid Scheme Names are: 'Average Acceleration', "
164 <<
"'Linear Acceleration', \n"
165 <<
"'Central Difference' and 'User Defined'.\n");
168 this->isInitialized_ =
false;
171 template <
class Scalar>
173 : out_(Teuchos::VerboseObjectBase::getDefaultOStream())
175 #ifdef VERBOSE_DEBUG_OUTPUT
176 *
out_ <<
"DEBUG: " << __PRETTY_FUNCTION__ <<
"\n";
191 template <
class Scalar>
195 bool useFSAL, std::string ICConsistency,
bool ICConsistencyCheck,
196 bool zeroInitialGuess, std::string schemeName, Scalar beta, Scalar gamma,
199 : out_(Teuchos::VerboseObjectBase::getDefaultOStream())
213 if (appModel != Teuchos::null) {
219 template <
class Scalar>
223 #ifdef VERBOSE_DEBUG_OUTPUT
224 *out_ <<
"DEBUG: " << __PRETTY_FUNCTION__ <<
"\n";
227 this->wrapperModel_ =
229 appModel,
"Newmark Implicit d-Form"));
232 std::logic_error,
"Error - Solver is not set!\n");
233 this->getSolver()->setModel(this->wrapperModel_);
235 this->isInitialized_ =
false;
238 template <
class Scalar>
242 if (appAction == Teuchos::null) {
244 stepperNewmarkImpAppAction_ =
248 stepperNewmarkImpAppAction_ = appAction;
251 this->isInitialized_ =
false;
254 template <
class Scalar>
258 #ifdef VERBOSE_DEBUG_OUTPUT
259 *out_ <<
"DEBUG: " << __PRETTY_FUNCTION__ <<
"\n";
261 this->checkInitialized();
265 TEMPUS_FUNC_TIME_MONITOR(
"Tempus::StepperNewmarkImplicitDForm::takeStep()");
268 solutionHistory->getNumStates() < 2, std::logic_error,
269 "Error - StepperNewmarkImplicitDForm<Scalar>::takeStep(...)\n"
270 <<
"Need at least two SolutionStates for NewmarkImplicitDForm.\n"
271 <<
" Number of States = " << solutionHistory->getNumStates()
272 <<
"\nTry setting in \"Solution History\" \"Storage Type\" = "
274 <<
" or \"Storage Type\" = \"Static\" and \"Storage Limit\" = "
277 auto thisStepper = Teuchos::rcpFromRef(*
this);
278 stepperNewmarkImpAppAction_->execute(
279 solutionHistory, thisStepper,
281 Scalar>::ACTION_LOCATION::BEGIN_STEP);
283 RCP<SolutionState<Scalar>> workingState =
284 solutionHistory->getWorkingState();
285 RCP<SolutionState<Scalar>> currentState =
286 solutionHistory->getCurrentState();
290 this->wrapperModel_);
293 RCP<const Thyra::VectorBase<Scalar>> d_old = currentState->getX();
294 RCP<Thyra::VectorBase<Scalar>> v_old = currentState->getXDot();
295 RCP<Thyra::VectorBase<Scalar>> a_old = currentState->getXDotDot();
299 RCP<Thyra::VectorBase<Scalar>> d_new = workingState->getX();
300 RCP<Thyra::VectorBase<Scalar>> v_new = workingState->getXDot();
301 RCP<Thyra::VectorBase<Scalar>> a_new = workingState->getXDotDot();
304 const Scalar time = currentState->getTime();
305 const Scalar dt = workingState->getTimeStep();
307 Scalar t = time + dt;
312 *out_ <<
"\n*** d_old ***\n";
313 RTOpPack::ConstSubVectorView<Scalar> dov;
314 d_old->acquireDetachedView(range, &dov);
315 auto doa = dov.values();
316 for (
auto i = 0; i < doa.size(); ++i) *out_ << doa[i] <<
" ";
317 *out_ <<
"\n*** d_old ***\n";
319 *out_ <<
"\n*** v_old ***\n";
320 RTOpPack::ConstSubVectorView<Scalar> vov;
321 v_old->acquireDetachedView(range, &vov);
322 auto voa = vov.values();
323 for (
auto i = 0; i < voa.size(); ++i) *out_ << voa[i] <<
" ";
324 *out_ <<
"\n*** v_old ***\n";
326 *out_ <<
"\n*** a_old ***\n";
327 RTOpPack::ConstSubVectorView<Scalar> aov;
328 a_old->acquireDetachedView(range, &aov);
329 auto aoa = aov.values();
330 for (
auto i = 0; i < aoa.size(); ++i) *out_ << aoa[i] <<
" ";
331 *out_ <<
"\n*** a_old ***\n";
335 RCP<Thyra::VectorBase<Scalar>> d_pred = Thyra::createMember(d_old->space());
336 RCP<Thyra::VectorBase<Scalar>> v_pred = Thyra::createMember(v_old->space());
339 predictDisplacement(*d_pred, *d_old, *v_old, *a_old, dt);
340 predictVelocity(*v_pred, *v_old, *a_old, dt);
343 *out_ <<
"\n*** d_pred ***\n";
344 RTOpPack::ConstSubVectorView<Scalar> dpv;
345 d_pred->acquireDetachedView(range, &dpv);
346 auto dpa = dpv.values();
347 for (
auto i = 0; i < dpa.size(); ++i) *out_ << dpa[i] <<
" ";
348 *out_ <<
"\n*** d_pred ***\n";
350 *out_ <<
"\n*** v_pred ***\n";
351 RTOpPack::ConstSubVectorView<Scalar> vpv;
352 v_pred->acquireDetachedView(range, &vpv);
353 auto vpa = vpv.values();
354 for (
auto i = 0; i < vpa.size(); ++i) *out_ << vpa[i] <<
" ";
355 *out_ <<
"\n*** v_pred ***\n";
359 wrapperModel->initializeNewmark(v_pred, d_pred, dt, t, beta_, gamma_);
362 RCP<Thyra::VectorBase<Scalar>> initial_guess =
363 Thyra::createMember(d_pred->space());
364 if ((time == solutionHistory->minTime()) &&
365 (this->initialGuess_ != Teuchos::null)) {
370 (initial_guess->space())->isCompatible(*this->initialGuess_->space());
372 is_compatible !=
true, std::logic_error,
373 "Error in Tempus::NemwarkImplicitDForm takeStep(): user-provided "
375 <<
"for Newton is not compatible with solution vector!\n");
376 Thyra::copy(*this->initialGuess_, initial_guess.ptr());
380 Thyra::copy(*d_pred, initial_guess.ptr());
383 stepperNewmarkImpAppAction_->execute(
384 solutionHistory, thisStepper,
386 Scalar>::ACTION_LOCATION::BEFORE_SOLVE);
390 (*(this->solver_)).solve(&*initial_guess);
392 workingState->setSolutionStatus(sStatus);
394 stepperNewmarkImpAppAction_->execute(
395 solutionHistory, thisStepper,
397 Scalar>::ACTION_LOCATION::AFTER_SOLVE);
401 Thyra::copy(*initial_guess, d_new.ptr());
404 correctAcceleration(*a_new, *d_pred, *d_new, dt);
405 correctVelocity(*v_new, *v_pred, *a_new, dt);
408 *out_ <<
"\n*** d_new ***\n";
409 RTOpPack::ConstSubVectorView<Scalar> dnv;
410 d_new->acquireDetachedView(range, &dnv);
411 auto dna = dnv.values();
412 for (
auto i = 0; i < dna.size(); ++i) *out_ << dna[i] <<
" ";
413 *out_ <<
"\n*** d_new ***\n";
415 *out_ <<
"\n*** v_new ***\n";
416 RTOpPack::ConstSubVectorView<Scalar> vnv;
417 v_new->acquireDetachedView(range, &vnv);
418 auto vna = vnv.values();
419 for (
auto i = 0; i < vna.size(); ++i) *out_ << vna[i] <<
" ";
420 *out_ <<
"\n*** v_new ***\n";
422 *out_ <<
"\n*** a_new ***\n";
423 RTOpPack::ConstSubVectorView<Scalar> anv;
424 a_new->acquireDetachedView(range, &anv);
425 auto ana = anv.values();
426 for (
auto i = 0; i < ana.size(); ++i) *out_ << ana[i] <<
" ";
427 *out_ <<
"\n*** a_new ***\n";
430 workingState->setOrder(this->getOrder());
431 workingState->computeNorms(currentState);
433 stepperNewmarkImpAppAction_->execute(
434 solutionHistory, thisStepper,
436 Scalar>::ACTION_LOCATION::END_STEP);
447 template <
class Scalar>
451 #ifdef VERBOSE_DEBUG_OUTPUT
452 *out_ <<
"DEBUG: " << __PRETTY_FUNCTION__ <<
"\n";
459 template <
class Scalar>
463 #ifdef VERBOSE_DEBUG_OUTPUT
464 *out_ <<
"DEBUG: " << __PRETTY_FUNCTION__ <<
"\n";
472 out <<
"--- StepperNewmarkImplicitDForm ---\n";
473 out <<
" schemeName_ = " << schemeName_ << std::endl;
474 out <<
" beta_ = " << beta_ << std::endl;
475 out <<
" gamma_ = " << gamma_ << std::endl;
476 out <<
"-----------------------------------" << std::endl;
479 template <
class Scalar>
484 bool isValidSetup =
true;
489 if (this->wrapperModel_->getAppModel() == Teuchos::null) {
490 isValidSetup =
false;
491 out <<
"The application ModelEvaluator is not set!\n";
494 if (this->wrapperModel_ == Teuchos::null) {
495 isValidSetup =
false;
496 out <<
"The wrapper ModelEvaluator is not set!\n";
499 if (this->solver_ == Teuchos::null) {
500 isValidSetup =
false;
501 out <<
"The solver is not set!\n";
507 template <
class Scalar>
511 #ifdef VERBOSE_DEBUG_OUTPUT
512 *out_ <<
"DEBUG: " << __PRETTY_FUNCTION__ <<
"\n";
514 auto pl = this->getValidParametersBasicImplicit();
516 auto newmarkPL = Teuchos::parameterList(
"Newmark Parameters");
517 newmarkPL->set<std::string>(
"Scheme Name", schemeName_);
518 newmarkPL->set<
double>(
"Beta", beta_);
519 newmarkPL->set<
double>(
"Gamma", gamma_);
520 pl->set(
"Newmark Parameters", *newmarkPL);
527 template <
class Scalar>
534 stepper->setStepperImplicitValues(pl);
536 if (pl != Teuchos::null) {
537 if (pl->
isSublist(
"Newmark Parameters")) {
538 auto newmarkPL = pl->
sublist(
"Newmark Parameters",
true);
539 std::string schemeName =
540 newmarkPL.
get<std::string>(
"Scheme Name",
"Average Acceleration");
541 stepper->setSchemeName(schemeName);
542 if (schemeName ==
"User Defined") {
543 stepper->setBeta(newmarkPL.get<
double>(
"Beta", 0.25));
544 stepper->setGamma(newmarkPL.get<
double>(
"Gamma", 0.5));
548 stepper->setSchemeName(
"Average Acceleration");
552 if (model != Teuchos::null) {
553 stepper->setModel(model);
554 stepper->initialize();
561 #endif // Tempus_StepperNewmarkImplicitDForm_impl_hpp
virtual RCP< const VectorSpaceBase< Scalar > > space() const =0
virtual void setSolver(Teuchos::RCP< Thyra::NonlinearSolverBase< Scalar > > solver) override
Set solver.
T & get(const std::string &name, T def_value)
Teuchos::RCP< StepperNewmarkImplicitDForm< Scalar > > createStepperNewmarkImplicitDForm(const Teuchos::RCP< const Thyra::ModelEvaluator< Scalar >> &model, Teuchos::RCP< Teuchos::ParameterList > pl)
Nonmember constructor - ModelEvaluator and ParameterList.
#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.
virtual void describe(Teuchos::FancyOStream &out, const Teuchos::EVerbosityLevel verbLevel) const override
A ModelEvaluator for residual evaluations given a state. This ModelEvaluator takes a state...
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
bool isSublist(const std::string &name) const
virtual void setDefaultSolver()
void setICConsistencyCheck(bool c)
virtual void setZeroInitialGuess(bool zIG)
SolutionHistory is basically a container of SolutionStates. SolutionHistory maintains a collection of...
basic_FancyOStream & setOutputToRootOnly(const int rootRank)
void validSecondOrderODE_DAE(const Teuchos::RCP< const Thyra::ModelEvaluator< Scalar > > &model)
ParameterList & sublist(const std::string &name, bool mustAlreadyExist=false, const std::string &docString="")
virtual void setUseFSAL(bool a)
void setStepperType(std::string s)
Set the stepper type.
void setICConsistency(std::string s)