From a13b6901f6b52f7dc2425190efb4dc0057ae5bab Mon Sep 17 00:00:00 2001 From: ianchb Date: Sun, 26 Apr 2026 15:38:13 +0800 Subject: [PATCH] Fix Z4C C++ gauge damping ordering --- AMSS_NCKU_source/z4c_rhs_c.C | 192 +++++++++++++++++++++++------------ 1 file changed, 128 insertions(+), 64 deletions(-) diff --git a/AMSS_NCKU_source/z4c_rhs_c.C b/AMSS_NCKU_source/z4c_rhs_c.C index b56f821..b679c4c 100644 --- a/AMSS_NCKU_source/z4c_rhs_c.C +++ b/AMSS_NCKU_source/z4c_rhs_c.C @@ -221,6 +221,58 @@ extern "C" void f_z4c_rhs_point( double &rTheta, double &Theta); +static inline void z4c_contract_gamma( + const double gxx, const double gxy, const double gxz, + const double gyy, const double gyz, const double gzz, + const double gxxx, const double gxyx, const double gxzx, + const double gyyx, const double gyzx, const double gzzx, + const double gxxy, const double gxyy, const double gxzy, + const double gyyy, const double gyzy, const double gzzy, + const double gxxz, const double gxyz, const double gxzz, + const double gyyz, const double gyzz, const double gzzz, + double &Gamxa, double &Gamya, double &Gamza) +{ + double det = gxx * gyy * gzz + gxy * gyz * gxz + gxz * gxy * gyz - + gxz * gyy * gxz - gxy * gxy * gzz - gxx * gyz * gyz; + const double gupxx = (gyy * gzz - gyz * gyz) / det; + const double gupxy = -(gxy * gzz - gyz * gxz) / det; + const double gupxz = (gxy * gyz - gyy * gxz) / det; + const double gupyy = (gxx * gzz - gxz * gxz) / det; + const double gupyz = -(gxx * gyz - gxy * gxz) / det; + const double gupzz = (gxx * gyy - gxy * gxy) / det; + + const double Gamxxx = 0.5 * (gupxx * gxxx + gupxy * (2.0 * gxyx - gxxy) + gupxz * (2.0 * gxzx - gxxz)); + const double Gamyxx = 0.5 * (gupxy * gxxx + gupyy * (2.0 * gxyx - gxxy) + gupyz * (2.0 * gxzx - gxxz)); + const double Gamzxx = 0.5 * (gupxz * gxxx + gupyz * (2.0 * gxyx - gxxy) + gupzz * (2.0 * gxzx - gxxz)); + + const double Gamxyy = 0.5 * (gupxx * (2.0 * gxyy - gyyx) + gupxy * gyyy + gupxz * (2.0 * gyzy - gyyz)); + const double Gamyyy = 0.5 * (gupxy * (2.0 * gxyy - gyyx) + gupyy * gyyy + gupyz * (2.0 * gyzy - gyyz)); + const double Gamzyy = 0.5 * (gupxz * (2.0 * gxyy - gyyx) + gupyz * gyyy + gupzz * (2.0 * gyzy - gyyz)); + + const double Gamxzz = 0.5 * (gupxx * (2.0 * gxzz - gzzx) + gupxy * (2.0 * gyzz - gzzy) + gupxz * gzzz); + const double Gamyzz = 0.5 * (gupxy * (2.0 * gxzz - gzzx) + gupyy * (2.0 * gyzz - gzzy) + gupyz * gzzz); + const double Gamzzz = 0.5 * (gupxz * (2.0 * gxzz - gzzx) + gupyz * (2.0 * gyzz - gzzy) + gupzz * gzzz); + + const double Gamxxy = 0.5 * (gupxx * gxxy + gupxy * gyyx + gupxz * (gxzy + gyzx - gxyz)); + const double Gamyxy = 0.5 * (gupxy * gxxy + gupyy * gyyx + gupyz * (gxzy + gyzx - gxyz)); + const double Gamzxy = 0.5 * (gupxz * gxxy + gupyz * gyyx + gupzz * (gxzy + gyzx - gxyz)); + + const double Gamxxz = 0.5 * (gupxx * gxxz + gupxy * (gxyz + gyzx - gxzy) + gupxz * gzzx); + const double Gamyxz = 0.5 * (gupxy * gxxz + gupyy * (gxyz + gyzx - gxzy) + gupyz * gzzx); + const double Gamzxz = 0.5 * (gupxz * gxxz + gupyz * (gxyz + gyzx - gxzy) + gupzz * gzzx); + + const double Gamxyz = 0.5 * (gupxx * (gxyz + gxzy - gyzx) + gupxy * gyyz + gupxz * gzzy); + const double Gamyyz = 0.5 * (gupxy * (gxyz + gxzy - gyzx) + gupyy * gyyz + gupyz * gzzy); + const double Gamzyz = 0.5 * (gupxz * (gxyz + gxzy - gyzx) + gupyz * gyyz + gupzz * gzzy); + + Gamxa = gupxx * Gamxxx + gupyy * Gamxyy + gupzz * Gamxzz + + 2.0 * (gupxy * Gamxxy + gupxz * Gamxxz + gupyz * Gamxyz); + Gamya = gupxx * Gamyxx + gupyy * Gamyyy + gupzz * Gamyzz + + 2.0 * (gupxy * Gamyxy + gupxz * Gamyxz + gupyz * Gamyyz); + Gamza = gupxx * Gamzxx + gupyy * Gamzyy + gupzz * Gamzzz + + 2.0 * (gupxy * Gamzxy + gupxz * Gamzxz + gupyz * Gamzyz); +} + static int compute_rhs_z4c_cartesian( int *ex, double &T, double *X, double *Y, double *Z, double *chi_state, double *chi_constraints, double *trK, @@ -348,6 +400,7 @@ static int compute_rhs_z4c_cartesian( for (int idx = 0; idx < all; ++idx) { + double point_kappa1 = 0.0; f_z4c_rhs_point( Axx[idx], Axy[idx], Axz[idx], Ayy[idx], Ayz[idx], Azz[idx], alpn1[idx], dtSfx[idx], dtSfy[idx], dtSfz[idx], @@ -391,7 +444,7 @@ static int compute_rhs_z4c_cartesian( Gamx[idx], gxx[idx], gxy[idx], gxz[idx], Gamy[idx], gyy[idx], gyz[idx], Gamz[idx], gzz[idx], - kappa1, kappa2, + point_kappa1, kappa2, trK[idx], Axx_rhs[idx], Axy_rhs[idx], Axz_rhs[idx], Ayy_rhs[idx], Ayz_rhs[idx], Azz_rhs[idx], chi_rhs[idx], @@ -427,73 +480,84 @@ static int compute_rhs_z4c_cartesian( #error "z4c_rhs_c.C currently supports GAUGE == 0 or GAUGE == 1 for Z4C" #endif + lopsided(ex, X, Y, Z, gxx, gxx_rhs, betax, betay, betaz, Symmetry, SSS); + lopsided(ex, X, Y, Z, gxy, gxy_rhs, betax, betay, betaz, Symmetry, AAS); + lopsided(ex, X, Y, Z, gxz, gxz_rhs, betax, betay, betaz, Symmetry, ASA); + lopsided(ex, X, Y, Z, gyy, gyy_rhs, betax, betay, betaz, Symmetry, SSS); + lopsided(ex, X, Y, Z, gyz, gyz_rhs, betax, betay, betaz, Symmetry, SAA); + lopsided(ex, X, Y, Z, gzz, gzz_rhs, betax, betay, betaz, Symmetry, SSS); + + lopsided(ex, X, Y, Z, Axx, Axx_rhs, betax, betay, betaz, Symmetry, SSS); + lopsided(ex, X, Y, Z, Axy, Axy_rhs, betax, betay, betaz, Symmetry, AAS); + lopsided(ex, X, Y, Z, Axz, Axz_rhs, betax, betay, betaz, Symmetry, ASA); + lopsided(ex, X, Y, Z, Ayy, Ayy_rhs, betax, betay, betaz, Symmetry, SSS); + lopsided(ex, X, Y, Z, Ayz, Ayz_rhs, betax, betay, betaz, Symmetry, SAA); + lopsided(ex, X, Y, Z, Azz, Azz_rhs, betax, betay, betaz, Symmetry, SSS); + + lopsided(ex, X, Y, Z, chi_state, chi_rhs, betax, betay, betaz, Symmetry, SSS); + lopsided(ex, X, Y, Z, trK, trK_rhs, betax, betay, betaz, Symmetry, SSS); + + lopsided(ex, X, Y, Z, Gamx, Gamx_rhs, betax, betay, betaz, Symmetry, ASS); + lopsided(ex, X, Y, Z, Gamy, Gamy_rhs, betax, betay, betaz, Symmetry, SAS); + lopsided(ex, X, Y, Z, Gamz, Gamz_rhs, betax, betay, betaz, Symmetry, SSA); + + lopsided(ex, X, Y, Z, Lap, Lap_rhs, betax, betay, betaz, Symmetry, SSS); + lopsided(ex, X, Y, Z, betax, betax_rhs, betax, betay, betaz, Symmetry, ASS); + lopsided(ex, X, Y, Z, betay, betay_rhs, betax, betay, betaz, Symmetry, SAS); + lopsided(ex, X, Y, Z, betaz, betaz_rhs, betax, betay, betaz, Symmetry, SSA); +#if (GAUGE == 0) + lopsided(ex, X, Y, Z, dtSfx, dtSfx_rhs, betax, betay, betaz, Symmetry, ASS); + lopsided(ex, X, Y, Z, dtSfy, dtSfy_rhs, betax, betay, betaz, Symmetry, SAS); + lopsided(ex, X, Y, Z, dtSfz, dtSfz_rhs, betax, betay, betaz, Symmetry, SSA); +#endif + lopsided(ex, X, Y, Z, TZ, TZ_rhs, betax, betay, betaz, Symmetry, SSS); + + for (int idx = 0; idx < all; ++idx) + { + double Gamxa = 0.0, Gamya = 0.0, Gamza = 0.0; + z4c_contract_gamma( + gxx[idx], gxy[idx], gxz[idx], gyy[idx], gyz[idx], gzz[idx], + gxxx[idx], gxyx[idx], gxzx[idx], gyyx[idx], gyzx[idx], gzzx[idx], + gxxy[idx], gxyy[idx], gxzy[idx], gyyy[idx], gyzy[idx], gzzy[idx], + gxxz[idx], gxyz[idx], gxzz[idx], gyyz[idx], gyzz[idx], gzzz[idx], + Gamxa, Gamya, Gamza); + + TZ_rhs[idx] -= alpn1[idx] * (TWO + kappa2) * kappa1 * TZ[idx]; + trK_rhs[idx] += alpn1[idx] * kappa1 * (ONE - kappa2) * TZ[idx]; + Gamx_rhs[idx] -= TWO * alpn1[idx] * kappa1 * (Gamx[idx] - Gamxa); + Gamy_rhs[idx] -= TWO * alpn1[idx] * kappa1 * (Gamy[idx] - Gamya); + Gamz_rhs[idx] -= TWO * alpn1[idx] * kappa1 * (Gamz[idx] - Gamza); + } + if (eps > 0.0) { - lopsided_kodis(ex, X, Y, Z, gxx, gxx_rhs, betax, betay, betaz, Symmetry, SSS, eps); - lopsided_kodis(ex, X, Y, Z, gxy, gxy_rhs, betax, betay, betaz, Symmetry, AAS, eps); - lopsided_kodis(ex, X, Y, Z, gxz, gxz_rhs, betax, betay, betaz, Symmetry, ASA, eps); - lopsided_kodis(ex, X, Y, Z, gyy, gyy_rhs, betax, betay, betaz, Symmetry, SSS, eps); - lopsided_kodis(ex, X, Y, Z, gyz, gyz_rhs, betax, betay, betaz, Symmetry, SAA, eps); - lopsided_kodis(ex, X, Y, Z, gzz, gzz_rhs, betax, betay, betaz, Symmetry, SSS, eps); - - lopsided_kodis(ex, X, Y, Z, Axx, Axx_rhs, betax, betay, betaz, Symmetry, SSS, eps); - lopsided_kodis(ex, X, Y, Z, Axy, Axy_rhs, betax, betay, betaz, Symmetry, AAS, eps); - lopsided_kodis(ex, X, Y, Z, Axz, Axz_rhs, betax, betay, betaz, Symmetry, ASA, eps); - lopsided_kodis(ex, X, Y, Z, Ayy, Ayy_rhs, betax, betay, betaz, Symmetry, SSS, eps); - lopsided_kodis(ex, X, Y, Z, Ayz, Ayz_rhs, betax, betay, betaz, Symmetry, SAA, eps); - lopsided_kodis(ex, X, Y, Z, Azz, Azz_rhs, betax, betay, betaz, Symmetry, SSS, eps); - - lopsided_kodis(ex, X, Y, Z, chi_state, chi_rhs, betax, betay, betaz, Symmetry, SSS, eps); - lopsided_kodis(ex, X, Y, Z, trK, trK_rhs, betax, betay, betaz, Symmetry, SSS, eps); - - lopsided_kodis(ex, X, Y, Z, Gamx, Gamx_rhs, betax, betay, betaz, Symmetry, ASS, eps); - lopsided_kodis(ex, X, Y, Z, Gamy, Gamy_rhs, betax, betay, betaz, Symmetry, SAS, eps); - lopsided_kodis(ex, X, Y, Z, Gamz, Gamz_rhs, betax, betay, betaz, Symmetry, SSA, eps); - - lopsided_kodis(ex, X, Y, Z, Lap, Lap_rhs, betax, betay, betaz, Symmetry, SSS, eps); - lopsided_kodis(ex, X, Y, Z, betax, betax_rhs, betax, betay, betaz, Symmetry, ASS, eps); - lopsided_kodis(ex, X, Y, Z, betay, betay_rhs, betax, betay, betaz, Symmetry, SAS, eps); - lopsided_kodis(ex, X, Y, Z, betaz, betaz_rhs, betax, betay, betaz, Symmetry, SSA, eps); + kodis(ex, X, Y, Z, chi_state, chi_rhs, SSS, Symmetry, eps); + kodis(ex, X, Y, Z, trK, trK_rhs, SSS, Symmetry, eps); + kodis(ex, X, Y, Z, gxx, gxx_rhs, SSS, Symmetry, eps); + kodis(ex, X, Y, Z, gxy, gxy_rhs, AAS, Symmetry, eps); + kodis(ex, X, Y, Z, gxz, gxz_rhs, ASA, Symmetry, eps); + kodis(ex, X, Y, Z, gyy, gyy_rhs, SSS, Symmetry, eps); + kodis(ex, X, Y, Z, gyz, gyz_rhs, SAA, Symmetry, eps); + kodis(ex, X, Y, Z, gzz, gzz_rhs, SSS, Symmetry, eps); + kodis(ex, X, Y, Z, Axx, Axx_rhs, SSS, Symmetry, eps); + kodis(ex, X, Y, Z, Axy, Axy_rhs, AAS, Symmetry, eps); + kodis(ex, X, Y, Z, Axz, Axz_rhs, ASA, Symmetry, eps); + kodis(ex, X, Y, Z, Ayy, Ayy_rhs, SSS, Symmetry, eps); + kodis(ex, X, Y, Z, Ayz, Ayz_rhs, SAA, Symmetry, eps); + kodis(ex, X, Y, Z, Azz, Azz_rhs, SSS, Symmetry, eps); + kodis(ex, X, Y, Z, Gamx, Gamx_rhs, ASS, Symmetry, eps); + kodis(ex, X, Y, Z, Gamy, Gamy_rhs, SAS, Symmetry, eps); + kodis(ex, X, Y, Z, Gamz, Gamz_rhs, SSA, Symmetry, eps); + kodis(ex, X, Y, Z, Lap, Lap_rhs, SSS, Symmetry, eps); + kodis(ex, X, Y, Z, betax, betax_rhs, ASS, Symmetry, eps); + kodis(ex, X, Y, Z, betay, betay_rhs, SAS, Symmetry, eps); + kodis(ex, X, Y, Z, betaz, betaz_rhs, SSA, Symmetry, eps); #if (GAUGE == 0) - lopsided_kodis(ex, X, Y, Z, dtSfx, dtSfx_rhs, betax, betay, betaz, Symmetry, ASS, eps); - lopsided_kodis(ex, X, Y, Z, dtSfy, dtSfy_rhs, betax, betay, betaz, Symmetry, SAS, eps); - lopsided_kodis(ex, X, Y, Z, dtSfz, dtSfz_rhs, betax, betay, betaz, Symmetry, SSA, eps); + kodis(ex, X, Y, Z, dtSfx, dtSfx_rhs, ASS, Symmetry, eps); + kodis(ex, X, Y, Z, dtSfy, dtSfy_rhs, SAS, Symmetry, eps); + kodis(ex, X, Y, Z, dtSfz, dtSfz_rhs, SSA, Symmetry, eps); #endif - lopsided_kodis(ex, X, Y, Z, TZ, TZ_rhs, betax, betay, betaz, Symmetry, SSS, eps); - } - else - { - lopsided(ex, X, Y, Z, gxx, gxx_rhs, betax, betay, betaz, Symmetry, SSS); - lopsided(ex, X, Y, Z, gxy, gxy_rhs, betax, betay, betaz, Symmetry, AAS); - lopsided(ex, X, Y, Z, gxz, gxz_rhs, betax, betay, betaz, Symmetry, ASA); - lopsided(ex, X, Y, Z, gyy, gyy_rhs, betax, betay, betaz, Symmetry, SSS); - lopsided(ex, X, Y, Z, gyz, gyz_rhs, betax, betay, betaz, Symmetry, SAA); - lopsided(ex, X, Y, Z, gzz, gzz_rhs, betax, betay, betaz, Symmetry, SSS); - - lopsided(ex, X, Y, Z, Axx, Axx_rhs, betax, betay, betaz, Symmetry, SSS); - lopsided(ex, X, Y, Z, Axy, Axy_rhs, betax, betay, betaz, Symmetry, AAS); - lopsided(ex, X, Y, Z, Axz, Axz_rhs, betax, betay, betaz, Symmetry, ASA); - lopsided(ex, X, Y, Z, Ayy, Ayy_rhs, betax, betay, betaz, Symmetry, SSS); - lopsided(ex, X, Y, Z, Ayz, Ayz_rhs, betax, betay, betaz, Symmetry, SAA); - lopsided(ex, X, Y, Z, Azz, Azz_rhs, betax, betay, betaz, Symmetry, SSS); - - lopsided(ex, X, Y, Z, chi_state, chi_rhs, betax, betay, betaz, Symmetry, SSS); - lopsided(ex, X, Y, Z, trK, trK_rhs, betax, betay, betaz, Symmetry, SSS); - - lopsided(ex, X, Y, Z, Gamx, Gamx_rhs, betax, betay, betaz, Symmetry, ASS); - lopsided(ex, X, Y, Z, Gamy, Gamy_rhs, betax, betay, betaz, Symmetry, SAS); - lopsided(ex, X, Y, Z, Gamz, Gamz_rhs, betax, betay, betaz, Symmetry, SSA); - - lopsided(ex, X, Y, Z, Lap, Lap_rhs, betax, betay, betaz, Symmetry, SSS); - lopsided(ex, X, Y, Z, betax, betax_rhs, betax, betay, betaz, Symmetry, ASS); - lopsided(ex, X, Y, Z, betay, betay_rhs, betax, betay, betaz, Symmetry, SAS); - lopsided(ex, X, Y, Z, betaz, betaz_rhs, betax, betay, betaz, Symmetry, SSA); -#if (GAUGE == 0) - lopsided(ex, X, Y, Z, dtSfx, dtSfx_rhs, betax, betay, betaz, Symmetry, ASS); - lopsided(ex, X, Y, Z, dtSfy, dtSfy_rhs, betax, betay, betaz, Symmetry, SAS); - lopsided(ex, X, Y, Z, dtSfz, dtSfz_rhs, betax, betay, betaz, Symmetry, SSA); -#endif - lopsided(ex, X, Y, Z, TZ, TZ_rhs, betax, betay, betaz, Symmetry, SSS); + kodis(ex, X, Y, Z, TZ, TZ_rhs, SSS, Symmetry, eps); } if (co == 0)