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