/* * Minimal distortion -- template for the equations definitions * Copyright (C) 2016 Anton Khirnov * * 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 . */ #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; } } }