diff options
Diffstat (limited to 'eq_gen/constraint_gen.py')
| -rw-r--r-- | eq_gen/constraint_gen.py | 62 |
1 files changed, 30 insertions, 32 deletions
diff --git a/eq_gen/constraint_gen.py b/eq_gen/constraint_gen.py index b7033e8..64003f0 100644 --- a/eq_gen/constraint_gen.py +++ b/eq_gen/constraint_gen.py @@ -67,18 +67,18 @@ for diff0 in range(psi_funcs.shape[0]): 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 = sympy.Piecewise((a * kappa_p * (L ** 5) * ((1 - ((x / L) ** 2)) ** 6), sympy.Abs(x) < L), +# (0, True)) +Im2 = a * kappa_p * (L ** 5) * ((1 - ((x / L) ** 2)) ** 6) #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 +v = r - r_0 u = -r - r_0 Im2v = Im2.subs(x, v) @@ -93,8 +93,6 @@ 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() @@ -103,7 +101,6 @@ K_rtheta = ((Ip2v / (r ** 2) - 3 * Ip1v / (r ** 3) + 6 * I0v / (r ** 4) - 6 * Im #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], @@ -143,8 +140,6 @@ 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], @@ -168,29 +163,32 @@ def replace_derivatives(expr, replace): 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) +#sys.stderr.write('before\n') +#R_scalar_val = replace_derivatives(R_scalar, replace_diff).subs(replace_val).simplify() +#sys.stderr.write('1\n') +#Dcurv_m_r_val = replace_derivatives(Dcurv_m[0], replace_diff).subs(replace_val).simplify() +#sys.stderr.write('2\n') +#Dcurv_m_t_val = replace_derivatives(Dcurv_m[1], replace_diff).subs(replace_val).simplify() +#sys.stderr.write('3\n') + +#funcs = [] +# +#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_0, 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_0, r, theta) gen = cg.codegen(funcs, language = 'C', argument_sequence = argument_sequence) sys.stdout.write(gen[1][1]) sys.stdout.write(gen[0][1]) |