9 #ifndef TEMPUS_TEST_VANDERPOL_IMEX_EXPLICITMODEL_IMPL_HPP
10 #define TEMPUS_TEST_VANDERPOL_IMEX_EXPLICITMODEL_IMPL_HPP
12 #include "Teuchos_StandardParameterEntryValidators.hpp"
14 #include "Thyra_DefaultSpmdVectorSpace.hpp"
15 #include "Thyra_DetachedVectorView.hpp"
16 #include "Thyra_DetachedMultiVectorView.hpp"
17 #include "Thyra_DefaultSerialDenseLinearOpWithSolveFactory.hpp"
18 #include "Thyra_DefaultMultiVectorLinearOpWithSolve.hpp"
19 #include "Thyra_DefaultLinearOpSource.hpp"
20 #include "Thyra_VectorStdOps.hpp"
21 #include "Thyra_MultiVectorStdOps.hpp"
22 #include "Thyra_DefaultMultiVectorProductVector.hpp"
23 #include "Thyra_DefaultBlockedLinearOp.hpp"
28 namespace Tempus_Test {
30 template<
class Scalar>
33 Teuchos::RCP<Teuchos::ParameterList> pList_,
bool useProductVector)
35 useProductVector_ = useProductVector;
36 isInitialized_ =
false;
42 acceptModelParams_ =
false;
49 if (useProductVector_ ==
false) {
50 x_space_ = Thyra::defaultSpmdVectorSpace<Scalar>(dim_);
51 f_space_ = Thyra::defaultSpmdVectorSpace<Scalar>(dim_);
55 Teuchos::Array<RCP<const Thyra::VectorSpaceBase<Scalar> > > yxSpaceArray;
56 xSpace_ = Thyra::defaultSpmdVectorSpace<Scalar>(1);
57 for (
int i=0; i < dim_; ++i) yxSpaceArray.push_back(xSpace_);
58 x_space_ = Thyra::productVectorSpace<Scalar>(yxSpaceArray);
59 f_space_ = Thyra::productVectorSpace<Scalar>(yxSpaceArray);
63 p_space_ = Thyra::defaultSpmdVectorSpace<Scalar>(np_);
64 g_space_ = Thyra::defaultSpmdVectorSpace<Scalar>(ng_);
65 dxdp_space_ = Thyra::multiVectorProductVectorSpace(x_space_, np_);
67 setParameterList(pList_);
70 template<
class Scalar>
71 Teuchos::RCP<const Thyra::VectorSpaceBase<Scalar> >
79 template<
class Scalar>
80 Teuchos::RCP<const Thyra::VectorSpaceBase<Scalar> >
88 template<
class Scalar>
89 Thyra::ModelEvaluatorBase::InArgs<Scalar>
93 TEUCHOS_TEST_FOR_EXCEPTION( !isInitialized_, std::logic_error,
94 "Error, setupInOutArgs_ must be called first!\n");
95 return nominalValues_;
99 template<
class Scalar>
100 Teuchos::RCP<Thyra::LinearOpBase<Scalar> >
104 Teuchos::RCP<Thyra::LinearOpBase<Scalar> > W_op;
105 if (useProductVector_ ==
false)
106 W_op = Thyra::createMembers(x_space_, dim_);
110 Teuchos::RCP<Thyra::LinearOpBase<Scalar> > A00 =
111 Thyra::createMembers(xSpace_, 1);
112 Teuchos::RCP<Thyra::LinearOpBase<Scalar> > A01 =
113 Thyra::createMembers(xSpace_, 1);
114 Teuchos::RCP<Thyra::LinearOpBase<Scalar> > A10 =
115 Thyra::createMembers(xSpace_, 1);
116 Teuchos::RCP<Thyra::LinearOpBase<Scalar> > A11 =
117 Thyra::createMembers(xSpace_, 1);
118 W_op = Thyra::nonconstBlock2x2(A00, A01, A10, A11);
124 template<
class Scalar>
125 Thyra::ModelEvaluatorBase::InArgs<Scalar>
134 template<
class Scalar>
135 Thyra::ModelEvaluatorBase::OutArgs<Scalar>
144 template<
class Scalar>
149 if (isInitialized_) {
155 Thyra::ModelEvaluatorBase::InArgsSetup<Scalar> inArgs;
156 inArgs.setModelEvalDescription(this->description());
157 inArgs.setSupports( Thyra::ModelEvaluatorBase::IN_ARG_x_dot );
158 inArgs.setSupports( Thyra::ModelEvaluatorBase::IN_ARG_x );
159 inArgs.setSupports( Thyra::ModelEvaluatorBase::IN_ARG_t );
160 inArgs.setSupports( Thyra::ModelEvaluatorBase::IN_ARG_alpha );
161 inArgs.setSupports( Thyra::ModelEvaluatorBase::IN_ARG_beta );
162 if (acceptModelParams_) {
170 Thyra::ModelEvaluatorBase::OutArgsSetup<Scalar> outArgs;
171 outArgs.setModelEvalDescription(this->description());
172 outArgs.setSupports( Thyra::ModelEvaluatorBase::OUT_ARG_f );
173 outArgs.setSupports( Thyra::ModelEvaluatorBase::OUT_ARG_W_op );
174 if (acceptModelParams_) {
175 outArgs.set_Np_Ng(Np_,Ng_);
176 outArgs.setSupports( Thyra::ModelEvaluatorBase::OUT_ARG_DfDp,0,
177 Thyra::ModelEvaluatorBase::DERIV_MV_BY_COL );
183 nominalValues_ = inArgs_;
188 nominalValues_.set_t(t0_ic_);
189 const RCP<Thyra::VectorBase<Scalar> > x_ic = createMember(x_space_);
191 Thyra::DetachedVectorView<Scalar> x_ic_view( *x_ic );
192 x_ic_view[0] = x0_ic_;
193 x_ic_view[1] = x1_ic_;
195 nominalValues_.set_x(x_ic);
197 if (acceptModelParams_) {
198 const RCP<Thyra::VectorBase<Scalar> > p_ic = createMember(p_space_);
200 Thyra::DetachedVectorView<Scalar> p_ic_view( *p_ic );
201 p_ic_view[0] = epsilon_;
203 nominalValues_.set_p(0,p_ic);
205 const RCP<Thyra::VectorBase<Scalar> > x_dot_ic = createMember(x_space_);
207 Thyra::DetachedVectorView<Scalar> x_dot_ic_view( *x_dot_ic );
208 x_dot_ic_view[0] = x1_ic_;
209 x_dot_ic_view[1] = ((1.0-x0_ic_*x0_ic_)*x1_ic_-x0_ic_)/epsilon_;
211 nominalValues_.set_x_dot(x_dot_ic);
214 isInitialized_ =
true;
217 template<
class Scalar>
223 using Teuchos::ParameterList;
224 Teuchos::RCP<ParameterList> tmpPL = Teuchos::rcp(
new ParameterList(
"VanDerPol_IMEX_ExplicitModel"));
225 if (paramList != Teuchos::null) tmpPL = paramList;
226 tmpPL->validateParametersAndSetDefaults(*this->getValidParameters());
227 this->setMyParamList(tmpPL);
228 Teuchos::RCP<ParameterList> pl = this->getMyNonconstParamList();
229 bool acceptModelParams = get<bool>(*pl,
"Accept model parameters");
230 bool haveIC = get<bool>(*pl,
"Provide nominal values");
231 bool useDfDpAsTangent = get<bool>(*pl,
"Use DfDp as Tangent");
232 if ( (acceptModelParams != acceptModelParams_) ||
235 isInitialized_ =
false;
237 acceptModelParams_ = acceptModelParams;
239 useDfDpAsTangent_ = useDfDpAsTangent;
240 epsilon_ = get<Scalar>(*pl,
"Coeff epsilon");
241 x0_ic_ = get<Scalar>(*pl,
"IC x0");
242 x1_ic_ = get<Scalar>(*pl,
"IC x1");
243 t0_ic_ = get<Scalar>(*pl,
"IC t0");
247 template<
class Scalar>
248 Teuchos::RCP<const Teuchos::ParameterList>
252 static Teuchos::RCP<const Teuchos::ParameterList> validPL;
253 if (is_null(validPL)) {
254 Teuchos::RCP<Teuchos::ParameterList> pl = Teuchos::parameterList();
255 pl->set(
"Accept model parameters",
false);
256 pl->set(
"Provide nominal values",
true);
257 pl->set(
"Use DfDp as Tangent",
false);
258 Teuchos::setDoubleParameter(
259 "Coeff epsilon", 1.0e-06,
"Coefficient a in model", &*pl);
260 Teuchos::setDoubleParameter(
261 "IC x0", 2.0,
"Initial Condition for x0", &*pl);
262 Teuchos::setDoubleParameter(
263 "IC x1", 0.0,
"Initial Condition for x1", &*pl);
264 Teuchos::setDoubleParameter(
265 "IC t0", 0.0,
"Initial time t0", &*pl);
271 template<
class Scalar>
272 Teuchos::RCP<const Thyra::VectorSpaceBase<Scalar> >
276 if (!acceptModelParams_) {
277 return Teuchos::null;
279 TEUCHOS_ASSERT_IN_RANGE_UPPER_EXCLUSIVE( l, 0, Np_ );
282 else if (l == 1 || l == 2)
284 return Teuchos::null;
287 template<
class Scalar>
288 Teuchos::RCP<const Teuchos::Array<std::string> >
292 if (!acceptModelParams_) {
293 return Teuchos::null;
295 TEUCHOS_ASSERT_IN_RANGE_UPPER_EXCLUSIVE( l, 0, Np_ );
296 Teuchos::RCP<Teuchos::Array<std::string> > p_strings =
297 Teuchos::rcp(
new Teuchos::Array<std::string>());
299 p_strings->push_back(
"Model Coefficient: epsilon");
301 p_strings->push_back(
"DxDp");
303 p_strings->push_back(
"Dx_dotDp");
307 template<
class Scalar>
308 Teuchos::RCP<const Thyra::VectorSpaceBase<Scalar> >
312 TEUCHOS_ASSERT_IN_RANGE_UPPER_EXCLUSIVE( j, 0, Ng_ );
316 template<
class Scalar>
320 const Thyra::ModelEvaluatorBase::OutArgs<Scalar> &outArgs)
const
322 typedef Thyra::DefaultMultiVectorProductVector<Scalar> DMVPV;
324 using Teuchos::rcp_dynamic_cast;
325 TEUCHOS_TEST_FOR_EXCEPTION( !isInitialized_, std::logic_error,
326 "Error, setupInOutArgs_ must be called first!\n");
328 const RCP<const Thyra::VectorBase<Scalar> > x_in =
329 inArgs.get_x().assert_not_null();
330 Thyra::ConstDetachedVectorView<Scalar> x_in_view( *x_in );
333 Scalar beta = inArgs.get_beta();
334 Scalar eps = epsilon_;
335 RCP<const Thyra::MultiVectorBase<Scalar> > dxdp_in, dxdotdp_in;
336 if (acceptModelParams_) {
337 const RCP<const Thyra::VectorBase<Scalar> > p_in =
339 if (p_in != Teuchos::null) {
340 Thyra::ConstDetachedVectorView<Scalar> p_in_view( *p_in );
343 if (inArgs.get_p(1) != Teuchos::null) {
345 rcp_dynamic_cast<
const DMVPV>(inArgs.get_p(1),
true)->getMultiVector();
347 if (inArgs.get_p(2) != Teuchos::null) {
349 rcp_dynamic_cast<
const DMVPV>(inArgs.get_p(2),
true)->getMultiVector();
353 const RCP<Thyra::VectorBase<Scalar> > f_out = outArgs.get_f();
354 const RCP<Thyra::LinearOpBase<Scalar> > W_out = outArgs.get_W_op();
355 RCP<Thyra::MultiVectorBase<Scalar> > DfDp_out;
356 if (acceptModelParams_) {
357 Thyra::ModelEvaluatorBase::Derivative<Scalar> DfDp = outArgs.get_DfDp(0);
358 DfDp_out = DfDp.getMultiVector();
361 if (inArgs.get_x_dot().is_null()) {
364 if (!is_null(f_out)) {
365 Thyra::DetachedVectorView<Scalar> f_out_view( *f_out );
366 f_out_view[0] = x_in_view[1];
367 f_out_view[1] = -x_in_view[0]/eps;
369 if (!is_null(DfDp_out)) {
370 Thyra::DetachedMultiVectorView<Scalar> DfDp_out_view( *DfDp_out );
371 DfDp_out_view(0,0) = 0.0;
372 DfDp_out_view(1,0) = x_in_view[0]/(eps*eps);
375 if (useDfDpAsTangent_ && !is_null(dxdp_in)) {
376 Thyra::ConstDetachedMultiVectorView<Scalar> dxdp( *dxdp_in );
377 DfDp_out_view(0,0) += dxdp(1,0);
378 DfDp_out_view(1,0) += -dxdp(0,0)/eps;
381 if (!is_null(W_out)) {
382 if (useProductVector_ ==
false) {
383 RCP<Thyra::MultiVectorBase<Scalar> > W =
384 Teuchos::rcp_dynamic_cast<Thyra::MultiVectorBase<Scalar> >(W_out,
true);
385 Thyra::DetachedMultiVectorView<Scalar> W_view( *W );
388 W_view(1,0) = -beta/eps;
393 RCP<Thyra::DefaultBlockedLinearOp<Scalar> > W_block =
394 Teuchos::rcp_dynamic_cast<Thyra::DefaultBlockedLinearOp<Scalar> >(W_out,
true);
395 RCP<Thyra::MultiVectorBase<Scalar> > W00 =
396 Teuchos::rcp_dynamic_cast<Thyra::MultiVectorBase<Scalar> >(W_block->getNonconstBlock(0,0),
true);
397 RCP<Thyra::MultiVectorBase<Scalar> > W01 =
398 Teuchos::rcp_dynamic_cast<Thyra::MultiVectorBase<Scalar> >(W_block->getNonconstBlock(0,1),
true);
399 RCP<Thyra::MultiVectorBase<Scalar> > W10 =
400 Teuchos::rcp_dynamic_cast<Thyra::MultiVectorBase<Scalar> >(W_block->getNonconstBlock(1,0),
true);
401 RCP<Thyra::MultiVectorBase<Scalar> > W11 =
402 Teuchos::rcp_dynamic_cast<Thyra::MultiVectorBase<Scalar> >(W_block->getNonconstBlock(1,1),
true);
403 Thyra::DetachedMultiVectorView<Scalar> W00_view( *W00 );
404 Thyra::DetachedMultiVectorView<Scalar> W01_view( *W01 );
405 Thyra::DetachedMultiVectorView<Scalar> W10_view( *W10 );
406 Thyra::DetachedMultiVectorView<Scalar> W11_view( *W11 );
408 W01_view(0,0) = beta;
409 W10_view(0,0) = -beta/eps;
417 RCP<const Thyra::VectorBase<Scalar> > x_dot_in;
418 x_dot_in = inArgs.get_x_dot().assert_not_null();
419 Scalar alpha = inArgs.get_alpha();
421 if (!is_null(f_out)) {
422 Thyra::DetachedVectorView<Scalar> f_out_view( *f_out );
423 Thyra::ConstDetachedVectorView<Scalar> x_dot_in_view( *x_dot_in );
424 f_out_view[0] = x_dot_in_view[0] - x_in_view[1];
425 f_out_view[1] = x_dot_in_view[1] + x_in_view[0]/eps;
427 if (!is_null(DfDp_out)) {
428 Thyra::DetachedMultiVectorView<Scalar> DfDp_out_view( *DfDp_out );
429 DfDp_out_view(0,0) = 0.0;
430 DfDp_out_view(1,0) = -x_in_view[0]/(eps*eps);
433 if (useDfDpAsTangent_ && !is_null(dxdotdp_in)) {
434 Thyra::ConstDetachedMultiVectorView<Scalar> dxdotdp( *dxdotdp_in );
435 DfDp_out_view(0,0) += dxdotdp(0,0);
436 DfDp_out_view(1,0) += dxdotdp(1,0);
438 if (useDfDpAsTangent_ && !is_null(dxdp_in)) {
439 Thyra::ConstDetachedMultiVectorView<Scalar> dxdp( *dxdp_in );
440 DfDp_out_view(0,0) += -dxdp(1,0);
441 DfDp_out_view(1,0) += dxdp(0,0)/eps;
444 if (!is_null(W_out)) {
445 if (useProductVector_ ==
false) {
446 RCP<Thyra::MultiVectorBase<Scalar> > W =
447 Teuchos::rcp_dynamic_cast<Thyra::MultiVectorBase<Scalar> >(W_out,
true);
448 Thyra::DetachedMultiVectorView<Scalar> W_view( *W );
451 W_view(1,0) = beta/eps;
456 RCP<Thyra::DefaultBlockedLinearOp<Scalar> > W_block =
457 Teuchos::rcp_dynamic_cast<Thyra::DefaultBlockedLinearOp<Scalar> >(W_out,
true);
458 RCP<Thyra::MultiVectorBase<Scalar> > W00 =
459 Teuchos::rcp_dynamic_cast<Thyra::MultiVectorBase<Scalar> >(W_block->getNonconstBlock(0,0),
true);
460 RCP<Thyra::MultiVectorBase<Scalar> > W01 =
461 Teuchos::rcp_dynamic_cast<Thyra::MultiVectorBase<Scalar> >(W_block->getNonconstBlock(0,1),
true);
462 RCP<Thyra::MultiVectorBase<Scalar> > W10 =
463 Teuchos::rcp_dynamic_cast<Thyra::MultiVectorBase<Scalar> >(W_block->getNonconstBlock(1,0),
true);
464 RCP<Thyra::MultiVectorBase<Scalar> > W11 =
465 Teuchos::rcp_dynamic_cast<Thyra::MultiVectorBase<Scalar> >(W_block->getNonconstBlock(1,1),
true);
466 Thyra::DetachedMultiVectorView<Scalar> W00_view( *W00 );
467 Thyra::DetachedMultiVectorView<Scalar> W01_view( *W01 );
468 Thyra::DetachedMultiVectorView<Scalar> W10_view( *W10 );
469 Thyra::DetachedMultiVectorView<Scalar> W11_view( *W11 );
470 W00_view(0,0) = alpha;
471 W01_view(0,0) = -beta;
472 W10_view(0,0) = beta/eps;
473 W11_view(0,0) = alpha;
482 #endif // TEMPUS_TEST_VANDERPOL_IMEX_EXPLICITMODEL_IMPL_HPP
Teuchos::RCP< const Teuchos::ParameterList > getValidParameters() const
void setupInOutArgs_() const
Teuchos::RCP< const Thyra::VectorSpaceBase< Scalar > > get_f_space() const
Teuchos::RCP< const Thyra::VectorSpaceBase< Scalar > > get_x_space() const
Teuchos::RCP< const Thyra::VectorSpaceBase< Scalar > > get_p_space(int l) const
void setParameterList(Teuchos::RCP< Teuchos::ParameterList > const ¶mList)
Teuchos::RCP< Thyra::LinearOpBase< Scalar > > create_W_op() const
Thyra::ModelEvaluatorBase::InArgs< Scalar > getNominalValues() const
Thyra::ModelEvaluatorBase::InArgs< Scalar > createInArgs() const
Thyra::ModelEvaluatorBase::OutArgs< Scalar > createOutArgsImpl() const
VanDerPol_IMEX_ExplicitModel(Teuchos::RCP< Teuchos::ParameterList > pList=Teuchos::null, bool useProductVector=false)
void evalModelImpl(const Thyra::ModelEvaluatorBase::InArgs< Scalar > &inArgs_bar, const Thyra::ModelEvaluatorBase::OutArgs< Scalar > &outArgs_bar) const
Teuchos::RCP< const Teuchos::Array< std::string > > get_p_names(int l) const
Teuchos::RCP< const Thyra::VectorSpaceBase< Scalar > > get_g_space(int j) const