diff --git a/.gitignore b/.gitignore index 399f2f6b..dc59879f 100644 --- a/.gitignore +++ b/.gitignore @@ -123,3 +123,9 @@ make-bin-BUCK.sh ##################### lib/qcd/spin/gamma-gen/*.h lib/qcd/spin/gamma-gen/*.cc + +# vs code editor files # +######################## +.vscode/ +.vscode/settings.json +settings.json diff --git a/.travis.yml b/.travis.yml index 7d8203ce..ad4e5b73 100644 --- a/.travis.yml +++ b/.travis.yml @@ -44,3 +44,4 @@ script: - make -j4 - ./benchmarks/Benchmark_dwf --threads 1 --debug-signals - make check + diff --git a/extras/Hadrons/Application.cc b/extras/Hadrons/Application.cc index 9a3366d4..7ba98ade 100644 --- a/extras/Hadrons/Application.cc +++ b/extras/Hadrons/Application.cc @@ -4,8 +4,7 @@ Grid physics library, www.github.com/paboyle/Grid Source file: extras/Hadrons/Application.cc -Copyright (C) 2015 -Copyright (C) 2016 +Copyright (C) 2015-2018 Author: Antonin Portelli @@ -43,6 +42,7 @@ using namespace Hadrons; // constructors //////////////////////////////////////////////////////////////// Application::Application(void) { + initLogger(); LOG(Message) << "Modules available:" << std::endl; auto list = ModuleFactory::getInstance().getBuilderList(); for (auto &m: list) diff --git a/extras/Hadrons/Application.hpp b/extras/Hadrons/Application.hpp index 4b2ce77b..8cd15433 100644 --- a/extras/Hadrons/Application.hpp +++ b/extras/Hadrons/Application.hpp @@ -4,8 +4,7 @@ Grid physics library, www.github.com/paboyle/Grid Source file: extras/Hadrons/Application.hpp -Copyright (C) 2015 -Copyright (C) 2016 +Copyright (C) 2015-2018 Author: Antonin Portelli diff --git a/extras/Hadrons/Environment.cc b/extras/Hadrons/Environment.cc index 6de13e86..82b0dda1 100644 --- a/extras/Hadrons/Environment.cc +++ b/extras/Hadrons/Environment.cc @@ -4,8 +4,7 @@ Grid physics library, www.github.com/paboyle/Grid Source file: extras/Hadrons/Environment.cc -Copyright (C) 2015 -Copyright (C) 2016 +Copyright (C) 2015-2018 Author: Antonin Portelli diff --git a/extras/Hadrons/Environment.hpp b/extras/Hadrons/Environment.hpp index 7f1bc26d..e9bfffe1 100644 --- a/extras/Hadrons/Environment.hpp +++ b/extras/Hadrons/Environment.hpp @@ -4,8 +4,7 @@ Grid physics library, www.github.com/paboyle/Grid Source file: extras/Hadrons/Environment.hpp -Copyright (C) 2015 -Copyright (C) 2016 +Copyright (C) 2015-2018 Author: Antonin Portelli diff --git a/extras/Hadrons/Exceptions.cc b/extras/Hadrons/Exceptions.cc index bf532c21..eedc03b1 100644 --- a/extras/Hadrons/Exceptions.cc +++ b/extras/Hadrons/Exceptions.cc @@ -4,7 +4,7 @@ Grid physics library, www.github.com/paboyle/Grid Source file: extras/Hadrons/Exceptions.cc -Copyright (C) 2017 +Copyright (C) 2015-2018 Author: Antonin Portelli @@ -54,4 +54,4 @@ CONST_EXC(Io, Runtime("IO error: " + msg, loc)) CONST_EXC(Memory, Runtime("memory error: " + msg, loc)) CONST_EXC(Parsing, Runtime("parsing error: " + msg, loc)) CONST_EXC(Program, Runtime("program error: " + msg, loc)) -CONST_EXC(System, Runtime("system error: " + msg, loc)) \ No newline at end of file +CONST_EXC(System, Runtime("system error: " + msg, loc)) diff --git a/extras/Hadrons/Exceptions.hpp b/extras/Hadrons/Exceptions.hpp index 8f04ab41..ab588e5e 100644 --- a/extras/Hadrons/Exceptions.hpp +++ b/extras/Hadrons/Exceptions.hpp @@ -4,7 +4,7 @@ Grid physics library, www.github.com/paboyle/Grid Source file: extras/Hadrons/Exceptions.hpp -Copyright (C) 2017 +Copyright (C) 2015-2018 Author: Antonin Portelli diff --git a/extras/Hadrons/Factory.hpp b/extras/Hadrons/Factory.hpp index 65ce03ca..705a639e 100644 --- a/extras/Hadrons/Factory.hpp +++ b/extras/Hadrons/Factory.hpp @@ -4,8 +4,7 @@ Grid physics library, www.github.com/paboyle/Grid Source file: extras/Hadrons/Factory.hpp -Copyright (C) 2015 -Copyright (C) 2016 +Copyright (C) 2015-2018 Author: Antonin Portelli diff --git a/extras/Hadrons/GeneticScheduler.hpp b/extras/Hadrons/GeneticScheduler.hpp index f199f1ed..9a6476c3 100644 --- a/extras/Hadrons/GeneticScheduler.hpp +++ b/extras/Hadrons/GeneticScheduler.hpp @@ -4,8 +4,7 @@ Grid physics library, www.github.com/paboyle/Grid Source file: extras/Hadrons/GeneticScheduler.hpp -Copyright (C) 2015 -Copyright (C) 2016 +Copyright (C) 2015-2018 Author: Antonin Portelli diff --git a/extras/Hadrons/Global.cc b/extras/Hadrons/Global.cc index 130ede96..9a90a08c 100644 --- a/extras/Hadrons/Global.cc +++ b/extras/Hadrons/Global.cc @@ -4,8 +4,7 @@ Grid physics library, www.github.com/paboyle/Grid Source file: extras/Hadrons/Global.cc -Copyright (C) 2015 -Copyright (C) 2016 +Copyright (C) 2015-2018 Author: Antonin Portelli @@ -39,6 +38,21 @@ HadronsLogger Hadrons::HadronsLogMessage(1,"Message"); HadronsLogger Hadrons::HadronsLogIterative(1,"Iterative"); HadronsLogger Hadrons::HadronsLogDebug(1,"Debug"); +void Hadrons::initLogger(void) +{ + auto w = std::string("Hadrons").length(); + GridLogError.setTopWidth(w); + GridLogWarning.setTopWidth(w); + GridLogMessage.setTopWidth(w); + GridLogIterative.setTopWidth(w); + GridLogDebug.setTopWidth(w); + HadronsLogError.Active(GridLogError.isActive()); + HadronsLogWarning.Active(GridLogWarning.isActive()); + HadronsLogMessage.Active(GridLogMessage.isActive()); + HadronsLogIterative.Active(GridLogIterative.isActive()); + HadronsLogDebug.Active(GridLogDebug.isActive()); +} + // type utilities ////////////////////////////////////////////////////////////// constexpr unsigned int maxNameSize = 1024u; @@ -53,3 +67,10 @@ std::string Hadrons::typeName(const std::type_info *info) return name; } + +// default writers/readers ///////////////////////////////////////////////////// +#ifdef HAVE_HDF5 +const std::string Hadrons::resultFileExt = "h5"; +#else +const std::string Hadrons::resultFileExt = "xml"; +#endif diff --git a/extras/Hadrons/Global.hpp b/extras/Hadrons/Global.hpp index ebfe94dc..fc069ed6 100644 --- a/extras/Hadrons/Global.hpp +++ b/extras/Hadrons/Global.hpp @@ -4,10 +4,10 @@ Grid physics library, www.github.com/paboyle/Grid Source file: extras/Hadrons/Global.hpp -Copyright (C) 2015 -Copyright (C) 2016 +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 @@ -61,6 +61,9 @@ using Grid::operator<<; #ifndef SIMPL #define SIMPL ScalarImplCR #endif +#ifndef GIMPL +#define GIMPL GimplTypesR +#endif BEGIN_HADRONS_NAMESPACE @@ -84,7 +87,8 @@ typedef std::function SolverFn##suffix; #define SINK_TYPE_ALIASES(suffix)\ -typedef std::function SinkFn##suffix; +typedef std::function SinkFn##suffix; #define FGS_TYPE_ALIASES(FImpl, suffix)\ FERM_TYPE_ALIASES(FImpl, suffix)\ @@ -108,6 +112,8 @@ extern HadronsLogger HadronsLogMessage; extern HadronsLogger HadronsLogIterative; extern HadronsLogger HadronsLogDebug; +void initLogger(void); + // singleton pattern #define SINGLETON(name)\ public:\ @@ -161,14 +167,19 @@ std::string typeName(void) } // default writers/readers +extern const std::string resultFileExt; + #ifdef HAVE_HDF5 -typedef Hdf5Reader CorrReader; -typedef Hdf5Writer CorrWriter; +typedef Hdf5Reader ResultReader; +typedef Hdf5Writer ResultWriter; #else -typedef XmlReader CorrReader; -typedef XmlWriter CorrWriter; +typedef XmlReader ResultReader; +typedef XmlWriter ResultWriter; #endif +#define RESULT_FILE_NAME(name) \ +name + "." + std::to_string(vm().getTrajectory()) + "." + resultFileExt + END_HADRONS_NAMESPACE #include diff --git a/extras/Hadrons/Graph.hpp b/extras/Hadrons/Graph.hpp index a9c240fa..67694aa8 100644 --- a/extras/Hadrons/Graph.hpp +++ b/extras/Hadrons/Graph.hpp @@ -4,8 +4,7 @@ Grid physics library, www.github.com/paboyle/Grid Source file: extras/Hadrons/Graph.hpp -Copyright (C) 2015 -Copyright (C) 2016 +Copyright (C) 2015-2018 Author: Antonin Portelli diff --git a/extras/Hadrons/HadronsXmlRun.cc b/extras/Hadrons/HadronsXmlRun.cc index 0dff8f9a..680f234b 100644 --- a/extras/Hadrons/HadronsXmlRun.cc +++ b/extras/Hadrons/HadronsXmlRun.cc @@ -4,8 +4,7 @@ Grid physics library, www.github.com/paboyle/Grid Source file: extras/Hadrons/HadronsXmlRun.cc -Copyright (C) 2015 -Copyright (C) 2016 +Copyright (C) 2015-2018 Author: Antonin Portelli @@ -55,12 +54,6 @@ 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; // execution Application application(parameterFileName); diff --git a/extras/Hadrons/HadronsXmlSchedule.cc b/extras/Hadrons/HadronsXmlSchedule.cc index a8ca9a63..55f3b231 100644 --- a/extras/Hadrons/HadronsXmlSchedule.cc +++ b/extras/Hadrons/HadronsXmlSchedule.cc @@ -4,8 +4,7 @@ Grid physics library, www.github.com/paboyle/Grid Source file: extras/Hadrons/HadronsXmlSchedule.cc -Copyright (C) 2015 -Copyright (C) 2016 +Copyright (C) 2015-2018 Author: Antonin Portelli @@ -49,12 +48,6 @@ 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; // execution Application application; diff --git a/extras/Hadrons/Module.cc b/extras/Hadrons/Module.cc index e5ef0fe4..54978f93 100644 --- a/extras/Hadrons/Module.cc +++ b/extras/Hadrons/Module.cc @@ -4,8 +4,7 @@ Grid physics library, www.github.com/paboyle/Grid Source file: extras/Hadrons/Module.cc -Copyright (C) 2015 -Copyright (C) 2016 +Copyright (C) 2015-2018 Author: Antonin Portelli diff --git a/extras/Hadrons/Module.hpp b/extras/Hadrons/Module.hpp index b71f779d..2ba425e4 100644 --- a/extras/Hadrons/Module.hpp +++ b/extras/Hadrons/Module.hpp @@ -4,8 +4,7 @@ Grid physics library, www.github.com/paboyle/Grid Source file: extras/Hadrons/Module.hpp -Copyright (C) 2015 -Copyright (C) 2016 +Copyright (C) 2015-2018 Author: Antonin Portelli diff --git a/extras/Hadrons/ModuleFactory.hpp b/extras/Hadrons/ModuleFactory.hpp index 48ab305c..d5c703fa 100644 --- a/extras/Hadrons/ModuleFactory.hpp +++ b/extras/Hadrons/ModuleFactory.hpp @@ -4,8 +4,7 @@ Grid physics library, www.github.com/paboyle/Grid Source file: extras/Hadrons/ModuleFactory.hpp -Copyright (C) 2015 -Copyright (C) 2016 +Copyright (C) 2015-2018 Author: Antonin Portelli diff --git a/extras/Hadrons/Modules.hpp b/extras/Hadrons/Modules.hpp index cf381d0f..c42e5e88 100644 --- a/extras/Hadrons/Modules.hpp +++ b/extras/Hadrons/Modules.hpp @@ -4,11 +4,10 @@ Grid physics library, www.github.com/paboyle/Grid Source file: extras/Hadrons/Modules.hpp -Copyright (C) 2015 -Copyright (C) 2016 -Copyright (C) 2017 +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 @@ -27,7 +26,6 @@ with this program; if not, write to the Free Software Foundation, Inc., See the full license in the file "LICENSE" in the top level distribution directory *************************************************************************************/ /* END LEGAL */ - #include #include #include @@ -46,10 +44,10 @@ 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 @@ -58,3 +56,9 @@ See the full license in the file "LICENSE" in the top level distribution directo #include #include #include +#include +#include +#include +#include +#include +#include diff --git a/extras/Hadrons/Modules/MAction/DWF.hpp b/extras/Hadrons/Modules/MAction/DWF.hpp index d99f1165..4dfd06cf 100644 --- a/extras/Hadrons/Modules/MAction/DWF.hpp +++ b/extras/Hadrons/Modules/MAction/DWF.hpp @@ -4,10 +4,10 @@ Grid physics library, www.github.com/paboyle/Grid Source file: extras/Hadrons/Modules/MAction/DWF.hpp -Copyright (C) 2015 -Copyright (C) 2016 +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 diff --git a/extras/Hadrons/Modules/MAction/Wilson.hpp b/extras/Hadrons/Modules/MAction/Wilson.hpp index 8ef755bb..6467b3ee 100644 --- a/extras/Hadrons/Modules/MAction/Wilson.hpp +++ b/extras/Hadrons/Modules/MAction/Wilson.hpp @@ -4,10 +4,10 @@ Grid physics library, www.github.com/paboyle/Grid Source file: extras/Hadrons/Modules/MAction/Wilson.hpp -Copyright (C) 2015 -Copyright (C) 2016 +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 diff --git a/extras/Hadrons/Modules/MAction/WilsonClover.hpp b/extras/Hadrons/Modules/MAction/WilsonClover.hpp new file mode 100644 index 00000000..c369f086 --- /dev/null +++ b/extras/Hadrons/Modules/MAction/WilsonClover.hpp @@ -0,0 +1,153 @@ +/************************************************************************************* + +Grid physics library, www.github.com/paboyle/Grid + +Source file: extras/Hadrons/Modules/MAction/Wilson.hpp + +Copyright (C) 2015 +Copyright (C) 2016 + +Author: Antonin Portelli + +This program is free software; you can redistribute it and/or modify +it under the terms of the GNU General Public License as published by +the Free Software Foundation; either version 2 of the License, or +(at your option) any later version. + +This program is distributed in the hope that it will be useful, +but WITHOUT ANY WARRANTY; without even the implied warranty of +MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +GNU General Public License for more details. + +You should have received a copy of the GNU General Public License along +with this program; if not, write to the Free Software Foundation, Inc., +51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. + +See the full license in the file "LICENSE" in the top level distribution directory +*************************************************************************************/ +/* END LEGAL */ + +#ifndef Hadrons_MAction_WilsonClover_hpp_ +#define Hadrons_MAction_WilsonClover_hpp_ + +#include +#include +#include + +BEGIN_HADRONS_NAMESPACE + +/****************************************************************************** + * TWilson quark action * + ******************************************************************************/ +BEGIN_MODULE_NAMESPACE(MAction) + +class WilsonCloverPar: Serializable +{ +public: + GRID_SERIALIZABLE_CLASS_MEMBERS(WilsonCloverPar, + std::string, gauge, + double , mass, + double , csw_r, + double , csw_t, + WilsonAnisotropyCoefficients ,clover_anisotropy, + std::string, boundary + ); +}; + +template +class TWilsonClover: public Module +{ +public: + FGS_TYPE_ALIASES(FImpl,); +public: + // constructor + TWilsonClover(const std::string name); + // destructor + virtual ~TWilsonClover(void) = default; + // dependencies/products + virtual std::vector getInput(void); + virtual std::vector getOutput(void); + // setup + virtual void setup(void); + // execution + virtual void execute(void); +}; + +MODULE_REGISTER_NS(WilsonClover, TWilsonClover, MAction); + +/****************************************************************************** + * TWilsonClover template implementation * + ******************************************************************************/ +// constructor ///////////////////////////////////////////////////////////////// +template +TWilsonClover::TWilsonClover(const std::string name) +: Module(name) +{} + +// dependencies/products /////////////////////////////////////////////////////// +template +std::vector TWilsonClover::getInput(void) +{ + std::vector in = {par().gauge}; + + return in; +} + +template +std::vector TWilsonClover::getOutput(void) +{ + std::vector out = {getName()}; + + return out; +} + +// setup /////////////////////////////////////////////////////////////////////// +template +void TWilsonClover::setup(void) +{ + //unsigned int size; + + // size = 2*env().template lattice4dSize(); + // env().registerObject(getName(), size); + + + LOG(Message) << "Setting up TWilsonClover fermion matrix with m= " << par().mass + << " using gauge field '" << par().gauge << "'" << std::endl; + LOG(Message) << "Fermion boundary conditions: " << par().boundary + << std::endl; + LOG(Message) << "Clover term csw_r: " << par().csw_r + << " csw_t: " << par().csw_t + << std::endl; + auto &U = envGet(LatticeGaugeField, par().gauge); + auto &grid = *env().getGrid(); + auto &gridRb = *env().getRbGrid(); + std::vector boundary = strToVec(par().boundary); + typename WilsonCloverFermion::ImplParams implParams(boundary); + envCreateDerived(FMat, WilsonCloverFermion, getName(), 1, U, grid, gridRb, par().mass, + par().csw_r, + par().csw_t, + par().clover_anisotropy, + implParams); + + + //FMat *fMatPt = new WilsonCloverFermion(U, grid, gridRb, par().mass, + // par().csw_r, + // par().csw_t, + // par().clover_anisotropy, + // implParams); + //env().setObject(getName(), fMatPt); + +} + +// execution /////////////////////////////////////////////////////////////////// +template +void TWilsonClover::execute() +{ + +} + +END_MODULE_NAMESPACE + +END_HADRONS_NAMESPACE + +#endif // Hadrons_WilsonClover_hpp_ diff --git a/extras/Hadrons/Modules/MContraction/Baryon.hpp b/extras/Hadrons/Modules/MContraction/Baryon.hpp index 1ef2e257..8966d95b 100644 --- a/extras/Hadrons/Modules/MContraction/Baryon.hpp +++ b/extras/Hadrons/Modules/MContraction/Baryon.hpp @@ -4,10 +4,10 @@ Grid physics library, www.github.com/paboyle/Grid Source file: extras/Hadrons/Modules/MContraction/Baryon.hpp -Copyright (C) 2015 -Copyright (C) 2016 +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 @@ -122,7 +122,7 @@ void TBaryon::execute(void) << " quarks '" << par().q1 << "', '" << par().q2 << "', and '" << par().q3 << "'" << std::endl; - CorrWriter writer(par().output); + ResultWriter writer(RESULT_FILE_NAME(par().output)); auto &q1 = envGet(PropagatorField1, par().q1); auto &q2 = envGet(PropagatorField2, par().q2); auto &q3 = envGet(PropagatorField3, par().q2); diff --git a/extras/Hadrons/Modules/MContraction/DiscLoop.hpp b/extras/Hadrons/Modules/MContraction/DiscLoop.hpp index ef50061c..539abbbb 100644 --- a/extras/Hadrons/Modules/MContraction/DiscLoop.hpp +++ b/extras/Hadrons/Modules/MContraction/DiscLoop.hpp @@ -4,9 +4,10 @@ Grid physics library, www.github.com/paboyle/Grid Source file: extras/Hadrons/Modules/MContraction/DiscLoop.hpp -Copyright (C) 2017 +Copyright (C) 2015-2018 -Author: Andrew Lawson +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 @@ -118,7 +119,7 @@ void TDiscLoop::execute(void) << "' using '" << par().q_loop << "' with " << par().gamma << " insertion." << std::endl; - CorrWriter writer(par().output); + ResultWriter writer(RESULT_FILE_NAME(par().output)); auto &q_loop = envGet(PropagatorField, par().q_loop); Gamma gamma(par().gamma); std::vector buf; diff --git a/extras/Hadrons/Modules/MContraction/Gamma3pt.hpp b/extras/Hadrons/Modules/MContraction/Gamma3pt.hpp index fb9a9d4b..b4327a13 100644 --- a/extras/Hadrons/Modules/MContraction/Gamma3pt.hpp +++ b/extras/Hadrons/Modules/MContraction/Gamma3pt.hpp @@ -4,9 +4,10 @@ Grid physics library, www.github.com/paboyle/Grid Source file: extras/Hadrons/Modules/MContraction/Gamma3pt.hpp -Copyright (C) 2017 +Copyright (C) 2015-2018 -Author: Andrew Lawson +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 @@ -152,7 +153,7 @@ void TGamma3pt::execute(void) // Initialise variables. q2 and q3 are normal propagators, q1 may be // sink smeared. - CorrWriter writer(par().output); + ResultWriter writer(RESULT_FILE_NAME(par().output)); auto &q1 = envGet(SlicedPropagator1, par().q1); auto &q2 = envGet(PropagatorField2, par().q2); auto &q3 = envGet(PropagatorField2, par().q3); diff --git a/extras/Hadrons/Modules/MContraction/Meson.hpp b/extras/Hadrons/Modules/MContraction/Meson.hpp index 46bbdb2e..b0988a5c 100644 --- a/extras/Hadrons/Modules/MContraction/Meson.hpp +++ b/extras/Hadrons/Modules/MContraction/Meson.hpp @@ -4,12 +4,10 @@ Grid physics library, www.github.com/paboyle/Grid Source file: extras/Hadrons/Modules/MContraction/Meson.hpp -Copyright (C) 2015 -Copyright (C) 2016 -Copyright (C) 2017 +Copyright (C) 2015-2018 Author: Antonin Portelli - Andrew Lawson +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 @@ -173,8 +171,8 @@ void TMeson::execute(void) LOG(Message) << "Computing meson contractions '" << getName() << "' using" << " quarks '" << par().q1 << "' and '" << par().q2 << "'" << std::endl; - - CorrWriter writer(par().output); + + ResultWriter writer(RESULT_FILE_NAME(par().output)); std::vector buf; std::vector result; Gamma g5(Gamma::Algebra::Gamma5); diff --git a/extras/Hadrons/Modules/MContraction/WardIdentity.hpp b/extras/Hadrons/Modules/MContraction/WardIdentity.hpp index 556450a8..2801d88c 100644 --- a/extras/Hadrons/Modules/MContraction/WardIdentity.hpp +++ b/extras/Hadrons/Modules/MContraction/WardIdentity.hpp @@ -4,9 +4,10 @@ Grid physics library, www.github.com/paboyle/Grid Source file: extras/Hadrons/Modules/MContraction/WardIdentity.hpp -Copyright (C) 2017 +Copyright (C) 2015-2018 -Author: Andrew Lawson +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 diff --git a/extras/Hadrons/Modules/MContraction/WeakHamiltonian.hpp b/extras/Hadrons/Modules/MContraction/WeakHamiltonian.hpp index 7df40370..9d8ada98 100644 --- a/extras/Hadrons/Modules/MContraction/WeakHamiltonian.hpp +++ b/extras/Hadrons/Modules/MContraction/WeakHamiltonian.hpp @@ -4,9 +4,10 @@ Grid physics library, www.github.com/paboyle/Grid Source file: extras/Hadrons/Modules/MContraction/WeakHamiltonian.hpp -Copyright (C) 2017 +Copyright (C) 2015-2018 -Author: Andrew Lawson +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 diff --git a/extras/Hadrons/Modules/MContraction/WeakHamiltonianEye.cc b/extras/Hadrons/Modules/MContraction/WeakHamiltonianEye.cc index 43dfa609..1d257fc7 100644 --- a/extras/Hadrons/Modules/MContraction/WeakHamiltonianEye.cc +++ b/extras/Hadrons/Modules/MContraction/WeakHamiltonianEye.cc @@ -4,9 +4,10 @@ Grid physics library, www.github.com/paboyle/Grid Source file: extras/Hadrons/Modules/MContraction/WeakHamiltonianEye.cc -Copyright (C) 2017 +Copyright (C) 2015-2018 -Author: Andrew Lawson +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 @@ -103,7 +104,7 @@ void TWeakHamiltonianEye::execute(void) << par().q2 << ", '" << par().q3 << "' and '" << par().q4 << "'." << std::endl; - CorrWriter writer(par().output); + ResultWriter writer(RESULT_FILE_NAME(par().output)); auto &q1 = envGet(SlicedPropagator, par().q1); auto &q2 = envGet(PropagatorField, par().q2); auto &q3 = envGet(PropagatorField, par().q3); diff --git a/extras/Hadrons/Modules/MContraction/WeakHamiltonianEye.hpp b/extras/Hadrons/Modules/MContraction/WeakHamiltonianEye.hpp index 3a2b9309..24f39f6c 100644 --- a/extras/Hadrons/Modules/MContraction/WeakHamiltonianEye.hpp +++ b/extras/Hadrons/Modules/MContraction/WeakHamiltonianEye.hpp @@ -4,9 +4,10 @@ Grid physics library, www.github.com/paboyle/Grid Source file: extras/Hadrons/Modules/MContraction/WeakHamiltonianEye.hpp -Copyright (C) 2017 +Copyright (C) 2015-2018 -Author: Andrew Lawson +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 diff --git a/extras/Hadrons/Modules/MContraction/WeakHamiltonianNonEye.cc b/extras/Hadrons/Modules/MContraction/WeakHamiltonianNonEye.cc index 8a7113e3..2ad2e7dc 100644 --- a/extras/Hadrons/Modules/MContraction/WeakHamiltonianNonEye.cc +++ b/extras/Hadrons/Modules/MContraction/WeakHamiltonianNonEye.cc @@ -4,9 +4,10 @@ Grid physics library, www.github.com/paboyle/Grid Source file: extras/Hadrons/Modules/MContraction/WeakHamiltonianNonEye.cc -Copyright (C) 2017 +Copyright (C) 2015-2018 -Author: Andrew Lawson +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 @@ -103,7 +104,7 @@ void TWeakHamiltonianNonEye::execute(void) << par().q2 << ", '" << par().q3 << "' and '" << par().q4 << "'." << std::endl; - CorrWriter writer(par().output); + ResultWriter writer(RESULT_FILE_NAME(par().output)); auto &q1 = envGet(PropagatorField, par().q1); auto &q2 = envGet(PropagatorField, par().q2); auto &q3 = envGet(PropagatorField, par().q3); diff --git a/extras/Hadrons/Modules/MContraction/WeakHamiltonianNonEye.hpp b/extras/Hadrons/Modules/MContraction/WeakHamiltonianNonEye.hpp index eb5abe3c..c4cd66f1 100644 --- a/extras/Hadrons/Modules/MContraction/WeakHamiltonianNonEye.hpp +++ b/extras/Hadrons/Modules/MContraction/WeakHamiltonianNonEye.hpp @@ -4,9 +4,10 @@ Grid physics library, www.github.com/paboyle/Grid Source file: extras/Hadrons/Modules/MContraction/WeakHamiltonianNonEye.hpp -Copyright (C) 2017 +Copyright (C) 2015-2018 -Author: Andrew Lawson +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 diff --git a/extras/Hadrons/Modules/MContraction/WeakNeutral4ptDisc.cc b/extras/Hadrons/Modules/MContraction/WeakNeutral4ptDisc.cc index 18423f3e..2c94b2ba 100644 --- a/extras/Hadrons/Modules/MContraction/WeakNeutral4ptDisc.cc +++ b/extras/Hadrons/Modules/MContraction/WeakNeutral4ptDisc.cc @@ -4,9 +4,10 @@ Grid physics library, www.github.com/paboyle/Grid Source file: extras/Hadrons/Modules/MContraction/WeakNeutral4ptDisc.cc -Copyright (C) 2017 +Copyright (C) 2015-2018 -Author: Andrew Lawson +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 @@ -103,7 +104,7 @@ void TWeakNeutral4ptDisc::execute(void) << par().q2 << ", '" << par().q3 << "' and '" << par().q4 << "'." << std::endl; - CorrWriter writer(par().output); + ResultWriter writer(RESULT_FILE_NAME(par().output)); auto &q1 = envGet(PropagatorField, par().q1); auto &q2 = envGet(PropagatorField, par().q2); auto &q3 = envGet(PropagatorField, par().q3); diff --git a/extras/Hadrons/Modules/MContraction/WeakNeutral4ptDisc.hpp b/extras/Hadrons/Modules/MContraction/WeakNeutral4ptDisc.hpp index f26d4636..5de2a751 100644 --- a/extras/Hadrons/Modules/MContraction/WeakNeutral4ptDisc.hpp +++ b/extras/Hadrons/Modules/MContraction/WeakNeutral4ptDisc.hpp @@ -4,9 +4,10 @@ Grid physics library, www.github.com/paboyle/Grid Source file: extras/Hadrons/Modules/MContraction/WeakNeutral4ptDisc.hpp -Copyright (C) 2017 +Copyright (C) 2015-2018 -Author: Andrew Lawson +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 diff --git a/extras/Hadrons/Modules/MFermion/GaugeProp.hpp b/extras/Hadrons/Modules/MFermion/GaugeProp.hpp index 05b3d17a..412e92d5 100644 --- a/extras/Hadrons/Modules/MFermion/GaugeProp.hpp +++ b/extras/Hadrons/Modules/MFermion/GaugeProp.hpp @@ -4,12 +4,10 @@ Grid physics library, www.github.com/paboyle/Grid Source file: extras/Hadrons/Modules/MFermion/GaugeProp.hpp -Copyright (C) 2015 -Copyright (C) 2016 -Copyright (C) 2017 +Copyright (C) 2015-2018 Author: Antonin Portelli - Andrew Lawson +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 @@ -96,7 +94,6 @@ private: }; MODULE_REGISTER_NS(GaugeProp, TGaugeProp, MFermion); - /****************************************************************************** * TGaugeProp implementation * ******************************************************************************/ @@ -156,7 +153,7 @@ void TGaugeProp::execute(void) LOG(Message) << "Inverting using solver '" << par().solver << "' on source '" << par().source << "'" << std::endl; for (unsigned int s = 0; s < Ns; ++s) - for (unsigned int c = 0; c < Nc; ++c) + for (unsigned int c = 0; c < FImpl::Dimension; ++c) { LOG(Message) << "Inversion for spin= " << s << ", color= " << c << std::endl; @@ -165,11 +162,11 @@ void TGaugeProp::execute(void) { if (Ls_ == 1) { - PropToFerm(source, fullSrc, s, c); + PropToFerm(source, fullSrc, s, c); } else { - PropToFerm(tmp, fullSrc, s, c); + PropToFerm(tmp, fullSrc, s, c); make_5D(tmp, source, Ls_); } } @@ -182,18 +179,18 @@ void TGaugeProp::execute(void) } else { - PropToFerm(source, fullSrc, s, c); + PropToFerm(source, fullSrc, s, c); } } sol = zero; solver(sol, source); - FermToProp(prop, sol, s, c); + FermToProp(prop, sol, s, c); // create 4D propagators from 5D one if necessary if (Ls_ > 1) { PropagatorField &p4d = envGet(PropagatorField, getName()); make_4D(sol, tmp, Ls_); - FermToProp(p4d, tmp, s, c); + FermToProp(p4d, tmp, s, c); } } } diff --git a/extras/Hadrons/Modules/MGauge/FundtoHirep.cc b/extras/Hadrons/Modules/MGauge/FundtoHirep.cc new file mode 100644 index 00000000..f15a3b7c --- /dev/null +++ b/extras/Hadrons/Modules/MGauge/FundtoHirep.cc @@ -0,0 +1,75 @@ +/************************************************************************************* + +Grid physics library, www.github.com/paboyle/Grid + +Source file: extras/Hadrons/Modules/MGauge/FundtoHirep.cc + +Copyright (C) 2015 +Copyright (C) 2016 + + +This program is free software; you can redistribute it and/or modify +it under the terms of the GNU General Public License as published by +the Free Software Foundation; either version 2 of the License, or +(at your option) any later version. + +This program is distributed in the hope that it will be useful, +but WITHOUT ANY WARRANTY; without even the implied warranty of +MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +GNU General Public License for more details. + +You should have received a copy of the GNU General Public License along +with this program; if not, write to the Free Software Foundation, Inc., +51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. + +See the full license in the file "LICENSE" in the top level distribution directory +*************************************************************************************/ +/* END LEGAL */ + +#include + +using namespace Grid; +using namespace Hadrons; +using namespace MGauge; + +// constructor ///////////////////////////////////////////////////////////////// +template +TFundtoHirep::TFundtoHirep(const std::string name) +: Module(name) +{} + +// dependencies/products /////////////////////////////////////////////////////// +template +std::vector TFundtoHirep::getInput(void) +{ + std::vector in; + return in; +} + +template +std::vector TFundtoHirep::getOutput(void) +{ + std::vector out = {getName()}; + return out; +} + +// setup /////////////////////////////////////////////////////////////////////// +template +void TFundtoHirep::setup(void) +{ + env().template registerLattice(getName()); +} + +// execution /////////////////////////////////////////////////////////////////// +template +void TFundtoHirep::execute(void) +{ + auto &U = *env().template getObject(par().gaugeconf); + LOG(Message) << "Transforming Representation" << std::endl; + + Rep TargetRepresentation(U._grid); + TargetRepresentation.update_representation(U); + + typename Rep::LatticeField &URep = *env().template createLattice(getName()); + URep = TargetRepresentation.U; +} diff --git a/extras/Hadrons/Modules/MGauge/FundtoHirep.hpp b/extras/Hadrons/Modules/MGauge/FundtoHirep.hpp new file mode 100644 index 00000000..6f072783 --- /dev/null +++ b/extras/Hadrons/Modules/MGauge/FundtoHirep.hpp @@ -0,0 +1,77 @@ +/************************************************************************************* + +Grid physics library, www.github.com/paboyle/Grid + +Source file: extras/Hadrons/Modules/MGauge/FundtoHirep.hpp + +Copyright (C) 2015 +Copyright (C) 2016 + +Author: David Preti + Guido Cossu + +This program is free software; you can redistribute it and/or modify +it under the terms of the GNU General Public License as published by +the Free Software Foundation; either version 2 of the License, or +(at your option) any later version. + +This program is distributed in the hope that it will be useful, +but WITHOUT ANY WARRANTY; without even the implied warranty of +MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +GNU General Public License for more details. + +You should have received a copy of the GNU General Public License along +with this program; if not, write to the Free Software Foundation, Inc., +51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. + +See the full license in the file "LICENSE" in the top level distribution directory +*************************************************************************************/ +/* END LEGAL */ + +#ifndef Hadrons_MGauge_FundtoHirep_hpp_ +#define Hadrons_MGauge_FundtoHirep_hpp_ + +#include +#include +#include + +BEGIN_HADRONS_NAMESPACE + +/****************************************************************************** + * Load a NERSC configuration * + ******************************************************************************/ +BEGIN_MODULE_NAMESPACE(MGauge) + +class FundtoHirepPar: Serializable +{ +public: + GRID_SERIALIZABLE_CLASS_MEMBERS(FundtoHirepPar, + std::string, gaugeconf); +}; + +template +class TFundtoHirep: public Module +{ +public: + // constructor + TFundtoHirep(const std::string name); + // destructor + virtual ~TFundtoHirep(void) = default; + // dependency relation + virtual std::vector getInput(void); + virtual std::vector getOutput(void); + // setup + void setup(void); + // execution + void execute(void); +}; + +//MODULE_REGISTER_NS(FundtoAdjoint, TFundtoHirep, MGauge); +//MODULE_REGISTER_NS(FundtoTwoIndexSym, TFundtoHirep, MGauge); +//MODULE_REGISTER_NS(FundtoTwoIndexAsym, TFundtoHirep, MGauge); + +END_MODULE_NAMESPACE + +END_HADRONS_NAMESPACE + +#endif // Hadrons_MGauge_FundtoHirep_hpp_ diff --git a/extras/Hadrons/Modules/MGauge/Random.cc b/extras/Hadrons/Modules/MGauge/Random.cc index 97afd338..962fc243 100644 --- a/extras/Hadrons/Modules/MGauge/Random.cc +++ b/extras/Hadrons/Modules/MGauge/Random.cc @@ -4,8 +4,7 @@ Grid physics library, www.github.com/paboyle/Grid Source file: extras/Hadrons/Modules/MGauge/Random.cc -Copyright (C) 2015 -Copyright (C) 2016 +Copyright (C) 2015-2018 Author: Antonin Portelli diff --git a/extras/Hadrons/Modules/MGauge/Random.hpp b/extras/Hadrons/Modules/MGauge/Random.hpp index a07130e4..51a08dbb 100644 --- a/extras/Hadrons/Modules/MGauge/Random.hpp +++ b/extras/Hadrons/Modules/MGauge/Random.hpp @@ -4,8 +4,7 @@ Grid physics library, www.github.com/paboyle/Grid Source file: extras/Hadrons/Modules/MGauge/Random.hpp -Copyright (C) 2015 -Copyright (C) 2016 +Copyright (C) 2015-2018 Author: Antonin Portelli diff --git a/extras/Hadrons/Modules/MGauge/StochEm.cc b/extras/Hadrons/Modules/MGauge/StochEm.cc index c5318573..21b7f626 100644 --- a/extras/Hadrons/Modules/MGauge/StochEm.cc +++ b/extras/Hadrons/Modules/MGauge/StochEm.cc @@ -4,9 +4,9 @@ Grid physics library, www.github.com/paboyle/Grid Source file: extras/Hadrons/Modules/MGauge/StochEm.cc -Copyright (C) 2015 -Copyright (C) 2016 +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/MGauge/StochEm.hpp b/extras/Hadrons/Modules/MGauge/StochEm.hpp index bacb5172..87b70880 100644 --- a/extras/Hadrons/Modules/MGauge/StochEm.hpp +++ b/extras/Hadrons/Modules/MGauge/StochEm.hpp @@ -4,9 +4,9 @@ Grid physics library, www.github.com/paboyle/Grid Source file: extras/Hadrons/Modules/MGauge/StochEm.hpp -Copyright (C) 2015 -Copyright (C) 2016 +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/MGauge/Unit.cc b/extras/Hadrons/Modules/MGauge/Unit.cc index 8bee1ecc..38b5f3aa 100644 --- a/extras/Hadrons/Modules/MGauge/Unit.cc +++ b/extras/Hadrons/Modules/MGauge/Unit.cc @@ -4,8 +4,7 @@ Grid physics library, www.github.com/paboyle/Grid Source file: extras/Hadrons/Modules/MGauge/Unit.cc -Copyright (C) 2015 -Copyright (C) 2016 +Copyright (C) 2015-2018 Author: Antonin Portelli diff --git a/extras/Hadrons/Modules/MGauge/Unit.hpp b/extras/Hadrons/Modules/MGauge/Unit.hpp index c1650cc7..d6ce5a6b 100644 --- a/extras/Hadrons/Modules/MGauge/Unit.hpp +++ b/extras/Hadrons/Modules/MGauge/Unit.hpp @@ -4,8 +4,7 @@ Grid physics library, www.github.com/paboyle/Grid Source file: extras/Hadrons/Modules/MGauge/Unit.hpp -Copyright (C) 2015 -Copyright (C) 2016 +Copyright (C) 2015-2018 Author: Antonin Portelli diff --git a/extras/Hadrons/Modules/MIO/LoadBinary.hpp b/extras/Hadrons/Modules/MIO/LoadBinary.hpp new file mode 100644 index 00000000..d9a8b5f8 --- /dev/null +++ b/extras/Hadrons/Modules/MIO/LoadBinary.hpp @@ -0,0 +1,140 @@ +/************************************************************************************* + +Grid physics library, www.github.com/paboyle/Grid + +Source file: extras/Hadrons/Modules/MIO/LoadBinary.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_MIO_LoadBinary_hpp_ +#define Hadrons_MIO_LoadBinary_hpp_ + +#include +#include +#include + +BEGIN_HADRONS_NAMESPACE + +/****************************************************************************** + * Load a binary configurations * + ******************************************************************************/ +BEGIN_MODULE_NAMESPACE(MIO) + +class LoadBinaryPar: Serializable +{ +public: + GRID_SERIALIZABLE_CLASS_MEMBERS(LoadBinaryPar, + std::string, file, + std::string, format); +}; + +template +class TLoadBinary: public Module +{ +public: + typedef typename Impl::Field Field; + typedef typename Impl::Simd Simd; + typedef typename Field::vector_object vobj; + typedef typename vobj::scalar_object sobj; + typedef typename sobj::DoublePrecision sobj_double; + typedef BinarySimpleMunger Munger; +public: + // constructor + TLoadBinary(const std::string name); + // destructor + virtual ~TLoadBinary(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(LoadBinary, TLoadBinary, MIO); +MODULE_REGISTER_NS(LoadBinaryScalarSU2, TLoadBinary>, MIO); +MODULE_REGISTER_NS(LoadBinaryScalarSU3, TLoadBinary>, MIO); +MODULE_REGISTER_NS(LoadBinaryScalarSU4, TLoadBinary>, MIO); +MODULE_REGISTER_NS(LoadBinaryScalarSU5, TLoadBinary>, MIO); +MODULE_REGISTER_NS(LoadBinaryScalarSU6, TLoadBinary>, MIO); + +/****************************************************************************** + * TLoadBinary implementation * + ******************************************************************************/ +// constructor ///////////////////////////////////////////////////////////////// +template +TLoadBinary::TLoadBinary(const std::string name) +: Module(name) +{} + +// dependencies/products /////////////////////////////////////////////////////// +template +std::vector TLoadBinary::getInput(void) +{ + std::vector in; + + return in; +} + +template +std::vector TLoadBinary::getOutput(void) +{ + std::vector out = {getName()}; + + return out; +} + +// setup /////////////////////////////////////////////////////////////////////// +template +void TLoadBinary::setup(void) +{ + envCreateLat(Field, getName()); +} + +// execution /////////////////////////////////////////////////////////////////// +template +void TLoadBinary::execute(void) +{ + Munger munge; + uint32_t nersc_csum, scidac_csuma, scidac_csumb; + auto &U = envGet(Field, getName()); + std::string filename = par().file + "." + + std::to_string(vm().getTrajectory()); + + LOG(Message) << "Loading " << par().format + << " binary configuration from file '" << filename + << "'" << std::endl; + BinaryIO::readLatticeObject(U, filename, munge, 0, + par().format, nersc_csum, + scidac_csuma, scidac_csumb); + LOG(Message) << "Checksums:" << std::endl; + LOG(Message) << " NERSC " << nersc_csum << std::endl; + LOG(Message) << " SciDAC A " << scidac_csuma << std::endl; + LOG(Message) << " SciDAC B " << scidac_csumb << std::endl; +} + +END_MODULE_NAMESPACE + +END_HADRONS_NAMESPACE + +#endif // Hadrons_MIO_LoadBinary_hpp_ diff --git a/extras/Hadrons/Modules/MGauge/Load.cc b/extras/Hadrons/Modules/MIO/LoadNersc.cc similarity index 82% rename from extras/Hadrons/Modules/MGauge/Load.cc rename to extras/Hadrons/Modules/MIO/LoadNersc.cc index b168a010..f20606fc 100644 --- a/extras/Hadrons/Modules/MGauge/Load.cc +++ b/extras/Hadrons/Modules/MIO/LoadNersc.cc @@ -2,10 +2,9 @@ Grid physics library, www.github.com/paboyle/Grid -Source file: extras/Hadrons/Modules/MGauge/Load.cc +Source file: extras/Hadrons/Modules/MIO/LoadNersc.cc -Copyright (C) 2015 -Copyright (C) 2016 +Copyright (C) 2015-2018 Author: Antonin Portelli @@ -26,30 +25,29 @@ with this program; if not, write to the Free Software Foundation, Inc., See the full license in the file "LICENSE" in the top level distribution directory *************************************************************************************/ /* END LEGAL */ - -#include +#include using namespace Grid; using namespace Hadrons; -using namespace MGauge; +using namespace MIO; /****************************************************************************** -* TLoad implementation * +* TLoadNersc implementation * ******************************************************************************/ // constructor ///////////////////////////////////////////////////////////////// -TLoad::TLoad(const std::string name) -: Module(name) +TLoadNersc::TLoadNersc(const std::string name) +: Module(name) {} // dependencies/products /////////////////////////////////////////////////////// -std::vector TLoad::getInput(void) +std::vector TLoadNersc::getInput(void) { std::vector in; return in; } -std::vector TLoad::getOutput(void) +std::vector TLoadNersc::getOutput(void) { std::vector out = {getName()}; @@ -57,13 +55,13 @@ std::vector TLoad::getOutput(void) } // setup /////////////////////////////////////////////////////////////////////// -void TLoad::setup(void) +void TLoadNersc::setup(void) { envCreateLat(LatticeGaugeField, getName()); } // execution /////////////////////////////////////////////////////////////////// -void TLoad::execute(void) +void TLoadNersc::execute(void) { FieldMetaData header; std::string fileName = par().file + "." diff --git a/extras/Hadrons/Modules/MGauge/Load.hpp b/extras/Hadrons/Modules/MIO/LoadNersc.hpp similarity index 76% rename from extras/Hadrons/Modules/MGauge/Load.hpp rename to extras/Hadrons/Modules/MIO/LoadNersc.hpp index a338af79..d6742e1e 100644 --- a/extras/Hadrons/Modules/MGauge/Load.hpp +++ b/extras/Hadrons/Modules/MIO/LoadNersc.hpp @@ -2,10 +2,9 @@ Grid physics library, www.github.com/paboyle/Grid -Source file: extras/Hadrons/Modules/MGauge/Load.hpp +Source file: extras/Hadrons/Modules/MIO/LoadNersc.hpp -Copyright (C) 2015 -Copyright (C) 2016 +Copyright (C) 2015-2018 Author: Antonin Portelli @@ -26,9 +25,8 @@ with this program; if not, write to the Free Software Foundation, Inc., See the full license in the file "LICENSE" in the top level distribution directory *************************************************************************************/ /* END LEGAL */ - -#ifndef Hadrons_MGauge_Load_hpp_ -#define Hadrons_MGauge_Load_hpp_ +#ifndef Hadrons_MIO_LoadNersc_hpp_ +#define Hadrons_MIO_LoadNersc_hpp_ #include #include @@ -37,38 +35,37 @@ See the full license in the file "LICENSE" in the top level distribution directo BEGIN_HADRONS_NAMESPACE /****************************************************************************** - * Load a NERSC configuration * + * Load a NERSC configuration * ******************************************************************************/ -BEGIN_MODULE_NAMESPACE(MGauge) +BEGIN_MODULE_NAMESPACE(MIO) -class LoadPar: Serializable +class LoadNerscPar: Serializable { public: - GRID_SERIALIZABLE_CLASS_MEMBERS(LoadPar, + GRID_SERIALIZABLE_CLASS_MEMBERS(LoadNerscPar, std::string, file); }; -class TLoad: public Module +class TLoadNersc: public Module { public: // constructor - TLoad(const std::string name); + TLoadNersc(const std::string name); // destructor - virtual ~TLoad(void) = default; + virtual ~TLoadNersc(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(Load, TLoad, MGauge); +MODULE_REGISTER_NS(LoadNersc, TLoadNersc, MIO); END_MODULE_NAMESPACE END_HADRONS_NAMESPACE -#endif // Hadrons_MGauge_Load_hpp_ +#endif // Hadrons_MIO_LoadNersc_hpp_ diff --git a/extras/Hadrons/Modules/MLoop/NoiseLoop.hpp b/extras/Hadrons/Modules/MLoop/NoiseLoop.hpp index 512c731a..e61bf163 100644 --- a/extras/Hadrons/Modules/MLoop/NoiseLoop.hpp +++ b/extras/Hadrons/Modules/MLoop/NoiseLoop.hpp @@ -4,9 +4,10 @@ Grid physics library, www.github.com/paboyle/Grid Source file: extras/Hadrons/Modules/MLoop/NoiseLoop.hpp -Copyright (C) 2016 +Copyright (C) 2015-2018 -Author: Andrew Lawson +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 diff --git a/extras/Hadrons/Modules/MScalar/ChargedProp.cc b/extras/Hadrons/Modules/MScalar/ChargedProp.cc index 6cb75a28..1470f1ad 100644 --- a/extras/Hadrons/Modules/MScalar/ChargedProp.cc +++ b/extras/Hadrons/Modules/MScalar/ChargedProp.cc @@ -1,3 +1,31 @@ +/************************************************************************************* + +Grid physics library, www.github.com/paboyle/Grid + +Source file: extras/Hadrons/Modules/MScalar/ChargedProp.cc + +Copyright (C) 2015-2018 + +Author: Antonin Portelli +Author: James Harrison + +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 @@ -105,7 +133,7 @@ void TChargedProp::execute(void) LOG(Message) << "Saving zero-momentum projection to '" << filename << "'..." << std::endl; - CorrWriter writer(filename); + ResultWriter writer(RESULT_FILE_NAME(par().output)); std::vector vecBuf; std::vector result; diff --git a/extras/Hadrons/Modules/MScalar/ChargedProp.hpp b/extras/Hadrons/Modules/MScalar/ChargedProp.hpp index cfcce28e..4d43aec2 100644 --- a/extras/Hadrons/Modules/MScalar/ChargedProp.hpp +++ b/extras/Hadrons/Modules/MScalar/ChargedProp.hpp @@ -1,3 +1,30 @@ +/************************************************************************************* + +Grid physics library, www.github.com/paboyle/Grid + +Source file: extras/Hadrons/Modules/MScalar/ChargedProp.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_MScalar_ChargedProp_hpp_ #define Hadrons_MScalar_ChargedProp_hpp_ diff --git a/extras/Hadrons/Modules/MScalar/FreeProp.cc b/extras/Hadrons/Modules/MScalar/FreeProp.cc index 924db288..ee86b9db 100644 --- a/extras/Hadrons/Modules/MScalar/FreeProp.cc +++ b/extras/Hadrons/Modules/MScalar/FreeProp.cc @@ -1,3 +1,30 @@ +/************************************************************************************* + +Grid physics library, www.github.com/paboyle/Grid + +Source file: extras/Hadrons/Modules/MScalar/FreeProp.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 #include diff --git a/extras/Hadrons/Modules/MScalar/FreeProp.hpp b/extras/Hadrons/Modules/MScalar/FreeProp.hpp index 6b956134..df17f44e 100644 --- a/extras/Hadrons/Modules/MScalar/FreeProp.hpp +++ b/extras/Hadrons/Modules/MScalar/FreeProp.hpp @@ -1,3 +1,30 @@ +/************************************************************************************* + +Grid physics library, www.github.com/paboyle/Grid + +Source file: extras/Hadrons/Modules/MScalar/FreeProp.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_MScalar_FreeProp_hpp_ #define Hadrons_MScalar_FreeProp_hpp_ diff --git a/extras/Hadrons/Modules/MScalar/Scalar.hpp b/extras/Hadrons/Modules/MScalar/Scalar.hpp index db702ff2..7272f1b3 100644 --- a/extras/Hadrons/Modules/MScalar/Scalar.hpp +++ b/extras/Hadrons/Modules/MScalar/Scalar.hpp @@ -1,3 +1,30 @@ +/************************************************************************************* + +Grid physics library, www.github.com/paboyle/Grid + +Source file: extras/Hadrons/Modules/MScalar/Scalar.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_Scalar_hpp_ #define Hadrons_Scalar_hpp_ diff --git a/extras/Hadrons/Modules/MScalarSUN/TrMag.hpp b/extras/Hadrons/Modules/MScalarSUN/TrMag.hpp new file mode 100644 index 00000000..96eb794e --- /dev/null +++ b/extras/Hadrons/Modules/MScalarSUN/TrMag.hpp @@ -0,0 +1,146 @@ +/************************************************************************************* + +Grid physics library, www.github.com/paboyle/Grid + +Source file: extras/Hadrons/Modules/MScalarSUN/TrMag.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_TrMag_hpp_ +#define Hadrons_MScalarSUN_TrMag_hpp_ + +#include +#include +#include + +BEGIN_HADRONS_NAMESPACE + +/****************************************************************************** + * Module to compute tr(mag^n) * + ******************************************************************************/ +BEGIN_MODULE_NAMESPACE(MScalarSUN) + +class TrMagPar: Serializable +{ +public: + GRID_SERIALIZABLE_CLASS_MEMBERS(TrMagPar, + std::string, field, + unsigned int, maxPow, + std::string, output); +}; + +template +class TTrMag: 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, + Real, value); + }; +public: + // constructor + TTrMag(const std::string name); + // destructor + virtual ~TTrMag(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(TrMagSU2, TTrMag>, MScalarSUN); +MODULE_REGISTER_NS(TrMagSU3, TTrMag>, MScalarSUN); +MODULE_REGISTER_NS(TrMagSU4, TTrMag>, MScalarSUN); +MODULE_REGISTER_NS(TrMagSU5, TTrMag>, MScalarSUN); +MODULE_REGISTER_NS(TrMagSU6, TTrMag>, MScalarSUN); + +/****************************************************************************** + * TTrMag implementation * + ******************************************************************************/ +// constructor ///////////////////////////////////////////////////////////////// +template +TTrMag::TTrMag(const std::string name) +: Module(name) +{} + +// dependencies/products /////////////////////////////////////////////////////// +template +std::vector TTrMag::getInput(void) +{ + std::vector in = {par().field}; + + return in; +} + +template +std::vector TTrMag::getOutput(void) +{ + std::vector out = {}; + + return out; +} + +// setup /////////////////////////////////////////////////////////////////////// +template +void TTrMag::setup(void) +{} + +// execution /////////////////////////////////////////////////////////////////// +template +void TTrMag::execute(void) +{ + LOG(Message) << "Computing tr(mag^n) for n even up to " << par().maxPow + << "..." << std::endl; + + std::vector result; + ResultWriter writer(RESULT_FILE_NAME(par().output)); + auto &phi = envGet(Field, par().field); + + auto m2 = sum(phi), mn = m2; + + m2 = -m2*m2; + mn = 1.; + for (unsigned int n = 2; n <= par().maxPow; n += 2) + { + Result r; + + mn = mn*m2; + r.op = "tr(mag^" + std::to_string(n) + ")"; + r.value = TensorRemove(trace(mn)).real(); + result.push_back(r); + } + write(writer, "trmag", result); +} + +END_MODULE_NAMESPACE + +END_HADRONS_NAMESPACE + +#endif // Hadrons_MScalarSUN_TrMag_hpp_ diff --git a/extras/Hadrons/Modules/MScalarSUN/TrPhi.hpp b/extras/Hadrons/Modules/MScalarSUN/TrPhi.hpp new file mode 100644 index 00000000..4586663d --- /dev/null +++ b/extras/Hadrons/Modules/MScalarSUN/TrPhi.hpp @@ -0,0 +1,182 @@ +/************************************************************************************* + +Grid physics library, www.github.com/paboyle/Grid + +Source file: extras/Hadrons/Modules/MScalarSUN/TrPhi.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_TrPhi_hpp_ +#define Hadrons_MScalarSUN_TrPhi_hpp_ + +#include +#include +#include + +BEGIN_HADRONS_NAMESPACE + +/****************************************************************************** + * Module to compute tr(phi^n) * + ******************************************************************************/ +BEGIN_MODULE_NAMESPACE(MScalarSUN) + +class TrPhiPar: Serializable +{ +public: + GRID_SERIALIZABLE_CLASS_MEMBERS(TrPhiPar, + std::string, field, + unsigned int, maxPow, + std::string, output); +}; + +template +class TTrPhi: 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, + Real, value); + }; +public: + // constructor + TTrPhi(const std::string name); + // destructor + virtual ~TTrPhi(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: + // output name generator + std::string outName(const unsigned int n); +}; + +MODULE_REGISTER_NS(TrPhiSU2, TTrPhi>, MScalarSUN); +MODULE_REGISTER_NS(TrPhiSU3, TTrPhi>, MScalarSUN); +MODULE_REGISTER_NS(TrPhiSU4, TTrPhi>, MScalarSUN); +MODULE_REGISTER_NS(TrPhiSU5, TTrPhi>, MScalarSUN); +MODULE_REGISTER_NS(TrPhiSU6, TTrPhi>, MScalarSUN); + +/****************************************************************************** + * TTrPhi implementation * + ******************************************************************************/ +// constructor ///////////////////////////////////////////////////////////////// +template +TTrPhi::TTrPhi(const std::string name) +: Module(name) +{} + +// dependencies/products /////////////////////////////////////////////////////// +template +std::vector TTrPhi::getInput(void) +{ + std::vector in = {par().field}; + + return in; +} + +template +std::vector TTrPhi::getOutput(void) +{ + std::vector out; + + for (unsigned int n = 2; n <= par().maxPow; n += 2) + { + out.push_back(outName(n)); + } + + return out; +} + +// setup /////////////////////////////////////////////////////////////////////// +template +void TTrPhi::setup(void) +{ + if (par().maxPow < 2) + { + HADRON_ERROR(Size, "'maxPow' should be at least equal to 2"); + } + envTmpLat(Field, "phi2"); + envTmpLat(Field, "buf"); + for (unsigned int n = 2; n <= par().maxPow; n += 2) + { + envCreateLat(ComplexField, outName(n)); + } +} + +// execution /////////////////////////////////////////////////////////////////// +template +void TTrPhi::execute(void) +{ + LOG(Message) << "Computing tr(phi^n) for n even up to " << par().maxPow + << "..." << std::endl; + + std::vector result; + auto &phi = envGet(Field, par().field); + + envGetTmp(Field, phi2); + envGetTmp(Field, buf); + buf = 1.; + phi2 = -phi*phi; + for (unsigned int n = 2; n <= par().maxPow; n += 2) + { + auto &phin = envGet(ComplexField, outName(n)); + + buf = buf*phi2; + phin = trace(buf); + if (!par().output.empty()) + { + Result r; + + r.op = "tr(phi^" + std::to_string(n) + ")"; + r.value = TensorRemove(sum(phin)).real(); + result.push_back(r); + } + } + if (result.size() > 0) + { + ResultWriter writer(RESULT_FILE_NAME(par().output)); + + write(writer, "trphi", result); + } +} + +// output name generator /////////////////////////////////////////////////////// +template +std::string TTrPhi::outName(const unsigned int n) +{ + return getName() + "_" + std::to_string(n); +} + +END_MODULE_NAMESPACE + +END_HADRONS_NAMESPACE + +#endif // Hadrons_MScalarSUN_TrPhi_hpp_ diff --git a/extras/Hadrons/Modules/MScalarSUN/TwoPoint.hpp b/extras/Hadrons/Modules/MScalarSUN/TwoPoint.hpp new file mode 100644 index 00000000..abfbf609 --- /dev/null +++ b/extras/Hadrons/Modules/MScalarSUN/TwoPoint.hpp @@ -0,0 +1,184 @@ +/************************************************************************************* + +Grid physics library, www.github.com/paboyle/Grid + +Source file: extras/Hadrons/Modules/MScalarSUN/TwoPoint.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_TwoPoint_hpp_ +#define Hadrons_MScalarSUN_TwoPoint_hpp_ + +#include +#include +#include + +BEGIN_HADRONS_NAMESPACE + +/****************************************************************************** + * 2-pt functions for a given set of operators * + ******************************************************************************/ +BEGIN_MODULE_NAMESPACE(MScalarSUN) + +class TwoPointPar: Serializable +{ +public: + GRID_SERIALIZABLE_CLASS_MEMBERS(TwoPointPar, + std::vector, op, + std::string, output); +}; + +template +class TTwoPoint: public Module +{ +public: + typedef typename SImpl::Field Field; + typedef typename SImpl::ComplexField ComplexField; + class Result: Serializable + { + public: + GRID_SERIALIZABLE_CLASS_MEMBERS(Result, + std::string, sink, + std::string, source, + std::vector, data); + }; +public: + // constructor + TTwoPoint(const std::string name); + // destructor + virtual ~TTwoPoint(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: + // make 2-pt function + template + std::vector makeTwoPoint(const std::vector &sink, + const std::vector &source); +}; + +MODULE_REGISTER_NS(TwoPointSU2, TTwoPoint>, MScalarSUN); +MODULE_REGISTER_NS(TwoPointSU3, TTwoPoint>, MScalarSUN); +MODULE_REGISTER_NS(TwoPointSU4, TTwoPoint>, MScalarSUN); +MODULE_REGISTER_NS(TwoPointSU5, TTwoPoint>, MScalarSUN); +MODULE_REGISTER_NS(TwoPointSU6, TTwoPoint>, MScalarSUN); + +/****************************************************************************** + * TTwoPoint implementation * + ******************************************************************************/ +// constructor ///////////////////////////////////////////////////////////////// +template +TTwoPoint::TTwoPoint(const std::string name) +: Module(name) +{} + +// dependencies/products /////////////////////////////////////////////////////// +template +std::vector TTwoPoint::getInput(void) +{ + return par().op; +} + +template +std::vector TTwoPoint::getOutput(void) +{ + std::vector out = {}; + + return out; +} + +// setup /////////////////////////////////////////////////////////////////////// +template +void TTwoPoint::setup(void) +{ + const unsigned int nt = env().getDim().back(); + envTmp(std::vector>, "slicedOp", 1, par().op.size(), + std::vector(nt)); +} + +// execution /////////////////////////////////////////////////////////////////// +template +void TTwoPoint::execute(void) +{ + LOG(Message) << "Computing 2-point functions for operators:" << std::endl; + for (auto &o: par().op) + { + LOG(Message) << " '" << o << "'" << std::endl; + } + + ResultWriter writer(RESULT_FILE_NAME(par().output)); + const unsigned int nd = env().getDim().size(); + std::vector result; + + envGetTmp(std::vector>, slicedOp); + for (unsigned int i = 0; i < par().op.size(); ++i) + { + auto &op = envGet(ComplexField, par().op[i]); + + sliceSum(op, slicedOp[i], nd - 1); + } + for (unsigned int i = 0; i < par().op.size(); ++i) + for (unsigned int j = 0; j < par().op.size(); ++j) + { + Result r; + + r.sink = par().op[i]; + r.source = par().op[j]; + r.data = makeTwoPoint(slicedOp[i], slicedOp[j]); + result.push_back(r); + } + write(writer, "twopt", result); +} + +// make 2-pt function ////////////////////////////////////////////////////////// +template +template +std::vector TTwoPoint::makeTwoPoint( + const std::vector &sink, + const std::vector &source) +{ + assert(sink.size() == source.size()); + + unsigned int nt = sink.size(); + std::vector res(nt, 0.); + + for (unsigned int dt = 0; dt < nt; ++dt) + { + for (unsigned int t = 0; t < nt; ++t) + { + res[dt] += TensorRemove(trace(sink[(t+dt)%nt]*source[t])); + } + res[dt] *= 1./static_cast(nt); + } + + return res; +} + +END_MODULE_NAMESPACE + +END_HADRONS_NAMESPACE + +#endif // Hadrons_MScalarSUN_TwoPoint_hpp_ diff --git a/extras/Hadrons/Modules/MSink/Point.hpp b/extras/Hadrons/Modules/MSink/Point.hpp index 43be3009..c5f6eff0 100644 --- a/extras/Hadrons/Modules/MSink/Point.hpp +++ b/extras/Hadrons/Modules/MSink/Point.hpp @@ -4,9 +4,10 @@ Grid physics library, www.github.com/paboyle/Grid Source file: extras/Hadrons/Modules/MSink/Point.hpp -Copyright (C) 2017 +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 diff --git a/extras/Hadrons/Modules/MSink/Smear.hpp b/extras/Hadrons/Modules/MSink/Smear.hpp index e85ab263..e72dece0 100644 --- a/extras/Hadrons/Modules/MSink/Smear.hpp +++ b/extras/Hadrons/Modules/MSink/Smear.hpp @@ -4,9 +4,10 @@ Grid physics library, www.github.com/paboyle/Grid Source file: extras/Hadrons/Modules/MSink/Smear.hpp -Copyright (C) 2017 +Copyright (C) 2015-2018 -Author: Andrew Lawson +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 diff --git a/extras/Hadrons/Modules/MSolver/RBPrecCG.hpp b/extras/Hadrons/Modules/MSolver/RBPrecCG.hpp index bb4f3f62..54c0f2d8 100644 --- a/extras/Hadrons/Modules/MSolver/RBPrecCG.hpp +++ b/extras/Hadrons/Modules/MSolver/RBPrecCG.hpp @@ -4,8 +4,7 @@ Grid physics library, www.github.com/paboyle/Grid Source file: extras/Hadrons/Modules/MSolver/RBPrecCG.hpp -Copyright (C) 2015 -Copyright (C) 2016 +Copyright (C) 2015-2018 Author: Antonin Portelli diff --git a/extras/Hadrons/Modules/MSource/Point.hpp b/extras/Hadrons/Modules/MSource/Point.hpp index 1d8241cf..ac6df252 100644 --- a/extras/Hadrons/Modules/MSource/Point.hpp +++ b/extras/Hadrons/Modules/MSource/Point.hpp @@ -4,10 +4,10 @@ Grid physics library, www.github.com/paboyle/Grid Source file: extras/Hadrons/Modules/MSource/Point.hpp -Copyright (C) 2015 -Copyright (C) 2016 +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 diff --git a/extras/Hadrons/Modules/MSource/SeqConserved.hpp b/extras/Hadrons/Modules/MSource/SeqConserved.hpp index 3e8ef457..ee8d8d56 100644 --- a/extras/Hadrons/Modules/MSource/SeqConserved.hpp +++ b/extras/Hadrons/Modules/MSource/SeqConserved.hpp @@ -2,11 +2,12 @@ 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) 2017 +Copyright (C) 2015-2018 -Author: Andrew Lawson +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 diff --git a/extras/Hadrons/Modules/MSource/SeqGamma.hpp b/extras/Hadrons/Modules/MSource/SeqGamma.hpp index abad5ace..40eda29f 100644 --- a/extras/Hadrons/Modules/MSource/SeqGamma.hpp +++ b/extras/Hadrons/Modules/MSource/SeqGamma.hpp @@ -4,11 +4,10 @@ Grid physics library, www.github.com/paboyle/Grid Source file: extras/Hadrons/Modules/MSource/SeqGamma.hpp -Copyright (C) 2015 -Copyright (C) 2016 -Copyright (C) 2017 +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 diff --git a/extras/Hadrons/Modules/MSource/Wall.hpp b/extras/Hadrons/Modules/MSource/Wall.hpp index 9d5f1f46..5853b11a 100644 --- a/extras/Hadrons/Modules/MSource/Wall.hpp +++ b/extras/Hadrons/Modules/MSource/Wall.hpp @@ -4,9 +4,10 @@ Grid physics library, www.github.com/paboyle/Grid Source file: extras/Hadrons/Modules/MSource/Wall.hpp -Copyright (C) 2017 +Copyright (C) 2015-2018 -Author: Andrew Lawson +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 diff --git a/extras/Hadrons/Modules/MSource/Z2.hpp b/extras/Hadrons/Modules/MSource/Z2.hpp index 3593cb34..4414e37f 100644 --- a/extras/Hadrons/Modules/MSource/Z2.hpp +++ b/extras/Hadrons/Modules/MSource/Z2.hpp @@ -4,8 +4,7 @@ Grid physics library, www.github.com/paboyle/Grid Source file: extras/Hadrons/Modules/MSource/Z2.hpp -Copyright (C) 2015 -Copyright (C) 2016 +Copyright (C) 2015-2018 Author: Antonin Portelli diff --git a/extras/Hadrons/Modules/MUtilities/TestSeqConserved.hpp b/extras/Hadrons/Modules/MUtilities/TestSeqConserved.hpp index 0647884c..6ee1e3c2 100644 --- a/extras/Hadrons/Modules/MUtilities/TestSeqConserved.hpp +++ b/extras/Hadrons/Modules/MUtilities/TestSeqConserved.hpp @@ -4,9 +4,10 @@ Grid physics library, www.github.com/paboyle/Grid Source file: extras/Hadrons/Modules/MUtilities/TestSeqConserved.hpp -Copyright (C) 2017 +Copyright (C) 2015-2018 -Author: Andrew Lawson +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 diff --git a/extras/Hadrons/Modules/MUtilities/TestSeqGamma.hpp b/extras/Hadrons/Modules/MUtilities/TestSeqGamma.hpp index fd53eab8..df35d887 100644 --- a/extras/Hadrons/Modules/MUtilities/TestSeqGamma.hpp +++ b/extras/Hadrons/Modules/MUtilities/TestSeqGamma.hpp @@ -4,9 +4,10 @@ Grid physics library, www.github.com/paboyle/Grid Source file: extras/Hadrons/Modules/MUtilities/TestSeqGamma.hpp -Copyright (C) 2017 +Copyright (C) 2015-2018 -Author: Andrew Lawson +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 diff --git a/extras/Hadrons/VirtualMachine.cc b/extras/Hadrons/VirtualMachine.cc index e0035bc1..d47bafb7 100644 --- a/extras/Hadrons/VirtualMachine.cc +++ b/extras/Hadrons/VirtualMachine.cc @@ -4,7 +4,7 @@ Grid physics library, www.github.com/paboyle/Grid Source file: extras/Hadrons/VirtualMachine.cc -Copyright (C) 2017 +Copyright (C) 2015-2018 Author: Antonin Portelli diff --git a/extras/Hadrons/VirtualMachine.hpp b/extras/Hadrons/VirtualMachine.hpp index 3af7d914..19a74f94 100644 --- a/extras/Hadrons/VirtualMachine.hpp +++ b/extras/Hadrons/VirtualMachine.hpp @@ -4,7 +4,7 @@ Grid physics library, www.github.com/paboyle/Grid Source file: extras/Hadrons/VirtualMachine.hpp -Copyright (C) 2017 +Copyright (C) 2015-2018 Author: Antonin Portelli diff --git a/extras/Hadrons/make_module_list.sh b/extras/Hadrons/make_module_list.sh index 8c6fa4da..ddc56ff6 100755 --- a/extras/Hadrons/make_module_list.sh +++ b/extras/Hadrons/make_module_list.sh @@ -7,36 +7,6 @@ echo 'modules_hpp =\' >> modules.inc find Modules -name '*.hpp' -type f -print | sed 's/^/ /;$q;s/$/ \\/' >> modules.inc echo '' >> modules.inc rm -f Modules.hpp -echo "/************************************************************************************* - -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 */ -" > Modules.hpp for f in `find Modules -name '*.hpp'`; do echo "#include " >> Modules.hpp done diff --git a/extras/Hadrons/modules.inc b/extras/Hadrons/modules.inc index 199bb5cd..3e62c91f 100644 --- a/extras/Hadrons/modules.inc +++ b/extras/Hadrons/modules.inc @@ -1,13 +1,16 @@ modules_cc =\ + Modules/MScalar/ChargedProp.cc \ + Modules/MScalar/FreeProp.cc \ Modules/MContraction/WeakHamiltonianEye.cc \ Modules/MContraction/WeakNeutral4ptDisc.cc \ Modules/MContraction/WeakHamiltonianNonEye.cc \ - Modules/MGauge/Load.cc \ Modules/MGauge/Unit.cc \ Modules/MGauge/StochEm.cc \ Modules/MGauge/Random.cc \ + Modules/MGauge/FundtoHirep.cc \ Modules/MScalar/FreeProp.cc \ - Modules/MScalar/ChargedProp.cc + Modules/MScalar/ChargedProp.cc \ + Modules/MIO/LoadNersc.cc modules_hpp =\ Modules/MContraction/Baryon.hpp \ @@ -28,10 +31,10 @@ modules_hpp =\ Modules/MSink/Smear.hpp \ Modules/MSink/Point.hpp \ Modules/MSolver/RBPrecCG.hpp \ - Modules/MGauge/Load.hpp \ Modules/MGauge/Unit.hpp \ Modules/MGauge/Random.hpp \ Modules/MGauge/StochEm.hpp \ + Modules/MGauge/FundtoHirep.hpp \ Modules/MUtilities/TestSeqGamma.hpp \ Modules/MUtilities/TestSeqConserved.hpp \ Modules/MLoop/NoiseLoop.hpp \ @@ -39,5 +42,11 @@ modules_hpp =\ Modules/MScalar/Scalar.hpp \ Modules/MScalar/ChargedProp.hpp \ Modules/MAction/DWF.hpp \ - Modules/MAction/Wilson.hpp + Modules/MAction/Wilson.hpp \ + Modules/MAction/WilsonClover.hpp \ + Modules/MScalarSUN/TrMag.hpp \ + Modules/MScalarSUN/TwoPoint.hpp \ + Modules/MScalarSUN/TrPhi.hpp \ + Modules/MIO/LoadNersc.hpp \ + Modules/MIO/LoadBinary.hpp diff --git a/lib/algorithms/LinearOperator.h b/lib/algorithms/LinearOperator.h index 26746e6e..38444c04 100644 --- a/lib/algorithms/LinearOperator.h +++ b/lib/algorithms/LinearOperator.h @@ -183,11 +183,13 @@ namespace Grid { virtual RealD Mpc (const Field &in, Field &out) =0; virtual RealD MpcDag (const Field &in, Field &out) =0; virtual void MpcDagMpc(const Field &in, Field &out,RealD &ni,RealD &no) { - Field tmp(in._grid); + Field tmp(in._grid); + tmp.checkerboard = in.checkerboard; ni=Mpc(in,tmp); no=MpcDag(tmp,out); } virtual void HermOpAndNorm(const Field &in, Field &out,RealD &n1,RealD &n2){ + out.checkerboard = in.checkerboard; MpcDagMpc(in,out,n1,n2); } virtual void HermOp(const Field &in, Field &out){ @@ -215,13 +217,15 @@ namespace Grid { public: SchurDiagMooeeOperator (Matrix &Mat): _Mat(Mat){}; virtual RealD Mpc (const Field &in, Field &out) { - Field tmp(in._grid); -// std::cout <<"grid pointers: in._grid="<< in._grid << " out._grid=" << out._grid << " _Mat.Grid=" << _Mat.Grid() << " _Mat.RedBlackGrid=" << _Mat.RedBlackGrid() << std::endl; + Field tmp(in._grid); + tmp.checkerboard = !in.checkerboard; + //std::cout <<"grid pointers: in._grid="<< in._grid << " out._grid=" << out._grid << " _Mat.Grid=" << _Mat.Grid() << " _Mat.RedBlackGrid=" << _Mat.RedBlackGrid() << std::endl; _Mat.Meooe(in,tmp); _Mat.MooeeInv(tmp,out); _Mat.Meooe(out,tmp); + //std::cout << "cb in " << in.checkerboard << " cb out " << out.checkerboard << std::endl; _Mat.Mooee(in,out); return axpy_norm(out,-1.0,tmp,out); } diff --git a/lib/cartesian/Cartesian_base.h b/lib/cartesian/Cartesian_base.h index acc870de..05a8a3da 100644 --- a/lib/cartesian/Cartesian_base.h +++ b/lib/cartesian/Cartesian_base.h @@ -59,6 +59,7 @@ public: virtual ~GridBase() = default; + // Physics Grid information. std::vector _simd_layout;// Which dimensions get relayed out over simd lanes. std::vector _fdimensions;// (full) Global dimensions of array prior to cb removal @@ -78,6 +79,8 @@ public: std::vector _lstart; // local start of array in gcoors _processor_coor[d]*_ldimensions[d] std::vector _lend ; // local end of array in gcoors _processor_coor[d]*_ldimensions[d]+_ldimensions_[d]-1 + bool _isCheckerBoarded; + public: //////////////////////////////////////////////////////////////// diff --git a/lib/cartesian/Cartesian_full.h b/lib/cartesian/Cartesian_full.h index 9273abf3..b6297d3d 100644 --- a/lib/cartesian/Cartesian_full.h +++ b/lib/cartesian/Cartesian_full.h @@ -97,6 +97,7 @@ public: /////////////////////// // Grid information /////////////////////// + _isCheckerBoarded = false; _ndimension = dimensions.size(); _fdimensions.resize(_ndimension); @@ -122,6 +123,7 @@ public: // Use a reduced simd grid _ldimensions[d] = _gdimensions[d] / _processors[d]; //local dimensions + //std::cout << _ldimensions[d] << " " << _gdimensions[d] << " " << _processors[d] << std::endl; assert(_ldimensions[d] * _processors[d] == _gdimensions[d]); _rdimensions[d] = _ldimensions[d] / _simd_layout[d]; //overdecomposition @@ -166,6 +168,7 @@ public: block = block * _rdimensions[d]; } }; + }; } #endif diff --git a/lib/cartesian/Cartesian_red_black.h b/lib/cartesian/Cartesian_red_black.h index ee424385..5a041f65 100644 --- a/lib/cartesian/Cartesian_red_black.h +++ b/lib/cartesian/Cartesian_red_black.h @@ -171,9 +171,8 @@ public: const std::vector &checker_dim_mask, int checker_dim) { - /////////////////////// - // Grid information - /////////////////////// + + _isCheckerBoarded = true; _checker_dim = checker_dim; assert(checker_dim_mask[checker_dim] == 1); _ndimension = dimensions.size(); diff --git a/lib/communicator/Communicator_base.h b/lib/communicator/Communicator_base.h index 359846c9..7d6911d3 100644 --- a/lib/communicator/Communicator_base.h +++ b/lib/communicator/Communicator_base.h @@ -82,6 +82,7 @@ public: void InitFromMPICommunicator(const std::vector &processors, Grid_MPI_Comm communicator_base); public: + //////////////////////////////////////////////////////////////////////////////////////// // Wraps MPI_Cart routines, or implements equivalent on other impls diff --git a/lib/log/Log.h b/lib/log/Log.h index ddff4c1d..011a7250 100644 --- a/lib/log/Log.h +++ b/lib/log/Log.h @@ -86,6 +86,7 @@ protected: Colours &Painter; int active; int timing_mode; + int topWidth{-1}; static int timestamp; std::string name, topName; std::string COLOUR; @@ -124,11 +125,17 @@ public: Reset(); } } + void setTopWidth(const int w) {topWidth = w;} friend std::ostream& operator<< (std::ostream& stream, Logger& log){ if ( log.active ) { - stream << log.background()<< std::left << log.topName << log.background()<< " : "; + stream << log.background()<< std::left; + if (log.topWidth > 0) + { + stream << std::setw(log.topWidth); + } + stream << log.topName << log.background()<< " : "; stream << log.colour() << std::left << log.name << log.background() << " : "; if ( log.timestamp ) { log.StopWatch->Stop(); diff --git a/lib/qcd/QCD.h b/lib/qcd/QCD.h index 531b71bd..e1e38e82 100644 --- a/lib/qcd/QCD.h +++ b/lib/qcd/QCD.h @@ -39,6 +39,7 @@ namespace QCD { static const int Zdir = 2; static const int Tdir = 3; + static const int Xp = 0; static const int Yp = 1; static const int Zp = 2; @@ -420,15 +421,16 @@ namespace QCD { ////////////////////////////////////////////// // Fermion <-> propagator assignements ////////////////////////////////////////////// - template - void FermToProp(Prop &p, const Ferm &f, const int s, const int c) + //template + template + void FermToProp(typename Fimpl::PropagatorField &p, const typename Fimpl::FermionField &f, const int s, const int c) { - for(int j = 0; j < Ns; ++j) + for(int j = 0; j < Ns; ++j) { auto pjs = peekSpin(p, j, s); auto fj = peekSpin(f, j); - for(int i = 0; i < Nc; ++i) + for(int i = 0; i < Fimpl::Dimension; ++i) { pokeColour(pjs, peekColour(fj, i), i, c); } @@ -436,15 +438,16 @@ namespace QCD { } } - template - void PropToFerm(Ferm &f, const Prop &p, const int s, const int c) + //template + template + void PropToFerm(typename Fimpl::FermionField &f, const typename Fimpl::PropagatorField &p, const int s, const int c) { for(int j = 0; j < Ns; ++j) { auto pjs = peekSpin(p, j, s); auto fj = peekSpin(f, j); - for(int i = 0; i < Nc; ++i) + for(int i = 0; i < Fimpl::Dimension; ++i) { pokeColour(fj, peekColour(pjs, i, c), i); } @@ -503,38 +506,6 @@ namespace QCD { } //namespace QCD } // Grid -/* -<<<<<<< HEAD -#include -#include -#include -#include -#include - -// Include representations -#include -#include -#include -#include - -// Scalar field -#include - -#include - -#include - -#include -#include -#include -#include - - -//#include -======= - ->>>>>>> develop -*/ #endif diff --git a/lib/qcd/action/fermion/Fermion.h b/lib/qcd/action/fermion/Fermion.h index ad2f383d..2a008cb7 100644 --- a/lib/qcd/action/fermion/Fermion.h +++ b/lib/qcd/action/fermion/Fermion.h @@ -50,11 +50,13 @@ Author: Peter Boyle //////////////////////////////////////////// #include // 4d wilson like -#include // 4d wilson like +#include // 4d wilson like +#include // 4d wilson clover fermions #include // 5d base used by all 5d overlap types -//#include + #include #include + #include // Cayley types #include #include @@ -104,10 +106,33 @@ typedef WilsonFermion WilsonTwoIndexSymmetricFermi typedef WilsonFermion WilsonTwoIndexSymmetricFermionF; typedef WilsonFermion WilsonTwoIndexSymmetricFermionD; +typedef WilsonFermion WilsonTwoIndexAntiSymmetricFermionR; +typedef WilsonFermion WilsonTwoIndexAntiSymmetricFermionF; +typedef WilsonFermion WilsonTwoIndexAntiSymmetricFermionD; + +// Twisted mass fermion typedef WilsonTMFermion WilsonTMFermionR; typedef WilsonTMFermion WilsonTMFermionF; typedef WilsonTMFermion WilsonTMFermionD; +// Clover fermions +typedef WilsonCloverFermion WilsonCloverFermionR; +typedef WilsonCloverFermion WilsonCloverFermionF; +typedef WilsonCloverFermion WilsonCloverFermionD; + +typedef WilsonCloverFermion WilsonCloverAdjFermionR; +typedef WilsonCloverFermion WilsonCloverAdjFermionF; +typedef WilsonCloverFermion WilsonCloverAdjFermionD; + +typedef WilsonCloverFermion WilsonCloverTwoIndexSymmetricFermionR; +typedef WilsonCloverFermion WilsonCloverTwoIndexSymmetricFermionF; +typedef WilsonCloverFermion WilsonCloverTwoIndexSymmetricFermionD; + +typedef WilsonCloverFermion WilsonCloverTwoIndexAntiSymmetricFermionR; +typedef WilsonCloverFermion WilsonCloverTwoIndexAntiSymmetricFermionF; +typedef WilsonCloverFermion WilsonCloverTwoIndexAntiSymmetricFermionD; + +// Domain Wall fermions typedef DomainWallFermion DomainWallFermionR; typedef DomainWallFermion DomainWallFermionF; typedef DomainWallFermion DomainWallFermionD; diff --git a/lib/qcd/action/fermion/FermionCore.h b/lib/qcd/action/fermion/FermionCore.h index 17006961..60632c3a 100644 --- a/lib/qcd/action/fermion/FermionCore.h +++ b/lib/qcd/action/fermion/FermionCore.h @@ -70,7 +70,9 @@ Author: Peter Boyle #define TwoIndexFermOpTemplateInstantiate(A) \ template class A; \ - template class A; + template class A; \ + template class A; \ + template class A; #define FermOp5dVecTemplateInstantiate(A) \ template class A; \ diff --git a/lib/qcd/action/fermion/FermionOperatorImpl.h b/lib/qcd/action/fermion/FermionOperatorImpl.h index a9510222..2b900569 100644 --- a/lib/qcd/action/fermion/FermionOperatorImpl.h +++ b/lib/qcd/action/fermion/FermionOperatorImpl.h @@ -261,8 +261,22 @@ namespace QCD { GaugeLinkField link(mat._grid); link = TraceIndex(outerProduct(Btilde,A)); PokeIndex(mat,link,mu); - } + } + + inline void outerProductImpl(PropagatorField &mat, const FermionField &B, const FermionField &A){ + mat = outerProduct(B,A); + } + + inline void TraceSpinImpl(GaugeLinkField &mat, PropagatorField&P) { + mat = TraceIndex(P); + } + inline void extractLinkField(std::vector &mat, DoubledGaugeField &Uds){ + for (int mu = 0; mu < Nd; mu++) + mat[mu] = PeekIndex(Uds, mu); + } + + inline void InsertForce5D(GaugeField &mat, FermionField &Btilde, FermionField Ã,int mu){ int Ls=Btilde._grid->_fdimensions[0]; @@ -393,6 +407,19 @@ class DomainWallVec5dImpl : public PeriodicGaugeImpl< GaugeImplTypes< S,Nrepres assert(0); } + inline void outerProductImpl(PropagatorField &mat, const FermionField &Btilde, const FermionField &A){ + assert(0); + } + + inline void TraceSpinImpl(GaugeLinkField &mat, PropagatorField&P) { + assert(0); + } + + inline void extractLinkField(std::vector &mat, DoubledGaugeField &Uds){ + assert(0); + } + + inline void InsertForce5D(GaugeField &mat, FermionField &Btilde, FermionField Ã, int mu) { assert(0); @@ -636,6 +663,25 @@ class GparityWilsonImpl : public ConjugateGaugeImpl(P); + parallel_for(auto ss = tmp.begin(); ss < tmp.end(); ss++) { + mat[ss]() = tmp[ss](0, 0) + conjugate(tmp[ss](1, 1)); + } + */ + } + + inline void extractLinkField(std::vector &mat, DoubledGaugeField &Uds){ + assert(0); + } + inline void InsertForce5D(GaugeField &mat, FermionField &Btilde, FermionField Ã, int mu) { int Ls = Btilde._grid->_fdimensions[0]; @@ -776,8 +822,8 @@ class StaggeredImpl : public PeriodicGaugeImpl(outerProduct(Btilde,A)); PokeIndex(mat,link,mu); - } - + } + inline void InsertForce5D(GaugeField &mat, FermionField &Btilde, FermionField Ã,int mu){ assert (0); // Must never hit @@ -983,6 +1029,10 @@ typedef WilsonImpl Wilso typedef WilsonImpl WilsonTwoIndexSymmetricImplF; // Float typedef WilsonImpl WilsonTwoIndexSymmetricImplD; // Double +typedef WilsonImpl WilsonTwoIndexAntiSymmetricImplR; // Real.. whichever prec +typedef WilsonImpl WilsonTwoIndexAntiSymmetricImplF; // Float +typedef WilsonImpl WilsonTwoIndexAntiSymmetricImplD; // Double + typedef DomainWallVec5dImpl DomainWallVec5dImplR; // Real.. whichever prec typedef DomainWallVec5dImpl DomainWallVec5dImplF; // Float typedef DomainWallVec5dImpl DomainWallVec5dImplD; // Double diff --git a/lib/qcd/action/fermion/WilsonCloverFermion.cc b/lib/qcd/action/fermion/WilsonCloverFermion.cc new file mode 100644 index 00000000..3c082446 --- /dev/null +++ b/lib/qcd/action/fermion/WilsonCloverFermion.cc @@ -0,0 +1,243 @@ +/************************************************************************************* + + Grid physics library, www.github.com/paboyle/Grid + + Source file: ./lib/qcd/action/fermion/WilsonCloverFermion.cc + + Copyright (C) 2017 + + Author: paboyle + Author: Guido Cossu + + 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 + +namespace Grid +{ +namespace QCD +{ + +// *NOT* EO +template +RealD WilsonCloverFermion::M(const FermionField &in, FermionField &out) +{ + FermionField temp(out._grid); + + // Wilson term + out.checkerboard = in.checkerboard; + this->Dhop(in, out, DaggerNo); + + // Clover term + Mooee(in, temp); + + out += temp; + return norm2(out); +} + +template +RealD WilsonCloverFermion::Mdag(const FermionField &in, FermionField &out) +{ + FermionField temp(out._grid); + + // Wilson term + out.checkerboard = in.checkerboard; + this->Dhop(in, out, DaggerYes); + + // Clover term + MooeeDag(in, temp); + + out += temp; + return norm2(out); +} + +template +void WilsonCloverFermion::ImportGauge(const GaugeField &_Umu) +{ + WilsonFermion::ImportGauge(_Umu); + GridBase *grid = _Umu._grid; + typename Impl::GaugeLinkField Bx(grid), By(grid), Bz(grid), Ex(grid), Ey(grid), Ez(grid); + + // Compute the field strength terms mu>nu + WilsonLoops::FieldStrength(Bx, _Umu, Zdir, Ydir); + WilsonLoops::FieldStrength(By, _Umu, Zdir, Xdir); + WilsonLoops::FieldStrength(Bz, _Umu, Ydir, Xdir); + WilsonLoops::FieldStrength(Ex, _Umu, Tdir, Xdir); + WilsonLoops::FieldStrength(Ey, _Umu, Tdir, Ydir); + WilsonLoops::FieldStrength(Ez, _Umu, Tdir, Zdir); + + // Compute the Clover Operator acting on Colour and Spin + // multiply here by the clover coefficients for the anisotropy + CloverTerm = fillCloverYZ(Bx) * csw_r; + CloverTerm += fillCloverXZ(By) * csw_r; + CloverTerm += fillCloverXY(Bz) * csw_r; + CloverTerm += fillCloverXT(Ex) * csw_t; + CloverTerm += fillCloverYT(Ey) * csw_t; + CloverTerm += fillCloverZT(Ez) * csw_t; + CloverTerm += diag_mass; + + int lvol = _Umu._grid->lSites(); + int DimRep = Impl::Dimension; + + Eigen::MatrixXcd EigenCloverOp = Eigen::MatrixXcd::Zero(Ns * DimRep, Ns * DimRep); + Eigen::MatrixXcd EigenInvCloverOp = Eigen::MatrixXcd::Zero(Ns * DimRep, Ns * DimRep); + + std::vector lcoor; + typename SiteCloverType::scalar_object Qx = zero, Qxinv = zero; + + for (int site = 0; site < lvol; site++) + { + grid->LocalIndexToLocalCoor(site, lcoor); + EigenCloverOp = Eigen::MatrixXcd::Zero(Ns * DimRep, Ns * DimRep); + peekLocalSite(Qx, CloverTerm, lcoor); + Qxinv = zero; + //if (csw!=0){ + for (int j = 0; j < Ns; j++) + for (int k = 0; k < Ns; k++) + for (int a = 0; a < DimRep; a++) + for (int b = 0; b < DimRep; b++) + EigenCloverOp(a + j * DimRep, b + k * DimRep) = Qx()(j, k)(a, b); + // if (site==0) std::cout << "site =" << site << "\n" << EigenCloverOp << std::endl; + + EigenInvCloverOp = EigenCloverOp.inverse(); + //std::cout << EigenInvCloverOp << std::endl; + for (int j = 0; j < Ns; j++) + for (int k = 0; k < Ns; k++) + for (int a = 0; a < DimRep; a++) + for (int b = 0; b < DimRep; b++) + Qxinv()(j, k)(a, b) = EigenInvCloverOp(a + j * DimRep, b + k * DimRep); + // if (site==0) std::cout << "site =" << site << "\n" << EigenInvCloverOp << std::endl; + // } + pokeLocalSite(Qxinv, CloverTermInv, lcoor); + } + + // Separate the even and odd parts + pickCheckerboard(Even, CloverTermEven, CloverTerm); + pickCheckerboard(Odd, CloverTermOdd, CloverTerm); + + pickCheckerboard(Even, CloverTermDagEven, adj(CloverTerm)); + pickCheckerboard(Odd, CloverTermDagOdd, adj(CloverTerm)); + + pickCheckerboard(Even, CloverTermInvEven, CloverTermInv); + pickCheckerboard(Odd, CloverTermInvOdd, CloverTermInv); + + pickCheckerboard(Even, CloverTermInvDagEven, adj(CloverTermInv)); + pickCheckerboard(Odd, CloverTermInvDagOdd, adj(CloverTermInv)); +} + +template +void WilsonCloverFermion::Mooee(const FermionField &in, FermionField &out) +{ + this->MooeeInternal(in, out, DaggerNo, InverseNo); +} + +template +void WilsonCloverFermion::MooeeDag(const FermionField &in, FermionField &out) +{ + this->MooeeInternal(in, out, DaggerYes, InverseNo); +} + +template +void WilsonCloverFermion::MooeeInv(const FermionField &in, FermionField &out) +{ + this->MooeeInternal(in, out, DaggerNo, InverseYes); +} + +template +void WilsonCloverFermion::MooeeInvDag(const FermionField &in, FermionField &out) +{ + this->MooeeInternal(in, out, DaggerYes, InverseYes); +} + +template +void WilsonCloverFermion::MooeeInternal(const FermionField &in, FermionField &out, int dag, int inv) +{ + out.checkerboard = in.checkerboard; + CloverFieldType *Clover; + assert(in.checkerboard == Odd || in.checkerboard == Even); + + if (dag) + { + if (in._grid->_isCheckerBoarded) + { + if (in.checkerboard == Odd) + { + Clover = (inv) ? &CloverTermInvDagOdd : &CloverTermDagOdd; + } + else + { + Clover = (inv) ? &CloverTermInvDagEven : &CloverTermDagEven; + } + out = *Clover * in; + } + else + { + Clover = (inv) ? &CloverTermInv : &CloverTerm; + out = adj(*Clover) * in; + } + } + else + { + if (in._grid->_isCheckerBoarded) + { + + if (in.checkerboard == Odd) + { + // std::cout << "Calling clover term Odd" << std::endl; + Clover = (inv) ? &CloverTermInvOdd : &CloverTermOdd; + } + else + { + // std::cout << "Calling clover term Even" << std::endl; + Clover = (inv) ? &CloverTermInvEven : &CloverTermEven; + } + out = *Clover * in; + // std::cout << GridLogMessage << "*Clover.checkerboard " << (*Clover).checkerboard << std::endl; + } + else + { + Clover = (inv) ? &CloverTermInv : &CloverTerm; + out = *Clover * in; + } + } + +} // MooeeInternal + + +// Derivative parts +template +void WilsonCloverFermion::MooDeriv(GaugeField &mat, const FermionField &X, const FermionField &Y, int dag) +{ + assert(0); +} + +// Derivative parts +template +void WilsonCloverFermion::MeeDeriv(GaugeField &mat, const FermionField &U, const FermionField &V, int dag) +{ + assert(0); // not implemented yet +} + +FermOpTemplateInstantiate(WilsonCloverFermion); +AdjointFermOpTemplateInstantiate(WilsonCloverFermion); +TwoIndexFermOpTemplateInstantiate(WilsonCloverFermion); +//GparityFermOpTemplateInstantiate(WilsonCloverFermion); +} +} diff --git a/lib/qcd/action/fermion/WilsonCloverFermion.h b/lib/qcd/action/fermion/WilsonCloverFermion.h new file mode 100644 index 00000000..268564c0 --- /dev/null +++ b/lib/qcd/action/fermion/WilsonCloverFermion.h @@ -0,0 +1,366 @@ +/************************************************************************************* + + Grid physics library, www.github.com/paboyle/Grid + + Source file: ./lib/qcd/action/fermion/WilsonCloverFermion.h + + Copyright (C) 2017 + + Author: Guido Cossu + Author: David Preti <> + + 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 GRID_QCD_WILSON_CLOVER_FERMION_H +#define GRID_QCD_WILSON_CLOVER_FERMION_H + +#include + +namespace Grid +{ +namespace QCD +{ + +/////////////////////////////////////////////////////////////////// +// Wilson Clover +// +// Operator ( with anisotropy coefficients): +// +// Q = 1 + (Nd-1)/xi_0 + m +// + W_t + (nu/xi_0) * W_s +// - 1/2*[ csw_t * sum_s (sigma_ts F_ts) + (csw_s/xi_0) * sum_ss (sigma_ss F_ss) ] +// +// s spatial, t temporal directions. +// where W_t and W_s are the temporal and spatial components of the +// Wilson Dirac operator +// +// csw_r = csw_t to recover the isotropic version +////////////////////////////////////////////////////////////////// + +template +class WilsonCloverFermion : public WilsonFermion +{ +public: + // Types definitions + INHERIT_IMPL_TYPES(Impl); + template + using iImplClover = iScalar, Ns>>; + typedef iImplClover SiteCloverType; + typedef Lattice CloverFieldType; + +public: + typedef WilsonFermion WilsonBase; + + virtual void Instantiatable(void){}; + // Constructors + WilsonCloverFermion(GaugeField &_Umu, GridCartesian &Fgrid, + GridRedBlackCartesian &Hgrid, + const RealD _mass, + const RealD _csw_r = 0.0, + const RealD _csw_t = 0.0, + const WilsonAnisotropyCoefficients &clover_anisotropy = WilsonAnisotropyCoefficients(), + const ImplParams &impl_p = ImplParams()) : WilsonFermion(_Umu, + Fgrid, + Hgrid, + _mass, impl_p, clover_anisotropy), + CloverTerm(&Fgrid), + CloverTermInv(&Fgrid), + CloverTermEven(&Hgrid), + CloverTermOdd(&Hgrid), + CloverTermInvEven(&Hgrid), + CloverTermInvOdd(&Hgrid), + CloverTermDagEven(&Hgrid), + CloverTermDagOdd(&Hgrid), + CloverTermInvDagEven(&Hgrid), + CloverTermInvDagOdd(&Hgrid) + { + assert(Nd == 4); // require 4 dimensions + + if (clover_anisotropy.isAnisotropic) + { + csw_r = _csw_r * 0.5 / clover_anisotropy.xi_0; + diag_mass = _mass + 1.0 + (Nd - 1) * (clover_anisotropy.nu / clover_anisotropy.xi_0); + } + else + { + csw_r = _csw_r * 0.5; + diag_mass = 4.0 + _mass; + } + csw_t = _csw_t * 0.5; + + if (csw_r == 0) + std::cout << GridLogWarning << "Initializing WilsonCloverFermion with csw_r = 0" << std::endl; + if (csw_t == 0) + std::cout << GridLogWarning << "Initializing WilsonCloverFermion with csw_t = 0" << std::endl; + + ImportGauge(_Umu); + } + + virtual RealD M(const FermionField &in, FermionField &out); + virtual RealD Mdag(const FermionField &in, FermionField &out); + + virtual void Mooee(const FermionField &in, FermionField &out); + virtual void MooeeDag(const FermionField &in, FermionField &out); + virtual void MooeeInv(const FermionField &in, FermionField &out); + virtual void MooeeInvDag(const FermionField &in, FermionField &out); + virtual void MooeeInternal(const FermionField &in, FermionField &out, int dag, int inv); + + //virtual void MDeriv(GaugeField &mat, const FermionField &U, const FermionField &V, int dag); + virtual void MooDeriv(GaugeField &mat, const FermionField &U, const FermionField &V, int dag); + virtual void MeeDeriv(GaugeField &mat, const FermionField &U, const FermionField &V, int dag); + + void ImportGauge(const GaugeField &_Umu); + + // Derivative parts unpreconditioned pseudofermions + void MDeriv(GaugeField &force, const FermionField &X, const FermionField &Y, int dag) + { + conformable(X._grid, Y._grid); + conformable(X._grid, force._grid); + GaugeLinkField force_mu(force._grid), lambda(force._grid); + GaugeField clover_force(force._grid); + PropagatorField Lambda(force._grid); + + // Guido: Here we are hitting some performance issues: + // need to extract the components of the DoubledGaugeField + // for each call + // Possible solution + // Create a vector object to store them? (cons: wasting space) + std::vector U(Nd, this->Umu._grid); + + Impl::extractLinkField(U, this->Umu); + + force = zero; + // Derivative of the Wilson hopping term + this->DhopDeriv(force, X, Y, dag); + + /////////////////////////////////////////////////////////// + // Clover term derivative + /////////////////////////////////////////////////////////// + Impl::outerProductImpl(Lambda, X, Y); + //std::cout << "Lambda:" << Lambda << std::endl; + + Gamma::Algebra sigma[] = { + Gamma::Algebra::SigmaXY, + Gamma::Algebra::SigmaXZ, + Gamma::Algebra::SigmaXT, + Gamma::Algebra::MinusSigmaXY, + Gamma::Algebra::SigmaYZ, + Gamma::Algebra::SigmaYT, + Gamma::Algebra::MinusSigmaXZ, + Gamma::Algebra::MinusSigmaYZ, + Gamma::Algebra::SigmaZT, + Gamma::Algebra::MinusSigmaXT, + Gamma::Algebra::MinusSigmaYT, + Gamma::Algebra::MinusSigmaZT}; + + /* + sigma_{\mu \nu}= + | 0 sigma[0] sigma[1] sigma[2] | + | sigma[3] 0 sigma[4] sigma[5] | + | sigma[6] sigma[7] 0 sigma[8] | + | sigma[9] sigma[10] sigma[11] 0 | + */ + + int count = 0; + clover_force = zero; + for (int mu = 0; mu < 4; mu++) + { + force_mu = zero; + for (int nu = 0; nu < 4; nu++) + { + if (mu == nu) + continue; + + RealD factor; + if (nu == 4 || mu == 4) + { + factor = 2.0 * csw_t; + } + else + { + factor = 2.0 * csw_r; + } + PropagatorField Slambda = Gamma(sigma[count]) * Lambda; // sigma checked + Impl::TraceSpinImpl(lambda, Slambda); // traceSpin ok + force_mu -= factor*Cmunu(U, lambda, mu, nu); // checked + count++; + } + + pokeLorentz(clover_force, U[mu] * force_mu, mu); + } + //clover_force *= csw; + force += clover_force; + } + + // Computing C_{\mu \nu}(x) as in Eq.(B.39) in Zbigniew Sroczynski's PhD thesis + GaugeLinkField Cmunu(std::vector &U, GaugeLinkField &lambda, int mu, int nu) + { + conformable(lambda._grid, U[0]._grid); + GaugeLinkField out(lambda._grid), tmp(lambda._grid); + // insertion in upper staple + // please check redundancy of shift operations + + // C1+ + tmp = lambda * U[nu]; + out = Impl::ShiftStaple(Impl::CovShiftForward(tmp, nu, Impl::CovShiftBackward(U[mu], mu, Impl::CovShiftIdentityBackward(U[nu], nu))), mu); + + // C2+ + tmp = U[mu] * Impl::ShiftStaple(adj(lambda), mu); + out += Impl::ShiftStaple(Impl::CovShiftForward(U[nu], nu, Impl::CovShiftBackward(tmp, mu, Impl::CovShiftIdentityBackward(U[nu], nu))), mu); + + // C3+ + tmp = U[nu] * Impl::ShiftStaple(adj(lambda), nu); + out += Impl::ShiftStaple(Impl::CovShiftForward(U[nu], nu, Impl::CovShiftBackward(U[mu], mu, Impl::CovShiftIdentityBackward(tmp, nu))), mu); + + // C4+ + out += Impl::ShiftStaple(Impl::CovShiftForward(U[nu], nu, Impl::CovShiftBackward(U[mu], mu, Impl::CovShiftIdentityBackward(U[nu], nu))), mu) * lambda; + + // insertion in lower staple + // C1- + out -= Impl::ShiftStaple(lambda, mu) * Impl::ShiftStaple(Impl::CovShiftBackward(U[nu], nu, Impl::CovShiftBackward(U[mu], mu, U[nu])), mu); + + // C2- + tmp = adj(lambda) * U[nu]; + out -= Impl::ShiftStaple(Impl::CovShiftBackward(tmp, nu, Impl::CovShiftBackward(U[mu], mu, U[nu])), mu); + + // C3- + tmp = lambda * U[nu]; + out -= Impl::ShiftStaple(Impl::CovShiftBackward(U[nu], nu, Impl::CovShiftBackward(U[mu], mu, tmp)), mu); + + // C4- + out -= Impl::ShiftStaple(Impl::CovShiftBackward(U[nu], nu, Impl::CovShiftBackward(U[mu], mu, U[nu])), mu) * lambda; + + return out; + } + +private: + // here fixing the 4 dimensions, make it more general? + + RealD csw_r; // Clover coefficient - spatial + RealD csw_t; // Clover coefficient - temporal + RealD diag_mass; // Mass term + CloverFieldType CloverTerm, CloverTermInv; // Clover term + CloverFieldType CloverTermEven, CloverTermOdd; // Clover term EO + CloverFieldType CloverTermInvEven, CloverTermInvOdd; // Clover term Inv EO + CloverFieldType CloverTermDagEven, CloverTermDagOdd; // Clover term Dag EO + CloverFieldType CloverTermInvDagEven, CloverTermInvDagOdd; // Clover term Inv Dag EO + + // eventually these can be compressed into 6x6 blocks instead of the 12x12 + // using the DeGrand-Rossi basis for the gamma matrices + CloverFieldType fillCloverYZ(const GaugeLinkField &F) + { + CloverFieldType T(F._grid); + T = zero; + PARALLEL_FOR_LOOP + for (int i = 0; i < CloverTerm._grid->oSites(); i++) + { + T._odata[i]()(0, 1) = timesMinusI(F._odata[i]()()); + T._odata[i]()(1, 0) = timesMinusI(F._odata[i]()()); + T._odata[i]()(2, 3) = timesMinusI(F._odata[i]()()); + T._odata[i]()(3, 2) = timesMinusI(F._odata[i]()()); + } + + return T; + } + + CloverFieldType fillCloverXZ(const GaugeLinkField &F) + { + CloverFieldType T(F._grid); + T = zero; + PARALLEL_FOR_LOOP + for (int i = 0; i < CloverTerm._grid->oSites(); i++) + { + T._odata[i]()(0, 1) = -F._odata[i]()(); + T._odata[i]()(1, 0) = F._odata[i]()(); + T._odata[i]()(2, 3) = -F._odata[i]()(); + T._odata[i]()(3, 2) = F._odata[i]()(); + } + + return T; + } + + CloverFieldType fillCloverXY(const GaugeLinkField &F) + { + CloverFieldType T(F._grid); + T = zero; + PARALLEL_FOR_LOOP + for (int i = 0; i < CloverTerm._grid->oSites(); i++) + { + + T._odata[i]()(0, 0) = timesMinusI(F._odata[i]()()); + T._odata[i]()(1, 1) = timesI(F._odata[i]()()); + T._odata[i]()(2, 2) = timesMinusI(F._odata[i]()()); + T._odata[i]()(3, 3) = timesI(F._odata[i]()()); + } + + return T; + } + + CloverFieldType fillCloverXT(const GaugeLinkField &F) + { + CloverFieldType T(F._grid); + T = zero; + PARALLEL_FOR_LOOP + for (int i = 0; i < CloverTerm._grid->oSites(); i++) + { + T._odata[i]()(0, 1) = timesI(F._odata[i]()()); + T._odata[i]()(1, 0) = timesI(F._odata[i]()()); + T._odata[i]()(2, 3) = timesMinusI(F._odata[i]()()); + T._odata[i]()(3, 2) = timesMinusI(F._odata[i]()()); + } + + return T; + } + + CloverFieldType fillCloverYT(const GaugeLinkField &F) + { + CloverFieldType T(F._grid); + T = zero; + PARALLEL_FOR_LOOP + for (int i = 0; i < CloverTerm._grid->oSites(); i++) + { + T._odata[i]()(0, 1) = -(F._odata[i]()()); + T._odata[i]()(1, 0) = (F._odata[i]()()); + T._odata[i]()(2, 3) = (F._odata[i]()()); + T._odata[i]()(3, 2) = -(F._odata[i]()()); + } + + return T; + } + + CloverFieldType fillCloverZT(const GaugeLinkField &F) + { + CloverFieldType T(F._grid); + T = zero; + PARALLEL_FOR_LOOP + for (int i = 0; i < CloverTerm._grid->oSites(); i++) + { + T._odata[i]()(0, 0) = timesI(F._odata[i]()()); + T._odata[i]()(1, 1) = timesMinusI(F._odata[i]()()); + T._odata[i]()(2, 2) = timesMinusI(F._odata[i]()()); + T._odata[i]()(3, 3) = timesI(F._odata[i]()()); + } + + return T; + } +}; +} +} + +#endif // GRID_QCD_WILSON_CLOVER_FERMION_H diff --git a/lib/qcd/action/fermion/WilsonFermion.cc b/lib/qcd/action/fermion/WilsonFermion.cc index 1a020e8a..2336f4b6 100644 --- a/lib/qcd/action/fermion/WilsonFermion.cc +++ b/lib/qcd/action/fermion/WilsonFermion.cc @@ -47,7 +47,8 @@ int WilsonFermionStatic::HandOptDslash; template WilsonFermion::WilsonFermion(GaugeField &_Umu, GridCartesian &Fgrid, GridRedBlackCartesian &Hgrid, RealD _mass, - const ImplParams &p) + const ImplParams &p, + const WilsonAnisotropyCoefficients &anis) : Kernels(p), _grid(&Fgrid), _cbgrid(&Hgrid), @@ -60,16 +61,41 @@ WilsonFermion::WilsonFermion(GaugeField &_Umu, GridCartesian &Fgrid, Umu(&Fgrid), UmuEven(&Hgrid), UmuOdd(&Hgrid), - _tmp(&Hgrid) + _tmp(&Hgrid), + anisotropyCoeff(anis) { // Allocate the required comms buffer ImportGauge(_Umu); + if (anisotropyCoeff.isAnisotropic){ + diag_mass = mass + 1.0 + (Nd-1)*(anisotropyCoeff.nu / anisotropyCoeff.xi_0); + } else { + diag_mass = 4.0 + mass; + } + + } template void WilsonFermion::ImportGauge(const GaugeField &_Umu) { GaugeField HUmu(_Umu._grid); - HUmu = _Umu * (-0.5); + + //Here multiply the anisotropy coefficients + if (anisotropyCoeff.isAnisotropic) + { + + for (int mu = 0; mu < Nd; mu++) + { + GaugeLinkField U_dir = (-0.5)*PeekIndex(_Umu, mu); + if (mu != anisotropyCoeff.t_direction) + U_dir *= (anisotropyCoeff.nu / anisotropyCoeff.xi_0); + + PokeIndex(HUmu, U_dir, mu); + } + } + else + { + HUmu = _Umu * (-0.5); + } Impl::DoubleStore(GaugeGrid(), Umu, HUmu); pickCheckerboard(Even, UmuEven, Umu); pickCheckerboard(Odd, UmuOdd, Umu); @@ -83,14 +109,14 @@ template RealD WilsonFermion::M(const FermionField &in, FermionField &out) { out.checkerboard = in.checkerboard; Dhop(in, out, DaggerNo); - return axpy_norm(out, 4 + mass, in, out); + return axpy_norm(out, diag_mass, in, out); } template RealD WilsonFermion::Mdag(const FermionField &in, FermionField &out) { out.checkerboard = in.checkerboard; Dhop(in, out, DaggerYes); - return axpy_norm(out, 4 + mass, in, out); + return axpy_norm(out, diag_mass, in, out); } template @@ -114,7 +140,7 @@ void WilsonFermion::MeooeDag(const FermionField &in, FermionField &out) { template void WilsonFermion::Mooee(const FermionField &in, FermionField &out) { out.checkerboard = in.checkerboard; - typename FermionField::scalar_type scal(4.0 + mass); + typename FermionField::scalar_type scal(diag_mass); out = scal * in; } @@ -127,7 +153,7 @@ void WilsonFermion::MooeeDag(const FermionField &in, FermionField &out) { template void WilsonFermion::MooeeInv(const FermionField &in, FermionField &out) { out.checkerboard = in.checkerboard; - out = (1.0/(4.0+mass))*in; + out = (1.0/(diag_mass))*in; } template @@ -204,7 +230,7 @@ void WilsonFermion::DerivInternal(StencilImpl &st, DoubledGaugeField &U, FermionField Btilde(B._grid); FermionField Atilde(B._grid); - Atilde = A; + Atilde = A;//redundant st.HaloExchange(B, compressor); diff --git a/lib/qcd/action/fermion/WilsonFermion.h b/lib/qcd/action/fermion/WilsonFermion.h index feba40ed..c0db827e 100644 --- a/lib/qcd/action/fermion/WilsonFermion.h +++ b/lib/qcd/action/fermion/WilsonFermion.h @@ -44,6 +44,21 @@ class WilsonFermionStatic { static const int npoint = 8; }; + struct WilsonAnisotropyCoefficients: Serializable + { + GRID_SERIALIZABLE_CLASS_MEMBERS(WilsonAnisotropyCoefficients, + bool, isAnisotropic, + int, t_direction, + double, xi_0, + double, nu); + + WilsonAnisotropyCoefficients(): + isAnisotropic(false), + t_direction(Nd-1), + xi_0(1.0), + nu(1.0){} +}; + template class WilsonFermion : public WilsonKernels, public WilsonFermionStatic { public: @@ -65,8 +80,8 @@ class WilsonFermion : public WilsonKernels, public WilsonFermionStatic { // override multiply; cut number routines if pass dagger argument // and also make interface more uniformly consistent ////////////////////////////////////////////////////////////////// - RealD M(const FermionField &in, FermionField &out); - RealD Mdag(const FermionField &in, FermionField &out); + virtual RealD M(const FermionField &in, FermionField &out); + virtual RealD Mdag(const FermionField &in, FermionField &out); ///////////////////////////////////////////////////////// // half checkerboard operations @@ -117,8 +132,9 @@ class WilsonFermion : public WilsonKernels, public WilsonFermionStatic { // Constructor WilsonFermion(GaugeField &_Umu, GridCartesian &Fgrid, - GridRedBlackCartesian &Hgrid, RealD _mass, - const ImplParams &p = ImplParams()); + GridRedBlackCartesian &Hgrid, RealD _mass, + const ImplParams &p = ImplParams(), + const WilsonAnisotropyCoefficients &anis = WilsonAnisotropyCoefficients() ); // DoubleStore impl dependent void ImportGauge(const GaugeField &_Umu); @@ -130,6 +146,7 @@ class WilsonFermion : public WilsonKernels, public WilsonFermionStatic { // protected: public: RealD mass; + RealD diag_mass; GridBase *_grid; GridBase *_cbgrid; @@ -146,6 +163,8 @@ class WilsonFermion : public WilsonKernels, public WilsonFermionStatic { LebesgueOrder Lebesgue; LebesgueOrder LebesgueEvenOdd; + + WilsonAnisotropyCoefficients anisotropyCoeff; /////////////////////////////////////////////////////////////// // Conserved current utilities diff --git a/lib/qcd/action/fermion/WilsonKernelsHand.cc b/lib/qcd/action/fermion/WilsonKernelsHand.cc index 80b81714..aa6b5f6b 100644 --- a/lib/qcd/action/fermion/WilsonKernelsHand.cc +++ b/lib/qcd/action/fermion/WilsonKernelsHand.cc @@ -946,5 +946,6 @@ INSTANTIATE_THEM(DomainWallVec5dImplFH); INSTANTIATE_THEM(DomainWallVec5dImplDF); INSTANTIATE_THEM(ZDomainWallVec5dImplFH); INSTANTIATE_THEM(ZDomainWallVec5dImplDF); - +INSTANTIATE_THEM(WilsonTwoIndexAntiSymmetricImplF); +INSTANTIATE_THEM(WilsonTwoIndexAntiSymmetricImplD); }} diff --git a/lib/qcd/action/scalar/ScalarImpl.h b/lib/qcd/action/scalar/ScalarImpl.h index f85ab840..55f5049d 100644 --- a/lib/qcd/action/scalar/ScalarImpl.h +++ b/lib/qcd/action/scalar/ScalarImpl.h @@ -16,12 +16,12 @@ class ScalarImplTypes { typedef iImplField SiteField; typedef SiteField SitePropagator; typedef SiteField SiteComplex; - + typedef Lattice Field; typedef Field ComplexField; typedef Field FermionField; typedef Field PropagatorField; - + static inline void generate_momenta(Field& P, GridParallelRNG& pRNG){ gaussian(pRNG, P); } @@ -47,54 +47,60 @@ class ScalarImplTypes { static inline void ColdConfiguration(GridParallelRNG &pRNG, Field &U) { U = 1.0; } - + static void MomentumSpacePropagator(Field &out, RealD m) { GridBase *grid = out._grid; Field kmu(grid), one(grid); const unsigned int nd = grid->_ndimension; std::vector &l = grid->_fdimensions; - + one = Complex(1.0,0.0); out = m*m; for(int mu = 0; mu < nd; mu++) { Real twoPiL = M_PI*2./l[mu]; - + LatticeCoordinate(kmu,mu); kmu = 2.*sin(.5*twoPiL*kmu); out = out + kmu*kmu; } out = one/out; } - + static void FreePropagator(const Field &in, Field &out, const Field &momKernel) { FFT fft((GridCartesian *)in._grid); Field inFT(in._grid); - + fft.FFT_all_dim(inFT, in, FFT::forward); inFT = inFT*momKernel; fft.FFT_all_dim(out, inFT, FFT::backward); } - + static void FreePropagator(const Field &in, Field &out, RealD m) { Field momKernel(in._grid); - + MomentumSpacePropagator(momKernel, m); FreePropagator(in, out, momKernel); } - + }; + #ifdef USE_FFT_ACCELERATION + #ifndef FFT_MASS + #error "USE_FFT_ACCELERATION is defined but not FFT_MASS" + #endif + #endif + template class ScalarAdjMatrixImplTypes { public: typedef S Simd; typedef QCD::SU Group; - + template using iImplField = iScalar>>; template @@ -103,24 +109,119 @@ class ScalarImplTypes { typedef iImplField SiteField; typedef SiteField SitePropagator; typedef iImplComplex SiteComplex; - + typedef Lattice Field; typedef Lattice ComplexField; typedef Field FermionField; typedef Field PropagatorField; - static inline void generate_momenta(Field& P, GridParallelRNG& pRNG) { + static void MomentaSquare(ComplexField &out) + { + GridBase *grid = out._grid; + const std::vector &l = grid->FullDimensions(); + ComplexField kmu(grid); + + for (int mu = 0; mu < grid->Nd(); mu++) + { + Real twoPiL = M_PI * 2.0 / l[mu]; + LatticeCoordinate(kmu, mu); + kmu = 2.0 * sin(0.5 * twoPiL * kmu); + out += kmu * kmu; + } + } + + static void MomentumSpacePropagator(ComplexField &out, RealD m) + { + GridBase *grid = out._grid; + ComplexField one(grid); + one = Complex(1.0, 0.0); + out = m * m; + MomentaSquare(out); + out = one / out; + } + + static inline void generate_momenta(Field &P, GridParallelRNG &pRNG) + { +#ifndef USE_FFT_ACCELERATION Group::GaussianFundamentalLieAlgebraMatrix(pRNG, P); +#else + + Field Pgaussian(P._grid), Pp(P._grid); + ComplexField p2(P._grid); p2 = zero; + RealD M = FFT_MASS; + + Group::GaussianFundamentalLieAlgebraMatrix(pRNG, Pgaussian); + + FFT theFFT((GridCartesian*)P._grid); + theFFT.FFT_all_dim(Pp, Pgaussian, FFT::forward); + MomentaSquare(p2); + p2 += M * M; + p2 = sqrt(p2); + Pp *= p2; + theFFT.FFT_all_dim(P, Pp, FFT::backward); + +#endif //USE_FFT_ACCELERATION } static inline Field projectForce(Field& P) {return P;} - static inline void update_field(Field& P, Field& U, double ep) { - U += P*ep; + static inline void update_field(Field &P, Field &U, double ep) + { +#ifndef USE_FFT_ACCELERATION + double t0=usecond(); + U += P * ep; + double t1=usecond(); + double total_time = (t1-t0)/1e6; + std::cout << GridLogIntegrator << "Total time for updating field (s) : " << total_time << std::endl; +#else + // FFT transform P(x) -> P(p) + // divide by (M^2+p^2) M external parameter (how to pass?) + // P'(p) = P(p)/(M^2+p^2) + // Transform back -> P'(x) + // U += P'(x)*ep + + Field Pp(U._grid), P_FFT(U._grid); + static ComplexField p2(U._grid); + RealD M = FFT_MASS; + + FFT theFFT((GridCartesian*)U._grid); + theFFT.FFT_all_dim(Pp, P, FFT::forward); + + static bool first_call = true; + if (first_call) + { + // avoid recomputing + MomentumSpacePropagator(p2, M); + first_call = false; + } + Pp *= p2; + theFFT.FFT_all_dim(P_FFT, Pp, FFT::backward); + U += P_FFT * ep; + +#endif //USE_FFT_ACCELERATION } - static inline RealD FieldSquareNorm(Field& U) { - return (TensorRemove(sum(trace(U*U))).real()); + static inline RealD FieldSquareNorm(Field &U) + { +#ifndef USE_FFT_ACCELERATION + return (TensorRemove(sum(trace(U * U))).real()); +#else + // In case of Fourier acceleration we have to: + // compute U(p)*U(p)/(M^2+p^2)) Parseval theorem + // 1 FFT needed U(x) -> U(p) + // M to be passed + + FFT theFFT((GridCartesian*)U._grid); + Field Up(U._grid); + + theFFT.FFT_all_dim(Up, U, FFT::forward); + RealD M = FFT_MASS; + ComplexField p2(U._grid); + MomentumSpacePropagator(p2, M); + Field Up2 = Up * p2; + // from the definition of the DFT we need to divide by the volume + return (-TensorRemove(sum(trace(adj(Up) * Up2))).real() / U._grid->gSites()); +#endif //USE_FFT_ACCELERATION } static inline void HotConfiguration(GridParallelRNG &pRNG, Field &U) { @@ -146,7 +247,7 @@ class ScalarImplTypes { typedef ScalarImplTypes ScalarImplCR; typedef ScalarImplTypes ScalarImplCF; typedef ScalarImplTypes ScalarImplCD; - + // Hardcoding here the size of the matrices typedef ScalarAdjMatrixImplTypes ScalarAdjImplR; typedef ScalarAdjMatrixImplTypes ScalarAdjImplF; @@ -155,7 +256,7 @@ class ScalarImplTypes { template using ScalarNxNAdjImplR = ScalarAdjMatrixImplTypes; template using ScalarNxNAdjImplF = ScalarAdjMatrixImplTypes; template using ScalarNxNAdjImplD = ScalarAdjMatrixImplTypes; - + //} } diff --git a/lib/qcd/action/scalar/ScalarInteractionAction.h b/lib/qcd/action/scalar/ScalarInteractionAction.h index 4d189352..8738b647 100644 --- a/lib/qcd/action/scalar/ScalarInteractionAction.h +++ b/lib/qcd/action/scalar/ScalarInteractionAction.h @@ -30,119 +30,179 @@ directory #ifndef SCALAR_INT_ACTION_H #define SCALAR_INT_ACTION_H - // Note: this action can completely absorb the ScalarAction for real float fields // use the scalarObjs to generalise the structure -namespace Grid { - // FIXME drop the QCD namespace everywhere here +namespace Grid +{ +// FIXME drop the QCD namespace everywhere here - template - class ScalarInteractionAction : public QCD::Action { - public: - INHERIT_FIELD_TYPES(Impl); - private: - RealD mass_square; - RealD lambda; +template +class ScalarInteractionAction : public QCD::Action +{ +public: + INHERIT_FIELD_TYPES(Impl); +private: + RealD mass_square; + RealD lambda; + RealD g; + const unsigned int N = Impl::Group::Dimension; - typedef typename Field::vector_object vobj; - typedef CartesianStencil Stencil; + typedef typename Field::vector_object vobj; + typedef CartesianStencil Stencil; - SimpleCompressor compressor; - int npoint = 2*Ndim; - std::vector directions;// = {0,1,2,3,0,1,2,3}; // forcing 4 dimensions - std::vector displacements;// = {1,1,1,1, -1,-1,-1,-1}; + SimpleCompressor compressor; + int npoint = 2 * Ndim; + std::vector directions; // + std::vector displacements; // - - public: - - ScalarInteractionAction(RealD ms, RealD l) : mass_square(ms), lambda(l), displacements(2*Ndim,0), directions(2*Ndim,0){ - for (int mu = 0 ; mu < Ndim; mu++){ - directions[mu] = mu; directions[mu+Ndim] = mu; - displacements[mu] = 1; displacements[mu+Ndim] = -1; - } +public: + ScalarInteractionAction(RealD ms, RealD l, RealD gval) : mass_square(ms), lambda(l), g(gval), displacements(2 * Ndim, 0), directions(2 * Ndim, 0) + { + for (int mu = 0; mu < Ndim; mu++) + { + directions[mu] = mu; + directions[mu + Ndim] = mu; + displacements[mu] = 1; + displacements[mu + Ndim] = -1; } + } - virtual std::string LogParameters() { - std::stringstream sstream; - sstream << GridLogMessage << "[ScalarAction] lambda : " << lambda << std::endl; - sstream << GridLogMessage << "[ScalarAction] mass_square : " << mass_square << std::endl; - return sstream.str(); - } + virtual std::string LogParameters() + { + std::stringstream sstream; + sstream << GridLogMessage << "[ScalarAction] lambda : " << lambda << std::endl; + sstream << GridLogMessage << "[ScalarAction] mass_square : " << mass_square << std::endl; + sstream << GridLogMessage << "[ScalarAction] g : " << g << std::endl; + return sstream.str(); + } - virtual std::string action_name() {return "ScalarAction";} + virtual std::string action_name() { return "ScalarAction"; } - virtual void refresh(const Field &U, GridParallelRNG &pRNG) {} + virtual void refresh(const Field &U, GridParallelRNG &pRNG) {} - virtual RealD S(const Field &p) { - assert(p._grid->Nd() == Ndim); - static Stencil phiStencil(p._grid, npoint, 0, directions, displacements); - phiStencil.HaloExchange(p, compressor); - Field action(p._grid), pshift(p._grid), phisquared(p._grid); - phisquared = p*p; - action = (2.0*Ndim + mass_square)*phisquared - lambda/24.*phisquared*phisquared; - for (int mu = 0; mu < Ndim; mu++) { - // pshift = Cshift(p, mu, +1); // not efficient, implement with stencils - parallel_for (int i = 0; i < p._grid->oSites(); i++) { - int permute_type; - StencilEntry *SE; - vobj temp2; - const vobj *temp, *t_p; - - SE = phiStencil.GetEntry(permute_type, mu, i); - t_p = &p._odata[i]; - if ( SE->_is_local ) { - temp = &p._odata[SE->_offset]; - if ( SE->_permute ) { - permute(temp2, *temp, permute_type); - action._odata[i] -= temp2*(*t_p) + (*t_p)*temp2; - } else { - action._odata[i] -= (*temp)*(*t_p) + (*t_p)*(*temp); - } - } else { - action._odata[i] -= phiStencil.CommBuf()[SE->_offset]*(*t_p) + (*t_p)*phiStencil.CommBuf()[SE->_offset]; - } - } - // action -= pshift*p + p*pshift; - } - // NB the trace in the algebra is normalised to 1/2 - // minus sign coming from the antihermitian fields - return -(TensorRemove(sum(trace(action)))).real(); - }; - - virtual void deriv(const Field &p, Field &force) { - assert(p._grid->Nd() == Ndim); - force = (2.0*Ndim + mass_square)*p - lambda/12.*p*p*p; - // move this outside - static Stencil phiStencil(p._grid, npoint, 0, directions, displacements); - phiStencil.HaloExchange(p, compressor); - - //for (int mu = 0; mu < QCD::Nd; mu++) force -= Cshift(p, mu, -1) + Cshift(p, mu, 1); - for (int point = 0; point < npoint; point++) { - parallel_for (int i = 0; i < p._grid->oSites(); i++) { - const vobj *temp; - vobj temp2; - int permute_type; - StencilEntry *SE; - SE = phiStencil.GetEntry(permute_type, point, i); - - if ( SE->_is_local ) { - temp = &p._odata[SE->_offset]; - if ( SE->_permute ) { - permute(temp2, *temp, permute_type); - force._odata[i] -= temp2; - } else { - force._odata[i] -= *temp; - } - } else { - force._odata[i] -= phiStencil.CommBuf()[SE->_offset]; - } - } + virtual RealD S(const Field &p) + { + assert(p._grid->Nd() == Ndim); + static Stencil phiStencil(p._grid, npoint, 0, directions, displacements); + phiStencil.HaloExchange(p, compressor); + Field action(p._grid), pshift(p._grid), phisquared(p._grid); + phisquared = p * p; + action = (2.0 * Ndim + mass_square) * phisquared - lambda * phisquared * phisquared; + for (int mu = 0; mu < Ndim; mu++) + { + // pshift = Cshift(p, mu, +1); // not efficient, implement with stencils + parallel_for(int i = 0; i < p._grid->oSites(); i++) + { + int permute_type; + StencilEntry *SE; + vobj temp2; + const vobj *temp, *t_p; + + SE = phiStencil.GetEntry(permute_type, mu, i); + t_p = &p._odata[i]; + if (SE->_is_local) + { + temp = &p._odata[SE->_offset]; + if (SE->_permute) + { + permute(temp2, *temp, permute_type); + action._odata[i] -= temp2 * (*t_p) + (*t_p) * temp2; + } + else + { + action._odata[i] -= (*temp) * (*t_p) + (*t_p) * (*temp); + } + } + else + { + action._odata[i] -= phiStencil.CommBuf()[SE->_offset] * (*t_p) + (*t_p) * phiStencil.CommBuf()[SE->_offset]; + } } + // action -= pshift*p + p*pshift; } + // NB the trace in the algebra is normalised to 1/2 + // minus sign coming from the antihermitian fields + return -(TensorRemove(sum(trace(action)))).real() * N / g; }; - -} // namespace Grid -#endif // SCALAR_INT_ACTION_H + virtual void deriv(const Field &p, Field &force) + { + double t0 = usecond(); + assert(p._grid->Nd() == Ndim); + force = (2. * Ndim + mass_square) * p - 2. * lambda * p * p * p; + double interm_t = usecond(); + + // move this outside + static Stencil phiStencil(p._grid, npoint, 0, directions, displacements); + + phiStencil.HaloExchange(p, compressor); + double halo_t = usecond(); + int chunk = 128; + //for (int mu = 0; mu < QCD::Nd; mu++) force -= Cshift(p, mu, -1) + Cshift(p, mu, 1); + + // inverting the order of the loops slows down the code(! g++ 7) + // cannot try to reduce the number of force writes by factor npoint... + // use cache blocking + for (int point = 0; point < npoint; point++) + { + +#pragma omp parallel +{ + int permute_type; + StencilEntry *SE; + const vobj *temp; + +#pragma omp for schedule(static, chunk) + for (int i = 0; i < p._grid->oSites(); i++) + { + SE = phiStencil.GetEntry(permute_type, point, i); + // prefetch next p? + + if (SE->_is_local) + { + temp = &p._odata[SE->_offset]; + + if (SE->_permute) + { + vobj temp2; + permute(temp2, *temp, permute_type); + force._odata[i] -= temp2; + } + else + { + force._odata[i] -= *temp; // slow part. Dominated by this read/write (BW) + } + } + else + { + force._odata[i] -= phiStencil.CommBuf()[SE->_offset]; + } + } + + } + } + force *= N / g; + + double t1 = usecond(); + double total_time = (t1 - t0) / 1e6; + double interm_time = (interm_t - t0) / 1e6; + double halo_time = (halo_t - interm_t) / 1e6; + double stencil_time = (t1 - halo_t) / 1e6; + std::cout << GridLogIntegrator << "Total time for force computation (s) : " << total_time << std::endl; + std::cout << GridLogIntegrator << "Intermediate time for force computation (s): " << interm_time << std::endl; + std::cout << GridLogIntegrator << "Halo time in force computation (s) : " << halo_time << std::endl; + std::cout << GridLogIntegrator << "Stencil time in force computation (s) : " << stencil_time << std::endl; + double flops = p._grid->gSites() * (14 * N * N * N + 18 * N * N + 2); + double flops_no_stencil = p._grid->gSites() * (14 * N * N * N + 6 * N * N + 2); + double Gflops = flops / (total_time * 1e9); + double Gflops_no_stencil = flops_no_stencil / (interm_time * 1e9); + std::cout << GridLogIntegrator << "Flops: " << flops << " - Gflop/s : " << Gflops << std::endl; + std::cout << GridLogIntegrator << "Flops NS: " << flops_no_stencil << " - Gflop/s NS: " << Gflops_no_stencil << std::endl; +} +}; + +} // namespace Grid + +#endif // SCALAR_INT_ACTION_H diff --git a/lib/qcd/hmc/GenericHMCrunner.h b/lib/qcd/hmc/GenericHMCrunner.h index 4f6c1af0..26fec3d5 100644 --- a/lib/qcd/hmc/GenericHMCrunner.h +++ b/lib/qcd/hmc/GenericHMCrunner.h @@ -211,7 +211,7 @@ typedef HMCWrapperTemplate ScalarAdjGenericHMCRunner; template -using ScalarNxNAdjGenericHMCRunner = HMCWrapperTemplate < ScalarNxNAdjImplR, MinimumNorm2, ScalarNxNMatrixFields >; +using ScalarNxNAdjGenericHMCRunner = HMCWrapperTemplate < ScalarNxNAdjImplR, ForceGradient, ScalarNxNMatrixFields >; } // namespace QCD } // namespace Grid diff --git a/lib/qcd/modules/ObservableModules.h b/lib/qcd/modules/ObservableModules.h index 24511617..fbffc236 100644 --- a/lib/qcd/modules/ObservableModules.h +++ b/lib/qcd/modules/ObservableModules.h @@ -92,6 +92,19 @@ class PlaquetteMod: public ObservableModule, NoParameters> PlaquetteMod(): ObsBase(NoParameters()){} }; +template < class Impl > +class PolyakovMod: public ObservableModule, NoParameters>{ + typedef ObservableModule, NoParameters> ObsBase; + using ObsBase::ObsBase; // for constructors + + // acquire resource + virtual void initialize(){ + this->ObservablePtr.reset(new PolyakovLogger()); + } + public: + PolyakovMod(): ObsBase(NoParameters()){} +}; + template < class Impl > class TopologicalChargeMod: public ObservableModule, TopologyObsParameters>{ diff --git a/lib/qcd/observables/hmc_observable.h b/lib/qcd/observables/hmc_observable.h index db629ce7..fcf11774 100644 --- a/lib/qcd/observables/hmc_observable.h +++ b/lib/qcd/observables/hmc_observable.h @@ -45,5 +45,7 @@ class HmcObservable { #include "plaquette.h" #include "topological_charge.h" +#include "polyakov_loop.h" + #endif // HMC_OBSERVABLE_H diff --git a/lib/qcd/observables/polyakov_loop.h b/lib/qcd/observables/polyakov_loop.h new file mode 100644 index 00000000..d708b474 --- /dev/null +++ b/lib/qcd/observables/polyakov_loop.h @@ -0,0 +1,68 @@ +/************************************************************************************* + +Grid physics library, www.github.com/paboyle/Grid + +Source file: ./lib/qcd/modules/polyakov_line.h + +Copyright (C) 2017 + +Author: David Preti + +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 HMC_POLYAKOV_H +#define HMC_POLYAKOV_H + +namespace Grid { +namespace QCD { + +// this is only defined for a gauge theory +template +class PolyakovLogger : public HmcObservable { + public: + // here forces the Impl to be of gauge fields + // if not the compiler will complain + INHERIT_GIMPL_TYPES(Impl); + + // necessary for HmcObservable compatibility + typedef typename Impl::Field Field; + + void TrajectoryComplete(int traj, + Field &U, + GridSerialRNG &sRNG, + GridParallelRNG &pRNG) { + + ComplexD polyakov = WilsonLoops::avgPolyakovLoop(U); + + int def_prec = std::cout.precision(); + + std::cout << GridLogMessage + << std::setprecision(std::numeric_limits::digits10 + 1) + << "Polyakov Loop: [ " << traj << " ] "<< polyakov << std::endl; + + std::cout.precision(def_prec); + + } +}; + +} // namespace QCD +} // namespace Grid + +#endif // HMC_POLYAKOV_H diff --git a/lib/qcd/utils/WilsonLoops.h b/lib/qcd/utils/WilsonLoops.h index 90d905d3..cdd76ecc 100644 --- a/lib/qcd/utils/WilsonLoops.h +++ b/lib/qcd/utils/WilsonLoops.h @@ -123,6 +123,28 @@ public: return sumplaq / vol / faces / Nc; // Nd , Nc dependent... FIXME } + + ////////////////////////////////////////////////// + // average over all x,y,z the temporal loop + ////////////////////////////////////////////////// + static ComplexD avgPolyakovLoop(const GaugeField &Umu) { //assume Nd=4 + GaugeMat Ut(Umu._grid), P(Umu._grid); + ComplexD out; + int T = Umu._grid->GlobalDimensions()[3]; + int X = Umu._grid->GlobalDimensions()[0]; + int Y = Umu._grid->GlobalDimensions()[1]; + int Z = Umu._grid->GlobalDimensions()[2]; + + Ut = peekLorentz(Umu,3); //Select temporal direction + P = Ut; + for (int t=1;t--+ +-->--+(x) + // (x)+-->--+ +-->--+(x) - h.c. // | | | | // +--<--+ +--<--+ @@ -335,7 +359,9 @@ static void StapleMult(GaugeMat &staple, const GaugeLorentz &Umu, int mu) { GaugeMat v = Vup - Vdn; GaugeMat u = PeekIndex(Umu, mu); // some redundant copies GaugeMat vu = v*u; - FS = 0.25*Ta(u*v + Cshift(vu, mu, -1)); + //FS = 0.25*Ta(u*v + Cshift(vu, mu, -1)); + FS = (u*v + Cshift(vu, mu, -1)); + FS = 0.125*(FS - adj(FS)); } static Real TopologicalCharge(GaugeLorentz &U){ @@ -360,6 +386,7 @@ static void StapleMult(GaugeMat &staple, const GaugeLorentz &Umu, int mu) { return TensorRemove(Tq).real(); } + ////////////////////////////////////////////////////// // Similar to above for rectangle is required ////////////////////////////////////////////////////// diff --git a/scripts/copyright b/scripts/copyright index cc9ed6e5..a461b54c 100755 --- a/scripts/copyright +++ b/scripts/copyright @@ -11,8 +11,7 @@ Grid physics library, www.github.com/paboyle/Grid Source file: $1 -Copyright (C) 2015 -Copyright (C) 2016 +Copyright (C) 2015-2018 EOF @@ -60,4 +59,4 @@ shift done - +rm message tmp.fil diff --git a/tests/core/Test_wilson_clover.cc b/tests/core/Test_wilson_clover.cc new file mode 100644 index 00000000..9281e298 --- /dev/null +++ b/tests/core/Test_wilson_clover.cc @@ -0,0 +1,357 @@ +/************************************************************************************* + + Grid physics library, www.github.com/paboyle/Grid + + Source file: ./benchmarks/Benchmark_wilson.cc + + Copyright (C) 2015 + + Author: Guido Cossu + + 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 std; +using namespace Grid; +using namespace Grid::QCD; + +int main(int argc, char **argv) +{ + Grid_init(&argc, &argv); + + std::vector latt_size = GridDefaultLatt(); + std::vector simd_layout = GridDefaultSimd(Nd, vComplex::Nsimd()); + std::vector mpi_layout = GridDefaultMpi(); + GridCartesian Grid(latt_size, simd_layout, mpi_layout); + GridRedBlackCartesian RBGrid(&Grid); + + int threads = GridThread::GetThreads(); + std::cout << GridLogMessage << "Grid is setup to use " << threads << " threads" << std::endl; + std::cout << GridLogMessage << "Grid floating point word size is REALF" << sizeof(RealF) << std::endl; + std::cout << GridLogMessage << "Grid floating point word size is REALD" << sizeof(RealD) << std::endl; + std::cout << GridLogMessage << "Grid floating point word size is REAL" << sizeof(Real) << std::endl; + + std::vector seeds({1, 2, 3, 4}); + GridParallelRNG pRNG(&Grid); + pRNG.SeedFixedIntegers(seeds); + // pRNG.SeedFixedIntegers(std::vector({45,12,81,9}); + + typedef typename WilsonCloverFermionR::FermionField FermionField; + typename WilsonCloverFermionR::ImplParams params; + WilsonAnisotropyCoefficients anis; + + FermionField src(&Grid); + random(pRNG, src); + FermionField result(&Grid); + result = zero; + FermionField result2(&Grid); + result2 = zero; + FermionField ref(&Grid); + ref = zero; + FermionField tmp(&Grid); + tmp = zero; + FermionField err(&Grid); + err = zero; + FermionField err2(&Grid); + err2 = zero; + FermionField phi(&Grid); + random(pRNG, phi); + FermionField chi(&Grid); + random(pRNG, chi); + LatticeGaugeField Umu(&Grid); + SU3::HotConfiguration(pRNG, Umu); + std::vector U(4, &Grid); + + double volume = 1; + for (int mu = 0; mu < Nd; mu++) + { + volume = volume * latt_size[mu]; + } + + RealD mass = 0.1; + RealD csw_r = 1.0; + RealD csw_t = 1.0; + + WilsonCloverFermionR Dwc(Umu, Grid, RBGrid, mass, csw_r, csw_t, anis, params); + //Dwc.ImportGauge(Umu); // not necessary, included in the constructor + + std::cout << GridLogMessage << "==========================================================" << std::endl; + std::cout << GridLogMessage << "= Testing that Deo + Doe = Dunprec " << std::endl; + std::cout << GridLogMessage << "==========================================================" << std::endl; + + FermionField src_e(&RBGrid); + FermionField src_o(&RBGrid); + FermionField r_e(&RBGrid); + FermionField r_o(&RBGrid); + FermionField r_eo(&Grid); + pickCheckerboard(Even, src_e, src); + pickCheckerboard(Odd, src_o, src); + + Dwc.Meooe(src_e, r_o); + std::cout << GridLogMessage << "Applied Meo" << std::endl; + Dwc.Meooe(src_o, r_e); + std::cout << GridLogMessage << "Applied Moe" << std::endl; + Dwc.Dhop(src, ref, DaggerNo); + + setCheckerboard(r_eo, r_o); + setCheckerboard(r_eo, r_e); + + err = ref - r_eo; + std::cout << GridLogMessage << "EO norm diff " << norm2(err) << " " << norm2(ref) << " " << norm2(r_eo) << std::endl; + + std::cout << GridLogMessage << "==============================================================" << std::endl; + std::cout << GridLogMessage << "= Test Ddagger is the dagger of D by requiring " << std::endl; + std::cout << GridLogMessage << "= < phi | Deo | chi > * = < chi | Deo^dag| phi> " << std::endl; + std::cout << GridLogMessage << "==============================================================" << std::endl; + + FermionField chi_e(&RBGrid); + FermionField chi_o(&RBGrid); + + FermionField dchi_e(&RBGrid); + FermionField dchi_o(&RBGrid); + + FermionField phi_e(&RBGrid); + FermionField phi_o(&RBGrid); + + FermionField dphi_e(&RBGrid); + FermionField dphi_o(&RBGrid); + + pickCheckerboard(Even, chi_e, chi); + pickCheckerboard(Odd, chi_o, chi); + pickCheckerboard(Even, phi_e, phi); + pickCheckerboard(Odd, phi_o, phi); + + Dwc.Meooe(chi_e, dchi_o); + Dwc.Meooe(chi_o, dchi_e); + Dwc.MeooeDag(phi_e, dphi_o); + Dwc.MeooeDag(phi_o, dphi_e); + + ComplexD pDce = innerProduct(phi_e, dchi_e); + ComplexD pDco = innerProduct(phi_o, dchi_o); + ComplexD cDpe = innerProduct(chi_e, dphi_e); + ComplexD cDpo = innerProduct(chi_o, dphi_o); + + std::cout << GridLogMessage << "e " << pDce << " " << cDpe << std::endl; + std::cout << GridLogMessage << "o " << pDco << " " << cDpo << std::endl; + + std::cout << GridLogMessage << "pDce - conj(cDpo) " << pDce - conj(cDpo) << std::endl; + std::cout << GridLogMessage << "pDco - conj(cDpe) " << pDco - conj(cDpe) << std::endl; + + std::cout << GridLogMessage << "==============================================================" << std::endl; + std::cout << GridLogMessage << "= Test MeeInv Mee = 1 (if csw!=0) " << std::endl; + std::cout << GridLogMessage << "==============================================================" << std::endl; + + pickCheckerboard(Even, chi_e, chi); + pickCheckerboard(Odd, chi_o, chi); + + Dwc.Mooee(chi_e, src_e); + Dwc.MooeeInv(src_e, phi_e); + + Dwc.Mooee(chi_o, src_o); + Dwc.MooeeInv(src_o, phi_o); + + setCheckerboard(phi, phi_e); + setCheckerboard(phi, phi_o); + + err = phi - chi; + std::cout << GridLogMessage << "norm diff " << norm2(err) << std::endl; + + std::cout << GridLogMessage << "==============================================================" << std::endl; + std::cout << GridLogMessage << "= Test MeeDag MeeInvDag = 1 (if csw!=0) " << std::endl; + std::cout << GridLogMessage << "==============================================================" << std::endl; + + pickCheckerboard(Even, chi_e, chi); + pickCheckerboard(Odd, chi_o, chi); + + Dwc.MooeeDag(chi_e, src_e); + Dwc.MooeeInvDag(src_e, phi_e); + + Dwc.MooeeDag(chi_o, src_o); + Dwc.MooeeInvDag(src_o, phi_o); + + setCheckerboard(phi, phi_e); + setCheckerboard(phi, phi_o); + + err = phi - chi; + std::cout << GridLogMessage << "norm diff " << norm2(err) << std::endl; + + std::cout << GridLogMessage << "==============================================================" << std::endl; + std::cout << GridLogMessage << "= Test MeeInv MeeDag = 1 (if csw!=0) " << std::endl; + std::cout << GridLogMessage << "==============================================================" << std::endl; + + pickCheckerboard(Even, chi_e, chi); + pickCheckerboard(Odd, chi_o, chi); + + Dwc.MooeeDag(chi_e, src_e); + Dwc.MooeeInv(src_e, phi_e); + + Dwc.MooeeDag(chi_o, src_o); + Dwc.MooeeInv(src_o, phi_o); + + setCheckerboard(phi, phi_e); + setCheckerboard(phi, phi_o); + + err = phi - chi; + std::cout << GridLogMessage << "norm diff " << norm2(err) << std::endl; + + std::cout << GridLogMessage << "================================================================" << std::endl; + std::cout << GridLogMessage << "= Testing gauge covariance Clover term with EO preconditioning " << std::endl; + std::cout << GridLogMessage << "================================================================" << std::endl; + + chi = zero; + phi = zero; + tmp = zero; + pickCheckerboard(Even, chi_e, chi); + pickCheckerboard(Odd, chi_o, chi); + pickCheckerboard(Even, phi_e, phi); + pickCheckerboard(Odd, phi_o, phi); + + Dwc.Mooee(src_e, chi_e); + Dwc.Mooee(src_o, chi_o); + setCheckerboard(chi, chi_e); + setCheckerboard(chi, chi_o); + setCheckerboard(src, src_e); + setCheckerboard(src, src_o); + + ////////////////////// Gauge Transformation + std::vector seeds2({5, 6, 7, 8}); + GridParallelRNG pRNG2(&Grid); + pRNG2.SeedFixedIntegers(seeds2); + LatticeColourMatrix Omega(&Grid); + LatticeColourMatrix ShiftedOmega(&Grid); + LatticeGaugeField U_prime(&Grid); + U_prime = zero; + LatticeColourMatrix U_prime_mu(&Grid); + U_prime_mu = zero; + SU::LieRandomize(pRNG2, Omega, 1.0); + for (int mu = 0; mu < Nd; mu++) + { + U[mu] = peekLorentz(Umu, mu); + ShiftedOmega = Cshift(Omega, mu, 1); + U_prime_mu = Omega * U[mu] * adj(ShiftedOmega); + pokeLorentz(U_prime, U_prime_mu, mu); + } + ///////////////// + + WilsonCloverFermionR Dwc_prime(U_prime, Grid, RBGrid, mass, csw_r, csw_t, anis, params); + Dwc_prime.ImportGauge(U_prime); + + tmp = Omega * src; + pickCheckerboard(Even, src_e, tmp); + pickCheckerboard(Odd, src_o, tmp); + + Dwc_prime.Mooee(src_e, phi_e); + Dwc_prime.Mooee(src_o, phi_o); + + setCheckerboard(phi, phi_e); + setCheckerboard(phi, phi_o); + + err = chi - adj(Omega) * phi; + std::cout << GridLogMessage << "norm diff " << norm2(err) << std::endl; + + std::cout << GridLogMessage << "=================================================================" << std::endl; + std::cout << GridLogMessage << "= Testing gauge covariance Clover term w/o EO preconditioning " << std::endl; + std::cout << GridLogMessage << "================================================================" << std::endl; + + chi = zero; + phi = zero; + + WilsonFermionR Dw(Umu, Grid, RBGrid, mass, params); + Dw.ImportGauge(Umu); + + Dw.M(src, result); + Dwc.M(src, chi); + + Dwc_prime.M(Omega * src, phi); + + WilsonFermionR Dw_prime(U_prime, Grid, RBGrid, mass, params); + Dw_prime.ImportGauge(U_prime); + Dw_prime.M(Omega * src, result2); + + err = chi - adj(Omega) * phi; + err2 = result - adj(Omega) * result2; + std::cout << GridLogMessage << "norm diff Wilson " << norm2(err) << std::endl; + std::cout << GridLogMessage << "norm diff WilsonClover " << norm2(err2) << std::endl; + + std::cout << GridLogMessage << "==========================================================" << std::endl; + std::cout << GridLogMessage << "= Testing Mooee(csw=0) Clover to reproduce Mooee Wilson " << std::endl; + std::cout << GridLogMessage << "==========================================================" << std::endl; + + chi = zero; + phi = zero; + err = zero; + WilsonCloverFermionR Dwc_csw0(Umu, Grid, RBGrid, mass, 0.0, 0.0, anis, params); // <-- Notice: csw=0 + Dwc_csw0.ImportGauge(Umu); + + pickCheckerboard(Even, phi_e, phi); + pickCheckerboard(Odd, phi_o, phi); + pickCheckerboard(Even, chi_e, chi); + pickCheckerboard(Odd, chi_o, chi); + + Dw.Mooee(src_e, chi_e); + Dw.Mooee(src_o, chi_o); + Dwc_csw0.Mooee(src_e, phi_e); + Dwc_csw0.Mooee(src_o, phi_o); + + setCheckerboard(chi, chi_e); + setCheckerboard(chi, chi_o); + setCheckerboard(phi, phi_e); + setCheckerboard(phi, phi_o); + setCheckerboard(src, src_e); + setCheckerboard(src, src_o); + + err = chi - phi; + std::cout << GridLogMessage << "norm diff " << norm2(err) << std::endl; + + std::cout << GridLogMessage << "==========================================================" << std::endl; + std::cout << GridLogMessage << "= Testing EO operator is equal to the unprec " << std::endl; + std::cout << GridLogMessage << "==========================================================" << std::endl; + + chi = zero; + phi = zero; + err = zero; + + pickCheckerboard(Even, phi_e, phi); + pickCheckerboard(Odd, phi_o, phi); + pickCheckerboard(Even, chi_e, chi); + pickCheckerboard(Odd, chi_o, chi); + + // M phi = (Mooee src_e + Meooe src_o , Meooe src_e + Mooee src_o) + + Dwc.M(src, ref); // Reference result from the unpreconditioned operator + + // EO matrix + Dwc.Mooee(src_e, chi_e); + Dwc.Mooee(src_o, chi_o); + Dwc.Meooe(src_o, phi_e); + Dwc.Meooe(src_e, phi_o); + + phi_o += chi_o; + phi_e += chi_e; + + setCheckerboard(phi, phi_e); + setCheckerboard(phi, phi_o); + + err = ref - phi; + std::cout << GridLogMessage << "ref (unpreconditioned operator) diff :" << norm2(ref) << std::endl; + std::cout << GridLogMessage << "phi (EO decomposition) diff :" << norm2(phi) << std::endl; + std::cout << GridLogMessage << "norm diff :" << norm2(err) << std::endl; + + Grid_finalize(); +} diff --git a/tests/forces/Test_wilson_force.cc b/tests/forces/Test_wilson_force.cc index 1f34a48a..f834726b 100644 --- a/tests/forces/Test_wilson_force.cc +++ b/tests/forces/Test_wilson_force.cc @@ -50,7 +50,12 @@ int main (int argc, char ** argv) std::vector seeds({1,2,3,4}); GridParallelRNG pRNG(&Grid); - pRNG.SeedFixedIntegers(std::vector({45,12,81,9})); + std::vector vrand(4); + std::srand(std::time(0)); + std::generate(vrand.begin(), vrand.end(), std::rand); + std::cout << GridLogMessage << vrand << std::endl; + pRNG.SeedFixedIntegers(vrand); + //pRNG.SeedFixedIntegers(std::vector({45,12,81,9})); LatticeFermion phi (&Grid); gaussian(pRNG,phi); LatticeFermion Mphi (&Grid); diff --git a/tests/forces/Test_wilsonclover_force.cc b/tests/forces/Test_wilsonclover_force.cc new file mode 100644 index 00000000..f7090845 --- /dev/null +++ b/tests/forces/Test_wilsonclover_force.cc @@ -0,0 +1,194 @@ +/************************************************************************************* + + Grid physics library, www.github.com/paboyle/Grid + + Source file: ./tests/Test_wilson_force.cc + + Copyright (C) 2015 + +Author: Peter Boyle + + This program is free software; you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation; either version 2 of the License, or + (at your option) any later version. + + This program is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License along + with this program; if not, write to the Free Software Foundation, Inc., + 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. + + See the full license in the file "LICENSE" in the top level distribution directory + *************************************************************************************/ +/* END LEGAL */ +#include + +using namespace std; +using namespace Grid; +using namespace Grid::QCD; + +int main(int argc, char **argv) +{ + Grid_init(&argc, &argv); + + std::vector latt_size = GridDefaultLatt(); + std::vector simd_layout = GridDefaultSimd(Nd, vComplex::Nsimd()); + std::vector mpi_layout = GridDefaultMpi(); + + GridCartesian Grid(latt_size, simd_layout, mpi_layout); + GridRedBlackCartesian RBGrid(&Grid); + + int threads = GridThread::GetThreads(); + std::cout << GridLogMessage << "Grid is setup to use " << threads << " threads" << std::endl; + + std::vector seeds({1, 2, 30, 50}); + + GridParallelRNG pRNG(&Grid); + + std::vector vrand(4); + std::srand(std::time(0)); + std::generate(vrand.begin(), vrand.end(), std::rand); + std::cout << GridLogMessage << vrand << std::endl; + pRNG.SeedFixedIntegers(vrand); + //pRNG.SeedFixedIntegers(seeds); + + LatticeFermion phi(&Grid); + gaussian(pRNG, phi); + LatticeFermion Mphi(&Grid); + LatticeFermion MphiPrime(&Grid); + + LatticeGaugeField U(&Grid); + + std::vector site = {0, 0, 0, 0}; + SU3::HotConfiguration(pRNG, U); + //SU3::ColdConfiguration(pRNG, U);// Clover term zero + + //////////////////////////////////// + // Unmodified matrix element + //////////////////////////////////// + RealD mass = 0.1; + Real csw = 1.0; + WilsonCloverFermionR Dw(U, Grid, RBGrid, mass, csw, csw); + Dw.ImportGauge(U); + Dw.M(phi, Mphi); + ComplexD S = innerProduct(Mphi, Mphi); // Action : pdag MdagM p + + // get the deriv of phidag MdagM phi with respect to "U" + LatticeGaugeField UdSdU(&Grid); + LatticeGaugeField tmp(&Grid); + + //////////////////////////////////////////// + Dw.MDeriv(tmp, Mphi, phi, DaggerNo); + UdSdU = tmp; + Dw.MDeriv(tmp, phi, Mphi, DaggerYes); + UdSdU += tmp; + ///////////////////////////////////////////// + + //////////////////////////////////// + // Modify the gauge field a little + //////////////////////////////////// + RealD dt = 0.00005; + RealD Hmom = 0.0; + RealD Hmomprime = 0.0; + RealD Hmompp = 0.0; + LatticeColourMatrix mommu(&Grid); + LatticeColourMatrix forcemu(&Grid); + LatticeGaugeField mom(&Grid); + LatticeGaugeField Uprime(&Grid); + + for (int mu = 0; mu < Nd; mu++) + { + // Traceless antihermitian momentum; gaussian in lie alg + SU3::GaussianFundamentalLieAlgebraMatrix(pRNG, mommu); + Hmom -= real(sum(trace(mommu * mommu))); + PokeIndex(mom, mommu, mu); + + parallel_for(int ss = 0; ss < mom._grid->oSites(); ss++) + { + Uprime[ss]._internal[mu] = ProjectOnGroup(Exponentiate(mom[ss]._internal[mu], dt, 12) * U[ss]._internal[mu]); + } + } + + std::cout << GridLogMessage << "Initial mom hamiltonian is " << Hmom << std::endl; + + // New action + Dw.ImportGauge(Uprime); + Dw.M(phi, MphiPrime); + ComplexD Sprime = innerProduct(MphiPrime, MphiPrime); + + ////////////////////////////////////////////// + // Use derivative to estimate dS + ////////////////////////////////////////////// + + LatticeComplex dS(&Grid); + dS = zero; + LatticeComplex dSmom(&Grid); + dSmom = zero; + LatticeComplex dSmom2(&Grid); + dSmom2 = zero; + + for (int mu = 0; mu < Nd; mu++) + { + mommu = PeekIndex(UdSdU, mu); // P_mu = + mommu = Ta(mommu) * 2.0; // Mom = (P_mu - P_mu^dag) - trace(P_mu - P_mu^dag) + PokeIndex(UdSdU, mommu, mu); // UdSdU_mu = Mom + } + + std::cout << GridLogMessage << "Antihermiticity tests" << std::endl; + for (int mu = 0; mu < Nd; mu++) + { + mommu = PeekIndex(mom, mu); + std::cout << GridLogMessage << " Mommu " << norm2(mommu) << std::endl; + mommu = mommu + adj(mommu); + std::cout << GridLogMessage << " Mommu + Mommudag " << norm2(mommu) << std::endl; + mommu = PeekIndex(UdSdU, mu); + std::cout << GridLogMessage << " dsdumu " << norm2(mommu) << std::endl; + mommu = mommu + adj(mommu); + std::cout << GridLogMessage << " dsdumu + dag " << norm2(mommu) << std::endl; + std::cout << "" << std::endl; + } + ///////////////////////////////////////////////////// + + for (int mu = 0; mu < Nd; mu++) + { + forcemu = PeekIndex(UdSdU, mu); + mommu = PeekIndex(mom, mu); + + // Update PF action density + dS = dS + trace(mommu * forcemu) * dt; + + dSmom = dSmom - trace(mommu * forcemu) * dt; + dSmom2 = dSmom2 - trace(forcemu * forcemu) * (0.25 * dt * dt); + + // Update mom action density + mommu = mommu + forcemu * (dt * 0.5); + + Hmomprime -= real(sum(trace(mommu * mommu))); + } + + ComplexD dSpred = sum(dS); + ComplexD dSm = sum(dSmom); + ComplexD dSm2 = sum(dSmom2); + + std::cout << GridLogMessage << "Initial mom hamiltonian is " << Hmom << std::endl; + std::cout << GridLogMessage << "Final mom hamiltonian is " << Hmomprime << std::endl; + std::cout << GridLogMessage << "Delta mom hamiltonian is " << Hmomprime - Hmom << std::endl; + + std::cout << GridLogMessage << " S " << S << std::endl; + std::cout << GridLogMessage << " Sprime " << Sprime << std::endl; + std::cout << GridLogMessage << "dS (S' - S) :" << Sprime - S << std::endl; + std::cout << GridLogMessage << "predict dS (force) :" << dSpred << std::endl; + std::cout << GridLogMessage << "dSm " << dSm << std::endl; + std::cout << GridLogMessage << "dSm2" << dSm2 << std::endl; + + std::cout << GridLogMessage << "Total dS " << Hmomprime - Hmom + Sprime - S << std::endl; + + assert(fabs(real(Sprime - S - dSpred)) < 1.0); + + std::cout << GridLogMessage << "Done" << std::endl; + Grid_finalize(); +} diff --git a/tests/hadrons/Test_hadrons_2AS_spectrum.cc b/tests/hadrons/Test_hadrons_2AS_spectrum.cc new file mode 100644 index 00000000..2f519834 --- /dev/null +++ b/tests/hadrons/Test_hadrons_2AS_spectrum.cc @@ -0,0 +1,168 @@ +/******************************************************************************* + Grid physics library, www.github.com/paboyle/Grid + + Source file: tests/hadrons/Test_hadrons_spectrum.cc + + Copyright (C) 2015 + + 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. + *******************************************************************************/ + +#include + +using namespace Grid; +using namespace Hadrons; + + + BEGIN_HADRONS_NAMESPACE + BEGIN_MODULE_NAMESPACE(MFermion) + MODULE_REGISTER_NS(GaugeProp2AS, TGaugeProp, MFermion); + END_MODULE_NAMESPACE + BEGIN_MODULE_NAMESPACE(MSource) + MODULE_REGISTER_NS(Point2AS, TPoint, MSource); + END_MODULE_NAMESPACE + BEGIN_MODULE_NAMESPACE(MContraction) + MODULE_REGISTER_NS(Meson2AS, ARG(TMeson), MContraction); +// MODULE_REGISTER_NS(BaryonMultirep, ARG(TBaryon), MContraction); + END_MODULE_NAMESPACE + BEGIN_MODULE_NAMESPACE(MSink) + MODULE_REGISTER_NS(ScalarPoint2AS, TPoint, MSink); + END_MODULE_NAMESPACE + BEGIN_MODULE_NAMESPACE(MSolver) + MODULE_REGISTER_NS(RBPrecCG2AS, TRBPrecCG, MSolver); + END_MODULE_NAMESPACE + BEGIN_MODULE_NAMESPACE(MAction) + MODULE_REGISTER_NS(WilsonClover2AS, TWilsonClover, MAction); + END_MODULE_NAMESPACE + END_HADRONS_NAMESPACE + + +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 = {"l", "s"}; + std::vector mass = {-0.01, -0.04}; + double csw = 1.0; + // 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"); + MSource::Point2AS::Par ptPar; + ptPar.position = "0 0 0 0"; + application.createModule("pt", ptPar); + // sink + MSink::ScalarPoint2AS::Par sinkPar; + sinkPar.mom = "0 0 0"; + application.createModule("sink", sinkPar); + + // set fermion boundary conditions to be periodic space, antiperiodic time. + std::string boundary = "1 1 1 -1"; + + for (unsigned int i = 0; i < flavour.size(); ++i) + { + // actions + MAction::WilsonClover2AS::Par actionPar; + actionPar.gauge = "gauge"; + actionPar.mass = mass[i]; + actionPar.csw_r = csw; + actionPar.csw_t = csw; + actionPar.clover_anisotropy.isAnisotropic= false; + actionPar.clover_anisotropy.t_direction = Nd-1 ; + actionPar.clover_anisotropy.xi_0 = 1.0 ; + actionPar.clover_anisotropy.nu = 1.0 ; + actionPar.boundary = boundary; + application.createModule("WilsonClover2AS_" + flavour[i], actionPar); + + // solvers + MSolver::RBPrecCG2AS::Par solverPar; + solverPar.action = "WilsonClover2AS_" + flavour[i]; + solverPar.residual = 1.0e-8; + application.createModule("CG_" + flavour[i], + solverPar); + + // propagators + MFermion::GaugeProp2AS::Par quarkPar; + quarkPar.solver = "CG_" + flavour[i]; + quarkPar.source = "pt"; + application.createModule("Qpt_" + flavour[i], quarkPar); + quarkPar.source = "z2"; + application.createModule("QZ2_" + flavour[i], quarkPar); + } + for (unsigned int i = 0; i < flavour.size(); ++i) + for (unsigned int j = i; j < flavour.size(); ++j) + { + MContraction::Meson2AS::Par mesPar; + + mesPar.output = "mesons2AS/pt_" + flavour[i] + flavour[j]; + mesPar.q1 = "Qpt_" + flavour[i]; + mesPar.q2 = "Qpt_" + flavour[j]; + mesPar.gammas = "all"; + mesPar.sink = "sink"; + application.createModule("meson_pt_" + + flavour[i] + flavour[j], + mesPar); + + // mesPar.output = "mesons2AS/Z2_" + flavour[i] + flavour[j]; + // mesPar.q1 = "QZ2_" + flavour[i]; + // mesPar.q2 = "QZ2_" + flavour[j]; + // mesPar.gammas = "all"; + // mesPar.sink = "sink"; + // application.createModule("meson_Z2_" + // + flavour[i] + flavour[j], + // mesPar); + } + for (unsigned int i = 0; i < flavour.size(); ++i) + for (unsigned int j = i; j < flavour.size(); ++j) + for (unsigned int k = j; k < flavour.size(); ++k) + { + MContraction::Baryon::Par barPar; + + barPar.output = "baryons/pt_" + flavour[i] + flavour[j] + flavour[k]; + barPar.q1 = "Qpt_" + flavour[i]; + barPar.q2 = "Qpt_" + flavour[j]; + barPar.q3 = "Qpt_" + flavour[k]; + application.createModule( + "baryon_pt_" + flavour[i] + flavour[j] + flavour[k], barPar); + } + + // execution + application.saveParameterFile("spectrum.xml"); + application.run(); + + // epilogue + LOG(Message) << "Grid is finalizing now" << std::endl; + Grid_finalize(); + + return EXIT_SUCCESS; +} diff --git a/tests/hadrons/Test_hadrons_quark.cc b/tests/hadrons/Test_hadrons_quark.cc index eac065e9..cf3d2590 100644 --- a/tests/hadrons/Test_hadrons_quark.cc +++ b/tests/hadrons/Test_hadrons_quark.cc @@ -130,8 +130,8 @@ int main(int argc, char **argv) for (int c = 0; c < Nc; ++c) { ref = prop; - PropToFerm(ferm, prop, s, c); - FermToProp(prop, ferm, s, c); + PropToFerm(ferm, prop, s, c); + FermToProp(prop, ferm, s, c); std::cout << "Spin = " << s << ", Colour = " << c << std::endl; ref -= prop; diff --git a/tests/hadrons/Test_hadrons_wilsonFund.cc b/tests/hadrons/Test_hadrons_wilsonFund.cc new file mode 100644 index 00000000..083d3b8c --- /dev/null +++ b/tests/hadrons/Test_hadrons_wilsonFund.cc @@ -0,0 +1,157 @@ +/******************************************************************************* + Grid physics library, www.github.com/paboyle/Grid + + Source file: tests/hadrons/Test_hadrons_spectrum.cc + + Copyright (C) 2015 + + 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. + *******************************************************************************/ + +#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 = {"l"}; + std::vector mass = {-0.1}; + double csw = 0.0; + + // global parameters + Application::GlobalPar globalPar; + + globalPar.trajCounter.start = 309; + globalPar.trajCounter.end = 310; + globalPar.trajCounter.step = 1; + + globalPar.seed = "1 2 3 4"; + + application.setPar(globalPar); + // gauge field + application.createModule("gauge"); + + // sources + //MSource::Z2::Par z2Par; + //z2Par.tA = 0; + //z2Par.tB = 0; + //application.createModule("z2", z2Par); + MSource::Point::Par ptPar; + ptPar.position = "0 0 0 0"; + application.createModule("pt", ptPar); + // sink + MSink::Point::Par sinkPar; + sinkPar.mom = "0 0 0"; + application.createModule("sink", sinkPar); + + // set fermion boundary conditions to be periodic space, antiperiodic time. + std::string boundary = "1 1 1 -1"; + + for (unsigned int i = 0; i < flavour.size(); ++i) + { + // actions + MAction::WilsonClover::Par actionPar; + actionPar.gauge = "gauge"; + actionPar.mass = mass[i]; + actionPar.boundary = boundary; + actionPar.csw_r = csw; + actionPar.csw_t = csw; + + // !!!!! Check if Anisotropy works !!!!! + actionPar.clover_anisotropy.isAnisotropic= false; + actionPar.clover_anisotropy.t_direction = 3 ; // Explicit for D=4 + actionPar.clover_anisotropy.xi_0 = 1.0 ; + actionPar.clover_anisotropy.nu = 1.0 ; + + application.createModule("WilsonClover_" + flavour[i], actionPar); + + // solvers + MSolver::RBPrecCG::Par solverPar; + solverPar.action = "WilsonClover_" + 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); + // quarkPar.source = "z2"; + // application.createModule("QZ2_" + flavour[i], quarkPar); + } + for (unsigned int i = 0; i < flavour.size(); ++i) + for (unsigned int j = i; j < flavour.size(); ++j) + { + MContraction::Meson::Par mesPar; + + mesPar.output = "Fund_mesons/pt_" + flavour[i] + flavour[j]; + mesPar.q1 = "Qpt_" + flavour[i]; + mesPar.q2 = "Qpt_" + flavour[j]; + mesPar.gammas = "all"; + mesPar.sink = "sink"; + application.createModule("meson_pt_" + + flavour[i] + flavour[j], + mesPar); + // mesPar.output = "mesons/Z2_" + flavour[i] + flavour[j]; + // mesPar.q1 = "QZ2_" + flavour[i]; + // mesPar.q2 = "QZ2_" + flavour[j]; + // mesPar.gammas = "all"; + // mesPar.sink = "sink"; + // application.createModule("meson_Z2_" + // + flavour[i] + flavour[j], + // mesPar); + } + for (unsigned int i = 0; i < flavour.size(); ++i) + for (unsigned int j = i; j < flavour.size(); ++j) + for (unsigned int k = j; k < flavour.size(); ++k) + { + MContraction::Baryon::Par barPar; + + barPar.output = "Fund_baryons/pt_" + flavour[i] + flavour[j] + flavour[k]; + barPar.q1 = "Qpt_" + flavour[i]; + barPar.q2 = "Qpt_" + flavour[j]; + barPar.q3 = "Qpt_" + flavour[k]; + application.createModule( + "baryon_pt_" + flavour[i] + flavour[j] + flavour[k], barPar); + } + + // execution + application.saveParameterFile("WilsonClover_spectrum.xml"); + application.run(); + + // epilogue + LOG(Message) << "Grid is finalizing now" << std::endl; + Grid_finalize(); + + return EXIT_SUCCESS; +} diff --git a/tests/hmc/Test_hmc_EOWilsonCloverFermionGauge.cc b/tests/hmc/Test_hmc_EOWilsonCloverFermionGauge.cc new file mode 100644 index 00000000..fbb8a4a3 --- /dev/null +++ b/tests/hmc/Test_hmc_EOWilsonCloverFermionGauge.cc @@ -0,0 +1,139 @@ +/************************************************************************************* + +Grid physics library, www.github.com/paboyle/Grid + +Source file: ./tests/Test_hmc_WilsonFermionGauge.cc + +Copyright (C) 2016 + +Author: Guido Cossu + +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 + +int main(int argc, char **argv) { + using namespace Grid; + using namespace Grid::QCD; + + Grid_init(&argc, &argv); + int threads = GridThread::GetThreads(); + // here make a routine to print all the relevant information on the run + std::cout << GridLogMessage << "Grid is setup to use " << threads << " threads" << std::endl; + + // Typedefs to simplify notation + typedef GenericHMCRunner HMCWrapper; // Uses the default minimum norm + typedef WilsonImplR FermionImplPolicy; + typedef WilsonCloverFermionR FermionAction; + typedef typename FermionAction::FermionField FermionField; + + + //:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: + HMCWrapper TheHMC; + + // Grid from the command line + TheHMC.Resources.AddFourDimGrid("gauge"); + // Possibile to create the module by hand + // hardcoding parameters or using a Reader + + + // Checkpointer definition + CheckpointerParameters CPparams; + CPparams.config_prefix = "ckpoint_lat"; + CPparams.rng_prefix = "ckpoint_rng"; + CPparams.saveInterval = 5; + CPparams.format = "IEEE64BIG"; + + TheHMC.Resources.LoadNerscCheckpointer(CPparams); + + RNGModuleParameters RNGpar; + RNGpar.serial_seeds = "1 2 3 4 5"; + RNGpar.parallel_seeds = "6 7 8 9 10"; + TheHMC.Resources.SetRNGSeeds(RNGpar); + + // Construct observables + // here there is too much indirection + typedef PlaquetteMod PlaqObs; + TheHMC.Resources.AddObservable(); + ////////////////////////////////////////////// + + ///////////////////////////////////////////////////////////// + // Collect actions, here use more encapsulation + // need wrappers of the fermionic classes + // that have a complex construction + // standard + RealD beta = 5.6 ; + WilsonGaugeActionR Waction(beta); + + // temporarily need a gauge field + auto GridPtr = TheHMC.Resources.GetCartesian(); + auto GridRBPtr = TheHMC.Resources.GetRBCartesian(); + + LatticeGaugeField U(GridPtr); + + Real mass = 0.01; + Real csw = 1.0; + + FermionAction FermOp(U, *GridPtr, *GridRBPtr, mass, csw); + + ConjugateGradient CG(1.0e-8, 2000); + + TwoFlavourEvenOddPseudoFermionAction Nf2(FermOp, CG, CG); + + // Set smearing (true/false), default: false + Nf2.is_smeared = false; + + + // Collect actions + ActionLevel Level1(1); + Level1.push_back(&Nf2); + + ActionLevel Level2(4); + Level2.push_back(&Waction); + + TheHMC.TheAction.push_back(Level1); + TheHMC.TheAction.push_back(Level2); + ///////////////////////////////////////////////////////////// + + /* + double rho = 0.1; // smearing parameter + int Nsmear = 2; // number of smearing levels + Smear_Stout Stout(rho); + SmearedConfiguration SmearingPolicy( + UGrid, Nsmear, Stout); + */ + + // HMC parameters are serialisable + TheHMC.Parameters.MD.MDsteps = 20; + TheHMC.Parameters.MD.trajL = 1.0; + + TheHMC.ReadCommandLine(argc, argv); // these can be parameters from file + TheHMC.Run(); // no smearing + // TheHMC.Run(SmearingPolicy); // for smearing + + Grid_finalize(); + +} // main + + + + + + + diff --git a/tests/hmc/Test_hmc_ScalarActionNxN.cc b/tests/hmc/Test_hmc_ScalarActionNxN.cc index a4dad1a3..9e40beac 100644 --- a/tests/hmc/Test_hmc_ScalarActionNxN.cc +++ b/tests/hmc/Test_hmc_ScalarActionNxN.cc @@ -31,7 +31,8 @@ class ScalarActionParameters : Serializable { public: GRID_SERIALIZABLE_CLASS_MEMBERS(ScalarActionParameters, double, mass_squared, - double, lambda); + double, lambda, + double, g); template ScalarActionParameters(Reader& Reader){ @@ -140,7 +141,7 @@ int main(int argc, char **argv) { // Scalar action in adjoint representation ScalarActionParameters SPar(Reader); - ScalarAction Saction(SPar.mass_squared, SPar.lambda); + ScalarAction Saction(SPar.mass_squared, SPar.lambda, SPar.g); // Collect actions ActionLevel> Level1(1); diff --git a/tests/hmc/Test_hmc_WC2ASFG_Production.cc b/tests/hmc/Test_hmc_WC2ASFG_Production.cc new file mode 100644 index 00000000..d255ab5d --- /dev/null +++ b/tests/hmc/Test_hmc_WC2ASFG_Production.cc @@ -0,0 +1,213 @@ +/************************************************************************************* + +Grid physics library, www.github.com/paboyle/Grid + +Source file: ./tests/Test_hmc_WilsonFermionGauge.cc + +Copyright (C) 2017 + +Author: Guido Cossu + +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 + + +namespace Grid{ + struct FermionParameters: Serializable { + GRID_SERIALIZABLE_CLASS_MEMBERS(FermionParameters, + double, mass, + double, csw, + double, StoppingCondition, + int, MaxCGIterations, + bool, ApplySmearing); + }; + + + struct WilsonCloverHMCParameters: Serializable { + GRID_SERIALIZABLE_CLASS_MEMBERS(WilsonCloverHMCParameters, + double, gauge_beta, + FermionParameters, WilsonClover) + + template + WilsonCloverHMCParameters(Reader& Reader){ + read(Reader, "Action", *this); + } + }; + + struct SmearingParameters: Serializable { + GRID_SERIALIZABLE_CLASS_MEMBERS(SmearingParameters, + double, rho, + Integer, Nsmear) + + template + SmearingParameters(Reader& Reader){ + read(Reader, "StoutSmearing", *this); + } + + }; + + +} + +int main(int argc, char **argv) +{ + using namespace Grid; + using namespace Grid::QCD; + + typedef Representations< FundamentalRepresentation, TwoIndexAntiSymmetricRepresentation > TheRepresentations; + + Grid_init(&argc, &argv); + int threads = GridThread::GetThreads(); + // here make a routine to print all the relevant information on the run + std::cout << GridLogMessage << "Grid is setup to use " << threads << " threads" << std::endl; + + // Typedefs to simplify notation + typedef GenericHMCRunnerHirep HMCWrapper; // Uses the default minimum norm + typedef WilsonTwoIndexAntiSymmetricImplR FermionImplPolicy; // gauge field implemetation for the pseudofermions + typedef WilsonCloverTwoIndexAntiSymmetricFermionR FermionAction; // type of lattice fermions (Wilson, DW, ...) + typedef typename FermionAction::FermionField FermionField; + typedef Grid::JSONReader Serialiser; + + //:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: + HMCWrapper TheHMC; + + // Grid from the command line + TheHMC.ReadCommandLine(argc, argv); + if (TheHMC.ParameterFile.empty()){ + std::cout << "Input file not specified." + << "Use --ParameterFile option in the command line.\nAborting" + << std::endl; + exit(1); + } + Serialiser Reader(TheHMC.ParameterFile); + WilsonCloverHMCParameters MyParams(Reader); + + // Apply smearing to the fermionic action + bool ApplySmearing = MyParams.WilsonClover.ApplySmearing; + + TheHMC.Resources.AddFourDimGrid("gauge"); + + // Checkpointer definition + CheckpointerParameters CPparams(Reader); + + /* + CPparams.config_prefix = "ckpoint_lat"; + CPparams.rng_prefix = "ckpoint_rng"; + CPparams.saveInterval = 5; + CPparams.format = "IEEE64BIG"; + */ + + TheHMC.Resources.LoadNerscCheckpointer(CPparams); + + RNGModuleParameters RNGpar(Reader); + /* + RNGpar.serial_seeds = "1 2 3 4 5"; + RNGpar.parallel_seeds = "6 7 8 9 10"; + TheHMC.Resources.SetRNGSeeds(RNGpar); + */ + TheHMC.Resources.SetRNGSeeds(RNGpar); + + // Construct observables + typedef PlaquetteMod PlaqObs; + TheHMC.Resources.AddObservable(); + + typedef PolyakovMod PolyakovObs; + TheHMC.Resources.AddObservable(); + + //typedef TopologicalChargeMod QObs; + //TopologyObsParameters TopParams(Reader); + //TheHMC.Resources.AddObservable(TopParams); + ////////////////////////////////////////////// + + ///////////////////////////////////////////////////////////// + // Collect actions, here use more encapsulation + // need wrappers of the fermionic classes + // that have a complex construction + // standard + + //RealD beta = 5.6; + WilsonGaugeActionR Waction(MyParams.gauge_beta); + + auto GridPtr = TheHMC.Resources.GetCartesian(); + auto GridRBPtr = TheHMC.Resources.GetRBCartesian(); + + // temporarily need a gauge field + TwoIndexAntiSymmetricRepresentation::LatticeField U(GridPtr); + + //Real mass = 0.01; + //Real csw = 1.0; + + Real mass = MyParams.WilsonClover.mass; + Real csw = MyParams.WilsonClover.csw; + + std::cout << "mass and csw" << mass << " and " << csw << std::endl; + + FermionAction FermOp(U, *GridPtr, *GridRBPtr, mass, csw, csw); + ConjugateGradient CG(MyParams.WilsonClover.StoppingCondition, MyParams.WilsonClover.MaxCGIterations); + TwoFlavourPseudoFermionAction Nf2(FermOp, CG, CG); + + // Set smearing (true/false), default: false + Nf2.is_smeared = ApplySmearing; + + // Collect actions + ActionLevel Level1(1); + Level1.push_back(&Nf2); + + ActionLevel Level2(4); + Level2.push_back(&Waction); + + TheHMC.TheAction.push_back(Level1); + TheHMC.TheAction.push_back(Level2); + ///////////////////////////////////////////////////////////// + + + /* + double rho = 0.1; // smearing parameter + int Nsmear = 2; // number of smearing levels + Smear_Stout Stout(rho); + SmearedConfiguration SmearingPolicy( + UGrid, Nsmear, Stout); + */ + + // HMC parameters are serialisable + + TheHMC.Parameters.initialize(Reader); + //TheHMC.Parameters.MD.MDsteps = 20; + //TheHMC.Parameters.MD.trajL = 1.0; + + if (ApplySmearing){ + SmearingParameters SmPar(Reader); + //double rho = 0.1; // smearing parameter + //int Nsmear = 3; // number of smearing levels + Smear_Stout Stout(SmPar.rho); + SmearedConfiguration SmearingPolicy(GridPtr, SmPar.Nsmear, Stout); + TheHMC.Run(SmearingPolicy); // for smearing + } else { + TheHMC.Run(); // no smearing + } + + //TheHMC.ReadCommandLine(argc, argv); // these can be parameters from file + //TheHMC.Run(); // no smearing + // TheHMC.Run(SmearingPolicy); // for smearing + + Grid_finalize(); + +} // main + diff --git a/tests/hmc/Test_hmc_WC2SFG_Production.cc b/tests/hmc/Test_hmc_WC2SFG_Production.cc new file mode 100644 index 00000000..8d5fc458 --- /dev/null +++ b/tests/hmc/Test_hmc_WC2SFG_Production.cc @@ -0,0 +1,212 @@ +/************************************************************************************* + +Grid physics library, www.github.com/paboyle/Grid + +Source file: ./tests/Test_hmc_WilsonFermionGauge.cc + +Copyright (C) 2017 + +Author: Guido Cossu + +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 + + +namespace Grid{ + struct FermionParameters: Serializable { + GRID_SERIALIZABLE_CLASS_MEMBERS(FermionParameters, + double, mass, + double, csw, + double, StoppingCondition, + int, MaxCGIterations, + bool, ApplySmearing); + }; + + + struct WilsonCloverHMCParameters: Serializable { + GRID_SERIALIZABLE_CLASS_MEMBERS(WilsonCloverHMCParameters, + double, gauge_beta, + FermionParameters, WilsonClover) + + template + WilsonCloverHMCParameters(Reader& Reader){ + read(Reader, "Action", *this); + } + }; + + struct SmearingParameters: Serializable { + GRID_SERIALIZABLE_CLASS_MEMBERS(SmearingParameters, + double, rho, + Integer, Nsmear) + + template + SmearingParameters(Reader& Reader){ + read(Reader, "StoutSmearing", *this); + } + + }; + + +} + +int main(int argc, char **argv) +{ + using namespace Grid; + using namespace Grid::QCD; + + typedef Representations< FundamentalRepresentation, TwoIndexSymmetricRepresentation > TheRepresentations; + + Grid_init(&argc, &argv); + int threads = GridThread::GetThreads(); + // here make a routine to print all the relevant information on the run + std::cout << GridLogMessage << "Grid is setup to use " << threads << " threads" << std::endl; + + // Typedefs to simplify notation + typedef GenericHMCRunnerHirep HMCWrapper; // Uses the default minimum norm + typedef WilsonTwoIndexSymmetricImplR FermionImplPolicy; // gauge field implemetation for the pseudofermions + typedef WilsonCloverTwoIndexSymmetricFermionR FermionAction; // type of lattice fermions (Wilson, DW, ...) + typedef typename FermionAction::FermionField FermionField; + typedef Grid::JSONReader Serialiser; + + //:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: + HMCWrapper TheHMC; + + // Grid from the command line + TheHMC.ReadCommandLine(argc, argv); + if (TheHMC.ParameterFile.empty()){ + std::cout << "Input file not specified." + << "Use --ParameterFile option in the command line.\nAborting" + << std::endl; + exit(1); + } + Serialiser Reader(TheHMC.ParameterFile); + WilsonCloverHMCParameters MyParams(Reader); + + // Apply smearing to the fermionic action + bool ApplySmearing = MyParams.WilsonClover.ApplySmearing; + + TheHMC.Resources.AddFourDimGrid("gauge"); + + // Checkpointer definition + CheckpointerParameters CPparams(Reader); + + /* + CPparams.config_prefix = "ckpoint_lat"; + CPparams.rng_prefix = "ckpoint_rng"; + CPparams.saveInterval = 5; + CPparams.format = "IEEE64BIG"; + */ + + TheHMC.Resources.LoadNerscCheckpointer(CPparams); + + RNGModuleParameters RNGpar(Reader); + /* + RNGpar.serial_seeds = "1 2 3 4 5"; + RNGpar.parallel_seeds = "6 7 8 9 10"; + TheHMC.Resources.SetRNGSeeds(RNGpar); + */ + TheHMC.Resources.SetRNGSeeds(RNGpar); + + // Construct observables + typedef PlaquetteMod PlaqObs; + TheHMC.Resources.AddObservable(); + + typedef PolyakovMod PolyakovObs; + TheHMC.Resources.AddObservable(); + + //typedef TopologicalChargeMod QObs; + //TopologyObsParameters TopParams(Reader); + //TheHMC.Resources.AddObservable(TopParams); + ////////////////////////////////////////////// + + ///////////////////////////////////////////////////////////// + // Collect actions, here use more encapsulation + // need wrappers of the fermionic classes + // that have a complex construction + // standard + + //RealD beta = 5.6; + WilsonGaugeActionR Waction(MyParams.gauge_beta); + + auto GridPtr = TheHMC.Resources.GetCartesian(); + auto GridRBPtr = TheHMC.Resources.GetRBCartesian(); + + // temporarily need a gauge field + TwoIndexSymmetricRepresentation::LatticeField U(GridPtr); + + //Real mass = 0.01; + //Real csw = 1.0; + + Real mass = MyParams.WilsonClover.mass; + Real csw = MyParams.WilsonClover.csw; + + std::cout << "mass and csw" << mass << " and " << csw << std::endl; + + FermionAction FermOp(U, *GridPtr, *GridRBPtr, mass, csw, csw); + ConjugateGradient CG(MyParams.WilsonClover.StoppingCondition, MyParams.WilsonClover.MaxCGIterations); + TwoFlavourPseudoFermionAction Nf2(FermOp, CG, CG); + + // Set smearing (true/false), default: false + Nf2.is_smeared = ApplySmearing; + + // Collect actions + ActionLevel Level1(1); + Level1.push_back(&Nf2); + + ActionLevel Level2(4); + Level2.push_back(&Waction); + + TheHMC.TheAction.push_back(Level1); + TheHMC.TheAction.push_back(Level2); + ///////////////////////////////////////////////////////////// + + + /* + double rho = 0.1; // smearing parameter + int Nsmear = 2; // number of smearing levels + Smear_Stout Stout(rho); + SmearedConfiguration SmearingPolicy( + UGrid, Nsmear, Stout); + */ + + // HMC parameters are serialisable + + TheHMC.Parameters.initialize(Reader); + //TheHMC.Parameters.MD.MDsteps = 20; + //TheHMC.Parameters.MD.trajL = 1.0; + + if (ApplySmearing){ + SmearingParameters SmPar(Reader); + //double rho = 0.1; // smearing parameter + //int Nsmear = 3; // number of smearing levels + Smear_Stout Stout(SmPar.rho); + SmearedConfiguration SmearingPolicy(GridPtr, SmPar.Nsmear, Stout); + TheHMC.Run(SmearingPolicy); // for smearing + } else { + TheHMC.Run(); // no smearing + } + + //TheHMC.ReadCommandLine(argc, argv); // these can be parameters from file + //TheHMC.Run(); // no smearing + // TheHMC.Run(SmearingPolicy); // for smearing + + Grid_finalize(); + +} // main diff --git a/tests/hmc/Test_hmc_WCFG_Production.cc b/tests/hmc/Test_hmc_WCFG_Production.cc new file mode 100644 index 00000000..895e4f81 --- /dev/null +++ b/tests/hmc/Test_hmc_WCFG_Production.cc @@ -0,0 +1,210 @@ +/************************************************************************************* + +Grid physics library, www.github.com/paboyle/Grid + +Source file: ./tests/Test_hmc_WilsonFermionGauge.cc + +Copyright (C) 2017 + +Author: Guido Cossu + +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 + + +namespace Grid{ + struct FermionParameters: Serializable { + GRID_SERIALIZABLE_CLASS_MEMBERS(FermionParameters, + double, mass, + double, csw, + double, StoppingCondition, + int, MaxCGIterations, + bool, ApplySmearing); + }; + + + struct WilsonCloverHMCParameters: Serializable { + GRID_SERIALIZABLE_CLASS_MEMBERS(WilsonCloverHMCParameters, + double, gauge_beta, + FermionParameters, WilsonClover) + + template + WilsonCloverHMCParameters(Reader& Reader){ + read(Reader, "Action", *this); + } + }; + + struct SmearingParameters: Serializable { + GRID_SERIALIZABLE_CLASS_MEMBERS(SmearingParameters, + double, rho, + Integer, Nsmear) + + template + SmearingParameters(Reader& Reader){ + read(Reader, "StoutSmearing", *this); + } + + }; + + +} + +int main(int argc, char **argv) +{ + using namespace Grid; + using namespace Grid::QCD; + + Grid_init(&argc, &argv); + int threads = GridThread::GetThreads(); + // here make a routine to print all the relevant information on the run + std::cout << GridLogMessage << "Grid is setup to use " << threads << " threads" << std::endl; + + // Typedefs to simplify notation + typedef GenericHMCRunner HMCWrapper; // Uses the default minimum norm + typedef WilsonImplR FermionImplPolicy; + typedef WilsonCloverFermionR FermionAction; + typedef typename FermionAction::FermionField FermionField; + typedef Grid::JSONReader Serialiser; + + //:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: + HMCWrapper TheHMC; + + // Grid from the command line + TheHMC.ReadCommandLine(argc, argv); + if (TheHMC.ParameterFile.empty()){ + std::cout << "Input file not specified." + << "Use --ParameterFile option in the command line.\nAborting" + << std::endl; + exit(1); + } + Serialiser Reader(TheHMC.ParameterFile); + WilsonCloverHMCParameters MyParams(Reader); + + // Apply smearing to the fermionic action + bool ApplySmearing = MyParams.WilsonClover.ApplySmearing; + + TheHMC.Resources.AddFourDimGrid("gauge"); + + // Checkpointer definition + CheckpointerParameters CPparams(Reader); + + /* + CPparams.config_prefix = "ckpoint_lat"; + CPparams.rng_prefix = "ckpoint_rng"; + CPparams.saveInterval = 5; + CPparams.format = "IEEE64BIG"; + */ + + TheHMC.Resources.LoadNerscCheckpointer(CPparams); + + RNGModuleParameters RNGpar(Reader); + /* + RNGpar.serial_seeds = "1 2 3 4 5"; + RNGpar.parallel_seeds = "6 7 8 9 10"; + TheHMC.Resources.SetRNGSeeds(RNGpar); + */ + TheHMC.Resources.SetRNGSeeds(RNGpar); + + // Construct observables + typedef PlaquetteMod PlaqObs; + TheHMC.Resources.AddObservable(); + + typedef PolyakovMod PolyakovObs; + TheHMC.Resources.AddObservable(); + + //typedef TopologicalChargeMod QObs; + //TopologyObsParameters TopParams(Reader); + //TheHMC.Resources.AddObservable(TopParams); + ////////////////////////////////////////////// + + ///////////////////////////////////////////////////////////// + // Collect actions, here use more encapsulation + // need wrappers of the fermionic classes + // that have a complex construction + // standard + + //RealD beta = 5.6; + WilsonGaugeActionR Waction(MyParams.gauge_beta); + + auto GridPtr = TheHMC.Resources.GetCartesian(); + auto GridRBPtr = TheHMC.Resources.GetRBCartesian(); + + // temporarily need a gauge field + LatticeGaugeField U(GridPtr); + + //Real mass = 0.01; + //Real csw = 1.0; + + Real mass = MyParams.WilsonClover.mass; + Real csw = MyParams.WilsonClover.csw; + + std::cout << "mass and csw" << mass << " and " << csw << std::endl; + + FermionAction FermOp(U, *GridPtr, *GridRBPtr, mass, csw, csw); + ConjugateGradient CG(MyParams.WilsonClover.StoppingCondition, MyParams.WilsonClover.MaxCGIterations); + TwoFlavourPseudoFermionAction Nf2(FermOp, CG, CG); + + // Set smearing (true/false), default: false + Nf2.is_smeared = ApplySmearing; + + // Collect actions + ActionLevel Level1(1); + Level1.push_back(&Nf2); + + ActionLevel Level2(4); + Level2.push_back(&Waction); + + TheHMC.TheAction.push_back(Level1); + TheHMC.TheAction.push_back(Level2); + ///////////////////////////////////////////////////////////// + + + /* + double rho = 0.1; // smearing parameter + int Nsmear = 2; // number of smearing levels + Smear_Stout Stout(rho); + SmearedConfiguration SmearingPolicy( + UGrid, Nsmear, Stout); + */ + + // HMC parameters are serialisable + + TheHMC.Parameters.initialize(Reader); + //TheHMC.Parameters.MD.MDsteps = 20; + //TheHMC.Parameters.MD.trajL = 1.0; + + if (ApplySmearing){ + SmearingParameters SmPar(Reader); + //double rho = 0.1; // smearing parameter + //int Nsmear = 3; // number of smearing levels + Smear_Stout Stout(SmPar.rho); + SmearedConfiguration SmearingPolicy(GridPtr, SmPar.Nsmear, Stout); + TheHMC.Run(SmearingPolicy); // for smearing + } else { + TheHMC.Run(); // no smearing + } + + //TheHMC.ReadCommandLine(argc, argv); // these can be parameters from file + //TheHMC.Run(); // no smearing + // TheHMC.Run(SmearingPolicy); // for smearing + + Grid_finalize(); + +} // main diff --git a/tests/hmc/Test_hmc_WCMixedRepFG_Production.cc b/tests/hmc/Test_hmc_WCMixedRepFG_Production.cc new file mode 100644 index 00000000..aa5cce85 --- /dev/null +++ b/tests/hmc/Test_hmc_WCMixedRepFG_Production.cc @@ -0,0 +1,224 @@ +/************************************************************************************* + +Grid physics library, www.github.com/paboyle/Grid + +Source file: ./tests/Test_hmc_WilsonAdjointFermionGauge.cc + +Copyright (C) 2015 + +Author: Peter Boyle +Author: Peter Boyle +Author: neo +Author: paboyle + +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 "Grid/Grid.h" + + +namespace Grid{ + struct FermionParameters: Serializable { + GRID_SERIALIZABLE_CLASS_MEMBERS(FermionParameters, + double, mass, + double, csw, + double, StoppingCondition, + int, MaxCGIterations, + bool, ApplySmearing); + }; + + struct WilsonCloverHMCParameters: Serializable { + GRID_SERIALIZABLE_CLASS_MEMBERS(WilsonCloverHMCParameters, + double, gauge_beta, + FermionParameters, WilsonCloverFund, + FermionParameters, WilsonCloverAS) + + template + WilsonCloverHMCParameters(Reader& Reader){ + read(Reader, "Action", *this); + } + }; + + struct SmearingParameters: Serializable { + GRID_SERIALIZABLE_CLASS_MEMBERS(SmearingParameters, + double, rho, + Integer, Nsmear) + + template + SmearingParameters(Reader& Reader){ + read(Reader, "StoutSmearing", *this); + } + + }; +} + + +int main(int argc, char **argv) { + using namespace Grid; + using namespace Grid::QCD; + + // Here change the allowed (higher) representations + typedef Representations< FundamentalRepresentation, TwoIndexAntiSymmetricRepresentation> TheRepresentations; + + Grid_init(&argc, &argv); + int threads = GridThread::GetThreads(); + // here make a routine to print all the relevant information on the run + std::cout << GridLogMessage << "Grid is setup to use " << threads << " threads" << std::endl; + + // Typedefs to simplify notation + typedef GenericHMCRunnerHirep HMCWrapper; + + typedef WilsonImplR FundImplPolicy; + typedef WilsonCloverFermionR FundFermionAction; + typedef typename FundFermionAction::FermionField FundFermionField; + + typedef WilsonTwoIndexAntiSymmetricImplR ASymmImplPolicy; + typedef WilsonCloverTwoIndexAntiSymmetricFermionR ASymmFermionAction; + typedef typename ASymmFermionAction::FermionField ASymmFermionField; + + typedef Grid::JSONReader Serialiser; + //:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: + HMCWrapper TheHMC; + + // Grid from the command line + TheHMC.ReadCommandLine(argc, argv); + if (TheHMC.ParameterFile.empty()){ + std::cout << "Input file not specified." + << "Use --ParameterFile option in the command line.\nAborting" + << std::endl; + exit(1); + } + Serialiser Reader(TheHMC.ParameterFile); + WilsonCloverHMCParameters MyParams(Reader); + + // Apply smearing to the fermionic action + bool ApplySmearingFund = MyParams.WilsonCloverFund.ApplySmearing; + bool ApplySmearingAS = MyParams.WilsonCloverAS.ApplySmearing; + + + TheHMC.Resources.AddFourDimGrid("gauge"); + + // Checkpointer definition + CheckpointerParameters CPparams(Reader); + + /* + CPparams.config_prefix = "ckpoint_lat"; + CPparams.rng_prefix = "ckpoint_rng"; + CPparams.saveInterval = 5; + CPparams.format = "IEEE64BIG"; + */ + + TheHMC.Resources.LoadNerscCheckpointer(CPparams); + + RNGModuleParameters RNGpar(Reader); + /* + RNGpar.serial_seeds = "1 2 3 4 5"; + RNGpar.parallel_seeds = "6 7 8 9 10"; + TheHMC.Resources.SetRNGSeeds(RNGpar); + */ + TheHMC.Resources.SetRNGSeeds(RNGpar); + + // Construct observables + typedef PlaquetteMod PlaqObs; + TheHMC.Resources.AddObservable(); + + typedef PolyakovMod PolyakovObs; + TheHMC.Resources.AddObservable(); + + typedef TopologicalChargeMod QObs; + TopologyObsParameters TopParams(Reader); + TheHMC.Resources.AddObservable(TopParams); + ////////////////////////////////////////////// + + ///////////////////////////////////////////////////////////// + // Collect actions, here use more encapsulation + // need wrappers of the fermionic classes + // that have a complex construction + // standard + + //RealD beta = 5.6; + WilsonGaugeActionR Waction(MyParams.gauge_beta); + + auto GridPtr = TheHMC.Resources.GetCartesian(); + auto GridRBPtr = TheHMC.Resources.GetRBCartesian(); + + // temporarily need a gauge field + FundamentalRepresentation::LatticeField UF(GridPtr); + TwoIndexAntiSymmetricRepresentation::LatticeField UAS(GridPtr); + + + Real Fundmass = MyParams.WilsonCloverFund.mass; + Real Fundcsw = MyParams.WilsonCloverFund.csw; + Real ASmass = MyParams.WilsonCloverAS.mass; + Real AScsw = MyParams.WilsonCloverAS.csw; + + + + std::cout << "Fund: mass and csw" << Fundmass << " and " << Fundcsw << std::endl; + std::cout << "AS : mass and csw" << ASmass << " and " << AScsw << std::endl; + + + FundFermionAction FundFermOp(UF, *GridPtr, *GridRBPtr, Fundmass, Fundcsw, Fundcsw); + ConjugateGradient CG_Fund(MyParams.WilsonCloverFund.StoppingCondition, MyParams.WilsonCloverFund.MaxCGIterations); + TwoFlavourPseudoFermionAction Nf2_Fund(FundFermOp, CG_Fund, CG_Fund); + + ASymmFermionAction ASFermOp(UAS, *GridPtr, *GridRBPtr, ASmass, AScsw, AScsw); + ConjugateGradient CG_AS(MyParams.WilsonCloverAS.StoppingCondition, MyParams.WilsonCloverAS.MaxCGIterations); + TwoFlavourPseudoFermionAction Nf2_AS(ASFermOp, CG_AS, CG_AS); + + Nf2_Fund.is_smeared = ApplySmearingFund; + Nf2_AS.is_smeared = ApplySmearingAS; + + + // Collect actions + ActionLevel Level1(1); + Level1.push_back(&Nf2_Fund); + Level1.push_back(&Nf2_AS); + + + ActionLevel Level2(4); + Level2.push_back(&Waction); + + TheHMC.TheAction.push_back(Level1); + TheHMC.TheAction.push_back(Level2); + + TheHMC.Parameters.initialize(Reader); + //TheHMC.Parameters.MD.MDsteps = 20; + //TheHMC.Parameters.MD.trajL = 1.0; +/* + if (ApplySmearingFund || ApplySmearingAS){ + SmearingParameters SmPar(Reader); + //double rho = 0.1; // smearing parameter + //int Nsmear = 3; // number of smearing levels + Smear_Stout Stout(SmPar.rho); + SmearedConfiguration SmearingPolicy(GridPtr, SmPar.Nsmear, Stout); + TheHMC.Run(SmearingPolicy); // for smearing + } else { + TheHMC.Run(); // no smearing + } +*/ + TheHMC.Run(); + + + //TheHMC.ReadCommandLine(argc, argv); // these can be parameters from file + //TheHMC.Run(); // no smearing + // TheHMC.Run(SmearingPolicy); // for smearing + + Grid_finalize(); + +} // main diff --git a/tests/hmc/Test_hmc_WCadjFG_Production.cc b/tests/hmc/Test_hmc_WCadjFG_Production.cc new file mode 100644 index 00000000..48cea756 --- /dev/null +++ b/tests/hmc/Test_hmc_WCadjFG_Production.cc @@ -0,0 +1,213 @@ +/************************************************************************************* + +Grid physics library, www.github.com/paboyle/Grid + +Source file: ./tests/Test_hmc_WilsonFermionGauge.cc + +Copyright (C) 2017 + +Author: Guido Cossu + +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 + + +namespace Grid{ + struct FermionParameters: Serializable { + GRID_SERIALIZABLE_CLASS_MEMBERS(FermionParameters, + double, mass, + double, csw, + double, StoppingCondition, + int, MaxCGIterations, + bool, ApplySmearing); + }; + + + struct WilsonCloverHMCParameters: Serializable { + GRID_SERIALIZABLE_CLASS_MEMBERS(WilsonCloverHMCParameters, + double, gauge_beta, + FermionParameters, WilsonClover) + + template + WilsonCloverHMCParameters(Reader& Reader){ + read(Reader, "Action", *this); + } + }; + + struct SmearingParameters: Serializable { + GRID_SERIALIZABLE_CLASS_MEMBERS(SmearingParameters, + double, rho, + Integer, Nsmear) + + template + SmearingParameters(Reader& Reader){ + read(Reader, "StoutSmearing", *this); + } + + }; + + +} + +int main(int argc, char **argv) +{ + using namespace Grid; + using namespace Grid::QCD; + + typedef Representations< FundamentalRepresentation, AdjointRepresentation > TheRepresentations; + + Grid_init(&argc, &argv); + int threads = GridThread::GetThreads(); + // here make a routine to print all the relevant information on the run + std::cout << GridLogMessage << "Grid is setup to use " << threads << " threads" << std::endl; + + // Typedefs to simplify notation + typedef GenericHMCRunnerHirep HMCWrapper; // Uses the default minimum norm + typedef WilsonAdjImplR FermionImplPolicy; // gauge field implemetation for the pseudofermions + typedef WilsonCloverAdjFermionR FermionAction; // type of lattice fermions (Wilson, DW, ...) + typedef typename FermionAction::FermionField FermionField; + typedef Grid::JSONReader Serialiser; + + //:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: + HMCWrapper TheHMC; + + // Grid from the command line + TheHMC.ReadCommandLine(argc, argv); + if (TheHMC.ParameterFile.empty()){ + std::cout << "Input file not specified." + << "Use --ParameterFile option in the command line.\nAborting" + << std::endl; + exit(1); + } + Serialiser Reader(TheHMC.ParameterFile); + WilsonCloverHMCParameters MyParams(Reader); + + // Apply smearing to the fermionic action + bool ApplySmearing = MyParams.WilsonClover.ApplySmearing; + + TheHMC.Resources.AddFourDimGrid("gauge"); + + // Checkpointer definition + CheckpointerParameters CPparams(Reader); + + /* + CPparams.config_prefix = "ckpoint_lat"; + CPparams.rng_prefix = "ckpoint_rng"; + CPparams.saveInterval = 5; + CPparams.format = "IEEE64BIG"; + */ + + TheHMC.Resources.LoadNerscCheckpointer(CPparams); + + RNGModuleParameters RNGpar(Reader); + /* + RNGpar.serial_seeds = "1 2 3 4 5"; + RNGpar.parallel_seeds = "6 7 8 9 10"; + TheHMC.Resources.SetRNGSeeds(RNGpar); + */ + TheHMC.Resources.SetRNGSeeds(RNGpar); + + // Construct observables + typedef PlaquetteMod PlaqObs; + TheHMC.Resources.AddObservable(); + + typedef PolyakovMod PolyakovObs; + TheHMC.Resources.AddObservable(); + + typedef TopologicalChargeMod QObs; + TopologyObsParameters TopParams(Reader); + TheHMC.Resources.AddObservable(TopParams); + ////////////////////////////////////////////// + + ///////////////////////////////////////////////////////////// + // Collect actions, here use more encapsulation + // need wrappers of the fermionic classes + // that have a complex construction + // standard + + //RealD beta = 5.6; + WilsonGaugeActionR Waction(MyParams.gauge_beta); + + auto GridPtr = TheHMC.Resources.GetCartesian(); + auto GridRBPtr = TheHMC.Resources.GetRBCartesian(); + + // temporarily need a gauge field + AdjointRepresentation::LatticeField U(GridPtr); + + //Real mass = 0.01; + //Real csw = 1.0; + + Real mass = MyParams.WilsonClover.mass; + Real csw = MyParams.WilsonClover.csw; + + std::cout << "mass and csw" << mass << " and " << csw << std::endl; + + FermionAction FermOp(U, *GridPtr, *GridRBPtr, mass, csw, csw); + ConjugateGradient CG(MyParams.WilsonClover.StoppingCondition, MyParams.WilsonClover.MaxCGIterations); + TwoFlavourPseudoFermionAction Nf2(FermOp, CG, CG); + + // Set smearing (true/false), default: false + Nf2.is_smeared = ApplySmearing; + + // Collect actions + ActionLevel Level1(1); + Level1.push_back(&Nf2); + + ActionLevel Level2(4); + Level2.push_back(&Waction); + + TheHMC.TheAction.push_back(Level1); + TheHMC.TheAction.push_back(Level2); + ///////////////////////////////////////////////////////////// + + + /* + double rho = 0.1; // smearing parameter + int Nsmear = 2; // number of smearing levels + Smear_Stout Stout(rho); + SmearedConfiguration SmearingPolicy( + UGrid, Nsmear, Stout); + */ + + // HMC parameters are serialisable + + TheHMC.Parameters.initialize(Reader); + //TheHMC.Parameters.MD.MDsteps = 20; + //TheHMC.Parameters.MD.trajL = 1.0; + + if (ApplySmearing){ + SmearingParameters SmPar(Reader); + //double rho = 0.1; // smearing parameter + //int Nsmear = 3; // number of smearing levels + Smear_Stout Stout(SmPar.rho); + SmearedConfiguration SmearingPolicy(GridPtr, SmPar.Nsmear, Stout); + TheHMC.Run(SmearingPolicy); // for smearing + } else { + TheHMC.Run(); // no smearing + } + + //TheHMC.ReadCommandLine(argc, argv); // these can be parameters from file + //TheHMC.Run(); // no smearing + // TheHMC.Run(SmearingPolicy); // for smearing + + Grid_finalize(); + +} // main + diff --git a/tests/hmc/Test_hmc_WG_Production.cc b/tests/hmc/Test_hmc_WG_Production.cc new file mode 100644 index 00000000..b99446d5 --- /dev/null +++ b/tests/hmc/Test_hmc_WG_Production.cc @@ -0,0 +1,114 @@ +/************************************************************************************* + +Grid physics library, www.github.com/paboyle/Grid + +Source file: ./tests/Test_hmc_WilsonFermionGauge.cc + +Copyright (C) 2015 + +Author: Guido Cossu + +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 + +namespace Grid{ + struct ActionParameters: Serializable { + GRID_SERIALIZABLE_CLASS_MEMBERS(ActionParameters, + double, beta) + + + template + ActionParameters(Reader& Reader){ + read(Reader, "Action", *this); + } + + }; + +} + + +int main(int argc, char **argv) { + using namespace Grid; + using namespace Grid::QCD; + + Grid_init(&argc, &argv); + GridLogLayout(); + + // Typedefs to simplify notation + typedef GenericHMCRunner HMCWrapper; // Uses the default minimum norm + HMCWrapper TheHMC; + typedef Grid::JSONReader Serialiser; + + // Grid from the command line + TheHMC.Resources.AddFourDimGrid("gauge"); + TheHMC.ReadCommandLine(argc, argv); // these can be parameters from file + // Reader, file should come from command line + if (TheHMC.ParameterFile.empty()){ + std::cout << "Input file not specified." + << "Use --ParameterFile option in the command line.\nAborting" + << std::endl; + exit(1); + } + Serialiser Reader(TheHMC.ParameterFile); + + + + // Checkpointer definition + CheckpointerParameters CPparams(Reader); + TheHMC.Resources.LoadNerscCheckpointer(CPparams); + + RNGModuleParameters RNGpar(Reader); + TheHMC.Resources.SetRNGSeeds(RNGpar); + + // Construct observables + // here there is too much indirection + typedef PlaquetteMod PlaqObs; + typedef TopologicalChargeMod QObs; + TheHMC.Resources.AddObservable(); + TopologyObsParameters TopParams(Reader); + TheHMC.Resources.AddObservable(TopParams); + ////////////////////////////////////////////// + + ///////////////////////////////////////////////////////////// + // Collect actions, here use more encapsulation + // need wrappers of the fermionic classes + // that have a complex construction + // standard + ActionParameters WilsonPar(Reader); + //RealD beta = 6.4 ; + WilsonGaugeActionR Waction(WilsonPar.beta); + + ActionLevel Level1(1); + Level1.push_back(&Waction); + //Level1.push_back(WGMod.getPtr()); + TheHMC.TheAction.push_back(Level1); + ///////////////////////////////////////////////////////////// + + // HMC parameters are serialisable + TheHMC.Parameters.initialize(Reader); + + //TheHMC.Parameters.MD.MDsteps = 17; + //TheHMC.Parameters.MD.trajL = 1.0; + + TheHMC.Run(); // no smearing + + Grid_finalize(); + +} // main diff --git a/tests/hmc/Test_hmc_WilsonCloverFermionGauge.cc b/tests/hmc/Test_hmc_WilsonCloverFermionGauge.cc new file mode 100644 index 00000000..8d5701b8 --- /dev/null +++ b/tests/hmc/Test_hmc_WilsonCloverFermionGauge.cc @@ -0,0 +1,129 @@ +/************************************************************************************* + +Grid physics library, www.github.com/paboyle/Grid + +Source file: ./tests/Test_hmc_WilsonFermionGauge.cc + +Copyright (C) 2017 + +Author: Guido Cossu + +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 + +int main(int argc, char **argv) +{ + using namespace Grid; + using namespace Grid::QCD; + + Grid_init(&argc, &argv); + int threads = GridThread::GetThreads(); + // here make a routine to print all the relevant information on the run + std::cout << GridLogMessage << "Grid is setup to use " << threads << " threads" << std::endl; + + // Typedefs to simplify notation + typedef GenericHMCRunner HMCWrapper; // Uses the default minimum norm + typedef WilsonImplR FermionImplPolicy; + typedef WilsonCloverFermionR FermionAction; + typedef typename FermionAction::FermionField FermionField; + + //:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: + HMCWrapper TheHMC; + + // Grid from the command line + TheHMC.Resources.AddFourDimGrid("gauge"); + + // Checkpointer definition + CheckpointerParameters CPparams; + CPparams.config_prefix = "ckpoint_lat"; + CPparams.rng_prefix = "ckpoint_rng"; + CPparams.saveInterval = 5; + CPparams.format = "IEEE64BIG"; + + TheHMC.Resources.LoadNerscCheckpointer(CPparams); + + RNGModuleParameters RNGpar; + RNGpar.serial_seeds = "1 2 3 4 5"; + RNGpar.parallel_seeds = "6 7 8 9 10"; + TheHMC.Resources.SetRNGSeeds(RNGpar); + + // Construct observables + typedef PlaquetteMod PlaqObs; + TheHMC.Resources.AddObservable(); + + typedef PolyakovMod PolyakovObs; + TheHMC.Resources.AddObservable(); + ////////////////////////////////////////////// + + ///////////////////////////////////////////////////////////// + // Collect actions, here use more encapsulation + // need wrappers of the fermionic classes + // that have a complex construction + // standard + RealD beta = 5.6; + WilsonGaugeActionR Waction(beta); + + auto GridPtr = TheHMC.Resources.GetCartesian(); + auto GridRBPtr = TheHMC.Resources.GetRBCartesian(); + + // temporarily need a gauge field + LatticeGaugeField U(GridPtr); + + Real mass = 0.01; + Real csw = 1.0; + + FermionAction FermOp(U, *GridPtr, *GridRBPtr, mass, csw); + ConjugateGradient CG(1.0e-8, 5000); + + TwoFlavourPseudoFermionAction Nf2(FermOp, CG, CG); + + // Set smearing (true/false), default: false + Nf2.is_smeared = false; + + // Collect actions + ActionLevel Level1(1); + Level1.push_back(&Nf2); + + ActionLevel Level2(4); + Level2.push_back(&Waction); + + TheHMC.TheAction.push_back(Level1); + TheHMC.TheAction.push_back(Level2); + ///////////////////////////////////////////////////////////// + + /* + double rho = 0.1; // smearing parameter + int Nsmear = 2; // number of smearing levels + Smear_Stout Stout(rho); + SmearedConfiguration SmearingPolicy( + UGrid, Nsmear, Stout); + */ + + // HMC parameters are serialisable + TheHMC.Parameters.MD.MDsteps = 20; + TheHMC.Parameters.MD.trajL = 1.0; + + TheHMC.ReadCommandLine(argc, argv); // these can be parameters from file + TheHMC.Run(); // no smearing + // TheHMC.Run(SmearingPolicy); // for smearing + + Grid_finalize(); + +} // main diff --git a/tests/lanczos/Test_WCMultiRep_lanczos.cc b/tests/lanczos/Test_WCMultiRep_lanczos.cc new file mode 100644 index 00000000..b6d69aee --- /dev/null +++ b/tests/lanczos/Test_WCMultiRep_lanczos.cc @@ -0,0 +1,178 @@ +/************************************************************************************* + +Grid physics library, www.github.com/paboyle/Grid + +Source file: ./tests/Test_dwf_lanczos.cc + +Copyright (C) 2015 + +Author: Peter Boyle + +This program is free software; you can redistribute it and/or modify +it under the terms of the GNU General Public License as published by +the Free Software Foundation; either version 2 of the License, or +(at your option) any later version. + +This program is distributed in the hope that it will be useful, +but WITHOUT ANY WARRANTY; without even the implied warranty of +MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +GNU General Public License for more details. + +You should have received a copy of the GNU General Public License along +with this program; if not, write to the Free Software Foundation, Inc., +51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. + +See the full license in the file "LICENSE" in the top level distribution +directory +*************************************************************************************/ +/* END LEGAL */ +#include + +using namespace std; +using namespace Grid; +using namespace Grid::QCD; + + +//typedef WilsonCloverFermionR FermionOp; +//typedef typename WilsonFermionR::FermionField FermionField; + +typedef WilsonImplR FundImplPolicy; +typedef WilsonCloverFermionR FundFermionAction; +typedef typename FundFermionAction::FermionField FundFermionField; + +typedef WilsonTwoIndexAntiSymmetricImplR ASymmImplPolicy; +typedef WilsonCloverTwoIndexAntiSymmetricFermionR ASymmFermionAction; +typedef typename ASymmFermionAction::FermionField ASymmFermionField; + + +RealD AllZero(RealD x) { return 0.; } + +int main(int argc, char** argv) { + Grid_init(&argc, &argv); + + GridCartesian* UGrid = SpaceTimeGrid::makeFourDimGrid( + GridDefaultLatt(), GridDefaultSimd(Nd, vComplex::Nsimd()), + GridDefaultMpi()); + GridRedBlackCartesian* UrbGrid = + SpaceTimeGrid::makeFourDimRedBlackGrid(UGrid); + GridCartesian* FGrid = UGrid; + GridRedBlackCartesian* FrbGrid = UrbGrid; + printf("UGrid=%p UrbGrid=%p FGrid=%p FrbGrid=%p\n", UGrid, UrbGrid, FGrid, + FrbGrid); + + std::vector seeds4({1, 2, 3, 4}); + std::vector seeds5({5, 6, 7, 8}); + GridParallelRNG RNG5(FGrid); + RNG5.SeedFixedIntegers(seeds5); + GridParallelRNG RNG4(UGrid); + RNG4.SeedFixedIntegers(seeds4); + GridParallelRNG RNG5rb(FrbGrid); + RNG5.SeedFixedIntegers(seeds5); + + GridParallelRNG pRNG(UGrid); + GridSerialRNG sRNG; + + FundamentalRepresentation::LatticeField Umu(UGrid); + + TwoIndexAntiSymmetricRepresentation HiRep(UGrid); + TwoIndexAntiSymmetricRepresentation::LatticeField UmuAS(UGrid); + + + CheckpointerParameters CPparams; + + CPparams.config_prefix = "ckpoint_lat"; + CPparams.rng_prefix = "ckpoint_rng"; + CPparams.format = "IEEE64BIG"; + +//NerscHmcCheckpointer Checkpoint(std::string("ckpoint_lat"), + // std::string("ckpoint_rng"), 1); + +NerscHmcCheckpointer Checkpoint(CPparams); + + int CNFGSTART=1; + int CNFGEND=2; + int CNFGSTEP=1; + + Real Fundmass = -0.1; + Real Fundcsw = 1.0; + Real ASmass = -0.1; + Real AScsw = 1.0; + + std::cout << "Fund: mass and csw" << Fundmass << " and " << Fundcsw << std::endl; + std::cout << "AS : mass and csw" << ASmass << " and " << AScsw << std::endl; + + const int Nstop = 30; + const int Nk = 40; + const int Np = 40; + const int Nm = Nk + Np; + const int MaxIt = 10000; + RealD resid = 1.0e-8; + + for (int cnfg=CNFGSTART;cnfg<=CNFGEND;cnfg+=CNFGSTEP){ + Checkpoint.CheckpointRestore(cnfg,Umu, sRNG, pRNG); + + //SU4::HotConfiguration(RNG4, Umu); // temporary, then read. + + HiRep.update_representation(Umu); + UmuAS = HiRep.U; + + FundFermionAction FundFermOp(Umu,*FGrid,*FrbGrid, Fundmass, Fundcsw, Fundcsw); + MdagMLinearOperator HermOpFund(FundFermOp); /// <----- + + ASymmFermionAction ASFermOp(UmuAS,*FGrid,*FrbGrid, ASmass, AScsw, AScsw); + MdagMLinearOperator HermOpAS(ASFermOp); /// <----- + + std::vector Coeffs{0, -1.}; + Polynomial FundPolyX(Coeffs); + //Chebyshev FundCheb(0.0, 10., 12); + + FunctionHermOp FundPolyXOp(FundPolyX,HermOpFund); + PlainHermOp FundOp (HermOpFund); + + ImplicitlyRestartedLanczos IRL_Fund(FundOp, FundPolyXOp, Nstop, Nk, Nm, + resid, MaxIt); + + Polynomial ASPolyX(Coeffs); + //Chebyshev ASCheb(0.0, 10., 12); + + FunctionHermOp ASPolyXOp(ASPolyX,HermOpAS); + PlainHermOp ASOp (HermOpAS); + + ImplicitlyRestartedLanczos IRL_AS(ASOp, ASPolyXOp, Nstop, Nk, Nm, + resid, MaxIt); + + std::vector Fundeval(Nm); + std::vector ASeval(Nm); + + FundFermionField Fundsrc(FGrid); + ASymmFermionField ASsrc(FGrid); + + gaussian(RNG5, Fundsrc); + gaussian(RNG5, ASsrc); + + std::vector Fundevec(Nm, FGrid); + std::vector ASevec(Nm, FGrid); + + for (int i = 0; i < 1; i++) { + std::cout << i << " / " << Nm << "Fund: grid pointer " << Fundevec[i]._grid + << std::endl; + }; + for (int i = 0; i < 1; i++) { + std::cout << i << " / " << Nm << "AS: grid pointer " << ASevec[i]._grid + << std::endl; + }; + + int FundNconv, ASNconv; + IRL_Fund.calc(Fundeval, Fundevec, Fundsrc, FundNconv); + IRL_AS.calc(ASeval, ASevec, ASsrc, ASNconv); + + for (int i=0;i +#include +#include +#include + +// Mass +double mq = 0.1; + +// Define Wilson Types +typedef Grid::QCD::WilsonImplR::FermionField FermionField; +typedef Grid::QCD::LatticeGaugeField GaugeField; + +enum ChromaAction +{ + Wilson, // Wilson + WilsonClover // CloverFermions +}; + +namespace Chroma +{ + +class ChromaWrapper +{ +public: + typedef multi1d U; + typedef LatticeFermion T4; + + static void ImportGauge(GaugeField &gr, + QDP::multi1d &ch) + { + Grid::QCD::LorentzColourMatrix LCM; + Grid::Complex cc; + QDP::ColorMatrix cm; + QDP::Complex c; + + std::vector x(4); + QDP::multi1d cx(4); + std::vector gd = gr._grid->GlobalDimensions(); + + for (x[0] = 0; x[0] < gd[0]; x[0]++) + { + for (x[1] = 0; x[1] < gd[1]; x[1]++) + { + for (x[2] = 0; x[2] < gd[2]; x[2]++) + { + for (x[3] = 0; x[3] < gd[3]; x[3]++) + { + cx[0] = x[0]; + cx[1] = x[1]; + cx[2] = x[2]; + cx[3] = x[3]; + Grid::peekSite(LCM, gr, x); + + for (int mu = 0; mu < 4; mu++) + { + for (int i = 0; i < 3; i++) + { + for (int j = 0; j < 3; j++) + { + cc = LCM(mu)()(i, j); + c = QDP::cmplx(QDP::Real(real(cc)), QDP::Real(imag(cc))); + QDP::pokeColor(cm, c, i, j); + } + } + QDP::pokeSite(ch[mu], cm, cx); + } + } + } + } + } + } + + static void ExportGauge(GaugeField &gr, + QDP::multi1d &ch) + { + Grid::QCD::LorentzColourMatrix LCM; + Grid::Complex cc; + QDP::ColorMatrix cm; + QDP::Complex c; + + std::vector x(4); + QDP::multi1d cx(4); + std::vector gd = gr._grid->GlobalDimensions(); + + for (x[0] = 0; x[0] < gd[0]; x[0]++) + { + for (x[1] = 0; x[1] < gd[1]; x[1]++) + { + for (x[2] = 0; x[2] < gd[2]; x[2]++) + { + for (x[3] = 0; x[3] < gd[3]; x[3]++) + { + cx[0] = x[0]; + cx[1] = x[1]; + cx[2] = x[2]; + cx[3] = x[3]; + + for (int mu = 0; mu < 4; mu++) + { + for (int i = 0; i < 3; i++) + { + for (int j = 0; j < 3; j++) + { + cm = QDP::peekSite(ch[mu], cx); + c = QDP::peekColor(cm, i, j); + cc = Grid::Complex(toDouble(real(c)), toDouble(imag(c))); + LCM(mu) + ()(i, j) = cc; + } + } + } + Grid::pokeSite(LCM, gr, x); + } + } + } + } + } + + // Specific for Wilson Fermions + static void ImportFermion(Grid::QCD::LatticeFermion &gr, + QDP::LatticeFermion &ch) + { + Grid::QCD::SpinColourVector F; + Grid::Complex c; + + QDP::Fermion cF; + QDP::SpinVector cS; + QDP::Complex cc; + + std::vector x(4); // explicit 4d fermions in Grid + QDP::multi1d cx(4); + std::vector gd = gr._grid->GlobalDimensions(); + + for (x[0] = 0; x[0] < gd[0]; x[0]++) + { + for (x[1] = 0; x[1] < gd[1]; x[1]++) + { + for (x[2] = 0; x[2] < gd[2]; x[2]++) + { + for (x[3] = 0; x[3] < gd[3]; x[3]++) + { + cx[0] = x[0]; + cx[1] = x[1]; + cx[2] = x[2]; + cx[3] = x[3]; + + Grid::peekSite(F, gr, x); + + for (int j = 0; j < 3; j++) + { + for (int sp = 0; sp < 4; sp++) + { + + c = F()(sp)(j); + + cc = QDP::cmplx(QDP::Real(real(c)), QDP::Real(imag(c))); + + QDP::pokeSpin(cS, cc, sp); + } + QDP::pokeColor(cF, cS, j); + } + QDP::pokeSite(ch, cF, cx); + } + } + } + } + } + + // Specific for 4d Wilson fermions + static void ExportFermion(Grid::QCD::LatticeFermion &gr, + QDP::LatticeFermion &ch) + { + Grid::QCD::SpinColourVector F; + Grid::Complex c; + + QDP::Fermion cF; + QDP::SpinVector cS; + QDP::Complex cc; + + std::vector x(4); // 4d fermions + QDP::multi1d cx(4); + std::vector gd = gr._grid->GlobalDimensions(); + + for (x[0] = 0; x[0] < gd[0]; x[0]++) + { + for (x[1] = 0; x[1] < gd[1]; x[1]++) + { + for (x[2] = 0; x[2] < gd[2]; x[2]++) + { + for (x[3] = 0; x[3] < gd[3]; x[3]++) + { + cx[0] = x[0]; + cx[1] = x[1]; + cx[2] = x[2]; + cx[3] = x[3]; + + cF = QDP::peekSite(ch, cx); + for (int sp = 0; sp < 4; sp++) + { + for (int j = 0; j < 3; j++) + { + cS = QDP::peekColor(cF, j); + cc = QDP::peekSpin(cS, sp); + c = Grid::Complex(QDP::toDouble(QDP::real(cc)), + QDP::toDouble(QDP::imag(cc))); + F() + (sp)(j) = c; + } + } + Grid::pokeSite(F, gr, x); + } + } + } + } + } + + static Handle> GetLinOp(U &u, ChromaAction params) + { + QDP::Real _mq(mq); + QDP::multi1d bcs(QDP::Nd); + + // Boundary conditions + bcs[0] = bcs[1] = bcs[2] = bcs[3] = 1; + + if (params == Wilson) + { + + Chroma::WilsonFermActParams p; + p.Mass = _mq; + AnisoParam_t _apar; + _apar.anisoP = true; + _apar.t_dir = 3; // in 4d + _apar.xi_0 = 2.0; + _apar.nu = 1.0; + p.anisoParam = _apar; + + Chroma::Handle> fbc(new Chroma::SimpleFermBC(bcs)); + Chroma::Handle> cfs(new Chroma::CreateSimpleFermState(fbc)); + Chroma::UnprecWilsonFermAct S_f(cfs, p); + Chroma::Handle> ffs(S_f.createState(u)); + return S_f.linOp(ffs); + } + + if (params == WilsonClover) + { + Chroma::CloverFermActParams p; + p.Mass = _mq; + p.clovCoeffR = QDP::Real(1.0); + p.clovCoeffT = QDP::Real(2.0); + p.u0 = QDP::Real(1.0); + AnisoParam_t _apar; + _apar.anisoP = true; + _apar.t_dir = 3; // in 4d + _apar.xi_0 = 2.0; + _apar.nu = 1.0; + p.anisoParam = _apar; + + Chroma::Handle> fbc(new Chroma::SimpleFermBC(bcs)); + Chroma::Handle> cfs(new Chroma::CreateSimpleFermState(fbc)); + Chroma::UnprecCloverFermAct S_f(cfs, p); + Chroma::Handle> ffs(S_f.createState(u)); + return S_f.linOp(ffs); + } + } +}; +} // namespace Chroma + +void calc_chroma(ChromaAction action, GaugeField &lat, FermionField &src, FermionField &res, int dag) +{ + QDP::multi1d u(4); + Chroma::ChromaWrapper::ImportGauge(lat, u); + + QDP::LatticeFermion check; + QDP::LatticeFermion result; + QDP::LatticeFermion psi; + + Chroma::ChromaWrapper::ImportFermion(src, psi); + + for (int mu = 0; mu < 4; mu++) + { + std::cout << "Imported Gauge norm [" << mu << "] " << QDP::norm2(u[mu]) << std::endl; + } + std::cout << "Imported Fermion norm " << QDP::norm2(psi) << std::endl; + + typedef QDP::LatticeFermion T; + typedef QDP::multi1d U; + + auto linop = Chroma::ChromaWrapper::GetLinOp(u, action); + + printf("Calling Chroma Linop\n"); + fflush(stdout); + + if (dag) + (*linop)(check, psi, Chroma::MINUS); + else + (*linop)(check, psi, Chroma::PLUS); + + printf("Called Chroma Linop\n"); + fflush(stdout); + + // std::cout << "Calling Chroma Linop " << std::endl; + // linop->evenEvenLinOp(tmp, psi, isign); + // check[rb[0]] = tmp; + // linop->oddOddLinOp(tmp, psi, isign); + // check[rb[1]] = tmp; + // linop->evenOddLinOp(tmp, psi, isign); + // check[rb[0]] += tmp; + // linop->oddEvenLinOp(tmp, psi, isign); + // check[rb[1]] += tmp; + + Chroma::ChromaWrapper::ExportFermion(res, check); +} + +void make_gauge(GaugeField &Umu, FermionField &src) +{ + using namespace Grid; + using namespace Grid::QCD; + + std::vector seeds4({1, 2, 3, 4}); + + Grid::GridCartesian *UGrid = (Grid::GridCartesian *)Umu._grid; + Grid::GridParallelRNG RNG4(UGrid); + RNG4.SeedFixedIntegers(seeds4); + Grid::QCD::SU3::HotConfiguration(RNG4, Umu); + + // Fermion field + Grid::gaussian(RNG4, src); + /* + Grid::QCD::SpinColourVector F; + Grid::Complex c; + + + + std::vector x(4); // 4d fermions + std::vector gd = src._grid->GlobalDimensions(); + + for (x[0] = 0; x[0] < gd[0]; x[0]++) + { + for (x[1] = 0; x[1] < gd[1]; x[1]++) + { + for (x[2] = 0; x[2] < gd[2]; x[2]++) + { + for (x[3] = 0; x[3] < gd[3]; x[3]++) + { + for (int sp = 0; sp < 4; sp++) + { + for (int j = 0; j < 3; j++) // colours + { + F()(sp)(j) = Grid::Complex(0.0,0.0); + if (((sp == 0)|| (sp==3)) && (j==2)) + { + c = Grid::Complex(1.0, 0.0); + F()(sp)(j) = c; + } + } + } + Grid::pokeSite(F, src, x); + + } + } + } + } + */ +} + +void calc_grid(ChromaAction action, Grid::QCD::LatticeGaugeField &Umu, Grid::QCD::LatticeFermion &src, Grid::QCD::LatticeFermion &res, int dag) +{ + using namespace Grid; + using namespace Grid::QCD; + + Grid::GridCartesian *UGrid = (Grid::GridCartesian *)Umu._grid; + Grid::GridRedBlackCartesian *UrbGrid = Grid::QCD::SpaceTimeGrid::makeFourDimRedBlackGrid(UGrid); + + Grid::RealD _mass = mq; + + if (action == Wilson) + { + WilsonAnisotropyCoefficients anis; + anis.isAnisotropic = true; + anis.t_direction = 3; + anis.xi_0 = 2.0; + anis.nu = 1.0; + WilsonImplParams iParam; + Grid::QCD::WilsonFermionR Wf(Umu, *UGrid, *UrbGrid, _mass, iParam, anis); + + std::cout << Grid::GridLogMessage << " Calling Grid Wilson Fermion multiply " << std::endl; + + if (dag) + Wf.Mdag(src, res); + else + Wf.M(src, res); + return; + } + + if (action == WilsonClover) + { + Grid::RealD _csw_r = 1.0; + Grid::RealD _csw_t = 2.0; + WilsonAnisotropyCoefficients anis; + anis.isAnisotropic = true; + anis.t_direction = 3; + anis.xi_0 = 2.0; + anis.nu = 1.0; + WilsonImplParams CloverImplParam; + Grid::QCD::WilsonCloverFermionR Wf(Umu, *UGrid, *UrbGrid, _mass, _csw_r, _csw_t, anis, CloverImplParam); + Wf.ImportGauge(Umu); + + std::cout << Grid::GridLogMessage << " Calling Grid Wilson Clover Fermion multiply " << std::endl; + + if (dag) + Wf.Mdag(src, res); + else + Wf.M(src, res); + return; + } + + assert(0); +} + +int main(int argc, char **argv) +{ + + /******************************************************** + * Setup QDP + *********************************************************/ + Chroma::initialize(&argc, &argv); + Chroma::WilsonTypeFermActs4DEnv::registerAll(); + + /******************************************************** + * Setup Grid + *********************************************************/ + Grid::Grid_init(&argc, &argv); + Grid::GridCartesian *UGrid = Grid::QCD::SpaceTimeGrid::makeFourDimGrid(Grid::GridDefaultLatt(), + Grid::GridDefaultSimd(Grid::QCD::Nd, Grid::vComplex::Nsimd()), + Grid::GridDefaultMpi()); + + std::vector gd = UGrid->GlobalDimensions(); + QDP::multi1d nrow(QDP::Nd); + for (int mu = 0; mu < 4; mu++) + nrow[mu] = gd[mu]; + + QDP::Layout::setLattSize(nrow); + QDP::Layout::create(); + + GaugeField Ug(UGrid); + FermionField src(UGrid); + FermionField res_chroma(UGrid); + FermionField res_grid(UGrid); + FermionField only_wilson(UGrid); + FermionField difference(UGrid); + + std::vector ActionList({Wilson, WilsonClover}); + std::vector ActionName({"Wilson", "WilsonClover"}); + + { + + for (int i = 0; i < ActionList.size(); i++) + { + std::cout << "*****************************" << std::endl; + std::cout << "Action " << ActionName[i] << std::endl; + std::cout << "*****************************" << std::endl; + make_gauge(Ug, src); // fills the gauge field and the fermion field with random numbers + + for (int dag = 0; dag < 2; dag++) + { + + { + + std::cout << "Dag = " << dag << std::endl; + + calc_chroma(ActionList[i], Ug, src, res_chroma, dag); + + // Remove the normalisation of Chroma Gauge links ???????? + std::cout << "Norm of Chroma " << ActionName[i] << " multiply " << Grid::norm2(res_chroma) << std::endl; + calc_grid(ActionList[i], Ug, src, res_grid, dag); + + std::cout << "Norm of gauge " << Grid::norm2(Ug) << std::endl; + + std::cout << "Norm of Grid " << ActionName[i] << " multiply " << Grid::norm2(res_grid) << std::endl; + + difference = res_chroma - res_grid; + std::cout << "Norm of difference " << Grid::norm2(difference) << std::endl; + } + } + + std::cout << "Finished test " << std::endl; + + Chroma::finalize(); + } + } +}