On Sat, Apr 17, 2021 at 8:32 PM Konrad Topolski k.topolski2@student.uw.edu.pl wrote:

Changing the shiftout did not help unfortunately.

Here's the link to the output, errors, parameter file and my running script: Let me know if that's a preferable format.

https://gist.github.com/konrad-topolski/57af20b18539cd7bcb4ae2d11993144f

If it helps troubleshooting in any way, when I set up fully nonsymmetric version (i.e. zmin=-120) (no reflection symmetry in z), the error persists but now the last coordinate is affected, as in:

CCTK_InterpLocalUniform(): interpolation point is either outside the grid, or inside but too close to the grid boundary! (this may be caused by a global interpolation with driver::ghost_size too small) 0-origin interpolation point number pt=2963 of N_interp_points=2964 interpolation point (x,y,z)=(2.63796e+243,1.58133e+243,3.77779e+243) grid x_min(delta_x)x_max = -122.5(1.25)122.5 grid y_min(delta_y)y_max = -122.5(1.25)122.5 grid z_min(delta_z)z_max = -122.5(1.25)3.75

Ah, these coordinates come from the interpolator! I thought this was the grid structure description output by Carpet. That's fine, then; the local interpolator sees only part of the grid at a time.

In the output I see this:

INFO (AHFinderDirect): AH 1/3: r=0.105016 at (2.999966,-0.002037,0.000000) INFO (AHFinderDirect): AH 1/3: area=3.141592763 m_irreducible=0.2500000044

which means that the first AH is found fine; it has the right mass. The coordinate radius is 0.1, which is a pointer as to what grid resolution you need there. For starters, I'd say the resolution there needs to be 0.1/15 (or 0.1/20); you can then fine-tune this later depending on how things converge.

Information about the second horizon is not visible since it's calculated on the second MPI process. You can either use the "-roe" command-line option to Cactus to see that output, or look for "BH_diagnostics" files. If you have 2 such files, then both horizons are found, and the horizon properties are described there.

When the interpolator reports an error, it is called from the function "QuasiLocalMeasures::qlm_calculate". This thorn examines horizons after they are found, and e.g. calculates their spin. (It also calculates many other quantities, but most people are only interested in the spin.) It shouldn't be called when a horizon is not found. Maybe that's the error – a miscommunication between AHFinderDirect, SphericalSurface, and QuasiLocalMeasures. I'll check the parameter file, brb...

... am back.

There are 3 surfaces (good); one for each horizon. There are 3 horizons; each is stored into the correct surface. (Horizons are numbered starting from 1, surfaces starting from 0.). Also correct. QuasiLocalMeasures examines 3 surfaces. Also correct.

This looks wrong: INFO (SphericalSurface): SphericalSurface[0]: finest refinement-level that fully contains this surface: sf_maxreflevel = 1920103026 I don't know where this comes from. Maybe someone else can explain?

This line AHF BH_diagnostics::save[2] found_flag=0 indicates that the second horizon wasn't found. I would start debugging there. You set the AH finder to be quite verbose (good for debugging!), and output such as INFO (AHFinderDirect): proc 1/horizon 2:it 19 r_grid=0.634 ||Theta||=2.1e+00 indicates that the AH finder is currently examining a surface at r=0.634 (this looks very large; the horizon should be much smaller), and the expansion of that surface is 2.1 (that's very large; a horizon has expansion 0). You need to find a better initial guess.

I looked at the AH finding iterations: $ grep 'proc 1/horizon 2' output.out INFO (AHFinderDirect): proc 1/horizon 2:it 1 r_grid=0.125 ||Theta||=1.1e+01 INFO (AHFinderDirect): proc 1/horizon 2:it 2 r_grid=0.137 ||Theta||=9.9e+00 INFO (AHFinderDirect): proc 1/horizon 2:it 3 r_grid=0.151 ||Theta||=9.0e+00 INFO (AHFinderDirect): proc 1/horizon 2:it 4 r_grid=0.166 ||Theta||=8.2e+00 INFO (AHFinderDirect): proc 1/horizon 2:it 5 r_grid=0.182 ||Theta||=7.5e+00 INFO (AHFinderDirect): proc 1/horizon 2:it 6 r_grid=0.199 ||Theta||=6.8e+00 INFO (AHFinderDirect): proc 1/horizon 2:it 7 r_grid=0.219 ||Theta||=6.2e+00 INFO (AHFinderDirect): proc 1/horizon 2:it 8 r_grid=0.240 ||Theta||=5.7e+00 INFO (AHFinderDirect): proc 1/horizon 2:it 9 r_grid=0.263 ||Theta||=5.2e+00 INFO (AHFinderDirect): proc 1/horizon 2:it 10 r_grid=0.288 ||Theta||=4.7e+00 INFO (AHFinderDirect): proc 1/horizon 2:it 11 r_grid=0.315 ||Theta||=4.3e+00 INFO (AHFinderDirect): proc 1/horizon 2:it 12 r_grid=0.344 ||Theta||=3.9e+00 INFO (AHFinderDirect): proc 1/horizon 2:it 13 r_grid=0.377 ||Theta||=3.6e+00 INFO (AHFinderDirect): proc 1/horizon 2:it 14 r_grid=0.412 ||Theta||=3.3e+00 INFO (AHFinderDirect): proc 1/horizon 2:it 15 r_grid=0.450 ||Theta||=3.0e+00 INFO (AHFinderDirect): proc 1/horizon 2:it 16 r_grid=0.491 ||Theta||=2.8e+00 INFO (AHFinderDirect): proc 1/horizon 2:it 17 r_grid=0.535 ||Theta||=2.5e+00 INFO (AHFinderDirect): proc 1/horizon 2:it 18 r_grid=0.583 ||Theta||=2.3e+00 INFO (AHFinderDirect): proc 1/horizon 2:it 19 r_grid=0.634 ||Theta||=2.1e+00 INFO (AHFinderDirect): proc 1/horizon 2:it 20 r_grid=0.689 ||Theta||=2.0e+00

The initial guess is probably too small (maybe r=0.3 is better?), and the AH finder never latches onto a horizon surface. Mind that a surface can have an arbitrary shape, and the output shows only the average radius. If the surface becomes distorted, but the actual horizon is spherical, then the AH finder won't find it.

If you find it difficult to find the horizon initially (this can be a tedious process), then you can set up a run and look for 50 horizons, all with different initial radii, or with different initial shapes (see "coordinate ellipsoid"). Look also at https://einsteintoolkit.org/thornguide/EinsteinAnalysis/AHFinderDirect/documentation.html for the gory details.

Good luck hunting for the horizon!

-erik