Tempus Version of the Day
Time Integration
Loading...
Searching...
No Matches
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
16namespace Tempus {
17
18template <typename Scalar>
21 : timeDer_(Teuchos::null), numExplicitOnlyBlocks_(0),
22 parameterIndex_(-1), useImplicitModel_(false)
23{}
24
25template <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
40template <typename Scalar>
41void
43setup(
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
58template <typename Scalar>
59void
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_ &&
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
100template <typename Scalar>
101void
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
112template <typename Scalar>
113Teuchos::RCP<const Thyra::ModelEvaluator<Scalar> >
115getAppModel() 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
124template <typename Scalar>
125Teuchos::RCP<const Thyra::VectorSpaceBase<Scalar> >
127get_x_space() const
128{
129 if (useImplicitModel_ == true) return implicitModel_->get_x_space();
130
131 return explicitModel_->get_x_space();
132}
133
134template <typename Scalar>
135Teuchos::RCP<const Thyra::VectorSpaceBase<Scalar> >
137get_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
144template <typename Scalar>
145Teuchos::RCP<const Thyra::VectorSpaceBase<Scalar> >
147get_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
154template <typename Scalar>
155void
158 const Teuchos::RCP<const Thyra::ModelEvaluator<Scalar> > & model )
159{
160 implicitModel_ = model;
161}
162
163template <typename Scalar>
164Teuchos::RCP<Thyra::VectorBase<Scalar> >
166getIMEXVector(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
201template <typename Scalar>
202Teuchos::RCP<const Thyra::VectorBase<Scalar> >
204getIMEXVector(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
241template <typename Scalar>
242Teuchos::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
276template <typename Scalar>
277Teuchos::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
312template <typename Scalar>
313void
315setParameterIndex(int parameterIndex)
316{
317 if (implicitModel_->Np() == 0) {
318 if (parameterIndex >= 0) {
319 Teuchos::RCP<Teuchos::FancyOStream> out = this->getOStream();
320 out->setOutputToRootOnly(0);
321 Teuchos::OSTab ostab(out,1,
322 "WrapperModelEvaluatorPairPartIMEX_Basic::setParameterIndex()");
323 *out << "Warning -- \n"
324 << " Invalid input parameter index = " << parameterIndex << "\n"
325 << " Should not be set since Np = "<<implicitModel_->Np() << "\n"
326 << " Not setting parameter index." << std::endl;
327 }
328 if (parameterIndex_ >= 0) {
329 Teuchos::RCP<Teuchos::FancyOStream> out = this->getOStream();
330 out->setOutputToRootOnly(0);
331 Teuchos::OSTab ostab(out,1,
332 "WrapperModelEvaluatorPairPartIMEX_Basic::setParameterIndex()");
333 *out << "Warning -- \n"
334 << " Invalid parameter index = " << parameterIndex_ << "\n"
335 << " Should not be set since Np = "<<implicitModel_->Np() << "\n"
336 << " Resetting to parameter index to unset state." << std::endl;
337 parameterIndex_ = -1;
338 }
339 } else {
340 if(parameterIndex >= 0) {
341 parameterIndex_ = parameterIndex;
342 } else if (parameterIndex_ < 0) {
343 parameterIndex_ = 0;
344 for(int i=0; i<implicitModel_->Np(); i++) {
345 if((*implicitModel_->get_p_names(i))[0] == "EXPLICIT_ONLY_VECTOR") {
346 parameterIndex_ = i;
347 break;
348 }
349 }
350 }
351 }
352
353 return;
354}
355
356template <typename Scalar>
357Teuchos::RCP<const Thyra::VectorSpaceBase<Scalar> >
359get_f_space() const
360{
361 if (useImplicitModel_ == true) return implicitModel_->get_f_space();
362
363 return explicitModel_->get_f_space();
364}
365
366template <typename Scalar>
367Thyra::ModelEvaluatorBase::InArgs<Scalar>
369getNominalValues() const
370{
371 typedef Thyra::ModelEvaluatorBase MEB;
372 MEB::InArgsSetup<Scalar> inArgs = this->createInArgs();
373 inArgs.set_Np(1);
374 return std::move(inArgs);
375}
376
377template <typename Scalar>
378Thyra::ModelEvaluatorBase::InArgs<Scalar>
380createInArgs() const
381{
382 typedef Thyra::ModelEvaluatorBase MEB;
383 MEB::InArgs<Scalar> implicitInArgs = implicitModel_->getNominalValues();
384 MEB::InArgs<Scalar> explicitInArgs = explicitModel_->getNominalValues();
385 const int np = std::max(implicitInArgs.Np(), explicitInArgs.Np());
386 if (useImplicitModel_ == true) {
387 MEB::InArgsSetup<Scalar> inArgs(implicitInArgs);
388 inArgs.setModelEvalDescription(this->description());
389 inArgs.set_Np(np);
390 return std::move(inArgs);
391 }
392
393 MEB::InArgsSetup<Scalar> inArgs(explicitInArgs);
394 inArgs.setModelEvalDescription(this->description());
395 inArgs.set_Np(np);
396 return std::move(inArgs);
397}
398
399template <typename Scalar>
400Thyra::ModelEvaluatorBase::OutArgs<Scalar>
402createOutArgsImpl() const
403{
404 typedef Thyra::ModelEvaluatorBase MEB;
405 MEB::OutArgs<Scalar> implicitOutArgs = implicitModel_->createOutArgs();
406 MEB::OutArgs<Scalar> explicitOutArgs = explicitModel_->createOutArgs();
407 const int np = std::max(implicitOutArgs.Np(), explicitOutArgs.Np());
408 if (useImplicitModel_ == true) {
409 MEB::OutArgsSetup<Scalar> outArgs(implicitOutArgs);
410 outArgs.setModelEvalDescription(this->description());
411 outArgs.set_Np_Ng(np,implicitOutArgs.Ng());
412 return std::move(outArgs);
413 }
414
415 MEB::OutArgsSetup<Scalar> outArgs(explicitOutArgs);
416 outArgs.setModelEvalDescription(this->description());
417 outArgs.set_Np_Ng(np,explicitOutArgs.Ng());
418 return std::move(outArgs);
419}
420
421template <typename Scalar>
423evalModelImpl(const Thyra::ModelEvaluatorBase::InArgs<Scalar> & inArgs,
424 const Thyra::ModelEvaluatorBase::OutArgs<Scalar> & outArgs) const
425{
426 typedef Thyra::ModelEvaluatorBase MEB;
427 using Teuchos::RCP;
428
429 RCP<const Thyra::VectorBase<Scalar> > x = inArgs.get_x();
430 RCP<Thyra::VectorBase<Scalar> > x_dot =
431 Thyra::createMember(implicitModel_->get_x_space());
432 timeDer_->compute(x, x_dot);
433
434 MEB::InArgs<Scalar> appImplicitInArgs (wrapperImplicitInArgs_);
435 MEB::OutArgs<Scalar> appImplicitOutArgs(wrapperImplicitOutArgs_);
436 appImplicitInArgs.set_x(x);
437 appImplicitInArgs.set_x_dot(x_dot);
438 for (int i=0; i<implicitModel_->Np(); ++i) {
439 // Copy over parameters except for the parameter for explicit-only vector!
440 if ((inArgs.get_p(i) != Teuchos::null) && (i != parameterIndex_))
441 appImplicitInArgs.set_p(i, inArgs.get_p(i));
442 }
443
444 appImplicitOutArgs.set_f(outArgs.get_f());
445 appImplicitOutArgs.set_W_op(outArgs.get_W_op());
446
447 implicitModel_->evalModel(appImplicitInArgs,appImplicitOutArgs);
448}
449
450} // end namespace Tempus
451
452#endif // Tempus_ModelEvaluatorPairPartIMEX_Basic_impl_hpp
virtual Teuchos::RCP< Thyra::VectorBase< Scalar > > getIMEXVector(const Teuchos::RCP< Thyra::VectorBase< Scalar > > &full) const
Extract IMEX vector from a full solution vector.
virtual void setParameterIndex(int parameterIndex=-1)
Set the parameter index for explicit-only vector.
virtual Teuchos::RCP< Thyra::VectorBase< Scalar > > getExplicitOnlyVector(const Teuchos::RCP< Thyra::VectorBase< Scalar > > &full) const
Extract explicit-only vector from a full solution vector.
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 void evalModelImpl(const Thyra::ModelEvaluatorBase::InArgs< Scalar > &in, const Thyra::ModelEvaluatorBase::OutArgs< Scalar > &out) const
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 Teuchos::RCP< const Thyra::VectorSpaceBase< Scalar > > get_g_space(int i) const
Get the g space.
virtual Thyra::ModelEvaluatorBase::InArgs< Scalar > getNominalValues() const
virtual Teuchos::RCP< const Thyra::ModelEvaluator< Scalar > > getAppModel() const
Get the underlying application ModelEvaluator.
virtual Teuchos::RCP< const Thyra::VectorSpaceBase< Scalar > > get_x_space() const
Get the x-solution space.
virtual Thyra::ModelEvaluatorBase::OutArgs< Scalar > createOutArgsImpl() const
virtual void setAppModel(const Teuchos::RCP< const Thyra::ModelEvaluator< Scalar > > &me)
Set the underlying application ModelEvaluator.
WrapperModelEvaluatorPairPartIMEX_Basic()
Default constructor – Still requires setting the models and running initialize.
virtual void setExplicitModel(const Teuchos::RCP< const Thyra::ModelEvaluator< Scalar > > &model)
virtual void setImplicitModel(const Teuchos::RCP< const Thyra::ModelEvaluator< Scalar > > &model)
virtual Thyra::ModelEvaluatorBase::InArgs< Scalar > createInArgs() const