From 4b1cf580e076c0a7f96e32d531dbeab0e7e02c52 Mon Sep 17 00:00:00 2001 From: Guido Cossu Date: Tue, 5 Apr 2016 16:19:30 +0900 Subject: [PATCH] Debugging the Smearing routines --- lib/qcd/action/ActionBase.h | 1 + lib/qcd/hmc/HMC.h | 2 +- lib/qcd/hmc/HmcRunner.h | 31 +++++--- lib/qcd/hmc/integrators/Integrator.h | 27 +++++-- .../hmc/integrators/Integrator_algorithm.h | 23 ++++-- lib/qcd/smearing/APEsmearing.h | 11 ++- lib/qcd/smearing/BaseSmearing.h | 2 +- lib/qcd/smearing/GaugeConfiguration.h | 78 +++++++++++-------- lib/qcd/smearing/StoutSmearing.h | 8 +- tests/Test_main.cc | 2 +- 10 files changed, 118 insertions(+), 67 deletions(-) diff --git a/lib/qcd/action/ActionBase.h b/lib/qcd/action/ActionBase.h index 8b91fb79..7c9f6f1c 100644 --- a/lib/qcd/action/ActionBase.h +++ b/lib/qcd/action/ActionBase.h @@ -35,6 +35,7 @@ template class Action { public: + bool is_smeared = false; // Boundary conditions? // Heatbath? virtual void refresh(const GaugeField &U, GridParallelRNG& pRNG) = 0;// refresh pseudofermions virtual RealD S (const GaugeField &U) = 0; // evaluate the action diff --git a/lib/qcd/hmc/HMC.h b/lib/qcd/hmc/HMC.h index 0b87b03e..b471eb3a 100644 --- a/lib/qcd/hmc/HMC.h +++ b/lib/qcd/hmc/HMC.h @@ -160,7 +160,7 @@ namespace Grid{ ///////////////////////////////////////// // Constructor ///////////////////////////////////////// - HybridMonteCarlo(HMCparameters Pms, IntegratorType &_Int, GridSerialRNG &_sRNG, GridParallelRNG &_pRNG, GaugeField &_U ) : + HybridMonteCarlo(HMCparameters Pms, IntegratorType &_Int, GridSerialRNG &_sRNG, GridParallelRNG &_pRNG, GaugeField &_U) : Params(Pms), TheIntegrator(_Int), sRNG(_sRNG), diff --git a/lib/qcd/hmc/HmcRunner.h b/lib/qcd/hmc/HmcRunner.h index ce92b1d2..6dac0d53 100644 --- a/lib/qcd/hmc/HmcRunner.h +++ b/lib/qcd/hmc/HmcRunner.h @@ -81,11 +81,22 @@ public: NumTraj = ivec[0]; } - // Create integrator - typedef MinimumNorm2 IntegratorType;// change here to change the algorithm - IntegratorParameters MDpar(20); - IntegratorType MDynamics(UGrid,MDpar, TheAction); + GridSerialRNG sRNG; + GridParallelRNG pRNG(UGrid); + LatticeGaugeField U(UGrid); // change this to an extended field (smearing class) + + std::vector SerSeed({1,2,3,4,5}); + std::vector ParSeed({6,7,8,9,10}); + + + // Create integrator, including the smearing policy + SmearedConfiguration SmearingPolicy; // simplest empty smearer, construct here more complex smearers + typedef MinimumNorm2 > IntegratorType;// change here to change the algorithm + IntegratorParameters MDpar(20); + IntegratorType MDynamics(UGrid, MDpar, TheAction, SmearingPolicy); + + // Checkpoint strategy NerscHmcCheckpointer Checkpoint(std::string("ckpoint_lat"),std::string("ckpoint_rng"),1); PlaquetteLogger PlaqLog(std::string("plaq")); @@ -94,12 +105,6 @@ public: HMCpar.StartTrajectory = StartTraj; HMCpar.Trajectories = NumTraj; - GridSerialRNG sRNG; - GridParallelRNG pRNG(UGrid); - LatticeGaugeField U(UGrid); // change this to an extended field (smearing class) - - std::vector SerSeed({1,2,3,4,5}); - std::vector ParSeed({6,7,8,9,10}); if ( StartType == HotStart ) { // Hot start @@ -129,7 +134,11 @@ public: Checkpoint.CheckpointRestore(StartTraj, U, sRNG, pRNG); } - HybridMonteCarlo HMC(HMCpar, MDynamics,sRNG,pRNG,U); // pass the extended field + // Attach the gauge field to the smearing Policy and create the fill the smeared set + // notice that the unit configuration is singular in this procedure + SmearingPolicy.set_GaugeField(U); + + HybridMonteCarlo HMC(HMCpar, MDynamics,sRNG,pRNG,U); HMC.AddObservable(&Checkpoint); HMC.AddObservable(&PlaqLog); diff --git a/lib/qcd/hmc/integrators/Integrator.h b/lib/qcd/hmc/integrators/Integrator.h index e6f6c1cd..5a70a1ec 100644 --- a/lib/qcd/hmc/integrators/Integrator.h +++ b/lib/qcd/hmc/integrators/Integrator.h @@ -68,7 +68,7 @@ namespace Grid{ }; /*! @brief Class for Molecular Dynamics management */ - template + template class Integrator { protected: @@ -85,6 +85,8 @@ namespace Grid{ GaugeField P; + SmearingPolicy &Smearer; + // Should match any legal (SU(n)) gauge field // Need to use this template to match Ncol to pass to SU class template void generate_momenta(Lattice< iVector< iScalar< iMatrix >, Nd> > & P,GridParallelRNG& pRNG){ @@ -113,7 +115,9 @@ namespace Grid{ void update_P(GaugeField &Mom,GaugeField&U, int level,double ep){ for(int a=0; aderiv(U,force); + as[level].actions.at(a)->deriv(U,force); // deriv should not include Ta + if (as[level].actions.at(a)->is_smeared) Smearer.smeared_force(force); + force = Ta(force); Mom = Mom - force*ep; } } @@ -137,6 +141,8 @@ namespace Grid{ ProjectOnGroup(Umu); PokeIndex(U, Umu, mu); } + // Update the smeared fields, can be implemented as observer + Smearer.set_GaugeField(U); } virtual void step (GaugeField& U,int level, int first,int last)=0; @@ -145,14 +151,17 @@ namespace Grid{ Integrator(GridBase* grid, IntegratorParameters Par, - ActionSet & Aset): + ActionSet & Aset, + SmearingPolicy &Sm): Params(Par), as(Aset), P(grid), - levels(Aset.size()) + levels(Aset.size()), + Smearer(Sm) { t_P.resize(levels,0.0); t_U=0.0; + // initialization of smearer delegated outside of Integrator }; virtual ~Integrator(){} @@ -163,7 +172,10 @@ namespace Grid{ generate_momenta(P,pRNG); for(int level=0; level< as.size(); ++level){ for(int actionID=0; actionIDrefresh(U, pRNG); + // get gauge field from the SmearingPolicy and + // based on the boolean is_smeared in actionID + GaugeField& Us = Smearer.get_U(as[level].actions.at(actionID)->is_smeared); + as[level].actions.at(actionID)->refresh(Us, pRNG); } } } @@ -186,7 +198,10 @@ namespace Grid{ // Actions for(int level=0; levelS(U); + // get gauge field from the SmearingPolicy and + // based on the boolean is_smeared in actionID + GaugeField& Us = Smearer.get_U(as[level].actions.at(actionID)->is_smeared); + Hterm = as[level].actions.at(actionID)->S(Us); std::cout< 0){ + GaugeField previous_u(U._grid); + std::cout<< GridLogDebug << "[SmearedConfiguration] Filling SmearedSet\n"; + + previous_u = *ThinLinks; + for(int smearLvl = 0; smearLvl < smearingLevels; ++smearLvl){ + StoutSmearing.smear(SmearedSet[smearLvl],previous_u); + previous_u = SmearedSet[smearLvl]; + } + } - } //==================================================================== GaugeField AnalyticSmearedForce(const GaugeField& SigmaKPrime, @@ -86,9 +92,9 @@ namespace Grid { LatticeReal u(grid), w(grid); LatticeComplex f0(grid), f1(grid), f2(grid); - LatticeReal xi0(grid), xi1(grid), fden(grid), tmp(grid); + LatticeReal xi0(grid), xi1(grid), tmp(grid); LatticeReal u2(grid), w2(grid), cosw(grid); - LatticeComplex emiu(grid), e2iu(grid), qt(grid); + LatticeComplex emiu(grid), e2iu(grid), qt(grid), fden(grid); LatticeComplex r01(grid), r11(grid), r21(grid), r02(grid), r12(grid); LatticeComplex r22(grid), tr1(grid), tr2(grid); LatticeComplex b10(grid), b11(grid), b12(grid), b20(grid), b21(grid), b22(grid); @@ -173,66 +179,70 @@ namespace Grid { links configuration */ /*! @brief Standard constructor */ - SmearedConfiguration(GridCartesian * UGrid, - unsigned int Nsmear, - Smear_Stout& Stout): - smearingLevels(Nsmear), + SmearedConfiguration(GridCartesian * UGrid, + unsigned int Nsmear, + Smear_Stout& Stout): + smearingLevels(Nsmear), StoutSmearing(Stout), - ThinLinks(new GaugeField(UGrid)){ + ThinLinks(NULL){ for (unsigned int i=0; i< smearingLevels; ++i) SmearedSet.push_back(*(new GaugeField(UGrid))); } /*! For just thin links */ - SmearedConfiguration(GridCartesian * UGrid): - smearingLevels(0), + SmearedConfiguration(): + smearingLevels(0), StoutSmearing(), - SmearedSet(0), - ThinLinks(new GaugeField(UGrid)){} + SmearedSet(), + ThinLinks(NULL){} + + + // attach the smeared routines to the thin links U and fill the smeared set + void set_GaugeField(GaugeField& U){ fill_smearedSet(U);} - void set_GaugeField(){ fill_smearedSet(); } void smeared_force(GaugeField& SigmaTilde) const{ GaugeField force = SigmaTilde;//actually = U*SigmaTilde, check this for Grid GaugeLinkField tmp_mu(SigmaTilde._grid); for (int mu = 0; mu < Nd; mu++){ tmp_mu = adj(peekLorentz(SmearedSet[smearingLevels-1], mu)) * peekLorentz(force,mu); - pokeLorentz(force, tmp_mu, mu); + pokeLorentz(force, tmp_mu, mu); } for(int ismr = smearingLevels - 1; ismr > 0; --ismr) force = AnalyticSmearedForce(force,get_smeared_conf(ismr-1)); - + force = AnalyticSmearedForce(force,*ThinLinks); - + for (int mu = 0; mu < Nd; mu++){ tmp_mu = peekLorentz(*ThinLinks, mu) * peekLorentz(force, mu); pokeLorentz(SigmaTilde, tmp_mu, mu); } } - - - GaugeField* get_SmearedU() const{ - return const_cast(&(SmearedSet[smearingLevels-1])); + + + GaugeField& get_SmearedU(){ + return SmearedSet[smearingLevels-1]; } - - GaugeField* get_U(bool smeared=false) const { + + GaugeField& get_U(bool smeared=false) { // get the config, thin links by default if (smeared){ if (smearingLevels) return get_SmearedU(); - else return ThinLinks; + else return *ThinLinks; } - else return ThinLinks; + else return *ThinLinks; } }; } - + } + #endif diff --git a/lib/qcd/smearing/StoutSmearing.h b/lib/qcd/smearing/StoutSmearing.h index 3538417a..dc4b29be 100644 --- a/lib/qcd/smearing/StoutSmearing.h +++ b/lib/qcd/smearing/StoutSmearing.h @@ -20,10 +20,14 @@ namespace Grid { public: INHERIT_GIMPL_TYPES(Gimpl) - Smear_Stout(Smear < Gimpl >* base):SmearBase(base){} + Smear_Stout(Smear < Gimpl >* base):SmearBase(base){ + static_assert(Nc==3, "Stout smearing currently implemented only for Nc==3"); + } /*! Default constructor */ - Smear_Stout():SmearBase(new Smear_APE < Gimpl > ()){} + Smear_Stout():SmearBase(new Smear_APE < Gimpl > ()){ + static_assert(Nc==3, "Stout smearing currently implemented only for Nc==3"); + } ~Smear_Stout(){} diff --git a/tests/Test_main.cc b/tests/Test_main.cc index 467c8a5e..9236fae1 100644 --- a/tests/Test_main.cc +++ b/tests/Test_main.cc @@ -608,7 +608,7 @@ int main (int argc, char ** argv) - // Testing Smearing routine compilation + // Testing Smearing routine compilation, separate in a different file GridCartesian Fine(latt_size,simd_layout,mpi_layout); Smear_APE< PeriodicGimplR > APEsmearing; // periodic gauge implemetation Smear_Stout< PeriodicGimplR > StoutSmearing(&APEsmearing);