9 #ifndef Tempus_ModelEvaluatorPairPartIMEX_Basic_impl_hpp
10 #define Tempus_ModelEvaluatorPairPartIMEX_Basic_impl_hpp
12 #include "Thyra_ProductVectorBase.hpp"
13 #include "Thyra_ProductVectorSpaceBase.hpp"
17 template <
typename Scalar>
18 WrapperModelEvaluatorPairPartIMEX_Basic<
20 : timeDer_(Teuchos::null),
21 numExplicitOnlyBlocks_(0),
23 useImplicitModel_(false)
27 template <
typename Scalar>
32 int numExplicitOnlyBlocks,
int parameterIndex)
33 : timeDer_(Teuchos::null),
34 numExplicitOnlyBlocks_(numExplicitOnlyBlocks),
35 parameterIndex_(parameterIndex),
36 useImplicitModel_(false)
44 template <
typename Scalar>
48 int numExplicitOnlyBlocks,
int parameterIndex)
50 timeDer_ = Teuchos::null;
51 numExplicitOnlyBlocks_ = numExplicitOnlyBlocks;
52 parameterIndex_ = parameterIndex;
53 useImplicitModel_ =
false;
54 setExplicitModel(explicitModel);
55 setImplicitModel(implicitModel);
60 template <
typename Scalar>
65 useImplicitModel_ =
true;
66 wrapperImplicitInArgs_ = this->createInArgs();
67 wrapperImplicitOutArgs_ = this->createOutArgs();
68 useImplicitModel_ =
false;
70 RCP<const Thyra::VectorBase<Scalar> > z =
71 explicitModel_->getNominalValues().get_x();
75 !(getIMEXVector(z)->space()->isCompatible(
76 *(implicitModel_->get_x_space()))),
78 "Error - WrapperModelEvaluatorPairIMEX_Basic::initialize()\n"
79 <<
" Explicit and Implicit vector x spaces are incompatible!\n"
80 <<
" Explicit vector x space = "
81 << *(getIMEXVector(z)->space())
82 <<
"\n Implicit vector x space = "
83 << *(implicitModel_->get_x_space()) <<
"\n");
86 RCP<const Thyra::ProductVectorBase<Scalar> > zPVector =
89 zPVector == Teuchos::null, std::logic_error,
90 "Error - WrapperModelEvaluatorPairPartIMEX_Basic::initialize()\n"
91 " was given a VectorBase that could not be cast to a\n"
92 " ProductVectorBase!\n");
94 int numBlocks = zPVector->productSpace()->numBlocks();
96 TEUCHOS_TEST_FOR_EXCEPTION(
97 !(0 <= numExplicitOnlyBlocks_ && numExplicitOnlyBlocks_ < numBlocks),
99 "Error - WrapperModelEvaluatorPairPartIMEX_Basic::initialize()\n"
100 <<
"Invalid number of explicit-only blocks = "
101 << numExplicitOnlyBlocks_
102 <<
"\nShould be in the interval [0, numBlocks) = [0, "
103 << numBlocks <<
")\n");
106 template <
typename Scalar>
111 true, std::logic_error,
112 "Error - WrapperModelEvaluatorPairPartIMEX_Basic<Scalar>::setAppModel\n"
113 " should not be used. One should instead use setExplicitModel,\n"
114 " setImplicitModel, or create a new "
115 "WrapperModelEvaluatorPairPartIMEX.\n");
118 template <
typename Scalar>
123 true, std::logic_error,
124 "Error - WrapperModelEvaluatorPairPartIMEX_Basic<Scalar>::getAppModel\n"
125 " should not be used. One should instead use getExplicitModel,\n"
126 " getImplicitModel, or directly use WrapperModelEvaluatorPairPartIMEX\n"
130 template <
typename Scalar>
134 if (useImplicitModel_ ==
true)
return implicitModel_->get_x_space();
136 return explicitModel_->get_x_space();
139 template <
typename Scalar>
143 if (useImplicitModel_ ==
true)
return implicitModel_->get_g_space(i);
145 return explicitModel_->get_g_space(i);
148 template <
typename Scalar>
152 if (useImplicitModel_ ==
true)
return implicitModel_->get_p_space(i);
154 return explicitModel_->get_p_space(i);
157 template <
typename Scalar>
161 implicitModel_ = model;
164 template <
typename Scalar>
170 using Teuchos::rcp_dynamic_cast;
173 if (full == Teuchos::null) {
174 vector = Teuchos::null;
176 else if (numExplicitOnlyBlocks_ == 0) {
180 RCP<Thyra::ProductVectorBase<Scalar> > blk_full =
183 blk_full == Teuchos::null, std::logic_error,
184 "Error - WrapperModelEvaluatorPairPartIMEX_Basic::getIMEXVector()\n"
185 " was given a VectorBase that could not be cast to a\n"
186 " ProductVectorBase!\n");
187 int numBlocks = blk_full->productSpace()->numBlocks();
190 if (numBlocks == numExplicitOnlyBlocks_ + 1)
191 vector = blk_full->getNonconstVectorBlock(numExplicitOnlyBlocks_);
193 TEUCHOS_TEST_FOR_EXCEPTION(
194 !(numExplicitOnlyBlocks_ == 0 || full == Teuchos::null ||
195 numBlocks == numExplicitOnlyBlocks_ + 1),
196 std::logic_error,
"Error - Invalid values!\n");
202 template <
typename Scalar>
208 using Teuchos::rcp_dynamic_cast;
211 if (full == Teuchos::null) {
212 vector = Teuchos::null;
214 else if (numExplicitOnlyBlocks_ == 0) {
220 RCP<const Thyra::ProductVectorBase<Scalar> > blk_full =
223 blk_full == Teuchos::null, std::logic_error,
224 "Error - WrapperModelEvaluatorPairPartIMEX_Basic::getIMEXVector()\n"
225 " was given a VectorBase that could not be cast to a\n"
226 " ProductVectorBase!\n");
227 int numBlocks = blk_full->productSpace()->numBlocks();
230 if (numBlocks == numExplicitOnlyBlocks_ + 1)
231 vector = blk_full->getVectorBlock(numExplicitOnlyBlocks_);
233 TEUCHOS_TEST_FOR_EXCEPTION(
234 !(numExplicitOnlyBlocks_ == 0 || full == Teuchos::null ||
235 numBlocks == numExplicitOnlyBlocks_ + 1),
236 std::logic_error,
"Error - Invalid values!\n");
242 template <
typename Scalar>
248 using Teuchos::rcp_dynamic_cast;
251 if (numExplicitOnlyBlocks_ == 0 || full == Teuchos::null) {
252 vector = Teuchos::null;
254 else if (numExplicitOnlyBlocks_ == 1) {
257 RCP<Thyra::ProductVectorBase<Scalar> > blk_full =
260 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");
265 vector = blk_full->getNonconstVectorBlock(0);
269 !((numExplicitOnlyBlocks_ == 0 || full == Teuchos::null) ||
270 (numExplicitOnlyBlocks_ == 1)),
271 std::logic_error,
"Error - Invalid values!\n");
276 template <
typename Scalar>
282 using Teuchos::rcp_dynamic_cast;
284 RCP<const Thyra::VectorBase<Scalar> > vector;
285 if (numExplicitOnlyBlocks_ == 0 || full == Teuchos::null) {
286 vector = Teuchos::null;
288 else if (numExplicitOnlyBlocks_ == 1) {
291 RCP<const Thyra::ProductVectorBase<Scalar> > blk_full =
294 blk_full == Teuchos::null, std::logic_error,
295 "Error - WrapperModelEvaluatorPairPartIMEX_Basic::getExplicitOnlyVector\n"
296 <<
" given a VectorBase that could not be cast to a ProductVectorBase!\n"
297 <<
" full = " << *full <<
"\n");
299 vector = blk_full->getVectorBlock(0);
303 !((numExplicitOnlyBlocks_ == 0 || full == Teuchos::null) ||
304 (numExplicitOnlyBlocks_ == 1)),
305 std::logic_error,
"Error - Invalid values!\n");
310 template <
typename Scalar>
314 if (implicitModel_->Np() == 0) {
315 if (parameterIndex >= 0) {
320 "WrapperModelEvaluatorPairPartIMEX_Basic::setParameterIndex()");
321 *out <<
"Warning -- \n"
322 <<
" Invalid input parameter index = " << parameterIndex <<
"\n"
323 <<
" Should not be set since Np = " << implicitModel_->Np() <<
"\n"
324 <<
" Not setting parameter index." << std::endl;
326 if (parameterIndex_ >= 0) {
331 "WrapperModelEvaluatorPairPartIMEX_Basic::setParameterIndex()");
332 *out <<
"Warning -- \n"
333 <<
" Invalid parameter index = " << parameterIndex_ <<
"\n"
334 <<
" Should not be set since Np = " << implicitModel_->Np() <<
"\n"
335 <<
" Resetting to parameter index to unset state." << std::endl;
336 parameterIndex_ = -1;
340 if (parameterIndex >= 0) {
341 parameterIndex_ = parameterIndex;
343 else if (parameterIndex_ < 0) {
345 for (
int i = 0; i < implicitModel_->Np(); i++) {
346 if ((*implicitModel_->get_p_names(i))[0] ==
"EXPLICIT_ONLY_VECTOR") {
357 template <
typename Scalar>
361 if (useImplicitModel_ ==
true)
return implicitModel_->get_f_space();
363 return explicitModel_->get_f_space();
366 template <
typename Scalar>
371 MEB::InArgsSetup<Scalar> inArgs = this->createInArgs();
372 return std::move(inArgs);
375 template <
typename Scalar>
380 MEB::InArgs<Scalar> implicitInArgs = implicitModel_->getNominalValues();
381 MEB::InArgs<Scalar> explicitInArgs = explicitModel_->getNominalValues();
382 const int np = std::max(implicitInArgs.Np(), explicitInArgs.Np());
383 if (useImplicitModel_ ==
true) {
384 MEB::InArgsSetup<Scalar> inArgs(implicitInArgs);
385 inArgs.setModelEvalDescription(this->description());
387 return std::move(inArgs);
390 MEB::InArgsSetup<Scalar> inArgs(explicitInArgs);
391 inArgs.setModelEvalDescription(this->description());
393 return std::move(inArgs);
396 template <
typename Scalar>
401 MEB::OutArgs<Scalar> implicitOutArgs = implicitModel_->createOutArgs();
402 MEB::OutArgs<Scalar> explicitOutArgs = explicitModel_->createOutArgs();
403 const int np = std::max(implicitOutArgs.Np(), explicitOutArgs.Np());
404 if (useImplicitModel_ ==
true) {
405 MEB::OutArgsSetup<Scalar> outArgs(implicitOutArgs);
406 outArgs.setModelEvalDescription(this->description());
407 outArgs.set_Np_Ng(np, implicitOutArgs.Ng());
408 return std::move(outArgs);
411 MEB::OutArgsSetup<Scalar> outArgs(explicitOutArgs);
412 outArgs.setModelEvalDescription(this->description());
413 outArgs.set_Np_Ng(np, explicitOutArgs.Ng());
414 return std::move(outArgs);
417 template <
typename Scalar>
425 RCP<const Thyra::VectorBase<Scalar> > x = inArgs.
get_x();
426 RCP<Thyra::VectorBase<Scalar> > x_dot =
427 Thyra::createMember(implicitModel_->get_x_space());
428 timeDer_->compute(x, x_dot);
430 MEB::InArgs<Scalar> appImplicitInArgs(wrapperImplicitInArgs_);
431 MEB::OutArgs<Scalar> appImplicitOutArgs(wrapperImplicitOutArgs_);
432 appImplicitInArgs.set_x(x);
433 appImplicitInArgs.set_x_dot(x_dot);
434 for (
int i = 0; i < implicitModel_->Np(); ++i) {
436 if ((inArgs.
get_p(i) != Teuchos::null) && (i != parameterIndex_))
437 appImplicitInArgs.set_p(i, inArgs.
get_p(i));
440 appImplicitOutArgs.set_f(outArgs.
get_f());
441 appImplicitOutArgs.set_W_op(outArgs.
get_W_op());
443 implicitModel_->evalModel(appImplicitInArgs, appImplicitOutArgs);
448 #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)
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
Evaluation< VectorBase< Scalar > > get_f() 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 Thyra::ModelEvaluatorBase::InArgs< Scalar > createInArgs() const
virtual void setParameterIndex(int parameterIndex=-1)
Set the parameter index for explicit-only vector.
virtual void initialize()
Initialize after setting member data.
virtual Teuchos::RCP< const Thyra::VectorSpaceBase< Scalar > > get_x_space() const
Get the x-solution space.
ModelEvaluator pair for implicit and explicit (IMEX) evaulations.
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)
basic_FancyOStream & setOutputToRootOnly(const int rootRank)
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()
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.
RCP< LinearOpBase< Scalar > > get_W_op() const
RCP< const VectorBase< Scalar > > get_x() const
virtual void evalModelImpl(const Thyra::ModelEvaluatorBase::InArgs< Scalar > &in, const Thyra::ModelEvaluatorBase::OutArgs< Scalar > &out) const
RCP< const VectorBase< Scalar > > get_p(int l) const