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->setStepperType(
"StaggeredForwardSensitivity");
32 this->setParams(Teuchos::null, Teuchos::null);
36 template<
class Scalar>
39 const Teuchos::RCP<
const Thyra::ModelEvaluator<Scalar> >& appModel,
40 const Teuchos::RCP<Teuchos::ParameterList>& pList,
41 const Teuchos::RCP<Teuchos::ParameterList>& sens_pList)
44 this->setStepperType(
"StaggeredForwardSensitivity");
45 this->setParams(pList, sens_pList);
46 this->setModel(appModel);
51 template<
class Scalar>
54 const Teuchos::RCP<
const Thyra::ModelEvaluator<Scalar> >& appModel)
58 using Teuchos::ParameterList;
61 Teuchos::RCP<Teuchos::ParameterList> spl = Teuchos::parameterList();
63 spl->remove(
"Reuse State Linear Solver");
64 spl->remove(
"Force W Update");
75 if (stateStepper_ == Teuchos::null)
76 stateStepper_ = sf->createStepper(stepperPL_, appModel);
78 stateStepper_->setModel(appModel);
80 if (sensitivityStepper_ == Teuchos::null)
81 sensitivityStepper_ = sf->createStepper(stepperPL_, fsa_model_);
83 sensitivityStepper_->setModel(fsa_model_);
85 this->isInitialized_ =
false;
89 template<
class Scalar>
90 Teuchos::RCP<const Thyra::ModelEvaluator<Scalar> >
94 return combined_fsa_model_;
98 template<
class Scalar>
101 Teuchos::RCP<Thyra::NonlinearSolverBase<Scalar> > solver)
103 stateStepper_->setSolver(solver);
104 sensitivityStepper_->setSolver(solver);
106 this->isInitialized_ =
false;
110 template<
class Scalar>
115 this->checkInitialized();
119 using Teuchos::rcp_dynamic_cast;
120 using Thyra::VectorBase;
121 using Thyra::MultiVectorBase;
123 using Thyra::createMember;
124 using Thyra::multiVectorProductVector;
125 using Thyra::multiVectorProductVectorSpace;
126 typedef Thyra::DefaultMultiVectorProductVector<Scalar> DMVPV;
127 typedef Thyra::DefaultMultiVectorProductVectorSpace<Scalar> DMVPVS;
133 if (stateSolutionHistory_ == Teuchos::null) {
134 RCP<Teuchos::ParameterList> shPL =
135 solutionHistory->getNonconstParameterList();
138 RCP<SolutionState<Scalar> > state = solutionHistory->getCurrentState();
139 RCP<DMVPV> X, XDot, XDotDot;
140 X = rcp_dynamic_cast<DMVPV>(state->getX(),
true);
141 XDot = rcp_dynamic_cast<DMVPV>(state->getXDot(),
true);
142 if (state->getXDotDot() != Teuchos::null)
143 XDotDot = rcp_dynamic_cast<DMVPV>(state->getXDotDot(),
true);
147 RCP<VectorBase<Scalar> > x, xdot, xdotdot;
148 x = X->getNonconstMultiVector()->col(0);
149 xdot = XDot->getNonconstMultiVector()->col(0);
150 if (XDotDot != Teuchos::null)
151 xdotdot = XDotDot->getNonconstMultiVector()->col(0);
154 RCP<SolutionState<Scalar> > state_state = state->clone();
155 state_state->setX(x);
156 state_state->setXDot(xdot);
157 state_state->setXDotDot(xdotdot);
159 stateSolutionHistory_->addState(state_state);
161 const int num_param = X->getMultiVector()->domain()->dim()-1;
162 TEUCHOS_ASSERT(num_param > 0);
163 const Teuchos::Range1D rng(1,num_param);
166 RCP<MultiVectorBase<Scalar> > dxdp, dxdotdp, dxdotdotdp;
167 dxdp = X->getNonconstMultiVector()->subView(rng);
168 dxdotdp = XDot->getNonconstMultiVector()->subView(rng);
169 if (XDotDot != Teuchos::null)
170 dxdotdotdp = XDotDot->getNonconstMultiVector()->subView(rng);
173 RCP<VectorBase<Scalar> > dxdp_vec, dxdotdp_vec, dxdotdotdp_vec;
174 RCP<const DMVPVS> prod_space =
175 multiVectorProductVectorSpace(X->getMultiVector()->range(), num_param);
176 dxdp_vec = multiVectorProductVector(prod_space, dxdp);
177 dxdotdp_vec = multiVectorProductVector(prod_space, dxdotdp);
178 if (XDotDot != Teuchos::null)
179 dxdotdotdp_vec = multiVectorProductVector(prod_space, dxdotdotdp);
182 RCP<SolutionState<Scalar> > sens_state = state->clone();
183 sens_state->setX(dxdp_vec);
184 sens_state->setXDot(dxdotdp_vec);
185 sens_state->setXDotDot(dxdotdotdp_vec);
187 sensSolutionHistory_->addState(sens_state);
191 RCP<SolutionState<Scalar> > prod_state = solutionHistory->getWorkingState();
192 RCP<DMVPV> X, XDot, XDotDot;
193 X = rcp_dynamic_cast<DMVPV>(prod_state->getX(),
true);
194 XDot = rcp_dynamic_cast<DMVPV>(prod_state->getXDot(),
true);
195 if (prod_state->getXDotDot() != Teuchos::null)
196 XDotDot = rcp_dynamic_cast<DMVPV>(prod_state->getXDotDot(),
true);
199 stateSolutionHistory_->initWorkingState();
200 RCP<SolutionState<Scalar> > state = stateSolutionHistory_->getWorkingState();
201 state->getMetaData()->copy(prod_state->getMetaData());
202 stateStepper_->takeStep(stateSolutionHistory_);
205 assign(X->getNonconstMultiVector()->col(0).ptr(), *(state->getX()));
206 assign(XDot->getNonconstMultiVector()->col(0).ptr(), *(state->getXDot()));
207 if (XDotDot != Teuchos::null)
208 assign(XDotDot->getNonconstMultiVector()->col(0).ptr(),
209 *(state->getXDotDot()));
210 prod_state->setOrder(state->getOrder());
217 stateSolutionHistory_->promoteWorkingState();
220 fsa_model_->setForwardSolutionHistory(stateSolutionHistory_);
221 if (reuse_solver_ && stateStepper_->getSolver() != Teuchos::null)
222 fsa_model_->setSolver(stateStepper_->getSolver(), force_W_update_);
225 sensSolutionHistory_->initWorkingState();
226 RCP<SolutionState<Scalar> > sens_state =
227 sensSolutionHistory_->getWorkingState();
228 sens_state->getMetaData()->copy(prod_state->getMetaData());
229 sensitivityStepper_->takeStep(sensSolutionHistory_);
232 RCP<const MultiVectorBase<Scalar> > dxdp =
233 rcp_dynamic_cast<
const DMVPV>(sens_state->getX(),
true)->getMultiVector();
234 const int num_param = dxdp->domain()->dim();
235 const Teuchos::Range1D rng(1,num_param);
236 assign(X->getNonconstMultiVector()->subView(rng).ptr(), *dxdp);
237 RCP<const MultiVectorBase<Scalar> > dxdotdp =
238 rcp_dynamic_cast<
const DMVPV>(sens_state->getXDot(),
true)->getMultiVector();
239 assign(XDot->getNonconstMultiVector()->subView(rng).ptr(), *dxdotdp);
240 if (sens_state->getXDotDot() != Teuchos::null) {
241 RCP<const MultiVectorBase<Scalar> > dxdotdotdp =
242 rcp_dynamic_cast<
const DMVPV>(sens_state->getXDotDot(),
true)->getMultiVector();
243 assign(XDotDot->getNonconstMultiVector()->subView(rng).ptr(), *dxdotdotdp);
245 prod_state->setOrder(std::min(state->getOrder(), sens_state->getOrder()));
252 sensSolutionHistory_->promoteWorkingState();
258 template<
class Scalar>
259 Teuchos::RCP<Tempus::StepperState<Scalar> >
265 Teuchos::RCP<Tempus::StepperState<Scalar> > stepperState =
271 template<
class Scalar>
274 Teuchos::FancyOStream &out,
275 const Teuchos::EVerbosityLevel verbLevel)
const
280 out <<
"--- StepperStaggeredForwardSensitivity ---\n";
281 out <<
" stateStepper_ = " << stateStepper_ <<std::endl;
282 out <<
" sensitivityStepper_ = " << sensitivityStepper_ <<std::endl;
283 out <<
" combined_fsa_model_ = " << combined_fsa_model_ <<std::endl;
284 out <<
" fsa_model_ = " << fsa_model_ <<std::endl;
285 out <<
" stateSolutionHistory_ = " << stateSolutionHistory_ <<std::endl;
286 out <<
" sensSolutionHistory_ = " << sensSolutionHistory_ <<std::endl;
287 out <<
"------------------------------------------" << std::endl;
289 out << description() <<
"::describe:" << std::endl;
290 stateStepper_->describe(out, verbLevel);
292 sensitivityStepper_->describe(out, verbLevel);
296 template<
class Scalar>
299 bool isValidSetup =
true;
303 if (stateStepper_ == Teuchos::null) {
304 isValidSetup =
false;
305 out <<
"The state stepper is not set!\n";
308 if (sensitivityStepper_ == Teuchos::null) {
309 isValidSetup =
false;
310 out <<
"The sensitivity stepper is not set!\n";
313 if (combined_fsa_model_ == Teuchos::null) {
314 isValidSetup =
false;
315 out <<
"The combined FSA model is not set!\n";
318 if (fsa_model_ == Teuchos::null) {
319 isValidSetup =
false;
320 out <<
"The FSA model is not set!\n";
327 template <
class Scalar>
330 Teuchos::RCP<Teuchos::ParameterList>
const& pList)
332 if (pList == Teuchos::null) {
334 if (this->stepperPL_ == Teuchos::null) this->stepperPL_ =
335 Teuchos::rcp_const_cast<Teuchos::ParameterList>(this->getValidParameters());
337 this->stepperPL_ = pList;
344 template<
class Scalar>
345 Teuchos::RCP<const Teuchos::ParameterList>
349 return stateStepper_->getValidParameters();
353 template <
class Scalar>
354 Teuchos::RCP<Teuchos::ParameterList>
362 template <
class Scalar>
363 Teuchos::RCP<Teuchos::ParameterList>
367 Teuchos::RCP<Teuchos::ParameterList> temp_plist = stepperPL_;
368 stepperPL_ = Teuchos::null;
373 template <
class Scalar>
376 Teuchos::RCP<Teuchos::ParameterList>
const& pList,
377 Teuchos::RCP<Teuchos::ParameterList>
const& spList)
379 if (pList == Teuchos::null)
380 stepperPL_ = Teuchos::rcp_const_cast<Teuchos::ParameterList>(this->getValidParameters());
384 if (spList == Teuchos::null)
385 sensPL_ = Teuchos::parameterList();
389 reuse_solver_ = sensPL_->get(
"Reuse State Linear Solver",
false);
390 force_W_update_ = sensPL_->get(
"Force W Update",
true);
397 template <
class Scalar>
398 Teuchos::RCP<const Thyra::VectorSpaceBase<Scalar> >
403 using Teuchos::rcp_dynamic_cast;
404 typedef Thyra::DefaultMultiVectorProductVectorSpace<Scalar> DMVPVS;
406 RCP<const Thyra::VectorSpaceBase<Scalar> > x_space =
407 stateStepper_->getModel()->get_x_space();
408 RCP<const DMVPVS> dxdp_space =
409 rcp_dynamic_cast<
const DMVPVS>(sensitivityStepper_->getModel()->get_x_space(),
true);
410 const int num_param = dxdp_space->numBlocks();
411 RCP<const Thyra::VectorSpaceBase<Scalar> > prod_space =
412 multiVectorProductVectorSpace(x_space, num_param+1);
418 #endif // Tempus_StepperStaggeredForwardSensitivity_impl_hpp
Teuchos::RCP< Teuchos::ParameterList > getNonconstParameterList()
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.
Thyra Base interface for time steppers.
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()
virtual void describe(Teuchos::FancyOStream &out, const Teuchos::EVerbosityLevel verbLevel) const
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 bool isValidSetup(Teuchos::FancyOStream &out) const
virtual void setSolver(Teuchos::RCP< Thyra::NonlinearSolverBase< Scalar > > solver=Teuchos::null)
Set solver.
void setParameterList(const Teuchos::RCP< Teuchos::ParameterList > &pl)
virtual void describe(Teuchos::FancyOStream &out, const Teuchos::EVerbosityLevel verbLevel) const