portelli/Grid
1
0
Fork 0
mirror of https://github.com/paboyle/Grid.git synced 2024-11-10 07:55:35 +00:00
Code Releases Activity

Merge pull request #133 from pretidav/feature/clover

Browse Source 
Feature/clover
This commit is contained in:
Guido Cossu 2017-10-24 15:15:21 +01:00 committed by GitHub
commit 0bc381f982
9 changed files with 837 additions and 134 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

1
.gitignore vendored
View File

@ -123,5 +123,6 @@ make-bin-BUCK.sh
lib/qcd/spin/gamma-gen/*.h
lib/qcd/spin/gamma-gen/*.cc
.vscode/
.vscode/settings.json
settings.json

3
.vscode/settings.json vendored
View File

@ -45,6 +45,7 @@
"istream": "cpp",
"ostream": "cpp",
"sstream": "cpp",
"streambuf": "cpp"
"streambuf": "cpp",
"algorithm": "cpp"
}
}

186
lib/qcd/action/fermion/WilsonCloverFermion.cc
View File

@ -39,30 +39,33 @@ namespace QCD
template <class Impl>
RealD WilsonCloverFermion<Impl>::M(const FermionField &in, FermionField &out)
{
FermionField temp(out._grid);
// Wilson term
out.checkerboard = in.checkerboard;
this->Dhop(in, out, DaggerNo);
// Clover term
// apply the sigma and Fmunu
FermionField temp(out._grid);
Mooee(in, temp);
// overall factor
out += temp;
return axpy_norm(out, 4 + this->mass, in, out);
return norm2(out);
}
template <class Impl>
RealD WilsonCloverFermion<Impl>::Mdag(const FermionField &in, FermionField &out)
{
FermionField temp(out._grid);
// Wilson term
out.checkerboard = in.checkerboard;
this->Dhop(in, out, DaggerYes);
// Clover term
// apply the sigma and Fmunu
FermionField temp(out._grid);
MooeeDag(in, temp);
out+=temp;
return axpy_norm(out, 4 + this->mass, in, out);
out += temp;
return norm2(out);
}
template <class Impl>
@ -72,22 +75,23 @@ void WilsonCloverFermion<Impl>::ImportGauge(const GaugeField &_Umu)
GridBase *grid = _Umu._grid;
typename Impl::GaugeLinkField Bx(grid), By(grid), Bz(grid), Ex(grid), Ey(grid), Ez(grid);
// Compute the field strength terms
WilsonLoops<Impl>::FieldStrength(Bx, _Umu, Ydir, Zdir);
// Compute the field strength terms mu>nu
WilsonLoops<Impl>::FieldStrength(Bx, _Umu, Zdir, Ydir);
WilsonLoops<Impl>::FieldStrength(By, _Umu, Zdir, Xdir);
WilsonLoops<Impl>::FieldStrength(Bz, _Umu, Xdir, Ydir);
WilsonLoops<Impl>::FieldStrength(Bz, _Umu, Ydir, Xdir);
WilsonLoops<Impl>::FieldStrength(Ex, _Umu, Tdir, Xdir);
WilsonLoops<Impl>::FieldStrength(Ey, _Umu, Tdir, Ydir);
WilsonLoops<Impl>::FieldStrength(Ez, _Umu, Tdir, Zdir);
// Compute the Clover Operator acting on Colour and Spin
CloverTerm = fillCloverYZ(Bx);
CloverTerm = fillCloverYZ(Bx);
CloverTerm += fillCloverXZ(By);
CloverTerm += fillCloverXY(Bz);
CloverTerm += fillCloverXT(Ex);
CloverTerm += fillCloverYT(Ey);
CloverTerm += fillCloverZT(Ez) ;
CloverTerm *= csw;
CloverTerm += fillCloverZT(Ez);
CloverTerm *= (0.5) * csw;
CloverTerm += (4.0 + this->mass);
int lvol = _Umu._grid->lSites();
int DimRep = Impl::Dimension;
@ -104,13 +108,13 @@ void WilsonCloverFermion<Impl>::ImportGauge(const GaugeField &_Umu)
EigenCloverOp = Eigen::MatrixXcd::Zero(Ns * DimRep, Ns * DimRep);
peekLocalSite(Qx, CloverTerm, lcoor);
Qxinv = zero;
//if (csw!=0){
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);
//std::cout << EigenCloverOp << std::endl;
// if (site==0) std::cout << "site =" << site << "\n" << EigenCloverOp << std::endl;
EigenInvCloverOp = EigenCloverOp.inverse();
//std::cout << EigenInvCloverOp << std::endl;
@ -119,30 +123,36 @@ void WilsonCloverFermion<Impl>::ImportGauge(const GaugeField &_Umu)
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);
// if (site==0) std::cout << "site =" << site << "\n" << EigenInvCloverOp << std::endl;
// }
pokeLocalSite(Qxinv, CloverTermInv, lcoor);
}
// Separate the even and odd parts.
// Separate the even and odd parts
pickCheckerboard(Even, CloverTermEven, CloverTerm);
pickCheckerboard( Odd, CloverTermOdd, CloverTerm);
pickCheckerboard(Odd, CloverTermOdd, CloverTerm);
pickCheckerboard(Even, CloverTermDagEven, adj(CloverTerm));
pickCheckerboard(Odd, CloverTermDagOdd, adj(CloverTerm));
pickCheckerboard(Even, CloverTermInvEven, CloverTermInv);
pickCheckerboard( Odd, CloverTermInvOdd, CloverTermInv);
pickCheckerboard(Odd, CloverTermInvOdd, CloverTermInv);
pickCheckerboard(Even, CloverTermInvDagEven, adj(CloverTermInv));
pickCheckerboard(Odd, CloverTermInvDagOdd, adj(CloverTermInv));
}
template <class Impl>
void WilsonCloverFermion<Impl>::Mooee(const FermionField &in, FermionField &out)
{
conformable(in,out);
conformable(in, out);
this->MooeeInternal(in, out, DaggerNo, InverseNo);
}
template <class Impl>
void WilsonCloverFermion<Impl>::MooeeDag(const FermionField &in, FermionField &out)
{
this->MooeeInternal(in, out, DaggerNo, InverseYes);
this->MooeeInternal(in, out, DaggerYes, InverseNo);
}
template <class Impl>
@ -154,7 +164,7 @@ void WilsonCloverFermion<Impl>::MooeeInv(const FermionField &in, FermionField &o
template <class Impl>
void WilsonCloverFermion<Impl>::MooeeInvDag(const FermionField &in, FermionField &out)
{
this->MooeeInternal(in, out, DaggerNo, InverseYes);
this->MooeeInternal(in, out, DaggerYes, InverseYes);
}
template <class Impl>
@ -164,32 +174,58 @@ void WilsonCloverFermion<Impl>::MooeeInternal(const FermionField &in, FermionFie
CloverFieldType *Clover;
assert(in.checkerboard == Odd || in.checkerboard == Even);
if (in._grid->_isCheckerBoarded)
if (dag)
{
if (in.checkerboard == Odd)
if (in._grid->_isCheckerBoarded)
{
std::cout << "Calling clover term Odd" << std::endl;
Clover = (inv) ? &CloverTermInvOdd : &CloverTermOdd;
if (in.checkerboard == Odd)
{
Clover = (inv) ? &CloverTermInvDagOdd : &CloverTermDagOdd;
}
else
{
Clover = (inv) ? &CloverTermInvDagEven : &CloverTermDagEven;
}
out = *Clover * in;
}
else
{
std::cout << "Calling clover term Even" << std::endl;
Clover = (inv) ? &CloverTermInvEven : &CloverTermEven;
Clover = (inv) ? &CloverTermInv : &CloverTerm;
out = adj(*Clover) * in;
}
}
else
{
Clover = (inv) ? &CloverTermInv : &CloverTerm;
if (in._grid->_isCheckerBoarded)
{
if (in.checkerboard == Odd)
{
// std::cout << "Calling clover term Odd" << std::endl;
Clover = (inv) ? &CloverTermInvOdd : &CloverTermOdd;
}
else
{
// std::cout << "Calling clover term Even" << std::endl;
Clover = (inv) ? &CloverTermInvEven : &CloverTermEven;
}
out = *Clover * in;
// std::cout << GridLogMessage << "*Clover.checkerboard " << (*Clover).checkerboard << std::endl;
}
else
{
Clover = (inv) ? &CloverTermInv : &CloverTerm;
out = *Clover * in;
}
}
std::cout << GridLogMessage << "*Clover.checkerboard " << (*Clover).checkerboard << std::endl;
if (dag){ out = adj(*Clover) * in;} else { out = *Clover * in;}
} // MooeeInternal
// Derivative parts
template <class Impl>
void WilsonCloverFermion<Impl>::MDeriv(GaugeField &mat, const FermionField &U, const FermionField &V, int dag)
{
GaugeField tmp(mat._grid);
conformable(U._grid, V._grid);
@ -205,10 +241,90 @@ void WilsonCloverFermion<Impl>::MDeriv(GaugeField &mat, const FermionField &U, c
// Derivative parts
template <class Impl>
void WilsonCloverFermion<Impl>::MooDeriv(GaugeField &mat, const FermionField &U, const FermionField &V, int dag)
void WilsonCloverFermion<Impl>::MooDeriv(GaugeField &mat, const FermionField &X, const FermionField &Y, int dag)
{
// Compute the 8 terms of the derivative
assert(0); // not implemented yet
GridBase *grid = mat._grid;
//GaugeLinkField Lambdaodd(grid), Lambdaeven(grid), tmp(grid);
//Lambdaodd = zero; //Yodd*dag(Xodd)+Xodd*dag(Yodd); // I have to peek spin and decide the color structure
//Lambdaeven = zero; //Teven*dag(Xeven)+Xeven*dag(Yeven) + 2*(Dee^-1)
GaugeLinkField Lambda(grid), tmp(grid);
Lambda = zero;
conformable(mat._grid, X._grid);
conformable(Y._grid, X._grid);
std::vector<GaugeLinkField> C1p(Nd, grid), C2p(Nd, grid), C3p(Nd, grid), C4p(Nd, grid);
std::vector<GaugeLinkField> C1m(Nd, grid), C2m(Nd, grid), C3m(Nd, grid), C4m(Nd, grid);
std::vector<GaugeLinkField> U(Nd, mat._grid);
for (int mu = 0; mu < Nd; mu++)
{
U[mu] = PeekIndex<LorentzIndex>(mat, mu);
C1p[mu] = zero;
C2p[mu] = zero;
C3p[mu] = zero;
C4p[mu] = zero;
C1m[mu] = zero;
C2m[mu] = zero;
C3m[mu] = zero;
C4m[mu] = zero;
}
/*
PARALLEL_FOR_LOOP
for (int i = 0; i < CloverTerm._grid->oSites(); i++)
{
T._odata[i]()(0, 1) = timesMinusI(F._odata[i]()());
T._odata[i]()(1, 0) = timesMinusI(F._odata[i]()());
T._odata[i]()(2, 3) = timesMinusI(F._odata[i]()());
T._odata[i]()(3, 2) = timesMinusI(F._odata[i]()());
}
*/
for (int i = 0; i < 4; i++)
{ //spin
for (int j = 0; j < 4; j++)
{ //spin
for (int mu = 0; mu < 4; mu++)
{ //color
for (int nu = 0; nu < 4; nu++)
{ //color
// insertion in upper staple
tmp = Lambda * U[nu];
C1p[mu] += Impl::ShiftStaple(Impl::CovShiftForward(tmp, nu, Impl::CovShiftBackward(U[mu], mu, Impl::CovShiftIdentityBackward(U[nu], nu))), mu);
tmp = Lambda * U[mu];
C2p[mu] += Impl::ShiftStaple(Impl::CovShiftForward(U[nu], nu, Impl::CovShiftBackward(tmp, mu, Impl::CovShiftIdentityBackward(U[nu], nu))), mu);
tmp = Impl::CovShiftIdentityForward(Lambda, nu) * U[nu];
C3p[mu] += Impl::ShiftStaple(Impl::CovShiftForward(U[nu], nu, Impl::CovShiftBackward(U[mu], mu, Impl::CovShiftIdentityBackward(tmp, nu))), mu);
tmp = Lambda;
C4p[mu] += Impl::ShiftStaple(Impl::CovShiftForward(U[nu], nu, Impl::CovShiftBackward(U[mu], mu, Impl::CovShiftIdentityBackward(U[nu], nu))), mu) * tmp;
// insertion in lower staple
tmp = Lambda * U[nu];
C1m[mu] += Impl::ShiftStaple(Impl::CovShiftBackward(tmp, nu, Impl::CovShiftBackward(U[mu], mu, U[nu])), mu);
tmp = Lambda * U[mu];
C2m[mu] += Impl::ShiftStaple(Impl::CovShiftBackward(U[nu], nu, Impl::CovShiftBackward(tmp, mu, U[nu])), mu);
tmp = Lambda * U[nu];
C3m[mu] += Impl::ShiftStaple(Impl::CovShiftBackward(U[nu], nu, Impl::CovShiftBackward(U[mu], mu, tmp)), mu);
tmp = Lambda;
C4m[mu] += Impl::ShiftStaple(Impl::CovShiftBackward(U[nu], nu, Impl::CovShiftBackward(U[mu], mu, U[nu])), mu) * tmp;
}
}
}
}
//Still implementing. Have to be tested, and understood how to project EO
}
// Derivative parts

99
lib/qcd/action/fermion/WilsonCloverFermion.h
View File

@ -26,6 +26,7 @@
See the full license in the file "LICENSE" in the top level distribution directory
*************************************************************************************/
/* END LEGAL */
#ifndef GRID_QCD_WILSON_CLOVER_FERMION_H
#define GRID_QCD_WILSON_CLOVER_FERMION_H
@ -42,9 +43,11 @@ class WilsonCloverFermion : public WilsonFermion<Impl>
public:
// Types definitions
INHERIT_IMPL_TYPES(Impl);
template <typename vtype> using iImplClover = iScalar<iMatrix<iMatrix<vtype, Impl::Dimension>, Ns> >;
typedef iImplClover<Simd> SiteCloverType;
typedef Lattice<SiteCloverType> CloverFieldType;
template <typename vtype>
using iImplClover = iScalar<iMatrix<iMatrix<vtype, Impl::Dimension>, Ns>>;
typedef iImplClover<Simd> SiteCloverType;
typedef Lattice<SiteCloverType> CloverFieldType;
public:
typedef WilsonFermion<Impl> WilsonBase;
@ -58,15 +61,21 @@ public:
Fgrid,
Hgrid,
_mass, p),
CloverTerm(&Fgrid),
CloverTermInv(&Fgrid),
CloverTermEven(&Hgrid),
CloverTermOdd(&Hgrid),
CloverTermInvEven(&Hgrid),
CloverTermInvOdd(&Hgrid)
CloverTerm(&Fgrid),
CloverTermInv(&Fgrid),
CloverTermEven(&Hgrid),
CloverTermOdd(&Hgrid),
CloverTermInvEven(&Hgrid),
CloverTermInvOdd(&Hgrid),
CloverTermDagEven(&Hgrid),
CloverTermDagOdd(&Hgrid),
CloverTermInvDagEven(&Hgrid),
CloverTermInvDagOdd(&Hgrid)
{
csw = _csw;
assert(Nd == 4); // require 4 dimensions
if (csw == 0) std::cout << GridLogWarning << "Initializing WilsonCloverFermion with csw = 0" << std::endl;
}
virtual RealD M(const FermionField &in, FermionField &out);
@ -87,10 +96,13 @@ public:
private:
// here fixing the 4 dimensions, make it more general?
RealD csw; // Clover coefficient
CloverFieldType CloverTerm, CloverTermInv; // Clover term
CloverFieldType CloverTermEven, CloverTermOdd;
CloverFieldType CloverTermInvEven, CloverTermInvOdd; // Clover term
RealD csw; // Clover coefficient
CloverFieldType CloverTerm, CloverTermInv; // Clover term
CloverFieldType CloverTermEven, CloverTermOdd; // Clover term EO
CloverFieldType CloverTermInvEven, CloverTermInvOdd; // Clover term Inv EO
CloverFieldType CloverTermDagEven, CloverTermDagOdd; // Clover term Dag EO
CloverFieldType CloverTermInvDagEven, CloverTermInvDagOdd; // Clover term Inv Dag EO
// eventually these two can be compressed into 6x6 blocks instead of the 12x12
// using the DeGrand-Rossi basis for the gamma matrices
@ -106,9 +118,9 @@ private:
T._odata[i]()(2, 3) = timesMinusI(F._odata[i]()());
T._odata[i]()(3, 2) = timesMinusI(F._odata[i]()());
}
return T;
}
return T;
}
CloverFieldType fillCloverXZ(const GaugeLinkField &F)
{
@ -122,9 +134,9 @@ private:
T._odata[i]()(2, 3) = -F._odata[i]()();
T._odata[i]()(3, 2) = F._odata[i]()();
}
return T;
}
return T;
}
CloverFieldType fillCloverXY(const GaugeLinkField &F)
{
@ -138,9 +150,9 @@ private:
T._odata[i]()(2, 2) = timesMinusI(F._odata[i]()());
T._odata[i]()(3, 3) = timesI(F._odata[i]()());
}
return T;
}
return T;
}
CloverFieldType fillCloverXT(const GaugeLinkField &F)
{
@ -149,14 +161,14 @@ private:
PARALLEL_FOR_LOOP
for (int i = 0; i < CloverTerm._grid->oSites(); i++)
{
T._odata[i]()(0, 1) = timesMinusI(F._odata[i]()());
T._odata[i]()(1, 0) = timesMinusI(F._odata[i]()());
T._odata[i]()(2, 3) = timesI(F._odata[i]()());
T._odata[i]()(3, 2) = timesI(F._odata[i]()());
T._odata[i]()(0, 1) = timesI(F._odata[i]()());
T._odata[i]()(1, 0) = timesI(F._odata[i]()());
T._odata[i]()(2, 3) = timesMinusI(F._odata[i]()());
T._odata[i]()(3, 2) = timesMinusI(F._odata[i]()());
}
return T;
}
return T;
}
CloverFieldType fillCloverYT(const GaugeLinkField &F)
{
@ -165,14 +177,14 @@ private:
PARALLEL_FOR_LOOP
for (int i = 0; i < CloverTerm._grid->oSites(); i++)
{
T._odata[i]()(0, 1) = (F._odata[i]()());
T._odata[i]()(1, 0) = -(F._odata[i]()());
T._odata[i]()(2, 3) = -(F._odata[i]()());
T._odata[i]()(3, 2) = (F._odata[i]()());
T._odata[i]()(0, 1) = -(F._odata[i]()());
T._odata[i]()(1, 0) = (F._odata[i]()());
T._odata[i]()(2, 3) = (F._odata[i]()());
T._odata[i]()(3, 2) = -(F._odata[i]()());
}
return T;
}
return T;
}
CloverFieldType fillCloverZT(const GaugeLinkField &F)
{
@ -181,17 +193,16 @@ private:
PARALLEL_FOR_LOOP
for (int i = 0; i < CloverTerm._grid->oSites(); i++)
{
T._odata[i]()(0, 0) = timesMinusI(F._odata[i]()());
T._odata[i]()(1, 1) = timesI(F._odata[i]()());
T._odata[i]()(2, 2) = timesI(F._odata[i]()());
T._odata[i]()(3, 3) = timesMinusI(F._odata[i]()());
T._odata[i]()(0, 0) = timesI(F._odata[i]()());
T._odata[i]()(1, 1) = timesMinusI(F._odata[i]()());
T._odata[i]()(2, 2) = timesMinusI(F._odata[i]()());
T._odata[i]()(3, 3) = timesI(F._odata[i]()());
}
return T;
}
return T;
}
};
}
}
#endif // GRID_QCD_WILSON_CLOVER_FERMION_H
#endif // GRID_QCD_WILSON_CLOVER_FERMION_H

4
lib/qcd/utils/WilsonLoops.h
View File

@ -337,7 +337,9 @@ static void StapleMult(GaugeMat &staple, const GaugeLorentz &Umu, int mu) {
GaugeMat v = Vup - Vdn;
GaugeMat u = PeekIndex<LorentzIndex>(Umu, mu); // some redundant copies
GaugeMat vu = v*u;
FS = 0.25*Ta(u*v + Cshift(vu, mu, -1));
//FS = 0.25*Ta(u*v + Cshift(vu, mu, -1));
FS = (u*v + Cshift(vu, mu, -1));
FS = 0.125*(FS - adj(FS));
}
static Real TopologicalCharge(GaugeLorentz &U){

176
tests/core/Test_wilson_clover.cc
View File

@ -55,13 +55,15 @@ int main (int argc, char ** argv)
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);
FermionField src (&Grid); random(pRNG,src);
FermionField result(&Grid); result=zero;
FermionField result2(&Grid); result2=zero;
FermionField ref(&Grid); ref=zero;
FermionField tmp(&Grid); tmp=zero;
FermionField err(&Grid); err=zero;
FermionField err2(&Grid); err2=zero;
FermionField phi (&Grid); random(pRNG,phi);
FermionField chi (&Grid); random(pRNG,chi);
LatticeGaugeField Umu(&Grid); SU3::HotConfiguration(pRNG,Umu);
std::vector<LatticeColourMatrix> U(4,&Grid);
@ -71,24 +73,9 @@ int main (int argc, char ** argv)
volume=volume*latt_size[mu];
}
// Only one non-zero (y)
for(int mu=0;mu<Nd;mu++){
U[mu] = PeekIndex<LorentzIndex>(Umu,mu);
/* Debug force unit
U[mu] = 1.0;
PokeIndex<LorentzIndex>(Umu,U[mu],mu);
*/
}
ref = zero;
RealD mass=0.1;
RealD mass= 0.1;
RealD csw = 1.0;
{ // Simple clover implementation
// ref = ref + mass * src;
}
WilsonCloverFermionR Dwc(Umu,Grid,RBGrid,mass,csw,params);
Dwc.ImportGauge(Umu);
@ -172,35 +159,30 @@ int main (int argc, char ** argv)
std::cout<<GridLogMessage <<"e "<<pDce<<" "<<cDpe <<std::endl;
std::cout<<GridLogMessage <<"o "<<pDco<<" "<<cDpo <<std::endl;
std::cout<<GridLogMessage <<"pDce - conj(cDpo) "<< pDce-conj(cDpo) <<std::endl;
std::cout<<GridLogMessage <<"pDco - conj(cDpe) "<< pDco-conj(cDpe) <<std::endl;
std::cout<<GridLogMessage <<"e "<<pDce<<" "<<cDpe <<std::endl;
std::cout<<GridLogMessage <<"o "<<pDco<<" "<<cDpo <<std::endl;
std::cout<<GridLogMessage <<"pDce - conj(cDpo) "<< pDce-conj(cDpo) <<std::endl;
std::cout<<GridLogMessage <<"pDco - conj(cDpe) "<< pDco-conj(cDpe) <<std::endl;
std::cout<<GridLogMessage<<"=============================================================="<<std::endl;
std::cout<<GridLogMessage<<"= Test MeeInv Mee = 1 "<<std::endl;
std::cout<<GridLogMessage<<"= Test MeeInv Mee = 1 (if csw!=0) "<<std::endl;
std::cout<<GridLogMessage<<"=============================================================="<<std::endl;
pickCheckerboard(Even,chi_e,chi);
pickCheckerboard(Odd ,chi_o,chi);
Dwc.Mooee(chi_e,src_e);
Dwc.MooeeInv(src_e,phi_e);
Dwc.MooeeInv(src_e,phi_e);
Dwc.Mooee(chi_o,src_o);
Dwc.MooeeInv(src_o,phi_o);
Dwc.MooeeInv(src_o,phi_o);
setCheckerboard(phi,phi_e);
setCheckerboard(phi,phi_o);
err = phi-chi;
std::cout<<GridLogMessage << "norm diff "<< norm2(err)<< std::endl;
std::cout<<GridLogMessage<<"=============================================================="<<std::endl;
std::cout<<GridLogMessage<<"= Test MeeInvDag MeeDag = 1 "<<std::endl;
std::cout<<GridLogMessage<<"= Test MeeDag MeeInvDag = 1 (if csw!=0) "<<std::endl;
std::cout<<GridLogMessage<<"=============================================================="<<std::endl;
pickCheckerboard(Even,chi_e,chi);
@ -219,33 +201,125 @@ int main (int argc, char ** argv)
std::cout<<GridLogMessage << "norm diff "<< norm2(err)<< std::endl;
std::cout<<GridLogMessage<<"=============================================================="<<std::endl;
std::cout<<GridLogMessage<<"= Test MpcDagMpc is Hermitian "<<std::endl;
std::cout<<GridLogMessage<<"= Test MeeInv MeeDag = 1 (if csw!=0) "<<std::endl;
std::cout<<GridLogMessage<<"=============================================================="<<std::endl;
random(pRNG,phi);
random(pRNG,chi);
pickCheckerboard(Even,chi_e,chi);
pickCheckerboard(Odd ,chi_o,chi);
Dwc.MooeeDag(chi_e,src_e);
Dwc.MooeeInv(src_e,phi_e);
Dwc.MooeeDag(chi_o,src_o);
Dwc.MooeeInv(src_o,phi_o);
setCheckerboard(phi,phi_e);
setCheckerboard(phi,phi_o);
err = phi-chi;
std::cout<<GridLogMessage << "norm diff "<< norm2(err)<< std::endl;
std::cout<<GridLogMessage<<"================================================================"<<std::endl;
std::cout<<GridLogMessage<<"= Testing gauge covariance Clover term with EO preconditioning "<<std::endl;
std::cout<<GridLogMessage<<"================================================================"<<std::endl;
chi=zero; phi=zero; tmp=zero;
pickCheckerboard(Even,chi_e,chi);
pickCheckerboard(Odd ,chi_o,chi);
pickCheckerboard(Even,phi_e,phi);
pickCheckerboard(Odd ,phi_o,phi);
SchurDiagMooeeOperator<WilsonCloverFermionR,FermionField> HermOpEO(Dwc);
HermOpEO.MpcDagMpc(chi_e,dchi_e,t1,t2);
HermOpEO.MpcDagMpc(chi_o,dchi_o,t1,t2);
Dwc.Mooee(src_e,chi_e);
Dwc.Mooee(src_o,chi_o);
setCheckerboard(chi,chi_e);
setCheckerboard(chi,chi_o);
setCheckerboard(src,src_e);
setCheckerboard(src,src_o);
HermOpEO.MpcDagMpc(phi_e,dphi_e,t1,t2);
HermOpEO.MpcDagMpc(phi_o,dphi_o,t1,t2);
////////////////////// Gauge Transformation
std::vector<int> seeds2({5,6,7,8});
GridParallelRNG pRNG2(&Grid); pRNG2.SeedFixedIntegers(seeds2);
LatticeColourMatrix Omega(&Grid);
LatticeColourMatrix ShiftedOmega(&Grid);
LatticeGaugeField U_prime(&Grid); U_prime=zero;
LatticeColourMatrix U_prime_mu(&Grid); U_prime_mu=zero;
SU<Nc>::LieRandomize(pRNG2, Omega, 1.0);
for (int mu=0;mu<Nd;mu++){
U[mu]=peekLorentz(Umu,mu);
ShiftedOmega=Cshift(Omega,mu,1);
U_prime_mu=Omega*U[mu]*adj(ShiftedOmega);
pokeLorentz(U_prime,U_prime_mu,mu);
}
/////////////////
WilsonCloverFermionR Dwc_prime(U_prime,Grid,RBGrid,mass,csw,params);
Dwc_prime.ImportGauge(U_prime);
pDce = innerProduct(phi_e,dchi_e);
pDco = innerProduct(phi_o,dchi_o);
cDpe = innerProduct(chi_e,dphi_e);
cDpo = innerProduct(chi_o,dphi_o);
tmp=Omega*src;
pickCheckerboard(Even,src_e,tmp);
pickCheckerboard(Odd ,src_o,tmp);
std::cout<<GridLogMessage <<"e "<<pDce<<" "<<cDpe <<std::endl;
std::cout<<GridLogMessage <<"o "<<pDco<<" "<<cDpo <<std::endl;
Dwc_prime.Mooee(src_e,phi_e);
Dwc_prime.Mooee(src_o,phi_o);
std::cout<<GridLogMessage <<"pDce - conj(cDpo) "<< pDco-conj(cDpo) <<std::endl;
std::cout<<GridLogMessage <<"pDco - conj(cDpe) "<< pDce-conj(cDpe) <<std::endl;
setCheckerboard(phi,phi_e);
setCheckerboard(phi,phi_o);
Grid_finalize();
err = chi - adj(Omega)*phi;
std::cout<<GridLogMessage << "norm diff "<< norm2(err)<< std::endl;
std::cout<<GridLogMessage<<"================================================================="<<std::endl;
std::cout<<GridLogMessage<<"= Testing gauge covariance Clover term w/o EO preconditioning "<<std::endl;
std::cout<<GridLogMessage<<"================================================================"<<std::endl;
chi=zero; phi=zero;
WilsonFermionR Dw(Umu,Grid,RBGrid,mass,params);
Dw.ImportGauge(Umu);
Dw.M(src,result);
Dwc.M(src,chi);
Dwc_prime.M(Omega*src,phi);
WilsonFermionR Dw_prime(U_prime,Grid,RBGrid,mass,params);
Dw_prime.ImportGauge(U_prime);
Dw_prime.M(Omega*src,result2);
err = chi-adj(Omega)*phi;
err2 = result-adj(Omega)*result2;
std::cout<<GridLogMessage << "norm diff Wilson "<< norm2(err)<< std::endl;
std::cout<<GridLogMessage << "norm diff WilsonClover "<< norm2(err2)<< std::endl;
std::cout<<GridLogMessage<<"=========================================================="<<std::endl;
std::cout<<GridLogMessage<<"= Testing Mooee(csw=0) Clover to reproduce Mooee Wilson "<<std::endl;
std::cout<<GridLogMessage<<"=========================================================="<<std::endl;
chi=zero; phi=zero; err=zero;
WilsonCloverFermionR Dwc_csw0(Umu,Grid,RBGrid,mass,0.0,params); // <-- Notice: csw=0
Dwc_csw0.ImportGauge(Umu);
pickCheckerboard(Even,phi_e,phi);
pickCheckerboard(Odd ,phi_o,phi);
pickCheckerboard(Even,chi_e,chi);
pickCheckerboard(Odd ,chi_o,chi);
Dw.Mooee(src_e,chi_e);
Dw.Mooee(src_o,chi_o);
Dwc_csw0.Mooee(src_e,phi_e);
Dwc_csw0.Mooee(src_o,phi_o);
setCheckerboard(chi,chi_e);
setCheckerboard(chi,chi_o);
setCheckerboard(phi,phi_e);
setCheckerboard(phi,phi_o);
setCheckerboard(src,src_e);
setCheckerboard(src,src_o);
err = chi - phi;
std::cout<<GridLogMessage << "norm diff "<< norm2(err)<< std::endl;
Grid_finalize();
}

3
tests/qdpxx/Makefile.am
View File

@ -1,4 +1,5 @@
AM_CXXFLAGS += `chroma-config --cxxflags`
AM_LDFLAGS += `chroma-config --ldflags` `chroma-config --libs`
AM_LDFLAGS += `chroma-config --ldflags`
LIBS += `chroma-config --libs`
include Make.inc

3
tests/qdpxx/Makefile.am.qdpxx
View File

@ -1,6 +1,7 @@
# additional include paths necessary to compile the C++ library
AM_CXXFLAGS = -I$(top_srcdir)/include `chroma-config --cxxflags`
AM_LDFLAGS = -L$(top_builddir)/lib `chroma-config --ldflags` `chroma-config --libs`
AM_LDFLAGS = -L$(top_builddir)/lib `chroma-config --ldflags`
AM_LIBS = `chroma-config --libs`
include Make.inc

496
tests/qdpxx/Test_qdpxx_wilson.cc Normal file
View File

@ -0,0 +1,496 @@
/*************************************************************************************
Grid physics library, www.github.com/paboyle/Grid
Source file: ./tests/qdpxx/Test_qdpxx_wilson.cc
Copyright (C) 2017
Author: Guido Cossu <guido.cossu@ed.ac.uk>
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 <Grid/Grid.h>
#include <chroma.h>
#include <actions/ferm/invert/syssolver_linop_cg_array.h>
#include <actions/ferm/invert/syssolver_linop_aggregate.h>
// Mass
double mq = 0.01;
// Define Wilson Types
typedef Grid::QCD::WilsonImplR::FermionField FermionField;
typedef Grid::QCD::LatticeGaugeField GaugeField;
enum ChromaAction
{
Wilson, // Wilson
WilsonClover // CloverFermions
};
namespace Chroma
{
class ChromaWrapper
{
public:
typedef multi1d<LatticeColorMatrix> U;
typedef LatticeFermion T4;
static void ImportGauge(GaugeField &gr,
QDP::multi1d<QDP::LatticeColorMatrix> &ch)
{
Grid::QCD::LorentzColourMatrix LCM;
Grid::Complex cc;
QDP::ColorMatrix cm;
QDP::Complex c;
std::vector<int> x(4);
QDP::multi1d<int> cx(4);
std::vector<int> gd = gr._grid->GlobalDimensions();
for (x[0] = 0; x[0] < gd[0]; x[0]++)
{
for (x[1] = 0; x[1] < gd[1]; x[1]++)
{
for (x[2] = 0; x[2] < gd[2]; x[2]++)
{
for (x[3] = 0; x[3] < gd[3]; x[3]++)
{
cx[0] = x[0];
cx[1] = x[1];
cx[2] = x[2];
cx[3] = x[3];
Grid::peekSite(LCM, gr, x);
for (int mu = 0; mu < 4; mu++)
{
for (int i = 0; i < 3; i++)
{
for (int j = 0; j < 3; j++)
{
cc = LCM(mu)()(i, j);
c = QDP::cmplx(QDP::Real(real(cc)), QDP::Real(imag(cc)));
QDP::pokeColor(cm, c, i, j);
}
}
QDP::pokeSite(ch[mu], cm, cx);
}
}
}
}
}
}
static void ExportGauge(GaugeField &gr,
QDP::multi1d<QDP::LatticeColorMatrix> &ch)
{
Grid::QCD::LorentzColourMatrix LCM;
Grid::Complex cc;
QDP::ColorMatrix cm;
QDP::Complex c;
std::vector<int> x(4);
QDP::multi1d<int> cx(4);
std::vector<int> gd = gr._grid->GlobalDimensions();
for (x[0] = 0; x[0] < gd[0]; x[0]++)
{
for (x[1] = 0; x[1] < gd[1]; x[1]++)
{
for (x[2] = 0; x[2] < gd[2]; x[2]++)
{
for (x[3] = 0; x[3] < gd[3]; x[3]++)
{
cx[0] = x[0];
cx[1] = x[1];
cx[2] = x[2];
cx[3] = x[3];
for (int mu = 0; mu < 4; mu++)
{
for (int i = 0; i < 3; i++)
{
for (int j = 0; j < 3; j++)
{
cm = QDP::peekSite(ch[mu], cx);
c = QDP::peekColor(cm, i, j);
cc = Grid::Complex(toDouble(real(c)), toDouble(imag(c)));
LCM(mu)
()(i, j) = cc;
}
}
}
Grid::pokeSite(LCM, gr, x);
}
}
}
}
}
// Specific for Wilson Fermions
static void ImportFermion(Grid::QCD::LatticeFermion &gr,
QDP::LatticeFermion &ch)
{
Grid::QCD::SpinColourVector F;
Grid::Complex c;
QDP::Fermion cF;
QDP::SpinVector cS;
QDP::Complex cc;
std::vector<int> x(4); // explicit 4d fermions in Grid
QDP::multi1d<int> cx(4);
std::vector<int> gd = gr._grid->GlobalDimensions();
for (x[0] = 0; x[0] < gd[0]; x[0]++)
{
for (x[1] = 0; x[1] < gd[1]; x[1]++)
{
for (x[2] = 0; x[2] < gd[2]; x[2]++)
{
for (x[3] = 0; x[3] < gd[3]; x[3]++)
{
cx[0] = x[0];
cx[1] = x[1];
cx[2] = x[2];
cx[3] = x[3];
Grid::peekSite(F, gr, x);
for (int j = 0; j < 3; j++)
{
for (int sp = 0; sp < 4; sp++)
{
c = F()(sp)(j);
cc = QDP::cmplx(QDP::Real(real(c)), QDP::Real(imag(c)));
QDP::pokeSpin(cS, cc, sp);
}
QDP::pokeColor(cF, cS, j);
}
QDP::pokeSite(ch, cF, cx);
}
}
}
}
}
// Specific for 4d Wilson fermions
static void ExportFermion(Grid::QCD::LatticeFermion &gr,
QDP::LatticeFermion &ch)
{
Grid::QCD::SpinColourVector F;
Grid::Complex c;
QDP::Fermion cF;
QDP::SpinVector cS;
QDP::Complex cc;
std::vector<int> x(4); // 4d fermions
QDP::multi1d<int> cx(4);
std::vector<int> gd = gr._grid->GlobalDimensions();
for (x[0] = 0; x[0] < gd[0]; x[0]++)
{
for (x[1] = 0; x[1] < gd[1]; x[1]++)
{
for (x[2] = 0; x[2] < gd[2]; x[2]++)
{
for (x[3] = 0; x[3] < gd[3]; x[3]++)
{
cx[0] = x[0];
cx[1] = x[1];
cx[2] = x[2];
cx[3] = x[3];
cF = QDP::peekSite(ch, cx);
for (int sp = 0; sp < 4; sp++)
{
for (int j = 0; j < 3; j++)
{
cS = QDP::peekColor(cF, j);
cc = QDP::peekSpin(cS, sp);
c = Grid::Complex(QDP::toDouble(QDP::real(cc)),
QDP::toDouble(QDP::imag(cc)));
F()
(sp)(j) = c;
}
}
Grid::pokeSite(F, gr, x);
}
}
}
}
}
static Handle<Chroma::UnprecLinearOperator<T4, U, U>> GetLinOp(U &u, ChromaAction params)
{
QDP::Real _mq(mq);
QDP::multi1d<int> bcs(QDP::Nd);
// Boundary conditions
bcs[0] = bcs[1] = bcs[2] = bcs[3] = 1;
if (params == Wilson)
{
Chroma::WilsonFermActParams p;
p.Mass = _mq;
Chroma::Handle<Chroma::FermBC<T4, U, U>> fbc(new Chroma::SimpleFermBC<T4, U, U>(bcs));
Chroma::Handle<Chroma::CreateFermState<T4, U, U>> cfs(new Chroma::CreateSimpleFermState<T4, U, U>(fbc));
Chroma::UnprecWilsonFermAct S_f(cfs, p);
Chroma::Handle<Chroma::FermState<T4, U, U>> ffs(S_f.createState(u));
return S_f.linOp(ffs);
}
if (params == WilsonClover)
{
Chroma::CloverFermActParams p;
p.Mass = _mq;
p.clovCoeffR = QDP::Real(1.0);
p.clovCoeffT = QDP::Real(1.0);
Real u0 = QDP::Real(1.0);
Chroma::Handle<Chroma::FermBC<T4, U, U>> fbc(new Chroma::SimpleFermBC<T4, U, U>(bcs));
Chroma::Handle<Chroma::CreateFermState<T4, U, U>> cfs(new Chroma::CreateSimpleFermState<T4, U, U>(fbc));
Chroma::UnprecCloverFermAct S_f(cfs, p);
Chroma::Handle<Chroma::FermState<T4, U, U>> ffs(S_f.createState(u));
return S_f.linOp(ffs);
}
}
};
} // namespace Chroma
void calc_chroma(ChromaAction action, GaugeField &lat, FermionField &src, FermionField &res, int dag)
{
QDP::multi1d<QDP::LatticeColorMatrix> u(4);
Chroma::ChromaWrapper::ImportGauge(lat, u);
QDP::LatticeFermion check;
QDP::LatticeFermion result;
QDP::LatticeFermion psi;
Chroma::ChromaWrapper::ImportFermion(src, psi);
for (int mu = 0; mu < 4; mu++)
{
std::cout << "Imported Gauge norm [" << mu << "] " << QDP::norm2(u[mu]) << std::endl;
}
std::cout << "Imported Fermion norm " << QDP::norm2(psi) << std::endl;
typedef QDP::LatticeFermion T;
typedef QDP::multi1d<QDP::LatticeColorMatrix> U;
auto linop = Chroma::ChromaWrapper::GetLinOp(u, action);
printf("Calling Chroma Linop\n");
fflush(stdout);
if (dag)
(*linop)(check, psi, Chroma::MINUS);
else
(*linop)(check, psi, Chroma::PLUS);
printf("Called Chroma Linop\n");
fflush(stdout);
// std::cout << "Calling Chroma Linop " << std::endl;
// linop->evenEvenLinOp(tmp, psi, isign);
// check[rb[0]] = tmp;
// linop->oddOddLinOp(tmp, psi, isign);
// check[rb[1]] = tmp;
// linop->evenOddLinOp(tmp, psi, isign);
// check[rb[0]] += tmp;
// linop->oddEvenLinOp(tmp, psi, isign);
// check[rb[1]] += tmp;
Chroma::ChromaWrapper::ExportFermion(res, check);
}
void make_gauge(GaugeField &Umu, FermionField &src)
{
using namespace Grid;
using namespace Grid::QCD;
std::vector<int> seeds4({1, 2, 3, 4});
Grid::GridCartesian *UGrid = (Grid::GridCartesian *)Umu._grid;
Grid::GridParallelRNG RNG4(UGrid);
RNG4.SeedFixedIntegers(seeds4);
Grid::QCD::SU3::HotConfiguration(RNG4, Umu);
// Fermion field
Grid::gaussian(RNG4, src);
/*
Grid::QCD::SpinColourVector F;
Grid::Complex c;
std::vector<int> x(4); // 4d fermions
std::vector<int> gd = src._grid->GlobalDimensions();
for (x[0] = 0; x[0] < gd[0]; x[0]++)
{
for (x[1] = 0; x[1] < gd[1]; x[1]++)
{
for (x[2] = 0; x[2] < gd[2]; x[2]++)
{
for (x[3] = 0; x[3] < gd[3]; x[3]++)
{
for (int sp = 0; sp < 4; sp++)
{
for (int j = 0; j < 3; j++) // colours
{
F()(sp)(j) = Grid::Complex(0.0,0.0);
if (((sp == 0)|| (sp==3)) && (j==2))
{
c = Grid::Complex(1.0, 0.0);
F()(sp)(j) = c;
}
}
}
Grid::pokeSite(F, src, x);
}
}
}
}
*/
}
void calc_grid(ChromaAction action, Grid::QCD::LatticeGaugeField &Umu, Grid::QCD::LatticeFermion &src, Grid::QCD::LatticeFermion &res, int dag)
{
using namespace Grid;
using namespace Grid::QCD;
Grid::GridCartesian *UGrid = (Grid::GridCartesian *)Umu._grid;
Grid::GridRedBlackCartesian *UrbGrid = Grid::QCD::SpaceTimeGrid::makeFourDimRedBlackGrid(UGrid);
Grid::RealD _mass = mq;
if (action == Wilson)
{
Grid::QCD::WilsonFermionR Wf(Umu, *UGrid, *UrbGrid, _mass);
std::cout << Grid::GridLogMessage << " Calling Grid Wilson Fermion multiply " << std::endl;
if (dag)
Wf.Mdag(src, res);
else
Wf.M(src, res);
return;
}
if (action == WilsonClover)
{
Grid::RealD _csw = 1.0;
Grid::QCD::WilsonCloverFermionR Wf(Umu, *UGrid, *UrbGrid, _mass, _csw);
Wf.ImportGauge(Umu);
std::cout << Grid::GridLogMessage << " Calling Grid Wilson Clover Fermion multiply " << std::endl;
if (dag)
Wf.Mdag(src, res);
else
Wf.M(src, res);
return;
}
assert(0);
}
int main(int argc, char **argv)
{
/********************************************************
* Setup QDP
*********************************************************/
Chroma::initialize(&argc, &argv);
Chroma::WilsonTypeFermActs4DEnv::registerAll();
/********************************************************
* Setup Grid
*********************************************************/
Grid::Grid_init(&argc, &argv);
Grid::GridCartesian *UGrid = Grid::QCD::SpaceTimeGrid::makeFourDimGrid(Grid::GridDefaultLatt(),
Grid::GridDefaultSimd(Grid::QCD::Nd, Grid::vComplex::Nsimd()),
Grid::GridDefaultMpi());
std::vector<int> gd = UGrid->GlobalDimensions();
QDP::multi1d<int> nrow(QDP::Nd);
for (int mu = 0; mu < 4; mu++)
nrow[mu] = gd[mu];
QDP::Layout::setLattSize(nrow);
QDP::Layout::create();
GaugeField Ug(UGrid);
FermionField src(UGrid);
FermionField res_chroma(UGrid);
FermionField res_grid(UGrid);
FermionField only_wilson(UGrid);
FermionField difference(UGrid);
std::vector<ChromaAction> ActionList({Wilson, WilsonClover});
std::vector<std::string> ActionName({"Wilson", "WilsonClover"});
{
for (int i = 0; i < ActionList.size(); i++)
{
std::cout << "*****************************" << std::endl;
std::cout << "Action " << ActionName[i] << std::endl;
std::cout << "*****************************" << std::endl;
make_gauge(Ug, src); // fills the gauge field and the fermion field with random numbers
for (int dag = 0; dag < 2; dag++)
{
{
std::cout << "Dag = " << dag << std::endl;
calc_chroma(ActionList[i], Ug, src, res_chroma, dag);
// Remove the normalisation of Chroma Gauge links ????????
std::cout << "Norm of Chroma " << ActionName[i] << " multiply " << Grid::norm2(res_chroma) << std::endl;
calc_grid(ActionList[i], Ug, src, res_grid, dag);
std::cout << "Norm of gauge " << Grid::norm2(Ug) << std::endl;
std::cout << "Norm of Grid " << ActionName[i] << " multiply " << Grid::norm2(res_grid) << std::endl;
difference = res_chroma - res_grid;
std::cout << "Norm of difference " << Grid::norm2(difference) << std::endl;
}
}
std::cout << "Finished test " << std::endl;
Chroma::finalize();
}
}
}