From 90381f1d94c437e58d233edeb43d0a1bb5c923f8 Mon Sep 17 00:00:00 2001 From: Anton Khirnov Date: Sat, 7 Apr 2018 15:25:59 +0200 Subject: Handle the upper solution branch. --- init.c | 103 ++++++++++++++++++++++++++++++++++++++++++++++-------- teukolsky_data.h | 2 ++ teukolsky_data.py | 1 + 3 files changed, 92 insertions(+), 14 deletions(-) diff --git a/init.c b/init.c index 4fdece4..b1495a6 100644 --- a/init.c +++ b/init.c @@ -51,6 +51,7 @@ typedef struct TDPriv { TDLogger logger; double *coeffs; + double *coeffs_tmp; ptrdiff_t coeffs_stride; int cpu_flags; @@ -191,6 +192,14 @@ static int teukolsky_init_check_options(TDContext *td) return -ENOMEM; memset(s->coeffs, 0, sizeof(*s->coeffs) * NB_EQUATIONS * td->nb_coeffs[0] * td->nb_coeffs[1]); + if (td->solution_branch > 0) { + ret = posix_memalign((void**)&s->coeffs_tmp, 32, + sizeof(*s->coeffs_tmp) * NB_EQUATIONS * td->nb_coeffs[0] * td->nb_coeffs[1]); + if (ret) + return -ENOMEM; + memset(s->coeffs_tmp, 0, sizeof(*s->coeffs_tmp) * NB_EQUATIONS * td->nb_coeffs[0] * td->nb_coeffs[1]); + } + for (int i = 0; i < NB_EQUATIONS; i++) td->coeffs[i] = s->coeffs + i * td->nb_coeffs[0] * td->nb_coeffs[1]; @@ -245,17 +254,26 @@ static int teukolsky_constraint_eval(void *opaque, unsigned int eq_idx, } #define AMPLITUDE_DIRECT_SOLVE 1e-3 +#define AMPLITUDE_CONVERGE 0.00920788374206543 +#define AMPLITUDE_DIVERGE 0.00920788374267578 int td_solve(TDContext *td, double *coeffs_init[3]) { TDPriv *s = td->priv; - double cur_amplitude = td->amplitude;//MIN(td->amplitude, AMPLITUDE_DIRECT_SOLVE); int ret; ret = teukolsky_init_check_options(td); if (ret < 0) return ret; + if (td->amplitude > AMPLITUDE_CONVERGE) { + tdi_log(&s->logger, 0, + "Amplitude A=%16.16g is above the branch point, 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); + } + if (coeffs_init) { for (int i = 0; i < 3; i++) { memcpy(td->coeffs[i], coeffs_init[i], sizeof(*td->coeffs[i]) * @@ -263,28 +281,84 @@ int td_solve(TDContext *td, double *coeffs_init[3]) } } - while (1) { - double scale_factor; + if (td->solution_branch == 0 || coeffs_init) { + // direct solve with default (flat space) or user-provided initial guess + ret = tdi_nlsolve_solve(s->solver, teukolsky_constraint_eval, NULL, + &td->amplitude, s->coeffs); + if (ret < 0) { + tdi_log(&s->logger, 0, "tdi_nlsolve_solve() failed: %d", ret); + return ret; + } + } else { + // second branch requested and no user-provided initial guess + // 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; + + double cur_amplitude, new_amplitude, inverse_step; + int dir; ret = tdi_nlsolve_solve(s->solver, teukolsky_constraint_eval, NULL, - &cur_amplitude, s->coeffs); + &a0, s->coeffs); + if (ret < 0) + return ret; + ret = tdi_nlsolve_solve(s->solver, teukolsky_constraint_eval, NULL, + &a1, s->coeffs_tmp); if (ret < 0) return ret; - if (fabs(cur_amplitude - td->amplitude) < DBL_EPSILON) - goto finish; + cblas_daxpy(N, -1.0, s->coeffs, 1, s->coeffs_tmp, 1); + cblas_daxpy(N, -1.0, s->coeffs_tmp, 1, s->coeffs, 1); - scale_factor = cur_amplitude; - cur_amplitude = MIN(td->amplitude, 2 * cur_amplitude); - scale_factor = cur_amplitude / scale_factor; + // obtain solution for a1 in the upper branch + ret = tdi_nlsolve_solve(s->solver, teukolsky_constraint_eval, NULL, + &a1, s->coeffs); + if (ret < 0) { + tdi_log(&s->logger, 0, "Failed to get into the upper branch\n"); + return ret; + } + + cur_amplitude = a1; + dir = SGN(td->amplitude - cur_amplitude); + inverse_step = 1.0 / td->amplitude - 1.0 / cur_amplitude; + while (fabs(cur_amplitude - td->amplitude) > DBL_EPSILON) { + //double scale_factor; + //scale_factor = cur_amplitude; + //cur_amplitude = MIN(td->amplitude, 2 * cur_amplitude); + //scale_factor = cur_amplitude / scale_factor; + + new_amplitude = 1.0 / ((1.0 / cur_amplitude) + inverse_step); + if (dir * (td->amplitude - new_amplitude) < 0) + new_amplitude = td->amplitude; + tdi_log(&s->logger, 2, "Trying amplitude %g\n", new_amplitude); + 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); + if (ret == -EDOM) { + inverse_step = 0.5 * inverse_step; + if (fabs(inverse_step) < 1e-2) + return ret; + continue; + } else if (ret < 0) + return ret; + cur_amplitude = new_amplitude; + memcpy(s->coeffs, s->coeffs_tmp, sizeof(*s->coeffs) * N); + } + //while (1) { - cblas_dscal(td->nb_coeffs[0] * td->nb_coeffs[1], SQR(scale_factor), td->coeffs[0], 1); - cblas_dscal(td->nb_coeffs[0] * td->nb_coeffs[1], scale_factor, td->coeffs[1], 1); - cblas_dscal(td->nb_coeffs[0] * td->nb_coeffs[1], scale_factor, td->coeffs[2], 1); - tdi_log(&s->logger, 2, "Trying amplitude %g %g\n", cur_amplitude, scale_factor); - } + // if (fabs(cur_amplitude - td->amplitude) < DBL_EPSILON) + // goto finish; + + // cblas_dscal(td->nb_coeffs[0] * td->nb_coeffs[1], SQR(scale_factor), td->coeffs[0], 1); + // cblas_dscal(td->nb_coeffs[0] * td->nb_coeffs[1], scale_factor, td->coeffs[1], 1); + // cblas_dscal(td->nb_coeffs[0] * td->nb_coeffs[1], scale_factor, td->coeffs[2], 1); + //} + + } finish: tdi_solver_print_stats(s->solver); @@ -347,6 +421,7 @@ void td_context_free(TDContext **ptd) tdi_basis_free(&s->basis[i][j]); free(s->coeffs); + free(s->coeffs_tmp); free(td->priv); free(td); diff --git a/teukolsky_data.h b/teukolsky_data.h index 78bf824..1668bff 100644 --- a/teukolsky_data.h +++ b/teukolsky_data.h @@ -103,6 +103,8 @@ typedef struct TDContext { unsigned int nb_threads; double *coeffs[3]; + + unsigned int solution_branch; } TDContext; /** diff --git a/teukolsky_data.py b/teukolsky_data.py index 5bc9364..c9150de 100644 --- a/teukolsky_data.py +++ b/teukolsky_data.py @@ -37,6 +37,7 @@ class TeukolskyData(object): ("atol", ctypes.c_double), ("nb_threads", ctypes.c_uint), ("coeffs", ctypes.POINTER(ctypes.c_double) * 3), + ("solution_branch", ctypes.c_uint), ] def __init__(self, **kwargs): -- cgit v1.2.3