Tempus Version of the Day
Time Integration
Loading...
Searching...
No Matches
Tempus_IntegratorAdjointSensitivity_impl.hpp
Go to the documentation of this file.
1// @HEADER
2// ****************************************************************************
3// Tempus: Copyright (2017) Sandia Corporation
4//
5// Distributed under BSD 3-clause license (See accompanying file Copyright.txt)
6// ****************************************************************************
7// @HEADER
8
9#ifndef Tempus_IntegratorAdjointSensitivity_impl_hpp
10#define Tempus_IntegratorAdjointSensitivity_impl_hpp
11
12#include "Teuchos_ParameterList.hpp"
13#include "Thyra_DefaultMultiVectorProductVectorSpace.hpp"
14#include "Thyra_DefaultMultiVectorProductVector.hpp"
15#include "Thyra_DefaultProductVectorSpace.hpp"
16#include "Thyra_DefaultProductVector.hpp"
17#include "Thyra_VectorStdOps.hpp"
18#include "Thyra_MultiVectorStdOps.hpp"
20#include "Thyra_LinearOpWithSolveFactoryHelpers.hpp"
21
22namespace Tempus {
23
24template <class Scalar>
26 const Teuchos::RCP<Thyra::ModelEvaluator<Scalar> > &model,
27 const Teuchos::RCP<IntegratorBasic<Scalar> > &state_integrator,
28 const Teuchos::RCP<Thyra::ModelEvaluator<Scalar> > &adjoint_model,
29 const Teuchos::RCP<AdjointAuxSensitivityModelEvaluator<Scalar> > &adjoint_aux_model,
30 const Teuchos::RCP<IntegratorBasic<Scalar> > &adjoint_integrator,
31 const Teuchos::RCP<SolutionHistory<Scalar> > &solutionHistory,
32 const int p_index,
33 const int g_index,
34 const bool g_depends_on_p,
35 const bool f_depends_on_p,
36 const bool ic_depends_on_p,
37 const bool mass_matrix_is_identity)
38 : model_(model)
39 , state_integrator_(state_integrator)
40 , adjoint_model_(adjoint_model)
41 , adjoint_aux_model_(adjoint_aux_model)
42 , adjoint_integrator_(adjoint_integrator)
43 , solutionHistory_(solutionHistory)
44 , p_index_(p_index)
45 , g_index_(g_index)
46 , g_depends_on_p_(g_depends_on_p)
47 , f_depends_on_p_(f_depends_on_p)
48 , ic_depends_on_p_(ic_depends_on_p)
49 , mass_matrix_is_identity_(mass_matrix_is_identity)
50 , stepMode_(SensitivityStepMode::Forward)
51{
52
53 TEUCHOS_TEST_FOR_EXCEPTION(
54 getStepper()->getUseFSAL(), std::logic_error,
55 "Error - IntegratorAdjointSensitivity(): Cannot use FSAL with\n"
56 " IntegratorAdjointSensitivity, because the state and adjoint\n"
57 " integrators require ModelEvaluator evaluation in the\n"
58 " constructor to make the initial conditions consistent.\n"
59 " For the adjoint integrator, this requires special construction\n"
60 " which has not been implemented yet.\n");
61}
62
63template<class Scalar>
66{
67 state_integrator_ = createIntegratorBasic<Scalar>();
68 adjoint_integrator_ = createIntegratorBasic<Scalar>();
70}
71
72template<class Scalar>
73bool
76{
77 const Scalar tfinal =
78 state_integrator_->getTimeStepControl()->getFinalTime();
79 return advanceTime(tfinal);
80}
81
82template<class Scalar>
83bool
85advanceTime(const Scalar timeFinal)
86{
87 TEMPUS_FUNC_TIME_MONITOR_DIFF("Tempus::IntegratorAdjointSensitivity::advanceTime()", TEMPUS_AS_AT);
88
89 using Teuchos::RCP;
90 using Teuchos::rcp_dynamic_cast;
92 using Thyra::VectorSpaceBase;
99 using Thyra::createMember;
100 using Thyra::createMembers;
101 using Thyra::assign;
102 typedef Thyra::ModelEvaluatorBase MEB;
103 typedef Thyra::DefaultMultiVectorProductVector<Scalar> DMVPV;
104 typedef Thyra::DefaultProductVector<Scalar> DPV;
105
106 // Get initial state for later
107 RCP<const SolutionHistory<Scalar> > state_solution_history =
108 state_integrator_->getSolutionHistory();
109 RCP<const SolutionState<Scalar> > initial_state =
110 (*state_solution_history)[0];
111
112 // Run state integrator and get solution
113 bool state_status = true;
114 {
115 TEMPUS_FUNC_TIME_MONITOR_DIFF("Tempus::IntegratorAdjointSensitivity::advanceTime::state", TEMPUS_AS_AT_FWD);
117 state_status = state_integrator_->advanceTime(timeFinal);
118 }
119
120 // For at least some time-stepping methods, the time of the last time step
121 // may not be timeFinal (e.g., it may be greater by at most delta_t).
122 // But since the adjoint model requires timeFinal in its formulation, reset
123 // it to the achieved final time.
124 adjoint_aux_model_->setFinalTime(state_integrator_->getTime());
125
126 // Set solution history in adjoint stepper
127 adjoint_aux_model_->setForwardSolutionHistory(state_solution_history);
128
129 // Compute dg/dx
130 RCP<const VectorSpaceBase<Scalar> > g_space = model_->get_g_space(g_index_);
131 RCP<const VectorSpaceBase<Scalar> > x_space = model_->get_x_space();
132 const int num_g = g_space->dim();
133 RCP<MultiVectorBase<Scalar> > dgdx = createMembers(x_space, num_g);
134 MEB::InArgs<Scalar> inargs = model_->getNominalValues();
135 RCP<const SolutionState<Scalar> > state =
136 state_solution_history->getCurrentState();
137 inargs.set_t(state->getTime());
138 inargs.set_x(state->getX());
139 inargs.set_x_dot(state->getXDot());
140 MEB::OutArgs<Scalar> outargs = model_->createOutArgs();
141 MEB::OutArgs<Scalar> adj_outargs = adjoint_model_->createOutArgs();
142 outargs.set_DgDx(g_index_,
143 MEB::Derivative<Scalar>(dgdx, MEB::DERIV_MV_GRADIENT_FORM));
144 model_->evalModel(inargs, outargs);
145 outargs.set_DgDx(g_index_, MEB::Derivative<Scalar>());
146
147 // Compute ICs == [ (df/dx_dot)^{-T} (dg/dx)^T; 0 ]
148 // For explicit form, we are relying on the user to inform us the
149 // the mass matrix is the identity. It would be nice to be able to determine
150 // somehow automatically that we are using an explicit stepper.
151 RCP<DPV> adjoint_init =
152 rcp_dynamic_cast<DPV>(Thyra::createMember(adjoint_aux_model_->get_x_space()));
153 RCP<MultiVectorBase<Scalar> > adjoint_init_mv =
154 rcp_dynamic_cast<DMVPV>(adjoint_init->getNonconstVectorBlock(0))->getNonconstMultiVector();
155 assign(adjoint_init->getNonconstVectorBlock(1).ptr(),
156 Teuchos::ScalarTraits<Scalar>::zero());
157 if (mass_matrix_is_identity_)
158 assign(adjoint_init_mv.ptr(), *dgdx);
159 else {
160 inargs.set_alpha(1.0);
161 inargs.set_beta(0.0);
162 RCP<LinearOpWithSolveBase<Scalar> > W;
163 if (adj_outargs.supports(MEB::OUT_ARG_W)) {
164 // Model supports W
165 W = adjoint_model_->create_W();
166 adj_outargs.set_W(W);
167 adjoint_model_->evalModel(inargs, adj_outargs);
168 adj_outargs.set_W(Teuchos::null);
169 }
170 else {
171 // Otherwise model must support a W_op and W factory
172 RCP<const LinearOpWithSolveFactoryBase<Scalar> > lowsfb =
173 adjoint_model_->get_W_factory();
174 TEUCHOS_TEST_FOR_EXCEPTION(
175 lowsfb == Teuchos::null, std::logic_error,
176 "Adjoint ME must support W out-arg or provide a W_factory for non-identity mass matrix");
177
178 // Compute W_op (and W_prec if supported)
179 RCP<LinearOpBase<Scalar> > W_op = adjoint_model_->create_W_op();
180 adj_outargs.set_W_op(W_op);
181 RCP<PreconditionerFactoryBase<Scalar> > prec_factory =
182 lowsfb->getPreconditionerFactory();
183 RCP<PreconditionerBase<Scalar> > W_prec;
184 if (prec_factory != Teuchos::null)
185 W_prec = prec_factory->createPrec();
186 else if (adj_outargs.supports(MEB::OUT_ARG_W_prec)) {
187 W_prec = adjoint_model_->create_W_prec();
188 adj_outargs.set_W_prec(W_prec);
189 }
190 adjoint_model_->evalModel(inargs, adj_outargs);
191 adj_outargs.set_W_op(Teuchos::null);
192 if (adj_outargs.supports(MEB::OUT_ARG_W_prec))
193 adj_outargs.set_W_prec(Teuchos::null);
194
195 // Create and initialize W
196 W = lowsfb->createOp();
197 if (W_prec != Teuchos::null) {
198 if (prec_factory != Teuchos::null)
199 prec_factory->initializePrec(
200 Thyra::defaultLinearOpSource<Scalar>(W_op), W_prec.get());
201 Thyra::initializePreconditionedOp<Scalar>(
202 *lowsfb, W_op, W_prec, W.ptr());
203 }
204 else
205 Thyra::initializeOp<Scalar>(*lowsfb, W_op, W.ptr());
206 }
207 TEUCHOS_TEST_FOR_EXCEPTION(
208 W == Teuchos::null, std::logic_error,
209 "A null W has been encountered in Tempus::IntegratorAdjointSensitivity::advanceTime!\n");
210 // Initialize adjoint_init_mv to zero before solve for linear solvers that
211 // use what is passed in as the initial guess
212 assign(adjoint_init_mv.ptr(), Teuchos::ScalarTraits<Scalar>::zero());
213 W->solve(Thyra::NOTRANS, *dgdx, adjoint_init_mv.ptr());
214 }
215
216 // Run sensitivity integrator and get solution
217 bool sens_status = true;
218 {
219 TEMPUS_FUNC_TIME_MONITOR_DIFF("Tempus::IntegratorAdjointSensitivity::advanceTime::adjoint", TEMPUS_AS_AT_ADJ);
221 const Scalar tinit = adjoint_integrator_->getTimeStepControl()->getInitTime();
222 adjoint_integrator_->initializeSolutionHistory(tinit, adjoint_init);
223 sens_status = adjoint_integrator_->advanceTime(timeFinal);
224 }
225 RCP<const SolutionHistory<Scalar> > adjoint_solution_history =
226 adjoint_integrator_->getSolutionHistory();
227
228 // Compute dg/dp at final time T
229 RCP<const VectorSpaceBase<Scalar> > p_space = model_->get_p_space(p_index_);
230 dgdp_ = createMembers(p_space, num_g);
231 if (g_depends_on_p_) {
232 MEB::DerivativeSupport dgdp_support =
233 outargs.supports(MEB::OUT_ARG_DgDp, g_index_, p_index_);
234 if (dgdp_support.supports(MEB::DERIV_MV_GRADIENT_FORM)) {
235 outargs.set_DgDp(g_index_, p_index_,
236 MEB::Derivative<Scalar>(dgdp_,
237 MEB::DERIV_MV_GRADIENT_FORM));
238 model_->evalModel(inargs, outargs);
239 }
240 else if (dgdp_support.supports(MEB::DERIV_MV_JACOBIAN_FORM)) {
241 const int num_p = p_space->dim();
242 RCP<MultiVectorBase<Scalar> > dgdp_trans =
243 createMembers(g_space, num_p);
244 outargs.set_DgDp(g_index_, p_index_,
245 MEB::Derivative<Scalar>(dgdp_trans,
246 MEB::DERIV_MV_JACOBIAN_FORM));
247 model_->evalModel(inargs, outargs);
248 Thyra::DetachedMultiVectorView<Scalar> dgdp_view(*dgdp_);
249 Thyra::DetachedMultiVectorView<Scalar> dgdp_trans_view(*dgdp_trans);
250 for (int i=0; i<num_p; ++i)
251 for (int j=0; j<num_g; ++j)
252 dgdp_view(i,j) = dgdp_trans_view(j,i);
253 }
254 else
255 TEUCHOS_TEST_FOR_EXCEPTION(true, std::logic_error,
256 "Invalid dg/dp support");
257 outargs.set_DgDp(g_index_, p_index_, MEB::Derivative<Scalar>());
258 }
259 else
260 assign(dgdp_.ptr(), Scalar(0.0));
261
262 // Add in initial condition term = (dx/dp^T(0))*(df/dx_dot^T(0))*y(0)
263 // If dxdp_init_ is null, assume it is zero
264 if (ic_depends_on_p_ && dxdp_init_ != Teuchos::null) {
265 RCP<const SolutionState<Scalar> > adjoint_state =
266 adjoint_solution_history->getCurrentState();
267 RCP<const VectorBase<Scalar> > adjoint_x =
268 rcp_dynamic_cast<const DPV>(adjoint_state->getX())->getVectorBlock(0);
269 RCP<const MultiVectorBase<Scalar> > adjoint_mv =
270 rcp_dynamic_cast<const DMVPV>(adjoint_x)->getMultiVector();
271 if (mass_matrix_is_identity_)
272 dxdp_init_->apply(Thyra::CONJTRANS, *adjoint_mv, dgdp_.ptr(), Scalar(1.0),
273 Scalar(1.0));
274 else {
275 inargs.set_t(initial_state->getTime());
276 inargs.set_x(initial_state->getX());
277 inargs.set_x_dot(initial_state->getXDot());
278 inargs.set_alpha(1.0);
279 inargs.set_beta(0.0);
280 RCP<LinearOpBase<Scalar> > W_op = adjoint_model_->create_W_op();
281 adj_outargs.set_W_op(W_op);
282 adjoint_model_->evalModel(inargs, adj_outargs);
283 adj_outargs.set_W_op(Teuchos::null);
284 RCP<MultiVectorBase<Scalar> > tmp = createMembers(x_space, num_g);
285 W_op->apply(Thyra::NOTRANS, *adjoint_mv, tmp.ptr(), Scalar(1.0),
286 Scalar(0.0));
287 dxdp_init_->apply(Thyra::CONJTRANS, *tmp, dgdp_.ptr(), Scalar(1.0),
288 Scalar(1.0));
289 }
290 }
291
292 // Add in model parameter term = \int_0^T( (df/dp^T(t)*y(t) )dt which
293 // is computed during the adjoint integration as an auxiliary integral
294 // (2nd block of the solution vector)
295 if (f_depends_on_p_) {
296 RCP<const SolutionState<Scalar> > adjoint_state =
297 adjoint_solution_history->getCurrentState();
298 RCP<const VectorBase<Scalar> > z =
299 rcp_dynamic_cast<const DPV>(adjoint_state->getX())->getVectorBlock(1);
300 RCP<const MultiVectorBase<Scalar> > z_mv =
301 rcp_dynamic_cast<const DMVPV>(z)->getMultiVector();
302 Thyra::V_VmV(dgdp_.ptr(), *dgdp_, *z_mv);
303 }
304
305 buildSolutionHistory(state_solution_history, adjoint_solution_history);
306
307 return state_status && sens_status;
308}
309
310template<class Scalar>
311Scalar
313getTime() const
314{
315 return solutionHistory_->getCurrentTime();
316}
317
318template<class Scalar>
319int
321getIndex() const
322{
323 return solutionHistory_->getCurrentIndex();
324}
325
326template<class Scalar>
327Status
329getStatus() const
330{
331 Status state_status = state_integrator_->getStatus();
332 Status sens_status = adjoint_integrator_->getStatus();
333 if (state_status == FAILED || sens_status == FAILED)
334 return FAILED;
335 if (state_status == WORKING || sens_status == WORKING)
336 return WORKING;
337 return PASSED;
338}
339
340template <class Scalar>
342 state_integrator_->setStatus(st);
343 adjoint_integrator_->setStatus(st);
344}
345
346template<class Scalar>
347Teuchos::RCP<Stepper<Scalar> >
349getStepper() const
350{
351 return state_integrator_->getStepper();
352}
353
354template<class Scalar>
355Teuchos::RCP<const SolutionHistory<Scalar> >
357getSolutionHistory() const
358{
359 return solutionHistory_;
360}
361
362template<class Scalar>
363Teuchos::RCP<const SolutionHistory<Scalar> >
366{
367 return state_integrator_->getSolutionHistory();
368}
369
370template<class Scalar>
371Teuchos::RCP<const SolutionHistory<Scalar> >
374{
375 return adjoint_integrator_->getSolutionHistory();
376}
377
378template<class Scalar>
379Teuchos::RCP<SolutionHistory<Scalar> >
382{
383 return solutionHistory_;
384}
385
386template<class Scalar>
387Teuchos::RCP<const TimeStepControl<Scalar> >
389getTimeStepControl() const
390{
391 return state_integrator_->getTimeStepControl();
392}
393
394template<class Scalar>
395Teuchos::RCP<TimeStepControl<Scalar> >
398{
399 return state_integrator_->getNonConstTimeStepControl();
400}
401
402template<class Scalar>
403Teuchos::RCP<TimeStepControl<Scalar> >
406{
407 return state_integrator_->getNonConstTimeStepControl();
408}
409
410template<class Scalar>
411Teuchos::RCP<TimeStepControl<Scalar> >
414{
415 return adjoint_integrator_->getNonConstTimeStepControl();
416}
417
418template<class Scalar>
421 Teuchos::RCP<const Thyra::VectorBase<Scalar> > x0,
422 Teuchos::RCP<const Thyra::VectorBase<Scalar> > xdot0,
423 Teuchos::RCP<const Thyra::VectorBase<Scalar> > xdotdot0,
424 Teuchos::RCP<const Thyra::MultiVectorBase<Scalar> > DxDp0,
425 Teuchos::RCP<const Thyra::MultiVectorBase<Scalar> > /* DxdotDp0 */,
426 Teuchos::RCP<const Thyra::MultiVectorBase<Scalar> > /* DxdotdotDp0 */)
427{
428 state_integrator_->initializeSolutionHistory(t0, x0, xdot0, xdotdot0);
429 dxdp_init_ = DxDp0;
430}
431
432template<class Scalar>
433Teuchos::RCP<IntegratorObserver<Scalar> >
436{
437 return state_integrator_->getObserver();
438}
439
440template<class Scalar>
442setObserver(Teuchos::RCP<IntegratorObserver<Scalar> > obs)
443{
444 state_integrator_->setObserver(obs);
445 // ETP 1/12/22 Disabling passing of the observer to the adjoint
446 // integrator to work around issues in Piro
447 //adjoint_integrator_->setObserver(obs);
448}
449
450template<class Scalar>
453{
454 state_integrator_->initialize();
455 adjoint_integrator_->initialize();
456}
457
458template<class Scalar>
459Teuchos::RCP<const Thyra::VectorBase<Scalar> >
461getX() const
462{
463 return state_integrator_->getX();
464}
465
466template<class Scalar>
467Teuchos::RCP<const Thyra::VectorBase<Scalar> >
469getXDot() const
470{
471 return state_integrator_->getXDot();
472}
473
474template<class Scalar>
475Teuchos::RCP<const Thyra::VectorBase<Scalar> >
477getXDotDot() const
478{
479 return state_integrator_->getXDotDot();
480}
481
482template<class Scalar>
483Teuchos::RCP<const Thyra::MultiVectorBase<Scalar> >
485getY() const
486{
487 using Teuchos::RCP;
488 using Teuchos::rcp_dynamic_cast;
489 typedef Thyra::DefaultProductVector<Scalar> DPV;
490 typedef Thyra::DefaultMultiVectorProductVector<Scalar> DMVPV;
491 RCP<const DPV> pv =
492 rcp_dynamic_cast<const DPV>(adjoint_integrator_->getX());
493 RCP<const DMVPV> mvpv =
494 rcp_dynamic_cast<const DMVPV>(pv->getVectorBlock(0));
495 return mvpv->getMultiVector();
496}
497
498template<class Scalar>
499Teuchos::RCP<const Thyra::MultiVectorBase<Scalar> >
501getYDot() const
502{
503 using Teuchos::RCP;
504 using Teuchos::rcp_dynamic_cast;
505 typedef Thyra::DefaultProductVector<Scalar> DPV;
506 typedef Thyra::DefaultMultiVectorProductVector<Scalar> DMVPV;
507 RCP<const DPV> pv =
508 rcp_dynamic_cast<const DPV>(adjoint_integrator_->getXDot());
509 RCP<const DMVPV> mvpv =
510 rcp_dynamic_cast<const DMVPV>(pv->getVectorBlock(0));
511 return mvpv->getMultiVector();
512}
513
514template<class Scalar>
515Teuchos::RCP<const Thyra::MultiVectorBase<Scalar> >
517getYDotDot() const
518{
519 using Teuchos::RCP;
520 using Teuchos::rcp_dynamic_cast;
521 typedef Thyra::DefaultProductVector<Scalar> DPV;
522 typedef Thyra::DefaultMultiVectorProductVector<Scalar> DMVPV;
523 RCP<const DPV> pv =
524 rcp_dynamic_cast<const DPV>(adjoint_integrator_->getXDotDot());
525 RCP<const DMVPV> mvpv =
526 rcp_dynamic_cast<const DMVPV>(pv->getVectorBlock(0));
527 return mvpv->getMultiVector();
528}
529
530template<class Scalar>
531Teuchos::RCP<const Thyra::MultiVectorBase<Scalar> >
533getDgDp() const
534{
535 return dgdp_;
536}
537
538template<class Scalar>
539std::string
541description() const
542{
543 std::string name = "Tempus::IntegratorAdjointSensitivity";
544 return(name);
545}
546
547template<class Scalar>
548void
551 Teuchos::FancyOStream &out,
552 const Teuchos::EVerbosityLevel verbLevel) const
553{
554 auto l_out = Teuchos::fancyOStream( out.getOStream() );
555 Teuchos::OSTab ostab(*l_out, 2, this->description());
556 l_out->setOutputToRootOnly(0);
557
558 *l_out << description() << "::describe" << std::endl;
559 state_integrator_->describe(*l_out, verbLevel);
560 adjoint_integrator_->describe(*l_out, verbLevel);
561}
562
563template<class Scalar>
566getStepMode() const
567{
568 return stepMode_;
569}
570
571template <class Scalar>
572Teuchos::RCP<AdjointAuxSensitivityModelEvaluator<Scalar> >
575 const Teuchos::RCP<Thyra::ModelEvaluator<Scalar> >& model,
576 const Teuchos::RCP<Thyra::ModelEvaluator<Scalar> >& adjoint_model,
577 const Teuchos::RCP<Teuchos::ParameterList>& inputPL)
578{
579 using Teuchos::ParameterList;
580 using Teuchos::RCP;
581 using Teuchos::rcp;
582
583 RCP<ParameterList> spl = Teuchos::parameterList();
584 if (inputPL != Teuchos::null)
585 {
586 *spl = inputPL->sublist("Sensitivities");
587 }
588 if (spl->isParameter("Response Depends on Parameters"))
589 spl->remove("Response Depends on Parameters");
590 if (spl->isParameter("Residual Depends on Parameters"))
591 spl->remove("Residual Depends on Parameters");
592 if (spl->isParameter("IC Depends on Parameters"))
593 spl->remove("IC Depends on Parameters");
594
595 const Scalar tinit = state_integrator_->getTimeStepControl()->getInitTime();
596 const Scalar tfinal = state_integrator_->getTimeStepControl()->getFinalTime();
597 return rcp(new AdjointAuxSensitivityModelEvaluator<Scalar>(
598 model, adjoint_model, tinit, tfinal, spl));
599}
600
601template<class Scalar>
602void
605 const Teuchos::RCP<const SolutionHistory<Scalar> >& state_solution_history,
606 const Teuchos::RCP<const SolutionHistory<Scalar> >& adjoint_solution_history)
607{
608 using Teuchos::RCP;
609 using Teuchos::rcp;
610 using Teuchos::rcp_dynamic_cast;
611 using Teuchos::ParameterList;
612 using Thyra::VectorBase;
614 using Thyra::VectorSpaceBase;
615 using Thyra::createMembers;
616 using Thyra::multiVectorProductVector;
617 using Thyra::assign;
618 typedef Thyra::DefaultProductVectorSpace<Scalar> DPVS;
619 typedef Thyra::DefaultProductVector<Scalar> DPV;
620
621 RCP<const VectorSpaceBase<Scalar> > x_space = model_->get_x_space();
622 RCP<const VectorSpaceBase<Scalar> > adjoint_space =
623 rcp_dynamic_cast<const DPVS>(adjoint_aux_model_->get_x_space())->getBlock(0);
624 Teuchos::Array< RCP<const VectorSpaceBase<Scalar> > > spaces(2);
625 spaces[0] = x_space;
626 spaces[1] = adjoint_space;
627 RCP<const DPVS > prod_space = Thyra::productVectorSpace(spaces());
628
629 int num_states = state_solution_history->getNumStates();
630 const Scalar t_init = state_integrator_->getTimeStepControl()->getInitTime();
631 const Scalar t_final = state_integrator_->getTime();
632 for (int i=0; i<num_states; ++i) {
633 RCP<const SolutionState<Scalar> > forward_state =
634 (*state_solution_history)[i];
635 RCP<const SolutionState<Scalar> > adjoint_state =
636 adjoint_solution_history->findState(t_final+t_init-forward_state->getTime());
637
638 // X
639 RCP<DPV> x = Thyra::defaultProductVector(prod_space);
640 RCP<const VectorBase<Scalar> > adjoint_x =
641 rcp_dynamic_cast<const DPV>(adjoint_state->getX())->getVectorBlock(0);
642 assign(x->getNonconstVectorBlock(0).ptr(), *(forward_state->getX()));
643 assign(x->getNonconstVectorBlock(1).ptr(), *(adjoint_x));
644 RCP<VectorBase<Scalar> > x_b = x;
645
646 // X-Dot
647 RCP<DPV> x_dot = Thyra::defaultProductVector(prod_space);
648 RCP<const VectorBase<Scalar> > adjoint_x_dot =
649 rcp_dynamic_cast<const DPV>(adjoint_state->getXDot())->getVectorBlock(0);
650 assign(x_dot->getNonconstVectorBlock(0).ptr(), *(forward_state->getXDot()));
651 assign(x_dot->getNonconstVectorBlock(1).ptr(), *(adjoint_x_dot));
652 RCP<VectorBase<Scalar> > x_dot_b = x_dot;
653
654 // X-Dot-Dot
655 RCP<DPV> x_dot_dot;
656 if (forward_state->getXDotDot() != Teuchos::null) {
657 x_dot_dot = Thyra::defaultProductVector(prod_space);
658 RCP<const VectorBase<Scalar> > adjoint_x_dot_dot =
659 rcp_dynamic_cast<const DPV>(
660 adjoint_state->getXDotDot())->getVectorBlock(0);
661 assign(x_dot_dot->getNonconstVectorBlock(0).ptr(),
662 *(forward_state->getXDotDot()));
663 assign(x_dot_dot->getNonconstVectorBlock(1).ptr(),
664 *(adjoint_x_dot_dot));
665 }
666 RCP<VectorBase<Scalar> > x_dot_dot_b = x_dot_dot;
667
668 RCP<SolutionState<Scalar> > prod_state = forward_state->clone();
669 prod_state->setX(x_b);
670 prod_state->setXDot(x_dot_b);
671 prod_state->setXDotDot(x_dot_dot_b);
672 prod_state->setPhysicsState(Teuchos::null);
673 solutionHistory_->addState(prod_state);
674 }
675}
676
678template<class Scalar>
679Teuchos::RCP<IntegratorAdjointSensitivity<Scalar> >
681 Teuchos::RCP<Teuchos::ParameterList> inputPL,
682 const Teuchos::RCP<Thyra::ModelEvaluator<Scalar> >& model
683 ,const Teuchos::RCP<Thyra::ModelEvaluator<Scalar> >& adjoint_model)
684{
685 // set the parameters
686 Teuchos::RCP<Teuchos::ParameterList> spl = Teuchos::parameterList();
687 if (inputPL != Teuchos::null)
688 *spl = inputPL->sublist("Sensitivities");
689
690 int p_index = spl->get<int>("Sensitivity Parameter Index", 0);
691 int g_index = spl->get<int>("Response Function Index", 0);
692 bool g_depends_on_p = spl->get<bool>("Response Depends on Parameters", true);
693 bool f_depends_on_p = spl->get<bool>("Residual Depends on Parameters", true);
694 bool ic_depends_on_p = spl->get<bool>("IC Depends on Parameters", true);
695 bool mass_matrix_is_identity = spl->get<bool>("Mass Matrix Is Identity", false);
696
697 auto state_integrator = createIntegratorBasic<Scalar>(inputPL, model);
698
699 // createAdjointModel
700 if (spl->isParameter("Response Depends on Parameters"))
701 spl->remove("Response Depends on Parameters");
702 if (spl->isParameter("Residual Depends on Parameters"))
703 spl->remove("Residual Depends on Parameters");
704 if (spl->isParameter("IC Depends on Parameters"))
705 spl->remove("IC Depends on Parameters");
706
707 const Scalar tinit = state_integrator->getTimeStepControl()->getInitTime();
708 const Scalar tfinal = state_integrator->getTimeStepControl()->getFinalTime();
709 //auto adjoint_model = Teuchos::rcp(new AdjointAuxSensitivityModelEvaluator<Scalar>(model, tfinal, spl));
710 //TODO: where is the adjoint ME coming from?
711
712 Teuchos::RCP<Thyra::ModelEvaluator<Scalar>> adjt_model = adjoint_model;
713 if (adjoint_model == Teuchos::null)
714 adjt_model = Thyra::implicitAdjointModelEvaluator(model);
715
716 auto adjoint_aux_model = Teuchos::rcp(new AdjointAuxSensitivityModelEvaluator<Scalar>(model, adjt_model, tinit, tfinal, spl));
717
718 // Create combined solution histories combining the forward and adjoint
719 // solutions. We do not include the auxiliary part from the adjoint solution.
720 auto integrator_name = inputPL->get<std::string>("Integrator Name");
721 auto integratorPL = Teuchos::sublist(inputPL, integrator_name, true);
722 auto shPL = Teuchos::sublist(integratorPL, "Solution History", true);
723 auto combined_solution_History = createSolutionHistoryPL<Scalar>(shPL);
724
725 auto adjoint_integrator = createIntegratorBasic<Scalar>(inputPL, adjoint_aux_model);
726
727 Teuchos::RCP<IntegratorAdjointSensitivity<Scalar>> integrator = Teuchos::rcp(new IntegratorAdjointSensitivity<Scalar>(
728 model, state_integrator, adjt_model, adjoint_aux_model, adjoint_integrator, combined_solution_History, p_index, g_index, g_depends_on_p,
729 f_depends_on_p, ic_depends_on_p, mass_matrix_is_identity));
730
731 return (integrator);
732}
733
735template<class Scalar>
736Teuchos::RCP<IntegratorAdjointSensitivity<Scalar> >
738{
739 Teuchos::RCP<IntegratorAdjointSensitivity<Scalar> > integrator =
740 Teuchos::rcp(new IntegratorAdjointSensitivity<Scalar>());
741 return(integrator);
742}
743
744} // namespace Tempus
745#endif // Tempus_IntegratorAdjointSensitivity_impl_hpp
Teuchos::RCP< const SolutionHistory< Scalar > > getSensSolutionHistory() const
virtual Teuchos::RCP< SolutionHistory< Scalar > > getNonConstSolutionHistory() override
Get the SolutionHistory.
virtual Teuchos::RCP< const Thyra::MultiVectorBase< Scalar > > getYDot() const
Get the current time derivative of the adjoint solution, ydot.
virtual void setObserver(Teuchos::RCP< IntegratorObserver< Scalar > > obs=Teuchos::null)
Set the Observer.
virtual Teuchos::RCP< const Thyra::VectorBase< Scalar > > getXDotDot() const
Get the current second time derivative of the solution, xdotdot.
virtual Scalar getTime() const override
Get current time.
virtual int getIndex() const override
Get current index.
virtual Teuchos::RCP< const SolutionHistory< Scalar > > getSolutionHistory() const override
Get the SolutionHistory.
void buildSolutionHistory(const Teuchos::RCP< const SolutionHistory< Scalar > > &state_solution_history, const Teuchos::RCP< const SolutionHistory< Scalar > > &adjoint_solution_history)
Teuchos::RCP< TimeStepControl< Scalar > > getSensNonConstTimeStepControl()
virtual void initialize()
Initializes the Integrator after set* function calls.
IntegratorAdjointSensitivity()
Constructor that requires a subsequent setParameterList, setStepper, and initialize calls.
virtual Teuchos::RCP< const Thyra::VectorBase< Scalar > > getX() const
Get the current solution, x.
virtual Teuchos::RCP< const Thyra::MultiVectorBase< Scalar > > getYDotDot() const
Get the current second time derivative of the adjoint solution, ydotdot.
virtual bool advanceTime()
Advance the solution to timeMax, and return true if successful.
virtual Teuchos::RCP< const Thyra::VectorBase< Scalar > > getXDot() const
Get the current time derivative of the solution, xdot.
virtual void setStatus(const Status st) override
Set Status.
SensitivityStepMode getStepMode() const
What mode the current time integration step is in.
void describe(Teuchos::FancyOStream &out, const Teuchos::EVerbosityLevel verbLevel) const override
virtual Teuchos::RCP< const TimeStepControl< Scalar > > getTimeStepControl() const override
Get the TimeStepControl.
virtual Teuchos::RCP< TimeStepControl< Scalar > > getNonConstTimeStepControl() override
virtual void initializeSolutionHistory(Scalar t0, Teuchos::RCP< const Thyra::VectorBase< Scalar > > x0, Teuchos::RCP< const Thyra::VectorBase< Scalar > > xdot0=Teuchos::null, Teuchos::RCP< const Thyra::VectorBase< Scalar > > xdotdot0=Teuchos::null, Teuchos::RCP< const Thyra::MultiVectorBase< Scalar > > DxDp0=Teuchos::null, Teuchos::RCP< const Thyra::MultiVectorBase< Scalar > > DxdotDp0=Teuchos::null, Teuchos::RCP< const Thyra::MultiVectorBase< Scalar > > DxdotdotDp0=Teuchos::null)
Set the initial state from Thyra::VectorBase(s)
virtual Teuchos::RCP< Stepper< Scalar > > getStepper() const override
Get the Stepper.
Teuchos::RCP< const SolutionHistory< Scalar > > getStateSolutionHistory() const
virtual Teuchos::RCP< const Thyra::MultiVectorBase< Scalar > > getDgDp() const
Return adjoint sensitivity stored in gradient format.
virtual Teuchos::RCP< const Thyra::MultiVectorBase< Scalar > > getY() const
Get the current adjoint solution, y.
virtual Teuchos::RCP< IntegratorObserver< Scalar > > getObserver()
Get the Observer.
Teuchos::RCP< TimeStepControl< Scalar > > getStateNonConstTimeStepControl()
Teuchos::RCP< AdjointAuxSensitivityModelEvaluator< Scalar > > createAdjointModel(const Teuchos::RCP< Thyra::ModelEvaluator< Scalar > > &model, const Teuchos::RCP< Thyra::ModelEvaluator< Scalar > > &adjoint_model, const Teuchos::RCP< Teuchos::ParameterList > &inputPL)
Status
Status for the Integrator, the Stepper and the SolutionState.
Teuchos::RCP< IntegratorAdjointSensitivity< Scalar > > createIntegratorAdjointSensitivity()
Nonmember constructor.
RCP< ImplicitAdjointModelEvaluator< Scalar > > implicitAdjointModelEvaluator(const RCP< const ModelEvaluator< Scalar > > &model)