portelli/Grid
1
0
Fork 0
mirror of https://github.com/paboyle/Grid.git synced 2025-06-20 16:56:55 +01:00
Code Releases Activity

Added new StapleAll and RectStapleAll functions that return the staples for all mu as an array

Browse Source 
Modified plaq+rectangle gauge actions to use the above
Added a test code to confirm the above changes
This commit is contained in:
Christopher Kelly
2023-06-26 15:48:47 -04:00
parent 4c6613d72c
commit 6f6844ccf1
4 changed files with 225 additions and 26 deletions
Download Patch File Download Diff File Expand all files Collapse all files

18
Grid/qcd/action/gauge/PlaqPlusRectangleAction.h
View File

@ -79,27 +79,19 @@ public:
GridBase *grid = Umu.Grid();
std::vector<GaugeLinkField> U (Nd,grid);
std::vector<GaugeLinkField> U2(Nd,grid);
for(int mu=0;mu<Nd;mu++){
U[mu] = PeekIndex<LorentzIndex>(Umu,mu);
WilsonLoops<Gimpl>::RectStapleDouble(U2[mu],U[mu],mu);
}
std::vector<GaugeLinkField> RectStaple(Nd,grid), Staple(Nd,grid);
WilsonLoops<Gimpl>::StapleAll(Staple, U);
WilsonLoops<Gimpl>::RectStapleAll(RectStaple, U);
GaugeLinkField dSdU_mu(grid);
GaugeLinkField staple(grid);
for (int mu=0; mu < Nd; mu++){
// Staple in direction mu
WilsonLoops<Gimpl>::Staple(staple,Umu,mu);
dSdU_mu = Ta(U[mu]*staple)*factor_p;
WilsonLoops<Gimpl>::RectStaple(Umu,staple,U2,U,mu);
dSdU_mu = dSdU_mu + Ta(U[mu]*staple)*factor_r;
dSdU_mu = Ta(U[mu]*Staple[mu])*factor_p;
dSdU_mu = dSdU_mu + Ta(U[mu]*RectStaple[mu])*factor_r;
PokeIndex<LorentzIndex>(dSdU, dSdU_mu, mu);
}

49
Grid/qcd/utils/WilsonLoops.h
View File

@ -290,7 +290,7 @@ public:
}
*/
//////////////////////////////////////////////////
// the sum over all staples on each site
// the sum over all nu-oriented staples for nu != mu on each site
//////////////////////////////////////////////////
static void Staple(GaugeMat &staple, const GaugeLorentz &Umu, int mu) {
@ -300,6 +300,10 @@ public:
for (int d = 0; d < Nd; d++) {
U[d] = PeekIndex<LorentzIndex>(Umu, d);
}
Staple(staple, U, mu);
}
static void Staple(GaugeMat &staple, const std::vector<GaugeMat> &U, int mu) {
staple = Zero();
for (int nu = 0; nu < Nd; nu++) {
@ -335,6 +339,16 @@ public:
}
}
/////////////
//Staples for each direction mu, summed over nu != mu
//staple: output staples for each mu (Nd)
//U: link array (Nd)
/////////////
static void StapleAll(std::vector<GaugeMat> &staple, const std::vector<GaugeMat> &U) {
assert(staple.size() == Nd); assert(U.size() == Nd);
for(int mu=0;mu<Nd;mu++) Staple(staple[mu], U, mu);
}
//////////////////////////////////////////////////
// the sum over all staples on each site in direction mu,nu, upper part
//////////////////////////////////////////////////
@ -713,8 +727,8 @@ public:
////////////////////////////////////////////////////////////////////////////
// Hop by two optimisation strategy. Use RectStapleDouble to obtain 'U2'
////////////////////////////////////////////////////////////////////////////
static void RectStapleOptimised(GaugeMat &Stap, std::vector<GaugeMat> &U2,
std::vector<GaugeMat> &U, int mu) {
static void RectStapleOptimised(GaugeMat &Stap, const std::vector<GaugeMat> &U2,
const std::vector<GaugeMat> &U, int mu) {
Stap = Zero();
@ -780,15 +794,6 @@ public:
}
}
static void RectStaple(GaugeMat &Stap, const GaugeLorentz &Umu, int mu) {
RectStapleUnoptimised(Stap, Umu, mu);
}
static void RectStaple(const GaugeLorentz &Umu, GaugeMat &Stap,
std::vector<GaugeMat> &U2, std::vector<GaugeMat> &U,
int mu) {
RectStapleOptimised(Stap, U2, U, mu);
}
static void RectStapleUnoptimised(GaugeMat &Stap, const GaugeLorentz &Umu,
int mu) {
GridBase *grid = Umu.Grid();
@ -887,6 +892,26 @@ public:
}
}
static void RectStaple(GaugeMat &Stap, const GaugeLorentz &Umu, int mu) {
RectStapleUnoptimised(Stap, Umu, mu);
}
static void RectStaple(const GaugeLorentz &Umu, GaugeMat &Stap,
std::vector<GaugeMat> &U2, std::vector<GaugeMat> &U,
int mu) {
RectStapleOptimised(Stap, U2, U, mu);
}
//////////////////////////////////////////////////////
//Compute the rectangular staples for all orientations
//Stap : Array of staples (Nd)
//U: Gauge links in each direction (Nd)
/////////////////////////////////////////////////////
static void RectStapleAll(std::vector<GaugeMat> &Stap, const std::vector<GaugeMat> &U){
assert(Stap.size() == Nd); assert(U.size() == Nd);
std::vector<GaugeMat> U2(Nd,U[0].Grid());
for(int mu=0;mu<Nd;mu++) RectStapleDouble(U2[mu], U[mu], mu);
for(int mu=0;mu<Nd;mu++) RectStapleOptimised(Stap[mu], U2, U, mu);
}
//////////////////////////////////////////////////
// Wilson loop of size (R1, R2), oriented in mu,nu plane
//////////////////////////////////////////////////