Rythmos - Transient Integration for Differential Equations Version of the Day
Loading...
Searching...
No Matches
Rythmos_StepperHelpers_def.hpp
1//@HEADER
2// ***********************************************************************
3//
4// Rythmos Package
5// Copyright (2006) Sandia Corporation
6//
7// Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
8// license for use of this work by or on behalf of the U.S. Government.
9//
10// This library is free software; you can redistribute it and/or modify
11// it under the terms of the GNU Lesser General Public License as
12// published by the Free Software Foundation; either version 2.1 of the
13// License, or (at your option) any later version.
14//
15// This library is distributed in the hope that it will be useful, but
16// WITHOUT ANY WARRANTY; without even the implied warranty of
17// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
18// Lesser General Public License for more details.
19//
20// You should have received a copy of the GNU Lesser General Public
21// License along with this library; if not, write to the Free Software
22// Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301
23// USA
24// Questions? Contact Todd S. Coffey (tscoffe@sandia.gov)
25//
26// ***********************************************************************
27//@HEADER
28
29#ifndef RYTHMOS_STEPPER_HELPERS_DEF_HPP
30#define RYTHMOS_STEPPER_HELPERS_DEF_HPP
31
32#include "Rythmos_StepperHelpers_decl.hpp"
33#include "Rythmos_InterpolationBufferHelpers.hpp"
34#include "Rythmos_InterpolatorBaseHelpers.hpp"
35#include "Teuchos_Assert.hpp"
36#include "Thyra_AssertOp.hpp"
37#include "Thyra_VectorStdOps.hpp"
38
39
40namespace Rythmos {
41
42using Teuchos::ConstNonconstObjectContainer;
43
44template<class Scalar>
45void assertValidModel(
46 const StepperBase<Scalar>& stepper,
47 const Thyra::ModelEvaluator<Scalar>& model
48 )
49{
50
51 typedef Thyra::ModelEvaluatorBase MEB;
52
53 TEUCHOS_ASSERT(stepper.acceptsModel());
54
55 const MEB::InArgs<Scalar> inArgs = model.createInArgs();
56 const MEB::OutArgs<Scalar> outArgs = model.createOutArgs();
57
58 //TEUCHOS_ASSERT(inArgs.supports(MEB::IN_ARG_t));
59 TEUCHOS_ASSERT(inArgs.supports(MEB::IN_ARG_x));
60 TEUCHOS_ASSERT(outArgs.supports(MEB::OUT_ARG_f));
61
62 if (stepper.isImplicit()) { // implicit stepper
63 TEUCHOS_ASSERT( inArgs.supports(MEB::IN_ARG_x_dot) );
64 TEUCHOS_ASSERT( inArgs.supports(MEB::IN_ARG_alpha) );
65 TEUCHOS_ASSERT( inArgs.supports(MEB::IN_ARG_beta) );
66 TEUCHOS_ASSERT( outArgs.supports(MEB::OUT_ARG_W) );
67 }
68 //else { // explicit stepper
69 // TEUCHOS_ASSERT( !inArgs.supports(MEB::IN_ARG_x_dot) );
70 // TEUCHOS_ASSERT( !inArgs.supports(MEB::IN_ARG_alpha) );
71 // TEUCHOS_ASSERT( !inArgs.supports(MEB::IN_ARG_beta) );
72 // TEUCHOS_ASSERT( !outArgs.supports(MEB::OUT_ARG_W) );
73 //}
74
75}
76
77
78template<class Scalar>
79bool setDefaultInitialConditionFromNominalValues(
80 const Thyra::ModelEvaluator<Scalar>& model,
81 const Ptr<StepperBase<Scalar> >& stepper
82 )
83{
84
85 typedef ScalarTraits<Scalar> ST;
86 typedef Thyra::ModelEvaluatorBase MEB;
87
88 if (isInitialized(*stepper))
89 return false; // Already has an initial condition
90
91 MEB::InArgs<Scalar> initCond = model.getNominalValues();
92
93 if (!is_null(initCond.get_x())) {
94 // IC has x, we will assume that initCont.get_t() is the valid start time.
95 // Therefore, we just need to check that x_dot is also set or we will
96 // create a zero x_dot
97#ifdef HAVE_RYTHMOS_DEBUG
98 THYRA_ASSERT_VEC_SPACES( "setInitialConditionIfExists(...)",
99 *model.get_x_space(), *initCond.get_x()->space() );
100#endif
101 if (initCond.supports(MEB::IN_ARG_x_dot)) {
102 if (is_null(initCond.get_x_dot())) {
103 const RCP<Thyra::VectorBase<Scalar> > x_dot =
104 createMember(model.get_x_space());
105 assign(x_dot.ptr(), ST::zero());
106 }
107 else {
108#ifdef HAVE_RYTHMOS_DEBUG
109 THYRA_ASSERT_VEC_SPACES( "setInitialConditionIfExists(...)",
110 *model.get_x_space(), *initCond.get_x_dot()->space() );
111#endif
112 }
113 }
114 stepper->setInitialCondition(initCond);
115 return true;
116 }
117
118 // The model has not nominal values for which to set the initial
119 // conditions so wo don't do anything! The stepper will still have not
120 return false;
121
122}
123
124
125template<class Scalar>
126void restart( StepperBase<Scalar> *stepper )
127{
128#ifdef HAVE_RYTHMOS_DEBUG
129 TEUCHOS_TEST_FOR_EXCEPT(0==stepper);
130#endif // HAVE_RYTHMOS_DEBUG
131 typedef Thyra::ModelEvaluatorBase MEB;
132 const Rythmos::StepStatus<double>
133 stepStatus = stepper->getStepStatus();
134 const RCP<const Thyra::ModelEvaluator<Scalar> >
135 model = stepper->getModel();
136 // First, copy all of the model's state, including parameter values etc.
137 MEB::InArgs<double> initialCondition = model->createInArgs();
138 initialCondition.setArgs(model->getNominalValues());
139 // Set the current values of the state and time
140 RCP<const Thyra::VectorBase<double> > x, x_dot;
141 Rythmos::get_x_and_x_dot(*stepper,stepStatus.time,&x,&x_dot);
142 initialCondition.set_x(x);
143 initialCondition.set_x_dot(x_dot);
144 initialCondition.set_t(stepStatus.time);
145 // Set the new initial condition back on the stepper. This will effectively
146 // reset the stepper to think that it is starting over again (which it is).
147 stepper->setInitialCondition(initialCondition);
148}
149
150template<class Scalar>
151void eval_model_explicit(
152 const Thyra::ModelEvaluator<Scalar> &model,
153 Thyra::ModelEvaluatorBase::InArgs<Scalar> &basePoint,
154 const VectorBase<Scalar>& x_in,
155 const typename Thyra::ModelEvaluatorBase::InArgs<Scalar>::ScalarMag &t_in,
156 const Ptr<VectorBase<Scalar> >& f_out,
157 const Scalar scaled_dt,
158 const Scalar stage_point
159 )
160{
161 typedef Thyra::ModelEvaluatorBase MEB;
162 MEB::InArgs<Scalar> inArgs = model.createInArgs();
163 MEB::OutArgs<Scalar> outArgs = model.createOutArgs();
164 inArgs.setArgs(basePoint);
165 inArgs.set_x(Teuchos::rcp(&x_in,false));
166 if (inArgs.supports(MEB::IN_ARG_t)) {
167 inArgs.set_t(t_in);
168 }
169 // For model evaluators whose state function f(x, x_dot, t) describes
170 // an implicit ODE, and which accept an optional x_dot input argument,
171 // make sure the latter is set to null in order to request the evaluation
172 // of a state function corresponding to the explicit ODE formulation
173 // x_dot = f(x, t)
174 if (inArgs.supports(MEB::IN_ARG_x_dot)) {
175 inArgs.set_x_dot(Teuchos::null);
176 }
177 if (inArgs.supports(MEB::IN_ARG_step_size)) {
178 inArgs.set_step_size(scaled_dt);
179 }
180 if (inArgs.supports(MEB::IN_ARG_stage_number)) {
181 inArgs.set_stage_number(stage_point);
182 }
183 outArgs.set_f(Teuchos::rcp(&*f_out,false));
184 model.evalModel(inArgs,outArgs);
185
186 //inArgs.set_x_dot(Teuchos::null);
187}
188
189
190#ifdef HAVE_THYRA_ME_POLYNOMIAL
191
192
193template<class Scalar>
194void eval_model_explicit_poly(
195 const Thyra::ModelEvaluator<Scalar> &model,
196 Thyra::ModelEvaluatorBase::InArgs<Scalar> &basePoint,
197 const Teuchos::Polynomial< VectorBase<Scalar> > &x_poly,
198 const typename Thyra::ModelEvaluatorBase::InArgs<Scalar>::ScalarMag &t,
199 const Ptr<Teuchos::Polynomial<VectorBase<Scalar> > >& f_poly
200 )
201{
202 typedef Thyra::ModelEvaluatorBase MEB;
203 MEB::InArgs<Scalar> inArgs = model.createInArgs();
204 MEB::OutArgs<Scalar> outArgs = model.createOutArgs();
205 inArgs.setArgs(basePoint);
206 inArgs.set_x_poly(Teuchos::rcp(&x_poly,false));
207 if (inArgs.supports(MEB::IN_ARG_t)) {
208 inArgs.set_t(t);
209 }
210 outArgs.set_f_poly(Teuchos::rcp(&*f_poly,false));
211
212 model.evalModel(inArgs,outArgs);
213}
214
215
216#endif // HAVE_THYRA_ME_POLYNOMIAL
217
218
219template<class Scalar>
220void defaultGetPoints(
221 const Scalar& t_old, // required inArg
222 const Ptr<const VectorBase<Scalar> >& x_old, // optional inArg
223 const Ptr<const VectorBase<Scalar> >& xdot_old, // optional inArg
224 const Scalar& t, // required inArg
225 const Ptr<const VectorBase<Scalar> >& x, // optional inArg
226 const Ptr<const VectorBase<Scalar> >& xdot, // optional inArg
227 const Array<Scalar>& time_vec, // required inArg
228 const Ptr<Array<Teuchos::RCP<const Thyra::VectorBase<Scalar> > > >& x_vec, // optional outArg
229 const Ptr<Array<Teuchos::RCP<const Thyra::VectorBase<Scalar> > > >& xdot_vec, // optional outArg
230 const Ptr<Array<typename Teuchos::ScalarTraits<Scalar>::magnitudeType> >& accuracy_vec, // optional outArg
231 const Ptr<InterpolatorBase<Scalar> > interpolator // optional inArg (note: not const)
232 )
233{
234 typedef Teuchos::ScalarTraits<Scalar> ST;
235 assertTimePointsAreSorted(time_vec);
236 TimeRange<Scalar> tr(t_old, t);
237 TEUCHOS_ASSERT( tr.isValid() );
238 if (!is_null(x_vec)) {
239 x_vec->clear();
240 }
241 if (!is_null(xdot_vec)) {
242 xdot_vec->clear();
243 }
244 if (!is_null(accuracy_vec)) {
245 accuracy_vec->clear();
246 }
247 typename Array<Scalar>::const_iterator time_it = time_vec.begin();
248 RCP<const VectorBase<Scalar> > tmpVec;
249 RCP<const VectorBase<Scalar> > tmpVecDot;
250 for (; time_it != time_vec.end() ; time_it++) {
251 Scalar time = *time_it;
252 asssertInTimeRange(tr, time);
253 Scalar accuracy = ST::zero();
254 if (compareTimeValues(time,t_old)==0) {
255 if (!is_null(x_old)) {
256 tmpVec = x_old->clone_v();
257 }
258 if (!is_null(xdot_old)) {
259 tmpVecDot = xdot_old->clone_v();
260 }
261 } else if (compareTimeValues(time,t)==0) {
262 if (!is_null(x)) {
263 tmpVec = x->clone_v();
264 }
265 if (!is_null(xdot)) {
266 tmpVecDot = xdot->clone_v();
267 }
268 } else {
269 TEUCHOS_TEST_FOR_EXCEPTION(
270 is_null(interpolator), std::logic_error,
271 "Error, getPoints: This stepper only supports time values on the boundaries!\n"
272 );
273 // At this point, we know time != t_old, time != t, interpolator != null,
274 // and time in [t_old,t], therefore, time in (t_old,t).
275 // t_old != t at this point because otherwise it would have been caught above.
276 // Now use the interpolator to pass out the interior points
277 typename DataStore<Scalar>::DataStoreVector_t ds_nodes;
278 typename DataStore<Scalar>::DataStoreVector_t ds_out;
279 {
280 // t_old
281 DataStore<Scalar> ds;
282 ds.time = t_old;
283 ds.x = rcp(x_old.get(),false);
284 ds.xdot = rcp(xdot_old.get(),false);
285 ds_nodes.push_back(ds);
286 }
287 {
288 // t
289 DataStore<Scalar> ds;
290 ds.time = t;
291 ds.x = rcp(x.get(),false);
292 ds.xdot = rcp(xdot.get(),false);
293 ds_nodes.push_back(ds);
294 }
295 Array<Scalar> time_vec_in;
296 time_vec_in.push_back(time);
297 interpolate<Scalar>(*interpolator,rcp(&ds_nodes,false),time_vec_in,&ds_out);
298 Array<Scalar> time_vec_out;
299 Array<RCP<const VectorBase<Scalar> > > x_vec_out;
300 Array<RCP<const VectorBase<Scalar> > > xdot_vec_out;
301 Array<typename Teuchos::ScalarTraits<Scalar>::magnitudeType> accuracy_vec_out;
302 dataStoreVectorToVector(ds_out,&time_vec_out,&x_vec_out,&xdot_vec_out,&accuracy_vec_out);
303 TEUCHOS_ASSERT( time_vec_out.length()==1 );
304 tmpVec = x_vec_out[0];
305 tmpVecDot = xdot_vec_out[0];
306 accuracy = accuracy_vec_out[0];
307 }
308 if (!is_null(x_vec)) {
309 x_vec->push_back(tmpVec);
310 }
311 if (!is_null(xdot_vec)) {
312 xdot_vec->push_back(tmpVecDot);
313 }
314 if (!is_null(accuracy_vec)) {
315 accuracy_vec->push_back(accuracy);
316 }
317 tmpVec = Teuchos::null;
318 tmpVecDot = Teuchos::null;
319 }
320}
321
322
323template<class Scalar>
324 void setStepperModel(
325 const Ptr<StepperBase<Scalar> >& stepper,
326 const RCP<const Thyra::ModelEvaluator<Scalar> >& model
327 )
328{
329 stepper->setModel(model);
330}
331
332template<class Scalar>
333 void setStepperModel(
334 const Ptr<StepperBase<Scalar> >& stepper,
335 const RCP<Thyra::ModelEvaluator<Scalar> >& model
336 )
337{
338 stepper->setNonconstModel(model);
339}
340
341template<class Scalar>
342 void setStepperModel(
343 const Ptr<StepperBase<Scalar> >& stepper,
344 ConstNonconstObjectContainer<Thyra::ModelEvaluator<Scalar> >& model
345 )
346{
347 if (model.isConst()) {
348 stepper->setModel(model.getConstObj());
349 }
350 else {
351 stepper->setNonconstModel(model.getNonconstObj());
352 }
353}
354
355
356//
357// Explicit Instantiation macro
358//
359// Must be expanded from within the Rythmos namespace!
360//
361
362
363#ifdef HAVE_THYRA_ME_POLYNOMIAL
364
365#define RYTHMOS_STEPPER_HELPERS_POLY_INSTANT(SCALAR) \
366 template void eval_model_explicit_poly( \
367 const Thyra::ModelEvaluator< SCALAR > &model, \
368 Thyra::ModelEvaluatorBase::InArgs< SCALAR > &basePoint, \
369 const Teuchos::Polynomial< VectorBase< SCALAR > > &x_poly, \
370 const Thyra::ModelEvaluatorBase::InArgs< SCALAR >::ScalarMag &t, \
371 const Ptr<Teuchos::Polynomial<VectorBase< SCALAR > > >& f_poly \
372 );
373
374#else // HAVE_THYRA_ME_POLYNOMIAL
375
376#define RYTHMOS_STEPPER_HELPERS_POLY_INSTANT(SCALAR)
377
378#endif // HAVE_THYRA_ME_POLYNOMIAL
379
380
381#define RYTHMOS_STEPPER_HELPERS_INSTANT(SCALAR) \
382 \
383 template void assertValidModel( \
384 const StepperBase< SCALAR >& stepper, \
385 const Thyra::ModelEvaluator< SCALAR >& model \
386 ); \
387 template bool setDefaultInitialConditionFromNominalValues( \
388 const Thyra::ModelEvaluator< SCALAR >& model, \
389 const Ptr<StepperBase< SCALAR > >& stepper \
390 ); \
391 template void restart( StepperBase< SCALAR > *stepper ); \
392 \
393 template void eval_model_explicit( \
394 const Thyra::ModelEvaluator< SCALAR > &model, \
395 Thyra::ModelEvaluatorBase::InArgs< SCALAR > &basePoint, \
396 const VectorBase< SCALAR >& x_in, \
397 const Thyra::ModelEvaluatorBase::InArgs< SCALAR >::ScalarMag &t_in, \
398 const Ptr<VectorBase< SCALAR > >& f_out, \
399 const SCALAR scaled_dt, \
400 const SCALAR stage_point\
401 ); \
402 \
403 RYTHMOS_STEPPER_HELPERS_POLY_INSTANT(SCALAR) \
404 \
405 template void defaultGetPoints( \
406 const SCALAR & t_old, \
407 const Ptr<const VectorBase< SCALAR > >& x_old, \
408 const Ptr<const VectorBase< SCALAR > >& xdot_old, \
409 const SCALAR & t, \
410 const Ptr<const VectorBase< SCALAR > >& x, \
411 const Ptr<const VectorBase< SCALAR > >& xdot, \
412 const Array< SCALAR >& time_vec, \
413 const Ptr<Array<Teuchos::RCP<const Thyra::VectorBase< SCALAR > > > >& x_vec, \
414 const Ptr<Array<Teuchos::RCP<const Thyra::VectorBase< SCALAR > > > >& xdot_vec, \
415 const Ptr<Array<Teuchos::ScalarTraits< SCALAR >::magnitudeType> >& accuracy_vec, \
416 const Ptr<InterpolatorBase< SCALAR > > interpolator \
417 ); \
418 \
419 template void setStepperModel( \
420 const Ptr<StepperBase< SCALAR > >& stepper, \
421 const RCP<const Thyra::ModelEvaluator< SCALAR > >& model \
422 ); \
423 \
424 template void setStepperModel( \
425 const Ptr<StepperBase< SCALAR > >& stepper, \
426 const RCP<Thyra::ModelEvaluator< SCALAR > >& model \
427 ); \
428 \
429 template void setStepperModel( \
430 const Ptr<StepperBase< SCALAR > >& stepper, \
431 Teuchos::ConstNonconstObjectContainer<Thyra::ModelEvaluator< SCALAR > >& model \
432 );
433
434} // namespace Rythmos
435
436
437#endif // RYTHMOS_STEPPER_HELPERS_DEF_HPP