From 8ff6316488d8fced54356be7ee3f10d8b133e0ab Mon Sep 17 00:00:00 2001 From: Thomas Robitaille Date: Tue, 21 Jul 2026 21:48:43 +0000 Subject: [PATCH 1/3] Only reset the frequency-resolved specific energy in cells that exceed the sublimation threshold when the sublimation mode is cap, rather than wiping it in every cell --- .../tests/test_specific_energy_spectrum.py | 24 +++++++++++++++++++ src/grid/grid_physics_3d.f90 | 13 +++++----- 2 files changed, 30 insertions(+), 7 deletions(-) diff --git a/hyperion/model/tests/test_specific_energy_spectrum.py b/hyperion/model/tests/test_specific_energy_spectrum.py index 78a639b4..61ccc638 100644 --- a/hyperion/model/tests/test_specific_energy_spectrum.py +++ b/hyperion/model/tests/test_specific_energy_spectrum.py @@ -162,6 +162,30 @@ def test_specific_energy_spectrum_custom_frequency_grid(tmpdir): _assert_spectrum_sums_to_specific_energy(out.filename) +@pytest.mark.requires_hyperion_binaries +def test_specific_energy_spectrum_with_sublimation_cap(tmpdir): + # With a dust type using the 'cap' sublimation mode, cells below the + # sublimation threshold must keep their frequency-resolved spectrum. This + # guards against the bug where sublimate_dust reset the spectrum in every + # cell rather than only in cells above the threshold, wiping the output. + m = Model() + m.set_cartesian_grid([-1., 0., 1.], [-1., 0., 1.], [-1., 0., 1.]) + dust = get_test_dust() + dust.set_sublimation_temperature('cap', temperature=1600.) + m.add_density_grid(np.ones((2, 2, 2)) * 1.e-16, dust) + s = m.add_point_source() + s.luminosity = 1. + s.temperature = 6000. + m.set_n_initial_iterations(3) + m.set_n_photons(initial=100000, imaging=0) + m.set_seed(-12345) + m.conf.output.output_specific_energy = 'last' + m.conf.output.output_specific_energy_spectrum = 'last' + m.write(tmpdir.join(random_id()).strpath) + out = m.run(tmpdir.join(random_id()).strpath) + _assert_spectrum_sums_to_specific_energy(out.filename) + + def test_specific_energy_spectrum_frequencies_roundtrip(tmpdir): # The custom frequency grid should survive a write/read round-trip. from ...conf.conf_files import RunConf diff --git a/src/grid/grid_physics_3d.f90 b/src/grid/grid_physics_3d.f90 index c3ce4948..3c969080 100644 --- a/src/grid/grid_physics_3d.f90 +++ b/src/grid/grid_physics_3d.f90 @@ -389,15 +389,14 @@ subroutine sublimate_dust() if(specific_energy(ic, id) > d(id)%sublimation_specific_energy) then specific_energy(ic, id) = d(id)%sublimation_specific_energy reset = reset + 1 - end if - - if (compute_specific_energy_spectrum) then - do idx=1,n_nu_bins - specific_energy_spectrum(ic,id,idx) = minimum_specific_energy(id) - end do - end if + if (compute_specific_energy_spectrum) then + do idx=1,n_nu_bins + specific_energy_spectrum(ic,id,idx) = minimum_specific_energy(id) + end do + end if + end if end do if(reset > 0) write(*,'(" [sublimate_dust] capping dust specific_energy in ",I0," cells")') reset From 34c5029769ec2fff18146ecfad991c404222f396 Mon Sep 17 00:00:00 2001 From: Thomas Robitaille Date: Tue, 21 Jul 2026 21:56:26 +0000 Subject: [PATCH 2/3] Only allocate the frequency-resolved specific energy arrays when they are requested, guarding whole-array operations with allocated checks, which also avoids referencing the first dust type in models with no dust --- src/grid/grid_physics_3d.f90 | 52 ++++++++++++++++++++++++------------ src/mpi/mpi_routines.f90 | 16 ++++++----- 2 files changed, 44 insertions(+), 24 deletions(-) diff --git a/src/grid/grid_physics_3d.f90 b/src/grid/grid_physics_3d.f90 index 3c969080..13ac5982 100644 --- a/src/grid/grid_physics_3d.f90 +++ b/src/grid/grid_physics_3d.f90 @@ -108,23 +108,35 @@ subroutine setup_grid_physics(group, use_mrw, use_pda, compute_specific_energy_s integer :: n_nu_bins ! specific_energy_spectrum is binned onto a user-specified frequency grid if one was given, - ! otherwise onto the frequency grid of the first dust type. - if (allocated(specific_energy_spectrum_frequencies)) then - n_nu_bins = size(specific_energy_spectrum_frequencies) + ! otherwise onto the frequency grid of the first dust type. The spectrum + ! arrays are only allocated when the frequency-resolved specific energy is + ! actually requested (they can be large), so all whole-array operations on + ! them are guarded by allocated() or compute_specific_energy_spectrum + ! checks. + if (compute_specific_energy_spectrum) then + if (allocated(specific_energy_spectrum_frequencies)) then + n_nu_bins = size(specific_energy_spectrum_frequencies) + else if (n_dust > 0) then + n_nu_bins = d(1)%n_nu + else + n_nu_bins = 0 + end if else - n_nu_bins = d(1)%n_nu + n_nu_bins = 0 end if ! Density allocate(density(geo%n_cells, n_dust)) allocate(specific_energy(geo%n_cells, n_dust)) - allocate(specific_energy_spectrum(geo%n_cells,n_dust,n_nu_bins)) - ! initialize: this array is only fully populated later (from the - ! photon deposits, or from the minimum specific energy), but with - ! specific_energy_type='additional' it is COPIED into - ! specific_energy_additional_spectrum before that happens -- copying - ! an uninitialized allocation injects heap garbage into the - ! frequency-resolved output at every iteration. - specific_energy_spectrum = 0._dp + if (compute_specific_energy_spectrum) then + allocate(specific_energy_spectrum(geo%n_cells,n_dust,n_nu_bins)) + ! initialize: this array is only fully populated later (from the + ! photon deposits, or from the minimum specific energy), but with + ! specific_energy_type='additional' it is COPIED into + ! specific_energy_additional_spectrum before that happens -- copying + ! an uninitialized allocation injects heap garbage into the + ! frequency-resolved output at every iteration. + specific_energy_spectrum = 0._dp + end if if(n_dust > 0) then @@ -203,7 +215,9 @@ subroutine setup_grid_physics(group, use_mrw, use_pda, compute_specific_energy_s ! energy. After the first iteration, specific_energy will get ! re-calculated and we will then add specific_energy_additional specific_energy_additional = specific_energy - specific_energy_additional_spectrum = specific_energy_spectrum + if (compute_specific_energy_spectrum) then + specific_energy_additional_spectrum = specific_energy_spectrum + end if do id=1,n_dust specific_energy(:,id) = minimum_specific_energy(id) @@ -247,8 +261,10 @@ subroutine setup_grid_physics(group, use_mrw, use_pda, compute_specific_energy_s specific_energy_sum = 0._dp ! Set up basics for the frequency-resolved specific energy - allocate(specific_energy_sum_spectrum(geo%n_cells, n_dust, n_nu_bins)) - specific_energy_sum_spectrum = 0._dp + if (compute_specific_energy_spectrum) then + allocate(specific_energy_sum_spectrum(geo%n_cells, n_dust, n_nu_bins)) + specific_energy_sum_spectrum = 0._dp + end if ! Cache the frequency grid (and its log) once so they do not have to be ! rebuilt for every photon. Photons are binned to the nearest grid point in @@ -334,7 +350,8 @@ subroutine sublimate_dust() integer :: reset integer :: n_nu_bins - n_nu_bins = size(specific_energy_spectrum, 3) + n_nu_bins = 0 + if (compute_specific_energy_spectrum) n_nu_bins = size(specific_energy_spectrum, 3) reset = 0 @@ -420,7 +437,8 @@ subroutine update_energy_abs(scale) integer :: id,idx integer :: n_nu_bins - n_nu_bins = size(specific_energy_spectrum, 3) + n_nu_bins = 0 + if (compute_specific_energy_spectrum) n_nu_bins = size(specific_energy_spectrum, 3) if(main_process()) write(*,'(" [grid_physics] updating energy_abs")') diff --git a/src/mpi/mpi_routines.f90 b/src/mpi/mpi_routines.f90 index bb585432..0069b4c4 100644 --- a/src/mpi/mpi_routines.f90 +++ b/src/mpi/mpi_routines.f90 @@ -280,13 +280,15 @@ subroutine mp_collect_physical_arrays() call mpi_reduce(specific_energy_sum, dummy_dp, size(specific_energy_sum), mpi_real8, mpi_sum, rank_main, mpi_comm_world, ierr) end if - if(main_process()) then - allocate(tmp_3d(size(specific_energy_sum_spectrum,1),size(specific_energy_sum_spectrum,2),size(specific_energy_sum_spectrum,3))) - call mpi_reduce(specific_energy_sum_spectrum, tmp_3d, size(specific_energy_sum_spectrum), mpi_real8, mpi_sum, rank_main, mpi_comm_world, ierr) - specific_energy_sum_spectrum = tmp_3d - deallocate(tmp_3d) - else - call mpi_reduce(specific_energy_sum_spectrum, dummy_dp, size(specific_energy_sum_spectrum), mpi_real8, mpi_sum, rank_main, mpi_comm_world, ierr) + if(allocated(specific_energy_sum_spectrum)) then + if(main_process()) then + allocate(tmp_3d(size(specific_energy_sum_spectrum,1),size(specific_energy_sum_spectrum,2),size(specific_energy_sum_spectrum,3))) + call mpi_reduce(specific_energy_sum_spectrum, tmp_3d, size(specific_energy_sum_spectrum), mpi_real8, mpi_sum, rank_main, mpi_comm_world, ierr) + specific_energy_sum_spectrum = tmp_3d + deallocate(tmp_3d) + else + call mpi_reduce(specific_energy_sum_spectrum, dummy_dp, size(specific_energy_sum_spectrum), mpi_real8, mpi_sum, rank_main, mpi_comm_world, ierr) + end if end if if(allocated(n_photons)) then From ec28788deec03002685646a6f26163ddea4e2f1c Mon Sep 17 00:00:00 2001 From: Thomas Robitaille Date: Tue, 21 Jul 2026 22:07:39 +0000 Subject: [PATCH 3/3] Find the specific_energy_spectrum frequency bin once at the start of grid_integrate rather than at every cell crossing since the photon frequency cannot change during the call --- src/grid/grid_propagate_3d.f90 | 10 ++++++---- 1 file changed, 6 insertions(+), 4 deletions(-) diff --git a/src/grid/grid_propagate_3d.f90 b/src/grid/grid_propagate_3d.f90 index 41549b4b..d2a3cc8c 100644 --- a/src/grid/grid_propagate_3d.f90 +++ b/src/grid/grid_propagate_3d.f90 @@ -59,6 +59,12 @@ subroutine grid_integrate(p,tau_required,tau_achieved) integer :: idx + ! The photon frequency does not change within a single grid_integrate call + ! (re-emission and scattering happen between calls, which is also when the + ! cached opacities in p%current_kappa are refreshed), so the frequency bin + ! only needs to be found once per call rather than at every cell crossing. + if (compute_specific_energy_spectrum) idx = minloc(abs(log_nu_bins - log10(p%nu)), DIM=1) + radial = (p%r .dot. p%v) > 0. if(debug) write(*,'(" [debug] start grid_integrate")') @@ -132,9 +138,6 @@ subroutine grid_integrate(p,tau_required,tau_achieved) p%r = p%r + tmin * p%v tau_achieved = tau_achieved + tau_cell - - if (compute_specific_energy_spectrum) idx = minloc(abs(log_nu_bins - log10(p%nu)), DIM=1) - do id=1,n_dust if(density(p%icell%ic, id) > 0._dp) then specific_energy_sum(p%icell%ic, id) = & @@ -202,7 +205,6 @@ subroutine grid_integrate(p,tau_required,tau_achieved) ! the spectrum low in any cell optically thick enough for ! photons to interact within a single crossing. if (compute_specific_energy_spectrum) then - idx = minloc(abs(log_nu_bins - log10(p%nu)), DIM=1) do id=1,n_dust if(density(p%icell%ic, id) > 0._dp) then specific_energy_sum_spectrum(p%icell%ic, id, idx) = &