Tempus Version of the Day
Time Integration
Loading...
Searching...
No Matches
Tempus_StepperStaggeredForwardSensitivity_impl.hpp
Go to the documentation of this file.
1// @HEADER
2// ****************************************************************************
3// Tempus: Copyright (2017) Sandia Corporation
4//
5// Distributed under BSD 3-clause license (See accompanying file Copyright.txt)
6// ****************************************************************************
7// @HEADER
8
9#ifndef Tempus_StepperStaggeredForwardSensitivity_impl_hpp
10#define Tempus_StepperStaggeredForwardSensitivity_impl_hpp
11
12#include "Thyra_VectorStdOps.hpp"
13#include "Thyra_MultiVectorStdOps.hpp"
14#include "Thyra_DefaultMultiVectorProductVectorSpace.hpp"
15#include "Thyra_DefaultMultiVectorProductVector.hpp"
16
17#include "Tempus_StepperFactory.hpp"
20
21
22namespace Tempus {
23
24
25template<class Scalar>
28 : stepMode_(SensitivityStepMode::Forward)
29{
30 this->setStepperName( "StaggeredForwardSensitivity");
31 this->setStepperType( "StaggeredForwardSensitivity");
32 this->setParams(Teuchos::null, Teuchos::null);
33}
34
35
36template<class Scalar>
39 const Teuchos::RCP<const Thyra::ModelEvaluator<Scalar> >& appModel,
40 const Teuchos::RCP<const Thyra::ModelEvaluator<Scalar> >& sens_residual_model,
41 const Teuchos::RCP<const Thyra::ModelEvaluator<Scalar> >& sens_solve_model,
42 const Teuchos::RCP<Teuchos::ParameterList>& pList,
43 const Teuchos::RCP<Teuchos::ParameterList>& sens_pList)
44 : stepMode_(SensitivityStepMode::Forward)
45{
46 // Set all the input parameters and call initialize
47 this->setStepperName( "StaggeredForwardSensitivity");
48 this->setStepperType( "StaggeredForwardSensitivity");
49 this->setParams(pList, sens_pList);
50 this->setModel(appModel, sens_residual_model, sens_solve_model);
51 this->initialize();
52}
53
54
55template<class Scalar>
58 const Teuchos::RCP<const Thyra::ModelEvaluator<Scalar> >& appModel,
59 const Teuchos::RCP<const Thyra::ModelEvaluator<Scalar> >& sens_residual_model,
60 const Teuchos::RCP<const Thyra::ModelEvaluator<Scalar> >& sens_solve_model)
61{
62 using Teuchos::RCP;
63 using Teuchos::rcp;
64 using Teuchos::ParameterList;
65
66 // Create forward sensitivity model evaluator wrapper
67 Teuchos::RCP<Teuchos::ParameterList> spl = Teuchos::parameterList();
68 *spl = *sensPL_;
69 spl->remove("Reuse State Linear Solver");
70 spl->remove("Force W Update");
72 appModel, sens_residual_model, sens_solve_model, false, spl);
73
74 // Create combined FSA ME which serves as "the" ME for this stepper,
75 // so that getModel() has a ME consistent the FSA problem (including both
76 // state and sensitivity components), e.g., the integrator may call
77 // getModel()->getNominalValues(), which needs to be consistent.
78 combined_fsa_model_ = wrapCombinedFSAModelEvaluator(
79 appModel, sens_residual_model, sens_solve_model, spl);
80
81 // Create state and sensitivity steppers
82 RCP<StepperFactory<Scalar> > sf =Teuchos::rcp(new StepperFactory<Scalar>());
83 if (stateStepper_ == Teuchos::null)
84 stateStepper_ = sf->createStepper(stepperPL_, appModel);
85 else
86 stateStepper_->setModel(appModel);
87
88 if (sensitivityStepper_ == Teuchos::null)
89 sensitivityStepper_ = sf->createStepper(stepperPL_, fsa_model_);
90 else
91 sensitivityStepper_->setModel(fsa_model_);
92
93 this->isInitialized_ = false;
94}
95
96
97template<class Scalar>
98Teuchos::RCP<const Thyra::ModelEvaluator<Scalar> >
100getModel() const
101{
102 return combined_fsa_model_;
103}
104
105
106template<class Scalar>
109 Teuchos::RCP<Thyra::NonlinearSolverBase<Scalar> > solver)
110{
111 stateStepper_->setSolver(solver);
112 sensitivityStepper_->setSolver(solver);
113
114 this->isInitialized_ = false;
115}
116
117
118template<class Scalar>
121 const Teuchos::RCP<SolutionHistory<Scalar> >& solutionHistory)
122{
123 this->checkInitialized();
124
125 using Teuchos::RCP;
126 using Teuchos::rcp;
127 using Teuchos::rcp_dynamic_cast;
128 using Thyra::VectorBase;
130 using Thyra::assign;
131 using Thyra::createMember;
132 using Thyra::multiVectorProductVector;
133 using Thyra::multiVectorProductVectorSpace;
134 typedef Thyra::DefaultMultiVectorProductVector<Scalar> DMVPV;
135 typedef Thyra::DefaultMultiVectorProductVectorSpace<Scalar> DMVPVS;
136
137 // Initialize state, sensitivity solution histories if necessary.
138 // We only need to split the solution history into state and sensitivity
139 // components for the first step, otherwise the state and sensitivity
140 // histories are updated from the previous step.
141 if (stateSolutionHistory_ == Teuchos::null) {
142 RCP<Teuchos::ParameterList> shPL =
143 solutionHistory->getNonconstParameterList();
144
145 // Get product X, XDot, XDotDot
146 RCP<SolutionState<Scalar> > state = solutionHistory->getCurrentState();
147 RCP<DMVPV> X, XDot, XDotDot;
148 X = rcp_dynamic_cast<DMVPV>(state->getX(),true);
149
150 XDot = rcp_dynamic_cast<DMVPV>(state->getXDot(),true);
151 if (state->getXDotDot() != Teuchos::null)
152 XDotDot = rcp_dynamic_cast<DMVPV>(state->getXDotDot(),true);
153
154 // Pull out state components (has to be non-const because of SolutionState
155 // constructor)
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);
161
162 // Create state solution history
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);
169
170 const int num_param = X->getMultiVector()->domain()->dim()-1;
171 TEUCHOS_ASSERT(num_param > 0);
172 const Teuchos::Range1D rng(1,num_param);
173
174 // Pull out sensitivity components
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);
180
181 // Create sensitivity product vectors
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);
189
190 // Create sensitivity solution history
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);
197 }
198
199 // Get our working 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);
206
207 // Take step for state equations
209 stateSolutionHistory_->initWorkingState();
210 RCP<SolutionState<Scalar> > state = stateSolutionHistory_->getWorkingState();
211 state->getMetaData()->copy(prod_state->getMetaData());
212 stateStepper_->takeStep(stateSolutionHistory_);
213
214 // Set state components of product state
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());
221
222 // If step passed promote the state, otherwise fail and stop
223 if (state->getSolutionStatus() == Status::FAILED) {
224 prod_state->setSolutionStatus(Status::FAILED);
225 return;
226 }
227 stateSolutionHistory_->promoteWorkingState();
228
229 // Get forward state in sensitivity model evaluator
230 fsa_model_->setForwardSolutionHistory(stateSolutionHistory_);
231 if (reuse_solver_ && stateStepper_->getSolver() != Teuchos::null)
232 fsa_model_->setSolver(stateStepper_->getSolver(), force_W_update_);
233
234 // Take step in sensitivity equations
236 sensSolutionHistory_->initWorkingState();
237 RCP<SolutionState<Scalar> > sens_state =
238 sensSolutionHistory_->getWorkingState();
239 sens_state->getMetaData()->copy(prod_state->getMetaData());
240 sensitivityStepper_->takeStep(sensSolutionHistory_);
241
242 // Set sensitivity components of product state
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);
255 }
256 prod_state->setOrder(std::min(state->getOrder(), sens_state->getOrder()));
257
258 // If step passed promote the state, otherwise fail and stop
259 if (sens_state->getSolutionStatus() == Status::FAILED) {
260 prod_state->setSolutionStatus(Status::FAILED);
261 }
262 else {
263 sensSolutionHistory_->promoteWorkingState();
264 prod_state->setSolutionStatus(Status::PASSED);
265 }
266}
267
268
269template<class Scalar>
270Teuchos::RCP<Tempus::StepperState<Scalar> >
273{
274 // ETP: Note, maybe this should be a special stepper state that combines
275 // states from both state and sensitivity steppers?
276 Teuchos::RCP<Tempus::StepperState<Scalar> > stepperState =
277 rcp(new StepperState<Scalar>(description()));
278 return stepperState;
279}
280
281
282template<class Scalar>
285 Teuchos::FancyOStream &out,
286 const Teuchos::EVerbosityLevel verbLevel) const
287{
288 out.setOutputToRootOnly(0);
289 out << std::endl;
290 Stepper<Scalar>::describe(out, verbLevel);
291
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;
300
301 out << description() << "::describe:" << std::endl;
302 stateStepper_->describe(out, verbLevel);
303 out << std::endl;
304 sensitivityStepper_->describe(out, verbLevel);
305}
306
307
308template<class Scalar>
309bool StepperStaggeredForwardSensitivity<Scalar>::isValidSetup(Teuchos::FancyOStream & out) const
310{
311 out.setOutputToRootOnly(0);
312 bool isValidSetup = true;
313
314 if ( !Stepper<Scalar>::isValidSetup(out) ) isValidSetup = false;
315
316 if (stateStepper_ == Teuchos::null) {
317 isValidSetup = false;
318 out << "The state stepper is not set!\n";
319 }
320
321 if (sensitivityStepper_ == Teuchos::null) {
322 isValidSetup = false;
323 out << "The sensitivity stepper is not set!\n";
324 }
325
326 if (combined_fsa_model_ == Teuchos::null) {
327 isValidSetup = false;
328 out << "The combined FSA model is not set!\n";
329 }
330
331 if (fsa_model_ == Teuchos::null) {
332 isValidSetup = false;
333 out << "The FSA model is not set!\n";
334 }
335
336 return isValidSetup;
337}
338
339
340template <class Scalar>
343 Teuchos::RCP<Teuchos::ParameterList> const& pList)
344{
345 if (pList == Teuchos::null) {
346 // Create default parameters if null, otherwise keep current parameters.
347 if (this->stepperPL_ == Teuchos::null) this->stepperPL_ =
348 Teuchos::rcp_const_cast<Teuchos::ParameterList>(this->getValidParameters());
349 } else {
350 this->stepperPL_ = pList;
351 }
352 // Can not validate because of optional Parameters (e.g., Solver Name).
353 //stepperPL_->validateParametersAndSetDefaults(*this->getValidParameters());
354}
355
356
357template<class Scalar>
358Teuchos::RCP<const Teuchos::ParameterList>
360getValidParameters() const
361{
362 return stateStepper_->getValidParameters();
363}
364
365
366template <class Scalar>
367Teuchos::RCP<Teuchos::ParameterList>
373
374
375template <class Scalar>
376Teuchos::RCP<Teuchos::ParameterList>
379{
380 Teuchos::RCP<Teuchos::ParameterList> temp_plist = stepperPL_;
381 stepperPL_ = Teuchos::null;
382 return temp_plist;
383}
384
385
386template <class Scalar>
389 Teuchos::RCP<Teuchos::ParameterList> const& pList,
390 Teuchos::RCP<Teuchos::ParameterList> const& spList)
391{
392 if (pList == Teuchos::null)
393 stepperPL_ = Teuchos::rcp_const_cast<Teuchos::ParameterList>(this->getValidParameters());
394 else
395 stepperPL_ = pList;
396
397 if (spList == Teuchos::null)
398 sensPL_ = Teuchos::parameterList();
399 else
400 sensPL_ = spList;
401
402 reuse_solver_ = sensPL_->get("Reuse State Linear Solver", false);
403 force_W_update_ = sensPL_->get("Force W Update", true);
404
405 // Can not validate because of optional Parameters (e.g., Solver Name).
406 //stepperPL_->validateParametersAndSetDefaults(*this->getValidParameters());
407}
408
409
410template <class Scalar>
411Teuchos::RCP<const Thyra::VectorSpaceBase<Scalar> >
413get_x_space() const
414{
415 using Teuchos::RCP;
416 using Teuchos::rcp_dynamic_cast;
417 typedef Thyra::DefaultMultiVectorProductVectorSpace<Scalar> DMVPVS;
418
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);
426 return prod_space;
427}
428
429} // namespace Tempus
430
431#endif // Tempus_StepperStaggeredForwardSensitivity_impl_hpp
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
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.
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.
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)