diff --git a/benchmarks/Benchmark_comms.cc b/benchmarks/Benchmark_comms.cc index 99ab190b..ce881ef6 100644 --- a/benchmarks/Benchmark_comms.cc +++ b/benchmarks/Benchmark_comms.cc @@ -31,6 +31,32 @@ using namespace std; using namespace Grid; using namespace Grid::QCD; +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(), [=](double x) { return x - mean; }); + double sq_sum = std::inner_product(diff.begin(), diff.end(), diff.begin(), 0.0); + 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 <1) nmu++; + std::cout << GridLogMessage << "Number of iterations to average: "<< Nloop << std::endl; + std::vector t_time(Nloop); + time_statistics timestat; + std::cout< requests; @@ -102,18 +132,24 @@ int main (int argc, char ** argv) } Grid.SendToRecvFromComplete(requests); Grid.Barrier(); - + double stop=usecond(); + t_time[i] = stop-start; // microseconds } - double stop=usecond(); + + timestat.statistics(t_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; @@ -258,28 +300,34 @@ int main (int argc, char ** argv) } Grid.StencilSendToRecvFromComplete(requests); Grid.Barrier(); + double stop=usecond(); + t_time[i] = stop-start; // microseconds } - double stop=usecond(); + + timestat.statistics(t_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; @@ -341,19 +389,27 @@ int main (int argc, char ** argv) } } - Grid.Barrier(); + Grid.Barrier(); + double stop=usecond(); + t_time[i] = stop-start; // microseconds } - double stop=usecond(); + + timestat.statistics(t_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<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 ebf93283..371256e8 100644 --- a/extras/Hadrons/Global.hpp +++ b/extras/Hadrons/Global.hpp @@ -51,22 +51,41 @@ 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 ScalarImplCR +#endif BEGIN_HADRONS_NAMESPACE // type aliases -#define TYPE_ALIASES(FImpl, suffix)\ -typedef FermionOperator FMat##suffix; \ -typedef typename FImpl::FermionField FermionField##suffix; \ -typedef typename FImpl::PropagatorField PropagatorField##suffix; \ -typedef typename FImpl::SitePropagator::scalar_object SitePropagator##suffix; \ -typedef typename FImpl::DoubledGaugeField DoubledGaugeField##suffix;\ -typedef std::function SolverFn##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::scalar_object SitePropagator##suffix; \ +typedef std::vector SlicedPropagator##suffix; + +#define GAUGE_TYPE_ALIASES(FImpl, suffix)\ +typedef typename FImpl::DoubledGaugeField DoubledGaugeField##suffix; + +#define SCALAR_TYPE_ALIASES(SImpl, suffix)\ +typedef typename SImpl::Field ScalarField##suffix;\ +typedef typename SImpl::Field PropagatorField##suffix; + +#define SOLVER_TYPE_ALIASES(FImpl, suffix)\ +typedef std::function SolverFn##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) // logger class HadronsLogger: public Logger diff --git a/extras/Hadrons/Modules.hpp b/extras/Hadrons/Modules.hpp index 53ec346c..5d2be8ee 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 @@ -39,8 +11,13 @@ See the full license in the file "LICENSE" in the top level distribution directo #include #include #include +#include #include #include +#include +#include +#include +#include #include #include #include diff --git a/extras/Hadrons/Modules/MAction/DWF.hpp b/extras/Hadrons/Modules/MAction/DWF.hpp index 49861e3e..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 @@ -48,14 +48,15 @@ public: std::string, gauge, unsigned int, Ls, double , mass, - double , M5); + double , M5, + std::string , boundary); }; template class TDWF: public Module { public: - TYPE_ALIASES(FImpl,); + FGS_TYPE_ALIASES(FImpl,); public: // constructor TDWF(const std::string name); @@ -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); } @@ -131,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 6ffa997d..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 @@ -46,14 +46,15 @@ class WilsonPar: Serializable public: GRID_SERIALIZABLE_CLASS_MEMBERS(WilsonPar, std::string, gauge, - double , mass); + double , mass, + std::string, boundary); }; template class TWilson: public Module { public: - TYPE_ALIASES(FImpl,); + FGS_TYPE_ALIASES(FImpl,); public: // constructor TWilson(const std::string name); @@ -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/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 4cbe1ac4..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; } @@ -131,12 +133,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) @@ -155,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) { @@ -162,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); @@ -208,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.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..12ce9fdc --- /dev/null +++ b/extras/Hadrons/Modules/MGauge/StochEm.hpp @@ -0,0 +1,75 @@ +/************************************************************************************* + +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_MGauge_StochEm_hpp_ +#define Hadrons_MGauge_StochEm_hpp_ + +#include +#include +#include + +BEGIN_HADRONS_NAMESPACE + +/****************************************************************************** + * StochEm * + ******************************************************************************/ +BEGIN_MODULE_NAMESPACE(MGauge) + +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_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.cc b/extras/Hadrons/Modules/MScalar/ChargedProp.cc new file mode 100644 index 00000000..cd8dc244 --- /dev/null +++ b/extras/Hadrons/Modules/MScalar/ChargedProp.cc @@ -0,0 +1,226 @@ +#include +#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 = {par().source, par().emField}; + + return in; +} + +std::vector TChargedProp::getOutput(void) +{ + std::vector out = {getName()}; + + return out; +} + +// setup /////////////////////////////////////////////////////////////////////// +void TChargedProp::setup(void) +{ + freeMomPropName_ = FREEMOMPROP(par().mass); + phaseName_.clear(); + for (unsigned int mu = 0; mu < env().getNd(); ++mu) + { + phaseName_.push_back("_shiftphase_" + std::to_string(mu)); + } + GFSrcName_ = "_" + getName() + "_DinvSrc"; + if (!env().hasRegisteredObject(freeMomPropName_)) + { + env().registerLattice(freeMomPropName_); + } + if (!env().hasRegisteredObject(phaseName_[0])) + { + for (unsigned int mu = 0; mu < env().getNd(); ++mu) + { + env().registerLattice(phaseName_[mu]); + } + } + if (!env().hasRegisteredObject(GFSrcName_)) + { + env().registerLattice(GFSrcName_); + } + env().registerLattice(getName()); +} + +// execution /////////////////////////////////////////////////////////////////// +void TChargedProp::execute(void) +{ + // CACHING ANALYTIC EXPRESSIONS + ScalarField &source = *env().getObject(par().source); + 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_); + SIMPL::MomentumSpacePropagator(*freeMomProp_, par().mass); + } + else + { + 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 shift phases..." << std::endl; + for (unsigned int mu = 0; mu < env().getNd(); ++mu) + { + Real twoPiL = M_PI*2./l[mu]; + + 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) + { + phase_.push_back(env().getObject(phaseName_[mu])); + } + } + + // 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; + + // 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 (momD2 = F*D2*Finv) + buf = GFSrc; + momD2(buf, fft); + 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; + + CorrWriter 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, "charge", q); + write(writer, "prop", result); + } +} + +void TChargedProp::momD1(ScalarField &s, FFT &fft) +{ + EmField &A = *env().getObject(par().emField); + ScalarField buf(env().getGrid()), result(env().getGrid()), + Amu(env().getGrid()); + Complex ci(0.0,1.0); + + result = zero; + + 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); + 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; +} + +void TChargedProp::momD2(ScalarField &s, FFT &fft) +{ + EmField &A = *env().getObject(par().emField); + ScalarField buf(env().getGrid()), result(env().getGrid()), + Amu(env().getGrid()); + + result = zero; + + 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); + 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; +} diff --git a/extras/Hadrons/Modules/MScalar/ChargedProp.hpp b/extras/Hadrons/Modules/MScalar/ChargedProp.hpp new file mode 100644 index 00000000..fbe75c05 --- /dev/null +++ b/extras/Hadrons/Modules/MScalar/ChargedProp.hpp @@ -0,0 +1,61 @@ +#ifndef Hadrons_MScalar_ChargedProp_hpp_ +#define Hadrons_MScalar_ChargedProp_hpp_ + +#include +#include +#include + +BEGIN_HADRONS_NAMESPACE + +/****************************************************************************** + * Charged scalar propagator * + ******************************************************************************/ +BEGIN_MODULE_NAMESPACE(MScalar) + +class ChargedPropPar: Serializable +{ +public: + GRID_SERIALIZABLE_CLASS_MEMBERS(ChargedPropPar, + std::string, emField, + std::string, source, + double, mass, + double, charge, + std::string, output); +}; + +class TChargedProp: public Module +{ +public: + SCALAR_TYPE_ALIASES(SIMPL,); + typedef PhotonR::GaugeField EmField; + typedef PhotonR::GaugeLinkField EmComp; +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); +private: + 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); + +END_MODULE_NAMESPACE + +END_HADRONS_NAMESPACE + +#endif // Hadrons_MScalar_ChargedProp_hpp_ diff --git a/extras/Hadrons/Modules/MScalar/FreeProp.cc b/extras/Hadrons/Modules/MScalar/FreeProp.cc new file mode 100644 index 00000000..674867e3 --- /dev/null +++ b/extras/Hadrons/Modules/MScalar/FreeProp.cc @@ -0,0 +1,79 @@ +#include +#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) +{ + freeMomPropName_ = FREEMOMPROP(par().mass); + + if (!env().hasRegisteredObject(freeMomPropName_)) + { + env().registerLattice(freeMomPropName_); + } + env().registerLattice(getName()); +} + +// execution /////////////////////////////////////////////////////////////////// +void TFreeProp::execute(void) +{ + ScalarField &prop = *env().createLattice(getName()); + ScalarField &source = *env().getObject(par().source); + ScalarField *freeMomProp; + + if (!env().hasCreatedObject(freeMomPropName_)) + { + LOG(Message) << "Caching momentum space free scalar propagator" + << " (mass= " << par().mass << ")..." << std::endl; + freeMomProp = env().createLattice(freeMomPropName_); + SIMPL::MomentumSpacePropagator(*freeMomProp, par().mass); + } + else + { + freeMomProp = env().getObject(freeMomPropName_); + } + LOG(Message) << "Computing free scalar propagator..." << std::endl; + SIMPL::FreePropagator(source, prop, *freeMomProp); + + 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 new file mode 100644 index 00000000..97cf288a --- /dev/null +++ b/extras/Hadrons/Modules/MScalar/FreeProp.hpp @@ -0,0 +1,50 @@ +#ifndef Hadrons_MScalar_FreeProp_hpp_ +#define Hadrons_MScalar_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, + std::string, output); +}; + +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); +private: + std::string freeMomPropName_; +}; + +MODULE_REGISTER_NS(FreeProp, TFreeProp, MScalar); + +END_MODULE_NAMESPACE + +END_HADRONS_NAMESPACE + +#endif // Hadrons_MScalar_FreeProp_hpp_ 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/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 36e1cc5b..7815e5c1 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 @@ -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 * @@ -132,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 cd5727be..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 @@ -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 * @@ -148,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 fff34edf..e22c3b1f 100644 --- a/extras/Hadrons/Modules/Quark.hpp +++ b/extras/Hadrons/Modules/Quark.hpp @@ -72,7 +72,7 @@ template class TQuark: public Module { public: - TYPE_ALIASES(FImpl,); + FGS_TYPE_ALIASES(FImpl,); public: // constructor TQuark(const std::string name); @@ -154,7 +154,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)) { 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 b57aa577..9fae0fe8 100644 --- a/extras/Hadrons/modules.inc +++ b/extras/Hadrons/modules.inc @@ -4,7 +4,10 @@ modules_cc =\ Modules/MContraction/WeakNeutral4ptDisc.cc \ Modules/MGauge/Load.cc \ Modules/MGauge/Random.cc \ - Modules/MGauge/Unit.cc + Modules/MGauge/StochEm.cc \ + Modules/MGauge/Unit.cc \ + Modules/MScalar/ChargedProp.cc \ + Modules/MScalar/FreeProp.cc modules_hpp =\ Modules/MAction/DWF.hpp \ @@ -20,8 +23,13 @@ modules_hpp =\ Modules/MContraction/WeakNeutral4ptDisc.hpp \ Modules/MGauge/Load.hpp \ Modules/MGauge/Random.hpp \ + Modules/MGauge/StochEm.hpp \ Modules/MGauge/Unit.hpp \ Modules/MLoop/NoiseLoop.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/SeqConserved.hpp \ diff --git a/extras/qed-fvol/Global.cc b/extras/qed-fvol/Global.cc new file mode 100644 index 00000000..57ed97cc --- /dev/null +++ b/extras/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/extras/qed-fvol/Global.hpp b/extras/qed-fvol/Global.hpp new file mode 100644 index 00000000..7f07200d --- /dev/null +++ b/extras/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/extras/qed-fvol/Makefile.am b/extras/qed-fvol/Makefile.am new file mode 100644 index 00000000..0a9030c7 --- /dev/null +++ b/extras/qed-fvol/Makefile.am @@ -0,0 +1,9 @@ +AM_CXXFLAGS += -I$(top_srcdir)/extras + +bin_PROGRAMS = qed-fvol + +qed_fvol_SOURCES = \ + qed-fvol.cc \ + Global.cc + +qed_fvol_LDADD = -lGrid diff --git a/extras/qed-fvol/WilsonLoops.h b/extras/qed-fvol/WilsonLoops.h new file mode 100644 index 00000000..98db6b7a --- /dev/null +++ b/extras/qed-fvol/WilsonLoops.h @@ -0,0 +1,265 @@ +#ifndef QEDFVOL_WILSONLOOPS_H +#define QEDFVOL_WILSONLOOPS_H + +#include + +BEGIN_QEDFVOL_NAMESPACE + +template class NewWilsonLoops : 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 < U[0]._grid->_ndimension; 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 Real sumPlaquette(const GaugeLorentz &Umu) { + std::vector U(4, Umu._grid); + + for (int mu = 0; mu < Umu._grid->_ndimension; 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 Real avgPlaquette(const GaugeLorentz &Umu) { + int ndim = Umu._grid->_ndimension; + Real sumplaq = sumPlaquette(Umu); + Real vol = Umu._grid->gSites(); + Real faces = (1.0 * ndim * (ndim - 1)) / 2.0; + return sumplaq / vol / faces / Nc; // 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 < 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; + 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 Real sumWilsonLoop(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); + + siteWilsonLoop(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 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); + 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 + } +}; + +END_QEDFVOL_NAMESPACE + +#endif // QEDFVOL_WILSONLOOPS_H \ No newline at end of file diff --git a/extras/qed-fvol/qed-fvol.cc b/extras/qed-fvol/qed-fvol.cc new file mode 100644 index 00000000..3ecac2fc --- /dev/null +++ b/extras/qed-fvol/qed-fvol.cc @@ -0,0 +1,88 @@ +#include +#include + +using namespace Grid; +using namespace QCD; +using namespace QedFVol; + +typedef PeriodicGaugeImpl QedPeriodicGimplR; +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 + 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; + + // 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); + EmField expA(&grid); + + Complex imag_unit(0, 1); + + Real wlA; + std::vector logWlAvg(NWILSON, 0.0), logWlTime(NWILSON, 0.0), logWlSpace(NWILSON, 0.0); + + pRNG.SeedRandomDevice(); + + 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 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 + LOG(Message) << "Grid is finalizing now" << std::endl; + Grid_finalize(); + + return EXIT_SUCCESS; +} diff --git a/lib/lattice/Lattice_unary.h b/lib/lattice/Lattice_unary.h index f5b324ec..44b7b4f1 100644 --- a/lib/lattice/Lattice_unary.h +++ b/lib/lattice/Lattice_unary.h @@ -62,14 +62,20 @@ 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); parallel_for(int ss=0;ssoSites();ss++){ ret._odata[ss]=Exponentiate(rhs._odata[ss],alpha, Nexp); } + return ret; + + + + + } diff --git a/lib/qcd/action/fermion/Fermion.h b/lib/qcd/action/fermion/Fermion.h index 55e71d20..0f803f44 100644 --- a/lib/qcd/action/fermion/Fermion.h +++ b/lib/qcd/action/fermion/Fermion.h @@ -237,4 +237,11 @@ typedef ImprovedStaggeredFermion5D ImprovedStaggeredFermion }} +//////////////////// +// Scalar QED actions +// TODO: this needs to move to another header after rename to Fermion.h +//////////////////// +#include +#include + #endif 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/action/gauge/Photon.h b/lib/qcd/action/gauge/Photon.h new file mode 100644 index 00000000..1512d4e3 --- /dev/null +++ b/lib/qcd/action/gauge/Photon.h @@ -0,0 +1,284 @@ +/************************************************************************************* + + Grid physics library, www.github.com/paboyle/Grid + + Source file: ./lib/qcd/action/gauge/Photon.h + + Copyright (C) 2015 + + Author: Peter Boyle + + 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 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 SiteLink; + typedef iImplGaugeField SiteField; + + typedef Lattice LinkField; + typedef Lattice Field; + }; + + typedef QedGimpl QedGimplR; + + template + class Photon + { + public: + INHERIT_GIMPL_TYPES(Gimpl); + 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; + 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); + private: + Gauge gauge_; + ZmScheme zmScheme_; + }; + + typedef Photon PhotonR; + + template + Photon::Photon(Gauge gauge, ZmScheme zmScheme) + : gauge_(gauge), zmScheme_(zmScheme) + {} + + template + void Photon::FreePropagator (const GaugeField &in,GaugeField &out) + { + FFT theFFT(in._grid); + + GaugeField in_k(in._grid); + GaugeField prop_k(in._grid); + + theFFT.FFT_all_dim(in_k,in,FFT::forward); + MomentumSpacePropagator(prop_k,in_k); + theFFT.FFT_all_dim(out,prop_k,FFT::backward); + } + + template + void Photon::invKHatSquared(GaugeLinkField &out) + { + GridBase *grid = out._grid; + GaugeLinkField kmu(grid), one(grid); + const unsigned int nd = grid->_ndimension; + std::vector &l = grid->_fdimensions; + std::vector zm(nd,0); + TComplex Tone = Complex(1.0,0.0); + TComplex Tzero= Complex(0.0,0.0); + + one = Complex(1.0,0.0); + out = zero; + 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; + } + pokeSite(Tone, out, zm); + out = one/out; + pokeSite(Tzero, out, zm); + } + + template + void Photon::zmSub(GaugeLinkField &out) + { + GridBase *grid = out._grid; + const unsigned int nd = grid->_ndimension; + + switch (zmScheme_) + { + case ZmScheme::qedTL: + { + std::vector zm(nd,0); + TComplex Tzero = Complex(0.0,0.0); + + pokeSite(Tzero, out, zm); + + break; + } + case ZmScheme::qedL: + { + LatticeInteger spNrm(grid), coor(grid); + GaugeLinkField z(grid); + + spNrm = zero; + for(int d = 0; d < grid->_ndimension - 1; d++) + { + LatticeCoordinate(coor,d); + spNrm = spNrm + coor*coor; + } + out = where(spNrm == Integer(0), 0.*out, out); + + break; + } + default: + break; + } + } + + template + void Photon::MomentumSpacePropagator(const GaugeField &in, + GaugeField &out) + { + GridBase *grid = out._grid; + LatticeComplex k2Inv(grid); + + invKHatSquared(k2Inv); + zmSub(k2Inv); + + out = in*k2Inv; + } + + template + void Photon::StochasticWeight(GaugeLinkField &weight) + { + 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(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 = weight*r; + pokeLorentz(aTilde, r, mu); + } + fft.FFT_all_dim(out, aTilde, FFT::backward); + + out = real(out); + } +// template +// void Photon::FeynmanGaugeMomentumSpacePropagator_L(GaugeField &out, +// const GaugeField &in) +// { +// +// FeynmanGaugeMomentumSpacePropagator_TL(out,in); +// +// GridBase *grid = out._grid; +// LatticeInteger coor(grid); +// GaugeField zz(grid); zz=zero; +// +// // xyzt +// for(int d = 0; d < grid->_ndimension-1;d++){ +// LatticeCoordinate(coor,d); +// out = where(coor==Integer(0),zz,out); +// } +// } +// +// template +// void Photon::FeynmanGaugeMomentumSpacePropagator_TL(GaugeField &out, +// const GaugeField &in) +// { +// +// // what type LatticeComplex +// GridBase *grid = out._grid; +// int nd = grid->_ndimension; +// +// typedef typename GaugeField::vector_type vector_type; +// typedef typename GaugeField::scalar_type ScalComplex; +// typedef Lattice > LatComplex; +// +// std::vector latt_size = grid->_fdimensions; +// +// LatComplex denom(grid); denom= zero; +// LatComplex one(grid); one = ScalComplex(1.0,0.0); +// LatComplex kmu(grid); +// +// ScalComplex ci(0.0,1.0); +// // momphase = n * 2pi / L +// for(int mu=0;mu zero_mode(nd,0); +// TComplexD Tone = ComplexD(1.0,0.0); +// TComplexD Tzero= ComplexD(0.0,0.0); +// +// pokeSite(Tone,denom,zero_mode); +// +// denom= one/denom; +// +// pokeSite(Tzero,denom,zero_mode); +// +// out = zero; +// out = in*denom; +// }; + +}} +#endif 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/action/scalar/ScalarImpl.h b/lib/qcd/action/scalar/ScalarImpl.h index ee2d2fb8..1b5fefea 100644 --- a/lib/qcd/action/scalar/ScalarImpl.h +++ b/lib/qcd/action/scalar/ScalarImpl.h @@ -14,9 +14,11 @@ namespace Grid { using iImplField = iScalar > >; typedef iImplField SiteField; - + typedef SiteField SitePropagator; typedef Lattice Field; + typedef Field FermionField; + typedef Field PropagatorField; static inline void generate_momenta(Field& P, GridParallelRNG& pRNG){ gaussian(pRNG, P); @@ -44,6 +46,45 @@ namespace Grid { U = 1.0; } + static void MomentumSpacePropagator(Field &out, RealD m) + { + GridBase *grid = out._grid; + Field 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++) + { + Real twoPiL = M_PI*2./l[mu]; + + LatticeCoordinate(kmu,mu); + kmu = 2.*sin(.5*twoPiL*kmu); + out = out + kmu*kmu; + } + out = one/out; + } + + static void FreePropagator(const Field &in, Field &out, + const Field &momKernel) + { + FFT fft((GridCartesian *)in._grid); + Field inFT(in._grid); + + fft.FFT_all_dim(inFT, in, FFT::forward); + inFT = inFT*momKernel; + fft.FFT_all_dim(out, inFT, FFT::backward); + } + + static void FreePropagator(const Field &in, Field &out, RealD m) + { + Field momKernel(in._grid); + + MomentumSpacePropagator(momKernel, m); + FreePropagator(in, out, momKernel); + } + }; template @@ -93,6 +134,9 @@ namespace Grid { typedef ScalarImplTypes ScalarImplR; typedef ScalarImplTypes ScalarImplF; typedef ScalarImplTypes ScalarImplD; + typedef ScalarImplTypes ScalarImplCR; + typedef ScalarImplTypes ScalarImplCF; + typedef ScalarImplTypes ScalarImplCD; //} } 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/qcd/smearing/WilsonFlow.h b/lib/qcd/smearing/WilsonFlow.h index 8b8f9a81..5e9f2d95 100644 --- a/lib/qcd/smearing/WilsonFlow.h +++ b/lib/qcd/smearing/WilsonFlow.h @@ -36,20 +36,23 @@ namespace QCD { template class WilsonFlow: public Smear{ unsigned int Nstep; - RealD epsilon; + unsigned int measure_interval; + mutable RealD epsilon, taus; + mutable WilsonGaugeAction SG; void evolve_step(typename Gimpl::GaugeField&) const; + void evolve_step_adaptive(typename Gimpl::GaugeField&, RealD); RealD tau(unsigned int t)const {return epsilon*(t+1.0); } - 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); @@ -72,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; }; @@ -98,23 +103,111 @@ 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, 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); + 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); 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 + + + 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(); + 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; + #ifdef WF_TIMING + std::cout << "Time to evolve " << diff.count() << " s\n"; + #endif 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; + } } } +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 diff --git a/lib/qcd/utils/WilsonLoops.h b/lib/qcd/utils/WilsonLoops.h index 5382882e..3c0e9ee4 100644 --- a/lib/qcd/utils/WilsonLoops.h +++ b/lib/qcd/utils/WilsonLoops.h @@ -188,6 +188,32 @@ public: } } + +// 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 + U[d] = PeekIndex(Umu, d); + } + staple = zero; + 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 ////////////////////////////////////////////////// @@ -200,7 +226,6 @@ public: U[d] = PeekIndex(Umu, d); } staple = zero; - GaugeMat tmp(grid); for (int nu = 0; nu < Nd; nu++) { @@ -214,7 +239,7 @@ public: // | // __| // - + staple += Gimpl::ShiftStaple( Gimpl::CovShiftForward( U[nu], nu, @@ -227,6 +252,7 @@ public: // |__ // // + staple += Gimpl::ShiftStaple( Gimpl::CovShiftBackward(U[nu], nu, Gimpl::CovShiftBackward(U[mu], mu, U[nu])), mu); @@ -289,8 +315,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 +332,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/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 } ////////////////////////////////////////////////////////////////////////////////////// diff --git a/lib/tensors/Tensor_exp.h b/lib/tensors/Tensor_exp.h index ad0d2071..e18fed70 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/RealD(i); + temp = unit + temp*arg; + } + return temp; + + } + + + + } #endif 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 3a642f24..a85beead 100644 --- a/tests/hadrons/Test_hadrons_rarekaon.cc +++ b/tests/hadrons/Test_hadrons_rarekaon.cc @@ -87,6 +87,10 @@ int main(int argc, char *argv[]) gaugePar.file = configStem; application.createModule(gaugeField, gaugePar); } + + // set fermion boundary conditions to be periodic space, antiperiodic time. + std::string boundary = "1 1 1 -1"; + for (unsigned int i = 0; i < flavour.size(); ++i) { // actions diff --git a/tests/hadrons/Test_hadrons_spectrum.cc b/tests/hadrons/Test_hadrons_spectrum.cc index 2d731ff4..8f7b30c8 100644 --- a/tests/hadrons/Test_hadrons_spectrum.cc +++ b/tests/hadrons/Test_hadrons_spectrum.cc @@ -63,6 +63,14 @@ 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"; + for (unsigned int i = 0; i < flavour.size(); ++i) { // actions @@ -71,6 +79,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 @@ -93,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); diff --git a/tests/smearing/Test_WilsonFlow.cc b/tests/smearing/Test_WilsonFlow.cc index 4e6bd0af..5db00d5d 100644 --- a/tests/smearing/Test_WilsonFlow.cc +++ b/tests/smearing/Test_WilsonFlow.cc @@ -28,6 +28,38 @@ directory /* END LEGAL */ #include +namespace Grid{ + struct WFParameters: Serializable { + GRID_SERIALIZABLE_CLASS_MEMBERS(WFParameters, + int, steps, + double, step_size, + int, meas_interval, + double, maxTau); // for the adaptive algorithm + + + 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; @@ -42,22 +74,38 @@ 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); + + 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); + + 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); + 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); - std::cout << GridLogMessage << "Plaquette: "<< WFlow_plaq << std::endl; + RealD WFlow_TC = WilsonLoops::TopologicalCharge(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; std::cout<< GridLogMessage << " Admissibility check:\n"; const double sp_adm = 0.067; // admissible threshold @@ -73,6 +121,32 @@ int main(int argc, char **argv) { std::cout<< GridLogMessage << " (sp_admissible = "<< sp_adm <<")\n"; //std::cout<< GridLogMessage << " sp_admissible - sp_max = "<