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