diff --git a/extras/Hadrons/Modules.hpp b/extras/Hadrons/Modules.hpp index c27254aa..262795e8 100644 --- a/extras/Hadrons/Modules.hpp +++ b/extras/Hadrons/Modules.hpp @@ -1,25 +1,27 @@ -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include #include #include -#include #include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include #include #include -#include -#include -#include -#include +#include +#include +#include +#include +#include +#include +#include +#include diff --git a/extras/Hadrons/Modules/MAction/WilsonClover.hpp b/extras/Hadrons/Modules/MAction/WilsonClover.hpp new file mode 100644 index 00000000..44b1f0b7 --- /dev/null +++ b/extras/Hadrons/Modules/MAction/WilsonClover.hpp @@ -0,0 +1,142 @@ +/************************************************************************************* + +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); +} + +// execution /////////////////////////////////////////////////////////////////// +template +void TWilsonClover::execute() +{ + 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 = *env().template getObject(par().gauge); + auto &grid = *env().getGrid(); + auto &gridRb = *env().getRbGrid(); + std::vector boundary = strToVec(par().boundary); + typename WilsonCloverFermion::ImplParams implParams(boundary); + FMat *fMatPt = new WilsonCloverFermion(U, grid, gridRb, par().mass, + par().csw_r, + par().csw_t, + par().clover_anisotropy, + implParams); + env().setObject(getName(), fMatPt); +} + +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 78bde5a2..358f7025 100644 --- a/extras/Hadrons/Modules/MContraction/Baryon.hpp +++ b/extras/Hadrons/Modules/MContraction/Baryon.hpp @@ -112,7 +112,9 @@ void TBaryon::execute(void) << " quarks '" << par().q1 << "', '" << par().q2 << "', and '" << par().q3 << "'" << std::endl; - CorrWriter writer(par().output); + std::string output_name = par().output + "." + std::to_string(env().getTrajectory()); + + CorrWriter writer(output_name); PropagatorField1 &q1 = *env().template getObject(par().q1); PropagatorField2 &q2 = *env().template getObject(par().q2); PropagatorField3 &q3 = *env().template getObject(par().q2); diff --git a/extras/Hadrons/Modules/MContraction/Meson.hpp b/extras/Hadrons/Modules/MContraction/Meson.hpp index 7810326a..5355bace 100644 --- a/extras/Hadrons/Modules/MContraction/Meson.hpp +++ b/extras/Hadrons/Modules/MContraction/Meson.hpp @@ -165,8 +165,10 @@ void TMeson::execute(void) LOG(Message) << "Computing meson contractions '" << getName() << "' using" << " quarks '" << par().q1 << "' and '" << par().q2 << "'" << std::endl; + + std::string output_name = par().output + "." + std::to_string(env().getTrajectory()); - CorrWriter writer(par().output); + CorrWriter writer(output_name); std::vector buf; std::vector result; Gamma g5(Gamma::Algebra::Gamma5); diff --git a/extras/Hadrons/Modules/MFermion/GaugeProp.hpp b/extras/Hadrons/Modules/MFermion/GaugeProp.hpp index b4f9edcc..4e802710 100644 --- a/extras/Hadrons/Modules/MFermion/GaugeProp.hpp +++ b/extras/Hadrons/Modules/MFermion/GaugeProp.hpp @@ -43,7 +43,6 @@ private: }; MODULE_REGISTER_NS(GaugeProp, TGaugeProp, MFermion); - /****************************************************************************** * TGaugeProp implementation * ******************************************************************************/ @@ -103,7 +102,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; @@ -112,12 +111,12 @@ void TGaugeProp::execute(void) { if (Ls_ == 1) { - PropToFerm(source, fullSrc, s, c); + PropToFerm(source, fullSrc, s, c); } else { source = zero; - PropToFerm(tmp, fullSrc, s, c); + PropToFerm(tmp, fullSrc, s, c); InsertSlice(tmp, source, 0, 0); InsertSlice(tmp, source, Ls_-1, 0); axpby_ssp_pplus(source, 0., source, 1., source, 0, 0); @@ -133,12 +132,12 @@ 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) { @@ -148,7 +147,7 @@ void TGaugeProp::execute(void) axpby_ssp_pminus(sol, 0., sol, 1., sol, 0, 0); axpby_ssp_pplus(sol, 1., sol, 1., sol, 0, Ls_-1); ExtractSlice(tmp, sol, 0, 0); - 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/Load.cc b/extras/Hadrons/Modules/MGauge/Load.cc index 062e7e98..97be9539 100644 --- a/extras/Hadrons/Modules/MGauge/Load.cc +++ b/extras/Hadrons/Modules/MGauge/Load.cc @@ -66,7 +66,8 @@ void TLoad::setup(void) void TLoad::execute(void) { FieldMetaData header; - std::string fileName = par().file + "." + + std::string fileName = par().file + "ckpoint_lat." + std::to_string(env().getTrajectory()); LOG(Message) << "Loading NERSC configuration from file '" << fileName diff --git a/extras/Hadrons/modules.inc b/extras/Hadrons/modules.inc index 669b08ba..089341c1 100644 --- a/extras/Hadrons/modules.inc +++ b/extras/Hadrons/modules.inc @@ -1,38 +1,41 @@ modules_cc =\ + Modules/MScalar/ChargedProp.cc \ + Modules/MScalar/FreeProp.cc \ Modules/MContraction/WeakHamiltonianEye.cc \ Modules/MContraction/WeakHamiltonianNonEye.cc \ Modules/MContraction/WeakNeutral4ptDisc.cc \ - Modules/MGauge/Load.cc \ - Modules/MGauge/Random.cc \ Modules/MGauge/StochEm.cc \ Modules/MGauge/Unit.cc \ - Modules/MScalar/ChargedProp.cc \ - Modules/MScalar/FreeProp.cc + Modules/MGauge/Load.cc \ + Modules/MGauge/FundtoHirep.cc \ + Modules/MGauge/Random.cc modules_hpp =\ - Modules/MAction/DWF.hpp \ - Modules/MAction/Wilson.hpp \ - Modules/MContraction/Baryon.hpp \ - Modules/MContraction/DiscLoop.hpp \ - Modules/MContraction/Gamma3pt.hpp \ - Modules/MContraction/Meson.hpp \ - Modules/MContraction/WeakHamiltonian.hpp \ - Modules/MContraction/WeakHamiltonianEye.hpp \ - Modules/MContraction/WeakHamiltonianNonEye.hpp \ - Modules/MContraction/WeakNeutral4ptDisc.hpp \ - Modules/MFermion/GaugeProp.hpp \ - Modules/MGauge/Load.hpp \ - Modules/MGauge/Random.hpp \ - Modules/MGauge/StochEm.hpp \ - Modules/MGauge/Unit.hpp \ Modules/MLoop/NoiseLoop.hpp \ Modules/MScalar/ChargedProp.hpp \ - Modules/MScalar/FreeProp.hpp \ Modules/MScalar/Scalar.hpp \ + Modules/MScalar/FreeProp.hpp \ + Modules/MSource/Wall.hpp \ + Modules/MSource/SeqGamma.hpp \ + Modules/MSource/Point.hpp \ + Modules/MSource/Z2.hpp \ + Modules/MFermion/GaugeProp.hpp \ + Modules/MContraction/Meson.hpp \ + Modules/MContraction/WeakHamiltonianNonEye.hpp \ + Modules/MContraction/WeakHamiltonianEye.hpp \ + Modules/MContraction/DiscLoop.hpp \ + Modules/MContraction/Baryon.hpp \ + Modules/MContraction/Gamma3pt.hpp \ + Modules/MContraction/WeakNeutral4ptDisc.hpp \ + Modules/MContraction/WeakHamiltonian.hpp \ Modules/MSink/Point.hpp \ Modules/MSolver/RBPrecCG.hpp \ - Modules/MSource/Point.hpp \ - Modules/MSource/SeqGamma.hpp \ - Modules/MSource/Wall.hpp \ - Modules/MSource/Z2.hpp + Modules/MGauge/StochEm.hpp \ + Modules/MGauge/FundtoHirep.hpp \ + Modules/MGauge/Unit.hpp \ + Modules/MGauge/Load.hpp \ + Modules/MGauge/Random.hpp \ + Modules/MAction/WilsonClover.hpp \ + Modules/MAction/DWF.hpp \ + Modules/MAction/Wilson.hpp diff --git a/lib/qcd/QCD.h b/lib/qcd/QCD.h index 9913a071..9c6d54d4 100644 --- a/lib/qcd/QCD.h +++ b/lib/qcd/QCD.h @@ -421,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); } @@ -437,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); } diff --git a/lib/qcd/action/fermion/Fermion.h b/lib/qcd/action/fermion/Fermion.h index bc8397ba..2a008cb7 100644 --- a/lib/qcd/action/fermion/Fermion.h +++ b/lib/qcd/action/fermion/Fermion.h @@ -106,6 +106,10 @@ 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; @@ -116,6 +120,19 @@ 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 89bd9a15..85d6ffea 100644 --- a/lib/qcd/action/fermion/FermionOperatorImpl.h +++ b/lib/qcd/action/fermion/FermionOperatorImpl.h @@ -1004,6 +1004,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 index 3ec90e06..3c082446 100644 --- a/lib/qcd/action/fermion/WilsonCloverFermion.cc +++ b/lib/qcd/action/fermion/WilsonCloverFermion.cc @@ -235,9 +235,9 @@ void WilsonCloverFermion::MeeDeriv(GaugeField &mat, const FermionField &U, assert(0); // not implemented yet } -FermOpTemplateInstantiate(WilsonCloverFermion); // now only for the fundamental representation -//AdjointFermOpTemplateInstantiate(WilsonCloverFermion); -//TwoIndexFermOpTemplateInstantiate(WilsonCloverFermion); +FermOpTemplateInstantiate(WilsonCloverFermion); +AdjointFermOpTemplateInstantiate(WilsonCloverFermion); +TwoIndexFermOpTemplateInstantiate(WilsonCloverFermion); //GparityFermOpTemplateInstantiate(WilsonCloverFermion); } } diff --git a/lib/qcd/action/fermion/WilsonFermion.h b/lib/qcd/action/fermion/WilsonFermion.h index ca5eba8b..0aea4b68 100644 --- a/lib/qcd/action/fermion/WilsonFermion.h +++ b/lib/qcd/action/fermion/WilsonFermion.h @@ -44,7 +44,8 @@ class WilsonFermionStatic { static const int npoint = 8; }; -struct WilsonAnisotropyCoefficients{ + struct WilsonAnisotropyCoefficients: Serializable + { GRID_SERIALIZABLE_CLASS_MEMBERS(WilsonAnisotropyCoefficients, bool, isAnisotropic, int, t_direction, 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/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 86609ffc..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 + + 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_wilsonFund.cc b/tests/hadrons/Test_hadrons_wilsonFund.cc new file mode 100644 index 00000000..aff8a670 --- /dev/null +++ b/tests/hadrons/Test_hadrons_wilsonFund.cc @@ -0,0 +1,160 @@ +/******************************************************************************* + 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 = 1; + globalPar.trajCounter.end = 2; + globalPar.trajCounter.step = 1; + + 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"); + //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_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_WilsonCloverFermionGauge.cc b/tests/hmc/Test_hmc_WilsonCloverFermionGauge.cc index 322bb304..8d5701b8 100644 --- a/tests/hmc/Test_hmc_WilsonCloverFermionGauge.cc +++ b/tests/hmc/Test_hmc_WilsonCloverFermionGauge.cc @@ -67,6 +67,9 @@ int main(int argc, char **argv) // Construct observables typedef PlaquetteMod PlaqObs; TheHMC.Resources.AddObservable(); + + typedef PolyakovMod PolyakovObs; + TheHMC.Resources.AddObservable(); ////////////////////////////////////////////// ///////////////////////////////////////////////////////////// diff --git a/tests/lanczos/Test_WCMultiRep_lanczos.cc b/tests/lanczos/Test_WCMultiRep_lanczos.cc new file mode 100644 index 00000000..98180db1 --- /dev/null +++ b/tests/lanczos/Test_WCMultiRep_lanczos.cc @@ -0,0 +1,170 @@ +/************************************************************************************* + +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); + ImplicitlyRestartedLanczos IRL_Fund(HermOpFund, FundPolyX, Nstop, Nk, Nm, + resid, MaxIt); + + Polynomial ASPolyX(Coeffs); + Chebyshev ASCheb(0.0, 10., 12); + ImplicitlyRestartedLanczos IRL_AS(HermOpAS, ASPolyX, 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