From fe05bf48b1534f30ed275281195a65e61361fed7 Mon Sep 17 00:00:00 2001 From: Mashy Green Date: Mon, 17 Feb 2025 18:52:04 +0000 Subject: [PATCH] Improvements to WilsonGaugeAction deriv function (#16) * patched version + modifications to deriv -> staple in qcd/gauge * Cleaning up and aligning variable naming between action deriv versions * Removing the regresion test files that were also in this branch for a clean PR * Reverting whitespace changes * Fixing after revering too much! --------- Co-authored-by: Mashy Green --- .../action/gauge/PlaqPlusRectangleAction.h | 16 +++++++------- Grid/qcd/action/gauge/WilsonGaugeAction.h | 17 ++++++++------ Grid/qcd/utils/WilsonLoops.h | 22 +++++++++---------- HMC/FTHMC2p1f.cc | 10 +++++++-- HMC/FTHMC2p1f_3GeV.cc | 9 ++++++++ HMC/HMC2p1f_3GeV.cc | 8 +++++++ 6 files changed, 53 insertions(+), 29 deletions(-) diff --git a/Grid/qcd/action/gauge/PlaqPlusRectangleAction.h b/Grid/qcd/action/gauge/PlaqPlusRectangleAction.h index b9d6ac16..e6ead936 100644 --- a/Grid/qcd/action/gauge/PlaqPlusRectangleAction.h +++ b/Grid/qcd/action/gauge/PlaqPlusRectangleAction.h @@ -71,27 +71,27 @@ public: return action; }; - virtual void deriv(const GaugeField &Umu,GaugeField & dSdU) { + virtual void deriv(const GaugeField &U, GaugeField &dSdU) { //extend Ta to include Lorentz indexes RealD factor_p = c_plaq/RealD(Nc)*0.5; RealD factor_r = c_rect/RealD(Nc)*0.5; - GridBase *grid = Umu.Grid(); + GridBase *grid = U.Grid(); - std::vector U (Nd,grid); + std::vector Umu (Nd,grid); for(int mu=0;mu(Umu,mu); + Umu[mu] = PeekIndex(U,mu); } std::vector RectStaple(Nd,grid), Staple(Nd,grid); - WilsonLoops::StapleAndRectStapleAll(Staple, RectStaple, U, workspace); + WilsonLoops::StapleAndRectStapleAll(Staple, RectStaple, Umu, workspace); GaugeLinkField dSdU_mu(grid); GaugeLinkField staple(grid); for (int mu=0; mu < Nd; mu++){ - dSdU_mu = Ta(U[mu]*Staple[mu])*factor_p; - dSdU_mu = dSdU_mu + Ta(U[mu]*RectStaple[mu])*factor_r; - + dSdU_mu = Ta(Umu[mu]*Staple[mu])*factor_p; + dSdU_mu = dSdU_mu + Ta(Umu[mu]*RectStaple[mu])*factor_r; + PokeIndex(dSdU, dSdU_mu, mu); } diff --git a/Grid/qcd/action/gauge/WilsonGaugeAction.h b/Grid/qcd/action/gauge/WilsonGaugeAction.h index f535b54f..952536b3 100644 --- a/Grid/qcd/action/gauge/WilsonGaugeAction.h +++ b/Grid/qcd/action/gauge/WilsonGaugeAction.h @@ -68,20 +68,23 @@ public: // extend Ta to include Lorentz indexes RealD factor = 0.5 * beta / RealD(Nc); + GridBase *grid = U.Grid(); - GaugeLinkField Umu(U.Grid()); - GaugeLinkField dSdU_mu(U.Grid()); + GaugeLinkField dSdU_mu(grid); + std::vector Umu(Nd, grid); for (int mu = 0; mu < Nd; mu++) { + Umu[mu] = PeekIndex(U, mu); + } - Umu = PeekIndex(U, mu); - + for (int mu = 0; mu < Nd; mu++) { // Staple in direction mu - WilsonLoops::Staple(dSdU_mu, U, mu); - dSdU_mu = Ta(Umu * dSdU_mu) * factor; - + WilsonLoops::Staple(dSdU_mu, Umu, mu); + dSdU_mu = Ta(Umu[mu] * dSdU_mu) * factor; + PokeIndex(dSdU, dSdU_mu, mu); } } + private: RealD beta; }; diff --git a/Grid/qcd/utils/WilsonLoops.h b/Grid/qcd/utils/WilsonLoops.h index 851ba172..655255d6 100644 --- a/Grid/qcd/utils/WilsonLoops.h +++ b/Grid/qcd/utils/WilsonLoops.h @@ -292,18 +292,16 @@ public: ////////////////////////////////////////////////// // the sum over all nu-oriented staples for nu != mu on each site ////////////////////////////////////////////////// - static void Staple(GaugeMat &staple, const GaugeLorentz &Umu, int mu) { + static void Staple(GaugeMat &staple, const GaugeLorentz &U, int mu) { - GridBase *grid = Umu.Grid(); - - std::vector U(Nd, grid); + std::vector Umu(Nd, U.grid()); for (int d = 0; d < Nd; d++) { - U[d] = PeekIndex(Umu, d); + Umu[d] = PeekIndex(U, d); } - Staple(staple, U, mu); + Staple(staple, Umu, mu); } - static void Staple(GaugeMat &staple, const std::vector &U, int mu) { + static void Staple(GaugeMat &staple, const std::vector &Umu, int mu) { staple = Zero(); for (int nu = 0; nu < Nd; nu++) { @@ -318,12 +316,12 @@ public: // | // __| // - + staple += Gimpl::ShiftStaple( Gimpl::CovShiftForward( - U[nu], nu, + Umu[nu], nu, Gimpl::CovShiftBackward( - U[mu], mu, Gimpl::CovShiftIdentityBackward(U[nu], nu))), + Umu[mu], mu, Gimpl::CovShiftIdentityBackward(Umu[nu], nu))), mu); // __ @@ -333,8 +331,8 @@ public: // staple += Gimpl::ShiftStaple( - Gimpl::CovShiftBackward(U[nu], nu, - Gimpl::CovShiftBackward(U[mu], mu, U[nu])), mu); + Gimpl::CovShiftBackward(Umu[nu], nu, + Gimpl::CovShiftBackward(Umu[mu], mu, Umu[nu])), mu); } } } diff --git a/HMC/FTHMC2p1f.cc b/HMC/FTHMC2p1f.cc index 7d93d168..1e914e87 100644 --- a/HMC/FTHMC2p1f.cc +++ b/HMC/FTHMC2p1f.cc @@ -25,13 +25,20 @@ directory *************************************************************************************/ /* END LEGAL */ #include + +#if Nc == 3 #include #include +#endif using namespace Grid; int main(int argc, char **argv) { +#if Nc != 3 +#warning FTHMC2p1f will not work for Nc != 3 + std::cout << "This program will currently only work for Nc == 3." << std::endl; +#else std::cout << std::setprecision(12); Grid_init(&argc, &argv); @@ -220,7 +227,6 @@ int main(int argc, char **argv) TheHMC.Run(SmearingPolicy); // for smearing Grid_finalize(); +#endif } // main - - diff --git a/HMC/FTHMC2p1f_3GeV.cc b/HMC/FTHMC2p1f_3GeV.cc index a8aa67f8..36d5caa3 100644 --- a/HMC/FTHMC2p1f_3GeV.cc +++ b/HMC/FTHMC2p1f_3GeV.cc @@ -24,14 +24,22 @@ See the full license in the file "LICENSE" in the top level distribution directory *************************************************************************************/ /* END LEGAL */ + #include + +#if Nc == 3 #include #include +#endif using namespace Grid; int main(int argc, char **argv) { +#if Nc != 3 +#warning FTHMC2p1f_3GeV will not work for Nc != 3 + std::cout << "This program will currently only work for Nc == 3." << std::endl; +#else std::cout << std::setprecision(12); Grid_init(&argc, &argv); @@ -220,6 +228,7 @@ int main(int argc, char **argv) TheHMC.Run(SmearingPolicy); // for smearing Grid_finalize(); +#endif } // main diff --git a/HMC/HMC2p1f_3GeV.cc b/HMC/HMC2p1f_3GeV.cc index 4bf088d7..199d4be8 100644 --- a/HMC/HMC2p1f_3GeV.cc +++ b/HMC/HMC2p1f_3GeV.cc @@ -25,13 +25,20 @@ directory *************************************************************************************/ /* END LEGAL */ #include + +#if Nc == 3 #include #include +#endif using namespace Grid; int main(int argc, char **argv) { +#if Nc != 3 +#warning HMC2p1f_3GeV will not work for Nc != 3 + std::cout << "This program will currently only work for Nc == 3." << std::endl; +#else std::cout << std::setprecision(12); Grid_init(&argc, &argv); @@ -220,6 +227,7 @@ int main(int argc, char **argv) TheHMC.Run(SmearingPolicy); // for smearing Grid_finalize(); +#endif } // main