portelli/Grid
1
0
Fork 0
mirror of https://github.com/paboyle/Grid.git synced 2025-04-09 21:50:45 +01:00
Code Releases Activity

Merge pull request #121 from Lanny91/feature/hadrons

Browse Source 
Feature/hadrons
This commit is contained in:
Peter Boyle 2017-08-24 12:59:08 +01:00 committed by GitHub
commit 2bcb704af2
40 changed files with 3075 additions and 54 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

11
extras/Hadrons/Global.hpp
View File

@ -62,12 +62,11 @@ BEGIN_HADRONS_NAMESPACE
// type aliases // type aliases
#define FERM_TYPE_ALIASES(FImpl, suffix)\ #define FERM_TYPE_ALIASES(FImpl, suffix)\
typedef FermionOperator<FImpl> FMat##suffix; \ typedef FermionOperator<FImpl> FMat##suffix; \
typedef typename FImpl::FermionField FermionField##suffix; \ typedef typename FImpl::FermionField FermionField##suffix; \
typedef typename FImpl::PropagatorField PropagatorField##suffix; \ typedef typename FImpl::PropagatorField PropagatorField##suffix; \
typedef typename FImpl::SitePropagator SitePropagator##suffix; \ typedef typename FImpl::SitePropagator::scalar_object SitePropagator##suffix; \
typedef std::vector<typename FImpl::SitePropagator::scalar_object> \ typedef std::vector<SitePropagator##suffix> SlicedPropagator##suffix;
SlicedPropagator##suffix;
#define GAUGE_TYPE_ALIASES(FImpl, suffix)\ #define GAUGE_TYPE_ALIASES(FImpl, suffix)\
typedef typename FImpl::DoubledGaugeField DoubledGaugeField##suffix; typedef typename FImpl::DoubledGaugeField DoubledGaugeField##suffix;

35
extras/Hadrons/Modules.hpp
View File

@ -1,9 +1,40 @@
/*************************************************************************************
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 <antonin.portelli@me.com>
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 <Grid/Hadrons/Modules/MAction/DWF.hpp> #include <Grid/Hadrons/Modules/MAction/DWF.hpp>
#include <Grid/Hadrons/Modules/MAction/Wilson.hpp> #include <Grid/Hadrons/Modules/MAction/Wilson.hpp>
#include <Grid/Hadrons/Modules/MContraction/Baryon.hpp> #include <Grid/Hadrons/Modules/MContraction/Baryon.hpp>
#include <Grid/Hadrons/Modules/MContraction/DiscLoop.hpp> #include <Grid/Hadrons/Modules/MContraction/DiscLoop.hpp>
#include <Grid/Hadrons/Modules/MContraction/Gamma3pt.hpp> #include <Grid/Hadrons/Modules/MContraction/Gamma3pt.hpp>
#include <Grid/Hadrons/Modules/MContraction/Meson.hpp> #include <Grid/Hadrons/Modules/MContraction/Meson.hpp>
#include <Grid/Hadrons/Modules/MContraction/WardIdentity.hpp>
#include <Grid/Hadrons/Modules/MContraction/WeakHamiltonian.hpp> #include <Grid/Hadrons/Modules/MContraction/WeakHamiltonian.hpp>
#include <Grid/Hadrons/Modules/MContraction/WeakHamiltonianEye.hpp> #include <Grid/Hadrons/Modules/MContraction/WeakHamiltonianEye.hpp>
#include <Grid/Hadrons/Modules/MContraction/WeakHamiltonianNonEye.hpp> #include <Grid/Hadrons/Modules/MContraction/WeakHamiltonianNonEye.hpp>
@ -18,8 +49,12 @@
#include <Grid/Hadrons/Modules/MScalar/FreeProp.hpp> #include <Grid/Hadrons/Modules/MScalar/FreeProp.hpp>
#include <Grid/Hadrons/Modules/MScalar/Scalar.hpp> #include <Grid/Hadrons/Modules/MScalar/Scalar.hpp>
#include <Grid/Hadrons/Modules/MSink/Point.hpp> #include <Grid/Hadrons/Modules/MSink/Point.hpp>
#include <Grid/Hadrons/Modules/MSink/Smear.hpp>
#include <Grid/Hadrons/Modules/MSolver/RBPrecCG.hpp> #include <Grid/Hadrons/Modules/MSolver/RBPrecCG.hpp>
#include <Grid/Hadrons/Modules/MSource/Point.hpp> #include <Grid/Hadrons/Modules/MSource/Point.hpp>
#include <Grid/Hadrons/Modules/MSource/SeqConserved.hpp>
#include <Grid/Hadrons/Modules/MSource/SeqGamma.hpp> #include <Grid/Hadrons/Modules/MSource/SeqGamma.hpp>
#include <Grid/Hadrons/Modules/MSource/Wall.hpp> #include <Grid/Hadrons/Modules/MSource/Wall.hpp>
#include <Grid/Hadrons/Modules/MSource/Z2.hpp> #include <Grid/Hadrons/Modules/MSource/Z2.hpp>
#include <Grid/Hadrons/Modules/MUtilities/TestSeqConserved.hpp>
#include <Grid/Hadrons/Modules/MUtilities/TestSeqGamma.hpp>

22
extras/Hadrons/Modules/MContraction/Gamma3pt.hpp
View File

@ -51,6 +51,14 @@ BEGIN_HADRONS_NAMESPACE
* q1 * q1
* *
* trace(g5*q1*adj(q2)*g5*gamma*q3) * 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, q2,
std::string, q3, std::string, q3,
Gamma::Algebra, gamma, Gamma::Algebra, gamma,
unsigned int, tSnk,
std::string, output); std::string, output);
}; };
@ -140,17 +149,22 @@ void TGamma3pt<FImpl1, FImpl2, FImpl3>::execute(void)
<< par().q3 << "', with " << par().gamma << " insertion." << par().q3 << "', with " << par().gamma << " insertion."
<< std::endl; << std::endl;
// Initialise variables. q2 and q3 are normal propagators, q1 may be
// sink smeared.
CorrWriter writer(par().output); CorrWriter writer(par().output);
PropagatorField1 &q1 = *env().template getObject<PropagatorField1>(par().q1); SlicedPropagator1 &q1 = *env().template getObject<SlicedPropagator1>(par().q1);
PropagatorField2 &q2 = *env().template getObject<PropagatorField2>(par().q2); PropagatorField2 &q2 = *env().template getObject<PropagatorField2>(par().q2);
PropagatorField3 &q3 = *env().template getObject<PropagatorField3>(par().q3); PropagatorField3 &q3 = *env().template getObject<PropagatorField2>(par().q3);
LatticeComplex c(env().getGrid()); LatticeComplex c(env().getGrid());
Gamma g5(Gamma::Algebra::Gamma5); Gamma g5(Gamma::Algebra::Gamma5);
Gamma gamma(par().gamma); Gamma gamma(par().gamma);
std::vector<TComplex> buf; std::vector<TComplex> buf;
Result result; 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); sliceSum(c, buf, Tp);
result.gamma = par().gamma; result.gamma = par().gamma;

3
extras/Hadrons/Modules/MContraction/Meson.hpp
View File

@ -51,8 +51,7 @@ BEGIN_HADRONS_NAMESPACE
in a sequence (e.g. "<Gamma5 Gamma5><Gamma5 GammaT>"). in a sequence (e.g. "<Gamma5 Gamma5><Gamma5 GammaT>").
Special values: "all" - perform all possible contractions. Special values: "all" - perform all possible contractions.
- mom: momentum insertion, space-separated float sequence (e.g ".1 .2 1. 0."), - sink: module to compute the sink to use in contraction (string).
given as multiples of (2*pi) / L.
*/ */
/****************************************************************************** /******************************************************************************

211
extras/Hadrons/Modules/MContraction/WardIdentity.hpp Normal file
View File

@ -0,0 +1,211 @@
/*************************************************************************************
Grid physics library, www.github.com/paboyle/Grid
Source file: extras/Hadrons/Modules/MContraction/WardIdentity.hpp
Copyright (C) 2017
Author: Andrew Lawson <andrew.lawson1991@gmail.com>
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_MContraction_WardIdentity_hpp_
#define Hadrons_MContraction_WardIdentity_hpp_
#include <Grid/Hadrons/Global.hpp>
#include <Grid/Hadrons/Module.hpp>
#include <Grid/Hadrons/ModuleFactory.hpp>
BEGIN_HADRONS_NAMESPACE
/*
Ward Identity contractions
-----------------------------
* options:
- 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.
*/
/******************************************************************************
* WardIdentity *
******************************************************************************/
BEGIN_MODULE_NAMESPACE(MContraction)
class WardIdentityPar: Serializable
{
public:
GRID_SERIALIZABLE_CLASS_MEMBERS(WardIdentityPar,
std::string, q,
std::string, action,
double, mass,
bool, test_axial);
};
template <typename FImpl>
class TWardIdentity: public Module<WardIdentityPar>
{
public:
FERM_TYPE_ALIASES(FImpl,);
public:
// constructor
TWardIdentity(const std::string name);
// destructor
virtual ~TWardIdentity(void) = default;
// dependency relation
virtual std::vector<std::string> getInput(void);
virtual std::vector<std::string> getOutput(void);
// setup
virtual void setup(void);
// execution
virtual void execute(void);
private:
unsigned int Ls_;
};
MODULE_REGISTER_NS(WardIdentity, TWardIdentity<FIMPL>, MContraction);
/******************************************************************************
* TWardIdentity implementation *
******************************************************************************/
// constructor /////////////////////////////////////////////////////////////////
template <typename FImpl>
TWardIdentity<FImpl>::TWardIdentity(const std::string name)
: Module<WardIdentityPar>(name)
{}
// dependencies/products ///////////////////////////////////////////////////////
template <typename FImpl>
std::vector<std::string> TWardIdentity<FImpl>::getInput(void)
{
std::vector<std::string> in = {par().q, par().action};
return in;
}
template <typename FImpl>
std::vector<std::string> TWardIdentity<FImpl>::getOutput(void)
{
std::vector<std::string> out = {getName()};
return out;
}
// setup ///////////////////////////////////////////////////////////////////////
template <typename FImpl>
void TWardIdentity<FImpl>::setup(void)
{
Ls_ = env().getObjectLs(par().q);
if (Ls_ != env().getObjectLs(par().action))
{
HADRON_ERROR("Ls mismatch between quark action and propagator");
}
}
// execution ///////////////////////////////////////////////////////////////////
template <typename FImpl>
void TWardIdentity<FImpl>::execute(void)
{
LOG(Message) << "Performing Ward Identity checks for quark '" << par().q
<< "'." << std::endl;
PropagatorField tmp(env().getGrid()), vector_WI(env().getGrid());
PropagatorField &q = *env().template getObject<PropagatorField>(par().q);
FMat &act = *(env().template getObject<FMat>(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)
{
act.ContractConservedCurrent(q, q, tmp, Current::Vector, mu);
tmp -= Cshift(tmp, mu, -1);
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;
if (par().test_axial)
{
PropagatorField psi(env().getGrid());
LatticeComplex PP(env().getGrid()), axial_defect(env().getGrid()),
PJ5q(env().getGrid());
std::vector<TComplex> axial_buf;
// Compute <P|D_mu A_mu>, D is backwards derivative.
axial_defect = zero;
for (unsigned int mu = 0; mu < Nd; ++mu)
{
act.ContractConservedCurrent(q, q, tmp, Current::Axial, mu);
tmp -= Cshift(tmp, mu, -1);
axial_defect += trace(g5*tmp);
}
// Get <P|J5q> for 5D (zero for 4D) and <P|P>.
PJ5q = zero;
if (Ls_ > 1)
{
// <P|P>
ExtractSlice(tmp, q, 0, 0);
psi = 0.5 * (tmp - g5*tmp);
ExtractSlice(tmp, q, Ls_ - 1, 0);
psi += 0.5 * (tmp + g5*tmp);
PP = trace(adj(psi)*psi);
// <P|5Jq>
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 identity <P|D_mu A_mu> = 2m<P|P> + 2<P|J5q>
LOG(Message) << "|D_mu A_mu|^2 = " << norm2(axial_defect) << std::endl;
LOG(Message) << "|PP|^2 = " << norm2(PP) << std::endl;
LOG(Message) << "|PJ5q|^2 = " << norm2(PJ5q) << std::endl;
LOG(Message) << "Axial Ward Identity defect Delta_mu A_mu = "
<< norm2(axial_defect) << std::endl;
// 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;
}
}
}
END_MODULE_NAMESPACE
END_HADRONS_NAMESPACE
#endif // Hadrons_WardIdentity_hpp_

1
extras/Hadrons/Modules/MContraction/WeakHamiltonian.hpp
View File

@ -76,6 +76,7 @@ public:
std::string, q2, std::string, q2,
std::string, q3, std::string, q3,
std::string, q4, std::string, q4,
unsigned int, tSnk,
std::string, output); std::string, output);
}; };

25
extras/Hadrons/Modules/MContraction/WeakHamiltonianEye.cc
View File

@ -54,6 +54,8 @@ using namespace MContraction;
* *
* S: trace(q3*g5*q1*adj(q2)*g5*gL[mu][p_1]*q4*gL[mu][p_2]) * 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]) * 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; << "'." << std::endl;
CorrWriter writer(par().output); CorrWriter writer(par().output);
PropagatorField &q1 = *env().template getObject<PropagatorField>(par().q1); SlicedPropagator &q1 = *env().template getObject<SlicedPropagator>(par().q1);
PropagatorField &q2 = *env().template getObject<PropagatorField>(par().q2); PropagatorField &q2 = *env().template getObject<PropagatorField>(par().q2);
PropagatorField &q3 = *env().template getObject<PropagatorField>(par().q3); PropagatorField &q3 = *env().template getObject<PropagatorField>(par().q3);
PropagatorField &q4 = *env().template getObject<PropagatorField>(par().q4); PropagatorField &q4 = *env().template getObject<PropagatorField>(par().q4);
Gamma g5 = Gamma(Gamma::Algebra::Gamma5); Gamma g5 = Gamma(Gamma::Algebra::Gamma5);
LatticeComplex expbuf(env().getGrid()); LatticeComplex expbuf(env().getGrid());
std::vector<TComplex> corrbuf; std::vector<TComplex> corrbuf;
std::vector<Result> result(n_eye_diag); std::vector<Result> result(n_eye_diag);
unsigned int ndim = env().getNd(); unsigned int ndim = env().getNd();
PropagatorField tmp1(env().getGrid()); PropagatorField tmp1(env().getGrid());
LatticeComplex tmp2(env().getGrid()); LatticeComplex tmp2(env().getGrid());
@ -111,10 +113,13 @@ void TWeakHamiltonianEye::execute(void)
std::vector<LatticeComplex> E_body(ndim, tmp2); std::vector<LatticeComplex> E_body(ndim, tmp2);
std::vector<LatticeComplex> E_loop(ndim, tmp2); std::vector<LatticeComplex> E_loop(ndim, tmp2);
// Get sink timeslice of q1.
SitePropagator q1Snk = q1[par().tSnk];
// Setup for S-type contractions. // Setup for S-type contractions.
for (int mu = 0; mu < ndim; ++mu) 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])); S_loop[mu] = MAKE_SE_LOOP(q4, GammaL(Gamma::gmu[mu]));
} }

65
extras/Hadrons/Modules/MFermion/GaugeProp.hpp
View File

@ -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 <antonin.portelli@me.com>
Andrew Lawson <andrew.lawson1991@gmail.com>
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_ #ifndef Hadrons_MFermion_GaugeProp_hpp_
#define Hadrons_MFermion_GaugeProp_hpp_ #define Hadrons_MFermion_GaugeProp_hpp_
@ -7,6 +38,27 @@
BEGIN_HADRONS_NAMESPACE BEGIN_HADRONS_NAMESPACE
/******************************************************************************
* 5D -> 4D and 4D -> 5D conversions. *
******************************************************************************/
template<class vobj> // Note that 5D object is modified.
inline void make_4D(Lattice<vobj> &in_5d, Lattice<vobj> &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<class vobj>
inline void make_5D(Lattice<vobj> &in_4d, Lattice<vobj> &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);
}
/****************************************************************************** /******************************************************************************
* GaugeProp * * GaugeProp *
******************************************************************************/ ******************************************************************************/
@ -116,12 +168,8 @@ void TGaugeProp<FImpl>::execute(void)
} }
else else
{ {
source = zero;
PropToFerm(tmp, fullSrc, s, c); PropToFerm(tmp, fullSrc, s, c);
InsertSlice(tmp, source, 0, 0); make_5D(tmp, source, Ls_);
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);
} }
} }
// source conversion for 5D sources // source conversion for 5D sources
@ -143,11 +191,8 @@ void TGaugeProp<FImpl>::execute(void)
if (Ls_ > 1) if (Ls_ > 1)
{ {
PropagatorField &p4d = PropagatorField &p4d =
*env().template getObject<PropagatorField>(getName()); *env().template getObject<PropagatorField>(getName());
make_4D(sol, tmp, Ls_);
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);
FermToProp(p4d, tmp, s, c); FermToProp(p4d, tmp, s, c);
} }
} }

28
extras/Hadrons/Modules/MSink/Point.hpp
View File

@ -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 <antonin.portelli@me.com>
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_ #ifndef Hadrons_MSink_Point_hpp_
#define Hadrons_MSink_Point_hpp_ #define Hadrons_MSink_Point_hpp_

127
extras/Hadrons/Modules/MSink/Smear.hpp Normal file
View File

@ -0,0 +1,127 @@
/*************************************************************************************
Grid physics library, www.github.com/paboyle/Grid
Source file: extras/Hadrons/Modules/MSink/Smear.hpp
Copyright (C) 2017
Author: Andrew Lawson <andrew.lawson1991@gmail.com>
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_
#include <Grid/Hadrons/Global.hpp>
#include <Grid/Hadrons/Module.hpp>
#include <Grid/Hadrons/ModuleFactory.hpp>
BEGIN_HADRONS_NAMESPACE
/******************************************************************************
* Smear *
******************************************************************************/
BEGIN_MODULE_NAMESPACE(MSink)
class SmearPar: Serializable
{
public:
GRID_SERIALIZABLE_CLASS_MEMBERS(SmearPar,
std::string, q,
std::string, sink);
};
template <typename FImpl>
class TSmear: public Module<SmearPar>
{
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<std::string> getInput(void);
virtual std::vector<std::string> getOutput(void);
// setup
virtual void setup(void);
// execution
virtual void execute(void);
};
MODULE_REGISTER_NS(Smear, TSmear<FIMPL>, MSink);
/******************************************************************************
* TSmear implementation *
******************************************************************************/
// constructor /////////////////////////////////////////////////////////////////
template <typename FImpl>
TSmear<FImpl>::TSmear(const std::string name)
: Module<SmearPar>(name)
{}
// dependencies/products ///////////////////////////////////////////////////////
template <typename FImpl>
std::vector<std::string> TSmear<FImpl>::getInput(void)
{
std::vector<std::string> in = {par().q, par().sink};
return in;
}
template <typename FImpl>
std::vector<std::string> TSmear<FImpl>::getOutput(void)
{
std::vector<std::string> out = {getName()};
return out;
}
// setup ///////////////////////////////////////////////////////////////////////
template <typename FImpl>
void TSmear<FImpl>::setup(void)
{
unsigned int nt = env().getDim(Tp);
unsigned int size = nt * sizeof(SitePropagator);
env().registerObject(getName(), size);
}
// execution ///////////////////////////////////////////////////////////////////
template <typename FImpl>
void TSmear<FImpl>::execute(void)
{
LOG(Message) << "Sink smearing propagator '" << par().q
<< "' using sink function '" << par().sink << "'."
<< std::endl;
SinkFn &sink = *env().template getObject<SinkFn>(par().sink);
PropagatorField &q = *env().template getObject<PropagatorField>(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_

2
extras/Hadrons/Modules/MSource/Point.hpp
View File

@ -119,7 +119,7 @@ template <typename FImpl>
void TPoint<FImpl>::execute(void) void TPoint<FImpl>::execute(void)
{ {
std::vector<int> position = strToVec<int>(par().position); std::vector<int> position = strToVec<int>(par().position);
typename SitePropagator::scalar_object id; SitePropagator id;
LOG(Message) << "Creating point source at position [" << par().position LOG(Message) << "Creating point source at position [" << par().position
<< "]" << std::endl; << "]" << std::endl;

158
extras/Hadrons/Modules/MSource/SeqConserved.hpp Normal file
View File

@ -0,0 +1,158 @@
/*************************************************************************************
Grid physics library, www.github.com/paboyle/Grid
Source file: extras/Hadrons/Modules/MContraction/SeqConserved.hpp
Copyright (C) 2017
Author: Andrew Lawson <andrew.lawson1991@gmail.com>
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_MSource_SeqConserved_hpp_
#define Hadrons_MSource_SeqConserved_hpp_
#include <Grid/Hadrons/Global.hpp>
#include <Grid/Hadrons/Module.hpp>
#include <Grid/Hadrons/ModuleFactory.hpp>
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 <typename FImpl>
class TSeqConserved: public Module<SeqConservedPar>
{
public:
FERM_TYPE_ALIASES(FImpl,);
public:
// constructor
TSeqConserved(const std::string name);
// destructor
virtual ~TSeqConserved(void) = default;
// dependency relation
virtual std::vector<std::string> getInput(void);
virtual std::vector<std::string> getOutput(void);
// setup
virtual void setup(void);
// execution
virtual void execute(void);
};
MODULE_REGISTER_NS(SeqConserved, TSeqConserved<FIMPL>, MSource);
/******************************************************************************
* TSeqConserved implementation *
******************************************************************************/
// constructor /////////////////////////////////////////////////////////////////
template <typename FImpl>
TSeqConserved<FImpl>::TSeqConserved(const std::string name)
: Module<SeqConservedPar>(name)
{}
// dependencies/products ///////////////////////////////////////////////////////
template <typename FImpl>
std::vector<std::string> TSeqConserved<FImpl>::getInput(void)
{
std::vector<std::string> in = {par().q, par().action};
return in;
}
template <typename FImpl>
std::vector<std::string> TSeqConserved<FImpl>::getOutput(void)
{
std::vector<std::string> out = {getName()};
return out;
}
// setup ///////////////////////////////////////////////////////////////////////
template <typename FImpl>
void TSeqConserved<FImpl>::setup(void)
{
auto Ls_ = env().getObjectLs(par().action);
env().template registerLattice<PropagatorField>(getName(), Ls_);
}
// execution ///////////////////////////////////////////////////////////////////
template <typename FImpl>
void TSeqConserved<FImpl>::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<PropagatorField>(getName());
PropagatorField &q = *env().template getObject<PropagatorField>(par().q);
FMat &mat = *(env().template getObject<FMat>(par().action));
std::vector<Real> mom = strToVec<Real>(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_

183
extras/Hadrons/Modules/MUtilities/TestSeqConserved.hpp Normal file
View File

@ -0,0 +1,183 @@
/*************************************************************************************
Grid physics library, www.github.com/paboyle/Grid
Source file: extras/Hadrons/Modules/MUtilities/TestSeqConserved.hpp
Copyright (C) 2017
Author: Andrew Lawson <andrew.lawson1991@gmail.com>
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_TestSeqConserved_hpp_
#define Hadrons_MUtilities_TestSeqConserved_hpp_
#include <Grid/Hadrons/Global.hpp>
#include <Grid/Hadrons/Module.hpp>
#include <Grid/Hadrons/ModuleFactory.hpp>
BEGIN_HADRONS_NAMESPACE
/*
Ward Identity contractions using sequential propagators.
-----------------------------
* options:
- q: point source propagator, 5D if available (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, qSeq,
std::string, action,
std::string, origin,
unsigned int, t_J,
unsigned int, mu,
Current, curr);
};
template <typename FImpl>
class TTestSeqConserved: public Module<TestSeqConservedPar>
{
public:
FERM_TYPE_ALIASES(FImpl,);
public:
// constructor
TTestSeqConserved(const std::string name);
// destructor
virtual ~TTestSeqConserved(void) = default;
// dependency relation
virtual std::vector<std::string> getInput(void);
virtual std::vector<std::string> getOutput(void);
// setup
virtual void setup(void);
// execution
virtual void execute(void);
};
MODULE_REGISTER_NS(TestSeqConserved, TTestSeqConserved<FIMPL>, MUtilities);
/******************************************************************************
* TTestSeqConserved implementation *
******************************************************************************/
// constructor /////////////////////////////////////////////////////////////////
template <typename FImpl>
TTestSeqConserved<FImpl>::TTestSeqConserved(const std::string name)
: Module<TestSeqConservedPar>(name)
{}
// dependencies/products ///////////////////////////////////////////////////////
template <typename FImpl>
std::vector<std::string> TTestSeqConserved<FImpl>::getInput(void)
{
std::vector<std::string> in = {par().q, par().qSeq, par().action};
return in;
}
template <typename FImpl>
std::vector<std::string> TTestSeqConserved<FImpl>::getOutput(void)
{
std::vector<std::string> out = {getName()};
return out;
}
// setup ///////////////////////////////////////////////////////////////////////
template <typename FImpl>
void TTestSeqConserved<FImpl>::setup(void)
{
auto Ls = env().getObjectLs(par().q);
if (Ls != env().getObjectLs(par().action))
{
HADRON_ERROR("Ls mismatch between quark action and propagator");
}
}
// execution ///////////////////////////////////////////////////////////////////
template <typename FImpl>
void TTestSeqConserved<FImpl>::execute(void)
{
PropagatorField tmp(env().getGrid());
PropagatorField &q = *env().template getObject<PropagatorField>(par().q);
PropagatorField &qSeq = *env().template getObject<PropagatorField>(par().qSeq);
FMat &act = *(env().template getObject<FMat>(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;
Complex test_S, test_V, check_S, check_V;
std::vector<TComplex> 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.
std::vector<int> siteCoord;
siteCoord = strToVec<int>(par().origin);
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*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;
LOG(Message) << "Check S = " << abs(check_S) << std::endl;
LOG(Message) << "Check V = " << abs(check_V) << std::endl;
// Check difference = 0
check_S -= test_S;
check_V -= test_V;
LOG(Message) << "Consistency check for sequential conserved "
<< par().curr << " current insertion: " << std::endl;
LOG(Message) << "Diff S = " << abs(check_S) << std::endl;
LOG(Message) << "Diff V = " << abs(check_V) << std::endl;
}
END_MODULE_NAMESPACE
END_HADRONS_NAMESPACE
#endif // Hadrons_TestSeqConserved_hpp_

147
extras/Hadrons/Modules/MUtilities/TestSeqGamma.hpp Normal file
View File

@ -0,0 +1,147 @@
/*************************************************************************************
Grid physics library, www.github.com/paboyle/Grid
Source file: extras/Hadrons/Modules/MUtilities/TestSeqGamma.hpp
Copyright (C) 2017
Author: Andrew Lawson <andrew.lawson1991@gmail.com>
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_
#include <Grid/Hadrons/Global.hpp>
#include <Grid/Hadrons/Module.hpp>
#include <Grid/Hadrons/ModuleFactory.hpp>
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 <typename FImpl>
class TTestSeqGamma: public Module<TestSeqGammaPar>
{
public:
FERM_TYPE_ALIASES(FImpl,);
public:
// constructor
TTestSeqGamma(const std::string name);
// destructor
virtual ~TTestSeqGamma(void) = default;
// dependency relation
virtual std::vector<std::string> getInput(void);
virtual std::vector<std::string> getOutput(void);
// setup
virtual void setup(void);
// execution
virtual void execute(void);
};
MODULE_REGISTER_NS(TestSeqGamma, TTestSeqGamma<FIMPL>, MUtilities);
/******************************************************************************
* TTestSeqGamma implementation *
******************************************************************************/
// constructor /////////////////////////////////////////////////////////////////
template <typename FImpl>
TTestSeqGamma<FImpl>::TTestSeqGamma(const std::string name)
: Module<TestSeqGammaPar>(name)
{}
// dependencies/products ///////////////////////////////////////////////////////
template <typename FImpl>
std::vector<std::string> TTestSeqGamma<FImpl>::getInput(void)
{
std::vector<std::string> in = {par().q, par().qSeq};
return in;
}
template <typename FImpl>
std::vector<std::string> TTestSeqGamma<FImpl>::getOutput(void)
{
std::vector<std::string> out = {getName()};
return out;
}
// setup ///////////////////////////////////////////////////////////////////////
template <typename FImpl>
void TTestSeqGamma<FImpl>::setup(void)
{
}
// execution ///////////////////////////////////////////////////////////////////
template <typename FImpl>
void TTestSeqGamma<FImpl>::execute(void)
{
PropagatorField &q = *env().template getObject<PropagatorField>(par().q);
PropagatorField &qSeq = *env().template getObject<PropagatorField>(par().qSeq);
LatticeComplex c(env().getGrid());
Gamma g5(Gamma::Algebra::Gamma5);
Gamma g(par().gamma);
SitePropagator qSite;
Complex test, check;
std::vector<TComplex> 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<int> siteCoord;
siteCoord = strToVec<int>(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_

8
extras/Hadrons/modules.inc
View File

@ -16,6 +16,7 @@ modules_hpp =\
Modules/MContraction/DiscLoop.hpp \ Modules/MContraction/DiscLoop.hpp \
Modules/MContraction/Gamma3pt.hpp \ Modules/MContraction/Gamma3pt.hpp \
Modules/MContraction/Meson.hpp \ Modules/MContraction/Meson.hpp \
Modules/MContraction/WardIdentity.hpp \
Modules/MContraction/WeakHamiltonian.hpp \ Modules/MContraction/WeakHamiltonian.hpp \
Modules/MContraction/WeakHamiltonianEye.hpp \ Modules/MContraction/WeakHamiltonianEye.hpp \
Modules/MContraction/WeakHamiltonianNonEye.hpp \ Modules/MContraction/WeakHamiltonianNonEye.hpp \
@ -30,9 +31,12 @@ modules_hpp =\
Modules/MScalar/FreeProp.hpp \ Modules/MScalar/FreeProp.hpp \
Modules/MScalar/Scalar.hpp \ Modules/MScalar/Scalar.hpp \
Modules/MSink/Point.hpp \ Modules/MSink/Point.hpp \
Modules/MSink/Smear.hpp \
Modules/MSolver/RBPrecCG.hpp \ Modules/MSolver/RBPrecCG.hpp \
Modules/MSource/Point.hpp \ Modules/MSource/Point.hpp \
Modules/MSource/SeqConserved.hpp \
Modules/MSource/SeqGamma.hpp \ Modules/MSource/SeqGamma.hpp \
Modules/MSource/Wall.hpp \ Modules/MSource/Wall.hpp \
Modules/MSource/Z2.hpp Modules/MSource/Z2.hpp \
Modules/MUtilities/TestSeqConserved.hpp \
Modules/MUtilities/TestSeqGamma.hpp

8
lib/qcd/QCD.h
View File

@ -492,6 +492,14 @@ namespace QCD {
return traceIndex<ColourIndex>(lhs); return traceIndex<ColourIndex>(lhs);
} }
//////////////////////////////////////////
// Current types
//////////////////////////////////////////
GRID_SERIALIZABLE_ENUM(Current, undef,
Vector, 0,
Axial, 1,
Tadpole, 2);
} //namespace QCD } //namespace QCD
} // Grid } // Grid

15
lib/qcd/action/fermion/FermionOperator.h
View File

@ -112,6 +112,21 @@ namespace Grid {
/////////////////////////////////////////////// ///////////////////////////////////////////////
virtual void ImportGauge(const GaugeField & _U)=0; 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<Real> mom,
unsigned int tmin,
unsigned int tmax)=0;
}; };
} }

29
lib/qcd/action/fermion/FermionOperatorImpl.h
View File

@ -212,6 +212,13 @@ namespace QCD {
StencilImpl &St) { StencilImpl &St) {
mult(&phi(), &U(mu), &chi()); mult(&phi(), &U(mu), &chi());
} }
inline void multLinkProp(SitePropagator &phi,
const SiteDoubledGaugeField &U,
const SitePropagator &chi,
int mu) {
mult(&phi(), &U(mu), &chi());
}
template <class ref> template <class ref>
inline void loadLinkElement(Simd &reg, ref &memory) { inline void loadLinkElement(Simd &reg, ref &memory) {
@ -340,7 +347,20 @@ class DomainWallVec5dImpl : public PeriodicGaugeImpl< GaugeImplTypes< S,Nrepres
} }
mult(&phi(), &UU(), &chi()); 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) inline void DoubleStore(GridBase *GaugeGrid, DoubledGaugeField &Uds,const GaugeField &Umu)
{ {
SiteScalarGaugeField ScalarUmu; SiteScalarGaugeField ScalarUmu;
@ -538,6 +558,13 @@ class GparityWilsonImpl : public ConjugateGaugeImpl<GaugeImplTypes<S, Nrepresent
} }
// Fixme: Gparity prop * link
inline void multLinkProp(SitePropagator &phi, const SiteDoubledGaugeField &U,
const SitePropagator &chi, int mu)
{
assert(0);
}
inline void DoubleStore(GridBase *GaugeGrid,DoubledGaugeField &Uds,const GaugeField &Umu) inline void DoubleStore(GridBase *GaugeGrid,DoubledGaugeField &Uds,const GaugeField &Umu)
{ {
conformable(Uds._grid,GaugeGrid); conformable(Uds._grid,GaugeGrid);

25
lib/qcd/action/fermion/ImprovedStaggeredFermion.cc
View File

@ -393,6 +393,31 @@ void ImprovedStaggeredFermion<Impl>::DhopInternal(StencilImpl &st, LebesgueOrder
} }
}; };
////////////////////////////////////////////////////////
// Conserved current - not yet implemented.
////////////////////////////////////////////////////////
template <class Impl>
void ImprovedStaggeredFermion<Impl>::ContractConservedCurrent(PropagatorField &q_in_1,
PropagatorField &q_in_2,
PropagatorField &q_out,
Current curr_type,
unsigned int mu)
{
assert(0);
}
template <class Impl>
void ImprovedStaggeredFermion<Impl>::SeqConservedCurrent(PropagatorField &q_in,
PropagatorField &q_out,
Current curr_type,
unsigned int mu,
std::vector<Real> mom,
unsigned int tmin,
unsigned int tmax)
{
assert(0);
}
FermOpStaggeredTemplateInstantiate(ImprovedStaggeredFermion); FermOpStaggeredTemplateInstantiate(ImprovedStaggeredFermion);
//AdjointFermOpTemplateInstantiate(ImprovedStaggeredFermion); //AdjointFermOpTemplateInstantiate(ImprovedStaggeredFermion);

16
lib/qcd/action/fermion/ImprovedStaggeredFermion.h
View File

@ -157,6 +157,22 @@ class ImprovedStaggeredFermion : public StaggeredKernels<Impl>, public ImprovedS
LebesgueOrder Lebesgue; LebesgueOrder Lebesgue;
LebesgueOrder LebesgueEvenOdd; 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<Real> mom,
unsigned int tmin,
unsigned int tmax);
}; };
typedef ImprovedStaggeredFermion<StaggeredImplF> ImprovedStaggeredFermionF; typedef ImprovedStaggeredFermion<StaggeredImplF> ImprovedStaggeredFermionF;

24
lib/qcd/action/fermion/ImprovedStaggeredFermion5D.cc
View File

@ -345,6 +345,30 @@ void ImprovedStaggeredFermion5D<Impl>::MooeeInvDag(const FermionField &in,
MooeeInv(in, out); MooeeInv(in, out);
} }
////////////////////////////////////////////////////////
// Conserved current - not yet implemented.
////////////////////////////////////////////////////////
template <class Impl>
void ImprovedStaggeredFermion5D<Impl>::ContractConservedCurrent(PropagatorField &q_in_1,
PropagatorField &q_in_2,
PropagatorField &q_out,
Current curr_type,
unsigned int mu)
{
assert(0);
}
template <class Impl>
void ImprovedStaggeredFermion5D<Impl>::SeqConservedCurrent(PropagatorField &q_in,
PropagatorField &q_out,
Current curr_type,
unsigned int mu,
std::vector<Real> mom,
unsigned int tmin,
unsigned int tmax)
{
assert(0);
}
FermOpStaggeredTemplateInstantiate(ImprovedStaggeredFermion5D); FermOpStaggeredTemplateInstantiate(ImprovedStaggeredFermion5D);
FermOpStaggeredVec5dTemplateInstantiate(ImprovedStaggeredFermion5D); FermOpStaggeredVec5dTemplateInstantiate(ImprovedStaggeredFermion5D);

15
lib/qcd/action/fermion/ImprovedStaggeredFermion5D.h
View File

@ -160,6 +160,21 @@ namespace QCD {
// Comms buffer // Comms buffer
std::vector<SiteHalfSpinor,alignedAllocator<SiteHalfSpinor> > comm_buf; std::vector<SiteHalfSpinor,alignedAllocator<SiteHalfSpinor> > 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<Real> mom,
unsigned int tmin,
unsigned int tmax);
}; };
}} }}

100
lib/qcd/action/fermion/WilsonFermion.cc
View File

@ -345,6 +345,106 @@ void WilsonFermion<Impl>::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 <class Impl>
void WilsonFermion<Impl>::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);
PropagatorField tmp1(_grid), tmp2(_grid);
q_out = zero;
// 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(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);
}
}
template <class Impl>
void WilsonFermion<Impl>::SeqConservedCurrent(PropagatorField &q_in,
PropagatorField &q_out,
Current curr_type,
unsigned int mu,
std::vector<Real> mom,
unsigned int tmin,
unsigned int tmax)
{
conformable(_grid, q_in._grid);
conformable(_grid, q_out._grid);
Lattice<iSinglet<Simd>> 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;
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);
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); FermOpTemplateInstantiate(WilsonFermion);
AdjointFermOpTemplateInstantiate(WilsonFermion); AdjointFermOpTemplateInstantiate(WilsonFermion);
TwoIndexFermOpTemplateInstantiate(WilsonFermion); TwoIndexFermOpTemplateInstantiate(WilsonFermion);

16
lib/qcd/action/fermion/WilsonFermion.h
View File

@ -146,6 +146,22 @@ class WilsonFermion : public WilsonKernels<Impl>, public WilsonFermionStatic {
LebesgueOrder Lebesgue; LebesgueOrder Lebesgue;
LebesgueOrder LebesgueEvenOdd; 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<Real> mom,
unsigned int tmin,
unsigned int tmax);
}; };
typedef WilsonFermion<WilsonImplF> WilsonFermionF; typedef WilsonFermion<WilsonImplF> WilsonFermionF;

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

@ -12,6 +12,7 @@ Author: Peter Boyle <paboyle@ph.ed.ac.uk>
Author: Peter Boyle <peterboyle@Peters-MacBook-Pro-2.local> Author: Peter Boyle <peterboyle@Peters-MacBook-Pro-2.local>
Author: paboyle <paboyle@ph.ed.ac.uk> Author: paboyle <paboyle@ph.ed.ac.uk>
Author: Guido Cossu <guido.cossu@ed.ac.uk> Author: Guido Cossu <guido.cossu@ed.ac.uk>
Author: Andrew Lawson <andrew.lawson1991@gmail.com>
This program is free software; you can redistribute it and/or modify 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 it under the terms of the GNU General Public License as published by
@ -671,6 +672,162 @@ void WilsonFermion5D<Impl>::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.
******************************************************************************/
// Helper macro to reverse Simd vector. Fixme: slow, generic implementation.
#define REVERSE_LS(qSite, qSiteRev, Nsimd) \
{ \
std::vector<typename SitePropagator::scalar_object> 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 <class Impl>
void WilsonFermion5D<Impl>::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 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),
// 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 * LLs;
unsigned int sF2 = (sU + 1) * LLs - 1;
for (unsigned int s = 0; s < LLs; ++s)
{
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],
qSite2,
q_out._odata[sU],
Umu, sU, mu, axial_sign);
Kernels::ContractConservedCurrentSiteBwd(q_in_1._odata[sF1],
qmuSite2,
q_out._odata[sU],
Umu, sU, mu, axial_sign);
sF1++;
sF2--;
}
}
}
template <class Impl>
void WilsonFermion5D<Impl>::SeqConservedCurrent(PropagatorField &q_in,
PropagatorField &q_out,
Current curr_type,
unsigned int mu,
std::vector<Real> mom,
unsigned int tmin,
unsigned int tmax)
{
conformable(q_in._grid, FermionGrid());
conformable(q_in._grid, q_out._grid);
Lattice<iSinglet<Simd>> ph(FermionGrid()), coor(FermionGrid());
PropagatorField tmpFwd(FermionGrid()), tmpBwd(FermionGrid()),
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;
for(unsigned int nu = 0; nu < Nd - 1; 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);
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)
{
// 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)
{
unsigned int sF = sU * LLs;
for (unsigned int s = 0; s < LLs; ++s)
{
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);
++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 * LLs;
for (unsigned int s = 0; s < LLs; ++s)
{
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);
++sF;
}
}
}
}
FermOpTemplateInstantiate(WilsonFermion5D); FermOpTemplateInstantiate(WilsonFermion5D);
GparityFermOpTemplateInstantiate(WilsonFermion5D); GparityFermOpTemplateInstantiate(WilsonFermion5D);

15
lib/qcd/action/fermion/WilsonFermion5D.h
View File

@ -214,6 +214,21 @@ namespace QCD {
// Comms buffer // Comms buffer
std::vector<SiteHalfSpinor,alignedAllocator<SiteHalfSpinor> > comm_buf; std::vector<SiteHalfSpinor,alignedAllocator<SiteHalfSpinor> > 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<Real> mom,
unsigned int tmin,
unsigned int tmax);
}; };
}} }}

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

@ -281,6 +281,172 @@ void WilsonKernels<Impl>::DhopDir( StencilImpl &st, DoubledGaugeField &U,SiteHal
vstream(out._odata[sF], result); 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.
******************************************************************************/
// 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))
/*******************************************************************************
* 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<class Impl>
void WilsonKernels<Impl>::ContractConservedCurrentSiteFwd(
const SitePropagator &q_in_1,
const SitePropagator &q_in_2,
SitePropagator &q_out,
DoubledGaugeField &U,
unsigned int sU,
unsigned int mu,
bool switch_sign)
{
SitePropagator result, tmp;
Gamma g5(Gamma::Algebra::Gamma5);
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 -= result;
}
else
{
q_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<class Impl>
void WilsonKernels<Impl>::ContractConservedCurrentSiteBwd(
const SitePropagator &q_in_1,
const SitePropagator &q_in_2,
SitePropagator &q_out,
DoubledGaugeField &U,
unsigned int sU,
unsigned int mu,
bool switch_sign)
{
SitePropagator result, tmp;
Gamma g5(Gamma::Algebra::Gamma5);
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 += result;
}
else
{
q_out -= result;
}
}
// G-parity requires more specialised implementation.
#define NO_CURR_SITE(Impl) \
template <> \
void WilsonKernels<Impl>::ContractConservedCurrentSiteFwd( \
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<Impl>::ContractConservedCurrentSiteBwd( \
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); \
}
NO_CURR_SITE(GparityWilsonImplF);
NO_CURR_SITE(GparityWilsonImplD);
NO_CURR_SITE(GparityWilsonImplFH);
NO_CURR_SITE(GparityWilsonImplDF);
/*******************************************************************************
* 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<class Impl>
void WilsonKernels<Impl>::SeqConservedCurrentSiteFwd(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);
result = WilsonCurrentFwd(result, mu);
// Zero any unwanted timeslice entries.
result = predicatedWhere(t_mask, result, 0.*result);
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<class Impl>
void WilsonKernels<Impl>::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);
// Zero any unwanted timeslice entries.
result = predicatedWhere(t_mask, result, 0.*result);
if (switch_sign)
{
q_out += result;
}
else
{
q_out -= result;
}
}
FermOpTemplateInstantiate(WilsonKernels); FermOpTemplateInstantiate(WilsonKernels);
AdjointFermOpTemplateInstantiate(WilsonKernels); AdjointFermOpTemplateInstantiate(WilsonKernels);
TwoIndexFermOpTemplateInstantiate(WilsonKernels); TwoIndexFermOpTemplateInstantiate(WilsonKernels);

32
lib/qcd/action/fermion/WilsonKernels.h
View File

@ -180,6 +180,38 @@ public:
void DhopDir(StencilImpl &st, DoubledGaugeField &U,SiteHalfSpinor * buf, void DhopDir(StencilImpl &st, DoubledGaugeField &U,SiteHalfSpinor * buf,
int sF, int sU, const FermionField &in, FermionField &out, int dirdisp, int gamma); int sF, int sU, const FermionField &in, FermionField &out, int dirdisp, int gamma);
//////////////////////////////////////////////////////////////////////////////
// Utilities for inserting Wilson conserved current.
//////////////////////////////////////////////////////////////////////////////
void ContractConservedCurrentSiteFwd(const SitePropagator &q_in_1,
const SitePropagator &q_in_2,
SitePropagator &q_out,
DoubledGaugeField &U,
unsigned int sU,
unsigned int mu,
bool switch_sign = false);
void ContractConservedCurrentSiteBwd(const SitePropagator &q_in_1,
const SitePropagator &q_in_2,
SitePropagator &q_out,
DoubledGaugeField &U,
unsigned int sU,
unsigned int mu,
bool switch_sign = false);
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: private:
// Specialised variants // Specialised variants
void GenericDhopSite(StencilImpl &st, LebesgueOrder &lo, DoubledGaugeField &U, SiteHalfSpinor * buf, void GenericDhopSite(StencilImpl &st, LebesgueOrder &lo, DoubledGaugeField &U, SiteHalfSpinor * buf,

25
lib/simd/Grid_avx.h
View File

@ -701,9 +701,28 @@ namespace Optimization {
//Integer Reduce //Integer Reduce
template<> template<>
inline Integer Reduce<Integer, __m256i>::operator()(__m256i in){ inline Integer Reduce<Integer, __m256i>::operator()(__m256i in){
// FIXME unimplemented __m128i ret;
printf("Reduce : Missing integer implementation -> FIX\n"); #if defined (AVX2)
assert(0); // 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 = _mm_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);
} }
} }

22
lib/simd/Grid_avx512.h
View File

@ -543,6 +543,24 @@ namespace Optimization {
u512d conv; conv.v = v1; u512d conv; conv.v = v1;
return conv.f[0]; return conv.f[0];
} }
//Integer Reduce
template<>
inline Integer Reduce<Integer, __m512i>::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 = _mm_add_epi32(u1, u2);
return _mm_cvtsi128_si32(ret);
}
#else #else
//Complex float Reduce //Complex float Reduce
template<> template<>
@ -570,9 +588,7 @@ namespace Optimization {
//Integer Reduce //Integer Reduce
template<> template<>
inline Integer Reduce<Integer, __m512i>::operator()(__m512i in){ inline Integer Reduce<Integer, __m512i>::operator()(__m512i in){
// FIXME unimplemented return _mm512_reduce_add_epi32(in);
printf("Reduce : Missing integer implementation -> FIX\n");
assert(0);
} }
#endif #endif

4
lib/simd/Grid_imci.h
View File

@ -401,9 +401,7 @@ namespace Optimization {
//Integer Reduce //Integer Reduce
template<> template<>
inline Integer Reduce<Integer, __m512i>::operator()(__m512i in){ inline Integer Reduce<Integer, __m512i>::operator()(__m512i in){
// FIXME unimplemented return _mm512_reduce_add_epi32(in);
printf("Reduce : Missing integer implementation -> FIX\n");
assert(0);
} }

90
lib/simd/Grid_qpx.h
View File

@ -374,6 +374,84 @@ namespace Optimization {
// Complex float // Complex float
FLOAT_WRAP_2(operator(), inline) 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
#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{ struct Permute{
//Complex double //Complex double
@ -497,15 +575,19 @@ namespace Optimization {
//Integer Reduce //Integer Reduce
template<> template<>
inline Integer Reduce<Integer, int>::operator()(int in){ inline Integer Reduce<Integer, veci>::operator()(veci in){
// FIXME unimplemented Integer a = 0;
printf("Reduce : Missing integer implementation -> FIX\n"); for (unsigned int i = 0; i < W<Integer>::r; ++i)
assert(0); {
a += in.v[i];
}
return a;
} }
} }
//////////////////////////////////////////////////////////////////////////////// ////////////////////////////////////////////////////////////////////////////////
// Here assign types // Here assign types
typedef Optimization::vech SIMD_Htype; // Half precision type
typedef Optimization::vector4float SIMD_Ftype; // Single precision type typedef Optimization::vector4float SIMD_Ftype; // Single precision type
typedef vector4double SIMD_Dtype; // Double precision type typedef vector4double SIMD_Dtype; // Double precision type
typedef Optimization::veci SIMD_Itype; // Integer type typedef Optimization::veci SIMD_Itype; // Integer type

6
lib/simd/Grid_sse4.h
View File

@ -570,9 +570,9 @@ namespace Optimization {
//Integer Reduce //Integer Reduce
template<> template<>
inline Integer Reduce<Integer, __m128i>::operator()(__m128i in){ inline Integer Reduce<Integer, __m128i>::operator()(__m128i in){
// FIXME unimplemented __m128i v1 = _mm_hadd_epi32(in, in);
printf("Reduce : Missing integer implementation -> FIX\n"); __m128i v2 = _mm_hadd_epi32(v1, v1);
assert(0); return _mm_cvtsi128_si32(v2);
} }
} }

47
tests/Test_simd.cc
View File

@ -183,8 +183,6 @@ void IntTester(const functor &func)
{ {
typedef Integer scal; typedef Integer scal;
typedef vInteger vec; typedef vInteger vec;
GridSerialRNG sRNG;
sRNG.SeedFixedIntegers(std::vector<int>({45,12,81,9}));
int Nsimd = vec::Nsimd(); int Nsimd = vec::Nsimd();
@ -287,6 +285,50 @@ void ReductionTester(const functor &func)
} }
template<class reduced,class scal, class vec,class functor >
void IntReductionTester(const functor &func)
{
int Nsimd = vec::Nsimd();
std::vector<scal> input1(Nsimd);
std::vector<scal> input2(Nsimd);
reduced result(0);
reduced reference(0);
reduced tmp;
std::vector<vec,alignedAllocator<vec> > buf(3);
vec & v_input1 = buf[0];
vec & v_input2 = buf[1];
for(int i=0;i<Nsimd;i++){
input1[i] = (i + 1) * 30;
input2[i] = (i + 1) * 20;
}
merge<vec,scal>(v_input1,input1);
merge<vec,scal>(v_input2,input2);
func.template vfunc<reduced,vec>(result,v_input1,v_input2);
for(int i=0;i<Nsimd;i++) {
func.template sfunc<reduced,scal>(tmp,input1[i],input2[i]);
reference+=tmp;
}
std::cout<<GridLogMessage << " " << func.name()<<std::endl;
int ok=0;
if ( reference-result != 0 ){
std::cout<<GridLogMessage<< "*****" << std::endl;
std::cout<<GridLogMessage<< reference-result << " " <<reference<< " " << result<<std::endl;
ok++;
}
if ( ok==0 ) {
std::cout<<GridLogMessage << " OK!" <<std::endl;
}
assert(ok==0);
}
class funcPermute { class funcPermute {
public: public:
@ -691,6 +733,7 @@ int main (int argc, char ** argv)
IntTester(funcPlus()); IntTester(funcPlus());
IntTester(funcMinus()); IntTester(funcMinus());
IntTester(funcTimes()); IntTester(funcTimes());
IntReductionTester<Integer, Integer, vInteger>(funcReduce());
std::cout<<GridLogMessage << "==================================="<< std::endl; std::cout<<GridLogMessage << "==================================="<< std::endl;
std::cout<<GridLogMessage << "Testing precisionChange "<< std::endl; std::cout<<GridLogMessage << "Testing precisionChange "<< std::endl;

659
tests/hadrons/Test_hadrons.hpp Normal file
View File

@ -0,0 +1,659 @@
/*******************************************************************************
Grid physics library, www.github.com/paboyle/Grid
Source file: tests/hadrons/Test_hadrons.hpp
Copyright (C) 2017
Author: Andrew Lawson <andrew.lawson1991@gmail.com>
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 <Grid/Hadrons/Application.hpp>
using namespace Grid;
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))
#define ADD_INDEX(s, n) (s + "_" + std::to_string(n))
#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)
#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);\
makeWallSource(application, srcName, tW, mom);\
makePropagator(application, propName, srcName, solver);\
}
#define MAKE_WALL_PROP(tW, propName, solver)\
MAKE_3MOM_WALL_PROP(tW, ZERO_MOM, propName, solver)
// Sequential source macros
#define MAKE_SEQUENTIAL_PROP(tS, qSrc, mom, seqPropName, solver, gamma)\
{\
std::string srcName = seqPropName + "_src";\
makeSequentialSource(application, srcName, qSrc, tS, gamma, mom);\
makePropagator(application, seqPropName, srcName, solver);\
}
// Point source macros
#define MAKE_POINT_PROP(pos, propName, solver)\
{\
std::string srcName = NAME_POINT_SOURCE(pos);\
makePointSource(application, srcName, pos);\
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.
* 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 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<MAction::Wilson>(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.
* 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, std::string boundary = "1 1 1 -1")
{
if (!(Environment::getInstance().hasModule(actionName)))
{
MAction::DWF::Par actionPar;
actionPar.gauge = gaugeField;
actionPar.Ls = Ls;
actionPar.M5 = M5;
actionPar.mass = mass;
actionPar.boundary = boundary;
application.createModule<MAction::DWF>(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<MSolver::RBPrecCG>(solverName,
solverPar);
}
}
/*******************************************************************************
* Name: makePointSource
* Purpose: Construct point source and add to application module.
* Parameters: application - main application that stores modules.
* srcName - name of source module to create.
* pos - Position of point source.
* Returns: None.
******************************************************************************/
inline void makePointSource(Application &application, std::string srcName,
std::string pos)
{
// If the source already exists, don't make the module again.
if (!(Environment::getInstance().hasModule(srcName)))
{
MSource::Point::Par pointPar;
pointPar.position = pos;
application.createModule<MSource::Point>(srcName, pointPar);
}
}
/*******************************************************************************
* Name: makeSequentialSource
* Purpose: Construct sequential source 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.
* tS - sequential source timeslice.
* mom - momentum insertion (default is zero).
* Returns: None.
******************************************************************************/
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.
if (!(Environment::getInstance().hasModule(srcName)))
{
MSource::SeqGamma::Par seqPar;
seqPar.q = qSrc;
seqPar.tA = tS;
seqPar.tB = tS;
seqPar.mom = mom;
seqPar.gamma = gamma;
application.createModule<MSource::SeqGamma>(srcName, seqPar);
}
}
/*******************************************************************************
* 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<MSource::SeqConserved>(srcName, seqPar);
}
}
/*******************************************************************************
* 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<MSource::Z2>(srcName, noisePar);
}
}
/*******************************************************************************
* Name: makeWallSource
* Purpose: Construct wall source and add to application module.
* Parameters: application - main application that stores modules.
* srcName - name of source module to create.
* tW - wall source timeslice.
* mom - momentum insertion (default is zero).
* Returns: None.
******************************************************************************/
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.
if (!(Environment::getInstance().hasModule(srcName)))
{
MSource::Wall::Par wallPar;
wallPar.tW = tW;
wallPar.mom = mom;
application.createModule<MSource::Wall>(srcName, wallPar);
}
}
/*******************************************************************************
* Name: makePointSink
* Purpose: Create function for point sink smearing of a propagator.
* Parameters: application - main application that stores modules.
* propName - name of input propagator.
* sinkFnct - name of output sink smearing module.
* mom - momentum insertion (default is zero).
* Returns: None.
******************************************************************************/
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<MSink::Point>(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.
if (!(Environment::getInstance().hasModule(smearedProp)))
{
MSink::Smear::Par smearPar;
smearPar.q = propName;
smearPar.sink = sinkFnct;
application.createModule<MSink::Smear>(smearedProp, smearPar);
}
}
/*******************************************************************************
* Name: makePropagator
* Purpose: Construct source and propagator then add to application module.
* Parameters: application - main application that stores modules.
* propName - name of propagator module to create.
* srcName - name of source module to use.
* solver - solver to use (default is CG).
* Returns: None.
******************************************************************************/
inline void makePropagator(Application &application, std::string &propName,
std::string &srcName, std::string &solver)
{
// If the propagator already exists, don't make the module again.
if (!(Environment::getInstance().hasModule(propName)))
{
MFermion::GaugeProp::Par quarkPar;
quarkPar.source = srcName;
quarkPar.solver = solver;
application.createModule<MFermion::GaugeProp>(propName, quarkPar);
}
}
/*******************************************************************************
* Name: makeLoop
* Purpose: Use noise source and inversion result to make loop propagator, then
* add to application module.
* Parameters: application - main application that stores modules.
* propName - name of propagator module to create.
* srcName - name of noise source module to use.
* resName - name of inversion result on given noise source.
* Returns: None.
******************************************************************************/
inline void makeLoop(Application &application, std::string &propName,
std::string &srcName, std::string &resName)
{
// If the loop propagator already exists, don't make the module again.
if (!(Environment::getInstance().hasModule(propName)))
{
MLoop::NoiseLoop::Par loopPar;
loopPar.q = resName;
loopPar.eta = srcName;
application.createModule<MLoop::NoiseLoop>(propName, loopPar);
}
}
/*******************************************************************************
* Contraction module creation.
******************************************************************************/
/*******************************************************************************
* Name: mesonContraction
* Purpose: Create meson contraction module and add to application module.
* Parameters: application - main application that stores modules.
* modName - unique module name.
* output - name of output files.
* q1 - quark propagator 1.
* q2 - quark propagator 2.
* sink - sink smearing module.
* gammas - gamma insertions at source and sink.
* Returns: None.
******************************************************************************/
inline void mesonContraction(Application &application,
std::string &modName, std::string &output,
std::string &q1, std::string &q2,
std::string &sink,
std::string gammas = "<Gamma5 Gamma5>")
{
if (!(Environment::getInstance().hasModule(modName)))
{
MContraction::Meson::Par mesPar;
mesPar.output = output;
mesPar.q1 = q1;
mesPar.q2 = q2;
mesPar.sink = sink;
mesPar.gammas = gammas;
application.createModule<MContraction::Meson>(modName, mesPar);
}
}
/*******************************************************************************
* Name: gamma3ptContraction
* 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, 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,
unsigned int tSnk = 0,
Gamma::Algebra gamma = Gamma::Algebra::Identity)
{
std::string modName = std::to_string(npt) + "pt_" + label;
if (!(Environment::getInstance().hasModule(modName)))
{
MContraction::Gamma3pt::Par gamma3ptPar;
gamma3ptPar.output = std::to_string(npt) + "pt/" + label;
gamma3ptPar.q1 = q1;
gamma3ptPar.q2 = q2;
gamma3ptPar.q3 = q3;
gamma3ptPar.tSnk = tSnk;
gamma3ptPar.gamma = gamma;
application.createModule<MContraction::Gamma3pt>(modName, gamma3ptPar);
}
}
/*******************************************************************************
* Name: weakContraction[Eye,NonEye]
* Purpose: Create Weak Hamiltonian contraction module for Eye/NonEye topology
* and add to application module.
* Parameters: application - main application that stores modules.
* npt - specify n-point correlator (for labelling).
* q1 - quark propagator 1.
* q2 - quark propagator 2.
* 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, unsigned int tSnk = 0)\
{\
std::string modName = std::to_string(npt) + "pt_" + label;\
if (!(Environment::getInstance().hasModule(modName)))\
{\
MContraction::WeakHamiltonian##top::Par weakPar;\
weakPar.output = std::to_string(npt) + "pt/" + label;\
weakPar.q1 = q1;\
weakPar.q2 = q2;\
weakPar.q3 = q3;\
weakPar.q4 = q4;\
weakPar.tSnk = tSnk;\
application.createModule<MContraction::WeakHamiltonian##top>(modName, weakPar);\
}\
}
HW_CONTRACTION(Eye) // weakContractionEye
HW_CONTRACTION(NonEye) // weakContractionNonEye
/*******************************************************************************
* Name: disc0Contraction
* Purpose: Create contraction module for 4pt Weak Hamiltonian + current
* disconnected topology for neutral mesons and add to application
* module.
* Parameters: application - main application that stores modules.
* q1 - quark propagator 1.
* q2 - quark propagator 2.
* q3 - quark propagator 3.
* q4 - quark propagator 4.
* label - unique label to construct module name.
* Returns: None.
******************************************************************************/
inline void disc0Contraction(Application &application,
std::string &q1, std::string &q2,
std::string &q3, std::string &q4,
std::string &label)
{
std::string modName = "4pt_" + label;
if (!(Environment::getInstance().hasModule(modName)))
{
MContraction::WeakNeutral4ptDisc::Par disc0Par;
disc0Par.output = "4pt/" + label;
disc0Par.q1 = q1;
disc0Par.q2 = q2;
disc0Par.q3 = q3;
disc0Par.q4 = q4;
application.createModule<MContraction::WeakNeutral4ptDisc>(modName, disc0Par);
}
}
/*******************************************************************************
* Name: discLoopContraction
* Purpose: Create contraction module for disconnected loop and add to
* application module.
* Parameters: application - main application that stores modules.
* q_loop - loop quark propagator.
* modName - unique module name.
* gamma - gamma matrix to use in contraction.
* Returns: None.
******************************************************************************/
inline void discLoopContraction(Application &application,
std::string &q_loop, std::string &modName,
Gamma::Algebra gamma = Gamma::Algebra::Identity)
{
if (!(Environment::getInstance().hasModule(modName)))
{
MContraction::DiscLoop::Par discPar;
discPar.output = "disc/" + modName;
discPar.q_loop = q_loop;
discPar.gamma = gamma;
application.createModule<MContraction::DiscLoop>(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).
* 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, bool test_axial = false)
{
if (!(Environment::getInstance().hasModule(modName)))
{
MContraction::WardIdentity::Par wiPar;
if (Ls > 1)
{
wiPar.q = LABEL_5D(propName);
}
else
{
wiPar.q = propName;
}
wiPar.action = actionName;
wiPar.mass = mass;
wiPar.test_axial = test_axial;
application.createModule<MContraction::WardIdentity>(modName, wiPar);
}
}
/*******************************************************************************
* 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 - 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.
* Returns: None.
******************************************************************************/
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;
seqPar.q = propName;
seqPar.qSeq = seqName;
seqPar.action = actionName;
seqPar.origin = origin;
seqPar.t_J = t_J;
seqPar.mu = mu;
seqPar.curr = curr;
application.createModule<MUtilities::TestSeqConserved>(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<MUtilities::TestSeqGamma>(modName, seqPar);
}
}

122
tests/hadrons/Test_hadrons_3pt_contractions.cc Normal file
View File

@ -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 <andrew.lawson1991@gmail.com>
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<MGauge::Unit>(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;
}

150
tests/hadrons/Test_hadrons_conserved_current.cc Normal file
View File

@ -0,0 +1,150 @@
/*******************************************************************************
Grid physics library, www.github.com/paboyle/Grid
Source file: tests/hadrons/Test_hadrons_conserved_current.cc
Copyright (C) 2017
Author: Andrew Lawson <andrew.lawson1991@gmail.com>
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;
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,
perform_axial_tests);
/***************************************************************************
* 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 //////////////////////////////////////////////////////////
HADRONS_DEFAULT_INIT;
// run setup ///////////////////////////////////////////////////////////////
Application application;
double mass = 0.04;
double M5 = 1.8;
unsigned int Ls = 12;
// global parameters
HADRONS_DEFAULT_GLOBALS(application);
// gauge field
std::string gaugeField = "gauge";
application.createModule<MGauge::Unit>(gaugeField);
// Setup each action and the conserved current tests relevant to it.
std::string actionName = "DWF";
makeDWFAction(application, actionName, gaugeField, mass, M5, Ls);
setupWardIdentityTests(application, actionName, mass, Ls, true);
actionName = "Wilson";
makeWilsonAction(application, actionName, gaugeField, mass);
setupWardIdentityTests(application, actionName, mass);
// execution
application.saveParameterFile("ConservedCurrentTest.xml");
application.run();
// epilogue
LOG(Message) << "Grid is finalizing now" << std::endl;
Grid_finalize();
return EXIT_SUCCESS;
}

115
tests/hadrons/Test_hadrons_meson_conserved_3pt.cc Normal file
View File

@ -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 <andrew.lawson1991@gmail.com>
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<MGauge::Unit>(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<MSink::ScalarPoint>(sink_0, sinkPar);
sinkPar.mom = pmom;
application.createModule<MSink::ScalarPoint>(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;
}

156
tests/hadrons/Test_hadrons_quark.cc Normal file
View File

@ -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 <andrew.lawson1991@gmail.com>
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 <Grid/Hadrons/Modules/MFermion/GaugeProp.hpp>
using namespace Grid;
using namespace QCD;
using namespace Hadrons;
/*******************************************************************************
* Unit test functions within Quark module.
******************************************************************************/
// Alternative 4D & 5D projections
template<class vobj>
inline void make_4D_with_gammas(Lattice<vobj> &in_5d, Lattice<vobj> &out_4d, int Ls)
{
GridBase *_grid(out_4d._grid);
Lattice<vobj> 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<class vobj>
inline void make_5D_with_gammas(Lattice<vobj> &in_4d, Lattice<vobj> &out_5d, int Ls)
{
out_5d = zero;
Gamma G5(Gamma::Algebra::Gamma5);
GridBase *_grid(in_4d._grid);
Lattice<vobj> 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<int> latt_size = GridDefaultLatt();
std::vector<int> simd_layout = GridDefaultSimd(Nd,vComplex::Nsimd());
std::vector<int> 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<int> seeds4({1,2,3,4});
std::vector<int> 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<int> 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;
}

89
tests/hadrons/Test_hadrons_seq_gamma.cc Normal file
View File

@ -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 <andrew.lawson1991@gmail.com>
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<MGauge::Unit>(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;
}