From 6a0395d7e63372be04d20ee447b5d88af4f09f69 Mon Sep 17 00:00:00 2001 From: jthorn Date: Thu, 14 Jun 2001 11:53:21 +0000 Subject: This commit was generated by cvs2svn to compensate for changes in r2, which included commits to RCS files with non-trunk default branches. git-svn-id: http://svn.einsteintoolkit.org/cactus/EinsteinAnalysis/AHFinderDirect/trunk@3 f88db872-0e4f-0410-b76b-b9085cfa78c5 --- src/gr/H.c | 159 ++++++++++++++++++++++++++++++++++++++++++ src/gr/Makefile | 29 ++++++++ src/gr/gr_gfas.minc | 25 +++++++ src/gr/horizon_function.maple | 127 +++++++++++++++++++++++++++++++++ src/gr/setup_gfas.maple | 49 +++++++++++++ 5 files changed, 389 insertions(+) create mode 100644 src/gr/H.c create mode 100644 src/gr/Makefile create mode 100644 src/gr/gr_gfas.minc create mode 100644 src/gr/horizon_function.maple create mode 100644 src/gr/setup_gfas.maple (limited to 'src/gr') diff --git a/src/gr/H.c b/src/gr/H.c new file mode 100644 index 0000000..c345f78 --- /dev/null +++ b/src/gr/H.c @@ -0,0 +1,159 @@ + t1 = g_uu[2][2]; + t2 = 1/r; + t4 = X_ud[0][2]; + t5 = Diff(h,rho); + t7 = X_ud[1][2]; + t8 = Diff(h,sigma); + t10 = zz*t2-t4*t5-t7*t8; + t11 = t1*t10; + t12 = g_uu[1][1]; + t14 = X_ud[0][1]; + t16 = X_ud[1][1]; + t18 = yy*t2-t14*t5-t16*t8; + t19 = t12*t18; + t21 = r*r; + t23 = 1/r/t21; + t30 = Diff(h,rho,rho); + t33 = Diff(h,rho,sigma); + t38 = Diff(h,sigma,sigma); + t40 = -yy*zz*t23-X_udd[0][1][2]*t5-X_udd[1][1][2]*t8-t14*t4*t30-t16*t4* +t33-t14*t7*t33-t16*t7*t38; + t44 = g_uu[0][2]; + t45 = t44*t44; + t46 = t10*t10; + t48 = yy*yy; + t49 = zz*zz; + t56 = X_ud[0][0]; + t57 = t56*t56; + t59 = X_ud[1][0]; + t63 = t59*t59; + t65 = (t48+t49)*t23-X_udd[0][0][0]*t5-X_udd[1][0][0]*t8-t57*t30-2.0*t59* +t56*t33-t63*t38; + t70 = xx*t2-t56*t5-t59*t8; + t71 = t70*t70; + t72 = t70*t71; + t74 = Diff(g_uu[0][0],zz); + t77 = g_uu[1][2]; + t78 = t77*t10; + t79 = xx*xx; + t86 = t14*t14; + t91 = t16*t16; + t93 = (t79+t49)*t23-X_udd[0][1][1]*t5-X_udd[1][1][1]*t8-t86*t30-2.0*t16* +t14*t33-t91*t38; + t97 = g_uu[0][1]; + t99 = Diff(g_uu[0][0],yy); + t102 = t77*t18; + t109 = t4*t4; + t114 = t7*t7; + t116 = (t79+t48)*t23-X_udd[0][2][2]*t5-X_udd[1][2][2]*t8-t109*t30-2.0*t7* +t4*t33-t114*t38; + t120 = t44*t10; + t135 = -xx*zz*t23-X_udd[0][0][2]*t5-X_udd[1][0][2]*t8-t56*t4*t30-t59*t4* +t33-t56*t7*t33-t59*t7*t38; + t139 = g_uu[0][0]; + t140 = t139*t70; + t141 = Diff(g_uu[1][2],xx); + t142 = t141*t10; + t145 = t44*t70; + t146 = Diff(g_uu[1][2],zz); + t147 = t146*t10; + t151 = Diff(g_uu[0][0],xx); + t154 = t10*t46; + t156 = Diff(g_uu[2][2],zz); + t159 = t97*t70; + t160 = Diff(g_uu[1][2],yy); + t161 = t160*t10; + t164 = -2.0*t11*t19*t40-t45*t46*t65-t44*t72*t74/2.0-2.0*t78*t19*t93-t97* +t72*t99/2.0-2.0*t11*t102*t116-2.0*t120*t102*t135-t140*t142*t18-t145*t147*t18- +t139*t72*t151/2.0-t1*t154*t156/2.0-t159*t161*t18; + t166 = Diff(g_uu[2][2],yy); + t170 = Diff(g_uu[2][2],xx); + t173 = t77*t46; + t188 = -xx*yy*t23-X_udd[0][0][1]*t5-X_udd[1][0][1]*t8-t56*t14*t30-t59*t14 +*t33-t56*t16*t33-t59*t16*t38; + t192 = t18*t18; + t193 = t18*t192; + t195 = Diff(g_uu[1][1],zz); + t199 = Diff(g_uu[1][1],xx); + t202 = t195*t192; + t205 = t199*t192; + t209 = Diff(g_uu[1][1],yy); + t212 = t1*t1; + t215 = t1*t46; + t216 = t77*t40; + t219 = t97*t192; + t221 = t12*t192; + t223 = -t77*t154*t166/2.0-t44*t154*t170/2.0-2.0*t173*t44*t188-t77*t193* +t195/2.0-t97*t193*t199/2.0-t145*t202/2.0-t140*t205/2.0-t12*t193*t209/2.0-t212* +t46*t116-2.0*t215*t216-t219*t142-t221*t161; + t225 = t74*t71; + t228 = t99*t71; + t235 = t151*t71; + t238 = t156*t46; + t241 = t77*t192; + t243 = t209*t192; + t246 = t77*t77; + t249 = t97*t18; + t252 = t170*t46; + t255 = t166*t46; + t258 = -t102*t225/2.0-t19*t228/2.0-t11*t225/2.0-t78*t228/2.0-t120*t235/ +2.0-t145*t238/2.0-t241*t147-t159*t243/2.0-t246*t46*t93-t249*t235/2.0-t140*t252/ +2.0-t159*t255/2.0; + t263 = Diff(g_uu[0][1],zz); + t266 = Diff(g_uu[0][1],yy); + t267 = t266*t70; + t269 = t139*t71; + t270 = Diff(g_uu[0][1],xx); + t271 = t270*t18; + t273 = t97*t71; + t274 = Diff(g_uu[0][2],yy); + t275 = t274*t10; + t277 = Diff(g_uu[0][2],xx); + t278 = t277*t10; + t282 = t266*t18; + t285 = t263*t18; + t290 = t44*t71; + t291 = Diff(g_uu[0][2],zz); + t292 = t291*t10; + t297 = -t249*t252/2.0-t102*t238/2.0-t241*t263*t70-t221*t267-t269*t271- +t273*t275-t269*t278-t19*t255/2.0-t78*t282*t70-t11*t285*t70-t120*t271*t70-t290* +t292-2.0*t273*t139*t188; + t303 = t97*t97; + t313 = t97*t188; + t316 = t12*t12; + t332 = -2.0*t290*t97*t40-t303*t71*t93-2.0*t290*t139*t135-t45*t71*t116- +t303*t192*t65-2.0*t221*t313-t316*t192*t93-2.0*t241*t97*t135-2.0*t241*t12*t40 +-2.0*t45*t10*t70*t135-t246*t192*t116-t273*t282; + t333 = t139*t139; + t336 = t145*t40; + t351 = t44*t46; + t357 = t291*t70; + t362 = t159*t40; + t365 = -t333*t71*t65-2.0*t78*t336-t249*t278*t70-t102*t292*t70-t19*t275* +t70-2.0*t11*t249*t135-2.0*t78*t249*t188-t351*t277*t70-t173*t274*t70-t290*t285- +t215*t357-2.0*t120*t19*t188-2.0*t11*t362; + t371 = t140*t65; + t374 = t140*t135; + t379 = t140*t188; + t382 = t159*t93; + t401 = -2.0*t246*t10*t18*t40-2.0*t120*t371-2.0*t11*t374-t219*t270*t70-2.0 +*t78*t379-2.0*t78*t382-2.0*t120*t159*t188-2.0*t19*t382-2.0*t303*t18*t70*t188 +-2.0*t19*t379-2.0*t249*t371-2.0*t249*t145*t135; + t406 = t145*t116; + t422 = t44*t135; + t427 = t146*t18; + t431 = -2.0*t102*t362-2.0*t102*t374-2.0*t11*t406-2.0*t102*t406-2.0*t19* +t336-t11*t202/2.0-t78*t243/2.0-2.0*t120*t249*t65-t120*t205/2.0-2.0*t215*t422- +t351*t141*t18-t215*t427-t173*t160*t18; + t441 = t269+2.0*t249*t70+2.0*t120*t70+t221+2.0*t78*t18+t215; + t442 = sqrt(t441); + t452 = t151*t70+t267+t357+t271+t209*t18+t427+t278+t161+t156*t10+t139*t65+ +2.0*t313+2.0*t422; + t456 = partial_d_ln_sqrt_g[0]; + t459 = partial_d_ln_sqrt_g[1]; + t462 = partial_d_ln_sqrt_g[2]; + t477 = t12*t93+2.0*t216+t1*t116+t139*t456*t70+t97*t459*t70+t44*t462*t70+ +t97*t456*t18+t12*t459*t18+t77*t462*t18+t44*t456*t10+t77*t459*t10+t1*t462*t10; + H = (t164+t223+t258+t297+t332+t365+t401+t431)/t442/t441+(t452+t477)/t442+ +(K_uu[0][0]*t71+2.0*K_uu[0][1]*t18*t70+2.0*K_uu[0][2]*t10*t70+K_uu[1][1]*t192+ +2.0*K_uu[1][2]*t10*t18+K_uu[2][2]*t46)/t441-K; diff --git a/src/gr/Makefile b/src/gr/Makefile new file mode 100644 index 0000000..892450c --- /dev/null +++ b/src/gr/Makefile @@ -0,0 +1,29 @@ +# Makefile for GR code +# $Id: Makefile,v 1.1.1.1 2001-06-14 11:53:21 jthorn Exp $ + +# +# Environment Variables: +# MAPLE_VERSION used via @ifdef for version control in Maple code; +# typically set to something like MAPLE_V_RELEASE_4 +# +# Targets: +# mm ==> preprocess all *.maple files to produce *.mm files +# clean ==> delete *.mm (mpp Maple preprocessor output) +# run ==> load all the code into Maple +# + +############################################################################### + +.PHONY : mm +mm : $(patsubst %.maple, %.mm, $(wildcard *.maple)) + +%.mm : %.maple $(wildcard *.minc $(gfa_dir)/*.minc $(gfa_dir)/*.maple) + mpp -D$(MAPLE_VERSION) <$< >$@ + +.PHONY : clean +clean : + -rm *.mm + +.PHONY : run +run : mm + maple -f maple.log diff --git a/src/gr/gr_gfas.minc b/src/gr/gr_gfas.minc new file mode 100644 index 0000000..92e781e --- /dev/null +++ b/src/gr/gr_gfas.minc @@ -0,0 +1,25 @@ +# Maple include file listing all gfa functional-dependence forms + +g_dd, +K_dd, + +g_uu, g_uu__fnd, +K_uu, K_uu__fnd, + +K, K__fnd, + +partial_d_ln_sqrt_g, partial_d_ln_sqrt_g__fnd, + +h, h__fnd, + +s_d, s_d__fnd, +partial_d_s_d, partial_d_s_d__fnd, + +n_u, n_u__fnd, + +H, H__fnd, + +HA, HA__fnd, +HB, HB__fnd, +HC, HC__fnd, +HD, HD__fnd # no comma diff --git a/src/gr/horizon_function.maple b/src/gr/horizon_function.maple new file mode 100644 index 0000000..c1750ea --- /dev/null +++ b/src/gr/horizon_function.maple @@ -0,0 +1,127 @@ +# horizon.maple - evaluate LHS of horizon function +# +# non_unit_normal - compute s_d +# non_unit_normal_gradient - compute partial_d_s_d +# horizon_function - compute HA,HB,HC,HD,H = fn(s_d) +# + +################################################################################ +################################################################################ +################################################################################ + +# +# This function computes the non-unit outward-pointing normal vector +# (field) to the horizon, using the equation on p.7.1 of my apparent +# horizon finding notes. +# +non_unit_normal := +proc() +global + @include "../maple/coords.minc", + @include "../gr/gr_gfas.minc"; +local i, u; + +assert_fnd_exists(h); +assert_fnd_exists(X_ud); +assert_fnd_exists(s_d, fnd); + + for i from 1 to N + do + s_d__fnd[i] := x_xyz[i]/r + - sum('X_ud[u,i]*Diff(h,x_rrs[u])', 'u'=A_min..N); + end do; + +NULL; +end proc; + + +################################################################################ + +# +# This function computes the gradient of the non-unit outward-pointing +# normal vector (field) to the horizon, using the equation on p.7.2 of +# my apparent horizon finding notes. +# +non_unit_normal_gradient := +proc() +global + @include "../maple/coords.minc", + @include "../gr/gr_gfas.minc"; +local temp, i, j, k, u, v; + +assert_fnd_exists(h); +assert_fnd_exists(X_ud); +assert_fnd_exists(X_udd); +assert_fnd_exists(partial_d_s_d, fnd); + + for i from 1 to N + do + for j from 1 to N + do + temp := `if`(i=j, sum('x_xyz[k]^2','k'=1..3), 0); + # simplify likes to put things over a common denominator, + # so we only apply it to the individual terms, not to + # the whole experssion + partial_d_s_d__fnd[i,j] + := + simplify( (temp-x_xyz[i]*x_xyz[j])/r^3 ) + - simplify( sum('X_udd[u,i,j]*Diff(h,x_rrs[u])', 'u'=A_min..N) ) + - simplify( msum('X_ud[u,i]*X_ud[v,j]*Diff(h,x_rrs[u],x_rrs[v])', + 'u'=A_min..N, 'v'=A_min..N) ); + end do; + end do; + +NULL; +end proc; + +################################################################################ + +# +# This function computes the LHS of the apparent horizon equation +# as a function of 1st and 2nd angular partial derivatives of the +# apparent horizon radius $h$, in the Maple gridfn array H.fnd , +# using equation (15) in my 1996 apparent horizon finding paper. +# +horizon_function := +proc() +global + @include "../maple/coords.minc", + @include "../gr/gr_gfas.minc"; +local i,j,k,l; + +assert_fnd_exists(g_uu); +assert_fnd_exists(K_uu); +assert_fnd_exists(partial_d_ln_sqrt_g); +assert_fnd_exists(s_d, fnd); +assert_fnd_exists(partial_d_s_d, fnd); + +# +# n.b. no simplify() in these subexpressions; empirically it makes +# them much *more* complicated (measured by codegen(...,optimized)) +# + +# (15a) in my paper +HA__fnd := - msum('g_uu[i,k]*s_d__fnd[k] * g_uu[j,l]*s_d__fnd[l] + * partial_d_s_d__fnd[i,j]', + 'i'=1..N, 'j'=1..N, 'k'=1..N, 'l'=1..N) + - (1/2)*msum('g_uu[i,j]*s_d__fnd[j] + * Diff(g_uu[k,l],x_xyz[i]) * s_d__fnd[k]*s_d__fnd[l]', + 'i'=1..N, 'j'=1..N, 'k'=1..N, 'l'=1..N); + +# (15b) in my paper +HB__fnd := + msum('Diff(g_uu[i,j],x_xyz[i])*s_d__fnd[j]', 'i'=1..N, 'j'=1..N) + + msum('g_uu[i,j]*partial_d_s_d__fnd[i,j]', 'i'=1..N, 'j'=1..N) + + msum('g_uu[i,j]*partial_d_ln_sqrt_g[i]*s_d__fnd[j]', + 'i'=1..N, 'j'=1..N); + +# (15c) in my paper +HC__fnd := msum('K_uu[i,j]*s_d__fnd[i]*s_d__fnd[j]', 'i'=1..N, 'j'=1..N); + +# (15d) in my paper +HD__fnd := msum('g_uu[i,j]*s_d__fnd[i]*s_d__fnd[j]', 'i'=1..N, 'j'=1..N); + +# n.b. no simplify() here, it would try to put things over a common +# denominator, which would make the equation much more complicated +H__fnd := HA__fnd/HD__fnd^(3/2) + HB__fnd/HD__fnd^(1/2) + HC__fnd/HD__fnd - K; + +NULL; +end proc; diff --git a/src/gr/setup_gfas.maple b/src/gr/setup_gfas.maple new file mode 100644 index 0000000..5da990b --- /dev/null +++ b/src/gr/setup_gfas.maple @@ -0,0 +1,49 @@ +# Maple code to set up all gridfn arrays +# $Id$ + +# +# setup_gr_gfas - setup all the GR gfas +# + +################################################################################ + +setup_gr_gfas := +proc() +global + @include "../gr/gr_gfas.minc"; + +# basic metric & extrinsic curvature fields +make_gfa('g_dd', {inert}, [1..N, 1..N], symmetric); +make_gfa('K_dd', {inert}, [1..N, 1..N], symmetric); + +make_gfa('g_uu', {inert,fnd}, [1..N, 1..N], symmetric); +make_gfa('K_uu', {inert,fnd}, [1..N, 1..N], symmetric); + +make_gfa('K', {inert, fnd}, [], none); + +# partial derivatives of metric determinant +make_gfa('partial_d_ln_sqrt_g', {inert,fnd}, [1..N], none); + +# radius of horizon +make_gfa('h', {inert, fnd}, [], none); + +# outward-pointing *non*-unit normal to horizon +# and it's xyz-coordinate partial derivatives +make_gfa('s_d', {inert,fnd}, [1..N], none); +make_gfa('partial_d_s_d', {inert,fnd}, [1..N, 1..N], none); + +# outward-pointing unit normal to horizon +make_gfa('n_u', {inert,fnd}, [1..N], none); + +# LHS of apparent horizon equation +make_gfa('H', {inert, fnd}, [], none); + +# subexpressions for computing LHS of apparent horizon equation +# ... these are defined by (15) in my 1996 apparent horizon finding paper +make_gfa('HA', {inert, fnd}, [], none); +make_gfa('HB', {inert, fnd}, [], none); +make_gfa('HC', {inert, fnd}, [], none); +make_gfa('HD', {inert, fnd}, [], none); + +NULL; +end proc; -- cgit v1.2.3