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