Tempus  Version of the Day
Time Integration
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups 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 namespace Tempus_Test {
27 
28 template <class Scalar>
31 {
32  isInitialized_ = false;
33  dim_ = 1;
34  Np_ = 5; // Number of parameter vectors (p, dx/dp, dx_dot/dp, dydp, y)
35  // Due to deficiencies in WrapperModelEvaluatorPairPartIMEX, y must
36  // be last
37  np_ = 1; // Number of parameters in this vector (for p and y) (1)
38  Ng_ = 0; // Number of observation functions (0)
39  ng_ = 0; // Number of elements in this observation function (1)
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 parameter spaces
51  y_space_ = Thyra::defaultSpmdVectorSpace<Scalar>(dim_);
52  p_space_ = Thyra::defaultSpmdVectorSpace<Scalar>(np_);
53  dxdp_space_ = Thyra::multiVectorProductVectorSpace(x_space_, np_);
54  dydp_space_ = Thyra::multiVectorProductVectorSpace(y_space_, np_);
55 
56  // g_space_ = Thyra::defaultSpmdVectorSpace<Scalar>(ng_);
57 
58  setParameterList(pList_);
59 }
60 
61 template <class Scalar>
64 {
65  return x_space_;
66 }
67 
68 template <class Scalar>
71 {
72  return f_space_;
73 }
74 
75 template <class Scalar>
78 {
79  TEUCHOS_TEST_FOR_EXCEPTION(!isInitialized_, std::logic_error,
80  "Error, setupInOutArgs_ must be called first!\n");
81  return nominalValues_;
82 }
83 
84 template <class Scalar>
87 {
88  using Teuchos::RCP;
90  this->get_W_factory();
91  RCP<Thyra::LinearOpBase<Scalar> > matrix = this->create_W_op();
92  {
93  // 01/20/09 tscoffe: This is a total hack to provide a full rank matrix to
94  // linearOpWithSolve because it ends up factoring the matrix during
95  // initialization, which it really shouldn't do, or I'm doing something
96  // wrong here. The net effect is that I get exceptions thrown in
97  // optimized mode due to the matrix being rank deficient unless I do this.
99  Teuchos::rcp_dynamic_cast<Thyra::MultiVectorBase<Scalar> >(matrix,
100  true);
101  {
102  RCP<Thyra::VectorBase<Scalar> > vec = Thyra::createMember(x_space_);
103  {
104  Thyra::DetachedVectorView<Scalar> vec_view(*vec);
105  vec_view[0] = 1.0;
106  }
107  V_V(multivec->col(0).ptr(), *vec);
108  }
109  }
111  Thyra::linearOpWithSolve<Scalar>(*W_factory, matrix);
112 
113  return W;
114 }
115 
116 template <class Scalar>
119 {
121  Thyra::createMembers(x_space_, dim_);
122  return (matrix);
123 }
124 
125 template <class Scalar>
128 {
130  Thyra::defaultSerialDenseLinearOpWithSolveFactory<Scalar>();
131  return W_factory;
132 }
133 
134 template <class Scalar>
137 {
138  setupInOutArgs_();
139  return inArgs_;
140 }
141 
142 // Private functions overridden from ModelEvaluatorDefaultBase
143 
144 template <class Scalar>
147 {
148  setupInOutArgs_();
149  return outArgs_;
150 }
151 
152 template <class Scalar>
155  const Thyra::ModelEvaluatorBase::OutArgs<Scalar> &outArgs) const
156 {
158  using Teuchos::RCP;
159  using Teuchos::rcp_dynamic_cast;
160  TEUCHOS_TEST_FOR_EXCEPTION(!isInitialized_, std::logic_error,
161  "Error, setupInOutArgs_ must be called first!\n");
162 
164  inArgs.get_x().assert_not_null();
166  Scalar x1 = x_in_view[0];
167 
168  // double t = inArgs.get_t();
169  Scalar beta = inArgs.get_beta();
170  Scalar eps = epsilon_;
171 
172  const RCP<const Thyra::VectorBase<Scalar> > p_in = inArgs.get_p(0);
173  if (p_in != Teuchos::null) {
175  eps = p_in_view[0];
176  }
177 
178  RCP<const Thyra::MultiVectorBase<Scalar> > dxdp_in, dxdotdp_in, dydp_in;
179  if (inArgs.get_p(1) != Teuchos::null) {
180  dxdp_in =
181  rcp_dynamic_cast<const DMVPV>(inArgs.get_p(1), true)->getMultiVector();
182  }
183  if (inArgs.get_p(2) != Teuchos::null) {
184  dxdotdp_in =
185  rcp_dynamic_cast<const DMVPV>(inArgs.get_p(2), true)->getMultiVector();
186  }
187  if (inArgs.get_p(3) != Teuchos::null) {
188  dydp_in =
189  rcp_dynamic_cast<const DMVPV>(inArgs.get_p(3), true)->getMultiVector();
190  }
191 
193  inArgs.get_p(4).assert_not_null();
195  Scalar x0 = y_in_view[0];
196 
197  const RCP<Thyra::VectorBase<Scalar> > f_out = outArgs.get_f();
198  const RCP<Thyra::LinearOpBase<Scalar> > W_out = outArgs.get_W_op();
199  const RCP<Thyra::MultiVectorBase<Scalar> > DfDp_out =
200  outArgs.get_DfDp(0).getMultiVector();
201  const RCP<Thyra::MultiVectorBase<Scalar> > DfDx0_out =
202  outArgs.get_DfDp(4).getMultiVector();
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] = (1.0 - x0 * x0) * x1 / eps;
209  }
210  if (!is_null(DfDp_out)) {
211  Thyra::DetachedMultiVectorView<Scalar> DfDp_out_view(*DfDp_out);
212  DfDp_out_view(0, 0) = -(1.0 - x0 * x0) * x1 / (eps * eps);
213 
214  // Compute df/dp + (df/dx) * (dx/dp) + (df/dy)*dy/dp
215  if (useDfDpAsTangent_ && !is_null(dxdp_in)) {
217  DfDp_out_view(0, 0) += (1.0 - x0 * x0) / eps * dxdp(0, 0);
218  }
219  if (useDfDpAsTangent_ && !is_null(dydp_in)) {
221  DfDp_out_view(0, 0) += -2.0 * x0 * x1 / eps * dydp(0, 0);
222  }
223  }
224  if (!is_null(DfDx0_out)) {
225  Thyra::DetachedMultiVectorView<Scalar> DfDx0_out_view(*DfDx0_out);
226  DfDx0_out_view(0, 0) = -2.0 * x0 * x1 / eps; // d(f0)/d(x0_n);
227  }
228  if (!is_null(W_out)) {
230  Teuchos::rcp_dynamic_cast<Thyra::MultiVectorBase<Scalar> >(W_out,
231  true);
233  W_view(0, 0) = beta * (1.0 - x0 * x0) / eps; // d(f0)/d(x1_n)
234  // Note: alpha = d(xdot)/d(x_n) and beta = d(x)/d(x_n)
235  }
236  }
237  else {
238  // Evaluate the implicit ODE f(xdot, x, t) [= 0]
240  x_dot_in = inArgs.get_x_dot().assert_not_null();
241  if (!is_null(f_out)) {
242  Thyra::DetachedVectorView<Scalar> f_out_view(*f_out);
243  Thyra::ConstDetachedVectorView<Scalar> x_dot_in_view(*x_dot_in);
244  f_out_view[0] = x_dot_in_view[0] - (1.0 - x0 * x0) * x1 / eps;
245  }
246  if (!is_null(DfDp_out)) {
247  Thyra::DetachedMultiVectorView<Scalar> DfDp_out_view(*DfDp_out);
248  DfDp_out_view(0, 0) = (1.0 - x0 * x0) * x1 / (eps * eps);
249 
250  // Compute df/dp + (df/dx_dot)*(dx_dot/dp) + (df/dx)*(dx/dp) +
251  // (df/dy)*dy/dp
252  if (useDfDpAsTangent_ && !is_null(dxdotdp_in)) {
253  Thyra::ConstDetachedMultiVectorView<Scalar> dxdotdp(*dxdotdp_in);
254  DfDp_out_view(0, 0) += dxdotdp(0, 0);
255  }
256  if (useDfDpAsTangent_ && !is_null(dxdp_in)) {
258  DfDp_out_view(0, 0) += -(1.0 - x0 * x0) / eps * dxdp(0, 0);
259  }
260  if (useDfDpAsTangent_ && !is_null(dydp_in)) {
262  DfDp_out_view(0, 0) += 2.0 * x0 * x1 / eps * dydp(0, 0);
263  }
264  }
265  if (!is_null(DfDx0_out)) {
266  Thyra::DetachedMultiVectorView<Scalar> DfDx0_out_view(*DfDx0_out);
267  DfDx0_out_view(0, 0) = 2.0 * x0 * x1 / eps; // d(f0)/d(x0_n);
268  }
269  if (!is_null(W_out)) {
270  Scalar alpha = inArgs.get_alpha();
272  Teuchos::rcp_dynamic_cast<Thyra::MultiVectorBase<Scalar> >(W_out,
273  true);
275  W_view(0, 0) = alpha - beta * (1.0 - x0 * x0) / eps; // d(f0)/d(x1_n)
276  // Note: alpha = d(xdot)/d(x_n) and beta = d(x)/d(x_n)
277  }
278  }
279 }
280 
281 template <class Scalar>
284 {
286  if (l == 0)
287  return p_space_;
288  else if (l == 1 || l == 2)
289  return dxdp_space_;
290  else if (l == 3)
291  return dydp_space_;
292  else if (l == 4)
293  return y_space_;
294  return Teuchos::null;
295 }
296 
297 template <class Scalar>
300 {
304  if (l == 0)
305  p_strings->push_back("Model Coefficient: epsilon");
306  else if (l == 1)
307  p_strings->push_back("DxDp");
308  else if (l == 2)
309  p_strings->push_back("Dx_dotDp");
310  else if (l == 3)
311  p_strings->push_back("DyDp");
312  else if (l == 4)
313  p_strings->push_back("EXPLICIT_ONLY_VECTOR");
314  return p_strings;
315 }
316 
317 template <class Scalar>
320 {
322  return g_space_;
323 }
324 
325 template <class Scalar>
327 {
328  if (isInitialized_) return;
329 
330  {
331  // Set up prototypical InArgs
333  inArgs.setModelEvalDescription(this->description());
334  inArgs.setSupports(Thyra::ModelEvaluatorBase::IN_ARG_x_dot);
335  inArgs.setSupports(Thyra::ModelEvaluatorBase::IN_ARG_x);
336  inArgs.setSupports(Thyra::ModelEvaluatorBase::IN_ARG_t);
337  inArgs.setSupports(Thyra::ModelEvaluatorBase::IN_ARG_alpha);
338  inArgs.setSupports(Thyra::ModelEvaluatorBase::IN_ARG_beta);
339  inArgs.set_Np(Np_);
340  inArgs_ = inArgs;
341  }
342 
343  {
344  // Set up prototypical OutArgs
346  outArgs.setModelEvalDescription(this->description());
347  outArgs.setSupports(Thyra::ModelEvaluatorBase::OUT_ARG_f);
348  outArgs.setSupports(Thyra::ModelEvaluatorBase::OUT_ARG_W_op);
349  outArgs.set_Np_Ng(Np_, Ng_);
350  outArgs.setSupports(Thyra::ModelEvaluatorBase::OUT_ARG_DfDp, 0,
351  Thyra::ModelEvaluatorBase::DERIV_MV_JACOBIAN_FORM);
352  outArgs.setSupports(Thyra::ModelEvaluatorBase::OUT_ARG_DfDp, 4,
353  Thyra::ModelEvaluatorBase::DERIV_MV_JACOBIAN_FORM);
354  outArgs_ = outArgs;
355  }
356 
357  // Set up nominal values
358  nominalValues_ = inArgs_;
359  if (haveIC_) {
360  using Teuchos::RCP;
361  nominalValues_.set_t(t0_ic_);
362  const RCP<Thyra::VectorBase<Scalar> > x_ic = createMember(x_space_);
363  { // scope to delete DetachedVectorView
364  Thyra::DetachedVectorView<Scalar> x_ic_view(*x_ic);
365  x_ic_view[0] = x1_ic_;
366  }
367  nominalValues_.set_x(x_ic);
368 
369  const RCP<Thyra::VectorBase<Scalar> > p_ic = createMember(p_space_);
370  const RCP<Thyra::VectorBase<Scalar> > y_ic = createMember(y_space_);
371  {
372  Thyra::DetachedVectorView<Scalar> p_ic_view(*p_ic);
373  Thyra::DetachedVectorView<Scalar> y_ic_view(*y_ic);
374  p_ic_view[0] = epsilon_;
375  y_ic_view[0] = x0_ic_;
376  }
377  nominalValues_.set_p(0, p_ic);
378  nominalValues_.set_p(4, y_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] = ((1.0 - x0_ic_ * x0_ic_) * x1_ic_ - x0_ic_) / epsilon_;
384  }
385  nominalValues_.set_x_dot(x_dot_ic);
386  }
387 
388  isInitialized_ = true;
389 }
390 
391 template <class Scalar>
393  Teuchos::RCP<Teuchos::ParameterList> const &paramList)
394 {
395  using Teuchos::get;
398  Teuchos::rcp(new ParameterList("VanDerPol_IMEXPart_ImplicitModel"));
399  if (paramList != Teuchos::null) tmpPL = paramList;
400  tmpPL->validateParametersAndSetDefaults(*this->getValidParameters());
401  this->setMyParamList(tmpPL);
402  Teuchos::RCP<ParameterList> pl = this->getMyNonconstParamList();
403  bool haveIC = get<bool>(*pl, "Provide nominal values");
404  bool useDfDpAsTangent = get<bool>(*pl, "Use DfDp as Tangent");
405  if (haveIC != haveIC_) isInitialized_ = false;
406  haveIC_ = haveIC;
407  useDfDpAsTangent_ = useDfDpAsTangent;
408  epsilon_ = get<Scalar>(*pl, "Coeff epsilon");
409  x0_ic_ = get<Scalar>(*pl, "IC x0");
410  x1_ic_ = get<Scalar>(*pl, "IC x1");
411  t0_ic_ = get<Scalar>(*pl, "IC t0");
412  setupInOutArgs_();
413 }
414 
415 template <class Scalar>
418 {
420  if (is_null(validPL)) {
421  Teuchos::RCP<Teuchos::ParameterList> pl = Teuchos::parameterList();
422  pl->set("Accept model parameters", false);
423  pl->set("Provide nominal values", true);
424  pl->set("Use DfDp as Tangent", false);
425  Teuchos::setDoubleParameter("Coeff epsilon", 1.0e-06,
426  "Coefficient a in model", &*pl);
427  Teuchos::setDoubleParameter("IC x0", 2.0, "Initial Condition for x0", &*pl);
428  Teuchos::setDoubleParameter("IC x1", 0.0, "Initial Condition for x1", &*pl);
429  Teuchos::setDoubleParameter("IC t0", 0.0, "Initial time t0", &*pl);
430  validPL = pl;
431  }
432  return validPL;
433 }
434 
435 } // namespace Tempus_Test
436 #endif // TEMPUS_TEST_VANDERPOL_IMEXPart_IMPLICITMODEL_IMPL_HPP
Thyra::ModelEvaluatorBase::InArgs< Scalar > getNominalValues() const
bool is_null(const boost::shared_ptr< T > &p)
Teuchos::RCP< const Teuchos::ParameterList > getValidParameters() const
Derivative< Scalar > get_DfDp(int l) const
RCP< const VectorBase< Scalar > > get_x_dot() const
void setModelEvalDescription(const std::string &modelEvalDescription)
Teuchos::RCP< Thyra::LinearOpWithSolveBase< Scalar > > create_W() const
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)
VanDerPol_IMEXPart_ImplicitModel(Teuchos::RCP< Teuchos::ParameterList > pList=Teuchos::null)
Evaluation< VectorBase< Scalar > > get_f() const
Teuchos::RCP< const Thyra::VectorSpaceBase< Scalar > > get_g_space(int j) const
#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)
Ptr< T > ptr() const
Thyra::ModelEvaluatorBase::OutArgs< Scalar > createOutArgsImpl() const
Teuchos::RCP< Thyra::LinearOpBase< Scalar > > create_W_op() const
void validateParametersAndSetDefaults(ParameterList const &validParamList, int const depth=1000)
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
void setSupports(EOutArgsMembers arg, bool supports=true)
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
void setModelEvalDescription(const std::string &modelEvalDescription)
Teuchos::RCP< const Thyra::LinearOpWithSolveFactoryBase< Scalar > > get_W_factory() const
RCP< LinearOpBase< Scalar > > get_W_op() const
RCP< const VectorBase< Scalar > > get_x() const
void evalModelImpl(const Thyra::ModelEvaluatorBase::InArgs< Scalar > &inArgs_bar, const Thyra::ModelEvaluatorBase::OutArgs< Scalar > &outArgs_bar) const
Thyra::ModelEvaluatorBase::InArgs< Scalar > createInArgs() const
RCP< const VectorBase< Scalar > > get_p(int l) const