Tempus  Version of the Day
Time Integration
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
SteadyQuadraticModel_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_STEADY_QUADRATIC_MODEL_IMPL_HPP
10 #define TEMPUS_TEST_STEADY_QUADRATIC_MODEL_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_MultiVectorStdOps.hpp"
21 #include "Thyra_VectorStdOps.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_ = 3; // Number of parameter vectors (p, dx/dp, dx_dot/dp)
35  np_ = 1; // Number of parameters in this vector (3)
36  Ng_ = 1; // Number of observation functions (1)
37  ng_ = 1; // Number of elements in this observation function (1)
38  useDfDpAsTangent_ = false;
39  b_ = 1.0;
40 
41  // Create x_space and f_space
42  x_space_ = Thyra::defaultSpmdVectorSpace<Scalar>(dim_);
43  f_space_ = Thyra::defaultSpmdVectorSpace<Scalar>(dim_);
44  // Create p_space and g_space
45  p_space_ = Thyra::defaultSpmdVectorSpace<Scalar>(np_);
46  g_space_ = Thyra::defaultSpmdVectorSpace<Scalar>(ng_);
47 
48  setParameterList(pList_);
49 
50  // Create DxDp product space
51  DxDp_space_ = Thyra::multiVectorProductVectorSpace(x_space_, np_);
52 }
53 
54 template <class Scalar>
56 {
57  if (b_ < 0.0) return b_;
58  return -b_;
59 }
60 
61 template <class Scalar>
63 {
64  if (b_ < 0.0) return 1.0;
65  return -1.0;
66 }
67 
68 template <class Scalar>
71 {
72  return x_space_;
73 }
74 
75 template <class Scalar>
78 {
79  return f_space_;
80 }
81 
82 template <class Scalar>
85 {
86  TEUCHOS_TEST_FOR_EXCEPTION(!isInitialized_, std::logic_error,
87  "Error, setupInOutArgs_ must be called first!\n");
88  return nominalValues_;
89 }
90 
91 template <class Scalar>
94 {
95  using Teuchos::RCP;
97  this->get_W_factory();
98  RCP<Thyra::LinearOpBase<Scalar> > matrix = this->create_W_op();
99  {
100  // 01/20/09 tscoffe: This is a total hack to provide a full rank matrix to
101  // linearOpWithSolve because it ends up factoring the matrix during
102  // initialization, which it really shouldn't do, or I'm doing something
103  // wrong here. The net effect is that I get exceptions thrown in
104  // optimized mode due to the matrix being rank deficient unless I do this.
106  Teuchos::rcp_dynamic_cast<Thyra::MultiVectorBase<Scalar> >(matrix,
107  true);
108  {
109  RCP<Thyra::VectorBase<Scalar> > vec = Thyra::createMember(x_space_);
110  {
111  Thyra::DetachedVectorView<Scalar> vec_view(*vec);
112  vec_view[0] = 1.0;
113  }
114  V_V(multivec->col(0).ptr(), *vec);
115  }
116  }
118  Thyra::linearOpWithSolve<Scalar>(*W_factory, matrix);
119  return W;
120 }
121 
122 template <class Scalar>
125 {
127  Thyra::createMembers(x_space_, dim_);
128  return (matrix);
129 }
130 
131 template <class Scalar>
134 {
136  Thyra::defaultSerialDenseLinearOpWithSolveFactory<Scalar>();
137  return W_factory;
138 }
139 
140 template <class Scalar>
143 {
144  setupInOutArgs_();
145  return inArgs_;
146 }
147 
148 // Private functions overridden from ModelEvaluatorDefaultBase
149 
150 template <class Scalar>
153 {
154  setupInOutArgs_();
155  return outArgs_;
156 }
157 
158 template <class Scalar>
161  const Thyra::ModelEvaluatorBase::OutArgs<Scalar> &outArgs) const
162 {
164  using Teuchos::RCP;
165  using Teuchos::rcp_dynamic_cast;
167  using Thyra::VectorBase;
168  TEUCHOS_TEST_FOR_EXCEPTION(!isInitialized_, std::logic_error,
169  "Error, setupInOutArgs_ must be called first!\n");
170 
171  const RCP<const VectorBase<Scalar> > x_in = inArgs.get_x().assert_not_null();
173 
174  Scalar b = b_;
175  const RCP<const VectorBase<Scalar> > p_in = inArgs.get_p(0);
176  if (p_in != Teuchos::null) {
178  b = p_in_view[0];
179  }
180 
181  RCP<const MultiVectorBase<Scalar> > DxDp_in, DxdotDp_in;
182  if (inArgs.get_p(1) != Teuchos::null)
183  DxDp_in = rcp_dynamic_cast<const DMVPV>(inArgs.get_p(1))->getMultiVector();
184  if (inArgs.get_p(2) != Teuchos::null)
185  DxdotDp_in =
186  rcp_dynamic_cast<const DMVPV>(inArgs.get_p(2))->getMultiVector();
187 
188  Scalar beta = inArgs.get_beta();
189 
190  const RCP<VectorBase<Scalar> > f_out = outArgs.get_f();
191  const RCP<Thyra::LinearOpBase<Scalar> > W_out = outArgs.get_W_op();
194  DfDp_out = DfDp.getMultiVector();
195 
196  if (inArgs.get_x_dot().is_null()) {
197  // Evaluate the Explicit ODE f(x,t) [= 0]
198  if (!is_null(f_out)) {
199  Thyra::DetachedVectorView<Scalar> f_out_view(*f_out);
200  f_out_view[0] = x_in_view[0] * x_in_view[0] - b * b;
201  }
202  if (!is_null(W_out)) {
204  Teuchos::rcp_dynamic_cast<Thyra::MultiVectorBase<Scalar> >(W_out,
205  true);
206  Thyra::DetachedMultiVectorView<Scalar> matrix_view(*matrix);
207  matrix_view(0, 0) = beta * 2.0 * x_in_view[0];
208  }
209  if (!is_null(DfDp_out)) {
210  Thyra::DetachedMultiVectorView<Scalar> DfDp_out_view(*DfDp_out);
211  DfDp_out_view(0, 0) = -2.0 * b;
212 
213  // Compute df/dp + (df/dx) * (dx/dp)
214  if (useDfDpAsTangent_ && !is_null(DxDp_in)) {
216  DfDp_out_view(0, 0) += 2.0 * x_in_view[0] * DxDp(0, 0);
217  }
218  }
219  }
220  else {
221  // Evaluate the implicit ODE f(xdot, x, t) [= 0]
222  RCP<const VectorBase<Scalar> > x_dot_in;
223  x_dot_in = inArgs.get_x_dot().assert_not_null();
224  Scalar alpha = inArgs.get_alpha();
225  if (!is_null(f_out)) {
226  Thyra::DetachedVectorView<Scalar> f_out_view(*f_out);
227  Thyra::ConstDetachedVectorView<Scalar> x_dot_in_view(*x_dot_in);
228  f_out_view[0] = x_dot_in_view[0] - (x_in_view[0] * x_in_view[0] - b * b);
229  }
230  if (!is_null(W_out)) {
232  Teuchos::rcp_dynamic_cast<Thyra::MultiVectorBase<Scalar> >(W_out,
233  true);
234  Thyra::DetachedMultiVectorView<Scalar> matrix_view(*matrix);
235  matrix_view(0, 0) = alpha - beta * 2.0 * x_in_view[0];
236  }
237  if (!is_null(DfDp_out)) {
238  Thyra::DetachedMultiVectorView<Scalar> DfDp_out_view(*DfDp_out);
239  DfDp_out_view(0, 0) = 2.0 * b;
240 
241  // Compute df/dp + (df/dx_dot) * (dx_dot/dp) + (df/dx) * (dx/dp)
242  if (useDfDpAsTangent_ && !is_null(DxdotDp_in)) {
243  Thyra::ConstDetachedMultiVectorView<Scalar> DxdotDp(*DxdotDp_in);
244  DfDp_out_view(0, 0) += DxdotDp(0, 0);
245  }
246  if (useDfDpAsTangent_ && !is_null(DxDp_in)) {
248  DfDp_out_view(0, 0) += -2.0 * x_in_view[0] * DxDp(0, 0);
249  }
250  }
251  }
252 
253  // Responses: g = x
254  RCP<VectorBase<Scalar> > g_out = outArgs.get_g(0);
255  if (g_out != Teuchos::null) Thyra::assign(g_out.ptr(), *x_in);
256 
258  outArgs.get_DgDp(0, 0).getMultiVector();
259  if (DgDp_out != Teuchos::null) Thyra::assign(DgDp_out.ptr(), Scalar(0.0));
260 
262  outArgs.get_DgDx(0).getMultiVector();
263  if (DgDx_out != Teuchos::null) {
264  Thyra::DetachedMultiVectorView<Scalar> DgDx_out_view(*DgDx_out);
265  DgDx_out_view(0, 0) = 1.0;
266  }
267 }
268 
269 template <class Scalar>
272 {
274  if (l == 0)
275  return p_space_;
276  else if (l == 1 || l == 2)
277  return DxDp_space_;
278  return Teuchos::null;
279 }
280 
281 template <class Scalar>
284 {
288  if (l == 0) {
289  p_strings->push_back("Model Coefficient: b");
290  }
291  else if (l == 1)
292  p_strings->push_back("DxDp");
293  else if (l == 2)
294  p_strings->push_back("Dx_dotDp");
295  return p_strings;
296 }
297 
298 template <class Scalar>
301 {
303  return g_space_;
304 }
305 
306 // private
307 
308 template <class Scalar>
310 {
311  if (isInitialized_) {
312  return;
313  }
314 
315  {
316  // Set up prototypical InArgs
318  inArgs.setModelEvalDescription(this->description());
319  inArgs.setSupports(Thyra::ModelEvaluatorBase::IN_ARG_t);
320  inArgs.setSupports(Thyra::ModelEvaluatorBase::IN_ARG_x);
321  inArgs.setSupports(Thyra::ModelEvaluatorBase::IN_ARG_beta);
322  inArgs.setSupports(Thyra::ModelEvaluatorBase::IN_ARG_x_dot);
323  inArgs.setSupports(Thyra::ModelEvaluatorBase::IN_ARG_alpha);
324  inArgs.set_Np(Np_);
325  inArgs_ = inArgs;
326  }
327 
328  {
329  // Set up prototypical OutArgs
331  outArgs.setModelEvalDescription(this->description());
332  outArgs.setSupports(Thyra::ModelEvaluatorBase::OUT_ARG_f);
333  outArgs.setSupports(Thyra::ModelEvaluatorBase::OUT_ARG_W_op);
334  outArgs.set_Np_Ng(Np_, Ng_);
335  outArgs.setSupports(Thyra::ModelEvaluatorBase::OUT_ARG_DfDp, 0,
336  Thyra::ModelEvaluatorBase::DERIV_MV_JACOBIAN_FORM);
337  outArgs.setSupports(Thyra::ModelEvaluatorBase::OUT_ARG_DgDx, 0,
338  Thyra::ModelEvaluatorBase::DERIV_MV_GRADIENT_FORM);
339  outArgs.setSupports(Thyra::ModelEvaluatorBase::OUT_ARG_DgDp, 0, 0,
340  Thyra::ModelEvaluatorBase::DERIV_MV_GRADIENT_FORM);
341  outArgs_ = outArgs;
342  }
343 
344  // Set up nominal values
345  nominalValues_ = inArgs_;
346  nominalValues_.set_t(0.0);
347  const Teuchos::RCP<Thyra::VectorBase<Scalar> > x_ic = createMember(x_space_);
348  { // scope to delete DetachedVectorView
349  Thyra::DetachedVectorView<Scalar> x_ic_view(*x_ic);
350  x_ic_view[0] = 0.0;
351  }
352  nominalValues_.set_x(x_ic);
353  const Teuchos::RCP<Thyra::VectorBase<Scalar> > p_ic = createMember(p_space_);
354  {
355  Thyra::DetachedVectorView<Scalar> p_ic_view(*p_ic);
356  p_ic_view[0] = b_;
357  }
358  nominalValues_.set_p(0, p_ic);
359 
360  const Teuchos::RCP<Thyra::VectorBase<Scalar> > x_dot_ic =
361  createMember(x_space_);
362  { // scope to delete DetachedVectorView
363  Thyra::DetachedVectorView<Scalar> x_dot_ic_view(*x_dot_ic);
364  x_dot_ic_view[0] = 0.0;
365  }
366  nominalValues_.set_x_dot(x_dot_ic);
367 
368  isInitialized_ = true;
369 }
370 
371 template <class Scalar>
373  Teuchos::RCP<Teuchos::ParameterList> const &paramList)
374 {
375  using Teuchos::get;
378  Teuchos::rcp(new ParameterList("SteadyQuadraticModel"));
379  if (paramList != Teuchos::null) tmpPL = paramList;
380  tmpPL->validateParametersAndSetDefaults(*this->getValidParameters());
381  this->setMyParamList(tmpPL);
382  Teuchos::RCP<ParameterList> pl = this->getMyNonconstParamList();
383  useDfDpAsTangent_ = get<bool>(*pl, "Use DfDp as Tangent");
384  b_ = get<Scalar>(*pl, "Coeff b");
385  isInitialized_ = false; // For setup of new in/out args
386  setupInOutArgs_();
387 }
388 
389 template <class Scalar>
392 {
394  if (is_null(validPL)) {
395  Teuchos::RCP<Teuchos::ParameterList> pl = Teuchos::parameterList();
396  pl->set("Use DfDp as Tangent", false);
397  Teuchos::setDoubleParameter("Coeff b", 1.0, "Coefficient b in model", &*pl);
398  validPL = pl;
399  }
400  return validPL;
401 }
402 
403 } // namespace Tempus_Test
404 #endif // TEMPUS_TEST_STEADY_QUADRATIC_MODEL_IMPL_HPP
Teuchos::RCP< Thyra::LinearOpBase< Scalar > > create_W_op() const
void evalModelImpl(const Thyra::ModelEvaluatorBase::InArgs< Scalar > &inArgs_bar, const Thyra::ModelEvaluatorBase::OutArgs< Scalar > &outArgs_bar) const
bool is_null(const boost::shared_ptr< T > &p)
Derivative< Scalar > get_DfDp(int l) const
RCP< const VectorBase< Scalar > > get_x_dot() const
void setModelEvalDescription(const std::string &modelEvalDescription)
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)
Teuchos::RCP< const Thyra::VectorSpaceBase< Scalar > > get_p_space(int l) const
Evaluation< VectorBase< Scalar > > get_g(int j) const
Evaluation< VectorBase< Scalar > > get_f() const
Thyra::ModelEvaluatorBase::InArgs< Scalar > createInArgs() const
Derivative< Scalar > get_DgDp(int j, int l) 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::RCP< const Teuchos::Array< std::string > > get_p_names(int l) const
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
void setParameterList(Teuchos::RCP< Teuchos::ParameterList > const &paramList)
Ptr< T > ptr() const
Teuchos::RCP< Thyra::LinearOpWithSolveBase< Scalar > > create_W() const
void validateParametersAndSetDefaults(ParameterList const &validParamList, int const depth=1000)
Thyra::ModelEvaluatorBase::InArgs< Scalar > getNominalValues() const
void setSupports(EOutArgsMembers arg, bool supports=true)
Teuchos::RCP< const Teuchos::ParameterList > getValidParameters() const
Teuchos::RCP< const Thyra::VectorSpaceBase< Scalar > > get_f_space() const
Teuchos::RCP< const Thyra::VectorSpaceBase< Scalar > > get_x_space() const
Teuchos::RCP< const Thyra::LinearOpWithSolveFactoryBase< Scalar > > get_W_factory() const
void setModelEvalDescription(const std::string &modelEvalDescription)
Thyra::ModelEvaluatorBase::OutArgs< Scalar > createOutArgsImpl() const
RCP< LinearOpBase< Scalar > > get_W_op() const
RCP< const VectorBase< Scalar > > get_x() const
SteadyQuadraticModel(Teuchos::RCP< Teuchos::ParameterList > pList=Teuchos::null)
RCP< const VectorBase< Scalar > > get_p(int l) const
RCP< MultiVectorBase< Scalar > > getMultiVector() const
Derivative< Scalar > get_DgDx(int j) const