# # 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])