Tempus_StepperRKButcherTableau.hpp

Go to the documentation of this file.

// @HEADER
// ****************************************************************************
//                Tempus: Copyright (2017) Sandia Corporation
//
// Distributed under BSD 3-clause license (See accompanying file Copyright.txt)
// ****************************************************************************
// @HEADER

#ifndef Tempus_StepperRKButcherTableau_hpp
#define Tempus_StepperRKButcherTableau_hpp

// disable clang warnings
#ifdef __clang__
#pragma clang system_header
#endif

#include "Tempus_config.hpp"
#include "Tempus_StepperExplicitRK.hpp"
#include "Tempus_StepperDIRK.hpp"
#include "Tempus_RKButcherTableau.hpp"


namespace Tempus {


// ----------------------------------------------------------------------------
template<class Scalar>
class StepperERK_ForwardEuler :
  virtual public StepperExplicitRK<Scalar>
{
public:
  StepperERK_ForwardEuler()
  {
    this->setStepperName("RK Forward Euler");
    this->setStepperType("RK Forward Euler");
    this->setupTableau();
    this->setupDefault();
    this->setUseFSAL(            true);
    this->setICConsistency(      "Consistent");
    this->setICConsistencyCheck( false);
  }

  StepperERK_ForwardEuler(
    const Teuchos::RCP<const Thyra::ModelEvaluator<Scalar> >& appModel,
    bool useFSAL,
    std::string ICConsistency,
    bool ICConsistencyCheck,
    bool useEmbedded,
    const Teuchos::RCP<StepperRKAppAction<Scalar> >& stepperRKAppAction)
  {
    this->setStepperName("RK Forward Euler");
    this->setStepperType("RK Forward Euler");
    this->setupTableau();
    this->setup(appModel, useFSAL, ICConsistency,
                ICConsistencyCheck, useEmbedded, stepperRKAppAction);
  }

  std::string getDescription() const
  {
    std::ostringstream Description;
    Description << this->getStepperType() << "\n"
                << "c = [ 0 ]'\n"
                << "A = [ 0 ]\n"
                << "b = [ 1 ]'";
    return Description.str();
  }

  void setUseFSAL(bool a) { this->useFSAL_ = a; this->isInitialized_ = false; }

protected:

  void setupTableau()
  {
    typedef Teuchos::ScalarTraits<Scalar> ST;
    Teuchos::SerialDenseMatrix<int,Scalar> A(1,1);
    Teuchos::SerialDenseVector<int,Scalar> b(1);
    Teuchos::SerialDenseVector<int,Scalar> c(1);
    A(0,0) = ST::zero();
    b(0) = ST::one();
    c(0) = ST::zero();
    int order = 1;

    this->tableau_ = Teuchos::rcp(new RKButcherTableau<Scalar>(
      this->getStepperType(),A,b,c,order,order,order));
    this->tableau_->setTVD(true);
    this->tableau_->setTVDCoeff(1.0);
  }
};


// ------------------------------------------------------------------------
template<class Scalar>
Teuchos::RCP<StepperERK_ForwardEuler<Scalar> >
createStepperERK_ForwardEuler(
  const Teuchos::RCP<const Thyra::ModelEvaluator<Scalar> >& model,
  Teuchos::RCP<Teuchos::ParameterList> pl)
{
  auto stepper = Teuchos::rcp(new StepperERK_ForwardEuler<Scalar>());

  // Test for aliases.
  if (pl != Teuchos::null) {
    auto stepperType =
      pl->get<std::string>("Stepper Type", stepper->getStepperType());

    TEUCHOS_TEST_FOR_EXCEPTION(
      stepperType != stepper->getStepperType() &&
      stepperType != "RK1", std::logic_error,
      "  ParameterList 'Stepper Type' (='" + stepperType +"')\n"
      "  does not match type for this Stepper (='" + stepper->getStepperType() +
      "')\n  or one of its aliases ('RK1').\n");

    // Reset default StepperType.
    pl->set<std::string>("Stepper Type", stepper->getStepperType());
  }

  stepper->setStepperRKValues(pl);

  if (model != Teuchos::null) {
    stepper->setModel(model);
    stepper->initialize();
  }

  return stepper;
}


// ----------------------------------------------------------------------------
template<class Scalar>
class StepperERK_4Stage4thOrder :
  virtual public StepperExplicitRK<Scalar>
{
public:
  StepperERK_4Stage4thOrder()
  {
    this->setStepperName("RK Explicit 4 Stage");
    this->setStepperType("RK Explicit 4 Stage");
    this->setupTableau();
    this->setupDefault();
    this->setUseFSAL(            false);
    this->setICConsistency(      "None");
    this->setICConsistencyCheck( false);
  }

  StepperERK_4Stage4thOrder(
    const Teuchos::RCP<const Thyra::ModelEvaluator<Scalar> >& appModel,
    bool useFSAL,
    std::string ICConsistency,
    bool ICConsistencyCheck,
    bool useEmbedded,
    const Teuchos::RCP<StepperRKAppAction<Scalar> >& stepperRKAppAction)
  {
    this->setStepperName("RK Explicit 4 Stage");
    this->setStepperType("RK Explicit 4 Stage");
    this->setupTableau();
    this->setup(appModel, useFSAL, ICConsistency,
                ICConsistencyCheck, useEmbedded, stepperRKAppAction);
  }

  std::string getDescription() const
  {
    std::ostringstream Description;
    Description << this->getStepperType() << "\n"
                << "\"The\" Runge-Kutta Method (explicit):\n"
                << "Solving Ordinary Differential Equations I:\n"
                << "Nonstiff Problems, 2nd Revised Edition\n"
                << "E. Hairer, S.P. Norsett, G. Wanner\n"
                << "Table 1.2, pg 138\n"
                << "c = [  0  1/2 1/2  1  ]'\n"
                << "A = [  0              ] \n"
                << "    [ 1/2  0          ]\n"
                << "    [  0  1/2  0      ]\n"
                << "    [  0   0   1   0  ]\n"
                << "b = [ 1/6 1/3 1/3 1/6 ]'";
    return Description.str();
  }

  void setupTableau()
  {
    typedef Teuchos::ScalarTraits<Scalar> ST;
    const Scalar one = ST::one();
    const Scalar zero = ST::zero();
    const Scalar onehalf = one/(2*one);
    const Scalar onesixth = one/(6*one);
    const Scalar onethird = one/(3*one);

    int NumStages = 4;
    Teuchos::SerialDenseMatrix<int,Scalar> A(NumStages,NumStages);
    Teuchos::SerialDenseVector<int,Scalar> b(NumStages);
    Teuchos::SerialDenseVector<int,Scalar> c(NumStages);

    // Fill A:
    A(0,0) =    zero; A(0,1) =    zero; A(0,2) = zero; A(0,3) = zero;
    A(1,0) = onehalf; A(1,1) =    zero; A(1,2) = zero; A(1,3) = zero;
    A(2,0) =    zero; A(2,1) = onehalf; A(2,2) = zero; A(2,3) = zero;
    A(3,0) =    zero; A(3,1) =    zero; A(3,2) =  one; A(3,3) = zero;

    // Fill b:
    b(0) = onesixth; b(1) = onethird; b(2) = onethird; b(3) = onesixth;

    // fill c:
    c(0) = zero; c(1) = onehalf; c(2) = onehalf; c(3) = one;

    int order = 4;

    this->tableau_ = Teuchos::rcp(new RKButcherTableau<Scalar>(
      this->getStepperType(),A,b,c,order,order,order));
  }
};


// ------------------------------------------------------------------------
template<class Scalar>
Teuchos::RCP<StepperERK_4Stage4thOrder<Scalar> >
createStepperERK_4Stage4thOrder(
  const Teuchos::RCP<const Thyra::ModelEvaluator<Scalar> >& model,
  Teuchos::RCP<Teuchos::ParameterList> pl)
{
  auto stepper = Teuchos::rcp(new StepperERK_4Stage4thOrder<Scalar>());
  stepper->setStepperRKValues(pl);

  if (model != Teuchos::null) {
    stepper->setModel(model);
    stepper->initialize();
  }

  return stepper;
}


// ----------------------------------------------------------------------------
template<class Scalar>
class StepperERK_BogackiShampine32 :
  virtual public StepperExplicitRK<Scalar>
{
public:
  StepperERK_BogackiShampine32()
  {
    this->setStepperName("Bogacki-Shampine 3(2) Pair");
    this->setStepperType("Bogacki-Shampine 3(2) Pair");
    this->setupTableau();
    this->setupDefault();
    this->setUseFSAL(            true);
    this->setICConsistency(      "Consistent");
    this->setICConsistencyCheck( false);
  }

  StepperERK_BogackiShampine32(
    const Teuchos::RCP<const Thyra::ModelEvaluator<Scalar> >& appModel,
    bool useFSAL,
    std::string ICConsistency,
    bool ICConsistencyCheck,
    bool useEmbedded,
    const Teuchos::RCP<StepperRKAppAction<Scalar> >& stepperRKAppAction)
  {
    this->setStepperName("Bogacki-Shampine 3(2) Pair");
    this->setStepperType("Bogacki-Shampine 3(2) Pair");
    this->setupTableau();
    this->setup(appModel, useFSAL, ICConsistency,
                ICConsistencyCheck, useEmbedded, stepperRKAppAction);
  }

  std::string getDescription() const
  {
    std::ostringstream Description;
    Description << this->getStepperType() << "\n"
                << "P. Bogacki and L.F. Shampine.\n"
                << "A 3(2) pair of Runge–Kutta formulas.\n"
                << "Applied Mathematics Letters, 2(4):321 – 325, 1989.\n"
                << "c =     [ 0     1/2  3/4   1  ]'\n"
                << "A =     [ 0                   ]\n"
                << "        [ 1/2    0            ]\n"
                << "        [  0    3/4   0       ]\n"
                << "        [ 2/9   1/3  4/9   0  ]\n"
                << "b     = [ 2/9   1/3  4/9   0  ]'\n"
                << "bstar = [ 7/24  1/4  1/3  1/8 ]'";
    return Description.str();
  }

  void setUseFSAL(bool a) { this->useFSAL_ = a; this->isInitialized_ = false; }

protected:

  void setupTableau()
  {
    typedef Teuchos::ScalarTraits<Scalar> ST;
    using Teuchos::as;
    int NumStages = 4;
    Teuchos::SerialDenseMatrix<int,Scalar> A(NumStages,NumStages);
    Teuchos::SerialDenseVector<int,Scalar> b(NumStages);
    Teuchos::SerialDenseVector<int,Scalar> c(NumStages);
    Teuchos::SerialDenseVector<int,Scalar> bstar(NumStages);

    const Scalar one = ST::one();
    const Scalar zero = ST::zero();
    const Scalar onehalf = one/(2*one);
    const Scalar onethird = one/(3*one);
    const Scalar threefourths = (3*one)/(4*one);
    const Scalar twoninths = (2*one)/(9*one);
    const Scalar fourninths = (4*one)/(9*one);

    // Fill A:
    A(0,0) =     zero; A(0,1) =        zero; A(0,2) =      zero; A(0,3) = zero;
    A(1,0) =  onehalf; A(1,1) =        zero; A(1,2) =      zero; A(1,3) = zero;
    A(2,0) =     zero; A(2,1) =threefourths; A(2,2) =      zero; A(2,3) = zero;
    A(3,0) =twoninths; A(3,1) =    onethird; A(3,2) =fourninths; A(3,3) = zero;

    // Fill b:
    b(0) = A(3,0); b(1) = A(3,1); b(2) = A(3,2); b(3) = A(3,3);

    // Fill c:
    c(0) = zero; c(1) = onehalf; c(2) = threefourths; c(3) = one;

    // Fill bstar
    bstar(0) = as<Scalar>(7*one/(24*one));
    bstar(1) = as<Scalar>(1*one/(4*one));
    bstar(2) = as<Scalar>(1*one/(3*one));
    bstar(3) = as<Scalar>(1*one/(8*one));
    int order = 3;

    this->tableau_ = Teuchos::rcp(new RKButcherTableau<Scalar>(
      this->getStepperType(),A,b,c,order,order,order,bstar));
  }
};


// ------------------------------------------------------------------------
template<class Scalar>
Teuchos::RCP<StepperERK_BogackiShampine32<Scalar> >
createStepperERK_BogackiShampine32(
  const Teuchos::RCP<const Thyra::ModelEvaluator<Scalar> >& model,
  Teuchos::RCP<Teuchos::ParameterList> pl)
{
  auto stepper = Teuchos::rcp(new StepperERK_BogackiShampine32<Scalar>());
  stepper->setStepperRKValues(pl);

  if (model != Teuchos::null) {
    stepper->setModel(model);
    stepper->initialize();
  }

  return stepper;
}


// ----------------------------------------------------------------------------
template<class Scalar>
class StepperERK_Merson45 :
  virtual public StepperExplicitRK<Scalar>
{
public:
  StepperERK_Merson45()
  {
    this->setStepperName("Merson 4(5) Pair");
    this->setStepperType("Merson 4(5) Pair");
    this->setupTableau();
    this->setupDefault();
    this->setUseFSAL(            false);
    this->setICConsistency(      "None");
    this->setICConsistencyCheck( false);
  }

  StepperERK_Merson45(
    const Teuchos::RCP<const Thyra::ModelEvaluator<Scalar> >& appModel,
    bool useFSAL,
    std::string ICConsistency,
    bool ICConsistencyCheck,
    bool useEmbedded,
    const Teuchos::RCP<StepperRKAppAction<Scalar> >& stepperRKAppAction)
  {
    this->setStepperName("Merson 4(5) Pair");
    this->setStepperType("Merson 4(5) Pair");
    this->setupTableau();
    this->setup(appModel, useFSAL, ICConsistency,
                ICConsistencyCheck, useEmbedded, stepperRKAppAction);
  }

  std::string getDescription() const
  {
    std::ostringstream Description;
    Description << this->getStepperType() << "\n"
                << "Solving Ordinary Differential Equations I:\n"
                << "Nonstiff Problems, 2nd Revised Edition\n"
                << "E. Hairer, S.P. Norsett, G. Wanner\n"
                << "Table 4.1, pg 167\n"
                << "c =     [  0    1/3  1/3  1/2   1  ]'\n"
                << "A =     [  0                       ]\n"
                << "        [ 1/3    0                 ]\n"
                << "        [ 1/6   1/6   0            ]\n"
                << "        [ 1/8    0   3/8   0       ]\n"
                << "        [ 1/2    0  -3/2   2    0  ]\n"
                << "b     = [ 1/6    0    0   2/3  1/6 ]'\n"
                << "bstar = [ 1/10   0  3/10  2/5  1/5 ]'";
    return Description.str();
  }


protected:

  void setupTableau()
  {
    typedef Teuchos::ScalarTraits<Scalar> ST;
    using Teuchos::as;
    int NumStages = 5;
    Teuchos::SerialDenseMatrix<int,Scalar> A(NumStages,NumStages, true);
    Teuchos::SerialDenseVector<int,Scalar> b(NumStages, true);
    Teuchos::SerialDenseVector<int,Scalar> c(NumStages, true);
    Teuchos::SerialDenseVector<int,Scalar> bstar(NumStages, true);

    const Scalar one = ST::one();
    const Scalar zero = ST::zero();

    // Fill A:
    A(1,0) = as<Scalar>(one/(3*one));;

    A(2,0) = as<Scalar>(one/(6*one));;
    A(2,1) = as<Scalar>(one/(6*one));;

    A(3,0) = as<Scalar>(one/(8*one));;
    A(3,2) = as<Scalar>(3*one/(8*one));;

    A(4,0) = as<Scalar>(one/(2*one));;
    A(4,2) = as<Scalar>(-3*one/(2*one));;
    A(4,3) = 2*one;

    // Fill b:
    b(0) = as<Scalar>(one/(6*one));
    b(3) = as<Scalar>(2*one/(3*one));
    b(4) = as<Scalar>(one/(6*one));

    // Fill c:
    c(0) = zero;
    c(1) = as<Scalar>(1*one/(3*one));
    c(2) = as<Scalar>(1*one/(3*one));
    c(3) = as<Scalar>(1*one/(2*one));
    c(4) = one;

    // Fill bstar
    bstar(0) = as<Scalar>(1*one/(10*one));
    bstar(2) = as<Scalar>(3*one/(10*one));
    bstar(3) = as<Scalar>(2*one/(5*one));
    bstar(4) = as<Scalar>(1*one/(5*one));
    int order = 4;

    this->tableau_ = Teuchos::rcp(new RKButcherTableau<Scalar>(
      this->getStepperType(),A,b,c,order,order,order,bstar));
  }
};


// ------------------------------------------------------------------------
template<class Scalar>
Teuchos::RCP<StepperERK_Merson45<Scalar> >
createStepperERK_Merson45(
  const Teuchos::RCP<const Thyra::ModelEvaluator<Scalar> >& model,
  Teuchos::RCP<Teuchos::ParameterList> pl)
{
  auto stepper = Teuchos::rcp(new StepperERK_Merson45<Scalar>());
  stepper->setStepperRKValues(pl);

  if (model != Teuchos::null) {
    stepper->setModel(model);
    stepper->initialize();
  }

  return stepper;
}


// ----------------------------------------------------------------------------
template<class Scalar>
class StepperERK_3_8Rule :
  virtual public StepperExplicitRK<Scalar>
{
public:

  StepperERK_3_8Rule()
  {
    this->setStepperName("RK Explicit 3/8 Rule");
    this->setStepperType("RK Explicit 3/8 Rule");
    this->setupTableau();
    this->setupDefault();
    this->setUseFSAL(            false);
    this->setICConsistency(      "None");
    this->setICConsistencyCheck( false);
  }

  StepperERK_3_8Rule(
    const Teuchos::RCP<const Thyra::ModelEvaluator<Scalar> >& appModel,
    bool useFSAL,
    std::string ICConsistency,
    bool ICConsistencyCheck,
    bool useEmbedded,
    const Teuchos::RCP<StepperRKAppAction<Scalar> >& stepperRKAppAction)
  {
    this->setStepperName("RK Explicit 3/8 Rule");
    this->setStepperType("RK Explicit 3/8 Rule");
    this->setupTableau();
    this->setup(appModel, useFSAL, ICConsistency,
                ICConsistencyCheck, useEmbedded, stepperRKAppAction);
  }

  std::string getDescription() const
  {
    std::ostringstream Description;
    Description << this->getStepperType() << "\n"
                << "Solving Ordinary Differential Equations I:\n"
                << "Nonstiff Problems, 2nd Revised Edition\n"
                << "E. Hairer, S.P. Norsett, G. Wanner\n"
                << "Table 1.2, pg 138\n"
                << "c = [  0  1/3 2/3  1  ]'\n"
                << "A = [  0              ]\n"
                << "    [ 1/3  0          ]\n"
                << "    [-1/3  1   0      ]\n"
                << "    [  1  -1   1   0  ]\n"
                << "b = [ 1/8 3/8 3/8 1/8 ]'";
    return Description.str();
  }


protected:

  void setupTableau()
  {
    typedef Teuchos::ScalarTraits<Scalar> ST;
    using Teuchos::as;
    int NumStages = 4;
    Teuchos::SerialDenseMatrix<int,Scalar> A(NumStages,NumStages);
    Teuchos::SerialDenseVector<int,Scalar> b(NumStages);
    Teuchos::SerialDenseVector<int,Scalar> c(NumStages);

    const Scalar one = ST::one();
    const Scalar zero = ST::zero();
    const Scalar onethird     = as<Scalar>(one/(3*one));
    const Scalar twothirds    = as<Scalar>(2*one/(3*one));
    const Scalar oneeighth    = as<Scalar>(one/(8*one));
    const Scalar threeeighths = as<Scalar>(3*one/(8*one));

    // Fill A:
    A(0,0) =      zero; A(0,1) = zero; A(0,2) = zero; A(0,3) = zero;
    A(1,0) =  onethird; A(1,1) = zero; A(1,2) = zero; A(1,3) = zero;
    A(2,0) = -onethird; A(2,1) =  one; A(2,2) = zero; A(2,3) = zero;
    A(3,0) =       one; A(3,1) = -one; A(3,2) =  one; A(3,3) = zero;

    // Fill b:
    b(0) =oneeighth; b(1) =threeeighths; b(2) =threeeighths; b(3) =oneeighth;

    // Fill c:
    c(0) = zero; c(1) = onethird; c(2) = twothirds; c(3) = one;

    int order = 4;

    this->tableau_ = Teuchos::rcp(new RKButcherTableau<Scalar>(
      this->getStepperType(),A,b,c,order,order,order));
  }
};


// ------------------------------------------------------------------------
template<class Scalar>
Teuchos::RCP<StepperERK_3_8Rule<Scalar> >
createStepperERK_3_8Rule(
  const Teuchos::RCP<const Thyra::ModelEvaluator<Scalar> >& model,
  Teuchos::RCP<Teuchos::ParameterList> pl)
{
  auto stepper = Teuchos::rcp(new StepperERK_3_8Rule<Scalar>());
  stepper->setStepperRKValues(pl);

  if (model != Teuchos::null) {
    stepper->setModel(model);
    stepper->initialize();
  }

  return stepper;
}


// ----------------------------------------------------------------------------
template<class Scalar>
class StepperERK_4Stage3rdOrderRunge :
  virtual public StepperExplicitRK<Scalar>
{
public:
  StepperERK_4Stage3rdOrderRunge()
  {
    this->setStepperName("RK Explicit 4 Stage 3rd order by Runge");
    this->setStepperType("RK Explicit 4 Stage 3rd order by Runge");
    this->setupTableau();
    this->setupDefault();
    this->setUseFSAL(            false);
    this->setICConsistency(      "None");
    this->setICConsistencyCheck( false);
  }

  StepperERK_4Stage3rdOrderRunge(
    const Teuchos::RCP<const Thyra::ModelEvaluator<Scalar> >& appModel,
    bool useFSAL,
    std::string ICConsistency,
    bool ICConsistencyCheck,
    bool useEmbedded,
    const Teuchos::RCP<StepperRKAppAction<Scalar> >& stepperRKAppAction)
  {
    this->setStepperName("RK Explicit 4 Stage 3rd order by Runge");
    this->setStepperType("RK Explicit 4 Stage 3rd order by Runge");
    this->setupTableau();
    this->setup(appModel, useFSAL, ICConsistency,
                ICConsistencyCheck, useEmbedded, stepperRKAppAction);
  }

  std::string getDescription() const
  {
    std::ostringstream Description;
    Description << this->getStepperType() << "\n"
                << "Solving Ordinary Differential Equations I:\n"
                << "Nonstiff Problems, 2nd Revised Edition\n"
                << "E. Hairer, S.P. Norsett, G. Wanner\n"
                << "Table 1.1, pg 135\n"
                << "c = [  0  1/2  1   1  ]'\n"
                << "A = [  0              ]\n"
                << "    [ 1/2  0          ]\n"
                << "    [  0   1   0      ]\n"
                << "    [  0   0   1   0  ]\n"
                << "b = [ 1/6 2/3  0  1/6 ]'";
    return Description.str();
  }
protected:

  void setupTableau()
  {
    typedef Teuchos::ScalarTraits<Scalar> ST;
    int NumStages = 4;
    Teuchos::SerialDenseMatrix<int,Scalar> A(NumStages,NumStages);
    Teuchos::SerialDenseVector<int,Scalar> b(NumStages);
    Teuchos::SerialDenseVector<int,Scalar> c(NumStages);

    const Scalar one = ST::one();
    const Scalar onehalf = one/(2*one);
    const Scalar onesixth = one/(6*one);
    const Scalar twothirds = 2*one/(3*one);
    const Scalar zero = ST::zero();

    // Fill A:
    A(0,0) =    zero; A(0,1) = zero; A(0,2) = zero; A(0,3) = zero;
    A(1,0) = onehalf; A(1,1) = zero; A(1,2) = zero; A(1,3) = zero;
    A(2,0) =    zero; A(2,1) =  one; A(2,2) = zero; A(2,3) = zero;
    A(3,0) =    zero; A(3,1) = zero; A(3,2) =  one; A(3,3) = zero;

    // Fill b:
    b(0) = onesixth; b(1) = twothirds; b(2) = zero; b(3) = onesixth;

    // Fill c:
    c(0) = zero; c(1) = onehalf; c(2) = one; c(3) = one;

    int order = 3;

    this->tableau_ = Teuchos::rcp(new RKButcherTableau<Scalar>(
      this->getStepperType(),A,b,c,order,order,order));
  }
};


// ------------------------------------------------------------------------
template<class Scalar>
Teuchos::RCP<StepperERK_4Stage3rdOrderRunge<Scalar> >
createStepperERK_4Stage3rdOrderRunge(
  const Teuchos::RCP<const Thyra::ModelEvaluator<Scalar> >& model,
  Teuchos::RCP<Teuchos::ParameterList> pl)
{
  auto stepper = Teuchos::rcp(new StepperERK_4Stage3rdOrderRunge<Scalar>());
  stepper->setStepperRKValues(pl);

  if (model != Teuchos::null) {
    stepper->setModel(model);
    stepper->initialize();
  }

  return stepper;
}


// ----------------------------------------------------------------------------
template<class Scalar>
class StepperERK_5Stage3rdOrderKandG :
  virtual public StepperExplicitRK<Scalar>
{
public:
  StepperERK_5Stage3rdOrderKandG()
  {
    this->setStepperName("RK Explicit 5 Stage 3rd order by Kinnmark and Gray");
    this->setStepperType("RK Explicit 5 Stage 3rd order by Kinnmark and Gray");
    this->setupTableau();
    this->setupDefault();
    this->setUseFSAL(            false);
    this->setICConsistency(      "None");
    this->setICConsistencyCheck( false);
  }

  StepperERK_5Stage3rdOrderKandG(
    const Teuchos::RCP<const Thyra::ModelEvaluator<Scalar> >& appModel,
    bool useFSAL,
    std::string ICConsistency,
    bool ICConsistencyCheck,
    bool useEmbedded,
    const Teuchos::RCP<StepperRKAppAction<Scalar> >& stepperRKAppAction)
  {
    this->setStepperName("RK Explicit 5 Stage 3rd order by Kinnmark and Gray");
    this->setStepperType("RK Explicit 5 Stage 3rd order by Kinnmark and Gray");
    this->setupTableau();
    this->setup(appModel, useFSAL, ICConsistency,
                ICConsistencyCheck, useEmbedded,stepperRKAppAction);
  }

  std::string getDescription() const
  {
    std::ostringstream Description;
    Description << this->getStepperType() << "\n"
                << "Kinnmark & Gray 5 stage, 3rd order scheme \n"
                << "Modified by P. Ullrich.  From the prim_advance_mod.F90 \n"
                << "routine in the HOMME atmosphere model code.\n"
                << "c = [  0  1/5  1/5  1/3  2/3  ]'\n"
                << "A = [  0                      ]\n"
                << "    [ 1/5  0                  ]\n"
                << "    [  0  1/5   0             ]\n"
                << "    [  0   0   1/3   0        ]\n"
                << "    [  0   0    0   2/3   0   ]\n"
                << "b = [ 1/4  0    0    0   3/4  ]'";
    return Description.str();
  }

protected:

  void setupTableau()
  {
    typedef Teuchos::ScalarTraits<Scalar> ST;
    int NumStages = 5;
    Teuchos::SerialDenseMatrix<int,Scalar> A(NumStages,NumStages);
    Teuchos::SerialDenseVector<int,Scalar> b(NumStages);
    Teuchos::SerialDenseVector<int,Scalar> c(NumStages);

    const Scalar one = ST::one();
    const Scalar onefifth = one/(5*one);
    const Scalar onefourth = one/(4*one);
    const Scalar onethird = one/(3*one);
    const Scalar twothirds = 2*one/(3*one);
    const Scalar threefourths = 3*one/(4*one);
    const Scalar zero = ST::zero();

    // Fill A:
    A(0,0) =     zero; A(0,1) =     zero; A(0,2) =     zero; A(0,3) =      zero; A(0,4) = zero;
    A(1,0) = onefifth; A(1,1) =     zero; A(1,2) =     zero; A(1,3) =      zero; A(1,4) = zero;
    A(2,0) =     zero; A(2,1) = onefifth; A(2,2) =     zero; A(2,3) =      zero; A(2,4) = zero;
    A(3,0) =     zero; A(3,1) =     zero; A(3,2) = onethird; A(3,3) =      zero; A(3,4) = zero;
    A(4,0) =     zero; A(4,1) =     zero; A(4,2) =     zero; A(4,3) = twothirds; A(4,4) = zero;

    // Fill b:
    b(0) =onefourth; b(1) =zero; b(2) =zero; b(3) =zero; b(4) =threefourths;

    // Fill c:
    c(0) =zero; c(1) =onefifth; c(2) =onefifth; c(3) =onethird; c(4) =twothirds;

    int order = 3;

    this->tableau_ = Teuchos::rcp(new RKButcherTableau<Scalar>(
      this->getStepperType(),A,b,c,order,order,order));
  }
};


// ------------------------------------------------------------------------
template<class Scalar>
Teuchos::RCP<StepperERK_5Stage3rdOrderKandG<Scalar> >
createStepperERK_5Stage3rdOrderKandG(
  const Teuchos::RCP<const Thyra::ModelEvaluator<Scalar> >& model,
  Teuchos::RCP<Teuchos::ParameterList> pl)
{
  auto stepper = Teuchos::rcp(new StepperERK_5Stage3rdOrderKandG<Scalar>());
  stepper->setStepperRKValues(pl);

  if (model != Teuchos::null) {
    stepper->setModel(model);
    stepper->initialize();
  }

  return stepper;
}


// ----------------------------------------------------------------------------
template<class Scalar>
class StepperERK_3Stage3rdOrder :
  virtual public StepperExplicitRK<Scalar>
{
public:
  StepperERK_3Stage3rdOrder()
  {
    this->setStepperName("RK Explicit 3 Stage 3rd order");
    this->setStepperType("RK Explicit 3 Stage 3rd order");
    this->setupTableau();
    this->setupDefault();
    this->setUseFSAL(            false);
    this->setICConsistency(      "None");
    this->setICConsistencyCheck( false);
  }

  StepperERK_3Stage3rdOrder(
    const Teuchos::RCP<const Thyra::ModelEvaluator<Scalar> >& appModel,
    bool useFSAL,
    std::string ICConsistency,
    bool ICConsistencyCheck,
    bool useEmbedded,
    const Teuchos::RCP<StepperRKAppAction<Scalar> >& stepperRKAppAction)
  {
    this->setStepperName("RK Explicit 3 Stage 3rd order");
    this->setStepperType("RK Explicit 3 Stage 3rd order");
    this->setupTableau();
    this->setup(appModel, useFSAL, ICConsistency,
                ICConsistencyCheck, useEmbedded, stepperRKAppAction);
  }

  std::string getDescription() const
  {
    std::ostringstream Description;
    Description << this->getStepperType() << "\n"
                << "c = [  0  1/2  1  ]'\n"
                << "A = [  0          ]\n"
                << "    [ 1/2  0      ]\n"
                << "    [ -1   2   0  ]\n"
                << "b = [ 1/6 4/6 1/6 ]'";
    return Description.str();
  }

protected:

  void setupTableau()
  {
    typedef Teuchos::ScalarTraits<Scalar> ST;
    const Scalar one = ST::one();
    const Scalar two = Teuchos::as<Scalar>(2*one);
    const Scalar zero = ST::zero();
    const Scalar onehalf = one/(2*one);
    const Scalar onesixth = one/(6*one);
    const Scalar foursixth = 4*one/(6*one);

    int NumStages = 3;
    Teuchos::SerialDenseMatrix<int,Scalar> A(NumStages,NumStages);
    Teuchos::SerialDenseVector<int,Scalar> b(NumStages);
    Teuchos::SerialDenseVector<int,Scalar> c(NumStages);

    // Fill A:
    A(0,0) =    zero; A(0,1) = zero; A(0,2) = zero;
    A(1,0) = onehalf; A(1,1) = zero; A(1,2) = zero;
    A(2,0) =    -one; A(2,1) =  two; A(2,2) = zero;

    // Fill b:
    b(0) = onesixth; b(1) = foursixth; b(2) = onesixth;

    // fill c:
    c(0) = zero; c(1) = onehalf; c(2) = one;

    int order = 3;

    this->tableau_ = Teuchos::rcp(new RKButcherTableau<Scalar>(
      this->getStepperType(),A,b,c,order,order,order));
  }
};


// ------------------------------------------------------------------------
template<class Scalar>
Teuchos::RCP<StepperERK_3Stage3rdOrder<Scalar> >
createStepperERK_3Stage3rdOrder(
  const Teuchos::RCP<const Thyra::ModelEvaluator<Scalar> >& model,
  Teuchos::RCP<Teuchos::ParameterList> pl)
{
  auto stepper = Teuchos::rcp(new StepperERK_3Stage3rdOrder<Scalar>());
  stepper->setStepperRKValues(pl);

  if (model != Teuchos::null) {
    stepper->setModel(model);
    stepper->initialize();
  }

  return stepper;
}


// ----------------------------------------------------------------------------
template<class Scalar>
class StepperERK_3Stage3rdOrderTVD :
  virtual public StepperExplicitRK<Scalar>
{
public:
  StepperERK_3Stage3rdOrderTVD()
  {
    this->setStepperName("RK Explicit 3 Stage 3rd order TVD");
    this->setStepperType("RK Explicit 3 Stage 3rd order TVD");
    this->setupTableau();
    this->setupDefault();
    this->setUseFSAL(            false);
    this->setICConsistency(      "None");
    this->setICConsistencyCheck( false);
  }

  StepperERK_3Stage3rdOrderTVD(
    const Teuchos::RCP<const Thyra::ModelEvaluator<Scalar> >& appModel,
    bool useFSAL,
    std::string ICConsistency,
    bool ICConsistencyCheck,
    bool useEmbedded,
    const Teuchos::RCP<StepperRKAppAction<Scalar> >& stepperRKAppAction)
  {
    this->setStepperName("RK Explicit 3 Stage 3rd order TVD");
    this->setStepperType("RK Explicit 3 Stage 3rd order TVD");
    this->setupTableau();
    this->setup(appModel, useFSAL, ICConsistency,
                ICConsistencyCheck, useEmbedded, stepperRKAppAction);
  }

  std::string getDescription() const
  {
    std::ostringstream Description;
    Description << this->getStepperType() << "\n"
                  << "This Stepper is known as 'RK Explicit 3 Stage 3rd order TVD' or 'SSPERK33'.\n"
                  << "Sigal Gottlieb and Chi-Wang Shu\n"
                  << "`Total Variation Diminishing Runge-Kutta Schemes'\n"
                  << "Mathematics of Computation\n"
                  << "Volume 67, Number 221, January 1998, pp. 73-85\n"
                  << "c = [  0   1  1/2 ]'\n"
                  << "A = [  0          ]\n"
                  << "    [  1   0      ]\n"
                  << "    [ 1/4 1/4  0  ]\n"
                  << "b = [ 1/6 1/6 4/6 ]'\n"
                  << "This is also written in the following set of updates.\n"
                  << "u1 = u^n + dt L(u^n)\n"
                  << "u2 = 3 u^n/4 + u1/4 + dt L(u1)/4\n"
                  << "u^(n+1) = u^n/3 + 2 u2/2 + 2 dt L(u2)/3";
    return Description.str();
  }

protected:

  void setupTableau()
  {
    typedef Teuchos::ScalarTraits<Scalar> ST;
    using Teuchos::as;
    const Scalar one = ST::one();
    const Scalar zero = ST::zero();
    const Scalar onehalf = one/(2*one);
    const Scalar onefourth = one/(4*one);
    const Scalar onesixth = one/(6*one);
    const Scalar foursixth = 4*one/(6*one);

    int NumStages = 3;
    Teuchos::SerialDenseMatrix<int,Scalar> A(NumStages,NumStages);
    Teuchos::SerialDenseVector<int,Scalar> b(NumStages);
    Teuchos::SerialDenseVector<int,Scalar> c(NumStages);
    Teuchos::SerialDenseVector<int,Scalar> bstar(NumStages);

    // Fill A:
    A(0,0) =      zero; A(0,1) =      zero; A(0,2) = zero;
    A(1,0) =       one; A(1,1) =      zero; A(1,2) = zero;
    A(2,0) = onefourth; A(2,1) = onefourth; A(2,2) = zero;

    // Fill b:
    b(0) = onesixth; b(1) = onesixth; b(2) = foursixth;

    // fill c:
    c(0) = zero; c(1) = one; c(2) = onehalf;

    // Fill bstar:
    bstar(0) = as<Scalar>(0.291485418878409);
    bstar(1) = as<Scalar>(0.291485418878409);
    bstar(2) = as<Scalar>(0.417029162243181);

    int order = 3;

    this->tableau_ = Teuchos::rcp(new RKButcherTableau<Scalar>(
      this->getStepperType(),A,b,c,order,order,order,bstar));
    this->tableau_->setTVD(true);
    this->tableau_->setTVDCoeff(1.0);
  }
};


// ------------------------------------------------------------------------
template<class Scalar>
Teuchos::RCP<StepperERK_3Stage3rdOrderTVD<Scalar> >
createStepperERK_3Stage3rdOrderTVD(
  const Teuchos::RCP<const Thyra::ModelEvaluator<Scalar> >& model,
  Teuchos::RCP<Teuchos::ParameterList> pl)
{
  auto stepper = Teuchos::rcp(new StepperERK_3Stage3rdOrderTVD<Scalar>());

  // Test for aliases.
  if (pl != Teuchos::null) {
    auto stepperType =
      pl->get<std::string>("Stepper Type", stepper->getStepperType());

    TEUCHOS_TEST_FOR_EXCEPTION(
      stepperType != stepper->getStepperType() &&
      stepperType != "SSPERK33" && stepperType != "SSPRK3", std::logic_error,
      "  ParameterList 'Stepper Type' (='" + stepperType +"')\n"
      "  does not match type for this Stepper (='" + stepper->getStepperType() +
      "')\n  or one of its aliases ('SSPERK33' or 'SSPRK3').\n");

    // Reset default StepperType.
    pl->set<std::string>("Stepper Type", stepper->getStepperType());
  }

  stepper->setStepperRKValues(pl);

  if (model != Teuchos::null) {
    stepper->setModel(model);
    stepper->initialize();
  }

  return stepper;
}


// ----------------------------------------------------------------------------
template<class Scalar>
class StepperERK_3Stage3rdOrderHeun :
  virtual public StepperExplicitRK<Scalar>
{
public:
  StepperERK_3Stage3rdOrderHeun()
  {
    this->setStepperName("RK Explicit 3 Stage 3rd order by Heun");
    this->setStepperType("RK Explicit 3 Stage 3rd order by Heun");
    this->setupTableau();
    this->setupDefault();
    this->setUseFSAL(            false);
    this->setICConsistency(      "None");
    this->setICConsistencyCheck( false);
  }

  StepperERK_3Stage3rdOrderHeun(
    const Teuchos::RCP<const Thyra::ModelEvaluator<Scalar> >& appModel,
    bool useFSAL,
    std::string ICConsistency,
    bool ICConsistencyCheck,
    bool useEmbedded,
    const Teuchos::RCP<StepperRKAppAction<Scalar> >& stepperRKAppAction)
  {
    this->setStepperName("RK Explicit 3 Stage 3rd order by Heun");
    this->setStepperType("RK Explicit 3 Stage 3rd order by Heun");
    this->setupTableau();
    this->setup(appModel, useFSAL, ICConsistency,
                ICConsistencyCheck, useEmbedded, stepperRKAppAction);
  }

  std::string getDescription() const
  {
    std::ostringstream Description;
    Description << this->getStepperType() << "\n"
                << "Solving Ordinary Differential Equations I:\n"
                << "Nonstiff Problems, 2nd Revised Edition\n"
                << "E. Hairer, S.P. Norsett, G. Wanner\n"
                << "Table 1.1, pg 135\n"
                << "c = [  0  1/3 2/3 ]'\n"
                << "A = [  0          ] \n"
                << "    [ 1/3  0      ]\n"
                << "    [  0  2/3  0  ]\n"
                << "b = [ 1/4  0  3/4 ]'";
    return Description.str();
  }

protected:

  void setupTableau()
  {
    typedef Teuchos::ScalarTraits<Scalar> ST;
    const Scalar one = ST::one();
    const Scalar zero = ST::zero();
    const Scalar onethird = one/(3*one);
    const Scalar twothirds = 2*one/(3*one);
    const Scalar onefourth = one/(4*one);
    const Scalar threefourths = 3*one/(4*one);

    int NumStages = 3;
    Teuchos::SerialDenseMatrix<int,Scalar> A(NumStages,NumStages);
    Teuchos::SerialDenseVector<int,Scalar> b(NumStages);
    Teuchos::SerialDenseVector<int,Scalar> c(NumStages);

    // Fill A:
    A(0,0) =     zero; A(0,1) =      zero; A(0,2) = zero;
    A(1,0) = onethird; A(1,1) =      zero; A(1,2) = zero;
    A(2,0) =     zero; A(2,1) = twothirds; A(2,2) = zero;

    // Fill b:
    b(0) = onefourth; b(1) = zero; b(2) = threefourths;

    // fill c:
    c(0) = zero; c(1) = onethird; c(2) = twothirds;

    int order = 3;

    this->tableau_ = Teuchos::rcp(new RKButcherTableau<Scalar>(
      this->getStepperType(),A,b,c,order,order,order));
  }
};


// ------------------------------------------------------------------------
template<class Scalar>
Teuchos::RCP<StepperERK_3Stage3rdOrderHeun<Scalar> >
createStepperERK_3Stage3rdOrderHeun(
  const Teuchos::RCP<const Thyra::ModelEvaluator<Scalar> >& model,
  Teuchos::RCP<Teuchos::ParameterList> pl)
{
  auto stepper = Teuchos::rcp(new StepperERK_3Stage3rdOrderHeun<Scalar>());
  stepper->setStepperRKValues(pl);

  if (model != Teuchos::null) {
    stepper->setModel(model);
    stepper->initialize();
  }

  return stepper;
}


// ----------------------------------------------------------------------------
template<class Scalar>
class StepperERK_Midpoint :
  virtual public StepperExplicitRK<Scalar>
{
public:
  StepperERK_Midpoint()
  {
    this->setStepperName("RK Explicit Midpoint");
    this->setStepperType("RK Explicit Midpoint");
    this->setupTableau();
    this->setupDefault();
    this->setUseFSAL(            false);
    this->setICConsistency(      "None");
    this->setICConsistencyCheck( false);
  }

  StepperERK_Midpoint(
    const Teuchos::RCP<const Thyra::ModelEvaluator<Scalar> >& appModel,
    bool useFSAL,
    std::string ICConsistency,
    bool ICConsistencyCheck,
    bool useEmbedded,
    const Teuchos::RCP<StepperRKAppAction<Scalar> >& stepperRKAppAction)
  {
    this->setStepperName("RK Explicit Midpoint");
    this->setStepperType("RK Explicit Midpoint");
    this->setupTableau();
    this->setup(appModel, useFSAL, ICConsistency,
                ICConsistencyCheck, useEmbedded, stepperRKAppAction);
  }

  std::string getDescription() const
  {
    std::ostringstream Description;
    Description << this->getStepperType() << "\n"
                << "Solving Ordinary Differential Equations I:\n"
                << "Nonstiff Problems, 2nd Revised Edition\n"
                << "E. Hairer, S.P. Norsett, G. Wanner\n"
                << "Table 1.1, pg 135\n"
                << "c = [  0  1/2 ]'\n"
                << "A = [  0      ]\n"
                << "    [ 1/2  0  ]\n"
                << "b = [  0   1  ]'";
    return Description.str();
  }

protected:

  void setupTableau()
  {
    typedef Teuchos::ScalarTraits<Scalar> ST;
    const Scalar one = ST::one();
    const Scalar zero = ST::zero();
    const Scalar onehalf = one/(2*one);

    int NumStages = 2;
    Teuchos::SerialDenseMatrix<int,Scalar> A(NumStages,NumStages);
    Teuchos::SerialDenseVector<int,Scalar> b(NumStages);
    Teuchos::SerialDenseVector<int,Scalar> c(NumStages);

    // Fill A:
    A(0,0) =    zero; A(0,1) = zero;
    A(1,0) = onehalf; A(1,1) = zero;

    // Fill b:
    b(0) = zero; b(1) = one;

    // fill c:
    c(0) = zero; c(1) = onehalf;

    int order = 2;

    this->tableau_ = Teuchos::rcp(new RKButcherTableau<Scalar>(
      this->getStepperType(),A,b,c,order,order,order));
  }
};


// ------------------------------------------------------------------------
template<class Scalar>
Teuchos::RCP<StepperERK_Midpoint<Scalar> >
createStepperERK_Midpoint(
  const Teuchos::RCP<const Thyra::ModelEvaluator<Scalar> >& model,
  Teuchos::RCP<Teuchos::ParameterList> pl)
{
  auto stepper = Teuchos::rcp(new StepperERK_Midpoint<Scalar>());
  stepper->setStepperRKValues(pl);

  if (model != Teuchos::null) {
    stepper->setModel(model);
    stepper->initialize();
  }

  return stepper;
}


// ----------------------------------------------------------------------------
template<class Scalar>
class StepperERK_Ralston :
  virtual public StepperExplicitRK<Scalar>
{
public:
  StepperERK_Ralston()
  {
    this->setStepperName("RK Explicit Ralston");
    this->setStepperType("RK Explicit Ralston");
    this->setupTableau();
    this->setupDefault();
    this->setUseFSAL(            false);
    this->setICConsistency(      "None");
    this->setICConsistencyCheck( false);
  }

  StepperERK_Ralston(
    const Teuchos::RCP<const Thyra::ModelEvaluator<Scalar> >& appModel,
    bool useFSAL,
    std::string ICConsistency,
    bool ICConsistencyCheck,
    bool useEmbedded,
    const Teuchos::RCP<StepperRKAppAction<Scalar> >& stepperRKAppAction)
  {
    this->setStepperName("RK Explicit Ralston");
    this->setStepperType("RK Explicit Ralston");
    this->setupTableau();
    this->setup(appModel, useFSAL, ICConsistency,
                ICConsistencyCheck, useEmbedded, stepperRKAppAction);
  }

  std::string getDescription() const
  {
    std::ostringstream Description;
    Description << this->getStepperType() << "\n"
                << "This Stepper is known as 'RK Explicit Ralston' or 'RK2'.\n"
                << "c = [   0   2/3 ]'\n"
                << "A = [   0       ]\n"
                << "    [  2/3   0  ]\n"
                << "b = [  1/4  3/4 ]'";
    return Description.str();
  }

protected:

  void setupTableau()
  {
   typedef Teuchos::ScalarTraits<Scalar> ST;
    const Scalar one = ST::one();
    const Scalar zero = ST::zero();

    const int NumStages = 2;
    const int order = 2;
    Teuchos::SerialDenseMatrix<int,Scalar> A(NumStages,NumStages);
    Teuchos::SerialDenseVector<int,Scalar> b(NumStages);
    Teuchos::SerialDenseVector<int,Scalar> c(NumStages);

    // Fill A:
    A(0,0) = zero; A(0,1) = zero; A(1,1) = zero;
    A(1,0) = (2*one)/(3*one);

    // Fill b:
    b(0) = (1*one)/(4*one);
    b(1) = (3*one)/(4*one);

    // fill c:
    c(0) = zero;
    c(1) = (2*one)/(3*one);


    this->tableau_ = Teuchos::rcp(new RKButcherTableau<Scalar>(
      this->getStepperType(),A,b,c,order,order,order));
    this->tableau_->setTVD(true);
    this->tableau_->setTVDCoeff(0.5);
  }
};


// ------------------------------------------------------------------------
template<class Scalar>
Teuchos::RCP<StepperERK_Ralston<Scalar> >
createStepperERK_Ralston(
  const Teuchos::RCP<const Thyra::ModelEvaluator<Scalar> >& model,
  Teuchos::RCP<Teuchos::ParameterList> pl)
{
  auto stepper = Teuchos::rcp(new StepperERK_Ralston<Scalar>());

  // Test for aliases.
  if (pl != Teuchos::null) {
    auto stepperType =
      pl->get<std::string>("Stepper Type", stepper->getStepperType());

    TEUCHOS_TEST_FOR_EXCEPTION(
      stepperType != stepper->getStepperType() &&
      stepperType != "RK2", std::logic_error,
      "  ParameterList 'Stepper Type' (='" + stepperType +"')\n"
      "  does not match type for this Stepper (='" + stepper->getStepperType() +
      "')\n  or one of its aliases ('RK2').\n");

    // Reset default StepperType.
    pl->set<std::string>("Stepper Type", stepper->getStepperType());
  }

  stepper->setStepperRKValues(pl);

  if (model != Teuchos::null) {
    stepper->setModel(model);
    stepper->initialize();
  }

  return stepper;
}


// ----------------------------------------------------------------------------
template<class Scalar>
class StepperERK_Trapezoidal :
  virtual public StepperExplicitRK<Scalar>
{
public:
  StepperERK_Trapezoidal()
  {
    this->setStepperName("RK Explicit Trapezoidal");
    this->setStepperType("RK Explicit Trapezoidal");
    this->setupTableau();
    this->setupDefault();
    this->setUseFSAL(            false);
    this->setICConsistency(      "None");
    this->setICConsistencyCheck( false);
  }

  StepperERK_Trapezoidal(
    const Teuchos::RCP<const Thyra::ModelEvaluator<Scalar> >& appModel,
    bool useFSAL,
    std::string ICConsistency,
    bool ICConsistencyCheck,
    bool useEmbedded,
    const Teuchos::RCP<StepperRKAppAction<Scalar> >& stepperRKAppAction)
  {
    this->setStepperName("RK Explicit Trapezoidal");
    this->setStepperType("RK Explicit Trapezoidal");
    this->setupTableau();
    this->setup(appModel, useFSAL, ICConsistency,
                ICConsistencyCheck, useEmbedded, stepperRKAppAction);
  }

  std::string getDescription() const
  {
    std::ostringstream Description;
    Description << this->getStepperType() << "\n"
                << "This Stepper is known as 'RK Explicit Trapezoidal' or 'Heuns Method' or 'SSPERK22'.\n"
                << "c      = [  0    1  ]'\n"
                << "A      = [  0       ]\n"
                << "         [  1    0  ]\n"
                << "b      = [ 1/2  1/2 ]\n"
                << "bstart = [ 3/4  1/4 ]'";
    return Description.str();
  }

protected:

  void setupTableau()
  {
   typedef Teuchos::ScalarTraits<Scalar> ST;
    using Teuchos::as;
    const Scalar one = ST::one();
    const Scalar zero = ST::zero();
    const Scalar onehalf = one/(2*one);

    int NumStages = 2;
    Teuchos::SerialDenseMatrix<int,Scalar> A(NumStages,NumStages);
    Teuchos::SerialDenseVector<int,Scalar> b(NumStages);
    Teuchos::SerialDenseVector<int,Scalar> c(NumStages);
    Teuchos::SerialDenseVector<int,Scalar> bstar(NumStages);

    // Fill A:
    A(0,0) = zero; A(0,1) = zero;
    A(1,0) =  one; A(1,1) = zero;

    // Fill b:
    b(0) = onehalf; b(1) = onehalf;

    // fill c:
    c(0) = zero; c(1) = one;

    // Fill bstar
    bstar(0) = as<Scalar>(3*one/(4*one));
    bstar(1) = as<Scalar>(1*one/(4*one));

    int order = 2;

    this->tableau_ = Teuchos::rcp(new RKButcherTableau<Scalar>(
      this->getStepperType(),A,b,c,order,order,order,bstar));
    this->tableau_->setTVD(true);
    this->tableau_->setTVDCoeff(1.0);
  }
};


// ------------------------------------------------------------------------
template<class Scalar>
Teuchos::RCP<StepperERK_Trapezoidal<Scalar> >
createStepperERK_Trapezoidal(
  const Teuchos::RCP<const Thyra::ModelEvaluator<Scalar> >& model,
  Teuchos::RCP<Teuchos::ParameterList> pl)
{
  auto stepper = Teuchos::rcp(new StepperERK_Trapezoidal<Scalar>());

  // Test for aliases.
  if (pl != Teuchos::null) {
    auto stepperType =
      pl->get<std::string>("Stepper Type", stepper->getStepperType());

    TEUCHOS_TEST_FOR_EXCEPTION(
      stepperType != stepper->getStepperType() &&
      stepperType != "Heuns Method" && stepperType != "SSPERK22" &&
      stepperType != "SSPRK2", std::logic_error,
      "  ParameterList 'Stepper Type' (='" + stepperType +"')\n"
      "  does not match type for this Stepper (='" + stepper->getStepperType() +
      "')\n  or one of its aliases ('Heuns Method' or 'SSPERK22' or 'SSPRK2').\n");

    // Reset default StepperType.
    pl->set<std::string>("Stepper Type", stepper->getStepperType());
  }

  stepper->setStepperRKValues(pl);

  if (model != Teuchos::null) {
    stepper->setModel(model);
    stepper->initialize();
  }

  return stepper;
}


// ----------------------------------------------------------------------------
template<class Scalar>
class StepperERK_SSPERK54 :
  virtual public StepperExplicitRK<Scalar>
{
  public:
  StepperERK_SSPERK54()
  {
    this->setStepperName("SSPERK54");
    this->setStepperType("SSPERK54");
    this->setupTableau();
    this->setupDefault();
    this->setUseFSAL(            false);
    this->setICConsistency(      "None");
    this->setICConsistencyCheck( false);
  }

  StepperERK_SSPERK54(
    const Teuchos::RCP<const Thyra::ModelEvaluator<Scalar> >& appModel,
    bool useFSAL,
    std::string ICConsistency,
    bool ICConsistencyCheck,
    bool useEmbedded,
    const Teuchos::RCP<StepperRKAppAction<Scalar> >& stepperRKAppAction)
  {
    this->setStepperName("SSPERK54");
    this->setStepperType("SSPERK54");
    this->setupTableau();
    this->setup(appModel, useFSAL, ICConsistency,
                ICConsistencyCheck, useEmbedded, stepperRKAppAction);
  }

  std::string getDescription() const
  {
    std::ostringstream Description;
    Description << this->getStepperType() << "\n"
                << "Strong Stability Preserving Explicit RK (stage=5, order=4)\n"
                << std::endl;
    return Description.str();
  }

protected:

  void setupTableau()
  {

    typedef Teuchos::ScalarTraits<Scalar> ST;
    using Teuchos::as;
    const int NumStages = 5;
    const int order     = 4;
    const Scalar sspcoef = 1.5082;
    Teuchos::SerialDenseMatrix<int,Scalar> A(NumStages,NumStages);
    Teuchos::SerialDenseVector<int,Scalar> b(NumStages);
    Teuchos::SerialDenseVector<int,Scalar> c(NumStages);
    Teuchos::SerialDenseVector<int,Scalar> bstar(NumStages);
    const Scalar zero = ST::zero();

    // Fill A:
    A(0,0) = A(0,1) =  A(0,2) = A(0,3) = A(0,4) = zero;

    A(1,0) = as<Scalar>(0.391752226571889);
    A(1,1) = A(1,2) = A(1,3) = A(0,4) = zero;

    A(2,0) = as<Scalar>(0.217669096261169);
    A(2,1) = as<Scalar>(0.368410593050372);
    A(2,2) = A(2,3) = A(2,4) = zero;

    A(3,0) = as<Scalar>(0.082692086657811);
    A(3,1) = as<Scalar>(0.139958502191896);
    A(3,2) = as<Scalar>(0.251891774271693);
    A(3,3) = A(2,4) = zero;

    A(4,0) = as<Scalar>(0.067966283637115);
    A(4,1) = as<Scalar>(0.115034698504632);
    A(4,2) = as<Scalar>(0.207034898597385);
    A(4,3) = as<Scalar>(0.544974750228520);
    A(4,4) = zero;

    // Fill b:
    b(0) = as<Scalar>(0.146811876084786);
    b(1) = as<Scalar>(0.248482909444976);
    b(2) = as<Scalar>(0.104258830331980);
    b(3) = as<Scalar>(0.274438900901350);
    b(4) = as<Scalar>(0.226007483236908);

    // fill c:
    c(0) = zero;
    c(1) = A(1,0);
    c(2) = A(2,0) + A(2,1);
    c(3) = A(3,0) + A(3,1) + A(3,2);
    c(4) = A(4,0) + A(4,1) + A(4,2) + A(4,3);

    // Fill bstar:
    bstar(0) = as<Scalar>(0.130649104813131);
    bstar(1) = as<Scalar>(0.317716031201302);
    bstar(2) = as<Scalar>(0.000000869337261);
    bstar(3) = as<Scalar>(0.304581512634772);
    bstar(4) = as<Scalar>(0.247052482013534);

    this->tableau_ = Teuchos::rcp(new RKButcherTableau<Scalar>(
      this->getStepperType(),A,b,c,order,order,order,bstar));
    this->tableau_->setTVD(true);
    this->tableau_->setTVDCoeff(sspcoef);
  }
};


// ------------------------------------------------------------------------
template<class Scalar>
Teuchos::RCP<StepperERK_SSPERK54<Scalar> >
createStepperERK_SSPERK54(
  const Teuchos::RCP<const Thyra::ModelEvaluator<Scalar> >& model,
  Teuchos::RCP<Teuchos::ParameterList> pl)
{
  auto stepper = Teuchos::rcp(new StepperERK_SSPERK54<Scalar>());
  stepper->setStepperRKValues(pl);

  if (model != Teuchos::null) {
    stepper->setModel(model);
    stepper->initialize();
  }

  return stepper;
}


// ----------------------------------------------------------------------------
template<class Scalar>
class StepperERK_General :
  virtual public StepperExplicitRK<Scalar>
{
public:
  StepperERK_General()
  {
    this->setStepperName("General ERK");
    this->setStepperType("General ERK");
    this->setupTableau();
    this->setupDefault();
    this->setUseFSAL(            false);
    this->setICConsistency(      "None");
    this->setICConsistencyCheck( false);
  }

  StepperERK_General(
    const Teuchos::RCP<const Thyra::ModelEvaluator<Scalar> >& appModel,
    bool useFSAL,
    std::string ICConsistency,
    bool ICConsistencyCheck,
    bool useEmbedded,
    const Teuchos::SerialDenseMatrix<int,Scalar>& A,
    const Teuchos::SerialDenseVector<int,Scalar>& b,
    const Teuchos::SerialDenseVector<int,Scalar>& c,
    const int order,
    const int orderMin,
    const int orderMax,
    const Teuchos::SerialDenseVector<int,Scalar>& bstar,
    const Teuchos::RCP<StepperRKAppAction<Scalar> >& stepperRKAppAction=Teuchos::null)
  {
    this->setStepperName("General ERK");
    this->setStepperType("General ERK");
    this->setTableau(A,b,c,order,orderMin,orderMax,bstar);

    TEUCHOS_TEST_FOR_EXCEPTION(
      this->tableau_->isImplicit() == true, std::logic_error,
      "Error - General ERK received an implicit Butcher Tableau!\n");

    this->setup(appModel, useFSAL, ICConsistency,
                ICConsistencyCheck, useEmbedded, stepperRKAppAction);
  }

  virtual std::string getDescription() const
  {
    std::stringstream Description;
    Description << this->getStepperType() << "\n"
      << "The format of the Butcher Tableau parameter list is\n"
      << "  <Parameter name=\"A\" type=\"string\" value=\"# # # ;\n"
      << "                                           # # # ;\n"
      << "                                           # # #\"/>\n"
      << "  <Parameter name=\"b\" type=\"string\" value=\"# # #\"/>\n"
      << "  <Parameter name=\"c\" type=\"string\" value=\"# # #\"/>\n\n"
      << "Note the number of stages is implicit in the number of entries.\n"
      << "The number of stages must be consistent.\n"
      << "\n"
      << "Default tableau is RK4 (order=4):\n"
      << "c = [  0  1/2 1/2  1  ]'\n"
      << "A = [  0              ]\n"
      << "    [ 1/2  0          ]\n"
      << "    [  0  1/2  0      ]\n"
      << "    [  0   0   1   0  ]\n"
      << "b = [ 1/6 1/3 1/3 1/6 ]'";
    return Description.str();
  }

  void setupTableau()
  {
    if (this->tableau_ == Teuchos::null) {
      // Set tableau to the default if null, otherwise keep current tableau.
      auto stepper = Teuchos::rcp(new StepperERK_4Stage4thOrder<Scalar>());
      auto t = stepper->getTableau();
      this->tableau_ = Teuchos::rcp(new RKButcherTableau<Scalar>(
                                 this->getStepperType(),
                                 t->A(),t->b(),t->c(),
                                 t->order(),t->orderMin(),t->orderMax(),
                                 t->bstar()));
    }
  }

  void setTableau(const Teuchos::SerialDenseMatrix<int,Scalar>& A,
                  const Teuchos::SerialDenseVector<int,Scalar>& b,
                  const Teuchos::SerialDenseVector<int,Scalar>& c,
                  const int order,
                  const int orderMin,
                  const int orderMax,
                  const Teuchos::SerialDenseVector<int,Scalar>&
                    bstar = Teuchos::SerialDenseVector<int,Scalar>())
  {
    this->tableau_ = Teuchos::rcp(new RKButcherTableau<Scalar>(
      this->getStepperType(),A,b,c,order,orderMin,orderMax,bstar));
    this->isInitialized_ = false;
  }

  virtual std::string getDefaultICConsistency() const { return "Consistent"; }

  Teuchos::RCP<const Teuchos::ParameterList>
  getValidParameters() const
  {
    auto pl = this->getValidParametersBasicERK();

    // Tableau ParameterList
    Teuchos::SerialDenseMatrix<int,Scalar> A = this->tableau_->A();
    Teuchos::SerialDenseVector<int,Scalar> b = this->tableau_->b();
    Teuchos::SerialDenseVector<int,Scalar> c = this->tableau_->c();
    Teuchos::SerialDenseVector<int,Scalar> bstar = this->tableau_->bstar();

    Teuchos::RCP<Teuchos::ParameterList> tableauPL = Teuchos::parameterList();

    std::ostringstream Astream;
    Astream.precision(15);
    for(int i = 0; i < A.numRows(); i++) {
      for(int j = 0; j < A.numCols()-1; j++) {
        Astream << A(i,j) << " ";
      }
      Astream << A(i,A.numCols()-1);
      if ( i != A.numRows()-1 ) Astream << "; ";
    }
    tableauPL->set<std::string>("A", Astream.str());

    std::ostringstream bstream;
    bstream.precision(15);
    for(int i = 0; i < b.length()-1; i++) {
        bstream << b(i) << " ";
    }
    bstream << b(b.length()-1);
    tableauPL->set<std::string>("b", bstream.str());

    std::ostringstream cstream;
    cstream.precision(15);
    for(int i = 0; i < c.length()-1; i++) {
        cstream << c(i) << " ";
    }
    cstream << c(c.length()-1);
    tableauPL->set<std::string>("c", cstream.str());

    tableauPL->set<int>("order", this->tableau_->order());

    if ( bstar.length() == 0 ) {
      tableauPL->set("bstar", "");
    } else {
      std::ostringstream bstarstream;
      bstarstream.precision(15);
      for(int i = 0; i < bstar.length()-1; i++) {
          bstarstream << bstar(i) << " ";
      }
      bstarstream << bstar(bstar.length()-1);
      tableauPL->set<std::string>("bstar", bstarstream.str());
    }

    pl->set("Tableau", *tableauPL);

    return pl;
  }
};


// ------------------------------------------------------------------------
template<class Scalar>
Teuchos::RCP<StepperERK_General<Scalar> >
createStepperERK_General(
  const Teuchos::RCP<const Thyra::ModelEvaluator<Scalar> >& model,
  Teuchos::RCP<Teuchos::ParameterList> pl)
{
  auto stepper = Teuchos::rcp(new StepperERK_General<Scalar>());
  stepper->setStepperRKValues(pl);

  if (pl != Teuchos::null) {
    if (pl->isParameter("Tableau")) {
      auto t = stepper->createTableau(pl);
      stepper->setTableau( t->A(),t->b(),t->c(),
                           t->order(),t->orderMin(),t->orderMax(),
                           t->bstar() );
    }
  }
  TEUCHOS_TEST_FOR_EXCEPTION(
    stepper->getTableau()->isImplicit() == true, std::logic_error,
    "Error - General ERK received an implicit Butcher Tableau!\n");

  if (model != Teuchos::null) {
    stepper->setModel(model);
    stepper->initialize();
  }

  return stepper;
}


// ----------------------------------------------------------------------------
template<class Scalar>
class StepperDIRK_BackwardEuler :
  virtual public StepperDIRK<Scalar>
{
public:
  StepperDIRK_BackwardEuler()
  {
    this->setStepperName("RK Backward Euler");
    this->setStepperType("RK Backward Euler");
    this->setupTableau();
    this->setupDefault();
    this->setUseFSAL(            false);
    this->setICConsistency(      "None");
    this->setICConsistencyCheck( false);
  }

  StepperDIRK_BackwardEuler(
    const Teuchos::RCP<const Thyra::ModelEvaluator<Scalar> >& appModel,
    const Teuchos::RCP<Thyra::NonlinearSolverBase<Scalar> >& solver,
    bool useFSAL,
    std::string ICConsistency,
    bool ICConsistencyCheck,
    bool useEmbedded,
    bool zeroInitialGuess,
    const Teuchos::RCP<StepperRKAppAction<Scalar> >& stepperRKAppAction)
  {
    this->setStepperName("RK Backward Euler");
    this->setStepperType("RK Backward Euler");
    this->setupTableau();
    this->setup(appModel, solver, useFSAL, ICConsistency, ICConsistencyCheck,
                useEmbedded, zeroInitialGuess, stepperRKAppAction);
  }

  std::string getDescription() const
  {
    std::ostringstream Description;
    Description << this->getStepperType() << "\n"
                << "c = [ 1 ]'\n"
                << "A = [ 1 ]\n"
                << "b = [ 1 ]'";
    return Description.str();
  }

protected:

  void setupTableau()
  {
    typedef Teuchos::ScalarTraits<Scalar> ST;
    const Scalar sspcoef = std::numeric_limits<Scalar>::max();
    int NumStages = 1;
    Teuchos::SerialDenseMatrix<int,Scalar> A(NumStages,NumStages);
    Teuchos::SerialDenseVector<int,Scalar> b(NumStages);
    Teuchos::SerialDenseVector<int,Scalar> c(NumStages);

    // Fill A:
    A(0,0) = ST::one();

    // Fill b:
    b(0) = ST::one();

    // Fill c:
    c(0) = ST::one();

    int order = 1;

    this->tableau_ = Teuchos::rcp(new RKButcherTableau<Scalar>(
      this->getStepperType(),A,b,c,order,order,order));
    this->tableau_->setTVD(true);
    this->tableau_->setTVDCoeff(sspcoef);
  }
};


// ------------------------------------------------------------------------
template<class Scalar>
Teuchos::RCP<StepperDIRK_BackwardEuler<Scalar> >
createStepperDIRK_BackwardEuler(
  const Teuchos::RCP<const Thyra::ModelEvaluator<Scalar> >& model,
  Teuchos::RCP<Teuchos::ParameterList> pl)
{
  auto stepper = Teuchos::rcp(new StepperDIRK_BackwardEuler<Scalar>());
  stepper->setStepperDIRKValues(pl);

  if (model != Teuchos::null) {
    stepper->setModel(model);
    stepper->initialize();
  }

  return stepper;
}


// ----------------------------------------------------------------------------
template<class Scalar>
class StepperSDIRK_2Stage2ndOrder :
  virtual public StepperDIRK<Scalar>
{
public:
  StepperSDIRK_2Stage2ndOrder()
  {
    typedef Teuchos::ScalarTraits<Scalar> ST;
    const Scalar one = ST::one();
    gammaDefault_ = Teuchos::as<Scalar>((2*one-ST::squareroot(2*one))/(2*one));

    this->setStepperName("SDIRK 2 Stage 2nd order");
    this->setStepperType("SDIRK 2 Stage 2nd order");
    this->setGamma(gammaDefault_);
    this->setupTableau();
    this->setupDefault();
    this->setUseFSAL(            false);
    this->setICConsistency(      "None");
    this->setICConsistencyCheck( false);
  }

  StepperSDIRK_2Stage2ndOrder(
    const Teuchos::RCP<const Thyra::ModelEvaluator<Scalar> >& appModel,
    const Teuchos::RCP<Thyra::NonlinearSolverBase<Scalar> >& solver,
    bool useFSAL,
    std::string ICConsistency,
    bool ICConsistencyCheck,
    bool useEmbedded,
    bool zeroInitialGuess,
    const Teuchos::RCP<StepperRKAppAction<Scalar> >& stepperRKAppAction,
    Scalar gamma = Scalar(0.2928932188134524))
  {
    typedef Teuchos::ScalarTraits<Scalar> ST;
    const Scalar one = ST::one();
    gammaDefault_ = Teuchos::as<Scalar>((2*one-ST::squareroot(2*one))/(2*one));

    this->setStepperName("SDIRK 2 Stage 2nd order");
    this->setStepperType("SDIRK 2 Stage 2nd order");
    this->setGamma(gamma);
    this->setupTableau();
    this->setup(appModel, solver, useFSAL, ICConsistency, ICConsistencyCheck,
                useEmbedded, zeroInitialGuess, stepperRKAppAction);
  }

  void setGamma(Scalar gamma)
  {
    gamma_ = gamma;
    this->isInitialized_ = false;
    this->setupTableau();
  }

  Scalar getGamma() const { return gamma_; }

  std::string getDescription() const
  {
    std::ostringstream Description;
    Description << this->getStepperType() << "\n"
                << "Computer Methods for ODEs and DAEs\n"
                << "U. M. Ascher and L. R. Petzold\n"
                << "p. 106\n"
                << "gamma = (2+-sqrt(2))/2\n"
                << "c = [  gamma   1     ]'\n"
                << "A = [  gamma   0     ]\n"
                << "    [ 1-gamma  gamma ]\n"
                << "b = [ 1-gamma  gamma ]'";
    return Description.str();
  }

  Teuchos::RCP<const Teuchos::ParameterList>
  getValidParameters() const
  {
    auto pl = this->getValidParametersBasicDIRK();
    pl->template set<double>("gamma", this->getGamma(),
      "The default value is gamma = (2-sqrt(2))/2. "
      "This will produce an L-stable 2nd order method with the stage "
      "times within the timestep.  Other values of gamma will still "
      "produce an L-stable scheme, but will only be 1st order accurate.");

    return pl;
  }

protected:

  void setupTableau()
  {
    typedef Teuchos::ScalarTraits<Scalar> ST;
    int NumStages = 2;
    Teuchos::SerialDenseMatrix<int,Scalar> A(NumStages,NumStages);
    Teuchos::SerialDenseVector<int,Scalar> b(NumStages);
    Teuchos::SerialDenseVector<int,Scalar> c(NumStages);

    const Scalar one = ST::one();
    const Scalar zero = ST::zero();

    // Fill A:
    A(0,0) =                              gamma_; A(0,1) = zero;
    A(1,0) = Teuchos::as<Scalar>( one - gamma_ ); A(1,1) = gamma_;

    // Fill b:
    b(0) = Teuchos::as<Scalar>( one - gamma_ ); b(1) = gamma_;

    // Fill c:
    c(0) = gamma_; c(1) = one;

    int order = 1;
    if ( std::abs((gamma_-gammaDefault_)/gamma_) < 1.0e-08 ) order = 2;

    this->tableau_ = Teuchos::rcp(new RKButcherTableau<Scalar>(
      this->getStepperType(),A,b,c,order,order,order));
  }

  private:
    Scalar gammaDefault_;
    Scalar gamma_;
};


// ------------------------------------------------------------------------
template<class Scalar>
Teuchos::RCP<StepperSDIRK_2Stage2ndOrder<Scalar> >
createStepperSDIRK_2Stage2ndOrder(
  const Teuchos::RCP<const Thyra::ModelEvaluator<Scalar> >& model,
  Teuchos::RCP<Teuchos::ParameterList> pl)
{
  auto stepper = Teuchos::rcp(new StepperSDIRK_2Stage2ndOrder<Scalar>());
  stepper->setStepperDIRKValues(pl);

  if (pl != Teuchos::null)
    stepper->setGamma(pl->get<double>("gamma", 0.2928932188134524));

  if (model != Teuchos::null) {
    stepper->setModel(model);
    stepper->initialize();
  }

  return stepper;
}


// ----------------------------------------------------------------------------
template<class Scalar>
class StepperSDIRK_3Stage2ndOrder :
  virtual public StepperDIRK<Scalar>
{
public:
  StepperSDIRK_3Stage2ndOrder()
  {
    this->setStepperName("SDIRK 3 Stage 2nd order");
    this->setStepperType("SDIRK 3 Stage 2nd order");
    this->setupTableau();
    this->setupDefault();
    this->setUseFSAL(            false);
    this->setICConsistency(      "None");
    this->setICConsistencyCheck( false);
  }

  StepperSDIRK_3Stage2ndOrder(
    const Teuchos::RCP<const Thyra::ModelEvaluator<Scalar> >& appModel,
    const Teuchos::RCP<Thyra::NonlinearSolverBase<Scalar> >& solver,
    bool useFSAL,
    std::string ICConsistency,
    bool ICConsistencyCheck,
    bool useEmbedded,
    bool zeroInitialGuess,
    const Teuchos::RCP<StepperRKAppAction<Scalar> >& stepperRKAppAction)
  {

    this->setStepperName("SDIRK 3 Stage 2nd order");
    this->setStepperType("SDIRK 3 Stage 2nd order");
    this->setupTableau();
    this->setup(appModel, solver, useFSAL, ICConsistency, ICConsistencyCheck,
                useEmbedded, zeroInitialGuess, stepperRKAppAction);
  }

  std::string getDescription() const
  {
    std::ostringstream Description;
    Description << this->getStepperType() << "\n"
                << "Implicit-explicit Runge-Kutta schemes and applications to\n"
                << "hyperbolic systems with relaxation\n"
                << "L Pareschi, G Russo\n"
                << "Journal of Scientific computing, 2005 - Springer\n"
                << "Table 5\n"
                << "gamma = 1/(2+sqrt(2))\n"
                << "c = [  gamma   (1-gamma)   1/2  ]'\n"
                << "A = [  gamma      0         0   ]\n"
                << "    [ 1-2gamma   gamma      0   ]\n"
                << "    [ 1/2-gamma   0      gamma  ]\n"
                << "b = [  1/6       1/6       2/3  ]'";
    return Description.str();
  }

protected:

  void setupTableau()
  {
    typedef Teuchos::ScalarTraits<Scalar> ST;
    using Teuchos::as;
    const int NumStages = 3;
    const int order = 2;
    const Scalar sspcoef = 1.0529;
    Teuchos::SerialDenseMatrix<int,Scalar> A(NumStages,NumStages);
    Teuchos::SerialDenseVector<int,Scalar> b(NumStages);
    Teuchos::SerialDenseVector<int,Scalar> c(NumStages);
    const Scalar one = ST::one();
    const Scalar zero = ST::zero();
    const Scalar gamma = as<Scalar>(one - ( one / ST::squareroot(2*one) ) );

    // Fill A:
    A(0,0) = A(1,1) = A(2,2) = gamma;
    A(0,1) = A(0,2) = A(1,2) = A(2,1) = zero;
    A(1,0) = as<Scalar>(one - 2*gamma);
    A(2,0) = as<Scalar>(  ( one/ (2.*one)) - gamma );

    // Fill b:
    b(0) = b(1) = ( one / (6*one) );
    b(2) = (2*one)/(3*one);

    // Fill c:
    c(0) = gamma;
    c(1) = one - gamma;
    c(2) = one / (2*one);

    this->tableau_ = Teuchos::rcp(new RKButcherTableau<Scalar>(
      this->getStepperType(),A,b,c,order,order,order));
    this->tableau_->setTVD(true);
    this->tableau_->setTVDCoeff(sspcoef);
  }

};


// ------------------------------------------------------------------------
template<class Scalar>
Teuchos::RCP<StepperSDIRK_3Stage2ndOrder<Scalar> >
createStepperSDIRK_3Stage2ndOrder(
  const Teuchos::RCP<const Thyra::ModelEvaluator<Scalar> >& model,
  Teuchos::RCP<Teuchos::ParameterList> pl)
{
  auto stepper = Teuchos::rcp(new StepperSDIRK_3Stage2ndOrder<Scalar>());
  stepper->setStepperDIRKValues(pl);

  if (model != Teuchos::null) {
    stepper->setModel(model);
    stepper->initialize();
  }

  return stepper;
}


// ----------------------------------------------------------------------------
template<class Scalar>
class StepperSDIRK_2Stage3rdOrder :
  virtual public StepperDIRK<Scalar>
{
public:
  StepperSDIRK_2Stage3rdOrder()
  {
    typedef Teuchos::ScalarTraits<Scalar> ST;
    using Teuchos::as;
    const Scalar one = ST::one();
    gammaDefault_ = as<Scalar>((3*one+ST::squareroot(3*one))/(6*one));
    gammaTypeDefault_ = "3rd Order A-stable";

    this->setStepperName("SDIRK 2 Stage 3rd order");
    this->setStepperType("SDIRK 2 Stage 3rd order");
    this->setGammaType(gammaTypeDefault_);
    this->setGamma(gammaDefault_);
    this->setupTableau();
    this->setupDefault();
    this->setUseFSAL(            false);
    this->setICConsistency(      "None");
    this->setICConsistencyCheck( false);
  }

  StepperSDIRK_2Stage3rdOrder(
    const Teuchos::RCP<const Thyra::ModelEvaluator<Scalar> >& appModel,
    const Teuchos::RCP<Thyra::NonlinearSolverBase<Scalar> >& solver,
    bool useFSAL,
    std::string ICConsistency,
    bool ICConsistencyCheck,
    bool useEmbedded,
    bool zeroInitialGuess,
    const Teuchos::RCP<StepperRKAppAction<Scalar> >& stepperRKAppAction,
    std::string gammaType = "3rd Order A-stable",
    Scalar gamma = Scalar(0.7886751345948128))
  {
    typedef Teuchos::ScalarTraits<Scalar> ST;
    using Teuchos::as;
    const Scalar one = ST::one();
    gammaDefault_ = as<Scalar>((3*one+ST::squareroot(3*one))/(6*one));
    gammaTypeDefault_ = "3rd Order A-stable";

    this->setStepperName("SDIRK 2 Stage 3rd order");
    this->setStepperType("SDIRK 2 Stage 3rd order");
    this->setGammaType(gammaType);
    this->setGamma(gamma);
    this->setupTableau();
    this->setup(appModel, solver, useFSAL, ICConsistency, ICConsistencyCheck,
                useEmbedded, zeroInitialGuess, stepperRKAppAction);
  }

  void setGammaType(std::string gammaType)
  {
    TEUCHOS_TEST_FOR_EXCEPTION(
      !(gammaType == "3rd Order A-stable" ||
        gammaType == "2nd Order L-stable" ||
        gammaType == "gamma"), std::logic_error,
      "gammaType needs to be '3rd Order A-stable', '2nd Order L-stable' "
      "or 'gamma'.");

    gammaType_ = gammaType;
    this->isInitialized_ = false;
    this->setupTableau();
  }

  std::string getGammaType() const { return gammaType_; }

  void setGamma(Scalar gamma)
  {
    if ( gammaType_ == "gamma" ) {
      gamma_ = gamma;
      this->setupTableau();
    }
    this->isInitialized_ = false;
  }

  Scalar getGamma() const { return gamma_; }

  std::string getDescription() const
  {
    std::ostringstream Description;
    Description << this->getStepperType() << "\n"
                << "Solving Ordinary Differential Equations I:\n"
                << "Nonstiff Problems, 2nd Revised Edition\n"
                << "E. Hairer, S. P. Norsett, and G. Wanner\n"
                << "Table 7.2, pg 207\n"
                << "gamma = (3+sqrt(3))/6 -> 3rd order and A-stable\n"
                << "gamma = (2-sqrt(2))/2 -> 2nd order and L-stable\n"
                << "c = [  gamma     1-gamma  ]'\n"
                << "A = [  gamma     0        ]\n"
                << "    [ 1-2*gamma  gamma    ]\n"
                << "b = [ 1/2        1/2      ]'";
    return Description.str();
  }

  Teuchos::RCP<const Teuchos::ParameterList>
  getValidParameters() const
  {
    auto pl = this->getValidParametersBasicDIRK();

    pl->template set<std::string>("Gamma Type", this->getGammaType(),
      "Valid values are '3rd Order A-stable' ((3+sqrt(3))/6.), "
      "'2nd Order L-stable' ((2-sqrt(2))/2) and 'gamma' (user defined).  "
      "The default value is '3rd Order A-stable'.");
    pl->template set<double>("gamma", this->getGamma(),
      "Equal to (3+sqrt(3))/6 if 'Gamma Type' = '3rd Order A-stable', or "
      "(2-sqrt(2))/2 if 'Gamma Type' = '2nd Order L-stable', or "
      "user-defined gamma value if 'Gamma Type = 'gamma'.  "
      "The default value is gamma = (3+sqrt(3))/6, which matches "
      "the default 'Gamma Type' = '3rd Order A-stable'.");

    return pl;
  }

protected:

  void setupTableau()
  {
    typedef Teuchos::ScalarTraits<Scalar> ST;
    using Teuchos::as;
    int NumStages = 2;
    Teuchos::SerialDenseMatrix<int,Scalar> A(NumStages,NumStages);
    Teuchos::SerialDenseVector<int,Scalar> b(NumStages);
    Teuchos::SerialDenseVector<int,Scalar> c(NumStages);
    const Scalar one = ST::one();
    const Scalar zero = ST::zero();

    int order = 0;
    if (gammaType_ == "3rd Order A-stable") {
      order = 3;
      gamma_ = gammaDefault_;
    } else if (gammaType_ == "2nd Order L-stable") {
      order = 2;
      gamma_ = as<Scalar>( (2*one - ST::squareroot(2*one))/(2*one) );
    } else if (gammaType_ == "gamma") {
      order = 2;
    }

    // Fill A:
    A(0,0) =                     gamma_; A(0,1) = zero;
    A(1,0) = as<Scalar>(one - 2*gamma_); A(1,1) = gamma_;

    // Fill b:
    b(0) = as<Scalar>( one/(2*one) ); b(1) = as<Scalar>( one/(2*one) );

    // Fill c:
    c(0) = gamma_; c(1) = as<Scalar>( one - gamma_ );

    this->tableau_ = Teuchos::rcp(new RKButcherTableau<Scalar>(
      this->getStepperType(),A,b,c,order,2,3));
  }

  private:
    std::string gammaTypeDefault_;
    std::string gammaType_;
    Scalar      gammaDefault_;
    Scalar      gamma_;
};


// ------------------------------------------------------------------------
template<class Scalar>
Teuchos::RCP<StepperSDIRK_2Stage3rdOrder<Scalar> >
createStepperSDIRK_2Stage3rdOrder(
  const Teuchos::RCP<const Thyra::ModelEvaluator<Scalar> >& model,
  Teuchos::RCP<Teuchos::ParameterList> pl)
{
  auto stepper = Teuchos::rcp(new StepperSDIRK_2Stage3rdOrder<Scalar>());
  stepper->setStepperDIRKValues(pl);

  if (pl != Teuchos::null) {
    stepper->setGammaType(
      pl->get<std::string>("Gamma Type", "3rd Order A-stable"));
    stepper->setGamma(pl->get<double>("gamma", 0.7886751345948128));
  }

  if (model != Teuchos::null) {
    stepper->setModel(model);
    stepper->initialize();
  }

  return stepper;
}


// ----------------------------------------------------------------------------
template<class Scalar>
class StepperEDIRK_2Stage3rdOrder :
  virtual public StepperDIRK<Scalar>
{
public:
  StepperEDIRK_2Stage3rdOrder()
  {
    this->setStepperName("EDIRK 2 Stage 3rd order");
    this->setStepperType("EDIRK 2 Stage 3rd order");
    this->setupTableau();
    this->setupDefault();
    this->setUseFSAL(            false);
    this->setICConsistency(      "None");
    this->setICConsistencyCheck( false);
  }

  StepperEDIRK_2Stage3rdOrder(
    const Teuchos::RCP<const Thyra::ModelEvaluator<Scalar> >& appModel,
    const Teuchos::RCP<Thyra::NonlinearSolverBase<Scalar> >& solver,
    bool useFSAL,
    std::string ICConsistency,
    bool ICConsistencyCheck,
    bool useEmbedded,
    bool zeroInitialGuess,
    const Teuchos::RCP<StepperRKAppAction<Scalar> >& stepperRKAppAction)
  {
    this->setStepperName("EDIRK 2 Stage 3rd order");
    this->setStepperType("EDIRK 2 Stage 3rd order");
    this->setupTableau();
    this->setup(appModel, solver, useFSAL, ICConsistency, ICConsistencyCheck,
                useEmbedded, zeroInitialGuess, stepperRKAppAction);
  }

  std::string getDescription() const
  {
    std::ostringstream Description;
    Description << this->getStepperType() << "\n"
                << "Hammer & Hollingsworth method\n"
                << "Solving Ordinary Differential Equations I:\n"
                << "Nonstiff Problems, 2nd Revised Edition\n"
                << "E. Hairer, S. P. Norsett, and G. Wanner\n"
                << "Table 7.1, pg 205\n"
                << "c = [  0   2/3 ]'\n"
                << "A = [  0    0  ]\n"
                << "    [ 1/3  1/3 ]\n"
                << "b = [ 1/4  3/4 ]'";
    return Description.str();
  }

protected:

  void setupTableau()
  {
    typedef Teuchos::ScalarTraits<Scalar> ST;
    using Teuchos::as;
    int NumStages = 2;
    Teuchos::SerialDenseMatrix<int,Scalar> A(NumStages,NumStages);
    Teuchos::SerialDenseVector<int,Scalar> b(NumStages);
    Teuchos::SerialDenseVector<int,Scalar> c(NumStages);
    const Scalar one = ST::one();
    const Scalar zero = ST::zero();

    // Fill A:
    A(0,0) =                      zero; A(0,1) =                      zero;
    A(1,0) = as<Scalar>( one/(3*one) ); A(1,1) = as<Scalar>( one/(3*one) );

    // Fill b:
    b(0) = as<Scalar>( one/(4*one) ); b(1) = as<Scalar>( 3*one/(4*one) );

    // Fill c:
    c(0) = zero; c(1) = as<Scalar>( 2*one/(3*one) );
    int order = 3;

    this->tableau_ = Teuchos::rcp(new RKButcherTableau<Scalar>(
      this->getStepperType(),A,b,c,order,order,order));
    this->tableau_->setTVD(true);
    this->tableau_->setTVDCoeff(1.5);
  }
};


template<class Scalar>
Teuchos::RCP<StepperEDIRK_2Stage3rdOrder<Scalar> >
createStepperEDIRK_2Stage3rdOrder(
  const Teuchos::RCP<const Thyra::ModelEvaluator<Scalar> >& model,
  Teuchos::RCP<Teuchos::ParameterList> pl)
{
  auto stepper = Teuchos::rcp(new StepperEDIRK_2Stage3rdOrder<Scalar>());
  stepper->setStepperDIRKValues(pl);

  if (model != Teuchos::null) {
    stepper->setModel(model);
    stepper->initialize();
  }

  return stepper;
}


// ----------------------------------------------------------------------------
template<class Scalar>
class StepperDIRK_1StageTheta :
  virtual public StepperDIRK<Scalar>
{
public:
  StepperDIRK_1StageTheta()
  {
    typedef Teuchos::ScalarTraits<Scalar> ST;
    thetaDefault_ = ST::one()/(2*ST::one());

    this->setStepperName("DIRK 1 Stage Theta Method");
    this->setStepperType("DIRK 1 Stage Theta Method");
    this->setTheta(thetaDefault_);
    this->setupTableau();
    this->setupDefault();
    this->setUseFSAL(            false);
    this->setICConsistency(      "None");
    this->setICConsistencyCheck( false);
  }

  StepperDIRK_1StageTheta(
    const Teuchos::RCP<const Thyra::ModelEvaluator<Scalar> >& appModel,
    const Teuchos::RCP<Thyra::NonlinearSolverBase<Scalar> >& solver,
    bool useFSAL,
    std::string ICConsistency,
    bool ICConsistencyCheck,
    bool useEmbedded,
    bool zeroInitialGuess,
    const Teuchos::RCP<StepperRKAppAction<Scalar> >& stepperRKAppAction,
    Scalar theta = Scalar(0.5))
  {
    typedef Teuchos::ScalarTraits<Scalar> ST;
    thetaDefault_ = ST::one()/(2*ST::one());

    this->setStepperName("DIRK 1 Stage Theta Method");
    this->setStepperType("DIRK 1 Stage Theta Method");
    this->setTheta(theta);
    this->setupTableau();
    this->setup(appModel, solver, useFSAL, ICConsistency, ICConsistencyCheck,
                useEmbedded, zeroInitialGuess, stepperRKAppAction);
  }

  void setTheta(Scalar theta)
  {
    TEUCHOS_TEST_FOR_EXCEPTION(
      theta == Teuchos::ScalarTraits<Scalar>::zero(), std::logic_error,
      "'theta' can not be zero, as it makes this stepper explicit. \n"
      "Try using the 'RK Forward Euler' stepper.\n");
    theta_ = theta;
    this->setupTableau();
    this->isInitialized_ = false;
  }

  Scalar getTheta() const { return theta_; }

  std::string getDescription() const
  {
    std::ostringstream Description;
    Description << this->getStepperType() << "\n"
                << "Non-standard finite-difference methods\n"
                << "in dynamical systems, P. Kama,\n"
                << "Dissertation, University of Pretoria, pg. 49.\n"
                << "Comment:  Generalized Implicit Midpoint Method\n"
                << "c = [ theta ]'\n"
                << "A = [ theta ]\n"
                << "b = [  1  ]'";
    return Description.str();
  }

  Teuchos::RCP<const Teuchos::ParameterList>
  getValidParameters() const
  {
    auto pl = this->getValidParametersBasicDIRK();

    pl->template set<double>("theta", getTheta(),
      "Valid values are 0 <= theta <= 1, where theta = 0 "
      "implies Forward Euler, theta = 1/2 implies implicit midpoint "
      "method (default), and theta = 1 implies Backward Euler. "
      "For theta != 1/2, this method is first-order accurate, "
      "and with theta = 1/2, it is second-order accurate.  "
      "This method is A-stable, but becomes L-stable with theta=1.");

    return pl;
  }

protected:

  void setupTableau()
  {
    typedef Teuchos::ScalarTraits<Scalar> ST;
    int NumStages = 1;
    Teuchos::SerialDenseMatrix<int,Scalar> A(NumStages,NumStages);
    Teuchos::SerialDenseVector<int,Scalar> b(NumStages);
    Teuchos::SerialDenseVector<int,Scalar> c(NumStages);
    A(0,0) = theta_;
    b(0) = ST::one();
    c(0) = theta_;

    int order = 1;
    if ( std::abs((theta_-thetaDefault_)/theta_) < 1.0e-08 ) order = 2;

    this->tableau_ = Teuchos::rcp(new RKButcherTableau<Scalar>(
      this->getStepperType(),A,b,c,order,1,2));
    this->tableau_->setTVD(true);
    this->tableau_->setTVDCoeff(2.0);
  }

  private:
    Scalar thetaDefault_;
    Scalar theta_;
};


// ------------------------------------------------------------------------
template<class Scalar>
Teuchos::RCP<StepperDIRK_1StageTheta<Scalar> >
createStepperDIRK_1StageTheta(
  const Teuchos::RCP<const Thyra::ModelEvaluator<Scalar> >& model,
  Teuchos::RCP<Teuchos::ParameterList> pl)
{
  auto stepper = Teuchos::rcp(new StepperDIRK_1StageTheta<Scalar>());
  stepper->setStepperDIRKValues(pl);

  if (pl != Teuchos::null) {
    stepper->setTheta(pl->get<double>("theta", 0.5));
  }

  if (model != Teuchos::null) {
    stepper->setModel(model);
    stepper->initialize();
  }

  return stepper;
}


// ----------------------------------------------------------------------------
template<class Scalar>
class StepperEDIRK_2StageTheta :
  virtual public StepperDIRK<Scalar>
{
public:
  StepperEDIRK_2StageTheta()
  {
    typedef Teuchos::ScalarTraits<Scalar> ST;
    thetaDefault_ = ST::one()/(2*ST::one());

    this->setStepperName("EDIRK 2 Stage Theta Method");
    this->setStepperType("EDIRK 2 Stage Theta Method");
    this->setTheta(thetaDefault_);
    this->setupTableau();
    this->setupDefault();
    this->setUseFSAL(            true);
    this->setICConsistency(      "Consistent");
    this->setICConsistencyCheck( false);
  }

  StepperEDIRK_2StageTheta(
    const Teuchos::RCP<const Thyra::ModelEvaluator<Scalar> >& appModel,
    const Teuchos::RCP<Thyra::NonlinearSolverBase<Scalar> >& solver,
    bool useFSAL,
    std::string ICConsistency,
    bool ICConsistencyCheck,
    bool useEmbedded,
    bool zeroInitialGuess,
    const Teuchos::RCP<StepperRKAppAction<Scalar> >& stepperRKAppAction,
    Scalar theta = Scalar(0.5))
  {
    typedef Teuchos::ScalarTraits<Scalar> ST;
    thetaDefault_ = ST::one()/(2*ST::one());

    this->setStepperName("EDIRK 2 Stage Theta Method");
    this->setStepperType("EDIRK 2 Stage Theta Method");
    this->setTheta(theta);
    this->setupTableau();
    this->setup(appModel, solver, useFSAL, ICConsistency, ICConsistencyCheck,
                useEmbedded, zeroInitialGuess, stepperRKAppAction);
  }

  void setTheta(Scalar theta)
  {
    TEUCHOS_TEST_FOR_EXCEPTION(
      theta == Teuchos::ScalarTraits<Scalar>::zero(), std::logic_error,
      "'theta' can not be zero, as it makes this stepper explicit. \n"
      "Try using the 'RK Forward Euler' stepper.\n");
    theta_ = theta;
    this->isInitialized_ = false;
    this->setupTableau();
  }

  Scalar getTheta() const { return theta_; }

  std::string getDescription() const
  {
    std::ostringstream Description;
    Description << this->getStepperType() << "\n"
                << "Computer Methods for ODEs and DAEs\n"
                << "U. M. Ascher and L. R. Petzold\n"
                << "p. 113\n"
                << "c = [  0       1     ]'\n"
                << "A = [  0       0     ]\n"
                << "    [ 1-theta  theta ]\n"
                << "b = [ 1-theta  theta ]'";
    return Description.str();
  }

  Teuchos::RCP<const Teuchos::ParameterList>
  getValidParameters() const
  {
    auto pl = this->getValidParametersBasicDIRK();

    pl->template set<double>("theta", getTheta(),
      "Valid values are 0 < theta <= 1, where theta = 0 "
      "implies Forward Euler, theta = 1/2 implies trapezoidal "
      "method (default), and theta = 1 implies Backward Euler. "
      "For theta != 1/2, this method is first-order accurate, "
      "and with theta = 1/2, it is second-order accurate.  "
      "This method is A-stable, but becomes L-stable with theta=1.");

    return pl;
  }

  void setUseFSAL(bool a) { this->useFSAL_ = a; this->isInitialized_ = false; }

protected:

  void setupTableau()
  {
    typedef Teuchos::ScalarTraits<Scalar> ST;
    const Scalar one = ST::one();
    const Scalar zero = ST::zero();

    int NumStages = 2;
    Teuchos::SerialDenseMatrix<int,Scalar> A(NumStages,NumStages);
    Teuchos::SerialDenseVector<int,Scalar> b(NumStages);
    Teuchos::SerialDenseVector<int,Scalar> c(NumStages);

    // Fill A:
    A(0,0) =                               zero; A(0,1) =  zero;
    A(1,0) = Teuchos::as<Scalar>( one - theta_ ); A(1,1) = theta_;

    // Fill b:
    b(0) = Teuchos::as<Scalar>( one - theta_ );
    b(1) = theta_;

    // Fill c:
    c(0) = zero;
    c(1) = one;

    int order = 1;
    if ( std::abs((theta_-thetaDefault_)/theta_) < 1.0e-08 ) order = 2;

    this->tableau_ = Teuchos::rcp(new RKButcherTableau<Scalar>(
      this->getStepperType(),A,b,c,order,1,2));
    this->tableau_->setTVD(true);
    this->tableau_->setTVDCoeff(2.0);
  }

  private:
    Scalar thetaDefault_;
    Scalar theta_;
};


// ------------------------------------------------------------------------
template<class Scalar>
Teuchos::RCP<StepperEDIRK_2StageTheta<Scalar> >
createStepperEDIRK_2StageTheta(
  const Teuchos::RCP<const Thyra::ModelEvaluator<Scalar> >& model,
  Teuchos::RCP<Teuchos::ParameterList> pl)
{
  auto stepper = Teuchos::rcp(new StepperEDIRK_2StageTheta<Scalar>());
  stepper->setStepperDIRKValues(pl);

  if (pl != Teuchos::null) {
    stepper->setTheta(pl->get<double>("theta", 0.5));
  }

  if (model != Teuchos::null) {
    stepper->setModel(model);
    stepper->initialize();
  }

  return stepper;
}


// ----------------------------------------------------------------------------
template<class Scalar>
class StepperEDIRK_TrapezoidalRule :
  virtual public StepperDIRK<Scalar>
{
public:
  StepperEDIRK_TrapezoidalRule()
  {
    this->setStepperName("RK Trapezoidal Rule");
    this->setStepperType("RK Trapezoidal Rule");
    this->setupTableau();
    this->setupDefault();
    this->setUseFSAL(            true);
    this->setICConsistency(      "Consistent");
    this->setICConsistencyCheck( false);
  }

  StepperEDIRK_TrapezoidalRule(
    const Teuchos::RCP<const Thyra::ModelEvaluator<Scalar> >& appModel,
    const Teuchos::RCP<Thyra::NonlinearSolverBase<Scalar> >& solver,
    bool useFSAL,
    std::string ICConsistency,
    bool ICConsistencyCheck,
    bool useEmbedded,
    bool zeroInitialGuess,
    const Teuchos::RCP<StepperRKAppAction<Scalar> >& stepperRKAppAction)
  {
    this->setStepperName("RK Trapezoidal Rule");
    this->setStepperType("RK Trapezoidal Rule");
    this->setupTableau();
    this->setup(appModel, solver, useFSAL, ICConsistency, ICConsistencyCheck,
                useEmbedded, zeroInitialGuess, stepperRKAppAction);
  }

  std::string getDescription() const
  {
    std::ostringstream Description;
    Description << this->getStepperType() << "\n"
                << "Also known as Crank-Nicolson Method.\n"
                << "c = [  0   1   ]'\n"
                << "A = [  0   0   ]\n"
                << "    [ 1/2  1/2 ]\n"
                << "b = [ 1/2  1/2 ]'";
    return Description.str();
  }

  void setUseFSAL(bool a) { this->useFSAL_ = a; this->isInitialized_ = false; }

protected:

  void setupTableau()
  {
    typedef Teuchos::ScalarTraits<Scalar> ST;
    const Scalar one = ST::one();
    const Scalar zero = ST::zero();
    const Scalar onehalf = ST::one()/(2*ST::one());

    int NumStages = 2;
    Teuchos::SerialDenseMatrix<int,Scalar> A(NumStages,NumStages);
    Teuchos::SerialDenseVector<int,Scalar> b(NumStages);
    Teuchos::SerialDenseVector<int,Scalar> c(NumStages);

    // Fill A:
    A(0,0) =    zero; A(0,1) =    zero;
    A(1,0) = onehalf; A(1,1) = onehalf;

    // Fill b:
    b(0) = onehalf;
    b(1) = onehalf;

    // Fill c:
    c(0) = zero;
    c(1) = one;

    int order = 2;

    this->tableau_ = Teuchos::rcp(new RKButcherTableau<Scalar>(
      this->getStepperType(),A,b,c,order,order,order));
    this->tableau_->setTVD(true);
    this->tableau_->setTVDCoeff(2.0);
  }
};


// ------------------------------------------------------------------------
template<class Scalar>
Teuchos::RCP<StepperEDIRK_TrapezoidalRule<Scalar> >
createStepperEDIRK_TrapezoidalRule(
  const Teuchos::RCP<const Thyra::ModelEvaluator<Scalar> >& model,
  Teuchos::RCP<Teuchos::ParameterList> pl)
{
  auto stepper = Teuchos::rcp(new StepperEDIRK_TrapezoidalRule<Scalar>());

  // Test for aliases.
  if (pl != Teuchos::null) {
    auto stepperType =
      pl->get<std::string>("Stepper Type", stepper->getStepperType());

    TEUCHOS_TEST_FOR_EXCEPTION(
      stepperType != stepper->getStepperType() &&
      stepperType != "RK Crank-Nicolson", std::logic_error,
      "  ParameterList 'Stepper Type' (='" + stepperType +"')\n"
      "  does not match type for this Stepper (='" + stepper->getStepperType() +
      "')\n  or one of its aliases ('RK Crank-Nicolson').\n");

    // Reset default StepperType.
    pl->set<std::string>("Stepper Type", stepper->getStepperType());
  }

  stepper->setStepperDIRKValues(pl);

  if (model != Teuchos::null) {
    stepper->setModel(model);
    stepper->initialize();
  }

  return stepper;
}


// ----------------------------------------------------------------------------
template<class Scalar>
class StepperSDIRK_ImplicitMidpoint :
  virtual public StepperDIRK<Scalar>
{
public:
  StepperSDIRK_ImplicitMidpoint()
  {
    this->setStepperName("RK Implicit Midpoint");
    this->setStepperType("RK Implicit Midpoint");
    this->setupTableau();
    this->setupDefault();
    this->setUseFSAL(            false);
    this->setICConsistency(      "None");
    this->setICConsistencyCheck( false);
  }

  StepperSDIRK_ImplicitMidpoint(
    const Teuchos::RCP<const Thyra::ModelEvaluator<Scalar> >& appModel,
    const Teuchos::RCP<Thyra::NonlinearSolverBase<Scalar> >& solver,
    bool useFSAL,
    std::string ICConsistency,
    bool ICConsistencyCheck,
    bool useEmbedded,
    bool zeroInitialGuess,
    const Teuchos::RCP<StepperRKAppAction<Scalar> >& stepperRKAppAction)
  {
    this->setStepperName("RK Implicit Midpoint");
    this->setStepperType("RK Implicit Midpoint");
    this->setupTableau();
    this->setup(appModel, solver, useFSAL, ICConsistency, ICConsistencyCheck,
                useEmbedded, zeroInitialGuess, stepperRKAppAction);
  }

  std::string getDescription() const
  {
    std::ostringstream Description;
    Description << this->getStepperType() << "\n"
                << "A-stable\n"
                << "Solving Ordinary Differential Equations II:\n"
                << "Stiff and Differential-Algebraic Problems,\n"
                << "2nd Revised Edition\n"
                << "E. Hairer and G. Wanner\n"
                << "Table 5.2, pg 72\n"
                << "Solving Ordinary Differential Equations I:\n"
                << "Nonstiff Problems, 2nd Revised Edition\n"
                << "E. Hairer, S. P. Norsett, and G. Wanner\n"
                << "Table 7.1, pg 205\n"
                << "c = [ 1/2 ]'\n"
                << "A = [ 1/2 ]\n"
                << "b = [  1  ]'";
    return Description.str();
  }

protected:

  void setupTableau()
  {
    typedef Teuchos::ScalarTraits<Scalar> ST;
    int NumStages = 1;
    Teuchos::SerialDenseMatrix<int,Scalar> A(NumStages,NumStages);
    Teuchos::SerialDenseVector<int,Scalar> b(NumStages);
    Teuchos::SerialDenseVector<int,Scalar> c(NumStages);
    const Scalar onehalf = ST::one()/(2*ST::one());
    const Scalar one = ST::one();

    // Fill A:
    A(0,0) = onehalf;

    // Fill b:
    b(0) = one;

    // Fill c:
    c(0) = onehalf;

    int order = 2;

    this->tableau_ = Teuchos::rcp(new RKButcherTableau<Scalar>(
      this->getStepperType(),A,b,c,order,order,order));
    this->tableau_->setTVD(true);
    this->tableau_->setTVDCoeff(2.0);
  }
};


// ------------------------------------------------------------------------
template<class Scalar>
Teuchos::RCP<StepperSDIRK_ImplicitMidpoint<Scalar> >
createStepperSDIRK_ImplicitMidpoint(
  const Teuchos::RCP<const Thyra::ModelEvaluator<Scalar> >& model,
  Teuchos::RCP<Teuchos::ParameterList> pl)
{
  auto stepper = Teuchos::rcp(new StepperSDIRK_ImplicitMidpoint<Scalar>());
  stepper->setStepperDIRKValues(pl);

  if (model != Teuchos::null) {
    stepper->setModel(model);
    stepper->initialize();
  }

  return stepper;
}


// ----------------------------------------------------------------------------
template<class Scalar>
class StepperSDIRK_SSPDIRK22 :
  virtual public StepperDIRK<Scalar>
{
  public:
  StepperSDIRK_SSPDIRK22()
  {
    this->setStepperName("SSPDIRK22");
    this->setStepperType("SSPDIRK22");
    this->setupTableau();
    this->setupDefault();
    this->setUseFSAL(            false);
    this->setICConsistency(      "None");
    this->setICConsistencyCheck( false);
  }

  StepperSDIRK_SSPDIRK22(
    const Teuchos::RCP<const Thyra::ModelEvaluator<Scalar> >& appModel,
    const Teuchos::RCP<Thyra::NonlinearSolverBase<Scalar> >& solver,
    bool useFSAL,
    std::string ICConsistency,
    bool ICConsistencyCheck,
    bool useEmbedded,
    bool zeroInitialGuess,
    const Teuchos::RCP<StepperRKAppAction<Scalar> >& stepperRKAppAction)
  {
    this->setStepperName("SSPDIRK22");
    this->setStepperType("SSPDIRK22");
    this->setupTableau();
    this->setup(appModel, solver, useFSAL, ICConsistency, ICConsistencyCheck,
                useEmbedded, zeroInitialGuess, stepperRKAppAction);
  }

  std::string getDescription() const
  {
    std::ostringstream Description;
    Description << this->getStepperType() << "\n"
      << "Strong Stability Preserving Diagonally-Implicit RK (stage=2, order=2)\n"
      << "SSP-Coef = 4\n"
      << "c =     [ 1/4   3/4 ]'\n"
      << "A =     [ 1/4       ]\n"
      << "        [ 1/2   1/4 ]\n"
      << "b     = [ 1/2   1/2 ]\n" << std::endl;
    return Description.str();
  }

protected:

  void setupTableau()
  {
    typedef Teuchos::ScalarTraits<Scalar> ST;
    using Teuchos::as;
    const int NumStages = 2;
    const int order     = 2;
    Teuchos::SerialDenseMatrix<int,Scalar> A(NumStages,NumStages);
    Teuchos::SerialDenseVector<int,Scalar> b(NumStages);
    Teuchos::SerialDenseVector<int,Scalar> c(NumStages);

    const Scalar one       = ST::one();
    const Scalar zero      = ST::zero();
    const Scalar onehalf   = one/(2*one);
    const Scalar onefourth = one/(4*one);

    // Fill A:
    A(0,0) = A(1,1) = onefourth;
    A(0,1) = zero;
    A(1,0) = onehalf;

    // Fill b:
    b(0) = b(1) = onehalf;

    // Fill c:
    c(0) = A(0,0);
    c(1) = A(1,0) + A(1,1);

    this->tableau_ = Teuchos::rcp(new RKButcherTableau<Scalar>(
      this->getStepperType(),A,b,c,order,order,order));
    this->tableau_->setTVD(true);
    this->tableau_->setTVDCoeff(4.0);
  }
};


// ------------------------------------------------------------------------
template<class Scalar>
Teuchos::RCP<StepperSDIRK_SSPDIRK22<Scalar> >
createStepperSDIRK_SSPDIRK22(
  const Teuchos::RCP<const Thyra::ModelEvaluator<Scalar> >& model,
  Teuchos::RCP<Teuchos::ParameterList> pl)
{
  auto stepper = Teuchos::rcp(new StepperSDIRK_SSPDIRK22<Scalar>());
  stepper->setStepperDIRKValues(pl);

  if (model != Teuchos::null) {
    stepper->setModel(model);
    stepper->initialize();
  }

  return stepper;
}


// ----------------------------------------------------------------------------
template<class Scalar>
class StepperSDIRK_SSPDIRK32 :
  virtual public StepperDIRK<Scalar>
{
  public:
  StepperSDIRK_SSPDIRK32()
  {
    this->setStepperName("SSPDIRK32");
    this->setStepperType("SSPDIRK32");
    this->setupTableau();
    this->setupDefault();
    this->setUseFSAL(            false);
    this->setICConsistency(      "None");
    this->setICConsistencyCheck( false);
  }

  StepperSDIRK_SSPDIRK32(
    const Teuchos::RCP<const Thyra::ModelEvaluator<Scalar> >& appModel,
    const Teuchos::RCP<Thyra::NonlinearSolverBase<Scalar> >& solver,
    bool useFSAL,
    std::string ICConsistency,
    bool ICConsistencyCheck,
    bool useEmbedded,
    bool zeroInitialGuess,
    const Teuchos::RCP<StepperRKAppAction<Scalar> >& stepperRKAppAction)
  {
    this->setStepperName("SSPDIRK32");
    this->setStepperType("SSPDIRK32");
    this->setupTableau();
    this->setup(appModel, solver, useFSAL, ICConsistency, ICConsistencyCheck,
                useEmbedded, zeroInitialGuess, stepperRKAppAction);
  }

  std::string getDescription() const
  {
    std::ostringstream Description;
    Description << this->getStepperType() << "\n"
                << "Strong Stability Preserving Diagonally-Implicit RK (stage=3, order=2)\n"
                << "SSP-Coef = 6\n"
                << "c = [ 1/6   1/2   5/6 ]'\n"
                << "A = [ 1/6             ]\n"
                << "    [ 1/3   1/6       ]\n"
                << "    [ 1/3   1/3   1/6 ]\n"
                << "b = [ 1/3   1/3   1/3 ]\n" << std::endl;
    return Description.str();
  }

protected:

  void setupTableau()
  {

    typedef Teuchos::ScalarTraits<Scalar> ST;
    using Teuchos::as;
    const int NumStages = 3;
    const int order     = 2;
    Teuchos::SerialDenseMatrix<int,Scalar> A(NumStages,NumStages);
    Teuchos::SerialDenseVector<int,Scalar> b(NumStages);
    Teuchos::SerialDenseVector<int,Scalar> c(NumStages);

    const Scalar one      = ST::one();
    const Scalar zero     = ST::zero();
    const Scalar onethird = one/(3*one);
    const Scalar onesixth = one/(6*one);

    // Fill A:
    A(0,0) = A(1,1) = A(2,2) = onesixth;
    A(1,0) = A(2,0) = A(2,1) = onethird;
    A(0,1) = A(0,2) = A(1,2) = zero;

    // Fill b:
    b(0) = b(1) = b(2) = onethird;

    // Fill c:
    c(0) = A(0,0);
    c(1) = A(1,0) + A(1,1);
    c(2) = A(2,0) + A(2,1) + A(2,2);

    this->tableau_ = Teuchos::rcp(new RKButcherTableau<Scalar>(
      this->getStepperType(),A,b,c,order,order,order));
    this->tableau_->setTVD(true);
    this->tableau_->setTVDCoeff(6.0);
  }
};


// ------------------------------------------------------------------------
template<class Scalar>
Teuchos::RCP<StepperSDIRK_SSPDIRK32<Scalar> >
createStepperSDIRK_SSPDIRK32(
  const Teuchos::RCP<const Thyra::ModelEvaluator<Scalar> >& model,
  Teuchos::RCP<Teuchos::ParameterList> pl)
{
  auto stepper = Teuchos::rcp(new StepperSDIRK_SSPDIRK32<Scalar>());
  stepper->setStepperDIRKValues(pl);

  if (model != Teuchos::null) {
    stepper->setModel(model);
    stepper->initialize();
  }

  return stepper;
}


// ----------------------------------------------------------------------------
template<class Scalar>
class StepperSDIRK_SSPDIRK23 :
  virtual public StepperDIRK<Scalar>
{
  public:
  StepperSDIRK_SSPDIRK23()
  {
    this->setStepperName("SSPDIRK23");
    this->setStepperType("SSPDIRK23");
    this->setupTableau();
    this->setupDefault();
    this->setUseFSAL(            false);
    this->setICConsistency(      "None");
    this->setICConsistencyCheck( false);
  }

  StepperSDIRK_SSPDIRK23(
    const Teuchos::RCP<const Thyra::ModelEvaluator<Scalar> >& appModel,
    const Teuchos::RCP<Thyra::NonlinearSolverBase<Scalar> >& solver,
    bool useFSAL,
    std::string ICConsistency,
    bool ICConsistencyCheck,
    bool useEmbedded,
    bool zeroInitialGuess,
    const Teuchos::RCP<StepperRKAppAction<Scalar> >& stepperRKAppAction)
  {
    this->setStepperName("SSPDIRK23");
    this->setStepperType("SSPDIRK23");
    this->setupTableau();
    this->setup(appModel, solver, useFSAL, ICConsistency, ICConsistencyCheck,
                useEmbedded, zeroInitialGuess, stepperRKAppAction);
  }

  std::string getDescription() const
  {
    std::ostringstream Description;
    Description << this->getStepperType() << "\n"
      << "Strong Stability Preserving Diagonally-Implicit RK (stage=2, order=3)\n"
      << "SSP-Coef = 1 + sqrt( 3 )\n"
      << "c =     [ 1/(3 + sqrt( 3 ))  (1/6)(3 + sqrt( 3 )) ] '\n"
      << "A =     [ 1/(3 + sqrt( 3 ))                       ] \n"
      << "        [ 1/sqrt( 3 )        1/(3 + sqrt( 3 ))    ] \n"
      << "b     = [ 1/2                   1/2               ] \n" << std::endl;
    return Description.str();
  }

protected:

  void setupTableau()
  {

    typedef Teuchos::ScalarTraits<Scalar> ST;
    using Teuchos::as;
    const int NumStages = 2;
    const int order     = 3;
    const Scalar sspcoef = 2.7321;
    Teuchos::SerialDenseMatrix<int,Scalar> A(NumStages,NumStages);
    Teuchos::SerialDenseVector<int,Scalar> b(NumStages);
    Teuchos::SerialDenseVector<int,Scalar> c(NumStages);

    const Scalar one       = ST::one();
    const Scalar zero      = ST::zero();
    const Scalar onehalf   = one/(2*one);
    const Scalar rootthree = ST::squareroot(3*one);

    // Fill A:
    A(0,0) = A(1,1) = one/(3*one + rootthree);
    A(1,0) = one/rootthree;
    A(0,1) = zero;

    // Fill b:
    b(0) = b(1) = onehalf;

    // Fill c:
    c(0) = A(0,0);
    c(1) = A(1,0) + A(1,1);

    this->tableau_ = Teuchos::rcp(new RKButcherTableau<Scalar>(
      this->getStepperType(),A,b,c,order,order,order));
    this->tableau_->setTVD(true);
    this->tableau_->setTVDCoeff(sspcoef);
  }
};


// ------------------------------------------------------------------------
template<class Scalar>
Teuchos::RCP<StepperSDIRK_SSPDIRK23<Scalar> >
createStepperSDIRK_SSPDIRK23(
  const Teuchos::RCP<const Thyra::ModelEvaluator<Scalar> >& model,
  Teuchos::RCP<Teuchos::ParameterList> pl)
{
  auto stepper = Teuchos::rcp(new StepperSDIRK_SSPDIRK23<Scalar>());
  stepper->setStepperDIRKValues(pl);

  if (model != Teuchos::null) {
    stepper->setModel(model);
    stepper->initialize();
  }

  return stepper;
}


// ----------------------------------------------------------------------------
template<class Scalar>
class StepperSDIRK_SSPDIRK33 :
  virtual public StepperDIRK<Scalar>
{
  public:
  StepperSDIRK_SSPDIRK33()
  {
    this->setStepperName("SSPDIRK33");
    this->setStepperType("SSPDIRK33");
    this->setupTableau();
    this->setupDefault();
    this->setUseFSAL(            false);
    this->setICConsistency(      "None");
    this->setICConsistencyCheck( false);
  }

  StepperSDIRK_SSPDIRK33(
    const Teuchos::RCP<const Thyra::ModelEvaluator<Scalar> >& appModel,
    const Teuchos::RCP<Thyra::NonlinearSolverBase<Scalar> >& solver,
    bool useFSAL,
    std::string ICConsistency,
    bool ICConsistencyCheck,
    bool useEmbedded,
    bool zeroInitialGuess,
    const Teuchos::RCP<StepperRKAppAction<Scalar> >& stepperRKAppAction)
  {
    this->setStepperName("SSPDIRK33");
    this->setStepperType("SSPDIRK33");
    this->setupTableau();
    this->setup(appModel, solver, useFSAL, ICConsistency, ICConsistencyCheck,
                useEmbedded, zeroInitialGuess, stepperRKAppAction);
  }

  std::string getDescription() const
  {
    std::ostringstream Description;
    Description << this->getStepperType() << "\n"
      << "Strong Stability Preserving Diagonally-Implicit RK (stage=3, order=3)\n"
      << "SSP-Coef = 2 + 2 sqrt(2)\n"
      << "c =     [ 1/( 4 + 2 sqrt(2)      1/2            (1/4)(2 + sqrt(2) ] '\n"
      << "A =     [ 1/( 4 + 2 sqrt(2)                                       ] \n"
      << "        [ 1/(2 sqrt(2)       1/( 4 + 2 sqrt(2)                    ] \n"
      << "        [ 1/(2 sqrt(2)        1/(2 sqrt(2)      1/( 4 + 2 sqrt(2) ] \n"
      << "b     = [ 1/3                    1/3                1/3           ] \n"
      << std::endl;
    return Description.str();
  }

protected:

  void setupTableau()
  {

    typedef Teuchos::ScalarTraits<Scalar> ST;
    using Teuchos::as;
    const int NumStages = 3;
    const int order     = 3;
    const Scalar sspcoef= 4.8284;
    Teuchos::SerialDenseMatrix<int,Scalar> A(NumStages,NumStages);
    Teuchos::SerialDenseVector<int,Scalar> b(NumStages);
    Teuchos::SerialDenseVector<int,Scalar> c(NumStages);

    const Scalar one      = ST::one();
    const Scalar zero     = ST::zero();
    const Scalar onethird = one/(3*one);
    const Scalar rootwo   = ST::squareroot(2*one);

    // Fill A:
    A(0,0) = A(1,1) = A(2,2) = one / (4*one + 2*rootwo);
    A(1,0) = A(2,0) = A(2,1) = one / (2*rootwo);
    A(0,1) = A(0,2) = A(1,2) = zero;

    // Fill b:
    b(0) = b(1) = b(2) = onethird;

    // Fill c:
    c(0) = A(0,0);
    c(1) = A(1,0) + A(1,1);
    c(2) = A(2,0) + A(2,1) + A(2,2);

    this->tableau_ = Teuchos::rcp(new RKButcherTableau<Scalar>(
      this->getStepperType(),A,b,c,order,order,order));
    this->tableau_->setTVD(true);
    this->tableau_->setTVDCoeff(sspcoef);
  }
};


// ------------------------------------------------------------------------
template<class Scalar>
Teuchos::RCP<StepperSDIRK_SSPDIRK33<Scalar> >
createStepperSDIRK_SSPDIRK33(
  const Teuchos::RCP<const Thyra::ModelEvaluator<Scalar> >& model,
  Teuchos::RCP<Teuchos::ParameterList> pl)
{
  auto stepper = Teuchos::rcp(new StepperSDIRK_SSPDIRK33<Scalar>());
  stepper->setStepperDIRKValues(pl);

  if (model != Teuchos::null) {
    stepper->setModel(model);
    stepper->initialize();
  }

  return stepper;
}


// ----------------------------------------------------------------------------
template<class Scalar>
class StepperDIRK_1Stage1stOrderRadauIA :
  virtual public StepperDIRK<Scalar>
{
public:
  StepperDIRK_1Stage1stOrderRadauIA()
  {
    this->setStepperName("RK Implicit 1 Stage 1st order Radau IA");
    this->setStepperType("RK Implicit 1 Stage 1st order Radau IA");
    this->setupTableau();
    this->setupDefault();
    this->setUseFSAL(            false);
    this->setICConsistency(      "None");
    this->setICConsistencyCheck( false);
  }

  StepperDIRK_1Stage1stOrderRadauIA(
    const Teuchos::RCP<const Thyra::ModelEvaluator<Scalar> >& appModel,
    const Teuchos::RCP<Thyra::NonlinearSolverBase<Scalar> >& solver,
    bool useFSAL,
    std::string ICConsistency,
    bool ICConsistencyCheck,
    bool useEmbedded,
    bool zeroInitialGuess,
    const Teuchos::RCP<StepperRKAppAction<Scalar> >& stepperRKAppAction)
  {
    this->setStepperName("RK Implicit 1 Stage 1st order Radau IA");
    this->setStepperType("RK Implicit 1 Stage 1st order Radau IA");
    this->setupTableau();
    this->setup(appModel, solver, useFSAL, ICConsistency, ICConsistencyCheck,
                useEmbedded, zeroInitialGuess, stepperRKAppAction);
  }

  std::string getDescription() const
  {
    std::ostringstream Description;
    Description << this->getStepperType() << "\n"
                << "A-stable\n"
                << "Solving Ordinary Differential Equations II:\n"
                << "Stiff and Differential-Algebraic Problems,\n"
                << "2nd Revised Edition\n"
                << "E. Hairer and G. Wanner\n"
                << "Table 5.3, pg 73\n"
                << "c = [ 0 ]'\n"
                << "A = [ 1 ]\n"
                << "b = [ 1 ]'";
    return Description.str();
  }

protected:

  void setupTableau()
  {
    typedef Teuchos::ScalarTraits<Scalar> ST;
    int NumStages = 1;
    Teuchos::SerialDenseMatrix<int,Scalar> A(NumStages,NumStages);
    Teuchos::SerialDenseVector<int,Scalar> b(NumStages);
    Teuchos::SerialDenseVector<int,Scalar> c(NumStages);
    const Scalar one = ST::one();
    const Scalar zero = ST::zero();
    A(0,0) = one;
    b(0) = one;
    c(0) = zero;
    int order = 1;

    auto emptyBStar = Teuchos::SerialDenseVector<int,Scalar>();
    this->tableau_ = Teuchos::rcp(new RKButcherTableau<Scalar>(
      this->getStepperType(),A,b,c,order,order,order,emptyBStar,false));
  }
};


// ------------------------------------------------------------------------
template<class Scalar>
Teuchos::RCP<StepperDIRK_1Stage1stOrderRadauIA<Scalar> >
createStepperDIRK_1Stage1stOrderRadauIA(
  const Teuchos::RCP<const Thyra::ModelEvaluator<Scalar> >& model,
  Teuchos::RCP<Teuchos::ParameterList> pl)
{
  auto stepper = Teuchos::rcp(new StepperDIRK_1Stage1stOrderRadauIA<Scalar>());
  stepper->setStepperDIRKValues(pl);

  if (model != Teuchos::null) {
    stepper->setModel(model);
    stepper->initialize();
  }

  return stepper;
}


// ----------------------------------------------------------------------------
template<class Scalar>
class StepperDIRK_2Stage2ndOrderLobattoIIIB :
  virtual public StepperDIRK<Scalar>
{
public:
  StepperDIRK_2Stage2ndOrderLobattoIIIB()
  {
    this->setStepperName("RK Implicit 2 Stage 2nd order Lobatto IIIB");
    this->setStepperType("RK Implicit 2 Stage 2nd order Lobatto IIIB");
    this->setupTableau();
    this->setupDefault();
    this->setUseFSAL(            false);
    this->setICConsistency(      "None");
    this->setICConsistencyCheck( false);
  }

  StepperDIRK_2Stage2ndOrderLobattoIIIB(
    const Teuchos::RCP<const Thyra::ModelEvaluator<Scalar> >& appModel,
    const Teuchos::RCP<Thyra::NonlinearSolverBase<Scalar> >& solver,
    bool useFSAL,
    std::string ICConsistency,
    bool ICConsistencyCheck,
    bool useEmbedded,
    bool zeroInitialGuess,
    const Teuchos::RCP<StepperRKAppAction<Scalar> >& stepperRKAppAction)
  {
    this->setStepperName("RK Implicit 2 Stage 2nd order Lobatto IIIB");
    this->setStepperType("RK Implicit 2 Stage 2nd order Lobatto IIIB");
    this->setupTableau();
    this->setup(appModel, solver, useFSAL, ICConsistency, ICConsistencyCheck,
                useEmbedded, zeroInitialGuess, stepperRKAppAction);
  }

  std::string getDescription() const
  {
    std::ostringstream Description;
    Description << this->getStepperType() << "\n"
                << "A-stable\n"
                << "Solving Ordinary Differential Equations II:\n"
                << "Stiff and Differential-Algebraic Problems,\n"
                << "2nd Revised Edition\n"
                << "E. Hairer and G. Wanner\n"
                << "Table 5.9, pg 76\n"
                << "c = [  0    1   ]'\n"
                << "A = [ 1/2   0   ]\n"
                << "    [ 1/2   0   ]\n"
                << "b = [ 1/2  1/2  ]'";
    return Description.str();
  }

protected:

  void setupTableau()
  {
    typedef Teuchos::ScalarTraits<Scalar> ST;
    using Teuchos::as;
    int NumStages = 2;
    Teuchos::SerialDenseMatrix<int,Scalar> A(NumStages,NumStages);
    Teuchos::SerialDenseVector<int,Scalar> b(NumStages);
    Teuchos::SerialDenseVector<int,Scalar> c(NumStages);
    const Scalar zero = ST::zero();
    const Scalar one = ST::one();

    // Fill A:
    A(0,0) = as<Scalar>( one/(2*one) );
    A(0,1) = zero;
    A(1,0) = as<Scalar>( one/(2*one) );
    A(1,1) = zero;

    // Fill b:
    b(0) = as<Scalar>( one/(2*one) );
    b(1) = as<Scalar>( one/(2*one) );

    // Fill c:
    c(0) = zero;
    c(1) = one;
    int order = 2;

    auto emptyBStar = Teuchos::SerialDenseVector<int,Scalar>();
    this->tableau_ = Teuchos::rcp(new RKButcherTableau<Scalar>(
      this->getStepperType(),A,b,c,order,order,order,emptyBStar,false));
    this->tableau_->setTVD(true);
    this->tableau_->setTVDCoeff(2.0);
    //TODO: fix this
  }

};


// ------------------------------------------------------------------------
template<class Scalar>
Teuchos::RCP<StepperDIRK_2Stage2ndOrderLobattoIIIB<Scalar> >
createStepperDIRK_2Stage2ndOrderLobattoIIIB(
  const Teuchos::RCP<const Thyra::ModelEvaluator<Scalar> >& model,
  Teuchos::RCP<Teuchos::ParameterList> pl)
{
  auto stepper = Teuchos::rcp(new StepperDIRK_2Stage2ndOrderLobattoIIIB<Scalar>());
  stepper->setStepperDIRKValues(pl);

  if (model != Teuchos::null) {
    stepper->setModel(model);
    stepper->initialize();
  }

  return stepper;
}


// ----------------------------------------------------------------------------
template<class Scalar>
class StepperSDIRK_5Stage4thOrder :
  virtual public StepperDIRK<Scalar>
{
public:
  StepperSDIRK_5Stage4thOrder()
  {
    this->setStepperName("SDIRK 5 Stage 4th order");
    this->setStepperType("SDIRK 5 Stage 4th order");
    this->setupTableau();
    this->setupDefault();
    this->setUseFSAL(            false);
    this->setICConsistency(      "None");
    this->setICConsistencyCheck( false);
  }

  StepperSDIRK_5Stage4thOrder(
    const Teuchos::RCP<const Thyra::ModelEvaluator<Scalar> >& appModel,
    const Teuchos::RCP<Thyra::NonlinearSolverBase<Scalar> >& solver,
    bool useFSAL,
    std::string ICConsistency,
    bool ICConsistencyCheck,
    bool useEmbedded,
    bool zeroInitialGuess,
    const Teuchos::RCP<StepperRKAppAction<Scalar> >& stepperRKAppAction)
  {
    this->setStepperName("SDIRK 5 Stage 4th order");
    this->setStepperType("SDIRK 5 Stage 4th order");
    this->setupTableau();
    this->setup(appModel, solver, useFSAL, ICConsistency, ICConsistencyCheck,
                useEmbedded, zeroInitialGuess, stepperRKAppAction);
  }

  std::string getDescription() const
  {
    std::ostringstream Description;
    Description << this->getStepperType() << "\n"
      << "L-stable\n"
      << "Solving Ordinary Differential Equations II:\n"
      << "Stiff and Differential-Algebraic Problems,\n"
      << "2nd Revised Edition\n"
      << "E. Hairer and G. Wanner\n"
      << "pg100 \n"
      << "c     = [ 1/4       3/4        11/20   1/2     1   ]'\n"
      << "A     = [ 1/4                                      ]\n"
      << "        [ 1/2       1/4                            ]\n"
      << "        [ 17/50     -1/25      1/4                 ]\n"
      << "        [ 371/1360  -137/2720  15/544  1/4         ]\n"
      << "        [ 25/24     -49/48     125/16  -85/12  1/4 ]\n"
      << "b     = [ 25/24     -49/48     125/16  -85/12  1/4 ]'";
      // << "b     = [ 59/48     -17/96     225/32  -85/12  0   ]'";
    return Description.str();
  }

protected:

  void setupTableau()
  {
    typedef Teuchos::ScalarTraits<Scalar> ST;
    using Teuchos::as;
    int NumStages = 5;
    Teuchos::SerialDenseMatrix<int,Scalar> A(NumStages,NumStages);
    Teuchos::SerialDenseVector<int,Scalar> b(NumStages);
    Teuchos::SerialDenseVector<int,Scalar> c(NumStages);
    const Scalar zero = ST::zero();
    const Scalar one = ST::one();
    const Scalar onequarter = as<Scalar>( one/(4*one) );

    // Fill A:
    A(0,0) = onequarter;
    A(0,1) = zero;
    A(0,2) = zero;
    A(0,3) = zero;
    A(0,4) = zero;

    A(1,0) = as<Scalar>( one / (2*one) );
    A(1,1) = onequarter;
    A(1,2) = zero;
    A(1,3) = zero;
    A(1,4) = zero;

    A(2,0) = as<Scalar>( 17*one/(50*one) );
    A(2,1) = as<Scalar>( -one/(25*one) );
    A(2,2) = onequarter;
    A(2,3) = zero;
    A(2,4) = zero;

    A(3,0) = as<Scalar>( 371*one/(1360*one) );
    A(3,1) = as<Scalar>( -137*one/(2720*one) );
    A(3,2) = as<Scalar>( 15*one/(544*one) );
    A(3,3) = onequarter;
    A(3,4) = zero;

    A(4,0) = as<Scalar>( 25*one/(24*one) );
    A(4,1) = as<Scalar>( -49*one/(48*one) );
    A(4,2) = as<Scalar>( 125*one/(16*one) );
    A(4,3) = as<Scalar>( -85*one/(12*one) );
    A(4,4) = onequarter;

    // Fill b:
    b(0) = as<Scalar>( 25*one/(24*one) );
    b(1) = as<Scalar>( -49*one/(48*one) );
    b(2) = as<Scalar>( 125*one/(16*one) );
    b(3) = as<Scalar>( -85*one/(12*one) );
    b(4) = onequarter;

    /*
    // Alternate version
    b(0) = as<Scalar>( 59*one/(48*one) );
    b(1) = as<Scalar>( -17*one/(96*one) );
    b(2) = as<Scalar>( 225*one/(32*one) );
    b(3) = as<Scalar>( -85*one/(12*one) );
    b(4) = zero;
    */

    // Fill c:
    c(0) = onequarter;
    c(1) = as<Scalar>( 3*one/(4*one) );
    c(2) = as<Scalar>( 11*one/(20*one) );
    c(3) = as<Scalar>( one/(2*one) );
    c(4) = one;

    int order = 4;

    this->tableau_ = Teuchos::rcp(new RKButcherTableau<Scalar>(
      this->getStepperType(),A,b,c,order,order,order));
  }
};


// ------------------------------------------------------------------------
template<class Scalar>
Teuchos::RCP<StepperSDIRK_5Stage4thOrder<Scalar> >
createStepperSDIRK_5Stage4thOrder(
  const Teuchos::RCP<const Thyra::ModelEvaluator<Scalar> >& model,
  Teuchos::RCP<Teuchos::ParameterList> pl)
{
  auto stepper = Teuchos::rcp(new StepperSDIRK_5Stage4thOrder<Scalar>());
  stepper->setStepperDIRKValues(pl);

  if (model != Teuchos::null) {
    stepper->setModel(model);
    stepper->initialize();
  }

  return stepper;
}


// ----------------------------------------------------------------------------
template<class Scalar>
class StepperSDIRK_3Stage4thOrder :
  virtual public StepperDIRK<Scalar>
{
public:
  StepperSDIRK_3Stage4thOrder()
  {
    this->setStepperName("SDIRK 3 Stage 4th order");
    this->setStepperType("SDIRK 3 Stage 4th order");
    this->setupTableau();
    this->setupDefault();
    this->setUseFSAL(            false);
    this->setICConsistency(      "None");
    this->setICConsistencyCheck( false);
  }

  StepperSDIRK_3Stage4thOrder(
    const Teuchos::RCP<const Thyra::ModelEvaluator<Scalar> >& appModel,
    const Teuchos::RCP<Thyra::NonlinearSolverBase<Scalar> >& solver,
    bool useFSAL,
    std::string ICConsistency,
    bool ICConsistencyCheck,
    bool useEmbedded,
    bool zeroInitialGuess,
    const Teuchos::RCP<StepperRKAppAction<Scalar> >& stepperRKAppAction)
  {
    this->setStepperName("SDIRK 3 Stage 4th order");
    this->setStepperType("SDIRK 3 Stage 4th order");
    this->setupTableau();
    this->setup(appModel, solver, useFSAL, ICConsistency, ICConsistencyCheck,
                useEmbedded, zeroInitialGuess, stepperRKAppAction);
  }

  std::string getDescription() const
  {
    std::ostringstream Description;
    Description << this->getStepperType() << "\n"
                << "A-stable\n"
                << "Solving Ordinary Differential Equations II:\n"
                << "Stiff and Differential-Algebraic Problems,\n"
                << "2nd Revised Edition\n"
                << "E. Hairer and G. Wanner\n"
                << "p. 100 \n"
                << "gamma = (1/sqrt(3))*cos(pi/18)+1/2\n"
                << "delta = 1/(6*(2*gamma-1)^2)\n"
                << "c = [ gamma      1/2        1-gamma ]'\n"
                << "A = [ gamma                         ]\n"
                << "    [ 1/2-gamma  gamma              ]\n"
                << "    [ 2*gamma    1-4*gamma  gamma   ]\n"
                << "b = [ delta      1-2*delta  delta   ]'";
    return Description.str();
  }

protected:

  void setupTableau()
  {
    typedef Teuchos::ScalarTraits<Scalar> ST;
    using Teuchos::as;
    int NumStages = 3;
    Teuchos::SerialDenseMatrix<int,Scalar> A(NumStages,NumStages);
    Teuchos::SerialDenseVector<int,Scalar> b(NumStages);
    Teuchos::SerialDenseVector<int,Scalar> c(NumStages);
    const Scalar zero = ST::zero();
    const Scalar one = ST::one();
    const Scalar pi = as<Scalar>(4*one)*std::atan(one);
    const Scalar gamma = as<Scalar>( one/ST::squareroot(3*one)*std::cos(pi/(18*one))+one/(2*one) );
    const Scalar delta = as<Scalar>( one/(6*one*std::pow(2*gamma-one,2*one)) );

    // Fill A:
    A(0,0) = gamma;
    A(0,1) = zero;
    A(0,2) = zero;

    A(1,0) = as<Scalar>( one/(2*one) - gamma );
    A(1,1) = gamma;
    A(1,2) = zero;

    A(2,0) = as<Scalar>( 2*gamma );
    A(2,1) = as<Scalar>( one - 4*gamma );
    A(2,2) = gamma;

    // Fill b:
    b(0) = delta;
    b(1) = as<Scalar>( one-2*delta );
    b(2) = delta;

    // Fill c:
    c(0) = gamma;
    c(1) = as<Scalar>( one/(2*one) );
    c(2) = as<Scalar>( one - gamma );

    int order = 4;

    this->tableau_ = Teuchos::rcp(new RKButcherTableau<Scalar>(
      this->getStepperType(),A,b,c,order,order,order));
  }
};


// ------------------------------------------------------------------------
template<class Scalar>
Teuchos::RCP<StepperSDIRK_3Stage4thOrder<Scalar> >
createStepperSDIRK_3Stage4thOrder(
  const Teuchos::RCP<const Thyra::ModelEvaluator<Scalar> >& model,
  Teuchos::RCP<Teuchos::ParameterList> pl)
{
  auto stepper = Teuchos::rcp(new StepperSDIRK_3Stage4thOrder<Scalar>());
  stepper->setStepperDIRKValues(pl);

  if (model != Teuchos::null) {
    stepper->setModel(model);
    stepper->initialize();
  }

  return stepper;
}


// ----------------------------------------------------------------------------
template<class Scalar>
class StepperSDIRK_5Stage5thOrder :
  virtual public StepperDIRK<Scalar>
{
public:
  StepperSDIRK_5Stage5thOrder()
  {
    this->setStepperName("SDIRK 5 Stage 5th order");
    this->setStepperType("SDIRK 5 Stage 5th order");
    this->setupTableau();
    this->setupDefault();
    this->setUseFSAL(            false);
    this->setICConsistency(      "None");
    this->setICConsistencyCheck( false);
  }

  StepperSDIRK_5Stage5thOrder(
    const Teuchos::RCP<const Thyra::ModelEvaluator<Scalar> >& appModel,
    const Teuchos::RCP<Thyra::NonlinearSolverBase<Scalar> >& solver,
    bool useFSAL,
    std::string ICConsistency,
    bool ICConsistencyCheck,
    bool useEmbedded,
    bool zeroInitialGuess,
    const Teuchos::RCP<StepperRKAppAction<Scalar> >& stepperRKAppAction)
  {
    this->setStepperName("SDIRK 5 Stage 5th order");
    this->setStepperType("SDIRK 5 Stage 5th order");
    this->setupTableau();
    this->setup(appModel, solver, useFSAL, ICConsistency, ICConsistencyCheck,
                useEmbedded, zeroInitialGuess, stepperRKAppAction);
  }

  std::string getDescription() const
  {
    std::ostringstream Description;
    Description << this->getStepperType() << "\n"
      << "Solving Ordinary Differential Equations II:\n"
      << "Stiff and Differential-Algebraic Problems,\n"
      << "2nd Revised Edition\n"
      << "E. Hairer and G. Wanner\n"
      << "pg101 \n"
      << "c = [ (6-sqrt(6))/10   ]\n"
      << "    [ (6+9*sqrt(6))/35 ]\n"
      << "    [ 1                ]\n"
      << "    [ (4-sqrt(6))/10   ]\n"
      << "    [ (4+sqrt(6))/10   ]\n"
      << "A = [ A1 A2 A3 A4 A5 ]\n"
      << "      A1 = [ (6-sqrt(6))/10               ]\n"
      << "           [ (-6+5*sqrt(6))/14            ]\n"
      << "           [ (888+607*sqrt(6))/2850       ]\n"
      << "           [ (3153-3082*sqrt(6))/14250    ]\n"
      << "           [ (-32583+14638*sqrt(6))/71250 ]\n"
      << "      A2 = [ 0                           ]\n"
      << "           [ (6-sqrt(6))/10              ]\n"
      << "           [ (126-161*sqrt(6))/1425      ]\n"
      << "           [ (3213+1148*sqrt(6))/28500   ]\n"
      << "           [ (-17199+364*sqrt(6))/142500 ]\n"
      << "      A3 = [ 0                       ]\n"
      << "           [ 0                       ]\n"
      << "           [ (6-sqrt(6))/10          ]\n"
      << "           [ (-267+88*sqrt(6))/500   ]\n"
      << "           [ (1329-544*sqrt(6))/2500 ]\n"
      << "      A4 = [ 0                     ]\n"
      << "           [ 0                     ]\n"
      << "           [ 0                     ]\n"
      << "           [ (6-sqrt(6))/10        ]\n"
      << "           [ (-96+131*sqrt(6))/625 ]\n"
      << "      A5 = [ 0              ]\n"
      << "           [ 0              ]\n"
      << "           [ 0              ]\n"
      << "           [ 0              ]\n"
      << "           [ (6-sqrt(6))/10 ]\n"
      << "b = [               0 ]\n"
      << "    [               0 ]\n"
      << "    [             1/9 ]\n"
      << "    [ (16-sqrt(6))/36 ]\n"
      << "    [ (16+sqrt(6))/36 ]'";
    return Description.str();
  }

protected:

  void setupTableau()
  {
    typedef Teuchos::ScalarTraits<Scalar> ST;
    using Teuchos::as;
    int NumStages = 5;
    Teuchos::SerialDenseMatrix<int,Scalar> A(NumStages,NumStages);
    Teuchos::SerialDenseVector<int,Scalar> b(NumStages);
    Teuchos::SerialDenseVector<int,Scalar> c(NumStages);
    const Scalar zero = ST::zero();
    const Scalar one = ST::one();
    const Scalar sqrt6 = ST::squareroot(as<Scalar>(6*one));
    const Scalar gamma = as<Scalar>( (6*one - sqrt6) / (10*one) ); // diagonal

    // Fill A:
    A(0,0) = gamma;
    A(0,1) = zero;
    A(0,2) = zero;
    A(0,3) = zero;
    A(0,4) = zero;

    A(1,0) = as<Scalar>( (-6*one+5*one*sqrt6)/(14*one) );
    A(1,1) = gamma;
    A(1,2) = zero;
    A(1,3) = zero;
    A(1,4) = zero;

    A(2,0) = as<Scalar>( (888*one+607*one*sqrt6)/(2850*one) );
    A(2,1) = as<Scalar>( (126*one-161*one*sqrt6)/(1425*one) );
    A(2,2) = gamma;
    A(2,3) = zero;
    A(2,4) = zero;

    A(3,0) = as<Scalar>( (3153*one-3082*one*sqrt6)/(14250*one) );
    A(3,1) = as<Scalar>( (3213*one+1148*one*sqrt6)/(28500*one) );
    A(3,2) = as<Scalar>( (-267*one+88*one*sqrt6)/(500*one) );
    A(3,3) = gamma;
    A(3,4) = zero;

    A(4,0) = as<Scalar>( (-32583*one+14638*one*sqrt6)/(71250*one) );
    A(4,1) = as<Scalar>( (-17199*one+364*one*sqrt6)/(142500*one) );
    A(4,2) = as<Scalar>( (1329*one-544*one*sqrt6)/(2500*one) );
    A(4,3) = as<Scalar>( (-96*one+131*sqrt6)/(625*one) );
    A(4,4) = gamma;

    // Fill b:
    b(0) = zero;
    b(1) = zero;
    b(2) = as<Scalar>( one/(9*one) );
    b(3) = as<Scalar>( (16*one-sqrt6)/(36*one) );
    b(4) = as<Scalar>( (16*one+sqrt6)/(36*one) );

    // Fill c:
    c(0) = gamma;
    c(1) = as<Scalar>( (6*one+9*one*sqrt6)/(35*one) );
    c(2) = one;
    c(3) = as<Scalar>( (4*one-sqrt6)/(10*one) );
    c(4) = as<Scalar>( (4*one+sqrt6)/(10*one) );

    int order = 5;

    this->tableau_ = Teuchos::rcp(new RKButcherTableau<Scalar>(
      this->getStepperType(),A,b,c,order,order,order));
  }
};


// ------------------------------------------------------------------------
template<class Scalar>
Teuchos::RCP<StepperSDIRK_5Stage5thOrder<Scalar> >
createStepperSDIRK_5Stage5thOrder(
  const Teuchos::RCP<const Thyra::ModelEvaluator<Scalar> >& model,
  Teuchos::RCP<Teuchos::ParameterList> pl)
{
  auto stepper = Teuchos::rcp(new StepperSDIRK_5Stage5thOrder<Scalar>());
  stepper->setStepperDIRKValues(pl);

  if (model != Teuchos::null) {
    stepper->setModel(model);
    stepper->initialize();
  }

  return stepper;
}


// ----------------------------------------------------------------------------
template<class Scalar>
class StepperSDIRK_21Pair :
  virtual public StepperDIRK<Scalar>
{
public:
  StepperSDIRK_21Pair()
  {
    this->setStepperName("SDIRK 2(1) Pair");
    this->setStepperType("SDIRK 2(1) Pair");
    this->setupTableau();
    this->setupDefault();
    this->setUseFSAL(            false);
    this->setICConsistency(      "None");
    this->setICConsistencyCheck( false);
  }

  StepperSDIRK_21Pair(
    const Teuchos::RCP<const Thyra::ModelEvaluator<Scalar> >& appModel,
    const Teuchos::RCP<Thyra::NonlinearSolverBase<Scalar> >& solver,
    bool useFSAL,
    std::string ICConsistency,
    bool ICConsistencyCheck,
    bool useEmbedded,
    bool zeroInitialGuess,
    const Teuchos::RCP<StepperRKAppAction<Scalar> >& stepperRKAppAction)
  {
    this->setStepperName("SDIRK 2(1) Pair");
    this->setStepperType("SDIRK 2(1) Pair");
    this->setupTableau();
    this->setup(appModel, solver, useFSAL, ICConsistency, ICConsistencyCheck,
                useEmbedded, zeroInitialGuess, stepperRKAppAction);
  }

  std::string getDescription() const
  {
    std::ostringstream Description;
    Description << this->getStepperType() << "\n"
                << "c =     [  1  0   ]'\n"
                << "A =     [  1      ]\n"
                << "        [ -1  1   ]\n"
                << "b     = [ 1/2 1/2 ]'\n"
                << "bstar = [  1  0   ]'";
    return Description.str();
  }

protected:

  void setupTableau()
  {
    typedef Teuchos::ScalarTraits<Scalar> ST;
    using Teuchos::as;
    int NumStages = 2;
    Teuchos::SerialDenseMatrix<int,Scalar> A(NumStages,NumStages);
    Teuchos::SerialDenseVector<int,Scalar> b(NumStages);
    Teuchos::SerialDenseVector<int,Scalar> c(NumStages);
    Teuchos::SerialDenseVector<int,Scalar> bstar(NumStages);

    const Scalar one = ST::one();
    const Scalar zero = ST::zero();

    // Fill A:
    A(0,0) =  one; A(0,1) = zero;
    A(1,0) = -one; A(1,1) =  one;

    // Fill b:
    b(0) = as<Scalar>(one/(2*one));
    b(1) = as<Scalar>(one/(2*one));

    // Fill c:
    c(0) = one;
    c(1) = zero;

    // Fill bstar
    bstar(0) = one;
    bstar(1) = zero;
    int order = 2;

    this->tableau_ = Teuchos::rcp(new RKButcherTableau<Scalar>(
      this->getStepperType(),A,b,c,order,order,order,bstar));
  }
};


// ------------------------------------------------------------------------
template<class Scalar>
Teuchos::RCP<StepperSDIRK_21Pair<Scalar> >
createStepperSDIRK_21Pair(
  const Teuchos::RCP<const Thyra::ModelEvaluator<Scalar> >& model,
  Teuchos::RCP<Teuchos::ParameterList> pl)
{
  auto stepper = Teuchos::rcp(new StepperSDIRK_21Pair<Scalar>());
  stepper->setStepperDIRKValues(pl);

  if (model != Teuchos::null) {
    stepper->setModel(model);
    stepper->initialize();
  }

  return stepper;
}


// ----------------------------------------------------------------------------
template<class Scalar>
class StepperDIRK_General :
  virtual public StepperDIRK<Scalar>
{
public:
  StepperDIRK_General()
  {
    this->setStepperName("General DIRK");
    this->setStepperType("General DIRK");
    this->setupTableau();
    this->setupDefault();
    this->setUseFSAL(            false);
    this->setICConsistency(      "None");
    this->setICConsistencyCheck( false);
  }

  StepperDIRK_General(
    const Teuchos::RCP<const Thyra::ModelEvaluator<Scalar> >& appModel,
    const Teuchos::RCP<Thyra::NonlinearSolverBase<Scalar> >& solver,
    bool useFSAL,
    std::string ICConsistency,
    bool ICConsistencyCheck,
    bool useEmbedded,
    bool zeroInitialGuess,
    const Teuchos::RCP<StepperRKAppAction<Scalar> >& stepperRKAppAction,
    const Teuchos::SerialDenseMatrix<int,Scalar>& A,
    const Teuchos::SerialDenseVector<int,Scalar>& b,
    const Teuchos::SerialDenseVector<int,Scalar>& c,
    const int order,
    const int orderMin,
    const int orderMax,
    const Teuchos::SerialDenseVector<int,Scalar>& bstar)
  {
    this->setStepperName("General DIRK");
    this->setStepperType("General DIRK");
    this->setTableau(A,b,c,order,orderMin,orderMax,bstar);

    TEUCHOS_TEST_FOR_EXCEPTION(
      this->tableau_->isImplicit() != true, std::logic_error,
      "Error - General DIRK did not receive a DIRK Butcher Tableau!\n");

    this->setup(appModel, solver, useFSAL, ICConsistency, ICConsistencyCheck,
                useEmbedded, zeroInitialGuess, stepperRKAppAction);
  }

  std::string getDescription() const
  {
    std::stringstream Description;
    Description << this->getStepperType() << "\n"
      << "The format of the Butcher Tableau parameter list is\n"
      << "  <Parameter name=\"A\" type=\"string\" value=\"# # # ;\n"
      << "                                           # # # ;\n"
      << "                                           # # #\"/>\n"
      << "  <Parameter name=\"b\" type=\"string\" value=\"# # #\"/>\n"
      << "  <Parameter name=\"c\" type=\"string\" value=\"# # #\"/>\n\n"
      << "Note the number of stages is implicit in the number of entries.\n"
      << "The number of stages must be consistent.\n"
      << "\n"
      << "Default tableau is 'SDIRK 2 Stage 2nd order':\n"
      << "  Computer Methods for ODEs and DAEs\n"
      << "  U. M. Ascher and L. R. Petzold\n"
      << "  p. 106\n"
      << "  gamma = (2-sqrt(2))/2\n"
      << "  c = [  gamma   1     ]'\n"
      << "  A = [  gamma   0     ]\n"
      << "      [ 1-gamma  gamma ]\n"
      << "  b = [ 1-gamma  gamma ]'";
    return Description.str();
  }

  void setupTableau()
  {
    if (this->tableau_ == Teuchos::null) {
      // Set tableau to the default if null, otherwise keep current tableau.
      auto stepper = Teuchos::rcp(new StepperSDIRK_2Stage2ndOrder<Scalar>());
      auto t = stepper->getTableau();
      this->tableau_ = Teuchos::rcp(new RKButcherTableau<Scalar>(
                                 this->getStepperType(),
                                 t->A(),t->b(),t->c(),
                                 t->order(),t->orderMin(),t->orderMax(),
                                 t->bstar()));
      this->isInitialized_ = false;
    }
  }

  void setTableau(const Teuchos::SerialDenseMatrix<int,Scalar>& A,
                  const Teuchos::SerialDenseVector<int,Scalar>& b,
                  const Teuchos::SerialDenseVector<int,Scalar>& c,
                  const int order,
                  const int orderMin,
                  const int orderMax,
                  const Teuchos::SerialDenseVector<int,Scalar>&
                    bstar = Teuchos::SerialDenseVector<int,Scalar>())
  {
    this->tableau_ = Teuchos::rcp(new RKButcherTableau<Scalar>(
      this->getStepperType(),A,b,c,order,orderMin,orderMax,bstar));
    this->isInitialized_ = false;
  }

  Teuchos::RCP<const Teuchos::ParameterList>
  getValidParameters() const
  {
    auto pl = this->getValidParametersBasicDIRK();

    // Tableau ParameterList
    Teuchos::SerialDenseMatrix<int,Scalar> A = this->tableau_->A();
    Teuchos::SerialDenseVector<int,Scalar> b = this->tableau_->b();
    Teuchos::SerialDenseVector<int,Scalar> c = this->tableau_->c();
    Teuchos::SerialDenseVector<int,Scalar> bstar = this->tableau_->bstar();

    Teuchos::RCP<Teuchos::ParameterList> tableauPL = Teuchos::parameterList();

    std::ostringstream Astream;
    Astream.precision(15);
    for(int i = 0; i < A.numRows(); i++) {
      for(int j = 0; j < A.numCols()-1; j++) {
        Astream << A(i,j) << " ";
      }
      Astream << A(i,A.numCols()-1);
      if ( i != A.numRows()-1 ) Astream << "; ";
    }
    tableauPL->set<std::string>("A", Astream.str());

    std::ostringstream bstream;
    bstream.precision(15);
    for(int i = 0; i < b.length()-1; i++) {
        bstream << b(i) << " ";
    }
    bstream << b(b.length()-1);
    tableauPL->set<std::string>("b", bstream.str());

    std::ostringstream cstream;
    cstream.precision(15);
    for(int i = 0; i < c.length()-1; i++) {
        cstream << c(i) << " ";
    }
    cstream << c(c.length()-1);
    tableauPL->set<std::string>("c", cstream.str());

    tableauPL->set<int>("order", this->tableau_->order());

    if ( bstar.length() == 0 ) {
      tableauPL->set("bstar", "");
    } else {
      std::ostringstream bstarstream;
      bstarstream.precision(15);
      for(int i = 0; i < bstar.length()-1; i++) {
          bstarstream << bstar(i) << " ";
      }
      bstarstream << bstar(bstar.length()-1);
      tableauPL->set<std::string>("bstar", bstarstream.str());
    }

    pl->set("Tableau", *tableauPL);

    return pl;
  }
};


// ------------------------------------------------------------------------
template<class Scalar>
Teuchos::RCP<StepperDIRK_General<Scalar> >
createStepperDIRK_General(
  const Teuchos::RCP<const Thyra::ModelEvaluator<Scalar> >& model,
  Teuchos::RCP<Teuchos::ParameterList> pl)
{
  auto stepper = Teuchos::rcp(new StepperDIRK_General<Scalar>());
  stepper->setStepperDIRKValues(pl);

  if (pl != Teuchos::null) {
    if (pl->isParameter("Tableau")) {
      auto t = stepper->createTableau(pl);
      stepper->setTableau( t->A(),t->b(),t->c(),
                           t->order(),t->orderMin(),t->orderMax(),
                           t->bstar() );
    }
  }
  TEUCHOS_TEST_FOR_EXCEPTION(
    stepper->getTableau()->isDIRK() != true, std::logic_error,
    "Error - General DIRK did not receive a DIRK Butcher Tableau!\n");

  if (model != Teuchos::null) {
    stepper->setModel(model);
    stepper->initialize();
  }

  return stepper;
}


} // namespace Tempus


#endif // Tempus_StepperRKButcherTableau_hpp