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);
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;
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_JACOBIAN_FORM );
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;
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");