97 Thyra::ModelEvaluator<Scalar>& ,
98 const Ptr<IntegratorBase<Scalar> >& forwardIntegrator
102 using Teuchos::describe;
103 using Teuchos::outArg;
104 using Teuchos::rcp_dynamic_cast;
105 using Teuchos::dyn_cast;
106 using Teuchos::dyn_cast;
107 using Teuchos::OSTab;
108 using Teuchos::sublist;
109 typedef Thyra::ModelEvaluatorBase MEB;
110 using Thyra::createMember;
111 namespace BDASTU = BasicDiscreteAdjointStepperTesterUtils;
113 const RCP<Teuchos::FancyOStream> out = this->getOStream();
114 const Teuchos::EVerbosityLevel verbLevel = this->getVerbLevel();
116 const RCP<ParameterList> paramList = this->getMyNonconstParamList();
123 *out <<
"\n*** Entering BasicDiscreteAdjointStepperTester<Scalar>::testAdjointStepper(...) ...\n";
126 const TimeRange<Scalar> fwdTimeRange = forwardIntegrator->getFwdTimeRange();
127 const RCP<Rythmos::StepperBase<Scalar> > fwdStepper =
128 forwardIntegrator->getNonconstStepper();
129 const RCP<const Thyra::ModelEvaluator<Scalar> > fwdModel =
130 fwdStepper->getModel();
131 const RCP<const Thyra::VectorSpaceBase<Scalar> > f_space = fwdModel->get_f_space();
134 *out <<
"\nA) Construct the IC basis B ...\n";
137 const RCP<Thyra::MultiVectorBase<Scalar> > B =
138 createMembers(fwdModel->get_x_space(), 1,
"B");
139 Thyra::seed_randomize<Scalar>(0);
140 Thyra::randomize<Scalar>( -1.0, +1.0, B.ptr() );
142 *out <<
"\nB: " << describe(*B, verbLevel);
145 *out <<
"\nB) Construct the forward sensitivity integrator ...\n";
148 const RCP<ForwardSensitivityStepper<Scalar> > fwdSensStepper =
149 forwardSensitivityStepper<Scalar>();
150 fwdSensStepper->initializeSyncedSteppersInitCondOnly(
153 fwdStepper->getInitialCondition(),
155 dyn_cast<SolverAcceptingStepperBase<Scalar> >(*fwdStepper).getNonconstSolver()
157 *out <<
"\nfwdSensStepper: " << describe(*fwdSensStepper, verbLevel);
159 const MEB::InArgs<Scalar> fwdIC =
160 forwardIntegrator->getStepper()->getInitialCondition();
162 MEB::InArgs<Scalar> fwdSensIC =
163 createStateAndSensInitialCondition<Scalar>(*fwdSensStepper, fwdIC, B);
164 fwdSensStepper->setInitialCondition(fwdSensIC);
166 const RCP<IntegratorBase<Scalar> > fwdSensIntegrator =
167 forwardIntegrator->cloneIntegrator();
168 fwdSensIntegrator->setStepper(fwdSensStepper,
169 forwardIntegrator->getFwdTimeRange().upper());
170 *out <<
"\nfwdSensIntegrator: " << describe(*fwdSensIntegrator, verbLevel);
173 *out <<
"\nC) Construct the adjoint sensitivity integrator ...\n";
176 RCP<AdjointModelEvaluator<Scalar> > adjModel =
177 adjointModelEvaluator<Scalar>(fwdModel, fwdTimeRange);
178 adjModel->setFwdStateSolutionBuffer(
179 dyn_cast<DefaultIntegrator<Scalar> >(*forwardIntegrator).getTrailingInterpolationBuffer()
186 *out <<
"\nadjModel: " << describe(*adjModel, verbLevel);
189 const RCP<Thyra::VectorBase<Scalar> > lambda_ic = createMember(f_space);
190 V_S( lambda_ic.ptr(), 0.0 );
193 const RCP<Thyra::VectorBase<Scalar> > lambda_dot_ic = createMember(f_space);
194 Thyra::V_S( lambda_dot_ic.ptr(), 0.0 );
196 MEB::InArgs<Scalar> adj_ic = adjModel->getNominalValues();
197 adj_ic.set_x(lambda_ic);
198 adj_ic.set_x_dot(lambda_dot_ic);
199 *out <<
"\nadj_ic: " << describe(adj_ic, Teuchos::VERB_EXTREME);
201 const RCP<Thyra::LinearNonlinearSolver<Scalar> > adjTimeStepSolver =
202 Thyra::linearNonlinearSolver<Scalar>();
203 const RCP<Rythmos::StepperBase<Scalar> > adjStepper =
204 forwardIntegrator->getStepper()->cloneStepperAlgorithm();
205 *out <<
"\nadjStepper: " << describe(*adjStepper, verbLevel);
207 adjStepper->setInitialCondition(adj_ic);
209 const RCP<IntegratorBase<Scalar> > adjIntegrator = forwardIntegrator->cloneIntegrator();
210 adjIntegrator->setStepper(adjStepper, forwardIntegrator->getFwdTimeRange().length());
213 *out <<
"\nD) Solve the forawrd problem storing the state along the way ...\n";
216 const double t_final = fwdTimeRange.upper();
217 RCP<const Thyra::VectorBase<Scalar> > x_final, x_dot_final;
218 get_fwd_x_and_x_dot( *forwardIntegrator, t_final, outArg(x_final), outArg(x_dot_final) );
220 *out <<
"\nt_final = " << t_final <<
"\n";
221 *out <<
"\nx_final: " << *x_final;
222 *out <<
"\nx_dot_final: " << *x_dot_final;
225 *out <<
"\nE) Solve the forawrd sensitivity equations and compute fwd reduced sens ...\n";
229 TEUCHOS_TEST_FOR_EXCEPT(
true);
232 TEUCHOS_TEST_FOR_EXCEPT(
true);
235 TEUCHOS_TEST_FOR_EXCEPT(
true);
238 *out <<
"\nF) Solve the adjoint equations and compute adj reduced sens ...\n";
242 TEUCHOS_TEST_FOR_EXCEPT(
true);
245 TEUCHOS_TEST_FOR_EXCEPT(
true);
248 TEUCHOS_TEST_FOR_EXCEPT(
true);
251 *out <<
"\nG) Compare forward and adjoint reduced sens ...\n";
254 TEUCHOS_TEST_FOR_EXCEPT(
true);
257 *out <<
"\n*** Leaving BasicDiscreteAdjointStepperTester<Scalar>::testAdjointStepper(...) ...\n";