From ec8cd11c1f7fce1c3deee79977745ba4f6c9776c Mon Sep 17 00:00:00 2001 From: Guido Cossu Date: Tue, 24 Oct 2017 13:21:17 +0100 Subject: [PATCH] Cleanup and prepare for pull request --- lib/qcd/action/fermion/WilsonCloverFermion.cc | 279 ++++++++---------- lib/qcd/action/fermion/WilsonCloverFermion.h | 104 +++---- tests/core/Test_wilson_clover.cc | 10 +- tests/qdpxx/Test_qdpxx_wilson.cc | 172 +++++------ 4 files changed, 258 insertions(+), 307 deletions(-) diff --git a/lib/qcd/action/fermion/WilsonCloverFermion.cc b/lib/qcd/action/fermion/WilsonCloverFermion.cc index 5e7e0034..73e2bf69 100644 --- a/lib/qcd/action/fermion/WilsonCloverFermion.cc +++ b/lib/qcd/action/fermion/WilsonCloverFermion.cc @@ -39,29 +39,33 @@ namespace QCD template RealD WilsonCloverFermion::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); + out += temp; - return axpy_norm(out, 4 + this->mass, in, out); + return norm2(out); } template RealD WilsonCloverFermion::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 @@ -80,14 +84,14 @@ void WilsonCloverFermion::ImportGauge(const GaugeField &_Umu) WilsonLoops::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 *= (0.5) * csw; - + CloverTerm *= (0.5) * csw; + CloverTerm += (4.0 + this->mass); int lvol = _Umu._grid->lSites(); int DimRep = Impl::Dimension; @@ -98,21 +102,20 @@ void WilsonCloverFermion::ImportGauge(const GaugeField &_Umu) 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; -//if (csw!=0){ + //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); - // 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(); //std::cout << EigenInvCloverOp << std::endl; for (int j = 0; j < Ns; j++) @@ -120,35 +123,29 @@ void WilsonCloverFermion::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; -// } + // 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(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)); - + pickCheckerboard(Odd, CloverTermInvDagOdd, adj(CloverTermInv)); } template void WilsonCloverFermion::Mooee(const FermionField &in, FermionField &out) { - conformable(in,out); + conformable(in, out); this->MooeeInternal(in, out, DaggerNo, InverseNo); } @@ -177,85 +174,50 @@ void WilsonCloverFermion::MooeeInternal(const FermionField &in, FermionFie CloverFieldType *Clover; assert(in.checkerboard == Odd || in.checkerboard == Even); - - - - if (dag){ - if (in._grid->_isCheckerBoarded){ - if (in.checkerboard == Odd){ -// std::cout << "Calling clover term adj Odd" << std::endl; - Clover = (inv) ? &CloverTermInvDagOdd : &CloverTermDagOdd; - -/* test - int DimRep = Impl::Dimension; - Eigen::MatrixXcd A = Eigen::MatrixXcd::Zero(Ns * DimRep, Ns * DimRep); - std::vector 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 */ - - - - } 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 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 */ - - + if (dag) + { + if (in._grid->_isCheckerBoarded) + { + if (in.checkerboard == Odd) + { + Clover = (inv) ? &CloverTermInvDagOdd : &CloverTermDagOdd; + } + else + { + Clover = (inv) ? &CloverTermInvDagEven : &CloverTermDagEven; + } + out = *Clover * in; + } + else + { + Clover = (inv) ? &CloverTermInv : &CloverTerm; + out = adj(*Clover) * in; } - // std::cout << GridLogMessage << "*Clover.checkerboard " << (*Clover).checkerboard << std::endl; - out = *Clover * in; - } else { - Clover = (inv) ? &CloverTermInv : &CloverTerm; - //out = adj(*Clover) * in; - out = adj(CloverTerm) * in; } + else + { + if (in._grid->_isCheckerBoarded) + { - - - - } else { - 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; + 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; + } } - } } // MooeeInternal @@ -264,7 +226,6 @@ template void WilsonCloverFermion::MDeriv(GaugeField &mat, const FermionField &U, const FermionField &V, int dag) { - GaugeField tmp(mat._grid); conformable(U._grid, V._grid); @@ -282,30 +243,37 @@ void WilsonCloverFermion::MDeriv(GaugeField &mat, const FermionField &U, c template void WilsonCloverFermion::MooDeriv(GaugeField &mat, const FermionField &X, const FermionField &Y, int dag) { - -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) + GridBase *grid = mat._grid; -GaugeLinkField Lambda(grid), tmp(grid); -Lambda=zero; + //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) -conformable(mat._grid, X._grid); -conformable(Y._grid, X._grid); + GaugeLinkField Lambda(grid), tmp(grid); + Lambda = zero; -std::vector C1p(Nd,grid), C2p(Nd,grid), C3p(Nd,grid), C4p(Nd,grid); -std::vector C1m(Nd,grid), C2m(Nd,grid), C3m(Nd,grid), C4m(Nd,grid); -std::vector U(Nd, mat._grid); + conformable(mat._grid, X._grid); + conformable(Y._grid, X._grid); -for (int mu = 0; mu < Nd; mu++) { - U[mu] = PeekIndex(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; -} + std::vector C1p(Nd, grid), C2p(Nd, grid), C3p(Nd, grid), C4p(Nd, grid); + std::vector C1m(Nd, grid), C2m(Nd, grid), C3m(Nd, grid), C4m(Nd, grid); + std::vector U(Nd, mat._grid); -/* + for (int mu = 0; mu < Nd; mu++) + { + U[mu] = PeekIndex(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++) { @@ -314,50 +282,49 @@ for (int mu = 0; mu < Nd; mu++) { 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 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 + 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); + // 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 = 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; + 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); -// 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; + C4p[mu] += Impl::ShiftStaple(Impl::CovShiftForward(U[nu], nu, Impl::CovShiftBackward(U[mu], mu, Impl::CovShiftIdentityBackward(U[nu], nu))), mu) * tmp; - tmp = Lambda * U[mu]; - C2m[mu]+= Impl::ShiftStaple(Impl::CovShiftBackward(U[nu], nu, Impl::CovShiftBackward(tmp, mu, U[nu])), mu); + // 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[nu]; - C3m[mu]+= Impl::ShiftStaple(Impl::CovShiftBackward(U[nu], nu, Impl::CovShiftBackward(U[mu], mu, tmp)), mu); + tmp = Lambda * U[mu]; + C2m[mu] += Impl::ShiftStaple(Impl::CovShiftBackward(U[nu], nu, Impl::CovShiftBackward(tmp, mu, U[nu])), mu); - tmp = Lambda; - C4m[mu]+= Impl::ShiftStaple(Impl::CovShiftBackward(U[nu], nu, Impl::CovShiftBackward(U[mu], mu, U[nu])), mu)* tmp; + 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 - - - + //Still implementing. Have to be tested, and understood how to project EO } // Derivative parts diff --git a/lib/qcd/action/fermion/WilsonCloverFermion.h b/lib/qcd/action/fermion/WilsonCloverFermion.h index 18386485..34482941 100644 --- a/lib/qcd/action/fermion/WilsonCloverFermion.h +++ b/lib/qcd/action/fermion/WilsonCloverFermion.h @@ -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 public: // Types definitions INHERIT_IMPL_TYPES(Impl); - template using iImplClover = iScalar, Ns> >; - typedef iImplClover SiteCloverType; - typedef Lattice CloverFieldType; + template + using iImplClover = iScalar, Ns>>; + typedef iImplClover SiteCloverType; + typedef Lattice CloverFieldType; + public: typedef WilsonFermion WilsonBase; @@ -58,19 +61,21 @@ public: Fgrid, Hgrid, _mass, p), - CloverTerm(&Fgrid), - CloverTermInv(&Fgrid), - CloverTermEven(&Hgrid), - CloverTermOdd(&Hgrid), - CloverTermInvEven(&Hgrid), - CloverTermInvOdd(&Hgrid), - CloverTermDagEven(&Hgrid), //test - CloverTermDagOdd(&Hgrid), //test - CloverTermInvDagEven(&Hgrid), //test - CloverTermInvDagOdd(&Hgrid) //test + 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); @@ -91,12 +96,12 @@ public: private: // here fixing the 4 dimensions, make it more general? - RealD csw; // Clover coefficient - CloverFieldType CloverTerm=zero, CloverTermInv=zero; // Clover term - CloverFieldType CloverTermEven=zero, CloverTermOdd=zero; // Clover term EO - CloverFieldType CloverTermInvEven=zero, CloverTermInvOdd=zero; // Clover term Inv EO - CloverFieldType CloverTermDagEven=zero, CloverTermDagOdd=zero; // Clover term Dag EO - CloverFieldType CloverTermInvDagEven=zero, CloverTermInvDagOdd=zero; // Clover term Inv Dag EO + 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 @@ -113,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) { @@ -129,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) { @@ -145,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) { @@ -156,14 +161,14 @@ private: PARALLEL_FOR_LOOP for (int i = 0; i < CloverTerm._grid->oSites(); 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]()()); + 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) { @@ -172,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) { @@ -188,17 +193,16 @@ private: PARALLEL_FOR_LOOP for (int i = 0; i < CloverTerm._grid->oSites(); 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]()()); + 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 diff --git a/tests/core/Test_wilson_clover.cc b/tests/core/Test_wilson_clover.cc index 9a5fffe2..9a55f6b2 100644 --- a/tests/core/Test_wilson_clover.cc +++ b/tests/core/Test_wilson_clover.cc @@ -237,7 +237,7 @@ int main (int argc, char ** argv) setCheckerboard(src,src_o); - //Gauge Transformation + ////////////////////// Gauge Transformation std::vector seeds2({5,6,7,8}); GridParallelRNG pRNG2(&Grid); pRNG2.SeedFixedIntegers(seeds2); LatticeColourMatrix Omega(&Grid); @@ -251,7 +251,7 @@ int main (int argc, char ** argv) 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); @@ -298,7 +298,7 @@ int main (int argc, char ** argv) std::cout< +#include +#include +#include // Mass -double mq = 0.0; +double mq = 0.01; // Define Wilson Types typedef Grid::QCD::WilsonImplR::FermionField FermionField; typedef Grid::QCD::LatticeGaugeField GaugeField; -#include -#include -#include - enum ChromaAction { Wilson, // Wilson 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 { @@ -286,91 +281,6 @@ public: }; } // 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 gd = UGrid->GlobalDimensions(); - QDP::multi1d 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 ActionList({Wilson, WilsonClover}); - std::vector 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) { QDP::multi1d 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) @@ -512,3 +421,76 @@ void calc_grid(ChromaAction action, Grid::QCD::LatticeGaugeField &Umu, Grid::QCD 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 gd = UGrid->GlobalDimensions(); + QDP::multi1d 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 ActionList({Wilson, WilsonClover}); + std::vector 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(); + } + } +}