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