Rythmos - Transient Integration for Differential Equations Version of the Day
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   TEUCHOS_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   TEUCHOS_TEST_FOR_EXCEPT(is_null(model));
00220   assertValidModel( *this, *model );
00221   model_ = model;
00222 }
00223 
00224 
00225 template<class Scalar>
00226 void ImplicitRKStepper<Scalar>::setNonconstModel(
00227   const RCP<Thyra::ModelEvaluator<Scalar> >& model
00228   )
00229 {
00230   this->setModel(model); // TODO 09/09/09 tscoffe:  use ConstNonconstObjectContainer!
00231 }
00232 
00233 
00234 template<class Scalar>
00235 RCP<const Thyra::ModelEvaluator<Scalar> >
00236 ImplicitRKStepper<Scalar>::getModel() const
00237 {
00238   return model_;
00239 }
00240 
00241 
00242 template<class Scalar>
00243 RCP<Thyra::ModelEvaluator<Scalar> >
00244 ImplicitRKStepper<Scalar>::getNonconstModel() 
00245 {
00246   return Teuchos::null;
00247 }
00248 
00249 
00250 template<class Scalar>
00251 void ImplicitRKStepper<Scalar>::setInitialCondition(
00252   const Thyra::ModelEvaluatorBase::InArgs<Scalar> &initialCondition
00253   )
00254 {
00255 
00256   typedef ScalarTraits<Scalar> ST;
00257   typedef Thyra::ModelEvaluatorBase MEB;
00258 
00259   basePoint_ = initialCondition;
00260 
00261   // x
00262 
00263   RCP<const Thyra::VectorBase<Scalar> >
00264     x_init = initialCondition.get_x();
00265 
00266 #ifdef HAVE_RYTHMOS_DEBUG
00267   TEUCHOS_TEST_FOR_EXCEPTION(
00268     is_null(x_init), std::logic_error,
00269     "Error, if the client passes in an intial condition to setInitialCondition(...),\n"
00270     "then x can not be null!" );
00271 #endif
00272 
00273   x_ = x_init->clone_v();
00274 
00275   // x_dot
00276 
00277   x_dot_ = createMember(x_->space());
00278 
00279   RCP<const Thyra::VectorBase<Scalar> >
00280     x_dot_init = initialCondition.get_x_dot();
00281 
00282   if (!is_null(x_dot_init))
00283     assign(x_dot_.ptr(),*x_dot_init);
00284   else
00285     assign(x_dot_.ptr(),ST::zero());
00286   
00287   // t
00288 
00289   const Scalar t =
00290     (
00291       initialCondition.supports(MEB::IN_ARG_t)
00292       ? initialCondition.get_t()
00293       : ST::zero()
00294       );
00295 
00296   timeRange_ = timeRange(t,t);
00297 
00298   // x_old
00299   x_old_ = x_->clone_v();
00300 
00301   haveInitialCondition_ = true;
00302 
00303 }
00304 
00305 
00306 template<class Scalar>
00307 Thyra::ModelEvaluatorBase::InArgs<Scalar> 
00308 ImplicitRKStepper<Scalar>::getInitialCondition() const
00309 {
00310   return basePoint_;
00311 }
00312 
00313 
00314 template<class Scalar>
00315 Scalar ImplicitRKStepper<Scalar>::takeStep(Scalar dt, StepSizeType stepSizeType)
00316 {
00317 
00318 
00319   using Teuchos::as;
00320   using Teuchos::incrVerbLevel;
00321   typedef ScalarTraits<Scalar> ST;
00322   typedef Thyra::NonlinearSolverBase<Scalar> NSB;
00323   typedef Teuchos::VerboseObjectTempState<NSB> VOTSNSB;
00324 
00325   RCP<FancyOStream> out = this->getOStream();
00326   Teuchos::EVerbosityLevel verbLevel = this->getVerbLevel();
00327   Teuchos::OSTab ostab(out,1,"takeStep");
00328   VOTSNSB solver_outputTempState(solver_,out,incrVerbLevel(verbLevel,-1));
00329 
00330   if ( !is_null(out) && as<int>(verbLevel) >= as<int>(Teuchos::VERB_LOW) ) {
00331     *out
00332       << "\nEntering " << Teuchos::TypeNameTraits<ImplicitRKStepper<Scalar> >::name()
00333       << "::takeStep("<<dt<<","<<toString(stepSizeType)<<") ...\n"; 
00334   }
00335 
00336   if (!isInitialized_) {
00337     initialize_();
00338   }
00339 
00340   TEUCHOS_TEST_FOR_EXCEPT( stepSizeType != STEP_TYPE_FIXED ); // ToDo: Handle variable case later
00341 
00342   // A) Set up the IRK ModelEvaluator so that it can represent the time step
00343   // equation to be solved.
00344 
00345   // Set irkModel_ with x_old_, t_old_, and dt
00346   V_V( x_old_.ptr(), *x_ );
00347   Scalar current_dt = dt;
00348   Scalar t = timeRange_.upper();
00349 
00350   // B) Solve the timestep equation
00351 
00352   // Set the guess for the stage derivatives to zero (unless we can think of
00353   // something better)
00354   V_S( Teuchos::rcp_dynamic_cast<Thyra::VectorBase<Scalar> >(x_stage_bar_).ptr(), ST::zero() );
00355 
00356   if (!isDirk_) { // General Implicit RK Case:
00357     RCP<ImplicitRKModelEvaluator<Scalar> > firkModel_ = 
00358       Teuchos::rcp_dynamic_cast<ImplicitRKModelEvaluator<Scalar> >(irkModel_,true);
00359     firkModel_->setTimeStepPoint( x_old_, t, current_dt );
00360 
00361     // Solve timestep equation
00362     solver_->solve( &*x_stage_bar_ );
00363 
00364   } else { // Diagonal Implicit RK Case:
00365 
00366     RCP<DiagonalImplicitRKModelEvaluator<Scalar> > dirkModel_ = 
00367       Teuchos::rcp_dynamic_cast<DiagonalImplicitRKModelEvaluator<Scalar> >(irkModel_,true);
00368     dirkModel_->setTimeStepPoint( x_old_, t, current_dt );
00369     int numStages = irkButcherTableau_->numStages();
00370     for (int stage=0 ; stage < numStages ; ++stage) {
00371       dirkModel_->setCurrentStage(stage);
00372       solver_->solve( &*(x_stage_bar_->getNonconstVectorBlock(stage)) );
00373       dirkModel_->setStageSolution( stage, *(x_stage_bar_->getVectorBlock(stage)) );
00374     }
00375 
00376   }
00377 
00378   // C) Complete the step ...
00379   
00380   // Combine the stage derivatives with the Butcher tableau "b" vector to obtain the solution at the final time.
00381   // x_{k+1} = x_k + dt*sum_{i}^{p}(b_i*x_stage_bar_[i])
00382 
00383   assembleIRKSolution( irkButcherTableau_->b(), current_dt, *x_old_, *x_stage_bar_,
00384     outArg(*x_)
00385     );
00386 
00387   // Update time range
00388   timeRange_ = timeRange(t,t+current_dt);
00389   numSteps_++;
00390 
00391   return current_dt;
00392 
00393 }
00394 
00395 
00396 template<class Scalar>
00397 const StepStatus<Scalar> ImplicitRKStepper<Scalar>::getStepStatus() const
00398 {
00399   StepStatus<Scalar> stepStatus;
00400 
00401   if (!isInitialized_) {
00402     stepStatus.stepStatus = STEP_STATUS_UNINITIALIZED;
00403     stepStatus.message = "This stepper is uninitialized.";
00404     return stepStatus;
00405   } 
00406   if (numSteps_ > 0) {
00407     stepStatus.stepStatus = STEP_STATUS_CONVERGED;
00408   } 
00409   else {
00410     stepStatus.stepStatus = STEP_STATUS_UNKNOWN;
00411   }
00412   stepStatus.stepSize = timeRange_.length();
00413   stepStatus.order = irkButcherTableau_->order();
00414   stepStatus.time = timeRange_.upper();
00415   stepStatus.solution = x_;
00416   return(stepStatus);
00417 }
00418 
00419 
00420 // Overridden from InterpolationBufferBase
00421 
00422 
00423 template<class Scalar>
00424 RCP<const Thyra::VectorSpaceBase<Scalar> >
00425 ImplicitRKStepper<Scalar>::get_x_space() const
00426 {
00427   return ( !is_null(model_) ? model_->get_x_space() : Teuchos::null );
00428 }
00429 
00430 
00431 template<class Scalar>
00432 void ImplicitRKStepper<Scalar>::addPoints(
00433     const Array<Scalar>& time_vec
00434     ,const Array<RCP<const Thyra::VectorBase<Scalar> > >& x_vec
00435     ,const Array<RCP<const Thyra::VectorBase<Scalar> > >& xdot_vec
00436     )
00437 {
00438   TEUCHOS_TEST_FOR_EXCEPT(true);
00439 }
00440 
00441 
00442 template<class Scalar>
00443 TimeRange<Scalar> ImplicitRKStepper<Scalar>::getTimeRange() const
00444 {
00445   return timeRange_;
00446 }
00447 
00448 
00449 template<class Scalar>
00450 void ImplicitRKStepper<Scalar>::getPoints(
00451   const Array<Scalar>& time_vec
00452   ,Array<RCP<const Thyra::VectorBase<Scalar> > >* x_vec
00453   ,Array<RCP<const Thyra::VectorBase<Scalar> > >* xdot_vec
00454   ,Array<ScalarMag>* accuracy_vec) const
00455 {
00456   using Teuchos::constOptInArg;
00457   using Teuchos::null;
00458   TEUCHOS_ASSERT(haveInitialCondition_);
00459   defaultGetPoints<Scalar>(
00460     timeRange_.lower(), constOptInArg(*x_old_),
00461     Ptr<const VectorBase<Scalar> >(null), // Sun
00462     timeRange_.upper(), constOptInArg(*x_),
00463     Ptr<const VectorBase<Scalar> >(null), // Sun
00464     time_vec,
00465     ptr(x_vec), ptr(xdot_vec), ptr(accuracy_vec),
00466     Ptr<InterpolatorBase<Scalar> >(null) // For Sun
00467     );
00468   // 04/17/09 tscoffe:  Currently, we don't have x_dot to pass out (TODO)
00469 }
00470 
00471 
00472 template<class Scalar>
00473 void ImplicitRKStepper<Scalar>::getNodes(Array<Scalar>* time_vec) const
00474 {
00475   TEUCHOS_ASSERT( time_vec != NULL );
00476   time_vec->clear();
00477   if (!haveInitialCondition_) {
00478     return;
00479   }
00480   time_vec->push_back(timeRange_.lower());
00481   if (numSteps_ > 0) {
00482     time_vec->push_back(timeRange_.upper());
00483   }
00484 }
00485 
00486 
00487 template<class Scalar>
00488 void ImplicitRKStepper<Scalar>::removeNodes(Array<Scalar>& time_vec) 
00489 {
00490   TEUCHOS_TEST_FOR_EXCEPT(true);
00491 }
00492 
00493 
00494 template<class Scalar>
00495 int ImplicitRKStepper<Scalar>::getOrder() const
00496 {
00497   return irkButcherTableau_->order();
00498 }
00499 
00500 
00501 // Overridden from Teuchos::ParameterListAcceptor
00502 
00503 
00504 template <class Scalar>
00505 void ImplicitRKStepper<Scalar>::setParameterList(
00506   RCP<ParameterList> const& paramList
00507   )
00508 {
00509   TEUCHOS_TEST_FOR_EXCEPT(is_null(paramList));
00510   paramList->validateParametersAndSetDefaults(*this->getValidParameters());
00511   paramList_ = paramList;
00512   Teuchos::readVerboseObjectSublist(&*paramList_,this);
00513 }
00514 
00515 
00516 template <class Scalar>
00517 RCP<ParameterList>
00518 ImplicitRKStepper<Scalar>::getNonconstParameterList()
00519 {
00520   return(paramList_);
00521 }
00522 
00523 
00524 template <class Scalar>
00525 RCP<ParameterList>
00526 ImplicitRKStepper<Scalar>::unsetParameterList()
00527 {
00528   RCP<ParameterList>
00529     temp_param_list = paramList_;
00530   paramList_ = Teuchos::null;
00531   return(temp_param_list);
00532 }
00533 
00534 
00535 template<class Scalar>
00536 RCP<const ParameterList>
00537 ImplicitRKStepper<Scalar>::getValidParameters() const
00538 {
00539   static RCP<const ParameterList> validPL;
00540   if (is_null(validPL)) {
00541     RCP<ParameterList> pl = Teuchos::parameterList();
00542     Teuchos::setupVerboseObjectSublist(&*pl);
00543     validPL = pl;
00544   }
00545   return validPL;
00546 }
00547 
00548 
00549 // Overridden from Teuchos::Describable
00550 
00551 
00552 template<class Scalar>
00553 void ImplicitRKStepper<Scalar>::describe(
00554   FancyOStream &out,
00555   const Teuchos::EVerbosityLevel verbLevel
00556   ) const
00557 {
00558   using std::endl;
00559   using Teuchos::as;
00560   if (!isInitialized_) {
00561     out << this->description() << " : This stepper is not initialized yet" << std::endl;
00562     return;
00563   }
00564   if (
00565     as<int>(verbLevel) == as<int>(Teuchos::VERB_DEFAULT)
00566     ||
00567     as<int>(verbLevel) >= as<int>(Teuchos::VERB_LOW)
00568     )
00569   {
00570     out << this->description() << ":" << endl;
00571     Teuchos::OSTab tab(out);
00572     out << "model = " << Teuchos::describe(*model_,verbLevel);
00573     out << "solver = " << Teuchos::describe(*solver_,verbLevel);
00574     out << "irk_W_factory = " << Teuchos::describe(*irk_W_factory_,verbLevel);
00575   }
00576 }
00577 
00578 
00579 // private
00580 
00581 
00582 template <class Scalar>
00583 void ImplicitRKStepper<Scalar>::initialize_()
00584 {
00585 
00586   typedef ScalarTraits<Scalar> ST;
00587   using Teuchos::rcp_dynamic_cast;
00588 
00589   TEUCHOS_TEST_FOR_EXCEPT(is_null(model_));
00590   TEUCHOS_TEST_FOR_EXCEPT(is_null(solver_));
00591   TEUCHOS_TEST_FOR_EXCEPT(irkButcherTableau_->numStages() == 0);
00592   TEUCHOS_ASSERT(haveInitialCondition_);
00593 
00594 #ifdef HAVE_RYTHMOS_DEBUG
00595   THYRA_ASSERT_VEC_SPACES(
00596     "Rythmos::ImplicitRKStepper::initialize_(...)",
00597     *x_->space(), *model_->get_x_space() );
00598 #endif
00599 
00600 
00601   // Set up the IRK mdoel
00602 
00603   if (!isDirk_) { // General Implicit RK 
00604     TEUCHOS_TEST_FOR_EXCEPT(is_null(irk_W_factory_));
00605     irkModel_ = implicitRKModelEvaluator(
00606       model_,basePoint_,irk_W_factory_,irkButcherTableau_);
00607   } else { // Diagonal Implicit RK
00608     irkModel_ = diagonalImplicitRKModelEvaluator(
00609       model_,basePoint_,irk_W_factory_,irkButcherTableau_);
00610   }
00611 
00612   solver_->setModel(irkModel_);
00613 
00614   // Set up the vector of stage derivatives ...
00615   const int numStages = irkButcherTableau_->numStages();
00616   RCP<const Thyra::ProductVectorSpaceBase<Scalar> > pvs = productVectorSpace(model_->get_x_space(),numStages);
00617   RCP<const Thyra::VectorSpaceBase<Scalar> > vs = rcp_dynamic_cast<const Thyra::VectorSpaceBase<Scalar> >(pvs,true);
00618   x_stage_bar_ = rcp_dynamic_cast<Thyra::ProductVectorBase<Scalar> >(createMember(vs),true);
00619 //  x_stage_bar_ = rcp_dynamic_cast<Thyra::ProductVectorBase<Scalar> >(
00620 //    createMember(irkModel_->get_x_space())
00621 //    );
00622 
00623   isInitialized_ = true;
00624 
00625 }
00626 
00627 template <class Scalar>
00628 void ImplicitRKStepper<Scalar>::setRKButcherTableau( const RCP<const RKButcherTableauBase<Scalar> > &rkButcherTableau )
00629 {
00630   TEUCHOS_ASSERT( !is_null(rkButcherTableau) );
00631   TEUCHOS_TEST_FOR_EXCEPTION( isInitialized_, std::logic_error,
00632       "Error!  The RK Butcher Tableau cannot be changed after internal initialization!"
00633       );
00634   validateIRKButcherTableau(*rkButcherTableau);
00635   irkButcherTableau_ = rkButcherTableau;
00636   E_RKButcherTableauTypes rkType = determineRKBTType<Scalar>(*irkButcherTableau_);
00637   if (
00638          (rkType == RYTHMOS_RK_BUTCHER_TABLEAU_TYPE_DIRK) 
00639       || (rkType == RYTHMOS_RK_BUTCHER_TABLEAU_TYPE_SDIRK) 
00640       || (irkButcherTableau_->numStages() == 1)
00641      ) 
00642   {
00643     isDirk_ = true;
00644   } 
00645 }
00646 
00647 template <class Scalar>
00648 RCP<const RKButcherTableauBase<Scalar> > ImplicitRKStepper<Scalar>::getRKButcherTableau() const
00649 {
00650   return irkButcherTableau_;
00651 }
00652 
00653 template<class Scalar>
00654 void ImplicitRKStepper<Scalar>::setDirk(bool isDirk)
00655 {
00656   TEUCHOS_TEST_FOR_EXCEPTION(isInitialized_, std::logic_error,
00657       "Error!  Cannot change DIRK flag after internal initialization is completed\n"
00658       );
00659   if (isDirk == true) {
00660     E_RKButcherTableauTypes rkType = determineRKBTType<Scalar>(*irkButcherTableau_);
00661     bool RKBT_is_DIRK = (
00662          (rkType == RYTHMOS_RK_BUTCHER_TABLEAU_TYPE_DIRK) 
00663       || (rkType == RYTHMOS_RK_BUTCHER_TABLEAU_TYPE_SDIRK) 
00664       || (irkButcherTableau_->numStages() == 1)
00665       );
00666     TEUCHOS_TEST_FOR_EXCEPTION( !RKBT_is_DIRK, std::logic_error,
00667         "Error!  Cannot set DIRK flag on a non-DIRK RK Butcher Tableau\n"
00668         );
00669   } else { // isDirk = false;
00670     isDirk_ = isDirk;
00671   }
00672 }
00673 
00674 // 
00675 // Explicit Instantiation macro
00676 //
00677 // Must be expanded from within the Rythmos namespace!
00678 //
00679 
00680 #define RYTHMOS_IMPLICIT_RK_STEPPER_INSTANT(SCALAR) \
00681   \
00682   template class ImplicitRKStepper< SCALAR >; \
00683   \
00684   template RCP< ImplicitRKStepper< SCALAR > > \
00685   implicitRKStepper();  \
00686   \
00687   template RCP< ImplicitRKStepper< SCALAR > > \
00688   implicitRKStepper( \
00689     const RCP<const Thyra::ModelEvaluator< SCALAR > >& model, \
00690     const RCP<Thyra::NonlinearSolverBase< SCALAR > >& solver, \
00691     const RCP<Thyra::LinearOpWithSolveFactoryBase< SCALAR > >& irk_W_factory, \
00692     const RCP<const RKButcherTableauBase< SCALAR > >& irkbt \
00693       ); 
00694 
00695 } // namespace Rythmos
00696 
00697 #endif //Rythmos_IMPLICIT_RK_STEPPER_DEF_H
 All Classes Functions Variables Typedefs Friends