Tempus Version of the Day
Time Integration
Loading...
Searching...
No Matches
VanDerPolModel_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_VANDERPOL_MODEL_IMPL_HPP
10#define TEMPUS_TEST_VANDERPOL_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
22#include <iostream>
23
24
25namespace Tempus_Test {
26
27template<class Scalar>
29VanDerPolModel(Teuchos::RCP<Teuchos::ParameterList> pList_)
30{
31 isInitialized_ = false;
32 dim_ = 2;
33 Np_ = 1; // Number of parameter vectors (1)
34 np_ = 1; // Number of parameters in this vector (1)
35 Ng_ = 0; // Number of observation functions (0)
36 ng_ = 0; // Number of elements in this observation function (0)
37 acceptModelParams_ = false;
38 haveIC_ = true;
39 epsilon_ = 1.0e-06;
40 x0_ic_ = 2.0;
41 x1_ic_ = 0.0;
42 t0_ic_ = 0.0;
43
44 // Create x_space and f_space
45 x_space_ = Thyra::defaultSpmdVectorSpace<Scalar>(dim_);
46 f_space_ = Thyra::defaultSpmdVectorSpace<Scalar>(dim_);
47 // Create p_space and g_space
48 p_space_ = Thyra::defaultSpmdVectorSpace<Scalar>(np_);
49 g_space_ = Thyra::defaultSpmdVectorSpace<Scalar>(ng_);
50
51 setParameterList(pList_);
52}
53
54template<class Scalar>
55Thyra::ModelEvaluatorBase::InArgs<Scalar>
57getExactSolution(double t) const
58{
59 TEUCHOS_TEST_FOR_EXCEPTION( true, std::logic_error,
60 "Error - No exact solution for van der Pol problem!\n");
61 return(inArgs_);
62}
63
64template<class Scalar>
65Thyra::ModelEvaluatorBase::InArgs<Scalar>
67getExactSensSolution(int j, double t) const
68{
69 TEUCHOS_TEST_FOR_EXCEPTION( !isInitialized_, std::logic_error,
70 "Error - No exact sensitivities for van der Pol problem!\n");
71 return(inArgs_);
72}
73
74template<class Scalar>
75Teuchos::RCP<const Thyra::VectorSpaceBase<Scalar> >
77get_x_space() const
78{
79 return x_space_;
80}
81
82
83template<class Scalar>
84Teuchos::RCP<const Thyra::VectorSpaceBase<Scalar> >
86get_f_space() const
87{
88 return f_space_;
89}
90
91
92template<class Scalar>
93Thyra::ModelEvaluatorBase::InArgs<Scalar>
95getNominalValues() const
96{
97 TEUCHOS_TEST_FOR_EXCEPTION( !isInitialized_, std::logic_error,
98 "Error, setupInOutArgs_ must be called first!\n");
99 return nominalValues_;
100}
101
102
103template<class Scalar>
104Teuchos::RCP<Thyra::LinearOpWithSolveBase<Scalar> >
106create_W() const
107{
108 using Teuchos::RCP;
109 RCP<const Thyra::LinearOpWithSolveFactoryBase<Scalar> > W_factory = this->get_W_factory();
110 RCP<Thyra::LinearOpBase<Scalar> > matrix = this->create_W_op();
111 {
112 // 01/20/09 tscoffe: This is a total hack to provide a full rank matrix to
113 // linearOpWithSolve because it ends up factoring the matrix during
114 // initialization, which it really shouldn't do, or I'm doing something
115 // wrong here. The net effect is that I get exceptions thrown in
116 // optimized mode due to the matrix being rank deficient unless I do this.
117 RCP<Thyra::MultiVectorBase<Scalar> > multivec = Teuchos::rcp_dynamic_cast<Thyra::MultiVectorBase<Scalar> >(matrix,true);
118 {
119 RCP<Thyra::VectorBase<Scalar> > vec = Thyra::createMember(x_space_);
120 {
121 Thyra::DetachedVectorView<Scalar> vec_view( *vec );
122 vec_view[0] = 0.0;
123 vec_view[1] = 1.0;
124 }
125 V_V(multivec->col(0).ptr(),*vec);
126 {
127 Thyra::DetachedVectorView<Scalar> vec_view( *vec );
128 vec_view[0] = 1.0;
129 vec_view[1] = 0.0;
130 }
131 V_V(multivec->col(1).ptr(),*vec);
132 }
133 }
134 RCP<Thyra::LinearOpWithSolveBase<Scalar> > W =
135 Thyra::linearOpWithSolve<Scalar>(*W_factory, matrix);
136
137 return W;
138}
139
140
141template<class Scalar>
142Teuchos::RCP<Thyra::LinearOpBase<Scalar> >
144create_W_op() const
145{
146 Teuchos::RCP<Thyra::MultiVectorBase<Scalar> > matrix = Thyra::createMembers(x_space_, dim_);
147 return(matrix);
148}
149
150
151template<class Scalar>
152Teuchos::RCP<const Thyra::LinearOpWithSolveFactoryBase<Scalar> >
154get_W_factory() const
155{
156 Teuchos::RCP<Thyra::LinearOpWithSolveFactoryBase<Scalar> > W_factory =
157 Thyra::defaultSerialDenseLinearOpWithSolveFactory<Scalar>();
158 return W_factory;
159}
160
161
162template<class Scalar>
163Thyra::ModelEvaluatorBase::InArgs<Scalar>
165createInArgs() const
166{
167 setupInOutArgs_();
168 return inArgs_;
169}
170
171
172// Private functions overridden from ModelEvaluatorDefaultBase
173
174
175template<class Scalar>
176Thyra::ModelEvaluatorBase::OutArgs<Scalar>
178createOutArgsImpl() const
179{
180 setupInOutArgs_();
181 return outArgs_;
182}
183
184
185template<class Scalar>
186void
189 const Thyra::ModelEvaluatorBase::InArgs<Scalar> &inArgs,
190 const Thyra::ModelEvaluatorBase::OutArgs<Scalar> &outArgs
191 ) const
192{
193 using Teuchos::RCP;
194 TEUCHOS_TEST_FOR_EXCEPTION( !isInitialized_, std::logic_error,
195 "Error, setupInOutArgs_ must be called first!\n");
196
197 const RCP<const Thyra::VectorBase<Scalar> > x_in =
198 inArgs.get_x().assert_not_null();
199 Thyra::ConstDetachedVectorView<Scalar> x_in_view( *x_in );
200
201 //double t = inArgs.get_t();
202 Scalar epsilon = epsilon_;
203 if (acceptModelParams_) {
204 const RCP<const Thyra::VectorBase<Scalar> > p_in =
205 inArgs.get_p(0).assert_not_null();
206 Thyra::ConstDetachedVectorView<Scalar> p_in_view( *p_in );
207 epsilon = p_in_view[0];
208 }
209
210 Scalar beta = inArgs.get_beta();
211
212 const RCP<Thyra::VectorBase<Scalar> > f_out = outArgs.get_f();
213 const RCP<Thyra::LinearOpBase<Scalar> > W_out = outArgs.get_W_op();
214 RCP<Thyra::MultiVectorBase<Scalar> > DfDp_out;
215 if (acceptModelParams_) {
216 Thyra::ModelEvaluatorBase::Derivative<Scalar> DfDp = outArgs.get_DfDp(0);
217 DfDp_out = DfDp.getMultiVector();
218 }
219
220 if (inArgs.get_x_dot().is_null()) {
221
222 // Evaluate the Explicit ODE f(x,t) [= xdot]
223 if (!is_null(f_out)) {
224 Thyra::DetachedVectorView<Scalar> f_out_view( *f_out );
225 f_out_view[0] = x_in_view[1];
226 f_out_view[1] =
227 ((1.0-x_in_view[0]*x_in_view[0])*x_in_view[1]-x_in_view[0])/epsilon;
228 }
229 if (!is_null(W_out)) {
230 RCP<Thyra::MultiVectorBase<Scalar> > matrix =
231 Teuchos::rcp_dynamic_cast<Thyra::MultiVectorBase<Scalar> >(W_out,true);
232 Thyra::DetachedMultiVectorView<Scalar> matrix_view( *matrix );
233 matrix_view(0,0) = 0.0; // d(f0)/d(x0_n)
234 matrix_view(0,1) = +beta; // d(f0)/d(x1_n)
235 matrix_view(1,0) =
236 -beta*(2.0*x_in_view[0]*x_in_view[1]+1.0)/epsilon; // d(f1)/d(x0_n)
237 matrix_view(1,1) =
238 beta*(1.0 - x_in_view[0]*x_in_view[0])/epsilon; // d(f1)/d(x1_n)
239 // Note: alpha = d(xdot)/d(x_n) and beta = d(x)/d(x_n)
240 }
241 if (!is_null(DfDp_out)) {
242 Thyra::DetachedMultiVectorView<Scalar> DfDp_out_view( *DfDp_out );
243 DfDp_out_view(0,0) = 0.0;
244 DfDp_out_view(1,0) =
245 -((1.0-x_in_view[0]*x_in_view[0])*x_in_view[1]-x_in_view[0])
246 /(epsilon*epsilon);
247 }
248 } else {
249
250 // Evaluate the implicit ODE f(xdot, x, t) [= 0]
251 RCP<const Thyra::VectorBase<Scalar> > x_dot_in;
252 x_dot_in = inArgs.get_x_dot().assert_not_null();
253 Scalar alpha = inArgs.get_alpha();
254 if (!is_null(f_out)) {
255 Thyra::DetachedVectorView<Scalar> f_out_view( *f_out );
256 Thyra::ConstDetachedVectorView<Scalar> x_dot_in_view( *x_dot_in );
257 f_out_view[0] = x_dot_in_view[0] - x_in_view[1];
258 f_out_view[1] = x_dot_in_view[1]
259 - ((1.0-x_in_view[0]*x_in_view[0])*x_in_view[1]-x_in_view[0])/epsilon;
260 }
261 if (!is_null(W_out)) {
262 RCP<Thyra::MultiVectorBase<Scalar> > matrix =
263 Teuchos::rcp_dynamic_cast<Thyra::MultiVectorBase<Scalar> >(W_out,true);
264 Thyra::DetachedMultiVectorView<Scalar> matrix_view( *matrix );
265 matrix_view(0,0) = alpha; // d(f0)/d(x0_n)
266 matrix_view(0,1) = -beta; // d(f0)/d(x1_n)
267 matrix_view(1,0) =
268 beta*(2.0*x_in_view[0]*x_in_view[1]+1.0)/epsilon; // d(f1)/d(x0_n)
269 matrix_view(1,1) = alpha
270 - beta*(1.0 - x_in_view[0]*x_in_view[0])/epsilon; // d(f1)/d(x1_n)
271 // Note: alpha = d(xdot)/d(x_n) and beta = d(x)/d(x_n)
272 }
273 if (!is_null(DfDp_out)) {
274 Thyra::DetachedMultiVectorView<Scalar> DfDp_out_view( *DfDp_out );
275 DfDp_out_view(0,0) = 0.0;
276 DfDp_out_view(1,0) =
277 ((1.0-x_in_view[0]*x_in_view[0])*x_in_view[1]-x_in_view[0])
278 /(epsilon*epsilon);
279 }
280 }
281}
282
283template<class Scalar>
284Teuchos::RCP<const Thyra::VectorSpaceBase<Scalar> >
286get_p_space(int l) const
287{
288 if (!acceptModelParams_) {
289 return Teuchos::null;
290 }
291 TEUCHOS_ASSERT_IN_RANGE_UPPER_EXCLUSIVE( l, 0, Np_ );
292 return p_space_;
293}
294
295template<class Scalar>
296Teuchos::RCP<const Teuchos::Array<std::string> >
298get_p_names(int l) const
299{
300 if (!acceptModelParams_) {
301 return Teuchos::null;
302 }
303 TEUCHOS_ASSERT_IN_RANGE_UPPER_EXCLUSIVE( l, 0, Np_ );
304 Teuchos::RCP<Teuchos::Array<std::string> > p_strings =
305 Teuchos::rcp(new Teuchos::Array<std::string>());
306 p_strings->push_back("Model Coefficient: epsilon");
307 return p_strings;
308}
309
310template<class Scalar>
311Teuchos::RCP<const Thyra::VectorSpaceBase<Scalar> >
313get_g_space(int j) const
314{
315 TEUCHOS_ASSERT_IN_RANGE_UPPER_EXCLUSIVE( j, 0, Ng_ );
316 return g_space_;
317}
318
319// private
320
321template<class Scalar>
322void
324setupInOutArgs_() const
325{
326 if (isInitialized_) {
327 return;
328 }
329
330 {
331 // Set up prototypical InArgs
332 Thyra::ModelEvaluatorBase::InArgsSetup<Scalar> inArgs;
333 inArgs.setModelEvalDescription(this->description());
334 inArgs.setSupports( Thyra::ModelEvaluatorBase::IN_ARG_t );
335 inArgs.setSupports( Thyra::ModelEvaluatorBase::IN_ARG_x );
336 inArgs.setSupports( Thyra::ModelEvaluatorBase::IN_ARG_beta );
337 inArgs.setSupports( Thyra::ModelEvaluatorBase::IN_ARG_x_dot );
338 inArgs.setSupports( Thyra::ModelEvaluatorBase::IN_ARG_alpha );
339 if (acceptModelParams_) {
340 inArgs.set_Np(Np_);
341 }
342 inArgs_ = inArgs;
343 }
344
345 {
346 // Set up prototypical OutArgs
347 Thyra::ModelEvaluatorBase::OutArgsSetup<Scalar> outArgs;
348 outArgs.setModelEvalDescription(this->description());
349 outArgs.setSupports( Thyra::ModelEvaluatorBase::OUT_ARG_f );
350 outArgs.setSupports( Thyra::ModelEvaluatorBase::OUT_ARG_W_op );
351 if (acceptModelParams_) {
352 outArgs.set_Np_Ng(Np_,Ng_);
353 outArgs.setSupports( Thyra::ModelEvaluatorBase::OUT_ARG_DfDp,0,
354 Thyra::ModelEvaluatorBase::DERIV_MV_JACOBIAN_FORM );
355 }
356 outArgs_ = outArgs;
357 }
358
359 // Set up nominal values
360 nominalValues_ = inArgs_;
361 if (haveIC_)
362 {
363 using Teuchos::RCP;
364 nominalValues_.set_t(t0_ic_);
365 const RCP<Thyra::VectorBase<Scalar> > x_ic = createMember(x_space_);
366 { // scope to delete DetachedVectorView
367 Thyra::DetachedVectorView<Scalar> x_ic_view( *x_ic );
368 x_ic_view[0] = x0_ic_;
369 x_ic_view[1] = x1_ic_;
370 }
371 nominalValues_.set_x(x_ic);
372 if (acceptModelParams_) {
373 const RCP<Thyra::VectorBase<Scalar> > p_ic = createMember(p_space_);
374 {
375 Thyra::DetachedVectorView<Scalar> p_ic_view( *p_ic );
376 p_ic_view[0] = epsilon_;
377 }
378 nominalValues_.set_p(0,p_ic);
379 }
380 const RCP<Thyra::VectorBase<Scalar> > x_dot_ic = createMember(x_space_);
381 { // scope to delete DetachedVectorView
382 Thyra::DetachedVectorView<Scalar> x_dot_ic_view( *x_dot_ic );
383 x_dot_ic_view[0] = x1_ic_;
384 x_dot_ic_view[1] = ((1.0-x0_ic_*x0_ic_)*x1_ic_-x0_ic_)/epsilon_;
385 }
386 nominalValues_.set_x_dot(x_dot_ic);
387 }
388
389 isInitialized_ = true;
390
391}
392
393template<class Scalar>
394void
396setParameterList(Teuchos::RCP<Teuchos::ParameterList> const& paramList)
397{
398 using Teuchos::get;
399 using Teuchos::ParameterList;
400 Teuchos::RCP<ParameterList> tmpPL = Teuchos::rcp(new ParameterList("VanDerPolModel"));
401 if (paramList != Teuchos::null) tmpPL = paramList;
402 tmpPL->validateParametersAndSetDefaults(*this->getValidParameters());
403 this->setMyParamList(tmpPL);
404 Teuchos::RCP<ParameterList> pl = this->getMyNonconstParamList();
405 bool acceptModelParams = get<bool>(*pl,"Accept model parameters");
406 bool haveIC = get<bool>(*pl,"Provide nominal values");
407 if ( (acceptModelParams != acceptModelParams_) ||
408 (haveIC != haveIC_)
409 ) {
410 isInitialized_ = false;
411 }
412 acceptModelParams_ = acceptModelParams;
413 haveIC_ = haveIC;
414 epsilon_ = get<Scalar>(*pl,"Coeff epsilon");
415 x0_ic_ = get<Scalar>(*pl,"IC x0");
416 x1_ic_ = get<Scalar>(*pl,"IC x1");
417 t0_ic_ = get<Scalar>(*pl,"IC t0");
418 setupInOutArgs_();
419}
420
421template<class Scalar>
422Teuchos::RCP<const Teuchos::ParameterList>
424getValidParameters() const
425{
426 static Teuchos::RCP<const Teuchos::ParameterList> validPL;
427 if (is_null(validPL)) {
428 Teuchos::RCP<Teuchos::ParameterList> pl = Teuchos::parameterList();
429 pl->set("Accept model parameters", false);
430 pl->set("Provide nominal values", true);
431 Teuchos::setDoubleParameter(
432 "Coeff epsilon", 1.0e-06, "Coefficient a in model", &*pl);
433 Teuchos::setDoubleParameter(
434 "IC x0", 2.0, "Initial Condition for x0", &*pl);
435 Teuchos::setDoubleParameter(
436 "IC x1", 0.0, "Initial Condition for x1", &*pl);
437 Teuchos::setDoubleParameter(
438 "IC t0", 0.0, "Initial time t0", &*pl);
439 Teuchos::setIntParameter(
440 "Number of Time Step Sizes", 1, "Number time step sizes for convergence study", &*pl);
441 validPL = pl;
442 }
443 return validPL;
444}
445
446} // namespace Tempus_Test
447#endif // TEMPUS_TEST_VANDERPOL_MODEL_IMPL_HPP
van der Pol model problem for nonlinear electrical circuit.
VanDerPolModel(Teuchos::RCP< Teuchos::ParameterList > pList=Teuchos::null)