Tempus_UnitTest_ERK.cpp

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

#include "Tempus_UnitTest_RK_Utils.hpp"

#include "Tempus_StepperRKButcherTableau.hpp"
#include "Tempus_StepperRKBase.hpp"
#include "Tempus_StepperRKModifierBase.hpp"
#include "Tempus_StepperRKModifierXBase.hpp"

#include "../TestModels/DahlquistTestModel.hpp"

namespace Tempus_Unit_Test {

using Teuchos::RCP;
using Teuchos::rcp;
using Teuchos::rcp_const_cast;
using Teuchos::rcp_dynamic_cast;

class StepperRKModifierBogackiShampineTest
  : virtual public Tempus::StepperRKModifierBase<double> {
 public:
  StepperRKModifierBogackiShampineTest(Teuchos::FancyOStream &Out,
                                       bool &Success)
    : out(Out), success(Success)
  {
  }

  virtual ~StepperRKModifierBogackiShampineTest() {}

  virtual void modify(
      Teuchos::RCP<Tempus::SolutionHistory<double> > sh,
      Teuchos::RCP<Tempus::StepperRKBase<double> > stepper,
      const typename Tempus::StepperRKAppAction<double>::ACTION_LOCATION actLoc)
  {
    const double relTol                       = 1.0e-14;
    auto stageNumber                          = stepper->getStageNumber();
    Teuchos::SerialDenseVector<int, double> c = stepper->getTableau()->c();

    auto currentState = sh->getCurrentState();
    auto workingState = sh->getWorkingState();
    const double dt   = workingState->getTimeStep();
    double time       = currentState->getTime();
    if (stageNumber >= 0) time += c(stageNumber) * dt;

    auto x    = workingState->getX();
    auto xDot = workingState->getXDot();
    if (xDot == Teuchos::null) xDot = stepper->getStepperXDot();

    if (actLoc == StepperRKAppAction<double>::BEGIN_STEP) {
      {
        auto DME = Teuchos::rcp_dynamic_cast<
            const Tempus_Test::DahlquistTestModel<double> >(
            stepper->getModel());
        TEST_FLOATING_EQUALITY(DME->getLambda(), -1.0, relTol);
      }
      TEST_FLOATING_EQUALITY(dt, 1.0, relTol);

      const double x_0    = get_ele(*(x), 0);
      const double xDot_0 = get_ele(*(xDot), 0);
      TEST_FLOATING_EQUALITY(x_0, 1.0, relTol);      // Should be x_0
      TEST_FLOATING_EQUALITY(xDot_0, -1.0, relTol);  // Should be xDot_0
      TEST_ASSERT(std::abs(time) < relTol);
      TEST_FLOATING_EQUALITY(dt, 1.0, relTol);
      TEST_COMPARE(stageNumber, ==, -1);
    }
    else if (actLoc == StepperRKAppAction<double>::BEGIN_STAGE ||
             actLoc == StepperRKAppAction<double>::BEFORE_SOLVE ||
             actLoc == StepperRKAppAction<double>::AFTER_SOLVE ||
             actLoc == StepperRKAppAction<double>::BEFORE_EXPLICIT_EVAL ||
             actLoc == StepperRKAppAction<double>::END_STAGE) {
      const double X_i = get_ele(*(x), 0);
      const double f_i = get_ele(*(xDot), 0);
      if (stageNumber == 0) {
        TEST_FLOATING_EQUALITY(X_i, 1.0, relTol);  // Should be X_1
        TEST_ASSERT(std::abs(time) < relTol);
        if (actLoc == StepperRKAppAction<double>::END_STAGE ||
            !stepper->getUseFSAL()) {
          TEST_FLOATING_EQUALITY(f_i, -1.0, relTol);  // Should be \bar{f}_1
        }
        else {
          TEST_ASSERT(std::abs(f_i) < relTol);  // Not set yet.
        }
      }
      else if (stageNumber == 1) {
        TEST_FLOATING_EQUALITY(X_i, 0.5, relTol);  // Should be X_2
        TEST_FLOATING_EQUALITY(time, 0.5, relTol);
        if (actLoc == StepperRKAppAction<double>::END_STAGE) {
          TEST_FLOATING_EQUALITY(f_i, -0.5, relTol);  // Should be \bar{f}_2
        }
        else {
          TEST_ASSERT(std::abs(f_i) < relTol);  // Not set yet.
        }
      }
      else if (stageNumber == 2) {
        TEST_FLOATING_EQUALITY(X_i, 5.0 / 8.0, relTol);  // Should be X_3
        TEST_FLOATING_EQUALITY(time, 0.75, relTol);
        if (actLoc == StepperRKAppAction<double>::END_STAGE) {
          TEST_FLOATING_EQUALITY(f_i, -5.0 / 8.0,
                                 relTol);  // Should be \bar{f}_3
        }
        else {
          TEST_ASSERT(std::abs(f_i) < relTol);  // Not set yet.
        }
      }
      else if (stageNumber == 3) {
        TEST_FLOATING_EQUALITY(X_i, 1.0 / 3.0, relTol);  // Should be X_4
        TEST_FLOATING_EQUALITY(time, 1.0, relTol);
        if (actLoc == StepperRKAppAction<double>::END_STAGE) {
          TEST_FLOATING_EQUALITY(f_i, -1.0 / 3.0,
                                 relTol);  // Should be \bar{f}_4
        }
        else if (workingState->getNConsecutiveFailures() > 0) {
          TEST_FLOATING_EQUALITY(f_i, -1.0, relTol);  // From IC and FSAL swap.
        }
        else {
          TEST_ASSERT(std::abs(f_i) < relTol);  // Not set yet.
        }
      }
      else {
        TEUCHOS_TEST_FOR_EXCEPT(!(-1 < stageNumber && stageNumber < 4));
      }
      TEST_FLOATING_EQUALITY(dt, 1.0, relTol);
    }
    else if (actLoc == StepperRKAppAction<double>::END_STEP) {
      const double x_1 = get_ele(*(x), 0);
      time             = workingState->getTime();
      TEST_FLOATING_EQUALITY(x_1, 1.0 / 3.0, relTol);  // Should be x_1
      TEST_FLOATING_EQUALITY(time, 1.0, relTol);
      TEST_FLOATING_EQUALITY(dt, 1.0, relTol);
      TEST_COMPARE(stageNumber, ==, -1);

      if (stepper->getUseEmbedded() == true) {
        TEST_FLOATING_EQUALITY(workingState->getTolRel(), 1.0, relTol);
        TEST_ASSERT(std::abs(workingState->getTolAbs()) < relTol);
        // e = 0 from doxygen above.
        TEST_ASSERT(std::abs(workingState->getErrorRel()) < relTol);
      }
    }
    else {
      TEUCHOS_TEST_FOR_EXCEPTION(true, std::logic_error,
                                 "Error - unknown action location.\n");
    }
  }

 private:
  Teuchos::FancyOStream &out;
  bool &success;
};

// ************************************************************
// ************************************************************
TEUCHOS_UNIT_TEST(ERK, Modifier)
{
  auto stepper = rcp(new Tempus::StepperERK_BogackiShampine32<double>());
  Teuchos::RCP<const Thyra::ModelEvaluator<double> > model =
      rcp(new Tempus_Test::DahlquistTestModel<double>());
  auto modifier = rcp(new StepperRKModifierBogackiShampineTest(out, success));

  stepper->setModel(model);
  stepper->setAppAction(modifier);
  stepper->setICConsistency("Consistent");
  stepper->setUseFSAL(false);
  stepper->initialize();

  // Create a SolutionHistory.
  auto solutionHistory = Tempus::createSolutionHistoryME(model);

  // Take one time step.
  stepper->setInitialConditions(solutionHistory);
  solutionHistory->initWorkingState();
  double dt = 1.0;
  solutionHistory->getWorkingState()->setTimeStep(dt);
  solutionHistory->getWorkingState()->setTime(dt);
  stepper->takeStep(solutionHistory);  // Primary testing occurs in
                                       // modifierX during takeStep().
  // Test stepper properties.
  TEUCHOS_ASSERT(stepper->getOrder() == 3);
}

// ************************************************************
// ************************************************************
TEUCHOS_UNIT_TEST(ERK, Modifier_FSAL)
{
  auto stepper = rcp(new Tempus::StepperERK_BogackiShampine32<double>());
  Teuchos::RCP<const Thyra::ModelEvaluator<double> > model =
      rcp(new Tempus_Test::DahlquistTestModel<double>());
  auto modifier = rcp(new StepperRKModifierBogackiShampineTest(out, success));

  stepper->setModel(model);
  stepper->setAppAction(modifier);
  stepper->setICConsistency("Consistent");
  stepper->setUseFSAL(true);
  stepper->initialize();

  // Create a SolutionHistory.
  auto solutionHistory = Tempus::createSolutionHistoryME(model);

  // Take one time step.
  stepper->setInitialConditions(solutionHistory);
  solutionHistory->initWorkingState();
  double dt = 1.0;
  solutionHistory->getWorkingState()->setTimeStep(dt);
  solutionHistory->getWorkingState()->setTime(dt);
  stepper->takeStep(solutionHistory);  // Primary testing occurs in
                                       // modifierX during takeStep().
  // Test stepper properties.
  TEUCHOS_ASSERT(stepper->getOrder() == 3);
}

// ************************************************************
// ************************************************************
TEUCHOS_UNIT_TEST(ERK, Modifier_FSAL_Failure)
{
  auto stepper = rcp(new Tempus::StepperERK_BogackiShampine32<double>());
  Teuchos::RCP<const Thyra::ModelEvaluator<double> > model =
      rcp(new Tempus_Test::DahlquistTestModel<double>());
  auto modifier = rcp(new StepperRKModifierBogackiShampineTest(out, success));

  stepper->setModel(model);
  stepper->setAppAction(modifier);
  stepper->setICConsistency("Consistent");
  stepper->setUseFSAL(true);
  stepper->initialize();

  // Create a SolutionHistory.
  auto solutionHistory = Tempus::createSolutionHistoryME(model);

  // Take one time step.
  stepper->setInitialConditions(solutionHistory);
  solutionHistory->initWorkingState();
  double dt = 1.0;
  solutionHistory->getWorkingState()->setTimeStep(dt);
  solutionHistory->getWorkingState()->setTime(dt);
  solutionHistory->getWorkingState()->setNConsecutiveFailures(1);
  stepper->takeStep(solutionHistory);  // Primary testing occurs in
                                       // modifierX during takeStep().
  // Test stepper properties.
  TEUCHOS_ASSERT(stepper->getOrder() == 3);
}

// ************************************************************
// ************************************************************
TEUCHOS_UNIT_TEST(ERK, Modifier_Embedded)
{
  auto stepper = rcp(new Tempus::StepperERK_BogackiShampine32<double>());
  Teuchos::RCP<const Thyra::ModelEvaluator<double> > model =
      rcp(new Tempus_Test::DahlquistTestModel<double>());
  auto modifier = rcp(new StepperRKModifierBogackiShampineTest(out, success));

  stepper->setModel(model);
  stepper->setAppAction(modifier);
  stepper->setICConsistency("Consistent");
  stepper->setUseFSAL(false);
  stepper->setUseEmbedded(true);
  stepper->initialize();

  // Create a SolutionHistory.
  auto solutionHistory = Tempus::createSolutionHistoryME(model);

  // Take one time step.
  stepper->setInitialConditions(solutionHistory);
  solutionHistory->initWorkingState();
  double dt = 1.0;
  solutionHistory->getWorkingState()->setTimeStep(dt);
  solutionHistory->getWorkingState()->setTime(dt);
  solutionHistory->getWorkingState()->setTolRel(1.0);
  solutionHistory->getWorkingState()->setTolAbs(0.0);
  stepper->takeStep(solutionHistory);  // Primary testing occurs in
                                       // modifierX during takeStep().
  // Test stepper properties.
  TEUCHOS_ASSERT(stepper->getOrder() == 3);
}

class StepperRKModifierXBogackiShampineTest
  : virtual public Tempus::StepperRKModifierXBase<double> {
 public:
  StepperRKModifierXBogackiShampineTest(Teuchos::FancyOStream &Out,
                                        bool &Success)
    : out(Out), success(Success)
  {
  }

  virtual ~StepperRKModifierXBogackiShampineTest() {}

  virtual void modify(
      Teuchos::RCP<Thyra::VectorBase<double> > x, const double time,
      const double dt, const int stageNumber,
      const typename Tempus::StepperRKModifierXBase<double>::MODIFIER_TYPE
          modType)
  {
    const double relTol = 1.0e-14;
    if (modType == StepperRKModifierXBase<double>::X_BEGIN_STEP) {
      const double x_0 = get_ele(*(x), 0);  // Should be x_0
      TEST_FLOATING_EQUALITY(x_0, 1.0, relTol);
      TEST_FLOATING_EQUALITY(time, 0.0, relTol);
      TEST_FLOATING_EQUALITY(dt, 1.0, relTol);
      TEST_COMPARE(stageNumber, ==, -1);
    }
    else if (modType == StepperRKModifierXBase<double>::X_BEGIN_STAGE ||
             modType == StepperRKModifierXBase<double>::X_BEFORE_SOLVE ||
             modType == StepperRKModifierXBase<double>::X_AFTER_SOLVE ||
             modType ==
                 StepperRKModifierXBase<double>::X_BEFORE_EXPLICIT_EVAL ||
             modType == StepperRKModifierXBase<double>::X_END_STAGE) {
      const double X_i = get_ele(*(x), 0);
      if (stageNumber == 0) {
        TEST_FLOATING_EQUALITY(X_i, 1.0, relTol);  // Should be X_1
        TEST_FLOATING_EQUALITY(time, 0.0, relTol);
      }
      else if (stageNumber == 1) {
        TEST_FLOATING_EQUALITY(X_i, 0.5, relTol);  // Should be X_2
        TEST_FLOATING_EQUALITY(time, 0.5, relTol);
      }
      else if (stageNumber == 2) {
        TEST_FLOATING_EQUALITY(X_i, 5.0 / 8.0, relTol);  // Should be X_3
        TEST_FLOATING_EQUALITY(time, 0.75, relTol);
      }
      else if (stageNumber == 3) {
        TEST_FLOATING_EQUALITY(X_i, 1.0 / 3.0, relTol);  // Should be X_4
        TEST_FLOATING_EQUALITY(time, 1.0, relTol);
      }
      else {
        TEUCHOS_TEST_FOR_EXCEPT(!(-1 < stageNumber && stageNumber < 4));
      }
      TEST_FLOATING_EQUALITY(dt, 1.0, relTol);
    }
    else if (modType == StepperRKModifierXBase<double>::X_END_STEP) {
      const double x_1 = get_ele(*(x), 0);
      TEST_FLOATING_EQUALITY(x_1, 1.0 / 3.0, relTol);  // Should be x_1
      TEST_FLOATING_EQUALITY(time, 1.0, relTol);
      TEST_FLOATING_EQUALITY(dt, 1.0, relTol);
      TEST_COMPARE(stageNumber, ==, -1);
    }
    else {
      TEUCHOS_TEST_FOR_EXCEPTION(true, std::logic_error,
                                 "Error - unknown action location.\n");
    }
  }

 private:
  Teuchos::FancyOStream &out;
  bool &success;
};

// ************************************************************
// ************************************************************
TEUCHOS_UNIT_TEST(ERK, ModifierX)
{
  auto stepper = rcp(new Tempus::StepperERK_BogackiShampine32<double>());
  Teuchos::RCP<const Thyra::ModelEvaluator<double> > model =
      rcp(new Tempus_Test::DahlquistTestModel<double>());
  auto modifierX = rcp(new StepperRKModifierXBogackiShampineTest(out, success));

  stepper->setModel(model);
  stepper->setAppAction(modifierX);
  stepper->setICConsistency("Consistent");
  stepper->setUseFSAL(false);
  stepper->initialize();

  // Create a SolutionHistory.
  auto solutionHistory = Tempus::createSolutionHistoryME(model);

  // Take one time step.
  stepper->setInitialConditions(solutionHistory);
  solutionHistory->initWorkingState();
  double dt = 1.0;
  solutionHistory->getWorkingState()->setTimeStep(dt);
  solutionHistory->getWorkingState()->setTime(dt);
  stepper->takeStep(solutionHistory);  // Primary testing occurs in
                                       // modifierX during takeStep().
  // Test stepper properties.
  TEUCHOS_ASSERT(stepper->getOrder() == 3);
}

// ************************************************************
// ************************************************************
TEUCHOS_UNIT_TEST(ERK, ModifierX_FSAL)
{
  auto stepper = rcp(new Tempus::StepperERK_BogackiShampine32<double>());
  Teuchos::RCP<const Thyra::ModelEvaluator<double> > model =
      rcp(new Tempus_Test::DahlquistTestModel<double>());
  auto modifierX = rcp(new StepperRKModifierXBogackiShampineTest(out, success));

  stepper->setModel(model);
  stepper->setAppAction(modifierX);
  stepper->setICConsistency("Consistent");
  stepper->setUseFSAL(true);
  stepper->initialize();

  // Create a SolutionHistory.
  auto solutionHistory = Tempus::createSolutionHistoryME(model);

  // Take one time step.
  stepper->setInitialConditions(solutionHistory);
  solutionHistory->initWorkingState();
  double dt = 1.0;
  solutionHistory->getWorkingState()->setTimeStep(dt);
  solutionHistory->getWorkingState()->setTime(dt);
  stepper->takeStep(solutionHistory);  // Primary testing occurs in
                                       // modifierX during takeStep().
  // Test stepper properties.
  TEUCHOS_ASSERT(stepper->getOrder() == 3);
}

// ************************************************************
// ************************************************************
TEUCHOS_UNIT_TEST(ERK, ModifierX_FSAL_Failure)
{
  auto stepper = rcp(new Tempus::StepperERK_BogackiShampine32<double>());
  Teuchos::RCP<const Thyra::ModelEvaluator<double> > model =
      rcp(new Tempus_Test::DahlquistTestModel<double>());
  auto modifierX = rcp(new StepperRKModifierXBogackiShampineTest(out, success));

  stepper->setModel(model);
  stepper->setAppAction(modifierX);
  stepper->setICConsistency("Consistent");
  stepper->setUseFSAL(true);
  stepper->initialize();

  // Create a SolutionHistory.
  auto solutionHistory = Tempus::createSolutionHistoryME(model);

  // Take one time step.
  stepper->setInitialConditions(solutionHistory);
  solutionHistory->initWorkingState();
  double dt = 1.0;
  solutionHistory->getWorkingState()->setTimeStep(dt);
  solutionHistory->getWorkingState()->setTime(dt);
  solutionHistory->getWorkingState()->setNConsecutiveFailures(1);
  stepper->takeStep(solutionHistory);  // Primary testing occurs in
                                       // modifierX during takeStep().
  // Test stepper properties.
  TEUCHOS_ASSERT(stepper->getOrder() == 3);
}

}  // namespace Tempus_Unit_Test