Rythmos - Transient Integration for Differential Equations  Version of the Day
 All Classes Functions Variables Typedefs Pages
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;
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 
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  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 
193 template<class Scalar>
194 void 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 
219 template<class Scalar>
220 void 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 
323 template<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 
332 template<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 
341 template<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
Generated by   doxygen 1.8.5