portelli/Grid
1
0
Fork 0
mirror of https://github.com/paboyle/Grid.git synced 2025-04-09 21:50:45 +01:00
Code Releases Activity

Namespace

Browse Source
This commit is contained in:
paboyle 2018-01-14 22:08:25 +00:00
parent 1f49f781bf
commit f4c06ed8c0
1 changed files with 117 additions and 119 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

236
lib/qcd/smearing/WilsonFlow.h
View File

@ -25,59 +25,58 @@ with this program; if not, write to the Free Software Foundation, Inc.,
See the full license in the file "LICENSE" in the top level distribution See the full license in the file "LICENSE" in the top level distribution
directory directory
*************************************************************************************/ *************************************************************************************/
/* END LEGAL */ /* END LEGAL */
#ifndef WILSONFLOW_H #ifndef WILSONFLOW_H
#define WILSONFLOW_H #define WILSONFLOW_H
namespace Grid { NAMESPACE_BEGIN(Grid);
namespace QCD {
template <class Gimpl> template <class Gimpl>
class WilsonFlow: public Smear<Gimpl>{ class WilsonFlow: public Smear<Gimpl>{
unsigned int Nstep; unsigned int Nstep;
unsigned int measure_interval; unsigned int measure_interval;
mutable RealD epsilon, taus; mutable RealD epsilon, taus;
mutable WilsonGaugeAction<Gimpl> SG; mutable WilsonGaugeAction<Gimpl> SG;
void evolve_step(typename Gimpl::GaugeField&) const; void evolve_step(typename Gimpl::GaugeField&) const;
void evolve_step_adaptive(typename Gimpl::GaugeField&, RealD); void evolve_step_adaptive(typename Gimpl::GaugeField&, RealD);
RealD tau(unsigned int t)const {return epsilon*(t+1.0); } RealD tau(unsigned int t)const {return epsilon*(t+1.0); }
public: public:
INHERIT_GIMPL_TYPES(Gimpl) 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):
Nstep(Nstep), Nstep(Nstep),
epsilon(epsilon), epsilon(epsilon),
measure_interval(interval), measure_interval(interval),
SG(WilsonGaugeAction<Gimpl>(3.0)) { SG(WilsonGaugeAction<Gimpl>(3.0)) {
// WilsonGaugeAction with beta 3.0 // WilsonGaugeAction with beta 3.0
assert(epsilon > 0.0); assert(epsilon > 0.0);
LogMessage(); LogMessage();
} }
void LogMessage() { void LogMessage() {
std::cout << GridLogMessage std::cout << GridLogMessage
<< "[WilsonFlow] Nstep : " << Nstep << std::endl; << "[WilsonFlow] Nstep : " << Nstep << std::endl;
std::cout << GridLogMessage std::cout << GridLogMessage
<< "[WilsonFlow] epsilon : " << epsilon << std::endl; << "[WilsonFlow] epsilon : " << epsilon << std::endl;
std::cout << GridLogMessage std::cout << GridLogMessage
<< "[WilsonFlow] full trajectory : " << Nstep * epsilon << std::endl; << "[WilsonFlow] full trajectory : " << Nstep * epsilon << std::endl;
} }
virtual void smear(GaugeField&, const GaugeField&) const; virtual void smear(GaugeField&, const GaugeField&) const;
virtual void derivative(GaugeField&, const 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); void smear_adaptive(GaugeField&, const GaugeField&, RealD maxTau);
RealD energyDensityPlaquette(unsigned int step, const GaugeField& U) const; RealD energyDensityPlaquette(unsigned int step, const GaugeField& U) const;
RealD energyDensityPlaquette(const GaugeField& U) const; RealD energyDensityPlaquette(const GaugeField& U) const;
}; };
@ -86,74 +85,74 @@ class WilsonFlow: public Smear<Gimpl>{
//////////////////////////////////////////////////////////////////////////////// ////////////////////////////////////////////////////////////////////////////////
template <class Gimpl> template <class Gimpl>
void WilsonFlow<Gimpl>::evolve_step(typename Gimpl::GaugeField &U) const{ void WilsonFlow<Gimpl>::evolve_step(typename Gimpl::GaugeField &U) const{
GaugeField Z(U._grid); GaugeField Z(U._grid);
GaugeField tmp(U._grid); GaugeField tmp(U._grid);
SG.deriv(U, Z); SG.deriv(U, Z);
Z *= 0.25; // Z0 = 1/4 * F(U) Z *= 0.25; // Z0 = 1/4 * F(U)
Gimpl::update_field(Z, U, -2.0*epsilon); // U = W1 = exp(ep*Z0)*W0 Gimpl::update_field(Z, U, -2.0*epsilon); // U = W1 = exp(ep*Z0)*W0
Z *= -17.0/8.0; Z *= -17.0/8.0;
SG.deriv(U, tmp); Z += tmp; // -17/32*Z0 +Z1 SG.deriv(U, tmp); Z += tmp; // -17/32*Z0 +Z1
Z *= 8.0/9.0; // Z = -17/36*Z0 +8/9*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 Gimpl::update_field(Z, U, -2.0*epsilon); // U_= W2 = exp(ep*Z)*W1
Z *= -4.0/3.0; Z *= -4.0/3.0;
SG.deriv(U, tmp); Z += tmp; // 4/3*(17/36*Z0 -8/9*Z1) +Z2 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 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 Gimpl::update_field(Z, U, -2.0*epsilon); // V(t+e) = exp(ep*Z)*W2
} }
template <class Gimpl> template <class Gimpl>
void WilsonFlow<Gimpl>::evolve_step_adaptive(typename Gimpl::GaugeField &U, RealD maxTau) { void WilsonFlow<Gimpl>::evolve_step_adaptive(typename Gimpl::GaugeField &U, RealD maxTau) {
if (maxTau - taus < epsilon){ if (maxTau - taus < epsilon){
epsilon = maxTau-taus; epsilon = maxTau-taus;
} }
//std::cout << GridLogMessage << "Integration epsilon : " << epsilon << std::endl; //std::cout << GridLogMessage << "Integration epsilon : " << epsilon << std::endl;
GaugeField Z(U._grid); GaugeField Z(U._grid);
GaugeField Zprime(U._grid); GaugeField Zprime(U._grid);
GaugeField tmp(U._grid), Uprime(U._grid); GaugeField tmp(U._grid), Uprime(U._grid);
Uprime = U; Uprime = U;
SG.deriv(U, Z); SG.deriv(U, Z);
Zprime = -Z; Zprime = -Z;
Z *= 0.25; // Z0 = 1/4 * F(U) Z *= 0.25; // Z0 = 1/4 * F(U)
Gimpl::update_field(Z, U, -2.0*epsilon); // U = W1 = exp(ep*Z0)*W0 Gimpl::update_field(Z, U, -2.0*epsilon); // U = W1 = exp(ep*Z0)*W0
Z *= -17.0/8.0; Z *= -17.0/8.0;
SG.deriv(U, tmp); Z += tmp; // -17/32*Z0 +Z1 SG.deriv(U, tmp); Z += tmp; // -17/32*Z0 +Z1
Zprime += 2.0*tmp; Zprime += 2.0*tmp;
Z *= 8.0/9.0; // Z = -17/36*Z0 +8/9*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 Gimpl::update_field(Z, U, -2.0*epsilon); // U_= W2 = exp(ep*Z)*W1
Z *= -4.0/3.0; Z *= -4.0/3.0;
SG.deriv(U, tmp); Z += tmp; // 4/3*(17/36*Z0 -8/9*Z1) +Z2 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 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 Gimpl::update_field(Z, U, -2.0*epsilon); // V(t+e) = exp(ep*Z)*W2
// Ramos // Ramos
Gimpl::update_field(Zprime, Uprime, -2.0*epsilon); // V'(t+e) = exp(ep*Z')*W0 Gimpl::update_field(Zprime, Uprime, -2.0*epsilon); // V'(t+e) = exp(ep*Z')*W0
// Compute distance as norm^2 of the difference // Compute distance as norm^2 of the difference
GaugeField diffU = U - Uprime; GaugeField diffU = U - Uprime;
RealD diff = norm2(diffU); RealD diff = norm2(diffU);
// adjust integration step // adjust integration step
taus += epsilon; taus += epsilon;
//std::cout << GridLogMessage << "Adjusting integration step with distance: " << diff << std::endl; //std::cout << GridLogMessage << "Adjusting integration step with distance: " << diff << std::endl;
epsilon = epsilon*0.95*std::pow(1e-4/diff,1./3.); epsilon = epsilon*0.95*std::pow(1e-4/diff,1./3.);
//std::cout << GridLogMessage << "New epsilon : " << epsilon << std::endl; //std::cout << GridLogMessage << "New epsilon : " << epsilon << std::endl;
} }
template <class Gimpl> template <class Gimpl>
RealD WilsonFlow<Gimpl>::energyDensityPlaquette(unsigned int step, const GaugeField& U) const { RealD WilsonFlow<Gimpl>::energyDensityPlaquette(unsigned int step, const GaugeField& U) const {
RealD td = tau(step); RealD td = tau(step);
return 2.0 * td * td * SG.S(U)/U._grid->gSites(); return 2.0 * td * td * SG.S(U)/U._grid->gSites();
} }
template <class Gimpl> template <class Gimpl>
RealD WilsonFlow<Gimpl>::energyDensityPlaquette(const GaugeField& U) const { RealD WilsonFlow<Gimpl>::energyDensityPlaquette(const GaugeField& U) const {
return 2.0 * taus * taus * SG.S(U)/U._grid->gSites(); return 2.0 * taus * taus * SG.S(U)/U._grid->gSites();
} }
@ -163,51 +162,50 @@ RealD WilsonFlow<Gimpl>::energyDensityPlaquette(const GaugeField& U) const {
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 {
out = in; out = in;
for (unsigned int step = 1; step <= Nstep; step++) { for (unsigned int step = 1; step <= Nstep; step++) {
auto start = std::chrono::high_resolution_clock::now(); auto start = std::chrono::high_resolution_clock::now();
evolve_step(out); evolve_step(out);
auto end = std::chrono::high_resolution_clock::now(); auto end = std::chrono::high_resolution_clock::now();
std::chrono::duration<double> diff = end - start; std::chrono::duration<double> diff = end - start;
#ifdef WF_TIMING #ifdef WF_TIMING
std::cout << "Time to evolve " << diff.count() << " s\n"; std::cout << "Time to evolve " << diff.count() << " s\n";
#endif #endif
std::cout << GridLogMessage << "[WilsonFlow] Energy density (plaq) : " std::cout << GridLogMessage << "[WilsonFlow] Energy density (plaq) : "
<< step << " " << step << " "
<< energyDensityPlaquette(step,out) << std::endl; << energyDensityPlaquette(step,out) << std::endl;
if( step % measure_interval == 0){ if( step % measure_interval == 0){
std::cout << GridLogMessage << "[WilsonFlow] Top. charge : " std::cout << GridLogMessage << "[WilsonFlow] Top. charge : "
<< step << " " << step << " "
<< WilsonLoops<PeriodicGimplR>::TopologicalCharge(out) << std::endl; << WilsonLoops<PeriodicGimplR>::TopologicalCharge(out) << std::endl;
}
} }
}
} }
template <class Gimpl> template <class Gimpl>
void WilsonFlow<Gimpl>::smear_adaptive(GaugeField& out, const GaugeField& in, RealD maxTau){ void WilsonFlow<Gimpl>::smear_adaptive(GaugeField& out, const GaugeField& in, RealD maxTau){
out = in; out = in;
taus = epsilon; taus = epsilon;
unsigned int step = 0; unsigned int step = 0;
do{ do{
step++; step++;
//std::cout << GridLogMessage << "Evolution time :"<< taus << std::endl; //std::cout << GridLogMessage << "Evolution time :"<< taus << std::endl;
evolve_step_adaptive(out, maxTau); evolve_step_adaptive(out, maxTau);
std::cout << GridLogMessage << "[WilsonFlow] Energy density (plaq) : " std::cout << GridLogMessage << "[WilsonFlow] Energy density (plaq) : "
<< step << " " << step << " "
<< energyDensityPlaquette(out) << std::endl; << energyDensityPlaquette(out) << std::endl;
if( step % measure_interval == 0){ if( step % measure_interval == 0){
std::cout << GridLogMessage << "[WilsonFlow] Top. charge : " std::cout << GridLogMessage << "[WilsonFlow] Top. charge : "
<< step << " " << step << " "
<< WilsonLoops<PeriodicGimplR>::TopologicalCharge(out) << std::endl; << WilsonLoops<PeriodicGimplR>::TopologicalCharge(out) << std::endl;
} }
} while (taus < maxTau); } while (taus < maxTau);
} }
} // namespace QCD NAMESPACE_END(Grid);
} // namespace Grid
#endif // WILSONFLOW_H #endif // WILSONFLOW_H