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/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/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/Global.cc b/extras/Hadrons/Global.cc index 9a90a08c..0164a1ab 100644 --- a/extras/Hadrons/Global.cc +++ b/extras/Hadrons/Global.cc @@ -37,20 +37,38 @@ 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) { - 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/extras/Hadrons/Global.hpp b/extras/Hadrons/Global.hpp index fc069ed6..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 @@ -111,6 +114,7 @@ extern HadronsLogger HadronsLogWarning; extern HadronsLogger HadronsLogMessage; extern HadronsLogger HadronsLogIterative; extern HadronsLogger HadronsLogDebug; +extern HadronsLogger HadronsLogIRL; void initLogger(void); @@ -180,6 +184,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/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/LanczosUtils.hpp b/extras/Hadrons/LanczosUtils.hpp new file mode 100644 index 00000000..a080da4b --- /dev/null +++ b/extras/Hadrons/LanczosUtils.hpp @@ -0,0 +1,115 @@ +/************************************************************************************* + +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_ + +#include +#include + +BEGIN_HADRONS_NAMESPACE + +// Lanczos type +#ifndef HADRONS_DEFAULT_LANCZOS_NBASIS +#define HADRONS_DEFAULT_LANCZOS_NBASIS 60 +#endif + +template +struct EigenPack +{ + typedef T VectorType; + 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< + typename LocalCoherenceLanczos::CoarseField>; + +END_HADRONS_NAMESPACE + +#endif // Hadrons_LanczosUtils_hpp_ diff --git a/extras/Hadrons/Makefile.am b/extras/Hadrons/Makefile.am index 3d07679a..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 @@ -21,6 +21,7 @@ nobase_libHadrons_a_HEADERS = \ GeneticScheduler.hpp \ Global.hpp \ Graph.hpp \ + LanczosUtils.hpp \ Module.hpp \ Modules.hpp \ ModuleFactory.hpp \ @@ -28,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 diff --git a/extras/Hadrons/Modules.hpp b/extras/Hadrons/Modules.hpp index 351860fb..04cc7a62 100644 --- a/extras/Hadrons/Modules.hpp +++ b/extras/Hadrons/Modules.hpp @@ -7,7 +7,9 @@ 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 @@ -43,12 +45,13 @@ See the full license in the file "LICENSE" in the top level distribution directo #include #include #include +#include #include #include #include +#include #include #include -#include #include #include #include @@ -60,9 +63,14 @@ 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 #include diff --git a/extras/Hadrons/Modules/MAction/WilsonClover.hpp b/extras/Hadrons/Modules/MAction/WilsonClover.hpp index c369f086..2967931c 100644 --- a/extras/Hadrons/Modules/MAction/WilsonClover.hpp +++ b/extras/Hadrons/Modules/MAction/WilsonClover.hpp @@ -2,12 +2,13 @@ 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/MAction/ZMobiusDWF.hpp b/extras/Hadrons/Modules/MAction/ZMobiusDWF.hpp new file mode 100644 index 00000000..518b88f2 --- /dev/null +++ b/extras/Hadrons/Modules/MAction/ZMobiusDWF.hpp @@ -0,0 +1,143 @@ +/************************************************************************************* + +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_ + +#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= " << par().Ls + << ", b= " << par().b << ", c= " << par().c + << " 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/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..508145bc 100644 --- a/extras/Hadrons/Modules/MGauge/FundtoHirep.hpp +++ b/extras/Hadrons/Modules/MGauge/FundtoHirep.hpp @@ -4,11 +4,10 @@ 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: Antonin Portelli +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/StochEm.hpp b/extras/Hadrons/Modules/MGauge/StochEm.hpp index 8ed58c4b..fc09a0ee 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/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/MScalarSUN/Div.hpp b/extras/Hadrons/Modules/MScalarSUN/Div.hpp index cab6e88c..0ecb1dbe 100644 --- a/extras/Hadrons/Modules/MScalarSUN/Div.hpp +++ b/extras/Hadrons/Modules/MScalarSUN/Div.hpp @@ -31,18 +31,18 @@ See the full license in the file "LICENSE" in the top level distribution directo #include #include #include +#include BEGIN_HADRONS_NAMESPACE /****************************************************************************** - * Div * + * Divergence of a vector field * ******************************************************************************/ 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 @@ -83,7 +83,7 @@ MODULE_REGISTER_NS(DivSU5, TDiv>, MScalarSUN); MODULE_REGISTER_NS(DivSU6, TDiv>, MScalarSUN); /****************************************************************************** - * TDiv implementation * + * TDiv implementation * ******************************************************************************/ // constructor ///////////////////////////////////////////////////////////////// template @@ -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/EMT.hpp b/extras/Hadrons/Modules/MScalarSUN/EMT.hpp new file mode 100644 index 00000000..8c1239df --- /dev/null +++ b/extras/Hadrons/Modules/MScalarSUN/EMT.hpp @@ -0,0 +1,181 @@ +/************************************************************************************* + +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_ + +#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; + + 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) +{ + 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"); +} + +// 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/MScalarSUN/TrKinetic.hpp b/extras/Hadrons/Modules/MScalarSUN/TrKinetic.hpp new file mode 100644 index 00000000..f1b520de --- /dev/null +++ b/extras/Hadrons/Modules/MScalarSUN/TrKinetic.hpp @@ -0,0 +1,170 @@ +/************************************************************************************* + +Grid physics library, www.github.com/paboyle/Grid + +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 +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 +#include + +BEGIN_HADRONS_NAMESPACE + +/****************************************************************************** + * Trace of kinetic term * + ******************************************************************************/ +BEGIN_MODULE_NAMESPACE(MScalarSUN) + +class TrKineticPar: Serializable +{ +public: + 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); +}; + +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(varName(getName(), 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, varName(getName(), 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) + { + 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, varName(getName(), 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); + } +} + +END_MODULE_NAMESPACE + +END_HADRONS_NAMESPACE + +#endif // Hadrons_MScalarSUN_TrKinetic_hpp_ diff --git a/extras/Hadrons/Modules/MScalarSUN/TrMag.hpp b/extras/Hadrons/Modules/MScalarSUN/TrMag.hpp index 2858984c..879951d5 100644 --- a/extras/Hadrons/Modules/MScalarSUN/TrMag.hpp +++ b/extras/Hadrons/Modules/MScalarSUN/TrMag.hpp @@ -31,11 +31,12 @@ See the full license in the file "LICENSE" in the top level distribution directo #include #include #include +#include BEGIN_HADRONS_NAMESPACE /****************************************************************************** - * Module to compute tr(mag^n) * + * Trace of powers of the magnetisation * ******************************************************************************/ BEGIN_MODULE_NAMESPACE(MScalarSUN) @@ -117,7 +118,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..aa98bfb3 100644 --- a/extras/Hadrons/Modules/MScalarSUN/TrPhi.hpp +++ b/extras/Hadrons/Modules/MScalarSUN/TrPhi.hpp @@ -31,11 +31,12 @@ See the full license in the file "LICENSE" in the top level distribution directo #include #include #include +#include BEGIN_HADRONS_NAMESPACE /****************************************************************************** - * Module to compute tr(phi^n) * + * Trace of powers of a scalar field * ******************************************************************************/ BEGIN_MODULE_NAMESPACE(MScalarSUN) @@ -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)); } } @@ -136,7 +134,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); @@ -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/TransProj.hpp b/extras/Hadrons/Modules/MScalarSUN/TransProj.hpp new file mode 100644 index 00000000..4532b4ab --- /dev/null +++ b/extras/Hadrons/Modules/MScalarSUN/TransProj.hpp @@ -0,0 +1,185 @@ +/************************************************************************************* + +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_ + +#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; + + 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/MScalarSUN/TwoPoint.hpp b/extras/Hadrons/Modules/MScalarSUN/TwoPoint.hpp index 58570dd5..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 @@ -87,7 +88,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/MScalarSUN/Utils.hpp b/extras/Hadrons/Modules/MScalarSUN/Utils.hpp new file mode 100644 index 00000000..1a17282a --- /dev/null +++ b/extras/Hadrons/Modules/MScalarSUN/Utils.hpp @@ -0,0 +1,107 @@ +/************************************************************************************* + +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_ + +#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/MSolver/LocalCoherenceLanczos.hpp b/extras/Hadrons/Modules/MSolver/LocalCoherenceLanczos.hpp new file mode 100644 index 00000000..6e2103ce --- /dev/null +++ b/extras/Hadrons/Modules/MSolver/LocalCoherenceLanczos.hpp @@ -0,0 +1,249 @@ +/************************************************************************************* + +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_ + +#include +#include +#include +#include + +BEGIN_HADRONS_NAMESPACE + +/****************************************************************************** + * Local coherence Lanczos eigensolver * + *****************************************************************************/ +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, + std::string, output); +}; + +template +class TLocalCoherenceLanczos: public Module +{ +public: + FERM_TYPE_ALIASES(FImpl,); + typedef LocalCoherenceLanczos LCL; + typedef FineEigenPack FinePack; + typedef CoarseEigenPack CoarsePack; + typedef HADRONS_DEFAULT_SCHUR_OP 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 fineName_, coarseName_; +}; + +MODULE_REGISTER_NS(LocalCoherenceLanczos, + ARG(TLocalCoherenceLanczos), + MSolver); +MODULE_REGISTER_NS(ZLocalCoherenceLanczos, + ARG(TLocalCoherenceLanczos), + MSolver); + +/****************************************************************************** + * TLocalCoherenceLanczos implementation * + ******************************************************************************/ +// constructor ///////////////////////////////////////////////////////////////// +template +TLocalCoherenceLanczos::TLocalCoherenceLanczos(const std::string name) +: Module(name) +{ + fineName_ = getName() + "_fine"; + coarseName_ = getName() + "_coarse"; +} + +// dependencies/products /////////////////////////////////////////////////////// +template +std::vector TLocalCoherenceLanczos::getInput(void) +{ + std::vector in = {par().action}; + + return in; +} + +template +std::vector TLocalCoherenceLanczos::getOutput(void) +{ + std::vector out = {fineName_, coarseName_}; + + 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(); + } + 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_); + 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, + Odd, fine.evec, coarse.evec, fine.eval, coarse.eval); +} + +// execution /////////////////////////////////////////////////////////////////// +template +void TLocalCoherenceLanczos::execute(void) +{ + auto &finePar = par().fineParams; + auto &coarsePar = par().coarseParams; + auto &fine = envGet(FinePack, fineName_); + auto &coarse = envGet(CoarsePack, coarseName_); + + envGetTmp(LCL, solver); + if (par().doFine) + { + LOG(Message) << "Performing fine grid IRL -- Nstop= " + << 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); + solver.testFine(finePar.resid*100.0); + 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= " + << 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); + solver.testCoarse(coarsePar.resid*100.0, par().smoother, + par().coarseRelaxTol); + if (!par().output.empty()) + { + coarse.write(par().output + "_coarse"); + } + } +} + +END_MODULE_NAMESPACE + +END_HADRONS_NAMESPACE + +#endif // Hadrons_MSolver_LocalCoherenceLanczos_hpp_ diff --git a/extras/Hadrons/Modules/MSolver/RBPrecCG.hpp b/extras/Hadrons/Modules/MSolver/RBPrecCG.hpp index 54c0f2d8..2b914625 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 @@ -69,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 * @@ -117,14 +119,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) diff --git a/extras/Hadrons/Modules/MSource/SeqConserved.hpp b/extras/Hadrons/Modules/MSource/SeqConserved.hpp index ee8d8d56..d7780528 100644 --- a/extras/Hadrons/Modules/MSource/SeqConserved.hpp +++ b/extras/Hadrons/Modules/MSource/SeqConserved.hpp @@ -8,6 +8,7 @@ Copyright (C) 2015-2018 Author: Antonin Portelli 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 @@ -38,9 +39,11 @@ 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 +51,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 +72,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 +83,8 @@ class TSeqConserved: public Module { public: FERM_TYPE_ALIASES(FImpl,); +public: + typedef PhotonR::GaugeField EmField; public: // constructor TSeqConserved(const std::string name); @@ -89,10 +98,14 @@ protected: virtual void setup(void); // execution virtual void execute(void); +private: + bool SeqhasPhase_{false}; + std::string SeqmomphName_; }; MODULE_REGISTER_NS(SeqConserved, TSeqConserved, MSource); + /****************************************************************************** * TSeqConserved implementation * ******************************************************************************/ @@ -100,6 +113,7 @@ MODULE_REGISTER_NS(SeqConserved, TSeqConserved, MSource); template TSeqConserved::TSeqConserved(const std::string name) : Module(name) +, SeqmomphName_ (name + "_Seqmomph") {} // dependencies/products /////////////////////////////////////////////////////// @@ -107,7 +121,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 +131,7 @@ std::vector TSeqConserved::getOutput(void) { std::vector out = {getName()}; - return out; + return out; } // setup /////////////////////////////////////////////////////////////////////// @@ -125,6 +140,10 @@ void TSeqConserved::setup(void) { auto Ls_ = env().getObjectLs(par().action); envCreateLat(PropagatorField, getName(), Ls_); + envTmpLat(PropagatorField, "src_tmp"); + envCacheLat(LatticeComplex, SeqmomphName_); + envTmpLat(LatticeComplex, "coor"); + envTmpLat(LatticeComplex, "latt_compl"); } // execution /////////////////////////////////////////////////////////////////// @@ -134,27 +153,79 @@ 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 << " 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 insertion (mu = " - << par().mu << ") for " << par().tA << " <= t <= " - << par().tB << std::endl; + << par().curr_type << " current for " + << par().tA << " <= t <= " + << 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, latt_compl); - std::vector mom = strToVec(par().mom); - mat.SeqConservedCurrent(q, src, par().curr_type, par().mu, - mom, par().tA, par().tB); + src = zero; + + //exp(ipx) + 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; + } + LOG(Message) << "Inserting momentum " << strToVec(par().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/VirtualMachine.cc b/extras/Hadrons/VirtualMachine.cc index d47bafb7..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 { @@ -296,12 +297,65 @@ 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; } +// 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) { @@ -424,11 +478,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); @@ -532,7 +592,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) { 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 diff --git a/extras/Hadrons/modules.inc b/extras/Hadrons/modules.inc index 9f642c5b..f7c6c414 100644 --- a/extras/Hadrons/modules.inc +++ b/extras/Hadrons/modules.inc @@ -31,12 +31,13 @@ modules_hpp =\ Modules/MSource/SeqConserved.hpp \ Modules/MSink/Smear.hpp \ Modules/MSink/Point.hpp \ + Modules/MSolver/LocalCoherenceLanczos.hpp \ Modules/MSolver/RBPrecCG.hpp \ Modules/MGauge/Random.hpp \ + Modules/MGauge/FundtoHirep.hpp \ Modules/MGauge/StochEm.hpp \ Modules/MGauge/Unit.hpp \ Modules/MGauge/UnitEm.hpp \ - Modules/MGauge/FundtoHirep.hpp \ Modules/MUtilities/TestSeqGamma.hpp \ Modules/MUtilities/TestSeqConserved.hpp \ Modules/MLoop/NoiseLoop.hpp \ @@ -48,10 +49,15 @@ 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/EMT.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 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 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 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); 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(); diff --git a/lib/qcd/action/fermion/FermionOperator.h b/lib/qcd/action/fermion/FermionOperator.h index 1d395d53..9eb74baa 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, + ComplexField &lattice_cmplx)=0; }; } diff --git a/lib/qcd/action/fermion/ImprovedStaggeredFermion.cc b/lib/qcd/action/fermion/ImprovedStaggeredFermion.cc index 5005b1a2..ed07fb5b 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, + ComplexField &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..d2426ef7 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, + ComplexField &lattice_cmplx); }; typedef ImprovedStaggeredFermion ImprovedStaggeredFermionF; diff --git a/lib/qcd/action/fermion/ImprovedStaggeredFermion5D.cc b/lib/qcd/action/fermion/ImprovedStaggeredFermion5D.cc index 7308ff30..3510b8f1 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, + ComplexField &lattice_cmplx) { assert(0); + } FermOpStaggeredTemplateInstantiate(ImprovedStaggeredFermion5D); diff --git a/lib/qcd/action/fermion/ImprovedStaggeredFermion5D.h b/lib/qcd/action/fermion/ImprovedStaggeredFermion5D.h index 22fb9e7d..e9cf2013 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, + ComplexField &lattice_cmplx); }; }} diff --git a/lib/qcd/action/fermion/WilsonFermion.cc b/lib/qcd/action/fermion/WilsonFermion.cc index dfaa6758..2d9cf22d 100644 --- a/lib/qcd/action/fermion/WilsonFermion.cc +++ b/lib/qcd/action/fermion/WilsonFermion.cc @@ -407,40 +407,30 @@ 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) + unsigned int tmax, + ComplexField &lattice_cmplx) { conformable(_grid, q_in._grid); conformable(_grid, q_out._grid); - Lattice> ph(_grid), coor(_grid); - ComplexD 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((RealD)(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) @@ -475,6 +465,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 c0db827e..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, - 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, + ComplexField &lattice_cmplx); }; typedef WilsonFermion WilsonFermionF; diff --git a/lib/qcd/action/fermion/WilsonFermion5D.cc b/lib/qcd/action/fermion/WilsonFermion5D.cc index 3e58fed6..6f82aad2 100644 --- a/lib/qcd/action/fermion/WilsonFermion5D.cc +++ b/lib/qcd/action/fermion/WilsonFermion5D.cc @@ -779,92 +779,89 @@ 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) + unsigned int tmax, + ComplexField &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()); - ComplexD i(0.0, 1.0); + PropagatorField tmp(GaugeGrid()),tmp2(GaugeGrid()); 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((RealD)(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) + for (unsigned int s = 0; s < LLs; ++s) { - // 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))); - Kernels::SeqConservedCurrentSiteFwd(tmpFwd._odata[sF], - q_out._odata[sF], Umu, sU, - mu, t_mask, axial_sign); - ++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); } - } + } } } + + + + FermOpTemplateInstantiate(WilsonFermion5D); GparityFermOpTemplateInstantiate(WilsonFermion5D); diff --git a/lib/qcd/action/fermion/WilsonFermion5D.h b/lib/qcd/action/fermion/WilsonFermion5D.h index d66f4a1d..21da4c31 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, + ComplexField &lattice_cmplx); }; }} 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; +} 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); }