Tempus_ConvergenceTestUtils.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_CONVERGENCE_TEST_UTILS_HPP
#define TEMPUS_CONVERGENCE_TEST_UTILS_HPP

#include <vector>
#include "Teuchos_as.hpp"
#include "Tempus_Stepper.hpp"
#include <fstream>


namespace Tempus_Test {

/** \brief Linear regression class.
 *  Copied and modified from Rythmos.
 */
template<class Scalar>
class LinearRegression
{
  public:
    LinearRegression();
    void setData(std::vector<Scalar>& x, std::vector<Scalar>& y);
    Scalar getSlope() const;
    Scalar getYIntercept() const;
  private:
    // Private functions
    void compute_();
    void validateXYData_(std::vector<Scalar>& x, std::vector<Scalar>& y);

    // Private data
    std::vector<Scalar> x_;
    std::vector<Scalar> y_;
    Scalar slope_;
    Scalar yIntercept_;
    bool isInitialized_;
};

// Member functions:

template<class Scalar>
LinearRegression<Scalar>::LinearRegression()
{
  isInitialized_ = false;
}

template<class Scalar>
void LinearRegression<Scalar>::setData(
  std::vector<Scalar>& x, std::vector<Scalar>& y)
{
  validateXYData_(x,y);
  x_ = x; // copy x data
  y_ = y; // copy y data
  isInitialized_ = true;
  compute_();
}

template<class Scalar>
void LinearRegression<Scalar>::validateXYData_(
  std::vector<Scalar>& x, std::vector<Scalar>& y)
{
  TEUCHOS_TEST_FOR_EXCEPTION(x.size() != y.size(), std::logic_error,
    "x and y data are note the same size for linear regression\n");
  TEUCHOS_TEST_FOR_EXCEPTION(x.size() < 2, std::logic_error,
    "Not enough data points for linear regression!\n");
  int N = Teuchos::as<int>(x.size());
  // There must be at least two unique x values
  Scalar alpha = x[0];
  int numUnique = 1;
  for (int i=1; i<N ; ++i) {
    if (x[i] != alpha) {
      numUnique++;
    }
  }
  TEUCHOS_TEST_FOR_EXCEPT(numUnique==1);
}

template<class Scalar>
Scalar LinearRegression<Scalar>::getSlope() const
{
  TEUCHOS_TEST_FOR_EXCEPT(!isInitialized_);
  return slope_;
}

template<class Scalar>
Scalar LinearRegression<Scalar>::getYIntercept() const
{
  TEUCHOS_TEST_FOR_EXCEPT(!isInitialized_);
  return yIntercept_;
}

template<class Scalar>
void LinearRegression<Scalar>::compute_()
{
  TEUCHOS_TEST_FOR_EXCEPT(!isInitialized_);
  typedef Teuchos::ScalarTraits<Scalar> ST;

  int N = Teuchos::as<int>(x_.size());

  Scalar sum1 = ST::zero();
  Scalar sum2 = ST::zero();
  for (int i=0 ; i<N ; ++i) {
    sum1 += x_[i]*y_[i];
    sum2 += x_[i]*x_[i];
  }
  sum1 *= Scalar(-2*ST::one());
  sum2 *= Scalar(-2*ST::one());

  Scalar sum3 = ST::zero();
  Scalar sum4 = ST::zero();
  for (int i=0 ; i<N ; ++i) {
    for (int j=0 ; j<N ; ++j) {
      sum3 += x_[i]*y_[j];
      sum4 += x_[i]*x_[j];
    }
  }
  sum3 *= Scalar(2*ST::one()/Scalar(N));
  sum4 *= Scalar(2*ST::one()/Scalar(N));

  slope_ = ( sum3 + sum1 ) / ( sum4 + sum2 );

  yIntercept_ = ST::zero();
  for (int i=0 ; i<N ; ++i ) {
    yIntercept_ += y_[i]-slope_*x_[i];
  }
  yIntercept_ *= Scalar(ST::one()/Scalar(N));
}

// Nonmember helper functions:
template<class Scalar>
Scalar computeLinearRegression(
  std::vector<Scalar>& x, std::vector<Scalar>& y)
{
  LinearRegression<Scalar> lr;
  lr.setData(x,y);
  return(lr.getSlope());
}

template<class Scalar>
void computeLinearRegression(
  std::vector<Scalar>& x, std::vector<Scalar>& y,
  Scalar & slope, Scalar & yIntercept)
{
  LinearRegression<Scalar> lr;
  lr.setData(x,y);
  slope = lr.getSlope();
  yIntercept = lr.getYIntercept();
  return;
}

template<class Scalar>
Scalar computeLinearRegressionLogLog(
  std::vector<Scalar>& x, std::vector<Scalar>& y)
{
  TEUCHOS_TEST_FOR_EXCEPT(x.size() != y.size());
  int N = Teuchos::as<int>(x.size());
  std::vector<Scalar> xlog;
  std::vector<Scalar> ylog;

  for (int i=0 ; i<N ; ++i) {
    xlog.push_back(log(x[i]));
    ylog.push_back(log(y[i]));
  }

  LinearRegression<Scalar> lr;
  lr.setData(xlog,ylog);
  return(lr.getSlope());
}

template<class Scalar>
Teuchos::RCP<LinearRegression<Scalar> > linearRegression()
{
  Teuchos::RCP<LinearRegression<Scalar> > lr =
    Teuchos::rcp(new LinearRegression<Scalar>());
  return lr;
}

template<class Scalar>
void writeOrderError(
  const std::string filename,
  Teuchos::RCP<Tempus::Stepper<Scalar> > stepper,
  std::vector<Scalar> & StepSize,
  std::vector<Teuchos::RCP<Thyra::VectorBase<Scalar>>> & solutions,
  std::vector<Scalar> & xErrorNorm,
  Scalar & xSlope,
  std::vector<Teuchos::RCP<Thyra::VectorBase<Scalar>>> & solutionsDot,
  std::vector<Scalar> & xDotErrorNorm,
  Scalar & xDotSlope,
  std::vector<Teuchos::RCP<Thyra::VectorBase<Scalar>>> & solutionsDotDot,
  std::vector<Scalar> & xDotDotErrorNorm,
  Scalar & xDotDotSlope)
{
  if ( solutions.empty() ) {
    std::cout << "No solutions to write out!\n" << std::endl;
    return;
  }
  std::size_t lastElem = solutions.size()-1;  // Last element is the ref soln.
  std::vector<Scalar> dt;

  auto ref_solution = solutions[lastElem];
  for (std::size_t i=0; i < lastElem; ++i) {
    dt.push_back(StepSize[i]);
    auto tmp = solutions[i];
    Thyra::Vp_StV(tmp.ptr(), -1.0, *ref_solution);
    Scalar L2norm = Thyra::norm_2(*tmp);
    xErrorNorm.push_back(L2norm);
  }
  xSlope = computeLinearRegressionLogLog<Scalar>(dt, xErrorNorm);

  if (!solutionsDot.empty()) {
    auto ref_solutionDot = solutionsDot[lastElem];
    for (std::size_t i=0; i < lastElem; ++i) {
      auto tmp = solutionsDot[i];
      Thyra::Vp_StV(tmp.ptr(), -1.0, *ref_solutionDot);
      Scalar L2norm = Thyra::norm_2(*tmp);
      xDotErrorNorm.push_back(L2norm);
    }
    xDotSlope = computeLinearRegressionLogLog<Scalar>(dt, xDotErrorNorm);
  }

  if (!solutionsDotDot.empty()) {
    auto ref_solutionDotDot = solutionsDotDot[solutions.size()-1];
    for (std::size_t i=0; i < lastElem; ++i) {
      auto tmp = solutionsDotDot[i];
      Thyra::Vp_StV(tmp.ptr(), -1.0, *ref_solutionDotDot);
      Scalar L2norm = Thyra::norm_2(*tmp);
      xDotDotErrorNorm.push_back(L2norm);
    }
    xDotDotSlope = computeLinearRegressionLogLog<Scalar>(dt, xDotDotErrorNorm);
  }

  Scalar order = stepper->getOrder();
  std::cout << "  Stepper = " << stepper->description() << std::endl;
  std::cout << "  =================================" << std::endl;
  std::cout << "  Expected Order = " << order        << std::endl;
  std::cout << "  x        Order = " << xSlope       << std::endl;
  if (!solutionsDot.empty())
    std::cout<<"  xDot     Order = " << xDotSlope    << std::endl;
  if (!solutionsDotDot.empty())
    std::cout<<"  xDotDot  Order = " << xDotDotSlope << std::endl;
  std::cout << "  =================================" << std::endl;

  std::ofstream ftmp(filename);
  Scalar factor = 0.0;
  for (std::size_t n=0; n<lastElem; n++) {
    factor = 0.8*(pow(dt[n]/dt[0], order));
    ftmp << dt[n]
         << "   " << xErrorNorm[n]       << "   " << factor*xErrorNorm[0];
    if (xDotErrorNorm.size() == lastElem)
      ftmp <<"   "<< xDotErrorNorm[n]    << "   " << factor*xDotErrorNorm[0];
    if (xDotDotErrorNorm.size() == lastElem)
      ftmp <<"   "<< xDotDotErrorNorm[n] << "   " << factor*xDotDotErrorNorm[0];
    ftmp << std::endl;
  }
  ftmp.close();
}

template<class Scalar>
void writeOrderError(
  const std::string filename,
  Teuchos::RCP<Tempus::Stepper<Scalar> > stepper,
  std::vector<Scalar> & StepSize,
  std::vector<Teuchos::RCP<Thyra::VectorBase<Scalar>>> & solutions,
  std::vector<Scalar> & xErrorNorm,
  Scalar & xSlope,
  std::vector<Teuchos::RCP<Thyra::VectorBase<Scalar>>> & solutionsDot,
  std::vector<Scalar> & xDotErrorNorm,
  Scalar & xDotSlope)
{
  std::vector<Teuchos::RCP<Thyra::VectorBase<Scalar>>> solutionsDotDot;
  std::vector<Scalar> xDotDotErrorNorm;
  Scalar xDotDotSlope = Scalar(0.0);
  writeOrderError(filename, stepper, StepSize,
                  solutions,       xErrorNorm,       xSlope,
                  solutionsDot,    xDotErrorNorm,    xDotSlope,
                  solutionsDotDot, xDotDotErrorNorm, xDotDotSlope);
}

template<class Scalar>
void writeOrderError(
  const std::string filename,
  Teuchos::RCP<Tempus::Stepper<Scalar> > stepper,
  std::vector<Scalar> & StepSize,
  std::vector<Teuchos::RCP<Thyra::VectorBase<Scalar>>> & solutions,
  std::vector<Scalar> & xErrorNorm,
  Scalar & xSlope)
{
  std::vector<Teuchos::RCP<Thyra::VectorBase<Scalar>>> solutionsDot;
  std::vector<Scalar> xDotErrorNorm;
  Scalar xDotSlope    = Scalar(0.0);
  std::vector<Teuchos::RCP<Thyra::VectorBase<Scalar>>> solutionsDotDot;
  std::vector<Scalar> xDotDotErrorNorm;
  Scalar xDotDotSlope = Scalar(0.0);
  writeOrderError(filename, stepper, StepSize,
                  solutions,       xErrorNorm,       xSlope,
                  solutionsDot,    xDotErrorNorm,    xDotSlope,
                  solutionsDotDot, xDotDotErrorNorm, xDotDotSlope);
}

template<class Scalar>
void writeSolution(
  const std::string filename,
  Teuchos::RCP<const Tempus::SolutionHistory<Scalar> > solutionHistory)
{
  std::ofstream ftmp(filename);
  Teuchos::RCP<const Thyra::VectorBase<Scalar> > x;
  Teuchos::RCP<const Thyra::VectorBase<Scalar> > xDot;
  Teuchos::RCP<const Thyra::VectorBase<Scalar> > xDotDot;
  for (int i=0; i<solutionHistory->getNumStates(); i++) {
    Teuchos::RCP<const Tempus::SolutionState<Scalar> > solutionState =
      (*solutionHistory)[i];
    Scalar time = solutionState->getTime();
    x       = solutionState->getX();
    xDot    = solutionState->getXDot();
    xDotDot = solutionState->getXDotDot();
    int J = x->space()->dim();

    ftmp << time;
    for (int j=0; j<J; j++)   ftmp << "   " << get_ele(*x, j);
    if (xDot != Teuchos::null)
      for (int j=0; j<J; j++) ftmp << "   " << get_ele(*xDot, j);
    if (xDotDot != Teuchos::null)
      for (int j=0; j<J; j++) ftmp << "   " << get_ele(*xDotDot, j);
    ftmp << std::endl;
  }
  ftmp.close();
}

template<class Scalar>
void writeSolution(
  const std::string filename,
  Teuchos::RCP<Tempus::SolutionHistory<Scalar> > solutionHistory)
{
  Teuchos::RCP<const Tempus::SolutionHistory<Scalar> > sh(solutionHistory);
  writeSolution(filename, sh);
}

} // namespace Tempus_Test

#endif // TEMPUS_CONVERGENCE_TEST_UTILS_HPP