Tempus  Version of the Day
Time Integration
SinCosModel_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_SINCOS_MODEL_IMPL_HPP
10 #define TEMPUS_TEST_SINCOS_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_VectorStdOps.hpp"
21 #include "Thyra_MultiVectorStdOps.hpp"
22 #include "Thyra_DefaultMultiVectorProductVector.hpp"
23 
24 #include <iostream>
25 
26 
27 namespace Tempus_Test {
28 
29 template<class Scalar>
30 SinCosModel<Scalar>::
31 SinCosModel(Teuchos::RCP<Teuchos::ParameterList> pList_)
32 {
33  isInitialized_ = false;
34  dim_ = 2;
35  Np_ = 3; // Number of parameter vectors (p, dx/dp, dx_dot/dp)
36  np_ = 3; // Number of parameters in this vector (3)
37  Ng_ = 1; // Number of observation functions (1)
38  ng_ = dim_; // Number of elements in this observation function ( == x )
39  acceptModelParams_ = false;
40  useDfDpAsTangent_ = false;
41  haveIC_ = true;
42  a_ = 0.0;
43  f_ = 1.0;
44  L_ = 1.0;
45  x0_ic_ = 0.0;
46  x1_ic_ = 1.0;
47  t0_ic_ = 0.0;
48 
49  // Create x_space and f_space
50  x_space_ = Thyra::defaultSpmdVectorSpace<Scalar>(dim_);
51  f_space_ = Thyra::defaultSpmdVectorSpace<Scalar>(dim_);
52  // Create p_space and g_space
53  p_space_ = Thyra::defaultSpmdVectorSpace<Scalar>(np_);
54  g_space_ = Thyra::defaultSpmdVectorSpace<Scalar>(ng_);
55 
56  setParameterList(pList_);
57 
58  // Create DxDp product space
59  DxDp_space_ = Thyra::multiVectorProductVectorSpace(x_space_, np_);
60 }
61 
62 template<class Scalar>
63 Thyra::ModelEvaluatorBase::InArgs<Scalar>
64 SinCosModel<Scalar>::
65 getExactSolution(double t) const
66 {
67  TEUCHOS_TEST_FOR_EXCEPTION( !isInitialized_, std::logic_error,
68  "Error, setupInOutArgs_ must be called first!\n");
69  Thyra::ModelEvaluatorBase::InArgs<Scalar> inArgs = inArgs_;
70  double exact_t = t;
71  inArgs.set_t(exact_t);
72  Teuchos::RCP<Thyra::VectorBase<Scalar> > exact_x = createMember(x_space_);
73  { // scope to delete DetachedVectorView
74  Thyra::DetachedVectorView<Scalar> exact_x_view(*exact_x);
75  exact_x_view[0] = a_+b_*sin((f_/L_)*t+phi_);
76  exact_x_view[1] = b_*(f_/L_)*cos((f_/L_)*t+phi_);
77  }
78  inArgs.set_x(exact_x);
79  Teuchos::RCP<Thyra::VectorBase<Scalar> > exact_x_dot = createMember(x_space_);
80  { // scope to delete DetachedVectorView
81  Thyra::DetachedVectorView<Scalar> exact_x_dot_view(*exact_x_dot);
82  exact_x_dot_view[0] = b_*(f_/L_)*cos((f_/L_)*t+phi_);
83  exact_x_dot_view[1] = -b_*(f_/L_)*(f_/L_)*sin((f_/L_)*t+phi_);
84  }
85  inArgs.set_x_dot(exact_x_dot);
86  return(inArgs);
87 }
88 
89 //
90 // 06/24/09 tscoffe:
91 // These are the exact sensitivities for the problem assuming the initial conditions are specified as:
92 // s(0) = [1;0] s(1) = [0;b/L] s(2) = [0;-b*f/(L*L)]
93 // sdot(0) = [0;0] sdot(1) = [0;-3*f*f*b/(L*L*L)] sdot(2) = [0;3*f*f*f*b/(L*L*L*L)]
94 //
95 template<class Scalar>
96 Thyra::ModelEvaluatorBase::InArgs<Scalar>
97 SinCosModel<Scalar>::
98 getExactSensSolution(int j, double t) const
99 {
100  TEUCHOS_TEST_FOR_EXCEPTION( !isInitialized_, std::logic_error,
101  "Error, setupInOutArgs_ must be called first!\n");
102  Thyra::ModelEvaluatorBase::InArgs<Scalar> inArgs = inArgs_;
103  if (!acceptModelParams_) {
104  return inArgs;
105  }
106  TEUCHOS_ASSERT_IN_RANGE_UPPER_EXCLUSIVE( j, 0, np_ );
107  double exact_t = t;
108  Scalar b = b_;
109  Scalar f = f_;
110  Scalar L = L_;
111  Scalar phi = phi_;
112  inArgs.set_t(exact_t);
113  Teuchos::RCP<Thyra::VectorBase<Scalar> > exact_s = createMember(x_space_);
114  { // scope to delete DetachedVectorView
115  Thyra::DetachedVectorView<Scalar> exact_s_view(*exact_s);
116  if (j == 0) { // dx/da
117  exact_s_view[0] = 1.0;
118  exact_s_view[1] = 0.0;
119  } else if (j == 1) { // dx/df
120  exact_s_view[0] = (b/L)*t*cos((f/L)*t+phi);
121  exact_s_view[1] = (b/L)*cos((f/L)*t+phi)-(b*f*t/(L*L))*sin((f/L)*t+phi);
122  } else { // dx/dL
123  exact_s_view[0] = -(b*f*t/(L*L))*cos((f/L)*t+phi);
124  exact_s_view[1] = -(b*f/(L*L))*cos((f/L)*t+phi)+(b*f*f*t/(L*L*L))*sin((f/L)*t+phi);
125  }
126  }
127  inArgs.set_x(exact_s);
128  Teuchos::RCP<Thyra::VectorBase<Scalar> > exact_s_dot = createMember(x_space_);
129  { // scope to delete DetachedVectorView
130  Thyra::DetachedVectorView<Scalar> exact_s_dot_view(*exact_s_dot);
131  if (j == 0) { // dxdot/da
132  exact_s_dot_view[0] = 0.0;
133  exact_s_dot_view[1] = 0.0;
134  } else if (j == 1) { // dxdot/df
135  exact_s_dot_view[0] = (b/L)*cos((f/L)*t+phi)-(b*f*t/(L*L))*sin((f/L)*t+phi);
136  exact_s_dot_view[1] = -(2.0*b*f/(L*L))*sin((f/L)*t+phi)-(b*f*f*t/(L*L*L))*cos((f/L)*t+phi);
137  } else { // dxdot/dL
138  exact_s_dot_view[0] = -(b*f/(L*L))*cos((f/L)*t+phi)+(b*f*f*t/(L*L*L))*sin((f/L)*t+phi);
139  exact_s_dot_view[1] = (2.0*b*f*f/(L*L*L))*sin((f/L)*t+phi)+(b*f*f*f*t/(L*L*L*L))*cos((f/L)*t+phi);
140  }
141  }
142  inArgs.set_x_dot(exact_s_dot);
143  return(inArgs);
144 }
145 
146 template<class Scalar>
147 Teuchos::RCP<const Thyra::VectorSpaceBase<Scalar> >
148 SinCosModel<Scalar>::
149 get_x_space() const
150 {
151  return x_space_;
152 }
153 
154 
155 template<class Scalar>
156 Teuchos::RCP<const Thyra::VectorSpaceBase<Scalar> >
157 SinCosModel<Scalar>::
158 get_f_space() const
159 {
160  return f_space_;
161 }
162 
163 
164 template<class Scalar>
165 Thyra::ModelEvaluatorBase::InArgs<Scalar>
166 SinCosModel<Scalar>::
167 getNominalValues() const
168 {
169  TEUCHOS_TEST_FOR_EXCEPTION( !isInitialized_, std::logic_error,
170  "Error, setupInOutArgs_ must be called first!\n");
171  return nominalValues_;
172 }
173 
174 
175 template<class Scalar>
176 Teuchos::RCP<Thyra::LinearOpWithSolveBase<Scalar> >
177 SinCosModel<Scalar>::
178 create_W() const
179 {
180  using Teuchos::RCP;
181  RCP<const Thyra::LinearOpWithSolveFactoryBase<Scalar> > W_factory = this->get_W_factory();
182  RCP<Thyra::LinearOpBase<Scalar> > matrix = this->create_W_op();
183  {
184  // 01/20/09 tscoffe: This is a total hack to provide a full rank matrix to
185  // linearOpWithSolve because it ends up factoring the matrix during
186  // initialization, which it really shouldn't do, or I'm doing something
187  // wrong here. The net effect is that I get exceptions thrown in
188  // optimized mode due to the matrix being rank deficient unless I do this.
189  RCP<Thyra::MultiVectorBase<Scalar> > multivec = Teuchos::rcp_dynamic_cast<Thyra::MultiVectorBase<Scalar> >(matrix,true);
190  {
191  RCP<Thyra::VectorBase<Scalar> > vec = Thyra::createMember(x_space_);
192  {
193  Thyra::DetachedVectorView<Scalar> vec_view( *vec );
194  vec_view[0] = 0.0;
195  vec_view[1] = 1.0;
196  }
197  V_V(multivec->col(0).ptr(),*vec);
198  {
199  Thyra::DetachedVectorView<Scalar> vec_view( *vec );
200  vec_view[0] = 1.0;
201  vec_view[1] = 0.0;
202  }
203  V_V(multivec->col(1).ptr(),*vec);
204  }
205  }
206  RCP<Thyra::LinearOpWithSolveBase<Scalar> > W =
207  Thyra::linearOpWithSolve<Scalar>(
208  *W_factory,
209  matrix
210  );
211  return W;
212 }
213 //Teuchos::RCP<Thyra::LinearOpWithSolveBase<Scalar> >
214 //SinCosModel<Scalar>::
215 //create_W() const
216 //{
217 // return Thyra::multiVectorLinearOpWithSolve<Scalar>();
218 //}
219 
220 
221 template<class Scalar>
222 Teuchos::RCP<Thyra::LinearOpBase<Scalar> >
223 SinCosModel<Scalar>::
224 create_W_op() const
225 {
226  Teuchos::RCP<Thyra::MultiVectorBase<Scalar> > matrix = Thyra::createMembers(x_space_, dim_);
227  return(matrix);
228 }
229 
230 
231 template<class Scalar>
232 Teuchos::RCP<const Thyra::LinearOpWithSolveFactoryBase<Scalar> >
233 SinCosModel<Scalar>::
234 get_W_factory() const
235 {
236  Teuchos::RCP<Thyra::LinearOpWithSolveFactoryBase<Scalar> > W_factory =
237  Thyra::defaultSerialDenseLinearOpWithSolveFactory<Scalar>();
238  return W_factory;
239 }
240 // Teuchos::RCP<const Thyra::LinearOpBase<Scalar> > fwdOp = this->create_W_op();
241 // Teuchos::RCP<Thyra::LinearOpWithSolveBase<Scalar> > W =
242 // Thyra::linearOpWithSolve<Scalar>(
243 // *W_factory,
244 // fwdOp
245 // );
246 // W_factory->initializeOp(
247 // Thyra::defaultLinearOpSource<Scalar>(fwdOp),
248 // &*W,
249 // Thyra::SUPPORT_SOLVE_UNSPECIFIED
250 // );
251 
252 
253 template<class Scalar>
254 Thyra::ModelEvaluatorBase::InArgs<Scalar>
255 SinCosModel<Scalar>::
256 createInArgs() const
257 {
258  setupInOutArgs_();
259  return inArgs_;
260 }
261 
262 
263 // Private functions overridden from ModelEvaluatorDefaultBase
264 
265 
266 template<class Scalar>
267 Thyra::ModelEvaluatorBase::OutArgs<Scalar>
268 SinCosModel<Scalar>::
269 createOutArgsImpl() const
270 {
271  setupInOutArgs_();
272  return outArgs_;
273 }
274 
275 
276 template<class Scalar>
277 void
278 SinCosModel<Scalar>::
279 evalModelImpl(
280  const Thyra::ModelEvaluatorBase::InArgs<Scalar> &inArgs,
281  const Thyra::ModelEvaluatorBase::OutArgs<Scalar> &outArgs
282  ) const
283 {
284  typedef Thyra::DefaultMultiVectorProductVector<Scalar> DMVPV;
285  using Teuchos::RCP;
286  using Thyra::VectorBase;
287  using Thyra::MultiVectorBase;
288  using Teuchos::rcp_dynamic_cast;
289  TEUCHOS_TEST_FOR_EXCEPTION( !isInitialized_, std::logic_error,
290  "Error, setupInOutArgs_ must be called first!\n");
291 
292  const RCP<const VectorBase<Scalar> > x_in = inArgs.get_x().assert_not_null();
293  Thyra::ConstDetachedVectorView<Scalar> x_in_view( *x_in );
294 
295  //double t = inArgs.get_t();
296  Scalar a = a_;
297  Scalar f = f_;
298  Scalar L = L_;
299  if (acceptModelParams_) {
300  const RCP<const VectorBase<Scalar> > p_in =
301  inArgs.get_p(0).assert_not_null();
302  Thyra::ConstDetachedVectorView<Scalar> p_in_view( *p_in );
303  a = p_in_view[0];
304  f = p_in_view[1];
305  L = p_in_view[2];
306  }
307 
308  RCP<const MultiVectorBase<Scalar> > DxDp_in, DxdotDp_in;
309  if (acceptModelParams_) {
310  if (inArgs.get_p(1) != Teuchos::null)
311  DxDp_in =
312  rcp_dynamic_cast<const DMVPV>(inArgs.get_p(1))->getMultiVector();
313  if (inArgs.get_p(2) != Teuchos::null)
314  DxdotDp_in =
315  rcp_dynamic_cast<const DMVPV>(inArgs.get_p(2))->getMultiVector();
316  }
317 
318  Scalar beta = inArgs.get_beta();
319 
320  const RCP<VectorBase<Scalar> > f_out = outArgs.get_f();
321  const RCP<Thyra::LinearOpBase<Scalar> > W_out = outArgs.get_W_op();
322  RCP<Thyra::MultiVectorBase<Scalar> > DfDp_out;
323  if (acceptModelParams_) {
324  Thyra::ModelEvaluatorBase::Derivative<Scalar> DfDp = outArgs.get_DfDp(0);
325  DfDp_out = DfDp.getMultiVector();
326  }
327 
328  if (inArgs.get_x_dot().is_null()) {
329 
330  // Evaluate the Explicit ODE f(x,t) [= 0]
331  if (!is_null(f_out)) {
332  Thyra::DetachedVectorView<Scalar> f_out_view( *f_out );
333  f_out_view[0] = x_in_view[1];
334  f_out_view[1] = (f/L)*(f/L)*(a-x_in_view[0]);
335  }
336  if (!is_null(W_out)) {
337  RCP<Thyra::MultiVectorBase<Scalar> > matrix =
338  Teuchos::rcp_dynamic_cast<Thyra::MultiVectorBase<Scalar> >(W_out,true);
339  Thyra::DetachedMultiVectorView<Scalar> matrix_view( *matrix );
340  matrix_view(0,0) = 0.0; // d(f0)/d(x0_n)
341  matrix_view(0,1) = +beta; // d(f0)/d(x1_n)
342  matrix_view(1,0) = -beta*(f/L)*(f/L); // d(f1)/d(x0_n)
343  matrix_view(1,1) = 0.0; // d(f1)/d(x1_n)
344  // Note: alpha = d(xdot)/d(x_n) and beta = d(x)/d(x_n)
345  }
346  if (!is_null(DfDp_out)) {
347  Thyra::DetachedMultiVectorView<Scalar> DfDp_out_view( *DfDp_out );
348  DfDp_out_view(0,0) = 0.0;
349  DfDp_out_view(0,1) = 0.0;
350  DfDp_out_view(0,2) = 0.0;
351  DfDp_out_view(1,0) = (f/L)*(f/L);
352  DfDp_out_view(1,1) = (2.0*f/(L*L))*(a-x_in_view[0]);
353  DfDp_out_view(1,2) = -(2.0*f*f/(L*L*L))*(a-x_in_view[0]);
354 
355  // Compute df/dp + (df/dx) * (dx/dp)
356  if (useDfDpAsTangent_ && !is_null(DxDp_in)) {
357  Thyra::ConstDetachedMultiVectorView<Scalar> DxDp( *DxDp_in );
358  DfDp_out_view(0,0) += DxDp(1,0);
359  DfDp_out_view(0,1) += DxDp(1,1);
360  DfDp_out_view(0,2) += DxDp(1,2);
361  DfDp_out_view(1,0) += -(f/L)*(f/L) * DxDp(0,0);
362  DfDp_out_view(1,1) += -(f/L)*(f/L) * DxDp(0,1);
363  DfDp_out_view(1,2) += -(f/L)*(f/L) * DxDp(0,2);
364  }
365  }
366  } else {
367 
368  // Evaluate the implicit ODE f(xdot, x, t) [= 0]
369  RCP<const VectorBase<Scalar> > x_dot_in;
370  x_dot_in = inArgs.get_x_dot().assert_not_null();
371  Scalar alpha = inArgs.get_alpha();
372  if (!is_null(f_out)) {
373  Thyra::DetachedVectorView<Scalar> f_out_view( *f_out );
374  Thyra::ConstDetachedVectorView<Scalar> x_dot_in_view( *x_dot_in );
375  f_out_view[0] = x_dot_in_view[0] - x_in_view[1];
376  f_out_view[1] = x_dot_in_view[1] - (f/L)*(f/L)*(a-x_in_view[0]);
377  }
378  if (!is_null(W_out)) {
379  RCP<Thyra::MultiVectorBase<Scalar> > matrix =
380  Teuchos::rcp_dynamic_cast<Thyra::MultiVectorBase<Scalar> >(W_out,true);
381  Thyra::DetachedMultiVectorView<Scalar> matrix_view( *matrix );
382  matrix_view(0,0) = alpha; // d(f0)/d(x0_n)
383  matrix_view(0,1) = -beta; // d(f0)/d(x1_n)
384  matrix_view(1,0) = +beta*(f/L)*(f/L); // d(f1)/d(x0_n)
385  matrix_view(1,1) = alpha; // d(f1)/d(x1_n)
386  // Note: alpha = d(xdot)/d(x_n) and beta = d(x)/d(x_n)
387  }
388  if (!is_null(DfDp_out)) {
389  Thyra::DetachedMultiVectorView<Scalar> DfDp_out_view( *DfDp_out );
390  DfDp_out_view(0,0) = 0.0;
391  DfDp_out_view(0,1) = 0.0;
392  DfDp_out_view(0,2) = 0.0;
393  DfDp_out_view(1,0) = -(f/L)*(f/L);
394  DfDp_out_view(1,1) = -(2.0*f/(L*L))*(a-x_in_view[0]);
395  DfDp_out_view(1,2) = +(2.0*f*f/(L*L*L))*(a-x_in_view[0]);
396 
397  // Compute df/dp + (df/dx_dot) * (dx_dot/dp) + (df/dx) * (dx/dp)
398  if (useDfDpAsTangent_ && !is_null(DxdotDp_in)) {
399  Thyra::ConstDetachedMultiVectorView<Scalar> DxdotDp( *DxdotDp_in );
400  DfDp_out_view(0,0) += DxdotDp(0,0);
401  DfDp_out_view(0,1) += DxdotDp(0,1);
402  DfDp_out_view(0,2) += DxdotDp(0,2);
403  DfDp_out_view(1,0) += DxdotDp(1,0);
404  DfDp_out_view(1,1) += DxdotDp(1,1);
405  DfDp_out_view(1,2) += DxdotDp(1,2);
406  }
407  if (useDfDpAsTangent_ && !is_null(DxDp_in)) {
408  Thyra::ConstDetachedMultiVectorView<Scalar> DxDp( *DxDp_in );
409  DfDp_out_view(0,0) += -DxDp(1,0);
410  DfDp_out_view(0,1) += -DxDp(1,1);
411  DfDp_out_view(0,2) += -DxDp(1,2);
412  DfDp_out_view(1,0) += (f/L)*(f/L) * DxDp(0,0);
413  DfDp_out_view(1,1) += (f/L)*(f/L) * DxDp(0,1);
414  DfDp_out_view(1,2) += (f/L)*(f/L) * DxDp(0,2);
415  }
416  }
417  }
418 
419  // Responses: g = x
420  if (acceptModelParams_) {
421  RCP<VectorBase<Scalar> > g_out = outArgs.get_g(0);
422  if (g_out != Teuchos::null)
423  Thyra::assign(g_out.ptr(), *x_in);
424 
425  RCP<Thyra::MultiVectorBase<Scalar> > DgDp_out =
426  outArgs.get_DgDp(0,0).getMultiVector();
427  if (DgDp_out != Teuchos::null)
428  Thyra::assign(DgDp_out.ptr(), Scalar(0.0));
429 
430  RCP<Thyra::MultiVectorBase<Scalar> > DgDx_out =
431  outArgs.get_DgDx(0).getMultiVector();
432  if (DgDx_out != Teuchos::null) {
433  Thyra::DetachedMultiVectorView<Scalar> DgDx_out_view( *DgDx_out );
434  DgDx_out_view(0,0) = 1.0;
435  DgDx_out_view(0,1) = 0.0;
436  DgDx_out_view(1,0) = 0.0;
437  DgDx_out_view(1,1) = 1.0;
438  }
439  }
440 }
441 
442 template<class Scalar>
443 Teuchos::RCP<const Thyra::VectorSpaceBase<Scalar> >
444 SinCosModel<Scalar>::
445 get_p_space(int l) const
446 {
447  if (!acceptModelParams_) {
448  return Teuchos::null;
449  }
450  TEUCHOS_ASSERT_IN_RANGE_UPPER_EXCLUSIVE( l, 0, Np_ );
451  if (l == 0)
452  return p_space_;
453  else if (l == 1 || l == 2)
454  return DxDp_space_;
455  return Teuchos::null;
456 }
457 
458 template<class Scalar>
459 Teuchos::RCP<const Teuchos::Array<std::string> >
460 SinCosModel<Scalar>::
461 get_p_names(int l) const
462 {
463  if (!acceptModelParams_) {
464  return Teuchos::null;
465  }
466  TEUCHOS_ASSERT_IN_RANGE_UPPER_EXCLUSIVE( l, 0, Np_ );
467  Teuchos::RCP<Teuchos::Array<std::string> > p_strings =
468  Teuchos::rcp(new Teuchos::Array<std::string>());
469  if (l == 0) {
470  p_strings->push_back("Model Coefficient: a");
471  p_strings->push_back("Model Coefficient: f");
472  p_strings->push_back("Model Coefficient: L");
473  }
474  else if (l == 1)
475  p_strings->push_back("DxDp");
476  else if (l == 2)
477  p_strings->push_back("Dx_dotDp");
478  return p_strings;
479 }
480 
481 template<class Scalar>
482 Teuchos::RCP<const Thyra::VectorSpaceBase<Scalar> >
483 SinCosModel<Scalar>::
484 get_g_space(int j) const
485 {
486  TEUCHOS_ASSERT_IN_RANGE_UPPER_EXCLUSIVE( j, 0, Ng_ );
487  return g_space_;
488 }
489 
490 // private
491 
492 template<class Scalar>
493 void
494 SinCosModel<Scalar>::
495 setupInOutArgs_() const
496 {
497  if (isInitialized_) {
498  return;
499  }
500 
501  {
502  // Set up prototypical InArgs
503  Thyra::ModelEvaluatorBase::InArgsSetup<Scalar> inArgs;
504  inArgs.setModelEvalDescription(this->description());
505  inArgs.setSupports( Thyra::ModelEvaluatorBase::IN_ARG_t );
506  inArgs.setSupports( Thyra::ModelEvaluatorBase::IN_ARG_x );
507  inArgs.setSupports( Thyra::ModelEvaluatorBase::IN_ARG_beta );
508  inArgs.setSupports( Thyra::ModelEvaluatorBase::IN_ARG_x_dot );
509  inArgs.setSupports( Thyra::ModelEvaluatorBase::IN_ARG_alpha );
510  if (acceptModelParams_) {
511  inArgs.set_Np(Np_);
512  }
513  inArgs_ = inArgs;
514  }
515 
516  {
517  // Set up prototypical OutArgs
518  Thyra::ModelEvaluatorBase::OutArgsSetup<Scalar> outArgs;
519  outArgs.setModelEvalDescription(this->description());
520  outArgs.setSupports( Thyra::ModelEvaluatorBase::OUT_ARG_f );
521  outArgs.setSupports( Thyra::ModelEvaluatorBase::OUT_ARG_W_op );
522  if (acceptModelParams_) {
523  outArgs.set_Np_Ng(Np_,Ng_);
524  outArgs.setSupports( Thyra::ModelEvaluatorBase::OUT_ARG_DfDp,0,
525  Thyra::ModelEvaluatorBase::DERIV_MV_BY_COL );
526  outArgs.setSupports( Thyra::ModelEvaluatorBase::OUT_ARG_DgDp,0,0,
527  Thyra::ModelEvaluatorBase::DERIV_MV_BY_COL );
528  outArgs.setSupports( Thyra::ModelEvaluatorBase::OUT_ARG_DgDx,0,
529  Thyra::ModelEvaluatorBase::DERIV_TRANS_MV_BY_ROW );
530  }
531  outArgs_ = outArgs;
532  }
533 
534  // Set up nominal values
535  nominalValues_ = inArgs_;
536  if (haveIC_)
537  {
538  nominalValues_.set_t(t0_ic_);
539  const Teuchos::RCP<Thyra::VectorBase<Scalar> > x_ic = createMember(x_space_);
540  { // scope to delete DetachedVectorView
541  Thyra::DetachedVectorView<Scalar> x_ic_view( *x_ic );
542  x_ic_view[0] = a_+b_*sin((f_/L_)*t0_ic_+phi_);
543  x_ic_view[1] = b_*(f_/L_)*cos((f_/L_)*t0_ic_+phi_);
544  }
545  nominalValues_.set_x(x_ic);
546  if (acceptModelParams_) {
547  const Teuchos::RCP<Thyra::VectorBase<Scalar> > p_ic = createMember(p_space_);
548  {
549  Thyra::DetachedVectorView<Scalar> p_ic_view( *p_ic );
550  p_ic_view[0] = a_;
551  p_ic_view[1] = f_;
552  p_ic_view[2] = L_;
553  }
554  nominalValues_.set_p(0,p_ic);
555  }
556  const Teuchos::RCP<Thyra::VectorBase<Scalar> > x_dot_ic = createMember(x_space_);
557  { // scope to delete DetachedVectorView
558  Thyra::DetachedVectorView<Scalar> x_dot_ic_view( *x_dot_ic );
559  x_dot_ic_view[0] = b_*(f_/L_)*cos((f_/L_)*t0_ic_+phi_);
560  x_dot_ic_view[1] = -b_*(f_/L_)*(f_/L_)*sin((f_/L_)*t0_ic_+phi_);
561  }
562  nominalValues_.set_x_dot(x_dot_ic);
563  }
564 
565  isInitialized_ = true;
566 
567 }
568 
569 template<class Scalar>
570 void
571 SinCosModel<Scalar>::
572 setParameterList(Teuchos::RCP<Teuchos::ParameterList> const& paramList)
573 {
574  using Teuchos::get;
575  using Teuchos::ParameterList;
576  Teuchos::RCP<ParameterList> tmpPL = Teuchos::rcp(new ParameterList("SinCosModel"));
577  if (paramList != Teuchos::null) tmpPL = paramList;
578  tmpPL->validateParametersAndSetDefaults(*this->getValidParameters());
579  this->setMyParamList(tmpPL);
580  Teuchos::RCP<ParameterList> pl = this->getMyNonconstParamList();
581  bool acceptModelParams = get<bool>(*pl,"Accept model parameters");
582  bool haveIC = get<bool>(*pl,"Provide nominal values");
583  bool useDfDpAsTangent = get<bool>(*pl, "Use DfDp as Tangent");
584  if ( (acceptModelParams != acceptModelParams_) ||
585  (haveIC != haveIC_)
586  ) {
587  isInitialized_ = false;
588  }
589  acceptModelParams_ = acceptModelParams;
590  haveIC_ = haveIC;
591  useDfDpAsTangent_ = useDfDpAsTangent;
592  a_ = get<Scalar>(*pl,"Coeff a");
593  f_ = get<Scalar>(*pl,"Coeff f");
594  L_ = get<Scalar>(*pl,"Coeff L");
595  x0_ic_ = get<Scalar>(*pl,"IC x0");
596  x1_ic_ = get<Scalar>(*pl,"IC x1");
597  t0_ic_ = get<Scalar>(*pl,"IC t0");
598  calculateCoeffFromIC_();
599  setupInOutArgs_();
600 }
601 
602 template<class Scalar>
603 Teuchos::RCP<const Teuchos::ParameterList>
604 SinCosModel<Scalar>::
605 getValidParameters() const
606 {
607  static Teuchos::RCP<const Teuchos::ParameterList> validPL;
608  if (is_null(validPL)) {
609  Teuchos::RCP<Teuchos::ParameterList> pl = Teuchos::parameterList();
610  pl->set("Accept model parameters", false);
611  pl->set("Provide nominal values", true);
612  pl->set("Use DfDp as Tangent", false);
613  pl->set<std::string>("Output File Name", "Tempus_BDF2_SinCos");
614  Teuchos::setDoubleParameter(
615  "Coeff a", 0.0, "Coefficient a in model", &*pl);
616  Teuchos::setDoubleParameter(
617  "Coeff f", 1.0, "Coefficient f in model", &*pl);
618  Teuchos::setDoubleParameter(
619  "Coeff L", 1.0, "Coefficient L in model", &*pl);
620  Teuchos::setDoubleParameter(
621  "IC x0", 0.0, "Initial Condition for x0", &*pl);
622  Teuchos::setDoubleParameter(
623  "IC x1", 1.0, "Initial Condition for x1", &*pl);
624  Teuchos::setDoubleParameter(
625  "IC t0", 0.0, "Initial time t0", &*pl);
626  Teuchos::setIntParameter(
627  "Number of Time Step Sizes", 1, "Number time step sizes for convergence study", &*pl);
628  validPL = pl;
629  }
630  return validPL;
631 }
632 
633 template<class Scalar>
634 void
635 SinCosModel<Scalar>::
636 calculateCoeffFromIC_()
637 {
638  phi_ = atan(((f_/L_)/x1_ic_)*(x0_ic_-a_))-(f_/L_)*t0_ic_;
639  b_ = x1_ic_/((f_/L_)*cos((f_/L_)*t0_ic_+phi_));
640 }
641 
642 } // namespace Tempus_Test
643 #endif // TEMPUS_TEST_SINCOS_MODEL_IMPL_HPP
Tempus_Test
Definition: Tempus_BackwardEuler_ASA.cpp:34
Tempus_Test::SinCosModel::get_p_space
Teuchos::RCP< const Thyra::VectorSpaceBase< Scalar > > get_p_space(int l) const
Definition: SinCosModel_impl.hpp:445
Tempus_Test::SinCosModel::create_W_op
Teuchos::RCP< Thyra::LinearOpBase< Scalar > > create_W_op() const
Definition: SinCosModel_impl.hpp:224
Tempus_Test::SinCosModel::get_x_space
Teuchos::RCP< const Thyra::VectorSpaceBase< Scalar > > get_x_space() const
Definition: SinCosModel_impl.hpp:149
Tempus_Test::SinCosModel::getNominalValues
Thyra::ModelEvaluatorBase::InArgs< Scalar > getNominalValues() const
Definition: SinCosModel_impl.hpp:167
Tempus_Test::SinCosModel::getExactSensSolution
Thyra::ModelEvaluatorBase::InArgs< Scalar > getExactSensSolution(int j, double t) const
Definition: SinCosModel_impl.hpp:98
Tempus_Test::SinCosModel::setParameterList
void setParameterList(Teuchos::RCP< Teuchos::ParameterList > const &paramList)
Definition: SinCosModel_impl.hpp:572
Tempus_Test::SinCosModel::get_W_factory
Teuchos::RCP< const Thyra::LinearOpWithSolveFactoryBase< Scalar > > get_W_factory() const
Definition: SinCosModel_impl.hpp:234
Tempus_Test::SinCosModel::createOutArgsImpl
Thyra::ModelEvaluatorBase::OutArgs< Scalar > createOutArgsImpl() const
Definition: SinCosModel_impl.hpp:269
Tempus_Test::SinCosModel::create_W
Teuchos::RCP< Thyra::LinearOpWithSolveBase< Scalar > > create_W() const
Definition: SinCosModel_impl.hpp:178
Tempus_Test::SinCosModel::get_f_space
Teuchos::RCP< const Thyra::VectorSpaceBase< Scalar > > get_f_space() const
Definition: SinCosModel_impl.hpp:158
Tempus_Test::SinCosModel::evalModelImpl
void evalModelImpl(const Thyra::ModelEvaluatorBase::InArgs< Scalar > &inArgs_bar, const Thyra::ModelEvaluatorBase::OutArgs< Scalar > &outArgs_bar) const
Definition: SinCosModel_impl.hpp:279
Tempus_Test::SinCosModel::SinCosModel
SinCosModel(Teuchos::RCP< Teuchos::ParameterList > pList=Teuchos::null)
Definition: SinCosModel_impl.hpp:31
Tempus_Test::SinCosModel::getValidParameters
Teuchos::RCP< const Teuchos::ParameterList > getValidParameters() const
Definition: SinCosModel_impl.hpp:605
Tempus_Test::SinCosModel::createInArgs
Thyra::ModelEvaluatorBase::InArgs< Scalar > createInArgs() const
Definition: SinCosModel_impl.hpp:256
Tempus_Test::SinCosModel::get_g_space
Teuchos::RCP< const Thyra::VectorSpaceBase< Scalar > > get_g_space(int j) const
Definition: SinCosModel_impl.hpp:484
Tempus_Test::SinCosModel::getExactSolution
Thyra::ModelEvaluatorBase::InArgs< Scalar > getExactSolution(double t) const
Definition: SinCosModel_impl.hpp:65
Tempus_Test::SinCosModel::setupInOutArgs_
void setupInOutArgs_() const
Definition: SinCosModel_impl.hpp:495
Tempus_Test::SinCosModel::get_p_names
Teuchos::RCP< const Teuchos::Array< std::string > > get_p_names(int l) const
Definition: SinCosModel_impl.hpp:461
Tempus_Test::SinCosModel::calculateCoeffFromIC_
void calculateCoeffFromIC_()
Definition: SinCosModel_impl.hpp:636