97 const Ptr<IntegratorBase<Scalar> > &fwdSensIntegrator
101 using Teuchos::outArg;
102 using Teuchos::rcp_dynamic_cast;
103 using Teuchos::OSTab;
104 using Teuchos::sublist;
105 typedef Thyra::ModelEvaluatorBase MEB;
106 namespace FSSTU = ForwardSensitivityStepperTesterUtils;
108 const Teuchos::EVerbosityLevel verbLevel = this->getVerbLevel();
109 const RCP<Teuchos::FancyOStream> out = this->getOStream();
111 const RCP<ParameterList> paramList = this->getMyNonconstParamList();
116 RCP<ForwardSensitivityStepper<Scalar> > fwdSensStepper =
117 rcp_dynamic_cast<ForwardSensitivityStepper<Scalar> >(
118 fwdSensIntegrator->getNonconstStepper()
122 RCP<StepperBase<Scalar> > stateStepper = fwdSensStepper->getNonconstStateStepper();
126 MEB::InArgs<Scalar> state_and_sens_ic = fwdSensStepper->getInitialCondition();
127 MEB::InArgs<Scalar> state_ic =
128 extractStateInitialCondition(*fwdSensStepper, state_and_sens_ic);
132 RCP<StepperAsModelEvaluator<Scalar> >
133 stateIntegratorAsModel = stepperAsModelEvaluator(
134 stateStepper, fwdSensIntegrator->cloneIntegrator(), state_ic
140 const Scalar finalTime = fwdSensIntegrator->getFwdTimeRange().upper();
142 *out <<
"\nCompute discrete forward sensitivities using fwdSensIntegrator ...\n";
144 RCP<const VectorBase<Scalar> > x_bar_final, x_bar_dot_final;
147 get_fwd_x_and_x_dot( *fwdSensIntegrator, finalTime,
148 outArg(x_bar_final), outArg(x_bar_dot_final) );
150 <<
"\nx_bar_final = x_bar(p,finalTime) evaluated using stateAndSensIntegratorAsModel:\n"
151 << describe(*x_bar_final, verbLevel);
156 *out <<
"\nCompute discrete forward sensitivities by finite differences ...\n";
158 RCP<Thyra::MultiVectorBase<Scalar> > DxDp_fd_final;
162 Thyra::DirectionalFiniteDiffCalculator<Scalar> fdCalc;
163 if (nonnull(paramList))
164 fdCalc.setParameterList(sublist(paramList, FSSTU::FdCalc_name));
165 fdCalc.setOStream(out);
166 fdCalc.setVerbLevel(Teuchos::VERB_NONE);
169 fdBasePoint = stateIntegratorAsModel->createInArgs();
171 fdBasePoint.setArgs(state_ic,
true);
172 fdBasePoint.set_t(finalTime);
174 const int g_index = 0;
175 const int p_index = fwdSensStepper->getFwdSensModel()->get_p_index();
179 stateIntegratorAsModel->get_g_space(g_index),
180 stateIntegratorAsModel->get_p_space(p_index)->dim()
183 typedef Thyra::DirectionalFiniteDiffCalculatorTypes::SelectedDerivatives
186 MEB::OutArgs<Scalar> fdOutArgs =
187 fdCalc.createOutArgs(
188 *stateIntegratorAsModel,
189 SelectedDerivatives().supports(MEB::OUT_ARG_DgDp, g_index, p_index)
191 fdOutArgs.set_DgDp(g_index, p_index, DxDp_fd_final);
195 stateStepper->setVerbLevel(Teuchos::VERB_NONE);
196 stateIntegratorAsModel->setVerbLevel(Teuchos::VERB_NONE);
198 fdCalc.calcDerivatives(
199 *stateIntegratorAsModel, fdBasePoint,
200 stateIntegratorAsModel->createOutArgs(),
205 <<
"\nFinite difference DxDp_fd_final = DxDp(p,finalTime): "
206 << describe(*DxDp_fd_final, verbLevel);
212 *out <<
"\nChecking that integrated DxDp(p,finalTime) and finite-diff DxDp(p,finalTime) are similar ...\n";
216 Teuchos::OSTab tab(out);
218 RCP<const Thyra::VectorBase<Scalar> >
219 DxDp_vec_final = Thyra::productVectorBase<Scalar>(x_bar_final)->getVectorBlock(1);
221 RCP<const Thyra::VectorBase<Scalar> >
222 DxDp_fd_vec_final = Thyra::multiVectorProductVector(
223 rcp_dynamic_cast<
const Thyra::DefaultMultiVectorProductVectorSpace<Scalar> >(
224 DxDp_vec_final->range()
229 const bool result = Thyra::testRelNormDiffErr(
230 "DxDp_vec_final", *DxDp_vec_final,
231 "DxDp_fd_vec_final", *DxDp_fd_vec_final,
232 "maxSensError", errorTol_,