portelli/Grid
1
0
Fork 0
mirror of https://github.com/paboyle/Grid.git synced 2024-11-10 07:55:35 +00:00
Code Releases Activity

Hadrons: most modules ported to the new interface, compiles but untested

Browse Source
This commit is contained in:
Antonin Portelli 2017-12-13 19:41:41 +00:00
parent 0887566134
commit 842754bea9
33 changed files with 504 additions and 273 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

4
extras/Hadrons/Module.hpp
View File

@ -92,8 +92,8 @@ static ns##mod##ModuleRegistrar ns##mod##ModuleRegistrarInstance;
#define envGet(type, name)\
*env().template getObject<type>(name)
#define envGetTmp(type, name)\
*env().template getObject<type>(getName() + "_tmp_" + name)
#define envGetTmp(type, var)\
type &var = *env().template getObject<type>(getName() + "_tmp_" + #var)
#define envHasType(type, name)\
env().template isObjectOfType<type>(name)

38
extras/Hadrons/Modules.hpp
View File

@ -30,31 +30,31 @@ See the full license in the file "LICENSE" in the top level distribution directo
#include <Grid/Hadrons/Modules/MAction/DWF.hpp>
#include <Grid/Hadrons/Modules/MAction/Wilson.hpp>
// #include <Grid/Hadrons/Modules/MContraction/Baryon.hpp>
// #include <Grid/Hadrons/Modules/MContraction/DiscLoop.hpp>
// #include <Grid/Hadrons/Modules/MContraction/Gamma3pt.hpp>
#include <Grid/Hadrons/Modules/MContraction/Baryon.hpp>
#include <Grid/Hadrons/Modules/MContraction/DiscLoop.hpp>
#include <Grid/Hadrons/Modules/MContraction/Gamma3pt.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/WeakHamiltonianEye.hpp>
// #include <Grid/Hadrons/Modules/MContraction/WeakHamiltonianNonEye.hpp>
// #include <Grid/Hadrons/Modules/MContraction/WeakNeutral4ptDisc.hpp>
#include <Grid/Hadrons/Modules/MContraction/WardIdentity.hpp>
#include <Grid/Hadrons/Modules/MContraction/WeakHamiltonian.hpp>
#include <Grid/Hadrons/Modules/MContraction/WeakHamiltonianEye.hpp>
#include <Grid/Hadrons/Modules/MContraction/WeakHamiltonianNonEye.hpp>
#include <Grid/Hadrons/Modules/MContraction/WeakNeutral4ptDisc.hpp>
#include <Grid/Hadrons/Modules/MFermion/GaugeProp.hpp>
// #include <Grid/Hadrons/Modules/MGauge/Load.hpp>
// #include <Grid/Hadrons/Modules/MGauge/Random.hpp>
// #include <Grid/Hadrons/Modules/MGauge/StochEm.hpp>
#include <Grid/Hadrons/Modules/MGauge/Load.hpp>
#include <Grid/Hadrons/Modules/MGauge/Random.hpp>
#include <Grid/Hadrons/Modules/MGauge/StochEm.hpp>
#include <Grid/Hadrons/Modules/MGauge/Unit.hpp>
// #include <Grid/Hadrons/Modules/MLoop/NoiseLoop.hpp>
#include <Grid/Hadrons/Modules/MLoop/NoiseLoop.hpp>
// #include <Grid/Hadrons/Modules/MScalar/ChargedProp.hpp>
// #include <Grid/Hadrons/Modules/MScalar/FreeProp.hpp>
// #include <Grid/Hadrons/Modules/MScalar/Scalar.hpp>
#include <Grid/Hadrons/Modules/MSink/Point.hpp>
// #include <Grid/Hadrons/Modules/MSink/Smear.hpp>
#include <Grid/Hadrons/Modules/MSink/Smear.hpp>
#include <Grid/Hadrons/Modules/MSolver/RBPrecCG.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/Wall.hpp>
// #include <Grid/Hadrons/Modules/MSource/Z2.hpp>
// #include <Grid/Hadrons/Modules/MUtilities/TestSeqConserved.hpp>
// #include <Grid/Hadrons/Modules/MUtilities/TestSeqGamma.hpp>
#include <Grid/Hadrons/Modules/MSource/SeqConserved.hpp>
#include <Grid/Hadrons/Modules/MSource/SeqGamma.hpp>
#include <Grid/Hadrons/Modules/MSource/Wall.hpp>
#include <Grid/Hadrons/Modules/MSource/Z2.hpp>
#include <Grid/Hadrons/Modules/MUtilities/TestSeqConserved.hpp>
#include <Grid/Hadrons/Modules/MUtilities/TestSeqGamma.hpp>

1
extras/Hadrons/Modules/MAction/DWF.hpp
View File

@ -119,6 +119,7 @@ void TDWF<FImpl>::setup(void)
<< std::endl;
LOG(Message) << "Fermion boundary conditions: " << par().boundary
<< std::endl;
env().createGrid(par().Ls);
auto &U = envGet(LatticeGaugeField, par().gauge);
auto &g4 = *env().getGrid();

1
extras/Hadrons/Modules/MAction/Wilson.hpp
View File

@ -115,6 +115,7 @@ void TWilson<FImpl>::setup(void)
<< " using gauge field '" << par().gauge << "'" << std::endl;
LOG(Message) << "Fermion boundary conditions: " << par().boundary
<< std::endl;
auto &U = envGet(LatticeGaugeField, par().gauge);
auto &grid = *env().getGrid();
auto &gridRb = *env().getRbGrid();

28
extras/Hadrons/Modules/MContraction/Baryon.hpp
View File

@ -71,8 +71,11 @@ public:
virtual ~TBaryon(void) = default;
// dependency relation
virtual std::vector<std::string> getInput(void);
virtual std::vector<std::string> getReference(void);
virtual std::vector<std::string> getOutput(void);
protected:
// setup
virtual void setup(void);
// execution
virtual void execute(void);
};
@ -97,14 +100,29 @@ std::vector<std::string> TBaryon<FImpl1, FImpl2, FImpl3>::getInput(void)
return input;
}
template <typename FImpl1, typename FImpl2, typename FImpl3>
std::vector<std::string> TBaryon<FImpl1, FImpl2, FImpl3>::getReference(void)
{
std::vector<std::string> ref = {};
return ref;
}
template <typename FImpl1, typename FImpl2, typename FImpl3>
std::vector<std::string> TBaryon<FImpl1, FImpl2, FImpl3>::getOutput(void)
{
std::vector<std::string> out = {getName()};
std::vector<std::string> out = {};
return out;
}
// setup ///////////////////////////////////////////////////////////////////////
template <typename FImpl1, typename FImpl2, typename FImpl3>
void TBaryon<FImpl1, FImpl2, FImpl3>::setup(void)
{
envTmpLat(LatticeComplex, "c");
}
// execution ///////////////////////////////////////////////////////////////////
template <typename FImpl1, typename FImpl2, typename FImpl3>
void TBaryon<FImpl1, FImpl2, FImpl3>::execute(void)
@ -114,10 +132,10 @@ void TBaryon<FImpl1, FImpl2, FImpl3>::execute(void)
<< par().q3 << "'" << std::endl;
CorrWriter writer(par().output);
PropagatorField1 &q1 = *env().template getObject<PropagatorField1>(par().q1);
PropagatorField2 &q2 = *env().template getObject<PropagatorField2>(par().q2);
PropagatorField3 &q3 = *env().template getObject<PropagatorField3>(par().q2);
LatticeComplex c(env().getGrid());
auto &q1 = envGet(PropagatorField1, par().q1);
auto &q2 = envGet(PropagatorField2, par().q2);
auto &q3 = envGet(PropagatorField3, par().q2);
envGetTmp(LatticeComplex, c);
Result result;
// FIXME: do contractions

17
extras/Hadrons/Modules/MContraction/DiscLoop.hpp
View File

@ -67,6 +67,7 @@ public:
virtual ~TDiscLoop(void) = default;
// dependency relation
virtual std::vector<std::string> getInput(void);
virtual std::vector<std::string> getReference(void);
virtual std::vector<std::string> getOutput(void);
protected:
// setup
@ -95,10 +96,18 @@ std::vector<std::string> TDiscLoop<FImpl>::getInput(void)
return in;
}
template <typename FImpl>
std::vector<std::string> TDiscLoop<FImpl>::getReference(void)
{
std::vector<std::string> out = {};
return out;
}
template <typename FImpl>
std::vector<std::string> TDiscLoop<FImpl>::getOutput(void)
{
std::vector<std::string> out = {getName()};
std::vector<std::string> out = {};
return out;
}
@ -107,7 +116,7 @@ std::vector<std::string> TDiscLoop<FImpl>::getOutput(void)
template <typename FImpl>
void TDiscLoop<FImpl>::setup(void)
{
envTmpLat(LatticeComplex, "c");
}
// execution ///////////////////////////////////////////////////////////////////
@ -119,12 +128,12 @@ void TDiscLoop<FImpl>::execute(void)
<< " insertion." << std::endl;
CorrWriter writer(par().output);
PropagatorField &q_loop = *env().template getObject<PropagatorField>(par().q_loop);
LatticeComplex c(env().getGrid());
auto &q_loop = envGet(PropagatorField, par().q_loop);
Gamma gamma(par().gamma);
std::vector<TComplex> buf;
Result result;
envGetTmp(LatticeComplex, c);
c = trace(gamma*q_loop);
sliceSum(c, buf, Tp);

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

@ -98,6 +98,7 @@ public:
virtual ~TGamma3pt(void) = default;
// dependency relation
virtual std::vector<std::string> getInput(void);
virtual std::vector<std::string> getReference(void);
virtual std::vector<std::string> getOutput(void);
protected:
// setup
@ -126,10 +127,18 @@ std::vector<std::string> TGamma3pt<FImpl1, FImpl2, FImpl3>::getInput(void)
return in;
}
template <typename FImpl1, typename FImpl2, typename FImpl3>
std::vector<std::string> TGamma3pt<FImpl1, FImpl2, FImpl3>::getReference(void)
{
std::vector<std::string> ref = {};
return ref;
}
template <typename FImpl1, typename FImpl2, typename FImpl3>
std::vector<std::string> TGamma3pt<FImpl1, FImpl2, FImpl3>::getOutput(void)
{
std::vector<std::string> out = {getName()};
std::vector<std::string> out = {};
return out;
}
@ -138,7 +147,7 @@ std::vector<std::string> TGamma3pt<FImpl1, FImpl2, FImpl3>::getOutput(void)
template <typename FImpl1, typename FImpl2, typename FImpl3>
void TGamma3pt<FImpl1, FImpl2, FImpl3>::setup(void)
{
envTmpLat(LatticeComplex, "c");
}
// execution ///////////////////////////////////////////////////////////////////
@ -153,10 +162,9 @@ void TGamma3pt<FImpl1, FImpl2, FImpl3>::execute(void)
// Initialise variables. q2 and q3 are normal propagators, q1 may be
// sink smeared.
CorrWriter writer(par().output);
SlicedPropagator1 &q1 = *env().template getObject<SlicedPropagator1>(par().q1);
PropagatorField2 &q2 = *env().template getObject<PropagatorField2>(par().q2);
PropagatorField3 &q3 = *env().template getObject<PropagatorField2>(par().q3);
LatticeComplex c(env().getGrid());
auto &q1 = envGet(SlicedPropagator1, par().q1);
auto &q2 = envGet(PropagatorField2, par().q2);
auto &q3 = envGet(PropagatorField2, par().q3);
Gamma g5(Gamma::Algebra::Gamma5);
Gamma gamma(par().gamma);
std::vector<TComplex> buf;
@ -165,6 +173,7 @@ void TGamma3pt<FImpl1, FImpl2, FImpl3>::execute(void)
// Extract relevant timeslice of sinked propagator q1, then contract &
// sum over all spacial positions of gamma insertion.
SitePropagator1 q1Snk = q1[par().tSnk];
envGetTmp(LatticeComplex, c);
c = trace(g5*q1Snk*adj(q2)*(g5*gamma)*q3);
sliceSum(c, buf, Tp);

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

@ -161,6 +161,7 @@ void TMeson<FImpl1, FImpl2>::parseGammaString(std::vector<GammaPair> &gammaList)
// Parse individual contractions from input string.
gammaList = strToVec<GammaPair>(par().gammas);
}
envTmpLat(LatticeComplex, "c");
}
// execution ///////////////////////////////////////////////////////////////////
@ -192,8 +193,8 @@ void TMeson<FImpl1, FImpl2>::execute(void)
if (envHasType(SlicedPropagator1, par().q1) and
envHasType(SlicedPropagator2, par().q2))
{
SlicedPropagator1 &q1 = envGet(SlicedPropagator1, par().q1);
SlicedPropagator2 &q2 = envGet(SlicedPropagator2, par().q2);
auto &q1 = envGet(SlicedPropagator1, par().q1);
auto &q2 = envGet(SlicedPropagator2, par().q2);
LOG(Message) << "(propagator already sinked)" << std::endl;
for (unsigned int i = 0; i < result.size(); ++i)
@ -209,10 +210,10 @@ void TMeson<FImpl1, FImpl2>::execute(void)
}
else
{
PropagatorField1 &q1 = envGet(PropagatorField1, par().q1);
PropagatorField2 &q2 = envGet(PropagatorField2, par().q2);
LatticeComplex c(env().getGrid());
auto &q1 = envGet(PropagatorField1, par().q1);
auto &q2 = envGet(PropagatorField2, par().q2);
envGetTmp(LatticeComplex, c);
LOG(Message) << "(using sink '" << par().sink << "')" << std::endl;
for (unsigned int i = 0; i < result.size(); ++i)
{

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

@ -73,6 +73,7 @@ public:
virtual ~TWardIdentity(void) = default;
// dependency relation
virtual std::vector<std::string> getInput(void);
virtual std::vector<std::string> getReference(void);
virtual std::vector<std::string> getOutput(void);
protected:
// setup
@ -103,10 +104,18 @@ std::vector<std::string> TWardIdentity<FImpl>::getInput(void)
return in;
}
template <typename FImpl>
std::vector<std::string> TWardIdentity<FImpl>::getReference(void)
{
std::vector<std::string> ref = {};
return ref;
}
template <typename FImpl>
std::vector<std::string> TWardIdentity<FImpl>::getOutput(void)
{
std::vector<std::string> out = {getName()};
std::vector<std::string> out = {};
return out;
}
@ -120,6 +129,15 @@ void TWardIdentity<FImpl>::setup(void)
{
HADRON_ERROR(Size, "Ls mismatch between quark action and propagator");
}
envTmpLat(PropagatorField, "tmp");
envTmpLat(PropagatorField, "vector_WI");
if (par().test_axial)
{
envTmpLat(PropagatorField, "psi");
envTmpLat(LatticeComplex, "PP");
envTmpLat(LatticeComplex, "axial_defect");
envTmpLat(LatticeComplex, "PJ5q");
}
}
// execution ///////////////////////////////////////////////////////////////////
@ -129,12 +147,13 @@ 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));
auto &q = envGet(PropagatorField, par().q);
auto &act = envGet(FMat, par().action);
Gamma g5(Gamma::Algebra::Gamma5);
// Compute D_mu V_mu, D here is backward derivative.
envGetTmp(PropagatorField, tmp);
envGetTmp(PropagatorField, vector_WI);
vector_WI = zero;
for (unsigned int mu = 0; mu < Nd; ++mu)
{
@ -149,9 +168,10 @@ void TWardIdentity<FImpl>::execute(void)
if (par().test_axial)
{
PropagatorField psi(env().getGrid());
LatticeComplex PP(env().getGrid()), axial_defect(env().getGrid()),
PJ5q(env().getGrid());
envGetTmp(PropagatorField, psi);
envGetTmp(LatticeComplex, PP);
envGetTmp(LatticeComplex, axial_defect);
envGetTmp(LatticeComplex, PJ5q);
std::vector<TComplex> axial_buf;
// Compute <P|D_mu A_mu>, D is backwards derivative.

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

@ -99,6 +99,7 @@ public:\
virtual ~T##modname(void) = default;\
/* dependency relation */ \
virtual std::vector<std::string> getInput(void);\
virtual std::vector<std::string> getReference(void);\
virtual std::vector<std::string> getOutput(void);\
public:\
std::vector<std::string> VA_label = {"V", "A"};\

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

@ -74,9 +74,16 @@ std::vector<std::string> TWeakHamiltonianEye::getInput(void)
return in;
}
std::vector<std::string> TWeakHamiltonianEye::getReference(void)
{
std::vector<std::string> out = {};
return out;
}
std::vector<std::string> TWeakHamiltonianEye::getOutput(void)
{
std::vector<std::string> out = {getName()};
std::vector<std::string> out = {};
return out;
}
@ -84,7 +91,15 @@ std::vector<std::string> TWeakHamiltonianEye::getOutput(void)
// setup ///////////////////////////////////////////////////////////////////////
void TWeakHamiltonianEye::setup(void)
{
unsigned int ndim = env().getNd();
envTmpLat(LatticeComplex, "expbuf");
envTmpLat(PropagatorField, "tmp1");
envTmpLat(LatticeComplex, "tmp2");
envTmp(std::vector<PropagatorField>, "S_body", 1, ndim, PropagatorField(env().getGrid()));
envTmp(std::vector<PropagatorField>, "S_loop", 1, ndim, PropagatorField(env().getGrid()));
envTmp(std::vector<LatticeComplex>, "E_body", 1, ndim, LatticeComplex(env().getGrid()));
envTmp(std::vector<LatticeComplex>, "E_loop", 1, ndim, LatticeComplex(env().getGrid()));
}
// execution ///////////////////////////////////////////////////////////////////
@ -96,22 +111,22 @@ void TWeakHamiltonianEye::execute(void)
<< "'." << std::endl;
CorrWriter writer(par().output);
SlicedPropagator &q1 = *env().template getObject<SlicedPropagator>(par().q1);
PropagatorField &q2 = *env().template getObject<PropagatorField>(par().q2);
PropagatorField &q3 = *env().template getObject<PropagatorField>(par().q3);
PropagatorField &q4 = *env().template getObject<PropagatorField>(par().q4);
auto &q1 = envGet(SlicedPropagator, par().q1);
auto &q2 = envGet(PropagatorField, par().q2);
auto &q3 = envGet(PropagatorField, par().q3);
auto &q4 = envGet(PropagatorField, par().q4);
Gamma g5 = Gamma(Gamma::Algebra::Gamma5);
LatticeComplex expbuf(env().getGrid());
std::vector<TComplex> corrbuf;
std::vector<Result> result(n_eye_diag);
unsigned int ndim = env().getNd();
PropagatorField tmp1(env().getGrid());
LatticeComplex tmp2(env().getGrid());
std::vector<PropagatorField> S_body(ndim, tmp1);
std::vector<PropagatorField> S_loop(ndim, tmp1);
std::vector<LatticeComplex> E_body(ndim, tmp2);
std::vector<LatticeComplex> E_loop(ndim, tmp2);
envGetTmp(LatticeComplex, expbuf);
envGetTmp(PropagatorField, tmp1);
envGetTmp(LatticeComplex, tmp2);
envGetTmp(std::vector<PropagatorField>, S_body);
envGetTmp(std::vector<PropagatorField>, S_loop);
envGetTmp(std::vector<LatticeComplex>, E_body);
envGetTmp(std::vector<LatticeComplex>, E_loop);
// Get sink timeslice of q1.
SitePropagator q1Snk = q1[par().tSnk];

38
extras/Hadrons/Modules/MContraction/WeakHamiltonianNonEye.cc
View File

@ -74,9 +74,15 @@ std::vector<std::string> TWeakHamiltonianNonEye::getInput(void)
return in;
}
std::vector<std::string> TWeakHamiltonianNonEye::getReference(void)
{
std::vector<std::string> out = {};
return out;
}
std::vector<std::string> TWeakHamiltonianNonEye::getOutput(void)
{
std::vector<std::string> out = {getName()};
std::vector<std::string> out = {};
return out;
}
@ -84,7 +90,15 @@ std::vector<std::string> TWeakHamiltonianNonEye::getOutput(void)
// setup ///////////////////////////////////////////////////////////////////////
void TWeakHamiltonianNonEye::setup(void)
{
unsigned int ndim = env().getNd();
envTmpLat(LatticeComplex, "expbuf");
envTmpLat(PropagatorField, "tmp1");
envTmpLat(LatticeComplex, "tmp2");
envTmp(std::vector<PropagatorField>, "C_i_side_loop", 1, ndim, PropagatorField(env().getGrid()));
envTmp(std::vector<PropagatorField>, "C_f_side_loop", 1, ndim, PropagatorField(env().getGrid()));
envTmp(std::vector<LatticeComplex>, "W_i_side_loop", 1, ndim, LatticeComplex(env().getGrid()));
envTmp(std::vector<LatticeComplex>, "W_f_side_loop", 1, ndim, LatticeComplex(env().getGrid()));
}
// execution ///////////////////////////////////////////////////////////////////
@ -96,22 +110,22 @@ void TWeakHamiltonianNonEye::execute(void)
<< "'." << std::endl;
CorrWriter writer(par().output);
PropagatorField &q1 = *env().template getObject<PropagatorField>(par().q1);
PropagatorField &q2 = *env().template getObject<PropagatorField>(par().q2);
PropagatorField &q3 = *env().template getObject<PropagatorField>(par().q3);
PropagatorField &q4 = *env().template getObject<PropagatorField>(par().q4);
auto &q1 = envGet(PropagatorField, par().q1);
auto &q2 = envGet(PropagatorField, par().q2);
auto &q3 = envGet(PropagatorField, par().q3);
auto &q4 = envGet(PropagatorField, par().q4);
Gamma g5 = Gamma(Gamma::Algebra::Gamma5);
LatticeComplex expbuf(env().getGrid());
std::vector<TComplex> corrbuf;
std::vector<Result> result(n_noneye_diag);
unsigned int ndim = env().getNd();
PropagatorField tmp1(env().getGrid());
LatticeComplex tmp2(env().getGrid());
std::vector<PropagatorField> C_i_side_loop(ndim, tmp1);
std::vector<PropagatorField> C_f_side_loop(ndim, tmp1);
std::vector<LatticeComplex> W_i_side_loop(ndim, tmp2);
std::vector<LatticeComplex> W_f_side_loop(ndim, tmp2);
envGetTmp(LatticeComplex, expbuf);
envGetTmp(PropagatorField, tmp1);
envGetTmp(LatticeComplex, tmp2);
envGetTmp(std::vector<PropagatorField>, C_i_side_loop);
envGetTmp(std::vector<PropagatorField>, C_f_side_loop);
envGetTmp(std::vector<LatticeComplex>, W_i_side_loop);
envGetTmp(std::vector<LatticeComplex>, W_f_side_loop);
// Setup for C-type contractions.
for (int mu = 0; mu < ndim; ++mu)

33
extras/Hadrons/Modules/MContraction/WeakNeutral4ptDisc.cc
View File

@ -76,9 +76,16 @@ std::vector<std::string> TWeakNeutral4ptDisc::getInput(void)
return in;
}
std::vector<std::string> TWeakNeutral4ptDisc::getReference(void)
{
std::vector<std::string> ref = {};
return ref;
}
std::vector<std::string> TWeakNeutral4ptDisc::getOutput(void)
{
std::vector<std::string> out = {getName()};
std::vector<std::string> out = {};
return out;
}
@ -86,7 +93,13 @@ std::vector<std::string> TWeakNeutral4ptDisc::getOutput(void)
// setup ///////////////////////////////////////////////////////////////////////
void TWeakNeutral4ptDisc::setup(void)
{
unsigned int ndim = env().getNd();
envTmpLat(LatticeComplex, "expbuf");
envTmpLat(PropagatorField, "tmp");
envTmpLat(LatticeComplex, "curr");
envTmp(std::vector<PropagatorField>, "meson", 1, ndim, PropagatorField(env().getGrid()));
envTmp(std::vector<PropagatorField>, "loop", 1, ndim, PropagatorField(env().getGrid()));
}
// execution ///////////////////////////////////////////////////////////////////
@ -98,20 +111,20 @@ void TWeakNeutral4ptDisc::execute(void)
<< "'." << std::endl;
CorrWriter writer(par().output);
PropagatorField &q1 = *env().template getObject<PropagatorField>(par().q1);
PropagatorField &q2 = *env().template getObject<PropagatorField>(par().q2);
PropagatorField &q3 = *env().template getObject<PropagatorField>(par().q3);
PropagatorField &q4 = *env().template getObject<PropagatorField>(par().q4);
auto &q1 = envGet(PropagatorField, par().q1);
auto &q2 = envGet(PropagatorField, par().q2);
auto &q3 = envGet(PropagatorField, par().q3);
auto &q4 = envGet(PropagatorField, par().q4);
Gamma g5 = Gamma(Gamma::Algebra::Gamma5);
LatticeComplex expbuf(env().getGrid());
std::vector<TComplex> corrbuf;
std::vector<Result> result(n_neut_disc_diag);
unsigned int ndim = env().getNd();
PropagatorField tmp(env().getGrid());
std::vector<PropagatorField> meson(ndim, tmp);
std::vector<PropagatorField> loop(ndim, tmp);
LatticeComplex curr(env().getGrid());
envGetTmp(LatticeComplex, expbuf);
envGetTmp(PropagatorField, tmp);
envGetTmp(LatticeComplex, curr);
envGetTmp(std::vector<PropagatorField>, meson);
envGetTmp(std::vector<PropagatorField>, loop);
// Setup for type 1 contractions.
for (int mu = 0; mu < ndim; ++mu)

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

@ -154,14 +154,14 @@ void TGaugeProp<FImpl>::execute(void)
LOG(Message) << "Computing quark propagator '" << getName() << "'"
<< std::endl;
FermionField &source = envGetTmp(FermionField, "source");
FermionField &sol = envGetTmp(FermionField, "sol");
FermionField &tmp = envGetTmp(FermionField, "tmp");
std::string propName = (Ls_ == 1) ? getName() : (getName() + "_5d");
PropagatorField &prop = envGet(PropagatorField, propName);
PropagatorField &fullSrc = envGet(PropagatorField, par().source);
SolverFn &solver = envGet(SolverFn, par().solver);
auto &prop = envGet(PropagatorField, propName);
auto &fullSrc = envGet(PropagatorField, par().source);
auto &solver = envGet(SolverFn, par().solver);
envGetTmp(FermionField, source);
envGetTmp(FermionField, sol);
envGetTmp(FermionField, tmp);
LOG(Message) << "Inverting using solver '" << par().solver
<< "' on source '" << par().source << "'" << std::endl;
for (unsigned int s = 0; s < Ns; ++s)

15
extras/Hadrons/Modules/MGauge/Load.cc
View File

@ -49,6 +49,13 @@ std::vector<std::string> TLoad::getInput(void)
return in;
}
std::vector<std::string> TLoad::getReference(void)
{
std::vector<std::string> ref;
return ref;
}
std::vector<std::string> TLoad::getOutput(void)
{
std::vector<std::string> out = {getName()};
@ -59,7 +66,7 @@ std::vector<std::string> TLoad::getOutput(void)
// setup ///////////////////////////////////////////////////////////////////////
void TLoad::setup(void)
{
env().registerLattice<LatticeGaugeField>(getName());
envCreateLat(LatticeGaugeField, getName());
}
// execution ///////////////////////////////////////////////////////////////////
@ -67,11 +74,11 @@ void TLoad::execute(void)
{
FieldMetaData header;
std::string fileName = par().file + "."
+ std::to_string(env().getTrajectory());
+ std::to_string(vm().getTrajectory());
LOG(Message) << "Loading NERSC configuration from file '" << fileName
<< "'" << std::endl;
LatticeGaugeField &U = *env().createLattice<LatticeGaugeField>(getName());
auto &U = envGet(LatticeGaugeField, getName());
NerscIO::readConfiguration(U, header, fileName);
LOG(Message) << "NERSC header:" << std::endl;
dump_meta_data(header, LOG(Message));

1
extras/Hadrons/Modules/MGauge/Load.hpp
View File

@ -57,6 +57,7 @@ public:
virtual ~TLoad(void) = default;
// dependency relation
virtual std::vector<std::string> getInput(void);
virtual std::vector<std::string> getReference(void);
virtual std::vector<std::string> getOutput(void);
protected:
// setup

16
extras/Hadrons/Modules/MGauge/Random.cc
View File

@ -44,7 +44,16 @@ TRandom::TRandom(const std::string name)
// dependencies/products ///////////////////////////////////////////////////////
std::vector<std::string> TRandom::getInput(void)
{
return std::vector<std::string>();
std::vector<std::string> in;
return in;
}
std::vector<std::string> TRandom::getReference(void)
{
std::vector<std::string> ref;
return ref;
}
std::vector<std::string> TRandom::getOutput(void)
@ -57,13 +66,14 @@ std::vector<std::string> TRandom::getOutput(void)
// setup ///////////////////////////////////////////////////////////////////////
void TRandom::setup(void)
{
env().registerLattice<LatticeGaugeField>(getName());
envCreateLat(LatticeGaugeField, getName());
}
// execution ///////////////////////////////////////////////////////////////////
void TRandom::execute(void)
{
LOG(Message) << "Generating random gauge configuration" << std::endl;
LatticeGaugeField &U = *env().createLattice<LatticeGaugeField>(getName());
auto &U = envGet(LatticeGaugeField, getName());
SU3::HotConfiguration(*env().get4dRng(), U);
}

1
extras/Hadrons/Modules/MGauge/Random.hpp
View File

@ -50,6 +50,7 @@ public:
virtual ~TRandom(void) = default;
// dependency relation
virtual std::vector<std::string> getInput(void);
virtual std::vector<std::string> getReference(void);
virtual std::vector<std::string> getOutput(void);
protected:
// setup

29
extras/Hadrons/Modules/MGauge/StochEm.cc
View File

@ -47,6 +47,13 @@ std::vector<std::string> TStochEm::getInput(void)
return in;
}
std::vector<std::string> TStochEm::getReference(void)
{
std::vector<std::string> ref = {};
return ref;
}
std::vector<std::string> TStochEm::getOutput(void)
{
std::vector<std::string> out = {getName()};
@ -57,32 +64,28 @@ std::vector<std::string> TStochEm::getOutput(void)
// setup ///////////////////////////////////////////////////////////////////////
void TStochEm::setup(void)
{
if (!env().hasRegisteredObject("_" + getName() + "_weight"))
if (!env().hasCreatedObject("_" + getName() + "_weight"))
{
env().registerLattice<EmComp>("_" + getName() + "_weight");
envCacheLat(EmComp, "_" + getName() + "_weight");
}
env().registerLattice<EmField>(getName());
envCreateLat(EmField, getName());
}
// execution ///////////////////////////////////////////////////////////////////
void TStochEm::execute(void)
{
LOG(Message) << "Generating stochatic EM potential..." << std::endl;
PhotonR photon(par().gauge, par().zmScheme);
EmField &a = *env().createLattice<EmField>(getName());
EmComp *w;
auto &a = envGet(EmField, getName());
auto &w = envGet(EmComp, "_" + getName() + "_weight");
if (!env().hasCreatedObject("_" + getName() + "_weight"))
{
LOG(Message) << "Caching stochatic EM potential weight (gauge: "
<< par().gauge << ", zero-mode scheme: "
<< par().zmScheme << ")..." << std::endl;
w = env().createLattice<EmComp>("_" + getName() + "_weight");
photon.StochasticWeight(*w);
photon.StochasticWeight(w);
}
else
{
w = env().getObject<EmComp>("_" + getName() + "_weight");
}
LOG(Message) << "Generating stochatic EM potential..." << std::endl;
photon.StochasticField(a, *env().get4dRng(), *w);
photon.StochasticField(a, *env().get4dRng(), w);
}

1
extras/Hadrons/Modules/MGauge/StochEm.hpp
View File

@ -59,6 +59,7 @@ public:
virtual ~TStochEm(void) = default;
// dependency relation
virtual std::vector<std::string> getInput(void);
virtual std::vector<std::string> getReference(void);
virtual std::vector<std::string> getOutput(void);
protected:
// setup

1
extras/Hadrons/Modules/MGauge/Unit.cc
View File

@ -71,6 +71,7 @@ void TUnit::setup(void)
void TUnit::execute(void)
{
LOG(Message) << "Creating unit gauge configuration" << std::endl;
auto &U = envGet(LatticeGaugeField, getName());
SU3::ColdConfiguration(*env().get4dRng(), U);
}

18
extras/Hadrons/Modules/MLoop/NoiseLoop.hpp
View File

@ -73,6 +73,7 @@ public:
virtual ~TNoiseLoop(void) = default;
// dependency relation
virtual std::vector<std::string> getInput(void);
virtual std::vector<std::string> getReference(void);
virtual std::vector<std::string> getOutput(void);
protected:
// setup
@ -101,6 +102,15 @@ std::vector<std::string> TNoiseLoop<FImpl>::getInput(void)
return in;
}
template <typename FImpl>
std::vector<std::string> TNoiseLoop<FImpl>::getReference(void)
{
std::vector<std::string> out = {};
return out;
}
template <typename FImpl>
std::vector<std::string> TNoiseLoop<FImpl>::getOutput(void)
{
@ -113,16 +123,16 @@ std::vector<std::string> TNoiseLoop<FImpl>::getOutput(void)
template <typename FImpl>
void TNoiseLoop<FImpl>::setup(void)
{
env().template registerLattice<PropagatorField>(getName());
envCreateLat(PropagatorField, getName());
}
// execution ///////////////////////////////////////////////////////////////////
template <typename FImpl>
void TNoiseLoop<FImpl>::execute(void)
{
PropagatorField &loop = *env().template createLattice<PropagatorField>(getName());
PropagatorField &q = *env().template getObject<PropagatorField>(par().q);
PropagatorField &eta = *env().template getObject<PropagatorField>(par().eta);
auto &loop = envGet(PropagatorField, getName());
auto &q = envGet(PropagatorField, par().q);
auto &eta = envGet(PropagatorField, par().eta);
loop = q*adj(eta);
}

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

@ -123,17 +123,18 @@ void TPoint<FImpl>::setup(void)
template <typename FImpl>
void TPoint<FImpl>::execute(void)
{
std::vector<Real> p = strToVec<Real>(par().mom);
LatticeComplex &ph = envGet(LatticeComplex, momphName_);
Complex i(0.0,1.0);
LOG(Message) << "Setting up point sink function for momentum ["
<< par().mom << "]" << std::endl;
auto &ph = envGet(LatticeComplex, momphName_);
if (!hasPhase_)
{
LatticeComplex &coor = envGetTmp(LatticeComplex, "coor");
Complex i(0.0,1.0);
std::vector<Real> p;
envGetTmp(LatticeComplex, coor);
p = strToVec<Real>(par().mom);
ph = zero;
for(unsigned int mu = 0; mu < env().getNd(); mu++)
{

23
extras/Hadrons/Modules/MSink/Smear.hpp
View File

@ -61,6 +61,7 @@ public:
virtual ~TSmear(void) = default;
// dependency relation
virtual std::vector<std::string> getInput(void);
virtual std::vector<std::string> getReference(void);
virtual std::vector<std::string> getOutput(void);
protected:
// setup
@ -89,6 +90,14 @@ std::vector<std::string> TSmear<FImpl>::getInput(void)
return in;
}
template <typename FImpl>
std::vector<std::string> TSmear<FImpl>::getReference(void)
{
std::vector<std::string> ref = {};
return ref;
}
template <typename FImpl>
std::vector<std::string> TSmear<FImpl>::getOutput(void)
{
@ -101,9 +110,7 @@ std::vector<std::string> TSmear<FImpl>::getOutput(void)
template <typename FImpl>
void TSmear<FImpl>::setup(void)
{
unsigned int nt = env().getDim(Tp);
unsigned int size = nt * sizeof(SitePropagator);
env().registerObject(getName(), size);
envCreate(SlicedPropagator, getName(), 1, env().getDim(Tp));
}
// execution ///////////////////////////////////////////////////////////////////
@ -114,11 +121,11 @@ void TSmear<FImpl>::execute(void)
<< "' 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);
auto &sink = envGet(SinkFn, par().sink);
auto &q = envGet(PropagatorField, par().q);
auto &out = envGet(SlicedPropagator, getName());
out = sink(q);
}
END_MODULE_NAMESPACE

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

@ -128,12 +128,13 @@ void TPoint<FImpl>::setup(void)
template <typename FImpl>
void TPoint<FImpl>::execute(void)
{
std::vector<int> position = strToVec<int>(par().position);
SitePropagator id;
LOG(Message) << "Creating point source at position [" << par().position
<< "]" << std::endl;
PropagatorField &src = envGet(PropagatorField, getName());
std::vector<int> position = strToVec<int>(par().position);
auto &src = envGet(PropagatorField, getName());
SitePropagator id;
id = 1.;
src = zero;
pokeSite(id, src, position);

17
extras/Hadrons/Modules/MSource/SeqConserved.hpp
View File

@ -82,6 +82,7 @@ public:
virtual ~TSeqConserved(void) = default;
// dependency relation
virtual std::vector<std::string> getInput(void);
virtual std::vector<std::string> getReference(void);
virtual std::vector<std::string> getOutput(void);
protected:
// setup
@ -110,6 +111,14 @@ std::vector<std::string> TSeqConserved<FImpl>::getInput(void)
return in;
}
template <typename FImpl>
std::vector<std::string> TSeqConserved<FImpl>::getReference(void)
{
std::vector<std::string> ref = {};
return ref;
}
template <typename FImpl>
std::vector<std::string> TSeqConserved<FImpl>::getOutput(void)
{
@ -123,7 +132,7 @@ template <typename FImpl>
void TSeqConserved<FImpl>::setup(void)
{
auto Ls_ = env().getObjectLs(par().action);
env().template registerLattice<PropagatorField>(getName(), Ls_);
envCreateLat(PropagatorField, getName(), Ls_);
}
// execution ///////////////////////////////////////////////////////////////////
@ -143,9 +152,9 @@ void TSeqConserved<FImpl>::execute(void)
<< 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));
auto &src = envGet(PropagatorField, getName());
auto &q = envGet(PropagatorField, par().q);
auto &mat = envGet(FMat, par().action);
std::vector<Real> mom = strToVec<Real>(par().mom);
mat.SeqConservedCurrent(q, src, par().curr_type, par().mu,

39
extras/Hadrons/Modules/MSource/SeqGamma.hpp
View File

@ -80,12 +80,16 @@ public:
virtual ~TSeqGamma(void) = default;
// dependency relation
virtual std::vector<std::string> getInput(void);
virtual std::vector<std::string> getReference(void);
virtual std::vector<std::string> getOutput(void);
protected:
// setup
virtual void setup(void);
// execution
virtual void execute(void);
private:
bool hasPhase_{false};
std::string momphName_, tName_;
};
MODULE_REGISTER_NS(SeqGamma, TSeqGamma<FIMPL>, MSource);
@ -97,6 +101,8 @@ MODULE_REGISTER_NS(SeqGamma, TSeqGamma<FIMPL>, MSource);
template <typename FImpl>
TSeqGamma<FImpl>::TSeqGamma(const std::string name)
: Module<SeqGammaPar>(name)
, momphName_ (name + "_momph")
, tName_ (name + "_t")
{}
// dependencies/products ///////////////////////////////////////////////////////
@ -108,6 +114,14 @@ std::vector<std::string> TSeqGamma<FImpl>::getInput(void)
return in;
}
template <typename FImpl>
std::vector<std::string> TSeqGamma<FImpl>::getReference(void)
{
std::vector<std::string> ref = {};
return ref;
}
template <typename FImpl>
std::vector<std::string> TSeqGamma<FImpl>::getOutput(void)
{
@ -120,7 +134,10 @@ std::vector<std::string> TSeqGamma<FImpl>::getOutput(void)
template <typename FImpl>
void TSeqGamma<FImpl>::setup(void)
{
env().template registerLattice<PropagatorField>(getName());
envCreateLat(PropagatorField, getName());
envCacheLat(Lattice<iScalar<vInteger>>, tName_);
envCacheLat(LatticeComplex, momphName_);
envTmpLat(LatticeComplex, "coor");
}
// execution ///////////////////////////////////////////////////////////////////
@ -138,23 +155,29 @@ void TSeqGamma<FImpl>::execute(void)
<< " sequential source for "
<< par().tA << " <= t <= " << par().tB << std::endl;
}
PropagatorField &src = *env().template createLattice<PropagatorField>(getName());
PropagatorField &q = *env().template getObject<PropagatorField>(par().q);
Lattice<iScalar<vInteger>> t(env().getGrid());
LatticeComplex ph(env().getGrid()), coor(env().getGrid());
auto &src = envGet(PropagatorField, getName());
auto &q = envGet(PropagatorField, par().q);
auto &ph = envGet(LatticeComplex, momphName_);
auto &t = envGet(Lattice<iScalar<vInteger>>, tName_);
Gamma g(par().gamma);
std::vector<Real> p;
Complex i(0.0,1.0);
if (!hasPhase_)
{
Complex i(0.0,1.0);
std::vector<Real> p;
envGetTmp(LatticeComplex, coor);
p = strToVec<Real>(par().mom);
ph = zero;
for(unsigned int mu = 0; mu < env().getNd(); mu++)
{
LatticeCoordinate(coor, mu);
ph = ph + p[mu]*coor*((1./(env().getGrid()->_fdimensions[mu])));
ph = ph + (p[mu]/env().getGrid()->_fdimensions[mu])*coor;
}
ph = exp((Real)(2*M_PI)*i*ph);
LatticeCoordinate(t, Tp);
hasPhase_ = true;
}
src = where((t >= par().tA) and (t <= par().tB), ph*(g*q), 0.*q);
}

39
extras/Hadrons/Modules/MSource/Wall.hpp
View File

@ -72,12 +72,16 @@ public:
virtual ~TWall(void) = default;
// dependency relation
virtual std::vector<std::string> getInput(void);
virtual std::vector<std::string> getReference(void);
virtual std::vector<std::string> getOutput(void);
protected:
// setup
virtual void setup(void);
// execution
virtual void execute(void);
private:
bool hasPhase_{false};
std::string momphName_, tName_;
};
MODULE_REGISTER_NS(Wall, TWall<FIMPL>, MSource);
@ -89,17 +93,27 @@ MODULE_REGISTER_NS(Wall, TWall<FIMPL>, MSource);
template <typename FImpl>
TWall<FImpl>::TWall(const std::string name)
: Module<WallPar>(name)
, momphName_ (name + "_momph")
, tName_ (name + "_t")
{}
// dependencies/products ///////////////////////////////////////////////////////
template <typename FImpl>
std::vector<std::string> TWall<FImpl>::getInput(void)
{
std::vector<std::string> in;
std::vector<std::string> in = {};
return in;
}
template <typename FImpl>
std::vector<std::string> TWall<FImpl>::getReference(void)
{
std::vector<std::string> ref = {};
return ref;
}
template <typename FImpl>
std::vector<std::string> TWall<FImpl>::getOutput(void)
{
@ -112,7 +126,7 @@ std::vector<std::string> TWall<FImpl>::getOutput(void)
template <typename FImpl>
void TWall<FImpl>::setup(void)
{
env().template registerLattice<PropagatorField>(getName());
envCreateLat(PropagatorField, getName());
}
// execution ///////////////////////////////////////////////////////////////////
@ -122,21 +136,28 @@ void TWall<FImpl>::execute(void)
LOG(Message) << "Generating wall source at t = " << par().tW
<< " with momentum " << par().mom << std::endl;
PropagatorField &src = *env().template createLattice<PropagatorField>(getName());
Lattice<iScalar<vInteger>> t(env().getGrid());
LatticeComplex ph(env().getGrid()), coor(env().getGrid());
std::vector<Real> p;
Complex i(0.0,1.0);
auto &src = envGet(PropagatorField, getName());
auto &ph = envGet(LatticeComplex, momphName_);
auto &t = envGet(Lattice<iScalar<vInteger>>, tName_);
if (!hasPhase_)
{
Complex i(0.0,1.0);
std::vector<Real> p;
envGetTmp(LatticeComplex, coor);
p = strToVec<Real>(par().mom);
ph = zero;
for(unsigned int mu = 0; mu < Nd; mu++)
for(unsigned int mu = 0; mu < env().getNd(); mu++)
{
LatticeCoordinate(coor, mu);
ph = ph + p[mu]*coor*((1./(env().getGrid()->_fdimensions[mu])));
ph = ph + (p[mu]/env().getGrid()->_fdimensions[mu])*coor;
}
ph = exp((Real)(2*M_PI)*i*ph);
LatticeCoordinate(t, Tp);
hasPhase_ = true;
}
src = 1.;
src = where((t == par().tW), src*ph, 0.*src);
}

32
extras/Hadrons/Modules/MSource/Z2.hpp
View File

@ -75,12 +75,16 @@ public:
virtual ~TZ2(void) = default;
// dependency relation
virtual std::vector<std::string> getInput(void);
virtual std::vector<std::string> getReference(void);
virtual std::vector<std::string> getOutput(void);
protected:
// setup
virtual void setup(void);
// execution
virtual void execute(void);
private:
bool hasT_{false};
std::string tName_;
};
MODULE_REGISTER_NS(Z2, TZ2<FIMPL>, MSource);
@ -93,6 +97,7 @@ MODULE_REGISTER_NS(ScalarZ2, TZ2<ScalarImplCR>, MSource);
template <typename FImpl>
TZ2<FImpl>::TZ2(const std::string name)
: Module<Z2Par>(name)
, tName_ (name + "_t")
{}
// dependencies/products ///////////////////////////////////////////////////////
@ -104,6 +109,14 @@ std::vector<std::string> TZ2<FImpl>::getInput(void)
return in;
}
template <typename FImpl>
std::vector<std::string> TZ2<FImpl>::getReference(void)
{
std::vector<std::string> ref = {};
return ref;
}
template <typename FImpl>
std::vector<std::string> TZ2<FImpl>::getOutput(void)
{
@ -116,17 +129,15 @@ std::vector<std::string> TZ2<FImpl>::getOutput(void)
template <typename FImpl>
void TZ2<FImpl>::setup(void)
{
env().template registerLattice<PropagatorField>(getName());
envCreateLat(PropagatorField, getName());
envCacheLat(Lattice<iScalar<vInteger>>, tName_);
envTmpLat(LatticeComplex, "eta");
}
// execution ///////////////////////////////////////////////////////////////////
template <typename FImpl>
void TZ2<FImpl>::execute(void)
{
Lattice<iScalar<vInteger>> t(env().getGrid());
LatticeComplex eta(env().getGrid());
Complex shift(1., 1.);
if (par().tA == par().tB)
{
LOG(Message) << "Generating Z_2 wall source at t= " << par().tA
@ -137,8 +148,17 @@ void TZ2<FImpl>::execute(void)
LOG(Message) << "Generating Z_2 band for " << par().tA << " <= t <= "
<< par().tB << std::endl;
}
PropagatorField &src = *env().template createLattice<PropagatorField>(getName());
auto &src = envGet(PropagatorField, getName());
auto &t = envGet(Lattice<iScalar<vInteger>>, getName());
Complex shift(1., 1.);
if (!hasT_)
{
LatticeCoordinate(t, Tp);
hasT_ = true;
}
envGetTmp(LatticeComplex, eta);
bernoulli(*env().get4dRng(), eta);
eta = (2.*eta - shift)*(1./::sqrt(2.));
eta = where((t >= par().tA) and (t <= par().tB), eta, 0.*eta);

30
extras/Hadrons/Modules/MUtilities/TestSeqConserved.hpp
View File

@ -79,6 +79,7 @@ public:
virtual ~TTestSeqConserved(void) = default;
// dependency relation
virtual std::vector<std::string> getInput(void);
virtual std::vector<std::string> getReference(void);
virtual std::vector<std::string> getOutput(void);
protected:
// setup
@ -107,6 +108,14 @@ std::vector<std::string> TTestSeqConserved<FImpl>::getInput(void)
return in;
}
template <typename FImpl>
std::vector<std::string> TTestSeqConserved<FImpl>::getReference(void)
{
std::vector<std::string> ref = {};
return ref;
}
template <typename FImpl>
std::vector<std::string> TTestSeqConserved<FImpl>::getOutput(void)
{
@ -124,16 +133,20 @@ void TTestSeqConserved<FImpl>::setup(void)
{
HADRON_ERROR(Size, "Ls mismatch between quark action and propagator");
}
envTmpLat(PropagatorField, "tmp");
envTmpLat(LatticeComplex, "c");
}
// 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));
// Check sequential insertion of current gives same result as conserved
// current sink upon contraction. Assume q uses a point source.
auto &q = envGet(PropagatorField, par().q);
auto &qSeq = envGet(PropagatorField, par().qSeq);
auto &act = envGet(FMat, par().action);
Gamma g5(Gamma::Algebra::Gamma5);
Gamma::Algebra gA = (par().curr == Current::Axial) ?
Gamma::Algebra::Gamma5 :
@ -142,18 +155,15 @@ void TTestSeqConserved<FImpl>::execute(void)
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;
envGetTmp(PropagatorField, tmp);
envGetTmp(LatticeComplex, c);
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]);

20
extras/Hadrons/Modules/MUtilities/TestSeqGamma.hpp
View File

@ -63,6 +63,7 @@ public:
virtual ~TTestSeqGamma(void) = default;
// dependency relation
virtual std::vector<std::string> getInput(void);
virtual std::vector<std::string> getReference(void);
virtual std::vector<std::string> getOutput(void);
protected:
// setup
@ -91,6 +92,14 @@ std::vector<std::string> TTestSeqGamma<FImpl>::getInput(void)
return in;
}
template <typename FImpl>
std::vector<std::string> TTestSeqGamma<FImpl>::getReference(void)
{
std::vector<std::string> ref = {};
return ref;
}
template <typename FImpl>
std::vector<std::string> TTestSeqGamma<FImpl>::getOutput(void)
{
@ -103,26 +112,27 @@ std::vector<std::string> TTestSeqGamma<FImpl>::getOutput(void)
template <typename FImpl>
void TTestSeqGamma<FImpl>::setup(void)
{
envTmpLat(LatticeComplex, "c");
}
// 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());
auto &q = envGet(PropagatorField, par().q);
auto &qSeq = envGet(PropagatorField, par().qSeq);
Gamma g5(Gamma::Algebra::Gamma5);
Gamma g(par().gamma);
SitePropagator qSite;
Complex test, check;
std::vector<TComplex> check_buf;
std::vector<int> siteCoord;
// 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;
envGetTmp(LatticeComplex, c);
siteCoord = strToVec<int>(par().origin);
peekSite(qSite, qSeq, siteCoord);
test = trace(g*qSite);

63
extras/Hadrons/modules.inc
View File

@ -1,45 +1,38 @@
modules_cc =\
Modules/MGauge/Unit.cc
# Modules/MContraction/WeakHamiltonianEye.cc \
# Modules/MContraction/WeakHamiltonianNonEye.cc \
# Modules/MContraction/WeakNeutral4ptDisc.cc \
# Modules/MGauge/Load.cc \
# Modules/MGauge/Random.cc \
# Modules/MGauge/StochEm.cc \
# Modules/MScalar/ChargedProp.cc \
# Modules/MScalar/FreeProp.cc
Modules/MGauge/Unit.cc \
Modules/MContraction/WeakHamiltonianEye.cc \
Modules/MContraction/WeakHamiltonianNonEye.cc \
Modules/MContraction/WeakNeutral4ptDisc.cc \
Modules/MGauge/Load.cc \
Modules/MGauge/Random.cc \
Modules/MGauge/StochEm.cc
modules_hpp =\
Modules/MAction/DWF.hpp \
Modules/MAction/Wilson.hpp \
Modules/MSink/Point.hpp \
Modules/MSource/Point.hpp \
Modules/MGauge/Load.hpp \
Modules/MGauge/Random.hpp \
Modules/MGauge/StochEm.hpp \
Modules/MGauge/Unit.hpp \
Modules/MSolver/RBPrecCG.hpp \
Modules/MFermion/GaugeProp.hpp \
Modules/MContraction/Meson.hpp
# Modules/MContraction/Baryon.hpp \
# Modules/MContraction/DiscLoop.hpp \
# Modules/MContraction/Gamma3pt.hpp \
# Modules/MContraction/WardIdentity.hpp \
# Modules/MContraction/WeakHamiltonian.hpp \
# Modules/MContraction/WeakHamiltonianEye.hpp \
# Modules/MContraction/WeakHamiltonianNonEye.hpp \
# Modules/MContraction/WeakNeutral4ptDisc.hpp \
# Modules/MFermion/GaugeProp.hpp \
# Modules/MGauge/Load.hpp \
# Modules/MGauge/Random.hpp \
# Modules/MGauge/StochEm.hpp \
# Modules/MLoop/NoiseLoop.hpp \
# Modules/MScalar/ChargedProp.hpp \
# Modules/MScalar/FreeProp.hpp \
# Modules/MScalar/Scalar.hpp \
# Modules/MSink/Smear.hpp \
# Modules/MSolver/RBPrecCG.hpp \
# Modules/MSource/SeqConserved.hpp \
# Modules/MSource/SeqGamma.hpp \
# Modules/MSource/Wall.hpp \
# Modules/MSource/Z2.hpp \
# Modules/MUtilities/TestSeqConserved.hpp \
# Modules/MUtilities/TestSeqGamma.hpp
Modules/MContraction/Baryon.hpp \
Modules/MContraction/DiscLoop.hpp \
Modules/MContraction/Gamma3pt.hpp \
Modules/MContraction/Meson.hpp \
Modules/MContraction/WardIdentity.hpp \
Modules/MContraction/WeakHamiltonian.hpp \
Modules/MContraction/WeakHamiltonianEye.hpp \
Modules/MContraction/WeakHamiltonianNonEye.hpp \
Modules/MContraction/WeakNeutral4ptDisc.hpp \
Modules/MLoop/NoiseLoop.hpp \
Modules/MSink/Smear.hpp \
Modules/MSolver/RBPrecCG.hpp \
Modules/MSource/SeqConserved.hpp \
Modules/MSource/SeqGamma.hpp \
Modules/MSource/Wall.hpp \
Modules/MSource/Z2.hpp \
Modules/MUtilities/TestSeqConserved.hpp \
Modules/MUtilities/TestSeqGamma.hpp

40
tests/hadrons/Test_hadrons.hpp
View File

@ -118,7 +118,7 @@ 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)))
if (!(VirtualMachine::getInstance().hasModule(actionName)))
{
MAction::Wilson::Par actionPar;
actionPar.gauge = gaugeField;
@ -144,7 +144,7 @@ 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)))
if (!(VirtualMachine::getInstance().hasModule(actionName)))
{
MAction::DWF::Par actionPar;
actionPar.gauge = gaugeField;
@ -173,7 +173,7 @@ inline void makeDWFAction(Application &application, std::string actionName,
inline void makeRBPrecCGSolver(Application &application, std::string &solverName,
std::string &actionName, double residual = 1e-8)
{
if (!(Environment::getInstance().hasModule(solverName)))
if (!(VirtualMachine::getInstance().hasModule(solverName)))
{
MSolver::RBPrecCG::Par solverPar;
solverPar.action = actionName;
@ -195,7 +195,7 @@ 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)))
if (!(VirtualMachine::getInstance().hasModule(srcName)))
{
MSource::Point::Par pointPar;
pointPar.position = pos;
@ -219,7 +219,7 @@ inline void makeSequentialSource(Application &application, std::string srcName,
std::string mom = ZERO_MOM)
{
// If the source already exists, don't make the module again.
if (!(Environment::getInstance().hasModule(srcName)))
if (!(VirtualMachine::getInstance().hasModule(srcName)))
{
MSource::SeqGamma::Par seqPar;
seqPar.q = qSrc;
@ -255,7 +255,7 @@ inline void makeConservedSequentialSource(Application &application,
std::string mom = ZERO_MOM)
{
// If the source already exists, don't make the module again.
if (!(Environment::getInstance().hasModule(srcName)))
if (!(VirtualMachine::getInstance().hasModule(srcName)))
{
MSource::SeqConserved::Par seqPar;
seqPar.q = qSrc;
@ -280,7 +280,7 @@ inline void makeConservedSequentialSource(Application &application,
inline void makeNoiseSource(Application &application, std::string &srcName,
unsigned int tA, unsigned int tB)
{
if (!(Environment::getInstance().hasModule(srcName)))
if (!(VirtualMachine::getInstance().hasModule(srcName)))
{
MSource::Z2::Par noisePar;
noisePar.tA = tA;
@ -302,7 +302,7 @@ 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)))
if (!(VirtualMachine::getInstance().hasModule(srcName)))
{
MSource::Wall::Par wallPar;
wallPar.tW = tW;
@ -324,7 +324,7 @@ 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)))
if (!(VirtualMachine::getInstance().hasModule(sinkFnct)))
{
MSink::Point::Par pointPar;
pointPar.mom = mom;
@ -345,7 +345,7 @@ 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)))
if (!(VirtualMachine::getInstance().hasModule(smearedProp)))
{
MSink::Smear::Par smearPar;
smearPar.q = propName;
@ -367,7 +367,7 @@ 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)))
if (!(VirtualMachine::getInstance().hasModule(propName)))
{
MFermion::GaugeProp::Par quarkPar;
quarkPar.source = srcName;
@ -390,7 +390,7 @@ 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)))
if (!(VirtualMachine::getInstance().hasModule(propName)))
{
MLoop::NoiseLoop::Par loopPar;
loopPar.q = resName;
@ -421,7 +421,7 @@ inline void mesonContraction(Application &application,
std::string &sink,
std::string gammas = "<Gamma5 Gamma5>")
{
if (!(Environment::getInstance().hasModule(modName)))
if (!(VirtualMachine::getInstance().hasModule(modName)))
{
MContraction::Meson::Par mesPar;
mesPar.output = output;
@ -453,7 +453,7 @@ inline void gamma3ptContraction(Application &application, unsigned int npt,
Gamma::Algebra gamma = Gamma::Algebra::Identity)
{
std::string modName = std::to_string(npt) + "pt_" + label;
if (!(Environment::getInstance().hasModule(modName)))
if (!(VirtualMachine::getInstance().hasModule(modName)))
{
MContraction::Gamma3pt::Par gamma3ptPar;
gamma3ptPar.output = std::to_string(npt) + "pt/" + label;
@ -487,7 +487,7 @@ inline void weakContraction##top(Application &application, unsigned int npt,\
std::string &label, unsigned int tSnk = 0)\
{\
std::string modName = std::to_string(npt) + "pt_" + label;\
if (!(Environment::getInstance().hasModule(modName)))\
if (!(VirtualMachine::getInstance().hasModule(modName)))\
{\
MContraction::WeakHamiltonian##top::Par weakPar;\
weakPar.output = std::to_string(npt) + "pt/" + label;\
@ -521,7 +521,7 @@ inline void disc0Contraction(Application &application,
std::string &label)
{
std::string modName = "4pt_" + label;
if (!(Environment::getInstance().hasModule(modName)))
if (!(VirtualMachine::getInstance().hasModule(modName)))
{
MContraction::WeakNeutral4ptDisc::Par disc0Par;
disc0Par.output = "4pt/" + label;
@ -547,7 +547,7 @@ inline void discLoopContraction(Application &application,
std::string &q_loop, std::string &modName,
Gamma::Algebra gamma = Gamma::Algebra::Identity)
{
if (!(Environment::getInstance().hasModule(modName)))
if (!(VirtualMachine::getInstance().hasModule(modName)))
{
MContraction::DiscLoop::Par discPar;
discPar.output = "disc/" + modName;
@ -574,7 +574,7 @@ 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)))
if (!(VirtualMachine::getInstance().hasModule(modName)))
{
MContraction::WardIdentity::Par wiPar;
if (Ls > 1)
@ -613,7 +613,7 @@ inline void makeSeqCurrComparison(Application &application, std::string &modName
std::string &actionName, std::string &origin,
unsigned int t_J, unsigned int mu, Current curr)
{
if (!(Environment::getInstance().hasModule(modName)))
if (!(VirtualMachine::getInstance().hasModule(modName)))
{
MUtilities::TestSeqConserved::Par seqPar;
seqPar.q = propName;
@ -646,7 +646,7 @@ inline void makeSeqGamComparison(Application &application, std::string &modName,
std::string &origin, Gamma::Algebra gamma,
unsigned int t_g)
{
if (!(Environment::getInstance().hasModule(modName)))
if (!(VirtualMachine::getInstance().hasModule(modName)))
{
MUtilities::TestSeqGamma::Par seqPar;
seqPar.q = propName;