Tempus Version of the Day
Time Integration
Loading...
Searching...
No Matches
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
27namespace Tempus_Test {
28
29template<class Scalar>
31SinCosModel(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
62template<class Scalar>
63Thyra::ModelEvaluatorBase::InArgs<Scalar>
65getExactSolution(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//
95template<class Scalar>
96Thyra::ModelEvaluatorBase::InArgs<Scalar>
98getExactSensSolution(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
146template<class Scalar>
147Teuchos::RCP<const Thyra::VectorSpaceBase<Scalar> >
149get_x_space() const
150{
151 return x_space_;
152}
153
154
155template<class Scalar>
156Teuchos::RCP<const Thyra::VectorSpaceBase<Scalar> >
158get_f_space() const
159{
160 return f_space_;
161}
162
163
164template<class Scalar>
165Thyra::ModelEvaluatorBase::InArgs<Scalar>
167getNominalValues() 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
175template<class Scalar>
176Teuchos::RCP<Thyra::LinearOpWithSolveBase<Scalar> >
178create_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
221template<class Scalar>
222Teuchos::RCP<Thyra::LinearOpBase<Scalar> >
224create_W_op() const
225{
226 Teuchos::RCP<Thyra::MultiVectorBase<Scalar> > matrix = Thyra::createMembers(x_space_, dim_);
227 return(matrix);
228}
229
230
231template<class Scalar>
232Teuchos::RCP<const Thyra::LinearOpWithSolveFactoryBase<Scalar> >
234get_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
253template<class Scalar>
254Thyra::ModelEvaluatorBase::InArgs<Scalar>
256createInArgs() const
257{
258 setupInOutArgs_();
259 return inArgs_;
260}
261
262
263// Private functions overridden from ModelEvaluatorDefaultBase
264
265
266template<class Scalar>
267Thyra::ModelEvaluatorBase::OutArgs<Scalar>
269createOutArgsImpl() const
270{
271 setupInOutArgs_();
272 return outArgs_;
273}
274
275
276template<class Scalar>
277void
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;
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
442template<class Scalar>
443Teuchos::RCP<const Thyra::VectorSpaceBase<Scalar> >
445get_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
458template<class Scalar>
459Teuchos::RCP<const Teuchos::Array<std::string> >
461get_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
481template<class Scalar>
482Teuchos::RCP<const Thyra::VectorSpaceBase<Scalar> >
484get_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
492template<class Scalar>
493void
495setupInOutArgs_() const
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
571template<class Scalar>
572void
574setParameterList(Teuchos::RCP<Teuchos::ParameterList> const& paramList)
575{
576 using Teuchos::get;
577 using Teuchos::ParameterList;
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
604template<class Scalar>
605Teuchos::RCP<const Teuchos::ParameterList>
607getValidParameters() const
608{
609 static Teuchos::RCP<const Teuchos::ParameterList> validPL;
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
635template<class Scalar>
636void
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
644template<class Scalar>
645Teuchos::RCP<Thyra::LinearOpWithSolveBase<Scalar> >
647create_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 }
675 RCP<Thyra::LinearOpWithSolveBase<Scalar> > W =
676 Thyra::linearOpWithSolve<Scalar>(
677 *W_factory,
678 matrix
679 );
680 return W;
681}
682
683template<class Scalar>
684Teuchos::RCP<Thyra::LinearOpBase<Scalar> >
686create_W_op() const
687{
688 Teuchos::RCP<Thyra::MultiVectorBase<Scalar> > matrix = Thyra::createMembers(this->f_space_, this->dim_);
689 return(matrix);
690}
691
692template<class Scalar>
693Thyra::ModelEvaluatorBase::InArgs<Scalar>
695createInArgs() const
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
708template<class Scalar>
709Thyra::ModelEvaluatorBase::OutArgs<Scalar>
711createOutArgsImpl() const
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
722template<class Scalar>
723void
726 const Thyra::ModelEvaluatorBase::InArgs<Scalar> &inArgs,
727 const Thyra::ModelEvaluatorBase::OutArgs<Scalar> &outArgs
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)) {
760 RCP<Thyra::MultiVectorBase<Scalar> > matrix =
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)) {
776 RCP<Thyra::MultiVectorBase<Scalar> > matrix =
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
Sine-Cosine model problem from Rythmos. This is a canonical Sine-Cosine differential equation.
SinCosModel(Teuchos::RCP< Teuchos::ParameterList > pList=Teuchos::null)