37 const Teuchos::RCP<SolutionHistory<Scalar> >& solutionHistory)
41 int numStates = solutionHistory->getNumStates();
43 TEUCHOS_TEST_FOR_EXCEPTION(numStates < 1, std::logic_error,
44 "Error - setInitialConditions() needs at least one SolutionState\n"
45 " to set the initial condition. Number of States = " << numStates);
47 RCP<SolutionState<Scalar> > initialState = solutionHistory->getCurrentState();
48 RCP<Thyra::VectorBase<Scalar> > x = initialState->getX();
49 RCP<Thyra::VectorBase<Scalar> > xDot = initialState->getXDot();
50 if (xDot() == Teuchos::null) xDot = this->getStepperXDot();
53 auto inArgs = this->wrapperModel_->getNominalValues();
55 RCP<Thyra::VectorBase<Scalar> > initialXDot = initialState->getXDot();
58 TEUCHOS_TEST_FOR_EXCEPTION(
59 !((x != Teuchos::null && initialXDot != Teuchos::null) ||
60 (inArgs.get_x() != Teuchos::null &&
61 inArgs.get_x_dot() != Teuchos::null)), std::logic_error,
62 "Error - We need to set the initial conditions for x and xDot from\n"
63 " either initialState or appModel_->getNominalValues::InArgs\n"
64 " (but not from a mixture of the two).\n");
69 if (x == Teuchos::null) {
70 TEUCHOS_TEST_FOR_EXCEPTION( (x == Teuchos::null) &&
71 (inArgs.get_x() == Teuchos::null), std::logic_error,
72 "Error - setInitialConditions needs the ICs from the SolutionHistory\n"
73 " or getNominalValues()!\n");
75 x = Teuchos::rcp_const_cast<Thyra::VectorBase<Scalar> >(inArgs.get_x());
76 initialState->setX(x);
81 using Teuchos::rcp_const_cast;
82 if ( x == Teuchos::null || xDot == Teuchos::null ) {
83 TEUCHOS_TEST_FOR_EXCEPTION( (inArgs.get_x() == Teuchos::null) ||
84 (inArgs.get_x_dot() == Teuchos::null), std::logic_error,
85 "Error - setInitialConditions() needs the ICs from the initialState\n"
86 " or getNominalValues()!\n");
87 x = rcp_const_cast<Thyra::VectorBase<Scalar> >(inArgs.get_x());
88 initialState->setX(x);
89 RCP<Thyra::VectorBase<Scalar> > x_dot =
90 rcp_const_cast<Thyra::VectorBase<Scalar> >(inArgs.get_x_dot());
91 initialState->setXDot(x_dot);
97 std::string icConsistency = this->getICConsistency();
98 if (icConsistency ==
"None") {
100 if (initialState->getXDot() == Teuchos::null)
101 Thyra::assign(xDot.ptr(), Scalar(0.0));
104 if (initialState->getXDotDot() == Teuchos::null)
105 Thyra::assign(this->getStepperXDotDot(initialState).ptr(), Scalar(0.0));
108 else if (icConsistency ==
"Zero") {
110 Thyra::assign(xDot.ptr(), Scalar(0.0));
112 Thyra::assign(this->getStepperXDotDot(initialState).ptr(), Scalar(0.0));
114 else if (icConsistency ==
"App") {
116 auto x_dot = Teuchos::rcp_const_cast<Thyra::VectorBase<Scalar> >(
118 TEUCHOS_TEST_FOR_EXCEPTION(x_dot == Teuchos::null, std::logic_error,
119 "Error - setInitialConditions() requested 'App' for IC consistency,\n"
120 " but 'App' returned a null pointer for xDot!\n");
121 Thyra::assign(xDot.ptr(), *x_dot);
124 auto xDotDot = Teuchos::rcp_const_cast<Thyra::VectorBase<Scalar> >(
125 inArgs.get_x_dot_dot());
126 TEUCHOS_TEST_FOR_EXCEPTION(xDotDot == Teuchos::null, std::logic_error,
127 "Error - setInitialConditions() requested 'App' for IC consistency,\n"
128 " but 'App' returned a null pointer for xDotDot!\n");
129 Thyra::assign(this->getStepperXDotDot(initialState).ptr(), *xDotDot);
132 else if (icConsistency ==
"Consistent") {
135 const Scalar time = initialState->getTime();
136 const Scalar dt = initialState->getTimeStep();
137 RCP<TimeDerivative<Scalar> > timeDer = Teuchos::null;
138 const Scalar alpha = Scalar(1.0);
139 const Scalar beta = Scalar(0.0);
140 auto p = Teuchos::rcp(
new ImplicitODEParameters<Scalar>(
143 const Thyra::SolveStatus<Scalar> sStatus =
144 this->solveImplicitODE(x, xDot, time, p);
146 TEUCHOS_TEST_FOR_EXCEPTION(
147 sStatus.solveStatus != Thyra::SOLVE_STATUS_CONVERGED, std::logic_error,
148 "Error - Solver failed while determining the initial conditions.\n"
149 " Solver status is "<<Thyra::toString(sStatus.solveStatus)<<
".\n");
153 TEUCHOS_TEST_FOR_EXCEPTION(
true, std::logic_error,
154 "Error - setInitialConditions(): 'Consistent' for second-order ODE\n"
155 " has not been implemented.\n");
159 TEUCHOS_TEST_FOR_EXCEPTION(
true, std::logic_error,
160 "Error - setInitialConditions() invalid IC consistency, "
161 << icConsistency <<
".\n");
166 initialState->setIsSynced(
true);
169 if (this->getICConsistencyCheck()) {
170 auto f = initialState->getX()->clone_v();
172 const Scalar time = initialState->getTime();
173 const Scalar dt = initialState->getTimeStep();
174 RCP<TimeDerivative<Scalar> > timeDer = Teuchos::null;
175 const Scalar alpha = Scalar(0.0);
176 const Scalar beta = Scalar(0.0);
177 auto p = Teuchos::rcp(
new ImplicitODEParameters<Scalar>(
180 this->evaluateImplicitODE(f, x, xDot, time, p);
182 Scalar normX = Thyra::norm(*x);
183 Scalar reldiff = Scalar(0.0);
184 if (normX == Scalar(0.0)) reldiff = Thyra::norm(*f);
185 else reldiff = Thyra::norm(*f)/normX;
187 Scalar eps = Scalar(100.0)*std::abs(Teuchos::ScalarTraits<Scalar>::eps());
188 RCP<Teuchos::FancyOStream> out = this->getOStream();
189 Teuchos::OSTab ostab(out,1,
"StepperImplicit::setInitialConditions()");
191 *out <<
"\n---------------------------------------------------\n"
192 <<
"Info -- Stepper = " << this->getStepperType() <<
"\n"
193 <<
" Initial condition PASSED consistency check!\n"
194 <<
" (||f(x,xDot,t)||/||x|| = " << reldiff <<
") < "
195 <<
"(eps = " << eps <<
")" << std::endl
196 <<
"---------------------------------------------------\n"<<std::endl;
198 *out <<
"\n---------------------------------------------------\n"
199 <<
"Info -- Stepper = " << this->getStepperType() <<
"\n"
200 <<
" Initial condition FAILED consistency check but continuing!\n"
201 <<
" (||f(x,xDot,t)||/||x|| = " << reldiff <<
") > "
202 <<
"(eps = " << eps <<
")" << std::endl
203 <<
" ||f(x,xDot,t)|| = " << Thyra::norm(*f) << std::endl
204 <<
" ||x|| = " << Thyra::norm(*x) << std::endl
205 <<
"---------------------------------------------------\n"<<std::endl;
246 const Teuchos::RCP<ImplicitODEParameters<Scalar> > & p,
250 typedef Thyra::ModelEvaluatorBase MEB;
251 MEB::InArgs<Scalar> inArgs = wrapperModel_->getInArgs();
252 MEB::OutArgs<Scalar> outArgs = wrapperModel_->getOutArgs();
254 if ( y != Teuchos::null ) inArgs.set_p(index, y);
255 if (inArgs.supports(MEB::IN_ARG_x_dot )) inArgs.set_x_dot (xDot);
256 if (inArgs.supports(MEB::IN_ARG_t )) inArgs.set_t (time);
257 if (inArgs.supports(MEB::IN_ARG_step_size))
258 inArgs.set_step_size(p->timeStepSize_);
259 if (inArgs.supports(MEB::IN_ARG_alpha )) inArgs.set_alpha (p->alpha_);
260 if (inArgs.supports(MEB::IN_ARG_beta )) inArgs.set_beta (p->beta_);
261 if (inArgs.supports(MEB::IN_ARG_stage_number))
262 inArgs.set_stage_number(p->stageNumber_);
264 wrapperModel_->setForSolve(p->timeDer_, inArgs, outArgs, p->evaluationType_);
266 Thyra::SolveStatus<Scalar> sStatus;
267 typedef Teuchos::ScalarTraits<Scalar> ST;
268 switch (p->evaluationType_)
271 if (getZeroInitialGuess()) Thyra::assign(x.ptr(), ST::zero());
272 sStatus = (*solver_).solve(&*x);
277 sStatus = (*solver_).solve(&*xDot);
281 TEUCHOS_TEST_FOR_EXCEPT(
"Invalid EVALUATION_TYPE!");