fopr_Chebyshev.cpp

Go to the documentation of this file.

#include "fopr_Chebyshev.h"

#ifdef USE_PARAMETERS_FACTORY
#include "parameters_factory.h"
#endif

//- parameter entries
namespace {
  void append_entry(Parameters& param)
  {
    param.Register_int("degree_of_polynomial", 0);
    param.Register_double("threshold_value", 0.0);
    param.Register_double("upper_bound", 0.0);

    param.Register_string("verbose_level", "NULL");
  }


#ifdef USE_PARAMETERS_FACTORY
  bool init_param = ParametersFactory::Register("Fopr.Chebyshev", append_entry);
#endif
}
//- end

//- parameters class
Parameters_Fopr_Chebyshev::Parameters_Fopr_Chebyshev() { append_entry(*this); }
//- end

//====================================================================
void Fopr_Chebyshev::set_parameters(const Parameters& params)
{
  const string str_vlevel = params.get_string("verbose_level");

  m_vl = vout.set_verbose_level(str_vlevel);

  //- fetch and check input parameters
  int    Np;
  double v_threshold, v_max;

  int err = 0;
  err += params.fetch_int("degree_of_polynomial", Np);
  err += params.fetch_double("threshold_value", v_threshold);
  err += params.fetch_double("upper_bound", v_max);

  if (err) {
    vout.crucial(m_vl, "Fopr_Chebyshev: fetch error, input parameter not found.\n");
    abort();
  }


  set_parameters(Np, v_threshold, v_max);
}


//====================================================================
void Fopr_Chebyshev::set_parameters(int Np, double v_threshold, double v_max)
{
  //- print input parameters
  vout.general(m_vl, "Parameters of Fopr_Chebyshev:\n");
  vout.general(m_vl, "  Np          = %d\n", Np);
  vout.general(m_vl, "  v_threshold = %16.8e\n", v_threshold);
  vout.general(m_vl, "  v_max       = %16.8e\n", v_max);


  //- range check
  int err = 0;
  err += ParameterCheck::non_negative(Np);
  // NB. v_threshold,v_max == 0 is allowed.

  if (err) {
    vout.crucial(m_vl, "Fopr_Chebyshev: parameter range check failed.\n");
    abort();
  }

  //- store values
  m_Npcb = Np;

  double b_max = v_max / v_threshold;
  double r     = 2.0 / (b_max * b_max - 1.0);
  double s     = v_threshold / sqrt(0.5 * r);

  m_Fcb1 = 2.0 / (s * s);
  m_Fcb2 = -(1.0 + r);

  vout.general(m_vl, "  Fcb1        = %16.8e\n", m_Fcb1);
  vout.general(m_vl, "  Fcb2        = %16.8e\n", m_Fcb2);
}


//====================================================================
const Field Fopr_Chebyshev::mult(const Field& w)
{
  std::valarray<Field> dj(3);
  int Nin  = w.nin();
  int Nvol = w.nvol();
  int Nex  = w.nex();

  Field v(Nin, Nvol, Nex);

  for (int k = 0; k < 3; ++k) {
    dj[k].reset(Nin, Nvol, Nex);
  }

  dj[0] = (-1.0) * w;
  dj[1] = 0.0;

  int jn  = 2;
  int jp1 = 1;
  int jp2 = 0;

  for (int j = m_Npcb; j >= 2; --j) {
    dj[jn]  = m_fopr->mult(dj[jp1]);
    dj[jn] *= m_Fcb1;
    dj[jn] += m_Fcb2 * dj[jp1];

    dj[jn] *= 2.0;
    dj[jn] -= 1.0 * dj[jp2];

    jn  = (jn + 1) % 3;
    jp1 = (jp1 + 1) % 3;
    jp2 = (jp2 + 1) % 3;
  }

  v  = m_fopr->mult(dj[jp1]);
  v *= m_Fcb1;
  v += m_Fcb2 * dj[jp1];
  v -= dj[jp2];


  return v;
}


//====================================================================
double Fopr_Chebyshev::mult(double x)
{
  std::valarray<double> dj(3);

  dj[0] = -1.0;
  dj[1] = 0.0;

  int jn  = 2;
  int jp1 = 1;
  int jp2 = 0;

  for (int j = m_Npcb; j >= 2; --j) {
    dj[jn]  = x * dj[jp1];
    dj[jn] *= m_Fcb1;
    dj[jn] += m_Fcb2 * dj[jp1];

    dj[jn] *= 2.0;
    dj[jn] -= 1.0 * dj[jp2];

    jn  = (jn + 1) % 3;
    jp1 = (jp1 + 1) % 3;
    jp2 = (jp2 + 1) % 3;
  }

  double v = x * dj[jp1];
  v *= m_Fcb1;
  v += m_Fcb2 * dj[jp1];
  v -= dj[jp2];

  return v;
}


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