portelli/Grid
1
0
Fork 0
mirror of https://github.com/paboyle/Grid.git synced 2024-11-09 23:45:36 +00:00
Code Releases Activity

Checking in fixed adaptive WilsonFlow

Browse Source
This commit is contained in:
Chulwoo Jung 2021-06-07 14:20:27 -04:00
parent 92def28bd3
commit b284d50863
2 changed files with 64 additions and 20 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

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

@ -6,7 +6,8 @@ Source file: ./lib/qcd/modules/plaquette.h
Copyright (C) 2017
Author: Guido Cossu <guido.cossu@ed.ac.uk>
Author: Guido Cossu <guido.cossu@ed.ac.uk>
Author: Chulwoo Jung <chulwoo@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
@ -35,7 +36,7 @@ template <class Gimpl>
class WilsonFlow: public Smear<Gimpl>{
unsigned int Nstep;
unsigned int measure_interval;
mutable RealD epsilon, taus;
mutable RealD epsilon, taus,tolerance;
mutable WilsonGaugeAction<Gimpl> SG;
@ -47,13 +48,15 @@ class WilsonFlow: public Smear<Gimpl>{
public:
INHERIT_GIMPL_TYPES(Gimpl)
explicit WilsonFlow(unsigned int Nstep, RealD epsilon, unsigned int interval = 1):
explicit WilsonFlow(unsigned int Nstep, RealD epsilon, unsigned int interval = 1, RealD tol = 1e-3):
Nstep(Nstep),
epsilon(epsilon),
tolerance(tol),
measure_interval(interval),
SG(WilsonGaugeAction<Gimpl>(3.0)) {
// WilsonGaugeAction with beta 3.0
assert(epsilon > 0.0);
assert(tolerance > 0.0);
LogMessage();
}
@ -64,6 +67,8 @@ public:
<< "[WilsonFlow] epsilon : " << epsilon << std::endl;
std::cout << GridLogMessage
<< "[WilsonFlow] full trajectory : " << Nstep * epsilon << std::endl;
std::cout << GridLogMessage
<< "[WilsonFlow] tolerance : " << tolerance << std::endl;
}
virtual void smear(GaugeField&, const GaugeField&) const;
@ -106,11 +111,14 @@ void WilsonFlow<Gimpl>::evolve_step_adaptive(typename Gimpl::GaugeField &U, Real
if (maxTau - taus < epsilon){
epsilon = maxTau-taus;
}
//std::cout << GridLogMessage << "Integration epsilon : " << epsilon << std::endl;
std::cout << GridLogMessage << "Integration epsilon : " << epsilon << std::endl;
GaugeField Z(U.Grid());
GaugeField Zprime(U.Grid());
GaugeField tmp(U.Grid()), Uprime(U.Grid());
GaugeField tmp(U.Grid()), Uprime(U.Grid()),Usave(U.Grid());
Uprime = U;
Usave = U;
SG.deriv(U, Z);
Zprime = -Z;
Z *= 0.25; // Z0 = 1/4 * F(U)
@ -128,18 +136,33 @@ void WilsonFlow<Gimpl>::evolve_step_adaptive(typename Gimpl::GaugeField &U, Real
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
// Ramos arXiv:1301.4388
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
// Wrong
// RealD diff = norm2(diffU);
// std::cout << GridLogMessage << "norm2: " << diff << std::endl;
// RealD tol=1e-3;
RealD diff = real(rankInnerMax(diffU,diffU));
diff = sqrt(diff)/18.; // distance defined in Ramos
GridBase *grid = diffU.Grid();
std::cout << GridLogMessage << "max: " << diff << std::endl;
grid->GlobalMax(diff);
std::cout << GridLogMessage << "max: " << diff << std::endl;
if(diff < tolerance) {
taus += epsilon;
//std::cout << GridLogMessage << "Adjusting integration step with distance: " << diff << std::endl;
// std::cout << GridLogMessage << "Adjusting integration step with distance: " << diff << std::endl;
} else {
U = Usave;
}
epsilon = epsilon*0.95*std::pow(1e-4/diff,1./3.);
//std::cout << GridLogMessage << "New epsilon : " << epsilon << std::endl;
epsilon = epsilon*0.95*std::pow(tolerance/diff,1./3.);
std::cout << GridLogMessage << "Distance : "<<diff<<"New epsilon : " << epsilon << std::endl;
}
@ -184,8 +207,11 @@ void WilsonFlow<Gimpl>::smear(GaugeField& out, const GaugeField& in) const {
template <class Gimpl>
void WilsonFlow<Gimpl>::smear_adaptive(GaugeField& out, const GaugeField& in, RealD maxTau){
out = in;
taus = epsilon;
// taus = epsilon;
taus = 0.;
unsigned int step = 0;
double measTau = epsilon*measure_interval;
std::cout << GridLogMessage << "measTau :"<< measTau << std::endl;
do{
step++;
//std::cout << GridLogMessage << "Evolution time :"<< taus << std::endl;
@ -193,10 +219,12 @@ void WilsonFlow<Gimpl>::smear_adaptive(GaugeField& out, const GaugeField& in, Re
std::cout << GridLogMessage << "[WilsonFlow] Energy density (plaq) : "
<< step << " " << taus << " "
<< energyDensityPlaquette(out) << std::endl;
if( step % measure_interval == 0){
// if( step % measure_interval == 0){
if( taus > measTau ) {
std::cout << GridLogMessage << "[WilsonFlow] Top. charge : "
<< step << " "
<< WilsonLoops<PeriodicGimplR>::TopologicalCharge(out) << std::endl;
measTau += epsilon*measure_interval;
}
} while (taus < maxTau);

30
tests/smearing/Test_WilsonFlow.cc
View File

@ -33,6 +33,7 @@ namespace Grid{
GRID_SERIALIZABLE_CLASS_MEMBERS(WFParameters,
int, steps,
double, step_size,
double, tol,
int, meas_interval,
double, maxTau); // for the adaptive algorithm
@ -82,13 +83,27 @@ int main(int argc, char **argv) {
SU<Nc>::HotConfiguration(pRNG, Umu);
typedef Grid::XmlReader Serialiser;
Serialiser Reader("input.xml");
WFParameters WFPar(Reader);
ConfParameters CPar(Reader);
CheckpointerParameters CPPar(CPar.conf_prefix, CPar.rng_prefix);
// Serialiser Reader("input.xml");
// WFParameters WFPar(Reader);
// ConfParameters CPar(Reader);
// WFParameters WFPar;
int steps = 800;
double step_size=0.02;
double tol=1e-4;
int meas_interval=50;
double maxTau = 16;
// ConfParameters CPar;
// CPar. conf_prefix="configurations/ckpoint_lat";
// CPar. rng_prefix="rngs/ckpoint_rng";
// CPar. StartConfiguration=100,
// CPar. EndConfiguration=110,
// CPar. Skip=1;
// CheckpointerParameters CPPar(CPar.conf_prefix, CPar.rng_prefix);
CheckpointerParameters CPPar("configurations/ckpoint_lat","rngs/ckpoint_rng");
BinaryHmcCheckpointer<PeriodicGimplR> CPBin(CPPar);
for (int conf = CPar.StartConfiguration; conf <= CPar.EndConfiguration; conf+= CPar.Skip){
// for (int conf = CPar.StartConfiguration; conf <= CPar.EndConfiguration; conf+= CPar.Skip){
for (int conf = 100; conf <= 110; conf+= 1){
CPBin.CheckpointRestore(conf, Umu, sRNG, pRNG);
@ -96,9 +111,10 @@ int main(int argc, char **argv) {
std::cout << GridLogMessage << "Initial plaquette: "
<< WilsonLoops<PeriodicGimplR>::avgPlaquette(Umu) << std::endl;
WilsonFlow<PeriodicGimplR> WF(WFPar.steps, WFPar.step_size, WFPar.meas_interval);
WilsonFlow<PeriodicGimplR> WF(steps, step_size, meas_interval);
WF.smear_adaptive(Uflow, Umu, WFPar.maxTau);
// WF.smear_adaptive(Uflow, Umu, maxTau);
WF.smear(Uflow, Umu);
RealD WFlow_plaq = WilsonLoops<PeriodicGimplR>::avgPlaquette(Uflow);
RealD WFlow_TC = WilsonLoops<PeriodicGimplR>::TopologicalCharge(Uflow);