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