From fa5e4add47bfdb2d0c54486ebb99236b9db11326 Mon Sep 17 00:00:00 2001 From: Guido Cossu Date: Tue, 31 Oct 2017 18:20:38 +0000 Subject: [PATCH] Added support for anisotropy to the WilsonFermion class --- lib/qcd/action/fermion/WilsonFermion.cc | 42 ++++++++++++++++++++----- lib/qcd/action/fermion/WilsonFermion.h | 21 +++++++++++-- tests/qdpxx/Test_qdpxx_wilson.cc | 28 ++++++++++++++--- 3 files changed, 76 insertions(+), 15 deletions(-) diff --git a/lib/qcd/action/fermion/WilsonFermion.cc b/lib/qcd/action/fermion/WilsonFermion.cc index 19f9674d..55ef5a51 100644 --- a/lib/qcd/action/fermion/WilsonFermion.cc +++ b/lib/qcd/action/fermion/WilsonFermion.cc @@ -47,7 +47,8 @@ int WilsonFermionStatic::HandOptDslash; template WilsonFermion::WilsonFermion(GaugeField &_Umu, GridCartesian &Fgrid, GridRedBlackCartesian &Hgrid, RealD _mass, - const ImplParams &p) + const ImplParams &p, + const WilsonAnisotropyCoefficients &anis) : Kernels(p), _grid(&Fgrid), _cbgrid(&Hgrid), @@ -60,16 +61,41 @@ WilsonFermion::WilsonFermion(GaugeField &_Umu, GridCartesian &Fgrid, Umu(&Fgrid), UmuEven(&Hgrid), UmuOdd(&Hgrid), - _tmp(&Hgrid) + _tmp(&Hgrid), + anisotropyCoeff(anis) { // Allocate the required comms buffer ImportGauge(_Umu); + if (anisotropyCoeff.isAnisotropic){ + diag_mass = mass + 1.0 + (Nd-1)*(anisotropyCoeff.nu / anisotropyCoeff.xi_0); + } else { + diag_mass = 4.0 + mass; + } + + } template void WilsonFermion::ImportGauge(const GaugeField &_Umu) { GaugeField HUmu(_Umu._grid); - HUmu = _Umu * (-0.5); + + //Here multiply the anisotropy coefficients + if (anisotropyCoeff.isAnisotropic) + { + + for (int mu = 0; mu < Nd; mu++) + { + GaugeLinkField U_dir = (-0.5)*PeekIndex(_Umu, mu); + if (mu != anisotropyCoeff.t_direction) + U_dir *= (anisotropyCoeff.nu / anisotropyCoeff.xi_0); + + PokeIndex(HUmu, U_dir, mu); + } + } + else + { + HUmu = _Umu * (-0.5); + } Impl::DoubleStore(GaugeGrid(), Umu, HUmu); pickCheckerboard(Even, UmuEven, Umu); pickCheckerboard(Odd, UmuOdd, Umu); @@ -83,14 +109,14 @@ template RealD WilsonFermion::M(const FermionField &in, FermionField &out) { out.checkerboard = in.checkerboard; Dhop(in, out, DaggerNo); - return axpy_norm(out, 4 + mass, in, out); + return axpy_norm(out, diag_mass, in, out); } template RealD WilsonFermion::Mdag(const FermionField &in, FermionField &out) { out.checkerboard = in.checkerboard; Dhop(in, out, DaggerYes); - return axpy_norm(out, 4 + mass, in, out); + return axpy_norm(out, diag_mass, in, out); } template @@ -114,7 +140,7 @@ void WilsonFermion::MeooeDag(const FermionField &in, FermionField &out) { template void WilsonFermion::Mooee(const FermionField &in, FermionField &out) { out.checkerboard = in.checkerboard; - typename FermionField::scalar_type scal(4.0 + mass); + typename FermionField::scalar_type scal(diag_mass); out = scal * in; } @@ -127,7 +153,7 @@ void WilsonFermion::MooeeDag(const FermionField &in, FermionField &out) { template void WilsonFermion::MooeeInv(const FermionField &in, FermionField &out) { out.checkerboard = in.checkerboard; - out = (1.0/(4.0+mass))*in; + out = (1.0/(diag_mass))*in; } template @@ -204,7 +230,7 @@ void WilsonFermion::DerivInternal(StencilImpl &st, DoubledGaugeField &U, FermionField Btilde(B._grid); FermionField Atilde(B._grid); - Atilde = A; + Atilde = A;//redundant st.HaloExchange(B, compressor); diff --git a/lib/qcd/action/fermion/WilsonFermion.h b/lib/qcd/action/fermion/WilsonFermion.h index 50f4f884..7b2b5206 100644 --- a/lib/qcd/action/fermion/WilsonFermion.h +++ b/lib/qcd/action/fermion/WilsonFermion.h @@ -44,6 +44,19 @@ class WilsonFermionStatic { static const int npoint = 8; }; +struct WilsonAnisotropyCoefficients{ + bool isAnisotropic; + int t_direction; + double xi_0; + double nu; + + WilsonAnisotropyCoefficients(): + isAnisotropic(false), + t_direction(Nd-1), + xi_0(1.0), + nu(1.0){} +}; + template class WilsonFermion : public WilsonKernels, public WilsonFermionStatic { public: @@ -117,8 +130,9 @@ class WilsonFermion : public WilsonKernels, public WilsonFermionStatic { // Constructor WilsonFermion(GaugeField &_Umu, GridCartesian &Fgrid, - GridRedBlackCartesian &Hgrid, RealD _mass, - const ImplParams &p = ImplParams()); + GridRedBlackCartesian &Hgrid, RealD _mass, + const ImplParams &p = ImplParams(), + const WilsonAnisotropyCoefficients &anis = WilsonAnisotropyCoefficients() ); // DoubleStore impl dependent void ImportGauge(const GaugeField &_Umu); @@ -130,6 +144,7 @@ class WilsonFermion : public WilsonKernels, public WilsonFermionStatic { // protected: public: RealD mass; + RealD diag_mass; GridBase *_grid; GridBase *_cbgrid; @@ -146,6 +161,8 @@ class WilsonFermion : public WilsonKernels, public WilsonFermionStatic { LebesgueOrder Lebesgue; LebesgueOrder LebesgueEvenOdd; + + WilsonAnisotropyCoefficients anisotropyCoeff; }; typedef WilsonFermion WilsonFermionF; diff --git a/tests/qdpxx/Test_qdpxx_wilson.cc b/tests/qdpxx/Test_qdpxx_wilson.cc index 6e6eb213..a084ebad 100644 --- a/tests/qdpxx/Test_qdpxx_wilson.cc +++ b/tests/qdpxx/Test_qdpxx_wilson.cc @@ -32,7 +32,7 @@ #include // Mass -double mq = 0.02; +double mq = 0.1; // Define Wilson Types typedef Grid::QCD::WilsonImplR::FermionField FermionField; @@ -255,6 +255,12 @@ public: Chroma::WilsonFermActParams p; p.Mass = _mq; + AnisoParam_t _apar; + _apar.anisoP = true; + _apar.t_dir = 3; // in 4d + _apar.xi_0 = 2.0; + _apar.nu = 1.0; + p.anisoParam = _apar; Chroma::Handle> fbc(new Chroma::SimpleFermBC(bcs)); Chroma::Handle> cfs(new Chroma::CreateSimpleFermState(fbc)); @@ -269,7 +275,13 @@ public: p.Mass = _mq; p.clovCoeffR = QDP::Real(1.0); p.clovCoeffT = QDP::Real(1.0); - Real u0 = QDP::Real(1.0); + p.u0 = QDP::Real(1.0); + AnisoParam_t _apar; + _apar.anisoP = false; + _apar.t_dir = 3; // in 4d + _apar.xi_0 = 2.0; + _apar.nu = 1.0; + p.anisoParam = _apar; Chroma::Handle> fbc(new Chroma::SimpleFermBC(bcs)); Chroma::Handle> cfs(new Chroma::CreateSimpleFermState(fbc)); @@ -391,8 +403,13 @@ void calc_grid(ChromaAction action, Grid::QCD::LatticeGaugeField &Umu, Grid::QCD if (action == Wilson) { - - Grid::QCD::WilsonFermionR Wf(Umu, *UGrid, *UrbGrid, _mass); + WilsonAnisotropyCoefficients anis; + anis.isAnisotropic = true; + anis.t_direction = 3; + anis.xi_0 = 2.0; + anis.nu = 1.0; + WilsonImplParams iParam; + Grid::QCD::WilsonFermionR Wf(Umu, *UGrid, *UrbGrid, _mass, iParam, anis); std::cout << Grid::GridLogMessage << " Calling Grid Wilson Fermion multiply " << std::endl; @@ -406,7 +423,8 @@ void calc_grid(ChromaAction action, Grid::QCD::LatticeGaugeField &Umu, Grid::QCD if (action == WilsonClover) { Grid::RealD _csw = 1.0; - + WilsonAnisotropyCoefficients anis; + WilsonImplParams implParam; Grid::QCD::WilsonCloverFermionR Wf(Umu, *UGrid, *UrbGrid, _mass, _csw); Wf.ImportGauge(Umu);