Rythmos - Transient Integration for Differential Equations  Version of the Day
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 
40 namespace Rythmos {
41 
42 using Teuchos::ConstNonconstObjectContainer;
43 
44 template<class Scalar>
45 void 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 
78 template<class Scalar>
79 bool 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 
125 template<class Scalar>
126 void 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;
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 
150 template<class Scalar>
151 void 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  )
158 {
159  typedef Thyra::ModelEvaluatorBase MEB;
160  MEB::InArgs<Scalar> inArgs = model.createInArgs();
161  MEB::OutArgs<Scalar> outArgs = model.createOutArgs();
162  inArgs.setArgs(basePoint);
163  inArgs.set_x(Teuchos::rcp(&x_in,false));
164  if (inArgs.supports(MEB::IN_ARG_t)) {
165  inArgs.set_t(t_in);
166  }
167  // For model evaluators whose state function f(x, x_dot, t) describes
168  // an implicit ODE, and which accept an optional x_dot input argument,
169  // make sure the latter is set to null in order to request the evaluation
170  // of a state function corresponding to the explicit ODE formulation
171  // x_dot = f(x, t)
172  if (inArgs.supports(MEB::IN_ARG_x_dot)) {
173  inArgs.set_x_dot(Teuchos::null);
174  }
175  outArgs.set_f(Teuchos::rcp(&*f_out,false));
176  model.evalModel(inArgs,outArgs);
177 }
178 
179 
180 #ifdef HAVE_THYRA_ME_POLYNOMIAL
181 
182 
183 template<class Scalar>
184 void eval_model_explicit_poly(
185  const Thyra::ModelEvaluator<Scalar> &model,
186  Thyra::ModelEvaluatorBase::InArgs<Scalar> &basePoint,
187  const Teuchos::Polynomial< VectorBase<Scalar> > &x_poly,
188  const typename Thyra::ModelEvaluatorBase::InArgs<Scalar>::ScalarMag &t,
189  const Ptr<Teuchos::Polynomial<VectorBase<Scalar> > >& f_poly
190  )
191 {
192  typedef Thyra::ModelEvaluatorBase MEB;
193  MEB::InArgs<Scalar> inArgs = model.createInArgs();
194  MEB::OutArgs<Scalar> outArgs = model.createOutArgs();
195  inArgs.setArgs(basePoint);
196  inArgs.set_x_poly(Teuchos::rcp(&x_poly,false));
197  if (inArgs.supports(MEB::IN_ARG_t)) {
198  inArgs.set_t(t);
199  }
200  outArgs.set_f_poly(Teuchos::rcp(&*f_poly,false));
201 
202  model.evalModel(inArgs,outArgs);
203 }
204 
205 
206 #endif // HAVE_THYRA_ME_POLYNOMIAL
207 
208 
209 template<class Scalar>
210 void defaultGetPoints(
211  const Scalar& t_old, // required inArg
212  const Ptr<const VectorBase<Scalar> >& x_old, // optional inArg
213  const Ptr<const VectorBase<Scalar> >& xdot_old, // optional inArg
214  const Scalar& t, // required inArg
215  const Ptr<const VectorBase<Scalar> >& x, // optional inArg
216  const Ptr<const VectorBase<Scalar> >& xdot, // optional inArg
217  const Array<Scalar>& time_vec, // required inArg
218  const Ptr<Array<Teuchos::RCP<const Thyra::VectorBase<Scalar> > > >& x_vec, // optional outArg
219  const Ptr<Array<Teuchos::RCP<const Thyra::VectorBase<Scalar> > > >& xdot_vec, // optional outArg
220  const Ptr<Array<typename Teuchos::ScalarTraits<Scalar>::magnitudeType> >& accuracy_vec, // optional outArg
221  const Ptr<InterpolatorBase<Scalar> > interpolator // optional inArg (note: not const)
222  )
223 {
224  typedef Teuchos::ScalarTraits<Scalar> ST;
225  assertTimePointsAreSorted(time_vec);
226  TimeRange<Scalar> tr(t_old, t);
227  TEUCHOS_ASSERT( tr.isValid() );
228  if (!is_null(x_vec)) {
229  x_vec->clear();
230  }
231  if (!is_null(xdot_vec)) {
232  xdot_vec->clear();
233  }
234  if (!is_null(accuracy_vec)) {
235  accuracy_vec->clear();
236  }
237  typename Array<Scalar>::const_iterator time_it = time_vec.begin();
238  RCP<const VectorBase<Scalar> > tmpVec;
239  RCP<const VectorBase<Scalar> > tmpVecDot;
240  for (; time_it != time_vec.end() ; time_it++) {
241  Scalar time = *time_it;
242  asssertInTimeRange(tr, time);
243  Scalar accuracy = ST::zero();
244  if (compareTimeValues(time,t_old)==0) {
245  if (!is_null(x_old)) {
246  tmpVec = x_old->clone_v();
247  }
248  if (!is_null(xdot_old)) {
249  tmpVecDot = xdot_old->clone_v();
250  }
251  } else if (compareTimeValues(time,t)==0) {
252  if (!is_null(x)) {
253  tmpVec = x->clone_v();
254  }
255  if (!is_null(xdot)) {
256  tmpVecDot = xdot->clone_v();
257  }
258  } else {
259  TEUCHOS_TEST_FOR_EXCEPTION(
260  is_null(interpolator), std::logic_error,
261  "Error, getPoints: This stepper only supports time values on the boundaries!\n"
262  );
263  // At this point, we know time != t_old, time != t, interpolator != null,
264  // and time in [t_old,t], therefore, time in (t_old,t).
265  // t_old != t at this point because otherwise it would have been caught above.
266  // Now use the interpolator to pass out the interior points
267  typename DataStore<Scalar>::DataStoreVector_t ds_nodes;
268  typename DataStore<Scalar>::DataStoreVector_t ds_out;
269  {
270  // t_old
271  DataStore<Scalar> ds;
272  ds.time = t_old;
273  ds.x = rcp(x_old.get(),false);
274  ds.xdot = rcp(xdot_old.get(),false);
275  ds_nodes.push_back(ds);
276  }
277  {
278  // t
279  DataStore<Scalar> ds;
280  ds.time = t;
281  ds.x = rcp(x.get(),false);
282  ds.xdot = rcp(xdot.get(),false);
283  ds_nodes.push_back(ds);
284  }
285  Array<Scalar> time_vec_in;
286  time_vec_in.push_back(time);
287  interpolate<Scalar>(*interpolator,rcp(&ds_nodes,false),time_vec_in,&ds_out);
288  Array<Scalar> time_vec_out;
289  Array<RCP<const VectorBase<Scalar> > > x_vec_out;
290  Array<RCP<const VectorBase<Scalar> > > xdot_vec_out;
291  Array<typename Teuchos::ScalarTraits<Scalar>::magnitudeType> accuracy_vec_out;
292  dataStoreVectorToVector(ds_out,&time_vec_out,&x_vec_out,&xdot_vec_out,&accuracy_vec_out);
293  TEUCHOS_ASSERT( time_vec_out.length()==1 );
294  tmpVec = x_vec_out[0];
295  tmpVecDot = xdot_vec_out[0];
296  accuracy = accuracy_vec_out[0];
297  }
298  if (!is_null(x_vec)) {
299  x_vec->push_back(tmpVec);
300  }
301  if (!is_null(xdot_vec)) {
302  xdot_vec->push_back(tmpVecDot);
303  }
304  if (!is_null(accuracy_vec)) {
305  accuracy_vec->push_back(accuracy);
306  }
307  tmpVec = Teuchos::null;
308  tmpVecDot = Teuchos::null;
309  }
310 }
311 
312 
313 template<class Scalar>
314  void setStepperModel(
315  const Ptr<StepperBase<Scalar> >& stepper,
316  const RCP<const Thyra::ModelEvaluator<Scalar> >& model
317  )
318 {
319  stepper->setModel(model);
320 }
321 
322 template<class Scalar>
323  void setStepperModel(
324  const Ptr<StepperBase<Scalar> >& stepper,
325  const RCP<Thyra::ModelEvaluator<Scalar> >& model
326  )
327 {
328  stepper->setNonconstModel(model);
329 }
330 
331 template<class Scalar>
332  void setStepperModel(
333  const Ptr<StepperBase<Scalar> >& stepper,
334  ConstNonconstObjectContainer<Thyra::ModelEvaluator<Scalar> >& model
335  )
336 {
337  if (model.isConst()) {
338  stepper->setModel(model.getConstObj());
339  }
340  else {
341  stepper->setNonconstModel(model.getNonconstObj());
342  }
343 }
344 
345 
346 //
347 // Explicit Instantiation macro
348 //
349 // Must be expanded from within the Rythmos namespace!
350 //
351 
352 
353 #ifdef HAVE_THYRA_ME_POLYNOMIAL
354 
355 #define RYTHMOS_STEPPER_HELPERS_POLY_INSTANT(SCALAR) \
356  template void eval_model_explicit_poly( \
357  const Thyra::ModelEvaluator< SCALAR > &model, \
358  Thyra::ModelEvaluatorBase::InArgs< SCALAR > &basePoint, \
359  const Teuchos::Polynomial< VectorBase< SCALAR > > &x_poly, \
360  const Thyra::ModelEvaluatorBase::InArgs< SCALAR >::ScalarMag &t, \
361  const Ptr<Teuchos::Polynomial<VectorBase< SCALAR > > >& f_poly \
362  );
363 
364 #else // HAVE_THYRA_ME_POLYNOMIAL
365 
366 #define RYTHMOS_STEPPER_HELPERS_POLY_INSTANT(SCALAR)
367 
368 #endif // HAVE_THYRA_ME_POLYNOMIAL
369 
370 
371 #define RYTHMOS_STEPPER_HELPERS_INSTANT(SCALAR) \
372  \
373  template void assertValidModel( \
374  const StepperBase< SCALAR >& stepper, \
375  const Thyra::ModelEvaluator< SCALAR >& model \
376  ); \
377  template bool setDefaultInitialConditionFromNominalValues( \
378  const Thyra::ModelEvaluator< SCALAR >& model, \
379  const Ptr<StepperBase< SCALAR > >& stepper \
380  ); \
381  template void restart( StepperBase< SCALAR > *stepper ); \
382  \
383  template void eval_model_explicit( \
384  const Thyra::ModelEvaluator< SCALAR > &model, \
385  Thyra::ModelEvaluatorBase::InArgs< SCALAR > &basePoint, \
386  const VectorBase< SCALAR >& x_in, \
387  const Thyra::ModelEvaluatorBase::InArgs< SCALAR >::ScalarMag &t_in, \
388  const Ptr<VectorBase< SCALAR > >& f_out \
389  ); \
390  \
391  RYTHMOS_STEPPER_HELPERS_POLY_INSTANT(SCALAR) \
392  \
393  template void defaultGetPoints( \
394  const SCALAR & t_old, \
395  const Ptr<const VectorBase< SCALAR > >& x_old, \
396  const Ptr<const VectorBase< SCALAR > >& xdot_old, \
397  const SCALAR & t, \
398  const Ptr<const VectorBase< SCALAR > >& x, \
399  const Ptr<const VectorBase< SCALAR > >& xdot, \
400  const Array< SCALAR >& time_vec, \
401  const Ptr<Array<Teuchos::RCP<const Thyra::VectorBase< SCALAR > > > >& x_vec, \
402  const Ptr<Array<Teuchos::RCP<const Thyra::VectorBase< SCALAR > > > >& xdot_vec, \
403  const Ptr<Array<Teuchos::ScalarTraits< SCALAR >::magnitudeType> >& accuracy_vec, \
404  const Ptr<InterpolatorBase< SCALAR > > interpolator \
405  ); \
406  \
407  template void setStepperModel( \
408  const Ptr<StepperBase< SCALAR > >& stepper, \
409  const RCP<const Thyra::ModelEvaluator< SCALAR > >& model \
410  ); \
411  \
412  template void setStepperModel( \
413  const Ptr<StepperBase< SCALAR > >& stepper, \
414  const RCP<Thyra::ModelEvaluator< SCALAR > >& model \
415  ); \
416  \
417  template void setStepperModel( \
418  const Ptr<StepperBase< SCALAR > >& stepper, \
419  Teuchos::ConstNonconstObjectContainer<Thyra::ModelEvaluator< SCALAR > >& model \
420  );
421 
422 } // namespace Rythmos
423 
424 
425 #endif // RYTHMOS_STEPPER_HELPERS_DEF_HPP