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 
27 namespace Tempus_Test {
28 
29 template<class Scalar>
32 {
33  isInitialized_ = false;
34  dim_ = 1;
35  Np_ = 3; // Number of parameter vectors (p, dx/dp, dx_dot/dp)
36  np_ = 1; // Number of parameters in this vector (3)
37  Ng_ = 1; // Number of observation functions (1)
38  ng_ = 1; // Number of elements in this observation function (1)
39  useDfDpAsTangent_ = false;
40  b_ = 1.0;
41 
42  // Create x_space and f_space
43  x_space_ = Thyra::defaultSpmdVectorSpace<Scalar>(dim_);
44  f_space_ = Thyra::defaultSpmdVectorSpace<Scalar>(dim_);
45  // Create p_space and g_space
46  p_space_ = Thyra::defaultSpmdVectorSpace<Scalar>(np_);
47  g_space_ = Thyra::defaultSpmdVectorSpace<Scalar>(ng_);
48 
49  setParameterList(pList_);
50 
51  // Create DxDp product space
52  DxDp_space_ = Thyra::multiVectorProductVectorSpace(x_space_, np_);
53 }
54 
55 template<class Scalar>
56 Scalar
59 {
60  if (b_ < 0.0)
61  return b_;
62  return -b_;
63 }
64 
65 template<class Scalar>
66 Scalar
69 {
70  if (b_ < 0.0)
71  return 1.0;
72  return -1.0;
73 }
74 
75 template<class Scalar>
78 get_x_space() const
79 {
80  return x_space_;
81 }
82 
83 
84 template<class Scalar>
87 get_f_space() const
88 {
89  return f_space_;
90 }
91 
92 
93 template<class Scalar>
97 {
98  TEUCHOS_TEST_FOR_EXCEPTION( !isInitialized_, std::logic_error,
99  "Error, setupInOutArgs_ must be called first!\n");
100  return nominalValues_;
101 }
102 
103 
104 template<class Scalar>
107 create_W() const
108 {
109  using Teuchos::RCP;
110  RCP<const Thyra::LinearOpWithSolveFactoryBase<Scalar> > W_factory = this->get_W_factory();
111  RCP<Thyra::LinearOpBase<Scalar> > matrix = this->create_W_op();
112  {
113  // 01/20/09 tscoffe: This is a total hack to provide a full rank matrix to
114  // linearOpWithSolve because it ends up factoring the matrix during
115  // initialization, which it really shouldn't do, or I'm doing something
116  // wrong here. The net effect is that I get exceptions thrown in
117  // optimized mode due to the matrix being rank deficient unless I do this.
119  Teuchos::rcp_dynamic_cast<Thyra::MultiVectorBase<Scalar> >(matrix,true);
120  {
121  RCP<Thyra::VectorBase<Scalar> > vec = Thyra::createMember(x_space_);
122  {
123  Thyra::DetachedVectorView<Scalar> vec_view( *vec );
124  vec_view[0] = 1.0;
125  }
126  V_V(multivec->col(0).ptr(),*vec);
127  }
128  }
130  Thyra::linearOpWithSolve<Scalar>(
131  *W_factory,
132  matrix
133  );
134  return W;
135 }
136 
137 template<class Scalar>
140 create_W_op() const
141 {
143  Thyra::createMembers(x_space_, dim_);
144  return(matrix);
145 }
146 
147 
148 template<class Scalar>
152 {
154  Thyra::defaultSerialDenseLinearOpWithSolveFactory<Scalar>();
155  return W_factory;
156 }
157 
158 
159 template<class Scalar>
163 {
164  setupInOutArgs_();
165  return inArgs_;
166 }
167 
168 
169 // Private functions overridden from ModelEvaluatorDefaultBase
170 
171 
172 template<class Scalar>
176 {
177  setupInOutArgs_();
178  return outArgs_;
179 }
180 
181 
182 template<class Scalar>
183 void
188  ) const
189 {
191  using Teuchos::RCP;
192  using Thyra::VectorBase;
194  using Teuchos::rcp_dynamic_cast;
195  TEUCHOS_TEST_FOR_EXCEPTION( !isInitialized_, std::logic_error,
196  "Error, setupInOutArgs_ must be called first!\n");
197 
198  const RCP<const VectorBase<Scalar> > x_in = inArgs.get_x().assert_not_null();
199  Thyra::ConstDetachedVectorView<Scalar> x_in_view( *x_in );
200 
201  Scalar b = b_;
202  const RCP<const VectorBase<Scalar> > p_in = inArgs.get_p(0);
203  if (p_in != Teuchos::null) {
204  Thyra::ConstDetachedVectorView<Scalar> p_in_view( *p_in );
205  b = p_in_view[0];
206  }
207 
208  RCP<const MultiVectorBase<Scalar> > DxDp_in, DxdotDp_in;
209  if (inArgs.get_p(1) != Teuchos::null)
210  DxDp_in =
211  rcp_dynamic_cast<const DMVPV>(inArgs.get_p(1))->getMultiVector();
212  if (inArgs.get_p(2) != Teuchos::null)
213  DxdotDp_in =
214  rcp_dynamic_cast<const DMVPV>(inArgs.get_p(2))->getMultiVector();
215 
216  Scalar beta = inArgs.get_beta();
217 
218  const RCP<VectorBase<Scalar> > f_out = outArgs.get_f();
219  const RCP<Thyra::LinearOpBase<Scalar> > W_out = outArgs.get_W_op();
222  DfDp_out = DfDp.getMultiVector();
223 
224  if (inArgs.get_x_dot().is_null()) {
225 
226  // Evaluate the Explicit ODE f(x,t) [= 0]
227  if (!is_null(f_out)) {
228  Thyra::DetachedVectorView<Scalar> f_out_view( *f_out );
229  f_out_view[0] = x_in_view[0]*x_in_view[0] - b*b;
230  }
231  if (!is_null(W_out)) {
233  Teuchos::rcp_dynamic_cast<Thyra::MultiVectorBase<Scalar> >(W_out,true);
234  Thyra::DetachedMultiVectorView<Scalar> matrix_view( *matrix );
235  matrix_view(0,0) = 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) * (dx/dp)
242  if (useDfDpAsTangent_ && !is_null(DxDp_in)) {
244  DfDp_out_view(0,0) += 2.0*x_in_view[0]*DxDp(0,0);
245  }
246  }
247  } else {
248 
249  // Evaluate the implicit ODE f(xdot, x, t) [= 0]
250  RCP<const VectorBase<Scalar> > x_dot_in;
251  x_dot_in = inArgs.get_x_dot().assert_not_null();
252  Scalar alpha = inArgs.get_alpha();
253  if (!is_null(f_out)) {
254  Thyra::DetachedVectorView<Scalar> f_out_view( *f_out );
255  Thyra::ConstDetachedVectorView<Scalar> x_dot_in_view( *x_dot_in );
256  f_out_view[0] = x_dot_in_view[0] - (x_in_view[0]*x_in_view[0] - b*b);
257  }
258  if (!is_null(W_out)) {
260  Teuchos::rcp_dynamic_cast<Thyra::MultiVectorBase<Scalar> >(W_out,true);
261  Thyra::DetachedMultiVectorView<Scalar> matrix_view( *matrix );
262  matrix_view(0,0) = alpha - beta*2.0*x_in_view[0];
263  }
264  if (!is_null(DfDp_out)) {
265  Thyra::DetachedMultiVectorView<Scalar> DfDp_out_view( *DfDp_out );
266  DfDp_out_view(0,0) = 2.0*b;
267 
268  // Compute df/dp + (df/dx_dot) * (dx_dot/dp) + (df/dx) * (dx/dp)
269  if (useDfDpAsTangent_ && !is_null(DxdotDp_in)) {
270  Thyra::ConstDetachedMultiVectorView<Scalar> DxdotDp( *DxdotDp_in );
271  DfDp_out_view(0,0) += DxdotDp(0,0);
272  }
273  if (useDfDpAsTangent_ && !is_null(DxDp_in)) {
275  DfDp_out_view(0,0) += -2.0*x_in_view[0]*DxDp(0,0);
276  }
277  }
278  }
279 
280  // Responses: g = x
281  RCP<VectorBase<Scalar> > g_out = outArgs.get_g(0);
282  if (g_out != Teuchos::null)
283  Thyra::assign(g_out.ptr(), *x_in);
284 
286  outArgs.get_DgDp(0,0).getMultiVector();
287  if (DgDp_out != Teuchos::null)
288  Thyra::assign(DgDp_out.ptr(), Scalar(0.0));
289 
291  outArgs.get_DgDx(0).getMultiVector();
292  if (DgDx_out != Teuchos::null) {
293  Thyra::DetachedMultiVectorView<Scalar> DgDx_out_view( *DgDx_out );
294  DgDx_out_view(0,0) = 1.0;
295  }
296 }
297 
298 template<class Scalar>
301 get_p_space(int l) const
302 {
304  if (l == 0)
305  return p_space_;
306  else if (l == 1 || l == 2)
307  return DxDp_space_;
308  return Teuchos::null;
309 }
310 
311 template<class Scalar>
314 get_p_names(int l) const
315 {
319  if (l == 0) {
320  p_strings->push_back("Model Coefficient: b");
321  }
322  else if (l == 1)
323  p_strings->push_back("DxDp");
324  else if (l == 2)
325  p_strings->push_back("Dx_dotDp");
326  return p_strings;
327 }
328 
329 template<class Scalar>
332 get_g_space(int j) const
333 {
335  return g_space_;
336 }
337 
338 // private
339 
340 template<class Scalar>
341 void
344 {
345  if (isInitialized_) {
346  return;
347  }
348 
349  {
350  // Set up prototypical InArgs
352  inArgs.setModelEvalDescription(this->description());
353  inArgs.setSupports( Thyra::ModelEvaluatorBase::IN_ARG_t );
354  inArgs.setSupports( Thyra::ModelEvaluatorBase::IN_ARG_x );
355  inArgs.setSupports( Thyra::ModelEvaluatorBase::IN_ARG_beta );
356  inArgs.setSupports( Thyra::ModelEvaluatorBase::IN_ARG_x_dot );
357  inArgs.setSupports( Thyra::ModelEvaluatorBase::IN_ARG_alpha );
358  inArgs.set_Np(Np_);
359  inArgs_ = inArgs;
360  }
361 
362  {
363  // Set up prototypical OutArgs
365  outArgs.setModelEvalDescription(this->description());
366  outArgs.setSupports( Thyra::ModelEvaluatorBase::OUT_ARG_f );
367  outArgs.setSupports( Thyra::ModelEvaluatorBase::OUT_ARG_W_op );
368  outArgs.set_Np_Ng(Np_,Ng_);
369  outArgs.setSupports( Thyra::ModelEvaluatorBase::OUT_ARG_DfDp,0,
370  Thyra::ModelEvaluatorBase::DERIV_MV_JACOBIAN_FORM );
371  outArgs.setSupports( Thyra::ModelEvaluatorBase::OUT_ARG_DgDx,0,
372  Thyra::ModelEvaluatorBase::DERIV_MV_GRADIENT_FORM );
373  outArgs.setSupports( Thyra::ModelEvaluatorBase::OUT_ARG_DgDp,0,0,
374  Thyra::ModelEvaluatorBase::DERIV_MV_GRADIENT_FORM );
375  outArgs_ = outArgs;
376  }
377 
378  // Set up nominal values
379  nominalValues_ = inArgs_;
380  nominalValues_.set_t(0.0);
381  const Teuchos::RCP<Thyra::VectorBase<Scalar> > x_ic = createMember(x_space_);
382  { // scope to delete DetachedVectorView
383  Thyra::DetachedVectorView<Scalar> x_ic_view( *x_ic );
384  x_ic_view[0] = 0.0;
385  }
386  nominalValues_.set_x(x_ic);
387  const Teuchos::RCP<Thyra::VectorBase<Scalar> > p_ic = createMember(p_space_);
388  {
389  Thyra::DetachedVectorView<Scalar> p_ic_view( *p_ic );
390  p_ic_view[0] = b_;
391  }
392  nominalValues_.set_p(0,p_ic);
393 
394  const Teuchos::RCP<Thyra::VectorBase<Scalar> > x_dot_ic =
395  createMember(x_space_);
396  { // scope to delete DetachedVectorView
397  Thyra::DetachedVectorView<Scalar> x_dot_ic_view( *x_dot_ic );
398  x_dot_ic_view[0] = 0.0;
399  }
400  nominalValues_.set_x_dot(x_dot_ic);
401 
402  isInitialized_ = true;
403 
404 }
405 
406 template<class Scalar>
407 void
410 {
411  using Teuchos::get;
414  Teuchos::rcp(new ParameterList("SteadyQuadraticModel"));
415  if (paramList != Teuchos::null) tmpPL = paramList;
416  tmpPL->validateParametersAndSetDefaults(*this->getValidParameters());
417  this->setMyParamList(tmpPL);
418  Teuchos::RCP<ParameterList> pl = this->getMyNonconstParamList();
419  useDfDpAsTangent_ = get<bool>(*pl, "Use DfDp as Tangent");
420  b_ = get<Scalar>(*pl,"Coeff b");
421  isInitialized_ = false; // For setup of new in/out args
422  setupInOutArgs_();
423 }
424 
425 template<class Scalar>
429 {
431  if (is_null(validPL)) {
432  Teuchos::RCP<Teuchos::ParameterList> pl = Teuchos::parameterList();
433  pl->set("Use DfDp as Tangent", false);
434  Teuchos::setDoubleParameter(
435  "Coeff b", 1.0, "Coefficient b in model", &*pl);
436  validPL = pl;
437  }
438  return validPL;
439 }
440 
441 } // namespace Tempus_Test
442 #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