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