From b2f2c352a2ecf89b32ad7a90b2f5441830cae509 Mon Sep 17 00:00:00 2001 From: Thomas Robitaille Date: Fri, 24 Jul 2026 06:42:15 +0000 Subject: [PATCH 1/3] Return invalid_cell from the AMR find_cell routines when a position lies outside all level 1 grids or maps to a ghost index with no target grid, instead of indexing grids(-1) or building a cell from out-of-range indices --- src/grid/grid_geometry_amr.f90 | 12 ++++++++++++ 1 file changed, 12 insertions(+) diff --git a/src/grid/grid_geometry_amr.f90 b/src/grid/grid_geometry_amr.f90 index 30eaacde..eafde7b9 100644 --- a/src/grid/grid_geometry_amr.f90 +++ b/src/grid/grid_geometry_amr.f90 @@ -531,6 +531,13 @@ recursive type(grid_cell) function find_position_in_grid(r, ilevel, igrid) resul ilevel_new = grid%goto_level(i1, i2, i3) igrid_new = grid%goto_grid(i1, i2, i3) if(ilevel_new == 0) then + if(i1 < 1 .or. i1 > grid%n1 .or. & + & i2 < 1 .or. i2 > grid%n2 .or. & + & i3 < 1 .or. i3 > grid%n3) then + call warn("find_position_in_grid","position is outside the grid") + cell = invalid_cell + return + end if cell = new_grid_cell(i1, i2, i3, ilevel, igrid, geo) else cell = find_position_in_grid(r, ilevel_new, igrid_new) @@ -556,6 +563,11 @@ type(grid_cell) function find_cell_position(r) result(icell) integer :: igrid if(debug) write(*,'(" [debug] find_cell")') igrid = locate_grid(geo%levels(1), r) ! Find grid in level 1 + if(igrid == -1) then + call warn("find_cell","photon not in any level 1 grid") + icell = invalid_cell + return + end if icell = find_position_in_grid(r, 1, igrid) ! Refine position end function find_cell_position From 46e2fa9d776e66ca1c684ae2fc31515b7aa4a0a2 Mon Sep 17 00:00:00 2001 From: Thomas Robitaille Date: Fri, 24 Jul 2026 06:48:21 +0000 Subject: [PATCH 2/3] Require zero azimuthal velocity in the spherical grid theta wall riding test so photons with matching polar angle but nonzero relative azimuth fall through to the cone intersection branch --- src/grid/grid_geometry_spherical_3d.f90 | 22 ++++++++++++++++++---- 1 file changed, 18 insertions(+), 4 deletions(-) diff --git a/src/grid/grid_geometry_spherical_3d.f90 b/src/grid/grid_geometry_spherical_3d.f90 index 1ad35c06..9c8d8a22 100644 --- a/src/grid/grid_geometry_spherical_3d.f90 +++ b/src/grid/grid_geometry_spherical_3d.f90 @@ -839,8 +839,15 @@ subroutine find_wall(p,radial,tnearest,id_min) ! Check if photon is on wall, and if so, whether it is moving along ! wall. If so, don't check for intersections and just set iext, which ! is used to specify walls that the photon is on even without an - ! intersection. - if(p%on_wall_id%w2 == -1 .and. equal_nulp(geo%wtant(p%icell%i2), sqrt(v2_xy) / p%v%z, 10)) then + ! intersection. Moving along the wall requires the direction of + ! motion to be parallel to the cone generator through the current + ! position, so in addition to the polar angle of the velocity + ! matching the cone angle, the azimuthal component of the velocity + ! has to vanish, i.e. rv_xy has to match w * v_z * tan(theta) (the + ! same azimuth-aware term as used in adjust_wall). + if(p%on_wall_id%w2 == -1 & + & .and. equal_nulp(geo%wtant(p%icell%i2), sqrt(v2_xy) / p%v%z, 10) & + & .and. equal_nulp(sqrt(r2_xy) * p%v%z * geo%wtant(p%icell%i2), rv_xy, 10)) then iext%w2 = -1 @@ -910,8 +917,15 @@ subroutine find_wall(p,radial,tnearest,id_min) ! Check if photon is on wall, and if so, whether it is moving along ! wall. If so, don't check for intersections and just set iext, which ! is used to specify walls that the photon is on even without an - ! intersection. - if(p%on_wall_id%w2 == +1 .and. equal_nulp(geo%wtant(p%icell%i2 + 1), sqrt(v2_xy) / p%v%z, 10)) then + ! intersection. Moving along the wall requires the direction of + ! motion to be parallel to the cone generator through the current + ! position, so in addition to the polar angle of the velocity + ! matching the cone angle, the azimuthal component of the velocity + ! has to vanish, i.e. rv_xy has to match w * v_z * tan(theta) (the + ! same azimuth-aware term as used in adjust_wall). + if(p%on_wall_id%w2 == +1 & + & .and. equal_nulp(geo%wtant(p%icell%i2 + 1), sqrt(v2_xy) / p%v%z, 10) & + & .and. equal_nulp(sqrt(r2_xy) * p%v%z * geo%wtant(p%icell%i2 + 1), rv_xy, 10)) then iext%w2 = +1 From b0684c798e90539c16d8ad78d0d68661309516bd Mon Sep 17 00:00:00 2001 From: Thomas Robitaille Date: Tue, 28 Jul 2026 13:37:59 +0000 Subject: [PATCH 3/3] Add regression tests for the AMR outside-grid guard and the theta wall riding condition --- hyperion/model/tests/test_grid_geometry.py | 103 +++++++++++++++++++++ 1 file changed, 103 insertions(+) create mode 100644 hyperion/model/tests/test_grid_geometry.py diff --git a/hyperion/model/tests/test_grid_geometry.py b/hyperion/model/tests/test_grid_geometry.py new file mode 100644 index 00000000..f3426c78 --- /dev/null +++ b/hyperion/model/tests/test_grid_geometry.py @@ -0,0 +1,103 @@ +from __future__ import print_function, division + +import os +import tempfile + +import numpy as np +import pytest + +from .. import Model +from ...grid import AMRGrid +from ...util.functions import random_id +from .test_helpers import get_test_dust + + +@pytest.mark.requires_hyperion_binaries +def test_amr_source_outside_grid(): + # Regression test: for AMR grids, a photon emitted at a position outside + # all level 1 grids used to index grids(-1) (undefined behavior, in + # practice a segmentation fault). find_cell now returns invalid_cell, so + # the run stops with a controlled, informative error instead. + + amr = AMRGrid() + level = amr.add_level() + grid = level.add_grid() + grid.xmin, grid.xmax = -1., 1. + grid.ymin, grid.ymax = -1., 1. + grid.zmin, grid.zmax = -1., 1. + grid.nx = grid.ny = grid.nz = 4 + grid.quantities['density'] = np.ones((4, 4, 4)) * 1.e-30 + + m = Model() + m.set_grid(amr) + m.add_density_grid(amr['density'], get_test_dust()) + + s = m.add_point_source() + s.position = (2., 0., 0.) # outside the level 1 grid + s.luminosity = 1. + s.temperature = 6000. + + m.set_n_initial_iterations(1) + m.set_n_photons(initial=100, imaging=0) + + tmpdir = tempfile.mkdtemp() + m.write(os.path.join(tmpdir, random_id())) + logfile = os.path.join(tmpdir, 'log') + + with pytest.raises(SystemExit): + m.run(os.path.join(tmpdir, 'out.rtout'), logfile=logfile) + + # The run should have failed through the controlled error path (with the + # warning from find_cell and the explanatory error from emit), not + # through a crash. + log = open(logfile).read() + assert 'photon not in any level 1 grid' in log + assert 'not emitted inside a cell' in log + + +@pytest.mark.requires_hyperion_binaries +def test_theta_wall_viewing_angle(): + # Regression test for the theta cone wall "riding" condition: peeloff + # rays whose direction has exactly the same polar angle as a theta wall + # (e.g. an observer at 60 degrees with a wall at 60 degrees) but a + # different azimuth used to be frozen on the wall when crossing it, and + # were then killed, so the SED at that exact viewing angle collapsed to + # zero. The source position and observer azimuth here are chosen so that + # every direct peeloff ray at 60 degrees crosses the 60 degree cone. + + m = Model() + r = np.array([0., 0.5, 1.0]) + t = np.array([0., np.radians(60.), np.pi]) # theta wall exactly at 60 deg + p = np.array([0., 2. * np.pi]) + m.set_spherical_polar_grid(r, t, p) + m.add_density_grid(np.ones((1, 2, 2)) * 1.e-8, get_test_dust()) + + s = m.add_point_source() + s.position = (0.3, 0., 0.05) # equator side of the cone, azimuth 0 + s.luminosity = 1. + s.temperature = 6000. + + sed = m.add_peeled_images(sed=True, image=False) + # opposing azimuth ensures the peeloff rays cross the cone + sed.set_viewing_angles([59., 60., 61.], [170., 170., 170.]) + sed.set_wavelength_range(10, 0.1, 100.) + + m.set_n_initial_iterations(0) + # check the cell consistency at every step so that mis-binned photons + # cannot go unnoticed + m.set_propagation_check_frequency(1.) + m.set_n_photons(imaging=5000) + + tmpdir = tempfile.mkdtemp() + m.write(os.path.join(tmpdir, random_id())) + mo = m.run(os.path.join(tmpdir, 'out.rtout')) + + wav, nufnu = mo.get_sed() + nufnu = np.array(nufnu) + tot = np.array([np.nansum(nufnu[i]) for i in range(3)]) + + # all three flux totals should be positive and essentially identical: + # the model is optically thin, so nothing special may happen at the + # viewing angle that exactly matches the theta wall + assert np.all(tot > 0.) + np.testing.assert_allclose(tot[1], 0.5 * (tot[0] + tot[2]), rtol=0.1)