101 const Teuchos::RCP<SolutionHistory<Scalar> >& solutionHistory)
105 int numStates = solutionHistory->getNumStates();
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);
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;
119 RCP<SolutionState<Scalar> > initialState = solutionHistory->getCurrentState();
120 RCP<Thyra::VectorBase<Scalar> > x = initialState->getX();
121 RCP<Thyra::VectorBase<Scalar> > xDot = initialState->getXDot();
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");
134 this->inArgs_ = this->appModel_->getNominalValues();
135 using Teuchos::rcp_const_cast;
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);
148 this->inArgs_.set_x(x);
149 this->inArgs_.set_x_dot(xDot);
153 if (initialState->getXDotDot() == Teuchos::null)
154 initialState->setXDotDot(initialState->getX()->clone_v());
156 this->setStepperXDotDot(initialState->getXDotDot());
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);
171 initialState->setIsSynced(
true);
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);
184 else if (icConsistency ==
"Consistent") {
186 auto p = Teuchos::rcp(
new ExplicitODEParameters<Scalar>(0.0));
187 this->evaluateExplicitODE(initialState->getXDotDot(), x, xDot,
188 initialState->getTime(), p);
192 initialState->setIsSynced(
true);
195 TEUCHOS_TEST_FOR_EXCEPTION(
true, std::logic_error,
196 "Error - setInitialConditions() invalid IC consistency, "
197 << icConsistency <<
".\n");
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);
210 Scalar eps = Scalar(100.0)*std::abs(Teuchos::ScalarTraits<Scalar>::eps());
211 if (normxDotDot > eps*reldiff) reldiff /= normxDotDot;
213 RCP<Teuchos::FancyOStream> out = this->getOStream();
214 Teuchos::OSTab ostab(out,1,
"StepperNewmarkExplicitAForm::setInitialConditions()");
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;
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;
238 const Teuchos::RCP<SolutionHistory<Scalar> >& solutionHistory)
240 this->checkInitialized();
244 TEMPUS_FUNC_TIME_MONITOR(
"Tempus::StepperNewmarkExplicitAForm::takeStep()");
246 TEUCHOS_TEST_FOR_EXCEPTION(solutionHistory->getNumStates() < 2,
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");
254 auto thisStepper = Teuchos::rcpFromRef(*
this);
255 stepperNewmarkExpAppAction_->execute(solutionHistory, thisStepper,
258 RCP<SolutionState<Scalar> > currentState=solutionHistory->getCurrentState();
259 RCP<SolutionState<Scalar> > workingState=solutionHistory->getWorkingState();
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();
267 const Scalar dt = workingState->getTimeStep();
268 const Scalar time_old = currentState->getTime();
270 auto p = Teuchos::rcp(
new ExplicitODEParameters<Scalar>(dt));
271 if (!(this->getUseFSAL()) || workingState->getNConsecutiveFailures() != 0) {
273 this->evaluateExplicitODE(a_old, d_old, v_old, time_old, p);
277 currentState->setIsSynced(
true);
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();
286 predictDisplacement(*d_new, *d_old, *v_old, *a_old, dt);
287 predictVelocity(*v_new, *v_old, *a_old, dt);
289 stepperNewmarkExpAppAction_->execute(solutionHistory, thisStepper,
293 this->evaluateExplicitODE(a_new, d_new, v_new, time_old, p);
295 stepperNewmarkExpAppAction_->execute(solutionHistory, thisStepper,
299 correctVelocity(*v_new, *v_new, *a_new, dt);
301 if ( this->getUseFSAL() ) {
303 const Scalar time_new = workingState->getTime();
304 this->evaluateExplicitODE(a_new, d_new, v_new, time_new, p);
308 workingState->setIsSynced(
true);
310 assign(a_new.ptr(), Teuchos::ScalarTraits<Scalar>::zero());
311 workingState->setIsSynced(
false);
315 workingState->setOrder(this->getOrder());
316 workingState->computeNorms(currentState);
318 stepperNewmarkExpAppAction_->execute(solutionHistory, thisStepper,