From 00f24f87653d92276f9ffeeaf403f78b5adef4cc Mon Sep 17 00:00:00 2001 From: david clarke Date: Mon, 22 Jan 2024 05:50:16 -0700 Subject: [PATCH] already found some bugs in projection, still needs testing --- Grid/qcd/smearing/HISQSmearing.h | 38 ++++++++++++++++++-------------- 1 file changed, 22 insertions(+), 16 deletions(-) diff --git a/Grid/qcd/smearing/HISQSmearing.h b/Grid/qcd/smearing/HISQSmearing.h index 43b06534..17959495 100644 --- a/Grid/qcd/smearing/HISQSmearing.h +++ b/Grid/qcd/smearing/HISQSmearing.h @@ -93,6 +93,7 @@ public: INHERIT_GIMPL_TYPES(Gimpl); typedef typename Gimpl::GaugeField GF; typedef typename Gimpl::GaugeLinkField LF; + typedef typename Gimpl::ComplexField CF; // Don't allow default values here. Smear_HISQ(GridCartesian* grid, Real c1, Real cnaik, Real c3, Real c5, Real c7, Real clp) @@ -329,31 +330,36 @@ public: // IN--u_mu void projectU3(GF& u_proj, GF& u_mu) const { - LF V, Q, sqrtQinv; - Real c1, c2, g0, g1, g2, S, R, theta, u, v, w, den, f0, f1, f2; + LF V(u_mu.Grid()), Q(u_mu.Grid()), sqrtQinv(u_mu.Grid()), id_3(u_mu.Grid()); + CF c0(u_mu.Grid()), c1(u_mu.Grid()), c2(u_mu.Grid()), g0(u_mu.Grid()), g1(u_mu.Grid()), + g2(u_mu.Grid()), S(u_mu.Grid()), R(u_mu.Grid()), theta(u_mu.Grid()), u(u_mu.Grid()), + v(u_mu.Grid()), w(u_mu.Grid()), den(u_mu.Grid()), f0(u_mu.Grid()), f1(u_mu.Grid()), + f2(u_mu.Grid()); // Follow MILC 10.1103/PhysRevD.82.074501, eqs (B2-B3) and (C1-C8) for (int mu = 0; mu < Nd; mu++) { V = PeekIndex(u_mu, mu); - Q = adj(V[mu])*V[mu]; - c1 = trace(Q*Q)/2.; // SU(N) matrices are traceless, so c0=0. - c2 = trace(Q*Q*Q)/3.; - S = c1/3.; - R = c2/2.; - theta = std::acos(R/std::pow(S,1.5)); - g0 = 2.*std::sqrt(S)*std::cos(theta/3.-2*M_PI/3.); - g1 = 2.*std::sqrt(S)*std::cos(theta/3. ); - g2 = 2.*std::sqrt(S)*std::cos(theta/3.+2*M_PI/3.); + Q = adj(V)*V; + c0 = real(trace(Q)); + c1 = (1/2.)*real(trace(Q*Q)); + c2 = (1/3.)*real(trace(Q*Q*Q)); + S = (1/3.)*c1-(1/18.)*c0*c0; + R = (1/2.)*c2-(1/3. )*c0*c1+(1/27.)*c0*c0*c0; + theta = acos(R*pow(S,-1.5)); + g0 = (1/3.)*c0+2.*sqrt(S)*cos((1/3.)*theta-2*M_PI/3.); + g1 = (1/3.)*c0+2.*sqrt(S)*cos((1/3.)*theta ); + g2 = (1/3.)*c0+2.*sqrt(S)*cos((1/3.)*theta+2*M_PI/3.); // if (fabs(Q.determinant()/(g0*g1*g2)-1.0) > 1e-5) {} - u = std::sqrt(g0) + std::sqrt(g1) + std::sqrt(g2); - v = std::sqrt(g0*g1) + std::sqrt(g0*g2) + std::sqrt(g1*g2); - w = std::sqrt(g0*g1*g2); + u = sqrt(g0) + sqrt(g1) + sqrt(g2); + v = sqrt(g0*g1) + sqrt(g0*g2) + sqrt(g1*g2); + w = sqrt(g0*g1*g2); den = w*(u*v-w); f0 = (-w*(u*u+v)+u*v*v)/den; - f1 = (-w-u*u*u+2*u*v)/den; + f1 = (-w-u*u*u+2.*u*v)/den; f2 = u/den; + id_3 = 1.; - sqrtQinv = f0 + f1*Q + f2*Q*Q; + sqrtQinv = f0*id_3 + f1*Q + f2*Q*Q; PokeIndex(u_proj, V*sqrtQinv, mu); } };