Tempus  Version of the Day
Time Integration
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Macros Pages
VanDerPolModel_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_VANDERPOL_MODEL_IMPL_HPP
10 #define TEMPUS_TEST_VANDERPOL_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 
22 #include <iostream>
23 
24 
25 namespace Tempus_Test {
26 
27 template<class Scalar>
29 VanDerPolModel(Teuchos::RCP<Teuchos::ParameterList> pList_)
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>
55 Thyra::ModelEvaluatorBase::InArgs<Scalar>
57 getExactSolution(double t) const
58 {
59  TEUCHOS_TEST_FOR_EXCEPTION( true, std::logic_error,
60  "Error - No exact solution for van der Pol problem!\n");
61  return(inArgs_);
62 }
63 
64 template<class Scalar>
65 Thyra::ModelEvaluatorBase::InArgs<Scalar>
67 getExactSensSolution(int j, double t) const
68 {
69  TEUCHOS_TEST_FOR_EXCEPTION( !isInitialized_, std::logic_error,
70  "Error - No exact sensitivities for van der Pol problem!\n");
71  return(inArgs_);
72 }
73 
74 template<class Scalar>
75 Teuchos::RCP<const Thyra::VectorSpaceBase<Scalar> >
77 get_x_space() const
78 {
79  return x_space_;
80 }
81 
82 
83 template<class Scalar>
84 Teuchos::RCP<const Thyra::VectorSpaceBase<Scalar> >
86 get_f_space() const
87 {
88  return f_space_;
89 }
90 
91 
92 template<class Scalar>
93 Thyra::ModelEvaluatorBase::InArgs<Scalar>
96 {
97  TEUCHOS_TEST_FOR_EXCEPTION( !isInitialized_, std::logic_error,
98  "Error, setupInOutArgs_ must be called first!\n");
99  return nominalValues_;
100 }
101 
102 
103 template<class Scalar>
104 Teuchos::RCP<Thyra::LinearOpWithSolveBase<Scalar> >
106 create_W() const
107 {
108  using Teuchos::RCP;
109  RCP<const Thyra::LinearOpWithSolveFactoryBase<Scalar> > W_factory = this->get_W_factory();
110  RCP<Thyra::LinearOpBase<Scalar> > matrix = this->create_W_op();
111  {
112  // 01/20/09 tscoffe: This is a total hack to provide a full rank matrix to
113  // linearOpWithSolve because it ends up factoring the matrix during
114  // initialization, which it really shouldn't do, or I'm doing something
115  // wrong here. The net effect is that I get exceptions thrown in
116  // optimized mode due to the matrix being rank deficient unless I do this.
117  RCP<Thyra::MultiVectorBase<Scalar> > multivec = Teuchos::rcp_dynamic_cast<Thyra::MultiVectorBase<Scalar> >(matrix,true);
118  {
119  RCP<Thyra::VectorBase<Scalar> > vec = Thyra::createMember(x_space_);
120  {
121  Thyra::DetachedVectorView<Scalar> vec_view( *vec );
122  vec_view[0] = 0.0;
123  vec_view[1] = 1.0;
124  }
125  V_V(multivec->col(0).ptr(),*vec);
126  {
127  Thyra::DetachedVectorView<Scalar> vec_view( *vec );
128  vec_view[0] = 1.0;
129  vec_view[1] = 0.0;
130  }
131  V_V(multivec->col(1).ptr(),*vec);
132  }
133  }
134  RCP<Thyra::LinearOpWithSolveBase<Scalar> > W =
135  Thyra::linearOpWithSolve<Scalar>(*W_factory, matrix);
136 
137  return W;
138 }
139 
140 
141 template<class Scalar>
142 Teuchos::RCP<Thyra::LinearOpBase<Scalar> >
144 create_W_op() const
145 {
146  Teuchos::RCP<Thyra::MultiVectorBase<Scalar> > matrix = Thyra::createMembers(x_space_, dim_);
147  return(matrix);
148 }
149 
150 
151 template<class Scalar>
152 Teuchos::RCP<const Thyra::LinearOpWithSolveFactoryBase<Scalar> >
155 {
156  Teuchos::RCP<Thyra::LinearOpWithSolveFactoryBase<Scalar> > W_factory =
157  Thyra::defaultSerialDenseLinearOpWithSolveFactory<Scalar>();
158  return W_factory;
159 }
160 
161 
162 template<class Scalar>
163 Thyra::ModelEvaluatorBase::InArgs<Scalar>
166 {
167  setupInOutArgs_();
168  return inArgs_;
169 }
170 
171 
172 // Private functions overridden from ModelEvaluatorDefaultBase
173 
174 
175 template<class Scalar>
176 Thyra::ModelEvaluatorBase::OutArgs<Scalar>
179 {
180  setupInOutArgs_();
181  return outArgs_;
182 }
183 
184 
185 template<class Scalar>
186 void
189  const Thyra::ModelEvaluatorBase::InArgs<Scalar> &inArgs,
190  const Thyra::ModelEvaluatorBase::OutArgs<Scalar> &outArgs
191  ) const
192 {
193  using Teuchos::RCP;
194  TEUCHOS_TEST_FOR_EXCEPTION( !isInitialized_, std::logic_error,
195  "Error, setupInOutArgs_ must be called first!\n");
196 
197  const RCP<const Thyra::VectorBase<Scalar> > x_in =
198  inArgs.get_x().assert_not_null();
199  Thyra::ConstDetachedVectorView<Scalar> x_in_view( *x_in );
200 
201  //double t = inArgs.get_t();
202  Scalar epsilon = epsilon_;
203  if (acceptModelParams_) {
204  const RCP<const Thyra::VectorBase<Scalar> > p_in =
205  inArgs.get_p(0).assert_not_null();
206  Thyra::ConstDetachedVectorView<Scalar> p_in_view( *p_in );
207  epsilon = p_in_view[0];
208  }
209 
210  Scalar beta = inArgs.get_beta();
211 
212  const RCP<Thyra::VectorBase<Scalar> > f_out = outArgs.get_f();
213  const RCP<Thyra::LinearOpBase<Scalar> > W_out = outArgs.get_W_op();
214  RCP<Thyra::MultiVectorBase<Scalar> > DfDp_out;
215  if (acceptModelParams_) {
216  Thyra::ModelEvaluatorBase::Derivative<Scalar> DfDp = outArgs.get_DfDp(0);
217  DfDp_out = DfDp.getMultiVector();
218  }
219 
220  if (inArgs.get_x_dot().is_null()) {
221 
222  // Evaluate the Explicit ODE f(x,t) [= xdot]
223  if (!is_null(f_out)) {
224  Thyra::DetachedVectorView<Scalar> f_out_view( *f_out );
225  f_out_view[0] = x_in_view[1];
226  f_out_view[1] =
227  ((1.0-x_in_view[0]*x_in_view[0])*x_in_view[1]-x_in_view[0])/epsilon;
228  }
229  if (!is_null(W_out)) {
230  RCP<Thyra::MultiVectorBase<Scalar> > matrix =
231  Teuchos::rcp_dynamic_cast<Thyra::MultiVectorBase<Scalar> >(W_out,true);
232  Thyra::DetachedMultiVectorView<Scalar> matrix_view( *matrix );
233  matrix_view(0,0) = 0.0; // d(f0)/d(x0_n)
234  matrix_view(0,1) = +beta; // d(f0)/d(x1_n)
235  matrix_view(1,0) =
236  -beta*(2.0*x_in_view[0]*x_in_view[1]+1.0)/epsilon; // d(f1)/d(x0_n)
237  matrix_view(1,1) =
238  beta*(1.0 - x_in_view[0]*x_in_view[0])/epsilon; // d(f1)/d(x1_n)
239  // Note: alpha = d(xdot)/d(x_n) and beta = d(x)/d(x_n)
240  }
241  if (!is_null(DfDp_out)) {
242  Thyra::DetachedMultiVectorView<Scalar> DfDp_out_view( *DfDp_out );
243  DfDp_out_view(0,0) = 0.0;
244  DfDp_out_view(1,0) =
245  -((1.0-x_in_view[0]*x_in_view[0])*x_in_view[1]-x_in_view[0])
246  /(epsilon*epsilon);
247  }
248  } else {
249 
250  // Evaluate the implicit ODE f(xdot, x, t) [= 0]
251  RCP<const Thyra::VectorBase<Scalar> > x_dot_in;
252  x_dot_in = inArgs.get_x_dot().assert_not_null();
253  Scalar alpha = inArgs.get_alpha();
254  if (!is_null(f_out)) {
255  Thyra::DetachedVectorView<Scalar> f_out_view( *f_out );
256  Thyra::ConstDetachedVectorView<Scalar> x_dot_in_view( *x_dot_in );
257  f_out_view[0] = x_dot_in_view[0] - x_in_view[1];
258  f_out_view[1] = x_dot_in_view[1]
259  - ((1.0-x_in_view[0]*x_in_view[0])*x_in_view[1]-x_in_view[0])/epsilon;
260  }
261  if (!is_null(W_out)) {
262  RCP<Thyra::MultiVectorBase<Scalar> > matrix =
263  Teuchos::rcp_dynamic_cast<Thyra::MultiVectorBase<Scalar> >(W_out,true);
264  Thyra::DetachedMultiVectorView<Scalar> matrix_view( *matrix );
265  matrix_view(0,0) = alpha; // d(f0)/d(x0_n)
266  matrix_view(0,1) = -beta; // d(f0)/d(x1_n)
267  matrix_view(1,0) =
268  beta*(2.0*x_in_view[0]*x_in_view[1]+1.0)/epsilon; // d(f1)/d(x0_n)
269  matrix_view(1,1) = alpha
270  - beta*(1.0 - x_in_view[0]*x_in_view[0])/epsilon; // d(f1)/d(x1_n)
271  // Note: alpha = d(xdot)/d(x_n) and beta = d(x)/d(x_n)
272  }
273  if (!is_null(DfDp_out)) {
274  Thyra::DetachedMultiVectorView<Scalar> DfDp_out_view( *DfDp_out );
275  DfDp_out_view(0,0) = 0.0;
276  DfDp_out_view(1,0) =
277  ((1.0-x_in_view[0]*x_in_view[0])*x_in_view[1]-x_in_view[0])
278  /(epsilon*epsilon);
279  }
280  }
281 }
282 
283 template<class Scalar>
284 Teuchos::RCP<const Thyra::VectorSpaceBase<Scalar> >
286 get_p_space(int l) const
287 {
288  if (!acceptModelParams_) {
289  return Teuchos::null;
290  }
291  TEUCHOS_ASSERT_IN_RANGE_UPPER_EXCLUSIVE( l, 0, Np_ );
292  return p_space_;
293 }
294 
295 template<class Scalar>
296 Teuchos::RCP<const Teuchos::Array<std::string> >
298 get_p_names(int l) const
299 {
300  if (!acceptModelParams_) {
301  return Teuchos::null;
302  }
303  TEUCHOS_ASSERT_IN_RANGE_UPPER_EXCLUSIVE( l, 0, Np_ );
304  Teuchos::RCP<Teuchos::Array<std::string> > p_strings =
305  Teuchos::rcp(new Teuchos::Array<std::string>());
306  p_strings->push_back("Model Coefficient: epsilon");
307  return p_strings;
308 }
309 
310 template<class Scalar>
311 Teuchos::RCP<const Thyra::VectorSpaceBase<Scalar> >
313 get_g_space(int j) const
314 {
315  TEUCHOS_ASSERT_IN_RANGE_UPPER_EXCLUSIVE( j, 0, Ng_ );
316  return g_space_;
317 }
318 
319 // private
320 
321 template<class Scalar>
322 void
325 {
326  if (isInitialized_) {
327  return;
328  }
329 
330  {
331  // Set up prototypical InArgs
332  Thyra::ModelEvaluatorBase::InArgsSetup<Scalar> inArgs;
333  inArgs.setModelEvalDescription(this->description());
334  inArgs.setSupports( Thyra::ModelEvaluatorBase::IN_ARG_t );
335  inArgs.setSupports( Thyra::ModelEvaluatorBase::IN_ARG_x );
336  inArgs.setSupports( Thyra::ModelEvaluatorBase::IN_ARG_beta );
337  inArgs.setSupports( Thyra::ModelEvaluatorBase::IN_ARG_x_dot );
338  inArgs.setSupports( Thyra::ModelEvaluatorBase::IN_ARG_alpha );
339  if (acceptModelParams_) {
340  inArgs.set_Np(Np_);
341  }
342  inArgs_ = inArgs;
343  }
344 
345  {
346  // Set up prototypical OutArgs
347  Thyra::ModelEvaluatorBase::OutArgsSetup<Scalar> outArgs;
348  outArgs.setModelEvalDescription(this->description());
349  outArgs.setSupports( Thyra::ModelEvaluatorBase::OUT_ARG_f );
350  outArgs.setSupports( Thyra::ModelEvaluatorBase::OUT_ARG_W_op );
351  if (acceptModelParams_) {
352  outArgs.set_Np_Ng(Np_,Ng_);
353  outArgs.setSupports( Thyra::ModelEvaluatorBase::OUT_ARG_DfDp,0,
354  Thyra::ModelEvaluatorBase::DERIV_MV_BY_COL );
355  }
356  outArgs_ = outArgs;
357  }
358 
359  // Set up nominal values
360  nominalValues_ = inArgs_;
361  if (haveIC_)
362  {
363  using Teuchos::RCP;
364  nominalValues_.set_t(t0_ic_);
365  const RCP<Thyra::VectorBase<Scalar> > x_ic = createMember(x_space_);
366  { // scope to delete DetachedVectorView
367  Thyra::DetachedVectorView<Scalar> x_ic_view( *x_ic );
368  x_ic_view[0] = x0_ic_;
369  x_ic_view[1] = x1_ic_;
370  }
371  nominalValues_.set_x(x_ic);
372  if (acceptModelParams_) {
373  const RCP<Thyra::VectorBase<Scalar> > p_ic = createMember(p_space_);
374  {
375  Thyra::DetachedVectorView<Scalar> p_ic_view( *p_ic );
376  p_ic_view[0] = epsilon_;
377  }
378  nominalValues_.set_p(0,p_ic);
379  }
380  const RCP<Thyra::VectorBase<Scalar> > x_dot_ic = createMember(x_space_);
381  { // scope to delete DetachedVectorView
382  Thyra::DetachedVectorView<Scalar> x_dot_ic_view( *x_dot_ic );
383  x_dot_ic_view[0] = x1_ic_;
384  x_dot_ic_view[1] = ((1.0-x0_ic_*x0_ic_)*x1_ic_-x0_ic_)/epsilon_;
385  }
386  nominalValues_.set_x_dot(x_dot_ic);
387  }
388 
389  isInitialized_ = true;
390 
391 }
392 
393 template<class Scalar>
394 void
396 setParameterList(Teuchos::RCP<Teuchos::ParameterList> const& paramList)
397 {
398  using Teuchos::get;
399  using Teuchos::ParameterList;
400  Teuchos::RCP<ParameterList> tmpPL = Teuchos::rcp(new ParameterList("VanDerPolModel"));
401  if (paramList != Teuchos::null) tmpPL = paramList;
402  tmpPL->validateParametersAndSetDefaults(*this->getValidParameters());
403  this->setMyParamList(tmpPL);
404  Teuchos::RCP<ParameterList> pl = this->getMyNonconstParamList();
405  bool acceptModelParams = get<bool>(*pl,"Accept model parameters");
406  bool haveIC = get<bool>(*pl,"Provide nominal values");
407  if ( (acceptModelParams != acceptModelParams_) ||
408  (haveIC != haveIC_)
409  ) {
410  isInitialized_ = false;
411  }
412  acceptModelParams_ = acceptModelParams;
413  haveIC_ = haveIC;
414  epsilon_ = get<Scalar>(*pl,"Coeff epsilon");
415  x0_ic_ = get<Scalar>(*pl,"IC x0");
416  x1_ic_ = get<Scalar>(*pl,"IC x1");
417  t0_ic_ = get<Scalar>(*pl,"IC t0");
418  setupInOutArgs_();
419 }
420 
421 template<class Scalar>
422 Teuchos::RCP<const Teuchos::ParameterList>
425 {
426  static Teuchos::RCP<const Teuchos::ParameterList> validPL;
427  if (is_null(validPL)) {
428  Teuchos::RCP<Teuchos::ParameterList> pl = Teuchos::parameterList();
429  pl->set("Accept model parameters", false);
430  pl->set("Provide nominal values", true);
431  Teuchos::setDoubleParameter(
432  "Coeff epsilon", 1.0e-06, "Coefficient a in model", &*pl);
433  Teuchos::setDoubleParameter(
434  "IC x0", 2.0, "Initial Condition for x0", &*pl);
435  Teuchos::setDoubleParameter(
436  "IC x1", 0.0, "Initial Condition for x1", &*pl);
437  Teuchos::setDoubleParameter(
438  "IC t0", 0.0, "Initial time t0", &*pl);
439  Teuchos::setIntParameter(
440  "Number of Time Step Sizes", 1, "Number time step sizes for convergence study", &*pl);
441  validPL = pl;
442  }
443  return validPL;
444 }
445 
446 } // namespace Tempus_Test
447 #endif // TEMPUS_TEST_VANDERPOL_MODEL_IMPL_HPP
Teuchos::RCP< const Teuchos::ParameterList > getValidParameters() 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)
Teuchos::RCP< Thyra::LinearOpBase< Scalar > > create_W_op() const
Teuchos::RCP< Thyra::LinearOpWithSolveBase< Scalar > > create_W() const
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 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
Thyra::ModelEvaluatorBase::OutArgs< Scalar > createOutArgsImpl() const
Teuchos::RCP< const Thyra::VectorSpaceBase< Scalar > > get_f_space() const