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: 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
27namespace Tempus_Test {
28
29template <class Scalar>
30SinCosModel<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
61template <class Scalar>
62Thyra::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//
96template <class Scalar>
97Thyra::ModelEvaluatorBase::InArgs<Scalar>
98SinCosModel<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
159template <class Scalar>
160Teuchos::RCP<const Thyra::VectorSpaceBase<Scalar> >
161SinCosModel<Scalar>::get_x_space() const
162{
163 return x_space_;
164}
165
166template <class Scalar>
167Teuchos::RCP<const Thyra::VectorSpaceBase<Scalar> >
168SinCosModel<Scalar>::get_f_space() const
169{
170 return f_space_;
171}
172
173template <class Scalar>
174Thyra::ModelEvaluatorBase::InArgs<Scalar>
175SinCosModel<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
182template <class Scalar>
183Teuchos::RCP<Thyra::LinearOpWithSolveBase<Scalar> >
184SinCosModel<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
226template <class Scalar>
227Teuchos::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
235template <class Scalar>
236Teuchos::RCP<const Thyra::LinearOpWithSolveFactoryBase<Scalar> >
237SinCosModel<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
255template <class Scalar>
256Thyra::ModelEvaluatorBase::InArgs<Scalar> SinCosModel<Scalar>::createInArgs()
257 const
258{
259 setupInOutArgs_();
260 return inArgs_;
261}
262
263// Private functions overridden from ModelEvaluatorDefaultBase
264
265template <class Scalar>
266Thyra::ModelEvaluatorBase::OutArgs<Scalar>
267SinCosModel<Scalar>::createOutArgsImpl() const
268{
269 setupInOutArgs_();
270 return outArgs_;
271}
272
273template <class Scalar>
274void 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
435template <class Scalar>
436Teuchos::RCP<const Thyra::VectorSpaceBase<Scalar> >
437SinCosModel<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
450template <class Scalar>
451Teuchos::RCP<const Teuchos::Array<std::string> >
452SinCosModel<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
472template <class Scalar>
473Teuchos::RCP<const Thyra::VectorSpaceBase<Scalar> >
474SinCosModel<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
482template <class Scalar>
483void 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
555template <class Scalar>
556void 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
586template <class Scalar>
587Teuchos::RCP<const Teuchos::ParameterList>
588SinCosModel<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
611template <class Scalar>
612void 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
618template <class Scalar>
619Teuchos::RCP<Thyra::LinearOpWithSolveBase<Scalar> >
620SinCosModelAdjoint<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
656template <class Scalar>
657Teuchos::RCP<Thyra::LinearOpBase<Scalar> >
658SinCosModelAdjoint<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
665template <class Scalar>
666Thyra::ModelEvaluatorBase::InArgs<Scalar>
667SinCosModelAdjoint<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
680template <class Scalar>
681Thyra::ModelEvaluatorBase::OutArgs<Scalar>
682SinCosModelAdjoint<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
694template <class Scalar>
695void 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
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::SinCosModelAdjoint::createOutArgsImpl
Thyra::ModelEvaluatorBase::OutArgs< Scalar > createOutArgsImpl() const
Definition SinCosModel_impl.hpp:682
Tempus_Test::SinCosModelAdjoint::create_W
Teuchos::RCP< Thyra::LinearOpWithSolveBase< Scalar > > create_W() const
Definition SinCosModel_impl.hpp:620
Tempus_Test::SinCosModelAdjoint::createInArgs
Thyra::ModelEvaluatorBase::InArgs< Scalar > createInArgs() const
Definition SinCosModel_impl.hpp:667
Tempus_Test::SinCosModelAdjoint::create_W_op
Teuchos::RCP< Thyra::LinearOpBase< Scalar > > create_W_op() const
Definition SinCosModel_impl.hpp:658
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::get_f_space
Teuchos::RCP< const Thyra::VectorSpaceBase< Scalar > > get_f_space() const
Definition SinCosModel_impl.hpp:168
Tempus_Test::SinCosModel::createInArgs
Thyra::ModelEvaluatorBase::InArgs< Scalar > createInArgs() const
Definition SinCosModel_impl.hpp:256
Tempus_Test::SinCosModel::create_W_op
Teuchos::RCP< Thyra::LinearOpBase< Scalar > > create_W_op() const
Definition SinCosModel_impl.hpp:227
Tempus_Test::SinCosModel::get_W_factory
Teuchos::RCP< const Thyra::LinearOpWithSolveFactoryBase< Scalar > > get_W_factory() const
Definition SinCosModel_impl.hpp:237
Tempus_Test::SinCosModel::getExactSolution
Thyra::ModelEvaluatorBase::InArgs< Scalar > getExactSolution(double t) const
Definition SinCosModel_impl.hpp:62
Tempus_Test::SinCosModel::get_x_space
Teuchos::RCP< const Thyra::VectorSpaceBase< Scalar > > get_x_space() const
Definition SinCosModel_impl.hpp:161
Tempus_Test::SinCosModel::getValidParameters
Teuchos::RCP< const Teuchos::ParameterList > getValidParameters() const
Definition SinCosModel_impl.hpp:588
Tempus_Test::SinCosModel::create_W
Teuchos::RCP< Thyra::LinearOpWithSolveBase< Scalar > > create_W() const
Definition SinCosModel_impl.hpp:184
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::getNominalValues
Thyra::ModelEvaluatorBase::InArgs< Scalar > getNominalValues() const
Definition SinCosModel_impl.hpp:175
Tempus_Test::SinCosModel::getExactSensSolution
Thyra::ModelEvaluatorBase::InArgs< Scalar > getExactSensSolution(int j, double t) const
Definition SinCosModel_impl.hpp:98
Tempus_Test::SinCosModel::SinCosModel
SinCosModel(Teuchos::RCP< Teuchos::ParameterList > pList=Teuchos::null)
Definition SinCosModel_impl.hpp:30
Tempus_Test::SinCosModel::setParameterList
void setParameterList(Teuchos::RCP< Teuchos::ParameterList > const &paramList)
Definition SinCosModel_impl.hpp:556
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::setupInOutArgs_
void setupInOutArgs_() const
Definition SinCosModel_impl.hpp:483
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
Definition Tempus_BackwardEuler_ASA.cpp:32
Generated by doxygen 1.9.8