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 
00038 
00039 namespace Rythmos {
00040 
00041 
00042 template<class Scalar>
00043 void assertValidModel(
00044   const StepperBase<Scalar>& stepper,
00045   const Thyra::ModelEvaluator<Scalar>& model
00046   )
00047 {
00048 
00049   typedef Thyra::ModelEvaluatorBase MEB;
00050 
00051   TEUCHOS_ASSERT(stepper.acceptsModel());
00052 
00053   const MEB::InArgs<Scalar> inArgs = model.createInArgs();
00054   const MEB::OutArgs<Scalar> outArgs = model.createOutArgs();
00055 
00056   //TEUCHOS_ASSERT(inArgs.supports(MEB::IN_ARG_t));
00057   TEUCHOS_ASSERT(inArgs.supports(MEB::IN_ARG_x));
00058   TEUCHOS_ASSERT(outArgs.supports(MEB::OUT_ARG_f));
00059   
00060   if (stepper.isImplicit()) { // implicit stepper
00061     TEUCHOS_ASSERT( inArgs.supports(MEB::IN_ARG_x_dot) );
00062     TEUCHOS_ASSERT( inArgs.supports(MEB::IN_ARG_alpha) );
00063     TEUCHOS_ASSERT( inArgs.supports(MEB::IN_ARG_beta) );
00064     TEUCHOS_ASSERT( outArgs.supports(MEB::OUT_ARG_W) );
00065   } 
00066   //else { // explicit stepper
00067   //  TEUCHOS_ASSERT( !inArgs.supports(MEB::IN_ARG_x_dot) );
00068   //  TEUCHOS_ASSERT( !inArgs.supports(MEB::IN_ARG_alpha) );
00069   //  TEUCHOS_ASSERT( !inArgs.supports(MEB::IN_ARG_beta) );
00070   //  TEUCHOS_ASSERT( !outArgs.supports(MEB::OUT_ARG_W) );
00071   //}
00072 
00073 }
00074 
00075 
00076 template<class Scalar>
00077 bool setDefaultInitialConditionFromNominalValues(
00078   const Thyra::ModelEvaluator<Scalar>& model,
00079   const Ptr<StepperBase<Scalar> >& stepper
00080   )
00081 {
00082 
00083   typedef ScalarTraits<Scalar> ST;
00084   typedef Thyra::ModelEvaluatorBase MEB;
00085 
00086   if (isInitialized(*stepper))
00087     return false;  // Already has an initial condition
00088   
00089   MEB::InArgs<Scalar> initCond = model.getNominalValues();
00090 
00091   if (!is_null(initCond.get_x())) {
00092     // IC has x, we will assume that initCont.get_t() is the valid start time.
00093     // Therefore, we just need to check that x_dot is also set or we will
00094     // create a zero x_dot
00095 #ifdef RYTHMOS_DEBUG
00096     THYRA_ASSERT_VEC_SPACES( "setInitialConditionIfExists(...)", 
00097       *model.get_x_space(), *initCond.get_x()->space() );
00098 #endif
00099     if (initCond.supports(MEB::IN_ARG_x_dot)) {
00100       if (is_null(initCond.get_x_dot())) {
00101         const RCP<Thyra::VectorBase<Scalar> > x_dot =
00102           createMember(model.get_x_space());
00103         assign(x_dot.ptr(), ST::zero());
00104       }
00105       else {
00106 #ifdef RYTHMOS_DEBUG
00107         THYRA_ASSERT_VEC_SPACES( "setInitialConditionIfExists(...)", 
00108           *model.get_x_space(), *initCond.get_x_dot()->space() );
00109 #endif
00110       }
00111     }
00112     stepper->setInitialCondition(initCond);
00113     return true;
00114   }
00115 
00116   // The model has not nominal values for which to set the initial
00117   // conditions so wo don't do anything!  The stepper will still have not
00118   return false;
00119 
00120 }
00121 
00122 
00123 template<class Scalar>
00124 void restart( StepperBase<Scalar> *stepper )
00125 {
00126 #ifdef RYTHMOS_DEBUG
00127   TEST_FOR_EXCEPT(0==stepper);
00128 #endif // RYTHMOS_DEBUG
00129   typedef Thyra::ModelEvaluatorBase MEB;
00130   const Rythmos::StepStatus<double>
00131     stepStatus = stepper->getStepStatus();
00132   const RCP<const Thyra::ModelEvaluator<Scalar> >
00133     model = stepper->getModel();
00134   // First, copy all of the model's state, including parameter values etc.
00135   MEB::InArgs<double> initialCondition = model->createInArgs();
00136   initialCondition.setArgs(model->getNominalValues());
00137   // Set the current values of the state and time
00138   RCP<const Thyra::VectorBase<double> > x, x_dot;
00139   Rythmos::get_x_and_x_dot(*stepper,stepStatus.time,&x,&x_dot);
00140   initialCondition.set_x(x);
00141   initialCondition.set_x_dot(x_dot);
00142   initialCondition.set_t(stepStatus.time);
00143   // Set the new initial condition back on the stepper.  This will effectively
00144   // reset the stepper to think that it is starting over again (which it is).
00145   stepper->setInitialCondition(initialCondition);
00146 }
00147 
00148 template<class Scalar>
00149 void eval_model_explicit(
00150     const Thyra::ModelEvaluator<Scalar> &model,
00151     Thyra::ModelEvaluatorBase::InArgs<Scalar> &basePoint,
00152     const VectorBase<Scalar>& x_in,
00153     const typename Thyra::ModelEvaluatorBase::InArgs<Scalar>::ScalarMag &t_in,
00154     const Ptr<VectorBase<Scalar> >& f_out
00155     )
00156 {
00157   typedef Thyra::ModelEvaluatorBase MEB;
00158   MEB::InArgs<Scalar> inArgs = model.createInArgs();
00159   MEB::OutArgs<Scalar> outArgs = model.createOutArgs();
00160   inArgs.setArgs(basePoint);
00161   inArgs.set_x(Teuchos::rcp(&x_in,false));
00162   if (inArgs.supports(MEB::IN_ARG_t)) {
00163     inArgs.set_t(t_in);
00164   }
00165   outArgs.set_f(Teuchos::rcp(&*f_out,false));
00166   model.evalModel(inArgs,outArgs);
00167 }
00168 
00169 template<class Scalar>
00170 void eval_model_explicit_poly(
00171     const Thyra::ModelEvaluator<Scalar> &model,
00172     Thyra::ModelEvaluatorBase::InArgs<Scalar> &basePoint,
00173     const Teuchos::Polynomial< VectorBase<Scalar> > &x_poly,
00174     const typename Thyra::ModelEvaluatorBase::InArgs<Scalar>::ScalarMag &t,
00175     const Ptr<Teuchos::Polynomial<VectorBase<Scalar> > >& f_poly
00176     )
00177 {
00178   typedef Thyra::ModelEvaluatorBase MEB;
00179   MEB::InArgs<Scalar> inArgs = model.createInArgs();
00180   MEB::OutArgs<Scalar> outArgs = model.createOutArgs();
00181   inArgs.setArgs(basePoint);
00182   inArgs.set_x_poly(Teuchos::rcp(&x_poly,false));
00183   if (inArgs.supports(MEB::IN_ARG_t)) {
00184     inArgs.set_t(t);
00185   }
00186   outArgs.set_f_poly(Teuchos::rcp(&*f_poly,false));
00187 
00188   model.evalModel(inArgs,outArgs);
00189 }
00190 
00191 
00192 template<class Scalar>
00193 void defaultGetPoints(
00194     const Scalar& t_old, // required inArg
00195     const Ptr<const VectorBase<Scalar> >& x_old, // optional inArg
00196     const Ptr<const VectorBase<Scalar> >& xdot_old, // optional inArg
00197     const Scalar& t, // required inArg
00198     const Ptr<const VectorBase<Scalar> >& x, // optional inArg
00199     const Ptr<const VectorBase<Scalar> >& xdot, // optional inArg
00200     const Array<Scalar>& time_vec, // required inArg
00201     const Ptr<Array<Teuchos::RCP<const Thyra::VectorBase<Scalar> > > >& x_vec, // optional outArg
00202     const Ptr<Array<Teuchos::RCP<const Thyra::VectorBase<Scalar> > > >& xdot_vec, // optional outArg
00203     const Ptr<Array<typename Teuchos::ScalarTraits<Scalar>::magnitudeType> >& accuracy_vec, // optional outArg
00204     const Ptr<InterpolatorBase<Scalar> > interpolator // optional inArg (note:  not const)
00205     ) 
00206 {
00207   typedef Teuchos::ScalarTraits<Scalar> ST;
00208   assertTimePointsAreSorted(time_vec);
00209   TimeRange<Scalar> tr(t_old, t);
00210   TEUCHOS_ASSERT( tr.isValid() );
00211   if (!is_null(x_vec)) {
00212     x_vec->clear();
00213   }
00214   if (!is_null(xdot_vec)) {
00215     xdot_vec->clear();
00216   }
00217   if (!is_null(accuracy_vec)) {
00218     accuracy_vec->clear();
00219   }
00220   typename Array<Scalar>::const_iterator time_it = time_vec.begin();
00221   RCP<const VectorBase<Scalar> > tmpVec;
00222   RCP<const VectorBase<Scalar> > tmpVecDot;
00223   for (; time_it != time_vec.end() ; time_it++) {
00224     Scalar time = *time_it;
00225     asssertInTimeRange(tr, time);
00226     Scalar accuracy = ST::zero();
00227     if (compareTimeValues(time,t_old)==0) {
00228       if (!is_null(x_old)) {
00229         tmpVec = x_old->clone_v();
00230       }
00231       if (!is_null(xdot_old)) {
00232         tmpVecDot = xdot_old->clone_v();
00233       }
00234     } else if (compareTimeValues(time,t)==0) {
00235       if (!is_null(x)) {
00236         tmpVec = x->clone_v();
00237       }
00238       if (!is_null(xdot)) {
00239         tmpVecDot = xdot->clone_v();
00240       }
00241     } else {
00242       TEST_FOR_EXCEPTION(
00243           is_null(interpolator), std::logic_error,
00244           "Error, getPoints:  This stepper only supports time values on the boundaries!\n"
00245           );
00246       // At this point, we know time != t_old, time != t, interpolator != null, 
00247       // and time in [t_old,t], therefore, time in (t_old,t).  
00248       // t_old != t at this point because otherwise it would have been caught above.
00249       // Now use the interpolator to pass out the interior points
00250       typename DataStore<Scalar>::DataStoreVector_t ds_nodes;
00251       typename DataStore<Scalar>::DataStoreVector_t ds_out;
00252       {
00253         // t_old
00254         DataStore<Scalar> ds;
00255         ds.time = t_old;
00256         ds.x = rcp(x_old.get(),false);
00257         ds.xdot = rcp(xdot_old.get(),false);
00258         ds_nodes.push_back(ds);
00259       }
00260       {
00261         // t
00262         DataStore<Scalar> ds;
00263         ds.time = t;
00264         ds.x = rcp(x.get(),false);
00265         ds.xdot = rcp(xdot.get(),false);
00266         ds_nodes.push_back(ds);
00267       }
00268       Array<Scalar> time_vec_in;
00269       time_vec_in.push_back(time);
00270       interpolate<Scalar>(*interpolator,rcp(&ds_nodes,false),time_vec_in,&ds_out);
00271       Array<Scalar> time_vec_out;
00272       Array<RCP<const VectorBase<Scalar> > > x_vec_out;
00273       Array<RCP<const VectorBase<Scalar> > > xdot_vec_out;
00274       Array<typename Teuchos::ScalarTraits<Scalar>::magnitudeType> accuracy_vec_out;
00275       dataStoreVectorToVector(ds_out,&time_vec_out,&x_vec_out,&xdot_vec_out,&accuracy_vec_out);
00276       TEUCHOS_ASSERT( time_vec_out.length()==1 );
00277       tmpVec = x_vec_out[0];
00278       tmpVecDot = xdot_vec_out[0];
00279       accuracy = accuracy_vec_out[0];
00280     }
00281     if (!is_null(x_vec)) {
00282       x_vec->push_back(tmpVec);
00283     }
00284     if (!is_null(xdot_vec)) {
00285       xdot_vec->push_back(tmpVecDot);
00286     }
00287     if (!is_null(accuracy_vec)) {
00288       accuracy_vec->push_back(accuracy);
00289     }
00290     tmpVec = Teuchos::null;
00291     tmpVecDot = Teuchos::null;
00292   }
00293 }
00294 
00295 // 
00296 // Explicit Instantiation macro
00297 //
00298 // Must be expanded from within the Rythmos namespace!
00299 //
00300 
00301 #define RYTHMOS_STEPPER_HELPERS_INSTANT(SCALAR) \
00302   \
00303   template void assertValidModel( \
00304     const StepperBase< SCALAR >& stepper, \
00305     const Thyra::ModelEvaluator< SCALAR >& model \
00306     ); \
00307   template bool setDefaultInitialConditionFromNominalValues( \
00308     const Thyra::ModelEvaluator< SCALAR >& model, \
00309     const Ptr<StepperBase< SCALAR > >& stepper \
00310     ); \
00311   template void restart( StepperBase< SCALAR > *stepper ); \
00312   \
00313   template void eval_model_explicit( \
00314       const Thyra::ModelEvaluator< SCALAR > &model, \
00315       Thyra::ModelEvaluatorBase::InArgs< SCALAR > &basePoint, \
00316       const VectorBase< SCALAR >& x_in, \
00317       const Thyra::ModelEvaluatorBase::InArgs< SCALAR >::ScalarMag &t_in, \
00318       const Ptr<VectorBase< SCALAR > >& f_out \
00319       ); \
00320   \
00321   template void eval_model_explicit_poly( \
00322       const Thyra::ModelEvaluator< SCALAR > &model, \
00323       Thyra::ModelEvaluatorBase::InArgs< SCALAR > &basePoint, \
00324       const Teuchos::Polynomial< VectorBase< SCALAR > > &x_poly, \
00325       const Thyra::ModelEvaluatorBase::InArgs< SCALAR >::ScalarMag &t, \
00326       const Ptr<Teuchos::Polynomial<VectorBase< SCALAR > > >& f_poly \
00327       ); \
00328   \
00329   template void defaultGetPoints( \
00330       const  SCALAR & t_old, \
00331       const Ptr<const VectorBase< SCALAR > >& x_old, \
00332       const Ptr<const VectorBase< SCALAR > >& xdot_old, \
00333       const  SCALAR & t, \
00334       const Ptr<const VectorBase< SCALAR > >& x, \
00335       const Ptr<const VectorBase< SCALAR > >& xdot, \
00336       const Array< SCALAR >& time_vec, \
00337       const Ptr<Array<Teuchos::RCP<const Thyra::VectorBase< SCALAR > > > >& x_vec, \
00338       const Ptr<Array<Teuchos::RCP<const Thyra::VectorBase< SCALAR > > > >& xdot_vec, \
00339       const Ptr<Array<Teuchos::ScalarTraits< SCALAR >::magnitudeType> >& accuracy_vec, \
00340       const Ptr<InterpolatorBase< SCALAR > > interpolator  \
00341       ); 
00342 
00343 } // namespace Rythmos
00344 
00345 
00346 #endif // RYTHMOS_STEPPER_HELPERS_DEF_HPP

Generated on Wed May 12 21:25:44 2010 for Rythmos - Transient Integration for Differential Equations by  doxygen 1.4.7