From ac1253bb764fe3424f1afa87dea4ed349ef24dd9 Mon Sep 17 00:00:00 2001 From: Lanny91 Date: Mon, 10 Apr 2017 17:42:55 +0100 Subject: [PATCH 01/50] Corrected solver in rare kaon test --- tests/hadrons/Test_hadrons_rarekaon.cc | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/tests/hadrons/Test_hadrons_rarekaon.cc b/tests/hadrons/Test_hadrons_rarekaon.cc index 26622525..7c76312d 100644 --- a/tests/hadrons/Test_hadrons_rarekaon.cc +++ b/tests/hadrons/Test_hadrons_rarekaon.cc @@ -109,11 +109,10 @@ int main(int argc, char *argv[]) application.createModule("DWF_" + flavour[i], actionPar); // solvers - // RBPrecCG -> CG MSolver::CG::Par solverPar; solverPar.action = "DWF_" + flavour[i]; solverPar.residual = 1.0e-8; - application.createModule(solvers[i], + application.createModule(solvers[i], solverPar); } From af2d6ce2e08d54ea7a46e5eeb56e0374905d9343 Mon Sep 17 00:00:00 2001 From: Lanny91 Date: Fri, 10 Mar 2017 14:59:11 +0000 Subject: [PATCH 02/50] Encapsulated 4D->5D and 5D->4D conversions in separate functions & added corresponding tests. --- extras/Hadrons/Modules/Quark.hpp | 32 ++++-- tests/hadrons/Test_hadrons_quark.cc | 156 ++++++++++++++++++++++++++++ 2 files changed, 179 insertions(+), 9 deletions(-) create mode 100644 tests/hadrons/Test_hadrons_quark.cc diff --git a/extras/Hadrons/Modules/Quark.hpp b/extras/Hadrons/Modules/Quark.hpp index be7426ab..c08e0192 100644 --- a/extras/Hadrons/Modules/Quark.hpp +++ b/extras/Hadrons/Modules/Quark.hpp @@ -36,6 +36,27 @@ See the full license in the file "LICENSE" in the top level distribution directo BEGIN_HADRONS_NAMESPACE +/****************************************************************************** + * 5D -> 4D and 4D -> 5D conversions. * + ******************************************************************************/ +template // Note that 5D object is modified. +inline void make_4D(Lattice &in_5d, Lattice &out_4d, int Ls) +{ + axpby_ssp_pminus(in_5d, 0., in_5d, 1., in_5d, 0, 0); + axpby_ssp_pplus(in_5d, 1., in_5d, 1., in_5d, 0, Ls-1); + ExtractSlice(out_4d, in_5d, 0, 0); +} + +template +inline void make_5D(const Lattice &in_4d, Lattice &out_5d, int Ls) +{ + out_5d = zero; + InsertSlice(in_4d, out_5d, 0, 0); + InsertSlice(in_4d, out_5d, Ls-1, 0); + axpby_ssp_pplus(out_5d, 0., out_5d, 1., out_5d, 0, 0); + axpby_ssp_pminus(out_5d, 0., out_5d, 1., out_5d, Ls-1, Ls-1); +} + /****************************************************************************** * TQuark * ******************************************************************************/ @@ -143,12 +164,8 @@ void TQuark::execute(void) } else { - source = zero; PropToFerm(tmp, fullSrc, s, c); - InsertSlice(tmp, source, 0, 0); - InsertSlice(tmp, source, Ls_-1, 0); - axpby_ssp_pplus(source, 0., source, 1., source, 0, 0); - axpby_ssp_pminus(source, 0., source, 1., source, Ls_-1, Ls_-1); + make_5D(tmp, source, Ls_); } } // source conversion for 5D sources @@ -171,10 +188,7 @@ void TQuark::execute(void) { PropagatorField &p4d = *env().template getObject(getName()); - - axpby_ssp_pminus(sol, 0., sol, 1., sol, 0, 0); - axpby_ssp_pplus(sol, 1., sol, 1., sol, 0, Ls_-1); - ExtractSlice(tmp, sol, 0, 0); + make_4D(sol, tmp, Ls_); FermToProp(p4d, tmp, s, c); } } diff --git a/tests/hadrons/Test_hadrons_quark.cc b/tests/hadrons/Test_hadrons_quark.cc new file mode 100644 index 00000000..6a142ff6 --- /dev/null +++ b/tests/hadrons/Test_hadrons_quark.cc @@ -0,0 +1,156 @@ +/******************************************************************************* + Grid physics library, www.github.com/paboyle/Grid + + Source file: tests/hadrons/Test_hadrons_quark.cc + + Copyright (C) 2017 + + Author: Andrew Lawson + + This program is free software; you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation; either version 2 of the License, or + (at your option) any later version. + + This program is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License along + with this program; if not, write to the Free Software Foundation, Inc., + 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. + + See the full license in the file "LICENSE" in the top level distribution + directory. + *******************************************************************************/ + +#include "Test_hadrons.hpp" +#include + +using namespace Grid; +using namespace QCD; +using namespace Hadrons; + +/******************************************************************************* + * Unit test functions within Quark module. + ******************************************************************************/ + +// Alternative 4D & 5D projections +template +inline void make_4D_with_gammas(Lattice &in_5d, Lattice &out_4d, int Ls) +{ + GridBase *_grid(out_4d._grid); + Lattice tmp(_grid); + Gamma G5(Gamma::Algebra::Gamma5); + + ExtractSlice(tmp, in_5d, 0, 0); + out_4d = 0.5 * (tmp - G5*tmp); + ExtractSlice(tmp, in_5d, Ls - 1, 0); + out_4d += 0.5 * (tmp + G5*tmp); +} + +template +inline void make_5D_with_gammas(Lattice &in_4d, Lattice &out_5d, int Ls) +{ + out_5d = zero; + Gamma G5(Gamma::Algebra::Gamma5); + GridBase *_grid(in_4d._grid); + Lattice tmp(_grid); + + tmp = 0.5 * (in_4d + G5*in_4d); + InsertSlice(tmp, out_5d, 0, 0); + tmp = 0.5 * (in_4d - G5*in_4d); + InsertSlice(tmp, out_5d, Ls - 1, 0); +} + +int main(int argc, char **argv) +{ + /*************************************************************************** + * Initialisation. + **************************************************************************/ + Grid_init(&argc, &argv); + + std::vector latt_size = GridDefaultLatt(); + std::vector simd_layout = GridDefaultSimd(Nd,vComplex::Nsimd()); + std::vector mpi_layout = GridDefaultMpi(); + + const int Ls = 8; + + GridCartesian UGrid(latt_size,simd_layout,mpi_layout); + GridCartesian *FGrid = SpaceTimeGrid::makeFiveDimGrid(Ls, &UGrid); + GridSerialRNG sRNG; + GridParallelRNG pRNG(&UGrid); + + std::vector seeds4({1,2,3,4}); + std::vector seeds5({5,6,7,8}); + GridParallelRNG rng4(&UGrid); + GridParallelRNG rng5(FGrid); + rng4.SeedFixedIntegers(seeds4); + rng5.SeedFixedIntegers(seeds5); + + /*************************************************************************** + * Build a 4D random source, and convert it to 5D. + **************************************************************************/ + LatticeFermion test4(&UGrid); + LatticeFermion test5(FGrid); + LatticeFermion check5(FGrid); + + gaussian(rng4, test4); + make_5D(test4, test5, Ls); + make_5D_with_gammas(test4, check5, Ls); + test5 -= check5; + std::cout << "4D -> 5D comparison, diff = " << Grid::sqrt(norm2(test5)) << std::endl; + + /*************************************************************************** + * Build a 5D random source, and project down to 4D. + **************************************************************************/ + LatticeFermion check4(&UGrid); + gaussian(rng5, test5); + check5 = test5; + + make_4D(test5, test4, Ls); + make_4D_with_gammas(check5, check4, Ls); + test4 -= check4; + std::cout << "5D -> 4D comparison, diff = " << Grid::sqrt(norm2(test4)) << std::endl; + + /*************************************************************************** + * Convert a propagator to a fermion & back. + **************************************************************************/ + LatticeFermion ferm(&UGrid); + LatticePropagator prop(&UGrid), ref(&UGrid); + gaussian(rng4, prop); + + // Define variables for sanity checking a single site. + typename SpinColourVector::scalar_object fermSite; + typename SpinColourMatrix::scalar_object propSite; + std::vector site(Nd, 0); + + for (int s = 0; s < Ns; ++s) + for (int c = 0; c < Nc; ++c) + { + ref = prop; + PropToFerm(ferm, prop, s, c); + FermToProp(prop, ferm, s, c); + + std::cout << "Spin = " << s << ", Colour = " << c << std::endl; + ref -= prop; + std::cout << "Prop->Ferm->Prop test, diff = " << Grid::sqrt(norm2(ref)) << std::endl; + + peekSite(fermSite, ferm, site); + peekSite(propSite, prop, site); + for (int s2 = 0; s2 < Ns; ++s2) + for (int c2 = 0; c2 < Nc; ++c2) + { + //if (propSite()(s2, s)(c2, c) != fermSite()(s2)(c2)) + //{ + std::cout << propSite()(s2, s)(c2, c) << " != " + << fermSite()(s2)(c2) << " for spin = " << s2 + << ", col = " << c2 << std::endl; + //} + } + } + + Grid_finalize(); + return EXIT_SUCCESS; +} From c382c351a5d088f6af5734ab5dbc11aec53e3cf4 Mon Sep 17 00:00:00 2001 From: Lanny91 Date: Fri, 10 Mar 2017 15:05:59 +0000 Subject: [PATCH 03/50] Quark test output correction. --- tests/hadrons/Test_hadrons_quark.cc | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/tests/hadrons/Test_hadrons_quark.cc b/tests/hadrons/Test_hadrons_quark.cc index 6a142ff6..5b9d0ce1 100644 --- a/tests/hadrons/Test_hadrons_quark.cc +++ b/tests/hadrons/Test_hadrons_quark.cc @@ -142,12 +142,12 @@ int main(int argc, char **argv) for (int s2 = 0; s2 < Ns; ++s2) for (int c2 = 0; c2 < Nc; ++c2) { - //if (propSite()(s2, s)(c2, c) != fermSite()(s2)(c2)) - //{ + if (propSite()(s2, s)(c2, c) != fermSite()(s2)(c2)) + { std::cout << propSite()(s2, s)(c2, c) << " != " << fermSite()(s2)(c2) << " for spin = " << s2 << ", col = " << c2 << std::endl; - //} + } } } From 1425afc72feae364a2629f30eb6c34783d6374eb Mon Sep 17 00:00:00 2001 From: Lanny91 Date: Tue, 25 Apr 2017 17:26:56 +0100 Subject: [PATCH 04/50] Rare Kaon test fix --- tests/hadrons/Test_hadrons_rarekaon.cc | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tests/hadrons/Test_hadrons_rarekaon.cc b/tests/hadrons/Test_hadrons_rarekaon.cc index 7c76312d..9d35c1bc 100644 --- a/tests/hadrons/Test_hadrons_rarekaon.cc +++ b/tests/hadrons/Test_hadrons_rarekaon.cc @@ -109,7 +109,7 @@ int main(int argc, char *argv[]) application.createModule("DWF_" + flavour[i], actionPar); // solvers - MSolver::CG::Par solverPar; + MSolver::RBPrecCG::Par solverPar; solverPar.action = "DWF_" + flavour[i]; solverPar.residual = 1.0e-8; application.createModule(solvers[i], From 44260643f6b2d61d8ea6a2543c67e08a76dab748 Mon Sep 17 00:00:00 2001 From: Lanny91 Date: Tue, 25 Apr 2017 18:00:24 +0100 Subject: [PATCH 05/50] First conserved current implementation for Wilson fermions only. Not implemented for Gparity or 5D-vectorised Wilson fermions. --- lib/qcd/QCD.h | 8 ++ lib/qcd/action/fermion/FermionOperator.h | 15 ++++ .../fermion/ImprovedStaggeredFermion.cc | 25 ++++++ .../action/fermion/ImprovedStaggeredFermion.h | 16 ++++ .../fermion/ImprovedStaggeredFermion5D.cc | 24 +++++ .../fermion/ImprovedStaggeredFermion5D.h | 15 ++++ lib/qcd/action/fermion/WilsonFermion.cc | 47 ++++++++++ lib/qcd/action/fermion/WilsonFermion.h | 16 ++++ lib/qcd/action/fermion/WilsonFermion5D.cc | 77 ++++++++++++++++ lib/qcd/action/fermion/WilsonFermion5D.h | 15 ++++ lib/qcd/action/fermion/WilsonKernels.cc | 89 +++++++++++++++++++ lib/qcd/action/fermion/WilsonKernels.h | 18 ++++ 12 files changed, 365 insertions(+) diff --git a/lib/qcd/QCD.h b/lib/qcd/QCD.h index 6e6144da..c66c7b13 100644 --- a/lib/qcd/QCD.h +++ b/lib/qcd/QCD.h @@ -489,6 +489,14 @@ namespace QCD { return traceIndex(lhs); } + ////////////////////////////////////////// + // Current types + ////////////////////////////////////////// + GRID_SERIALIZABLE_ENUM(Current, undef, + Vector, 0, + Axial, 1, + Tadpole, 2); + } //namespace QCD } // Grid diff --git a/lib/qcd/action/fermion/FermionOperator.h b/lib/qcd/action/fermion/FermionOperator.h index 676a0e83..144b70f6 100644 --- a/lib/qcd/action/fermion/FermionOperator.h +++ b/lib/qcd/action/fermion/FermionOperator.h @@ -112,6 +112,21 @@ namespace Grid { /////////////////////////////////////////////// virtual void ImportGauge(const GaugeField & _U)=0; + ////////////////////////////////////////////////////////////////////// + // Conserved currents, either contract at sink or insert sequentially. + ////////////////////////////////////////////////////////////////////// + virtual void ContractConservedCurrent(PropagatorField &q_in_1, + PropagatorField &q_in_2, + PropagatorField &q_out, + Current curr_type, + unsigned int mu)=0; + virtual void SeqConservedCurrent(PropagatorField &q_in, + PropagatorField &q_out, + Current curr_type, + unsigned int mu, + std::vector mom, + unsigned int tmin, + unsigned int tmax)=0; }; } diff --git a/lib/qcd/action/fermion/ImprovedStaggeredFermion.cc b/lib/qcd/action/fermion/ImprovedStaggeredFermion.cc index 2ba4f4af..ef8c79bd 100644 --- a/lib/qcd/action/fermion/ImprovedStaggeredFermion.cc +++ b/lib/qcd/action/fermion/ImprovedStaggeredFermion.cc @@ -395,6 +395,31 @@ void ImprovedStaggeredFermion::DhopInternal(StencilImpl &st, LebesgueOrder } }; +//////////////////////////////////////////////////////// +// Conserved current - not yet implemented. +//////////////////////////////////////////////////////// +template +void ImprovedStaggeredFermion::ContractConservedCurrent(PropagatorField &q_in_1, + PropagatorField &q_in_2, + PropagatorField &q_out, + Current curr_type, + unsigned int mu) +{ + assert(0); +} + +template +void ImprovedStaggeredFermion::SeqConservedCurrent(PropagatorField &q_in, + PropagatorField &q_out, + Current curr_type, + unsigned int mu, + std::vector mom, + unsigned int tmin, + unsigned int tmax) +{ + assert(0); +} + FermOpStaggeredTemplateInstantiate(ImprovedStaggeredFermion); //AdjointFermOpTemplateInstantiate(ImprovedStaggeredFermion); diff --git a/lib/qcd/action/fermion/ImprovedStaggeredFermion.h b/lib/qcd/action/fermion/ImprovedStaggeredFermion.h index 7d1f2996..9d5270c6 100644 --- a/lib/qcd/action/fermion/ImprovedStaggeredFermion.h +++ b/lib/qcd/action/fermion/ImprovedStaggeredFermion.h @@ -157,6 +157,22 @@ class ImprovedStaggeredFermion : public StaggeredKernels, public ImprovedS LebesgueOrder Lebesgue; LebesgueOrder LebesgueEvenOdd; + + /////////////////////////////////////////////////////////////// + // Conserved current utilities + /////////////////////////////////////////////////////////////// + void ContractConservedCurrent(PropagatorField &q_in_1, + PropagatorField &q_in_2, + PropagatorField &q_out, + Current curr_type, + unsigned int mu); + void SeqConservedCurrent(PropagatorField &q_in, + PropagatorField &q_out, + Current curr_type, + unsigned int mu, + std::vector mom, + unsigned int tmin, + unsigned int tmax); }; typedef ImprovedStaggeredFermion ImprovedStaggeredFermionF; diff --git a/lib/qcd/action/fermion/ImprovedStaggeredFermion5D.cc b/lib/qcd/action/fermion/ImprovedStaggeredFermion5D.cc index 61a3c559..293077f7 100644 --- a/lib/qcd/action/fermion/ImprovedStaggeredFermion5D.cc +++ b/lib/qcd/action/fermion/ImprovedStaggeredFermion5D.cc @@ -345,6 +345,30 @@ void ImprovedStaggeredFermion5D::MooeeInvDag(const FermionField &in, MooeeInv(in, out); } +//////////////////////////////////////////////////////// +// Conserved current - not yet implemented. +//////////////////////////////////////////////////////// +template +void ImprovedStaggeredFermion5D::ContractConservedCurrent(PropagatorField &q_in_1, + PropagatorField &q_in_2, + PropagatorField &q_out, + Current curr_type, + unsigned int mu) +{ + assert(0); +} + +template +void ImprovedStaggeredFermion5D::SeqConservedCurrent(PropagatorField &q_in, + PropagatorField &q_out, + Current curr_type, + unsigned int mu, + std::vector mom, + unsigned int tmin, + unsigned int tmax) +{ + assert(0); +} FermOpStaggeredTemplateInstantiate(ImprovedStaggeredFermion5D); FermOpStaggeredVec5dTemplateInstantiate(ImprovedStaggeredFermion5D); diff --git a/lib/qcd/action/fermion/ImprovedStaggeredFermion5D.h b/lib/qcd/action/fermion/ImprovedStaggeredFermion5D.h index 4961da49..1c540892 100644 --- a/lib/qcd/action/fermion/ImprovedStaggeredFermion5D.h +++ b/lib/qcd/action/fermion/ImprovedStaggeredFermion5D.h @@ -160,6 +160,21 @@ namespace QCD { // Comms buffer std::vector > comm_buf; + /////////////////////////////////////////////////////////////// + // Conserved current utilities + /////////////////////////////////////////////////////////////// + void ContractConservedCurrent(PropagatorField &q_in_1, + PropagatorField &q_in_2, + PropagatorField &q_out, + Current curr_type, + unsigned int mu); + void SeqConservedCurrent(PropagatorField &q_in, + PropagatorField &q_out, + Current curr_type, + unsigned int mu, + std::vector mom, + unsigned int tmin, + unsigned int tmax); }; }} diff --git a/lib/qcd/action/fermion/WilsonFermion.cc b/lib/qcd/action/fermion/WilsonFermion.cc index 32083d5e..839f5215 100644 --- a/lib/qcd/action/fermion/WilsonFermion.cc +++ b/lib/qcd/action/fermion/WilsonFermion.cc @@ -347,6 +347,53 @@ void WilsonFermion::DhopInternal(StencilImpl &st, LebesgueOrder &lo, } }; +/******************************************************************************* + * Conserved current utilities for Wilson fermions, for contracting propagators + * to make a conserved current sink or inserting the conserved current + * sequentially. + ******************************************************************************/ +template +void WilsonFermion::ContractConservedCurrent(PropagatorField &q_in_1, + PropagatorField &q_in_2, + PropagatorField &q_out, + Current curr_type, + unsigned int mu) +{ + Gamma g5(Gamma::Algebra::Gamma5); + conformable(_grid, q_in_1._grid); + conformable(_grid, q_in_2._grid); + conformable(_grid, q_out._grid); + Kernels::ContractConservedCurrentInternal(q_in_1, q_in_2, q_out, + Umu, curr_type, mu); +} + +template +void WilsonFermion::SeqConservedCurrent(PropagatorField &q_in, + PropagatorField &q_out, + Current curr_type, + unsigned int mu, + std::vector mom, + unsigned int tmin, + unsigned int tmax) +{ + conformable(_grid, q_in._grid); + conformable(_grid, q_out._grid); + Lattice> ph(_grid), coor(_grid); + Complex i(0.0,1.0); + + // Momentum projection + ph = zero; + for(unsigned int mu = 0; mu < Nd - 1; mu++) + { + LatticeCoordinate(coor, mu); + ph = ph + mom[mu]*coor*((1./(_grid->_fdimensions[mu]))); + } + ph = exp((Real)(2*M_PI)*i*ph); + + Kernels::SeqConservedCurrentInternal(q_in, q_out, Umu, curr_type, mu, ph, + tmin, tmax); +} + FermOpTemplateInstantiate(WilsonFermion); AdjointFermOpTemplateInstantiate(WilsonFermion); TwoIndexFermOpTemplateInstantiate(WilsonFermion); diff --git a/lib/qcd/action/fermion/WilsonFermion.h b/lib/qcd/action/fermion/WilsonFermion.h index 933be732..feba40ed 100644 --- a/lib/qcd/action/fermion/WilsonFermion.h +++ b/lib/qcd/action/fermion/WilsonFermion.h @@ -146,6 +146,22 @@ class WilsonFermion : public WilsonKernels, public WilsonFermionStatic { LebesgueOrder Lebesgue; LebesgueOrder LebesgueEvenOdd; + + /////////////////////////////////////////////////////////////// + // Conserved current utilities + /////////////////////////////////////////////////////////////// + void ContractConservedCurrent(PropagatorField &q_in_1, + PropagatorField &q_in_2, + PropagatorField &q_out, + Current curr_type, + unsigned int mu); + void SeqConservedCurrent(PropagatorField &q_in, + PropagatorField &q_out, + Current curr_type, + unsigned int mu, + std::vector mom, + unsigned int tmin, + unsigned int tmax); }; typedef WilsonFermion WilsonFermionF; diff --git a/lib/qcd/action/fermion/WilsonFermion5D.cc b/lib/qcd/action/fermion/WilsonFermion5D.cc index 88bc425a..d0d3d055 100644 --- a/lib/qcd/action/fermion/WilsonFermion5D.cc +++ b/lib/qcd/action/fermion/WilsonFermion5D.cc @@ -679,6 +679,83 @@ void WilsonFermion5D::MomentumSpacePropagatorHw(FermionField &out,const Fe } +/******************************************************************************* + * Conserved current utilities for Wilson fermions, for contracting propagators + * to make a conserved current sink or inserting the conserved current + * sequentially. + ******************************************************************************/ +template +void WilsonFermion5D::ContractConservedCurrent(PropagatorField &q_in_1, + PropagatorField &q_in_2, + PropagatorField &q_out, + Current curr_type, + unsigned int mu) +{ + conformable(q_in_1._grid, FermionGrid()); + conformable(q_in_1._grid, q_in_2._grid); + conformable(_FourDimGrid, q_out._grid); + + PropagatorField q1_s(_FourDimGrid); + PropagatorField q2_s(_FourDimGrid); + PropagatorField tmp(_FourDimGrid); + + // Contract across 5th dimension. + q_out = zero; + for (int s = 0; s < Ls; ++s) + { + ExtractSlice(q1_s, q_in_1, 0, s); + ExtractSlice(q2_s, q_in_2, 0, Ls - s - 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; + } +} + + +template +void WilsonFermion5D::SeqConservedCurrent(PropagatorField &q_in, + PropagatorField &q_out, + Current curr_type, + unsigned int mu, + std::vector mom, + unsigned int tmin, + unsigned int tmax) +{ + conformable(q_in._grid, FermionGrid()); + conformable(q_in._grid, q_out._grid); + Lattice> ph(_FourDimGrid), coor(_FourDimGrid); + Complex i(0.0, 1.0); + + // Momentum projection + ph = zero; + for(unsigned int nu = 0; nu < Nd - 1; nu++) + { + LatticeCoordinate(coor, nu); + ph = ph + mom[nu]*coor*((1./(_FourDimGrid->_fdimensions[nu]))); + } + ph = exp((Real)(2*M_PI)*i*ph); + + // Sequential insertion + Kernels::SeqConservedCurrentInternal(q_in, q_out, Umu, curr_type, + mu, ph, tmin, tmax); + + // Axial current sign. + if (curr_type == Current::Axial) + { + SitePropagator result; + parallel_for(int sU = 0; sU < Umu._grid->oSites(); sU++) + { + int sF = sU * Ls; + for (int s = 0; s < Ls/2; s++) + { + vstream(q_out._odata[sF], -q_out._odata[sF]); + sF++; + } + } + } +} FermOpTemplateInstantiate(WilsonFermion5D); GparityFermOpTemplateInstantiate(WilsonFermion5D); diff --git a/lib/qcd/action/fermion/WilsonFermion5D.h b/lib/qcd/action/fermion/WilsonFermion5D.h index e87e927e..d66f4a1d 100644 --- a/lib/qcd/action/fermion/WilsonFermion5D.h +++ b/lib/qcd/action/fermion/WilsonFermion5D.h @@ -214,6 +214,21 @@ namespace QCD { // Comms buffer std::vector > comm_buf; + /////////////////////////////////////////////////////////////// + // Conserved current utilities + /////////////////////////////////////////////////////////////// + void ContractConservedCurrent(PropagatorField &q_in_1, + PropagatorField &q_in_2, + PropagatorField &q_out, + Current curr_type, + unsigned int mu); + void SeqConservedCurrent(PropagatorField &q_in, + PropagatorField &q_out, + Current curr_type, + unsigned int mu, + std::vector mom, + unsigned int tmin, + unsigned int tmax); }; }} diff --git a/lib/qcd/action/fermion/WilsonKernels.cc b/lib/qcd/action/fermion/WilsonKernels.cc index 6e72e089..fbf8dc00 100644 --- a/lib/qcd/action/fermion/WilsonKernels.cc +++ b/lib/qcd/action/fermion/WilsonKernels.cc @@ -567,6 +567,95 @@ void WilsonKernels::DhopDir( StencilImpl &st, DoubledGaugeField &U,SiteHal vstream(out._odata[sF], result); } +/******************************************************************************* + * Conserved current utilities for Wilson fermions, for contracting propagators + * 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)) + +template +void WilsonKernels::ContractConservedCurrentInternal(const PropagatorField &q_in_1, + const PropagatorField &q_in_2, + PropagatorField &q_out, + DoubledGaugeField &U, + Current curr_type, + unsigned int mu) +{ + Gamma g5(Gamma::Algebra::Gamma5); + PropagatorField tmp(q_out._grid); + GaugeLinkField Umu(U._grid); + Umu = PeekIndex(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); +} + + +template +void WilsonKernels::SeqConservedCurrentInternal(const PropagatorField &q_in, + PropagatorField &q_out, + DoubledGaugeField &U, + Current curr_type, + unsigned int mu, + Lattice> &ph, + unsigned int tmin, + unsigned int tmax) +{ + int tshift = (mu == Nd - 1) ? 1 : 0; + Real G_T = (curr_type == Current::Tadpole) ? -1. : 1.; + PropagatorField tmp(q_in._grid); + GaugeLinkField Umu(U._grid); + Umu = PeekIndex(U, mu); + Lattice> t(q_in._grid); + + tmp = this->CovShiftForward(Umu, mu, q_in)*ph; + 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); +} + + +// GParity, (Z)DomainWallVec5D -> require special implementation +#define NO_CURR(Impl) \ +template <> void \ +WilsonKernels::ContractConservedCurrentInternal(const PropagatorField &q_in_1, \ + const PropagatorField &q_in_2, \ + PropagatorField &q_out, \ + DoubledGaugeField &U, \ + Current curr_type, \ + unsigned int mu) \ +{ \ + assert(0); \ +} \ +template <> void \ +WilsonKernels::SeqConservedCurrentInternal(const PropagatorField &q_in, \ + PropagatorField &q_out, \ + DoubledGaugeField &U, \ + Current curr_type, \ + unsigned int mu, \ + Lattice> &ph, \ + unsigned int tmin, \ + unsigned int tmax) \ +{ \ + assert(0); \ +} + +NO_CURR(GparityWilsonImplF); +NO_CURR(GparityWilsonImplD); +NO_CURR(DomainWallVec5dImplF); +NO_CURR(DomainWallVec5dImplD); +NO_CURR(ZDomainWallVec5dImplF); +NO_CURR(ZDomainWallVec5dImplD); + FermOpTemplateInstantiate(WilsonKernels); AdjointFermOpTemplateInstantiate(WilsonKernels); TwoIndexFermOpTemplateInstantiate(WilsonKernels); diff --git a/lib/qcd/action/fermion/WilsonKernels.h b/lib/qcd/action/fermion/WilsonKernels.h index 20ee87f2..34820274 100644 --- a/lib/qcd/action/fermion/WilsonKernels.h +++ b/lib/qcd/action/fermion/WilsonKernels.h @@ -166,6 +166,24 @@ public: void DhopDir(StencilImpl &st, DoubledGaugeField &U,SiteHalfSpinor * buf, int sF, int sU, const FermionField &in, FermionField &out, int dirdisp, int gamma); + ////////////////////////////////////////////////////////////////////////////// + // Utilities for inserting Wilson conserved current. + ////////////////////////////////////////////////////////////////////////////// + void ContractConservedCurrentInternal(const PropagatorField &q_in_1, + const PropagatorField &q_in_2, + PropagatorField &q_out, + DoubledGaugeField &U, + Current curr_type, + unsigned int mu); + void SeqConservedCurrentInternal(const PropagatorField &q_in, + PropagatorField &q_out, + DoubledGaugeField &U, + Current curr_type, + unsigned int mu, + Lattice> &ph, + unsigned int tmin, + unsigned int tmax); + private: // Specialised variants void GenericDhopSite(StencilImpl &st, LebesgueOrder &lo, DoubledGaugeField &U, SiteHalfSpinor * buf, From dc5a6404eadc237985a1b4ffac7b8a51760e6bc4 Mon Sep 17 00:00:00 2001 From: Lanny91 Date: Tue, 25 Apr 2017 22:08:33 +0100 Subject: [PATCH 06/50] Hadrons: modules for testing conserved current contractions and sequential insertion. --- extras/Hadrons/Modules.hpp | 3 + .../Modules/MContraction/WardIdentity.hpp | 151 ++++++++++++++++++ .../Modules/MContraction/WardIdentitySeq.hpp | 117 ++++++++++++++ .../Hadrons/Modules/MSource/SeqConserved.hpp | 129 +++++++++++++++ extras/Hadrons/modules.inc | 3 + 5 files changed, 403 insertions(+) create mode 100644 extras/Hadrons/Modules/MContraction/WardIdentity.hpp create mode 100644 extras/Hadrons/Modules/MContraction/WardIdentitySeq.hpp create mode 100644 extras/Hadrons/Modules/MSource/SeqConserved.hpp diff --git a/extras/Hadrons/Modules.hpp b/extras/Hadrons/Modules.hpp index 05ad1697..67762246 100644 --- a/extras/Hadrons/Modules.hpp +++ b/extras/Hadrons/Modules.hpp @@ -32,6 +32,8 @@ See the full license in the file "LICENSE" in the top level distribution directo #include #include #include +#include +#include #include #include #include @@ -42,6 +44,7 @@ See the full license in the file "LICENSE" in the top level distribution directo #include #include #include +#include #include #include #include diff --git a/extras/Hadrons/Modules/MContraction/WardIdentity.hpp b/extras/Hadrons/Modules/MContraction/WardIdentity.hpp new file mode 100644 index 00000000..39221148 --- /dev/null +++ b/extras/Hadrons/Modules/MContraction/WardIdentity.hpp @@ -0,0 +1,151 @@ +#ifndef Hadrons_WardIdentity_hpp_ +#define Hadrons_WardIdentity_hpp_ + +#include +#include +#include + +BEGIN_HADRONS_NAMESPACE + +/* + Ward Identity contractions + ----------------------------- + + * 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) +*/ + +/****************************************************************************** + * WardIdentity * + ******************************************************************************/ +BEGIN_MODULE_NAMESPACE(MContraction) + +class WardIdentityPar: Serializable +{ +public: + GRID_SERIALIZABLE_CLASS_MEMBERS(WardIdentityPar, + std::string, q, + std::string, q4d, + std::string, action, + double, mass); +}; + +template +class TWardIdentity: public Module +{ +public: + TYPE_ALIASES(FImpl,); +public: + // constructor + TWardIdentity(const std::string name); + // destructor + virtual ~TWardIdentity(void) = default; + // dependency relation + virtual std::vector getInput(void); + virtual std::vector getOutput(void); + // setup + virtual void setup(void); + // execution + virtual void execute(void); +private: + unsigned int Ls_; +}; + +MODULE_REGISTER_NS(WardIdentity, TWardIdentity, MContraction); + +/****************************************************************************** + * TWardIdentity implementation * + ******************************************************************************/ +// constructor ///////////////////////////////////////////////////////////////// +template +TWardIdentity::TWardIdentity(const std::string name) +: Module(name) +{} + +// dependencies/products /////////////////////////////////////////////////////// +template +std::vector TWardIdentity::getInput(void) +{ + std::vector in = {par().q, par().q4d, par().action}; + + return in; +} + +template +std::vector TWardIdentity::getOutput(void) +{ + std::vector out = {getName()}; + + return out; +} + +// setup /////////////////////////////////////////////////////////////////////// +template +void TWardIdentity::setup(void) +{ + Ls_ = env().getObjectLs(par().q); +} + +// execution /////////////////////////////////////////////////////////////////// +template +void TWardIdentity::execute(void) +{ + LOG(Message) << "Performing Ward Identity checks for quark '" << par().q + << "'." << std::endl; + + PropagatorField psi(env().getGrid()), tmp(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); + for (unsigned int mu = 0; mu < Nd; ++mu) + { + act.ContractConservedCurrent(q, q, tmp, Current::Vector, mu); + Vmu[mu] = trace(g5*tmp); + act.ContractConservedCurrent(q, q, tmp, Current::Axial, mu); + Amu[mu] = trace(g5*tmp); + } + + // Get PJ5q for 5D (zero for 4D). + if (Ls_ > 1) + { + ExtractSlice(psi, q, 0, Ls_/2 - 1); + psi = 0.5 * (psi + g5*psi); + ExtractSlice(tmp, q, 0, Ls_/2); + 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; +} + +END_MODULE_NAMESPACE + +END_HADRONS_NAMESPACE + +#endif // Hadrons_WardIdentity_hpp_ diff --git a/extras/Hadrons/Modules/MContraction/WardIdentitySeq.hpp b/extras/Hadrons/Modules/MContraction/WardIdentitySeq.hpp new file mode 100644 index 00000000..3e72c11e --- /dev/null +++ b/extras/Hadrons/Modules/MContraction/WardIdentitySeq.hpp @@ -0,0 +1,117 @@ +#ifndef Hadrons_WardIdentitySeq_hpp_ +#define Hadrons_WardIdentitySeq_hpp_ + +#include +#include +#include + +BEGIN_HADRONS_NAMESPACE + +/* + Ward Identity contractions using sequential propagators. + ----------------------------- + + * options: + - q_x: propagator, mu = x current insertion (string). + - q_y: propagator, mu = y current insertion (string). + - q_z: propagator, mu = z current insertion (string). + - q_t: propagator, mu = t current insertion (string). +*/ + +/****************************************************************************** + * WardIdentitySeq * + ******************************************************************************/ +BEGIN_MODULE_NAMESPACE(MContraction) + +class WardIdentitySeqPar: Serializable +{ +public: + GRID_SERIALIZABLE_CLASS_MEMBERS(WardIdentitySeqPar, + std::string, q_x, + std::string, q_y, + std::string, q_z, + std::string, q_t); +}; + +template +class TWardIdentitySeq: public Module +{ +public: + TYPE_ALIASES(FImpl,); +public: + // constructor + TWardIdentitySeq(const std::string name); + // destructor + virtual ~TWardIdentitySeq(void) = default; + // dependency relation + virtual std::vector getInput(void); + virtual std::vector getOutput(void); + // setup + virtual void setup(void); + // execution + virtual void execute(void); +}; + +MODULE_REGISTER_NS(WardIdentitySeq, TWardIdentitySeq, MContraction); + +/****************************************************************************** + * TWardIdentitySeq implementation * + ******************************************************************************/ +// constructor ///////////////////////////////////////////////////////////////// +template +TWardIdentitySeq::TWardIdentitySeq(const std::string name) +: Module(name) +{} + +// dependencies/products /////////////////////////////////////////////////////// +template +std::vector TWardIdentitySeq::getInput(void) +{ + std::vector in = {par().q_x, par().q_y, par().q_z, par().q_t}; + + return in; +} + +template +std::vector TWardIdentitySeq::getOutput(void) +{ + std::vector out = {getName()}; + + return out; +} + +// setup /////////////////////////////////////////////////////////////////////// +template +void TWardIdentitySeq::setup(void) +{ + +} + +// execution /////////////////////////////////////////////////////////////////// +template +void TWardIdentitySeq::execute(void) +{ + LatticeComplex vector_WI(env().getGrid()), c(env().getGrid()); + PropagatorField q_x = *env().template getObject(par().q_x); + PropagatorField q_y = *env().template getObject(par().q_y); + PropagatorField q_z = *env().template getObject(par().q_z); + PropagatorField q_t = *env().template getObject(par().q_t); + PropagatorField *q[Nd] = {&q_x, &q_y, &q_z, &q_t}; + Gamma g5(Gamma::Algebra::Gamma5); + + // Check D_mu V_mu = 0 + for (unsigned int mu = 0; mu < Nd; ++mu) + { + c = trace(g5*(*q[mu])); + vector_WI += c - Cshift(c, mu, -1); + } + + LOG(Message) << "Ward Identity checks for sequential vector current " + << "insertion = " << norm2(vector_WI) << std::endl; +} + +END_MODULE_NAMESPACE + +END_HADRONS_NAMESPACE + +#endif // Hadrons_WardIdentitySeq_hpp_ diff --git a/extras/Hadrons/Modules/MSource/SeqConserved.hpp b/extras/Hadrons/Modules/MSource/SeqConserved.hpp new file mode 100644 index 00000000..7d4974f4 --- /dev/null +++ b/extras/Hadrons/Modules/MSource/SeqConserved.hpp @@ -0,0 +1,129 @@ +#ifndef Hadrons_SeqConserved_hpp_ +#define Hadrons_SeqConserved_hpp_ + +#include +#include +#include + +BEGIN_HADRONS_NAMESPACE + +/* + + Sequential source + ----------------------------- + * src_x = q_x * theta(x_3 - tA) * theta(tB - x_3) * J_mu * exp(i x.mom) + + * options: + - q: input propagator (string) + - action: fermion action used for propagator q (string) + - tA: begin timeslice (integer) + - tB: end timesilce (integer) + - curr_type: type of conserved current to insert (Current) + - mu: Lorentz index of current to insert (integer) + - mom: momentum insertion, space-separated float sequence (e.g ".1 .2 1. 0.") + + */ + +/****************************************************************************** + * SeqConserved * + ******************************************************************************/ +BEGIN_MODULE_NAMESPACE(MSource) + +class SeqConservedPar: Serializable +{ +public: + GRID_SERIALIZABLE_CLASS_MEMBERS(SeqConservedPar, + std::string, q, + std::string, action, + unsigned int, tA, + unsigned int, tB, + Current, curr_type, + unsigned int, mu, + std::string, mom); +}; + +template +class TSeqConserved: public Module +{ +public: + TYPE_ALIASES(FImpl,); +public: + // constructor + TSeqConserved(const std::string name); + // destructor + virtual ~TSeqConserved(void) = default; + // dependency relation + virtual std::vector getInput(void); + virtual std::vector getOutput(void); + // setup + virtual void setup(void); + // execution + virtual void execute(void); +}; + +MODULE_REGISTER_NS(SeqConserved, TSeqConserved, MSource); + +/****************************************************************************** + * TSeqConserved implementation * + ******************************************************************************/ +// constructor ///////////////////////////////////////////////////////////////// +template +TSeqConserved::TSeqConserved(const std::string name) +: Module(name) +{} + +// dependencies/products /////////////////////////////////////////////////////// +template +std::vector TSeqConserved::getInput(void) +{ + std::vector in; + + return in; +} + +template +std::vector TSeqConserved::getOutput(void) +{ + std::vector out = {getName()}; + + return out; +} + +// setup /////////////////////////////////////////////////////////////////////// +template +void TSeqConserved::setup(void) +{ + +} + +// execution /////////////////////////////////////////////////////////////////// +template +void TSeqConserved::execute(void) +{ + if (par().tA == par().tB) + { + LOG(Message) << "Generating sequential source with conserved " + << par().curr_type << " current insertion (mu = " + << par().mu << ") at " << "t = " << par().tA << std::endl; + } + else + { + LOG(Message) << "Generating sequential source with conserved " + << par().curr_type << " current insertion (mu = " + << par().mu << ") for " << par().tA << " <= t <= " + << par().tB << std::endl; + } + PropagatorField &src = *env().template createLattice(getName()); + PropagatorField &q = *env().template getObject(par().q); + FMat &mat = *(env().template getObject(par().action)); + + std::vector mom = strToVec(par().mom); + mat.SeqConservedCurrent(q, src, par().curr_type, par().mu, + mom, par().tA, par().tB); +} + +END_MODULE_NAMESPACE + +END_HADRONS_NAMESPACE + +#endif // Hadrons_SeqConserved_hpp_ diff --git a/extras/Hadrons/modules.inc b/extras/Hadrons/modules.inc index af291631..32655c3b 100644 --- a/extras/Hadrons/modules.inc +++ b/extras/Hadrons/modules.inc @@ -13,6 +13,8 @@ modules_hpp =\ Modules/MContraction/DiscLoop.hpp \ Modules/MContraction/Gamma3pt.hpp \ Modules/MContraction/Meson.hpp \ + Modules/MContraction/WardIdentity.hpp \ + Modules/MContraction/WardIdentitySeq.hpp \ Modules/MContraction/WeakHamiltonian.hpp \ Modules/MContraction/WeakHamiltonianEye.hpp \ Modules/MContraction/WeakHamiltonianNonEye.hpp \ @@ -23,6 +25,7 @@ modules_hpp =\ Modules/MLoop/NoiseLoop.hpp \ Modules/MSolver/RBPrecCG.hpp \ Modules/MSource/Point.hpp \ + Modules/MSource/SeqConserved.hpp \ Modules/MSource/SeqGamma.hpp \ Modules/MSource/Wall.hpp \ Modules/MSource/Z2.hpp \ From 159770e21ba2515d145bb331305593474ce33b01 Mon Sep 17 00:00:00 2001 From: Lanny91 Date: Wed, 26 Apr 2017 09:32:57 +0100 Subject: [PATCH 07/50] Legal Banners added --- .../Modules/MContraction/WardIdentity.hpp | 28 +++++++++++++++++++ .../Modules/MContraction/WardIdentitySeq.hpp | 28 +++++++++++++++++++ .../Hadrons/Modules/MSource/SeqConserved.hpp | 28 +++++++++++++++++++ 3 files changed, 84 insertions(+) diff --git a/extras/Hadrons/Modules/MContraction/WardIdentity.hpp b/extras/Hadrons/Modules/MContraction/WardIdentity.hpp index 39221148..355126da 100644 --- a/extras/Hadrons/Modules/MContraction/WardIdentity.hpp +++ b/extras/Hadrons/Modules/MContraction/WardIdentity.hpp @@ -1,3 +1,31 @@ +/************************************************************************************* + +Grid physics library, www.github.com/paboyle/Grid + +Source file: extras/Hadrons/Modules/MContraction/WardIdentity.hpp + +Copyright (C) 2017 + +Author: Andrew Lawson + +This program is free software; you can redistribute it and/or modify +it under the terms of the GNU General Public License as published by +the Free Software Foundation; either version 2 of the License, or +(at your option) any later version. + +This program is distributed in the hope that it will be useful, +but WITHOUT ANY WARRANTY; without even the implied warranty of +MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +GNU General Public License for more details. + +You should have received a copy of the GNU General Public License along +with this program; if not, write to the Free Software Foundation, Inc., +51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. + +See the full license in the file "LICENSE" in the top level distribution directory +*************************************************************************************/ +/* END LEGAL */ + #ifndef Hadrons_WardIdentity_hpp_ #define Hadrons_WardIdentity_hpp_ diff --git a/extras/Hadrons/Modules/MContraction/WardIdentitySeq.hpp b/extras/Hadrons/Modules/MContraction/WardIdentitySeq.hpp index 3e72c11e..31409925 100644 --- a/extras/Hadrons/Modules/MContraction/WardIdentitySeq.hpp +++ b/extras/Hadrons/Modules/MContraction/WardIdentitySeq.hpp @@ -1,3 +1,31 @@ +/************************************************************************************* + +Grid physics library, www.github.com/paboyle/Grid + +Source file: extras/Hadrons/Modules/MContraction/WardIdentitySeq.hpp + +Copyright (C) 2017 + +Author: Andrew Lawson + +This program is free software; you can redistribute it and/or modify +it under the terms of the GNU General Public License as published by +the Free Software Foundation; either version 2 of the License, or +(at your option) any later version. + +This program is distributed in the hope that it will be useful, +but WITHOUT ANY WARRANTY; without even the implied warranty of +MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +GNU General Public License for more details. + +You should have received a copy of the GNU General Public License along +with this program; if not, write to the Free Software Foundation, Inc., +51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. + +See the full license in the file "LICENSE" in the top level distribution directory +*************************************************************************************/ +/* END LEGAL */ + #ifndef Hadrons_WardIdentitySeq_hpp_ #define Hadrons_WardIdentitySeq_hpp_ diff --git a/extras/Hadrons/Modules/MSource/SeqConserved.hpp b/extras/Hadrons/Modules/MSource/SeqConserved.hpp index 7d4974f4..ccfb68f4 100644 --- a/extras/Hadrons/Modules/MSource/SeqConserved.hpp +++ b/extras/Hadrons/Modules/MSource/SeqConserved.hpp @@ -1,3 +1,31 @@ +/************************************************************************************* + +Grid physics library, www.github.com/paboyle/Grid + +Source file: extras/Hadrons/Modules/MContraction/SeqConserved.hpp + +Copyright (C) 2017 + +Author: Andrew Lawson + +This program is free software; you can redistribute it and/or modify +it under the terms of the GNU General Public License as published by +the Free Software Foundation; either version 2 of the License, or +(at your option) any later version. + +This program is distributed in the hope that it will be useful, +but WITHOUT ANY WARRANTY; without even the implied warranty of +MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +GNU General Public License for more details. + +You should have received a copy of the GNU General Public License along +with this program; if not, write to the Free Software Foundation, Inc., +51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. + +See the full license in the file "LICENSE" in the top level distribution directory +*************************************************************************************/ +/* END LEGAL */ + #ifndef Hadrons_SeqConserved_hpp_ #define Hadrons_SeqConserved_hpp_ From a39daecb62a0bf8a128f4650d311badaa1659fda Mon Sep 17 00:00:00 2001 From: Lanny91 Date: Wed, 26 Apr 2017 12:39:07 +0100 Subject: [PATCH 08/50] Removed make_5D const declaration to avoid compilation error --- extras/Hadrons/Modules/Quark.hpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/extras/Hadrons/Modules/Quark.hpp b/extras/Hadrons/Modules/Quark.hpp index c08e0192..fff34edf 100644 --- a/extras/Hadrons/Modules/Quark.hpp +++ b/extras/Hadrons/Modules/Quark.hpp @@ -48,7 +48,7 @@ inline void make_4D(Lattice &in_5d, Lattice &out_4d, int Ls) } template -inline void make_5D(const Lattice &in_4d, Lattice &out_5d, int Ls) +inline void make_5D(Lattice &in_4d, Lattice &out_5d, int Ls) { out_5d = zero; InsertSlice(in_4d, out_5d, 0, 0); From 6299dd35f57b03131447d70f7b8e7f002dc4cdf9 Mon Sep 17 00:00:00 2001 From: Lanny91 Date: Wed, 26 Apr 2017 12:41:39 +0100 Subject: [PATCH 09/50] Hadrons: Added test of conserved current code. Tests Ward identities for conserved vector and partially conserved axial currents. --- .../hadrons/Test_hadrons_conserved_current.cc | 127 ++++++++++++++++++ 1 file changed, 127 insertions(+) create mode 100644 tests/hadrons/Test_hadrons_conserved_current.cc diff --git a/tests/hadrons/Test_hadrons_conserved_current.cc b/tests/hadrons/Test_hadrons_conserved_current.cc new file mode 100644 index 00000000..df774ac0 --- /dev/null +++ b/tests/hadrons/Test_hadrons_conserved_current.cc @@ -0,0 +1,127 @@ +/******************************************************************************* + Grid physics library, www.github.com/paboyle/Grid + + Source file: tests/hadrons/Test_hadrons_conserved_current.cc + + Copyright (C) 2017 + + Author: Andrew Lawson + + This program is free software; you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation; either version 2 of the License, or + (at your option) any later version. + + This program is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License along + with this program; if not, write to the Free Software Foundation, Inc., + 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. + + See the full license in the file "LICENSE" in the top level distribution + directory. + *******************************************************************************/ + +#include "Test_hadrons.hpp" + +using namespace Grid; +using namespace Hadrons; + +int main(int argc, char *argv[]) +{ + // initialization ////////////////////////////////////////////////////////// + Grid_init(&argc, &argv); + HadronsLogError.Active(GridLogError.isActive()); + HadronsLogWarning.Active(GridLogWarning.isActive()); + HadronsLogMessage.Active(GridLogMessage.isActive()); + HadronsLogIterative.Active(GridLogIterative.isActive()); + HadronsLogDebug.Active(GridLogDebug.isActive()); + LOG(Message) << "Grid initialized" << std::endl; + + // run setup /////////////////////////////////////////////////////////////// + Application application; + unsigned int nt = GridDefaultLatt()[Tp]; + double mass = 0.04; + + // global parameters + Application::GlobalPar globalPar; + globalPar.trajCounter.start = 1500; + globalPar.trajCounter.end = 1520; + globalPar.trajCounter.step = 20; + globalPar.seed = "1 2 3 4"; + globalPar.genetic.maxGen = 1000; + globalPar.genetic.maxCstGen = 200; + globalPar.genetic.popSize = 20; + globalPar.genetic.mutationRate = .1; + application.setPar(globalPar); + + // gauge field + application.createModule("gauge"); + + // action + std::string actionName = "DWF"; + MAction::DWF::Par actionPar; + actionPar.gauge = "gauge"; + actionPar.Ls = 12; + actionPar.M5 = 1.8; + actionPar.mass = mass; + application.createModule(actionName, actionPar); + + // solver + std::string solverName = "CG"; + MSolver::RBPrecCG::Par solverPar; + solverPar.action = actionName; + solverPar.residual = 1.0e-8; + application.createModule(solverName, + solverPar); + + // Conserved current sink contractions: use a single point propagator. + std::string pointProp = "q_0"; + std::string pos = "0 0 0 0"; + std::string modName = "Ward Identity Test"; + MAKE_POINT_PROP(pos, pointProp, solverName); + if (!(Environment::getInstance().hasModule(modName))) + { + MContraction::WardIdentity::Par wiPar; + wiPar.q = pointProp + "_5d"; + wiPar.q4d = pointProp; + wiPar.action = actionName; + wiPar.mass = mass; + application.createModule(modName, wiPar); + } + + // Conserved current contractions with sequential insertion of vector + // current. + std::string q_x = "q_x"; + std::string q_y = "q_y"; + std::string q_z = "q_z"; + std::string q_t = "q_t"; + std::string mom = ZERO_MOM; + modName = "Sequential Ward Identity Test"; + MAKE_SEQUENTIAL_PROP(nt/2, pointProp, mom, q_x, solverName); + MAKE_SEQUENTIAL_PROP(nt/2, pointProp, mom, q_y, solverName); + MAKE_SEQUENTIAL_PROP(nt/2, pointProp, mom, q_z, solverName); + MAKE_SEQUENTIAL_PROP(nt/2, pointProp, mom, q_t, solverName); + if (!(Environment::getInstance().hasModule(modName))) + { + MContraction::WardIdentitySeq::Par wiPar; + wiPar.q_x = q_x; + wiPar.q_y = q_y; + wiPar.q_z = q_z; + wiPar.q_t = q_t; + application.createModule(modName, wiPar); + } + + // execution + application.saveParameterFile("ConservedCurrentTest.xml"); + application.run(); + + // epilogue + LOG(Message) << "Grid is finalizing now" << std::endl; + Grid_finalize(); + + return EXIT_SUCCESS; +} \ No newline at end of file From d2003f24f49b68f49fd0cb5764ef32f6ef4cd498 Mon Sep 17 00:00:00 2001 From: Lanny91 Date: Wed, 26 Apr 2017 17:25:28 +0100 Subject: [PATCH 10/50] Corrected incorrect usage of ExtractSlice for conserved current code. --- extras/Hadrons/Modules/MContraction/WardIdentity.hpp | 4 ++-- lib/qcd/action/fermion/WilsonFermion5D.cc | 4 ++-- 2 files changed, 4 insertions(+), 4 deletions(-) diff --git a/extras/Hadrons/Modules/MContraction/WardIdentity.hpp b/extras/Hadrons/Modules/MContraction/WardIdentity.hpp index 355126da..41d8c6d1 100644 --- a/extras/Hadrons/Modules/MContraction/WardIdentity.hpp +++ b/extras/Hadrons/Modules/MContraction/WardIdentity.hpp @@ -150,9 +150,9 @@ void TWardIdentity::execute(void) // Get PJ5q for 5D (zero for 4D). if (Ls_ > 1) { - ExtractSlice(psi, q, 0, Ls_/2 - 1); + ExtractSlice(psi, q, Ls_/2 - 1, 0); psi = 0.5 * (psi + g5*psi); - ExtractSlice(tmp, q, 0, Ls_/2); + ExtractSlice(tmp, q, Ls_/2, 0); psi += 0.5 * (tmp - g5*tmp); PJ5q = trace(adj(psi)*psi); } diff --git a/lib/qcd/action/fermion/WilsonFermion5D.cc b/lib/qcd/action/fermion/WilsonFermion5D.cc index d0d3d055..99ff0dc1 100644 --- a/lib/qcd/action/fermion/WilsonFermion5D.cc +++ b/lib/qcd/action/fermion/WilsonFermion5D.cc @@ -703,8 +703,8 @@ void WilsonFermion5D::ContractConservedCurrent(PropagatorField &q_in_1, q_out = zero; for (int s = 0; s < Ls; ++s) { - ExtractSlice(q1_s, q_in_1, 0, s); - ExtractSlice(q2_s, q_in_2, 0, Ls - s - 1); + ExtractSlice(q1_s, q_in_1, s, 0); + ExtractSlice(q2_s, q_in_2, Ls - s - 1, 0); Kernels::ContractConservedCurrentInternal(q1_s, q2_s, tmp, Umu, curr_type, mu); // Axial current sign From a6ccbbe1080f8c7d2bd193e08ad38b97d0f5af1d Mon Sep 17 00:00:00 2001 From: Lanny91 Date: Fri, 28 Apr 2017 10:43:47 +0100 Subject: [PATCH 11/50] Conserved current sequential source now registered properly and fixed module inputs. --- extras/Hadrons/Modules/MSource/SeqConserved.hpp | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/extras/Hadrons/Modules/MSource/SeqConserved.hpp b/extras/Hadrons/Modules/MSource/SeqConserved.hpp index ccfb68f4..6e5fb197 100644 --- a/extras/Hadrons/Modules/MSource/SeqConserved.hpp +++ b/extras/Hadrons/Modules/MSource/SeqConserved.hpp @@ -104,7 +104,7 @@ TSeqConserved::TSeqConserved(const std::string name) template std::vector TSeqConserved::getInput(void) { - std::vector in; + std::vector in = {par().q, par().action}; return in; } @@ -121,7 +121,8 @@ std::vector TSeqConserved::getOutput(void) template void TSeqConserved::setup(void) { - + auto Ls_ = env().getObjectLs(par().action); + env().template registerLattice(getName(), Ls_); } // execution /////////////////////////////////////////////////////////////////// From f302eea91e00e1b5e47190253708f9fd1ac36ff7 Mon Sep 17 00:00:00 2001 From: Lanny91 Date: Fri, 28 Apr 2017 15:27:49 +0100 Subject: [PATCH 12/50] SitePropagator redefined to be a scalar object in TYPE_ALIASES. --- extras/Hadrons/Global.hpp | 14 +++++++------- extras/Hadrons/Modules/MSource/Point.hpp | 2 +- 2 files changed, 8 insertions(+), 8 deletions(-) diff --git a/extras/Hadrons/Global.hpp b/extras/Hadrons/Global.hpp index 3e11ddf8..ebf93283 100644 --- a/extras/Hadrons/Global.hpp +++ b/extras/Hadrons/Global.hpp @@ -60,13 +60,13 @@ BEGIN_HADRONS_NAMESPACE // type aliases #define TYPE_ALIASES(FImpl, suffix)\ -typedef FermionOperator FMat##suffix; \ -typedef typename FImpl::FermionField FermionField##suffix; \ -typedef typename FImpl::PropagatorField PropagatorField##suffix; \ -typedef typename FImpl::SitePropagator SitePropagator##suffix; \ -typedef typename FImpl::DoubledGaugeField DoubledGaugeField##suffix;\ -typedef std::function SolverFn##suffix; +typedef FermionOperator FMat##suffix; \ +typedef typename FImpl::FermionField FermionField##suffix; \ +typedef typename FImpl::PropagatorField PropagatorField##suffix; \ +typedef typename FImpl::SitePropagator::scalar_object SitePropagator##suffix; \ +typedef typename FImpl::DoubledGaugeField DoubledGaugeField##suffix;\ +typedef std::function SolverFn##suffix; // logger class HadronsLogger: public Logger diff --git a/extras/Hadrons/Modules/MSource/Point.hpp b/extras/Hadrons/Modules/MSource/Point.hpp index a0ecbc2a..36e1cc5b 100644 --- a/extras/Hadrons/Modules/MSource/Point.hpp +++ b/extras/Hadrons/Modules/MSource/Point.hpp @@ -118,7 +118,7 @@ template void TPoint::execute(void) { std::vector position = strToVec(par().position); - typename SitePropagator::scalar_object id; + SitePropagator id; LOG(Message) << "Creating point source at position [" << par().position << "]" << std::endl; From b9356d38662144acfa5b3f1d5184123b9769598d Mon Sep 17 00:00:00 2001 From: Lanny91 Date: Fri, 28 Apr 2017 16:46:40 +0100 Subject: [PATCH 13/50] Added more complete test of sequential insertion of conserved current. --- extras/Hadrons/Modules.hpp | 1 + .../Modules/MUtilities/TestSeqConserved.hpp | 166 ++++++++++++++++++ extras/Hadrons/modules.inc | 1 + 3 files changed, 168 insertions(+) create mode 100644 extras/Hadrons/Modules/MUtilities/TestSeqConserved.hpp diff --git a/extras/Hadrons/Modules.hpp b/extras/Hadrons/Modules.hpp index 67762246..0286333c 100644 --- a/extras/Hadrons/Modules.hpp +++ b/extras/Hadrons/Modules.hpp @@ -48,4 +48,5 @@ See the full license in the file "LICENSE" in the top level distribution directo #include #include #include +#include #include diff --git a/extras/Hadrons/Modules/MUtilities/TestSeqConserved.hpp b/extras/Hadrons/Modules/MUtilities/TestSeqConserved.hpp new file mode 100644 index 00000000..0730b8ed --- /dev/null +++ b/extras/Hadrons/Modules/MUtilities/TestSeqConserved.hpp @@ -0,0 +1,166 @@ +/************************************************************************************* + +Grid physics library, www.github.com/paboyle/Grid + +Source file: extras/Hadrons/Modules/MUtilities/TestSeqConserved.hpp + +Copyright (C) 2017 + +Author: Andrew Lawson + +This program is free software; you can redistribute it and/or modify +it under the terms of the GNU General Public License as published by +the Free Software Foundation; either version 2 of the License, or +(at your option) any later version. + +This program is distributed in the hope that it will be useful, +but WITHOUT ANY WARRANTY; without even the implied warranty of +MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +GNU General Public License for more details. + +You should have received a copy of the GNU General Public License along +with this program; if not, write to the Free Software Foundation, Inc., +51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. + +See the full license in the file "LICENSE" in the top level distribution directory +*************************************************************************************/ +/* END LEGAL */ + +#ifndef Hadrons_TestSeqConserved_hpp_ +#define Hadrons_TestSeqConserved_hpp_ + +#include +#include +#include + +BEGIN_HADRONS_NAMESPACE + +/* + Ward Identity contractions using sequential propagators. + ----------------------------- + + * options: + - q: point source propagator, 5D if available (string) + - q4d: 4D point source propagator, duplicate of q if q is 4D (string) + - qSeq: result of sequential insertion of conserved current using q (string) + - action: action used for computation of q (string) + - origin: string giving point source origin of q (string) + - t_J: time at which sequential current is inserted (int) + - mu: Lorentz index of current inserted (int) + - curr: current type, e.g. vector/axial (Current) +*/ + +/****************************************************************************** + * TestSeqConserved * + ******************************************************************************/ +BEGIN_MODULE_NAMESPACE(MUtilities) + +class TestSeqConservedPar: Serializable +{ +public: + GRID_SERIALIZABLE_CLASS_MEMBERS(TestSeqConservedPar, + std::string, q, + std::string, q4d, + std::string, qSeq, + std::string, action, + std::string, origin, + unsigned int, t_J, + unsigned int, mu, + Current, curr); +}; + +template +class TTestSeqConserved: public Module +{ +public: + TYPE_ALIASES(FImpl,); +public: + // constructor + TTestSeqConserved(const std::string name); + // destructor + virtual ~TTestSeqConserved(void) = default; + // dependency relation + virtual std::vector getInput(void); + virtual std::vector getOutput(void); + // setup + virtual void setup(void); + // execution + virtual void execute(void); +}; + +MODULE_REGISTER_NS(TestSeqConserved, TTestSeqConserved, MUtilities); + +/****************************************************************************** + * TTestSeqConserved implementation * + ******************************************************************************/ +// constructor ///////////////////////////////////////////////////////////////// +template +TTestSeqConserved::TTestSeqConserved(const std::string name) +: Module(name) +{} + +// dependencies/products /////////////////////////////////////////////////////// +template +std::vector TTestSeqConserved::getInput(void) +{ + std::vector in = {par().q, par().q4d, + par().qSeq, par().action}; + + return in; +} + +template +std::vector TTestSeqConserved::getOutput(void) +{ + std::vector out = {getName()}; + + return out; +} + +// setup /////////////////////////////////////////////////////////////////////// +template +void TTestSeqConserved::setup(void) +{ + +} + +// execution /////////////////////////////////////////////////////////////////// +template +void TTestSeqConserved::execute(void) +{ + PropagatorField tmp(env().getGrid()); + PropagatorField &q = *env().template getObject(par().q); + PropagatorField &q4d = *env().template getObject(par().q4d); + PropagatorField &qSeq = *env().template getObject(par().qSeq); + FMat &act = *(env().template getObject(par().action)); + Gamma g5(Gamma::Algebra::Gamma5); + SitePropagator qSite; + LatticeComplex c(env().getGrid()); + Complex seq_res, check_res; + std::vector check_buf; + + // Check sequential insertion of current gives same result as conserved + // current sink upon contraction. Assume q uses a point source. + std::vector siteCoord; + siteCoord = strToVec(par().origin); + peekSite(qSite, q, siteCoord); + seq_res = trace(g5*qSite); + + act.ContractConservedCurrent(q, q, tmp, par().curr, par().mu); + c = trace(tmp); + sliceSum(c, check_buf, Tp); + check_res = TensorRemove(check_buf[par().t_J]); + + // Check difference = 0 + check_res -= seq_res; + + LOG(Message) << "Consistency check for sequential conserved " + << par().curr << " current insertion = " << abs(check_res) + << std::endl; +} + +END_MODULE_NAMESPACE + +END_HADRONS_NAMESPACE + +#endif // Hadrons_TestSeqConserved_hpp_ diff --git a/extras/Hadrons/modules.inc b/extras/Hadrons/modules.inc index 32655c3b..4ab51ce0 100644 --- a/extras/Hadrons/modules.inc +++ b/extras/Hadrons/modules.inc @@ -29,5 +29,6 @@ modules_hpp =\ Modules/MSource/SeqGamma.hpp \ Modules/MSource/Wall.hpp \ Modules/MSource/Z2.hpp \ + Modules/MUtilities/TestSeqConserved.hpp \ Modules/Quark.hpp From db14fb30df901c5e6bb445a2c59a55bf28ac75d5 Mon Sep 17 00:00:00 2001 From: Lanny91 Date: Fri, 28 Apr 2017 16:48:00 +0100 Subject: [PATCH 14/50] Hadrons: overhaul of conserved current test --- tests/hadrons/Test_hadrons.hpp | 125 +++++++++++++++++- .../hadrons/Test_hadrons_conserved_current.cc | 54 ++++---- 2 files changed, 145 insertions(+), 34 deletions(-) diff --git a/tests/hadrons/Test_hadrons.hpp b/tests/hadrons/Test_hadrons.hpp index 26d02a5c..c4dcedaf 100644 --- a/tests/hadrons/Test_hadrons.hpp +++ b/tests/hadrons/Test_hadrons.hpp @@ -40,6 +40,7 @@ using namespace Hadrons; #define LABEL_3PT(s, t1, t2) ADD_INDEX(INIT_INDEX(s, t1), t2) #define LABEL_4PT(s, t1, t2, t3) ADD_INDEX(ADD_INDEX(INIT_INDEX(s, t1), t2), t3) #define LABEL_4PT_NOISE(s, t1, t2, t3, nn) ADD_INDEX(ADD_INDEX(ADD_INDEX(INIT_INDEX(s, t1), t2), t3), nn) +#define LABEL_5D(s) s + "_5d"; // Wall source/sink macros #define NAME_3MOM_WALL_SOURCE(t, mom) ("wall_" + std::to_string(t) + "_" + mom) @@ -123,6 +124,44 @@ inline void makeSequentialSource(Application &application, std::string srcName, } } +/******************************************************************************* + * Name: makeConservedSequentialSource + * Purpose: Construct sequential source with conserved current insertion and + * add to application module. + * Parameters: application - main application that stores modules. + * srcName - name of source module to create. + * qSrc - Input quark for sequential inversion. + * actionName - action corresponding to quark. + * tS - sequential source timeslice. + * curr - conserved current type to insert. + * mu - Lorentz index of current to insert. + * mom - momentum insertion (default is zero). + * Returns: None. + ******************************************************************************/ +inline void makeConservedSequentialSource(Application &application, + std::string &srcName, + std::string &qSrc, + std::string &actionName, + unsigned int tS, + Current curr, + unsigned int mu, + std::string mom = ZERO_MOM) +{ + // If the source already exists, don't make the module again. + if (!(Environment::getInstance().hasModule(srcName))) + { + MSource::SeqConserved::Par seqPar; + seqPar.q = qSrc; + seqPar.action = actionName; + seqPar.tA = tS; + seqPar.tB = tS; + seqPar.curr_type = curr; + seqPar.mu = mu; + seqPar.mom = mom; + application.createModule(srcName, seqPar); + } +} + /******************************************************************************* * Name: makeWallSource * Purpose: Construct wall source and add to application module. @@ -132,7 +171,7 @@ inline void makeSequentialSource(Application &application, std::string srcName, * mom - momentum insertion (default is zero). * Returns: None. ******************************************************************************/ -inline void makeWallSource(Application &application, std::string srcName, +inline void makeWallSource(Application &application, std::string &srcName, unsigned int tW, std::string mom = ZERO_MOM) { // If the source already exists, don't make the module again. @@ -154,8 +193,8 @@ inline void makeWallSource(Application &application, std::string srcName, * mom - momentum insertion (default is zero). * Returns: None. ******************************************************************************/ -inline void makeWallSink(Application &application, std::string propName, - std::string wallName, std::string mom = ZERO_MOM) +inline void makeWallSink(Application &application, std::string &propName, + std::string &wallName, std::string mom = ZERO_MOM) { // If the propagator has already been smeared, don't smear it again. // Temporarily removed, strategy for sink smearing likely to change. @@ -365,4 +404,82 @@ inline void discLoopContraction(Application &application, discPar.gamma = gamma; application.createModule(modName, discPar); } - } +} + +/******************************************************************************* + * Name: makeWITest + * Purpose: Create module to test Ward Identities for conserved current + * contractions and add to application module. + * Parameters: application - main application that stores modules. + * modName - name of module to create. + * propName - 4D quark propagator. + * actionName - action used to compute quark propagator. + * mass - mass of quark. + * Ls - length of 5th dimension (default = 1). + * Returns: None. + ******************************************************************************/ +inline void makeWITest(Application &application, std::string &modName, + std::string &propName, std::string &actionName, + double mass, unsigned int Ls = 1) +{ + if (!(Environment::getInstance().hasModule(modName))) + { + MContraction::WardIdentity::Par wiPar; + if (Ls > 1) + { + wiPar.q = LABEL_5D(propName); + } + else + { + wiPar.q = propName; + } + wiPar.q4d = propName; + wiPar.action = actionName; + wiPar.mass = mass; + application.createModule(modName, wiPar); + } +} + +/******************************************************************************* + * Name: makeSeqTest + * Purpose: Create module to test sequential insertion of conserved current + * and add to application module. + * Parameters: application - main application that stores modules. + * modName - name of module to create. + * propName - 4D quark propagator. + * seqProp - 4D quark propagator with sequential insertion of + * conserved current. + * actionName - action used to compute quark propagators. + * t_J - time at which sequential current is inserted. + * mu - Lorentz index of sequential current. + * curr - type of conserved current inserted. + * Ls - length of 5th dimension (default = 1). + * Returns: None. + ******************************************************************************/ +inline void makeSeqTest(Application &application, std::string &modName, + std::string &propName, std::string &seqName, + std::string &actionName, std::string &origin, + unsigned int t_J, unsigned int mu, Current curr, + unsigned int Ls = 1) +{ + if (!(Environment::getInstance().hasModule(modName))) + { + MUtilities::TestSeqConserved::Par seqPar; + if (Ls > 1) + { + seqPar.q = LABEL_5D(propName); + } + else + { + seqPar.q = propName; + } + seqPar.q4d = propName; + seqPar.qSeq = seqName; + seqPar.action = actionName; + seqPar.origin = origin; + seqPar.t_J = t_J; + seqPar.mu = mu; + seqPar.curr = curr; + application.createModule(modName, seqPar); + } +} diff --git a/tests/hadrons/Test_hadrons_conserved_current.cc b/tests/hadrons/Test_hadrons_conserved_current.cc index df774ac0..a11a3530 100644 --- a/tests/hadrons/Test_hadrons_conserved_current.cc +++ b/tests/hadrons/Test_hadrons_conserved_current.cc @@ -45,6 +45,7 @@ int main(int argc, char *argv[]) Application application; unsigned int nt = GridDefaultLatt()[Tp]; double mass = 0.04; + unsigned int Ls = 12; // global parameters Application::GlobalPar globalPar; @@ -65,7 +66,7 @@ int main(int argc, char *argv[]) std::string actionName = "DWF"; MAction::DWF::Par actionPar; actionPar.gauge = "gauge"; - actionPar.Ls = 12; + actionPar.Ls = Ls; actionPar.M5 = 1.8; actionPar.mass = mass; application.createModule(actionName, actionPar); @@ -83,37 +84,30 @@ int main(int argc, char *argv[]) std::string pos = "0 0 0 0"; std::string modName = "Ward Identity Test"; MAKE_POINT_PROP(pos, pointProp, solverName); - if (!(Environment::getInstance().hasModule(modName))) - { - MContraction::WardIdentity::Par wiPar; - wiPar.q = pointProp + "_5d"; - wiPar.q4d = pointProp; - wiPar.action = actionName; - wiPar.mass = mass; - application.createModule(modName, wiPar); - } + makeWITest(application, modName, pointProp, actionName, mass, Ls); - // Conserved current contractions with sequential insertion of vector + // Conserved current contractions with sequential insertion of vector/axial // current. - std::string q_x = "q_x"; - std::string q_y = "q_y"; - std::string q_z = "q_z"; - std::string q_t = "q_t"; - std::string mom = ZERO_MOM; - modName = "Sequential Ward Identity Test"; - MAKE_SEQUENTIAL_PROP(nt/2, pointProp, mom, q_x, solverName); - MAKE_SEQUENTIAL_PROP(nt/2, pointProp, mom, q_y, solverName); - MAKE_SEQUENTIAL_PROP(nt/2, pointProp, mom, q_z, solverName); - MAKE_SEQUENTIAL_PROP(nt/2, pointProp, mom, q_t, solverName); - if (!(Environment::getInstance().hasModule(modName))) - { - MContraction::WardIdentitySeq::Par wiPar; - wiPar.q_x = q_x; - wiPar.q_y = q_y; - wiPar.q_z = q_z; - wiPar.q_t = q_t; - application.createModule(modName, wiPar); - } + std::string mom = ZERO_MOM; + unsigned int t_J = nt/2; + std::string seqPropA = ADD_INDEX(pointProp + "_seq_A", t_J); + std::string seqPropV = ADD_INDEX(pointProp + "_seq_V", t_J); + std::string seqSrcA = seqPropA + "_src"; + std::string seqSrcV = seqPropV + "_src"; + std::string point5d = LABEL_5D(pointProp); + makeConservedSequentialSource(application, seqSrcA, point5d, + actionName, t_J, Current::Axial, Tp, mom); + makePropagator(application, seqPropA, seqSrcA, solverName); + makeConservedSequentialSource(application, seqSrcV, point5d, + actionName, t_J, Current::Vector, Tp, mom); + makePropagator(application, seqPropV, seqSrcV, solverName); + + std::string modNameA = "Axial Sequential Test"; + std::string modNameV = "Vector Sequential Test"; + makeSeqTest(application, modNameA, pointProp, seqPropA, + actionName, pos, t_J, Tp, Current::Axial, Ls); + makeSeqTest(application, modNameV, pointProp, seqPropV, + actionName, pos, t_J, Tp, Current::Vector, Ls); // execution application.saveParameterFile("ConservedCurrentTest.xml"); From 51d84ec057a80b9898369c7181dacce9979a945d Mon Sep 17 00:00:00 2001 From: Lanny91 Date: Fri, 28 Apr 2017 16:49:14 +0100 Subject: [PATCH 15/50] Bugfixes in Wilson 5D sequential conserved current insertion --- lib/qcd/action/fermion/WilsonFermion5D.cc | 24 ++++++++++------------- 1 file changed, 10 insertions(+), 14 deletions(-) diff --git a/lib/qcd/action/fermion/WilsonFermion5D.cc b/lib/qcd/action/fermion/WilsonFermion5D.cc index 99ff0dc1..bae5ae70 100644 --- a/lib/qcd/action/fermion/WilsonFermion5D.cc +++ b/lib/qcd/action/fermion/WilsonFermion5D.cc @@ -726,6 +726,8 @@ void WilsonFermion5D::SeqConservedCurrent(PropagatorField &q_in, conformable(q_in._grid, FermionGrid()); conformable(q_in._grid, q_out._grid); Lattice> ph(_FourDimGrid), coor(_FourDimGrid); + PropagatorField q_in_s(_FourDimGrid); + PropagatorField q_out_s(_FourDimGrid); Complex i(0.0, 1.0); // Momentum projection @@ -737,23 +739,17 @@ void WilsonFermion5D::SeqConservedCurrent(PropagatorField &q_in, } ph = exp((Real)(2*M_PI)*i*ph); - // Sequential insertion - Kernels::SeqConservedCurrentInternal(q_in, q_out, Umu, curr_type, - mu, ph, tmin, tmax); - - // Axial current sign. - if (curr_type == Current::Axial) + // Sequential insertion across 5th dimension + for (int s = 0; s < Ls; s++) { - SitePropagator result; - parallel_for(int sU = 0; sU < Umu._grid->oSites(); sU++) + ExtractSlice(q_in_s, q_in, s, 0); + Kernels::SeqConservedCurrentInternal(q_in_s, q_out_s, Umu, curr_type, + mu, ph, tmin, tmax); + if ((curr_type == Current::Axial) && (s < Ls/2)) { - int sF = sU * Ls; - for (int s = 0; s < Ls/2; s++) - { - vstream(q_out._odata[sF], -q_out._odata[sF]); - sF++; - } + q_out_s = -q_out_s; } + InsertSlice(q_out_s, q_out, s, 0); } } From 49331a3e72b6b644cd48424b8566440a7338fa64 Mon Sep 17 00:00:00 2001 From: Lanny91 Date: Fri, 28 Apr 2017 16:50:17 +0100 Subject: [PATCH 16/50] Minor improvements to Ward Identity checks --- .../Hadrons/Modules/MContraction/WardIdentity.hpp | 13 ++++++++++--- 1 file changed, 10 insertions(+), 3 deletions(-) diff --git a/extras/Hadrons/Modules/MContraction/WardIdentity.hpp b/extras/Hadrons/Modules/MContraction/WardIdentity.hpp index 41d8c6d1..d312bd4d 100644 --- a/extras/Hadrons/Modules/MContraction/WardIdentity.hpp +++ b/extras/Hadrons/Modules/MContraction/WardIdentity.hpp @@ -115,6 +115,10 @@ template void TWardIdentity::setup(void) { Ls_ = env().getObjectLs(par().q); + if (Ls_ != env().getObjectLs(par().action)) + { + HADRON_ERROR("Ls mismatch between quark action and propagator"); + } } // execution /////////////////////////////////////////////////////////////////// @@ -125,8 +129,8 @@ void TWardIdentity::execute(void) << "'." << std::endl; PropagatorField psi(env().getGrid()), tmp(env().getGrid()); - PropagatorField q = *env().template getObject(par().q); - PropagatorField q4d = *env().template getObject(par().q4d); + 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()), @@ -142,7 +146,7 @@ void TWardIdentity::execute(void) for (unsigned int mu = 0; mu < Nd; ++mu) { act.ContractConservedCurrent(q, q, tmp, Current::Vector, mu); - Vmu[mu] = trace(g5*tmp); + Vmu[mu] = trace(tmp); act.ContractConservedCurrent(q, q, tmp, Current::Axial, mu); Amu[mu] = trace(g5*tmp); } @@ -170,6 +174,9 @@ void TWardIdentity::execute(void) << 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; } END_MODULE_NAMESPACE From 77e0af9c2eca8ce816b6c6a54bc2d0edef26e213 Mon Sep 17 00:00:00 2001 From: Lanny91 Date: Fri, 5 May 2017 12:27:50 +0100 Subject: [PATCH 17/50] Compilation fix after merge - conserved current code not yet operational for vectorised 5D or Gparity Impl. --- lib/qcd/action/fermion/WilsonKernels.cc | 6 ++++++ 1 file changed, 6 insertions(+) diff --git a/lib/qcd/action/fermion/WilsonKernels.cc b/lib/qcd/action/fermion/WilsonKernels.cc index 5deec27c..8dc6bd8c 100644 --- a/lib/qcd/action/fermion/WilsonKernels.cc +++ b/lib/qcd/action/fermion/WilsonKernels.cc @@ -365,10 +365,16 @@ WilsonKernels::SeqConservedCurrentInternal(const PropagatorField &q_in, NO_CURR(GparityWilsonImplF); NO_CURR(GparityWilsonImplD); +NO_CURR(GparityWilsonImplFH); +NO_CURR(GparityWilsonImplDF); NO_CURR(DomainWallVec5dImplF); NO_CURR(DomainWallVec5dImplD); +NO_CURR(DomainWallVec5dImplFH); +NO_CURR(DomainWallVec5dImplDF); NO_CURR(ZDomainWallVec5dImplF); NO_CURR(ZDomainWallVec5dImplD); +NO_CURR(ZDomainWallVec5dImplFH); +NO_CURR(ZDomainWallVec5dImplDF); FermOpTemplateInstantiate(WilsonKernels); AdjointFermOpTemplateInstantiate(WilsonKernels); From d44cc204d166e19c0198ab9176061d3d8595930a Mon Sep 17 00:00:00 2001 From: Lanny91 Date: Fri, 12 May 2017 14:58:17 +0100 Subject: [PATCH 18/50] Added test module for sequential gamma matrix insertion --- extras/Hadrons/Modules.hpp | 1 + .../Modules/MUtilities/TestSeqGamma.hpp | 119 ++++++++++++++++++ extras/Hadrons/modules.inc | 1 + 3 files changed, 121 insertions(+) create mode 100644 extras/Hadrons/Modules/MUtilities/TestSeqGamma.hpp diff --git a/extras/Hadrons/Modules.hpp b/extras/Hadrons/Modules.hpp index 0286333c..dd6a6010 100644 --- a/extras/Hadrons/Modules.hpp +++ b/extras/Hadrons/Modules.hpp @@ -49,4 +49,5 @@ See the full license in the file "LICENSE" in the top level distribution directo #include #include #include +#include #include diff --git a/extras/Hadrons/Modules/MUtilities/TestSeqGamma.hpp b/extras/Hadrons/Modules/MUtilities/TestSeqGamma.hpp new file mode 100644 index 00000000..b3e99617 --- /dev/null +++ b/extras/Hadrons/Modules/MUtilities/TestSeqGamma.hpp @@ -0,0 +1,119 @@ +#ifndef Hadrons_TestSeqGamma_hpp_ +#define Hadrons_TestSeqGamma_hpp_ + +#include +#include +#include + +BEGIN_HADRONS_NAMESPACE + +/****************************************************************************** + * TestSeqGamma * + ******************************************************************************/ +BEGIN_MODULE_NAMESPACE(MUtilities) + +class TestSeqGammaPar: Serializable +{ +public: + GRID_SERIALIZABLE_CLASS_MEMBERS(TestSeqGammaPar, + std::string, q, + std::string, qSeq, + std::string, origin, + Gamma::Algebra, gamma, + unsigned int, t_g); +}; + +template +class TTestSeqGamma: public Module +{ +public: + TYPE_ALIASES(FImpl,); +public: + // constructor + TTestSeqGamma(const std::string name); + // destructor + virtual ~TTestSeqGamma(void) = default; + // dependency relation + virtual std::vector getInput(void); + virtual std::vector getOutput(void); + // setup + virtual void setup(void); + // execution + virtual void execute(void); +}; + +MODULE_REGISTER_NS(TestSeqGamma, TTestSeqGamma, MUtilities); + +/****************************************************************************** + * TTestSeqGamma implementation * + ******************************************************************************/ +// constructor ///////////////////////////////////////////////////////////////// +template +TTestSeqGamma::TTestSeqGamma(const std::string name) +: Module(name) +{} + +// dependencies/products /////////////////////////////////////////////////////// +template +std::vector TTestSeqGamma::getInput(void) +{ + std::vector in = {par().q, par().qSeq}; + + return in; +} + +template +std::vector TTestSeqGamma::getOutput(void) +{ + std::vector out = {getName()}; + + return out; +} + +// setup /////////////////////////////////////////////////////////////////////// +template +void TTestSeqGamma::setup(void) +{ + +} + +// execution /////////////////////////////////////////////////////////////////// +template +void TTestSeqGamma::execute(void) +{ + PropagatorField &q = *env().template getObject(par().q); + PropagatorField &qSeq = *env().template getObject(par().qSeq); + LatticeComplex c(env().getGrid()); + Gamma g5(Gamma::Algebra::Gamma5); + Gamma g(par().gamma); + SitePropagator qSite; + Complex test, check; + std::vector check_buf; + + // Check sequential insertion of gamma matrix gives same result as + // insertion of gamma at sink upon contraction. Assume q uses a point + // source. + std::vector siteCoord; + siteCoord = strToVec(par().origin); + peekSite(qSite, qSeq, siteCoord); + test = trace(g*qSite); + + c = trace(adj(g)*g5*adj(q)*g5*g*q); + sliceSum(c, check_buf, Tp); + check = TensorRemove(check_buf[par().t_g]); + + LOG(Message) << "Seq Result = " << abs(test) << std::endl; + LOG(Message) << "Reference = " << abs(check) << std::endl; + + // Check difference = 0 + check -= test; + + LOG(Message) << "Consistency check for sequential " << par().gamma + << " insertion = " << abs(check) << std::endl; +} + +END_MODULE_NAMESPACE + +END_HADRONS_NAMESPACE + +#endif // Hadrons_TestSeqGamma_hpp_ diff --git a/extras/Hadrons/modules.inc b/extras/Hadrons/modules.inc index 4ab51ce0..0364502a 100644 --- a/extras/Hadrons/modules.inc +++ b/extras/Hadrons/modules.inc @@ -30,5 +30,6 @@ modules_hpp =\ Modules/MSource/Wall.hpp \ Modules/MSource/Z2.hpp \ Modules/MUtilities/TestSeqConserved.hpp \ + Modules/MUtilities/TestSeqGamma.hpp \ Modules/Quark.hpp From 98f610ce5384f75883d7e4a31be54e55c6251410 Mon Sep 17 00:00:00 2001 From: Lanny91 Date: Fri, 12 May 2017 16:15:26 +0100 Subject: [PATCH 19/50] Reduced code duplication in hadron tests --- tests/hadrons/Test_hadrons.hpp | 96 ++++++++++++++++++++++++++ tests/hadrons/Test_hadrons_rarekaon.cc | 41 +++-------- 2 files changed, 107 insertions(+), 30 deletions(-) diff --git a/tests/hadrons/Test_hadrons.hpp b/tests/hadrons/Test_hadrons.hpp index c4dcedaf..c0a596b5 100644 --- a/tests/hadrons/Test_hadrons.hpp +++ b/tests/hadrons/Test_hadrons.hpp @@ -33,6 +33,30 @@ using namespace Hadrons; /******************************************************************************* * Macros to reduce code duplication. ******************************************************************************/ +// Common initialisation +#define HADRONS_DEFAULT_INIT \ + Grid_init(&argc, &argv); \ + HadronsLogError.Active(GridLogError.isActive()); \ + HadronsLogWarning.Active(GridLogWarning.isActive()); \ + HadronsLogMessage.Active(GridLogMessage.isActive()); \ + HadronsLogIterative.Active(GridLogIterative.isActive()); \ + HadronsLogDebug.Active(GridLogDebug.isActive()); \ + LOG(Message) << "Grid initialized" << std::endl; + +#define HADRONS_DEFAULT_GLOBALS(application) \ +{ \ + Application::GlobalPar globalPar; \ + globalPar.trajCounter.start = 1500; \ + globalPar.trajCounter.end = 1520; \ + globalPar.trajCounter.step = 20; \ + globalPar.seed = "1 2 3 4"; \ + globalPar.genetic.maxGen = 1000; \ + globalPar.genetic.maxCstGen = 200; \ + globalPar.genetic.popSize = 20; \ + globalPar.genetic.mutationRate = .1; \ + application.setPar(globalPar); \ +} + // Useful definitions #define ZERO_MOM "0. 0. 0. 0." #define INIT_INDEX(s, n) (std::string(s) + "_" + std::to_string(n)) @@ -73,10 +97,82 @@ using namespace Hadrons; makePropagator(application, propName, srcName, solver);\ } +/******************************************************************************* + * Action setups. + ******************************************************************************/ + +/******************************************************************************* + * Name: makeWilsonAction + * Parameters: application - main application that stores modules. + * actionName - name of action module to create. + * gaugeField - gauge field module. + * mass - quark mass. + * Returns: None. + ******************************************************************************/ +inline void makeWilsonAction(Application &application, std::string actionName, + std::string &gaugeField, double mass) +{ + if (!(Environment::getInstance().hasModule(actionName))) + { + MAction::Wilson::Par actionPar; + actionPar.gauge = gaugeField; + actionPar.mass = mass; + application.createModule(actionName, actionPar); + } +} + +/******************************************************************************* + * Name: makeDWFAction + * Parameters: application - main application that stores modules. + * actionName - name of action module to create. + * gaugeField - gauge field module. + * mass - quark mass. + * M5 - domain wall height. + * Ls - fifth dimension extent. + * Returns: None. + ******************************************************************************/ +inline void makeDWFAction(Application &application, std::string actionName, + std::string &gaugeField, double mass, double M5, + unsigned int Ls) +{ + if (!(Environment::getInstance().hasModule(actionName))) + { + MAction::DWF::Par actionPar; + actionPar.gauge = gaugeField; + actionPar.Ls = Ls; + actionPar.M5 = M5; + actionPar.mass = mass; + application.createModule(actionName, actionPar); + } +} + /******************************************************************************* * Functions for propagator construction. ******************************************************************************/ +/******************************************************************************* + * Name: makeRBPrecCGSolver + * Purpose: Make RBPrecCG solver module for specified action. + * Parameters: application - main application that stores modules. + * solverName - name of solver module to create. + * actionName - action module corresponding to propagators to be + * computed. + * residual - CG target residual. + * Returns: None. + ******************************************************************************/ +inline void makeRBPrecCGSolver(Application &application, std::string &solverName, + std::string &actionName, double residual = 1e-8) +{ + if (!(Environment::getInstance().hasModule(solverName))) + { + MSolver::RBPrecCG::Par solverPar; + solverPar.action = actionName; + solverPar.residual = residual; + application.createModule(solverName, + solverPar); + } +} + /******************************************************************************* * Name: makePointSource * Purpose: Construct point source and add to application module. diff --git a/tests/hadrons/Test_hadrons_rarekaon.cc b/tests/hadrons/Test_hadrons_rarekaon.cc index 9d35c1bc..1b5a45d9 100644 --- a/tests/hadrons/Test_hadrons_rarekaon.cc +++ b/tests/hadrons/Test_hadrons_rarekaon.cc @@ -51,13 +51,7 @@ int main(int argc, char *argv[]) configStem = argv[1]; // initialization ////////////////////////////////////////////////////////// - Grid_init(&argc, &argv); - HadronsLogError.Active(GridLogError.isActive()); - HadronsLogWarning.Active(GridLogWarning.isActive()); - HadronsLogMessage.Active(GridLogMessage.isActive()); - HadronsLogIterative.Active(GridLogIterative.isActive()); - HadronsLogDebug.Active(GridLogDebug.isActive()); - LOG(Message) << "Grid initialized" << std::endl; + HADRONS_DEFAULT_INIT; // run setup /////////////////////////////////////////////////////////////// Application application; @@ -74,46 +68,33 @@ int main(int argc, char *argv[]) unsigned int n_noise = 1; unsigned int nt = 32; bool do_disconnected(false); + Gamma::Algebra gT = Gamma::Algebra::GammaT; + unsigned int Ls = 16; + double M5 = 1.8; // Global parameters. - Application::GlobalPar globalPar; - globalPar.trajCounter.start = 1500; - globalPar.trajCounter.end = 1520; - globalPar.trajCounter.step = 20; - globalPar.seed = "1 2 3 4"; - globalPar.genetic.maxGen = 1000; - globalPar.genetic.maxCstGen = 200; - globalPar.genetic.popSize = 20; - globalPar.genetic.mutationRate = .1; - application.setPar(globalPar); + HADRONS_DEFAULT_GLOBALS(application); // gauge field + std::string gaugeField = "gauge"; if (configStem == "None") { - application.createModule("gauge"); + application.createModule(gaugeField); } else { MGauge::Load::Par gaugePar; gaugePar.file = configStem; - application.createModule("gauge", gaugePar); + application.createModule(gaugeField, gaugePar); } for (unsigned int i = 0; i < flavour.size(); ++i) { // actions - MAction::DWF::Par actionPar; - actionPar.gauge = "gauge"; - actionPar.Ls = 16; - actionPar.M5 = 1.8; - actionPar.mass = mass[i]; - application.createModule("DWF_" + flavour[i], actionPar); + std::string actionName = "DWF_" + flavour[i]; + makeDWFAction(application, actionName, gaugeField, mass[i], M5, Ls); // solvers - MSolver::RBPrecCG::Par solverPar; - solverPar.action = "DWF_" + flavour[i]; - solverPar.residual = 1.0e-8; - application.createModule(solvers[i], - solverPar); + makeRBPrecCGSolver(application, solvers[i], actionName); } // Create noise propagators for loops. From c2010f21aba12b3f2fd7166211b6c3243b428ed9 Mon Sep 17 00:00:00 2001 From: Lanny91 Date: Fri, 12 May 2017 16:23:01 +0100 Subject: [PATCH 20/50] Added sequential propagator test for gamma matrix insertion --- tests/hadrons/Test_hadrons.hpp | 42 ++++++++++-- tests/hadrons/Test_hadrons_rarekaon.cc | 10 +-- tests/hadrons/Test_hadrons_seq_gamma.cc | 89 +++++++++++++++++++++++++ 3 files changed, 131 insertions(+), 10 deletions(-) create mode 100644 tests/hadrons/Test_hadrons_seq_gamma.cc diff --git a/tests/hadrons/Test_hadrons.hpp b/tests/hadrons/Test_hadrons.hpp index c0a596b5..61e90bac 100644 --- a/tests/hadrons/Test_hadrons.hpp +++ b/tests/hadrons/Test_hadrons.hpp @@ -82,11 +82,11 @@ using namespace Hadrons; MAKE_3MOM_WALL_PROP(tW, ZERO_MOM, propName, solver) // Sequential source macros -#define MAKE_SEQUENTIAL_PROP(tS, qSrc, mom, propName, solver)\ +#define MAKE_SEQUENTIAL_PROP(tS, qSrc, mom, seqPropName, solver, gamma)\ {\ - std::string srcName = ADD_INDEX(qSrc + "_seq", tS);\ - makeSequentialSource(application, srcName, qSrc, tS, mom);\ - makePropagator(application, propName, srcName, solver);\ + std::string srcName = seqPropName + "_src";\ + makeSequentialSource(application, srcName, qSrc, tS, gamma, mom);\ + makePropagator(application, seqPropName, srcName, solver);\ } // Point source macros @@ -205,6 +205,7 @@ inline void makePointSource(Application &application, std::string srcName, ******************************************************************************/ inline void makeSequentialSource(Application &application, std::string srcName, std::string qSrc, unsigned int tS, + Gamma::Algebra gamma = Gamma::Algebra::GammaT, std::string mom = ZERO_MOM) { // If the source already exists, don't make the module again. @@ -215,7 +216,7 @@ inline void makeSequentialSource(Application &application, std::string srcName, seqPar.tA = tS; seqPar.tB = tS; seqPar.mom = mom; - seqPar.gamma = Gamma::Algebra::GammaT; + seqPar.gamma = gamma; application.createModule(srcName, seqPar); } } @@ -579,3 +580,34 @@ inline void makeSeqTest(Application &application, std::string &modName, application.createModule(modName, seqPar); } } + +/******************************************************************************* + * Name: makeSeqGamComparison + * Purpose: Create module to compare sequential insertion of gamma matrix + * against sink contraction and add to application module. + * Parameters: application - main application that stores modules. + * modName - name of module to create. + * propName - 4D quark propagator. + * seqProp - 4D quark propagator with sequential insertion of + * gamma matrix. + * gamma - Inserted gamma matrix. + * t_g - time at which gamma matrix is inserted + * sequentially. + * Returns: None. + ******************************************************************************/ +inline void makeSeqGamComparison(Application &application, std::string &modName, + std::string &propName, std::string &seqProp, + std::string &origin, Gamma::Algebra gamma, + unsigned int t_g) +{ + if (!(Environment::getInstance().hasModule(modName))) + { + MUtilities::TestSeqGamma::Par seqPar; + seqPar.q = propName; + seqPar.qSeq = seqProp; + seqPar.origin = origin; + seqPar.t_g = t_g; + seqPar.gamma = gamma; + application.createModule(modName, seqPar); + } +} diff --git a/tests/hadrons/Test_hadrons_rarekaon.cc b/tests/hadrons/Test_hadrons_rarekaon.cc index 1b5a45d9..3a642f24 100644 --- a/tests/hadrons/Test_hadrons_rarekaon.cc +++ b/tests/hadrons/Test_hadrons_rarekaon.cc @@ -212,10 +212,10 @@ int main(int argc, char *argv[]) std::string q_KsCs_mq = LABEL_3PT("Q_KsCs_mq", tK, tJ); std::string q_pilCl_q = LABEL_3PT("Q_pilCl_q", tpi, tJ); std::string q_pilCl_mq = LABEL_3PT("Q_pilCl_mq", tpi, tJ); - MAKE_SEQUENTIAL_PROP(tJ, q_Kl_0, qmom, q_KlCl_q, solvers[light]); - MAKE_SEQUENTIAL_PROP(tJ, q_Ks_k, mqmom, q_KsCs_mq, solvers[strange]); - MAKE_SEQUENTIAL_PROP(tJ, q_pil_p, qmom, q_pilCl_q, solvers[light]); - MAKE_SEQUENTIAL_PROP(tJ, q_pil_0, mqmom, q_pilCl_mq, solvers[light]); + MAKE_SEQUENTIAL_PROP(tJ, q_Kl_0, qmom, q_KlCl_q, solvers[light], gT); + MAKE_SEQUENTIAL_PROP(tJ, q_Ks_k, mqmom, q_KsCs_mq, solvers[strange], gT); + MAKE_SEQUENTIAL_PROP(tJ, q_pil_p, qmom, q_pilCl_q, solvers[light], gT); + MAKE_SEQUENTIAL_PROP(tJ, q_pil_0, mqmom, q_pilCl_mq, solvers[light], gT); /******************************************************************* * CONTRACTIONS: pi and K 3pt contractions with current insertion. @@ -271,7 +271,7 @@ int main(int argc, char *argv[]) std::string loop_qCq = LABEL_3PT(loop_stem + flavour[f], tJ, nn); std::string loop_qCq_res = loop_qCq + "_res"; MAKE_SEQUENTIAL_PROP(tJ, noiseRes[f][nn], qmom, - loop_qCq_res, solvers[f]); + loop_qCq_res, solvers[f], gT); makeLoop(application, loop_qCq, eta, loop_qCq_res); /******************************************************* diff --git a/tests/hadrons/Test_hadrons_seq_gamma.cc b/tests/hadrons/Test_hadrons_seq_gamma.cc new file mode 100644 index 00000000..22c35ecb --- /dev/null +++ b/tests/hadrons/Test_hadrons_seq_gamma.cc @@ -0,0 +1,89 @@ +/******************************************************************************* + Grid physics library, www.github.com/paboyle/Grid + + Source file: tests/hadrons/Test_hadrons_seq_gamma.cc + + Copyright (C) 2017 + + Author: Andrew Lawson + + This program is free software; you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation; either version 2 of the License, or + (at your option) any later version. + + This program is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License along + with this program; if not, write to the Free Software Foundation, Inc., + 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. + + See the full license in the file "LICENSE" in the top level distribution + directory. + *******************************************************************************/ + +#include "Test_hadrons.hpp" + +using namespace Grid; +using namespace QCD; +using namespace Hadrons; + +/******************************************************************************* + * Consistency test for sequential gamma insertion. + ******************************************************************************/ + +int main(int argc, char *argv[]) +{ + // initialization ////////////////////////////////////////////////////////// + HADRONS_DEFAULT_INIT; + + // run setup /////////////////////////////////////////////////////////////// + Application application; + unsigned int nt = GridDefaultLatt()[Tp]; + unsigned int tS = nt / 2; + unsigned int Ls = 12; + double mass = 0.04; + double M5 = 1.8; + + // global parameters + HADRONS_DEFAULT_GLOBALS(application); + + // gauge field + std::string gaugeField = "gauge"; + application.createModule(gaugeField); + + // action + std::string actionName = "DWF"; + makeDWFAction(application, actionName, gaugeField, mass, M5, Ls); + + // solver + std::string solverName = "CG"; + makeRBPrecCGSolver(application, solverName, actionName); + + // test sequential propagator, with g5 insertion. + Gamma::Algebra g = Gamma::Algebra::Gamma5; + std::string pointProp = "q_0"; + std::string point5d = LABEL_5D(pointProp); + std::string origin = "0 0 0 0"; + MAKE_POINT_PROP(origin, pointProp, solverName); + + std::string seqProp = ADD_INDEX(pointProp + "_seqg5", tS); + std::string seqSrc = seqProp + "_src"; + MAKE_SEQUENTIAL_PROP(tS, pointProp, ZERO_MOM, seqProp, solverName, g); + + std::string modName = "Test g5 sequential insertion"; + makeSeqGamComparison(application, modName, pointProp, seqProp, origin, g, tS); + + // execution + application.saveParameterFile("SeqGamma5Test.xml"); + application.run(); + + // epilogue + LOG(Message) << "Grid is finalizing now" << std::endl; + Grid_finalize(); + + return EXIT_SUCCESS; +} \ No newline at end of file From 34332fe3934754025446cec92b6c099c6828df9f Mon Sep 17 00:00:00 2001 From: Lanny91 Date: Fri, 12 May 2017 16:30:43 +0100 Subject: [PATCH 21/50] Improvement to sequential conserved current insertion tests --- .../Modules/MUtilities/TestSeqConserved.hpp | 45 +++-- tests/hadrons/Test_hadrons.hpp | 31 ++-- .../hadrons/Test_hadrons_conserved_current.cc | 156 +++++++++++------- 3 files changed, 131 insertions(+), 101 deletions(-) diff --git a/extras/Hadrons/Modules/MUtilities/TestSeqConserved.hpp b/extras/Hadrons/Modules/MUtilities/TestSeqConserved.hpp index 0730b8ed..3ae1b8b0 100644 --- a/extras/Hadrons/Modules/MUtilities/TestSeqConserved.hpp +++ b/extras/Hadrons/Modules/MUtilities/TestSeqConserved.hpp @@ -41,7 +41,6 @@ BEGIN_HADRONS_NAMESPACE * options: - q: point source propagator, 5D if available (string) - - q4d: 4D point source propagator, duplicate of q if q is 4D (string) - qSeq: result of sequential insertion of conserved current using q (string) - action: action used for computation of q (string) - origin: string giving point source origin of q (string) @@ -60,7 +59,6 @@ class TestSeqConservedPar: Serializable public: GRID_SERIALIZABLE_CLASS_MEMBERS(TestSeqConservedPar, std::string, q, - std::string, q4d, std::string, qSeq, std::string, action, std::string, origin, @@ -103,8 +101,7 @@ TTestSeqConserved::TTestSeqConserved(const std::string name) template std::vector TTestSeqConserved::getInput(void) { - std::vector in = {par().q, par().q4d, - par().qSeq, par().action}; + std::vector in = {par().q, par().qSeq, par().action}; return in; } @@ -121,7 +118,11 @@ std::vector TTestSeqConserved::getOutput(void) template void TTestSeqConserved::setup(void) { - + auto Ls = env().getObjectLs(par().q); + if (Ls != env().getObjectLs(par().action)) + { + HADRON_ERROR("Ls mismatch between quark action and propagator"); + } } // execution /////////////////////////////////////////////////////////////////// @@ -130,33 +131,43 @@ void TTestSeqConserved::execute(void) { PropagatorField tmp(env().getGrid()); PropagatorField &q = *env().template getObject(par().q); - PropagatorField &q4d = *env().template getObject(par().q4d); PropagatorField &qSeq = *env().template getObject(par().qSeq); FMat &act = *(env().template getObject(par().action)); Gamma g5(Gamma::Algebra::Gamma5); + Gamma::Algebra gA = (par().curr == Current::Axial) ? + Gamma::Algebra::Gamma5 : + Gamma::Algebra::Identity; + Gamma g(gA); SitePropagator qSite; - LatticeComplex c(env().getGrid()); - Complex seq_res, check_res; - std::vector check_buf; + Complex test_S, test_V, check_S, check_V; + std::vector check_buf; // Check sequential insertion of current gives same result as conserved // current sink upon contraction. Assume q uses a point source. std::vector siteCoord; siteCoord = strToVec(par().origin); - peekSite(qSite, q, siteCoord); - seq_res = trace(g5*qSite); + peekSite(qSite, qSeq, siteCoord); + test_S = trace(qSite*g); + test_V = trace(qSite*g*Gamma::gmu[par().mu]); act.ContractConservedCurrent(q, q, tmp, par().curr, par().mu); - c = trace(tmp); - sliceSum(c, check_buf, Tp); - check_res = TensorRemove(check_buf[par().t_J]); + sliceSum(tmp, check_buf, Tp); + check_S = TensorRemove(trace(check_buf[par().t_J]*g)); + check_V = TensorRemove(trace(check_buf[par().t_J]*g*Gamma::gmu[par().mu])); + + LOG(Message) << "Test S = " << abs(test_S) << std::endl; + LOG(Message) << "Test V = " << abs(test_V) << std::endl; + LOG(Message) << "Check S = " << abs(check_S) << std::endl; + LOG(Message) << "Check V = " << abs(check_V) << std::endl; // Check difference = 0 - check_res -= seq_res; + check_S -= test_S; + check_V -= test_V; LOG(Message) << "Consistency check for sequential conserved " - << par().curr << " current insertion = " << abs(check_res) - << std::endl; + << par().curr << " current insertion: " << std::endl; + LOG(Message) << "Check S = " << abs(check_S) << std::endl; + LOG(Message) << "Check V = " << abs(check_V) << std::endl; } END_MODULE_NAMESPACE diff --git a/tests/hadrons/Test_hadrons.hpp b/tests/hadrons/Test_hadrons.hpp index 61e90bac..1b038388 100644 --- a/tests/hadrons/Test_hadrons.hpp +++ b/tests/hadrons/Test_hadrons.hpp @@ -538,39 +538,30 @@ inline void makeWITest(Application &application, std::string &modName, } /******************************************************************************* - * Name: makeSeqTest - * Purpose: Create module to test sequential insertion of conserved current - * and add to application module. + * Name: makeSeqCurrComparison + * Purpose: Create module to compare sequential insertion of conserved current + * against sink contraction and add to application module. * Parameters: application - main application that stores modules. * modName - name of module to create. - * propName - 4D quark propagator. - * seqProp - 4D quark propagator with sequential insertion of + * propName - quark propagator (point source), 5D if available. + * seqName - 4D quark propagator with sequential insertion of * conserved current. * actionName - action used to compute quark propagators. + * origin - origin of point source propagator. * t_J - time at which sequential current is inserted. * mu - Lorentz index of sequential current. * curr - type of conserved current inserted. - * Ls - length of 5th dimension (default = 1). * Returns: None. ******************************************************************************/ -inline void makeSeqTest(Application &application, std::string &modName, - std::string &propName, std::string &seqName, - std::string &actionName, std::string &origin, - unsigned int t_J, unsigned int mu, Current curr, - unsigned int Ls = 1) +inline void makeSeqCurrComparison(Application &application, std::string &modName, + std::string &propName, std::string &seqName, + std::string &actionName, std::string &origin, + unsigned int t_J, unsigned int mu, Current curr) { if (!(Environment::getInstance().hasModule(modName))) { MUtilities::TestSeqConserved::Par seqPar; - if (Ls > 1) - { - seqPar.q = LABEL_5D(propName); - } - else - { - seqPar.q = propName; - } - seqPar.q4d = propName; + seqPar.q = propName; seqPar.qSeq = seqName; seqPar.action = actionName; seqPar.origin = origin; diff --git a/tests/hadrons/Test_hadrons_conserved_current.cc b/tests/hadrons/Test_hadrons_conserved_current.cc index a11a3530..080fef73 100644 --- a/tests/hadrons/Test_hadrons_conserved_current.cc +++ b/tests/hadrons/Test_hadrons_conserved_current.cc @@ -30,84 +30,112 @@ using namespace Grid; using namespace Hadrons; +inline void setupSeqCurrTests(Application &application, std::string modStem, + std::string &pointProp, std::string &seqStem, + std::string &actionName, std::string &solverName, + std::string &origin, Current curr, + unsigned int t_J, unsigned int mu, + unsigned int Ls = 1) +{ + std::string modName = ADD_INDEX(modStem, mu); + std::string seqProp = ADD_INDEX(seqStem, mu); + std::string seqSrc = seqProp + "_src"; + + // 5D actions require 5D propagator as input for conserved current + // insertions. + std::string propIn; + if (Ls > 1) + { + propIn = LABEL_5D(pointProp); + } + else + { + propIn = pointProp; + } + + makeConservedSequentialSource(application, seqSrc, propIn, + actionName, t_J, curr, mu); + makePropagator(application, seqProp, seqSrc, solverName); + makeSeqCurrComparison(application, modName, propIn, seqProp, + actionName, origin, t_J, mu, curr); +} + +inline void setupWardIdentityTests(Application &application, + std::string &actionName, + double mass, + unsigned int Ls = 1, + bool perform_axial_tests = false) +{ + // solver + std::string solverName = actionName + "_CG"; + makeRBPrecCGSolver(application, solverName, actionName); + + unsigned int nt = GridDefaultLatt()[Tp]; + unsigned int t_J = nt/2; + + /*************************************************************************** + * Conserved current sink contractions: use a single point propagator for + * the Ward Identity test. + **************************************************************************/ + std::string pointProp = actionName + "_q_0"; + 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); + + /*************************************************************************** + * Conserved current tests with sequential insertion of vector/axial + * current. If above Ward Identity passes, sufficient to test sequential + * insertion of conserved current agrees with contracted version. + **************************************************************************/ + // Compare sequential insertion to contraction. Should be enough to perform + // for time and one space component. + std::string seqStem = ADD_INDEX(pointProp + "seq_V", t_J); + std::string modStem = actionName + " Vector Sequential Test mu"; + setupSeqCurrTests(application, modStem, pointProp, seqStem, actionName, + solverName, origin, Current::Vector, t_J, Tp, Ls); + setupSeqCurrTests(application, modStem, pointProp, seqStem, actionName, + solverName, origin, Current::Vector, t_J, Xp, Ls); + + // Perform axial tests only if partially-conserved axial current exists for + // the action. + if (perform_axial_tests) + { + seqStem = ADD_INDEX(pointProp + "seq_A", t_J); + modStem = actionName + " Axial Sequential Test mu"; + setupSeqCurrTests(application, modStem, pointProp, seqStem, actionName, + solverName, origin, Current::Axial, t_J, Tp, Ls); + setupSeqCurrTests(application, modStem, pointProp, seqStem, actionName, + solverName, origin, Current::Axial, t_J, Xp, Ls); + } +} + int main(int argc, char *argv[]) { // initialization ////////////////////////////////////////////////////////// - Grid_init(&argc, &argv); - HadronsLogError.Active(GridLogError.isActive()); - HadronsLogWarning.Active(GridLogWarning.isActive()); - HadronsLogMessage.Active(GridLogMessage.isActive()); - HadronsLogIterative.Active(GridLogIterative.isActive()); - HadronsLogDebug.Active(GridLogDebug.isActive()); - LOG(Message) << "Grid initialized" << std::endl; - + HADRONS_DEFAULT_INIT; + // run setup /////////////////////////////////////////////////////////////// Application application; - unsigned int nt = GridDefaultLatt()[Tp]; double mass = 0.04; + double M5 = 1.8; unsigned int Ls = 12; // global parameters - Application::GlobalPar globalPar; - globalPar.trajCounter.start = 1500; - globalPar.trajCounter.end = 1520; - globalPar.trajCounter.step = 20; - globalPar.seed = "1 2 3 4"; - globalPar.genetic.maxGen = 1000; - globalPar.genetic.maxCstGen = 200; - globalPar.genetic.popSize = 20; - globalPar.genetic.mutationRate = .1; - application.setPar(globalPar); + HADRONS_DEFAULT_GLOBALS(application); // gauge field - application.createModule("gauge"); + std::string gaugeField = "gauge"; + application.createModule(gaugeField); - // action + // Setup each action and the conserved current tests relevant to it. std::string actionName = "DWF"; - MAction::DWF::Par actionPar; - actionPar.gauge = "gauge"; - actionPar.Ls = Ls; - actionPar.M5 = 1.8; - actionPar.mass = mass; - application.createModule(actionName, actionPar); + makeDWFAction(application, actionName, gaugeField, mass, M5, Ls); + setupWardIdentityTests(application, actionName, mass, Ls, true); - // solver - std::string solverName = "CG"; - MSolver::RBPrecCG::Par solverPar; - solverPar.action = actionName; - solverPar.residual = 1.0e-8; - application.createModule(solverName, - solverPar); - - // Conserved current sink contractions: use a single point propagator. - std::string pointProp = "q_0"; - std::string pos = "0 0 0 0"; - std::string modName = "Ward Identity Test"; - MAKE_POINT_PROP(pos, pointProp, solverName); - makeWITest(application, modName, pointProp, actionName, mass, Ls); - - // Conserved current contractions with sequential insertion of vector/axial - // current. - std::string mom = ZERO_MOM; - unsigned int t_J = nt/2; - std::string seqPropA = ADD_INDEX(pointProp + "_seq_A", t_J); - std::string seqPropV = ADD_INDEX(pointProp + "_seq_V", t_J); - std::string seqSrcA = seqPropA + "_src"; - std::string seqSrcV = seqPropV + "_src"; - std::string point5d = LABEL_5D(pointProp); - makeConservedSequentialSource(application, seqSrcA, point5d, - actionName, t_J, Current::Axial, Tp, mom); - makePropagator(application, seqPropA, seqSrcA, solverName); - makeConservedSequentialSource(application, seqSrcV, point5d, - actionName, t_J, Current::Vector, Tp, mom); - makePropagator(application, seqPropV, seqSrcV, solverName); - - std::string modNameA = "Axial Sequential Test"; - std::string modNameV = "Vector Sequential Test"; - makeSeqTest(application, modNameA, pointProp, seqPropA, - actionName, pos, t_J, Tp, Current::Axial, Ls); - makeSeqTest(application, modNameV, pointProp, seqPropV, - actionName, pos, t_J, Tp, Current::Vector, Ls); + actionName = "Wilson"; + makeWilsonAction(application, actionName, gaugeField, mass); + setupWardIdentityTests(application, actionName, mass); // execution application.saveParameterFile("ConservedCurrentTest.xml"); From 08b314fd0fadd012492710bf57fd17b37ea9cf54 Mon Sep 17 00:00:00 2001 From: Lanny91 Date: Thu, 18 May 2017 13:16:14 +0100 Subject: [PATCH 22/50] 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 From eec79e0a1e8cdaf58c55547d520e2b9e9a894898 Mon Sep 17 00:00:00 2001 From: Lanny91 Date: Mon, 5 Jun 2017 11:55:41 +0100 Subject: [PATCH 23/50] Ward Identity test improvements and conserved current bug fixes --- .../Modules/MContraction/WardIdentity.hpp | 35 ++++++++++++------- lib/qcd/action/fermion/WilsonFermion5D.cc | 10 ++++-- lib/qcd/action/fermion/WilsonKernels.cc | 21 +++++------ 3 files changed, 42 insertions(+), 24 deletions(-) diff --git a/extras/Hadrons/Modules/MContraction/WardIdentity.hpp b/extras/Hadrons/Modules/MContraction/WardIdentity.hpp index fa51ce95..7fc7d15d 100644 --- a/extras/Hadrons/Modules/MContraction/WardIdentity.hpp +++ b/extras/Hadrons/Modules/MContraction/WardIdentity.hpp @@ -143,6 +143,7 @@ void TWardIdentity::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::execute(void) { LatticeComplex PP(env().getGrid()), axial_defect(env().getGrid()), PJ5q(env().getGrid()); + std::vector axial_buf; - // Compute D_mu A_mu, D is backwards derivative. + // Compute , 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 for 5D (zero for 4D) and . PJ5q = zero; if (Ls_ > 1) { - // PP + // 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 + // 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::execute(void) 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; + // Test ward identity = 2m + 2 + 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; + } } } diff --git a/lib/qcd/action/fermion/WilsonFermion5D.cc b/lib/qcd/action/fermion/WilsonFermion5D.cc index f616a080..3bbc03b4 100644 --- a/lib/qcd/action/fermion/WilsonFermion5D.cc +++ b/lib/qcd/action/fermion/WilsonFermion5D.cc @@ -700,8 +700,14 @@ void WilsonFermion5D::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; + } } } diff --git a/lib/qcd/action/fermion/WilsonKernels.cc b/lib/qcd/action/fermion/WilsonKernels.cc index 8dc6bd8c..802c0940 100644 --- a/lib/qcd/action/fermion/WilsonKernels.cc +++ b/lib/qcd/action/fermion/WilsonKernels.cc @@ -286,8 +286,9 @@ void WilsonKernels::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 void WilsonKernels::ContractConservedCurrentInternal(const PropagatorField &q_in_1, @@ -300,13 +301,13 @@ void WilsonKernels::ContractConservedCurrentInternal(const PropagatorField Gamma g5(Gamma::Algebra::Gamma5); PropagatorField tmp(q_out._grid); GaugeLinkField Umu(U._grid); - Umu = PeekIndex(U, mu); + Umu = PeekIndex(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::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(U, mu); + Umu = PeekIndex(U, mu); Lattice> 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); } From 622a21bec673ccf3a3b895584678afacd1f59c4b Mon Sep 17 00:00:00 2001 From: Lanny91 Date: Mon, 5 Jun 2017 15:55:32 +0100 Subject: [PATCH 24/50] Improvements to sequential conserved current test and small bugfix. --- .../Modules/MUtilities/TestSeqConserved.hpp | 18 ++++++++++++------ lib/qcd/action/fermion/WilsonKernels.cc | 1 + 2 files changed, 13 insertions(+), 6 deletions(-) diff --git a/extras/Hadrons/Modules/MUtilities/TestSeqConserved.hpp b/extras/Hadrons/Modules/MUtilities/TestSeqConserved.hpp index 3ae1b8b0..eccb00cc 100644 --- a/extras/Hadrons/Modules/MUtilities/TestSeqConserved.hpp +++ b/extras/Hadrons/Modules/MUtilities/TestSeqConserved.hpp @@ -140,7 +140,8 @@ void TTestSeqConserved::execute(void) Gamma g(gA); SitePropagator qSite; Complex test_S, test_V, check_S, check_V; - std::vector check_buf; + std::vector check_buf; + LatticeComplex c(env().getGrid()); // Check sequential insertion of current gives same result as conserved // current sink upon contraction. Assume q uses a point source. @@ -151,9 +152,14 @@ void TTestSeqConserved::execute(void) test_V = trace(qSite*g*Gamma::gmu[par().mu]); act.ContractConservedCurrent(q, q, tmp, par().curr, par().mu); - sliceSum(tmp, check_buf, Tp); - check_S = TensorRemove(trace(check_buf[par().t_J]*g)); - check_V = TensorRemove(trace(check_buf[par().t_J]*g*Gamma::gmu[par().mu])); + + c = trace(tmp*g); + sliceSum(c, check_buf, Tp); + check_S = TensorRemove(check_buf[par().t_J]); + + c = trace(tmp*g*Gamma::gmu[par().mu]); + sliceSum(c, check_buf, Tp); + check_V = TensorRemove(check_buf[par().t_J]); LOG(Message) << "Test S = " << abs(test_S) << std::endl; LOG(Message) << "Test V = " << abs(test_V) << std::endl; @@ -166,8 +172,8 @@ void TTestSeqConserved::execute(void) LOG(Message) << "Consistency check for sequential conserved " << par().curr << " current insertion: " << std::endl; - LOG(Message) << "Check S = " << abs(check_S) << std::endl; - LOG(Message) << "Check V = " << abs(check_V) << std::endl; + LOG(Message) << "Diff S = " << abs(check_S) << std::endl; + LOG(Message) << "Diff V = " << abs(check_V) << std::endl; } END_MODULE_NAMESPACE diff --git a/lib/qcd/action/fermion/WilsonKernels.cc b/lib/qcd/action/fermion/WilsonKernels.cc index 802c0940..8d5406f4 100644 --- a/lib/qcd/action/fermion/WilsonKernels.cc +++ b/lib/qcd/action/fermion/WilsonKernels.cc @@ -327,6 +327,7 @@ void WilsonKernels::SeqConservedCurrentInternal(const PropagatorField &q_i GaugeLinkField Umu(U._grid); Umu = PeekIndex(U, mu); Lattice> t(q_in._grid); + LatticeCoordinate(t, mu); tmp = this->CovShiftForward(Umu, mu, q_in)*ph; tmp = where((t >= tmin) and (t <= tmax), tmp, 0.*tmp); From c504b4dbad611b8f36599fb5d6202a85b465134d Mon Sep 17 00:00:00 2001 From: Lanny91 Date: Mon, 5 Jun 2017 15:56:43 +0100 Subject: [PATCH 25/50] Code cleaning --- extras/Hadrons/Modules/MContraction/WardIdentity.hpp | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/extras/Hadrons/Modules/MContraction/WardIdentity.hpp b/extras/Hadrons/Modules/MContraction/WardIdentity.hpp index 7fc7d15d..fb2ea173 100644 --- a/extras/Hadrons/Modules/MContraction/WardIdentity.hpp +++ b/extras/Hadrons/Modules/MContraction/WardIdentity.hpp @@ -128,12 +128,11 @@ void TWardIdentity::execute(void) LOG(Message) << "Performing Ward Identity checks for quark '" << par().q << "'." << std::endl; - PropagatorField psi(env().getGrid()), tmp(env().getGrid()), - vector_WI(env().getGrid()); + PropagatorField tmp(env().getGrid()), vector_WI(env().getGrid()); PropagatorField &q = *env().template getObject(par().q); FMat &act = *(env().template getObject(par().action)); Gamma g5(Gamma::Algebra::Gamma5); - + // Compute D_mu V_mu, D here is backward derivative. vector_WI = zero; for (unsigned int mu = 0; mu < Nd; ++mu) @@ -149,6 +148,7 @@ void TWardIdentity::execute(void) if (par().test_axial) { + PropagatorField psi(env().getGrid()); LatticeComplex PP(env().getGrid()), axial_defect(env().getGrid()), PJ5q(env().getGrid()); std::vector axial_buf; @@ -159,7 +159,7 @@ void TWardIdentity::execute(void) { act.ContractConservedCurrent(q, q, tmp, Current::Axial, mu); tmp -= Cshift(tmp, mu, -1); - axial_defect += 2.*trace(g5*tmp); + axial_defect += trace(g5*tmp); } // Get for 5D (zero for 4D) and . @@ -191,7 +191,7 @@ void TWardIdentity::execute(void) LOG(Message) << "|PJ5q|^2 = " << norm2(PJ5q) << std::endl; LOG(Message) << "Axial Ward Identity defect Delta_mu A_mu = " << norm2(axial_defect) << std::endl; - + // Axial defect by timeslice. axial_defect -= 2.*(par().mass*PP + PJ5q); LOG(Message) << "Check Axial defect by timeslice" << std::endl; From e5c8b7369e2cb259379d987260cf21f2b96e404f Mon Sep 17 00:00:00 2001 From: Lanny91 Date: Tue, 6 Jun 2017 14:19:10 +0100 Subject: [PATCH 26/50] Boundary condition option in quark actions for hadrons tests. --- tests/hadrons/Test_hadrons.hpp | 11 +++++++++-- 1 file changed, 9 insertions(+), 2 deletions(-) diff --git a/tests/hadrons/Test_hadrons.hpp b/tests/hadrons/Test_hadrons.hpp index 6dbe3425..a554425d 100644 --- a/tests/hadrons/Test_hadrons.hpp +++ b/tests/hadrons/Test_hadrons.hpp @@ -107,16 +107,20 @@ using namespace Hadrons; * actionName - name of action module to create. * gaugeField - gauge field module. * mass - quark mass. + * boundary - fermion boundary conditions (default to periodic + * space, antiperiodic time). * Returns: None. ******************************************************************************/ inline void makeWilsonAction(Application &application, std::string actionName, - std::string &gaugeField, double mass) + std::string &gaugeField, double mass, + std::string boundary = "1 1 1 -1") { if (!(Environment::getInstance().hasModule(actionName))) { MAction::Wilson::Par actionPar; actionPar.gauge = gaugeField; actionPar.mass = mass; + actionPar.boundary = boundary; application.createModule(actionName, actionPar); } } @@ -129,11 +133,13 @@ inline void makeWilsonAction(Application &application, std::string actionName, * mass - quark mass. * M5 - domain wall height. * Ls - fifth dimension extent. + * boundary - fermion boundary conditions (default to periodic + * space, antiperiodic time). * Returns: None. ******************************************************************************/ inline void makeDWFAction(Application &application, std::string actionName, std::string &gaugeField, double mass, double M5, - unsigned int Ls) + unsigned int Ls, std::string boundary = "1 1 1 -1") { if (!(Environment::getInstance().hasModule(actionName))) { @@ -142,6 +148,7 @@ inline void makeDWFAction(Application &application, std::string actionName, actionPar.Ls = Ls; actionPar.M5 = M5; actionPar.mass = mass; + actionPar.boundary = boundary; application.createModule(actionName, actionPar); } } From 8d442b502dc59f7fe4407b02142677299ac63740 Mon Sep 17 00:00:00 2001 From: Lanny91 Date: Tue, 6 Jun 2017 17:06:40 +0100 Subject: [PATCH 27/50] Sequential current fix for spacial indices. --- lib/qcd/action/fermion/WilsonKernels.cc | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/lib/qcd/action/fermion/WilsonKernels.cc b/lib/qcd/action/fermion/WilsonKernels.cc index 8d5406f4..62ae93fa 100644 --- a/lib/qcd/action/fermion/WilsonKernels.cc +++ b/lib/qcd/action/fermion/WilsonKernels.cc @@ -327,7 +327,7 @@ void WilsonKernels::SeqConservedCurrentInternal(const PropagatorField &q_i GaugeLinkField Umu(U._grid); Umu = PeekIndex(U, mu); Lattice> t(q_in._grid); - LatticeCoordinate(t, mu); + LatticeCoordinate(t, Tp); tmp = this->CovShiftForward(Umu, mu, q_in)*ph; tmp = where((t >= tmin) and (t <= tmax), tmp, 0.*tmp); From 60f11bfd72f2c74cfdb0b91eaa9f44d80dd9946c Mon Sep 17 00:00:00 2001 From: Lanny91 Date: Wed, 7 Jun 2017 12:34:47 +0100 Subject: [PATCH 28/50] Removed redundant test module --- extras/Hadrons/Modules.hpp | 1 - .../Modules/MContraction/WardIdentitySeq.hpp | 145 ------------------ extras/Hadrons/modules.inc | 1 - 3 files changed, 147 deletions(-) delete mode 100644 extras/Hadrons/Modules/MContraction/WardIdentitySeq.hpp diff --git a/extras/Hadrons/Modules.hpp b/extras/Hadrons/Modules.hpp index dd6a6010..53ec346c 100644 --- a/extras/Hadrons/Modules.hpp +++ b/extras/Hadrons/Modules.hpp @@ -33,7 +33,6 @@ See the full license in the file "LICENSE" in the top level distribution directo #include #include #include -#include #include #include #include diff --git a/extras/Hadrons/Modules/MContraction/WardIdentitySeq.hpp b/extras/Hadrons/Modules/MContraction/WardIdentitySeq.hpp deleted file mode 100644 index 31409925..00000000 --- a/extras/Hadrons/Modules/MContraction/WardIdentitySeq.hpp +++ /dev/null @@ -1,145 +0,0 @@ -/************************************************************************************* - -Grid physics library, www.github.com/paboyle/Grid - -Source file: extras/Hadrons/Modules/MContraction/WardIdentitySeq.hpp - -Copyright (C) 2017 - -Author: Andrew Lawson - -This program is free software; you can redistribute it and/or modify -it under the terms of the GNU General Public License as published by -the Free Software Foundation; either version 2 of the License, or -(at your option) any later version. - -This program is distributed in the hope that it will be useful, -but WITHOUT ANY WARRANTY; without even the implied warranty of -MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the -GNU General Public License for more details. - -You should have received a copy of the GNU General Public License along -with this program; if not, write to the Free Software Foundation, Inc., -51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. - -See the full license in the file "LICENSE" in the top level distribution directory -*************************************************************************************/ -/* END LEGAL */ - -#ifndef Hadrons_WardIdentitySeq_hpp_ -#define Hadrons_WardIdentitySeq_hpp_ - -#include -#include -#include - -BEGIN_HADRONS_NAMESPACE - -/* - Ward Identity contractions using sequential propagators. - ----------------------------- - - * options: - - q_x: propagator, mu = x current insertion (string). - - q_y: propagator, mu = y current insertion (string). - - q_z: propagator, mu = z current insertion (string). - - q_t: propagator, mu = t current insertion (string). -*/ - -/****************************************************************************** - * WardIdentitySeq * - ******************************************************************************/ -BEGIN_MODULE_NAMESPACE(MContraction) - -class WardIdentitySeqPar: Serializable -{ -public: - GRID_SERIALIZABLE_CLASS_MEMBERS(WardIdentitySeqPar, - std::string, q_x, - std::string, q_y, - std::string, q_z, - std::string, q_t); -}; - -template -class TWardIdentitySeq: public Module -{ -public: - TYPE_ALIASES(FImpl,); -public: - // constructor - TWardIdentitySeq(const std::string name); - // destructor - virtual ~TWardIdentitySeq(void) = default; - // dependency relation - virtual std::vector getInput(void); - virtual std::vector getOutput(void); - // setup - virtual void setup(void); - // execution - virtual void execute(void); -}; - -MODULE_REGISTER_NS(WardIdentitySeq, TWardIdentitySeq, MContraction); - -/****************************************************************************** - * TWardIdentitySeq implementation * - ******************************************************************************/ -// constructor ///////////////////////////////////////////////////////////////// -template -TWardIdentitySeq::TWardIdentitySeq(const std::string name) -: Module(name) -{} - -// dependencies/products /////////////////////////////////////////////////////// -template -std::vector TWardIdentitySeq::getInput(void) -{ - std::vector in = {par().q_x, par().q_y, par().q_z, par().q_t}; - - return in; -} - -template -std::vector TWardIdentitySeq::getOutput(void) -{ - std::vector out = {getName()}; - - return out; -} - -// setup /////////////////////////////////////////////////////////////////////// -template -void TWardIdentitySeq::setup(void) -{ - -} - -// execution /////////////////////////////////////////////////////////////////// -template -void TWardIdentitySeq::execute(void) -{ - LatticeComplex vector_WI(env().getGrid()), c(env().getGrid()); - PropagatorField q_x = *env().template getObject(par().q_x); - PropagatorField q_y = *env().template getObject(par().q_y); - PropagatorField q_z = *env().template getObject(par().q_z); - PropagatorField q_t = *env().template getObject(par().q_t); - PropagatorField *q[Nd] = {&q_x, &q_y, &q_z, &q_t}; - Gamma g5(Gamma::Algebra::Gamma5); - - // Check D_mu V_mu = 0 - for (unsigned int mu = 0; mu < Nd; ++mu) - { - c = trace(g5*(*q[mu])); - vector_WI += c - Cshift(c, mu, -1); - } - - LOG(Message) << "Ward Identity checks for sequential vector current " - << "insertion = " << norm2(vector_WI) << std::endl; -} - -END_MODULE_NAMESPACE - -END_HADRONS_NAMESPACE - -#endif // Hadrons_WardIdentitySeq_hpp_ diff --git a/extras/Hadrons/modules.inc b/extras/Hadrons/modules.inc index 0364502a..b57aa577 100644 --- a/extras/Hadrons/modules.inc +++ b/extras/Hadrons/modules.inc @@ -14,7 +14,6 @@ modules_hpp =\ Modules/MContraction/Gamma3pt.hpp \ Modules/MContraction/Meson.hpp \ Modules/MContraction/WardIdentity.hpp \ - Modules/MContraction/WardIdentitySeq.hpp \ Modules/MContraction/WeakHamiltonian.hpp \ Modules/MContraction/WeakHamiltonianEye.hpp \ Modules/MContraction/WeakHamiltonianNonEye.hpp \ From b8e45ae490729a9ed79983974e1eeec1778a1e8d Mon Sep 17 00:00:00 2001 From: Lanny91 Date: Wed, 7 Jun 2017 16:26:22 +0100 Subject: [PATCH 29/50] Fixed remaining fermion type aliases after merge. --- extras/Hadrons/Modules/MContraction/WardIdentity.hpp | 2 +- extras/Hadrons/Modules/MSource/SeqConserved.hpp | 2 +- extras/Hadrons/Modules/MUtilities/TestSeqConserved.hpp | 2 +- extras/Hadrons/Modules/MUtilities/TestSeqGamma.hpp | 2 +- 4 files changed, 4 insertions(+), 4 deletions(-) diff --git a/extras/Hadrons/Modules/MContraction/WardIdentity.hpp b/extras/Hadrons/Modules/MContraction/WardIdentity.hpp index fb2ea173..82b0317a 100644 --- a/extras/Hadrons/Modules/MContraction/WardIdentity.hpp +++ b/extras/Hadrons/Modules/MContraction/WardIdentity.hpp @@ -65,7 +65,7 @@ template class TWardIdentity: public Module { public: - TYPE_ALIASES(FImpl,); + FERM_TYPE_ALIASES(FImpl,); public: // constructor TWardIdentity(const std::string name); diff --git a/extras/Hadrons/Modules/MSource/SeqConserved.hpp b/extras/Hadrons/Modules/MSource/SeqConserved.hpp index 6e5fb197..67086f11 100644 --- a/extras/Hadrons/Modules/MSource/SeqConserved.hpp +++ b/extras/Hadrons/Modules/MSource/SeqConserved.hpp @@ -74,7 +74,7 @@ template class TSeqConserved: public Module { public: - TYPE_ALIASES(FImpl,); + FERM_TYPE_ALIASES(FImpl,); public: // constructor TSeqConserved(const std::string name); diff --git a/extras/Hadrons/Modules/MUtilities/TestSeqConserved.hpp b/extras/Hadrons/Modules/MUtilities/TestSeqConserved.hpp index eccb00cc..faebab0a 100644 --- a/extras/Hadrons/Modules/MUtilities/TestSeqConserved.hpp +++ b/extras/Hadrons/Modules/MUtilities/TestSeqConserved.hpp @@ -71,7 +71,7 @@ template class TTestSeqConserved: public Module { public: - TYPE_ALIASES(FImpl,); + FERM_TYPE_ALIASES(FImpl,); public: // constructor TTestSeqConserved(const std::string name); diff --git a/extras/Hadrons/Modules/MUtilities/TestSeqGamma.hpp b/extras/Hadrons/Modules/MUtilities/TestSeqGamma.hpp index b3e99617..1b057c29 100644 --- a/extras/Hadrons/Modules/MUtilities/TestSeqGamma.hpp +++ b/extras/Hadrons/Modules/MUtilities/TestSeqGamma.hpp @@ -27,7 +27,7 @@ template class TTestSeqGamma: public Module { public: - TYPE_ALIASES(FImpl,); + FERM_TYPE_ALIASES(FImpl,); public: // constructor TTestSeqGamma(const std::string name); From 2d433ba30720f621b4d0a1bae91434aa4a42fe36 Mon Sep 17 00:00:00 2001 From: Lanny91 Date: Mon, 12 Jun 2017 10:32:14 +0100 Subject: [PATCH 30/50] Changed header include guards to match new convention --- extras/Hadrons/Modules/MContraction/WardIdentity.hpp | 4 ++-- extras/Hadrons/Modules/MSource/SeqConserved.hpp | 4 ++-- extras/Hadrons/Modules/MUtilities/TestSeqConserved.hpp | 4 ++-- extras/Hadrons/Modules/MUtilities/TestSeqGamma.hpp | 4 ++-- 4 files changed, 8 insertions(+), 8 deletions(-) diff --git a/extras/Hadrons/Modules/MContraction/WardIdentity.hpp b/extras/Hadrons/Modules/MContraction/WardIdentity.hpp index 82b0317a..8a56e0eb 100644 --- a/extras/Hadrons/Modules/MContraction/WardIdentity.hpp +++ b/extras/Hadrons/Modules/MContraction/WardIdentity.hpp @@ -26,8 +26,8 @@ See the full license in the file "LICENSE" in the top level distribution directo *************************************************************************************/ /* END LEGAL */ -#ifndef Hadrons_WardIdentity_hpp_ -#define Hadrons_WardIdentity_hpp_ +#ifndef Hadrons_MContraction_WardIdentity_hpp_ +#define Hadrons_MContraction_WardIdentity_hpp_ #include #include diff --git a/extras/Hadrons/Modules/MSource/SeqConserved.hpp b/extras/Hadrons/Modules/MSource/SeqConserved.hpp index 67086f11..86a7dfb9 100644 --- a/extras/Hadrons/Modules/MSource/SeqConserved.hpp +++ b/extras/Hadrons/Modules/MSource/SeqConserved.hpp @@ -26,8 +26,8 @@ See the full license in the file "LICENSE" in the top level distribution directo *************************************************************************************/ /* END LEGAL */ -#ifndef Hadrons_SeqConserved_hpp_ -#define Hadrons_SeqConserved_hpp_ +#ifndef Hadrons_MSource_SeqConserved_hpp_ +#define Hadrons_MSource_SeqConserved_hpp_ #include #include diff --git a/extras/Hadrons/Modules/MUtilities/TestSeqConserved.hpp b/extras/Hadrons/Modules/MUtilities/TestSeqConserved.hpp index faebab0a..b085eb8c 100644 --- a/extras/Hadrons/Modules/MUtilities/TestSeqConserved.hpp +++ b/extras/Hadrons/Modules/MUtilities/TestSeqConserved.hpp @@ -26,8 +26,8 @@ See the full license in the file "LICENSE" in the top level distribution directo *************************************************************************************/ /* END LEGAL */ -#ifndef Hadrons_TestSeqConserved_hpp_ -#define Hadrons_TestSeqConserved_hpp_ +#ifndef Hadrons_MUtilities_TestSeqConserved_hpp_ +#define Hadrons_MUtilities_TestSeqConserved_hpp_ #include #include diff --git a/extras/Hadrons/Modules/MUtilities/TestSeqGamma.hpp b/extras/Hadrons/Modules/MUtilities/TestSeqGamma.hpp index 1b057c29..3dbd7d63 100644 --- a/extras/Hadrons/Modules/MUtilities/TestSeqGamma.hpp +++ b/extras/Hadrons/Modules/MUtilities/TestSeqGamma.hpp @@ -1,5 +1,5 @@ -#ifndef Hadrons_TestSeqGamma_hpp_ -#define Hadrons_TestSeqGamma_hpp_ +#ifndef Hadrons_MUtilities_TestSeqGamma_hpp_ +#define Hadrons_MUtilities_TestSeqGamma_hpp_ #include #include From 5633a2db20e99cec2b5f11906632beb20eaadb31 Mon Sep 17 00:00:00 2001 From: Lanny91 Date: Mon, 12 Jun 2017 10:41:02 +0100 Subject: [PATCH 31/50] Faster implementation of conserved current site contraction. Added 5D vectorised support, but not G-parity. --- lib/qcd/action/fermion/FermionOperatorImpl.h | 29 ++++- lib/qcd/action/fermion/WilsonFermion.cc | 20 ++- lib/qcd/action/fermion/WilsonFermion5D.cc | 48 ++++--- lib/qcd/action/fermion/WilsonKernels.cc | 124 +++++++++++++++---- lib/qcd/action/fermion/WilsonKernels.h | 26 +++- 5 files changed, 198 insertions(+), 49 deletions(-) diff --git a/lib/qcd/action/fermion/FermionOperatorImpl.h b/lib/qcd/action/fermion/FermionOperatorImpl.h index 20458b6d..f330fb0d 100644 --- a/lib/qcd/action/fermion/FermionOperatorImpl.h +++ b/lib/qcd/action/fermion/FermionOperatorImpl.h @@ -212,6 +212,13 @@ namespace QCD { StencilImpl &St) { mult(&phi(), &U(mu), &chi()); } + + inline void multLinkProp(SitePropagator &phi, + const SiteDoubledGaugeField &U, + const SitePropagator &chi, + int mu) { + mult(&phi(), &U(mu), &chi()); + } template inline void loadLinkElement(Simd ®, ref &memory) { @@ -340,7 +347,20 @@ class DomainWallVec5dImpl : public PeriodicGaugeImpl< GaugeImplTypes< S,Nrepres } mult(&phi(), &UU(), &chi()); } - + + inline void multLinkProp(SitePropagator &phi, + const SiteDoubledGaugeField &U, + const SitePropagator &chi, + int mu) { + SiteGaugeLink UU; + for (int i = 0; i < Nrepresentation; i++) { + for (int j = 0; j < Nrepresentation; j++) { + vsplat(UU()()(i, j), U(mu)()(i, j)); + } + } + mult(&phi(), &UU(), &chi()); + } + inline void DoubleStore(GridBase *GaugeGrid, DoubledGaugeField &Uds,const GaugeField &Umu) { SiteScalarGaugeField ScalarUmu; @@ -538,6 +558,13 @@ class GparityWilsonImpl : public ConjugateGaugeImpl::ContractConservedCurrent(PropagatorField &q_in_1, conformable(_grid, q_in_1._grid); conformable(_grid, q_in_2._grid); conformable(_grid, q_out._grid); - Kernels::ContractConservedCurrentInternal(q_in_1, q_in_2, q_out, - Umu, curr_type, mu); + PropagatorField tmp(_grid); + q_out = zero; + + // Forward, need q1(x + mu), q2(x) + tmp = Cshift(q_in_1, mu, 1); + parallel_for (unsigned int sU = 0; sU < Umu._grid->oSites(); ++sU) + { + Kernels::ContractConservedCurrentSiteFwd(tmp, q_in_2, q_out, Umu, + mu, sU, sU, sU, sU); + } + + // Backward, need q1(x), q2(x + mu) + tmp = Cshift(q_in_2, mu, 1); + parallel_for (unsigned int sU = 0; sU < Umu._grid->oSites(); ++sU) + { + Kernels::ContractConservedCurrentSiteBwd(q_in_1, tmp, q_out, Umu, + mu, sU, sU, sU, sU); + } } template diff --git a/lib/qcd/action/fermion/WilsonFermion5D.cc b/lib/qcd/action/fermion/WilsonFermion5D.cc index 3bbc03b4..b69a18ba 100644 --- a/lib/qcd/action/fermion/WilsonFermion5D.cc +++ b/lib/qcd/action/fermion/WilsonFermion5D.cc @@ -687,26 +687,44 @@ void WilsonFermion5D::ContractConservedCurrent(PropagatorField &q_in_1, conformable(q_in_1._grid, q_in_2._grid); conformable(_FourDimGrid, q_out._grid); - PropagatorField q1_s(_FourDimGrid); - PropagatorField q2_s(_FourDimGrid); - PropagatorField tmp(_FourDimGrid); - - // Contract across 5th dimension. + PropagatorField tmp(FermionGrid()); q_out = zero; - for (int s = 0; s < Ls; ++s) - { - ExtractSlice(q1_s, q_in_1, s, 0); - ExtractSlice(q2_s, q_in_2, Ls - s - 1, 0); - Kernels::ContractConservedCurrentInternal(q1_s, q2_s, tmp, Umu, curr_type, mu); - // Axial current sign - if ((curr_type == Current::Axial) && (s < (Ls / 2))) + // Forward, need q1(x + mu, s), q2(x, Ls - 1 - s). 5D lattice so shift + // 4D coordinate mu by one. + tmp = Cshift(q_in_1, mu + 1, 1); + parallel_for (unsigned int sU = 0; sU < Umu._grid->oSites(); ++sU) + { + unsigned int sF1 = sU * Ls; + unsigned int sF2 = (sU + 1) * Ls - 1; + for (int s = 0; s < Ls; ++s) { - q_out -= tmp; + bool axial_sign = ((curr_type == Current::Axial) && (s < (Ls / 2))) ? \ + true : false; + Kernels::ContractConservedCurrentSiteFwd(tmp, q_in_2, q_out, Umu, + mu, sF1, sF2, sU, sU, + axial_sign); + sF1++; + sF2--; } - else + } + + // Backward, need q1(x, s), q2(x + mu, Ls - 1 - s). 5D lattice so shift + // 4D coordinate mu by one. + tmp = Cshift(q_in_2, mu + 1, 1); + parallel_for (unsigned int sU = 0; sU < Umu._grid->oSites(); ++sU) + { + unsigned int sF1 = sU * Ls; + unsigned int sF2 = (sU + 1) * Ls - 1; + for (int s = 0; s < Ls; ++s) { - q_out += tmp; + bool axial_sign = ((curr_type == Current::Axial) && (s < (Ls / 2))) ? \ + true : false; + Kernels::ContractConservedCurrentSiteBwd(q_in_1, tmp, q_out, Umu, + mu, sF1, sF2, sU, sU, + axial_sign); + sF1++; + sF2--; } } } diff --git a/lib/qcd/action/fermion/WilsonKernels.cc b/lib/qcd/action/fermion/WilsonKernels.cc index 62ae93fa..c519dc56 100644 --- a/lib/qcd/action/fermion/WilsonKernels.cc +++ b/lib/qcd/action/fermion/WilsonKernels.cc @@ -290,26 +290,110 @@ void WilsonKernels::DhopDir( StencilImpl &st, DoubledGaugeField &U,SiteHal #define WilsonCurrentFwd(expr, mu) ((expr - Gamma::gmu[mu]*expr)) #define WilsonCurrentBwd(expr, mu) ((expr + Gamma::gmu[mu]*expr)) +/******************************************************************************* + * Name: ContractConservedCurrentSiteFwd + * Operation: (1/2) * q2[x] * U(x) * (g[mu] - 1) * q1[x + mu] + * Notes: - DoubledGaugeField U assumed to contain -1/2 factor. + * - Pass in q_in_1 shifted in +ve mu direction. + ******************************************************************************/ template -void WilsonKernels::ContractConservedCurrentInternal(const PropagatorField &q_in_1, - const PropagatorField &q_in_2, - PropagatorField &q_out, - DoubledGaugeField &U, - Current curr_type, - unsigned int mu) +void WilsonKernels::ContractConservedCurrentSiteFwd( + const PropagatorField &q_in_1, + const PropagatorField &q_in_2, + PropagatorField &q_out, + DoubledGaugeField &U, + unsigned int mu, + unsigned int sF_in_1, + unsigned int sF_in_2, + unsigned int sF_out, + unsigned int sU, + bool switch_sign) { + SitePropagator result, tmp; Gamma g5(Gamma::Algebra::Gamma5); - PropagatorField tmp(q_out._grid); - GaugeLinkField Umu(U._grid); - Umu = PeekIndex(U, mu); - - tmp = this->CovShiftForward(Umu, mu, q_in_1); - q_out = (g5*adj(q_in_2)*g5)*WilsonCurrentFwd(tmp, mu); - - tmp = this->CovShiftForward(Umu, mu, q_in_2); - q_out -= (g5*adj(tmp)*g5)*WilsonCurrentBwd(q_in_1, mu); + multLinkProp(tmp, U._odata[sU], q_in_1._odata[sF_in_1], mu); + result = g5 * adj(q_in_2._odata[sF_in_2]) * g5 * WilsonCurrentFwd(tmp, mu); + if (switch_sign) + { + q_out._odata[sF_out] -= result; + } + else + { + q_out._odata[sF_out] += result; + } } +/******************************************************************************* + * Name: ContractConservedCurrentSiteBwd + * Operation: (1/2) * q2[x + mu] * U^dag(x) * (g[mu] + 1) * q1[x] + * Notes: - DoubledGaugeField U assumed to contain -1/2 factor. + * - Pass in q_in_2 shifted in +ve mu direction. + ******************************************************************************/ +template +void WilsonKernels::ContractConservedCurrentSiteBwd( + const PropagatorField &q_in_1, + const PropagatorField &q_in_2, + PropagatorField &q_out, + DoubledGaugeField &U, + unsigned int mu, + unsigned int sF_in_1, + unsigned int sF_in_2, + unsigned int sF_out, + unsigned int sU, + bool switch_sign) +{ + SitePropagator result, tmp; + Gamma g5(Gamma::Algebra::Gamma5); + multLinkProp(tmp, U._odata[sU], q_in_1._odata[sF_in_1], mu + Nd); + result = g5 * adj(q_in_2._odata[sF_in_2]) * g5 * WilsonCurrentBwd(tmp, mu); + if (switch_sign) + { + q_out._odata[sF_out] += result; + } + else + { + q_out._odata[sF_out] -= result; + } +} + +// G-parity requires more specialised implementation. +#define NO_CURR_SITE(Impl) \ +template <> \ +void WilsonKernels::ContractConservedCurrentSiteFwd( \ + const PropagatorField &q_in_1, \ + const PropagatorField &q_in_2, \ + PropagatorField &q_out, \ + DoubledGaugeField &U, \ + unsigned int mu, \ + unsigned int sF_in_1, \ + unsigned int sF_in_2, \ + unsigned int sF_out, \ + unsigned int sU, \ + bool switch_sign) \ +{ \ + assert(0); \ +} \ +template <> \ +void WilsonKernels::ContractConservedCurrentSiteBwd( \ + const PropagatorField &q_in_1, \ + const PropagatorField &q_in_2, \ + PropagatorField &q_out, \ + DoubledGaugeField &U, \ + unsigned int mu, \ + unsigned int sF_in_1, \ + unsigned int sF_in_2, \ + unsigned int sF_out, \ + unsigned int sU, \ + bool switch_sign) \ +{ \ + assert(0); \ +} + +NO_CURR_SITE(GparityWilsonImplF); +NO_CURR_SITE(GparityWilsonImplD); +NO_CURR_SITE(GparityWilsonImplFH); +NO_CURR_SITE(GparityWilsonImplDF); + template void WilsonKernels::SeqConservedCurrentInternal(const PropagatorField &q_in, @@ -342,16 +426,6 @@ void WilsonKernels::SeqConservedCurrentInternal(const PropagatorField &q_i // GParity, (Z)DomainWallVec5D -> require special implementation #define NO_CURR(Impl) \ -template <> void \ -WilsonKernels::ContractConservedCurrentInternal(const PropagatorField &q_in_1, \ - const PropagatorField &q_in_2, \ - PropagatorField &q_out, \ - DoubledGaugeField &U, \ - Current curr_type, \ - unsigned int mu) \ -{ \ - assert(0); \ -} \ template <> void \ WilsonKernels::SeqConservedCurrentInternal(const PropagatorField &q_in, \ PropagatorField &q_out, \ diff --git a/lib/qcd/action/fermion/WilsonKernels.h b/lib/qcd/action/fermion/WilsonKernels.h index 25c956ef..95155ccc 100644 --- a/lib/qcd/action/fermion/WilsonKernels.h +++ b/lib/qcd/action/fermion/WilsonKernels.h @@ -183,12 +183,26 @@ public: ////////////////////////////////////////////////////////////////////////////// // Utilities for inserting Wilson conserved current. ////////////////////////////////////////////////////////////////////////////// - void ContractConservedCurrentInternal(const PropagatorField &q_in_1, - const PropagatorField &q_in_2, - PropagatorField &q_out, - DoubledGaugeField &U, - Current curr_type, - unsigned int mu); + void ContractConservedCurrentSiteFwd(const PropagatorField &q_in_1, + const PropagatorField &q_in_2, + PropagatorField &q_out, + DoubledGaugeField &U, + unsigned int mu, + unsigned int sF_in_1, + unsigned int sF_in_2, + unsigned int sF_out, + unsigned int sU, + bool switch_sign = false); + void ContractConservedCurrentSiteBwd(const PropagatorField &q_in_1, + const PropagatorField &q_in_2, + PropagatorField &q_out, + DoubledGaugeField &U, + unsigned int mu, + unsigned int sF_in_1, + unsigned int sF_in_2, + unsigned int sF_out, + unsigned int sU, + bool switch_sign = false); void SeqConservedCurrentInternal(const PropagatorField &q_in, PropagatorField &q_out, DoubledGaugeField &U, From 2ad54c5a0222946571f133837f1a2cfaef807181 Mon Sep 17 00:00:00 2001 From: Lanny91 Date: Wed, 14 Jun 2017 10:53:39 +0100 Subject: [PATCH 32/50] QPX exchange support --- lib/simd/Grid_qpx.h | 43 +++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 43 insertions(+) diff --git a/lib/simd/Grid_qpx.h b/lib/simd/Grid_qpx.h index cbca9118..9fc8ef3c 100644 --- a/lib/simd/Grid_qpx.h +++ b/lib/simd/Grid_qpx.h @@ -375,6 +375,49 @@ namespace Optimization { FLOAT_WRAP_2(operator(), inline) }; + ////////////////////////////////////////////// + // Exchange support +#define FLOAT_WRAP_EXCHANGE(fn) \ + static inline void fn(vector4float &out1, vector4float &out2, \ + vector4float in1, vector4float in2) \ + { \ + vector4double out1d, out2d, in1d, in2d; \ + in1d = Vset()(in1); \ + in2d = Vset()(in2); \ + fn(out1d, out2d, in1d, in2d); \ + Vstore()(out1d, out1); \ + Vstore()(out2d, out2); \ + } + + struct Exchange{ + + // double precision + static inline void Exchange0(vector4double &out1, vector4double &out2, + vector4double in1, vector4double in2) { + out1 = vec_perm(in1, in2, vec_gpci(0145)); + out2 = vec_perm(in1, in2, vec_gpci(02367)); + } + static inline void Exchange1(vector4double &out1, vector4double &out2, + vector4double in1, vector4double in2) { + out1 = vec_perm(in1, in2, vec_gpci(0426)); + out2 = vec_perm(in1, in2, vec_gpci(01537)); + } + static inline void Exchange2(vector4double &out1, vector4double &out2, + vector4double in1, vector4double in2) { + assert(0); + } + static inline void Exchange3(vector4double &out1, vector4double &out2, + vector4double in1, vector4double in2) { + assert(0); + } + + // single precision + FLOAT_WRAP_EXCHANGE(Exchange0); + FLOAT_WRAP_EXCHANGE(Exchange1); + FLOAT_WRAP_EXCHANGE(Exchange2); + FLOAT_WRAP_EXCHANGE(Exchange3); + }; + struct Permute{ //Complex double static inline vector4double Permute0(vector4double v){ //0123 -> 2301 From 735cbdb983703fd3ffadc6133d792b4d058a897b Mon Sep 17 00:00:00 2001 From: Lanny91 Date: Wed, 14 Jun 2017 10:55:10 +0100 Subject: [PATCH 33/50] QPX Integer reduction (+ integer reduction test) --- lib/simd/Grid_qpx.h | 11 +++++++---- tests/Test_simd.cc | 47 +++++++++++++++++++++++++++++++++++++++++++-- 2 files changed, 52 insertions(+), 6 deletions(-) diff --git a/lib/simd/Grid_qpx.h b/lib/simd/Grid_qpx.h index 9fc8ef3c..00dbace5 100644 --- a/lib/simd/Grid_qpx.h +++ b/lib/simd/Grid_qpx.h @@ -540,10 +540,13 @@ namespace Optimization { //Integer Reduce template<> - inline Integer Reduce::operator()(int in){ - // FIXME unimplemented - printf("Reduce : Missing integer implementation -> FIX\n"); - assert(0); + inline Integer Reduce::operator()(veci in){ + Integer a = 0; + for (unsigned int i = 0; i < W::r; ++i) + { + a += in.v[i]; + } + return a; } } diff --git a/tests/Test_simd.cc b/tests/Test_simd.cc index c0bbef1d..b2e8d68e 100644 --- a/tests/Test_simd.cc +++ b/tests/Test_simd.cc @@ -183,8 +183,6 @@ void IntTester(const functor &func) { typedef Integer scal; typedef vInteger vec; - GridSerialRNG sRNG; - sRNG.SeedFixedIntegers(std::vector({45,12,81,9})); int Nsimd = vec::Nsimd(); @@ -287,6 +285,50 @@ void ReductionTester(const functor &func) } +template +void IntReductionTester(const functor &func) +{ + int Nsimd = vec::Nsimd(); + + std::vector input1(Nsimd); + std::vector input2(Nsimd); + reduced result(0); + reduced reference(0); + reduced tmp; + + std::vector > buf(3); + vec & v_input1 = buf[0]; + vec & v_input2 = buf[1]; + + for(int i=0;i(v_input1,input1); + merge(v_input2,input2); + + func.template vfunc(result,v_input1,v_input2); + + for(int i=0;i(tmp,input1[i],input2[i]); + reference+=tmp; + } + + std::cout<(funcReduce()); std::cout< Date: Fri, 16 Jun 2017 15:04:26 +0100 Subject: [PATCH 34/50] Placeholder precision change functions to allow Grid to compile with QPX (warning: no actual functionality) --- lib/simd/Grid_qpx.h | 36 ++++++++++++++++++++++++++++++++++++ 1 file changed, 36 insertions(+) diff --git a/lib/simd/Grid_qpx.h b/lib/simd/Grid_qpx.h index 00dbace5..8de7bde8 100644 --- a/lib/simd/Grid_qpx.h +++ b/lib/simd/Grid_qpx.h @@ -374,6 +374,41 @@ namespace Optimization { // Complex float FLOAT_WRAP_2(operator(), inline) }; +#define USE_FP16 + struct PrecisionChange { + static inline vech StoH (const vector4float &a, const vector4float &b) { + vech ret; + std::cout << GridLogError << "QPX single to half precision conversion not yet supported." << std::endl; + assert(0); + return ret; + } + static inline void HtoS (vech h, vector4float &sa, vector4float &sb) { + std::cout << GridLogError << "QPX half to single precision conversion not yet supported." << std::endl; + assert(0); + } + static inline vector4float DtoS (vector4double a, vector4double b) { + vector4float ret; + std::cout << GridLogError << "QPX double to single precision conversion not yet supported." << std::endl; + assert(0); + return ret; + } + static inline void StoD (vector4float s, vector4double &a, vector4double &b) { + std::cout << GridLogError << "QPX single to double precision conversion not yet supported." << std::endl; + assert(0); + } + static inline vech DtoH (vector4double a, vector4double b, + vector4double c, vector4double d) { + vech ret; + std::cout << GridLogError << "QPX double to half precision conversion not yet supported." << std::endl; + assert(0); + return ret; + } + static inline void HtoD (vech h, vector4double &a, vector4double &b, + vector4double &c, vector4double &d) { + std::cout << GridLogError << "QPX half to double precision conversion not yet supported." << std::endl; + assert(0); + } + }; ////////////////////////////////////////////// // Exchange support @@ -552,6 +587,7 @@ namespace Optimization { //////////////////////////////////////////////////////////////////////////////// // Here assign types +typedef Optimization::vech SIMD_Htype; // Half precision type typedef Optimization::vector4float SIMD_Ftype; // Single precision type typedef vector4double SIMD_Dtype; // Double precision type typedef Optimization::veci SIMD_Itype; // Integer type From a833f88c3237f9c941e9eb79ad459d0e260d2a2b Mon Sep 17 00:00:00 2001 From: Lanny91 Date: Fri, 16 Jun 2017 15:58:47 +0100 Subject: [PATCH 35/50] Added missing SIMD integer reduction implementation for AVX, AVX-512, SSE4, IMCI --- lib/simd/Grid_avx.h | 25 ++++++++++++++++++++++--- lib/simd/Grid_avx512.h | 22 +++++++++++++++++++--- lib/simd/Grid_imci.h | 4 +--- lib/simd/Grid_sse4.h | 6 +++--- 4 files changed, 45 insertions(+), 12 deletions(-) diff --git a/lib/simd/Grid_avx.h b/lib/simd/Grid_avx.h index 52be9c05..57d9064d 100644 --- a/lib/simd/Grid_avx.h +++ b/lib/simd/Grid_avx.h @@ -701,9 +701,28 @@ namespace Optimization { //Integer Reduce template<> inline Integer Reduce::operator()(__m256i in){ - // FIXME unimplemented - printf("Reduce : Missing integer implementation -> FIX\n"); - assert(0); + __m128i ret; +#if defined (AVX2) + // AVX2 horizontal adds within upper and lower halves of register; use + // SSE to add upper and lower halves for result. + __m256i v1, v2; + __m128i u1, u2; + v1 = _mm256_hadd_epi32(in, in); + v2 = _mm256_hadd_epi32(v1, v1); + u1 = _mm256_castsi256_si128(v2); // upper half + u2 = _mm256_extracti128_si256(v2, 1); // lower half + ret = _mm256_add_epi32(u1, u2); +#else + // No AVX horizontal add; extract upper and lower halves of register & use + // SSE intrinsics. + __m128i u1, u2, u3; + u1 = _mm256_extractf128_si256(in, 0); // upper half + u2 = _mm256_extractf128_si256(in, 1); // lower half + u3 = _mm_add_epi32(u1, u2); + u1 = _mm_hadd_epi32(u3, u3); + ret = _mm_hadd_epi32(u1, u1); +#endif + return _mm_cvtsi128_si32(ret); } } diff --git a/lib/simd/Grid_avx512.h b/lib/simd/Grid_avx512.h index ba054665..458a8f7c 100644 --- a/lib/simd/Grid_avx512.h +++ b/lib/simd/Grid_avx512.h @@ -543,6 +543,24 @@ namespace Optimization { u512d conv; conv.v = v1; return conv.f[0]; } + + //Integer Reduce + template<> + inline Integer Reduce::operator()(__m512i in){ + // No full vector reduce, use AVX to add upper and lower halves of register + // and perform AVX reduction. + __m256i v1, v2, v3; + __m128i u1, u2, ret; + v1 = _mm512_castsi512_si256(in); // upper half + v2 = _mm512_extracti32x8_epi32(in, 1); // lower half + v3 = _mm256_add_epi32(v1, v2); + v1 = _mm256_hadd_epi32(v3, v3); + v2 = _mm256_hadd_epi32(v1, v1); + u1 = _mm256_castsi256_si128(v2) // upper half + u2 = _mm256_extracti128_si256(v2, 1); // lower half + ret = _mm256_add_epi32(u1, u2); + return _mm_cvtsi128_si32(ret); + } #else //Complex float Reduce template<> @@ -570,9 +588,7 @@ namespace Optimization { //Integer Reduce template<> inline Integer Reduce::operator()(__m512i in){ - // FIXME unimplemented - printf("Reduce : Missing integer implementation -> FIX\n"); - assert(0); + return _mm512_reduce_add_epi32(in); } #endif diff --git a/lib/simd/Grid_imci.h b/lib/simd/Grid_imci.h index 173e57d8..a1dae565 100644 --- a/lib/simd/Grid_imci.h +++ b/lib/simd/Grid_imci.h @@ -401,9 +401,7 @@ namespace Optimization { //Integer Reduce template<> inline Integer Reduce::operator()(__m512i in){ - // FIXME unimplemented - printf("Reduce : Missing integer implementation -> FIX\n"); - assert(0); + return _mm512_reduce_add_epi32(in); } diff --git a/lib/simd/Grid_sse4.h b/lib/simd/Grid_sse4.h index 2fb2df76..0b1f9ffb 100644 --- a/lib/simd/Grid_sse4.h +++ b/lib/simd/Grid_sse4.h @@ -570,9 +570,9 @@ namespace Optimization { //Integer Reduce template<> inline Integer Reduce::operator()(__m128i in){ - // FIXME unimplemented - printf("Reduce : Missing integer implementation -> FIX\n"); - assert(0); + __m128i v1 = _mm_hadd_epi32(in, in); + __m128i v2 = _mm_hadd_epi32(v1, v1); + return _mm_cvtsi128_si32(v2); } } From 41af8c12d70145320a1f2fd924464802f26cffff Mon Sep 17 00:00:00 2001 From: Lanny91 Date: Fri, 16 Jun 2017 16:38:59 +0100 Subject: [PATCH 36/50] Code cleaning for conserved current contractions. Will now be easier to implement mobius conserved current. --- lib/qcd/action/fermion/WilsonFermion.cc | 26 ++++---- lib/qcd/action/fermion/WilsonFermion5D.cc | 40 ++++--------- lib/qcd/action/fermion/WilsonKernels.cc | 72 ++++++++++------------- lib/qcd/action/fermion/WilsonKernels.h | 22 +++---- 4 files changed, 64 insertions(+), 96 deletions(-) diff --git a/lib/qcd/action/fermion/WilsonFermion.cc b/lib/qcd/action/fermion/WilsonFermion.cc index 28842cdd..eff7d958 100644 --- a/lib/qcd/action/fermion/WilsonFermion.cc +++ b/lib/qcd/action/fermion/WilsonFermion.cc @@ -361,23 +361,23 @@ void WilsonFermion::ContractConservedCurrent(PropagatorField &q_in_1, conformable(_grid, q_in_1._grid); conformable(_grid, q_in_2._grid); conformable(_grid, q_out._grid); - PropagatorField tmp(_grid); + PropagatorField tmp1(_grid), tmp2(_grid); q_out = zero; - // Forward, need q1(x + mu), q2(x) - tmp = Cshift(q_in_1, mu, 1); + // Forward, need q1(x + mu), q2(x). Backward, need q1(x), q2(x + mu). + // Inefficient comms method but not performance critical. + tmp1 = Cshift(q_in_1, mu, 1); + tmp2 = Cshift(q_in_2, mu, 1); parallel_for (unsigned int sU = 0; sU < Umu._grid->oSites(); ++sU) { - Kernels::ContractConservedCurrentSiteFwd(tmp, q_in_2, q_out, Umu, - mu, sU, sU, sU, sU); - } - - // Backward, need q1(x), q2(x + mu) - tmp = Cshift(q_in_2, mu, 1); - parallel_for (unsigned int sU = 0; sU < Umu._grid->oSites(); ++sU) - { - Kernels::ContractConservedCurrentSiteBwd(q_in_1, tmp, q_out, Umu, - mu, sU, sU, sU, sU); + Kernels::ContractConservedCurrentSiteFwd(tmp1._odata[sU], + q_in_2._odata[sU], + q_out._odata[sU], + Umu, sU, mu); + Kernels::ContractConservedCurrentSiteBwd(q_in_1._odata[sU], + tmp2._odata[sU], + q_out._odata[sU], + Umu, sU, mu); } } diff --git a/lib/qcd/action/fermion/WilsonFermion5D.cc b/lib/qcd/action/fermion/WilsonFermion5D.cc index b69a18ba..76218098 100644 --- a/lib/qcd/action/fermion/WilsonFermion5D.cc +++ b/lib/qcd/action/fermion/WilsonFermion5D.cc @@ -686,13 +686,13 @@ void WilsonFermion5D::ContractConservedCurrent(PropagatorField &q_in_1, conformable(q_in_1._grid, FermionGrid()); conformable(q_in_1._grid, q_in_2._grid); conformable(_FourDimGrid, q_out._grid); - - PropagatorField tmp(FermionGrid()); + PropagatorField tmp1(FermionGrid()), tmp2(FermionGrid()); q_out = zero; - // Forward, need q1(x + mu, s), q2(x, Ls - 1 - s). 5D lattice so shift - // 4D coordinate mu by one. - tmp = Cshift(q_in_1, mu + 1, 1); + // Forward, need q1(x + mu, s), q2(x, Ls - 1 - s). Backward, need q1(x, s), + // q2(x + mu, Ls - 1 - s). 5D lattice so shift 4D coordinate mu by one. + tmp1 = Cshift(q_in_1, mu + 1, 1); + tmp2 = Cshift(q_in_2, mu + 1, 1); parallel_for (unsigned int sU = 0; sU < Umu._grid->oSites(); ++sU) { unsigned int sF1 = sU * Ls; @@ -701,28 +701,14 @@ void WilsonFermion5D::ContractConservedCurrent(PropagatorField &q_in_1, { bool axial_sign = ((curr_type == Current::Axial) && (s < (Ls / 2))) ? \ true : false; - Kernels::ContractConservedCurrentSiteFwd(tmp, q_in_2, q_out, Umu, - mu, sF1, sF2, sU, sU, - axial_sign); - sF1++; - sF2--; - } - } - - // Backward, need q1(x, s), q2(x + mu, Ls - 1 - s). 5D lattice so shift - // 4D coordinate mu by one. - tmp = Cshift(q_in_2, mu + 1, 1); - parallel_for (unsigned int sU = 0; sU < Umu._grid->oSites(); ++sU) - { - unsigned int sF1 = sU * Ls; - unsigned int sF2 = (sU + 1) * Ls - 1; - for (int s = 0; s < Ls; ++s) - { - bool axial_sign = ((curr_type == Current::Axial) && (s < (Ls / 2))) ? \ - true : false; - Kernels::ContractConservedCurrentSiteBwd(q_in_1, tmp, q_out, Umu, - mu, sF1, sF2, sU, sU, - axial_sign); + Kernels::ContractConservedCurrentSiteFwd(tmp1._odata[sF1], + q_in_2._odata[sF2], + q_out._odata[sU], + Umu, sU, mu, axial_sign); + Kernels::ContractConservedCurrentSiteBwd(q_in_1._odata[sF1], + tmp2._odata[sF2], + q_out._odata[sU], + Umu, sU, mu, axial_sign); sF1++; sF2--; } diff --git a/lib/qcd/action/fermion/WilsonKernels.cc b/lib/qcd/action/fermion/WilsonKernels.cc index c519dc56..6b193766 100644 --- a/lib/qcd/action/fermion/WilsonKernels.cc +++ b/lib/qcd/action/fermion/WilsonKernels.cc @@ -298,28 +298,25 @@ void WilsonKernels::DhopDir( StencilImpl &st, DoubledGaugeField &U,SiteHal ******************************************************************************/ template void WilsonKernels::ContractConservedCurrentSiteFwd( - const PropagatorField &q_in_1, - const PropagatorField &q_in_2, - PropagatorField &q_out, + const SitePropagator &q_in_1, + const SitePropagator &q_in_2, + SitePropagator &q_out, DoubledGaugeField &U, - unsigned int mu, - unsigned int sF_in_1, - unsigned int sF_in_2, - unsigned int sF_out, unsigned int sU, + unsigned int mu, bool switch_sign) { SitePropagator result, tmp; Gamma g5(Gamma::Algebra::Gamma5); - multLinkProp(tmp, U._odata[sU], q_in_1._odata[sF_in_1], mu); - result = g5 * adj(q_in_2._odata[sF_in_2]) * g5 * WilsonCurrentFwd(tmp, mu); + Impl::multLinkProp(tmp, U._odata[sU], q_in_1, mu); + result = g5 * adj(q_in_2) * g5 * WilsonCurrentFwd(tmp, mu); if (switch_sign) { - q_out._odata[sF_out] -= result; + q_out -= result; } else { - q_out._odata[sF_out] += result; + q_out += result; } } @@ -331,28 +328,25 @@ void WilsonKernels::ContractConservedCurrentSiteFwd( ******************************************************************************/ template void WilsonKernels::ContractConservedCurrentSiteBwd( - const PropagatorField &q_in_1, - const PropagatorField &q_in_2, - PropagatorField &q_out, + const SitePropagator &q_in_1, + const SitePropagator &q_in_2, + SitePropagator &q_out, DoubledGaugeField &U, - unsigned int mu, - unsigned int sF_in_1, - unsigned int sF_in_2, - unsigned int sF_out, unsigned int sU, + unsigned int mu, bool switch_sign) { SitePropagator result, tmp; Gamma g5(Gamma::Algebra::Gamma5); - multLinkProp(tmp, U._odata[sU], q_in_1._odata[sF_in_1], mu + Nd); - result = g5 * adj(q_in_2._odata[sF_in_2]) * g5 * WilsonCurrentBwd(tmp, mu); + Impl::multLinkProp(tmp, U._odata[sU], q_in_1, mu + Nd); + result = g5 * adj(q_in_2) * g5 * WilsonCurrentBwd(tmp, mu); if (switch_sign) { - q_out._odata[sF_out] += result; + q_out += result; } else { - q_out._odata[sF_out] -= result; + q_out -= result; } } @@ -360,31 +354,25 @@ void WilsonKernels::ContractConservedCurrentSiteBwd( #define NO_CURR_SITE(Impl) \ template <> \ void WilsonKernels::ContractConservedCurrentSiteFwd( \ - const PropagatorField &q_in_1, \ - const PropagatorField &q_in_2, \ - PropagatorField &q_out, \ - DoubledGaugeField &U, \ - unsigned int mu, \ - unsigned int sF_in_1, \ - unsigned int sF_in_2, \ - unsigned int sF_out, \ - unsigned int sU, \ - bool switch_sign) \ + const SitePropagator &q_in_1, \ + const SitePropagator &q_in_2, \ + SitePropagator &q_out, \ + DoubledGaugeField &U, \ + unsigned int sU, \ + unsigned int mu, \ + bool switch_sign) \ { \ assert(0); \ } \ template <> \ void WilsonKernels::ContractConservedCurrentSiteBwd( \ - const PropagatorField &q_in_1, \ - const PropagatorField &q_in_2, \ - PropagatorField &q_out, \ - DoubledGaugeField &U, \ - unsigned int mu, \ - unsigned int sF_in_1, \ - unsigned int sF_in_2, \ - unsigned int sF_out, \ - unsigned int sU, \ - bool switch_sign) \ + const SitePropagator &q_in_1, \ + const SitePropagator &q_in_2, \ + SitePropagator &q_out, \ + DoubledGaugeField &U, \ + unsigned int mu, \ + unsigned int sU, \ + bool switch_sign) \ { \ assert(0); \ } diff --git a/lib/qcd/action/fermion/WilsonKernels.h b/lib/qcd/action/fermion/WilsonKernels.h index 95155ccc..0294c740 100644 --- a/lib/qcd/action/fermion/WilsonKernels.h +++ b/lib/qcd/action/fermion/WilsonKernels.h @@ -183,25 +183,19 @@ public: ////////////////////////////////////////////////////////////////////////////// // Utilities for inserting Wilson conserved current. ////////////////////////////////////////////////////////////////////////////// - void ContractConservedCurrentSiteFwd(const PropagatorField &q_in_1, - const PropagatorField &q_in_2, - PropagatorField &q_out, + void ContractConservedCurrentSiteFwd(const SitePropagator &q_in_1, + const SitePropagator &q_in_2, + SitePropagator &q_out, DoubledGaugeField &U, - unsigned int mu, - unsigned int sF_in_1, - unsigned int sF_in_2, - unsigned int sF_out, unsigned int sU, + unsigned int mu, bool switch_sign = false); - void ContractConservedCurrentSiteBwd(const PropagatorField &q_in_1, - const PropagatorField &q_in_2, - PropagatorField &q_out, + void ContractConservedCurrentSiteBwd(const SitePropagator &q_in_1, + const SitePropagator &q_in_2, + SitePropagator &q_out, DoubledGaugeField &U, - unsigned int mu, - unsigned int sF_in_1, - unsigned int sF_in_2, - unsigned int sF_out, unsigned int sU, + unsigned int mu, bool switch_sign = false); void SeqConservedCurrentInternal(const PropagatorField &q_in, PropagatorField &q_out, From 1bd311ba9ccd8506d13064cb6f6829515a0f0240 Mon Sep 17 00:00:00 2001 From: Lanny91 Date: Fri, 16 Jun 2017 16:43:15 +0100 Subject: [PATCH 37/50] Faster sequential conserved current implementation, now compatible with 5D vectorisation & G-parity. --- lib/qcd/action/fermion/WilsonFermion.cc | 41 ++++++++- lib/qcd/action/fermion/WilsonFermion5D.cc | 68 +++++++++++--- lib/qcd/action/fermion/WilsonKernels.cc | 105 ++++++++++++---------- lib/qcd/action/fermion/WilsonKernels.h | 22 +++-- 4 files changed, 164 insertions(+), 72 deletions(-) diff --git a/lib/qcd/action/fermion/WilsonFermion.cc b/lib/qcd/action/fermion/WilsonFermion.cc index eff7d958..b986edd7 100644 --- a/lib/qcd/action/fermion/WilsonFermion.cc +++ b/lib/qcd/action/fermion/WilsonFermion.cc @@ -394,6 +394,8 @@ void WilsonFermion::SeqConservedCurrent(PropagatorField &q_in, conformable(_grid, q_out._grid); Lattice> ph(_grid), coor(_grid); Complex i(0.0,1.0); + PropagatorField tmpFwd(_grid), tmpBwd(_grid), tmp(_grid); + int tshift = (mu == Tp) ? 1 : 0; // Momentum projection ph = zero; @@ -404,8 +406,43 @@ void WilsonFermion::SeqConservedCurrent(PropagatorField &q_in, } ph = exp((Real)(2*M_PI)*i*ph); - Kernels::SeqConservedCurrentInternal(q_in, q_out, Umu, curr_type, mu, ph, - tmin, tmax); + q_out = zero; + LatticeInteger coords(_grid); + LatticeCoordinate(coords, Tp); + + // Need q(x + mu) and q(x - mu). + tmp = Cshift(q_in, mu, 1); + tmpFwd = tmp*ph; + tmp = ph*q_in; + tmpBwd = Cshift(tmp, mu, -1); + + parallel_for (unsigned int sU = 0; sU < Umu._grid->oSites(); ++sU) + { + // Compute the sequential conserved current insertion only if our simd + // object contains a timeslice we need. + vInteger t_mask = ((coords._odata[sU] >= tmin) && + (coords._odata[sU] <= tmax)); + Integer timeSlices = Reduce(t_mask); + + if (timeSlices > 0) + { + Kernels::SeqConservedCurrentSiteFwd(tmpFwd._odata[sU], + q_out._odata[sU], + Umu, sU, mu, t_mask); + } + + // Repeat for backward direction. + t_mask = ((coords._odata[sU] >= (tmin + tshift)) && + (coords._odata[sU] <= (tmax + tshift))); + timeSlices = Reduce(t_mask); + + if (timeSlices > 0) + { + Kernels::SeqConservedCurrentSiteBwd(tmpBwd._odata[sU], + q_out._odata[sU], + Umu, sU, mu, t_mask); + } + } } FermOpTemplateInstantiate(WilsonFermion); diff --git a/lib/qcd/action/fermion/WilsonFermion5D.cc b/lib/qcd/action/fermion/WilsonFermion5D.cc index 76218098..5daed3de 100644 --- a/lib/qcd/action/fermion/WilsonFermion5D.cc +++ b/lib/qcd/action/fermion/WilsonFermion5D.cc @@ -727,31 +727,73 @@ void WilsonFermion5D::SeqConservedCurrent(PropagatorField &q_in, { conformable(q_in._grid, FermionGrid()); conformable(q_in._grid, q_out._grid); - Lattice> ph(_FourDimGrid), coor(_FourDimGrid); - PropagatorField q_in_s(_FourDimGrid); - PropagatorField q_out_s(_FourDimGrid); + Lattice> ph(FermionGrid()), coor(FermionGrid()); + PropagatorField tmpFwd(FermionGrid()), tmpBwd(FermionGrid()), + tmp(FermionGrid()); Complex i(0.0, 1.0); + int tshift = (mu == Tp) ? 1 : 0; - // Momentum projection + // Momentum projection. ph = zero; for(unsigned int nu = 0; nu < Nd - 1; nu++) { - LatticeCoordinate(coor, nu); + // Shift coordinate lattice index by 1 to account for 5th dimension. + LatticeCoordinate(coor, nu + 1); ph = ph + mom[nu]*coor*((1./(_FourDimGrid->_fdimensions[nu]))); } ph = exp((Real)(2*M_PI)*i*ph); - // Sequential insertion across 5th dimension - for (int s = 0; s < Ls; s++) + q_out = zero; + LatticeInteger coords(_FourDimGrid); + LatticeCoordinate(coords, Tp); + + // Need q(x + mu, s) and q(x - mu, s). 5D lattice so shift 4D coordinate mu + // by one. + tmp = Cshift(q_in, mu + 1, 1); + tmpFwd = tmp*ph; + tmp = ph*q_in; + tmpBwd = Cshift(tmp, mu + 1, -1); + + parallel_for (unsigned int sU = 0; sU < Umu._grid->oSites(); ++sU) { - ExtractSlice(q_in_s, q_in, s, 0); - Kernels::SeqConservedCurrentInternal(q_in_s, q_out_s, Umu, curr_type, - mu, ph, tmin, tmax); - if ((curr_type == Current::Axial) && (s < Ls/2)) + // Compute the sequential conserved current insertion only if our simd + // object contains a timeslice we need. + vInteger t_mask = ((coords._odata[sU] >= tmin) && + (coords._odata[sU] <= tmax)); + Integer timeSlices = Reduce(t_mask); + + if (timeSlices > 0) { - q_out_s = -q_out_s; + unsigned int sF = sU * Ls; + for (unsigned int s = 0; s < Ls; ++s) + { + bool axial_sign = ((curr_type == Current::Axial) && (s < (Ls / 2))) ? \ + true : false; + Kernels::SeqConservedCurrentSiteFwd(tmpFwd._odata[sF], + q_out._odata[sF], Umu, sU, + mu, t_mask, axial_sign); + ++sF; + } + } + + // Repeat for backward direction. + t_mask = ((coords._odata[sU] >= (tmin + tshift)) && + (coords._odata[sU] <= (tmax + tshift))); + timeSlices = Reduce(t_mask); + + if (timeSlices > 0) + { + unsigned int sF = sU * Ls; + for (unsigned int s = 0; s < Ls; ++s) + { + bool axial_sign = ((curr_type == Current::Axial) && (s < (Ls / 2))) ? \ + true : false; + Kernels::SeqConservedCurrentSiteBwd(tmpBwd._odata[sF], + q_out._odata[sF], Umu, sU, + mu, t_mask, axial_sign); + ++sF; + } } - InsertSlice(q_out_s, q_out, s, 0); } } diff --git a/lib/qcd/action/fermion/WilsonKernels.cc b/lib/qcd/action/fermion/WilsonKernels.cc index 6b193766..dc66db23 100644 --- a/lib/qcd/action/fermion/WilsonKernels.cc +++ b/lib/qcd/action/fermion/WilsonKernels.cc @@ -383,63 +383,70 @@ NO_CURR_SITE(GparityWilsonImplFH); NO_CURR_SITE(GparityWilsonImplDF); -template -void WilsonKernels::SeqConservedCurrentInternal(const PropagatorField &q_in, - PropagatorField &q_out, - DoubledGaugeField &U, - Current curr_type, - unsigned int mu, - Lattice> &ph, - unsigned int tmin, - unsigned int tmax) +/******************************************************************************* + * Name: SeqConservedCurrentSiteFwd + * Operation: (1/2) * U(x) * (g[mu] - 1) * q[x + mu] + * Notes: - DoubledGaugeField U assumed to contain -1/2 factor. + * - Pass in q_in shifted in +ve mu direction. + ******************************************************************************/ +template +void WilsonKernels::SeqConservedCurrentSiteFwd(const SitePropagator &q_in, + SitePropagator &q_out, + DoubledGaugeField &U, + unsigned int sU, + unsigned int mu, + vInteger t_mask, + bool switch_sign) { - 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(U, mu); - Lattice> t(q_in._grid); - LatticeCoordinate(t, Tp); + SitePropagator result; + Impl::multLinkProp(result, U._odata[sU], q_in, mu); + result = WilsonCurrentFwd(result, mu); - tmp = this->CovShiftForward(Umu, mu, q_in)*ph; - tmp = where((t >= tmin) and (t <= tmax), tmp, 0.*tmp); - q_out = G_T*WilsonCurrentFwd(tmp, mu); + // Zero any unwanted timeslice entries. + result = predicatedWhere(t_mask, result, 0.*result); - tmp = q_in*ph; - tmp = this->CovShiftBackward(Umu, mu, tmp); - tmp = where((t >= tmin + tshift) and (t <= tmax + tshift), tmp, 0.*tmp); - q_out -= WilsonCurrentBwd(tmp, mu); + if (switch_sign) + { + q_out -= result; + } + else + { + q_out += result; + } } +/******************************************************************************* + * Name: SeqConservedCurrentSiteFwd + * Operation: (1/2) * U^dag(x) * (g[mu] + 1) * q[x - mu] + * Notes: - DoubledGaugeField U assumed to contain -1/2 factor. + * - Pass in q_in shifted in -ve mu direction. + ******************************************************************************/ +template +void WilsonKernels::SeqConservedCurrentSiteBwd(const SitePropagator &q_in, + SitePropagator &q_out, + DoubledGaugeField &U, + unsigned int sU, + unsigned int mu, + vInteger t_mask, + bool switch_sign) +{ + SitePropagator result; + Impl::multLinkProp(result, U._odata[sU], q_in, mu + Nd); + result = WilsonCurrentBwd(result, mu); -// GParity, (Z)DomainWallVec5D -> require special implementation -#define NO_CURR(Impl) \ -template <> void \ -WilsonKernels::SeqConservedCurrentInternal(const PropagatorField &q_in, \ - PropagatorField &q_out, \ - DoubledGaugeField &U, \ - Current curr_type, \ - unsigned int mu, \ - Lattice> &ph, \ - unsigned int tmin, \ - unsigned int tmax) \ -{ \ - assert(0); \ + // Zero any unwanted timeslice entries. + result = predicatedWhere(t_mask, result, 0.*result); + + if (switch_sign) + { + q_out += result; + } + else + { + q_out -= result; + } } -NO_CURR(GparityWilsonImplF); -NO_CURR(GparityWilsonImplD); -NO_CURR(GparityWilsonImplFH); -NO_CURR(GparityWilsonImplDF); -NO_CURR(DomainWallVec5dImplF); -NO_CURR(DomainWallVec5dImplD); -NO_CURR(DomainWallVec5dImplFH); -NO_CURR(DomainWallVec5dImplDF); -NO_CURR(ZDomainWallVec5dImplF); -NO_CURR(ZDomainWallVec5dImplD); -NO_CURR(ZDomainWallVec5dImplFH); -NO_CURR(ZDomainWallVec5dImplDF); - FermOpTemplateInstantiate(WilsonKernels); AdjointFermOpTemplateInstantiate(WilsonKernels); TwoIndexFermOpTemplateInstantiate(WilsonKernels); diff --git a/lib/qcd/action/fermion/WilsonKernels.h b/lib/qcd/action/fermion/WilsonKernels.h index 0294c740..ed8d6be9 100644 --- a/lib/qcd/action/fermion/WilsonKernels.h +++ b/lib/qcd/action/fermion/WilsonKernels.h @@ -197,14 +197,20 @@ public: unsigned int sU, unsigned int mu, bool switch_sign = false); - void SeqConservedCurrentInternal(const PropagatorField &q_in, - PropagatorField &q_out, - DoubledGaugeField &U, - Current curr_type, - unsigned int mu, - Lattice> &ph, - unsigned int tmin, - unsigned int tmax); + void SeqConservedCurrentSiteFwd(const SitePropagator &q_in, + SitePropagator &q_out, + DoubledGaugeField &U, + unsigned int sU, + unsigned int mu, + vInteger t_mask, + bool switch_sign = false); + void SeqConservedCurrentSiteBwd(const SitePropagator &q_in, + SitePropagator &q_out, + DoubledGaugeField &U, + unsigned int sU, + unsigned int mu, + vInteger t_mask, + bool switch_sign = false); private: // Specialised variants From 863bb2ad1007d7a00c8aedac932571d56d64a75f Mon Sep 17 00:00:00 2001 From: Lanny91 Date: Thu, 22 Jun 2017 16:02:15 +0200 Subject: [PATCH 38/50] Moving overly-specialised code out of Grid --- tests/hadrons/Test_hadrons_rarekaon.cc | 321 ------------------------- 1 file changed, 321 deletions(-) delete mode 100644 tests/hadrons/Test_hadrons_rarekaon.cc diff --git a/tests/hadrons/Test_hadrons_rarekaon.cc b/tests/hadrons/Test_hadrons_rarekaon.cc deleted file mode 100644 index a85beead..00000000 --- a/tests/hadrons/Test_hadrons_rarekaon.cc +++ /dev/null @@ -1,321 +0,0 @@ -/******************************************************************************* - Grid physics library, www.github.com/paboyle/Grid - - Source file: tests/hadrons/Test_hadrons_rarekaon.cc - - Copyright (C) 2017 - - Author: Andrew Lawson - - This program is free software; you can redistribute it and/or modify - it under the terms of the GNU General Public License as published by - the Free Software Foundation; either version 2 of the License, or - (at your option) any later version. - - This program is distributed in the hope that it will be useful, - but WITHOUT ANY WARRANTY; without even the implied warranty of - MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the - GNU General Public License for more details. - - You should have received a copy of the GNU General Public License along - with this program; if not, write to the Free Software Foundation, Inc., - 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. - - See the full license in the file "LICENSE" in the top level distribution - directory. - *******************************************************************************/ - -#include "Test_hadrons.hpp" - -using namespace Grid; -using namespace Hadrons; - -enum quarks -{ - light = 0, - strange = 1, - charm = 2 -}; - -int main(int argc, char *argv[]) -{ - // parse command line ////////////////////////////////////////////////////// - std::string configStem; - - if (argc < 2) - { - std::cerr << "usage: " << argv[0] << " [Grid options]"; - std::cerr << std::endl; - std::exit(EXIT_FAILURE); - } - configStem = argv[1]; - - // initialization ////////////////////////////////////////////////////////// - HADRONS_DEFAULT_INIT; - - // run setup /////////////////////////////////////////////////////////////// - Application application; - std::vector mass = {.01, .04, .2}; - std::vector flavour = {"l", "s", "c"}; - std::vector solvers = {"CG_l", "CG_s", "CG_c"}; - std::string kmom = "0. 0. 0. 0."; - std::string pmom = "1. 0. 0. 0."; - std::string qmom = "-1. 0. 0. 0."; - std::string mqmom = "1. 0. 0. 0."; - std::vector tKs = {0}; - unsigned int dt_pi = 16; - std::vector tJs = {8}; - unsigned int n_noise = 1; - unsigned int nt = 32; - bool do_disconnected(false); - Gamma::Algebra gT = Gamma::Algebra::GammaT; - unsigned int Ls = 16; - double M5 = 1.8; - - // Global parameters. - HADRONS_DEFAULT_GLOBALS(application); - - // gauge field - std::string gaugeField = "gauge"; - if (configStem == "None") - { - application.createModule(gaugeField); - } - else - { - MGauge::Load::Par gaugePar; - gaugePar.file = configStem; - application.createModule(gaugeField, gaugePar); - } - - // set fermion boundary conditions to be periodic space, antiperiodic time. - std::string boundary = "1 1 1 -1"; - - for (unsigned int i = 0; i < flavour.size(); ++i) - { - // actions - std::string actionName = "DWF_" + flavour[i]; - makeDWFAction(application, actionName, gaugeField, mass[i], M5, Ls); - - // solvers - makeRBPrecCGSolver(application, solvers[i], actionName); - } - - // Create noise propagators for loops. - std::vector noiseSrcs; - std::vector> noiseRes; - std::vector> noiseProps; - if (n_noise > 0) - { - MSource::Z2::Par noisePar; - noisePar.tA = 0; - noisePar.tB = nt - 1; - std::string loop_stem = "loop_"; - - noiseRes.resize(flavour.size()); - noiseProps.resize(flavour.size()); - for (unsigned int nn = 0; nn < n_noise; ++nn) - { - std::string eta = INIT_INDEX("noise", nn); - application.createModule(eta, noisePar); - noiseSrcs.push_back(eta); - - for (unsigned int f = 0; f < flavour.size(); ++f) - { - std::string loop_prop = INIT_INDEX(loop_stem + flavour[f], nn); - std::string loop_res = loop_prop + "_res"; - makePropagator(application, loop_res, eta, solvers[f]); - makeLoop(application, loop_prop, eta, loop_res); - noiseRes[f].push_back(loop_res); - noiseProps[f].push_back(loop_prop); - } - } - } - - // Translate rare kaon decay across specified timeslices. - for (unsigned int i = 0; i < tKs.size(); ++i) - { - // Zero-momentum wall source propagators for kaon and pion. - unsigned int tK = tKs[i]; - unsigned int tpi = (tK + dt_pi) % nt; - std::string q_Kl_0 = INIT_INDEX("Q_l_0", tK); - std::string q_pil_0 = INIT_INDEX("Q_l_0", tpi); - MAKE_WALL_PROP(tK, q_Kl_0, solvers[light]); - MAKE_WALL_PROP(tpi, q_pil_0, solvers[light]); - - // Wall sources for kaon and pion with momentum insertion. If either - // p or k are zero, or p = k, re-use the existing name to avoid - // duplicating a propagator. - std::string q_Ks_k = INIT_INDEX("Q_Ks_k", tK); - std::string q_Ks_p = INIT_INDEX((kmom == pmom) ? "Q_Ks_k" : "Q_Ks_p", tK); - std::string q_pil_k = INIT_INDEX((kmom == ZERO_MOM) ? "Q_l_0" : "Q_l_k", tpi); - std::string q_pil_p = INIT_INDEX((pmom == kmom) ? q_pil_k : ((pmom == ZERO_MOM) ? "Q_l_0" : "Q_l_p"), tpi); - MAKE_3MOM_WALL_PROP(tK, kmom, q_Ks_k, solvers[strange]); - MAKE_3MOM_WALL_PROP(tK, pmom, q_Ks_p, solvers[strange]); - MAKE_3MOM_WALL_PROP(tpi, kmom, q_pil_k, solvers[light]); - MAKE_3MOM_WALL_PROP(tpi, pmom, q_pil_p, solvers[light]); - - /*********************************************************************** - * CONTRACTIONS: pi and K 2pt contractions with mom = p, k. - **********************************************************************/ - // Wall-Point - std::string PW_K_k = INIT_INDEX("PW_K_k", tK); - std::string PW_K_p = INIT_INDEX("PW_K_p", tK); - std::string PW_pi_k = INIT_INDEX("PW_pi_k", tpi); - std::string PW_pi_p = INIT_INDEX("PW_pi_p", tpi); - mesonContraction(application, 2, q_Kl_0, q_Ks_k, PW_K_k, kmom); - mesonContraction(application, 2, q_Kl_0, q_Ks_p, PW_K_p, pmom); - mesonContraction(application, 2, q_pil_k, q_pil_0, PW_pi_k, kmom); - mesonContraction(application, 2, q_pil_p, q_pil_0, PW_pi_p, pmom); - // Wall-Wall, to be done - requires modification of meson module. - - /*********************************************************************** - * CONTRACTIONS: 3pt Weak Hamiltonian, C & W (non-Eye type) classes. - **********************************************************************/ - std::string HW_CW_k = LABEL_3PT("HW_CW_k", tK, tpi); - std::string HW_CW_p = LABEL_3PT("HW_CW_p", tK, tpi); - weakContractionNonEye(application, 3, q_Kl_0, q_Ks_k, q_pil_k, q_pil_0, HW_CW_k); - weakContractionNonEye(application, 3, q_Kl_0, q_Ks_p, q_pil_p, q_pil_0, HW_CW_p); - - /*********************************************************************** - * CONTRACTIONS: 3pt sd insertion. - **********************************************************************/ - // Note: eventually will use wall sink smeared q_Kl_0 instead. - std::string sd_k = LABEL_3PT("sd_k", tK, tpi); - std::string sd_p = LABEL_3PT("sd_p", tK, tpi); - gamma3ptContraction(application, 3, q_Kl_0, q_Ks_k, q_pil_k, sd_k); - gamma3ptContraction(application, 3, q_Kl_0, q_Ks_p, q_pil_p, sd_p); - - for (unsigned int nn = 0; nn < n_noise; ++nn) - { - /******************************************************************* - * CONTRACTIONS: 3pt Weak Hamiltonian, S and E (Eye type) classes. - ******************************************************************/ - // Note: eventually will use wall sink smeared q_Kl_0 instead. - for (unsigned int f = 0; f < flavour.size(); ++f) - { - if ((f != strange) || do_disconnected) - { - std::string HW_SE_k = LABEL_3PT("HW_SE_k_" + flavour[f], tK, tpi); - std::string HW_SE_p = LABEL_3PT("HW_SE_p_" + flavour[f], tK, tpi); - std::string loop_q = noiseProps[f][nn]; - weakContractionEye(application, 3, q_Kl_0, q_Ks_k, q_pil_k, loop_q, HW_CW_k); - weakContractionEye(application, 3, q_Kl_0, q_Ks_p, q_pil_p, loop_q, HW_CW_p); - } - } - } - - // Perform separate contractions for each t_J position. - for (unsigned int j = 0; j < tJs.size(); ++j) - { - // Sequential sources for current insertions. Local for now, - // gamma_0 only. - unsigned int tJ = (tJs[j] + tK) % nt; - MSource::SeqGamma::Par seqPar; - std::string q_KlCl_q = LABEL_3PT("Q_KlCl_q", tK, tJ); - std::string q_KsCs_mq = LABEL_3PT("Q_KsCs_mq", tK, tJ); - std::string q_pilCl_q = LABEL_3PT("Q_pilCl_q", tpi, tJ); - std::string q_pilCl_mq = LABEL_3PT("Q_pilCl_mq", tpi, tJ); - MAKE_SEQUENTIAL_PROP(tJ, q_Kl_0, qmom, q_KlCl_q, solvers[light], gT); - MAKE_SEQUENTIAL_PROP(tJ, q_Ks_k, mqmom, q_KsCs_mq, solvers[strange], gT); - MAKE_SEQUENTIAL_PROP(tJ, q_pil_p, qmom, q_pilCl_q, solvers[light], gT); - MAKE_SEQUENTIAL_PROP(tJ, q_pil_0, mqmom, q_pilCl_mq, solvers[light], gT); - - /******************************************************************* - * CONTRACTIONS: pi and K 3pt contractions with current insertion. - ******************************************************************/ - // Wall-Point - std::string C_PW_Kl = LABEL_3PT("C_PW_Kl", tK, tJ); - std::string C_PW_Ksb = LABEL_3PT("C_PW_Ksb", tK, tJ); - std::string C_PW_pilb = LABEL_3PT("C_PW_pilb", tK, tJ); - std::string C_PW_pil = LABEL_3PT("C_PW_pil", tK, tJ); - mesonContraction(application, 3, q_KlCl_q, q_Ks_k, C_PW_Kl, pmom); - mesonContraction(application, 3, q_Kl_0, q_KsCs_mq, C_PW_Ksb, pmom); - mesonContraction(application, 3, q_pil_0, q_pilCl_q, C_PW_pilb, kmom); - mesonContraction(application, 3, q_pilCl_mq, q_pil_p, C_PW_pil, kmom); - // Wall-Wall, to be done. - - /******************************************************************* - * CONTRACTIONS: 4pt contractions, C & W classes. - ******************************************************************/ - std::string CW_Kl = LABEL_4PT("CW_Kl", tK, tJ, tpi); - std::string CW_Ksb = LABEL_4PT("CW_Ksb", tK, tJ, tpi); - std::string CW_pilb = LABEL_4PT("CW_pilb", tK, tJ, tpi); - std::string CW_pil = LABEL_4PT("CW_pil", tK, tJ, tpi); - weakContractionNonEye(application, 4, q_KlCl_q, q_Ks_k, q_pil_p, q_pil_0, CW_Kl); - weakContractionNonEye(application, 4, q_Kl_0, q_KsCs_mq, q_pil_p, q_pil_0, CW_Ksb); - weakContractionNonEye(application, 4, q_Kl_0, q_Ks_k, q_pilCl_q, q_pil_0, CW_pilb); - weakContractionNonEye(application, 4, q_Kl_0, q_Ks_k, q_pil_p, q_pilCl_mq, CW_pil); - - /******************************************************************* - * CONTRACTIONS: 4pt contractions, sd insertions. - ******************************************************************/ - // Note: eventually will use wall sink smeared q_Kl_0/q_KlCl_q instead. - std::string sd_Kl = LABEL_4PT("sd_Kl", tK, tJ, tpi); - std::string sd_Ksb = LABEL_4PT("sd_Ksb", tK, tJ, tpi); - std::string sd_pilb = LABEL_4PT("sd_pilb", tK, tJ, tpi); - gamma3ptContraction(application, 4, q_KlCl_q, q_Ks_k, q_pil_p, sd_Kl); - gamma3ptContraction(application, 4, q_Kl_0, q_KsCs_mq, q_pil_p, sd_Ksb); - gamma3ptContraction(application, 4, q_Kl_0, q_Ks_k, q_pilCl_q, sd_pilb); - - // Sequential sources for each noise propagator. - for (unsigned int nn = 0; nn < n_noise; ++nn) - { - std::string loop_stem = "loop_"; - - // Contraction required for each quark flavour - alternatively - // drop the strange loop if not performing disconnected - // contractions or neglecting H_W operators Q_3 -> Q_10. - for (unsigned int f = 0; f < flavour.size(); ++f) - { - if ((f != strange) || do_disconnected) - { - std::string eta = noiseSrcs[nn]; - std::string loop_q = noiseProps[f][nn]; - std::string loop_qCq = LABEL_3PT(loop_stem + flavour[f], tJ, nn); - std::string loop_qCq_res = loop_qCq + "_res"; - MAKE_SEQUENTIAL_PROP(tJ, noiseRes[f][nn], qmom, - loop_qCq_res, solvers[f], gT); - makeLoop(application, loop_qCq, eta, loop_qCq_res); - - /******************************************************* - * CONTRACTIONS: 4pt contractions, S & E classes. - ******************************************************/ - // Note: eventually will use wall sink smeared q_Kl_0/q_KlCl_q instead. - std::string SE_Kl = LABEL_4PT_NOISE("SE_Kl", tK, tJ, tpi, nn); - std::string SE_Ksb = LABEL_4PT_NOISE("SE_Ksb", tK, tJ, tpi, nn); - std::string SE_pilb = LABEL_4PT_NOISE("SE_pilb", tK, tJ, tpi, nn); - std::string SE_loop = LABEL_4PT_NOISE("SE_loop", tK, tJ, tpi, nn); - weakContractionEye(application, 4, q_KlCl_q, q_Ks_k, q_pil_p, loop_q, SE_Kl); - weakContractionEye(application, 4, q_Kl_0, q_KsCs_mq, q_pil_p, loop_q, SE_Ksb); - weakContractionEye(application, 4, q_Kl_0, q_Ks_k, q_pilCl_q, loop_q, SE_pilb); - weakContractionEye(application, 4, q_Kl_0, q_Ks_k, q_pil_p, loop_qCq, SE_loop); - - /******************************************************* - * CONTRACTIONS: 4pt contractions, pi0 disconnected - * loop. - ******************************************************/ - std::string disc0 = LABEL_4PT_NOISE("disc0", tK, tJ, tpi, nn); - disc0Contraction(application, q_Kl_0, q_Ks_k, q_pilCl_q, loop_q, disc0); - - /******************************************************* - * CONTRACTIONS: Disconnected loop. - ******************************************************/ - std::string discLoop = "disc_" + loop_qCq; - discLoopContraction(application, loop_qCq, discLoop); - } - } - } - } - } - // execution - std::string par_file_name = "rarekaon_000_100_tK0_tpi16_tJ8_noloop_mc0.2.xml"; - application.saveParameterFile(par_file_name); - application.run(); - - // epilogue - LOG(Message) << "Grid is finalizing now" << std::endl; - Grid_finalize(); - - return EXIT_SUCCESS; -} From 18211eb5b13fa2738061b4c48a7518c7368c2645 Mon Sep 17 00:00:00 2001 From: Lanny91 Date: Thu, 22 Jun 2017 16:03:59 +0200 Subject: [PATCH 39/50] Hadrons: Fixed test to use new implementation of meson module. --- tests/hadrons/Test_hadrons.hpp | 20 ++++++----- tests/hadrons/Test_hadrons_meson_3pt.cc | 47 ++++++++++--------------- 2 files changed, 30 insertions(+), 37 deletions(-) diff --git a/tests/hadrons/Test_hadrons.hpp b/tests/hadrons/Test_hadrons.hpp index a554425d..3492816d 100644 --- a/tests/hadrons/Test_hadrons.hpp +++ b/tests/hadrons/Test_hadrons.hpp @@ -71,6 +71,9 @@ using namespace Hadrons; #define NAME_WALL_SOURCE(t) NAME_3MOM_WALL_SOURCE(t, ZERO_MOM) #define NAME_POINT_SOURCE(pos) ("point_" + pos) +// Meson module "gammas" special values +#define ALL_GAMMAS "all" + #define MAKE_3MOM_WALL_PROP(tW, mom, propName, solver)\ {\ std::string srcName = NAME_3MOM_WALL_SOURCE(tW, mom);\ @@ -364,28 +367,27 @@ inline void makeLoop(Application &application, std::string &propName, * Name: mesonContraction * Purpose: Create meson contraction module and add to application module. * Parameters: application - main application that stores modules. - * npt - specify n-point correlator (for labelling). + * modName - unique module name. + * output - name of output files. * q1 - quark propagator 1. * q2 - quark propagator 2. - * label - unique label to construct module name. - * mom - momentum to project (default is zero) + * sink - sink smearing module. * gammas - gamma insertions at source and sink. * Returns: None. ******************************************************************************/ -inline void mesonContraction(Application &application, unsigned int npt, +inline void mesonContraction(Application &application, + std::string &modName, std::string &output, std::string &q1, std::string &q2, - std::string &label, - std::string mom = ZERO_MOM, + std::string &sink, std::string gammas = "") { - std::string modName = std::to_string(npt) + "pt_" + label; if (!(Environment::getInstance().hasModule(modName))) { MContraction::Meson::Par mesPar; - mesPar.output = std::to_string(npt) + "pt/" + label; + mesPar.output = output; mesPar.q1 = q1; mesPar.q2 = q2; - mesPar.mom = mom; + mesPar.sink = sink; mesPar.gammas = gammas; application.createModule(modName, mesPar); } diff --git a/tests/hadrons/Test_hadrons_meson_3pt.cc b/tests/hadrons/Test_hadrons_meson_3pt.cc index 7e487153..1cbb866d 100644 --- a/tests/hadrons/Test_hadrons_meson_3pt.cc +++ b/tests/hadrons/Test_hadrons_meson_3pt.cc @@ -25,7 +25,7 @@ directory. *******************************************************************************/ -#include +#include "Test_hadrons.hpp" using namespace Grid; using namespace Hadrons; @@ -127,43 +127,34 @@ int main(int argc, char *argv[]) } } + // Point sink. + std::string sink = "sink"; + MSink::Point::Par sinkPar; + sinkPar.mom = ZERO_MOM; + application.createModule(sink, sinkPar); + // contractions MContraction::Meson::Par mesPar; for (unsigned int i = 0; i < flavour.size(); ++i) for (unsigned int j = i; j < flavour.size(); ++j) { - mesPar.output = "mesons/Z2_" + flavour[i] + flavour[j]; - mesPar.q1 = qName[i]; - mesPar.q2 = qName[j]; - mesPar.gammas = "all"; - mesPar.mom = "0. 0. 0. 0."; - application.createModule("meson_Z2_" - + std::to_string(t) - + "_" - + flavour[i] - + flavour[j], - mesPar); + std::string modName = "meson_Z2_" + std::to_string(t) + "_" + \ + flavour[i] + flavour[j]; + std::string output = "mesons/Z2_" + flavour[i] + flavour[j]; + mesonContraction(application, modName, output, qName[i], qName[j], + sink, ALL_GAMMAS); } for (unsigned int i = 0; i < flavour.size(); ++i) for (unsigned int j = 0; j < flavour.size(); ++j) for (unsigned int mu = 0; mu < Nd; ++mu) { - MContraction::Meson::Par mesPar; - - mesPar.output = "3pt/Z2_" + flavour[i] + flavour[j] + "_" - + std::to_string(mu); - mesPar.q1 = qName[i]; - mesPar.q2 = seqName[j][mu]; - mesPar.gammas = "all"; - mesPar.mom = "0. 0. 0. 0."; - application.createModule("3pt_Z2_" - + std::to_string(t) - + "_" - + flavour[i] - + flavour[j] - + "_" - + std::to_string(mu), - mesPar); + std::string modName = "3pt_Z2_" + std::to_string(t) + "_" + \ + flavour[i] + flavour[j] + "_" + \ + std::to_string(mu); + std::string output = "3pt/Z2_" + flavour[i] + \ + flavour[j] + "_" + std::to_string(mu); + mesonContraction(application, modName, output, + qName[i], seqName[j][mu], sink, ALL_GAMMAS); } } From 7a3bd5c66c3e1f125801f17fcb0b3b4f1de5d274 Mon Sep 17 00:00:00 2001 From: Lanny91 Date: Thu, 22 Jun 2017 16:06:15 +0200 Subject: [PATCH 40/50] Hadrons: new conserved current contraction test (for regression testing) --- .../Test_hadrons_meson_conserved_3pt.cc | 115 ++++++++++++++++++ 1 file changed, 115 insertions(+) create mode 100644 tests/hadrons/Test_hadrons_meson_conserved_3pt.cc diff --git a/tests/hadrons/Test_hadrons_meson_conserved_3pt.cc b/tests/hadrons/Test_hadrons_meson_conserved_3pt.cc new file mode 100644 index 00000000..c9aeb2cc --- /dev/null +++ b/tests/hadrons/Test_hadrons_meson_conserved_3pt.cc @@ -0,0 +1,115 @@ +/******************************************************************************* + Grid physics library, www.github.com/paboyle/Grid + + Source file: tests/hadrons/Test_hadrons_meson_conserved_3pt.cc + + Copyright (C) 2017 + + Author: Andrew Lawson + + This program is free software; you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation; either version 2 of the License, or + (at your option) any later version. + + This program is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License along + with this program; if not, write to the Free Software Foundation, Inc., + 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. + + See the full license in the file "LICENSE" in the top level distribution + directory. + *******************************************************************************/ + +#include "Test_hadrons.hpp" + +using namespace Grid; +using namespace Hadrons; + +int main(int argc, char *argv[]) +{ + // initialization ////////////////////////////////////////////////////////// + HADRONS_DEFAULT_INIT; + + // run setup /////////////////////////////////////////////////////////////// + Application application; + + // actions parameters + double mass = 0.04; + unsigned int Ls = 16; + double M5 = 1.8; + + // kinematics + unsigned int nt = GridDefaultLatt()[Tp]; + unsigned int tSrc = 0; + unsigned int tJ = nt / 4; + std::string kmom = "0. 0. 0. 0."; + std::string pmom = "1. 0. 0. 0."; + + // Global parameters. + HADRONS_DEFAULT_GLOBALS(application); + + // Unit gauge field. + std::string gaugeField = "Unit gauge"; + application.createModule(gaugeField); + + // DWF action + std::string actionName = "DWF"; + makeDWFAction(application, actionName, gaugeField, mass, M5, Ls); + + // Solver + std::string solver = "CG"; + makeRBPrecCGSolver(application, solver, actionName); + + // main test body ////////////////////////////////////////////////////////// + // Point sink modules. + std::string sink_0 = "sink_0"; + std::string sink_p = "sink_p"; + MSink::Point::Par sinkPar; + sinkPar.mom = kmom; + application.createModule(sink_0, sinkPar); + sinkPar.mom = pmom; + application.createModule(sink_p, sinkPar); + + // 2pt pion contraction, zero momentum. + std::string q_0 = "Q_0"; + MAKE_WALL_PROP(tSrc, q_0, solver); + std::string modName = INIT_INDEX("2pt_pion_WP", tSrc); + std::string output = "2pt/pion_WP_0"; + mesonContraction(application, modName, output, q_0, q_0, sink_0); + + // 2pt pion contraction, with momentum p. + std::string q_p = "Q_p"; + MAKE_3MOM_WALL_PROP(tSrc, pmom, q_p, solver); + modName = INIT_INDEX("2pt_pion_WP_p", tSrc); + output = "2pt/pion_WP_p"; + mesonContraction(application, modName, output, q_0, q_p, sink_p); + + // 3pt pion(0) -> pion(p), with sequentially inserted vector current in + // time direction. + std::string qSeq = q_0 + INIT_INDEX("_seq_Vc3", tJ); + std::string q5d = LABEL_5D(q_0); // Need 5D prop for DWF conserved current. + std::string srcName = qSeq + "_src"; + modName = LABEL_3PT("3pt_pion_Vc3", tSrc, tJ); + output = "3pt/pion_Vc3_p"; + makeConservedSequentialSource(application, srcName, q5d, actionName, + tJ, Current::Vector, Tp, pmom); + makePropagator(application, qSeq, srcName, solver); + mesonContraction(application, modName, output, q_0, qSeq, sink_p); + + std::string par_file_name = "conserved_3pt.xml"; + application.saveParameterFile(par_file_name); + application.run(); + + // epilogue + LOG(Message) << "Grid is finalizing now" << std::endl; + Grid_finalize(); + + return EXIT_SUCCESS; +} + + \ No newline at end of file From dc6b2d30d2dea8ede405aef3f1a753eb7f5127f1 Mon Sep 17 00:00:00 2001 From: Lanny91 Date: Thu, 22 Jun 2017 16:09:45 +0200 Subject: [PATCH 41/50] Documentation fix --- extras/Hadrons/Modules/MContraction/Meson.hpp | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/extras/Hadrons/Modules/MContraction/Meson.hpp b/extras/Hadrons/Modules/MContraction/Meson.hpp index 7810326a..b71f7c08 100644 --- a/extras/Hadrons/Modules/MContraction/Meson.hpp +++ b/extras/Hadrons/Modules/MContraction/Meson.hpp @@ -51,8 +51,7 @@ BEGIN_HADRONS_NAMESPACE in a sequence (e.g. ""). Special values: "all" - perform all possible contractions. - - mom: momentum insertion, space-separated float sequence (e.g ".1 .2 1. 0."), - given as multiples of (2*pi) / L. + - sink: module to compute the sink to use in contraction (string). */ /****************************************************************************** From 08b0e472aa46eb0cef6ad00eaab46cce35357781 Mon Sep 17 00:00:00 2001 From: Lanny91 Date: Thu, 22 Jun 2017 16:34:33 +0200 Subject: [PATCH 42/50] Fixed hadrons tests after merge --- tests/hadrons/Test_hadrons.hpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/tests/hadrons/Test_hadrons.hpp b/tests/hadrons/Test_hadrons.hpp index 3492816d..6ea51d72 100644 --- a/tests/hadrons/Test_hadrons.hpp +++ b/tests/hadrons/Test_hadrons.hpp @@ -329,10 +329,10 @@ inline void makePropagator(Application &application, std::string &propName, // If the propagator already exists, don't make the module again. if (!(Environment::getInstance().hasModule(propName))) { - Quark::Par quarkPar; + MFermion::GaugeProp::Par quarkPar; quarkPar.source = srcName; quarkPar.solver = solver; - application.createModule(propName, quarkPar); + application.createModule(propName, quarkPar); } } From af71c63f4ce48ccbe9bfdaf40d4171913483add7 Mon Sep 17 00:00:00 2001 From: Lanny91 Date: Fri, 23 Jun 2017 11:03:12 +0200 Subject: [PATCH 43/50] AVX2 fix --- lib/simd/Grid_avx.h | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/lib/simd/Grid_avx.h b/lib/simd/Grid_avx.h index 57d9064d..f4634432 100644 --- a/lib/simd/Grid_avx.h +++ b/lib/simd/Grid_avx.h @@ -711,7 +711,7 @@ namespace Optimization { v2 = _mm256_hadd_epi32(v1, v1); u1 = _mm256_castsi256_si128(v2); // upper half u2 = _mm256_extracti128_si256(v2, 1); // lower half - ret = _mm256_add_epi32(u1, u2); + ret = _mm_add_epi32(u1, u2); #else // No AVX horizontal add; extract upper and lower halves of register & use // SSE intrinsics. From 56abbdf4c2fa3848fe9037cf95cf5e4930631d3a Mon Sep 17 00:00:00 2001 From: Lanny91 Date: Fri, 23 Jun 2017 11:09:14 +0200 Subject: [PATCH 44/50] AVX512 integer reduce fix (for non-intel compiler) --- lib/simd/Grid_avx512.h | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/lib/simd/Grid_avx512.h b/lib/simd/Grid_avx512.h index 458a8f7c..85d27421 100644 --- a/lib/simd/Grid_avx512.h +++ b/lib/simd/Grid_avx512.h @@ -558,7 +558,7 @@ namespace Optimization { v2 = _mm256_hadd_epi32(v1, v1); u1 = _mm256_castsi256_si128(v2) // upper half u2 = _mm256_extracti128_si256(v2, 1); // lower half - ret = _mm256_add_epi32(u1, u2); + ret = _mm_add_epi32(u1, u2); return _mm_cvtsi128_si32(ret); } #else From 852ade029a64c6376d391205e607ae655c6d1c80 Mon Sep 17 00:00:00 2001 From: Lanny91 Date: Sun, 16 Jul 2017 13:41:47 +0100 Subject: [PATCH 45/50] Hadrons: Added module to sink a propagator --- extras/Hadrons/Modules.hpp | 1 + extras/Hadrons/Modules/MSink/Smear.hpp | 99 ++++++++++++++++++++++++++ extras/Hadrons/modules.inc | 1 + 3 files changed, 101 insertions(+) create mode 100644 extras/Hadrons/Modules/MSink/Smear.hpp diff --git a/extras/Hadrons/Modules.hpp b/extras/Hadrons/Modules.hpp index 6e1b8823..d0d0d80d 100644 --- a/extras/Hadrons/Modules.hpp +++ b/extras/Hadrons/Modules.hpp @@ -19,6 +19,7 @@ #include #include #include +#include #include #include #include diff --git a/extras/Hadrons/Modules/MSink/Smear.hpp b/extras/Hadrons/Modules/MSink/Smear.hpp new file mode 100644 index 00000000..9327001f --- /dev/null +++ b/extras/Hadrons/Modules/MSink/Smear.hpp @@ -0,0 +1,99 @@ +#ifndef Hadrons_MSink_Smear_hpp_ +#define Hadrons_MSink_Smear_hpp_ + +#include +#include +#include + +BEGIN_HADRONS_NAMESPACE + +/****************************************************************************** + * Smear * + ******************************************************************************/ +BEGIN_MODULE_NAMESPACE(MSink) + +class SmearPar: Serializable +{ +public: + GRID_SERIALIZABLE_CLASS_MEMBERS(SmearPar, + std::string, q, + std::string, sink); +}; + +template +class TSmear: public Module +{ +public: + FERM_TYPE_ALIASES(FImpl,); + SINK_TYPE_ALIASES(); +public: + // constructor + TSmear(const std::string name); + // destructor + virtual ~TSmear(void) = default; + // dependency relation + virtual std::vector getInput(void); + virtual std::vector getOutput(void); + // setup + virtual void setup(void); + // execution + virtual void execute(void); +}; + +MODULE_REGISTER_NS(Smear, TSmear, MSink); + +/****************************************************************************** + * TSmear implementation * + ******************************************************************************/ +// constructor ///////////////////////////////////////////////////////////////// +template +TSmear::TSmear(const std::string name) +: Module(name) +{} + +// dependencies/products /////////////////////////////////////////////////////// +template +std::vector TSmear::getInput(void) +{ + std::vector in = {par().q, par().sink}; + + return in; +} + +template +std::vector TSmear::getOutput(void) +{ + std::vector out = {getName()}; + + return out; +} + +// setup /////////////////////////////////////////////////////////////////////// +template +void TSmear::setup(void) +{ + unsigned int nt = env().getDim(Tp); + unsigned int size = nt * sizeof(SitePropagator); + env().registerObject(getName(), size); +} + +// execution /////////////////////////////////////////////////////////////////// +template +void TSmear::execute(void) +{ + LOG(Message) << "Sink smearing propagator '" << par().q + << "' using sink function '" << par().sink << "'." + << std::endl; + + SinkFn &sink = *env().template getObject(par().sink); + PropagatorField &q = *env().template getObject(par().q); + SlicedPropagator *out = new SlicedPropagator(env().getDim(Tp)); + *out = sink(q); + env().setObject(getName(), out); +} + +END_MODULE_NAMESPACE + +END_HADRONS_NAMESPACE + +#endif // Hadrons_MSink_Smear_hpp_ diff --git a/extras/Hadrons/modules.inc b/extras/Hadrons/modules.inc index 91d0bbe1..fbbb2eb9 100644 --- a/extras/Hadrons/modules.inc +++ b/extras/Hadrons/modules.inc @@ -31,6 +31,7 @@ modules_hpp =\ Modules/MScalar/FreeProp.hpp \ Modules/MScalar/Scalar.hpp \ Modules/MSink/Point.hpp \ + Modules/MSink/Smear.hpp \ Modules/MSolver/RBPrecCG.hpp \ Modules/MSource/Point.hpp \ Modules/MSource/SeqConserved.hpp \ From 6293d438cd6ff2201300298ba29b985962991202 Mon Sep 17 00:00:00 2001 From: Lanny91 Date: Sun, 16 Jul 2017 13:43:25 +0100 Subject: [PATCH 46/50] Hadrons: sink smearing compatibility for 3pt contraction modules. --- .../Hadrons/Modules/MContraction/Gamma3pt.hpp | 22 +++++++++++++--- .../Modules/MContraction/WeakHamiltonian.hpp | 1 + .../MContraction/WeakHamiltonianEye.cc | 25 +++++++++++-------- 3 files changed, 34 insertions(+), 14 deletions(-) diff --git a/extras/Hadrons/Modules/MContraction/Gamma3pt.hpp b/extras/Hadrons/Modules/MContraction/Gamma3pt.hpp index 7f643d49..162ab786 100644 --- a/extras/Hadrons/Modules/MContraction/Gamma3pt.hpp +++ b/extras/Hadrons/Modules/MContraction/Gamma3pt.hpp @@ -51,6 +51,14 @@ BEGIN_HADRONS_NAMESPACE * q1 * * trace(g5*q1*adj(q2)*g5*gamma*q3) + * + * options: + * - q1: sink smeared propagator, source at i + * - q2: propagator, source at i + * - q3: propagator, source at f + * - gamma: gamma matrix to insert + * - tSnk: sink position for propagator q1. + * */ /****************************************************************************** @@ -66,6 +74,7 @@ public: std::string, q2, std::string, q3, Gamma::Algebra, gamma, + unsigned int, tSnk, std::string, output); }; @@ -140,17 +149,22 @@ void TGamma3pt::execute(void) << par().q3 << "', with " << par().gamma << " insertion." << std::endl; + // Initialise variables. q2 and q3 are normal propagators, q1 may be + // sink smeared. CorrWriter writer(par().output); - PropagatorField1 &q1 = *env().template getObject(par().q1); + SlicedPropagator1 &q1 = *env().template getObject(par().q1); PropagatorField2 &q2 = *env().template getObject(par().q2); - PropagatorField3 &q3 = *env().template getObject(par().q3); + PropagatorField3 &q3 = *env().template getObject(par().q3); LatticeComplex c(env().getGrid()); Gamma g5(Gamma::Algebra::Gamma5); Gamma gamma(par().gamma); std::vector buf; Result result; - - c = trace(g5*q1*adj(q2)*(g5*gamma)*q3); + + // Extract relevant timeslice of sinked propagator q1, then contract & + // sum over all spacial positions of gamma insertion. + SitePropagator1 q1Snk = q1[par().tSnk]; + c = trace(g5*q1Snk*adj(q2)*(g5*gamma)*q3); sliceSum(c, buf, Tp); result.gamma = par().gamma; diff --git a/extras/Hadrons/Modules/MContraction/WeakHamiltonian.hpp b/extras/Hadrons/Modules/MContraction/WeakHamiltonian.hpp index 0a3c2e31..302b207e 100644 --- a/extras/Hadrons/Modules/MContraction/WeakHamiltonian.hpp +++ b/extras/Hadrons/Modules/MContraction/WeakHamiltonian.hpp @@ -76,6 +76,7 @@ public: std::string, q2, std::string, q3, std::string, q4, + unsigned int, tSnk, std::string, output); }; diff --git a/extras/Hadrons/Modules/MContraction/WeakHamiltonianEye.cc b/extras/Hadrons/Modules/MContraction/WeakHamiltonianEye.cc index a44c2534..314b080a 100644 --- a/extras/Hadrons/Modules/MContraction/WeakHamiltonianEye.cc +++ b/extras/Hadrons/Modules/MContraction/WeakHamiltonianEye.cc @@ -54,6 +54,8 @@ using namespace MContraction; * * S: trace(q3*g5*q1*adj(q2)*g5*gL[mu][p_1]*q4*gL[mu][p_2]) * E: trace(q3*g5*q1*adj(q2)*g5*gL[mu][p_1])*trace(q4*gL[mu][p_2]) + * + * Note q1 must be sink smeared. */ /****************************************************************************** @@ -94,15 +96,15 @@ void TWeakHamiltonianEye::execute(void) << "'." << std::endl; CorrWriter writer(par().output); - PropagatorField &q1 = *env().template getObject(par().q1); - PropagatorField &q2 = *env().template getObject(par().q2); - PropagatorField &q3 = *env().template getObject(par().q3); - PropagatorField &q4 = *env().template getObject(par().q4); - Gamma g5 = Gamma(Gamma::Algebra::Gamma5); - LatticeComplex expbuf(env().getGrid()); - std::vector corrbuf; - std::vector result(n_eye_diag); - unsigned int ndim = env().getNd(); + SlicedPropagator &q1 = *env().template getObject(par().q1); + PropagatorField &q2 = *env().template getObject(par().q2); + PropagatorField &q3 = *env().template getObject(par().q3); + PropagatorField &q4 = *env().template getObject(par().q4); + Gamma g5 = Gamma(Gamma::Algebra::Gamma5); + LatticeComplex expbuf(env().getGrid()); + std::vector corrbuf; + std::vector result(n_eye_diag); + unsigned int ndim = env().getNd(); PropagatorField tmp1(env().getGrid()); LatticeComplex tmp2(env().getGrid()); @@ -111,10 +113,13 @@ void TWeakHamiltonianEye::execute(void) std::vector E_body(ndim, tmp2); std::vector E_loop(ndim, tmp2); + // Get sink timeslice of q1. + SitePropagator q1Snk = q1[par().tSnk]; + // Setup for S-type contractions. for (int mu = 0; mu < ndim; ++mu) { - S_body[mu] = MAKE_SE_BODY(q1, q2, q3, GammaL(Gamma::gmu[mu])); + S_body[mu] = MAKE_SE_BODY(q1Snk, q2, q3, GammaL(Gamma::gmu[mu])); S_loop[mu] = MAKE_SE_LOOP(q4, GammaL(Gamma::gmu[mu])); } From 0366288b1c8a42ff15eacb0f5e23eee2e89fb50f Mon Sep 17 00:00:00 2001 From: Lanny91 Date: Sun, 16 Jul 2017 13:45:55 +0100 Subject: [PATCH 47/50] Hadrons: added tests for 3pt contractions. --- tests/hadrons/Test_hadrons.hpp | 75 ++++++++--- .../hadrons/Test_hadrons_3pt_contractions.cc | 122 ++++++++++++++++++ 2 files changed, 182 insertions(+), 15 deletions(-) create mode 100644 tests/hadrons/Test_hadrons_3pt_contractions.cc diff --git a/tests/hadrons/Test_hadrons.hpp b/tests/hadrons/Test_hadrons.hpp index 6ea51d72..9bd3ee0a 100644 --- a/tests/hadrons/Test_hadrons.hpp +++ b/tests/hadrons/Test_hadrons.hpp @@ -269,6 +269,26 @@ inline void makeConservedSequentialSource(Application &application, } } +/******************************************************************************* + * Name: makeNoiseSource + * Parameters: application - main application that stores modules. + * srcName - name of source module to create. + * tA - lower source timeslice limit. + * tB - upper source timeslice limit. + * Returns: None. + ******************************************************************************/ +inline void makeNoiseSource(Application &application, std::string &srcName, + unsigned int tA, unsigned int tB) +{ + if (!(Environment::getInstance().hasModule(srcName))) + { + MSource::Z2::Par noisePar; + noisePar.tA = tA; + noisePar.tB = tB; + application.createModule(srcName, noisePar); + } + } + /******************************************************************************* * Name: makeWallSource * Purpose: Construct wall source and add to application module. @@ -292,26 +312,46 @@ inline void makeWallSource(Application &application, std::string &srcName, } /******************************************************************************* - * Name: makeWallSink - * Purpose: Wall sink smearing of a propagator. + * Name: makePointSink + * Purpose: Create function for point sink smearing of a propagator. * Parameters: application - main application that stores modules. * propName - name of input propagator. - * wallName - name of smeared propagator. + * sinkFnct - name of output sink smearing module. * mom - momentum insertion (default is zero). * Returns: None. ******************************************************************************/ -inline void makeWallSink(Application &application, std::string &propName, - std::string &wallName, std::string mom = ZERO_MOM) +inline void makePointSink(Application &application, std::string &sinkFnct, + std::string mom = ZERO_MOM) +{ + // If the sink function already exists, don't make it again. + if (!(Environment::getInstance().hasModule(sinkFnct))) + { + MSink::Point::Par pointPar; + pointPar.mom = mom; + application.createModule(sinkFnct, pointPar); + } +} + +/******************************************************************************* + * Name: sinkSmear + * Purpose: Perform sink smearing of a propagator. + * Parameters: application - main application that stores modules. + * sinkFnct - sink smearing module. + * propName - propagator to smear. + * smearedProp - name of output smeared propagator. + * Returns: None. + ******************************************************************************/ +inline void sinkSmear(Application &application, std::string &sinkFnct, + std::string &propName, std::string &smearedProp) { // If the propagator has already been smeared, don't smear it again. - // Temporarily removed, strategy for sink smearing likely to change. - /*if (!(Environment::getInstance().hasModule(wallName))) + if (!(Environment::getInstance().hasModule(smearedProp))) { - MSink::Wall::Par wallPar; - wallPar.q = propName; - wallPar.mom = mom; - application.createModule(wallName, wallPar); - }*/ + MSink::Smear::Par smearPar; + smearPar.q = propName; + smearPar.sink = sinkFnct; + application.createModule(smearedProp, smearPar); + } } /******************************************************************************* @@ -398,16 +438,18 @@ inline void mesonContraction(Application &application, * Purpose: Create gamma3pt contraction module and add to application module. * Parameters: application - main application that stores modules. * npt - specify n-point correlator (for labelling). - * q1 - quark propagator 1. + * q1 - quark propagator 1, sink smeared. * q2 - quark propagator 2. * q3 - quark propagator 3. * label - unique label to construct module name. + * tSnk - sink position of sink for q1. * gamma - gamma insertions between q2 and q3. * Returns: None. ******************************************************************************/ inline void gamma3ptContraction(Application &application, unsigned int npt, std::string &q1, std::string &q2, - std::string &q3, std::string &label, + std::string &q3, std::string &label, + unsigned int tSnk = 0, Gamma::Algebra gamma = Gamma::Algebra::Identity) { std::string modName = std::to_string(npt) + "pt_" + label; @@ -418,6 +460,7 @@ inline void gamma3ptContraction(Application &application, unsigned int npt, gamma3ptPar.q1 = q1; gamma3ptPar.q2 = q2; gamma3ptPar.q3 = q3; + gamma3ptPar.tSnk = tSnk; gamma3ptPar.gamma = gamma; application.createModule(modName, gamma3ptPar); } @@ -434,13 +477,14 @@ inline void gamma3ptContraction(Application &application, unsigned int npt, * q3 - quark propagator 3. * q4 - quark propagator 4. * label - unique label to construct module name. + * tSnk - time position of sink (for sink smearing). * Returns: None. ******************************************************************************/ #define HW_CONTRACTION(top) \ inline void weakContraction##top(Application &application, unsigned int npt,\ std::string &q1, std::string &q2, \ std::string &q3, std::string &q4, \ - std::string &label)\ + std::string &label, unsigned int tSnk = 0)\ {\ std::string modName = std::to_string(npt) + "pt_" + label;\ if (!(Environment::getInstance().hasModule(modName)))\ @@ -451,6 +495,7 @@ inline void weakContraction##top(Application &application, unsigned int npt,\ weakPar.q2 = q2;\ weakPar.q3 = q3;\ weakPar.q4 = q4;\ + weakPar.tSnk = tSnk;\ application.createModule(modName, weakPar);\ }\ } diff --git a/tests/hadrons/Test_hadrons_3pt_contractions.cc b/tests/hadrons/Test_hadrons_3pt_contractions.cc new file mode 100644 index 00000000..452fc34d --- /dev/null +++ b/tests/hadrons/Test_hadrons_3pt_contractions.cc @@ -0,0 +1,122 @@ +/******************************************************************************* + Grid physics library, www.github.com/paboyle/Grid + + Source file: tests/hadrons/Test_hadrons_3pt_contractions.cc + + Copyright (C) 2017 + + Author: Andrew Lawson + + This program is free software; you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation; either version 2 of the License, or + (at your option) any later version. + + This program is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License along + with this program; if not, write to the Free Software Foundation, Inc., + 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. + + See the full license in the file "LICENSE" in the top level distribution + directory. + *******************************************************************************/ + +#include "Test_hadrons.hpp" + +using namespace Grid; +using namespace Hadrons; + +int main(int argc, char *argv[]) +{ + // initialization ////////////////////////////////////////////////////////// + HADRONS_DEFAULT_INIT; + + // run setup /////////////////////////////////////////////////////////////// + Application application; + double mass = 0.04; + double M5 = 1.8; + unsigned int Ls = 12; + unsigned int nt = GridDefaultLatt()[Tp]; + unsigned int t_i = 0; + unsigned int t_f = nt / 2; + std::string mom = "1. 0. 0. 0."; + + // global parameters + HADRONS_DEFAULT_GLOBALS(application); + + // gauge field + std::string gaugeField = "gauge"; + application.createModule(gaugeField); + + // Action & solver setup. + std::string action = "DWF"; + std::string solver = "CG"; + makeDWFAction(application, action, gaugeField, mass, M5, Ls); + makeRBPrecCGSolver(application, solver, action); + + /*************************************************************************** + * Weak Contraction test: Non-Eye class. + **************************************************************************/ + // Make wall source propagators for each leg of 4-quark vertex. + std::string q_i_0 = "q_i_0"; + std::string q_i_p = "q_i_p"; + std::string q_f_0 = "q_f_0"; + std::string q_f_p = "q_f_p"; + MAKE_WALL_PROP(t_i, q_i_0, solver); + MAKE_WALL_PROP(t_f, q_f_0, solver); + MAKE_3MOM_WALL_PROP(t_i, mom, q_i_p, solver); + MAKE_3MOM_WALL_PROP(t_f, mom, q_f_p, solver); + + // Perform contractions, zero and non-zero momentum. + std::string HW_CW_0 = LABEL_3PT("HW_CW_0", t_i, t_f); + std::string HW_CW_p = LABEL_3PT("HW_CW_p", t_i, t_f); + weakContractionNonEye(application, 3, q_i_0, q_i_0, q_f_0, q_f_0, HW_CW_0); + weakContractionNonEye(application, 3, q_i_0, q_i_p, q_f_p, q_f_0, HW_CW_p); + + /*************************************************************************** + * Weak Contraction test: Eye-class. + **************************************************************************/ + // Create random propagator for loop. + std::string eta = "noise_source"; + makeNoiseSource(application, eta, 0, nt - 1); + std::string loopProp = "loop"; + std::string loopRes = loopProp + "_res"; + makePropagator(application, loopRes, eta, solver); + makeLoop(application, loopProp, eta, loopRes); + + // Wall sink smear the propagator directly connecting the source & sink. + // (i.e. make point sink but smear before the contraction) + std::string wallSink = "wall_sink"; + std::string qWall = "q_wall"; + makePointSink(application, wallSink); + sinkSmear(application, wallSink, q_i_0, qWall); + + // Perform contractions, zero and non-zero momentum. + std::string HW_SE_0 = LABEL_3PT("HW_SE_0", t_i, t_f); + std::string HW_SE_p = LABEL_3PT("HW_SE_p", t_i, t_f); + weakContractionEye(application, 3, qWall, q_i_0, q_f_p, loopProp, HW_SE_0, t_f); + weakContractionEye(application, 3, qWall, q_i_p, q_f_p, loopProp, HW_SE_p, t_f); + + /*************************************************************************** + * Gamma insertion test. + **************************************************************************/ + Gamma::Algebra gamma = Gamma::Algebra::GammaT; + std::string sd_0 = LABEL_3PT("sd_0", t_i, t_f); + std::string sd_p = LABEL_3PT("sd_p", t_i, t_f); + gamma3ptContraction(application, 3, qWall, q_i_0, q_f_0, sd_0, t_f, gamma); + gamma3ptContraction(application, 3, qWall, q_i_p, q_f_p, sd_p, t_f, gamma); + + // execution + application.saveParameterFile("ContractionTest3pt.xml"); + application.run(); + + // epilogue + LOG(Message) << "Grid is finalizing now" << std::endl; + Grid_finalize(); + + return EXIT_SUCCESS; +} \ No newline at end of file From 875e1a841f24166084cc26e16aea363c1200070c Mon Sep 17 00:00:00 2001 From: Lanny91 Date: Sun, 16 Jul 2017 13:47:00 +0100 Subject: [PATCH 48/50] Hadrons: updated Quark -> MFermion/GaugeProp module name in test. --- tests/hadrons/Test_hadrons_quark.cc | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tests/hadrons/Test_hadrons_quark.cc b/tests/hadrons/Test_hadrons_quark.cc index 5b9d0ce1..eac065e9 100644 --- a/tests/hadrons/Test_hadrons_quark.cc +++ b/tests/hadrons/Test_hadrons_quark.cc @@ -26,7 +26,7 @@ *******************************************************************************/ #include "Test_hadrons.hpp" -#include +#include using namespace Grid; using namespace QCD; From 67b34e5789aec1b39d34c2cdbedb156ff9509e11 Mon Sep 17 00:00:00 2001 From: Lanny91 Date: Mon, 31 Jul 2017 11:35:01 +0100 Subject: [PATCH 49/50] Modified conserved current 5th dimension loop for compatibility with 5D vectorisation. --- lib/qcd/action/fermion/WilsonFermion5D.cc | 61 +++++++++++++++++------ 1 file changed, 46 insertions(+), 15 deletions(-) diff --git a/lib/qcd/action/fermion/WilsonFermion5D.cc b/lib/qcd/action/fermion/WilsonFermion5D.cc index 5daed3de..5ddfde9a 100644 --- a/lib/qcd/action/fermion/WilsonFermion5D.cc +++ b/lib/qcd/action/fermion/WilsonFermion5D.cc @@ -12,6 +12,7 @@ Author: Peter Boyle Author: Peter Boyle Author: paboyle Author: Guido Cossu +Author: Andrew Lawson This program is free software; you can redistribute it and/or modify it under the terms of the GNU General Public License as published by @@ -676,6 +677,21 @@ void WilsonFermion5D::MomentumSpacePropagatorHw(FermionField &out,const Fe * to make a conserved current sink or inserting the conserved current * sequentially. ******************************************************************************/ + +// Helper macro to reverse Simd vector. Fixme: slow, generic implementation. +#define REVERSE_LS(qSite, qSiteRev, Nsimd) \ +{ \ + std::vector qSiteVec(Nsimd); \ + extract(qSite, qSiteVec); \ + for (int i = 0; i < Nsimd / 2; ++i) \ + { \ + typename SitePropagator::scalar_object tmp = qSiteVec[i]; \ + qSiteVec[i] = qSiteVec[Nsimd - i - 1]; \ + qSiteVec[Nsimd - i - 1] = tmp; \ + } \ + merge(qSiteRev, qSiteVec); \ +} + template void WilsonFermion5D::ContractConservedCurrent(PropagatorField &q_in_1, PropagatorField &q_in_2, @@ -687,6 +703,7 @@ void WilsonFermion5D::ContractConservedCurrent(PropagatorField &q_in_1, conformable(q_in_1._grid, q_in_2._grid); conformable(_FourDimGrid, q_out._grid); PropagatorField tmp1(FermionGrid()), tmp2(FermionGrid()); + unsigned int LLs = q_in_1._grid->_rdimensions[0]; q_out = zero; // Forward, need q1(x + mu, s), q2(x, Ls - 1 - s). Backward, need q1(x, s), @@ -695,18 +712,33 @@ void WilsonFermion5D::ContractConservedCurrent(PropagatorField &q_in_1, tmp2 = Cshift(q_in_2, mu + 1, 1); parallel_for (unsigned int sU = 0; sU < Umu._grid->oSites(); ++sU) { - unsigned int sF1 = sU * Ls; - unsigned int sF2 = (sU + 1) * Ls - 1; - for (int s = 0; s < Ls; ++s) + unsigned int sF1 = sU * LLs; + unsigned int sF2 = (sU + 1) * LLs - 1; + + for (unsigned int s = 0; s < LLs; ++s) { - bool axial_sign = ((curr_type == Current::Axial) && (s < (Ls / 2))) ? \ - true : false; + bool axial_sign = ((curr_type == Current::Axial) && \ + (s < (LLs / 2))); + SitePropagator qSite2, qmuSite2; + + // If vectorised in 5th dimension, reverse q2 vector to match up + // sites correctly. + if (Impl::LsVectorised) + { + REVERSE_LS(q_in_2._odata[sF2], qSite2, Ls / LLs); + REVERSE_LS(tmp2._odata[sF2], qmuSite2, Ls / LLs); + } + else + { + qSite2 = q_in_2._odata[sF2]; + qmuSite2 = tmp2._odata[sF2]; + } Kernels::ContractConservedCurrentSiteFwd(tmp1._odata[sF1], - q_in_2._odata[sF2], + qSite2, q_out._odata[sU], Umu, sU, mu, axial_sign); Kernels::ContractConservedCurrentSiteBwd(q_in_1._odata[sF1], - tmp2._odata[sF2], + qmuSite2, q_out._odata[sU], Umu, sU, mu, axial_sign); sF1++; @@ -732,6 +764,7 @@ void WilsonFermion5D::SeqConservedCurrent(PropagatorField &q_in, tmp(FermionGrid()); Complex i(0.0, 1.0); int tshift = (mu == Tp) ? 1 : 0; + unsigned int LLs = q_in._grid->_rdimensions[0]; // Momentum projection. ph = zero; @@ -764,11 +797,10 @@ void WilsonFermion5D::SeqConservedCurrent(PropagatorField &q_in, if (timeSlices > 0) { - unsigned int sF = sU * Ls; - for (unsigned int s = 0; s < Ls; ++s) + unsigned int sF = sU * LLs; + for (unsigned int s = 0; s < LLs; ++s) { - bool axial_sign = ((curr_type == Current::Axial) && (s < (Ls / 2))) ? \ - true : false; + bool axial_sign = ((curr_type == Current::Axial) && (s < (LLs / 2))); Kernels::SeqConservedCurrentSiteFwd(tmpFwd._odata[sF], q_out._odata[sF], Umu, sU, mu, t_mask, axial_sign); @@ -783,11 +815,10 @@ void WilsonFermion5D::SeqConservedCurrent(PropagatorField &q_in, if (timeSlices > 0) { - unsigned int sF = sU * Ls; - for (unsigned int s = 0; s < Ls; ++s) + unsigned int sF = sU * LLs; + for (unsigned int s = 0; s < LLs; ++s) { - bool axial_sign = ((curr_type == Current::Axial) && (s < (Ls / 2))) ? \ - true : false; + bool axial_sign = ((curr_type == Current::Axial) && (s < (LLs / 2))); Kernels::SeqConservedCurrentSiteBwd(tmpBwd._odata[sF], q_out._odata[sF], Umu, sU, mu, t_mask, axial_sign); From 323e9c439ab0889d69c60b4736e1bd07d7724c06 Mon Sep 17 00:00:00 2001 From: Lanny91 Date: Mon, 31 Jul 2017 12:26:34 +0100 Subject: [PATCH 50/50] Hadrons: Legal banner fixes --- extras/Hadrons/Modules.hpp | 30 ++++++++++++++++++ extras/Hadrons/Modules/MFermion/GaugeProp.hpp | 31 +++++++++++++++++++ extras/Hadrons/Modules/MSink/Point.hpp | 28 +++++++++++++++++ extras/Hadrons/Modules/MSink/Smear.hpp | 28 +++++++++++++++++ .../Modules/MUtilities/TestSeqGamma.hpp | 28 +++++++++++++++++ 5 files changed, 145 insertions(+) diff --git a/extras/Hadrons/Modules.hpp b/extras/Hadrons/Modules.hpp index d0d0d80d..e1f06f32 100644 --- a/extras/Hadrons/Modules.hpp +++ b/extras/Hadrons/Modules.hpp @@ -1,3 +1,33 @@ +/************************************************************************************* + +Grid physics library, www.github.com/paboyle/Grid + +Source file: extras/Hadrons/Modules.hpp + +Copyright (C) 2015 +Copyright (C) 2016 +Copyright (C) 2017 + +Author: Antonin Portelli + +This program is free software; you can redistribute it and/or modify +it under the terms of the GNU General Public License as published by +the Free Software Foundation; either version 2 of the License, or +(at your option) any later version. + +This program is distributed in the hope that it will be useful, +but WITHOUT ANY WARRANTY; without even the implied warranty of +MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +GNU General Public License for more details. + +You should have received a copy of the GNU General Public License along +with this program; if not, write to the Free Software Foundation, Inc., +51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. + +See the full license in the file "LICENSE" in the top level distribution directory +*************************************************************************************/ +/* END LEGAL */ + #include #include #include diff --git a/extras/Hadrons/Modules/MFermion/GaugeProp.hpp b/extras/Hadrons/Modules/MFermion/GaugeProp.hpp index 050f8381..8add9a00 100644 --- a/extras/Hadrons/Modules/MFermion/GaugeProp.hpp +++ b/extras/Hadrons/Modules/MFermion/GaugeProp.hpp @@ -1,3 +1,34 @@ +/************************************************************************************* + +Grid physics library, www.github.com/paboyle/Grid + +Source file: extras/Hadrons/Modules/MFermion/GaugeProp.hpp + +Copyright (C) 2015 +Copyright (C) 2016 +Copyright (C) 2017 + +Author: Antonin Portelli + Andrew Lawson + +This program is free software; you can redistribute it and/or modify +it under the terms of the GNU General Public License as published by +the Free Software Foundation; either version 2 of the License, or +(at your option) any later version. + +This program is distributed in the hope that it will be useful, +but WITHOUT ANY WARRANTY; without even the implied warranty of +MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +GNU General Public License for more details. + +You should have received a copy of the GNU General Public License along +with this program; if not, write to the Free Software Foundation, Inc., +51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. + +See the full license in the file "LICENSE" in the top level distribution directory +*************************************************************************************/ +/* END LEGAL */ + #ifndef Hadrons_MFermion_GaugeProp_hpp_ #define Hadrons_MFermion_GaugeProp_hpp_ diff --git a/extras/Hadrons/Modules/MSink/Point.hpp b/extras/Hadrons/Modules/MSink/Point.hpp index 7b3aa9de..0761c4c4 100644 --- a/extras/Hadrons/Modules/MSink/Point.hpp +++ b/extras/Hadrons/Modules/MSink/Point.hpp @@ -1,3 +1,31 @@ +/************************************************************************************* + +Grid physics library, www.github.com/paboyle/Grid + +Source file: extras/Hadrons/Modules/MSink/Point.hpp + +Copyright (C) 2017 + +Author: Antonin Portelli + +This program is free software; you can redistribute it and/or modify +it under the terms of the GNU General Public License as published by +the Free Software Foundation; either version 2 of the License, or +(at your option) any later version. + +This program is distributed in the hope that it will be useful, +but WITHOUT ANY WARRANTY; without even the implied warranty of +MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +GNU General Public License for more details. + +You should have received a copy of the GNU General Public License along +with this program; if not, write to the Free Software Foundation, Inc., +51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. + +See the full license in the file "LICENSE" in the top level distribution directory +*************************************************************************************/ +/* END LEGAL */ + #ifndef Hadrons_MSink_Point_hpp_ #define Hadrons_MSink_Point_hpp_ diff --git a/extras/Hadrons/Modules/MSink/Smear.hpp b/extras/Hadrons/Modules/MSink/Smear.hpp index 9327001f..c3973d2b 100644 --- a/extras/Hadrons/Modules/MSink/Smear.hpp +++ b/extras/Hadrons/Modules/MSink/Smear.hpp @@ -1,3 +1,31 @@ +/************************************************************************************* + +Grid physics library, www.github.com/paboyle/Grid + +Source file: extras/Hadrons/Modules/MSink/Smear.hpp + +Copyright (C) 2017 + +Author: Andrew Lawson + +This program is free software; you can redistribute it and/or modify +it under the terms of the GNU General Public License as published by +the Free Software Foundation; either version 2 of the License, or +(at your option) any later version. + +This program is distributed in the hope that it will be useful, +but WITHOUT ANY WARRANTY; without even the implied warranty of +MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +GNU General Public License for more details. + +You should have received a copy of the GNU General Public License along +with this program; if not, write to the Free Software Foundation, Inc., +51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. + +See the full license in the file "LICENSE" in the top level distribution directory +*************************************************************************************/ +/* END LEGAL */ + #ifndef Hadrons_MSink_Smear_hpp_ #define Hadrons_MSink_Smear_hpp_ diff --git a/extras/Hadrons/Modules/MUtilities/TestSeqGamma.hpp b/extras/Hadrons/Modules/MUtilities/TestSeqGamma.hpp index 3dbd7d63..2799e5d0 100644 --- a/extras/Hadrons/Modules/MUtilities/TestSeqGamma.hpp +++ b/extras/Hadrons/Modules/MUtilities/TestSeqGamma.hpp @@ -1,3 +1,31 @@ +/************************************************************************************* + +Grid physics library, www.github.com/paboyle/Grid + +Source file: extras/Hadrons/Modules/MUtilities/TestSeqGamma.hpp + +Copyright (C) 2017 + +Author: Andrew Lawson + +This program is free software; you can redistribute it and/or modify +it under the terms of the GNU General Public License as published by +the Free Software Foundation; either version 2 of the License, or +(at your option) any later version. + +This program is distributed in the hope that it will be useful, +but WITHOUT ANY WARRANTY; without even the implied warranty of +MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +GNU General Public License for more details. + +You should have received a copy of the GNU General Public License along +with this program; if not, write to the Free Software Foundation, Inc., +51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. + +See the full license in the file "LICENSE" in the top level distribution directory +*************************************************************************************/ +/* END LEGAL */ + #ifndef Hadrons_MUtilities_TestSeqGamma_hpp_ #define Hadrons_MUtilities_TestSeqGamma_hpp_