9#ifndef Tempus_StepperStaggeredForwardSensitivity_impl_hpp
10#define Tempus_StepperStaggeredForwardSensitivity_impl_hpp
12#include "Thyra_VectorStdOps.hpp"
13#include "Thyra_MultiVectorStdOps.hpp"
14#include "Thyra_DefaultMultiVectorProductVectorSpace.hpp"
15#include "Thyra_DefaultMultiVectorProductVector.hpp"
17#include "Tempus_StepperFactory.hpp"
32 this->
setParams(Teuchos::null, Teuchos::null);
42 const Teuchos::RCP<Teuchos::ParameterList>& pList,
43 const Teuchos::RCP<Teuchos::ParameterList>& sens_pList)
50 this->
setModel(appModel, sens_residual_model, sens_solve_model);
64 using Teuchos::ParameterList;
67 Teuchos::RCP<Teuchos::ParameterList> spl = Teuchos::parameterList();
69 spl->remove(
"Reuse State Linear Solver");
70 spl->remove(
"Force W Update");
72 appModel, sens_residual_model, sens_solve_model,
false, spl);
79 appModel, sens_residual_model, sens_solve_model, spl);
83 if (stateStepper_ == Teuchos::null)
84 stateStepper_ = sf->createStepper(stepperPL_, appModel);
86 stateStepper_->setModel(appModel);
88 if (sensitivityStepper_ == Teuchos::null)
89 sensitivityStepper_ = sf->createStepper(stepperPL_, fsa_model_);
91 sensitivityStepper_->setModel(fsa_model_);
93 this->isInitialized_ =
false;
98Teuchos::RCP<const Thyra::ModelEvaluator<Scalar> >
102 return combined_fsa_model_;
106template<
class Scalar>
109 Teuchos::RCP<Thyra::NonlinearSolverBase<Scalar> > solver)
111 stateStepper_->setSolver(solver);
112 sensitivityStepper_->setSolver(solver);
114 this->isInitialized_ =
false;
118template<
class Scalar>
123 this->checkInitialized();
127 using Teuchos::rcp_dynamic_cast;
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 =
143 solutionHistory->getNonconstParameterList();
146 RCP<SolutionState<Scalar> > state = solutionHistory->getCurrentState();
147 RCP<DMVPV> X, XDot, XDotDot;
148 X = rcp_dynamic_cast<DMVPV>(state->getX(),
true);
150 XDot = rcp_dynamic_cast<DMVPV>(state->getXDot(),
true);
151 if (state->getXDotDot() != Teuchos::null)
152 XDotDot = rcp_dynamic_cast<DMVPV>(state->getXDotDot(),
true);
156 RCP<VectorBase<Scalar> > x, xdot, xdotdot;
157 x = X->getNonconstMultiVector()->col(0);
158 xdot = XDot->getNonconstMultiVector()->col(0);
159 if (XDotDot != Teuchos::null)
160 xdotdot = XDotDot->getNonconstMultiVector()->col(0);
163 RCP<SolutionState<Scalar> > state_state = state->clone();
164 state_state->setX(x);
165 state_state->setXDot(xdot);
166 state_state->setXDotDot(xdotdot);
167 stateSolutionHistory_ = createSolutionHistoryPL<Scalar>(shPL);
168 stateSolutionHistory_->addState(state_state);
170 const int num_param = X->getMultiVector()->domain()->dim()-1;
171 TEUCHOS_ASSERT(num_param > 0);
172 const Teuchos::Range1D rng(1,num_param);
175 RCP<MultiVectorBase<Scalar> > dxdp, dxdotdp, dxdotdotdp;
176 dxdp = X->getNonconstMultiVector()->subView(rng);
177 dxdotdp = XDot->getNonconstMultiVector()->subView(rng);
178 if (XDotDot != Teuchos::null)
179 dxdotdotdp = XDotDot->getNonconstMultiVector()->subView(rng);
182 RCP<VectorBase<Scalar> > dxdp_vec, dxdotdp_vec, dxdotdotdp_vec;
183 RCP<const DMVPVS> prod_space =
184 multiVectorProductVectorSpace(X->getMultiVector()->range(), num_param);
185 dxdp_vec = multiVectorProductVector(prod_space, dxdp);
186 dxdotdp_vec = multiVectorProductVector(prod_space, dxdotdp);
187 if (XDotDot != Teuchos::null)
188 dxdotdotdp_vec = multiVectorProductVector(prod_space, dxdotdotdp);
191 RCP<SolutionState<Scalar> > sens_state = state->clone();
192 sens_state->setX(dxdp_vec);
193 sens_state->setXDot(dxdotdp_vec);
194 sens_state->setXDotDot(dxdotdotdp_vec);
195 sensSolutionHistory_ = createSolutionHistoryPL<Scalar>(shPL);
196 sensSolutionHistory_->addState(sens_state);
200 RCP<SolutionState<Scalar> > prod_state = solutionHistory->getWorkingState();
201 RCP<DMVPV> X, XDot, XDotDot;
202 X = rcp_dynamic_cast<DMVPV>(prod_state->getX(),
true);
203 XDot = rcp_dynamic_cast<DMVPV>(prod_state->getXDot(),
true);
204 if (prod_state->getXDotDot() != Teuchos::null)
205 XDotDot = rcp_dynamic_cast<DMVPV>(prod_state->getXDotDot(),
true);
209 stateSolutionHistory_->initWorkingState();
210 RCP<SolutionState<Scalar> > state = stateSolutionHistory_->getWorkingState();
211 state->getMetaData()->copy(prod_state->getMetaData());
212 stateStepper_->takeStep(stateSolutionHistory_);
215 assign(X->getNonconstMultiVector()->col(0).ptr(), *(state->getX()));
216 assign(XDot->getNonconstMultiVector()->col(0).ptr(), *(state->getXDot()));
217 if (XDotDot != Teuchos::null)
218 assign(XDotDot->getNonconstMultiVector()->col(0).ptr(),
219 *(state->getXDotDot()));
220 prod_state->setOrder(state->getOrder());
227 stateSolutionHistory_->promoteWorkingState();
230 fsa_model_->setForwardSolutionHistory(stateSolutionHistory_);
231 if (reuse_solver_ && stateStepper_->getSolver() != Teuchos::null)
232 fsa_model_->setSolver(stateStepper_->getSolver(), force_W_update_);
236 sensSolutionHistory_->initWorkingState();
237 RCP<SolutionState<Scalar> > sens_state =
238 sensSolutionHistory_->getWorkingState();
239 sens_state->getMetaData()->copy(prod_state->getMetaData());
240 sensitivityStepper_->takeStep(sensSolutionHistory_);
243 RCP<const MultiVectorBase<Scalar> > dxdp =
244 rcp_dynamic_cast<const DMVPV>(sens_state->getX(),
true)->getMultiVector();
245 const int num_param = dxdp->domain()->dim();
246 const Teuchos::Range1D rng(1,num_param);
247 assign(X->getNonconstMultiVector()->subView(rng).ptr(), *dxdp);
248 RCP<const MultiVectorBase<Scalar> > dxdotdp =
249 rcp_dynamic_cast<const DMVPV>(sens_state->getXDot(),
true)->getMultiVector();
250 assign(XDot->getNonconstMultiVector()->subView(rng).ptr(), *dxdotdp);
251 if (sens_state->getXDotDot() != Teuchos::null) {
252 RCP<const MultiVectorBase<Scalar> > dxdotdotdp =
253 rcp_dynamic_cast<const DMVPV>(sens_state->getXDotDot(),
true)->getMultiVector();
254 assign(XDotDot->getNonconstMultiVector()->subView(rng).ptr(), *dxdotdotdp);
256 prod_state->setOrder(std::min(state->getOrder(), sens_state->getOrder()));
263 sensSolutionHistory_->promoteWorkingState();
269template<
class Scalar>
270Teuchos::RCP<Tempus::StepperState<Scalar> >
276 Teuchos::RCP<Tempus::StepperState<Scalar> > stepperState =
282template<
class Scalar>
285 Teuchos::FancyOStream &out,
286 const Teuchos::EVerbosityLevel verbLevel)
const
288 out.setOutputToRootOnly(0);
292 out <<
"--- StepperStaggeredForwardSensitivity ---\n";
293 out <<
" stateStepper_ = " << stateStepper_ <<std::endl;
294 out <<
" sensitivityStepper_ = " << sensitivityStepper_ <<std::endl;
295 out <<
" combined_fsa_model_ = " << combined_fsa_model_ <<std::endl;
296 out <<
" fsa_model_ = " << fsa_model_ <<std::endl;
297 out <<
" stateSolutionHistory_ = " << stateSolutionHistory_ <<std::endl;
298 out <<
" sensSolutionHistory_ = " << sensSolutionHistory_ <<std::endl;
299 out <<
"------------------------------------------" << std::endl;
301 out << description() <<
"::describe:" << std::endl;
302 stateStepper_->describe(out, verbLevel);
304 sensitivityStepper_->describe(out, verbLevel);
308template<
class Scalar>
311 out.setOutputToRootOnly(0);
312 bool isValidSetup =
true;
316 if (stateStepper_ == Teuchos::null) {
317 isValidSetup =
false;
318 out <<
"The state stepper is not set!\n";
321 if (sensitivityStepper_ == Teuchos::null) {
322 isValidSetup =
false;
323 out <<
"The sensitivity stepper is not set!\n";
326 if (combined_fsa_model_ == Teuchos::null) {
327 isValidSetup =
false;
328 out <<
"The combined FSA model is not set!\n";
331 if (fsa_model_ == Teuchos::null) {
332 isValidSetup =
false;
333 out <<
"The FSA model is not set!\n";
340template <
class Scalar>
343 Teuchos::RCP<Teuchos::ParameterList>
const& pList)
345 if (pList == Teuchos::null) {
347 if (this->stepperPL_ == Teuchos::null) this->stepperPL_ =
348 Teuchos::rcp_const_cast<Teuchos::ParameterList>(this->getValidParameters());
350 this->stepperPL_ = pList;
357template<
class Scalar>
358Teuchos::RCP<const Teuchos::ParameterList>
362 return stateStepper_->getValidParameters();
366template <
class Scalar>
367Teuchos::RCP<Teuchos::ParameterList>
375template <
class Scalar>
376Teuchos::RCP<Teuchos::ParameterList>
380 Teuchos::RCP<Teuchos::ParameterList> temp_plist = stepperPL_;
381 stepperPL_ = Teuchos::null;
386template <
class Scalar>
389 Teuchos::RCP<Teuchos::ParameterList>
const& pList,
390 Teuchos::RCP<Teuchos::ParameterList>
const& spList)
392 if (pList == Teuchos::null)
393 stepperPL_ = Teuchos::rcp_const_cast<Teuchos::ParameterList>(this->getValidParameters());
397 if (spList == Teuchos::null)
398 sensPL_ = Teuchos::parameterList();
402 reuse_solver_ = sensPL_->get(
"Reuse State Linear Solver",
false);
403 force_W_update_ = sensPL_->get(
"Force W Update",
true);
410template <
class Scalar>
411Teuchos::RCP<const Thyra::VectorSpaceBase<Scalar> >
416 using Teuchos::rcp_dynamic_cast;
417 typedef Thyra::DefaultMultiVectorProductVectorSpace<Scalar> DMVPVS;
419 RCP<const Thyra::VectorSpaceBase<Scalar> > x_space =
420 stateStepper_->getModel()->get_x_space();
421 RCP<const DMVPVS> dxdp_space =
422 rcp_dynamic_cast<const DMVPVS>(sensitivityStepper_->getModel()->get_x_space(),
true);
423 const int num_param = dxdp_space->numBlocks();
424 RCP<const Thyra::VectorSpaceBase<Scalar> > prod_space =
425 multiVectorProductVectorSpace(x_space, num_param+1);
SolutionHistory is basically a container of SolutionStates. SolutionHistory maintains a collection of...
Teuchos::RCP< Teuchos::ParameterList > getNonconstParameterList()
virtual bool isValidSetup(Teuchos::FancyOStream &out) const
Teuchos::RCP< const Teuchos::ParameterList > getValidParameters() const
virtual Teuchos::RCP< const Thyra::ModelEvaluator< Scalar > > getModel() const
virtual Teuchos::RCP< Tempus::StepperState< Scalar > > getDefaultStepperState()
Get a default (initial) StepperState.
virtual void describe(Teuchos::FancyOStream &out, const Teuchos::EVerbosityLevel verbLevel) const
Teuchos::RCP< Teuchos::ParameterList > unsetParameterList()
void setParams(const Teuchos::RCP< Teuchos::ParameterList > &pl, const Teuchos::RCP< Teuchos::ParameterList > &spl)
void setParameterList(const Teuchos::RCP< Teuchos::ParameterList > &pl)
virtual void takeStep(const Teuchos::RCP< SolutionHistory< Scalar > > &solutionHistory)
Take the specified timestep, dt, and return true if successful.
StepperStaggeredForwardSensitivity()
Default constructor.
virtual void setModel(const Teuchos::RCP< const Thyra::ModelEvaluator< Scalar > > &appModel)
Teuchos::RCP< const Thyra::VectorSpaceBase< Scalar > > get_x_space() const
virtual void setSolver(Teuchos::RCP< Thyra::NonlinearSolverBase< Scalar > > solver=Teuchos::null)
Set solver.
StepperState is a simple class to hold state information about the stepper.
Thyra Base interface for time steppers.
void setStepperName(std::string s)
Set the stepper name.
virtual void initialize()
Initialize after construction and changing input parameters.
void setStepperType(std::string s)
Set the stepper type.
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)
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)