portelli/Grid
1
0
Fork 0
mirror of https://github.com/paboyle/Grid.git synced 2024-11-10 07:55:35 +00:00
Code Releases Activity

Merge pull request #412 from giltirn/patch/adaptive-wflow

Browse Source 
Patch/adaptive wflow
This commit is contained in:
Peter Boyle 2022-10-04 17:23:19 -04:00 committed by GitHub
commit 584a3ee45c
No known key found for this signature in database
GPG Key ID: 4AEE18F83AFDEB23
3 changed files with 336 additions and 134 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

15
Grid/qcd/observables/topological_charge.h
View File

@ -31,15 +31,16 @@ directory
NAMESPACE_BEGIN(Grid); NAMESPACE_BEGIN(Grid);
struct TopologySmearingParameters : Serializable { struct TopologySmearingParameters : Serializable {
GRID_SERIALIZABLE_CLASS_MEMBERS(TopologySmearingParameters, GRID_SERIALIZABLE_CLASS_MEMBERS(TopologySmearingParameters,
int, steps,
float, step_size,
int, meas_interval, 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): TopologySmearingParameters(float ss = 0.0f, int mi = 0, float mT = 0.0f, float tol = 1e-4):
steps(s), step_size(ss), meas_interval(mi), maxTau(mT){} init_step_size(ss), meas_interval(mi), maxTau(mT), tolerance(tol){}
template < class ReaderClass > template < class ReaderClass >
TopologySmearingParameters(Reader<ReaderClass>& Reader){ TopologySmearingParameters(Reader<ReaderClass>& Reader){
@ -97,8 +98,8 @@ public:
if (Pars.do_smearing){ if (Pars.do_smearing){
// using wilson flow by default here // using wilson flow by default here
WilsonFlow<PeriodicGimplR> WF(Pars.Smearing.steps, Pars.Smearing.step_size, Pars.Smearing.meas_interval); WilsonFlowAdaptive<PeriodicGimplR> WF(Pars.Smearing.init_step_size, Pars.Smearing.maxTau, Pars.Smearing.tolerance, Pars.Smearing.meas_interval);
WF.smear_adaptive(Usmear, U, Pars.Smearing.maxTau); WF.smear(Usmear, U);
Real T0 = WF.energyDensityPlaquette(Pars.Smearing.maxTau, Usmear); Real T0 = WF.energyDensityPlaquette(Pars.Smearing.maxTau, Usmear);
std::cout << GridLogMessage << std::setprecision(std::numeric_limits<Real>::digits10 + 1) std::cout << GridLogMessage << std::setprecision(std::numeric_limits<Real>::digits10 + 1)
<< "T0 : [ " << traj << " ] "<< T0 << std::endl; << "T0 : [ " << traj << " ] "<< T0 << std::endl;

288
Grid/qcd/smearing/WilsonFlow.h
View File

@ -33,27 +33,25 @@ directory
NAMESPACE_BEGIN(Grid); NAMESPACE_BEGIN(Grid);
template <class Gimpl> template <class Gimpl>
class WilsonFlow: public Smear<Gimpl>{ class WilsonFlowBase: public Smear<Gimpl>{
public: public:
//Store generic measurements to take during smearing process using std::function //Store generic measurements to take during smearing process using std::function
typedef std::function<void(int, RealD, const typename Gimpl::GaugeField &)> FunctionType; //int: step, RealD: flow time, GaugeField : the gauge field typedef std::function<void(int, RealD, const typename Gimpl::GaugeField &)> FunctionType; //int: step, RealD: flow time, GaugeField : the gauge field
private: protected:
unsigned int Nstep;
RealD epsilon; //for regular smearing this is the time step, for adaptive it is the initial time step
std::vector< std::pair<int, FunctionType> > functions; //The int maps to the measurement frequency std::vector< std::pair<int, FunctionType> > functions; //The int maps to the measurement frequency
mutable WilsonGaugeAction<Gimpl> SG; mutable WilsonGaugeAction<Gimpl> 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: public:
INHERIT_GIMPL_TYPES(Gimpl) INHERIT_GIMPL_TYPES(Gimpl)
explicit WilsonFlowBase(unsigned int meas_interval =1):
SG(WilsonGaugeAction<Gimpl>(3.0)) {
// WilsonGaugeAction with beta 3.0
setDefaultMeasurements(meas_interval);
}
void resetActions(){ functions.clear(); } void resetActions(){ functions.clear(); }
void addMeasurement(int meas_interval, FunctionType meas){ functions.push_back({meas_interval, meas}); } void addMeasurement(int meas_interval, FunctionType meas){ functions.push_back({meas_interval, meas}); }
@ -64,34 +62,11 @@ public:
//and output to stdout //and output to stdout
void setDefaultMeasurements(int topq_meas_interval = 1); void setDefaultMeasurements(int topq_meas_interval = 1);
explicit WilsonFlow(unsigned int Nstep, RealD epsilon, unsigned int interval = 1): void derivative(GaugeField&, const GaugeField&, const GaugeField&) const override{
Nstep(Nstep),
epsilon(epsilon),
SG(WilsonGaugeAction<Gimpl>(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 {
assert(0); assert(0);
// undefined for WilsonFlow // undefined for WilsonFlow
} }
void smear_adaptive(GaugeField&, const GaugeField&, RealD maxTau) const;
//Compute t^2 <E(t)> for time t from the plaquette //Compute t^2 <E(t)> for time t from the plaquette
static RealD energyDensityPlaquette(const RealD t, const GaugeField& U); static RealD energyDensityPlaquette(const RealD t, const GaugeField& U);
@ -115,82 +90,63 @@ public:
std::vector<RealD> flowMeasureEnergyDensityCloverleaf(const GaugeField& U, int measure_interval = 1); std::vector<RealD> flowMeasureEnergyDensityCloverleaf(const GaugeField& U, int measure_interval = 1);
}; };
//Basic iterative Wilson flow
template <class Gimpl>
class WilsonFlow: public WilsonFlowBase<Gimpl>{
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<Gimpl>(meas_interval), Nstep(Nstep), epsilon(epsilon){}
void smear(GaugeField& out, const GaugeField& in) const override;
};
//Wilson flow with adaptive step size
template <class Gimpl>
class WilsonFlowAdaptive: public WilsonFlowBase<Gimpl>{
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<Gimpl>(meas_interval), init_epsilon(init_epsilon), maxTau(maxTau), tolerance(tolerance){}
void smear(GaugeField& out, const GaugeField& in) const override;
};
//////////////////////////////////////////////////////////////////////////////// ////////////////////////////////////////////////////////////////////////////////
// Implementations // Implementations
//////////////////////////////////////////////////////////////////////////////// ////////////////////////////////////////////////////////////////////////////////
template <class Gimpl> template <class Gimpl>
void WilsonFlow<Gimpl>::evolve_step(typename Gimpl::GaugeField &U, RealD &tau) const{ RealD WilsonFlowBase<Gimpl>::energyDensityPlaquette(const RealD t, const GaugeField& U){
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 <class Gimpl>
void WilsonFlow<Gimpl>::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 <class Gimpl>
RealD WilsonFlow<Gimpl>::energyDensityPlaquette(const RealD t, const GaugeField& U){
static WilsonGaugeAction<Gimpl> SG(3.0); static WilsonGaugeAction<Gimpl> SG(3.0);
return 2.0 * t * t * SG.S(U)/U.Grid()->gSites(); return 2.0 * t * t * SG.S(U)/U.Grid()->gSites();
} }
//Compute t^2 <E(t)> for time from the 1x1 cloverleaf form //Compute t^2 <E(t)> for time from the 1x1 cloverleaf form
template <class Gimpl> template <class Gimpl>
RealD WilsonFlow<Gimpl>::energyDensityCloverleaf(const RealD t, const GaugeField& U){ RealD WilsonFlowBase<Gimpl>::energyDensityCloverleaf(const RealD t, const GaugeField& U){
typedef typename Gimpl::GaugeLinkField GaugeMat; typedef typename Gimpl::GaugeLinkField GaugeMat;
typedef typename Gimpl::GaugeField GaugeLorentz; typedef typename Gimpl::GaugeField GaugeLorentz;
@ -215,7 +171,7 @@ RealD WilsonFlow<Gimpl>::energyDensityCloverleaf(const RealD t, const GaugeField
template <class Gimpl> template <class Gimpl>
std::vector<RealD> WilsonFlow<Gimpl>::flowMeasureEnergyDensityPlaquette(GaugeField &V, const GaugeField& U, int measure_interval){ std::vector<RealD> WilsonFlowBase<Gimpl>::flowMeasureEnergyDensityPlaquette(GaugeField &V, const GaugeField& U, int measure_interval){
std::vector<RealD> out; std::vector<RealD> out;
resetActions(); resetActions();
addMeasurement(measure_interval, [&out](int step, RealD t, const typename Gimpl::GaugeField &U){ addMeasurement(measure_interval, [&out](int step, RealD t, const typename Gimpl::GaugeField &U){
@ -227,13 +183,13 @@ std::vector<RealD> WilsonFlow<Gimpl>::flowMeasureEnergyDensityPlaquette(GaugeFie
} }
template <class Gimpl> template <class Gimpl>
std::vector<RealD> WilsonFlow<Gimpl>::flowMeasureEnergyDensityPlaquette(const GaugeField& U, int measure_interval){ std::vector<RealD> WilsonFlowBase<Gimpl>::flowMeasureEnergyDensityPlaquette(const GaugeField& U, int measure_interval){
GaugeField V(U); GaugeField V(U);
return flowMeasureEnergyDensityPlaquette(V,U, measure_interval); return flowMeasureEnergyDensityPlaquette(V,U, measure_interval);
} }
template <class Gimpl> template <class Gimpl>
std::vector<RealD> WilsonFlow<Gimpl>::flowMeasureEnergyDensityCloverleaf(GaugeField &V, const GaugeField& U, int measure_interval){ std::vector<RealD> WilsonFlowBase<Gimpl>::flowMeasureEnergyDensityCloverleaf(GaugeField &V, const GaugeField& U, int measure_interval){
std::vector<RealD> out; std::vector<RealD> out;
resetActions(); resetActions();
addMeasurement(measure_interval, [&out](int step, RealD t, const typename Gimpl::GaugeField &U){ addMeasurement(measure_interval, [&out](int step, RealD t, const typename Gimpl::GaugeField &U){
@ -245,16 +201,52 @@ std::vector<RealD> WilsonFlow<Gimpl>::flowMeasureEnergyDensityCloverleaf(GaugeFi
} }
template <class Gimpl> template <class Gimpl>
std::vector<RealD> WilsonFlow<Gimpl>::flowMeasureEnergyDensityCloverleaf(const GaugeField& U, int measure_interval){ std::vector<RealD> WilsonFlowBase<Gimpl>::flowMeasureEnergyDensityCloverleaf(const GaugeField& U, int measure_interval){
GaugeField V(U); GaugeField V(U);
return flowMeasureEnergyDensityCloverleaf(V,U, measure_interval); return flowMeasureEnergyDensityCloverleaf(V,U, measure_interval);
} }
template <class Gimpl>
void WilsonFlowBase<Gimpl>::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<Gimpl>::TopologicalCharge(U) << std::endl;
});
}
//#define WF_TIMING
template <class Gimpl>
void WilsonFlow<Gimpl>::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 <class Gimpl> template <class Gimpl>
void WilsonFlow<Gimpl>::smear(GaugeField& out, const GaugeField& in) const{ void WilsonFlow<Gimpl>::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; out = in;
RealD taus = 0.; RealD taus = 0.;
for (unsigned int step = 1; step <= Nstep; step++) { //step indicates the number of smearing steps applied at the time of measurement for (unsigned int step = 1; step <= Nstep; step++) { //step indicates the number of smearing steps applied at the time of measurement
@ -266,36 +258,92 @@ void WilsonFlow<Gimpl>::smear(GaugeField& out, const GaugeField& in) const{
std::cout << "Time to evolve " << diff.count() << " s\n"; std::cout << "Time to evolve " << diff.count() << " s\n";
#endif #endif
//Perform measurements //Perform measurements
for(auto const &meas : functions) for(auto const &meas : this->functions)
if( step % meas.first == 0 ) meas.second(step,taus,out); if( step % meas.first == 0 ) meas.second(step,taus,out);
} }
} }
template <class Gimpl> template <class Gimpl>
void WilsonFlow<Gimpl>::smear_adaptive(GaugeField& out, const GaugeField& in, RealD maxTau) const{ int WilsonFlowAdaptive<Gimpl>::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<Nd;mu++){
typename Gimpl::GaugeLinkField diffU_mu = PeekIndex<LorentzIndex>(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 <class Gimpl>
void WilsonFlowAdaptive<Gimpl>::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; out = in;
RealD taus = 0.; RealD taus = 0.;
RealD eps = epsilon; RealD eps = init_epsilon;
unsigned int step = 0; unsigned int step = 0;
do{ do{
step++; int step_success = evolve_step_adaptive(out, taus, eps);
//std::cout << GridLogMessage << "Evolution time :"<< taus << std::endl; step += step_success; //step will not be incremented if the integration step fails
evolve_step_adaptive(out, taus, eps, maxTau);
//Perform measurements //Perform measurements
for(auto const &meas : functions) if(step_success)
for(auto const &meas : this->functions)
if( step % meas.first == 0 ) meas.second(step,taus,out); if( step % meas.first == 0 ) meas.second(step,taus,out);
} while (taus < maxTau); } while (taus < maxTau);
} }
template <class Gimpl>
void WilsonFlow<Gimpl>::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<Gimpl>::TopologicalCharge(U) << std::endl;
});
}
NAMESPACE_END(Grid); NAMESPACE_END(Grid);

153
tests/smearing/Test_WilsonFlow_adaptive.cc Normal file
View File

@ -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 <ckelly@bnl.gov>
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 <Grid/Grid.h>
using namespace Grid;
//Linearly interpolate between two nearest times
RealD interpolate(const RealD t_int, const std::vector<std::pair<RealD,RealD> > &data){
RealD tdiff1=1e32; int t1_idx=-1;
RealD tdiff2=1e32; int t2_idx=-1;
for(int i=0;i<data.size();i++){
RealD diff = fabs(data[i].first-t_int);
//std::cout << "targ " << t_int << " cur " << data[i].first << " diff " << diff << " best diff1 " << tdiff1 << " diff2 " << tdiff2 << std::endl;
if(diff < tdiff1){
if(tdiff1 < tdiff2){ //swap out tdiff2
tdiff2 = tdiff1; t2_idx = t1_idx;
}
tdiff1 = diff; t1_idx = i;
}
else if(diff < tdiff2){ tdiff2 = diff; t2_idx = i; }
}
assert(t1_idx != -1 && t2_idx != -1);
RealD t2 = data[t2_idx].first, v2 = data[t2_idx].second;
RealD t1 = data[t1_idx].first, v1 = data[t1_idx].second;
//v = a + bt
//v2-v1 = b(t2-t1)
RealD b = (v2-v1)/(t2-t1);
RealD a = v1 - b*t1;
RealD vout = a + b*t_int;
//std::cout << "Interpolate to " << t_int << " two closest points " << t1 << " " << t2
//<< " with values " << v1 << " "<< v2 << " : got " << vout << std::endl;
return vout;
}
int main(int argc, char **argv) {
Grid_init(&argc, &argv);
GridLogLayout();
auto latt_size = GridDefaultLatt();
auto simd_layout = GridDefaultSimd(Nd, vComplex::Nsimd());
auto mpi_layout = GridDefaultMpi();
GridCartesian Grid(latt_size, simd_layout, mpi_layout);
GridRedBlackCartesian RBGrid(&Grid);
std::vector<int> seeds({1, 2, 3, 4, 5});
GridSerialRNG sRNG;
GridParallelRNG pRNG(&Grid);
pRNG.SeedFixedIntegers(seeds);
LatticeGaugeField U(&Grid);
SU<Nc>::HotConfiguration(pRNG, U);
int Nstep = 300;
RealD epsilon = 0.01;
RealD maxTau = Nstep*epsilon;
RealD tolerance = 1e-4;
for(int i=1;i<argc;i++){
std::string sarg(argv[i]);
if(sarg == "--tolerance"){
std::stringstream ss; ss << argv[i+1]; ss >> tolerance;
}
}
std::cout << "Adaptive smear tolerance " << tolerance << std::endl;
//Setup iterative Wilson flow
WilsonFlow<PeriodicGimplD> wflow(epsilon,Nstep);
wflow.resetActions();
std::vector<std::pair<RealD, RealD> > 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<PeriodicGimplD> wflow_ad(epsilon,maxTau,tolerance);
wflow_ad.resetActions();
std::vector<std::pair<RealD, RealD> > 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<meas_adaptive.size();i++){
RealD t = meas_adaptive[i].first;
RealD v_adaptive = meas_adaptive[i].second;
RealD v_orig = interpolate(t, meas_orig); //should be very precise due to fine timestep
std::cout << t << " orig: " << v_orig << " adaptive: " << v_adaptive << " reldiff: " << (v_adaptive-v_orig)/v_orig << std::endl;
}
std::cout << GridLogMessage << "Done" << std::endl;
Grid_finalize();
}