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