From 55e39df30f1953e0d54b79d2b1c6061093468e76 Mon Sep 17 00:00:00 2001 From: Vera Guelpers Date: Fri, 22 Dec 2017 11:36:31 +0000 Subject: [PATCH 01/38] tadpole insertion for DWF --- lib/qcd/action/fermion/WilsonFermion5D.cc | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/lib/qcd/action/fermion/WilsonFermion5D.cc b/lib/qcd/action/fermion/WilsonFermion5D.cc index 393ee7f3..bc703187 100644 --- a/lib/qcd/action/fermion/WilsonFermion5D.cc +++ b/lib/qcd/action/fermion/WilsonFermion5D.cc @@ -833,9 +833,12 @@ void WilsonFermion5D::SeqConservedCurrent(PropagatorField &q_in, for (unsigned int s = 0; s < LLs; ++s) { bool axial_sign = ((curr_type == Current::Axial) && (s < (LLs / 2))); + bool tadpole_sign = (curr_type == Current::Tadpole); + bool switch_sgn = tadpole_sign || axial_sign; + Kernels::SeqConservedCurrentSiteFwd(tmpFwd._odata[sF], q_out._odata[sF], Umu, sU, - mu, t_mask, axial_sign); + mu, t_mask, switch_sgn); ++sF; } } From 935cd1e173e6f50f5185efb30e6b9c3b92790c3d Mon Sep 17 00:00:00 2001 From: Vera Guelpers Date: Fri, 22 Dec 2017 11:38:45 +0000 Subject: [PATCH 02/38] conserved current insertion summed over Lorentzindex --- extras/Hadrons/Modules.hpp | 67 ++----- .../Modules/MSource/SeqConservedSummed.hpp | 165 ++++++++++++++++++ extras/Hadrons/modules.inc | 46 ++--- 3 files changed, 208 insertions(+), 70 deletions(-) create mode 100644 extras/Hadrons/Modules/MSource/SeqConservedSummed.hpp diff --git a/extras/Hadrons/Modules.hpp b/extras/Hadrons/Modules.hpp index e1f06f32..a52b9a6e 100644 --- a/extras/Hadrons/Modules.hpp +++ b/extras/Hadrons/Modules.hpp @@ -1,60 +1,31 @@ -/************************************************************************************* - -Grid physics library, www.github.com/paboyle/Grid - -Source file: extras/Hadrons/Modules.hpp - -Copyright (C) 2015 -Copyright (C) 2016 -Copyright (C) 2017 - -Author: Antonin Portelli - -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 -#include -#include -#include +#include #include -#include +#include +#include +#include #include -#include +#include +#include #include -#include -#include -#include -#include -#include +#include +#include +#include #include -#include #include +#include #include #include #include -#include +#include +#include +#include #include #include #include -#include #include -#include -#include +#include +#include +#include +#include +#include +#include diff --git a/extras/Hadrons/Modules/MSource/SeqConservedSummed.hpp b/extras/Hadrons/Modules/MSource/SeqConservedSummed.hpp new file mode 100644 index 00000000..2296701a --- /dev/null +++ b/extras/Hadrons/Modules/MSource/SeqConservedSummed.hpp @@ -0,0 +1,165 @@ +/************************************************************************************* + +Grid physics library, www.github.com/paboyle/Grid + +Source file: extras/Hadrons/Modules/MContraction/SeqConservedSummed.hpp + +Copyright (C) 2017 + +Author: Andrew Lawson +Author: Vera Guelpers + +This program is free software; you can redistribute it and/or modify +it under the terms of the GNU General Public License as published by +the Free Software Foundation; either version 2 of the License, or +(at your option) any later version. + +This program is distributed in the hope that it will be useful, +but WITHOUT ANY WARRANTY; without even the implied warranty of +MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +GNU General Public License for more details. + +You should have received a copy of the GNU General Public License along +with this program; if not, write to the Free Software Foundation, Inc., +51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. + +See the full license in the file "LICENSE" in the top level distribution directory +*************************************************************************************/ +/* END LEGAL */ + +#ifndef Hadrons_MSource_SeqConservedSummed_hpp_ +#define Hadrons_MSource_SeqConservedSummed_hpp_ + +#include +#include +#include + +BEGIN_HADRONS_NAMESPACE + +/* + + Sequential source summed over the Lorentz index of current + ----------------------------- + * src_x = sum_mu q_x * theta(x_3 - tA) * theta(tB - x_3) * J_mu * exp(i x.mom) + + * options: + - q: input propagator (string) + - action: fermion action used for propagator q (string) + - tA: begin timeslice (integer) + - tB: end timesilce (integer) + - curr_type: type of conserved current to insert (Current) + - mom: momentum insertion, space-separated float sequence (e.g ".1 .2 1. 0.") + + */ + +/****************************************************************************** + * SeqConservedSummed * + ******************************************************************************/ +BEGIN_MODULE_NAMESPACE(MSource) + +class SeqConservedSummedPar: Serializable +{ +public: + GRID_SERIALIZABLE_CLASS_MEMBERS(SeqConservedSummedPar, + std::string, q, + std::string, action, + unsigned int, tA, + unsigned int, tB, + Current, curr_type, + std::string, mom); +}; + +template +class TSeqConservedSummed: public Module +{ +public: + FERM_TYPE_ALIASES(FImpl,); +public: + // constructor + TSeqConservedSummed(const std::string name); + // destructor + virtual ~TSeqConservedSummed(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(SeqConservedSummed, TSeqConservedSummed, MSource); + +/****************************************************************************** + * TSeqConservedSummed implementation * + ******************************************************************************/ +// constructor ///////////////////////////////////////////////////////////////// +template +TSeqConservedSummed::TSeqConservedSummed(const std::string name) +: Module(name) +{} + +// dependencies/products /////////////////////////////////////////////////////// +template +std::vector TSeqConservedSummed::getInput(void) +{ + std::vector in = {par().q, par().action}; + + return in; +} + +template +std::vector TSeqConservedSummed::getOutput(void) +{ + std::vector out = {getName()}; + + return out; +} + +// setup /////////////////////////////////////////////////////////////////////// +template +void TSeqConservedSummed::setup(void) +{ + auto Ls_ = env().getObjectLs(par().action); + env().template registerLattice(getName(), Ls_); +} + +// execution /////////////////////////////////////////////////////////////////// +template +void TSeqConservedSummed::execute(void) +{ + if (par().tA == par().tB) + { + LOG(Message) << "Generating sequential source with conserved " + << par().curr_type << " current insertion summed over mu at " + << "t = " << par().tA << std::endl; + } + else + { + LOG(Message) << "Generating sequential source with conserved " + << par().curr_type << " current insertion summed over mu for " + << par().tA << " <= t <= " + << par().tB << std::endl; + } + PropagatorField &src = *env().template createLattice(getName()); + auto src_buf = src; + PropagatorField &q = *env().template getObject(par().q); + FMat &mat = *(env().template getObject(par().action)); + + std::vector mom = strToVec(par().mom); + src = zero; + for(int mu=0;mu<=3;mu++) + { + mat.SeqConservedCurrent(q, src_buf, par().curr_type, mu, + mom, par().tA, par().tB); + src += src_buf; + + } + +} + +END_MODULE_NAMESPACE + +END_HADRONS_NAMESPACE + +#endif // Hadrons_SeqConservedSummed_hpp_ diff --git a/extras/Hadrons/modules.inc b/extras/Hadrons/modules.inc index fbbb2eb9..1ad4750a 100644 --- a/extras/Hadrons/modules.inc +++ b/extras/Hadrons/modules.inc @@ -1,42 +1,44 @@ modules_cc =\ Modules/MContraction/WeakHamiltonianEye.cc \ - Modules/MContraction/WeakHamiltonianNonEye.cc \ Modules/MContraction/WeakNeutral4ptDisc.cc \ - Modules/MGauge/Load.cc \ + Modules/MContraction/WeakHamiltonianNonEye.cc \ + Modules/MScalar/FreeProp.cc \ + Modules/MScalar/ChargedProp.cc \ Modules/MGauge/Random.cc \ Modules/MGauge/StochEm.cc \ Modules/MGauge/Unit.cc \ - Modules/MScalar/ChargedProp.cc \ - Modules/MScalar/FreeProp.cc + Modules/MGauge/Load.cc modules_hpp =\ - Modules/MAction/DWF.hpp \ - Modules/MAction/Wilson.hpp \ - Modules/MContraction/Baryon.hpp \ - Modules/MContraction/DiscLoop.hpp \ - Modules/MContraction/Gamma3pt.hpp \ - Modules/MContraction/Meson.hpp \ + Modules/MSolver/RBPrecCG.hpp \ Modules/MContraction/WardIdentity.hpp \ - Modules/MContraction/WeakHamiltonian.hpp \ + Modules/MContraction/Meson.hpp \ + Modules/MContraction/Gamma3pt.hpp \ + Modules/MContraction/DiscLoop.hpp \ Modules/MContraction/WeakHamiltonianEye.hpp \ - Modules/MContraction/WeakHamiltonianNonEye.hpp \ + Modules/MContraction/Baryon.hpp \ + Modules/MContraction/WeakHamiltonian.hpp \ Modules/MContraction/WeakNeutral4ptDisc.hpp \ - Modules/MFermion/GaugeProp.hpp \ - Modules/MGauge/Load.hpp \ - Modules/MGauge/Random.hpp \ - Modules/MGauge/StochEm.hpp \ - Modules/MGauge/Unit.hpp \ + Modules/MContraction/WeakHamiltonianNonEye.hpp \ + Modules/MUtilities/TestSeqGamma.hpp \ + Modules/MUtilities/TestSeqConserved.hpp \ Modules/MLoop/NoiseLoop.hpp \ - Modules/MScalar/ChargedProp.hpp \ Modules/MScalar/FreeProp.hpp \ + Modules/MScalar/ChargedProp.hpp \ Modules/MScalar/Scalar.hpp \ Modules/MSink/Point.hpp \ Modules/MSink/Smear.hpp \ - Modules/MSolver/RBPrecCG.hpp \ + Modules/MFermion/GaugeProp.hpp \ + Modules/MSource/SeqConservedSummed.hpp \ + Modules/MSource/Wall.hpp \ Modules/MSource/Point.hpp \ Modules/MSource/SeqConserved.hpp \ Modules/MSource/SeqGamma.hpp \ - Modules/MSource/Wall.hpp \ Modules/MSource/Z2.hpp \ - Modules/MUtilities/TestSeqConserved.hpp \ - Modules/MUtilities/TestSeqGamma.hpp + Modules/MGauge/Load.hpp \ + Modules/MGauge/StochEm.hpp \ + Modules/MGauge/Random.hpp \ + Modules/MGauge/Unit.hpp \ + Modules/MAction/DWF.hpp \ + Modules/MAction/Wilson.hpp + From 389731d373e5f348c80a29520a4328414470c0a2 Mon Sep 17 00:00:00 2001 From: Vera Guelpers Date: Tue, 23 Jan 2018 10:11:33 +0000 Subject: [PATCH 03/38] changed SeqConservedSummed.hpp to work with new hadrons interface --- .../Modules/MSource/SeqConservedSummed.hpp | 20 +++++++++++-------- 1 file changed, 12 insertions(+), 8 deletions(-) diff --git a/extras/Hadrons/Modules/MSource/SeqConservedSummed.hpp b/extras/Hadrons/Modules/MSource/SeqConservedSummed.hpp index 2296701a..243d05f2 100644 --- a/extras/Hadrons/Modules/MSource/SeqConservedSummed.hpp +++ b/extras/Hadrons/Modules/MSource/SeqConservedSummed.hpp @@ -4,8 +4,9 @@ Grid physics library, www.github.com/paboyle/Grid Source file: extras/Hadrons/Modules/MContraction/SeqConservedSummed.hpp -Copyright (C) 2017 +Copyright (C) 2015-2018 +Author: Antonin Portelli Author: Andrew Lawson Author: Vera Guelpers @@ -82,6 +83,7 @@ public: // dependency relation virtual std::vector getInput(void); virtual std::vector getOutput(void); +protected: // setup virtual void setup(void); // execution @@ -121,7 +123,8 @@ template void TSeqConservedSummed::setup(void) { auto Ls_ = env().getObjectLs(par().action); - env().template registerLattice(getName(), Ls_); + envCreateLat(PropagatorField, getName(), Ls_); + envTmpLat(PropagatorField, "src_tmp"); } // execution /////////////////////////////////////////////////////////////////// @@ -141,18 +144,19 @@ void TSeqConservedSummed::execute(void) << par().tA << " <= t <= " << par().tB << std::endl; } - PropagatorField &src = *env().template createLattice(getName()); - auto src_buf = src; - PropagatorField &q = *env().template getObject(par().q); - FMat &mat = *(env().template getObject(par().action)); + auto &src = envGet(PropagatorField, getName()); + envGetTmp(PropagatorField, src_tmp); + src_tmp = src; + auto &q = envGet(PropagatorField, par().q); + auto &mat = envGet(FMat, par().action); std::vector mom = strToVec(par().mom); src = zero; for(int mu=0;mu<=3;mu++) { - mat.SeqConservedCurrent(q, src_buf, par().curr_type, mu, + mat.SeqConservedCurrent(q, src_tmp, par().curr_type, mu, mom, par().tA, par().tB); - src += src_buf; + src += src_tmp; } From b6fe03eb26a8725d4f52d3215104ea9743c071b0 Mon Sep 17 00:00:00 2001 From: Vera Guelpers Date: Fri, 2 Feb 2018 15:29:38 +0000 Subject: [PATCH 04/38] BugFix: Now the stochatic EM potential weight is generated when calling for the first time --- extras/Hadrons/Modules/MGauge/StochEm.cc | 6 ++++-- extras/Hadrons/Modules/MGauge/StochEm.hpp | 2 ++ 2 files changed, 6 insertions(+), 2 deletions(-) diff --git a/extras/Hadrons/Modules/MGauge/StochEm.cc b/extras/Hadrons/Modules/MGauge/StochEm.cc index 21b7f626..cfbd28ef 100644 --- a/extras/Hadrons/Modules/MGauge/StochEm.cc +++ b/extras/Hadrons/Modules/MGauge/StochEm.cc @@ -57,9 +57,11 @@ std::vector TStochEm::getOutput(void) // setup /////////////////////////////////////////////////////////////////////// void TStochEm::setup(void) { + create_weight = false; if (!env().hasCreatedObject("_" + getName() + "_weight")) { envCacheLat(EmComp, "_" + getName() + "_weight"); + create_weight = true; } envCreateLat(EmField, getName()); } @@ -67,13 +69,13 @@ void TStochEm::setup(void) // execution /////////////////////////////////////////////////////////////////// void TStochEm::execute(void) { - LOG(Message) << "Generating stochatic EM potential..." << std::endl; + LOG(Message) << "Generating stochastic EM potential..." << std::endl; PhotonR photon(par().gauge, par().zmScheme); auto &a = envGet(EmField, getName()); auto &w = envGet(EmComp, "_" + getName() + "_weight"); - if (!env().hasCreatedObject("_" + getName() + "_weight")) + if (create_weight) { LOG(Message) << "Caching stochatic EM potential weight (gauge: " << par().gauge << ", zero-mode scheme: " diff --git a/extras/Hadrons/Modules/MGauge/StochEm.hpp b/extras/Hadrons/Modules/MGauge/StochEm.hpp index 87b70880..05162443 100644 --- a/extras/Hadrons/Modules/MGauge/StochEm.hpp +++ b/extras/Hadrons/Modules/MGauge/StochEm.hpp @@ -60,6 +60,8 @@ public: // dependency relation virtual std::vector getInput(void); virtual std::vector getOutput(void); +private: + bool create_weight; protected: // setup virtual void setup(void); From 49a0ae73eb097fc9e7de02897bcde6ff6a3db170 Mon Sep 17 00:00:00 2001 From: Vera Guelpers Date: Mon, 12 Feb 2018 09:36:08 +0000 Subject: [PATCH 05/38] Insertion of photon field in seqential conserved current --- extras/Hadrons/Modules.hpp | 1 - .../Hadrons/Modules/MSource/SeqConserved.hpp | 93 ++++++++-- .../Modules/MSource/SeqConservedSummed.hpp | 169 ------------------ extras/Hadrons/modules.inc | 1 - lib/qcd/action/fermion/FermionOperator.h | 4 +- .../fermion/ImprovedStaggeredFermion.cc | 16 +- .../action/fermion/ImprovedStaggeredFermion.h | 12 +- .../fermion/ImprovedStaggeredFermion5D.cc | 15 +- .../fermion/ImprovedStaggeredFermion5D.h | 12 +- lib/qcd/action/fermion/WilsonFermion.cc | 90 ++++++++-- lib/qcd/action/fermion/WilsonFermion.h | 14 +- lib/qcd/action/fermion/WilsonFermion5D.cc | 127 +++++++++++-- lib/qcd/action/fermion/WilsonFermion5D.h | 10 +- tests/hadrons/Test_hadrons.hpp | 3 +- 14 files changed, 310 insertions(+), 257 deletions(-) delete mode 100644 extras/Hadrons/Modules/MSource/SeqConservedSummed.hpp diff --git a/extras/Hadrons/Modules.hpp b/extras/Hadrons/Modules.hpp index ed6254f4..734ff0fc 100644 --- a/extras/Hadrons/Modules.hpp +++ b/extras/Hadrons/Modules.hpp @@ -28,7 +28,6 @@ See the full license in the file "LICENSE" in the top level distribution directo #include #include #include -#include #include #include #include diff --git a/extras/Hadrons/Modules/MSource/SeqConserved.hpp b/extras/Hadrons/Modules/MSource/SeqConserved.hpp index ee8d8d56..0fb23cc4 100644 --- a/extras/Hadrons/Modules/MSource/SeqConserved.hpp +++ b/extras/Hadrons/Modules/MSource/SeqConserved.hpp @@ -2,12 +2,13 @@ Grid physics library, www.github.com/paboyle/Grid -Source file: extras/Hadrons/Modules/MSource/SeqConserved.hpp +Source file: extras/Hadrons/Modules/MContraction/SeqConserved.hpp Copyright (C) 2015-2018 Author: Antonin Portelli -Author: Lanny91 +Author: Andrew Lawson +Author: Vera Guelpers 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 @@ -34,13 +35,17 @@ See the full license in the file "LICENSE" in the top level distribution directo #include #include +#include + BEGIN_HADRONS_NAMESPACE /* - Sequential source + Sequential source with insertion of conserved current. + Additionally optional insertion of a photon field A_\mu(x). ----------------------------- - * src_x = q_x * theta(x_3 - tA) * theta(tB - x_3) * J_mu * exp(i x.mom) + * src_x = sum_{mu=mu_min}^{mu_max} + q_x * theta(x_3 - tA) * theta(tB - x_3) * J_mu * exp(i x.mom) (* A_\mu(x)) * options: - q: input propagator (string) @@ -48,8 +53,10 @@ BEGIN_HADRONS_NAMESPACE - tA: begin timeslice (integer) - tB: end timesilce (integer) - curr_type: type of conserved current to insert (Current) - - mu: Lorentz index of current to insert (integer) + - mu_min: begin Lorentz Index (integer) + - mu_max: end Lorentz Index (integer) - mom: momentum insertion, space-separated float sequence (e.g ".1 .2 1. 0.") + - photon: optional photon field (string) */ @@ -67,8 +74,10 @@ public: unsigned int, tA, unsigned int, tB, Current, curr_type, - unsigned int, mu, - std::string, mom); + unsigned int, mu_min, + unsigned int, mu_max, + std::string, mom, + std::string, photon); }; template @@ -76,6 +85,8 @@ class TSeqConserved: public Module { public: FERM_TYPE_ALIASES(FImpl,); +public: + typedef PhotonR::GaugeField EmField; public: // constructor TSeqConserved(const std::string name); @@ -93,6 +104,7 @@ protected: MODULE_REGISTER_NS(SeqConserved, TSeqConserved, MSource); + /****************************************************************************** * TSeqConserved implementation * ******************************************************************************/ @@ -107,7 +119,8 @@ template std::vector TSeqConserved::getInput(void) { std::vector in = {par().q, par().action}; - + if (!par().photon.empty()) in.push_back(par().photon); + return in; } @@ -116,7 +129,7 @@ std::vector TSeqConserved::getOutput(void) { std::vector out = {getName()}; - return out; + return out; } // setup /////////////////////////////////////////////////////////////////////// @@ -125,6 +138,10 @@ void TSeqConserved::setup(void) { auto Ls_ = env().getObjectLs(par().action); envCreateLat(PropagatorField, getName(), Ls_); + envTmpLat(PropagatorField, "src_tmp"); + envTmpLat(LatticeComplex, "mom_phase"); + envTmpLat(LatticeComplex, "coor"); + envTmpLat(LatticeComplex, "latt_compl"); } // execution /////////////////////////////////////////////////////////////////// @@ -134,27 +151,71 @@ void TSeqConserved::execute(void) if (par().tA == par().tB) { LOG(Message) << "Generating sequential source with conserved " - << par().curr_type << " current insertion (mu = " - << par().mu << ") at " << "t = " << par().tA << std::endl; + << par().curr_type << " current at " + << "t = " << par().tA << std::endl; } else { LOG(Message) << "Generating sequential source with conserved " - << par().curr_type << " current insertion (mu = " - << par().mu << ") for " << par().tA << " <= t <= " + << par().curr_type << " current for " + << par().tA << " <= t <= " << par().tB << std::endl; } auto &src = envGet(PropagatorField, getName()); + envGetTmp(PropagatorField, src_tmp); + src_tmp = src; auto &q = envGet(PropagatorField, par().q); auto &mat = envGet(FMat, par().action); + envGetTmp(LatticeComplex, mom_phase); + envGetTmp(LatticeComplex, coor); + envGetTmp(LatticeComplex, latt_compl); + src = zero; + + //exp(ipx) std::vector mom = strToVec(par().mom); - mat.SeqConservedCurrent(q, src, par().curr_type, par().mu, - mom, par().tA, par().tB); + mom_phase = zero; + Complex i(0.0,1.0); + for(unsigned int mu = 0; mu < env().getNd(); mu++) + { + LatticeCoordinate(coor, mu); + mom_phase = mom_phase + (mom[mu]/env().getGrid()->_fdimensions[mu])*coor; + } + mom_phase = exp((Real)(2*M_PI)*i*mom_phase); + LOG(Message) << "Inserting momentum " << mom << std::endl; + + + + if (!par().photon.empty()) + { + LOG(Message) << "Inserting the stochastic photon field " << par().photon << std::endl; + } + + for(unsigned int mu=par().mu_min;mu<=par().mu_max;mu++) + { + if (!par().photon.empty()) + { + //Get the stochastic photon field, if required + auto &stoch_photon = envGet(EmField, par().photon); + latt_compl = PeekIndex(stoch_photon, mu) * mom_phase; + } + else + { + latt_compl = mom_phase; + } + + mat.SeqConservedCurrent(q, src_tmp, par().curr_type, mu, + par().tA, par().tB, latt_compl); + src += src_tmp; + + } + + } + END_MODULE_NAMESPACE END_HADRONS_NAMESPACE -#endif // Hadrons_SeqConserved_hpp_ +#endif // Hadrons_MSource_SeqConserved_hpp_ diff --git a/extras/Hadrons/Modules/MSource/SeqConservedSummed.hpp b/extras/Hadrons/Modules/MSource/SeqConservedSummed.hpp deleted file mode 100644 index 243d05f2..00000000 --- a/extras/Hadrons/Modules/MSource/SeqConservedSummed.hpp +++ /dev/null @@ -1,169 +0,0 @@ -/************************************************************************************* - -Grid physics library, www.github.com/paboyle/Grid - -Source file: extras/Hadrons/Modules/MContraction/SeqConservedSummed.hpp - -Copyright (C) 2015-2018 - -Author: Antonin Portelli -Author: Andrew Lawson -Author: Vera Guelpers - -This program is free software; you can redistribute it and/or modify -it under the terms of the GNU General Public License as published by -the Free Software Foundation; either version 2 of the License, or -(at your option) any later version. - -This program is distributed in the hope that it will be useful, -but WITHOUT ANY WARRANTY; without even the implied warranty of -MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the -GNU General Public License for more details. - -You should have received a copy of the GNU General Public License along -with this program; if not, write to the Free Software Foundation, Inc., -51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. - -See the full license in the file "LICENSE" in the top level distribution directory -*************************************************************************************/ -/* END LEGAL */ - -#ifndef Hadrons_MSource_SeqConservedSummed_hpp_ -#define Hadrons_MSource_SeqConservedSummed_hpp_ - -#include -#include -#include - -BEGIN_HADRONS_NAMESPACE - -/* - - Sequential source summed over the Lorentz index of current - ----------------------------- - * src_x = sum_mu q_x * theta(x_3 - tA) * theta(tB - x_3) * J_mu * exp(i x.mom) - - * options: - - q: input propagator (string) - - action: fermion action used for propagator q (string) - - tA: begin timeslice (integer) - - tB: end timesilce (integer) - - curr_type: type of conserved current to insert (Current) - - mom: momentum insertion, space-separated float sequence (e.g ".1 .2 1. 0.") - - */ - -/****************************************************************************** - * SeqConservedSummed * - ******************************************************************************/ -BEGIN_MODULE_NAMESPACE(MSource) - -class SeqConservedSummedPar: Serializable -{ -public: - GRID_SERIALIZABLE_CLASS_MEMBERS(SeqConservedSummedPar, - std::string, q, - std::string, action, - unsigned int, tA, - unsigned int, tB, - Current, curr_type, - std::string, mom); -}; - -template -class TSeqConservedSummed: public Module -{ -public: - FERM_TYPE_ALIASES(FImpl,); -public: - // constructor - TSeqConservedSummed(const std::string name); - // destructor - virtual ~TSeqConservedSummed(void) = default; - // dependency relation - virtual std::vector getInput(void); - virtual std::vector getOutput(void); -protected: - // setup - virtual void setup(void); - // execution - virtual void execute(void); -}; - -MODULE_REGISTER_NS(SeqConservedSummed, TSeqConservedSummed, MSource); - -/****************************************************************************** - * TSeqConservedSummed implementation * - ******************************************************************************/ -// constructor ///////////////////////////////////////////////////////////////// -template -TSeqConservedSummed::TSeqConservedSummed(const std::string name) -: Module(name) -{} - -// dependencies/products /////////////////////////////////////////////////////// -template -std::vector TSeqConservedSummed::getInput(void) -{ - std::vector in = {par().q, par().action}; - - return in; -} - -template -std::vector TSeqConservedSummed::getOutput(void) -{ - std::vector out = {getName()}; - - return out; -} - -// setup /////////////////////////////////////////////////////////////////////// -template -void TSeqConservedSummed::setup(void) -{ - auto Ls_ = env().getObjectLs(par().action); - envCreateLat(PropagatorField, getName(), Ls_); - envTmpLat(PropagatorField, "src_tmp"); -} - -// execution /////////////////////////////////////////////////////////////////// -template -void TSeqConservedSummed::execute(void) -{ - if (par().tA == par().tB) - { - LOG(Message) << "Generating sequential source with conserved " - << par().curr_type << " current insertion summed over mu at " - << "t = " << par().tA << std::endl; - } - else - { - LOG(Message) << "Generating sequential source with conserved " - << par().curr_type << " current insertion summed over mu for " - << par().tA << " <= t <= " - << par().tB << std::endl; - } - auto &src = envGet(PropagatorField, getName()); - envGetTmp(PropagatorField, src_tmp); - src_tmp = src; - auto &q = envGet(PropagatorField, par().q); - auto &mat = envGet(FMat, par().action); - - std::vector mom = strToVec(par().mom); - src = zero; - for(int mu=0;mu<=3;mu++) - { - mat.SeqConservedCurrent(q, src_tmp, par().curr_type, mu, - mom, par().tA, par().tB); - src += src_tmp; - - } - -} - -END_MODULE_NAMESPACE - -END_HADRONS_NAMESPACE - -#endif // Hadrons_SeqConservedSummed_hpp_ diff --git a/extras/Hadrons/modules.inc b/extras/Hadrons/modules.inc index d005caec..63b674d4 100644 --- a/extras/Hadrons/modules.inc +++ b/extras/Hadrons/modules.inc @@ -20,7 +20,6 @@ modules_hpp =\ Modules/MContraction/WardIdentity.hpp \ Modules/MContraction/WeakHamiltonianEye.hpp \ Modules/MFermion/GaugeProp.hpp \ - Modules/MSource/SeqConservedSummed.hpp \ Modules/MSource/SeqGamma.hpp \ Modules/MSource/Point.hpp \ Modules/MSource/Wall.hpp \ diff --git a/lib/qcd/action/fermion/FermionOperator.h b/lib/qcd/action/fermion/FermionOperator.h index 1d395d53..94e065cc 100644 --- a/lib/qcd/action/fermion/FermionOperator.h +++ b/lib/qcd/action/fermion/FermionOperator.h @@ -125,9 +125,9 @@ namespace Grid { PropagatorField &q_out, Current curr_type, unsigned int mu, - std::vector mom, unsigned int tmin, - unsigned int tmax)=0; + unsigned int tmax, + Lattice> &lattice_cmplx)=0; }; } diff --git a/lib/qcd/action/fermion/ImprovedStaggeredFermion.cc b/lib/qcd/action/fermion/ImprovedStaggeredFermion.cc index 5005b1a2..3a296e52 100644 --- a/lib/qcd/action/fermion/ImprovedStaggeredFermion.cc +++ b/lib/qcd/action/fermion/ImprovedStaggeredFermion.cc @@ -407,17 +407,19 @@ void ImprovedStaggeredFermion::ContractConservedCurrent(PropagatorField &q } template -void ImprovedStaggeredFermion::SeqConservedCurrent(PropagatorField &q_in, - PropagatorField &q_out, - Current curr_type, - unsigned int mu, - std::vector mom, - unsigned int tmin, - unsigned int tmax) +void ImprovedStaggeredFermion::SeqConservedCurrent(PropagatorField &q_in, + PropagatorField &q_out, + Current curr_type, + unsigned int mu, + unsigned int tmin, + unsigned int tmax, + Lattice> &lattice_cmplx) { assert(0); + } + FermOpStaggeredTemplateInstantiate(ImprovedStaggeredFermion); //AdjointFermOpTemplateInstantiate(ImprovedStaggeredFermion); diff --git a/lib/qcd/action/fermion/ImprovedStaggeredFermion.h b/lib/qcd/action/fermion/ImprovedStaggeredFermion.h index 9d5270c6..015eb6bb 100644 --- a/lib/qcd/action/fermion/ImprovedStaggeredFermion.h +++ b/lib/qcd/action/fermion/ImprovedStaggeredFermion.h @@ -166,13 +166,13 @@ class ImprovedStaggeredFermion : public StaggeredKernels, public ImprovedS PropagatorField &q_out, Current curr_type, unsigned int mu); - void SeqConservedCurrent(PropagatorField &q_in, + void SeqConservedCurrent(PropagatorField &q_in, PropagatorField &q_out, - Current curr_type, - unsigned int mu, - std::vector mom, - unsigned int tmin, - unsigned int tmax); + Current curr_type, + unsigned int mu, + unsigned int tmin, + unsigned int tmax, + Lattice> &lattice_cmplx); }; typedef ImprovedStaggeredFermion ImprovedStaggeredFermionF; diff --git a/lib/qcd/action/fermion/ImprovedStaggeredFermion5D.cc b/lib/qcd/action/fermion/ImprovedStaggeredFermion5D.cc index 7308ff30..3640a26b 100644 --- a/lib/qcd/action/fermion/ImprovedStaggeredFermion5D.cc +++ b/lib/qcd/action/fermion/ImprovedStaggeredFermion5D.cc @@ -419,15 +419,16 @@ void ImprovedStaggeredFermion5D::ContractConservedCurrent(PropagatorField } template -void ImprovedStaggeredFermion5D::SeqConservedCurrent(PropagatorField &q_in, - PropagatorField &q_out, - Current curr_type, - unsigned int mu, - std::vector mom, - unsigned int tmin, - unsigned int tmax) +void ImprovedStaggeredFermion5D::SeqConservedCurrent(PropagatorField &q_in, + PropagatorField &q_out, + Current curr_type, + unsigned int mu, + unsigned int tmin, + unsigned int tmax, + Lattice> &lattice_cmplx) { assert(0); + } FermOpStaggeredTemplateInstantiate(ImprovedStaggeredFermion5D); diff --git a/lib/qcd/action/fermion/ImprovedStaggeredFermion5D.h b/lib/qcd/action/fermion/ImprovedStaggeredFermion5D.h index 22fb9e7d..1c211938 100644 --- a/lib/qcd/action/fermion/ImprovedStaggeredFermion5D.h +++ b/lib/qcd/action/fermion/ImprovedStaggeredFermion5D.h @@ -178,13 +178,13 @@ namespace QCD { PropagatorField &q_out, Current curr_type, unsigned int mu); - void SeqConservedCurrent(PropagatorField &q_in, + void SeqConservedCurrent(PropagatorField &q_in, PropagatorField &q_out, - Current curr_type, - unsigned int mu, - std::vector mom, - unsigned int tmin, - unsigned int tmax); + Current curr_type, + unsigned int mu, + unsigned int tmin, + unsigned int tmax, + Lattice> &lattice_cmplx); }; }} diff --git a/lib/qcd/action/fermion/WilsonFermion.cc b/lib/qcd/action/fermion/WilsonFermion.cc index 1a020e8a..560c371d 100644 --- a/lib/qcd/action/fermion/WilsonFermion.cc +++ b/lib/qcd/action/fermion/WilsonFermion.cc @@ -381,40 +381,100 @@ void WilsonFermion::ContractConservedCurrent(PropagatorField &q_in_1, } } +//template +//void WilsonFermion::SeqConservedCurrent(PropagatorField &q_in, +// PropagatorField &q_out, +// Current curr_type, +// unsigned int mu, +// std::vector mom, +// unsigned int tmin, +// unsigned int tmax) +//{ +// conformable(_grid, q_in._grid); +// conformable(_grid, q_out._grid); +// Lattice> ph(_grid), coor(_grid); +// Complex i(0.0,1.0); +// PropagatorField tmpFwd(_grid), tmpBwd(_grid), tmp(_grid); +// unsigned int tshift = (mu == Tp) ? 1 : 0; +// unsigned int LLt = GridDefaultLatt()[Tp]; +// +// // Momentum projection +// ph = zero; +// for(unsigned int mu = 0; mu < Nd - 1; mu++) +// { +// LatticeCoordinate(coor, mu); +// ph = ph + mom[mu]*coor*((1./(_grid->_fdimensions[mu]))); +// } +// ph = exp((Real)(2*M_PI)*i*ph); +// +// q_out = zero; +// LatticeInteger coords(_grid); +// LatticeCoordinate(coords, Tp); +// +// // Need q(x + mu) and q(x - mu). +// tmp = Cshift(q_in, mu, 1); +// tmpFwd = tmp*ph; +// tmp = ph*q_in; +// tmpBwd = Cshift(tmp, mu, -1); +// +// parallel_for (unsigned int sU = 0; sU < Umu._grid->oSites(); ++sU) +// { +// // Compute the sequential conserved current insertion only if our simd +// // object contains a timeslice we need. +// vInteger t_mask = ((coords._odata[sU] >= tmin) && +// (coords._odata[sU] <= tmax)); +// Integer timeSlices = Reduce(t_mask); +// +// if (timeSlices > 0) +// { +// Kernels::SeqConservedCurrentSiteFwd(tmpFwd._odata[sU], +// q_out._odata[sU], +// Umu, sU, mu, t_mask); +// } +// +// // Repeat for backward direction. +// t_mask = ((coords._odata[sU] >= (tmin + tshift)) && +// (coords._odata[sU] <= (tmax + tshift))); +// +// //if tmax = LLt-1 (last timeslice) include timeslice 0 if the time is shifted (mu=3) +// unsigned int t0 = 0; +// if((tmax==LLt-1) && (tshift==1)) t_mask = (t_mask || (coords._odata[sU] == t0 )); +// +// timeSlices = Reduce(t_mask); +// +// if (timeSlices > 0) +// { +// Kernels::SeqConservedCurrentSiteBwd(tmpBwd._odata[sU], +// q_out._odata[sU], +// Umu, sU, mu, t_mask); +// } +// } +//} + template void WilsonFermion::SeqConservedCurrent(PropagatorField &q_in, PropagatorField &q_out, Current curr_type, unsigned int mu, - std::vector mom, unsigned int tmin, - unsigned int tmax) + unsigned int tmax, + Lattice> &lattice_cmplx) { conformable(_grid, q_in._grid); conformable(_grid, q_out._grid); - Lattice> ph(_grid), coor(_grid); Complex i(0.0,1.0); PropagatorField tmpFwd(_grid), tmpBwd(_grid), tmp(_grid); unsigned int tshift = (mu == Tp) ? 1 : 0; unsigned int LLt = GridDefaultLatt()[Tp]; - // Momentum projection - ph = zero; - for(unsigned int mu = 0; mu < Nd - 1; mu++) - { - LatticeCoordinate(coor, mu); - ph = ph + mom[mu]*coor*((1./(_grid->_fdimensions[mu]))); - } - ph = exp((Real)(2*M_PI)*i*ph); - q_out = zero; LatticeInteger coords(_grid); LatticeCoordinate(coords, Tp); // Need q(x + mu) and q(x - mu). tmp = Cshift(q_in, mu, 1); - tmpFwd = tmp*ph; - tmp = ph*q_in; + tmpFwd = tmp*lattice_cmplx; + tmp = lattice_cmplx*q_in; tmpBwd = Cshift(tmp, mu, -1); parallel_for (unsigned int sU = 0; sU < Umu._grid->oSites(); ++sU) @@ -449,6 +509,8 @@ void WilsonFermion::SeqConservedCurrent(PropagatorField &q_in, Umu, sU, mu, t_mask); } } + + } FermOpTemplateInstantiate(WilsonFermion); diff --git a/lib/qcd/action/fermion/WilsonFermion.h b/lib/qcd/action/fermion/WilsonFermion.h index feba40ed..db32529d 100644 --- a/lib/qcd/action/fermion/WilsonFermion.h +++ b/lib/qcd/action/fermion/WilsonFermion.h @@ -155,13 +155,13 @@ class WilsonFermion : public WilsonKernels, public WilsonFermionStatic { PropagatorField &q_out, Current curr_type, unsigned int mu); - void SeqConservedCurrent(PropagatorField &q_in, - PropagatorField &q_out, - Current curr_type, - unsigned int mu, - std::vector mom, - unsigned int tmin, - unsigned int tmax); + void SeqConservedCurrent(PropagatorField &q_in, + PropagatorField &q_out, + Current curr_type, + unsigned int mu, + unsigned int tmin, + unsigned int tmax, + Lattice> &lattice_cmplx); }; typedef WilsonFermion WilsonFermionF; diff --git a/lib/qcd/action/fermion/WilsonFermion5D.cc b/lib/qcd/action/fermion/WilsonFermion5D.cc index bc703187..92280c13 100644 --- a/lib/qcd/action/fermion/WilsonFermion5D.cc +++ b/lib/qcd/action/fermion/WilsonFermion5D.cc @@ -779,18 +779,110 @@ void WilsonFermion5D::ContractConservedCurrent(PropagatorField &q_in_1, } +//template +//void WilsonFermion5D::SeqConservedCurrent(PropagatorField &q_in, +// PropagatorField &q_out, +// Current curr_type, +// unsigned int mu, +// std::vector mom, +// unsigned int tmin, +// unsigned int tmax) +//{ +// conformable(q_in._grid, FermionGrid()); +// conformable(q_in._grid, q_out._grid); +// Lattice> ph(FermionGrid()), coor(FermionGrid()); +// PropagatorField tmpFwd(FermionGrid()), tmpBwd(FermionGrid()), +// tmp(FermionGrid()); +// Complex i(0.0, 1.0); +// unsigned int tshift = (mu == Tp) ? 1 : 0; +// unsigned int LLs = q_in._grid->_rdimensions[0]; +// unsigned int LLt = GridDefaultLatt()[Tp]; +// +// // Momentum projection. +// ph = zero; +// for(unsigned int nu = 0; nu < Nd - 1; nu++) +// { +// // Shift coordinate lattice index by 1 to account for 5th dimension. +// LatticeCoordinate(coor, nu + 1); +// ph = ph + mom[nu]*coor*((1./(_FourDimGrid->_fdimensions[nu]))); +// } +// ph = exp((Real)(2*M_PI)*i*ph); +// +// q_out = zero; +// LatticeInteger coords(_FourDimGrid); +// LatticeCoordinate(coords, Tp); +// +// +// // Need q(x + mu, s) and q(x - mu, s). 5D lattice so shift 4D coordinate mu +// // by one. +// tmp = Cshift(q_in, mu + 1, 1); +// tmpFwd = tmp*ph; +// tmp = ph*q_in; +// tmpBwd = Cshift(tmp, mu + 1, -1); +// +// parallel_for (unsigned int sU = 0; sU < Umu._grid->oSites(); ++sU) +// { +// // Compute the sequential conserved current insertion only if our simd +// // object contains a timeslice we need. +// vInteger t_mask = ((coords._odata[sU] >= tmin) && +// (coords._odata[sU] <= tmax)); +// Integer timeSlices = Reduce(t_mask); +// +// if (timeSlices > 0) +// { +// unsigned int sF = sU * LLs; +// for (unsigned int s = 0; s < LLs; ++s) +// { +// bool axial_sign = ((curr_type == Current::Axial) && (s < (LLs / 2))); +// bool tadpole_sign = (curr_type == Current::Tadpole); +// bool switch_sgn = tadpole_sign || axial_sign; +// +// Kernels::SeqConservedCurrentSiteFwd(tmpFwd._odata[sF], +// q_out._odata[sF], Umu, sU, +// mu, t_mask, switch_sgn); +// ++sF; +// } +// } +// +// // Repeat for backward direction. +// t_mask = ((coords._odata[sU] >= (tmin + tshift)) && +// (coords._odata[sU] <= (tmax + tshift))); +// +// //if tmax = LLt-1 (last timeslice) include timeslice 0 if the time is shifted (mu=3) +// unsigned int t0 = 0; +// if((tmax==LLt-1) && (tshift==1)) t_mask = (t_mask || (coords._odata[sU] == t0 )); +// +// timeSlices = Reduce(t_mask); +// +// if (timeSlices > 0) +// { +// unsigned int sF = sU * LLs; +// for (unsigned int s = 0; s < LLs; ++s) +// { +// bool axial_sign = ((curr_type == Current::Axial) && (s < (LLs / 2))); +// Kernels::SeqConservedCurrentSiteBwd(tmpBwd._odata[sF], +// q_out._odata[sF], Umu, sU, +// mu, t_mask, axial_sign); +// ++sF; +// } +// } +// } +//} + + + + template void WilsonFermion5D::SeqConservedCurrent(PropagatorField &q_in, PropagatorField &q_out, Current curr_type, unsigned int mu, - std::vector mom, unsigned int tmin, - unsigned int tmax) + unsigned int tmax, + Lattice> &lattice_cmplx) { conformable(q_in._grid, FermionGrid()); conformable(q_in._grid, q_out._grid); - Lattice> ph(FermionGrid()), coor(FermionGrid()); PropagatorField tmpFwd(FermionGrid()), tmpBwd(FermionGrid()), tmp(FermionGrid()); Complex i(0.0, 1.0); @@ -798,25 +890,26 @@ void WilsonFermion5D::SeqConservedCurrent(PropagatorField &q_in, unsigned int LLs = q_in._grid->_rdimensions[0]; unsigned int LLt = GridDefaultLatt()[Tp]; - // Momentum projection. - ph = zero; - for(unsigned int nu = 0; nu < Nd - 1; nu++) - { - // Shift coordinate lattice index by 1 to account for 5th dimension. - LatticeCoordinate(coor, nu + 1); - ph = ph + mom[nu]*coor*((1./(_FourDimGrid->_fdimensions[nu]))); - } - ph = exp((Real)(2*M_PI)*i*ph); - q_out = zero; LatticeInteger coords(_FourDimGrid); LatticeCoordinate(coords, Tp); + + //QED: photon field is 4dim, but need a 5dim object to multiply to + // DWF PropagatorField + Lattice> lattice_cmplx_5d(FermionGrid()); + for (unsigned int s = 0; s < LLs; ++s) + { + InsertSlice(lattice_cmplx,lattice_cmplx_5d, s, 0); + } + + + // Need q(x + mu, s) and q(x - mu, s). 5D lattice so shift 4D coordinate mu // by one. tmp = Cshift(q_in, mu + 1, 1); - tmpFwd = tmp*ph; - tmp = ph*q_in; + tmpFwd = tmp*lattice_cmplx_5d; + tmp = lattice_cmplx_5d*q_in; tmpBwd = Cshift(tmp, mu + 1, -1); parallel_for (unsigned int sU = 0; sU < Umu._grid->oSites(); ++sU) @@ -868,6 +961,10 @@ void WilsonFermion5D::SeqConservedCurrent(PropagatorField &q_in, } } + + + + FermOpTemplateInstantiate(WilsonFermion5D); GparityFermOpTemplateInstantiate(WilsonFermion5D); diff --git a/lib/qcd/action/fermion/WilsonFermion5D.h b/lib/qcd/action/fermion/WilsonFermion5D.h index d66f4a1d..6a7c6c7d 100644 --- a/lib/qcd/action/fermion/WilsonFermion5D.h +++ b/lib/qcd/action/fermion/WilsonFermion5D.h @@ -222,13 +222,13 @@ namespace QCD { PropagatorField &q_out, Current curr_type, unsigned int mu); - void SeqConservedCurrent(PropagatorField &q_in, + void SeqConservedCurrent(PropagatorField &q_in, PropagatorField &q_out, - Current curr_type, + Current curr_type, unsigned int mu, - std::vector mom, - unsigned int tmin, - unsigned int tmax); + unsigned int tmin, + unsigned int tmax, + Lattice> &lattice_cmplx); }; }} diff --git a/tests/hadrons/Test_hadrons.hpp b/tests/hadrons/Test_hadrons.hpp index 0265f5a6..6188a580 100644 --- a/tests/hadrons/Test_hadrons.hpp +++ b/tests/hadrons/Test_hadrons.hpp @@ -263,7 +263,8 @@ inline void makeConservedSequentialSource(Application &application, seqPar.tA = tS; seqPar.tB = tS; seqPar.curr_type = curr; - seqPar.mu = mu; + seqPar.mu_min = mu; + seqPar.mu_min = mu; seqPar.mom = mom; application.createModule(srcName, seqPar); } From c6e1f645737d6785d513246db47323afaa5850e2 Mon Sep 17 00:00:00 2001 From: Vera Guelpers Date: Tue, 13 Feb 2018 09:30:23 +0000 Subject: [PATCH 06/38] Test for QED --- tests/hadrons/Test_QED.cc | 260 ++++++++++++++++++++++++++++++++++++++ 1 file changed, 260 insertions(+) create mode 100644 tests/hadrons/Test_QED.cc diff --git a/tests/hadrons/Test_QED.cc b/tests/hadrons/Test_QED.cc new file mode 100644 index 00000000..44dd0578 --- /dev/null +++ b/tests/hadrons/Test_QED.cc @@ -0,0 +1,260 @@ +/******************************************************************************* + Grid physics library, www.github.com/paboyle/Grid + + Source file: tests/hadrons/Test_hadrons_spectrum.cc + + Copyright (C) 2015-2018 + + Author: Antonin Portelli + Author: Vera Guelpers + + This program is free software; you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation; either version 2 of the License, or + (at your option) any later version. + + This program is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License along + with this program; if not, write to the Free Software Foundation, Inc., + 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. + + See the full license in the file "LICENSE" in the top level distribution + directory. + *******************************************************************************/ + +#include + +using namespace Grid; +using namespace Hadrons; + +int main(int argc, char *argv[]) +{ + // initialization ////////////////////////////////////////////////////////// + Grid_init(&argc, &argv); + HadronsLogError.Active(GridLogError.isActive()); + HadronsLogWarning.Active(GridLogWarning.isActive()); + HadronsLogMessage.Active(GridLogMessage.isActive()); + HadronsLogIterative.Active(GridLogIterative.isActive()); + HadronsLogDebug.Active(GridLogDebug.isActive()); + LOG(Message) << "Grid initialized" << std::endl; + + + // run setup /////////////////////////////////////////////////////////////// + Application application; + std::vector flavour = {"h"}; //{"l", "s", "c1", "c2", "c3"}; + std::vector mass = {.2}; //{.01, .04, .2 , .25 , .3 }; + + unsigned int nt = GridDefaultLatt()[Tp]; + + // global parameters + Application::GlobalPar globalPar; + globalPar.trajCounter.start = 1500; + globalPar.trajCounter.end = 1520; + globalPar.trajCounter.step = 20; + globalPar.seed = "1 2 3 4"; + application.setPar(globalPar); + // gauge field + application.createModule("gauge"); + // pt source + 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"; + + //stochastic photon field + MGauge::StochEm::Par photonPar; + photonPar.gauge = PhotonR::Gauge::feynman; + photonPar.zmScheme = PhotonR::ZmScheme::qedL; + application.createModule("ph_field", photonPar); + + + + for (unsigned int i = 0; i < flavour.size(); ++i) + { + // actions + MAction::DWF::Par actionPar; + actionPar.gauge = "gauge"; + actionPar.Ls = 8; + actionPar.M5 = 1.8; + actionPar.mass = mass[i]; + actionPar.boundary = boundary; + application.createModule("DWF_" + flavour[i], actionPar); + + + // solvers + MSolver::RBPrecCG::Par solverPar; + solverPar.action = "DWF_" + flavour[i]; + solverPar.residual = 1.0e-8; + application.createModule("CG_" + flavour[i], + solverPar); + + // propagators + MFermion::GaugeProp::Par quarkPar; + quarkPar.solver = "CG_" + flavour[i]; + quarkPar.source = "pt"; + application.createModule("Qpt_" + flavour[i], + quarkPar); + + + //seq sources with tadpole insertion + MSource::SeqConserved::Par seqPar_T; + seqPar_T.q = "Qpt_" + flavour[i] + "_5d"; + seqPar_T.action = "DWF_" + flavour[i]; + seqPar_T.tA = 0; + seqPar_T.tB = nt-1; + seqPar_T.curr_type = Current::Tadpole; + seqPar_T.mu_min = 0; + seqPar_T.mu_max = 3; + seqPar_T.mom = "0. 0. 0. 0."; + application.createModule("Qpt_" + flavour[i] + + "_seq_T", seqPar_T); + // seq propagator with tadpole insertion + MFermion::GaugeProp::Par quarkPar_seq_T; + quarkPar_seq_T.solver = "CG_" + flavour[i]; + quarkPar_seq_T.source = "Qpt_" + flavour[i] + "_seq_T"; + application.createModule("Qpt_" + flavour[i] + + "_seq_T" + flavour[i], + quarkPar_seq_T); + + + + //seq sources with conserved vector and photon insertion + MSource::SeqConserved::Par seqPar_V; + seqPar_V.q = "Qpt_" + flavour[i] + "_5d"; + seqPar_V.action = "DWF_" + flavour[i]; + seqPar_V.tA = 0; + seqPar_V.tB = nt-1; + seqPar_V.curr_type = Current::Vector; + seqPar_V.mu_min = 0; + seqPar_V.mu_max = 3; + seqPar_V.mom = "0. 0. 0. 0."; + seqPar_V.photon = "ph_field"; + application.createModule("Qpt_" + flavour[i] + + "_seq_V_ph", seqPar_V); + // seq propagator with conserved vector and photon insertion + MFermion::GaugeProp::Par quarkPar_seq_V; + quarkPar_seq_V.solver = "CG_" + flavour[i]; + quarkPar_seq_V.source = "Qpt_" + flavour[i] + "_seq_V_ph"; + application.createModule("Qpt_" + flavour[i] + + "_seq_V_ph_" + flavour[i], + quarkPar_seq_V); + + + + //double seq sources with conserved vector and photon insertion + //(for self energy) + MSource::SeqConserved::Par seqPar_VV; + seqPar_VV.q = "Qpt_" + flavour[i] + "_seq_V_ph_" + + flavour[i] + "_5d"; + seqPar_VV.action = "DWF_" + flavour[i]; + seqPar_VV.tA = 0; + seqPar_VV.tB = nt-1; + seqPar_VV.curr_type = Current::Vector; + seqPar_VV.mu_min = 0; + seqPar_VV.mu_max = 3; + seqPar_VV.mom = "0. 0. 0. 0."; + seqPar_VV.photon = "ph_field"; + application.createModule("Qpt_" + flavour[i] + + "_seq_V_ph" + flavour[i] + + "_seq_V_ph", seqPar_VV); + //double seq propagator with conserved vector and photon insertion + MFermion::GaugeProp::Par quarkPar_seq_VV; + quarkPar_seq_VV.solver = "CG_" + flavour[i]; + quarkPar_seq_VV.source = "Qpt_" + flavour[i] + "_seq_V_ph" + + flavour[i] + "_seq_V_ph"; + application.createModule("Qpt_" + flavour[i] + + "_seq_V_ph_" + flavour[i] + + "_seq_V_ph_" + flavour[i], + quarkPar_seq_VV); + + + + } + for (unsigned int i = 0; i < flavour.size(); ++i) + for (unsigned int j = i; j < flavour.size(); ++j) + { + //2pt function contraction + MContraction::Meson::Par mesPar; + mesPar.output = "QED/pt_" + flavour[i] + flavour[j]; + mesPar.q1 = "Qpt_" + flavour[i]; + mesPar.q2 = "Qpt_" + flavour[j]; + mesPar.gammas = ""; + mesPar.sink = "sink"; + application.createModule("meson_pt_" + + flavour[i] + flavour[j], + mesPar); + + + + //tadpole contraction + MContraction::Meson::Par mesPar_seq_T; + mesPar_seq_T.output = "QED/tadpole_pt_" + flavour[i] + "_T_" + + flavour[i] + "__" + flavour[j]; + mesPar_seq_T.q1 = "Qpt_" + flavour[i] + "_seq_T" + flavour[i]; + mesPar_seq_T.q2 = "Qpt_" + flavour[j]; + mesPar_seq_T.gammas = ""; + mesPar_seq_T.sink = "sink"; + application.createModule("meson_tadpole_pt_" + + flavour[i] + "_seq_T" + + flavour[i] + flavour[j], + mesPar_seq_T); + + + + //photon exchange contraction + MContraction::Meson::Par mesPar_seq_E; + mesPar_seq_E.output = "QED/exchange_pt_" + flavour[i] + "_V_ph_" + + flavour[i] + "__" + flavour[j] + "_V_ph_" + + flavour[j]; + mesPar_seq_E.q1 = "Qpt_" + flavour[i] + "_seq_V_ph_" + flavour[i]; + mesPar_seq_E.q2 = "Qpt_" + flavour[j] + "_seq_V_ph_" + flavour[j]; + mesPar_seq_E.gammas = ""; + mesPar_seq_E.sink = "sink"; + application.createModule("meson_exchange_pt_" + + flavour[i] + "_seq_V_ph_" + flavour[i] + + flavour[j] + "_seq_V_ph_" + flavour[j], + mesPar_seq_E); + + + + //self energy contraction + MContraction::Meson::Par mesPar_seq_S; + mesPar_seq_S.output = "QED/selfenergy_pt_" + flavour[i] + "_V_ph_" + + flavour[i] + "_V_ph_" + flavour[i] + "__" + + flavour[j]; + mesPar_seq_S.q1 = "Qpt_" + flavour[i] + "_seq_V_ph_" + flavour[i] + + "_seq_V_ph_" + flavour[i]; + mesPar_seq_S.q2 = "Qpt_" + flavour[j]; + mesPar_seq_S.gammas = ""; + mesPar_seq_S.sink = "sink"; + application.createModule("meson_selfenergy_pt_" + + flavour[i] + "_seq_V_ph_" + + flavour[i] + "_seq_V_ph_" + + flavour[i] + flavour[j], + mesPar_seq_S); + + } + + + + + // execution + application.saveParameterFile("QED.xml"); + application.run(); + + // epilogue + LOG(Message) << "Grid is finalizing now" << std::endl; + Grid_finalize(); + + return EXIT_SUCCESS; +} From d55212c998bba55b78bb676eb8a33f86fe6825cc Mon Sep 17 00:00:00 2001 From: Vera Guelpers Date: Wed, 14 Feb 2018 10:45:18 +0000 Subject: [PATCH 07/38] restructure SeqConservedCurrent for DWF to need less memory --- lib/qcd/action/fermion/WilsonFermion5D.cc | 96 ++++++++++------------- 1 file changed, 43 insertions(+), 53 deletions(-) diff --git a/lib/qcd/action/fermion/WilsonFermion5D.cc b/lib/qcd/action/fermion/WilsonFermion5D.cc index 92280c13..6b830a03 100644 --- a/lib/qcd/action/fermion/WilsonFermion5D.cc +++ b/lib/qcd/action/fermion/WilsonFermion5D.cc @@ -883,8 +883,7 @@ void WilsonFermion5D::SeqConservedCurrent(PropagatorField &q_in, { conformable(q_in._grid, FermionGrid()); conformable(q_in._grid, q_out._grid); - PropagatorField tmpFwd(FermionGrid()), tmpBwd(FermionGrid()), - tmp(FermionGrid()); + PropagatorField tmp(GaugeGrid()),tmp2(GaugeGrid()); Complex i(0.0, 1.0); unsigned int tshift = (mu == Tp) ? 1 : 0; unsigned int LLs = q_in._grid->_rdimensions[0]; @@ -895,69 +894,60 @@ void WilsonFermion5D::SeqConservedCurrent(PropagatorField &q_in, LatticeCoordinate(coords, Tp); - //QED: photon field is 4dim, but need a 5dim object to multiply to - // DWF PropagatorField - Lattice> lattice_cmplx_5d(FermionGrid()); for (unsigned int s = 0; s < LLs; ++s) - { - InsertSlice(lattice_cmplx,lattice_cmplx_5d, s, 0); - } - - - - // Need q(x + mu, s) and q(x - mu, s). 5D lattice so shift 4D coordinate mu - // by one. - tmp = Cshift(q_in, mu + 1, 1); - tmpFwd = tmp*lattice_cmplx_5d; - tmp = lattice_cmplx_5d*q_in; - tmpBwd = Cshift(tmp, mu + 1, -1); - - parallel_for (unsigned int sU = 0; sU < Umu._grid->oSites(); ++sU) { - // Compute the sequential conserved current insertion only if our simd - // object contains a timeslice we need. - vInteger t_mask = ((coords._odata[sU] >= tmin) && - (coords._odata[sU] <= tmax)); - Integer timeSlices = Reduce(t_mask); + bool axial_sign = ((curr_type == Current::Axial) && (s < (LLs / 2))); + bool tadpole_sign = (curr_type == Current::Tadpole); + bool switch_sgn = tadpole_sign || axial_sign; - if (timeSlices > 0) - { - unsigned int sF = sU * LLs; - for (unsigned int s = 0; s < LLs; ++s) + + //forward direction: Need q(x + mu, s)*A(x) + ExtractSlice(tmp2, q_in, s, 0); //q(x,s) + tmp = Cshift(tmp2, mu, 1); //q(x+mu,s) + tmp2 = tmp*lattice_cmplx; //q(x+mu,s)*A(x) + + parallel_for (unsigned int sU = 0; sU < Umu._grid->oSites(); ++sU) + { + // Compute the sequential conserved current insertion only if our simd + // object contains a timeslice we need. + vInteger t_mask = ((coords._odata[sU] >= tmin) && + (coords._odata[sU] <= tmax)); + Integer timeSlices = Reduce(t_mask); + + if (timeSlices > 0) { - bool axial_sign = ((curr_type == Current::Axial) && (s < (LLs / 2))); - bool tadpole_sign = (curr_type == Current::Tadpole); - bool switch_sgn = tadpole_sign || axial_sign; - - Kernels::SeqConservedCurrentSiteFwd(tmpFwd._odata[sF], - q_out._odata[sF], Umu, sU, - mu, t_mask, switch_sgn); - ++sF; + unsigned int sF = sU * LLs + s; + Kernels::SeqConservedCurrentSiteFwd(tmp2._odata[sU], + q_out._odata[sF], Umu, sU, + mu, t_mask, switch_sgn); } + } - // Repeat for backward direction. - t_mask = ((coords._odata[sU] >= (tmin + tshift)) && - (coords._odata[sU] <= (tmax + tshift))); + //backward direction: Need q(x - mu, s)*A(x-mu) + ExtractSlice(tmp2, q_in, s, 0); //q(x,s) + tmp = lattice_cmplx*tmp2; //q(x,s)*A(x) + tmp2 = Cshift(tmp, mu, -1); //q(x-mu,s)*A(x-mu,s) - //if tmax = LLt-1 (last timeslice) include timeslice 0 if the time is shifted (mu=3) - unsigned int t0 = 0; - if((tmax==LLt-1) && (tshift==1)) t_mask = (t_mask || (coords._odata[sU] == t0 )); + parallel_for (unsigned int sU = 0; sU < Umu._grid->oSites(); ++sU) + { + vInteger t_mask = ((coords._odata[sU] >= (tmin + tshift)) && + (coords._odata[sU] <= (tmax + tshift))); - timeSlices = Reduce(t_mask); + //if tmax = LLt-1 (last timeslice) include timeslice 0 if the time is shifted (mu=3) + unsigned int t0 = 0; + if((tmax==LLt-1) && (tshift==1)) t_mask = (t_mask || (coords._odata[sU] == t0 )); - if (timeSlices > 0) - { - unsigned int sF = sU * LLs; - for (unsigned int s = 0; s < LLs; ++s) + Integer timeSlices = Reduce(t_mask); + + if (timeSlices > 0) { - bool axial_sign = ((curr_type == Current::Axial) && (s < (LLs / 2))); - Kernels::SeqConservedCurrentSiteBwd(tmpBwd._odata[sF], - q_out._odata[sF], Umu, sU, - mu, t_mask, axial_sign); - ++sF; + unsigned int sF = sU * LLs + s; + Kernels::SeqConservedCurrentSiteBwd(tmp2._odata[sU], + q_out._odata[sF], Umu, sU, + mu, t_mask, axial_sign); } - } + } } } From fdbd42e542dee144c0fd238917baaee90399ffac Mon Sep 17 00:00:00 2001 From: Antonin Portelli Date: Mon, 26 Feb 2018 19:22:43 +0000 Subject: [PATCH 08/38] Hadrons: first implementation of local coherence Lanczos --- extras/Hadrons/Environment.cc | 2 +- extras/Hadrons/Modules.hpp | 31 +-- .../Modules/MSolver/LocalCoherenceLanczos.hpp | 211 ++++++++++++++++++ 3 files changed, 214 insertions(+), 30 deletions(-) create mode 100644 extras/Hadrons/Modules/MSolver/LocalCoherenceLanczos.hpp diff --git a/extras/Hadrons/Environment.cc b/extras/Hadrons/Environment.cc index 82b0dda1..6554122e 100644 --- a/extras/Hadrons/Environment.cc +++ b/extras/Hadrons/Environment.cc @@ -270,7 +270,7 @@ int Environment::getObjectModule(const std::string name) const unsigned int Environment::getObjectLs(const unsigned int address) const { - if (hasObject(address)) + if (hasCreatedObject(address)) { return object_[address].Ls; } diff --git a/extras/Hadrons/Modules.hpp b/extras/Hadrons/Modules.hpp index 135934fb..d1f947a2 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-2018 - -Author: Antonin Portelli -Author: Lanny91 - -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 @@ -43,11 +15,12 @@ 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/MSolver/LocalCoherenceLanczos.hpp b/extras/Hadrons/Modules/MSolver/LocalCoherenceLanczos.hpp new file mode 100644 index 00000000..1ce87c2b --- /dev/null +++ b/extras/Hadrons/Modules/MSolver/LocalCoherenceLanczos.hpp @@ -0,0 +1,211 @@ +#ifndef Hadrons_MSolver_LocalCoherenceLanczos_hpp_ +#define Hadrons_MSolver_LocalCoherenceLanczos_hpp_ + +#include +#include +#include +#include + +#ifndef DEFAULT_LANCZOS_NBASIS +#define DEFAULT_LANCZOS_NBASIS 60 +#endif + +BEGIN_HADRONS_NAMESPACE + +/****************************************************************************** + * LocalCoherenceLanczos * + ******************************************************************************/ +BEGIN_MODULE_NAMESPACE(MSolver) + +class LocalCoherenceLanczosPar: Serializable +{ +public: + GRID_SERIALIZABLE_CLASS_MEMBERS(LocalCoherenceLanczosPar, + std::string, action, + int, doFine, + int, doCoarse, + LanczosParams, fineParams, + LanczosParams, coarseParams, + ChebyParams, smoother, + RealD, coarseRelaxTol, + std::string, blockSize); +}; + +template +class TLocalCoherenceLanczos: public Module +{ +public: + FERM_TYPE_ALIASES(FImpl,); + typedef LocalCoherenceLanczos LCL; + typedef typename LCL::FineField FineField; + typedef typename LCL::CoarseField CoarseField; + typedef SchurDiagMooeeOperator SchurFMat; +public: + // constructor + TLocalCoherenceLanczos(const std::string name); + // destructor + virtual ~TLocalCoherenceLanczos(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 makeCoarseGrid(void); +private: + std::vector coarseDim_; + int Ls_, cLs_{1}; + std::unique_ptr coarseGrid4_{nullptr}; + std::unique_ptr coarseGrid_{nullptr}; + std::unique_ptr coarseGrid4Rb_{nullptr}; + std::unique_ptr coarseGridRb_{nullptr}; + std::string fevecName_, cevecName_, fevalName_, cevalName_; +}; + +MODULE_REGISTER_NS(LocalCoherenceLanczos, + ARG(TLocalCoherenceLanczos), + MSolver); + +/****************************************************************************** + * TLocalCoherenceLanczos implementation * + ******************************************************************************/ +// constructor ///////////////////////////////////////////////////////////////// +template +TLocalCoherenceLanczos::TLocalCoherenceLanczos(const std::string name) +: Module(name) +{ + fevecName_ = getName() + "_fineEvec"; + cevecName_ = getName() + "_coarseEvec"; + fevalName_ = getName() + "_fineEval"; + cevalName_ = getName() + "_coarseEval"; +} + +// dependencies/products /////////////////////////////////////////////////////// +template +std::vector TLocalCoherenceLanczos::getInput(void) +{ + std::vector in = {par().action}; + + return in; +} + +template +std::vector TLocalCoherenceLanczos::getOutput(void) +{ + std::vector out = {fevecName_, cevecName_, fevalName_, cevalName_}; + + return out; +} + +// setup /////////////////////////////////////////////////////////////////////// +template +void TLocalCoherenceLanczos::makeCoarseGrid(void) +{ + int nd = env().getNd(); + std::vector blockSize = strToVec(par().blockSize); + auto fineDim = env().getDim(); + + Ls_ = env().getObjectLs(par().action); + env().createGrid(Ls_); + coarseDim_.resize(nd); + for (int d = 0; d < coarseDim_.size(); d++) + { + coarseDim_[d] = fineDim[d]/blockSize[d]; + if (coarseDim_[d]*blockSize[d] != fineDim[d]) + { + HADRON_ERROR(Size, "Fine dimension " + std::to_string(d) + + " (" + std::to_string(fineDim[d]) + + ") not divisible by coarse dimension (" + + std::to_string(coarseDim_[d]) + ")"); + } + } + if (blockSize.size() > nd) + { + cLs_ = Ls_/blockSize[nd]; + if (cLs_*blockSize[nd] != Ls_) + { + HADRON_ERROR(Size, "Fine Ls (" + std::to_string(Ls_) + + ") not divisible by coarse Ls (" + + std::to_string(cLs_) + ")"); + } + } + if (Ls_ > 1) + { + coarseGrid4_.reset(SpaceTimeGrid::makeFourDimGrid( + coarseDim_, GridDefaultSimd(nd, vComplex::Nsimd()), + GridDefaultMpi())); + coarseGrid4Rb_.reset(SpaceTimeGrid::makeFourDimRedBlackGrid(coarseGrid4_.get())); + coarseGrid_.reset(SpaceTimeGrid::makeFiveDimGrid(cLs_, coarseGrid4_.get())); + coarseGridRb_.reset(SpaceTimeGrid::makeFiveDimRedBlackGrid(cLs_, coarseGrid4_.get())); + } + else + { + coarseGrid_.reset(SpaceTimeGrid::makeFourDimGrid( + coarseDim_, GridDefaultSimd(nd, vComplex::Nsimd()), + GridDefaultMpi())); + coarseGridRb_.reset(SpaceTimeGrid::makeFourDimRedBlackGrid(coarseGrid_.get())); + } +} + +template +void TLocalCoherenceLanczos::setup(void) +{ + LOG(Message) << "Setting up local coherence Lanczos eigensolver for" + << " action '" << par().action << "' (" << nBasis + << " eigenvectors)..." << std::endl; + + if (!coarseGrid_) + { + makeCoarseGrid(); + } + envCreate(std::vector, fevecName_, Ls_, par().fineParams.Nm, + env().getRbGrid(Ls_)); + envCreate(std::vector, cevecName_, Ls_, par().coarseParams.Nm, + coarseGridRb_.get()); + envCreate(std::vector, fevalName_, Ls_, par().fineParams.Nm); + envCreate(std::vector, cevalName_, Ls_, par().coarseParams.Nm); + envTmp(SchurFMat, "mat", Ls_, envGet(FMat, par().action)); + envGetTmp(SchurFMat, mat); + envTmp(LCL, "solver", Ls_, env().getRbGrid(Ls_), coarseGridRb_.get(), mat, Odd, + envGet(std::vector, fevecName_), + envGet(std::vector, cevecName_), + envGet(std::vector, fevalName_), + envGet(std::vector, cevalName_)); +} + +// execution /////////////////////////////////////////////////////////////////// +template +void TLocalCoherenceLanczos::execute(void) +{ + auto &fine = par().fineParams; + auto &coarse = par().coarseParams; + + envGetTmp(LCL, solver); + if (par().doFine) + { + LOG(Message) << "Performing fine grid IRL -- Nstop= " + << fine.Nstop << ", Nk= " << fine.Nk << ", Nm= " + << fine.Nm << std::endl; + solver.calcFine(fine.Cheby, fine.Nstop, fine.Nk, fine.Nm, fine.resid, + fine.MaxIt, fine.betastp, fine.MinRes); + LOG(Message) << "Orthogonalising" << std::endl; + solver.Orthogonalise(); + } + if (par().doCoarse) + { + LOG(Message) << "Performing coarse grid IRL -- Nstop= " + << fine.Nstop << ", Nk= " << fine.Nk << ", Nm= " + << fine.Nm << std::endl; + solver.calcCoarse(coarse.Cheby, par().smoother, par().coarseRelaxTol, + coarse.Nstop, coarse.Nk, coarse.Nm, coarse.resid, + coarse.MaxIt, coarse.betastp, coarse.MinRes); + } +} + +END_MODULE_NAMESPACE + +END_HADRONS_NAMESPACE + +#endif // Hadrons_MSolver_LocalCoherenceLanczos_hpp_ From 43578a3eb496fe97cfc791eb6a44584c0e1f738a Mon Sep 17 00:00:00 2001 From: Antonin Portelli Date: Mon, 26 Feb 2018 19:24:19 +0000 Subject: [PATCH 09/38] Hadrons: copyright update --- extras/Hadrons/Modules.hpp | 30 +++++++++++++++++++ .../Hadrons/Modules/MAction/WilsonClover.hpp | 8 ++--- extras/Hadrons/Modules/MFermion/GaugeProp.hpp | 2 ++ extras/Hadrons/Modules/MGauge/FundtoHirep.hpp | 6 ++-- .../Modules/MSolver/LocalCoherenceLanczos.hpp | 27 +++++++++++++++++ extras/Hadrons/modules.inc | 5 ++-- 6 files changed, 67 insertions(+), 11 deletions(-) diff --git a/extras/Hadrons/Modules.hpp b/extras/Hadrons/Modules.hpp index d1f947a2..1334a9b4 100644 --- a/extras/Hadrons/Modules.hpp +++ b/extras/Hadrons/Modules.hpp @@ -1,3 +1,33 @@ +/************************************************************************************* + +Grid physics library, www.github.com/paboyle/Grid + +Source file: extras/Hadrons/Modules.hpp + +Copyright (C) 2015-2018 + +Author: Antonin Portelli +Author: Guido Cossu +Author: Lanny91 +Author: pretidav + +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 diff --git a/extras/Hadrons/Modules/MAction/WilsonClover.hpp b/extras/Hadrons/Modules/MAction/WilsonClover.hpp index c369f086..d5174c6d 100644 --- a/extras/Hadrons/Modules/MAction/WilsonClover.hpp +++ b/extras/Hadrons/Modules/MAction/WilsonClover.hpp @@ -2,12 +2,12 @@ Grid physics library, www.github.com/paboyle/Grid -Source file: extras/Hadrons/Modules/MAction/Wilson.hpp +Source file: extras/Hadrons/Modules/MAction/WilsonClover.hpp -Copyright (C) 2015 -Copyright (C) 2016 +Copyright (C) 2015-2018 -Author: Antonin Portelli +Author: Guido Cossu +Author: pretidav 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 diff --git a/extras/Hadrons/Modules/MFermion/GaugeProp.hpp b/extras/Hadrons/Modules/MFermion/GaugeProp.hpp index 412e92d5..45b9de6c 100644 --- a/extras/Hadrons/Modules/MFermion/GaugeProp.hpp +++ b/extras/Hadrons/Modules/MFermion/GaugeProp.hpp @@ -7,7 +7,9 @@ Source file: extras/Hadrons/Modules/MFermion/GaugeProp.hpp Copyright (C) 2015-2018 Author: Antonin Portelli +Author: Guido Cossu Author: Lanny91 +Author: pretidav 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 diff --git a/extras/Hadrons/Modules/MGauge/FundtoHirep.hpp b/extras/Hadrons/Modules/MGauge/FundtoHirep.hpp index 6f072783..4fba59f0 100644 --- a/extras/Hadrons/Modules/MGauge/FundtoHirep.hpp +++ b/extras/Hadrons/Modules/MGauge/FundtoHirep.hpp @@ -4,11 +4,9 @@ Grid physics library, www.github.com/paboyle/Grid Source file: extras/Hadrons/Modules/MGauge/FundtoHirep.hpp -Copyright (C) 2015 -Copyright (C) 2016 +Copyright (C) 2015-2018 -Author: David Preti - Guido Cossu +Author: pretidav 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 diff --git a/extras/Hadrons/Modules/MSolver/LocalCoherenceLanczos.hpp b/extras/Hadrons/Modules/MSolver/LocalCoherenceLanczos.hpp index 1ce87c2b..5fd486fd 100644 --- a/extras/Hadrons/Modules/MSolver/LocalCoherenceLanczos.hpp +++ b/extras/Hadrons/Modules/MSolver/LocalCoherenceLanczos.hpp @@ -1,3 +1,30 @@ +/************************************************************************************* + +Grid physics library, www.github.com/paboyle/Grid + +Source file: extras/Hadrons/Modules/MSolver/LocalCoherenceLanczos.hpp + +Copyright (C) 2015-2018 + +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 */ #ifndef Hadrons_MSolver_LocalCoherenceLanczos_hpp_ #define Hadrons_MSolver_LocalCoherenceLanczos_hpp_ diff --git a/extras/Hadrons/modules.inc b/extras/Hadrons/modules.inc index 90602275..4a9b2630 100644 --- a/extras/Hadrons/modules.inc +++ b/extras/Hadrons/modules.inc @@ -1,6 +1,4 @@ modules_cc =\ - Modules/MScalar/ChargedProp.cc \ - Modules/MScalar/FreeProp.cc \ Modules/MContraction/WeakHamiltonianEye.cc \ Modules/MContraction/WeakNeutral4ptDisc.cc \ Modules/MContraction/WeakHamiltonianNonEye.cc \ @@ -30,11 +28,12 @@ modules_hpp =\ Modules/MSource/SeqConserved.hpp \ Modules/MSink/Smear.hpp \ Modules/MSink/Point.hpp \ + Modules/MSolver/LocalCoherenceLanczos.hpp \ Modules/MSolver/RBPrecCG.hpp \ Modules/MGauge/Unit.hpp \ Modules/MGauge/Random.hpp \ - Modules/MGauge/StochEm.hpp \ Modules/MGauge/FundtoHirep.hpp \ + Modules/MGauge/StochEm.hpp \ Modules/MUtilities/TestSeqGamma.hpp \ Modules/MUtilities/TestSeqConserved.hpp \ Modules/MLoop/NoiseLoop.hpp \ From 8a049f27b840cad2e54f822a9c3c28cf423d996a Mon Sep 17 00:00:00 2001 From: Antonin Portelli Date: Tue, 27 Feb 2018 13:46:59 +0000 Subject: [PATCH 10/38] Hadrons: Lanczos code improvement --- extras/Hadrons/Global.cc | 3 + extras/Hadrons/Global.hpp | 1 + extras/Hadrons/LanczosUtils.hpp | 88 +++++++++++++++++++ extras/Hadrons/Makefile.am | 1 + .../Modules/MSolver/LocalCoherenceLanczos.hpp | 80 +++++++++-------- 5 files changed, 135 insertions(+), 38 deletions(-) create mode 100644 extras/Hadrons/LanczosUtils.hpp diff --git a/extras/Hadrons/Global.cc b/extras/Hadrons/Global.cc index 9a90a08c..b121434f 100644 --- a/extras/Hadrons/Global.cc +++ b/extras/Hadrons/Global.cc @@ -37,6 +37,7 @@ HadronsLogger Hadrons::HadronsLogWarning(1,"Warning"); HadronsLogger Hadrons::HadronsLogMessage(1,"Message"); HadronsLogger Hadrons::HadronsLogIterative(1,"Iterative"); HadronsLogger Hadrons::HadronsLogDebug(1,"Debug"); +HadronsLogger Hadrons::HadronsLogIRL(1,"IRL"); void Hadrons::initLogger(void) { @@ -46,11 +47,13 @@ void Hadrons::initLogger(void) GridLogMessage.setTopWidth(w); GridLogIterative.setTopWidth(w); GridLogDebug.setTopWidth(w); + GridLogIRL.setTopWidth(w); HadronsLogError.Active(GridLogError.isActive()); HadronsLogWarning.Active(GridLogWarning.isActive()); HadronsLogMessage.Active(GridLogMessage.isActive()); HadronsLogIterative.Active(GridLogIterative.isActive()); HadronsLogDebug.Active(GridLogDebug.isActive()); + HadronsLogIRL.Active(GridLogIRL.isActive()); } // type utilities ////////////////////////////////////////////////////////////// diff --git a/extras/Hadrons/Global.hpp b/extras/Hadrons/Global.hpp index fc069ed6..7f60ebcb 100644 --- a/extras/Hadrons/Global.hpp +++ b/extras/Hadrons/Global.hpp @@ -111,6 +111,7 @@ extern HadronsLogger HadronsLogWarning; extern HadronsLogger HadronsLogMessage; extern HadronsLogger HadronsLogIterative; extern HadronsLogger HadronsLogDebug; +extern HadronsLogger HadronsLogIRL; void initLogger(void); diff --git a/extras/Hadrons/LanczosUtils.hpp b/extras/Hadrons/LanczosUtils.hpp new file mode 100644 index 00000000..e6a78753 --- /dev/null +++ b/extras/Hadrons/LanczosUtils.hpp @@ -0,0 +1,88 @@ +#ifndef Hadrons_LanczosUtils_hpp_ +#define Hadrons_LanczosUtils_hpp_ + +#include +#include + +BEGIN_HADRONS_NAMESPACE + +// Lanczos type +#ifndef HADRONS_DEFAULT_LANCZOS_NBASIS +#define HADRONS_DEFAULT_LANCZOS_NBASIS 60 +#endif + +template +using LCL = LocalCoherenceLanczos; + +template +struct EigenPack +{ + std::vector eval; + std::vector evec; + + EigenPack(void) = default; + + EigenPack(const size_t size, GridBase *grid) + { + resize(size, grid); + } + + void resize(const size_t size, GridBase *grid) + { + eval.resize(size); + evec.resize(size, grid); + } + + void read(const std::string fileStem) + { + std::string evecFilename = fileStem + "_evec.bin"; + std::string evalFilename = fileStem + "_eval.xml"; + emptyUserRecord record; + ScidacReader binReader; + XmlReader xmlReader(evalFilename); + + LOG(Message) << "Reading " << evec.size() << " eigenvectors from '" + << evecFilename << "'" << std::endl; + binReader.open(evecFilename); + for(int k = 0; k < evec.size(); ++k) + { + binReader.readScidacFieldRecord(evec[k], record); + } + binReader.close(); + LOG(Message) << "Reading " << eval.size() << " eigenvalues from '" + << evalFilename << "'" << std::endl; + Grid::read(xmlReader, "evals", eval); + } + + void write(const std::string fileStem) + { + std::string evecFilename = fileStem + "_evec.bin"; + std::string evalFilename = fileStem + "_eval.xml"; + emptyUserRecord record; + ScidacWriter binWriter; + XmlWriter xmlWriter(evalFilename); + + LOG(Message) << "Writing " << evec.size() << " eigenvectors to '" + << evecFilename << "'" << std::endl; + binWriter.open(fileStem + "_evec.bin"); + for(int k = 0; k < evec.size(); ++k) + { + binWriter.writeScidacFieldRecord(evec[k], record); + } + binWriter.close(); + LOG(Message) << "Writing " << eval.size() << " eigenvalues to '" + << evalFilename << "'" << std::endl; + Grid::write(xmlWriter, "evals", eval); + } +}; + +template +using FineEigenPack = EigenPack; + +template +using CoarseEigenPack = EigenPack::CoarseField>; + +END_HADRONS_NAMESPACE + +#endif // Hadrons_LanczosUtils_hpp_ \ No newline at end of file diff --git a/extras/Hadrons/Makefile.am b/extras/Hadrons/Makefile.am index 3d07679a..477b93e4 100644 --- a/extras/Hadrons/Makefile.am +++ b/extras/Hadrons/Makefile.am @@ -21,6 +21,7 @@ nobase_libHadrons_a_HEADERS = \ GeneticScheduler.hpp \ Global.hpp \ Graph.hpp \ + LanczosUtils.hpp \ Module.hpp \ Modules.hpp \ ModuleFactory.hpp \ diff --git a/extras/Hadrons/Modules/MSolver/LocalCoherenceLanczos.hpp b/extras/Hadrons/Modules/MSolver/LocalCoherenceLanczos.hpp index 5fd486fd..a41a7839 100644 --- a/extras/Hadrons/Modules/MSolver/LocalCoherenceLanczos.hpp +++ b/extras/Hadrons/Modules/MSolver/LocalCoherenceLanczos.hpp @@ -31,11 +31,7 @@ See the full license in the file "LICENSE" in the top level distribution directo #include #include #include -#include - -#ifndef DEFAULT_LANCZOS_NBASIS -#define DEFAULT_LANCZOS_NBASIS 60 -#endif +#include BEGIN_HADRONS_NAMESPACE @@ -55,7 +51,8 @@ public: LanczosParams, coarseParams, ChebyParams, smoother, RealD, coarseRelaxTol, - std::string, blockSize); + std::string, blockSize, + std::string, output); }; template @@ -63,9 +60,9 @@ class TLocalCoherenceLanczos: public Module { public: FERM_TYPE_ALIASES(FImpl,); - typedef LocalCoherenceLanczos LCL; - typedef typename LCL::FineField FineField; - typedef typename LCL::CoarseField CoarseField; + typedef LCL LCL; + typedef FineEigenPack FineEigenPack; + typedef CoarseEigenPack CoarseEigenPack; typedef SchurDiagMooeeOperator SchurFMat; public: // constructor @@ -88,12 +85,12 @@ private: std::unique_ptr coarseGrid_{nullptr}; std::unique_ptr coarseGrid4Rb_{nullptr}; std::unique_ptr coarseGridRb_{nullptr}; - std::string fevecName_, cevecName_, fevalName_, cevalName_; + std::string fineName_, coarseName_; }; MODULE_REGISTER_NS(LocalCoherenceLanczos, - ARG(TLocalCoherenceLanczos), - MSolver); + ARG(TLocalCoherenceLanczos), + MSolver); /****************************************************************************** * TLocalCoherenceLanczos implementation * @@ -103,10 +100,8 @@ template TLocalCoherenceLanczos::TLocalCoherenceLanczos(const std::string name) : Module(name) { - fevecName_ = getName() + "_fineEvec"; - cevecName_ = getName() + "_coarseEvec"; - fevalName_ = getName() + "_fineEval"; - cevalName_ = getName() + "_coarseEval"; + fineName_ = getName() + "_fine"; + coarseName_ = getName() + "_coarse"; } // dependencies/products /////////////////////////////////////////////////////// @@ -121,7 +116,7 @@ std::vector TLocalCoherenceLanczos::getInput(void) template std::vector TLocalCoherenceLanczos::getOutput(void) { - std::vector out = {fevecName_, cevecName_, fevalName_, cevalName_}; + std::vector out = {fineName_, coarseName_}; return out; } @@ -187,47 +182,56 @@ void TLocalCoherenceLanczos::setup(void) { makeCoarseGrid(); } - envCreate(std::vector, fevecName_, Ls_, par().fineParams.Nm, + envCreate(FineEigenPack, fineName_, Ls_, par().fineParams.Nm, env().getRbGrid(Ls_)); - envCreate(std::vector, cevecName_, Ls_, par().coarseParams.Nm, + envCreate(CoarseEigenPack, coarseName_, Ls_, par().coarseParams.Nm, coarseGridRb_.get()); - envCreate(std::vector, fevalName_, Ls_, par().fineParams.Nm); - envCreate(std::vector, cevalName_, Ls_, par().coarseParams.Nm); + auto &fine = envGet(FineEigenPack, fineName_); + auto &coarse = envGet(CoarseEigenPack, coarseName_); envTmp(SchurFMat, "mat", Ls_, envGet(FMat, par().action)); envGetTmp(SchurFMat, mat); - envTmp(LCL, "solver", Ls_, env().getRbGrid(Ls_), coarseGridRb_.get(), mat, Odd, - envGet(std::vector, fevecName_), - envGet(std::vector, cevecName_), - envGet(std::vector, fevalName_), - envGet(std::vector, cevalName_)); + envTmp(LCL, "solver", Ls_, env().getRbGrid(Ls_), coarseGridRb_.get(), mat, + Odd, fine.evec, coarse.evec, fine.eval, coarse.eval); } // execution /////////////////////////////////////////////////////////////////// template void TLocalCoherenceLanczos::execute(void) { - auto &fine = par().fineParams; - auto &coarse = par().coarseParams; - + auto &finePar = par().fineParams; + auto &coarsePar = par().coarseParams; + auto &fine = envGet(FineEigenPack, fineName_); + auto &coarse = envGet(CoarseEigenPack, coarseName_); + envGetTmp(LCL, solver); if (par().doFine) { LOG(Message) << "Performing fine grid IRL -- Nstop= " - << fine.Nstop << ", Nk= " << fine.Nk << ", Nm= " - << fine.Nm << std::endl; - solver.calcFine(fine.Cheby, fine.Nstop, fine.Nk, fine.Nm, fine.resid, - fine.MaxIt, fine.betastp, fine.MinRes); + << finePar.Nstop << ", Nk= " << finePar.Nk << ", Nm= " + << finePar.Nm << std::endl; + solver.calcFine(finePar.Cheby, finePar.Nstop, finePar.Nk, finePar.Nm, + finePar.resid,finePar.MaxIt, finePar.betastp, + finePar.MinRes); LOG(Message) << "Orthogonalising" << std::endl; solver.Orthogonalise(); + if (!par().output.empty()) + { + fine.write(par().output + "_fine"); + } } if (par().doCoarse) { LOG(Message) << "Performing coarse grid IRL -- Nstop= " - << fine.Nstop << ", Nk= " << fine.Nk << ", Nm= " - << fine.Nm << std::endl; - solver.calcCoarse(coarse.Cheby, par().smoother, par().coarseRelaxTol, - coarse.Nstop, coarse.Nk, coarse.Nm, coarse.resid, - coarse.MaxIt, coarse.betastp, coarse.MinRes); + << coarsePar.Nstop << ", Nk= " << coarsePar.Nk << ", Nm= " + << coarsePar.Nm << std::endl; + solver.calcCoarse(coarsePar.Cheby, par().smoother, par().coarseRelaxTol, + coarsePar.Nstop, coarsePar.Nk, coarsePar.Nm, + coarsePar.resid, coarsePar.MaxIt, coarsePar.betastp, + coarsePar.MinRes); + if (!par().output.empty()) + { + coarse.write(par().output + "_coarse"); + } } } From 48282260956b914eae53ab0b32774bbbfc4e4334 Mon Sep 17 00:00:00 2001 From: Antonin Portelli Date: Tue, 27 Feb 2018 14:43:51 +0000 Subject: [PATCH 11/38] Hadrons: prettier log --- extras/Hadrons/Global.cc | 17 ++++++++++++++++- lib/log/Log.h | 10 ++++++++-- 2 files changed, 24 insertions(+), 3 deletions(-) diff --git a/extras/Hadrons/Global.cc b/extras/Hadrons/Global.cc index b121434f..0164a1ab 100644 --- a/extras/Hadrons/Global.cc +++ b/extras/Hadrons/Global.cc @@ -41,19 +41,34 @@ HadronsLogger Hadrons::HadronsLogIRL(1,"IRL"); void Hadrons::initLogger(void) { - auto w = std::string("Hadrons").length(); + auto w = std::string("Hadrons").length(); + int cw = 8; + + GridLogError.setTopWidth(w); GridLogWarning.setTopWidth(w); GridLogMessage.setTopWidth(w); GridLogIterative.setTopWidth(w); GridLogDebug.setTopWidth(w); GridLogIRL.setTopWidth(w); + GridLogError.setChanWidth(cw); + GridLogWarning.setChanWidth(cw); + GridLogMessage.setChanWidth(cw); + GridLogIterative.setChanWidth(cw); + GridLogDebug.setChanWidth(cw); + GridLogIRL.setChanWidth(cw); HadronsLogError.Active(GridLogError.isActive()); HadronsLogWarning.Active(GridLogWarning.isActive()); HadronsLogMessage.Active(GridLogMessage.isActive()); HadronsLogIterative.Active(GridLogIterative.isActive()); HadronsLogDebug.Active(GridLogDebug.isActive()); HadronsLogIRL.Active(GridLogIRL.isActive()); + HadronsLogError.setChanWidth(cw); + HadronsLogWarning.setChanWidth(cw); + HadronsLogMessage.setChanWidth(cw); + HadronsLogIterative.setChanWidth(cw); + HadronsLogDebug.setChanWidth(cw); + HadronsLogIRL.setChanWidth(cw); } // type utilities ////////////////////////////////////////////////////////////// diff --git a/lib/log/Log.h b/lib/log/Log.h index 011a7250..b58c5d16 100644 --- a/lib/log/Log.h +++ b/lib/log/Log.h @@ -86,7 +86,7 @@ protected: Colours &Painter; int active; int timing_mode; - int topWidth{-1}; + int topWidth{-1}, chanWidth{-1}; static int timestamp; std::string name, topName; std::string COLOUR; @@ -126,6 +126,7 @@ public: } } void setTopWidth(const int w) {topWidth = w;} + void setChanWidth(const int w) {chanWidth = w;} friend std::ostream& operator<< (std::ostream& stream, Logger& log){ @@ -136,7 +137,12 @@ public: stream << std::setw(log.topWidth); } stream << log.topName << log.background()<< " : "; - stream << log.colour() << std::left << log.name << log.background() << " : "; + stream << log.colour() << std::left; + if (log.chanWidth > 0) + { + stream << std::setw(log.chanWidth); + } + stream << log.name << log.background() << " : "; if ( log.timestamp ) { log.StopWatch->Stop(); GridTime now = log.StopWatch->Elapsed(); From 16ebbfff2991b5bbaee5d2fa337dce0540ebb978 Mon Sep 17 00:00:00 2001 From: Antonin Portelli Date: Tue, 27 Feb 2018 18:45:23 +0000 Subject: [PATCH 12/38] Hadrons: Schur convention globally defined through a macro --- extras/Hadrons/Global.hpp | 12 ++++++++++++ .../Modules/MSolver/LocalCoherenceLanczos.hpp | 8 ++++---- extras/Hadrons/Modules/MSolver/RBPrecCG.hpp | 13 ++++++++----- 3 files changed, 24 insertions(+), 9 deletions(-) diff --git a/extras/Hadrons/Global.hpp b/extras/Hadrons/Global.hpp index 7f60ebcb..675a97c2 100644 --- a/extras/Hadrons/Global.hpp +++ b/extras/Hadrons/Global.hpp @@ -181,6 +181,18 @@ typedef XmlWriter ResultWriter; #define RESULT_FILE_NAME(name) \ name + "." + std::to_string(vm().getTrajectory()) + "." + resultFileExt +// default Schur convention + +#ifndef HADRONS_DEFAULT_SCHUR +#define HADRONS_DEFAULT_SCHUR DiagMooee +#endif +#define _HADRONS_SCHUR_OP_(conv) Schur##conv##Operator +#define HADRONS_SCHUR_OP(conv) _HADRONS_SCHUR_OP_(conv) +#define HADRONS_DEFAULT_SCHUR_OP HADRONS_SCHUR_OP(HADRONS_DEFAULT_SCHUR) +#define _HADRONS_SCHUR_SOLVE_(conv) SchurRedBlack##conv##Solve +#define HADRONS_SCHUR_SOLVE(conv) _HADRONS_SCHUR_SOLVE_(conv) +#define HADRONS_DEFAULT_SCHUR_SOLVE HADRONS_SCHUR_SOLVE(HADRONS_DEFAULT_SCHUR) + END_HADRONS_NAMESPACE #include diff --git a/extras/Hadrons/Modules/MSolver/LocalCoherenceLanczos.hpp b/extras/Hadrons/Modules/MSolver/LocalCoherenceLanczos.hpp index a41a7839..387b9455 100644 --- a/extras/Hadrons/Modules/MSolver/LocalCoherenceLanczos.hpp +++ b/extras/Hadrons/Modules/MSolver/LocalCoherenceLanczos.hpp @@ -60,10 +60,10 @@ class TLocalCoherenceLanczos: public Module { public: FERM_TYPE_ALIASES(FImpl,); - typedef LCL LCL; - typedef FineEigenPack FineEigenPack; - typedef CoarseEigenPack CoarseEigenPack; - typedef SchurDiagMooeeOperator SchurFMat; + typedef LCL LCL; + typedef FineEigenPack FineEigenPack; + typedef CoarseEigenPack CoarseEigenPack; + typedef HADRONS_DEFAULT_SCHUR_OP SchurFMat; public: // constructor TLocalCoherenceLanczos(const std::string name); diff --git a/extras/Hadrons/Modules/MSolver/RBPrecCG.hpp b/extras/Hadrons/Modules/MSolver/RBPrecCG.hpp index 54c0f2d8..77ce6cd4 100644 --- a/extras/Hadrons/Modules/MSolver/RBPrecCG.hpp +++ b/extras/Hadrons/Modules/MSolver/RBPrecCG.hpp @@ -43,9 +43,10 @@ BEGIN_MODULE_NAMESPACE(MSolver) class RBPrecCGPar: Serializable { public: - GRID_SERIALIZABLE_CLASS_MEMBERS(RBPrecCGPar, - std::string, action, - double , residual); + GRID_SERIALIZABLE_CLASS_MEMBERS(RBPrecCGPar , + std::string , action, + unsigned int , maxIteration, + double , residual); }; template @@ -117,14 +118,16 @@ void TRBPrecCG::setup(void) auto &mat = envGet(FMat, par().action); auto solver = [&mat, this](FermionField &sol, const FermionField &source) { - ConjugateGradient cg(par().residual, 10000); - SchurRedBlackDiagMooeeSolve schurSolver(cg); + ConjugateGradient cg(par().residual, + par().maxIteration); + HADRONS_DEFAULT_SCHUR_SOLVE schurSolver(cg); schurSolver(mat, source, sol); }; envCreate(SolverFn, getName(), Ls, solver); } + // execution /////////////////////////////////////////////////////////////////// template void TRBPrecCG::execute(void) From abb7d4d2f5026c348342e2d19c2f479a8a5b0755 Mon Sep 17 00:00:00 2001 From: Antonin Portelli Date: Tue, 27 Feb 2018 19:32:19 +0000 Subject: [PATCH 13/38] Hadrons: z-Mobius action --- extras/Hadrons/Global.hpp | 3 + extras/Hadrons/Modules.hpp | 31 +---- extras/Hadrons/Modules/MAction/ZMobiusDWF.hpp | 116 ++++++++++++++++++ .../Modules/MSolver/LocalCoherenceLanczos.hpp | 3 + extras/Hadrons/Modules/MSolver/RBPrecCG.hpp | 3 +- extras/Hadrons/modules.inc | 1 + 6 files changed, 126 insertions(+), 31 deletions(-) create mode 100644 extras/Hadrons/Modules/MAction/ZMobiusDWF.hpp diff --git a/extras/Hadrons/Global.hpp b/extras/Hadrons/Global.hpp index 675a97c2..ed8f4f32 100644 --- a/extras/Hadrons/Global.hpp +++ b/extras/Hadrons/Global.hpp @@ -58,6 +58,9 @@ using Grid::operator<<; #ifndef FIMPL #define FIMPL WilsonImplR #endif +#ifndef ZFIMPL +#define ZFIMPL ZWilsonImplR +#endif #ifndef SIMPL #define SIMPL ScalarImplCR #endif diff --git a/extras/Hadrons/Modules.hpp b/extras/Hadrons/Modules.hpp index 1334a9b4..1b78a85d 100644 --- a/extras/Hadrons/Modules.hpp +++ b/extras/Hadrons/Modules.hpp @@ -1,33 +1,3 @@ -/************************************************************************************* - -Grid physics library, www.github.com/paboyle/Grid - -Source file: extras/Hadrons/Modules.hpp - -Copyright (C) 2015-2018 - -Author: Antonin Portelli -Author: Guido Cossu -Author: Lanny91 -Author: pretidav - -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 @@ -60,6 +30,7 @@ See the full license in the file "LICENSE" in the top level distribution directo #include #include #include +#include #include #include #include diff --git a/extras/Hadrons/Modules/MAction/ZMobiusDWF.hpp b/extras/Hadrons/Modules/MAction/ZMobiusDWF.hpp new file mode 100644 index 00000000..b0c614df --- /dev/null +++ b/extras/Hadrons/Modules/MAction/ZMobiusDWF.hpp @@ -0,0 +1,116 @@ +#ifndef Hadrons_MAction_ZMobiusDWF_hpp_ +#define Hadrons_MAction_ZMobiusDWF_hpp_ + +#include +#include +#include + +BEGIN_HADRONS_NAMESPACE + +/****************************************************************************** + * ZMobiusDWF * + ******************************************************************************/ +BEGIN_MODULE_NAMESPACE(MAction) + +class ZMobiusDWFPar: Serializable +{ +public: + GRID_SERIALIZABLE_CLASS_MEMBERS(ZMobiusDWFPar, + std::string , gauge, + unsigned int , Ls, + double , mass, + double , M5, + double , b, + double , c, + std::vector>, omega, + std::string , boundary); +}; + +template +class TZMobiusDWF: public Module +{ +public: + FGS_TYPE_ALIASES(FImpl,); +public: + // constructor + TZMobiusDWF(const std::string name); + // destructor + virtual ~TZMobiusDWF(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(ZMobiusDWF, TZMobiusDWF, MAction); + +/****************************************************************************** + * TZMobiusDWF implementation * + ******************************************************************************/ +// constructor ///////////////////////////////////////////////////////////////// +template +TZMobiusDWF::TZMobiusDWF(const std::string name) +: Module(name) +{} + +// dependencies/products /////////////////////////////////////////////////////// +template +std::vector TZMobiusDWF::getInput(void) +{ + std::vector in = {par().gauge}; + + return in; +} + +template +std::vector TZMobiusDWF::getOutput(void) +{ + std::vector out = {getName()}; + + return out; +} + +// setup /////////////////////////////////////////////////////////////////////// +template +void TZMobiusDWF::setup(void) +{ + LOG(Message) << "Setting up z-Mobius domain wall fermion matrix with m= " + << par().mass << ", M5= " << par().M5 << ", Ls= " + << ", b= " << par().b << ", c= " << par().c + << par().Ls << " using gauge field '" << par().gauge << "'" + << std::endl; + LOG(Message) << "Omegas: " << std::endl; + for (unsigned int i = 0; i < par().omega.size(); ++i) + { + LOG(Message) << " omega[" << i << "]= " << par().omega[i] << std::endl; + } + LOG(Message) << "Fermion boundary conditions: " << par().boundary + << std::endl; + + env().createGrid(par().Ls); + auto &U = envGet(LatticeGaugeField, par().gauge); + auto &g4 = *env().getGrid(); + auto &grb4 = *env().getRbGrid(); + auto &g5 = *env().getGrid(par().Ls); + auto &grb5 = *env().getRbGrid(par().Ls); + auto omega = par().omega; + std::vector boundary = strToVec(par().boundary); + typename ZMobiusFermion::ImplParams implParams(boundary); + envCreateDerived(FMat, ZMobiusFermion, getName(), par().Ls, U, g5, + grb5, g4, grb4, par().mass, par().M5, omega, + par().b, par().c, implParams); +} + +// execution /////////////////////////////////////////////////////////////////// +template +void TZMobiusDWF::execute(void) +{} + +END_MODULE_NAMESPACE + +END_HADRONS_NAMESPACE + +#endif // Hadrons_MAction_ZMobiusDWF_hpp_ diff --git a/extras/Hadrons/Modules/MSolver/LocalCoherenceLanczos.hpp b/extras/Hadrons/Modules/MSolver/LocalCoherenceLanczos.hpp index 387b9455..feb1d0be 100644 --- a/extras/Hadrons/Modules/MSolver/LocalCoherenceLanczos.hpp +++ b/extras/Hadrons/Modules/MSolver/LocalCoherenceLanczos.hpp @@ -91,6 +91,9 @@ private: MODULE_REGISTER_NS(LocalCoherenceLanczos, ARG(TLocalCoherenceLanczos), MSolver); +MODULE_REGISTER_NS(ZLocalCoherenceLanczos, + ARG(TLocalCoherenceLanczos), + MSolver); /****************************************************************************** * TLocalCoherenceLanczos implementation * diff --git a/extras/Hadrons/Modules/MSolver/RBPrecCG.hpp b/extras/Hadrons/Modules/MSolver/RBPrecCG.hpp index 77ce6cd4..2b914625 100644 --- a/extras/Hadrons/Modules/MSolver/RBPrecCG.hpp +++ b/extras/Hadrons/Modules/MSolver/RBPrecCG.hpp @@ -70,7 +70,8 @@ protected: virtual void execute(void); }; -MODULE_REGISTER_NS(RBPrecCG, TRBPrecCG, MSolver); +MODULE_REGISTER_NS(RBPrecCG, TRBPrecCG, MSolver); +MODULE_REGISTER_NS(ZRBPrecCG, TRBPrecCG, MSolver); /****************************************************************************** * TRBPrecCG template implementation * diff --git a/extras/Hadrons/modules.inc b/extras/Hadrons/modules.inc index 4a9b2630..79db9396 100644 --- a/extras/Hadrons/modules.inc +++ b/extras/Hadrons/modules.inc @@ -43,6 +43,7 @@ modules_hpp =\ Modules/MAction/DWF.hpp \ Modules/MAction/Wilson.hpp \ Modules/MAction/WilsonClover.hpp \ + Modules/MAction/ZMobiusDWF.hpp \ Modules/MScalarSUN/Div.hpp \ Modules/MScalarSUN/TrMag.hpp \ Modules/MScalarSUN/TwoPoint.hpp \ From 6ec42b4b82311b6eaa20e37dec49a17af6df0f95 Mon Sep 17 00:00:00 2001 From: Antonin Portelli Date: Thu, 1 Mar 2018 12:27:29 +0000 Subject: [PATCH 14/38] LCL: external storage fix --- lib/algorithms/iterative/LocalCoherenceLanczos.h | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/lib/algorithms/iterative/LocalCoherenceLanczos.h b/lib/algorithms/iterative/LocalCoherenceLanczos.h index b8348c0c..3f997a01 100644 --- a/lib/algorithms/iterative/LocalCoherenceLanczos.h +++ b/lib/algorithms/iterative/LocalCoherenceLanczos.h @@ -373,14 +373,14 @@ public: RealD MaxIt, RealD betastp, int MinRes) { Chebyshev Cheby(cheby_op); - ProjectedHermOp Op(_FineOp,_subspace); - ProjectedFunctionHermOp ChebyOp (Cheby,_FineOp,_subspace); + ProjectedHermOp Op(_FineOp,subspace); + ProjectedFunctionHermOp ChebyOp (Cheby,_FineOp,subspace); ////////////////////////////////////////////////////////////////////////////////////////////////// // create a smoother and see if we can get a cheap convergence test and smooth inside the IRL ////////////////////////////////////////////////////////////////////////////////////////////////// Chebyshev ChebySmooth(cheby_smooth); - ImplicitlyRestartedLanczosSmoothedTester ChebySmoothTester(ChebyOp,ChebySmooth,_FineOp,_subspace,relax); + ImplicitlyRestartedLanczosSmoothedTester ChebySmoothTester(ChebyOp,ChebySmooth,_FineOp,subspace,relax); evals_coarse.resize(Nm); evec_coarse.resize(Nm,_CoarseGrid); From b8b05f143fa7bdf66d2b5bbeed193b74b32037fa Mon Sep 17 00:00:00 2001 From: Antonin Portelli Date: Thu, 1 Mar 2018 12:53:16 +0000 Subject: [PATCH 15/38] Hadrons: Lanczos more conservative type names --- extras/Hadrons/LanczosUtils.hpp | 12 +++++----- .../Modules/MSolver/LocalCoherenceLanczos.hpp | 22 +++++++++---------- 2 files changed, 17 insertions(+), 17 deletions(-) diff --git a/extras/Hadrons/LanczosUtils.hpp b/extras/Hadrons/LanczosUtils.hpp index e6a78753..ab3ddbd5 100644 --- a/extras/Hadrons/LanczosUtils.hpp +++ b/extras/Hadrons/LanczosUtils.hpp @@ -11,16 +11,13 @@ BEGIN_HADRONS_NAMESPACE #define HADRONS_DEFAULT_LANCZOS_NBASIS 60 #endif -template -using LCL = LocalCoherenceLanczos; - template struct EigenPack { + typedef T VectorType; std::vector eval; std::vector evec; - + EigenPack(void) = default; EigenPack(const size_t size, GridBase *grid) @@ -81,7 +78,10 @@ template using FineEigenPack = EigenPack; template -using CoarseEigenPack = EigenPack::CoarseField>; +using CoarseEigenPack = EigenPack< + typename LocalCoherenceLanczos::CoarseField>; END_HADRONS_NAMESPACE diff --git a/extras/Hadrons/Modules/MSolver/LocalCoherenceLanczos.hpp b/extras/Hadrons/Modules/MSolver/LocalCoherenceLanczos.hpp index feb1d0be..b708d59d 100644 --- a/extras/Hadrons/Modules/MSolver/LocalCoherenceLanczos.hpp +++ b/extras/Hadrons/Modules/MSolver/LocalCoherenceLanczos.hpp @@ -60,9 +60,11 @@ class TLocalCoherenceLanczos: public Module { public: FERM_TYPE_ALIASES(FImpl,); - typedef LCL LCL; - typedef FineEigenPack FineEigenPack; - typedef CoarseEigenPack CoarseEigenPack; + typedef LocalCoherenceLanczos LCL; + typedef FineEigenPack FinePack; + typedef CoarseEigenPack CoarsePack; typedef HADRONS_DEFAULT_SCHUR_OP SchurFMat; public: // constructor @@ -185,12 +187,10 @@ void TLocalCoherenceLanczos::setup(void) { makeCoarseGrid(); } - envCreate(FineEigenPack, fineName_, Ls_, par().fineParams.Nm, - env().getRbGrid(Ls_)); - envCreate(CoarseEigenPack, coarseName_, Ls_, par().coarseParams.Nm, - coarseGridRb_.get()); - auto &fine = envGet(FineEigenPack, fineName_); - auto &coarse = envGet(CoarseEigenPack, coarseName_); + envCreate(FinePack, fineName_, Ls_, par().fineParams.Nm, env().getRbGrid(Ls_)); + envCreate(CoarsePack, coarseName_, Ls_, par().coarseParams.Nm, coarseGridRb_.get()); + auto &fine = envGet(FinePack, fineName_); + auto &coarse = envGet(CoarsePack, coarseName_); envTmp(SchurFMat, "mat", Ls_, envGet(FMat, par().action)); envGetTmp(SchurFMat, mat); envTmp(LCL, "solver", Ls_, env().getRbGrid(Ls_), coarseGridRb_.get(), mat, @@ -203,8 +203,8 @@ void TLocalCoherenceLanczos::execute(void) { auto &finePar = par().fineParams; auto &coarsePar = par().coarseParams; - auto &fine = envGet(FineEigenPack, fineName_); - auto &coarse = envGet(CoarseEigenPack, coarseName_); + auto &fine = envGet(FinePack, fineName_); + auto &coarse = envGet(CoarsePack, coarseName_); envGetTmp(LCL, solver); if (par().doFine) From e418b044f747f0375cdf420699a3c026f5f82d1b Mon Sep 17 00:00:00 2001 From: Antonin Portelli Date: Thu, 1 Mar 2018 12:57:28 +0000 Subject: [PATCH 16/38] Hadrons: code cleaning --- extras/Hadrons/LanczosUtils.hpp | 27 +++++++++++++++++ extras/Hadrons/Modules.hpp | 30 +++++++++++++++++++ .../Hadrons/Modules/MAction/WilsonClover.hpp | 1 + extras/Hadrons/Modules/MAction/ZMobiusDWF.hpp | 27 +++++++++++++++++ extras/Hadrons/Modules/MGauge/FundtoHirep.hpp | 1 + extras/Hadrons/Modules/MIO/LoadNersc.cc | 2 -- .../Modules/MSolver/LocalCoherenceLanczos.hpp | 18 +++++------ 7 files changed, 95 insertions(+), 11 deletions(-) diff --git a/extras/Hadrons/LanczosUtils.hpp b/extras/Hadrons/LanczosUtils.hpp index ab3ddbd5..f8e7693a 100644 --- a/extras/Hadrons/LanczosUtils.hpp +++ b/extras/Hadrons/LanczosUtils.hpp @@ -1,3 +1,30 @@ +/************************************************************************************* + +Grid physics library, www.github.com/paboyle/Grid + +Source file: extras/Hadrons/LanczosUtils.hpp + +Copyright (C) 2015-2018 + +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 */ #ifndef Hadrons_LanczosUtils_hpp_ #define Hadrons_LanczosUtils_hpp_ diff --git a/extras/Hadrons/Modules.hpp b/extras/Hadrons/Modules.hpp index 1b78a85d..57fa93a3 100644 --- a/extras/Hadrons/Modules.hpp +++ b/extras/Hadrons/Modules.hpp @@ -1,3 +1,33 @@ +/************************************************************************************* + +Grid physics library, www.github.com/paboyle/Grid + +Source file: extras/Hadrons/Modules.hpp + +Copyright (C) 2015-2018 + +Author: Antonin Portelli +Author: Guido Cossu +Author: Lanny91 +Author: pretidav + +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 diff --git a/extras/Hadrons/Modules/MAction/WilsonClover.hpp b/extras/Hadrons/Modules/MAction/WilsonClover.hpp index d5174c6d..2967931c 100644 --- a/extras/Hadrons/Modules/MAction/WilsonClover.hpp +++ b/extras/Hadrons/Modules/MAction/WilsonClover.hpp @@ -6,6 +6,7 @@ Source file: extras/Hadrons/Modules/MAction/WilsonClover.hpp Copyright (C) 2015-2018 +Author: Antonin Portelli Author: Guido Cossu Author: pretidav diff --git a/extras/Hadrons/Modules/MAction/ZMobiusDWF.hpp b/extras/Hadrons/Modules/MAction/ZMobiusDWF.hpp index b0c614df..fd812a2a 100644 --- a/extras/Hadrons/Modules/MAction/ZMobiusDWF.hpp +++ b/extras/Hadrons/Modules/MAction/ZMobiusDWF.hpp @@ -1,3 +1,30 @@ +/************************************************************************************* + +Grid physics library, www.github.com/paboyle/Grid + +Source file: extras/Hadrons/Modules/MAction/ZMobiusDWF.hpp + +Copyright (C) 2015-2018 + +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 */ #ifndef Hadrons_MAction_ZMobiusDWF_hpp_ #define Hadrons_MAction_ZMobiusDWF_hpp_ diff --git a/extras/Hadrons/Modules/MGauge/FundtoHirep.hpp b/extras/Hadrons/Modules/MGauge/FundtoHirep.hpp index 4fba59f0..508145bc 100644 --- a/extras/Hadrons/Modules/MGauge/FundtoHirep.hpp +++ b/extras/Hadrons/Modules/MGauge/FundtoHirep.hpp @@ -6,6 +6,7 @@ Source file: extras/Hadrons/Modules/MGauge/FundtoHirep.hpp Copyright (C) 2015-2018 +Author: Antonin Portelli Author: pretidav This program is free software; you can redistribute it and/or modify diff --git a/extras/Hadrons/Modules/MIO/LoadNersc.cc b/extras/Hadrons/Modules/MIO/LoadNersc.cc index f20606fc..36a56d54 100644 --- a/extras/Hadrons/Modules/MIO/LoadNersc.cc +++ b/extras/Hadrons/Modules/MIO/LoadNersc.cc @@ -71,6 +71,4 @@ void TLoadNersc::execute(void) auto &U = envGet(LatticeGaugeField, getName()); NerscIO::readConfiguration(U, header, fileName); - LOG(Message) << "NERSC header:" << std::endl; - dump_meta_data(header, LOG(Message)); } diff --git a/extras/Hadrons/Modules/MSolver/LocalCoherenceLanczos.hpp b/extras/Hadrons/Modules/MSolver/LocalCoherenceLanczos.hpp index b708d59d..6ab0412a 100644 --- a/extras/Hadrons/Modules/MSolver/LocalCoherenceLanczos.hpp +++ b/extras/Hadrons/Modules/MSolver/LocalCoherenceLanczos.hpp @@ -44,15 +44,15 @@ class LocalCoherenceLanczosPar: Serializable { public: GRID_SERIALIZABLE_CLASS_MEMBERS(LocalCoherenceLanczosPar, - std::string, action, - int, doFine, - int, doCoarse, - LanczosParams, fineParams, - LanczosParams, coarseParams, - ChebyParams, smoother, - RealD, coarseRelaxTol, - std::string, blockSize, - std::string, output); + std::string, action, + int, doFine, + int, doCoarse, + LanczosParams, fineParams, + LanczosParams, coarseParams, + ChebyParams, smoother, + RealD, coarseRelaxTol, + std::string, blockSize, + std::string, output); }; template From 5b937e3644dba95b6e1c3e868fc567ccf991c2ea Mon Sep 17 00:00:00 2001 From: Antonin Portelli Date: Thu, 1 Mar 2018 17:28:38 +0000 Subject: [PATCH 17/38] Hadrons: VM memory profiling fix --- extras/Hadrons/VirtualMachine.cc | 10 ++++++++-- 1 file changed, 8 insertions(+), 2 deletions(-) diff --git a/extras/Hadrons/VirtualMachine.cc b/extras/Hadrons/VirtualMachine.cc index d47bafb7..c91667e7 100644 --- a/extras/Hadrons/VirtualMachine.cc +++ b/extras/Hadrons/VirtualMachine.cc @@ -424,11 +424,17 @@ void VirtualMachine::memoryProfile(const unsigned int address) cleanEnvironment(); for (auto &in: m->getInput()) { - memoryProfile(env().getObjectModule(in)); + if (!env().hasCreatedObject(in)) + { + memoryProfile(env().getObjectModule(in)); + } } for (auto &ref: m->getReference()) { - memoryProfile(env().getObjectModule(ref)); + if (!env().hasCreatedObject(ref)) + { + memoryProfile(env().getObjectModule(ref)); + } } m->setup(); updateProfile(address); From 4d2a32ae7ace84c8a762d5f97cd34146c8edbf7d Mon Sep 17 00:00:00 2001 From: Antonin Portelli Date: Thu, 1 Mar 2018 18:03:44 +0000 Subject: [PATCH 18/38] Hadrons: z-Mobius message fix --- extras/Hadrons/Modules/MAction/ZMobiusDWF.hpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/extras/Hadrons/Modules/MAction/ZMobiusDWF.hpp b/extras/Hadrons/Modules/MAction/ZMobiusDWF.hpp index fd812a2a..518b88f2 100644 --- a/extras/Hadrons/Modules/MAction/ZMobiusDWF.hpp +++ b/extras/Hadrons/Modules/MAction/ZMobiusDWF.hpp @@ -105,9 +105,9 @@ template void TZMobiusDWF::setup(void) { LOG(Message) << "Setting up z-Mobius domain wall fermion matrix with m= " - << par().mass << ", M5= " << par().M5 << ", Ls= " + << par().mass << ", M5= " << par().M5 << ", Ls= " << par().Ls << ", b= " << par().b << ", c= " << par().c - << par().Ls << " using gauge field '" << par().gauge << "'" + << " using gauge field '" << par().gauge << "'" << std::endl; LOG(Message) << "Omegas: " << std::endl; for (unsigned int i = 0; i < par().omega.size(); ++i) From 15767a1491a36edd02ee29eeef4efd7ac1fd85fe Mon Sep 17 00:00:00 2001 From: Antonin Portelli Date: Thu, 1 Mar 2018 18:04:08 +0000 Subject: [PATCH 19/38] Hadrons: LCL fine convergence test --- extras/Hadrons/Modules/MSolver/LocalCoherenceLanczos.hpp | 1 + 1 file changed, 1 insertion(+) diff --git a/extras/Hadrons/Modules/MSolver/LocalCoherenceLanczos.hpp b/extras/Hadrons/Modules/MSolver/LocalCoherenceLanczos.hpp index 6ab0412a..e30788bd 100644 --- a/extras/Hadrons/Modules/MSolver/LocalCoherenceLanczos.hpp +++ b/extras/Hadrons/Modules/MSolver/LocalCoherenceLanczos.hpp @@ -215,6 +215,7 @@ void TLocalCoherenceLanczos::execute(void) solver.calcFine(finePar.Cheby, finePar.Nstop, finePar.Nk, finePar.Nm, finePar.resid,finePar.MaxIt, finePar.betastp, finePar.MinRes); + solver.testFine(finePar.resid*100.0); LOG(Message) << "Orthogonalising" << std::endl; solver.Orthogonalise(); if (!par().output.empty()) From e5ea84d5311e86aab789e771fb4557d6b87ef3cc Mon Sep 17 00:00:00 2001 From: Antonin Portelli Date: Thu, 1 Mar 2018 19:33:11 +0000 Subject: [PATCH 20/38] Hadrons: LCL: orthogonalise coarse evec --- extras/Hadrons/Modules/MSolver/LocalCoherenceLanczos.hpp | 2 ++ 1 file changed, 2 insertions(+) diff --git a/extras/Hadrons/Modules/MSolver/LocalCoherenceLanczos.hpp b/extras/Hadrons/Modules/MSolver/LocalCoherenceLanczos.hpp index e30788bd..c9adc2a7 100644 --- a/extras/Hadrons/Modules/MSolver/LocalCoherenceLanczos.hpp +++ b/extras/Hadrons/Modules/MSolver/LocalCoherenceLanczos.hpp @@ -232,6 +232,8 @@ void TLocalCoherenceLanczos::execute(void) coarsePar.Nstop, coarsePar.Nk, coarsePar.Nm, coarsePar.resid, coarsePar.MaxIt, coarsePar.betastp, coarsePar.MinRes); + LOG(Message) << "Orthogonalising" << std::endl; + solver.Orthogonalise(); if (!par().output.empty()) { coarse.write(par().output + "_coarse"); From ba6db55cb0c9dd95e009420634860ad01003170c Mon Sep 17 00:00:00 2001 From: Antonin Portelli Date: Thu, 1 Mar 2018 23:30:58 +0000 Subject: [PATCH 21/38] Hadrons: reverse last commit --- extras/Hadrons/Modules/MSolver/LocalCoherenceLanczos.hpp | 1 - 1 file changed, 1 deletion(-) diff --git a/extras/Hadrons/Modules/MSolver/LocalCoherenceLanczos.hpp b/extras/Hadrons/Modules/MSolver/LocalCoherenceLanczos.hpp index c9adc2a7..ae289ddf 100644 --- a/extras/Hadrons/Modules/MSolver/LocalCoherenceLanczos.hpp +++ b/extras/Hadrons/Modules/MSolver/LocalCoherenceLanczos.hpp @@ -233,7 +233,6 @@ void TLocalCoherenceLanczos::execute(void) coarsePar.resid, coarsePar.MaxIt, coarsePar.betastp, coarsePar.MinRes); LOG(Message) << "Orthogonalising" << std::endl; - solver.Orthogonalise(); if (!par().output.empty()) { coarse.write(par().output + "_coarse"); From c4274e1660ea529e5c9e8c396677ca61c8b11537 Mon Sep 17 00:00:00 2001 From: Antonin Portelli Date: Fri, 2 Mar 2018 10:18:33 +0000 Subject: [PATCH 22/38] Hadrons: LCL cleaning --- extras/Hadrons/Modules/MSolver/LocalCoherenceLanczos.hpp | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/extras/Hadrons/Modules/MSolver/LocalCoherenceLanczos.hpp b/extras/Hadrons/Modules/MSolver/LocalCoherenceLanczos.hpp index ae289ddf..3a299bd6 100644 --- a/extras/Hadrons/Modules/MSolver/LocalCoherenceLanczos.hpp +++ b/extras/Hadrons/Modules/MSolver/LocalCoherenceLanczos.hpp @@ -174,6 +174,7 @@ void TLocalCoherenceLanczos::makeCoarseGrid(void) GridDefaultMpi())); coarseGridRb_.reset(SpaceTimeGrid::makeFourDimRedBlackGrid(coarseGrid_.get())); } + LOG(Message) << "Coarse grid: " << coarseGrid_->GlobalDimensions() << std::endl; } template @@ -232,7 +233,8 @@ void TLocalCoherenceLanczos::execute(void) coarsePar.Nstop, coarsePar.Nk, coarsePar.Nm, coarsePar.resid, coarsePar.MaxIt, coarsePar.betastp, coarsePar.MinRes); - LOG(Message) << "Orthogonalising" << std::endl; + solver.testCoarse(coarsePar.resid*100.0, par().smoother, + par().coarseRelaxTol); if (!par().output.empty()) { coarse.write(par().output + "_coarse"); From 83a101db8332a0b26958635614293be20a661991 Mon Sep 17 00:00:00 2001 From: Antonin Portelli Date: Fri, 2 Mar 2018 11:05:02 +0000 Subject: [PATCH 23/38] Hadrons: more LCL fixes --- extras/Hadrons/Modules/MSolver/LocalCoherenceLanczos.hpp | 2 +- lib/algorithms/iterative/LocalCoherenceLanczos.h | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/extras/Hadrons/Modules/MSolver/LocalCoherenceLanczos.hpp b/extras/Hadrons/Modules/MSolver/LocalCoherenceLanczos.hpp index 3a299bd6..284b233b 100644 --- a/extras/Hadrons/Modules/MSolver/LocalCoherenceLanczos.hpp +++ b/extras/Hadrons/Modules/MSolver/LocalCoherenceLanczos.hpp @@ -174,7 +174,6 @@ void TLocalCoherenceLanczos::makeCoarseGrid(void) GridDefaultMpi())); coarseGridRb_.reset(SpaceTimeGrid::makeFourDimRedBlackGrid(coarseGrid_.get())); } - LOG(Message) << "Coarse grid: " << coarseGrid_->GlobalDimensions() << std::endl; } template @@ -188,6 +187,7 @@ void TLocalCoherenceLanczos::setup(void) { makeCoarseGrid(); } + LOG(Message) << "Coarse grid: " << coarseGrid_->GlobalDimensions() << std::endl; envCreate(FinePack, fineName_, Ls_, par().fineParams.Nm, env().getRbGrid(Ls_)); envCreate(CoarsePack, coarseName_, Ls_, par().coarseParams.Nm, coarseGridRb_.get()); auto &fine = envGet(FinePack, fineName_); diff --git a/lib/algorithms/iterative/LocalCoherenceLanczos.h b/lib/algorithms/iterative/LocalCoherenceLanczos.h index 3f997a01..0c7f4f65 100644 --- a/lib/algorithms/iterative/LocalCoherenceLanczos.h +++ b/lib/algorithms/iterative/LocalCoherenceLanczos.h @@ -332,7 +332,7 @@ public: // create a smoother and see if we can get a cheap convergence test and smooth inside the IRL ////////////////////////////////////////////////////////////////////////////////////////////////// Chebyshev ChebySmooth(cheby_smooth); - ProjectedFunctionHermOp ChebyOp (ChebySmooth,_FineOp,_subspace); + ProjectedFunctionHermOp ChebyOp (ChebySmooth,_FineOp,subspace); ImplicitlyRestartedLanczosSmoothedTester ChebySmoothTester(ChebyOp,ChebySmooth,_FineOp,subspace,relax); for(int k=0;k Date: Fri, 2 Mar 2018 11:47:38 +0000 Subject: [PATCH 24/38] changes requested for the pull request --- .../Hadrons/Modules/MSource/SeqConserved.hpp | 44 +++++---- extras/Hadrons/modules.inc | 2 +- lib/qcd/action/fermion/FermionOperator.h | 2 +- .../fermion/ImprovedStaggeredFermion.cc | 2 +- .../action/fermion/ImprovedStaggeredFermion.h | 2 +- .../fermion/ImprovedStaggeredFermion5D.cc | 2 +- .../fermion/ImprovedStaggeredFermion5D.h | 2 +- lib/qcd/action/fermion/WilsonFermion.cc | 71 +------------- lib/qcd/action/fermion/WilsonFermion.h | 4 +- lib/qcd/action/fermion/WilsonFermion5D.cc | 94 +------------------ lib/qcd/action/fermion/WilsonFermion5D.h | 2 +- 11 files changed, 38 insertions(+), 189 deletions(-) diff --git a/extras/Hadrons/Modules/MSource/SeqConserved.hpp b/extras/Hadrons/Modules/MSource/SeqConserved.hpp index 0fb23cc4..5c4a320b 100644 --- a/extras/Hadrons/Modules/MSource/SeqConserved.hpp +++ b/extras/Hadrons/Modules/MSource/SeqConserved.hpp @@ -2,7 +2,7 @@ Grid physics library, www.github.com/paboyle/Grid -Source file: extras/Hadrons/Modules/MContraction/SeqConserved.hpp +Source file: extras/Hadrons/Modules/MSource/SeqConserved.hpp Copyright (C) 2015-2018 @@ -35,8 +35,6 @@ See the full license in the file "LICENSE" in the top level distribution directo #include #include -#include - BEGIN_HADRONS_NAMESPACE /* @@ -100,6 +98,9 @@ protected: virtual void setup(void); // execution virtual void execute(void); +private: + bool SeqhasPhase_{false}; + std::string SeqmomphName_; }; MODULE_REGISTER_NS(SeqConserved, TSeqConserved, MSource); @@ -112,6 +113,7 @@ MODULE_REGISTER_NS(SeqConserved, TSeqConserved, MSource); template TSeqConserved::TSeqConserved(const std::string name) : Module(name) +, SeqmomphName_ (name + "_Seqmomph") {} // dependencies/products /////////////////////////////////////////////////////// @@ -139,7 +141,7 @@ void TSeqConserved::setup(void) auto Ls_ = env().getObjectLs(par().action); envCreateLat(PropagatorField, getName(), Ls_); envTmpLat(PropagatorField, "src_tmp"); - envTmpLat(LatticeComplex, "mom_phase"); + envCacheLat(LatticeComplex, SeqmomphName_); envTmpLat(LatticeComplex, "coor"); envTmpLat(LatticeComplex, "latt_compl"); } @@ -152,37 +154,45 @@ void TSeqConserved::execute(void) { LOG(Message) << "Generating sequential source with conserved " << par().curr_type << " current at " - << "t = " << par().tA << std::endl; + << "t = " << par().tA << " summed over the indices " + << par().mu_min << " <= mu <= " << par().mu_max + << std::endl; } else { LOG(Message) << "Generating sequential source with conserved " << par().curr_type << " current for " << par().tA << " <= t <= " - << par().tB << std::endl; + << par().tB << " summed over the indices " + << par().mu_min << " <= mu <= " << par().mu_max + << std::endl; } auto &src = envGet(PropagatorField, getName()); envGetTmp(PropagatorField, src_tmp); src_tmp = src; auto &q = envGet(PropagatorField, par().q); auto &mat = envGet(FMat, par().action); - envGetTmp(LatticeComplex, mom_phase); - envGetTmp(LatticeComplex, coor); envGetTmp(LatticeComplex, latt_compl); src = zero; //exp(ipx) - std::vector mom = strToVec(par().mom); - mom_phase = zero; - Complex i(0.0,1.0); - for(unsigned int mu = 0; mu < env().getNd(); mu++) - { - LatticeCoordinate(coor, mu); - mom_phase = mom_phase + (mom[mu]/env().getGrid()->_fdimensions[mu])*coor; + auto &mom_phase = envGet(LatticeComplex, SeqmomphName_); + if (!SeqhasPhase_) + { + std::vector mom = strToVec(par().mom); + mom_phase = zero; + Complex i(0.0,1.0); + envGetTmp(LatticeComplex, coor); + for(unsigned int mu = 0; mu < env().getNd(); mu++) + { + LatticeCoordinate(coor, mu); + mom_phase = mom_phase + (mom[mu]/env().getGrid()->_fdimensions[mu])*coor; + } + mom_phase = exp((Real)(2*M_PI)*i*mom_phase); + SeqhasPhase_ = true; } - mom_phase = exp((Real)(2*M_PI)*i*mom_phase); - LOG(Message) << "Inserting momentum " << mom << std::endl; + LOG(Message) << "Inserting momentum " << strToVec(par().mom) << std::endl; diff --git a/extras/Hadrons/modules.inc b/extras/Hadrons/modules.inc index 8bf22cc0..90602275 100644 --- a/extras/Hadrons/modules.inc +++ b/extras/Hadrons/modules.inc @@ -49,5 +49,5 @@ modules_hpp =\ Modules/MScalarSUN/TwoPoint.hpp \ Modules/MScalarSUN/TrPhi.hpp \ Modules/MIO/LoadNersc.hpp \ -Modules/MIO/LoadBinary.hpp + Modules/MIO/LoadBinary.hpp diff --git a/lib/qcd/action/fermion/FermionOperator.h b/lib/qcd/action/fermion/FermionOperator.h index 94e065cc..9eb74baa 100644 --- a/lib/qcd/action/fermion/FermionOperator.h +++ b/lib/qcd/action/fermion/FermionOperator.h @@ -127,7 +127,7 @@ namespace Grid { unsigned int mu, unsigned int tmin, unsigned int tmax, - Lattice> &lattice_cmplx)=0; + ComplexField &lattice_cmplx)=0; }; } diff --git a/lib/qcd/action/fermion/ImprovedStaggeredFermion.cc b/lib/qcd/action/fermion/ImprovedStaggeredFermion.cc index 3a296e52..ed07fb5b 100644 --- a/lib/qcd/action/fermion/ImprovedStaggeredFermion.cc +++ b/lib/qcd/action/fermion/ImprovedStaggeredFermion.cc @@ -413,7 +413,7 @@ void ImprovedStaggeredFermion::SeqConservedCurrent(PropagatorField &q_in, unsigned int mu, unsigned int tmin, unsigned int tmax, - Lattice> &lattice_cmplx) + ComplexField &lattice_cmplx) { assert(0); diff --git a/lib/qcd/action/fermion/ImprovedStaggeredFermion.h b/lib/qcd/action/fermion/ImprovedStaggeredFermion.h index 015eb6bb..d2426ef7 100644 --- a/lib/qcd/action/fermion/ImprovedStaggeredFermion.h +++ b/lib/qcd/action/fermion/ImprovedStaggeredFermion.h @@ -172,7 +172,7 @@ class ImprovedStaggeredFermion : public StaggeredKernels, public ImprovedS unsigned int mu, unsigned int tmin, unsigned int tmax, - Lattice> &lattice_cmplx); + ComplexField &lattice_cmplx); }; typedef ImprovedStaggeredFermion ImprovedStaggeredFermionF; diff --git a/lib/qcd/action/fermion/ImprovedStaggeredFermion5D.cc b/lib/qcd/action/fermion/ImprovedStaggeredFermion5D.cc index 3640a26b..3510b8f1 100644 --- a/lib/qcd/action/fermion/ImprovedStaggeredFermion5D.cc +++ b/lib/qcd/action/fermion/ImprovedStaggeredFermion5D.cc @@ -425,7 +425,7 @@ void ImprovedStaggeredFermion5D::SeqConservedCurrent(PropagatorField &q_in unsigned int mu, unsigned int tmin, unsigned int tmax, - Lattice> &lattice_cmplx) + ComplexField &lattice_cmplx) { assert(0); diff --git a/lib/qcd/action/fermion/ImprovedStaggeredFermion5D.h b/lib/qcd/action/fermion/ImprovedStaggeredFermion5D.h index 1c211938..e9cf2013 100644 --- a/lib/qcd/action/fermion/ImprovedStaggeredFermion5D.h +++ b/lib/qcd/action/fermion/ImprovedStaggeredFermion5D.h @@ -184,7 +184,7 @@ namespace QCD { unsigned int mu, unsigned int tmin, unsigned int tmax, - Lattice> &lattice_cmplx); + ComplexField &lattice_cmplx); }; }} diff --git a/lib/qcd/action/fermion/WilsonFermion.cc b/lib/qcd/action/fermion/WilsonFermion.cc index 5b14651b..2d9cf22d 100644 --- a/lib/qcd/action/fermion/WilsonFermion.cc +++ b/lib/qcd/action/fermion/WilsonFermion.cc @@ -407,75 +407,6 @@ void WilsonFermion::ContractConservedCurrent(PropagatorField &q_in_1, } } -//template -//void WilsonFermion::SeqConservedCurrent(PropagatorField &q_in, -// PropagatorField &q_out, -// Current curr_type, -// unsigned int mu, -// std::vector mom, -// unsigned int tmin, -// unsigned int tmax) -//{ -// conformable(_grid, q_in._grid); -// conformable(_grid, q_out._grid); -// Lattice> ph(_grid), coor(_grid); -// Complex i(0.0,1.0); -// PropagatorField tmpFwd(_grid), tmpBwd(_grid), tmp(_grid); -// unsigned int tshift = (mu == Tp) ? 1 : 0; -// unsigned int LLt = GridDefaultLatt()[Tp]; -// -// // Momentum projection -// ph = zero; -// for(unsigned int mu = 0; mu < Nd - 1; mu++) -// { -// LatticeCoordinate(coor, mu); -// ph = ph + mom[mu]*coor*((1./(_grid->_fdimensions[mu]))); -// } -// ph = exp((Real)(2*M_PI)*i*ph); -// -// q_out = zero; -// LatticeInteger coords(_grid); -// LatticeCoordinate(coords, Tp); -// -// // Need q(x + mu) and q(x - mu). -// tmp = Cshift(q_in, mu, 1); -// tmpFwd = tmp*ph; -// tmp = ph*q_in; -// tmpBwd = Cshift(tmp, mu, -1); -// -// parallel_for (unsigned int sU = 0; sU < Umu._grid->oSites(); ++sU) -// { -// // Compute the sequential conserved current insertion only if our simd -// // object contains a timeslice we need. -// vInteger t_mask = ((coords._odata[sU] >= tmin) && -// (coords._odata[sU] <= tmax)); -// Integer timeSlices = Reduce(t_mask); -// -// if (timeSlices > 0) -// { -// Kernels::SeqConservedCurrentSiteFwd(tmpFwd._odata[sU], -// q_out._odata[sU], -// Umu, sU, mu, t_mask); -// } -// -// // Repeat for backward direction. -// t_mask = ((coords._odata[sU] >= (tmin + tshift)) && -// (coords._odata[sU] <= (tmax + tshift))); -// -// //if tmax = LLt-1 (last timeslice) include timeslice 0 if the time is shifted (mu=3) -// unsigned int t0 = 0; -// if((tmax==LLt-1) && (tshift==1)) t_mask = (t_mask || (coords._odata[sU] == t0 )); -// -// timeSlices = Reduce(t_mask); -// -// if (timeSlices > 0) -// { -// Kernels::SeqConservedCurrentSiteBwd(tmpBwd._odata[sU], -// q_out._odata[sU], -// Umu, sU, mu, t_mask); -// } -// } -//} template void WilsonFermion::SeqConservedCurrent(PropagatorField &q_in, @@ -484,7 +415,7 @@ void WilsonFermion::SeqConservedCurrent(PropagatorField &q_in, unsigned int mu, unsigned int tmin, unsigned int tmax, - Lattice> &lattice_cmplx) + ComplexField &lattice_cmplx) { conformable(_grid, q_in._grid); conformable(_grid, q_out._grid); diff --git a/lib/qcd/action/fermion/WilsonFermion.h b/lib/qcd/action/fermion/WilsonFermion.h index 1ef6451e..ea25ed7f 100644 --- a/lib/qcd/action/fermion/WilsonFermion.h +++ b/lib/qcd/action/fermion/WilsonFermion.h @@ -174,13 +174,13 @@ class WilsonFermion : public WilsonKernels, public WilsonFermionStatic { PropagatorField &q_out, Current curr_type, unsigned int mu); - void SeqConservedCurrent(PropagatorField &q_in, + void SeqConservedCurrent(PropagatorField &q_in, PropagatorField &q_out, Current curr_type, unsigned int mu, unsigned int tmin, unsigned int tmax, - Lattice> &lattice_cmplx); + ComplexField &lattice_cmplx); }; typedef WilsonFermion WilsonFermionF; diff --git a/lib/qcd/action/fermion/WilsonFermion5D.cc b/lib/qcd/action/fermion/WilsonFermion5D.cc index 8a795a1e..6f82aad2 100644 --- a/lib/qcd/action/fermion/WilsonFermion5D.cc +++ b/lib/qcd/action/fermion/WilsonFermion5D.cc @@ -779,98 +779,6 @@ void WilsonFermion5D::ContractConservedCurrent(PropagatorField &q_in_1, } -//template -//void WilsonFermion5D::SeqConservedCurrent(PropagatorField &q_in, -// PropagatorField &q_out, -// Current curr_type, -// unsigned int mu, -// std::vector mom, -// unsigned int tmin, -// unsigned int tmax) -//{ -// conformable(q_in._grid, FermionGrid()); -// conformable(q_in._grid, q_out._grid); -// Lattice> ph(FermionGrid()), coor(FermionGrid()); -// PropagatorField tmpFwd(FermionGrid()), tmpBwd(FermionGrid()), -// tmp(FermionGrid()); -// Complex i(0.0, 1.0); -// unsigned int tshift = (mu == Tp) ? 1 : 0; -// unsigned int LLs = q_in._grid->_rdimensions[0]; -// unsigned int LLt = GridDefaultLatt()[Tp]; -// -// // Momentum projection. -// ph = zero; -// for(unsigned int nu = 0; nu < Nd - 1; nu++) -// { -// // Shift coordinate lattice index by 1 to account for 5th dimension. -// LatticeCoordinate(coor, nu + 1); -// ph = ph + mom[nu]*coor*((1./(_FourDimGrid->_fdimensions[nu]))); -// } -// ph = exp((Real)(2*M_PI)*i*ph); -// -// q_out = zero; -// LatticeInteger coords(_FourDimGrid); -// LatticeCoordinate(coords, Tp); -// -// -// // Need q(x + mu, s) and q(x - mu, s). 5D lattice so shift 4D coordinate mu -// // by one. -// tmp = Cshift(q_in, mu + 1, 1); -// tmpFwd = tmp*ph; -// tmp = ph*q_in; -// tmpBwd = Cshift(tmp, mu + 1, -1); -// -// parallel_for (unsigned int sU = 0; sU < Umu._grid->oSites(); ++sU) -// { -// // Compute the sequential conserved current insertion only if our simd -// // object contains a timeslice we need. -// vInteger t_mask = ((coords._odata[sU] >= tmin) && -// (coords._odata[sU] <= tmax)); -// Integer timeSlices = Reduce(t_mask); -// -// if (timeSlices > 0) -// { -// unsigned int sF = sU * LLs; -// for (unsigned int s = 0; s < LLs; ++s) -// { -// bool axial_sign = ((curr_type == Current::Axial) && (s < (LLs / 2))); -// bool tadpole_sign = (curr_type == Current::Tadpole); -// bool switch_sgn = tadpole_sign || axial_sign; -// -// Kernels::SeqConservedCurrentSiteFwd(tmpFwd._odata[sF], -// q_out._odata[sF], Umu, sU, -// mu, t_mask, switch_sgn); -// ++sF; -// } -// } -// -// // Repeat for backward direction. -// t_mask = ((coords._odata[sU] >= (tmin + tshift)) && -// (coords._odata[sU] <= (tmax + tshift))); -// -// //if tmax = LLt-1 (last timeslice) include timeslice 0 if the time is shifted (mu=3) -// unsigned int t0 = 0; -// if((tmax==LLt-1) && (tshift==1)) t_mask = (t_mask || (coords._odata[sU] == t0 )); -// -// timeSlices = Reduce(t_mask); -// -// if (timeSlices > 0) -// { -// unsigned int sF = sU * LLs; -// for (unsigned int s = 0; s < LLs; ++s) -// { -// bool axial_sign = ((curr_type == Current::Axial) && (s < (LLs / 2))); -// Kernels::SeqConservedCurrentSiteBwd(tmpBwd._odata[sF], -// q_out._odata[sF], Umu, sU, -// mu, t_mask, axial_sign); -// ++sF; -// } -// } -// } -//} - - - template void WilsonFermion5D::SeqConservedCurrent(PropagatorField &q_in, @@ -879,7 +787,7 @@ void WilsonFermion5D::SeqConservedCurrent(PropagatorField &q_in, unsigned int mu, unsigned int tmin, unsigned int tmax, - Lattice> &lattice_cmplx) + ComplexField &lattice_cmplx) { conformable(q_in._grid, FermionGrid()); conformable(q_in._grid, q_out._grid); diff --git a/lib/qcd/action/fermion/WilsonFermion5D.h b/lib/qcd/action/fermion/WilsonFermion5D.h index 6a7c6c7d..21da4c31 100644 --- a/lib/qcd/action/fermion/WilsonFermion5D.h +++ b/lib/qcd/action/fermion/WilsonFermion5D.h @@ -228,7 +228,7 @@ namespace QCD { unsigned int mu, unsigned int tmin, unsigned int tmax, - Lattice> &lattice_cmplx); + ComplexField &lattice_cmplx); }; }} From 37fe9442248f6a894220e6f0d0314fc247d69b25 Mon Sep 17 00:00:00 2001 From: Antonin Portelli Date: Fri, 2 Mar 2018 14:14:11 +0000 Subject: [PATCH 25/38] Hadrons: scalar kinetic term --- extras/Hadrons/Modules.hpp | 31 +-- .../Hadrons/Modules/MScalarSUN/TrKinetic.hpp | 198 ++++++++++++++++++ extras/Hadrons/modules.inc | 1 + 3 files changed, 200 insertions(+), 30 deletions(-) create mode 100644 extras/Hadrons/Modules/MScalarSUN/TrKinetic.hpp diff --git a/extras/Hadrons/Modules.hpp b/extras/Hadrons/Modules.hpp index 57fa93a3..274a6d86 100644 --- a/extras/Hadrons/Modules.hpp +++ b/extras/Hadrons/Modules.hpp @@ -1,33 +1,3 @@ -/************************************************************************************* - -Grid physics library, www.github.com/paboyle/Grid - -Source file: extras/Hadrons/Modules.hpp - -Copyright (C) 2015-2018 - -Author: Antonin Portelli -Author: Guido Cossu -Author: Lanny91 -Author: pretidav - -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 @@ -65,5 +35,6 @@ See the full license in the file "LICENSE" in the top level distribution directo #include #include #include +#include #include #include diff --git a/extras/Hadrons/Modules/MScalarSUN/TrKinetic.hpp b/extras/Hadrons/Modules/MScalarSUN/TrKinetic.hpp new file mode 100644 index 00000000..83e0c29e --- /dev/null +++ b/extras/Hadrons/Modules/MScalarSUN/TrKinetic.hpp @@ -0,0 +1,198 @@ +/************************************************************************************* + +Grid physics library, www.github.com/paboyle/Grid + +Source file: extras/Hadrons/Modules/MScalarSUN/TrKinetic.hpp + +Copyright (C) 2015-2018 + + +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_MScalarSUN_TrKinetic_hpp_ +#define Hadrons_MScalarSUN_TrKinetic_hpp_ + +#include +#include +#include + +BEGIN_HADRONS_NAMESPACE + +/****************************************************************************** + * TrKinetic * + ******************************************************************************/ +BEGIN_MODULE_NAMESPACE(MScalarSUN) + +class TrKineticPar: Serializable +{ +public: + GRID_SERIALIZABLE_ENUM(DiffType, undef, forward, 1, backward, 2, central, 3); + GRID_SERIALIZABLE_CLASS_MEMBERS(TrKineticPar, + std::string, field, + DiffType, type, + std::string, output); +}; + +template +class TTrKinetic: public Module +{ +public: + typedef typename SImpl::Field Field; + typedef typename SImpl::ComplexField ComplexField; + class Result: Serializable + { + public: + GRID_SERIALIZABLE_CLASS_MEMBERS(Result, + std::string, op, + Complex , value); + }; +public: + // constructor + TTrKinetic(const std::string name); + // destructor + virtual ~TTrKinetic(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 outName(const unsigned int mu, const unsigned int nu); + std::string bufName(const unsigned int mu); +}; + +MODULE_REGISTER_NS(TrKineticSU2, TTrKinetic>, MScalarSUN); +MODULE_REGISTER_NS(TrKineticSU3, TTrKinetic>, MScalarSUN); +MODULE_REGISTER_NS(TrKineticSU4, TTrKinetic>, MScalarSUN); +MODULE_REGISTER_NS(TrKineticSU5, TTrKinetic>, MScalarSUN); +MODULE_REGISTER_NS(TrKineticSU6, TTrKinetic>, MScalarSUN); + +/****************************************************************************** + * TTrKinetic implementation * + ******************************************************************************/ +// constructor ///////////////////////////////////////////////////////////////// +template +TTrKinetic::TTrKinetic(const std::string name) +: Module(name) +{} + +// dependencies/products /////////////////////////////////////////////////////// +template +std::vector TTrKinetic::getInput(void) +{ + std::vector in = {par().field}; + + return in; +} + +template +std::vector TTrKinetic::getOutput(void) +{ + std::vector out ; + + for (unsigned int mu = 0; mu < env().getNd(); ++mu) + for (unsigned int nu = mu; nu < env().getNd(); ++nu) + { + out.push_back(outName(mu, nu)); + } + + return out; +} + +// setup /////////////////////////////////////////////////////////////////////// +template +void TTrKinetic::setup(void) +{ + for (unsigned int mu = 0; mu < env().getNd(); ++mu) + for (unsigned int nu = mu; nu < env().getNd(); ++nu) + { + envCreateLat(ComplexField, outName(mu, nu)); + } + envTmp(std::vector, "der", 1, env().getNd(), env().getGrid()); +} + +// execution /////////////////////////////////////////////////////////////////// +template +void TTrKinetic::execute(void) +{ + LOG(Message) << "Computing tr(d_mu phi*d_nu phi) using " << par().type + << " derivative" << std::endl; + + std::vector result; + auto &phi = envGet(Field, par().field); + + envGetTmp(std::vector, der); + for (unsigned int mu = 0; mu < env().getNd(); ++mu) + { + switch(par().type) + { + case TrKineticPar::DiffType::backward: + der[mu] = phi - Cshift(phi, mu, -1); + break; + case TrKineticPar::DiffType::forward: + der[mu] = Cshift(phi, mu, 1) - phi; + break; + case TrKineticPar::DiffType::central: + der[mu] = 0.5*(Cshift(phi, mu, 1) - Cshift(phi, mu, -1)); + break; + } + } + for (unsigned int mu = 0; mu < env().getNd(); ++mu) + for (unsigned int nu = mu; nu < env().getNd(); ++nu) + { + auto &out = envGet(ComplexField, outName(mu, nu)); + + out = -trace(der[mu]*der[nu]); + if (!par().output.empty()) + { + Result r; + + r.op = "tr(d_" + std::to_string(mu) + "phi*d_" + + std::to_string(nu) + "phi)"; + r.value = TensorRemove(sum(out)); + result.push_back(r); + } + } + if (result.size() > 0) + { + saveResult(par().output, "trkinetic", result); + } +} + +// variable name generators //////////////////////////////////////////////////// +template +std::string TTrKinetic::outName(const unsigned int mu, + const unsigned int nu) +{ + return getName() + "_" + std::to_string(mu) + "_" + std::to_string(nu); +} + +template +std::string TTrKinetic::bufName(const unsigned int mu) +{ + return "d_" + std::to_string(mu); +} + + +END_MODULE_NAMESPACE + +END_HADRONS_NAMESPACE + +#endif // Hadrons_MScalarSUN_TrKinetic_hpp_ diff --git a/extras/Hadrons/modules.inc b/extras/Hadrons/modules.inc index 79db9396..3ee3fb69 100644 --- a/extras/Hadrons/modules.inc +++ b/extras/Hadrons/modules.inc @@ -48,6 +48,7 @@ modules_hpp =\ Modules/MScalarSUN/TrMag.hpp \ Modules/MScalarSUN/TwoPoint.hpp \ Modules/MScalarSUN/TrPhi.hpp \ + Modules/MScalarSUN/TrKinetic.hpp \ Modules/MIO/LoadNersc.hpp \ Modules/MIO/LoadBinary.hpp From c0a929aef7e72bcc68e8a355c758d34cbe94f90c Mon Sep 17 00:00:00 2001 From: Antonin Portelli Date: Fri, 2 Mar 2018 14:29:54 +0000 Subject: [PATCH 26/38] Hadrons: code cleaning --- extras/Hadrons/LanczosUtils.hpp | 2 +- extras/Hadrons/Modules.hpp | 30 +++++++++++++++++++ extras/Hadrons/Modules/MScalarSUN/Div.hpp | 4 +-- .../Hadrons/Modules/MScalarSUN/TrKinetic.hpp | 1 + extras/Hadrons/Modules/MScalarSUN/TrMag.hpp | 4 +-- extras/Hadrons/Modules/MScalarSUN/TrPhi.hpp | 4 +-- .../Hadrons/Modules/MScalarSUN/TwoPoint.hpp | 2 +- .../Modules/MSolver/LocalCoherenceLanczos.hpp | 6 ++-- 8 files changed, 42 insertions(+), 11 deletions(-) diff --git a/extras/Hadrons/LanczosUtils.hpp b/extras/Hadrons/LanczosUtils.hpp index f8e7693a..a080da4b 100644 --- a/extras/Hadrons/LanczosUtils.hpp +++ b/extras/Hadrons/LanczosUtils.hpp @@ -112,4 +112,4 @@ using CoarseEigenPack = EigenPack< END_HADRONS_NAMESPACE -#endif // Hadrons_LanczosUtils_hpp_ \ No newline at end of file +#endif // Hadrons_LanczosUtils_hpp_ diff --git a/extras/Hadrons/Modules.hpp b/extras/Hadrons/Modules.hpp index 274a6d86..dc0d2420 100644 --- a/extras/Hadrons/Modules.hpp +++ b/extras/Hadrons/Modules.hpp @@ -1,3 +1,33 @@ +/************************************************************************************* + +Grid physics library, www.github.com/paboyle/Grid + +Source file: extras/Hadrons/Modules.hpp + +Copyright (C) 2015-2018 + +Author: Antonin Portelli +Author: Guido Cossu +Author: Lanny91 +Author: pretidav + +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 diff --git a/extras/Hadrons/Modules/MScalarSUN/Div.hpp b/extras/Hadrons/Modules/MScalarSUN/Div.hpp index cab6e88c..7d5b08b2 100644 --- a/extras/Hadrons/Modules/MScalarSUN/Div.hpp +++ b/extras/Hadrons/Modules/MScalarSUN/Div.hpp @@ -35,7 +35,7 @@ See the full license in the file "LICENSE" in the top level distribution directo BEGIN_HADRONS_NAMESPACE /****************************************************************************** - * Div * + * Divergence of a vector field * ******************************************************************************/ BEGIN_MODULE_NAMESPACE(MScalarSUN) @@ -83,7 +83,7 @@ MODULE_REGISTER_NS(DivSU5, TDiv>, MScalarSUN); MODULE_REGISTER_NS(DivSU6, TDiv>, MScalarSUN); /****************************************************************************** - * TDiv implementation * + * TDiv implementation * ******************************************************************************/ // constructor ///////////////////////////////////////////////////////////////// template diff --git a/extras/Hadrons/Modules/MScalarSUN/TrKinetic.hpp b/extras/Hadrons/Modules/MScalarSUN/TrKinetic.hpp index 83e0c29e..75318488 100644 --- a/extras/Hadrons/Modules/MScalarSUN/TrKinetic.hpp +++ b/extras/Hadrons/Modules/MScalarSUN/TrKinetic.hpp @@ -6,6 +6,7 @@ Source file: extras/Hadrons/Modules/MScalarSUN/TrKinetic.hpp Copyright (C) 2015-2018 +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 diff --git a/extras/Hadrons/Modules/MScalarSUN/TrMag.hpp b/extras/Hadrons/Modules/MScalarSUN/TrMag.hpp index 2858984c..d7d1cce7 100644 --- a/extras/Hadrons/Modules/MScalarSUN/TrMag.hpp +++ b/extras/Hadrons/Modules/MScalarSUN/TrMag.hpp @@ -35,7 +35,7 @@ See the full license in the file "LICENSE" in the top level distribution directo BEGIN_HADRONS_NAMESPACE /****************************************************************************** - * Module to compute tr(mag^n) * + * Trace of powers of the magnetisation * ******************************************************************************/ BEGIN_MODULE_NAMESPACE(MScalarSUN) @@ -117,7 +117,7 @@ template void TTrMag::execute(void) { LOG(Message) << "Computing tr(mag^n) for n even up to " << par().maxPow - << "..." << std::endl; + << std::endl; std::vector result; auto &phi = envGet(Field, par().field); diff --git a/extras/Hadrons/Modules/MScalarSUN/TrPhi.hpp b/extras/Hadrons/Modules/MScalarSUN/TrPhi.hpp index 4b305128..e5e2b0ef 100644 --- a/extras/Hadrons/Modules/MScalarSUN/TrPhi.hpp +++ b/extras/Hadrons/Modules/MScalarSUN/TrPhi.hpp @@ -35,7 +35,7 @@ See the full license in the file "LICENSE" in the top level distribution directo BEGIN_HADRONS_NAMESPACE /****************************************************************************** - * Module to compute tr(phi^n) * + * Trace of powers of a scalar field * ******************************************************************************/ BEGIN_MODULE_NAMESPACE(MScalarSUN) @@ -136,7 +136,7 @@ template void TTrPhi::execute(void) { LOG(Message) << "Computing tr(phi^n) for n even up to " << par().maxPow - << "..." << std::endl; + << std::endl; std::vector result; auto &phi = envGet(Field, par().field); diff --git a/extras/Hadrons/Modules/MScalarSUN/TwoPoint.hpp b/extras/Hadrons/Modules/MScalarSUN/TwoPoint.hpp index 58570dd5..8772cfbe 100644 --- a/extras/Hadrons/Modules/MScalarSUN/TwoPoint.hpp +++ b/extras/Hadrons/Modules/MScalarSUN/TwoPoint.hpp @@ -87,7 +87,7 @@ MODULE_REGISTER_NS(TwoPointSU5, TTwoPoint>, MScalarSUN); MODULE_REGISTER_NS(TwoPointSU6, TTwoPoint>, MScalarSUN); /****************************************************************************** - * TTwoPoint implementation * + * TTwoPoint implementation * ******************************************************************************/ // constructor ///////////////////////////////////////////////////////////////// template diff --git a/extras/Hadrons/Modules/MSolver/LocalCoherenceLanczos.hpp b/extras/Hadrons/Modules/MSolver/LocalCoherenceLanczos.hpp index 284b233b..6e2103ce 100644 --- a/extras/Hadrons/Modules/MSolver/LocalCoherenceLanczos.hpp +++ b/extras/Hadrons/Modules/MSolver/LocalCoherenceLanczos.hpp @@ -36,8 +36,8 @@ See the full license in the file "LICENSE" in the top level distribution directo BEGIN_HADRONS_NAMESPACE /****************************************************************************** - * LocalCoherenceLanczos * - ******************************************************************************/ + * Local coherence Lanczos eigensolver * + *****************************************************************************/ BEGIN_MODULE_NAMESPACE(MSolver) class LocalCoherenceLanczosPar: Serializable @@ -98,7 +98,7 @@ MODULE_REGISTER_NS(ZLocalCoherenceLanczos, MSolver); /****************************************************************************** - * TLocalCoherenceLanczos implementation * + * TLocalCoherenceLanczos implementation * ******************************************************************************/ // constructor ///////////////////////////////////////////////////////////////// template From 550142bd6a2a4326b50871da51f362bc56fe10f7 Mon Sep 17 00:00:00 2001 From: Antonin Portelli Date: Fri, 2 Mar 2018 14:30:45 +0000 Subject: [PATCH 27/38] Hadrons: more code cleaning --- extras/Hadrons/Modules/MScalarSUN/TrKinetic.hpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/extras/Hadrons/Modules/MScalarSUN/TrKinetic.hpp b/extras/Hadrons/Modules/MScalarSUN/TrKinetic.hpp index 75318488..19776bef 100644 --- a/extras/Hadrons/Modules/MScalarSUN/TrKinetic.hpp +++ b/extras/Hadrons/Modules/MScalarSUN/TrKinetic.hpp @@ -35,7 +35,7 @@ See the full license in the file "LICENSE" in the top level distribution directo BEGIN_HADRONS_NAMESPACE /****************************************************************************** - * TrKinetic * + * Trace of kinetic term * ******************************************************************************/ BEGIN_MODULE_NAMESPACE(MScalarSUN) @@ -86,7 +86,7 @@ MODULE_REGISTER_NS(TrKineticSU5, TTrKinetic>, MScalarSUN); MODULE_REGISTER_NS(TrKineticSU6, TTrKinetic>, MScalarSUN); /****************************************************************************** - * TTrKinetic implementation * + * TTrKinetic implementation * ******************************************************************************/ // constructor ///////////////////////////////////////////////////////////////// template From 614a0e8277f23c713ae6898b4744cc9c411d907f Mon Sep 17 00:00:00 2001 From: Antonin Portelli Date: Fri, 2 Mar 2018 17:34:23 +0000 Subject: [PATCH 28/38] Hadrons: Scalar SU(N) utility functions --- extras/Hadrons/Modules.hpp | 1 + extras/Hadrons/Modules/MScalarSUN/Div.hpp | 19 +---- .../Hadrons/Modules/MScalarSUN/TrKinetic.hpp | 39 ++------- extras/Hadrons/Modules/MScalarSUN/TrMag.hpp | 1 + extras/Hadrons/Modules/MScalarSUN/TrPhi.hpp | 17 +--- .../Hadrons/Modules/MScalarSUN/TwoPoint.hpp | 1 + extras/Hadrons/Modules/MScalarSUN/Utils.hpp | 80 +++++++++++++++++++ extras/Hadrons/modules.inc | 1 + 8 files changed, 97 insertions(+), 62 deletions(-) create mode 100644 extras/Hadrons/Modules/MScalarSUN/Utils.hpp diff --git a/extras/Hadrons/Modules.hpp b/extras/Hadrons/Modules.hpp index dc0d2420..68b96822 100644 --- a/extras/Hadrons/Modules.hpp +++ b/extras/Hadrons/Modules.hpp @@ -65,6 +65,7 @@ See the full license in the file "LICENSE" in the top level distribution directo #include #include #include +#include #include #include #include diff --git a/extras/Hadrons/Modules/MScalarSUN/Div.hpp b/extras/Hadrons/Modules/MScalarSUN/Div.hpp index 7d5b08b2..0ecb1dbe 100644 --- a/extras/Hadrons/Modules/MScalarSUN/Div.hpp +++ b/extras/Hadrons/Modules/MScalarSUN/Div.hpp @@ -31,6 +31,7 @@ See the full license in the file "LICENSE" in the top level distribution directo #include #include #include +#include BEGIN_HADRONS_NAMESPACE @@ -42,7 +43,6 @@ BEGIN_MODULE_NAMESPACE(MScalarSUN) class DivPar: Serializable { public: - GRID_SERIALIZABLE_ENUM(DiffType, undef, forward, 1, backward, 2, central, 3); GRID_SERIALIZABLE_CLASS_MEMBERS(DivPar, std::vector, op, DiffType, type, @@ -59,8 +59,8 @@ public: { public: GRID_SERIALIZABLE_CLASS_MEMBERS(Result, - DivPar::DiffType, type, - Complex, value); + DiffType, type, + Complex, value); }; public: // constructor @@ -135,18 +135,7 @@ void TDiv::execute(void) for (unsigned int mu = 0; mu < nd; ++mu) { auto &op = envGet(ComplexField, par().op[mu]); - switch(par().type) - { - case DivPar::DiffType::backward: - div += op - Cshift(op, mu, -1); - break; - case DivPar::DiffType::forward: - div += Cshift(op, mu, 1) - op; - break; - case DivPar::DiffType::central: - div += 0.5*(Cshift(op, mu, 1) - Cshift(op, mu, -1)); - break; - } + dmuAcc(div, op, mu, par().type); } if (!par().output.empty()) { diff --git a/extras/Hadrons/Modules/MScalarSUN/TrKinetic.hpp b/extras/Hadrons/Modules/MScalarSUN/TrKinetic.hpp index 19776bef..f1b520de 100644 --- a/extras/Hadrons/Modules/MScalarSUN/TrKinetic.hpp +++ b/extras/Hadrons/Modules/MScalarSUN/TrKinetic.hpp @@ -31,6 +31,7 @@ See the full license in the file "LICENSE" in the top level distribution directo #include #include #include +#include BEGIN_HADRONS_NAMESPACE @@ -42,7 +43,6 @@ BEGIN_MODULE_NAMESPACE(MScalarSUN) class TrKineticPar: Serializable { public: - GRID_SERIALIZABLE_ENUM(DiffType, undef, forward, 1, backward, 2, central, 3); GRID_SERIALIZABLE_CLASS_MEMBERS(TrKineticPar, std::string, field, DiffType, type, @@ -74,9 +74,6 @@ public: virtual void setup(void); // execution virtual void execute(void); -private: - std::string outName(const unsigned int mu, const unsigned int nu); - std::string bufName(const unsigned int mu); }; MODULE_REGISTER_NS(TrKineticSU2, TTrKinetic>, MScalarSUN); @@ -111,7 +108,7 @@ std::vector TTrKinetic::getOutput(void) for (unsigned int mu = 0; mu < env().getNd(); ++mu) for (unsigned int nu = mu; nu < env().getNd(); ++nu) { - out.push_back(outName(mu, nu)); + out.push_back(varName(getName(), mu, nu)); } return out; @@ -124,7 +121,7 @@ void TTrKinetic::setup(void) for (unsigned int mu = 0; mu < env().getNd(); ++mu) for (unsigned int nu = mu; nu < env().getNd(); ++nu) { - envCreateLat(ComplexField, outName(mu, nu)); + envCreateLat(ComplexField, varName(getName(), mu, nu)); } envTmp(std::vector, "der", 1, env().getNd(), env().getGrid()); } @@ -142,23 +139,12 @@ void TTrKinetic::execute(void) envGetTmp(std::vector, der); for (unsigned int mu = 0; mu < env().getNd(); ++mu) { - switch(par().type) - { - case TrKineticPar::DiffType::backward: - der[mu] = phi - Cshift(phi, mu, -1); - break; - case TrKineticPar::DiffType::forward: - der[mu] = Cshift(phi, mu, 1) - phi; - break; - case TrKineticPar::DiffType::central: - der[mu] = 0.5*(Cshift(phi, mu, 1) - Cshift(phi, mu, -1)); - break; - } + dmu(der[mu], phi, mu, par().type); } for (unsigned int mu = 0; mu < env().getNd(); ++mu) for (unsigned int nu = mu; nu < env().getNd(); ++nu) { - auto &out = envGet(ComplexField, outName(mu, nu)); + auto &out = envGet(ComplexField, varName(getName(), mu, nu)); out = -trace(der[mu]*der[nu]); if (!par().output.empty()) @@ -177,21 +163,6 @@ void TTrKinetic::execute(void) } } -// variable name generators //////////////////////////////////////////////////// -template -std::string TTrKinetic::outName(const unsigned int mu, - const unsigned int nu) -{ - return getName() + "_" + std::to_string(mu) + "_" + std::to_string(nu); -} - -template -std::string TTrKinetic::bufName(const unsigned int mu) -{ - return "d_" + std::to_string(mu); -} - - END_MODULE_NAMESPACE END_HADRONS_NAMESPACE diff --git a/extras/Hadrons/Modules/MScalarSUN/TrMag.hpp b/extras/Hadrons/Modules/MScalarSUN/TrMag.hpp index d7d1cce7..879951d5 100644 --- a/extras/Hadrons/Modules/MScalarSUN/TrMag.hpp +++ b/extras/Hadrons/Modules/MScalarSUN/TrMag.hpp @@ -31,6 +31,7 @@ See the full license in the file "LICENSE" in the top level distribution directo #include #include #include +#include BEGIN_HADRONS_NAMESPACE diff --git a/extras/Hadrons/Modules/MScalarSUN/TrPhi.hpp b/extras/Hadrons/Modules/MScalarSUN/TrPhi.hpp index e5e2b0ef..aa98bfb3 100644 --- a/extras/Hadrons/Modules/MScalarSUN/TrPhi.hpp +++ b/extras/Hadrons/Modules/MScalarSUN/TrPhi.hpp @@ -31,6 +31,7 @@ See the full license in the file "LICENSE" in the top level distribution directo #include #include #include +#include BEGIN_HADRONS_NAMESPACE @@ -73,9 +74,6 @@ public: virtual void setup(void); // execution virtual void execute(void); -private: - // output name generator - std::string outName(const unsigned int n); }; MODULE_REGISTER_NS(TrPhiSU2, TTrPhi>, MScalarSUN); @@ -109,7 +107,7 @@ std::vector TTrPhi::getOutput(void) for (unsigned int n = 2; n <= par().maxPow; n += 2) { - out.push_back(outName(n)); + out.push_back(varName(getName(), n)); } return out; @@ -127,7 +125,7 @@ void TTrPhi::setup(void) envTmpLat(Field, "buf"); for (unsigned int n = 2; n <= par().maxPow; n += 2) { - envCreateLat(ComplexField, outName(n)); + envCreateLat(ComplexField, varName(getName(), n)); } } @@ -147,7 +145,7 @@ void TTrPhi::execute(void) phi2 = -phi*phi; for (unsigned int n = 2; n <= par().maxPow; n += 2) { - auto &phin = envGet(ComplexField, outName(n)); + auto &phin = envGet(ComplexField, varName(getName(), n)); buf = buf*phi2; phin = trace(buf); @@ -166,13 +164,6 @@ void TTrPhi::execute(void) } } -// output name generator /////////////////////////////////////////////////////// -template -std::string TTrPhi::outName(const unsigned int n) -{ - return getName() + "_" + std::to_string(n); -} - END_MODULE_NAMESPACE END_HADRONS_NAMESPACE diff --git a/extras/Hadrons/Modules/MScalarSUN/TwoPoint.hpp b/extras/Hadrons/Modules/MScalarSUN/TwoPoint.hpp index 8772cfbe..5227e95a 100644 --- a/extras/Hadrons/Modules/MScalarSUN/TwoPoint.hpp +++ b/extras/Hadrons/Modules/MScalarSUN/TwoPoint.hpp @@ -31,6 +31,7 @@ See the full license in the file "LICENSE" in the top level distribution directo #include #include #include +#include BEGIN_HADRONS_NAMESPACE diff --git a/extras/Hadrons/Modules/MScalarSUN/Utils.hpp b/extras/Hadrons/Modules/MScalarSUN/Utils.hpp new file mode 100644 index 00000000..1a13f8f5 --- /dev/null +++ b/extras/Hadrons/Modules/MScalarSUN/Utils.hpp @@ -0,0 +1,80 @@ +#ifndef Hadrons_MScalarSUN_Utils_hpp_ +#define Hadrons_MScalarSUN_Utils_hpp_ + +#include +#include + +BEGIN_HADRONS_NAMESPACE + +BEGIN_MODULE_NAMESPACE(MScalarSUN) + +GRID_SERIALIZABLE_ENUM(DiffType, undef, forward, 1, backward, 2, central, 3); + +template +inline void dmu(Field &out, const Field &in, const unsigned int mu, const DiffType type) +{ + auto & env = Environment::getInstance(); + + if (mu >= env.getNd()) + { + HADRON_ERROR(Range, "Derivative direction out of range"); + } + switch(type) + { + case DiffType::backward: + out = in - Cshift(in, mu, -1); + break; + case DiffType::forward: + out = Cshift(in, mu, 1) - in; + break; + case DiffType::central: + out = 0.5*(Cshift(in, mu, 1) - Cshift(in, mu, -1)); + break; + default: + HADRON_ERROR(Argument, "Derivative type invalid"); + break; + } +} + +template +inline void dmuAcc(Field &out, const Field &in, const unsigned int mu, const DiffType type) +{ + auto & env = Environment::getInstance(); + + if (mu >= env.getNd()) + { + HADRON_ERROR(Range, "Derivative direction out of range"); + } + switch(type) + { + case DiffType::backward: + out += in - Cshift(in, mu, -1); + break; + case DiffType::forward: + out += Cshift(in, mu, 1) - in; + break; + case DiffType::central: + out += 0.5*(Cshift(in, mu, 1) - Cshift(in, mu, -1)); + break; + default: + HADRON_ERROR(Argument, "Derivative type invalid"); + break; + } +} + +inline std::string varName(const std::string name, const unsigned int mu) +{ + return name + "_" + std::to_string(mu); +} + +inline std::string varName(const std::string name, const unsigned int mu, + const unsigned int nu) +{ + return name + "_" + std::to_string(mu) + "_" + std::to_string(nu); +} + +END_MODULE_NAMESPACE + +END_HADRONS_NAMESPACE + +#endif // Hadrons_MScalarSUN_Utils_hpp_ diff --git a/extras/Hadrons/modules.inc b/extras/Hadrons/modules.inc index 3ee3fb69..ea94179e 100644 --- a/extras/Hadrons/modules.inc +++ b/extras/Hadrons/modules.inc @@ -48,6 +48,7 @@ modules_hpp =\ Modules/MScalarSUN/TrMag.hpp \ Modules/MScalarSUN/TwoPoint.hpp \ Modules/MScalarSUN/TrPhi.hpp \ + Modules/MScalarSUN/Utils.hpp \ Modules/MScalarSUN/TrKinetic.hpp \ Modules/MIO/LoadNersc.hpp \ Modules/MIO/LoadBinary.hpp From d9c435e282e3bce3ea42027f96475ceb74cdc959 Mon Sep 17 00:00:00 2001 From: Antonin Portelli Date: Fri, 2 Mar 2018 17:35:12 +0000 Subject: [PATCH 29/38] Hadrons: Scalar SU(N) transverse projection module --- extras/Hadrons/Modules.hpp | 1 + .../Hadrons/Modules/MScalarSUN/TransProj.hpp | 158 ++++++++++++++++++ extras/Hadrons/modules.inc | 1 + 3 files changed, 160 insertions(+) create mode 100644 extras/Hadrons/Modules/MScalarSUN/TransProj.hpp diff --git a/extras/Hadrons/Modules.hpp b/extras/Hadrons/Modules.hpp index 68b96822..d92f2661 100644 --- a/extras/Hadrons/Modules.hpp +++ b/extras/Hadrons/Modules.hpp @@ -66,6 +66,7 @@ See the full license in the file "LICENSE" in the top level distribution directo #include #include #include +#include #include #include #include diff --git a/extras/Hadrons/Modules/MScalarSUN/TransProj.hpp b/extras/Hadrons/Modules/MScalarSUN/TransProj.hpp new file mode 100644 index 00000000..59b7097c --- /dev/null +++ b/extras/Hadrons/Modules/MScalarSUN/TransProj.hpp @@ -0,0 +1,158 @@ +#ifndef Hadrons_MScalarSUN_TransProj_hpp_ +#define Hadrons_MScalarSUN_TransProj_hpp_ + +#include +#include +#include +#include + +BEGIN_HADRONS_NAMESPACE + +/****************************************************************************** + * Transverse projection * + ******************************************************************************/ +BEGIN_MODULE_NAMESPACE(MScalarSUN) + +class TransProjPar: Serializable +{ +public: + GRID_SERIALIZABLE_CLASS_MEMBERS(TransProjPar, + std::string, op, + DiffType, type, + std::string, output); +}; + +template +class TTransProj: public Module +{ +public: + typedef typename SImpl::Field Field; + typedef typename SImpl::ComplexField ComplexField; + class Result: Serializable + { + public: + GRID_SERIALIZABLE_CLASS_MEMBERS(Result, + std::string, op, + Complex , value); + }; +public: + // constructor + TTransProj(const std::string name); + // destructor + virtual ~TTransProj(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(TransProjSU2, TTransProj>, MScalarSUN); +MODULE_REGISTER_NS(TransProjSU3, TTransProj>, MScalarSUN); +MODULE_REGISTER_NS(TransProjSU4, TTransProj>, MScalarSUN); +MODULE_REGISTER_NS(TransProjSU5, TTransProj>, MScalarSUN); +MODULE_REGISTER_NS(TransProjSU6, TTransProj>, MScalarSUN); + +/****************************************************************************** + * TTransProj implementation * + ******************************************************************************/ +// constructor ///////////////////////////////////////////////////////////////// +template +TTransProj::TTransProj(const std::string name) +: Module(name) +{} + +// dependencies/products /////////////////////////////////////////////////////// +template +std::vector TTransProj::getInput(void) +{ + std::vector in = {par().op}; + + return in; +} + +template +std::vector TTransProj::getOutput(void) +{ + std::vector out = {getName()}; + + for (unsigned int mu = 0; mu < env().getNd(); ++mu) + for (unsigned int nu = mu; nu < env().getNd(); ++nu) + { + out.push_back(varName(getName(), mu, nu)); + } + + return out; +} + +// setup /////////////////////////////////////////////////////////////////////// +template +void TTransProj::setup(void) +{ + for (unsigned int mu = 0; mu < env().getNd(); ++mu) + for (unsigned int nu = mu; nu < env().getNd(); ++nu) + { + envCreateLat(ComplexField, varName(getName(), mu, nu)); + } + envTmpLat(ComplexField, "buf1"); + envTmpLat(ComplexField, "buf2"); + envTmpLat(ComplexField, "lap"); +} + +// execution /////////////////////////////////////////////////////////////////// +template +void TTransProj::execute(void) +{ + LOG(Message) << "Computing (delta_mu,nu d^2 - d_mu*d_nu)*op using " + << par().type << " derivatives and op= '" << par().op + << "'" << std::endl; + + std::vector result; + auto &op = envGet(ComplexField, par().op); + + envGetTmp(ComplexField, buf1); + envGetTmp(ComplexField, buf2); + envGetTmp(ComplexField, lap); + lap = zero; + for (unsigned int mu = 0; mu < env().getNd(); ++mu) + { + dmu(buf1, op, mu, par().type); + dmu(buf2, buf1, mu, par().type); + lap += buf2; + } + for (unsigned int mu = 0; mu < env().getNd(); ++mu) + for (unsigned int nu = mu; nu < env().getNd(); ++nu) + { + auto &out = envGet(ComplexField, varName(getName(), mu, nu)); + dmu(buf1, op, mu, par().type); + dmu(buf2, buf1, nu, par().type); + out = -buf2; + if (mu == nu) + { + out += lap; + } + if (!par().output.empty()) + { + Result r; + + r.op = "(delta_" + std::to_string(mu) + "," + std::to_string(nu) + + " d^2 - d_" + std::to_string(mu) + "*d_" + + std::to_string(nu) + ")*op"; + r.value = TensorRemove(sum(out)); + result.push_back(r); + } + } + if (result.size() > 0) + { + saveResult(par().output, "transproj", result); + } +} + +END_MODULE_NAMESPACE + +END_HADRONS_NAMESPACE + +#endif // Hadrons_MScalarSUN_TransProj_hpp_ diff --git a/extras/Hadrons/modules.inc b/extras/Hadrons/modules.inc index ea94179e..dbfa359c 100644 --- a/extras/Hadrons/modules.inc +++ b/extras/Hadrons/modules.inc @@ -49,6 +49,7 @@ modules_hpp =\ Modules/MScalarSUN/TwoPoint.hpp \ Modules/MScalarSUN/TrPhi.hpp \ Modules/MScalarSUN/Utils.hpp \ + Modules/MScalarSUN/TransProj.hpp \ Modules/MScalarSUN/TrKinetic.hpp \ Modules/MIO/LoadNersc.hpp \ Modules/MIO/LoadBinary.hpp From 1d83521daadf983e6bb599982c2d5552eaadfc1d Mon Sep 17 00:00:00 2001 From: Antonin Portelli Date: Fri, 2 Mar 2018 17:55:18 +0000 Subject: [PATCH 30/38] Hadrons: scalar SU(N) EMT --- extras/Hadrons/Modules.hpp | 1 + extras/Hadrons/Modules/MScalarSUN/EMT.hpp | 150 ++++++++++++++++++++++ extras/Hadrons/modules.inc | 1 + 3 files changed, 152 insertions(+) create mode 100644 extras/Hadrons/Modules/MScalarSUN/EMT.hpp diff --git a/extras/Hadrons/Modules.hpp b/extras/Hadrons/Modules.hpp index d92f2661..695dfd02 100644 --- a/extras/Hadrons/Modules.hpp +++ b/extras/Hadrons/Modules.hpp @@ -63,6 +63,7 @@ See the full license in the file "LICENSE" in the top level distribution directo #include #include #include +#include #include #include #include diff --git a/extras/Hadrons/Modules/MScalarSUN/EMT.hpp b/extras/Hadrons/Modules/MScalarSUN/EMT.hpp new file mode 100644 index 00000000..7f1a612d --- /dev/null +++ b/extras/Hadrons/Modules/MScalarSUN/EMT.hpp @@ -0,0 +1,150 @@ +#ifndef Hadrons_MScalarSUN_EMT_hpp_ +#define Hadrons_MScalarSUN_EMT_hpp_ + +#include +#include +#include +#include + +BEGIN_HADRONS_NAMESPACE + +/****************************************************************************** + * Energy-momentum tensor * + ******************************************************************************/ +BEGIN_MODULE_NAMESPACE(MScalarSUN) + +class EMTPar: Serializable +{ +public: + GRID_SERIALIZABLE_CLASS_MEMBERS(EMTPar, + std::string, kinetic, + std::string, phiPow, + std::string, improvement, + double , m2, + double , lambda, + double , g, + double , xi, + std::string, output); +}; + +template +class TEMT: public Module +{ +public: + typedef typename SImpl::Field Field; + typedef typename SImpl::ComplexField ComplexField; +public: + // constructor + TEMT(const std::string name); + // destructor + virtual ~TEMT(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(EMTSU2, TEMT>, MScalarSUN); +MODULE_REGISTER_NS(EMTSU3, TEMT>, MScalarSUN); +MODULE_REGISTER_NS(EMTSU4, TEMT>, MScalarSUN); +MODULE_REGISTER_NS(EMTSU5, TEMT>, MScalarSUN); +MODULE_REGISTER_NS(EMTSU6, TEMT>, MScalarSUN); + +/****************************************************************************** + * TEMT implementation * + ******************************************************************************/ +// constructor ///////////////////////////////////////////////////////////////// +template +TEMT::TEMT(const std::string name) +: Module(name) +{} + +// dependencies/products /////////////////////////////////////////////////////// +template +std::vector TEMT::getInput(void) +{ + std::vector in; + + for (unsigned int mu = 0; mu < env().getNd(); ++mu) + for (unsigned int nu = mu; nu < env().getNd(); ++nu) + { + in.push_back(varName(par().kinetic, mu, nu)); + in.push_back(varName(par().improvement, mu, nu)); + } + in.push_back(varName(par().phiPow, 2)); + in.push_back(varName(par().phiPow, 4)); + + return in; +} + +template +std::vector TEMT::getOutput(void) +{ + std::vector out = {getName()}; + + for (unsigned int mu = 0; mu < env().getNd(); ++mu) + for (unsigned int nu = mu; nu < env().getNd(); ++nu) + { + out.push_back(varName(getName(), mu, nu)); + } + + return out; +} + +// setup /////////////////////////////////////////////////////////////////////// +template +void TEMT::setup(void) +{ + envCreateLat(ComplexField, getName()); + envTmpLat(ComplexField, "sumkin"); +} + +// execution /////////////////////////////////////////////////////////////////// +template +void TEMT::execute(void) +{ + LOG(Message) << "Computing energy-momentum tensor" << std::endl; + LOG(Message) << " kinetic terms: '" << par().kinetic << "'" << std::endl; + LOG(Message) << " tr(phi^n): '" << par().phiPow << "'" << std::endl; + LOG(Message) << " improvement: '" << par().improvement << "'" << std::endl; + LOG(Message) << " m^2= " << par().m2 << std::endl; + LOG(Message) << " lambda= " << par().lambda << std::endl; + LOG(Message) << " g= " << par().g << std::endl; + LOG(Message) << " xi= " << par().xi << std::endl; + + const unsigned int N = SImpl::Group::Dimension; + auto &trphi2 = envGet(ComplexField, varName(par().phiPow, 2)); + auto &trphi4 = envGet(ComplexField, varName(par().phiPow, 4)); + + envGetTmp(ComplexField, sumkin); + sumkin = zero; + for (unsigned int mu = 0; mu < env().getNd(); ++mu) + { + auto &trkin = envGet(ComplexField, varName(par().kinetic, mu, mu)); + + sumkin += trkin; + } + for (unsigned int mu = 0; mu < env().getNd(); ++mu) + for (unsigned int nu = mu; nu < env().getNd(); ++nu) + { + auto &out = envGet(ComplexField, varName(getName(), mu, nu)); + auto &trkin = envGet(ComplexField, varName(par().kinetic, mu, nu)); + auto &imp = envGet(ComplexField, varName(par().improvement, mu, nu)); + + out = 2.*trkin + par().xi*imp; + if (mu == nu) + { + out -= sumkin + par().m2*trphi2 + par().lambda*trphi4; + } + out *= N/par().g; + } +} + +END_MODULE_NAMESPACE + +END_HADRONS_NAMESPACE + +#endif // Hadrons_MScalarSUN_EMT_hpp_ diff --git a/extras/Hadrons/modules.inc b/extras/Hadrons/modules.inc index dbfa359c..44df3091 100644 --- a/extras/Hadrons/modules.inc +++ b/extras/Hadrons/modules.inc @@ -46,6 +46,7 @@ modules_hpp =\ Modules/MAction/ZMobiusDWF.hpp \ Modules/MScalarSUN/Div.hpp \ Modules/MScalarSUN/TrMag.hpp \ + Modules/MScalarSUN/EMT.hpp \ Modules/MScalarSUN/TwoPoint.hpp \ Modules/MScalarSUN/TrPhi.hpp \ Modules/MScalarSUN/Utils.hpp \ From 018801d9731e8a01ecdef8c17fa9b773b5c85260 Mon Sep 17 00:00:00 2001 From: Antonin Portelli Date: Fri, 2 Mar 2018 17:56:00 +0000 Subject: [PATCH 31/38] Hadrons: legal update --- extras/Hadrons/Modules/MScalarSUN/EMT.hpp | 27 +++++++++++++++++++ .../Hadrons/Modules/MScalarSUN/TransProj.hpp | 27 +++++++++++++++++++ extras/Hadrons/Modules/MScalarSUN/Utils.hpp | 27 +++++++++++++++++++ 3 files changed, 81 insertions(+) diff --git a/extras/Hadrons/Modules/MScalarSUN/EMT.hpp b/extras/Hadrons/Modules/MScalarSUN/EMT.hpp index 7f1a612d..850efe6a 100644 --- a/extras/Hadrons/Modules/MScalarSUN/EMT.hpp +++ b/extras/Hadrons/Modules/MScalarSUN/EMT.hpp @@ -1,3 +1,30 @@ +/************************************************************************************* + +Grid physics library, www.github.com/paboyle/Grid + +Source file: extras/Hadrons/Modules/MScalarSUN/EMT.hpp + +Copyright (C) 2015-2018 + +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 */ #ifndef Hadrons_MScalarSUN_EMT_hpp_ #define Hadrons_MScalarSUN_EMT_hpp_ diff --git a/extras/Hadrons/Modules/MScalarSUN/TransProj.hpp b/extras/Hadrons/Modules/MScalarSUN/TransProj.hpp index 59b7097c..09386cea 100644 --- a/extras/Hadrons/Modules/MScalarSUN/TransProj.hpp +++ b/extras/Hadrons/Modules/MScalarSUN/TransProj.hpp @@ -1,3 +1,30 @@ +/************************************************************************************* + +Grid physics library, www.github.com/paboyle/Grid + +Source file: extras/Hadrons/Modules/MScalarSUN/TransProj.hpp + +Copyright (C) 2015-2018 + +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 */ #ifndef Hadrons_MScalarSUN_TransProj_hpp_ #define Hadrons_MScalarSUN_TransProj_hpp_ diff --git a/extras/Hadrons/Modules/MScalarSUN/Utils.hpp b/extras/Hadrons/Modules/MScalarSUN/Utils.hpp index 1a13f8f5..1a17282a 100644 --- a/extras/Hadrons/Modules/MScalarSUN/Utils.hpp +++ b/extras/Hadrons/Modules/MScalarSUN/Utils.hpp @@ -1,3 +1,30 @@ +/************************************************************************************* + +Grid physics library, www.github.com/paboyle/Grid + +Source file: extras/Hadrons/Modules/MScalarSUN/Utils.hpp + +Copyright (C) 2015-2018 + +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 */ #ifndef Hadrons_MScalarSUN_Utils_hpp_ #define Hadrons_MScalarSUN_Utils_hpp_ From 2f4dac35313d342a36eb4a09ed3a04db717ff8c4 Mon Sep 17 00:00:00 2001 From: Antonin Portelli Date: Fri, 2 Mar 2018 18:10:58 +0000 Subject: [PATCH 32/38] Hadrons: legal update --- extras/Hadrons/Modules/MGauge/StochEm.hpp | 1 + extras/Hadrons/Modules/MSource/SeqConserved.hpp | 4 ++-- 2 files changed, 3 insertions(+), 2 deletions(-) diff --git a/extras/Hadrons/Modules/MGauge/StochEm.hpp b/extras/Hadrons/Modules/MGauge/StochEm.hpp index 05162443..9e2ba1fb 100644 --- a/extras/Hadrons/Modules/MGauge/StochEm.hpp +++ b/extras/Hadrons/Modules/MGauge/StochEm.hpp @@ -7,6 +7,7 @@ Source file: extras/Hadrons/Modules/MGauge/StochEm.hpp Copyright (C) 2015-2018 Author: Antonin Portelli +Author: Vera Guelpers 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 diff --git a/extras/Hadrons/Modules/MSource/SeqConserved.hpp b/extras/Hadrons/Modules/MSource/SeqConserved.hpp index 5c4a320b..d7780528 100644 --- a/extras/Hadrons/Modules/MSource/SeqConserved.hpp +++ b/extras/Hadrons/Modules/MSource/SeqConserved.hpp @@ -7,8 +7,8 @@ Source file: extras/Hadrons/Modules/MSource/SeqConserved.hpp Copyright (C) 2015-2018 Author: Antonin Portelli -Author: Andrew Lawson -Author: Vera Guelpers +Author: Lanny91 +Author: Vera Guelpers 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 From c4baf876d4896318817dbcc7d5b03680c8142bc9 Mon Sep 17 00:00:00 2001 From: Antonin Portelli Date: Fri, 2 Mar 2018 18:40:18 +0000 Subject: [PATCH 33/38] Hadrons: graph consistency check --- extras/Hadrons/VirtualMachine.cc | 13 ++++++++++++- 1 file changed, 12 insertions(+), 1 deletion(-) diff --git a/extras/Hadrons/VirtualMachine.cc b/extras/Hadrons/VirtualMachine.cc index c91667e7..1bd18100 100644 --- a/extras/Hadrons/VirtualMachine.cc +++ b/extras/Hadrons/VirtualMachine.cc @@ -296,7 +296,18 @@ void VirtualMachine::makeModuleGraph(void) { for (auto &in: module_[m].input) { - graph.addEdge(env().getObjectModule(in), m); + int min = env().getObjectModule(in); + + if (min < 0) + { + HADRON_ERROR(Definition, "object with address " + + std::to_string(in) + + " is not produced by any module"); + } + else + { + graph.addEdge(min, m); + } } } graph_ = graph; From 480708b9a0aac92e00a34bf1ba988f6e0994c215 Mon Sep 17 00:00:00 2001 From: Antonin Portelli Date: Fri, 2 Mar 2018 19:19:37 +0000 Subject: [PATCH 34/38] Hadrons: safer error handling for HadronsXmlRun --- extras/Hadrons/Exceptions.hpp | 1 - extras/Hadrons/HadronsXmlRun.cc | 24 +++++++--- extras/Hadrons/HadronsXmlSchedule.cc | 65 ---------------------------- extras/Hadrons/Makefile.am | 5 +-- 4 files changed, 19 insertions(+), 76 deletions(-) delete mode 100644 extras/Hadrons/HadronsXmlSchedule.cc diff --git a/extras/Hadrons/Exceptions.hpp b/extras/Hadrons/Exceptions.hpp index ab588e5e..ce11618a 100644 --- a/extras/Hadrons/Exceptions.hpp +++ b/extras/Hadrons/Exceptions.hpp @@ -37,7 +37,6 @@ See the full license in the file "LICENSE" in the top level distribution directo #define SRC_LOC std::string(__FUNCTION__) + " at " + std::string(__FILE__) + ":"\ + std::to_string(__LINE__) #define HADRON_ERROR(exc, msg)\ -LOG(Error) << msg << std::endl;\ throw(Exceptions::exc(msg, SRC_LOC)); #define DECL_EXC(name, base) \ diff --git a/extras/Hadrons/HadronsXmlRun.cc b/extras/Hadrons/HadronsXmlRun.cc index 680f234b..4fcfabf2 100644 --- a/extras/Hadrons/HadronsXmlRun.cc +++ b/extras/Hadrons/HadronsXmlRun.cc @@ -56,14 +56,26 @@ int main(int argc, char *argv[]) Grid_init(&argc, &argv); // execution - Application application(parameterFileName); - - application.parseParameterFile(parameterFileName); - if (!scheduleFileName.empty()) + try { - application.loadSchedule(scheduleFileName); + Application application(parameterFileName); + + application.parseParameterFile(parameterFileName); + if (!scheduleFileName.empty()) + { + application.loadSchedule(scheduleFileName); + } + application.run(); + } + catch (const std::exception& e) + { + LOG(Error) << "FATAL ERROR -- Exception " << typeName(&typeid(e)) << std::endl; + LOG(Error) << e.what() << std::endl; + LOG(Error) << "Aborting program" << std::endl; + Grid_finalize(); + + return EXIT_FAILURE; } - application.run(); // epilogue LOG(Message) << "Grid is finalizing now" << std::endl; diff --git a/extras/Hadrons/HadronsXmlSchedule.cc b/extras/Hadrons/HadronsXmlSchedule.cc deleted file mode 100644 index 55f3b231..00000000 --- a/extras/Hadrons/HadronsXmlSchedule.cc +++ /dev/null @@ -1,65 +0,0 @@ -/************************************************************************************* - -Grid physics library, www.github.com/paboyle/Grid - -Source file: extras/Hadrons/HadronsXmlSchedule.cc - -Copyright (C) 2015-2018 - -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 - -using namespace Grid; -using namespace QCD; -using namespace Hadrons; - -int main(int argc, char *argv[]) -{ - // parse command line - std::string parameterFileName, scheduleFileName; - - if (argc < 3) - { - std::cerr << "usage: " << argv[0] << " [Grid options]"; - std::cerr << std::endl; - std::exit(EXIT_FAILURE); - } - parameterFileName = argv[1]; - scheduleFileName = argv[2]; - - // initialization - Grid_init(&argc, &argv); - - // execution - Application application; - - application.parseParameterFile(parameterFileName); - application.schedule(); - application.printSchedule(); - application.saveSchedule(scheduleFileName); - - // epilogue - LOG(Message) << "Grid is finalizing now" << std::endl; - Grid_finalize(); - - return EXIT_SUCCESS; -} diff --git a/extras/Hadrons/Makefile.am b/extras/Hadrons/Makefile.am index 477b93e4..3b393cd1 100644 --- a/extras/Hadrons/Makefile.am +++ b/extras/Hadrons/Makefile.am @@ -1,5 +1,5 @@ lib_LIBRARIES = libHadrons.a -bin_PROGRAMS = HadronsXmlRun HadronsXmlSchedule +bin_PROGRAMS = HadronsXmlRun include modules.inc @@ -29,6 +29,3 @@ nobase_libHadrons_a_HEADERS = \ HadronsXmlRun_SOURCES = HadronsXmlRun.cc HadronsXmlRun_LDADD = libHadrons.a -lGrid - -HadronsXmlSchedule_SOURCES = HadronsXmlSchedule.cc -HadronsXmlSchedule_LDADD = libHadrons.a -lGrid From 90f40009358acf4c8e75888c0425369c3dea1d4d Mon Sep 17 00:00:00 2001 From: Antonin Portelli Date: Fri, 2 Mar 2018 19:20:01 +0000 Subject: [PATCH 35/38] Hadrons: scheduler debug less verbose --- extras/Hadrons/GeneticScheduler.hpp | 8 ++++---- extras/Hadrons/VirtualMachine.cc | 2 +- 2 files changed, 5 insertions(+), 5 deletions(-) diff --git a/extras/Hadrons/GeneticScheduler.hpp b/extras/Hadrons/GeneticScheduler.hpp index 9a6476c3..d6f9bdee 100644 --- a/extras/Hadrons/GeneticScheduler.hpp +++ b/extras/Hadrons/GeneticScheduler.hpp @@ -130,7 +130,7 @@ void GeneticScheduler::nextGeneration(void) { initPopulation(); } - LOG(Debug) << "Starting population:\n" << *this << std::endl; + //LOG(Debug) << "Starting population:\n" << *this << std::endl; // random mutations //PARALLEL_FOR_LOOP @@ -138,7 +138,7 @@ void GeneticScheduler::nextGeneration(void) { doMutation(); } - LOG(Debug) << "After mutations:\n" << *this << std::endl; + //LOG(Debug) << "After mutations:\n" << *this << std::endl; // mating //PARALLEL_FOR_LOOP @@ -146,14 +146,14 @@ void GeneticScheduler::nextGeneration(void) { doCrossover(); } - LOG(Debug) << "After mating:\n" << *this << std::endl; + //LOG(Debug) << "After mating:\n" << *this << std::endl; // grim reaper auto it = population_.begin(); std::advance(it, par_.popSize); population_.erase(it, population_.end()); - LOG(Debug) << "After grim reaper:\n" << *this << std::endl; + //LOG(Debug) << "After grim reaper:\n" << *this << std::endl; } // evolution steps ///////////////////////////////////////////////////////////// diff --git a/extras/Hadrons/VirtualMachine.cc b/extras/Hadrons/VirtualMachine.cc index 1bd18100..1c60aa5b 100644 --- a/extras/Hadrons/VirtualMachine.cc +++ b/extras/Hadrons/VirtualMachine.cc @@ -549,7 +549,7 @@ VirtualMachine::Program VirtualMachine::schedule(const GeneticPar &par) gen = 0; do { - LOG(Debug) << "Generation " << gen << ":" << std::endl; + //LOG(Debug) << "Generation " << gen << ":" << std::endl; scheduler.nextGeneration(); if (gen != 0) { From fcac5c07722a3fb200a6e2ae62e998c526f44473 Mon Sep 17 00:00:00 2001 From: Antonin Portelli Date: Fri, 2 Mar 2018 19:20:23 +0000 Subject: [PATCH 36/38] Hadrons: scalar SU(N) fixes --- extras/Hadrons/Modules/MScalarSUN/EMT.hpp | 8 ++++++-- extras/Hadrons/Modules/MScalarSUN/TransProj.hpp | 2 +- 2 files changed, 7 insertions(+), 3 deletions(-) diff --git a/extras/Hadrons/Modules/MScalarSUN/EMT.hpp b/extras/Hadrons/Modules/MScalarSUN/EMT.hpp index 850efe6a..8c1239df 100644 --- a/extras/Hadrons/Modules/MScalarSUN/EMT.hpp +++ b/extras/Hadrons/Modules/MScalarSUN/EMT.hpp @@ -110,7 +110,7 @@ std::vector TEMT::getInput(void) template std::vector TEMT::getOutput(void) { - std::vector out = {getName()}; + std::vector out; for (unsigned int mu = 0; mu < env().getNd(); ++mu) for (unsigned int nu = mu; nu < env().getNd(); ++nu) @@ -125,7 +125,11 @@ std::vector TEMT::getOutput(void) template void TEMT::setup(void) { - envCreateLat(ComplexField, getName()); + for (unsigned int mu = 0; mu < env().getNd(); ++mu) + for (unsigned int nu = mu; nu < env().getNd(); ++nu) + { + envCreateLat(ComplexField, varName(getName(), mu, nu)); + } envTmpLat(ComplexField, "sumkin"); } diff --git a/extras/Hadrons/Modules/MScalarSUN/TransProj.hpp b/extras/Hadrons/Modules/MScalarSUN/TransProj.hpp index 09386cea..4532b4ab 100644 --- a/extras/Hadrons/Modules/MScalarSUN/TransProj.hpp +++ b/extras/Hadrons/Modules/MScalarSUN/TransProj.hpp @@ -104,7 +104,7 @@ std::vector TTransProj::getInput(void) template std::vector TTransProj::getOutput(void) { - std::vector out = {getName()}; + std::vector out; for (unsigned int mu = 0; mu < env().getNd(); ++mu) for (unsigned int nu = mu; nu < env().getNd(); ++nu) From 2e88408f5ce1bc1ba4052be07c4b1e94f0a99f5a Mon Sep 17 00:00:00 2001 From: Fionn O hOgain Date: Fri, 2 Mar 2018 22:27:41 +0000 Subject: [PATCH 37/38] Some changes needed for deflation interface --- lib/algorithms/iterative/Deflation.h | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/lib/algorithms/iterative/Deflation.h b/lib/algorithms/iterative/Deflation.h index b6aa0d3d..b2239c55 100644 --- a/lib/algorithms/iterative/Deflation.h +++ b/lib/algorithms/iterative/Deflation.h @@ -59,7 +59,7 @@ public: assert(evec.size()==eval.size()); auto N = evec.size(); for (int i=0;i Date: Sat, 3 Mar 2018 13:42:36 +0000 Subject: [PATCH 38/38] Hadrons: basic GraphViz visualisation --- extras/Hadrons/VirtualMachine.cc | 43 +++++++++++++++++++++++++++++++ extras/Hadrons/VirtualMachine.hpp | 6 ++++- 2 files changed, 48 insertions(+), 1 deletion(-) diff --git a/extras/Hadrons/VirtualMachine.cc b/extras/Hadrons/VirtualMachine.cc index 1c60aa5b..eb1f68ba 100644 --- a/extras/Hadrons/VirtualMachine.cc +++ b/extras/Hadrons/VirtualMachine.cc @@ -111,6 +111,7 @@ void VirtualMachine::pushModule(VirtualMachine::ModPt &pt) { // output does not exists, add it env().addObject(out, address); + module_[address].output.push_back(env().getObjectAddress(out)); } else { @@ -313,6 +314,48 @@ void VirtualMachine::makeModuleGraph(void) graph_ = graph; } +// dump GraphViz graph ///////////////////////////////////////////////////////// +void VirtualMachine::dumpModuleGraph(std::ostream &out) +{ + makeModuleGraph(); + out << "digraph hadrons {" << std::endl; + out << "node [shape=record, fontname=\"Courier\", fontsize=\"11\"];" << std::endl; + out << "graph [fontname = \"Courier\", fontsize=\"11\"];" << std::endl; + out << "edge [fontname = \"Courier\", fontsize=\"11\"];"<< std::endl; + for (unsigned int m = 0; m < module_.size(); ++m) + { + + } + for (unsigned int m = 0; m < module_.size(); ++m) + { + for (auto &in: module_[m].input) + { + int min = env().getObjectModule(in); + + out << min << " -> " << m << " [ label = \"" + << env().getObjectName(in) << "\" ];" << std::endl; + } + } + for (unsigned int m = 0; m < module_.size(); ++m) + { + out << m << " [ label = \"{ " << getModule(m)->getRegisteredName() + << " | " << getModuleName(m) << "}\" ];" << std::endl; + } + out << "}\n" << std::endl; +} + +void VirtualMachine::dumpModuleGraph(void) +{ + dumpModuleGraph(std::cout); +} + +void VirtualMachine::dumpModuleGraph(const std::string filename) +{ + std::ofstream f(filename); + + dumpModuleGraph(f); +} + // memory profile ////////////////////////////////////////////////////////////// const VirtualMachine::MemoryProfile & VirtualMachine::getMemoryProfile(void) { diff --git a/extras/Hadrons/VirtualMachine.hpp b/extras/Hadrons/VirtualMachine.hpp index 19a74f94..68eeb0c0 100644 --- a/extras/Hadrons/VirtualMachine.hpp +++ b/extras/Hadrons/VirtualMachine.hpp @@ -84,7 +84,7 @@ private: const std::type_info *type{nullptr}; std::string name; ModPt data{nullptr}; - std::vector input; + std::vector input, output; size_t maxAllocated; }; public: @@ -120,6 +120,10 @@ public: void printContent(void) const; // module graph (could be a const reference if topoSort was const) Graph getModuleGraph(void); + // dump GraphViz graph + void dumpModuleGraph(std::ostream &out); + void dumpModuleGraph(void); + void dumpModuleGraph(const std::string filename); // memory profile const MemoryProfile &getMemoryProfile(void); // garbage collector