Tempus_PhysicsStateTest.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 "Teuchos_UnitTestHarness.hpp"
#include "Teuchos_XMLParameterListHelpers.hpp"
#include "Teuchos_TimeMonitor.hpp"

#include "Thyra_VectorStdOps.hpp"

#include "Tempus_IntegratorBasic.hpp"
#include "Tempus_PhysicsStateTest_StepperForwardEuler.hpp"
#include "Tempus_PhysicsStateCounter.hpp"

#include "../TestModels/SinCosModel.hpp"
#include "../TestModels/VanDerPolModel.hpp"
#include "../TestUtils/Tempus_ConvergenceTestUtils.hpp"

#include <fstream>
#include <vector>

namespace Tempus_Test {

using Teuchos::RCP;
using Teuchos::ParameterList;
using Teuchos::sublist;
using Teuchos::getParametersFromXmlFile;

using Tempus::IntegratorBasic;
using Tempus::SolutionHistory;
using Tempus::SolutionState;


// ************************************************************
// ************************************************************
TEUCHOS_UNIT_TEST(PhysicsState, SinCos)
{
  std::vector<double> StepSize;
  std::vector<double> ErrorNorm;
  const int nTimeStepSizes = 7;
  double dt = 0.2;
  double order = 0.0;
  for (int n=0; n<nTimeStepSizes; n++) {

    // Read params from .xml file
    RCP<ParameterList> pList =
      getParametersFromXmlFile("Tempus_PhysicsState_SinCos.xml");

    //std::ofstream ftmp("PL.txt");
    //pList->print(ftmp);
    //ftmp.close();

    // Setup the SinCosModel
    RCP<ParameterList> scm_pl = sublist(pList, "SinCosModel", true);
    //RCP<SinCosModel<double> > model = sineCosineModel(scm_pl);
    RCP<SinCosModel<double> > model =
      Teuchos::rcp(new SinCosModel<double> (scm_pl));

    dt /= 2;

    // Setup the Integrator and reset initial time step
    RCP<ParameterList> pl = sublist(pList, "Tempus", true);
    pl->sublist("Demo Integrator")
       .sublist("Time Step Control").set("Initial Time Step", dt);
    RCP<Tempus::IntegratorBasic<double> > integrator =
      Tempus::integratorBasic<double>(pl, model);

    // Replace Tempus::StepperForwardEuler with
    // Tempus_Test::PhysicsStateTest_StepperForwardEuler
    RCP<ParameterList> stepperPL =
      integrator->getStepper()->getNonconstParameterList();
    Teuchos::RCP<Tempus::Stepper<double> > physicsStepper = Teuchos::rcp(
      new PhysicsStateTest_StepperForwardEuler<double>(model, stepperPL));
    integrator->setStepperWStepper(physicsStepper);
    order = integrator->getStepper()->getOrder();

    // Initial Conditions
    // During the Integrator construction, the initial SolutionState
    // is set by default to model->getNominalVales().get_x().  However,
    // the application can set it also by integrator->initializeSolutionHistory.
    RCP<Thyra::VectorBase<double> > x0 =
      model->getNominalValues().get_x()->clone_v();
    integrator->initializeSolutionHistory(0.0, x0);


    // Replace Tempus::PhysicsState with
    // Tempus_Test::PhysicsStateCounter
    RCP<PhysicsStateCounter<double> > pSC = Teuchos::rcp(
      new PhysicsStateCounter<double> ("PhysicsStateTest", 0));
    integrator->getSolutionHistory()->getCurrentState()->setPhysicsState(pSC);

    // Integrate to timeMax
    bool integratorStatus = integrator->advanceTime();
    TEST_ASSERT(integratorStatus)

    // Test PhysicsState
    Teuchos::RCP<Tempus::PhysicsState<double> > pS =
      integrator->getSolutionHistory()->getCurrentState()->getPhysicsState();
    TEST_EQUALITY(pS->getName(), "PhysicsStateTest");
    pSC = Teuchos::rcp_dynamic_cast<PhysicsStateCounter<double> >(pS);
    //std::cout << " Counter = " << pSC->getCounter() << std::endl;
    TEST_EQUALITY(pSC->getCounter(), (int)(1.0/dt));


    // Test if at 'Final Time'
    double time = integrator->getTime();
    double timeFinal = pl->sublist("Demo Integrator")
      .sublist("Time Step Control").get<double>("Final Time");
    TEST_FLOATING_EQUALITY(time, timeFinal, 1.0e-14);

    // Time-integrated solution and the exact solution
    RCP<Thyra::VectorBase<double> > x = integrator->getX();
    RCP<const Thyra::VectorBase<double> > x_exact =
      model->getExactSolution(time).get_x();

    // Plot sample solution and exact solution
    if (n == 0) {
      std::ofstream ftmp("Tempus_ForwardEuler_SinCos.dat");
      RCP<const SolutionHistory<double> > solutionHistory =
        integrator->getSolutionHistory();
      RCP<const Thyra::VectorBase<double> > x_exact_plot;
      for (int i=0; i<solutionHistory->getNumStates(); i++) {
        RCP<const SolutionState<double> > solutionState = (*solutionHistory)[i];
        double time_i = solutionState->getTime();
        RCP<const Thyra::VectorBase<double> > x_plot = solutionState->getX();
        x_exact_plot = model->getExactSolution(time_i).get_x();
        ftmp << time_i << "   "
             << Thyra::get_ele(*(x_plot), 0) << "   "
             << Thyra::get_ele(*(x_plot), 1) << "   "
             << Thyra::get_ele(*(x_exact_plot), 0) << "   "
             << Thyra::get_ele(*(x_exact_plot), 1) << std::endl;
      }
      ftmp.close();
    }

    // Calculate the error
    RCP<Thyra::VectorBase<double> > xdiff = x->clone_v();
    Thyra::V_StVpStV(xdiff.ptr(), 1.0, *x_exact, -1.0, *(x));
    StepSize.push_back(dt);
    const double L2norm = Thyra::norm_2(*xdiff);
    ErrorNorm.push_back(L2norm);
  }

  // Check the order and intercept
  double slope = computeLinearRegressionLogLog<double>(StepSize, ErrorNorm);
  std::cout << "  Stepper = ForwardEuler" << std::endl;
  std::cout << "  =========================" << std::endl;
  std::cout << "  Expected order: " << order << std::endl;
  std::cout << "  Observed order: " << slope << std::endl;
  std::cout << "  =========================" << std::endl;
  TEST_FLOATING_EQUALITY( slope, order, 0.01 );
  TEST_FLOATING_EQUALITY( ErrorNorm[0], 0.051123, 1.0e-4 );

  std::ofstream ftmp("Tempus_ForwardEuler_SinCos-Error.dat");
  double error0 = 0.8*ErrorNorm[0];
  for (int n=0; n<nTimeStepSizes; n++) {
    ftmp << StepSize[n]  << "   " << ErrorNorm[n] << "   "
         << error0*(pow(StepSize[n]/StepSize[0],order)) << std::endl;
  }
  ftmp.close();

  Teuchos::TimeMonitor::summarize();
}

} // namespace Tempus_Test