170 const Thyra::ModelEvaluatorBase::InArgs<Scalar> &inArgs,
171 const Thyra::ModelEvaluatorBase::OutArgs<Scalar> &outArgs
174 typedef Thyra::DefaultMultiVectorProductVector<Scalar> DMVPV;
176 using Teuchos::rcp_dynamic_cast;
177 TEUCHOS_TEST_FOR_EXCEPTION( !isInitialized_, std::logic_error,
178 "Error, setupInOutArgs_ must be called first!\n");
180 const RCP<const Thyra::VectorBase<Scalar> > x_in =
181 inArgs.get_x().assert_not_null();
182 Thyra::ConstDetachedVectorView<Scalar> x_in_view( *x_in );
183 Scalar x1 = x_in_view[0];
186 Scalar beta = inArgs.get_beta();
187 Scalar eps = epsilon_;
189 const RCP<const Thyra::VectorBase<Scalar> > p_in = inArgs.get_p(0);
190 if (p_in != Teuchos::null) {
191 Thyra::ConstDetachedVectorView<Scalar> p_in_view( *p_in );
195 RCP<const Thyra::MultiVectorBase<Scalar> > dxdp_in, dxdotdp_in, dydp_in;
196 if (inArgs.get_p(1) != Teuchos::null) {
198 rcp_dynamic_cast<const DMVPV>(inArgs.get_p(1),
true)->getMultiVector();
200 if (inArgs.get_p(2) != Teuchos::null) {
202 rcp_dynamic_cast<const DMVPV>(inArgs.get_p(2),
true)->getMultiVector();
204 if (inArgs.get_p(3) != Teuchos::null) {
206 rcp_dynamic_cast<const DMVPV>(inArgs.get_p(3),
true)->getMultiVector();
209 const RCP<const Thyra::VectorBase<Scalar> > y_in =
210 inArgs.get_p(4).assert_not_null();
211 Thyra::ConstDetachedVectorView<Scalar> y_in_view( *y_in );
212 Scalar x0 = y_in_view[0];
214 const RCP<Thyra::VectorBase<Scalar> > f_out = outArgs.get_f();
215 const RCP<Thyra::LinearOpBase<Scalar> > W_out = outArgs.get_W_op();
216 const RCP<Thyra::MultiVectorBase<Scalar> > DfDp_out =
217 outArgs.get_DfDp(0).getMultiVector();
218 const RCP<Thyra::MultiVectorBase<Scalar> > DfDx0_out =
219 outArgs.get_DfDp(4).getMultiVector();
221 if (inArgs.get_x_dot().is_null()) {
224 if (!is_null(f_out)) {
225 Thyra::DetachedVectorView<Scalar> f_out_view( *f_out );
226 f_out_view[0] = (1.0-x0*x0)*x1/eps;
228 if (!is_null(DfDp_out)) {
229 Thyra::DetachedMultiVectorView<Scalar> DfDp_out_view( *DfDp_out );
230 DfDp_out_view(0,0) = -(1.0-x0*x0)*x1/(eps*eps);
233 if (useDfDpAsTangent_ && !is_null(dxdp_in)) {
234 Thyra::ConstDetachedMultiVectorView<Scalar> dxdp( *dxdp_in );
235 DfDp_out_view(0,0) += (1.0 - x0*x0)/eps * dxdp(0,0);
237 if (useDfDpAsTangent_ && !is_null(dydp_in)) {
238 Thyra::ConstDetachedMultiVectorView<Scalar> dydp( *dydp_in );
239 DfDp_out_view(0,0) += -2.0*x0*x1/eps * dydp(0,0);
242 if (!is_null(DfDx0_out)) {
243 Thyra::DetachedMultiVectorView<Scalar> DfDx0_out_view( *DfDx0_out );
244 DfDx0_out_view(0,0) = -2.0*x0*x1/eps;
246 if (!is_null(W_out)) {
247 RCP<Thyra::MultiVectorBase<Scalar> > W =
248 Teuchos::rcp_dynamic_cast<Thyra::MultiVectorBase<Scalar> >(W_out,
true);
249 Thyra::DetachedMultiVectorView<Scalar> W_view( *W );
250 W_view(0,0) = beta*(1.0 - x0*x0)/eps;
256 RCP<const Thyra::VectorBase<Scalar> > x_dot_in;
257 x_dot_in = inArgs.get_x_dot().assert_not_null();
258 if (!is_null(f_out)) {
259 Thyra::DetachedVectorView<Scalar> f_out_view( *f_out );
260 Thyra::ConstDetachedVectorView<Scalar> x_dot_in_view( *x_dot_in );
261 f_out_view[0] = x_dot_in_view[0] - (1.0-x0*x0)*x1/eps;
263 if (!is_null(DfDp_out)) {
264 Thyra::DetachedMultiVectorView<Scalar> DfDp_out_view( *DfDp_out );
265 DfDp_out_view(0,0) = (1.0-x0*x0)*x1/(eps*eps);
268 if (useDfDpAsTangent_ && !is_null(dxdotdp_in)) {
269 Thyra::ConstDetachedMultiVectorView<Scalar> dxdotdp( *dxdotdp_in );
270 DfDp_out_view(0,0) += dxdotdp(0,0);
272 if (useDfDpAsTangent_ && !is_null(dxdp_in)) {
273 Thyra::ConstDetachedMultiVectorView<Scalar> dxdp( *dxdp_in );
274 DfDp_out_view(0,0) += -(1.0 - x0*x0)/eps * dxdp(0,0);
276 if (useDfDpAsTangent_ && !is_null(dydp_in)) {
277 Thyra::ConstDetachedMultiVectorView<Scalar> dydp( *dydp_in );
278 DfDp_out_view(0,0) += 2.0*x0*x1/eps * dydp(0,0);
281 if (!is_null(DfDx0_out)) {
282 Thyra::DetachedMultiVectorView<Scalar> DfDx0_out_view( *DfDx0_out );
283 DfDx0_out_view(0,0) = 2.0*x0*x1/eps;
285 if (!is_null(W_out)) {
286 Scalar alpha = inArgs.get_alpha();
287 RCP<Thyra::MultiVectorBase<Scalar> > W =
288 Teuchos::rcp_dynamic_cast<Thyra::MultiVectorBase<Scalar> >(W_out,
true);
289 Thyra::DetachedMultiVectorView<Scalar> W_view( *W );
290 W_view(0,0) = alpha - beta*(1.0 - x0*x0)/eps;
348 if (isInitialized_)
return;
352 Thyra::ModelEvaluatorBase::InArgsSetup<Scalar> inArgs;
353 inArgs.setModelEvalDescription(this->description());
354 inArgs.setSupports( Thyra::ModelEvaluatorBase::IN_ARG_x_dot );
355 inArgs.setSupports( Thyra::ModelEvaluatorBase::IN_ARG_x );
356 inArgs.setSupports( Thyra::ModelEvaluatorBase::IN_ARG_t );
357 inArgs.setSupports( Thyra::ModelEvaluatorBase::IN_ARG_alpha );
358 inArgs.setSupports( Thyra::ModelEvaluatorBase::IN_ARG_beta );
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 outArgs.set_Np_Ng(Np_,Ng_);
370 outArgs.setSupports( Thyra::ModelEvaluatorBase::OUT_ARG_DfDp,0,
371 Thyra::ModelEvaluatorBase::DERIV_MV_JACOBIAN_FORM );
372 outArgs.setSupports( Thyra::ModelEvaluatorBase::OUT_ARG_DfDp,4,
373 Thyra::ModelEvaluatorBase::DERIV_MV_JACOBIAN_FORM );
378 nominalValues_ = inArgs_;
382 nominalValues_.set_t(t0_ic_);
383 const RCP<Thyra::VectorBase<Scalar> > x_ic = createMember(x_space_);
385 Thyra::DetachedVectorView<Scalar> x_ic_view( *x_ic );
386 x_ic_view[0] = x1_ic_;
388 nominalValues_.set_x(x_ic);
390 const RCP<Thyra::VectorBase<Scalar> > p_ic = createMember(p_space_);
391 const RCP<Thyra::VectorBase<Scalar> > y_ic = createMember(y_space_);
393 Thyra::DetachedVectorView<Scalar> p_ic_view( *p_ic );
394 Thyra::DetachedVectorView<Scalar> y_ic_view( *y_ic );
395 p_ic_view[0] = epsilon_;
396 y_ic_view[0] = x0_ic_;
398 nominalValues_.set_p(0,p_ic);
399 nominalValues_.set_p(4,y_ic);
401 const RCP<Thyra::VectorBase<Scalar> > x_dot_ic = createMember(x_space_);
403 Thyra::DetachedVectorView<Scalar> x_dot_ic_view( *x_dot_ic );
404 x_dot_ic_view[0] = ((1.0-x0_ic_*x0_ic_)*x1_ic_-x0_ic_)/epsilon_;
406 nominalValues_.set_x_dot(x_dot_ic);
409 isInitialized_ =
true;
419 using Teuchos::ParameterList;
420 Teuchos::RCP<ParameterList> tmpPL =
421 Teuchos::rcp(
new ParameterList(
"VanDerPol_IMEXPart_ImplicitModel"));
422 if (paramList != Teuchos::null) tmpPL = paramList;
423 tmpPL->validateParametersAndSetDefaults(*this->getValidParameters());
424 this->setMyParamList(tmpPL);
425 Teuchos::RCP<ParameterList> pl = this->getMyNonconstParamList();
426 bool haveIC = get<bool>(*pl,
"Provide nominal values");
427 bool useDfDpAsTangent = get<bool>(*pl,
"Use DfDp as Tangent");
428 if (haveIC != haveIC_) isInitialized_ =
false;
430 useDfDpAsTangent_ = useDfDpAsTangent;
431 epsilon_ = get<Scalar>(*pl,
"Coeff epsilon");
432 x0_ic_ = get<Scalar>(*pl,
"IC x0");
433 x1_ic_ = get<Scalar>(*pl,
"IC x1");
434 t0_ic_ = get<Scalar>(*pl,
"IC t0");