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_JACOBIAN_FORM );
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;
319evalModelImpl(
const Thyra::ModelEvaluatorBase::InArgs<Scalar> &inArgs,
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;