273 using Teuchos::rcp_dynamic_cast;
274 using Thyra::VectorSpaceBase;
276 using Thyra::createMember;
277 typedef Thyra::DefaultMultiVectorProductVector<Scalar> DMVPV;
282 RCP< const VectorSpaceBase<Scalar> > space = sens_model_->get_x_space();
283 RCP<DMVPV> X = rcp_dynamic_cast<DMVPV>(createMember(space));
284 RCP<DMVPV> Xdot = rcp_dynamic_cast<DMVPV>(createMember(space));
285 RCP<DMVPV> Xdotdot = rcp_dynamic_cast<DMVPV>(createMember(space));
286 const Scalar zero = Teuchos::ScalarTraits<Scalar>::zero();
289 if (DxDp0 == Teuchos::null)
290 assign(X->getNonconstMultiVector().ptr(), zero);
292 assign(X->getNonconstMultiVector().ptr(), *DxDp0);
295 if (DxdotDp0 == Teuchos::null)
296 assign(Xdot->getNonconstMultiVector().ptr(), zero);
298 assign(Xdot->getNonconstMultiVector().ptr(), *DxdotDp0);
301 if (DxdotDp0 == Teuchos::null)
302 assign(Xdotdot->getNonconstMultiVector().ptr(), zero);
304 assign(Xdotdot->getNonconstMultiVector().ptr(), *DxdotdotDp0);
306 state_integrator_->initializeSolutionHistory(t0, x0, xdot0, xdotdot0);
307 sens_integrator_->initializeSolutionHistory(t0, X, Xdot, Xdotdot);
470 using Teuchos::rcp_dynamic_cast;
471 using Teuchos::ParameterList;
474 using Thyra::VectorSpaceBase;
475 using Thyra::createMembers;
476 using Thyra::multiVectorProductVector;
478 typedef Thyra::DefaultMultiVectorProductVector<Scalar> DMVPV;
484 RCP<ParameterList> shPL;
486 solutionHistory_ = createSolutionHistoryPL<Scalar>(shPL);
488 const int num_param =
489 rcp_dynamic_cast<const DMVPV>(sens_integrator_->getX())->getMultiVector()->domain()->dim();
490 RCP<const VectorSpaceBase<Scalar> > x_space = model_->get_x_space();
491 RCP<const Thyra::DefaultMultiVectorProductVectorSpace<Scalar> > prod_space =
492 Thyra::multiVectorProductVectorSpace(x_space, num_param+1);
493 const Teuchos::Range1D rng(1,num_param);
494 const Scalar zero = Teuchos::ScalarTraits<Scalar>::zero();
496 RCP<const SolutionHistory<Scalar> > state_solution_history =
497 state_integrator_->getSolutionHistory();
498 int num_states = state_solution_history->getNumStates();
499 for (
int i=0; i<num_states; ++i) {
500 RCP<const SolutionState<Scalar> > state = (*state_solution_history)[i];
503 RCP< MultiVectorBase<Scalar> > x_mv =
504 createMembers(x_space, num_param+1);
505 assign(x_mv->col(0).ptr(), *(state->getX()));
506 assign(x_mv->subView(rng).ptr(), zero);
507 RCP<VectorBase<Scalar> > x = multiVectorProductVector(prod_space, x_mv);
510 RCP<VectorBase<Scalar> > x_dot;
511 if (state->getXDot() != Teuchos::null) {
512 RCP< MultiVectorBase<Scalar> > x_dot_mv =
513 createMembers(x_space, num_param+1);
514 assign(x_dot_mv->col(0).ptr(), *(state->getXDot()));
515 assign(x_dot_mv->subView(rng).ptr(), zero);
516 x_dot = multiVectorProductVector(prod_space, x_dot_mv);
520 RCP<VectorBase<Scalar> > x_dot_dot;
521 if (state->getXDotDot() != Teuchos::null) {
522 RCP< MultiVectorBase<Scalar> > x_dot_dot_mv =
523 createMembers(x_space, num_param+1);
524 assign(x_dot_dot_mv->col(0).ptr(), *(state->getXDotDot()));
525 assign(x_dot_dot_mv->subView(rng).ptr(), zero);
526 x_dot_dot = multiVectorProductVector(prod_space, x_dot_dot_mv);
529 RCP<SolutionState<Scalar> > prod_state = state->clone();
531 prod_state->setXDot(x_dot);
532 prod_state->setXDotDot(x_dot_dot);
533 solutionHistory_->addState(prod_state);
536 RCP<const VectorBase<Scalar> > frozen_x =
537 state_solution_history->getCurrentState()->getX();
538 RCP<const VectorBase<Scalar> > frozen_x_dot =
539 state_solution_history->getCurrentState()->getXDot();
540 RCP<const VectorBase<Scalar> > frozen_x_dot_dot =
541 state_solution_history->getCurrentState()->getXDotDot();
542 RCP<const SolutionHistory<Scalar> > sens_solution_history =
543 sens_integrator_->getSolutionHistory();
544 num_states = sens_solution_history->getNumStates();
545 for (
int i=0; i<num_states; ++i) {
546 RCP<const SolutionState<Scalar> > state = (*sens_solution_history)[i];
549 RCP< MultiVectorBase<Scalar> > x_mv =
550 createMembers(x_space, num_param+1);
551 RCP<const MultiVectorBase<Scalar> > dxdp =
552 rcp_dynamic_cast<const DMVPV>(state->getX())->getMultiVector();
553 assign(x_mv->col(0).ptr(), *(frozen_x));
554 assign(x_mv->subView(rng).ptr(), *dxdp);
555 RCP<VectorBase<Scalar> > x = multiVectorProductVector(prod_space, x_mv);
558 RCP<VectorBase<Scalar> > x_dot;
559 if (state->getXDot() != Teuchos::null) {
560 RCP< MultiVectorBase<Scalar> > x_dot_mv =
561 createMembers(x_space, num_param+1);
562 RCP<const MultiVectorBase<Scalar> > dxdotdp =
563 rcp_dynamic_cast<const DMVPV>(state->getXDot())->getMultiVector();
564 assign(x_dot_mv->col(0).ptr(), *(frozen_x_dot));
565 assign(x_dot_mv->subView(rng).ptr(), *dxdotdp);
566 x_dot = multiVectorProductVector(prod_space, x_dot_mv);
570 RCP<VectorBase<Scalar> > x_dot_dot;
571 if (state->getXDotDot() != Teuchos::null) {
572 RCP< MultiVectorBase<Scalar> > x_dot_dot_mv =
573 createMembers(x_space, num_param+1);
574 RCP<const MultiVectorBase<Scalar> > dxdotdotdp =
575 rcp_dynamic_cast<const DMVPV>(state->getXDotDot())->getMultiVector();
576 assign(x_dot_dot_mv->col(0).ptr(), *(frozen_x_dot_dot));
577 assign(x_dot_dot_mv->subView(rng).ptr(), *dxdotdotdp);
578 x_dot_dot = multiVectorProductVector(prod_space, x_dot_dot_mv);
581 RCP<SolutionState<Scalar> > prod_state = state->clone();
583 prod_state->setXDot(x_dot);
584 prod_state->setXDotDot(x_dot_dot);
585 solutionHistory_->addState(prod_state,
false);
593 Teuchos::RCP<Teuchos::ParameterList> pList,
599 auto fwd_integrator = createIntegratorBasic<Scalar>(pList, model);
600 Teuchos::RCP<SensitivityModelEvaluatorBase<Scalar> > sens_model;
601 Teuchos::RCP<IntegratorBasic<Scalar> > sens_integrator;
604 Teuchos::RCP<Teuchos::ParameterList> pl = Teuchos::rcp(
new Teuchos::ParameterList);
605 Teuchos::RCP<const Teuchos::ParameterList> integrator_pl = fwd_integrator->getValidParameters();
606 Teuchos::RCP<const Teuchos::ParameterList> sensitivity_pl =
608 pl->setParameters(*integrator_pl);
609 pl->sublist(
"Sensitivities").setParameters(*sensitivity_pl);
610 pl->sublist(
"Sensitivities").set(
"Reuse State Linear Solver",
false);
611 pl->sublist(
"Sensitivities").set(
"Force W Update",
false);
612 pl->sublist(
"Sensitivities").set(
"Cache Matrices",
false);
613 pList->setParametersNotAlreadySet(*pl);
616 bool reuse_solver = pList->sublist(
"Sensitivities").get(
"Reuse State Linear Solver",
false);
617 bool force_W_update = pList->sublist(
"Sensitivities").get(
"Force W Update",
false);
618 bool cache_matrices = pList->sublist(
"Sensitivities").get(
"Cache Matrices",
false);
621 Teuchos::RCP<Teuchos::ParameterList> pl = Teuchos::parameterList();
622 if (pList!= Teuchos::null)
624 *pl = pList->sublist(
"Sensitivities");
626 pl->remove(
"Reuse State Linear Solver");
627 pl->remove(
"Force W Update");
628 pl->remove(
"Cache Matrices");
630 model, sens_residual_model, sens_solve_model, cache_matrices, pl);
631 sens_integrator = createIntegratorBasic<Scalar>(pList, sens_model);
634 Teuchos::RCP<Tempus::IntegratorPseudoTransientForwardSensitivity<Scalar>> integrator =