From 19da647e3c54a38ff65ea123fbb3aa5981dba2b3 Mon Sep 17 00:00:00 2001 From: Christopher Kelly Date: Fri, 9 Sep 2022 12:47:09 -0400 Subject: [PATCH 1/2] Added support for non-periodic gauge field implementations in the random gauge shift performed at the start of the HMC trajectory (The above required exposing the gauge implementation to the HMC class through the Integrator class) Made the random shift optional (default on) through a parameter in HMCparameters Modified ConjugateBC::CshiftLink such that it supports any shift in -L < shift < L rather than just +-1 Added a tester for the BC-respecting Cshift Fixed a missing system header include in SSE4 intrinsics wrapper Fixed sumD_cpu for single-prec types performing an incorrect conversion to a single-prec data type at the end, that fails to compile on some systems --- Grid/lattice/Lattice_reduction.h | 5 +- Grid/qcd/hmc/HMC.h | 43 ++-- Grid/qcd/hmc/integrators/Integrator.h | 4 +- .../hmc/integrators/Integrator_algorithm.h | 15 +- Grid/qcd/utils/CovariantCshift.h | 30 ++- Grid/qcd/utils/GaugeFix.h | 4 +- Grid/simd/Grid_sse4.h | 2 +- tests/Test_gfield_shift.cc | 183 ++++++++++++++++++ 8 files changed, 248 insertions(+), 38 deletions(-) create mode 100644 tests/Test_gfield_shift.cc diff --git a/Grid/lattice/Lattice_reduction.h b/Grid/lattice/Lattice_reduction.h index bcd09c04..af6f8c00 100644 --- a/Grid/lattice/Lattice_reduction.h +++ b/Grid/lattice/Lattice_reduction.h @@ -91,10 +91,7 @@ inline typename vobj::scalar_objectD sumD_cpu(const vobj *arg, Integer osites) for(int i=0;i * > ObsListType; //pass these from the resource manager @@ -138,26 +142,37 @@ private: GridBase *Grid = U.Grid(); - ////////////////////////////////////////////////////////////////////////////////////////////////////// - // Mainly for DDHMC perform a random translation of U modulo volume - ////////////////////////////////////////////////////////////////////////////////////////////////////// - std::cout << GridLogMessage << "--------------------------------------------------\n"; - std::cout << GridLogMessage << "Random shifting gauge field by ["; - for(int d=0;dNd();d++) { + if(Params.PerformRandomShift){ + ////////////////////////////////////////////////////////////////////////////////////////////////////// + // Mainly for DDHMC perform a random translation of U modulo volume + ////////////////////////////////////////////////////////////////////////////////////////////////////// + std::cout << GridLogMessage << "--------------------------------------------------\n"; + std::cout << GridLogMessage << "Random shifting gauge field by ["; - int L = Grid->GlobalDimensions()[d]; + std::vector Umu(Grid->Nd(), U.Grid()); + for(int mu=0;muNd();mu++) Umu[mu] = PeekIndex(U, mu); - RealD rn_uniform; random(sRNG, rn_uniform); + for(int d=0;dNd();d++) { - int shift = (int) (rn_uniform*L); + int L = Grid->GlobalDimensions()[d]; - std::cout << shift; - if(dNd()-1) std::cout <<","; - else std::cout <<"]\n"; + RealD rn_uniform; random(sRNG, rn_uniform); + + int shift = (int) (rn_uniform*L); + + std::cout << shift; + if(dNd()-1) std::cout <<","; + else std::cout <<"]\n"; - U = Cshift(U,d,shift); + //shift all fields together in a way that respects the gauge BCs + for(int mu=0; mu < Grid->Nd(); mu++) + Umu[mu] = FieldImplementation::CshiftLink(Umu[mu],d,shift); + } + + for(int mu=0;muNd();mu++) PokeIndex(U,Umu[mu],mu); + + std::cout << GridLogMessage << "--------------------------------------------------\n"; } - std::cout << GridLogMessage << "--------------------------------------------------\n"; TheIntegrator.reset_timer(); diff --git a/Grid/qcd/hmc/integrators/Integrator.h b/Grid/qcd/hmc/integrators/Integrator.h index 7f676ce7..f36e1f64 100644 --- a/Grid/qcd/hmc/integrators/Integrator.h +++ b/Grid/qcd/hmc/integrators/Integrator.h @@ -63,10 +63,10 @@ public: }; /*! @brief Class for Molecular Dynamics management */ -template +template class Integrator { protected: - + typedef FieldImplementation_ FieldImplementation; typedef typename FieldImplementation::Field MomentaField; //for readability typedef typename FieldImplementation::Field Field; diff --git a/Grid/qcd/hmc/integrators/Integrator_algorithm.h b/Grid/qcd/hmc/integrators/Integrator_algorithm.h index b05c4ea8..d8b34571 100644 --- a/Grid/qcd/hmc/integrators/Integrator_algorithm.h +++ b/Grid/qcd/hmc/integrators/Integrator_algorithm.h @@ -92,10 +92,11 @@ NAMESPACE_BEGIN(Grid); * P 1/2 P 1/2 */ -template > -class LeapFrog : public Integrator +template > +class LeapFrog : public Integrator { public: + typedef FieldImplementation_ FieldImplementation; typedef LeapFrog Algorithm; INHERIT_FIELD_TYPES(FieldImplementation); @@ -135,13 +136,14 @@ public: } }; -template > -class MinimumNorm2 : public Integrator +template > +class MinimumNorm2 : public Integrator { private: const RealD lambda = 0.1931833275037836; public: + typedef FieldImplementation_ FieldImplementation; INHERIT_FIELD_TYPES(FieldImplementation); MinimumNorm2(GridBase* grid, IntegratorParameters Par, ActionSet& Aset, SmearingPolicy& Sm) @@ -192,8 +194,8 @@ public: } }; -template > -class ForceGradient : public Integrator +template > +class ForceGradient : public Integrator { private: const RealD lambda = 1.0 / 6.0; @@ -202,6 +204,7 @@ private: const RealD theta = 0.0; public: + typedef FieldImplementation_ FieldImplementation; INHERIT_FIELD_TYPES(FieldImplementation); // Looks like dH scales as dt^4. tested wilson/wilson 2 level. diff --git a/Grid/qcd/utils/CovariantCshift.h b/Grid/qcd/utils/CovariantCshift.h index 79cf8e0f..23984145 100644 --- a/Grid/qcd/utils/CovariantCshift.h +++ b/Grid/qcd/utils/CovariantCshift.h @@ -227,26 +227,38 @@ namespace ConjugateBC { //shift = -1 //Out(x) = U_\mu(x-mu) | x_\mu != 0 // = U*_\mu(L-1) | x_\mu == 0 + //shift = 2 + //Out(x) = U_\mu(x+2\hat\mu) | x_\mu < L-2 + // = U*_\mu(1) | x_\mu == L-1 + // = U*_\mu(0) | x_\mu == L-2 + //shift = -2 + //Out(x) = U_\mu(x-2mu) | x_\mu > 1 + // = U*_\mu(L-2) | x_\mu == 0 + // = U*_\mu(L-1) | x_\mu == 1 + //etc template Lattice CshiftLink(const Lattice &Link, int mu, int shift) { GridBase *grid = Link.Grid(); - int Lmu = grid->GlobalDimensions()[mu] - 1; + int Lmu = grid->GlobalDimensions()[mu]; + assert(abs(shift) < Lmu && "Invalid shift value"); Lattice> coor(grid); LatticeCoordinate(coor, mu); Lattice tmp(grid); - if(shift == 1){ - tmp = Cshift(Link, mu, 1); - tmp = where(coor == Lmu, conjugate(tmp), tmp); + if(shift > 0){ + tmp = Cshift(Link, mu, shift); + tmp = where(coor >= Lmu-shift, conjugate(tmp), tmp); return tmp; - }else if(shift == -1){ + }else if(shift < 0){ tmp = Link; - tmp = where(coor == Lmu, conjugate(tmp), tmp); - return Cshift(tmp, mu, -1); - }else assert(0 && "Invalid shift value"); - return tmp; //shuts up the compiler fussing about the return type + tmp = where(coor >= Lmu+shift, conjugate(tmp), tmp); + return Cshift(tmp, mu, shift); + } + + //shift == 0 + return Link; } } diff --git a/Grid/qcd/utils/GaugeFix.h b/Grid/qcd/utils/GaugeFix.h index d9d03c54..fc723fe3 100644 --- a/Grid/qcd/utils/GaugeFix.h +++ b/Grid/qcd/utils/GaugeFix.h @@ -72,12 +72,12 @@ public: //Fix the gauge field Umu //0 < alpha < 1 is related to the step size, cf https://arxiv.org/pdf/1405.5812.pdf - static void SteepestDescentGaugeFix(GaugeLorentz &Umu,Real & alpha,int maxiter,Real Omega_tol, Real Phi_tol,bool Fourier=false,int orthog=-1,bool err_on_no_converge=true) { + static void SteepestDescentGaugeFix(GaugeLorentz &Umu,Real alpha,int maxiter,Real Omega_tol, Real Phi_tol,bool Fourier=false,int orthog=-1,bool err_on_no_converge=true) { GridBase *grid = Umu.Grid(); GaugeMat xform(grid); SteepestDescentGaugeFix(Umu,xform,alpha,maxiter,Omega_tol,Phi_tol,Fourier,orthog,err_on_no_converge); } - static void SteepestDescentGaugeFix(GaugeLorentz &Umu,GaugeMat &xform,Real & alpha,int maxiter,Real Omega_tol, Real Phi_tol,bool Fourier=false,int orthog=-1,bool err_on_no_converge=true) { + static void SteepestDescentGaugeFix(GaugeLorentz &Umu,GaugeMat &xform,Real alpha,int maxiter,Real Omega_tol, Real Phi_tol,bool Fourier=false,int orthog=-1,bool err_on_no_converge=true) { //Fix the gauge field Umu and also return the gauge transformation from the original gauge field, xform GridBase *grid = Umu.Grid(); diff --git a/Grid/simd/Grid_sse4.h b/Grid/simd/Grid_sse4.h index eb76427e..2b2a80fd 100644 --- a/Grid/simd/Grid_sse4.h +++ b/Grid/simd/Grid_sse4.h @@ -35,7 +35,7 @@ Author: neo */ // Time-stamp: <2015-06-16 23:27:54 neo> //---------------------------------------------------------------------- - +#include #include NAMESPACE_BEGIN(Grid); diff --git a/tests/Test_gfield_shift.cc b/tests/Test_gfield_shift.cc new file mode 100644 index 00000000..cf450d3c --- /dev/null +++ b/tests/Test_gfield_shift.cc @@ -0,0 +1,183 @@ + /************************************************************************************* + + Grid physics library, www.github.com/paboyle/Grid + + Source file: ./tests/Test_gfield_shift.cc + + Copyright (C) 2015 + +Author: Christopher Kelly +Author: Azusa Yamaguchi +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 */ + +//Test the shifting of the gauge field that respects the boundary conditions +#include + +using namespace Grid; + ; + +typedef ConjugateGimplR Gimpl; //can choose periodic / charge conjugate directions at wil +typedef Gimpl::GaugeField GaugeField; +typedef Gimpl::GaugeLinkField GaugeLinkField; +typedef Gimpl::SiteGaugeField SiteGaugeField; +typedef Gimpl::SiteGaugeLink SiteGaugeLink; + +GaugeField CshiftGaugeField(const GaugeField &U, const int dir, const int shift){ + GridBase *Grid = U.Grid(); + + GaugeField out(Grid); + GaugeLinkField Umu(Grid); + for(int mu=0;muNd();mu++){ + Umu = PeekIndex(U, mu); + Umu = Gimpl::CshiftLink(Umu,dir,shift); + PokeIndex(out,Umu,mu); + } + return out; +} + + +int main (int argc, char ** argv) +{ + Grid_init(&argc,&argv); + + auto latt_size = GridDefaultLatt(); + auto simd_layout = GridDefaultSimd(4,vComplex::Nsimd()); + auto mpi_layout = GridDefaultMpi(); + + std::vector conj_dirs = {1,1,0,0}; + Gimpl::setDirections(conj_dirs); + + GridCartesian Fine(latt_size,simd_layout,mpi_layout); + + GridParallelRNG FineRNG(&Fine); FineRNG.SeedFixedIntegers(std::vector({45,12,81,9})); + + + GaugeField U(&Fine); + GaugeField ShiftU(&Fine); + + GaugeLinkField link_field(&Fine), link_field_2(&Fine); + + //Like Test_cshift we put the lex coordinate index on each site but make it imaginary + //so we can tell when it was complex conjugated + LatticeComplex lex(&Fine); + lex=Zero(); + U = Zero(); + { + LatticeComplex coor(&Fine); + Integer stride =1; + for(int d=0;d<4;d++){ + LatticeCoordinate(coor,d); + lex = lex + coor*stride; + stride=stride*latt_size[d]; + } + PokeIndex(link_field, lex, 0,0); //place on 0,0 element of link + + for(int mu=0;mu(U, link_field_2, mu); + } + } + + std::stringstream ss; + ss<<"error"; + for(int d=0;d 0 && coor[dir] >= latt_size[dir]-shift && conj_dirs[dir] ) + || + ( shift < 0 && coor[dir] <= -shift-1 && conj_dirs[dir] ) + ) + scm = conjugate(scm); //CC if pulled over boundary + + cm = um(mu)()(0,0); + + RealD nrm = abs(scm-cm()()()); + //std::cout << cm << " " << scm << std::endl; + + Coordinate peer(4); + Complex tmp =cm; + Integer index=real(tmp); + + Integer cm_mu = index / vol4d; + index = index % vol4d; + Lexicographic::CoorFromIndex(peer,index,latt_size); + + if (nrm > 0){ + ferr<<"FAIL mu " << mu << " shift "<< shift<<" in dir "<< dir<<" ["< Date: Mon, 3 Oct 2022 10:59:38 -0400 Subject: [PATCH 2/2] Refactored Wilson flow class; previously the class implemented both iterative and adaptive smearing, but only the iterative method was accessible through the Smearing base class. The implementation of Smearing also forced a clunky need to pass iterative smearing parameters through the constructor but adaptive smearing parameters through the function call. Now there is a WilsonFlowBase class that implements common functionality, and separate WilsonFlow (iterative) and WilsonFlowAdaptive (adaptive) classes, both of which implement Smearing virtual functions. Modified the Wilson flow adaptive smearing step size update to implement the original Ramos definition of the distance, where previously it used the norm of a difference which scales with the volume and so would choose too coarse or too fine steps depending on the volume. This is based on Chulwoo's code. Added a test comparing adaptive (with tuneable tolerance) to iterative Wilson flow smearing on a random gauge configuration. --- Grid/qcd/observables/topological_charge.h | 15 +- Grid/qcd/smearing/WilsonFlow.h | 302 ++++++++++++--------- tests/smearing/Test_WilsonFlow_adaptive.cc | 153 +++++++++++ 3 files changed, 336 insertions(+), 134 deletions(-) create mode 100644 tests/smearing/Test_WilsonFlow_adaptive.cc diff --git a/Grid/qcd/observables/topological_charge.h b/Grid/qcd/observables/topological_charge.h index 7c09a180..220ed738 100644 --- a/Grid/qcd/observables/topological_charge.h +++ b/Grid/qcd/observables/topological_charge.h @@ -31,15 +31,16 @@ directory NAMESPACE_BEGIN(Grid); + struct TopologySmearingParameters : Serializable { GRID_SERIALIZABLE_CLASS_MEMBERS(TopologySmearingParameters, - int, steps, - float, step_size, int, meas_interval, - float, maxTau); + float, init_step_size, + float, maxTau, + float, tolerance); - TopologySmearingParameters(int s = 0, float ss = 0.0f, int mi = 0, float mT = 0.0f): - steps(s), step_size(ss), meas_interval(mi), maxTau(mT){} + TopologySmearingParameters(float ss = 0.0f, int mi = 0, float mT = 0.0f, float tol = 1e-4): + init_step_size(ss), meas_interval(mi), maxTau(mT), tolerance(tol){} template < class ReaderClass > TopologySmearingParameters(Reader& Reader){ @@ -97,8 +98,8 @@ public: if (Pars.do_smearing){ // using wilson flow by default here - WilsonFlow WF(Pars.Smearing.steps, Pars.Smearing.step_size, Pars.Smearing.meas_interval); - WF.smear_adaptive(Usmear, U, Pars.Smearing.maxTau); + WilsonFlowAdaptive WF(Pars.Smearing.init_step_size, Pars.Smearing.maxTau, Pars.Smearing.tolerance, Pars.Smearing.meas_interval); + WF.smear(Usmear, U); Real T0 = WF.energyDensityPlaquette(Pars.Smearing.maxTau, Usmear); std::cout << GridLogMessage << std::setprecision(std::numeric_limits::digits10 + 1) << "T0 : [ " << traj << " ] "<< T0 << std::endl; diff --git a/Grid/qcd/smearing/WilsonFlow.h b/Grid/qcd/smearing/WilsonFlow.h index 0d1ee5d2..f169d02b 100644 --- a/Grid/qcd/smearing/WilsonFlow.h +++ b/Grid/qcd/smearing/WilsonFlow.h @@ -33,27 +33,25 @@ directory NAMESPACE_BEGIN(Grid); template -class WilsonFlow: public Smear{ +class WilsonFlowBase: public Smear{ public: //Store generic measurements to take during smearing process using std::function typedef std::function FunctionType; //int: step, RealD: flow time, GaugeField : the gauge field - -private: - unsigned int Nstep; - RealD epsilon; //for regular smearing this is the time step, for adaptive it is the initial time step - + +protected: std::vector< std::pair > functions; //The int maps to the measurement frequency mutable WilsonGaugeAction SG; - - //Evolve the gauge field by 1 step and update tau - void evolve_step(typename Gimpl::GaugeField &U, RealD &tau) const; - //Evolve the gauge field by 1 step and update tau and the current time step eps - void evolve_step_adaptive(typename Gimpl::GaugeField&U, RealD &tau, RealD &eps, RealD maxTau) const; - + public: INHERIT_GIMPL_TYPES(Gimpl) + explicit WilsonFlowBase(unsigned int meas_interval =1): + SG(WilsonGaugeAction(3.0)) { + // WilsonGaugeAction with beta 3.0 + setDefaultMeasurements(meas_interval); + } + void resetActions(){ functions.clear(); } void addMeasurement(int meas_interval, FunctionType meas){ functions.push_back({meas_interval, meas}); } @@ -64,34 +62,11 @@ public: //and output to stdout void setDefaultMeasurements(int topq_meas_interval = 1); - explicit WilsonFlow(unsigned int Nstep, RealD epsilon, unsigned int interval = 1): - Nstep(Nstep), - epsilon(epsilon), - SG(WilsonGaugeAction(3.0)) { - // WilsonGaugeAction with beta 3.0 - assert(epsilon > 0.0); - LogMessage(); - setDefaultMeasurements(interval); - } - - void LogMessage() { - std::cout << GridLogMessage - << "[WilsonFlow] Nstep : " << Nstep << std::endl; - std::cout << GridLogMessage - << "[WilsonFlow] epsilon : " << epsilon << std::endl; - std::cout << GridLogMessage - << "[WilsonFlow] full trajectory : " << Nstep * epsilon << std::endl; - } - - virtual void smear(GaugeField&, const GaugeField&) const; - - virtual void derivative(GaugeField&, const GaugeField&, const GaugeField&) const { + void derivative(GaugeField&, const GaugeField&, const GaugeField&) const override{ assert(0); // undefined for WilsonFlow } - void smear_adaptive(GaugeField&, const GaugeField&, RealD maxTau) const; - //Compute t^2 for time t from the plaquette static RealD energyDensityPlaquette(const RealD t, const GaugeField& U); @@ -115,82 +90,63 @@ public: std::vector flowMeasureEnergyDensityCloverleaf(const GaugeField& U, int measure_interval = 1); }; +//Basic iterative Wilson flow +template +class WilsonFlow: public WilsonFlowBase{ +private: + int Nstep; //number of steps + RealD epsilon; //step size + + //Evolve the gauge field by 1 step of size eps and update tau + void evolve_step(typename Gimpl::GaugeField &U, RealD &tau) const; + +public: + INHERIT_GIMPL_TYPES(Gimpl) + + //Integrate the Wilson flow for Nstep steps of size epsilon + WilsonFlow(const RealD epsilon, const int Nstep, unsigned int meas_interval = 1): WilsonFlowBase(meas_interval), Nstep(Nstep), epsilon(epsilon){} + + void smear(GaugeField& out, const GaugeField& in) const override; +}; + +//Wilson flow with adaptive step size +template +class WilsonFlowAdaptive: public WilsonFlowBase{ +private: + RealD init_epsilon; //initial step size + RealD maxTau; //integrate to t=maxTau + RealD tolerance; //integration error tolerance + + //Evolve the gauge field by 1 step and update tau and the current time step eps + // + //If the step size eps is too large that a significant integration error results, + //the gauge field (U) and tau will not be updated and the function will return 0; eps will be adjusted to a smaller + //value for the next iteration. + // + //For a successful integration step the function will return 1 + int evolve_step_adaptive(typename Gimpl::GaugeField&U, RealD &tau, RealD &eps) const; + +public: + INHERIT_GIMPL_TYPES(Gimpl) + + WilsonFlowAdaptive(const RealD init_epsilon, const RealD maxTau, const RealD tolerance, unsigned int meas_interval = 1): + WilsonFlowBase(meas_interval), init_epsilon(init_epsilon), maxTau(maxTau), tolerance(tolerance){} + + void smear(GaugeField& out, const GaugeField& in) const override; +}; //////////////////////////////////////////////////////////////////////////////// // Implementations //////////////////////////////////////////////////////////////////////////////// template -void WilsonFlow::evolve_step(typename Gimpl::GaugeField &U, RealD &tau) const{ - GaugeField Z(U.Grid()); - GaugeField tmp(U.Grid()); - SG.deriv(U, Z); - Z *= 0.25; // Z0 = 1/4 * F(U) - Gimpl::update_field(Z, U, -2.0*epsilon); // U = W1 = exp(ep*Z0)*W0 - - Z *= -17.0/8.0; - SG.deriv(U, tmp); Z += tmp; // -17/32*Z0 +Z1 - Z *= 8.0/9.0; // Z = -17/36*Z0 +8/9*Z1 - Gimpl::update_field(Z, U, -2.0*epsilon); // U_= W2 = exp(ep*Z)*W1 - - Z *= -4.0/3.0; - SG.deriv(U, tmp); Z += tmp; // 4/3*(17/36*Z0 -8/9*Z1) +Z2 - Z *= 3.0/4.0; // Z = 17/36*Z0 -8/9*Z1 +3/4*Z2 - Gimpl::update_field(Z, U, -2.0*epsilon); // V(t+e) = exp(ep*Z)*W2 - tau += epsilon; -} - -template -void WilsonFlow::evolve_step_adaptive(typename Gimpl::GaugeField &U, RealD &tau, RealD &eps, RealD maxTau) const{ - if (maxTau - tau < eps){ - eps = maxTau-tau; - } - //std::cout << GridLogMessage << "Integration epsilon : " << epsilon << std::endl; - GaugeField Z(U.Grid()); - GaugeField Zprime(U.Grid()); - GaugeField tmp(U.Grid()), Uprime(U.Grid()); - Uprime = U; - SG.deriv(U, Z); - Zprime = -Z; - Z *= 0.25; // Z0 = 1/4 * F(U) - Gimpl::update_field(Z, U, -2.0*eps); // U = W1 = exp(ep*Z0)*W0 - - Z *= -17.0/8.0; - SG.deriv(U, tmp); Z += tmp; // -17/32*Z0 +Z1 - Zprime += 2.0*tmp; - Z *= 8.0/9.0; // Z = -17/36*Z0 +8/9*Z1 - Gimpl::update_field(Z, U, -2.0*eps); // U_= W2 = exp(ep*Z)*W1 - - - Z *= -4.0/3.0; - SG.deriv(U, tmp); Z += tmp; // 4/3*(17/36*Z0 -8/9*Z1) +Z2 - Z *= 3.0/4.0; // Z = 17/36*Z0 -8/9*Z1 +3/4*Z2 - Gimpl::update_field(Z, U, -2.0*eps); // V(t+e) = exp(ep*Z)*W2 - - // Ramos - Gimpl::update_field(Zprime, Uprime, -2.0*eps); // V'(t+e) = exp(ep*Z')*W0 - // Compute distance as norm^2 of the difference - GaugeField diffU = U - Uprime; - RealD diff = norm2(diffU); - // adjust integration step - - tau += eps; - //std::cout << GridLogMessage << "Adjusting integration step with distance: " << diff << std::endl; - - eps = eps*0.95*std::pow(1e-4/diff,1./3.); - //std::cout << GridLogMessage << "New epsilon : " << epsilon << std::endl; - -} - - -template -RealD WilsonFlow::energyDensityPlaquette(const RealD t, const GaugeField& U){ +RealD WilsonFlowBase::energyDensityPlaquette(const RealD t, const GaugeField& U){ static WilsonGaugeAction SG(3.0); return 2.0 * t * t * SG.S(U)/U.Grid()->gSites(); } //Compute t^2 for time from the 1x1 cloverleaf form template -RealD WilsonFlow::energyDensityCloverleaf(const RealD t, const GaugeField& U){ +RealD WilsonFlowBase::energyDensityCloverleaf(const RealD t, const GaugeField& U){ typedef typename Gimpl::GaugeLinkField GaugeMat; typedef typename Gimpl::GaugeField GaugeLorentz; @@ -215,7 +171,7 @@ RealD WilsonFlow::energyDensityCloverleaf(const RealD t, const GaugeField template -std::vector WilsonFlow::flowMeasureEnergyDensityPlaquette(GaugeField &V, const GaugeField& U, int measure_interval){ +std::vector WilsonFlowBase::flowMeasureEnergyDensityPlaquette(GaugeField &V, const GaugeField& U, int measure_interval){ std::vector out; resetActions(); addMeasurement(measure_interval, [&out](int step, RealD t, const typename Gimpl::GaugeField &U){ @@ -227,13 +183,13 @@ std::vector WilsonFlow::flowMeasureEnergyDensityPlaquette(GaugeFie } template -std::vector WilsonFlow::flowMeasureEnergyDensityPlaquette(const GaugeField& U, int measure_interval){ +std::vector WilsonFlowBase::flowMeasureEnergyDensityPlaquette(const GaugeField& U, int measure_interval){ GaugeField V(U); return flowMeasureEnergyDensityPlaquette(V,U, measure_interval); } template -std::vector WilsonFlow::flowMeasureEnergyDensityCloverleaf(GaugeField &V, const GaugeField& U, int measure_interval){ +std::vector WilsonFlowBase::flowMeasureEnergyDensityCloverleaf(GaugeField &V, const GaugeField& U, int measure_interval){ std::vector out; resetActions(); addMeasurement(measure_interval, [&out](int step, RealD t, const typename Gimpl::GaugeField &U){ @@ -245,16 +201,52 @@ std::vector WilsonFlow::flowMeasureEnergyDensityCloverleaf(GaugeFi } template -std::vector WilsonFlow::flowMeasureEnergyDensityCloverleaf(const GaugeField& U, int measure_interval){ +std::vector WilsonFlowBase::flowMeasureEnergyDensityCloverleaf(const GaugeField& U, int measure_interval){ GaugeField V(U); return flowMeasureEnergyDensityCloverleaf(V,U, measure_interval); } +template +void WilsonFlowBase::setDefaultMeasurements(int topq_meas_interval){ + addMeasurement(1, [](int step, RealD t, const typename Gimpl::GaugeField &U){ + std::cout << GridLogMessage << "[WilsonFlow] Energy density (plaq) : " << step << " " << t << " " << energyDensityPlaquette(t,U) << std::endl; + }); + addMeasurement(topq_meas_interval, [](int step, RealD t, const typename Gimpl::GaugeField &U){ + std::cout << GridLogMessage << "[WilsonFlow] Top. charge : " << step << " " << WilsonLoops::TopologicalCharge(U) << std::endl; + }); +} -//#define WF_TIMING + +template +void WilsonFlow::evolve_step(typename Gimpl::GaugeField &U, RealD &tau) const{ + GaugeField Z(U.Grid()); + GaugeField tmp(U.Grid()); + this->SG.deriv(U, Z); + Z *= 0.25; // Z0 = 1/4 * F(U) + Gimpl::update_field(Z, U, -2.0*epsilon); // U = W1 = exp(ep*Z0)*W0 + + Z *= -17.0/8.0; + this->SG.deriv(U, tmp); Z += tmp; // -17/32*Z0 +Z1 + Z *= 8.0/9.0; // Z = -17/36*Z0 +8/9*Z1 + Gimpl::update_field(Z, U, -2.0*epsilon); // U_= W2 = exp(ep*Z)*W1 + + Z *= -4.0/3.0; + this->SG.deriv(U, tmp); Z += tmp; // 4/3*(17/36*Z0 -8/9*Z1) +Z2 + Z *= 3.0/4.0; // Z = 17/36*Z0 -8/9*Z1 +3/4*Z2 + Gimpl::update_field(Z, U, -2.0*epsilon); // V(t+e) = exp(ep*Z)*W2 + tau += epsilon; +} + template void WilsonFlow::smear(GaugeField& out, const GaugeField& in) const{ + std::cout << GridLogMessage + << "[WilsonFlow] Nstep : " << Nstep << std::endl; + std::cout << GridLogMessage + << "[WilsonFlow] epsilon : " << epsilon << std::endl; + std::cout << GridLogMessage + << "[WilsonFlow] full trajectory : " << Nstep * epsilon << std::endl; + out = in; RealD taus = 0.; for (unsigned int step = 1; step <= Nstep; step++) { //step indicates the number of smearing steps applied at the time of measurement @@ -266,37 +258,93 @@ void WilsonFlow::smear(GaugeField& out, const GaugeField& in) const{ std::cout << "Time to evolve " << diff.count() << " s\n"; #endif //Perform measurements - for(auto const &meas : functions) + for(auto const &meas : this->functions) if( step % meas.first == 0 ) meas.second(step,taus,out); } } + + template -void WilsonFlow::smear_adaptive(GaugeField& out, const GaugeField& in, RealD maxTau) const{ - out = in; - RealD taus = 0.; - RealD eps = epsilon; - unsigned int step = 0; - do{ - step++; - //std::cout << GridLogMessage << "Evolution time :"<< taus << std::endl; - evolve_step_adaptive(out, taus, eps, maxTau); - //Perform measurements - for(auto const &meas : functions) - if( step % meas.first == 0 ) meas.second(step,taus,out); - } while (taus < maxTau); +int WilsonFlowAdaptive::evolve_step_adaptive(typename Gimpl::GaugeField &U, RealD &tau, RealD &eps) const{ + if (maxTau - tau < eps){ + eps = maxTau-tau; + } + //std::cout << GridLogMessage << "Integration epsilon : " << epsilon << std::endl; + GaugeField Z(U.Grid()); + GaugeField Zprime(U.Grid()); + GaugeField tmp(U.Grid()), Uprime(U.Grid()), Usave(U.Grid()); + Uprime = U; + Usave = U; + + this->SG.deriv(U, Z); + Zprime = -Z; + Z *= 0.25; // Z0 = 1/4 * F(U) + Gimpl::update_field(Z, U, -2.0*eps); // U = W1 = exp(ep*Z0)*W0 + + Z *= -17.0/8.0; + this->SG.deriv(U, tmp); Z += tmp; // -17/32*Z0 +Z1 + Zprime += 2.0*tmp; + Z *= 8.0/9.0; // Z = -17/36*Z0 +8/9*Z1 + Gimpl::update_field(Z, U, -2.0*eps); // U_= W2 = exp(ep*Z)*W1 + + + Z *= -4.0/3.0; + this->SG.deriv(U, tmp); Z += tmp; // 4/3*(17/36*Z0 -8/9*Z1) +Z2 + Z *= 3.0/4.0; // Z = 17/36*Z0 -8/9*Z1 +3/4*Z2 + Gimpl::update_field(Z, U, -2.0*eps); // V(t+e) = exp(ep*Z)*W2 + + // Ramos arXiv:1301.4388 + Gimpl::update_field(Zprime, Uprime, -2.0*eps); // V'(t+e) = exp(ep*Z')*W0 + + // Compute distance using Ramos' definition + GaugeField diffU = U - Uprime; + RealD max_dist = 0; + + for(int mu=0;mu(diffU, mu); + RealD dist_mu = sqrt( maxLocalNorm2(diffU_mu) ) /Nc/Nc; //maximize over sites + max_dist = std::max(max_dist, dist_mu); //maximize over mu + } + + int ret; + if(max_dist < tolerance) { + tau += eps; + ret = 1; + } else { + U = Usave; + ret = 0; + } + eps = eps*0.95*std::pow(tolerance/max_dist,1./3.); + std::cout << GridLogMessage << "Adaptive smearing : Distance: "<< max_dist <<" Step successful: " << ret << " New epsilon: " << eps << std::endl; + + return ret; } template -void WilsonFlow::setDefaultMeasurements(int topq_meas_interval){ - addMeasurement(1, [](int step, RealD t, const typename Gimpl::GaugeField &U){ - std::cout << GridLogMessage << "[WilsonFlow] Energy density (plaq) : " << step << " " << t << " " << energyDensityPlaquette(t,U) << std::endl; - }); - addMeasurement(topq_meas_interval, [](int step, RealD t, const typename Gimpl::GaugeField &U){ - std::cout << GridLogMessage << "[WilsonFlow] Top. charge : " << step << " " << WilsonLoops::TopologicalCharge(U) << std::endl; - }); +void WilsonFlowAdaptive::smear(GaugeField& out, const GaugeField& in) const{ + std::cout << GridLogMessage + << "[WilsonFlow] initial epsilon : " << init_epsilon << std::endl; + std::cout << GridLogMessage + << "[WilsonFlow] full trajectory : " << maxTau << std::endl; + std::cout << GridLogMessage + << "[WilsonFlow] tolerance : " << tolerance << std::endl; + out = in; + RealD taus = 0.; + RealD eps = init_epsilon; + unsigned int step = 0; + do{ + int step_success = evolve_step_adaptive(out, taus, eps); + step += step_success; //step will not be incremented if the integration step fails + + //Perform measurements + if(step_success) + for(auto const &meas : this->functions) + if( step % meas.first == 0 ) meas.second(step,taus,out); + } while (taus < maxTau); } + NAMESPACE_END(Grid); diff --git a/tests/smearing/Test_WilsonFlow_adaptive.cc b/tests/smearing/Test_WilsonFlow_adaptive.cc new file mode 100644 index 00000000..23123eb9 --- /dev/null +++ b/tests/smearing/Test_WilsonFlow_adaptive.cc @@ -0,0 +1,153 @@ +/************************************************************************************* + +Grid physics library, www.github.com/paboyle/Grid + +Source file: ./tests/hmc/Test_WilsonFlow_adaptive.cc + +Copyright (C) 2017 + +Author: Christopher Kelly + +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; + +//Linearly interpolate between two nearest times +RealD interpolate(const RealD t_int, const std::vector > &data){ + RealD tdiff1=1e32; int t1_idx=-1; + RealD tdiff2=1e32; int t2_idx=-1; + + for(int i=0;i seeds({1, 2, 3, 4, 5}); + GridSerialRNG sRNG; + GridParallelRNG pRNG(&Grid); + pRNG.SeedFixedIntegers(seeds); + + LatticeGaugeField U(&Grid); + SU::HotConfiguration(pRNG, U); + + int Nstep = 300; + RealD epsilon = 0.01; + RealD maxTau = Nstep*epsilon; + RealD tolerance = 1e-4; + + for(int i=1;i> tolerance; + } + } + std::cout << "Adaptive smear tolerance " << tolerance << std::endl; + + //Setup iterative Wilson flow + WilsonFlow wflow(epsilon,Nstep); + wflow.resetActions(); + + std::vector > meas_orig; + + wflow.addMeasurement(1, [&wflow,&meas_orig](int step, RealD t, const LatticeGaugeField &U){ + std::cout << GridLogMessage << "[WilsonFlow] Computing Cloverleaf energy density for step " << step << std::endl; + meas_orig.push_back( {t, wflow.energyDensityCloverleaf(t,U)} ); + }); + + //Setup adaptive Wilson flow + WilsonFlowAdaptive wflow_ad(epsilon,maxTau,tolerance); + wflow_ad.resetActions(); + + std::vector > meas_adaptive; + + wflow_ad.addMeasurement(1, [&wflow_ad,&meas_adaptive](int step, RealD t, const LatticeGaugeField &U){ + std::cout << GridLogMessage << "[WilsonFlow] Computing Cloverleaf energy density for step " << step << std::endl; + meas_adaptive.push_back( {t, wflow_ad.energyDensityCloverleaf(t,U)} ); + }); + + //Run + LatticeGaugeFieldD Vtmp(U.Grid()); + wflow.smear(Vtmp, U); //basic smear + + Vtmp = Zero(); + wflow_ad.smear(Vtmp, U); + + //Output values for plotting + { + std::ofstream out("wflow_t2E_orig.dat"); + out.precision(16); + for(auto const &e: meas_orig){ + out << e.first << " " << e.second << std::endl; + } + } + { + std::ofstream out("wflow_t2E_adaptive.dat"); + out.precision(16); + for(auto const &e: meas_adaptive){ + out << e.first << " " << e.second << std::endl; + } + } + + //Compare at times available with adaptive smearing + for(int i=0;i