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