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);
  }
};

// ------------------------------------------------------------------------
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