Tempus  Version of the Day
Time Integration
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
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>
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>
65 getExactSolution(double t) const
66 {
67  TEUCHOS_TEST_FOR_EXCEPTION( !isInitialized_, std::logic_error,
68  "Error, setupInOutArgs_ must be called first!\n");
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>
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");
103  if (!acceptModelParams_) {
104  return inArgs;
105  }
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>
149 get_x_space() const
150 {
151  return x_space_;
152 }
153 
154 
155 template<class Scalar>
158 get_f_space() const
159 {
160  return f_space_;
161 }
162 
163 
164 template<class Scalar>
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>
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  }
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>
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>
235 {
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>
257 {
258  setupInOutArgs_();
259  return inArgs_;
260 }
261 
262 
263 // Private functions overridden from ModelEvaluatorDefaultBase
264 
265 
266 template<class Scalar>
270 {
271  setupInOutArgs_();
272  return outArgs_;
273 }
274 
275 
276 template<class Scalar>
277 void
282  ) const
283 {
285  using Teuchos::RCP;
286  using Thyra::VectorBase;
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();
323  if (acceptModelParams_) {
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)) {
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)) {
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)) {
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)) {
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 
426  outArgs.get_DgDp(0,0).getMultiVector();
427  if (DgDp_out != Teuchos::null)
428  Thyra::assign(DgDp_out.ptr(), Scalar(0.0));
429 
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>
445 get_p_space(int l) const
446 {
447  if (!acceptModelParams_) {
448  return Teuchos::null;
449  }
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>
461 get_p_names(int l) const
462 {
463  if (!acceptModelParams_) {
464  return Teuchos::null;
465  }
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>
484 get_g_space(int j) const
485 {
487  return g_space_;
488 }
489 
490 // private
491 
492 template<class Scalar>
493 void
496 {
497  if (isInitialized_) {
498  return;
499  }
500 
501  using Teuchos::RCP;
502  typedef Thyra::ModelEvaluatorBase MEB;
503  {
504  // Set up prototypical InArgs
505  MEB::InArgsSetup<Scalar> inArgs;
506  inArgs.setModelEvalDescription(this->description());
507  inArgs.setSupports( MEB::IN_ARG_t );
508  inArgs.setSupports( MEB::IN_ARG_x );
509  inArgs.setSupports( MEB::IN_ARG_beta );
510  inArgs.setSupports( MEB::IN_ARG_x_dot );
511  inArgs.setSupports( MEB::IN_ARG_alpha );
512  if (acceptModelParams_) {
513  inArgs.set_Np(Np_);
514  }
515  inArgs_ = inArgs;
516  }
517 
518  {
519  // Set up prototypical OutArgs
520  MEB::OutArgsSetup<Scalar> outArgs;
521  outArgs.setModelEvalDescription(this->description());
522  outArgs.setSupports( MEB::OUT_ARG_f );
523  outArgs.setSupports( MEB::OUT_ARG_W_op );
524  if (acceptModelParams_) {
525  outArgs.set_Np_Ng(Np_,Ng_);
526  outArgs.setSupports( MEB::OUT_ARG_DfDp,0,
527  MEB::DERIV_MV_JACOBIAN_FORM );
528  outArgs.setSupports( MEB::OUT_ARG_DgDp,0,0,
529  MEB::DERIV_MV_JACOBIAN_FORM );
530  outArgs.setSupports( MEB::OUT_ARG_DgDx,0,
531  MEB::DERIV_MV_GRADIENT_FORM );
532  }
533  outArgs_ = outArgs;
534  }
535 
536  // Set up nominal values
537  nominalValues_ = inArgs_;
538  if (haveIC_)
539  {
540  nominalValues_.set_t(t0_ic_);
541  const RCP<Thyra::VectorBase<Scalar> > x_ic = createMember(x_space_);
542  { // scope to delete DetachedVectorView
543  Thyra::DetachedVectorView<Scalar> x_ic_view( *x_ic );
544  x_ic_view[0] = a_+b_*sin((f_/L_)*t0_ic_+phi_);
545  x_ic_view[1] = b_*(f_/L_)*cos((f_/L_)*t0_ic_+phi_);
546  }
547  nominalValues_.set_x(x_ic);
548  if (acceptModelParams_) {
549  const RCP<Thyra::VectorBase<Scalar> > p_ic = createMember(p_space_);
550  {
551  Thyra::DetachedVectorView<Scalar> p_ic_view( *p_ic );
552  p_ic_view[0] = a_;
553  p_ic_view[1] = f_;
554  p_ic_view[2] = L_;
555  }
556  nominalValues_.set_p(0,p_ic);
557  }
558  const RCP<Thyra::VectorBase<Scalar> > x_dot_ic = createMember(x_space_);
559  { // scope to delete DetachedVectorView
560  Thyra::DetachedVectorView<Scalar> x_dot_ic_view( *x_dot_ic );
561  x_dot_ic_view[0] = b_*(f_/L_)*cos((f_/L_)*t0_ic_+phi_);
562  x_dot_ic_view[1] = -b_*(f_/L_)*(f_/L_)*sin((f_/L_)*t0_ic_+phi_);
563  }
564  nominalValues_.set_x_dot(x_dot_ic);
565  }
566 
567  isInitialized_ = true;
568 
569 }
570 
571 template<class Scalar>
572 void
575 {
576  using Teuchos::get;
578  Teuchos::RCP<ParameterList> tmpPL = Teuchos::rcp(new ParameterList("SinCosModel"));
579  if (paramList != Teuchos::null) tmpPL = paramList;
580  tmpPL->validateParametersAndSetDefaults(*this->getValidParameters());
581  this->setMyParamList(tmpPL);
582  Teuchos::RCP<ParameterList> pl = this->getMyNonconstParamList();
583  bool acceptModelParams = get<bool>(*pl,"Accept model parameters");
584  bool haveIC = get<bool>(*pl,"Provide nominal values");
585  bool useDfDpAsTangent = get<bool>(*pl, "Use DfDp as Tangent");
586  if ( (acceptModelParams != acceptModelParams_) ||
587  (haveIC != haveIC_)
588  ) {
589  isInitialized_ = false;
590  }
591  acceptModelParams_ = acceptModelParams;
592  haveIC_ = haveIC;
593  useDfDpAsTangent_ = useDfDpAsTangent;
594  a_ = get<Scalar>(*pl,"Coeff a");
595  f_ = get<Scalar>(*pl,"Coeff f");
596  L_ = get<Scalar>(*pl,"Coeff L");
597  x0_ic_ = get<Scalar>(*pl,"IC x0");
598  x1_ic_ = get<Scalar>(*pl,"IC x1");
599  t0_ic_ = get<Scalar>(*pl,"IC t0");
600  calculateCoeffFromIC_();
601  setupInOutArgs_();
602 }
603 
604 template<class Scalar>
608 {
610  if (is_null(validPL)) {
611  Teuchos::RCP<Teuchos::ParameterList> pl = Teuchos::parameterList();
612  pl->set("Accept model parameters", false);
613  pl->set("Provide nominal values", true);
614  pl->set("Use DfDp as Tangent", false);
615  pl->set<std::string>("Output File Name", "Tempus_BDF2_SinCos");
616  Teuchos::setDoubleParameter(
617  "Coeff a", 0.0, "Coefficient a in model", &*pl);
618  Teuchos::setDoubleParameter(
619  "Coeff f", 1.0, "Coefficient f in model", &*pl);
620  Teuchos::setDoubleParameter(
621  "Coeff L", 1.0, "Coefficient L in model", &*pl);
622  Teuchos::setDoubleParameter(
623  "IC x0", 0.0, "Initial Condition for x0", &*pl);
624  Teuchos::setDoubleParameter(
625  "IC x1", 1.0, "Initial Condition for x1", &*pl);
626  Teuchos::setDoubleParameter(
627  "IC t0", 0.0, "Initial time t0", &*pl);
628  Teuchos::setIntParameter(
629  "Number of Time Step Sizes", 1, "Number time step sizes for convergence study", &*pl);
630  validPL = pl;
631  }
632  return validPL;
633 }
634 
635 template<class Scalar>
636 void
639 {
640  phi_ = atan(((f_/L_)/x1_ic_)*(x0_ic_-a_))-(f_/L_)*t0_ic_;
641  b_ = x1_ic_/((f_/L_)*cos((f_/L_)*t0_ic_+phi_));
642 }
643 
644 template<class Scalar>
647 create_W() const
648 {
649  using Teuchos::RCP;
650  RCP<const Thyra::LinearOpWithSolveFactoryBase<Scalar> > W_factory = this->get_W_factory();
651  RCP<Thyra::LinearOpBase<Scalar> > matrix = this->create_W_op();
652  {
653  // 01/20/09 tscoffe: This is a total hack to provide a full rank matrix to
654  // linearOpWithSolve because it ends up factoring the matrix during
655  // initialization, which it really shouldn't do, or I'm doing something
656  // wrong here. The net effect is that I get exceptions thrown in
657  // optimized mode due to the matrix being rank deficient unless I do this.
658  RCP<Thyra::MultiVectorBase<Scalar> > multivec = Teuchos::rcp_dynamic_cast<Thyra::MultiVectorBase<Scalar> >(matrix,true);
659  {
660  RCP<Thyra::VectorBase<Scalar> > vec = Thyra::createMember(this->f_space_);
661  {
662  Thyra::DetachedVectorView<Scalar> vec_view( *vec );
663  vec_view[0] = 0.0;
664  vec_view[1] = 1.0;
665  }
666  V_V(multivec->col(0).ptr(),*vec);
667  {
668  Thyra::DetachedVectorView<Scalar> vec_view( *vec );
669  vec_view[0] = 1.0;
670  vec_view[1] = 0.0;
671  }
672  V_V(multivec->col(1).ptr(),*vec);
673  }
674  }
676  Thyra::linearOpWithSolve<Scalar>(
677  *W_factory,
678  matrix
679  );
680  return W;
681 }
682 
683 template<class Scalar>
686 create_W_op() const
687 {
688  Teuchos::RCP<Thyra::MultiVectorBase<Scalar> > matrix = Thyra::createMembers(this->f_space_, this->dim_);
689  return(matrix);
690 }
691 
692 template<class Scalar>
696 {
697  // This ME should use the same InArgs as the base SinCosModel. However
698  // we can't just use it's InArgs directly because the description won't
699  // match (which is checked in debug builds). Instead create a new InArgsSetup
700  // initialized by SinCosModel::createInArgs() and set the description
701  // appropriately.
702  typedef Thyra::ModelEvaluatorBase MEB;
703  MEB::InArgsSetup<Scalar> inArgs = SinCosModel<Scalar>::createInArgs();
704  inArgs.setModelEvalDescription(this->description());
705  return inArgs;
706 }
707 
708 template<class Scalar>
712 {
713  typedef Thyra::ModelEvaluatorBase MEB;
714  MEB::OutArgsSetup<Scalar> outArgs;
715  outArgs.setModelEvalDescription(this->description());
716  outArgs.setSupports( MEB::OUT_ARG_f ); // Apparently all models have to support f
717  outArgs.setSupports( MEB::OUT_ARG_W_op );
718  outArgs.set_Np_Ng(this->Np_,0);
719  return outArgs;
720 }
721 
722 template<class Scalar>
723 void
728  ) const
729 {
730  using Teuchos::RCP;
731  using Thyra::VectorBase;
733  using Teuchos::rcp_dynamic_cast;
734  TEUCHOS_TEST_FOR_EXCEPTION( !this->isInitialized_, std::logic_error,
735  "Error, setupInOutArgs_ must be called first!\n");
736 
737  const RCP<const VectorBase<Scalar> > x_in = inArgs.get_x().assert_not_null();
738  Thyra::ConstDetachedVectorView<Scalar> x_in_view( *x_in );
739 
740  //double t = inArgs.get_t();
741  //Scalar a = this->a_;
742  Scalar f = this->f_;
743  Scalar L = this->L_;
744  if (this->acceptModelParams_) {
745  const RCP<const VectorBase<Scalar> > p_in =
746  inArgs.get_p(0).assert_not_null();
747  Thyra::ConstDetachedVectorView<Scalar> p_in_view( *p_in );
748  //a = p_in_view[0];
749  f = p_in_view[1];
750  L = p_in_view[2];
751  }
752 
753  Scalar beta = inArgs.get_beta();
754 
755  const RCP<Thyra::LinearOpBase<Scalar> > W_out = outArgs.get_W_op();
756  if (inArgs.get_x_dot().is_null()) {
757 
758  // Evaluate the Explicit ODE f(x,t) [= 0]
759  if (!is_null(W_out)) {
761  Teuchos::rcp_dynamic_cast<Thyra::MultiVectorBase<Scalar> >(W_out,true);
762  Thyra::DetachedMultiVectorView<Scalar> matrix_view( *matrix );
763  matrix_view(0,0) = 0.0; // d(f0)/d(x0_n)
764  matrix_view(1,0) = +beta; // d(f0)/d(x1_n)
765  matrix_view(0,1) = -beta*(f/L)*(f/L); // d(f1)/d(x0_n)
766  matrix_view(1,1) = 0.0; // d(f1)/d(x1_n)
767  // Note: alpha = d(xdot)/d(x_n) and beta = d(x)/d(x_n)
768  }
769  } else {
770 
771  // Evaluate the implicit ODE f(xdot, x, t) [= 0]
772  RCP<const VectorBase<Scalar> > x_dot_in;
773  x_dot_in = inArgs.get_x_dot().assert_not_null();
774  Scalar alpha = inArgs.get_alpha();
775  if (!is_null(W_out)) {
777  Teuchos::rcp_dynamic_cast<Thyra::MultiVectorBase<Scalar> >(W_out,true);
778  Thyra::DetachedMultiVectorView<Scalar> matrix_view( *matrix );
779  matrix_view(0,0) = alpha; // d(f0)/d(x0_n)
780  matrix_view(1,0) = -beta; // d(f0)/d(x1_n)
781  matrix_view(0,1) = +beta*(f/L)*(f/L); // d(f1)/d(x0_n)
782  matrix_view(1,1) = alpha; // d(f1)/d(x1_n)
783  // Note: alpha = d(xdot)/d(x_n) and beta = d(x)/d(x_n)
784  }
785  }
786 }
787 
788 } // namespace Tempus_Test
789 #endif // TEMPUS_TEST_SINCOS_MODEL_IMPL_HPP
bool is_null(const boost::shared_ptr< T > &p)
Derivative< Scalar > get_DfDp(int l) const
RCP< const VectorBase< Scalar > > get_x_dot() 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)
Thyra::ModelEvaluatorBase::InArgs< Scalar > createInArgs() const
Evaluation< VectorBase< Scalar > > get_g(int j) const
Evaluation< VectorBase< Scalar > > get_f() const
Teuchos::RCP< Thyra::LinearOpBase< Scalar > > create_W_op() const
Teuchos::RCP< Thyra::LinearOpWithSolveBase< Scalar > > create_W() const
Thyra::ModelEvaluatorBase::InArgs< Scalar > getNominalValues() const
Derivative< Scalar > get_DgDp(int j, int l) const
Teuchos::RCP< const Thyra::VectorSpaceBase< Scalar > > get_f_space() const
Teuchos::RCP< Thyra::LinearOpBase< Scalar > > create_W_op() const
#define TEUCHOS_ASSERT_IN_RANGE_UPPER_EXCLUSIVE(index, lower_inclusive, upper_exclusive)
void set_x(const RCP< const VectorBase< Scalar > > &x)
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
void evalModelImpl(const Thyra::ModelEvaluatorBase::InArgs< Scalar > &inArgs_bar, const Thyra::ModelEvaluatorBase::OutArgs< Scalar > &outArgs_bar) const
Teuchos::RCP< const Thyra::VectorSpaceBase< Scalar > > get_x_space() const
Thyra::ModelEvaluatorBase::InArgs< Scalar > createInArgs() const
Thyra::ModelEvaluatorBase::InArgs< Scalar > getExactSolution(double t) const
Ptr< T > ptr() const
SinCosModel(Teuchos::RCP< Teuchos::ParameterList > pList=Teuchos::null)
void validateParametersAndSetDefaults(ParameterList const &validParamList, int const depth=1000)
Teuchos::RCP< const Thyra::VectorSpaceBase< Scalar > > get_p_space(int l) const
Thyra::ModelEvaluatorBase::OutArgs< Scalar > createOutArgsImpl() const
Teuchos::RCP< const Teuchos::ParameterList > getValidParameters() const
void set_x_dot(const RCP< const VectorBase< Scalar > > &x_dot)
void evalModelImpl(const Thyra::ModelEvaluatorBase::InArgs< Scalar > &inArgs_bar, const Thyra::ModelEvaluatorBase::OutArgs< Scalar > &outArgs_bar) const
Teuchos::RCP< const Thyra::VectorSpaceBase< Scalar > > get_g_space(int j) const
Teuchos::RCP< const Teuchos::Array< std::string > > get_p_names(int l) const
Thyra::ModelEvaluatorBase::OutArgs< Scalar > createOutArgsImpl() const
Thyra::ModelEvaluatorBase::InArgs< Scalar > getExactSensSolution(int j, double t) const
Teuchos::RCP< const Thyra::LinearOpWithSolveFactoryBase< Scalar > > get_W_factory() const
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
Teuchos::RCP< Thyra::LinearOpWithSolveBase< Scalar > > create_W() const
void setParameterList(Teuchos::RCP< Teuchos::ParameterList > const &paramList)
Derivative< Scalar > get_DgDx(int j) const