Tempus Version of the Day
Time Integration
Loading...
Searching...
No Matches
Tempus_WrapperModelEvaluatorPairPartIMEX_CombinedFSA_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_ModelEvaluatorPairPartIMEX_CombinedFSA_impl_hpp
10#define Tempus_ModelEvaluatorPairPartIMEX_CombinedFSA_impl_hpp
11
12#include "Thyra_VectorStdOps.hpp"
13#include "Thyra_MultiVectorStdOps.hpp"
14
15namespace Tempus {
16
17template <typename Scalar>
20 const Teuchos::RCP<const WrapperModelEvaluatorPairPartIMEX_Basic<Scalar> >& forwardModel,
21 const Teuchos::RCP<const Teuchos::ParameterList>& pList) :
22 forwardModel_(forwardModel),
23 use_dfdp_as_tangent_(false),
24 y_tangent_index_(3)
25{
26 using Teuchos::RCP;
27 using Teuchos::rcp;
28 using Thyra::multiVectorProductVectorSpace;
29
30 RCP<Teuchos::ParameterList> pl = Teuchos::parameterList();
31 if (pList != Teuchos::null)
32 *pl = *pList;
33 pl->validateParametersAndSetDefaults(*this->getValidParameters());
34 use_dfdp_as_tangent_ = pl->get<bool>("Use DfDp as Tangent");
35 y_tangent_index_ = pl->get<int>("Sensitivity Y Tangent Index");
36 pl->remove("Sensitivity Y Tangent Index");
37
38 appExplicitModel_ = forwardModel_->getExplicitModel();
39 appImplicitModel_ = forwardModel_->getImplicitModel();
42
43 const int y_param_index = forwardModel_->getParameterIndex();
44 const int sens_param_index = pl->get<int>("Sensitivity Parameter Index");
45 const int num_sens_param =
46 appImplicitModel_->get_p_space(sens_param_index)->dim();
47 RCP<const Thyra::VectorSpaceBase<Scalar> > explicit_y_space =
48 appImplicitModel_->get_p_space(y_param_index);
49 RCP<const Thyra::VectorSpaceBase<Scalar> > implicit_x_space =
50 appImplicitModel_->get_x_space();
52 multiVectorProductVectorSpace(explicit_y_space, num_sens_param+1);
54 multiVectorProductVectorSpace(explicit_y_space, num_sens_param);
56 multiVectorProductVectorSpace(implicit_x_space, num_sens_param+1);
57
59 forwardModel_->getNumExplicitOnlyBlocks(),
60 y_param_index);
61}
62
63template <typename Scalar>
64void
67{
68 using Teuchos::RCP;
69 using Teuchos::rcp_dynamic_cast;
70
71 this->useImplicitModel_ = true;
72 this->wrapperImplicitInArgs_ = this->createInArgs();
73 this->wrapperImplicitOutArgs_ = this->createOutArgs();
74 this->useImplicitModel_ = false;
75
76 RCP<const Thyra::VectorBase<Scalar> > z =
77 this->explicitModel_->getNominalValues().get_x();
78
79 // A Thyra::VectorSpace requirement
80 TEUCHOS_TEST_FOR_EXCEPTION( !(getIMEXVector(z)->space()->isCompatible(
81 *(this->implicitModel_->get_x_space()))),
82 std::logic_error,
83 "Error - WrapperModelEvaluatorPairIMEX_CombinedFSA::initialize()\n"
84 " Explicit and Implicit vector x spaces are incompatible!\n"
85 " Explicit vector x space = " << *(getIMEXVector(z)->space()) << "\n"
86 " Implicit vector x space = " << *(this->implicitModel_->get_x_space()) <<
87 "\n");
88
89 // Valid number of blocks?
90 const RCP<const DMVPV> z_dmvpv = rcp_dynamic_cast<const DMVPV>(z,true);
91 const RCP<const Thyra::MultiVectorBase<Scalar> > z_mv =
92 z_dmvpv->getMultiVector();
93 RCP<const PMVB> zPVector = rcp_dynamic_cast<const PMVB>(z_mv);
94 TEUCHOS_TEST_FOR_EXCEPTION( zPVector == Teuchos::null, std::logic_error,
95 "Error - WrapperModelEvaluatorPairPartIMEX_CombinedFSA::initialize()\n"
96 " was given a VectorBase that could not be cast to a\n"
97 " ProductVectorBase!\n");
98
99 int numBlocks = zPVector->productSpace()->numBlocks();
100
101 TEUCHOS_TEST_FOR_EXCEPTION( !(0 <= this->numExplicitOnlyBlocks_ &&
102 this->numExplicitOnlyBlocks_ < numBlocks),
103 std::logic_error,
104 "Error - WrapperModelEvaluatorPairPartIMEX_CombinedFSA::initialize()\n"
105 "Invalid number of explicit-only blocks = " <<
106 this->numExplicitOnlyBlocks_ << "\n"
107 "Should be in the interval [0, numBlocks) = [0, " << numBlocks << ")\n");
108}
109
110template <typename Scalar>
111Teuchos::RCP<const Thyra::VectorSpaceBase<Scalar> >
113get_p_space(int i) const
114{
115 if (this->useImplicitModel_) {
116 if (i == this->parameterIndex_)
117 return explicit_y_dydp_prod_space_;
118 else
119 return appImplicitModel_->get_p_space(i);
120 }
121
122 return appExplicitModel_->get_p_space(i);
123}
124
125template <typename Scalar>
126Teuchos::RCP<Thyra::VectorBase<Scalar> >
128getIMEXVector(const Teuchos::RCP<Thyra::VectorBase<Scalar> > & full) const
129{
130 using Teuchos::RCP;
131 using Teuchos::rcp_dynamic_cast;
133 using Thyra::VectorBase;
134 using Thyra::multiVectorProductVector;
135
136 // CombinedFSA ME stores vectors as DMVPV's. To extract the implicit
137 // part of the vector, cast it to DMVPV, extract the multi-vector,
138 // cast it to a product multi-vector, extract the IMEX block, then
139 // create a DMVPV from it.
140
141 if(full == Teuchos::null)
142 return Teuchos::null;
143
144 if (this->numExplicitOnlyBlocks_==0)
145 return full;
146
147 const RCP<DMVPV> full_dmvpv = rcp_dynamic_cast<DMVPV>(full,true);
148 const RCP<MultiVectorBase<Scalar> > full_mv =
149 full_dmvpv->getNonconstMultiVector();
150 const RCP<PMVB> blk_full_mv = rcp_dynamic_cast<PMVB>(full_mv,true);
151
152 // special case where the implicit terms are not blocked
153 const int numBlocks = blk_full_mv->productSpace()->numBlocks();
154 const int numExplicitBlocks = this->numExplicitOnlyBlocks_;
155 if (numBlocks == numExplicitBlocks+1) {
156 const RCP<MultiVectorBase<Scalar> > imex_mv =
157 blk_full_mv->getNonconstMultiVectorBlock(numExplicitBlocks);
158 return multiVectorProductVector(imex_x_dxdp_prod_space_, imex_mv);
159 }
160
161 // Not supposed to get here, apparently
162 TEUCHOS_ASSERT(false);
163 return Teuchos::null;
164}
165
166template <typename Scalar>
167Teuchos::RCP<const Thyra::VectorBase<Scalar> >
169getIMEXVector(const Teuchos::RCP<const Thyra::VectorBase<Scalar> > & full) const
170{
171 using Teuchos::RCP;
172 using Teuchos::rcp_dynamic_cast;
174 using Thyra::VectorBase;
175 using Thyra::multiVectorProductVector;
176
177 // CombinedFSA ME stores vectors as DMVPV's. To extract the implicit
178 // part of the vector, cast it to DMVPV, extract the multi-vector,
179 // cast it to a product multi-vector, extract the IMEX block, then
180 // create a DMVPV from it.
181
182 if(full == Teuchos::null)
183 return Teuchos::null;
184
185 if (this->numExplicitOnlyBlocks_==0)
186 return full;
187
188 const RCP<const DMVPV> full_dmvpv = rcp_dynamic_cast<const DMVPV>(full,true);
189 const RCP<const MultiVectorBase<Scalar> > full_mv =
190 full_dmvpv->getMultiVector();
191 const RCP<const PMVB> blk_full_mv =
192 rcp_dynamic_cast<const PMVB>(full_mv,true);
193
194 // special case where the implicit terms are not blocked
195 const int numBlocks = blk_full_mv->productSpace()->numBlocks();
196 const int numExplicitBlocks = this->numExplicitOnlyBlocks_;
197 if (numBlocks == numExplicitBlocks+1) {
198 const RCP<const MultiVectorBase<Scalar> > imex_mv =
199 blk_full_mv->getMultiVectorBlock(numExplicitBlocks);
200 return multiVectorProductVector(imex_x_dxdp_prod_space_, imex_mv);
201 }
202
203 // Not supposed to get here, apparently
204 TEUCHOS_ASSERT(false);
205 return Teuchos::null;
206}
207
208template <typename Scalar>
209Teuchos::RCP<Thyra::VectorBase<Scalar> >
212 const Teuchos::RCP<Thyra::VectorBase<Scalar> > & full) const
213{
214 using Teuchos::RCP;
215 using Teuchos::rcp_dynamic_cast;
217 using Thyra::VectorBase;
218 using Thyra::multiVectorProductVectorSpace;
219 using Thyra::multiVectorProductVector;
220
221 // CombinedFSA ME stores vectors as DMVPV's. To extract the explicit
222 // part of the vector, cast it to DMVPV, extract the multi-vector,
223 // cast it to a product multi-vector, extract the explicit block, then
224 // create a DMVPV from it.
225
226 if(full == Teuchos::null)
227 return Teuchos::null;
228
229 if (this->numExplicitOnlyBlocks_==0)
230 return full;
231
232 const RCP<DMVPV> full_dmvpv = rcp_dynamic_cast<DMVPV>(full,true);
233 const RCP<MultiVectorBase<Scalar> > full_mv =
234 full_dmvpv->getNonconstMultiVector();
235 const RCP<PMVB> blk_full_mv = rcp_dynamic_cast<PMVB>(full_mv,true);
236
237 // special case where the explicit terms are not blocked
238 const int numExplicitBlocks = this->numExplicitOnlyBlocks_;
239 if (numExplicitBlocks == 1) {
240 const RCP<MultiVectorBase<Scalar> > explicit_mv =
241 blk_full_mv->getNonconstMultiVectorBlock(0);
242 return multiVectorProductVector(explicit_y_dydp_prod_space_, explicit_mv);
243 }
244
245 // Not supposed to get here, apparently
246 TEUCHOS_ASSERT(false);
247 return Teuchos::null;
248}
249
250template <typename Scalar>
251Teuchos::RCP<const Thyra::VectorBase<Scalar> >
254 const Teuchos::RCP<const Thyra::VectorBase<Scalar> > & full) const
255{
256 using Teuchos::RCP;
257 using Teuchos::rcp_dynamic_cast;
259 using Thyra::VectorBase;
260 using Thyra::multiVectorProductVectorSpace;
261 using Thyra::multiVectorProductVector;
262
263 // CombinedFSA ME stores vectors as DMVPV's. To extract the explicit
264 // part of the vector, cast it to DMVPV, extract the multi-vector,
265 // cast it to a product multi-vector, extract the explicit block, then
266 // create a DMVPV from it.
267
268 if(full == Teuchos::null)
269 return Teuchos::null;
270
271 if (this->numExplicitOnlyBlocks_==0)
272 return full;
273
274 const RCP<const DMVPV> full_dmvpv = rcp_dynamic_cast<const DMVPV>(full,true);
275 const RCP<const MultiVectorBase<Scalar> > full_mv =
276 full_dmvpv->getMultiVector();
277 const RCP<const PMVB> blk_full_mv =
278 rcp_dynamic_cast<const PMVB>(full_mv,true);
279
280 // special case where the explicit terms are not blocked
281 const int numExplicitBlocks = this->numExplicitOnlyBlocks_;
282 if (numExplicitBlocks == 1) {
283 const RCP<const MultiVectorBase<Scalar> > explicit_mv =
284 blk_full_mv->getMultiVectorBlock(0);
285 return multiVectorProductVector(explicit_y_dydp_prod_space_, explicit_mv);
286 }
287
288 // Not supposed to get here, apparently
289 TEUCHOS_ASSERT(false);
290 return Teuchos::null;
291}
292
293template <typename Scalar>
294Teuchos::RCP<const Thyra::ModelEvaluator<Scalar> >
296getForwardModel() const
297{
298 return forwardModel_;
299}
300
301template <typename Scalar>
302Thyra::ModelEvaluatorBase::InArgs<Scalar>
304createInArgs() const
305{
306 using Teuchos::RCP;
307 using Teuchos::rcp_dynamic_cast;
308 using Thyra::createMember;
309
310 Thyra::ModelEvaluatorBase::InArgs<Scalar> inArgs = Base::createInArgs();
311
312 // Set p to be the correct product vector form for the explicit only vector y
313 if (this->useImplicitModel_ == true) {
314 RCP< const Thyra::VectorBase<Scalar> > y =
315 inArgs.get_p(this->parameterIndex_);
316 if (y != Teuchos::null) {
317 RCP<DMVPV> y_dydp =
318 rcp_dynamic_cast<DMVPV>(createMember(*explicit_y_dydp_prod_space_));
319 Thyra::assign(y_dydp->getNonconstMultiVector().ptr(), Scalar(0.0));
320 Thyra::assign(y_dydp->getNonconstMultiVector()->col(0).ptr(), *y);
321 inArgs.set_p(this->parameterIndex_, y_dydp);
322 }
323 }
324 return inArgs;
325}
326
327template <typename Scalar>
328void
330evalModelImpl(const Thyra::ModelEvaluatorBase::InArgs<Scalar> & inArgs,
331 const Thyra::ModelEvaluatorBase::OutArgs<Scalar> & outArgs) const
332{
333 typedef Thyra::ModelEvaluatorBase MEB;
334 using Teuchos::RCP;
335 using Teuchos::rcp_dynamic_cast;
336 using Teuchos::Range1D;
337
338 const int p_index = this->parameterIndex_;
339
340 //
341 // From Base::evalModelImpl()
342 //
343 RCP<const Thyra::VectorBase<Scalar> > x = inArgs.get_x();
344 RCP<Thyra::VectorBase<Scalar> > x_dot =
345 Thyra::createMember(fsaImplicitModel_->get_x_space());
346 this->timeDer_->compute(x, x_dot);
347
348 MEB::InArgs<Scalar> fsaImplicitInArgs (this->wrapperImplicitInArgs_);
349 MEB::OutArgs<Scalar> fsaImplicitOutArgs(this->wrapperImplicitOutArgs_);
350 fsaImplicitInArgs.set_x(x);
351 fsaImplicitInArgs.set_x_dot(x_dot);
352 for (int i=0; i<fsaImplicitModel_->Np(); ++i) {
353 // Copy over parameters except for the parameter for explicit-only vector!
354 if ((inArgs.get_p(i) != Teuchos::null) && (i != p_index))
355 fsaImplicitInArgs.set_p(i, inArgs.get_p(i));
356 }
357
358 // p-vector for index parameterIndex_ is part of the IMEX solution vector,
359 // and therefore is an n+1 column multi-vector where n is the number of
360 // sensitivity parameters. Pull out the sensitivity components before
361 // passing along to the ME, then use them for adding in dg/dy*dy/dp term.
362 RCP<const Thyra::VectorBase<Scalar> > y;
363 RCP<const Thyra::MultiVectorBase<Scalar> > dydp;
364 if (fsaImplicitInArgs.get_p(p_index) != Teuchos::null) {
365 RCP<const Thyra::VectorBase<Scalar> > p =
366 fsaImplicitInArgs.get_p(p_index);
367 RCP<const Thyra::MultiVectorBase<Scalar> > p_mv =
368 rcp_dynamic_cast<const DMVPV>(p,true)->getMultiVector();
369 const int num_param = p_mv->domain()->dim()-1;
370 y = p_mv->col(0);
371 dydp = p_mv->subView(Range1D(1,num_param));
372 fsaImplicitInArgs.set_p(p_index, y);
373 }
374 if (use_dfdp_as_tangent_) {
375 RCP< const Thyra::VectorBase<Scalar> > dydp_vec =
376 Thyra::multiVectorProductVector(explicit_dydp_prod_space_, dydp);
377 fsaImplicitInArgs.set_p(y_tangent_index_, dydp_vec);
378 }
379
380 fsaImplicitOutArgs.set_f(outArgs.get_f());
381 fsaImplicitOutArgs.set_W_op(outArgs.get_W_op());
382
383 fsaImplicitModel_->evalModel(fsaImplicitInArgs,fsaImplicitOutArgs);
384
385 // Compute derivative of implicit residual with respect to explicit only
386 // vector y, which is passed as a parameter
387 if (!use_dfdp_as_tangent_ && outArgs.get_f() != Teuchos::null) {
388 MEB::InArgs<Scalar> appImplicitInArgs =
389 appImplicitModel_->getNominalValues();
390 RCP< const Thyra::VectorBase<Scalar> > app_x =
391 rcp_dynamic_cast<const DMVPV>(x,true)->getMultiVector()->col(0);
392 RCP< const Thyra::VectorBase<Scalar> > app_x_dot =
393 rcp_dynamic_cast<const DMVPV>(x_dot,true)->getMultiVector()->col(0);
394 appImplicitInArgs.set_x(app_x);
395 appImplicitInArgs.set_x_dot(app_x_dot);
396 for (int i=0; i<appImplicitModel_->Np(); ++i) {
397 if (i != p_index)
398 appImplicitInArgs.set_p(i, inArgs.get_p(i));
399 }
400 appImplicitInArgs.set_p(p_index, y);
401 if (appImplicitInArgs.supports(MEB::IN_ARG_t))
402 appImplicitInArgs.set_t(inArgs.get_t());
403 MEB::OutArgs<Scalar> appImplicitOutArgs =
404 appImplicitModel_->createOutArgs();
405 MEB::DerivativeSupport dfdp_support =
406 appImplicitOutArgs.supports(MEB::OUT_ARG_DfDp, p_index);
407 Thyra::EOpTransp trans = Thyra::NOTRANS;
408 if (dfdp_support.supports(MEB::DERIV_LINEAR_OP)) {
409 if (my_dfdp_op_ == Teuchos::null)
410 my_dfdp_op_ = appImplicitModel_->create_DfDp_op(p_index);
411 appImplicitOutArgs.set_DfDp(p_index,
412 MEB::Derivative<Scalar>(my_dfdp_op_));
413 trans = Thyra::NOTRANS;
414 }
415 else if (dfdp_support.supports(MEB::DERIV_MV_JACOBIAN_FORM)) {
416 if (my_dfdp_mv_ == Teuchos::null)
417 my_dfdp_mv_ = Thyra::createMembers(
418 appImplicitModel_->get_f_space(),
419 appImplicitModel_->get_p_space(p_index)->dim());
420 appImplicitOutArgs.set_DfDp(
421 p_index, MEB::Derivative<Scalar>(my_dfdp_mv_,
422 MEB::DERIV_MV_JACOBIAN_FORM));
423 my_dfdp_op_ = my_dfdp_mv_;
424 trans = Thyra::NOTRANS;
425 }
426 else if (dfdp_support.supports(MEB::DERIV_MV_GRADIENT_FORM)) {
427 if (my_dfdp_mv_ == Teuchos::null)
428 my_dfdp_mv_ = Thyra::createMembers(
429 appImplicitModel_->get_p_space(p_index),
430 appImplicitModel_->get_f_space()->dim());
431 appImplicitOutArgs.set_DfDp(
432 p_index, MEB::Derivative<Scalar>(my_dfdp_mv_,
433 MEB::DERIV_MV_GRADIENT_FORM));
434 my_dfdp_op_ = my_dfdp_mv_;
435 trans = Thyra::TRANS;
436 }
437 else
438 TEUCHOS_TEST_FOR_EXCEPTION(
439 true, std::logic_error, "Invalid df/dp support");
440
441 appImplicitModel_->evalModel(appImplicitInArgs, appImplicitOutArgs);
442
443 // Add df/dy*dy/dp term to residual
444 RCP<DMVPV> f_dfdp = rcp_dynamic_cast<DMVPV>(outArgs.get_f(),true);
445 const int num_param = f_dfdp->getNonconstMultiVector()->domain()->dim()-1;
446 RCP<Thyra::MultiVectorBase<Scalar> > dfdp =
447 f_dfdp->getNonconstMultiVector()->subView(Range1D(1,num_param));
448 my_dfdp_op_->apply(trans, *dydp, dfdp.ptr(), Scalar(1.0), Scalar(1.0));
449 }
450}
451
452template <typename Scalar>
453Teuchos::RCP<const Teuchos::ParameterList>
455getValidParameters() const
456{
457 Teuchos::RCP<const Teuchos::ParameterList> fsa_pl =
459 Teuchos::RCP<Teuchos::ParameterList> pl = Teuchos::parameterList(*fsa_pl);
460 pl->set<int>("Sensitivity Y Tangent Index", 3);
461 return pl;
462}
463
464} // namespace Tempus
465
466#endif // Tempus_ModelEvaluatorPairPartIMEX_CombinedFSA_impl_hpp
ModelEvaluator pair for implicit and explicit (IMEX) evaulations.
void setup(const Teuchos::RCP< const Thyra::ModelEvaluator< Scalar > > &explicitModel, const Teuchos::RCP< const Thyra::ModelEvaluator< Scalar > > &implicitModel, int numExplicitOnlyBlocks=0, int parameterIndex=-1)
Setup ME when using default constructor – for derived classes.
Teuchos::RCP< const WrapperModelEvaluatorPairPartIMEX_Basic< Scalar > > forwardModel_
WrapperModelEvaluatorPairPartIMEX_CombinedFSA(const Teuchos::RCP< const WrapperModelEvaluatorPairPartIMEX_Basic< Scalar > > &forwardModel, const Teuchos::RCP< const Teuchos::ParameterList > &pList=Teuchos::null)
Constructor.
virtual void evalModelImpl(const Thyra::ModelEvaluatorBase::InArgs< Scalar > &inArgs, const Thyra::ModelEvaluatorBase::OutArgs< Scalar > &outArgs) const
virtual Teuchos::RCP< Thyra::VectorBase< Scalar > > getExplicitOnlyVector(const Teuchos::RCP< Thyra::VectorBase< Scalar > > &full) const
Extract explicit-only vector from a full solution vector.
virtual Teuchos::RCP< Thyra::VectorBase< Scalar > > getIMEXVector(const Teuchos::RCP< Thyra::VectorBase< Scalar > > &full) const
Extract IMEX vector from a full solution vector.
virtual Teuchos::RCP< const Thyra::VectorSpaceBase< Scalar > > get_p_space(int i) const
Get the p space.
virtual Teuchos::RCP< const Thyra::ModelEvaluator< Scalar > > getForwardModel() const
Get the underlying forward model.