135 if (isInitialized_) {
141 Thyra::ModelEvaluatorBase::InArgsSetup<Scalar> inArgs;
142 inArgs.setModelEvalDescription(this->description());
143 inArgs.setSupports(Thyra::ModelEvaluatorBase::IN_ARG_x_dot);
144 inArgs.setSupports(Thyra::ModelEvaluatorBase::IN_ARG_x);
145 inArgs.setSupports(Thyra::ModelEvaluatorBase::IN_ARG_t);
146 inArgs.setSupports(Thyra::ModelEvaluatorBase::IN_ARG_alpha);
147 inArgs.setSupports(Thyra::ModelEvaluatorBase::IN_ARG_beta);
148 if (acceptModelParams_) {
156 Thyra::ModelEvaluatorBase::OutArgsSetup<Scalar> outArgs;
157 outArgs.setModelEvalDescription(this->description());
158 outArgs.setSupports(Thyra::ModelEvaluatorBase::OUT_ARG_f);
159 outArgs.setSupports(Thyra::ModelEvaluatorBase::OUT_ARG_W_op);
160 if (acceptModelParams_) {
161 outArgs.set_Np_Ng(Np_, Ng_);
162 outArgs.setSupports(Thyra::ModelEvaluatorBase::OUT_ARG_DfDp, 0,
163 Thyra::ModelEvaluatorBase::DERIV_MV_JACOBIAN_FORM);
169 nominalValues_ = inArgs_;
173 nominalValues_.set_t(t0_ic_);
174 const RCP<Thyra::VectorBase<Scalar> > x_ic = createMember(x_space_);
176 Thyra::DetachedVectorView<Scalar> x_ic_view(*x_ic);
177 x_ic_view[0] = x0_ic_;
178 x_ic_view[1] = x1_ic_;
180 nominalValues_.set_x(x_ic);
182 if (acceptModelParams_) {
183 const RCP<Thyra::VectorBase<Scalar> > p_ic = createMember(p_space_);
185 Thyra::DetachedVectorView<Scalar> p_ic_view(*p_ic);
186 p_ic_view[0] = epsilon_;
188 nominalValues_.set_p(0, p_ic);
190 const RCP<Thyra::VectorBase<Scalar> > x_dot_ic = createMember(x_space_);
192 Thyra::DetachedVectorView<Scalar> x_dot_ic_view(*x_dot_ic);
193 x_dot_ic_view[0] = x1_ic_;
194 x_dot_ic_view[1] = ((1.0 - x0_ic_ * x0_ic_) * x1_ic_ - x0_ic_) / epsilon_;
196 nominalValues_.set_x_dot(x_dot_ic);
199 isInitialized_ =
true;
294 const Thyra::ModelEvaluatorBase::InArgs<Scalar> &inArgs,
295 const Thyra::ModelEvaluatorBase::OutArgs<Scalar> &outArgs)
const
297 typedef Thyra::DefaultMultiVectorProductVector<Scalar> DMVPV;
299 using Teuchos::rcp_dynamic_cast;
300 TEUCHOS_TEST_FOR_EXCEPTION(!isInitialized_, std::logic_error,
301 "Error, setupInOutArgs_ must be called first!\n");
303 const RCP<const Thyra::VectorBase<Scalar> > x_in =
304 inArgs.get_x().assert_not_null();
305 Thyra::ConstDetachedVectorView<Scalar> x_in_view(*x_in);
308 Scalar beta = inArgs.get_beta();
309 Scalar eps = epsilon_;
310 RCP<const Thyra::MultiVectorBase<Scalar> > dxdp_in, dxdotdp_in;
311 if (acceptModelParams_) {
312 const RCP<const Thyra::VectorBase<Scalar> > p_in = inArgs.get_p(0);
313 if (p_in != Teuchos::null) {
314 Thyra::ConstDetachedVectorView<Scalar> p_in_view(*p_in);
317 if (inArgs.get_p(1) != Teuchos::null) {
318 dxdp_in = rcp_dynamic_cast<const DMVPV>(inArgs.get_p(1),
true)
321 if (inArgs.get_p(2) != Teuchos::null) {
322 dxdotdp_in = rcp_dynamic_cast<const DMVPV>(inArgs.get_p(2),
true)
327 const RCP<Thyra::VectorBase<Scalar> > f_out = outArgs.get_f();
328 const RCP<Thyra::LinearOpBase<Scalar> > W_out = outArgs.get_W_op();
329 RCP<Thyra::MultiVectorBase<Scalar> > DfDp_out;
330 if (acceptModelParams_) {
331 Thyra::ModelEvaluatorBase::Derivative<Scalar> DfDp = outArgs.get_DfDp(0);
332 DfDp_out = DfDp.getMultiVector();
335 if (inArgs.get_x_dot().is_null()) {
337 if (!is_null(f_out)) {
338 Thyra::DetachedVectorView<Scalar> f_out_view(*f_out);
339 f_out_view[0] = x_in_view[1];
340 f_out_view[1] = -x_in_view[0] / eps;
342 if (!is_null(DfDp_out)) {
343 Thyra::DetachedMultiVectorView<Scalar> DfDp_out_view(*DfDp_out);
344 DfDp_out_view(0, 0) = 0.0;
345 DfDp_out_view(1, 0) = x_in_view[0] / (eps * eps);
348 if (useDfDpAsTangent_ && !is_null(dxdp_in)) {
349 Thyra::ConstDetachedMultiVectorView<Scalar> dxdp(*dxdp_in);
350 DfDp_out_view(0, 0) += dxdp(1, 0);
351 DfDp_out_view(1, 0) += -dxdp(0, 0) / eps;
354 if (!is_null(W_out)) {
355 if (useProductVector_ ==
false) {
356 RCP<Thyra::MultiVectorBase<Scalar> > W =
357 Teuchos::rcp_dynamic_cast<Thyra::MultiVectorBase<Scalar> >(W_out,
359 Thyra::DetachedMultiVectorView<Scalar> W_view(*W);
362 W_view(1, 0) = -beta / eps;
367 RCP<Thyra::DefaultBlockedLinearOp<Scalar> > W_block =
368 Teuchos::rcp_dynamic_cast<Thyra::DefaultBlockedLinearOp<Scalar> >(
370 RCP<Thyra::MultiVectorBase<Scalar> > W00 =
371 Teuchos::rcp_dynamic_cast<Thyra::MultiVectorBase<Scalar> >(
372 W_block->getNonconstBlock(0, 0),
true);
373 RCP<Thyra::MultiVectorBase<Scalar> > W01 =
374 Teuchos::rcp_dynamic_cast<Thyra::MultiVectorBase<Scalar> >(
375 W_block->getNonconstBlock(0, 1),
true);
376 RCP<Thyra::MultiVectorBase<Scalar> > W10 =
377 Teuchos::rcp_dynamic_cast<Thyra::MultiVectorBase<Scalar> >(
378 W_block->getNonconstBlock(1, 0),
true);
379 RCP<Thyra::MultiVectorBase<Scalar> > W11 =
380 Teuchos::rcp_dynamic_cast<Thyra::MultiVectorBase<Scalar> >(
381 W_block->getNonconstBlock(1, 1),
true);
382 Thyra::DetachedMultiVectorView<Scalar> W00_view(*W00);
383 Thyra::DetachedMultiVectorView<Scalar> W01_view(*W01);
384 Thyra::DetachedMultiVectorView<Scalar> W10_view(*W10);
385 Thyra::DetachedMultiVectorView<Scalar> W11_view(*W11);
386 W00_view(0, 0) = 0.0;
387 W01_view(0, 0) = beta;
388 W10_view(0, 0) = -beta / eps;
389 W11_view(0, 0) = 0.0;
396 RCP<const Thyra::VectorBase<Scalar> > x_dot_in;
397 x_dot_in = inArgs.get_x_dot().assert_not_null();
398 Scalar alpha = inArgs.get_alpha();
400 if (!is_null(f_out)) {
401 Thyra::DetachedVectorView<Scalar> f_out_view(*f_out);
402 Thyra::ConstDetachedVectorView<Scalar> x_dot_in_view(*x_dot_in);
403 f_out_view[0] = x_dot_in_view[0] - x_in_view[1];
404 f_out_view[1] = x_dot_in_view[1] + x_in_view[0] / eps;
406 if (!is_null(DfDp_out)) {
407 Thyra::DetachedMultiVectorView<Scalar> DfDp_out_view(*DfDp_out);
408 DfDp_out_view(0, 0) = 0.0;
409 DfDp_out_view(1, 0) = -x_in_view[0] / (eps * eps);
412 if (useDfDpAsTangent_ && !is_null(dxdotdp_in)) {
413 Thyra::ConstDetachedMultiVectorView<Scalar> dxdotdp(*dxdotdp_in);
414 DfDp_out_view(0, 0) += dxdotdp(0, 0);
415 DfDp_out_view(1, 0) += dxdotdp(1, 0);
417 if (useDfDpAsTangent_ && !is_null(dxdp_in)) {
418 Thyra::ConstDetachedMultiVectorView<Scalar> dxdp(*dxdp_in);
419 DfDp_out_view(0, 0) += -dxdp(1, 0);
420 DfDp_out_view(1, 0) += dxdp(0, 0) / eps;
423 if (!is_null(W_out)) {
424 if (useProductVector_ ==
false) {
425 RCP<Thyra::MultiVectorBase<Scalar> > W =
426 Teuchos::rcp_dynamic_cast<Thyra::MultiVectorBase<Scalar> >(W_out,
428 Thyra::DetachedMultiVectorView<Scalar> W_view(*W);
429 W_view(0, 0) = alpha;
430 W_view(0, 1) = -beta;
431 W_view(1, 0) = beta / eps;
432 W_view(1, 1) = alpha;
436 RCP<Thyra::DefaultBlockedLinearOp<Scalar> > W_block =
437 Teuchos::rcp_dynamic_cast<Thyra::DefaultBlockedLinearOp<Scalar> >(
439 RCP<Thyra::MultiVectorBase<Scalar> > W00 =
440 Teuchos::rcp_dynamic_cast<Thyra::MultiVectorBase<Scalar> >(
441 W_block->getNonconstBlock(0, 0),
true);
442 RCP<Thyra::MultiVectorBase<Scalar> > W01 =
443 Teuchos::rcp_dynamic_cast<Thyra::MultiVectorBase<Scalar> >(
444 W_block->getNonconstBlock(0, 1),
true);
445 RCP<Thyra::MultiVectorBase<Scalar> > W10 =
446 Teuchos::rcp_dynamic_cast<Thyra::MultiVectorBase<Scalar> >(
447 W_block->getNonconstBlock(1, 0),
true);
448 RCP<Thyra::MultiVectorBase<Scalar> > W11 =
449 Teuchos::rcp_dynamic_cast<Thyra::MultiVectorBase<Scalar> >(
450 W_block->getNonconstBlock(1, 1),
true);
451 Thyra::DetachedMultiVectorView<Scalar> W00_view(*W00);
452 Thyra::DetachedMultiVectorView<Scalar> W01_view(*W01);
453 Thyra::DetachedMultiVectorView<Scalar> W10_view(*W10);
454 Thyra::DetachedMultiVectorView<Scalar> W11_view(*W11);
455 W00_view(0, 0) = alpha;
456 W01_view(0, 0) = -beta;
457 W10_view(0, 0) = beta / eps;
458 W11_view(0, 0) = alpha;