10 #ifndef Tempus_StepperStaggeredForwardSensitivity_impl_hpp
11 #define Tempus_StepperStaggeredForwardSensitivity_impl_hpp
13 #include "Thyra_VectorStdOps.hpp"
14 #include "Thyra_MultiVectorStdOps.hpp"
15 #include "Thyra_DefaultMultiVectorProductVectorSpace.hpp"
16 #include "Thyra_DefaultMultiVectorProductVector.hpp"
18 #include "Tempus_StepperFactory.hpp"
24 template <
class Scalar>
30 this->
setParams(Teuchos::null, Teuchos::null);
33 template <
class Scalar>
47 this->
setModel(appModel, sens_residual_model, sens_solve_model);
51 template <
class Scalar>
65 spl->
remove(
"Reuse State Linear Solver");
66 spl->
remove(
"Force W Update");
68 sens_solve_model,
false, spl);
75 appModel, sens_residual_model, sens_solve_model, spl);
79 if (stateStepper_ == Teuchos::null)
80 stateStepper_ = sf->createStepper(stepperPL_, appModel);
82 stateStepper_->setModel(appModel);
84 if (sensitivityStepper_ == Teuchos::null)
85 sensitivityStepper_ = sf->createStepper(stepperPL_, fsa_model_);
87 sensitivityStepper_->setModel(fsa_model_);
89 this->isInitialized_ =
false;
92 template <
class Scalar>
96 return combined_fsa_model_;
99 template <
class Scalar>
103 stateStepper_->setSolver(solver);
104 sensitivityStepper_->setSolver(solver);
106 this->isInitialized_ =
false;
109 template <
class Scalar>
113 this->checkInitialized();
117 using Teuchos::rcp_dynamic_cast;
119 using Thyra::createMember;
121 using Thyra::multiVectorProductVector;
122 using Thyra::multiVectorProductVectorSpace;
131 if (stateSolutionHistory_ == Teuchos::null) {
132 RCP<Teuchos::ParameterList> shPL =
133 solutionHistory->getNonconstParameterList();
136 RCP<SolutionState<Scalar> > state = solutionHistory->getCurrentState();
137 RCP<DMVPV> X, XDot, XDotDot;
138 X = rcp_dynamic_cast<DMVPV>(state->getX(),
true);
140 XDot = rcp_dynamic_cast<DMVPV>(state->getXDot(),
true);
141 if (state->getXDotDot() != Teuchos::null)
142 XDotDot = rcp_dynamic_cast<DMVPV>(state->getXDotDot(),
true);
146 RCP<VectorBase<Scalar> > x, xdot, xdotdot;
147 x = X->getNonconstMultiVector()->col(0);
148 xdot = XDot->getNonconstMultiVector()->col(0);
149 if (XDotDot != Teuchos::null)
150 xdotdot = XDotDot->getNonconstMultiVector()->col(0);
153 RCP<SolutionState<Scalar> > state_state = state->
clone();
154 state_state->setX(x);
155 state_state->setXDot(xdot);
156 state_state->setXDotDot(xdotdot);
157 stateSolutionHistory_ = createSolutionHistoryPL<Scalar>(shPL);
158 stateSolutionHistory_->addState(state_state);
160 const int num_param = X->getMultiVector()->domain()->dim() - 1;
165 RCP<MultiVectorBase<Scalar> > dxdp, dxdotdp, dxdotdotdp;
166 dxdp = X->getNonconstMultiVector()->subView(rng);
167 dxdotdp = XDot->getNonconstMultiVector()->subView(rng);
168 if (XDotDot != Teuchos::null)
169 dxdotdotdp = XDotDot->getNonconstMultiVector()->subView(rng);
172 RCP<VectorBase<Scalar> > dxdp_vec, dxdotdp_vec, dxdotdotdp_vec;
173 RCP<const DMVPVS> prod_space =
174 multiVectorProductVectorSpace(X->getMultiVector()->range(), num_param);
175 dxdp_vec = multiVectorProductVector(prod_space, dxdp);
176 dxdotdp_vec = multiVectorProductVector(prod_space, dxdotdp);
177 if (XDotDot != Teuchos::null)
178 dxdotdotdp_vec = multiVectorProductVector(prod_space, dxdotdotdp);
181 RCP<SolutionState<Scalar> > sens_state = state->clone();
182 sens_state->setX(dxdp_vec);
183 sens_state->setXDot(dxdotdp_vec);
184 sens_state->setXDotDot(dxdotdotdp_vec);
185 sensSolutionHistory_ = createSolutionHistoryPL<Scalar>(shPL);
186 sensSolutionHistory_->addState(sens_state);
190 RCP<SolutionState<Scalar> > prod_state = solutionHistory->getWorkingState();
191 RCP<DMVPV> X, XDot, XDotDot;
192 X = rcp_dynamic_cast<DMVPV>(prod_state->getX(),
true);
193 XDot = rcp_dynamic_cast<DMVPV>(prod_state->getXDot(),
true);
194 if (prod_state->getXDotDot() != Teuchos::null)
195 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_);
226 sensSolutionHistory_->initWorkingState();
227 RCP<SolutionState<Scalar> > sens_state =
228 sensSolutionHistory_->getWorkingState();
229 sens_state->getMetaData()->copy(prod_state->getMetaData());
230 sensitivityStepper_->takeStep(sensSolutionHistory_);
233 RCP<const MultiVectorBase<Scalar> > dxdp =
234 rcp_dynamic_cast<
const DMVPV>(sens_state->getX(),
true)->getMultiVector();
235 const int num_param = dxdp->domain()->dim();
237 assign(X->getNonconstMultiVector()->subView(rng).ptr(), *dxdp);
238 RCP<const MultiVectorBase<Scalar> > dxdotdp =
239 rcp_dynamic_cast<
const DMVPV>(sens_state->getXDot(),
true)
241 assign(XDot->getNonconstMultiVector()->subView(rng).ptr(), *dxdotdp);
242 if (sens_state->getXDotDot() != Teuchos::null) {
243 RCP<const MultiVectorBase<Scalar> > dxdotdotdp =
244 rcp_dynamic_cast<
const DMVPV>(sens_state->getXDotDot(),
true)
246 assign(XDotDot->getNonconstMultiVector()->subView(rng).ptr(), *dxdotdotdp);
248 prod_state->setOrder(std::min(state->getOrder(), sens_state->getOrder()));
255 sensSolutionHistory_->promoteWorkingState();
260 template <
class Scalar>
271 template <
class Scalar>
279 out <<
"--- StepperStaggeredForwardSensitivity ---\n";
280 out <<
" stateStepper_ = " << stateStepper_ << std::endl;
281 out <<
" sensitivityStepper_ = " << sensitivityStepper_ << std::endl;
282 out <<
" combined_fsa_model_ = " << combined_fsa_model_ << std::endl;
283 out <<
" fsa_model_ = " << fsa_model_ << std::endl;
284 out <<
" stateSolutionHistory_ = " << stateSolutionHistory_ << std::endl;
285 out <<
" sensSolutionHistory_ = " << sensSolutionHistory_ << std::endl;
286 out <<
"------------------------------------------" << std::endl;
288 out << description() <<
"::describe:" << std::endl;
289 stateStepper_->describe(out, verbLevel);
291 sensitivityStepper_->describe(out, verbLevel);
294 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";
326 template <
class Scalar>
330 if (pList == Teuchos::null) {
332 if (this->stepperPL_ == Teuchos::null)
334 this->getValidParameters());
337 this->stepperPL_ = pList;
343 template <
class Scalar>
347 return stateStepper_->getValidParameters();
350 template <
class Scalar>
357 template <
class Scalar>
362 stepperPL_ = Teuchos::null;
366 template <
class Scalar>
371 if (pList == Teuchos::null)
373 this->getValidParameters());
377 if (spList == Teuchos::null)
378 sensPL_ = Teuchos::parameterList();
382 reuse_solver_ = sensPL_->
get(
"Reuse State Linear Solver",
false);
383 force_W_update_ = sensPL_->get(
"Force W Update",
true);
389 template <
class Scalar>
394 using Teuchos::rcp_dynamic_cast;
397 RCP<const Thyra::VectorSpaceBase<Scalar> > x_space =
398 stateStepper_->getModel()->get_x_space();
399 RCP<const DMVPVS> dxdp_space = rcp_dynamic_cast<
const DMVPVS>(
400 sensitivityStepper_->getModel()->get_x_space(),
true);
401 const int num_param = dxdp_space->
numBlocks();
402 RCP<const Thyra::VectorSpaceBase<Scalar> > prod_space =
403 multiVectorProductVectorSpace(x_space, num_param + 1);
409 #endif // Tempus_StepperStaggeredForwardSensitivity_impl_hpp
Teuchos::RCP< Teuchos::ParameterList > getNonconstParameterList()
RCP< const VectorSpaceBase< Scalar > > clone() const
T & get(const std::string &name, T def_value)
Teuchos::RCP< SensitivityModelEvaluatorBase< Scalar > > wrapStaggeredFSAModelEvaluator(const Teuchos::RCP< const Thyra::ModelEvaluator< Scalar > > &model, const Teuchos::RCP< const Thyra::ModelEvaluator< Scalar > > &sens_residual_model, const Teuchos::RCP< const Thyra::ModelEvaluator< Scalar > > &sens_solve_model, const bool is_pseudotransient, const Teuchos::RCP< const Teuchos::ParameterList > &pList=Teuchos::null)
void setStepperName(std::string s)
Set the stepper name.
virtual void takeStep(const Teuchos::RCP< SolutionHistory< Scalar > > &solutionHistory)
Take the specified timestep, dt, and return true if successful.
StepperStaggeredForwardSensitivity()
Default constructor.
Teuchos::RCP< const Thyra::VectorSpaceBase< Scalar > > get_x_space() const
virtual void initialize()
Initialize after construction and changing input parameters.
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)
bool remove(std::string const &name, bool throwIfNotExists=true)
StepperState is a simple class to hold state information about the stepper.
Teuchos::RCP< Teuchos::ParameterList > unsetParameterList()
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
virtual void describe(Teuchos::FancyOStream &out, const Teuchos::EVerbosityLevel verbLevel) const
Teuchos::RCP< SensitivityModelEvaluatorBase< Scalar > > wrapCombinedFSAModelEvaluator(const Teuchos::RCP< const Thyra::ModelEvaluator< Scalar > > &model, const Teuchos::RCP< const Thyra::ModelEvaluator< Scalar > > &sens_residual_model, const Teuchos::RCP< const Thyra::ModelEvaluator< Scalar > > &sens_solve_model, const Teuchos::RCP< const Teuchos::ParameterList > &pList=Teuchos::null)
SolutionHistory is basically a container of SolutionStates. SolutionHistory maintains a collection of...
basic_FancyOStream & setOutputToRootOnly(const int rootRank)
Teuchos::RCP< const Teuchos::ParameterList > getValidParameters() const
#define TEUCHOS_ASSERT(assertion_test)
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 setStepperType(std::string s)
Set the stepper type.
virtual Teuchos::RCP< const Thyra::ModelEvaluator< Scalar > > getModel() const
void setParameterList(const Teuchos::RCP< Teuchos::ParameterList > &pl)
virtual void describe(Teuchos::FancyOStream &out, const Teuchos::EVerbosityLevel verbLevel) const