32 const Teuchos::RCP<SolutionHistory<Scalar> >& solutionHistory)
36 int numStates = solutionHistory->getNumStates();
38 TEUCHOS_TEST_FOR_EXCEPTION(numStates < 1, std::logic_error,
39 "Error - setInitialConditions() needs at least one SolutionState\n"
40 " to set the initial condition. Number of States = " << numStates);
43 RCP<Teuchos::FancyOStream> out = this->getOStream();
44 Teuchos::OSTab ostab(out,1,
"StepperExplicit::setInitialConditions()");
45 *out <<
"Warning -- SolutionHistory has more than one state!\n"
46 <<
"Setting the initial conditions on the currentState.\n"<<std::endl;
49 RCP<SolutionState<Scalar> > initialState = solutionHistory->getCurrentState();
50 RCP<Thyra::VectorBase<Scalar> > x = initialState->getX();
51 RCP<Thyra::VectorBase<Scalar> > xDot = initialState->getXDot();
52 if (xDot == Teuchos::null) xDot = this->getStepperXDot();
55 inArgs_ = appModel_->getNominalValues();
57 RCP<Thyra::VectorBase<Scalar> > initialXDot = initialState->getXDot();
60 TEUCHOS_TEST_FOR_EXCEPTION(
61 !((x != Teuchos::null && initialXDot != Teuchos::null) ||
62 (inArgs_.get_x() != Teuchos::null &&
63 inArgs_.get_x_dot() != Teuchos::null)), std::logic_error,
64 "Error - We need to set the initial conditions for x and xDot from\n"
65 " either initialState or appModel_->getNominalValues::InArgs\n"
66 " (but not from a mixture of the two).\n");
70 if (x == Teuchos::null) {
72 TEUCHOS_TEST_FOR_EXCEPTION( (x == Teuchos::null) &&
73 (inArgs_.get_x() == Teuchos::null), std::logic_error,
74 "Error - setInitialConditions needs the ICs from the SolutionHistory\n"
75 " or getNominalValues()!\n");
77 x = Teuchos::rcp_const_cast<Thyra::VectorBase<Scalar> >(inArgs_.get_x());
78 initialState->setX(x);
82 using Teuchos::rcp_const_cast;
84 if ( initialState->getX() == Teuchos::null ||
85 initialState->getXDot() == Teuchos::null ) {
86 TEUCHOS_TEST_FOR_EXCEPTION( (inArgs_.get_x() == Teuchos::null) ||
87 (inArgs_.get_x_dot() == Teuchos::null), std::logic_error,
88 "Error - setInitialConditions() needs the ICs from the initialState\n"
89 " or getNominalValues()!\n");
90 x = rcp_const_cast<Thyra::VectorBase<Scalar> >(inArgs_.get_x());
91 initialState->setX(x);
92 RCP<Thyra::VectorBase<Scalar> > x_dot
93 = rcp_const_cast<Thyra::VectorBase<Scalar> >(inArgs_.get_x_dot());
94 initialState->setXDot(x_dot);
99 std::string icConsistency = this->getICConsistency();
100 if (icConsistency ==
"None") {
102 if (initialState->getXDot() == Teuchos::null) {
103 RCP<Teuchos::FancyOStream> out = this->getOStream();
104 Teuchos::OSTab ostab(out,1,
"StepperExplicit::setInitialConditions()");
105 *out <<
"Warning -- Requested IC consistency of 'None' but\n"
106 <<
" initialState does not have an xDot.\n"
107 <<
" Setting a 'Consistent' xDot!\n" << std::endl;
108 auto p = Teuchos::rcp(
new ExplicitODEParameters<Scalar>(0.0));
109 evaluateExplicitODE(xDot, x, initialState->getTime(), p);
110 initialState->setIsSynced(
true);
114 if (initialState->getXDotDot() == Teuchos::null) {
115 RCP<Teuchos::FancyOStream> out = this->getOStream();
116 Teuchos::OSTab ostab(out,1,
"StepperExplicit::setInitialConditions()");
117 *out <<
"Warning -- Requested IC consistency of 'None' but\n"
118 <<
" initialState does not have an xDotDot.\n"
119 <<
" Setting a 'Consistent' xDotDot!\n" << std::endl;
120 auto p = Teuchos::rcp(
new ExplicitODEParameters<Scalar>(0.0));
121 this->evaluateExplicitODE(this->getStepperXDotDot(initialState), x,
122 initialState->getXDot(),
123 initialState->getTime(), p);
124 initialState->setIsSynced(
true);
128 else if (icConsistency ==
"Zero") {
130 Thyra::assign(xDot.ptr(), Scalar(0.0));
132 Thyra::assign(this->getStepperXDotDot(initialState).ptr(), Scalar(0.0));
134 else if (icConsistency ==
"App") {
136 auto x_dot = Teuchos::rcp_const_cast<Thyra::VectorBase<Scalar> >(
137 inArgs_.get_x_dot());
138 TEUCHOS_TEST_FOR_EXCEPTION(xDot == Teuchos::null, std::logic_error,
139 "Error - setInitialConditions() requested 'App' for IC consistency,\n"
140 " but 'App' returned a null pointer for xDot!\n");
141 Thyra::assign(xDot.ptr(), *x_dot);
144 auto xDotDot = Teuchos::rcp_const_cast<Thyra::VectorBase<Scalar> >(
145 inArgs_.get_x_dot_dot());
146 TEUCHOS_TEST_FOR_EXCEPTION(xDotDot == Teuchos::null, std::logic_error,
147 "Error - setInitialConditions() requested 'App' for IC consistency,\n"
148 " but 'App' returned a null pointer for xDotDot!\n");
149 Thyra::assign(this->getStepperXDotDot(initialState).ptr(), *xDotDot);
152 else if (icConsistency ==
"Consistent") {
155 auto p = Teuchos::rcp(
new ExplicitODEParameters<Scalar>(0.0));
156 evaluateExplicitODE(xDot, x, initialState->getTime(), p);
160 auto p = Teuchos::rcp(
new ExplicitODEParameters<Scalar>(0.0));
161 this->evaluateExplicitODE(initialState->getXDotDot(), x,
162 initialState->getXDot(),
163 initialState->getTime(), p);
167 TEUCHOS_TEST_FOR_EXCEPTION(
true, std::logic_error,
168 "Error - setInitialConditions() invalid IC consistency, '"
169 << icConsistency <<
"'.\n");
174 initialState->setIsSynced(
true);
177 if (this->getICConsistencyCheck()) {
179 auto f = initialState->getX()->clone_v();
180 auto p = Teuchos::rcp(
new ExplicitODEParameters<Scalar>(0.0));
181 evaluateExplicitODE(f, x, initialState->getTime(), p);
182 Thyra::Vp_StV(f.ptr(), Scalar(-1.0), *(xDot));
183 Scalar normX = Thyra::norm(*x);
184 Scalar reldiff = Scalar(0.0);
185 if (normX == Scalar(0.0)) reldiff = Thyra::norm(*f);
186 else reldiff = Thyra::norm(*f)/normX;
188 Scalar eps = Scalar(100.0)*std::abs(Teuchos::ScalarTraits<Scalar>::eps());
189 RCP<Teuchos::FancyOStream> out = this->getOStream();
190 Teuchos::OSTab ostab(out,1,
"StepperExplicit::setInitialConditions()");
192 *out <<
"\n---------------------------------------------------\n"
193 <<
"Info -- Stepper = " << this->getStepperType() <<
"\n"
194 <<
" Initial condition PASSED consistency check!\n"
195 <<
" (||xDot-f(x,t)||/||x|| = " << reldiff <<
") < "
196 <<
"(eps = " << eps <<
")" << std::endl
197 <<
"---------------------------------------------------\n"<<std::endl;
199 *out <<
"\n---------------------------------------------------\n"
200 <<
"Info -- Stepper = " << this->getStepperType() <<
"\n"
201 <<
" Initial condition FAILED consistency check but continuing!\n"
202 <<
" (||xDot-f(x,t)||/||x|| = " << reldiff <<
") > "
203 <<
"(eps = " << eps <<
")" << std::endl
204 <<
" ||xDot-f(x,t)|| = " << Thyra::norm(*f) << std::endl
205 <<
" ||x|| = " << Thyra::norm(*x) << std::endl
206 <<
"---------------------------------------------------\n"<<std::endl;
210 auto xDotDot = initialState->getXDotDot();
211 auto f = initialState->getX()->clone_v();
212 auto p = Teuchos::rcp(
new ExplicitODEParameters<Scalar>(0.0));
213 this->evaluateExplicitODE(f, x, initialState->getXDot(),
214 initialState->getTime(), p);
215 Thyra::Vp_StV(f.ptr(), Scalar(-1.0), *(xDotDot));
216 Scalar normX = Thyra::norm(*x);
217 Scalar reldiff = Scalar(0.0);
218 if (normX == Scalar(0.0)) reldiff = Thyra::norm(*f);
219 else reldiff = Thyra::norm(*f)/normX;
221 Scalar eps = Scalar(100.0)*std::abs(Teuchos::ScalarTraits<Scalar>::eps());
222 RCP<Teuchos::FancyOStream> out = this->getOStream();
223 Teuchos::OSTab ostab(out,1,
"StepperExplicit::setInitialConditions()");
225 *out <<
"\n---------------------------------------------------\n"
226 <<
"Info -- Stepper = " << this->getStepperType() <<
"\n"
227 <<
" Initial condition PASSED consistency check!\n"
228 <<
" (||xDotDot-f(x,xDot,t)||/||x|| = " << reldiff <<
") > "
229 <<
"(eps = " << eps <<
")" << std::endl
230 <<
"---------------------------------------------------\n"<<std::endl;
232 *out <<
"\n---------------------------------------------------\n"
233 <<
"Info -- Stepper = " << this->getStepperType() <<
"\n"
234 <<
"Initial condition FAILED consistency check but continuing!\n"
235 <<
" (||xDotDot-f(x,xDot,t)||/||x|| = " << reldiff <<
") > "
236 <<
"(eps = " << eps <<
")" << std::endl
237 <<
" ||xDotDot-f(x,xDot,t)|| = " << Thyra::norm(*f) << std::endl
238 <<
" ||x|| = " << Thyra::norm(*x) << std::endl
239 <<
"---------------------------------------------------\n"<<std::endl;