path: root/src/fill_eq_coeffs_template.c

#define LOAD(arr, idx) (arr[idx])

static inline double fd1_4(const double *arr, const ptrdiff_t stride, const double h_inv)
{
    const double ap2 = LOAD(arr,  2 * stride);
    const double ap1 = LOAD(arr,  1 * stride);
    const double am1 = LOAD(arr, -1 * stride);
    const double am2 = LOAD(arr, -2 * stride);
    return (-ap2 + 8.0 * ap1 - 8.0 * am1 + am2) * h_inv / 12.0;
}

static inline double fd2_4(const double *arr, const ptrdiff_t stride, const double h_inv)
{
    const double ap2 = LOAD(arr,  2 * stride);
    const double ap1 = LOAD(arr,  1 * stride);
    const double a0  = LOAD(arr,  0 * stride);
    const double am1 = LOAD(arr, -1 * stride);
    const double am2 = LOAD(arr, -2 * stride);
    return (-ap2 + 16.0 * ap1 - 30.0 * a0 + 16.0 * am1 - am2) * SQR(h_inv) / 12.0;
}

static inline double fd11_4(const double *arr, ptrdiff_t stride_x, ptrdiff_t stride_z,
                            double dx_inv, double dz_inv)
{
    const double ap2p2 = LOAD(arr,  2 * stride_x + 2 * stride_z);
    const double ap2p1 = LOAD(arr,  2 * stride_x + 1 * stride_z);
    const double ap2m1 = LOAD(arr,  2 * stride_x - 1 * stride_z);
    const double ap2m2 = LOAD(arr,  2 * stride_x - 2 * stride_z);

    const double ap1p2 = LOAD(arr,  1 * stride_x + 2 * stride_z);
    const double ap1p1 = LOAD(arr,  1 * stride_x + 1 * stride_z);
    const double ap1m1 = LOAD(arr,  1 * stride_x - 1 * stride_z);
    const double ap1m2 = LOAD(arr,  1 * stride_x - 2 * stride_z);

    const double am1p2 = LOAD(arr, -1 * stride_x + 2 * stride_z);
    const double am1p1 = LOAD(arr, -1 * stride_x + 1 * stride_z);
    const double am1m1 = LOAD(arr, -1 * stride_x - 1 * stride_z);
    const double am1m2 = LOAD(arr, -1 * stride_x - 2 * stride_z);

    const double am2p2 = LOAD(arr, -2 * stride_x + 2 * stride_z);
    const double am2p1 = LOAD(arr, -2 * stride_x + 1 * stride_z);
    const double am2m1 = LOAD(arr, -2 * stride_x - 1 * stride_z);
    const double am2m2 = LOAD(arr, -2 * stride_x - 2 * stride_z);

    return (-1.0 * (-1.0 * ap2p2 + 8.0 * ap1p2 - 8.0 * am1p2 + 1.0 * am2p2)
            +8.0 * (-1.0 * ap2p1 + 8.0 * ap1p1 - 8.0 * am1p1 + 1.0 * am2p1)
            -8.0 * (-1.0 * ap2m1 + 8.0 * ap1m1 - 8.0 * am1m1 + 1.0 * am2m1)
            +1.0 * (-1.0 * ap2m2 + 8.0 * ap1m2 - 8.0 * am1m2 + 1.0 * am2m2)) * dx_inv * dz_inv / 144.0;
}

#undef LOAD

#if 0
#define FD8(arr, stride, h_inv) \
    ((-3.0 * arr[4 * stride] + 32.0 * arr[3 * stride] - 168.0 * arr[2 * stride] + 672.0 * arr[1 * stride] \
      - 672.0 * arr[-1 * stride] + 168.0 * arr[-2 * stride] - 32.0 * arr[-3 * stride] + 3.0 * arr[-4 * stride]) * h_inv / 840.0)
#define FD28(arr, stride, h_inv) \
    ((-9.0 * arr[ 4 * stride] + 128.0 * arr[ 3 * stride] - 1008.0 * arr[ 2 * stride] + 8064.0 * arr[ 1 * stride]    \
      -9.0 * arr[-4 * stride] + 128.0 * arr[-3 * stride] - 1008.0 * arr[-2 * stride] + 8064.0 * arr[-1 * stride]    \
      - 14350.0 * arr[0]) * SQR(h_inv) / 5040.0)

static double diff8_dx(const double *arr, ptrdiff_t stride, double h_inv)
{
    return FD8(arr, 1, h_inv);
}
static double diff8_dz(const double *arr, ptrdiff_t stride, double h_inv)
{
    return FD8(arr, stride, h_inv);
}

static double diff8_dxx(const double *arr, ptrdiff_t stride, double h_inv)
{
    return FD28(arr, 1, h_inv);
}
static double diff8_dzz(const double *arr, ptrdiff_t stride, double h_inv)
{
    return FD28(arr, stride, h_inv);
}
static double diff8_dxz(const double *arr, ptrdiff_t stride, double dx_inv, double dz_inv)
{
    static const double fd_coeffs[] = {  3.0, -32.0,   168.0, -672.0,      0.0,  672.0,  -168.0,  32.0, -3.0 };
    static const int stencil = ARRAY_ELEMS(fd_coeffs) / 2;
    double val = 0.0;
    for (int j = -stencil; j <= stencil; j++) {
        double tmp = 0.0;
        for (int i = -stencil; i <= stencil; i++)
            tmp += fd_coeffs[i + stencil] * arr[j * stride + i];
        val += fd_coeffs[j + stencil] * tmp;
    }
    return val * dx_inv * dz_inv / SQR(840.0);
}

#define diff1   fd1_8
#define diff2   fd2_8
#define diff11 fd11_8
#else
#define diff1   fd1_4
#define diff2   fd2_4
#define diff11 fd11_4
#endif

#define diff_dx(arr, stride, dx_inv, dz_inv)    (diff1 (arr,      1,         dx_inv))
#define diff_dz(arr, stride, dx_inv, dz_inv)    (diff1 (arr, stride,         dz_inv))
#define diff_dxx(arr, stride, dx_inv, dz_inv)   (diff2 (arr,      1,         dx_inv))
#define diff_dzz(arr, stride, dx_inv, dz_inv)   (diff2 (arr, stride,         dz_inv))
#define diff_dxz(arr, stride, dx_inv, dz_inv)   (diff11(arr,      1, stride, dx_inv, dz_inv))

static void fill_eq_coeffs_line(const cGH *gh, const CoordPatch *cp, MG2DContext *solver,
                                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];

#define LOAD_VAR(var)  (        vars[RHS_VAR_ ## var][idx_src])
#define DX(var)        (diff1 (&vars[RHS_VAR_ ## var][idx_src],        1,           dx_inv        ))
#define DZ(var)        (diff1 (&vars[RHS_VAR_ ## var][idx_src],           stride_z,         dz_inv))
#define DXX(var)       (diff2 (&vars[RHS_VAR_ ## var][idx_src],        1,           dx_inv))
#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] = {
                {
                    { 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] = {
                {
                    {
                        { 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 },
                    },
                    {
                        { 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 },
                    },

                },
                {
                    {
                        { 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 },
                    },

                },
            };

            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
#undef DZ
#undef DXX
#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;
            }

            // ∂_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 γ^{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 = γ^{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 β_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) + (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 += gu[j][l] * K[l][k];
                    Km[j][k] = val;
                }

            // 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;
                }

            // ∂_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];

                    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];
                    }
                }


            // \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]);

                        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 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];

                    D2alpha[j][k] = d2alpha[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];
                }

            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];
                }

            // 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 += 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;
                }

            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];

                        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];

            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];
            }
        }
}