From 9dc345e8e8858cf3eec2131efd8f46fedd8208df Mon Sep 17 00:00:00 2001 From: Guido Cossu Date: Wed, 13 Jul 2016 17:51:18 +0100 Subject: [PATCH] Debugged smearing and adding HMC functions for hirep --- lib/qcd/QCD.h | 7 +-- lib/qcd/action/ActionBase.h | 11 +++-- lib/qcd/hmc/integrators/Integrator.h | 48 +++++++++++++++++-- .../hmc/integrators/Integrator_algorithm.h | 26 ++++++---- lib/qcd/representations/adjoint.h | 22 +++++---- lib/qcd/representations/fundamental.h | 39 +++++++++++++++ lib/qcd/representations/hmc_types.h | 19 +++++++- lib/qcd/utils/SUn.h | 4 ++ lib/qcd/utils/SUnAdjoint.h | 4 +- tests/Test_hmc_EOWilsonFermionGauge.cc | 2 +- tests/Test_hmc_WilsonAdjointFermionGauge.cc | 4 +- tests/Test_lie_generators.cc | 2 +- 12 files changed, 150 insertions(+), 38 deletions(-) create mode 100644 lib/qcd/representations/fundamental.h diff --git a/lib/qcd/QCD.h b/lib/qcd/QCD.h index c475e74d..f154fa6a 100644 --- a/lib/qcd/QCD.h +++ b/lib/qcd/QCD.h @@ -496,11 +496,12 @@ namespace QCD { #include +#include + +#include + #include #include #include -#include - -#include #endif diff --git a/lib/qcd/action/ActionBase.h b/lib/qcd/action/ActionBase.h index 7c9f6f1c..83bb3ff7 100644 --- a/lib/qcd/action/ActionBase.h +++ b/lib/qcd/action/ActionBase.h @@ -66,13 +66,13 @@ template struct ActionLevel{ public: typedef Action* ActPtr; // now force the same colours as the rest of the code - - int multiplier; + + unsigned int multiplier; std::vector actions; - ActionLevel(int mul = 1) : multiplier(mul) { - assert (mul > 0); + ActionLevel(unsigned int mul = 1) : actions(0), multiplier(mul) { + assert (mul >= 1); }; void push_back(ActPtr ptr){ @@ -83,5 +83,6 @@ public: template using ActionSet = std::vector >; -}} +} +} #endif diff --git a/lib/qcd/hmc/integrators/Integrator.h b/lib/qcd/hmc/integrators/Integrator.h index a8e10f55..e9f0c0c7 100644 --- a/lib/qcd/hmc/integrators/Integrator.h +++ b/lib/qcd/hmc/integrators/Integrator.h @@ -64,7 +64,7 @@ struct IntegratorParameters { }; /*! @brief Class for Molecular Dynamics management */ -template +template class Integrator { protected: typedef IntegratorParameters ParameterType; @@ -81,6 +81,8 @@ class Integrator { SmearingPolicy& Smearer; + RepresentationPolicy Representations; + // Should match any legal (SU(n)) gauge field // Need to use this template to match Ncol to pass to SU class template @@ -108,8 +110,28 @@ class Integrator { << " dt " << ep << " : t_P " << t_P[level] << std::endl; } + // to be used by the actionlevel class to iterate + // over the representations + template + void update_P_core(Level repr_level, GaugeField& Mom, GaugeField& U, + double ep) { + typedef typename Level::LatticeField FieldType; + FieldType Ur = repr_level->getRepresentation();// update U is better + for (int a = 0; a < repr_level.size(); ++a) { + FieldType forceR(U._grid); + // Implement smearing only for the fundamental representation now + repr_level.at(a)->deriv(Ur, forceR); + GaugeField force = repr_level.at(a)->RtoFundamentalProject(forceR); + std::cout << GridLogIntegrator + << "Hirep Force average: " << norm2(force) / (U._grid->gSites()) + << std::endl; + Mom -= force * ep; + } + } + // Add the specialized class for the fundamental case + void update_P(GaugeField& Mom, GaugeField& U, int level, double ep) { - // input U actually not used... + // input U actually not used in the fundamental case for (int a = 0; a < as[level].actions.size(); ++a) { GaugeField force(U._grid); GaugeField& Us = Smearer.get_U(as[level].actions.at(a)->is_smeared); @@ -125,6 +147,8 @@ class Integrator { << std::endl; Mom -= force * ep; } + // Add here the other representations + // as[level].apply(update_P_hireps, Args...) } void update_U(GaugeField& U, double ep) { @@ -147,6 +171,8 @@ class Integrator { } // Update the smeared fields, can be implemented as observer Smearer.set_GaugeField(U); + // Update the higher representations fields + //Representations.update(U);// void functions if fundamental representation } virtual void step(GaugeField& U, int level, int first, int last) = 0; @@ -154,7 +180,7 @@ class Integrator { public: Integrator(GridBase* grid, IntegratorParameters Par, ActionSet& Aset, SmearingPolicy& Sm) - : Params(Par), as(Aset), P(grid), levels(Aset.size()), Smearer(Sm) { + : Params(Par), as(Aset), P(grid), levels(Aset.size()), Smearer(Sm), Representations(grid) { t_P.resize(levels, 0.0); t_U = 0.0; // initialization of smearer delegated outside of Integrator @@ -166,8 +192,16 @@ class Integrator { void refresh(GaugeField& U, GridParallelRNG& pRNG) { std::cout << GridLogIntegrator << "Integrator refresh\n"; generate_momenta(P, pRNG); + + // Update the smeared fields, can be implemented as observer + // necessary to keep the fields updated even after a reject + // of the Metropolis + Smearer.set_GaugeField(U); + // Set the (eventual) representations gauge fields + // Representations.update(U); + // The Smearer is attached to a pointer of the gauge field - // automatically gets the updated field + // automatically gets the correct field // whether or not has been accepted in the previous sweep for (int level = 0; level < as.size(); ++level) { for (int actionID = 0; actionID < as[level].actions.size(); ++actionID) { @@ -178,6 +212,9 @@ class Integrator { as[level].actions.at(actionID)->refresh(Us, pRNG); } } + + + } // Calculate action @@ -220,7 +257,8 @@ class Integrator { t_P[level] = 0; } - for (int step = 0; step < Params.MDsteps; ++step) { // MD step + + for (int step = 0; step < Params.MDsteps; ++step) { // MD step int first_step = (step == 0); int last_step = (step == Params.MDsteps - 1); this->step(U, 0, first_step, last_step); diff --git a/lib/qcd/hmc/integrators/Integrator_algorithm.h b/lib/qcd/hmc/integrators/Integrator_algorithm.h index 66390742..9e12c0c3 100644 --- a/lib/qcd/hmc/integrators/Integrator_algorithm.h +++ b/lib/qcd/hmc/integrators/Integrator_algorithm.h @@ -91,17 +91,19 @@ namespace Grid{ * P 1/2 P 1/2 */ - template class LeapFrog : - public Integrator { + template > class LeapFrog : + public Integrator { public: - typedef LeapFrog Algorithm; + typedef LeapFrog Algorithm; LeapFrog(GridBase* grid, IntegratorParameters Par, ActionSet & Aset, SmearingPolicy & Sm): - Integrator(grid,Par,Aset,Sm) {}; + Integrator(grid,Par,Aset,Sm) {}; void step (GaugeField& U, int level,int _first, int _last){ @@ -138,8 +140,10 @@ namespace Grid{ } }; - template class MinimumNorm2 : - public Integrator { + template > class MinimumNorm2 : + public Integrator { private: const RealD lambda = 0.1931833275037836; @@ -149,7 +153,7 @@ namespace Grid{ IntegratorParameters Par, ActionSet & Aset, SmearingPolicy& Sm): - Integrator(grid,Par,Aset,Sm) {}; + Integrator(grid,Par,Aset,Sm) {}; void step (GaugeField& U, int level, int _first,int _last){ @@ -197,8 +201,10 @@ namespace Grid{ }; - template class ForceGradient : - public Integrator { + template > class ForceGradient : + public Integrator { private: const RealD lambda = 1.0/6.0;; const RealD chi = 1.0/72.0; @@ -211,7 +217,7 @@ namespace Grid{ IntegratorParameters Par, ActionSet & Aset, SmearingPolicy &Sm): - Integrator(grid,Par,Aset, Sm) {}; + Integrator(grid,Par,Aset, Sm) {}; void FG_update_P(GaugeField&U, int level,double fg_dt,double ep){ diff --git a/lib/qcd/representations/adjoint.h b/lib/qcd/representations/adjoint.h index 93f8a27c..b211e122 100644 --- a/lib/qcd/representations/adjoint.h +++ b/lib/qcd/representations/adjoint.h @@ -19,11 +19,16 @@ namespace QCD { template class AdjointRep { public: - typename SU_Adjoint::LatticeAdjMatrix U; + // typdef to be used by the Representations class in HMC to get the + // types for the higher representation fields + typedef typename SU_Adjoint::LatticeAdjMatrix LatticeField; const int Dimension = ncolour * ncolour - 1; - explicit AdjointRep(GridBase* grid):U(grid) {} - void update_representation(const LatticeGaugeField& Uin) { + LatticeField U; + + + explicit AdjointRep(GridBase* grid) : U(grid) {} + LatticeField update_representation(const LatticeGaugeField& Uin) { // Uin is in the fundamental representation // get the U in AdjointRep // (U_adj)_B = tr[e^a U e^b U^dag] @@ -34,19 +39,20 @@ class AdjointRep { U = zero; LatticeGaugeField tmp(Uin._grid); - Vector::Matrix > ta(ncolour * ncolour - 1); + Vector::Matrix> ta(ncolour * ncolour - 1); - // FIXME probably not very efficient to get all the generators everytime + // FIXME probably not very efficient to get all the generators + // everytime for (int a = 0; a < Dimension; a++) SU::generator(a, ta[a]); for (int a = 0; a < Dimension; a++) { - tmp = 2.0 * adj(Uin) * ta[a] * Uin; + tmp = 2.0 * adj(Uin) * ta[a] * Uin; for (int b = 0; b < (ncolour * ncolour - 1); b++) { auto Tr = TensorRemove(trace(tmp * ta[b])); - pokeColour(U, Tr, a,b); + pokeColour(U, Tr, a, b); } } - + } }; diff --git a/lib/qcd/representations/fundamental.h b/lib/qcd/representations/fundamental.h new file mode 100644 index 00000000..886dd399 --- /dev/null +++ b/lib/qcd/representations/fundamental.h @@ -0,0 +1,39 @@ +/* + * Policy classes for the HMC + * Author: Guido Cossu +*/ + +#ifndef FUNDAMENTAL_H +#define FUNDAMENTAL_H + + +namespace Grid { +namespace QCD { + +/* +* This is an helper class for the HMC +* Empty since HMC updates already the fundamental representation +*/ + +template +class FundamentalRep { + public: + const int Dimension = ncolour; + + // typdef to be used by the Representations class in HMC to get the + // types for the higher representation fields + typedef typename SU::LatticeMatrix LatticeField; + + explicit FundamentalRep(GridBase* grid) {} //do nothing + void update_representation(const LatticeGaugeField& Uin) {} // do nothing +}; + +typedef FundamentalRep FundamentalRepresentation; + +} +} + + + + +#endif diff --git a/lib/qcd/representations/hmc_types.h b/lib/qcd/representations/hmc_types.h index b4e92c6b..4a281e75 100644 --- a/lib/qcd/representations/hmc_types.h +++ b/lib/qcd/representations/hmc_types.h @@ -4,16 +4,27 @@ #include #include #include +#include namespace Grid { namespace QCD { -// Utility to add support for representations other than the fundamental +// Supported types +//enum {Fundamental, Adjoint} repr_type; + +// Utility to add support to the HMC for representations other than the fundamental template class Representations{ public: typedef std::tuple Representation_type; + + // To access the Reptypes (FundamentalRepresentation, AdjointRepresentation) + template + using repr_type = typename std::tuple_element::type; + // in order to get the typename of the field use + // type repr_type::LatticeField + Representation_type rep; // Multiple types constructor @@ -29,12 +40,16 @@ public: template inline typename std::enable_if < - I::type update(LatticeGaugeField& U) { + I::type update(LatticeGaugeField& U) { std::get(rep).update_representation(U); update(U); } }; + +typedef Representations JustTheFundamental; + + } } diff --git a/lib/qcd/utils/SUn.h b/lib/qcd/utils/SUn.h index 4d07cf5f..c3e8295f 100644 --- a/lib/qcd/utils/SUn.h +++ b/lib/qcd/utils/SUn.h @@ -728,6 +728,10 @@ typedef SU<2> SU2; typedef SU<3> SU3; typedef SU<4> SU4; typedef SU<5> SU5; + + +typedef SU FundamentalMatrices; + } } #endif diff --git a/lib/qcd/utils/SUnAdjoint.h b/lib/qcd/utils/SUnAdjoint.h index 87b64e93..18c21e0b 100644 --- a/lib/qcd/utils/SUnAdjoint.h +++ b/lib/qcd/utils/SUnAdjoint.h @@ -137,12 +137,14 @@ class SU_Adjoint : public SU { - +// Some useful type names typedef SU_Adjoint<2> SU2Adjoint; typedef SU_Adjoint<3> SU3Adjoint; typedef SU_Adjoint<4> SU4Adjoint; typedef SU_Adjoint<5> SU5Adjoint; + +typedef SU_Adjoint AdjointMatrices; } } diff --git a/tests/Test_hmc_EOWilsonFermionGauge.cc b/tests/Test_hmc_EOWilsonFermionGauge.cc index a322fa2d..6065884a 100644 --- a/tests/Test_hmc_EOWilsonFermionGauge.cc +++ b/tests/Test_hmc_EOWilsonFermionGauge.cc @@ -67,7 +67,7 @@ public: TwoFlavourEvenOddPseudoFermionAction Nf2(FermOp,CG,CG); //Set smearing (true/false), default: false - Nf2.is_smeared=false; + Nf2.is_smeared=true; //Collect actions ActionLevel Level1(1); diff --git a/tests/Test_hmc_WilsonAdjointFermionGauge.cc b/tests/Test_hmc_WilsonAdjointFermionGauge.cc index 3f5bf82c..b0a5cd2c 100644 --- a/tests/Test_hmc_WilsonAdjointFermionGauge.cc +++ b/tests/Test_hmc_WilsonAdjointFermionGauge.cc @@ -43,8 +43,8 @@ class HmcRunner : public NerscHmcRunner { void BuildTheAction(int argc, char **argv) { - typedef WilsonImplR ImplPolicy; // gauge field implemetation - typedef WilsonFermionR FermionAction; // type of lattice fermions + typedef WilsonImplR ImplPolicy; // gauge field implemetation for the pseudofermions + typedef WilsonFermionR FermionAction; // type of lattice fermions (Wilson, DW, ...) typedef typename FermionAction::FermionField FermionField; UGrid = SpaceTimeGrid::makeFourDimGrid( diff --git a/tests/Test_lie_generators.cc b/tests/Test_lie_generators.cc index 13f711eb..b30be869 100644 --- a/tests/Test_lie_generators.cc +++ b/tests/Test_lie_generators.cc @@ -112,7 +112,7 @@ int main(int argc, char** argv) { AdjointRep<3> AdjRep(grid); // AdjointRepresentation has the predefined number of colours Nc - Representations RepresentationTypes(grid); + Representations RepresentationTypes(grid); Grid_finalize(); }