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 --- eq_gen/constraint_gen.py | 196 +++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 196 insertions(+) create mode 100644 eq_gen/constraint_gen.py (limited to 'eq_gen') 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]) -- cgit v1.2.3