From 08b314fd0fadd012492710bf57fd17b37ea9cf54 Mon Sep 17 00:00:00 2001 From: Lanny91 Date: Thu, 18 May 2017 13:16:14 +0100 Subject: [PATCH] Hadrons: conserved current test fixes. Axial current tests now also optional. --- .../Modules/MContraction/WardIdentity.hpp | 106 ++++++++++-------- tests/hadrons/Test_hadrons.hpp | 13 ++- .../hadrons/Test_hadrons_conserved_current.cc | 3 +- 3 files changed, 69 insertions(+), 53 deletions(-) diff --git a/extras/Hadrons/Modules/MContraction/WardIdentity.hpp b/extras/Hadrons/Modules/MContraction/WardIdentity.hpp index d312bd4d..fa51ce95 100644 --- a/extras/Hadrons/Modules/MContraction/WardIdentity.hpp +++ b/extras/Hadrons/Modules/MContraction/WardIdentity.hpp @@ -40,10 +40,10 @@ BEGIN_HADRONS_NAMESPACE ----------------------------- * options: - - q: propagator, 5D if available (string) - - q4d: 4D propagator, duplicate of q if q is not 5D (string) - - action: action module used for propagator solution (string) - - mass: mass of quark (double) + - q: propagator, 5D if available (string) + - action: action module used for propagator solution (string) + - mass: mass of quark (double) + - test_axial: whether or not to test PCAC relation. */ /****************************************************************************** @@ -56,9 +56,9 @@ class WardIdentityPar: Serializable public: GRID_SERIALIZABLE_CLASS_MEMBERS(WardIdentityPar, std::string, q, - std::string, q4d, std::string, action, - double, mass); + double, mass, + bool, test_axial); }; template @@ -97,7 +97,7 @@ TWardIdentity::TWardIdentity(const std::string name) template std::vector TWardIdentity::getInput(void) { - std::vector in = {par().q, par().q4d, par().action}; + std::vector in = {par().q, par().action}; return in; } @@ -128,55 +128,69 @@ void TWardIdentity::execute(void) LOG(Message) << "Performing Ward Identity checks for quark '" << par().q << "'." << std::endl; - PropagatorField psi(env().getGrid()), tmp(env().getGrid()); + PropagatorField psi(env().getGrid()), tmp(env().getGrid()), + vector_WI(env().getGrid()); PropagatorField &q = *env().template getObject(par().q); - PropagatorField &q4d = *env().template getObject(par().q4d); FMat &act = *(env().template getObject(par().action)); Gamma g5(Gamma::Algebra::Gamma5); - LatticeComplex PP(env().getGrid()), PA(env().getGrid()), - c(env().getGrid()), PJ5q(env().getGrid()), - vector_WI(env().getGrid()), defect(env().getGrid()); - c = zero; PJ5q = zero; vector_WI = zero; defect = zero; - std::vector Vmu(Nd, c); - std::vector Amu(Nd, c); - - // Get PP, PA, V_mu, A_mu for 4D. - PP = trace(adj(q4d)*q4d); - PA = trace(adj(q4d)*g5*q4d); + + // Compute D_mu V_mu, D here is backward derivative. + vector_WI = zero; for (unsigned int mu = 0; mu < Nd; ++mu) { act.ContractConservedCurrent(q, q, tmp, Current::Vector, mu); - Vmu[mu] = trace(tmp); - act.ContractConservedCurrent(q, q, tmp, Current::Axial, mu); - Amu[mu] = trace(g5*tmp); + tmp -= Cshift(tmp, mu, -1); + vector_WI += tmp; } - // Get PJ5q for 5D (zero for 4D). - if (Ls_ > 1) - { - ExtractSlice(psi, q, Ls_/2 - 1, 0); - psi = 0.5 * (psi + g5*psi); - ExtractSlice(tmp, q, Ls_/2, 0); - psi += 0.5 * (tmp - g5*tmp); - PJ5q = trace(adj(psi)*psi); - } - - // Test ward identities, D_mu V_mu = 0; D_mu A_mu = 2m + 2 PJ5q - for (unsigned int mu = 0; mu < Nd; ++mu) - { - vector_WI += Vmu[mu] - Cshift(Vmu[mu], mu, -1); - defect += Amu[mu] - Cshift(Amu[mu], mu, -1); - } - defect -= 2.*PJ5q; - defect -= 2.*(par().mass)*PP; - LOG(Message) << "Vector Ward Identity check Delta_mu V_mu = " << norm2(vector_WI) << std::endl; - LOG(Message) << "Axial Ward Identity defect Delta_mu A_mu = " - << norm2(defect) << std::endl; - LOG(Message) << "norm2(PP) = " << norm2(PP) << std::endl; - LOG(Message) << "norm2(PA) = " << norm2(PA) << std::endl; - LOG(Message) << "norm2(PJ5q) = " << norm2(PJ5q) << std::endl; + + if (par().test_axial) + { + LatticeComplex PP(env().getGrid()), axial_defect(env().getGrid()), + PJ5q(env().getGrid()); + + // Compute D_mu A_mu, D is backwards derivative. + axial_defect = zero; + for (unsigned int mu = 0; mu < Nd; ++mu) + { + act.ContractConservedCurrent(q, q, tmp, Current::Axial, mu); + tmp -= Cshift(tmp, mu, -1); + axial_defect += trace(g5*tmp); + } + + // Get PJ5q for 5D (zero for 4D) and PP. + PJ5q = zero; + if (Ls_ > 1) + { + // PP + ExtractSlice(tmp, q, 0, 0); + psi = (tmp - g5*tmp); + ExtractSlice(tmp, q, Ls_ - 1, 0); + psi += (tmp + g5*tmp); + PP = trace(adj(psi)*psi); + + // P5Jq + ExtractSlice(tmp, q, Ls_/2 - 1, 0); + psi = 0.5 * (tmp + g5*tmp); + ExtractSlice(tmp, q, Ls_/2, 0); + psi += 0.5 * (tmp - g5*tmp); + PJ5q = trace(adj(psi)*psi); + } + else + { + PP = trace(adj(q)*q); + } + + // Test ward identities, D_mu V_mu = 0; D_mu A_mu = 2m + 2 PJ5q + axial_defect -= 2.*PJ5q; + axial_defect -= 2.*(par().mass)*PP; + LOG(Message) << "Axial Ward Identity defect Delta_mu A_mu = " + << norm2(axial_defect) << std::endl; + LOG(Message) << "norm2(PP) = " << norm2(PP) << std::endl; + LOG(Message) << "norm2(PJ5q) = " << norm2(PJ5q) << std::endl; + } } END_MODULE_NAMESPACE diff --git a/tests/hadrons/Test_hadrons.hpp b/tests/hadrons/Test_hadrons.hpp index 1b038388..6dbe3425 100644 --- a/tests/hadrons/Test_hadrons.hpp +++ b/tests/hadrons/Test_hadrons.hpp @@ -513,26 +513,27 @@ inline void discLoopContraction(Application &application, * actionName - action used to compute quark propagator. * mass - mass of quark. * Ls - length of 5th dimension (default = 1). + * test_axial - whether or not to check PCAC relation. * Returns: None. ******************************************************************************/ inline void makeWITest(Application &application, std::string &modName, std::string &propName, std::string &actionName, - double mass, unsigned int Ls = 1) + double mass, unsigned int Ls = 1, bool test_axial = false) { if (!(Environment::getInstance().hasModule(modName))) { MContraction::WardIdentity::Par wiPar; if (Ls > 1) { - wiPar.q = LABEL_5D(propName); + wiPar.q = LABEL_5D(propName); } else { - wiPar.q = propName; + wiPar.q = propName; } - wiPar.q4d = propName; - wiPar.action = actionName; - wiPar.mass = mass; + wiPar.action = actionName; + wiPar.mass = mass; + wiPar.test_axial = test_axial; application.createModule(modName, wiPar); } } diff --git a/tests/hadrons/Test_hadrons_conserved_current.cc b/tests/hadrons/Test_hadrons_conserved_current.cc index 080fef73..37ef30d9 100644 --- a/tests/hadrons/Test_hadrons_conserved_current.cc +++ b/tests/hadrons/Test_hadrons_conserved_current.cc @@ -81,7 +81,8 @@ inline void setupWardIdentityTests(Application &application, std::string origin = "0 0 0 0"; std::string modName = actionName + " Ward Identity Test"; MAKE_POINT_PROP(origin, pointProp, solverName); - makeWITest(application, modName, pointProp, actionName, mass, Ls); + makeWITest(application, modName, pointProp, actionName, mass, Ls, + perform_axial_tests); /*************************************************************************** * Conserved current tests with sequential insertion of vector/axial