Tempus Version of the Day
Time Integration
Loading...
Searching...
No Matches
VanDerPol_IMEX_ExplicitModel_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_EXPLICITMODEL_IMPL_HPP
10#define TEMPUS_TEST_VANDERPOL_IMEX_EXPLICITMODEL_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#include "Thyra_DefaultBlockedLinearOp.hpp"
24
25#include <iostream>
26
27
28namespace Tempus_Test {
29
30template<class Scalar>
33 Teuchos::RCP<Teuchos::ParameterList> pList_, bool useProductVector)
34{
35 useProductVector_ = useProductVector;
36 isInitialized_ = false;
37 dim_ = 2;
38 Np_ = 3; // Number of parameter vectors (p, dx/dp, dx_dot/dp)
39 np_ = 1; // Number of parameters in this vector (1)
40 Ng_ = 0; // Number of observation functions (0)
41 ng_ = 0; // Number of elements in this observation function (0)
42 acceptModelParams_ = false;
43 haveIC_ = true;
44 epsilon_ = 1.0e-06;
45 x0_ic_ = 2.0;
46 x1_ic_ = 0.0;
47 t0_ic_ = 0.0;
48
49 if (useProductVector_ == false) {
50 x_space_ = Thyra::defaultSpmdVectorSpace<Scalar>(dim_);
51 f_space_ = Thyra::defaultSpmdVectorSpace<Scalar>(dim_);
52 } else {
53 // Create using ProductVector so we can use in partitioned IMEX-RK Stepper.
54 using Teuchos::RCP;
55 Teuchos::Array<RCP<const Thyra::VectorSpaceBase<Scalar> > > yxSpaceArray;
56 xSpace_ = Thyra::defaultSpmdVectorSpace<Scalar>(1);
57 for (int i=0; i < dim_; ++i) yxSpaceArray.push_back(xSpace_);
58 x_space_ = Thyra::productVectorSpace<Scalar>(yxSpaceArray);
59 f_space_ = Thyra::productVectorSpace<Scalar>(yxSpaceArray);
60 }
61
62 // Create p_space and g_space
63 p_space_ = Thyra::defaultSpmdVectorSpace<Scalar>(np_);
64 g_space_ = Thyra::defaultSpmdVectorSpace<Scalar>(ng_);
65 dxdp_space_ = Thyra::multiVectorProductVectorSpace(x_space_, np_);
66
67 setParameterList(pList_);
68}
69
70template<class Scalar>
71Teuchos::RCP<const Thyra::VectorSpaceBase<Scalar> >
73get_x_space() const
74{
75 return x_space_;
76}
77
78
79template<class Scalar>
80Teuchos::RCP<const Thyra::VectorSpaceBase<Scalar> >
82get_f_space() const
83{
84 return f_space_;
85}
86
87
88template<class Scalar>
89Thyra::ModelEvaluatorBase::InArgs<Scalar>
91getNominalValues() const
92{
93 TEUCHOS_TEST_FOR_EXCEPTION( !isInitialized_, std::logic_error,
94 "Error, setupInOutArgs_ must be called first!\n");
95 return nominalValues_;
96}
97
98
99template<class Scalar>
100Teuchos::RCP<Thyra::LinearOpBase<Scalar> >
102create_W_op() const
103{
104 Teuchos::RCP<Thyra::LinearOpBase<Scalar> > W_op;
105 if (useProductVector_ == false)
106 W_op = Thyra::createMembers(x_space_, dim_);
107 else {
108 // When using product vector formulation, we have to create a block
109 // operator so that its domain space is consistent with x_space_
110 Teuchos::RCP<Thyra::LinearOpBase<Scalar> > A00 =
111 Thyra::createMembers(xSpace_, 1);
112 Teuchos::RCP<Thyra::LinearOpBase<Scalar> > A01 =
113 Thyra::createMembers(xSpace_, 1);
114 Teuchos::RCP<Thyra::LinearOpBase<Scalar> > A10 =
115 Thyra::createMembers(xSpace_, 1);
116 Teuchos::RCP<Thyra::LinearOpBase<Scalar> > A11 =
117 Thyra::createMembers(xSpace_, 1);
118 W_op = Thyra::nonconstBlock2x2(A00, A01, A10, A11);
119 }
120 return W_op;
121}
122
123
124template<class Scalar>
125Thyra::ModelEvaluatorBase::InArgs<Scalar>
127createInArgs() const
128{
129 setupInOutArgs_();
130 return inArgs_;
131}
132
133
134template<class Scalar>
135Thyra::ModelEvaluatorBase::OutArgs<Scalar>
137createOutArgsImpl() const
138{
139 setupInOutArgs_();
140 return outArgs_;
141}
142
143
144template<class Scalar>
145void
147setupInOutArgs_() const
148{
149 if (isInitialized_) {
150 return;
151 }
152
153 {
154 // Set up prototypical InArgs
155 Thyra::ModelEvaluatorBase::InArgsSetup<Scalar> inArgs;
156 inArgs.setModelEvalDescription(this->description());
157 inArgs.setSupports( Thyra::ModelEvaluatorBase::IN_ARG_x_dot );
158 inArgs.setSupports( Thyra::ModelEvaluatorBase::IN_ARG_x );
159 inArgs.setSupports( Thyra::ModelEvaluatorBase::IN_ARG_t );
160 inArgs.setSupports( Thyra::ModelEvaluatorBase::IN_ARG_alpha );
161 inArgs.setSupports( Thyra::ModelEvaluatorBase::IN_ARG_beta );
162 if (acceptModelParams_) {
163 inArgs.set_Np(Np_);
164 }
165 inArgs_ = inArgs;
166 }
167
168 {
169 // Set up prototypical OutArgs
170 Thyra::ModelEvaluatorBase::OutArgsSetup<Scalar> outArgs;
171 outArgs.setModelEvalDescription(this->description());
172 outArgs.setSupports( Thyra::ModelEvaluatorBase::OUT_ARG_f );
173 outArgs.setSupports( Thyra::ModelEvaluatorBase::OUT_ARG_W_op );
174 if (acceptModelParams_) {
175 outArgs.set_Np_Ng(Np_,Ng_);
176 outArgs.setSupports( Thyra::ModelEvaluatorBase::OUT_ARG_DfDp,0,
177 Thyra::ModelEvaluatorBase::DERIV_MV_JACOBIAN_FORM );
178 }
179 outArgs_ = outArgs;
180 }
181
182 // Set up nominal values
183 nominalValues_ = inArgs_;
184 if (haveIC_)
185 {
186 using Teuchos::RCP;
187
188 nominalValues_.set_t(t0_ic_);
189 const RCP<Thyra::VectorBase<Scalar> > x_ic = createMember(x_space_);
190 { // scope to delete DetachedVectorView
191 Thyra::DetachedVectorView<Scalar> x_ic_view( *x_ic );
192 x_ic_view[0] = x0_ic_;
193 x_ic_view[1] = x1_ic_;
194 }
195 nominalValues_.set_x(x_ic);
196
197 if (acceptModelParams_) {
198 const RCP<Thyra::VectorBase<Scalar> > p_ic = createMember(p_space_);
199 {
200 Thyra::DetachedVectorView<Scalar> p_ic_view( *p_ic );
201 p_ic_view[0] = epsilon_;
202 }
203 nominalValues_.set_p(0,p_ic);
204 }
205 const RCP<Thyra::VectorBase<Scalar> > x_dot_ic = createMember(x_space_);
206 { // scope to delete DetachedVectorView
207 Thyra::DetachedVectorView<Scalar> x_dot_ic_view( *x_dot_ic );
208 x_dot_ic_view[0] = x1_ic_;
209 x_dot_ic_view[1] = ((1.0-x0_ic_*x0_ic_)*x1_ic_-x0_ic_)/epsilon_;
210 }
211 nominalValues_.set_x_dot(x_dot_ic);
212 }
213
214 isInitialized_ = true;
215}
216
217template<class Scalar>
218void
220setParameterList(Teuchos::RCP<Teuchos::ParameterList> const& paramList)
221{
222 using Teuchos::get;
223 using Teuchos::ParameterList;
224 Teuchos::RCP<ParameterList> tmpPL = Teuchos::rcp(new ParameterList("VanDerPol_IMEX_ExplicitModel"));
225 if (paramList != Teuchos::null) tmpPL = paramList;
226 tmpPL->validateParametersAndSetDefaults(*this->getValidParameters());
227 this->setMyParamList(tmpPL);
228 Teuchos::RCP<ParameterList> pl = this->getMyNonconstParamList();
229 bool acceptModelParams = get<bool>(*pl,"Accept model parameters");
230 bool haveIC = get<bool>(*pl,"Provide nominal values");
231 bool useDfDpAsTangent = get<bool>(*pl, "Use DfDp as Tangent");
232 if ( (acceptModelParams != acceptModelParams_) ||
233 (haveIC != haveIC_)
234 ) {
235 isInitialized_ = false;
236 }
237 acceptModelParams_ = acceptModelParams;
238 haveIC_ = haveIC;
239 useDfDpAsTangent_ = useDfDpAsTangent;
240 epsilon_ = get<Scalar>(*pl,"Coeff epsilon");
241 x0_ic_ = get<Scalar>(*pl,"IC x0");
242 x1_ic_ = get<Scalar>(*pl,"IC x1");
243 t0_ic_ = get<Scalar>(*pl,"IC t0");
244 setupInOutArgs_();
245}
246
247template<class Scalar>
248Teuchos::RCP<const Teuchos::ParameterList>
250getValidParameters() const
251{
252 static Teuchos::RCP<const Teuchos::ParameterList> validPL;
253 if (is_null(validPL)) {
254 Teuchos::RCP<Teuchos::ParameterList> pl = Teuchos::parameterList();
255 pl->set("Accept model parameters", false);
256 pl->set("Provide nominal values", true);
257 pl->set("Use DfDp as Tangent", false);
258 Teuchos::setDoubleParameter(
259 "Coeff epsilon", 1.0e-06, "Coefficient a in model", &*pl);
260 Teuchos::setDoubleParameter(
261 "IC x0", 2.0, "Initial Condition for x0", &*pl);
262 Teuchos::setDoubleParameter(
263 "IC x1", 0.0, "Initial Condition for x1", &*pl);
264 Teuchos::setDoubleParameter(
265 "IC t0", 0.0, "Initial time t0", &*pl);
266 validPL = pl;
267 }
268 return validPL;
269}
270
271template<class Scalar>
272Teuchos::RCP<const Thyra::VectorSpaceBase<Scalar> >
274get_p_space(int l) const
275{
276 if (!acceptModelParams_) {
277 return Teuchos::null;
278 }
279 TEUCHOS_ASSERT_IN_RANGE_UPPER_EXCLUSIVE( l, 0, Np_ );
280 if (l == 0)
281 return p_space_;
282 else if (l == 1 || l == 2)
283 return dxdp_space_;
284 return Teuchos::null;
285}
286
287template<class Scalar>
288Teuchos::RCP<const Teuchos::Array<std::string> >
290get_p_names(int l) const
291{
292 if (!acceptModelParams_) {
293 return Teuchos::null;
294 }
295 TEUCHOS_ASSERT_IN_RANGE_UPPER_EXCLUSIVE( l, 0, Np_ );
296 Teuchos::RCP<Teuchos::Array<std::string> > p_strings =
297 Teuchos::rcp(new Teuchos::Array<std::string>());
298 if (l == 0)
299 p_strings->push_back("Model Coefficient: epsilon");
300 else if (l == 1)
301 p_strings->push_back("DxDp");
302 else if (l == 2)
303 p_strings->push_back("Dx_dotDp");
304 return p_strings;
305}
306
307template<class Scalar>
308Teuchos::RCP<const Thyra::VectorSpaceBase<Scalar> >
310get_g_space(int j) const
311{
312 TEUCHOS_ASSERT_IN_RANGE_UPPER_EXCLUSIVE( j, 0, Ng_ );
313 return g_space_;
314}
315
316template<class Scalar>
317void
319evalModelImpl(const Thyra::ModelEvaluatorBase::InArgs<Scalar> &inArgs,
320 const Thyra::ModelEvaluatorBase::OutArgs<Scalar> &outArgs) const
321{
322 typedef Thyra::DefaultMultiVectorProductVector<Scalar> DMVPV;
323 using Teuchos::RCP;
324 using Teuchos::rcp_dynamic_cast;
325 TEUCHOS_TEST_FOR_EXCEPTION( !isInitialized_, std::logic_error,
326 "Error, setupInOutArgs_ must be called first!\n");
327
328 const RCP<const Thyra::VectorBase<Scalar> > x_in =
329 inArgs.get_x().assert_not_null();
330 Thyra::ConstDetachedVectorView<Scalar> x_in_view( *x_in );
331
332 //double t = inArgs.get_t();
333 Scalar beta = inArgs.get_beta();
334 Scalar eps = epsilon_;
335 RCP<const Thyra::MultiVectorBase<Scalar> > dxdp_in, dxdotdp_in;
336 if (acceptModelParams_) {
337 const RCP<const Thyra::VectorBase<Scalar> > p_in =
338 inArgs.get_p(0);
339 if (p_in != Teuchos::null) {
340 Thyra::ConstDetachedVectorView<Scalar> p_in_view( *p_in );
341 eps = p_in_view[0];
342 }
343 if (inArgs.get_p(1) != Teuchos::null) {
344 dxdp_in =
345 rcp_dynamic_cast<const DMVPV>(inArgs.get_p(1),true)->getMultiVector();
346 }
347 if (inArgs.get_p(2) != Teuchos::null) {
348 dxdotdp_in =
349 rcp_dynamic_cast<const DMVPV>(inArgs.get_p(2),true)->getMultiVector();
350 }
351 }
352
353 const RCP<Thyra::VectorBase<Scalar> > f_out = outArgs.get_f();
354 const RCP<Thyra::LinearOpBase<Scalar> > W_out = outArgs.get_W_op();
355 RCP<Thyra::MultiVectorBase<Scalar> > DfDp_out;
356 if (acceptModelParams_) {
357 Thyra::ModelEvaluatorBase::Derivative<Scalar> DfDp = outArgs.get_DfDp(0);
358 DfDp_out = DfDp.getMultiVector();
359 }
360
361 if (inArgs.get_x_dot().is_null()) {
362
363 // Evaluate the Explicit ODE f(x,t) [= xdot]
364 if (!is_null(f_out)) {
365 Thyra::DetachedVectorView<Scalar> f_out_view( *f_out );
366 f_out_view[0] = x_in_view[1];
367 f_out_view[1] = -x_in_view[0]/eps;
368 }
369 if (!is_null(DfDp_out)) {
370 Thyra::DetachedMultiVectorView<Scalar> DfDp_out_view( *DfDp_out );
371 DfDp_out_view(0,0) = 0.0;
372 DfDp_out_view(1,0) = x_in_view[0]/(eps*eps);
373
374 // Compute df/dp + (df/dx) * (dx/dp)
375 if (useDfDpAsTangent_ && !is_null(dxdp_in)) {
376 Thyra::ConstDetachedMultiVectorView<Scalar> dxdp( *dxdp_in );
377 DfDp_out_view(0,0) += dxdp(1,0);
378 DfDp_out_view(1,0) += -dxdp(0,0)/eps;
379 }
380 }
381 if (!is_null(W_out)) {
382 if (useProductVector_ == false) {
383 RCP<Thyra::MultiVectorBase<Scalar> > W =
384 Teuchos::rcp_dynamic_cast<Thyra::MultiVectorBase<Scalar> >(W_out,true);
385 Thyra::DetachedMultiVectorView<Scalar> W_view( *W );
386 W_view(0,0) = 0.0; // d(f0)/d(x0_n)
387 W_view(0,1) = beta; // d(f0)/d(x1_n)
388 W_view(1,0) = -beta/eps; // d(f1)/d(x0_n)
389 W_view(1,1) = 0.0; // d(f1)/d(x1_n)
390 // Note: alpha = d(xdot)/d(x_n) and beta = d(x)/d(x_n)
391 }
392 else {
393 RCP<Thyra::DefaultBlockedLinearOp<Scalar> > W_block =
394 Teuchos::rcp_dynamic_cast<Thyra::DefaultBlockedLinearOp<Scalar> >(W_out, true);
395 RCP<Thyra::MultiVectorBase<Scalar> > W00 =
396 Teuchos::rcp_dynamic_cast<Thyra::MultiVectorBase<Scalar> >(W_block->getNonconstBlock(0,0),true);
397 RCP<Thyra::MultiVectorBase<Scalar> > W01 =
398 Teuchos::rcp_dynamic_cast<Thyra::MultiVectorBase<Scalar> >(W_block->getNonconstBlock(0,1),true);
399 RCP<Thyra::MultiVectorBase<Scalar> > W10 =
400 Teuchos::rcp_dynamic_cast<Thyra::MultiVectorBase<Scalar> >(W_block->getNonconstBlock(1,0),true);
401 RCP<Thyra::MultiVectorBase<Scalar> > W11 =
402 Teuchos::rcp_dynamic_cast<Thyra::MultiVectorBase<Scalar> >(W_block->getNonconstBlock(1,1),true);
403 Thyra::DetachedMultiVectorView<Scalar> W00_view( *W00 );
404 Thyra::DetachedMultiVectorView<Scalar> W01_view( *W01 );
405 Thyra::DetachedMultiVectorView<Scalar> W10_view( *W10 );
406 Thyra::DetachedMultiVectorView<Scalar> W11_view( *W11 );
407 W00_view(0,0) = 0.0; // d(f0)/d(x0_n)
408 W01_view(0,0) = beta; // d(f0)/d(x1_n)
409 W10_view(0,0) = -beta/eps; // d(f1)/d(x0_n)
410 W11_view(0,0) = 0.0; // d(f1)/d(x1_n)
411 // Note: alpha = d(xdot)/d(x_n) and beta = d(x)/d(x_n)
412 }
413 }
414 } else {
415
416 // Evaluate the implicit ODE f(xdot, x, t) [= 0]
417 RCP<const Thyra::VectorBase<Scalar> > x_dot_in;
418 x_dot_in = inArgs.get_x_dot().assert_not_null();
419 Scalar alpha = inArgs.get_alpha();
420
421 if (!is_null(f_out)) {
422 Thyra::DetachedVectorView<Scalar> f_out_view( *f_out );
423 Thyra::ConstDetachedVectorView<Scalar> x_dot_in_view( *x_dot_in );
424 f_out_view[0] = x_dot_in_view[0] - x_in_view[1];
425 f_out_view[1] = x_dot_in_view[1] + x_in_view[0]/eps;
426 }
427 if (!is_null(DfDp_out)) {
428 Thyra::DetachedMultiVectorView<Scalar> DfDp_out_view( *DfDp_out );
429 DfDp_out_view(0,0) = 0.0;
430 DfDp_out_view(1,0) = -x_in_view[0]/(eps*eps);
431
432 // Compute df/dp + (df/dx_dot)*(dx_dot/dp) + (df/dx) * (dx/dp)
433 if (useDfDpAsTangent_ && !is_null(dxdotdp_in)) {
434 Thyra::ConstDetachedMultiVectorView<Scalar> dxdotdp( *dxdotdp_in );
435 DfDp_out_view(0,0) += dxdotdp(0,0);
436 DfDp_out_view(1,0) += dxdotdp(1,0);
437 }
438 if (useDfDpAsTangent_ && !is_null(dxdp_in)) {
439 Thyra::ConstDetachedMultiVectorView<Scalar> dxdp( *dxdp_in );
440 DfDp_out_view(0,0) += -dxdp(1,0);
441 DfDp_out_view(1,0) += dxdp(0,0)/eps;
442 }
443 }
444 if (!is_null(W_out)) {
445 if (useProductVector_ == false) {
446 RCP<Thyra::MultiVectorBase<Scalar> > W =
447 Teuchos::rcp_dynamic_cast<Thyra::MultiVectorBase<Scalar> >(W_out,true);
448 Thyra::DetachedMultiVectorView<Scalar> W_view( *W );
449 W_view(0,0) = alpha; // d(f0)/d(x0_n)
450 W_view(0,1) = -beta; // d(f0)/d(x1_n)
451 W_view(1,0) = beta/eps; // d(f1)/d(x0_n)
452 W_view(1,1) = alpha; // d(f1)/d(x1_n)
453 // Note: alpha = d(xdot)/d(x_n) and beta = d(x)/d(x_n)
454 }
455 else {
456 RCP<Thyra::DefaultBlockedLinearOp<Scalar> > W_block =
457 Teuchos::rcp_dynamic_cast<Thyra::DefaultBlockedLinearOp<Scalar> >(W_out, true);
458 RCP<Thyra::MultiVectorBase<Scalar> > W00 =
459 Teuchos::rcp_dynamic_cast<Thyra::MultiVectorBase<Scalar> >(W_block->getNonconstBlock(0,0),true);
460 RCP<Thyra::MultiVectorBase<Scalar> > W01 =
461 Teuchos::rcp_dynamic_cast<Thyra::MultiVectorBase<Scalar> >(W_block->getNonconstBlock(0,1),true);
462 RCP<Thyra::MultiVectorBase<Scalar> > W10 =
463 Teuchos::rcp_dynamic_cast<Thyra::MultiVectorBase<Scalar> >(W_block->getNonconstBlock(1,0),true);
464 RCP<Thyra::MultiVectorBase<Scalar> > W11 =
465 Teuchos::rcp_dynamic_cast<Thyra::MultiVectorBase<Scalar> >(W_block->getNonconstBlock(1,1),true);
466 Thyra::DetachedMultiVectorView<Scalar> W00_view( *W00 );
467 Thyra::DetachedMultiVectorView<Scalar> W01_view( *W01 );
468 Thyra::DetachedMultiVectorView<Scalar> W10_view( *W10 );
469 Thyra::DetachedMultiVectorView<Scalar> W11_view( *W11 );
470 W00_view(0,0) = alpha; // d(f0)/d(x0_n)
471 W01_view(0,0) = -beta; // d(f0)/d(x1_n)
472 W10_view(0,0) = beta/eps; // d(f1)/d(x0_n)
473 W11_view(0,0) = alpha; // d(f1)/d(x1_n)
474 // Note: alpha = d(xdot)/d(x_n) and beta = d(x)/d(x_n)
475 }
476 }
477 }
478}
479
480
481} // namespace Tempus_Test
482#endif // TEMPUS_TEST_VANDERPOL_IMEX_EXPLICITMODEL_IMPL_HPP
Teuchos::RCP< const Teuchos::ParameterList > getValidParameters() const
Teuchos::RCP< Thyra::LinearOpBase< Scalar > > create_W_op() const
void setParameterList(Teuchos::RCP< Teuchos::ParameterList > const &paramList)
void evalModelImpl(const Thyra::ModelEvaluatorBase::InArgs< Scalar > &inArgs_bar, const Thyra::ModelEvaluatorBase::OutArgs< Scalar > &outArgs_bar) const
Teuchos::RCP< const Teuchos::Array< std::string > > get_p_names(int l) const
Teuchos::RCP< const Thyra::VectorSpaceBase< Scalar > > get_g_space(int j) const
Teuchos::RCP< const Thyra::VectorSpaceBase< Scalar > > get_x_space() const
Thyra::ModelEvaluatorBase::OutArgs< Scalar > createOutArgsImpl() const
Teuchos::RCP< const Thyra::VectorSpaceBase< Scalar > > get_f_space() const
VanDerPol_IMEX_ExplicitModel(Teuchos::RCP< Teuchos::ParameterList > pList=Teuchos::null, bool useProductVector=false)
Thyra::ModelEvaluatorBase::InArgs< Scalar > createInArgs() const
Teuchos::RCP< const Thyra::VectorSpaceBase< Scalar > > get_p_space(int l) const
Thyra::ModelEvaluatorBase::InArgs< Scalar > getNominalValues() const