Tempus  Version of the Day
Time Integration
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  if(full == Teuchos::null)
172  return Teuchos::null;
173 
174  if(numExplicitOnlyBlocks_==0)
175  return full;
176 
177  RCP<Thyra::ProductVectorBase<Scalar> > blk_full =
178  rcp_dynamic_cast<Thyra::ProductVectorBase<Scalar> >(full);
179  TEUCHOS_TEST_FOR_EXCEPTION( blk_full == Teuchos::null, std::logic_error,
180  "Error - WrapperModelEvaluatorPairPartIMEX_Basic::getIMEXVector()\n"
181  " was given a VectorBase that could not be cast to a\n"
182  " ProductVectorBase!\n");
183  int numBlocks = blk_full->productSpace()->numBlocks();
184 
185  // special case where the implicit terms are not blocked
186  if(numBlocks==numExplicitOnlyBlocks_+1)
187  return blk_full->getNonconstVectorBlock(numExplicitOnlyBlocks_);
188 
189  TEUCHOS_ASSERT(false);
190  return Teuchos::null;
191 }
192 
193 template <typename Scalar>
194 Teuchos::RCP<const Thyra::VectorBase<Scalar> >
196 getIMEXVector(const Teuchos::RCP<const Thyra::VectorBase<Scalar> > & full) const
197 {
198  using Teuchos::RCP;
199  using Teuchos::rcp_dynamic_cast;
200 
201  if(full == Teuchos::null)
202  return Teuchos::null;
203 
204  if(numExplicitOnlyBlocks_==0)
205  return full;
206 
207  RCP<const Thyra::ProductVectorBase<Scalar> > blk_full =
208  rcp_dynamic_cast<const Thyra::ProductVectorBase<Scalar> >(full);
209  TEUCHOS_TEST_FOR_EXCEPTION( blk_full == Teuchos::null, std::logic_error,
210  "Error - WrapperModelEvaluatorPairPartIMEX_Basic::getIMEXVector()\n"
211  " was given a VectorBase that could not be cast to a\n"
212  " ProductVectorBase!\n");
213  int numBlocks = blk_full->productSpace()->numBlocks();
214 
215  // special case where the implicit terms are not blocked
216  if(numBlocks==numExplicitOnlyBlocks_+1)
217  return blk_full->getVectorBlock(numExplicitOnlyBlocks_);
218 
219  TEUCHOS_ASSERT(false);
220  return Teuchos::null;
221 }
222 
223 template <typename Scalar>
224 Teuchos::RCP<Thyra::VectorBase<Scalar> >
227  const Teuchos::RCP<Thyra::VectorBase<Scalar> > & full) const
228 {
229  using Teuchos::RCP;
230  using Teuchos::rcp_dynamic_cast;
231 
232  if(numExplicitOnlyBlocks_ == 0 || full == Teuchos::null)
233  return Teuchos::null;
234 
235  RCP<Thyra::ProductVectorBase<Scalar> > blk_full =
236  rcp_dynamic_cast<Thyra::ProductVectorBase<Scalar> >(full);
237  TEUCHOS_TEST_FOR_EXCEPTION( blk_full == Teuchos::null, std::logic_error,
238  "Error - WrapperModelEvaluatorPairPartIMEX_Basic::getExplicitOnlyVector()\n"
239  " was given a VectorBase that could not be cast to a ProductVectorBase!\n"
240  " full = " << *full << "\n");
241 
242  // special case where the explicit terms are not blocked
243  if(numExplicitOnlyBlocks_==1)
244  return blk_full->getNonconstVectorBlock(0);
245 
246  TEUCHOS_ASSERT(false);
247  return Teuchos::null;
248 
249 }
250 
251 template <typename Scalar>
252 Teuchos::RCP<const Thyra::VectorBase<Scalar> >
255  const Teuchos::RCP<const Thyra::VectorBase<Scalar> > & full) const
256 {
257  using Teuchos::RCP;
258  using Teuchos::rcp_dynamic_cast;
259 
260  if(numExplicitOnlyBlocks_ == 0 || full == Teuchos::null)
261  return Teuchos::null;
262 
263  RCP<const Thyra::ProductVectorBase<Scalar> > blk_full =
264  rcp_dynamic_cast<const Thyra::ProductVectorBase<Scalar> >(full);
265  TEUCHOS_TEST_FOR_EXCEPTION( blk_full == Teuchos::null, std::logic_error,
266  "Error - WrapperModelEvaluatorPairPartIMEX_Basic::getExplicitOnlyVector()\n"
267  " was given a VectorBase that could not be cast to a ProductVectorBase!\n"
268  " full = " << *full << "\n");
269 
270  // special case where the explicit terms are not blocked
271  if(numExplicitOnlyBlocks_==1)
272  return blk_full->getVectorBlock(0);
273 
274  TEUCHOS_ASSERT(false);
275  return Teuchos::null;
276 
277 }
278 
279 template <typename Scalar>
280 void
282 setParameterIndex(int parameterIndex)
283 {
284  if (implicitModel_->Np() == 0) {
285  if (parameterIndex >= 0) {
286  Teuchos::RCP<Teuchos::FancyOStream> out = this->getOStream();
287  Teuchos::OSTab ostab(out,1,
288  "WrapperModelEvaluatorPairPartIMEX_Basic::setParameterIndex()");
289  *out << "Warning -- \n"
290  << " Invalid input parameter index = " << parameterIndex << "\n"
291  << " Should not be set since Np = "<<implicitModel_->Np() << "\n"
292  << " Not setting parameter index." << std::endl;
293  }
294  if (parameterIndex_ >= 0) {
295  Teuchos::RCP<Teuchos::FancyOStream> out = this->getOStream();
296  Teuchos::OSTab ostab(out,1,
297  "WrapperModelEvaluatorPairPartIMEX_Basic::setParameterIndex()");
298  *out << "Warning -- \n"
299  << " Invalid parameter index = " << parameterIndex_ << "\n"
300  << " Should not be set since Np = "<<implicitModel_->Np() << "\n"
301  << " Resetting to parameter index to unset state." << std::endl;
302  parameterIndex_ = -1;
303  }
304  } else {
305  if(parameterIndex >= 0) {
306  parameterIndex_ = parameterIndex;
307  } else if (parameterIndex_ < 0) {
308  parameterIndex_ = 0;
309  for(int i=0; i<implicitModel_->Np(); i++) {
310  if((*implicitModel_->get_p_names(i))[0] == "EXPLICIT_ONLY_VECTOR") {
311  parameterIndex_ = i;
312  break;
313  }
314  }
315  }
316  }
317 
318  return;
319 }
320 
321 template <typename Scalar>
322 Teuchos::RCP<const Thyra::VectorSpaceBase<Scalar> >
324 get_f_space() const
325 {
326  if (useImplicitModel_ == true) return implicitModel_->get_f_space();
327 
328  return explicitModel_->get_f_space();
329 }
330 
331 
332 template <typename Scalar>
333 Thyra::ModelEvaluatorBase::InArgs<Scalar>
336 {
337  typedef Thyra::ModelEvaluatorBase MEB;
338  MEB::InArgsSetup<Scalar> inArgs = this->createInArgs();
339  inArgs.set_Np(1);
340  return inArgs;
341 }
342 
343 template <typename Scalar>
344 Thyra::ModelEvaluatorBase::InArgs<Scalar>
347 {
348  typedef Thyra::ModelEvaluatorBase MEB;
349  MEB::InArgs<Scalar> implicitInArgs = implicitModel_->getNominalValues();
350  MEB::InArgs<Scalar> explicitInArgs = explicitModel_->getNominalValues();
351  const int np = std::max(implicitInArgs.Np(), explicitInArgs.Np());
352  if (useImplicitModel_ == true) {
353  MEB::InArgsSetup<Scalar> inArgs(implicitInArgs);
354  inArgs.setModelEvalDescription(this->description());
355  inArgs.set_Np(np);
356  return inArgs;
357  }
358 
359  MEB::InArgsSetup<Scalar> inArgs(explicitInArgs);
360  inArgs.setModelEvalDescription(this->description());
361  inArgs.set_Np(np);
362  return inArgs;
363 }
364 
365 template <typename Scalar>
366 Thyra::ModelEvaluatorBase::OutArgs<Scalar>
369 {
370  typedef Thyra::ModelEvaluatorBase MEB;
371  MEB::OutArgs<Scalar> implicitOutArgs = implicitModel_->createOutArgs();
372  MEB::OutArgs<Scalar> explicitOutArgs = explicitModel_->createOutArgs();
373  const int np = std::max(implicitOutArgs.Np(), explicitOutArgs.Np());
374  if (useImplicitModel_ == true) {
375  MEB::OutArgsSetup<Scalar> outArgs(implicitOutArgs);
376  outArgs.setModelEvalDescription(this->description());
377  outArgs.set_Np_Ng(np,implicitOutArgs.Ng());
378  return outArgs;
379  }
380 
381  MEB::OutArgsSetup<Scalar> outArgs(explicitOutArgs);
382  outArgs.setModelEvalDescription(this->description());
383  outArgs.set_Np_Ng(np,explicitOutArgs.Ng());
384  return outArgs;
385 }
386 
387 template <typename Scalar>
389 evalModelImpl(const Thyra::ModelEvaluatorBase::InArgs<Scalar> & inArgs,
390  const Thyra::ModelEvaluatorBase::OutArgs<Scalar> & outArgs) const
391 {
392  typedef Thyra::ModelEvaluatorBase MEB;
393  using Teuchos::RCP;
394 
395  RCP<const Thyra::VectorBase<Scalar> > x = inArgs.get_x();
396  RCP<Thyra::VectorBase<Scalar> > x_dot =
397  Thyra::createMember(implicitModel_->get_x_space());
398  timeDer_->compute(x, x_dot);
399 
400  MEB::InArgs<Scalar> appImplicitInArgs (wrapperImplicitInArgs_);
401  MEB::OutArgs<Scalar> appImplicitOutArgs(wrapperImplicitOutArgs_);
402  appImplicitInArgs.set_x(x);
403  appImplicitInArgs.set_x_dot(x_dot);
404  for (int i=0; i<implicitModel_->Np(); ++i) {
405  // Copy over parameters except for the parameter for explicit-only vector!
406  if ((inArgs.get_p(i) != Teuchos::null) and (i != parameterIndex_))
407  appImplicitInArgs.set_p(i, inArgs.get_p(i));
408  }
409 
410  appImplicitOutArgs.set_f(outArgs.get_f());
411  appImplicitOutArgs.set_W_op(outArgs.get_W_op());
412 
413  implicitModel_->evalModel(appImplicitInArgs,appImplicitOutArgs);
414 }
415 
416 } // end namespace Tempus
417 
418 #endif // Tempus_ModelEvaluatorPairPartIMEX_Basic_impl_hpp
Tempus::WrapperModelEvaluatorPairPartIMEX_Basic::WrapperModelEvaluatorPairPartIMEX_Basic
WrapperModelEvaluatorPairPartIMEX_Basic()
Default constructor – Still requires setting the models and running initialize.
Definition: Tempus_WrapperModelEvaluatorPairPartIMEX_Basic_impl.hpp:20
Tempus::WrapperModelEvaluatorPairPartIMEX_Basic::setup
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.
Definition: Tempus_WrapperModelEvaluatorPairPartIMEX_Basic_impl.hpp:43
Tempus::WrapperModelEvaluatorPairPartIMEX_Basic::createOutArgsImpl
virtual Thyra::ModelEvaluatorBase::OutArgs< Scalar > createOutArgsImpl() const
Definition: Tempus_WrapperModelEvaluatorPairPartIMEX_Basic_impl.hpp:368
Tempus::WrapperModelEvaluatorPairPartIMEX_Basic::get_p_space
virtual Teuchos::RCP< const Thyra::VectorSpaceBase< Scalar > > get_p_space(int i) const
Get the p space.
Definition: Tempus_WrapperModelEvaluatorPairPartIMEX_Basic_impl.hpp:147
Tempus::WrapperModelEvaluatorPairPartIMEX_Basic::getExplicitOnlyVector
virtual Teuchos::RCP< Thyra::VectorBase< Scalar > > getExplicitOnlyVector(const Teuchos::RCP< Thyra::VectorBase< Scalar > > &full) const
Extract explicit-only vector from a full solution vector.
Definition: Tempus_WrapperModelEvaluatorPairPartIMEX_Basic_impl.hpp:226
Tempus
Definition: Tempus_AdjointAuxSensitivityModelEvaluator_decl.hpp:20
Tempus::WrapperModelEvaluatorPairPartIMEX_Basic::get_g_space
virtual Teuchos::RCP< const Thyra::VectorSpaceBase< Scalar > > get_g_space(int i) const
Get the g space.
Definition: Tempus_WrapperModelEvaluatorPairPartIMEX_Basic_impl.hpp:137
Tempus::WrapperModelEvaluatorPairPartIMEX_Basic::initialize
virtual void initialize()
Initialize after setting member data.
Definition: Tempus_WrapperModelEvaluatorPairPartIMEX_Basic_impl.hpp:61
Tempus::WrapperModelEvaluatorPairPartIMEX_Basic::get_f_space
virtual Teuchos::RCP< const Thyra::VectorSpaceBase< Scalar > > get_f_space() const
Definition: Tempus_WrapperModelEvaluatorPairPartIMEX_Basic_impl.hpp:324
Tempus::WrapperModelEvaluatorPairPartIMEX_Basic::setExplicitModel
virtual void setExplicitModel(const Teuchos::RCP< const Thyra::ModelEvaluator< Scalar > > &model)
Definition: Tempus_WrapperModelEvaluatorPairPartIMEX_Basic_decl.hpp:111
Tempus::WrapperModelEvaluatorPairPartIMEX_Basic::getAppModel
virtual Teuchos::RCP< const Thyra::ModelEvaluator< Scalar > > getAppModel() const
Get the underlying application ModelEvaluator.
Definition: Tempus_WrapperModelEvaluatorPairPartIMEX_Basic_impl.hpp:115
Tempus::WrapperModelEvaluatorPairPartIMEX_Basic::setAppModel
virtual void setAppModel(const Teuchos::RCP< const Thyra::ModelEvaluator< Scalar > > &me)
Set the underlying application ModelEvaluator.
Definition: Tempus_WrapperModelEvaluatorPairPartIMEX_Basic_impl.hpp:103
Tempus::WrapperModelEvaluatorPairPartIMEX_Basic::setParameterIndex
virtual void setParameterIndex(int parameterIndex=-1)
Set the parameter index for explicit-only vector.
Definition: Tempus_WrapperModelEvaluatorPairPartIMEX_Basic_impl.hpp:282
Tempus::WrapperModelEvaluatorPairPartIMEX_Basic::getNominalValues
virtual Thyra::ModelEvaluatorBase::InArgs< Scalar > getNominalValues() const
Definition: Tempus_WrapperModelEvaluatorPairPartIMEX_Basic_impl.hpp:335
Tempus::WrapperModelEvaluatorPairPartIMEX_Basic::createInArgs
virtual Thyra::ModelEvaluatorBase::InArgs< Scalar > createInArgs() const
Definition: Tempus_WrapperModelEvaluatorPairPartIMEX_Basic_impl.hpp:346
Tempus::WrapperModelEvaluatorPairPartIMEX_Basic::getIMEXVector
virtual Teuchos::RCP< Thyra::VectorBase< Scalar > > getIMEXVector(const Teuchos::RCP< Thyra::VectorBase< Scalar > > &full) const
Extract IMEX vector from a full solution vector.
Definition: Tempus_WrapperModelEvaluatorPairPartIMEX_Basic_impl.hpp:166
Teuchos
Definition: Tempus_Integrator.hpp:19
Tempus::WrapperModelEvaluatorPairPartIMEX_Basic::get_x_space
virtual Teuchos::RCP< const Thyra::VectorSpaceBase< Scalar > > get_x_space() const
Get the x-solution space.
Definition: Tempus_WrapperModelEvaluatorPairPartIMEX_Basic_impl.hpp:127
Tempus::WrapperModelEvaluatorPairPartIMEX_Basic::setImplicitModel
virtual void setImplicitModel(const Teuchos::RCP< const Thyra::ModelEvaluator< Scalar > > &model)
Definition: Tempus_WrapperModelEvaluatorPairPartIMEX_Basic_impl.hpp:157
Tempus::WrapperModelEvaluatorPairPartIMEX_Basic::evalModelImpl
virtual void evalModelImpl(const Thyra::ModelEvaluatorBase::InArgs< Scalar > &in, const Thyra::ModelEvaluatorBase::OutArgs< Scalar > &out) const
Definition: Tempus_WrapperModelEvaluatorPairPartIMEX_Basic_impl.hpp:389