portelli/Grid
1
0
Fork 0
mirror of https://github.com/paboyle/Grid.git synced 2024-09-21 01:25:48 +01:00
Code Releases Activity
145a295231
Grid/lib/qcd/utils/WilsonLoops.h

379 lines
10 KiB
C
Raw Normal View History

Adding some wilson loop support
2015-06-05 10:02:36 +01:00
#ifndef QCD_UTILS_WILSON_LOOPS_H
#define QCD_UTILS_WILSON_LOOPS_H
namespace Grid {
namespace QCD {
// Common wilson loop observables
Global edit on HMC sector -- making GaugeField a template parameter and preparing to pass integrator, smearing, bc's as policy classes to hmc. Propose to unify "integrator" and integrator algorithm in a base/derived way to override step. Want to read through ForceGradient to ensure that abstraction covers the force gradient case.
2015-08-30 12:18:34 +01:00
template<class GaugeLorentz>
Adding some wilson loop support
2015-06-05 10:02:36 +01:00
class WilsonLoops {
public:
Global edit on HMC sector -- making GaugeField a template parameter and preparing to pass integrator, smearing, bc's as policy classes to hmc. Propose to unify "integrator" and integrator algorithm in a base/derived way to override step. Want to read through ForceGradient to ensure that abstraction covers the force gradient case.
2015-08-30 12:18:34 +01:00
typedef LorentzScalar<GaugeLorentz> GaugeMat;
Adding some wilson loop support
2015-06-05 10:02:36 +01:00
//////////////////////////////////////////////////
// directed plaquette oriented in mu,nu plane
//////////////////////////////////////////////////
static void dirPlaquette(GaugeMat &plaq,const std::vector<GaugeMat> &U, const int mu, const int nu)
{
plaq=CovShiftForward(U[mu],mu,U[nu])*adj(CovShiftForward(U[nu],nu,U[mu]));
}
//////////////////////////////////////////////////
// trace of directed plaquette oriented in mu,nu plane
//////////////////////////////////////////////////
static void traceDirPlaquette(LatticeComplex &plaq, const std::vector<GaugeMat> &U, const int mu, const int nu)
{
GaugeMat sp(U[0]._grid);
dirPlaquette(sp,U,mu,nu);
plaq=trace(sp);
}
//////////////////////////////////////////////////
// sum over all planes of plaquette
//////////////////////////////////////////////////
static void sitePlaquette(LatticeComplex &Plaq,const std::vector<GaugeMat> &U)
{
LatticeComplex sitePlaq(U[0]._grid);
uninitialised bug fix
2015-06-16 14:07:05 +01:00
Plaq=zero;
Adding some wilson loop support
2015-06-05 10:02:36 +01:00
for(int mu=1;mu<Nd;mu++){
for(int nu=0;nu<mu;nu++){
traceDirPlaquette(sitePlaq,U,mu,nu);
Plaq = Plaq + sitePlaq;
}
}
}
//////////////////////////////////////////////////
// sum over all x,y,z,t and over all planes of plaquette
//////////////////////////////////////////////////
static RealD sumPlaquette(const GaugeLorentz &Umu){
std::vector<GaugeMat> U(4,Umu._grid);
uninitialised bug fix
2015-06-16 14:07:05 +01:00
Adding some wilson loop support
2015-06-05 10:02:36 +01:00
for(int mu=0;mu<Nd;mu++){
big commit fixing nocompiles in defective C++11 compilers (gcc, icpc). stared getting to near the bleeding edge I guess
2015-06-30 15:03:11 +01:00
U[mu] = PeekIndex<LorentzIndex>(Umu,mu);
Adding some wilson loop support
2015-06-05 10:02:36 +01:00
}
LatticeComplex Plaq(Umu._grid);
sitePlaquette(Plaq,U);
TComplex Tp = sum(Plaq);
Complex p = TensorRemove(Tp);
return p.real();
}
//////////////////////////////////////////////////
// average over all x,y,z,t and over all planes of plaquette
//////////////////////////////////////////////////
static RealD avgPlaquette(const GaugeLorentz &Umu){
RealD sumplaq = sumPlaquette(Umu);
double vol = Umu._grid->gSites();
double faces = (1.0*Nd*(Nd-1))/2.0;
return sumplaq/vol/faces/Nc; // Nd , Nc dependent... FIXME
}
Binary IO file for generic Grid array parallel I/O. Number of IO MPI tasks can be varied by selecting which dimensions use parallel IO and which dimensions use Serial send to boss I/O. Thus can neck down from, say 1024 nodes = 4x4x8x8 to {1,8,32,64,128,256,1024} nodes doing the I/O. Interpolates nicely between ALL nodes write their data, a single boss per time-plane in processor space [old UKQCD fortran code did this], and a single node doing all I/O. Not sure I have the transfer sizes big enough and am not overly convinced fstream is guaranteed to not give buffer inconsistencies unless I set streambuf size to zero. Practically it has worked on 8 tasks, 2x1x2x2 writing /cloning NERSC configurations on my MacOS + OpenMPI and Clang environment. It is VERY easy to switch to pwrite at a later date, and also easy to send x-strips around from each node in order to gather bigger chunks at the syscall level. That would push us up to the circa 8x 18*4*8 == 4KB size write chunk, and by taking, say, x/y non parallel we get to 16MB contiguous chunks written in multi 4KB transactions per IOnode in 64^3 lattices for configuration I/O. I suspect this is fine for system performance.
2015-08-26 13:40:29 +01:00
static RealD linkTrace(const GaugeLorentz &Umu){
std::vector<GaugeMat> U(4,Umu._grid);
LatticeComplex Tr(Umu._grid); Tr=zero;
for(int mu=0;mu<Nd;mu++){
U[mu] = PeekIndex<LorentzIndex>(Umu,mu);
Tr = Tr+trace(U[mu]);
}
TComplex Tp = sum(Tr);
Complex p = TensorRemove(Tp);
double vol = Umu._grid->gSites();
return p.real()/vol/4.0/3.0;
};
Adding some wilson loop support
2015-06-05 10:02:36 +01:00
//////////////////////////////////////////////////
// the sum over all staples on each site
//////////////////////////////////////////////////
const safety
2015-06-14 00:59:50 +01:00
static void Staple(GaugeMat &staple,const GaugeLorentz &Umu,int mu){
Adding some wilson loop support
2015-06-05 10:02:36 +01:00
uninitialised bug fix
2015-06-16 14:07:05 +01:00
GridBase *grid = Umu._grid;
std::vector<GaugeMat> U(4,grid);
Adding some wilson loop support
2015-06-05 10:02:36 +01:00
for(int d=0;d<Nd;d++){
big commit fixing nocompiles in defective C++11 compilers (gcc, icpc). stared getting to near the bleeding edge I guess
2015-06-30 15:03:11 +01:00
U[d] = PeekIndex<LorentzIndex>(Umu,d);
Adding some wilson loop support
2015-06-05 10:02:36 +01:00
}
staple = zero;
uninitialised bug fix
2015-06-16 14:07:05 +01:00
GaugeMat tmp(grid);
Adding some wilson loop support
2015-06-05 10:02:36 +01:00
for(int nu=0;nu<Nd;nu++){
if(nu != mu) {
// mu
// ^
Plaq, Rectangle, Iwasaki, Symanzik and DBW2 workign and HMC regresses to http://arxiv.org/pdf/hep-lat/0610075.pdf
2015-12-28 16:38:31 +00:00
// |__> nu
Adding some wilson loop support
2015-06-05 10:02:36 +01:00
// __ __
// | |
// __| = __| *
//
Plaq, Rectangle, Iwasaki, Symanzik and DBW2 workign and HMC regresses to http://arxiv.org/pdf/hep-lat/0610075.pdf
2015-12-28 16:38:31 +00:00
// staple += CovShiftForward(U[nu],nu,U[mu])*Cshift(adj(U[nu]),mu,+1);
staple +=
Cshift(CovShiftForward (U[nu],nu,
CovShiftBackward(U[mu],mu,Cshift(adj(U[nu]),nu,-1))),mu,1);
// Unu(x) Umu(x+nu) Unu^dag(x+mu) ; left mult by Umu^dag(x) to close ring
Adding some wilson loop support
2015-06-05 10:02:36 +01:00
//
// __ __
// | |
// |__ = | * __
//
//
Plaq, Rectangle, Iwasaki, Symanzik and DBW2 workign and HMC regresses to http://arxiv.org/pdf/hep-lat/0610075.pdf
2015-12-28 16:38:31 +00:00
staple +=
Cshift(
CovShiftBackward(U[nu],nu,
CovShiftBackward(U[mu],mu,U[nu])),mu,1);
// tmp = CovShiftForward (U[mu],mu,U[nu]);
// staple+= CovShiftBackward(U[nu],nu,tmp);
// Unu^dag(x-nu) Umu(x-nu) Unu(x+mu-nu) ; left mult by Umu^dag(x) to close ring.
Adding some wilson loop support
2015-06-05 10:02:36 +01:00
}
}
}
//////////////////////////////////////////////////////
// Similar to above for rectangle is required
//////////////////////////////////////////////////////
Plaq, Rectangle, Iwasaki, Symanzik and DBW2 workign and HMC regresses to http://arxiv.org/pdf/hep-lat/0610075.pdf
2015-12-28 16:38:31 +00:00
static void dirRectangle(GaugeMat &rect,const std::vector<GaugeMat> &U, const int mu, const int nu)
{
rect = CovShiftForward(U[mu],mu,CovShiftForward(U[mu],mu,U[nu]))* // ->->|
adj(CovShiftForward(U[nu],nu,CovShiftForward(U[mu],mu,U[mu]))) ;
rect = rect +
CovShiftForward(U[mu],mu,CovShiftForward(U[nu],nu,U[nu]))* // ->||
adj(CovShiftForward(U[nu],nu,CovShiftForward(U[nu],nu,U[mu]))) ;
}
static void traceDirRectangle(LatticeComplex &rect, const std::vector<GaugeMat> &U, const int mu, const int nu)
{
GaugeMat sp(U[0]._grid);
dirRectangle(sp,U,mu,nu);
rect=trace(sp);
}
static void siteRectangle(LatticeComplex &Rect,const std::vector<GaugeMat> &U)
{
LatticeComplex siteRect(U[0]._grid);
Rect=zero;
for(int mu=1;mu<Nd;mu++){
for(int nu=0;nu<mu;nu++){
traceDirRectangle(siteRect,U,mu,nu);
Rect = Rect + siteRect;
}
}
}
//////////////////////////////////////////////////
// sum over all x,y,z,t and over all planes of plaquette
//////////////////////////////////////////////////
static RealD sumRectangle(const GaugeLorentz &Umu){
std::vector<GaugeMat> U(4,Umu._grid);
for(int mu=0;mu<Nd;mu++){
U[mu] = PeekIndex<LorentzIndex>(Umu,mu);
}
LatticeComplex Rect(Umu._grid);
siteRectangle(Rect,U);
TComplex Tp = sum(Rect);
Complex p = TensorRemove(Tp);
return p.real();
}
//////////////////////////////////////////////////
// average over all x,y,z,t and over all planes of plaquette
//////////////////////////////////////////////////
static RealD avgRectangle(const GaugeLorentz &Umu){
RealD sumrect = sumRectangle(Umu);
double vol = Umu._grid->gSites();
double faces = (1.0*Nd*(Nd-1)); // 2 distinct orientations summed
return sumrect/vol/faces/Nc; // Nd , Nc dependent... FIXME
}
//////////////////////////////////////////////////
// the sum over all staples on each site
//////////////////////////////////////////////////
Optimised version of rectangle term staples. ~3.4x faster than the naive.
2015-12-29 19:22:59 +00:00
static void RectStapleDouble(GaugeMat &U2,const GaugeMat & U,int mu){
U2 = U * Cshift(U,mu,1);
}
static void RectStapleOptimised(GaugeMat &Stap,std::vector<GaugeMat> &U2,std::vector<GaugeMat> &U,int mu){
Stap = zero;
GridBase *grid = U[0]._grid;
GaugeMat Staple2x1 (grid);
GaugeMat tmp (grid);
for(int nu=0;nu<Nd;nu++){
if ( nu!=mu) {
// Up staple ___ ___
// | |
tmp = Cshift(adj(U[nu]),nu,-1);
tmp = adj(U2[mu])*tmp;
tmp = Cshift(tmp,mu,-2);
Staple2x1 = CovShiftForward (U[nu],nu,tmp);
// Down staple
// |___ ___|
//
tmp = adj(U2[mu])*U[nu];
Staple2x1+= CovShiftBackward(U[nu],nu,Cshift(tmp,mu,-2));
// ___ ___
// | ___|
// |___ ___|
//
Stap+= Cshift(CovShiftForward (U[mu],mu,Staple2x1),mu,1);
// ___ ___
// |___ |
// |___ ___|
//
// tmp= Staple2x1* Cshift(U[mu],mu,-2);
// Stap+= Cshift(tmp,mu,1) ;
Stap+= Cshift(Staple2x1,mu,1)*Cshift(U[mu],mu,-1); ;
// --
// | |
//
// | |
tmp = Cshift(adj(U2[nu]),nu,-2);
tmp = CovShiftBackward(U[mu],mu,tmp);
tmp = U2[nu]*Cshift(tmp,nu,2);
Stap+= Cshift(tmp, mu, 1);
// | |
//
// | |
// --
tmp = CovShiftBackward(U[mu],mu,U2[nu]);
tmp = adj(U2[nu])*tmp;
tmp = Cshift(tmp,nu,-2);
Stap+=Cshift(tmp, mu, 1);
}}
}
Plaq, Rectangle, Iwasaki, Symanzik and DBW2 workign and HMC regresses to http://arxiv.org/pdf/hep-lat/0610075.pdf
2015-12-28 16:38:31 +00:00
static void RectStaple(GaugeMat &Stap,const GaugeLorentz &Umu,int mu){
Optimised version of rectangle term staples. ~3.4x faster than the naive.
2015-12-29 19:22:59 +00:00
#if 1
RectStapleUnoptimised(Stap,Umu,mu);
#else
GridBase *grid = Umu._grid;
std::vector<GaugeMat> U (Nd,grid);
std::vector<GaugeMat> U2(Nd,grid);
for(int mu=0;mu<Nd;mu++){
U[mu] = PeekIndex<LorentzIndex>(Umu,mu);
RectStapleDouble(U2[mu],U[mu],mu);
}
RectStapleOptimised(Stap,U2,U,mu);
#endif
}
static void RectStapleUnoptimised(GaugeMat &Stap,const GaugeLorentz &Umu,int mu){
Plaq, Rectangle, Iwasaki, Symanzik and DBW2 workign and HMC regresses to http://arxiv.org/pdf/hep-lat/0610075.pdf
2015-12-28 16:38:31 +00:00
GridBase *grid = Umu._grid;
std::vector<GaugeMat> U(4,grid);
for(int d=0;d<Nd;d++){
U[d] = PeekIndex<LorentzIndex>(Umu,d);
}
Stap=zero;
for(int nu=0;nu<Nd;nu++){
if ( nu!=mu) {
// __ ___
// | __ |
//
Stap+= Cshift(CovShiftForward (U[mu],mu,
CovShiftForward (U[nu],nu,
CovShiftBackward(U[mu],mu,
CovShiftBackward(U[mu],mu,
Cshift(adj(U[nu]),nu,-1))))) , mu, 1);
// __
// |__ __ |
Stap+= Cshift(CovShiftForward (U[mu],mu,
CovShiftBackward(U[nu],nu,
CovShiftBackward(U[mu],mu,
CovShiftBackward(U[mu],mu, U[nu])))) , mu, 1);
// __
// |__ __ |
Stap+= Cshift(CovShiftBackward(U[nu],nu,
CovShiftBackward(U[mu],mu,
CovShiftBackward(U[mu],mu,
CovShiftForward(U[nu],nu,U[mu])))) , mu, 1);
// __ ___
// |__ |
Stap+= Cshift(CovShiftForward (U[nu],nu,
CovShiftBackward(U[mu],mu,
CovShiftBackward(U[mu],mu,
CovShiftBackward(U[nu],nu,U[mu])))) , mu, 1);
// --
// | |
//
// | |
Stap+= Cshift(CovShiftForward(U[nu],nu,
CovShiftForward(U[nu],nu,
CovShiftBackward(U[mu],mu,
CovShiftBackward(U[nu],nu,
Cshift(adj(U[nu]),nu,-1))))) , mu, 1);
// | |
//
// | |
// --
Stap+= Cshift(CovShiftBackward(U[nu],nu,
CovShiftBackward(U[nu],nu,
CovShiftBackward(U[mu],mu,
CovShiftForward (U[nu],nu,U[nu])))) , mu, 1);
}}
}
Adding some wilson loop support
2015-06-05 10:02:36 +01:00
};
Compile fix
2015-06-05 10:29:42 +01:00
Global edit on HMC sector -- making GaugeField a template parameter and preparing to pass integrator, smearing, bc's as policy classes to hmc. Propose to unify "integrator" and integrator algorithm in a base/derived way to override step. Want to read through ForceGradient to ensure that abstraction covers the force gradient case.
2015-08-30 12:18:34 +01:00
typedef WilsonLoops<LatticeGaugeField> ColourWilsonLoops;
typedef WilsonLoops<LatticeGaugeField> U1WilsonLoops;
typedef WilsonLoops<LatticeGaugeField> SU2WilsonLoops;
typedef WilsonLoops<LatticeGaugeField> SU3WilsonLoops;
Compile fix
2015-06-05 10:29:42 +01:00
Adding some wilson loop support
2015-06-05 10:02:36 +01:00
}}
#endif
Copy Permalink