solve: add the shift terms

1 files changed, 537 insertions, 54 deletions

diff --git a/src/qms_solve.c b/src/qms_solve.c
index 7136f38..c84c1d8 100644
--- a/src/qms_solve.c
+++ b/src/qms_solve.c
@@ -73,7 +73,6 @@ enum MetricVars {
     XTDOT_X,
     XTDOT_Y,
     XTDOT_Z,
-    PHIDOT,
     NB_METRIC_VARS,
 };
 
@@ -85,23 +84,66 @@ enum InterpMetricVars {
     I_GTXY,
     I_GTXZ,
     I_GTYZ,
+    I_GTXX_DX,
+    I_GTYY_DX,
+    I_GTZZ_DX,
+    I_GTXZ_DX,
+    I_GTXX_DZ,
+    I_GTYY_DZ,
+    I_GTZZ_DZ,
+    I_GTXZ_DZ,
+    I_GTXX_DXX,
+    I_GTYY_DXX,
+    I_GTZZ_DXX,
+    I_GTXZ_DXX,
+    I_GTXX_DXZ,
+    I_GTYY_DXZ,
+    I_GTZZ_DXZ,
+    I_GTXZ_DXZ,
+    I_GTXX_DZZ,
+    I_GTYY_DZZ,
+    I_GTZZ_DZZ,
+    I_GTXZ_DZZ,
     I_PHI,
     I_PHI_DX,
     I_PHI_DY,
     I_PHI_DZ,
+    I_PHI_DXX,
+    I_PHI_DZZ,
+    I_PHI_DXZ,
     I_ATXX,
     I_ATYY,
     I_ATZZ,
     I_ATXY,
     I_ATXZ,
     I_ATYZ,
+    I_ATXX_DX,
+    I_ATYY_DX,
+    I_ATZZ_DX,
+    I_ATXZ_DX,
+    I_ATXX_DZ,
+    I_ATYY_DZ,
+    I_ATZZ_DZ,
+    I_ATXZ_DZ,
     I_K,
+    I_TRK_DX,
+    I_TRK_DZ,
     I_XTX,
     I_XTY,
     I_XTZ,
     I_BETAX,
     I_BETAY,
     I_BETAZ,
+    I_BETAX_DX,
+    I_BETAX_DZ,
+    I_BETAX_DXX,
+    I_BETAX_DXZ,
+    I_BETAX_DZZ,
+    I_BETAZ_DX,
+    I_BETAZ_DZ,
+    I_BETAZ_DXX,
+    I_BETAZ_DXZ,
+    I_BETAZ_DZZ,
     I_ALPHA,
     I_ALPHA_DX,
     I_ALPHA_DY,
@@ -121,10 +163,6 @@ enum InterpMetricVars {
     I_XTDOT_X,
     I_XTDOT_Y,
     I_XTDOT_Z,
-    I_PHIDOT,
-    I_PHIDOT_DX,
-    I_PHIDOT_DY,
-    I_PHIDOT_DZ,
     NB_INTERP_VARS,
 };
 
@@ -199,7 +237,6 @@ static const char *metric_vars[] = {
     [XTDOT_X] = "ML_CCZ4::Xtdot1",
     [XTDOT_Y] = "ML_CCZ4::Xtdot2",
     [XTDOT_Z] = "ML_CCZ4::Xtdot3",
-    [PHIDOT]  = "ML_CCZ4::phidot",
 #else
     [GTXX]  = "ML_BSSN::gt11",
     [GTYY]  = "ML_BSSN::gt22",
@@ -232,7 +269,6 @@ static const char *metric_vars[] = {
     [XTDOT_X] = "ML_BSSN::Xtdot1",
     [XTDOT_Y] = "ML_BSSN::Xtdot2",
     [XTDOT_Z] = "ML_BSSN::Xtdot3",
-    [PHIDOT]  = "ML_BSSN::phidot",
 #endif
 };
 
@@ -244,23 +280,66 @@ static const CCTK_INT interp_operation_indices[] = {
     [I_GTXY]     = GTXY,
     [I_GTXZ]     = GTXZ,
     [I_GTYZ]     = GTYZ,
+    [I_GTXX_DX]  = GTXX,
+    [I_GTYY_DX]  = GTYY,
+    [I_GTZZ_DX]  = GTZZ,
+    [I_GTXZ_DX]  = GTXZ,
+    [I_GTXX_DZ]  = GTXX,
+    [I_GTYY_DZ]  = GTYY,
+    [I_GTZZ_DZ]  = GTZZ,
+    [I_GTXZ_DZ]  = GTXZ,
+    [I_GTXX_DXX] = GTXX,
+    [I_GTYY_DXX] = GTYY,
+    [I_GTZZ_DXX] = GTZZ,
+    [I_GTXZ_DXX] = GTXZ,
+    [I_GTXX_DXZ] = GTXX,
+    [I_GTYY_DXZ] = GTYY,
+    [I_GTZZ_DXZ] = GTZZ,
+    [I_GTXZ_DXZ] = GTXZ,
+    [I_GTXX_DZZ] = GTXX,
+    [I_GTYY_DZZ] = GTYY,
+    [I_GTZZ_DZZ] = GTZZ,
+    [I_GTXZ_DZZ] = GTXZ,
     [I_PHI]      = PHI,
     [I_PHI_DX]   = PHI,
     [I_PHI_DY]   = PHI,
     [I_PHI_DZ]   = PHI,
+    [I_PHI_DXX]  = PHI,
+    [I_PHI_DZZ]  = PHI,
+    [I_PHI_DXZ]  = PHI,
     [I_ATXX]     = ATXX,
     [I_ATYY]     = ATYY,
     [I_ATZZ]     = ATZZ,
     [I_ATXY]     = ATXY,
     [I_ATXZ]     = ATXZ,
     [I_ATYZ]     = ATYZ,
+    [I_ATXX_DX]  = ATXX,
+    [I_ATYY_DX]  = ATYY,
+    [I_ATZZ_DX]  = ATZZ,
+    [I_ATXZ_DX]  = ATXZ,
+    [I_ATXX_DZ]  = ATXX,
+    [I_ATYY_DZ]  = ATYY,
+    [I_ATZZ_DZ]  = ATZZ,
+    [I_ATXZ_DZ]  = ATXZ,
     [I_K]        = K,
+    [I_TRK_DX]   = K,
+    [I_TRK_DZ]   = K,
     [I_XTX]      = XTX,
     [I_XTY]      = XTY,
     [I_XTZ]      = XTZ,
-    [I_BETAX]    = BETAX,
-    [I_BETAY]    = BETAY,
-    [I_BETAZ]    = BETAZ,
+    [I_BETAX]     = BETAX,
+    [I_BETAY]     = BETAY,
+    [I_BETAZ]     = BETAZ,
+    [I_BETAX_DX]  = BETAX,
+    [I_BETAX_DZ]  = BETAX,
+    [I_BETAX_DXX] = BETAX,
+    [I_BETAX_DXZ] = BETAX,
+    [I_BETAX_DZZ] = BETAX,
+    [I_BETAZ_DX]  = BETAZ,
+    [I_BETAZ_DZ]  = BETAZ,
+    [I_BETAZ_DXX] = BETAZ,
+    [I_BETAZ_DXZ] = BETAZ,
+    [I_BETAZ_DZZ] = BETAZ,
     [I_ALPHA]    = ALPHA,
     [I_ALPHA_DX] = ALPHA,
     [I_ALPHA_DY] = ALPHA,
@@ -280,10 +359,6 @@ static const CCTK_INT interp_operation_indices[] = {
     [I_XTDOT_X]   = XTDOT_X,
     [I_XTDOT_Y]   = XTDOT_Y,
     [I_XTDOT_Z]   = XTDOT_Z,
-    [I_PHIDOT]    = PHIDOT,
-    [I_PHIDOT_DX] = PHIDOT,
-    [I_PHIDOT_DY] = PHIDOT,
-    [I_PHIDOT_DZ] = PHIDOT,
 };
 
 /* the operation (plain value or x/y/z-derivative) to apply during interpolation */
@@ -294,23 +369,66 @@ static const CCTK_INT interp_operation_codes[] = {
     [I_GTXY]     = 0,
     [I_GTXZ]     = 0,
     [I_GTYZ]     = 0,
+    [I_GTXX_DX]  = 1,
+    [I_GTYY_DX]  = 1,
+    [I_GTZZ_DX]  = 1,
+    [I_GTXZ_DX]  = 1,
+    [I_GTXX_DZ]  = 3,
+    [I_GTYY_DZ]  = 3,
+    [I_GTZZ_DZ]  = 3,
+    [I_GTXZ_DZ]  = 3,
+    [I_GTXX_DXX] = 11,
+    [I_GTYY_DXX] = 11,
+    [I_GTZZ_DXX] = 11,
+    [I_GTXZ_DXX] = 11,
+    [I_GTXX_DXZ] = 13,
+    [I_GTYY_DXZ] = 13,
+    [I_GTZZ_DXZ] = 13,
+    [I_GTXZ_DXZ] = 13,
+    [I_GTXX_DZZ] = 33,
+    [I_GTYY_DZZ] = 33,
+    [I_GTZZ_DZZ] = 33,
+    [I_GTXZ_DZZ] = 33,
     [I_PHI]      = 0,
     [I_PHI_DX]   = 1,
     [I_PHI_DY]   = 2,
     [I_PHI_DZ]   = 3,
+    [I_PHI_DXX]  = 11,
+    [I_PHI_DZZ]  = 33,
+    [I_PHI_DXZ]  = 13,
     [I_ATXX]     = 0,
     [I_ATYY]     = 0,
     [I_ATZZ]     = 0,
     [I_ATXY]     = 0,
     [I_ATXZ]     = 0,
     [I_ATYZ]     = 0,
+    [I_ATXX_DX]  = 1,
+    [I_ATYY_DX]  = 1,
+    [I_ATZZ_DX]  = 1,
+    [I_ATXZ_DX]  = 1,
+    [I_ATXX_DZ]  = 3,
+    [I_ATYY_DZ]  = 3,
+    [I_ATZZ_DZ]  = 3,
+    [I_ATXZ_DZ]  = 3,
     [I_K]        = 0,
+    [I_TRK_DX]   = 1,
+    [I_TRK_DZ]   = 3,
     [I_XTX]      = 0,
     [I_XTY]      = 0,
     [I_XTZ]      = 0,
-    [I_BETAX]    = 0,
-    [I_BETAY]    = 0,
-    [I_BETAZ]    = 0,
+    [I_BETAX]     = 0,
+    [I_BETAY]     = 0,
+    [I_BETAZ]     = 0,
+    [I_BETAX_DX]  = 1,
+    [I_BETAX_DZ]  = 3,
+    [I_BETAX_DXX] = 11,
+    [I_BETAX_DXZ] = 13,
+    [I_BETAX_DZZ] = 33,
+    [I_BETAZ_DX]  = 1,
+    [I_BETAZ_DZ]  = 3,
+    [I_BETAZ_DXX] = 11,
+    [I_BETAZ_DXZ] = 13,
+    [I_BETAZ_DZZ] = 33,
     [I_ALPHA]    = 0,
     [I_ALPHA_DX] = 1,
     [I_ALPHA_DY] = 2,
@@ -330,10 +448,6 @@ static const CCTK_INT interp_operation_codes[] = {
     [I_XTDOT_X]   = 0,
     [I_XTDOT_Y]   = 0,
     [I_XTDOT_Z]   = 0,
-    [I_PHIDOT]    = 0,
-    [I_PHIDOT_DX] = 1,
-    [I_PHIDOT_DY] = 2,
-    [I_PHIDOT_DZ] = 3,
 };
 
 /* interpolate the cactus gridfunctions onto the pseudospectral grid */
@@ -363,8 +477,16 @@ static int calc_eq_coeffs(QMSSolver *ctx)
         const double z = s->interp_coords[2][i];
         const int zaxis = x <= EPS;
 
-        double Am[3][3], K[3][3], Km[3][3], Ku[3][3], gtu[3][3];
-        double k2, kij_dij_alpha, k_kdot, k3;
+        double K[3][3], Km[3][3], Ku[3][3], dK[3][3][3];
+        double gtu[3][3], g[3][3], gu[3][3];
+        double dg[3][3][3], dgu[3][3][3], d2g[3][3][3][3], gudot[3][3];
+        double G[3][3][3], dG[3][3][3][3], Gdot[3][3][3];
+        double Ric[3][3], Ricm[3][3];
+        double dgdot[3][3][3];
+        double k2, Kij_Dij_alpha, k_kdot, k3;
+        double Gammadot_term, beta_term;
+
+        double betal[3], dbetal[3][3], d2betal[3][3][3];
 
         const double gtxx = s->interp_values[I_GTXX][i];
         const double gtyy = s->interp_values[I_GTYY][i];
@@ -377,6 +499,111 @@ static int calc_eq_coeffs(QMSSolver *ctx)
                                  { gtxy, gtyy, gtyz },
                                  { gtxz, gtyz, gtzz }};
 
+        const double dx_gt11 = s->interp_values[I_GTXX_DX][i];
+        const double dx_gt22 = s->interp_values[I_GTYY_DX][i];
+        const double dx_gt33 = s->interp_values[I_GTZZ_DX][i];
+        const double dx_gt13 = s->interp_values[I_GTXZ_DX][i];
+
+        const double dz_gt11 = s->interp_values[I_GTXX_DZ][i];
+        const double dz_gt22 = s->interp_values[I_GTYY_DZ][i];
+        const double dz_gt33 = s->interp_values[I_GTZZ_DZ][i];
+        const double dz_gt13 = s->interp_values[I_GTXZ_DZ][i];
+
+        const double dgt[3][3][3] = {
+            {
+                { dx_gt11,     0.0, dx_gt13 },
+                {     0.0, dx_gt22,     0.0 },
+                { dx_gt13,     0.0, dx_gt33 },
+            },
+            {
+                {     0.0, zaxis ? dx_gt11 - dx_gt22 : (gtxx - gtyy) / x, 0.0 },
+                { zaxis ? dx_gt11 - dx_gt22 : (gtxx - gtyy) / x, 0.0, zaxis ? dx_gt13 : gtxz / x },
+                { 0.0, zaxis ? dx_gt13 : gtxz / x, 0.0 },
+            },
+            {
+                { dz_gt11,     0.0, dz_gt13 },
+                {     0.0, dz_gt22,     0.0 },
+                { dz_gt13,     0.0, dz_gt33 },
+            },
+        };
+
+        const double dxx_gt11 = s->interp_values[I_GTXX_DXX][i];
+        const double dxx_gt22 = s->interp_values[I_GTYY_DXX][i];
+        const double dxx_gt33 = s->interp_values[I_GTZZ_DXX][i];
+        const double dxx_gt13 = s->interp_values[I_GTXZ_DXX][i];
+
+        const double dxz_gt11 = s->interp_values[I_GTXX_DXZ][i];
+        const double dxz_gt22 = s->interp_values[I_GTYY_DXZ][i];
+        const double dxz_gt33 = s->interp_values[I_GTZZ_DXZ][i];
+        const double dxz_gt13 = s->interp_values[I_GTXZ_DXZ][i];
+
+        const double dzz_gt11 = s->interp_values[I_GTXX_DZZ][i];
+        const double dzz_gt22 = s->interp_values[I_GTYY_DZZ][i];
+        const double dzz_gt33 = s->interp_values[I_GTZZ_DZZ][i];
+        const double dzz_gt13 = s->interp_values[I_GTXZ_DZZ][i];
+
+        const double d2gt[3][3][3][3] = {
+            {
+                {
+                    { dxx_gt11,      0.0, dxx_gt13 },
+                    { 0.0,      dxx_gt22,      0.0 },
+                    { dxx_gt13,      0.0, dxx_gt33 },
+                },
+                {
+                    {      0.0, zaxis ? 0.5 * (dxx_gt11 - dxx_gt22) : (dx_gt11 - dx_gt22) / x - (gtxx - gtyy) / SQR(x), 0.0 },
+                    { zaxis ? 0.5 * (dxx_gt11 - dxx_gt22) : (dx_gt11 - dx_gt22) / x - (gtxx - gtyy) / SQR(x), 0.0,
+                        zaxis ? 0.5 * dxx_gt13 : dx_gt13 / x - gtxz / SQR(x) },
+                    { 0.0, zaxis ? 0.5 * dxx_gt13 : dx_gt13 / x - gtxz / SQR(x), 0.0 },
+                },
+                {
+                    { dxz_gt11,     0.0, dxz_gt13 },
+                    { 0.0,     dxz_gt22,      0.0 },
+                    { dxz_gt13,     0.0, dxz_gt33 },
+                },
+
+            },
+            {
+                {
+                    { 0.0, zaxis ? 0.5 * (dxx_gt11 - dxx_gt22) : (dx_gt11 - dx_gt22) / x - (gtxx - gtyy) / SQR(x), 0.0 },
+                    { zaxis ? 0.5 * (dxx_gt11 - dxx_gt22) : (dx_gt11 - dx_gt22) / x - (gtxx - gtyy) / SQR(x), 0.0,
+                        zaxis ? 0.5 * dxx_gt13 : dx_gt13 / x - gtxz / SQR(x) },
+                    { 0.0, zaxis ? 0.5 * dxx_gt13 : dx_gt13 / x - gtxz / SQR(x), 0.0 },
+                },
+                {
+                    { zaxis ? dxx_gt22 : dx_gt11 / x - 2 * (gtxx - gtyy) / SQR(x), 0.0,
+                       zaxis ? 0.5 * dxx_gt13 : dx_gt13 / x - gtxz / SQR(x) },
+                    { 0.0, zaxis ? dxx_gt11 : dx_gt22 / x + 2.0 * (gtxx - gtyy) / SQR(x), 0.0 },
+                    { zaxis ? 0.5 * dxx_gt13 : dx_gt13 / x - gtxz / SQR(x), 0.0, zaxis ? dxx_gt33 : dx_gt33 / x },
+                },
+                {
+                    { 0.0, zaxis ? dxz_gt11 - dxz_gt22 : (dz_gt11 - dz_gt22) / x, 0.0 },
+                    { zaxis ? dxz_gt11 - dxz_gt22 : (dz_gt11 - dz_gt22) / x, 0.0,
+                       zaxis ? dxz_gt13 : dz_gt13 / x },
+                    { 0.0, zaxis ? dxz_gt13 : dz_gt13 / x, 0.0 },
+                },
+
+            },
+            {
+                {
+                    { dxz_gt11,      0.0, dxz_gt13 },
+                    {      0.0, dxz_gt22,      0.0 },
+                    { dxz_gt13,      0.0, dxz_gt33 },
+                },
+                {
+                    { 0.0, zaxis ? dxz_gt11 - dxz_gt22 : (dz_gt11 - dz_gt22) / x, 0.0 },
+                    { zaxis ? dxz_gt11 - dxz_gt22 : (dz_gt11 - dz_gt22) / x, 0.0,
+                       zaxis ? dxz_gt13 : dz_gt13 / x },
+                    { 0.0, zaxis ? dxz_gt13 : dz_gt13 / x, 0.0 },
+                },
+                {
+                    { dzz_gt11,      0.0, dzz_gt13 },
+                    {      0.0, dzz_gt22,      0.0 },
+                    { dzz_gt13,      0.0, dzz_gt33 },
+                },
+
+            },
+        };
+
         const double Atxx = s->interp_values[I_ATXX][i];
         const double Atyy = s->interp_values[I_ATYY][i];
         const double Atzz = s->interp_values[I_ATZZ][i];
@@ -384,17 +611,62 @@ static int calc_eq_coeffs(QMSSolver *ctx)
         const double Atxz = s->interp_values[I_ATXZ][i];
         const double Atyz = s->interp_values[I_ATYZ][i];
 
-        const double phi  = s->interp_values[I_PHI][i];
-
-        const double phidot    = s->interp_values[I_PHIDOT][i];
-        const double phidot_dx = s->interp_values[I_PHIDOT_DX][i];
-        const double phidot_dz = s->interp_values[I_PHIDOT_DZ][i];
-
         const double At[3][3] = {{ Atxx, Atxy, Atxz },
                                  { Atxy, Atyy, Atyz },
                                  { Atxz, Atyz, Atzz }};
 
+        const double dx_At11 = s->interp_values[I_ATXX_DX][i];
+        const double dx_At22 = s->interp_values[I_ATYY_DX][i];
+        const double dx_At33 = s->interp_values[I_ATZZ_DX][i];
+        const double dx_At13 = s->interp_values[I_ATXZ_DX][i];
+
+        const double dz_At11 = s->interp_values[I_ATXX_DZ][i];
+        const double dz_At22 = s->interp_values[I_ATYY_DZ][i];
+        const double dz_At33 = s->interp_values[I_ATZZ_DZ][i];
+        const double dz_At13 = s->interp_values[I_ATXZ_DZ][i];
+
+        const double dAt[3][3][3] = {
+            {
+                { dx_At11,     0.0, dx_At13 },
+                {     0.0, dx_At22,     0.0 },
+                { dx_At13,     0.0, dx_At33 },
+            },
+            {
+                {     0.0, zaxis ? dx_At11 - dx_At22 : (Atxx - Atyy) / x, 0.0 },
+                { zaxis ? dx_At11 - dx_At22 : (Atxx - Atyy) / x, 0.0, zaxis ? dx_At13 : Atxz / x },
+                { 0.0, zaxis ? dx_At13 : Atxz / x, 0.0 },
+            },
+            {
+                { dz_At11,     0.0, dz_At13 },
+                {     0.0, dz_At22,     0.0 },
+                { dz_At13,     0.0, dz_At33 },
+            },
+        };
+
+        const double phi  = s->interp_values[I_PHI][i];
+
+        const double dx_phi = s->interp_values[I_PHI_DX][i];
+        const double dz_phi = s->interp_values[I_PHI_DZ][i];
+
+        const double dphi[3] = { dx_phi, 0.0, dz_phi };
+
+        const double phi_dxx = s->interp_values[I_PHI_DXX][i];
+        const double phi_dzz = s->interp_values[I_PHI_DZZ][i];
+        const double phi_dxz = s->interp_values[I_PHI_DXZ][i];
+
+        const double d2phi[3][3] = {
+            { phi_dxx, 0.0, phi_dxz },
+            {     0.0, zaxis ? phi_dxx : dx_phi / x, 0.0 },
+            { phi_dxz,    0.0,  phi_dzz },
+        };
+
         const double trK     = s->interp_values[I_K][i];
+
+        const double dx_trK = s->interp_values[I_TRK_DX][i];
+        const double dz_trK = s->interp_values[I_TRK_DZ][i];
+
+        const double dtrK[3] = { dx_trK, 0.0, dz_trK };
+
         const double kdot_xx = s->interp_values[I_KDOT_XX][i];
         const double kdot_yy = s->interp_values[I_KDOT_YY][i];
         const double kdot_zz = s->interp_values[I_KDOT_ZZ][i];
@@ -409,13 +681,57 @@ static int calc_eq_coeffs(QMSSolver *ctx)
         const double alpha     = s->interp_values[I_ALPHA][i];
         const double dx_alpha  = s->interp_values[I_ALPHA_DX][i];
         const double dz_alpha  = s->interp_values[I_ALPHA_DZ][i];
+
+        const double dalpha[3] = { dx_alpha, 0.0, dz_alpha };
+
         const double dxx_alpha = s->interp_values[I_ALPHA_DXX][i];
         const double dzz_alpha = s->interp_values[I_ALPHA_DZZ][i];
         const double dxz_alpha = s->interp_values[I_ALPHA_DXZ][i];
 
-        const double dij_alpha[3][3] = {{ dxx_alpha,                                0, dxz_alpha },
-                                        {         0, zaxis ? dxx_alpha : dx_alpha / x,         0 },
-                                        { dxz_alpha,                                0, dzz_alpha }};
+        const double d2alpha[3][3] = {{ dxx_alpha,                                0, dxz_alpha },
+                                      {         0, zaxis ? dxx_alpha : dx_alpha / x,         0 },
+                                      { dxz_alpha,                                0, dzz_alpha }};
+
+        const double betax  = s->interp_values[I_BETAX][i];
+        const double betaz  = s->interp_values[I_BETAZ][i];
+
+        const double beta[3] = { betax, 0.0, betaz };
+
+        const double dx_betax = s->interp_values[I_BETAX_DX][i];
+        const double dz_betax = s->interp_values[I_BETAX_DZ][i];
+
+        const double dx_betaz = s->interp_values[I_BETAZ_DX][i];
+        const double dz_betaz = s->interp_values[I_BETAZ_DZ][i];
+
+        const double dbeta[3][3] = {{ dx_betax,      0.0, dx_betaz },
+                                    {      0.0, zaxis ? dx_betax : betax / x,      0.0 },
+                                    { dz_betax,      0.0, dz_betaz }};
+
+        const double dxx_betax = s->interp_values[I_BETAX_DXX][i];
+        const double dxz_betax = s->interp_values[I_BETAX_DXZ][i];
+        const double dzz_betax = s->interp_values[I_BETAX_DZZ][i];
+
+        const double dxx_betaz = s->interp_values[I_BETAZ_DXX][i];
+        const double dxz_betaz = s->interp_values[I_BETAZ_DXZ][i];
+        const double dzz_betaz = s->interp_values[I_BETAZ_DZZ][i];
+
+        const double d2beta[3][3][3] = {
+            {
+                { dxx_betax, 0.0, dxx_betaz },
+                { 0.0, zaxis ? 0.5 * dxx_betax : dx_betax / x - betax / SQR(x), 0.0 },
+                { dxz_betax, 0.0, dxz_betaz },
+            },
+            {
+                { 0.0, zaxis ? 0.5 * dxx_betax : dx_betax / x - betax / SQR(x), 0.0 },
+                { zaxis ? 0.5 * dxx_betax : dx_betax / x - betax / SQR(x), 0.0, zaxis ? dxx_betaz : dx_betaz / x },
+                { 0.0, zaxis ? dxz_betax : dz_betax / x, 0.0 },
+            },
+            {
+                { dxz_betax, 0.0, dxz_betaz },
+                { 0.0, zaxis ? dxz_betax : dz_betax / x, 0.0 },
+                { dzz_betax, 0.0, dzz_betaz },
+            },
+        };
 
         const double Xtx  = s->interp_values[I_XTX][i];
         const double Xtz  = s->interp_values[I_XTZ][i];
@@ -436,16 +752,107 @@ static int calc_eq_coeffs(QMSSolver *ctx)
         gtu[2][0] = gtu[0][2];
         gtu[2][1] = gtu[1][2];
 
+        //　γ_{jk}/^{jk}
+        for (int j = 0; j < 3; j++)
+            for (int k = 0; k < 3; k++) {
+                gu[j][k] = SQR(phi) * gtu[j][k];
+                g[j][k]  = gt[j][k] / SQR(phi);
+            }
+
+        // β_j
+        for (int j = 0; j < 3; j++) {
+            double val = 0.0;
+            for (int k = 0; k < 3; k++)
+                val += g[j][k] * beta[k];
+            betal[j] = val;
+        }
+
+        // ∂_j γ_{kl}
+        for (int j = 0; j < 3; j++)
+            for (int k = 0; k < 3; k++)
+                for (int l = 0; l < 3; l++) {
+                    dg[j][k][l] = -2.0 * dphi[j] * gt[k][l] / (phi * SQR(phi)) + dgt[j][k][l] / SQR(phi);
+                    dK[j][k][l] = -2.0 * dphi[j] * At[k][l] / (phi * SQR(phi)) + dAt[j][k][l] / SQR(phi) +
+                                  (1.0 / 3.0) * (dtrK[j] * g[k][l] + trK * dg[j][k][l]);
+                }
+
+        // ∂_{jk} g_{lm}
+        for (int j = 0; j < 3; j++)
+            for (int k = 0; k < 3; k++)
+                for (int l = 0; l < 3; l++)
+                    for (int m = 0; m < 3; m++) {
+                        d2g[j][k][l][m] = 6.0 *  gt      [l][m] * dphi[j] * dphi[k] / SQR(SQR(phi))    -
+                                          2.0 *  gt      [l][m] * d2phi[j][k]       / (phi * SQR(phi)) -
+                                          2.0 * dgt   [j][l][m] * dphi[k]           / (phi * SQR(phi)) -
+                                          2.0 * dgt   [k][l][m] * dphi[j]           / (phi * SQR(phi)) +
+                                               d2gt[j][k][l][m]                     / SQR(phi);
+                    }
+
+        // ∂_j γ^{kl}
+        for (int j = 0; j < 3; j++)
+            for (int k = 0; k < 3; k++)
+                for (int l = 0; l < 3; l++) {
+                    double val = 0.0;
+                    for (int m = 0; m < 3; m++)
+                        for (int n = 0; n < 3; n++)
+                            val += -gu[k][m] * gu[l][n] * dg[j][m][n];
+                    dgu[j][k][l] = val;
+                }
+
+        // Γ^j_{kl}
+        for (int j = 0; j < 3; j++)
+            for (int k = 0; k < 3; k++)
+                for (int l = 0; l < 3; l++) {
+                    double val = 0.0;
+                    for (int m = 0; m < 3; m++)
+                        val += 0.5 * gu[j][m] * (dg[k][l][m] + dg[l][k][m] - dg[m][k][l]);
+                    G[j][k][l] = val;
+                }
+
+        // ∂_j Γ^k_{lm}
+        for (int j = 0; j < 3; j++)
+            for (int k = 0; k < 3; k++)
+                for (int l = 0; l < 3; l++)
+                    for (int m = 0; m < 3; m++) {
+                        double val = 0.0;
+                        for (int n = 0; n < 3; n++) {
+                            val += dgu[j][k][n] * (dg    [l][m][n] +  dg   [m][l][n] -  dg   [n][l][m]) +
+                                    gu   [k][n] * (d2g[j][l][m][n] + d2g[j][m][l][n] - d2g[j][n][l][m]);
+                        }
+                        dG[j][k][l][m] = 0.5 * val;
+                    }
+
+        // ∂_j β_k
+        for (int j = 0; j < 3; j++)
+            for (int k = 0; k < 3; k++) {
+                double val = 0.0;
+                for (int l = 0; l < 3; l++)
+                    val += g[k][l] * dbeta[j][l] + dg[j][k][l] * beta[l];
+                dbetal[j][k] = val;
+            }
+
+        // ∂_{jk} β_l
+        for (int j = 0; j < 3; j++)
+            for (int k = 0; k < 3; k++)
+                for (int l = 0; l < 3; l++) {
+                    double val = 0.0;
+                    for (int m = 0; m < 3; m++)
+                        val += d2g[j][k][l][m] * beta[m] + dg[k][l][m] * dbeta[j][m] + dg[j][l][m] * dbeta[k][m] +
+                               g[l][m] * d2beta[j][k][m];
+                    d2betal[j][k][l] = val;
+                }
+
+
         // K_{ij}
         for (int j = 0; j < 3; j++)
             for (int k = 0; k < 3; k++)
-                K[j][k] = At[j][k] / SQR(phi) + gt[j][k] * trK;
+                K[j][k] = At[j][k] / SQR(phi) + (1.0 / 3.0) * g[j][k] * trK;
 
         for (int j = 0; j < 3; j++)
             for (int k = 0; k < 3; k++) {
                 double val = 0.0;
                 for (int l = 0; l < 3; l++)
-                    val += SQR(phi) * gtu[j][l] * K[l][k];
+                    val += gu[j][l] * K[l][k];
                 Km[j][k] = val;
             }
 
@@ -454,23 +861,65 @@ static int calc_eq_coeffs(QMSSolver *ctx)
             for (int k = 0; k < 3; k++) {
                 double val = 0.0;
                 for (int l = 0; l < 3; l++)
-                    val += SQR(phi) * gtu[j][l] * Km[k][l];
+                    val += gu[j][l] * Km[k][l];
                 Ku[j][k] = val;
             }
 
-        // \tilde{A}_{i}^j
+        // ∂_j \dot{γ_kl}
+        for (int j = 0; j < 3; j++)
+            for (int k = 0; k < 3; k++)
+                for (int l = 0; l < 3; l++) {
+                    dgdot[j][k][l] = -2.0 * (dalpha[j] * K[k][l] + alpha * dK[j][k][l]) +
+                                     d2betal[j][k][l] + d2betal[j][l][k];
+                    for (int m = 0; m < 3; m++)
+                        dgdot[j][k][l] += -2.0 * (dG[j][m][k][l] * betal[m] + G[m][k][l] * dbetal[j][m]);
+                }
+
+        // \dot{γ^{jk}
         for (int j = 0; j < 3; j++)
             for (int k = 0; k < 3; k++) {
-                double val = 0.0;
-                for (int l = 0; l < 3; l++)
-                    val += gtu[j][l] * At[l][k];
-                Am[j][k] = val;
+                gudot[j][k] = 2.0 * alpha * Ku[j][k];
+
+                for (int l = 0; l < 3; l++) {
+                    gudot[j][k] += -gu[j][l] * dbeta[l][k] - gu[k][l] * dbeta[l][j];
+
+                    for (int m = 0; m < 3; m++)
+                        gudot[j][k] += -gu[j][l] * G[k][l][m] * beta[m] - gu[k][l] * G[j][l][m] * beta[m];
+                }
             }
 
-        kij_dij_alpha = 0.0;
+        // \dot{Γ}^j_{kl}
         for (int j = 0; j < 3; j++)
             for (int k = 0; k < 3; k++)
-                kij_dij_alpha += Ku[j][k] * dij_alpha[j][k];
+                for (int l = 0; l < 3; l++) {
+                    double val = 0.0;
+                    for (int m = 0; m < 3; m++)
+                        val += gudot[j][m] * (dg   [k][l][m] + dg   [l][k][m] - dg   [m][k][l]) +
+                               gu   [j][m] * (dgdot[k][l][m] + dgdot[l][k][m] - dgdot[m][k][l]);
+
+                    Gdot[j][k][l] = 0.5 * val;
+                }
+
+        Gammadot_term = 0.0;
+        for (int j = 0; j < 3; j++) {
+            double val = 0.0;
+            for (int k = 0; k < 3; k++)
+                for (int l = 0; l < 3; l++) {
+                    val += gu[k][l] * Gdot[j][k][l];
+                }
+
+            Gammadot_term += val * dalpha[j];
+        }
+
+        Kij_Dij_alpha = 0.0;
+        for (int j = 0; j < 3; j++)
+            for (int k = 0; k < 3; k++) {
+                double val = 0.0;
+                for (int l = 0; l < 3; l++)
+                    val += G[l][j][k] * dalpha[l];
+
+                Kij_Dij_alpha += Ku[j][k] * (d2alpha[j][k] - val);
+            }
 
         k_kdot = 0.0;
         for (int j = 0; j < 3; j++)
@@ -482,8 +931,8 @@ static int calc_eq_coeffs(QMSSolver *ctx)
             for (int k = 0; k < 3; k++) {
                 double val = 0.0;
                 for (int l = 0; l < 3; l++)
-                    val += Km[k][l] * Ku[l][j];
-                k3 += val * K[j][k];
+                    val += Km[k][l] * Km[l][j];
+                k3 += val * Km[j][k];
             }
 
         // K_{ij} K^{ij}
@@ -492,6 +941,48 @@ static int calc_eq_coeffs(QMSSolver *ctx)
             for (int k = 0; k < 3; k++)
                 k2 += Km[j][k] * Km[k][j];
 
+        // Ric_{jk}
+        for (int j = 0; j < 3; j++)
+            for (int k = 0; k < 3; k++) {
+                double val = 0.0;
+                for (int m = 0; m < 3; m++)
+                    val += dG[m][m][j][k] - dG[k][m][j][m];
+                for (int m = 0; m < 3; m++)
+                    for (int l = 0; l < 3; l++)
+                        val += G[l][l][m] * G[m][j][k] - G[l][k][m] * G[m][j][l];
+                Ric[j][k] = val;
+            }
+
+        // Ric^j_k
+        for (int j = 0; j < 3; j++)
+            for (int k = 0; k < 3; k++) {
+                double val = 0.0;
+                for (int l = 0; l < 3; l++)
+                    val += gu[j][l] * Ric[l][k];
+                Ricm[j][k] = val;
+            }
+
+        beta_term = 0.0;
+        for (int j = 0; j < 3; j++)
+            for (int k = 0; k < 3; k++)
+                for (int l = 0; l < 3; l++) {
+                    double D2beta = 0.0;
+
+                    D2beta += d2beta[j][k][l];
+
+                    for (int m = 0; m < 3; m++) {
+                        D2beta += dG[j][l][k][m] * beta[m] + G[l][k][m] * dbeta[j][m] + G[l][j][m] * dbeta[k][m] - G[m][j][k] * dbeta[m][l];
+
+                        for (int n = 0; n < 3; n++)
+                            D2beta += G[l][j][m] * G[m][k][n] * beta[n] - G[m][j][k] * G[l][m][n] * beta[n];
+                    }
+
+                    beta_term += -gu[j][k] * D2beta * dalpha[l];
+                }
+        for (int j = 0; j < 3; j++)
+            for (int k = 0; k < 3; k++)
+                beta_term += -beta[k] * Ricm[j][k] * dalpha[j];
+
         {
             const double gtuxx         = gtu[0][0];
             const double gtuyy         = gtu[1][1];
@@ -504,17 +995,9 @@ static int calc_eq_coeffs(QMSSolver *ctx)
             const double Xtx           = s->interp_values[I_XTX][i];
             const double Xtz           = s->interp_values[I_XTZ][i];
 
-            const double betax         = s->interp_values[I_BETAX][i];
-            const double betaz         = s->interp_values[I_BETAZ][i];
-
             const double Xx = SQR(phi) * (Xtx + (phi_dx * gtuxx + phi_dz * gtuxz) / phi);
             const double Xz = SQR(phi) * (Xtz + (phi_dx * gtuxz + phi_dz * gtuzz) / phi);
 
-            const double Xdot_x = 2 * phi * phidot * Xtx + SQR(phi) * Xtdot_x + phi * (phidot_dx * gtuxx + phidot_dz * gtuxz) -
-                                  phidot * (phi_dx * gtuxx + phi_dz * gtuxz) + 2 * alpha * (phi_dx * Ku[0][0] + phi_dz * Ku[0][2]) / phi;
-            const double Xdot_z = 2 * phi * phidot * Xtz + SQR(phi) * Xtdot_z + phi * (phidot_dz * gtuzz + phidot_dx * gtuxz) -
-                                  phidot * (phi_dz * gtuzz + phi_dx * gtuxz) + 2 * alpha * (phi_dz * Ku[2][2] + phi_dx * Ku[0][2]) / phi;
-
             s->eq_coeffs[PSSOLVE_DIFF_ORDER_20][i] = SQR(phi) * (gtuxx + ((x <= EPS) ? gtuyy : 0.0));
             s->eq_coeffs[PSSOLVE_DIFF_ORDER_02][i] = SQR(phi) * gtuzz;
             s->eq_coeffs[PSSOLVE_DIFF_ORDER_11][i] = SQR(phi) * gtuxz * 2;
@@ -522,8 +1005,8 @@ static int calc_eq_coeffs(QMSSolver *ctx)
             s->eq_coeffs[PSSOLVE_DIFF_ORDER_01][i] = -Xz;
             s->eq_coeffs[PSSOLVE_DIFF_ORDER_00][i] = -k2;
 
-            s->rhs[i] = -2 * alpha * kij_dij_alpha + Xdot_x * dx_alpha + Xdot_z * dz_alpha +
-                        2 * (k_kdot + 2 * alpha * k3) * alpha;
+            s->rhs[i] = -2 * alpha * Kij_Dij_alpha + Gammadot_term +
+                        2 * alpha * (k_kdot + 2 * alpha * k3) + beta_term;
         }
     }