From a0d399e5ce32d8bf0f51cceaa01cf920b907bde8 Mon Sep 17 00:00:00 2001 From: Antonin Portelli Date: Fri, 18 May 2018 20:49:26 +0100 Subject: [PATCH] Hadrons: yet other attempts at EMT NPR --- extras/Hadrons/Modules.hpp | 2 + .../Modules/MScalarSUN/TimeMomProbe.cc | 11 + .../Modules/MScalarSUN/TimeMomProbe.hpp | 241 ++++++++++++++++++ .../Hadrons/Modules/MScalarSUN/TwoPoint.hpp | 31 --- .../Hadrons/Modules/MScalarSUN/TwoPointNPR.cc | 11 + .../Modules/MScalarSUN/TwoPointNPR.hpp | 182 +++++++++++++ extras/Hadrons/Modules/MScalarSUN/Utils.hpp | 22 ++ extras/Hadrons/modules.inc | 4 + 8 files changed, 473 insertions(+), 31 deletions(-) create mode 100644 extras/Hadrons/Modules/MScalarSUN/TimeMomProbe.cc create mode 100644 extras/Hadrons/Modules/MScalarSUN/TimeMomProbe.hpp create mode 100644 extras/Hadrons/Modules/MScalarSUN/TwoPointNPR.cc create mode 100644 extras/Hadrons/Modules/MScalarSUN/TwoPointNPR.hpp diff --git a/extras/Hadrons/Modules.hpp b/extras/Hadrons/Modules.hpp index 1a18eb3a..3b33eb6e 100644 --- a/extras/Hadrons/Modules.hpp +++ b/extras/Hadrons/Modules.hpp @@ -32,8 +32,10 @@ #include #include #include +#include #include #include +#include #include #include #include diff --git a/extras/Hadrons/Modules/MScalarSUN/TimeMomProbe.cc b/extras/Hadrons/Modules/MScalarSUN/TimeMomProbe.cc new file mode 100644 index 00000000..f22d08af --- /dev/null +++ b/extras/Hadrons/Modules/MScalarSUN/TimeMomProbe.cc @@ -0,0 +1,11 @@ +#include + +using namespace Grid; +using namespace Hadrons; +using namespace MScalarSUN; + +template class Grid::Hadrons::MScalarSUN::TTimeMomProbe>; +template class Grid::Hadrons::MScalarSUN::TTimeMomProbe>; +template class Grid::Hadrons::MScalarSUN::TTimeMomProbe>; +template class Grid::Hadrons::MScalarSUN::TTimeMomProbe>; +template class Grid::Hadrons::MScalarSUN::TTimeMomProbe>; diff --git a/extras/Hadrons/Modules/MScalarSUN/TimeMomProbe.hpp b/extras/Hadrons/Modules/MScalarSUN/TimeMomProbe.hpp new file mode 100644 index 00000000..4807e865 --- /dev/null +++ b/extras/Hadrons/Modules/MScalarSUN/TimeMomProbe.hpp @@ -0,0 +1,241 @@ +#ifndef Hadrons_MScalarSUN_TimeMomProbe_hpp_ +#define Hadrons_MScalarSUN_TimeMomProbe_hpp_ + +#include +#include +#include +#include + +BEGIN_HADRONS_NAMESPACE + +/****************************************************************************** + * n-point functions O(t,p)*tr(phi(t_1,p_1)*...*phi(t_n,p_n)) * + ******************************************************************************/ +BEGIN_MODULE_NAMESPACE(MScalarSUN) + +class TimeMomProbePar: Serializable +{ +public: + GRID_SERIALIZABLE_CLASS_MEMBERS(TimeMomProbePar, + std::string, field, + std::vector, op, + std::vector>, timeMom, + std::string, output); +}; + +class TimeMomProbeResult: Serializable +{ +public: + GRID_SERIALIZABLE_CLASS_MEMBERS(TimeMomProbeResult, + std::string, op, + std::vector>, timeMom, + std::vector, data); +}; + +template +class TTimeMomProbe: public Module +{ +public: + typedef typename SImpl::Field Field; + typedef typename SImpl::SiteField::scalar_object Site; + typedef typename SImpl::ComplexField ComplexField; + typedef std::vector SlicedOp; +public: + // constructor + TTimeMomProbe(const std::string name); + // destructor + virtual ~TTimeMomProbe(void) {}; + // dependency relation + virtual std::vector getInput(void); + virtual std::vector getOutput(void); + // setup + virtual void setup(void); + // execution + virtual void execute(void); +private: + void vectorModulo(std::vector &v); +}; + +MODULE_REGISTER_TMP(TimeMomProbeSU2, TTimeMomProbe>, MScalarSUN); +MODULE_REGISTER_TMP(TimeMomProbeSU3, TTimeMomProbe>, MScalarSUN); +MODULE_REGISTER_TMP(TimeMomProbeSU4, TTimeMomProbe>, MScalarSUN); +MODULE_REGISTER_TMP(TimeMomProbeSU5, TTimeMomProbe>, MScalarSUN); +MODULE_REGISTER_TMP(TimeMomProbeSU6, TTimeMomProbe>, MScalarSUN); + +/****************************************************************************** + * TTimeMomProbe implementation * + ******************************************************************************/ +// constructor ///////////////////////////////////////////////////////////////// +template +TTimeMomProbe::TTimeMomProbe(const std::string name) +: Module(name) +{} + +// dependencies/products /////////////////////////////////////////////////////// +template +std::vector TTimeMomProbe::getInput(void) +{ + std::vector in = par().op; + + in.push_back(par().field); + + return in; +} + +template +std::vector TTimeMomProbe::getOutput(void) +{ + std::vector out; + + return out; +} + +// setup /////////////////////////////////////////////////////////////////////// +template +void TTimeMomProbe::setup(void) +{ + envTmpLat(ComplexField, "ftBuf"); + envTmpLat(Field, "ftMatBuf"); +} + +// execution /////////////////////////////////////////////////////////////////// +// NB: time is direction 0 +template +void TTimeMomProbe::vectorModulo(std::vector &v) +{ + for (unsigned int mu = 0; mu < env().getNd(); ++mu) + { + auto d = env().getDim(mu); + v[mu] = ((v[mu] % d) + d) % d; + } +} + +template +void TTimeMomProbe::execute(void) +{ + const unsigned int nd = env().getNd(); + const unsigned int nt = env().getDim(0); + double partVol = 1.; + std::set> timeMomSet; + std::vector>> timeMom; + std::vector> transferMom; + FFT fft(env().getGrid()); + std::vector dMask(nd, 1); + std::vector result; + std::map> slicedOp; + std::vector slicedProbe; + auto &phi = envGet(Field, par().field); + + envGetTmp(ComplexField, ftBuf); + envGetTmp(Field, ftMatBuf); + dMask[0] = 0; + for (unsigned int mu = 1; mu < nd; ++mu) + { + partVol *= env().getDim(mu); + } + timeMom.resize(par().timeMom.size()); + for (unsigned int p = 0; p < timeMom.size(); ++p) + { + for (auto &tms: par().timeMom[p]) + { + std::vector tm = strToVec(tms); + + timeMom[p].push_back(tm); + timeMomSet.insert(tm); + } + transferMom.push_back(std::vector(nd - 1, 0)); + for (auto &tm: timeMom[p]) + { + for (unsigned int j = 1; j < nd; ++j) + { + transferMom[p][j - 1] -= tm[j]; + } + } + LOG(Message) << "Probe " << p << " (" << timeMom[p].size() << " points) : " << std::endl; + LOG(Message) << " phi(t_i, p_i) for (t_i, p_i) in " << timeMom[p] << std::endl; + LOG(Message) << " operator with momentum " << transferMom[p] << std::endl; + } + LOG(Message) << "FFT: field '" << par().field << "'" << std::endl; + fft.FFT_dim_mask(ftMatBuf, phi, dMask, FFT::forward); + slicedProbe.resize(timeMom.size()); + for (unsigned int p = 0; p < timeMom.size(); ++p) + { + std::vector qt; + + LOG(Message) << "Making probe " << p << std::endl; + slicedProbe[p].resize(nt); + for (unsigned int t = 0; t < nt; ++t) + { + Site acc; + + for (unsigned int i = 0; i < timeMom[p].size(); ++i) + { + Site buf; + + qt = timeMom[p][i]; + qt[0] += t; + vectorModulo(qt); + peekSite(buf, ftMatBuf, qt); + if (i == 0) + { + acc = buf; + } + else + { + acc *= buf; + } + } + slicedProbe[p][t] = TensorRemove(trace(acc)); + } + //std::cout << slicedProbe[p]<< std::endl; + } + for (auto &o: par().op) + { + auto &op = envGet(ComplexField, o); + + slicedOp[o].resize(transferMom.size()); + LOG(Message) << "FFT: operator '" << o << "'" << std::endl; + fft.FFT_dim_mask(ftBuf, op, dMask, FFT::forward); + //std::cout << ftBuf << std::endl; + for (unsigned int p = 0; p < transferMom.size(); ++p) + { + std::vector qt(nd, 0); + + for (unsigned int j = 1; j < nd; ++j) + { + qt[j] = transferMom[p][j - 1]; + } + slicedOp[o][p].resize(nt); + for (unsigned int t = 0; t < nt; ++t) + { + TComplex buf; + + qt[0] = t; + vectorModulo(qt); + peekSite(buf, ftBuf, qt); + slicedOp[o][p][t] = TensorRemove(buf); + } + //std::cout << ftBuf << std::endl; + //std::cout << slicedOp[o][p] << std::endl; + } + } + LOG(Message) << "Making correlators" << std::endl; + for (auto &o: par().op) + for (unsigned int p = 0; p < timeMom.size(); ++p) + { + TimeMomProbeResult r; + + LOG(Message) << " <" << o << " probe_" << p << ">" << std::endl; + r.op = o; + r.timeMom = timeMom[p]; + r.data = makeTwoPoint(slicedOp[o][p], slicedProbe[p], 1./partVol); + result.push_back(r); + } + saveResult(par().output, "timemomprobe", result); +} + +END_MODULE_NAMESPACE + +END_HADRONS_NAMESPACE + +#endif // Hadrons_MScalarSUN_TimeMomProbe_hpp_ diff --git a/extras/Hadrons/Modules/MScalarSUN/TwoPoint.hpp b/extras/Hadrons/Modules/MScalarSUN/TwoPoint.hpp index 8ee80242..2407f26f 100644 --- a/extras/Hadrons/Modules/MScalarSUN/TwoPoint.hpp +++ b/extras/Hadrons/Modules/MScalarSUN/TwoPoint.hpp @@ -79,12 +79,6 @@ public: virtual void setup(void); // execution virtual void execute(void); -private: - // make 2-pt function - template - std::vector makeTwoPoint(const std::vector &sink, - const std::vector &source, - const double factor = 1.); private: std::vector> mom_; }; @@ -221,31 +215,6 @@ void TTwoPoint::execute(void) saveResult(par().output, "twopt", result); } -// make 2-pt function ////////////////////////////////////////////////////////// -template -template -std::vector TTwoPoint::makeTwoPoint( - const std::vector &sink, - const std::vector &source, - const double factor) -{ - assert(sink.size() == source.size()); - - unsigned int nt = sink.size(); - std::vector res(nt, 0.); - - for (unsigned int dt = 0; dt < nt; ++dt) - { - for (unsigned int t = 0; t < nt; ++t) - { - res[dt] += sink[(t+dt)%nt]*adj(source[t]); - } - res[dt] *= factor/static_cast(nt); - } - - return res; -} - END_MODULE_NAMESPACE END_HADRONS_NAMESPACE diff --git a/extras/Hadrons/Modules/MScalarSUN/TwoPointNPR.cc b/extras/Hadrons/Modules/MScalarSUN/TwoPointNPR.cc new file mode 100644 index 00000000..70fb2445 --- /dev/null +++ b/extras/Hadrons/Modules/MScalarSUN/TwoPointNPR.cc @@ -0,0 +1,11 @@ +#include + +using namespace Grid; +using namespace Hadrons; +using namespace MScalarSUN; + +template class Grid::Hadrons::MScalarSUN::TTwoPointNPR>; +template class Grid::Hadrons::MScalarSUN::TTwoPointNPR>; +template class Grid::Hadrons::MScalarSUN::TTwoPointNPR>; +template class Grid::Hadrons::MScalarSUN::TTwoPointNPR>; +template class Grid::Hadrons::MScalarSUN::TTwoPointNPR>; diff --git a/extras/Hadrons/Modules/MScalarSUN/TwoPointNPR.hpp b/extras/Hadrons/Modules/MScalarSUN/TwoPointNPR.hpp new file mode 100644 index 00000000..4c2cdb0d --- /dev/null +++ b/extras/Hadrons/Modules/MScalarSUN/TwoPointNPR.hpp @@ -0,0 +1,182 @@ +#ifndef Hadrons_MScalarSUN_TwoPointNPR_hpp_ +#define Hadrons_MScalarSUN_TwoPointNPR_hpp_ + +#include +#include +#include +#include + +BEGIN_HADRONS_NAMESPACE + +/****************************************************************************** + * TwoPointNPR * + ******************************************************************************/ +BEGIN_MODULE_NAMESPACE(MScalarSUN) + +class TwoPointNPRPar: Serializable +{ +public: + GRID_SERIALIZABLE_CLASS_MEMBERS(TwoPointNPRPar, + std::vector, op, + std::string, field, + std::string, output); +}; + +class TwoPointNPRResult: Serializable +{ +public: + GRID_SERIALIZABLE_CLASS_MEMBERS(TwoPointNPRResult, + std::string, op, + std::vector, data); +}; + +template +class TTwoPointNPR: public Module +{ +public: + typedef typename SImpl::Field Field; + typedef typename SImpl::SiteField::scalar_object Site; + typedef typename SImpl::ComplexField ComplexField; +public: + // constructor + TTwoPointNPR(const std::string name); + // destructor + virtual ~TTwoPointNPR(void) {}; + // dependency relation + virtual std::vector getInput(void); + virtual std::vector getOutput(void); + // setup + virtual void setup(void); + // execution + virtual void execute(void); +}; + +MODULE_REGISTER_TMP(TwoPointNPRSU2, TTwoPointNPR>, MScalarSUN); +MODULE_REGISTER_TMP(TwoPointNPRSU3, TTwoPointNPR>, MScalarSUN); +MODULE_REGISTER_TMP(TwoPointNPRSU4, TTwoPointNPR>, MScalarSUN); +MODULE_REGISTER_TMP(TwoPointNPRSU5, TTwoPointNPR>, MScalarSUN); +MODULE_REGISTER_TMP(TwoPointNPRSU6, TTwoPointNPR>, MScalarSUN); + +/****************************************************************************** + * TTwoPointNPR implementation * + ******************************************************************************/ +// constructor ///////////////////////////////////////////////////////////////// +template +TTwoPointNPR::TTwoPointNPR(const std::string name) +: Module(name) +{} + +// dependencies/products /////////////////////////////////////////////////////// +template +std::vector TTwoPointNPR::getInput(void) +{ + std::vector in = par().op; + + in.push_back(par().field); + + return in; +} + +template +std::vector TTwoPointNPR::getOutput(void) +{ + std::vector out; + + return out; +} + +// setup /////////////////////////////////////////////////////////////////////// +template +void TTwoPointNPR::setup(void) +{ + const unsigned int nl = env().getDim(0); + + for (unsigned int mu = 1; mu < env().getNd(); ++mu) + { + if (nl != env().getDim(mu)) + { + HADRONS_ERROR(Size, "non-cubic grid"); + } + } + envTmpLat(ComplexField, "ftBuf"); + envTmpLat(Field, "ftMatBuf"); +} + +// execution /////////////////////////////////////////////////////////////////// +template +void TTwoPointNPR::execute(void) +{ + const unsigned int nd = env().getNd(); + const unsigned int nl = env().getDim(0); + FFT fft(env().getGrid()); + std::vector result; + TwoPointNPRResult twoPtp1, twoPtp2; + auto &phi = envGet(Field, par().field); + bool doTwoPt = true; + + envGetTmp(ComplexField, ftBuf); + envGetTmp(Field, ftMatBuf); + LOG(Message) << "FFT: field '" << par().field << "'" << std::endl; + fft.FFT_all_dim(ftMatBuf, phi, FFT::forward); + for (auto &opName: par().op) + { + auto &op = envGet(ComplexField, opName); + std::vector p1, p2, p; + Site phip1, phip2; + TComplex opp; + TwoPointNPRResult r; + + LOG(Message) << "FFT: operator '" << opName << "'" << std::endl; + fft.FFT_all_dim(ftBuf, op, FFT::forward); + LOG(Message) << "Generating vertex function" << std::endl; + r.op = opName; + r.data.resize(nl); + if (doTwoPt) + { + twoPtp1.op = "phi_prop_p1"; + twoPtp1.data.resize(nl); + twoPtp2.op = "phi_prop_p2"; + twoPtp2.data.resize(nl); + } + for (unsigned int n = 0; n < nl; ++n) + { + p1.assign(nd, 0); + p2.assign(nd, 0); + p.assign(nd, 0); + // non-exceptional RI/SMOM kinematic + // p1 = mu*(1,1,0): in mom + // p2 = mu*(0,1,1): out mom + // p = p1 - p2 = mu*(1,0,-1) + // mu = 2*n*pi/L + p1[0] = n; + p1[1] = n; + p2[1] = n; + p2[2] = n; + p[0] = n; + p[2] = (nl - n) % nl; + peekSite(phip1, ftMatBuf, p1); + peekSite(phip2, ftMatBuf, p2); + peekSite(opp, ftBuf, p); + if (doTwoPt) + { + twoPtp1.data[n] = TensorRemove(trace(phip1*adj(phip1))); + twoPtp2.data[n] = TensorRemove(trace(phip2*adj(phip2))); + } + r.data[n] = TensorRemove(trace(phip2*adj(phip1))*opp); + } + if (doTwoPt) + { + result.push_back(twoPtp1); + result.push_back(twoPtp2); + } + result.push_back(r); + doTwoPt = false; + } + saveResult(par().output, "twoptnpr", result); +} + +END_MODULE_NAMESPACE + +END_HADRONS_NAMESPACE + +#endif // Hadrons_MScalarSUN_TwoPointNPR_hpp_ diff --git a/extras/Hadrons/Modules/MScalarSUN/Utils.hpp b/extras/Hadrons/Modules/MScalarSUN/Utils.hpp index 37a9e137..68009f5e 100644 --- a/extras/Hadrons/Modules/MScalarSUN/Utils.hpp +++ b/extras/Hadrons/Modules/MScalarSUN/Utils.hpp @@ -89,6 +89,28 @@ inline void dmuAcc(Field &out, const Field &in, const unsigned int mu, const Dif } } +template +std::vector makeTwoPoint(const std::vector &sink, + const std::vector &source, + const double factor = 1.) +{ + assert(sink.size() == source.size()); + + unsigned int nt = sink.size(); + std::vector res(nt, 0.); + + for (unsigned int dt = 0; dt < nt; ++dt) + { + for (unsigned int t = 0; t < nt; ++t) + { + res[dt] += trace(sink[(t+dt)%nt]*source[t]); + } + res[dt] *= factor/static_cast(nt); + } + + return res; +} + inline std::string varName(const std::string name, const std::string suf) { return name + "_" + suf; diff --git a/extras/Hadrons/modules.inc b/extras/Hadrons/modules.inc index c2e84086..a1f34b95 100644 --- a/extras/Hadrons/modules.inc +++ b/extras/Hadrons/modules.inc @@ -32,6 +32,7 @@ modules_cc =\ Modules/MAction/DWF.cc \ Modules/MScalarSUN/TrPhi.cc \ Modules/MScalarSUN/Grad.cc \ + Modules/MScalarSUN/TimeMomProbe.cc \ Modules/MScalarSUN/TrMag.cc \ Modules/MScalarSUN/TrKinetic.cc \ Modules/MScalarSUN/EMT.cc \ @@ -39,6 +40,7 @@ modules_cc =\ Modules/MScalarSUN/TransProj.cc \ Modules/MScalarSUN/StochFreeField.cc \ Modules/MScalarSUN/TwoPoint.cc \ + Modules/MScalarSUN/TwoPointNPR.cc \ Modules/MScalarSUN/Div.cc \ Modules/MIO/LoadEigenPack.cc \ Modules/MIO/LoadBinary.cc \ @@ -80,8 +82,10 @@ modules_hpp =\ Modules/MAction/WilsonClover.hpp \ Modules/MAction/ZMobiusDWF.hpp \ Modules/MScalarSUN/StochFreeField.hpp \ + Modules/MScalarSUN/TwoPointNPR.hpp \ Modules/MScalarSUN/ShiftProbe.hpp \ Modules/MScalarSUN/Div.hpp \ + Modules/MScalarSUN/TimeMomProbe.hpp \ Modules/MScalarSUN/TrMag.hpp \ Modules/MScalarSUN/EMT.hpp \ Modules/MScalarSUN/TwoPoint.hpp \