Tempus Version of the Day
Time Integration
Loading...
Searching...
No Matches
DahlquistTestModel_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_DAHLQUIST_TEST_MODEL_IMPL_HPP
10#define TEMPUS_TEST_DAHLQUIST_TEST_MODEL_IMPL_HPP
11
12#include "Thyra_DefaultSpmdVectorSpace.hpp"
13#include "Thyra_DetachedVectorView.hpp"
14#include "Thyra_DetachedMultiVectorView.hpp"
15#include "Thyra_DefaultSerialDenseLinearOpWithSolveFactory.hpp"
16#include "Thyra_DefaultMultiVectorLinearOpWithSolve.hpp"
17#include "Thyra_DefaultLinearOpSource.hpp"
18#include "Thyra_VectorStdOps.hpp"
19#include "Thyra_MultiVectorStdOps.hpp"
20#include "Thyra_DefaultMultiVectorProductVector.hpp"
21
22//#include <iostream>
23
24
25namespace Tempus_Test {
26
27template<class Scalar>
29{ constructDahlquistTestModel(-1.0, false); }
30
31
32template<class Scalar>
34DahlquistTestModel(Scalar lambda, bool includeXDot)
35{ constructDahlquistTestModel(lambda, includeXDot); }
36
37
38template<class Scalar>
39void
41constructDahlquistTestModel(Scalar lambda, bool includeXDot)
42{
43 lambda_ = lambda;
44 includeXDot_ = includeXDot;
45 isInitialized_ = false;
46 xIC_ = Scalar( 1.0);
47 xDotIC_ = Scalar(lambda);
48 dim_ = 1;
49 haveIC_ = true;
50 Np_ = 1;
51 np_ = 1;
52 Ng_ = 1;
53 ng_ = dim_;
54 acceptModelParams_ = false;
55
56 // Create x_space and f_space
57 x_space_ = Thyra::defaultSpmdVectorSpace<Scalar>(dim_);
58 f_space_ = Thyra::defaultSpmdVectorSpace<Scalar>(dim_);
59
60 using Teuchos::RCP;
61 typedef Thyra::ModelEvaluatorBase MEB;
62 {
63 // Set up prototypical InArgs
64 MEB::InArgsSetup<Scalar> inArgs;
65 inArgs.setModelEvalDescription(this->description());
66 inArgs.setSupports( MEB::IN_ARG_t );
67 inArgs.setSupports( MEB::IN_ARG_x );
68 inArgs.setSupports( MEB::IN_ARG_x_dot );
69
70 inArgs.setSupports( MEB::IN_ARG_beta );
71 inArgs.setSupports( MEB::IN_ARG_alpha );
72 if (acceptModelParams_) inArgs.set_Np( Np_ );
73
74 inArgs_ = inArgs;
75 }
76
77 {
78 // Set up prototypical OutArgs
79 MEB::OutArgsSetup<Scalar> outArgs;
80 outArgs.setModelEvalDescription(this->description());
81 outArgs.setSupports( MEB::OUT_ARG_f );
82
83 outArgs.setSupports( MEB::OUT_ARG_W_op );
84 if (acceptModelParams_) {
85 outArgs.set_Np_Ng(Np_,Ng_);
86 outArgs.setSupports( MEB::OUT_ARG_DfDp,0,
87 MEB::DERIV_MV_JACOBIAN_FORM );
88 outArgs.setSupports( MEB::OUT_ARG_DgDp,0,0,
89 MEB::DERIV_MV_JACOBIAN_FORM );
90 outArgs.setSupports( MEB::OUT_ARG_DgDx,0,
91 MEB::DERIV_MV_GRADIENT_FORM );
92 }
93 outArgs_ = outArgs;
94 }
95
96 // Set up nominal values (initial conditions)
97 nominalValues_ = inArgs_;
98 {
99 nominalValues_.set_t(Scalar(0.0));
100 const RCP<Thyra::VectorBase<Scalar> > x_ic = createMember(x_space_);
101 { // scope to delete DetachedVectorView
102 Thyra::DetachedVectorView<Scalar> x_ic_view( *x_ic );
103 x_ic_view[0] = xIC_;
104 }
105 nominalValues_.set_x(x_ic);
106 }
107
108 if (includeXDot_) {
109 const RCP<Thyra::VectorBase<Scalar> > x_dot_ic = createMember(x_space_);
110 { // scope to delete DetachedVectorView
111 Thyra::DetachedVectorView<Scalar> x_dot_ic_view( *x_dot_ic );
112 x_dot_ic_view[0] = xDotIC_;
113 }
114 nominalValues_.set_x_dot(x_dot_ic);
115 }
116
117 isInitialized_ = true;
118}
119
120
121template<class Scalar>
122Thyra::ModelEvaluatorBase::InArgs<Scalar>
124getExactSolution(double t) const
125{
126 Thyra::ModelEvaluatorBase::InArgs<Scalar> inArgs = inArgs_;
127 double exact_t = t;
128 inArgs.set_t(exact_t);
129
130 // Set the exact solution, x.
131 Teuchos::RCP<Thyra::VectorBase<Scalar> > exact_x = createMember(x_space_);
132 { // scope to delete DetachedVectorView
133 Thyra::DetachedVectorView<Scalar> exact_x_view(*exact_x);
134 exact_x_view[0] = exp(lambda_*exact_t);
135 }
136 inArgs.set_x(exact_x);
137
138 // Set the exact solution time derivative, xDot.
139 if (includeXDot_) {
140 Teuchos::RCP<Thyra::VectorBase<Scalar> > exact_x_dot = createMember(x_space_);
141 { // scope to delete DetachedVectorView
142 Thyra::DetachedVectorView<Scalar> exact_x_dot_view(*exact_x_dot);
143 exact_x_dot_view[0] = lambda_ * exp(lambda_*exact_t);
144 }
145 inArgs.set_x_dot(exact_x_dot);
146 }
147
148 return inArgs;
149}
150
151
152template<class Scalar>
153Teuchos::RCP<const Thyra::VectorSpaceBase<Scalar> >
155get_x_space() const
156{
157 return x_space_;
158}
159
160
161template<class Scalar>
162Teuchos::RCP<const Thyra::VectorSpaceBase<Scalar> >
164get_f_space() const
165{
166 return f_space_;
167}
168
169
170template<class Scalar>
171Thyra::ModelEvaluatorBase::InArgs<Scalar>
173getNominalValues() const
174{
175 return nominalValues_;
176}
177
178
179template<class Scalar>
180Thyra::ModelEvaluatorBase::InArgs<Scalar>
182createInArgs() const
183{
184 return inArgs_;
185}
186
187
188template<class Scalar>
189Thyra::ModelEvaluatorBase::OutArgs<Scalar>
191createOutArgsImpl() const
192{
193 return outArgs_;
194}
195
196
197template<class Scalar>
198void
201 const Thyra::ModelEvaluatorBase::InArgs<Scalar> &inArgs,
202 const Thyra::ModelEvaluatorBase::OutArgs<Scalar> &outArgs
203 ) const
204{
205 using Teuchos::RCP;
206 using Thyra::VectorBase;
208 using Teuchos::rcp_dynamic_cast;
209 TEUCHOS_TEST_FOR_EXCEPTION( !isInitialized_, std::logic_error,
210 "Error, setupInOutArgs_ must be called first!\n");
211
212 const RCP<const VectorBase<Scalar> > x_in = inArgs.get_x().assert_not_null();
213 Thyra::ConstDetachedVectorView<Scalar> x_in_view( *x_in );
214 const RCP<VectorBase<Scalar> > f_out = outArgs.get_f();
215 const RCP<Thyra::LinearOpBase<Scalar> > W_out = outArgs.get_W_op();
216
217 if (inArgs.get_x_dot().is_null()) {
218
219 // Evaluate the Explicit ODE f(x,t) [= 0]
220 if (!is_null(f_out)) {
221 Thyra::DetachedVectorView<Scalar> f_out_view( *f_out );
222 f_out_view[0] = lambda_*x_in_view[0];
223 } else {
224 TEUCHOS_TEST_FOR_EXCEPTION( true, std::logic_error,
225 "Error -- Dahlquist Test Model requires f_out!\n");
226 }
227
228 if (!is_null(W_out)) {
229 RCP<Thyra::MultiVectorBase<Scalar> > matrix =
230 Teuchos::rcp_dynamic_cast<Thyra::MultiVectorBase<Scalar> >(W_out,true);
231 Thyra::DetachedMultiVectorView<Scalar> matrix_view( *matrix );
232 matrix_view(0,0) = lambda_;
233 }
234
235 } else {
236
237 // Evaluate the implicit ODE f(xdot, x, t) [=0]
238 RCP<const VectorBase<Scalar> > x_dot_in;
239 x_dot_in = inArgs.get_x_dot().assert_not_null();
240 Scalar alpha = inArgs.get_alpha();
241 Scalar beta = inArgs.get_beta();
242
243 if (!is_null(f_out)) {
244 Thyra::DetachedVectorView<Scalar> f_out_view( *f_out );
245 Thyra::ConstDetachedVectorView<Scalar> x_dot_in_view( *x_dot_in );
246 f_out_view[0] = x_dot_in_view[0] - lambda_*x_in_view[0];
247 }
248
249 if (!is_null(W_out)) {
250 RCP<Thyra::MultiVectorBase<Scalar> > matrix =
251 Teuchos::rcp_dynamic_cast<Thyra::MultiVectorBase<Scalar> >(W_out,true);
252 Thyra::DetachedMultiVectorView<Scalar> matrix_view( *matrix );
253 matrix_view(0,0) = alpha - beta*lambda_; // d(f0)/d(x0_n)
254 // Note: alpha = d(xdot)/d(x_n) and beta = d(x)/d(x_n)
255 }
256
257 }
258
259}
260
261template<class Scalar>
262Teuchos::RCP<Thyra::LinearOpWithSolveBase<Scalar> >
264create_W() const
265{
266 using Teuchos::RCP;
267 RCP<const Thyra::LinearOpWithSolveFactoryBase<Scalar> > W_factory = this->get_W_factory();
268 RCP<Thyra::LinearOpBase<Scalar> > matrix = this->create_W_op();
269 {
270 RCP<Thyra::MultiVectorBase<Scalar> > multivec = Teuchos::rcp_dynamic_cast<Thyra::MultiVectorBase<Scalar> >(matrix,true);
271 {
272 RCP<Thyra::VectorBase<Scalar> > vec = Thyra::createMember(x_space_);
273 {
274 Thyra::DetachedVectorView<Scalar> vec_view( *vec );
275 vec_view[0] = lambda_;
276 }
277 V_V(multivec->col(0).ptr(),*vec);
278 }
279 }
280 RCP<Thyra::LinearOpWithSolveBase<Scalar> > W =
281 Thyra::linearOpWithSolve<Scalar>(
282 *W_factory,
283 matrix
284 );
285 return W;
286}
287
288template<class Scalar>
289Teuchos::RCP<Thyra::LinearOpBase<Scalar> >
291create_W_op() const
292{
293 //const int dim_ = 1;
294 Teuchos::RCP<Thyra::MultiVectorBase<Scalar> > matrix = Thyra::createMembers(x_space_, dim_);
295 return(matrix);
296}
297
298
299template<class Scalar>
300Teuchos::RCP<const Thyra::LinearOpWithSolveFactoryBase<Scalar> >
302get_W_factory() const
303{
304 Teuchos::RCP<Thyra::LinearOpWithSolveFactoryBase<Scalar> > W_factory =
305 Thyra::defaultSerialDenseLinearOpWithSolveFactory<Scalar>();
306 return W_factory;
307}
308
309template<class Scalar>
310Teuchos::RCP<const Thyra::VectorSpaceBase<Scalar> >
312get_p_space(int l) const
313{
314 if (!acceptModelParams_) {
315 return Teuchos::null;
316 }
317 TEUCHOS_ASSERT_IN_RANGE_UPPER_EXCLUSIVE( l, 0, Np_ );
318 if (l == 0)
319 return p_space_;
320 else if (l == 1 || l == 2)
321 return DxDp_space_;
322 return Teuchos::null;
323}
324
325
326template<class Scalar>
327Teuchos::RCP<const Teuchos::Array<std::string> >
329get_p_names(int l) const
330{
331 if (!acceptModelParams_) {
332 return Teuchos::null;
333 }
334 TEUCHOS_ASSERT_IN_RANGE_UPPER_EXCLUSIVE( l, 0, Np_ );
335 Teuchos::RCP<Teuchos::Array<std::string> > p_strings =
336 Teuchos::rcp(new Teuchos::Array<std::string>());
337 if (l == 0) {
338 p_strings->push_back("Model Coefficient: a");
339 //p_strings->push_back("Model Coefficient: f");
340 //p_strings->push_back("Model Coefficient: L");
341 }
342 else if (l == 1)
343 p_strings->push_back("DxDp");
344 else if (l == 2)
345 p_strings->push_back("Dx_dotDp");
346 return p_strings;
347}
348
349template<class Scalar>
350Teuchos::RCP<const Thyra::VectorSpaceBase<Scalar> >
352get_g_space(int j) const
353{
354 TEUCHOS_ASSERT_IN_RANGE_UPPER_EXCLUSIVE( j, 0, Ng_ );
355 return g_space_;
356}
357
358} // namespace Tempus_Test
359#endif // TEMPUS_TEST_DAHLQUIST_TEST_MODEL_IMPL_HPP
Teuchos::RCP< const Thyra::VectorSpaceBase< Scalar > > get_g_space(int j) const
Thyra::ModelEvaluatorBase::OutArgs< Scalar > createOutArgsImpl() const
Teuchos::RCP< const Thyra::VectorSpaceBase< Scalar > > get_f_space() const
void constructDahlquistTestModel(Scalar lambda, bool includeXDot)
Teuchos::RCP< Thyra::LinearOpWithSolveBase< Scalar > > create_W() const
Teuchos::RCP< const Thyra::VectorSpaceBase< Scalar > > get_x_space() const
Teuchos::RCP< Thyra::LinearOpBase< Scalar > > create_W_op() const
Thyra::ModelEvaluatorBase::InArgs< Scalar > getExactSolution(double t) const
Exact solution.
Teuchos::RCP< const Thyra::LinearOpWithSolveFactoryBase< Scalar > > get_W_factory() const
Teuchos::RCP< const Thyra::VectorSpaceBase< Scalar > > get_p_space(int l) const
void evalModelImpl(const Thyra::ModelEvaluatorBase::InArgs< Scalar > &inArgs_bar, const Thyra::ModelEvaluatorBase::OutArgs< Scalar > &outArgs_bar) const
Thyra::ModelEvaluatorBase::InArgs< Scalar > createInArgs() const
Thyra::ModelEvaluatorBase::InArgs< Scalar > getNominalValues() const
Teuchos::RCP< const Teuchos::Array< std::string > > get_p_names(int l) const