aboutsummaryrefslogtreecommitdiff
path: root/td_constraints_template.c
diff options
context:
space:
mode:
Diffstat (limited to 'td_constraints_template.c')
-rw-r--r--td_constraints_template.c198
1 files changed, 0 insertions, 198 deletions
diff --git a/td_constraints_template.c b/td_constraints_template.c
deleted file mode 100644
index df09290..0000000
--- a/td_constraints_template.c
+++ /dev/null
@@ -1,198 +0,0 @@
-/*
- * Teukolsky data -- template for the constraint equations definitions
- * Copyright (C) 2017 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 L 1.0
-#define R_0 0.0
-
-static void EQNAME(double a,
- unsigned int nb_coords[2], double *coords[2],
- double *vars[TD_CONSTRAINT_VAR_NB][PSSOLVE_DIFF_ORDER_NB],
- double *out)
-{
- for (int idx1 = 0; idx1 < nb_coords[1]; idx1++) {
- const double coord1 = coords[1][idx1];
-
- for (int idx0 = 0; idx0 < nb_coords[0]; idx0++) {
- const double coord0 = coords[0][idx0];
-
- const double r = coord0;
- const double theta = coord1;
-
- const int idx = idx1 * nb_coords[0] + idx0;
-
- const double psi = vars[TD_CONSTRAINT_VAR_CONFFAC][PSSOLVE_DIFF_ORDER_00][idx] + 1.0;
- const double K_rr = vars[TD_CONSTRAINT_VAR_CURV_RR][PSSOLVE_DIFF_ORDER_00][idx];
- const double K_phiphi = vars[TD_CONSTRAINT_VAR_CURV_PHIPHI][PSSOLVE_DIFF_ORDER_00][idx];
-
- const double dpsi_10 = vars[TD_CONSTRAINT_VAR_CONFFAC][PSSOLVE_DIFF_ORDER_10][idx];
- const double dpsi_01 = vars[TD_CONSTRAINT_VAR_CONFFAC][PSSOLVE_DIFF_ORDER_01][idx];
- const double dpsi_11 = vars[TD_CONSTRAINT_VAR_CONFFAC][PSSOLVE_DIFF_ORDER_11][idx];
- const double dpsi_20 = vars[TD_CONSTRAINT_VAR_CONFFAC][PSSOLVE_DIFF_ORDER_20][idx];
- const double dpsi_02 = vars[TD_CONSTRAINT_VAR_CONFFAC][PSSOLVE_DIFF_ORDER_02][idx];
-
- const double dK_rr_10 = vars[TD_CONSTRAINT_VAR_CURV_RR][PSSOLVE_DIFF_ORDER_10][idx];
- const double dK_rr_01 = vars[TD_CONSTRAINT_VAR_CURV_RR][PSSOLVE_DIFF_ORDER_01][idx];
- const double dK_phiphi_10 = vars[TD_CONSTRAINT_VAR_CURV_PHIPHI][PSSOLVE_DIFF_ORDER_10][idx];
- const double dK_phiphi_01 = vars[TD_CONSTRAINT_VAR_CURV_PHIPHI][PSSOLVE_DIFF_ORDER_01][idx];
-
- //const double K_rtheta_val = K_rtheta(a, r, theta);
- //const double dK_rtheta_r_val = K_rtheta_dr(a, r, theta);
- //const double dK_rtheta_theta_val = K_rtheta_dtheta(a, r, theta);
- const double K_rtheta_val = eval_K_rtheta(a, L, R_0, coord0, coord1);
- const double dK_rtheta_r_val = eval_dK_rtheta_r(a, L, R_0, coord0, coord1);
- const double dK_rtheta_theta_val = eval_dK_rtheta_theta(a, L, R_0, coord0, coord1);
-
- const double Km[3][3] = {{ K_rr, K_rtheta_val, 0 },
- { K_rtheta_val / SQR(r), -(K_rr + K_phiphi), 0 },
- { 0, 0, K_phiphi }};
-
-#if 0
- double gu[3][3];
- double dgu[3][3][3], G[3][3][3], dG[3][3][3][3];
- double Ric[3][3];
-
- double g[3][3], dg[3][3][3], d2g[3][3][3][3];
-
- eval_metric (a, L, R_0, coord0, coord1, psi, dpsi_01, dpsi_02, dpsi_10, dpsi_11, dpsi_20, g);
- for (int i = 0; i < 3; i++)
- eval_dmetric[i] (a, L, R_0, coord0, coord1, psi, dpsi_01, dpsi_02, dpsi_10, dpsi_11, dpsi_20, dg[i]);
- for (int i = 0; i < 3; i++)
- for (int j = 0; j < 3; j++)
- eval_d2metric[i][j](a, L, R_0, coord0, coord1, psi, dpsi_01, dpsi_02, dpsi_10, dpsi_11, dpsi_20, d2g[i][j]);
-
- const double det = g[0][0] * g[1][1] * g[2][2] + 2 * g[0][1] * g[1][2] * g[0][2] - g[2][2] * SQR(g[0][1]) - SQR(g[0][2]) * g[1][1] - g[0][0] * SQR(g[1][2]);
-
- // γ^{ij}
- gu[0][0] = (g[1][1] * g[2][2] - SQR(g[1][2])) / det;
- gu[1][1] = (g[0][0] * g[2][2] - SQR(g[0][2])) / det;
- gu[2][2] = (g[0][0] * g[1][1] - SQR(g[0][1])) / det;
- gu[0][1] = -(g[0][1] * g[2][2] - g[1][2] * g[0][2]) / det;
- gu[0][2] = (g[0][1] * g[1][2] - g[1][1] * g[0][2]) / det;
- gu[1][2] = -(g[0][0] * g[1][2] - g[0][1] * g[0][2]) / det;
- gu[1][0] = gu[0][1];
- gu[2][0] = gu[0][2];
- gu[2][1] = gu[1][2];
-
- // ∂_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;
- }
-
- // 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;
- }
-
- double Rscal = 0.0;
- for (int i = 0; i < 3; i++)
- for (int j = 0; j < 3; j++)
- Rscal += gu[i][j] * Ric[i][j];
-
- const double dKm[3][3][3] = {
- {
- { dK_rr_10, dK_rtheta_r_val, 0 },
- { gu[1][1] / gu[0][0] * dK_rtheta_r_val + dgu[0][1][1] / gu[0][0] * K_rtheta_val - dgu[0][0][0] / SQR(gu[0][0]) * gu[1][1] * K_rtheta_val,
- -(dK_rr_10 + dK_phiphi_10), 0 },
- { 0, 0, dK_phiphi_10 },
- },
- {
- { dK_rr_01, dK_rtheta_theta_val, 0 },
- { gu[1][1] / gu[0][0] * dK_rtheta_theta_val + dgu[1][1][1] / gu[0][0] * K_rtheta_val - dgu[1][0][0] / SQR(gu[0][0]) * gu[1][1] * K_rtheta_val,
- -(dK_rr_01 + dK_phiphi_01), 0 },
- { 0, 0, dK_phiphi_01 },
- },
- { 0 },
- };
-#endif
-
- if (EQUATION == TD_CONSTRAINT_EQ_HAM) {
- double Rscal = eval_R_scalar(a, L, coord0, coord1, psi, dpsi_01, dpsi_02, dpsi_10, dpsi_11, dpsi_20,
- K_rr, dK_rr_01, dK_rr_10, K_phiphi, dK_phiphi_01, dK_phiphi_10);
- double K2 = 0.0;
- for (int i = 0; i < 3; i++)
- for (int j = 0; j < 3; j++)
- K2 += Km[i][j] * Km[j][i];
-
- out[idx] = Rscal - K2;
- } else if (EQUATION == TD_CONSTRAINT_EQ_MOM_0) {
- out[idx] = eval_Dcurv_m_r(a, L, coord0, coord1, psi, dpsi_01, dpsi_02, dpsi_10, dpsi_11, dpsi_20,
- K_rr, dK_rr_01, dK_rr_10, K_phiphi, dK_phiphi_01, dK_phiphi_10);
- } else {
- out[idx] = eval_Dcurv_m_t(a, L, coord0, coord1, psi, dpsi_01, dpsi_02, dpsi_10, dpsi_11, dpsi_20,
- K_rr, dK_rr_01, dK_rr_10, K_phiphi, dK_phiphi_01, dK_phiphi_10);
- }
-#if 0
- } else {
- int mom_comp = (EQUATION == TD_CONSTRAINT_EQ_MOM_0) ? 0 : 1;
-
- double dKm_sum = 0.0;
- for (int i = 0; i < 3; i++)
- dKm_sum += dKm[i][i][mom_comp];
-
- double Gterm1 = 0.0;
- for (int i = 0; i < 3; i++)
- for (int j = 0; j < 3; j++)
- Gterm1 += G[i][i][j] * Km[j][mom_comp];
-
- double Gterm2 = 0.0;
- for (int i = 0; i < 3; i++)
- for (int j = 0; j < 3; j++)
- Gterm2 += G[j][i][mom_comp] * Km[i][j];
-
- out[idx] = dKm_sum + Gterm1 - Gterm2;
- }
-#endif
- }
- }
-}
generated by cgit v1.2.3 (git 2.39.1) at 2024-05-07 14:03:38 +0000