From 3910a4ee954866524bab129dfa4a0403f2ebe164 Mon Sep 17 00:00:00 2001 From: Anton Khirnov Date: Sun, 22 Apr 2018 12:58:22 +0200 Subject: Rewrite constraint evaluation code. Make it much more efficient and easier to add other seed functions. --- init.c | 90 ++++++++++++-- nlsolve.c | 4 +- nlsolve.h | 11 ++ td_constraints.c | 294 ++++++++++++++++++++++++++++++++++++++++------ td_constraints.h | 51 +++++++- td_constraints_template.c | 198 ------------------------------- 6 files changed, 397 insertions(+), 251 deletions(-) delete mode 100644 td_constraints_template.c diff --git a/init.c b/init.c index e90bbb8..d93d5a3 100644 --- a/init.c +++ b/init.c @@ -36,6 +36,7 @@ #include "cpu.h" #include "log.h" #include "nlsolve.h" +#include "td_constraints.h" #include "teukolsky_data.h" #include "threadpool.h" @@ -241,16 +242,47 @@ static int teukolsky_init_check_options(TDContext *td) return 0; } +static int constraint_eval_alloc(TDContext *td, double amplitude, + TDConstraintEvalContext **pce) +{ + TDPriv *priv = td->priv; + TDConstraintEvalContext *ce; + int ret; + + ret = tdi_constraint_eval_alloc(&ce, TD_FAMILY_AE_TIME_ANTISYM); + if (ret < 0) { + tdi_log(&priv->logger, 0, "Error allocating the constraints evaluator\n"); + return ret; + } + + ce->logger = priv->logger; + ce->amplitude = amplitude; + ce->nb_coords[0] = td->nb_coeffs[0]; + ce->nb_coords[1] = td->nb_coeffs[1]; + ce->coords[0] = priv->solver->colloc_grid[0][0]; + ce->coords[1] = priv->solver->colloc_grid[0][1]; + + ret = tdi_constraint_eval_init(ce); + if (ret < 0) { + tdi_constraint_eval_free(&ce); + tdi_log(&priv->logger, 0, "Error initializing the constraints evaluator\n"); + return ret; + } + + *pce = ce; + return 0; +} + static int teukolsky_constraint_eval(void *opaque, unsigned int eq_idx, const unsigned int *colloc_grid_order, const double * const *colloc_grid, const double * const (*vars)[PSSOLVE_DIFF_ORDER_NB], double *dst) { + TDConstraintEvalContext *ce = opaque; const double *amplitude = opaque; - return tdi_constraints_eval(eq_idx, *amplitude, colloc_grid_order, colloc_grid, - vars, dst); + return tdi_constraint_eval(ce, eq_idx, vars, dst); } #define AMPLITUDE_DIRECT_SOLVE 1e-3 @@ -260,19 +292,26 @@ static int teukolsky_constraint_eval(void *opaque, unsigned int eq_idx, int td_solve(TDContext *td, double *coeffs_init[3]) { TDPriv *s = td->priv; + TDConstraintEvalContext *ce; + double a0; int ret; ret = teukolsky_init_check_options(td); if (ret < 0) return ret; - if (td->amplitude > AMPLITUDE_CONVERGE) { + ret = constraint_eval_alloc(td, 0.0, &ce); + if (ret < 0) + return ret; + if (fabs(td->amplitude) >= ce->a_diverge) { tdi_log(&s->logger, 0, - "Amplitude A=%16.16g is above the branch point, no solutions " + "Amplitude A=%16.16g is above the point A_{max}=%g, no solutions " "are known to exist there. Set solution_branch=1 to get to the " "second branch where mass increases with decreasing amplitude\n", - td->amplitude); + td->amplitude, ce->a_converge); } + a0 = ce->a_converge; + tdi_constraint_eval_free(&ce); if (coeffs_init) { for (int i = 0; i < 3; i++) { @@ -283,8 +322,13 @@ int td_solve(TDContext *td, double *coeffs_init[3]) if (td->solution_branch == 0 || coeffs_init) { // direct solve with default (flat space) or user-provided initial guess + ret = constraint_eval_alloc(td, td->amplitude, &ce); + if (ret < 0) + return ret; + ret = tdi_nlsolve_solve(s->solver, teukolsky_constraint_eval, NULL, - &td->amplitude, s->coeffs, 0); + ce, s->coeffs, 0); + tdi_constraint_eval_free(&ce); if (ret < 0) { tdi_log(&s->logger, 0, "tdi_nlsolve_solve() failed: %d", ret); return ret; @@ -294,18 +338,27 @@ int td_solve(TDContext *td, double *coeffs_init[3]) // execute two lower-branch solutions and extrapolate to get to the upper branch const int N = NB_EQUATIONS * td->nb_coeffs[0] * td->nb_coeffs[1]; const double delta = 1e-5; - const double a0 = AMPLITUDE_CONVERGE; - const double a1 = AMPLITUDE_CONVERGE - delta; + const double a1 = a0 - delta; double cur_amplitude, new_amplitude, inverse_step; int dir; + ret = constraint_eval_alloc(td, a0, &ce); + if (ret < 0) + return ret; + ret = tdi_nlsolve_solve(s->solver, teukolsky_constraint_eval, NULL, - &a0, s->coeffs, 0); + ce, s->coeffs, 0); + tdi_constraint_eval_free(&ce); + if (ret < 0) + return ret; + + ret = constraint_eval_alloc(td, a1, &ce); if (ret < 0) return ret; ret = tdi_nlsolve_solve(s->solver, teukolsky_constraint_eval, NULL, - &a1, s->coeffs_tmp, 0); + ce, s->coeffs_tmp, 0); + tdi_constraint_eval_free(&ce); if (ret < 0) return ret; @@ -313,8 +366,13 @@ int td_solve(TDContext *td, double *coeffs_init[3]) cblas_daxpy(N, -1.0, s->coeffs_tmp, 1, s->coeffs, 1); // obtain solution for a1 in the upper branch + ret = constraint_eval_alloc(td, a1, &ce); + if (ret < 0) + return ret; + ret = tdi_nlsolve_solve(s->solver, teukolsky_constraint_eval, NULL, - &a1, s->coeffs, 0); + ce, s->coeffs, 0); + tdi_constraint_eval_free(&ce); if (ret < 0) { tdi_log(&s->logger, 0, "Failed to get into the upper branch\n"); return ret; @@ -333,9 +391,15 @@ int td_solve(TDContext *td, double *coeffs_init[3]) if (dir * (td->amplitude - new_amplitude) < 0) new_amplitude = td->amplitude; tdi_log(&s->logger, 2, "Trying amplitude %g\n", new_amplitude); + + ret = constraint_eval_alloc(td, new_amplitude, &ce); + if (ret < 0) + return ret; + memcpy(s->coeffs_tmp, s->coeffs, sizeof(*s->coeffs) * N); ret = tdi_nlsolve_solve(s->solver, teukolsky_constraint_eval, NULL, - &new_amplitude, s->coeffs_tmp, 1); + ce, s->coeffs_tmp, 1); + tdi_constraint_eval_free(&ce); if (ret == -EDOM) { inverse_step = 0.5 * inverse_step; if (fabs(inverse_step) < 1e-2) @@ -360,7 +424,7 @@ int td_solve(TDContext *td, double *coeffs_init[3]) } finish: - tdi_solver_print_stats(s->solver); + tdi_nlsolve_print_stats(s->solver); return 0; } diff --git a/nlsolve.c b/nlsolve.c index 708cd08..2137767 100644 --- a/nlsolve.c +++ b/nlsolve.c @@ -235,7 +235,7 @@ finish: return ret; } -void tdi_solver_print_stats(NLSolveContext *ctx) +void tdi_nlsolve_print_stats(NLSolveContext *ctx) { NLSolvePriv *s = ctx->priv; @@ -497,6 +497,8 @@ int tdi_nlsolve_context_init(NLSolveContext *ctx) return ret; } + ctx->colloc_grid = s->ps_ctx->colloc_grid; + /* init the per-equation state */ for (int i = 0; i < s->nb_equations; i++) { ret = eq_ctx_init(ctx, i); diff --git a/nlsolve.h b/nlsolve.h index f9b69fe..426cb24 100644 --- a/nlsolve.h +++ b/nlsolve.h @@ -72,6 +72,17 @@ typedef struct NLSolveContext { */ unsigned int (*solve_order)[2]; + /** + * Locations of the collocation points. The equation coefficients in NLEqCallback and + * NLEqJacobianCallback should be evaluated at those grid positions. + * + * colloc_grid[i][j] is an array of length solve_order[i][j] and contains + * the collocation points for the i-th equation in the j-th direction. + * + * Set by the solver after tdi_nlsolve_context_init(). + */ + double *(*colloc_grid)[2]; + #if HAVE_OPENCL cl_context ocl_ctx; cl_command_queue ocl_queue; diff --git a/td_constraints.c b/td_constraints.c index 6795891..fe6bc18 100644 --- a/td_constraints.c +++ b/td_constraints.c @@ -15,49 +15,275 @@ * along with this program. If not, see . */ +#include #include #include "common.h" #include "pssolve.h" #include "td_constraints.h" -#include "eval_constraints.c" -#include "eval_k_rtheta.c" - -#define EQNAME constraint_eq_ham -#define EQUATION TD_CONSTRAINT_EQ_HAM -#include "td_constraints_template.c" -#undef EQNAME -#undef EQUATION - -#define EQNAME constraint_eq_mom_0 -#define EQUATION TD_CONSTRAINT_EQ_MOM_0 -#include "td_constraints_template.c" -#undef EQNAME -#undef EQUATION - -#define EQNAME constraint_eq_mom_1 -#define EQUATION TD_CONSTRAINT_EQ_MOM_1 -#include "td_constraints_template.c" -#undef EQNAME -#undef EQUATION - -void -(*constraint_funcs[TD_CONSTRAINT_EQ_NB])(const double a, - unsigned int nb_coords[2], double *coords[2], +struct TDConstraintEvalPriv { + double *k_rtheta; + double *dk_rtheta[2]; +}; + +typedef struct TDFamilyDef { + double (*eval_krt) (TDConstraintEvalContext *ctx, double r, double theta); + double (*eval_krt_dr)(TDConstraintEvalContext *ctx, double r, double theta); + double (*eval_krt_dt)(TDConstraintEvalContext *ctx, double r, double theta); + double a_converge; + double a_diverge; +} TDFamilyDef; + +static double k_rtheta_eval_ae_time_antisym(TDConstraintEvalContext *ctx, + double r, double theta) +{ + const double r2 = SQR(r); + const double r3 = r * r2; + return -60.0 * sqrt(2) * (3.0 * r - 2.0 * r3) * exp(-r2) * sin(2.0 * theta) * ctx->amplitude / sqrt(M_PI); +} + +static double k_rtheta_eval_dr_ae_time_antisym(TDConstraintEvalContext *ctx, + double r, double theta) +{ + const double r2 = SQR(r); + const double r4 = SQR(r2); + return -60.0 * sqrt(2) * (3.0 - 12.0 * r2 + 4.0 * r4) * exp(-r2) * sin(2.0 * theta) * ctx->amplitude / sqrt(M_PI); +} + +static double k_rtheta_eval_dt_ae_time_antisym(TDConstraintEvalContext *ctx, + double r, double theta) +{ + const double r2 = SQR(r); + const double r3 = r * r2; + return -120.0 * sqrt(2) * (3.0 * r - 2.0 * r3) * exp(-r2) * cos(2.0 * theta) * ctx->amplitude / sqrt(M_PI); +} + +static const TDFamilyDef ae_time_antisym = { + .eval_krt = k_rtheta_eval_ae_time_antisym, + .eval_krt_dr = k_rtheta_eval_dr_ae_time_antisym, + .eval_krt_dt = k_rtheta_eval_dt_ae_time_antisym, + .a_converge = 0.0092078837420, + .a_diverge = 0.0092078837427, +}; + +#if 0 +static const TDFamilyDef simple_time_antisym = { + .eval_krt = k_rtheta_eval_simple_time_antisym, + .eval_krt_dr = k_rtheta_eval_dr_simple_time_antisym, + .eval_krt_dt = k_rtheta_eval_dt_simple_time_antisym, + .a_converge = 1.018918628476058, + .a_diverge = 1.018918628476968, +}; +#endif + +static const TDFamilyDef *td_families[] = { + [TD_FAMILY_AE_TIME_ANTISYM] = &ae_time_antisym, +}; + +static void +constraint_eval_ham_confflat(TDConstraintEvalContext *ctx, + double *vars[TD_CONSTRAINT_VAR_NB][PSSOLVE_DIFF_ORDER_NB], + double *out) +{ + TDConstraintEvalPriv *priv = ctx->priv; + + for (int idx1 = 0; idx1 < ctx->nb_coords[1]; idx1++) { + const double theta = ctx->coords[1][idx1]; + const double tt = tan(theta); + + for (int idx0 = 0; idx0 < ctx->nb_coords[0]; idx0++) { + const double r = ctx->coords[0][idx0]; + + const int idx = idx1 * ctx->nb_coords[0] + idx0; + + const double K_rt = priv->k_rtheta[idx]; + + const double K_rr = vars[TD_CONSTRAINT_VAR_CURV_RR][PSSOLVE_DIFF_ORDER_00][idx]; + const double K_pp = vars[TD_CONSTRAINT_VAR_CURV_PHIPHI][PSSOLVE_DIFF_ORDER_00][idx]; + const double psi = vars[TD_CONSTRAINT_VAR_CONFFAC][PSSOLVE_DIFF_ORDER_00][idx] + 1.0; + const double dpsi_r = vars[TD_CONSTRAINT_VAR_CONFFAC][PSSOLVE_DIFF_ORDER_10][idx]; + const double dpsi_t = vars[TD_CONSTRAINT_VAR_CONFFAC][PSSOLVE_DIFF_ORDER_01][idx]; + const double d2psi_rr = vars[TD_CONSTRAINT_VAR_CONFFAC][PSSOLVE_DIFF_ORDER_20][idx]; + const double d2psi_tt = vars[TD_CONSTRAINT_VAR_CONFFAC][PSSOLVE_DIFF_ORDER_02][idx]; + + const double Km[3][3] = {{ K_rr, K_rt, 0 }, + { K_rt / SQR(r), -(K_rr + K_pp), 0 }, + { 0, 0, K_pp }}; + + const double r2 = SQR(r); + const double psi2 = SQR(psi); + const double psi4 = SQR(psi2); + const double Rscal = -8.0 * (d2psi_rr + 2.0 * dpsi_r / r + d2psi_tt / r2 + dpsi_t / (r2 * tt)) / (psi4 * psi); + + 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; + } + } +} + +static void +constraint_eval_mom_r_confflat(TDConstraintEvalContext *ctx, + double *vars[TD_CONSTRAINT_VAR_NB][PSSOLVE_DIFF_ORDER_NB], + double *out) +{ + TDConstraintEvalPriv *priv = ctx->priv; + + for (int idx1 = 0; idx1 < ctx->nb_coords[1]; idx1++) { + const double theta = ctx->coords[1][idx1]; + const double tt = tan(theta); + + for (int idx0 = 0; idx0 < ctx->nb_coords[0]; idx0++) { + const double r = ctx->coords[0][idx0]; + + const int idx = idx1 * ctx->nb_coords[0] + idx0; + + const double K_rt = priv->k_rtheta[idx]; + const double dK_rt_t = priv->dk_rtheta[1][idx]; + + const double K_rr = vars[TD_CONSTRAINT_VAR_CURV_RR][PSSOLVE_DIFF_ORDER_00][idx]; + const double dK_rr_r = vars[TD_CONSTRAINT_VAR_CURV_RR][PSSOLVE_DIFF_ORDER_10][idx]; + const double psi = vars[TD_CONSTRAINT_VAR_CONFFAC][PSSOLVE_DIFF_ORDER_00][idx] + 1.0; + const double dpsi_r = vars[TD_CONSTRAINT_VAR_CONFFAC][PSSOLVE_DIFF_ORDER_10][idx]; + const double dpsi_t = vars[TD_CONSTRAINT_VAR_CONFFAC][PSSOLVE_DIFF_ORDER_01][idx]; + + const double r2 = SQR(r); + + out[idx] = 6.0 * (K_rr * dpsi_r + K_rt * dpsi_t / r2) / psi + + dK_rr_r + 3.0 * K_rr / r + K_rt / (r2 * tt) + dK_rt_t / r2; + } + } +} + +static void +constraint_eval_mom_t_confflat(TDConstraintEvalContext *ctx, + double *vars[TD_CONSTRAINT_VAR_NB][PSSOLVE_DIFF_ORDER_NB], + double *out) +{ + TDConstraintEvalPriv *priv = ctx->priv; + + for (int idx1 = 0; idx1 < ctx->nb_coords[1]; idx1++) { + const double theta = ctx->coords[1][idx1]; + const double tt = tan(theta); + + for (int idx0 = 0; idx0 < ctx->nb_coords[0]; idx0++) { + const double r = ctx->coords[0][idx0]; + + const int idx = idx1 * ctx->nb_coords[0] + idx0; + + const double K_rt = priv->k_rtheta[idx]; + const double dK_rt_r = priv->dk_rtheta[0][idx]; + + const double K_rr = vars[TD_CONSTRAINT_VAR_CURV_RR][PSSOLVE_DIFF_ORDER_00][idx]; + const double dK_rr_t = vars[TD_CONSTRAINT_VAR_CURV_RR][PSSOLVE_DIFF_ORDER_01][idx]; + const double K_pp = vars[TD_CONSTRAINT_VAR_CURV_PHIPHI][PSSOLVE_DIFF_ORDER_00][idx]; + const double dK_pp_t = vars[TD_CONSTRAINT_VAR_CURV_PHIPHI][PSSOLVE_DIFF_ORDER_01][idx]; + const double psi = vars[TD_CONSTRAINT_VAR_CONFFAC][PSSOLVE_DIFF_ORDER_00][idx] + 1.0; + const double dpsi_r = vars[TD_CONSTRAINT_VAR_CONFFAC][PSSOLVE_DIFF_ORDER_10][idx]; + const double dpsi_t = vars[TD_CONSTRAINT_VAR_CONFFAC][PSSOLVE_DIFF_ORDER_01][idx]; + + out[idx] = -(K_rr + 2.0 * K_pp) / tt - 6.0 * (K_rr + K_pp) * dpsi_t / psi + + 6.0 * K_rt * dpsi_r / psi - dK_pp_t - dK_rr_t + dK_rt_r + 2.0 * K_rt / r; + } + } +} + +int tdi_constraint_eval_init(TDConstraintEvalContext *ctx) +{ + TDConstraintEvalPriv *priv = ctx->priv; + const TDFamilyDef *fd = td_families[ctx->family]; + int ret; + + ret = posix_memalign((void**)&priv->k_rtheta, 32, + sizeof(*priv->k_rtheta) * ctx->nb_coords[0] * ctx->nb_coords[1]); + if (ret) + return -ENOMEM; + + for (int i = 0; i < ARRAY_ELEMS(priv->dk_rtheta); i++) { + ret = posix_memalign((void**)&priv->dk_rtheta[i], 32, + sizeof(*priv->dk_rtheta[i]) * ctx->nb_coords[0] * ctx->nb_coords[1]); + if (ret) + return -ENOMEM; + } + + ctx->a_converge = fd->a_converge; + ctx->a_diverge = fd->a_diverge; + + for (int idx1 = 0; idx1 < ctx->nb_coords[1]; idx1++) { + const double theta = ctx->coords[1][idx1]; + + for (int idx0 = 0; idx0 < ctx->nb_coords[0]; idx0++) { + const double r = ctx->coords[0][idx0]; + + const int idx = idx1 * ctx->nb_coords[0] + idx0; + + priv->k_rtheta[idx] = fd->eval_krt (ctx, r, theta); + priv->dk_rtheta[0][idx] = fd->eval_krt_dr(ctx, r, theta); + priv->dk_rtheta[1][idx] = fd->eval_krt_dt(ctx, r, theta); + } + } + + return 0; +} + +int tdi_constraint_eval_alloc(TDConstraintEvalContext **pctx, enum TDFamily family) +{ + TDConstraintEvalContext *ctx; + int ret; + + if (family < 0 || family >= ARRAY_ELEMS(td_families)) + return -EINVAL; + + ctx = calloc(1, sizeof(*ctx)); + if (!ctx) + return -ENOMEM; + + ctx->priv = calloc(1, sizeof(*ctx->priv)); + if (!ctx->priv) { + tdi_constraint_eval_free(&ctx); + return -ENOMEM; + } + + ctx->family = family; + + *pctx = ctx; + return 0; +} + +void tdi_constraint_eval_free(TDConstraintEvalContext **pctx) +{ + TDConstraintEvalContext *ctx = *pctx; + + if (!ctx) + return; + + if (ctx->priv) { + free(ctx->priv->k_rtheta); + free(ctx->priv->dk_rtheta[0]); + free(ctx->priv->dk_rtheta[1]); + } + free(ctx->priv); + + free(ctx); + *pctx = NULL; +} + +static const void +(*constraint_funcs[TD_CONSTRAINT_EQ_NB])(TDConstraintEvalContext *ctx, double *vars[TD_CONSTRAINT_VAR_NB][PSSOLVE_DIFF_ORDER_NB], double *out) = { - [TD_CONSTRAINT_EQ_HAM] = constraint_eq_ham, - [TD_CONSTRAINT_EQ_MOM_0] = constraint_eq_mom_0, - [TD_CONSTRAINT_EQ_MOM_1] = constraint_eq_mom_1, + [TD_CONSTRAINT_EQ_HAM] = constraint_eval_ham_confflat, + [TD_CONSTRAINT_EQ_MOM_0] = constraint_eval_mom_r_confflat, + [TD_CONSTRAINT_EQ_MOM_1] = constraint_eval_mom_t_confflat, }; -int tdi_constraints_eval(enum TDConstraintEq eq, - const double a, - unsigned int nb_coords[2], double *coords[2], - double *vars[TD_CONSTRAINT_VAR_NB][PSSOLVE_DIFF_ORDER_NB], - double *out) +int tdi_constraint_eval(TDConstraintEvalContext *ctx, enum TDConstraintEq eq, + const double *vars[TD_CONSTRAINT_VAR_NB][PSSOLVE_DIFF_ORDER_NB], + double *out) { - constraint_funcs[eq](a, nb_coords, coords, vars, out); - return 0; + constraint_funcs[eq](ctx, vars, out); } diff --git a/td_constraints.h b/td_constraints.h index 44f5257..f0a892c 100644 --- a/td_constraints.h +++ b/td_constraints.h @@ -20,6 +20,20 @@ #include "pssolve.h" +/** + * Identifies the specific initial data family to use. + */ +enum TDFamily { + /** + * The time-antisymmetric initial data used in Abrahams&Evans PhysRevD v49,n8 (1994). + * Conformally flat spatial metric. + * r / x x 3 \ / x 2 \ + * K = -60√(2/π) a | --- - (---) | * exp| - (---) | sin(2θ) + * θ \ L L / \ L / + */ + TD_FAMILY_AE_TIME_ANTISYM = 0, +}; + enum TDConstraintEq { TD_CONSTRAINT_EQ_HAM = 0, TD_CONSTRAINT_EQ_MOM_0, @@ -34,10 +48,37 @@ enum TDConstraintVar { TD_CONSTRAINT_VAR_NB, }; -int tdi_constraints_eval(enum TDConstraintEq eq, - const double a, - unsigned int nb_coords[2], double *coords[2], - double *vars[TD_CONSTRAINT_VAR_NB][PSSOLVE_DIFF_ORDER_NB], - double *out); + +typedef struct TDConstraintEvalPriv TDConstraintEvalPriv; + +typedef struct TDConstraintEvalContext { + /** + * Private data, not to be touched by the caller. + */ + TDConstraintEvalPriv *priv; + + /** + * The logging context. + * Set by the caller before tdi_constraint_eval_init(). + */ + TDLogger logger; + + unsigned int nb_coords[2]; + double *coords[2]; + + enum TDFamily family; + + double amplitude; + double a_converge; + double a_diverge; +} TDConstraintEvalContext; + +int tdi_constraint_eval_alloc(TDConstraintEvalContext **ctx, enum TDFamily family); +int tdi_constraint_eval_init (TDConstraintEvalContext *ctx); +void tdi_constraint_eval_free(TDConstraintEvalContext **ctx); + +int tdi_constraint_eval(TDConstraintEvalContext *ctx, enum TDConstraintEq eq, + const double *vars[TD_CONSTRAINT_VAR_NB][PSSOLVE_DIFF_ORDER_NB], + double *out); #endif // TEUKOLSKY_DATA_CONSTRAINTS_H 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 - * - * 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 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 - } - } -} -- cgit v1.2.3