Fixes related to ISRF changes

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

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
 
@@ -389,15 +406,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
@@ -421,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")')
 

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) = &

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