From a877278dc70a0a5c0c172e163ba5b2e9d37e2870 Mon Sep 17 00:00:00 2001 From: Anton Khirnov Date: Thu, 3 May 2018 17:02:38 +0200 Subject: Implement evaluating maximal lapse for the initial data. --- init.c | 291 ++++++++++++++++++++++++++++++++++++++++++++++++++++- libteukolskydata.v | 2 +- teukolsky_data.h | 4 + teukolsky_data.py | 2 + 4 files changed, 294 insertions(+), 5 deletions(-) diff --git a/init.c b/init.c index 33e5587..e5ad3d8 100644 --- a/init.c +++ b/init.c @@ -53,6 +53,8 @@ typedef struct TDPriv { double *coeffs_tmp; ptrdiff_t coeffs_stride; + double *coeffs_lapse; + int cpu_flags; double (*scalarproduct_metric)(size_t len1, size_t len2, double *mat, @@ -500,12 +502,240 @@ void td_context_free(TDContext **ptd) free(s->coeffs); free(s->coeffs_tmp); + free(s->coeffs_lapse); free(td->priv); free(td); *ptd = NULL; } +typedef struct MaximalSlicingEval { + double *psi; + double *dpsi_r; + double *dpsi_t; + double *krr; + double *kpp; + double *krt; +} MaximalSlicingEval; + +static int maximal_slicing_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) +{ + MaximalSlicingEval *max = opaque; + + for (int idx1 = 0; idx1 < colloc_grid_order[1]; idx1++) { + const double theta = colloc_grid[1][idx1]; + const double st = sin(theta); + const double st2 = SQR(st); + const double ct = cos(theta); + + for (int idx0 = 0; idx0 < colloc_grid_order[0]; idx0++) { + const double r = colloc_grid[0][idx0]; + const double r2 = SQR(r); + + const int idx = idx1 * colloc_grid_order[0] + idx0; + + const double alpha = vars[0][PSSOLVE_DIFF_ORDER_00][idx] + 1.0; + const double dalpha_r = vars[0][PSSOLVE_DIFF_ORDER_10][idx]; + const double dalpha_t = vars[0][PSSOLVE_DIFF_ORDER_01][idx]; + const double d2alpha_rr = vars[0][PSSOLVE_DIFF_ORDER_20][idx]; + const double d2alpha_tt = vars[0][PSSOLVE_DIFF_ORDER_02][idx]; + const double d2alpha_rt = vars[0][PSSOLVE_DIFF_ORDER_11][idx]; + + const double d2alpha[3][3] = { + { d2alpha_rr, d2alpha_rt, 0.0 }, + { d2alpha_rt, d2alpha_tt, 0.0 }, + { 0.0, 0.0, 0.0 }, + }; + + const double psi = max->psi[idx]; + const double psi2 = SQR(psi); + const double psi3 = psi * psi2; + const double psi4 = SQR(psi2); + + const double dpsi_r = max->dpsi_r[idx]; + const double dpsi_t = max->dpsi_t[idx]; + + const double g[3][3] = { + { psi4, 0.0, 0.0 }, + { 0.0, psi4 * r2, 0.0 }, + { 0.0, 0.0, psi4 * r2 * st2}, + }; + const double gu[3][3] = { + { 1.0 / psi4, 0.0, 0.0 }, + { 0.0, 1.0 / (psi4 * r2), 0.0 }, + { 0.0, 0.0, 1.0 / (psi4 * r2 * st2) }, + }; + + const double dg[3][3][3] = { + { + { 4.0 * dpsi_r * psi3, 0.0, 0.0 }, + { 0.0, 4.0 * dpsi_r * psi3 * r2 + 2.0 * r * psi4, 0.0 }, + { 0.0, 0.0, 4.0 * dpsi_r * psi3 * r2 * st2 + 2.0 * r * psi4 * st2 }, + }, + { + { 4.0 * dpsi_t * psi3, 0.0, 0.0 }, + { 0.0, 4.0 * dpsi_t * psi3 * r2, 0.0 }, + { 0.0, 0.0, 4.0 * dpsi_t * psi3 * r2 * st2 + 2.0 * psi4 * r2 * st * ct }, + + }, + { + { 0.0, 0.0, 0.0 }, + { 0.0, 0.0, 0.0 }, + { 0.0, 0.0, 0.0 }, + }, + }; + + const double krr = max->krr[idx]; + const double kpp = max->kpp[idx]; + const double krt = max->krt[idx]; + + const double Km[3][3] = { + { krr, krt, 0.0 }, + { krt / r2, -(krr + kpp), 0.0 }, + { 0.0, 0.0, kpp }, + }; + + double G[3][3][3], X[3]; + double laplace_alpha, k2; + + for (int i = 0; i < 3; i++) + 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[i][l] * (dg[k][j][l] + dg[j][k][l] - dg[l][j][k]); + + G[i][j][k] = 0.5 * val; + } + + for (int i = 0; i < 3; i++) { + double val = 0.0; + for (int j = 0; j < 3; j++) + for (int k = 0; k < 3; k++) + val += gu[j][k] * G[i][j][k]; + X[i] = val; + } + + laplace_alpha = 0.0; + for (int i = 0; i < 3; i++) + for (int j = 0; j < 3; j++) + laplace_alpha += gu[i][j] * d2alpha[i][j]; + + laplace_alpha -= X[0] * dalpha_r + X[1] * dalpha_t; + + k2 = 0.0; + for (int i = 0; i < 3; i++) + for (int j = 0; j < 3; j++) + k2 += Km[i][j] * Km[j][i]; + + dst[idx] = laplace_alpha - alpha * k2; + } + } + + return 0; +} + +static int lapse_solve_max(const TDContext *td) +{ + MaximalSlicingEval max = { NULL }; + TDPriv *priv = td->priv; + NLSolveContext *nl; + double *coords_r = NULL; + double *coords_t = NULL; + double *coeffs; + int ret = 0; + + ret = tdi_nlsolve_context_alloc(&nl, 1); + + if (ret < 0) { + tdi_log(&priv->logger, 0, "Error allocating the non-linear solver\n"); + return ret; + } + + nl->logger = priv->logger; + nl->tp = priv->tp; + nl->maxiter = td->max_iter; + nl->atol = td->atol; + + nl->basis[0][0] = priv->basis[0][0]; + nl->basis[0][1] = priv->basis[0][1]; + nl->solve_order[0][0] = priv->basis_order[0][0]; + nl->solve_order[0][1] = priv->basis_order[0][1]; + +#if HAVE_OPENCL + nl->ocl_ctx = priv->ocl_ctx; + nl->ocl_queue = priv->ocl_queue; +#endif + + ret = tdi_nlsolve_context_init(nl); + if (ret < 0) { + tdi_log(&priv->logger, 0, "Error initializing the non-linear solver\n"); + goto fail; + } + + ret = posix_memalign((void**)&max.psi, 32, sizeof(*max.psi) * td->nb_coeffs[0] * td->nb_coeffs[1]); + ret |= posix_memalign((void**)&max.dpsi_r, 32, sizeof(*max.psi) * td->nb_coeffs[0] * td->nb_coeffs[1]); + ret |= posix_memalign((void**)&max.dpsi_t, 32, sizeof(*max.psi) * td->nb_coeffs[0] * td->nb_coeffs[1]); + ret |= posix_memalign((void**)&max.krr, 32, sizeof(*max.krr) * td->nb_coeffs[0] * td->nb_coeffs[1]); + ret |= posix_memalign((void**)&max.kpp, 32, sizeof(*max.kpp) * td->nb_coeffs[0] * td->nb_coeffs[1]); + ret |= posix_memalign((void**)&max.krt, 32, sizeof(*max.krt) * td->nb_coeffs[0] * td->nb_coeffs[1]); + ret |= posix_memalign((void**)&coords_r, 32, sizeof(*coords_r) * td->nb_coeffs[0] * td->nb_coeffs[1]); + ret |= posix_memalign((void**)&coords_t, 32, sizeof(*coords_t) * td->nb_coeffs[0] * td->nb_coeffs[1]); + ret |= posix_memalign((void**)&coeffs, 32, sizeof(*coeffs) * td->nb_coeffs[0] * td->nb_coeffs[1]); + if (ret) { + ret = -ENOMEM; + goto fail; + } + + for (int j = 0; j < td->nb_coeffs[1]; j++) { + for (int i = 0; i < td->nb_coeffs[0]; i++) { + coords_r[j * td->nb_coeffs[0] + i] = nl->colloc_grid[0][0][i]; + coords_t[j * td->nb_coeffs[0] + i] = nl->colloc_grid[0][1][j]; + } + } + + td_eval_psi(td, td->nb_coeffs[0] * td->nb_coeffs[1], coords_r, coords_t, + (const unsigned int [2]){ 0, 0 }, max.psi); + td_eval_psi(td, td->nb_coeffs[0] * td->nb_coeffs[1], coords_r, coords_t, + (const unsigned int [2]){ 1, 0 }, max.dpsi_r); + td_eval_psi(td, td->nb_coeffs[0] * td->nb_coeffs[1], coords_r, coords_t, + (const unsigned int [2]){ 0, 1 }, max.dpsi_t); + td_eval_krr(td, td->nb_coeffs[0] * td->nb_coeffs[1], coords_r, coords_t, + (const unsigned int [2]){ 0, 0 }, max.krr); + td_eval_kpp(td, td->nb_coeffs[0] * td->nb_coeffs[1], coords_r, coords_t, + (const unsigned int [2]){ 0, 0 }, max.kpp); + td_eval_krt(td, td->nb_coeffs[0] * td->nb_coeffs[1], coords_r, coords_t, + (const unsigned int [2]){ 0, 0 }, max.krt); + + ret = tdi_nlsolve_solve(nl, maximal_slicing_eval, NULL, &max, coeffs, 0); + if (ret < 0) { + tdi_log(&priv->logger, 0, "Error solving the maximal slicing equation\n"); + goto fail; + } + + priv->coeffs_lapse = coeffs; + coeffs = NULL; + +fail: + free(coords_r); + free(coords_t); + free(coeffs); + free(max.psi); + free(max.dpsi_r); + free(max.dpsi_t); + free(max.krr); + free(max.kpp); + free(max.krt); + + tdi_nlsolve_context_free(&nl); + return ret; +} + static double scalarproduct_metric_c(size_t len1, size_t len2, double *mat, double *vec1, double *vec2) { @@ -528,7 +758,7 @@ static int eval_var(const TDContext *td, unsigned int var_idx, double *basis_val[2] = { NULL }; int ret = 0; - if (diff_order[0] > 0 || diff_order[1] > 0) { + if (diff_order[0] > 2 || diff_order[1] > 2) { tdi_log(&s->logger, 0, "Derivatives of higher order than 2 are not supported.\n"); return -ENOSYS; } @@ -546,12 +776,12 @@ static int eval_var(const TDContext *td, unsigned int var_idx, double theta_val = theta[i]; double r_val = r[i]; - double val = (var_idx == 0) ? 1.0 : 0.0; + double val = (var_idx == 0 && diff_order[0] == 0 && diff_order[1] == 0) ? 1.0 : 0.0; for (int k = 0; k < td->nb_coeffs[0]; k++) - basis_val[0][k] = tdi_basis_eval(s->basis[var_idx][0], BS_EVAL_TYPE_VALUE, r_val, k); + basis_val[0][k] = tdi_basis_eval(s->basis[var_idx][0], diff_order[0], r_val, k); for (int k = 0; k < td->nb_coeffs[1]; k++) - basis_val[1][k] = tdi_basis_eval(s->basis[var_idx][1], BS_EVAL_TYPE_VALUE, theta_val, k); + basis_val[1][k] = tdi_basis_eval(s->basis[var_idx][1], diff_order[1], theta_val, k); val += scalarproduct_metric_c(td->nb_coeffs[0], td->nb_coeffs[1], td->coeffs[var_idx], basis_val[0], basis_val[1]); @@ -621,3 +851,56 @@ int td_eval_krt(const TDContext *td, return 0; } + +int td_eval_lapse(const TDContext *td, + size_t nb_coords, const double *r, const double *theta, + const unsigned int diff_order[2], + double *out) +{ + TDPriv *priv = td->priv; + double *basis_val[2] = { NULL }; + int ret; + + if (!priv->coeffs_lapse) { + ret = lapse_solve_max(td); + if (ret < 0) + return ret; + } + + if (diff_order[0] > 2 || diff_order[1] > 2) { + tdi_log(&priv->logger, 0, "Derivatives of higher order than 2 are not supported.\n"); + return -ENOSYS; + } + + posix_memalign((void**)&basis_val[0], 32, sizeof(*basis_val[0]) * td->nb_coeffs[0]); + posix_memalign((void**)&basis_val[1], 32, sizeof(*basis_val[1]) * td->nb_coeffs[1]); + if (!basis_val[0] || !basis_val[1]) { + ret = -ENOMEM; + goto fail; + } + memset(basis_val[0], 0, sizeof(*basis_val[0]) * td->nb_coeffs[0]); + memset(basis_val[1], 0, sizeof(*basis_val[1]) * td->nb_coeffs[1]); + + for (int i = 0; i < nb_coords; i++) { + double theta_val = theta[i]; + double r_val = r[i]; + + double val = (diff_order[0] == 0 && diff_order[1] == 0) ? 1.0 : 0.0; + + for (int k = 0; k < td->nb_coeffs[0]; k++) + basis_val[0][k] = tdi_basis_eval(priv->basis[0][0], diff_order[0], r_val, k); + for (int k = 0; k < td->nb_coeffs[1]; k++) + basis_val[1][k] = tdi_basis_eval(priv->basis[0][1], diff_order[1], theta_val, k); + + val += scalarproduct_metric_c(td->nb_coeffs[0], td->nb_coeffs[1], priv->coeffs_lapse, + basis_val[0], basis_val[1]); + + out[i] = val; + } + +fail: + free(basis_val[0]); + free(basis_val[1]); + + return 0; +} diff --git a/libteukolskydata.v b/libteukolskydata.v index 50bd496..2ecfdde 100644 --- a/libteukolskydata.v +++ b/libteukolskydata.v @@ -1,4 +1,4 @@ -LIBTEUKOLSKYDATA_1 { +LIBTEUKOLSKYDATA_2 { global: td_*; local: *; }; diff --git a/teukolsky_data.h b/teukolsky_data.h index d95dbfe..4967f8f 100644 --- a/teukolsky_data.h +++ b/teukolsky_data.h @@ -191,5 +191,9 @@ int td_eval_krt(const TDContext *td, size_t nb_coords, const double *r, const double *theta, const unsigned int diff_order[2], double *out); +int td_eval_lapse(const TDContext *td, + size_t nb_coords, const double *r, const double *theta, + const unsigned int diff_order[2], + double *out); #endif /* TEUKOLSKY_DATA_H */ diff --git a/teukolsky_data.py b/teukolsky_data.py index 2458bf7..86db668 100644 --- a/teukolsky_data.py +++ b/teukolsky_data.py @@ -124,6 +124,8 @@ class TeukolskyData(object): return self._eval_var(self._libtd.td_eval_kpp, r, theta, diff_order) def eval_krt(self, r, theta, diff_order = None): return self._eval_var(self._libtd.td_eval_krt, r, theta, diff_order) + def eval_lapse(self, r, theta, diff_order = None): + return self._eval_var(self._libtd.td_eval_lapse, r, theta, diff_order) @property -- cgit v1.2.3