diff --git a/.gitignore b/.gitignore index d743ee06..6b92b1a4 100644 --- a/.gitignore +++ b/.gitignore @@ -123,3 +123,5 @@ make-bin-BUCK.sh lib/qcd/spin/gamma-gen/*.h lib/qcd/spin/gamma-gen/*.cc +.vscode/settings.json +settings.json diff --git a/.vscode/settings.json b/.vscode/settings.json index dd8f0473..3e49029b 100644 --- a/.vscode/settings.json +++ b/.vscode/settings.json @@ -36,6 +36,15 @@ "tuple": "cpp", "type_traits": "cpp", "typeinfo": "cpp", - "utility": "cpp" + "utility": "cpp", + "iostream": "cpp", + "strstream": "cpp", + "complex": "cpp", + "fstream": "cpp", + "iomanip": "cpp", + "istream": "cpp", + "ostream": "cpp", + "sstream": "cpp", + "streambuf": "cpp" } } \ No newline at end of file diff --git a/lib/cartesian/Cartesian_base.h b/lib/cartesian/Cartesian_base.h index 0db6ce0d..e435bbba 100644 --- a/lib/cartesian/Cartesian_base.h +++ b/lib/cartesian/Cartesian_base.h @@ -69,6 +69,8 @@ public: std::vector _lstart; // local start of array in gcoors _processor_coor[d]*_ldimensions[d] std::vector _lend ; // local end of array in gcoors _processor_coor[d]*_ldimensions[d]+_ldimensions_[d]-1 + bool _isCheckerBoarded; + public: //////////////////////////////////////////////////////////////// diff --git a/lib/cartesian/Cartesian_full.h b/lib/cartesian/Cartesian_full.h index b0e47fa4..3be3e8cd 100644 --- a/lib/cartesian/Cartesian_full.h +++ b/lib/cartesian/Cartesian_full.h @@ -69,6 +69,7 @@ public: /////////////////////// // Grid information /////////////////////// + _isCheckerBoarded = false; _ndimension = dimensions.size(); _fdimensions.resize(_ndimension); @@ -76,8 +77,8 @@ public: _ldimensions.resize(_ndimension); _rdimensions.resize(_ndimension); _simd_layout.resize(_ndimension); - _lstart.resize(_ndimension); - _lend.resize(_ndimension); + _lstart.resize(_ndimension); + _lend.resize(_ndimension); _ostride.resize(_ndimension); _istride.resize(_ndimension); diff --git a/lib/cartesian/Cartesian_red_black.h b/lib/cartesian/Cartesian_red_black.h index 3037de00..a440538a 100644 --- a/lib/cartesian/Cartesian_red_black.h +++ b/lib/cartesian/Cartesian_red_black.h @@ -139,6 +139,7 @@ public: /////////////////////// // Grid information /////////////////////// + _isCheckerBoarded = true; _checker_dim = checker_dim; assert(checker_dim_mask[checker_dim]==1); _ndimension = dimensions.size(); diff --git a/lib/qcd/action/fermion/WilsonCloverFermion.cc b/lib/qcd/action/fermion/WilsonCloverFermion.cc index e1900830..f8b62ba4 100644 --- a/lib/qcd/action/fermion/WilsonCloverFermion.cc +++ b/lib/qcd/action/fermion/WilsonCloverFermion.cc @@ -35,27 +35,6 @@ namespace Grid namespace QCD { -//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) @@ -89,10 +68,10 @@ RealD WilsonCloverFermion::Mdag(const FermionField &in, FermionField &out) template void WilsonCloverFermion::ImportGauge(const GaugeField &_Umu) { - this->ImportGauge(_Umu); + WilsonFermion::ImportGauge(_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::FieldStrength(Bx, _Umu, Ydir, Zdir); WilsonLoops::FieldStrength(By, _Umu, Zdir, Xdir); @@ -102,12 +81,12 @@ void WilsonCloverFermion::ImportGauge(const GaugeField &_Umu) WilsonLoops::FieldStrength(Ez, _Umu, Tdir, Zdir); // 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 = fillCloverYZ(Bx); + CloverTerm += fillCloverXZ(By); + CloverTerm += fillCloverXY(Bz); + CloverTerm += fillCloverXT(Ex); + CloverTerm += fillCloverYT(Ey); + CloverTerm += fillCloverZT(Ez) ; CloverTerm *= csw; int lvol = _Umu._grid->lSites(); @@ -130,8 +109,11 @@ void WilsonCloverFermion::ImportGauge(const GaugeField &_Umu) 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; + EigenInvCloverOp = EigenCloverOp.inverse(); + //std::cout << EigenInvCloverOp << std::endl; for (int j = 0; j < Ns; j++) for (int k = 0; k < Ns; k++) for (int a = 0; a < DimRep; a++) @@ -139,17 +121,21 @@ void WilsonCloverFermion::ImportGauge(const GaugeField &_Umu) Qxinv()(j, k)(a, b) = EigenInvCloverOp(a + j * DimRep, b + k * DimRep); pokeLocalSite(Qxinv, CloverTermInv, lcoor); - // Separate the even and odd parts. - pickCheckerboard(Even, CloverTermEven, CloverTerm); - pickCheckerboard( Odd, CloverTermOdd, CloverTerm); - pickCheckerboard(Even, CloverTermInvEven, CloverTermInv); - pickCheckerboard( Odd, CloverTermInvOdd, CloverTermInv); } + + // Separate the even and odd parts. + pickCheckerboard(Even, CloverTermEven, CloverTerm); + pickCheckerboard( Odd, CloverTermOdd, CloverTerm); + + pickCheckerboard(Even, CloverTermInvEven, CloverTermInv); + pickCheckerboard( Odd, CloverTermInvOdd, CloverTermInv); + } template void WilsonCloverFermion::Mooee(const FermionField &in, FermionField &out) { + conformable(in,out); this->MooeeInternal(in, out, DaggerNo, InverseNo); } @@ -176,15 +162,27 @@ void WilsonCloverFermion::MooeeInternal(const FermionField &in, FermionFie { out.checkerboard = in.checkerboard; CloverFieldType *Clover; - if (in.checkerboard == Odd){ - std::cout << "Calling clover term Odd" << std::endl; - Clover = (inv) ? &CloverTermInvOdd : &CloverTermOdd; + assert(in.checkerboard == Odd || in.checkerboard == Even); + + 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; + } } - if (in.checkerboard == Even){ - std::cout << "Calling clover term Even" << std::endl; - Clover = (inv) ? &CloverTermInvEven : &CloverTermEven; + else + { + Clover = (inv) ? &CloverTermInv : &CloverTerm; } + std::cout << GridLogMessage << "*Clover.checkerboard " << (*Clover).checkerboard << std::endl; if (dag){ out = adj(*Clover) * in;} else { out = *Clover * in;} } // MooeeInternal diff --git a/lib/qcd/action/fermion/WilsonCloverFermion.h b/lib/qcd/action/fermion/WilsonCloverFermion.h index c9e7be39..fd9d1f60 100644 --- a/lib/qcd/action/fermion/WilsonCloverFermion.h +++ b/lib/qcd/action/fermion/WilsonCloverFermion.h @@ -94,17 +94,102 @@ private: // eventually these two can be compressed into 6x6 blocks instead of the 12x12 // using the DeGrand-Rossi basis for the gamma matrices - CloverFieldType fillClover(const GaugeLinkField& F){ + CloverFieldType fillCloverYZ(const GaugeLinkField &F) + { CloverFieldType T(F._grid); + T = zero; 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]()(); + 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]()()); } - return T; - } + return T; +} + + CloverFieldType fillCloverXZ(const GaugeLinkField &F) + { + CloverFieldType T(F._grid); + T = zero; + 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]()(); + } + + return T; +} + + CloverFieldType fillCloverXY(const GaugeLinkField &F) + { + CloverFieldType T(F._grid); + T = zero; + 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) = timesMinusI(F._odata[i]()()); + T._odata[i]()(3, 3) = timesI(F._odata[i]()()); + } + + return T; +} + + CloverFieldType fillCloverXT(const GaugeLinkField &F) + { + CloverFieldType T(F._grid); + T = 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) = timesI(F._odata[i]()()); + T._odata[i]()(3, 2) = timesI(F._odata[i]()()); + } + + return T; +} + + CloverFieldType fillCloverYT(const GaugeLinkField &F) + { + CloverFieldType T(F._grid); + T = zero; + 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]()()); + } + + return T; +} + + CloverFieldType fillCloverZT(const GaugeLinkField &F) + { + CloverFieldType T(F._grid); + T = zero; + 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]()()); + } + + return T; +} + }; } } diff --git a/tests/core/Test_wilson_clover.cc b/tests/core/Test_wilson_clover.cc index 3df69e3b..1b208e2d 100644 --- a/tests/core/Test_wilson_clover.cc +++ b/tests/core/Test_wilson_clover.cc @@ -91,7 +91,7 @@ int main (int argc, char ** argv) } WilsonCloverFermionR Dwc(Umu,Grid,RBGrid,mass,csw,params); - + Dwc.ImportGauge(Umu); std::cout<