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

Cleanup and prepare for pull request

Browse Source
This commit is contained in:
Guido Cossu 2017-10-24 13:21:17 +01:00
parent cbda4f66e0
commit ec8cd11c1f
4 changed files with 258 additions and 307 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

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

@ -39,29 +39,33 @@ namespace QCD
template <class Impl> template <class Impl>
RealD WilsonCloverFermion<Impl>::M(const FermionField &in, FermionField &out) RealD WilsonCloverFermion<Impl>::M(const FermionField &in, FermionField &out)
{ {
FermionField temp(out._grid);
// Wilson term // Wilson term
out.checkerboard = in.checkerboard; out.checkerboard = in.checkerboard;
this->Dhop(in, out, DaggerNo); this->Dhop(in, out, DaggerNo);
// Clover term // Clover term
// apply the sigma and Fmunu
FermionField temp(out._grid);
Mooee(in, temp); Mooee(in, temp);
out += temp; out += temp;
return axpy_norm(out, 4 + this->mass, in, out); return norm2(out);
} }
template <class Impl> template <class Impl>
RealD WilsonCloverFermion<Impl>::Mdag(const FermionField &in, FermionField &out) RealD WilsonCloverFermion<Impl>::Mdag(const FermionField &in, FermionField &out)
{ {
FermionField temp(out._grid);
// Wilson term // Wilson term
out.checkerboard = in.checkerboard; out.checkerboard = in.checkerboard;
this->Dhop(in, out, DaggerYes); this->Dhop(in, out, DaggerYes);
// Clover term // Clover term
// apply the sigma and Fmunu
FermionField temp(out._grid);
MooeeDag(in, temp); MooeeDag(in, temp);
out+=temp;
return axpy_norm(out, 4 + this->mass, in, out); out += temp;
return norm2(out);
} }
template <class Impl> template <class Impl>
@ -80,14 +84,14 @@ void WilsonCloverFermion<Impl>::ImportGauge(const GaugeField &_Umu)
WilsonLoops<Impl>::FieldStrength(Ez, _Umu, Tdir, Zdir); WilsonLoops<Impl>::FieldStrength(Ez, _Umu, Tdir, Zdir);
// Compute the Clover Operator acting on Colour and Spin // Compute the Clover Operator acting on Colour and Spin
CloverTerm = fillCloverYZ(Bx); CloverTerm = fillCloverYZ(Bx);
CloverTerm += fillCloverXZ(By); CloverTerm += fillCloverXZ(By);
CloverTerm += fillCloverXY(Bz); CloverTerm += fillCloverXY(Bz);
CloverTerm += fillCloverXT(Ex); CloverTerm += fillCloverXT(Ex);
CloverTerm += fillCloverYT(Ey); CloverTerm += fillCloverYT(Ey);
CloverTerm += fillCloverZT(Ez); CloverTerm += fillCloverZT(Ez);
CloverTerm *= (0.5) * csw; CloverTerm *= (0.5) * csw;
CloverTerm += (4.0 + this->mass);
int lvol = _Umu._grid->lSites(); int lvol = _Umu._grid->lSites();
int DimRep = Impl::Dimension; int DimRep = Impl::Dimension;
@ -98,20 +102,19 @@ void WilsonCloverFermion<Impl>::ImportGauge(const GaugeField &_Umu)
std::vector<int> lcoor; std::vector<int> lcoor;
typename SiteCloverType::scalar_object Qx = zero, Qxinv = zero; typename SiteCloverType::scalar_object Qx = zero, Qxinv = zero;
for (int site = 0; site < lvol; site++) for (int site = 0; site < lvol; site++)
{ {
grid->LocalIndexToLocalCoor(site, lcoor); grid->LocalIndexToLocalCoor(site, lcoor);
EigenCloverOp = Eigen::MatrixXcd::Zero(Ns * DimRep, Ns * DimRep); EigenCloverOp = Eigen::MatrixXcd::Zero(Ns * DimRep, Ns * DimRep);
peekLocalSite(Qx, CloverTerm, lcoor); peekLocalSite(Qx, CloverTerm, lcoor);
Qxinv = zero; Qxinv = zero;
//if (csw!=0){ //if (csw!=0){
for (int j = 0; j < Ns; j++) for (int j = 0; j < Ns; j++)
for (int k = 0; k < Ns; k++) for (int k = 0; k < Ns; k++)
for (int a = 0; a < DimRep; a++) for (int a = 0; a < DimRep; a++)
for (int b = 0; b < DimRep; b++) for (int b = 0; b < DimRep; b++)
EigenCloverOp(a + j * DimRep, b + k * DimRep) = Qx()(j, k)(a, b); EigenCloverOp(a + j * DimRep, b + k * DimRep) = Qx()(j, k)(a, b);
// if (site==0) std::cout << "site =" << site << "\n" << EigenCloverOp << std::endl; // if (site==0) std::cout << "site =" << site << "\n" << EigenCloverOp << std::endl;
EigenInvCloverOp = EigenCloverOp.inverse(); EigenInvCloverOp = EigenCloverOp.inverse();
//std::cout << EigenInvCloverOp << std::endl; //std::cout << EigenInvCloverOp << std::endl;
@ -120,35 +123,29 @@ void WilsonCloverFermion<Impl>::ImportGauge(const GaugeField &_Umu)
for (int a = 0; a < DimRep; a++) for (int a = 0; a < DimRep; a++)
for (int b = 0; b < DimRep; b++) for (int b = 0; b < DimRep; b++)
Qxinv()(j, k)(a, b) = EigenInvCloverOp(a + j * DimRep, b + k * DimRep); Qxinv()(j, k)(a, b) = EigenInvCloverOp(a + j * DimRep, b + k * DimRep);
// if (site==0) std::cout << "site =" << site << "\n" << EigenInvCloverOp << std::endl; // if (site==0) std::cout << "site =" << site << "\n" << EigenInvCloverOp << std::endl;
// } // }
pokeLocalSite(Qxinv, CloverTermInv, lcoor); pokeLocalSite(Qxinv, CloverTermInv, lcoor);
} }
// Separate the even and odd parts
// Separate the even and odd parts.
pickCheckerboard(Even, CloverTermEven, CloverTerm); pickCheckerboard(Even, CloverTermEven, CloverTerm);
pickCheckerboard( Odd, CloverTermOdd, CloverTerm); pickCheckerboard(Odd, CloverTermOdd, CloverTerm);
pickCheckerboard(Even, CloverTermDagEven, adj(CloverTerm)); pickCheckerboard(Even, CloverTermDagEven, adj(CloverTerm));
pickCheckerboard( Odd, CloverTermDagOdd, adj(CloverTerm)); pickCheckerboard(Odd, CloverTermDagOdd, adj(CloverTerm));
pickCheckerboard(Even, CloverTermInvEven, CloverTermInv); pickCheckerboard(Even, CloverTermInvEven, CloverTermInv);
pickCheckerboard( Odd, CloverTermInvOdd, CloverTermInv); pickCheckerboard(Odd, CloverTermInvOdd, CloverTermInv);
pickCheckerboard(Even, CloverTermInvDagEven, adj(CloverTermInv)); pickCheckerboard(Even, CloverTermInvDagEven, adj(CloverTermInv));
pickCheckerboard( Odd, CloverTermInvDagOdd, adj(CloverTermInv)); pickCheckerboard(Odd, CloverTermInvDagOdd, adj(CloverTermInv));
} }
template <class Impl> template <class Impl>
void WilsonCloverFermion<Impl>::Mooee(const FermionField &in, FermionField &out) void WilsonCloverFermion<Impl>::Mooee(const FermionField &in, FermionField &out)
{ {
conformable(in,out); conformable(in, out);
this->MooeeInternal(in, out, DaggerNo, InverseNo); this->MooeeInternal(in, out, DaggerNo, InverseNo);
} }
@ -177,85 +174,50 @@ void WilsonCloverFermion<Impl>::MooeeInternal(const FermionField &in, FermionFie
CloverFieldType *Clover; CloverFieldType *Clover;
assert(in.checkerboard == Odd || in.checkerboard == Even); assert(in.checkerboard == Odd || in.checkerboard == Even);
if (dag)
{
if (in._grid->_isCheckerBoarded)
if (dag){ {
if (in._grid->_isCheckerBoarded){ if (in.checkerboard == Odd)
if (in.checkerboard == Odd){ {
// std::cout << "Calling clover term adj Odd" << std::endl; Clover = (inv) ? &CloverTermInvDagOdd : &CloverTermDagOdd;
Clover = (inv) ? &CloverTermInvDagOdd : &CloverTermDagOdd; }
else
/* test {
int DimRep = Impl::Dimension; Clover = (inv) ? &CloverTermInvDagEven : &CloverTermDagEven;
Eigen::MatrixXcd A = Eigen::MatrixXcd::Zero(Ns * DimRep, Ns * DimRep); }
std::vector<int> lcoor; out = *Clover * in;
typename SiteCloverType::scalar_object Qx2 = zero;
GridBase *grid = in._grid;
int site = 0 ;
grid->LocalIndexToLocalCoor(site, lcoor);
peekLocalSite(Qx2, *Clover, lcoor);
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++)
A(a + j * DimRep, b + k * DimRep) = Qx2()(j, k)(a, b);
std::cout << "adj Odd =" << site << "\n" << A << std::endl;
end test */
} else {
// std::cout << "Calling clover term adj Even" << std::endl;
Clover = (inv) ? &CloverTermInvDagEven : &CloverTermDagEven;
/* test
int DimRep = Impl::Dimension;
Eigen::MatrixXcd A = Eigen::MatrixXcd::Zero(Ns * DimRep, Ns * DimRep);
std::vector<int> lcoor;
typename SiteCloverType::scalar_object Qx2 = zero;
GridBase *grid = in._grid;
int site = 0 ;
grid->LocalIndexToLocalCoor(site, lcoor);
peekLocalSite(Qx2, *Clover, lcoor);
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++)
A(a + j * DimRep, b + k * DimRep) = Qx2()(j, k)(a, b);
std::cout << "adj Odd =" << site << "\n" << A << std::endl;
end test */
} }
// std::cout << GridLogMessage << "*Clover.checkerboard " << (*Clover).checkerboard << std::endl; else
out = *Clover * in; {
} else { Clover = (inv) ? &CloverTermInv : &CloverTerm;
Clover = (inv) ? &CloverTermInv : &CloverTerm; out = adj(*Clover) * in;
//out = adj(*Clover) * in; }
out = adj(CloverTerm) * in; }
} else
{
if (in._grid->_isCheckerBoarded)
{
} else { if (in.checkerboard == Odd)
if (in._grid->_isCheckerBoarded){ {
// std::cout << "Calling clover term Odd" << std::endl;
if (in.checkerboard == Odd){ Clover = (inv) ? &CloverTermInvOdd : &CloverTermOdd;
// std::cout << "Calling clover term Odd" << std::endl; }
Clover = (inv) ? &CloverTermInvOdd : &CloverTermOdd; else
} else { {
// std::cout << "Calling clover term Even" << std::endl; // std::cout << "Calling clover term Even" << std::endl;
Clover = (inv) ? &CloverTermInvEven : &CloverTermEven; Clover = (inv) ? &CloverTermInvEven : &CloverTermEven;
}
out = *Clover * in;
// std::cout << GridLogMessage << "*Clover.checkerboard " << (*Clover).checkerboard << std::endl;
}
else
{
Clover = (inv) ? &CloverTermInv : &CloverTerm;
out = *Clover * in;
} }
out = *Clover * in;
// std::cout << GridLogMessage << "*Clover.checkerboard " << (*Clover).checkerboard << std::endl;
} else {
Clover = (inv) ? &CloverTermInv : &CloverTerm;
out = *Clover * in;
} }
}
} // MooeeInternal } // MooeeInternal
@ -264,7 +226,6 @@ template <class Impl>
void WilsonCloverFermion<Impl>::MDeriv(GaugeField &mat, const FermionField &U, const FermionField &V, int dag) void WilsonCloverFermion<Impl>::MDeriv(GaugeField &mat, const FermionField &U, const FermionField &V, int dag)
{ {
GaugeField tmp(mat._grid); GaugeField tmp(mat._grid);
conformable(U._grid, V._grid); conformable(U._grid, V._grid);
@ -283,29 +244,36 @@ template <class Impl>
void WilsonCloverFermion<Impl>::MooDeriv(GaugeField &mat, const FermionField &X, const FermionField &Y, int dag) void WilsonCloverFermion<Impl>::MooDeriv(GaugeField &mat, const FermionField &X, const FermionField &Y, int dag)
{ {
GridBase *grid = mat._grid; GridBase *grid = mat._grid;
//GaugeLinkField Lambdaodd(grid), Lambdaeven(grid), tmp(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 //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) //Lambdaeven = zero; //Teven*dag(Xeven)+Xeven*dag(Yeven) + 2*(Dee^-1)
GaugeLinkField Lambda(grid), tmp(grid); GaugeLinkField Lambda(grid), tmp(grid);
Lambda=zero; Lambda = zero;
conformable(mat._grid, X._grid); conformable(mat._grid, X._grid);
conformable(Y._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> 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> C1m(Nd, grid), C2m(Nd, grid), C3m(Nd, grid), C4m(Nd, grid);
std::vector<GaugeLinkField> U(Nd, mat._grid); std::vector<GaugeLinkField> U(Nd, mat._grid);
for (int mu = 0; mu < Nd; mu++) { 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; U[mu] = PeekIndex<LorentzIndex>(mat, mu);
C1m[mu]=zero; C2m[mu]=zero; C3m[mu]=zero; C4m[mu]=zero; 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 PARALLEL_FOR_LOOP
for (int i = 0; i < CloverTerm._grid->oSites(); i++) for (int i = 0; i < CloverTerm._grid->oSites(); i++)
{ {
@ -316,48 +284,47 @@ for (int mu = 0; mu < Nd; mu++) {
} }
*/ */
for (int i=0;i<4;i++){ //spin for (int i = 0; i < 4; i++)
for(int j=0;j<4;j++){ //spin { //spin
for (int j = 0; j < 4; j++)
{ //spin
for (int mu=0;mu<4;mu++){ //color for (int mu = 0; mu < 4; mu++)
for (int nu=0;nu<4;nu++){ //color { //color
for (int nu = 0; nu < 4; nu++)
{ //color
// insertion in upper staple // insertion in upper staple
tmp = Lambda * U[nu]; tmp = Lambda * U[nu];
C1p[mu]+=Impl::ShiftStaple(Impl::CovShiftForward(tmp, nu, Impl::CovShiftBackward(U[mu], mu, Impl::CovShiftIdentityBackward(U[nu], nu))), mu); C1p[mu] += Impl::ShiftStaple(Impl::CovShiftForward(tmp, nu, Impl::CovShiftBackward(U[mu], mu, Impl::CovShiftIdentityBackward(U[nu], nu))), mu);
tmp = Lambda * U[mu]; tmp = Lambda * U[mu];
C2p[mu]+= Impl::ShiftStaple(Impl::CovShiftForward(U[nu], nu, Impl::CovShiftBackward(tmp, mu, Impl::CovShiftIdentityBackward(U[nu], nu))), 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]; 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); C3p[mu] += Impl::ShiftStaple(Impl::CovShiftForward(U[nu], nu, Impl::CovShiftBackward(U[mu], mu, Impl::CovShiftIdentityBackward(tmp, nu))), mu);
tmp = Lambda; tmp = Lambda;
C4p[mu]+= Impl::ShiftStaple(Impl::CovShiftForward(U[nu], nu, Impl::CovShiftBackward(U[mu], mu, Impl::CovShiftIdentityBackward(U[nu], nu))),mu) * tmp; 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 // insertion in lower staple
tmp = Lambda * U[nu]; tmp = Lambda * U[nu];
C1m[mu]+= Impl::ShiftStaple(Impl::CovShiftBackward(tmp, nu, Impl::CovShiftBackward(U[mu], mu, U[nu])), mu); C1m[mu] += Impl::ShiftStaple(Impl::CovShiftBackward(tmp, nu, Impl::CovShiftBackward(U[mu], mu, U[nu])), mu);
tmp = Lambda * U[mu]; tmp = Lambda * U[mu];
C2m[mu]+= Impl::ShiftStaple(Impl::CovShiftBackward(U[nu], nu, Impl::CovShiftBackward(tmp, mu, U[nu])), mu); C2m[mu] += Impl::ShiftStaple(Impl::CovShiftBackward(U[nu], nu, Impl::CovShiftBackward(tmp, mu, U[nu])), mu);
tmp = Lambda * U[nu]; tmp = Lambda * U[nu];
C3m[mu]+= Impl::ShiftStaple(Impl::CovShiftBackward(U[nu], nu, Impl::CovShiftBackward(U[mu], mu, tmp)), mu); C3m[mu] += Impl::ShiftStaple(Impl::CovShiftBackward(U[nu], nu, Impl::CovShiftBackward(U[mu], mu, tmp)), mu);
tmp = Lambda; tmp = Lambda;
C4m[mu]+= Impl::ShiftStaple(Impl::CovShiftBackward(U[nu], nu, Impl::CovShiftBackward(U[mu], mu, U[nu])), mu)* tmp; 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
//Still implementing. Have to be tested, and understood how to project EO
} }
// Derivative parts // Derivative parts

70
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 See the full license in the file "LICENSE" in the top level distribution directory
*************************************************************************************/ *************************************************************************************/
/* END LEGAL */ /* END LEGAL */
#ifndef GRID_QCD_WILSON_CLOVER_FERMION_H #ifndef GRID_QCD_WILSON_CLOVER_FERMION_H
#define GRID_QCD_WILSON_CLOVER_FERMION_H #define GRID_QCD_WILSON_CLOVER_FERMION_H
@ -42,9 +43,11 @@ class WilsonCloverFermion : public WilsonFermion<Impl>
public: public:
// Types definitions // Types definitions
INHERIT_IMPL_TYPES(Impl); INHERIT_IMPL_TYPES(Impl);
template <typename vtype> using iImplClover = iScalar<iMatrix<iMatrix<vtype, Impl::Dimension>, Ns> >; template <typename vtype>
typedef iImplClover<Simd> SiteCloverType; using iImplClover = iScalar<iMatrix<iMatrix<vtype, Impl::Dimension>, Ns>>;
typedef Lattice<SiteCloverType> CloverFieldType; typedef iImplClover<Simd> SiteCloverType;
typedef Lattice<SiteCloverType> CloverFieldType;
public: public:
typedef WilsonFermion<Impl> WilsonBase; typedef WilsonFermion<Impl> WilsonBase;
@ -58,19 +61,21 @@ public:
Fgrid, Fgrid,
Hgrid, Hgrid,
_mass, p), _mass, p),
CloverTerm(&Fgrid), CloverTerm(&Fgrid),
CloverTermInv(&Fgrid), CloverTermInv(&Fgrid),
CloverTermEven(&Hgrid), CloverTermEven(&Hgrid),
CloverTermOdd(&Hgrid), CloverTermOdd(&Hgrid),
CloverTermInvEven(&Hgrid), CloverTermInvEven(&Hgrid),
CloverTermInvOdd(&Hgrid), CloverTermInvOdd(&Hgrid),
CloverTermDagEven(&Hgrid), //test CloverTermDagEven(&Hgrid),
CloverTermDagOdd(&Hgrid), //test CloverTermDagOdd(&Hgrid),
CloverTermInvDagEven(&Hgrid), //test CloverTermInvDagEven(&Hgrid),
CloverTermInvDagOdd(&Hgrid) //test CloverTermInvDagOdd(&Hgrid)
{ {
csw = _csw; csw = _csw;
assert(Nd == 4); // require 4 dimensions 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); virtual RealD M(const FermionField &in, FermionField &out);
@ -91,12 +96,12 @@ public:
private: private:
// here fixing the 4 dimensions, make it more general? // here fixing the 4 dimensions, make it more general?
RealD csw; // Clover coefficient RealD csw; // Clover coefficient
CloverFieldType CloverTerm=zero, CloverTermInv=zero; // Clover term CloverFieldType CloverTerm, CloverTermInv; // Clover term
CloverFieldType CloverTermEven=zero, CloverTermOdd=zero; // Clover term EO CloverFieldType CloverTermEven, CloverTermOdd; // Clover term EO
CloverFieldType CloverTermInvEven=zero, CloverTermInvOdd=zero; // Clover term Inv EO CloverFieldType CloverTermInvEven, CloverTermInvOdd; // Clover term Inv EO
CloverFieldType CloverTermDagEven=zero, CloverTermDagOdd=zero; // Clover term Dag EO CloverFieldType CloverTermDagEven, CloverTermDagOdd; // Clover term Dag EO
CloverFieldType CloverTermInvDagEven=zero, CloverTermInvDagOdd=zero; // Clover term Inv Dag EO CloverFieldType CloverTermInvDagEven, CloverTermInvDagOdd; // Clover term Inv Dag EO
// eventually these two can be compressed into 6x6 blocks instead of the 12x12 // eventually these two can be compressed into 6x6 blocks instead of the 12x12
// using the DeGrand-Rossi basis for the gamma matrices // using the DeGrand-Rossi basis for the gamma matrices
@ -114,8 +119,8 @@ private:
T._odata[i]()(3, 2) = timesMinusI(F._odata[i]()()); T._odata[i]()(3, 2) = timesMinusI(F._odata[i]()());
} }
return T; return T;
} }
CloverFieldType fillCloverXZ(const GaugeLinkField &F) CloverFieldType fillCloverXZ(const GaugeLinkField &F)
{ {
@ -130,8 +135,8 @@ private:
T._odata[i]()(3, 2) = F._odata[i]()(); T._odata[i]()(3, 2) = F._odata[i]()();
} }
return T; return T;
} }
CloverFieldType fillCloverXY(const GaugeLinkField &F) CloverFieldType fillCloverXY(const GaugeLinkField &F)
{ {
@ -146,8 +151,8 @@ private:
T._odata[i]()(3, 3) = timesI(F._odata[i]()()); T._odata[i]()(3, 3) = timesI(F._odata[i]()());
} }
return T; return T;
} }
CloverFieldType fillCloverXT(const GaugeLinkField &F) CloverFieldType fillCloverXT(const GaugeLinkField &F)
{ {
@ -162,8 +167,8 @@ private:
T._odata[i]()(3, 2) = timesMinusI(F._odata[i]()()); T._odata[i]()(3, 2) = timesMinusI(F._odata[i]()());
} }
return T; return T;
} }
CloverFieldType fillCloverYT(const GaugeLinkField &F) CloverFieldType fillCloverYT(const GaugeLinkField &F)
{ {
@ -178,8 +183,8 @@ private:
T._odata[i]()(3, 2) = -(F._odata[i]()()); T._odata[i]()(3, 2) = -(F._odata[i]()());
} }
return T; return T;
} }
CloverFieldType fillCloverZT(const GaugeLinkField &F) CloverFieldType fillCloverZT(const GaugeLinkField &F)
{ {
@ -194,11 +199,10 @@ private:
T._odata[i]()(3, 3) = timesI(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

10
tests/core/Test_wilson_clover.cc
View File

@ -237,7 +237,7 @@ int main (int argc, char ** argv)
setCheckerboard(src,src_o); setCheckerboard(src,src_o);
//Gauge Transformation ////////////////////// Gauge Transformation
std::vector<int> seeds2({5,6,7,8}); std::vector<int> seeds2({5,6,7,8});
GridParallelRNG pRNG2(&Grid); pRNG2.SeedFixedIntegers(seeds2); GridParallelRNG pRNG2(&Grid); pRNG2.SeedFixedIntegers(seeds2);
LatticeColourMatrix Omega(&Grid); LatticeColourMatrix Omega(&Grid);
@ -251,7 +251,7 @@ int main (int argc, char ** argv)
U_prime_mu=Omega*U[mu]*adj(ShiftedOmega); U_prime_mu=Omega*U[mu]*adj(ShiftedOmega);
pokeLorentz(U_prime,U_prime_mu,mu); pokeLorentz(U_prime,U_prime_mu,mu);
} }
/////////////////
WilsonCloverFermionR Dwc_prime(U_prime,Grid,RBGrid,mass,csw,params); WilsonCloverFermionR Dwc_prime(U_prime,Grid,RBGrid,mass,csw,params);
Dwc_prime.ImportGauge(U_prime); Dwc_prime.ImportGauge(U_prime);
@ -298,7 +298,7 @@ int main (int argc, char ** argv)
std::cout<<GridLogMessage<<"=========================================================="<<std::endl; std::cout<<GridLogMessage<<"=========================================================="<<std::endl;
chi=zero; phi=zero; err=zero; chi=zero; phi=zero; err=zero;
WilsonCloverFermionR Dwc_csw0(Umu,Grid,RBGrid,mass,0.0,params); // <-- csw=0 WilsonCloverFermionR Dwc_csw0(Umu,Grid,RBGrid,mass,0.0,params); // <-- Notice: csw=0
Dwc_csw0.ImportGauge(Umu); Dwc_csw0.ImportGauge(Umu);
pickCheckerboard(Even,phi_e,phi); pickCheckerboard(Even,phi_e,phi);
@ -318,9 +318,7 @@ int main (int argc, char ** argv)
setCheckerboard(src,src_e); setCheckerboard(src,src_e);
setCheckerboard(src,src_o); setCheckerboard(src,src_o);
FermionField::scalar_type scal(4.0 + mass); err = chi - phi;
err = chi - (phi + scal*src) ; // subtraction of the mass term (not present in Mooee Clover!)
std::cout<<GridLogMessage << "norm diff "<< norm2(err)<< std::endl; std::cout<<GridLogMessage << "norm diff "<< norm2(err)<< std::endl;
Grid_finalize(); Grid_finalize();

172
tests/qdpxx/Test_qdpxx_wilson.cc
View File

@ -27,28 +27,23 @@
/* END LEGAL */ /* END LEGAL */
#include <Grid/Grid.h> #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 // Mass
double mq = 0.0; double mq = 0.01;
// Define Wilson Types // Define Wilson Types
typedef Grid::QCD::WilsonImplR::FermionField FermionField; typedef Grid::QCD::WilsonImplR::FermionField FermionField;
typedef Grid::QCD::LatticeGaugeField GaugeField; typedef Grid::QCD::LatticeGaugeField GaugeField;
#include <chroma.h>
#include <actions/ferm/invert/syssolver_linop_cg_array.h>
#include <actions/ferm/invert/syssolver_linop_aggregate.h>
enum ChromaAction enum ChromaAction
{ {
Wilson, // Wilson Wilson, // Wilson
WilsonClover // CloverFermions WilsonClover // CloverFermions
}; };
void make_gauge(GaugeField &lat, FermionField &src);
void calc_grid(ChromaAction CA, GaugeField &lat, FermionField &src, FermionField &res, int dag);
void calc_chroma(ChromaAction CA, GaugeField &lat, FermionField &src, FermionField &res, int dag);
namespace Chroma namespace Chroma
{ {
@ -286,91 +281,6 @@ public:
}; };
} // namespace Chroma } // namespace Chroma
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;
// Isolate Clover term
/*
calc_grid(Wilson, Ug, src, only_wilson, dag); // Wilson term
res_grid -= only_wilson;
res_chroma -= only_wilson;
std::cout << "Chroma:" << res_chroma << std::endl;
std::cout << "Grid :" << res_grid << std::endl;
difference = (res_grid-res_chroma);
std::cout << "Difference :" << difference << std::endl;
*/
}
}
std::cout << "Finished test " << std::endl;
Chroma::finalize();
}
}
}
void calc_chroma(ChromaAction action, GaugeField &lat, FermionField &src, FermionField &res, int dag) void calc_chroma(ChromaAction action, GaugeField &lat, FermionField &src, FermionField &res, int dag)
{ {
QDP::multi1d<QDP::LatticeColorMatrix> u(4); QDP::multi1d<QDP::LatticeColorMatrix> u(4);
@ -467,7 +377,6 @@ void make_gauge(GaugeField &Umu, FermionField &src)
} }
} }
*/ */
} }
void calc_grid(ChromaAction action, Grid::QCD::LatticeGaugeField &Umu, Grid::QCD::LatticeFermion &src, Grid::QCD::LatticeFermion &res, int dag) void calc_grid(ChromaAction action, Grid::QCD::LatticeGaugeField &Umu, Grid::QCD::LatticeFermion &src, Grid::QCD::LatticeFermion &res, int dag)
@ -512,3 +421,76 @@ void calc_grid(ChromaAction action, Grid::QCD::LatticeGaugeField &Umu, Grid::QCD
assert(0); 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();
}
}
}