From 1fd6f3b91bebf86c5cc142351395888979c23f02 Mon Sep 17 00:00:00 2001 From: Thomas Robitaille Date: Fri, 24 Jul 2026 06:46:44 +0000 Subject: [PATCH 1/4] Rescale the specific energy spectrum in PDA cells so it stays consistent with the diffusion-updated specific energy --- .../tests/test_specific_energy_spectrum.py | 20 ++++++++++++++++++- src/grid/grid_pda_3d.f90 | 14 ++++++++++++- 2 files changed, 32 insertions(+), 2 deletions(-) diff --git a/hyperion/model/tests/test_specific_energy_spectrum.py b/hyperion/model/tests/test_specific_energy_spectrum.py index 61ccc638..6cd6c07e 100644 --- a/hyperion/model/tests/test_specific_energy_spectrum.py +++ b/hyperion/model/tests/test_specific_energy_spectrum.py @@ -25,17 +25,22 @@ def _read_dataset(filename, suffix): return np.asarray(f[matches[-1]][()], dtype=float) -def _assert_spectrum_sums_to_specific_energy(filename): +def _assert_spectrum_sums_to_specific_energy(filename, exclude_empty_spectrum=False): # Summed over frequency, specific_energy_spectrum must reproduce specific_energy. # Cells with no absorption are clamped up to the minimum specific energy # (which only affects specific_energy, not the unclamped specific_energy_spectrum), # so we compare only cells that were genuinely heated above that floor. + # With exclude_empty_spectrum, cells whose spectrum is identically zero are + # also skipped: the PDA can heat cells that received no photons at all, and + # for those there is no Monte-Carlo spectral shape to rescale. se = _read_dataset(filename, '/specific_energy') se_nu = _read_dataset(filename, '/specific_energy_spectrum') assert se_nu is not None nu_sum = se_nu.sum(axis=0) floor = se.min() heated = se > floor * (1. + 1.e-6) + if exclude_empty_spectrum: + heated &= nu_sum > 0. assert np.count_nonzero(heated) >= 1 np.testing.assert_allclose(nu_sum[heated], se[heated], rtol=1.e-6) @@ -85,6 +90,19 @@ def test_specific_energy_spectrum_sums_to_specific_energy(tmpdir): _assert_spectrum_sums_to_specific_energy(out.filename) +@pytest.mark.requires_hyperion_binaries +def test_specific_energy_spectrum_sums_to_specific_energy_with_pda(tmpdir): + # The partial diffusion approximation overwrites specific_energy in cells + # with few photon deposits, so the frequency-resolved values in those cells + # have to be rescaled to match, rather than keeping the noisy Monte-Carlo + # values. Few photons over many cells guarantees PDA cells exist. + m = _cartesian_model('last', n_cells=8, n_photons=2000) + m.set_pda(True) + m.write(tmpdir.join(random_id()).strpath) + out = m.run(tmpdir.join(random_id()).strpath) + _assert_spectrum_sums_to_specific_energy(out.filename, exclude_empty_spectrum=True) + + @pytest.mark.requires_hyperion_binaries def test_specific_energy_spectrum_frequencies_written(tmpdir): # specific_energy_spectrum_frequencies is a per-frequency (1-D) quantity, not diff --git a/src/grid/grid_pda_3d.f90 b/src/grid/grid_pda_3d.f90 index dbe237b9..e23dc839 100644 --- a/src/grid/grid_pda_3d.f90 +++ b/src/grid/grid_pda_3d.f90 @@ -11,6 +11,7 @@ module grid_pda use dust_main use grid_io use grid_pda_geometry + use settings, only : compute_specific_energy_spectrum implicit none save @@ -36,10 +37,11 @@ subroutine update_specific_energy(ic) implicit none integer,intent(in) :: ic integer :: id - real(dp) :: s_prev, s, smin, smax + real(dp) :: s_prev, s, s_old, smin, smax do id=1,n_dust s = specific_energy(ic, id) + s_old = s smin = d(id)%specific_energy(1) smax = d(id)%specific_energy(d(id)%n_e) @@ -58,6 +60,16 @@ subroutine update_specific_energy(ic) end do end if specific_energy(ic, id) = s + + ! Keep the frequency-resolved specific energy consistent with the new + ! scalar value by rescaling the spectrum bins, preserving their shape. + ! If the cell had no energy before, the bins are all zero and there is + ! no shape to rescale, so they are left unchanged. + if(compute_specific_energy_spectrum) then + if(s_old > 0._dp) then + specific_energy_spectrum(ic, id, :) = specific_energy_spectrum(ic, id, :) * (s / s_old) + end if + end if end do end subroutine update_specific_energy From 9f69aa371439f92e7cb1b9826ecc9f7a38bb4ec7 Mon Sep 17 00:00:00 2001 From: Thomas Robitaille Date: Tue, 28 Jul 2026 21:29:09 +0000 Subject: [PATCH 2/4] Deposit MRW energy into the frequency-resolved specific energy spectrum by distributing it over the frequency bins according to the local emissivity --- .../tests/test_specific_energy_spectrum.py | 61 ++++++++++++++++++- src/dust/dust_type_4elem.f90 | 27 ++++++++ src/grid/grid_mrw_3d.f90 | 21 ++++++- src/grid/grid_physics_3d.f90 | 35 +++++++++++ 4 files changed, 141 insertions(+), 3 deletions(-) diff --git a/hyperion/model/tests/test_specific_energy_spectrum.py b/hyperion/model/tests/test_specific_energy_spectrum.py index 6cd6c07e..145fc545 100644 --- a/hyperion/model/tests/test_specific_energy_spectrum.py +++ b/hyperion/model/tests/test_specific_energy_spectrum.py @@ -10,7 +10,7 @@ from .. import Model from ...grid import AMRGrid from ...util.functions import random_id -from .test_helpers import get_test_dust +from .test_helpers import get_test_dust, get_realistic_test_dust def _read_dataset(filename, suffix): @@ -180,6 +180,33 @@ 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_mrw(tmpdir): + # Energy deposited during modified-random-walk steps must appear in the + # frequency-resolved spectrum as well as in the scalar specific energy. + # This guards against the bug where grid_do_mrw updated only + # specific_energy_sum, so the spectrum was biased low in the optically + # thick cells where the MRW handles most of the absorption. + m = Model() + m.set_cartesian_grid([-1., 0., 1.], [-1., 0., 1.], [-1., 0., 1.]) + # The density is chosen so that the cells are optically thick enough to + # the local Planck-mean opacity for MRW steps to actually occur. + m.add_density_grid(np.ones((2, 2, 2)) * 1.e5, get_realistic_test_dust()) + s = m.add_point_source() + s.luminosity = 1. + s.temperature = 6000. + m.set_n_initial_iterations(3) + m.set_n_photons(initial=1000, imaging=0) + m.set_seed(-12345) + m.set_mrw(True, gamma=2.) + m.set_max_interactions(1000000000) + 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) + + @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 @@ -310,3 +337,35 @@ def test_specific_energy_spectrum_mpi_matches_serial(tmpdir): total_serial = np.nansum(_read_dataset(out_serial.filename, '/specific_energy_spectrum')) total_mpi = np.nansum(_read_dataset(out_mpi.filename, '/specific_energy_spectrum')) np.testing.assert_allclose(total_mpi, total_serial, rtol=2.e-2) + + +@pytest.mark.requires_hyperion_binaries +def test_specific_energy_spectrum_sums_to_specific_energy_with_mrw(tmpdir): + # In cells heated via the Modified Random Walk, the deposited energy is + # distributed over the frequency bins according to the local emissivity + # (the radiation field is Planckian in the diffusion regime), so the sum + # over frequency must still reproduce the scalar specific energy. The + # high density here makes the model optically thick enough for the MRW + # to handle most of the energy deposition. + m = _cartesian_model('last', n_photons=2000, density=8.) + m.set_mrw(True, gamma=2.) + m.write(tmpdir.join(random_id()).strpath) + out = m.run(tmpdir.join(random_id()).strpath) + _assert_spectrum_sums_to_specific_energy(out.filename) + + +@pytest.mark.requires_hyperion_binaries +def test_specific_energy_spectrum_is_passive_with_mrw(tmpdir): + # The MRW spectrum deposit distributes energy deterministically instead + # of sampling a frequency, so enabling the spectrum must not consume + # random numbers and therefore must not change the specific energy at + # all, even in MRW-dominated models. + se = {} + for output in ('none', 'last'): + m = _cartesian_model(output_specific_energy_spectrum=output, + n_photons=2000, density=8.) + m.set_mrw(True, gamma=2.) + m.write(tmpdir.join(random_id()).strpath) + out = m.run(tmpdir.join(random_id()).strpath) + se[output] = _read_dataset(out.filename, '/specific_energy') + np.testing.assert_array_equal(se['last'], se['none']) diff --git a/src/dust/dust_type_4elem.f90 b/src/dust/dust_type_4elem.f90 index 5558f513..ff90ce89 100644 --- a/src/dust/dust_type_4elem.f90 +++ b/src/dust/dust_type_4elem.f90 @@ -20,6 +20,7 @@ module type_dust public :: get_j_nu_interp public :: get_j_nu_binned + public :: get_j_nu_bin_fractions public :: get_chi_nu_interp public :: get_chi_nu_binned @@ -748,6 +749,32 @@ function get_j_nu_binned(d, n_nu, nu_min, nu_max, jnu_var_id) result(j_nu) end function get_j_nu_binned + function get_j_nu_bin_fractions(d, jnu_var_id, nu_edges) result(frac) + + ! Fraction of the emissivity falling into each of the bins whose edges + ! are given by nu_edges (in ascending order). The fractions are + ! normalized so that they add up to exactly one. + + implicit none + + type(dust),intent(in) :: d + integer,intent(in) :: jnu_var_id + real(dp),intent(in) :: nu_edges(:) + real(dp) :: frac(size(nu_edges)-1) + + integer :: inu + real(dp) :: norm + + do inu=1, size(nu_edges) - 1 + frac(inu) = integral_loglog(d%j_nu(jnu_var_id)%x, d%j_nu(jnu_var_id)%pdf, & + & nu_edges(inu), nu_edges(inu+1)) + end do + + norm = sum(frac) + if(norm > 0._dp) frac = frac / norm + + end function get_j_nu_bin_fractions + function get_chi_nu_interp(d, nu) result(chi_nu) implicit none diff --git a/src/grid/grid_mrw_3d.f90 b/src/grid/grid_mrw_3d.f90 index 6a6e4d26..5f4e4375 100644 --- a/src/grid/grid_mrw_3d.f90 +++ b/src/grid/grid_mrw_3d.f90 @@ -8,6 +8,7 @@ module grid_mrw use grid_physics use dust_main, only : kappa_planck, chi_inv_planck, n_dust, d use type_dust, only : dust_sample_b_nu + use settings, only : compute_specific_energy_spectrum implicit none save @@ -60,9 +61,9 @@ subroutine grid_do_mrw(p) type(photon),intent(inout) :: p real(dp) :: R0 - real(dp) :: e,y,ct + real(dp) :: e,y,ct,fr type(vector3d_dp) :: dr - integer :: id + integer :: id,iv ! Find distance to closest wall R0 = distance_to_closest_wall(p) @@ -79,6 +80,22 @@ subroutine grid_do_mrw(p) ! Insert ct into (9), get energy deposited for Lucy method e = p%energy * ct * kappa_planck(id, specific_energy(p%icell%ic, id)) specific_energy_sum(p%icell%ic, id) = specific_energy_sum(p%icell%ic, id) + e + ! Deposit the same energy in the frequency-resolved spectrum. + ! In the diffusion regime the radiation field is Planckian, so + ! the absorbed energy is distributed in frequency as + ! kappa_nu * B_nu, which is the local emissivity. Distribute the + ! energy over the frequency bins according to that distribution + ! (pre-computed in j_nu_bin_frac), interpolating between the two + ! adjacent emissivity states. No frequency sampling is involved, + ! so computing the spectrum does not affect the random number + ! stream. + if (compute_specific_energy_spectrum) then + iv = jnu_var_id(p%icell%ic, id) + fr = jnu_var_frac(p%icell%ic, id) + specific_energy_sum_spectrum(p%icell%ic, id, :) = & + & specific_energy_sum_spectrum(p%icell%ic, id, :) & + & + e * ((1._dp - fr) * j_nu_bin_frac(:, iv, id) + fr * j_nu_bin_frac(:, iv + 1, id)) + end if end if end do diff --git a/src/grid/grid_physics_3d.f90 b/src/grid/grid_physics_3d.f90 index ffc014a4..1066cf66 100644 --- a/src/grid/grid_physics_3d.f90 +++ b/src/grid/grid_physics_3d.f90 @@ -42,6 +42,11 @@ module grid_physics real(dp),allocatable, public :: nu_bins(:) real(dp),allocatable, public :: log_nu_bins(:) + ! Fraction of the emissivity falling into each frequency bin, for each + ! emissivity state and dust type - used to distribute energy deposited by + ! the MRW over the specific energy spectrum. Indexed (bin, state, dust). + real(dp),allocatable, public :: j_nu_bin_frac(:,:,:) + real(dp),allocatable, public :: specific_energy_additional(:,:) real(dp),allocatable, public :: specific_energy_additional_spectrum(:,:,:) real(dp),allocatable, public :: energy_abs_tot(:) @@ -105,6 +110,8 @@ subroutine setup_grid_physics(group, use_mrw, use_pda, compute_specific_energy_s integer(hid_t),intent(in) :: group logical,intent(in) :: use_mrw, use_pda, compute_specific_energy_spectrum integer :: n_nu_bins + integer :: iv, n_jnu_max + real(dp),allocatable :: nu_bin_edges(:) ! 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. The spectrum @@ -279,6 +286,34 @@ subroutine setup_grid_physics(group, use_mrw, use_pda, compute_specific_energy_s end if allocate(log_nu_bins(n_nu_bins)) log_nu_bins = log10(nu_bins) + + ! Pre-compute, for each dust type and emissivity state, the fraction + ! of the emissivity falling into each frequency bin. This is used to + ! distribute the energy deposited by the MRW over the specific energy + ! spectrum without sampling. The bins are defined by the same + ! nearest-neighbour convention in log-frequency as the photon binning + ! in grid_propagate, so the edges are placed half-way (in log space) + ! between the bin frequencies, with the outer bins extending to all + ! lower/higher frequencies. + allocate(nu_bin_edges(n_nu_bins+1)) + nu_bin_edges(1) = 10._dp**(log_nu_bins(1) - 99._dp) + do idx=2,n_nu_bins + nu_bin_edges(idx) = 10._dp**(0.5_dp * (log_nu_bins(idx-1) + log_nu_bins(idx))) + end do + nu_bin_edges(n_nu_bins+1) = 10._dp**(log_nu_bins(n_nu_bins) + 99._dp) + + n_jnu_max = 0 + do id=1,n_dust + n_jnu_max = max(n_jnu_max, d(id)%n_jnu) + end do + allocate(j_nu_bin_frac(n_nu_bins, n_jnu_max, n_dust)) + j_nu_bin_frac = 0._dp + do id=1,n_dust + do iv=1,d(id)%n_jnu + j_nu_bin_frac(:, iv, id) = get_j_nu_bin_fractions(d(id), iv, nu_bin_edges) + end do + end do + deallocate(nu_bin_edges) end if ! Total energy absorbed From 307eb645a08a0055b5c9e7e78a5ecceb2b81c798 Mon Sep 17 00:00:00 2001 From: Thomas Robitaille Date: Fri, 24 Jul 2026 06:44:04 +0000 Subject: [PATCH 3/4] Rescale the specific energy spectrum by the sublimation cap ratio instead of resetting it to the minimum in sublimation modes 2 and 3 --- src/grid/grid_physics_3d.f90 | 19 +++++++++++-------- 1 file changed, 11 insertions(+), 8 deletions(-) diff --git a/src/grid/grid_physics_3d.f90 b/src/grid/grid_physics_3d.f90 index 1066cf66..12c9f60f 100644 --- a/src/grid/grid_physics_3d.f90 +++ b/src/grid/grid_physics_3d.f90 @@ -419,16 +419,17 @@ subroutine sublimate_dust() & * d(id)%sublimation_specific_energy / specific_energy(ic, id) & & * (chi_rosseland(id, specific_energy(ic,id)) & & / chi_rosseland(id, d(id)%sublimation_specific_energy))**2 - specific_energy(ic,id) = d(id)%sublimation_specific_energy - reset = reset + 1 - if (compute_specific_energy_spectrum) then do idx=1,n_nu_bins - specific_energy_spectrum(ic,id,idx) = minimum_specific_energy(id) + specific_energy_spectrum(ic,id,idx) = specific_energy_spectrum(ic,id,idx) & + & * d(id)%sublimation_specific_energy / specific_energy(ic, id) end do end if - + + specific_energy(ic,id) = d(id)%sublimation_specific_energy + reset = reset + 1 + end if end do @@ -438,15 +439,17 @@ subroutine sublimate_dust() do ic=1,geo%n_cells if(specific_energy(ic, id) > d(id)%sublimation_specific_energy) then - specific_energy(ic, id) = d(id)%sublimation_specific_energy - reset = reset + 1 if (compute_specific_energy_spectrum) then do idx=1,n_nu_bins - specific_energy_spectrum(ic,id,idx) = minimum_specific_energy(id) + specific_energy_spectrum(ic,id,idx) = specific_energy_spectrum(ic,id,idx) & + & * d(id)%sublimation_specific_energy / specific_energy(ic, id) end do end if + specific_energy(ic, id) = d(id)%sublimation_specific_energy + reset = reset + 1 + end if end do From 48b5b83fc5afa6fcb805a93f519e0ee8afad7029 Mon Sep 17 00:00:00 2001 From: Thomas Robitaille Date: Tue, 28 Jul 2026 22:07:43 +0000 Subject: [PATCH 4/4] Move the specific energy spectrum bin fraction setup, scaling, and MRW deposition into helper subroutines --- src/grid/grid_mrw_3d.f90 | 25 ++----- src/grid/grid_pda_3d.f90 | 7 +- src/grid/grid_physics_3d.f90 | 135 ++++++++++++++++++++++++----------- 3 files changed, 101 insertions(+), 66 deletions(-) diff --git a/src/grid/grid_mrw_3d.f90 b/src/grid/grid_mrw_3d.f90 index 5f4e4375..1ac3c1d8 100644 --- a/src/grid/grid_mrw_3d.f90 +++ b/src/grid/grid_mrw_3d.f90 @@ -8,7 +8,6 @@ module grid_mrw use grid_physics use dust_main, only : kappa_planck, chi_inv_planck, n_dust, d use type_dust, only : dust_sample_b_nu - use settings, only : compute_specific_energy_spectrum implicit none save @@ -61,9 +60,9 @@ subroutine grid_do_mrw(p) type(photon),intent(inout) :: p real(dp) :: R0 - real(dp) :: e,y,ct,fr + real(dp) :: e,y,ct type(vector3d_dp) :: dr - integer :: id,iv + integer :: id ! Find distance to closest wall R0 = distance_to_closest_wall(p) @@ -80,22 +79,10 @@ subroutine grid_do_mrw(p) ! Insert ct into (9), get energy deposited for Lucy method e = p%energy * ct * kappa_planck(id, specific_energy(p%icell%ic, id)) specific_energy_sum(p%icell%ic, id) = specific_energy_sum(p%icell%ic, id) + e - ! Deposit the same energy in the frequency-resolved spectrum. - ! In the diffusion regime the radiation field is Planckian, so - ! the absorbed energy is distributed in frequency as - ! kappa_nu * B_nu, which is the local emissivity. Distribute the - ! energy over the frequency bins according to that distribution - ! (pre-computed in j_nu_bin_frac), interpolating between the two - ! adjacent emissivity states. No frequency sampling is involved, - ! so computing the spectrum does not affect the random number - ! stream. - if (compute_specific_energy_spectrum) then - iv = jnu_var_id(p%icell%ic, id) - fr = jnu_var_frac(p%icell%ic, id) - specific_energy_sum_spectrum(p%icell%ic, id, :) = & - & specific_energy_sum_spectrum(p%icell%ic, id, :) & - & + e * ((1._dp - fr) * j_nu_bin_frac(:, iv, id) + fr * j_nu_bin_frac(:, iv + 1, id)) - end if + ! Deposit the same energy in the frequency-resolved spectrum, + ! distributed over the frequency bins according to the local + ! emissivity (see deposit_specific_energy_spectrum). + call deposit_specific_energy_spectrum(p%icell%ic, id, e) end if end do diff --git a/src/grid/grid_pda_3d.f90 b/src/grid/grid_pda_3d.f90 index e23dc839..925b781b 100644 --- a/src/grid/grid_pda_3d.f90 +++ b/src/grid/grid_pda_3d.f90 @@ -11,7 +11,6 @@ module grid_pda use dust_main use grid_io use grid_pda_geometry - use settings, only : compute_specific_energy_spectrum implicit none save @@ -65,11 +64,7 @@ subroutine update_specific_energy(ic) ! scalar value by rescaling the spectrum bins, preserving their shape. ! If the cell had no energy before, the bins are all zero and there is ! no shape to rescale, so they are left unchanged. - if(compute_specific_energy_spectrum) then - if(s_old > 0._dp) then - specific_energy_spectrum(ic, id, :) = specific_energy_spectrum(ic, id, :) * (s / s_old) - end if - end if + if(s_old > 0._dp) call scale_specific_energy_spectrum(ic, id, s / s_old) end do end subroutine update_specific_energy diff --git a/src/grid/grid_physics_3d.f90 b/src/grid/grid_physics_3d.f90 index 12c9f60f..3a486186 100644 --- a/src/grid/grid_physics_3d.f90 +++ b/src/grid/grid_physics_3d.f90 @@ -16,6 +16,8 @@ module grid_physics private public :: setup_grid_physics + public :: scale_specific_energy_spectrum + public :: deposit_specific_energy_spectrum public :: sublimate_dust public :: update_alpha_inv_planck public :: check_energy_abs @@ -110,8 +112,6 @@ subroutine setup_grid_physics(group, use_mrw, use_pda, compute_specific_energy_s integer(hid_t),intent(in) :: group logical,intent(in) :: use_mrw, use_pda, compute_specific_energy_spectrum integer :: n_nu_bins - integer :: iv, n_jnu_max - real(dp),allocatable :: nu_bin_edges(:) ! 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. The spectrum @@ -287,33 +287,7 @@ subroutine setup_grid_physics(group, use_mrw, use_pda, compute_specific_energy_s allocate(log_nu_bins(n_nu_bins)) log_nu_bins = log10(nu_bins) - ! Pre-compute, for each dust type and emissivity state, the fraction - ! of the emissivity falling into each frequency bin. This is used to - ! distribute the energy deposited by the MRW over the specific energy - ! spectrum without sampling. The bins are defined by the same - ! nearest-neighbour convention in log-frequency as the photon binning - ! in grid_propagate, so the edges are placed half-way (in log space) - ! between the bin frequencies, with the outer bins extending to all - ! lower/higher frequencies. - allocate(nu_bin_edges(n_nu_bins+1)) - nu_bin_edges(1) = 10._dp**(log_nu_bins(1) - 99._dp) - do idx=2,n_nu_bins - nu_bin_edges(idx) = 10._dp**(0.5_dp * (log_nu_bins(idx-1) + log_nu_bins(idx))) - end do - nu_bin_edges(n_nu_bins+1) = 10._dp**(log_nu_bins(n_nu_bins) + 99._dp) - - n_jnu_max = 0 - do id=1,n_dust - n_jnu_max = max(n_jnu_max, d(id)%n_jnu) - end do - allocate(j_nu_bin_frac(n_nu_bins, n_jnu_max, n_dust)) - j_nu_bin_frac = 0._dp - do id=1,n_dust - do iv=1,d(id)%n_jnu - j_nu_bin_frac(:, iv, id) = get_j_nu_bin_fractions(d(id), iv, nu_bin_edges) - end do - end do - deallocate(nu_bin_edges) + call setup_j_nu_bin_fractions() end if ! Total energy absorbed @@ -354,6 +328,93 @@ subroutine setup_grid_physics(group, use_mrw, use_pda, compute_specific_energy_s end subroutine setup_grid_physics + subroutine setup_j_nu_bin_fractions() + + ! Pre-compute, for each dust type and emissivity state, the fraction of + ! the emissivity falling into each frequency bin, stored in + ! j_nu_bin_frac. This is used to distribute the energy deposited by the + ! MRW over the specific energy spectrum without sampling. The bins are + ! defined by the same nearest-neighbour convention in log-frequency as + ! the photon binning in grid_propagate, so the edges are placed half-way + ! (in log space) between the bin frequencies, with the outer bins + ! extending to all lower/higher frequencies. + + implicit none + + integer :: inu, id, iv, n_nu_bins, n_jnu_max + real(dp),allocatable :: nu_bin_edges(:) + + n_nu_bins = size(log_nu_bins) + + allocate(nu_bin_edges(n_nu_bins+1)) + nu_bin_edges(1) = 10._dp**(log_nu_bins(1) - 99._dp) + do inu=2,n_nu_bins + nu_bin_edges(inu) = 10._dp**(0.5_dp * (log_nu_bins(inu-1) + log_nu_bins(inu))) + end do + nu_bin_edges(n_nu_bins+1) = 10._dp**(log_nu_bins(n_nu_bins) + 99._dp) + + n_jnu_max = 0 + do id=1,n_dust + n_jnu_max = max(n_jnu_max, d(id)%n_jnu) + end do + allocate(j_nu_bin_frac(n_nu_bins, n_jnu_max, n_dust)) + j_nu_bin_frac = 0._dp + do id=1,n_dust + do iv=1,d(id)%n_jnu + j_nu_bin_frac(:, iv, id) = get_j_nu_bin_fractions(d(id), iv, nu_bin_edges) + end do + end do + deallocate(nu_bin_edges) + + end subroutine setup_j_nu_bin_fractions + + subroutine scale_specific_energy_spectrum(ic, id, factor) + + ! Multiply the frequency-resolved specific energy of a cell and dust + ! type by the given factor, preserving the spectral shape. Does nothing + ! if the frequency-resolved specific energy is not being computed. + + implicit none + + integer,intent(in) :: ic, id + real(dp),intent(in) :: factor + + if(.not.compute_specific_energy_spectrum) return + + specific_energy_spectrum(ic, id, :) = specific_energy_spectrum(ic, id, :) * factor + + end subroutine scale_specific_energy_spectrum + + subroutine deposit_specific_energy_spectrum(ic, id, energy) + + ! Add energy to the frequency-resolved specific energy of a cell and + ! dust type, distributed over the frequency bins according to the local + ! emissivity (pre-computed in j_nu_bin_frac), interpolating between the + ! two adjacent emissivity states. This is used by the MRW, for which + ! the radiation field is Planckian, so the absorbed energy is + ! distributed in frequency as kappa_nu * B_nu, which is the local + ! emissivity. No frequency sampling is involved, so computing the + ! spectrum does not affect the random number stream. Does nothing if + ! the frequency-resolved specific energy is not being computed. + + implicit none + + integer,intent(in) :: ic, id + real(dp),intent(in) :: energy + + integer :: iv + real(dp) :: fr + + if(.not.compute_specific_energy_spectrum) return + + iv = jnu_var_id(ic, id) + fr = jnu_var_frac(ic, id) + + specific_energy_sum_spectrum(ic, id, :) = specific_energy_sum_spectrum(ic, id, :) & + & + energy * ((1._dp - fr) * j_nu_bin_frac(:, iv, id) + fr * j_nu_bin_frac(:, iv + 1, id)) + + end subroutine deposit_specific_energy_spectrum + subroutine update_alpha_inv_planck() ! Optimization: could pre-compute alpha_inv_planck just for masked (valid) cells @@ -420,12 +481,8 @@ subroutine sublimate_dust() & * (chi_rosseland(id, specific_energy(ic,id)) & & / chi_rosseland(id, d(id)%sublimation_specific_energy))**2 - if (compute_specific_energy_spectrum) then - do idx=1,n_nu_bins - specific_energy_spectrum(ic,id,idx) = specific_energy_spectrum(ic,id,idx) & - & * d(id)%sublimation_specific_energy / specific_energy(ic, id) - end do - end if + call scale_specific_energy_spectrum(ic, id, & + & d(id)%sublimation_specific_energy / specific_energy(ic, id)) specific_energy(ic,id) = d(id)%sublimation_specific_energy reset = reset + 1 @@ -440,12 +497,8 @@ subroutine sublimate_dust() do ic=1,geo%n_cells if(specific_energy(ic, id) > d(id)%sublimation_specific_energy) then - if (compute_specific_energy_spectrum) then - do idx=1,n_nu_bins - specific_energy_spectrum(ic,id,idx) = specific_energy_spectrum(ic,id,idx) & - & * d(id)%sublimation_specific_energy / specific_energy(ic, id) - end do - end if + call scale_specific_energy_spectrum(ic, id, & + & d(id)%sublimation_specific_energy / specific_energy(ic, id)) specific_energy(ic, id) = d(id)%sublimation_specific_energy reset = reset + 1