path: root/src/md_solve_template.c

/*
 * Minimal distortion -- template for the equations definitions
 * Copyright (C) 2016 Anton Khirnov <anton@khirnov.net>
 *
 * This program is free software: you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation, either version 3 of the License, or
 * (at your option) any later version.
 *
 * This program is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with this program.  If not, see <http://www.gnu.org/licenses/>.
 */

#define FUNC3(a, b) a ## _ ## b
#define FUNC2(a, b) FUNC3(a, b)
#define FUNC(name) FUNC2(name, EQUATION)

/**
 * A template for calculating the equation coefficients.
 */
static void FUNC(calc_eq_coeffs)(void *arg,
                                 unsigned int job_idx,    unsigned int nb_jobs,
                                 unsigned int thread_idx, unsigned int nb_threads)
{
    const MDCalcEqThread   *et = arg;
    const MDSolver        *ctx = et->ctx;
    MDEquationContext  *eq_ctx = et->eq_ctx;

    const int start = job_idx * et->block_size;
    const int end   = MIN((job_idx + 1) * et->block_size, NB_COLLOC_POINTS(ctx));

    for (int i = start; i < end; i++) {
        const double x = eq_ctx->interp_coords[0][i];
        const double z = eq_ctx->interp_coords[2][i];
        const int zaxis = x <= EPS;

        double c1o3 = (1.0 / 3.0);

        double gtu[3][3], g[3][3], gu[3][3];
        double dg[3][3][3], d2g[3][3][3][3], dgu[3][3][3], G[3][3][3], dG[3][3][3][3];
        double A[3][3], Au[3][3];
        double dA[3][3][3], dAu[3][3][3];
        double Ric[3][3], Ricm[3][3];
        double rhs_x, rhs_z;

        const double gtxx = eq_ctx->interp_values[I_GTXX][i];
        const double gtyy = eq_ctx->interp_values[I_GTYY][i];
        const double gtzz = eq_ctx->interp_values[I_GTZZ][i];
        const double gtxy = eq_ctx->interp_values[I_GTXY][i];
        const double gtxz = eq_ctx->interp_values[I_GTXZ][i];
        const double gtyz = eq_ctx->interp_values[I_GTYZ][i];

        const double gt[3][3] = {{ gtxx, gtxy, gtxz },
                                 { gtxy, gtyy, gtyz },
                                 { gtxz, gtyz, gtzz }};

        const double dx_gt11 = eq_ctx->interp_values[I_GTXX_DX][i];
        const double dx_gt22 = eq_ctx->interp_values[I_GTYY_DX][i];
        const double dx_gt33 = eq_ctx->interp_values[I_GTZZ_DX][i];
        const double dx_gt13 = eq_ctx->interp_values[I_GTXZ_DX][i];

        const double dz_gt11 = eq_ctx->interp_values[I_GTXX_DZ][i];
        const double dz_gt22 = eq_ctx->interp_values[I_GTYY_DZ][i];
        const double dz_gt33 = eq_ctx->interp_values[I_GTZZ_DZ][i];
        const double dz_gt13 = eq_ctx->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 = eq_ctx->interp_values[I_GTXX_DXX][i];
        const double dxx_gt22 = eq_ctx->interp_values[I_GTYY_DXX][i];
        const double dxx_gt33 = eq_ctx->interp_values[I_GTZZ_DXX][i];
        const double dxx_gt13 = eq_ctx->interp_values[I_GTXZ_DXX][i];

        const double dxz_gt11 = eq_ctx->interp_values[I_GTXX_DXZ][i];
        const double dxz_gt22 = eq_ctx->interp_values[I_GTYY_DXZ][i];
        const double dxz_gt33 = eq_ctx->interp_values[I_GTZZ_DXZ][i];
        const double dxz_gt13 = eq_ctx->interp_values[I_GTXZ_DXZ][i];

        const double dzz_gt11 = eq_ctx->interp_values[I_GTXX_DZZ][i];
        const double dzz_gt22 = eq_ctx->interp_values[I_GTYY_DZZ][i];
        const double dzz_gt33 = eq_ctx->interp_values[I_GTZZ_DZZ][i];
        const double dzz_gt13 = eq_ctx->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 = eq_ctx->interp_values[I_ATXX][i];
        const double Atyy = eq_ctx->interp_values[I_ATYY][i];
        const double Atzz = eq_ctx->interp_values[I_ATZZ][i];
        const double Atxy = eq_ctx->interp_values[I_ATXY][i];
        const double Atxz = eq_ctx->interp_values[I_ATXZ][i];
        const double Atyz = eq_ctx->interp_values[I_ATYZ][i];

        const double dx_At11 = eq_ctx->interp_values[I_ATXX_DX][i];
        const double dx_At22 = eq_ctx->interp_values[I_ATYY_DX][i];
        const double dx_At33 = eq_ctx->interp_values[I_ATZZ_DX][i];
        const double dx_At13 = eq_ctx->interp_values[I_ATXZ_DX][i];

        const double dz_At11 = eq_ctx->interp_values[I_ATXX_DZ][i];
        const double dz_At22 = eq_ctx->interp_values[I_ATYY_DZ][i];
        const double dz_At33 = eq_ctx->interp_values[I_ATZZ_DZ][i];
        const double dz_At13 = eq_ctx->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    = eq_ctx->interp_values[I_PHI][i];

        const double phi_dx = eq_ctx->interp_values[I_PHI_DX][i];
        const double phi_dz = eq_ctx->interp_values[I_PHI_DZ][i];

        const double dphi[3] = { phi_dx, 0.0, phi_dz };

        const double phi_dxx = eq_ctx->interp_values[I_PHI_DXX][i];
        const double phi_dzz = eq_ctx->interp_values[I_PHI_DZZ][i];
        const double phi_dxz = eq_ctx->interp_values[I_PHI_DXZ][i];

        const double d2phi[3][3] = {
            { phi_dxx, 0.0, phi_dxz },
            {     0.0, zaxis ? phi_dxx : phi_dx / x, 0.0 },
            { phi_dxz,    0.0,  phi_dzz },
        };

        const double At[3][3] = {{ Atxx, Atxy, Atxz },
                                 { Atxy, Atyy, Atyz },
                                 { Atxz, Atyz, Atzz }};

        const double alpha     = eq_ctx->interp_values[I_ALPHA][i];
        const double dx_alpha  = eq_ctx->interp_values[I_ALPHA_DX][i];
        const double dz_alpha  = eq_ctx->interp_values[I_ALPHA_DZ][i];

        const double dalpha[3] = { dx_alpha, 0.0, dz_alpha };

        const double Xtx  = eq_ctx->interp_values[I_XTX][i];
        const double Xtz  = eq_ctx->interp_values[I_XTZ][i];

        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}
        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);
                    dA[j][k][l] = -2.0 * dphi[j] * At[k][l] / (phi * SQR(phi)) + dAt[j][k][l] / 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;
                }

        // ∂_{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++)
                        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;
                    }

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

        // A_{jk}
        for (int j = 0; j < 3; j++)
            for (int k = 0; k < 3; k++) {
                A[j][k] = At[j][k] / SQR(phi);
            }

        // d_j A^{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 += dgu[j][k][m] * gu[l][n] * A[m][n] + gu[k][m] * dgu[j][l][n] * A[m][n] + gu[k][m] * gu[l][n] * dA[j][m][n];
                    dAu[j][k][l] = val;
                }

        // A^{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++)
                    for (int m = 0; m < 3; m++)
                        val += gu[j][l] * gu[k][m] * A[l][m];
                Au[j][k] = val;
            }

        rhs_x = 0.0;
        rhs_z = 0.0;
        for (int j = 0; j < 3; j++) {
            rhs_x += dalpha[j] * Au[0][j];
            rhs_z += dalpha[j] * Au[2][j];
        }
        for (int j = 0; j < 3; j++) {
            rhs_x += alpha * dAu[j][0][j];
            rhs_z += alpha * dAu[j][2][j];
        }
        for (int j = 0; j < 3; j++) {
            double val_x = 0.0;
            double val_z = 0.0;
            for (int k = 0; k < 3; k++) {
                val_x += G[0][j][k] * Au[k][j];
                val_z += G[2][j][k] * Au[k][j];
            }
            rhs_x += val_x * alpha;
            rhs_z += val_z * alpha;
        }
        for (int j = 0; j < 3; j++) {
            double val_x = 0.0;
            double val_z = 0.0;
            for (int k = 0; k < 3; k++) {
                val_x += G[j][j][k] * Au[0][k];
                val_z += G[j][j][k] * Au[2][k];
            }
            rhs_x += val_x * alpha;
            rhs_z += val_z * alpha;
        }

        rhs_x *= 2.0;
        rhs_z *= 2.0;

        double X[3] = { 0.0 };
        for (int j = 0; j < 3; j++)
            for (int k = 0; k < 3; k++) {
                X[0] += gu[j][k] * G[0][j][k];
                X[2] += gu[j][k] * G[2][j][k];
            }

        if (EQUATION == 0) {
            /* eq 0 */
            /* ∂_{xx}β^x */
            eq_ctx->eq_coeffs[0][PSSOLVE_DIFF_ORDER_20][i] = gu[0][0] + c1o3 * gu[0][0] + (zaxis ? 0.5 * (gu[1][1] + c1o3 * gu[0][0]) : 0.0);
            /* ∂_{xx}β^z */
            eq_ctx->eq_coeffs[1][PSSOLVE_DIFF_ORDER_20][i] = 0.0;
            /* ∂_{zz}β^x */
            eq_ctx->eq_coeffs[0][PSSOLVE_DIFF_ORDER_02][i] = gu[2][2];
            /* ∂_{zz}β^z */
            eq_ctx->eq_coeffs[1][PSSOLVE_DIFF_ORDER_02][i] = c1o3 * gu[0][2];

            /* ∂_{xz}β^x */
            eq_ctx->eq_coeffs[0][PSSOLVE_DIFF_ORDER_11][i] = 2.0 * gu[0][2] + c1o3 * gu[0][2] + (zaxis ? c1o3 * gu[0][2] : 0.0);
            /* ∂_{xz}β^z */
            eq_ctx->eq_coeffs[1][PSSOLVE_DIFF_ORDER_11][i] = c1o3 * gu[0][0];

            /* ∂_{x}β^x */
            {
                double t0 = 0.0;
                double t1 = 0.0;
                for (int j = 0; j < 3; j++) {
                    t0 += gu[0][j] * G[0][j][0];
                    t1 += G[j][j][0];
                }
                eq_ctx->eq_coeffs[0][PSSOLVE_DIFF_ORDER_10][i] = 2.0 * t0 - X[0] + c1o3 * gu[0][0] * t1 + (zaxis ? 2.0 * gu[1][1] * G[0][1][1] : (gu[1][1] + c1o3 * gu[0][0]) / x);
            }
            /* ∂_{x}β^z */
            {
                double t0 = 0.0;
                double t1 = 0.0;
                for (int j = 0; j < 3; j++) {
                    t0 += gu[0][j] * G[0][j][2];
                    t1 += G[j][j][2];
                }
                eq_ctx->eq_coeffs[1][PSSOLVE_DIFF_ORDER_10][i] = 2.0 * t0 + c1o3 * gu[0][0] * t1;
            }

            /* ∂_{z}β^x */
            {
                double t0 = 0.0;
                double t1 = 0.0;
                for (int j = 0; j < 3; j++) {
                    t0 += gu[2][j] * G[0][j][0];
                    t1 += G[j][j][0];
                }
                eq_ctx->eq_coeffs[0][PSSOLVE_DIFF_ORDER_01][i] = 2.0 * t0 - X[2] + c1o3 * gu[0][2] * t1 + (zaxis ? 0.0 : c1o3 * gu[0][2] / x);
            }
            /* ∂_{z}β^z */
            {
                double t0 = 0.0;
                double t1 = 0.0;
                for (int j = 0; j < 3; j++) {
                    t0 += gu[2][j] * G[0][j][2];
                    t1 += G[j][j][2];
                }
                eq_ctx->eq_coeffs[1][PSSOLVE_DIFF_ORDER_01][i] = 2.0 * t0 + c1o3 * gu[0][2] * t1;
            }

            /* β^x */
            {
                double t0 = 0.0;
                double t1 = 0.0;
                for (int k = 0; k < 3; k++) {
                    for (int l = 0; l < 3; l++) {
                        double val = 0.0;
                        for (int j = 0; j < 3; j++)
                            val += G[0][k][j] * G[j][l][0] - G[j][k][l] * G[0][0][j];
                        t0 += gu[k][l] * (dG[k][0][l][0] + val);
                        t1 += gu[0][k] * dG[k][l][l][0];
                    }
                }
                eq_ctx->eq_coeffs[0][PSSOLVE_DIFF_ORDER_00][i] = t0 + c1o3 * t1 + Ricm[0][0] + (zaxis ? 0.0 : 2.0 * gu[1][1] * G[0][1][1] / x - (gu[1][1] + c1o3 * gu[0][0]) / SQR(x));
            }

            /* β^z */
            {
                double t0 = 0.0;
                double t1 = 0.0;
                for (int k = 0; k < 3; k++) {
                    for (int l = 0; l < 3; l++) {
                        double val = 0.0;
                        for (int j = 0; j < 3; j++)
                            val += G[0][k][j] * G[j][l][2] - G[j][k][l] * G[0][2][j];
                        t0 += gu[k][l] * (dG[k][0][l][2] + val);
                        t1 += gu[0][k] * dG[k][l][l][2];
                    }
                }
                eq_ctx->eq_coeffs[1][PSSOLVE_DIFF_ORDER_00][i] = t0 + c1o3 * t1 + Ricm[0][2];
            }

            eq_ctx->rhs[i]                                 = rhs_x;
        } else {
            /* eq 1 */
            /* ∂_{xx}β^x */
            eq_ctx->eq_coeffs[0][PSSOLVE_DIFF_ORDER_20][i] = c1o3 * gu[2][0] + (zaxis ? c1o3 * 0.5 * gu[2][0] : 0.0);
            /* ∂_{xx}β^z */
            eq_ctx->eq_coeffs[1][PSSOLVE_DIFF_ORDER_20][i] = gu[0][0] + (zaxis ? gu[1][1] : 0.0);
            /* ∂_{zz}β^x */
            eq_ctx->eq_coeffs[0][PSSOLVE_DIFF_ORDER_02][i] = 0.0;
            /* ∂_{zz}β^z */
            eq_ctx->eq_coeffs[1][PSSOLVE_DIFF_ORDER_02][i] = gu[2][2] + c1o3 * gu[2][2];
            /* ∂_{xz}β^x */
            eq_ctx->eq_coeffs[0][PSSOLVE_DIFF_ORDER_11][i] = c1o3 * gu[2][2] + (zaxis ? c1o3 * gu[2][2] : 0.0);
            /* ∂_{xz}β^z */
            eq_ctx->eq_coeffs[1][PSSOLVE_DIFF_ORDER_11][i] = 2.0 * gu[0][2] + c1o3 * gu[0][2];

            /* ∂_{x}β^x */
            {
                double t0 = 0.0;
                double t1 = 0.0;
                for (int j = 0; j < 3; j++) {
                    t0 += gu[0][j] * G[2][j][0];
                    t1 += G[j][j][0];
                }
                eq_ctx->eq_coeffs[0][PSSOLVE_DIFF_ORDER_10][i] = 2.0 * t0 + c1o3 * gu[2][0] * t1 + (zaxis ? 2.0 * gu[1][1] * G[2][1][1] : c1o3 * gu[2][0] / x);
            }
            /* ∂_{x}β^z */
            {
                double t0 = 0.0;
                double t1 = 0.0;
                for (int j = 0; j < 3; j++) {
                    t0 += gu[0][j] * G[2][j][2];
                    t1 += G[j][j][2];
                }
                eq_ctx->eq_coeffs[1][PSSOLVE_DIFF_ORDER_10][i] = 2.0 * t0 - X[0] + c1o3 * gu[2][0] * t1 + (zaxis ? 0.0 : gu[1][1] / x);
            }
            /* ∂_{z}β^x */
            {
                double t0 = 0.0;
                double t1 = 0.0;
                for (int j = 0; j < 3; j++) {
                    t0 += gu[2][j] * G[2][j][0];
                    t1 += G[j][j][0];
                }
                eq_ctx->eq_coeffs[0][PSSOLVE_DIFF_ORDER_01][i] = 2.0 * t0 + c1o3 * gu[2][2] * t1 + (zaxis ? 0.0 : c1o3 * gu[2][2] / x);
            }
            /* ∂_{z}β^z */
            {
                double t0 = 0.0;
                double t1 = 0.0;
                for (int j = 0; j < 3; j++) {
                    t0 += gu[2][j] * G[2][j][2];
                    t1 += G[j][j][2];
                }
                eq_ctx->eq_coeffs[1][PSSOLVE_DIFF_ORDER_01][i] = 2.0 * t0 - X[2] + c1o3 * gu[2][2] * t1;
            }

            /* β^x */
            {
                double t0 = 0.0;
                double t1 = 0.0;
                for (int k = 0; k < 3; k++) {
                    for (int l = 0; l < 3; l++) {
                        double val = 0.0;
                        for (int j = 0; j < 3; j++)
                            val += G[2][k][j] * G[j][l][0] - G[j][k][l] * G[2][0][j];
                        t0 += gu[k][l] * (dG[k][2][l][0] + val);
                        t1 += gu[2][k] * dG[k][l][l][0];
                    }
                }
                eq_ctx->eq_coeffs[0][PSSOLVE_DIFF_ORDER_00][i] = t0 + c1o3 * t1 + Ricm[2][0] + (zaxis ? 0.0 : 2.0 * gu[1][1] * G[2][1][1] / x - c1o3 * gu[2][0] / SQR(x));
            }

            /* β^z */
            {
                double t0 = 0.0;
                double t1 = 0.0;
                for (int k = 0; k < 3; k++) {
                    for (int l = 0; l < 3; l++) {
                        double val = 0.0;
                        for (int j = 0; j < 3; j++)
                            val += G[2][k][j] * G[j][l][2] - G[j][k][l] * G[2][2][j];
                        t0 += gu[k][l] * (dG[k][2][l][2] + val);
                        t1 += gu[2][k] * dG[k][l][l][2];
                    }
                }
                eq_ctx->eq_coeffs[1][PSSOLVE_DIFF_ORDER_00][i] = t0 + c1o3 * t1 + Ricm[2][2];
            }

            eq_ctx->rhs[i] = rhs_z;
        }
    }
}