77 TEUCHOS_TEST_FOR_EXCEPTION( !isInitialized_, std::logic_error,
78 "Error, setupInOutArgs_ must be called first!\n");
79 Thyra::ModelEvaluatorBase::InArgs<Scalar> inArgs = inArgs_;
81 inArgs.set_t(exact_t);
82 Teuchos::RCP<VectorBase<Scalar> > exact_x = createMember(x_space_);
84 Thyra::DetachedVectorView<Scalar> exact_x_view(*exact_x);
87 exact_x_view[0] = t*(1.0+0.5*f_*t);
89 exact_x_view[0] = (c_-f_)/(c_*c_)*(1.0-exp(-c_*t))
93 exact_x_view[0] = 1.0/sqrt(k_)*sin(sqrt(k_)*t) + f_/k_*(1.0-cos(sqrt(k_)*t));
96 inArgs.set_x(exact_x);
97 Teuchos::RCP<VectorBase<Scalar> > exact_x_dot = createMember(x_space_);
99 Thyra::DetachedVectorView<Scalar> exact_x_dot_view(*exact_x_dot);
102 exact_x_dot_view[0] = 1.0+f_*t;
104 exact_x_dot_view[0] = (c_-f_)/c_*exp(-c_*t)+f_/c_;
107 exact_x_dot_view[0] = cos(sqrt(k_)*t) + f_/sqrt(k_)*sin(sqrt(k_)*t);
110 inArgs.set_x_dot(exact_x_dot);
111 Teuchos::RCP<VectorBase<Scalar> > exact_x_dot_dot = createMember(x_space_);
113 Thyra::DetachedVectorView<Scalar> exact_x_dot_dot_view(*exact_x_dot_dot);
116 exact_x_dot_dot_view[0] = f_;
118 exact_x_dot_dot_view[0] = (f_-c_)*exp(-c_*t);
121 exact_x_dot_dot_view[0] = f_*cos(sqrt(k_)*t) - sqrt(k_)*sin(sqrt(k_)*t);
124 inArgs.set_x_dot_dot(exact_x_dot_dot);
162 Teuchos::RCP<const Thyra::LinearOpWithSolveFactoryBase<Scalar> > W_factory = this->get_W_factory();
163 Teuchos::RCP<Thyra::LinearOpBase<Scalar> > matrix = this->create_W_op();
164 Teuchos::RCP<Thyra::MultiVectorBase<Scalar> > matrix_mv = Teuchos::rcp_dynamic_cast<Thyra::MultiVectorBase<Scalar> >(matrix,
true);
165 Thyra::DetachedMultiVectorView<Scalar> matrix_view( *matrix_mv );
167 matrix_view(0,0) = 1.0;
168 Teuchos::RCP<Thyra::LinearOpWithSolveBase<Scalar> > W =
169 Thyra::linearOpWithSolve<Scalar>(*W_factory, matrix );
222 const Thyra::ModelEvaluatorBase::InArgs<Scalar> &inArgs,
223 const Thyra::ModelEvaluatorBase::OutArgs<Scalar> &outArgs
228 TEUCHOS_TEST_FOR_EXCEPTION( !isInitialized_, std::logic_error,
229 "Error, setupInOutArgs_ must be called first!\n");
231 RCP<const VectorBase<Scalar> > x_in = inArgs.get_x();
232 double beta = inArgs.get_beta();
234 TEUCHOS_TEST_FOR_EXCEPTION(
true, std::logic_error,
235 "\n ERROR: HarmonicOscillatorModel requires x as InArgs.\n");
237 Thyra::ConstDetachedVectorView<Scalar> x_in_view( *x_in );
239 auto myVecLength = x_in_view.subDim();
241 RCP<const VectorBase<Scalar> > x_dot_in = inArgs.get_x_dot();
242 double alpha = inArgs.get_alpha();
244 RCP<const VectorBase<Scalar> > x_dotdot_in = inArgs.get_x_dot_dot();
245 double omega = inArgs.get_W_x_dot_dot_coeff();
248 RCP<VectorBase<Scalar> > f_out = outArgs.get_f();
249 RCP<VectorBase<Scalar> > g_out = outArgs.get_g(0);
250 const RCP<Thyra::LinearOpBase<Scalar> > W_out = outArgs.get_W_op();
252 Scalar neg_sign = 1.0;
254 if (inArgs.get_x_dot_dot().is_null()) neg_sign = -1.0;
257 if (f_out != Teuchos::null) {
258 Thyra::DetachedVectorView<Scalar> f_out_view( *f_out );
259 for (
int i=0; i<myVecLength; i++) {
262 if (x_dotdot_in != Teuchos::null) {
263 Thyra::ConstDetachedVectorView<Scalar> x_dotdot_in_view( *x_dotdot_in);
264 for (
int i=0; i<myVecLength; i++) {
265 f_out_view[i] = x_dotdot_in_view[i] - f_out_view[i];
268 if (x_dot_in != Teuchos::null) {
269 Thyra::ConstDetachedVectorView<Scalar> x_dot_in_view( *x_dot_in);
270 for (
int i=0; i<myVecLength; i++) {
271 f_out_view[i] += neg_sign*c_*x_dot_in_view[i];
274 if (x_in != Teuchos::null) {
275 for (
int i=0; i<myVecLength; i++) {
276 f_out_view[i] += neg_sign*k_*x_in_view[i];
282 if (W_out != Teuchos::null) {
283 RCP<Thyra::MultiVectorBase<Scalar> > matrix = Teuchos::rcp_dynamic_cast<Thyra::MultiVectorBase<Scalar> >(W_out,
true);
284 Thyra::DetachedMultiVectorView<Scalar> matrix_view( *matrix );
286 TEUCHOS_TEST_FOR_EXCEPTION(
true, std::logic_error,
287 "\n ERROR: omega = 0 in HarmonicOscillatorModel!\n");
289 matrix_view(0,0) = omega;
290 if (x_dot_in != Teuchos::null) {
291 matrix_view(0,0) += neg_sign*c_*alpha;
293 if (x_in != Teuchos::null) {
294 matrix_view(0,0) += neg_sign*k_*beta;
300 if (g_out != Teuchos::null) {
301 Thyra::DetachedVectorView<Scalar> g_out_view(*g_out);
302 g_out_view[0] = Thyra::sum(*x_in)/vecLength_;
345 if (isInitialized_)
return;
348 Thyra::ModelEvaluatorBase::InArgsSetup<Scalar> inArgs;
349 inArgs.setModelEvalDescription(this->description());
351 inArgs.setSupports(Thyra::ModelEvaluatorBase::IN_ARG_x);
352 inArgs.setSupports(Thyra::ModelEvaluatorBase::IN_ARG_x_dot);
353 inArgs.setSupports(Thyra::ModelEvaluatorBase::IN_ARG_x_dot_dot);
354 inArgs.setSupports(Thyra::ModelEvaluatorBase::IN_ARG_t);
355 inArgs.setSupports(Thyra::ModelEvaluatorBase::IN_ARG_W_x_dot_dot_coeff);
356 inArgs.setSupports(Thyra::ModelEvaluatorBase::IN_ARG_alpha);
357 inArgs.setSupports(Thyra::ModelEvaluatorBase::IN_ARG_beta);
361 Thyra::ModelEvaluatorBase::OutArgsSetup<Scalar> outArgs;
362 outArgs.setModelEvalDescription(this->description());
363 outArgs.set_Np_Ng(0, numResponses_);
365 outArgs.setSupports(Thyra::ModelEvaluatorBase::OUT_ARG_f);
366 outArgs.setSupports(Thyra::ModelEvaluatorBase::OUT_ARG_W_op);
374 nominalValues_ = inArgs_;
375 nominalValues_.set_t(0.0);
376 nominalValues_.set_x(x_vec_);
377 nominalValues_.set_x_dot(x_dot_vec_);
378 nominalValues_.set_x_dot_dot(x_dot_dot_vec_);
380 isInitialized_ =
true;
391 using Teuchos::ParameterList;
392 RCP<ParameterList> tmpPL = Teuchos::rcp(
new ParameterList(
"HarmonicOscillatorModel"));
393 if (paramList != Teuchos::null) tmpPL = paramList;
394 tmpPL->validateParametersAndSetDefaults(*this->getValidParameters());
395 this->setMyParamList(tmpPL);
396 RCP<ParameterList> pl = this->getMyNonconstParamList();
397 c_ = get<Scalar>(*pl,
"Damping coeff c");
398 f_ = get<Scalar>(*pl,
"Forcing coeff f");
399 k_ = get<Scalar>(*pl,
"x coeff k");
400 m_ = get<Scalar>(*pl,
"Mass coeff m");
402 TEUCHOS_TEST_FOR_EXCEPTION(
true, std::logic_error,
403 "Error: invalid value of Mass coeff m = " << m_ <<
"! Mass coeff m must be > 0.\n");
406 TEUCHOS_TEST_FOR_EXCEPTION(
true, std::logic_error,
407 "Error: invalid value of x coeff k = " << k_ <<
"! x coeff k must be >= 0.\n");
409 if ((k_ > 0.0) && (c_ != 0.0)) {
410 TEUCHOS_TEST_FOR_EXCEPTION(
true, std::logic_error,
411 "Error: HarmonicOscillator model only supports x coeff k > 0 when Damping coeff c = 0. You have "
412 <<
"specified x coeff k = " << k_ <<
" and Damping coeff c = " << c_ <<
".\n");