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

Ward Identity test improvements and conserved current bug fixes

Browse Source
This commit is contained in:
Lanny91 2017-06-05 11:55:41 +01:00
parent 23135aa58a
commit eec79e0a1e
3 changed files with 42 additions and 24 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

35
extras/Hadrons/Modules/MContraction/WardIdentity.hpp
View File

@ -143,6 +143,7 @@ void TWardIdentity<FImpl>::execute(void)
vector_WI += tmp;
}
// Test ward identity D_mu V_mu = 0;
LOG(Message) << "Vector Ward Identity check Delta_mu V_mu = "
<< norm2(vector_WI) << std::endl;
@ -150,28 +151,29 @@ void TWardIdentity<FImpl>::execute(void)
{
LatticeComplex PP(env().getGrid()), axial_defect(env().getGrid()),
PJ5q(env().getGrid());
std::vector<TComplex> axial_buf;
// Compute D_mu A_mu, D is backwards derivative.
// Compute <P|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);
axial_defect += 2.*trace(g5*tmp);
}
// Get PJ5q for 5D (zero for 4D) and PP.
// Get <P|J5q> for 5D (zero for 4D) and <P|P>.
PJ5q = zero;
if (Ls_ > 1)
{
// PP
// <P|P>
ExtractSlice(tmp, q, 0, 0);
psi = (tmp - g5*tmp);
psi = 0.5 * (tmp - g5*tmp);
ExtractSlice(tmp, q, Ls_ - 1, 0);
psi += (tmp + g5*tmp);
psi += 0.5 * (tmp + g5*tmp);
PP = trace(adj(psi)*psi);
// P5Jq
// <P|5Jq>
ExtractSlice(tmp, q, Ls_/2 - 1, 0);
psi = 0.5 * (tmp + g5*tmp);
ExtractSlice(tmp, q, Ls_/2, 0);
@ -183,13 +185,22 @@ void TWardIdentity<FImpl>::execute(void)
PP = trace(adj(q)*q);
}
// Test ward identities, D_mu V_mu = 0; D_mu A_mu = 2m<PP> + 2 PJ5q
axial_defect -= 2.*PJ5q;
axial_defect -= 2.*(par().mass)*PP;
// Test ward identity <P|D_mu A_mu> = 2m<P|P> + 2<P|J5q>
LOG(Message) << "|D_mu A_mu|^2 = " << norm2(axial_defect) << std::endl;
LOG(Message) << "|PP|^2 = " << norm2(PP) << std::endl;
LOG(Message) << "|PJ5q|^2 = " << norm2(PJ5q) << std::endl;
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;
// Axial defect by timeslice.
axial_defect -= 2.*(par().mass*PP + PJ5q);
LOG(Message) << "Check Axial defect by timeslice" << std::endl;
sliceSum(axial_defect, axial_buf, Tp);
for (int t = 0; t < axial_buf.size(); ++t)
{
LOG(Message) << "t = " << t << ": "
<< TensorRemove(axial_buf[t]) << std::endl;
}
}
}

10
lib/qcd/action/fermion/WilsonFermion5D.cc
View File

@ -700,8 +700,14 @@ void WilsonFermion5D<Impl>::ContractConservedCurrent(PropagatorField &q_in_1,
Kernels::ContractConservedCurrentInternal(q1_s, q2_s, tmp, Umu, curr_type, mu);
// Axial current sign
Real G_s = (curr_type == Current::Axial) ? ((s < Ls/2) ? -1. : 1.) : 1.;
q_out += G_s*tmp;
if ((curr_type == Current::Axial) && (s < (Ls / 2)))
{
q_out -= tmp;
}
else
{
q_out += tmp;
}
}
}

21
lib/qcd/action/fermion/WilsonKernels.cc
View File

@ -286,8 +286,9 @@ void WilsonKernels<Impl>::DhopDir( StencilImpl &st, DoubledGaugeField &U,SiteHal
* to make a conserved current sink or inserting the conserved current
* sequentially. Common to both 4D and 5D.
******************************************************************************/
#define WilsonCurrentFwd(expr, mu) (0.5*(Gamma::gmu[mu]*expr - expr))
#define WilsonCurrentBwd(expr, mu) (0.5*(Gamma::gmu[mu]*expr + expr))
// N.B. Functions below assume a -1/2 factor within U.
#define WilsonCurrentFwd(expr, mu) ((expr - Gamma::gmu[mu]*expr))
#define WilsonCurrentBwd(expr, mu) ((expr + Gamma::gmu[mu]*expr))
template<class Impl>
void WilsonKernels<Impl>::ContractConservedCurrentInternal(const PropagatorField &q_in_1,
@ -300,13 +301,13 @@ void WilsonKernels<Impl>::ContractConservedCurrentInternal(const PropagatorField
Gamma g5(Gamma::Algebra::Gamma5);
PropagatorField tmp(q_out._grid);
GaugeLinkField Umu(U._grid);
Umu = PeekIndex<LorentzIndex>(U, mu);
Umu = PeekIndex<LorentzIndex>(U, mu);
tmp = this->CovShiftForward(Umu, mu, q_in_1);
q_out = (g5*adj(q_in_2)*g5)*WilsonCurrentFwd(tmp, mu);
tmp = adj(Umu)*q_in_1;
q_out += (g5*adj(this->CovShiftForward(Umu, mu, q_in_2))*g5)*WilsonCurrentBwd(q_in_1, mu);
tmp = this->CovShiftForward(Umu, mu, q_in_2);
q_out -= (g5*adj(tmp)*g5)*WilsonCurrentBwd(q_in_1, mu);
}
@ -320,21 +321,21 @@ void WilsonKernels<Impl>::SeqConservedCurrentInternal(const PropagatorField &q_i
unsigned int tmin,
unsigned int tmax)
{
int tshift = (mu == Nd - 1) ? 1 : 0;
int tshift = (mu == Tp) ? 1 : 0;
Real G_T = (curr_type == Current::Tadpole) ? -1. : 1.;
PropagatorField tmp(q_in._grid);
GaugeLinkField Umu(U._grid);
Umu = PeekIndex<LorentzIndex>(U, mu);
Umu = PeekIndex<LorentzIndex>(U, mu);
Lattice<iScalar<vInteger>> t(q_in._grid);
tmp = this->CovShiftForward(Umu, mu, q_in)*ph;
where((t >= tmin) and (t <= tmax), tmp, 0.*tmp);
tmp = where((t >= tmin) and (t <= tmax), tmp, 0.*tmp);
q_out = G_T*WilsonCurrentFwd(tmp, mu);
tmp = q_in*ph;
tmp = this->CovShiftBackward(Umu, mu, tmp);
where((t >= tmin + tshift) and (t <= tmax + tshift), tmp, 0.*tmp);
q_out += WilsonCurrentBwd(tmp, mu);
tmp = where((t >= tmin + tshift) and (t <= tmax + tshift), tmp, 0.*tmp);
q_out -= WilsonCurrentBwd(tmp, mu);
}