Tempus  Version of the Day
Time Integration
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
HarmonicOscillatorModel_impl.hpp
Go to the documentation of this file.
1 // @HEADER
2 // ****************************************************************************
3 // Tempus: Copyright (2017) Sandia Corporation
4 //
5 // Distributed under BSD 3-clause license (See accompanying file Copyright.txt)
6 // ****************************************************************************
7 // @HEADER
8 
9 #ifndef TEMPUS_TEST_HARMONIC_OSCILLATOR_MODEL_IMPL_HPP
10 #define TEMPUS_TEST_HARMONIC_OSCILLATOR_MODEL_IMPL_HPP
11 
12 #include "Teuchos_StandardParameterEntryValidators.hpp"
13 
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"
21 #include <iostream>
22 
23 namespace Tempus_Test {
24 
25 template <class Scalar>
27  Teuchos::RCP<Teuchos::ParameterList> pList_, const bool use_accel_IC)
28  : out_(Teuchos::VerboseObjectBase::getDefaultOStream())
29 {
30  isInitialized_ = false;
31  setParameterList(pList_);
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";
36  // Divide all coefficients by m_
37  k_ /= m_;
38  f_ /= m_;
39  c_ /= m_;
40  m_ = 1.0;
41  // Set up space and initial guess for solution vector
42  vecLength_ = 1;
43  x_space_ = Thyra::defaultSpmdVectorSpace<Scalar>(vecLength_);
44  x_vec_ = createMember(x_space_);
45  Thyra::put_scalar(0.0, x_vec_.ptr());
46  x_dot_vec_ = createMember(x_space_);
47  Thyra::put_scalar(1.0, x_dot_vec_.ptr());
48  x_dot_dot_vec_ = createMember(x_space_);
49  // The following is the initial condition for the acceleration
50  // Commenting this out to check that IC for acceleration
51  // is computed correctly using displacement and velocity ICs
52  // inside 2nd order steppers.
53  // Thyra::put_scalar(f_-c_, x_dot_dot_vec_.ptr());
54  if (use_accel_IC == true) {
55  Thyra::put_scalar(-2.0, x_dot_dot_vec_.ptr());
56  }
57  else {
58  // Instead of real IC, putting arbitrary, incorrect IC to check correctness
59  // in stepper involving calculation of a IC.
60  Thyra::put_scalar(7.0, x_dot_dot_vec_.ptr());
61  }
62 
63  // Set up responses
64  numResponses_ = 1;
65  g_space_ = Thyra::defaultSpmdVectorSpace<Scalar>(numResponses_);
66 
68 }
69 
70 template <class Scalar>
73 {
74  using Thyra::VectorBase;
75  TEUCHOS_TEST_FOR_EXCEPTION(!isInitialized_, std::logic_error,
76  "Error, setupInOutArgs_ must be called first!\n");
78  double exact_t = t;
79  inArgs.set_t(exact_t);
80  Teuchos::RCP<VectorBase<Scalar> > exact_x = createMember(x_space_);
81  { // scope to delete DetachedVectorView
82  Thyra::DetachedVectorView<Scalar> exact_x_view(*exact_x);
83  if (k_ == 0) {
84  if (c_ == 0)
85  exact_x_view[0] = t * (1.0 + 0.5 * f_ * t);
86  else
87  exact_x_view[0] =
88  (c_ - f_) / (c_ * c_) * (1.0 - exp(-c_ * t)) + f_ * t / c_;
89  }
90  else {
91  exact_x_view[0] = 1.0 / sqrt(k_) * sin(sqrt(k_) * t) +
92  f_ / k_ * (1.0 - cos(sqrt(k_) * t));
93  }
94  }
95  inArgs.set_x(exact_x);
96  Teuchos::RCP<VectorBase<Scalar> > exact_x_dot = createMember(x_space_);
97  { // scope to delete DetachedVectorView
98  Thyra::DetachedVectorView<Scalar> exact_x_dot_view(*exact_x_dot);
99  if (k_ == 0) {
100  if (c_ == 0)
101  exact_x_dot_view[0] = 1.0 + f_ * t;
102  else
103  exact_x_dot_view[0] = (c_ - f_) / c_ * exp(-c_ * t) + f_ / c_;
104  }
105  else {
106  exact_x_dot_view[0] =
107  cos(sqrt(k_) * t) + f_ / sqrt(k_) * sin(sqrt(k_) * t);
108  }
109  }
110  inArgs.set_x_dot(exact_x_dot);
111  Teuchos::RCP<VectorBase<Scalar> > exact_x_dot_dot = createMember(x_space_);
112  { // scope to delete DetachedVectorView
113  Thyra::DetachedVectorView<Scalar> exact_x_dot_dot_view(*exact_x_dot_dot);
114  if (k_ == 0) {
115  if (c_ == 0)
116  exact_x_dot_dot_view[0] = f_;
117  else
118  exact_x_dot_dot_view[0] = (f_ - c_) * exp(-c_ * t);
119  }
120  else {
121  exact_x_dot_dot_view[0] =
122  f_ * cos(sqrt(k_) * t) - sqrt(k_) * sin(sqrt(k_) * t);
123  }
124  }
125  inArgs.set_x_dot_dot(exact_x_dot_dot);
126  return (inArgs);
127 }
128 
129 template <class Scalar>
132 {
133  return x_space_;
134 }
135 
136 template <class Scalar>
139 {
140  return x_space_;
141 }
142 
143 template <class Scalar>
146 {
147  TEUCHOS_TEST_FOR_EXCEPTION(!isInitialized_, std::logic_error,
148  "Error, setupInOutArgs_ must be called first!\n");
149  return nominalValues_;
150 }
151 
152 template <class Scalar>
155 {
157  this->get_W_factory();
158  Teuchos::RCP<Thyra::LinearOpBase<Scalar> > matrix = this->create_W_op();
160  Teuchos::rcp_dynamic_cast<Thyra::MultiVectorBase<Scalar> >(matrix, true);
161  Thyra::DetachedMultiVectorView<Scalar> matrix_view(*matrix_mv);
162  // IKT: is it necessary for W to be non-singular when initialized?
163  matrix_view(0, 0) = 1.0;
165  Thyra::linearOpWithSolve<Scalar>(*W_factory, matrix);
166  return W;
167 }
168 
169 template <class Scalar>
172 {
174  Thyra::createMembers(x_space_, vecLength_);
175  return (matrix);
176 }
177 
178 template <class Scalar>
181 {
183  Thyra::defaultSerialDenseLinearOpWithSolveFactory<Scalar>();
184  return W_factory;
185 }
186 
187 template <class Scalar>
190 {
191  setupInOutArgs_();
192  return inArgs_;
193 }
194 
195 // Private functions overridden from ModelEvaluatorDefaultBase
196 
197 template <class Scalar>
200 {
201  setupInOutArgs_();
202  return outArgs_;
203 }
204 
205 template <class Scalar>
208  const Thyra::ModelEvaluatorBase::OutArgs<Scalar> &outArgs) const
209 {
210  using Teuchos::RCP;
211  using Thyra::VectorBase;
212  TEUCHOS_TEST_FOR_EXCEPTION(!isInitialized_, std::logic_error,
213  "Error, setupInOutArgs_ must be called first!\n");
214 
215  RCP<const VectorBase<Scalar> > x_in = inArgs.get_x();
216  double beta = inArgs.get_beta();
217  if (!x_in.get()) {
219  true, std::logic_error,
220  "\n ERROR: HarmonicOscillatorModel requires x as InArgs.\n");
221  }
223  // IKT, FIXME: check that subDim() is the write routine to get local length of
224  // a Thyra::ConstDetachedVectorView
225  auto myVecLength = x_in_view.subDim();
226 
227  RCP<const VectorBase<Scalar> > x_dot_in = inArgs.get_x_dot();
228  double alpha = inArgs.get_alpha();
229 
230  RCP<const VectorBase<Scalar> > x_dotdot_in = inArgs.get_x_dot_dot();
231  double omega = inArgs.get_W_x_dot_dot_coeff();
232 
233  // Parse OutArgs
234  RCP<VectorBase<Scalar> > f_out = outArgs.get_f();
235  RCP<VectorBase<Scalar> > g_out = outArgs.get_g(0);
236  const RCP<Thyra::LinearOpBase<Scalar> > W_out = outArgs.get_W_op();
237 
238  Scalar neg_sign = 1.0;
239  // Explicit ODE
240  if (inArgs.get_x_dot_dot().is_null()) neg_sign = -1.0;
241 
242  // Populate residual and Jacobian
243  if (f_out != Teuchos::null) {
244  Thyra::DetachedVectorView<Scalar> f_out_view(*f_out);
245  for (int i = 0; i < myVecLength; i++) {
246  f_out_view[i] = f_;
247  }
248  if (x_dotdot_in != Teuchos::null) {
249  Thyra::ConstDetachedVectorView<Scalar> x_dotdot_in_view(*x_dotdot_in);
250  for (int i = 0; i < myVecLength; i++) {
251  f_out_view[i] = x_dotdot_in_view[i] - f_out_view[i];
252  }
253  }
254  if (x_dot_in != Teuchos::null) {
255  Thyra::ConstDetachedVectorView<Scalar> x_dot_in_view(*x_dot_in);
256  for (int i = 0; i < myVecLength; i++) {
257  f_out_view[i] += neg_sign * c_ * x_dot_in_view[i];
258  }
259  }
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];
263  }
264  }
265  }
266 
267  // Note: W = alpha*df/dxdot + beta*df/dx + omega*df/dxdotdot
268  if (W_out != Teuchos::null) {
270  Teuchos::rcp_dynamic_cast<Thyra::MultiVectorBase<Scalar> >(W_out, true);
271  Thyra::DetachedMultiVectorView<Scalar> matrix_view(*matrix);
272  if (omega == 0.0) {
274  true, std::logic_error,
275  "\n ERROR: omega = 0 in HarmonicOscillatorModel!\n");
276  }
277  matrix_view(0, 0) = omega;
278  if (x_dot_in != Teuchos::null) {
279  matrix_view(0, 0) += neg_sign * c_ * alpha;
280  }
281  if (x_in != Teuchos::null) {
282  matrix_view(0, 0) += neg_sign * k_ * beta;
283  }
284  }
285 
286  // Calculated response(s) g
287  // g = mean value of x
288  if (g_out != Teuchos::null) {
289  Thyra::DetachedVectorView<Scalar> g_out_view(*g_out);
290  g_out_view[0] = Thyra::sum(*x_in) / vecLength_;
291  }
292 }
293 
294 template <class Scalar>
297 {
299  true, std::logic_error,
300  "\n Error! HarmonicOscillatorModel::get_p_space() is not supported!\n");
301  return Teuchos::null;
302 }
303 
304 template <class Scalar>
307 {
309  true, std::logic_error,
310  "\n Error! HarmonicOscillatorModel::get_p_names() is not supported!\n");
311  return Teuchos::null;
312 }
313 
314 template <class Scalar>
317 {
319  j != 0, std::logic_error,
320  "\n Error! HarmonicOscillatorModel::get_g_space() only "
321  << " supports 1 parameter vector. Supplied index l = " << j << "\n");
322  return g_space_;
323 }
324 
325 // private
326 
327 template <class Scalar>
329 {
330  if (isInitialized_) return;
331 
332  // Set up InArgs
334  inArgs.setModelEvalDescription(this->description());
335  inArgs.set_Np(0);
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);
343  inArgs_ = inArgs;
344 
345  // Set up OutArgs
347  outArgs.setModelEvalDescription(this->description());
348  outArgs.set_Np_Ng(0, numResponses_);
349 
350  outArgs.setSupports(Thyra::ModelEvaluatorBase::OUT_ARG_f);
351  outArgs.setSupports(Thyra::ModelEvaluatorBase::OUT_ARG_W_op);
352  // outArgs.setSupports(OUT_ARG_W,true);
353  // IKT, what is the following supposed to do??
354  // outArgs.set_W_properties( DerivativeProperties(
355  // DERIV_LINEARITY_UNKNOWN, DERIV_RANK_FULL, true));
356  outArgs_ = outArgs;
357 
358  // Set up nominal values
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_);
364 
365  isInitialized_ = true;
366 }
367 
368 template <class Scalar>
370  Teuchos::RCP<Teuchos::ParameterList> const &paramList)
371 {
372  using Teuchos::get;
374  using Teuchos::RCP;
375  RCP<ParameterList> tmpPL =
376  Teuchos::rcp(new ParameterList("HarmonicOscillatorModel"));
377  if (paramList != Teuchos::null) tmpPL = paramList;
378  tmpPL->validateParametersAndSetDefaults(*this->getValidParameters());
379  this->setMyParamList(tmpPL);
380  RCP<ParameterList> pl = this->getMyNonconstParamList();
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");
385  if (m_ <= 0.0) {
386  TEUCHOS_TEST_FOR_EXCEPTION(true, std::logic_error,
387  "Error: invalid value of Mass coeff m = "
388  << m_ << "! Mass coeff m must be > 0.\n");
389  }
390  if (k_ < 0.0) {
391  TEUCHOS_TEST_FOR_EXCEPTION(true, std::logic_error,
392  "Error: invalid value of x coeff k = "
393  << k_ << "! x coeff k must be >= 0.\n");
394  }
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_
401  << ".\n");
402  }
403 }
404 
405 template <class Scalar>
408 {
410  if (is_null(validPL)) {
411  Teuchos::RCP<Teuchos::ParameterList> pl = Teuchos::parameterList();
412  validPL = pl;
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",
418  &*pl);
419  Teuchos::setDoubleParameter("Mass coeff m", 1.0,
420  "Mass coefficient in model", &*pl);
421  }
422  return validPL;
423 }
424 
425 } // namespace Tempus_Test
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_
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
T * get() 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)
void set_x_dot(const RCP< const VectorBase< Scalar > > &x_dot)
void setParameterList(Teuchos::RCP< Teuchos::ParameterList > const &paramList)
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 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