Add the remaining parts for teukolsky waves initial data

1 files changed, 196 insertions, 0 deletions

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 <anton@khirnov.net>
+#
+# 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 <http://www.gnu.org/licenses/>.
+#
+
+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])