9 #ifndef Tempus_StepperStaggeredForwardSensitivity_impl_hpp
10 #define Tempus_StepperStaggeredForwardSensitivity_impl_hpp
12 #include "Tempus_config.hpp"
16 #include "Teuchos_VerboseObjectParameterListHelpers.hpp"
18 #include "Thyra_DefaultMultiVectorProductVectorSpace.hpp"
19 #include "Thyra_DefaultMultiVectorProductVector.hpp"
24 template<
class Scalar>
class StepperFactory;
27 template<
class Scalar>
31 this->setParams(Teuchos::null, Teuchos::null);
35 template<
class Scalar>
38 const Teuchos::RCP<
const Thyra::ModelEvaluator<Scalar> >& appModel,
39 const Teuchos::RCP<Teuchos::ParameterList>& pList,
40 const Teuchos::RCP<Teuchos::ParameterList>& sens_pList)
43 this->setParams(pList, sens_pList);
44 this->setModel(appModel);
49 template<
class Scalar>
52 const Teuchos::RCP<
const Thyra::ModelEvaluator<Scalar> >& appModel)
56 using Teuchos::ParameterList;
59 Teuchos::RCP<Teuchos::ParameterList> spl = Teuchos::parameterList();
61 spl->remove(
"Reuse State Linear Solver");
62 spl->remove(
"Force W Update");
73 if (stateStepper_ == Teuchos::null)
74 stateStepper_ = sf->createStepper(stepperPL_, appModel);
76 stateStepper_->setModel(appModel);
78 if (sensitivityStepper_ == Teuchos::null)
79 sensitivityStepper_ = sf->createStepper(stepperPL_, fsa_model_);
81 sensitivityStepper_->setModel(fsa_model_);
85 template<
class Scalar>
88 const Teuchos::RCP<Thyra::ModelEvaluator<Scalar> >& appModel)
90 this->setModel(appModel);
94 template<
class Scalar>
95 Teuchos::RCP<const Thyra::ModelEvaluator<Scalar> >
99 return combined_fsa_model_;
103 template<
class Scalar>
106 Teuchos::RCP<Thyra::NonlinearSolverBase<Scalar> > solver)
108 stateStepper_->setSolver(solver);
109 sensitivityStepper_->setSolver(solver);
113 template<
class Scalar>
120 template<
class Scalar>
127 using Teuchos::rcp_dynamic_cast;
128 using Thyra::VectorBase;
129 using Thyra::MultiVectorBase;
131 using Thyra::createMember;
132 using Thyra::multiVectorProductVector;
133 using Thyra::multiVectorProductVectorSpace;
134 typedef Thyra::DefaultMultiVectorProductVector<Scalar> DMVPV;
135 typedef Thyra::DefaultMultiVectorProductVectorSpace<Scalar> DMVPVS;
141 if (stateSolutionHistory_ == Teuchos::null) {
142 RCP<Teuchos::ParameterList> shPL =
146 RCP<SolutionState<Scalar> > state =
solutionHistory->getCurrentState();
147 RCP<DMVPV> X, XDot, XDotDot;
148 X = rcp_dynamic_cast<DMVPV>(state->getX(),
true);
149 XDot = rcp_dynamic_cast<DMVPV>(state->getXDot(),
true);
150 if (state->getXDotDot() != Teuchos::null)
151 XDotDot = rcp_dynamic_cast<DMVPV>(state->getXDotDot(),
true);
155 RCP<VectorBase<Scalar> > x, xdot, xdotdot;
156 x = X->getNonconstMultiVector()->col(0);
157 xdot = XDot->getNonconstMultiVector()->col(0);
158 if (XDotDot != Teuchos::null)
159 xdotdot = XDotDot->getNonconstMultiVector()->col(0);
162 RCP<SolutionState<Scalar> > state_state =
165 state->getStepperState()->clone()));
167 stateSolutionHistory_->addState(state_state);
169 const int num_param = X->getMultiVector()->domain()->dim()-1;
170 TEUCHOS_ASSERT(num_param > 0);
171 const Teuchos::Range1D rng(1,num_param);
174 RCP<MultiVectorBase<Scalar> > dxdp, dxdotdp, dxdotdotdp;
175 dxdp = X->getNonconstMultiVector()->subView(rng);
176 dxdotdp = XDot->getNonconstMultiVector()->subView(rng);
177 if (XDotDot != Teuchos::null)
178 dxdotdotdp = XDotDot->getNonconstMultiVector()->subView(rng);
181 RCP<VectorBase<Scalar> > dxdp_vec, dxdotdp_vec, dxdotdotdp_vec;
182 RCP<const DMVPVS> prod_space =
183 multiVectorProductVectorSpace(X->getMultiVector()->range(), num_param);
184 dxdp_vec = multiVectorProductVector(prod_space, dxdp);
185 dxdotdp_vec = multiVectorProductVector(prod_space, dxdotdp);
186 if (XDotDot != Teuchos::null)
187 dxdotdotdp_vec = multiVectorProductVector(prod_space, dxdotdotdp);
190 RCP<SolutionState<Scalar> > sens_state =
192 dxdp_vec, dxdotdp_vec, dxdotdotdp_vec,
193 state->getStepperState()->clone()));
195 sensSolutionHistory_->addState(sens_state);
199 RCP<SolutionState<Scalar> > prod_state =
solutionHistory->getWorkingState();
200 RCP<DMVPV> X, XDot, XDotDot;
201 X = rcp_dynamic_cast<DMVPV>(prod_state->getX(),
true);
202 XDot = rcp_dynamic_cast<DMVPV>(prod_state->getXDot(),
true);
203 if (prod_state->getXDotDot() != Teuchos::null)
204 XDotDot = rcp_dynamic_cast<DMVPV>(prod_state->getXDotDot(),
true);
207 stateSolutionHistory_->initWorkingState();
208 RCP<SolutionState<Scalar> > state = stateSolutionHistory_->getWorkingState();
209 state->getMetaData()->copy(prod_state->getMetaData());
210 stateStepper_->takeStep(stateSolutionHistory_);
213 assign(X->getNonconstMultiVector()->col(0).ptr(), *(state->getX()));
214 assign(XDot->getNonconstMultiVector()->col(0).ptr(), *(state->getXDot()));
215 if (XDotDot != Teuchos::null)
216 assign(XDotDot->getNonconstMultiVector()->col(0).ptr(),
217 *(state->getXDotDot()));
218 prod_state->setOrder(state->getOrder());
225 stateSolutionHistory_->promoteWorkingState();
228 fsa_model_->setForwardSolutionHistory(stateSolutionHistory_);
229 if (reuse_solver_ && stateStepper_->getSolver() != Teuchos::null)
230 fsa_model_->setSolver(stateStepper_->getSolver(), force_W_update_);
233 sensSolutionHistory_->initWorkingState();
234 RCP<SolutionState<Scalar> > sens_state =
235 sensSolutionHistory_->getWorkingState();
236 sens_state->getMetaData()->copy(prod_state->getMetaData());
237 sensitivityStepper_->takeStep(sensSolutionHistory_);
240 RCP<const MultiVectorBase<Scalar> > dxdp =
241 rcp_dynamic_cast<
const DMVPV>(sens_state->getX(),
true)->getMultiVector();
242 const int num_param = dxdp->domain()->dim();
243 const Teuchos::Range1D rng(1,num_param);
244 assign(X->getNonconstMultiVector()->subView(rng).ptr(), *dxdp);
245 RCP<const MultiVectorBase<Scalar> > dxdotdp =
246 rcp_dynamic_cast<
const DMVPV>(sens_state->getXDot(),
true)->getMultiVector();
247 assign(XDot->getNonconstMultiVector()->subView(rng).ptr(), *dxdotdp);
248 if (sens_state->getXDotDot() != Teuchos::null) {
249 RCP<const MultiVectorBase<Scalar> > dxdotdotdp =
250 rcp_dynamic_cast<
const DMVPV>(sens_state->getXDotDot(),
true)->getMultiVector();
251 assign(XDotDot->getNonconstMultiVector()->subView(rng).ptr(), *dxdotdotdp);
253 prod_state->setOrder(std::min(state->getOrder(), sens_state->getOrder()));
260 sensSolutionHistory_->promoteWorkingState();
266 template<
class Scalar>
267 Teuchos::RCP<Tempus::StepperState<Scalar> >
273 Teuchos::RCP<Tempus::StepperState<Scalar> > stepperState =
279 template<
class Scalar>
282 Teuchos::FancyOStream &out,
283 const Teuchos::EVerbosityLevel verbLevel)
const
285 out << description() <<
"::describe:" << std::endl;
286 stateStepper_->describe(out, verbLevel);
288 sensitivityStepper_->describe(out, verbLevel);
292 template <
class Scalar>
295 Teuchos::RCP<Teuchos::ParameterList>
const& pList)
297 if (pList == Teuchos::null) {
299 if (this->stepperPL_ == Teuchos::null) this->stepperPL_ =
300 Teuchos::rcp_const_cast<Teuchos::ParameterList>(this->getValidParameters());
302 this->stepperPL_ = pList;
309 template<
class Scalar>
310 Teuchos::RCP<const Teuchos::ParameterList>
314 return stateStepper_->getValidParameters();
318 template <
class Scalar>
319 Teuchos::RCP<Teuchos::ParameterList>
327 template <
class Scalar>
328 Teuchos::RCP<Teuchos::ParameterList>
332 Teuchos::RCP<Teuchos::ParameterList> temp_plist = stepperPL_;
333 stepperPL_ = Teuchos::null;
338 template <
class Scalar>
341 Teuchos::RCP<Teuchos::ParameterList>
const& pList,
342 Teuchos::RCP<Teuchos::ParameterList>
const& spList)
344 if (pList == Teuchos::null)
345 stepperPL_ = Teuchos::rcp_const_cast<Teuchos::ParameterList>(this->getValidParameters());
349 if (spList == Teuchos::null)
350 sensPL_ = Teuchos::parameterList();
354 reuse_solver_ = sensPL_->get(
"Reuse State Linear Solver",
false);
355 force_W_update_ = sensPL_->get(
"Force W Update",
true);
362 template <
class Scalar>
363 Teuchos::RCP<const Thyra::VectorSpaceBase<Scalar> >
368 using Teuchos::rcp_dynamic_cast;
369 typedef Thyra::DefaultMultiVectorProductVectorSpace<Scalar> DMVPVS;
371 RCP<const Thyra::VectorSpaceBase<Scalar> > x_space =
372 stateStepper_->getModel()->get_x_space();
373 RCP<const DMVPVS> dxdp_space =
374 rcp_dynamic_cast<
const DMVPVS>(sensitivityStepper_->getModel()->get_x_space(),
true);
375 const int num_param = dxdp_space->numBlocks();
376 RCP<const Thyra::VectorSpaceBase<Scalar> > prod_space =
377 multiVectorProductVectorSpace(x_space, num_param+1);
383 #endif // Tempus_StepperStaggeredForwardSensitivity_impl_hpp
virtual void initialize()
Initialize during construction and after changing input parameters.
Teuchos::RCP< Teuchos::ParameterList > getNonconstParameterList()
virtual void setNonConstModel(const Teuchos::RCP< Thyra::ModelEvaluator< Scalar > > &appModel)
virtual Teuchos::RCP< const Thyra::ModelEvaluator< Scalar > > getModel()
virtual void takeStep(const Teuchos::RCP< SolutionHistory< Scalar > > &solutionHistory)
Take the specified timestep, dt, and return true if successful.
StepperStaggeredForwardSensitivity()
Default constructor.
Teuchos::RCP< SensitivityModelEvaluatorBase< Scalar > > wrapStaggeredFSAModelEvaluator(const Teuchos::RCP< const Thyra::ModelEvaluator< Scalar > > &model, const Teuchos::RCP< const Teuchos::ParameterList > &pList=Teuchos::null)
Teuchos::RCP< const Thyra::VectorSpaceBase< Scalar > > get_x_space() const
virtual Teuchos::RCP< Tempus::StepperState< Scalar > > getDefaultStepperState()
Get a default (initial) StepperState.
virtual void setModel(const Teuchos::RCP< const Thyra::ModelEvaluator< Scalar > > &appModel)
StepperState is a simple class to hold state information about the stepper.
Teuchos::RCP< Teuchos::ParameterList > unsetParameterList()
Teuchos::RCP< SolutionHistory< Scalar > > solutionHistory(Teuchos::RCP< Teuchos::ParameterList > pList=Teuchos::null)
Nonmember constructor.
SolutionHistory is basically a container of SolutionStates. SolutionHistory maintains a collection of...
Teuchos::RCP< const Teuchos::ParameterList > getValidParameters() const
Teuchos::RCP< SensitivityModelEvaluatorBase< Scalar > > wrapCombinedFSAModelEvaluator(const Teuchos::RCP< const Thyra::ModelEvaluator< Scalar > > &model, const Teuchos::RCP< const Teuchos::ParameterList > &pList=Teuchos::null)
void setParams(const Teuchos::RCP< Teuchos::ParameterList > &pl, const Teuchos::RCP< Teuchos::ParameterList > &spl)
virtual void setSolver(Teuchos::RCP< Thyra::NonlinearSolverBase< Scalar > > solver=Teuchos::null)
Set solver.
Solution state for integrators and steppers. SolutionState contains the metadata for solutions and th...
void setParameterList(const Teuchos::RCP< Teuchos::ParameterList > &pl)
virtual void describe(Teuchos::FancyOStream &out, const Teuchos::EVerbosityLevel verbLevel) const