From 234e722c5abfb15c9850b4224cabbb3e2d0c3657 Mon Sep 17 00:00:00 2001 From: Anton Khirnov Date: Fri, 9 Feb 2018 15:34:25 +0100 Subject: Add the remaining parts for teukolsky waves initial data --- Makefile | 2 + eq_gen/constraint_gen.py | 196 +++++++++++++++++++++ eval_constraints.c | 26 +++ eval_k_rtheta.c | 42 +++++ init.c | 438 ++++++++++++++++++++++++++++++++++++++++++++++ libteukolskydata.v | 4 + td_constraints.c | 195 +++++++++++++++++++++ td_constraints.h | 43 +++++ td_constraints_template.c | 198 +++++++++++++++++++++ teukolsky_data.h | 165 +++++++++++++++++ teukolsky_data.py | 128 ++++++++++++++ 11 files changed, 1437 insertions(+) create mode 100644 eq_gen/constraint_gen.py create mode 100644 eval_constraints.c create mode 100644 eval_k_rtheta.c create mode 100644 init.c create mode 100644 libteukolskydata.v create mode 100644 td_constraints.c create mode 100644 td_constraints.h create mode 100644 td_constraints_template.c create mode 100644 teukolsky_data.h create mode 100644 teukolsky_data.py diff --git a/Makefile b/Makefile index 14ad028..be6287a 100644 --- a/Makefile +++ b/Makefile @@ -8,9 +8,11 @@ OBJS = basis.o \ bicgstab.o \ cpu.o \ cpuid.o \ + init.o \ log.o \ nlsolve.o \ pssolve.o \ + td_constraints.o \ threadpool.o \ TESTPROGS = nlsolve \ diff --git a/eq_gen/constraint_gen.py b/eq_gen/constraint_gen.py new file mode 100644 index 0000000..b7033e8 --- /dev/null +++ b/eq_gen/constraint_gen.py @@ -0,0 +1,196 @@ +# +# Copyright 2018 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 . +# + +import sympy +import sympy.utilities.codegen as cg +import sys + +coords = sympy.symbols('x_0 x_1 x_2') + +psi_funcs = sympy.eye(3) +psi_vals = sympy.eye(3) +K_rr_funcs = sympy.eye(3) +K_rr_vals = sympy.eye(3) +K_phiphi_funcs = sympy.eye(3) +K_phiphi_vals = sympy.eye(3) +for diff0 in range(psi_funcs.shape[0]): + for diff1 in range(psi_funcs.shape[1]): + psi_funcs [diff0, diff1] = sympy.Function('psi_%d%d' % (diff0, diff1))(coords[0], coords[1]) + K_rr_funcs [diff0, diff1] = sympy.Function('K_rr_%d%d' % (diff0, diff1))(coords[0], coords[1]) + K_phiphi_funcs[diff0, diff1] = sympy.Function('K_phiphi_%d%d' % (diff0, diff1))(coords[0], coords[1]) + if diff0 + diff1 <= 2: + psi_vals [diff0, diff1] = sympy.symbols('psi_val_%d%d' % (diff0, diff1)) + else: + psi_vals[diff0, diff1] = None + + if diff0 + diff1 <= 1: + K_rr_vals [diff0, diff1] = sympy.symbols('K_rr_val_%d%d' % (diff0, diff1)) + K_phiphi_vals [diff0, diff1] = sympy.symbols('K_phiphi_val_%d%d' % (diff0, diff1)) + else: + K_rr_vals[diff0, diff1] = None + K_phiphi_vals[diff0, diff1] = None +psi = psi_funcs[0, 0] +K_rr = K_rr_funcs[0, 0] +K_phiphi = K_phiphi_funcs[0, 0] + +replace_diff = [] +replace_val = [] +for diff0 in range(psi_funcs.shape[0]): + for diff1 in range(psi_funcs.shape[1]): + replace_val.append((psi_funcs [diff0, diff1], psi_vals [diff0, diff1])) + replace_val.append((K_rr_funcs [diff0, diff1], K_rr_vals [diff0, diff1])) + replace_val.append((K_phiphi_funcs[diff0, diff1], K_phiphi_vals[diff0, diff1])) + + if diff0 + diff1 + 1 > 2: + continue + replace_diff.append((psi_funcs [diff0, diff1].diff(coords[0]), psi_funcs [diff0 + 1, diff1])) + replace_diff.append((psi_funcs [diff0, diff1].diff(coords[1]), psi_funcs [diff0, diff1 + 1])) + replace_diff.append((K_rr_funcs [diff0, diff1].diff(coords[0]), K_rr_funcs [diff0 + 1, diff1])) + replace_diff.append((K_rr_funcs [diff0, diff1].diff(coords[1]), K_rr_funcs [diff0, diff1 + 1])) + replace_diff.append((K_phiphi_funcs[diff0, diff1].diff(coords[0]), K_phiphi_funcs[diff0 + 1, diff1])) + replace_diff.append((K_phiphi_funcs[diff0, diff1].diff(coords[1]), K_phiphi_funcs[diff0, diff1 + 1])) + +r, theta, phi = coords + +x, a, L, r_0 = sympy.symbols('x a L r_0') +r_0 = 20 +kappa_p = sympy.sqrt(143 / sympy.pi) / 12288 +c_n = (((225 / (8 * sympy.pi)) ** 2)) ** (sympy.Rational(1, 4)) +Im2 = sympy.Piecewise((a * kappa_p * (L ** 5) * ((1 - ((x / L) ** 2)) ** 6), sympy.Abs(x) < L), + (0, True)) +#Im2 = a * c_n * (L ** 5) * sympy.exp(-(x ** 2) / (L ** 2)) +Im1 = Im2.diff(x) +I0 = Im1.diff(x) +Ip1 = I0.diff(x) +Ip2 = Ip1.diff(x) + +v = r# - r_0 +u = -r - r_0 + +Im2v = Im2.subs(x, v) +Im1v = Im1.subs(x, v) +I0v = I0.subs (x, v) +Ip1v = Ip1.subs(x, v) +Ip2v = Ip2.subs(x, v) + +Im2u = Im2.subs(x, u) +Im1u = Im1.subs(x, u) +I0u = I0.subs (x, u) +Ip1u = Ip1.subs(x, u) +Ip2u = Ip2.subs(x, u) + +Ip1v = sympy.Function('Ip1v')(v) +Ip2v = Ip1v.diff(v) +eta = ((Ip2v / r - 2 * Ip1v / (r ** 2)) * (sympy.sin(theta) ** 2)).simplify() +K_rtheta = ((Ip2v / (r ** 2) - 3 * Ip1v / (r ** 3) + 6 * I0v / (r ** 4) - 6 * Im1v / (r ** 5)) * r * sympy.sin(2 * theta)).simplify() + +#eta = (((Ip2u - Ip2v) / r + 2 * (Ip1u + Ip1v) / (r ** 2)) * (sympy.sin(theta) ** 2)).simplify() +#eta = 0 +#K_rtheta = (((Ip2u + Ip2v) / (r ** 2) + 3 * (Ip1u - Ip1v) / (r ** 3) + +# 6 * (I0u + I0v) / (r ** 4) + 6 * (Im1u - Im1v) / (r ** 5)) * r * sympy.sin(2 * theta)).simplify() + +eta = sympy.Function('eta')(r, theta) +e23eta = sympy.exp(2 * eta / 3) +metric = sympy.MutableDenseNDimArray([ + [e23eta, 0, 0], + [ 0, (r ** 2) * e23eta, 0], + [ 0, 0, (r * sympy.sin(theta)) ** 2 / (e23eta ** 2)], +]) * (psi ** 4) + +dmetric = sympy.derive_by_array(metric, coords) +d2metric = sympy.derive_by_array(dmetric, coords) + +metric_u = metric.applyfunc(lambda x: 0 if x == 0 else 1 / x) +dmetric_u = sympy.derive_by_array(metric_u, coords) + +Gamma = sympy.MutableDenseNDimArray.zeros(3, 3, 3) +for i in range(3): + for j in range(3): + for k in range(3): + val = sympy.Integer(0) + for l in range(3): + val += metric_u[i, l] * (dmetric[j, k, l] + dmetric[k, j, l] - dmetric[l, j, k]) + Gamma[i, j, k] = val / 2 + +dGamma = sympy.derive_by_array(Gamma, coords) + +Ricci = sympy.MutableDenseNDimArray.zeros(3, 3) +for j in range(3): + for k in range(3): + val = sympy.Integer(0) + for m in range(3): + val += dGamma[m, m, j, k] - dGamma[k, m, j, m] + for m in range(3): + for l in range(3): + val += Gamma[l, l, m] * Gamma[m, j, k] - Gamma[l, k, m] * Gamma[m, j, l] + Ricci[j, k] = val + +R_scalar = sympy.Integer(0) +for i in range(3): + for j in range(3): + R_scalar += metric_u[i, j] * Ricci[i, j] +sympy.pretty_print(R_scalar.simplify(), num_columns = 200) +sys.exit(0) + +curv_m = sympy.MutableDenseNDimArray([ + [ K_rr, K_rtheta, 0], + [ K_rtheta * metric_u[1, 1] / metric_u[0, 0], -(K_rr + K_phiphi), 0], + [ 0, 0, K_phiphi], +]) +dcurv_m = sympy.derive_by_array(curv_m, coords) + +Dcurv_m = sympy.MutableDenseNDimArray.zeros(3) +for i in range(3): + val = sympy.Integer(0) + for j in range(3): + val += dcurv_m[j, j, i] + for j in range(3): + for k in range(3): + val += Gamma[j, j, k] * curv_m[k, i] - Gamma[j, k, i] * curv_m[k, j] + Dcurv_m[i] = val + +def replace_derivatives(expr, replace): + while expr.find(sympy.Derivative): + expr = expr.subs(replace) + return expr + +R_scalar_val = replace_derivatives(R_scalar, replace_diff).subs(replace_val).simplify() +Dcurv_m_r_val = replace_derivatives(Dcurv_m[0], replace_diff).subs(replace_val).simplify() +Dcurv_m_t_val = replace_derivatives(Dcurv_m[1], replace_diff).subs(replace_val).simplify() + +funcs = [] + +res = sympy.MatrixSymbol('res', 3, 3) +funcs.append(('eval_R_scalar', R_scalar_val)) +funcs.append(('eval_Dcurv_m_r', Dcurv_m_r_val)) +funcs.append(('eval_Dcurv_m_t', Dcurv_m_t_val)) + +argument_sequence = ( + a, L, r, theta, + *filter(lambda x: x is not None, psi_vals), + *filter(lambda x: x is not None, K_rr_vals), + *filter(lambda x: x is not None, K_phiphi_vals), +) +#funcs = [ +# ('eval_K_rtheta', K_rtheta), +# ('eval_dK_rtheta_r', K_rtheta.diff(r)), +# ('eval_dK_rtheta_theta', K_rtheta.diff(theta)), +#] +#argument_sequence = (a, L, r, theta) +gen = cg.codegen(funcs, language = 'C', argument_sequence = argument_sequence) +sys.stdout.write(gen[1][1]) +sys.stdout.write(gen[0][1]) diff --git a/eval_constraints.c b/eval_constraints.c new file mode 100644 index 0000000..71d78bf --- /dev/null +++ b/eval_constraints.c @@ -0,0 +1,26 @@ +static inline +double eval_R_scalar(double a, double L, double x_0, double x_1, double psi_val_00, double psi_val_01, double psi_val_02, double psi_val_10, double psi_val_11, double psi_val_20, double K_rr_val_00, double K_rr_val_01, double K_rr_val_10, double K_phiphi_val_00, double K_phiphi_val_01, double K_phiphi_val_10) { + + double eval_R_scalar_result; + eval_R_scalar_result = -(8.0*psi_val_01/tan(x_1) + 8.0*psi_val_02 + 16.0*psi_val_10*x_0 + 8.0*psi_val_20*pow(x_0, 2))/(pow(psi_val_00, 5)*pow(x_0, 2)); + return eval_R_scalar_result; + +} + +static inline +double eval_Dcurv_m_r(double a, double L, double x_0, double x_1, double psi_val_00, double psi_val_01, double psi_val_02, double psi_val_10, double psi_val_11, double psi_val_20, double K_rr_val_00, double K_rr_val_01, double K_rr_val_10, double K_phiphi_val_00, double K_phiphi_val_01, double K_phiphi_val_10) { + + double eval_Dcurv_m_r_result; + eval_Dcurv_m_r_result = 3.0*K_rr_val_00/x_0 + 6.0*K_rr_val_00*psi_val_10/psi_val_00 + 1.0*K_rr_val_10 - 540.0*sqrt(2)*a*exp(-pow(x_0, 2)/pow(L, 2))*cos(2*x_1)/(sqrt(M_PI)*L*x_0) - 180.0*sqrt(2)*a*exp(-pow(x_0, 2)/pow(L, 2))/(sqrt(M_PI)*L*x_0) - 1080.0*sqrt(2)*a*psi_val_01*exp(-pow(x_0, 2)/pow(L, 2))*sin(2*x_1)/(sqrt(M_PI)*L*psi_val_00*x_0) + 360.0*sqrt(2)*a*x_0*exp(-pow(x_0, 2)/pow(L, 2))*cos(2*x_1)/(sqrt(M_PI)*pow(L, 3)) + 120.0*sqrt(2)*a*x_0*exp(-pow(x_0, 2)/pow(L, 2))/(sqrt(M_PI)*pow(L, 3)) + 720.0*sqrt(2)*a*psi_val_01*x_0*exp(-pow(x_0, 2)/pow(L, 2))*sin(2*x_1)/(sqrt(M_PI)*pow(L, 3)*psi_val_00); + return eval_Dcurv_m_r_result; + +} + +static inline +double eval_Dcurv_m_t(double a, double L, double x_0, double x_1, double psi_val_00, double psi_val_01, double psi_val_02, double psi_val_10, double psi_val_11, double psi_val_20, double K_rr_val_00, double K_rr_val_01, double K_rr_val_10, double K_phiphi_val_00, double K_phiphi_val_01, double K_phiphi_val_10) { + + double eval_Dcurv_m_t_result; + eval_Dcurv_m_t_result = -2.0*K_phiphi_val_00/tan(x_1) - 6.0*K_phiphi_val_00*psi_val_01/psi_val_00 - 1.0*K_phiphi_val_01 - 1.0*K_rr_val_00/tan(x_1) - 6.0*K_rr_val_00*psi_val_01/psi_val_00 - 1.0*K_rr_val_01 - 540.0*sqrt(2)*a*exp(-pow(x_0, 2)/pow(L, 2))*sin(2*x_1)/(sqrt(M_PI)*L) - 1080.0*sqrt(2)*a*psi_val_10*x_0*exp(-pow(x_0, 2)/pow(L, 2))*sin(2*x_1)/(sqrt(M_PI)*L*psi_val_00) + 960.0*sqrt(2)*a*pow(x_0, 2)*exp(-pow(x_0, 2)/pow(L, 2))*sin(2*x_1)/(sqrt(M_PI)*pow(L, 3)) + 720.0*sqrt(2)*a*psi_val_10*pow(x_0, 3)*exp(-pow(x_0, 2)/pow(L, 2))*sin(2*x_1)/(sqrt(M_PI)*pow(L, 3)*psi_val_00) - 240.0*sqrt(2)*a*pow(x_0, 4)*exp(-pow(x_0, 2)/pow(L, 2))*sin(2*x_1)/(sqrt(M_PI)*pow(L, 5)); + return eval_Dcurv_m_t_result; + +} diff --git a/eval_k_rtheta.c b/eval_k_rtheta.c new file mode 100644 index 0000000..0f55c9c --- /dev/null +++ b/eval_k_rtheta.c @@ -0,0 +1,42 @@ +/****************************************************************************** + * Code generated with sympy 1.1.1 * + * * + * See http://www.sympy.org/ for more information. * + * * + * This file is part of 'project' * + ******************************************************************************/ + + +#ifndef PROJECT__EVAL_K_RTHETA__H +#define PROJECT__EVAL_K_RTHETA__H + +double eval_K_rtheta(double a, double L, double r_0, double x_0, double x_1); +double eval_dK_rtheta_r(double a, double L, double r_0, double x_0, double x_1); +double eval_dK_rtheta_theta(double a, double L, double r_0, double x_0, double x_1); + +#endif + + +double eval_K_rtheta(double a, double L, double r_0, double x_0, double x_1) { + + double eval_K_rtheta_result; + eval_K_rtheta_result = 15*sqrt(2)*a*(-3*pow(L, 6)*r_0*exp(4*r_0*x_0/pow(L, 2)) + 3*pow(L, 6)*r_0 + 6*pow(L, 4)*pow(r_0, 2)*x_0*exp(4*r_0*x_0/pow(L, 2)) + 6*pow(L, 4)*pow(r_0, 2)*x_0 - 3*pow(L, 4)*r_0*pow(x_0, 2)*exp(4*r_0*x_0/pow(L, 2)) + 3*pow(L, 4)*r_0*pow(x_0, 2) - 6*pow(L, 2)*pow(r_0, 3)*pow(x_0, 2)*exp(4*r_0*x_0/pow(L, 2)) + 6*pow(L, 2)*pow(r_0, 3)*pow(x_0, 2) + 6*pow(L, 2)*pow(r_0, 2)*pow(x_0, 3)*exp(4*r_0*x_0/pow(L, 2)) + 6*pow(L, 2)*pow(r_0, 2)*pow(x_0, 3) + 6*pow(L, 2)*r_0*pow(x_0, 4)*exp(4*r_0*x_0/pow(L, 2)) - 6*pow(L, 2)*r_0*pow(x_0, 4) - 6*pow(L, 2)*pow(x_0, 5)*exp(4*r_0*x_0/pow(L, 2)) - 6*pow(L, 2)*pow(x_0, 5) + 4*pow(r_0, 4)*pow(x_0, 3)*exp(4*r_0*x_0/pow(L, 2)) + 4*pow(r_0, 4)*pow(x_0, 3) - 16*pow(r_0, 3)*pow(x_0, 4)*exp(4*r_0*x_0/pow(L, 2)) + 16*pow(r_0, 3)*pow(x_0, 4) + 24*pow(r_0, 2)*pow(x_0, 5)*exp(4*r_0*x_0/pow(L, 2)) + 24*pow(r_0, 2)*pow(x_0, 5) - 16*r_0*pow(x_0, 6)*exp(4*r_0*x_0/pow(L, 2)) + 16*r_0*pow(x_0, 6) + 4*pow(x_0, 7)*exp(4*r_0*x_0/pow(L, 2)) + 4*pow(x_0, 7))*exp(-(pow(r_0, 2) + 2*r_0*x_0 + pow(x_0, 2))/pow(L, 2))*sin(2*x_1)/(sqrt(M_PI)*pow(L, 3)*pow(x_0, 4)); + return eval_K_rtheta_result; + +} + +double eval_dK_rtheta_r(double a, double L, double r_0, double x_0, double x_1) { + + double eval_dK_rtheta_r_result; + eval_dK_rtheta_r_result = 15*sqrt(2)*a*(-6*pow(L, 4)*pow(r_0, 2)*exp(4*r_0*x_0/pow(L, 2)) + 6*pow(L, 4)*pow(r_0, 2) - 6*pow(L, 4)*r_0*x_0*exp(4*r_0*x_0/pow(L, 2)) + 6*pow(L, 4)*r_0*x_0 + 12*pow(L, 2)*pow(r_0, 3)*x_0*exp(4*r_0*x_0/pow(L, 2)) + 12*pow(L, 2)*pow(r_0, 3)*x_0 + 6*pow(L, 2)*pow(r_0, 2)*pow(x_0, 2)*exp(4*r_0*x_0/pow(L, 2)) + 18*pow(L, 2)*pow(r_0, 2)*pow(x_0, 2) + 24*pow(L, 2)*r_0*pow(x_0, 3)*exp(4*r_0*x_0/pow(L, 2)) - 24*pow(L, 2)*r_0*pow(x_0, 3) - 30*pow(L, 2)*pow(x_0, 4)*exp(4*r_0*x_0/pow(L, 2)) - 30*pow(L, 2)*pow(x_0, 4) - 12*pow(r_0, 4)*pow(x_0, 2)*exp(4*r_0*x_0/pow(L, 2)) + 12*pow(r_0, 4)*pow(x_0, 2) - 40*pow(r_0, 3)*pow(x_0, 3)*exp(4*r_0*x_0/pow(L, 2)) + 64*pow(r_0, 3)*pow(x_0, 3) + 144*pow(r_0, 2)*pow(x_0, 4)*exp(4*r_0*x_0/pow(L, 2)) + 120*pow(r_0, 2)*pow(x_0, 4) - 120*r_0*pow(x_0, 5)*exp(4*r_0*x_0/pow(L, 2)) + 96*r_0*pow(x_0, 5) + 28*pow(x_0, 6)*exp(4*r_0*x_0/pow(L, 2)) + 28*pow(x_0, 6) + 16*pow(r_0, 5)*pow(x_0, 3)*exp(4*r_0*x_0/pow(L, 2))/pow(L, 2) - 64*pow(r_0, 4)*pow(x_0, 4)*exp(4*r_0*x_0/pow(L, 2))/pow(L, 2) + 96*pow(r_0, 3)*pow(x_0, 5)*exp(4*r_0*x_0/pow(L, 2))/pow(L, 2) - 64*pow(r_0, 2)*pow(x_0, 6)*exp(4*r_0*x_0/pow(L, 2))/pow(L, 2) + 16*r_0*pow(x_0, 7)*exp(4*r_0*x_0/pow(L, 2))/pow(L, 2))*exp(-(pow(r_0, 2) + 2*r_0*x_0 + pow(x_0, 2))/pow(L, 2))*sin(2*x_1)/(sqrt(M_PI)*pow(L, 3)*pow(x_0, 4)) - 60*sqrt(2)*a*(-3*pow(L, 6)*r_0*exp(4*r_0*x_0/pow(L, 2)) + 3*pow(L, 6)*r_0 + 6*pow(L, 4)*pow(r_0, 2)*x_0*exp(4*r_0*x_0/pow(L, 2)) + 6*pow(L, 4)*pow(r_0, 2)*x_0 - 3*pow(L, 4)*r_0*pow(x_0, 2)*exp(4*r_0*x_0/pow(L, 2)) + 3*pow(L, 4)*r_0*pow(x_0, 2) - 6*pow(L, 2)*pow(r_0, 3)*pow(x_0, 2)*exp(4*r_0*x_0/pow(L, 2)) + 6*pow(L, 2)*pow(r_0, 3)*pow(x_0, 2) + 6*pow(L, 2)*pow(r_0, 2)*pow(x_0, 3)*exp(4*r_0*x_0/pow(L, 2)) + 6*pow(L, 2)*pow(r_0, 2)*pow(x_0, 3) + 6*pow(L, 2)*r_0*pow(x_0, 4)*exp(4*r_0*x_0/pow(L, 2)) - 6*pow(L, 2)*r_0*pow(x_0, 4) - 6*pow(L, 2)*pow(x_0, 5)*exp(4*r_0*x_0/pow(L, 2)) - 6*pow(L, 2)*pow(x_0, 5) + 4*pow(r_0, 4)*pow(x_0, 3)*exp(4*r_0*x_0/pow(L, 2)) + 4*pow(r_0, 4)*pow(x_0, 3) - 16*pow(r_0, 3)*pow(x_0, 4)*exp(4*r_0*x_0/pow(L, 2)) + 16*pow(r_0, 3)*pow(x_0, 4) + 24*pow(r_0, 2)*pow(x_0, 5)*exp(4*r_0*x_0/pow(L, 2)) + 24*pow(r_0, 2)*pow(x_0, 5) - 16*r_0*pow(x_0, 6)*exp(4*r_0*x_0/pow(L, 2)) + 16*r_0*pow(x_0, 6) + 4*pow(x_0, 7)*exp(4*r_0*x_0/pow(L, 2)) + 4*pow(x_0, 7))*exp(-(pow(r_0, 2) + 2*r_0*x_0 + pow(x_0, 2))/pow(L, 2))*sin(2*x_1)/(sqrt(M_PI)*pow(L, 3)*pow(x_0, 5)) - 15*sqrt(2)*a*(2*r_0 + 2*x_0)*(-3*pow(L, 6)*r_0*exp(4*r_0*x_0/pow(L, 2)) + 3*pow(L, 6)*r_0 + 6*pow(L, 4)*pow(r_0, 2)*x_0*exp(4*r_0*x_0/pow(L, 2)) + 6*pow(L, 4)*pow(r_0, 2)*x_0 - 3*pow(L, 4)*r_0*pow(x_0, 2)*exp(4*r_0*x_0/pow(L, 2)) + 3*pow(L, 4)*r_0*pow(x_0, 2) - 6*pow(L, 2)*pow(r_0, 3)*pow(x_0, 2)*exp(4*r_0*x_0/pow(L, 2)) + 6*pow(L, 2)*pow(r_0, 3)*pow(x_0, 2) + 6*pow(L, 2)*pow(r_0, 2)*pow(x_0, 3)*exp(4*r_0*x_0/pow(L, 2)) + 6*pow(L, 2)*pow(r_0, 2)*pow(x_0, 3) + 6*pow(L, 2)*r_0*pow(x_0, 4)*exp(4*r_0*x_0/pow(L, 2)) - 6*pow(L, 2)*r_0*pow(x_0, 4) - 6*pow(L, 2)*pow(x_0, 5)*exp(4*r_0*x_0/pow(L, 2)) - 6*pow(L, 2)*pow(x_0, 5) + 4*pow(r_0, 4)*pow(x_0, 3)*exp(4*r_0*x_0/pow(L, 2)) + 4*pow(r_0, 4)*pow(x_0, 3) - 16*pow(r_0, 3)*pow(x_0, 4)*exp(4*r_0*x_0/pow(L, 2)) + 16*pow(r_0, 3)*pow(x_0, 4) + 24*pow(r_0, 2)*pow(x_0, 5)*exp(4*r_0*x_0/pow(L, 2)) + 24*pow(r_0, 2)*pow(x_0, 5) - 16*r_0*pow(x_0, 6)*exp(4*r_0*x_0/pow(L, 2)) + 16*r_0*pow(x_0, 6) + 4*pow(x_0, 7)*exp(4*r_0*x_0/pow(L, 2)) + 4*pow(x_0, 7))*exp(-(pow(r_0, 2) + 2*r_0*x_0 + pow(x_0, 2))/pow(L, 2))*sin(2*x_1)/(sqrt(M_PI)*pow(L, 5)*pow(x_0, 4)); + return eval_dK_rtheta_r_result; + +} + +double eval_dK_rtheta_theta(double a, double L, double r_0, double x_0, double x_1) { + + double eval_dK_rtheta_theta_result; + eval_dK_rtheta_theta_result = 30*sqrt(2)*a*(-3*pow(L, 6)*r_0*exp(4*r_0*x_0/pow(L, 2)) + 3*pow(L, 6)*r_0 + 6*pow(L, 4)*pow(r_0, 2)*x_0*exp(4*r_0*x_0/pow(L, 2)) + 6*pow(L, 4)*pow(r_0, 2)*x_0 - 3*pow(L, 4)*r_0*pow(x_0, 2)*exp(4*r_0*x_0/pow(L, 2)) + 3*pow(L, 4)*r_0*pow(x_0, 2) - 6*pow(L, 2)*pow(r_0, 3)*pow(x_0, 2)*exp(4*r_0*x_0/pow(L, 2)) + 6*pow(L, 2)*pow(r_0, 3)*pow(x_0, 2) + 6*pow(L, 2)*pow(r_0, 2)*pow(x_0, 3)*exp(4*r_0*x_0/pow(L, 2)) + 6*pow(L, 2)*pow(r_0, 2)*pow(x_0, 3) + 6*pow(L, 2)*r_0*pow(x_0, 4)*exp(4*r_0*x_0/pow(L, 2)) - 6*pow(L, 2)*r_0*pow(x_0, 4) - 6*pow(L, 2)*pow(x_0, 5)*exp(4*r_0*x_0/pow(L, 2)) - 6*pow(L, 2)*pow(x_0, 5) + 4*pow(r_0, 4)*pow(x_0, 3)*exp(4*r_0*x_0/pow(L, 2)) + 4*pow(r_0, 4)*pow(x_0, 3) - 16*pow(r_0, 3)*pow(x_0, 4)*exp(4*r_0*x_0/pow(L, 2)) + 16*pow(r_0, 3)*pow(x_0, 4) + 24*pow(r_0, 2)*pow(x_0, 5)*exp(4*r_0*x_0/pow(L, 2)) + 24*pow(r_0, 2)*pow(x_0, 5) - 16*r_0*pow(x_0, 6)*exp(4*r_0*x_0/pow(L, 2)) + 16*r_0*pow(x_0, 6) + 4*pow(x_0, 7)*exp(4*r_0*x_0/pow(L, 2)) + 4*pow(x_0, 7))*exp(-(pow(r_0, 2) + 2*r_0*x_0 + pow(x_0, 2))/pow(L, 2))*cos(2*x_1)/(sqrt(M_PI)*pow(L, 3)*pow(x_0, 4)); + return eval_dK_rtheta_theta_result; + +} diff --git a/init.c b/init.c new file mode 100644 index 0000000..4fdece4 --- /dev/null +++ b/init.c @@ -0,0 +1,438 @@ +/* + * Copyright 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 . + */ + +#include "config.h" + +#include +#include +#include +#include +#include +#include + +#include + +#if HAVE_OPENCL +#include +#include +#endif + +#include "basis.h" +#include "common.h" +#include "cpu.h" +#include "log.h" +#include "nlsolve.h" +#include "teukolsky_data.h" +#include "threadpool.h" + +#define NB_EQUATIONS 3 + +typedef struct TDPriv { + unsigned int basis_order[NB_EQUATIONS][2]; + BasisSetContext *basis[NB_EQUATIONS][2]; + + NLSolveContext *solver; + + ThreadPoolContext *tp; + TDLogger logger; + + double *coeffs; + ptrdiff_t coeffs_stride; + + int cpu_flags; + + double (*scalarproduct_metric)(size_t len1, size_t len2, double *mat, + double *vec1, double *vec2); +} TDPriv; + +double tdi_scalarproduct_metric_fma3(size_t len1, size_t len2, double *mat, + double *vec1, double *vec2); +double tdi_scalarproduct_metric_avx(size_t len1, size_t len2, double *mat, + double *vec1, double *vec2); +double tdi_scalarproduct_metric_sse3(size_t len1, size_t len2, double *mat, + double *vec1, double *vec2); +double tdi_scalarproduct_metric_c(size_t len1, size_t len2, double *mat, + double *vec1, double *vec2); + +static void init_opencl(TDPriv *s) +#if HAVE_OPENCL +{ + int err, count; + cl_platform_id platform; + cl_context_properties props[3]; + cl_device_id ocl_device; + + err = clGetPlatformIDs(1, &platform, &count); + if (err != CL_SUCCESS || count < 1) { + tdi_log(&s->logger, 0, "Could not get an OpenCL platform ID\n"); + return; + } + + err = clGetDeviceIDs(platform, CL_DEVICE_TYPE_GPU, 1, &ocl_device, &count); + if (err != CL_SUCCESS || count < 1) { + tdi_log(&s->logger, 0, "Could not get an OpenCL device ID\n"); + return; + } + + props[0] = CL_CONTEXT_PLATFORM; + props[1] = (cl_context_properties)platform; + props[2] = 0; + + s->ocl_ctx = clCreateContext(props, 1, &ocl_device, NULL, NULL, &err); + if (err != CL_SUCCESS || !s->ocl_ctx) { + tdi_log(&s->logger, 0, "Could not create an OpenCL context\n"); + return; + } + + s->ocl_queue = clCreateCommandQueue(s->ocl_ctx, ocl_device, 0, &err); + if (err != CL_SUCCESS || !s->ocl_queue) { + tdi_log(&s->logger, 0, "Could not create an OpenCL command queue: %d\n", err); + goto fail; + } + + err = clblasSetup(); + if (err != CL_SUCCESS) { + tdi_log(&s->logger, 0, "Error setting up clBLAS\n"); + goto fail; + } + + return; +fail: + if (s->ocl_queue) + clReleaseCommandQueue(s->ocl_queue); + s->ocl_queue = 0; + + if (s->ocl_ctx) + clReleaseContext(s->ocl_ctx); + s->ocl_ctx = 0; +} +#else +{ +} +#endif + +static const enum BasisFamily basis_sets[NB_EQUATIONS][2] = { + { BASIS_FAMILY_SB_EVEN, BASIS_FAMILY_COS_EVEN }, + { BASIS_FAMILY_SB_EVEN, BASIS_FAMILY_COS_EVEN }, + { BASIS_FAMILY_SB_EVEN, BASIS_FAMILY_COS_EVEN }, +}; + +static int teukolsky_init_check_options(TDContext *td) +{ + TDPriv *s = td->priv; + int ret; + + s->cpu_flags = tdi_init_cpu_flags(); + + if (!td->nb_threads) { + td->nb_threads = tdi_cpu_count(); + if (!td->nb_threads) + td->nb_threads = 1; + } + + ret = tdi_threadpool_init(&s->tp, td->nb_threads); + if (ret < 0) + return ret; + + init_opencl(s); + + s->logger.log = tdi_log_default_callback; + + //s->scalarproduct_metric = tdi_scalarproduct_metric_c; + //if (EXTERNAL_SSE3(s->cpu_flags)) + // s->scalarproduct_metric = tdi_scalarproduct_metric_sse3; + //if (EXTERNAL_AVX(s->cpu_flags)) + // s->scalarproduct_metric = tdi_scalarproduct_metric_avx; + //if (EXTERNAL_FMA3(s->cpu_flags)) + // s->scalarproduct_metric = tdi_scalarproduct_metric_fma3; +// for (int i = 0; i < ctx->nb_equations; i++) +// for (int j = 0; j < 2; j++) { +// s->ps_ctx->basis[i][j] = ctx->basis[i][j]; +// s->ps_ctx->solve_order[i][j] = basis_order[i][j]; +// max_order = MAX(max_order, basis_order[i][j]); +// } +// + //for (int i = 0; i < max_order; i++) { + // tdi_log(&s->logger, 2, "%d ", i); + // for (int j = 0; j < ctx->nb_equations; j++) + // for (int k = 0; k < 2; k++) { + // if (i < s->ps_ctx->solve_order[j][k]) + // tdi_log(&s->logger, 2, "%8.8g\t", s->ps_ctx->colloc_grid[j][k][i]); + // else + // tdi_log(&s->logger, 2, " "); + // } + // tdi_log(&s->logger, 2, "\n"); + //} + + s->basis_order[0][0] = td->nb_coeffs[0]; + s->basis_order[0][1] = td->nb_coeffs[1]; + s->basis_order[1][0] = td->nb_coeffs[0]; + s->basis_order[1][1] = td->nb_coeffs[1]; + s->basis_order[2][0] = td->nb_coeffs[0]; + s->basis_order[2][1] = td->nb_coeffs[1]; + + ret = posix_memalign((void**)&s->coeffs, 32, + sizeof(*s->coeffs) * NB_EQUATIONS * td->nb_coeffs[0] * td->nb_coeffs[1]); + if (ret) + return -ENOMEM; + memset(s->coeffs, 0, sizeof(*s->coeffs) * 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]; + + for (int i = 0; i < ARRAY_ELEMS(basis_sets); i++) + for (int j = 0; j < ARRAY_ELEMS(basis_sets[i]); j++) { + double sf; + + ret = tdi_basis_init(&s->basis[i][j], basis_sets[i][j], td->basis_scale_factor[j]); + if (ret < 0) + return ret; + } + + + ret = tdi_nlsolve_context_alloc(&s->solver, ARRAY_ELEMS(basis_sets)); + if (ret < 0) { + tdi_log(&s->logger, 0, "Error allocating the non-linear solver\n"); + return ret; + } + + s->solver->logger = s->logger; + s->solver->tp = s->tp; + s->solver->maxiter = td->max_iter; + s->solver->atol = td->atol; + + memcpy(s->solver->basis, s->basis, sizeof(s->basis)); + memcpy(s->solver->solve_order, s->basis_order, sizeof(s->basis_order)); + +#if HAVE_OPENCL + s->solver->ocl_ctx = s->ocl_ctx; + s->solver->ocl_queue = s->ocl_queue; +#endif + + ret = tdi_nlsolve_context_init(s->solver); + if (ret < 0) { + tdi_log(&s->logger, 0, "Error initializing the non-linear solver\n"); + return ret; + } + + 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) +{ + const double *amplitude = opaque; + + return tdi_constraints_eval(eq_idx, *amplitude, colloc_grid_order, colloc_grid, + vars, dst); +} + +#define AMPLITUDE_DIRECT_SOLVE 1e-3 + +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 (coeffs_init) { + for (int i = 0; i < 3; i++) { + memcpy(td->coeffs[i], coeffs_init[i], sizeof(*td->coeffs[i]) * + td->nb_coeffs[0] * td->nb_coeffs[1]); + } + } + + while (1) { + double scale_factor; + + ret = tdi_nlsolve_solve(s->solver, teukolsky_constraint_eval, NULL, + &cur_amplitude, s->coeffs); + if (ret < 0) + return ret; + + if (fabs(cur_amplitude - td->amplitude) < DBL_EPSILON) + goto finish; + + scale_factor = cur_amplitude; + cur_amplitude = MIN(td->amplitude, 2 * cur_amplitude); + scale_factor = cur_amplitude / scale_factor; + + 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); + } + +finish: + tdi_solver_print_stats(s->solver); + + return 0; +} + +TDContext *td_context_alloc(void) +{ + TDContext *td = calloc(1, sizeof(*td)); + + if (!td) + return NULL; + + td->nb_threads = 1; + + td->amplitude = 1.0; + + td->nb_coeffs[0] = 32; + td->nb_coeffs[1] = 16; + + td->basis_scale_factor[0] = 3.0; + td->basis_scale_factor[1] = 3.0; + + td->max_iter = 16; + td->atol = 1e-12; + + td->log_callback = tdi_log_default_callback; + + td->priv = calloc(1, sizeof(TDPriv)); + if (!td->priv) { + free(td); + return NULL; + } + + return td; +} + +void td_context_free(TDContext **ptd) +{ + TDContext *td = *ptd; + TDPriv *s; + + if (!td) + return; + + s = td->priv; + + tdi_nlsolve_context_free(&s->solver); + tdi_threadpool_free(&s->tp); + +#if HAVE_OPENCL + if (s->ocl_queue) + clReleaseCommandQueue(s->ocl_queue); + if (s->ocl_ctx) + clReleaseContext(s->ocl_ctx); +#endif + + for (int i = 0; i < ARRAY_ELEMS(s->basis); i++) + for (int j = 0; j < ARRAY_ELEMS(s->basis[i]); j++) + tdi_basis_free(&s->basis[i][j]); + + free(s->coeffs); + + free(td->priv); + free(td); + *ptd = NULL; +} + +static double scalarproduct_metric_c(size_t len1, size_t len2, double *mat, + double *vec1, double *vec2) +{ + double val = 0.0; + for (int m = 0; m < len2; m++) { + double tmp = 0.0; + for (int n = 0; n < len1; n++) + tmp += mat[m * len1 + n] * vec1[n]; + val += tmp * vec2[m]; + } + return val; +} + +static int eval_var(const TDContext *td, unsigned int var_idx, + size_t nb_coords, const double *r, const double *theta, + const unsigned int diff_order[2], + double *out) +{ + TDPriv *s = td->priv; + double *basis_val[2] = { NULL }; + int ret = 0; + + if (diff_order[0] > 0 || diff_order[1] > 0) { + tdi_log(&s->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 = (var_idx == 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); + 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); + + val += scalarproduct_metric_c(td->nb_coeffs[0], td->nb_coeffs[1], td->coeffs[var_idx], + basis_val[0], basis_val[1]); + + out[i] = val; + } + + +fail: + free(basis_val[0]); + free(basis_val[1]); + + return ret; +} + +int td_eval_psi(const TDContext *td, + size_t nb_coords, const double *r, const double *theta, + const unsigned int diff_order[2], + double *out) +{ + return eval_var(td, 0, nb_coords, r, theta, diff_order, out); +} + +int td_eval_krr(const TDContext *td, + size_t nb_coords, const double *r, const double *theta, + const unsigned int diff_order[2], + double *out) +{ + return eval_var(td, 1, nb_coords, r, theta, diff_order, out); +} +int td_eval_kpp(const TDContext *td, + size_t nb_coords, const double *r, const double *theta, + const unsigned int diff_order[2], + double *out) +{ + return eval_var(td, 2, nb_coords, r, theta, diff_order, out); +} diff --git a/libteukolskydata.v b/libteukolskydata.v new file mode 100644 index 0000000..18295a2 --- /dev/null +++ b/libteukolskydata.v @@ -0,0 +1,4 @@ +LIBTEUKOLSKYDATA_0 { + global: td_*; + local: *; +}; diff --git a/td_constraints.c b/td_constraints.c new file mode 100644 index 0000000..9ef6a80 --- /dev/null +++ b/td_constraints.c @@ -0,0 +1,195 @@ +/* + * Copyright 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 . + */ + +#include + +#include "common.h" +#include "pssolve.h" +#include "td_constraints.h" + +#define I_C sqrt(225.0 / (8.0 * M_PI)) +static double I_m2(double x, double a) +{ + return a * I_C * exp(-SQR(x)); +} +static double I_m1(double x, double a) +{ + return a * I_C * (-2.0 * x * exp(-SQR(x))); +} +static double I_0(double x, double a) +{ + return a * I_C * (-2.0 + 4.0 * SQR(x)) * exp(-SQR(x)); +} +static double I_p1(double x, double a) +{ + return a * I_C * (12.0 * x - 8.0 * x * SQR(x)) * exp(-SQR(x)); +} +static double I_p2(double x, double a) +{ + return a * I_C * (12.0 - 48.0 * SQR(x) + 16.0 * SQR(SQR(x))) * exp(-SQR(x)); +} +static double I_p3(double x, double a) +{ + return a * I_C * (-120.0 * x + 144.0 * x * SQR(x) - 32.0 * x * SQR(SQR(x))) * exp(-SQR(x)); +} + +static double K_rtheta(double a, double r, double theta) +{ + double r2 = r * r; + double r3 = r2 * r; + + // equal to the below + return a * I_C * exp(-r2) * sin(2 * theta) * (32.0 * r3 - 48.0 * r); + + //return ((I_p2(-r, a) + I_p2(r, a)) * r3 + 3.0 * (I_p1(-r, a) - I_p1(r, a)) * r2 + + // 6.0 * (I_0(-r, a) + I_0(r, a)) * r + 6.0 * (I_m1(-r, a) - I_m1(r, a))) * sin(2.0 * theta) / r4; +} + +static double K_rtheta_dtheta(double a, double r, double theta) +{ + return 2.0 * cos(2.0 * theta) * K_rtheta(a, r, theta) / sin(2.0 * theta); +} + +static double K_rtheta_dr(double a, double r, double theta) +{ + double r2 = r * r; + return -2 * r * K_rtheta(a, r, theta) + a * I_C * exp(-r2) * sin(2 * theta) * (3.0 * 32.0 * r2 - 48.0); + //return ((-I_p3(-r, a) + I_p3(r, a)) / r + 3.0 * (-I_p2(-r, a) - I_p2(r, a)) / (SQR(r)) + + // 6.0 * (-I_p1(-r, a) + I_p1(r, a)) / (r * SQR(r)) + 6.0 * (-I_0(-r, a) - I_0(r, a)) / (SQR(SQR(r))) + + // (I_p2(-r, a) + I_p2(r, a)) / (-SQR(r)) + 3.0 * (I_p1(-r, a) - I_p1(r, a)) * (-2.0 / (r * SQR(r))) + + // 6.0 * (I_0(-r, a) + I_0(r, a)) * (-3.0 / (SQR(SQR(r)))) + 6.0 * (I_m1(-r, a) - I_m1(r, a)) * (-4.0 / (r * SQR(SQR(r))))) * sin(2.0 * theta); +} + +#include "eval_constraints.c" +#include "eval_k_rtheta.c" + +//typedef void (*eval_metric_func)(double a, double L, double r_0, double x_0, double x_1, double psi_val_00, double psi_val_01, double psi_val_02, double psi_val_10, double psi_val_11, double psi_val_20, double *res); +// +//static const eval_metric_func eval_dmetric[3] = { +// eval_dmetric_0, +// eval_dmetric_1, +// eval_dmetric_2, +//}; +// +//static const eval_metric_func eval_d2metric[3][3] = { +// { +// eval_d2metric_00, +// eval_d2metric_01, +// eval_d2metric_02, +// }, +// { +// eval_d2metric_10, +// eval_d2metric_11, +// eval_d2metric_12, +// }, +// { +// eval_d2metric_20, +// eval_d2metric_21, +// eval_d2metric_22, +// }, +//}; + +#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], + 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, +}; + +static double q_gundlach(const double a, double rho, double z) +{ + return a * SQR(rho) * exp(-(SQR(rho) + SQR(z))); +} + +static double dq_rho_gundlach(const double a, double rho, double z) +{ + return a * 2 * rho * exp(-SQR(rho) - SQR(z)) * (1.0 - rho * (rho)); +} + +static double d2q_rho_gundlach(const double a, double rho, double z) +{ + double rho2 = SQR(rho); + return a * 2 * exp(-SQR(rho) - SQR(z)) * (1.0 - 4.0 * rho * (rho) - rho2 + 2.0 * rho2 * SQR(rho)); +} + +static double d2q_rho_z_gundlach(const double a, double rho, double z) +{ + return a * 4 * z * rho * exp(-SQR(rho) - SQR(z)) * (rho * (rho) - 1.0); +} + +static double dq_z_gundlach(const double a, double rho, double z) +{ + return -a * 2 * z * SQR(rho) * exp(-SQR(rho) - SQR(z)); +} + +static double d2q_z_gundlach(const double a, double rho, double z) +{ + return a * 2 * SQR(rho) * exp(-SQR(rho) - SQR(z)) * (2 * SQR(z) - 1.0); +} + + +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) +{ + //for (int idx1 = 0; idx1 < nb_coords[1]; idx1++) { + // const double theta = coords[1][idx1]; + + // for (int idx0 = 0; idx0 < nb_coords[0]; idx0++) { + // const double r = coords[0][idx0]; + // const double rho = r * sin(theta); + // const double z = r * cos(theta); + + // const double ct = cos(theta); + // const double st = sin(theta); + + // const int idx = idx1 * nb_coords[0] + idx0; + + // const double laplace = vars[eq][PSSOLVE_DIFF_ORDER_20][idx] + + // vars[eq][PSSOLVE_DIFF_ORDER_10][idx] * (2.0 / r) + + // vars[eq][PSSOLVE_DIFF_ORDER_02][idx] * (1.0 / SQR(r)) + + // vars[eq][PSSOLVE_DIFF_ORDER_01][idx] * (ct / (SQR(r) * st)); + + // out[idx] = laplace + 0.25 * (d2q_z_gundlach(a, rho, z) + d2q_rho_gundlach(a, rho, z)) * (1.0 + vars[eq][PSSOLVE_DIFF_ORDER_00][idx]); + // } + //} + constraint_funcs[eq](a, nb_coords, coords, vars, out); + return 0; +} diff --git a/td_constraints.h b/td_constraints.h new file mode 100644 index 0000000..44f5257 --- /dev/null +++ b/td_constraints.h @@ -0,0 +1,43 @@ +/* + * Copyright 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 . + */ + +#ifndef TEUKOLSKY_DATA_CONSTRAINTS_H +#define TEUKOLSKY_DATA_CONSTRAINTS_H + +#include "pssolve.h" + +enum TDConstraintEq { + TD_CONSTRAINT_EQ_HAM = 0, + TD_CONSTRAINT_EQ_MOM_0, + TD_CONSTRAINT_EQ_MOM_1, + TD_CONSTRAINT_EQ_NB, +}; + +enum TDConstraintVar { + TD_CONSTRAINT_VAR_CONFFAC = 0, + TD_CONSTRAINT_VAR_CURV_RR, + TD_CONSTRAINT_VAR_CURV_PHIPHI, + 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); + +#endif // TEUKOLSKY_DATA_CONSTRAINTS_H diff --git a/td_constraints_template.c b/td_constraints_template.c new file mode 100644 index 0000000..df09290 --- /dev/null +++ b/td_constraints_template.c @@ -0,0 +1,198 @@ +/* + * 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 + } + } +} diff --git a/teukolsky_data.h b/teukolsky_data.h new file mode 100644 index 0000000..78bf824 --- /dev/null +++ b/teukolsky_data.h @@ -0,0 +1,165 @@ +/* + * Copyright 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 . + */ + +#ifndef TEUKOLSKY_DATA_H +#define TEUKOLSKY_DATA_H + +#include +#include +#include + +/** + * API usage: + * + * First, allocate the solver context with td_context_alloc(). All interaction + * with the solver is done through this context. + * + * Fill any fields in the context that are described as settable by the caller. + * Call td_solve() to solve the equation determined by the option values. + * + * The metric/extrinsic curvature can be evaluated with td_eval_metric()/td_eval_curv(). + * + * Finally, free the solver context with td_context_free(). + */ + +typedef struct TDContext { + /** + * Solver internals, not to be accessed by the caller + */ + void *priv; + + /******************************** + * options, set by the caller * + ********************************/ + + /** + * A callback that will be used to print diagnostic messages. + * + * Defaults to fprintf(stderr, ...) + */ + void (*log_callback)(const struct TDContext *td, int level, + const char *fmt, va_list); + + /** + * Arbitrary user data, e.g. to be used by the log callback. + */ + void *opaque; + + /******************************** + * initial data parameters * + ********************************/ + + /** + * The amplitude in the I function. + * Defaults to 1. + */ + double amplitude; + + /******************************** + * solver options * + ********************************/ + + /** + * The number of basis functions in each direction. + * [0] - radial, [1] - angular + */ + unsigned int nb_coeffs[2]; + + /** + * The scaling factor used in the basis functions. + * [0] - radial, [1] - angular + */ + double basis_scale_factor[2]; + + /** + * Maximum number of Newton iterations. + */ + unsigned int max_iter; + + /** + * Absolute tolerance. The solver is deemed to have converged + * after maximum difference between iterations is below this bound. + */ + double atol; + + /** + * Number of threads to use for parallelization. Set to 0 to automatically + * detect according to the number of CPU cores. + */ + unsigned int nb_threads; + + double *coeffs[3]; +} TDContext; + +/** + * Allocate and initialize the solver. + */ +TDContext *td_context_alloc(void); + +/** + * Free the solver and everything associated with it. + */ +void td_context_free(TDContext **td); + +/** + * Solve the equation for the conformal factor ψ and export the expansion coefficients in the + * context. + * + * @return >= 0 on success, a negative error code on failure + */ +int td_solve(TDContext *td, double *coeffs_init[3]); + +/** + * Evaluate the 3-metric γ_ij at the specified rectangular grid (in cylindrical + * coordinates { ρ, z, φ }). + * + * @param td the solver context + * @param rho the array of ρ coordinates. + * @param nb_coords_rho the number of elements in rho. + * @param z the array of z coordinates. + * @param nb_coords_z the number of elements in z. + * @param nb_comp number of needed components of the metric + * @param comp a nb_comp-sized array specifying the components of the metric to evaluate + * @param diff_order a nb_comp-sized array specifying the order of the derivatives of the + * metric to evaluate. The first element specifies the derivative wrt ρ, + * the second wrt z. I.e. diff_order[i] = { 0, 0 } evaluates the metric + * itself, diff_order[i] = { 0, 1 } evaluates ∂γ/∂z etc. + * @param out a nb_comp-sized array of pointers to the arrays into which the values of the + * metric will be written. Each requested component is evaluated on the grid + * formed by the outer product of the rho and z vectors. I.e. + * out[l][j * out_strides[l] + i] = γ_comp[l](rho[i], z[j]). The length of each + * array in out must be nb_coords_rho * nb_coords_z. + * @param out_strides a nb_comp-sized array of distances (in double-sized elements), for each + * array in out, between two elements corresponding to the same ρ but one + * step in z. Each element in out_strides must be at least nb_coords_rho. + * + * @return >= 0 on success, a negative error code on failure. + */ +int td_eval_psi(const TDContext *td, + size_t nb_coords, const double *r, const double *theta, + const unsigned int diff_order[2], + double *out); +int td_eval_krr(const TDContext *td, + size_t nb_coords, const double *r, const double *theta, + const unsigned int diff_order[2], + double *out); +int td_eval_kpp(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 new file mode 100644 index 0000000..5bc9364 --- /dev/null +++ b/teukolsky_data.py @@ -0,0 +1,128 @@ +# +# Copyright 2014 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 . +# + + +import ctypes +import numpy as np + +class TeukolskyData(object): + + coeffs = None + + _libtd = None + _tdctx = None + + class _TDContext(ctypes.Structure): + _fields_ = [("priv", ctypes.c_void_p), + ("log_callback", ctypes.c_void_p), + ("opaque", ctypes.c_void_p), + ("amplitude", ctypes.c_double), + ("nb_coeffs", ctypes.c_uint * 2), + ("basis_scale_factor", ctypes.c_double * 2), + ("max_iter", ctypes.c_uint), + ("atol", ctypes.c_double), + ("nb_threads", ctypes.c_uint), + ("coeffs", ctypes.POINTER(ctypes.c_double) * 3), + ] + + def __init__(self, **kwargs): + self._libtd = ctypes.CDLL('libteukolskydata.so') + tdctx_alloc = self._libtd.td_context_alloc + tdctx_alloc.restype = ctypes.POINTER(self._TDContext) + self._tdctx = self._libtd.td_context_alloc() + + coeffs_init = ctypes.c_void_p() + + for arg, value in kwargs.iteritems(): + if arg == 'coeffs_init': + coeffs_init = (ctypes.POINTER(ctypes.c_double) * 3)() + coeffs_init[0] = ctypes.cast(np.ctypeslib.as_ctypes(value[0]), ctypes.POINTER(ctypes.c_double)) + coeffs_init[1] = ctypes.cast(np.ctypeslib.as_ctypes(value[1]), ctypes.POINTER(ctypes.c_double)) + coeffs_init[2] = ctypes.cast(np.ctypeslib.as_ctypes(value[2]), ctypes.POINTER(ctypes.c_double)) + continue + + try: + self._tdctx.contents.__setattr__(arg, value) + except TypeError as e: + # try assigning items of an iterable + try: + for i, it in enumerate(value): + self._tdctx.contents.__getattribute__(arg)[i] = it + except: + raise e + + ret = self._libtd.td_solve(self._tdctx, coeffs_init) + if ret < 0: + raise RuntimeError('Error solving the equation') + + self.coeffs = [None] * 3 + for i in xrange(3): + self.coeffs[i] = np.copy(np.ctypeslib.as_array(self._tdctx.contents.coeffs[i], (self._tdctx.contents.nb_coeffs[1], self._tdctx.contents.nb_coeffs[0]))) + + def __del__(self): + if self._tdctx: + addr_tdctx = ctypes.c_void_p(ctypes.addressof(self._tdctx)) + self._libtd.td_context_free(addr_tdctx) + self._tdctx = None + + def _eval_var(self, eval_func, r, theta, diff_order = None): + if diff_order is None: + diff_order = [0, 0] + + c_diff_order = (ctypes.c_uint * 2)() + c_diff_order[0] = diff_order[0] + c_diff_order[1] = diff_order[1] + + if r.ndim == 2: + if r.shape != theta.shape: + raise TypeError('r and theta must be identically-shaped 2-dimensional arrays') + R, Theta = r.view(), theta.view() + elif r.ndim == 1: + if theta.ndim != 1: + raise TypeError('r and theta must both be 1-dimensional NumPy arrays') + R, Theta = np.meshgrid(r, theta) + else: + raise TypeError('invalid r/theta parameters') + + out = np.empty(R.shape[0] * R.shape[1]) + + R.shape = out.shape + Theta.shape = out.shape + + c_out = np.ctypeslib.as_ctypes(out) + c_r = np.ctypeslib.as_ctypes(R) + c_theta = np.ctypeslib.as_ctypes(Theta) + + ret = eval_func(self._tdctx, out.shape[0], c_r, c_theta, + c_diff_order, c_out, ctypes.c_long(r.shape[0])) + if ret < 0: + raise RuntimeError('Error evaluating the variable') + + out.shape = (theta.shape[0], r.shape[0]) + return out + + def eval_psi(self, r, theta, diff_order = None): + return self._eval_var(self._libtd.td_eval_psi, r, theta, diff_order) + def eval_krr(self, r, theta, diff_order = None): + return self._eval_var(self._libtd.td_eval_krr, r, theta, diff_order) + def eval_kpp(self, r, theta, diff_order = None): + return self._eval_var(self._libtd.td_eval_kpp, r, theta, diff_order) + + + @property + def amplitude(self): + return self._tdctx.contents.__getattribute__('amplitude') -- cgit v1.2.3