9 #ifndef TEMPUS_TEST_HARMONIC_OSCILLATOR_MODEL_IMPL_HPP
10 #define TEMPUS_TEST_HARMONIC_OSCILLATOR_MODEL_IMPL_HPP
12 #include "Teuchos_StandardParameterEntryValidators.hpp"
14 #include "Thyra_DefaultSpmdVectorSpace.hpp"
15 #include "Thyra_DetachedVectorView.hpp"
16 #include "Thyra_DetachedMultiVectorView.hpp"
17 #include "Thyra_DefaultSerialDenseLinearOpWithSolveFactory.hpp"
18 #include "Thyra_DefaultMultiVectorLinearOpWithSolve.hpp"
19 #include "Thyra_DefaultLinearOpSource.hpp"
20 #include "Thyra_VectorStdOps.hpp"
23 namespace Tempus_Test {
25 template <
class Scalar>
28 : out_(Teuchos::VerboseObjectBase::getDefaultOStream())
32 *
out_ <<
"\n\nDamping coeff c = " <<
c_ <<
"\n";
33 *
out_ <<
"Forcing coeff f = " <<
f_ <<
"\n";
34 *
out_ <<
"x coeff k = " <<
k_ <<
"\n";
35 *
out_ <<
"Mass coeff m = " <<
m_ <<
"\n";
45 Thyra::put_scalar(0.0, x_vec_.ptr());
47 Thyra::put_scalar(1.0, x_dot_vec_.ptr());
54 if (use_accel_IC ==
true) {
55 Thyra::put_scalar(-2.0, x_dot_dot_vec_.ptr());
60 Thyra::put_scalar(7.0, x_dot_dot_vec_.ptr());
70 template <
class Scalar>
76 "Error, setupInOutArgs_ must be called first!\n");
79 inArgs.
set_t(exact_t);
85 exact_x_view[0] = t * (1.0 + 0.5 * f_ * t);
88 (c_ - f_) / (c_ * c_) * (1.0 - exp(-c_ * t)) + f_ * t / c_;
91 exact_x_view[0] = 1.0 / sqrt(k_) * sin(sqrt(k_) * t) +
92 f_ / k_ * (1.0 - cos(sqrt(k_) * t));
95 inArgs.
set_x(exact_x);
101 exact_x_dot_view[0] = 1.0 + f_ * t;
103 exact_x_dot_view[0] = (c_ - f_) / c_ * exp(-c_ * t) + f_ / c_;
106 exact_x_dot_view[0] =
107 cos(sqrt(k_) * t) + f_ / sqrt(k_) * sin(sqrt(k_) * t);
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] =
122 f_ * cos(sqrt(k_) * t) - sqrt(k_) * sin(sqrt(k_) * t);
129 template <
class Scalar>
136 template <
class Scalar>
143 template <
class Scalar>
148 "Error, setupInOutArgs_ must be called first!\n");
149 return nominalValues_;
152 template <
class Scalar>
157 this->get_W_factory();
163 matrix_view(0, 0) = 1.0;
165 Thyra::linearOpWithSolve<Scalar>(*W_factory, matrix);
169 template <
class Scalar>
174 Thyra::createMembers(x_space_, vecLength_);
178 template <
class Scalar>
183 Thyra::defaultSerialDenseLinearOpWithSolveFactory<Scalar>();
187 template <
class Scalar>
197 template <
class Scalar>
205 template <
class Scalar>
213 "Error, setupInOutArgs_ must be called first!\n");
219 true, std::logic_error,
220 "\n ERROR: HarmonicOscillatorModel requires x as InArgs.\n");
225 auto myVecLength = x_in_view.
subDim();
238 Scalar neg_sign = 1.0;
243 if (f_out != Teuchos::null) {
245 for (
int i = 0; i < myVecLength; i++) {
248 if (x_dotdot_in != Teuchos::null) {
250 for (
int i = 0; i < myVecLength; i++) {
251 f_out_view[i] = x_dotdot_in_view[i] - f_out_view[i];
254 if (x_dot_in != Teuchos::null) {
256 for (
int i = 0; i < myVecLength; i++) {
257 f_out_view[i] += neg_sign * c_ * x_dot_in_view[i];
260 if (x_in != Teuchos::null) {
261 for (
int i = 0; i < myVecLength; i++) {
262 f_out_view[i] += neg_sign * k_ * x_in_view[i];
268 if (W_out != Teuchos::null) {
274 true, std::logic_error,
275 "\n ERROR: omega = 0 in HarmonicOscillatorModel!\n");
277 matrix_view(0, 0) = omega;
278 if (x_dot_in != Teuchos::null) {
279 matrix_view(0, 0) += neg_sign * c_ * alpha;
281 if (x_in != Teuchos::null) {
282 matrix_view(0, 0) += neg_sign * k_ * beta;
288 if (g_out != Teuchos::null) {
290 g_out_view[0] = Thyra::sum(*x_in) / vecLength_;
294 template <
class Scalar>
299 true, std::logic_error,
300 "\n Error! HarmonicOscillatorModel::get_p_space() is not supported!\n");
301 return Teuchos::null;
304 template <
class Scalar>
309 true, std::logic_error,
310 "\n Error! HarmonicOscillatorModel::get_p_names() is not supported!\n");
311 return Teuchos::null;
314 template <
class Scalar>
319 j != 0, std::logic_error,
320 "\n Error! HarmonicOscillatorModel::get_g_space() only "
321 <<
" supports 1 parameter vector. Supplied index l = " << j <<
"\n");
327 template <
class Scalar>
330 if (isInitialized_)
return;
336 inArgs.
setSupports(Thyra::ModelEvaluatorBase::IN_ARG_x);
337 inArgs.
setSupports(Thyra::ModelEvaluatorBase::IN_ARG_x_dot);
338 inArgs.
setSupports(Thyra::ModelEvaluatorBase::IN_ARG_x_dot_dot);
339 inArgs.
setSupports(Thyra::ModelEvaluatorBase::IN_ARG_t);
340 inArgs.
setSupports(Thyra::ModelEvaluatorBase::IN_ARG_W_x_dot_dot_coeff);
341 inArgs.
setSupports(Thyra::ModelEvaluatorBase::IN_ARG_alpha);
342 inArgs.
setSupports(Thyra::ModelEvaluatorBase::IN_ARG_beta);
350 outArgs.
setSupports(Thyra::ModelEvaluatorBase::OUT_ARG_f);
351 outArgs.
setSupports(Thyra::ModelEvaluatorBase::OUT_ARG_W_op);
359 nominalValues_ = inArgs_;
360 nominalValues_.set_t(0.0);
361 nominalValues_.set_x(x_vec_);
362 nominalValues_.set_x_dot(x_dot_vec_);
363 nominalValues_.set_x_dot_dot(x_dot_dot_vec_);
365 isInitialized_ =
true;
368 template <
class Scalar>
377 if (paramList != Teuchos::null) tmpPL = paramList;
379 this->setMyParamList(tmpPL);
381 c_ = get<Scalar>(*pl,
"Damping coeff c");
382 f_ = get<Scalar>(*pl,
"Forcing coeff f");
383 k_ = get<Scalar>(*pl,
"x coeff k");
384 m_ = get<Scalar>(*pl,
"Mass coeff m");
387 "Error: invalid value of Mass coeff m = "
388 << m_ <<
"! Mass coeff m must be > 0.\n");
392 "Error: invalid value of x coeff k = "
393 << k_ <<
"! x coeff k must be >= 0.\n");
395 if ((k_ > 0.0) && (c_ != 0.0)) {
397 true, std::logic_error,
398 "Error: HarmonicOscillator model only supports x coeff k > 0 when "
399 "Damping coeff c = 0. You have "
400 <<
"specified x coeff k = " << k_ <<
" and Damping coeff c = " << c_
405 template <
class Scalar>
413 Teuchos::setDoubleParameter(
"Damping coeff c", 0.0,
414 "Damping coefficient in model", &*pl);
415 Teuchos::setDoubleParameter(
"Forcing coeff f", -1.0,
416 "Forcing coefficient in model", &*pl);
417 Teuchos::setDoubleParameter(
"x coeff k", 0.0,
"x coefficient in model",
419 Teuchos::setDoubleParameter(
"Mass coeff m", 1.0,
420 "Mass coefficient in model", &*pl);
426 #endif // TEMPUS_TEST_HARMONIC_OSCILLATOR_MODEL_IMPL_HPP
Thyra::ModelEvaluatorBase::InArgs< Scalar > getNominalValues() const
Thyra::ModelEvaluatorBase::OutArgs< Scalar > createOutArgsImpl() const
bool is_null(const boost::shared_ptr< T > &p)
void set_x_dot_dot(const RCP< const VectorBase< Scalar > > &x_dot_dot)
RCP< const VectorBase< Scalar > > get_x_dot() const
void setModelEvalDescription(const std::string &modelEvalDescription)
Teuchos::RCP< Thyra::VectorBase< Scalar > > x_dot_vec_
void set_Np_Ng(int Np, int Ng)
Teuchos::RCP< const Thyra::VectorSpaceBase< Scalar > > g_space_
Teuchos::RCP< const Thyra::LinearOpWithSolveFactoryBase< Scalar > > get_W_factory() const
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
Teuchos::RCP< const Teuchos::ParameterList > getValidParameters() const
Evaluation< VectorBase< Scalar > > get_g(int j) const
Evaluation< VectorBase< Scalar > > get_f() const
Teuchos::RCP< const Teuchos::Array< std::string > > get_p_names(int l) const
Thyra::ModelEvaluatorBase::InArgs< Scalar > getExactSolution(double t) const
void evalModelImpl(const Thyra::ModelEvaluatorBase::InArgs< Scalar > &inArgs_bar, const Thyra::ModelEvaluatorBase::OutArgs< Scalar > &outArgs_bar) const
Teuchos::RCP< const Thyra::VectorSpaceBase< Scalar > > get_x_space() const
void setSupports(EInArgsMembers arg, bool supports=true)
Teuchos::RCP< const Thyra::VectorSpaceBase< Scalar > > get_g_space(int j) const
void set_x(const RCP< const VectorBase< Scalar > > &x)
Teuchos::RCP< Teuchos::FancyOStream > out_
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
Teuchos::RCP< Thyra::VectorBase< Scalar > > x_dot_dot_vec_
Teuchos_Ordinal subDim() const
Teuchos::RCP< const Thyra::VectorSpaceBase< Scalar > > get_p_space(int l) const
Teuchos::RCP< Thyra::LinearOpWithSolveBase< Scalar > > create_W() const
void validateParametersAndSetDefaults(ParameterList const &validParamList, int const depth=1000)
Scalar get_W_x_dot_dot_coeff() const
void set_x_dot(const RCP< const VectorBase< Scalar > > &x_dot)
void setParameterList(Teuchos::RCP< Teuchos::ParameterList > const ¶mList)
Teuchos::RCP< const Thyra::VectorSpaceBase< Scalar > > x_space_
void setSupports(EOutArgsMembers arg, bool supports=true)
Teuchos::RCP< Thyra::LinearOpBase< Scalar > > create_W_op() const
HarmonicOscillatorModel(Teuchos::RCP< Teuchos::ParameterList > pList=Teuchos::null, const bool use_accel_IC=false)
Teuchos::RCP< Thyra::VectorBase< Scalar > > x_vec_
RCP< const VectorBase< Scalar > > get_x_dot_dot() const
void setupInOutArgs_() const
void setModelEvalDescription(const std::string &modelEvalDescription)
Thyra::ModelEvaluatorBase::InArgs< Scalar > createInArgs() const
RCP< LinearOpBase< Scalar > > get_W_op() const
RCP< const VectorBase< Scalar > > get_x() const
Teuchos::RCP< const Thyra::VectorSpaceBase< Scalar > > get_f_space() const