37 const Teuchos::RCP<Thyra::NonlinearSolverBase<Scalar> >& solver,
39 std::string ICConsistency,
40 bool ICConsistencyCheck,
41 bool zeroInitialGuess,
42 const Teuchos::RCP<StepperTrapezoidalAppAction<Scalar> >& stepperTrapAppAction)
44 this->setStepperName(
"Trapezoidal Method");
45 this->setStepperType(
"Trapezoidal Method");
46 this->setUseFSAL( useFSAL);
47 this->setICConsistency( ICConsistency);
48 this->setICConsistencyCheck( ICConsistencyCheck);
49 this->setZeroInitialGuess( zeroInitialGuess);
51 this->setAppAction(stepperTrapAppAction);
52 this->setSolver(solver);
54 if (appModel != Teuchos::null) {
55 this->setModel(appModel);
78 const Teuchos::RCP<SolutionHistory<Scalar> >& solutionHistory)
82 RCP<SolutionState<Scalar> > initialState = solutionHistory->getCurrentState();
85 if (initialState->getXDot() == Teuchos::null)
86 this->setStepperXDot(initialState->getX()->clone_v());
88 this->setStepperXDot(initialState->getXDot());
92 TEUCHOS_TEST_FOR_EXCEPTION( !(this->getUseFSAL()), std::logic_error,
93 "Error - The First-Same-As-Last (FSAL) principle is required\n"
94 " for the Trapezoidal Stepper (i.e., useFSAL=true)!\n");
107 const Teuchos::RCP<SolutionHistory<Scalar> >& solutionHistory)
109 this->checkInitialized();
113 TEMPUS_FUNC_TIME_MONITOR(
"Tempus::StepperTrapezoidal::takeStep()");
115 TEUCHOS_TEST_FOR_EXCEPTION(solutionHistory->getNumStates() < 2,
117 "Error - StepperTrapezoidal<Scalar>::takeStep(...)\n"
118 "Need at least two SolutionStates for Trapezoidal.\n"
119 " Number of States = " << solutionHistory->getNumStates() <<
"\n"
120 "Try setting in \"Solution History\" \"Storage Type\" = \"Undo\"\n"
121 " or \"Storage Type\" = \"Static\" and \"Storage Limit\" = \"2\"\n");
122 RCP<StepperTrapezoidal<Scalar> > thisStepper = Teuchos::rcpFromRef(*
this);
123 stepperTrapAppAction_->execute(solutionHistory, thisStepper,
126 RCP<SolutionState<Scalar> > workingState=solutionHistory->getWorkingState();
127 RCP<SolutionState<Scalar> > currentState=solutionHistory->getCurrentState();
129 RCP<const Thyra::VectorBase<Scalar> > xOld = currentState->getX();
130 RCP<const Thyra::VectorBase<Scalar> > xDotOld = currentState->getXDot();
131 RCP<Thyra::VectorBase<Scalar> > x = workingState->getX();
132 if (workingState->getXDot() != Teuchos::null)
133 this->setStepperXDot(workingState->getXDot());
134 RCP<Thyra::VectorBase<Scalar> > xDot = this->getStepperXDot();
136 const Scalar time = workingState->getTime();
137 const Scalar dt = workingState->getTimeStep();
138 const Scalar alpha = getAlpha(dt);
139 const Scalar beta = getBeta (dt);
142 Teuchos::RCP<TimeDerivative<Scalar> > timeDer =
143 Teuchos::rcp(
new StepperTrapezoidalTimeDerivative<Scalar>(
144 alpha, xOld, xDotOld));
146 auto p = Teuchos::rcp(
new ImplicitODEParameters<Scalar>(
147 timeDer, dt, alpha, beta));
148 stepperTrapAppAction_->execute(solutionHistory, thisStepper,
151 const Thyra::SolveStatus<Scalar> sStatus =
152 this->solveImplicitODE(x, xDot, time, p);
153 stepperTrapAppAction_->execute(solutionHistory, thisStepper,
156 if (workingState->getXDot() != Teuchos::null)
157 timeDer->compute(x, xDot);
159 workingState->setSolutionStatus(sStatus);
160 workingState->setOrder(this->getOrder());
161 workingState->computeNorms(currentState);
162 stepperTrapAppAction_->execute(solutionHistory, thisStepper,