From 6e4a06e180f7500df13ceea362b71294b8da74ff Mon Sep 17 00:00:00 2001 From: Antonin Portelli Date: Thu, 20 Oct 2016 15:04:00 +0100 Subject: [PATCH 01/53] qed-fvol: initial commit --- Makefile.am | 2 +- configure.ac | 2 ++ programs/Makefile.am | 1 + programs/qed-fvol/Global.cc | 11 +++++++++ programs/qed-fvol/Global.hpp | 42 +++++++++++++++++++++++++++++++++++ programs/qed-fvol/Makefile.am | 9 ++++++++ programs/qed-fvol/qed-fvol.cc | 36 ++++++++++++++++++++++++++++++ 7 files changed, 102 insertions(+), 1 deletion(-) create mode 100644 programs/Makefile.am create mode 100644 programs/qed-fvol/Global.cc create mode 100644 programs/qed-fvol/Global.hpp create mode 100644 programs/qed-fvol/Makefile.am create mode 100644 programs/qed-fvol/qed-fvol.cc diff --git a/Makefile.am b/Makefile.am index 90c5cd71..8cc860a9 100644 --- a/Makefile.am +++ b/Makefile.am @@ -1,5 +1,5 @@ # additional include paths necessary to compile the C++ library -SUBDIRS = lib benchmarks tests +SUBDIRS = lib benchmarks tests programs AM_CXXFLAGS += -I$(top_builddir)/include ACLOCAL_AMFLAGS = -I m4 diff --git a/configure.ac b/configure.ac index 7bcdc49f..81ced467 100644 --- a/configure.ac +++ b/configure.ac @@ -326,6 +326,8 @@ AC_CONFIG_FILES(tests/hmc/Makefile) AC_CONFIG_FILES(tests/solver/Makefile) AC_CONFIG_FILES(tests/qdpxx/Makefile) AC_CONFIG_FILES(benchmarks/Makefile) +AC_CONFIG_FILES(programs/Makefile) +AC_CONFIG_FILES(programs/qed-fvol/Makefile) AC_OUTPUT echo " diff --git a/programs/Makefile.am b/programs/Makefile.am new file mode 100644 index 00000000..ff7f6584 --- /dev/null +++ b/programs/Makefile.am @@ -0,0 +1 @@ +SUBDIRS = qed-fvol diff --git a/programs/qed-fvol/Global.cc b/programs/qed-fvol/Global.cc new file mode 100644 index 00000000..57ed97cc --- /dev/null +++ b/programs/qed-fvol/Global.cc @@ -0,0 +1,11 @@ +#include + +using namespace Grid; +using namespace QCD; +using namespace QedFVol; + +QedFVolLogger QedFVol::QedFVolLogError(1,"Error"); +QedFVolLogger QedFVol::QedFVolLogWarning(1,"Warning"); +QedFVolLogger QedFVol::QedFVolLogMessage(1,"Message"); +QedFVolLogger QedFVol::QedFVolLogIterative(1,"Iterative"); +QedFVolLogger QedFVol::QedFVolLogDebug(1,"Debug"); diff --git a/programs/qed-fvol/Global.hpp b/programs/qed-fvol/Global.hpp new file mode 100644 index 00000000..7f07200d --- /dev/null +++ b/programs/qed-fvol/Global.hpp @@ -0,0 +1,42 @@ +#ifndef QedFVol_Global_hpp_ +#define QedFVol_Global_hpp_ + +#include + +#define BEGIN_QEDFVOL_NAMESPACE \ +namespace Grid {\ +using namespace QCD;\ +namespace QedFVol {\ +using Grid::operator<<; +#define END_QEDFVOL_NAMESPACE }} + +/* the 'using Grid::operator<<;' statement prevents a very nasty compilation + * error with GCC (clang compiles fine without it). + */ + +BEGIN_QEDFVOL_NAMESPACE + +class QedFVolLogger: public Logger +{ +public: + QedFVolLogger(int on, std::string nm): Logger("QedFVol", on, nm, + GridLogColours, "BLACK"){}; +}; + +#define LOG(channel) std::cout << QedFVolLog##channel +#define QEDFVOL_ERROR(msg)\ +LOG(Error) << msg << " (" << __FUNCTION__ << " at " << __FILE__ << ":"\ + << __LINE__ << ")" << std::endl;\ +abort(); + +#define DEBUG_VAR(var) LOG(Debug) << #var << "= " << (var) << std::endl; + +extern QedFVolLogger QedFVolLogError; +extern QedFVolLogger QedFVolLogWarning; +extern QedFVolLogger QedFVolLogMessage; +extern QedFVolLogger QedFVolLogIterative; +extern QedFVolLogger QedFVolLogDebug; + +END_QEDFVOL_NAMESPACE + +#endif // QedFVol_Global_hpp_ diff --git a/programs/qed-fvol/Makefile.am b/programs/qed-fvol/Makefile.am new file mode 100644 index 00000000..cd762e94 --- /dev/null +++ b/programs/qed-fvol/Makefile.am @@ -0,0 +1,9 @@ +AM_CXXFLAGS += -I$(top_srcdir)/programs -I../$(top_srcdir)/programs + +bin_PROGRAMS = qed-fvol + +qed_fvol_SOURCES = \ + qed-fvol.cc \ + Global.cc + +qed_fvol_LDADD = -lGrid diff --git a/programs/qed-fvol/qed-fvol.cc b/programs/qed-fvol/qed-fvol.cc new file mode 100644 index 00000000..bb3204c6 --- /dev/null +++ b/programs/qed-fvol/qed-fvol.cc @@ -0,0 +1,36 @@ +#include + +using namespace Grid; +using namespace QCD; +using namespace QedFVol; + +int main(int argc, char *argv[]) +{ + // parse command line + std::string parameterFileName; + + if (argc < 2) + { + std::cerr << "usage: " << argv[0] << " [Grid options]"; + std::cerr << std::endl; + std::exit(EXIT_FAILURE); + } + parameterFileName = argv[1]; + + // initialization + Grid_init(&argc, &argv); + QedFVolLogError.Active(GridLogError.isActive()); + QedFVolLogWarning.Active(GridLogWarning.isActive()); + QedFVolLogMessage.Active(GridLogMessage.isActive()); + QedFVolLogIterative.Active(GridLogIterative.isActive()); + QedFVolLogDebug.Active(GridLogDebug.isActive()); + LOG(Message) << "Grid initialized" << std::endl; + + + + // epilogue + LOG(Message) << "Grid is finalizing now" << std::endl; + Grid_finalize(); + + return EXIT_SUCCESS; +} From 0d889b70410bfdaf70b5cbffe2fb92157943cc03 Mon Sep 17 00:00:00 2001 From: Antonin Portelli Date: Fri, 21 Oct 2016 15:21:32 +0100 Subject: [PATCH 02/53] QedFVol: first attempt at generating a QED field --- programs/qed-fvol/qed-fvol.cc | 38 ++++++++++++++++++++++++++++++++++- 1 file changed, 37 insertions(+), 1 deletion(-) diff --git a/programs/qed-fvol/qed-fvol.cc b/programs/qed-fvol/qed-fvol.cc index bb3204c6..53e01de9 100644 --- a/programs/qed-fvol/qed-fvol.cc +++ b/programs/qed-fvol/qed-fvol.cc @@ -4,6 +4,31 @@ using namespace Grid; using namespace QCD; using namespace QedFVol; +template +class QedGimpl +{ +public: + typedef S Simd; + + template + using iImplGaugeLink = iScalar>>; + template + using iImplGaugeField = iVector>, Nd>; + + typedef iImplGaugeLink SiteGaugeLink; + typedef iImplGaugeField SiteGaugeField; + + typedef Lattice GaugeLinkField; // bit ugly naming; polarised + // gauge field, lorentz... all + // ugly + typedef Lattice GaugeField; +}; + +typedef QedGimpl QedGimplR; +typedef Photon PhotonR; +typedef PhotonR::GaugeField EmField; +typedef PhotonR::GaugeLinkField EmComp; + int main(int argc, char *argv[]) { // parse command line @@ -26,8 +51,19 @@ int main(int argc, char *argv[]) QedFVolLogDebug.Active(GridLogDebug.isActive()); LOG(Message) << "Grid initialized" << std::endl; + // QED stuff + std::vector latt_size = GridDefaultLatt(); + std::vector simd_layout = GridDefaultSimd(4, vComplex::Nsimd()); + std::vector mpi_layout = GridDefaultMpi(); + GridCartesian grid(latt_size,simd_layout,mpi_layout); + GridParallelRNG pRNG(&grid); + PhotonR photon(PhotonR::Gauge::Feynman, + PhotonR::ZmScheme::QedL); + EmField a(&grid); + + pRNG.SeedRandomDevice(); + photon.StochasticField(a, pRNG); - // epilogue LOG(Message) << "Grid is finalizing now" << std::endl; Grid_finalize(); From 3ab4c8c0bbde6a572d41074405c2baa8e9a0119c Mon Sep 17 00:00:00 2001 From: James Harrison Date: Tue, 25 Oct 2016 13:32:02 +0100 Subject: [PATCH 03/53] QedFVol: calculate plaquette and 2x2 Wilson loop of stochastic QED field --- programs/qed-fvol/WilsonLoops.h | 167 ++++++++++++++++++++++++++++++++ programs/qed-fvol/qed-fvol.cc | 45 +++++++++ 2 files changed, 212 insertions(+) create mode 100644 programs/qed-fvol/WilsonLoops.h diff --git a/programs/qed-fvol/WilsonLoops.h b/programs/qed-fvol/WilsonLoops.h new file mode 100644 index 00000000..610fdc7b --- /dev/null +++ b/programs/qed-fvol/WilsonLoops.h @@ -0,0 +1,167 @@ +#ifndef QEDFVOL_WILSONLOOPS_H +#define QEDFVOL_WILSONLOOPS_H + +#include + +BEGIN_QEDFVOL_NAMESPACE + +template class WilsonLoops : public Gimpl { +public: + INHERIT_GIMPL_TYPES(Gimpl); + + typedef typename Gimpl::GaugeLinkField GaugeMat; + typedef typename Gimpl::GaugeField GaugeLorentz; + + ////////////////////////////////////////////////// + // directed plaquette oriented in mu,nu plane + ////////////////////////////////////////////////// + static void dirPlaquette(GaugeMat &plaq, const std::vector &U, + const int mu, const int nu) { + // Annoyingly, must use either scope resolution to find dependent base + // class, + // or this-> ; there is no "this" in a static method. This forces explicit + // Gimpl scope + // resolution throughout the usage in this file, and rather defeats the + // purpose of deriving + // from Gimpl. + plaq = Gimpl::CovShiftBackward( + U[mu], mu, Gimpl::CovShiftBackward( + U[nu], nu, Gimpl::CovShiftForward(U[mu], mu, U[nu]))); + } + ////////////////////////////////////////////////// + // trace of directed plaquette oriented in mu,nu plane + ////////////////////////////////////////////////// + static void traceDirPlaquette(LatticeComplex &plaq, + const std::vector &U, const int mu, + const int nu) { + GaugeMat sp(U[0]._grid); + dirPlaquette(sp, U, mu, nu); + plaq = trace(sp); + } + ////////////////////////////////////////////////// + // sum over all planes of plaquette + ////////////////////////////////////////////////// + static void sitePlaquette(LatticeComplex &Plaq, + const std::vector &U) { + LatticeComplex sitePlaq(U[0]._grid); + Plaq = zero; + for (int mu = 1; mu < Nd; mu++) { + for (int nu = 0; nu < mu; nu++) { + traceDirPlaquette(sitePlaq, U, mu, nu); + Plaq = Plaq + sitePlaq; + } + } + } + ////////////////////////////////////////////////// + // sum over all x,y,z,t and over all planes of plaquette + ////////////////////////////////////////////////// + static RealD sumPlaquette(const GaugeLorentz &Umu) { + std::vector U(4, Umu._grid); + + for (int mu = 0; mu < Nd; mu++) { + U[mu] = PeekIndex(Umu, mu); + } + + LatticeComplex Plaq(Umu._grid); + + sitePlaquette(Plaq, U); + + TComplex Tp = sum(Plaq); + Complex p = TensorRemove(Tp); + return p.real(); + } + ////////////////////////////////////////////////// + // average over all x,y,z,t and over all planes of plaquette + ////////////////////////////////////////////////// + static RealD avgPlaquette(const GaugeLorentz &Umu) { + RealD sumplaq = sumPlaquette(Umu); + double vol = Umu._grid->gSites(); + double faces = (1.0 * Nd * (Nd - 1)) / 2.0; + return sumplaq / vol / faces / Nc; // Nd , Nc dependent... FIXME + } + + ////////////////////////////////////////////////// + // Wilson loop of size (R1, R2), oriented in mu,nu plane + ////////////////////////////////////////////////// + static void wilsonLoop(GaugeMat &wl, const std::vector &U, + const int Rmu, const int Rnu, + const int mu, const int nu) { + wl = U[nu]; + + for(int i = 0; i < Rnu-1; i++){ + wl = Gimpl::CovShiftForward(U[nu], nu, wl); + } + + for(int i = 0; i < Rmu; i++){ + wl = Gimpl::CovShiftForward(U[mu], mu, wl); + } + + for(int i = 0; i < Rnu; i++){ + wl = Gimpl::CovShiftBackward(U[nu], nu, wl); + } + + for(int i = 0; i < Rmu; i++){ + wl = Gimpl::CovShiftBackward(U[mu], mu, wl); + } + } + ////////////////////////////////////////////////// + // trace of Wilson Loop oriented in mu,nu plane + ////////////////////////////////////////////////// + static void traceWilsonLoop(LatticeComplex &wl, + const std::vector &U, + const int Rmu, const int Rnu, + const int mu, const int nu) { + GaugeMat sp(U[0]._grid); + WilsonLoop(sp, U, Rmu, Rnu, mu, nu); + wl = trace(sp); + } + ////////////////////////////////////////////////// + // sum over all planes of Wilson loop + ////////////////////////////////////////////////// + static void siteWilsonLoop(LatticeComplex &Wl, + const std::vector &U + const int R1, const int R2) { + LatticeComplex siteWl(U[0]._grid); + Wl = zero; + for (int mu = 1; mu < Nd; mu++) { + for (int nu = 0; nu < mu; nu++) { + traceWilsonLoop(siteWl, U, R1, R2, mu, nu); + Wl = Wl + siteWl; + traceWilsonLoop(siteWl, U, R2, R1, mu, nu); + Wl = Wl + siteWl; + } + } + } + ////////////////////////////////////////////////// + // sum over all x,y,z,t and over all planes of Wilson loop + ////////////////////////////////////////////////// + static RealD sumWilsonLoop(const GaugeLorentz &Umu, + const int R1, const int R2) { + std::vector U(4, Umu._grid); + + for (int mu = 0; mu < Nd; mu++) { + U[mu] = PeekIndex(Umu, mu); + } + + LatticeComplex Wl(Umu._grid); + + siteWilsonLoop(Wl, U, R1, R2); + + TComplex Tp = sum(Wl); + Complex p = TensorRemove(Tp); + return p.real(); + } + ////////////////////////////////////////////////// + // average over all x,y,z,t and over all planes of Wilson loop + ////////////////////////////////////////////////// + static RealD avgPlaquette(const GaugeLorentz &Umu, + const int R1, const int R2) { + RealD sumWl = sumWilsonLoop(Umu); + double vol = Umu._grid->gSites(); + double faces = 1.0 * Nd * (Nd - 1); + return sumWl / vol / faces / Nc; // Nd , Nc dependent... FIXME + } + +END_QEDFVOL_NAMESPACE + +#endif // QEDFVOL_WILSONLOOPS_H \ No newline at end of file diff --git a/programs/qed-fvol/qed-fvol.cc b/programs/qed-fvol/qed-fvol.cc index 53e01de9..02c36a67 100644 --- a/programs/qed-fvol/qed-fvol.cc +++ b/programs/qed-fvol/qed-fvol.cc @@ -64,6 +64,51 @@ int main(int argc, char *argv[]) pRNG.SeedRandomDevice(); photon.StochasticField(a, pRNG); + // Calculate log of plaquette + EmComp plaqA(&grid); + EmComp wlA(&grid); + EmComp tmp(&grid); + std::vector a_comp(4, &grid); + + for (int dir = 0; dir < Nd; dir++) { + a_comp[dir] = PeekIndex(a, dir); + } + + plaqA = zero; + wlA = zero; + + for(int mu = 1; mu < Nd; mu++) { + for(int nu = 0; nu < mu; nu++) { + tmp = a_comp[mu] + Cshift(a_comp[nu], mu, 1) - Cshift(a_comp[mu], nu, 1) - a_comp[nu]; + plaqA = plaqA + cos(tmp); + + tmp = a_comp[mu] + Cshift(a_comp[mu], mu, 1) + + Cshift(a_comp[nu], mu, 2) + Cshift(Cshift(a_comp[nu], mu, 2), nu, 1) + - Cshift(Cshift(a_comp[mu], nu, 2), mu, 1) - Cshift(a_comp[mu], nu, 2) + - Cshift(a_comp[nu], nu, 1) - a_comp[nu]; + wlA = wlA + cos(tmp); + } + } + + double vol = grid.gSites(); + double faces = (1.0 * Nd * (Nd - 1)) / 2.0; + + Complex avgPlaqA = sum(trace(plaqA)); + avgPlaqA = avgPlaqA / vol / faces; + + Complex avgWlA = sum(trace(wlA)); + avgWlA = avgWlA / vol / faces; + + TComplex tplaqsite; + LatticeComplex plaqtrace = trace(plaqA); + std::vector site0 = {0,0,0,0}; + peekSite(tplaqsite, plaqtrace, site0); + Complex plaqsite = TensorRemove(tplaqsite); + + LOG(Message) << "Plaquette average: " << avgPlaqA << std::endl; + LOG(Message) << "2x2 Wilson Loop average: " << avgWlA << std::endl; + LOG(Message) << "Plaquette (one site): " << plaqsite / faces << std::endl; + // epilogue LOG(Message) << "Grid is finalizing now" << std::endl; Grid_finalize(); From 78c7bcee36f7d937c8f5c6afe0f2088f85ebda51 Mon Sep 17 00:00:00 2001 From: James Harrison Date: Tue, 1 Nov 2016 13:30:11 +0000 Subject: [PATCH 04/53] QedFVol: Change variables of type "double" to type "Real". --- programs/qed-fvol/qed-fvol.cc | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/programs/qed-fvol/qed-fvol.cc b/programs/qed-fvol/qed-fvol.cc index 02c36a67..fd780edf 100644 --- a/programs/qed-fvol/qed-fvol.cc +++ b/programs/qed-fvol/qed-fvol.cc @@ -90,8 +90,8 @@ int main(int argc, char *argv[]) } } - double vol = grid.gSites(); - double faces = (1.0 * Nd * (Nd - 1)) / 2.0; + Real vol = grid.gSites(); + Real faces = (1.0 * Nd * (Nd - 1)) / 2.0; Complex avgPlaqA = sum(trace(plaqA)); avgPlaqA = avgPlaqA / vol / faces; From c30d96ea5097ab1760a3f0a6ea1aed8ac1e6142b Mon Sep 17 00:00:00 2001 From: James Harrison Date: Wed, 9 Nov 2016 11:06:20 +0000 Subject: [PATCH 05/53] QedFVol: x86intrin.h namespace fix --- lib/PerfCount.h | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/lib/PerfCount.h b/lib/PerfCount.h index 9ac58883..5ab07c02 100644 --- a/lib/PerfCount.h +++ b/lib/PerfCount.h @@ -43,6 +43,9 @@ Author: paboyle #else #include #endif +#ifdef __x86_64__ +#include +#endif namespace Grid { @@ -86,7 +89,6 @@ inline uint64_t cyclecount(void){ return tmp; } #elif defined __x86_64__ -#include inline uint64_t cyclecount(void){ return __rdtsc(); // unsigned int dummy; From cf167d0cd1c561bed3557eaf89350b8d8eb8d9b1 Mon Sep 17 00:00:00 2001 From: James Harrison Date: Mon, 14 Nov 2016 17:02:29 +0000 Subject: [PATCH 06/53] QedFVol: implement exponentiation of photon field --- programs/qed-fvol/WilsonLoops.h | 19 ++++++++++--------- programs/qed-fvol/qed-fvol.cc | 32 +++++++++++++++++++++++++------- 2 files changed, 35 insertions(+), 16 deletions(-) diff --git a/programs/qed-fvol/WilsonLoops.h b/programs/qed-fvol/WilsonLoops.h index 610fdc7b..c40fbaf3 100644 --- a/programs/qed-fvol/WilsonLoops.h +++ b/programs/qed-fvol/WilsonLoops.h @@ -5,7 +5,7 @@ BEGIN_QEDFVOL_NAMESPACE -template class WilsonLoops : public Gimpl { +template class NewWilsonLoops : public Gimpl { public: INHERIT_GIMPL_TYPES(Gimpl); @@ -55,7 +55,7 @@ public: ////////////////////////////////////////////////// // sum over all x,y,z,t and over all planes of plaquette ////////////////////////////////////////////////// - static RealD sumPlaquette(const GaugeLorentz &Umu) { + static Real sumPlaquette(const GaugeLorentz &Umu) { std::vector U(4, Umu._grid); for (int mu = 0; mu < Nd; mu++) { @@ -73,8 +73,8 @@ public: ////////////////////////////////////////////////// // average over all x,y,z,t and over all planes of plaquette ////////////////////////////////////////////////// - static RealD avgPlaquette(const GaugeLorentz &Umu) { - RealD sumplaq = sumPlaquette(Umu); + static Real avgPlaquette(const GaugeLorentz &Umu) { + Real sumplaq = sumPlaquette(Umu); double vol = Umu._grid->gSites(); double faces = (1.0 * Nd * (Nd - 1)) / 2.0; return sumplaq / vol / faces / Nc; // Nd , Nc dependent... FIXME @@ -112,14 +112,14 @@ public: const int Rmu, const int Rnu, const int mu, const int nu) { GaugeMat sp(U[0]._grid); - WilsonLoop(sp, U, Rmu, Rnu, mu, nu); + wilsonLoop(sp, U, Rmu, Rnu, mu, nu); wl = trace(sp); } ////////////////////////////////////////////////// // sum over all planes of Wilson loop ////////////////////////////////////////////////// static void siteWilsonLoop(LatticeComplex &Wl, - const std::vector &U + const std::vector &U, const int R1, const int R2) { LatticeComplex siteWl(U[0]._grid); Wl = zero; @@ -135,7 +135,7 @@ public: ////////////////////////////////////////////////// // sum over all x,y,z,t and over all planes of Wilson loop ////////////////////////////////////////////////// - static RealD sumWilsonLoop(const GaugeLorentz &Umu, + static Real sumWilsonLoop(const GaugeLorentz &Umu, const int R1, const int R2) { std::vector U(4, Umu._grid); @@ -154,13 +154,14 @@ public: ////////////////////////////////////////////////// // average over all x,y,z,t and over all planes of Wilson loop ////////////////////////////////////////////////// - static RealD avgPlaquette(const GaugeLorentz &Umu, + static Real avgWilsonLoop(const GaugeLorentz &Umu, const int R1, const int R2) { - RealD sumWl = sumWilsonLoop(Umu); + Real sumWl = sumWilsonLoop(Umu, R1, R2); double vol = Umu._grid->gSites(); double faces = 1.0 * Nd * (Nd - 1); return sumWl / vol / faces / Nc; // Nd , Nc dependent... FIXME } +}; END_QEDFVOL_NAMESPACE diff --git a/programs/qed-fvol/qed-fvol.cc b/programs/qed-fvol/qed-fvol.cc index fd780edf..68705b8f 100644 --- a/programs/qed-fvol/qed-fvol.cc +++ b/programs/qed-fvol/qed-fvol.cc @@ -1,4 +1,5 @@ #include +#include using namespace Grid; using namespace QCD; @@ -24,10 +25,11 @@ public: typedef Lattice GaugeField; }; -typedef QedGimpl QedGimplR; -typedef Photon PhotonR; -typedef PhotonR::GaugeField EmField; -typedef PhotonR::GaugeLinkField EmComp; +typedef QedGimpl QedGimplR; +typedef PeriodicGaugeImpl QedPeriodicGimplR; +typedef Photon PhotonR; +typedef PhotonR::GaugeField EmField; +typedef PhotonR::GaugeLinkField EmComp; int main(int argc, char *argv[]) { @@ -60,11 +62,18 @@ int main(int argc, char *argv[]) PhotonR photon(PhotonR::Gauge::Feynman, PhotonR::ZmScheme::QedL); EmField a(&grid); + EmField expA(&grid); + + Real avgPlaqAexp, avgWl2x2Aexp; pRNG.SeedRandomDevice(); photon.StochasticField(a, pRNG); - // Calculate log of plaquette + // Exponentiate photon field + Complex imag_unit(0, 1); + expA = exp(imag_unit*0.5*(a+conjugate(a))); + + // Calculate plaquette from photon field EmComp plaqA(&grid); EmComp wlA(&grid); EmComp tmp(&grid); @@ -105,8 +114,17 @@ int main(int argc, char *argv[]) peekSite(tplaqsite, plaqtrace, site0); Complex plaqsite = TensorRemove(tplaqsite); - LOG(Message) << "Plaquette average: " << avgPlaqA << std::endl; - LOG(Message) << "2x2 Wilson Loop average: " << avgWlA << std::endl; + // Calculate plaquette from exponentiated photon field + avgPlaqAexp = NewWilsonLoops::avgPlaquette(expA); + avgWl2x2Aexp = NewWilsonLoops::avgWilsonLoop(expA, 2, 2); + + avgPlaqAexp = avgPlaqAexp*3; + avgWl2x2Aexp = avgWl2x2Aexp*3; + + LOG(Message) << "Plaquette average (from A): " << avgPlaqA << std::endl; + LOG(Message) << "Plaquette average (from exp(A)): " << avgPlaqAexp << std::endl; + LOG(Message) << "2x2 Wilson Loop average (from A): " << avgWlA << std::endl; + LOG(Message) << "2x2 Wilson Loop average (from exp(A)): " << avgWl2x2Aexp << std::endl; LOG(Message) << "Plaquette (one site): " << plaqsite / faces << std::endl; // epilogue From f4ebea3381046026276864f3f908cb10b114d6a5 Mon Sep 17 00:00:00 2001 From: James Harrison Date: Mon, 14 Nov 2016 17:51:53 +0000 Subject: [PATCH 07/53] QedFVol: add functions for computing spatial and timelike Wilson loops --- programs/qed-fvol/WilsonLoops.h | 117 +++++++++++++++++++++++++++++--- programs/qed-fvol/qed-fvol.cc | 8 ++- 2 files changed, 114 insertions(+), 11 deletions(-) diff --git a/programs/qed-fvol/WilsonLoops.h b/programs/qed-fvol/WilsonLoops.h index c40fbaf3..98db6b7a 100644 --- a/programs/qed-fvol/WilsonLoops.h +++ b/programs/qed-fvol/WilsonLoops.h @@ -45,7 +45,7 @@ public: const std::vector &U) { LatticeComplex sitePlaq(U[0]._grid); Plaq = zero; - for (int mu = 1; mu < Nd; mu++) { + for (int mu = 1; mu < U[0]._grid->_ndimension; mu++) { for (int nu = 0; nu < mu; nu++) { traceDirPlaquette(sitePlaq, U, mu, nu); Plaq = Plaq + sitePlaq; @@ -58,7 +58,7 @@ public: static Real sumPlaquette(const GaugeLorentz &Umu) { std::vector U(4, Umu._grid); - for (int mu = 0; mu < Nd; mu++) { + for (int mu = 0; mu < Umu._grid->_ndimension; mu++) { U[mu] = PeekIndex(Umu, mu); } @@ -74,10 +74,11 @@ public: // average over all x,y,z,t and over all planes of plaquette ////////////////////////////////////////////////// static Real avgPlaquette(const GaugeLorentz &Umu) { + int ndim = Umu._grid->_ndimension; Real sumplaq = sumPlaquette(Umu); - double vol = Umu._grid->gSites(); - double faces = (1.0 * Nd * (Nd - 1)) / 2.0; - return sumplaq / vol / faces / Nc; // Nd , Nc dependent... FIXME + Real vol = Umu._grid->gSites(); + Real faces = (1.0 * ndim * (ndim - 1)) / 2.0; + return sumplaq / vol / faces / Nc; // Nc dependent... FIXME } ////////////////////////////////////////////////// @@ -123,7 +124,42 @@ public: const int R1, const int R2) { LatticeComplex siteWl(U[0]._grid); Wl = zero; - for (int mu = 1; mu < Nd; mu++) { + for (int mu = 1; mu < U[0]._grid->_ndimension; mu++) { + for (int nu = 0; nu < mu; nu++) { + traceWilsonLoop(siteWl, U, R1, R2, mu, nu); + Wl = Wl + siteWl; + traceWilsonLoop(siteWl, U, R2, R1, mu, nu); + Wl = Wl + siteWl; + } + } + } + ////////////////////////////////////////////////// + // sum over planes of Wilson loop with length R1 + // in the time direction + ////////////////////////////////////////////////// + static void siteTimelikeWilsonLoop(LatticeComplex &Wl, + const std::vector &U, + const int R1, const int R2) { + LatticeComplex siteWl(U[0]._grid); + + int ndim = U[0]._grid->_ndimension; + + Wl = zero; + for (int nu = 0; nu < ndim - 1; nu++) { + traceWilsonLoop(siteWl, U, R1, R2, ndim-1, nu); + Wl = Wl + siteWl; + } + } + ////////////////////////////////////////////////// + // sum Wilson loop over all planes orthogonal to the time direction + ////////////////////////////////////////////////// + static void siteSpatialWilsonLoop(LatticeComplex &Wl, + const std::vector &U, + const int R1, const int R2) { + LatticeComplex siteWl(U[0]._grid); + + Wl = zero; + for (int mu = 1; mu < U[0]._grid->_ndimension - 1; mu++) { for (int nu = 0; nu < mu; nu++) { traceWilsonLoop(siteWl, U, R1, R2, mu, nu); Wl = Wl + siteWl; @@ -139,7 +175,7 @@ public: const int R1, const int R2) { std::vector U(4, Umu._grid); - for (int mu = 0; mu < Nd; mu++) { + for (int mu = 0; mu < Umu._grid->_ndimension; mu++) { U[mu] = PeekIndex(Umu, mu); } @@ -152,14 +188,75 @@ public: return p.real(); } ////////////////////////////////////////////////// + // sum over all x,y,z,t and over all planes of timelike Wilson loop + ////////////////////////////////////////////////// + static Real sumTimelikeWilsonLoop(const GaugeLorentz &Umu, + const int R1, const int R2) { + std::vector U(4, Umu._grid); + + for (int mu = 0; mu < Umu._grid->_ndimension; mu++) { + U[mu] = PeekIndex(Umu, mu); + } + + LatticeComplex Wl(Umu._grid); + + siteTimelikeWilsonLoop(Wl, U, R1, R2); + + TComplex Tp = sum(Wl); + Complex p = TensorRemove(Tp); + return p.real(); + } + ////////////////////////////////////////////////// + // sum over all x,y,z,t and over all planes of spatial Wilson loop + ////////////////////////////////////////////////// + static Real sumSpatialWilsonLoop(const GaugeLorentz &Umu, + const int R1, const int R2) { + std::vector U(4, Umu._grid); + + for (int mu = 0; mu < Umu._grid->_ndimension; mu++) { + U[mu] = PeekIndex(Umu, mu); + } + + LatticeComplex Wl(Umu._grid); + + siteSpatialWilsonLoop(Wl, U, R1, R2); + + TComplex Tp = sum(Wl); + Complex p = TensorRemove(Tp); + return p.real(); + } + ////////////////////////////////////////////////// // average over all x,y,z,t and over all planes of Wilson loop ////////////////////////////////////////////////// static Real avgWilsonLoop(const GaugeLorentz &Umu, const int R1, const int R2) { + int ndim = Umu._grid->_ndimension; Real sumWl = sumWilsonLoop(Umu, R1, R2); - double vol = Umu._grid->gSites(); - double faces = 1.0 * Nd * (Nd - 1); - return sumWl / vol / faces / Nc; // Nd , Nc dependent... FIXME + Real vol = Umu._grid->gSites(); + Real faces = 1.0 * ndim * (ndim - 1); + return sumWl / vol / faces / Nc; // Nc dependent... FIXME + } + ////////////////////////////////////////////////// + // average over all x,y,z,t and over all planes of timelike Wilson loop + ////////////////////////////////////////////////// + static Real avgTimelikeWilsonLoop(const GaugeLorentz &Umu, + const int R1, const int R2) { + int ndim = Umu._grid->_ndimension; + Real sumWl = sumTimelikeWilsonLoop(Umu, R1, R2); + Real vol = Umu._grid->gSites(); + Real faces = 1.0 * (ndim - 1); + return sumWl / vol / faces / Nc; // Nc dependent... FIXME + } + ////////////////////////////////////////////////// + // average over all x,y,z,t and over all planes of spatial Wilson loop + ////////////////////////////////////////////////// + static Real avgSpatialWilsonLoop(const GaugeLorentz &Umu, + const int R1, const int R2) { + int ndim = Umu._grid->_ndimension; + Real sumWl = sumSpatialWilsonLoop(Umu, R1, R2); + Real vol = Umu._grid->gSites(); + Real faces = 1.0 * (ndim - 1) * (ndim - 2); + return sumWl / vol / faces / Nc; // Nc dependent... FIXME } }; diff --git a/programs/qed-fvol/qed-fvol.cc b/programs/qed-fvol/qed-fvol.cc index 68705b8f..d026057e 100644 --- a/programs/qed-fvol/qed-fvol.cc +++ b/programs/qed-fvol/qed-fvol.cc @@ -64,7 +64,7 @@ int main(int argc, char *argv[]) EmField a(&grid); EmField expA(&grid); - Real avgPlaqAexp, avgWl2x2Aexp; + Real avgPlaqAexp, avgWl2x2Aexp, avgWl2x2Aexp_time, avgWl2x2Aexp_space; pRNG.SeedRandomDevice(); photon.StochasticField(a, pRNG); @@ -117,14 +117,20 @@ int main(int argc, char *argv[]) // Calculate plaquette from exponentiated photon field avgPlaqAexp = NewWilsonLoops::avgPlaquette(expA); avgWl2x2Aexp = NewWilsonLoops::avgWilsonLoop(expA, 2, 2); + avgWl2x2Aexp_time = NewWilsonLoops::avgTimelikeWilsonLoop(expA, 2, 2); + avgWl2x2Aexp_space = NewWilsonLoops::avgSpatialWilsonLoop(expA, 2, 2); avgPlaqAexp = avgPlaqAexp*3; avgWl2x2Aexp = avgWl2x2Aexp*3; + avgWl2x2Aexp_time = avgWl2x2Aexp_time*3; + avgWl2x2Aexp_space = avgWl2x2Aexp_space*3; LOG(Message) << "Plaquette average (from A): " << avgPlaqA << std::endl; LOG(Message) << "Plaquette average (from exp(A)): " << avgPlaqAexp << std::endl; LOG(Message) << "2x2 Wilson Loop average (from A): " << avgWlA << std::endl; LOG(Message) << "2x2 Wilson Loop average (from exp(A)): " << avgWl2x2Aexp << std::endl; + LOG(Message) << "2x2 Wilson Loop timelike average (from exp(A)): " << avgWl2x2Aexp_time << std::endl; + LOG(Message) << "2x2 Wilson Loop spatial average (from exp(A)): " << avgWl2x2Aexp_space << std::endl; LOG(Message) << "Plaquette (one site): " << plaqsite / faces << std::endl; // epilogue From 92ec3404f8a404e7d6420ebfa0f113af5eb6ec6d Mon Sep 17 00:00:00 2001 From: James Harrison Date: Mon, 14 Nov 2016 17:59:02 +0000 Subject: [PATCH 08/53] Set imaginary part of stochastic QED field to zero after FFT into position space --- lib/qcd/action/gauge/Photon.h | 2 ++ 1 file changed, 2 insertions(+) diff --git a/lib/qcd/action/gauge/Photon.h b/lib/qcd/action/gauge/Photon.h index 852ecb3e..ca0a8d40 100644 --- a/lib/qcd/action/gauge/Photon.h +++ b/lib/qcd/action/gauge/Photon.h @@ -172,6 +172,8 @@ namespace QCD{ pokeLorentz(aTilde, r, mu); } fft.FFT_all_dim(out, aTilde, FFT::backward); + + out = 0.5*(out + conjugate(out)); } // template // void Photon::FeynmanGaugeMomentumSpacePropagator_L(GaugeField &out, From a71b69389b6fc7360dbebbf5ed8d4fa3a6952016 Mon Sep 17 00:00:00 2001 From: James Harrison Date: Mon, 14 Nov 2016 18:23:04 +0000 Subject: [PATCH 09/53] QedFVol: calculate square Wilson loops up to 10x10 --- programs/qed-fvol/qed-fvol.cc | 74 +++++++---------------------------- 1 file changed, 14 insertions(+), 60 deletions(-) diff --git a/programs/qed-fvol/qed-fvol.cc b/programs/qed-fvol/qed-fvol.cc index d026057e..31312b1e 100644 --- a/programs/qed-fvol/qed-fvol.cc +++ b/programs/qed-fvol/qed-fvol.cc @@ -64,75 +64,29 @@ int main(int argc, char *argv[]) EmField a(&grid); EmField expA(&grid); - Real avgPlaqAexp, avgWl2x2Aexp, avgWl2x2Aexp_time, avgWl2x2Aexp_space; + Real wlA, logWlA; pRNG.SeedRandomDevice(); photon.StochasticField(a, pRNG); // Exponentiate photon field Complex imag_unit(0, 1); - expA = exp(imag_unit*0.5*(a+conjugate(a))); + expA = exp(imag_unit*a); - // Calculate plaquette from photon field - EmComp plaqA(&grid); - EmComp wlA(&grid); - EmComp tmp(&grid); - std::vector a_comp(4, &grid); - - for (int dir = 0; dir < Nd; dir++) { - a_comp[dir] = PeekIndex(a, dir); + // Calculate Wilson loops + for(int i=1; i<=10; i++){ + LOG(Message) << i << 'x' << i << " Wilson loop" << std::endl; + wlA = NewWilsonLoops::avgWilsonLoop(expA, i, i) * 3; + logWlA = -2*log(wlA); + LOG(Message) << "-2log(W) average: " << logWlA << std::endl; + wlA = NewWilsonLoops::avgTimelikeWilsonLoop(expA, i, i) * 3; + logWlA = -2*log(wlA); + LOG(Message) << "-2log(W) timelike: " << logWlA << std::endl; + wlA = NewWilsonLoops::avgSpatialWilsonLoop(expA, i, i) * 3; + logWlA = -2*log(wlA); + LOG(Message) << "-2log(W) spatial: " << logWlA << std::endl; } - plaqA = zero; - wlA = zero; - - for(int mu = 1; mu < Nd; mu++) { - for(int nu = 0; nu < mu; nu++) { - tmp = a_comp[mu] + Cshift(a_comp[nu], mu, 1) - Cshift(a_comp[mu], nu, 1) - a_comp[nu]; - plaqA = plaqA + cos(tmp); - - tmp = a_comp[mu] + Cshift(a_comp[mu], mu, 1) - + Cshift(a_comp[nu], mu, 2) + Cshift(Cshift(a_comp[nu], mu, 2), nu, 1) - - Cshift(Cshift(a_comp[mu], nu, 2), mu, 1) - Cshift(a_comp[mu], nu, 2) - - Cshift(a_comp[nu], nu, 1) - a_comp[nu]; - wlA = wlA + cos(tmp); - } - } - - Real vol = grid.gSites(); - Real faces = (1.0 * Nd * (Nd - 1)) / 2.0; - - Complex avgPlaqA = sum(trace(plaqA)); - avgPlaqA = avgPlaqA / vol / faces; - - Complex avgWlA = sum(trace(wlA)); - avgWlA = avgWlA / vol / faces; - - TComplex tplaqsite; - LatticeComplex plaqtrace = trace(plaqA); - std::vector site0 = {0,0,0,0}; - peekSite(tplaqsite, plaqtrace, site0); - Complex plaqsite = TensorRemove(tplaqsite); - - // Calculate plaquette from exponentiated photon field - avgPlaqAexp = NewWilsonLoops::avgPlaquette(expA); - avgWl2x2Aexp = NewWilsonLoops::avgWilsonLoop(expA, 2, 2); - avgWl2x2Aexp_time = NewWilsonLoops::avgTimelikeWilsonLoop(expA, 2, 2); - avgWl2x2Aexp_space = NewWilsonLoops::avgSpatialWilsonLoop(expA, 2, 2); - - avgPlaqAexp = avgPlaqAexp*3; - avgWl2x2Aexp = avgWl2x2Aexp*3; - avgWl2x2Aexp_time = avgWl2x2Aexp_time*3; - avgWl2x2Aexp_space = avgWl2x2Aexp_space*3; - - LOG(Message) << "Plaquette average (from A): " << avgPlaqA << std::endl; - LOG(Message) << "Plaquette average (from exp(A)): " << avgPlaqAexp << std::endl; - LOG(Message) << "2x2 Wilson Loop average (from A): " << avgWlA << std::endl; - LOG(Message) << "2x2 Wilson Loop average (from exp(A)): " << avgWl2x2Aexp << std::endl; - LOG(Message) << "2x2 Wilson Loop timelike average (from exp(A)): " << avgWl2x2Aexp_time << std::endl; - LOG(Message) << "2x2 Wilson Loop spatial average (from exp(A)): " << avgWl2x2Aexp_space << std::endl; - LOG(Message) << "Plaquette (one site): " << plaqsite / faces << std::endl; - // epilogue LOG(Message) << "Grid is finalizing now" << std::endl; Grid_finalize(); From 739c2308b5ce9a9464dbbd9057dbe49f6b04cf59 Mon Sep 17 00:00:00 2001 From: James Harrison Date: Tue, 15 Nov 2016 13:07:52 +0000 Subject: [PATCH 10/53] Set imaginary part of stochastic QED field to zero using real() instead of conjugate(). --- lib/qcd/action/gauge/Photon.h | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/lib/qcd/action/gauge/Photon.h b/lib/qcd/action/gauge/Photon.h index ca0a8d40..b6c1b76f 100644 --- a/lib/qcd/action/gauge/Photon.h +++ b/lib/qcd/action/gauge/Photon.h @@ -173,7 +173,7 @@ namespace QCD{ } fft.FFT_all_dim(out, aTilde, FFT::backward); - out = 0.5*(out + conjugate(out)); + out = real(out); } // template // void Photon::FeynmanGaugeMomentumSpacePropagator_L(GaugeField &out, From 6ad73145bc9754a5f26093eee5a34473ba0cff82 Mon Sep 17 00:00:00 2001 From: James Harrison Date: Wed, 30 Nov 2016 15:17:22 +0000 Subject: [PATCH 11/53] Calculate Wilson loop average over multiple configurations. --- programs/qed-fvol/qed-fvol.cc | 47 +++++++++++++++++++++++------------ 1 file changed, 31 insertions(+), 16 deletions(-) diff --git a/programs/qed-fvol/qed-fvol.cc b/programs/qed-fvol/qed-fvol.cc index 31312b1e..f0f5079f 100644 --- a/programs/qed-fvol/qed-fvol.cc +++ b/programs/qed-fvol/qed-fvol.cc @@ -31,6 +31,9 @@ typedef Photon PhotonR; typedef PhotonR::GaugeField EmField; typedef PhotonR::GaugeLinkField EmComp; +const int NCONFIGS = 10; +const int NWILSON = 10; + int main(int argc, char *argv[]) { // parse command line @@ -64,27 +67,39 @@ int main(int argc, char *argv[]) EmField a(&grid); EmField expA(&grid); - Real wlA, logWlA; + Complex imag_unit(0, 1); + + Real wlA; + std::vector logWlAvg(NWILSON, 0.0), logWlTime(NWILSON, 0.0), logWlSpace(NWILSON, 0.0); pRNG.SeedRandomDevice(); - photon.StochasticField(a, pRNG); - // Exponentiate photon field - Complex imag_unit(0, 1); - expA = exp(imag_unit*a); + LOG(Message) << "Wilson loop calculation beginning" << std::endl; + for(int ic = 0; ic < NCONFIGS; ic++){ + LOG(Message) << "Configuration " << ic <::avgWilsonLoop(expA, iw, iw) * 3; + logWlAvg[iw-1] -= 2*log(wlA); + wlA = NewWilsonLoops::avgTimelikeWilsonLoop(expA, iw, iw) * 3; + logWlTime[iw-1] -= 2*log(wlA); + wlA = NewWilsonLoops::avgSpatialWilsonLoop(expA, iw, iw) * 3; + logWlSpace[iw-1] -= 2*log(wlA); + } + } + LOG(Message) << "Wilson loop calculation completed" << std::endl; + // Calculate Wilson loops - for(int i=1; i<=10; i++){ - LOG(Message) << i << 'x' << i << " Wilson loop" << std::endl; - wlA = NewWilsonLoops::avgWilsonLoop(expA, i, i) * 3; - logWlA = -2*log(wlA); - LOG(Message) << "-2log(W) average: " << logWlA << std::endl; - wlA = NewWilsonLoops::avgTimelikeWilsonLoop(expA, i, i) * 3; - logWlA = -2*log(wlA); - LOG(Message) << "-2log(W) timelike: " << logWlA << std::endl; - wlA = NewWilsonLoops::avgSpatialWilsonLoop(expA, i, i) * 3; - logWlA = -2*log(wlA); - LOG(Message) << "-2log(W) spatial: " << logWlA << std::endl; + for(int iw=1; iw<=10; iw++){ + LOG(Message) << iw << 'x' << iw << " Wilson loop" << std::endl; + LOG(Message) << "-2log(W) average: " << logWlAvg[iw-1]/NCONFIGS << std::endl; + LOG(Message) << "-2log(W) timelike: " << logWlTime[iw-1]/NCONFIGS << std::endl; + LOG(Message) << "-2log(W) spatial: " << logWlSpace[iw-1]/NCONFIGS << std::endl; } // epilogue From 2e3c5890b6035a4c9d661102c2117c53f93f00fd Mon Sep 17 00:00:00 2001 From: Antonin Portelli Date: Thu, 15 Dec 2016 20:06:46 +0000 Subject: [PATCH 12/53] qed-fvol: build fix --- extras/Makefile.am | 2 +- extras/qed-fvol/Makefile.am | 2 +- lib/qcd/action/Actions.h | 1 + 3 files changed, 3 insertions(+), 2 deletions(-) diff --git a/extras/Makefile.am b/extras/Makefile.am index d8c2b675..416a9fc8 100644 --- a/extras/Makefile.am +++ b/extras/Makefile.am @@ -1 +1 @@ -SUBDIRS = Hadrons \ No newline at end of file +SUBDIRS = Hadrons qed-fvol \ No newline at end of file diff --git a/extras/qed-fvol/Makefile.am b/extras/qed-fvol/Makefile.am index cd762e94..0a9030c7 100644 --- a/extras/qed-fvol/Makefile.am +++ b/extras/qed-fvol/Makefile.am @@ -1,4 +1,4 @@ -AM_CXXFLAGS += -I$(top_srcdir)/programs -I../$(top_srcdir)/programs +AM_CXXFLAGS += -I$(top_srcdir)/extras bin_PROGRAMS = qed-fvol diff --git a/lib/qcd/action/Actions.h b/lib/qcd/action/Actions.h index 4a30f8c3..fea75f8a 100644 --- a/lib/qcd/action/Actions.h +++ b/lib/qcd/action/Actions.h @@ -57,6 +57,7 @@ Author: paboyle //////////////////////////////////////////// // Gauge Actions //////////////////////////////////////////// +#include #include #include From 2af9ab903445291377bb323ee349ddf9c7e94abf Mon Sep 17 00:00:00 2001 From: Antonin Portelli Date: Tue, 20 Dec 2016 12:40:26 +0100 Subject: [PATCH 13/53] old Makefile cleaning --- programs/Makefile.am | 1 - 1 file changed, 1 deletion(-) delete mode 100644 programs/Makefile.am diff --git a/programs/Makefile.am b/programs/Makefile.am deleted file mode 100644 index ff7f6584..00000000 --- a/programs/Makefile.am +++ /dev/null @@ -1 +0,0 @@ -SUBDIRS = qed-fvol From 9ac3ac41df095e3208c126f4b52bdf9f1b58937a Mon Sep 17 00:00:00 2001 From: Antonin Portelli Date: Tue, 20 Dec 2016 12:41:01 +0100 Subject: [PATCH 14/53] serialisable Photon parameters --- lib/qcd/action/gauge/Photon.h | 9 +++++---- 1 file changed, 5 insertions(+), 4 deletions(-) diff --git a/lib/qcd/action/gauge/Photon.h b/lib/qcd/action/gauge/Photon.h index b6c1b76f..bbe3ebf7 100644 --- a/lib/qcd/action/gauge/Photon.h +++ b/lib/qcd/action/gauge/Photon.h @@ -28,6 +28,7 @@ #ifndef QCD_PHOTON_ACTION_H #define QCD_PHOTON_ACTION_H + namespace Grid{ namespace QCD{ @@ -36,8 +37,8 @@ namespace QCD{ { public: INHERIT_GIMPL_TYPES(Gimpl); - enum class Gauge {Feynman, Coulomb, Landau}; - enum class ZmScheme {QedL, QedTL}; + GRID_SERIALIZABLE_ENUM(Gauge, undef, feynman, 1, coulomb, 2, landau, 3); + GRID_SERIALIZABLE_ENUM(ZmScheme, undef, qedL, 1, qedTL, 2); public: Photon(Gauge gauge, ZmScheme zmScheme); virtual ~Photon(void) = default; @@ -104,7 +105,7 @@ namespace QCD{ switch (zmScheme_) { - case ZmScheme::QedTL: + case ZmScheme::qedTL: { std::vector zm(nd,0); TComplex Tzero = Complex(0.0,0.0); @@ -113,7 +114,7 @@ namespace QCD{ break; } - case ZmScheme::QedL: + case ZmScheme::qedL: { LatticeInteger spNrm(grid), coor(grid); GaugeLinkField z(grid); From db9c28a773c5d93d3c757ea4f8af75876106b948 Mon Sep 17 00:00:00 2001 From: Antonin Portelli Date: Tue, 20 Dec 2016 12:41:39 +0100 Subject: [PATCH 15/53] qed-fvol: Photon parameter name fix --- extras/qed-fvol/qed-fvol.cc | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/extras/qed-fvol/qed-fvol.cc b/extras/qed-fvol/qed-fvol.cc index f0f5079f..951c36ad 100644 --- a/extras/qed-fvol/qed-fvol.cc +++ b/extras/qed-fvol/qed-fvol.cc @@ -62,8 +62,8 @@ int main(int argc, char *argv[]) std::vector mpi_layout = GridDefaultMpi(); GridCartesian grid(latt_size,simd_layout,mpi_layout); GridParallelRNG pRNG(&grid); - PhotonR photon(PhotonR::Gauge::Feynman, - PhotonR::ZmScheme::QedL); + PhotonR photon(PhotonR::Gauge::feynman, + PhotonR::ZmScheme::qedL); EmField a(&grid); EmField expA(&grid); From 17b3a10d46e46823f0a380647c4953c2ffd74ea4 Mon Sep 17 00:00:00 2001 From: Antonin Portelli Date: Thu, 22 Dec 2016 00:29:19 +0100 Subject: [PATCH 16/53] stochastic QED: function to cache 1/sqrt(khat^2) --- lib/qcd/action/gauge/Photon.h | 47 +++++++++++++++++++++++++---------- 1 file changed, 34 insertions(+), 13 deletions(-) diff --git a/lib/qcd/action/gauge/Photon.h b/lib/qcd/action/gauge/Photon.h index bbe3ebf7..faa63b42 100644 --- a/lib/qcd/action/gauge/Photon.h +++ b/lib/qcd/action/gauge/Photon.h @@ -44,7 +44,10 @@ namespace QCD{ virtual ~Photon(void) = default; void FreePropagator(const GaugeField &in, GaugeField &out); void MomentumSpacePropagator(const GaugeField &in, GaugeField &out); + void StochasticWeight(GaugeLinkField &weight); void StochasticField(GaugeField &out, GridParallelRNG &rng); + void StochasticField(GaugeField &out, GridParallelRNG &rng, + const GaugeLinkField &weight); private: void invKHatSquared(GaugeLinkField &out); void zmSub(GaugeLinkField &out); @@ -148,32 +151,50 @@ namespace QCD{ } template - void Photon::StochasticField(GaugeField &out, GridParallelRNG &rng) + void Photon::StochasticWeight(GaugeLinkField &weight) { - auto *grid = dynamic_cast(out._grid); - const unsigned int nd = grid->_ndimension; - std::vector latt_size = grid->_fdimensions; - GaugeLinkField sqrtK2Inv(grid), r(grid); - GaugeField aTilde(grid); - FFT fft(grid); + auto *grid = dynamic_cast(weight._grid); + const unsigned int nd = grid->_ndimension; + std::vector latt_size = grid->_fdimensions; Integer vol = 1; for(int d = 0; d < nd; d++) { vol = vol * latt_size[d]; } - - invKHatSquared(sqrtK2Inv); - sqrtK2Inv = sqrt(vol*real(sqrtK2Inv)); - zmSub(sqrtK2Inv); + invKHatSquared(weight); + weight = sqrt(vol*real(weight)); + zmSub(weight); + } + + template + void Photon::StochasticField(GaugeField &out, GridParallelRNG &rng) + { + auto *grid = dynamic_cast(out._grid); + GaugeLinkField weight(grid); + + StochasticWeight(weight); + StochasticField(out, rng, weight); + } + + template + void Photon::StochasticField(GaugeField &out, GridParallelRNG &rng, + const GaugeLinkField &weight) + { + auto *grid = dynamic_cast(out._grid); + const unsigned int nd = grid->_ndimension; + GaugeLinkField r(grid); + GaugeField aTilde(grid); + FFT fft(grid); + for(int mu = 0; mu < nd; mu++) { gaussian(rng, r); - r = sqrtK2Inv*r; + r = weight*r; pokeLorentz(aTilde, r, mu); } fft.FFT_all_dim(out, aTilde, FFT::backward); - + out = real(out); } // template From 4c3fd9fa3f6976c1297715d1e2239797bb0dd45b Mon Sep 17 00:00:00 2001 From: Antonin Portelli Date: Thu, 22 Dec 2016 00:29:41 +0100 Subject: [PATCH 17/53] stochastic QED field module in Hadrons --- extras/Hadrons/Modules.hpp | 1 + extras/Hadrons/Modules/MGauge/StochEm.cc | 88 +++++++++++++++++++++ extras/Hadrons/Modules/MGauge/StochEm.hpp | 96 +++++++++++++++++++++++ extras/Hadrons/modules.inc | 2 + 4 files changed, 187 insertions(+) create mode 100644 extras/Hadrons/Modules/MGauge/StochEm.cc create mode 100644 extras/Hadrons/Modules/MGauge/StochEm.hpp diff --git a/extras/Hadrons/Modules.hpp b/extras/Hadrons/Modules.hpp index 77ae08b7..5d1a456c 100644 --- a/extras/Hadrons/Modules.hpp +++ b/extras/Hadrons/Modules.hpp @@ -32,6 +32,7 @@ See the full license in the file "LICENSE" in the top level distribution directo #include #include #include +#include #include #include #include diff --git a/extras/Hadrons/Modules/MGauge/StochEm.cc b/extras/Hadrons/Modules/MGauge/StochEm.cc new file mode 100644 index 00000000..c7a9fc4f --- /dev/null +++ b/extras/Hadrons/Modules/MGauge/StochEm.cc @@ -0,0 +1,88 @@ +/************************************************************************************* + +Grid physics library, www.github.com/paboyle/Grid + +Source file: extras/Hadrons/Modules/MGauge/StochEm.cc + +Copyright (C) 2015 +Copyright (C) 2016 + + +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 + +using namespace Grid; +using namespace Hadrons; +using namespace MGauge; + +/****************************************************************************** +* TStochEm implementation * +******************************************************************************/ +// constructor ///////////////////////////////////////////////////////////////// +TStochEm::TStochEm(const std::string name) +: Module(name) +{} + +// dependencies/products /////////////////////////////////////////////////////// +std::vector TStochEm::getInput(void) +{ + std::vector in; + + return in; +} + +std::vector TStochEm::getOutput(void) +{ + std::vector out = {getName()}; + + return out; +} + +// setup /////////////////////////////////////////////////////////////////////// +void TStochEm::setup(void) +{ + if (!env().hasRegisteredObject("_" + getName() + "_weight")) + { + env().registerLattice("_" + getName() + "_weight"); + } + env().registerLattice(getName()); +} + +// execution /////////////////////////////////////////////////////////////////// +void TStochEm::execute(void) +{ + PhotonR photon(par().gauge, par().zmScheme); + EmField &a = *env().createLattice(getName()); + EmComp *w; + + 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("_" + getName() + "_weight"); + photon.StochasticWeight(*w); + } + else + { + w = env().getObject("_" + getName() + "_weight"); + } + LOG(Message) << "Generating stochatic EM potential..." << std::endl; + photon.StochasticField(a, *env().get4dRng(), *w); +} diff --git a/extras/Hadrons/Modules/MGauge/StochEm.hpp b/extras/Hadrons/Modules/MGauge/StochEm.hpp new file mode 100644 index 00000000..04a7c48c --- /dev/null +++ b/extras/Hadrons/Modules/MGauge/StochEm.hpp @@ -0,0 +1,96 @@ +/************************************************************************************* + +Grid physics library, www.github.com/paboyle/Grid + +Source file: extras/Hadrons/Modules/MGauge/StochEm.hpp + +Copyright (C) 2015 +Copyright (C) 2016 + + +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_StochEm_hpp_ +#define Hadrons_StochEm_hpp_ + +#include +#include +#include + +BEGIN_HADRONS_NAMESPACE + +/****************************************************************************** + * StochEm * + ******************************************************************************/ +BEGIN_MODULE_NAMESPACE(MGauge) + +template +class QedGimpl +{ +public: + typedef S Simd; + + template + using iImplGaugeLink = iScalar>>; + template + using iImplGaugeField = iVector>, Nd>; + + typedef iImplGaugeLink SiteGaugeLink; + typedef iImplGaugeField SiteGaugeField; + + typedef Lattice GaugeLinkField; + typedef Lattice GaugeField; +}; + +typedef QedGimpl QedGimplR; +typedef Photon PhotonR; + +class StochEmPar: Serializable +{ +public: + GRID_SERIALIZABLE_CLASS_MEMBERS(StochEmPar, + PhotonR::Gauge, gauge, + PhotonR::ZmScheme, zmScheme); +}; + +class TStochEm: public Module +{ +public: + typedef PhotonR::GaugeField EmField; + typedef PhotonR::GaugeLinkField EmComp; +public: + // constructor + TStochEm(const std::string name); + // destructor + virtual ~TStochEm(void) = default; + // dependency relation + virtual std::vector getInput(void); + virtual std::vector getOutput(void); + // setup + virtual void setup(void); + // execution + virtual void execute(void); +}; + +MODULE_REGISTER_NS(StochEm, TStochEm, MGauge); + +END_MODULE_NAMESPACE + +END_HADRONS_NAMESPACE + +#endif // Hadrons_StochEm_hpp_ diff --git a/extras/Hadrons/modules.inc b/extras/Hadrons/modules.inc index 4251ffa3..8b559024 100644 --- a/extras/Hadrons/modules.inc +++ b/extras/Hadrons/modules.inc @@ -1,6 +1,7 @@ modules_cc =\ Modules/MGauge/Load.cc \ Modules/MGauge/Random.cc \ + Modules/MGauge/StochEm.cc \ Modules/MGauge/Unit.cc modules_hpp =\ @@ -10,6 +11,7 @@ modules_hpp =\ Modules/MContraction/Meson.hpp \ Modules/MGauge/Load.hpp \ Modules/MGauge/Random.hpp \ + Modules/MGauge/StochEm.hpp \ Modules/MGauge/Unit.hpp \ Modules/MSolver/RBPrecCG.hpp \ Modules/MSource/Point.hpp \ From 8c3cc3236447b4b8eef95a29da1b48166b5eb03d Mon Sep 17 00:00:00 2001 From: Antonin Portelli Date: Thu, 29 Dec 2016 22:42:58 +0100 Subject: [PATCH 18/53] Scalar action --- lib/qcd/action/Actions.h | 5 + lib/qcd/action/scalar/Scalar.h | 211 +++++++++++++++++++++++++++++++++ 2 files changed, 216 insertions(+) create mode 100644 lib/qcd/action/scalar/Scalar.h diff --git a/lib/qcd/action/Actions.h b/lib/qcd/action/Actions.h index fea75f8a..efd6a5bc 100644 --- a/lib/qcd/action/Actions.h +++ b/lib/qcd/action/Actions.h @@ -292,4 +292,9 @@ typedef MobiusFermion GparityMobiusFermionD; #include #include +//////////////////// +// Scalar actions +//////////////////// +#include + #endif diff --git a/lib/qcd/action/scalar/Scalar.h b/lib/qcd/action/scalar/Scalar.h new file mode 100644 index 00000000..194f6767 --- /dev/null +++ b/lib/qcd/action/scalar/Scalar.h @@ -0,0 +1,211 @@ +#ifndef QCD_SCALAR_ACTION_H +#define QCD_SCALAR_ACTION_H + +#define INHERIT_SIMPL_TYPES(Impl)\ +typedef typename Impl::SiteScalar SiteScalar; \ +typedef typename Impl::SiteSpinor SiteSpinor; \ +typedef typename Impl::SitePropagator SitePropagator; \ +typedef typename Impl::ScalarField ScalarField; \ +typedef typename Impl::FermionField FermionField; \ +typedef typename Impl::PropagatorField PropagatorField; \ +typedef typename Impl::StencilImpl StencilImpl; + +namespace Grid{ +namespace QCD{ + // Scalar implementation class /////////////////////////////////////////////// + // FIXME: it is not very nice to have the FImpl aliases + template , + class _Coeff_t = RealD> + class ScalarImpl: + public PeriodicGaugeImpl> + { + public: + static constexpr unsigned int rDim = Representation::Dimension; + public: + // gauge types + typedef PeriodicGaugeImpl> Gimpl; + INHERIT_GIMPL_TYPES(Gimpl); + // site types + // (using classes instead of aliases to allow for partial specialisation) + template + class iImplScalar + { + public: + typedef iScalar>> type; + }; + template + class iImplScalar + { + public: + typedef iScalar>> type; + }; + template + class iImplPropagator + { + public: + typedef iScalar>> type; + }; + template + class iImplPropagator + { + public: + typedef iScalar>> type; + }; + // type aliases + typedef typename iImplScalar::type SiteScalar; + typedef SiteScalar SiteSpinor; + typedef typename iImplPropagator::type SitePropagator; + typedef Lattice ScalarField; + typedef ScalarField FermionField; + typedef Lattice PropagatorField; + typedef CartesianStencil StencilImpl; + }; + + // single scalar implementation + typedef ScalarImpl ScalarImplR; + + // Scalar action ///////////////////////////////////////////////////////////// + template + class Scalar: + public CheckerBoardedSparseMatrixBase, + public SImpl + { + public: + INHERIT_GIMPL_TYPES(SImpl); + INHERIT_SIMPL_TYPES(SImpl); + public: + // constructor + Scalar(GaugeField &_Umu, GridCartesian &Sgrid, GridRedBlackCartesian &Hgrid, + RealD _mass) + : _grid(&Sgrid) + , _cbgrid(&Hgrid) + , mass(_mass) + , Lebesgue(_grid) + , LebesgueEvenOdd(_cbgrid) + , Umu(&Sgrid) + , UmuEven(&Hgrid) + , UmuOdd(&Hgrid) + { + Umu = _Umu; + pickCheckerboard(Even, UmuEven, Umu); + pickCheckerboard(Odd, UmuOdd, Umu); + } + // grid access + virtual GridBase *RedBlackGrid(void) {return _grid;} + // half checkerboard operations + // FIXME: do implementation + virtual void Meooe(const ScalarField &in, ScalarField &out) + { + assert(0); + } + virtual void Mooee(const ScalarField &in, ScalarField &out) + { + assert(0); + } + virtual void MooeeInv(const ScalarField &in, ScalarField &out) + { + assert(0); + } + virtual void MeooeDag(const ScalarField &in, ScalarField &out) + { + assert(0); + } + virtual void MooeeDag(const ScalarField &in, ScalarField &out) + { + assert(0); + } + virtual void MooeeInvDag(const ScalarField &in, ScalarField &out) + { + assert(0); + } + // free propagators + static void MomentumSpacePropagator(ScalarField &out, RealD m); + static void FreePropagator(const ScalarField &in, ScalarField &out, + const ScalarField &momKernel); + static void FreePropagator(const ScalarField &in, ScalarField &out, RealD m); + public: + RealD mass; + + GridBase *_grid; + GridBase *_cbgrid; + + // Defines the stencils for even and odd + StencilImpl Stencil; + StencilImpl StencilEven; + StencilImpl StencilOdd; + + // Copy of the gauge field, with even and odd subsets + GaugeField Umu; + GaugeField UmuEven; + GaugeField UmuOdd; + + LebesgueOrder Lebesgue; + LebesgueOrder LebesgueEvenOdd; + }; + + template + void Scalar::MomentumSpacePropagator(ScalarField &out, RealD m) + { + GridBase *grid = out._grid; + ScalarField kmu(grid); + const unsigned int nd = grid->_ndimension; + std::vector &l = grid->_fdimensions; + + out = m*m; + for(int mu = 0; mu < nd; mu++) + { + Real twoPiL = M_PI*2./l[mu]; + + LatticeCoordinate(kmu,mu); + kmu = 2.*sin(.5*twoPiL*kmu); + out = out + kmu*kmu; + } + } + + template + void Scalar::FreePropagator(const ScalarField &in, ScalarField &out, + const ScalarField &FTKernel) + { + FFT fft((GridCartesian *)in._grid); + ScalarField inFT(in._grid); + + fft.FFT_all_dim(inFT, in, FFT::forward); + inFT = inFT*FTKernel; + fft.FFT_all_dim(out, inFT, FFT::backward); + } + + template + void Scalar::FreePropagator(const ScalarField &in, ScalarField &out, + RealD m) + { + ScalarField FTKernel(in._grid); + + MomentumSpacePropagator(FTKernel, m); + FreePropagator(in, out, FTKernel); + } + + template + void ScalarToProp(typename SImpl::PropagatorField &p, + const typename SImpl::ScalarField &s, + const int c) + { + for(int i = 0; i < SImpl::rDim; ++i) + { + pokeColour(p, peekColour(s, i), i); + } + } + + template + void PropToScalar(typename SImpl::ScalarField &s, + const typename SImpl::PropagatorField &p, + const int c) + { + for(int i = 0; i < SImpl::rDim; ++i) + { + pokeColour(s, peekColour(p, i), i); + } + } +}} + +#endif From afbf7d4c37df8f134aa6bb191d1fd29d1709b16e Mon Sep 17 00:00:00 2001 From: Antonin Portelli Date: Thu, 29 Dec 2016 22:43:38 +0100 Subject: [PATCH 19/53] QED Gimpl moved in Photon.h --- extras/Hadrons/Modules/MGauge/StochEm.hpp | 21 --------------------- extras/qed-fvol/qed-fvol.cc | 22 ---------------------- lib/qcd/action/gauge/Photon.h | 22 +++++++++++++++++++++- 3 files changed, 21 insertions(+), 44 deletions(-) diff --git a/extras/Hadrons/Modules/MGauge/StochEm.hpp b/extras/Hadrons/Modules/MGauge/StochEm.hpp index 04a7c48c..50a77435 100644 --- a/extras/Hadrons/Modules/MGauge/StochEm.hpp +++ b/extras/Hadrons/Modules/MGauge/StochEm.hpp @@ -39,27 +39,6 @@ BEGIN_HADRONS_NAMESPACE ******************************************************************************/ BEGIN_MODULE_NAMESPACE(MGauge) -template -class QedGimpl -{ -public: - typedef S Simd; - - template - using iImplGaugeLink = iScalar>>; - template - using iImplGaugeField = iVector>, Nd>; - - typedef iImplGaugeLink SiteGaugeLink; - typedef iImplGaugeField SiteGaugeField; - - typedef Lattice GaugeLinkField; - typedef Lattice GaugeField; -}; - -typedef QedGimpl QedGimplR; -typedef Photon PhotonR; - class StochEmPar: Serializable { public: diff --git a/extras/qed-fvol/qed-fvol.cc b/extras/qed-fvol/qed-fvol.cc index 951c36ad..3ecac2fc 100644 --- a/extras/qed-fvol/qed-fvol.cc +++ b/extras/qed-fvol/qed-fvol.cc @@ -5,29 +5,7 @@ using namespace Grid; using namespace QCD; using namespace QedFVol; -template -class QedGimpl -{ -public: - typedef S Simd; - - template - using iImplGaugeLink = iScalar>>; - template - using iImplGaugeField = iVector>, Nd>; - - typedef iImplGaugeLink SiteGaugeLink; - typedef iImplGaugeField SiteGaugeField; - - typedef Lattice GaugeLinkField; // bit ugly naming; polarised - // gauge field, lorentz... all - // ugly - typedef Lattice GaugeField; -}; - -typedef QedGimpl QedGimplR; typedef PeriodicGaugeImpl QedPeriodicGimplR; -typedef Photon PhotonR; typedef PhotonR::GaugeField EmField; typedef PhotonR::GaugeLinkField EmComp; diff --git a/lib/qcd/action/gauge/Photon.h b/lib/qcd/action/gauge/Photon.h index faa63b42..73405297 100644 --- a/lib/qcd/action/gauge/Photon.h +++ b/lib/qcd/action/gauge/Photon.h @@ -28,9 +28,27 @@ #ifndef QCD_PHOTON_ACTION_H #define QCD_PHOTON_ACTION_H - namespace Grid{ namespace QCD{ + template + class QedGimpl + { + public: + typedef S Simd; + + template + using iImplGaugeLink = iScalar>>; + template + using iImplGaugeField = iVector>, Nd>; + + typedef iImplGaugeLink SiteGaugeLink; + typedef iImplGaugeField SiteGaugeField; + + typedef Lattice GaugeLinkField; + typedef Lattice GaugeField; + }; + + typedef QedGimpl QedGimplR; template class Photon @@ -56,6 +74,8 @@ namespace QCD{ ZmScheme zmScheme_; }; + typedef Photon PhotonR; + template Photon::Photon(Gauge gauge, ZmScheme zmScheme) : gauge_(gauge), zmScheme_(zmScheme) From 4c60e31070f3d05ba7e52c15c4bf6e59644a046a Mon Sep 17 00:00:00 2001 From: Antonin Portelli Date: Thu, 29 Dec 2016 22:44:08 +0100 Subject: [PATCH 20/53] Hadrons: code cleaning --- extras/Hadrons/Modules/Quark.hpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/extras/Hadrons/Modules/Quark.hpp b/extras/Hadrons/Modules/Quark.hpp index e441a096..0cf7314b 100644 --- a/extras/Hadrons/Modules/Quark.hpp +++ b/extras/Hadrons/Modules/Quark.hpp @@ -133,7 +133,7 @@ void TQuark::execute(void) for (unsigned int c = 0; c < Nc; ++c) { LOG(Message) << "Inversion for spin= " << s << ", color= " << c - << std::endl; + << std::endl; // source conversion for 4D sources if (!env().isObject5d(par().source)) { From bbc0eff078cfd331ca31f2dd0c95b3030ee8a261 Mon Sep 17 00:00:00 2001 From: Antonin Portelli Date: Thu, 29 Dec 2016 22:44:22 +0100 Subject: [PATCH 21/53] Hadrons: scalar sources --- extras/Hadrons/Modules/MSource/Point.hpp | 5 +++-- extras/Hadrons/Modules/MSource/Z2.hpp | 5 +++-- 2 files changed, 6 insertions(+), 4 deletions(-) diff --git a/extras/Hadrons/Modules/MSource/Point.hpp b/extras/Hadrons/Modules/MSource/Point.hpp index a0ecbc2a..8d0b4de8 100644 --- a/extras/Hadrons/Modules/MSource/Point.hpp +++ b/extras/Hadrons/Modules/MSource/Point.hpp @@ -63,7 +63,7 @@ template class TPoint: public Module { public: - TYPE_ALIASES(FImpl,); + FERM_TYPE_ALIASES(FImpl,); public: // constructor TPoint(const std::string name); @@ -78,7 +78,8 @@ public: virtual void execute(void); }; -MODULE_REGISTER_NS(Point, TPoint, MSource); +MODULE_REGISTER_NS(Point, TPoint, MSource); +MODULE_REGISTER_NS(ScalarPoint, TPoint, MSource); /****************************************************************************** * TPoint template implementation * diff --git a/extras/Hadrons/Modules/MSource/Z2.hpp b/extras/Hadrons/Modules/MSource/Z2.hpp index cd5727be..6fa49cfe 100644 --- a/extras/Hadrons/Modules/MSource/Z2.hpp +++ b/extras/Hadrons/Modules/MSource/Z2.hpp @@ -67,7 +67,7 @@ template class TZ2: public Module { public: - TYPE_ALIASES(FImpl,); + FERM_TYPE_ALIASES(FImpl,); public: // constructor TZ2(const std::string name); @@ -82,7 +82,8 @@ public: virtual void execute(void); }; -MODULE_REGISTER_NS(Z2, TZ2, MSource); +MODULE_REGISTER_NS(Z2, TZ2, MSource); +MODULE_REGISTER_NS(ScalarZ2, TZ2, MSource); /****************************************************************************** * TZ2 template implementation * From 673994b281e6c464b4021c62c80a9976e0035176 Mon Sep 17 00:00:00 2001 From: Antonin Portelli Date: Thu, 29 Dec 2016 22:44:58 +0100 Subject: [PATCH 22/53] Hadrons: modules for scalar propagators --- extras/Hadrons/Global.hpp | 25 ++++++-- extras/Hadrons/Modules.hpp | 30 +--------- extras/Hadrons/Modules/MScalar/ChargedProp.cc | 40 +++++++++++++ .../Hadrons/Modules/MScalar/ChargedProp.hpp | 44 ++++++++++++++ extras/Hadrons/Modules/MScalar/FreeProp.cc | 57 +++++++++++++++++++ extras/Hadrons/Modules/MScalar/FreeProp.hpp | 47 +++++++++++++++ extras/Hadrons/modules.inc | 6 +- 7 files changed, 215 insertions(+), 34 deletions(-) create mode 100644 extras/Hadrons/Modules/MScalar/ChargedProp.cc create mode 100644 extras/Hadrons/Modules/MScalar/ChargedProp.hpp create mode 100644 extras/Hadrons/Modules/MScalar/FreeProp.cc create mode 100644 extras/Hadrons/Modules/MScalar/FreeProp.hpp diff --git a/extras/Hadrons/Global.hpp b/extras/Hadrons/Global.hpp index 81afab13..bcb282fc 100644 --- a/extras/Hadrons/Global.hpp +++ b/extras/Hadrons/Global.hpp @@ -51,23 +51,38 @@ using Grid::operator<<; * error with GCC 5 (clang & GCC 6 compile fine without it). */ -// FIXME: find a way to do that in a more general fashion #ifndef FIMPL #define FIMPL WilsonImplR #endif +#ifndef SIMPL +#define SIMPL ScalarImplR +#endif BEGIN_HADRONS_NAMESPACE // type aliases -#define TYPE_ALIASES(FImpl, suffix)\ +#define FERM_TYPE_ALIASES(FImpl, suffix)\ typedef FermionOperator FMat##suffix; \ typedef typename FImpl::FermionField FermionField##suffix; \ typedef typename FImpl::PropagatorField PropagatorField##suffix; \ -typedef typename FImpl::SitePropagator SitePropagator##suffix; \ -typedef typename FImpl::DoubledGaugeField DoubledGaugeField##suffix;\ -typedef std::function SolverFn##suffix; +#define TYPE_ALIASES(FImpl, suffix)\ +FERM_TYPE_ALIASES(FImpl, suffix)\ +GAUGE_TYPE_ALIASES(FImpl, suffix)\ +SOLVER_TYPE_ALIASES(FImpl, suffix) + // logger class HadronsLogger: public Logger { diff --git a/extras/Hadrons/Modules.hpp b/extras/Hadrons/Modules.hpp index 5d1a456c..ad31d2a7 100644 --- a/extras/Hadrons/Modules.hpp +++ b/extras/Hadrons/Modules.hpp @@ -1,31 +1,3 @@ -/************************************************************************************* - -Grid physics library, www.github.com/paboyle/Grid - -Source file: extras/Hadrons/Modules.hpp - -Copyright (C) 2015 -Copyright (C) 2016 - -Author: Antonin Portelli - -This program is free software; you can redistribute it and/or modify -it under the terms of the GNU General Public License as published by -the Free Software Foundation; either version 2 of the License, or -(at your option) any later version. - -This program is distributed in the hope that it will be useful, -but WITHOUT ANY WARRANTY; without even the implied warranty of -MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the -GNU General Public License for more details. - -You should have received a copy of the GNU General Public License along -with this program; if not, write to the Free Software Foundation, Inc., -51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. - -See the full license in the file "LICENSE" in the top level distribution directory -*************************************************************************************/ -/* END LEGAL */ #include #include #include @@ -34,6 +6,8 @@ See the full license in the file "LICENSE" in the top level distribution directo #include #include #include +#include +#include #include #include #include diff --git a/extras/Hadrons/Modules/MScalar/ChargedProp.cc b/extras/Hadrons/Modules/MScalar/ChargedProp.cc new file mode 100644 index 00000000..1137c6f0 --- /dev/null +++ b/extras/Hadrons/Modules/MScalar/ChargedProp.cc @@ -0,0 +1,40 @@ +#include + +using namespace Grid; +using namespace Hadrons; +using namespace MScalar; + +/****************************************************************************** +* TChargedProp implementation * +******************************************************************************/ +// constructor ///////////////////////////////////////////////////////////////// +TChargedProp::TChargedProp(const std::string name) +: Module(name) +{} + +// dependencies/products /////////////////////////////////////////////////////// +std::vector TChargedProp::getInput(void) +{ + std::vector in; + + return in; +} + +std::vector TChargedProp::getOutput(void) +{ + std::vector out = {getName()}; + + return out; +} + +// setup /////////////////////////////////////////////////////////////////////// +void TChargedProp::setup(void) +{ + +} + +// execution /////////////////////////////////////////////////////////////////// +void TChargedProp::execute(void) +{ + +} diff --git a/extras/Hadrons/Modules/MScalar/ChargedProp.hpp b/extras/Hadrons/Modules/MScalar/ChargedProp.hpp new file mode 100644 index 00000000..7a60c2ad --- /dev/null +++ b/extras/Hadrons/Modules/MScalar/ChargedProp.hpp @@ -0,0 +1,44 @@ +#ifndef Hadrons_ChargedProp_hpp_ +#define Hadrons_ChargedProp_hpp_ + +#include +#include +#include + +BEGIN_HADRONS_NAMESPACE + +/****************************************************************************** + * ChargedProp * + ******************************************************************************/ +BEGIN_MODULE_NAMESPACE(MScalar) + +class ChargedPropPar: Serializable +{ +public: + GRID_SERIALIZABLE_CLASS_MEMBERS(ChargedPropPar, + unsigned int, i); +}; + +class TChargedProp: public Module +{ +public: + // constructor + TChargedProp(const std::string name); + // destructor + virtual ~TChargedProp(void) = default; + // dependency relation + virtual std::vector getInput(void); + virtual std::vector getOutput(void); + // setup + virtual void setup(void); + // execution + virtual void execute(void); +}; + +MODULE_REGISTER_NS(ChargedProp, TChargedProp, MScalar); + +END_MODULE_NAMESPACE + +END_HADRONS_NAMESPACE + +#endif // Hadrons_ChargedProp_hpp_ diff --git a/extras/Hadrons/Modules/MScalar/FreeProp.cc b/extras/Hadrons/Modules/MScalar/FreeProp.cc new file mode 100644 index 00000000..7419a954 --- /dev/null +++ b/extras/Hadrons/Modules/MScalar/FreeProp.cc @@ -0,0 +1,57 @@ +#include + +using namespace Grid; +using namespace Hadrons; +using namespace MScalar; + +/****************************************************************************** +* TFreeProp implementation * +******************************************************************************/ +// constructor ///////////////////////////////////////////////////////////////// +TFreeProp::TFreeProp(const std::string name) +: Module(name) +{} + +// dependencies/products /////////////////////////////////////////////////////// +std::vector TFreeProp::getInput(void) +{ + std::vector in = {par().source}; + + return in; +} + +std::vector TFreeProp::getOutput(void) +{ + std::vector out = {getName()}; + + return out; +} + +// setup /////////////////////////////////////////////////////////////////////// +void TFreeProp::setup(void) +{ + env().registerLattice(getName()); +} + +// execution /////////////////////////////////////////////////////////////////// +void TFreeProp::execute(void) +{ + ScalarField &prop = *env().createLattice(getName()); + ScalarField &source = *env().getObject(par().source); + ScalarField *momKernel; + std::string kerName = "_" + getName() + "_momKernel"; + + if (!env().hasCreatedObject(kerName)) + { + LOG(Message) << "Caching momentum space free scalar propagator" + << "(mass= " << par().mass << ")..." << std::endl; + momKernel = env().template createLattice(kerName); + Scalar::MomentumSpacePropagator(*momKernel, par().mass); + } + else + { + momKernel = env().getObject(kerName); + } + LOG(Message) << "Computing free scalar propagator..." << std::endl; + Scalar::FreePropagator(source, prop, *momKernel); +} diff --git a/extras/Hadrons/Modules/MScalar/FreeProp.hpp b/extras/Hadrons/Modules/MScalar/FreeProp.hpp new file mode 100644 index 00000000..6a0cd930 --- /dev/null +++ b/extras/Hadrons/Modules/MScalar/FreeProp.hpp @@ -0,0 +1,47 @@ +#ifndef Hadrons_FreeProp_hpp_ +#define Hadrons_FreeProp_hpp_ + +#include +#include +#include + +BEGIN_HADRONS_NAMESPACE + +/****************************************************************************** + * FreeProp * + ******************************************************************************/ +BEGIN_MODULE_NAMESPACE(MScalar) + +class FreePropPar: Serializable +{ +public: + GRID_SERIALIZABLE_CLASS_MEMBERS(FreePropPar, + std::string, source, + double, mass); +}; + +class TFreeProp: public Module +{ +public: + SCALAR_TYPE_ALIASES(SIMPL,); +public: + // constructor + TFreeProp(const std::string name); + // destructor + virtual ~TFreeProp(void) = default; + // dependency relation + virtual std::vector getInput(void); + virtual std::vector getOutput(void); + // setup + virtual void setup(void); + // execution + virtual void execute(void); +}; + +MODULE_REGISTER_NS(FreeProp, TFreeProp, MScalar); + +END_MODULE_NAMESPACE + +END_HADRONS_NAMESPACE + +#endif // Hadrons_FreeProp_hpp_ diff --git a/extras/Hadrons/modules.inc b/extras/Hadrons/modules.inc index 8b559024..b091c38b 100644 --- a/extras/Hadrons/modules.inc +++ b/extras/Hadrons/modules.inc @@ -2,7 +2,9 @@ modules_cc =\ Modules/MGauge/Load.cc \ Modules/MGauge/Random.cc \ Modules/MGauge/StochEm.cc \ - Modules/MGauge/Unit.cc + Modules/MGauge/Unit.cc \ + Modules/MScalar/ChargedProp.cc \ + Modules/MScalar/FreeProp.cc modules_hpp =\ Modules/MAction/DWF.hpp \ @@ -13,6 +15,8 @@ modules_hpp =\ Modules/MGauge/Random.hpp \ Modules/MGauge/StochEm.hpp \ Modules/MGauge/Unit.hpp \ + Modules/MScalar/ChargedProp.hpp \ + Modules/MScalar/FreeProp.hpp \ Modules/MSolver/RBPrecCG.hpp \ Modules/MSource/Point.hpp \ Modules/MSource/SeqGamma.hpp \ From 82b3f546970fde47b1e1220679f3c6c772d5eff0 Mon Sep 17 00:00:00 2001 From: Antonin Portelli Date: Thu, 5 Jan 2017 14:58:07 +0000 Subject: [PATCH 23/53] scalar free propagator fix --- lib/qcd/action/scalar/Scalar.h | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/lib/qcd/action/scalar/Scalar.h b/lib/qcd/action/scalar/Scalar.h index 194f6767..c053e15e 100644 --- a/lib/qcd/action/scalar/Scalar.h +++ b/lib/qcd/action/scalar/Scalar.h @@ -148,10 +148,11 @@ namespace QCD{ void Scalar::MomentumSpacePropagator(ScalarField &out, RealD m) { GridBase *grid = out._grid; - ScalarField kmu(grid); + ScalarField kmu(grid), one(grid); const unsigned int nd = grid->_ndimension; std::vector &l = grid->_fdimensions; + one = Complex(1.0,0.0); out = m*m; for(int mu = 0; mu < nd; mu++) { @@ -161,6 +162,7 @@ namespace QCD{ kmu = 2.*sin(.5*twoPiL*kmu); out = out + kmu*kmu; } + out = one/out; } template From 97843e2b5818667ab5f6802003bb4c04d6076503 Mon Sep 17 00:00:00 2001 From: Antonin Portelli Date: Thu, 5 Jan 2017 14:58:55 +0000 Subject: [PATCH 24/53] Hadrons: free scalar buffer fix and output --- extras/Hadrons/Modules/MScalar/FreeProp.cc | 34 ++++++++++++++++++--- extras/Hadrons/Modules/MScalar/FreeProp.hpp | 5 +-- 2 files changed, 32 insertions(+), 7 deletions(-) diff --git a/extras/Hadrons/Modules/MScalar/FreeProp.cc b/extras/Hadrons/Modules/MScalar/FreeProp.cc index 7419a954..ba85e041 100644 --- a/extras/Hadrons/Modules/MScalar/FreeProp.cc +++ b/extras/Hadrons/Modules/MScalar/FreeProp.cc @@ -1,11 +1,13 @@ #include +#define KERNAME "_" + getName() + "_momKernel" + using namespace Grid; using namespace Hadrons; using namespace MScalar; /****************************************************************************** -* TFreeProp implementation * +* TFreeProp implementation * ******************************************************************************/ // constructor ///////////////////////////////////////////////////////////////// TFreeProp::TFreeProp(const std::string name) @@ -30,6 +32,12 @@ std::vector TFreeProp::getOutput(void) // setup /////////////////////////////////////////////////////////////////////// void TFreeProp::setup(void) { + std::string kerName = KERNAME; + + if (!env().hasRegisteredObject(kerName)) + { + env().registerLattice(kerName); + } env().registerLattice(getName()); } @@ -39,13 +47,13 @@ void TFreeProp::execute(void) ScalarField &prop = *env().createLattice(getName()); ScalarField &source = *env().getObject(par().source); ScalarField *momKernel; - std::string kerName = "_" + getName() + "_momKernel"; - + std::string kerName = KERNAME; + if (!env().hasCreatedObject(kerName)) { LOG(Message) << "Caching momentum space free scalar propagator" - << "(mass= " << par().mass << ")..." << std::endl; - momKernel = env().template createLattice(kerName); + << " (mass= " << par().mass << ")..." << std::endl; + momKernel = env().createLattice(kerName); Scalar::MomentumSpacePropagator(*momKernel, par().mass); } else @@ -54,4 +62,20 @@ void TFreeProp::execute(void) } LOG(Message) << "Computing free scalar propagator..." << std::endl; Scalar::FreePropagator(source, prop, *momKernel); + + if (!par().output.empty()) + { + TextWriter writer(par().output + "." + + std::to_string(env().getTrajectory())); + std::vector buf; + std::vector result; + + sliceSum(prop, buf, Tp); + result.resize(buf.size()); + for (unsigned int t = 0; t < buf.size(); ++t) + { + result[t] = TensorRemove(buf[t]); + } + write(writer, "prop", result); + } } diff --git a/extras/Hadrons/Modules/MScalar/FreeProp.hpp b/extras/Hadrons/Modules/MScalar/FreeProp.hpp index 6a0cd930..81bb8121 100644 --- a/extras/Hadrons/Modules/MScalar/FreeProp.hpp +++ b/extras/Hadrons/Modules/MScalar/FreeProp.hpp @@ -8,7 +8,7 @@ BEGIN_HADRONS_NAMESPACE /****************************************************************************** - * FreeProp * + * FreeProp * ******************************************************************************/ BEGIN_MODULE_NAMESPACE(MScalar) @@ -17,7 +17,8 @@ class FreePropPar: Serializable public: GRID_SERIALIZABLE_CLASS_MEMBERS(FreePropPar, std::string, source, - double, mass); + double, mass, + std::string, output); }; class TFreeProp: public Module From fc760016b3e12b5fea25c8ca288525d4e2dad7c7 Mon Sep 17 00:00:00 2001 From: Antonin Portelli Date: Wed, 11 Jan 2017 18:39:58 +0000 Subject: [PATCH 25/53] More uniform cache name for scalar momentum propagators --- extras/Hadrons/Modules.hpp | 1 + extras/Hadrons/Modules/MScalar/FreeProp.cc | 22 ++++++++++----------- extras/Hadrons/Modules/MScalar/FreeProp.hpp | 2 ++ extras/Hadrons/Modules/MScalar/Scalar.hpp | 6 ++++++ extras/Hadrons/modules.inc | 1 + 5 files changed, 20 insertions(+), 12 deletions(-) create mode 100644 extras/Hadrons/Modules/MScalar/Scalar.hpp diff --git a/extras/Hadrons/Modules.hpp b/extras/Hadrons/Modules.hpp index ad31d2a7..a25419c5 100644 --- a/extras/Hadrons/Modules.hpp +++ b/extras/Hadrons/Modules.hpp @@ -8,6 +8,7 @@ #include #include #include +#include #include #include #include diff --git a/extras/Hadrons/Modules/MScalar/FreeProp.cc b/extras/Hadrons/Modules/MScalar/FreeProp.cc index ba85e041..f0a503ff 100644 --- a/extras/Hadrons/Modules/MScalar/FreeProp.cc +++ b/extras/Hadrons/Modules/MScalar/FreeProp.cc @@ -1,6 +1,5 @@ #include - -#define KERNAME "_" + getName() + "_momKernel" +#include using namespace Grid; using namespace Hadrons; @@ -32,11 +31,11 @@ std::vector TFreeProp::getOutput(void) // setup /////////////////////////////////////////////////////////////////////// void TFreeProp::setup(void) { - std::string kerName = KERNAME; + freeMomPropName_ = FREEMOMPROP(par().mass); - if (!env().hasRegisteredObject(kerName)) + if (!env().hasRegisteredObject(freeMomPropName_)) { - env().registerLattice(kerName); + env().registerLattice(freeMomPropName_); } env().registerLattice(getName()); } @@ -46,22 +45,21 @@ void TFreeProp::execute(void) { ScalarField &prop = *env().createLattice(getName()); ScalarField &source = *env().getObject(par().source); - ScalarField *momKernel; - std::string kerName = KERNAME; + ScalarField *freeMomProp; - if (!env().hasCreatedObject(kerName)) + if (!env().hasCreatedObject(freeMomPropName_)) { LOG(Message) << "Caching momentum space free scalar propagator" << " (mass= " << par().mass << ")..." << std::endl; - momKernel = env().createLattice(kerName); - Scalar::MomentumSpacePropagator(*momKernel, par().mass); + freeMomProp = env().createLattice(freeMomPropName_); + Scalar::MomentumSpacePropagator(*freeMomProp, par().mass); } else { - momKernel = env().getObject(kerName); + freeMomProp = env().getObject(freeMomPropName_); } LOG(Message) << "Computing free scalar propagator..." << std::endl; - Scalar::FreePropagator(source, prop, *momKernel); + Scalar::FreePropagator(source, prop, *freeMomProp); if (!par().output.empty()) { diff --git a/extras/Hadrons/Modules/MScalar/FreeProp.hpp b/extras/Hadrons/Modules/MScalar/FreeProp.hpp index 81bb8121..29f15eda 100644 --- a/extras/Hadrons/Modules/MScalar/FreeProp.hpp +++ b/extras/Hadrons/Modules/MScalar/FreeProp.hpp @@ -37,6 +37,8 @@ public: virtual void setup(void); // execution virtual void execute(void); +private: + std::string freeMomPropName_; }; MODULE_REGISTER_NS(FreeProp, TFreeProp, MScalar); diff --git a/extras/Hadrons/Modules/MScalar/Scalar.hpp b/extras/Hadrons/Modules/MScalar/Scalar.hpp new file mode 100644 index 00000000..db702ff2 --- /dev/null +++ b/extras/Hadrons/Modules/MScalar/Scalar.hpp @@ -0,0 +1,6 @@ +#ifndef Hadrons_Scalar_hpp_ +#define Hadrons_Scalar_hpp_ + +#define FREEMOMPROP(m) "_scalar_mom_prop_" + std::to_string(m) + +#endif // Hadrons_Scalar_hpp_ diff --git a/extras/Hadrons/modules.inc b/extras/Hadrons/modules.inc index b091c38b..dfbe85ff 100644 --- a/extras/Hadrons/modules.inc +++ b/extras/Hadrons/modules.inc @@ -17,6 +17,7 @@ modules_hpp =\ Modules/MGauge/Unit.hpp \ Modules/MScalar/ChargedProp.hpp \ Modules/MScalar/FreeProp.hpp \ + Modules/MScalar/Scalar.hpp \ Modules/MSolver/RBPrecCG.hpp \ Modules/MSource/Point.hpp \ Modules/MSource/SeqGamma.hpp \ From ad98b6193d4b08ad42c9da79370b4ccd7382b4cb Mon Sep 17 00:00:00 2001 From: Antonin Portelli Date: Wed, 11 Jan 2017 18:40:43 +0000 Subject: [PATCH 26/53] creating the necessary caches for the FFT EM scalar propagator --- extras/Hadrons/Modules/MScalar/ChargedProp.cc | 68 ++++++++++++++++++- .../Hadrons/Modules/MScalar/ChargedProp.hpp | 10 ++- 2 files changed, 74 insertions(+), 4 deletions(-) diff --git a/extras/Hadrons/Modules/MScalar/ChargedProp.cc b/extras/Hadrons/Modules/MScalar/ChargedProp.cc index 1137c6f0..1cd0cae6 100644 --- a/extras/Hadrons/Modules/MScalar/ChargedProp.cc +++ b/extras/Hadrons/Modules/MScalar/ChargedProp.cc @@ -1,4 +1,5 @@ #include +#include using namespace Grid; using namespace Hadrons; @@ -15,7 +16,7 @@ TChargedProp::TChargedProp(const std::string name) // dependencies/products /////////////////////////////////////////////////////// std::vector TChargedProp::getInput(void) { - std::vector in; + std::vector in = {par().source, par().emField}; return in; } @@ -30,11 +31,72 @@ std::vector TChargedProp::getOutput(void) // setup /////////////////////////////////////////////////////////////////////// void TChargedProp::setup(void) { - + freeMomPropName_ = FREEMOMPROP(par().mass); + shiftedMomPropName_.clear(); + for (unsigned int mu = 0; mu < env().getNd(); ++mu) + { + shiftedMomPropName_.push_back(freeMomPropName_ + "_" + + std::to_string(mu)); + } + if (!env().hasRegisteredObject(freeMomPropName_)) + { + env().registerLattice(freeMomPropName_); + } + if (!env().hasRegisteredObject(shiftedMomPropName_[0])) + { + for (unsigned int mu = 0; mu < env().getNd(); ++mu) + { + env().registerLattice(shiftedMomPropName_[mu]); + } + } + env().registerLattice(getName()); + } // execution /////////////////////////////////////////////////////////////////// void TChargedProp::execute(void) { - + ScalarField &prop = *env().createLattice(getName()); + ScalarField &source = *env().getObject(par().source); + ScalarField *freeMomProp; + std::vector shiftedMomProp; + Complex ci(0.0,1.0); + + if (!env().hasCreatedObject(freeMomPropName_)) + { + LOG(Message) << "Caching momentum space free scalar propagator" + << " (mass= " << par().mass << ")..." << std::endl; + freeMomProp = env().createLattice(freeMomPropName_); + Scalar::MomentumSpacePropagator(*freeMomProp, par().mass); + } + else + { + freeMomProp = env().getObject(freeMomPropName_); + } + if (!env().hasCreatedObject(shiftedMomPropName_[0])) + { + std::vector &l = env().getGrid()->_fdimensions; + + LOG(Message) << "Caching shifted momentum space free scalar propagator" + << " (mass= " << par().mass << ")..." << std::endl; + for (unsigned int mu = 0; mu < env().getNd(); ++mu) + { + Real twoPiL = M_PI*2./l[mu]; + + shiftedMomProp.push_back( + env().createLattice(shiftedMomPropName_[mu])); + LatticeCoordinate(*(shiftedMomProp[mu]), mu); + *(shiftedMomProp[mu]) = exp(ci*twoPiL*(*(shiftedMomProp[mu]))) + *(*freeMomProp); + } + } + else + { + for (unsigned int mu = 0; mu < env().getNd(); ++mu) + { + shiftedMomProp.push_back( + env().getObject(shiftedMomPropName_[mu])); + } + } + } diff --git a/extras/Hadrons/Modules/MScalar/ChargedProp.hpp b/extras/Hadrons/Modules/MScalar/ChargedProp.hpp index 7a60c2ad..91ea2355 100644 --- a/extras/Hadrons/Modules/MScalar/ChargedProp.hpp +++ b/extras/Hadrons/Modules/MScalar/ChargedProp.hpp @@ -16,11 +16,16 @@ class ChargedPropPar: Serializable { public: GRID_SERIALIZABLE_CLASS_MEMBERS(ChargedPropPar, - unsigned int, i); + std::string, emField, + std::string, source, + double, mass, + std::string, output); }; class TChargedProp: public Module { +public: + SCALAR_TYPE_ALIASES(SIMPL,); public: // constructor TChargedProp(const std::string name); @@ -33,6 +38,9 @@ public: virtual void setup(void); // execution virtual void execute(void); +private: + std::string freeMomPropName_; + std::vector shiftedMomPropName_; }; MODULE_REGISTER_NS(ChargedProp, TChargedProp, MScalar); From 889d828bc289d2ee4ea5939af29ff56fe7466db5 Mon Sep 17 00:00:00 2001 From: Antonin Portelli Date: Thu, 12 Jan 2017 18:17:44 +0000 Subject: [PATCH 27/53] Code cleaning --- extras/Hadrons/Modules/MScalar/ChargedProp.cc | 2 +- extras/Hadrons/Modules/MScalar/ChargedProp.hpp | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/extras/Hadrons/Modules/MScalar/ChargedProp.cc b/extras/Hadrons/Modules/MScalar/ChargedProp.cc index 1cd0cae6..ff53fa0b 100644 --- a/extras/Hadrons/Modules/MScalar/ChargedProp.cc +++ b/extras/Hadrons/Modules/MScalar/ChargedProp.cc @@ -6,7 +6,7 @@ using namespace Hadrons; using namespace MScalar; /****************************************************************************** -* TChargedProp implementation * +* TChargedProp implementation * ******************************************************************************/ // constructor ///////////////////////////////////////////////////////////////// TChargedProp::TChargedProp(const std::string name) diff --git a/extras/Hadrons/Modules/MScalar/ChargedProp.hpp b/extras/Hadrons/Modules/MScalar/ChargedProp.hpp index 91ea2355..001f6494 100644 --- a/extras/Hadrons/Modules/MScalar/ChargedProp.hpp +++ b/extras/Hadrons/Modules/MScalar/ChargedProp.hpp @@ -8,7 +8,7 @@ BEGIN_HADRONS_NAMESPACE /****************************************************************************** - * ChargedProp * + * Charged scalar propagator * ******************************************************************************/ BEGIN_MODULE_NAMESPACE(MScalar) From 65987a8a5810434b3f7ee54e8b6a4cf400108c74 Mon Sep 17 00:00:00 2001 From: Antonin Portelli Date: Thu, 12 Jan 2017 20:44:23 +0000 Subject: [PATCH 28/53] First implementation of the scalar QED propagator, runs but absolutely not checked --- extras/Hadrons/Modules/MScalar/ChargedProp.cc | 122 +++++++++++++++--- .../Hadrons/Modules/MScalar/ChargedProp.hpp | 13 +- 2 files changed, 115 insertions(+), 20 deletions(-) diff --git a/extras/Hadrons/Modules/MScalar/ChargedProp.cc b/extras/Hadrons/Modules/MScalar/ChargedProp.cc index ff53fa0b..dd260798 100644 --- a/extras/Hadrons/Modules/MScalar/ChargedProp.cc +++ b/extras/Hadrons/Modules/MScalar/ChargedProp.cc @@ -32,23 +32,28 @@ std::vector TChargedProp::getOutput(void) void TChargedProp::setup(void) { freeMomPropName_ = FREEMOMPROP(par().mass); - shiftedMomPropName_.clear(); + phaseName_.clear(); for (unsigned int mu = 0; mu < env().getNd(); ++mu) { - shiftedMomPropName_.push_back(freeMomPropName_ + "_" + phaseName_.push_back(freeMomPropName_ + "_" + std::to_string(mu)); } + GFSrcName_ = "_" + getName() + "_DinvSrc"; if (!env().hasRegisteredObject(freeMomPropName_)) { env().registerLattice(freeMomPropName_); } - if (!env().hasRegisteredObject(shiftedMomPropName_[0])) + if (!env().hasRegisteredObject(phaseName_[0])) { for (unsigned int mu = 0; mu < env().getNd(); ++mu) { - env().registerLattice(shiftedMomPropName_[mu]); + env().registerLattice(phaseName_[mu]); } } + if (!env().hasRegisteredObject(GFSrcName_)) + { + env().registerLattice(GFSrcName_); + } env().registerLattice(getName()); } @@ -56,24 +61,26 @@ void TChargedProp::setup(void) // execution /////////////////////////////////////////////////////////////////// void TChargedProp::execute(void) { + // CACHING ANALYTIC EXPRESSIONS ScalarField &prop = *env().createLattice(getName()); ScalarField &source = *env().getObject(par().source); - ScalarField *freeMomProp; - std::vector shiftedMomProp; - Complex ci(0.0,1.0); + Complex ci(0.0,1.0); + FFT fft(env().getGrid()); + // cache free scalar propagator if (!env().hasCreatedObject(freeMomPropName_)) { LOG(Message) << "Caching momentum space free scalar propagator" << " (mass= " << par().mass << ")..." << std::endl; - freeMomProp = env().createLattice(freeMomPropName_); - Scalar::MomentumSpacePropagator(*freeMomProp, par().mass); + freeMomProp_ = env().createLattice(freeMomPropName_); + Scalar::MomentumSpacePropagator(*freeMomProp_, par().mass); } else { - freeMomProp = env().getObject(freeMomPropName_); + freeMomProp_ = env().getObject(freeMomPropName_); } - if (!env().hasCreatedObject(shiftedMomPropName_[0])) + // cache phases + if (!env().hasCreatedObject(phaseName_[0])) { std::vector &l = env().getGrid()->_fdimensions; @@ -83,20 +90,99 @@ void TChargedProp::execute(void) { Real twoPiL = M_PI*2./l[mu]; - shiftedMomProp.push_back( - env().createLattice(shiftedMomPropName_[mu])); - LatticeCoordinate(*(shiftedMomProp[mu]), mu); - *(shiftedMomProp[mu]) = exp(ci*twoPiL*(*(shiftedMomProp[mu]))) - *(*freeMomProp); + phase_.push_back(env().createLattice(phaseName_[mu])); + LatticeCoordinate(*(phase_[mu]), mu); + *(phase_[mu]) = exp(ci*twoPiL*(*(phase_[mu]))); } } else { for (unsigned int mu = 0; mu < env().getNd(); ++mu) { - shiftedMomProp.push_back( - env().getObject(shiftedMomPropName_[mu])); + phase_.push_back(env().getObject(phaseName_[mu])); } } + // cache G*F*src + if (!env().hasCreatedObject(GFSrcName_)) + + { + GFSrc_ = env().createLattice(GFSrcName_); + fft.FFT_all_dim(*GFSrc_, source, FFT::forward); + *GFSrc_ = (*freeMomProp_)*(*GFSrc_); + } + else + { + GFSrc_ = env().getObject(GFSrcName_); + } + // PROPAGATOR CALCULATION + ScalarField buf(env().getGrid()); + ScalarField &GFSrc = *GFSrc_, &G = *freeMomProp_; + double q = par().charge; + + // G*F*Src + prop = GFSrc; + // - q*G*momD1*G*F*Src (momD1 = F*D1*Finv) + buf = GFSrc; + momD1(buf, fft); + buf = G*buf; + prop = prop - q*buf; + // + q^2*G*momD1*G*momD1*G*F*Src (here buf = G*momD1*G*F*Src) + momD1(buf, fft); + prop = prop + q*q*G*buf; + // + q^2*G*momD2*G*F*Src (momD1 = F*D2*Finv) + buf = GFSrc; + momD2(buf, fft); + prop = prop + q*q*G*buf; + // final FT + fft.FFT_all_dim(prop, prop, FFT::backward); +} + +void TChargedProp::momD1(ScalarField &s, FFT &fft) +{ + EmField &A = *env().getObject(par().emField); + ScalarField buf(env().getGrid()), Amu(env().getGrid()); + Complex ci(0.0,1.0); + + for (unsigned int mu = 0; mu < env().getNd(); ++mu) + { + Amu = peekLorentz(A, mu); + fft.FFT_all_dim(buf, s, FFT::backward); + buf = Amu*buf; + fft.FFT_all_dim(buf, buf, FFT::forward); + s = s + ci*adj(*phase_[mu])*buf; + } + for (unsigned int mu = 0; mu < env().getNd(); ++mu) + { + Amu = peekLorentz(A, mu); + buf = (*phase_[mu])*s; + fft.FFT_all_dim(buf, buf, FFT::backward); + buf = Amu*buf; + fft.FFT_all_dim(buf, buf, FFT::forward); + s = s - ci*buf; + } +} + +void TChargedProp::momD2(ScalarField &s, FFT &fft) +{ + EmField &A = *env().getObject(par().emField); + ScalarField buf(env().getGrid()), Amu(env().getGrid()); + + for (unsigned int mu = 0; mu < env().getNd(); ++mu) + { + Amu = peekLorentz(A, mu); + fft.FFT_all_dim(buf, s, FFT::backward); + buf = Amu*Amu*buf; + fft.FFT_all_dim(buf, buf, FFT::forward); + s = s + .5*adj(*phase_[mu])*buf; + } + for (unsigned int mu = 0; mu < env().getNd(); ++mu) + { + Amu = peekLorentz(A, mu); + buf = (*phase_[mu])*s; + fft.FFT_all_dim(buf, buf, FFT::backward); + buf = Amu*Amu*buf; + fft.FFT_all_dim(buf, buf, FFT::forward); + s = s + .5*buf; + } } diff --git a/extras/Hadrons/Modules/MScalar/ChargedProp.hpp b/extras/Hadrons/Modules/MScalar/ChargedProp.hpp index 001f6494..8bb5faa0 100644 --- a/extras/Hadrons/Modules/MScalar/ChargedProp.hpp +++ b/extras/Hadrons/Modules/MScalar/ChargedProp.hpp @@ -19,6 +19,7 @@ public: std::string, emField, std::string, source, double, mass, + double, charge, std::string, output); }; @@ -26,6 +27,8 @@ class TChargedProp: public Module { public: SCALAR_TYPE_ALIASES(SIMPL,); + typedef PhotonR::GaugeField EmField; + typedef PhotonR::GaugeLinkField EmComp; public: // constructor TChargedProp(const std::string name); @@ -39,8 +42,14 @@ public: // execution virtual void execute(void); private: - std::string freeMomPropName_; - std::vector shiftedMomPropName_; + void momD1(ScalarField &s, FFT &fft); + void momD2(ScalarField &s, FFT &fft); +private: + std::string freeMomPropName_, GFSrcName_; + std::vector phaseName_; + ScalarField *freeMomProp_, *GFSrc_; + std::vector phase_; + EmField *A; }; MODULE_REGISTER_NS(ChargedProp, TChargedProp, MScalar); From 92f8950a5658f75fa4e184fdffda492d5e45b200 Mon Sep 17 00:00:00 2001 From: Antonin Portelli Date: Fri, 13 Jan 2017 13:30:56 +0000 Subject: [PATCH 29/53] Charged scalar prop: cleaning and output --- extras/Hadrons/Modules/MScalar/ChargedProp.cc | 61 +++++++++++++------ 1 file changed, 42 insertions(+), 19 deletions(-) diff --git a/extras/Hadrons/Modules/MScalar/ChargedProp.cc b/extras/Hadrons/Modules/MScalar/ChargedProp.cc index dd260798..f8323705 100644 --- a/extras/Hadrons/Modules/MScalar/ChargedProp.cc +++ b/extras/Hadrons/Modules/MScalar/ChargedProp.cc @@ -35,8 +35,7 @@ void TChargedProp::setup(void) phaseName_.clear(); for (unsigned int mu = 0; mu < env().getNd(); ++mu) { - phaseName_.push_back(freeMomPropName_ + "_" - + std::to_string(mu)); + phaseName_.push_back("_shiftphase_" + std::to_string(mu)); } GFSrcName_ = "_" + getName() + "_DinvSrc"; if (!env().hasRegisteredObject(freeMomPropName_)) @@ -55,14 +54,12 @@ void TChargedProp::setup(void) env().registerLattice(GFSrcName_); } env().registerLattice(getName()); - } // execution /////////////////////////////////////////////////////////////////// void TChargedProp::execute(void) { // CACHING ANALYTIC EXPRESSIONS - ScalarField &prop = *env().createLattice(getName()); ScalarField &source = *env().getObject(par().source); Complex ci(0.0,1.0); FFT fft(env().getGrid()); @@ -79,13 +76,24 @@ void TChargedProp::execute(void) { freeMomProp_ = env().getObject(freeMomPropName_); } + // cache G*F*src + if (!env().hasCreatedObject(GFSrcName_)) + + { + GFSrc_ = env().createLattice(GFSrcName_); + fft.FFT_all_dim(*GFSrc_, source, FFT::forward); + *GFSrc_ = (*freeMomProp_)*(*GFSrc_); + } + else + { + GFSrc_ = env().getObject(GFSrcName_); + } // cache phases if (!env().hasCreatedObject(phaseName_[0])) { std::vector &l = env().getGrid()->_fdimensions; - LOG(Message) << "Caching shifted momentum space free scalar propagator" - << " (mass= " << par().mass << ")..." << std::endl; + LOG(Message) << "Caching shift phases..." << std::endl; for (unsigned int mu = 0; mu < env().getNd(); ++mu) { Real twoPiL = M_PI*2./l[mu]; @@ -102,20 +110,13 @@ void TChargedProp::execute(void) phase_.push_back(env().getObject(phaseName_[mu])); } } - // cache G*F*src - if (!env().hasCreatedObject(GFSrcName_)) - - { - GFSrc_ = env().createLattice(GFSrcName_); - fft.FFT_all_dim(*GFSrc_, source, FFT::forward); - *GFSrc_ = (*freeMomProp_)*(*GFSrc_); - } - else - { - GFSrc_ = env().getObject(GFSrcName_); - } - + // PROPAGATOR CALCULATION + LOG(Message) << "Computing charged scalar propagator" + << " (mass= " << par().mass + << ", charge= " << par().charge << ")..." << std::endl; + + ScalarField &prop = *env().createLattice(getName()); ScalarField buf(env().getGrid()); ScalarField &GFSrc = *GFSrc_, &G = *freeMomProp_; double q = par().charge; @@ -136,6 +137,28 @@ void TChargedProp::execute(void) prop = prop + q*q*G*buf; // final FT fft.FFT_all_dim(prop, prop, FFT::backward); + + // OUTPUT IF NECESSARY + if (!par().output.empty()) + { + std::string filename = par().output + "." + + std::to_string(env().getTrajectory()); + + LOG(Message) << "Saving zero-momentum projection to '" + << filename << "'..." << std::endl; + + TextWriter writer(filename); + std::vector vecBuf; + std::vector result; + + sliceSum(prop, vecBuf, Tp); + result.resize(vecBuf.size()); + for (unsigned int t = 0; t < vecBuf.size(); ++t) + { + result[t] = TensorRemove(vecBuf[t]); + } + write(writer, "prop", result); + } } void TChargedProp::momD1(ScalarField &s, FFT &fft) From ae99e99da235ebf5d9a149ef4a7ad0b93d1f7474 Mon Sep 17 00:00:00 2001 From: James Harrison Date: Mon, 23 Jan 2017 17:27:50 +0000 Subject: [PATCH 30/53] Fixed bug in ChargedProp --- extras/Hadrons/Modules/MScalar/ChargedProp.cc | 40 ++++++++++++------- 1 file changed, 26 insertions(+), 14 deletions(-) diff --git a/extras/Hadrons/Modules/MScalar/ChargedProp.cc b/extras/Hadrons/Modules/MScalar/ChargedProp.cc index f8323705..d88fdc45 100644 --- a/extras/Hadrons/Modules/MScalar/ChargedProp.cc +++ b/extras/Hadrons/Modules/MScalar/ChargedProp.cc @@ -123,18 +123,22 @@ void TChargedProp::execute(void) // G*F*Src prop = GFSrc; + // - q*G*momD1*G*F*Src (momD1 = F*D1*Finv) buf = GFSrc; momD1(buf, fft); buf = G*buf; prop = prop - q*buf; + // + q^2*G*momD1*G*momD1*G*F*Src (here buf = G*momD1*G*F*Src) momD1(buf, fft); prop = prop + q*q*G*buf; - // + q^2*G*momD2*G*F*Src (momD1 = F*D2*Finv) + + // - q^2*G*momD2*G*F*Src (momD2 = F*D2*Finv) buf = GFSrc; momD2(buf, fft); - prop = prop + q*q*G*buf; + prop = prop - q*q*G*buf; + // final FT fft.FFT_all_dim(prop, prop, FFT::backward); @@ -164,16 +168,18 @@ void TChargedProp::execute(void) void TChargedProp::momD1(ScalarField &s, FFT &fft) { EmField &A = *env().getObject(par().emField); - ScalarField buf(env().getGrid()), Amu(env().getGrid()); + ScalarField buf(env().getGrid()), fs(env().getGrid()), result(env().getGrid()), Amu(env().getGrid()); Complex ci(0.0,1.0); - + + result = zero; + + fft.FFT_all_dim(fs, s, FFT::backward); for (unsigned int mu = 0; mu < env().getNd(); ++mu) { Amu = peekLorentz(A, mu); - fft.FFT_all_dim(buf, s, FFT::backward); - buf = Amu*buf; + buf = Amu*fs; fft.FFT_all_dim(buf, buf, FFT::forward); - s = s + ci*adj(*phase_[mu])*buf; + result = result + ci*adj(*phase_[mu])*buf; } for (unsigned int mu = 0; mu < env().getNd(); ++mu) { @@ -182,22 +188,26 @@ void TChargedProp::momD1(ScalarField &s, FFT &fft) fft.FFT_all_dim(buf, buf, FFT::backward); buf = Amu*buf; fft.FFT_all_dim(buf, buf, FFT::forward); - s = s - ci*buf; + result = result - ci*buf; } + + s = result; } void TChargedProp::momD2(ScalarField &s, FFT &fft) { EmField &A = *env().getObject(par().emField); - ScalarField buf(env().getGrid()), Amu(env().getGrid()); + ScalarField buf(env().getGrid()), fs(env().getGrid()), result(env().getGrid()), Amu(env().getGrid()); + + result = zero; + fft.FFT_all_dim(fs, s, FFT::backward); for (unsigned int mu = 0; mu < env().getNd(); ++mu) { - Amu = peekLorentz(A, mu); - fft.FFT_all_dim(buf, s, FFT::backward); - buf = Amu*Amu*buf; + Amu = peekLorentz(A, mu); + buf = Amu*Amu*fs; fft.FFT_all_dim(buf, buf, FFT::forward); - s = s + .5*adj(*phase_[mu])*buf; + result = result + .5*adj(*phase_[mu])*buf; } for (unsigned int mu = 0; mu < env().getNd(); ++mu) { @@ -206,6 +216,8 @@ void TChargedProp::momD2(ScalarField &s, FFT &fft) fft.FFT_all_dim(buf, buf, FFT::backward); buf = Amu*Amu*buf; fft.FFT_all_dim(buf, buf, FFT::forward); - s = s + .5*buf; + result = result + .5*buf; } + + s = result; } From f65a585236f420d7ed966b8f9e0b7cbfb0857d8c Mon Sep 17 00:00:00 2001 From: James Harrison Date: Thu, 26 Jan 2017 15:02:30 +0000 Subject: [PATCH 31/53] ChargedProp: Switch to HDF5 output --- extras/Hadrons/Modules/MScalar/ChargedProp.cc | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/extras/Hadrons/Modules/MScalar/ChargedProp.cc b/extras/Hadrons/Modules/MScalar/ChargedProp.cc index d88fdc45..f2890b2a 100644 --- a/extras/Hadrons/Modules/MScalar/ChargedProp.cc +++ b/extras/Hadrons/Modules/MScalar/ChargedProp.cc @@ -151,7 +151,7 @@ void TChargedProp::execute(void) LOG(Message) << "Saving zero-momentum projection to '" << filename << "'..." << std::endl; - TextWriter writer(filename); + Hdf5Writer writer(filename); std::vector vecBuf; std::vector result; @@ -161,6 +161,7 @@ void TChargedProp::execute(void) { result[t] = TensorRemove(vecBuf[t]); } + write(writer, "charge", q); write(writer, "prop", result); } } From ee93f0218bebd84d98211d0ed87cff48951d76d7 Mon Sep 17 00:00:00 2001 From: James Harrison Date: Fri, 27 Jan 2017 12:22:48 +0000 Subject: [PATCH 32/53] ChargedProp: remove ScalarField fs --- extras/Hadrons/Modules/MScalar/ChargedProp.cc | 38 ++++++++++--------- 1 file changed, 20 insertions(+), 18 deletions(-) diff --git a/extras/Hadrons/Modules/MScalar/ChargedProp.cc b/extras/Hadrons/Modules/MScalar/ChargedProp.cc index f2890b2a..40d4504c 100644 --- a/extras/Hadrons/Modules/MScalar/ChargedProp.cc +++ b/extras/Hadrons/Modules/MScalar/ChargedProp.cc @@ -169,19 +169,12 @@ void TChargedProp::execute(void) void TChargedProp::momD1(ScalarField &s, FFT &fft) { EmField &A = *env().getObject(par().emField); - ScalarField buf(env().getGrid()), fs(env().getGrid()), result(env().getGrid()), Amu(env().getGrid()); + ScalarField buf(env().getGrid()), result(env().getGrid()), + Amu(env().getGrid()); Complex ci(0.0,1.0); result = zero; - fft.FFT_all_dim(fs, s, FFT::backward); - for (unsigned int mu = 0; mu < env().getNd(); ++mu) - { - Amu = peekLorentz(A, mu); - buf = Amu*fs; - fft.FFT_all_dim(buf, buf, FFT::forward); - result = result + ci*adj(*phase_[mu])*buf; - } for (unsigned int mu = 0; mu < env().getNd(); ++mu) { Amu = peekLorentz(A, mu); @@ -191,6 +184,14 @@ void TChargedProp::momD1(ScalarField &s, FFT &fft) fft.FFT_all_dim(buf, buf, FFT::forward); result = result - ci*buf; } + fft.FFT_all_dim(s, s, FFT::backward); + for (unsigned int mu = 0; mu < env().getNd(); ++mu) + { + Amu = peekLorentz(A, mu); + buf = Amu*s; + fft.FFT_all_dim(buf, buf, FFT::forward); + result = result + ci*adj(*phase_[mu])*buf; + } s = result; } @@ -198,18 +199,11 @@ void TChargedProp::momD1(ScalarField &s, FFT &fft) void TChargedProp::momD2(ScalarField &s, FFT &fft) { EmField &A = *env().getObject(par().emField); - ScalarField buf(env().getGrid()), fs(env().getGrid()), result(env().getGrid()), Amu(env().getGrid()); + ScalarField buf(env().getGrid()), result(env().getGrid()), + Amu(env().getGrid()); result = zero; - fft.FFT_all_dim(fs, s, FFT::backward); - for (unsigned int mu = 0; mu < env().getNd(); ++mu) - { - Amu = peekLorentz(A, mu); - buf = Amu*Amu*fs; - fft.FFT_all_dim(buf, buf, FFT::forward); - result = result + .5*adj(*phase_[mu])*buf; - } for (unsigned int mu = 0; mu < env().getNd(); ++mu) { Amu = peekLorentz(A, mu); @@ -219,6 +213,14 @@ void TChargedProp::momD2(ScalarField &s, FFT &fft) fft.FFT_all_dim(buf, buf, FFT::forward); result = result + .5*buf; } + fft.FFT_all_dim(s, s, FFT::backward); + for (unsigned int mu = 0; mu < env().getNd(); ++mu) + { + Amu = peekLorentz(A, mu); + buf = Amu*Amu*s; + fft.FFT_all_dim(buf, buf, FFT::forward); + result = result + .5*adj(*phase_[mu])*buf; + } s = result; } From b39f0d1fb675f453e710ed953583bb68bfe2b18f Mon Sep 17 00:00:00 2001 From: Antonin Portelli Date: Fri, 27 Jan 2017 18:12:35 -0800 Subject: [PATCH 33/53] Hadrons: default I/O to HDF5 if possible, XML otherwise --- extras/Hadrons/Global.hpp | 9 +++++++++ extras/Hadrons/Modules/MScalar/ChargedProp.cc | 2 +- 2 files changed, 10 insertions(+), 1 deletion(-) diff --git a/extras/Hadrons/Global.hpp b/extras/Hadrons/Global.hpp index bcb282fc..8dbb08ca 100644 --- a/extras/Hadrons/Global.hpp +++ b/extras/Hadrons/Global.hpp @@ -160,6 +160,15 @@ std::string typeName(void) return typeName(typeIdPt()); } +// default writers/readers +#ifdef HAVE_HDF5 +typedef Hdf5Reader CorrReader; +typedef Hdf5Writer CorrWriter; +#else +typedef XmlReader CorrReader; +typedef XmlWriter CorrWriter; +#endif + END_HADRONS_NAMESPACE #endif // Hadrons_Global_hpp_ diff --git a/extras/Hadrons/Modules/MScalar/ChargedProp.cc b/extras/Hadrons/Modules/MScalar/ChargedProp.cc index 40d4504c..dc6481f3 100644 --- a/extras/Hadrons/Modules/MScalar/ChargedProp.cc +++ b/extras/Hadrons/Modules/MScalar/ChargedProp.cc @@ -151,7 +151,7 @@ void TChargedProp::execute(void) LOG(Message) << "Saving zero-momentum projection to '" << filename << "'..." << std::endl; - Hdf5Writer writer(filename); + CorrWriter writer(filename); std::vector vecBuf; std::vector result; From 140741875555fb9c788e78bb8b6080e480776c0f Mon Sep 17 00:00:00 2001 From: Antonin Portelli Date: Thu, 13 Apr 2017 15:32:30 +0100 Subject: [PATCH 34/53] Old qed-fvol program build disabled --- extras/Makefile.am | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/extras/Makefile.am b/extras/Makefile.am index 416a9fc8..d8c2b675 100644 --- a/extras/Makefile.am +++ b/extras/Makefile.am @@ -1 +1 @@ -SUBDIRS = Hadrons qed-fvol \ No newline at end of file +SUBDIRS = Hadrons \ No newline at end of file From 10f2872aae87626df0d2ec2a591df88bfbfa68a0 Mon Sep 17 00:00:00 2001 From: Guido Cossu Date: Mon, 15 May 2017 15:51:16 +0100 Subject: [PATCH 35/53] Faster exponentiation for lattice fields --- lib/lattice/Lattice_unary.h | 13 ++++++++++++- 1 file changed, 12 insertions(+), 1 deletion(-) diff --git a/lib/lattice/Lattice_unary.h b/lib/lattice/Lattice_unary.h index f5b324ec..0afeaa32 100644 --- a/lib/lattice/Lattice_unary.h +++ b/lib/lattice/Lattice_unary.h @@ -66,10 +66,21 @@ namespace Grid { Lattice ret(rhs._grid); ret.checkerboard = rhs.checkerboard; conformable(ret,rhs); + obj unit(1.0); parallel_for(int ss=0;ssoSites();ss++){ - ret._odata[ss]=Exponentiate(rhs._odata[ss],alpha, Nexp); + //ret._odata[ss]=Exponentiate(rhs._odata[ss],alpha, Nexp); + ret._odata[ss] = unit; + for(int i=Nexp; i>=1;--i) + ret._odata[ss] = unit + ret._odata[ss]*rhs._odata[ss]*(alpha/RealD(i)); + } + return ret; + + + + + } From bc862ce3ab7b642fba1874f65d867fcbc61c8df2 Mon Sep 17 00:00:00 2001 From: Guido Cossu Date: Thu, 18 May 2017 14:44:56 +0100 Subject: [PATCH 36/53] Fixing an allocation issue in Benchmark_comms --- benchmarks/Benchmark_comms.cc | 19 ++++++++++--------- 1 file changed, 10 insertions(+), 9 deletions(-) diff --git a/benchmarks/Benchmark_comms.cc b/benchmarks/Benchmark_comms.cc index 99ab190b..a49e3033 100644 --- a/benchmarks/Benchmark_comms.cc +++ b/benchmarks/Benchmark_comms.cc @@ -40,15 +40,17 @@ int main (int argc, char ** argv) int threads = GridThread::GetThreads(); std::cout<1) nmu++; + + std::cout< xbuf(8); std::vector rbuf(8); + Grid.ShmInitGeneric(); Grid.ShmBufferFreeAll(); for(int d=0;d<8;d++){ xbuf[d] = (HalfSpinColourVectorD *)Grid.ShmBufferMalloc(lat*lat*lat*Ls*sizeof(HalfSpinColourVectorD)); @@ -269,13 +271,11 @@ int main (int argc, char ** argv) double time = stop-start; // microseconds - std::cout< xbuf(8); std::vector rbuf(8); + Grid.ShmInitGeneric(); Grid.ShmBufferFreeAll(); for(int d=0;d<8;d++){ xbuf[d] = (HalfSpinColourVectorD *)Grid.ShmBufferMalloc(lat*lat*lat*Ls*sizeof(HalfSpinColourVectorD)); @@ -353,7 +354,7 @@ int main (int argc, char ** argv) double time = stop-start; // microseconds - std::cout< Date: Thu, 18 May 2017 17:48:11 +0100 Subject: [PATCH 37/53] Reverting change in Bechmark_comms. Keeping 300 iterations --- benchmarks/Benchmark_comms.cc | 4 +--- 1 file changed, 1 insertion(+), 3 deletions(-) diff --git a/benchmarks/Benchmark_comms.cc b/benchmarks/Benchmark_comms.cc index a49e3033..319f01dc 100644 --- a/benchmarks/Benchmark_comms.cc +++ b/benchmarks/Benchmark_comms.cc @@ -40,7 +40,7 @@ int main (int argc, char ** argv) int threads = GridThread::GetThreads(); std::cout<1) nmu++; @@ -213,7 +213,6 @@ int main (int argc, char ** argv) std::vector xbuf(8); std::vector rbuf(8); - Grid.ShmInitGeneric(); Grid.ShmBufferFreeAll(); for(int d=0;d<8;d++){ xbuf[d] = (HalfSpinColourVectorD *)Grid.ShmBufferMalloc(lat*lat*lat*Ls*sizeof(HalfSpinColourVectorD)); @@ -293,7 +292,6 @@ int main (int argc, char ** argv) std::vector xbuf(8); std::vector rbuf(8); - Grid.ShmInitGeneric(); Grid.ShmBufferFreeAll(); for(int d=0;d<8;d++){ xbuf[d] = (HalfSpinColourVectorD *)Grid.ShmBufferMalloc(lat*lat*lat*Ls*sizeof(HalfSpinColourVectorD)); From 8e19c99c7d87bd4af0e9ecb9dd72c68f5368f891 Mon Sep 17 00:00:00 2001 From: Guido Cossu Date: Thu, 18 May 2017 19:07:35 +0100 Subject: [PATCH 38/53] Adding more statistical info in the Benchmark_comms --- benchmarks/Benchmark_comms.cc | 69 +++++++++++++++++++++++------------ 1 file changed, 46 insertions(+), 23 deletions(-) diff --git a/benchmarks/Benchmark_comms.cc b/benchmarks/Benchmark_comms.cc index 319f01dc..d4c59785 100644 --- a/benchmarks/Benchmark_comms.cc +++ b/benchmarks/Benchmark_comms.cc @@ -31,6 +31,17 @@ using namespace std; using namespace Grid; using namespace Grid::QCD; +void statistics(std::vector v, double &mean, double &std_err){ + double sum = std::accumulate(v.begin(), v.end(), 0.0); + mean = sum / v.size(); + + std::vector diff(v.size()); + std::transform(v.begin(), v.end(), diff.begin(), [mean](double x) { return x - mean; }); + double sq_sum = std::inner_product(diff.begin(), diff.end(), diff.begin(), 0.0); + std_err = std::sqrt(sq_sum / (v.size()*(v.size() - 1))); +} + + int main (int argc, char ** argv) { Grid_init(&argc,&argv); @@ -40,12 +51,14 @@ int main (int argc, char ** argv) int threads = GridThread::GetThreads(); std::cout<1) nmu++; - + std::cout << GridLogMessage << "Number of iterations to average: "<< Nloop << std::endl; + std::vector t_time(Nloop); + double mean_time, std_err_time; std::cout< requests; @@ -104,18 +117,21 @@ int main (int argc, char ** argv) } Grid.SendToRecvFromComplete(requests); Grid.Barrier(); - + double stop=usecond(); + t_time[i] = stop-start; // microseconds } - double stop=usecond(); + + statistics(t_time, mean_time, std_err_time); double dbytes = bytes; - double xbytes = Nloop*dbytes*2.0*ncomm; + double xbytes = dbytes*2.0*ncomm; double rbytes = xbytes; double bidibytes = xbytes+rbytes; - double time = stop-start; // microseconds + std::cout< requests; @@ -259,18 +278,19 @@ int main (int argc, char ** argv) } Grid.StencilSendToRecvFromComplete(requests); Grid.Barrier(); + double stop=usecond(); + t_time[i] = stop-start; // microseconds } - double stop=usecond(); double dbytes = bytes; double xbytes = Nloop*dbytes*2.0*ncomm; double rbytes = xbytes; double bidibytes = xbytes+rbytes; - double time = stop-start; // microseconds - - std::cout< requests; @@ -340,19 +360,22 @@ int main (int argc, char ** argv) } } - Grid.Barrier(); + Grid.Barrier(); + double stop=usecond(); + t_time[i] = stop-start; // microseconds } - double stop=usecond(); double dbytes = bytes; double xbytes = Nloop*dbytes*2.0*ncomm; double rbytes = xbytes; double bidibytes = xbytes+rbytes; - double time = stop-start; // microseconds - std::cout< Date: Fri, 19 May 2017 10:55:04 +0100 Subject: [PATCH 39/53] Adding more statistics to the Benchmark_comms. Min and max --- benchmarks/Benchmark_comms.cc | 88 ++++++++++++++++++++++++----------- 1 file changed, 61 insertions(+), 27 deletions(-) diff --git a/benchmarks/Benchmark_comms.cc b/benchmarks/Benchmark_comms.cc index d4c59785..ce881ef6 100644 --- a/benchmarks/Benchmark_comms.cc +++ b/benchmarks/Benchmark_comms.cc @@ -31,16 +31,31 @@ using namespace std; using namespace Grid; using namespace Grid::QCD; -void statistics(std::vector v, double &mean, double &std_err){ +struct time_statistics{ + double mean; + double err; + double min; + double max; + + void statistics(std::vector v){ double sum = std::accumulate(v.begin(), v.end(), 0.0); mean = sum / v.size(); std::vector diff(v.size()); - std::transform(v.begin(), v.end(), diff.begin(), [mean](double x) { return x - mean; }); + std::transform(v.begin(), v.end(), diff.begin(), [=](double x) { return x - mean; }); double sq_sum = std::inner_product(diff.begin(), diff.end(), diff.begin(), 0.0); - std_err = std::sqrt(sq_sum / (v.size()*(v.size() - 1))); -} + err = std::sqrt(sq_sum / (v.size()*(v.size() - 1))); + auto result = std::minmax_element(v.begin(), v.end()); + min = *result.first; + max = *result.second; +} +}; + +void header(){ + std::cout < t_time(Nloop); - double mean_time, std_err_time; + time_statistics timestat; std::cout< Date: Fri, 19 May 2017 16:39:36 +0100 Subject: [PATCH 40/53] Fixing a compilation error for generic SIMD --- lib/simd/Grid_generic.h | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/lib/simd/Grid_generic.h b/lib/simd/Grid_generic.h index bdf5aa01..e1d5f894 100644 --- a/lib/simd/Grid_generic.h +++ b/lib/simd/Grid_generic.h @@ -281,8 +281,8 @@ namespace Optimization { struct PrecisionChange { static inline vech StoH (const vecf &a,const vecf &b) { -#ifdef USE_FP16 vech ret; +#ifdef USE_FP16 vech *ha = (vech *)&a; vech *hb = (vech *)&b; const int nf = W::r; @@ -493,6 +493,8 @@ namespace Optimization { return a; } + + #undef acc // EIGEN compatibility } ////////////////////////////////////////////////////////////////////////////////////// From ab3596d4d368c9112a3d169d3cd8459d0f072f05 Mon Sep 17 00:00:00 2001 From: Guido Cossu Date: Thu, 25 May 2017 12:07:47 +0100 Subject: [PATCH 41/53] Using Cayley-Hamilton form for the exponential of SU(3) matrices --- lib/lattice/Lattice_unary.h | 9 +-- lib/qcd/action/gauge/GaugeImplTypes.h | 21 ++++-- lib/qcd/smearing/StoutSmearing.h | 2 + lib/tensors/Tensor_exp.h | 103 ++++++++++++++++++++++---- 4 files changed, 106 insertions(+), 29 deletions(-) diff --git a/lib/lattice/Lattice_unary.h b/lib/lattice/Lattice_unary.h index 0afeaa32..44b7b4f1 100644 --- a/lib/lattice/Lattice_unary.h +++ b/lib/lattice/Lattice_unary.h @@ -62,17 +62,12 @@ namespace Grid { return ret; } - template Lattice expMat(const Lattice &rhs, ComplexD alpha, Integer Nexp = DEFAULT_MAT_EXP){ + template Lattice expMat(const Lattice &rhs, RealD alpha, Integer Nexp = DEFAULT_MAT_EXP){ Lattice ret(rhs._grid); ret.checkerboard = rhs.checkerboard; conformable(ret,rhs); - obj unit(1.0); parallel_for(int ss=0;ssoSites();ss++){ - //ret._odata[ss]=Exponentiate(rhs._odata[ss],alpha, Nexp); - ret._odata[ss] = unit; - for(int i=Nexp; i>=1;--i) - ret._odata[ss] = unit + ret._odata[ss]*rhs._odata[ss]*(alpha/RealD(i)); - + ret._odata[ss]=Exponentiate(rhs._odata[ss],alpha, Nexp); } return ret; diff --git a/lib/qcd/action/gauge/GaugeImplTypes.h b/lib/qcd/action/gauge/GaugeImplTypes.h index 9d36eead..29e79581 100644 --- a/lib/qcd/action/gauge/GaugeImplTypes.h +++ b/lib/qcd/action/gauge/GaugeImplTypes.h @@ -59,7 +59,7 @@ public: typedef iImplGaugeLink SiteLink; typedef iImplGaugeField SiteField; - typedef Lattice LinkField; + typedef Lattice LinkField; typedef Lattice Field; // Guido: we can probably separate the types from the HMC functions @@ -80,7 +80,7 @@ public: /////////////////////////////////////////////////////////// // Move these to another class - // HMC auxiliary functions + // HMC auxiliary functions static inline void generate_momenta(Field &P, GridParallelRNG &pRNG) { // specific for SU gauge fields LinkField Pmu(P._grid); @@ -92,14 +92,19 @@ public: } static inline Field projectForce(Field &P) { return Ta(P); } - + static inline void update_field(Field& P, Field& U, double ep){ - for (int mu = 0; mu < Nd; mu++) { - auto Umu = PeekIndex(U, mu); - auto Pmu = PeekIndex(P, mu); - Umu = expMat(Pmu, ep, Nexp) * Umu; - PokeIndex(U, ProjectOnGroup(Umu), mu); + //static std::chrono::duration diff; + + //auto start = std::chrono::high_resolution_clock::now(); + parallel_for(int ss=0;ssoSites();ss++){ + for (int mu = 0; mu < Nd; mu++) + U[ss]._internal[mu] = ProjectOnGroup(Exponentiate(P[ss]._internal[mu], ep, Nexp) * U[ss]._internal[mu]); } + + //auto end = std::chrono::high_resolution_clock::now(); + // diff += end - start; + // std::cout << "Time to exponentiate matrix " << diff.count() << " s\n"; } static inline RealD FieldSquareNorm(Field& U){ diff --git a/lib/qcd/smearing/StoutSmearing.h b/lib/qcd/smearing/StoutSmearing.h index b50ffe21..bfc37d0d 100644 --- a/lib/qcd/smearing/StoutSmearing.h +++ b/lib/qcd/smearing/StoutSmearing.h @@ -58,6 +58,8 @@ class Smear_Stout : public Smear { SmearBase->smear(C, U); }; + + // Repetion of code here (use the Tensor_exp.h function) void exponentiate_iQ(GaugeLinkField& e_iQ, const GaugeLinkField& iQ) const { // Put this outside // only valid for SU(3) matrices diff --git a/lib/tensors/Tensor_exp.h b/lib/tensors/Tensor_exp.h index ad0d2071..64c06bb5 100644 --- a/lib/tensors/Tensor_exp.h +++ b/lib/tensors/Tensor_exp.h @@ -37,30 +37,105 @@ namespace Grid { /////////////////////////////////////////////// - template inline iScalar Exponentiate(const iScalar&r, ComplexD alpha , Integer Nexp = DEFAULT_MAT_EXP) + template inline iScalar Exponentiate(const iScalar&r, RealD alpha , Integer Nexp = DEFAULT_MAT_EXP) { iScalar ret; ret._internal = Exponentiate(r._internal, alpha, Nexp); return ret; } - - template::TensorLevel == 0 >::type * =nullptr> - inline iMatrix Exponentiate(const iMatrix &arg, ComplexD alpha , Integer Nexp = DEFAULT_MAT_EXP ) +template inline iVector Exponentiate(const iVector&r, RealD alpha , Integer Nexp = DEFAULT_MAT_EXP) { - iMatrix unit(1.0); - iMatrix temp(unit); - - for(int i=Nexp; i>=1;--i){ - temp *= alpha/ComplexD(i); - temp = unit + temp*arg; - } - - return temp; - + iVector ret; + for (int i = 0; i < N; i++) + ret._internal[i] = Exponentiate(r._internal[i], alpha, Nexp); + return ret; } + // Specialisation: Cayley-Hamilton exponential for SU(3) + template::TensorLevel == 0>::type * =nullptr> + inline iMatrix Exponentiate(const iMatrix &arg, RealD alpha , Integer Nexp = DEFAULT_MAT_EXP ) + { + // for SU(3) 2x faster than the std implementation using Nexp=12 + // notice that it actually computes + // exp ( input matrix ) + // the i sign is coming from outside + // input matrix is anti-hermitian NOT hermitian + typedef iMatrix mat; + typedef iScalar scalar; + mat unit(1.0); + mat temp(unit); + const Complex one_over_three = 1.0 / 3.0; + const Complex one_over_two = 1.0 / 2.0; + + scalar c0, c1, tmp, c0max, theta, u, w; + scalar xi0, u2, w2, cosw; + scalar fden, h0, h1, h2; + scalar e2iu, emiu, ixi0, qt; + scalar f0, f1, f2; + scalar unity(1.0); + + mat iQ2 = arg*arg*alpha*alpha; + mat iQ3 = arg*iQ2*alpha; + // sign in c0 from the conventions on the Ta + c0 = -imag( trace(iQ3) ) * one_over_three; + c1 = -real( trace(iQ2) ) * one_over_two; + + // Cayley Hamilton checks to machine precision, tested + tmp = c1 * one_over_three; + c0max = 2.0 * pow(tmp, 1.5); + + theta = acos(c0 / c0max) * one_over_three; + u = sqrt(tmp) * cos(theta); + w = sqrt(c1) * sin(theta); + + xi0 = sin(w) / w; + u2 = u * u; + w2 = w * w; + cosw = cos(w); + + ixi0 = timesI(xi0); + emiu = cos(u) - timesI(sin(u)); + e2iu = cos(2.0 * u) + timesI(sin(2.0 * u)); + + h0 = e2iu * (u2 - w2) + + emiu * ((8.0 * u2 * cosw) + (2.0 * u * (3.0 * u2 + w2) * ixi0)); + h1 = e2iu * (2.0 * u) - emiu * ((2.0 * u * cosw) - (3.0 * u2 - w2) * ixi0); + h2 = e2iu - emiu * (cosw + (3.0 * u) * ixi0); + + fden = unity / (9.0 * u2 - w2); // reals + f0 = h0 * fden; + f1 = h1 * fden; + f2 = h2 * fden; + + return (f0 * unit + timesMinusI(f1) * arg*alpha - f2 * iQ2); + } + + + +// General exponential +template::TensorLevel == 0 >::type * =nullptr> + inline iMatrix Exponentiate(const iMatrix &arg, RealD alpha , Integer Nexp = DEFAULT_MAT_EXP ) + { + // notice that it actually computes + // exp ( input matrix ) + // the i sign is coming from outside + // input matrix is anti-hermitian NOT hermitian + typedef iMatrix mat; + mat unit(1.0); + mat temp(unit); + for(int i=Nexp; i>=1;--i){ + temp *= alpha/ComplexD(i); + temp = unit + temp*arg; + } + return temp; + + } + + + + } #endif From 3c112a7a25eebf5d1ae611f900d0d65057336fb3 Mon Sep 17 00:00:00 2001 From: Guido Cossu Date: Thu, 25 May 2017 12:09:00 +0100 Subject: [PATCH 42/53] Small correction to the general exp definition --- lib/tensors/Tensor_exp.h | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/lib/tensors/Tensor_exp.h b/lib/tensors/Tensor_exp.h index 64c06bb5..2f584134 100644 --- a/lib/tensors/Tensor_exp.h +++ b/lib/tensors/Tensor_exp.h @@ -127,7 +127,7 @@ template::Tens mat unit(1.0); mat temp(unit); for(int i=Nexp; i>=1;--i){ - temp *= alpha/ComplexD(i); + temp *= alpha/ReadD(i); temp = unit + temp*arg; } return temp; From 75856f29454eaa3cfd7f002fa90ebd074ba663f4 Mon Sep 17 00:00:00 2001 From: Guido Cossu Date: Thu, 25 May 2017 12:44:56 +0100 Subject: [PATCH 43/53] Compilation fix in the Tensor_exp --- lib/tensors/Tensor_exp.h | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/lib/tensors/Tensor_exp.h b/lib/tensors/Tensor_exp.h index 2f584134..e18fed70 100644 --- a/lib/tensors/Tensor_exp.h +++ b/lib/tensors/Tensor_exp.h @@ -127,7 +127,7 @@ template::Tens mat unit(1.0); mat temp(unit); for(int i=Nexp; i>=1;--i){ - temp *= alpha/ReadD(i); + temp *= alpha/RealD(i); temp = unit + temp*arg; } return temp; From a74c34315c79c4b52871e87b312943174dbaf1e2 Mon Sep 17 00:00:00 2001 From: Antonin Portelli Date: Thu, 25 May 2017 14:27:49 +0100 Subject: [PATCH 44/53] Bootstrap script fix --- bootstrap.sh | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/bootstrap.sh b/bootstrap.sh index 66eb63ee..dfb6735d 100755 --- a/bootstrap.sh +++ b/bootstrap.sh @@ -1,4 +1,4 @@ -]#!/usr/bin/env bash +#!/usr/bin/env bash EIGEN_URL='http://bitbucket.org/eigen/eigen/get/3.3.3.tar.bz2' From f4e8bf28589d2629364b7e4a3279d76879ef818a Mon Sep 17 00:00:00 2001 From: Guido Cossu Date: Fri, 26 May 2017 12:45:59 +0100 Subject: [PATCH 45/53] Fixing the topological charge. Wilson Flow tested, ok --- lib/qcd/smearing/WilsonFlow.h | 18 +++++++++++++++--- lib/qcd/utils/WilsonLoops.h | 26 ++++++++++++++++++++------ tests/smearing/Test_WilsonFlow.cc | 11 +++++++++-- 3 files changed, 44 insertions(+), 11 deletions(-) diff --git a/lib/qcd/smearing/WilsonFlow.h b/lib/qcd/smearing/WilsonFlow.h index 8b8f9a81..34800d11 100644 --- a/lib/qcd/smearing/WilsonFlow.h +++ b/lib/qcd/smearing/WilsonFlow.h @@ -36,8 +36,10 @@ namespace QCD { template class WilsonFlow: public Smear{ unsigned int Nstep; + unsigned int measure_interval; RealD epsilon; + mutable WilsonGaugeAction SG; void evolve_step(typename Gimpl::GaugeField&) const; @@ -47,9 +49,10 @@ class WilsonFlow: public Smear{ public: INHERIT_GIMPL_TYPES(Gimpl) - explicit WilsonFlow(unsigned int Nstep, RealD epsilon): + explicit WilsonFlow(unsigned int Nstep, RealD epsilon, unsigned int interval = 1): Nstep(Nstep), epsilon(epsilon), + measure_interval(interval), SG(WilsonGaugeAction(3.0)) { // WilsonGaugeAction with beta 3.0 assert(epsilon > 0.0); @@ -107,11 +110,20 @@ RealD WilsonFlow::energyDensityPlaquette(unsigned int step, const GaugeFi template void WilsonFlow::smear(GaugeField& out, const GaugeField& in) const { out = in; - for (unsigned int step = 0; step < Nstep; step++) { + for (unsigned int step = 1; step <= Nstep; step++) { + auto start = std::chrono::high_resolution_clock::now(); evolve_step(out); + auto end = std::chrono::high_resolution_clock::now(); + std::chrono::duration diff = end - start; + std::cout << "Time to evolve " << diff.count() << " s\n"; std::cout << GridLogMessage << "[WilsonFlow] Energy density (plaq) : " - << step << " " + << step << " " << energyDensityPlaquette(step,out) << std::endl; + if( step % measure_interval == 0){ + std::cout << GridLogMessage << "[WilsonFlow] Top. charge : " + << step << " " + << WilsonLoops::TopologicalCharge(out) << std::endl; + } } } diff --git a/lib/qcd/utils/WilsonLoops.h b/lib/qcd/utils/WilsonLoops.h index 5382882e..77b5db95 100644 --- a/lib/qcd/utils/WilsonLoops.h +++ b/lib/qcd/utils/WilsonLoops.h @@ -197,10 +197,13 @@ public: std::vector U(Nd, grid); for (int d = 0; d < Nd; d++) { + // this operation is taking too much time U[d] = PeekIndex(Umu, d); } staple = zero; - GaugeMat tmp(grid); + GaugeMat tmp1(grid); + GaugeMat tmp2(grid); + for (int nu = 0; nu < Nd; nu++) { @@ -214,22 +217,34 @@ public: // | // __| // + tmp1 = Cshift(U[nu], mu, 1); + tmp2 = Cshift(U[mu], nu, 1); + staple += tmp1*adj(U[nu] * tmp2); + + +/* staple += Gimpl::ShiftStaple( Gimpl::CovShiftForward( U[nu], nu, Gimpl::CovShiftBackward( U[mu], mu, Gimpl::CovShiftIdentityBackward(U[nu], nu))), mu); - +*/ // __ // | // |__ // // + + tmp2 = adj(U[mu]*tmp1)*U[nu]; + staple += Cshift(tmp2, nu, -1); + + /* staple += Gimpl::ShiftStaple( Gimpl::CovShiftBackward(U[nu], nu, Gimpl::CovShiftBackward(U[mu], mu, U[nu])), mu); +*/ } } } @@ -289,8 +304,7 @@ public: // staple = Gimpl::ShiftStaple( Gimpl::CovShiftBackward(U[nu], nu, - Gimpl::CovShiftBackward(U[mu], mu, U[nu])), - mu); + Gimpl::CovShiftBackward(U[mu], mu, U[nu])), mu); } } @@ -307,10 +321,10 @@ public: GaugeMat Vup(Umu._grid), Vdn(Umu._grid); StapleUpper(Vup, Umu, mu, nu); StapleLower(Vdn, Umu, mu, nu); - GaugeMat v = adj(Vup) - adj(Vdn); + GaugeMat v = Vup - Vdn; GaugeMat u = PeekIndex(Umu, mu); // some redundant copies GaugeMat vu = v*u; - FS = 0.25*Ta(u*v + Cshift(vu, mu, +1)); + FS = 0.25*Ta(u*v + Cshift(vu, mu, -1)); } static Real TopologicalCharge(GaugeLorentz &U){ diff --git a/tests/smearing/Test_WilsonFlow.cc b/tests/smearing/Test_WilsonFlow.cc index 4e6bd0af..19ae5575 100644 --- a/tests/smearing/Test_WilsonFlow.cc +++ b/tests/smearing/Test_WilsonFlow.cc @@ -42,22 +42,29 @@ int main(int argc, char **argv) { GridRedBlackCartesian RBGrid(latt_size, simd_layout, mpi_layout); std::vector seeds({1, 2, 3, 4, 5}); + GridSerialRNG sRNG; GridParallelRNG pRNG(&Grid); pRNG.SeedFixedIntegers(seeds); LatticeGaugeField Umu(&Grid), Uflow(&Grid); SU::HotConfiguration(pRNG, Umu); + CheckpointerParameters CPPar("ckpoint_lat", "ckpoint_rng"); + BinaryHmcCheckpointer CPBin(CPPar); + + CPBin.CheckpointRestore(3000, Umu, sRNG, pRNG); std::cout << std::setprecision(15); std::cout << GridLogMessage << "Plaquette: " << WilsonLoops::avgPlaquette(Umu) << std::endl; - WilsonFlow WF(200, 0.01); + WilsonFlow WF(200, 0.01, 50); WF.smear(Uflow, Umu); RealD WFlow_plaq = WilsonLoops::avgPlaquette(Uflow); - std::cout << GridLogMessage << "Plaquette: "<< WFlow_plaq << std::endl; + RealD WFlow_TC = WilsonLoops::TopologicalCharge(Uflow); + std::cout << GridLogMessage << "Plaquette : "<< WFlow_plaq << std::endl; + std::cout << GridLogMessage << "TopologicalCharge : "<< WFlow_TC << std::endl; std::cout<< GridLogMessage << " Admissibility check:\n"; const double sp_adm = 0.067; // admissible threshold From 0de314870d5d926fb71d25e87dfdd0ba38f89742 Mon Sep 17 00:00:00 2001 From: Guido Cossu Date: Fri, 26 May 2017 14:31:49 +0100 Subject: [PATCH 46/53] Faster derivative for WilsonGauge --- lib/qcd/action/gauge/WilsonGaugeAction.h | 12 ++++-- lib/qcd/utils/WilsonLoops.h | 47 +++++++++++++++--------- 2 files changed, 37 insertions(+), 22 deletions(-) diff --git a/lib/qcd/action/gauge/WilsonGaugeAction.h b/lib/qcd/action/gauge/WilsonGaugeAction.h index 77c2424c..1ea780b7 100644 --- a/lib/qcd/action/gauge/WilsonGaugeAction.h +++ b/lib/qcd/action/gauge/WilsonGaugeAction.h @@ -71,14 +71,18 @@ class WilsonGaugeAction : public Action { RealD factor = 0.5 * beta / RealD(Nc); - GaugeLinkField Umu(U._grid); + //GaugeLinkField Umu(U._grid); GaugeLinkField dSdU_mu(U._grid); for (int mu = 0; mu < Nd; mu++) { - Umu = PeekIndex(U, mu); + //Umu = PeekIndex(U, mu); // Staple in direction mu - WilsonLoops::Staple(dSdU_mu, U, mu); - dSdU_mu = Ta(Umu * dSdU_mu) * factor; + //WilsonLoops::Staple(dSdU_mu, U, mu); + //dSdU_mu = Ta(Umu * dSdU_mu) * factor; + + + WilsonLoops::StapleMult(dSdU_mu, U, mu); + dSdU_mu = Ta(dSdU_mu) * factor; PokeIndex(dSdU, dSdU_mu, mu); } diff --git a/lib/qcd/utils/WilsonLoops.h b/lib/qcd/utils/WilsonLoops.h index 77b5db95..3c0e9ee4 100644 --- a/lib/qcd/utils/WilsonLoops.h +++ b/lib/qcd/utils/WilsonLoops.h @@ -188,13 +188,10 @@ public: } } - ////////////////////////////////////////////////// - // the sum over all staples on each site - ////////////////////////////////////////////////// - static void Staple(GaugeMat &staple, const GaugeLorentz &Umu, int mu) { +// For the force term +static void StapleMult(GaugeMat &staple, const GaugeLorentz &Umu, int mu) { GridBase *grid = Umu._grid; - std::vector U(Nd, grid); for (int d = 0; d < Nd; d++) { // this operation is taking too much time @@ -204,6 +201,31 @@ public: GaugeMat tmp1(grid); GaugeMat tmp2(grid); + for (int nu = 0; nu < Nd; nu++) { + if (nu != mu) { + // this is ~10% faster than the Staple + tmp1 = Cshift(U[nu], mu, 1); + tmp2 = Cshift(U[mu], nu, 1); + staple += tmp1* adj(U[nu]*tmp2); + tmp2 = adj(U[mu]*tmp1)*U[nu]; + staple += Cshift(tmp2, nu, -1); + } + } + staple = U[mu]*staple; +} + + ////////////////////////////////////////////////// + // the sum over all staples on each site + ////////////////////////////////////////////////// + static void Staple(GaugeMat &staple, const GaugeLorentz &Umu, int mu) { + + GridBase *grid = Umu._grid; + + std::vector U(Nd, grid); + for (int d = 0; d < Nd; d++) { + U[d] = PeekIndex(Umu, d); + } + staple = zero; for (int nu = 0; nu < Nd; nu++) { @@ -217,34 +239,23 @@ public: // | // __| // - tmp1 = Cshift(U[nu], mu, 1); - tmp2 = Cshift(U[mu], nu, 1); - - staple += tmp1*adj(U[nu] * tmp2); - - -/* + staple += Gimpl::ShiftStaple( Gimpl::CovShiftForward( U[nu], nu, Gimpl::CovShiftBackward( U[mu], mu, Gimpl::CovShiftIdentityBackward(U[nu], nu))), mu); -*/ + // __ // | // |__ // // - tmp2 = adj(U[mu]*tmp1)*U[nu]; - staple += Cshift(tmp2, nu, -1); - - /* staple += Gimpl::ShiftStaple( Gimpl::CovShiftBackward(U[nu], nu, Gimpl::CovShiftBackward(U[mu], mu, U[nu])), mu); -*/ } } } From 8e0ced627a9b134e8a99bc94b9b5032e8a8141ed Mon Sep 17 00:00:00 2001 From: Lanny91 Date: Fri, 26 May 2017 15:59:15 +0100 Subject: [PATCH 47/53] Hadrons: Fermion boundary conditions can now be set in measurement code. --- extras/Hadrons/Modules/MAction/DWF.hpp | 10 ++++++++-- extras/Hadrons/Modules/MAction/Wilson.hpp | 10 ++++++++-- tests/hadrons/Test_hadrons_meson_3pt.cc | 5 +++++ tests/hadrons/Test_hadrons_rarekaon.cc | 5 +++++ tests/hadrons/Test_hadrons_spectrum.cc | 5 +++++ 5 files changed, 31 insertions(+), 4 deletions(-) diff --git a/extras/Hadrons/Modules/MAction/DWF.hpp b/extras/Hadrons/Modules/MAction/DWF.hpp index 49861e3e..880fe7b9 100644 --- a/extras/Hadrons/Modules/MAction/DWF.hpp +++ b/extras/Hadrons/Modules/MAction/DWF.hpp @@ -48,7 +48,8 @@ public: std::string, gauge, unsigned int, Ls, double , mass, - double , M5); + double , M5, + std::string , boundary); }; template @@ -116,14 +117,19 @@ void TDWF::execute(void) << par().mass << ", M5= " << par().M5 << " and Ls= " << par().Ls << " using gauge field '" << par().gauge << "'" << std::endl; + LOG(Message) << "Fermion boundary conditions: " << par().boundary + << std::endl; env().createGrid(par().Ls); auto &U = *env().template getObject(par().gauge); auto &g4 = *env().getGrid(); auto &grb4 = *env().getRbGrid(); auto &g5 = *env().getGrid(par().Ls); auto &grb5 = *env().getRbGrid(par().Ls); + std::vector boundary = strToVec(par().boundary); + typename DomainWallFermion::ImplParams implParams(boundary); FMat *fMatPt = new DomainWallFermion(U, g5, grb5, g4, grb4, - par().mass, par().M5); + par().mass, par().M5, + implParams); env().setObject(getName(), fMatPt); } diff --git a/extras/Hadrons/Modules/MAction/Wilson.hpp b/extras/Hadrons/Modules/MAction/Wilson.hpp index 6ffa997d..4b84bda5 100644 --- a/extras/Hadrons/Modules/MAction/Wilson.hpp +++ b/extras/Hadrons/Modules/MAction/Wilson.hpp @@ -46,7 +46,8 @@ class WilsonPar: Serializable public: GRID_SERIALIZABLE_CLASS_MEMBERS(WilsonPar, std::string, gauge, - double , mass); + double , mass, + std::string, boundary); }; template @@ -112,10 +113,15 @@ void TWilson::execute() { LOG(Message) << "Setting up TWilson fermion matrix with m= " << par().mass << " using gauge field '" << par().gauge << "'" << std::endl; + LOG(Message) << "Fermion boundary conditions: " << par().boundary + << std::endl; auto &U = *env().template getObject(par().gauge); auto &grid = *env().getGrid(); auto &gridRb = *env().getRbGrid(); - FMat *fMatPt = new WilsonFermion(U, grid, gridRb, par().mass); + std::vector boundary = strToVec(par().boundary); + typename WilsonFermion::ImplParams implParams(boundary); + FMat *fMatPt = new WilsonFermion(U, grid, gridRb, par().mass, + implParams); env().setObject(getName(), fMatPt); } diff --git a/tests/hadrons/Test_hadrons_meson_3pt.cc b/tests/hadrons/Test_hadrons_meson_3pt.cc index efef6931..7e487153 100644 --- a/tests/hadrons/Test_hadrons_meson_3pt.cc +++ b/tests/hadrons/Test_hadrons_meson_3pt.cc @@ -61,6 +61,10 @@ int main(int argc, char *argv[]) // gauge field application.createModule("gauge"); + + // set fermion boundary conditions to be periodic space, antiperiodic time. + std::string boundary = "1 1 1 -1"; + for (unsigned int i = 0; i < flavour.size(); ++i) { // actions @@ -69,6 +73,7 @@ int main(int argc, char *argv[]) actionPar.Ls = 12; actionPar.M5 = 1.8; actionPar.mass = mass[i]; + actionPar.boundary = boundary; application.createModule("DWF_" + flavour[i], actionPar); // solvers diff --git a/tests/hadrons/Test_hadrons_rarekaon.cc b/tests/hadrons/Test_hadrons_rarekaon.cc index 89d7d501..ab4d3ef1 100644 --- a/tests/hadrons/Test_hadrons_rarekaon.cc +++ b/tests/hadrons/Test_hadrons_rarekaon.cc @@ -98,6 +98,10 @@ int main(int argc, char *argv[]) gaugePar.file = configStem; application.createModule("gauge", gaugePar); } + + // set fermion boundary conditions to be periodic space, antiperiodic time. + std::string boundary = "1 1 1 -1"; + for (unsigned int i = 0; i < flavour.size(); ++i) { // actions @@ -106,6 +110,7 @@ int main(int argc, char *argv[]) actionPar.Ls = 16; actionPar.M5 = 1.8; actionPar.mass = mass[i]; + actionPar.boundary = boundary; application.createModule("DWF_" + flavour[i], actionPar); // solvers diff --git a/tests/hadrons/Test_hadrons_spectrum.cc b/tests/hadrons/Test_hadrons_spectrum.cc index 2d731ff4..55f3346e 100644 --- a/tests/hadrons/Test_hadrons_spectrum.cc +++ b/tests/hadrons/Test_hadrons_spectrum.cc @@ -63,6 +63,10 @@ int main(int argc, char *argv[]) MSource::Point::Par ptPar; ptPar.position = "0 0 0 0"; application.createModule("pt", ptPar); + + // set fermion boundary conditions to be periodic space, antiperiodic time. + std::string boundary = "1 1 1 -1"; + for (unsigned int i = 0; i < flavour.size(); ++i) { // actions @@ -71,6 +75,7 @@ int main(int argc, char *argv[]) actionPar.Ls = 12; actionPar.M5 = 1.8; actionPar.mass = mass[i]; + actionPar.boundary = boundary; application.createModule("DWF_" + flavour[i], actionPar); // solvers From 7c6cc85df6b6f6a43f86099f033aa762777ede55 Mon Sep 17 00:00:00 2001 From: Guido Cossu Date: Sat, 27 May 2017 18:03:49 +0100 Subject: [PATCH 48/53] Updating WilsonFlow test --- lib/qcd/smearing/WilsonFlow.h | 5 +++ tests/smearing/Test_WilsonFlow.cc | 52 ++++++++++++++++++++++++++----- 2 files changed, 50 insertions(+), 7 deletions(-) diff --git a/lib/qcd/smearing/WilsonFlow.h b/lib/qcd/smearing/WilsonFlow.h index 34800d11..00b7622e 100644 --- a/lib/qcd/smearing/WilsonFlow.h +++ b/lib/qcd/smearing/WilsonFlow.h @@ -107,6 +107,9 @@ RealD WilsonFlow::energyDensityPlaquette(unsigned int step, const GaugeFi return 2.0 * td * td * SG.S(U)/U._grid->gSites(); } + +//#define WF_TIMING + template void WilsonFlow::smear(GaugeField& out, const GaugeField& in) const { out = in; @@ -115,7 +118,9 @@ void WilsonFlow::smear(GaugeField& out, const GaugeField& in) const { evolve_step(out); auto end = std::chrono::high_resolution_clock::now(); std::chrono::duration diff = end - start; + #ifdef WF_TIMING std::cout << "Time to evolve " << diff.count() << " s\n"; + #endif std::cout << GridLogMessage << "[WilsonFlow] Energy density (plaq) : " << step << " " << energyDensityPlaquette(step,out) << std::endl; diff --git a/tests/smearing/Test_WilsonFlow.cc b/tests/smearing/Test_WilsonFlow.cc index 19ae5575..c337bf29 100644 --- a/tests/smearing/Test_WilsonFlow.cc +++ b/tests/smearing/Test_WilsonFlow.cc @@ -28,6 +28,37 @@ directory /* END LEGAL */ #include +namespace Grid{ + struct WFParameters: Serializable { + GRID_SERIALIZABLE_CLASS_MEMBERS(WFParameters, + int, steps, + double, step_size, + int, meas_interval); + + + template + WFParameters(Reader& Reader){ + read(Reader, "WilsonFlow", *this); + } + + }; + + struct ConfParameters: Serializable { + GRID_SERIALIZABLE_CLASS_MEMBERS(ConfParameters, + std::string, conf_prefix, + std::string, rng_prefix, + int, StartConfiguration, + int, EndConfiguration, + int, Skip); + + template + ConfParameters(Reader& Reader){ + read(Reader, "Configurations", *this); + } + + }; +} + int main(int argc, char **argv) { using namespace Grid; using namespace Grid::QCD; @@ -48,23 +79,30 @@ int main(int argc, char **argv) { LatticeGaugeField Umu(&Grid), Uflow(&Grid); SU::HotConfiguration(pRNG, Umu); - CheckpointerParameters CPPar("ckpoint_lat", "ckpoint_rng"); + + typedef Grid::JSONReader Serialiser; + Serialiser Reader("input.json"); + WFParameters WFPar(Reader); + ConfParameters CPar(Reader); + CheckpointerParameters CPPar(CPar.conf_prefix, CPar.rng_prefix); BinaryHmcCheckpointer CPBin(CPPar); - CPBin.CheckpointRestore(3000, Umu, sRNG, pRNG); + for (int conf = CPar.StartConfiguration; conf <= CPar.EndConfiguration; conf+= CPar.Skip){ + + CPBin.CheckpointRestore(conf, Umu, sRNG, pRNG); std::cout << std::setprecision(15); - std::cout << GridLogMessage << "Plaquette: " + std::cout << GridLogMessage << "Initial plaquette: " << WilsonLoops::avgPlaquette(Umu) << std::endl; - WilsonFlow WF(200, 0.01, 50); + WilsonFlow WF(WFPar.steps, WFPar.step_size, WFPar.meas_interval); WF.smear(Uflow, Umu); RealD WFlow_plaq = WilsonLoops::avgPlaquette(Uflow); RealD WFlow_TC = WilsonLoops::TopologicalCharge(Uflow); - std::cout << GridLogMessage << "Plaquette : "<< WFlow_plaq << std::endl; - std::cout << GridLogMessage << "TopologicalCharge : "<< WFlow_TC << std::endl; + std::cout << GridLogMessage << "Plaquette "<< conf << " " << WFlow_plaq << std::endl; + std::cout << GridLogMessage << "TopologicalCharge "<< conf << " " << WFlow_TC << std::endl; std::cout<< GridLogMessage << " Admissibility check:\n"; const double sp_adm = 0.067; // admissible threshold @@ -80,6 +118,6 @@ int main(int argc, char **argv) { std::cout<< GridLogMessage << " (sp_admissible = "<< sp_adm <<")\n"; //std::cout<< GridLogMessage << " sp_admissible - sp_max = "< Date: Mon, 29 May 2017 12:57:33 +0100 Subject: [PATCH 49/53] Hadrons: mesons gamma list fix --- extras/Hadrons/Modules/MContraction/Meson.hpp | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/extras/Hadrons/Modules/MContraction/Meson.hpp b/extras/Hadrons/Modules/MContraction/Meson.hpp index 4cbe1ac4..09c2a6e1 100644 --- a/extras/Hadrons/Modules/MContraction/Meson.hpp +++ b/extras/Hadrons/Modules/MContraction/Meson.hpp @@ -131,12 +131,11 @@ std::vector TMeson::getOutput(void) template void TMeson::parseGammaString(std::vector &gammaList) { + gammaList.clear(); // Determine gamma matrices to insert at source/sink. if (par().gammas.compare("all") == 0) { // Do all contractions. - unsigned int n_gam = Ns * Ns; - gammaList.resize(n_gam*n_gam); for (unsigned int i = 1; i < Gamma::nGamma; i += 2) { for (unsigned int j = 1; j < Gamma::nGamma; j += 2) From aaf1e33a77a31ba0d06e69d1ed476722ce015f5c Mon Sep 17 00:00:00 2001 From: Guido Cossu Date: Fri, 2 Jun 2017 16:32:35 +0100 Subject: [PATCH 50/53] Adding adaptive integration in the WilsonFlow --- lib/qcd/smearing/WilsonFlow.h | 53 +++++++++++++++++++++++++++++-- tests/smearing/Test_WilsonFlow.cc | 27 ++++++++++++++++ 2 files changed, 78 insertions(+), 2 deletions(-) diff --git a/lib/qcd/smearing/WilsonFlow.h b/lib/qcd/smearing/WilsonFlow.h index 00b7622e..1053775d 100644 --- a/lib/qcd/smearing/WilsonFlow.h +++ b/lib/qcd/smearing/WilsonFlow.h @@ -37,15 +37,15 @@ template class WilsonFlow: public Smear{ unsigned int Nstep; unsigned int measure_interval; - RealD epsilon; + mutable RealD epsilon, taus; mutable WilsonGaugeAction SG; void evolve_step(typename Gimpl::GaugeField&) const; + void evolve_step_adaptive(typename Gimpl::GaugeField&) const; RealD tau(unsigned int t)const {return epsilon*(t+1.0); } - public: INHERIT_GIMPL_TYPES(Gimpl) @@ -101,20 +101,63 @@ void WilsonFlow::evolve_step(typename Gimpl::GaugeField &U) const{ Gimpl::update_field(Z, U, -2.0*epsilon); // V(t+e) = exp(ep*Z)*W2 } +template +void WilsonFlow::evolve_step_adaptive(typename Gimpl::GaugeField &U) const { + GaugeField Z(U._grid); + GaugeField Zprime(U._grid); + GaugeField tmp(U._grid), Uprime(U._grid); + Uprime = U; + SG.deriv(U, Z); + Zprime = -Z; + Z *= 0.25; // Z0 = 1/4 * F(U) + Gimpl::update_field(Z, U, -2.0*epsilon); // U = W1 = exp(ep*Z0)*W0 + + Z *= -17.0/8.0; + SG.deriv(U, tmp); Z += tmp; // -17/32*Z0 +Z1 + Zprime += 2.0*tmp; + Z *= 8.0/9.0; // Z = -17/36*Z0 +8/9*Z1 + Gimpl::update_field(Z, U, -2.0*epsilon); // U_= W2 = exp(ep*Z)*W1 + + + Z *= -4.0/3.0; + SG.deriv(U, tmp); Z += tmp; // 4/3*(17/36*Z0 -8/9*Z1) +Z2 + Z *= 3.0/4.0; // Z = 17/36*Z0 -8/9*Z1 +3/4*Z2 + Gimpl::update_field(Z, U, -2.0*epsilon); // V(t+e) = exp(ep*Z)*W2 + + // Ramos + Gimpl::update_field(Zprime, Uprime, -2.0*epsilon); // V'(t+e) = exp(ep*Z')*W0 + // Compute distance as norm^2 of the difference + GaugeField diffU = U - Uprime; + RealD diff = norm2(diffU); + // adjust integration step + + taus += epsilon; + std::cout << GridLogMessage << "Adjusting integration step with distance: " << diff << std::endl; + + epsilon = epsilon*0.95*std::pow(1e-4/diff,1./3.); + std::cout << GridLogMessage << "New epsilon : " << epsilon << std::endl; + +} + template RealD WilsonFlow::energyDensityPlaquette(unsigned int step, const GaugeField& U) const { RealD td = tau(step); + //RealD td = tau; return 2.0 * td * td * SG.S(U)/U._grid->gSites(); } //#define WF_TIMING + + template void WilsonFlow::smear(GaugeField& out, const GaugeField& in) const { out = in; + taus = epsilon; for (unsigned int step = 1; step <= Nstep; step++) { auto start = std::chrono::high_resolution_clock::now(); + std::cout << GridLogMessage << "Evolution time :"<< tau(step) << std::endl; evolve_step(out); auto end = std::chrono::high_resolution_clock::now(); std::chrono::duration diff = end - start; @@ -130,6 +173,12 @@ void WilsonFlow::smear(GaugeField& out, const GaugeField& in) const { << WilsonLoops::TopologicalCharge(out) << std::endl; } } + std::cout << GridLogMessage << "[WilsonFlow] Energy density (plaq) : " + << step << " " + << energyDensityPlaquette(Nstep,out) << std::endl; + std::cout << GridLogMessage << "[WilsonFlow] Top. charge : " + << step << " " + << WilsonLoops::TopologicalCharge(out) << std::endl; } } // namespace QCD diff --git a/tests/smearing/Test_WilsonFlow.cc b/tests/smearing/Test_WilsonFlow.cc index c337bf29..77e7c5d1 100644 --- a/tests/smearing/Test_WilsonFlow.cc +++ b/tests/smearing/Test_WilsonFlow.cc @@ -101,7 +101,9 @@ int main(int argc, char **argv) { RealD WFlow_plaq = WilsonLoops::avgPlaquette(Uflow); RealD WFlow_TC = WilsonLoops::TopologicalCharge(Uflow); + RealD WFlow_T0 = WF.energyDensityPlaquette(WFPar.steps, Uflow); std::cout << GridLogMessage << "Plaquette "<< conf << " " << WFlow_plaq << std::endl; + std::cout << GridLogMessage << "T0 "<< conf << " " << WFlow_T0 << std::endl; std::cout << GridLogMessage << "TopologicalCharge "<< conf << " " << WFlow_TC << std::endl; std::cout<< GridLogMessage << " Admissibility check:\n"; @@ -121,3 +123,28 @@ int main(int argc, char **argv) { } Grid_finalize(); } // main + + +/* +Input file example + + +JSON + +{ + "WilsonFlow":{ + "steps": 200, + "step_size": 0.01, + "meas_interval": 50 + }, + "Configurations":{ + "conf_prefix": "ckpoint_lat", + "rng_prefix": "ckpoint_rng", + "StartConfiguration": 3000, + "EndConfiguration": 3000, + "Skip": 5 + } +} + + +*/ \ No newline at end of file From 7da4856e8e90131870d29ab61bd16a7660e78d1c Mon Sep 17 00:00:00 2001 From: Guido Cossu Date: Fri, 2 Jun 2017 16:55:53 +0100 Subject: [PATCH 51/53] Wilson flow with adaptive steps --- lib/qcd/smearing/WilsonFlow.h | 47 +++++++++++++++++++++++++++-------- 1 file changed, 37 insertions(+), 10 deletions(-) diff --git a/lib/qcd/smearing/WilsonFlow.h b/lib/qcd/smearing/WilsonFlow.h index 1053775d..5e9f2d95 100644 --- a/lib/qcd/smearing/WilsonFlow.h +++ b/lib/qcd/smearing/WilsonFlow.h @@ -43,7 +43,7 @@ class WilsonFlow: public Smear{ mutable WilsonGaugeAction SG; void evolve_step(typename Gimpl::GaugeField&) const; - void evolve_step_adaptive(typename Gimpl::GaugeField&) const; + void evolve_step_adaptive(typename Gimpl::GaugeField&, RealD); RealD tau(unsigned int t)const {return epsilon*(t+1.0); } public: @@ -75,7 +75,9 @@ class WilsonFlow: public Smear{ // undefined for WilsonFlow } + void smear_adaptive(GaugeField&, const GaugeField&, RealD maxTau); RealD energyDensityPlaquette(unsigned int step, const GaugeField& U) const; + RealD energyDensityPlaquette(const GaugeField& U) const; }; @@ -102,7 +104,11 @@ void WilsonFlow::evolve_step(typename Gimpl::GaugeField &U) const{ } template -void WilsonFlow::evolve_step_adaptive(typename Gimpl::GaugeField &U) const { +void WilsonFlow::evolve_step_adaptive(typename Gimpl::GaugeField &U, RealD maxTau) { + if (maxTau - taus < epsilon){ + epsilon = maxTau-taus; + } + std::cout << GridLogMessage << "Integration epsilon : " << epsilon << std::endl; GaugeField Z(U._grid); GaugeField Zprime(U._grid); GaugeField tmp(U._grid), Uprime(U._grid); @@ -142,10 +148,14 @@ void WilsonFlow::evolve_step_adaptive(typename Gimpl::GaugeField &U) cons template RealD WilsonFlow::energyDensityPlaquette(unsigned int step, const GaugeField& U) const { RealD td = tau(step); - //RealD td = tau; return 2.0 * td * td * SG.S(U)/U._grid->gSites(); } +template +RealD WilsonFlow::energyDensityPlaquette(const GaugeField& U) const { + return 2.0 * taus * taus * SG.S(U)/U._grid->gSites(); +} + //#define WF_TIMING @@ -154,7 +164,6 @@ RealD WilsonFlow::energyDensityPlaquette(unsigned int step, const GaugeFi template void WilsonFlow::smear(GaugeField& out, const GaugeField& in) const { out = in; - taus = epsilon; for (unsigned int step = 1; step <= Nstep; step++) { auto start = std::chrono::high_resolution_clock::now(); std::cout << GridLogMessage << "Evolution time :"<< tau(step) << std::endl; @@ -173,14 +182,32 @@ void WilsonFlow::smear(GaugeField& out, const GaugeField& in) const { << WilsonLoops::TopologicalCharge(out) << std::endl; } } - std::cout << GridLogMessage << "[WilsonFlow] Energy density (plaq) : " - << step << " " - << energyDensityPlaquette(Nstep,out) << std::endl; - std::cout << GridLogMessage << "[WilsonFlow] Top. charge : " - << step << " " - << WilsonLoops::TopologicalCharge(out) << std::endl; } +template +void WilsonFlow::smear_adaptive(GaugeField& out, const GaugeField& in, RealD maxTau){ + out = in; + taus = epsilon; + unsigned int step = 0; + do{ + step++; + std::cout << GridLogMessage << "Evolution time :"<< taus << std::endl; + evolve_step_adaptive(out, maxTau); + std::cout << GridLogMessage << "[WilsonFlow] Energy density (plaq) : " + << step << " " + << energyDensityPlaquette(out) << std::endl; + if( step % measure_interval == 0){ + std::cout << GridLogMessage << "[WilsonFlow] Top. charge : " + << step << " " + << WilsonLoops::TopologicalCharge(out) << std::endl; + } + } while (taus < maxTau); + + + +} + + } // namespace QCD } // namespace Grid From 4a8c4ccfba1d05159348d21a9698028ea847e77b Mon Sep 17 00:00:00 2001 From: Guido Cossu Date: Fri, 2 Jun 2017 17:02:29 +0100 Subject: [PATCH 52/53] Test wilson flow, added maxTau for adaptive flow --- tests/smearing/Test_WilsonFlow.cc | 10 ++++++---- 1 file changed, 6 insertions(+), 4 deletions(-) diff --git a/tests/smearing/Test_WilsonFlow.cc b/tests/smearing/Test_WilsonFlow.cc index 77e7c5d1..5db00d5d 100644 --- a/tests/smearing/Test_WilsonFlow.cc +++ b/tests/smearing/Test_WilsonFlow.cc @@ -33,7 +33,8 @@ namespace Grid{ GRID_SERIALIZABLE_CLASS_MEMBERS(WFParameters, int, steps, double, step_size, - int, meas_interval); + int, meas_interval, + double, maxTau); // for the adaptive algorithm template @@ -97,11 +98,11 @@ int main(int argc, char **argv) { WilsonFlow WF(WFPar.steps, WFPar.step_size, WFPar.meas_interval); - WF.smear(Uflow, Umu); + WF.smear_adaptive(Uflow, Umu, WFPar.maxTau); RealD WFlow_plaq = WilsonLoops::avgPlaquette(Uflow); RealD WFlow_TC = WilsonLoops::TopologicalCharge(Uflow); - RealD WFlow_T0 = WF.energyDensityPlaquette(WFPar.steps, Uflow); + RealD WFlow_T0 = WF.energyDensityPlaquette(Uflow); std::cout << GridLogMessage << "Plaquette "<< conf << " " << WFlow_plaq << std::endl; std::cout << GridLogMessage << "T0 "<< conf << " " << WFlow_T0 << std::endl; std::cout << GridLogMessage << "TopologicalCharge "<< conf << " " << WFlow_TC << std::endl; @@ -135,7 +136,8 @@ JSON "WilsonFlow":{ "steps": 200, "step_size": 0.01, - "meas_interval": 50 + "meas_interval": 50, + "maxTau": 2.0 }, "Configurations":{ "conf_prefix": "ckpoint_lat", From 22749699a30da633f58a4d47642721c639048f31 Mon Sep 17 00:00:00 2001 From: Antonin Portelli Date: Tue, 6 Jun 2017 11:45:30 -0500 Subject: [PATCH 53/53] Fixes after merge and point sink module --- extras/Hadrons/Environment.cc | 56 ++++++++- extras/Hadrons/Environment.hpp | 44 ++++++- extras/Hadrons/Global.hpp | 9 +- extras/Hadrons/Modules.hpp | 1 + extras/Hadrons/Modules/MAction/DWF.hpp | 8 +- extras/Hadrons/Modules/MAction/Wilson.hpp | 6 +- .../Hadrons/Modules/MContraction/Baryon.hpp | 14 +-- .../Hadrons/Modules/MContraction/DiscLoop.hpp | 8 +- .../Hadrons/Modules/MContraction/Gamma3pt.hpp | 12 +- extras/Hadrons/Modules/MContraction/Meson.hpp | 102 ++++++++++------ .../Modules/MContraction/WeakHamiltonian.hpp | 8 +- .../MContraction/WeakHamiltonianEye.hpp | 6 +- .../MContraction/WeakHamiltonianNonEye.hpp | 6 +- .../MContraction/WeakNeutral4ptDisc.hpp | 6 +- extras/Hadrons/Modules/MGauge/Load.hpp | 6 +- extras/Hadrons/Modules/MGauge/Random.hpp | 6 +- extras/Hadrons/Modules/MGauge/StochEm.hpp | 6 +- extras/Hadrons/Modules/MGauge/Unit.hpp | 6 +- extras/Hadrons/Modules/MLoop/NoiseLoop.hpp | 8 +- .../Hadrons/Modules/MScalar/ChargedProp.hpp | 6 +- extras/Hadrons/Modules/MScalar/FreeProp.hpp | 6 +- extras/Hadrons/Modules/MSink/Point.hpp | 114 ++++++++++++++++++ extras/Hadrons/Modules/MSolver/RBPrecCG.hpp | 8 +- extras/Hadrons/Modules/MSource/Point.hpp | 6 +- extras/Hadrons/Modules/MSource/SeqGamma.hpp | 8 +- extras/Hadrons/Modules/MSource/Wall.hpp | 8 +- extras/Hadrons/Modules/MSource/Z2.hpp | 6 +- extras/Hadrons/Modules/Quark.hpp | 2 +- .../templates/Module_in_NS.hpp.template | 6 +- .../templates/Module_tmp_in_NS.hpp.template | 6 +- extras/Hadrons/modules.inc | 1 + tests/hadrons/Test_hadrons_spectrum.cc | 24 ++-- 32 files changed, 385 insertions(+), 134 deletions(-) create mode 100644 extras/Hadrons/Modules/MSink/Point.hpp diff --git a/extras/Hadrons/Environment.cc b/extras/Hadrons/Environment.cc index 37f2a3d7..0e7a4326 100644 --- a/extras/Hadrons/Environment.cc +++ b/extras/Hadrons/Environment.cc @@ -41,9 +41,10 @@ using namespace Hadrons; // constructor ///////////////////////////////////////////////////////////////// Environment::Environment(void) { - nd_ = GridDefaultLatt().size(); + dim_ = GridDefaultLatt(); + nd_ = dim_.size(); grid4d_.reset(SpaceTimeGrid::makeFourDimGrid( - GridDefaultLatt(), GridDefaultSimd(nd_, vComplex::Nsimd()), + dim_, GridDefaultSimd(nd_, vComplex::Nsimd()), GridDefaultMpi())); gridRb4d_.reset(SpaceTimeGrid::makeFourDimRedBlackGrid(grid4d_.get())); auto loc = getGrid()->LocalDimensions(); @@ -132,6 +133,16 @@ unsigned int Environment::getNd(void) const return nd_; } +std::vector Environment::getDim(void) const +{ + return dim_; +} + +int Environment::getDim(const unsigned int mu) const +{ + return dim_[mu]; +} + // random number generator ///////////////////////////////////////////////////// void Environment::setSeed(const std::vector &seed) { @@ -271,6 +282,21 @@ std::string Environment::getModuleType(const std::string name) const return getModuleType(getModuleAddress(name)); } +std::string Environment::getModuleNamespace(const unsigned int address) const +{ + std::string type = getModuleType(address), ns; + + auto pos2 = type.rfind("::"); + auto pos1 = type.rfind("::", pos2 - 2); + + return type.substr(pos1 + 2, pos2 - pos1 - 2); +} + +std::string Environment::getModuleNamespace(const std::string name) const +{ + return getModuleNamespace(getModuleAddress(name)); +} + bool Environment::hasModule(const unsigned int address) const { return (address < module_.size()); @@ -492,7 +518,14 @@ std::string Environment::getObjectType(const unsigned int address) const { if (hasRegisteredObject(address)) { - return typeName(object_[address].type); + if (object_[address].type) + { + return typeName(object_[address].type); + } + else + { + return ""; + } } else if (hasObject(address)) { @@ -532,6 +565,23 @@ Environment::Size Environment::getObjectSize(const std::string name) const return getObjectSize(getObjectAddress(name)); } +unsigned int Environment::getObjectModule(const unsigned int address) const +{ + if (hasObject(address)) + { + return object_[address].module; + } + else + { + HADRON_ERROR("no object with address " + std::to_string(address)); + } +} + +unsigned int Environment::getObjectModule(const std::string name) const +{ + return getObjectModule(getObjectAddress(name)); +} + unsigned int Environment::getObjectLs(const unsigned int address) const { if (hasRegisteredObject(address)) diff --git a/extras/Hadrons/Environment.hpp b/extras/Hadrons/Environment.hpp index 2628e5a0..13264bd5 100644 --- a/extras/Hadrons/Environment.hpp +++ b/extras/Hadrons/Environment.hpp @@ -106,6 +106,8 @@ public: void createGrid(const unsigned int Ls); GridCartesian * getGrid(const unsigned int Ls = 1) const; GridRedBlackCartesian * getRbGrid(const unsigned int Ls = 1) const; + std::vector getDim(void) const; + int getDim(const unsigned int mu) const; unsigned int getNd(void) const; // random number generator void setSeed(const std::vector &seed); @@ -131,6 +133,8 @@ public: std::string getModuleName(const unsigned int address) const; std::string getModuleType(const unsigned int address) const; std::string getModuleType(const std::string name) const; + std::string getModuleNamespace(const unsigned int address) const; + std::string getModuleNamespace(const std::string name) const; bool hasModule(const unsigned int address) const; bool hasModule(const std::string name) const; Graph makeModuleGraph(void) const; @@ -171,6 +175,8 @@ public: std::string getObjectType(const std::string name) const; Size getObjectSize(const unsigned int address) const; Size getObjectSize(const std::string name) const; + unsigned int getObjectModule(const unsigned int address) const; + unsigned int getObjectModule(const std::string name) const; unsigned int getObjectLs(const unsigned int address) const; unsigned int getObjectLs(const std::string name) const; bool hasObject(const unsigned int address) const; @@ -181,6 +187,10 @@ public: bool hasCreatedObject(const std::string name) const; bool isObject5d(const unsigned int address) const; bool isObject5d(const std::string name) const; + template + bool isObjectOfType(const unsigned int address) const; + template + bool isObjectOfType(const std::string name) const; Environment::Size getTotalSize(void) const; void addOwnership(const unsigned int owner, const unsigned int property); @@ -197,6 +207,7 @@ private: bool dryRun_{false}; unsigned int traj_, locVol_; // grids + std::vector dim_; GridPt grid4d_; std::map grid5d_; GridRbPt gridRb4d_; @@ -343,7 +354,7 @@ T * Environment::getObject(const unsigned int address) const else { HADRON_ERROR("object with address " + std::to_string(address) + - " does not have type '" + typeid(T).name() + + " does not have type '" + typeName(&typeid(T)) + "' (has type '" + getObjectType(address) + "')"); } } @@ -380,6 +391,37 @@ T * Environment::createLattice(const std::string name) return createLattice(getObjectAddress(name)); } +template +bool Environment::isObjectOfType(const unsigned int address) const +{ + if (hasRegisteredObject(address)) + { + if (auto h = dynamic_cast *>(object_[address].data.get())) + { + return true; + } + else + { + return false; + } + } + else if (hasObject(address)) + { + HADRON_ERROR("object with address " + std::to_string(address) + + " exists but is not registered"); + } + else + { + HADRON_ERROR("no object with address " + std::to_string(address)); + } +} + +template +bool Environment::isObjectOfType(const std::string name) const +{ + return isObjectOfType(getObjectAddress(name)); +} + END_HADRONS_NAMESPACE #endif // Hadrons_Environment_hpp_ diff --git a/extras/Hadrons/Global.hpp b/extras/Hadrons/Global.hpp index 3ff79ea3..9de01623 100644 --- a/extras/Hadrons/Global.hpp +++ b/extras/Hadrons/Global.hpp @@ -65,7 +65,9 @@ BEGIN_HADRONS_NAMESPACE typedef FermionOperator FMat##suffix; \ typedef typename FImpl::FermionField FermionField##suffix; \ typedef typename FImpl::PropagatorField PropagatorField##suffix; \ -typedef typename FImpl::SitePropagator SitePropagator##suffix; +typedef typename FImpl::SitePropagator SitePropagator##suffix; \ +typedef std::vector \ + SlicedPropagator##suffix; #define GAUGE_TYPE_ALIASES(FImpl, suffix)\ typedef typename FImpl::DoubledGaugeField DoubledGaugeField##suffix; @@ -78,7 +80,10 @@ typedef typename SImpl::Field PropagatorField##suffix; typedef std::function SolverFn##suffix; -#define TYPE_ALIASES(FImpl, suffix)\ +#define SINK_TYPE_ALIASES(suffix)\ +typedef std::function SinkFn##suffix; + +#define FGS_TYPE_ALIASES(FImpl, suffix)\ FERM_TYPE_ALIASES(FImpl, suffix)\ GAUGE_TYPE_ALIASES(FImpl, suffix)\ SOLVER_TYPE_ALIASES(FImpl, suffix) diff --git a/extras/Hadrons/Modules.hpp b/extras/Hadrons/Modules.hpp index 7155c02a..42a1f651 100644 --- a/extras/Hadrons/Modules.hpp +++ b/extras/Hadrons/Modules.hpp @@ -16,6 +16,7 @@ #include #include #include +#include #include #include #include diff --git a/extras/Hadrons/Modules/MAction/DWF.hpp b/extras/Hadrons/Modules/MAction/DWF.hpp index 880fe7b9..78e0916c 100644 --- a/extras/Hadrons/Modules/MAction/DWF.hpp +++ b/extras/Hadrons/Modules/MAction/DWF.hpp @@ -27,8 +27,8 @@ See the full license in the file "LICENSE" in the top level distribution directo *************************************************************************************/ /* END LEGAL */ -#ifndef Hadrons_DWF_hpp_ -#define Hadrons_DWF_hpp_ +#ifndef Hadrons_MAction_DWF_hpp_ +#define Hadrons_MAction_DWF_hpp_ #include #include @@ -56,7 +56,7 @@ template class TDWF: public Module { public: - TYPE_ALIASES(FImpl,); + FGS_TYPE_ALIASES(FImpl,); public: // constructor TDWF(const std::string name); @@ -137,4 +137,4 @@ END_MODULE_NAMESPACE END_HADRONS_NAMESPACE -#endif // Hadrons_DWF_hpp_ +#endif // Hadrons_MAction_DWF_hpp_ diff --git a/extras/Hadrons/Modules/MAction/Wilson.hpp b/extras/Hadrons/Modules/MAction/Wilson.hpp index 4b84bda5..aab54245 100644 --- a/extras/Hadrons/Modules/MAction/Wilson.hpp +++ b/extras/Hadrons/Modules/MAction/Wilson.hpp @@ -27,8 +27,8 @@ See the full license in the file "LICENSE" in the top level distribution directo *************************************************************************************/ /* END LEGAL */ -#ifndef Hadrons_Wilson_hpp_ -#define Hadrons_Wilson_hpp_ +#ifndef Hadrons_MAction_Wilson_hpp_ +#define Hadrons_MAction_Wilson_hpp_ #include #include @@ -54,7 +54,7 @@ template class TWilson: public Module { public: - TYPE_ALIASES(FImpl,); + FGS_TYPE_ALIASES(FImpl,); public: // constructor TWilson(const std::string name); diff --git a/extras/Hadrons/Modules/MContraction/Baryon.hpp b/extras/Hadrons/Modules/MContraction/Baryon.hpp index be7d919c..78bde5a2 100644 --- a/extras/Hadrons/Modules/MContraction/Baryon.hpp +++ b/extras/Hadrons/Modules/MContraction/Baryon.hpp @@ -27,8 +27,8 @@ See the full license in the file "LICENSE" in the top level distribution directo *************************************************************************************/ /* END LEGAL */ -#ifndef Hadrons_Baryon_hpp_ -#define Hadrons_Baryon_hpp_ +#ifndef Hadrons_MContraction_Baryon_hpp_ +#define Hadrons_MContraction_Baryon_hpp_ #include #include @@ -55,9 +55,9 @@ template class TBaryon: public Module { public: - TYPE_ALIASES(FImpl1, 1); - TYPE_ALIASES(FImpl2, 2); - TYPE_ALIASES(FImpl3, 3); + FERM_TYPE_ALIASES(FImpl1, 1); + FERM_TYPE_ALIASES(FImpl2, 2); + FERM_TYPE_ALIASES(FImpl3, 3); class Result: Serializable { public: @@ -121,11 +121,11 @@ void TBaryon::execute(void) // FIXME: do contractions - write(writer, "meson", result); + // write(writer, "meson", result); } END_MODULE_NAMESPACE END_HADRONS_NAMESPACE -#endif // Hadrons_Baryon_hpp_ +#endif // Hadrons_MContraction_Baryon_hpp_ diff --git a/extras/Hadrons/Modules/MContraction/DiscLoop.hpp b/extras/Hadrons/Modules/MContraction/DiscLoop.hpp index 4ad12e90..4f782cd3 100644 --- a/extras/Hadrons/Modules/MContraction/DiscLoop.hpp +++ b/extras/Hadrons/Modules/MContraction/DiscLoop.hpp @@ -26,8 +26,8 @@ See the full license in the file "LICENSE" in the top level distribution directo *************************************************************************************/ /* END LEGAL */ -#ifndef Hadrons_DiscLoop_hpp_ -#define Hadrons_DiscLoop_hpp_ +#ifndef Hadrons_MContraction_DiscLoop_hpp_ +#define Hadrons_MContraction_DiscLoop_hpp_ #include #include @@ -52,7 +52,7 @@ public: template class TDiscLoop: public Module { - TYPE_ALIASES(FImpl,); + FERM_TYPE_ALIASES(FImpl,); class Result: Serializable { public: @@ -141,4 +141,4 @@ END_MODULE_NAMESPACE END_HADRONS_NAMESPACE -#endif // Hadrons_DiscLoop_hpp_ +#endif // Hadrons_MContraction_DiscLoop_hpp_ diff --git a/extras/Hadrons/Modules/MContraction/Gamma3pt.hpp b/extras/Hadrons/Modules/MContraction/Gamma3pt.hpp index e5e73fa6..7f643d49 100644 --- a/extras/Hadrons/Modules/MContraction/Gamma3pt.hpp +++ b/extras/Hadrons/Modules/MContraction/Gamma3pt.hpp @@ -26,8 +26,8 @@ See the full license in the file "LICENSE" in the top level distribution directo *************************************************************************************/ /* END LEGAL */ -#ifndef Hadrons_Gamma3pt_hpp_ -#define Hadrons_Gamma3pt_hpp_ +#ifndef Hadrons_MContraction_Gamma3pt_hpp_ +#define Hadrons_MContraction_Gamma3pt_hpp_ #include #include @@ -72,9 +72,9 @@ public: template class TGamma3pt: public Module { - TYPE_ALIASES(FImpl1, 1); - TYPE_ALIASES(FImpl2, 2); - TYPE_ALIASES(FImpl3, 3); + FERM_TYPE_ALIASES(FImpl1, 1); + FERM_TYPE_ALIASES(FImpl2, 2); + FERM_TYPE_ALIASES(FImpl3, 3); class Result: Serializable { public: @@ -167,4 +167,4 @@ END_MODULE_NAMESPACE END_HADRONS_NAMESPACE -#endif // Hadrons_Gamma3pt_hpp_ +#endif // Hadrons_MContraction_Gamma3pt_hpp_ diff --git a/extras/Hadrons/Modules/MContraction/Meson.hpp b/extras/Hadrons/Modules/MContraction/Meson.hpp index 09c2a6e1..7810326a 100644 --- a/extras/Hadrons/Modules/MContraction/Meson.hpp +++ b/extras/Hadrons/Modules/MContraction/Meson.hpp @@ -29,8 +29,8 @@ See the full license in the file "LICENSE" in the top level distribution directo *************************************************************************************/ /* END LEGAL */ -#ifndef Hadrons_Meson_hpp_ -#define Hadrons_Meson_hpp_ +#ifndef Hadrons_MContraction_Meson_hpp_ +#define Hadrons_MContraction_Meson_hpp_ #include #include @@ -69,7 +69,7 @@ public: std::string, q1, std::string, q2, std::string, gammas, - std::string, mom, + std::string, sink, std::string, output); }; @@ -77,8 +77,10 @@ template class TMeson: public Module { public: - TYPE_ALIASES(FImpl1, 1); - TYPE_ALIASES(FImpl2, 2); + FERM_TYPE_ALIASES(FImpl1, 1); + FERM_TYPE_ALIASES(FImpl2, 2); + FERM_TYPE_ALIASES(ScalarImplCR, Scalar); + SINK_TYPE_ALIASES(Scalar); class Result: Serializable { public: @@ -115,7 +117,7 @@ TMeson::TMeson(const std::string name) template std::vector TMeson::getInput(void) { - std::vector input = {par().q1, par().q2}; + std::vector input = {par().q1, par().q2, par().sink}; return input; } @@ -154,6 +156,9 @@ void TMeson::parseGammaString(std::vector &gammaList) // execution /////////////////////////////////////////////////////////////////// +#define mesonConnected(q1, q2, gSnk, gSrc) \ +(g5*(gSnk))*(q1)*(adj(gSrc)*g5)*adj(q2) + template void TMeson::execute(void) { @@ -161,43 +166,72 @@ void TMeson::execute(void) << " quarks '" << par().q1 << "' and '" << par().q2 << "'" << std::endl; - CorrWriter writer(par().output); - PropagatorField1 &q1 = *env().template getObject(par().q1); - PropagatorField2 &q2 = *env().template getObject(par().q2); - LatticeComplex c(env().getGrid()); - Gamma g5(Gamma::Algebra::Gamma5); - std::vector gammaList; + CorrWriter writer(par().output); std::vector buf; std::vector result; - std::vector p; - - p = strToVec(par().mom); - LatticeComplex ph(env().getGrid()), coor(env().getGrid()); - Complex i(0.0,1.0); - 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 = exp((Real)(2*M_PI)*i*ph); + Gamma g5(Gamma::Algebra::Gamma5); + std::vector gammaList; + int nt = env().getDim(Tp); parseGammaString(gammaList); - result.resize(gammaList.size()); for (unsigned int i = 0; i < result.size(); ++i) { - Gamma gSnk(gammaList[i].first); - Gamma gSrc(gammaList[i].second); - c = trace((g5*gSnk)*q1*(adj(gSrc)*g5)*adj(q2))*ph; - sliceSum(c, buf, Tp); - result[i].gamma_snk = gammaList[i].first; result[i].gamma_src = gammaList[i].second; - result[i].corr.resize(buf.size()); - for (unsigned int t = 0; t < buf.size(); ++t) + result[i].corr.resize(nt); + } + if (env().template isObjectOfType(par().q1) and + env().template isObjectOfType(par().q2)) + { + SlicedPropagator1 &q1 = *env().template getObject(par().q1); + SlicedPropagator2 &q2 = *env().template getObject(par().q2); + + LOG(Message) << "(propagator already sinked)" << std::endl; + for (unsigned int i = 0; i < result.size(); ++i) { - result[i].corr[t] = TensorRemove(buf[t]); + Gamma gSnk(gammaList[i].first); + Gamma gSrc(gammaList[i].second); + + for (unsigned int t = 0; t < buf.size(); ++t) + { + result[i].corr[t] = TensorRemove(trace(mesonConnected(q1[t], q2[t], gSnk, gSrc))); + } + } + } + else + { + PropagatorField1 &q1 = *env().template getObject(par().q1); + PropagatorField2 &q2 = *env().template getObject(par().q2); + LatticeComplex c(env().getGrid()); + + LOG(Message) << "(using sink '" << par().sink << "')" << std::endl; + for (unsigned int i = 0; i < result.size(); ++i) + { + Gamma gSnk(gammaList[i].first); + Gamma gSrc(gammaList[i].second); + std::string ns; + + ns = env().getModuleNamespace(env().getObjectModule(par().sink)); + if (ns == "MSource") + { + PropagatorField1 &sink = + *env().template getObject(par().sink); + + c = trace(mesonConnected(q1, q2, gSnk, gSrc)*sink); + sliceSum(c, buf, Tp); + } + else if (ns == "MSink") + { + SinkFnScalar &sink = *env().template getObject(par().sink); + + c = trace(mesonConnected(q1, q2, gSnk, gSrc)); + buf = sink(c); + } + for (unsigned int t = 0; t < buf.size(); ++t) + { + result[i].corr[t] = TensorRemove(buf[t]); + } } } write(writer, "meson", result); @@ -207,4 +241,4 @@ END_MODULE_NAMESPACE END_HADRONS_NAMESPACE -#endif // Hadrons_Meson_hpp_ +#endif // Hadrons_MContraction_Meson_hpp_ diff --git a/extras/Hadrons/Modules/MContraction/WeakHamiltonian.hpp b/extras/Hadrons/Modules/MContraction/WeakHamiltonian.hpp index 23482feb..0a3c2e31 100644 --- a/extras/Hadrons/Modules/MContraction/WeakHamiltonian.hpp +++ b/extras/Hadrons/Modules/MContraction/WeakHamiltonian.hpp @@ -26,8 +26,8 @@ See the full license in the file "LICENSE" in the top level distribution directo *************************************************************************************/ /* END LEGAL */ -#ifndef Hadrons_WeakHamiltonian_hpp_ -#define Hadrons_WeakHamiltonian_hpp_ +#ifndef Hadrons_MContraction_WeakHamiltonian_hpp_ +#define Hadrons_MContraction_WeakHamiltonian_hpp_ #include #include @@ -83,7 +83,7 @@ public: class T##modname: public Module\ {\ public:\ - TYPE_ALIASES(FIMPL,)\ + FERM_TYPE_ALIASES(FIMPL,)\ class Result: Serializable\ {\ public:\ @@ -111,4 +111,4 @@ END_MODULE_NAMESPACE END_HADRONS_NAMESPACE -#endif // Hadrons_WeakHamiltonian_hpp_ +#endif // Hadrons_MContraction_WeakHamiltonian_hpp_ diff --git a/extras/Hadrons/Modules/MContraction/WeakHamiltonianEye.hpp b/extras/Hadrons/Modules/MContraction/WeakHamiltonianEye.hpp index 2ee87895..3a2b9309 100644 --- a/extras/Hadrons/Modules/MContraction/WeakHamiltonianEye.hpp +++ b/extras/Hadrons/Modules/MContraction/WeakHamiltonianEye.hpp @@ -26,8 +26,8 @@ See the full license in the file "LICENSE" in the top level distribution directo *************************************************************************************/ /* END LEGAL */ -#ifndef Hadrons_WeakHamiltonianEye_hpp_ -#define Hadrons_WeakHamiltonianEye_hpp_ +#ifndef Hadrons_MContraction_WeakHamiltonianEye_hpp_ +#define Hadrons_MContraction_WeakHamiltonianEye_hpp_ #include @@ -55,4 +55,4 @@ END_MODULE_NAMESPACE END_HADRONS_NAMESPACE -#endif // Hadrons_WeakHamiltonianEye_hpp_ +#endif // Hadrons_MContraction_WeakHamiltonianEye_hpp_ diff --git a/extras/Hadrons/Modules/MContraction/WeakHamiltonianNonEye.hpp b/extras/Hadrons/Modules/MContraction/WeakHamiltonianNonEye.hpp index 69bb8005..eb5abe3c 100644 --- a/extras/Hadrons/Modules/MContraction/WeakHamiltonianNonEye.hpp +++ b/extras/Hadrons/Modules/MContraction/WeakHamiltonianNonEye.hpp @@ -26,8 +26,8 @@ See the full license in the file "LICENSE" in the top level distribution directo *************************************************************************************/ /* END LEGAL */ -#ifndef Hadrons_WeakHamiltonianNonEye_hpp_ -#define Hadrons_WeakHamiltonianNonEye_hpp_ +#ifndef Hadrons_MContraction_WeakHamiltonianNonEye_hpp_ +#define Hadrons_MContraction_WeakHamiltonianNonEye_hpp_ #include @@ -54,4 +54,4 @@ END_MODULE_NAMESPACE END_HADRONS_NAMESPACE -#endif // Hadrons_WeakHamiltonianNonEye_hpp_ +#endif // Hadrons_MContraction_WeakHamiltonianNonEye_hpp_ diff --git a/extras/Hadrons/Modules/MContraction/WeakNeutral4ptDisc.hpp b/extras/Hadrons/Modules/MContraction/WeakNeutral4ptDisc.hpp index c0d8f829..f26d4636 100644 --- a/extras/Hadrons/Modules/MContraction/WeakNeutral4ptDisc.hpp +++ b/extras/Hadrons/Modules/MContraction/WeakNeutral4ptDisc.hpp @@ -26,8 +26,8 @@ See the full license in the file "LICENSE" in the top level distribution directo *************************************************************************************/ /* END LEGAL */ -#ifndef Hadrons_WeakNeutral4ptDisc_hpp_ -#define Hadrons_WeakNeutral4ptDisc_hpp_ +#ifndef Hadrons_MContraction_WeakNeutral4ptDisc_hpp_ +#define Hadrons_MContraction_WeakNeutral4ptDisc_hpp_ #include @@ -56,4 +56,4 @@ END_MODULE_NAMESPACE END_HADRONS_NAMESPACE -#endif // Hadrons_WeakNeutral4ptDisc_hpp_ +#endif // Hadrons_MContraction_WeakNeutral4ptDisc_hpp_ diff --git a/extras/Hadrons/Modules/MGauge/Load.hpp b/extras/Hadrons/Modules/MGauge/Load.hpp index c41f9b8c..5ff6da0f 100644 --- a/extras/Hadrons/Modules/MGauge/Load.hpp +++ b/extras/Hadrons/Modules/MGauge/Load.hpp @@ -27,8 +27,8 @@ See the full license in the file "LICENSE" in the top level distribution directo *************************************************************************************/ /* END LEGAL */ -#ifndef Hadrons_Load_hpp_ -#define Hadrons_Load_hpp_ +#ifndef Hadrons_MGauge_Load_hpp_ +#define Hadrons_MGauge_Load_hpp_ #include #include @@ -70,4 +70,4 @@ END_MODULE_NAMESPACE END_HADRONS_NAMESPACE -#endif // Hadrons_Load_hpp_ +#endif // Hadrons_MGauge_Load_hpp_ diff --git a/extras/Hadrons/Modules/MGauge/Random.hpp b/extras/Hadrons/Modules/MGauge/Random.hpp index e3fbcf1a..a97d25cf 100644 --- a/extras/Hadrons/Modules/MGauge/Random.hpp +++ b/extras/Hadrons/Modules/MGauge/Random.hpp @@ -27,8 +27,8 @@ See the full license in the file "LICENSE" in the top level distribution directo *************************************************************************************/ /* END LEGAL */ -#ifndef Hadrons_Random_hpp_ -#define Hadrons_Random_hpp_ +#ifndef Hadrons_MGauge_Random_hpp_ +#define Hadrons_MGauge_Random_hpp_ #include #include @@ -63,4 +63,4 @@ END_MODULE_NAMESPACE END_HADRONS_NAMESPACE -#endif // Hadrons_Random_hpp_ +#endif // Hadrons_MGauge_Random_hpp_ diff --git a/extras/Hadrons/Modules/MGauge/StochEm.hpp b/extras/Hadrons/Modules/MGauge/StochEm.hpp index 50a77435..12ce9fdc 100644 --- a/extras/Hadrons/Modules/MGauge/StochEm.hpp +++ b/extras/Hadrons/Modules/MGauge/StochEm.hpp @@ -25,8 +25,8 @@ with this program; if not, write to the Free Software Foundation, Inc., See the full license in the file "LICENSE" in the top level distribution directory *************************************************************************************/ /* END LEGAL */ -#ifndef Hadrons_StochEm_hpp_ -#define Hadrons_StochEm_hpp_ +#ifndef Hadrons_MGauge_StochEm_hpp_ +#define Hadrons_MGauge_StochEm_hpp_ #include #include @@ -72,4 +72,4 @@ END_MODULE_NAMESPACE END_HADRONS_NAMESPACE -#endif // Hadrons_StochEm_hpp_ +#endif // Hadrons_MGauge_StochEm_hpp_ diff --git a/extras/Hadrons/Modules/MGauge/Unit.hpp b/extras/Hadrons/Modules/MGauge/Unit.hpp index 2ff10bfd..7cd15ef7 100644 --- a/extras/Hadrons/Modules/MGauge/Unit.hpp +++ b/extras/Hadrons/Modules/MGauge/Unit.hpp @@ -27,8 +27,8 @@ See the full license in the file "LICENSE" in the top level distribution directo *************************************************************************************/ /* END LEGAL */ -#ifndef Hadrons_Unit_hpp_ -#define Hadrons_Unit_hpp_ +#ifndef Hadrons_MGauge_Unit_hpp_ +#define Hadrons_MGauge_Unit_hpp_ #include #include @@ -63,4 +63,4 @@ END_MODULE_NAMESPACE END_HADRONS_NAMESPACE -#endif // Hadrons_Unit_hpp_ +#endif // Hadrons_MGauge_Unit_hpp_ diff --git a/extras/Hadrons/Modules/MLoop/NoiseLoop.hpp b/extras/Hadrons/Modules/MLoop/NoiseLoop.hpp index 3d2850d1..5d2c4a13 100644 --- a/extras/Hadrons/Modules/MLoop/NoiseLoop.hpp +++ b/extras/Hadrons/Modules/MLoop/NoiseLoop.hpp @@ -26,8 +26,8 @@ See the full license in the file "LICENSE" in the top level distribution directo *************************************************************************************/ /* END LEGAL */ -#ifndef Hadrons_NoiseLoop_hpp_ -#define Hadrons_NoiseLoop_hpp_ +#ifndef Hadrons_MLoop_NoiseLoop_hpp_ +#define Hadrons_MLoop_NoiseLoop_hpp_ #include #include @@ -65,7 +65,7 @@ template class TNoiseLoop: public Module { public: - TYPE_ALIASES(FImpl,); + FERM_TYPE_ALIASES(FImpl,); public: // constructor TNoiseLoop(const std::string name); @@ -129,4 +129,4 @@ END_MODULE_NAMESPACE END_HADRONS_NAMESPACE -#endif // Hadrons_NoiseLoop_hpp_ +#endif // Hadrons_MLoop_NoiseLoop_hpp_ diff --git a/extras/Hadrons/Modules/MScalar/ChargedProp.hpp b/extras/Hadrons/Modules/MScalar/ChargedProp.hpp index 8bb5faa0..fbe75c05 100644 --- a/extras/Hadrons/Modules/MScalar/ChargedProp.hpp +++ b/extras/Hadrons/Modules/MScalar/ChargedProp.hpp @@ -1,5 +1,5 @@ -#ifndef Hadrons_ChargedProp_hpp_ -#define Hadrons_ChargedProp_hpp_ +#ifndef Hadrons_MScalar_ChargedProp_hpp_ +#define Hadrons_MScalar_ChargedProp_hpp_ #include #include @@ -58,4 +58,4 @@ END_MODULE_NAMESPACE END_HADRONS_NAMESPACE -#endif // Hadrons_ChargedProp_hpp_ +#endif // Hadrons_MScalar_ChargedProp_hpp_ diff --git a/extras/Hadrons/Modules/MScalar/FreeProp.hpp b/extras/Hadrons/Modules/MScalar/FreeProp.hpp index 29f15eda..97cf288a 100644 --- a/extras/Hadrons/Modules/MScalar/FreeProp.hpp +++ b/extras/Hadrons/Modules/MScalar/FreeProp.hpp @@ -1,5 +1,5 @@ -#ifndef Hadrons_FreeProp_hpp_ -#define Hadrons_FreeProp_hpp_ +#ifndef Hadrons_MScalar_FreeProp_hpp_ +#define Hadrons_MScalar_FreeProp_hpp_ #include #include @@ -47,4 +47,4 @@ END_MODULE_NAMESPACE END_HADRONS_NAMESPACE -#endif // Hadrons_FreeProp_hpp_ +#endif // Hadrons_MScalar_FreeProp_hpp_ diff --git a/extras/Hadrons/Modules/MSink/Point.hpp b/extras/Hadrons/Modules/MSink/Point.hpp new file mode 100644 index 00000000..7b3aa9de --- /dev/null +++ b/extras/Hadrons/Modules/MSink/Point.hpp @@ -0,0 +1,114 @@ +#ifndef Hadrons_MSink_Point_hpp_ +#define Hadrons_MSink_Point_hpp_ + +#include +#include +#include + +BEGIN_HADRONS_NAMESPACE + +/****************************************************************************** + * Point * + ******************************************************************************/ +BEGIN_MODULE_NAMESPACE(MSink) + +class PointPar: Serializable +{ +public: + GRID_SERIALIZABLE_CLASS_MEMBERS(PointPar, + std::string, mom); +}; + +template +class TPoint: public Module +{ +public: + FERM_TYPE_ALIASES(FImpl,); + SINK_TYPE_ALIASES(); +public: + // constructor + TPoint(const std::string name); + // destructor + virtual ~TPoint(void) = default; + // dependency relation + virtual std::vector getInput(void); + virtual std::vector getOutput(void); + // setup + virtual void setup(void); + // execution + virtual void execute(void); +}; + +MODULE_REGISTER_NS(Point, TPoint, MSink); +MODULE_REGISTER_NS(ScalarPoint, TPoint, MSink); + +/****************************************************************************** + * TPoint implementation * + ******************************************************************************/ +// constructor ///////////////////////////////////////////////////////////////// +template +TPoint::TPoint(const std::string name) +: Module(name) +{} + +// dependencies/products /////////////////////////////////////////////////////// +template +std::vector TPoint::getInput(void) +{ + std::vector in; + + return in; +} + +template +std::vector TPoint::getOutput(void) +{ + std::vector out = {getName()}; + + return out; +} + +// setup /////////////////////////////////////////////////////////////////////// +template +void TPoint::setup(void) +{ + unsigned int size; + + size = env().template lattice4dSize(); + env().registerObject(getName(), size); +} + +// execution /////////////////////////////////////////////////////////////////// +template +void TPoint::execute(void) +{ + std::vector p = strToVec(par().mom); + LatticeComplex ph(env().getGrid()), coor(env().getGrid()); + Complex i(0.0,1.0); + + LOG(Message) << "Setting up point sink function for momentum [" + << par().mom << "]" << std::endl; + ph = zero; + for(unsigned int mu = 0; mu < env().getNd(); mu++) + { + LatticeCoordinate(coor, mu); + ph = ph + (p[mu]/env().getGrid()->_fdimensions[mu])*coor; + } + ph = exp((Real)(2*M_PI)*i*ph); + auto sink = [ph](const PropagatorField &field) + { + SlicedPropagator res; + PropagatorField tmp = ph*field; + + sliceSum(tmp, res, Tp); + + return res; + }; + env().setObject(getName(), new SinkFn(sink)); +} + +END_MODULE_NAMESPACE + +END_HADRONS_NAMESPACE + +#endif // Hadrons_MSink_Point_hpp_ diff --git a/extras/Hadrons/Modules/MSolver/RBPrecCG.hpp b/extras/Hadrons/Modules/MSolver/RBPrecCG.hpp index d7220271..b1f63a5d 100644 --- a/extras/Hadrons/Modules/MSolver/RBPrecCG.hpp +++ b/extras/Hadrons/Modules/MSolver/RBPrecCG.hpp @@ -27,8 +27,8 @@ See the full license in the file "LICENSE" in the top level distribution directo *************************************************************************************/ /* END LEGAL */ -#ifndef Hadrons_RBPrecCG_hpp_ -#define Hadrons_RBPrecCG_hpp_ +#ifndef Hadrons_MSolver_RBPrecCG_hpp_ +#define Hadrons_MSolver_RBPrecCG_hpp_ #include #include @@ -53,7 +53,7 @@ template class TRBPrecCG: public Module { public: - TYPE_ALIASES(FImpl,); + FGS_TYPE_ALIASES(FImpl,); public: // constructor TRBPrecCG(const std::string name); @@ -129,4 +129,4 @@ END_MODULE_NAMESPACE END_HADRONS_NAMESPACE -#endif // Hadrons_RBPrecCG_hpp_ +#endif // Hadrons_MSolver_RBPrecCG_hpp_ diff --git a/extras/Hadrons/Modules/MSource/Point.hpp b/extras/Hadrons/Modules/MSource/Point.hpp index 3c0fc9a1..0c415807 100644 --- a/extras/Hadrons/Modules/MSource/Point.hpp +++ b/extras/Hadrons/Modules/MSource/Point.hpp @@ -27,8 +27,8 @@ See the full license in the file "LICENSE" in the top level distribution directo *************************************************************************************/ /* END LEGAL */ -#ifndef Hadrons_Point_hpp_ -#define Hadrons_Point_hpp_ +#ifndef Hadrons_MSource_Point_hpp_ +#define Hadrons_MSource_Point_hpp_ #include #include @@ -133,4 +133,4 @@ END_MODULE_NAMESPACE END_HADRONS_NAMESPACE -#endif // Hadrons_Point_hpp_ +#endif // Hadrons_MSource_Point_hpp_ diff --git a/extras/Hadrons/Modules/MSource/SeqGamma.hpp b/extras/Hadrons/Modules/MSource/SeqGamma.hpp index 366ebee7..e2129a46 100644 --- a/extras/Hadrons/Modules/MSource/SeqGamma.hpp +++ b/extras/Hadrons/Modules/MSource/SeqGamma.hpp @@ -28,8 +28,8 @@ See the full license in the file "LICENSE" in the top level distribution directo *************************************************************************************/ /* END LEGAL */ -#ifndef Hadrons_SeqGamma_hpp_ -#define Hadrons_SeqGamma_hpp_ +#ifndef Hadrons_MSource_SeqGamma_hpp_ +#define Hadrons_MSource_SeqGamma_hpp_ #include #include @@ -72,7 +72,7 @@ template class TSeqGamma: public Module { public: - TYPE_ALIASES(FImpl,); + FGS_TYPE_ALIASES(FImpl,); public: // constructor TSeqGamma(const std::string name); @@ -161,4 +161,4 @@ END_MODULE_NAMESPACE END_HADRONS_NAMESPACE -#endif // Hadrons_SeqGamma_hpp_ +#endif // Hadrons_MSource_SeqGamma_hpp_ diff --git a/extras/Hadrons/Modules/MSource/Wall.hpp b/extras/Hadrons/Modules/MSource/Wall.hpp index 8722876f..4de37e4d 100644 --- a/extras/Hadrons/Modules/MSource/Wall.hpp +++ b/extras/Hadrons/Modules/MSource/Wall.hpp @@ -26,8 +26,8 @@ See the full license in the file "LICENSE" in the top level distribution directo *************************************************************************************/ /* END LEGAL */ -#ifndef Hadrons_WallSource_hpp_ -#define Hadrons_WallSource_hpp_ +#ifndef Hadrons_MSource_WallSource_hpp_ +#define Hadrons_MSource_WallSource_hpp_ #include #include @@ -64,7 +64,7 @@ template class TWall: public Module { public: - TYPE_ALIASES(FImpl,); + FERM_TYPE_ALIASES(FImpl,); public: // constructor TWall(const std::string name); @@ -144,4 +144,4 @@ END_MODULE_NAMESPACE END_HADRONS_NAMESPACE -#endif // Hadrons_WallSource_hpp_ +#endif // Hadrons_MSource_WallSource_hpp_ diff --git a/extras/Hadrons/Modules/MSource/Z2.hpp b/extras/Hadrons/Modules/MSource/Z2.hpp index 761ae139..a7f7a3e6 100644 --- a/extras/Hadrons/Modules/MSource/Z2.hpp +++ b/extras/Hadrons/Modules/MSource/Z2.hpp @@ -27,8 +27,8 @@ See the full license in the file "LICENSE" in the top level distribution directo *************************************************************************************/ /* END LEGAL */ -#ifndef Hadrons_Z2_hpp_ -#define Hadrons_Z2_hpp_ +#ifndef Hadrons_MSource_Z2_hpp_ +#define Hadrons_MSource_Z2_hpp_ #include #include @@ -149,4 +149,4 @@ END_MODULE_NAMESPACE END_HADRONS_NAMESPACE -#endif // Hadrons_Z2_hpp_ +#endif // Hadrons_MSource_Z2_hpp_ diff --git a/extras/Hadrons/Modules/Quark.hpp b/extras/Hadrons/Modules/Quark.hpp index c0d1f65a..cf7d4c28 100644 --- a/extras/Hadrons/Modules/Quark.hpp +++ b/extras/Hadrons/Modules/Quark.hpp @@ -51,7 +51,7 @@ template class TQuark: public Module { public: - TYPE_ALIASES(FImpl,); + FGS_TYPE_ALIASES(FImpl,); public: // constructor TQuark(const std::string name); diff --git a/extras/Hadrons/Modules/templates/Module_in_NS.hpp.template b/extras/Hadrons/Modules/templates/Module_in_NS.hpp.template index ece2bb58..ea77b12a 100644 --- a/extras/Hadrons/Modules/templates/Module_in_NS.hpp.template +++ b/extras/Hadrons/Modules/templates/Module_in_NS.hpp.template @@ -1,5 +1,5 @@ -#ifndef Hadrons____FILEBASENAME____hpp_ -#define Hadrons____FILEBASENAME____hpp_ +#ifndef Hadrons____NAMESPACE_______FILEBASENAME____hpp_ +#define Hadrons____NAMESPACE_______FILEBASENAME____hpp_ #include #include @@ -41,4 +41,4 @@ END_MODULE_NAMESPACE END_HADRONS_NAMESPACE -#endif // Hadrons____FILEBASENAME____hpp_ +#endif // Hadrons____NAMESPACE_______FILEBASENAME____hpp_ diff --git a/extras/Hadrons/Modules/templates/Module_tmp_in_NS.hpp.template b/extras/Hadrons/Modules/templates/Module_tmp_in_NS.hpp.template index a330652d..b79c0ad3 100644 --- a/extras/Hadrons/Modules/templates/Module_tmp_in_NS.hpp.template +++ b/extras/Hadrons/Modules/templates/Module_tmp_in_NS.hpp.template @@ -1,5 +1,5 @@ -#ifndef Hadrons____FILEBASENAME____hpp_ -#define Hadrons____FILEBASENAME____hpp_ +#ifndef Hadrons____NAMESPACE_______FILEBASENAME____hpp_ +#define Hadrons____NAMESPACE_______FILEBASENAME____hpp_ #include #include @@ -82,4 +82,4 @@ END_MODULE_NAMESPACE END_HADRONS_NAMESPACE -#endif // Hadrons____FILEBASENAME____hpp_ +#endif // Hadrons____NAMESPACE_______FILEBASENAME____hpp_ diff --git a/extras/Hadrons/modules.inc b/extras/Hadrons/modules.inc index 3cf69144..f51ede5a 100644 --- a/extras/Hadrons/modules.inc +++ b/extras/Hadrons/modules.inc @@ -28,6 +28,7 @@ modules_hpp =\ Modules/MScalar/ChargedProp.hpp \ Modules/MScalar/FreeProp.hpp \ Modules/MScalar/Scalar.hpp \ + Modules/MSink/Point.hpp \ Modules/MSolver/RBPrecCG.hpp \ Modules/MSource/Point.hpp \ Modules/MSource/SeqGamma.hpp \ diff --git a/tests/hadrons/Test_hadrons_spectrum.cc b/tests/hadrons/Test_hadrons_spectrum.cc index 55f3346e..8f7b30c8 100644 --- a/tests/hadrons/Test_hadrons_spectrum.cc +++ b/tests/hadrons/Test_hadrons_spectrum.cc @@ -63,6 +63,10 @@ int main(int argc, char *argv[]) MSource::Point::Par ptPar; ptPar.position = "0 0 0 0"; application.createModule("pt", ptPar); + // sink + MSink::Point::Par sinkPar; + sinkPar.mom = "0 0 0"; + application.createModule("sink", sinkPar); // set fermion boundary conditions to be periodic space, antiperiodic time. std::string boundary = "1 1 1 -1"; @@ -98,19 +102,19 @@ int main(int argc, char *argv[]) { MContraction::Meson::Par mesPar; - mesPar.output = "mesons/pt_" + flavour[i] + flavour[j]; - mesPar.q1 = "Qpt_" + flavour[i]; - mesPar.q2 = "Qpt_" + flavour[j]; - mesPar.gammas = "all"; - mesPar.mom = "0. 0. 0. 0."; + mesPar.output = "mesons/pt_" + flavour[i] + flavour[j]; + mesPar.q1 = "Qpt_" + flavour[i]; + mesPar.q2 = "Qpt_" + flavour[j]; + mesPar.gammas = "all"; + mesPar.sink = "sink"; application.createModule("meson_pt_" + flavour[i] + flavour[j], mesPar); - mesPar.output = "mesons/Z2_" + flavour[i] + flavour[j]; - mesPar.q1 = "QZ2_" + flavour[i]; - mesPar.q2 = "QZ2_" + flavour[j]; - mesPar.gammas = "all"; - mesPar.mom = "0. 0. 0. 0."; + mesPar.output = "mesons/Z2_" + flavour[i] + flavour[j]; + mesPar.q1 = "QZ2_" + flavour[i]; + mesPar.q2 = "QZ2_" + flavour[j]; + mesPar.gammas = "all"; + mesPar.sink = "sink"; application.createModule("meson_Z2_" + flavour[i] + flavour[j], mesPar);