9 #ifndef TEMPUS_TEST_VANDERPOL_IMEX_IMPLICITMODEL_IMPL_HPP
10 #define TEMPUS_TEST_VANDERPOL_IMEX_IMPLICITMODEL_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"
27 namespace Tempus_Test {
29 template<
class Scalar>
33 isInitialized_ =
false;
39 acceptModelParams_ =
false;
40 useDfDpAsTangent_ =
false;
48 x_space_ = Thyra::defaultSpmdVectorSpace<Scalar>(dim_);
49 f_space_ = Thyra::defaultSpmdVectorSpace<Scalar>(dim_);
51 p_space_ = Thyra::defaultSpmdVectorSpace<Scalar>(np_);
52 dxdp_space_ = Thyra::multiVectorProductVectorSpace(x_space_, np_);
56 setParameterList(pList_);
59 template<
class Scalar>
60 Teuchos::RCP<const Thyra::VectorSpaceBase<Scalar> >
67 template<
class Scalar>
68 Teuchos::RCP<const Thyra::VectorSpaceBase<Scalar> >
75 template<
class Scalar>
76 Thyra::ModelEvaluatorBase::InArgs<Scalar>
80 TEUCHOS_TEST_FOR_EXCEPTION( !isInitialized_, std::logic_error,
81 "Error, setupInOutArgs_ must be called first!\n");
82 return nominalValues_;
85 template<
class Scalar>
86 Teuchos::RCP<Thyra::LinearOpWithSolveBase<Scalar> >
91 RCP<const Thyra::LinearOpWithSolveFactoryBase<Scalar> > W_factory = this->get_W_factory();
92 RCP<Thyra::LinearOpBase<Scalar> > matrix = this->create_W_op();
99 RCP<Thyra::MultiVectorBase<Scalar> > multivec = Teuchos::rcp_dynamic_cast<Thyra::MultiVectorBase<Scalar> >(matrix,
true);
101 RCP<Thyra::VectorBase<Scalar> > vec = Thyra::createMember(x_space_);
103 Thyra::DetachedVectorView<Scalar> vec_view( *vec );
107 V_V(multivec->col(0).ptr(),*vec);
109 Thyra::DetachedVectorView<Scalar> vec_view( *vec );
113 V_V(multivec->col(1).ptr(),*vec);
116 RCP<Thyra::LinearOpWithSolveBase<Scalar> > W =
117 Thyra::linearOpWithSolve<Scalar>(*W_factory, matrix);
123 template<
class Scalar>
124 Teuchos::RCP<Thyra::LinearOpBase<Scalar> >
128 Teuchos::RCP<Thyra::MultiVectorBase<Scalar> > matrix = Thyra::createMembers(x_space_, dim_);
133 template<
class Scalar>
134 Teuchos::RCP<const Thyra::LinearOpWithSolveFactoryBase<Scalar> >
138 Teuchos::RCP<Thyra::LinearOpWithSolveFactoryBase<Scalar> > W_factory =
139 Thyra::defaultSerialDenseLinearOpWithSolveFactory<Scalar>();
144 template<
class Scalar>
145 Thyra::ModelEvaluatorBase::InArgs<Scalar>
157 template<
class Scalar>
158 Thyra::ModelEvaluatorBase::OutArgs<Scalar>
167 template<
class Scalar>
171 const Thyra::ModelEvaluatorBase::InArgs<Scalar> &inArgs,
172 const Thyra::ModelEvaluatorBase::OutArgs<Scalar> &outArgs
175 typedef Thyra::DefaultMultiVectorProductVector<Scalar> DMVPV;
177 using Teuchos::rcp_dynamic_cast;
178 TEUCHOS_TEST_FOR_EXCEPTION( !isInitialized_, std::logic_error,
179 "Error, setupInOutArgs_ must be called first!\n");
181 const RCP<const Thyra::VectorBase<Scalar> > x_in =
182 inArgs.get_x().assert_not_null();
183 Thyra::ConstDetachedVectorView<Scalar> x_in_view( *x_in );
186 Scalar beta = inArgs.get_beta();
187 Scalar eps = epsilon_;
188 RCP<const Thyra::MultiVectorBase<Scalar> > dxdp_in, dxdotdp_in;
189 if (acceptModelParams_) {
190 const RCP<const Thyra::VectorBase<Scalar> > p_in = inArgs.get_p(0);
191 if (p_in != Teuchos::null) {
192 Thyra::ConstDetachedVectorView<Scalar> p_in_view( *p_in );
195 if (inArgs.get_p(1) != Teuchos::null) {
197 rcp_dynamic_cast<
const DMVPV>(inArgs.get_p(1),
true)->getMultiVector();
199 if (inArgs.get_p(2) != Teuchos::null) {
201 rcp_dynamic_cast<
const DMVPV>(inArgs.get_p(2),
true)->getMultiVector();
205 const RCP<Thyra::VectorBase<Scalar> > f_out = outArgs.get_f();
206 const RCP<Thyra::LinearOpBase<Scalar> > W_out = outArgs.get_W_op();
207 RCP<Thyra::MultiVectorBase<Scalar> > DfDp_out;
208 if (acceptModelParams_) {
209 Thyra::ModelEvaluatorBase::Derivative<Scalar> DfDp = outArgs.get_DfDp(0);
210 DfDp_out = DfDp.getMultiVector();
213 if (inArgs.get_x_dot().is_null()) {
216 if (!is_null(f_out)) {
217 Thyra::DetachedVectorView<Scalar> f_out_view( *f_out );
219 f_out_view[1] = (1.0-x_in_view[0]*x_in_view[0])*x_in_view[1]/eps;
221 if (!is_null(DfDp_out)) {
222 Thyra::DetachedMultiVectorView<Scalar> DfDp_out_view( *DfDp_out );
223 DfDp_out_view(0,0) = 0.0;
224 DfDp_out_view(1,0) = -(1.0-x_in_view[0]*x_in_view[0])*x_in_view[1]/(eps*eps);
227 if (useDfDpAsTangent_ && !is_null(dxdp_in)) {
228 Thyra::ConstDetachedMultiVectorView<Scalar> dxdp( *dxdp_in );
229 DfDp_out_view(1,0) +=
230 -2.0*x_in_view[0]*x_in_view[1]/eps * dxdp(0,0) +
231 (1.0 - x_in_view[0]*x_in_view[0])/eps * dxdp(1,0);
234 if (!is_null(W_out)) {
235 RCP<Thyra::MultiVectorBase<Scalar> > W =
236 Teuchos::rcp_dynamic_cast<Thyra::MultiVectorBase<Scalar> >(W_out,
true);
237 Thyra::DetachedMultiVectorView<Scalar> W_view( *W );
241 -2.0*beta*x_in_view[0]*x_in_view[1]/eps;
243 beta*(1.0 - x_in_view[0]*x_in_view[0])/eps;
249 RCP<const Thyra::VectorBase<Scalar> > x_dot_in;
250 x_dot_in = inArgs.get_x_dot().assert_not_null();
251 Scalar alpha = inArgs.get_alpha();
253 if (!is_null(f_out)) {
254 Thyra::DetachedVectorView<Scalar> f_out_view( *f_out );
255 Thyra::ConstDetachedVectorView<Scalar> x_dot_in_view( *x_dot_in );
256 f_out_view[0] = x_dot_in_view[0];
257 f_out_view[1] = x_dot_in_view[1]
258 - (1.0-x_in_view[0]*x_in_view[0])*x_in_view[1]/eps;;
260 if (!is_null(DfDp_out)) {
261 Thyra::DetachedMultiVectorView<Scalar> DfDp_out_view( *DfDp_out );
262 DfDp_out_view(0,0) = 0.0;
263 DfDp_out_view(1,0) = (1.0-x_in_view[0]*x_in_view[0])*x_in_view[1]/(eps*eps);
266 if (useDfDpAsTangent_ && !is_null(dxdotdp_in)) {
267 Thyra::ConstDetachedMultiVectorView<Scalar> dxdotdp( *dxdotdp_in );
268 DfDp_out_view(0,0) += dxdotdp(0,0);
269 DfDp_out_view(1,0) += dxdotdp(1,0);
271 if (useDfDpAsTangent_ && !is_null(dxdp_in)) {
272 Thyra::ConstDetachedMultiVectorView<Scalar> dxdp( *dxdp_in );
273 DfDp_out_view(1,0) +=
274 2.0*x_in_view[0]*x_in_view[1]/eps * dxdp(0,0) -
275 (1.0 - x_in_view[0]*x_in_view[0])/eps * dxdp(1,0);
278 if (!is_null(W_out)) {
279 RCP<Thyra::MultiVectorBase<Scalar> > W =
280 Teuchos::rcp_dynamic_cast<Thyra::MultiVectorBase<Scalar> >(W_out,
true);
281 Thyra::DetachedMultiVectorView<Scalar> W_view( *W );
285 2.0*beta*x_in_view[0]*x_in_view[1]/eps;
287 - beta*(1.0 - x_in_view[0]*x_in_view[0])/eps;
293 template<
class Scalar>
294 Teuchos::RCP<const Thyra::VectorSpaceBase<Scalar> >
298 if (!acceptModelParams_) {
299 return Teuchos::null;
301 TEUCHOS_ASSERT_IN_RANGE_UPPER_EXCLUSIVE( l, 0, Np_ );
304 else if (l == 1 || l == 2)
306 return Teuchos::null;
309 template<
class Scalar>
310 Teuchos::RCP<const Teuchos::Array<std::string> >
314 if (!acceptModelParams_) {
315 return Teuchos::null;
317 TEUCHOS_ASSERT_IN_RANGE_UPPER_EXCLUSIVE( l, 0, Np_ );
318 Teuchos::RCP<Teuchos::Array<std::string> > p_strings =
319 Teuchos::rcp(
new Teuchos::Array<std::string>());
320 p_strings->push_back(
"Model Coefficient: epsilon");
322 p_strings->push_back(
"Model Coefficient: epsilon");
324 p_strings->push_back(
"DxDp");
326 p_strings->push_back(
"Dx_dotDp");
330 template<
class Scalar>
331 Teuchos::RCP<const Thyra::VectorSpaceBase<Scalar> >
335 TEUCHOS_ASSERT_IN_RANGE_UPPER_EXCLUSIVE( j, 0, Ng_ );
339 template<
class Scalar>
344 if (isInitialized_) {
350 Thyra::ModelEvaluatorBase::InArgsSetup<Scalar> inArgs;
351 inArgs.setModelEvalDescription(this->description());
352 inArgs.setSupports( Thyra::ModelEvaluatorBase::IN_ARG_x_dot );
353 inArgs.setSupports( Thyra::ModelEvaluatorBase::IN_ARG_x );
354 inArgs.setSupports( Thyra::ModelEvaluatorBase::IN_ARG_t );
355 inArgs.setSupports( Thyra::ModelEvaluatorBase::IN_ARG_alpha );
356 inArgs.setSupports( Thyra::ModelEvaluatorBase::IN_ARG_beta );
357 if (acceptModelParams_) {
365 Thyra::ModelEvaluatorBase::OutArgsSetup<Scalar> outArgs;
366 outArgs.setModelEvalDescription(this->description());
367 outArgs.setSupports( Thyra::ModelEvaluatorBase::OUT_ARG_f );
368 outArgs.setSupports( Thyra::ModelEvaluatorBase::OUT_ARG_W_op );
369 if (acceptModelParams_) {
370 outArgs.set_Np_Ng(Np_,Ng_);
371 outArgs.setSupports( Thyra::ModelEvaluatorBase::OUT_ARG_DfDp,0,
372 Thyra::ModelEvaluatorBase::DERIV_MV_BY_COL );
378 nominalValues_ = inArgs_;
381 nominalValues_.set_t(t0_ic_);
382 const Teuchos::RCP<Thyra::VectorBase<Scalar> > x_ic = createMember(x_space_);
384 Thyra::DetachedVectorView<Scalar> x_ic_view( *x_ic );
385 x_ic_view[0] = x0_ic_;
386 x_ic_view[1] = x1_ic_;
388 nominalValues_.set_x(x_ic);
389 if (acceptModelParams_) {
390 const Teuchos::RCP<Thyra::VectorBase<Scalar> > p_ic = createMember(p_space_);
392 Thyra::DetachedVectorView<Scalar> p_ic_view( *p_ic );
393 p_ic_view[0] = epsilon_;
395 nominalValues_.set_p(0,p_ic);
397 const Teuchos::RCP<Thyra::VectorBase<Scalar> > x_dot_ic = createMember(x_space_);
399 Thyra::DetachedVectorView<Scalar> x_dot_ic_view( *x_dot_ic );
400 x_dot_ic_view[0] = x1_ic_;
401 x_dot_ic_view[1] = ((1.0-x0_ic_*x0_ic_)*x1_ic_-x0_ic_)/epsilon_;
403 nominalValues_.set_x_dot(x_dot_ic);
406 isInitialized_ =
true;
410 template<
class Scalar>
416 using Teuchos::ParameterList;
417 Teuchos::RCP<ParameterList> tmpPL =
418 Teuchos::rcp(
new ParameterList(
"VanDerPol_IMEX_ImplicitModel"));
419 if (paramList != Teuchos::null) tmpPL = paramList;
420 tmpPL->validateParametersAndSetDefaults(*this->getValidParameters());
421 this->setMyParamList(tmpPL);
422 Teuchos::RCP<ParameterList> pl = this->getMyNonconstParamList();
423 bool acceptModelParams = get<bool>(*pl,
"Accept model parameters");
424 bool haveIC = get<bool>(*pl,
"Provide nominal values");
425 bool useDfDpAsTangent = get<bool>(*pl,
"Use DfDp as Tangent");
426 if ( (acceptModelParams != acceptModelParams_) ||
429 isInitialized_ =
false;
431 acceptModelParams_ = acceptModelParams;
433 useDfDpAsTangent_ = useDfDpAsTangent;
434 epsilon_ = get<Scalar>(*pl,
"Coeff epsilon");
435 x0_ic_ = get<Scalar>(*pl,
"IC x0");
436 x1_ic_ = get<Scalar>(*pl,
"IC x1");
437 t0_ic_ = get<Scalar>(*pl,
"IC t0");
441 template<
class Scalar>
442 Teuchos::RCP<const Teuchos::ParameterList>
446 static Teuchos::RCP<const Teuchos::ParameterList> validPL;
447 if (is_null(validPL)) {
448 Teuchos::RCP<Teuchos::ParameterList> pl = Teuchos::parameterList();
449 pl->set(
"Accept model parameters",
false);
450 pl->set(
"Provide nominal values",
true);
451 pl->set(
"Use DfDp as Tangent",
false);
452 Teuchos::setDoubleParameter(
453 "Coeff epsilon", 1.0e-06,
"Coefficient a in model", &*pl);
454 Teuchos::setDoubleParameter(
455 "IC x0", 2.0,
"Initial Condition for x0", &*pl);
456 Teuchos::setDoubleParameter(
457 "IC x1", 0.0,
"Initial Condition for x1", &*pl);
458 Teuchos::setDoubleParameter(
459 "IC t0", 0.0,
"Initial time t0", &*pl);
466 #endif // TEMPUS_TEST_VANDERPOL_IMEX_IMPLICITMODEL_IMPL_HPP
Teuchos::RCP< const Thyra::VectorSpaceBase< Scalar > > get_p_space(int l) const
void setupInOutArgs_() const
Teuchos::RCP< Thyra::LinearOpWithSolveBase< Scalar > > create_W() const
Teuchos::RCP< const Teuchos::ParameterList > getValidParameters() const
VanDerPol_IMEX_ImplicitModel(Teuchos::RCP< Teuchos::ParameterList > pList=Teuchos::null)
Teuchos::RCP< const Thyra::LinearOpWithSolveFactoryBase< Scalar > > get_W_factory() const
Thyra::ModelEvaluatorBase::OutArgs< Scalar > createOutArgsImpl() const
Thyra::ModelEvaluatorBase::InArgs< Scalar > getNominalValues() const
Teuchos::RCP< const Thyra::VectorSpaceBase< Scalar > > get_g_space(int j) const
Teuchos::RCP< const Thyra::VectorSpaceBase< Scalar > > get_x_space() const
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
Thyra::ModelEvaluatorBase::InArgs< Scalar > createInArgs() const
void setParameterList(Teuchos::RCP< Teuchos::ParameterList > const ¶mList)
Teuchos::RCP< const Thyra::VectorSpaceBase< Scalar > > get_f_space() const
Teuchos::RCP< Thyra::LinearOpBase< Scalar > > create_W_op() const