diff --git a/lib/qcd/smearing/WilsonFlow.h b/lib/qcd/smearing/WilsonFlow.h index 00b7622e..5e9f2d95 100644 --- a/lib/qcd/smearing/WilsonFlow.h +++ b/lib/qcd/smearing/WilsonFlow.h @@ -37,15 +37,15 @@ template class WilsonFlow: public Smear{ unsigned int Nstep; unsigned int measure_interval; - RealD epsilon; + mutable RealD epsilon, taus; mutable WilsonGaugeAction SG; void evolve_step(typename Gimpl::GaugeField&) const; + void evolve_step_adaptive(typename Gimpl::GaugeField&, RealD); RealD tau(unsigned int t)const {return epsilon*(t+1.0); } - public: INHERIT_GIMPL_TYPES(Gimpl) @@ -75,7 +75,9 @@ class WilsonFlow: public Smear{ // undefined for WilsonFlow } + void smear_adaptive(GaugeField&, const GaugeField&, RealD maxTau); RealD energyDensityPlaquette(unsigned int step, const GaugeField& U) const; + RealD energyDensityPlaquette(const GaugeField& U) const; }; @@ -101,20 +103,70 @@ void WilsonFlow::evolve_step(typename Gimpl::GaugeField &U) const{ Gimpl::update_field(Z, U, -2.0*epsilon); // V(t+e) = exp(ep*Z)*W2 } +template +void WilsonFlow::evolve_step_adaptive(typename Gimpl::GaugeField &U, RealD maxTau) { + if (maxTau - taus < epsilon){ + epsilon = maxTau-taus; + } + 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*epsilon); // 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*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 + + // Ramos + Gimpl::update_field(Zprime, Uprime, -2.0*epsilon); // 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 + + taus += epsilon; + std::cout << GridLogMessage << "Adjusting integration step with distance: " << diff << std::endl; + + epsilon = epsilon*0.95*std::pow(1e-4/diff,1./3.); + std::cout << GridLogMessage << "New epsilon : " << epsilon << std::endl; + +} + template RealD WilsonFlow::energyDensityPlaquette(unsigned int step, const GaugeField& U) const { RealD td = tau(step); return 2.0 * td * td * SG.S(U)/U._grid->gSites(); } +template +RealD WilsonFlow::energyDensityPlaquette(const GaugeField& U) const { + return 2.0 * taus * taus * SG.S(U)/U._grid->gSites(); +} + //#define WF_TIMING + + template void WilsonFlow::smear(GaugeField& out, const GaugeField& in) const { out = in; for (unsigned int step = 1; step <= Nstep; step++) { auto start = std::chrono::high_resolution_clock::now(); + std::cout << GridLogMessage << "Evolution time :"<< tau(step) << std::endl; evolve_step(out); auto end = std::chrono::high_resolution_clock::now(); std::chrono::duration diff = end - start; @@ -132,6 +184,30 @@ void WilsonFlow::smear(GaugeField& out, const GaugeField& in) const { } } +template +void WilsonFlow::smear_adaptive(GaugeField& out, const GaugeField& in, RealD maxTau){ + out = in; + taus = epsilon; + unsigned int step = 0; + do{ + step++; + std::cout << GridLogMessage << "Evolution time :"<< taus << std::endl; + evolve_step_adaptive(out, maxTau); + std::cout << GridLogMessage << "[WilsonFlow] Energy density (plaq) : " + << step << " " + << energyDensityPlaquette(out) << std::endl; + if( step % measure_interval == 0){ + std::cout << GridLogMessage << "[WilsonFlow] Top. charge : " + << step << " " + << WilsonLoops::TopologicalCharge(out) << std::endl; + } + } while (taus < maxTau); + + + +} + + } // namespace QCD } // namespace Grid diff --git a/tests/smearing/Test_WilsonFlow.cc b/tests/smearing/Test_WilsonFlow.cc index c337bf29..5db00d5d 100644 --- a/tests/smearing/Test_WilsonFlow.cc +++ b/tests/smearing/Test_WilsonFlow.cc @@ -33,7 +33,8 @@ namespace Grid{ GRID_SERIALIZABLE_CLASS_MEMBERS(WFParameters, int, steps, double, step_size, - int, meas_interval); + int, meas_interval, + double, maxTau); // for the adaptive algorithm template @@ -97,11 +98,13 @@ int main(int argc, char **argv) { WilsonFlow WF(WFPar.steps, WFPar.step_size, WFPar.meas_interval); - WF.smear(Uflow, Umu); + WF.smear_adaptive(Uflow, Umu, WFPar.maxTau); RealD WFlow_plaq = WilsonLoops::avgPlaquette(Uflow); RealD WFlow_TC = WilsonLoops::TopologicalCharge(Uflow); + RealD WFlow_T0 = WF.energyDensityPlaquette(Uflow); std::cout << GridLogMessage << "Plaquette "<< conf << " " << WFlow_plaq << std::endl; + std::cout << GridLogMessage << "T0 "<< conf << " " << WFlow_T0 << std::endl; std::cout << GridLogMessage << "TopologicalCharge "<< conf << " " << WFlow_TC << std::endl; std::cout<< GridLogMessage << " Admissibility check:\n"; @@ -121,3 +124,29 @@ int main(int argc, char **argv) { } Grid_finalize(); } // main + + +/* +Input file example + + +JSON + +{ + "WilsonFlow":{ + "steps": 200, + "step_size": 0.01, + "meas_interval": 50, + "maxTau": 2.0 + }, + "Configurations":{ + "conf_prefix": "ckpoint_lat", + "rng_prefix": "ckpoint_rng", + "StartConfiguration": 3000, + "EndConfiguration": 3000, + "Skip": 5 + } +} + + +*/ \ No newline at end of file