Rythmos_ImplicitRKStepper_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_IMPLICIT_RK_STEPPER_DEF_H
00030 #define Rythmos_IMPLICIT_RK_STEPPER_DEF_H
00031 
00032 #include "Rythmos_ImplicitRKStepper_decl.hpp"
00033 
00034 #include "Rythmos_StepperHelpers.hpp"
00035 #include "Rythmos_SingleResidualModelEvaluator.hpp"
00036 #include "Rythmos_ImplicitRKModelEvaluator.hpp"
00037 #include "Rythmos_DiagonalImplicitRKModelEvaluator.hpp"
00038 #include "Rythmos_RKButcherTableau.hpp"
00039 #include "Rythmos_RKButcherTableauHelpers.hpp"
00040 
00041 #include "Thyra_ModelEvaluatorHelpers.hpp"
00042 #include "Thyra_ProductVectorSpaceBase.hpp"
00043 #include "Thyra_AssertOp.hpp"
00044 #include "Thyra_TestingTools.hpp"
00045 
00046 #include "Teuchos_VerboseObjectParameterListHelpers.hpp"
00047 #include "Teuchos_as.hpp"
00048 
00049 
00050 namespace Rythmos {
00051 
00056 template<class Scalar>
00057 RCP<ImplicitRKStepper<Scalar> >
00058 implicitRKStepper()
00059 {
00060   RCP<ImplicitRKStepper<Scalar> > stepper(new ImplicitRKStepper<Scalar>());
00061   return stepper;
00062 }
00063 
00064 template<class Scalar>
00065 RCP<ImplicitRKStepper<Scalar> >
00066 implicitRKStepper(
00067   const RCP<const Thyra::ModelEvaluator<Scalar> >  &model,
00068   const RCP<Thyra::NonlinearSolverBase<Scalar> >  &solver,
00069   const RCP<Thyra::LinearOpWithSolveFactoryBase<Scalar> > &irk_W_factory,
00070   const RCP<const RKButcherTableauBase<Scalar> > &irkbt
00071   )
00072 {
00073   RCP<ImplicitRKStepper<Scalar> > stepper(new ImplicitRKStepper<Scalar>());
00074 
00075   stepper->setModel(model);
00076   stepper->setSolver(solver);
00077   stepper->set_W_factory(irk_W_factory);
00078   stepper->setRKButcherTableau(irkbt);
00079   return stepper;
00080 }
00081 
00082 
00083 // ////////////////////////////
00084 // Defintions
00085 
00086 
00087 // Constructors, intializers, Misc.
00088 
00089 
00090 template<class Scalar>
00091 ImplicitRKStepper<Scalar>::ImplicitRKStepper()
00092 {
00093   this->defaultInitializeAll_();
00094   irkButcherTableau_ = rKButcherTableau<Scalar>();
00095   numSteps_ = 0;
00096 }
00097 
00098 template<class Scalar>
00099 void ImplicitRKStepper<Scalar>::defaultInitializeAll_()
00100 {
00101   isInitialized_ = false;
00102   model_ = Teuchos::null;
00103   solver_ = Teuchos::null;
00104   irk_W_factory_ = Teuchos::null;
00105   paramList_ = Teuchos::null;
00106   //basePoint_;
00107   x_ = Teuchos::null;
00108   x_old_ = Teuchos::null;
00109   x_dot_ = Teuchos::null;
00110   //timeRange_;
00111   irkModel_ = Teuchos::null;
00112   irkButcherTableau_ = Teuchos::null;
00113   isDirk_ = false;
00114   numSteps_ = -1;
00115   haveInitialCondition_ = false;
00116   x_stage_bar_ = Teuchos::null;
00117 }
00118 
00119 template<class Scalar>
00120 void ImplicitRKStepper<Scalar>::set_W_factory(
00121   const RCP<Thyra::LinearOpWithSolveFactoryBase<Scalar> > &irk_W_factory
00122   )
00123 {
00124   TEUCHOS_ASSERT( !is_null(irk_W_factory) ); 
00125   irk_W_factory_ = irk_W_factory;
00126 }
00127 
00128 template<class Scalar>
00129 RCP<const Thyra::LinearOpWithSolveFactoryBase<Scalar> > ImplicitRKStepper<Scalar>::get_W_factory() const
00130 {
00131   return irk_W_factory_;
00132 }
00133 
00134 // Overridden from SolverAcceptingStepperBase
00135 
00136 
00137 template<class Scalar>
00138 void ImplicitRKStepper<Scalar>::setSolver(
00139   const RCP<Thyra::NonlinearSolverBase<Scalar> > &solver
00140   )
00141 {
00142   TEST_FOR_EXCEPT(is_null(solver));
00143   solver_ = solver;
00144 }
00145 
00146 
00147 template<class Scalar>
00148 RCP<Thyra::NonlinearSolverBase<Scalar> >
00149 ImplicitRKStepper<Scalar>::getNonconstSolver()
00150 {
00151   return solver_;
00152 }
00153 
00154 
00155 template<class Scalar>
00156 RCP<const Thyra::NonlinearSolverBase<Scalar> >
00157 ImplicitRKStepper<Scalar>::getSolver() const
00158 {
00159   return solver_;
00160 }
00161 
00162 
00163 // Overridden from StepperBase
00164  
00165 
00166 template<class Scalar>
00167 bool ImplicitRKStepper<Scalar>::isImplicit() const
00168 {
00169   return true;
00170 }
00171 
00172 template<class Scalar>
00173 bool ImplicitRKStepper<Scalar>::supportsCloning() const
00174 {
00175   return true;
00176 }
00177 
00178 
00179 template<class Scalar>
00180 RCP<StepperBase<Scalar> >
00181 ImplicitRKStepper<Scalar>::cloneStepperAlgorithm() const
00182 {
00183   // Just use the interface to clone the algorithm in a basically
00184   // uninitialized state
00185   RCP<ImplicitRKStepper<Scalar> >
00186     stepper = Teuchos::rcp(new ImplicitRKStepper<Scalar>());
00187 
00188   if (!is_null(model_)) {
00189     stepper->setModel(model_); // Shallow copy is okay!
00190   }
00191 
00192   if (!is_null(irkButcherTableau_)) {
00193     // 06/16/09 tscoffe:  should we clone the RKBT here?
00194     stepper->setRKButcherTableau(irkButcherTableau_);
00195   }
00196 
00197   if (!is_null(solver_)) {
00198     stepper->setSolver(solver_->cloneNonlinearSolver().assert_not_null());
00199   }
00200 
00201   if (!is_null(irk_W_factory_)) {
00202     // 06/16/09 tscoffe:  should we clone the W_factory here?
00203     stepper->set_W_factory(irk_W_factory_);
00204   }
00205 
00206   if (!is_null(paramList_)) {
00207     stepper->setParameterList(Teuchos::parameterList(*paramList_));
00208   }
00209 
00210   return stepper;
00211 }
00212 
00213 
00214 template<class Scalar>
00215 void ImplicitRKStepper<Scalar>::setModel(
00216   const RCP<const Thyra::ModelEvaluator<Scalar> > & model
00217   )
00218 {
00219   TEST_FOR_EXCEPT(is_null(model));
00220   assertValidModel( *this, *model );
00221   model_ = model;
00222 }
00223 
00224 
00225 template<class Scalar>
00226 RCP<const Thyra::ModelEvaluator<Scalar> >
00227 ImplicitRKStepper<Scalar>::getModel() const
00228 {
00229   return model_;
00230 }
00231 
00232 
00233 template<class Scalar>
00234 void ImplicitRKStepper<Scalar>::setInitialCondition(
00235   const Thyra::ModelEvaluatorBase::InArgs<Scalar> &initialCondition
00236   )
00237 {
00238 
00239   typedef ScalarTraits<Scalar> ST;
00240   typedef Thyra::ModelEvaluatorBase MEB;
00241 
00242   basePoint_ = initialCondition;
00243 
00244   // x
00245 
00246   RCP<const Thyra::VectorBase<Scalar> >
00247     x_init = initialCondition.get_x();
00248 
00249 #ifdef RYTHMOS_DEBUG
00250   TEST_FOR_EXCEPTION(
00251     is_null(x_init), std::logic_error,
00252     "Error, if the client passes in an intial condition to setInitialCondition(...),\n"
00253     "then x can not be null!" );
00254 #endif
00255 
00256   x_ = x_init->clone_v();
00257 
00258   // x_dot
00259 
00260   x_dot_ = createMember(x_->space());
00261 
00262   RCP<const Thyra::VectorBase<Scalar> >
00263     x_dot_init = initialCondition.get_x_dot();
00264 
00265   if (!is_null(x_dot_init))
00266     assign(&*x_dot_,*x_dot_init);
00267   else
00268     assign(&*x_dot_,ST::zero());
00269   
00270   // t
00271 
00272   const Scalar t =
00273     (
00274       initialCondition.supports(MEB::IN_ARG_t)
00275       ? initialCondition.get_t()
00276       : ST::zero()
00277       );
00278 
00279   timeRange_ = timeRange(t,t);
00280 
00281   // x_old
00282   x_old_ = x_->clone_v();
00283 
00284   haveInitialCondition_ = true;
00285 
00286 }
00287 
00288 
00289 template<class Scalar>
00290 Scalar ImplicitRKStepper<Scalar>::takeStep(Scalar dt, StepSizeType stepSizeType)
00291 {
00292 
00293 
00294   using Teuchos::as;
00295   using Teuchos::incrVerbLevel;
00296   typedef ScalarTraits<Scalar> ST;
00297   typedef Thyra::NonlinearSolverBase<Scalar> NSB;
00298   typedef Teuchos::VerboseObjectTempState<NSB> VOTSNSB;
00299 
00300   RCP<FancyOStream> out = this->getOStream();
00301   Teuchos::EVerbosityLevel verbLevel = this->getVerbLevel();
00302   Teuchos::OSTab ostab(out,1,"takeStep");
00303   VOTSNSB solver_outputTempState(solver_,out,incrVerbLevel(verbLevel,-1));
00304 
00305   if ( !is_null(out) && as<int>(verbLevel) >= as<int>(Teuchos::VERB_LOW) ) {
00306     *out
00307       << "\nEntering " << Teuchos::TypeNameTraits<ImplicitRKStepper<Scalar> >::name()
00308       << "::takeStep("<<dt<<","<<toString(stepSizeType)<<") ...\n"; 
00309   }
00310 
00311   if (!isInitialized_) {
00312     initialize_();
00313   }
00314 
00315   TEST_FOR_EXCEPT( stepSizeType != STEP_TYPE_FIXED ); // ToDo: Handle variable case later
00316 
00317   // A) Set up the IRK ModelEvaluator so that it can represent the time step
00318   // equation to be solved.
00319 
00320   // Set irkModel_ with x_old_, t_old_, and dt
00321   V_V( &*x_old_, *x_ );
00322   Scalar current_dt = dt;
00323   Scalar t = timeRange_.upper();
00324 
00325   // B) Solve the timestep equation
00326 
00327   // Set the guess for the stage derivatives to zero (unless we can think of
00328   // something better)
00329   V_S( &*x_stage_bar_, ST::zero() );
00330 
00331   if (!isDirk_) { // General Implicit RK Case:
00332     RCP<ImplicitRKModelEvaluator<Scalar> > firkModel_ = 
00333       Teuchos::rcp_dynamic_cast<ImplicitRKModelEvaluator<Scalar> >(irkModel_,true);
00334     firkModel_->setTimeStepPoint( x_old_, t, current_dt );
00335 
00336     // Solve timestep equation
00337     solver_->solve( &*x_stage_bar_ );
00338 
00339   } else { // Diagonal Implicit RK Case:
00340 
00341     RCP<DiagonalImplicitRKModelEvaluator<Scalar> > dirkModel_ = 
00342       Teuchos::rcp_dynamic_cast<DiagonalImplicitRKModelEvaluator<Scalar> >(irkModel_,true);
00343     dirkModel_->setTimeStepPoint( x_old_, t, current_dt );
00344     int numStages = irkButcherTableau_->numStages();
00345     for (int stage=0 ; stage < numStages ; ++stage) {
00346       dirkModel_->setCurrentStage(stage);
00347       solver_->solve( &*(x_stage_bar_->getNonconstVectorBlock(stage)) );
00348       dirkModel_->setStageSolution( stage, *(x_stage_bar_->getVectorBlock(stage)) );
00349     }
00350 
00351   }
00352 
00353   // C) Complete the step ...
00354   
00355   // Combine the stage derivatives with the Butcher tableau "b" vector to obtain the solution at the final time.
00356   // x_{k+1} = x_k + dt*sum_{i}^{p}(b_i*x_stage_bar_[i])
00357 
00358   assembleIRKSolution( irkButcherTableau_->b(), current_dt, *x_old_, *x_stage_bar_,
00359     outArg(*x_)
00360     );
00361 
00362   // Update time range
00363   timeRange_ = timeRange(t,t+current_dt);
00364   numSteps_++;
00365 
00366   return current_dt;
00367 
00368 }
00369 
00370 
00371 template<class Scalar>
00372 const StepStatus<Scalar> ImplicitRKStepper<Scalar>::getStepStatus() const
00373 {
00374   StepStatus<Scalar> stepStatus;
00375 
00376   if (!isInitialized_) {
00377     stepStatus.stepStatus = STEP_STATUS_UNINITIALIZED;
00378     stepStatus.message = "This stepper is uninitialized.";
00379     return stepStatus;
00380   } 
00381   if (numSteps_ > 0) {
00382     stepStatus.stepStatus = STEP_STATUS_CONVERGED;
00383   } 
00384   else {
00385     stepStatus.stepStatus = STEP_STATUS_UNKNOWN;
00386   }
00387   stepStatus.stepSize = timeRange_.length();
00388   stepStatus.order = irkButcherTableau_->order();
00389   stepStatus.time = timeRange_.upper();
00390   stepStatus.solution = x_;
00391   return(stepStatus);
00392 }
00393 
00394 
00395 // Overridden from InterpolationBufferBase
00396 
00397 
00398 template<class Scalar>
00399 RCP<const Thyra::VectorSpaceBase<Scalar> >
00400 ImplicitRKStepper<Scalar>::get_x_space() const
00401 {
00402   return ( !is_null(model_) ? model_->get_x_space() : Teuchos::null );
00403 }
00404 
00405 
00406 template<class Scalar>
00407 void ImplicitRKStepper<Scalar>::addPoints(
00408     const Array<Scalar>& time_vec
00409     ,const Array<RCP<const Thyra::VectorBase<Scalar> > >& x_vec
00410     ,const Array<RCP<const Thyra::VectorBase<Scalar> > >& xdot_vec
00411     )
00412 {
00413   TEST_FOR_EXCEPT(true);
00414 }
00415 
00416 
00417 template<class Scalar>
00418 TimeRange<Scalar> ImplicitRKStepper<Scalar>::getTimeRange() const
00419 {
00420   return timeRange_;
00421 }
00422 
00423 
00424 template<class Scalar>
00425 void ImplicitRKStepper<Scalar>::getPoints(
00426   const Array<Scalar>& time_vec
00427   ,Array<RCP<const Thyra::VectorBase<Scalar> > >* x_vec
00428   ,Array<RCP<const Thyra::VectorBase<Scalar> > >* xdot_vec
00429   ,Array<ScalarMag>* accuracy_vec) const
00430 {
00431   using Teuchos::constOptInArg;
00432   using Teuchos::null;
00433   TEUCHOS_ASSERT(haveInitialCondition_);
00434   defaultGetPoints<Scalar>(
00435       timeRange_.lower(),constOptInArg(*x_old_),null,
00436       timeRange_.upper(),constOptInArg(*x_),null,
00437       time_vec,ptr(x_vec),ptr(xdot_vec),ptr(accuracy_vec),
00438       null
00439       );
00440   // 04/17/09 tscoffe:  Currently, we don't have x_dot to pass out (TODO)
00441 }
00442 
00443 
00444 template<class Scalar>
00445 void ImplicitRKStepper<Scalar>::getNodes(Array<Scalar>* time_vec) const
00446 {
00447   TEUCHOS_ASSERT( time_vec != NULL );
00448   time_vec->clear();
00449   if (!haveInitialCondition_) {
00450     return;
00451   }
00452   time_vec->push_back(timeRange_.lower());
00453   if (numSteps_ > 0) {
00454     time_vec->push_back(timeRange_.upper());
00455   }
00456 }
00457 
00458 
00459 template<class Scalar>
00460 void ImplicitRKStepper<Scalar>::removeNodes(Array<Scalar>& time_vec) 
00461 {
00462   TEST_FOR_EXCEPT(true);
00463 }
00464 
00465 
00466 template<class Scalar>
00467 int ImplicitRKStepper<Scalar>::getOrder() const
00468 {
00469   return irkButcherTableau_->order();
00470 }
00471 
00472 
00473 // Overridden from Teuchos::ParameterListAcceptor
00474 
00475 
00476 template <class Scalar>
00477 void ImplicitRKStepper<Scalar>::setParameterList(
00478   RCP<ParameterList> const& paramList
00479   )
00480 {
00481   TEST_FOR_EXCEPT(is_null(paramList));
00482   paramList->validateParametersAndSetDefaults(*this->getValidParameters());
00483   paramList_ = paramList;
00484   Teuchos::readVerboseObjectSublist(&*paramList_,this);
00485 }
00486 
00487 
00488 template <class Scalar>
00489 RCP<ParameterList>
00490 ImplicitRKStepper<Scalar>::getNonconstParameterList()
00491 {
00492   return(paramList_);
00493 }
00494 
00495 
00496 template <class Scalar>
00497 RCP<ParameterList>
00498 ImplicitRKStepper<Scalar>::unsetParameterList()
00499 {
00500   RCP<ParameterList>
00501     temp_param_list = paramList_;
00502   paramList_ = Teuchos::null;
00503   return(temp_param_list);
00504 }
00505 
00506 
00507 template<class Scalar>
00508 RCP<const ParameterList>
00509 ImplicitRKStepper<Scalar>::getValidParameters() const
00510 {
00511   static RCP<const ParameterList> validPL;
00512   if (is_null(validPL)) {
00513     RCP<ParameterList> pl = Teuchos::parameterList();
00514     Teuchos::setupVerboseObjectSublist(&*pl);
00515     validPL = pl;
00516   }
00517   return validPL;
00518 }
00519 
00520 
00521 // Overridden from Teuchos::Describable
00522 
00523 
00524 template<class Scalar>
00525 void ImplicitRKStepper<Scalar>::describe(
00526   FancyOStream &out,
00527   const Teuchos::EVerbosityLevel verbLevel
00528   ) const
00529 {
00530   using std::endl;
00531   using Teuchos::as;
00532   Teuchos::OSTab tab(out);
00533   if (!isInitialized_) {
00534     out << this->description() << " : This stepper is not initialized yet" << std::endl;
00535     return;
00536   }
00537   if (
00538     as<int>(verbLevel) == as<int>(Teuchos::VERB_DEFAULT)
00539     ||
00540     as<int>(verbLevel) >= as<int>(Teuchos::VERB_LOW)
00541     )
00542   {
00543     out << this->description() << ":" << endl;
00544     Teuchos::OSTab tab(out);
00545     out << "model = " << Teuchos::describe(*model_,verbLevel);
00546     out << "solver = " << Teuchos::describe(*solver_,verbLevel);
00547     out << "irk_W_factory = " << Teuchos::describe(*irk_W_factory_,verbLevel);
00548   }
00549 }
00550 
00551 
00552 // private
00553 
00554 
00555 template <class Scalar>
00556 void ImplicitRKStepper<Scalar>::initialize_()
00557 {
00558 
00559   typedef ScalarTraits<Scalar> ST;
00560   using Teuchos::rcp_dynamic_cast;
00561 
00562   TEST_FOR_EXCEPT(is_null(model_));
00563   TEST_FOR_EXCEPT(is_null(solver_));
00564   TEST_FOR_EXCEPT(irkButcherTableau_->numStages() == 0);
00565   TEUCHOS_ASSERT(haveInitialCondition_);
00566 
00567 #ifdef RYTHMOS_DEBUG
00568   THYRA_ASSERT_VEC_SPACES(
00569     "Rythmos::ImplicitRKStepper::initialize_(...)",
00570     *x_->space(), *model_->get_x_space() );
00571 #endif
00572 
00573 
00574   // Set up the IRK mdoel
00575 
00576   if (!isDirk_) { // General Implicit RK 
00577     TEST_FOR_EXCEPT(is_null(irk_W_factory_));
00578     irkModel_ = implicitRKModelEvaluator(
00579       model_,basePoint_,irk_W_factory_,irkButcherTableau_);
00580   } else { // Diagonal Implicit RK
00581     irkModel_ = diagonalImplicitRKModelEvaluator(
00582       model_,basePoint_,irk_W_factory_,irkButcherTableau_);
00583   }
00584 
00585   solver_->setModel(irkModel_);
00586 
00587   // Set up the vector of stage derivatives ...
00588   const int numStages = irkButcherTableau_->numStages();
00589   RCP<const Thyra::ProductVectorSpaceBase<Scalar> > pvs = productVectorSpace(model_->get_x_space(),numStages);
00590   RCP<const Thyra::VectorSpaceBase<Scalar> > vs = rcp_dynamic_cast<const Thyra::VectorSpaceBase<Scalar> >(pvs,true);
00591   x_stage_bar_ = rcp_dynamic_cast<Thyra::ProductVectorBase<Scalar> >(createMember(vs),true);
00592 //  x_stage_bar_ = rcp_dynamic_cast<Thyra::ProductVectorBase<Scalar> >(
00593 //    createMember(irkModel_->get_x_space())
00594 //    );
00595 
00596   isInitialized_ = true;
00597 
00598 }
00599 
00600 template <class Scalar>
00601 void ImplicitRKStepper<Scalar>::setRKButcherTableau( const RCP<const RKButcherTableauBase<Scalar> > &rkButcherTableau )
00602 {
00603   TEUCHOS_ASSERT( !is_null(rkButcherTableau) );
00604   TEST_FOR_EXCEPTION( isInitialized_, std::logic_error,
00605       "Error!  The RK Butcher Tableau cannot be changed after internal initialization!"
00606       );
00607   validateIRKButcherTableau(*rkButcherTableau);
00608   irkButcherTableau_ = rkButcherTableau;
00609   E_RKButcherTableauTypes rkType = determineRKBTType<Scalar>(*irkButcherTableau_);
00610   if (
00611          (rkType == RYTHMOS_RK_BUTCHER_TABLEAU_TYPE_DIRK) 
00612       || (rkType == RYTHMOS_RK_BUTCHER_TABLEAU_TYPE_SDIRK) 
00613       || (irkButcherTableau_->numStages() == 1)
00614      ) 
00615   {
00616     isDirk_ = true;
00617   } 
00618 }
00619 
00620 template <class Scalar>
00621 RCP<const RKButcherTableauBase<Scalar> > ImplicitRKStepper<Scalar>::getRKButcherTableau() const
00622 {
00623   return irkButcherTableau_;
00624 }
00625 
00626 template<class Scalar>
00627 void ImplicitRKStepper<Scalar>::setDirk(bool isDirk)
00628 {
00629   TEST_FOR_EXCEPTION(isInitialized_, std::logic_error,
00630       "Error!  Cannot change DIRK flag after internal initialization is completed\n"
00631       );
00632   if (isDirk == true) {
00633     E_RKButcherTableauTypes rkType = determineRKBTType<Scalar>(*irkButcherTableau_);
00634     bool RKBT_is_DIRK = (
00635          (rkType == RYTHMOS_RK_BUTCHER_TABLEAU_TYPE_DIRK) 
00636       || (rkType == RYTHMOS_RK_BUTCHER_TABLEAU_TYPE_SDIRK) 
00637       || (irkButcherTableau_->numStages() == 1)
00638       );
00639     TEST_FOR_EXCEPTION( !RKBT_is_DIRK, std::logic_error,
00640         "Error!  Cannot set DIRK flag on a non-DIRK RK Butcher Tableau\n"
00641         );
00642   } else { // isDirk = false;
00643     isDirk_ = isDirk;
00644   }
00645 }
00646 
00647 // 
00648 // Explicit Instantiation macro
00649 //
00650 // Must be expanded from within the Rythmos namespace!
00651 //
00652 
00653 #define RYTHMOS_IMPLICIT_RK_STEPPER_INSTANT(SCALAR) \
00654   \
00655   template class ImplicitRKStepper< SCALAR >; \
00656   \
00657   template RCP< ImplicitRKStepper< SCALAR > > \
00658   implicitRKStepper();  \
00659   \
00660   template RCP< ImplicitRKStepper< SCALAR > > \
00661   implicitRKStepper( \
00662     const RCP<const Thyra::ModelEvaluator< SCALAR > > &model, \
00663     const RCP<Thyra::NonlinearSolverBase< SCALAR > >  &solver, \
00664     const RCP<Thyra::LinearOpWithSolveFactoryBase< SCALAR > > &irk_W_factory, \
00665     const RCP<const RKButcherTableauBase< SCALAR > > &irkbt \
00666       ); 
00667 
00668 } // namespace Rythmos
00669 
00670 #endif //Rythmos_IMPLICIT_RK_STEPPER_DEF_H

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