gradientFlow_AdaptiveRungeKutta.cpp

Go to the documentation of this file.

#include "gradientFlow_AdaptiveRungeKutta.h"

const std::string GradientFlow_AdaptiveRungeKutta::class_name = "GradientFlow_AdaptiveRungeKutta";

//====================================================================
void GradientFlow_AdaptiveRungeKutta::flow(double& t, double& Estep, Field_G& U)
{
  while (true)
  {
    copy(m_U_rough, U);
    copy(m_U_fine, U);

    //- rough step
    double t2     = t;
    double estep2 = 2.0 * Estep;

    m_RK->flow(t2, estep2, m_U_rough);

    //- two original steps
    double t1     = t;
    double estep1 = Estep;

    m_RK->flow(t1, estep1, m_U_fine);
    m_RK->flow(t1, estep1, m_U_fine);

    double diff = max_diff_U(m_U_rough, m_U_fine);

    vout.general(m_vl, "  Estep,diff,m_tolerance = %e %e %e\n", Estep, diff, m_tolerance);

    if (diff < m_tolerance) {
      copy(U, m_U_fine);
      t += estep2;

      // extend stepsize
      Estep *= m_safety * pow(m_tolerance / diff, 1.0 / (Norder_RK() + 1));

      break;  // exit loop
    } else {
      // shrink stepsize
      Estep *= m_safety * pow(m_tolerance / diff, 1.0 / (Norder_RK() + 1));

      if (Estep < m_tolerance) {
        vout.crucial(m_vl, "Error at %s: too small Estep = %e\n", class_name.c_str(), Estep);
        exit(EXIT_FAILURE);
      }

      // proceed to next trial
    }
  }
}


//====================================================================
double GradientFlow_AdaptiveRungeKutta::max_diff_U(const Field_G& U1, const Field_G& U0) const
{
  const int Nvol = U1.nvol();
  const int Nex  = U0.nex();
  const int Nc   = CommonParameters::Nc();

  double max_norm = 0.0;

  for (int ex = 0; ex < Nex; ++ex) {
    for (int site = 0; site < Nvol; ++site) {
      Mat_SU_N u0(Nc);
      U1.mat(u0, site, ex);

      Mat_SU_N u1(Nc);
      U0.mat(u1, site, ex);

      Mat_SU_N u2(Nc);
      u2 = u1 - u0;

      double norm = sqrt(u2.norm2()) / Nc;

      if (norm > max_norm) max_norm = norm;
    }
  }

  max_norm = Communicator::reduce_max(max_norm);

  return max_norm;
}


//====================================================================
//============================================================END=====