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: 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_IMEXPart_IMPLICITMODEL_IMPL_HPP
11 #define TEMPUS_TEST_VANDERPOL_IMEXPart_IMPLICITMODEL_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 #include "Thyra_MultiVectorStdOps.hpp"
23 #include "Thyra_DefaultMultiVectorProductVector.hpp"
24 
25 #include <iostream>
26 
27 namespace Tempus_Test {
28 
29 template <class Scalar>
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>
65 {
66  return x_space_;
67 }
68 
69 template <class Scalar>
72 {
73  return f_space_;
74 }
75 
76 template <class 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>
88 {
89  using Teuchos::RCP;
91  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.
100  Teuchos::rcp_dynamic_cast<Thyra::MultiVectorBase<Scalar> >(matrix,
101  true);
102  {
103  RCP<Thyra::VectorBase<Scalar> > vec = Thyra::createMember(x_space_);
104  {
105  Thyra::DetachedVectorView<Scalar> vec_view(*vec);
106  vec_view[0] = 1.0;
107  }
108  V_V(multivec->col(0).ptr(), *vec);
109  }
110  }
112  Thyra::linearOpWithSolve<Scalar>(*W_factory, matrix);
113 
114  return W;
115 }
116 
117 template <class Scalar>
120 {
122  Thyra::createMembers(x_space_, dim_);
123  return (matrix);
124 }
125 
126 template <class Scalar>
129 {
131  Thyra::defaultSerialDenseLinearOpWithSolveFactory<Scalar>();
132  return W_factory;
133 }
134 
135 template <class Scalar>
138 {
139  setupInOutArgs_();
140  return inArgs_;
141 }
142 
143 // Private functions overridden from ModelEvaluatorDefaultBase
144 
145 template <class Scalar>
148 {
149  setupInOutArgs_();
150  return outArgs_;
151 }
152 
153 template <class Scalar>
156  const Thyra::ModelEvaluatorBase::OutArgs<Scalar> &outArgs) const
157 {
159  using Teuchos::RCP;
160  using Teuchos::rcp_dynamic_cast;
161  TEUCHOS_TEST_FOR_EXCEPTION(!isInitialized_, std::logic_error,
162  "Error, setupInOutArgs_ must be called first!\n");
163 
165  inArgs.get_x().assert_not_null();
167  Scalar x1 = x_in_view[0];
168 
169  // double t = inArgs.get_t();
170  Scalar beta = inArgs.get_beta();
171  Scalar eps = epsilon_;
172 
173  const RCP<const Thyra::VectorBase<Scalar> > p_in = inArgs.get_p(0);
174  if (p_in != Teuchos::null) {
176  eps = p_in_view[0];
177  }
178 
179  RCP<const Thyra::MultiVectorBase<Scalar> > dxdp_in, dxdotdp_in, dydp_in;
180  if (inArgs.get_p(1) != Teuchos::null) {
181  dxdp_in =
182  rcp_dynamic_cast<const DMVPV>(inArgs.get_p(1), true)->getMultiVector();
183  }
184  if (inArgs.get_p(2) != Teuchos::null) {
185  dxdotdp_in =
186  rcp_dynamic_cast<const DMVPV>(inArgs.get_p(2), true)->getMultiVector();
187  }
188  if (inArgs.get_p(3) != Teuchos::null) {
189  dydp_in =
190  rcp_dynamic_cast<const DMVPV>(inArgs.get_p(3), true)->getMultiVector();
191  }
192 
194  inArgs.get_p(4).assert_not_null();
196  Scalar x0 = y_in_view[0];
197 
198  const RCP<Thyra::VectorBase<Scalar> > f_out = outArgs.get_f();
199  const RCP<Thyra::LinearOpBase<Scalar> > W_out = outArgs.get_W_op();
200  const RCP<Thyra::MultiVectorBase<Scalar> > DfDp_out =
201  outArgs.get_DfDp(0).getMultiVector();
202  const RCP<Thyra::MultiVectorBase<Scalar> > DfDx0_out =
203  outArgs.get_DfDp(4).getMultiVector();
204 
205  if (inArgs.get_x_dot().is_null()) {
206  // Evaluate the Explicit ODE f(x,t) [= xdot]
207  if (!is_null(f_out)) {
208  Thyra::DetachedVectorView<Scalar> f_out_view(*f_out);
209  f_out_view[0] = (1.0 - x0 * x0) * x1 / eps;
210  }
211  if (!is_null(DfDp_out)) {
212  Thyra::DetachedMultiVectorView<Scalar> DfDp_out_view(*DfDp_out);
213  DfDp_out_view(0, 0) = -(1.0 - x0 * x0) * x1 / (eps * eps);
214 
215  // Compute df/dp + (df/dx) * (dx/dp) + (df/dy)*dy/dp
216  if (useDfDpAsTangent_ && !is_null(dxdp_in)) {
218  DfDp_out_view(0, 0) += (1.0 - x0 * x0) / eps * dxdp(0, 0);
219  }
220  if (useDfDpAsTangent_ && !is_null(dydp_in)) {
222  DfDp_out_view(0, 0) += -2.0 * x0 * x1 / eps * dydp(0, 0);
223  }
224  }
225  if (!is_null(DfDx0_out)) {
226  Thyra::DetachedMultiVectorView<Scalar> DfDx0_out_view(*DfDx0_out);
227  DfDx0_out_view(0, 0) = -2.0 * x0 * x1 / eps; // d(f0)/d(x0_n);
228  }
229  if (!is_null(W_out)) {
231  Teuchos::rcp_dynamic_cast<Thyra::MultiVectorBase<Scalar> >(W_out,
232  true);
234  W_view(0, 0) = beta * (1.0 - x0 * x0) / eps; // d(f0)/d(x1_n)
235  // Note: alpha = d(xdot)/d(x_n) and beta = d(x)/d(x_n)
236  }
237  }
238  else {
239  // Evaluate the implicit ODE f(xdot, x, t) [= 0]
241  x_dot_in = inArgs.get_x_dot().assert_not_null();
242  if (!is_null(f_out)) {
243  Thyra::DetachedVectorView<Scalar> f_out_view(*f_out);
244  Thyra::ConstDetachedVectorView<Scalar> x_dot_in_view(*x_dot_in);
245  f_out_view[0] = x_dot_in_view[0] - (1.0 - x0 * x0) * x1 / eps;
246  }
247  if (!is_null(DfDp_out)) {
248  Thyra::DetachedMultiVectorView<Scalar> DfDp_out_view(*DfDp_out);
249  DfDp_out_view(0, 0) = (1.0 - x0 * x0) * x1 / (eps * eps);
250 
251  // Compute df/dp + (df/dx_dot)*(dx_dot/dp) + (df/dx)*(dx/dp) +
252  // (df/dy)*dy/dp
253  if (useDfDpAsTangent_ && !is_null(dxdotdp_in)) {
254  Thyra::ConstDetachedMultiVectorView<Scalar> dxdotdp(*dxdotdp_in);
255  DfDp_out_view(0, 0) += dxdotdp(0, 0);
256  }
257  if (useDfDpAsTangent_ && !is_null(dxdp_in)) {
259  DfDp_out_view(0, 0) += -(1.0 - x0 * x0) / eps * dxdp(0, 0);
260  }
261  if (useDfDpAsTangent_ && !is_null(dydp_in)) {
263  DfDp_out_view(0, 0) += 2.0 * x0 * x1 / eps * dydp(0, 0);
264  }
265  }
266  if (!is_null(DfDx0_out)) {
267  Thyra::DetachedMultiVectorView<Scalar> DfDx0_out_view(*DfDx0_out);
268  DfDx0_out_view(0, 0) = 2.0 * x0 * x1 / eps; // d(f0)/d(x0_n);
269  }
270  if (!is_null(W_out)) {
271  Scalar alpha = inArgs.get_alpha();
273  Teuchos::rcp_dynamic_cast<Thyra::MultiVectorBase<Scalar> >(W_out,
274  true);
276  W_view(0, 0) = alpha - beta * (1.0 - x0 * x0) / eps; // d(f0)/d(x1_n)
277  // Note: alpha = d(xdot)/d(x_n) and beta = d(x)/d(x_n)
278  }
279  }
280 }
281 
282 template <class Scalar>
285 {
287  if (l == 0)
288  return p_space_;
289  else if (l == 1 || l == 2)
290  return dxdp_space_;
291  else if (l == 3)
292  return dydp_space_;
293  else if (l == 4)
294  return y_space_;
295  return Teuchos::null;
296 }
297 
298 template <class Scalar>
301 {
305  if (l == 0)
306  p_strings->push_back("Model Coefficient: epsilon");
307  else if (l == 1)
308  p_strings->push_back("DxDp");
309  else if (l == 2)
310  p_strings->push_back("Dx_dotDp");
311  else if (l == 3)
312  p_strings->push_back("DyDp");
313  else if (l == 4)
314  p_strings->push_back("EXPLICIT_ONLY_VECTOR");
315  return p_strings;
316 }
317 
318 template <class Scalar>
321 {
323  return g_space_;
324 }
325 
326 template <class Scalar>
328 {
329  if (isInitialized_) return;
330 
331  {
332  // Set up prototypical InArgs
334  inArgs.setModelEvalDescription(this->description());
335  inArgs.setSupports(Thyra::ModelEvaluatorBase::IN_ARG_x_dot);
336  inArgs.setSupports(Thyra::ModelEvaluatorBase::IN_ARG_x);
337  inArgs.setSupports(Thyra::ModelEvaluatorBase::IN_ARG_t);
338  inArgs.setSupports(Thyra::ModelEvaluatorBase::IN_ARG_alpha);
339  inArgs.setSupports(Thyra::ModelEvaluatorBase::IN_ARG_beta);
340  inArgs.set_Np(Np_);
341  inArgs_ = inArgs;
342  }
343 
344  {
345  // Set up prototypical OutArgs
347  outArgs.setModelEvalDescription(this->description());
348  outArgs.setSupports(Thyra::ModelEvaluatorBase::OUT_ARG_f);
349  outArgs.setSupports(Thyra::ModelEvaluatorBase::OUT_ARG_W_op);
350  outArgs.set_Np_Ng(Np_, Ng_);
351  outArgs.setSupports(Thyra::ModelEvaluatorBase::OUT_ARG_DfDp, 0,
352  Thyra::ModelEvaluatorBase::DERIV_MV_JACOBIAN_FORM);
353  outArgs.setSupports(Thyra::ModelEvaluatorBase::OUT_ARG_DfDp, 4,
354  Thyra::ModelEvaluatorBase::DERIV_MV_JACOBIAN_FORM);
355  outArgs_ = outArgs;
356  }
357 
358  // Set up nominal values
359  nominalValues_ = inArgs_;
360  if (haveIC_) {
361  using Teuchos::RCP;
362  nominalValues_.set_t(t0_ic_);
363  const RCP<Thyra::VectorBase<Scalar> > x_ic = createMember(x_space_);
364  { // scope to delete DetachedVectorView
365  Thyra::DetachedVectorView<Scalar> x_ic_view(*x_ic);
366  x_ic_view[0] = x1_ic_;
367  }
368  nominalValues_.set_x(x_ic);
369 
370  const RCP<Thyra::VectorBase<Scalar> > p_ic = createMember(p_space_);
371  const RCP<Thyra::VectorBase<Scalar> > y_ic = createMember(y_space_);
372  {
373  Thyra::DetachedVectorView<Scalar> p_ic_view(*p_ic);
374  Thyra::DetachedVectorView<Scalar> y_ic_view(*y_ic);
375  p_ic_view[0] = epsilon_;
376  y_ic_view[0] = x0_ic_;
377  }
378  nominalValues_.set_p(0, p_ic);
379  nominalValues_.set_p(4, y_ic);
380 
381  const RCP<Thyra::VectorBase<Scalar> > x_dot_ic = createMember(x_space_);
382  { // scope to delete DetachedVectorView
383  Thyra::DetachedVectorView<Scalar> x_dot_ic_view(*x_dot_ic);
384  x_dot_ic_view[0] = ((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 template <class Scalar>
394  Teuchos::RCP<Teuchos::ParameterList> const &paramList)
395 {
396  using Teuchos::get;
399  Teuchos::rcp(new ParameterList("VanDerPol_IMEXPart_ImplicitModel"));
400  if (paramList != Teuchos::null) tmpPL = paramList;
401  tmpPL->validateParametersAndSetDefaults(*this->getValidParameters());
402  this->setMyParamList(tmpPL);
403  Teuchos::RCP<ParameterList> pl = this->getMyNonconstParamList();
404  bool haveIC = get<bool>(*pl, "Provide nominal values");
405  bool useDfDpAsTangent = get<bool>(*pl, "Use DfDp as Tangent");
406  if (haveIC != haveIC_) isInitialized_ = false;
407  haveIC_ = haveIC;
408  useDfDpAsTangent_ = useDfDpAsTangent;
409  epsilon_ = get<Scalar>(*pl, "Coeff epsilon");
410  x0_ic_ = get<Scalar>(*pl, "IC x0");
411  x1_ic_ = get<Scalar>(*pl, "IC x1");
412  t0_ic_ = get<Scalar>(*pl, "IC t0");
413  setupInOutArgs_();
414 }
415 
416 template <class Scalar>
419 {
421  if (is_null(validPL)) {
422  Teuchos::RCP<Teuchos::ParameterList> pl = Teuchos::parameterList();
423  pl->set("Accept model parameters", false);
424  pl->set("Provide nominal values", true);
425  pl->set("Use DfDp as Tangent", false);
426  Teuchos::setDoubleParameter("Coeff epsilon", 1.0e-06,
427  "Coefficient a in model", &*pl);
428  Teuchos::setDoubleParameter("IC x0", 2.0, "Initial Condition for x0", &*pl);
429  Teuchos::setDoubleParameter("IC x1", 0.0, "Initial Condition for x1", &*pl);
430  Teuchos::setDoubleParameter("IC t0", 0.0, "Initial time t0", &*pl);
431  validPL = pl;
432  }
433  return validPL;
434 }
435 
436 } // namespace Tempus_Test
437 #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