Rythmos - Transient Integration for Differential Equations Version of the Day
Rythmos_StepperHelpers_def.hpp
00001 //@HEADER
00002 // ***********************************************************************
00003 //
00004 //                           Rythmos Package
00005 //                 Copyright (2006) Sandia Corporation
00006 //
00007 // Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
00008 // license for use of this work by or on behalf of the U.S. Government.
00009 //
00010 // This library is free software; you can redistribute it and/or modify
00011 // it under the terms of the GNU Lesser General Public License as
00012 // published by the Free Software Foundation; either version 2.1 of the
00013 // License, or (at your option) any later version.
00014 //
00015 // This library is distributed in the hope that it will be useful, but
00016 // WITHOUT ANY WARRANTY; without even the implied warranty of
00017 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
00018 // Lesser General Public License for more details.
00019 //
00020 // You should have received a copy of the GNU Lesser General Public
00021 // License along with this library; if not, write to the Free Software
00022 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307
00023 // USA
00024 // Questions? Contact Todd S. Coffey (tscoffe@sandia.gov)
00025 //
00026 // ***********************************************************************
00027 //@HEADER
00028 
00029 #ifndef RYTHMOS_STEPPER_HELPERS_DEF_HPP
00030 #define RYTHMOS_STEPPER_HELPERS_DEF_HPP
00031 
00032 #include "Rythmos_StepperHelpers_decl.hpp"
00033 #include "Rythmos_InterpolationBufferHelpers.hpp"
00034 #include "Rythmos_InterpolatorBaseHelpers.hpp"
00035 #include "Teuchos_Assert.hpp"
00036 #include "Thyra_AssertOp.hpp"
00037 #include "Thyra_VectorStdOps.hpp"
00038 
00039 
00040 namespace Rythmos {
00041 
00042 using Teuchos::ConstNonconstObjectContainer;
00043 
00044 template<class Scalar>
00045 void assertValidModel(
00046   const StepperBase<Scalar>& stepper,
00047   const Thyra::ModelEvaluator<Scalar>& model
00048   )
00049 {
00050 
00051   typedef Thyra::ModelEvaluatorBase MEB;
00052 
00053   TEUCHOS_ASSERT(stepper.acceptsModel());
00054 
00055   const MEB::InArgs<Scalar> inArgs = model.createInArgs();
00056   const MEB::OutArgs<Scalar> outArgs = model.createOutArgs();
00057 
00058   //TEUCHOS_ASSERT(inArgs.supports(MEB::IN_ARG_t));
00059   TEUCHOS_ASSERT(inArgs.supports(MEB::IN_ARG_x));
00060   TEUCHOS_ASSERT(outArgs.supports(MEB::OUT_ARG_f));
00061   
00062   if (stepper.isImplicit()) { // implicit stepper
00063     TEUCHOS_ASSERT( inArgs.supports(MEB::IN_ARG_x_dot) );
00064     TEUCHOS_ASSERT( inArgs.supports(MEB::IN_ARG_alpha) );
00065     TEUCHOS_ASSERT( inArgs.supports(MEB::IN_ARG_beta) );
00066     TEUCHOS_ASSERT( outArgs.supports(MEB::OUT_ARG_W) );
00067   } 
00068   //else { // explicit stepper
00069   //  TEUCHOS_ASSERT( !inArgs.supports(MEB::IN_ARG_x_dot) );
00070   //  TEUCHOS_ASSERT( !inArgs.supports(MEB::IN_ARG_alpha) );
00071   //  TEUCHOS_ASSERT( !inArgs.supports(MEB::IN_ARG_beta) );
00072   //  TEUCHOS_ASSERT( !outArgs.supports(MEB::OUT_ARG_W) );
00073   //}
00074 
00075 }
00076 
00077 
00078 template<class Scalar>
00079 bool setDefaultInitialConditionFromNominalValues(
00080   const Thyra::ModelEvaluator<Scalar>& model,
00081   const Ptr<StepperBase<Scalar> >& stepper
00082   )
00083 {
00084 
00085   typedef ScalarTraits<Scalar> ST;
00086   typedef Thyra::ModelEvaluatorBase MEB;
00087 
00088   if (isInitialized(*stepper))
00089     return false;  // Already has an initial condition
00090   
00091   MEB::InArgs<Scalar> initCond = model.getNominalValues();
00092 
00093   if (!is_null(initCond.get_x())) {
00094     // IC has x, we will assume that initCont.get_t() is the valid start time.
00095     // Therefore, we just need to check that x_dot is also set or we will
00096     // create a zero x_dot
00097 #ifdef HAVE_RYTHMOS_DEBUG
00098     THYRA_ASSERT_VEC_SPACES( "setInitialConditionIfExists(...)", 
00099       *model.get_x_space(), *initCond.get_x()->space() );
00100 #endif
00101     if (initCond.supports(MEB::IN_ARG_x_dot)) {
00102       if (is_null(initCond.get_x_dot())) {
00103         const RCP<Thyra::VectorBase<Scalar> > x_dot =
00104           createMember(model.get_x_space());
00105         assign(x_dot.ptr(), ST::zero());
00106       }
00107       else {
00108 #ifdef HAVE_RYTHMOS_DEBUG
00109         THYRA_ASSERT_VEC_SPACES( "setInitialConditionIfExists(...)", 
00110           *model.get_x_space(), *initCond.get_x_dot()->space() );
00111 #endif
00112       }
00113     }
00114     stepper->setInitialCondition(initCond);
00115     return true;
00116   }
00117 
00118   // The model has not nominal values for which to set the initial
00119   // conditions so wo don't do anything!  The stepper will still have not
00120   return false;
00121 
00122 }
00123 
00124 
00125 template<class Scalar>
00126 void restart( StepperBase<Scalar> *stepper )
00127 {
00128 #ifdef HAVE_RYTHMOS_DEBUG
00129   TEUCHOS_TEST_FOR_EXCEPT(0==stepper);
00130 #endif // HAVE_RYTHMOS_DEBUG
00131   typedef Thyra::ModelEvaluatorBase MEB;
00132   const Rythmos::StepStatus<double>
00133     stepStatus = stepper->getStepStatus();
00134   const RCP<const Thyra::ModelEvaluator<Scalar> >
00135     model = stepper->getModel();
00136   // First, copy all of the model's state, including parameter values etc.
00137   MEB::InArgs<double> initialCondition = model->createInArgs();
00138   initialCondition.setArgs(model->getNominalValues());
00139   // Set the current values of the state and time
00140   RCP<const Thyra::VectorBase<double> > x, x_dot;
00141   Rythmos::get_x_and_x_dot(*stepper,stepStatus.time,&x,&x_dot);
00142   initialCondition.set_x(x);
00143   initialCondition.set_x_dot(x_dot);
00144   initialCondition.set_t(stepStatus.time);
00145   // Set the new initial condition back on the stepper.  This will effectively
00146   // reset the stepper to think that it is starting over again (which it is).
00147   stepper->setInitialCondition(initialCondition);
00148 }
00149 
00150 template<class Scalar>
00151 void eval_model_explicit(
00152     const Thyra::ModelEvaluator<Scalar> &model,
00153     Thyra::ModelEvaluatorBase::InArgs<Scalar> &basePoint,
00154     const VectorBase<Scalar>& x_in,
00155     const typename Thyra::ModelEvaluatorBase::InArgs<Scalar>::ScalarMag &t_in,
00156     const Ptr<VectorBase<Scalar> >& f_out
00157     )
00158 {
00159   typedef Thyra::ModelEvaluatorBase MEB;
00160   MEB::InArgs<Scalar> inArgs = model.createInArgs();
00161   MEB::OutArgs<Scalar> outArgs = model.createOutArgs();
00162   inArgs.setArgs(basePoint);
00163   inArgs.set_x(Teuchos::rcp(&x_in,false));
00164   if (inArgs.supports(MEB::IN_ARG_t)) {
00165     inArgs.set_t(t_in);
00166   }
00167   // For model evaluators whose state function f(x, x_dot, t) describes
00168   // an implicit ODE, and which accept an optional x_dot input argument,
00169   // make sure the latter is set to null in order to request the evaluation
00170   // of a state function corresponding to the explicit ODE formulation
00171   // x_dot = f(x, t)
00172   if (inArgs.supports(MEB::IN_ARG_x_dot)) {
00173     inArgs.set_x_dot(Teuchos::null);
00174   }
00175   outArgs.set_f(Teuchos::rcp(&*f_out,false));
00176   model.evalModel(inArgs,outArgs);
00177 }
00178 
00179 
00180 #ifdef HAVE_THYRA_ME_POLYNOMIAL
00181 
00182 
00183 template<class Scalar>
00184 void eval_model_explicit_poly(
00185     const Thyra::ModelEvaluator<Scalar> &model,
00186     Thyra::ModelEvaluatorBase::InArgs<Scalar> &basePoint,
00187     const Teuchos::Polynomial< VectorBase<Scalar> > &x_poly,
00188     const typename Thyra::ModelEvaluatorBase::InArgs<Scalar>::ScalarMag &t,
00189     const Ptr<Teuchos::Polynomial<VectorBase<Scalar> > >& f_poly
00190     )
00191 {
00192   typedef Thyra::ModelEvaluatorBase MEB;
00193   MEB::InArgs<Scalar> inArgs = model.createInArgs();
00194   MEB::OutArgs<Scalar> outArgs = model.createOutArgs();
00195   inArgs.setArgs(basePoint);
00196   inArgs.set_x_poly(Teuchos::rcp(&x_poly,false));
00197   if (inArgs.supports(MEB::IN_ARG_t)) {
00198     inArgs.set_t(t);
00199   }
00200   outArgs.set_f_poly(Teuchos::rcp(&*f_poly,false));
00201 
00202   model.evalModel(inArgs,outArgs);
00203 }
00204 
00205 
00206 #endif // HAVE_THYRA_ME_POLYNOMIAL
00207 
00208 
00209 template<class Scalar>
00210 void defaultGetPoints(
00211     const Scalar& t_old, // required inArg
00212     const Ptr<const VectorBase<Scalar> >& x_old, // optional inArg
00213     const Ptr<const VectorBase<Scalar> >& xdot_old, // optional inArg
00214     const Scalar& t, // required inArg
00215     const Ptr<const VectorBase<Scalar> >& x, // optional inArg
00216     const Ptr<const VectorBase<Scalar> >& xdot, // optional inArg
00217     const Array<Scalar>& time_vec, // required inArg
00218     const Ptr<Array<Teuchos::RCP<const Thyra::VectorBase<Scalar> > > >& x_vec, // optional outArg
00219     const Ptr<Array<Teuchos::RCP<const Thyra::VectorBase<Scalar> > > >& xdot_vec, // optional outArg
00220     const Ptr<Array<typename Teuchos::ScalarTraits<Scalar>::magnitudeType> >& accuracy_vec, // optional outArg
00221     const Ptr<InterpolatorBase<Scalar> > interpolator // optional inArg (note:  not const)
00222     ) 
00223 {
00224   typedef Teuchos::ScalarTraits<Scalar> ST;
00225   assertTimePointsAreSorted(time_vec);
00226   TimeRange<Scalar> tr(t_old, t);
00227   TEUCHOS_ASSERT( tr.isValid() );
00228   if (!is_null(x_vec)) {
00229     x_vec->clear();
00230   }
00231   if (!is_null(xdot_vec)) {
00232     xdot_vec->clear();
00233   }
00234   if (!is_null(accuracy_vec)) {
00235     accuracy_vec->clear();
00236   }
00237   typename Array<Scalar>::const_iterator time_it = time_vec.begin();
00238   RCP<const VectorBase<Scalar> > tmpVec;
00239   RCP<const VectorBase<Scalar> > tmpVecDot;
00240   for (; time_it != time_vec.end() ; time_it++) {
00241     Scalar time = *time_it;
00242     asssertInTimeRange(tr, time);
00243     Scalar accuracy = ST::zero();
00244     if (compareTimeValues(time,t_old)==0) {
00245       if (!is_null(x_old)) {
00246         tmpVec = x_old->clone_v();
00247       }
00248       if (!is_null(xdot_old)) {
00249         tmpVecDot = xdot_old->clone_v();
00250       }
00251     } else if (compareTimeValues(time,t)==0) {
00252       if (!is_null(x)) {
00253         tmpVec = x->clone_v();
00254       }
00255       if (!is_null(xdot)) {
00256         tmpVecDot = xdot->clone_v();
00257       }
00258     } else {
00259       TEUCHOS_TEST_FOR_EXCEPTION(
00260           is_null(interpolator), std::logic_error,
00261           "Error, getPoints:  This stepper only supports time values on the boundaries!\n"
00262           );
00263       // At this point, we know time != t_old, time != t, interpolator != null, 
00264       // and time in [t_old,t], therefore, time in (t_old,t).  
00265       // t_old != t at this point because otherwise it would have been caught above.
00266       // Now use the interpolator to pass out the interior points
00267       typename DataStore<Scalar>::DataStoreVector_t ds_nodes;
00268       typename DataStore<Scalar>::DataStoreVector_t ds_out;
00269       {
00270         // t_old
00271         DataStore<Scalar> ds;
00272         ds.time = t_old;
00273         ds.x = rcp(x_old.get(),false);
00274         ds.xdot = rcp(xdot_old.get(),false);
00275         ds_nodes.push_back(ds);
00276       }
00277       {
00278         // t
00279         DataStore<Scalar> ds;
00280         ds.time = t;
00281         ds.x = rcp(x.get(),false);
00282         ds.xdot = rcp(xdot.get(),false);
00283         ds_nodes.push_back(ds);
00284       }
00285       Array<Scalar> time_vec_in;
00286       time_vec_in.push_back(time);
00287       interpolate<Scalar>(*interpolator,rcp(&ds_nodes,false),time_vec_in,&ds_out);
00288       Array<Scalar> time_vec_out;
00289       Array<RCP<const VectorBase<Scalar> > > x_vec_out;
00290       Array<RCP<const VectorBase<Scalar> > > xdot_vec_out;
00291       Array<typename Teuchos::ScalarTraits<Scalar>::magnitudeType> accuracy_vec_out;
00292       dataStoreVectorToVector(ds_out,&time_vec_out,&x_vec_out,&xdot_vec_out,&accuracy_vec_out);
00293       TEUCHOS_ASSERT( time_vec_out.length()==1 );
00294       tmpVec = x_vec_out[0];
00295       tmpVecDot = xdot_vec_out[0];
00296       accuracy = accuracy_vec_out[0];
00297     }
00298     if (!is_null(x_vec)) {
00299       x_vec->push_back(tmpVec);
00300     }
00301     if (!is_null(xdot_vec)) {
00302       xdot_vec->push_back(tmpVecDot);
00303     }
00304     if (!is_null(accuracy_vec)) {
00305       accuracy_vec->push_back(accuracy);
00306     }
00307     tmpVec = Teuchos::null;
00308     tmpVecDot = Teuchos::null;
00309   }
00310 }
00311 
00312 
00313 template<class Scalar>
00314   void setStepperModel(
00315       const Ptr<StepperBase<Scalar> >& stepper,
00316       const RCP<const Thyra::ModelEvaluator<Scalar> >& model
00317       )
00318 {
00319   stepper->setModel(model);
00320 }
00321 
00322 template<class Scalar>
00323   void setStepperModel(
00324       const Ptr<StepperBase<Scalar> >& stepper,
00325       const RCP<Thyra::ModelEvaluator<Scalar> >& model
00326       )
00327 {
00328   stepper->setNonconstModel(model);
00329 }
00330 
00331 template<class Scalar>
00332   void setStepperModel(
00333       const Ptr<StepperBase<Scalar> >& stepper,
00334       ConstNonconstObjectContainer<Thyra::ModelEvaluator<Scalar> >& model
00335       )
00336 {
00337   if (model.isConst()) {
00338     stepper->setModel(model.getConstObj());
00339   } 
00340   else {
00341     stepper->setNonconstModel(model.getNonconstObj());
00342   }
00343 }
00344 
00345 
00346 // 
00347 // Explicit Instantiation macro
00348 //
00349 // Must be expanded from within the Rythmos namespace!
00350 //
00351 
00352 
00353 #ifdef HAVE_THYRA_ME_POLYNOMIAL
00354 
00355 #define RYTHMOS_STEPPER_HELPERS_POLY_INSTANT(SCALAR) \
00356   template void eval_model_explicit_poly( \
00357       const Thyra::ModelEvaluator< SCALAR > &model, \
00358       Thyra::ModelEvaluatorBase::InArgs< SCALAR > &basePoint, \
00359       const Teuchos::Polynomial< VectorBase< SCALAR > > &x_poly, \
00360       const Thyra::ModelEvaluatorBase::InArgs< SCALAR >::ScalarMag &t, \
00361       const Ptr<Teuchos::Polynomial<VectorBase< SCALAR > > >& f_poly \
00362       );
00363 
00364 #else // HAVE_THYRA_ME_POLYNOMIAL
00365 
00366 #define RYTHMOS_STEPPER_HELPERS_POLY_INSTANT(SCALAR)
00367 
00368 #endif // HAVE_THYRA_ME_POLYNOMIAL
00369 
00370 
00371 #define RYTHMOS_STEPPER_HELPERS_INSTANT(SCALAR) \
00372   \
00373   template void assertValidModel( \
00374     const StepperBase< SCALAR >& stepper, \
00375     const Thyra::ModelEvaluator< SCALAR >& model \
00376     ); \
00377   template bool setDefaultInitialConditionFromNominalValues( \
00378     const Thyra::ModelEvaluator< SCALAR >& model, \
00379     const Ptr<StepperBase< SCALAR > >& stepper \
00380     ); \
00381   template void restart( StepperBase< SCALAR > *stepper ); \
00382   \
00383   template void eval_model_explicit( \
00384       const Thyra::ModelEvaluator< SCALAR > &model, \
00385       Thyra::ModelEvaluatorBase::InArgs< SCALAR > &basePoint, \
00386       const VectorBase< SCALAR >& x_in, \
00387       const Thyra::ModelEvaluatorBase::InArgs< SCALAR >::ScalarMag &t_in, \
00388       const Ptr<VectorBase< SCALAR > >& f_out \
00389       ); \
00390   \
00391   RYTHMOS_STEPPER_HELPERS_POLY_INSTANT(SCALAR) \
00392   \
00393   template void defaultGetPoints( \
00394       const  SCALAR & t_old, \
00395       const Ptr<const VectorBase< SCALAR > >& x_old, \
00396       const Ptr<const VectorBase< SCALAR > >& xdot_old, \
00397       const  SCALAR & t, \
00398       const Ptr<const VectorBase< SCALAR > >& x, \
00399       const Ptr<const VectorBase< SCALAR > >& xdot, \
00400       const Array< SCALAR >& time_vec, \
00401       const Ptr<Array<Teuchos::RCP<const Thyra::VectorBase< SCALAR > > > >& x_vec, \
00402       const Ptr<Array<Teuchos::RCP<const Thyra::VectorBase< SCALAR > > > >& xdot_vec, \
00403       const Ptr<Array<Teuchos::ScalarTraits< SCALAR >::magnitudeType> >& accuracy_vec, \
00404       const Ptr<InterpolatorBase< SCALAR > > interpolator  \
00405       );  \
00406   \
00407   template void setStepperModel( \
00408         const Ptr<StepperBase< SCALAR > >& stepper, \
00409         const RCP<const Thyra::ModelEvaluator< SCALAR > >& model \
00410         ); \
00411   \
00412   template void setStepperModel( \
00413         const Ptr<StepperBase< SCALAR > >& stepper, \
00414         const RCP<Thyra::ModelEvaluator< SCALAR > >& model \
00415         ); \
00416   \
00417   template void setStepperModel( \
00418         const Ptr<StepperBase< SCALAR > >& stepper, \
00419         Teuchos::ConstNonconstObjectContainer<Thyra::ModelEvaluator< SCALAR > >& model \
00420         );
00421 
00422 } // namespace Rythmos
00423 
00424 
00425 #endif // RYTHMOS_STEPPER_HELPERS_DEF_HPP
 All Classes Functions Variables Typedefs Friends