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