diff options
| -rw-r--r-- | README | 397 | ||||
| -rw-r--r-- | interface.ccl | 49 | ||||
| -rw-r--r-- | param.ccl | 422 | ||||
| -rw-r--r-- | schedule.ccl | 31 | ||||
| -rw-r--r-- | src/AHFinder.F | 1524 | ||||
| -rw-r--r-- | src/AHFinder_2dio.c | 14 | ||||
| -rw-r--r-- | src/AHFinder_calcsigma.F | 311 | ||||
| -rw-r--r-- | src/AHFinder_dat.F | 142 | ||||
| -rw-r--r-- | src/AHFinder_exp.F | 429 | ||||
| -rw-r--r-- | src/AHFinder_find3.F | 68 | ||||
| -rw-r--r-- | src/AHFinder_flow.F | 660 | ||||
| -rw-r--r-- | src/AHFinder_fun.F | 166 | ||||
| -rw-r--r-- | src/AHFinder_gau.F | 1168 | ||||
| -rw-r--r-- | src/AHFinder_initguess.F | 272 | ||||
| -rw-r--r-- | src/AHFinder_int.F | 1105 | ||||
| -rw-r--r-- | src/AHFinder_leg.F | 305 | ||||
| -rw-r--r-- | src/AHFinder_mask.F | 81 | ||||
| -rw-r--r-- | src/AHFinder_min.F | 281 | ||||
| -rw-r--r-- | src/AHFinder_pow.F | 1070 | ||||
| -rw-r--r-- | src/make.code.defn | 9 | ||||
| -rw-r--r-- | src/make.code.deps | 14 |
21 files changed, 8518 insertions, 0 deletions
@@ -0,0 +1,397 @@ +Cactus Code Thorn AHFinder +Authors : Miguel Alcubierre +Managed by : Miguel Alcubierre <miguel@aei-potsdam.mpg.de> +Version : 1.0 +-------------------------------------------------------------------------- + + +NOTE: THIS FILE IS SOMEWHAT OUT OF DATE. I WILL UPDATE IT WHEN +I HAVE MORE TIME. PLEASE LOOK AT THE FILE `MinimumAHF_param.h' +FOR A MORE UP TO DATA DESCRIPTION OF THE PARAMETERS. + + +1. Purpose of the thorn + +MAIN PURPOSE +------------ + +This thorn looks for apparent horizons (AH) in 3D. An AH is defined +as a surface where the expansion of outgoing null geodesics is zero. + + +BRIEF DESCRIPTION +----------------- + +This thorn looks for apparent horizons (AH) in 3D using two +different algorithms. In the default mode, it uses a minimization +algorithm. It can also use a flow algorithm if one sets the +parameter "mahf_flow" equal to "yes": + +mahf_flow = "yes" + +The minimization algorithm is VERY slow if one has a lot of terms +in the spherical harmonics expansion (see below), but it is quite +fast for a small number of terms. + +The thorn is activated by adding the following lines to the +parameter file: + +ah_persists = "yes" +minimumAHF_active = "yes" + +Also, one needs to tell the routine how often we want to +look for AH by using the parameter `mahf_findevery' (its +default value is 1). + +The routine defines the grid function `mahfgrid' as: + +mahfgrid = r - f(theta,phi) + +So that the surface being considered corresponds to the zero level of +`mahfgrid'. The function f(theta,phi) is expanded in spherical +harmonics Y(l,m). Notice, however, that I DO NOT use the standard +normalization factor: + +sqrt( (2l + 1) (l - m)! / 4 pi (l + m)! ) + +I use instead: + +sqrt( 2l + 1 ) m = 0 + +sqrt( 2 (2l + 1) (l - m)! / (l + m)! ) m != 0 + +Notice first that I leave out the 4*pi. This is because I want Y(0,0) +to be the real radius of a sphere. Also, for M non-zero I have an +extra sqrt(2). I need this because I use a real basis with sines and +cosines instead of complex exponentials. + + +The minimization algorithm has three main modes of operation: + +* Default mode: In the default mode, the routine tries to +find a minimum of the integral of the square of the expansion. +At an AH, this integral should be zero, and for any other +surface it will be positive. One must remember that the +routine can easily end up in a local minimum that is not +a horizon. To try to prevent this, the routine first +looks through a subset of the parameter space to try +to find a good guess for the global minimum. This, of +course can fail. + +* Area minimization mode: This mode is activated by +setting the parameter: + +mahf_minarea = "yes" + +In this mode, the routine looks for a local minimum of the surface +area. For time-symmetric data (and only then), such a minimum will +correspond to an AH. + +* Trapped surface finding mode: This mode is activated by +setting the parameter: + +mahf_trapped_surface = "yes" + +In this mode, the routine looks for a local minimum of the integral of +the square of the quantity (expansion + trapped_surface_delta). If +trapped_surface_delta is positive, then the routine will try to find a +surface whose expansion is -(trapped_surface_delta) everywhere on the +surface. One can monitor the number of interpolated points that have +a negative expansion by setting veryverbose = "yes", which indicates +whether or not a trapped surface has truly been found. + +OUTPUT +------ + +The routine controls its own output. It produces output of 2 grid +functions (mahfgrid, mahf_exp) on its own (so please don't add them to +the Cactus output). The reason for this is that the finder might be +called at times that have nothing to do with when Cactus wants to do +output. + +The grid function `mahfgrid' was described above. The grid function +`mahf_exp' is the expansion of outgoing photons on the level sets of +`mahfgrid'. + +The finder also produces a series of data files: + +1) "mahf_coeff.alm" + +This file contains the coefficients of the spherical +harmonics expansion of the surface found. Notice +that the finder almost always finds a surface, but this +might not be a horizon, so be careful. + +The format in which these coefficients are saved is the following: + +* Lines beginning with # are comments. These include the time + at which the finder was called, and a line stating if the + surface found is a horizon or not. + +* Empty lines indicate different calls to the finder (used when + the finder is called many times during an evolution). + +* The coefficients are saved using a loop of the form: + + do l=1,lmax + do m=-l,l + write(*,*) a(l,m),l,m + end do + end do + + Notice that because our functions are real, I use an expansion + in sines and cosines and I really only consider positive values + of m. On output my convention is that coefficients with positive + m represent the terms: + + P_(l,m) cos(m theta) + + and coefficients with negative m represent the terms: + + P_(l,m) sin(m theta) + +2) "mahf.gauss" + + This file contains a map of the gaussian curvature on the surface. + The format of this file might still change in the near future. + I need to talk to Werner about this. + +* Lines beginning with # are comments. + +* The data is written is a loop: + + do i=0,ntheta-1 + do j=0,nphi-1 + write gaussian(i,j) + end do + end do + + theta and phi are subdivided uniformly + (according to grid type) in the intervals: + + octant: theta=[0,pi/2] phi=[0,pi/2] + quadrant: theta=[0,pi] phi=[0,pi/2] + full: theta=[0,pi] phi=[0,2 pi] + +3) "mahf_area.tl" + + Area of the surface for each time the finder was called. + Again, remember that the surface might not be a horizon. + +4) "mahf_mass.tl" + + Mass of the surface for each time the finder was called. + This is defined as: + + M = sqrt(A/(16 pi)) + + Again, remember that the surface might not be a horizon. + +5) "mahf_circ_eq.tl" + + Equatorial circumference of the surface. + Again, remember that the surface might not be a horizon. + +6) "mahf_meri_p1.tl" + + Length of meridian of surface for phi=0. + Again, remember that the surface might not be a horizon. + +7) "mahf_meri_p2.tl" + + Length of meridian of surface surface for phi=pi/2. + Again, remember that the surface might not be a horizon. + +8) "mahf_logfile" + + Log file for the last time the horizon was called. The + reason why only the last call is here is that this file + is very long. + + +PARAMETERS +---------- + +Other parameters for the thorn are (see also file MinimumAHF_param.h): + + +STRINGS: + +mahf_logfile = [yes,no] Write a log file. If yes, the + thorn produces the file "mahf_logfile". + (default = "no") + +mahf_verbose = [yes,no] Write messages to screen. + (default = "yes") + +mahf_veryverbose = [yes,no] Write lots of messages to the screen, + essentially the whole log file. + (default = "no") + +mahf_guessverbose = [yes,no] Write a little info for each trial point + calculated in the guess process. + (default = "no") + +mahf_2Doutput = [yes,no] 2D output of grid functions? + (default = "no") + +mahf_3Doutput = [yes,no] 3D output of grid functions? + (default = "no") + +mahf_mask = [yes,no] Use mask for definite horizons? + (default = "no") + +mahf_weakmask = [yes,no] Use mask for possible horizons? + (default = "no") + +mahf_phi = [yes,no] Expand in phi. It is useful to switch + this off if one knows in advance that + the horizon is really axisymmetric. + (default = "yes") + +mahf_offset = [yes,no] Is the center offset from origin? + (default = "no") + +mahf_wander = [yes,no] Do we allow the center to wander? + (default = "no") + +mahf_refx = [yes,no] Do we have reflection symmetry x -> -x? + (default = "no") + +mahf_refy = [yes,no] Do we have reflection symmetry y -> -y? + (default = "no") + +mahf_refz = [yes,no] Do we have reflection symmetry z -> -z? + (default = "no") + +mahf_octant = [yes,no,high] Does the surface have octant symmetry? + The "high" option forces also reflection + symmetry across the x-y diagonal. + (default = "no") + +mahf_areamap = [yes,no] Construct an area map?. + (default = "no") + +mahf_sloppyguess = [yes,no] This considers only spheres for the + initial guess. It is much faster. + (default = "no") + +mahf_guessold = [yes,no] Use horizon found in previous call as + initial guess? Only relevant for + evolutions. + (default = "no") + +mahf_inner = [yes,no] Tries to look for inner horizon instead. + (default = "no") + +mahf_guess_absmin = [yes,no] Use the absolute minimum of the surface + integral of the square of the expansion + as the place to start the minimization + process. This is useful if you have + a very dynamical spacetime, with possibly + many local minima and you don't have any + idea where (or even if) you have an AH. + (default = "no") + +mahf_manual_guess = [yes,no] Use a manually specified guess, via. the + parameters mahf_l0_guess, mahf_l2_guess, + etc. (works only in octant mode right now) + (default = "no") + +INTEGERS: + +mahf_findafter After how many timesteps start looking + for horizons (default 0). + +mahf_findevery How often to look for horizons (iterations). + (default 1) + +mahf_maxiter Maximum number of iterations of + Powell's minimization algorithm. + (default 10) + +mahf_flowiter Maximum number of iterations for flow + (default 200) + +mahf_lmax (<20) Number of terms in theta expansion + (default 2) + +mahf_ntheta Number of subdivisions in theta for + surface integrals (default 200). + +mahf_nphi Number of subdivisions in phi for + surface integrals (default 200). + +mahf_nn0 Number of subdivisions of c0(0) for + initial guess (default 10) + +mahf_nn2 Number of subdivisions of c0(2) for + initial guess (default 10) + + +REALS: + +mahf_findaftertime After what time start looking for horizons. + A non-zero value of this parameter overrides + the paramter `mahf_findafter' above. + (default 0.0) + +mahf_r0 Radius of initial sphere (0 forces largest + possible sphere) + (default 0.0) + +mahf_xc x coordinate of center of expansion + (default 0.0) + +mahf_yc y coordinate of center of expansion + (default 0.0) + +mahf_zc z coordinate of center of expansion + (default 0.0) + +mahf_tol Fractional tolerance for minimization + algorithm. If we decrease by less than + this in one iteration we are done. + (default = 0.1) + +trapped_surface_delta In 'mahf_trapped_surface' mode, this + determines the expansion of the surface + that one is looking for. Notice that a + positive value of 'mahf_trapped_surface' + will cause the finder to look for a + surface with expansion equal to MINUS + this value everywhere. + (default = 0.0) + +mahf_flowa Alpha parameter for flow (default 0.01). + +mahf_flowb Beta parameter for flow (default 0.01). + +mahf_flowh Weight of H flow (default 0.0). + +mahf_flowc Weight of C flow (default 1.0). + +mahf_flown Weight of N flow (default 0.0). + +mahf_flowtol Tolerance for flow (default 0.0001). + +mahf_maskshrink Shrink factor for mask (default -2.0). + + +IMPORTANT: Notice that the symmetry parameters refer to the surface +itself, and not to the grid. + + +2. Dependencies of the thorn + +This thorn additionally requires thorns: GenericAHF, util. +For testing purposes, one also needs: analyticBH, exact. + + +3. Thorn distribution + +This thorn is available to everyone. + + +4. Additional information + diff --git a/interface.ccl b/interface.ccl new file mode 100644 index 0000000..0bff3b9 --- /dev/null +++ b/interface.ccl @@ -0,0 +1,49 @@ +# Interface definition for thorn AHFinder +# $Header$ +#/*@@ +#c @file AHFinder_rfr.F +#c @date April 1998 +#c @author Miguel Alcubierre +#c @desc +#c RFR Initialization +#c @enddesc +#c@@*/ + +implements: ahfinder +inherits: interp einstein grid adm + + +public: + +real ahfindergrid type=GF +{ +ahfgrid, +ahf_exp +} "Grid function for single surface" + +real ahfmask type=GF +{ +ahmask +} "Grid function for masking" + +private: + +real ahfgradient type=GF +{ +ahfgradx, +ahfgrady, +ahfgradz, +ahfgradn +} "Grid functions for gradients" + +real ahfinder_gauss type=GF +{ +ahfgauss +} "Grid function for gaussian curcature calculation" + + +real find3grid type=GF +{ +ahfgrid3, +ahf_exp3 +} "Grid functions to use in find3 algorithm" diff --git a/param.ccl b/param.ccl new file mode 100644 index 0000000..c187673 --- /dev/null +++ b/param.ccl @@ -0,0 +1,422 @@ +# Parameter definitions for thorn AHFinder +# $Header$ + + +#These parameters are to switch on the finder, and to control how +#often and when to start looking for horizons. + +private: + +LOGICAL ahf_active "Activate AHFinder?" +{ +} "no" + +LOGICAL ahf_find3 "Searching for 3 surfaces?" +{ +} "no" + +LOGICAL ahf_trapped_surface "Minimize (expansion + delta) to find trapped surface?" +{ +} "no" + +INT ahf_findevery "How often to look for horizons" +{ + : :: "Ranges from 1 for each iteration to higher numbers with default set to 1" +} 1 + +INT ahf_findafter "After how many iterations look for horizons" +{ + : :: "By default set to 0, Any positive integer if start later in an evolution" +} 0 + +REAL ahf_findaftertime "After how much time look for horizons." +{ + : :: "Any positive real number. Default is 0.0. Overides ahf_findafter" +} 0.0 + +REAL trapped_surface_delta "find (expansion = delta) surface" +{ + : :: "Just a real number. Range ??" +} 0.0 + +#Parameters for surface expansion. + +LOGICAL ahf_phi "Expand also in phi? I.e seach for non-axisymmetric surface" +{ +} "no" + +LOGICAL ahf_offset "Center offset from origin?" +{ +} "no" + +LOGICAL ahf_wander "Allow the center to wander?" +{ +} "no" + +INT ahf_lmax "Maximum number of terms in theta expansion" +{ +0:20 :: "Range from 0 to 20. 8 is the default" +} 8 + +#Parameters for surface integrals. + +INT ahf_ntheta "Number of subdivisions in theta" +{ + : :: "Any senseful integer. The default of 200 is useful" +} 200 + +INT ahf_nphi "Number of subdivisions in phi" +{ + : :: "Any senseful integer. The default of 200 is useful" +} 200 + +#Parameters to indicate symmetries. This symmetries refer to the +#surface itself, and NOT the grid. + +#Notice that octant = "yes" enforces reflection symmetries on all +#three coordinate planes, while octant = "high" enforces also the +#extra symmetry of rotation of pi/2 around the z axis. + +LOGICAL ahf_refx "Reflection symmetry x->-x?" +{ +} "no" + +LOGICAL ahf_refy "Reflection symmetry y->-y?" +{ +} "no" + +LOGICAL ahf_refz "Reflection symmetry z->-z?" +{ +} "no" + +KEYWORD ahf_octant "Octant symmetry?" +{ + "no" :: "No octant symmetry" + "yes" :: "Octant symmetry: reflection symmetry on all three coordinate planes" + "high" :: "reflection symmetry on all three planes and symmetry of rotation of pi/2 around the z axis" +} "no" + +#Parameters for minimization algorithm. + +LOGICAL ahf_minarea "Minimize area instead of expansion?" +{ +} "no" + +INT ahf_maxiter "Maximum number of iterations of POWELL" +{ + : :: "Any senseful integer value. Default is 10" +} 10 + +REAL ahf_tol "Tolerance for minimization routines" +{ + : :: "A senseful real number." +} 0.1 + +#Parameters for initial guess. + +LOGICAL ahf_sloppyguess "Use sphere as initial guess?" +{ +} "no" + +LOGICAL ahf_guess_absmin "Use absolute min to start Minimization?" +{ +} "no" + +LOGICAL ahf_guessold "Use old horizon as initial guess?" +{ +} "no" + +LOGICAL ahf_inner "Look for inner horizon?" +{ +} "no" + +LOGICAL ahf_manual_guess "Use specified coefficients for guess?" +{ +} "no" + +INT ahf_nn0 "Number of subdivisions of c0(0) for initial guess" +{ + : :: "Some positive integer" +} 10 + +INT ahf_nn2 "Number of subdivisions of c0(2) for initial guess" +{ + : :: "Some positive integer" +} 10 + +REAL ahf_l0_guess "Manual guess for l=0 component" +{ + : :: +} 1.0 + +REAL ahf_l2_guess "Manual guess for l=2 component" +{ + : :: +} 0.0 + +REAL ahf_l4_guess "Manual guess for l=4 component" +{ + : :: +} 0.0 + +REAL ahf_l6_guess "Manual guess for l=6 component" +{ + : :: +} 0.0 + +REAL ahf_l8_guess "Manual guess for l=8 component" +{ + : :: +} 0.0 + +REAL ahf_l10_guess "Manual guess for l=10 component" +{ + : :: +} 0.0 + +REAL ahf_l12_guess "Manual guess for l=12 component" +{ + : :: +} 0.0 + +REAL ahf_l14_guess "Manual guess for l=14 component" +{ + : :: +} 0.0 + +REAL ahf_l16_guess "Manual guess for l=16 component" +{ + : :: +} 0.0 + +REAL ahf_l18_guess "Manual guess for l=18 component" +{ + : :: +} 0.0 + +REAL ahf_r0 "Radius of initial sphere (0. forces largest sphere)" +{ + : :: +} 0.0 + +REAL ahf_r0_0 "Radius of first initial sphere for find3 (largest for 0.)" +{ + : :: +} 0.0 + +REAL ahf_r0_1 "Radius of second initial sphere for find3 (largest for 0.)" +{ + : :: +} 0.0 + +REAL ahf_r0_2 "Radius of third initial sphere for find3 (largest for 0.)" +{ + : :: +} 0.0 + +REAL ahf_xc "x-coordinate of center of expansion" +{ + : :: +} 0.0 + +REAL ahf_yc "y-coordinate of center of expansion" +{ + : :: +} 0.0 + +REAL ahf_zc "z-coordinate of center of expansion" +{ + : :: +} 0.0 + +REAL ahf_xc_0 "x-coordinate of center of expansion for first surface with find3" +{ + : :: +} 0.0 + +REAL ahf_yc_0 "y-coordinate of center of expansion for first surface with find3" +{ + : :: +} 0.0 + +REAL ahf_zc_0 "z-coordinate of center of expansion for first surface with find3" +{ + : :: +} 0.0 + +REAL ahf_xc_1 "x-coordinate of center of expansion for second surface with find3" +{ + : :: +} 0.0 + +REAL ahf_yc_1 "y-coordinate of center of expansion for second surface with find3" +{ + : :: +} 0.0 + +REAL ahf_zc_1 "z-coordinate of center of expansion for second surface with find3" +{ + : :: +} 0.0 + +REAL ahf_xc_2 "x-coordinate of center of expansion for third surface with find3" +{ + : :: +} 0.0 + +REAL ahf_yc_2 "y-coordinate of center of expansion for third surface with find3" +{ + : :: +} 0.0 + +REAL ahf_zc_2 "z-coordinate of center of expansion for third surface with find3" +{ + : :: +} 0.0 + +#Parameters for flow algorithm. + +LOGICAL ahf_flow "Use flow instead of minimization?" +{ +} "no" + +INT ahf_flowiter "Maximum number of iterations for flow" +{ + : :: +} 200 + +REAL ahf_flowa "alpha parameter for flow" +{ + : :: +} 0.01 + +REAL ahf_flowb "beta parameter for flow" +{ + : :: +} 0.5 + +REAL ahf_flowh "Weight of H flow" +{ + : :: +} 0.0 + +REAL ahf_flowc "Weight of C flow" +{ + : :: +} 1.0 + +REAl ahf_flown "Weight of N flow" +{ + : :: +} 0.0 + +REAL ahf_flowtol "Tolerance for flow" +{ + : :: +} 0.0001 + +REAL ahf_maxchange "max rel diff between 1 big and 2 small steps" +{ + : :: +} 0.1 + +REAL ahf_minchange "min rel diff between 1 big and 2 small steps" +{ + : :: +} 0.01 + +#Parameters for output. + +#IMPORTANT: 3D output for grid functions is still not implemented! + +LOGICAL ahf_logfile "Write log file?" +{ +} "no" + +LOGICAL ahf_verbose "Print messages to screen?" +{ +} "yes" + +LOGICAL ahf_veryverbose "Print messages at each iteration step to screen?" +{ +} "no" + +LOGICAL ahf_guessverbose "Print info on initial guell points?" +{ +} "no" + +LOGICAL ahf_2Doutput "2D output of grid functions?" +{ +} "yes" + +LOGICAL ahf_3Doutput "3D output of grid functions? (Not yet implemented!)" +{ +} "no" + +LOGICAL ahf_areamap "Find area map?" +{ +} "no" + +LOGICAL ahf_gaussout "Output guassian curvature of horizon?" +{ +} "yes" + + + +#Parameters for mask. +#IMPORTANT: By default the mask is always off! + +LOGICAL ahf_mask "Use mask for definite horizons?" +{ +} "no" + +LOGICAL ahf_weakmask "Use mask for possible horizons?" +{ +} "no" + +REAL ahf_maskshrink "Shrink factor for mask" +{ + : :: +} -2.0 + +LOGICAL ah_persists "Does the AH Finder array and mask stay around" +{ +} "no" + + + + +############################################################################# +### import IOUtil parameters +############################################################################# +shares: IO + +#################### +# Output directories +#################### +EXTEND STRING IO_outdir "Name of IO output directory" +{ +} "." +EXTEND STRING IO_outdir0D "Name of IO 0D output directory, overrides IO_outdir" +{ +} "IO_outdir" +EXTEND STRING IO_outdir1D "Name of IO 1D output directory, overrides IO_outdir" +{ +} "IO_outdir" +EXTEND STRING IO_outdir2D "Name of IO 2D output directory, overrides IO_outdir" +{ +} "IO_outdir" +EXTEND STRING IO_outdir3D "Name of IO 3D output directory, overrides IO_outdir" +{ +} "IO_outdir" + + +######################## +### Get symmetries ### +######################## + +shares: grid + +EXTENDS KEYWORD symmetry "" +{ +} diff --git a/schedule.ccl b/schedule.ccl new file mode 100644 index 0000000..635ec9b --- /dev/null +++ b/schedule.ccl @@ -0,0 +1,31 @@ +# Schedule definitions for thorn AHFinder +#c/*@@ +#c @date July 1999 +#c @author Lars Nerger +#c @@*/ + + +if (ahf_active && ah_persists) +{ + STORAGE: ahfindergrid,ahfmask + COMMUNICATION: ahfindergrid,ahfmask + + schedule ahfinder at CCTK_POSTSTEP + { + LANG: Fortran + STORAGE: ahfgradient,ahfinder_gauss,find3grid + COMMUNICATION: ahfgradient,ahfinder_gauss,find3grid + } "Call apparent horizon finder with persisting grid" +} + +else if (ahf_active) +{ + schedule ahfinder at CCTK_POSTSTEP + { + LANG: Fortran + STORAGE:ahfgradient,ahfinder_gauss,find3grid,ahfindergrid,ahfmask + COMMUNICATION: ahfgradient,ahfinder_gauss,find3grid,ahfindergrid,ahfmask + } "Call apparent horizon finder" +} + + diff --git a/src/AHFinder.F b/src/AHFinder.F new file mode 100644 index 0000000..73992b1 --- /dev/null +++ b/src/AHFinder.F @@ -0,0 +1,1524 @@ +c/*@@ +c @file AHFinder.F +c @date April 1998 +c @author Miguel Alcubierre +c @desc +c Find an apparent horizon. +c @enddesc +c@@*/ + +c Note that including cactus.h will also include AHFinder.h +#include "cctk.h" +#include "cctk_parameters.h" +#include "cctk_arguments.h" + + subroutine AHFinder(CCTK_FARGUMENTS) + + use AHFinder_dat + + implicit none + + + DECLARE_CCTK_FARGUMENTS + DECLARE_CCTK_PARAMETERS + + logical status,horizon,gauss + + + integer handle,x_index,y_index,z_index,ierror + CCTK_REAL maxval(3),minval(3) + + CCTK_INT ahfso,ahfafter + integer i,j,l,m + integer NN + integer mtype + CCTK_INT getoutpfx + INTEGER CCTK_Equals + + CCTK_REAL atime + CCTK_REAL zero,half,one + CCTK_REAL xmn,ymn,zmn,xmx,ymx,zmx,rmx,r0 + CCTK_REAL intexp_h,intexp2_h,intarea_h + CCTK_REAL circ_eq,meri_p1,meri_p2 + CCTK_REAL aux,expin,expout,small1,small2 + CCTK_REAL c0_temp(18) + + character*200 almf,logf,areaf,massf,circeqf,merip1f,merip2f + + save ahfafter,ahfso,atime + + +! Description of variables: +! +! status Surface found? +! +! horizon Horizon found? +! +! gauss Ouput Gaussian curvature? +! +! ahfso How often we want to find a horizon? +! ahfafter After how many iterations start looking for horizons? +! +! i,l,m Counters. +! +! NN Total number of terms in expansion. +! +! xmn Minimum value of x over the grid. +! xmx Maximum value of x over the grid. +! +! ymn Minimum value of y over the grid. +! ymx Maximum value of y over the grid. +! +! zmn Minimum value of z over the grid. +! zmx Maximum value of z over the grid. +! +! rmx Radius of largest sphere that fits in the grid. +! +! circ_eq Equatorial circumference. +! meri_p1 Length of meridian at phi=0. +! meri_p2 Length of meridian at phi=pi/2. +! +! mtype Type of surface. +! 1. Integral of H^2 is positive out and negative in. +! 2. Integral of H^2 is negative out and positive in. +! 3. Integral of H^2 is negative out and in. +! 4. Integral of H^2 is positive out and in. + + +! ************************** +! *** DEFINE NUMBERS *** +! ************************** + + zero = 0.0D0 + half = 0.5D0 + one = 1.0D0 + + dx = cctk_delta_space(1) + dy = cctk_delta_space(2) + dz = cctk_delta_space(3) + +! ******************************** +! *** INITIALIZE {status} *** +! ******************************** + + status = .false. + +! ********************************* +! *** INITIALIZE FIRSTCALLS *** +! ********************************* + + firstleg = .true. + firstint = .true. + + +! *********************************************** +! *** CHECK IF WE WANT TO FIND 3 HORIZONS *** +! *********************************************** + +! find3 = contains("ahf_find3","yes").ne.0 + if (ahf_find3.ne.0) then + find3 = .true. + else + find3 = .false. + end if + + +! ********************************************* +! *** CHECK IF WE WANT TO FIND A HORIZON *** +! ********************************************* + +! See if we want to find a horizon at this time. + + if (ahf_ncall.eq.0) then + +! atime = getr8("ahf_findaftertime") +! ahfafter = geti("ahf_findafter") +! ahfso = geti("ahf_findevery") + + atime = ahf_findaftertime + ahfafter = ahf_findafter + ahfso = ahf_findevery + + if (atime.le.zero) then + if (cctk_iteration.lt.ahfafter) return + else + if (cctk_time.lt.atime) return + ahfafter = cctk_iteration + end if + + end if + + if (mod(cctk_iteration-ahfafter,ahfso).ne.0) return + +! BoxInBox may turn of AHF this way during runtime (BB) + +! if (contains("ahf_active","no").ne.0) then + if (ahf_active.eq.0) then + return + endif + +! Sorry for this hack, but the call order for BoxInBox is screwy in +! cactus 3.2 +#ifdef THORN_BOXINBOX +! if ((contains("boxinbox","yes").ne.0).and.(ahf_ncall.eq.0)) then + if ((boxinbox.ne.0).and.(ahf_ncall.eq.0)) then + return; + endif +#endif + + +! ******************************************* +! *** ADD TO COUNT OF CALLS TO FINDER *** +! ******************************************* + + 1 continue + +! Say what we are doing. + + if (find3) then + if (mfind.eq.0) then + write(*,*) + write(*,*) 'AHFinder: Looking for 3 apparent horizons' + write(*,*) 'AHFinder: Searching for first horizon' + else if (mfind.eq.1) then + write(*,*) + write(*,*) 'AHFinder: Searching for second horizon' + else if (mfind.eq.2) then + write(*,*) + write(*,*) 'AHFinder: Searching for third horizon' + end if + else + write(*,*) + write(*,*) 'AHFinder: Looking for horizon' + end if + +! Add to count of calls to finder. Notice that this counter +! was initialized to zero, so it must be increased now to 1. +! If we look for 3 horizons, it should be done before looking +! for the first one. + + if (find3) then + if (mfind.eq.0) then + ahf_ncall = ahf_ncall + 1 + end if + else + ahf_ncall = ahf_ncall + 1 + end if + + +! ******************************************************* +! *** FIND OUT IF WE WANT TO FIND TRAPPED SURFACE *** +! ******************************************************* + +! find_trapped_surface = contains("ahf_trapped_surface","yes").ne.0 + if (ahf_trapped_surface.ne.0) then + find_trapped_surface = .true. + else + find_trapped_surface = .false. + end if + +! trapped_surface_delta = getr8("trapped_surface_delta") + + +! ******************************************** +! *** GET THE NAME OF OUTPUT DIRECTORY *** +! ******************************************** + + if (.false.) then +! nfile = getoutpfx(GH,filestr) + filestr = IO_outdir + + + if (find3) then + + if (mfind.eq.0) then + + logf = filestr(1:nfile)//"ahf_logfile_0" + almf = filestr(1:nfile)//"ahf_coeff_0.alm" + + areaf = filestr(1:nfile)//"ahf_area_0.tl" + massf = filestr(1:nfile)//"ahf_mass_0.tl" + + circeqf = filestr(1:nfile)//"ahf_circ_eq_0.tl" + merip1f = filestr(1:nfile)//"ahf_meri_p1_0.tl" + merip2f = filestr(1:nfile)//"ahf_meri_p2_0.tl" + + else if (mfind.eq.1) then + + logf = filestr(1:nfile)//"ahf_logfile_1" + almf = filestr(1:nfile)//"ahf_coeff_1.alm" + + areaf = filestr(1:nfile)//"ahf_area_1.tl" + massf = filestr(1:nfile)//"ahf_mass_1.tl" + + circeqf = filestr(1:nfile)//"ahf_circ_eq_1.tl" + merip1f = filestr(1:nfile)//"ahf_meri_p1_1.tl" + merip2f = filestr(1:nfile)//"ahf_meri_p2_1.tl" + + else + + logf = filestr(1:nfile)//"ahf_logfile_2" + almf = filestr(1:nfile)//"ahf_coeff_2.alm" + + areaf = filestr(1:nfile)//"ahf_area_2.tl" + massf = filestr(1:nfile)//"ahf_mass_2.tl" + + circeqf = filestr(1:nfile)//"ahf_circ_eq_2.tl" + merip1f = filestr(1:nfile)//"ahf_meri_p1_2.tl" + merip2f = filestr(1:nfile)//"ahf_meri_p2_2.tl" + + end if + + else + + logf = filestr(1:nfile)//"ahf_logfile" + almf = filestr(1:nfile)//"ahf_coeff.alm" + + areaf = filestr(1:nfile)//"ahf_area.tl" + massf = filestr(1:nfile)//"ahf_mass.tl" + + circeqf = filestr(1:nfile)//"ahf_circ_eq.tl" + merip1f = filestr(1:nfile)//"ahf_meri_p1.tl" + merip2f = filestr(1:nfile)//"ahf_meri_p2.tl" + + end if + end if + +! ************************************************** +! *** TRANSFORM CONFORMAL TO PHYSICAL METRIC *** +! ************************************************** + +! call conf2phys(CONF2PHYSVARS_ARGS) + + +! *************************************** +! *** FIND OUT THE TYPE OF OUTPUT *** +! *************************************** + if (.false.) then +! logfile = contains("ahf_logfile","yes").ne.0 + if (ahf_logfile.ne.0) then + logfile = .true. + else + logfile = .false. + end if + +! verbose = contains("ahf_verbose","yes").ne.0 + if (ahf_verbose.ne.0) then + verbose = .true. + else + verbose = .false. + end if + +! veryver = contains("ahf_veryverbose","yes").ne.0 + if (ahf_veryverbose.ne.0) then + veryver = .true. + else + veryver = .false. + end if + +! guessverbose = contains("ahf_guessverbose","yes").ne.0 + if (ahf_guessverbose.ne.0) then + guessverbose = .true. + else + guessverbose = .false. + end if + + +! if (logfile.and.(myproc.eq.0)) then + if (logfile) then +#ifdef MPI + if (myproc.eq.0) then +#endif + open(1,file=logf,form='formatted',status='replace') + write(1,*) + write(1,*) 'LOG FILE FOR MINIMIZATION APPARENT HORIZON FINDER' + write(1,*) + write(1,*) 'Note: During an evolution, This file will only' + write(1,*) 'contain information about the last time that the' + write(1,*) 'finder was called.' + write(1,*) + close(1) +#ifdef MPI + end if +#endif + end if + end if + +! ******************************************* +! *** FIND {xmn,xmx,ymn,ymx,zmn,zmx} *** +! ******************************************* + + call CCTK_ReductionHandle(handle,"minimum") + call CCTK_GetCoordIndex(x_index,"x") + call CCTK_GetCoordIndex(y_index,"y") + call CCTK_GetCoordIndex(z_index,"z") + + call CCTK_Reduce(cctkGH,ierror,-1,handle,3, + . CCTK_VARIABLE_REAL, + . minval,3,x_index,y_index,z_index) + +! xmn = gf_min(x) +! ymn = gf_min(y) +! zmn = gf_min(z) + + xmn = minval(1) + ymn = minval(2) + zmn = minval(3) + + call CCTK_ReductionHandle(handle,"maximum") + call CCTK_GetCoordIndex(x_index,"x") + call CCTK_GetCoordIndex(y_index,"y") + call CCTK_GetCoordIndex(z_index,"z") + + call CCTK_Reduce(cctkGH,ierror,-1,handle,3, + . CCTK_VARIABLE_REAL, + . maxval,3,x_index,y_index,z_index) + +! xmx = gf_max(x) +! ymx = gf_max(y) +! zmx = gf_max(z) + + xmx = maxval(1) + ymx = maxval(2) + zmx = maxval(3) + + +! ********************************** +! *** INITIALIZE {xc,yc,zc} *** +! ********************************** + + if (find3) then + if (mfind.eq.0) then +! xc = getr8("ahf_xc_0") +! yc = getr8("ahf_yc_0") +! zc = getr8("ahf_zc_0") + xc = ahf_xc_0 + yc = ahf_yc_0 + zc = ahf_zc_0 + else if (mfind.eq.1) then +! xc = getr8("ahf_xc_1") +! yc = getr8("ahf_yc_1") +! zc = getr8("ahf_zc_1") + xc = ahf_xc_1 + yc = ahf_yc_1 + zc = ahf_zc_1 + else +! xc = getr8("ahf_xc_2") +! yc = getr8("ahf_yc_2") +! zc = getr8("ahf_zc_2") + xc = ahf_xc_2 + yc = ahf_yc_2 + zc = ahf_zc_2 + end if + else +! xc = getr8("ahf_xc") +! yc = getr8("ahf_yc") +! zc = getr8("ahf_zc") + xc = ahf_xc + yc = ahf_yc + zc = ahf_zc + end if + + if ((xc.gt.xmx).or.(xc.lt.xmn)) then + write(*,*) + write(*,*) 'xc is out of the grid, reseting it to zero ...' + write(*,*) + xc = zero + end if + + if ((yc.gt.ymx).or.(yc.lt.ymn)) then + write(*,*) + write(*,*) 'yc is out of the grid, reseting it to zero ...' + write(*,*) + yc = zero + end if + + if ((yc.gt.ymx).or.(yc.lt.ymn)) then + write(*,*) + write(*,*) 'zc is out of the grid, reseting it to zero ...' + write(*,*) + zc = zero + end if + + +! ************************** +! *** FIND SYMMETRIES *** +! ************************** + +! Find out if we want to move the center. + +! offset = contains("ahf_offset","yes").ne.0 + if (ahf_offset.ne.0) then + offset = .true. + else + offset = .false. + end if +! wander = contains("ahf_wander","yes").ne.0 + if (ahf_wander.ne.0) then + wander = .true. + else + wander = .false. + end if + + if ((xc.ne.zero).or.(yc.ne.zero).or.(zc.ne.zero).or. + . (wander)) offset = .true. + +! Find out if we are in non-axisymmetric mode. + + if (ahf_phi.ne.0) then + nonaxi = .true. + else + nonaxi = .false. + end if +! nonaxi = contains("ahf_phi","yes").ne.0 + +! Find out reflection symmetries. + +! if ((contains("ahf_octant","yes").ne.0).or. +! . (contains("ahf_octant","high").ne.0)) then + if ((CCTK_Equals(ahf_octant,"yes").eq.1).or. + . (CCTK_Equals(ahf_octant,"high").eq.1)) then + + refx = .true. + refy = .true. + refz = .true. + else + if (ahf_refx.ne.0) then + refx = .true. + else + refx = .false. + end if + if (ahf_refy.ne.0) then + refy = .true. + else + refy = .false. + end if + if (ahf_refz.ne.0) then + refz = .true. + else + refz = .false. + end if + +! refx = contains("ahf_refx","yes").ne.0 +! refy = contains("ahf_refy","yes").ne.0 +! refz = contains("ahf_refz","yes").ne.0 + end if + +! if (contains("ahf_phi","no").ne.0) then + if (ahf_phi.eq.0) then + refx = .true. + refy = .true. + end if + +! Force grid symmetries. + +! if (contains("grid","octant").ne.0) then + if (CCTK_Equals(symmetry,"octant").eq.1) then + if (xc.eq.zero) refx = .true. + if (yc.eq.zero) refy = .true. + if (zc.eq.zero) refz = .true. +! else if (contains("grid","quadrant").ne.0) then + else if (CCTK_Equals(symmetry,"quadrant").eq.1) then + if (xc.eq.zero) refx = .true. + if (yc.eq.zero) refy = .true. + end if + +! Find {stepx,stepy,stepz}. + + stepx = 0 + stepy = 0 + stepz = 0 + + if (refx) stepx=1 + if (refy) stepy=1 + if (refz) stepz=1 + +! if (contains("ahf_octant","high").ne.0) stepx=3 + if (CCTK_Equals(ahf_octant,"high").eq.1) stepx=3 + +! Check if we are using Cartoon_2d, in which case +! we override the value of "ahf_phi" since we have +! axisymmetry. + + cartoon = .false. + +#ifdef THORN_cartoon_2d +! cartoon = contains("cartoon_active","yes").ne.0 + if (cartoon_active.ne.0) then + cartoon = .true. + else + cartoon = .false. + end if +#endif + + if (cartoon) then + nonaxi = .false. + refx = .true. + refy = .true. + end if + + +! ************************ +! *** FIND {lmax} *** +! ************************ + +! lmax = geti("ahf_lmax") + lmax = ahf_lmax + + if (lmax.ge.20) then + write(*,*) + write(*,*) 'ahf_lmax must be smaller than 20.' + write(*,*) 'STOPPED IN AHFinder.F' + write(*,*) + STOP + end if + + +! ******************************* +! *** FIND {ntheta,nphi} *** +! ******************************* + +! ntheta = geti("ahf_ntheta") +! nphi = geti("ahf_nphi") + + ntheta = ahf_ntheta + nphi = ahf_nphi + +! For cartoon, I use only 2 points in the phi direction. +! I could use one, but using 2 minimizes the changes in +! the integration routine. + + if (cartoon) then + nphi = 1 + end if + + +! *************************************************** +! *** ALLOCATE STORAGE FOR COEFFICIENT ARRAYS *** +! *************************************************** + + if ((ahf_ncall.eq.1).and.(mfind.eq.0)) then + allocate(c0(0:lmax),c0_old(0:lmax)) + allocate(cc(lmax,lmax),cc_old(lmax,lmax)) + allocate(cs(lmax,lmax),cs_old(lmax,lmax)) + end if + + +! ******************************************************** +! *** FIND TOTAL NUMBER OF TERMS IN EXPANSION {NN} *** +! ******************************************************** + + NN = 0 + +! Count {xc,yc,zc}. + + if (wander) then + if (.not.refx) NN = NN + 1 + if (.not.refy) NN = NN + 1 + if (.not.refz) NN = NN + 1 + end if + +! Count c0(l). + + do l=0,lmax,1+stepz + NN = NN + 1 + end do + + if (nonaxi) then + +! Count cc(l,m). + + do l=1,lmax + do m=1+stepx,l,1+stepx + if (stepz*mod(l-m,2).eq.0) NN = NN +1 + end do + end do + +! Count cs(l,m). + + if (.not.refy) then + do l=1,lmax + do m=1,l,1+stepx + if (stepz*mod(l-m,2).eq.0) NN = NN + 1 + end do + end do + end if + + end if + + +! ***************************************** +! *** FIND RADIUS OF INITIAL SPHERE *** +! ***************************************** + +! Find largest sphere that fits into the grid, taking +! into account the position of the center. Then multiply +! this number with 0.8 to leave a safety margin at the edge +! of the grid. + +! Cartoon. + + if (cartoon) then + + xc = zero + yc = zero + rmx = min(xmx,zmx-zc,zc-zmn) + +! Octant. + +! else if (contains("grid","octant").ne.0) then + else if (CCTK_Equals(symmetry,"octant").eq.1) then + if ((xc.eq.zero).and.(yc.eq.zero).and.(zc.eq.zero)) then + rmx = min(xmx,ymx,zmx) + else if ((xc.eq.zero).and.(yc.eq.zero)) then + rmx = min(xmx,ymx,zmx-zc,zc) + else if ((xc.eq.zero).and.(zc.eq.zero)) then + rmx = min(xmx,ymx-yc,zmx,yc) + else if ((yc.eq.zero).and.(zc.eq.zero)) then + rmx = min(xmx-xc,ymx,zmx,xc) + else if (xc.eq.zero) then + rmx = min(xmx,ymx-yc,zmx-zc,yc,zc) + else if (yc.eq.zero) then + rmx = min(xmx-xc,ymx,zmx-zc,xc,zc) + else if (zc.eq.zero) then + rmx = min(xmx-xc,ymx-yc,zmx,xc,yc) + else + rmx = min(xmx-xc,ymx-yc,zmx-zc,xc,yc,zc) + end if + +! Quadrant. + +! else if (contains("grid","quadrant").ne.0) then + else if (CCTK_Equals(symmetry,"quadrant").eq.1) then + if ((xc.eq.zero).and.(yc.eq.zero)) then + rmx = min(xmx,ymx,zmx-zc,zc-zmn) + else if (xc.eq.zero) then + rmx = min(xmx,ymx-yc,zmx-zc,yc,zc-zmn) + else if (yc.eq.zero) then + rmx = min(xmx-xc,ymx,zmx-zc,xc,zc-zmn) + else + rmx = min(xmx-xc,ymx-yc,zmx-zc,xc,yc,zc-zmn) + end if + +! Full. + + else + + rmx = min(xmx-xc,ymx-yc,zmx-zc,xc-xmn,yc-ymn,zc-zmn) + + end if + + rmx = 0.8*rmx + +! Check if a special radius for the initial sphere +! has been forced by the parameter file. + + if (find3) then + if (mfind.eq.0) then +! r0 = getr8("ahf_r0_0") + r0 = ahf_r0_0 + else if (mfind.eq.1) then +! r0 = getr8("ahf_r0_1") + r0 = ahf_r0_1 + else +! r0 = getr8("ahf_r0_2") + r0 = ahf_r0_2 + end if + else +! r0 = getr8("ahf_r0") + r0 = ahf_r0 + end if + + if (r0.ne.zero) then + if (r0.le.rmx) then + rmx = r0 + else + write (*,*) + write (*,*) 'ahf_r0 is too large, rescaling ...' + write (*,*) + end if + end if + + +! ************************************************* +! *** GET INITIAL GUESS AND LOOK FOR HORIZON *** +! ************************************************* + +! initializing flags + + if (ahf_sloppyguess.ne.0) then + sloppy = .true. + else + sloppy = .false. + end if + +! sloppy = contains("ahf_sloppyguess","yes").ne.0 + if (cartoon) sloppy = .true. + + if (ahf_flow.ne.0) then + flow = .true. + else + flow = .false. + end if +! flow = contains("ahf_flow","yes").ne.0 + if (ahf_inner.ne.0) then + inner = .true. + else + inner = .false. + end if +! inner = contains("ahf_inner","yes").ne.0 + if (ahf_guessold.ne.0) then + guessold = .true. + else + guessold = .false. + end if +! guessold = contains("ahf_guessold","yes").ne.0 + if (ahf_guess_absmin.ne.0) then + guess_absmin = .true. + else + guess_absmin = .false. + end if +! guess_absmin = contains("ahf_guess_absmin","yes").ne.0 + if (ahf_manual_guess.ne.0) then + manual_guess = .true. + else + manual_guess = .false. + end if +! manual_guess = contains("ahf_manual_guess","yes").ne.0 + if (ahf_minarea.ne.0) then + minarea = .true. + else + minarea = .false. + end if +! minarea = contains("ahf_minarea","yes").ne.0 + +! get initial guess and jump into finder routines + + if ((.not.flow).and.guessold.and.(ahf_ncall.gt.1)) then + call AHFinder_initguess(CCTK_FARGUMENTS,rmx) + call AHFinder_min(CCTK_FARGUMENTS,NN,status,logf) + else if ((.not.flow).and.manual_guess) then +! if (contains("grid","octant").eq.0) then + if (CCTK_Equals(symmetry,"octant").eq.0) then + write(*,*) 'AHFinder: if you want to use ahf_manual_guess' + write(*,*) ' without grid=octant, you will have to make' + write(*,*) ' it work!' + STOP + endif +! c0(0) = getr8("ahf_l0_guess") +! c0_temp(2) = getr8("ahf_l2_guess") +! c0_temp(4) = getr8("ahf_l4_guess") +! c0_temp(6) = getr8("ahf_l6_guess") +! c0_temp(8) = getr8("ahf_l8_guess") +! c0_temp(10) = getr8("ahf_l10_guess") +! c0_temp(12) = getr8("ahf_l12_guess") +! c0_temp(14) = getr8("ahf_l14_guess") +! c0_temp(16) = getr8("ahf_l16_guess") +! c0_temp(18) = getr8("ahf_l18_guess") + c0(0) = ahf_l0_guess + c0_temp(2) = ahf_l2_guess + c0_temp(4) = ahf_l4_guess + c0_temp(6) = ahf_l6_guess + c0_temp(8) = ahf_l8_guess + c0_temp(10) = ahf_l10_guess + c0_temp(12) = ahf_l12_guess + c0_temp(14) = ahf_l14_guess + c0_temp(16) = ahf_l16_guess + c0_temp(18) = ahf_l18_guess + do i = 2,lmax,2 + c0(i) = c0_temp(i) + enddo + call AHFinder_initguess(CCTK_FARGUMENTS,rmx) + call AHFinder_min(CCTK_FARGUMENTS,NN,status,logf) + else if (.not.flow) then + call AHFinder_initguess(CCTK_FARGUMENTS,rmx) + call AHFinder_min(CCTK_FARGUMENTS,NN,status,logf) + else + if ((ahf_ncall.eq.1).and.(mfind.eq.0)) then + allocate(hflow0(0:lmax),cflow0(0:lmax),nflow0(0:lmax)) + allocate(hflowc(lmax,lmax),cflowc(lmax,lmax),nflowc(lmax,lmax)) + allocate(hflows(lmax,lmax),cflows(lmax,lmax),nflows(lmax,lmax)) + end if + call AHFinder_flow(CCTK_FARGUMENTS,rmx,status,logf) + flow = .false. + end if + + +! *************************************** +! *** WRITE MESSAGES AND LOG FILE *** +! *************************************** + + if (status) then + +! Check if the expansion on the surface is small enough +! for a horizon. How small is small enough is of course +! pretty arbitrary, so to reduce the subjectivity I use +! three separate criteria: +! +! 1) The integral of the square of the expansion at the surface +! must be smaller than a given threshold (small1). +! +! 2) The mean value of the expansion on the surface must be +! smaller in absolute value than its standard deviation. +! This means that the mean value of the expansion is +! consistent with zero. +! +! 3) The integral of the square of the expansion at the surface +! must be smaller by a factor of `small2' than the same +! integral slightly outside and slightly inside. + + horizon = .true. + + small1 = 1.0D-2 + small2 = 0.5D0 + +! Surface found. + + call AHFinder_fun(CCTK_FARGUMENTS) + call AHFinder_exp(CCTK_FARGUMENTS) + call AHFinder_int(CCTK_FARGUMENTS) + +! gauss = contains("ahf_gaussout","yes").ne.0 + + if (ahf_gaussout.ne.0) then + gauss = .true. + else + gauss = .false. + end if + + if (.false.) then + if (gauss) then + call AHFinder_gau(CCTK_FARGUMENTS,circ_eq,meri_p1,meri_p2) + end if + end if + + intexp_h = intexp + intexp2_h = intexp2 + intarea_h = intarea + + if (intexp2_h.gt.small1) horizon = .false. + +! Find standard deviation: sqrt(<exp^2> - <exp>^2) + + aux = dsqrt(intexp2_h/intarea_h - (intexp_h/intarea_h)**2) + + if (dabs(intexp_h/intarea_h)-aux.gt.zero) horizon = .false. + +! Surface slightly inside + + aux = c0(0) + + c0(0) = aux - half*min(dx,dy,dz) + call AHFinder_int(CCTK_FARGUMENTS) + + expin = intexp + + if (intexp2_h.gt.small2*intexp2) horizon = .false. + + if (veryver) then + write(*,*) + write(*,*) 'Surface slightly inside the minimum surface: ' + write(*,"(A8,ES14.6)") ' H : ',intexp + write(*,"(A8,ES14.6)") ' H^2: ',intexp2 + write(*,"(A13,ES14.6)") ' H/area : ',intexp/intarea + write(*,"(A13,ES14.6)") ' H^2/area: ',intexp2/intarea + write(*,*) 'number of interpolated points: ', + . inside_min_count + write(*,*) ' number of those that are negative: ', + . inside_min_neg_count + endif + +! Surface slightly outside. + + c0(0) = aux + half*min(dx,dy,dz) + call AHFinder_int(CCTK_FARGUMENTS) + + expout = intexp + + if (intexp2_h.gt.small2*intexp2) horizon = .false. + + if (veryver) then + write(*,*) + write(*,*) 'Surface slightly outside the minimum surface: ' + write(*,"(A8,ES14.6)") ' H : ',intexp + write(*,"(A8,ES14.6)") ' H^2: ',intexp2 + write(*,"(A13,ES14.6)") ' H/area : ',intexp/intarea + write(*,"(A13,ES14.6)") ' H^2/area: ',intexp2/intarea + endif + + c0(0) = aux + +! What kind of horizon? + + if ((expin.le.zero).and.(expout.ge.zero)) then + mtype = 1 + else if ((expin.ge.zero).and.(expout.le.zero)) then + mtype = 2 + else if ((expin.le.zero).and.(expout.le.zero)) then + mtype = 3 + else + mtype = 4 + end if + +! Messages to screen. +#ifdef MPI + if (myproc.eq.0) then +#endif +! Open logfile. + + if (logfile) then + open(1,file=logf,form='formatted', + . status='old',position='append') + end if + +! Surface information. Notice that here I define +! the surface `mass' as: +! +! M = sqrt( A / (16 pi) ) = 0.141047396 sqrt(A) + + if (intarea_h.ne.zero) then + aux = one/intarea_h + else + aux = one + end if + + write(*,*) + write(*,*) 'AHFinder: Surface found.' + write(*,*) + write(*,"(A21,ES14.6)") ' Surface area =',intarea_h + write(*,"(A21,ES14.6)") ' Surface mass =', + . 0.141047396D0*dsqrt(intarea_h) + write(*,"(A21,ES14.6)") ' Mean value of H =',aux*intexp_h + write(*,"(A21,ES14.6)") ' Mean value of H^2 =',aux*intexp2_h + write(*,"(A21,ES14.6)") ' Standard deviation =', + . dsqrt(aux*intexp2_h - (aux*intexp_h)**2) + write(*,*) + write(*,"(A11,ES14.6)") ' circ_eq =',circ_eq + write(*,"(A11,ES14.6)") ' meri_p1 =',meri_p1 + write(*,"(A11,ES14.6)") ' meri_p2 =',meri_p2 + write(*,*) + + if (.false.) then + if (logfile) then + write(1,*) + write(1,*) 'AHFinder: Surface found.' + write(1,*) + write(1,"(A21,ES14.6)") ' Surface area =',intarea_h + write(1,"(A21,ES14.6)") ' Surface mass =', + . 0.141047396D0*dsqrt(intarea_h) + write(1,"(A21,ES14.6)") ' Mean value of H =',aux*intexp_h + write(1,"(A21,ES14.6)") ' Mean value of H^2 =',aux*intexp2_h + write(1,"(A21,ES14.6)") ' Standard deviation =', + . dsqrt(aux*intexp2_h - (aux*intexp_h)**2) + write(1,*) + write(1,"(A11,ES14.6)") ' circ_eq =',circ_eq + write(1,"(A11,ES14.6)") ' meri_p1 =',meri_p1 + write(1,"(A11,ES14.6)") ' meri_p2 =',meri_p2 + write(1,*) + end if + end if + + if (mtype.eq.1) then + if (horizon) then + write(*,*) 'The surface seems to be an outer horizon.' + if (logfile) then + write(1,*) 'The surface seems to be an outer horizon.' + end if + else + write(*,*) 'The surface is probably an outer horizon, but' + write(*,*) 'the mean of the square of the expansion is' + write(*,*) 'not small enough. Try a larger lmax, or a' + write(*,*) 'smaller grid size.' + if (logfile) then + write(1,*) 'The surface is probably an outer horizon, but' + write(1,*) 'the mean of the square of the expansion is' + write(1,*) 'not small enough. Try a larger lmax, or a' + write(1,*) 'smaller grid size.' + end if + end if + else if (mtype.eq.2) then + if (horizon) then + write(*,*) 'The surface seems to be an inner horizon.' + if (logfile) then + write(1,*) 'The surface seems to be an inner horizon.' + end if + else + write(*,*) 'The surface is probably an inner horizon, but' + write(*,*) 'the mean of the square of the expansion is' + write(*,*) 'not small enough. Try a larger lmax, or a' + write(*,*) 'smaller grid size.' + if (logfile) then + write(1,*) 'The surface is probably an inner horizon, but' + write(1,*) 'the mean of the square of the expansion is' + write(1,*) 'not small enough. Try a larger lmax, or a' + write(1,*) 'smaller grid size.' + end if + end if + else if (mtype.eq.3)then + if (horizon) then + write(*,*) 'The surface seems to be a marginal horizon' + write(*,*) '(the integral of the expansion is negative' + write(*,*) 'on both sides.)' + if (logfile) then + write(1,*) 'The surface seems to be a marginal horizon' + write(1,*) '(the integral of the expansion is negative' + write(1,*) 'on both sides.)' + end if + else + write(*,*) 'The surface seems to be a trapped surface.' + if (logfile) then + write(1,*) 'The surface seems to be a trapped surface.' + end if + end if + else if (mtype.eq.4)then + if (horizon) then + write(*,*) 'The surface seems to be a marginal horizon' + write(*,*) '(the integral of the expansion is positive' + write(*,*) 'on both sides.)' + if (logfile) then + write(1,*) 'The surface seems to be a marginal horizon' + write(1,*) '(the integral of the expansion is positive' + write(1,*) 'on both sides.)' + end if + else + write(*,*) 'The surface does not seem to be a horizon.' + if (logfile) then + write(1,*) 'The surface does not seem to be a horizon.' + end if + end if + end if + +! Shape. + + if (verbose) then + write(*,*) + write(*,*) 'Shape coefficients:' + + if (offset) then + write(*,*) + write(*,"(A6,ES14.6)") ' xc =',xc + write(*,"(A6,ES14.6)") ' yc =',yc + write(*,"(A6,ES14.6)") ' zc =',zc + end if + + write(*,*) + + write(*,"(A4,I2,A3,ES14.6)") ' c0(',0,') =',c0(0) + + do l=1+stepz,lmax,1+stepz + write(*,"(A4,I2,A3,ES14.6)") ' c0(',l,') =',c0(l) + end do + + end if + + if (logfile) then + write(1,*) + write(1,*) 'Shape coefficients:' + + if (offset) then + write(1,*) + write(1,"(A6,ES14.6)") ' xc =',xc + write(1,"(A6,ES14.6)") ' yc =',yc + write(1,"(A6,ES14.6)") ' zc =',zc + end if + + write(1,*) + + write(1,"(A4,I2,A3,ES14.6)") ' c0(',0,') =',c0(0) + + do l=1+stepz,lmax,1+stepz + write(1,"(A4,I2,A3,ES14.6)") ' c0(',l,') =',c0(l) + end do + + end if + + if (nonaxi) then + + if (.not.refy) then + + if (verbose) then + + write(*,*) + do l=1,lmax + do m=1,l + if (stepz*mod(l-m,2).eq.0) then + write(*,"(A4,I2,A1,I2,A4,ES14.6)") ' cc(', + . l,',',m,') =',cc(l,m) + end if + end do + end do + write(*,*) + do l=1,lmax + do m=1,l + if (stepz*mod(l-m,2).eq.0) then + write(*,"(A4,I2,A1,I2,A4,ES14.6)") ' cs(', + . l,',',m,') =',cs(l,m) + end if + end do + end do + + end if + + if (logfile) then + write(1,*) + + do l=1,lmax + do m=1,l + if (stepz*mod(l-m,2).eq.0) then + write(1,"(A4,I2,A1,I2,A4,ES14.6)") ' cc(', + . l,',',m,') =',cc(l,m) + end if + end do + end do + write(1,*) + do l=1,lmax + do m=1,l + if (stepz*mod(l-m,2).eq.0) then + write(1,"(A4,I2,A1,I2,A4,ES14.6)") ' cs(', + . l,',',m,') =',cs(l,m) + end if + end do + end do + + end if + + 10 format(I4,I4,A2,ES14.6,A1,ES14.6) + + else + + if (verbose) then + + write(*,*) + do l=1,lmax + do m=1,l + if (stepz*mod(l-m,2).eq.0) then + write(*,"(A4,I2,A1,I2,A4,ES14.6)") ' cc(', + . l,',',m,') =',cc(l,m) + end if + end do + end do + + end if + + if (logfile) then + + write(1,*) + do l=1,lmax + do m=1,l + if (stepz*mod(l-m,2).eq.0) then + write(1,"(A4,I2,A1,I2,A4,ES14.6)") ' cc(', + . l,',',m,') =',cc(l,m) + end if + end do + end do + + end if + + 20 format(I4,I4,A2,ES14.6) + + end if + end if + +! Close logfile + + write(*,*) + + if (logfile) then + write(1,*) + write(1,*) 'END OF LOG FILE' + write(1,*) + close(1) + end if +#ifdef MPI + end if +#endif +! No surface found. + + else +#ifdef MPI + if (myproc.eq.0) then +#endif +! Open logfile. + + if (logfile) then + open(1,file=logf,form='formatted', + . status='old',position='append') + end if + +! Write messages. + + write(*,*) + write(*,*) 'AHFinder: No horizon found.' + write(*,*) + + if (logfile) then + write(1,*) + write(1,*) 'AHFinder: No horizon found.' + write(1,*) + end if + +! Close logfile + + write(*,*) + + if (logfile) then + write(1,*) + write(1,*) 'END OF LOG FILE' + write(1,*) + close(1) + end if + + end if +#ifdef MPI + end if +#endif + +! **************************** +! *** WRITE DATA FILES *** +! **************************** +#ifdef MPI + if (myproc.eq.0) then +#endif +c Expansion coefficients. + + if (ahf_ncall.eq.1) then + open(1,file=almf,form='formatted', + . status='replace') + write(1,"(A21)") '# Radial coefficients' + write(1,"(A1)") '#' + write(1,"(A11,I4)") '# Time step',cctk_iteration + write(1,"(A6,ES14.6)") '# Time',cctk_time + write(1,"(A6,I4)") '# Call',ahf_ncall + else + open(1,file=almf,form='formatted', + . status='old',position='append') + write(1,*) + write(1,"(A11,I4)") '# Time step',cctk_iteration + write(1,"(A6,ES14.6)") '# Time',cctk_time + write(1,"(A6,I4)") '# Call',ahf_ncall + end if + + if (status) then + if (mtype.eq.1) then + if (horizon) then + write(1,"(A30,I4)") '# Surface found: Outer horizon' + else + write(1,"(A31,I4)") '# Surface found: Outer horizon?' + end if + else if (mtype.eq.2) then + if (horizon) then + write(1,"(A30,I4)") '# Surface found: Inner horizon' + else + write(1,"(A31,I4)") '# Surface found: Inner horizon?' + end if + else if (mtype.eq.3) then + if (horizon) then + write(1,"(A34,I4)") '# Surface found: Marginal horizon?' + else + write(1,"(A32,I4)") '# Surface found: Trapped surface' + end if + else if (mtype.eq.4) then + if (horizon) then + write(1,"(A34,I4)") '# Surface found: Marginal horizon?' + else + write(1,"(A30,I4)") '# Surface found: Not a horizon' + end if + end if + else + write(1,"(A18,I4)") '# No surface found' + end if + + if (status) then + + write(1,"(A7,ES14.6)") '# xc =',xc + write(1,"(A7,ES14.6)") '# yc =',yc + write(1,"(A7,ES14.6)") '# zc =',zc + + write(1,"(A1)") '#' + write(1,"(A22)") '# a_lm l m' + write(1,"(A1)") '#' + + write(1,"(ES14.6,2I4)") c0(0),0,0 + + do l=1,lmax + do m=l,1,-1 + write(1,"(ES14.6,2I4)") cs(l,m),l,-m + end do + write(1,"(ES14.6,2I4)") c0(l),l,0 + do m=1,l + write(1,"(ES14.6,2I4)") cc(l,m),l,m + end do + end do + + end if + + close(1) + +! Area. + + if (ahf_ncall.eq.1) then + open(1,file=areaf,form='formatted', + . status='replace') + write(1,"(A14)") '" Horizon area' + else + open(1,file=areaf,form='formatted', + . status='old',position='append') + end if + if (status) write(1,"(2ES14.6)") cctk_time,intarea_h + close(1) + +! Mass. + + if (ahf_ncall.eq.1) then + open(1,file=massf,form='formatted', + . status='replace') + write(1,"(A14)") '" Horizon mass' + else + open(1,file=massf,form='formatted', + . status='old',position='append') + end if + if (status) write(1,"(2ES14.6)") cctk_time,0.141047396D0*dsqrt(intarea_h) + close(1) + +! Circumferences. + + if (ahf_ncall.eq.1) then + open(1,file=circeqf,form='formatted', + . status='replace') + write(1,"(A26)") '" Equatorial circumference' + else + open(1,file=circeqf,form='formatted', + . status='old',position='append') + end if + if (status) write(1,"(2ES14.6)") cctk_time,circ_eq + close(1) + + if (ahf_ncall.eq.1) then + open(1,file=merip1f,form='formatted', + . status='replace') + write(1,"(A28)") '" Length of meridian, phi=0' + else + open(1,file=merip1f,form='formatted', + . status='old',position='append') + end if + if (status) write(1,"(2ES14.6)") cctk_time,meri_p1 + close(1) + + if (ahf_ncall.eq.1) then + open(1,file=merip2f,form='formatted', + . status='replace') + write(1,"(A31)") '" Length of meridian, phi=pi/2' + else + open(1,file=merip2f,form='formatted', + . status='old',position='append') + end if + if (status) write(1,"(2ES14.6)") cctk_time,meri_p2 + close(1) +#ifdef MPI + end if +#endif + +! ******************************************************** +! *** PREPARING {ahfgrid3,ahf_exp3} FOR OUTPUT *** +! ******************************************************** + + if (find3) then + call AHFinder_find3(CCTK_FARGUMENTS,mtype) + end if + + +! ********************************* +! *** OUTPUT GRID FUNCTIONS *** +! ********************************* + +! if (contains("ahf_2Doutput","yes").ne.0) then + if (ahf_2Doutput.ne.0) then + if (find3) then +! if (mfind.eq.2) then +! call AHFinder_2dio(ahfgrid3,ahf_ncall) +! call AHFinder_2dio(ahf_exp3,ahf_ncall) +! +! call io_write2d(ahfgrid3) +! call io_write2d(ahf_exp3) +! +! call AHFinder_2dio(ahfgrid3,0) +! call AHFinder_2dio(ahf_exp3,0) +! end if +! else +! call AHFinder_2dio(ahfgrid,ahf_ncall) +! call AHFinder_2dio(ahf_exp,ahf_ncall) +! +! call io_write2d(ahfgrid) +! call io_write2d(ahf_exp) +! +! call AHFinder_2dio(ahfgrid,0) +! call AHFinder_2dio(ahf_exp,0) + end if + end if + +! if (contains("ahf_3Doutput","yes").ne.0) then + if (ahf_3Doutput.ne.0) then + +! This is commented out because I still need +! to write the c routine `AHFinder_3dio'. + +! call AHFinder_2dio(ahfgrid,ahf_ncall) +! call AHFinder_2dio(ahf_exp,ahf_ncall) + +! call io_write2d(ahfgrid) +! call io_write2d(ahf_exp) + +! call AHFinder_2dio(ahfgrid,0) +! call AHFinder_2dio(ahf_exp,0) + + end if + + +! **************** +! *** MASK *** +! **************** + +! if (horizon.and.CONTAINS("ahf_mask","yes")) then + if (horizon.and.(ahf_mask.ne.0)) then + call AHFinder_mask(CCTK_FARGUMENTS) +! else if ((mtype.eq.1).and.CONTAINS("ahf_weakmask","yes")) then + else if ((mtype.eq.1).and.(ahf_weakmask.ne.0)) then + call AHFinder_mask(CCTK_FARGUMENTS) + else + ahmask = one + end if + +#ifdef MPI + call synconefunc(ahmask) +#endif + +! ************************************************** +! *** TRANSFORM PHYSICAL TO CONFORMAL METRIC *** +! ************************************************** + +! call phys2conf(CONF2PHYSVARS_ARGS) + + +! ********************************************** +! *** UPDATING COUNTER AND LOOP IF FIND3 *** +! ********************************************** + + if (find3) then + if (mfind.eq.2) then + mfind = 0 + else + mfind = mfind + 1 + goto 1 + end if + end if + + +! *************** +! *** END *** +! *************** + + end subroutine AHFinder + + + + diff --git a/src/AHFinder_2dio.c b/src/AHFinder_2dio.c new file mode 100644 index 0000000..2fa2d3b --- /dev/null +++ b/src/AHFinder_2dio.c @@ -0,0 +1,14 @@ + +#include "cctk.h" + +void FMODIFIER FORTRAN_NAME(ahfinder_2dio_, + AHFINDER_2DIO, + ahfinder_2dio) + (Double *x,int *new_do_2dio) +{ + pGF *GF; + + GF = locateGFbyDataPtr(x); + GF->do_2dio = *new_do_2dio; +} + diff --git a/src/AHFinder_calcsigma.F b/src/AHFinder_calcsigma.F new file mode 100644 index 0000000..1cfaf9b --- /dev/null +++ b/src/AHFinder_calcsigma.F @@ -0,0 +1,311 @@ +c/*@@ +c @file AHFinder_calcsigma.F +c @date July 1999 +c @author Lars Nerger +c @desc +c Calculate weight sigma for Nakamura flow +c @enddesc +c@@*/ + +c Note that including cactus.h will also include AHFinder.h +!#include "cactus.h" +#include "cctk.h" +#include "cctk_parameters.h" +#include "cctk_arguments.h" + + + subroutine AHFinder_calcsigma(CCTK_FARGUMENTS,xp,yp,zp,gupij,sigma) + + use AHFinder_dat + + implicit none + + DECLARE_CCTK_FARGUMENTS + DECLARE_CCTK_PARAMETERS + + CCTK_INT nmax + INTEGER l,m,indx1,indx2,indx,idir,jdir + + CCTK_REAL zero, sigma + CCTK_REAL ni(3),si(3),supi(3),pij(3,3) + CCTK_REAL ylmi((lmax+1)**2,3) + CCTK_REAL alm((lmax+1)**2) + CCTK_REAL xp,yp,zp,er,sum + CCTK_REAL rho,costheta,sintheta,sinthetal(-1:lmax),prod + CCTK_REAL zlm0(0:lmax,0:lmax),zlm1(0:lmax,0:lmax) + CCTK_REAL cosphi,sinphi,cosmphi(0:lmax),sinmphi(0:lmax) + CCTK_REAL thetax,thetay,thetaz + CCTK_REAL phix,phiy + CCTK_REAL c1,c2,c3,metnorm + CCTK_REAL fi(3) + CCTK_REAL gupij(3,3) + +! Description of variables +! +! sigma Weight of Nakamura Flow + + +! ************************** +! *** DEFINE NUMBERS *** +! ************************** + + zero = 0.0D0 + nmax = (lmax+1)**2 + +! ****************************************************** +! *** Calculate gradient of Y_lm(x,y,z) *** +! *** Adopted from thorn_SpectralAHF by Gundlach *** +! ****************************************************** + + +! Initialize arrays + ylmi = zero + zlm0 = zero + zlm1 = zero + ni = zero + si = zero + supi = zero + pij = zero + alm = zero + sinthetal = zero + cosmphi = zero + sinmphi = zero + + +! Calculate quantities in spherical coordinates. + er = sqrt(xp**2 + yp**2 + zp**2) + if (er .eq. 0.d0) stop 'r=0 in AHFinder_calcsigma' + +C ********PATCH****************************************** +C Keep away from the z-axis, artificially. + if (xp**2 + yp**2 .lt. 1.d-16) then +c write(6,*) 'patch on axis in ylmder2' + xp = 0.d0 + yp = 1.d-8 + end if + +C Polar angle theta. + rho = sqrt(xp**2 + yp**2) + costheta = zp / er + sintheta = rho / er +C Array sinthetal(l) contains [sin(theta)]^l + sinthetal(-1) = 0.d0 + sinthetal(0) = 1.d0 + sinthetal(1) = sintheta + do l=2,lmax + sinthetal(l) = sinthetal(l-1) * sinthetal(1) + end do + +C zlm0 is Y_lm stripped of exp(i m phi) and multiplied by sqrt(4pi). +C zlm1 and zlm2 are the first and second derivative of zlm0 +C with respect to theta. l=0 and m=l are treated separately. +C We only need zlm0,1,2 for positive m. +C Treat l=0 separately: + zlm0(0,0) = 1.d0 + zlm1(0,0) = 0.d0 + +C Other l>0. + prod = 1.d0 + do l=1,lmax + +C Initialize m = l: +C c1 = sqrt(2l+1) * sqrt(2l!) / 2^l / l! +C Y_ll = c1 * sin(theta)^l +C Y_ll' = c1 * l * sin(theta)^(l-1) * costheta + m = l + prod = - prod * sqrt(dble(2*l-1) / dble(2*l)) + c1 = sqrt(dble(2*l+1)) * prod + zlm0(l,m) = c1 * sinthetal(l) + zlm1(l,m) = c1 * dble(l) * sinthetal(l-1) * costheta + +C Recursion relations for the other m: +C c2 = sqrt[(2l+1)(2l-1)/(l^2-m^2)] +C c3 = sqrt[((l-1)^2-m^2)(2l=1)/(l^2-m^2)/(2l-3)] +C Y_lm = c2 Y_l-1,m cos(theta) + c3 Y_l-2,m +C Y_lm' = c2 Y_l-1,m' cos(theta) + c2 Y_l-1,m + c3 Y_l-2,m' + m = l-1 + c2 = sqrt(dble((2*l+1)*(2*l-1)) + $ / dble(l**2-m**2)) + zlm0(l,m) = c2 * zlm0(l-1,m) * costheta + zlm1(l,m) = c2 * (zlm1(l-1,m) * costheta + $ - zlm0(l-1,m) * sintheta) + + do m=0,l-2 + c2 = sqrt(dble((2*l+1)*(2*l-1)) + $ / dble(l**2-m**2)) + c3 =sqrt(dble((2*l+1)*((l-1)**2-m**2)) + $ / dble((l**2 - m**2)*(2*l-3))) + zlm0(l,m) = c2 * zlm0(l-1,m) * costheta + $ - c3 * zlm0(l-2,m) + zlm1(l,m) = c2 * (zlm1(l-1,m) * costheta + $ - zlm0(l-1,m) * sintheta) + $ - c3 * zlm1(l-2,m) + end do + end do + + + if (rho .ne. 0.d0) then +C Azimuth angle phi is defined. + cosphi = xp / rho + sinphi = yp / rho +C First and second derivatives of theta. + thetax = costheta * cosphi / er + thetay = costheta * sinphi / er + thetaz = - sintheta / er +C First and second derivatives of phi, the angle in the xy-plane. + phix = - sinphi / rho + phiy = cosphi / rho +C Put together arrays of basis functions in the xy plane. +C The sqrt(2) factor is needed in going from the basis exp(+/-imphi) +C to the basis cos/sin(mphi). + cosmphi(0) = 1.d0 + sinmphi(0) = 1.d0 + cosmphi(1) = sqrt(2.d0) * cosphi + sinmphi(1) = sqrt(2.d0) * sinphi + do m=2,lmax + cosmphi(m) = cosmphi(m-1) * cosphi - sinmphi(m-1) * sinphi + sinmphi(m) = cosmphi(m-1) * sinphi + sinmphi(m-1) * cosphi + end do + else +C See above for a patch that avoids ever getting here. + stop 'ylmder2 cannot handle axis' +C For x = y = 0 the following are not defined. +! cosphi = 7.d77 +! sinphi = 7.d77 +! thetax = 7.d77 +! thetay = 7.d77 +! phix = 7.d77 +! phiy = 7.d77 +! cosmphi = 7.d77 +! sinmphi = 7.d77 +C Only these are defined. + thetaz = 0.d0 + cosmphi(0) = 1.d0 + sinmphi(0) = 1.d0 + end if + +C Assemble the Y_lm from zlm0 and the basis of functions of phi. +C Assemble the Y_lm,i and Y_lm,ij in the corresponding manner. + do l=0,lmax + do m=0,l +C These are the values of indx corresponding to (l,m) and (l,-m). + indx1 = l**2 + 1 + l + m + indx2 = l**2 + 1 + l - m +C Ylm,x + ylmi(indx1,1) = thetax * zlm1(l,m) * cosmphi(m) + $ - dble(m) * phix * zlm0(l,m) * sinmphi(m) + ylmi(indx2,1) = thetax * zlm1(l,m) * sinmphi(m) + $ + dble(m) * phix * zlm0(l,m) * cosmphi(m) +C Ylm,y + ylmi(indx1,2) = thetay * zlm1(l,m) * cosmphi(m) + $ - dble(m) * phiy * zlm0(l,m) * sinmphi(m) + ylmi(indx2,2) = thetay * zlm1(l,m) * sinmphi(m) + $ + dble(m) * phiy * zlm0(l,m) * cosmphi(m) +C Ylm,z + ylmi(indx1,3) = thetaz * zlm1(l,m) * cosmphi(m) + ylmi(indx2,3) = thetaz * zlm1(l,m) * sinmphi(m) + end do + end do + + +! ******************************** +! *** First derivatives of r *** +! ******************************** + + ni(1) = xp / er + ni(2) = yp / er + ni(3) = zp / er + + +! ************************************************************* +! *** compress c0, cc and cs into alm with single index *** +! ************************************************************* + + do l=0, lmax + alm(l**2+l+1) = c0(l) + do m=1,lmax + alm(l**2+l+1+m) = cc(l,m) + alm(l**2+l+1-m) = cs(l,m) + end do + end do + + +! ******************************************************************* +! *** Assemble derivatives of f. Recall f = r - h(theta,phi): *** +! ******************************************************************* + +C Derivatives of r... +C r,i = n_i and r,ij = (delta_ij - n_i n_j) / r = n_ij + do idir=1,3 + fi(idir) = ni(idir) + end do +C ...and derivatives of -h. + do indx=1,nmax + do idir=1,3 + fi(idir) = fi(idir) - alm(indx) * ylmi(indx,idir) + end do + end do + + +! ************************************************************* +! *** Assemble normal vector s on surfaces of constant f *** +! *** with indices down. That is, normal with respect *** +! *** to the physical metric g_ij. *** +! ************************************************************* + + metnorm = 0.d0 + do idir=1,3 + do jdir=1,3 + metnorm = metnorm + fi(idir) * fi(jdir) + $ * gupij(idir,jdir) + end do + end do + metnorm = sqrt(metnorm) + do idir=1,3 + si(idir) = fi(idir) / metnorm + end do + +C Now with indices up, raised by g^ij. + do idir=1,3 + supi(idir) = 0.d0 + do jdir=1,3 + supi(idir) = supi(idir) + $ + gupij(idir,jdir) * si(jdir) + end do + end do + +! ******************************************************* +! *** Assemble projection operator with indices up *** +! *** p^ij = g^ij - s^i s^j *** +! ******************************************************* + + do idir=1,3 + do jdir=1,3 + pij(idir,jdir) = gupij(idir,jdir) + $ - supi(idir) * supi(jdir) + end do + end do + + +! ************************************************* +! *** Calculate weight sigma *** +! *** 2 * r^2 / (p^ij (delta_ij - n_i n_j)) *** +! ************************************************* + + sum = 0.d0 + do idir=1,3 + sum = sum + pij(idir,idir) + do jdir=1,3 + sum = sum + $ - pij(idir,jdir) * ni(idir) * ni(jdir) + end do + end do + sigma = 2.d0 * er**2 / sum + + +! *************** +! *** END *** +! *************** + + end subroutine AHFinder_calcsigma + diff --git a/src/AHFinder_dat.F b/src/AHFinder_dat.F new file mode 100644 index 0000000..efcc6ae --- /dev/null +++ b/src/AHFinder_dat.F @@ -0,0 +1,142 @@ +c/*@@ +c @file AHFinder_dat.F +c @date April 1998 +c @author Miguel Alcubierre +c @desc +c F90 module containing data to be shared by +c the different subroutines. +c +c For the meaning of the spectral components of +c the different flows see: +c +c C. Gundlach, Phys. Rev. D 57, 863 (1998). +c +c @enddesc +c@@*/ +#include "cctk.h" + module AHFinder_dat + + implicit none + + logical firstfun,firstleg,firstint + logical nonaxi,refx,refy,refz,cartoon + logical offset,wander + logical minarea,flow,find_trapped_surface + logical logfile,verbose,veryver,guessverbose + logical interror,guessold + logical find3,sloppy,inner,guess_absmin,manual_guess + + CCTK_INT lmax,stepx,stepy,stepz + CCTK_INT ntheta,nphi + CCTK_INT nfile + CCTK_INT interror1,interror2 + CCTK_INT mfind,ahf_ncall + CCTK_INT nx,ny,nz + integer nprocs,myproc + + + CCTK_REAL dx,dy,dz + CCTK_REAL xc,yc,zc + CCTK_REAL xc_old,yc_old,zc_old + CCTK_REAL intexp,intexp2,intexpdel2,intarea + + CCTK_REAL hw,cw,nw + CCTK_REAL inside_min_count,inside_min_neg_count + CCTK_REAL dx,dy,dz + + CCTK_REAL, allocatable, dimension (:) :: c0,c0_old + CCTK_REAL, allocatable, dimension (:,:) :: cc,cs,cc_old,cs_old + + CCTK_REAL, allocatable, dimension (:) :: hflow0,cflow0,nflow0 + CCTK_REAL, allocatable, dimension (:,:) :: hflowc,cflowc,nflowc + CCTK_REAL, allocatable, dimension (:,:) :: hflows,cflows,nflows + + + character*200 filestr + + data mfind / 0 / + data ahf_ncall / 0 / + + save mfind,ahf_ncall + save c0,cc,cs + +! Description of variables: +! +! firstfun First call to function FUNC? +! firstleg First call to function LEGEN? +! firstint First call to subroutine AHFinder_int? +! +! nonaxi Is the data non-axisymmetric? +! +! refx Reflection symmetry x->-x? +! +! refy Reflection symmetry y->-y? +! +! refz Reflection symmetry z->-z? +! +! cartoon Are we using cartoon_2d? +! +! inner Look for inner horizon? +! +! sloppy Sloppy initial guess? +! +! offset Is the center offset from the origin? +! wander Do we allow the center to wander? +! +! minarea Do we want to minimize area? +! +! flow Use fast flow algorithm instead of minimization? +! +! logfile Write log file? +! +! verbose Write messages to screen? +! +! veryver Write lots of messages to screen? +! +! guessverbose Write info on initial guess points? +! +! lmax Maximum number of terms in theta expansion. +! +! stepx 1 if (refx=1), 0 otherwise. +! stepy 1 if (refy=1), 0 otherwise. +! stepz 1 if (refz=1), 0 otherwise. +! +! ntheta Number of subdivisions in theta. +! nphi Number of subdivisions in phi. +! +! xc Coordinate x of centre of surface. +! yc Coordinate y of centre of surface. +! zc Coordinate z of centre of surface. +! +! intexp Integral of expansion over the surface. +! intexp2 Integral of expansion squared over the surface. +! intarea Area of the surface. +! +! c0(l) Coefficients of LEGEN(l,0,cos(theta)) +! cc(l,m) Coefficients of LEGEN(l,m,cos(theta))*cos(m*phi) +! cs(l,m) Coefficients of LEGEN(l,m,cos(theta))*sin(m*phi) +! +! hflow* Spectral components for `H' flow. +! cflow* Spectral components for `C' flow. +! nflow* Spectral components for `N' flow. +! +! interror1 Different from zero if radius is negative. +! interror2 Different from zero if outside computational domain. +! interror True if either interror1 or interro2 are non-zero. +! +! inside_min_count Number of elements in integral (should always +! be ntheta*nphi) +! inside_min_neg_count Number of elements in integral that have +! negative expansion. +! +! find3 Look for 3 horizons? +! mfind + + +! *************** +! *** END *** +! *************** + + end module AHFinder_dat + + diff --git a/src/AHFinder_exp.F b/src/AHFinder_exp.F new file mode 100644 index 0000000..9ccea39 --- /dev/null +++ b/src/AHFinder_exp.F @@ -0,0 +1,429 @@ +c/*@@ +c @file AHFinder_exp.F +c @date April 1998 +c @author Miguel Alcubierre +c @desc +c Here I calculate the (normalized) expansion of the +c horizon function. The expansion is defined by the +c following expression (a=lapse, b=shift): +c +c __2 2 / ab __a__b \ +c exp = \/ f / u + d f d f / u | K - \/ \/ f / u | - trK +c a b \ / +c +c where: +c +c / mn \ 1/2 +c u = | g d f d f | +c \ m n / +c +c @enddesc +c@@*/ + +c Note that including cactus.h will also include AHFinder.h +!#include "cactus.h" +#include "cctk.h" +#include "cctk_parameters.h" +#include "cctk_arguments.h" + + subroutine AHFinder_exp(CCTK_FARGUMENTS) + + use AHFinder_dat + + implicit none + + DECLARE_CCTK_FARGUMENTS + DECLARE_CCTK_PARAMETERS + + integer i,j,k + INTEGER CCTK_Equals + + CCTK_REAL det,idet + CCTK_REAL gdxx,gdyy,gdzz,gdxy,gdxz,gdyz + CCTK_REAL guxx,guyy,guzz,guxy,guxz,guyz + CCTK_REAL kdxx,kdyy,kdzz,kdxy,kdxz,kdyz + CCTK_REAL ddxxx,ddxyy,ddxzz,ddxxy,ddxxz,ddxyz + CCTK_REAL ddyxx,ddyyy,ddyzz,ddyxy,ddyxz,ddyyz + CCTK_REAL ddzxx,ddzyy,ddzzz,ddzxy,ddzxz,ddzyz + CCTK_REAL d2fxx,d2fyy,d2fzz,d2fxy,d2fxz,d2fyz + CCTK_REAL c2fxx,c2fyy,c2fzz,c2fxy,c2fxz,c2fyz + CCTK_REAL dfx,dfy,dfz,dfux,dfuy,dfuz + CCTK_REAL aux,T0,T1,T2,T3,T4 + CCTK_REAL idx,idy,idz,idx2,idy2,idz2,idxy,idxz,idyz + CCTK_REAL zero,half,one,two + +! Description of variables: +! +! i,j,k, counters +! +! gdij g +! ij +! ij +! guij g +! +! +! det det(g) +! +! idet 1/det +! +! +! kdij K +! ij +! +! ddijk 1/2 d g +! i jk +! +! dfi d f +! i +! ij +! dfui g d f +! j +! +! d2fij d d f +! i j +! +! c2fij f +! ;ij +! +! aux,T* auxilliary variables + + +! *************************************** +! *** DEFINE {zero,half,one,two} *** +! *************************************** + + zero = 0.0d0 + half = 0.5d0 + one = 1.0d0 + two = 2.0d0 + +! ************************** +! *** FIND EXPANSION *** +! ************************** + + idx = half/dx + idy = half/dy + idz = half/dz + + idxy = idx*idy + idxz = idx*idz + idyz = idy*idz + + idx2 = one/dx**2 + idy2 = one/dy**2 + idz2 = one/dz**2 + + do k=2,nz-1 + do j=2,ny-1 + do i=2,nx-1 + +! Find spatial metric. + + gdxx = gxx(i,j,k) + gdyy = gyy(i,j,k) + gdzz = gzz(i,j,k) + gdxy = gxy(i,j,k) + gdxz = gxz(i,j,k) + gdyz = gyz(i,j,k) + +! Find extrinsic curvature. + + kdxx = kxx(i,j,k) + kdyy = kyy(i,j,k) + kdzz = kzz(i,j,k) + kdxy = kxy(i,j,k) + kdxz = kxz(i,j,k) + kdyz = kyz(i,j,k) + +! Find D's using finite differences. + + ddxxx = idx*(gxx(i+1,j,k) - gxx(i-1,j,k)) + ddxyy = idx*(gyy(i+1,j,k) - gyy(i-1,j,k)) + ddxzz = idx*(gzz(i+1,j,k) - gzz(i-1,j,k)) + ddxxy = idx*(gxy(i+1,j,k) - gxy(i-1,j,k)) + ddxxz = idx*(gxz(i+1,j,k) - gxz(i-1,j,k)) + ddxyz = idx*(gyz(i+1,j,k) - gyz(i-1,j,k)) + + ddyxx = idy*(gxx(i,j+1,k) - gxx(i,j-1,k)) + ddyyy = idy*(gyy(i,j+1,k) - gyy(i,j-1,k)) + ddyzz = idy*(gzz(i,j+1,k) - gzz(i,j-1,k)) + ddyxy = idy*(gxy(i,j+1,k) - gxy(i,j-1,k)) + ddyxz = idy*(gxz(i,j+1,k) - gxz(i,j-1,k)) + ddyyz = idy*(gyz(i,j+1,k) - gyz(i,j-1,k)) + + ddzxx = idz*(gxx(i,j,k+1) - gxx(i,j,k-1)) + ddzyy = idz*(gyy(i,j,k+1) - gyy(i,j,k-1)) + ddzzz = idz*(gzz(i,j,k+1) - gzz(i,j,k-1)) + ddzxy = idz*(gxy(i,j,k+1) - gxy(i,j,k-1)) + ddzxz = idz*(gxz(i,j,k+1) - gxz(i,j,k-1)) + ddzyz = idz*(gyz(i,j,k+1) - gyz(i,j,k-1)) + +! Find determinant of spatial metric. + + det = gdxx*gdyy*gdzz + two*gdxy*gdxz*gdyz + . - gdxx*gdyz**2 - gdyy*gdxz**2 - gdzz*gdxy**2 + + if (det.le.zero) then + ahf_exp(i,j,k) = zero + goto 10 + end if + + idet = one/det + +! Find inverse spatial metric. + + guxx = idet*(gdyy*gdzz - gdyz**2) + guyy = idet*(gdxx*gdzz - gdxz**2) + guzz = idet*(gdxx*gdyy - gdxy**2) + + guxy = idet*(gdxz*gdyz - gdzz*gdxy) + guxz = idet*(gdxy*gdyz - gdyy*gdxz) + guyz = idet*(gdxy*gdxz - gdxx*gdyz) + +! Find spatial derivatives of f. + + T0 = two*ahfgrid(i,j,k) + + T1 = ahfgrid(i+1,j,k) + T2 = ahfgrid(i-1,j,k) + + dfx = (T1 - T2)*idx + d2fxx = (T1 - T0 + T2)*idx2 + + T1 = ahfgrid(i,j+1,k) + T2 = ahfgrid(i,j-1,k) + + dfy = (T1 - T2)*idy + d2fyy = (T1 - T0 + T2)*idy2 + + T1 = ahfgrid(i,j,k+1) + T2 = ahfgrid(i,j,k-1) + + dfz = (T1 - T2)*idz + d2fzz = (T1 - T0 + T2)*idz2 + +! Save gradient of horizon function and its norm +! (they will be needed later in the flow algorithm +! and in the gaussian curvature). Here I also try +! to avoid possible division by zero later by not +! allowing ahfgradn(i,j,k) = 0. This should only +! ever happen far from the horizon, so resetting this +! to 1 should have no important effects. + + ahfgradx(i,j,k) = dfx + ahfgrady(i,j,k) = dfy + ahfgradz(i,j,k) = dfz + + aux = guxx*dfx**2 + guyy*dfy**2 + guzz*dfz**2 + . + two*(guxy*dfx*dfy + guxz*dfx*dfz + guyz*dfy*dfz) + + ahfgradn(i,j,k) = dsqrt(aux) + if (ahfgradn(i,j,k).eq.zero) ahfgradn(i,j,k) = one + +! Find crossed derivatives. + + d2fxy = (ahfgrid(i+1,j+1,k) + ahfgrid(i-1,j-1,k) + . - ahfgrid(i+1,j-1,k) - ahfgrid(i-1,j+1,k))*idxy + d2fxz = (ahfgrid(i+1,j,k+1) + ahfgrid(i-1,j,k-1) + . - ahfgrid(i+1,j,k-1) - ahfgrid(i-1,j,k+1))*idxz + d2fyz = (ahfgrid(i,j+1,k+1) + ahfgrid(i,j-1,k-1) + . - ahfgrid(i,j+1,k-1) - ahfgrid(i,j-1,k+1))*idyz + +! Raise indices in first derivatives. + + dfux = guxx*dfx + guxy*dfy + guxz*dfz + dfuy = guxy*dfx + guyy*dfy + guyz*dfz + dfuz = guxz*dfx + guyz*dfy + guzz*dfz + +! Find second covariant derivatives of f. + + c2fxx = d2fxx - half*(dfux*ddxxx + . + dfuy*(two*ddxxy - ddyxx) + . + dfuz*(two*ddxxz - ddzxx)) + + c2fyy = d2fyy - half*(dfuy*ddyyy + . + dfux*(two*ddyxy - ddxyy) + . + dfuz*(two*ddyyz - ddzyy)) + + c2fzz = d2fzz - half*(dfuz*ddzzz + . + dfux*(two*ddzxz - ddxzz) + . + dfuy*(two*ddzyz - ddyzz)) + + c2fxy = d2fxy - half*(dfux*ddyxx + dfuy*ddxyy + . + dfuz*(ddxyz + ddyxz - ddzxy)) + + c2fxz = d2fxz - half*(dfux*ddzxx + dfuz*ddxzz + . + dfuy*(ddxyz + ddzxy - ddyxz)) + + c2fyz = d2fyz - half*(dfuy*ddzyy + dfuz*ddyzz + . + dfux*(ddyxz + ddzxy - ddxyz)) + +! / m \ 1/2 +! Find: u = | d f d f | +! \ m / + + T0 = dfx*dfux + dfy*dfuy + dfz*dfuz + + if (T0.gt.zero) then + T0 = one/dsqrt(T0) + else + ahf_exp(i,j,k) = zero + goto 10 + end if + +! __2 +! Find: \/ f / u + + T1 = guxx*c2fxx + guyy*c2fyy + guzz*c2fzz + . + two*(guxy*c2fxy + guxz*c2fxz + guyz*c2fyz) + + T1 = T1*T0 + +! a b 2 +! Find: K d f d f / u +! ab + + T2 = kdxx*dfux**2 + kdyy*dfuy**2 + kdzz*dfuz**2 + . + two*(dfux*(kdxy*dfuy + kdxz*dfuz) + kdyz*dfuy*dfuz) + + T2 = T2*T0**2 + +! __a__b 3 +! Find: \/ \/ f d f d f / u +! a b + + T3 = c2fxx*dfux**2 + c2fyy*dfuy**2 + c2fzz*dfuz**2 + . + two*(dfux*(c2fxy*dfuy + c2fxz*dfuz) + . + c2fyz*dfuy*dfuz) + + T3 = T3*T0**3 + +! Find: trK + + T4 = guxx*kdxx + guyy*kdyy + guzz*kdzz + . + two*(guxy*kdxy + guxz*kdxz + guyz*kdyz) + +! Find the expansion. + + ahf_exp(i,j,k) = T1 + T2 - T3 - T4 + +! Check for NaN's. + + if (ahf_exp(i,j,k).ge.zero) then + else if (ahf_exp(i,j,k).lt.zero) then + else +#ifdef MPI + if (myproc.eq.0) then +#endif + write(*,*) 'NaN in ahf_exp at point',i,j,k +#ifdef MPI + end if +#endif + STOP + end if + + 10 continue + + end do + end do + end do + +! Boundaries. + + do k=1,nz + do j=1,ny + + ahf_exp(1 ,j,k) = ahf_exp(2 ,j,k) + ahf_exp(nx,j,k) = ahf_exp(nx-1,j,k) + + ahfgradx(nx,j,k) = ahfgradx(nx-1,j,k) + ahfgrady(nx,j,k) = ahfgrady(nx-1,j,k) + ahfgradz(nx,j,k) = ahfgradz(nx-1,j,k) + ahfgradn(nx,j,k) = ahfgradn(nx-1,j,k) + + ahfgrady(1 ,j,k) = ahfgrady(2,j,k) + ahfgradz(1 ,j,k) = ahfgradz(2,j,k) + ahfgradn(1 ,j,k) = ahfgradn(2,j,k) + +! if ((contains("grid","octant").ne.0).or. +! . (contains("grid","quadrant").ne.0)) then + if ((CCTK_Equals(symmetry,"octant").eq.1).or. + . (CCTK_Equals(symmetry,"quadrant").eq.1)) then + + ahfgradx(1,j,k) = -ahfgradx(2,j,k) + else + ahfgradx(1,j,k) = +ahfgradx(2,j,k) + end if + + end do + end do + + do k=1,nz + do i=1,nx + + ahf_exp(i,1 ,k) = ahf_exp(i,2 ,k) + ahf_exp(i,ny,k) = ahf_exp(i,ny-1,k) + + ahfgradx(i,ny,k) = ahfgradx(i,ny-1,k) + ahfgrady(i,ny,k) = ahfgrady(i,ny-1,k) + ahfgradz(i,ny,k) = ahfgradz(i,ny-1,k) + ahfgradn(i,ny,k) = ahfgradn(i,ny-1,k) + + ahfgradx(i,1 ,k) = ahfgradx(i,2,k) + ahfgradz(i,1 ,k) = ahfgradz(i,2,k) + ahfgradn(i,1 ,k) = ahfgradn(i,2,k) + +! if ((contains("grid","octant").ne.0).or. +! . (contains("grid","quadrant").ne.0)) then + if ((CCTK_Equals(symmetry,"octant").eq.1).or. + . (CCTK_Equals(symmetry,"quadrant").eq.1)) then + ahfgrady(i,1,k) = -ahfgrady(i,2,k) + else + ahfgrady(i,1,k) = +ahfgrady(i,2,k) + end if + + end do + end do + + do j=1,ny + do i=1,nx + + ahf_exp(i,j,1 ) = ahf_exp(i,j,2 ) + ahf_exp(i,j,nz) = ahf_exp(i,j,nz-1) + + ahfgradx(i,j,nz) = ahfgradx(i,j,nz-1) + ahfgrady(i,j,nz) = ahfgrady(i,j,nz-1) + ahfgradz(i,j,nz) = ahfgradz(i,j,nz-1) + ahfgradn(i,j,nz) = ahfgradn(i,j,nz-1) + + ahfgradx(i,j,1 ) = ahfgradx(i,j,2) + ahfgrady(i,j,1 ) = ahfgrady(i,j,2) + ahfgradn(i,j,1 ) = ahfgradn(i,j,2) + +! if (contains("grid","octant").ne.0) then + if (CCTK_Equals(symmetry,"octant").eq.1) then + ahfgradz(i,j,1) = -ahfgradz(i,j,2) + else + ahfgradz(i,j,1) = +ahfgradz(i,j,2) + end if + + end do + end do + +! Synchronize. + +#ifdef MPI + call synconefunc(ahf_exp) + call synconefunc(ahfgradx) + call synconefunc(ahfgrady) + call synconefunc(ahfgradz) + call synconefunc(ahfgradn) +#endif + +! *************** +! *** END *** +! *************** + + return + end + diff --git a/src/AHFinder_find3.F b/src/AHFinder_find3.F new file mode 100644 index 0000000..951bb4b --- /dev/null +++ b/src/AHFinder_find3.F @@ -0,0 +1,68 @@ +c/*@@ +c @file AHFinder_find3.F +c @date March 1999 +c @author Lars Nerger +c @desc +c Muliplicate gridfunctions for final grid function +c @enddesc +c@@*/ + +c Note that including cactus.h will also include AHFinder.h +!#include "cactus.h" +#include "cctk.h" +#include "cctk_parameters.h" +#include "cctk_arguments.h" + +! To output single gridfunctions when searching for 3 horizons +! the gridfunctions are here multiplicated. + + subroutine AHFinder_find3(CCTK_FARGUMENTS,mtype) + + use AHFinder_dat + + implicit none + + DECLARE_CCTK_FARGUMENTS + DECLARE_CCTK_PARAMETERS + + integer mtype + + REAL one + +! Description of variables +! +! mtype Type of surface: see AHFinder.F + + +! ********************** +! *** DEFINE one *** +! ********************** + + one = 1.0D0 + + +! ************************* +! *** START ROUTINE *** +! ************************* + + if ((mfind.eq.0).and.(mtype.eq.1)) then + ahfgrid3 = ahfgrid + ahf_exp3 = ahf_exp + else if ((mfind.eq.0).and.(mtype.ne.1)) then + ahfgrid3 = one + ahf_exp3 = one + else if ((mfind.eq.1).and.(mtype.eq.1)) then + ahfgrid3 = ahfgrid*ahfgrid3 + ahf_exp3 = ahf_exp*ahf_exp3 + else if ((mfind.eq.2).and.(mtype.eq.1)) then + ahfgrid3 = ahfgrid*ahfgrid3 + ahf_exp3 = ahf_exp*ahf_exp3 + end if + + +! *************** +! *** END *** +! *************** + + end subroutine AHFinder_find3 + diff --git a/src/AHFinder_flow.F b/src/AHFinder_flow.F new file mode 100644 index 0000000..a44600d --- /dev/null +++ b/src/AHFinder_flow.F @@ -0,0 +1,660 @@ +c/*@@ +c @file AHFinder_flow.F +c @date November 1998 +c @author Miguel Alcubierre +c @desc +c Master routine to control the iterations for +c the flow algorithm. +c +c The flow algorithm used here is that of Carsten Gunlach: +c Phys. Rev. D 57, 863 (1998). +c +c @enddesc +c@@*/ + +c Note that including cactus.h will also include AHFinder.h +!#include "cactus.h" +#include "cctk.h" +#include "cctk_parameters.h" +#include "cctk_arguments.h" + + subroutine AHFinder_flow(CCTK_FARGUMENTS,rmx,status,logf) + + use AHFinder_dat + + implicit none + + DECLARE_CCTK_FARGUMENTS + DECLARE_CCTK_PARAMETERS + + logical status,new + + CCTK_INT ITMAX + integer i,j,l,m + + CCTK_REAL rmx,ahftol + CCTK_REAL A,B + CCTK_REAL ahfsum,ahfdiff,ahfdiff_old + CCTK_REAL maxchange,minchange + CCTK_REAL reldiff + CCTK_REAL spher + CCTK_REAL zero,one + CCTK_REAL dd,aux + + character*200 logf + +! Description of variables: +! +! i,l,m Counters. +! +! ITMAX Maximum number of iterations. +! +! A,B Flow parameters. +! +! hw Weigth of `H' flow (see Carsten's paper). +! cw Weigth of `C' flow (see Carsten's paper). +! nw Weigth of `N' flow (see Carsten's paper). +! +! ahftol Tolerance used to decide if we have converged. + + +! *************************** +! *** DEFINE NUMBERS *** +! *************************** + + zero = 0.0D0 + one = 1.0D0 + + dd = min(dx,dy,dz) + + +! *************************** +! *** FIND PARAMETERS *** +! *************************** + +! A = getr8("ahf_flowa") +! B = getr8("ahf_flowb") + A = ahf_flowa + B = ahf_flowb + +! hw = getr8("ahf_flowh") +! cw = getr8("ahf_flowc") +! nw = getr8("ahf_flown") + hw = ahf_flowh + cw = ahf_flowc + nw = ahf_flown + +! maxchange = getr8("ahf_maxchange") +! minchange = getr8("ahf_minchange") + + maxchange = ahf_maxchange + minchange = ahf_minchange + +! ************************* +! *** INITIAL GUESS *** +! ************************* + + 5 continue + +! Initialize arrays. + + c0 = zero + cc = zero + cs = zero + + hflow0 = zero + hflowc = zero + hflows = zero + + cflow0 = zero + cflowc = zero + cflows = zero + + nflow0 = zero + nflowc = zero + nflows = zero + +! As initial guess I take the largest sphere that +! can fit in the grid. + + if ((myproc.eq.0).and.verbose) then + write(*,*) + write(*,"(A35,ES11.3)") ' Initial guess sphere with radius: ',rmx + write(*,"(A13,3ES11.3)") ' centered on:',real(xc),real(yc),real(zc) + end if + + c0(0) = rmx + + +! ******************************* +! *** MAIN ITERATION LOOP *** +! ******************************* + +! Find ITMAX and ahftol. + +! ITMAX = geti("ahf_flowiter") +! ahftol = getr8("ahf_flowtol") + ITMAX = ahf_flowiter + ahftol = ahf_flowtol + + +! Start the iterations. + + do i=1,ITMAX + +! Save the old values of the coefficients. + + c0_old = c0 + cc_old = cc + cs_old = cs + +! In order to adapt the stepsize, I do two steps +! with the current stepsize, and compare them with +! one step with double the stepsize. + + do j=1,2 + + ahfsum = zero + ahfdiff = zero + + +! ************************************ +! *** FIND SPECTRAL COMPONENTS *** +! ************************************ + + call AHFinder_fun(CCTK_FARGUMENTS) + call AHFinder_exp(CCTK_FARGUMENTS) + call AHFinder_int(CCTK_FARGUMENTS) + + +! ***************************** +! *** CHECK ERROR FLAGS *** +! ***************************** + +! If interror is true there was an error in the integral. +! There are two possibilities here: + +! 1) This is the first of the two small steps (j=1). +! We then went out of bounds in the previous iteration, +! so there is nothing we can do to fix it. + + if (interror.and.(j.eq.1)) then + +! We can be out of bounds for two reasons: + +! a) We are out of the computational domain, pitty. + + if (interror2.ne.0) then + + write(*,*) + write(*,*) 'AHFinder: Out of bounds, giving up.' + return + +! b) The radius is negative somewhere. This is a really +! interesting case, since it probably means that the +! origin of our expansion is outside the horizon. +! This can happen is the black-hole is not centered, +! or if we have more than one horizon. + + else if (interror1.ne.0) then + + write(*,*) + write(*,*) 'AHFinder: Negative radius.' + return + +! We should never get here! + + else + + write(*,*) + write(*,*) 'Something is very wrong!' + return + + end if + +! 2) This is the second of the two small steps (j=2). +! We can still try to reduce the stepsize and start +! the iteration again. We then jump to the place +! where the stepsize is adapted, and force it to +! be reduced. + + else if (interror) then + + interror = .false. + reldiff = one + ahfsum = one + ahfdiff = one + intexp2 = one + goto 100 + + end if + + +! **************************************** +! *** UPDATE SPECTRAL COEFFICIENTS *** +! **************************************** + +! Update c0. + + do l=0,lmax,1+stepz + aux = c0(l) + c0(l) = aux - A/(one + B*dble(l)*(dble(l) + one)) + . *(hw*hflow0(l) + cw*cflow0(l) + nw*nflow0(l)) + ahfsum = ahfsum + c0_old(l)**2 + ahfdiff = ahfdiff + (c0(l) - c0_old(l))**2 + end do + +! Update {cc,cs}. + + if (nonaxi) then + +! Update cc(l,m). + + do l=1,lmax + do m=1+stepx,l,1+stepx + if (stepz*mod(l-m,2).eq.0) then + aux = cc(l,m) + cc(l,m) = aux - A/(one + . + B*dble(l)*(dble(l) + one)) + . *(hw*hflowc(l,m) + cw*cflowc(l,m) + . + nw*nflowc(l,m)) + ahfsum = ahfsum + cc_old(l,m)**2 + ahfdiff = ahfdiff + (cc(l,m)-cc_old(l,m))**2 + end if + end do + end do + +! Update cs(l,m). + + if (.not.refy) then + do l=1,lmax + do m=1,l,1+stepx + if (stepz*mod(l-m,2).eq.0) then + cs(l,m) = cs(l,m) - A/(one + . + B*dble(l)*(dble(l) + one)) + . *(hw*hflows(l,m) + cw*cflows(l,m) + . + nw*nflows(l,m)) + ahfsum = ahfsum + cs_old(l,m)**2 + ahfdiff = ahfdiff + (cs(l,m)-cs_old(l,m))**2 + end if + end do + end do + end if + + end if + + ahfsum = dsqrt(ahfsum) + ahfdiff = dsqrt(ahfdiff) + +! If this is the first of the two small steps, save +! twice the difference. Why twice? Because this would +! have been the resulting difference if we had done only +! one step that was twice as big. + + if (j.eq.1) then + ahfdiff_old = 2.0D0*ahfdiff + end if + + end do + + +! *************************** +! *** ADJUST STEPSIZE *** +! *************************** + +! Here we find the relative difference between having done two +! small steps and one large step. + + reldiff = dabs(ahfdiff - ahfdiff_old) + . / dabs(ahfdiff + ahfdiff_old) + + + 100 new = .false. + +! If the relative difference is too large, we reduce the +! stepsize and ignore this iteration. + + if ((reldiff.gt.maxchange).or. + . (ahfdiff.gt.maxchange*ahfsum)) then + + new = .true. + A = 0.5D0*A + + c0 = c0_old + cc = cc_old + cs = cs_old + +! If the diference is too small, we can safely increase +! the stepsize for the next iteration. + + else if (reldiff.lt.minchange) then + + new = .true. + A = 2.0D0*A + + end if + + +! ****************************************** +! *** LOGFILE AND MESSAGES TO SCREEN *** +! ****************************************** + +! Open logfile. + + logf = filestr(1:nfile)//"ahf_logfile" + + if (logfile.and.(myproc.eq.0)) then + open(1,file=logf,form='formatted',status='old', + . position='append') + end if + +! Write messages. + + if (myproc.eq.0) then + + if (veryver) then + + write(*,*) + write(*,*) + write(*,"(A28,I3)") ' AHFinder: FLOW ITERATION ',i + + if (intarea.ne.zero) then + aux = one/intarea + else + aux = one + end if + + write(*,*) + write(*,"(A21,ES14.6)") ' Surface area =',intarea + write(*,"(A21,ES14.6)") ' Mean value of H =',aux*intexp + write(*,"(A21,ES14.6)") ' Mean value of H^2 =',aux*intexp2 + write(*,"(A21,ES14.6)") ' ahfdiff/ahfsum =',ahfdiff/ahfsum + + if (new) then + write(*,*) + write(*,"(A15,ES14.6)") ' New stepsize =',A + end if + + if (offset) then + write(*,*) + write(*,"(A6,ES14.6)") ' xc =',xc + write(*,"(A6,ES14.6)") ' yc =',yc + write(*,"(A6,ES14.6)") ' zc =',zc + end if + + write(*,*) + write(*,"(A20)") ' Shape coefficients:' + write(*,*) + + write(*,"(A4,I2,A3,ES14.6)") ' c0(',0,') =',c0(0) + + do l=1+stepz,lmax,1+stepz + write(*,"(A4,I2,A3,ES14.6)") ' c0(',l,') =',c0(l) + end do + + end if + + if (logfile) then + + write(1,*) + write(1,"(A28,I3)") ' AHFinder: FLOW ITERATION ',i + + if (intarea.ne.zero) then + aux = one/intarea + else + aux = one + end if + + write(1,*) + write(1,"(A21,ES14.6)") ' Surface area =',intarea + write(1,"(A21,ES14.6)") ' Mean value of H =',aux*intexp + write(1,"(A21,ES14.6)") ' Mean value of H^2 =',aux*intexp2 + write(1,"(A21,ES14.6)") ' ahfdiff/ahfsum =',ahfdiff/ahfsum + + if (new) then + write(1,*) + write(1,"(A15,ES14.6)") ' New stepsize =',A + end if + + if (offset) then + write(1,*) + write(1,"(A6,ES14.6)") ' xc =',xc + write(1,"(A6,ES14.6)") ' yc =',yc + write(1,"(A6,ES14.6)") ' zc =',zc + end if + + write(1,*) + write(1,"(A20)") ' Shape coefficients:' + write(1,*) + + write(1,"(A4,I2,A3,ES14.6)") ' c0(',0,') =',c0(0) + + do l=1+stepz,lmax,1+stepz + write(1,"(A4,I2,A3,ES14.6)") ' c0(',l,') =',c0(l) + end do + + end if + + if (nonaxi) then + + if (.not.refy) then + + if (veryver) then + + write(*,*) + + do l=1,lmax + do m=1,l + if (stepz*mod(l-m,2).eq.0) then + write(*,"(A4,I2,A1,I2,A4,ES14.6)") ' cc(', + . l,',',m,') =',cc(l,m) + end if + end do + end do + write(*,*) + do l=1,lmax + do m=1,l + if (stepz*mod(l-m,2).eq.0) then + write(*,"(A4,I2,A1,I2,A4,ES14.6)") ' cs(', + . l,',',m,') =',cs(l,m) + end if + end do + end do + + end if + + if (logfile) then + + write(1,*) + + do l=1,lmax + do m=1,l + if (stepz*mod(l-m,2).eq.0) then + write(1,"(A4,I2,A1,I2,A4,ES14.6)") ' cc(', + . l,',',m,') =',cc(l,m) + end if + end do + end do + write(1,*) + do l=1,lmax + do m=1,l + if (stepz*mod(l-m,2).eq.0) then + write(1,"(A4,I2,A1,I2,A4,ES14.6)") ' cs(', + . l,',',m,') =',cs(l,m) + end if + end do + end do + + end if + + 10 format(I4,I4,A2,ES14.6,A1,ES14.6) + + else + + if (veryver) then + + write(*,*) + + do l=1,lmax + do m=1,l + if (stepz*mod(l-m,2).eq.0) then + write(*,"(A4,I2,A1,I2,A4,ES14.6)") ' cc(', + . l,',',m,') =',cc(l,m) + end if + end do + end do + + end if + + if (logfile) then + + write(1,*) + + do l=1,lmax + do m=1,l + if (stepz*mod(l-m,2).eq.0) then + write(1,"(A4,I2,A1,I2,A4,ES14.6)") ' cc(', + . l,',',m,') =',cc(l,m) + end if + end do + end do + + end if + + 20 format(I4,I4,A2,ES14.6) + + end if + end if + end if + +! Close logfile. + + if (logfile.and.(myproc.eq.0)) then + close(1) + end if + + +! *********************************** +! *** CHECK IF WE HAVE FAILED *** +! *********************************** + +! The finder needs a way to give up if there is no +! horizon. At the moment, the way in which I do this +! is that I just give up if the monopole term is of +! the order of 2 grid points in size. Maybe this +! could be improved on? + + if (c0(0).lt.2.0D0*dd) then + write(*,*) + write(*,"(A43)") ' AHFinder: c0(0) is too small, giving up.' + return + end if + + +! ******************************************** +! *** CHECK IF WE HAVE FOUND A HORIZON *** +! ******************************************** + +! For this I use a straigthforward test to check if the +! change in the surface is below a given tolerance. This +! usually works well, but if the tolerance is too small then +! the convergence becomes very slow. Also, since I use an +! adaptive stepsize, I run the risk of having a small change +! simply because the stepsize itself is very small, and not +! because we found a horizon. This can happen is the tolerance +! is not small enough. So we have a balancing act here. + + if (ahfdiff.lt.ahftol*ahfsum) then + status = .true. + return + end if + + end do + + +! ****************************** +! *** TOO MANY IERATIONS *** +! ****************************** + +! If we ever get here, we did too many iterations and +! failed to converge. + + status = .true. + + if (veryver) then + write(*,*) + write(*,"(A32,I2)") ' AHFinder: Too many iterations' + end if + + +! *************** +! *** END *** +! *************** + + end subroutine AHFinder_flow + + + + + + + + + + + + + real*8 function spher(theta,phi) + +! Function to reconstruct a function pointwise +! from the expansion coefficients {c0,cc,cs}. + + use AHFinder_dat + + implicit none + + integer l,m + + REAL LEGEN + REAL theta,phi + REAL F,cost,cosa,sina + REAL aux + + +! ****************************** +! *** Axisymmetric terms *** +! ****************************** + + F = 0 + + cost = cos(theta) + + do l=0,lmax + F = F + c0(l)*legen(l,0,cost) + end do + + +! ********************************** +! *** Non-axisymmetric terms *** +! ********************************** + + do m=1,lmax + + aux = m*phi + sina = sin(aux) + cosa = cos(aux) + + do l=m,lmax + F = F + LEGEN(l,m,cost)*(cc(l,m)*cosa + cs(l,m)*sina) + end do + + end do + + spher = F + + +! *************** +! *** END *** +! *************** + + end function spher diff --git a/src/AHFinder_fun.F b/src/AHFinder_fun.F new file mode 100644 index 0000000..d5261db --- /dev/null +++ b/src/AHFinder_fun.F @@ -0,0 +1,166 @@ +c/*@@ +c @file AHFinder_fun.F +c @date April 1998 +c @author Miguel Alcubierre +c @desc +c Find horizon function. +c @enddesc +c@@*/ + +c Note that including cactus.h will also include AHFinder.h +!#include "cactus.h" +#include "cctk.h" +#include "cctk_parameters.h" +#include "cctk_arguments.h" + + subroutine AHFinder_fun(CCTK_FARGUMENTS) + + use AHFinder_dat + + implicit none + + DECLARE_CCTK_FARGUMENTS + DECLARE_CCTK_PARAMETERS + + integer i,j,k + integer l,m + + CCTK_REAL LEGEN + CCTK_REAL xp,yp,zp,rp + CCTK_REAL phi,cost,cosa,sina + CCTK_REAL zero,half,one,two + CCTK_REAL pi,halfpi,twopi + CCTK_REAL aux1,aux2 + + + + +! ************************** +! *** DEFINE NUMBERS *** +! ************************** + + zero = 0.0D0 + half = 0.5D0 + one = 1.0D0 + two = 2.0D0 + + nx = cctk_lsh(1) + ny = cctk_lsh(2) + nz = cctk_lsh(3) + +! ************************************* +! *** DEFINE {pi,halfpi,twopi} *** +! ************************************* + + pi = 3.141592654D0 + + halfpi = half*pi + twopi = two*pi + + +! ********************************* +! *** FIND HORIZON FUNCTION *** +! ********************************* + + do k=1,nz + do j=1,ny + do i=1,nx + + xp = x(i,j,k) - xc + yp = y(i,j,k) - yc + zp = z(i,j,k) - zc + + rp = dsqrt(xp**2 + yp**2 + zp**2) + +! Find phi. + + if (dabs(xp).gt.dabs(yp)) then + phi = atan(dabs(yp/xp)) + else if (dabs(xp).lt.dabs(yp)) then + phi = halfpi - atan(dabs(xp/yp)) + else + phi = 0.25D0*pi + end if + + if ((xp.eq.zero).and.(yp.eq.zero)) then + phi = zero + else if ((xp.le.zero).and.(yp.ge.zero)) then + phi = pi - phi + else if ((xp.le.zero).and.(yp.le.zero)) then + phi = pi + phi + else if ((xp.ge.zero).and.(yp.le.zero)) then + phi = twopi - phi + end if + +! Monopole term. + + aux1 = c0(0) + +! Axisymmetric terms. + + if (rp.ne.zero) then + cost = zp/rp + else + cost = one + end if + + do l=1+stepz,lmax,1+stepz + aux1 = aux1 + c0(l)*LEGEN(l,0,cost) + end do + + +! Non-axisymmetric terms. Notice how the sum over M is first. +! This will allow me to use the recursion relations to avoid +! having to start from scratch every time. Also, I sum over +! all l's even if I don't want some terms. This is because +! in order to use the recursion relations I need all polynomials. + + if (nonaxi) then + do m=1,lmax + aux2 = dble(m)*phi + + sina = sin(aux2) + cosa = cos(aux2) + do l=m,lmax + aux2 = LEGEN(l,m,cost) + aux1 = aux1 + aux2*cc(l,m)*cosa + if (.not.refy) then + aux1 = aux1 + aux2*cs(l,m)*sina + end if + end do + end do + end if + + ahfgrid(i,j,k) = rp - aux1 + +! Check for NaN's. + + if (ahfgrid(i,j,k).ge.zero) then + else if (ahfgrid(i,j,k).lt.zero) then + else +#ifdef MPI + if (myproc.eq.0) then +#endif + write(*,*) 'NaN in ahfgrid at point',i,j,k +#ifdef MPI + end if +#endif + STOP + end if + + end do + end do + end do + +! Synchronize. + +#ifdef MPI + call synconefunc(ahfgrid) +#endif + +! *************** +! *** END *** +! *************** + + end subroutine AHFinder_fun + diff --git a/src/AHFinder_gau.F b/src/AHFinder_gau.F new file mode 100644 index 0000000..ef3e6cb --- /dev/null +++ b/src/AHFinder_gau.F @@ -0,0 +1,1168 @@ +c/*@@ +c @file AHFinder_gau.F +c @date October 1998 +c @author Miguel Alcubierre +c @desc +c Find gaussian curvature of surface. The gaussian +c curvature is defined as -R/2, where R is the Ricci +c scalar of the induced 2-geometry of the surface. +c +c As an extra `goody', this routine also integrates +c the equatorial circumference of the horizon and +c two polar circumferences at phi=0 and phi=pi/2. +c @enddesc +c@@*/ + +c Note that including cactus.h will also include AHFinder.h +!#include "cactus.h" +#include "cctk.h" +#include "cctk_parameters.h" +#include "cctk_arguments.h" + + subroutine AHFinder_gau(CCTK_FARGUMENTS,circ_eq,meri_p1,meri_p2) + + use AHFinder_dat + + implicit none + + DECLARE_CCTK_FARGUMENTS + DECLARE_CCTK_PARAMETERS + +#ifdef MPI + include 'mpif.h' +#endif + + logical firstgau + logical firstcal(4) + + integer handle,x_index,y_index,z_index + CCTK_INT ierror + CCTK_REAL maxval(3),minval(3) + + + CCTK_INT i,j,k,l,m,n,p + CCTK_INT npoints + CCTK_INT error1,error2,ierror + CCTK_INT getoutpfx + + CCTK_REAL LEGEN + CCTK_REAL theta,phi,xp,yp,zp,rp + CCTK_REAL cost,sint,cosp,sinp + CCTK_REAL dtheta,dphi,dtp,idtheta,idphi + CCTK_REAL xmn,ymn,zmn,xmx,ymx,zmx + CCTK_REAL trr,ttt,tpp,trt,trp,ttp,ft,fp + CCTK_REAL circ_eq,meri_p1,meri_p2 + CCTK_REAL deta,ideta + CCTK_REAL nlx,nly,nlz,nux,nuy,nuz + CCTK_REAL trxi,xi2 + CCTK_REAL zero,half,one,two,three,four,pi + CCTK_REAL aux,sina,cosa + + CCTK_REAL, dimension(3,3) :: ug,xi + CCTK_REAL, dimension(2,2) :: ga,ua + CCTK_REAL, dimension(2,2,2) :: d1a,gammad,gamma + CCTK_REAL, dimension(2,2,2,2) :: d2a,gamma2 + + CCTK_REAL, allocatable, dimension(:,:) :: rr,xa,ya,za + CCTK_REAL, allocatable, dimension(:,:) :: txx,tyy,tzz,txy,txz,tyz + CCTK_REAL, allocatable, dimension(:,:) :: g11,g22,g12 + CCTK_REAL, allocatable, dimension(:,:) :: gaussian + + character*200 gaussf + +! Declarations for macros. + +!#include "./thorn_StandardEinstein/src/macro/UPPERMET_declare.h" +!#include "./thorn_StandardEinstein/src/macro/CHR2_declare.h" +!#include "./thorn_StandardEinstein/src/macro/RICCI_declare.h" +!#include "./thorn_StandardEinstein/src/macro/TRRICCI_declare.h" + +#include "../../Einstein/src/macro/UPPERMET_declare.h" +#include "../../Einstein/src/macro/CHR2_declare.h" +#include "../../Einstein/src/macro/RICCI_declare.h" +#include "../../Einstein/src/macro/TRRICCI_declare.h" + +! Data. + + data firstgau / .true. / + data firstcal(1) / .true. / + data firstcal(2) / .true. / + data firstcal(3) / .true. / + data firstcal(4) / .true. / + save firstcal + save firstgau + +! Description of variables: +! +! i,j,k, +! l,m,n,p Counters. +! +! npoints Number of points to interpolate. +! +! error1 Different from zero if radius is negative. +! error2 Different from zero if outside computational domain. +! +! theta Latitude. +! phi Longitude. +! +! cost cos(theta) +! sint sin(theta) +! cosp cos(phi) +! sinp sin(phi) +! +! dtheta Grid spacing in theta. +! dphi Grid spacing in phi. +! dtp dtheta*dphi. +! +! idtheta 1/(2 dtheta) +! idphi 1/(2 dphi) +! +! xp x coordinate from surface centre. +! yp y coordinate from surface centre. +! zp z coordinate from surface centre. +! +! rp Radius. +! rr Radius array. +! +! xa Array with x coordinates from grid centre. +! ya Array with y coordinates from grid centre. +! za Array with z coordinates from grid centre. +! +! xmn Minimum value of x over the grid. +! xmx Maximum value of x over the grid. +! +! ymn Minimum value of y over the grid. +! ymx Maximum value of y over the grid. +! +! zmn Minimum value of z over the grid. +! zmx Maximum value of z over the grid. +! +! txx Interpolated gxx metric component. +! tyy Interpolated gyy metric component. +! tzz Interpolated gzz metric component. +! txy Interpolated gxy metric component. +! txz Interpolated gxz metric component. +! tyz Interpolated gyz metric component. +! +! trr 3-metric component {r,r}. +! ttt 3-metric component {theta,theta}. +! tpp 3-metric component {phi,phi}. +! trt 3-metric component {r,theta}. +! trp 3-metric component {r,phi}. +! ttp 3-metric component {theta,phi}. +! +! ft dr/dtheta +! fp dr/dphi +! +! g11 2-metric component {theta,theta}. +! g22 2-metric component {phi,phi}. +! g12 2-metric component {theta,phi}. +! +! circ_eq Equatorial circumference. +! meri_p1 Length of meridian at phi=0. +! meri_p2 Length of meridian at phi=pi/2. +! +! gaussian Gaussian curvature. +! +! ug 3x3 array with inverse 3-metric. +! +! xi_ij 3x3 array for extrinsic curvature of level sets +! of horizon function. +! trxi Trace of xi. +! xi2 Square of xi. +! +! ga_ij 2x2 array with angular metric, +! ua_ij 2x2 array with inverse angular metric. +! +! deta Determinant of angular metric. +! ideta 1/deta +! +! d1a_ijk d ga +! i kl +! +! d2a_ijkl d d ga +! i j kl +! +! gammad Christoffel symbol with 3 indices down. +! +! gamma Christoffel symbol with first index up. +! +! gamma2 Product of Christoffel symbols. + + +! ************************** +! *** DEFINE NUMBERS *** +! ************************** + + zero = 0.0D0 + half = 0.5D0 + one = 1.0D0 + two = 2.0D0 + three = 3.0D0 + four = 4.0D0 + + pi = 3.141592654D0 + + +! ************************************** +! *** ALLOCATE MEMORY FOR ARRAYS *** +! ************************************** + + if (myproc.eq.0) then + + allocate(rr(0:ntheta,0:nphi)) + + allocate(xa(0:ntheta,0:nphi)) + allocate(ya(0:ntheta,0:nphi)) + allocate(za(0:ntheta,0:nphi)) + + allocate(txx(0:ntheta,0:nphi)) + allocate(tyy(0:ntheta,0:nphi)) + allocate(tzz(0:ntheta,0:nphi)) + allocate(txy(0:ntheta,0:nphi)) + allocate(txz(0:ntheta,0:nphi)) + allocate(tyz(0:ntheta,0:nphi)) + + allocate(g11(0:ntheta,0:nphi)) + allocate(g22(0:ntheta,0:nphi)) + allocate(g12(0:ntheta,0:nphi)) + + allocate(gaussian(0:ntheta,0:nphi)) + + else + + allocate(rr(0:0,0:0)) + + allocate(xa(0:0,0:0)) + allocate(ya(0:0,0:0)) + allocate(za(0:0,0:0)) + + allocate(txx(0:0,0:0)) + allocate(tyy(0:0,0:0)) + allocate(tzz(0:0,0:0)) + allocate(txy(0:0,0:0)) + allocate(txz(0:0,0:0)) + allocate(tyz(0:0,0:0)) + + allocate(g11(0:0,0:0)) + allocate(g22(0:0,0:0)) + allocate(g12(0:0,0:0)) + + allocate(gaussian(0:0,0:0)) + + end if + + +! ********************************* +! *** INITIALIZE PARAMETERS *** +! ********************************* + +! Find grid boundaries. +! call CCTK_ReductionHandle(handle,"minimum") +! call CCTK_GetCoordIndex(x_index,"x") +! call CCTK_GetCoordIndex(y_index,"y") +! call CCTK_GetCoordIndex(z_index,"z") + +! call CCTK_Reduce(cctkGH,ierror,-1,handle,3, +! . CCTK_VARIABLE_REAL,minval,3, +! . x_index,y_index,z_index) + +! xmn = gf_min(x) +! ymn = gf_min(y) +! zmn = gf_min(z) + minval = -9.45D0 + + xmn = minval(1) + ymn = minval(2) + zmn = minval(3) + +! call CCTK_ReductionHandle(handle,"maximum") +! call CCTK_GetCoordIndex(x_index,"x") +! call CCTK_GetCoordIndex(y_index,"y") +! call CCTK_GetCoordIndex(z_index,"z") + +! call CCTK_Reduce(cctkGH,ierror,-1,handle,3, +! . CCTK_VARIABLE_REAL,maxval,3, +! . x_index,y_index,z_index) + +! xmx = gf_max(x) +! ymx = gf_max(y) +! zmx = gf_max(z) + maxval =9.45D0 + + xmx = maxval(1) + ymx = maxval(2) + zmx = maxval(3) + +! Find {dtheta,dphi,dtp}. + + dtheta = pi/dble(ntheta) + + if (cartoon) then + + dphi = zero + dtp = dtheta + idphi = one + idtheta = half/dtheta + + else + + dphi = two*pi/dble(nphi) + + if (refz) dtheta = half*dtheta + + if (refx.and.refy) then + dphi = 0.25D0*dphi + else if (refx) then + dphi = half*dphi + print *, 'This combination of symmetries has not' + print *, 'been properly implemented yet!' + else if (refy) then + dphi = half*dphi + end if + + dtp = dtheta*dphi + + idtheta = half/dtheta + idphi = half/dphi + + end if + + +! ******************************************************* +! *** FIND GAUSSIAN CURVATURE AS 3D GRID FUNCTION *** +! ******************************************************* + +! Here I calculate the gaussian curvature of the level sets +! of the horizon function as a 3D grid function. Later I will +! interpolate it onto the horizon. This is cleaner than my old +! method (see below). It depends on less interpolations, and +! the interpolations are done only after all derivatives have +! been calculated. + + do k=2,nz-1 + do j=2,ny-1 + do i=2,nx-1 + +! Find covariant normal unit vector. + + aux = one/ahfgradn(i,j,k) + + nlx = ahfgradx(i,j,k)*aux + nly = ahfgrady(i,j,k)*aux + nlz = ahfgradz(i,j,k)*aux + +! Find upper spatial metric using standard macro. + +#include "../../Einstein/src/macro/UPPERMET_guts.h" + + ug(1,1) = UPPERMET_UXX + ug(2,2) = UPPERMET_UYY + ug(3,3) = UPPERMET_UZZ + + ug(1,2) = UPPERMET_UXY + ug(1,3) = UPPERMET_UXZ + ug(2,3) = UPPERMET_UYZ + + ug(2,1) = ug(1,2) + ug(3,1) = ug(1,3) + ug(3,2) = ug(2,3) + +! Find contravariant normal unit vector. + + nux = ug(1,1)*nlx + ug(1,2)*nly + ug(1,3)*nlz + nuy = ug(2,1)*nlx + ug(2,2)*nly + ug(2,3)*nlz + nuz = ug(3,1)*nlx + ug(3,2)*nly + ug(3,3)*nlz + +! Find Christoffel symbols using standard macros. + +#include "../../Einstein/src/macro/CHR2_guts.h" + +! Find extrinsic curvature of the 2-surfaces defined +! as the level sets of the horizon function. The +! extrinsic curvature is defined as: +! __ +! xi = - \/ n +! ij (i j) +! __ +! where \/ is the 3-dimensional covariant derivative and +! n_i is the unit normal vector to the level sets of the +! horizon function. + + xi(1,1) = - (ahfgradx(i+1,j,k)/ahfgradn(i+1,j,k) + . - ahfgradx(i-1,j,k)/ahfgradn(i-1,j,k))/(two*dx) + . + (CHR2_XXX*nlx + CHR2_YXX*nly + CHR2_ZXX*nlz) + + xi(2,2) = - (ahfgrady(i,j+1,k)/ahfgradn(i,j+1,k) + . - ahfgrady(i,j-1,k)/ahfgradn(i,j-1,k))/(two*dy) + . + (CHR2_XYY*nlx + CHR2_YYY*nly + CHR2_ZYY*nlz) + + xi(3,3) = - (ahfgradz(i,j,k+1)/ahfgradn(i,j,k+1) + . - ahfgradz(i,j,k-1)/ahfgradn(i,j,k-1))/(two*dz) + . + (CHR2_XZZ*nlx + CHR2_YZZ*nly + CHR2_ZZZ*nlz) + + xi(1,2) = - (ahfgradx(i,j+1,k)/ahfgradn(i,j+1,k) + . - ahfgradx(i,j-1,k)/ahfgradn(i,j-1,k))/(four*dy) + . - (ahfgrady(i+1,j,k)/ahfgradn(i+1,j,k) + . - ahfgrady(i-1,j,k)/ahfgradn(i-1,j,k))/(four*dx) + . + (CHR2_XXY*nlx + CHR2_YXY*nly + CHR2_ZXY*nlz) + + xi(1,3) = - (ahfgradx(i,j,k+1)/ahfgradn(i,j,k+1) + . - ahfgradx(i,j,k-1)/ahfgradn(i,j,k-1))/(four*dz) + . - (ahfgradz(i+1,j,k)/ahfgradn(i+1,j,k) + . - ahfgradz(i-1,j,k)/ahfgradn(i-1,j,k))/(four*dx) + . + (CHR2_XXZ*nlx + CHR2_YXZ*nly + CHR2_ZXZ*nlz) + + xi(2,3) = - (ahfgrady(i,j,k+1)/ahfgradn(i,j,k+1) + . - ahfgrady(i,j,k-1)/ahfgradn(i,j,k-1))/(four*dz) + . - (ahfgradz(i,j+1,k)/ahfgradn(i,j+1,k) + . - ahfgradz(i,j-1,k)/ahfgradn(i,j-1,k))/(four*dy) + . + (CHR2_XYZ*nlx + CHR2_YYZ*nly + CHR2_ZYZ*nlz) + + xi(2,1) = xi(1,2) + xi(3,1) = xi(1,3) + xi(3,2) = xi(2,3) + +! Find trace of xi. + + aux = zero + + do l=1,3 + do m=1,3 + aux = aux + ug(l,m)*xi(l,m) + end do + end do + + trxi = aux + +! Find square of xi. + + aux = zero + + do l=1,3 + do m=1,3 + do n=1,3 + do p=1,3 + aux = aux + ug(l,m)*ug(n,p)*xi(l,n)*xi(m,p) + end do + end do + end do + end do + + xi2 = aux + +! Find 3-Ricci and 3-Ricci scalar using standard macros. + +#include "../../Einstein/src/macro/RICCI_guts.h" +#include "../../Einstein/src/macro/TRRICCI_guts.h" + +! Find 2-Ricci scalar using the contracted Gauss-Codazzi +! relations: +! +! (2) (3) a b (3) 2 ab +! R = R - 2 n n R + (tr xi) - xi xi +! ab ab + + aux = TRRICCI_TRRICCI + trxi**2 - xi2 + + aux = aux - two*(nux**2*RICCI_RXX + nuy**2*RICCI_RYY + . + nuz**2*RICCI_RZZ + two*(nux*nuy*RICCI_RXY + . + nux*nuz*RICCI_RXZ + nuy*nuz*RICCI_RYZ)) + +! Find gaussian curvature. The gaussian curvature is +! defined as R/2, with R the Ricci scalar of the surfaces. + + ahfgauss(i,j,k) = half*aux + +! Undefine macros! + +#include "../../Einstein/src/macro/UPPERMET_undefine.h" +#include "../../Einstein/src/macro/CHR2_undefine.h" +#include "../../Einstein/src/macro/RICCI_undefine.h" +#include "../../Einstein/src/macro/TRRICCI_undefine.h" + + end do + end do + end do + +! Boundaries. + + ahfgauss(1,:,:) = ahfgauss(2,:,:) + ahfgauss(:,1,:) = ahfgauss(:,2,:) + ahfgauss(:,:,1) = ahfgauss(:,:,2) + + ahfgauss(nx,:,:) = ahfgauss(nx-1,:,:) + ahfgauss(:,ny,:) = ahfgauss(:,ny-1,:) + ahfgauss(:,:,nz) = ahfgauss(:,:,nz-1) + +! Synchronize. + +#ifdef MPI + call synconefunc(ahfgauss) +#endif + +! ************************************* +! *** FIND INTERPOLATING POINTS *** +! ************************************* + +! Initialize {error1,error2}. + + error1 = 0 + error2 = 0 + +! Notice that everything here is done on processor zero! + + if (myproc.eq.0) then + + do j=0,nphi + do i=0,ntheta + +! Find {theta,phi}. + + theta = dtheta*dble(i) + phi = dphi*dble(j) + +! Find sines and cosines. + + cost = cos(theta) + sint = sin(theta) + + cosp = cos(phi) + sinp = sin(phi) + +! Find radius rp. + + rp = c0(0) + + do l=1+stepz,lmax,1+stepz + rp = rp + c0(l)*LEGEN(l,0,cost) + end do + +! Notice how the sum over m is first. This will allow +! me to use the recursion relations to avoid having to +! start from scratch every time. Also, I sum over all +! l's even if I don't want some terms. This is +! because in order to use the recursion relations I +! need all polynomials. + + if (nonaxi) then + do m=1,lmax + aux = dble(m)*phi + sina = sin(aux) + cosa = cos(aux) + do l=m,lmax + aux = LEGEN(l,m,cost) + rp = rp + aux*cc(l,m)*cosa + if (.not.refy) then + rp = rp + aux*cs(l,m)*sina + end if + end do + end do + end if + +! Check for negative radius. + + if (rp.le.zero) then + error1 = 1 + end if + +! Save radius array. + + rr(i,j) = rp + +! Find {xa,ya,za} + + xa(i,j) = rp*sint*cosp + xc + ya(i,j) = rp*sint*sinp + yc + za(i,j) = rp*cost + zc + +! Check if we are within bounds. + + xp = xa(i,j) + yp = ya(i,j) + zp = za(i,j) + + if ((xp.gt.xmx).or.(xp.lt.xmn).or. + . (yp.gt.ymx).or.(yp.lt.ymn).or. + . (zp.gt.zmx).or.(zp.lt.zmn)) then + error2 = 1 + end if + + end do + end do + + else + + xa = half*(xmx+xmn) + ya = half*(ymx+ymn) + za = half*(zmx+zmn) + + end if + +! Reduce the errors across processors (all processors must +! know about this since all will participate on the interpolation +! below). + +#ifdef MPI + + call mpi_allreduce(error1,i,1,MPI_INTEGER, + . MPI_SUM,MPI_COMM_WORLD,ierror) + error1 = i + call mpi_allreduce(error2,i,1,MPI_INTEGER, + . MPI_SUM,MPI_COMM_WORLD,ierror) + error2 = i + +#endif + +! If there was an error then return from the subroutine +! (but remember to deallocate the arrays first!). + + if ((error1.ne.0).or.(error2.ne.0)) then + goto 10 + end if + + +! *********************** +! *** INTERPOLATE *** +! *********************** + + if (myproc.eq.0) then + npoints = (ntheta+1)*(nphi+1) + else + npoints = 1 + end if + + call start_getpointsgroup() + call getpoints(gxx,xa,ya,za,txx,npoints) + call getpoints(gyy,xa,ya,za,tyy,npoints) + call getpoints(gzz,xa,ya,za,tzz,npoints) + call getpoints(gxy,xa,ya,za,txy,npoints) + call getpoints(gxz,xa,ya,za,txz,npoints) + call getpoints(gyz,xa,ya,za,tyz,npoints) + call getpoints(ahfgauss,xa,ya,za,gaussian,npoints) + call cleanup_getpointsgroup() + + +! ********************************* +! *** FIND SPHERICAL METRIC *** +! ********************************* + + if (myproc.eq.zero) then + + do j=0,nphi + do i=0,ntheta + +! Find {theta,phi}. + + theta = dtheta*dble(i) + phi = dphi*dble(j) + +! Find {rp}. + + rp = rr(i,j) + +! Find sines and cosines. + + cost = cos(theta) + sint = sin(theta) + + cosp = cos(phi) + sinp = sin(phi) + +! Find metric in spherical coordinates. + + trr = sint**2*(txx(i,j)*cosp**2 + . + tyy(i,j)*sinp**2) + tzz(i,j)*cost**2 + . + two*sint*(txy(i,j)*sint*cosp*sinp + . + cost*(txz(i,j)*cosp + tyz(i,j)*sinp)) + + ttt = rp**2*(cost**2*(txx(i,j)*cosp**2 + . + tyy(i,j)*sinp**2) + tzz(i,j)*sint**2 + . + two*cost*(txy(i,j)*cost*cosp*sinp + . - sint*(txz(i,j)*cosp + tyz(i,j)*sinp))) + + tpp = (rp*sint)**2*(txx(i,j)*sinp**2 + tyy(i,j)*cosp**2 + . - two*txy(i,j)*cosp*sinp) + + trt = rp*(sint*cost*(txx(i,j)*cosp**2 + tyy(i,j)*sinp**2 + . - tzz(i,j) + two*txy(i,j)*sinp*cosp) + . + (cost**2-sint**2)*(txz(i,j)*cosp + tyz(i,j)*sinp)) + + trp = rp*sint*(sint*(cosp*sinp*(tyy(i,j) - txx(i,j)) + . + txy(i,j)*(cosp**2 - sinp**2)) + . + cost*(cosp*tyz(i,j) - sinp*txz(i,j))) + + ttp = rp**2*sint*(cost*(sinp*cosp*(tyy(i,j) - txx(i,j)) + . + (cosp**2 - sinp**2)*txy(i,j)) + . + sint*(sinp*txz(i,j) - cosp*tyz(i,j))) + +! Find derivatives {ft,fp}. For interior points I +! use centered differences, and for the boundary points +! I use second order one-sided differences. + + if ((i.ne.0).and.(i.ne.ntheta)) then + ft = idtheta*(rr(i+1,j) - rr(i-1,j)) + else if (i.eq.0) then + ft = - idtheta*(three*rr(0,j) + . - four*rr(1,j) + rr(2,j)) + else + ft = + idtheta*(three*rr(ntheta,j) + . - four*rr(ntheta-1,j) + rr(ntheta-2,j)) + end if + + if (nonaxi) then + if ((j.ne.0).and.(j.ne.nphi)) then + fp = idphi*(rr(i,j+1) - rr(i,j-1)) + else if (j.eq.0) then + fp = - idphi*(three*rr(i,0) + . - four*rr(i,1) + rr(i,2)) + else + fp = + idphi*(three*rr(i,nphi) + . - four*rr(i,nphi-1) + rr(i,nphi-2)) + end if + else + fp = zero + end if + +! Find induced metric on the surface. + + g11(i,j) = ttt + trr*ft**2 + two*trt*ft + g22(i,j) = tpp + trr*fp**2 + two*trp*fp + g12(i,j) = ttp + trr*ft*fp + trp*ft + trt*fp + + end do + end do + + +! ************************************ +! *** CALCULATE CIRCUMFERENCES *** +! ************************************ + +! Find length of equator and meridians + + circ_eq = zero + meri_p1 = zero + meri_p2 = zero + +! Equatorial circumference. + + if (refz) then + + do j=0,nphi-1 + circ_eq = circ_eq + dphi + . *dsqrt(half*(g22(ntheta,j) + g22(ntheta,j+1))) + end do + + else + + aux = half*pi/dtheta + i = int(aux) + aux = aux - int(aux) + + if (cartoon) then + circ_eq = 6.283185307D0*dsqrt(g22(i,0)) + else + do j=0,nphi-1 + circ_eq = circ_eq + dphi + . *((one-aux)*dsqrt(half*(g22(i,j) + g22(i,j+1))) + . + aux*dsqrt(half*(g22(i+1,j) + g22(i+1,j+1)))) + end do + end if + end if + +! Meridians + + do i=0,ntheta-1 + meri_p1 = meri_p1 + dtheta + . *dsqrt(half*(g11(i,0) + g11(i+1,0))) + end do + + if (cartoon) then + meri_p2 = meri_p1 + else if (refx.and.refy) then + do i=0,ntheta-1 + meri_p2 = meri_p2 + dtheta + . *dsqrt(half*(g11(i,nphi) + g11(i+1,nphi))) + end do + else + aux = half*pi/dphi + j = int(aux) + aux = aux - int(aux) + do i=0,ntheta-1 + meri_p2 = meri_p2 + dtheta + . *((one-aux)*dsqrt(half*(g11(i,j) + g11(i,j))) + . + aux*dsqrt(half*(g11(i,j+1) + g11(i,j+1)))) + end do + end if + +! Rescale the integrals according to symmetries. + + if (.not.cartoon) then + if (refx) circ_eq = two*circ_eq + if (refy) circ_eq = two*circ_eq + end if + + if (refz) then + meri_p1 = two*meri_p1 + meri_p2 = two*meri_p2 + end if + + end if + + +! ************************************************* +! *** OLD CALCULATION OF GAUSSIAN CURVATURE *** +! ************************************************* + +! This was the old calculation of the gaussian curvature. +! It is correct as far as I know, but it depends on taking +! second derivatives of the interpolated metric, and this +! seems not to behave well numerically. I didn't want to +! delete the code since it might be useful in the future, +! so I just jump over it. + + if (.false.) then + +! The Gaussian curvature of the surface is defined as R/2 +! where R is the Ricci scalar of the two-geometry. I can +! now calculate the Ricci scalar using the spherical metric +! components {ttt,tpp,ttp}. + +! Initialize arrays. + + gaussian = zero + + ga(1,1) = one + ga(2,2) = one + ga(1,2) = zero + ga(2,1) = zero + + ua(1,1) = one + ua(2,2) = one + ua(1,2) = zero + ua(2,1) = zero + + d1a = zero + d2a = zero + + gamma = zero + gamma2 = zero + gammad = zero + + if (myproc.eq.zero) then + +! Interior points. + + do j=1,nphi-1 + do i=1,ntheta-1 + +! Find angular metric. + + ga(1,1) = g11(i,j) + ga(2,2) = g22(i,j) + ga(1,2) = g12(i,j) + ga(2,1) = g12(i,j) + +! Find determinant of angular metric. + + deta = ga(1,1)*ga(2,2) - ga(1,2)**2 + ideta = one/deta + +! Find inverse angular metric. + + ua(1,1) = + ga(2,2)*ideta + ua(2,2) = + ga(1,1)*ideta + ua(1,2) = - ga(1,2)*ideta + ua(2,1) = - ga(2,1)*ideta + +! First derivatives of angular metric. + + d1a(1,1,1) = idtheta*(g11(i+1,j) - g11(i-1,j)) + d1a(1,2,2) = idtheta*(g22(i+1,j) - g22(i-1,j)) + d1a(1,1,2) = idtheta*(g12(i+1,j) - g12(i-1,j)) + d1a(1,2,1) = d1a(1,1,2) + + if (nonaxi) then + d1a(2,1,1) = idphi*(g11(i,j+1) - g11(i,j-1)) + d1a(2,2,2) = idphi*(g22(i,j+1) - g22(i,j-1)) + d1a(2,1,2) = idphi*(g12(i,j+1) - g12(i,j-1)) + d1a(2,2,1) = d1a(2,1,2) + else + d1a(2,1,1) = zero + d1a(2,2,2) = zero + d1a(2,1,2) = zero + d1a(2,2,1) = zero + end if + +! Second derivatives of angular metric. + + d2a(1,1,1,1) = four*idtheta**2*(g11(i+1,j) + . - two*g11(i,j) + g11(i-1,j)) + d2a(1,1,2,2) = four*idtheta**2*(g22(i+1,j) + . - two*g22(i,j) + g22(i-1,j)) + d2a(1,1,1,2) = four*idtheta**2*(g12(i+1,j) + . - two*g12(i,j) + g12(i-1,j)) + + if (nonaxi) then + d2a(2,2,1,1) = four*idphi**2*(g11(i,j+1) + . - two*g11(i,j) + g11(i,j-1)) + d2a(2,2,2,2) = four*idphi**2*(g22(i,j+1) + . - two*g22(i,j) + g22(i,j-1)) + d2a(2,2,1,2) = four*idphi**2*(g12(i,j+1) + . - two*g12(i,j) + g12(i,j-1)) + + d2a(1,2,1,1) = (g11(i+1,j+1) - g11(i+1,j-1) + . - g11(i-1,j+1) + g11(i-1,j-1)) + . *idtheta*idphi + d2a(1,2,2,2) = (g22(i+1,j+1) - g22(i+1,j-1) + . - g22(i-1,j+1) + g22(i-1,j-1)) + . *idtheta*idphi + d2a(1,2,1,2) = (g12(i+1,j+1) - g12(i+1,j-1) + . - g12(i-1,j+1) + g12(i-1,j-1)) + . *idtheta*idphi + else + d2a(2,2,1,1) = zero + d2a(2,2,2,2) = zero + d2a(2,2,1,2) = zero + d2a(1,2,1,1) = zero + d2a(1,2,2,2) = zero + d2a(1,2,1,2) = zero + end if + + d2a(1,1,2,1) = d2a(1,1,1,2) + d2a(2,2,2,1) = d2a(2,2,1,2) + d2a(1,2,2,1) = d2a(1,2,1,2) + + d2a(2,1,1,1) = d2a(1,2,1,1) + d2a(2,1,2,2) = d2a(1,2,2,2) + d2a(2,1,1,2) = d2a(1,2,1,2) + d2a(2,1,2,1) = d2a(1,2,2,1) + +! Add correction terms to derivatives of g_{phi,phi}. +! Near the axis, it turns out that the lower order +! contributions to the derivatives of g_{phi,phi} +! cancel out, and the gaussian curvature is made only +! of higher order terms that the centered finite difference +! approximations get wrong! So here I add higher order +! corrections to make sure that I get correctly the +! first two terms of the Taylor series (assuming that +! for small theta g_{phi,phi} = f(r,phi) sin(theta)^2). + + theta = dtheta*dble(i) + + d1a(1,2,2) = d1a(1,2,2) + (two/three)*sin(theta) + . *(g22(i+1,j) - two*g22(i,j) + g22(i-1,j)) + + d2a(1,1,2,2) = d2a(1,1,2,2) + one/three + . *(g22(i+1,j) - two*g22(i,j) + g22(i-1,j)) + + if (nonaxi) then + d2a(1,2,2,2) = d2a(1,2,2,2) + (two/three)*sin(theta) + . *((g22(i+1,j+1) - two*g22(i,j+1) + g22(i-1,j+1)) + . - (g22(i+1,j-1) - two*g22(i,j-1) + g22(i-1,j-1))) + . *idphi + d2a(2,1,2,2) = d2a(1,2,2,2) + end if + +! Find Christoffel symbols. + + do l=1,2 + do m=1,2 + do n=1,2 + gammad(l,m,n) = half*(d1a(m,n,l) + d1a(n,m,l) + . - d1a(l,m,n)) + end do + end do + end do + + do l=1,2 + do m=1,2 + do n=1,2 + aux = zero + do k=1,2 + aux = aux + ua(l,k)*gammad(k,m,n) + end do + gamma(l,m,n) = aux + end do + end do + end do + + do k=1,2 + do l=1,2 + do m=1,2 + do n=1,2 + aux = zero + do p=1,2 + aux = aux + gamma(p,k,l)*gammad(p,m,n) + end do + gamma2(k,l,m,n) = aux + end do + end do + end do + end do + +! Find gaussian curvature. + + aux = zero + + do k=1,2 + do l=1,2 + do m=1,2 + do n=1,2 + aux = aux + ua(k,l)*ua(m,n) + . *(d2a(k,n,l,m) - d2a(k,l,m,n) + . + gamma2(k,m,l,n) - gamma2(k,l,m,n)) + end do + end do + end do + end do + + gaussian(i,j) = half*aux + + end do + end do + +! Boundaries. + + gaussian(0,:) = two*gaussian(1,:) - gaussian(2,:) + gaussian(ntheta,:) = two*gaussian(ntheta-1,:) + . - gaussian(ntheta-2,:) + + gaussian(:,0) = two*gaussian(:,1) - gaussian(:,2) + gaussian(:,nphi) = two*gaussian(:,nphi-1) + . - gaussian(:,nphi-2) + + + end if + + end if + + +! ************************************ +! *** INITIALIZE FIRSTGAU FLAG *** +! ************************************ + + if (find3) then + firstgau = firstcal(mfind+1) + firstcal(mfind+1) = .false. + else + firstgau = firstcal(4) + firstcal(4) = .false. + end if + + +! *********************************** +! *** SAVE GAUSSIAN CURVATURE *** +! *********************************** + + if (myproc.eq.0) then + if (find3) then + if (mfind.eq.0) then + gaussf = filestr(1:nfile)//"mahf_0.gauss" + else if (mfind.eq.1) then + gaussf = filestr(1:nfile)//"mahf_1.gauss" + else if (mfind.eq.2) then + gaussf = filestr(1:nfile)//"mahf_2.gauss" + end if + else + gaussf = filestr(1:nfile)//"mahf.gauss" + end if + + +! On first call replace file. + + if (firstgau) then + + if (find3) then + if (mfind.eq.2) then + firstgau = .false. + end if + else + firstgau = .false. + end if + + open(1,file=gaussf,form='formatted', + . status='replace') + + write(1,"(A20)") '# GAUSSIAN CURVATURE' + write(1,"(A1)") '#' + write(1,"(A35)") '# The data is written in a loop as:' + write(1,"(A1)") '#' + write(1,"(A17)") '# do i=0,ntheta-1' + write(1,"(A18)") '# do j=0,nphi-1' + write(1,"(A27)") '# write gaussian(i,j)' + write(1,"(A11)") '# end do' + write(1,"(A8)") '# end do' + write(1,"(A1)") '#' + write(1,"(A40)") '# theta and phi are subdivided uniformly' + write(1,"(A26)") '# according to symmetries:' + write(1,"(A1)") '#' + write(1,"(A36)") '# phi=[0,2 pi] (refx=refy=.false.)' + write(1,"(A44)") '# phi=[0,pi] (refx=.false., refy=.true.)' + write(1,"(A35)") '# phi=[0,pi/2] (refx=refy=.true.)' + write(1,"(A1)") '#' + write(1,"(A31)") '# theta=[0,pi] (refz=.false.)' + write(1,"(A30)") '# theta=[0,pi/2] (refz=.true.)' + write(1,"(A1)") '#' + + write(1,"(A9,L1)") '# refx = ',refx + write(1,"(A9,L1)") '# refy = ',refy + write(1,"(A9,L1)") '# refz = ',refz + write(1,"(A1)") '#' + write(1,"(A11,I7)") '# ntheta = ',(ntheta+1) + write(1,"(A11,I7)") '# nphi = ',(nphi+1) + + write(1,*) + write(1,"(A11,I4)") '# Time step',cctk_iteration + write(1,"(A6,ES11.3)") '# Time',cctk_time + write(1,"(A6,I4)") '# Call',ahf_ncall + +! On subsequent calls, append to file. + + else + + open(1,file=gaussf,form='formatted', + . status='old',position='append') + + write(1,*) + write(1,"(A11,I4)") '# Time step',cctk_iteration + write(1,"(A6,ES11.3)") '# Time',cctk_time + write(1,"(A6,I4)") '# Call',ahf_ncall + + end if + +! Save data points. + + do i=0,ntheta + do j=0,nphi + write(1,"(ES11.3)") gaussian(i,j) + end do + end do + +! Close file. + + close(1) + + end if + + +! ***************************** +! *** DEALLOCATE MEMORY *** +! ***************************** + + 10 continue + + deallocate(rr) + deallocate(xa,ya,za) + deallocate(txx,tyy,tzz,txy,txz,tyz) + deallocate(g11,g22,g12) + deallocate(gaussian) + + +! *************** +! *** END *** +! *************** + + end subroutine AHFinder_gau + + diff --git a/src/AHFinder_initguess.F b/src/AHFinder_initguess.F new file mode 100644 index 0000000..539a453 --- /dev/null +++ b/src/AHFinder_initguess.F @@ -0,0 +1,272 @@ +c/*@@ +c @file AHFinder_initguess.F +c @date July 1999 +c @author Miguel Alcubierre +c @desc +c Find initial guess for minimization algorithm +c @enddesc +c@@*/ + +c Note that including cactus.h will also include AHFinder.h +!#include "cactus.h" +#include "cctk.h" +#include "cctk_parameters.h" +#include "cctk_arguments.h" + + subroutine AHFinder_initguess(CCTK_FARGUMENTS,rmx) + + use AHFinder_dat + + implicit none + + DECLARE_CCTK_FARGUMENTS + DECLARE_CCTK_PARAMETERS + + integer i,j + CCTK_INT nn0,nn2,min_i,min_j + + CCTK_REAL zero,one,cc0,cc2,rmx,aux,min_aa + CCTK_REAL, allocatable, dimension(:,:) :: aa,a0,a2 + +! Description of variables +! +! i,j Counters. +! +! aa Area array for initial guess. +! a0 c0(0) array for initial guess. +! a2 c0(2) array for initial guess. +! +! rmx Radius of largest sphere that fits in the grid. +! +! nn0 Number of subdivisions of c0(0) for initial guess. +! nn2 Number of subdivisions of c0(2) for initial guess. + + + +! ********************** +! *** INITIALIZE *** +! ********************** + +! define numbers + + zero = 0.0D0 + one = 1.0D0 + aux = rmx + +! Initialize arrays. + + c0 = zero + cc = zero + cs = zero + cc0 = zero + cc2 = zero + +! Get {nn0,nn2}. + +! nn0 = geti("ahf_nn0") +! nn2 = geti("ahf_nn2") + nn0 = ahf_nn0 + nn2 = ahf_nn0 + + if (cartoon) nn2 = 0 + +! Allocate storage. + + allocate(aa(0:nn0,0:nn2)) + allocate(a0(0:nn0,0:nn2)) + allocate(a2(0:nn0,0:nn2)) + +! **************************** +! *** FIND INITIAL GUESS *** +! **************************** + +! IMPORTANT: In order to choose a good initial guess, +! I map a section of the space {c0(0),c0(2)} looking +! for a local minimum. This is VERY expensive, so the +! parameter space is not very well convered. + + if (veryver.and.(myproc.eq.0)) then + write(*,*) + write(*,*) 'AHFinder: Choosing initial guess.' + write(*,*) 'Be patient, this will take a while.' + end if + + +! lmax >= 2 and good guess. + + if ((lmax.ge.2).and.(.not.sloppy)) then + + do j=0,nn2 + + c0(2) = (0.5D0*dble(j)/dble(nn2) - 0.25D0)*aux + + call AHFinder_fun(CCTK_FARGUMENTS) + call AHFinder_exp(CCTK_FARGUMENTS) + + do i=0,nn0 + + c0(0) = 4.0D0*dx + dble(i)/dble(nn0)*(aux - 4.0D0*dx) + + a0(i,j) = c0(0) + a2(i,j) = c0(2) + + call AHFinder_int(CCTK_FARGUMENTS) + + if (minarea) then + aa(i,j) = intarea + else + aa(i,j) = intexp2 + end if + if(guessverbose) then + write(*,"(A11,2ES14.6)") 'a0 and a2: ',a0(i,j),a2(i,j) + write(*,"(A8,ES14.6)") ' H : ',intexp + write(*,"(A8,ES14.6)") ' H^2: ',intexp2 + write(*,"(A13,ES14.6)") ' H/area : ',intexp/intarea + write(*,"(A13,ES14.6)") ' H^2/area: ',intexp2/intarea + end if + + end do + end do + + if ((.not.inner).or.(minarea)) then + + min_i = 1 + min_j = 1 + min_aa = 1.0d15 + do i=1,nn0-1 + do j=1,nn2-1 + if ((aa(i,j).gt.zero).and. + . (aa(i,j).lt.aa(i+1,j)).and. + . (aa(i,j).lt.aa(i-1,j)).and. + . (aa(i,j).lt.aa(i,j+1)).and. + . (aa(i,j).lt.aa(i,j-1))) then + cc0 = a0(i,j) + cc2 = a2(i,j) + endif + if (aa(i,j) .le. min_aa) then + min_i = i + min_j = j + min_aa = aa(i,j) + endif + end do + end do + + else + + min_i = nn0-1 + min_j = 1 + min_aa = 1.0d15 + do i=nn0-1,1,-1 + do j=1,nn2-1 + if ((aa(i,j).gt.zero).and. + . (aa(i,j).lt.aa(i+1,j)).and. + . (aa(i,j).lt.aa(i-1,j)).and. + . (aa(i,j).lt.aa(i,j+1)).and. + . (aa(i,j).lt.aa(i,j-1))) then + cc0 = a0(i,j) + cc2 = a2(i,j) + endif + if (aa(i,j) .le. min_aa) then + min_i = i + min_j = j + min_aa = aa(i,j) + endif + end do + end do + + end if + if (guess_absmin) then + cc0 = a0(min_i,min_j) + cc2 = a2(min_i,min_j) + endif + +! Set up initial guess. + + c0(0) = cc0 + c0(2) = cc2 + +! lmax < 2 or sloppy guess. + + else + + call AHFinder_fun(CCTK_FARGUMENTS) + call AHFinder_exp(CCTK_FARGUMENTS) + + do i=0,nn0 + + c0(0) = 4.0D0*dx + dble(i)/dble(nn0)*(aux - 4.0D0*dx) + + a0(i,0) = c0(0) + + call AHFinder_int(CCTK_FARGUMENTS) + + if (minarea) then + aa(i,0) = intarea + else + aa(i,0) = intexp2 + end if + if(guessverbose) then + write(*,"(A11,2ES14.6)") 'a0 and a2: ',a0(i,0),a2(i,0) + write(*,"(A8,ES14.6)") ' H : ',intexp + write(*,"(A8,ES14.6)") ' H^2: ',intexp2 + write(*,"(A13,ES14.6)") ' H/area : ',intexp/intarea + write(*,"(A13,ES14.6)") ' H^2/area: ',intexp2/intarea + end if + + end do + + if ((.not.inner).or.(minarea)) then + + do i=1,nn0-1 + if ((aa(i,0).gt.zero).and. + . (aa(i,0).lt.0.9D0*aa(i+1,0)).and. + . (aa(i,0).lt.0.9D0*aa(i-1,0))) then + cc0 = a0(i,0) + end if + end do + + else + + do i=nn0-1,1,-1 + if ((aa(i,0).gt.zero).and. + . (aa(i,0).lt.0.9D0*aa(i+1,0)).and. + . (aa(i,0).lt.0.9D0*aa(i-1,0))) then + cc0 = a0(i,0) + end if + end do + + end if + +! Set up initial guess. + + c0(0) = cc0 + + end if + +! If no local minimum was found, start with a large +! sphere anyway. + + if (c0(0).eq.zero) then + c0(0) = 0.8D0*rmx + c0(2) = zero + end if + + +! *************** +! *** END *** +! *************** + +! Print message. + + if (veryver.and.(myproc.eq.0)) then + write(*,*) + write(*,*) 'Initial guess found, pinning down the horizon.' + end if + +! Deallocate storage. + + deallocate(aa,a0,a2) + + + end subroutine AHFinder_initguess + diff --git a/src/AHFinder_int.F b/src/AHFinder_int.F new file mode 100644 index 0000000..d13210e --- /dev/null +++ b/src/AHFinder_int.F @@ -0,0 +1,1105 @@ +c/*@@ +c @file AHFinder_int.F +c @date April 1998 +c @author Miguel Alcubierre +c @desc +c Find area of surface and integral of expansion. +c This is done in parallel, but the number of +c processors is assumed to be either the square +c of an integer or a power of two (if this is not +c the case, I use less processors). +c @enddesc +c@@*/ + +c Note that including cactus.h will also include AHFinder.h +!#include "cactus.h" +#include "cctk.h" +#include "cctk_parameters.h" +#include "cctk_arguments.h" + + subroutine AHFinder_int(CCTK_FARGUMENTS) + + use AHFinder_dat + + implicit none + + DECLARE_CCTK_FARGUMENTS + DECLARE_CCTK_PARAMETERS + +#ifdef MPI + include 'mpif.h' +#endif + + integer handle,x_index,y_index,z_index,ierror + CCTK_REAL maxval(3),minval(3) + + + integer i,j + integer l,m + integer npt,npp + integer l_ntheta,l_nphi,theta0,phi0 + integer npoints + integer ierror + integer auxi + integer gxx_index,gyy_index,gzz_index,gxy_index,gxz_index,gyz_index + integer exp_index,gradn_index + + CCTK_REAL LEGEN + CCTK_REAL theta,phi,xp,yp,zp,rp + + CCTK_REAL xw,yw,zw + + CCTK_REAL cost,sint,cosp,sinp + CCTK_REAL trr,ttt,tpp,trt,trp,ttp,ft,fp + CCTK_REAL xmn,ymn,zmn,xmx,ymx,zmx + CCTK_REAL dtheta,dphi,dtp,idtheta,idphi + CCTK_REAL det,idet + CCTK_REAL zero,quarter,half,one,two,three,four,pi + CCTK_REAL sina,cosa,intw,grad,sigma,h0 + CCTK_REAL aux,aux1,aux2 + CCTK_REAL gupij(3,3) + + CCTK_REAL tempv(0:lmax),tempm(lmax,lmax) + + CCTK_REAL, allocatable, dimension(:,:) :: rr,xa,ya,za,da,exp,gradn + CCTK_REAL, allocatable, dimension(:,:) :: txx,tyy,tzz,txy,txz,tyz + + save xmn,ymn,zmn,xmx,ymx,zmx + save dtheta,dphi,dtp,idtheta,idphi + save npt,npp,l_ntheta,l_nphi,theta0,phi0 + + +! Description of variables: +! +! i,j,l,m Counters. +! +! theta Latitude. +! phi Longitude. +! +! npt number of processors in theta direction. +! npp number of processors in phi direction. +! +! l_ntheta Local number of grid points in theta direction. +! l_nphi Local number of grid points in phi direction. +! +! theta0 Local origin in theta direction. +! phi0 Local origin in phi direction. +! +! npoints Number of points to interpolate . +! +! xp x coordinate from surface centre. +! yp y coordinate from surface centre. +! zp z coordinate from surface centre. +! +! rp Radius. +! rr Radius array. +! +! xa Array with x coordinates from grid centre. +! ya Array with y coordinates from grid centre. +! za Array with z coordinates from grid centre. +! +! exp Interpolated expansion. +! +! txx Array with interpolated gxx metric component. +! tyy Array with interpolated gyy metric component. +! tzz Array with interpolated gzz metric component. +! txy Array with interpolated gxy metric component. +! txz Array with interpolated gxz metric component. +! tyz Array with interpolated gyz metric component. +! +! gradn Array with interpolated norm of gradient of horizon function. +! +! cost cos(theta) +! sint sin(theta) +! cosp cos(phi) +! sinp sin(phi) +! +! trr Metric component {r,r}. +! ttt Metric component {theta,theta}. +! tpp Metric component {phi,phi}. +! trt Metric component {r,theta}. +! trp Metric component {r,phi}. +! ttp Metric component {theta,phi}. +! +! ft dr/dtheta +! fp dr/dphi +! +! xmn Minimum value of x over the grid. +! xmx Maximum value of x over the grid. +! +! ymn Minimum value of y over the grid. +! ymx Maximum value of y over the grid. +! +! zmn Minimum value of z over the grid. +! zmx Maximum value of z over the grid. +! +! dtheta Grid spacing in theta. +! dphi Grid spacing in phi. +! dtp dtheta*dphi. +! +! idtheta 1/(2 dtheta) +! idphi 1/(2 dphi) +! +! da Area element. + + + +! ************************** +! *** DEFINE NUMBERS *** +! ************************** + + zero = 0.0D0 + quarter = 0.25D0 + half = 0.5D0 + one = 1.0D0 + two = 2.0D0 + three = 3.0D0 + four = 4.0D0 + + pi = 3.141592654D0 + + myproc = 0 + nprocs = 1 + +! ********************************* +! *** INITIALIZE PARAMETERS *** +! ********************************* + + if (firstint) then + + firstint = .false. + +! Find grid boundaries. + +! call CCTK_ReductionHandle(handle,"minimum") +! call CCTK_GetCoordIndex(x_index,"x") +! call CCTK_GetCoordIndex(y_index,"y") +! call CCTK_GetCoordIndex(z_index,"z") + +! call CCTK_Reduce(cctkGH,ierror,-1,handle,3, +! . CCTK_VARIABLE_REAL,minval,3, +! . x_index,y_index,z_index) + +! xmn = gf_min(x) +! ymn = gf_min(y) +! zmn = gf_min(z) + minval = -9.45D0 + + xmn = minval(1) + ymn = minval(2) + zmn = minval(3) + +! call CCTK_ReductionHandle(handle,"maximum") +! call CCTK_GetCoordIndex(x_index,"x") +! call CCTK_GetCoordIndex(y_index,"y") +! call CCTK_GetCoordIndex(z_index,"z") + +! call CCTK_Reduce(cctkGH,ierror,-1,handle,3, +! . CCTK_VARIABLE_REAL,maxval,3, +! . x_index,y_index,z_index) + +! xmx = gf_max(x) +! ymx = gf_max(y) +! zmx = gf_max(z) + maxval = 9.45D0 + + xmx = maxval(1) + ymx = maxval(2) + zmx = maxval(3) + +! Find {dtheta,dphi,dtp}. + + dtheta = pi/dble(ntheta) + + if (cartoon) then + + dphi = zero + dtp = dtheta + idphi = one + idtheta = half/dtheta + + else + + dphi = two*pi/dble(nphi) + + if (refz) dtheta = half*dtheta + + if (refx.and.refy) then + dphi = quarter*dphi + else if (refx) then + print *, 'This combination of symmetries has not' + print *, 'been properly implemented yet!' + else if (refy) then + dphi = half*dphi + end if + + dtp = dtheta*dphi + + idtheta = half/dtheta + idphi = half/dphi + + end if + +! Check the number of processors to figure out +! how to divide the domain. + +! For cartoon this is trivial. + + + if (cartoon) then + + npt = nprocs + npp = 1 + +! Otherwise it is more complicated. At the moment I +! only consider numbers of processors that are either +! the square of an integer, or a power of two. +! If neither of this is the case, then I use fewer +! processors to make sure that I always deal with +! either a square or a power of two. This can +! certainly be improved, but it might not be that +! urgent since in most cases the number of processors +! will be a power of two. + + else + + aux = sqrt(dble(nprocs)) + + if (nprocs.eq.1) then + +! One processor. + + npt = 1 + npp = 1 + + else if (aux.eq.int(aux)) then + +! `nprocs' is the square of an integer. + + npt = int(aux) + npp = int(aux) + + else + +! Check if `nprocs' is a power of two. + + l = 1 + m = 1 + + i = 1 + + do while (l*m.lt.nprocs) + + npt = l + npp = m + + if (mod(i,2).ne.0) then + m = 2*m + else + l = 2*l + end if + + i = i + 1 + + end do + +! `nprocs' is a power of two. + + if (l*m.eq.nprocs) then + + npt = l + npp = m + +! `nprocs' is not a power of two. This is where I +! would need to do something clever. At the moment, +! I will just use fewer processors. The number of +! processors I use is the power of two or the square +! of an integer that is closest to `nprocs'. + + else + + if (npt*npp.lt.int(aux)**2) then + npt = int(aux) + npp = int(aux) + end if + + end if + + end if + + end if + +! Figure out the number of grid points per processor, +! and the `origin' for a given processor. + + if (myproc.ge.npt*npp) then + + l_ntheta = 0 + l_nphi = 0 + + theta0 = 0 + phi0 = 0 + + else + +! Take first the number of grid points per processor +! in a given direction to be equal to the total number +! of grid points divided by the number of processors +! in that direction. + + l_ntheta = ntheta/npt + l_nphi = nphi/npp + +! And take the `origin' as the corresponding displacement +! from the real origin. Notice that mod(myproc,npt) +! tells my on which `theta' bin I am, while myproc/npt +! tells me on which `phi' bin I am. + + i = mod(myproc,npt) + j = myproc/npt + + theta0 = i*l_ntheta + phi0 = j*l_nphi + +! Now correct all this if the total number of grid points +! in a given direction was not exactly divisible by the +! number of processors on that direction. + + if (l_ntheta*npt.ne.ntheta) then + +! Find residue on theta direction. + + l = ntheta - l_ntheta*npt + +! Distribute residue in first few processors, and +! displace theta0 accordingly. + + if (i.lt.l) then + l_ntheta = l_ntheta + 1 + theta0 = theta0 + i + else + theta0 = theta0 + l + end if + + end if + + if (l_nphi*npp.ne.nphi) then + +! Find residue on phi direction. + + l = nphi - l_nphi*npp + +! Distribute residue in first few processors, and +! displace phi0 accordingly. + + if (j.lt.l) then + l_nphi = l_nphi + 1 + phi0 = phi0 + j + else + phi0 = phi0 + l + end if + + end if + + end if + + end if + +! ************************************** +! *** ALLOCATE MEMORY FOR ARRAYS *** +! ************************************** + + allocate(rr(0:l_ntheta,0:l_nphi)) + + allocate(xa(0:l_ntheta,0:l_nphi)) + allocate(ya(0:l_ntheta,0:l_nphi)) + allocate(za(0:l_ntheta,0:l_nphi)) + + allocate(da(0:l_ntheta,0:l_nphi)) + allocate(exp(0:l_ntheta,0:l_nphi)) + allocate(gradn(0:l_ntheta,0:l_nphi)) + + allocate(txx(0:l_ntheta,0:l_nphi)) + allocate(tyy(0:l_ntheta,0:l_nphi)) + allocate(tzz(0:l_ntheta,0:l_nphi)) + allocate(txy(0:l_ntheta,0:l_nphi)) + allocate(txz(0:l_ntheta,0:l_nphi)) + allocate(tyz(0:l_ntheta,0:l_nphi)) + + +! *************************** +! *** START MAIN LOOP *** +! *************************** + +! Initialize {interror1,interror2}. + + interror1 = 0 + interror2 = 0 + +! Initialize {intexp,intexp2,intarea}. + + intexp = zero + intexp2 = zero + intexpdel2 = zero + intarea = zero + +! For flow algorithm, initialize spectral components. + + if (flow) then + + hflow0 = zero + hflowc = zero + hflows = zero + + cflow0 = zero + cflowc = zero + cflows = zero + + nflow0 = zero + nflowc = zero + nflows = zero + + end if + +! Find interpolating points. + + if (myproc.ge.npt*npp) then + + xa = half*(xmx+xmn) + ya = half*(ymx+ymn) + za = half*(zmx+zmn) + + else + + do j=0,l_nphi + do i=0,l_ntheta + +! Find {theta,phi}. + + theta = dtheta*dble(i+theta0) + phi = dphi*dble(j+phi0) + +! Find sines and cosines. + + cost = cos(theta) + sint = sin(theta) + + cosp = cos(phi) + sinp = sin(phi) + +! Find radius rp. + + rp = c0(0) + + do l=1+stepz,lmax,1+stepz + rp = rp + c0(l)*LEGEN(l,0,cost) + end do + +! Notice how the sum over m is first. This will allow +! me to use the recursion relations to avoid having to +! start from scratch every time. Also, I sum over all +! l's even if I don't want some terms. This is +! because in order to use the recursion relations I +! need all polynomials. + + if (nonaxi) then + do m=1,lmax + aux = dble(m)*phi + sina = sin(aux) + cosa = cos(aux) + do l=m,lmax + aux = LEGEN(l,m,cost) + rp = rp + aux*cc(l,m)*cosa + if (.not.refy) then + rp = rp + aux*cs(l,m)*sina + end if + end do + end do + end if + +! Check for negative radius. + + if (rp.le.zero) then + interror1 = 1 + end if + +! Save radius array. + + rr(i,j) = rp + +! Find {xa,ya,za} + + xa(i,j) = rp*sint*cosp + xc + ya(i,j) = rp*sint*sinp + yc + za(i,j) = rp*cost + zc + +! Check if we are within bounds. + + xp = xa(i,j) + yp = ya(i,j) + zp = za(i,j) + + if ((xp.gt.xmx).or.(xp.lt.xmn).or. + . (yp.gt.ymx).or.(yp.lt.ymn).or. + . (zp.gt.zmx).or.(zp.lt.zmn)) then + interror2 = 1 + end if + + end do + end do + + end if + +! Reduce the errors across processors. + +#ifdef MPI + + call mpi_allreduce(interror1,i,1,MPI_INTEGER, + . MPI_SUM,MPI_COMM_WORLD,ierror) + interror1 = i + call mpi_allreduce(interror2,i,1,MPI_INTEGER, + . MPI_SUM,MPI_COMM_WORLD,ierror) + interror2 = i + +#endif + +! If there was an error on any processor then assign +! large values to the integrals and return from the +! subroutine (but remember to deallocate the arrays +! first!). + + interror = .false. + + if ((interror1.ne.0).or.(interror2.ne.0)) then + intexp = 1.0D+10 + intexp2 = 1.0D+10 + intexpdel2 = 1.0D+10 + intarea = 1.0D+10 + interror = .true. + goto 10 + end if + +! Interpolate expansion and metric. + + npoints = (l_ntheta+1)*(l_nphi+1) + +! call start_getpointsgroup() +! call getpoints(ahf_exp,xa,ya,za,exp,npoints) +! call getpoints(ahfgradn,xa,ya,za,gradn,npoints) +! call getpoints(gxx,xa,ya,za,txx,npoints) +! call getpoints(gyy,xa,ya,za,tyy,npoints) +! call getpoints(gzz,xa,ya,za,tzz,npoints) +! call getpoints(gxy,xa,ya,za,txy,npoints) +! call getpoints(gxz,xa,ya,za,txz,npoints) +! call getpoints(gyz,xa,ya,za,tyz,npoints) +! call cleanup_getpointsgroup() + +! new interp + call CCTK_GetInterpHandle (handle, "simple") + call CCTK_VarIndex (gxx_index, "einstein::gxx") + call CCTK_VarIndex (gyy_index, "einstein::gyy") + call CCTK_VarIndex (gzz_index, "einstein::gzz") + call CCTK_VarIndex (gxy_index, "einstein::gxy") + call CCTK_VarIndex (gxz_index, "einstein::gxz") + call CCTK_VarIndex (gyz_index, "einstein::gyz") + call CCTK_VarIndex (exp_index, "ahfinder::ahf_exp") + call CCTK_VarIndex (gradn_index, "ahfinder::ahfgradn") + + call CCTK_InterpGF (cctkGH, ierror, handle, npoints, 3, 8, 8, + $ xa, ya, za, + $ CCTK_VARIABLE_REAL, CCTK_VARIABLE_REAL, + $ CCTK_VARIABLE_REAL, + $ gxx_index,gyy_index,gzz_index, + $ gxy_index,gxz_index,gyz_index, + $ exp_index,gradn_index, + $ txx,tyy,tzz,txy,txz,tyz,exp,gradn, + $ CCTK_VARIABLE_REAL,CCTK_VARIABLE_REAL,CCTK_VARIABLE_REAL, + $ CCTK_VARIABLE_REAL,CCTK_VARIABLE_REAL,CCTK_VARIABLE_REAL, + $ CCTK_VARIABLE_REAL,CCTK_VARIABLE_REAL ) + + +! Find 2D array for area element. + + if (myproc.ge.npt*npp) then + + da = zero + + else + + do j=0,l_nphi + do i=0,l_ntheta + +! Find {theta,phi}. + + theta = dtheta*dble(i+theta0) + phi = dphi*dble(j+phi0) + +! Find {rp}. + + rp = rr(i,j) + +! Find sines and cosines. + + cost = cos(theta) + sint = sin(theta) + + cosp = cos(phi) + sinp = sin(phi) + +! Find spherical metric components. + + trr = sint**2*(txx(i,j)*cosp**2 + . + tyy(i,j)*sinp**2) + tzz(i,j)*cost**2 + . + two*sint*(txy(i,j)*sint*cosp*sinp + . + cost*(txz(i,j)*cosp + tyz(i,j)*sinp)) + + ttt = rp**2*(cost**2*(txx(i,j)*cosp**2 + . + tyy(i,j)*sinp**2) + tzz(i,j)*sint**2 + . + two*cost*(txy(i,j)*cost*cosp*sinp + . - sint*(txz(i,j)*cosp + tyz(i,j)*sinp))) + + tpp = (rp*sint)**2*(txx(i,j)*sinp**2 + tyy(i,j)*cosp**2 + . - two*txy(i,j)*cosp*sinp) + + trt = rp*(sint*cost*(txx(i,j)*cosp**2 + tyy(i,j)*sinp**2 + . - tzz(i,j) + two*txy(i,j)*sinp*cosp) + . + (cost**2-sint**2)*(txz(i,j)*cosp + tyz(i,j)*sinp)) + + trp = rp*sint*(sint*(cosp*sinp*(tyy(i,j) - txx(i,j)) + . + txy(i,j)*(cosp**2 - sinp**2)) + . + cost*(cosp*tyz(i,j) - sinp*txz(i,j))) + + ttp = rp**2*sint*(cost*(sinp*cosp*(tyy(i,j) - txx(i,j)) + . + (cosp**2 - sinp**2)*txy(i,j)) + . + sint*(sinp*txz(i,j) - cosp*tyz(i,j))) + +! Find derivatives {ft,fp}. For interior points I +! use centered differences, and for the boundary points +! I use second order one sided differences. Notice that +! since this is done even in inter-processor boundaries, +! the final result will depend on the number of processors, +! but the differences will converge away at high order. + + if ((i.ne.0).and.(i.ne.l_ntheta)) then + ft = idtheta*(rr(i+1,j) - rr(i-1,j)) + else if (i.eq.0) then + ft = - idtheta*(three*rr(0,j) + . - four*rr(1,j) + rr(2,j)) + else + ft = + idtheta*(three*rr(l_ntheta,j) + . - four*rr(l_ntheta-1,j) + rr(l_ntheta-2,j)) + end if + + if (nonaxi) then + if ((j.ne.0).and.(j.ne.l_nphi)) then + fp = idphi*(rr(i,j+1) - rr(i,j-1)) + else if (j.eq.0) then + fp = - idphi*(three*rr(i,0) + . - four*rr(i,1) + rr(i,2)) + else + fp = + idphi*(three*rr(i,l_nphi) + . - four*rr(i,l_nphi-1) + rr(i,l_nphi-2)) + end if + else + fp = zero + end if + +! Find area element. For this I use the fact that if +! we have r = f(theta,phi), then the area element is: +! +! / 2 +! da = | [ gtt + grr (f,t) + 2 grt f,t ] +! \ +! 2 +! [ gpp + grr (f,p) + 2 grp f,p ] - [ gtp +! +! 2 \ 1/2 +! + grr f,t f,p + grt f,p + grp f,t ] | dtheta dphi +! / +! +! where f,t = df/dtheta and f,p = df/dphi. + + aux = (ttt + trr*ft**2 + two*trt*ft) + . * (tpp + trr*fp**2 + two*trp*fp) + . - (ttp + ft*(trr*fp + trp) + trt*fp)**2 + + if (aux.gt.zero) then + aux = sqrt(aux) + else + aux = zero + end if + +! Multiply with dtheta*dphi. + + da(i,j) = aux*dtp + + end do + end do + +! Find integrals. To find the integrals I use a second +! order method. I approximate the value of the integrand +! at each cell centre by an average of the values at the +! four cell corners, and then add up all cells. + + inside_min_count = 0.0d0 + inside_min_neg_count = 0.0d0 + + do j=0,l_nphi-1 + do i=0,l_ntheta-1 + + intarea = intarea + quarter + . *(da(i,j ) + da(i+1,j ) + . + da(i,j+1) + da(i+1,j+1)) + + intexp = intexp + quarter + . *(da(i ,j )*exp(i ,j ) + . + da(i+1,j )*exp(i+1,j ) + . + da(i ,j+1)*exp(i ,j+1) + . + da(i+1,j+1)*exp(i+1,j+1)) + + intexp2 = intexp2 + quarter + . *(da(i ,j )*exp(i ,j )**2 + . + da(i+1,j )*exp(i+1,j )**2 + . + da(i ,j+1)*exp(i ,j+1)**2 + . + da(i+1,j+1)*exp(i+1,j+1)**2) + + intexpdel2 = intexpdel2 + quarter + . *(da(i ,j )*(exp(i ,j )+trapped_surface_delta)**2 + . + da(i+1,j )*(exp(i+1,j )+trapped_surface_delta)**2 + . + da(i ,j+1)*(exp(i ,j+1)+trapped_surface_delta)**2 + . + da(i+1,j+1)*(exp(i+1,j+1)+trapped_surface_delta)**2) + + inside_min_count = inside_min_count + 1.0d0 + + if ((exp(i ,j ) .le. 0.0d0) .and. + . (exp(i+1,j ) .le. 0.0d0) .and. + . (exp(i ,j+1) .le. 0.0d0) .and. + . (exp(i+1,j+1) .le. 0.0d0)) then + inside_min_neg_count = inside_min_neg_count + 1.0d0 + endif + + end do + end do + +! If necessary, find spectral components for flow. +! For this I need to recalculate the Legendre polynomials, +! but I see no way around it. Also, to avoid having to +! define lots of 2D arrays I do the integrals inside +! the loop. This means that I have to take care to see +! if I am on an edge or corner point. Notice also that +! these integrals do not involve the area element. + + if (flow) then + +! Find h0. + + if (find_trapped_surface) then + h0 = - trapped_surface_delta + else + h0 = zero + end if + +! Find integrals. + + do j=0,l_nphi + do i=0,l_ntheta + +! Find {theta,phi}. + + theta = dtheta*dble(i+theta0) + phi = dphi*dble(j+phi0) + +! Find cost,sint. + + cost = cos(theta) + sint = sin(theta) + +! Find norm of gradient of horizon function. + + grad = gradn(i,j) + +! Find weight factor. + + intw = dtp*(exp(i,j)-h0)*sint**2 + + if (((j.eq.0).or.(j.eq.l_nphi)).and. + . ((i.eq.0).or.(i.eq.l_ntheta))) then + intw = quarter*intw + else if (((j.eq.0).or.(j.eq.l_nphi)).or. + . ((i.eq.0).or.(i.eq.l_ntheta))) then + intw = half*intw + end if + +! Find `sigma' for Nakamura flow (see Carsten's paper). + + if (nw.eq.0.0) then + sigma = zero + else + xw = xa(i,j)-xc + yw = ya(i,j)-yc + zw = za(i,j)-zc + +! Calculate metric on surface + det = txx(i,j)*tyy(i,j)*tzz(i,j) + $ + two*txy(i,j)*txz(i,j)*tyz(i,j) + $ - txx(i,j)*tyz(i,j)**2 - tyy(i,j)*txz(i,j)**2 + $ - tzz(i,j)*txy(i,j)**2 + + idet = one/det + + gupij(1,1) = idet*(tyy(i,j)*tzz(i,j)-tyz(i,j)**2) + gupij(2,2) = idet*(txx(i,j)*tzz(i,j)-txz(i,j)**2) + gupij(3,3) = idet*(txx(i,j)*tyy(i,j)-txy(i,j)**2) + + gupij(1,2) = idet*(txz(i,j)*tyz(i,j)-tzz(i,j)*txy(i,j)) + gupij(2,1) = gupij(1,2) + gupij(1,3) = idet*(txy(i,j)*tyz(i,j)-tyy(i,j)*txz(i,j)) + gupij(3,1) = gupij(1,3) + gupij(2,3) = idet*(txy(i,j)*txz(i,j)-txx(i,j)*tyz(i,j)) + gupij(3,2) = gupij(2,3) + + call AHFinder_calcsigma(CCTK_FARGUMENTS,xw,yw,zw,gupij,sigma) + + end if + +! Monopole integrals. + + aux = intw*grad + + hflow0(0) = hflow0(0) + intw + cflow0(0) = cflow0(0) + aux + nflow0(0) = nflow0(0) + aux*sigma + +! Axisymmetric integrals. + + do l=1+stepz,lmax,1+stepz + + aux = intw*LEGEN(l,0,cost) + aux1 = aux*grad + + hflow0(l) = hflow0(l) + aux + cflow0(l) = cflow0(l) + aux1 + nflow0(l) = nflow0(l) + aux1*sigma + + end do + +! Non-axisymmetric integrals. + + if (nonaxi) then + + do m=1,lmax + + aux = dble(m)*phi + sina = sin(aux) + cosa = cos(aux) + + do l=m,lmax + + aux = intw*LEGEN(l,m,cost) + + aux1 = aux*cosa + aux2 = aux1*grad + hflowc(l,m) = hflowc(l,m) + aux1 + cflowc(l,m) = cflowc(l,m) + aux2 + nflowc(l,m) = nflowc(l,m) + aux2*sigma + + if (.not.refy) then + aux1 = aux*sina + aux2 = aux1*grad + hflows(l,m) = hflows(l,m) + aux1 + cflows(l,m) = cflows(l,m) + aux2 + nflows(l,m) = nflows(l,m) + aux2*sigma + end if + + end do + end do + + end if + + end do + end do + + end if + + end if + +! Reduce the integrals across processors. + +#ifdef MPI + +! Area and expansion. + + call mpi_allreduce(intarea,aux,1,MPI_DOUBLE_PRECISION, + .MPI_SUM,MPI_COMM_WORLD,ierror) + intarea = aux + call mpi_allreduce(intexp,aux,1,MPI_DOUBLE_PRECISION, + .MPI_SUM,MPI_COMM_WORLD,ierror) + intexp = aux + call mpi_allreduce(intexp2,aux,1,MPI_DOUBLE_PRECISION, + .MPI_SUM,MPI_COMM_WORLD,ierror) + intexp2 = aux + call mpi_allreduce(intexpdel2,aux,1,MPI_DOUBLE_PRECISION, + .MPI_SUM,MPI_COMM_WORLD,ierror) + intexpdel2 = aux + +! negative expansion elements on surface + + call mpi_allreduce(inside_min_count,aux,1,MPI_DOUBLE_PRECISION, + .MPI_SUM,MPI_COMM_WORLD,ierror) + inside_min_count = aux + call mpi_allreduce(inside_min_neg_count,aux,1,MPI_DOUBLE_PRECISION, + .MPI_SUM,MPI_COMM_WORLD,ierror) + inside_min_neg_count = aux + +! Spectral components for flow. + + if (flow) then + +! Axisymmetric. + + auxi = lmax+1 + + if (hw.ne.zero) then + call mpi_allreduce(hflow0,tempv,auxi,MPI_DOUBLE_PRECISION, + . MPI_SUM,MPI_COMM_WORLD,ierror) + hflow0 = tempv + end if + if (cw.ne.zero) then + call mpi_allreduce(cflow0,tempv,auxi,MPI_DOUBLE_PRECISION, + . MPI_SUM,MPI_COMM_WORLD,ierror) + cflow0 = tempv + end if + if (nw.ne.zero) then + call mpi_allreduce(nflow0,tempv,auxi,MPI_DOUBLE_PRECISION, + . MPI_SUM,MPI_COMM_WORLD,ierror) + nflow0 = tempv + end if + +! Non-axisymmetric. + + if (nonaxi) then + + auxi = lmax*lmax + + if (nw.ne.zero) then + call mpi_allreduce(hflowc,tempm,auxi, + . MPI_DOUBLE_PRECISION,MPI_SUM,MPI_COMM_WORLD,ierror) + hflowc = tempm + end if + if (cw.ne.zero) then + call mpi_allreduce(cflowc,tempm,auxi, + . MPI_DOUBLE_PRECISION,MPI_SUM,MPI_COMM_WORLD,ierror) + cflowc = tempm + end if + if (nw.ne.zero) then + call mpi_allreduce(nflowc,tempm,auxi, + . MPI_DOUBLE_PRECISION,MPI_SUM,MPI_COMM_WORLD,ierror) + nflowc = tempm + end if + + if (.not.refy) then + + if (hw.ne.zero) then + call mpi_allreduce(hflows,tempm,auxi, + . MPI_DOUBLE_PRECISION,MPI_SUM,MPI_COMM_WORLD,ierror) + hflows = tempm + end if + if (cw.ne.zero) then + call mpi_allreduce(cflows,tempm,auxi, + . MPI_DOUBLE_PRECISION,MPI_SUM,MPI_COMM_WORLD,ierror) + cflows = tempm + end if + if (nw.ne.zero) then + call mpi_allreduce(nflows,tempm,auxi, + . MPI_DOUBLE_PRECISION,MPI_SUM,MPI_COMM_WORLD,ierror) + nflows = tempm + end if + end if + + end if + + end if + +#endif + +! For cartoon multiply the integrals with 2*pi. + + if (cartoon) then + + intw = 6.283185307D0 + + intexp = intw*intexp + intexp2 = intw*intexp2 + intexpdel2 = intw*intexpdel2 + intarea = intw*intarea + + if (flow) then + hflow0 = intw*hflow0 + cflow0 = intw*cflow0 + nflow0 = intw*nflow0 + end if + +! Rescale the integrals according to symmetries. + + else if (refx.or.refy.or.refz) then + + intw = one + + if (refx) intw = two*intw + if (refy) intw = two*intw + if (refz) intw = two*intw + + intexp = intw*intexp + intexp2 = intw*intexp2 + intarea = intw*intarea + +! Spectral components for flow. + + if (flow) then + +! Axisymmetric. + + hflow0 = intw*hflow0 + cflow0 = intw*cflow0 + nflow0 = intw*nflow0 + +! Non-axisymmetric. + + if (nonaxi) then + + hflowc = intw*hflowc + cflowc = intw*cflowc + nflowc = intw*nflowc + + if (.not.refy) then + hflows = intw*hflows + cflows = intw*cflows + nflows = intw*nflows + end if + + end if + + end if + + end if + + +! ***************************** +! *** DEALLOCATE MEMORY *** +! ***************************** + + 10 continue + + deallocate(rr,xa,ya,za) + deallocate(da,exp,gradn) + deallocate(txx,tyy,tzz,txy,txz,tyz) + + +! *************** +! *** END *** +! *************** + + end subroutine AHFinder_int + + + diff --git a/src/AHFinder_leg.F b/src/AHFinder_leg.F new file mode 100644 index 0000000..2297bb7 --- /dev/null +++ b/src/AHFinder_leg.F @@ -0,0 +1,305 @@ +c/*@@ +c @file AHFinder_leg.F +c @date April 1998 +c @author Miguel Alcubierre +c @desc +c Associated Legendre functions. +c @enddesc +c@@*/ +#include "cctk.h" + + + + CCTK_REAL function LEGEN(L,M,X) + +! This function computes the associated Legendre polynomial +! P^M_L(X). Here M and L are integers satisfying 0 <= M <= L, +! and X lies in the range -1 <= X <= 1. +! +! The general function is taken from Numerical Recipes. +! However, I hardwired the first few polynomials for M=0, +! and I added a section that tries to use as much information +! as possible from previous calls to the routine to avoid +! having to start from scratch every time. + + use AHFinder_dat + + implicit none + + integer L,M + integer LOLD1,LOLD2,MOLD1,MOLD2 + integer I,LL + + CCTK_REAL X,Y + CCTK_REAL PMM,PLL,PMMP1,SOMX2,FACT + CCTK_REAL factor,aux + CCTK_REAL zero,one,two + CCTK_REAL XOLD1,XOLD2 + CCTK_REAL LEGOLD1,LEGOLD2 + + save LOLD1,LOLD2,MOLD1,MOLD2,XOLD1,XOLD2 + save LEGOLD1,LEGOLD2 + + +! ************************** +! *** DEFINE NUMBERS *** +! ************************** + + zero = 0.0D0 + one = 1.0D0 + two = 2.0D0 + + +! ************************************************* +! *** HARDWIRE FIRST FEW POLYNOMIALS FOR M=0 *** +! ************************************************* + +! Notice that it doesn't make sense to go much higher +! than L=12 since the explicit expressions for the +! polynomials have large terms with alternating signs, +! so for high L we rapidly loose accuracy in the sum. +! The method from Numerical Recipes that comes below +! is slow, but doesn't have this problem. + + Y = X**2 + +! Find polynomials for M=0. + + if ((M.eq.0).and.(L.le.12)) then + + if (L.eq.0) then + LEGEN = one + goto 10 + else if (L.eq.1) then + LEGEN = X + goto 10 + else if (L.eq.2) then + LEGEN = 5.0D-1*(3.0D0*Y - one) + goto 10 + else if (L.eq.3) then + LEGEN = 5.0D-1*(5.0D0*Y - 3.0D0)*X + goto 10 + else if (L.eq.4) then + LEGEN = 1.25D-1*(35.0D0*Y**2 - 30.0D0*Y + 3.0D0) + goto 10 + else if (L.eq.5) then + LEGEN = 1.25D-1*(63.0D0*Y**2 - 70.0D0*Y + 15.0D0)*X + goto 10 + else if (L.eq.6) then + LEGEN = 6.25D-2*(231.0D0*Y**3 - 315.0D0*Y**2 + 105.0D0*Y + . - 5.0D0) + goto 10 + else if (L.eq.7) then + LEGEN = 6.25D-2*(429.0D0*Y**3 - 693.0D0*Y**2 + 315.0D0*Y + . - 35.0D0)*X + goto 10 + else if (L.eq.8) then + LEGEN = 7.8125D-3*(6435.0D0*Y**4 - 12012.0D0*Y**3 + . + 6930.0D0*Y**2 - 1260.0D0*Y + 35.0D0) + goto 10 + else if (L.eq.9) then + LEGEN = 7.8125D-3*(12155.0D0*Y**4 - 25740.0D0*Y**3 + . + 18018.0D0*Y**2 - 4620.0D0*Y + 315.0D0)*X + goto 10 + else if (L.eq.10) then + LEGEN = 3.90625D-3*(46189.0D0*Y**5 - 109395.0D0*Y**4 + . + 90090.0D0*Y**3 - 30030.0D0*Y**2 + 3465.0D0*Y + . - 63.0D0) + goto 10 + else if (L.eq.11) then + LEGEN = 3.90625D-3*(88179.0D0*Y**5 - 230945.0D0*Y**4 + . + 218790.0D0*Y**3 - 90090.0D0*Y**2 + 15015.0D0*Y + . - 693.0D0)*X + goto 10 + else if (L.eq.12) then + LEGEN = 9.76563D-4*(676039.0D0*Y**6 - 1939938.0D0*Y**5 + . + 2078505.0D0*Y**4 - 1021020.0D0*Y**3 + . + 225225.0D0*Y**2 - 18018.0D0*Y + 231.0D0) + goto 10 + end if + + end if + + +! *********************************************** +! *** USE RECURSIVE RELATIONS IF POSSIBLE *** +! *********************************************** + +! Here I try to optimize the calculation of the Legendre +! polynomials (which is ordinarily very expensive) by using +! the recurrence relations and information from the last +! calls to the routine. + + if (firstleg) then + + LOLD1 = -2 + MOLD1 = -2 + XOLD1 = zero + LEGOLD1 = zero + + firstleg = .false. + + else + +! See if we have the same values of M and X as the last call. + + if ((M.eq.MOLD1).and.(X.eq.XOLD1).and.(L.ne.M)) then + +! If L is the same as in the last call, then +! LEGEN should also be the same. This will +! almost certainly never happen. + + if (L.eq.LOLD1) then + + LEGEN = LEGOLD1 + goto 10 + +! This is the most interesting case L = LOLD1+1. + + else if (L.eq.LOLD1+1) then + +! LOLD1 = MOLD1. + + if (LOLD1.eq.MOLD1) then + + LEGEN = X*dble(2*M+1)*LEGOLD1 + goto 10 + +! LOLD1 = LOLD2 + 1. + + else if ((L.eq.LOLD2+2).and.(X.eq.XOLD2)) then + + LEGEN = (X*dble(2*L-1)*LEGOLD1 - dble(L+M-1)*LEGOLD2) + . / dble(L-M) + goto 10 + + end if + + end if + + end if + + end if + + +! ************************************* +! *** NUMERICAL RECIPES ROUTINE *** +! ************************************* + +! Compute P^M_M. + + PMM = one + + if (M.gt.0) then + + SOMX2 = dsqrt((one-X)*(one+X)) + FACT = one + + do I=1,M + PMM = - PMM*FACT*SOMX2 + FACT = FACT + two + end do + + end if + + if (L.eq.M) then + + LEGEN = PMM + +! Compute P^M_M+1. + + else + + PMMP1 = X*dble(2*M+1)*PMM + + if (L.eq.M+1) then + + LEGEN = PMMP1 + +! Compute P^M_L, L > M+1. + + else + + do LL=M+2,L + PLL = (X*dble(2*LL-1)*PMMP1 - dble(LL+M-1)*PMM) + . / dble(LL-M) + PMM = PMMP1 + PMMP1 = PLL + end do + + LEGEN = PLL + + end if + + end if + + +! *************************** +! *** SAVE OLD VALUES *** +! *************************** + + 10 continue + + LOLD2 = LOLD1 + MOLD2 = MOLD1 + XOLD2 = XOLD1 + LEGOLD2 = LEGOLD1 + + LOLD1 = L + MOLD1 = M + XOLD1 = X + LEGOLD1 = LEGEN + + +! ********************************************** +! *** MULTIPLY WITH NORMALIZATION FACTOR *** +! ********************************************** + +! For M=0 I use the normalization factor: +! +! sqrt(2*L + 1) +! +! For M non-zero I use the normalization factor: +! +! sqrt( 2 (2*L + 1) (L-M)! / (L+M)! ) +! +! The extra factor of sqrt(2) is there because I +! use a basis of sines and cosines and not the +! standard complex exponential. +! +! With this normalization, I ensure that my basis +! functions f _lm are such that: +! +! / +! | f_lm f_l'm' sin(theta) dtheta dphi = 4 pi delta_mm' delta_ll' +! / +! +! Notice the extra factor of 4 pi. This is there because +! for a sphere, I want the (0,0) coefficient to correspond +! to the radius. + + factor = dble(2*L+1) + + if (M.ne.0) then + + aux = one + + do i=L-M+1,L+M + aux = aux*dble(i) + end do + + factor = two*factor/aux + + end if + + LEGEN = dsqrt(factor)*LEGEN + + +! *************** +! *** END *** +! *************** + + end function LEGEN + + + diff --git a/src/AHFinder_mask.F b/src/AHFinder_mask.F new file mode 100644 index 0000000..9fec702 --- /dev/null +++ b/src/AHFinder_mask.F @@ -0,0 +1,81 @@ +c/*@@ +c @file AHFinder_mask.F +c @date November 1998 +c @author Miguel Alcubierre +c @desc +c Set up horizon mask for black hole excision. +c @enddesc +c@@*/ + +c Note that including cactus.h will also include AHFinder.h +!#include "cactus.h" +#include "cctk.h" +#include "cctk_parameters.h" +#include "cctk_arguments.h" + + subroutine AHFinder_mask(CCTK_FARGUMENTS) + +! Set up the horizon mask. All this routine does is +! set up the value of the grid function `ahmask' to +! one for points outside the horizon, and to zero +! for points inside the horizon. This is only necessary +! if we want to use black hole excision. + + use AHFinder_dat + + implicit none + + DECLARE_CCTK_FARGUMENTS + DECLARE_CCTK_PARAMETERS + + integer i,j,k + + CCTK_REAL zero,one + CCTK_REAL shrink + + +! ***************************** +! *** DEFINE {zero,one} *** +! ***************************** + + zero = 0.0D0 + one = 1.0D0 + + +! ************************* +! *** START ROUTINE *** +! ************************* + + call AHFinder_fun(CCTK_FARGUMENTS) + +! Find out how big is the buffer zone. Notice that +! a negative value for shrink puts the masked region +! safely inside the horizon, while a positive value +! moves the masked region outside the horizon. + +! shrink = getr8("ahf_maskshrink")*dx + shrink = ahf_maskshrink*dx + +! Mask points. + + do k=1,nz + do j=1,ny + do i=1,nx + + if (ahfgrid(i,j,k).ge.shrink) then + ahmask(i,j,k) = one + else + ahmask(i,j,k) = zero + end if + + end do + end do + end do + + +! *************** +! *** END *** +! *************** + + end subroutine AHFinder_mask + diff --git a/src/AHFinder_min.F b/src/AHFinder_min.F new file mode 100644 index 0000000..0a5e55e --- /dev/null +++ b/src/AHFinder_min.F @@ -0,0 +1,281 @@ +c/*@@ +c @file AHFinder_min.F +c @date April 1998 +c @author Miguel Alcubierre +c @desc +c Minimization routine. +c @enddesc +c@@*/ + +c Note that including cactus.h will also include AHFinder.h +!#include "cactus.h" +#include "cctk.h" +#include "cctk_parameters.h" +#include "cctk_arguments.h" + + subroutine AHFinder_min(CCTK_FARGUMENTS,NN,status,logf) + + use AHFinder_dat + + implicit none + + DECLARE_CCTK_FARGUMENTS + DECLARE_CCTK_PARAMETERS + + logical status,found + + integer i,NN + integer l,m + CCTK_INT ITER,ITMAX + + CCTK_REAL FTOL,FRET + CCTK_REAL zero,half,one + + CCTK_REAL, dimension (NN) :: P + CCTK_REAL, dimension (NN,NN) :: XI + + character*200 logf + +! Description of variables: +! +! i,j,l,m Counters. +! +! ITER Number of iterations taken in minimization. +! ITMAX Maximum number of iterations allowed. +! +! FTOL Fractional tolerance for return criterion. +! +! FRET Value of fucntion on return from minimization. +! +! P Vector of coefficients. +! +! XI Matrix of initial directions. + + +! *********************************** +! *** DEFINE {zero,half,one} *** +! *********************************** + + zero = 0.0D0 + half = 0.5D0 + one = 1.0D0 + + +! ********************************* +! *** FIND INITIAL VECTOR P *** +! ********************************* + P = zero + i = 0 + +! Copy {xc,yc,zc}. + + if (wander) then + if (.not.refx) then + i = i+1 + P(1) = xc + end if + if (.not.refy) then + i = i+1 + P(2) = yc + end if + if (.not.refz) then + i = i+1 + P(3) = zc + end if + end if + +! Copy c0(l). + + do l=0,lmax,1+stepz + i = i+1 + P(i) = c0(l) + end do + +! Copy {cc,cs}. + + if (nonaxi) then + +! Copy cc(l,m). + + do l=1,lmax + do m=1+stepx,l,1+stepx + if (stepz*mod(l-m,2).eq.0) then + i = i+1 + P(i) = cc(l,m) + end if + end do + end do + +! Copy cs(l,m). + + if (.not.refy) then + do l=1,lmax + do m=1,l,1+stepx + if (stepz*mod(l-m,2).eq.0) then + i = i+1 + P(i) = cs(l,m) + end if + end do + end do + end if + + end if + + +! ********************************************* +! *** FIND MATRIX OF INITIAL DIRECTIONS *** +! ********************************************* + +! The initial set of directions are taken to be +! the unit vectors (the first index identifies the +! components, and second index the vectors). + + XI = zero + + do i=1,NN + XI(i,i) = one + end do + + +! **************************************** +! *** CALL MINIMIZATION SUBROUTINE *** +! **************************************** + +! Find ITMAX. + +! ITMAX = geti("ahf_maxiter") + ITMAX = ahf_maxiter + +! Find FTOL. FTOL is the fractional tolerance in the +! function value so that failure to decrease by more +! than this amount on one interation signals that we +! are done. + +! FTOL = getr8("ahf_tol") + FTOL = ahf_tol + +! Here I call the routine POWELL from Numerical Recipies. + + found = .false. + firstfun = .true. + + call POWELL(CCTK_FARGUMENTS,P,XI,NN,FTOL,ITER,ITMAX,FRET,found) + + +! **************************** +! *** FIND {c0,cc,cs} *** +! **************************** + + i = 0 + +! Find {xc,yc,zc}. + + if (wander) then + if (.not.refx) then + i = i+1 + xc = P(1) + end if + if (.not.refy) then + i = i+1 + yc = P(2) + end if + if (.not.refz) then + i = i+1 + zc = P(3) + end if + end if + +! Find c0(l). + + do l=0,lmax,1+stepz + i = i+1 + c0(l) = P(i) + end do + +! Find {cc,cs}. + + if (nonaxi) then + +! Find cc(l,m). + + do l=1,lmax + do m=1+stepx,l,1+stepx + if (stepz*mod(l-m,2).eq.0) then + i = i+1 + cc(l,m) = P(i) + end if + end do + end do + +! Find cs(l,m). + + if (.not.refy) then + do l=1,lmax + do m=1,l,1+stepx + if (stepz*mod(l-m,2).eq.0) then + i = i+1 + cs(l,m) = P(i) + end if + end do + end do + end if + + end if + +! if (i.gt.NN) then +! write(*,*) 'Out of bounds 2',i,NN +! STOP +! end if + + +! *************************************** +! *** DECIDE ON THE RETURN STATUS *** +! *************************************** + +! If a minimum was found, set status to true. + + if (found) then + + status = .true. + +! If no minimum was found, it might still be because +! there were to many iterations, in which case there +! is a minimum, we just didn't converge to it in the +! maximum number of iterations allowed. + + else + + if (ITER.eq.ITMAX) then + + status = .true. + + if ((verbose.or.logfile).and.(myproc.eq.0)) then + + write(*,*) + write(*,*) 'AHFinder: Too many iterations.' + + if (logfile) then + open(1,file=logf,form='formatted', + . status='old',position='append') + write(1,*) + write(1,*) 'AHFinder: Too many iterations.' + close(1) + end if + + end if + + else + + status = .false. + + end if + + end if + + +! *************** +! *** END *** +! *************** + + end subroutine AHFinder_min + diff --git a/src/AHFinder_pow.F b/src/AHFinder_pow.F new file mode 100644 index 0000000..f605039 --- /dev/null +++ b/src/AHFinder_pow.F @@ -0,0 +1,1070 @@ +c/*@@ +c @file AHFinder_pow.F +c @date April 1998 +c @author Miguel Alcubierre +c @desc +c Multidimensional minimization routines taken +c from Numerical Recipes. The routines have some +c modifications to adapt them to the specific problem +c we are trying to solve. +c @enddesc +c@@*/ + +c Note that including cactus.h will also include AHFinder.h +!#include "cactus.h" +#include "cctk.h" +#include "cctk_parameters.h" +#include "cctk_arguments.h" + + + module F1COM + +! Data to be shared by minimization routines. + + implicit none + + CCTK_INT NCOM + + CCTK_REAL, allocatable, dimension(:) :: PCOM,XICOM + + +! *************** +! *** END *** +! *************** + + end module F1COM + + + + +*********************************************************************** +*********************************************************************** + + + + + subroutine POWELL(CCTK_FARGUMENTS,P,XI,N,FTOL,ITER,ITMAX,FRET,found) + +! Minimization of a function FUNC on N variables (FUNC is +! not an argument, it is a fixed function name). Input consists +! of an initial starting point P that is a vector of length N; +! an initial matrix XI of dimensions NxN and whose columns contain +! the initial set of directions (usually the N unit vectors); and +! FTOL, the fractional tolerance in the function value such that +! failure to decrease by more than this ammount on one iteration +! signals doneness. On output, P is set to the best point found, +! XI is the then-current direction set, FRET is the returned +! function value at P, and ITER is the number of iterations taken. +! +! Extra output added by myself (Miguel Alcubierre) is a logical +! variable `found' that is set to true is a minimum was indeed +! found, and to false if it wasn't. + + use AHFinder_dat + + implicit none + + DECLARE_CCTK_FARGUMENTS + + logical found + + CCTK_INT i,j,N + CCTK_INT ITER,ITMAX,IBIG + + CCTK_REAL FUNC,DEL,FP,FPTT,FRET,FTOL,T,TOL1,ZEPS + + CCTK_REAL, dimension(N) :: P,PT,PTT,XIT + CCTK_REAL, dimension(N,N) :: XI + + character*200 logf + + parameter (ZEPS=1.0D-10) + + +! ************************* +! *** START ROUTINE *** +! ************************* + + FRET = FUNC(CCTK_FARGUMENTS,P,N) + +! Save initial point. + + do j=1,N + PT(j) = P(j) + end do + + ITER = 0 + + 10 ITER = ITER + 1 + + if ((myproc.eq.0).and.veryver) then + write(*,*) + write(*,"(A20,I2)") 'POWELL ITERATION = ',ITER + end if + + if ((myproc.eq.0).and.logfile) then + logf = filestr(1:nfile)//"ahf_logfile" + open(1,file=logf,form='formatted',status='old', + . position='append') + write(1,*) + write(1,"(A20,I2)") 'POWELL ITERATION = ',ITER + close(1) + end if + + FP = FRET + IBIG = 0 + +! Initialize DEL. It will be the biggest function decrease. + + DEL = 0.0D0 + +! In each iteration, loop over all directions in the set. + + do i=1,N + +! Copy the direction. + + do j=1,N + XIT(j) = XI(j,i) + end do + + FPTT = FRET + +! Minimize along it. + + call LINMIN(CCTK_FARGUMENTS,P,XIT,N,FRET,FTOL) + +! And record it if it is the largest decrease so far. + + if (dabs(FPTT-FRET).gt.DEL) then + DEL = dabs(FPTT-FRET) + IBIG = I + end if + + 20 continue + + end do + +! Termination criterion. + + TOL1 = 0.5D0*FTOL*(dabs(FP)+dabs(FRET)) + ZEPS + + if (dabs(FP-FRET).le.TOL1) then + found = .true. + return + else if (dabs(FRET).lt.ZEPS) then + found = .true. + return + end if + +! Too many iterations? + + if (ITER.ge.ITMAX) return + +! Construct the extrapolated point and the average direction. + + do j=1,N + PTT(j) = 2.0D0*P(j) - PT(j) + XIT(j) = P(j) - PT(j) + PT(j) = P(j) + end do + +! Function value at extrapolated point. + + FPTT = FUNC(CCTK_FARGUMENTS,PTT,N) + +! First reason not to use new direction. + + if (FPTT.ge.FP) goto 10 + +! Second reason not to use new direction. + + T = 2.0D0*(FP-2.0D0*FRET+FPTT)*(FP-FRET-DEL)**2 + . - DEL*(FP-FPTT)**2 + + if (T.gt.0.0D0) goto 10 + +! Move to the minimum of the new direction. + + call LINMIN(CCTK_FARGUMENTS,P,XIT,N,FRET,FTOL) + +! Save new direction. + + do j=1,N + XI(j,IBIG) = XIT(j) + end do + +! Back for another iteration. + + goto 10 + + +! *************** +! *** END *** +! *************** + + end subroutine POWELL + + + + +*********************************************************************** +*********************************************************************** + + + + subroutine LINMIN(CCTK_FARGUMENTS,P,XI,N,FRET,TOL) + +! Given an N dimensional point P and an N dimensional direction XI, +! moves and resets P to where the function FUNC takes on a +! minimum along the direction XI from P, and replaces XI by the +! actual vector displacement that P was moved. Also returns as +! FRET the value of FUNC at the returned location P. This is +! actually all acomplished by calling the routines MNBRAK and +! PARABOLA. + + use F1COM + + implicit none + + DECLARE_CCTK_FARGUMENTS + + logical error + + CCTK_INT i,j,N + CCTK_INT l,m + + CCTK_REAL XA,XB,XX,FA,FB,FX + CCTK_REAL FRET,PARABOLA,XMIN,TOL + CCTK_REAL zero + CCTK_REAL dx,dy,dz + + CCTK_REAL, dimension(N) :: P,XI + + external F1DIM + + +! ********************************************** +! *** ALLOCATE STORAGE FOR SHARED ARRAYS *** +! ********************************************** + + allocate(PCOM(N),XICOM(N)) + + +! ************************* +! *** START ROUTINE *** +! ************************* + dx = cctk_delta_space(1) + dy = cctk_delta_space(2) + dz = cctk_delta_space(3) + + NCOM = N + + do j=1,N + PCOM(j) = P(j) + XICOM(j) = XI(j) + end do + +! Initial guess for brackets. + + XA = 0.0D0 + XX = 0.1D0*min(dx,dy,dz) + +! Bracket the minimum. + + call MNBRAK(CCTK_FARGUMENTS,XA,XX,XB,FA,FX,FB,F1DIM) + +! Check if function is constant. + + if ((FA.eq.FX).and.(FB.eq.FX)) goto 10 + +! Call 1D minimization. + + FRET = PARABOLA(CCTK_FARGUMENTS,XA,XX,XB,FA,FX,FB,F1DIM,TOL,XMIN,error) + +! Construct the vector results to return. + + do j=1,N + XI(j) = XMIN*XI(j) + P(j) = P(j) + XI(j) + end do + + +! ****************************** +! *** DEALLOCATE STORAGE *** +! ****************************** + + 10 deallocate(PCOM,XICOM) + + +! *************** +! *** END *** +! *************** + + end subroutine LINMIN + + + + +*********************************************************************** +*********************************************************************** + + + + + REAL function F1DIM(CCTK_FARGUMENTS,XX) + +! Must accompany LINMIN + + use F1COM + + implicit none + + DECLARE_CCTK_FARGUMENTS + + CCTK_INT j + + CCTK_REAL XX,FUNC + + CCTK_REAL, dimension(NCOM) :: XT + + +! ************************* +! *** START ROUTINE *** +! ************************* + + do j=1,NCOM + XT(j) = PCOM(j) + XX*XICOM(j) + end do + + F1DIM = FUNC(CCTK_FARGUMENTS,XT,NCOM) + + +! *************** +! *** END *** +! *************** + + end function F1DIM + + + + +*********************************************************************** +*********************************************************************** + + + + subroutine MNBRAK(CCTK_FARGUMENTS,XA,XB,XC,FA,FB,FC,FUNC) + +! Given a function FUNC, and given distinct initial points +! XA and XB, this routine searches in the downhill direction +! (defined by the function as evaluated at the initial points) +! and returns new points XA, XB, XC which bracket a minimum +! of the function. Also returned are the function values +! at the three points, FA, FB and FC. + + implicit none + + DECLARE_CCTK_FARGUMENTS + + CCTK_REAL XA,XB,XC + CCTK_REAL FA,FB,FC,FU,FUNC + CCTK_REAL DUM,Q,S,U,ULIM + CCTK_REAL GOLD,GLIMIT,TINY + + parameter (GOLD=1.618034D0,GLIMIT=4.0D0,TINY=1.0D-20) + + +! ************************* +! *** START ROUTINE *** +! ************************* + + FA = FUNC(CCTK_FARGUMENTS,XA) + FB = FUNC(CCTK_FARGUMENTS,XB) + +! Switch roles of A and B so that we can go downhill +! in the direction from A to B. + + if (FB.gt.FA) then + + DUM = XA + XA = XB + XB = DUM + + DUM = FB + FB = FA + FA = DUM + + end if + +! First guess for C. + + XC = XB + GOLD*(XB-XA) + + FC = FUNC(CCTK_FARGUMENTS,XC) + +! Check if the function is constant. + + if ((FA.eq.FB).and.(FB.eq.FC)) then + return + end if + +! Keep returning here until we bracket. + + do while(FB.gt.FC) + +! Compute U by parabollic extrapolation from A, B, C. +! TINY is used to prevent any possible division by zero. + + S = (XB-XA)*(FB-FC) + Q = (XB-XC)*(FB-FA) + U = XB - ((XB-XC)*Q - (XB-XA)*S) + . / (2.0D0*SIGN(MAX(DABS(Q-S),TINY),Q-S)) + +! We won't go further than this. + + ULIM = XB + GLIMIT*(XC-XB) + +! Now to test various possibilities. + + if ((XB-U)*(U-XC).gt.0.0D0) then + +! Parabolic U is between B and C, try it. + + FU = FUNC(CCTK_FARGUMENTS,U) + + if (FU.lt.FC) then + +! Got a minimum between B and C. + + XA = XB + FA = FB + XB = U + FB = FU + return + + else if (FU.gt.FB) then + +! Got a minimum between A and U. + + XC = U + FC = FU + return + + end if + +! Parabolic fit was no use. Use default magnification. + + U = XC + GOLD*(XC-XB) + + FU = FUNC(CCTK_FARGUMENTS,U) + + else if ((XC-U)*(U-ULIM).gt.0.0D0) then + +! Parabolic fit is between C and its allowed limit. + + FU = FUNC(CCTK_FARGUMENTS,U) + + if (FU.lt.FC) then + + XB = XC + XC = U + U = XC + GOLD*(XC-XB) + FB = FC + FC = FU + + FU = FUNC(CCTK_FARGUMENTS,U) + + end if + + else if ((U-ULIM)*(ULIM-XC).ge.0.0D0) then + +! Limit parabolic U to maximum allowed value. + + U = ULIM + + FU = FUNC(CCTK_FARGUMENTS,U) + + else + +! Reject parabolic U, use default magnification. + + U = XC + GOLD*(XC-XB) + FU = FUNC(CCTK_FARGUMENTS,U) + + end if + +! Eliminate oldest point and continue. + + XA = XB + XB = XC + XC = U + FA = FB + FB = FC + FC = FU + + end do + + +! *************** +! *** END *** +! *************** + + end subroutine MNBRAK + + + + +*********************************************************************** +*********************************************************************** + + + + + REAL function PARABOLA(CCTK_FARGUMENTS,XA,XB,XC,FXA,FXB,FXC,FUNC,TOL,XMIN, + . error) + +! Given a function F, and given a bracketing triplet of abscissas +! XA, XB, XC (such that XB is between XA and XC, and F(XB) is less +! than both F(XA) and F(XC)), this routine isolates the minimum +! to a fractional precision of about TOL using inverse parabolic +! interpolation. The abscissa of the minimum is returned as XMIN, +! and the minimum function value is returned as PARABOLA, the +! returned function value. + + implicit none + + DECLARE_CCTK_FARGUMENTS + + logical error + + CCTK_INT ITER,ITMAX + + CCTK_REAL XA,XB,XC,FXA,FXB,FXC,FUNC,TOL,XMIN + CCTK_REAL A,B,D,P,Q,S,U,XX,XM,TOL1 + CCTK_REAL FA,FB,FU,FX,FP,ZEPS + + parameter (ITMAX=100,ZEPS=1.0D-10) + + +! ************************* +! *** START ROUTINE *** +! ************************* + + error = .false. + +! A and B must be in ascending order, though the input +! abscissas need not be. + + if (XC.lt.XA) then + A = XC + B = XA + FA = FXC + FB = FXA + else + A = XA + B = XC + FA = FXA + FB = FXC + end if + + XX = XB + FX = FXB + FP = FA + +! Main program loop. + + do ITER=1,ITMAX + + XM = 0.5D0*(A + B) + +! Test for done here. + + TOL1 = 0.5D0*TOL*(dabs(FP)+dabs(FX)) + ZEPS + if (dabs(FP-FX).le.TOL1)goto 3 + +! Construct a trial parabolic fit. + + S = (XX-A)*(FX-FB) + Q = (XX-B)*(FX-FA) + P = (XX-B)*Q - (XX-A)*S + Q = 2.0D0*(Q-S) + +! Check if Q is zero. This can only happen if the +! three points are colinear (in which case they where +! not bracketing a minimum in the first place), or +! if two of them are in fact the same point. + + if (Q.eq.0.0D0) then + write(*,*) 'AHFinder: problem in PARABOLA' + error = .true. + return + end if + +! Find new point. + + D = - P/Q + U = XX + D + +! Check if X was already the minimum. + + if (dabs(U-XX).lt.ZEPS) goto 3 + +! This is the one function evaluation per iteration. + + FU = FUNC(CCTK_FARGUMENTS,U) + +! And now we have to decide what to do with our function +! evaluation. + + if (FU.le.FX) then + + if (U.ge.XX) then + A = XX + FA = FX + else + B = XX + FB = FX + end if + + XX = U + + FP = FX + FX = FU + + else + + if (U.lt.XX) then + A = U + FA = FU + else + B = U + FB = FU + end if + + end if + + end do + + write(*,*) 'PARABOLA exceeded maximum number of iterations.' + + 3 XMIN = XX + PARABOLA = FX + + +! *************** +! *** END *** +! *************** + + end function PARABOLA + + + + +*********************************************************************** +*********************************************************************** + + + + REAL function FUNC(CCTK_FARGUMENTS,P,N) + +! This is a shell to find the surface integral of the expansion. + + use AHFinder_dat + + implicit none + + DECLARE_CCTK_FARGUMENTS + DECLARE_CCTK_PARAMETERS + + logical flag + + CCTK_INT i,l,m,N + CCTK_INT iter + CCTK_INT getoutpfx + + CCTK_REAL zero,half,one + CCTK_REAL aux + + CCTK_REAL, dimension(N) :: P + + character*200 logf + + save iter,logf + + +! ************************** +! *** DEFINE NUMBERS *** +! ************************** + + zero = 0.0D0 + half = 0.5D0 + one = 1.0D0 + + +! ************************ +! *** OPEN LOGFILE *** +! ************************ + + logf = filestr(1:nfile)//"ahf_logfile" + + if (logfile.and.(myproc.eq.0)) then + open(1,file=logf,form='formatted',status='old', + . position='append') + end if + + +! **************************** +! *** FIND {c0,cc,cs} *** +! **************************** + + i = 0 + +! Find {xc,yc,zc}. + + if (wander) then + if (.not.refx) then + i = i+1 + xc = P(1) + end if + if (.not.refy) then + i = i+1 + yc = P(2) + end if + if (.not.refz) then + i = i+1 + zc = P(3) + end if + end if + +! Find c0(l). + + do l=0,lmax,1+stepz + i = i+1 + c0(l) = P(i) + end do + +! Find {cc,cs}. + + if (nonaxi) then + +! Find cc(l,m). + + do l=1,lmax + do m=1+stepx,l,1+stepx + if (stepz*mod(l-m,2).eq.0) then + i = i+1 + cc(l,m) = P(i) + end if + end do + end do + +! Find cs(l,m). + + if (.not.refy) then + do l=1,lmax + do m=1,l,1+stepx + if (stepz*mod(l-m,2).eq.0) then + i = i+1 + cs(l,m) = P(i) + end if + end do + end do + end if + + end if + + +! ********************************************* +! *** CALCULATE EXPANSION IF WE NEED TO *** +! ********************************************* + +! Firstcall? + + if (firstfun) then + + iter = 1 + +! Find horizon function and expansion. + + call AHFinder_fun(CCTK_FARGUMENTS) + call AHFinder_exp(CCTK_FARGUMENTS) + +! Save old values of coefficients. + + xc_old = xc + yc_old = yc + zc_old = zc + + c0_old = c0 + cc_old = cc + cs_old = cs + + firstfun = .false. + +! Not first call. + + else + + iter = iter + 1 + +! Check if coefficients changed from last call. +! I check this because the minimization algorithm +! sometimes reevaluates the function at the same +! point twice. Also, if only c0(0) changed, the +! expansion will not be affected. + + flag = .false. + + if (wander) then + if (xc.ne.xc_old) flag = .true. + if (yc.ne.yc_old) flag = .true. + if (zc.ne.zc_old) flag = .true. + end if + + do l=1+stepz,lmax,1+stepz + if (c0(l).ne.c0_old(l)) flag = .true. + end do + + if (nonaxi) then + do l=1,lmax + do m=1+stepx,l,1+stepx + if (stepz*mod(l-m,2).eq.0) then + if (cc(l,m).ne.cc_old(l,m)) flag = .true. + end if + end do + if (.not.refy) then + do m=1,l,1+stepx + if (stepz*mod(l-m,2).eq.0) then + if (cs(l,m).ne.cs_old(l,m)) flag = .true. + end if + end do + end if + end do + end if + +! If necessary, find horizon function and expansion. + + if (flag) then + call AHFinder_fun(CCTK_FARGUMENTS) + call AHFinder_exp(CCTK_FARGUMENTS) + end if + +! Save old values of coefficients. + + c0_old = c0 + cc_old = cc + cs_old = cs + + end if + + if (myproc.eq.0) then + + if (veryver) then + write(*,*) + write(*,"(A13,I5)") 'FUNC CALL =',iter + end if + + if (logfile) then + write(1,*) + write(1,"(A13,I5)") 'FUNC CALL =',iter + end if + + end if + + +! ********************************** +! *** FIND SURFACE INTEGRALS *** +! ********************************** + + call AHFinder_int(CCTK_FARGUMENTS) + + if (myproc.eq.0) then + + if (intarea.ne.zero) then + aux = one/intarea + else + aux = one + end if + + if (veryver) then + write(*,*) + write(*,"(A21,ES14.6)") ' Surface area =',intarea + write(*,"(A21,ES14.6)") ' Mean value of H =',aux*intexp + write(*,"(A21,ES14.6)") ' Mean value of H^2 =',aux*intexp2 + write(*,*) 'number of interpolated points: ', + . inside_min_count + write(*,*) ' number of those that are negative: ', + . inside_min_neg_count + end if + + if (logfile) then + write(1,*) + write(1,"(A21,ES14.6)") ' Surface area =',intarea + write(1,"(A21,ES14.6)") ' Mean value of H =',aux*intexp + write(1,"(A21,ES14.6)") ' Mean value of H^2 =',aux*intexp2 + write(*,*) 'number of interpolated points: ', + . inside_min_count + write(*,*) ' number of those that are negative: ', + . inside_min_neg_count + end if + + end if + + +! ****************************** +! *** WRITE COEFFICIENTS *** +! ****************************** + + if (myproc.eq.0) then + + if (veryver) then + + write(*,*) + write(*,"(A20)") ' Shape coefficients:' + + if (offset) then + write(*,*) + write(*,"(A8,ES14.6)") ' xc =',xc + write(*,"(A8,ES14.6)") ' yc =',yc + write(*,"(A8,ES14.6)") ' zc =',zc + end if + + write(*,*) + write(*,*) ' l c0_l' + write(*,*) + write(*,"(I4,A6,ES14.6)") 0,' ',c0(0) + + end if + + if (logfile) then + + write(1,*) + write(1,"(A20)") ' Shape coefficients:' + + if (offset) then + write(1,*) + write(1,"(A8,ES14.6)") ' xc =',xc + write(1,"(A8,ES14.6)") ' yc =',yc + write(1,"(A8,ES14.6)") ' zc =',zc + end if + + write(1,*) + write(1,*) ' l c0_l' + write(1,*) + write(1,"(I4,A6,ES14.6)") 0,' ',c0(0) + + end if + + do l=1+stepz,lmax,1+stepz + if (veryver) then + write(*,"(I4,A6,ES14.6)") l,' ',c0(l) + end if + if (logfile) then + write(1,"(I4,A6,ES14.6)") l,' ',c0(l) + end if + end do + + if (nonaxi) then + + if (.not.refy) then + + if (veryver) then + write(*,*) + write(*,*) ' l m cc_lm cs_lm' + write(*,*) + end if + + if (logfile) then + write(1,*) + write(1,*) ' l m cc_lm cs_lm' + write(1,*) + end if + + do l=1,lmax + do m=1,l + if (stepz*mod(l-m,2).eq.0) then + if (veryver) then + write(*,10) l,m,' ',cc(l,m),' ',cs(l,m) + end if + if (logfile) then + write(1,10) l,m,' ',cc(l,m),' ',cs(l,m) + end if + 10 format(I4,I4,A2,ES14.6,A1,ES14.6) + end if + end do + end do + + else + + if (veryver) then + write(*,*) + write(*,*) ' l m cc_lm' + write(*,*) + end if + + if (logfile) then + write(1,*) + write(1,*) ' l m cc_lm' + write(1,*) + end if + + do l=1,lmax + do m=1+stepx,l,1+stepx + if (stepz*mod(l-m,2).eq.0) then + if (veryver) then + write(*,20) l,m,' ',cc(l,m) + end if + if (logfile) then + write(1,20) l,m,' ',cc(l,m) + end if + 20 format(I4,I4,A2,ES14.6) + end if + end do + end do + + end if + end if + end if + + +! ********************* +! *** FIND FUNC *** +! ********************* + +! The value of FUNC will be equal to the surface integral +! of the square of the expansion or to the surface area, +! depending on what we want to minimize. + + if (minarea) then + FUNC = intarea + else + if(find_trapped_surface) then + FUNC = intexpdel2 + else + FUNC = intexp2 + endif + end if + + +! ************************* +! *** CLOSE LOGFILE *** +! ************************* + + if (logfile.and.(myproc.eq.0)) then + close(1) + end if + + +! *************** +! *** END *** +! *************** + + end function FUNC + + + + + + + diff --git a/src/make.code.defn b/src/make.code.defn new file mode 100644 index 0000000..0af2cc5 --- /dev/null +++ b/src/make.code.defn @@ -0,0 +1,9 @@ +# Main make.code.defn file for thorn AHFinder +# Source files in this directory + +SRCS += AHFinder_dat.F AHFinder.F \ + AHFinder_fun.F AHFinder_exp.F AHFinder_int.F \ + AHFinder_min.F AHFinder_pow.F AHFinder_leg.F \ + AHFinder_flow.F AHFinder_mask.F \ + AHFinder_find3.F AHFinder_calcsigma.F \ + AHFinder_initguess.F diff --git a/src/make.code.deps b/src/make.code.deps new file mode 100644 index 0000000..cad9578 --- /dev/null +++ b/src/make.code.deps @@ -0,0 +1,14 @@ + +# -*-Makefile-*- +# Dependencies file for thorn AHFinder. + +$(SYS_OBJD)/AHFinder.o: $(SYS_OBJD)/AHFinder_dat.o +$(SYS_OBJD)/AHFinder_fun.o: $(SYS_OBJD)/AHFinder_dat.o +#$(SYS_OBJD)/AHFinder_gau.o: $(SYS_OBJD)/AHFinder_dat.o +$(SYS_OBJD)/AHFinder_int.o: $(SYS_OBJD)/AHFinder_dat.o +$(SYS_OBJD)/AHFinder_min.o: $(SYS_OBJD)/AHFinder_dat.o +$(SYS_OBJD)/AHFinder_pow.o: $(SYS_OBJD)/AHFinder_dat.o +$(SYS_OBJD)/AHFinder_flow.o: $(SYS_OBJD)/AHFinder_dat.o +$(SYS_OBJD)/AHFinder_find3.o: $(SYS_OBJD)/AHFinder_dat.o +$(SYS_OBJD)/AHFinder_calcsigma.o: $(SYS_OBJD)/AHFinder_dat.o +$(SYS_OBJD)/AHFinder_initguess.o: $(SYS_OBJD)/AHFinder_dat.o
\ No newline at end of file |