Tempus  Version of the Day
Time Integration
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
VanDerPol_IMEX_ExplicitModel_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_IMEX_EXPLICITMODEL_IMPL_HPP
11 #define TEMPUS_TEST_VANDERPOL_IMEX_EXPLICITMODEL_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 #include "Thyra_DefaultBlockedLinearOp.hpp"
25 
26 #include <iostream>
27 
28 namespace Tempus_Test {
29 
30 template <class Scalar>
32  Teuchos::RCP<Teuchos::ParameterList> pList_, bool useProductVector)
33 {
34  useProductVector_ = useProductVector;
35  isInitialized_ = false;
36  dim_ = 2;
37  Np_ = 3; // Number of parameter vectors (p, dx/dp, dx_dot/dp)
38  np_ = 1; // Number of parameters in this vector (1)
39  Ng_ = 0; // Number of observation functions (0)
40  ng_ = 0; // Number of elements in this observation function (0)
41  acceptModelParams_ = 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  if (useProductVector_ == false) {
49  x_space_ = Thyra::defaultSpmdVectorSpace<Scalar>(dim_);
50  f_space_ = Thyra::defaultSpmdVectorSpace<Scalar>(dim_);
51  }
52  else {
53  // Create using ProductVector so we can use in partitioned IMEX-RK Stepper.
54  using Teuchos::RCP;
56  xSpace_ = Thyra::defaultSpmdVectorSpace<Scalar>(1);
57  for (int i = 0; i < dim_; ++i) yxSpaceArray.push_back(xSpace_);
58  x_space_ = Thyra::productVectorSpace<Scalar>(yxSpaceArray);
59  f_space_ = Thyra::productVectorSpace<Scalar>(yxSpaceArray);
60  }
61 
62  // Create p_space and g_space
63  p_space_ = Thyra::defaultSpmdVectorSpace<Scalar>(np_);
64  g_space_ = Thyra::defaultSpmdVectorSpace<Scalar>(ng_);
65  dxdp_space_ = Thyra::multiVectorProductVectorSpace(x_space_, np_);
66 
67  setParameterList(pList_);
68 }
69 
70 template <class Scalar>
73 {
74  return x_space_;
75 }
76 
77 template <class Scalar>
80 {
81  return f_space_;
82 }
83 
84 template <class Scalar>
87 {
88  TEUCHOS_TEST_FOR_EXCEPTION(!isInitialized_, std::logic_error,
89  "Error, setupInOutArgs_ must be called first!\n");
90  return nominalValues_;
91 }
92 
93 template <class Scalar>
96 {
98  if (useProductVector_ == false)
99  W_op = Thyra::createMembers(x_space_, dim_);
100  else {
101  // When using product vector formulation, we have to create a block
102  // operator so that its domain space is consistent with x_space_
104  Thyra::createMembers(xSpace_, 1);
106  Thyra::createMembers(xSpace_, 1);
108  Thyra::createMembers(xSpace_, 1);
110  Thyra::createMembers(xSpace_, 1);
111  W_op = Thyra::nonconstBlock2x2(A00, A01, A10, A11);
112  }
113  return W_op;
114 }
115 
116 template <class Scalar>
119 {
120  setupInOutArgs_();
121  return inArgs_;
122 }
123 
124 template <class Scalar>
127 {
128  setupInOutArgs_();
129  return outArgs_;
130 }
131 
132 template <class Scalar>
134 {
135  if (isInitialized_) {
136  return;
137  }
138 
139  {
140  // Set up prototypical InArgs
142  inArgs.setModelEvalDescription(this->description());
143  inArgs.setSupports(Thyra::ModelEvaluatorBase::IN_ARG_x_dot);
144  inArgs.setSupports(Thyra::ModelEvaluatorBase::IN_ARG_x);
145  inArgs.setSupports(Thyra::ModelEvaluatorBase::IN_ARG_t);
146  inArgs.setSupports(Thyra::ModelEvaluatorBase::IN_ARG_alpha);
147  inArgs.setSupports(Thyra::ModelEvaluatorBase::IN_ARG_beta);
148  if (acceptModelParams_) {
149  inArgs.set_Np(Np_);
150  }
151  inArgs_ = inArgs;
152  }
153 
154  {
155  // Set up prototypical OutArgs
157  outArgs.setModelEvalDescription(this->description());
158  outArgs.setSupports(Thyra::ModelEvaluatorBase::OUT_ARG_f);
159  outArgs.setSupports(Thyra::ModelEvaluatorBase::OUT_ARG_W_op);
160  if (acceptModelParams_) {
161  outArgs.set_Np_Ng(Np_, Ng_);
162  outArgs.setSupports(Thyra::ModelEvaluatorBase::OUT_ARG_DfDp, 0,
163  Thyra::ModelEvaluatorBase::DERIV_MV_JACOBIAN_FORM);
164  }
165  outArgs_ = outArgs;
166  }
167 
168  // Set up nominal values
169  nominalValues_ = inArgs_;
170  if (haveIC_) {
171  using Teuchos::RCP;
172 
173  nominalValues_.set_t(t0_ic_);
174  const RCP<Thyra::VectorBase<Scalar> > x_ic = createMember(x_space_);
175  { // scope to delete DetachedVectorView
176  Thyra::DetachedVectorView<Scalar> x_ic_view(*x_ic);
177  x_ic_view[0] = x0_ic_;
178  x_ic_view[1] = x1_ic_;
179  }
180  nominalValues_.set_x(x_ic);
181 
182  if (acceptModelParams_) {
183  const RCP<Thyra::VectorBase<Scalar> > p_ic = createMember(p_space_);
184  {
185  Thyra::DetachedVectorView<Scalar> p_ic_view(*p_ic);
186  p_ic_view[0] = epsilon_;
187  }
188  nominalValues_.set_p(0, p_ic);
189  }
190  const RCP<Thyra::VectorBase<Scalar> > x_dot_ic = createMember(x_space_);
191  { // scope to delete DetachedVectorView
192  Thyra::DetachedVectorView<Scalar> x_dot_ic_view(*x_dot_ic);
193  x_dot_ic_view[0] = x1_ic_;
194  x_dot_ic_view[1] = ((1.0 - x0_ic_ * x0_ic_) * x1_ic_ - x0_ic_) / epsilon_;
195  }
196  nominalValues_.set_x_dot(x_dot_ic);
197  }
198 
199  isInitialized_ = true;
200 }
201 
202 template <class Scalar>
204  Teuchos::RCP<Teuchos::ParameterList> const &paramList)
205 {
206  using Teuchos::get;
209  Teuchos::rcp(new ParameterList("VanDerPol_IMEX_ExplicitModel"));
210  if (paramList != Teuchos::null) tmpPL = paramList;
211  tmpPL->validateParametersAndSetDefaults(*this->getValidParameters());
212  this->setMyParamList(tmpPL);
213  Teuchos::RCP<ParameterList> pl = this->getMyNonconstParamList();
214  bool acceptModelParams = get<bool>(*pl, "Accept model parameters");
215  bool haveIC = get<bool>(*pl, "Provide nominal values");
216  bool useDfDpAsTangent = get<bool>(*pl, "Use DfDp as Tangent");
217  if ((acceptModelParams != acceptModelParams_) || (haveIC != haveIC_)) {
218  isInitialized_ = false;
219  }
220  acceptModelParams_ = acceptModelParams;
221  haveIC_ = haveIC;
222  useDfDpAsTangent_ = useDfDpAsTangent;
223  epsilon_ = get<Scalar>(*pl, "Coeff epsilon");
224  x0_ic_ = get<Scalar>(*pl, "IC x0");
225  x1_ic_ = get<Scalar>(*pl, "IC x1");
226  t0_ic_ = get<Scalar>(*pl, "IC t0");
227  setupInOutArgs_();
228 }
229 
230 template <class Scalar>
233 {
235  if (is_null(validPL)) {
236  Teuchos::RCP<Teuchos::ParameterList> pl = Teuchos::parameterList();
237  pl->set("Accept model parameters", false);
238  pl->set("Provide nominal values", true);
239  pl->set("Use DfDp as Tangent", false);
240  Teuchos::setDoubleParameter("Coeff epsilon", 1.0e-06,
241  "Coefficient a in model", &*pl);
242  Teuchos::setDoubleParameter("IC x0", 2.0, "Initial Condition for x0", &*pl);
243  Teuchos::setDoubleParameter("IC x1", 0.0, "Initial Condition for x1", &*pl);
244  Teuchos::setDoubleParameter("IC t0", 0.0, "Initial time t0", &*pl);
245  validPL = pl;
246  }
247  return validPL;
248 }
249 
250 template <class Scalar>
253 {
254  if (!acceptModelParams_) {
255  return Teuchos::null;
256  }
258  if (l == 0)
259  return p_space_;
260  else if (l == 1 || l == 2)
261  return dxdp_space_;
262  return Teuchos::null;
263 }
264 
265 template <class Scalar>
268 {
269  if (!acceptModelParams_) {
270  return Teuchos::null;
271  }
275  if (l == 0)
276  p_strings->push_back("Model Coefficient: epsilon");
277  else if (l == 1)
278  p_strings->push_back("DxDp");
279  else if (l == 2)
280  p_strings->push_back("Dx_dotDp");
281  return p_strings;
282 }
283 
284 template <class Scalar>
287 {
289  return g_space_;
290 }
291 
292 template <class Scalar>
295  const Thyra::ModelEvaluatorBase::OutArgs<Scalar> &outArgs) const
296 {
298  using Teuchos::RCP;
299  using Teuchos::rcp_dynamic_cast;
300  TEUCHOS_TEST_FOR_EXCEPTION(!isInitialized_, std::logic_error,
301  "Error, setupInOutArgs_ must be called first!\n");
302 
304  inArgs.get_x().assert_not_null();
306 
307  // double t = inArgs.get_t();
308  Scalar beta = inArgs.get_beta();
309  Scalar eps = epsilon_;
310  RCP<const Thyra::MultiVectorBase<Scalar> > dxdp_in, dxdotdp_in;
311  if (acceptModelParams_) {
312  const RCP<const Thyra::VectorBase<Scalar> > p_in = inArgs.get_p(0);
313  if (p_in != Teuchos::null) {
315  eps = p_in_view[0];
316  }
317  if (inArgs.get_p(1) != Teuchos::null) {
318  dxdp_in = rcp_dynamic_cast<const DMVPV>(inArgs.get_p(1), true)
319  ->getMultiVector();
320  }
321  if (inArgs.get_p(2) != Teuchos::null) {
322  dxdotdp_in = rcp_dynamic_cast<const DMVPV>(inArgs.get_p(2), true)
323  ->getMultiVector();
324  }
325  }
326 
327  const RCP<Thyra::VectorBase<Scalar> > f_out = outArgs.get_f();
328  const RCP<Thyra::LinearOpBase<Scalar> > W_out = outArgs.get_W_op();
330  if (acceptModelParams_) {
332  DfDp_out = DfDp.getMultiVector();
333  }
334 
335  if (inArgs.get_x_dot().is_null()) {
336  // Evaluate the Explicit ODE f(x,t) [= xdot]
337  if (!is_null(f_out)) {
338  Thyra::DetachedVectorView<Scalar> f_out_view(*f_out);
339  f_out_view[0] = x_in_view[1];
340  f_out_view[1] = -x_in_view[0] / eps;
341  }
342  if (!is_null(DfDp_out)) {
343  Thyra::DetachedMultiVectorView<Scalar> DfDp_out_view(*DfDp_out);
344  DfDp_out_view(0, 0) = 0.0;
345  DfDp_out_view(1, 0) = x_in_view[0] / (eps * eps);
346 
347  // Compute df/dp + (df/dx) * (dx/dp)
348  if (useDfDpAsTangent_ && !is_null(dxdp_in)) {
350  DfDp_out_view(0, 0) += dxdp(1, 0);
351  DfDp_out_view(1, 0) += -dxdp(0, 0) / eps;
352  }
353  }
354  if (!is_null(W_out)) {
355  if (useProductVector_ == false) {
357  Teuchos::rcp_dynamic_cast<Thyra::MultiVectorBase<Scalar> >(W_out,
358  true);
360  W_view(0, 0) = 0.0; // d(f0)/d(x0_n)
361  W_view(0, 1) = beta; // d(f0)/d(x1_n)
362  W_view(1, 0) = -beta / eps; // d(f1)/d(x0_n)
363  W_view(1, 1) = 0.0; // d(f1)/d(x1_n)
364  // Note: alpha = d(xdot)/d(x_n) and beta = d(x)/d(x_n)
365  }
366  else {
368  Teuchos::rcp_dynamic_cast<Thyra::DefaultBlockedLinearOp<Scalar> >(
369  W_out, true);
371  Teuchos::rcp_dynamic_cast<Thyra::MultiVectorBase<Scalar> >(
372  W_block->getNonconstBlock(0, 0), true);
374  Teuchos::rcp_dynamic_cast<Thyra::MultiVectorBase<Scalar> >(
375  W_block->getNonconstBlock(0, 1), true);
377  Teuchos::rcp_dynamic_cast<Thyra::MultiVectorBase<Scalar> >(
378  W_block->getNonconstBlock(1, 0), true);
380  Teuchos::rcp_dynamic_cast<Thyra::MultiVectorBase<Scalar> >(
381  W_block->getNonconstBlock(1, 1), true);
386  W00_view(0, 0) = 0.0; // d(f0)/d(x0_n)
387  W01_view(0, 0) = beta; // d(f0)/d(x1_n)
388  W10_view(0, 0) = -beta / eps; // d(f1)/d(x0_n)
389  W11_view(0, 0) = 0.0; // d(f1)/d(x1_n)
390  // Note: alpha = d(xdot)/d(x_n) and beta = d(x)/d(x_n)
391  }
392  }
393  }
394  else {
395  // Evaluate the implicit ODE f(xdot, x, t) [= 0]
397  x_dot_in = inArgs.get_x_dot().assert_not_null();
398  Scalar alpha = inArgs.get_alpha();
399 
400  if (!is_null(f_out)) {
401  Thyra::DetachedVectorView<Scalar> f_out_view(*f_out);
402  Thyra::ConstDetachedVectorView<Scalar> x_dot_in_view(*x_dot_in);
403  f_out_view[0] = x_dot_in_view[0] - x_in_view[1];
404  f_out_view[1] = x_dot_in_view[1] + x_in_view[0] / eps;
405  }
406  if (!is_null(DfDp_out)) {
407  Thyra::DetachedMultiVectorView<Scalar> DfDp_out_view(*DfDp_out);
408  DfDp_out_view(0, 0) = 0.0;
409  DfDp_out_view(1, 0) = -x_in_view[0] / (eps * eps);
410 
411  // Compute df/dp + (df/dx_dot)*(dx_dot/dp) + (df/dx) * (dx/dp)
412  if (useDfDpAsTangent_ && !is_null(dxdotdp_in)) {
413  Thyra::ConstDetachedMultiVectorView<Scalar> dxdotdp(*dxdotdp_in);
414  DfDp_out_view(0, 0) += dxdotdp(0, 0);
415  DfDp_out_view(1, 0) += dxdotdp(1, 0);
416  }
417  if (useDfDpAsTangent_ && !is_null(dxdp_in)) {
419  DfDp_out_view(0, 0) += -dxdp(1, 0);
420  DfDp_out_view(1, 0) += dxdp(0, 0) / eps;
421  }
422  }
423  if (!is_null(W_out)) {
424  if (useProductVector_ == false) {
426  Teuchos::rcp_dynamic_cast<Thyra::MultiVectorBase<Scalar> >(W_out,
427  true);
429  W_view(0, 0) = alpha; // d(f0)/d(x0_n)
430  W_view(0, 1) = -beta; // d(f0)/d(x1_n)
431  W_view(1, 0) = beta / eps; // d(f1)/d(x0_n)
432  W_view(1, 1) = alpha; // d(f1)/d(x1_n)
433  // Note: alpha = d(xdot)/d(x_n) and beta = d(x)/d(x_n)
434  }
435  else {
437  Teuchos::rcp_dynamic_cast<Thyra::DefaultBlockedLinearOp<Scalar> >(
438  W_out, true);
440  Teuchos::rcp_dynamic_cast<Thyra::MultiVectorBase<Scalar> >(
441  W_block->getNonconstBlock(0, 0), true);
443  Teuchos::rcp_dynamic_cast<Thyra::MultiVectorBase<Scalar> >(
444  W_block->getNonconstBlock(0, 1), true);
446  Teuchos::rcp_dynamic_cast<Thyra::MultiVectorBase<Scalar> >(
447  W_block->getNonconstBlock(1, 0), true);
449  Teuchos::rcp_dynamic_cast<Thyra::MultiVectorBase<Scalar> >(
450  W_block->getNonconstBlock(1, 1), true);
455  W00_view(0, 0) = alpha; // d(f0)/d(x0_n)
456  W01_view(0, 0) = -beta; // d(f0)/d(x1_n)
457  W10_view(0, 0) = beta / eps; // d(f1)/d(x0_n)
458  W11_view(0, 0) = alpha; // d(f1)/d(x1_n)
459  // Note: alpha = d(xdot)/d(x_n) and beta = d(x)/d(x_n)
460  }
461  }
462  }
463 }
464 
465 } // namespace Tempus_Test
466 #endif // TEMPUS_TEST_VANDERPOL_IMEX_EXPLICITMODEL_IMPL_HPP
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< const Thyra::VectorSpaceBase< Scalar > > get_f_space() const
ParameterList & set(std::string const &name, T const &value, std::string const &docString="", RCP< const ParameterEntryValidator > const &validator=null)
Teuchos::RCP< const Thyra::VectorSpaceBase< Scalar > > get_x_space() const
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
Teuchos::RCP< const Thyra::VectorSpaceBase< Scalar > > get_p_space(int l) const
Evaluation< VectorBase< Scalar > > get_f() const
void setParameterList(Teuchos::RCP< Teuchos::ParameterList > const &paramList)
#define TEUCHOS_ASSERT_IN_RANGE_UPPER_EXCLUSIVE(index, lower_inclusive, upper_exclusive)
void setSupports(EInArgsMembers arg, bool supports=true)
Teuchos::RCP< Thyra::LinearOpBase< Scalar > > create_W_op() const
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
Thyra::ModelEvaluatorBase::InArgs< Scalar > getNominalValues() const
Thyra::ModelEvaluatorBase::InArgs< Scalar > createInArgs() const
void validateParametersAndSetDefaults(ParameterList const &validParamList, int const depth=1000)
void push_back(const value_type &x)
void setSupports(EOutArgsMembers arg, bool supports=true)
Thyra::ModelEvaluatorBase::OutArgs< Scalar > createOutArgsImpl() const
VanDerPol_IMEX_ExplicitModel(Teuchos::RCP< Teuchos::ParameterList > pList=Teuchos::null, bool useProductVector=false)
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 setModelEvalDescription(const std::string &modelEvalDescription)
RCP< LinearOpBase< Scalar > > get_W_op() const
RCP< const VectorBase< Scalar > > get_x() const
RCP< const VectorBase< Scalar > > get_p(int l) const
RCP< MultiVectorBase< Scalar > > getMultiVector() const