337 const Teuchos::RCP<SolutionHistory<Scalar> >& solutionHistory)
339#ifdef VERBOSE_DEBUG_OUTPUT
340 *out_ <<
"DEBUG: " << __PRETTY_FUNCTION__ <<
"\n";
342 this->checkInitialized();
346 TEMPUS_FUNC_TIME_MONITOR(
"Tempus::StepperHHTAlpha::takeStep()");
348 TEUCHOS_TEST_FOR_EXCEPTION(solutionHistory->getNumStates() < 2,
350 "Error - StepperHHTAlpha<Scalar>::takeStep(...)\n"
351 "Need at least two SolutionStates for HHTAlpha.\n"
352 " Number of States = " << solutionHistory->getNumStates() <<
"\n"
353 "Try setting in \"Solution History\" \"Storage Type\" = \"Undo\"\n"
354 " or \"Storage Type\" = \"Static\" and \"Storage Limit\" = \"2\"\n");
356 RCP<StepperHHTAlpha<Scalar> > thisStepper = Teuchos::rcpFromRef(*
this);
357 stepperHHTAlphaAppAction_->execute(solutionHistory, thisStepper,
360 RCP<SolutionState<Scalar> > workingState=solutionHistory->getWorkingState();
361 RCP<SolutionState<Scalar> > currentState=solutionHistory->getCurrentState();
363 Teuchos::RCP<WrapperModelEvaluatorSecondOrder<Scalar> > wrapperModel =
364 Teuchos::rcp_dynamic_cast<WrapperModelEvaluatorSecondOrder<Scalar> >(
365 this->wrapperModel_);
368 RCP<const Thyra::VectorBase<Scalar> > d_old = currentState->getX();
369 RCP<const Thyra::VectorBase<Scalar> > v_old = currentState->getXDot();
370 RCP<Thyra::VectorBase<Scalar> > a_old = currentState->getXDotDot();
375 *out_ <<
"IKT d_old = " << Thyra::max(*d_old) <<
"\n";
376 *out_ <<
"IKT v_old = " << Thyra::max(*v_old) <<
"\n";
381 RCP<Thyra::VectorBase<Scalar> > d_new = workingState->getX();
382 RCP<Thyra::VectorBase<Scalar> > v_new = workingState->getXDot();
383 RCP<Thyra::VectorBase<Scalar> > a_new = workingState->getXDotDot();
386 const Scalar time = currentState->getTime();
387 const Scalar dt = workingState->getTimeStep();
393 if (time == solutionHistory->minTime()) {
394 RCP<Thyra::VectorBase<Scalar> > d_init = Thyra::createMember(d_old->space());
395 RCP<Thyra::VectorBase<Scalar> > v_init = Thyra::createMember(v_old->space());
396 RCP<Thyra::VectorBase<Scalar> > a_init = Thyra::createMember(a_old->space());
397 Thyra::copy(*d_old, d_init.ptr());
398 Thyra::copy(*v_old, v_init.ptr());
399 if (this->initialGuess_ != Teuchos::null) {
401 bool is_compatible = (a_init->space())->isCompatible(*this->initialGuess_->space());
402 TEUCHOS_TEST_FOR_EXCEPTION(
403 is_compatible !=
true, std::logic_error,
404 "Error in Tempus::NemwarkImplicitAForm takeStep(): user-provided initial guess'!\n"
405 <<
"for Newton is not compatible with solution vector!\n");
406 Thyra::copy(*this->initialGuess_, a_init.ptr());
409 Thyra::put_scalar(0.0, a_init.ptr());
411 wrapperModel->initializeNewmark(v_init,d_init,0.0,time,beta_,gamma_);
412 const Thyra::SolveStatus<Scalar> sStatus=(*(this->solver_)).solve(&*a_init);
414 workingState->setSolutionStatus(sStatus);
415 Thyra::copy(*a_init, a_old.ptr());
419 *out_ <<
"IKT a_old = " << Thyra::max(*a_old) <<
"\n";
423 RCP<Thyra::VectorBase<Scalar> > d_pred =Thyra::createMember(d_old->space());
424 RCP<Thyra::VectorBase<Scalar> > v_pred =Thyra::createMember(v_old->space());
427 predictDisplacement(*d_pred, *d_old, *v_old, *a_old, dt);
428 predictVelocity(*v_pred, *v_old, *a_old, dt);
431 predictDisplacement_alpha_f(*d_pred, *d_old);
432 predictVelocity_alpha_f(*v_pred, *v_old);
435 wrapperModel->initializeNewmark(v_pred,d_pred,dt,t,beta_,gamma_);
438 stepperHHTAlphaAppAction_->execute(solutionHistory, thisStepper,
442 const Thyra::SolveStatus<Scalar> sStatus = (*(this->solver_)).solve(&*a_new);
444 stepperHHTAlphaAppAction_->execute(solutionHistory, thisStepper,
448 correctAcceleration(*a_new, *a_old);
451 correctVelocity(*v_new, *v_pred, *a_new, dt);
452 correctDisplacement(*d_new, *d_pred, *a_new, dt);
454 workingState->setSolutionStatus(sStatus);
455 workingState->setOrder(this->getOrder());
456 workingState->computeNorms(currentState);
458 stepperHHTAlphaAppAction_->execute(solutionHistory, thisStepper,
568 Teuchos::RCP<Teuchos::ParameterList> pl)
570 auto stepper = Teuchos::rcp(
new StepperHHTAlpha<Scalar>());
571 stepper->setStepperImplicitValues(pl);
573 if (pl != Teuchos::null) {
574 if (pl->isSublist(
"HHT-Alpha Parameters")) {
575 auto hhtalphaPL = pl->sublist(
"HHT-Alpha Parameters",
true);
576 std::string schemeName =
577 hhtalphaPL.get<std::string>(
"Scheme Name",
"Newmark Beta Average Acceleration");
578 stepper->setSchemeName(schemeName);
579 if (schemeName ==
"Newmark Beta User Defined") {
580 stepper->setBeta (hhtalphaPL.get<
double>(
"Beta", 0.25));
581 stepper->setGamma(hhtalphaPL.get<
double>(
"Gamma", 0.5 ));
583 stepper->setAlphaF(hhtalphaPL.get<
double>(
"Alpha_f", 0.0));
584 stepper->setAlphaM(hhtalphaPL.get<
double>(
"Alpha_m", 0.0));
586 stepper->setSchemeName(
"Newmark Beta Average Acceleration");
587 stepper->setAlphaF(0.0);
588 stepper->setAlphaM(0.0);
592 if (model != Teuchos::null) {
593 stepper->setModel(model);
594 stepper->initialize();