Tempus  Version of the Day
Time Integration
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Macros Pages
Tempus_WrapperModelEvaluatorPairPartIMEX_Basic_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_Basic_impl_hpp
10 #define Tempus_ModelEvaluatorPairPartIMEX_Basic_impl_hpp
11 
12 #include "Thyra_ProductVectorBase.hpp"
13 #include "Thyra_ProductVectorSpaceBase.hpp"
14 
15 
16 namespace Tempus {
17 
18 template <typename Scalar>
21  : timeDer_(Teuchos::null), numExplicitOnlyBlocks_(0),
22  parameterIndex_(-1), useImplicitModel_(false)
23 {}
24 
25 template <typename Scalar>
28  const Teuchos::RCP<const Thyra::ModelEvaluator<Scalar> >& explicitModel,
29  const Teuchos::RCP<const Thyra::ModelEvaluator<Scalar> >& implicitModel,
30  int numExplicitOnlyBlocks, int parameterIndex)
31  : timeDer_(Teuchos::null), numExplicitOnlyBlocks_(numExplicitOnlyBlocks),
32  parameterIndex_(parameterIndex), useImplicitModel_(false)
33 {
34  setExplicitModel(explicitModel);
35  setImplicitModel(implicitModel);
37  initialize();
38 }
39 
40 template <typename Scalar>
41 void
44  const Teuchos::RCP<const Thyra::ModelEvaluator<Scalar> >& explicitModel,
45  const Teuchos::RCP<const Thyra::ModelEvaluator<Scalar> >& implicitModel,
46  int numExplicitOnlyBlocks, int parameterIndex)
47 {
48  timeDer_ = Teuchos::null;
49  numExplicitOnlyBlocks_ = numExplicitOnlyBlocks;
50  parameterIndex_ = parameterIndex;
51  useImplicitModel_ = false;
52  setExplicitModel(explicitModel);
53  setImplicitModel(implicitModel);
54  setParameterIndex();
55  initialize();
56 }
57 
58 template <typename Scalar>
59 void
62 {
63  using Teuchos::RCP;
64 
65  useImplicitModel_ = true;
66  wrapperImplicitInArgs_ = this->createInArgs();
67  wrapperImplicitOutArgs_ = this->createOutArgs();
68  useImplicitModel_ = false;
69 
70  RCP<const Thyra::VectorBase<Scalar> > z =
71  explicitModel_->getNominalValues().get_x();
72 
73  // A Thyra::VectorSpace requirement
74  TEUCHOS_TEST_FOR_EXCEPTION( !(getIMEXVector(z)->space()->isCompatible(
75  *(implicitModel_->get_x_space()))),
76  std::logic_error,
77  "Error - WrapperModelEvaluatorPairIMEX_Basic::initialize()\n"
78  " Explicit and Implicit vector x spaces are incompatible!\n"
79  " Explicit vector x space = " << *(getIMEXVector(z)->space()) << "\n"
80  " Implicit vector x space = " << *(implicitModel_->get_x_space())<< "\n");
81 
82  // Valid number of blocks?
83  RCP<const Thyra::ProductVectorBase<Scalar> > zPVector =
84  Teuchos::rcp_dynamic_cast<const Thyra::ProductVectorBase<Scalar> >(z);
85  TEUCHOS_TEST_FOR_EXCEPTION( zPVector == Teuchos::null, std::logic_error,
86  "Error - WrapperModelEvaluatorPairPartIMEX_Basic::initialize()\n"
87  " was given a VectorBase that could not be cast to a\n"
88  " ProductVectorBase!\n");
89 
90  int numBlocks = zPVector->productSpace()->numBlocks();
91 
92  TEUCHOS_TEST_FOR_EXCEPTION( !(0 <= numExplicitOnlyBlocks_ and
93  numExplicitOnlyBlocks_ < numBlocks),
94  std::logic_error,
95  "Error - WrapperModelEvaluatorPairPartIMEX_Basic::initialize()\n"
96  "Invalid number of explicit-only blocks = "<<numExplicitOnlyBlocks_<<"\n"
97  "Should be in the interval [0, numBlocks) = [0, "<<numBlocks<<")\n");
98 }
99 
100 template <typename Scalar>
101 void
104  const Teuchos::RCP<const Thyra::ModelEvaluator<Scalar> > & /* me */)
105 {
106  TEUCHOS_TEST_FOR_EXCEPTION( true, std::logic_error,
107  "Error - WrapperModelEvaluatorPairPartIMEX_Basic<Scalar>::setAppModel\n"
108  " should not be used. One should instead use setExplicitModel,\n"
109  " setImplicitModel, or create a new WrapperModelEvaluatorPairPartIMEX.\n");
110 }
111 
112 template <typename Scalar>
113 Teuchos::RCP<const Thyra::ModelEvaluator<Scalar> >
115 getAppModel() const
116 {
117  TEUCHOS_TEST_FOR_EXCEPTION( true, std::logic_error,
118  "Error - WrapperModelEvaluatorPairPartIMEX_Basic<Scalar>::getAppModel\n"
119  " should not be used. One should instead use getExplicitModel,\n"
120  " getImplicitModel, or directly use WrapperModelEvaluatorPairPartIMEX\n"
121  " object.\n");
122 }
123 
124 template <typename Scalar>
125 Teuchos::RCP<const Thyra::VectorSpaceBase<Scalar> >
127 get_x_space() const
128 {
129  if (useImplicitModel_ == true) return implicitModel_->get_x_space();
130 
131  return explicitModel_->get_x_space();
132 }
133 
134 template <typename Scalar>
135 Teuchos::RCP<const Thyra::VectorSpaceBase<Scalar> >
137 get_g_space(int i) const
138 {
139  if (useImplicitModel_ == true) return implicitModel_->get_g_space(i);
140 
141  return explicitModel_->get_g_space(i);
142 }
143 
144 template <typename Scalar>
145 Teuchos::RCP<const Thyra::VectorSpaceBase<Scalar> >
147 get_p_space(int i) const
148 {
149  if (useImplicitModel_ == true) return implicitModel_->get_p_space(i);
150 
151  return explicitModel_->get_p_space(i);
152 }
153 
154 template <typename Scalar>
155 void
158  const Teuchos::RCP<const Thyra::ModelEvaluator<Scalar> > & model )
159 {
160  implicitModel_ = model;
161 }
162 
163 template <typename Scalar>
164 Teuchos::RCP<Thyra::VectorBase<Scalar> >
166 getIMEXVector(const Teuchos::RCP<Thyra::VectorBase<Scalar> > & full) const
167 {
168  using Teuchos::RCP;
169  using Teuchos::rcp_dynamic_cast;
170 
171  Teuchos::RCP<Thyra::VectorBase<Scalar> > vector;
172  if(full == Teuchos::null) {
173  vector = Teuchos::null;
174  }
175  else if(numExplicitOnlyBlocks_ == 0) {
176  vector = full;
177  }
178  else {
179 
180  RCP<Thyra::ProductVectorBase<Scalar> > blk_full =
181  rcp_dynamic_cast<Thyra::ProductVectorBase<Scalar> >(full);
182  TEUCHOS_TEST_FOR_EXCEPTION( blk_full == Teuchos::null, std::logic_error,
183  "Error - WrapperModelEvaluatorPairPartIMEX_Basic::getIMEXVector()\n"
184  " was given a VectorBase that could not be cast to a\n"
185  " ProductVectorBase!\n");
186  int numBlocks = blk_full->productSpace()->numBlocks();
187 
188  // special case where the implicit terms are not blocked
189  if(numBlocks == numExplicitOnlyBlocks_+1)
190  vector = blk_full->getNonconstVectorBlock(numExplicitOnlyBlocks_);
191 
192  TEUCHOS_TEST_FOR_EXCEPTION(
193  !( numExplicitOnlyBlocks_ == 0 || full == Teuchos::null ||
194  numBlocks == numExplicitOnlyBlocks_+1 ),
195  std::logic_error, "Error - Invalid values!\n");
196  }
197 
198  return vector;
199 }
200 
201 template <typename Scalar>
202 Teuchos::RCP<const Thyra::VectorBase<Scalar> >
204 getIMEXVector(const Teuchos::RCP<const Thyra::VectorBase<Scalar> > & full) const
205 {
206  using Teuchos::RCP;
207  using Teuchos::rcp_dynamic_cast;
208 
209  Teuchos::RCP<const Thyra::VectorBase<Scalar> > vector;
210  if(full == Teuchos::null) {
211  vector = Teuchos::null;
212  }
213  else if(numExplicitOnlyBlocks_ == 0) {
214  vector = full;
215  }
216  else {
217 
218  // special case where the implicit terms are not blocked
219 
220  RCP<const Thyra::ProductVectorBase<Scalar> > blk_full =
221  rcp_dynamic_cast<const Thyra::ProductVectorBase<Scalar> >(full);
222  TEUCHOS_TEST_FOR_EXCEPTION( blk_full == Teuchos::null, std::logic_error,
223  "Error - WrapperModelEvaluatorPairPartIMEX_Basic::getIMEXVector()\n"
224  " was given a VectorBase that could not be cast to a\n"
225  " ProductVectorBase!\n");
226  int numBlocks = blk_full->productSpace()->numBlocks();
227 
228  // special case where the implicit terms are not blocked
229  if(numBlocks == numExplicitOnlyBlocks_+1)
230  vector = blk_full->getVectorBlock(numExplicitOnlyBlocks_);
231 
232  TEUCHOS_TEST_FOR_EXCEPTION(
233  !( numExplicitOnlyBlocks_ == 0 || full == Teuchos::null ||
234  numBlocks == numExplicitOnlyBlocks_+1 ),
235  std::logic_error, "Error - Invalid values!\n");
236  }
237 
238  return vector;
239 }
240 
241 template <typename Scalar>
242 Teuchos::RCP<Thyra::VectorBase<Scalar> >
245  const Teuchos::RCP<Thyra::VectorBase<Scalar> > & full) const
246 {
247  using Teuchos::RCP;
248  using Teuchos::rcp_dynamic_cast;
249 
250  Teuchos::RCP<Thyra::VectorBase<Scalar> > vector;
251  if(numExplicitOnlyBlocks_ == 0 || full == Teuchos::null) {
252  vector = Teuchos::null;
253  }
254  else if(numExplicitOnlyBlocks_ == 1) {
255 
256  // special case where the explicit terms are not blocked
257 
258  RCP<Thyra::ProductVectorBase<Scalar> > blk_full =
259  rcp_dynamic_cast<Thyra::ProductVectorBase<Scalar> >(full);
260  TEUCHOS_TEST_FOR_EXCEPTION( blk_full == Teuchos::null, std::logic_error,
261  "Error - WrapperModelEvaluatorPairPartIMEX_Basic::getExplicitOnlyVector\n"
262  " given a VectorBase that could not be cast to a ProductVectorBase!\n"
263  " full = " << *full << "\n");
264 
265  vector = blk_full->getNonconstVectorBlock(0);
266  }
267 
268  TEUCHOS_TEST_FOR_EXCEPTION(
269  !( (numExplicitOnlyBlocks_ == 0 || full == Teuchos::null) ||
270  (numExplicitOnlyBlocks_ == 1) ),
271  std::logic_error, "Error - Invalid values!\n");
272 
273  return vector;
274 }
275 
276 template <typename Scalar>
277 Teuchos::RCP<const Thyra::VectorBase<Scalar> >
280  const Teuchos::RCP<const Thyra::VectorBase<Scalar> > & full) const
281 {
282  using Teuchos::RCP;
283  using Teuchos::rcp_dynamic_cast;
284 
285  RCP<const Thyra::VectorBase<Scalar> > vector;
286  if(numExplicitOnlyBlocks_ == 0 || full == Teuchos::null) {
287  vector = Teuchos::null;
288  }
289  else if(numExplicitOnlyBlocks_ == 1) {
290 
291  // special case where the explicit terms are not blocked
292 
293  RCP<const Thyra::ProductVectorBase<Scalar> > blk_full =
294  rcp_dynamic_cast<const Thyra::ProductVectorBase<Scalar> >(full);
295  TEUCHOS_TEST_FOR_EXCEPTION( blk_full == Teuchos::null, std::logic_error,
296  "Error - WrapperModelEvaluatorPairPartIMEX_Basic::getExplicitOnlyVector\n"
297  " given a VectorBase that could not be cast to a ProductVectorBase!\n"
298  " full = " << *full << "\n");
299 
300  vector = blk_full->getVectorBlock(0);
301 
302  }
303 
304  TEUCHOS_TEST_FOR_EXCEPTION(
305  !( (numExplicitOnlyBlocks_ == 0 || full == Teuchos::null) ||
306  (numExplicitOnlyBlocks_ == 1) ),
307  std::logic_error, "Error - Invalid values!\n");
308 
309  return vector;
310 }
311 
312 template <typename Scalar>
313 void
315 setParameterIndex(int parameterIndex)
316 {
317  if (implicitModel_->Np() == 0) {
318  if (parameterIndex >= 0) {
319  Teuchos::RCP<Teuchos::FancyOStream> out = this->getOStream();
320  Teuchos::OSTab ostab(out,1,
321  "WrapperModelEvaluatorPairPartIMEX_Basic::setParameterIndex()");
322  *out << "Warning -- \n"
323  << " Invalid input parameter index = " << parameterIndex << "\n"
324  << " Should not be set since Np = "<<implicitModel_->Np() << "\n"
325  << " Not setting parameter index." << std::endl;
326  }
327  if (parameterIndex_ >= 0) {
328  Teuchos::RCP<Teuchos::FancyOStream> out = this->getOStream();
329  Teuchos::OSTab ostab(out,1,
330  "WrapperModelEvaluatorPairPartIMEX_Basic::setParameterIndex()");
331  *out << "Warning -- \n"
332  << " Invalid parameter index = " << parameterIndex_ << "\n"
333  << " Should not be set since Np = "<<implicitModel_->Np() << "\n"
334  << " Resetting to parameter index to unset state." << std::endl;
335  parameterIndex_ = -1;
336  }
337  } else {
338  if(parameterIndex >= 0) {
339  parameterIndex_ = parameterIndex;
340  } else if (parameterIndex_ < 0) {
341  parameterIndex_ = 0;
342  for(int i=0; i<implicitModel_->Np(); i++) {
343  if((*implicitModel_->get_p_names(i))[0] == "EXPLICIT_ONLY_VECTOR") {
344  parameterIndex_ = i;
345  break;
346  }
347  }
348  }
349  }
350 
351  return;
352 }
353 
354 template <typename Scalar>
355 Teuchos::RCP<const Thyra::VectorSpaceBase<Scalar> >
357 get_f_space() const
358 {
359  if (useImplicitModel_ == true) return implicitModel_->get_f_space();
360 
361  return explicitModel_->get_f_space();
362 }
363 
364 
365 template <typename Scalar>
366 Thyra::ModelEvaluatorBase::InArgs<Scalar>
369 {
370  typedef Thyra::ModelEvaluatorBase MEB;
371  MEB::InArgsSetup<Scalar> inArgs = this->createInArgs();
372  inArgs.set_Np(1);
373  return inArgs;
374 }
375 
376 template <typename Scalar>
377 Thyra::ModelEvaluatorBase::InArgs<Scalar>
380 {
381  typedef Thyra::ModelEvaluatorBase MEB;
382  MEB::InArgs<Scalar> implicitInArgs = implicitModel_->getNominalValues();
383  MEB::InArgs<Scalar> explicitInArgs = explicitModel_->getNominalValues();
384  const int np = std::max(implicitInArgs.Np(), explicitInArgs.Np());
385  if (useImplicitModel_ == true) {
386  MEB::InArgsSetup<Scalar> inArgs(implicitInArgs);
387  inArgs.setModelEvalDescription(this->description());
388  inArgs.set_Np(np);
389  return inArgs;
390  }
391 
392  MEB::InArgsSetup<Scalar> inArgs(explicitInArgs);
393  inArgs.setModelEvalDescription(this->description());
394  inArgs.set_Np(np);
395  return inArgs;
396 }
397 
398 template <typename Scalar>
399 Thyra::ModelEvaluatorBase::OutArgs<Scalar>
402 {
403  typedef Thyra::ModelEvaluatorBase MEB;
404  MEB::OutArgs<Scalar> implicitOutArgs = implicitModel_->createOutArgs();
405  MEB::OutArgs<Scalar> explicitOutArgs = explicitModel_->createOutArgs();
406  const int np = std::max(implicitOutArgs.Np(), explicitOutArgs.Np());
407  if (useImplicitModel_ == true) {
408  MEB::OutArgsSetup<Scalar> outArgs(implicitOutArgs);
409  outArgs.setModelEvalDescription(this->description());
410  outArgs.set_Np_Ng(np,implicitOutArgs.Ng());
411  return outArgs;
412  }
413 
414  MEB::OutArgsSetup<Scalar> outArgs(explicitOutArgs);
415  outArgs.setModelEvalDescription(this->description());
416  outArgs.set_Np_Ng(np,explicitOutArgs.Ng());
417  return outArgs;
418 }
419 
420 template <typename Scalar>
422 evalModelImpl(const Thyra::ModelEvaluatorBase::InArgs<Scalar> & inArgs,
423  const Thyra::ModelEvaluatorBase::OutArgs<Scalar> & outArgs) const
424 {
425  typedef Thyra::ModelEvaluatorBase MEB;
426  using Teuchos::RCP;
427 
428  RCP<const Thyra::VectorBase<Scalar> > x = inArgs.get_x();
429  RCP<Thyra::VectorBase<Scalar> > x_dot =
430  Thyra::createMember(implicitModel_->get_x_space());
431  timeDer_->compute(x, x_dot);
432 
433  MEB::InArgs<Scalar> appImplicitInArgs (wrapperImplicitInArgs_);
434  MEB::OutArgs<Scalar> appImplicitOutArgs(wrapperImplicitOutArgs_);
435  appImplicitInArgs.set_x(x);
436  appImplicitInArgs.set_x_dot(x_dot);
437  for (int i=0; i<implicitModel_->Np(); ++i) {
438  // Copy over parameters except for the parameter for explicit-only vector!
439  if ((inArgs.get_p(i) != Teuchos::null) and (i != parameterIndex_))
440  appImplicitInArgs.set_p(i, inArgs.get_p(i));
441  }
442 
443  appImplicitOutArgs.set_f(outArgs.get_f());
444  appImplicitOutArgs.set_W_op(outArgs.get_W_op());
445 
446  implicitModel_->evalModel(appImplicitInArgs,appImplicitOutArgs);
447 }
448 
449 } // end namespace Tempus
450 
451 #endif // Tempus_ModelEvaluatorPairPartIMEX_Basic_impl_hpp
virtual Teuchos::RCP< const Thyra::ModelEvaluator< Scalar > > getAppModel() const
Get the underlying application ModelEvaluator.
virtual Teuchos::RCP< const Thyra::VectorSpaceBase< Scalar > > get_g_space(int i) const
Get the g space.
virtual void setExplicitModel(const Teuchos::RCP< const Thyra::ModelEvaluator< Scalar > > &model)
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 Thyra::ModelEvaluatorBase::InArgs< Scalar > createInArgs() const
virtual void setParameterIndex(int parameterIndex=-1)
Set the parameter index for explicit-only vector.
virtual Teuchos::RCP< const Thyra::VectorSpaceBase< Scalar > > get_x_space() const
Get the x-solution space.
virtual Thyra::ModelEvaluatorBase::InArgs< Scalar > getNominalValues() const
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.
virtual Thyra::ModelEvaluatorBase::OutArgs< Scalar > createOutArgsImpl() const
virtual void setImplicitModel(const Teuchos::RCP< const Thyra::ModelEvaluator< Scalar > > &model)
virtual void setAppModel(const Teuchos::RCP< const Thyra::ModelEvaluator< Scalar > > &me)
Set the underlying application ModelEvaluator.
virtual Teuchos::RCP< Thyra::VectorBase< Scalar > > getIMEXVector(const Teuchos::RCP< Thyra::VectorBase< Scalar > > &full) const
Extract IMEX vector from a full solution vector.
WrapperModelEvaluatorPairPartIMEX_Basic()
Default constructor – Still requires setting the models and running initialize.
virtual Teuchos::RCP< const Thyra::VectorSpaceBase< Scalar > > get_f_space() const
virtual Teuchos::RCP< const Thyra::VectorSpaceBase< Scalar > > get_p_space(int i) const
Get the p space.
virtual void evalModelImpl(const Thyra::ModelEvaluatorBase::InArgs< Scalar > &in, const Thyra::ModelEvaluatorBase::OutArgs< Scalar > &out) const