9 #ifndef Tempus_ModelEvaluatorPairPartIMEX_Basic_impl_hpp
10 #define Tempus_ModelEvaluatorPairPartIMEX_Basic_impl_hpp
12 #include "Thyra_ProductVectorBase.hpp"
13 #include "Thyra_ProductVectorSpaceBase.hpp"
18 template <
typename Scalar>
21 : timeDer_(
Teuchos::null), numExplicitOnlyBlocks_(0),
22 parameterIndex_(-1), useImplicitModel_(false)
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)
40 template <
typename Scalar>
44 const Teuchos::RCP<
const Thyra::ModelEvaluator<Scalar> >& explicitModel,
45 const Teuchos::RCP<
const Thyra::ModelEvaluator<Scalar> >& implicitModel,
46 int numExplicitOnlyBlocks,
int parameterIndex)
48 timeDer_ = Teuchos::null;
49 numExplicitOnlyBlocks_ = numExplicitOnlyBlocks;
50 parameterIndex_ = parameterIndex;
51 useImplicitModel_ =
false;
52 setExplicitModel(explicitModel);
53 setImplicitModel(implicitModel);
58 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();
74 TEUCHOS_TEST_FOR_EXCEPTION( !(getIMEXVector(z)->space()->isCompatible(
75 *(implicitModel_->get_x_space()))),
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");
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");
90 int numBlocks = zPVector->productSpace()->numBlocks();
92 TEUCHOS_TEST_FOR_EXCEPTION( !(0 <= numExplicitOnlyBlocks_ and
93 numExplicitOnlyBlocks_ < numBlocks),
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");
100 template <
typename Scalar>
104 const Teuchos::RCP<
const Thyra::ModelEvaluator<Scalar> > & me)
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");
112 template <
typename Scalar>
113 Teuchos::RCP<const Thyra::ModelEvaluator<Scalar> >
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"
124 template <
typename Scalar>
125 Teuchos::RCP<const Thyra::VectorSpaceBase<Scalar> >
129 if (useImplicitModel_ ==
true)
return implicitModel_->get_x_space();
131 return explicitModel_->get_x_space();
134 template <
typename Scalar>
135 Teuchos::RCP<const Thyra::VectorSpaceBase<Scalar> >
139 if (useImplicitModel_ ==
true)
return implicitModel_->get_g_space(i);
141 return explicitModel_->get_g_space(i);
144 template <
typename Scalar>
145 Teuchos::RCP<const Thyra::VectorSpaceBase<Scalar> >
149 if (useImplicitModel_ ==
true)
return implicitModel_->get_p_space(i);
151 return explicitModel_->get_p_space(i);
154 template <
typename Scalar>
158 const Teuchos::RCP<
const Thyra::ModelEvaluator<Scalar> > & model )
160 implicitModel_ = model;
163 template <
typename Scalar>
164 Teuchos::RCP<Thyra::VectorBase<Scalar> >
169 using Teuchos::rcp_dynamic_cast;
171 if(full == Teuchos::null)
172 return Teuchos::null;
174 if(numExplicitOnlyBlocks_==0)
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();
186 if(numBlocks==numExplicitOnlyBlocks_+1)
187 return blk_full->getNonconstVectorBlock(numExplicitOnlyBlocks_);
189 TEUCHOS_ASSERT(
false);
190 return Teuchos::null;
193 template <
typename Scalar>
194 Teuchos::RCP<const Thyra::VectorBase<Scalar> >
196 getIMEXVector(
const Teuchos::RCP<
const Thyra::VectorBase<Scalar> > & full)
const
199 using Teuchos::rcp_dynamic_cast;
201 if(full == Teuchos::null)
202 return Teuchos::null;
204 if(numExplicitOnlyBlocks_==0)
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();
216 if(numBlocks==numExplicitOnlyBlocks_+1)
217 return blk_full->getVectorBlock(numExplicitOnlyBlocks_);
219 TEUCHOS_ASSERT(
false);
220 return Teuchos::null;
223 template <
typename Scalar>
224 Teuchos::RCP<Thyra::VectorBase<Scalar> >
227 const Teuchos::RCP<Thyra::VectorBase<Scalar> > & full)
const
230 using Teuchos::rcp_dynamic_cast;
232 if(numExplicitOnlyBlocks_ == 0 || full == Teuchos::null)
233 return Teuchos::null;
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");
243 if(numExplicitOnlyBlocks_==1)
244 return blk_full->getNonconstVectorBlock(0);
246 TEUCHOS_ASSERT(
false);
247 return Teuchos::null;
251 template <
typename Scalar>
252 Teuchos::RCP<const Thyra::VectorBase<Scalar> >
255 const Teuchos::RCP<
const Thyra::VectorBase<Scalar> > & full)
const
258 using Teuchos::rcp_dynamic_cast;
260 if(numExplicitOnlyBlocks_ == 0 || full == Teuchos::null)
261 return Teuchos::null;
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");
271 if(numExplicitOnlyBlocks_==1)
272 return blk_full->getVectorBlock(0);
274 TEUCHOS_ASSERT(
false);
275 return Teuchos::null;
279 template <
typename Scalar>
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;
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;
305 if(parameterIndex >= 0) {
306 parameterIndex_ = parameterIndex;
307 }
else if (parameterIndex_ < 0) {
309 for(
int i=0; i<implicitModel_->Np(); i++) {
310 if((*implicitModel_->get_p_names(i))[0] ==
"EXPLICIT_ONLY_VECTOR") {
321 template <
typename Scalar>
322 Teuchos::RCP<const Thyra::VectorSpaceBase<Scalar> >
326 if (useImplicitModel_ ==
true)
return implicitModel_->get_f_space();
328 return explicitModel_->get_f_space();
332 template <
typename Scalar>
333 Thyra::ModelEvaluatorBase::InArgs<Scalar>
337 typedef Thyra::ModelEvaluatorBase MEB;
338 MEB::InArgsSetup<Scalar> inArgs = this->createInArgs();
343 template <
typename Scalar>
344 Thyra::ModelEvaluatorBase::InArgs<Scalar>
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());
359 MEB::InArgsSetup<Scalar> inArgs(explicitInArgs);
360 inArgs.setModelEvalDescription(this->description());
365 template <
typename Scalar>
366 Thyra::ModelEvaluatorBase::OutArgs<Scalar>
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());
381 MEB::OutArgsSetup<Scalar> outArgs(explicitOutArgs);
382 outArgs.setModelEvalDescription(this->description());
383 outArgs.set_Np_Ng(np,explicitOutArgs.Ng());
387 template <
typename Scalar>
390 const Thyra::ModelEvaluatorBase::OutArgs<Scalar> & outArgs)
const
392 typedef Thyra::ModelEvaluatorBase MEB;
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);
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) {
406 if ((inArgs.get_p(i) != Teuchos::null) and (i != parameterIndex_))
407 appImplicitInArgs.set_p(i, inArgs.get_p(i));
410 appImplicitOutArgs.set_f(outArgs.get_f());
411 appImplicitOutArgs.set_W_op(outArgs.get_W_op());
413 implicitModel_->evalModel(appImplicitInArgs,appImplicitOutArgs);
418 #endif // Tempus_ModelEvaluatorPairPartIMEX_Basic_impl_hpp