fill_eq_coeffs_template: reindent

1 files changed, 467 insertions, 467 deletions

diff --git a/src/fill_eq_coeffs_template.c b/src/fill_eq_coeffs_template.c
index ae8fcd8..20706f6 100644
--- a/src/fill_eq_coeffs_template.c
+++ b/src/fill_eq_coeffs_template.c
@@ -110,26 +110,26 @@ static void fill_eq_coeffs_line(const cGH *gh, const CoordPatch *cp, MG2DContext
                                 const double * const vars[], const int idx_z,
                                 const ptrdiff_t stride_z, const double dx_inv, const double dz_inv)
 {
-        for (int idx_x = 0; idx_x < solver->local_size[0]; idx_x++) {
-            const int idx_src = CCTK_GFINDEX3D(gh, idx_x + cp->offset_left[0], cp->y_idx, idx_z + cp->offset_left[1]);
-            const int idx_dc  = idx_z * solver->diff_coeffs[0]->stride + idx_x;
-            const int idx_rhs = idx_z * solver->rhs_stride + idx_x;
-
-            const double x = vars[RHS_VAR_X][idx_src];
-            const int on_axis = fabs(x) < 1e-13;
-
-            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], X[3];
-            double Ric[3][3], Ricm[3][3];
-            double dgdot[3][3][3];
-            double D2alpha[3][3];
-            double Kdot[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];
+    for (int idx_x = 0; idx_x < solver->local_size[0]; idx_x++) {
+        const int idx_src = CCTK_GFINDEX3D(gh, idx_x + cp->offset_left[0], cp->y_idx, idx_z + cp->offset_left[1]);
+        const int idx_dc  = idx_z * solver->diff_coeffs[0]->stride + idx_x;
+        const int idx_rhs = idx_z * solver->rhs_stride + idx_x;
+
+        const double x = vars[RHS_VAR_X][idx_src];
+        const int on_axis = fabs(x) < 1e-13;
+
+        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], X[3];
+        double Ric[3][3], Ricm[3][3];
+        double dgdot[3][3][3];
+        double D2alpha[3][3];
+        double Kdot[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];
 
 #define LOAD_VAR(var)  (        vars[RHS_VAR_ ## var][idx_src])
 #define DX(var)        (diff1 (&vars[RHS_VAR_ ## var][idx_src],        1,           dx_inv        ))
@@ -138,217 +138,217 @@ static void fill_eq_coeffs_line(const cGH *gh, const CoordPatch *cp, MG2DContext
 #define DZZ(var)       (diff2 (&vars[RHS_VAR_ ## var][idx_src],           stride_z,         dz_inv))
 #define DXZ(var)       (diff11(&vars[RHS_VAR_ ## var][idx_src],        1, stride_z, dx_inv, dz_inv))
 
-            const double gtxx  = LOAD_VAR(GTXX);
-            const double gtyy  = LOAD_VAR(GTYY);
-            const double gtzz  = LOAD_VAR(GTZZ);
-            const double gtxy  = LOAD_VAR(GTXY);
-            const double gtxz  = LOAD_VAR(GTXZ);
-            const double gtyz  = LOAD_VAR(GTYZ);
-
-            const double gt[3][3] = {{ gtxx, gtxy, gtxz },
-                                     { gtxy, gtyy, gtyz },
-                                     { gtxz, gtyz, gtzz }};
-
-            const double dx_gt11 = DX(GTXX);
-            const double dx_gt22 = DX(GTYY);
-            const double dx_gt33 = DX(GTZZ);
-            const double dx_gt13 = DX(GTXZ);
-
-            const double dz_gt11 = DZ(GTXX);
-            const double dz_gt22 = DZ(GTYY);
-            const double dz_gt33 = DZ(GTZZ);
-            const double dz_gt13 = DZ(GTXZ);
-
-            const double dgt[3][3][3] = {
+        const double gtxx  = LOAD_VAR(GTXX);
+        const double gtyy  = LOAD_VAR(GTYY);
+        const double gtzz  = LOAD_VAR(GTZZ);
+        const double gtxy  = LOAD_VAR(GTXY);
+        const double gtxz  = LOAD_VAR(GTXZ);
+        const double gtyz  = LOAD_VAR(GTYZ);
+
+        const double gt[3][3] = {{ gtxx, gtxy, gtxz },
+                                 { gtxy, gtyy, gtyz },
+                                 { gtxz, gtyz, gtzz }};
+
+        const double dx_gt11 = DX(GTXX);
+        const double dx_gt22 = DX(GTYY);
+        const double dx_gt33 = DX(GTZZ);
+        const double dx_gt13 = DX(GTXZ);
+
+        const double dz_gt11 = DZ(GTXX);
+        const double dz_gt22 = DZ(GTYY);
+        const double dz_gt33 = DZ(GTZZ);
+        const double dz_gt13 = DZ(GTXZ);
+
+        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, on_axis ? dx_gt11 - dx_gt22 : (gtxx - gtyy) / x, 0.0 },
+                { on_axis ? dx_gt11 - dx_gt22 : (gtxx - gtyy) / x, 0.0, on_axis ? dx_gt13 : gtxz / x },
+                { 0.0, on_axis ? 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 = DXX(GTXX);
+        const double dxx_gt22 = DXX(GTYY);
+        const double dxx_gt33 = DXX(GTZZ);
+        const double dxx_gt13 = DXX(GTXZ);
+
+        const double dxz_gt11 = DXZ(GTXX);
+        const double dxz_gt22 = DXZ(GTYY);
+        const double dxz_gt33 = DXZ(GTZZ);
+        const double dxz_gt13 = DXZ(GTXZ);
+
+        const double dzz_gt11 = DZZ(GTXX);
+        const double dzz_gt22 = DZZ(GTYY);
+        const double dzz_gt33 = DZZ(GTZZ);
+        const double dzz_gt13 = DZZ(GTXZ);
+
+        const double d2gt[3][3][3][3] = {
+            {
                 {
-                    { dx_gt11,     0.0, dx_gt13 },
-                    {     0.0, dx_gt22,     0.0 },
-                    { dx_gt13,     0.0, dx_gt33 },
+                    { dxx_gt11,      0.0, dxx_gt13 },
+                    { 0.0,      dxx_gt22,      0.0 },
+                    { dxx_gt13,      0.0, dxx_gt33 },
                 },
                 {
-                    {     0.0, on_axis ? dx_gt11 - dx_gt22 : (gtxx - gtyy) / x, 0.0 },
-                    { on_axis ? dx_gt11 - dx_gt22 : (gtxx - gtyy) / x, 0.0, on_axis ? dx_gt13 : gtxz / x },
-                    { 0.0, on_axis ? dx_gt13 : gtxz / x, 0.0 },
+                    {      0.0, on_axis ? 0.5 * (dxx_gt11 - dxx_gt22) : (dx_gt11 - dx_gt22) / x - (gtxx - gtyy) / SQR(x), 0.0 },
+                    { on_axis ? 0.5 * (dxx_gt11 - dxx_gt22) : (dx_gt11 - dx_gt22) / x - (gtxx - gtyy) / SQR(x), 0.0,
+                        on_axis ? 0.5 * dxx_gt13 : dx_gt13 / x - gtxz / SQR(x) },
+                    { 0.0, on_axis ? 0.5 * dxx_gt13 : dx_gt13 / x - gtxz / SQR(x), 0.0 },
                 },
                 {
-                    { dz_gt11,     0.0, dz_gt13 },
-                    {     0.0, dz_gt22,     0.0 },
-                    { dz_gt13,     0.0, dz_gt33 },
+                    { dxz_gt11,     0.0, dxz_gt13 },
+                    { 0.0,     dxz_gt22,      0.0 },
+                    { dxz_gt13,     0.0, dxz_gt33 },
                 },
-            };
 
-            const double dxx_gt11 = DXX(GTXX);
-            const double dxx_gt22 = DXX(GTYY);
-            const double dxx_gt33 = DXX(GTZZ);
-            const double dxx_gt13 = DXX(GTXZ);
-
-            const double dxz_gt11 = DXZ(GTXX);
-            const double dxz_gt22 = DXZ(GTYY);
-            const double dxz_gt33 = DXZ(GTZZ);
-            const double dxz_gt13 = DXZ(GTXZ);
-
-            const double dzz_gt11 = DZZ(GTXX);
-            const double dzz_gt22 = DZZ(GTYY);
-            const double dzz_gt33 = DZZ(GTZZ);
-            const double dzz_gt13 = DZZ(GTXZ);
-
-            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, on_axis ? 0.5 * (dxx_gt11 - dxx_gt22) : (dx_gt11 - dx_gt22) / x - (gtxx - gtyy) / SQR(x), 0.0 },
-                        { on_axis ? 0.5 * (dxx_gt11 - dxx_gt22) : (dx_gt11 - dx_gt22) / x - (gtxx - gtyy) / SQR(x), 0.0,
-                            on_axis ? 0.5 * dxx_gt13 : dx_gt13 / x - gtxz / SQR(x) },
-                        { 0.0, on_axis ? 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, on_axis ? 0.5 * (dxx_gt11 - dxx_gt22) : (dx_gt11 - dx_gt22) / x - (gtxx - gtyy) / SQR(x), 0.0 },
+                    { on_axis ? 0.5 * (dxx_gt11 - dxx_gt22) : (dx_gt11 - dx_gt22) / x - (gtxx - gtyy) / SQR(x), 0.0,
+                        on_axis ? 0.5 * dxx_gt13 : dx_gt13 / x - gtxz / SQR(x) },
+                    { 0.0, on_axis ? 0.5 * dxx_gt13 : dx_gt13 / x - gtxz / SQR(x), 0.0 },
                 },
                 {
-                    {
-                        { 0.0, on_axis ? 0.5 * (dxx_gt11 - dxx_gt22) : (dx_gt11 - dx_gt22) / x - (gtxx - gtyy) / SQR(x), 0.0 },
-                        { on_axis ? 0.5 * (dxx_gt11 - dxx_gt22) : (dx_gt11 - dx_gt22) / x - (gtxx - gtyy) / SQR(x), 0.0,
-                            on_axis ? 0.5 * dxx_gt13 : dx_gt13 / x - gtxz / SQR(x) },
-                        { 0.0, on_axis ? 0.5 * dxx_gt13 : dx_gt13 / x - gtxz / SQR(x), 0.0 },
-                    },
-                    {
-                        { on_axis ? dxx_gt22 : dx_gt11 / x - 2 * (gtxx - gtyy) / SQR(x), 0.0,
-                           on_axis ? 0.5 * dxx_gt13 : dx_gt13 / x - gtxz / SQR(x) },
-                        { 0.0, on_axis ? dxx_gt11 : dx_gt22 / x + 2.0 * (gtxx - gtyy) / SQR(x), 0.0 },
-                        { on_axis ? 0.5 * dxx_gt13 : dx_gt13 / x - gtxz / SQR(x), 0.0, on_axis ? dxx_gt33 : dx_gt33 / x },
-                    },
-                    {
-                        { 0.0, on_axis ? dxz_gt11 - dxz_gt22 : (dz_gt11 - dz_gt22) / x, 0.0 },
-                        { on_axis ? dxz_gt11 - dxz_gt22 : (dz_gt11 - dz_gt22) / x, 0.0,
-                           on_axis ? dxz_gt13 : dz_gt13 / x },
-                        { 0.0, on_axis ? dxz_gt13 : dz_gt13 / x, 0.0 },
-                    },
-
+                    { on_axis ? dxx_gt22 : dx_gt11 / x - 2 * (gtxx - gtyy) / SQR(x), 0.0,
+                       on_axis ? 0.5 * dxx_gt13 : dx_gt13 / x - gtxz / SQR(x) },
+                    { 0.0, on_axis ? dxx_gt11 : dx_gt22 / x + 2.0 * (gtxx - gtyy) / SQR(x), 0.0 },
+                    { on_axis ? 0.5 * dxx_gt13 : dx_gt13 / x - gtxz / SQR(x), 0.0, on_axis ? dxx_gt33 : dx_gt33 / x },
                 },
                 {
-                    {
-                        { dxz_gt11,      0.0, dxz_gt13 },
-                        {      0.0, dxz_gt22,      0.0 },
-                        { dxz_gt13,      0.0, dxz_gt33 },
-                    },
-                    {
-                        { 0.0, on_axis ? dxz_gt11 - dxz_gt22 : (dz_gt11 - dz_gt22) / x, 0.0 },
-                        { on_axis ? dxz_gt11 - dxz_gt22 : (dz_gt11 - dz_gt22) / x, 0.0,
-                           on_axis ? dxz_gt13 : dz_gt13 / x },
-                        { 0.0, on_axis ? 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 },
-                    },
-
+                    { 0.0, on_axis ? dxz_gt11 - dxz_gt22 : (dz_gt11 - dz_gt22) / x, 0.0 },
+                    { on_axis ? dxz_gt11 - dxz_gt22 : (dz_gt11 - dz_gt22) / x, 0.0,
+                       on_axis ? dxz_gt13 : dz_gt13 / x },
+                    { 0.0, on_axis ? dxz_gt13 : dz_gt13 / x, 0.0 },
                 },
-            };
-
-            const double Atxx  = LOAD_VAR(ATXX);
-            const double Atyy  = LOAD_VAR(ATYY);
-            const double Atzz  = LOAD_VAR(ATZZ);
-            const double Atxy  = LOAD_VAR(ATXY);
-            const double Atxz  = LOAD_VAR(ATXZ);
-            const double Atyz  = LOAD_VAR(ATYZ);
-
-            const double At[3][3] = {
-                { Atxx, Atxy, Atxz },
-                { Atxy, Atyy, Atyz },
-                { Atxz, Atyz, Atzz }};
-
-            const double dx_At11 = DX(ATXX);
-            const double dx_At22 = DX(ATYY);
-            const double dx_At33 = DX(ATZZ);
-            const double dx_At13 = DX(ATXZ);
-
-            const double dz_At11 = DZ(ATXX);
-            const double dz_At22 = DZ(ATYY);
-            const double dz_At33 = DZ(ATZZ);
-            const double dz_At13 = DZ(ATXZ);
-
-            const double dAt[3][3][3] = {
+
+            },
+            {
                 {
-                    { dx_At11,     0.0, dx_At13 },
-                    {     0.0, dx_At22,     0.0 },
-                    { dx_At13,     0.0, dx_At33 },
+                    { dxz_gt11,      0.0, dxz_gt13 },
+                    {      0.0, dxz_gt22,      0.0 },
+                    { dxz_gt13,      0.0, dxz_gt33 },
                 },
                 {
-                    {     0.0, on_axis ? dx_At11 - dx_At22 : (Atxx - Atyy) / x, 0.0 },
-                    { on_axis ? dx_At11 - dx_At22 : (Atxx - Atyy) / x, 0.0, on_axis ? dx_At13 : Atxz / x },
-                    { 0.0, on_axis ? dx_At13 : Atxz / x, 0.0 },
+                    { 0.0, on_axis ? dxz_gt11 - dxz_gt22 : (dz_gt11 - dz_gt22) / x, 0.0 },
+                    { on_axis ? dxz_gt11 - dxz_gt22 : (dz_gt11 - dz_gt22) / x, 0.0,
+                       on_axis ? dxz_gt13 : dz_gt13 / x },
+                    { 0.0, on_axis ? dxz_gt13 : dz_gt13 / x, 0.0 },
                 },
                 {
-                    { dz_At11,     0.0, dz_At13 },
-                    {     0.0, dz_At22,     0.0 },
-                    { dz_At13,     0.0, dz_At33 },
+                    { dzz_gt11,      0.0, dzz_gt13 },
+                    {      0.0, dzz_gt22,      0.0 },
+                    { dzz_gt13,      0.0, dzz_gt33 },
                 },
-            };
-
-            const double phi     = LOAD_VAR(PHI);
-            const double dphi[3] = { DX(PHI), 0.0, DZ(PHI) };
-
-            const double dxx_phi = DXX(PHI);
-            const double dzz_phi = DZZ(PHI);
-            const double dxz_phi = DXZ(PHI);
-
-            const double d2phi[3][3] = {
-                { dxx_phi,                             0.0, dxz_phi },
-                {     0.0, on_axis ? dxx_phi : dphi[0] / x,     0.0 },
-                { dxz_phi,                             0.0, dzz_phi },
-            };
-
-            const double trK     = LOAD_VAR(TRK);
-            const double dtrK[3] = { DX(TRK), 0.0, DZ(TRK) };
 
-            const double  Xtx  = LOAD_VAR(XTX);
-            const double  Xtz  = LOAD_VAR(XTZ);
-
-            const double alpha     = LOAD_VAR(ALPHA);
-            const double dx_alpha  = DX(ALPHA);
-            const double dz_alpha  = DZ(ALPHA);
-
-            const double dalpha[3] = { dx_alpha, 0.0, dz_alpha };
-
-            const double dxx_alpha  = DXX(ALPHA);
-            const double dzz_alpha  = DZZ(ALPHA);
-            const double dxz_alpha  = DXZ(ALPHA);
-
-            const double d2alpha[3][3] = {{ dxx_alpha,                                0, dxz_alpha },
-                                          {         0, on_axis ? dxx_alpha : dx_alpha / x,         0 },
-                                          { dxz_alpha,                                0, dzz_alpha }};
-
-            const double betax  = LOAD_VAR(BETAX);
-            const double betaz  = LOAD_VAR(BETAX);
-
-            const double beta[3] = { betax, 0.0, betaz };
-
-            const double dx_betax  = DX(BETAX);
-            const double dz_betax  = DZ(BETAX);
-
-            const double dx_betaz  = DX(BETAZ);
-            const double dz_betaz  = DZ(BETAZ);
-
-            const double dbeta[3][3] = {{ dx_betax,      0.0, dx_betaz },
-                                        {      0.0, on_axis ? dx_betax : betax / x,      0.0 },
-                                        { dz_betax,      0.0, dz_betaz }};
-
-            const double dxx_betax = DXX(BETAX);
-            const double dxz_betax = DXZ(BETAX);
-            const double dzz_betax = DZZ(BETAX);
-
-            const double dxx_betaz = DXX(BETAZ);
-            const double dxz_betaz = DXZ(BETAZ);
-            const double dzz_betaz = DZZ(BETAZ);
+            },
+        };
+
+        const double Atxx  = LOAD_VAR(ATXX);
+        const double Atyy  = LOAD_VAR(ATYY);
+        const double Atzz  = LOAD_VAR(ATZZ);
+        const double Atxy  = LOAD_VAR(ATXY);
+        const double Atxz  = LOAD_VAR(ATXZ);
+        const double Atyz  = LOAD_VAR(ATYZ);
+
+        const double At[3][3] = {
+            { Atxx, Atxy, Atxz },
+            { Atxy, Atyy, Atyz },
+            { Atxz, Atyz, Atzz }};
+
+        const double dx_At11 = DX(ATXX);
+        const double dx_At22 = DX(ATYY);
+        const double dx_At33 = DX(ATZZ);
+        const double dx_At13 = DX(ATXZ);
+
+        const double dz_At11 = DZ(ATXX);
+        const double dz_At22 = DZ(ATYY);
+        const double dz_At33 = DZ(ATZZ);
+        const double dz_At13 = DZ(ATXZ);
+
+        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, on_axis ? dx_At11 - dx_At22 : (Atxx - Atyy) / x, 0.0 },
+                { on_axis ? dx_At11 - dx_At22 : (Atxx - Atyy) / x, 0.0, on_axis ? dx_At13 : Atxz / x },
+                { 0.0, on_axis ? 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     = LOAD_VAR(PHI);
+        const double dphi[3] = { DX(PHI), 0.0, DZ(PHI) };
+
+        const double dxx_phi = DXX(PHI);
+        const double dzz_phi = DZZ(PHI);
+        const double dxz_phi = DXZ(PHI);
+
+        const double d2phi[3][3] = {
+            { dxx_phi,                             0.0, dxz_phi },
+            {     0.0, on_axis ? dxx_phi : dphi[0] / x,     0.0 },
+            { dxz_phi,                             0.0, dzz_phi },
+        };
+
+        const double trK     = LOAD_VAR(TRK);
+        const double dtrK[3] = { DX(TRK), 0.0, DZ(TRK) };
+
+        const double  Xtx  = LOAD_VAR(XTX);
+        const double  Xtz  = LOAD_VAR(XTZ);
+
+        const double alpha     = LOAD_VAR(ALPHA);
+        const double dx_alpha  = DX(ALPHA);
+        const double dz_alpha  = DZ(ALPHA);
+
+        const double dalpha[3] = { dx_alpha, 0.0, dz_alpha };
+
+        const double dxx_alpha  = DXX(ALPHA);
+        const double dzz_alpha  = DZZ(ALPHA);
+        const double dxz_alpha  = DXZ(ALPHA);
+
+        const double d2alpha[3][3] = {{ dxx_alpha,                                0, dxz_alpha },
+                                      {         0, on_axis ? dxx_alpha : dx_alpha / x,         0 },
+                                      { dxz_alpha,                                0, dzz_alpha }};
+
+        const double betax  = LOAD_VAR(BETAX);
+        const double betaz  = LOAD_VAR(BETAX);
+
+        const double beta[3] = { betax, 0.0, betaz };
+
+        const double dx_betax  = DX(BETAX);
+        const double dz_betax  = DZ(BETAX);
+
+        const double dx_betaz  = DX(BETAZ);
+        const double dz_betaz  = DZ(BETAZ);
+
+        const double dbeta[3][3] = {{ dx_betax,      0.0, dx_betaz },
+                                    {      0.0, on_axis ? dx_betax : betax / x,      0.0 },
+                                    { dz_betax,      0.0, dz_betaz }};
+
+        const double dxx_betax = DXX(BETAX);
+        const double dxz_betax = DXZ(BETAX);
+        const double dzz_betax = DZZ(BETAX);
+
+        const double dxx_betaz = DXX(BETAZ);
+        const double dxz_betaz = DXZ(BETAZ);
+        const double dzz_betaz = DZZ(BETAZ);
 
 #undef LOAD_VAR
 #undef DX
@@ -357,308 +357,308 @@ static void fill_eq_coeffs_line(const cGH *gh, const CoordPatch *cp, MG2DContext
 #undef DZZ
 #undef DXZ
 
-            const double d2beta[3][3][3] = {
-                {
-                    { dxx_betax, 0.0, dxx_betaz },
-                    { 0.0, on_axis ? 0.5 * dxx_betax : dx_betax / x - betax / SQR(x), 0.0 },
-                    { dxz_betax, 0.0, dxz_betaz },
-                },
-                {
-                    { 0.0, on_axis ? 0.5 * dxx_betax : dx_betax / x - betax / SQR(x), 0.0 },
-                    { on_axis ? 0.5 * dxx_betax : dx_betax / x - betax / SQR(x), 0.0, on_axis ? dxx_betaz : dx_betaz / x },
-                    { 0.0, on_axis ? dxz_betax : dz_betax / x, 0.0 },
-                },
-                {
-                    { dxz_betax, 0.0, dxz_betaz },
-                    { 0.0, on_axis ? dxz_betax : dz_betax / x, 0.0 },
-                    { dzz_betax, 0.0, dzz_betaz },
-                },
-            };
-
-            const double det = gtxx * gtyy * gtzz + 2 * gtxy * gtyz * gtxz - gtzz * SQR(gtxy) - SQR(gtxz) * gtyy - gtxx * SQR(gtyz);
-
-            // \tilde{γ}^{ij}
-            gtu[0][0] =  (gtyy * gtzz - SQR(gtyz)) / det;
-            gtu[1][1] =  (gtxx * gtzz - SQR(gtxz)) / det;
-            gtu[2][2] =  (gtxx * gtyy - SQR(gtxy)) / det;
-            gtu[0][1] = -(gtxy * gtzz - gtyz * gtxz) / det;
-            gtu[0][2] =  (gtxy * gtyz - gtyy * gtxz) / det;
-            gtu[1][2] = -(gtxx * gtyz - gtxy * gtxz) / det;
-            gtu[1][0] = gtu[0][1];
-            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;
+        const double d2beta[3][3][3] = {
+            {
+                { dxx_betax, 0.0, dxx_betaz },
+                { 0.0, on_axis ? 0.5 * dxx_betax : dx_betax / x - betax / SQR(x), 0.0 },
+                { dxz_betax, 0.0, dxz_betaz },
+            },
+            {
+                { 0.0, on_axis ? 0.5 * dxx_betax : dx_betax / x - betax / SQR(x), 0.0 },
+                { on_axis ? 0.5 * dxx_betax : dx_betax / x - betax / SQR(x), 0.0, on_axis ? dxx_betaz : dx_betaz / x },
+                { 0.0, on_axis ? dxz_betax : dz_betax / x, 0.0 },
+            },
+            {
+                { dxz_betax, 0.0, dxz_betaz },
+                { 0.0, on_axis ? dxz_betax : dz_betax / x, 0.0 },
+                { dzz_betax, 0.0, dzz_betaz },
+            },
+        };
+
+        const double det = gtxx * gtyy * gtzz + 2 * gtxy * gtyz * gtxz - gtzz * SQR(gtxy) - SQR(gtxz) * gtyy - gtxx * SQR(gtyz);
+
+        // \tilde{γ}^{ij}
+        gtu[0][0] =  (gtyy * gtzz - SQR(gtyz)) / det;
+        gtu[1][1] =  (gtxx * gtzz - SQR(gtxz)) / det;
+        gtu[2][2] =  (gtxx * gtyy - SQR(gtxy)) / det;
+        gtu[0][1] = -(gtxy * gtzz - gtyz * gtxz) / det;
+        gtu[0][2] =  (gtxy * gtyz - gtyy * gtxz) / det;
+        gtu[1][2] = -(gtxx * gtyz - gtxy * gtxz) / det;
+        gtu[1][0] = gtu[0][1];
+        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 γ_{kl}
-            // ∂_j K_{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
+        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++) {
-                        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}
+        // ∂_j K_{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]);
+                }
 
-            // Γ^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;
+        // ∂_{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 Γ^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 = γ^{kl} Γ_{kl}^j
-            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] * G[j][k][l];
-                X[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++) {
+                    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 β_k
-            for (int j = 0; j < 3; j++)
-                for (int k = 0; k < 3; k++) {
+        // Γ^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 l = 0; l < 3; l++)
-                        val += g[k][l] * dbeta[j][l] + dg[j][k][l] * beta[l];
-                    dbetal[j][k] = val;
+                    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;
                 }
 
-            // ∂_{jk} β_l
-            for (int j = 0; j < 3; j++)
-                for (int k = 0; k < 3; k++)
-                    for (int l = 0; l < 3; l++) {
+        // ∂_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 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;
+                        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;
                     }
 
-            // K_{ij}
-            for (int j = 0; j < 3; j++)
-                for (int k = 0; k < 3; k++)
-                    K[j][k] = At[j][k] / SQR(phi) + (1.0 / 3.0) * g[j][k] * trK;
+        // Γ^j = γ^{kl} Γ_{kl}^j
+        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] * G[j][k][l];
+            X[j] = val;
+        }
 
-            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] * K[l][k];
-                    Km[j][k] = 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;
+            }
 
-            // K^{ij}
-            for (int j = 0; j < 3; j++)
-                for (int k = 0; k < 3; k++) {
+        // ∂_{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 l = 0; l < 3; l++)
-                        val += gu[j][l] * Km[k][l];
-                    Ku[j][k] = val;
+                    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;
                 }
 
-            // ∂_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]);
-                    }
+        // K_{ij}
+        for (int j = 0; j < 3; j++)
+            for (int k = 0; k < 3; k++)
+                K[j][k] = At[j][k] / SQR(phi) + (1.0 / 3.0) * g[j][k] * trK;
 
-            // \dot{γ^{jk}}
-            for (int j = 0; j < 3; j++)
-                for (int k = 0; k < 3; k++) {
-                    gudot[j][k] = 2.0 * alpha * Ku[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] * K[l][k];
+                Km[j][k] = val;
+            }
 
-                    for (int l = 0; l < 3; l++) {
-                        gudot[j][k] += -gu[j][l] * dbeta[l][k] - gu[k][l] * dbeta[l][j];
+        // K^{ij}
+        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] * Km[k][l];
+                Ku[j][k] = val;
+            }
 
-                        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];
-                    }
+        // ∂_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++) {
+                gudot[j][k] = 2.0 * alpha * Ku[j][k];
 
-            // \dot{Γ}^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 += 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]);
+                for (int l = 0; l < 3; l++) {
+                    gudot[j][k] += -gu[j][l] * dbeta[l][k] - gu[k][l] * dbeta[l][j];
 
-                        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];
+                    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];
+                }
             }
 
-            for (int j = 0; j < 3; j++)
-                for (int k = 0; k < 3; k++) {
+
+        // \dot{Γ}^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 l = 0; l < 3; l++)
-                        val += G[l][j][k] * dalpha[l];
+                    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]);
 
-                    D2alpha[j][k] = d2alpha[j][k] - val;
+                    Gdot[j][k][l] = 0.5 * val;
                 }
 
-            Kij_Dij_alpha = 0.0;
-            for (int j = 0; j < 3; j++)
-                for (int k = 0; k < 3; k++) {
-                    Kij_Dij_alpha += Ku[j][k] * D2alpha[j][k];
+        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];
                 }
 
-            k3 = 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 += Km[k][l] * Km[l][j];
-                    k3 += val * Km[j][k];
-                }
+            Gammadot_term += val * dalpha[j];
+        }
+
+        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];
 
-            // K_{ij} K^{ij}
-            k2 = 0.0;
-            for (int j = 0; j < 3; j++)
-                for (int k = 0; k < 3; k++)
-                    k2 += Km[j][k] * Km[k][j];
+                D2alpha[j][k] = d2alpha[j][k] - val;
+            }
 
-            // 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;
-                }
+        Kij_Dij_alpha = 0.0;
+        for (int j = 0; j < 3; j++)
+            for (int k = 0; k < 3; k++) {
+                Kij_Dij_alpha += Ku[j][k] * D2alpha[j][k];
+            }
 
-            // 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;
-                }
+        k3 = 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 += Km[k][l] * Km[l][j];
+                k3 += val * Km[j][k];
+            }
 
-            for (int j = 0; j < 3; j++)
-                for (int k = 0; k < 3; k++) {
-                    double val = 0.0;
+        // K_{ij} K^{ij}
+        k2 = 0.0;
+        for (int j = 0; j < 3; j++)
+            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 += K[j][l] * Km[l][k];
+                        val += G[l][l][m] * G[m][j][k] - G[l][k][m] * G[m][j][l];
+                Ric[j][k] = val;
+            }
 
-                    Kdot[j][k] = -D2alpha[j][k] + alpha * (Ric[j][k] + trK * K[j][k] - 2 * 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;
+            }
 
-            k_kdot = 0.0;
-            for (int j = 0; j < 3; j++)
-                for (int k = 0; k < 3; k++) {
-                    //k_kdot += kdot[j][k] * Ku[j][k];
-                    k_kdot += Kdot[j][k] * Ku[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 += K[j][l] * Km[l][k];
+
+                Kdot[j][k] = -D2alpha[j][k] + alpha * (Ric[j][k] + trK * K[j][k] - 2 * val);
+            }
 
+        k_kdot = 0.0;
+        for (int j = 0; j < 3; j++)
+            for (int k = 0; k < 3; k++) {
+                //k_kdot += kdot[j][k] * Ku[j][k];
+                k_kdot += Kdot[j][k] * Ku[j][k];
+            }
 
-            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];
+        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;
 
-                        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];
+                    D2beta += d2beta[j][k][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];
-                        }
+                    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];
 
-                        beta_term += -gu[j][k] * D2beta * dalpha[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];
                     }
-            for (int j = 0; j < 3; j++)
-                for (int k = 0; k < 3; k++)
-                    beta_term += -beta[k] * Ricm[j][k] * dalpha[j];
-
-            solver->diff_coeffs[MG2D_DIFF_COEFF_20]->data[idx_dc] = gu[0][0] + (on_axis ? gu[1][1] : 0.0);
-            solver->diff_coeffs[MG2D_DIFF_COEFF_02]->data[idx_dc] = gu[2][2];
-            solver->diff_coeffs[MG2D_DIFF_COEFF_11]->data[idx_dc] = 2.0 * gu[0][2];
-            solver->diff_coeffs[MG2D_DIFF_COEFF_10]->data[idx_dc] = -X[0] + (on_axis ? 0.0 :  gu[1][1] / x);
-            solver->diff_coeffs[MG2D_DIFF_COEFF_01]->data[idx_dc] = -X[2];
-            solver->diff_coeffs[MG2D_DIFF_COEFF_00]->data[idx_dc] = -k2;
-
-            solver->rhs[idx_rhs] = -2.0 * alpha * Kij_Dij_alpha + Gammadot_term +
-                                    2.0 * alpha * (k_kdot + 2.0 * alpha * k3) + beta_term;
-
-            if (on_axis) {
-                solver->diff_coeffs[MG2D_DIFF_COEFF_20]->boundaries[MG2D_BOUNDARY_0L].val[idx_z] = gu[1][1];
-                solver->diff_coeffs[MG2D_DIFF_COEFF_10]->boundaries[MG2D_BOUNDARY_0L].val[idx_z] = gu[1][1];
-            }
+
+                    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];
+
+        solver->diff_coeffs[MG2D_DIFF_COEFF_20]->data[idx_dc] = gu[0][0] + (on_axis ? gu[1][1] : 0.0);
+        solver->diff_coeffs[MG2D_DIFF_COEFF_02]->data[idx_dc] = gu[2][2];
+        solver->diff_coeffs[MG2D_DIFF_COEFF_11]->data[idx_dc] = 2.0 * gu[0][2];
+        solver->diff_coeffs[MG2D_DIFF_COEFF_10]->data[idx_dc] = -X[0] + (on_axis ? 0.0 :  gu[1][1] / x);
+        solver->diff_coeffs[MG2D_DIFF_COEFF_01]->data[idx_dc] = -X[2];
+        solver->diff_coeffs[MG2D_DIFF_COEFF_00]->data[idx_dc] = -k2;
+
+        solver->rhs[idx_rhs] = -2.0 * alpha * Kij_Dij_alpha + Gammadot_term +
+                                2.0 * alpha * (k_kdot + 2.0 * alpha * k3) + beta_term;
+
+        if (on_axis) {
+            solver->diff_coeffs[MG2D_DIFF_COEFF_20]->boundaries[MG2D_BOUNDARY_0L].val[idx_z] = gu[1][1];
+            solver->diff_coeffs[MG2D_DIFF_COEFF_10]->boundaries[MG2D_BOUNDARY_0L].val[idx_z] = gu[1][1];
         }
+    }
 }