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 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 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 RYTHMOS_DEBUG
00129   TEST_FOR_EXCEPT(0==stepper);
00130 #endif // 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   outArgs.set_f(Teuchos::rcp(&*f_out,false));
00168   model.evalModel(inArgs,outArgs);
00169 }
00170 
00171 
00172 #ifdef HAVE_THYRA_ME_POLYNOMIAL
00173 
00174 
00175 template<class Scalar>
00176 void eval_model_explicit_poly(
00177     const Thyra::ModelEvaluator<Scalar> &model,
00178     Thyra::ModelEvaluatorBase::InArgs<Scalar> &basePoint,
00179     const Teuchos::Polynomial< VectorBase<Scalar> > &x_poly,
00180     const typename Thyra::ModelEvaluatorBase::InArgs<Scalar>::ScalarMag &t,
00181     const Ptr<Teuchos::Polynomial<VectorBase<Scalar> > >& f_poly
00182     )
00183 {
00184   typedef Thyra::ModelEvaluatorBase MEB;
00185   MEB::InArgs<Scalar> inArgs = model.createInArgs();
00186   MEB::OutArgs<Scalar> outArgs = model.createOutArgs();
00187   inArgs.setArgs(basePoint);
00188   inArgs.set_x_poly(Teuchos::rcp(&x_poly,false));
00189   if (inArgs.supports(MEB::IN_ARG_t)) {
00190     inArgs.set_t(t);
00191   }
00192   outArgs.set_f_poly(Teuchos::rcp(&*f_poly,false));
00193 
00194   model.evalModel(inArgs,outArgs);
00195 }
00196 
00197 
00198 #endif // HAVE_THYRA_ME_POLYNOMIAL
00199 
00200 
00201 template<class Scalar>
00202 void defaultGetPoints(
00203     const Scalar& t_old, // required inArg
00204     const Ptr<const VectorBase<Scalar> >& x_old, // optional inArg
00205     const Ptr<const VectorBase<Scalar> >& xdot_old, // optional inArg
00206     const Scalar& t, // required inArg
00207     const Ptr<const VectorBase<Scalar> >& x, // optional inArg
00208     const Ptr<const VectorBase<Scalar> >& xdot, // optional inArg
00209     const Array<Scalar>& time_vec, // required inArg
00210     const Ptr<Array<Teuchos::RCP<const Thyra::VectorBase<Scalar> > > >& x_vec, // optional outArg
00211     const Ptr<Array<Teuchos::RCP<const Thyra::VectorBase<Scalar> > > >& xdot_vec, // optional outArg
00212     const Ptr<Array<typename Teuchos::ScalarTraits<Scalar>::magnitudeType> >& accuracy_vec, // optional outArg
00213     const Ptr<InterpolatorBase<Scalar> > interpolator // optional inArg (note:  not const)
00214     ) 
00215 {
00216   typedef Teuchos::ScalarTraits<Scalar> ST;
00217   assertTimePointsAreSorted(time_vec);
00218   TimeRange<Scalar> tr(t_old, t);
00219   TEUCHOS_ASSERT( tr.isValid() );
00220   if (!is_null(x_vec)) {
00221     x_vec->clear();
00222   }
00223   if (!is_null(xdot_vec)) {
00224     xdot_vec->clear();
00225   }
00226   if (!is_null(accuracy_vec)) {
00227     accuracy_vec->clear();
00228   }
00229   typename Array<Scalar>::const_iterator time_it = time_vec.begin();
00230   RCP<const VectorBase<Scalar> > tmpVec;
00231   RCP<const VectorBase<Scalar> > tmpVecDot;
00232   for (; time_it != time_vec.end() ; time_it++) {
00233     Scalar time = *time_it;
00234     asssertInTimeRange(tr, time);
00235     Scalar accuracy = ST::zero();
00236     if (compareTimeValues(time,t_old)==0) {
00237       if (!is_null(x_old)) {
00238         tmpVec = x_old->clone_v();
00239       }
00240       if (!is_null(xdot_old)) {
00241         tmpVecDot = xdot_old->clone_v();
00242       }
00243     } else if (compareTimeValues(time,t)==0) {
00244       if (!is_null(x)) {
00245         tmpVec = x->clone_v();
00246       }
00247       if (!is_null(xdot)) {
00248         tmpVecDot = xdot->clone_v();
00249       }
00250     } else {
00251       TEST_FOR_EXCEPTION(
00252           is_null(interpolator), std::logic_error,
00253           "Error, getPoints:  This stepper only supports time values on the boundaries!\n"
00254           );
00255       // At this point, we know time != t_old, time != t, interpolator != null, 
00256       // and time in [t_old,t], therefore, time in (t_old,t).  
00257       // t_old != t at this point because otherwise it would have been caught above.
00258       // Now use the interpolator to pass out the interior points
00259       typename DataStore<Scalar>::DataStoreVector_t ds_nodes;
00260       typename DataStore<Scalar>::DataStoreVector_t ds_out;
00261       {
00262         // t_old
00263         DataStore<Scalar> ds;
00264         ds.time = t_old;
00265         ds.x = rcp(x_old.get(),false);
00266         ds.xdot = rcp(xdot_old.get(),false);
00267         ds_nodes.push_back(ds);
00268       }
00269       {
00270         // t
00271         DataStore<Scalar> ds;
00272         ds.time = t;
00273         ds.x = rcp(x.get(),false);
00274         ds.xdot = rcp(xdot.get(),false);
00275         ds_nodes.push_back(ds);
00276       }
00277       Array<Scalar> time_vec_in;
00278       time_vec_in.push_back(time);
00279       interpolate<Scalar>(*interpolator,rcp(&ds_nodes,false),time_vec_in,&ds_out);
00280       Array<Scalar> time_vec_out;
00281       Array<RCP<const VectorBase<Scalar> > > x_vec_out;
00282       Array<RCP<const VectorBase<Scalar> > > xdot_vec_out;
00283       Array<typename Teuchos::ScalarTraits<Scalar>::magnitudeType> accuracy_vec_out;
00284       dataStoreVectorToVector(ds_out,&time_vec_out,&x_vec_out,&xdot_vec_out,&accuracy_vec_out);
00285       TEUCHOS_ASSERT( time_vec_out.length()==1 );
00286       tmpVec = x_vec_out[0];
00287       tmpVecDot = xdot_vec_out[0];
00288       accuracy = accuracy_vec_out[0];
00289     }
00290     if (!is_null(x_vec)) {
00291       x_vec->push_back(tmpVec);
00292     }
00293     if (!is_null(xdot_vec)) {
00294       xdot_vec->push_back(tmpVecDot);
00295     }
00296     if (!is_null(accuracy_vec)) {
00297       accuracy_vec->push_back(accuracy);
00298     }
00299     tmpVec = Teuchos::null;
00300     tmpVecDot = Teuchos::null;
00301   }
00302 }
00303 
00304 
00305 template<class Scalar>
00306   void setStepperModel(
00307       const Ptr<StepperBase<Scalar> >& stepper,
00308       const RCP<const Thyra::ModelEvaluator<Scalar> >& model
00309       )
00310 {
00311   stepper->setModel(model);
00312 }
00313 
00314 template<class Scalar>
00315   void setStepperModel(
00316       const Ptr<StepperBase<Scalar> >& stepper,
00317       const RCP<Thyra::ModelEvaluator<Scalar> >& model
00318       )
00319 {
00320   stepper->setNonconstModel(model);
00321 }
00322 
00323 template<class Scalar>
00324   void setStepperModel(
00325       const Ptr<StepperBase<Scalar> >& stepper,
00326       ConstNonconstObjectContainer<Thyra::ModelEvaluator<Scalar> >& model
00327       )
00328 {
00329   if (model.isConst()) {
00330     stepper->setModel(model.getConstObj());
00331   } 
00332   else {
00333     stepper->setNonconstModel(model.getNonconstObj());
00334   }
00335 }
00336 
00337 
00338 // 
00339 // Explicit Instantiation macro
00340 //
00341 // Must be expanded from within the Rythmos namespace!
00342 //
00343 
00344 
00345 #ifdef HAVE_THYRA_ME_POLYNOMIAL
00346 
00347 #define RYTHMOS_STEPPER_HELPERS_POLY_INSTANT(SCALAR) \
00348   template void eval_model_explicit_poly( \
00349       const Thyra::ModelEvaluator< SCALAR > &model, \
00350       Thyra::ModelEvaluatorBase::InArgs< SCALAR > &basePoint, \
00351       const Teuchos::Polynomial< VectorBase< SCALAR > > &x_poly, \
00352       const Thyra::ModelEvaluatorBase::InArgs< SCALAR >::ScalarMag &t, \
00353       const Ptr<Teuchos::Polynomial<VectorBase< SCALAR > > >& f_poly \
00354       );
00355 
00356 #else // HAVE_THYRA_ME_POLYNOMIAL
00357 
00358 #define RYTHMOS_STEPPER_HELPERS_POLY_INSTANT(SCALAR)
00359 
00360 #endif // HAVE_THYRA_ME_POLYNOMIAL
00361 
00362 
00363 #define RYTHMOS_STEPPER_HELPERS_INSTANT(SCALAR) \
00364   \
00365   template void assertValidModel( \
00366     const StepperBase< SCALAR >& stepper, \
00367     const Thyra::ModelEvaluator< SCALAR >& model \
00368     ); \
00369   template bool setDefaultInitialConditionFromNominalValues( \
00370     const Thyra::ModelEvaluator< SCALAR >& model, \
00371     const Ptr<StepperBase< SCALAR > >& stepper \
00372     ); \
00373   template void restart( StepperBase< SCALAR > *stepper ); \
00374   \
00375   template void eval_model_explicit( \
00376       const Thyra::ModelEvaluator< SCALAR > &model, \
00377       Thyra::ModelEvaluatorBase::InArgs< SCALAR > &basePoint, \
00378       const VectorBase< SCALAR >& x_in, \
00379       const Thyra::ModelEvaluatorBase::InArgs< SCALAR >::ScalarMag &t_in, \
00380       const Ptr<VectorBase< SCALAR > >& f_out \
00381       ); \
00382   \
00383   RYTHMOS_STEPPER_HELPERS_POLY_INSTANT(SCALAR) \
00384   \
00385   template void defaultGetPoints( \
00386       const  SCALAR & t_old, \
00387       const Ptr<const VectorBase< SCALAR > >& x_old, \
00388       const Ptr<const VectorBase< SCALAR > >& xdot_old, \
00389       const  SCALAR & t, \
00390       const Ptr<const VectorBase< SCALAR > >& x, \
00391       const Ptr<const VectorBase< SCALAR > >& xdot, \
00392       const Array< SCALAR >& time_vec, \
00393       const Ptr<Array<Teuchos::RCP<const Thyra::VectorBase< SCALAR > > > >& x_vec, \
00394       const Ptr<Array<Teuchos::RCP<const Thyra::VectorBase< SCALAR > > > >& xdot_vec, \
00395       const Ptr<Array<Teuchos::ScalarTraits< SCALAR >::magnitudeType> >& accuracy_vec, \
00396       const Ptr<InterpolatorBase< SCALAR > > interpolator  \
00397       );  \
00398   \
00399   template void setStepperModel( \
00400         const Ptr<StepperBase< SCALAR > >& stepper, \
00401         const RCP<const Thyra::ModelEvaluator< SCALAR > >& model \
00402         ); \
00403   \
00404   template void setStepperModel( \
00405         const Ptr<StepperBase< SCALAR > >& stepper, \
00406         const RCP<Thyra::ModelEvaluator< SCALAR > >& model \
00407         ); \
00408   \
00409   template void setStepperModel( \
00410         const Ptr<StepperBase< SCALAR > >& stepper, \
00411         Teuchos::ConstNonconstObjectContainer<Thyra::ModelEvaluator< SCALAR > >& model \
00412         );
00413 
00414 } // namespace Rythmos
00415 
00416 
00417 #endif // RYTHMOS_STEPPER_HELPERS_DEF_HPP

Generated on Tue Jul 13 09:23:54 2010 for Rythmos - Transient Integration for Differential Equations by  doxygen 1.4.7