Tempus  Version of the Day
Time Integration
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
VanDerPolModel_impl.hpp
Go to the documentation of this file.
1 //@HEADER
2 // *****************************************************************************
3 // Tempus: Time Integration and Sensitivity Analysis Package
4 //
5 // Copyright 2017 NTESS and the Tempus contributors.
6 // SPDX-License-Identifier: BSD-3-Clause
7 // *****************************************************************************
8 //@HEADER
9 
10 #ifndef TEMPUS_TEST_VANDERPOL_MODEL_IMPL_HPP
11 #define TEMPUS_TEST_VANDERPOL_MODEL_IMPL_HPP
12 
13 #include "Teuchos_StandardParameterEntryValidators.hpp"
14 
15 #include "Thyra_DefaultSpmdVectorSpace.hpp"
16 #include "Thyra_DetachedVectorView.hpp"
17 #include "Thyra_DetachedMultiVectorView.hpp"
18 #include "Thyra_DefaultSerialDenseLinearOpWithSolveFactory.hpp"
19 #include "Thyra_DefaultMultiVectorLinearOpWithSolve.hpp"
20 #include "Thyra_DefaultLinearOpSource.hpp"
21 #include "Thyra_VectorStdOps.hpp"
22 
23 #include <iostream>
24 
25 namespace Tempus_Test {
26 
27 template <class Scalar>
30 {
31  isInitialized_ = false;
32  dim_ = 2;
33  Np_ = 1; // Number of parameter vectors (1)
34  np_ = 1; // Number of parameters in this vector (1)
35  Ng_ = 0; // Number of observation functions (0)
36  ng_ = 0; // Number of elements in this observation function (0)
37  acceptModelParams_ = false;
38  haveIC_ = true;
39  epsilon_ = 1.0e-06;
40  x0_ic_ = 2.0;
41  x1_ic_ = 0.0;
42  t0_ic_ = 0.0;
43 
44  // Create x_space and f_space
45  x_space_ = Thyra::defaultSpmdVectorSpace<Scalar>(dim_);
46  f_space_ = Thyra::defaultSpmdVectorSpace<Scalar>(dim_);
47  // Create p_space and g_space
48  p_space_ = Thyra::defaultSpmdVectorSpace<Scalar>(np_);
49  g_space_ = Thyra::defaultSpmdVectorSpace<Scalar>(ng_);
50 
51  setParameterList(pList_);
52 }
53 
54 template <class Scalar>
57 {
59  true, std::logic_error,
60  "Error - No exact solution for van der Pol problem!\n");
61  return (inArgs_);
62 }
63 
64 template <class Scalar>
67 {
69  !isInitialized_, std::logic_error,
70  "Error - No exact sensitivities for van der Pol problem!\n");
71  return (inArgs_);
72 }
73 
74 template <class Scalar>
77 {
78  return x_space_;
79 }
80 
81 template <class Scalar>
84 {
85  return f_space_;
86 }
87 
88 template <class Scalar>
91 {
92  TEUCHOS_TEST_FOR_EXCEPTION(!isInitialized_, std::logic_error,
93  "Error, setupInOutArgs_ must be called first!\n");
94  return nominalValues_;
95 }
96 
97 template <class Scalar>
100 {
101  using Teuchos::RCP;
103  this->get_W_factory();
104  RCP<Thyra::LinearOpBase<Scalar> > matrix = this->create_W_op();
105  {
106  // 01/20/09 tscoffe: This is a total hack to provide a full rank matrix to
107  // linearOpWithSolve because it ends up factoring the matrix during
108  // initialization, which it really shouldn't do, or I'm doing something
109  // wrong here. The net effect is that I get exceptions thrown in
110  // optimized mode due to the matrix being rank deficient unless I do this.
112  Teuchos::rcp_dynamic_cast<Thyra::MultiVectorBase<Scalar> >(matrix,
113  true);
114  {
115  RCP<Thyra::VectorBase<Scalar> > vec = Thyra::createMember(x_space_);
116  {
117  Thyra::DetachedVectorView<Scalar> vec_view(*vec);
118  vec_view[0] = 0.0;
119  vec_view[1] = 1.0;
120  }
121  V_V(multivec->col(0).ptr(), *vec);
122  {
123  Thyra::DetachedVectorView<Scalar> vec_view(*vec);
124  vec_view[0] = 1.0;
125  vec_view[1] = 0.0;
126  }
127  V_V(multivec->col(1).ptr(), *vec);
128  }
129  }
131  Thyra::linearOpWithSolve<Scalar>(*W_factory, matrix);
132 
133  return W;
134 }
135 
136 template <class Scalar>
138  const
139 {
141  Thyra::createMembers(x_space_, dim_);
142  return (matrix);
143 }
144 
145 template <class Scalar>
148 {
150  Thyra::defaultSerialDenseLinearOpWithSolveFactory<Scalar>();
151  return W_factory;
152 }
153 
154 template <class Scalar>
156  const
157 {
158  setupInOutArgs_();
159  return inArgs_;
160 }
161 
162 // Private functions overridden from ModelEvaluatorDefaultBase
163 
164 template <class Scalar>
167 {
168  setupInOutArgs_();
169  return outArgs_;
170 }
171 
172 template <class Scalar>
175  const Thyra::ModelEvaluatorBase::OutArgs<Scalar> &outArgs) const
176 {
177  using Teuchos::RCP;
178  TEUCHOS_TEST_FOR_EXCEPTION(!isInitialized_, std::logic_error,
179  "Error, setupInOutArgs_ must be called first!\n");
180 
182  inArgs.get_x().assert_not_null();
184 
185  // double t = inArgs.get_t();
186  Scalar epsilon = epsilon_;
187  if (acceptModelParams_) {
189  inArgs.get_p(0).assert_not_null();
191  epsilon = p_in_view[0];
192  }
193 
194  Scalar beta = inArgs.get_beta();
195 
196  const RCP<Thyra::VectorBase<Scalar> > f_out = outArgs.get_f();
197  const RCP<Thyra::LinearOpBase<Scalar> > W_out = outArgs.get_W_op();
199  if (acceptModelParams_) {
201  DfDp_out = DfDp.getMultiVector();
202  }
203 
204  if (inArgs.get_x_dot().is_null()) {
205  // Evaluate the Explicit ODE f(x,t) [= xdot]
206  if (!is_null(f_out)) {
207  Thyra::DetachedVectorView<Scalar> f_out_view(*f_out);
208  f_out_view[0] = x_in_view[1];
209  f_out_view[1] =
210  ((1.0 - x_in_view[0] * x_in_view[0]) * x_in_view[1] - x_in_view[0]) /
211  epsilon;
212  }
213  if (!is_null(W_out)) {
215  Teuchos::rcp_dynamic_cast<Thyra::MultiVectorBase<Scalar> >(W_out,
216  true);
217  Thyra::DetachedMultiVectorView<Scalar> matrix_view(*matrix);
218  matrix_view(0, 0) = 0.0; // d(f0)/d(x0_n)
219  matrix_view(0, 1) = +beta; // d(f0)/d(x1_n)
220  matrix_view(1, 0) = -beta * (2.0 * x_in_view[0] * x_in_view[1] + 1.0) /
221  epsilon; // d(f1)/d(x0_n)
222  matrix_view(1, 1) = beta * (1.0 - x_in_view[0] * x_in_view[0]) /
223  epsilon; // d(f1)/d(x1_n)
224  // Note: alpha = d(xdot)/d(x_n) and beta = d(x)/d(x_n)
225  }
226  if (!is_null(DfDp_out)) {
227  Thyra::DetachedMultiVectorView<Scalar> DfDp_out_view(*DfDp_out);
228  DfDp_out_view(0, 0) = 0.0;
229  DfDp_out_view(1, 0) =
230  -((1.0 - x_in_view[0] * x_in_view[0]) * x_in_view[1] - x_in_view[0]) /
231  (epsilon * epsilon);
232  }
233  }
234  else {
235  // Evaluate the implicit ODE f(xdot, x, t) [= 0]
237  x_dot_in = inArgs.get_x_dot().assert_not_null();
238  Scalar alpha = inArgs.get_alpha();
239  if (!is_null(f_out)) {
240  Thyra::DetachedVectorView<Scalar> f_out_view(*f_out);
241  Thyra::ConstDetachedVectorView<Scalar> x_dot_in_view(*x_dot_in);
242  f_out_view[0] = x_dot_in_view[0] - x_in_view[1];
243  f_out_view[1] =
244  x_dot_in_view[1] -
245  ((1.0 - x_in_view[0] * x_in_view[0]) * x_in_view[1] - x_in_view[0]) /
246  epsilon;
247  }
248  if (!is_null(W_out)) {
250  Teuchos::rcp_dynamic_cast<Thyra::MultiVectorBase<Scalar> >(W_out,
251  true);
252  Thyra::DetachedMultiVectorView<Scalar> matrix_view(*matrix);
253  matrix_view(0, 0) = alpha; // d(f0)/d(x0_n)
254  matrix_view(0, 1) = -beta; // d(f0)/d(x1_n)
255  matrix_view(1, 0) = beta * (2.0 * x_in_view[0] * x_in_view[1] + 1.0) /
256  epsilon; // d(f1)/d(x0_n)
257  matrix_view(1, 1) = alpha - beta * (1.0 - x_in_view[0] * x_in_view[0]) /
258  epsilon; // d(f1)/d(x1_n)
259  // Note: alpha = d(xdot)/d(x_n) and beta = d(x)/d(x_n)
260  }
261  if (!is_null(DfDp_out)) {
262  Thyra::DetachedMultiVectorView<Scalar> DfDp_out_view(*DfDp_out);
263  DfDp_out_view(0, 0) = 0.0;
264  DfDp_out_view(1, 0) =
265  ((1.0 - x_in_view[0] * x_in_view[0]) * x_in_view[1] - x_in_view[0]) /
266  (epsilon * epsilon);
267  }
268  }
269 }
270 
271 template <class Scalar>
274 {
275  if (!acceptModelParams_) {
276  return Teuchos::null;
277  }
279  return p_space_;
280 }
281 
282 template <class Scalar>
285 {
286  if (!acceptModelParams_) {
287  return Teuchos::null;
288  }
292  p_strings->push_back("Model Coefficient: epsilon");
293  return p_strings;
294 }
295 
296 template <class Scalar>
299 {
301  return g_space_;
302 }
303 
304 // private
305 
306 template <class Scalar>
308 {
309  if (isInitialized_) {
310  return;
311  }
312 
313  {
314  // Set up prototypical InArgs
316  inArgs.setModelEvalDescription(this->description());
317  inArgs.setSupports(Thyra::ModelEvaluatorBase::IN_ARG_t);
318  inArgs.setSupports(Thyra::ModelEvaluatorBase::IN_ARG_x);
319  inArgs.setSupports(Thyra::ModelEvaluatorBase::IN_ARG_beta);
320  inArgs.setSupports(Thyra::ModelEvaluatorBase::IN_ARG_x_dot);
321  inArgs.setSupports(Thyra::ModelEvaluatorBase::IN_ARG_alpha);
322  if (acceptModelParams_) {
323  inArgs.set_Np(Np_);
324  }
325  inArgs_ = inArgs;
326  }
327 
328  {
329  // Set up prototypical OutArgs
331  outArgs.setModelEvalDescription(this->description());
332  outArgs.setSupports(Thyra::ModelEvaluatorBase::OUT_ARG_f);
333  outArgs.setSupports(Thyra::ModelEvaluatorBase::OUT_ARG_W_op);
334  if (acceptModelParams_) {
335  outArgs.set_Np_Ng(Np_, Ng_);
336  outArgs.setSupports(Thyra::ModelEvaluatorBase::OUT_ARG_DfDp, 0,
337  Thyra::ModelEvaluatorBase::DERIV_MV_JACOBIAN_FORM);
338  }
339  outArgs_ = outArgs;
340  }
341 
342  // Set up nominal values
343  nominalValues_ = inArgs_;
344  if (haveIC_) {
345  using Teuchos::RCP;
346  nominalValues_.set_t(t0_ic_);
347  const RCP<Thyra::VectorBase<Scalar> > x_ic = createMember(x_space_);
348  { // scope to delete DetachedVectorView
349  Thyra::DetachedVectorView<Scalar> x_ic_view(*x_ic);
350  x_ic_view[0] = x0_ic_;
351  x_ic_view[1] = x1_ic_;
352  }
353  nominalValues_.set_x(x_ic);
354  if (acceptModelParams_) {
355  const RCP<Thyra::VectorBase<Scalar> > p_ic = createMember(p_space_);
356  {
357  Thyra::DetachedVectorView<Scalar> p_ic_view(*p_ic);
358  p_ic_view[0] = epsilon_;
359  }
360  nominalValues_.set_p(0, p_ic);
361  }
362  const RCP<Thyra::VectorBase<Scalar> > x_dot_ic = createMember(x_space_);
363  { // scope to delete DetachedVectorView
364  Thyra::DetachedVectorView<Scalar> x_dot_ic_view(*x_dot_ic);
365  x_dot_ic_view[0] = x1_ic_;
366  x_dot_ic_view[1] = ((1.0 - x0_ic_ * x0_ic_) * x1_ic_ - x0_ic_) / epsilon_;
367  }
368  nominalValues_.set_x_dot(x_dot_ic);
369  }
370 
371  isInitialized_ = true;
372 }
373 
374 template <class Scalar>
376  Teuchos::RCP<Teuchos::ParameterList> const &paramList)
377 {
378  using Teuchos::get;
381  Teuchos::rcp(new ParameterList("VanDerPolModel"));
382  if (paramList != Teuchos::null) tmpPL = paramList;
383  tmpPL->validateParametersAndSetDefaults(*this->getValidParameters());
384  this->setMyParamList(tmpPL);
385  Teuchos::RCP<ParameterList> pl = this->getMyNonconstParamList();
386  bool acceptModelParams = get<bool>(*pl, "Accept model parameters");
387  bool haveIC = get<bool>(*pl, "Provide nominal values");
388  if ((acceptModelParams != acceptModelParams_) || (haveIC != haveIC_)) {
389  isInitialized_ = false;
390  }
391  acceptModelParams_ = acceptModelParams;
392  haveIC_ = haveIC;
393  epsilon_ = get<Scalar>(*pl, "Coeff epsilon");
394  x0_ic_ = get<Scalar>(*pl, "IC x0");
395  x1_ic_ = get<Scalar>(*pl, "IC x1");
396  t0_ic_ = get<Scalar>(*pl, "IC t0");
397  setupInOutArgs_();
398 }
399 
400 template <class Scalar>
403 {
405  if (is_null(validPL)) {
406  Teuchos::RCP<Teuchos::ParameterList> pl = Teuchos::parameterList();
407  pl->set("Accept model parameters", false);
408  pl->set("Provide nominal values", true);
409  Teuchos::setDoubleParameter("Coeff epsilon", 1.0e-06,
410  "Coefficient a in model", &*pl);
411  Teuchos::setDoubleParameter("IC x0", 2.0, "Initial Condition for x0", &*pl);
412  Teuchos::setDoubleParameter("IC x1", 0.0, "Initial Condition for x1", &*pl);
413  Teuchos::setDoubleParameter("IC t0", 0.0, "Initial time t0", &*pl);
414  Teuchos::setIntParameter("Number of Time Step Sizes", 1,
415  "Number time step sizes for convergence study",
416  &*pl);
417  validPL = pl;
418  }
419  return validPL;
420 }
421 
422 } // namespace Tempus_Test
423 #endif // TEMPUS_TEST_VANDERPOL_MODEL_IMPL_HPP
Teuchos::RCP< const Teuchos::ParameterList > getValidParameters() const
bool is_null(const boost::shared_ptr< T > &p)
Derivative< Scalar > get_DfDp(int l) const
RCP< const VectorBase< Scalar > > get_x_dot() const
void setModelEvalDescription(const std::string &modelEvalDescription)
ParameterList & set(std::string const &name, T const &value, std::string const &docString="", RCP< const ParameterEntryValidator > const &validator=null)
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
Evaluation< VectorBase< Scalar > > get_f() const
Thyra::ModelEvaluatorBase::InArgs< Scalar > getNominalValues() const
Thyra::ModelEvaluatorBase::InArgs< Scalar > getExactSolution(double t) const
Thyra::ModelEvaluatorBase::InArgs< Scalar > createInArgs() const
Teuchos::RCP< const Thyra::LinearOpWithSolveFactoryBase< Scalar > > get_W_factory() const
VanDerPolModel(Teuchos::RCP< Teuchos::ParameterList > pList=Teuchos::null)
#define TEUCHOS_ASSERT_IN_RANGE_UPPER_EXCLUSIVE(index, lower_inclusive, upper_exclusive)
void setSupports(EInArgsMembers arg, bool supports=true)
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
Teuchos::RCP< Thyra::LinearOpBase< Scalar > > create_W_op() const
Teuchos::RCP< Thyra::LinearOpWithSolveBase< Scalar > > create_W() const
Ptr< T > ptr() const
void validateParametersAndSetDefaults(ParameterList const &validParamList, int const depth=1000)
void evalModelImpl(const Thyra::ModelEvaluatorBase::InArgs< Scalar > &inArgs_bar, const Thyra::ModelEvaluatorBase::OutArgs< Scalar > &outArgs_bar) const
Teuchos::RCP< const Teuchos::Array< std::string > > get_p_names(int l) const
Teuchos::RCP< const Thyra::VectorSpaceBase< Scalar > > get_g_space(int j) const
void setSupports(EOutArgsMembers arg, bool supports=true)
void setParameterList(Teuchos::RCP< Teuchos::ParameterList > const &paramList)
Teuchos::RCP< const Thyra::VectorSpaceBase< Scalar > > get_x_space() const
Thyra::ModelEvaluatorBase::InArgs< Scalar > getExactSensSolution(int j, double t) const
Teuchos::RCP< const Thyra::VectorSpaceBase< Scalar > > get_p_space(int l) const
void setModelEvalDescription(const std::string &modelEvalDescription)
RCP< LinearOpBase< Scalar > > get_W_op() const
RCP< const VectorBase< Scalar > > get_x() const
Thyra::ModelEvaluatorBase::OutArgs< Scalar > createOutArgsImpl() const
RCP< const VectorBase< Scalar > > get_p(int l) const
Teuchos::RCP< const Thyra::VectorSpaceBase< Scalar > > get_f_space() const
RCP< MultiVectorBase< Scalar > > getMultiVector() const