diff --git a/lib/qcd/action/fermion/Fermion.h b/lib/qcd/action/fermion/Fermion.h index 2eaf42fa..9d999c6d 100644 --- a/lib/qcd/action/fermion/Fermion.h +++ b/lib/qcd/action/fermion/Fermion.h @@ -49,7 +49,7 @@ Author: Peter Boyle #include // 4d wilson like #include // 4d wilson like -#include // 4d wilson like +#include // 4d wilson clover fermions #include // 5d base used by all 5d overlap types #include diff --git a/lib/qcd/action/fermion/WilsonCloverFermion.cc b/lib/qcd/action/fermion/WilsonCloverFermion.cc index c7fd211d..ebea565b 100644 --- a/lib/qcd/action/fermion/WilsonCloverFermion.cc +++ b/lib/qcd/action/fermion/WilsonCloverFermion.cc @@ -27,28 +27,35 @@ *************************************************************************************/ /* END LEGAL */ #include +#include #include namespace Grid { namespace QCD { - template - void WilsonCloverFermion::AddCloverTerm(const FermionField& in, - FermionField& out){ - FermionField tmp(out._grid); - tmp = zero; - // the product sigma_munu Fmunu is hermitian - tmp += Bx*(Gamma(Gamma::Algebra::SigmaYZ)*in); - tmp += By*(Gamma(Gamma::Algebra::MinusSigmaXZ)*in); - tmp += Bz*(Gamma(Gamma::Algebra::SigmaXY)*in); - tmp += Ex*(Gamma(Gamma::Algebra::MinusSigmaXT)*in); - tmp += Ey*(Gamma(Gamma::Algebra::MinusSigmaYT)*in); - tmp += Ez*(Gamma(Gamma::Algebra::MinusSigmaZT)*in); - out += tmp*csw; // check signs - - } +//WilsonLoop::CloverPlaquette +///////////////////////////////////////////////////// +//// Clover plaquette combination in mu,nu plane with Double Stored U +//////////////////////////////////////////////////// +//static void CloverPlaquette(GaugeMat &Q, const std::vector &U, +// const int mu, const int nu){ +// Q = zero; +// Q += Gimpl::CovShiftBackward( +// U[mu], mu, Gimpl::CovShiftBackward( +// U[nu], nu, Gimpl::CovShiftForward(U[mu], mu, U[nu] ))); +// Q += Gimpl::CovShiftForward( +// U[mu], mu, Gimpl::CovShiftForward( +// U[nu], nu, Gimpl::CovShiftBackward(U[mu], mu, U[nu+Nd] ))); +// Q += Gimpl::CovShiftBackward( +// U[nu], nu, Gimpl::CovShiftForward( +// U[mu], mu, Gimpl::CovShiftForward(U[nu], nu, U[mu+Nd] ))); +// Q += Gimpl::CovShiftForward( +// U[mu], mu, Gimpl::CovShiftBackward( +// U[nu], nu, Gimpl::CovShiftBackward(U[mu], mu, U[nu] ))); +// } +// *NOT* EO template RealD WilsonCloverFermion::M(const FermionField& in, FermionField& out) { // Wilson term @@ -56,7 +63,7 @@ namespace QCD { this->Dhop(in, out, DaggerNo); // Clover term // apply the sigma and Fmunu - AddCloverTerm(in, out); + Mooee(in, out); // overall factor return axpy_norm(out, 4 + this->mass, in, out); } @@ -68,13 +75,16 @@ namespace QCD { this->Dhop(in, out, DaggerYes); // Clover term // apply the sigma and Fmunu - AddCloverTerm(in, out); + MooeeDag(in, out); return axpy_norm(out, 4 + this->mass, in, out); } template void WilsonCloverFermion::ImportGauge(const GaugeField& _Umu) { this->ImportGauge(_Umu); + GridBase* grid = _Umu._grid; + assert(Nd==4); //only works in 4 dim + typename Impl::GaugeLinkField Bx(grid), By(grid), Bz(grid), Ex(grid), Ey(grid), Ez(grid); // Compute the field strength terms WilsonLoops::FieldStrength(Bx, _Umu, Ydir, Zdir); WilsonLoops::FieldStrength(By, _Umu, Zdir, Xdir); @@ -82,31 +92,77 @@ namespace QCD { WilsonLoops::FieldStrength(Ex, _Umu, Tdir, Xdir); WilsonLoops::FieldStrength(Ey, _Umu, Tdir, Ydir); WilsonLoops::FieldStrength(Ez, _Umu, Tdir, Zdir); - // Save the contracted term with sigma - // into a dense matrix site by site - // Invert the Moo, Mee terms (using Eigen) + // Compute the Clover Operator acting on Colour and Spin + CloverTerm = fillClover(Bx)*(Gamma(Gamma::Algebra::SigmaYZ)); + CloverTerm += fillClover(By)*(Gamma(Gamma::Algebra::MinusSigmaXZ)); + CloverTerm += fillClover(Bz)*(Gamma(Gamma::Algebra::SigmaXY)); + CloverTerm += fillClover(Ex)*(Gamma(Gamma::Algebra::MinusSigmaXT)); + CloverTerm += fillClover(Ey)*(Gamma(Gamma::Algebra::MinusSigmaYT)); + CloverTerm += fillClover(Ez)*(Gamma(Gamma::Algebra::MinusSigmaZT)); + CloverTerm *= csw; + + + int lvol = _Umu._grid->lSites(); + int DimRep = Impl::Dimension; + + Eigen::MatrixXcd EigenCloverOp = Eigen::MatrixXcd::Zero(Ns*DimRep,Ns*DimRep); + Eigen::MatrixXcd EigenInvCloverOp = Eigen::MatrixXcd::Zero(Ns*DimRep,Ns*DimRep); + + std::vector lcoor; + typename SiteCloverType::scalar_object Qx = zero, Qxinv = zero; + + for (int site = 0; site < lvol; site++){ + grid->LocalIndexToLocalCoor(site,lcoor); + EigenCloverOp=Eigen::MatrixXcd::Zero(Ns*DimRep,Ns*DimRep); + peekLocalSite(Qx,CloverTerm,lcoor); + Qxinv = zero; + for(int j = 0; j < Ns; j++) + for (int k = 0; k < Ns; k++) + for(int a = 0; a < DimRep; a++) + for(int b = 0; b < DimRep; b++) + EigenCloverOp(a+j*DimRep,b+k*DimRep) = Qx()(j,k)(a,b); + + EigenInvCloverOp = EigenCloverOp.inverse(); + for(int j = 0; j < Ns; j++) + for (int k = 0; k < Ns; k++) + for(int a = 0; a < DimRep; a++) + for(int b = 0; b < DimRep; b++) + Qxinv()(j,k)(a,b) = EigenInvCloverOp(a+j*DimRep,b+k*DimRep); + + pokeLocalSite(Qxinv,CloverTermInv,lcoor); + } + } + + template + void WilsonCloverFermion::Mooee(const FermionField &in, FermionField &out){ + this -> MooeeInternal(in, out, DaggerNo, InverseNo); + } + + template + void WilsonCloverFermion::MooeeDag(const FermionField &in, FermionField &out){ + this -> MooeeInternal(in, out, DaggerNo, InverseYes); + } + + template + void WilsonCloverFermion::MooeeInv(const FermionField &in, FermionField &out){ + this -> MooeeInternal(in, out, DaggerNo, InverseYes); + } + + template + void WilsonCloverFermion::MooeeInvDag(const FermionField &in, FermionField &out){ + this -> MooeeInternal(in, out, DaggerNo, InverseYes); } template - void WilsonCloverFermion::Mooee(const FermionField &in, FermionField &out) { - out.checkerboard = in.checkerboard; - assert(0); // to be completed - } + void WilsonCloverFermion::MooeeInternal(const FermionField &in, FermionField &out, int dag, int inv){ + out.checkerboard = in.checkerboard; + CloverFieldType *Clover; - template - void WilsonCloverFermion::MooeeDag(const FermionField &in, FermionField &out) { - assert(0); // not implemented yet - } - template - void WilsonCloverFermion::MooeeInv(const FermionField &in, FermionField &out) { - assert(0); // not implemented yet - } - template - void WilsonCloverFermion::MooeeInvDag(const FermionField &in, FermionField &out) { - assert(0); // not implemented yet - } + Clover = (inv) ? &CloverTermInv : &CloverTerm; + if(dag){ out = adj(*Clover)*in;} else {out = *Clover*in;} + } // MooeeInternal // Derivative parts template @@ -128,17 +184,6 @@ namespace QCD { template void WilsonCloverFermion::MooDeriv(GaugeField&mat, const FermionField&U, const FermionField&V, int dag){ // Compute the 8 terms of the derivative - - // Pseudocode - // Using Chroma as a template - - // for loop on mu and nu, but upper matrix - // take the outer product factor * U x (sigma_mu_nu V) - - // derivative of loops - // end of loop - - assert(0); // not implemented yet } @@ -148,7 +193,10 @@ namespace QCD { assert(0); // not implemented yet } - FermOpTemplateInstantiate(WilsonCloverFermion); +FermOpTemplateInstantiate(WilsonCloverFermion); // now only for the fundamental representation +//AdjointFermOpTemplateInstantiate(WilsonCloverFermion); +//TwoIndexFermOpTemplateInstantiate(WilsonCloverFermion); +//GparityFermOpTemplateInstantiate(WilsonCloverFermion); } } diff --git a/lib/qcd/action/fermion/WilsonCloverFermion.h b/lib/qcd/action/fermion/WilsonCloverFermion.h index e942de1f..0fa0d57d 100644 --- a/lib/qcd/action/fermion/WilsonCloverFermion.h +++ b/lib/qcd/action/fermion/WilsonCloverFermion.h @@ -31,14 +31,20 @@ #include -namespace Grid { -namespace QCD { +namespace Grid +{ +namespace QCD +{ template -class WilsonCloverFermion : public WilsonFermion { +class WilsonCloverFermion : public WilsonFermion +{ public: + // Types definitions INHERIT_IMPL_TYPES(Impl); - + template using iImplClover = iScalar, Ns> >; + typedef iImplClover SiteCloverType; + typedef Lattice CloverFieldType; public: typedef WilsonFermion WilsonBase; @@ -51,43 +57,48 @@ public: const ImplParams &p = ImplParams()) : WilsonFermion(_Umu, Fgrid, Hgrid, - _mass, p), - Bx(_Umu._grid), - By(_Umu._grid), - Bz(_Umu._grid), - Ex(_Umu._grid), - Ey(_Umu._grid), - Ez(_Umu._grid) + _mass, p), + CloverTerm(&Fgrid), + CloverTermInv(&Fgrid) { csw = _csw; assert(Nd == 4); // require 4 dimensions } - virtual RealD M(const FermionField& in, FermionField& out); - virtual RealD Mdag(const FermionField& in, FermionField& out); + virtual RealD M(const FermionField &in, FermionField &out); + virtual RealD Mdag(const FermionField &in, FermionField &out); virtual void Mooee(const FermionField &in, FermionField &out); virtual void MooeeDag(const FermionField &in, FermionField &out); virtual void MooeeInv(const FermionField &in, FermionField &out); virtual void MooeeInvDag(const FermionField &in, FermionField &out); + virtual void MooeeInternal(const FermionField &in, FermionField &out, int dag, int inv); - virtual void MDeriv(GaugeField&mat, const FermionField&U, const FermionField&V, int dag); - virtual void MooDeriv(GaugeField&mat, const FermionField&U, const FermionField&V, int dag); - virtual void MeeDeriv(GaugeField&mat, const FermionField&U, const FermionField&V, int dag); - + virtual void MDeriv(GaugeField &mat, const FermionField &U, const FermionField &V, int dag); + virtual void MooDeriv(GaugeField &mat, const FermionField &U, const FermionField &V, int dag); + virtual void MeeDeriv(GaugeField &mat, const FermionField &U, const FermionField &V, int dag); void ImportGauge(const GaugeField &_Umu); + private: // here fixing the 4 dimensions, make it more general? - // Field strengths - GaugeLinkField Bx, By, Bz, Ex, Ey, Ez; + RealD csw; // Clover coefficient + CloverFieldType CloverTerm, CloverTermInv; // Clover term + // eventually these two can be compressed into 6x6 blocks instead of the 12x12 + // using the DeGrand-Rossi basis for the gamma matrices - RealD csw; // Clover coefficient - - - // Methods - void AddCloverTerm(const FermionField& in, FermionField& out); + CloverFieldType fillClover(const GaugeLinkField& F){ + CloverFieldType T(F._grid); + PARALLEL_FOR_LOOP + for (int i = 0; i < CloverTerm._grid->oSites(); i++){ + for (int s1 = 0; s1 < Nc; s1++) + for (int s2 = 0; s2 < Nc; s2++) + T._odata[i]()(s1,s2) = F._odata[i]()(); + } + return T; + } + }; } } diff --git a/tests/core/Test_wilson_clover.cc b/tests/core/Test_wilson_clover.cc new file mode 100644 index 00000000..08516d80 --- /dev/null +++ b/tests/core/Test_wilson_clover.cc @@ -0,0 +1,251 @@ + /************************************************************************************* + + Grid physics library, www.github.com/paboyle/Grid + + Source file: ./benchmarks/Benchmark_wilson.cc + + Copyright (C) 2015 + + Author: Guido Cossu + + This program is free software; you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation; either version 2 of the License, or + (at your option) any later version. + + This program is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License along + with this program; if not, write to the Free Software Foundation, Inc., + 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. + + See the full license in the file "LICENSE" in the top level distribution directory + *************************************************************************************/ + /* END LEGAL */ +#include + +using namespace std; +using namespace Grid; +using namespace Grid::QCD; + +int main (int argc, char ** argv) +{ + Grid_init(&argc,&argv); + + std::vector latt_size = GridDefaultLatt(); + std::vector simd_layout = GridDefaultSimd(Nd,vComplex::Nsimd()); + std::vector mpi_layout = GridDefaultMpi(); + GridCartesian Grid(latt_size,simd_layout,mpi_layout); + GridRedBlackCartesian RBGrid(latt_size,simd_layout,mpi_layout); + + int threads = GridThread::GetThreads(); + std::cout< seeds({1,2,3,4}); + GridParallelRNG pRNG(&Grid); + pRNG.SeedFixedIntegers(seeds); + // pRNG.SeedFixedIntegers(std::vector({45,12,81,9}); + + typedef typename WilsonCloverFermionR::FermionField FermionField; + typename WilsonCloverFermionR::ImplParams params; + + FermionField src (&Grid); random(pRNG,src); + FermionField result(&Grid); result=zero; + FermionField ref(&Grid); ref=zero; + FermionField tmp(&Grid); tmp=zero; + FermionField err(&Grid); tmp=zero; + FermionField phi (&Grid); random(pRNG,phi); + FermionField chi (&Grid); random(pRNG,chi); + LatticeGaugeField Umu(&Grid); SU3::HotConfiguration(pRNG,Umu); + std::vector U(4,&Grid); + + + double volume=1; + for(int mu=0;mu(Umu,mu); + /* Debug force unit + U[mu] = 1.0; + PokeIndex(Umu,U[mu],mu); + */ + } + + ref = zero; + + RealD mass=0.1; + RealD csw = 1.0; + + { // Simple clover implementation + + // ref = ref + mass * src; + } + + WilsonCloverFermionR Dwc(Umu,Grid,RBGrid,mass,csw,params); + + + std::cout< * = < chi | Deo^dag| phi> "< HermOpEO(Dwc); + HermOpEO.MpcDagMpc(chi_e,dchi_e,t1,t2); + HermOpEO.MpcDagMpc(chi_o,dchi_o,t1,t2); + + HermOpEO.MpcDagMpc(phi_e,dphi_e,t1,t2); + HermOpEO.MpcDagMpc(phi_o,dphi_o,t1,t2); + + pDce = innerProduct(phi_e,dchi_e); + pDco = innerProduct(phi_o,dchi_o); + cDpe = innerProduct(chi_e,dphi_e); + cDpo = innerProduct(chi_o,dphi_o); + + std::cout<