Initial commit.

1 files changed, 507 insertions, 0 deletions

diff --git a/src/gamma_freeze_template.c b/src/gamma_freeze_template.c
new file mode 100644
index 0000000..8edda4d
--- /dev/null
+++ b/src/gamma_freeze_template.c
@@ -0,0 +1,507 @@
+/*
+ * 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], dgtu[3][3][3], G[3][3][3], dG[3][3][3][3];
+        double Gt[3][3][3];
+        double dXt[3][3];
+        double A[3][3], Au[3][3], Atu[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 trK  = eq_ctx->interp_values[I_TRK][i];
+
+        const double dx_trK = eq_ctx->interp_values[I_TRK_DX][i];
+        const double dz_trK = eq_ctx->interp_values[I_TRK_DZ][i];
+
+        const double dtrK[3] = { dx_trK, 0.0, dz_trK };
+
+        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 Xt[3] = { Xtx, 0.0, Xtz };
+
+        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 \tilde{γ}^{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 += -gtu[k][m] * gtu[l][n] * dgt[j][m][n];
+                    dgtu[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++)
+                        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);
+                    }
+
+        // \tilde{Γ}^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 * gtu[j][m] * (dgt[k][l][m] + dgt[l][k][m] - dgt[m][k][l]);
+                    Gt[j][k][l] = val;
+                }
+
+        // Γ^j_{kl}
+        for (int j = 0; j < 3; j++)
+            for (int k = 0; k < 3; k++)
+                for (int l = 0; l < 3; l++) {
+                    double val = 0.0;
+                    for (int m = 0; m < 3; m++)
+                        val += 0.5 * gu[j][m] * (dg[k][l][m] + dg[l][k][m] - dg[m][k][l]);
+                    G[j][k][l] = val;
+                }
+
+        // ∂_j Γ^k_{lm}
+        for (int j = 0; j < 3; j++)
+            for (int k = 0; k < 3; k++)
+                for (int l = 0; l < 3; l++)
+                    for (int m = 0; m < 3; m++) {
+                        double val = 0.0;
+                        for (int n = 0; n < 3; n++) {
+                            val += dgu[j][k][n] * (dg    [l][m][n] +  dg   [m][l][n] -  dg   [n][l][m]) +
+                                    gu   [k][n] * (d2g[j][l][m][n] + d2g[j][m][l][n] - d2g[j][n][l][m]);
+                        }
+                        dG[j][k][l][m] = 0.5 * val;
+                    }
+
+        // ∂_j Γ^k
+        for (int j = 0; j < 3; j++)
+            for (int k = 0; k < 3; k++) {
+                double val = 0.0;
+                for (int l = 0; l < 3; l++)
+                    for (int m = 0; m < 3; m++)
+                        val += gtu[l][m] * dG[j][k][l][m] + dgtu[j][l][m] * G[k][l][m];
+                dXt[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;
+            }
+
+        // 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;
+                }
+
+        // \tilde{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 += gtu[j][l] * gtu[k][m] * At[l][m];
+                Atu[j][k] = 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] * Atu[0][j];
+            rhs_z += dalpha[j] * Atu[2][j];
+        }
+        double val_x = 0.0;
+        double val_z = 0.0;
+        for (int j = 0; j < 3; j++) {
+            for (int k = 0; k < 3; k++) {
+                val_x += -Gt[0][j][k] * Atu[j][k];
+                val_z += -Gt[2][j][k] * Atu[j][k];
+            }
+        }
+        rhs_x += val_x * alpha;
+        rhs_z += val_z * alpha;
+        for (int j = 0; j < 3; j++) {
+            rhs_x += alpha * (2.0 / 3.0) * gtu[0][j] * dtrK[j];
+            rhs_z += alpha * (2.0 / 3.0) * gtu[2][j] * dtrK[j];
+        }
+        for (int j = 0; j < 3; j++) {
+            rhs_x += alpha * 3.0 * Atu[0][j] * dphi[j]/ phi;
+            rhs_z += alpha * 3.0 * Atu[2][j] * dphi[j]/ phi;
+        }
+
+        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] = gtu[0][0] + c1o3 * gtu[0][0] + (zaxis ? 0.5 * (gtu[1][1] + c1o3 * gtu[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] = gtu[2][2];
+            /* ∂_{zz}β^z */
+            eq_ctx->eq_coeffs[1][PSSOLVE_DIFF_ORDER_02][i] = c1o3 * gtu[0][2];
+
+            /* ∂_{xz}β^x */
+            eq_ctx->eq_coeffs[0][PSSOLVE_DIFF_ORDER_11][i] = 2.0 * gtu[0][2] + c1o3 * gtu[0][2] + (zaxis ? c1o3 * gtu[0][2] : 0.0);
+            /* ∂_{xz}β^z */
+            eq_ctx->eq_coeffs[1][PSSOLVE_DIFF_ORDER_11][i] = c1o3 * gtu[0][0];
+
+            /* ∂_{x}β^x */
+            eq_ctx->eq_coeffs[0][PSSOLVE_DIFF_ORDER_10][i] = -Xt[0] + (2.0 / 3.0) * Xt[0] + (zaxis ? (2.0 / 3.0) * Xt[0] : (gtu[1][1] + c1o3 * gtu[0][0]) / x);
+            /* ∂_{x}β^z */
+            eq_ctx->eq_coeffs[1][PSSOLVE_DIFF_ORDER_10][i] = 0.0;
+
+            /* ∂_{z}β^x */
+            eq_ctx->eq_coeffs[0][PSSOLVE_DIFF_ORDER_01][i] = -Xt[2] + (zaxis ? 0.0 : c1o3 * gtu[0][2] / x);
+            /* ∂_{z}β^z */
+            eq_ctx->eq_coeffs[1][PSSOLVE_DIFF_ORDER_01][i] = (2.0 / 3.0) * Xt[0];
+
+            /* β^x */
+            eq_ctx->eq_coeffs[0][PSSOLVE_DIFF_ORDER_00][i] = dXt[0][0] + (zaxis ? 0.0 : (2.0 / 3.0) * Xt[0] / x - (gtu[1][1] + c1o3 * gtu[0][0]) / SQR(x));
+
+            /* β^z */
+            eq_ctx->eq_coeffs[1][PSSOLVE_DIFF_ORDER_00][i] = dXt[2][0];
+
+            eq_ctx->rhs[i]                                 = rhs_x;
+        } else {
+            /* eq 1 */
+            /* ∂_{xx}β^x */
+            eq_ctx->eq_coeffs[0][PSSOLVE_DIFF_ORDER_20][i] = c1o3 * gtu[2][0] + (zaxis ? 0.5 * c1o3 * gtu[2][0] : 0.0);
+            /* ∂_{xx}β^z */
+            eq_ctx->eq_coeffs[1][PSSOLVE_DIFF_ORDER_20][i] = gtu[0][0] + (zaxis ? gtu[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] = gtu[2][2] + c1o3 * gtu[2][2];
+            /* ∂_{xz}β^x */
+            eq_ctx->eq_coeffs[0][PSSOLVE_DIFF_ORDER_11][i] = c1o3 * gtu[2][2] + (zaxis ? c1o3 * gtu[2][2] : 0.0);
+            /* ∂_{xz}β^z */
+            eq_ctx->eq_coeffs[1][PSSOLVE_DIFF_ORDER_11][i] = 2.0 * gtu[0][2] + c1o3 * gtu[0][2];
+
+            /* ∂_{x}β^x */
+            eq_ctx->eq_coeffs[0][PSSOLVE_DIFF_ORDER_10][i] = (2.0 / 3.0) * Xt[2] + (zaxis ? (2.0 / 3.0) * Xt[2] : c1o3 * gtu[2][0] / x);
+            /* ∂_{x}β^z */
+            eq_ctx->eq_coeffs[1][PSSOLVE_DIFF_ORDER_10][i] = -Xt[0] + (zaxis ? 0.0 : gtu[1][1] / x);
+            /* ∂_{z}β^x */
+            eq_ctx->eq_coeffs[0][PSSOLVE_DIFF_ORDER_01][i] = (zaxis ? 0.0 : c1o3 * gtu[2][2] / x);
+            /* ∂_{z}β^z */
+            eq_ctx->eq_coeffs[1][PSSOLVE_DIFF_ORDER_01][i] = -Xt[2] + (2.0 / 3.0) * Xt[2];
+
+            /* β^x */
+            eq_ctx->eq_coeffs[0][PSSOLVE_DIFF_ORDER_00][i] = dXt[0][2] + (zaxis ? 0.0 : (2.0 / 3.0) * Xt[2] / x - c1o3 * gtu[2][0] / SQR(x));
+
+            /* β^z */
+            eq_ctx->eq_coeffs[1][PSSOLVE_DIFF_ORDER_00][i] = dXt[2][2];
+
+            eq_ctx->rhs[i] = rhs_z;
+        }
+    }
+}