aboutsummaryrefslogtreecommitdiff
diff options
context:
space:
mode:
authorAnton Khirnov <anton@khirnov.net>2018-04-07 15:25:59 +0200
committerAnton Khirnov <anton@khirnov.net>2018-04-07 15:25:59 +0200
commit90381f1d94c437e58d233edeb43d0a1bb5c923f8 (patch)
tree2098d4cdc3bdd4fca7192a10c4bf068c9d2361f5
parent234e722c5abfb15c9850b4224cabbb3e2d0c3657 (diff)
Handle the upper solution branch.
-rw-r--r--init.c103
-rw-r--r--teukolsky_data.h2
-rw-r--r--teukolsky_data.py1
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):