From 711b282b47a90bce30d2cfdb7971b42ed4326e78 Mon Sep 17 00:00:00 2001 From: dnarayanan Date: Mon, 4 Oct 2021 09:34:48 -0400 Subject: [PATCH 01/38] strawperson architecture in for computing specific energy as a function of frequency. note, this has a manual setting of frequencies and right now only computes, not saves --- src/grid/grid_physics_3d.f90 | 5 ++++ src/grid/grid_propagate_3d.f90 | 55 +++++++++++++++++++++++++++++++++- 2 files changed, 59 insertions(+), 1 deletion(-) diff --git a/src/grid/grid_physics_3d.f90 b/src/grid/grid_physics_3d.f90 index 90674116..429f041d 100644 --- a/src/grid/grid_physics_3d.f90 +++ b/src/grid/grid_physics_3d.f90 @@ -37,6 +37,8 @@ module grid_physics integer(idp),allocatable, public :: last_photon_id(:) real(dp),allocatable, public :: specific_energy(:,:) real(dp),allocatable, public :: specific_energy_sum(:,:) + real(dp),allocatable, public :: specific_energy_sum_nu(:,:,:) + real(dp),allocatable, public :: specific_energy_additional(:,:) real(dp),allocatable, public :: energy_abs_tot(:) real(dp),allocatable, public :: minimum_specific_energy(:) @@ -191,6 +193,9 @@ subroutine setup_grid_physics(group, use_mrw, use_pda) allocate(specific_energy_sum(geo%n_cells, n_dust)) specific_energy_sum = 0._dp + allocate(specific_energy_sum_nu(geo%n_cells, n_dust, 10)) + specific_energy_sum_nu = 0._dp + ! Total energy absorbed allocate(energy_abs_tot(n_dust)) energy_abs_tot = 0._dp diff --git a/src/grid/grid_propagate_3d.f90 b/src/grid/grid_propagate_3d.f90 index 3c345161..d2f2a723 100644 --- a/src/grid/grid_propagate_3d.f90 +++ b/src/grid/grid_propagate_3d.f90 @@ -5,7 +5,7 @@ module grid_propagate use type_grid_cell use dust_main, only : n_dust use grid_geometry, only : escaped, find_wall, in_correct_cell, next_cell, opposite_wall - use grid_physics, only : specific_energy_sum, density, n_photons, last_photon_id + use grid_physics, only : specific_energy_sum, specific_energy_sum_nu, density, n_photons, last_photon_id use sources use counters use settings, only : frac_check => propagation_check_frequency @@ -131,12 +131,30 @@ subroutine grid_integrate(p,tau_required,tau_achieved) tau_achieved = tau_achieved + tau_cell ! NEED INDIVIDUAL ALPHA HERE + + !print *,'[grid_propagate_3d] p%energy=',p%energy + !print *,'[grid_propagate_3d] p%nu=',p%nu + + !DN CRAZY ADDITIONS + !print *,'[grid_propagate_3d] p%nu=',p%nu + !print *,' ',minloc(abs(energy_frequency_bins-p%nu)) + + !DN CRAZY ADDITIONS + !idx = minloc(abs(energy_frequency_bins-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) = & & specific_energy_sum(p%icell%ic, id) + tmin * p%current_kappa(id) * p%energy end if + + !if(density(p%icell%ic,id) > 0._dp) then + + + !specific_energy_sum_nu(p%icell%ic,id,1) = & + ! & specific_energy_sum_nu(p%icell%ic,id,idx) + tmin * p%current_kappa(id) * p%energy + + !end if end do p%on_wall = .true. @@ -221,6 +239,22 @@ subroutine grid_integrate_noenergy(p,tau_required,tau_achieved) integer :: source_id + !DN CRAZY ADDITIONS + integer :: idx + !DN CRAZY ADDITIONS + real, dimension(10) :: energy_frequency_bins + energy_frequency_bins(1) = 10.**14.4768207 + energy_frequency_bins(2) = 10.**14.59237683 + energy_frequency_bins(3) = 10.**14.70793296 + energy_frequency_bins(4) = 10.**14.82348909 + energy_frequency_bins(5) = 10.**14.93904522 + energy_frequency_bins(6) = 10.**15.05460135 + energy_frequency_bins(7) = 10.**15.17015748 + energy_frequency_bins(8) = 10.**15.28571361 + energy_frequency_bins(9) = 10.**15.40126974 + energy_frequency_bins(10) = 10.**15.51682586 + + radial = (p%r .dot. p%v) > 0. if(debug) write(*,'(" [debug] start grid_integrate_noenergy")') @@ -268,9 +302,28 @@ subroutine grid_integrate_noenergy(p,tau_required,tau_achieved) end if chi_rho_total = 0._dp + + + + !DN CRAZY ADDITIONS + !idx = minloc(abs(energy_frequency_bins-p%nu),DIM=1) + !print *,'[grid_propagate_3d last iteration] p%energy=',p%energy + !print *,'[grid_propagate_3d last iteration] p%nu=',p%nu + idx = minloc(abs(energy_frequency_bins-p%nu),DIM=1) + do id=1,n_dust chi_rho_total = chi_rho_total + p%current_chi(id) * density(p%icell%ic, id) + + !DN CRAZY ADDITIONS + if(density(p%icell%ic,id) > 0._dp) then + specific_energy_sum_nu(p%icell%ic,id,1) = & + & specific_energy_sum_nu(p%icell%ic,id,idx) + tmin * p%current_kappa(id) * p%energy + print *,'[grid_propage_3d last iteration] specific_energy_sum_nu=',specific_energy_sum_nu + end if + end do + + tau_cell = chi_rho_total * tmin if(tau_cell < tau_needed) then From 258972c2d9b51c5dde9dfd52301690187c4f8128 Mon Sep 17 00:00:00 2001 From: dnarayanan Date: Wed, 10 Nov 2021 10:27:49 -0500 Subject: [PATCH 02/38] specific_energy_sum_nu is in all places specific_energy_sum is now i think --- src/grid/grid_generic.f90 | 24 +++++++++++++++++++++++- src/grid/grid_propagate_3d.f90 | 2 +- src/mpi/mpi_routines.f90 | 13 +++++++++++++ 3 files changed, 37 insertions(+), 2 deletions(-) diff --git a/src/grid/grid_generic.f90 b/src/grid/grid_generic.f90 index d050bca0..29e19179 100644 --- a/src/grid/grid_generic.f90 +++ b/src/grid/grid_generic.f90 @@ -6,7 +6,7 @@ module grid_generic use grid_io, only : write_grid_3d, write_grid_4d use grid_geometry, only : geo - use grid_physics, only : n_photons, last_photon_id, specific_energy_sum, specific_energy, density, density_original + use grid_physics, only : n_photons, last_photon_id, specific_energy_sum, specific_energy_sum_nu, specific_energy, density, density_original use settings, only : output_n_photons, output_specific_energy, output_density, output_density_diff, physics_io_type implicit none @@ -23,6 +23,7 @@ subroutine grid_reset_energy() if(allocated(n_photons)) n_photons = 0 if(allocated(last_photon_id)) last_photon_id = 0 if(allocated(specific_energy_sum)) specific_energy_sum = 0._dp + if(allocated(specific_energy_sum_nu)) specific_energy_sum_nu = 0._dp end subroutine grid_reset_energy subroutine output_grid(group, iter, n_iter) @@ -61,6 +62,27 @@ subroutine output_grid(group, iter, n_iter) end if end if + + + ! ENERGY/PATH LENGTHS + + if(trim(output_specific_energy)=='all' .or. (trim(output_specific_energy)=='last'.and.iter==n_iter)) then + if(allocated(specific_energy)) then + select case(physics_io_type) + case(sp) + call write_grid_4d(group, 'specific_energy_sum_nu', real(specific_energy_sum, sp), geo) + case(dp) + call write_grid_4d(group, 'specific_energy_sum_nu', real(specific_energy_sum, dp), geo) + case default + call error("output_grid","unexpected value of physics_io_type (should be sp or dp)") + end select + else + call warn("output_grid","specific_energy array is not allocated") + end if + end if + + + ! DENSITY if(trim(output_density)=='all' .or. (trim(output_density)=='last'.and.iter==n_iter)) then diff --git a/src/grid/grid_propagate_3d.f90 b/src/grid/grid_propagate_3d.f90 index d2f2a723..4a04f6f9 100644 --- a/src/grid/grid_propagate_3d.f90 +++ b/src/grid/grid_propagate_3d.f90 @@ -318,7 +318,7 @@ subroutine grid_integrate_noenergy(p,tau_required,tau_achieved) if(density(p%icell%ic,id) > 0._dp) then specific_energy_sum_nu(p%icell%ic,id,1) = & & specific_energy_sum_nu(p%icell%ic,id,idx) + tmin * p%current_kappa(id) * p%energy - print *,'[grid_propage_3d last iteration] specific_energy_sum_nu=',specific_energy_sum_nu + !print *,'[grid_propage_3d last iteration] specific_energy_sum_nu=',specific_energy_sum_nu end if end do diff --git a/src/mpi/mpi_routines.f90 b/src/mpi/mpi_routines.f90 index fb0d6e7b..8f4e8a15 100644 --- a/src/mpi/mpi_routines.f90 +++ b/src/mpi/mpi_routines.f90 @@ -266,6 +266,8 @@ subroutine mp_collect_physical_arrays() real(dp) :: tmp real(dp) :: dummy_dp integer(idp) :: dummy_idp + !DN crazy additions + real(dp),allocatable :: tmp_3d(:,:,:) real(dp),allocatable :: tmp_2d(:,:) integer(idp),allocatable :: tmp_int_1d(:) @@ -278,6 +280,17 @@ 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_nu,1),size(specific_energy_sum_nu,2),size(specific_energy_sum_nu,3))) + call mpi_reduce(specific_energy_sum_nu, tmp_3d, size(specific_energy_sum_nu), mpi_real8, mpi_sum, rank_main, mpi_comm_world, ierr) + specific_energy_sum_nu = tmp_3d + deallocate(tmp_3d) + else + call mpi_reduce(specific_energy_sum_nu, dummy_dp, size(specific_energy_sum_nu), mpi_real8, mpi_sum, rank_main, mpi_comm_world, ierr) + end if + + + if(allocated(n_photons)) then if(main_process()) then allocate(tmp_int_1d(size(n_photons,1))) From 8b91d3c6064919575f097669e0b233d3f03dfa17 Mon Sep 17 00:00:00 2001 From: dnarayanan Date: Tue, 23 Nov 2021 09:36:38 -0500 Subject: [PATCH 03/38] code now saves a ISRF in every cell. notably: AMR is still broken, and m.read is broken still with the new specific_energy_sum_nu array dimensions --- hyperion/conf/conf_files.py | 3 + hyperion/model/model.py | 4 + src/grid/grid_generic.f90 | 31 ++++- src/grid/grid_io.f90 | 212 ++++++++++++++++++++++++++++++ src/grid/grid_io_1d.f90 | 178 +++++++++++++++++++++++++ src/grid/grid_io_amr.f90 | 181 +++++++++++++++++++++++++ src/grid/grid_io_amr_template.f90 | 53 ++++++++ src/grid/grid_io_template.f90 | 66 ++++++++++ src/grid/grid_physics_3d.f90 | 2 +- src/grid/grid_propagate_3d.f90 | 45 +++++-- 10 files changed, 760 insertions(+), 15 deletions(-) diff --git a/hyperion/conf/conf_files.py b/hyperion/conf/conf_files.py index 59c62a64..dac010d7 100644 --- a/hyperion/conf/conf_files.py +++ b/hyperion/conf/conf_files.py @@ -18,6 +18,7 @@ def __init__(self): self.output_density = 'none' self.output_density_diff = 'none' self.output_specific_energy = 'last' + self.output_specific_energy_nu = 'last' self.output_n_photons = 'none' self._freeze() @@ -27,6 +28,7 @@ def read(cls, group): self.output_density = group.attrs['output_density'].decode('utf-8') self.output_density_diff = group.attrs['output_density_diff'].decode('utf-8') self.output_specific_energy = group.attrs['output_specific_energy'].decode('utf-8') + self.output_specific_energy_nu = group.attrs['output_specific_energy_nu'].decode('utf-8') self.output_n_photons = group.attrs['output_n_photons'].decode('utf-8') return self @@ -34,6 +36,7 @@ def write(self, group): group.attrs['output_density'] = np.bytes_(self.output_density.encode('utf-8')) group.attrs['output_density_diff'] = np.bytes_(self.output_density_diff.encode('utf-8')) group.attrs['output_specific_energy'] = np.bytes_(self.output_specific_energy.encode('utf-8')) + group.attrs['output_specific_energy_nu'] = np.bytes_(self.output_specific_energ_nu.encode('utf-8')) group.attrs['output_n_photons'] = np.bytes_(self.output_n_photons.encode('utf-8')) diff --git a/hyperion/model/model.py b/hyperion/model/model.py index c781ee7b..23a6edba 100644 --- a/hyperion/model/model.py +++ b/hyperion/model/model.py @@ -313,6 +313,10 @@ def use_quantities(self, filename, quantities=['density', 'specific_energy'], if 'specific_energy' in f['/Grid/Quantities']: quantities_path['specific_energy'] = '/Grid/Quantities' + if 'specific_energy_nu' in quantities: + if 'specific_energy_nu' in f['/Grid/Quantities']: + quantities_path['specific_energy_nu'] = '/Grid/Quantities' + # Minimum specific energy if use_minimum_specific_energy: minimum_specific_energy_path = '/Grid/Quantities' diff --git a/src/grid/grid_generic.f90 b/src/grid/grid_generic.f90 index 29e19179..14397b6a 100644 --- a/src/grid/grid_generic.f90 +++ b/src/grid/grid_generic.f90 @@ -4,7 +4,7 @@ module grid_generic use mpi_hdf5_io, only : mp_create_group use mpi_core, only : main_process - use grid_io, only : write_grid_3d, write_grid_4d + use grid_io, only : write_grid_3d, write_grid_4d, write_grid_5d use grid_geometry, only : geo use grid_physics, only : n_photons, last_photon_id, specific_energy_sum, specific_energy_sum_nu, specific_energy, density, density_original use settings, only : output_n_photons, output_specific_energy, output_density, output_density_diff, physics_io_type @@ -32,6 +32,9 @@ subroutine output_grid(group, iter, n_iter) integer(hid_t),intent(in) :: group integer,intent(in) :: iter, n_iter + !DN CRAZY ADDITION + integer :: n_cells, n_dust, n_isrf_lam + integer :: i,j,k if(main_process()) write(*,'(" [output_grid] outputting grid arrays for iteration")') @@ -64,20 +67,36 @@ subroutine output_grid(group, iter, n_iter) - ! ENERGY/PATH LENGTHS + ! DN Crazy Additions + n_cells = size(specific_energy_sum_nu, 1) + n_dust = size(specific_energy_sum_nu, 2) + n_isrf_lam = size(specific_energy_sum_nu,3) + do i=1,n_cells + do j=1,n_dust + !print *, "specific_energy_sum = ",specific_energy_sum(i,j) + do k=1,n_isrf_lam + print *, "specific_energy_sum_nu = ",specific_energy_sum_nu(i,j,k) + end do + end do + end do + + if(trim(output_specific_energy)=='all' .or. (trim(output_specific_energy)=='last'.and.iter==n_iter)) then - if(allocated(specific_energy)) then + if(allocated(specific_energy_sum_nu)) then + !print *, "shape(specific_energy_sum_nu", shape(specific_energy_sum_nu) + !print *, "specific_energy_sum_nu = ",specific_energy_sum_nu select case(physics_io_type) + case(sp) - call write_grid_4d(group, 'specific_energy_sum_nu', real(specific_energy_sum, sp), geo) + call write_grid_5d(group, 'specific_energy_sum_nu', real(specific_energy_sum_nu, sp), geo) case(dp) - call write_grid_4d(group, 'specific_energy_sum_nu', real(specific_energy_sum, dp), geo) + call write_grid_5d(group, 'specific_energy_sum_nu', real(specific_energy_sum_nu, dp), geo) case default call error("output_grid","unexpected value of physics_io_type (should be sp or dp)") end select else - call warn("output_grid","specific_energy array is not allocated") + call warn("output_grid","specific_energy_sum_nu array is not allocated") end if end if diff --git a/src/grid/grid_io.f90 b/src/grid/grid_io.f90 index 8cfc82f5..a4aaee28 100644 --- a/src/grid/grid_io.f90 +++ b/src/grid/grid_io.f90 @@ -12,8 +12,10 @@ module grid_io public :: grid_exists public :: read_grid_3d public :: read_grid_4d + public :: read_grid_5d public :: write_grid_3d public :: write_grid_4d + public :: write_grid_5d interface read_grid_3d module procedure read_grid_3d_sp @@ -29,6 +31,14 @@ module grid_io module procedure read_grid_4d_int8 end interface read_grid_4d + interface read_grid_5d + module procedure read_grid_5d_sp + module procedure read_grid_5d_dp + module procedure read_grid_5d_int + module procedure read_grid_5d_int8 + end interface read_grid_5d + + interface write_grid_3d module procedure write_grid_3d_sp module procedure write_grid_3d_dp @@ -43,6 +53,14 @@ module grid_io module procedure write_grid_4d_int8 end interface write_grid_4d + interface write_grid_5d + module procedure write_grid_5d_sp + module procedure write_grid_5d_dp + module procedure write_grid_5d_int + module procedure write_grid_5d_int8 + end interface write_grid_5d + + contains logical function grid_exists(group, name) @@ -52,6 +70,36 @@ logical function grid_exists(group, name) grid_exists = mp_path_exists(group, name) end function grid_exists + !DN CRAZY ADDITIONS + subroutine read_grid_5d_int8(group, path, array, geo) + + implicit none + + integer(hid_t), intent(in) :: group + character(len=*), intent(in) :: path + integer(idp), intent(out) :: array(:,:,:) + type(grid_geometry_desc),intent(in) :: geo + integer(idp), allocatable :: array5d(:,:,:,:,:) + integer :: n_cells, n_dust, n_isrf_lam + + character(len=32) :: geometry_id_check + + call mp_read_keyword(group,path, 'geometry', geometry_id_check) + if(geometry_id_check.ne.geo%id) then + call error("read_grid", "geometry IDs do not match") + end if + call mp_read_array_auto(group,path, array5d) + + if(any(is_nan(array5d))) call error("read_grid_5d", "NaN values in 5D array") + + n_cells = size(array, 1) + n_dust = size(array, 2) + n_isrf_lam = size(array,3) + + array = reshape(array5d, (/n_cells, n_dust, n_isrf_lam/)) + + end subroutine read_grid_5d_int8 + subroutine read_grid_4d_int8(group, path, array, geo) @@ -108,6 +156,26 @@ subroutine read_grid_3d_int8(group, path, array, geo) end subroutine read_grid_3d_int8 + + + !DN CRAZY ADDITIONS + subroutine write_grid_5d_int8(group, path, array, geo) + + implicit none + + integer(hid_t), intent(in) :: group + character(len=*), intent(in) :: path + integer(idp), intent(in) :: array(:,:,:) + type(grid_geometry_desc),intent(in) :: geo + + call mp_write_array(group, path, reshape(array, (/geo%n1, geo%n2, geo%n3, size(array,2), size(array,3)/))) + call mp_write_keyword(group, path, 'geometry', geo%id) + + end subroutine write_grid_5d_int8 + + + + subroutine write_grid_4d_int8(group, path, array, geo) implicit none @@ -136,6 +204,37 @@ subroutine write_grid_3d_int8(group, path, array, geo) end subroutine write_grid_3d_int8 + !DN CRAZY ADDITIONS + subroutine read_grid_5d_int(group, path, array, geo) + + implicit none + + integer(hid_t), intent(in) :: group + character(len=*), intent(in) :: path + integer, intent(out) :: array(:,:,:) + type(grid_geometry_desc),intent(in) :: geo + integer, allocatable :: array5d(:,:,:,:,:) + integer :: n_cells, n_dust, n_isrf_lam + + character(len=32) :: geometry_id_check + + call mp_read_keyword(group,path, 'geometry', geometry_id_check) + if(geometry_id_check.ne.geo%id) then + call error("read_grid", "geometry IDs do not match") + end if + call mp_read_array_auto(group,path, array5d) + + if(any(is_nan(array5d))) call error("read_grid_5d", "NaN values in 5D array") + + n_cells = size(array, 1) + n_dust = size(array, 2) + n_isrf_lam = size(array,3) + + array = reshape(array5d, (/n_cells, n_dust, n_isrf_lam/)) + + end subroutine read_grid_5d_int + + subroutine read_grid_4d_int(group, path, array, geo) @@ -192,6 +291,24 @@ subroutine read_grid_3d_int(group, path, array, geo) end subroutine read_grid_3d_int + + !DN CRAZY ADDITION + + subroutine write_grid_5d_int(group, path, array, geo) + + implicit none + + integer(hid_t), intent(in) :: group + character(len=*), intent(in) :: path + integer, intent(in) :: array(:,:,:) + type(grid_geometry_desc),intent(in) :: geo + + call mp_write_array(group, path, reshape(array, (/geo%n1, geo%n2, geo%n3, size(array,2), size(array,3)/))) + call mp_write_keyword(group, path, 'geometry', geo%id) + + end subroutine write_grid_5d_int + + subroutine write_grid_4d_int(group, path, array, geo) implicit none @@ -220,6 +337,36 @@ subroutine write_grid_3d_int(group, path, array, geo) end subroutine write_grid_3d_int + !DN CRAZY ADDITIONS + subroutine read_grid_5d_dp(group, path, array, geo) + + implicit none + + integer(hid_t), intent(in) :: group + character(len=*), intent(in) :: path + real(dp), intent(out) :: array(:,:,:) + type(grid_geometry_desc),intent(in) :: geo + real(dp), allocatable :: array5d(:,:,:,:,:) + integer :: n_cells, n_dust, n_isrf_lam + + character(len=32) :: geometry_id_check + + call mp_read_keyword(group,path, 'geometry', geometry_id_check) + if(geometry_id_check.ne.geo%id) then + call error("read_grid", "geometry IDs do not match") + end if + call mp_read_array_auto(group,path, array5d) + + if(any(is_nan(array5d))) call error("read_grid_5d", "NaN values in 5D array") + + n_cells = size(array, 1) + n_dust = size(array, 2) + n_isrf_lam = size(array,3) + + array = reshape(array5d, (/n_cells, n_dust, n_isrf_lam/)) + + end subroutine read_grid_5d_dp + subroutine read_grid_4d_dp(group, path, array, geo) @@ -276,6 +423,23 @@ subroutine read_grid_3d_dp(group, path, array, geo) end subroutine read_grid_3d_dp + + !CRAZY DN ADDITIONS + subroutine write_grid_5d_dp(group, path, array, geo) + + implicit none + + integer(hid_t), intent(in) :: group + character(len=*), intent(in) :: path + real(dp), intent(in) :: array(:,:,:) + type(grid_geometry_desc),intent(in) :: geo + + call mp_write_array(group, path, reshape(array, (/geo%n1, geo%n2, geo%n3, size(array,2), size(array,3)/))) + call mp_write_keyword(group, path, 'geometry', geo%id) + + end subroutine write_grid_5d_dp + + subroutine write_grid_4d_dp(group, path, array, geo) implicit none @@ -305,6 +469,37 @@ subroutine write_grid_3d_dp(group, path, array, geo) end subroutine write_grid_3d_dp + !DN CRAZY ADDITIONS + subroutine read_grid_5d_sp(group, path, array, geo) + + implicit none + + integer(hid_t), intent(in) :: group + character(len=*), intent(in) :: path + real(sp), intent(out) :: array(:,:,:) + type(grid_geometry_desc),intent(in) :: geo + real(sp), allocatable :: array5d(:,:,:,:,:) + integer :: n_cells, n_dust, n_isrf_lam + + character(len=32) :: geometry_id_check + + call mp_read_keyword(group,path, 'geometry', geometry_id_check) + if(geometry_id_check.ne.geo%id) then + call error("read_grid", "geometry IDs do not match") + end if + call mp_read_array_auto(group,path, array5d) + + if(any(is_nan(array5d))) call error("read_grid_5d", "NaN values in 5D array") + + n_cells = size(array, 1) + n_dust = size(array, 2) + n_isrf_lam = size(array,3) + + array = reshape(array5d, (/n_cells, n_dust, n_isrf_lam/)) + + end subroutine read_grid_5d_sp + + subroutine read_grid_4d_sp(group, path, array, geo) implicit none @@ -360,6 +555,23 @@ subroutine read_grid_3d_sp(group, path, array, geo) end subroutine read_grid_3d_sp + + !CRAZY DN ADDITIONS + subroutine write_grid_5d_sp(group, path, array, geo) + + implicit none + + integer(hid_t), intent(in) :: group + character(len=*), intent(in) :: path + real(sp), intent(in) :: array(:,:,:) + type(grid_geometry_desc),intent(in) :: geo + + call mp_write_array(group, path, reshape(array, (/geo%n1, geo%n2, geo%n3, size(array,2), size(array,3)/))) + call mp_write_keyword(group, path, 'geometry', geo%id) + + end subroutine write_grid_5d_sp + + subroutine write_grid_4d_sp(group, path, array, geo) implicit none diff --git a/src/grid/grid_io_1d.f90 b/src/grid/grid_io_1d.f90 index 5d87aaa2..98779162 100644 --- a/src/grid/grid_io_1d.f90 +++ b/src/grid/grid_io_1d.f90 @@ -12,8 +12,11 @@ module grid_io public :: grid_exists public :: read_grid_3d public :: read_grid_4d + public :: read_grid_5d public :: write_grid_3d public :: write_grid_4d + public :: write_grid_5d + interface read_grid_3d module procedure read_grid_3d_sp @@ -29,6 +32,14 @@ module grid_io module procedure read_grid_4d_int8 end interface read_grid_4d + interface read_grid_5d + module procedure read_grid_5d_sp + module procedure read_grid_5d_dp + module procedure read_grid_5d_int + module procedure read_grid_5d_int8 + end interface read_grid_5d + + interface write_grid_3d module procedure write_grid_3d_sp module procedure write_grid_3d_dp @@ -43,6 +54,14 @@ module grid_io module procedure write_grid_4d_int8 end interface write_grid_4d + interface write_grid_5d + module procedure write_grid_5d_sp + module procedure write_grid_5d_dp + module procedure write_grid_5d_int + module procedure write_grid_5d_int8 + end interface write_grid_5d + + contains logical function grid_exists(group, name) @@ -52,6 +71,28 @@ logical function grid_exists(group, name) grid_exists = mp_path_exists(group, name) end function grid_exists + !DN CRAZY ADDITIONS + subroutine read_grid_5d_int8(group, path, array, geo) + + implicit none + + integer(hid_t), intent(in) :: group + character(len=*), intent(in) :: path + integer(idp), intent(out) :: array(:,:,:) + type(grid_geometry_desc),intent(in) :: geo + + character(len=32) :: geometry_id_check + + call mp_read_keyword(group,path, 'geometry', geometry_id_check) + if(geometry_id_check.ne.geo%id) then + call error("read_grid", "geometry IDs do not match") + end if + call mp_read_array(group, path, array) + + if(any(is_nan(array))) call error("read_grid_5d", "NaN values in 5D array") + + end subroutine read_grid_5d_int8 + subroutine read_grid_4d_int8(group, path, array, geo) @@ -95,6 +136,23 @@ subroutine read_grid_3d_int8(group, path, array, geo) end subroutine read_grid_3d_int8 + + !DN CRAZY ADDITIONS + subroutine write_grid_5d_int8(group, path, array, geo) + + implicit none + + integer(hid_t), intent(in) :: group + character(len=*), intent(in) :: path + integer(idp), intent(in) :: array(:,:,:) + type(grid_geometry_desc),intent(in) :: geo + + call mp_write_array(group, path, array) + call mp_write_keyword(group, path, 'geometry', geo%id) + + end subroutine write_grid_5d_int8 + + subroutine write_grid_4d_int8(group, path, array, geo) implicit none @@ -123,6 +181,28 @@ subroutine write_grid_3d_int8(group, path, array, geo) end subroutine write_grid_3d_int8 + !DN CRAZY ADDITIONS + subroutine read_grid_5d_int(group, path, array, geo) + + implicit none + + integer(hid_t), intent(in) :: group + character(len=*), intent(in) :: path + integer, intent(out) :: array(:,:,:) + type(grid_geometry_desc),intent(in) :: geo + + character(len=32) :: geometry_id_check + + call mp_read_keyword(group,path, 'geometry', geometry_id_check) + if(geometry_id_check.ne.geo%id) then + call error("read_grid", "geometry IDs do not match") + end if + call mp_read_array(group, path, array) + + if(any(is_nan(array))) call error("read_grid_5d", "NaN values in 5D array") + + end subroutine read_grid_5d_int + subroutine read_grid_4d_int(group, path, array, geo) @@ -166,6 +246,25 @@ subroutine read_grid_3d_int(group, path, array, geo) end subroutine read_grid_3d_int + + + !DN CRAZY ADDITION + + subroutine write_grid_5d_int(group, path, array, geo) + + implicit none + + integer(hid_t), intent(in) :: group + character(len=*), intent(in) :: path + integer, intent(in) :: array(:,:,:) + type(grid_geometry_desc),intent(in) :: geo + + call mp_write_array(group, path, array) + call mp_write_keyword(group, path, 'geometry', geo%id) + + end subroutine write_grid_5d_int + + subroutine write_grid_4d_int(group, path, array, geo) implicit none @@ -194,6 +293,28 @@ subroutine write_grid_3d_int(group, path, array, geo) end subroutine write_grid_3d_int + !DN CRAZY ADDITIONS + subroutine read_grid_5d_dp(group, path, array, geo) + + implicit none + + integer(hid_t), intent(in) :: group + character(len=*), intent(in) :: path + real(dp), intent(out) :: array(:,:,:) + type(grid_geometry_desc),intent(in) :: geo + + character(len=32) :: geometry_id_check + + call mp_read_keyword(group,path, 'geometry', geometry_id_check) + if(geometry_id_check.ne.geo%id) then + call error("read_grid", "geometry IDs do not match") + end if + call mp_read_array(group, path, array) + + if(any(is_nan(array))) call error("read_grid_5d", "NaN values in 5D array") + + end subroutine read_grid_5d_dp + subroutine read_grid_4d_dp(group, path, array, geo) @@ -237,6 +358,23 @@ subroutine read_grid_3d_dp(group, path, array, geo) end subroutine read_grid_3d_dp + + !CRAZY DN ADDITIONS + subroutine write_grid_5d_dp(group, path, array, geo) + + implicit none + + integer(hid_t), intent(in) :: group + character(len=*), intent(in) :: path + real(dp), intent(in) :: array(:,:,:) + type(grid_geometry_desc),intent(in) :: geo + + call mp_write_array(group, path, array) + call mp_write_keyword(group, path, 'geometry', geo%id) + + end subroutine write_grid_5d_dp + + subroutine write_grid_4d_dp(group, path, array, geo) implicit none @@ -265,6 +403,29 @@ subroutine write_grid_3d_dp(group, path, array, geo) end subroutine write_grid_3d_dp + + !DN CRAZY ADDITIONS + subroutine read_grid_5d_sp(group, path, array, geo) + + implicit none + + integer(hid_t), intent(in) :: group + character(len=*), intent(in) :: path + real(sp), intent(out) :: array(:,:,:) + type(grid_geometry_desc),intent(in) :: geo + + character(len=32) :: geometry_id_check + + call mp_read_keyword(group,path, 'geometry', geometry_id_check) + if(geometry_id_check.ne.geo%id) then + call error("read_grid", "geometry IDs do not match") + end if + call mp_read_array(group, path, array) + + if(any(is_nan(array))) call error("read_grid_5d", "NaN values in 5D array") + + end subroutine read_grid_5d_sp + subroutine read_grid_4d_sp(group, path, array, geo) @@ -308,6 +469,23 @@ subroutine read_grid_3d_sp(group, path, array, geo) end subroutine read_grid_3d_sp + + !CRAZY DN ADDITIONS + subroutine write_grid_5d_sp(group, path, array, geo) + + implicit none + + integer(hid_t), intent(in) :: group + character(len=*), intent(in) :: path + real(sp), intent(in) :: array(:,:,:) + type(grid_geometry_desc),intent(in) :: geo + + call mp_write_array(group, path, array) + call mp_write_keyword(group, path, 'geometry', geo%id) + + end subroutine write_grid_5d_sp + + subroutine write_grid_4d_sp(group, path, array, geo) implicit none diff --git a/src/grid/grid_io_amr.f90 b/src/grid/grid_io_amr.f90 index 832e1725..6f6b1a48 100644 --- a/src/grid/grid_io_amr.f90 +++ b/src/grid/grid_io_amr.f90 @@ -14,6 +14,8 @@ module grid_io public :: read_grid_4d public :: write_grid_3d public :: write_grid_4d + public :: write_grid_5d + interface read_grid_3d module procedure read_grid_3d_sp @@ -43,6 +45,14 @@ module grid_io module procedure write_grid_4d_int8 end interface write_grid_4d + interface write_grid_5d + module procedure write_grid_5d_sp + module procedure write_grid_5d_dp + module procedure write_grid_5d_int + module procedure write_grid_5d_int8 + end interface write_grid_5d + + contains logical function grid_exists(group, name) @@ -146,6 +156,51 @@ end subroutine test end subroutine read_grid_3d_int8 + + !DN CRAZY ADDITIONS + subroutine write_grid_5d_int8(group, path, array, geo) + + implicit none + + integer(hid_t), intent(in) :: group + character(len=*), intent(in) :: path + integer(idp), intent(in) :: array(:,:,:) + type(grid_geometry_desc),intent(in),target :: geo + integer(hid_t) :: g_level, g_grid + character(len=100) :: name + integer :: ilevel, igrid + type(level_desc), pointer :: level + type(grid_desc), pointer :: grid + + do ilevel=1,size(geo%levels) + level => geo%levels(ilevel) + write(name, '("level_", I5.5)') ilevel + if(mp_path_exists(group, name)) then + g_level = mp_open_group(group, name) + else + g_level = mp_create_group(group, name) + end if + do igrid=1,size(level%grids) + grid => level%grids(igrid) + write(name, '("grid_", I5.5)') igrid + if(mp_path_exists(g_level, name)) then + g_grid = mp_open_group(g_level, name) + else + g_grid = mp_create_group(g_level, name) + end if + call mp_write_array(g_grid, path, reshape(array(grid%start_id:grid%start_id + grid%n_cells - 1, :), & + & (/grid%n1, grid%n2, grid%n3, size(array,2), size(array,3)/))) + call mp_close_group(g_grid) + end do + call mp_close_group(g_level) + + end do + + end subroutine write_grid_5d_int8 + + + + subroutine write_grid_4d_int8(group, path, array, geo) implicit none @@ -304,6 +359,47 @@ end subroutine test end subroutine read_grid_3d_int + + subroutine write_grid_5d_int(group, path, array, geo) + + implicit none + + integer(hid_t), intent(in) :: group + character(len=*), intent(in) :: path + integer, intent(in) :: array(:,:,:) + type(grid_geometry_desc),intent(in),target :: geo + integer(hid_t) :: g_level, g_grid + character(len=100) :: name + integer :: ilevel, igrid + type(level_desc), pointer :: level + type(grid_desc), pointer :: grid + + do ilevel=1,size(geo%levels) + level => geo%levels(ilevel) + write(name, '("level_", I5.5)') ilevel + if(mp_path_exists(group, name)) then + g_level = mp_open_group(group, name) + else + g_level = mp_create_group(group, name) + end if + do igrid=1,size(level%grids) + grid => level%grids(igrid) + write(name, '("grid_", I5.5)') igrid + if(mp_path_exists(g_level, name)) then + g_grid = mp_open_group(g_level, name) + else + g_grid = mp_create_group(g_level, name) + end if + call mp_write_array(g_grid, path, reshape(array(grid%start_id:grid%start_id + grid%n_cells - 1, :), & + & (/grid%n1, grid%n2, grid%n3, size(array,2), size(array,3)/))) + call mp_close_group(g_grid) + end do + call mp_close_group(g_level) + end do + + end subroutine write_grid_5d_int + + subroutine write_grid_4d_int(group, path, array, geo) implicit none @@ -462,6 +558,48 @@ end subroutine test end subroutine read_grid_3d_dp + +!DN CRAZY ADDITIONS + subroutine write_grid_5d_dp(group, path, array, geo) + + implicit none + + integer(hid_t), intent(in) :: group + character(len=*), intent(in) :: path + real(dp), intent(in) :: array(:,:,:) + type(grid_geometry_desc),intent(in),target :: geo + integer(hid_t) :: g_level, g_grid + character(len=100) :: name + integer :: ilevel, igrid + type(level_desc), pointer :: level + type(grid_desc), pointer :: grid + + do ilevel=1,size(geo%levels) + level => geo%levels(ilevel) + write(name, '("level_", I5.5)') ilevel + if(mp_path_exists(group, name)) then + g_level = mp_open_group(group, name) + else + g_level = mp_create_group(group, name) + end if + do igrid=1,size(level%grids) + grid => level%grids(igrid) + write(name, '("grid_", I5.5)') igrid + if(mp_path_exists(g_level, name)) then + g_grid = mp_open_group(g_level, name) + else + g_grid = mp_create_group(g_level, name) + end if + call mp_write_array(g_grid, path, reshape(array(grid%start_id:grid%start_id + grid%n_cells - 1, :), & + & (/grid%n1, grid%n2, grid%n3, size(array,2), size(array,3)/))) + call mp_close_group(g_grid) + end do + call mp_close_group(g_level) + end do + + end subroutine write_grid_5d_dp + + subroutine write_grid_4d_dp(group, path, array, geo) implicit none @@ -620,6 +758,49 @@ end subroutine test end subroutine read_grid_3d_sp + +!DN CRAZY ADDITIONS + subroutine write_grid_5d_sp(group, path, array, geo) + + implicit none + + integer(hid_t), intent(in) :: group + character(len=*), intent(in) :: path + real(sp), intent(in) :: array(:,:,:) + type(grid_geometry_desc),intent(in),target :: geo + integer(hid_t) :: g_level, g_grid + character(len=100) :: name + integer :: ilevel, igrid + type(level_desc), pointer :: level + type(grid_desc), pointer :: grid + + do ilevel=1,size(geo%levels) + level => geo%levels(ilevel) + write(name, '("level_", I5.5)') ilevel + if(mp_path_exists(group, name)) then + g_level = mp_open_group(group, name) + else + g_level = mp_create_group(group, name) + end if + do igrid=1,size(level%grids) + grid => level%grids(igrid) + write(name, '("grid_", I5.5)') igrid + if(mp_path_exists(g_level, name)) then + g_grid = mp_open_group(g_level, name) + else + g_grid = mp_create_group(g_level, name) + end if + call mp_write_array(g_grid, path, reshape(array(grid%start_id:grid%start_id + grid%n_cells - 1, :), & + & (/grid%n1, grid%n2, grid%n3, size(array,2),size(array,3)/))) + call mp_close_group(g_grid) + end do + call mp_close_group(g_level) + end do + + end subroutine write_grid_5d_sp + + + subroutine write_grid_4d_sp(group, path, array, geo) implicit none diff --git a/src/grid/grid_io_amr_template.f90 b/src/grid/grid_io_amr_template.f90 index f80eb96c..e3abfedc 100644 --- a/src/grid/grid_io_amr_template.f90 +++ b/src/grid/grid_io_amr_template.f90 @@ -13,6 +13,8 @@ module grid_io public :: read_grid_4d public :: write_grid_3d public :: write_grid_4d + public :: write_grid_5d + interface read_grid_3d module procedure read_grid_3d_sp @@ -42,6 +44,15 @@ module grid_io module procedure write_grid_4d_int8 end interface write_grid_4d + interface write_grid_5d + module procedure write_grid_5d_sp + module procedure write_grid_5d_dp + module procedure write_grid_5d_int + module procedure write_grid_5d_int8 + end interface write_grid_5d + + + contains logical function grid_exists(group, name) @@ -146,6 +157,48 @@ end subroutine test end subroutine read_grid_3d_ + + !DN CRAZY ADDITIONS + subroutine write_grid_5d_(group, path, array, geo) + + implicit none + + integer(hid_t), intent(in) :: group + character(len=*), intent(in) :: path + @T, intent(in) :: array(:,:,:) + type(grid_geometry_desc),intent(in),target :: geo + integer(hid_t) :: g_level, g_grid + character(len=100) :: name + integer :: ilevel, igrid + type(level_desc), pointer :: level + type(grid_desc), pointer :: grid + + do ilevel=1,size(geo%levels) + level => geo%levels(ilevel) + write(name, '("level_", I5.5)') ilevel + if(mp_path_exists(group, name)) then + g_level = mp_open_group(group, name) + else + g_level = mp_create_group(group, name) + end if + do igrid=1,size(level%grids) + grid => level%grids(igrid) + write(name, '("grid_", I5.5)') igrid + if(mp_path_exists(g_level, name)) then + g_grid = mp_open_group(g_level, name) + else + g_grid = mp_create_group(g_level, name) + end if + call mp_write_array(g_grid, path, reshape(array(grid%start_id:grid%start_id + grid%n_cells - 1, :), & + & (/grid%n1, grid%n2, grid%n3, size(array,2), size(array,3)/))) + call mp_close_group(g_grid) + end do + call mp_close_group(g_level) + + end do + + end subroutine write_grid_5d_ + subroutine write_grid_4d_(group, path, array, geo) implicit none diff --git a/src/grid/grid_io_template.f90 b/src/grid/grid_io_template.f90 index f42e6eab..439ca0d4 100644 --- a/src/grid/grid_io_template.f90 +++ b/src/grid/grid_io_template.f90 @@ -11,8 +11,11 @@ module grid_io public :: grid_exists public :: read_grid_3d public :: read_grid_4d + public :: read_grid_5d public :: write_grid_3d public :: write_grid_4d + public :: write_grid_5d + interface read_grid_3d module procedure read_grid_3d_sp @@ -28,6 +31,13 @@ module grid_io module procedure read_grid_4d_int8 end interface read_grid_4d + interface read_grid_5d + module procedure read_grid_5d_sp + module procedure read_grid_5d_dp + module procedure read_grid_5d_int + module procedure read_grid_5d_int8 + end interface read_grid_5d + interface write_grid_3d module procedure write_grid_3d_sp module procedure write_grid_3d_dp @@ -42,6 +52,14 @@ module grid_io module procedure write_grid_4d_int8 end interface write_grid_4d + interface write_grid_5d + module procedure write_grid_5d_sp + module procedure write_grid_5d_dp + module procedure write_grid_5d_int + module procedure write_grid_5d_int8 + end interface write_grid_5d + + contains logical function grid_exists(group, name) @@ -53,6 +71,38 @@ end function grid_exists !!@FOR real(sp):sp real(dp):dp integer:int integer(idp):int8 + + subroutine read_grid_5d_(group, path, array, geo) + + implicit none + + integer(hid_t), intent(in) :: group + character(len=*), intent(in) :: path + @T, intent(out) :: array(:,:,:) + type(grid_geometry_desc),intent(in) :: geo + @T, allocatable :: array5d(:,:,:,:,:) + integer :: n_cells, n_dust, n_isrf_lam + + character(len=32) :: geometry_id_check + + call mp_read_keyword(group,path, 'geometry', geometry_id_check) + if(geometry_id_check.ne.geo%id) then + call error("read_grid", "geometry IDs do not match") + end if + call mp_read_array_auto(group,path, array5d) + + if(any(is_nan(array5d))) call error("read_grid_5d", "NaN values in 5D array") + + n_cells = size(array, 1) + n_dust = size(array, 2) + n_isrf_lam = size(array,3) + + array = reshape(array5d, (/n_cells, n_dust, n_isrf_lam/)) + + end subroutine read_grid_5d_ + + + subroutine read_grid_4d_(group, path, array, geo) implicit none @@ -107,6 +157,22 @@ subroutine read_grid_3d_(group, path, array, geo) array = reshape(array3d, (/n_cells/)) end subroutine read_grid_3d_ + + subroutine write_grid_5d_(group, path, array, geo) + + implicit none + + integer(hid_t), intent(in) :: group + character(len=*), intent(in) :: path + @T, intent(in) :: array(:,:,:) + type(grid_geometry_desc),intent(in) :: geo + + call mp_write_array(group, path, reshape(array, (/geo%n1, geo%n2, geo%n3, size(array,2), size(array,3)/))) + call mp_write_keyword(group, path, 'geometry', geo%id) + + end subroutine write_grid_5d_ + + subroutine write_grid_4d_(group, path, array, geo) diff --git a/src/grid/grid_physics_3d.f90 b/src/grid/grid_physics_3d.f90 index 429f041d..6665f46e 100644 --- a/src/grid/grid_physics_3d.f90 +++ b/src/grid/grid_physics_3d.f90 @@ -195,7 +195,7 @@ subroutine setup_grid_physics(group, use_mrw, use_pda) allocate(specific_energy_sum_nu(geo%n_cells, n_dust, 10)) specific_energy_sum_nu = 0._dp - + ! Total energy absorbed allocate(energy_abs_tot(n_dust)) energy_abs_tot = 0._dp diff --git a/src/grid/grid_propagate_3d.f90 b/src/grid/grid_propagate_3d.f90 index 4a04f6f9..b5d3bc9f 100644 --- a/src/grid/grid_propagate_3d.f90 +++ b/src/grid/grid_propagate_3d.f90 @@ -10,6 +10,9 @@ module grid_propagate use counters use settings, only : frac_check => propagation_check_frequency + !DN CRAZY ADDITIONS + use grid_geometry, only : geo + implicit none save @@ -57,6 +60,23 @@ subroutine grid_integrate(p,tau_required,tau_achieved) integer :: source_id + + !DN CRAZY ADDITIONS + integer :: idx + real, dimension(10) :: energy_frequency_bins + energy_frequency_bins(1) = 10.**14.4768207 + energy_frequency_bins(2) = 10.**14.59237683 + energy_frequency_bins(3) = 10.**14.70793296 + energy_frequency_bins(4) = 10.**14.82348909 + energy_frequency_bins(5) = 10.**14.93904522 + energy_frequency_bins(6) = 10.**15.05460135 + energy_frequency_bins(7) = 10.**15.17015748 + energy_frequency_bins(8) = 10.**15.28571361 + energy_frequency_bins(9) = 10.**15.40126974 + energy_frequency_bins(10) = 10.**15.51682586 + + + radial = (p%r .dot. p%v) > 0. if(debug) write(*,'(" [debug] start grid_integrate")') @@ -140,21 +160,27 @@ subroutine grid_integrate(p,tau_required,tau_achieved) !print *,' ',minloc(abs(energy_frequency_bins-p%nu)) !DN CRAZY ADDITIONS - !idx = minloc(abs(energy_frequency_bins-p%nu),DIM=1) + idx = minloc(abs(energy_frequency_bins-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) = & & specific_energy_sum(p%icell%ic, id) + tmin * p%current_kappa(id) * p%energy end if - - !if(density(p%icell%ic,id) > 0._dp) then + + + !DN CRAZY ADDITIONS + if(density(p%icell%ic,id) > 0._dp) then + specific_energy_sum_nu(p%icell%ic,id,idx) = & + & specific_energy_sum_nu(p%icell%ic,id,idx) + tmin * p%current_kappa(id) * p%energy + end if - !specific_energy_sum_nu(p%icell%ic,id,1) = & - ! & specific_energy_sum_nu(p%icell%ic,id,idx) + tmin * p%current_kappa(id) * p%energy - !end if end do p%on_wall = .true. @@ -309,6 +335,9 @@ subroutine grid_integrate_noenergy(p,tau_required,tau_achieved) !idx = minloc(abs(energy_frequency_bins-p%nu),DIM=1) !print *,'[grid_propagate_3d last iteration] p%energy=',p%energy !print *,'[grid_propagate_3d last iteration] p%nu=',p%nu + !print *, "shape(specific_energy_sum_nu", shape(specific_energy_sum_nu) + + idx = minloc(abs(energy_frequency_bins-p%nu),DIM=1) do id=1,n_dust @@ -316,9 +345,9 @@ subroutine grid_integrate_noenergy(p,tau_required,tau_achieved) !DN CRAZY ADDITIONS if(density(p%icell%ic,id) > 0._dp) then - specific_energy_sum_nu(p%icell%ic,id,1) = & + specific_energy_sum_nu(p%icell%ic,id,idx) = & & specific_energy_sum_nu(p%icell%ic,id,idx) + tmin * p%current_kappa(id) * p%energy - !print *,'[grid_propage_3d last iteration] specific_energy_sum_nu=',specific_energy_sum_nu + !print *,'[grid_propage_3d last iteration] specific_energy_sum_nu=',specific_energy_sum_nu(p%icell%ic,id,idx) end if end do From 4ca14bd2c2cf3b600c5d623c77e9b5bab7248741 Mon Sep 17 00:00:00 2001 From: dnarayanan Date: Wed, 24 Nov 2021 15:54:52 -0500 Subject: [PATCH 04/38] specific_energy_sum_nu saves and reads. AMR still needs to be dealt with --- hyperion/conf/conf_files.py | 2 +- hyperion/grid/grid_helpers.py | 12 +++++-- src/grid/grid_io_1d_template.f90 | 56 ++++++++++++++++++++++++++++++++ src/grid/grid_physics_3d.f90 | 2 ++ 4 files changed, 69 insertions(+), 3 deletions(-) diff --git a/hyperion/conf/conf_files.py b/hyperion/conf/conf_files.py index dac010d7..6e66ab0f 100644 --- a/hyperion/conf/conf_files.py +++ b/hyperion/conf/conf_files.py @@ -36,7 +36,7 @@ def write(self, group): group.attrs['output_density'] = np.bytes_(self.output_density.encode('utf-8')) group.attrs['output_density_diff'] = np.bytes_(self.output_density_diff.encode('utf-8')) group.attrs['output_specific_energy'] = np.bytes_(self.output_specific_energy.encode('utf-8')) - group.attrs['output_specific_energy_nu'] = np.bytes_(self.output_specific_energ_nu.encode('utf-8')) + group.attrs['output_specific_energy_nu'] = np.bytes_(self.output_specific_energy_nu.encode('utf-8')) group.attrs['output_n_photons'] = np.bytes_(self.output_n_photons.encode('utf-8')) diff --git a/hyperion/grid/grid_helpers.py b/hyperion/grid/grid_helpers.py index cadb3496..55194e84 100644 --- a/hyperion/grid/grid_helpers.py +++ b/hyperion/grid/grid_helpers.py @@ -22,7 +22,7 @@ def single_grid_dims(data, ndim=3): ''' if type(data) in [list, tuple]: - + n_pop = len(data) shape = None for item in data: @@ -34,14 +34,22 @@ def single_grid_dims(data, ndim=3): if shape is not None and len(shape) != ndim: raise ValueError("Grids should be %i-dimensional" % ndim) - elif isinstance(data, np.ndarray): + + elif isinstance(data, np.ndarray): if data.ndim == ndim: n_pop = None shape = data.shape elif data.ndim == ndim + 1: n_pop = data.shape[0] shape = data[0].shape + + #DN CRAZY ADDITIONS + elif data.ndim == ndim+2: + n_pop = data.shape[1] + shape = data.shape[-1] + shape = tuple(map(int,str(shape).split(' '))) + else: raise Exception("Unexpected number of dimensions: %i" % data.ndim) diff --git a/src/grid/grid_io_1d_template.f90 b/src/grid/grid_io_1d_template.f90 index 1abbe776..f34ceec1 100644 --- a/src/grid/grid_io_1d_template.f90 +++ b/src/grid/grid_io_1d_template.f90 @@ -11,6 +11,7 @@ module grid_io public :: grid_exists public :: read_grid_3d public :: read_grid_4d + public :: read_grid_5d public :: write_grid_3d public :: write_grid_4d @@ -28,6 +29,13 @@ module grid_io module procedure read_grid_4d_int8 end interface read_grid_4d + interface read_grid_5d + module procedure read_grid_5d_sp + module procedure read_grid_5d_dp + module procedure read_grid_5d_int + module procedure read_grid_5d_int8 + end interface read_grid_5d + interface write_grid_3d module procedure write_grid_3d_sp module procedure write_grid_3d_dp @@ -41,6 +49,15 @@ module grid_io module procedure write_grid_4d_int module procedure write_grid_4d_int8 end interface write_grid_4d + + interface write_grid_5d + module procedure write_grid_5d_sp + module procedure write_grid_5d_dp + module procedure write_grid_5d_int + module procedure write_grid_5d_int8 + end interface write_grid_5d + + contains @@ -53,6 +70,29 @@ end function grid_exists !!@FOR real(sp):sp real(dp):dp integer:int integer(idp):int8 + + subroutine read_grid_5d_(group, path, array, geo) + + implicit none + + integer(hid_t), intent(in) :: group + character(len=*), intent(in) :: path + @T, intent(out) :: array(:,:,:) + type(grid_geometry_desc),intent(in) :: geo + + character(len=32) :: geometry_id_check + + call mp_read_keyword(group,path, 'geometry', geometry_id_check) + if(geometry_id_check.ne.geo%id) then + call error("read_grid", "geometry IDs do not match") + end if + call mp_read_array(group, path, array) + + if(any(is_nan(array))) call error("read_grid_5d", "NaN values in 5D array") + + end subroutine read_grid_5d_ + + subroutine read_grid_4d_(group, path, array, geo) implicit none @@ -95,6 +135,22 @@ subroutine read_grid_3d_(group, path, array, geo) end subroutine read_grid_3d_ + + subroutine write_grid_5d_(group, path, array, geo) + + implicit none + + integer(hid_t), intent(in) :: group + character(len=*), intent(in) :: path + @T, intent(in) :: array(:,:,:) + type(grid_geometry_desc),intent(in) :: geo + + call mp_write_array(group, path, array) + call mp_write_keyword(group, path, 'geometry', geo%id) + + end subroutine write_grid_5d_ + + subroutine write_grid_4d_(group, path, array, geo) implicit none diff --git a/src/grid/grid_physics_3d.f90 b/src/grid/grid_physics_3d.f90 index 6665f46e..2c1d600e 100644 --- a/src/grid/grid_physics_3d.f90 +++ b/src/grid/grid_physics_3d.f90 @@ -193,6 +193,8 @@ subroutine setup_grid_physics(group, use_mrw, use_pda) allocate(specific_energy_sum(geo%n_cells, n_dust)) specific_energy_sum = 0._dp + + !DN CRAZY ADDITION allocate(specific_energy_sum_nu(geo%n_cells, n_dust, 10)) specific_energy_sum_nu = 0._dp From 4f59a2b620991cd9651e07be487a6e1d6b820b06 Mon Sep 17 00:00:00 2001 From: dnarayanan Date: Thu, 2 Dec 2021 08:29:06 -0500 Subject: [PATCH 05/38] grid_io_amr_template updated but commetned out since i cant quite get it to compile --- src/grid/grid_io_amr.f90 | 188 +++++++++++++++++++++++++++++++++++++-- 1 file changed, 183 insertions(+), 5 deletions(-) diff --git a/src/grid/grid_io_amr.f90 b/src/grid/grid_io_amr.f90 index 6f6b1a48..0258db78 100644 --- a/src/grid/grid_io_amr.f90 +++ b/src/grid/grid_io_amr.f90 @@ -31,6 +31,21 @@ module grid_io module procedure read_grid_4d_int8 end interface read_grid_4d +! The 5d grid reading is commented out currently becuase the compilation crashes in read_grid_5d, owing to some issue in the where statement. we get an error: + +!src/grid/grid_io_amr.f90(113): error #6783: The variable being defined does not conform with the mask-expr of the where-construct or where-stmt. [ARRAY5D] +! array5d(:,:,:,:,:) = 0 + +!i think this could be due to something that needs to be updated in type_grid_amr.f90 + +! interface read_grid_5d +! module procedure read_grid_5d_sp +! module procedure read_grid_5d_dp +! module procedure read_grid_5d_int +! module procedure read_grid_5d_int8 +! end interface read_grid_5d + + interface write_grid_3d module procedure write_grid_3d_sp module procedure write_grid_3d_dp @@ -78,6 +93,47 @@ logical function grid_exists(group, name) end function grid_exists +! subroutine read_grid_5d_int8(group, path, array, geo) + +! implicit none + +! integer(hid_t), intent(in) :: group +! character(len=*), intent(in) :: path +! integer(idp), intent(out) :: array(:,:,:) +! type(grid_geometry_desc),intent(in),target :: geo +! integer(idp), allocatable :: array5d(:,:,:,:,:) +! character(len=100) :: full_path +! integer :: ilevel, igrid, idust +! type(level_desc), pointer :: level +! type(grid_desc), pointer :: grid + +! do ilevel=1,size(geo%levels) +! level => geo%levels(ilevel) +! do igrid=1,size(level%grids) +! grid => level%grids(igrid) +! write(full_path, '("level_", I5.5, "/grid_", I5.5,"/")') ilevel, igrid +! full_path = trim(full_path)//trim(path) +! call mp_read_array_auto(group, full_path, array5d) +! if(any(is_nan(array5d))) call error("read_grid_5d", "NaN values in 5D array") +! do idust=1,size(array5d, 5) +! where(grid%goto_grid(1:grid%n1,1:grid%n2,1:grid%n3) > 0) +! array5d(:,:,:,:,idust) = 0 +! end where +! end do +! array(grid%start_id:grid%start_id + grid%n_cells - 1, :, :) = reshape(array5d, (/grid%n_cells, size(array, 2), size(array,3)/)) +! end do +! end do + +! ! The following three lines provide a workaround for the PGI Fortran +! ! compiler, which otherwise crashes with the following error: +! ! 0: RESHAPE: result type != SOURCE type +! contains +! subroutine test() +! end subroutine test + +! end subroutine read_grid_5d_int8 + + subroutine read_grid_4d_int8(group, path, array, geo) implicit none @@ -188,7 +244,8 @@ subroutine write_grid_5d_int8(group, path, array, geo) else g_grid = mp_create_group(g_level, name) end if - call mp_write_array(g_grid, path, reshape(array(grid%start_id:grid%start_id + grid%n_cells - 1, :), & + + call mp_write_array(g_grid, path, reshape(array(grid%start_id:grid%start_id + grid%n_cells - 1, :,:), & & (/grid%n1, grid%n2, grid%n3, size(array,2), size(array,3)/))) call mp_close_group(g_grid) end do @@ -281,6 +338,48 @@ subroutine write_grid_3d_int8(group, path, array, geo) end subroutine write_grid_3d_int8 + +! subroutine read_grid_5d_int(group, path, array, geo) + +! implicit none + +! integer(hid_t), intent(in) :: group +! character(len=*), intent(in) :: path +! integer, intent(out) :: array(:,:,:) +! type(grid_geometry_desc),intent(in),target :: geo +! integer, allocatable :: array5d(:,:,:,:,:) +! character(len=100) :: full_path +! integer :: ilevel, igrid, idust +! type(level_desc), pointer :: level +! type(grid_desc), pointer :: grid + +! do ilevel=1,size(geo%levels) +! level => geo%levels(ilevel) +! do igrid=1,size(level%grids) +! grid => level%grids(igrid) +! write(full_path, '("level_", I5.5, "/grid_", I5.5,"/")') ilevel, igrid +! full_path = trim(full_path)//trim(path) +! call mp_read_array_auto(group, full_path, array5d) +! if(any(is_nan(array5d))) call error("read_grid_5d", "NaN values in 5D array") +! do idust=1,size(array5d, 5) +! where(grid%goto_grid(1:grid%n1,1:grid%n2,1:grid%n3) > 0) +! array5d(:,:,:,:,idust) = 0 +! end where +! end do +! array(grid%start_id:grid%start_id + grid%n_cells - 1, :, :) = reshape(array5d, (/grid%n_cells, size(array, 2), size(array,3)/)) +! end do +! end do + +! ! The following three lines provide a workaround for the PGI Fortran +! ! compiler, which otherwise crashes with the following error: +! ! 0: RESHAPE: result type != SOURCE type +! contains +! subroutine test() +! end subroutine test + +! end subroutine read_grid_5d_int + + subroutine read_grid_4d_int(group, path, array, geo) implicit none @@ -359,7 +458,7 @@ end subroutine test end subroutine read_grid_3d_int - + !DN CRAZY ADDITIONS subroutine write_grid_5d_int(group, path, array, geo) implicit none @@ -390,7 +489,7 @@ subroutine write_grid_5d_int(group, path, array, geo) else g_grid = mp_create_group(g_level, name) end if - call mp_write_array(g_grid, path, reshape(array(grid%start_id:grid%start_id + grid%n_cells - 1, :), & + call mp_write_array(g_grid, path, reshape(array(grid%start_id:grid%start_id + grid%n_cells - 1, :,:), & & (/grid%n1, grid%n2, grid%n3, size(array,2), size(array,3)/))) call mp_close_group(g_grid) end do @@ -480,6 +579,46 @@ subroutine write_grid_3d_int(group, path, array, geo) end subroutine write_grid_3d_int +! subroutine read_grid_5d_dp(group, path, array, geo) + +! implicit none + +! integer(hid_t), intent(in) :: group +! character(len=*), intent(in) :: path +! real(dp), intent(out) :: array(:,:,:) +! type(grid_geometry_desc),intent(in),target :: geo +! real(dp), allocatable :: array5d(:,:,:,:,:) +! character(len=100) :: full_path +! integer :: ilevel, igrid, idust +! type(level_desc), pointer :: level +! type(grid_desc), pointer :: grid + +! do ilevel=1,size(geo%levels) +! level => geo%levels(ilevel) +! do igrid=1,size(level%grids) +! grid => level%grids(igrid) +! write(full_path, '("level_", I5.5, "/grid_", I5.5,"/")') ilevel, igrid +! full_path = trim(full_path)//trim(path) +! call mp_read_array_auto(group, full_path, array5d) +! if(any(is_nan(array5d))) call error("read_grid_5d", "NaN values in 5D array") +! do idust=1,size(array5d, 5) +! where(grid%goto_grid(1:grid%n1,1:grid%n2,1:grid%n3) > 0) +! array5d(:,:,:,:,idust) = 0 +! end where +! end do +! array(grid%start_id:grid%start_id + grid%n_cells - 1, :, :) = reshape(array5d, (/grid%n_cells, size(array, 2), size(array,3)/)) +! end do +! end do + +! ! The following three lines provide a workaround for the PGI Fortran +! ! compiler, which otherwise crashes with the following error: +! ! 0: RESHAPE: result type != SOURCE type +! contains +! subroutine test() +! end subroutine test + +! end subroutine read_grid_5d_dp + subroutine read_grid_4d_dp(group, path, array, geo) implicit none @@ -590,7 +729,7 @@ subroutine write_grid_5d_dp(group, path, array, geo) else g_grid = mp_create_group(g_level, name) end if - call mp_write_array(g_grid, path, reshape(array(grid%start_id:grid%start_id + grid%n_cells - 1, :), & + call mp_write_array(g_grid, path, reshape(array(grid%start_id:grid%start_id + grid%n_cells - 1, :, :), & & (/grid%n1, grid%n2, grid%n3, size(array,2), size(array,3)/))) call mp_close_group(g_grid) end do @@ -679,6 +818,45 @@ subroutine write_grid_3d_dp(group, path, array, geo) end subroutine write_grid_3d_dp +! subroutine read_grid_5d_sp(group, path, array, geo) + +! implicit none + +! integer(hid_t), intent(in) :: group +! character(len=*), intent(in) :: path +! real(sp), intent(out) :: array(:,:,:) +! type(grid_geometry_desc),intent(in),target :: geo +! real(sp), allocatable :: array5d(:,:,:,:,:) +! character(len=100) :: full_path +! integer :: ilevel, igrid, idust +! type(level_desc), pointer :: level +! type(grid_desc), pointer :: grid + +! do ilevel=1,size(geo%levels) +! level => geo%levels(ilevel) +! do igrid=1,size(level%grids) +! grid => level%grids(igrid) +! write(full_path, '("level_", I5.5, "/grid_", I5.5,"/")') ilevel, igrid +! full_path = trim(full_path)//trim(path) +! call mp_read_array_auto(group, full_path, array5d) +! if(any(is_nan(array5d))) call error("read_grid_5d", "NaN values in 5D array") +! do idust=1,size(array5d, 5) +! where(grid%goto_grid(1:grid%n1,1:grid%n2,1:grid%n3) > 0) +! array5d(:,:,:,:,idust) = 0 +! end where +! end do +! array(grid%start_id:grid%start_id + grid%n_cells - 1, :,:) = reshape(array5d, (/grid%n_cells, size(array, 2),size(array,3)/)) +! end do +! end do + +! ! The following three lines provide a workaround for the PGI Fortran +! ! compiler, which otherwise crashes with the following error: +! ! 0: RESHAPE: result type != SOURCE type +! contains +! subroutine test() +! end subroutine test + +! end subroutine read_grid_5d_sp subroutine read_grid_4d_sp(group, path, array, geo) @@ -790,7 +968,7 @@ subroutine write_grid_5d_sp(group, path, array, geo) else g_grid = mp_create_group(g_level, name) end if - call mp_write_array(g_grid, path, reshape(array(grid%start_id:grid%start_id + grid%n_cells - 1, :), & + call mp_write_array(g_grid, path, reshape(array(grid%start_id:grid%start_id + grid%n_cells - 1, :, :), & & (/grid%n1, grid%n2, grid%n3, size(array,2),size(array,3)/))) call mp_close_group(g_grid) end do From 57a1914fceb35a73a9a9465b508bc4f9ccefcf73 Mon Sep 17 00:00:00 2001 From: dnarayanan Date: Thu, 9 Dec 2021 15:32:13 -0500 Subject: [PATCH 06/38] grid modules now takes the n_isrf_wavelengths from the dust files --- src/grid/grid_generic.f90 | 14 ++++++------ src/grid/grid_physics_3d.f90 | 6 ++++- src/grid/grid_propagate_3d.f90 | 41 ++++++++++++---------------------- 3 files changed, 26 insertions(+), 35 deletions(-) diff --git a/src/grid/grid_generic.f90 b/src/grid/grid_generic.f90 index 14397b6a..db6c6794 100644 --- a/src/grid/grid_generic.f90 +++ b/src/grid/grid_generic.f90 @@ -71,14 +71,14 @@ subroutine output_grid(group, iter, n_iter) n_cells = size(specific_energy_sum_nu, 1) n_dust = size(specific_energy_sum_nu, 2) n_isrf_lam = size(specific_energy_sum_nu,3) - do i=1,n_cells - do j=1,n_dust + !do i=1,n_cells + ! do j=1,n_dust !print *, "specific_energy_sum = ",specific_energy_sum(i,j) - do k=1,n_isrf_lam - print *, "specific_energy_sum_nu = ",specific_energy_sum_nu(i,j,k) - end do - end do - end do + !do k=1,n_isrf_lam + !print *, "specific_energy_sum_nu = ",specific_energy_sum_nu(i,j,k) + ! end do + ! end do + !end do diff --git a/src/grid/grid_physics_3d.f90 b/src/grid/grid_physics_3d.f90 index 2c1d600e..cfc3179b 100644 --- a/src/grid/grid_physics_3d.f90 +++ b/src/grid/grid_physics_3d.f90 @@ -59,6 +59,8 @@ module grid_physics type(pdf_discrete_dp) :: absorption + + contains real(dp) function tau_inv_planck_to_closest_wall(p) result(tau) @@ -97,6 +99,7 @@ subroutine setup_grid_physics(group, use_mrw, use_pda) integer(hid_t),intent(in) :: group logical,intent(in) :: use_mrw, use_pda + integer :: n_isrf_wavelengths ! Density allocate(density(geo%n_cells, n_dust)) @@ -195,7 +198,8 @@ subroutine setup_grid_physics(group, use_mrw, use_pda) !DN CRAZY ADDITION - allocate(specific_energy_sum_nu(geo%n_cells, n_dust, 10)) + n_isrf_wavelengths = d(1)%n_nu + allocate(specific_energy_sum_nu(geo%n_cells, n_dust, n_isrf_wavelengths)) specific_energy_sum_nu = 0._dp ! Total energy absorbed diff --git a/src/grid/grid_propagate_3d.f90 b/src/grid/grid_propagate_3d.f90 index b5d3bc9f..428f419f 100644 --- a/src/grid/grid_propagate_3d.f90 +++ b/src/grid/grid_propagate_3d.f90 @@ -3,7 +3,7 @@ module grid_propagate use core_lib use type_photon, only : photon use type_grid_cell - use dust_main, only : n_dust + use dust_main, only : n_dust,d use grid_geometry, only : escaped, find_wall, in_correct_cell, next_cell, opposite_wall use grid_physics, only : specific_energy_sum, specific_energy_sum_nu, density, n_photons, last_photon_id use sources @@ -63,19 +63,14 @@ subroutine grid_integrate(p,tau_required,tau_achieved) !DN CRAZY ADDITIONS integer :: idx - real, dimension(10) :: energy_frequency_bins - energy_frequency_bins(1) = 10.**14.4768207 - energy_frequency_bins(2) = 10.**14.59237683 - energy_frequency_bins(3) = 10.**14.70793296 - energy_frequency_bins(4) = 10.**14.82348909 - energy_frequency_bins(5) = 10.**14.93904522 - energy_frequency_bins(6) = 10.**15.05460135 - energy_frequency_bins(7) = 10.**15.17015748 - energy_frequency_bins(8) = 10.**15.28571361 - energy_frequency_bins(9) = 10.**15.40126974 - energy_frequency_bins(10) = 10.**15.51682586 + real, dimension(d(1)%n_nu) :: energy_frequency_bins + + do id=1,d(1)%n_nu + energy_frequency_bins(id) = d(1)%nu(id) + end do + radial = (p%r .dot. p%v) > 0. @@ -163,10 +158,8 @@ subroutine grid_integrate(p,tau_required,tau_achieved) idx = minloc(abs(energy_frequency_bins-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) = & & specific_energy_sum(p%icell%ic, id) + tmin * p%current_kappa(id) * p%energy @@ -266,20 +259,14 @@ subroutine grid_integrate_noenergy(p,tau_required,tau_achieved) integer :: source_id !DN CRAZY ADDITIONS + integer :: id integer :: idx - !DN CRAZY ADDITIONS - real, dimension(10) :: energy_frequency_bins - energy_frequency_bins(1) = 10.**14.4768207 - energy_frequency_bins(2) = 10.**14.59237683 - energy_frequency_bins(3) = 10.**14.70793296 - energy_frequency_bins(4) = 10.**14.82348909 - energy_frequency_bins(5) = 10.**14.93904522 - energy_frequency_bins(6) = 10.**15.05460135 - energy_frequency_bins(7) = 10.**15.17015748 - energy_frequency_bins(8) = 10.**15.28571361 - energy_frequency_bins(9) = 10.**15.40126974 - energy_frequency_bins(10) = 10.**15.51682586 + real, dimension(d(0)%n_nu) :: energy_frequency_bins + do id=1,d(0)%n_nu + energy_frequency_bins(id) = d(0)%nu(id) + print *,'[grid_propagate_3d last iteration] energy_frequency_bins(id) = ',energy_frequency_bins(id) + end do radial = (p%r .dot. p%v) > 0. From 964016eca9388a5bf3915f6130194742ab4d4748 Mon Sep 17 00:00:00 2001 From: dnarayanan Date: Mon, 13 Dec 2021 13:09:07 -0500 Subject: [PATCH 07/38] outputting the ISRF frequency bins --- src/grid/grid_generic.f90 | 29 ++++++++++++++++++----------- 1 file changed, 18 insertions(+), 11 deletions(-) diff --git a/src/grid/grid_generic.f90 b/src/grid/grid_generic.f90 index db6c6794..c80fca7c 100644 --- a/src/grid/grid_generic.f90 +++ b/src/grid/grid_generic.f90 @@ -3,12 +3,16 @@ module grid_generic use core_lib, only : sp, dp, hid_t, warn, error use mpi_hdf5_io, only : mp_create_group use mpi_core, only : main_process - + use grid_io, only : write_grid_3d, write_grid_4d, write_grid_5d use grid_geometry, only : geo use grid_physics, only : n_photons, last_photon_id, specific_energy_sum, specific_energy_sum_nu, specific_energy, density, density_original use settings, only : output_n_photons, output_specific_energy, output_density, output_density_diff, physics_io_type + !DN Crazy Additions + use dust_main, only: d + + implicit none save @@ -35,6 +39,7 @@ subroutine output_grid(group, iter, n_iter) !DN CRAZY ADDITION integer :: n_cells, n_dust, n_isrf_lam integer :: i,j,k + real, dimension(d(1)%n_nu) :: energy_frequency_bins if(main_process()) write(*,'(" [output_grid] outputting grid arrays for iteration")') @@ -68,24 +73,26 @@ subroutine output_grid(group, iter, n_iter) ! DN Crazy Additions + ! WRITE THE ISRF n_cells = size(specific_energy_sum_nu, 1) n_dust = size(specific_energy_sum_nu, 2) n_isrf_lam = size(specific_energy_sum_nu,3) - !do i=1,n_cells - ! do j=1,n_dust - !print *, "specific_energy_sum = ",specific_energy_sum(i,j) - !do k=1,n_isrf_lam - !print *, "specific_energy_sum_nu = ",specific_energy_sum_nu(i,j,k) - ! end do - ! end do - !end do + do i=1,d(1)%n_nu + energy_frequency_bins(i) = d(1)%nu(i) + end do + if(trim(output_specific_energy)=='all' .or. (trim(output_specific_energy)=='last'.and.iter==n_iter)) then + !if(allocated(energy_frequency_bins)) then + call write_grid_3d(group, 'ISRF_frequency_bins',energy_frequency_bins, geo) + else + call warn("output_grid","energy_frequency bins [ISRF wavelengths] array is not allocated") + !end if + end if if(trim(output_specific_energy)=='all' .or. (trim(output_specific_energy)=='last'.and.iter==n_iter)) then if(allocated(specific_energy_sum_nu)) then - !print *, "shape(specific_energy_sum_nu", shape(specific_energy_sum_nu) - !print *, "specific_energy_sum_nu = ",specific_energy_sum_nu + select case(physics_io_type) case(sp) From e6928a08c097f24c6bf82c1074399cf86e9bc8f2 Mon Sep 17 00:00:00 2001 From: dnarayanan Date: Tue, 14 Dec 2021 11:09:55 -0500 Subject: [PATCH 08/38] can now call whether ISRF is computed and saved via m.compute_isrf(True) --- hyperion/conf/conf_files.py | 26 ++++++++++++ hyperion/conf/tests/test_conf_io.py | 10 +++++ src/grid/grid_generic.f90 | 63 +++++++++++++++-------------- src/grid/grid_physics_3d.f90 | 4 +- src/grid/grid_propagate_3d.f90 | 44 +++++++++----------- src/main/settings.f90 | 2 +- src/main/setup_rt.f90 | 5 ++- 7 files changed, 94 insertions(+), 60 deletions(-) diff --git a/hyperion/conf/conf_files.py b/hyperion/conf/conf_files.py index 6e66ab0f..ffc9f144 100644 --- a/hyperion/conf/conf_files.py +++ b/hyperion/conf/conf_files.py @@ -55,6 +55,9 @@ def __init__(self): self.set_max_reabsorptions(1000000) self.set_pda(False) self.set_mrw(False) + #DN CRAZY ADDITION + self.compute_isrf(False) + self.set_convergence(False) self.set_kill_on_absorb(False) self.set_kill_on_scatter(False) @@ -394,6 +397,27 @@ def _read_pda(self, group): def _write_pda(self, group): group.attrs['pda'] = bool2str(self.pda) + def _read_isrf(self,group): + self.isrf = str2bool(group.attrs['isrf']) + + def _write_isrf(self,group): + group.attrs['isrf'] = bool2str(self.isrf) + + + #DN CRAZY ADDITIONS + def compute_isrf(self,isrf): + + ''' + + Set whether or not to compute and save the ISRF in each cell + + If enabled, the ISRF is computed in every cell at the + frequencies of the dust opacity tables. + + ''' + self.isrf = isrf + + def set_mrw(self, mrw, gamma=1.0, inter_max=1000, warn=True): ''' Set whether to use the Modified Random Walk (MRW) approximation @@ -722,6 +746,7 @@ def read_run_conf(self, group): # not a class method because inherited self._read_max_reabsorptions(group) self._read_pda(group) self._read_mrw(group) + self._reada_isrf(group) self._read_convergence(group) self._read_kill_on_absorb(group) self._read_kill_on_scatter(group) @@ -750,6 +775,7 @@ def write_run_conf(self, group): self._write_max_reabsorptions(group) self._write_pda(group) self._write_mrw(group) + self._write_isrf(group) self._write_convergence(group) self._write_kill_on_absorb(group) self._write_kill_on_scatter(group) diff --git a/hyperion/conf/tests/test_conf_io.py b/hyperion/conf/tests/test_conf_io.py index da688d32..738ef480 100644 --- a/hyperion/conf/tests/test_conf_io.py +++ b/hyperion/conf/tests/test_conf_io.py @@ -249,6 +249,16 @@ def test_io_run_conf_mrw(value): r2.read_run_conf(v) assert r2.mrw == r1.mrw +@pytest.mark.parametrize(('value'), [True, False]) +def test_io_run_conf_isrf(value): + r1 = RunConf() + r1.compute_isrf(value) + r1.set_n_photons(1, 2) + v = virtual_file() + r1.write_run_conf(v) + r2 = RunConf() + r2.read_run_conf(v) + assert r2.isrf == r1.isrf @pytest.mark.parametrize(('value'), [False, True]) def test_io_run_conf_convergence(value): diff --git a/src/grid/grid_generic.f90 b/src/grid/grid_generic.f90 index c80fca7c..31337337 100644 --- a/src/grid/grid_generic.f90 +++ b/src/grid/grid_generic.f90 @@ -7,7 +7,7 @@ module grid_generic use grid_io, only : write_grid_3d, write_grid_4d, write_grid_5d use grid_geometry, only : geo use grid_physics, only : n_photons, last_photon_id, specific_energy_sum, specific_energy_sum_nu, specific_energy, density, density_original - use settings, only : output_n_photons, output_specific_energy, output_density, output_density_diff, physics_io_type + use settings, only : output_n_photons, output_specific_energy, output_density, output_density_diff, physics_io_type, compute_isrf !DN Crazy Additions use dust_main, only: d @@ -74,39 +74,42 @@ subroutine output_grid(group, iter, n_iter) ! DN Crazy Additions ! WRITE THE ISRF - n_cells = size(specific_energy_sum_nu, 1) - n_dust = size(specific_energy_sum_nu, 2) - n_isrf_lam = size(specific_energy_sum_nu,3) - - do i=1,d(1)%n_nu - energy_frequency_bins(i) = d(1)%nu(i) - end do - - if(trim(output_specific_energy)=='all' .or. (trim(output_specific_energy)=='last'.and.iter==n_iter)) then - !if(allocated(energy_frequency_bins)) then + if (compute_isrf) then + + n_cells = size(specific_energy_sum_nu, 1) + n_dust = size(specific_energy_sum_nu, 2) + n_isrf_lam = size(specific_energy_sum_nu,3) + + do i=1,d(1)%n_nu + energy_frequency_bins(i) = d(1)%nu(i) + end do + + if(trim(output_specific_energy)=='all' .or. (trim(output_specific_energy)=='last'.and.iter==n_iter)) then + !if(allocated(energy_frequency_bins)) then call write_grid_3d(group, 'ISRF_frequency_bins',energy_frequency_bins, geo) else call warn("output_grid","energy_frequency bins [ISRF wavelengths] array is not allocated") - !end if - end if - - if(trim(output_specific_energy)=='all' .or. (trim(output_specific_energy)=='last'.and.iter==n_iter)) then - if(allocated(specific_energy_sum_nu)) then - - select case(physics_io_type) - - case(sp) - call write_grid_5d(group, 'specific_energy_sum_nu', real(specific_energy_sum_nu, sp), geo) - case(dp) - call write_grid_5d(group, 'specific_energy_sum_nu', real(specific_energy_sum_nu, dp), geo) - case default - call error("output_grid","unexpected value of physics_io_type (should be sp or dp)") - end select - else - call warn("output_grid","specific_energy_sum_nu array is not allocated") + !end if end if - end if - + + if(trim(output_specific_energy)=='all' .or. (trim(output_specific_energy)=='last'.and.iter==n_iter)) then + if(allocated(specific_energy_sum_nu)) then + + select case(physics_io_type) + + case(sp) + call write_grid_5d(group, 'specific_energy_sum_nu', real(specific_energy_sum_nu, sp), geo) + case(dp) + call write_grid_5d(group, 'specific_energy_sum_nu', real(specific_energy_sum_nu, dp), geo) + case default + call error("output_grid","unexpected value of physics_io_type (should be sp or dp)") + end select + else + call warn("output_grid","specific_energy_sum_nu array is not allocated") + end if + end if + + endif ! DENSITY diff --git a/src/grid/grid_physics_3d.f90 b/src/grid/grid_physics_3d.f90 index cfc3179b..641a1160 100644 --- a/src/grid/grid_physics_3d.f90 +++ b/src/grid/grid_physics_3d.f90 @@ -93,12 +93,12 @@ integer function select_dust_specific_energy_rho(icell) result(id_select) id_select = sample_pdf(absorption) end function select_dust_specific_energy_rho - subroutine setup_grid_physics(group, use_mrw, use_pda) + subroutine setup_grid_physics(group, use_mrw, use_pda, compute_isrf) implicit none integer(hid_t),intent(in) :: group - logical,intent(in) :: use_mrw, use_pda + logical,intent(in) :: use_mrw, use_pda, compute_isrf integer :: n_isrf_wavelengths ! Density diff --git a/src/grid/grid_propagate_3d.f90 b/src/grid/grid_propagate_3d.f90 index 428f419f..ca2979f8 100644 --- a/src/grid/grid_propagate_3d.f90 +++ b/src/grid/grid_propagate_3d.f90 @@ -8,7 +8,7 @@ module grid_propagate use grid_physics, only : specific_energy_sum, specific_energy_sum_nu, density, n_photons, last_photon_id use sources use counters - use settings, only : frac_check => propagation_check_frequency + use settings, only : frac_check => propagation_check_frequency, compute_isrf !DN CRAZY ADDITIONS use grid_geometry, only : geo @@ -63,7 +63,6 @@ subroutine grid_integrate(p,tau_required,tau_achieved) !DN CRAZY ADDITIONS integer :: idx - real, dimension(d(1)%n_nu) :: energy_frequency_bins @@ -145,36 +144,29 @@ subroutine grid_integrate(p,tau_required,tau_achieved) p%r = p%r + tmin * p%v tau_achieved = tau_achieved + tau_cell - ! NEED INDIVIDUAL ALPHA HERE - - !print *,'[grid_propagate_3d] p%energy=',p%energy - !print *,'[grid_propagate_3d] p%nu=',p%nu - - !DN CRAZY ADDITIONS - !print *,'[grid_propagate_3d] p%nu=',p%nu - !print *,' ',minloc(abs(energy_frequency_bins-p%nu)) - !DN CRAZY ADDITIONS - idx = minloc(abs(energy_frequency_bins-p%nu),DIM=1) + if (compute_isrf) then + idx = minloc(abs(energy_frequency_bins-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) = & - & specific_energy_sum(p%icell%ic, id) + tmin * p%current_kappa(id) * p%energy - end if - - !DN CRAZY ADDITIONS - if(density(p%icell%ic,id) > 0._dp) then - specific_energy_sum_nu(p%icell%ic,id,idx) = & - & specific_energy_sum_nu(p%icell%ic,id,idx) + tmin * p%current_kappa(id) * p%energy - end if - + do id=1,n_dust + + if(density(p%icell%ic, id) > 0._dp) then + specific_energy_sum(p%icell%ic, id) = & + & specific_energy_sum(p%icell%ic, id) + tmin * p%current_kappa(id) * p%energy + end if + + + !DN CRAZY ADDITIONS + if(density(p%icell%ic,id) > 0._dp) then + specific_energy_sum_nu(p%icell%ic,id,idx) = & + & specific_energy_sum_nu(p%icell%ic,id,idx) + tmin * p%current_kappa(id) * p%energy + end if + end do + endif - end do p%on_wall = .true. p%icell = next_cell(p%icell, id_min, intersection=p%r) diff --git a/src/main/settings.f90 b/src/main/settings.f90 index daf4ad23..d8038eef 100644 --- a/src/main/settings.f90 +++ b/src/main/settings.f90 @@ -20,7 +20,7 @@ module settings integer(idp),public :: n_last_photons_dust = 0 integer(idp),public :: n_raytracing_photons_sources = 0 integer(idp),public :: n_raytracing_photons_dust = 0 - logical,public :: use_raytracing, use_mrw, use_pda + logical,public :: use_raytracing, use_mrw, use_pda, compute_isrf logical, public :: kill_on_absorb, kill_on_scatter real(dp),public :: mrw_gamma logical, public :: forced_first_interaction diff --git a/src/main/setup_rt.f90 b/src/main/setup_rt.f90 index 649799fa..941b899e 100644 --- a/src/main/setup_rt.f90 +++ b/src/main/setup_rt.f90 @@ -75,6 +75,9 @@ subroutine setup_initial(input_handle) call mp_read_keyword(input_handle, '/', 'pda', use_pda) call mp_read_keyword(input_handle, '/', 'mrw', use_mrw) + !DN CRAZY ADDITIONS + call mp_read_keyword(input_handle, '/', 'isrf', compute_isrf) + if(use_mrw) then call mp_read_keyword(input_handle, '/', 'mrw_gamma', mrw_gamma) call mp_read_keyword(input_handle, '/', 'n_inter_mrw_max', n_mrw_max) @@ -173,7 +176,7 @@ subroutine setup_initial(input_handle) call mp_close_group(g_geometry) g_physics = mp_open_group(input_handle, '/Grid/Quantities') - call setup_grid_physics(g_physics, use_mrw, use_pda) + call setup_grid_physics(g_physics, use_mrw, use_pda, compute_isrf) call mp_close_group(g_physics) call mp_read_keyword(input_handle, '/', 'physics_io_bytes', physics_io_bytes) From 5bbc6100d1c005bd2d273690e568726f35812b33 Mon Sep 17 00:00:00 2001 From: dnarayanan Date: Tue, 14 Dec 2021 11:50:39 -0500 Subject: [PATCH 09/38] getting rid of internal comments showing where ive added new lines of code --- hyperion/conf/conf_files.py | 3 --- src/grid/grid_generic.f90 | 3 --- src/grid/grid_io.f90 | 12 ++++++------ src/grid/grid_io_1d.f90 | 16 ++++++++-------- src/grid/grid_io_amr.f90 | 8 ++++---- src/grid/grid_io_amr_template.f90 | 2 +- src/grid/grid_physics_3d.f90 | 3 +-- src/grid/grid_propagate_3d.f90 | 18 ++++-------------- src/main/setup_rt.f90 | 2 -- src/mpi/mpi_routines.f90 | 2 +- 10 files changed, 25 insertions(+), 44 deletions(-) diff --git a/hyperion/conf/conf_files.py b/hyperion/conf/conf_files.py index ffc9f144..e2a84faa 100644 --- a/hyperion/conf/conf_files.py +++ b/hyperion/conf/conf_files.py @@ -55,7 +55,6 @@ def __init__(self): self.set_max_reabsorptions(1000000) self.set_pda(False) self.set_mrw(False) - #DN CRAZY ADDITION self.compute_isrf(False) self.set_convergence(False) @@ -403,8 +402,6 @@ def _read_isrf(self,group): def _write_isrf(self,group): group.attrs['isrf'] = bool2str(self.isrf) - - #DN CRAZY ADDITIONS def compute_isrf(self,isrf): ''' diff --git a/src/grid/grid_generic.f90 b/src/grid/grid_generic.f90 index 31337337..1ee0f0c1 100644 --- a/src/grid/grid_generic.f90 +++ b/src/grid/grid_generic.f90 @@ -9,7 +9,6 @@ module grid_generic use grid_physics, only : n_photons, last_photon_id, specific_energy_sum, specific_energy_sum_nu, specific_energy, density, density_original use settings, only : output_n_photons, output_specific_energy, output_density, output_density_diff, physics_io_type, compute_isrf - !DN Crazy Additions use dust_main, only: d @@ -36,7 +35,6 @@ subroutine output_grid(group, iter, n_iter) integer(hid_t),intent(in) :: group integer,intent(in) :: iter, n_iter - !DN CRAZY ADDITION integer :: n_cells, n_dust, n_isrf_lam integer :: i,j,k real, dimension(d(1)%n_nu) :: energy_frequency_bins @@ -72,7 +70,6 @@ subroutine output_grid(group, iter, n_iter) - ! DN Crazy Additions ! WRITE THE ISRF if (compute_isrf) then diff --git a/src/grid/grid_io.f90 b/src/grid/grid_io.f90 index a4aaee28..382e4d7b 100644 --- a/src/grid/grid_io.f90 +++ b/src/grid/grid_io.f90 @@ -70,7 +70,7 @@ logical function grid_exists(group, name) grid_exists = mp_path_exists(group, name) end function grid_exists - !DN CRAZY ADDITIONS + subroutine read_grid_5d_int8(group, path, array, geo) implicit none @@ -158,7 +158,7 @@ end subroutine read_grid_3d_int8 - !DN CRAZY ADDITIONS + subroutine write_grid_5d_int8(group, path, array, geo) implicit none @@ -204,7 +204,7 @@ subroutine write_grid_3d_int8(group, path, array, geo) end subroutine write_grid_3d_int8 - !DN CRAZY ADDITIONS + subroutine read_grid_5d_int(group, path, array, geo) implicit none @@ -292,7 +292,7 @@ subroutine read_grid_3d_int(group, path, array, geo) end subroutine read_grid_3d_int - !DN CRAZY ADDITION + subroutine write_grid_5d_int(group, path, array, geo) @@ -337,7 +337,7 @@ subroutine write_grid_3d_int(group, path, array, geo) end subroutine write_grid_3d_int - !DN CRAZY ADDITIONS + subroutine read_grid_5d_dp(group, path, array, geo) implicit none @@ -424,7 +424,7 @@ subroutine read_grid_3d_dp(group, path, array, geo) end subroutine read_grid_3d_dp - !CRAZY DN ADDITIONS + subroutine write_grid_5d_dp(group, path, array, geo) implicit none diff --git a/src/grid/grid_io_1d.f90 b/src/grid/grid_io_1d.f90 index 98779162..0c642387 100644 --- a/src/grid/grid_io_1d.f90 +++ b/src/grid/grid_io_1d.f90 @@ -71,7 +71,7 @@ logical function grid_exists(group, name) grid_exists = mp_path_exists(group, name) end function grid_exists - !DN CRAZY ADDITIONS + subroutine read_grid_5d_int8(group, path, array, geo) implicit none @@ -137,7 +137,7 @@ subroutine read_grid_3d_int8(group, path, array, geo) end subroutine read_grid_3d_int8 - !DN CRAZY ADDITIONS + subroutine write_grid_5d_int8(group, path, array, geo) implicit none @@ -181,7 +181,7 @@ subroutine write_grid_3d_int8(group, path, array, geo) end subroutine write_grid_3d_int8 - !DN CRAZY ADDITIONS + subroutine read_grid_5d_int(group, path, array, geo) implicit none @@ -248,7 +248,7 @@ end subroutine read_grid_3d_int - !DN CRAZY ADDITION + subroutine write_grid_5d_int(group, path, array, geo) @@ -293,7 +293,7 @@ subroutine write_grid_3d_int(group, path, array, geo) end subroutine write_grid_3d_int - !DN CRAZY ADDITIONS + subroutine read_grid_5d_dp(group, path, array, geo) implicit none @@ -359,7 +359,7 @@ subroutine read_grid_3d_dp(group, path, array, geo) end subroutine read_grid_3d_dp - !CRAZY DN ADDITIONS + subroutine write_grid_5d_dp(group, path, array, geo) implicit none @@ -404,7 +404,7 @@ subroutine write_grid_3d_dp(group, path, array, geo) end subroutine write_grid_3d_dp - !DN CRAZY ADDITIONS + subroutine read_grid_5d_sp(group, path, array, geo) implicit none @@ -470,7 +470,7 @@ subroutine read_grid_3d_sp(group, path, array, geo) end subroutine read_grid_3d_sp - !CRAZY DN ADDITIONS + subroutine write_grid_5d_sp(group, path, array, geo) implicit none diff --git a/src/grid/grid_io_amr.f90 b/src/grid/grid_io_amr.f90 index 0258db78..0d3e3135 100644 --- a/src/grid/grid_io_amr.f90 +++ b/src/grid/grid_io_amr.f90 @@ -213,7 +213,7 @@ end subroutine test end subroutine read_grid_3d_int8 - !DN CRAZY ADDITIONS + subroutine write_grid_5d_int8(group, path, array, geo) implicit none @@ -458,7 +458,7 @@ end subroutine test end subroutine read_grid_3d_int - !DN CRAZY ADDITIONS + subroutine write_grid_5d_int(group, path, array, geo) implicit none @@ -698,7 +698,7 @@ end subroutine test end subroutine read_grid_3d_dp -!DN CRAZY ADDITIONS + subroutine write_grid_5d_dp(group, path, array, geo) implicit none @@ -937,7 +937,7 @@ end subroutine test end subroutine read_grid_3d_sp -!DN CRAZY ADDITIONS + subroutine write_grid_5d_sp(group, path, array, geo) implicit none diff --git a/src/grid/grid_io_amr_template.f90 b/src/grid/grid_io_amr_template.f90 index e3abfedc..5e06a8c1 100644 --- a/src/grid/grid_io_amr_template.f90 +++ b/src/grid/grid_io_amr_template.f90 @@ -158,7 +158,7 @@ end subroutine test end subroutine read_grid_3d_ - !DN CRAZY ADDITIONS + subroutine write_grid_5d_(group, path, array, geo) implicit none diff --git a/src/grid/grid_physics_3d.f90 b/src/grid/grid_physics_3d.f90 index 641a1160..938e54a4 100644 --- a/src/grid/grid_physics_3d.f90 +++ b/src/grid/grid_physics_3d.f90 @@ -196,8 +196,7 @@ subroutine setup_grid_physics(group, use_mrw, use_pda, compute_isrf) allocate(specific_energy_sum(geo%n_cells, n_dust)) specific_energy_sum = 0._dp - - !DN CRAZY ADDITION + ! Set up basics for ISRF calculation n_isrf_wavelengths = d(1)%n_nu allocate(specific_energy_sum_nu(geo%n_cells, n_dust, n_isrf_wavelengths)) specific_energy_sum_nu = 0._dp diff --git a/src/grid/grid_propagate_3d.f90 b/src/grid/grid_propagate_3d.f90 index ca2979f8..232a6ecd 100644 --- a/src/grid/grid_propagate_3d.f90 +++ b/src/grid/grid_propagate_3d.f90 @@ -10,7 +10,6 @@ module grid_propagate use counters use settings, only : frac_check => propagation_check_frequency, compute_isrf - !DN CRAZY ADDITIONS use grid_geometry, only : geo implicit none @@ -60,8 +59,6 @@ subroutine grid_integrate(p,tau_required,tau_achieved) integer :: source_id - - !DN CRAZY ADDITIONS integer :: idx real, dimension(d(1)%n_nu) :: energy_frequency_bins @@ -144,7 +141,8 @@ subroutine grid_integrate(p,tau_required,tau_achieved) p%r = p%r + tmin * p%v tau_achieved = tau_achieved + tau_cell - !DN CRAZY ADDITIONS + + ! Compute the ISRF if (compute_isrf) then idx = minloc(abs(energy_frequency_bins-p%nu),DIM=1) @@ -158,7 +156,6 @@ subroutine grid_integrate(p,tau_required,tau_achieved) end if - !DN CRAZY ADDITIONS if(density(p%icell%ic,id) > 0._dp) then specific_energy_sum_nu(p%icell%ic,id,idx) = & & specific_energy_sum_nu(p%icell%ic,id,idx) + tmin * p%current_kappa(id) * p%energy @@ -250,7 +247,6 @@ subroutine grid_integrate_noenergy(p,tau_required,tau_achieved) integer :: source_id - !DN CRAZY ADDITIONS integer :: id integer :: idx real, dimension(d(0)%n_nu) :: energy_frequency_bins @@ -309,20 +305,14 @@ subroutine grid_integrate_noenergy(p,tau_required,tau_achieved) chi_rho_total = 0._dp - - !DN CRAZY ADDITIONS - !idx = minloc(abs(energy_frequency_bins-p%nu),DIM=1) - !print *,'[grid_propagate_3d last iteration] p%energy=',p%energy - !print *,'[grid_propagate_3d last iteration] p%nu=',p%nu - !print *, "shape(specific_energy_sum_nu", shape(specific_energy_sum_nu) - + !Compute the ISRF + !Figure out what frequency bin in the ISRF calculation the current photon's frequency is closest to idx = minloc(abs(energy_frequency_bins-p%nu),DIM=1) do id=1,n_dust chi_rho_total = chi_rho_total + p%current_chi(id) * density(p%icell%ic, id) - !DN CRAZY ADDITIONS if(density(p%icell%ic,id) > 0._dp) then specific_energy_sum_nu(p%icell%ic,id,idx) = & & specific_energy_sum_nu(p%icell%ic,id,idx) + tmin * p%current_kappa(id) * p%energy diff --git a/src/main/setup_rt.f90 b/src/main/setup_rt.f90 index 941b899e..5be59417 100644 --- a/src/main/setup_rt.f90 +++ b/src/main/setup_rt.f90 @@ -74,8 +74,6 @@ subroutine setup_initial(input_handle) call mp_read_keyword(input_handle, '/', 'pda', use_pda) call mp_read_keyword(input_handle, '/', 'mrw', use_mrw) - - !DN CRAZY ADDITIONS call mp_read_keyword(input_handle, '/', 'isrf', compute_isrf) if(use_mrw) then diff --git a/src/mpi/mpi_routines.f90 b/src/mpi/mpi_routines.f90 index 8f4e8a15..32e1c777 100644 --- a/src/mpi/mpi_routines.f90 +++ b/src/mpi/mpi_routines.f90 @@ -266,7 +266,7 @@ subroutine mp_collect_physical_arrays() real(dp) :: tmp real(dp) :: dummy_dp integer(idp) :: dummy_idp - !DN crazy additions + real(dp),allocatable :: tmp_3d(:,:,:) real(dp),allocatable :: tmp_2d(:,:) integer(idp),allocatable :: tmp_int_1d(:) From 935b6d1d7a47d5fe506b1ed8fc3d97f640c2bbd7 Mon Sep 17 00:00:00 2001 From: dnarayanan Date: Mon, 9 May 2022 08:51:16 -0400 Subject: [PATCH 10/38] some bugs corrected fixing the ISRF calculation. but some issues remain now when running without ISRF calculation turned on --- hyperion/model/model.py | 11 ++++- src/grid/grid_generic.f90 | 28 +++++++---- src/grid/grid_io.f90 | 3 ++ src/grid/grid_physics_3d.f90 | 96 ++++++++++++++++++++++++++++++++++-- 4 files changed, 122 insertions(+), 16 deletions(-) diff --git a/hyperion/model/model.py b/hyperion/model/model.py index 23a6edba..fc07a678 100644 --- a/hyperion/model/model.py +++ b/hyperion/model/model.py @@ -223,7 +223,7 @@ def use_geometry(self, filename): # Close the file f.close() - def use_quantities(self, filename, quantities=['density', 'specific_energy'], + def use_quantities(self, filename, quantities=['density', 'specific_energy', 'specific_energy_nu'], use_minimum_specific_energy=True, use_dust=True, copy=True, only_initial=False): ''' @@ -297,6 +297,13 @@ def use_quantities(self, filename, quantities=['density', 'specific_energy'], else: quantities_path['specific_energy'] = last_iteration + if 'specific_energy_nu' in quantities: + if only_initial or last_iteration is None: + if 'specific_energy_nu' in f['/Input/Grid/Quantities']: + quantities_path['specific_energy_nu'] = '/Input/Grid/Quantities' + else: + quantities_path['specific_energy_nu'] = last_iteration + # Minimum specific energy if use_minimum_specific_energy: minimum_specific_energy_path = '/Input/Grid/Quantities' @@ -812,6 +819,7 @@ class (e.g. ``SphericalDust``, ``IsotropicDust``, ...). self.grid['density'] = [] if specific_energy is not None: self.grid['specific_energy'] = [] + self.grid['specific_energy_nu'] = [] # Check whether the density can be added to an existing one if merge_if_possible: @@ -853,6 +861,7 @@ class (e.g. ``SphericalDust``, ``IsotropicDust``, ...). # Set specific energy if specified if specific_energy is not None: self.grid['specific_energy'].append(specific_energy) + self.grid['specific_energy_nu'].append(specific_energy) def set_cartesian_grid(self, x_wall, y_wall, z_wall): self.set_grid(CartesianGrid(x_wall, y_wall, z_wall)) diff --git a/src/grid/grid_generic.f90 b/src/grid/grid_generic.f90 index 1ee0f0c1..710d70ed 100644 --- a/src/grid/grid_generic.f90 +++ b/src/grid/grid_generic.f90 @@ -6,7 +6,7 @@ module grid_generic use grid_io, only : write_grid_3d, write_grid_4d, write_grid_5d use grid_geometry, only : geo - use grid_physics, only : n_photons, last_photon_id, specific_energy_sum, specific_energy_sum_nu, specific_energy, density, density_original + use grid_physics, only : n_photons, last_photon_id, specific_energy_sum, specific_energy_sum_nu, specific_energy, specific_energy_nu, density, density_original use settings, only : output_n_photons, output_specific_energy, output_density, output_density_diff, physics_io_type, compute_isrf use dust_main, only: d @@ -39,6 +39,7 @@ subroutine output_grid(group, iter, n_iter) integer :: i,j,k real, dimension(d(1)%n_nu) :: energy_frequency_bins + if(main_process()) write(*,'(" [output_grid] outputting grid arrays for iteration")') ! NUMBER OF PHOTONS IN EACH CELL @@ -73,9 +74,9 @@ subroutine output_grid(group, iter, n_iter) ! WRITE THE ISRF if (compute_isrf) then - n_cells = size(specific_energy_sum_nu, 1) - n_dust = size(specific_energy_sum_nu, 2) - n_isrf_lam = size(specific_energy_sum_nu,3) + n_cells = size(specific_energy_nu, 1) + n_dust = size(specific_energy_nu, 2) + n_isrf_lam = size(specific_energy_nu,3) do i=1,d(1)%n_nu energy_frequency_bins(i) = d(1)%nu(i) @@ -84,25 +85,32 @@ subroutine output_grid(group, iter, n_iter) if(trim(output_specific_energy)=='all' .or. (trim(output_specific_energy)=='last'.and.iter==n_iter)) then !if(allocated(energy_frequency_bins)) then call write_grid_3d(group, 'ISRF_frequency_bins',energy_frequency_bins, geo) - else - call warn("output_grid","energy_frequency bins [ISRF wavelengths] array is not allocated") + !else + !call warn("output_grid","energy_frequency bins [ISRF wavelengths] array is not allocated") !end if end if if(trim(output_specific_energy)=='all' .or. (trim(output_specific_energy)=='last'.and.iter==n_iter)) then - if(allocated(specific_energy_sum_nu)) then + if(allocated(specific_energy_nu)) then + print *,'[GRID_GENERIC.F90] ENTERING BLOCK FOR WRITING SPECIFIC_ENERGY_NU' + select case(physics_io_type) case(sp) - call write_grid_5d(group, 'specific_energy_sum_nu', real(specific_energy_sum_nu, sp), geo) + print *,'[GRID_GENERIC.F90] diving into write_grid_5d sp' + call write_grid_5d(group, 'specific_energy_nu', real(specific_energy_nu, sp), geo) + case(dp) - call write_grid_5d(group, 'specific_energy_sum_nu', real(specific_energy_sum_nu, dp), geo) + print *,'[GRID_GENERIC.F90] diving into write_grid_5d dp' + call write_grid_5d(group, 'specific_energy_nu', real(specific_energy_nu, dp), geo) + print *,'[GRID_GENERIC.F90] finished with write_grid_5d dp' + case default call error("output_grid","unexpected value of physics_io_type (should be sp or dp)") end select else - call warn("output_grid","specific_energy_sum_nu array is not allocated") + call warn("output_grid","specific_energy_nu array is not allocated") end if end if diff --git a/src/grid/grid_io.f90 b/src/grid/grid_io.f90 index 382e4d7b..ea8ccb77 100644 --- a/src/grid/grid_io.f90 +++ b/src/grid/grid_io.f90 @@ -163,11 +163,13 @@ subroutine write_grid_5d_int8(group, path, array, geo) implicit none + integer(hid_t), intent(in) :: group character(len=*), intent(in) :: path integer(idp), intent(in) :: array(:,:,:) type(grid_geometry_desc),intent(in) :: geo + call mp_write_array(group, path, reshape(array, (/geo%n1, geo%n2, geo%n3, size(array,2), size(array,3)/))) call mp_write_keyword(group, path, 'geometry', geo%id) @@ -434,6 +436,7 @@ subroutine write_grid_5d_dp(group, path, array, geo) real(dp), intent(in) :: array(:,:,:) type(grid_geometry_desc),intent(in) :: geo + call mp_write_array(group, path, reshape(array, (/geo%n1, geo%n2, geo%n3, size(array,2), size(array,3)/))) call mp_write_keyword(group, path, 'geometry', geo%id) diff --git a/src/grid/grid_physics_3d.f90 b/src/grid/grid_physics_3d.f90 index 938e54a4..cd121b18 100644 --- a/src/grid/grid_physics_3d.f90 +++ b/src/grid/grid_physics_3d.f90 @@ -3,7 +3,7 @@ module grid_physics use core_lib use type_photon use type_grid_cell - use grid_io, only : read_grid_4d, grid_exists + use grid_io, only : read_grid_4d, grid_exists, read_grid_5d use dust_main ! many variables and routines use type_dust use grid_geometry @@ -36,10 +36,12 @@ module grid_physics integer(idp),allocatable, public :: n_photons(:) integer(idp),allocatable, public :: last_photon_id(:) real(dp),allocatable, public :: specific_energy(:,:) + real(dp),allocatable, public :: specific_energy_nu(:,:,:) real(dp),allocatable, public :: specific_energy_sum(:,:) real(dp),allocatable, public :: specific_energy_sum_nu(:,:,:) real(dp),allocatable, public :: specific_energy_additional(:,:) + real(dp),allocatable, public :: specific_energy_additional_nu(:,:,:) real(dp),allocatable, public :: energy_abs_tot(:) real(dp),allocatable, public :: minimum_specific_energy(:) @@ -53,7 +55,9 @@ module grid_physics ! Temporary variables (used for convenience in external modules) real(dp), allocatable, public :: tmp_column_density(:) + !Counting indices integer :: id + integer :: idx logical :: debug = .false. @@ -101,9 +105,11 @@ subroutine setup_grid_physics(group, use_mrw, use_pda, compute_isrf) logical,intent(in) :: use_mrw, use_pda, compute_isrf integer :: n_isrf_wavelengths + n_isrf_wavelengths = d(1)%n_nu ! Density allocate(density(geo%n_cells, n_dust)) allocate(specific_energy(geo%n_cells, n_dust)) + allocate(specific_energy_nu(geo%n_cells,n_dust,n_isrf_wavelengths)) if(n_dust > 0) then @@ -147,10 +153,13 @@ subroutine setup_grid_physics(group, use_mrw, use_pda, compute_isrf) ! Read in specific_energy call read_grid_4d(group, 'specific_energy', specific_energy, geo) + call read_grid_5d(group, 'specific_energy_nu', specific_energy_nu, geo) ! Check number of dust types for specific_energy if(size(specific_energy, 2).ne.n_dust) call error("setup_grid","specific_energy array has wrong number of dust types") + + ! Reset specific energy to zero in masked cells if(geo%masked) then if(main_process()) write(*, '(" [grid_physics] applying mask to specific_energy grid")') @@ -158,19 +167,39 @@ subroutine setup_grid_physics(group, use_mrw, use_pda, compute_isrf) where(.not.geo%mask) specific_energy(:, id) = 0. end where + + if (compute_isrf) then + do idx=1,n_isrf_wavelengths + where(.not.geo%mask) + specific_energy_nu(:, id, idx) = 0. + endwhere + end do + end if + end do end if if(trim(specific_energy_type) == 'additional') then allocate(specific_energy_additional(geo%n_cells, n_dust)) + if (compute_isrf) then + allocate(specific_energy_additional_nu(geo%n_cells,n_dust,n_isrf_wavelengths)) + end if ! We store a copy of the initial specific energy in a separate ! array, and we set the specific energy to the minimum specific ! 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_nu = specific_energy_nu do id=1,n_dust specific_energy(:,id) = minimum_specific_energy(id) - end do + + if (compute_isrf) then + do idx=1,n_isrf_wavelengths + specific_energy_nu(:,id,idx) = minimum_specific_energy(id) + end do + end if + + end do end if @@ -183,6 +212,13 @@ subroutine setup_grid_physics(group, use_mrw, use_pda, compute_isrf) ! Set all specific_energy to minimum requested do id=1,n_dust specific_energy(:,id) = minimum_specific_energy(id) + + if (compute_isrf) then + do idx = 1,n_isrf_wavelengths + specific_energy_nu(:,id,idx) = minimum_specific_energy(id) + end do + end if + end do end if @@ -197,7 +233,7 @@ subroutine setup_grid_physics(group, use_mrw, use_pda, compute_isrf) specific_energy_sum = 0._dp ! Set up basics for ISRF calculation - n_isrf_wavelengths = d(1)%n_nu + !n_isrf_wavelengths = d(1)%n_nu allocate(specific_energy_sum_nu(geo%n_cells, n_dust, n_isrf_wavelengths)) specific_energy_sum_nu = 0._dp @@ -267,6 +303,9 @@ subroutine sublimate_dust() implicit none integer :: ic, id integer :: reset + integer :: n_isrf_wavelengths + + n_isrf_wavelengths = d(1)%n_nu reset = 0 @@ -280,7 +319,15 @@ subroutine sublimate_dust() density(ic, id) = 0. specific_energy(ic, id) = minimum_specific_energy(id) reset = reset + 1 + + if (compute_isrf) then + do idx=1,n_isrf_wavelengths + specific_energy_nu(ic,id,idx) = minimum_specific_energy(id) + end do + end if + end if + end do if(reset > 0) write(*,'(" [sublimate_dust] dust removed in ",I0," cells")') reset @@ -294,6 +341,14 @@ subroutine sublimate_dust() & / chi_rosseland(id, d(id)%sublimation_specific_energy))**2 specific_energy(ic,id) = d(id)%sublimation_specific_energy reset = reset + 1 + + + if (compute_isrf) then + do idx=1,n_isrf_wavelengths + specific_energy_nu(ic,id,idx) = minimum_specific_energy(id) + end do + end if + end if end do @@ -306,6 +361,14 @@ subroutine sublimate_dust() specific_energy(ic, id) = d(id)%sublimation_specific_energy reset = reset + 1 end if + + + if (compute_isrf) then + do idx=1,n_isrf_wavelengths + specific_energy_nu(ic,id,idx) = minimum_specific_energy(id) + end do + end if + end do if(reset > 0) write(*,'(" [sublimate_dust] capping dust specific_energy in ",I0," cells")') reset @@ -326,12 +389,23 @@ subroutine update_energy_abs(scale) real(dp), intent(in) :: scale - integer :: id + integer :: id,idx + integer :: n_isrf_wavelengths + + n_isrf_wavelengths = d(1)%n_nu if(main_process()) write(*,'(" [grid_physics] updating energy_abs")') do id=1,n_dust specific_energy(:,id) = specific_energy_sum(:,id) * scale / geo%volume + + if (compute_isrf) then + do idx=1,n_isrf_wavelengths + specific_energy_nu(:,id,idx) = specific_energy_sum_nu(:,id,idx) * scale/geo%volume + end do + end if + + where(geo%volume == 0._dp) specific_energy(:,id) = 0._dp end where @@ -341,10 +415,17 @@ subroutine update_energy_abs(scale) write(*,'(" [update_energy_abs] ",I0," cells have no energy")') count(specific_energy==0.and.density>0.) end if + ! Add in additional source of heating if(trim(specific_energy_type) == 'additional') then if(main_process()) write(*,'(" [grid_physics] adding additional heating source")') specific_energy = specific_energy + specific_energy_additional + + + if (compute_isrf) then + specific_energy_nu = specific_energy_nu + specific_energy_additional_nu + end if + end if call update_energy_abs_tot() @@ -367,7 +448,10 @@ subroutine check_energy_abs() if(main_process()) call warn("update_energy_abs","specific_energy below minimum requested in some cells - resetting") where(specific_energy(:,id) < minimum_specific_energy(id)) specific_energy(:,id) = minimum_specific_energy(id) + + end where + end if if(any(specific_energy(:,id) < d(id)%specific_energy(1))) then @@ -375,7 +459,8 @@ subroutine check_energy_abs() if(main_process()) call warn("update_energy_abs","specific_energy below minimum allowed in some cells - resetting") where(specific_energy(:,id) < d(id)%specific_energy(1)) specific_energy(:,id) = d(id)%specific_energy(1) - end where + + end where else if(main_process()) call warn("update_energy_abs","specific_energy below minimum allowed in some cells - will pick closest emissivities") end if @@ -386,6 +471,7 @@ subroutine check_energy_abs() if(main_process()) call warn("update_energy_abs","specific_energy above maximum allowed in some cells - resetting") where(specific_energy(:,id) > d(id)%specific_energy(d(id)%n_e)) specific_energy(:,id) = d(id)%specific_energy(d(id)%n_e) + end where else if(main_process()) call warn("update_energy_abs","specific_energy above maximum allowed in some cells - will pick closest emissivities") From e2f49f87a86c71c05ee6859e695fa903b272d5f5 Mon Sep 17 00:00:00 2001 From: dnarayanan Date: Mon, 9 May 2022 12:22:59 -0400 Subject: [PATCH 11/38] bugs fixed that kept code from running when ISRF is turned off --- hyperion/model/model.py | 30 ++++++++++++++++++------------ src/grid/grid_generic.f90 | 36 ++++++++++++++++++++---------------- src/grid/grid_physics_3d.f90 | 1 - 3 files changed, 38 insertions(+), 29 deletions(-) diff --git a/hyperion/model/model.py b/hyperion/model/model.py index fc07a678..567dd9aa 100644 --- a/hyperion/model/model.py +++ b/hyperion/model/model.py @@ -297,12 +297,13 @@ def use_quantities(self, filename, quantities=['density', 'specific_energy', 'sp else: quantities_path['specific_energy'] = last_iteration - if 'specific_energy_nu' in quantities: - if only_initial or last_iteration is None: - if 'specific_energy_nu' in f['/Input/Grid/Quantities']: - quantities_path['specific_energy_nu'] = '/Input/Grid/Quantities' - else: - quantities_path['specific_energy_nu'] = last_iteration + if self.compute_isrf == True: + if 'specific_energy_nu' in quantities: + if only_initial or last_iteration is None: + if 'specific_energy_nu' in f['/Input/Grid/Quantities']: + quantities_path['specific_energy_nu'] = '/Input/Grid/Quantities' + else: + quantities_path['specific_energy_nu'] = last_iteration # Minimum specific energy if use_minimum_specific_energy: @@ -320,10 +321,11 @@ def use_quantities(self, filename, quantities=['density', 'specific_energy', 'sp if 'specific_energy' in f['/Grid/Quantities']: quantities_path['specific_energy'] = '/Grid/Quantities' - if 'specific_energy_nu' in quantities: - if 'specific_energy_nu' in f['/Grid/Quantities']: - quantities_path['specific_energy_nu'] = '/Grid/Quantities' - + if self.compute_isrf == True: + if 'specific_energy_nu' in quantities: + if 'specific_energy_nu' in f['/Grid/Quantities']: + quantities_path['specific_energy_nu'] = '/Grid/Quantities' + # Minimum specific energy if use_minimum_specific_energy: minimum_specific_energy_path = '/Grid/Quantities' @@ -819,7 +821,9 @@ class (e.g. ``SphericalDust``, ``IsotropicDust``, ...). self.grid['density'] = [] if specific_energy is not None: self.grid['specific_energy'] = [] - self.grid['specific_energy_nu'] = [] + + if self.compute_isrf == True: + self.grid['specific_energy_nu'] = [] # Check whether the density can be added to an existing one if merge_if_possible: @@ -861,7 +865,9 @@ class (e.g. ``SphericalDust``, ``IsotropicDust``, ...). # Set specific energy if specified if specific_energy is not None: self.grid['specific_energy'].append(specific_energy) - self.grid['specific_energy_nu'].append(specific_energy) + + if self.compute_isrf == True: + self.grid['specific_energy_nu'].append(specific_energy) def set_cartesian_grid(self, x_wall, y_wall, z_wall): self.set_grid(CartesianGrid(x_wall, y_wall, z_wall)) diff --git a/src/grid/grid_generic.f90 b/src/grid/grid_generic.f90 index 710d70ed..2dc2457c 100644 --- a/src/grid/grid_generic.f90 +++ b/src/grid/grid_generic.f90 @@ -93,24 +93,28 @@ subroutine output_grid(group, iter, n_iter) if(trim(output_specific_energy)=='all' .or. (trim(output_specific_energy)=='last'.and.iter==n_iter)) then if(allocated(specific_energy_nu)) then - print *,'[GRID_GENERIC.F90] ENTERING BLOCK FOR WRITING SPECIFIC_ENERGY_NU' - select case(physics_io_type) + if (compute_isrf) then + print *,'[GRID_GENERIC.F90] ENTERING BLOCK FOR WRITING SPECIFIC_ENERGY_NU' - case(sp) - print *,'[GRID_GENERIC.F90] diving into write_grid_5d sp' - call write_grid_5d(group, 'specific_energy_nu', real(specific_energy_nu, sp), geo) - - case(dp) - print *,'[GRID_GENERIC.F90] diving into write_grid_5d dp' - call write_grid_5d(group, 'specific_energy_nu', real(specific_energy_nu, dp), geo) - print *,'[GRID_GENERIC.F90] finished with write_grid_5d dp' - - case default - call error("output_grid","unexpected value of physics_io_type (should be sp or dp)") - end select - else - call warn("output_grid","specific_energy_nu array is not allocated") + select case(physics_io_type) + + case(sp) + print *,'[GRID_GENERIC.F90] diving into write_grid_5d sp' + call write_grid_5d(group, 'specific_energy_nu', real(specific_energy_nu, sp), geo) + + case(dp) + print *,'[GRID_GENERIC.F90] diving into write_grid_5d dp' + call write_grid_5d(group, 'specific_energy_nu', real(specific_energy_nu, dp), geo) + print *,'[GRID_GENERIC.F90] finished with write_grid_5d dp' + + case default + call error("output_grid","unexpected value of physics_io_type (should be sp or dp)") + end select + else + call warn("output_grid","specific_energy_nu array is not allocated") + + end if end if end if diff --git a/src/grid/grid_physics_3d.f90 b/src/grid/grid_physics_3d.f90 index cd121b18..1d8c804e 100644 --- a/src/grid/grid_physics_3d.f90 +++ b/src/grid/grid_physics_3d.f90 @@ -153,7 +153,6 @@ subroutine setup_grid_physics(group, use_mrw, use_pda, compute_isrf) ! Read in specific_energy call read_grid_4d(group, 'specific_energy', specific_energy, geo) - call read_grid_5d(group, 'specific_energy_nu', specific_energy_nu, geo) ! Check number of dust types for specific_energy if(size(specific_energy, 2).ne.n_dust) call error("setup_grid","specific_energy array has wrong number of dust types") From 7f86b0d3857a5c6dea383dafd603e54c9d512559 Mon Sep 17 00:00:00 2001 From: dnarayanan Date: Tue, 10 May 2022 12:14:36 -0400 Subject: [PATCH 12/38] getting rid of read_grid_5d since its not actually used --- src/grid/grid_io.f90 | 131 ------------------------- src/grid/grid_io_1d.f90 | 98 ------------------- src/grid/grid_io_amr.f90 | 176 ---------------------------------- src/grid/grid_io_template.f90 | 39 -------- src/grid/grid_physics_3d.f90 | 2 +- 5 files changed, 1 insertion(+), 445 deletions(-) diff --git a/src/grid/grid_io.f90 b/src/grid/grid_io.f90 index ea8ccb77..0cba3f1d 100644 --- a/src/grid/grid_io.f90 +++ b/src/grid/grid_io.f90 @@ -12,7 +12,6 @@ module grid_io public :: grid_exists public :: read_grid_3d public :: read_grid_4d - public :: read_grid_5d public :: write_grid_3d public :: write_grid_4d public :: write_grid_5d @@ -31,14 +30,6 @@ module grid_io module procedure read_grid_4d_int8 end interface read_grid_4d - interface read_grid_5d - module procedure read_grid_5d_sp - module procedure read_grid_5d_dp - module procedure read_grid_5d_int - module procedure read_grid_5d_int8 - end interface read_grid_5d - - interface write_grid_3d module procedure write_grid_3d_sp module procedure write_grid_3d_dp @@ -71,36 +62,6 @@ logical function grid_exists(group, name) end function grid_exists - subroutine read_grid_5d_int8(group, path, array, geo) - - implicit none - - integer(hid_t), intent(in) :: group - character(len=*), intent(in) :: path - integer(idp), intent(out) :: array(:,:,:) - type(grid_geometry_desc),intent(in) :: geo - integer(idp), allocatable :: array5d(:,:,:,:,:) - integer :: n_cells, n_dust, n_isrf_lam - - character(len=32) :: geometry_id_check - - call mp_read_keyword(group,path, 'geometry', geometry_id_check) - if(geometry_id_check.ne.geo%id) then - call error("read_grid", "geometry IDs do not match") - end if - call mp_read_array_auto(group,path, array5d) - - if(any(is_nan(array5d))) call error("read_grid_5d", "NaN values in 5D array") - - n_cells = size(array, 1) - n_dust = size(array, 2) - n_isrf_lam = size(array,3) - - array = reshape(array5d, (/n_cells, n_dust, n_isrf_lam/)) - - end subroutine read_grid_5d_int8 - - subroutine read_grid_4d_int8(group, path, array, geo) implicit none @@ -207,37 +168,6 @@ subroutine write_grid_3d_int8(group, path, array, geo) end subroutine write_grid_3d_int8 - subroutine read_grid_5d_int(group, path, array, geo) - - implicit none - - integer(hid_t), intent(in) :: group - character(len=*), intent(in) :: path - integer, intent(out) :: array(:,:,:) - type(grid_geometry_desc),intent(in) :: geo - integer, allocatable :: array5d(:,:,:,:,:) - integer :: n_cells, n_dust, n_isrf_lam - - character(len=32) :: geometry_id_check - - call mp_read_keyword(group,path, 'geometry', geometry_id_check) - if(geometry_id_check.ne.geo%id) then - call error("read_grid", "geometry IDs do not match") - end if - call mp_read_array_auto(group,path, array5d) - - if(any(is_nan(array5d))) call error("read_grid_5d", "NaN values in 5D array") - - n_cells = size(array, 1) - n_dust = size(array, 2) - n_isrf_lam = size(array,3) - - array = reshape(array5d, (/n_cells, n_dust, n_isrf_lam/)) - - end subroutine read_grid_5d_int - - - subroutine read_grid_4d_int(group, path, array, geo) implicit none @@ -340,36 +270,6 @@ subroutine write_grid_3d_int(group, path, array, geo) end subroutine write_grid_3d_int - subroutine read_grid_5d_dp(group, path, array, geo) - - implicit none - - integer(hid_t), intent(in) :: group - character(len=*), intent(in) :: path - real(dp), intent(out) :: array(:,:,:) - type(grid_geometry_desc),intent(in) :: geo - real(dp), allocatable :: array5d(:,:,:,:,:) - integer :: n_cells, n_dust, n_isrf_lam - - character(len=32) :: geometry_id_check - - call mp_read_keyword(group,path, 'geometry', geometry_id_check) - if(geometry_id_check.ne.geo%id) then - call error("read_grid", "geometry IDs do not match") - end if - call mp_read_array_auto(group,path, array5d) - - if(any(is_nan(array5d))) call error("read_grid_5d", "NaN values in 5D array") - - n_cells = size(array, 1) - n_dust = size(array, 2) - n_isrf_lam = size(array,3) - - array = reshape(array5d, (/n_cells, n_dust, n_isrf_lam/)) - - end subroutine read_grid_5d_dp - - subroutine read_grid_4d_dp(group, path, array, geo) implicit none @@ -472,37 +372,6 @@ subroutine write_grid_3d_dp(group, path, array, geo) end subroutine write_grid_3d_dp - !DN CRAZY ADDITIONS - subroutine read_grid_5d_sp(group, path, array, geo) - - implicit none - - integer(hid_t), intent(in) :: group - character(len=*), intent(in) :: path - real(sp), intent(out) :: array(:,:,:) - type(grid_geometry_desc),intent(in) :: geo - real(sp), allocatable :: array5d(:,:,:,:,:) - integer :: n_cells, n_dust, n_isrf_lam - - character(len=32) :: geometry_id_check - - call mp_read_keyword(group,path, 'geometry', geometry_id_check) - if(geometry_id_check.ne.geo%id) then - call error("read_grid", "geometry IDs do not match") - end if - call mp_read_array_auto(group,path, array5d) - - if(any(is_nan(array5d))) call error("read_grid_5d", "NaN values in 5D array") - - n_cells = size(array, 1) - n_dust = size(array, 2) - n_isrf_lam = size(array,3) - - array = reshape(array5d, (/n_cells, n_dust, n_isrf_lam/)) - - end subroutine read_grid_5d_sp - - subroutine read_grid_4d_sp(group, path, array, geo) implicit none diff --git a/src/grid/grid_io_1d.f90 b/src/grid/grid_io_1d.f90 index 0c642387..7388046d 100644 --- a/src/grid/grid_io_1d.f90 +++ b/src/grid/grid_io_1d.f90 @@ -12,7 +12,6 @@ module grid_io public :: grid_exists public :: read_grid_3d public :: read_grid_4d - public :: read_grid_5d public :: write_grid_3d public :: write_grid_4d public :: write_grid_5d @@ -32,14 +31,6 @@ module grid_io module procedure read_grid_4d_int8 end interface read_grid_4d - interface read_grid_5d - module procedure read_grid_5d_sp - module procedure read_grid_5d_dp - module procedure read_grid_5d_int - module procedure read_grid_5d_int8 - end interface read_grid_5d - - interface write_grid_3d module procedure write_grid_3d_sp module procedure write_grid_3d_dp @@ -72,28 +63,6 @@ logical function grid_exists(group, name) end function grid_exists - subroutine read_grid_5d_int8(group, path, array, geo) - - implicit none - - integer(hid_t), intent(in) :: group - character(len=*), intent(in) :: path - integer(idp), intent(out) :: array(:,:,:) - type(grid_geometry_desc),intent(in) :: geo - - character(len=32) :: geometry_id_check - - call mp_read_keyword(group,path, 'geometry', geometry_id_check) - if(geometry_id_check.ne.geo%id) then - call error("read_grid", "geometry IDs do not match") - end if - call mp_read_array(group, path, array) - - if(any(is_nan(array))) call error("read_grid_5d", "NaN values in 5D array") - - end subroutine read_grid_5d_int8 - - subroutine read_grid_4d_int8(group, path, array, geo) implicit none @@ -182,28 +151,6 @@ subroutine write_grid_3d_int8(group, path, array, geo) end subroutine write_grid_3d_int8 - subroutine read_grid_5d_int(group, path, array, geo) - - implicit none - - integer(hid_t), intent(in) :: group - character(len=*), intent(in) :: path - integer, intent(out) :: array(:,:,:) - type(grid_geometry_desc),intent(in) :: geo - - character(len=32) :: geometry_id_check - - call mp_read_keyword(group,path, 'geometry', geometry_id_check) - if(geometry_id_check.ne.geo%id) then - call error("read_grid", "geometry IDs do not match") - end if - call mp_read_array(group, path, array) - - if(any(is_nan(array))) call error("read_grid_5d", "NaN values in 5D array") - - end subroutine read_grid_5d_int - - subroutine read_grid_4d_int(group, path, array, geo) implicit none @@ -294,28 +241,6 @@ subroutine write_grid_3d_int(group, path, array, geo) end subroutine write_grid_3d_int - subroutine read_grid_5d_dp(group, path, array, geo) - - implicit none - - integer(hid_t), intent(in) :: group - character(len=*), intent(in) :: path - real(dp), intent(out) :: array(:,:,:) - type(grid_geometry_desc),intent(in) :: geo - - character(len=32) :: geometry_id_check - - call mp_read_keyword(group,path, 'geometry', geometry_id_check) - if(geometry_id_check.ne.geo%id) then - call error("read_grid", "geometry IDs do not match") - end if - call mp_read_array(group, path, array) - - if(any(is_nan(array))) call error("read_grid_5d", "NaN values in 5D array") - - end subroutine read_grid_5d_dp - - subroutine read_grid_4d_dp(group, path, array, geo) implicit none @@ -404,29 +329,6 @@ subroutine write_grid_3d_dp(group, path, array, geo) end subroutine write_grid_3d_dp - - subroutine read_grid_5d_sp(group, path, array, geo) - - implicit none - - integer(hid_t), intent(in) :: group - character(len=*), intent(in) :: path - real(sp), intent(out) :: array(:,:,:) - type(grid_geometry_desc),intent(in) :: geo - - character(len=32) :: geometry_id_check - - call mp_read_keyword(group,path, 'geometry', geometry_id_check) - if(geometry_id_check.ne.geo%id) then - call error("read_grid", "geometry IDs do not match") - end if - call mp_read_array(group, path, array) - - if(any(is_nan(array))) call error("read_grid_5d", "NaN values in 5D array") - - end subroutine read_grid_5d_sp - - subroutine read_grid_4d_sp(group, path, array, geo) implicit none diff --git a/src/grid/grid_io_amr.f90 b/src/grid/grid_io_amr.f90 index 0d3e3135..b49281c4 100644 --- a/src/grid/grid_io_amr.f90 +++ b/src/grid/grid_io_amr.f90 @@ -31,21 +31,6 @@ module grid_io module procedure read_grid_4d_int8 end interface read_grid_4d -! The 5d grid reading is commented out currently becuase the compilation crashes in read_grid_5d, owing to some issue in the where statement. we get an error: - -!src/grid/grid_io_amr.f90(113): error #6783: The variable being defined does not conform with the mask-expr of the where-construct or where-stmt. [ARRAY5D] -! array5d(:,:,:,:,:) = 0 - -!i think this could be due to something that needs to be updated in type_grid_amr.f90 - -! interface read_grid_5d -! module procedure read_grid_5d_sp -! module procedure read_grid_5d_dp -! module procedure read_grid_5d_int -! module procedure read_grid_5d_int8 -! end interface read_grid_5d - - interface write_grid_3d module procedure write_grid_3d_sp module procedure write_grid_3d_dp @@ -93,47 +78,6 @@ logical function grid_exists(group, name) end function grid_exists -! subroutine read_grid_5d_int8(group, path, array, geo) - -! implicit none - -! integer(hid_t), intent(in) :: group -! character(len=*), intent(in) :: path -! integer(idp), intent(out) :: array(:,:,:) -! type(grid_geometry_desc),intent(in),target :: geo -! integer(idp), allocatable :: array5d(:,:,:,:,:) -! character(len=100) :: full_path -! integer :: ilevel, igrid, idust -! type(level_desc), pointer :: level -! type(grid_desc), pointer :: grid - -! do ilevel=1,size(geo%levels) -! level => geo%levels(ilevel) -! do igrid=1,size(level%grids) -! grid => level%grids(igrid) -! write(full_path, '("level_", I5.5, "/grid_", I5.5,"/")') ilevel, igrid -! full_path = trim(full_path)//trim(path) -! call mp_read_array_auto(group, full_path, array5d) -! if(any(is_nan(array5d))) call error("read_grid_5d", "NaN values in 5D array") -! do idust=1,size(array5d, 5) -! where(grid%goto_grid(1:grid%n1,1:grid%n2,1:grid%n3) > 0) -! array5d(:,:,:,:,idust) = 0 -! end where -! end do -! array(grid%start_id:grid%start_id + grid%n_cells - 1, :, :) = reshape(array5d, (/grid%n_cells, size(array, 2), size(array,3)/)) -! end do -! end do - -! ! The following three lines provide a workaround for the PGI Fortran -! ! compiler, which otherwise crashes with the following error: -! ! 0: RESHAPE: result type != SOURCE type -! contains -! subroutine test() -! end subroutine test - -! end subroutine read_grid_5d_int8 - - subroutine read_grid_4d_int8(group, path, array, geo) implicit none @@ -339,47 +283,6 @@ end subroutine write_grid_3d_int8 -! subroutine read_grid_5d_int(group, path, array, geo) - -! implicit none - -! integer(hid_t), intent(in) :: group -! character(len=*), intent(in) :: path -! integer, intent(out) :: array(:,:,:) -! type(grid_geometry_desc),intent(in),target :: geo -! integer, allocatable :: array5d(:,:,:,:,:) -! character(len=100) :: full_path -! integer :: ilevel, igrid, idust -! type(level_desc), pointer :: level -! type(grid_desc), pointer :: grid - -! do ilevel=1,size(geo%levels) -! level => geo%levels(ilevel) -! do igrid=1,size(level%grids) -! grid => level%grids(igrid) -! write(full_path, '("level_", I5.5, "/grid_", I5.5,"/")') ilevel, igrid -! full_path = trim(full_path)//trim(path) -! call mp_read_array_auto(group, full_path, array5d) -! if(any(is_nan(array5d))) call error("read_grid_5d", "NaN values in 5D array") -! do idust=1,size(array5d, 5) -! where(grid%goto_grid(1:grid%n1,1:grid%n2,1:grid%n3) > 0) -! array5d(:,:,:,:,idust) = 0 -! end where -! end do -! array(grid%start_id:grid%start_id + grid%n_cells - 1, :, :) = reshape(array5d, (/grid%n_cells, size(array, 2), size(array,3)/)) -! end do -! end do - -! ! The following three lines provide a workaround for the PGI Fortran -! ! compiler, which otherwise crashes with the following error: -! ! 0: RESHAPE: result type != SOURCE type -! contains -! subroutine test() -! end subroutine test - -! end subroutine read_grid_5d_int - - subroutine read_grid_4d_int(group, path, array, geo) implicit none @@ -579,46 +482,6 @@ subroutine write_grid_3d_int(group, path, array, geo) end subroutine write_grid_3d_int -! subroutine read_grid_5d_dp(group, path, array, geo) - -! implicit none - -! integer(hid_t), intent(in) :: group -! character(len=*), intent(in) :: path -! real(dp), intent(out) :: array(:,:,:) -! type(grid_geometry_desc),intent(in),target :: geo -! real(dp), allocatable :: array5d(:,:,:,:,:) -! character(len=100) :: full_path -! integer :: ilevel, igrid, idust -! type(level_desc), pointer :: level -! type(grid_desc), pointer :: grid - -! do ilevel=1,size(geo%levels) -! level => geo%levels(ilevel) -! do igrid=1,size(level%grids) -! grid => level%grids(igrid) -! write(full_path, '("level_", I5.5, "/grid_", I5.5,"/")') ilevel, igrid -! full_path = trim(full_path)//trim(path) -! call mp_read_array_auto(group, full_path, array5d) -! if(any(is_nan(array5d))) call error("read_grid_5d", "NaN values in 5D array") -! do idust=1,size(array5d, 5) -! where(grid%goto_grid(1:grid%n1,1:grid%n2,1:grid%n3) > 0) -! array5d(:,:,:,:,idust) = 0 -! end where -! end do -! array(grid%start_id:grid%start_id + grid%n_cells - 1, :, :) = reshape(array5d, (/grid%n_cells, size(array, 2), size(array,3)/)) -! end do -! end do - -! ! The following three lines provide a workaround for the PGI Fortran -! ! compiler, which otherwise crashes with the following error: -! ! 0: RESHAPE: result type != SOURCE type -! contains -! subroutine test() -! end subroutine test - -! end subroutine read_grid_5d_dp - subroutine read_grid_4d_dp(group, path, array, geo) implicit none @@ -818,45 +681,6 @@ subroutine write_grid_3d_dp(group, path, array, geo) end subroutine write_grid_3d_dp -! subroutine read_grid_5d_sp(group, path, array, geo) - -! implicit none - -! integer(hid_t), intent(in) :: group -! character(len=*), intent(in) :: path -! real(sp), intent(out) :: array(:,:,:) -! type(grid_geometry_desc),intent(in),target :: geo -! real(sp), allocatable :: array5d(:,:,:,:,:) -! character(len=100) :: full_path -! integer :: ilevel, igrid, idust -! type(level_desc), pointer :: level -! type(grid_desc), pointer :: grid - -! do ilevel=1,size(geo%levels) -! level => geo%levels(ilevel) -! do igrid=1,size(level%grids) -! grid => level%grids(igrid) -! write(full_path, '("level_", I5.5, "/grid_", I5.5,"/")') ilevel, igrid -! full_path = trim(full_path)//trim(path) -! call mp_read_array_auto(group, full_path, array5d) -! if(any(is_nan(array5d))) call error("read_grid_5d", "NaN values in 5D array") -! do idust=1,size(array5d, 5) -! where(grid%goto_grid(1:grid%n1,1:grid%n2,1:grid%n3) > 0) -! array5d(:,:,:,:,idust) = 0 -! end where -! end do -! array(grid%start_id:grid%start_id + grid%n_cells - 1, :,:) = reshape(array5d, (/grid%n_cells, size(array, 2),size(array,3)/)) -! end do -! end do - -! ! The following three lines provide a workaround for the PGI Fortran -! ! compiler, which otherwise crashes with the following error: -! ! 0: RESHAPE: result type != SOURCE type -! contains -! subroutine test() -! end subroutine test - -! end subroutine read_grid_5d_sp subroutine read_grid_4d_sp(group, path, array, geo) diff --git a/src/grid/grid_io_template.f90 b/src/grid/grid_io_template.f90 index 439ca0d4..696575d0 100644 --- a/src/grid/grid_io_template.f90 +++ b/src/grid/grid_io_template.f90 @@ -11,7 +11,6 @@ module grid_io public :: grid_exists public :: read_grid_3d public :: read_grid_4d - public :: read_grid_5d public :: write_grid_3d public :: write_grid_4d public :: write_grid_5d @@ -31,13 +30,6 @@ module grid_io module procedure read_grid_4d_int8 end interface read_grid_4d - interface read_grid_5d - module procedure read_grid_5d_sp - module procedure read_grid_5d_dp - module procedure read_grid_5d_int - module procedure read_grid_5d_int8 - end interface read_grid_5d - interface write_grid_3d module procedure write_grid_3d_sp module procedure write_grid_3d_dp @@ -72,37 +64,6 @@ end function grid_exists !!@FOR real(sp):sp real(dp):dp integer:int integer(idp):int8 - subroutine read_grid_5d_(group, path, array, geo) - - implicit none - - integer(hid_t), intent(in) :: group - character(len=*), intent(in) :: path - @T, intent(out) :: array(:,:,:) - type(grid_geometry_desc),intent(in) :: geo - @T, allocatable :: array5d(:,:,:,:,:) - integer :: n_cells, n_dust, n_isrf_lam - - character(len=32) :: geometry_id_check - - call mp_read_keyword(group,path, 'geometry', geometry_id_check) - if(geometry_id_check.ne.geo%id) then - call error("read_grid", "geometry IDs do not match") - end if - call mp_read_array_auto(group,path, array5d) - - if(any(is_nan(array5d))) call error("read_grid_5d", "NaN values in 5D array") - - n_cells = size(array, 1) - n_dust = size(array, 2) - n_isrf_lam = size(array,3) - - array = reshape(array5d, (/n_cells, n_dust, n_isrf_lam/)) - - end subroutine read_grid_5d_ - - - subroutine read_grid_4d_(group, path, array, geo) implicit none diff --git a/src/grid/grid_physics_3d.f90 b/src/grid/grid_physics_3d.f90 index 1d8c804e..a9a28fdf 100644 --- a/src/grid/grid_physics_3d.f90 +++ b/src/grid/grid_physics_3d.f90 @@ -3,7 +3,7 @@ module grid_physics use core_lib use type_photon use type_grid_cell - use grid_io, only : read_grid_4d, grid_exists, read_grid_5d + use grid_io, only : read_grid_4d, grid_exists use dust_main ! many variables and routines use type_dust use grid_geometry From e1d312b5cead0bbb77ea63ba8b4635bf090ac70c Mon Sep 17 00:00:00 2001 From: dnarayanan Date: Tue, 24 May 2022 12:42:56 -0400 Subject: [PATCH 13/38] one more 5d file turned off --- src/grid/grid_io_1d_template.f90 | 30 ------------------------------ 1 file changed, 30 deletions(-) diff --git a/src/grid/grid_io_1d_template.f90 b/src/grid/grid_io_1d_template.f90 index f34ceec1..f36aa0b7 100644 --- a/src/grid/grid_io_1d_template.f90 +++ b/src/grid/grid_io_1d_template.f90 @@ -11,7 +11,6 @@ module grid_io public :: grid_exists public :: read_grid_3d public :: read_grid_4d - public :: read_grid_5d public :: write_grid_3d public :: write_grid_4d @@ -29,13 +28,6 @@ module grid_io module procedure read_grid_4d_int8 end interface read_grid_4d - interface read_grid_5d - module procedure read_grid_5d_sp - module procedure read_grid_5d_dp - module procedure read_grid_5d_int - module procedure read_grid_5d_int8 - end interface read_grid_5d - interface write_grid_3d module procedure write_grid_3d_sp module procedure write_grid_3d_dp @@ -71,28 +63,6 @@ end function grid_exists !!@FOR real(sp):sp real(dp):dp integer:int integer(idp):int8 - subroutine read_grid_5d_(group, path, array, geo) - - implicit none - - integer(hid_t), intent(in) :: group - character(len=*), intent(in) :: path - @T, intent(out) :: array(:,:,:) - type(grid_geometry_desc),intent(in) :: geo - - character(len=32) :: geometry_id_check - - call mp_read_keyword(group,path, 'geometry', geometry_id_check) - if(geometry_id_check.ne.geo%id) then - call error("read_grid", "geometry IDs do not match") - end if - call mp_read_array(group, path, array) - - if(any(is_nan(array))) call error("read_grid_5d", "NaN values in 5D array") - - end subroutine read_grid_5d_ - - subroutine read_grid_4d_(group, path, array, geo) implicit none From 7642fe84dd3f1f2c4b0f06e86d82ab53b87127ba Mon Sep 17 00:00:00 2001 From: dnarayanan Date: Thu, 11 Aug 2022 11:49:45 -0400 Subject: [PATCH 14/38] critical bug fix changing indexes starting from 0 to 1 [which caused, naturally, random segfaults] --- src/grid/grid_propagate_3d.f90 | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/src/grid/grid_propagate_3d.f90 b/src/grid/grid_propagate_3d.f90 index 232a6ecd..1f10ddcb 100644 --- a/src/grid/grid_propagate_3d.f90 +++ b/src/grid/grid_propagate_3d.f90 @@ -249,11 +249,11 @@ subroutine grid_integrate_noenergy(p,tau_required,tau_achieved) integer :: id integer :: idx - real, dimension(d(0)%n_nu) :: energy_frequency_bins + real, dimension(d(1)%n_nu) :: energy_frequency_bins - do id=1,d(0)%n_nu - energy_frequency_bins(id) = d(0)%nu(id) - print *,'[grid_propagate_3d last iteration] energy_frequency_bins(id) = ',energy_frequency_bins(id) + do id=1,d(1)%n_nu + energy_frequency_bins(id) = d(1)%nu(id) + !print *,'[grid_propagate_3d last iteration] energy_frequency_bins(id) = ',energy_frequency_bins(id) end do radial = (p%r .dot. p%v) > 0. From 65a48cb0fb2804a9cf923d722e84c2b4dab8a758 Mon Sep 17 00:00:00 2001 From: Thomas Robitaille Date: Wed, 10 Jun 2026 12:59:02 +0000 Subject: [PATCH 15/38] Always accumulate specific_energy_sum so non-ISRF runs are heated correctly --- src/grid/grid_propagate_3d.f90 | 28 ++++++++++------------------ 1 file changed, 10 insertions(+), 18 deletions(-) diff --git a/src/grid/grid_propagate_3d.f90 b/src/grid/grid_propagate_3d.f90 index 1f10ddcb..805b2104 100644 --- a/src/grid/grid_propagate_3d.f90 +++ b/src/grid/grid_propagate_3d.f90 @@ -142,26 +142,18 @@ subroutine grid_integrate(p,tau_required,tau_achieved) tau_achieved = tau_achieved + tau_cell - ! Compute the ISRF - if (compute_isrf) then + if (compute_isrf) idx = minloc(abs(energy_frequency_bins - p%nu), DIM=1) - idx = minloc(abs(energy_frequency_bins-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) = & - & specific_energy_sum(p%icell%ic, id) + tmin * p%current_kappa(id) * p%energy - end if - - - if(density(p%icell%ic,id) > 0._dp) then - specific_energy_sum_nu(p%icell%ic,id,idx) = & - & specific_energy_sum_nu(p%icell%ic,id,idx) + tmin * p%current_kappa(id) * p%energy + do id=1,n_dust + if(density(p%icell%ic, id) > 0._dp) then + specific_energy_sum(p%icell%ic, id) = & + & specific_energy_sum(p%icell%ic, id) + tmin * p%current_kappa(id) * p%energy + if (compute_isrf) then + specific_energy_sum_nu(p%icell%ic, id, idx) = & + & specific_energy_sum_nu(p%icell%ic, id, idx) + tmin * p%current_kappa(id) * p%energy end if - end do - endif + end if + end do From d4186f3fe9c01607df7936e4fbe6e83aa14924cc Mon Sep 17 00:00:00 2001 From: Thomas Robitaille Date: Wed, 10 Jun 2026 13:02:27 +0000 Subject: [PATCH 16/38] Guard ISRF accumulation behind compute_isrf and a dust check so dustless models no longer segfault --- src/grid/grid_propagate_3d.f90 | 39 +++++++++++++++++----------------- 1 file changed, 19 insertions(+), 20 deletions(-) diff --git a/src/grid/grid_propagate_3d.f90 b/src/grid/grid_propagate_3d.f90 index 805b2104..c4cf0a4c 100644 --- a/src/grid/grid_propagate_3d.f90 +++ b/src/grid/grid_propagate_3d.f90 @@ -60,13 +60,14 @@ subroutine grid_integrate(p,tau_required,tau_achieved) integer :: source_id integer :: idx - real, dimension(d(1)%n_nu) :: energy_frequency_bins + real, allocatable :: energy_frequency_bins(:) - - do id=1,d(1)%n_nu - energy_frequency_bins(id) = d(1)%nu(id) - end do - + if (compute_isrf .and. n_dust > 0) then + allocate(energy_frequency_bins(d(1)%n_nu)) + do id=1,d(1)%n_nu + energy_frequency_bins(id) = d(1)%nu(id) + end do + end if radial = (p%r .dot. p%v) > 0. @@ -241,12 +242,14 @@ subroutine grid_integrate_noenergy(p,tau_required,tau_achieved) integer :: id integer :: idx - real, dimension(d(1)%n_nu) :: energy_frequency_bins + real, allocatable :: energy_frequency_bins(:) - do id=1,d(1)%n_nu - energy_frequency_bins(id) = d(1)%nu(id) - !print *,'[grid_propagate_3d last iteration] energy_frequency_bins(id) = ',energy_frequency_bins(id) - end do + if (compute_isrf .and. n_dust > 0) then + allocate(energy_frequency_bins(d(1)%n_nu)) + do id=1,d(1)%n_nu + energy_frequency_bins(id) = d(1)%nu(id) + end do + end if radial = (p%r .dot. p%v) > 0. @@ -296,19 +299,15 @@ subroutine grid_integrate_noenergy(p,tau_required,tau_achieved) chi_rho_total = 0._dp - - !Compute the ISRF - - !Figure out what frequency bin in the ISRF calculation the current photon's frequency is closest to - idx = minloc(abs(energy_frequency_bins-p%nu),DIM=1) + ! Figure out which ISRF frequency bin the current photon is closest to + if (compute_isrf) idx = minloc(abs(energy_frequency_bins - p%nu), DIM=1) do id=1,n_dust chi_rho_total = chi_rho_total + p%current_chi(id) * density(p%icell%ic, id) - if(density(p%icell%ic,id) > 0._dp) then - specific_energy_sum_nu(p%icell%ic,id,idx) = & - & specific_energy_sum_nu(p%icell%ic,id,idx) + tmin * p%current_kappa(id) * p%energy - !print *,'[grid_propage_3d last iteration] specific_energy_sum_nu=',specific_energy_sum_nu(p%icell%ic,id,idx) + if (compute_isrf .and. density(p%icell%ic, id) > 0._dp) then + specific_energy_sum_nu(p%icell%ic, id, idx) = & + & specific_energy_sum_nu(p%icell%ic, id, idx) + tmin * p%current_kappa(id) * p%energy end if end do From cda3ea6cd6292941cfbecc4acc1b6fa8bde6e048 Mon Sep 17 00:00:00 2001 From: Thomas Robitaille Date: Wed, 10 Jun 2026 13:06:11 +0000 Subject: [PATCH 17/38] Write ISRF frequency bins as a 1D dataset instead of a grid array to avoid an out-of-bounds segfault --- src/grid/grid_generic.f90 | 73 +++++++++++++++------------------------ 1 file changed, 27 insertions(+), 46 deletions(-) diff --git a/src/grid/grid_generic.f90 b/src/grid/grid_generic.f90 index 2dc2457c..ecb716dd 100644 --- a/src/grid/grid_generic.f90 +++ b/src/grid/grid_generic.f90 @@ -1,7 +1,7 @@ module grid_generic use core_lib, only : sp, dp, hid_t, warn, error - use mpi_hdf5_io, only : mp_create_group + use mpi_hdf5_io, only : mp_create_group, mp_write_array use mpi_core, only : main_process use grid_io, only : write_grid_3d, write_grid_4d, write_grid_5d @@ -35,9 +35,8 @@ subroutine output_grid(group, iter, n_iter) integer(hid_t),intent(in) :: group integer,intent(in) :: iter, n_iter - integer :: n_cells, n_dust, n_isrf_lam - integer :: i,j,k - real, dimension(d(1)%n_nu) :: energy_frequency_bins + integer :: i + real, allocatable :: energy_frequency_bins(:) if(main_process()) write(*,'(" [output_grid] outputting grid arrays for iteration")') @@ -72,52 +71,34 @@ subroutine output_grid(group, iter, n_iter) ! WRITE THE ISRF - if (compute_isrf) then - - n_cells = size(specific_energy_nu, 1) - n_dust = size(specific_energy_nu, 2) - n_isrf_lam = size(specific_energy_nu,3) - - do i=1,d(1)%n_nu - energy_frequency_bins(i) = d(1)%nu(i) - end do - - if(trim(output_specific_energy)=='all' .or. (trim(output_specific_energy)=='last'.and.iter==n_iter)) then - !if(allocated(energy_frequency_bins)) then - call write_grid_3d(group, 'ISRF_frequency_bins',energy_frequency_bins, geo) - !else - !call warn("output_grid","energy_frequency bins [ISRF wavelengths] array is not allocated") - !end if - end if - + if (compute_isrf) then + if(trim(output_specific_energy)=='all' .or. (trim(output_specific_energy)=='last'.and.iter==n_iter)) then + + ! ISRF frequency bins are a per-frequency quantity, not per-cell, so + ! they are written as a plain 1-D dataset rather than a grid array. + allocate(energy_frequency_bins(d(1)%n_nu)) + do i=1,d(1)%n_nu + energy_frequency_bins(i) = d(1)%nu(i) + end do + call mp_write_array(group, 'ISRF_frequency_bins', energy_frequency_bins) + deallocate(energy_frequency_bins) + if(allocated(specific_energy_nu)) then - - - if (compute_isrf) then - print *,'[GRID_GENERIC.F90] ENTERING BLOCK FOR WRITING SPECIFIC_ENERGY_NU' - - select case(physics_io_type) - - case(sp) - print *,'[GRID_GENERIC.F90] diving into write_grid_5d sp' - call write_grid_5d(group, 'specific_energy_nu', real(specific_energy_nu, sp), geo) - - case(dp) - print *,'[GRID_GENERIC.F90] diving into write_grid_5d dp' - call write_grid_5d(group, 'specific_energy_nu', real(specific_energy_nu, dp), geo) - print *,'[GRID_GENERIC.F90] finished with write_grid_5d dp' - - case default - call error("output_grid","unexpected value of physics_io_type (should be sp or dp)") - end select - else - call warn("output_grid","specific_energy_nu array is not allocated") - - end if + select case(physics_io_type) + case(sp) + call write_grid_5d(group, 'specific_energy_nu', real(specific_energy_nu, sp), geo) + case(dp) + call write_grid_5d(group, 'specific_energy_nu', real(specific_energy_nu, dp), geo) + case default + call error("output_grid","unexpected value of physics_io_type (should be sp or dp)") + end select + else + call warn("output_grid","specific_energy_nu array is not allocated") end if + end if - + endif From f73ade86b1381f2f960af296f1cf7c76f45c7b3a Mon Sep 17 00:00:00 2001 From: Thomas Robitaille Date: Wed, 10 Jun 2026 13:10:58 +0000 Subject: [PATCH 18/38] Fix typo in ISRF config reader method name so models can be read back --- hyperion/conf/conf_files.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/hyperion/conf/conf_files.py b/hyperion/conf/conf_files.py index e2a84faa..c0c240a6 100644 --- a/hyperion/conf/conf_files.py +++ b/hyperion/conf/conf_files.py @@ -743,7 +743,7 @@ def read_run_conf(self, group): # not a class method because inherited self._read_max_reabsorptions(group) self._read_pda(group) self._read_mrw(group) - self._reada_isrf(group) + self._read_isrf(group) self._read_convergence(group) self._read_kill_on_absorb(group) self._read_kill_on_scatter(group) From 2332847957adade303f7ff1d998c7d8cbd46c4e2 Mon Sep 17 00:00:00 2001 From: Thomas Robitaille Date: Wed, 10 Jun 2026 13:28:25 +0000 Subject: [PATCH 19/38] Cache the ISRF frequency grid at setup and stop accumulating it during the final imaging iteration --- src/grid/grid_physics_3d.f90 | 10 +++++++-- src/grid/grid_propagate_3d.f90 | 37 +++------------------------------- 2 files changed, 11 insertions(+), 36 deletions(-) diff --git a/src/grid/grid_physics_3d.f90 b/src/grid/grid_physics_3d.f90 index a9a28fdf..913aafcf 100644 --- a/src/grid/grid_physics_3d.f90 +++ b/src/grid/grid_physics_3d.f90 @@ -39,6 +39,7 @@ module grid_physics real(dp),allocatable, public :: specific_energy_nu(:,:,:) real(dp),allocatable, public :: specific_energy_sum(:,:) real(dp),allocatable, public :: specific_energy_sum_nu(:,:,:) + real(dp),allocatable, public :: isrf_nu(:) real(dp),allocatable, public :: specific_energy_additional(:,:) real(dp),allocatable, public :: specific_energy_additional_nu(:,:,:) @@ -232,10 +233,15 @@ subroutine setup_grid_physics(group, use_mrw, use_pda, compute_isrf) specific_energy_sum = 0._dp ! Set up basics for ISRF calculation - !n_isrf_wavelengths = d(1)%n_nu allocate(specific_energy_sum_nu(geo%n_cells, n_dust, n_isrf_wavelengths)) specific_energy_sum_nu = 0._dp - + + ! Cache the ISRF frequency grid once so it does not have to be rebuilt for every photon + allocate(isrf_nu(n_isrf_wavelengths)) + do idx=1,n_isrf_wavelengths + isrf_nu(idx) = d(1)%nu(idx) + end do + ! Total energy absorbed allocate(energy_abs_tot(n_dust)) energy_abs_tot = 0._dp diff --git a/src/grid/grid_propagate_3d.f90 b/src/grid/grid_propagate_3d.f90 index c4cf0a4c..5aa7bb7d 100644 --- a/src/grid/grid_propagate_3d.f90 +++ b/src/grid/grid_propagate_3d.f90 @@ -3,15 +3,13 @@ module grid_propagate use core_lib use type_photon, only : photon use type_grid_cell - use dust_main, only : n_dust,d + use dust_main, only : n_dust use grid_geometry, only : escaped, find_wall, in_correct_cell, next_cell, opposite_wall - use grid_physics, only : specific_energy_sum, specific_energy_sum_nu, density, n_photons, last_photon_id + use grid_physics, only : specific_energy_sum, specific_energy_sum_nu, isrf_nu, density, n_photons, last_photon_id use sources use counters use settings, only : frac_check => propagation_check_frequency, compute_isrf - use grid_geometry, only : geo - implicit none save @@ -60,14 +58,6 @@ subroutine grid_integrate(p,tau_required,tau_achieved) integer :: source_id integer :: idx - real, allocatable :: energy_frequency_bins(:) - - if (compute_isrf .and. n_dust > 0) then - allocate(energy_frequency_bins(d(1)%n_nu)) - do id=1,d(1)%n_nu - energy_frequency_bins(id) = d(1)%nu(id) - end do - end if radial = (p%r .dot. p%v) > 0. @@ -143,7 +133,7 @@ subroutine grid_integrate(p,tau_required,tau_achieved) tau_achieved = tau_achieved + tau_cell - if (compute_isrf) idx = minloc(abs(energy_frequency_bins - p%nu), DIM=1) + if (compute_isrf) idx = minloc(abs(isrf_nu - p%nu), DIM=1) do id=1,n_dust if(density(p%icell%ic, id) > 0._dp) then @@ -240,17 +230,6 @@ subroutine grid_integrate_noenergy(p,tau_required,tau_achieved) integer :: source_id - integer :: id - integer :: idx - real, allocatable :: energy_frequency_bins(:) - - if (compute_isrf .and. n_dust > 0) then - allocate(energy_frequency_bins(d(1)%n_nu)) - do id=1,d(1)%n_nu - energy_frequency_bins(id) = d(1)%nu(id) - end do - end if - radial = (p%r .dot. p%v) > 0. if(debug) write(*,'(" [debug] start grid_integrate_noenergy")') @@ -298,18 +277,8 @@ subroutine grid_integrate_noenergy(p,tau_required,tau_achieved) end if chi_rho_total = 0._dp - - ! Figure out which ISRF frequency bin the current photon is closest to - if (compute_isrf) idx = minloc(abs(energy_frequency_bins - p%nu), DIM=1) - do id=1,n_dust chi_rho_total = chi_rho_total + p%current_chi(id) * density(p%icell%ic, id) - - if (compute_isrf .and. density(p%icell%ic, id) > 0._dp) then - specific_energy_sum_nu(p%icell%ic, id, idx) = & - & specific_energy_sum_nu(p%icell%ic, id, idx) + tmin * p%current_kappa(id) * p%energy - end if - end do From 7109c0680518b26e080ecc3a36d902f8040dd34f Mon Sep 17 00:00:00 2001 From: Thomas Robitaille Date: Wed, 10 Jun 2026 13:28:25 +0000 Subject: [PATCH 20/38] Remove a leftover comment and extra blank lines from the ISRF work --- src/grid/grid_io.f90 | 1 - src/mpi/mpi_routines.f90 | 2 -- 2 files changed, 3 deletions(-) diff --git a/src/grid/grid_io.f90 b/src/grid/grid_io.f90 index 0cba3f1d..18c4e343 100644 --- a/src/grid/grid_io.f90 +++ b/src/grid/grid_io.f90 @@ -428,7 +428,6 @@ subroutine read_grid_3d_sp(group, path, array, geo) end subroutine read_grid_3d_sp - !CRAZY DN ADDITIONS subroutine write_grid_5d_sp(group, path, array, geo) implicit none diff --git a/src/mpi/mpi_routines.f90 b/src/mpi/mpi_routines.f90 index 32e1c777..6adbf7e0 100644 --- a/src/mpi/mpi_routines.f90 +++ b/src/mpi/mpi_routines.f90 @@ -289,8 +289,6 @@ subroutine mp_collect_physical_arrays() call mpi_reduce(specific_energy_sum_nu, dummy_dp, size(specific_energy_sum_nu), mpi_real8, mpi_sum, rank_main, mpi_comm_world, ierr) end if - - if(allocated(n_photons)) then if(main_process()) then allocate(tmp_int_1d(size(n_photons,1))) From 364550056409134a0494ec3e57ade81cc9edc378 Mon Sep 17 00:00:00 2001 From: Thomas Robitaille Date: Wed, 10 Jun 2026 13:38:16 +0000 Subject: [PATCH 21/38] Only populate the cached ISRF frequency grid when ISRF is enabled so dustless models do not segfault --- src/grid/grid_physics_3d.f90 | 10 ++++++---- 1 file changed, 6 insertions(+), 4 deletions(-) diff --git a/src/grid/grid_physics_3d.f90 b/src/grid/grid_physics_3d.f90 index 913aafcf..1389b225 100644 --- a/src/grid/grid_physics_3d.f90 +++ b/src/grid/grid_physics_3d.f90 @@ -237,10 +237,12 @@ subroutine setup_grid_physics(group, use_mrw, use_pda, compute_isrf) specific_energy_sum_nu = 0._dp ! Cache the ISRF frequency grid once so it does not have to be rebuilt for every photon - allocate(isrf_nu(n_isrf_wavelengths)) - do idx=1,n_isrf_wavelengths - isrf_nu(idx) = d(1)%nu(idx) - end do + if (compute_isrf) then + allocate(isrf_nu(n_isrf_wavelengths)) + do idx=1,n_isrf_wavelengths + isrf_nu(idx) = d(1)%nu(idx) + end do + end if ! Total energy absorbed allocate(energy_abs_tot(n_dust)) From b020ca7848e1d33f8623d79ac8925b9e3bdfdb01 Mon Sep 17 00:00:00 2001 From: Thomas Robitaille Date: Wed, 10 Jun 2026 13:58:52 +0000 Subject: [PATCH 22/38] Fix the AMR write_grid_5d template slice so a regenerated grid_io_amr matches the working code --- src/grid/grid_io_amr_template.f90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/grid/grid_io_amr_template.f90 b/src/grid/grid_io_amr_template.f90 index 5e06a8c1..a37d6354 100644 --- a/src/grid/grid_io_amr_template.f90 +++ b/src/grid/grid_io_amr_template.f90 @@ -189,7 +189,7 @@ subroutine write_grid_5d_(group, path, array, geo) else g_grid = mp_create_group(g_level, name) end if - call mp_write_array(g_grid, path, reshape(array(grid%start_id:grid%start_id + grid%n_cells - 1, :), & + call mp_write_array(g_grid, path, reshape(array(grid%start_id:grid%start_id + grid%n_cells - 1, :, :), & & (/grid%n1, grid%n2, grid%n3, size(array,2), size(array,3)/))) call mp_close_group(g_grid) end do From f10e30ec5ebf7e54b66be444c6c2c5514a746231 Mon Sep 17 00:00:00 2001 From: Thomas Robitaille Date: Wed, 10 Jun 2026 14:08:09 +0000 Subject: [PATCH 23/38] Add regression tests for the ISRF covering heating, energy conservation, frequency bins, dustless models, AMR, and MPI --- hyperion/model/tests/test_isrf.py | 163 ++++++++++++++++++++++++++++++ 1 file changed, 163 insertions(+) create mode 100644 hyperion/model/tests/test_isrf.py diff --git a/hyperion/model/tests/test_isrf.py b/hyperion/model/tests/test_isrf.py new file mode 100644 index 00000000..ed988920 --- /dev/null +++ b/hyperion/model/tests/test_isrf.py @@ -0,0 +1,163 @@ +from __future__ import print_function, division + +import shutil + +import numpy as np +import pytest + +import h5py + +from .. import Model +from ...grid import AMRGrid +from ...util.functions import random_id +from .test_helpers import get_test_dust + + +def _read_dataset(filename, suffix): + """Return the last grid dataset whose HDF5 path ends with ``suffix``.""" + with h5py.File(filename, 'r') as f: + matches = [] + f.visititems(lambda name, obj: matches.append(name) + if name.endswith(suffix) and isinstance(obj, h5py.Dataset) + else None) + if not matches: + return None + return np.asarray(f[matches[-1]][()], dtype=float) + + +def _assert_isrf_reconstructs_specific_energy(filename): + # Summed over frequency, the ISRF must reproduce the specific energy. Cells + # with no absorption are clamped up to the minimum specific energy (which + # only affects specific_energy, not the unclamped specific_energy_nu), so we + # compare only cells that were genuinely heated above that floor. + se = _read_dataset(filename, '/specific_energy') + se_nu = _read_dataset(filename, '/specific_energy_nu') + assert se_nu is not None + nu_sum = se_nu.sum(axis=0) + floor = se.min() + heated = se > floor * (1. + 1.e-6) + assert np.count_nonzero(heated) >= 1 + np.testing.assert_allclose(nu_sum[heated], se[heated], rtol=1.e-6) + + +def _cartesian_isrf_model(compute_isrf, n_cells=2, n_photons=100000, density=1.e-18): + m = Model() + edges = np.linspace(-1., 1., n_cells + 1) + m.set_cartesian_grid(edges, edges, edges) + dust = get_test_dust() + m.add_density_grid(np.ones((n_cells, n_cells, n_cells)) * density, dust) + s = m.add_point_source() + s.luminosity = 1. + s.temperature = 6000. + m.set_n_initial_iterations(3) + m.set_n_photons(initial=n_photons, imaging=0) + m.set_seed(-12345) + m.conf.output.output_specific_energy = 'last' + m.compute_isrf(compute_isrf) + return m + + +@pytest.mark.requires_hyperion_binaries +def test_isrf_does_not_change_specific_energy(tmpdir): + # Computing the ISRF must not perturb the ordinary specific-energy + # (temperature) calculation. This guards against the regression where the + # specific_energy_sum accumulation was accidentally gated behind compute_isrf. + se = {} + for isrf in (False, True): + m = _cartesian_isrf_model(isrf) + m.write(tmpdir.join(random_id()).strpath) + out = m.run(tmpdir.join(random_id()).strpath) + se[isrf] = _read_dataset(out.filename, '/specific_energy') + # Same seed, and the ISRF code path does not touch the random stream, so + # the specific energy should be identical with and without the ISRF. + np.testing.assert_allclose(se[True], se[False], rtol=1.e-10) + + +@pytest.mark.requires_hyperion_binaries +def test_isrf_sums_to_specific_energy(tmpdir): + # The frequency-resolved ISRF, summed over frequency, must reproduce the + # total specific energy. + m = _cartesian_isrf_model(True) + m.write(tmpdir.join(random_id()).strpath) + out = m.run(tmpdir.join(random_id()).strpath) + _assert_isrf_reconstructs_specific_energy(out.filename) + + +@pytest.mark.requires_hyperion_binaries +def test_isrf_frequency_bins_written(tmpdir): + # ISRF_frequency_bins is a per-frequency (1-D) quantity, not per-cell. This + # guards against the regression where it was written as a grid array, which + # segfaulted whenever the number of frequencies was smaller than n_cells. + m = _cartesian_isrf_model(True, n_cells=8) # 512 cells, dust has 2 frequencies + m.write(tmpdir.join(random_id()).strpath) + out = m.run(tmpdir.join(random_id()).strpath) + bins = _read_dataset(out.filename, '/ISRF_frequency_bins') + assert bins is not None + assert bins.shape == (2,) # == number of dust frequencies, not n_cells + + +@pytest.mark.requires_hyperion_binaries +def test_isrf_dustless_model_runs(tmpdir): + # A model with no dust must still run (and image) without crashing. This + # guards against the regression where ISRF instrumentation dereferenced the + # first dust type unconditionally, segfaulting dustless models. + m = Model() + m.set_cartesian_grid([-1., 1.], [-1., 1.], [-1., 1.]) + s = m.add_point_source() + s.luminosity = 1. + s.temperature = 6000. + i = m.add_peeled_images(sed=True, image=False) + i.set_viewing_angles([45.], [45.]) + i.set_wavelength_range(5, 0.1, 100.) + m.set_n_initial_iterations(0) + m.set_n_photons(imaging=10000) + m.write(tmpdir.join(random_id()).strpath) + # Should complete without raising (previously segfaulted in the final iteration). + m.run(tmpdir.join(random_id()).strpath) + + +@pytest.mark.requires_hyperion_binaries +def test_isrf_amr(tmpdir): + # The ISRF must work on AMR grids: it should run and, summed over frequency, + # reproduce the specific energy in every cell. + m = Model() + 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, 4, 4 + grid.quantities['density'] = [np.ones((4, 4, 4)) * 1.e-16] + m.set_amr_grid(amr) + m.add_density_grid(amr['density'][0], get_test_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(-9) + m.conf.output.output_specific_energy = 'last' + m.compute_isrf(True) + m.write(tmpdir.join(random_id()).strpath) + out = m.run(tmpdir.join(random_id()).strpath) + _assert_isrf_reconstructs_specific_energy(out.filename) + + +@pytest.mark.requires_hyperion_binaries +def test_isrf_mpi_matches_serial(tmpdir): + # The saved ISRF must not depend on the number of MPI processes. We compare + # the total ISRF energy (which is conserved, hence robust to Monte Carlo + # noise) between a serial run and a 2-process MPI run. + if shutil.which('mpirun') is None and shutil.which('mpiexec') is None: + pytest.skip("no MPI launcher available") + + m = _cartesian_isrf_model(True, n_photons=200000) + m.write(tmpdir.join(random_id()).strpath) + + out_serial = m.run(tmpdir.join(random_id()).strpath) + out_mpi = m.run(tmpdir.join(random_id()).strpath, mpi=True, n_processes=2) + + total_serial = np.nansum(_read_dataset(out_serial.filename, '/specific_energy_nu')) + total_mpi = np.nansum(_read_dataset(out_mpi.filename, '/specific_energy_nu')) + np.testing.assert_allclose(total_mpi, total_serial, rtol=2.e-2) From 199093f8be63f36af7b7317bdaad22e68e20fec7 Mon Sep 17 00:00:00 2001 From: Thomas Robitaille Date: Wed, 10 Jun 2026 14:57:04 +0000 Subject: [PATCH 24/38] Document the compute_isrf option and the specific_energy_nu and ISRF_frequency_bins outputs --- docs/setup/setup_conf.rst | 29 +++++++++++++++++++++++++++++ hyperion/conf/conf_files.py | 19 ++++++++++++++----- 2 files changed, 43 insertions(+), 5 deletions(-) diff --git a/docs/setup/setup_conf.rst b/docs/setup/setup_conf.rst index b4b7908f..74992fa2 100644 --- a/docs/setup/setup_conf.rst +++ b/docs/setup/setup_conf.rst @@ -195,6 +195,35 @@ only after last iteration), or ``none`` (do not output). The default is to output only the last iteration of ``specific_energy``. To find out how to view these values, see :doc:`../postprocessing/postprocessing` +Interstellar radiation field +---------------------------- + +By default, Hyperion saves the total specific energy absorbed in each cell. If +you also want the frequency-resolved specific energy absorbed (the interstellar +radiation field, ISRF), you can enable this with:: + + m.compute_isrf(True) + +When enabled, two extra arrays are written out following the same +``output_specific_energy`` setting described above: + +* ``specific_energy_nu`` -- the specific energy absorbed in each cell as a + function of frequency, in erg/s/g. +* ``ISRF_frequency_bins`` -- the frequencies (in Hz) corresponding to the last + axis of ``specific_energy_nu``. + +The frequency grid is taken from the first dust type, so all dust types should +share the same frequency grid. Summed over frequency, ``specific_energy_nu`` +recovers the total ``specific_energy``. This option is disabled by default and +works for all grid types, including AMR. + +``specific_energy_nu`` is not currently returned by ``get_quantities`` and +should be read directly from the output file, e.g.:: + + import h5py + with h5py.File('output.rtout', 'r') as f: + senu = f['iteration_00005/specific_energy_nu'][()] + Advanced parameters ------------------- diff --git a/hyperion/conf/conf_files.py b/hyperion/conf/conf_files.py index c0c240a6..44313a4e 100644 --- a/hyperion/conf/conf_files.py +++ b/hyperion/conf/conf_files.py @@ -402,15 +402,24 @@ def _read_isrf(self,group): def _write_isrf(self,group): group.attrs['isrf'] = bool2str(self.isrf) - def compute_isrf(self,isrf): + def compute_isrf(self, isrf): ''' + Set whether to compute and save the interstellar radiation field (ISRF) + in each cell. - Set whether or not to compute and save the ISRF in each cell - - If enabled, the ISRF is computed in every cell at the - frequencies of the dust opacity tables. + If enabled, the specific energy absorbed is also accumulated as a + function of frequency and saved as the ``specific_energy_nu`` quantity + (in erg/s/g), alongside the ``ISRF_frequency_bins`` array giving the + corresponding frequencies (in Hz). The frequency grid is taken from the + first dust type, so all dust types should share the same frequency grid. + These arrays are written following the ``output_specific_energy`` + setting. This is disabled by default. + Parameters + ---------- + isrf : bool + Whether to compute and save the ISRF in each cell. ''' self.isrf = isrf From 7b83de8dff424ba7b4a59a41440ba96f9205c7e8 Mon Sep 17 00:00:00 2001 From: Thomas Robitaille Date: Wed, 10 Jun 2026 15:59:09 +0000 Subject: [PATCH 25/38] Allow the ISRF frequency grid to be set independently of the dust grid and bin photons in log space --- src/grid/grid_generic.f90 | 15 ++------------- src/grid/grid_physics_3d.f90 | 31 +++++++++++++++++++++++-------- src/grid/grid_propagate_3d.f90 | 4 ++-- src/main/settings.f90 | 4 ++++ src/main/setup_rt.f90 | 7 +++++++ 5 files changed, 38 insertions(+), 23 deletions(-) diff --git a/src/grid/grid_generic.f90 b/src/grid/grid_generic.f90 index ecb716dd..326560bc 100644 --- a/src/grid/grid_generic.f90 +++ b/src/grid/grid_generic.f90 @@ -6,12 +6,9 @@ module grid_generic use grid_io, only : write_grid_3d, write_grid_4d, write_grid_5d use grid_geometry, only : geo - use grid_physics, only : n_photons, last_photon_id, specific_energy_sum, specific_energy_sum_nu, specific_energy, specific_energy_nu, density, density_original + use grid_physics, only : n_photons, last_photon_id, specific_energy_sum, specific_energy_sum_nu, specific_energy, specific_energy_nu, isrf_nu, density, density_original use settings, only : output_n_photons, output_specific_energy, output_density, output_density_diff, physics_io_type, compute_isrf - use dust_main, only: d - - implicit none save @@ -35,9 +32,6 @@ subroutine output_grid(group, iter, n_iter) integer(hid_t),intent(in) :: group integer,intent(in) :: iter, n_iter - integer :: i - real, allocatable :: energy_frequency_bins(:) - if(main_process()) write(*,'(" [output_grid] outputting grid arrays for iteration")') @@ -77,12 +71,7 @@ subroutine output_grid(group, iter, n_iter) ! ISRF frequency bins are a per-frequency quantity, not per-cell, so ! they are written as a plain 1-D dataset rather than a grid array. - allocate(energy_frequency_bins(d(1)%n_nu)) - do i=1,d(1)%n_nu - energy_frequency_bins(i) = d(1)%nu(i) - end do - call mp_write_array(group, 'ISRF_frequency_bins', energy_frequency_bins) - deallocate(energy_frequency_bins) + call mp_write_array(group, 'ISRF_frequency_bins', isrf_nu) if(allocated(specific_energy_nu)) then select case(physics_io_type) diff --git a/src/grid/grid_physics_3d.f90 b/src/grid/grid_physics_3d.f90 index 1389b225..de743a5e 100644 --- a/src/grid/grid_physics_3d.f90 +++ b/src/grid/grid_physics_3d.f90 @@ -40,6 +40,7 @@ module grid_physics real(dp),allocatable, public :: specific_energy_sum(:,:) real(dp),allocatable, public :: specific_energy_sum_nu(:,:,:) real(dp),allocatable, public :: isrf_nu(:) + real(dp),allocatable, public :: log_isrf_nu(:) real(dp),allocatable, public :: specific_energy_additional(:,:) real(dp),allocatable, public :: specific_energy_additional_nu(:,:,:) @@ -106,7 +107,13 @@ subroutine setup_grid_physics(group, use_mrw, use_pda, compute_isrf) logical,intent(in) :: use_mrw, use_pda, compute_isrf integer :: n_isrf_wavelengths - n_isrf_wavelengths = d(1)%n_nu + ! The ISRF is binned onto a user-specified frequency grid if one was given, + ! otherwise onto the frequency grid of the first dust type. + if (allocated(isrf_frequencies)) then + n_isrf_wavelengths = size(isrf_frequencies) + else + n_isrf_wavelengths = d(1)%n_nu + end if ! Density allocate(density(geo%n_cells, n_dust)) allocate(specific_energy(geo%n_cells, n_dust)) @@ -236,12 +243,20 @@ subroutine setup_grid_physics(group, use_mrw, use_pda, compute_isrf) allocate(specific_energy_sum_nu(geo%n_cells, n_dust, n_isrf_wavelengths)) specific_energy_sum_nu = 0._dp - ! Cache the ISRF frequency grid once so it does not have to be rebuilt for every photon + ! Cache the ISRF 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 + ! log-frequency space (see grid_propagate). if (compute_isrf) then allocate(isrf_nu(n_isrf_wavelengths)) - do idx=1,n_isrf_wavelengths - isrf_nu(idx) = d(1)%nu(idx) - end do + if (allocated(isrf_frequencies)) then + isrf_nu = isrf_frequencies + else + do idx=1,n_isrf_wavelengths + isrf_nu(idx) = d(1)%nu(idx) + end do + end if + allocate(log_isrf_nu(n_isrf_wavelengths)) + log_isrf_nu = log10(isrf_nu) end if ! Total energy absorbed @@ -312,7 +327,7 @@ subroutine sublimate_dust() integer :: reset integer :: n_isrf_wavelengths - n_isrf_wavelengths = d(1)%n_nu + n_isrf_wavelengths = size(specific_energy_nu, 3) reset = 0 @@ -398,8 +413,8 @@ subroutine update_energy_abs(scale) integer :: id,idx integer :: n_isrf_wavelengths - - n_isrf_wavelengths = d(1)%n_nu + + n_isrf_wavelengths = size(specific_energy_nu, 3) if(main_process()) write(*,'(" [grid_physics] updating energy_abs")') diff --git a/src/grid/grid_propagate_3d.f90 b/src/grid/grid_propagate_3d.f90 index 5aa7bb7d..0f57ceb2 100644 --- a/src/grid/grid_propagate_3d.f90 +++ b/src/grid/grid_propagate_3d.f90 @@ -5,7 +5,7 @@ module grid_propagate use type_grid_cell use dust_main, only : n_dust use grid_geometry, only : escaped, find_wall, in_correct_cell, next_cell, opposite_wall - use grid_physics, only : specific_energy_sum, specific_energy_sum_nu, isrf_nu, density, n_photons, last_photon_id + use grid_physics, only : specific_energy_sum, specific_energy_sum_nu, log_isrf_nu, density, n_photons, last_photon_id use sources use counters use settings, only : frac_check => propagation_check_frequency, compute_isrf @@ -133,7 +133,7 @@ subroutine grid_integrate(p,tau_required,tau_achieved) tau_achieved = tau_achieved + tau_cell - if (compute_isrf) idx = minloc(abs(isrf_nu - p%nu), DIM=1) + if (compute_isrf) idx = minloc(abs(log_isrf_nu - log10(p%nu)), DIM=1) do id=1,n_dust if(density(p%icell%ic, id) > 0._dp) then diff --git a/src/main/settings.f90 b/src/main/settings.f90 index d8038eef..9e3102d0 100644 --- a/src/main/settings.f90 +++ b/src/main/settings.f90 @@ -30,6 +30,10 @@ module settings real(dp) :: monochromatic_energy_threshold real(dp),allocatable :: frequencies(:) + ! Optional user-specified frequency grid for the ISRF. If not allocated, the + ! frequency grid of the first dust type is used instead. + real(dp),allocatable :: isrf_frequencies(:) + integer :: physics_io_type character(len=4) :: output_density diff --git a/src/main/setup_rt.f90 b/src/main/setup_rt.f90 index 5be59417..462968b0 100644 --- a/src/main/setup_rt.f90 +++ b/src/main/setup_rt.f90 @@ -76,6 +76,13 @@ subroutine setup_initial(input_handle) call mp_read_keyword(input_handle, '/', 'mrw', use_mrw) call mp_read_keyword(input_handle, '/', 'isrf', compute_isrf) + ! Optional user-specified ISRF frequency grid (otherwise the dust grid is used) + if(compute_isrf) then + if(mp_path_exists(input_handle, 'isrf_frequencies')) then + call mp_table_read_column_auto(input_handle, 'isrf_frequencies', 'nu', isrf_frequencies) + end if + end if + if(use_mrw) then call mp_read_keyword(input_handle, '/', 'mrw_gamma', mrw_gamma) call mp_read_keyword(input_handle, '/', 'n_inter_mrw_max', n_mrw_max) From 2b28c2321ba125b58cf878962ea5fee0c22ff796 Mon Sep 17 00:00:00 2001 From: Thomas Robitaille Date: Wed, 10 Jun 2026 15:59:09 +0000 Subject: [PATCH 26/38] Add a frequencies argument to compute_isrf for a custom ISRF frequency grid --- hyperion/conf/conf_files.py | 36 ++++++++++++++++++++++++++++-------- 1 file changed, 28 insertions(+), 8 deletions(-) diff --git a/hyperion/conf/conf_files.py b/hyperion/conf/conf_files.py index 44313a4e..f5d26222 100644 --- a/hyperion/conf/conf_files.py +++ b/hyperion/conf/conf_files.py @@ -396,13 +396,21 @@ def _read_pda(self, group): def _write_pda(self, group): group.attrs['pda'] = bool2str(self.pda) - def _read_isrf(self,group): + def _read_isrf(self, group): self.isrf = str2bool(group.attrs['isrf']) - - def _write_isrf(self,group): + if 'isrf_frequencies' in group: + self.isrf_frequencies = np.array(group['isrf_frequencies']['nu']) + else: + self.isrf_frequencies = None + + def _write_isrf(self, group): group.attrs['isrf'] = bool2str(self.isrf) + if self.isrf and self.isrf_frequencies is not None: + group.create_dataset('isrf_frequencies', + data=np.array(list(zip(self.isrf_frequencies)), + dtype=[('nu', float)])) - def compute_isrf(self, isrf): + def compute_isrf(self, isrf, frequencies=None): ''' Set whether to compute and save the interstellar radiation field (ISRF) @@ -411,17 +419,29 @@ def compute_isrf(self, isrf): If enabled, the specific energy absorbed is also accumulated as a function of frequency and saved as the ``specific_energy_nu`` quantity (in erg/s/g), alongside the ``ISRF_frequency_bins`` array giving the - corresponding frequencies (in Hz). The frequency grid is taken from the - first dust type, so all dust types should share the same frequency grid. - These arrays are written following the ``output_specific_energy`` - setting. This is disabled by default. + corresponding frequencies (in Hz). Photons are binned to the nearest + frequency in log space. These arrays are written following the + ``output_specific_energy`` setting. This is disabled by default. Parameters ---------- isrf : bool Whether to compute and save the ISRF in each cell. + frequencies : iterable of float, optional + The frequencies (in Hz) onto which to bin the ISRF. If not + specified, the frequency grid of the first dust type is used. Only + used when ``isrf`` is `True`. ''' self.isrf = isrf + if frequencies is None: + self.isrf_frequencies = None + else: + frequencies = np.asarray(frequencies, dtype=float) + if frequencies.ndim != 1 or frequencies.size < 1: + raise ValueError("frequencies should be a 1-d array of at least one value") + if np.any(frequencies <= 0.): + raise ValueError("frequencies should be positive (in Hz)") + self.isrf_frequencies = np.sort(frequencies) def set_mrw(self, mrw, gamma=1.0, inter_max=1000, warn=True): From 0aede1508c5d699ba528a0d0e3f55b4f9d577848 Mon Sep 17 00:00:00 2001 From: Thomas Robitaille Date: Wed, 10 Jun 2026 15:59:09 +0000 Subject: [PATCH 27/38] Test and document the configurable ISRF frequency grid --- docs/setup/setup_conf.rst | 14 ++++++--- hyperion/model/tests/test_isrf.py | 47 +++++++++++++++++++++++++++++-- 2 files changed, 55 insertions(+), 6 deletions(-) diff --git a/docs/setup/setup_conf.rst b/docs/setup/setup_conf.rst index 74992fa2..e6cc80a3 100644 --- a/docs/setup/setup_conf.rst +++ b/docs/setup/setup_conf.rst @@ -212,10 +212,16 @@ When enabled, two extra arrays are written out following the same * ``ISRF_frequency_bins`` -- the frequencies (in Hz) corresponding to the last axis of ``specific_energy_nu``. -The frequency grid is taken from the first dust type, so all dust types should -share the same frequency grid. Summed over frequency, ``specific_energy_nu`` -recovers the total ``specific_energy``. This option is disabled by default and -works for all grid types, including AMR. +By default, the ISRF is binned onto the frequency grid of the first dust type. +You can instead specify your own frequency grid (in Hz), independent of the dust +properties, by passing the ``frequencies`` argument:: + + import numpy as np + m.compute_isrf(True, frequencies=np.logspace(11., 16., 100)) + +Photons are binned to the nearest frequency in log space. Summed over frequency, +``specific_energy_nu`` recovers the total ``specific_energy``. This option is +disabled by default and works for all grid types, including AMR. ``specific_energy_nu`` is not currently returned by ``get_quantities`` and should be read directly from the output file, e.g.:: diff --git a/hyperion/model/tests/test_isrf.py b/hyperion/model/tests/test_isrf.py index ed988920..df2c14fa 100644 --- a/hyperion/model/tests/test_isrf.py +++ b/hyperion/model/tests/test_isrf.py @@ -40,7 +40,8 @@ def _assert_isrf_reconstructs_specific_energy(filename): np.testing.assert_allclose(nu_sum[heated], se[heated], rtol=1.e-6) -def _cartesian_isrf_model(compute_isrf, n_cells=2, n_photons=100000, density=1.e-18): +def _cartesian_isrf_model(compute_isrf, n_cells=2, n_photons=100000, density=1.e-18, + frequencies=None): m = Model() edges = np.linspace(-1., 1., n_cells + 1) m.set_cartesian_grid(edges, edges, edges) @@ -53,7 +54,7 @@ def _cartesian_isrf_model(compute_isrf, n_cells=2, n_photons=100000, density=1.e m.set_n_photons(initial=n_photons, imaging=0) m.set_seed(-12345) m.conf.output.output_specific_energy = 'last' - m.compute_isrf(compute_isrf) + m.compute_isrf(compute_isrf, frequencies=frequencies) return m @@ -144,6 +145,48 @@ def test_isrf_amr(tmpdir): _assert_isrf_reconstructs_specific_energy(out.filename) +@pytest.mark.requires_hyperion_binaries +def test_isrf_custom_frequency_grid(tmpdir): + # A user-specified ISRF frequency grid (independent of the dust frequency + # grid) should be used for the output bins, and energy must still be + # conserved when summing over frequency. + frequencies = np.logspace(10., 16., 12) + m = _cartesian_isrf_model(True, density=1.e-16, frequencies=frequencies) + m.write(tmpdir.join(random_id()).strpath) + out = m.run(tmpdir.join(random_id()).strpath) + bins = _read_dataset(out.filename, '/ISRF_frequency_bins') + np.testing.assert_allclose(bins, frequencies, rtol=1.e-6) + se_nu = _read_dataset(out.filename, '/specific_energy_nu') + assert se_nu.shape[0] == len(frequencies) + _assert_isrf_reconstructs_specific_energy(out.filename) + + +def test_isrf_frequencies_roundtrip(tmpdir): + # The custom ISRF frequency grid should survive a write/read round-trip. + from ...conf.conf_files import RunConf + frequencies = np.logspace(10., 16., 8) + conf = RunConf() + conf.compute_isrf(True, frequencies=frequencies) + with h5py.File(tmpdir.join('conf.h5').strpath, 'w') as f: + conf._write_isrf(f) + conf2 = RunConf() + conf2._read_isrf(f) + assert conf2.isrf + np.testing.assert_allclose(conf2.isrf_frequencies, frequencies) + + +def test_isrf_frequencies_validation(): + from ...conf.conf_files import RunConf + conf = RunConf() + with pytest.raises(ValueError): + conf.compute_isrf(True, frequencies=[[1.e10, 1.e12], [1.e14, 1.e16]]) + with pytest.raises(ValueError): + conf.compute_isrf(True, frequencies=[1.e10, -1.e12]) + # Default (no custom grid) leaves isrf_frequencies unset + conf.compute_isrf(True) + assert conf.isrf_frequencies is None + + @pytest.mark.requires_hyperion_binaries def test_isrf_mpi_matches_serial(tmpdir): # The saved ISRF must not depend on the number of MPI processes. We compare From e0c8545a9f8c9d93f4d4ccc388995200c6cf32d6 Mon Sep 17 00:00:00 2001 From: Thomas Robitaille Date: Wed, 10 Jun 2026 16:30:18 +0000 Subject: [PATCH 28/38] Expose specific_energy_nu through get_quantities and skip the ISRF frequency metadata in grid quantity reads --- hyperion/grid/cartesian_grid.py | 2 ++ hyperion/grid/cylindrical_polar_grid.py | 2 ++ hyperion/grid/grid_helpers.py | 12 +++++++----- hyperion/grid/octree_grid.py | 2 ++ hyperion/grid/spherical_polar_grid.py | 2 ++ hyperion/grid/voronoi_grid.py | 2 ++ hyperion/model/model.py | 22 ++++++++++------------ 7 files changed, 27 insertions(+), 17 deletions(-) diff --git a/hyperion/grid/cartesian_grid.py b/hyperion/grid/cartesian_grid.py index 7676162f..ac280672 100644 --- a/hyperion/grid/cartesian_grid.py +++ b/hyperion/grid/cartesian_grid.py @@ -279,6 +279,8 @@ def read_quantities(self, group, quantities='all'): # Read in physical quantities if quantities is not None: for quantity in group: + if quantity == 'ISRF_frequency_bins': + continue # per-frequency metadata, not a per-cell grid quantity if quantities == 'all' or quantity in quantities: array = np.array(group[quantity]) if array.ndim == 4: # if array is 4D, it is a list of 3D arrays diff --git a/hyperion/grid/cylindrical_polar_grid.py b/hyperion/grid/cylindrical_polar_grid.py index f78342ba..78c11bc6 100644 --- a/hyperion/grid/cylindrical_polar_grid.py +++ b/hyperion/grid/cylindrical_polar_grid.py @@ -307,6 +307,8 @@ def read_quantities(self, group, quantities='all'): # Read in physical quantities if quantities is not None: for quantity in group: + if quantity == 'ISRF_frequency_bins': + continue # per-frequency metadata, not a per-cell grid quantity if quantities == 'all' or quantity in quantities: array = np.array(group[quantity]) if array.ndim == 4: # if array is 4D, it is a list of 3D arrays diff --git a/hyperion/grid/grid_helpers.py b/hyperion/grid/grid_helpers.py index 55194e84..a60dd8d8 100644 --- a/hyperion/grid/grid_helpers.py +++ b/hyperion/grid/grid_helpers.py @@ -43,12 +43,11 @@ def single_grid_dims(data, ndim=3): elif data.ndim == ndim + 1: n_pop = data.shape[0] shape = data[0].shape - - #DN CRAZY ADDITIONS - elif data.ndim == ndim+2: + + elif data.ndim == ndim + 2: + # Frequency-resolved quantity, stored as (n_freq, n_pop, *grid) n_pop = data.shape[1] - shape = data.shape[-1] - shape = tuple(map(int,str(shape).split(' '))) + shape = data.shape[-ndim:] else: raise Exception("Unexpected number of dimensions: %i" % data.ndim) @@ -64,6 +63,9 @@ def single_grid_dims(data, ndim=3): elif len(shape) == ndim + 1: n_pop = shape[0] shape = shape[1:] + elif len(shape) == ndim + 2: + n_pop = shape[1] + shape = shape[-ndim:] else: raise Exception("Unexpected number of dimensions: %i" % len(shape)) else: diff --git a/hyperion/grid/octree_grid.py b/hyperion/grid/octree_grid.py index 8396db9d..5bae928f 100644 --- a/hyperion/grid/octree_grid.py +++ b/hyperion/grid/octree_grid.py @@ -372,6 +372,8 @@ def read_quantities(self, group, quantities='all'): # Read in physical quantities if quantities is not None: for quantity in group: + if quantity == 'ISRF_frequency_bins': + continue # per-frequency metadata, not a per-cell grid quantity if quantities == 'all' or quantity in quantities: array = np.array(group[quantity]) if array.ndim == 2: # if array is 2D, it is a list of 1D arrays diff --git a/hyperion/grid/spherical_polar_grid.py b/hyperion/grid/spherical_polar_grid.py index 5d4e7eba..451e1f17 100644 --- a/hyperion/grid/spherical_polar_grid.py +++ b/hyperion/grid/spherical_polar_grid.py @@ -317,6 +317,8 @@ def read_quantities(self, group, quantities='all'): # Read in physical quantities if quantities is not None: for quantity in group: + if quantity == 'ISRF_frequency_bins': + continue # per-frequency metadata, not a per-cell grid quantity if quantities == 'all' or quantity in quantities: array = np.array(group[quantity]) if array.ndim == 4: # if array is 4D, it is a list of 3D arrays diff --git a/hyperion/grid/voronoi_grid.py b/hyperion/grid/voronoi_grid.py index 2b03a62f..87aff852 100644 --- a/hyperion/grid/voronoi_grid.py +++ b/hyperion/grid/voronoi_grid.py @@ -396,6 +396,8 @@ def read_quantities(self, group, quantities='all'): # Read in physical quantities if quantities is not None: for quantity in group: + if quantity == 'ISRF_frequency_bins': + continue # per-frequency metadata, not a per-cell grid quantity if quantities == 'all' or quantity in quantities: array = np.array(group[quantity]) if array.ndim == 2: # if array is 2D, it is a list of 1D arrays diff --git a/hyperion/model/model.py b/hyperion/model/model.py index 567dd9aa..ec26206b 100644 --- a/hyperion/model/model.py +++ b/hyperion/model/model.py @@ -297,13 +297,12 @@ def use_quantities(self, filename, quantities=['density', 'specific_energy', 'sp else: quantities_path['specific_energy'] = last_iteration - if self.compute_isrf == True: - if 'specific_energy_nu' in quantities: - if only_initial or last_iteration is None: - if 'specific_energy_nu' in f['/Input/Grid/Quantities']: - quantities_path['specific_energy_nu'] = '/Input/Grid/Quantities' - else: - quantities_path['specific_energy_nu'] = last_iteration + if 'specific_energy_nu' in quantities: + if only_initial or last_iteration is None: + if 'specific_energy_nu' in f['/Input/Grid/Quantities']: + quantities_path['specific_energy_nu'] = '/Input/Grid/Quantities' + elif 'specific_energy_nu' in f[last_iteration]: + quantities_path['specific_energy_nu'] = last_iteration # Minimum specific energy if use_minimum_specific_energy: @@ -321,11 +320,10 @@ def use_quantities(self, filename, quantities=['density', 'specific_energy', 'sp if 'specific_energy' in f['/Grid/Quantities']: quantities_path['specific_energy'] = '/Grid/Quantities' - if self.compute_isrf == True: - if 'specific_energy_nu' in quantities: - if 'specific_energy_nu' in f['/Grid/Quantities']: - quantities_path['specific_energy_nu'] = '/Grid/Quantities' - + if 'specific_energy_nu' in quantities: + if 'specific_energy_nu' in f['/Grid/Quantities']: + quantities_path['specific_energy_nu'] = '/Grid/Quantities' + # Minimum specific energy if use_minimum_specific_energy: minimum_specific_energy_path = '/Grid/Quantities' From 9546e73d1a8d3a9e00f6b24e23f5f9ff208e8415 Mon Sep 17 00:00:00 2001 From: Thomas Robitaille Date: Wed, 10 Jun 2026 16:30:18 +0000 Subject: [PATCH 29/38] Test and document retrieving specific_energy_nu via get_quantities --- docs/setup/setup_conf.rst | 8 +++----- hyperion/model/tests/test_isrf.py | 19 +++++++++++++++++++ 2 files changed, 22 insertions(+), 5 deletions(-) diff --git a/docs/setup/setup_conf.rst b/docs/setup/setup_conf.rst index e6cc80a3..767a6c4d 100644 --- a/docs/setup/setup_conf.rst +++ b/docs/setup/setup_conf.rst @@ -223,12 +223,10 @@ Photons are binned to the nearest frequency in log space. Summed over frequency, ``specific_energy_nu`` recovers the total ``specific_energy``. This option is disabled by default and works for all grid types, including AMR. -``specific_energy_nu`` is not currently returned by ``get_quantities`` and -should be read directly from the output file, e.g.:: +``specific_energy_nu`` can be retrieved like other grid quantities, as an array +with an extra leading frequency axis:: - import h5py - with h5py.File('output.rtout', 'r') as f: - senu = f['iteration_00005/specific_energy_nu'][()] + senu = m.get_quantities()['specific_energy_nu'] Advanced parameters ------------------- diff --git a/hyperion/model/tests/test_isrf.py b/hyperion/model/tests/test_isrf.py index df2c14fa..485d9bec 100644 --- a/hyperion/model/tests/test_isrf.py +++ b/hyperion/model/tests/test_isrf.py @@ -187,6 +187,25 @@ def test_isrf_frequencies_validation(): assert conf.isrf_frequencies is None +@pytest.mark.requires_hyperion_binaries +def test_isrf_get_quantities(tmpdir): + # specific_energy_nu should be retrievable through get_quantities, the + # per-frequency ISRF_frequency_bins metadata should not pollute the grid + # quantities, and the retrieved array should still conserve energy. + m = _cartesian_isrf_model(True, density=1.e-16) + m.write(tmpdir.join(random_id()).strpath) + out = m.run(tmpdir.join(random_id()).strpath) + g = out.get_quantities() + assert 'specific_energy_nu' in g.quantities + assert 'ISRF_frequency_bins' not in g.quantities + se = np.array(g.quantities['specific_energy']) # (dust, nz, ny, nx) + se_nu = np.array(g.quantities['specific_energy_nu']) # (nu, dust, nz, ny, nx) + assert se_nu.shape[0] == 2 # number of dust frequencies + nu_sum = se_nu.sum(axis=0) + heated = se > se.min() * (1. + 1.e-6) + np.testing.assert_allclose(nu_sum[heated], se[heated], rtol=1.e-6) + + @pytest.mark.requires_hyperion_binaries def test_isrf_mpi_matches_serial(tmpdir): # The saved ISRF must not depend on the number of MPI processes. We compare From 2392b8bbc32fb7329b50a05d988768431457443c Mon Sep 17 00:00:00 2001 From: Thomas Robitaille Date: Wed, 10 Jun 2026 16:34:22 +0000 Subject: [PATCH 30/38] Add a Voronoi ISRF regression test --- hyperion/model/tests/test_isrf.py | 25 +++++++++++++++++++++++++ 1 file changed, 25 insertions(+) diff --git a/hyperion/model/tests/test_isrf.py b/hyperion/model/tests/test_isrf.py index 485d9bec..b959b9dd 100644 --- a/hyperion/model/tests/test_isrf.py +++ b/hyperion/model/tests/test_isrf.py @@ -187,6 +187,31 @@ def test_isrf_frequencies_validation(): assert conf.isrf_frequencies is None +@pytest.mark.requires_hyperion_binaries +def test_isrf_voronoi(tmpdir): + # The ISRF must also work on (unstructured) Voronoi grids, where cells are a + # flat list rather than a structured grid. + rng = np.random.RandomState(12345) + n = 30 + x = rng.uniform(-1., 1., n) + y = rng.uniform(-1., 1., n) + z = rng.uniform(-1., 1., n) + m = Model() + m.set_voronoi_grid(x, y, z) + m.add_density_grid(np.ones(n) * 1.e-16, get_test_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.compute_isrf(True) + m.write(tmpdir.join(random_id()).strpath) + out = m.run(tmpdir.join(random_id()).strpath) + _assert_isrf_reconstructs_specific_energy(out.filename) + + @pytest.mark.requires_hyperion_binaries def test_isrf_get_quantities(tmpdir): # specific_energy_nu should be retrievable through get_quantities, the From f0bc62dd16bf3394a14d43370d7e2bea508ebebd Mon Sep 17 00:00:00 2001 From: Thomas Robitaille Date: Wed, 10 Jun 2026 16:47:55 +0000 Subject: [PATCH 31/38] Skip the frequency-resolved ISRF quantity when exporting a grid to yt --- hyperion/grid/yt3_wrappers.py | 9 +++++++++ hyperion/model/tests/test_isrf.py | 17 +++++++++++++++++ 2 files changed, 26 insertions(+) diff --git a/hyperion/grid/yt3_wrappers.py b/hyperion/grid/yt3_wrappers.py index 5c9aa426..ef2eb1ce 100644 --- a/hyperion/grid/yt3_wrappers.py +++ b/hyperion/grid/yt3_wrappers.py @@ -79,6 +79,9 @@ def amr_grid_to_yt_stream(levels, dust_id=0): grid_dict['level'] = ilevel for field in grid.quantities: + if not isinstance(grid.quantities[field], list): + logger.warning("Skipping frequency-resolved quantity '{0}' in yt export (select a single frequency first)".format(field)) + continue grid_dict[('gas', field)] = grid.quantities[field][dust_id].transpose() grid_data.append(grid_dict) @@ -158,6 +161,9 @@ def octree_grid_to_yt_stream(grid, dust_id=0): quantities = {} for field in grid.quantities: + if not isinstance(grid.quantities[field], list): + logger.warning("Skipping frequency-resolved quantity '{0}' in yt export (select a single frequency first)".format(field)) + continue quantities[('gas', field)] = np.atleast_2d(grid.quantities[field][dust_id][order][~refined]).transpose() bbox = np.array([[xmin, xmax], [ymin, ymax], [zmin, zmax]]) @@ -180,6 +186,9 @@ def cartesian_grid_to_yt_stream(grid, xmin, xmax, ymin, ymax, zmin, zmax, dust_i # Make data dict which should contain (array, unit) tuples data = {} for field in grid.quantities: + if not isinstance(grid.quantities[field], list): + logger.warning("Skipping frequency-resolved quantity '{0}' in yt export (select a single frequency first)".format(field)) + continue data[field] = grid.quantities[field][dust_id].transpose(), '' # Load cartesian grid into yt diff --git a/hyperion/model/tests/test_isrf.py b/hyperion/model/tests/test_isrf.py index b959b9dd..55b206ad 100644 --- a/hyperion/model/tests/test_isrf.py +++ b/hyperion/model/tests/test_isrf.py @@ -231,6 +231,23 @@ def test_isrf_get_quantities(tmpdir): np.testing.assert_allclose(nu_sum[heated], se[heated], rtol=1.e-6) +@pytest.mark.requires_hyperion_binaries +def test_isrf_yt_export_skips_frequency_resolved(tmpdir): + # The frequency-resolved specific_energy_nu is not a per-cell scalar, so it + # must be skipped (not exported as a malformed field) when converting to yt, + # while the ordinary scalar quantities are still exported. + pytest.importorskip("yt") + m = _cartesian_isrf_model(True, density=1.e-16) + m.write(tmpdir.join(random_id()).strpath) + out = m.run(tmpdir.join(random_id()).strpath) + g = out.get_quantities() + assert 'specific_energy_nu' in g.quantities + ds = g.to_yt() + fields = [str(f) for f in ds.field_list] + assert not any('specific_energy_nu' in f for f in fields) + assert any(f.endswith("'density')") for f in fields) + + @pytest.mark.requires_hyperion_binaries def test_isrf_mpi_matches_serial(tmpdir): # The saved ISRF must not depend on the number of MPI processes. We compare From 2e4176805809990f0bbfcd06ced4c447039f4540 Mon Sep 17 00:00:00 2001 From: Thomas Robitaille Date: Wed, 10 Jun 2026 17:23:02 +0000 Subject: [PATCH 32/38] Control the frequency-resolved specific energy via output_specific_energy_nu and drop the ISRF naming --- docs/setup/setup_conf.rst | 42 +++-- hyperion/conf/conf_files.py | 73 ++++----- hyperion/grid/cartesian_grid.py | 2 +- hyperion/grid/cylindrical_polar_grid.py | 2 +- hyperion/grid/octree_grid.py | 2 +- hyperion/grid/spherical_polar_grid.py | 2 +- hyperion/grid/voronoi_grid.py | 2 +- hyperion/model/model.py | 6 - ...est_isrf.py => test_specific_energy_nu.py} | 145 +++++++++--------- src/grid/grid_generic.f90 | 16 +- src/grid/grid_physics_3d.f90 | 84 +++++----- src/grid/grid_propagate_3d.f90 | 8 +- src/main/settings.f90 | 9 +- src/main/setup_rt.f90 | 29 +++- 14 files changed, 216 insertions(+), 206 deletions(-) rename hyperion/model/tests/{test_isrf.py => test_specific_energy_nu.py} (59%) diff --git a/docs/setup/setup_conf.rst b/docs/setup/setup_conf.rst index 767a6c4d..9b3682fa 100644 --- a/docs/setup/setup_conf.rst +++ b/docs/setup/setup_conf.rst @@ -195,33 +195,41 @@ only after last iteration), or ``none`` (do not output). The default is to output only the last iteration of ``specific_energy``. To find out how to view these values, see :doc:`../postprocessing/postprocessing` -Interstellar radiation field ----------------------------- +Frequency-resolved specific energy +---------------------------------- -By default, Hyperion saves the total specific energy absorbed in each cell. If -you also want the frequency-resolved specific energy absorbed (the interstellar -radiation field, ISRF), you can enable this with:: +In addition to the total specific energy absorbed in each cell +(``output_specific_energy``), Hyperion can save the specific energy absorbed as +a function of frequency. This is controlled in the same way as the other output +quantities -- there is no separate switch to enable it:: - m.compute_isrf(True) + m.conf.output.output_specific_energy_nu = 'all' # every iteration + # 'last' -> only final iteration + # 'none' -> not computed or saved (default) -When enabled, two extra arrays are written out following the same -``output_specific_energy`` setting described above: +When it is not ``'none'``, two extra arrays are written out: * ``specific_energy_nu`` -- the specific energy absorbed in each cell as a function of frequency, in erg/s/g. -* ``ISRF_frequency_bins`` -- the frequencies (in Hz) corresponding to the last - axis of ``specific_energy_nu``. +* ``specific_energy_nu_frequencies`` -- the frequencies (in Hz) corresponding to + the leading axis of ``specific_energy_nu``. -By default, the ISRF is binned onto the frequency grid of the first dust type. -You can instead specify your own frequency grid (in Hz), independent of the dust -properties, by passing the ``frequencies`` argument:: +This is the *absorbed* (deposited) energy spectrum, not the mean intensity: each +contribution is weighted by the dust absorption opacity, so it is proportional +to :math:`\kappa_\nu J_\nu`. To recover the mean intensity :math:`J_\nu` (the +radiation field), divide by the dust absorption opacity at each frequency. +Summed over frequency, ``specific_energy_nu`` recovers the total +``specific_energy``. + +By default the binning uses the frequency grid of the first dust type. You can +instead provide your own frequency grid (in Hz), independent of the dust +properties:: import numpy as np - m.compute_isrf(True, frequencies=np.logspace(11., 16., 100)) + m.set_specific_energy_nu_frequencies(np.logspace(11., 16., 100)) -Photons are binned to the nearest frequency in log space. Summed over frequency, -``specific_energy_nu`` recovers the total ``specific_energy``. This option is -disabled by default and works for all grid types, including AMR. +Photons are binned to the nearest frequency in log space. This works for all +grid types, including AMR and Voronoi. ``specific_energy_nu`` can be retrieved like other grid quantities, as an array with an extra leading frequency axis:: diff --git a/hyperion/conf/conf_files.py b/hyperion/conf/conf_files.py index f5d26222..6046828c 100644 --- a/hyperion/conf/conf_files.py +++ b/hyperion/conf/conf_files.py @@ -18,7 +18,7 @@ def __init__(self): self.output_density = 'none' self.output_density_diff = 'none' self.output_specific_energy = 'last' - self.output_specific_energy_nu = 'last' + self.output_specific_energy_nu = 'none' self.output_n_photons = 'none' self._freeze() @@ -28,7 +28,10 @@ def read(cls, group): self.output_density = group.attrs['output_density'].decode('utf-8') self.output_density_diff = group.attrs['output_density_diff'].decode('utf-8') self.output_specific_energy = group.attrs['output_specific_energy'].decode('utf-8') - self.output_specific_energy_nu = group.attrs['output_specific_energy_nu'].decode('utf-8') + if 'output_specific_energy_nu' in group.attrs: + self.output_specific_energy_nu = group.attrs['output_specific_energy_nu'].decode('utf-8') + else: + self.output_specific_energy_nu = 'none' self.output_n_photons = group.attrs['output_n_photons'].decode('utf-8') return self @@ -55,7 +58,7 @@ def __init__(self): self.set_max_reabsorptions(1000000) self.set_pda(False) self.set_mrw(False) - self.compute_isrf(False) + self.specific_energy_nu_frequencies = None self.set_convergence(False) self.set_kill_on_absorb(False) @@ -396,52 +399,40 @@ def _read_pda(self, group): def _write_pda(self, group): group.attrs['pda'] = bool2str(self.pda) - def _read_isrf(self, group): - self.isrf = str2bool(group.attrs['isrf']) - if 'isrf_frequencies' in group: - self.isrf_frequencies = np.array(group['isrf_frequencies']['nu']) + def _read_specific_energy_nu_frequencies(self, group): + if 'specific_energy_nu_frequencies' in group: + self.specific_energy_nu_frequencies = np.array(group['specific_energy_nu_frequencies']['nu']) else: - self.isrf_frequencies = None + self.specific_energy_nu_frequencies = None - def _write_isrf(self, group): - group.attrs['isrf'] = bool2str(self.isrf) - if self.isrf and self.isrf_frequencies is not None: - group.create_dataset('isrf_frequencies', - data=np.array(list(zip(self.isrf_frequencies)), + def _write_specific_energy_nu_frequencies(self, group): + if self.specific_energy_nu_frequencies is not None: + group.create_dataset('specific_energy_nu_frequencies', + data=np.array(list(zip(self.specific_energy_nu_frequencies)), dtype=[('nu', float)])) - def compute_isrf(self, isrf, frequencies=None): + def set_specific_energy_nu_frequencies(self, frequencies): ''' - Set whether to compute and save the interstellar radiation field (ISRF) - in each cell. + Set the frequency grid onto which the frequency-resolved specific energy + (``specific_energy_nu``) is binned. - If enabled, the specific energy absorbed is also accumulated as a - function of frequency and saved as the ``specific_energy_nu`` quantity - (in erg/s/g), alongside the ``ISRF_frequency_bins`` array giving the - corresponding frequencies (in Hz). Photons are binned to the nearest - frequency in log space. These arrays are written following the - ``output_specific_energy`` setting. This is disabled by default. + This is only relevant if ``conf.output.output_specific_energy_nu`` is set + to ``'all'`` or ``'last'``. If this method is not called, the frequency + grid of the first dust type is used. Photons are binned to the nearest + frequency in log space. Parameters ---------- - isrf : bool - Whether to compute and save the ISRF in each cell. - frequencies : iterable of float, optional - The frequencies (in Hz) onto which to bin the ISRF. If not - specified, the frequency grid of the first dust type is used. Only - used when ``isrf`` is `True`. - ''' - self.isrf = isrf - if frequencies is None: - self.isrf_frequencies = None - else: - frequencies = np.asarray(frequencies, dtype=float) - if frequencies.ndim != 1 or frequencies.size < 1: - raise ValueError("frequencies should be a 1-d array of at least one value") - if np.any(frequencies <= 0.): - raise ValueError("frequencies should be positive (in Hz)") - self.isrf_frequencies = np.sort(frequencies) + frequencies : iterable of float + The frequencies (in Hz) onto which to bin ``specific_energy_nu``. + ''' + frequencies = np.asarray(frequencies, dtype=float) + if frequencies.ndim != 1 or frequencies.size < 1: + raise ValueError("frequencies should be a 1-d array of at least one value") + if np.any(frequencies <= 0.): + raise ValueError("frequencies should be positive (in Hz)") + self.specific_energy_nu_frequencies = np.sort(frequencies) def set_mrw(self, mrw, gamma=1.0, inter_max=1000, warn=True): @@ -772,7 +763,7 @@ def read_run_conf(self, group): # not a class method because inherited self._read_max_reabsorptions(group) self._read_pda(group) self._read_mrw(group) - self._read_isrf(group) + self._read_specific_energy_nu_frequencies(group) self._read_convergence(group) self._read_kill_on_absorb(group) self._read_kill_on_scatter(group) @@ -801,7 +792,7 @@ def write_run_conf(self, group): self._write_max_reabsorptions(group) self._write_pda(group) self._write_mrw(group) - self._write_isrf(group) + self._write_specific_energy_nu_frequencies(group) self._write_convergence(group) self._write_kill_on_absorb(group) self._write_kill_on_scatter(group) diff --git a/hyperion/grid/cartesian_grid.py b/hyperion/grid/cartesian_grid.py index ac280672..dc9f71cb 100644 --- a/hyperion/grid/cartesian_grid.py +++ b/hyperion/grid/cartesian_grid.py @@ -279,7 +279,7 @@ def read_quantities(self, group, quantities='all'): # Read in physical quantities if quantities is not None: for quantity in group: - if quantity == 'ISRF_frequency_bins': + if quantity == 'specific_energy_nu_frequencies': continue # per-frequency metadata, not a per-cell grid quantity if quantities == 'all' or quantity in quantities: array = np.array(group[quantity]) diff --git a/hyperion/grid/cylindrical_polar_grid.py b/hyperion/grid/cylindrical_polar_grid.py index 78c11bc6..a6a8b7cf 100644 --- a/hyperion/grid/cylindrical_polar_grid.py +++ b/hyperion/grid/cylindrical_polar_grid.py @@ -307,7 +307,7 @@ def read_quantities(self, group, quantities='all'): # Read in physical quantities if quantities is not None: for quantity in group: - if quantity == 'ISRF_frequency_bins': + if quantity == 'specific_energy_nu_frequencies': continue # per-frequency metadata, not a per-cell grid quantity if quantities == 'all' or quantity in quantities: array = np.array(group[quantity]) diff --git a/hyperion/grid/octree_grid.py b/hyperion/grid/octree_grid.py index 5bae928f..31696047 100644 --- a/hyperion/grid/octree_grid.py +++ b/hyperion/grid/octree_grid.py @@ -372,7 +372,7 @@ def read_quantities(self, group, quantities='all'): # Read in physical quantities if quantities is not None: for quantity in group: - if quantity == 'ISRF_frequency_bins': + if quantity == 'specific_energy_nu_frequencies': continue # per-frequency metadata, not a per-cell grid quantity if quantities == 'all' or quantity in quantities: array = np.array(group[quantity]) diff --git a/hyperion/grid/spherical_polar_grid.py b/hyperion/grid/spherical_polar_grid.py index 451e1f17..9fa45bdc 100644 --- a/hyperion/grid/spherical_polar_grid.py +++ b/hyperion/grid/spherical_polar_grid.py @@ -317,7 +317,7 @@ def read_quantities(self, group, quantities='all'): # Read in physical quantities if quantities is not None: for quantity in group: - if quantity == 'ISRF_frequency_bins': + if quantity == 'specific_energy_nu_frequencies': continue # per-frequency metadata, not a per-cell grid quantity if quantities == 'all' or quantity in quantities: array = np.array(group[quantity]) diff --git a/hyperion/grid/voronoi_grid.py b/hyperion/grid/voronoi_grid.py index 87aff852..81b48d2d 100644 --- a/hyperion/grid/voronoi_grid.py +++ b/hyperion/grid/voronoi_grid.py @@ -396,7 +396,7 @@ def read_quantities(self, group, quantities='all'): # Read in physical quantities if quantities is not None: for quantity in group: - if quantity == 'ISRF_frequency_bins': + if quantity == 'specific_energy_nu_frequencies': continue # per-frequency metadata, not a per-cell grid quantity if quantities == 'all' or quantity in quantities: array = np.array(group[quantity]) diff --git a/hyperion/model/model.py b/hyperion/model/model.py index ec26206b..a9143c03 100644 --- a/hyperion/model/model.py +++ b/hyperion/model/model.py @@ -819,9 +819,6 @@ class (e.g. ``SphericalDust``, ``IsotropicDust``, ...). self.grid['density'] = [] if specific_energy is not None: self.grid['specific_energy'] = [] - - if self.compute_isrf == True: - self.grid['specific_energy_nu'] = [] # Check whether the density can be added to an existing one if merge_if_possible: @@ -863,9 +860,6 @@ class (e.g. ``SphericalDust``, ``IsotropicDust``, ...). # Set specific energy if specified if specific_energy is not None: self.grid['specific_energy'].append(specific_energy) - - if self.compute_isrf == True: - self.grid['specific_energy_nu'].append(specific_energy) def set_cartesian_grid(self, x_wall, y_wall, z_wall): self.set_grid(CartesianGrid(x_wall, y_wall, z_wall)) diff --git a/hyperion/model/tests/test_isrf.py b/hyperion/model/tests/test_specific_energy_nu.py similarity index 59% rename from hyperion/model/tests/test_isrf.py rename to hyperion/model/tests/test_specific_energy_nu.py index 55b206ad..d7483da0 100644 --- a/hyperion/model/tests/test_isrf.py +++ b/hyperion/model/tests/test_specific_energy_nu.py @@ -25,11 +25,11 @@ def _read_dataset(filename, suffix): return np.asarray(f[matches[-1]][()], dtype=float) -def _assert_isrf_reconstructs_specific_energy(filename): - # Summed over frequency, the ISRF must reproduce the specific energy. Cells - # with no absorption are clamped up to the minimum specific energy (which - # only affects specific_energy, not the unclamped specific_energy_nu), so we - # compare only cells that were genuinely heated above that floor. +def _assert_nu_sums_to_specific_energy(filename): + # Summed over frequency, specific_energy_nu 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_nu), + # so we compare only cells that were genuinely heated above that floor. se = _read_dataset(filename, '/specific_energy') se_nu = _read_dataset(filename, '/specific_energy_nu') assert se_nu is not None @@ -40,8 +40,8 @@ def _assert_isrf_reconstructs_specific_energy(filename): np.testing.assert_allclose(nu_sum[heated], se[heated], rtol=1.e-6) -def _cartesian_isrf_model(compute_isrf, n_cells=2, n_photons=100000, density=1.e-18, - frequencies=None): +def _cartesian_model(output_specific_energy_nu='last', n_cells=2, n_photons=100000, + density=1.e-18, frequencies=None): m = Model() edges = np.linspace(-1., 1., n_cells + 1) m.set_cartesian_grid(edges, edges, edges) @@ -54,53 +54,54 @@ def _cartesian_isrf_model(compute_isrf, n_cells=2, n_photons=100000, density=1.e m.set_n_photons(initial=n_photons, imaging=0) m.set_seed(-12345) m.conf.output.output_specific_energy = 'last' - m.compute_isrf(compute_isrf, frequencies=frequencies) + m.conf.output.output_specific_energy_nu = output_specific_energy_nu + if frequencies is not None: + m.set_specific_energy_nu_frequencies(frequencies) return m @pytest.mark.requires_hyperion_binaries -def test_isrf_does_not_change_specific_energy(tmpdir): - # Computing the ISRF must not perturb the ordinary specific-energy - # (temperature) calculation. This guards against the regression where the - # specific_energy_sum accumulation was accidentally gated behind compute_isrf. +def test_specific_energy_nu_does_not_change_specific_energy(tmpdir): + # Computing the frequency-resolved specific energy must not perturb the + # ordinary specific-energy (temperature) calculation. This guards against the + # regression where specific_energy_sum accumulation was gated on the option. se = {} - for isrf in (False, True): - m = _cartesian_isrf_model(isrf) + for output in ('none', 'last'): + m = _cartesian_model(output_specific_energy_nu=output) m.write(tmpdir.join(random_id()).strpath) out = m.run(tmpdir.join(random_id()).strpath) - se[isrf] = _read_dataset(out.filename, '/specific_energy') - # Same seed, and the ISRF code path does not touch the random stream, so - # the specific energy should be identical with and without the ISRF. - np.testing.assert_allclose(se[True], se[False], rtol=1.e-10) + se[output] = _read_dataset(out.filename, '/specific_energy') + # Same seed, and the extra code path does not touch the random stream, so + # the specific energy should be identical with and without the option. + np.testing.assert_allclose(se['last'], se['none'], rtol=1.e-10) @pytest.mark.requires_hyperion_binaries -def test_isrf_sums_to_specific_energy(tmpdir): - # The frequency-resolved ISRF, summed over frequency, must reproduce the - # total specific energy. - m = _cartesian_isrf_model(True) +def test_specific_energy_nu_sums_to_specific_energy(tmpdir): + # specific_energy_nu, summed over frequency, must reproduce specific_energy. + m = _cartesian_model('last') m.write(tmpdir.join(random_id()).strpath) out = m.run(tmpdir.join(random_id()).strpath) - _assert_isrf_reconstructs_specific_energy(out.filename) + _assert_nu_sums_to_specific_energy(out.filename) @pytest.mark.requires_hyperion_binaries -def test_isrf_frequency_bins_written(tmpdir): - # ISRF_frequency_bins is a per-frequency (1-D) quantity, not per-cell. This - # guards against the regression where it was written as a grid array, which - # segfaulted whenever the number of frequencies was smaller than n_cells. - m = _cartesian_isrf_model(True, n_cells=8) # 512 cells, dust has 2 frequencies +def test_specific_energy_nu_frequencies_written(tmpdir): + # specific_energy_nu_frequencies is a per-frequency (1-D) quantity, not + # per-cell. This guards against the regression where it was written as a grid + # array, which segfaulted whenever n_nu was smaller than n_cells. + m = _cartesian_model('last', n_cells=8) # 512 cells, dust has 2 frequencies m.write(tmpdir.join(random_id()).strpath) out = m.run(tmpdir.join(random_id()).strpath) - bins = _read_dataset(out.filename, '/ISRF_frequency_bins') + bins = _read_dataset(out.filename, '/specific_energy_nu_frequencies') assert bins is not None assert bins.shape == (2,) # == number of dust frequencies, not n_cells @pytest.mark.requires_hyperion_binaries -def test_isrf_dustless_model_runs(tmpdir): +def test_specific_energy_nu_dustless_model_runs(tmpdir): # A model with no dust must still run (and image) without crashing. This - # guards against the regression where ISRF instrumentation dereferenced the + # guards against the regression where the instrumentation dereferenced the # first dust type unconditionally, segfaulting dustless models. m = Model() m.set_cartesian_grid([-1., 1.], [-1., 1.], [-1., 1.]) @@ -118,9 +119,9 @@ def test_isrf_dustless_model_runs(tmpdir): @pytest.mark.requires_hyperion_binaries -def test_isrf_amr(tmpdir): - # The ISRF must work on AMR grids: it should run and, summed over frequency, - # reproduce the specific energy in every cell. +def test_specific_energy_nu_amr(tmpdir): + # specific_energy_nu must work on AMR grids: it should run and, summed over + # frequency, reproduce specific_energy in every cell. m = Model() amr = AMRGrid() level = amr.add_level() @@ -139,58 +140,56 @@ def test_isrf_amr(tmpdir): m.set_n_photons(initial=100000, imaging=0) m.set_seed(-9) m.conf.output.output_specific_energy = 'last' - m.compute_isrf(True) + m.conf.output.output_specific_energy_nu = 'last' m.write(tmpdir.join(random_id()).strpath) out = m.run(tmpdir.join(random_id()).strpath) - _assert_isrf_reconstructs_specific_energy(out.filename) + _assert_nu_sums_to_specific_energy(out.filename) @pytest.mark.requires_hyperion_binaries -def test_isrf_custom_frequency_grid(tmpdir): - # A user-specified ISRF frequency grid (independent of the dust frequency - # grid) should be used for the output bins, and energy must still be - # conserved when summing over frequency. +def test_specific_energy_nu_custom_frequency_grid(tmpdir): + # A user-specified frequency grid (independent of the dust frequency grid) + # should be used for the output bins, and energy must still be conserved + # when summing over frequency. frequencies = np.logspace(10., 16., 12) - m = _cartesian_isrf_model(True, density=1.e-16, frequencies=frequencies) + m = _cartesian_model('last', density=1.e-16, frequencies=frequencies) m.write(tmpdir.join(random_id()).strpath) out = m.run(tmpdir.join(random_id()).strpath) - bins = _read_dataset(out.filename, '/ISRF_frequency_bins') + bins = _read_dataset(out.filename, '/specific_energy_nu_frequencies') np.testing.assert_allclose(bins, frequencies, rtol=1.e-6) se_nu = _read_dataset(out.filename, '/specific_energy_nu') assert se_nu.shape[0] == len(frequencies) - _assert_isrf_reconstructs_specific_energy(out.filename) + _assert_nu_sums_to_specific_energy(out.filename) -def test_isrf_frequencies_roundtrip(tmpdir): - # The custom ISRF frequency grid should survive a write/read round-trip. +def test_specific_energy_nu_frequencies_roundtrip(tmpdir): + # The custom frequency grid should survive a write/read round-trip. from ...conf.conf_files import RunConf frequencies = np.logspace(10., 16., 8) conf = RunConf() - conf.compute_isrf(True, frequencies=frequencies) + conf.set_specific_energy_nu_frequencies(frequencies) with h5py.File(tmpdir.join('conf.h5').strpath, 'w') as f: - conf._write_isrf(f) + conf._write_specific_energy_nu_frequencies(f) conf2 = RunConf() - conf2._read_isrf(f) - assert conf2.isrf - np.testing.assert_allclose(conf2.isrf_frequencies, frequencies) + conf2._read_specific_energy_nu_frequencies(f) + np.testing.assert_allclose(conf2.specific_energy_nu_frequencies, frequencies) -def test_isrf_frequencies_validation(): +def test_specific_energy_nu_frequencies_validation(): from ...conf.conf_files import RunConf conf = RunConf() with pytest.raises(ValueError): - conf.compute_isrf(True, frequencies=[[1.e10, 1.e12], [1.e14, 1.e16]]) + conf.set_specific_energy_nu_frequencies([[1.e10, 1.e12], [1.e14, 1.e16]]) with pytest.raises(ValueError): - conf.compute_isrf(True, frequencies=[1.e10, -1.e12]) - # Default (no custom grid) leaves isrf_frequencies unset - conf.compute_isrf(True) - assert conf.isrf_frequencies is None + conf.set_specific_energy_nu_frequencies([1.e10, -1.e12]) + # Default: no custom grid set + assert conf.specific_energy_nu_frequencies is None @pytest.mark.requires_hyperion_binaries -def test_isrf_voronoi(tmpdir): - # The ISRF must also work on (unstructured) Voronoi grids, where cells are a - # flat list rather than a structured grid. +def test_specific_energy_nu_voronoi(tmpdir): + # specific_energy_nu must also work on (unstructured) Voronoi grids, where + # cells are a flat list rather than a structured grid. rng = np.random.RandomState(12345) n = 30 x = rng.uniform(-1., 1., n) @@ -206,23 +205,23 @@ def test_isrf_voronoi(tmpdir): m.set_n_photons(initial=100000, imaging=0) m.set_seed(-12345) m.conf.output.output_specific_energy = 'last' - m.compute_isrf(True) + m.conf.output.output_specific_energy_nu = 'last' m.write(tmpdir.join(random_id()).strpath) out = m.run(tmpdir.join(random_id()).strpath) - _assert_isrf_reconstructs_specific_energy(out.filename) + _assert_nu_sums_to_specific_energy(out.filename) @pytest.mark.requires_hyperion_binaries -def test_isrf_get_quantities(tmpdir): +def test_specific_energy_nu_get_quantities(tmpdir): # specific_energy_nu should be retrievable through get_quantities, the - # per-frequency ISRF_frequency_bins metadata should not pollute the grid - # quantities, and the retrieved array should still conserve energy. - m = _cartesian_isrf_model(True, density=1.e-16) + # per-frequency frequencies metadata should not pollute the grid quantities, + # and the retrieved array should still conserve energy. + m = _cartesian_model('last', density=1.e-16) m.write(tmpdir.join(random_id()).strpath) out = m.run(tmpdir.join(random_id()).strpath) g = out.get_quantities() assert 'specific_energy_nu' in g.quantities - assert 'ISRF_frequency_bins' not in g.quantities + assert 'specific_energy_nu_frequencies' not in g.quantities se = np.array(g.quantities['specific_energy']) # (dust, nz, ny, nx) se_nu = np.array(g.quantities['specific_energy_nu']) # (nu, dust, nz, ny, nx) assert se_nu.shape[0] == 2 # number of dust frequencies @@ -232,12 +231,12 @@ def test_isrf_get_quantities(tmpdir): @pytest.mark.requires_hyperion_binaries -def test_isrf_yt_export_skips_frequency_resolved(tmpdir): +def test_specific_energy_nu_yt_export_skips_frequency_resolved(tmpdir): # The frequency-resolved specific_energy_nu is not a per-cell scalar, so it # must be skipped (not exported as a malformed field) when converting to yt, # while the ordinary scalar quantities are still exported. pytest.importorskip("yt") - m = _cartesian_isrf_model(True, density=1.e-16) + m = _cartesian_model('last', density=1.e-16) m.write(tmpdir.join(random_id()).strpath) out = m.run(tmpdir.join(random_id()).strpath) g = out.get_quantities() @@ -249,14 +248,14 @@ def test_isrf_yt_export_skips_frequency_resolved(tmpdir): @pytest.mark.requires_hyperion_binaries -def test_isrf_mpi_matches_serial(tmpdir): - # The saved ISRF must not depend on the number of MPI processes. We compare - # the total ISRF energy (which is conserved, hence robust to Monte Carlo - # noise) between a serial run and a 2-process MPI run. +def test_specific_energy_nu_mpi_matches_serial(tmpdir): + # The saved spectrum must not depend on the number of MPI processes. We + # compare the total (which is conserved, hence robust to Monte Carlo noise) + # between a serial run and a 2-process MPI run. if shutil.which('mpirun') is None and shutil.which('mpiexec') is None: pytest.skip("no MPI launcher available") - m = _cartesian_isrf_model(True, n_photons=200000) + m = _cartesian_model('last', n_photons=200000) m.write(tmpdir.join(random_id()).strpath) out_serial = m.run(tmpdir.join(random_id()).strpath) diff --git a/src/grid/grid_generic.f90 b/src/grid/grid_generic.f90 index 326560bc..786296b5 100644 --- a/src/grid/grid_generic.f90 +++ b/src/grid/grid_generic.f90 @@ -6,8 +6,8 @@ module grid_generic use grid_io, only : write_grid_3d, write_grid_4d, write_grid_5d use grid_geometry, only : geo - use grid_physics, only : n_photons, last_photon_id, specific_energy_sum, specific_energy_sum_nu, specific_energy, specific_energy_nu, isrf_nu, density, density_original - use settings, only : output_n_photons, output_specific_energy, output_density, output_density_diff, physics_io_type, compute_isrf + use grid_physics, only : n_photons, last_photon_id, specific_energy_sum, specific_energy_sum_nu, specific_energy, specific_energy_nu, nu_bins, density, density_original + use settings, only : output_n_photons, output_specific_energy, output_specific_energy_nu, output_density, output_density_diff, physics_io_type, compute_specific_energy_nu implicit none save @@ -64,14 +64,14 @@ subroutine output_grid(group, iter, n_iter) - ! WRITE THE ISRF - if (compute_isrf) then + ! WRITE THE FREQUENCY-RESOLVED SPECIFIC ENERGY + if (compute_specific_energy_nu) then - if(trim(output_specific_energy)=='all' .or. (trim(output_specific_energy)=='last'.and.iter==n_iter)) then + if(trim(output_specific_energy_nu)=='all' .or. (trim(output_specific_energy_nu)=='last'.and.iter==n_iter)) then - ! ISRF frequency bins are a per-frequency quantity, not per-cell, so - ! they are written as a plain 1-D dataset rather than a grid array. - call mp_write_array(group, 'ISRF_frequency_bins', isrf_nu) + ! The frequencies are a per-frequency quantity, not per-cell, so they + ! are written as a plain 1-D dataset rather than a grid array. + call mp_write_array(group, 'specific_energy_nu_frequencies', nu_bins) if(allocated(specific_energy_nu)) then select case(physics_io_type) diff --git a/src/grid/grid_physics_3d.f90 b/src/grid/grid_physics_3d.f90 index de743a5e..e48a5958 100644 --- a/src/grid/grid_physics_3d.f90 +++ b/src/grid/grid_physics_3d.f90 @@ -39,8 +39,8 @@ module grid_physics real(dp),allocatable, public :: specific_energy_nu(:,:,:) real(dp),allocatable, public :: specific_energy_sum(:,:) real(dp),allocatable, public :: specific_energy_sum_nu(:,:,:) - real(dp),allocatable, public :: isrf_nu(:) - real(dp),allocatable, public :: log_isrf_nu(:) + real(dp),allocatable, public :: nu_bins(:) + real(dp),allocatable, public :: log_nu_bins(:) real(dp),allocatable, public :: specific_energy_additional(:,:) real(dp),allocatable, public :: specific_energy_additional_nu(:,:,:) @@ -99,25 +99,25 @@ integer function select_dust_specific_energy_rho(icell) result(id_select) id_select = sample_pdf(absorption) end function select_dust_specific_energy_rho - subroutine setup_grid_physics(group, use_mrw, use_pda, compute_isrf) + subroutine setup_grid_physics(group, use_mrw, use_pda, compute_specific_energy_nu) implicit none integer(hid_t),intent(in) :: group - logical,intent(in) :: use_mrw, use_pda, compute_isrf - integer :: n_isrf_wavelengths + logical,intent(in) :: use_mrw, use_pda, compute_specific_energy_nu + integer :: n_nu_bins - ! The ISRF is binned onto a user-specified frequency grid if one was given, + ! specific_energy_nu is binned onto a user-specified frequency grid if one was given, ! otherwise onto the frequency grid of the first dust type. - if (allocated(isrf_frequencies)) then - n_isrf_wavelengths = size(isrf_frequencies) + if (allocated(specific_energy_nu_frequencies)) then + n_nu_bins = size(specific_energy_nu_frequencies) else - n_isrf_wavelengths = d(1)%n_nu + n_nu_bins = d(1)%n_nu end if ! Density allocate(density(geo%n_cells, n_dust)) allocate(specific_energy(geo%n_cells, n_dust)) - allocate(specific_energy_nu(geo%n_cells,n_dust,n_isrf_wavelengths)) + allocate(specific_energy_nu(geo%n_cells,n_dust,n_nu_bins)) if(n_dust > 0) then @@ -175,8 +175,8 @@ subroutine setup_grid_physics(group, use_mrw, use_pda, compute_isrf) specific_energy(:, id) = 0. end where - if (compute_isrf) then - do idx=1,n_isrf_wavelengths + if (compute_specific_energy_nu) then + do idx=1,n_nu_bins where(.not.geo%mask) specific_energy_nu(:, id, idx) = 0. endwhere @@ -188,8 +188,8 @@ subroutine setup_grid_physics(group, use_mrw, use_pda, compute_isrf) if(trim(specific_energy_type) == 'additional') then allocate(specific_energy_additional(geo%n_cells, n_dust)) - if (compute_isrf) then - allocate(specific_energy_additional_nu(geo%n_cells,n_dust,n_isrf_wavelengths)) + if (compute_specific_energy_nu) then + allocate(specific_energy_additional_nu(geo%n_cells,n_dust,n_nu_bins)) end if ! We store a copy of the initial specific energy in a separate ! array, and we set the specific energy to the minimum specific @@ -200,8 +200,8 @@ subroutine setup_grid_physics(group, use_mrw, use_pda, compute_isrf) do id=1,n_dust specific_energy(:,id) = minimum_specific_energy(id) - if (compute_isrf) then - do idx=1,n_isrf_wavelengths + if (compute_specific_energy_nu) then + do idx=1,n_nu_bins specific_energy_nu(:,id,idx) = minimum_specific_energy(id) end do end if @@ -220,8 +220,8 @@ subroutine setup_grid_physics(group, use_mrw, use_pda, compute_isrf) do id=1,n_dust specific_energy(:,id) = minimum_specific_energy(id) - if (compute_isrf) then - do idx = 1,n_isrf_wavelengths + if (compute_specific_energy_nu) then + do idx = 1,n_nu_bins specific_energy_nu(:,id,idx) = minimum_specific_energy(id) end do end if @@ -239,24 +239,24 @@ subroutine setup_grid_physics(group, use_mrw, use_pda, compute_isrf) allocate(specific_energy_sum(geo%n_cells, n_dust)) specific_energy_sum = 0._dp - ! Set up basics for ISRF calculation - allocate(specific_energy_sum_nu(geo%n_cells, n_dust, n_isrf_wavelengths)) + ! Set up basics for the frequency-resolved specific energy + allocate(specific_energy_sum_nu(geo%n_cells, n_dust, n_nu_bins)) specific_energy_sum_nu = 0._dp - ! Cache the ISRF frequency grid (and its log) once so they do not have to be + ! 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 ! log-frequency space (see grid_propagate). - if (compute_isrf) then - allocate(isrf_nu(n_isrf_wavelengths)) - if (allocated(isrf_frequencies)) then - isrf_nu = isrf_frequencies + if (compute_specific_energy_nu) then + allocate(nu_bins(n_nu_bins)) + if (allocated(specific_energy_nu_frequencies)) then + nu_bins = specific_energy_nu_frequencies else - do idx=1,n_isrf_wavelengths - isrf_nu(idx) = d(1)%nu(idx) + do idx=1,n_nu_bins + nu_bins(idx) = d(1)%nu(idx) end do end if - allocate(log_isrf_nu(n_isrf_wavelengths)) - log_isrf_nu = log10(isrf_nu) + allocate(log_nu_bins(n_nu_bins)) + log_nu_bins = log10(nu_bins) end if ! Total energy absorbed @@ -325,9 +325,9 @@ subroutine sublimate_dust() implicit none integer :: ic, id integer :: reset - integer :: n_isrf_wavelengths + integer :: n_nu_bins - n_isrf_wavelengths = size(specific_energy_nu, 3) + n_nu_bins = size(specific_energy_nu, 3) reset = 0 @@ -342,8 +342,8 @@ subroutine sublimate_dust() specific_energy(ic, id) = minimum_specific_energy(id) reset = reset + 1 - if (compute_isrf) then - do idx=1,n_isrf_wavelengths + if (compute_specific_energy_nu) then + do idx=1,n_nu_bins specific_energy_nu(ic,id,idx) = minimum_specific_energy(id) end do end if @@ -365,8 +365,8 @@ subroutine sublimate_dust() reset = reset + 1 - if (compute_isrf) then - do idx=1,n_isrf_wavelengths + if (compute_specific_energy_nu) then + do idx=1,n_nu_bins specific_energy_nu(ic,id,idx) = minimum_specific_energy(id) end do end if @@ -385,8 +385,8 @@ subroutine sublimate_dust() end if - if (compute_isrf) then - do idx=1,n_isrf_wavelengths + if (compute_specific_energy_nu) then + do idx=1,n_nu_bins specific_energy_nu(ic,id,idx) = minimum_specific_energy(id) end do end if @@ -412,17 +412,17 @@ subroutine update_energy_abs(scale) real(dp), intent(in) :: scale integer :: id,idx - integer :: n_isrf_wavelengths + integer :: n_nu_bins - n_isrf_wavelengths = size(specific_energy_nu, 3) + n_nu_bins = size(specific_energy_nu, 3) if(main_process()) write(*,'(" [grid_physics] updating energy_abs")') do id=1,n_dust specific_energy(:,id) = specific_energy_sum(:,id) * scale / geo%volume - if (compute_isrf) then - do idx=1,n_isrf_wavelengths + if (compute_specific_energy_nu) then + do idx=1,n_nu_bins specific_energy_nu(:,id,idx) = specific_energy_sum_nu(:,id,idx) * scale/geo%volume end do end if @@ -444,7 +444,7 @@ subroutine update_energy_abs(scale) specific_energy = specific_energy + specific_energy_additional - if (compute_isrf) then + if (compute_specific_energy_nu) then specific_energy_nu = specific_energy_nu + specific_energy_additional_nu end if diff --git a/src/grid/grid_propagate_3d.f90 b/src/grid/grid_propagate_3d.f90 index 0f57ceb2..3ccac335 100644 --- a/src/grid/grid_propagate_3d.f90 +++ b/src/grid/grid_propagate_3d.f90 @@ -5,10 +5,10 @@ module grid_propagate use type_grid_cell use dust_main, only : n_dust use grid_geometry, only : escaped, find_wall, in_correct_cell, next_cell, opposite_wall - use grid_physics, only : specific_energy_sum, specific_energy_sum_nu, log_isrf_nu, density, n_photons, last_photon_id + use grid_physics, only : specific_energy_sum, specific_energy_sum_nu, log_nu_bins, density, n_photons, last_photon_id use sources use counters - use settings, only : frac_check => propagation_check_frequency, compute_isrf + use settings, only : frac_check => propagation_check_frequency, compute_specific_energy_nu implicit none save @@ -133,13 +133,13 @@ subroutine grid_integrate(p,tau_required,tau_achieved) tau_achieved = tau_achieved + tau_cell - if (compute_isrf) idx = minloc(abs(log_isrf_nu - log10(p%nu)), DIM=1) + if (compute_specific_energy_nu) 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) = & & specific_energy_sum(p%icell%ic, id) + tmin * p%current_kappa(id) * p%energy - if (compute_isrf) then + if (compute_specific_energy_nu) then specific_energy_sum_nu(p%icell%ic, id, idx) = & & specific_energy_sum_nu(p%icell%ic, id, idx) + tmin * p%current_kappa(id) * p%energy end if diff --git a/src/main/settings.f90 b/src/main/settings.f90 index 9e3102d0..faf0ccd7 100644 --- a/src/main/settings.f90 +++ b/src/main/settings.f90 @@ -20,7 +20,7 @@ module settings integer(idp),public :: n_last_photons_dust = 0 integer(idp),public :: n_raytracing_photons_sources = 0 integer(idp),public :: n_raytracing_photons_dust = 0 - logical,public :: use_raytracing, use_mrw, use_pda, compute_isrf + logical,public :: use_raytracing, use_mrw, use_pda, compute_specific_energy_nu logical, public :: kill_on_absorb, kill_on_scatter real(dp),public :: mrw_gamma logical, public :: forced_first_interaction @@ -30,15 +30,16 @@ module settings real(dp) :: monochromatic_energy_threshold real(dp),allocatable :: frequencies(:) - ! Optional user-specified frequency grid for the ISRF. If not allocated, the - ! frequency grid of the first dust type is used instead. - real(dp),allocatable :: isrf_frequencies(:) + ! Optional user-specified frequency grid for specific_energy_nu. If not + ! allocated, the frequency grid of the first dust type is used instead. + real(dp),allocatable :: specific_energy_nu_frequencies(:) integer :: physics_io_type character(len=4) :: output_density character(len=4) :: output_density_diff character(len=4) :: output_specific_energy + character(len=4) :: output_specific_energy_nu character(len=4) :: output_n_photons logical :: check_convergence = .false. diff --git a/src/main/setup_rt.f90 b/src/main/setup_rt.f90 index 462968b0..73ef5641 100644 --- a/src/main/setup_rt.f90 +++ b/src/main/setup_rt.f90 @@ -74,12 +74,29 @@ subroutine setup_initial(input_handle) call mp_read_keyword(input_handle, '/', 'pda', use_pda) call mp_read_keyword(input_handle, '/', 'mrw', use_mrw) - call mp_read_keyword(input_handle, '/', 'isrf', compute_isrf) + ! The frequency-resolved specific energy is computed whenever it is to be + ! output, i.e. unless output_specific_energy_nu is 'none'. This is read here + ! (ahead of the main OUTPUT block below) because it is needed to set up the + ! grid physics arrays. + g_output = mp_open_group(input_handle, '/Output') + if(mp_exists_keyword(g_output, '.', 'output_specific_energy_nu')) then + call mp_read_keyword(g_output, '.', 'output_specific_energy_nu', output_specific_energy_nu) + else + output_specific_energy_nu = 'none' + end if + call mp_close_group(g_output) + + if(trim(output_specific_energy_nu).ne.'all' & + & .and.trim(output_specific_energy_nu).ne.'last' & + & .and.trim(output_specific_energy_nu).ne.'none') & + & call error("setup_initial","output_specific_energy_nu should be one of all/last/none") + + compute_specific_energy_nu = trim(output_specific_energy_nu).ne.'none' - ! Optional user-specified ISRF frequency grid (otherwise the dust grid is used) - if(compute_isrf) then - if(mp_path_exists(input_handle, 'isrf_frequencies')) then - call mp_table_read_column_auto(input_handle, 'isrf_frequencies', 'nu', isrf_frequencies) + ! Optional user-specified frequency grid (otherwise the dust grid is used) + if(compute_specific_energy_nu) then + if(mp_path_exists(input_handle, 'specific_energy_nu_frequencies')) then + call mp_table_read_column_auto(input_handle, 'specific_energy_nu_frequencies', 'nu', specific_energy_nu_frequencies) end if end if @@ -181,7 +198,7 @@ subroutine setup_initial(input_handle) call mp_close_group(g_geometry) g_physics = mp_open_group(input_handle, '/Grid/Quantities') - call setup_grid_physics(g_physics, use_mrw, use_pda, compute_isrf) + call setup_grid_physics(g_physics, use_mrw, use_pda, compute_specific_energy_nu) call mp_close_group(g_physics) call mp_read_keyword(input_handle, '/', 'physics_io_bytes', physics_io_bytes) From 85141ac3c8040b91b3cf716df58fbba423a1dce8 Mon Sep 17 00:00:00 2001 From: Thomas Robitaille Date: Wed, 10 Jun 2026 17:31:31 +0000 Subject: [PATCH 33/38] Update config IO tests for output_specific_energy_nu and the frequency-grid setter --- hyperion/conf/tests/test_conf_io.py | 12 +++++++----- 1 file changed, 7 insertions(+), 5 deletions(-) diff --git a/hyperion/conf/tests/test_conf_io.py b/hyperion/conf/tests/test_conf_io.py index 738ef480..fb667e72 100644 --- a/hyperion/conf/tests/test_conf_io.py +++ b/hyperion/conf/tests/test_conf_io.py @@ -4,6 +4,7 @@ from itertools import product +import numpy as np from numpy.testing import assert_equal import pytest @@ -13,7 +14,8 @@ @pytest.mark.parametrize(('attribute', 'value'), list(product(['output_density', 'output_density_diff', - 'output_specific_energy', 'output_n_photons'], + 'output_specific_energy', 'output_specific_energy_nu', + 'output_n_photons'], ['none', 'last', 'all']))) def test_io_output_conf(attribute, value): o1 = OutputConf() @@ -249,16 +251,16 @@ def test_io_run_conf_mrw(value): r2.read_run_conf(v) assert r2.mrw == r1.mrw -@pytest.mark.parametrize(('value'), [True, False]) -def test_io_run_conf_isrf(value): +def test_io_run_conf_specific_energy_nu_frequencies(): r1 = RunConf() - r1.compute_isrf(value) + r1.set_specific_energy_nu_frequencies(np.logspace(11., 16., 10)) r1.set_n_photons(1, 2) v = virtual_file() r1.write_run_conf(v) r2 = RunConf() r2.read_run_conf(v) - assert r2.isrf == r1.isrf + np.testing.assert_allclose(r2.specific_energy_nu_frequencies, + r1.specific_energy_nu_frequencies) @pytest.mark.parametrize(('value'), [False, True]) def test_io_run_conf_convergence(value): From e873d63d2c0c2944b0d0f7333e6d1045017b62b6 Mon Sep 17 00:00:00 2001 From: Thomas Robitaille Date: Wed, 10 Jun 2026 18:42:41 +0000 Subject: [PATCH 34/38] Drop a redundant phrase from the frequency-resolved specific energy docs --- docs/setup/setup_conf.rst | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/docs/setup/setup_conf.rst b/docs/setup/setup_conf.rst index 9b3682fa..b0128798 100644 --- a/docs/setup/setup_conf.rst +++ b/docs/setup/setup_conf.rst @@ -201,7 +201,7 @@ Frequency-resolved specific energy In addition to the total specific energy absorbed in each cell (``output_specific_energy``), Hyperion can save the specific energy absorbed as a function of frequency. This is controlled in the same way as the other output -quantities -- there is no separate switch to enable it:: +quantities:: m.conf.output.output_specific_energy_nu = 'all' # every iteration # 'last' -> only final iteration From 0879991b9b3a5fa8402480d585c828c11f33a514 Mon Sep 17 00:00:00 2001 From: Thomas Robitaille Date: Wed, 10 Jun 2026 18:43:46 +0000 Subject: [PATCH 35/38] Clarify that specific_energy_nu frequencies are bin centers not edges --- docs/setup/setup_conf.rst | 6 ++++-- hyperion/conf/conf_files.py | 3 ++- 2 files changed, 6 insertions(+), 3 deletions(-) diff --git a/docs/setup/setup_conf.rst b/docs/setup/setup_conf.rst index b0128798..2be2b9db 100644 --- a/docs/setup/setup_conf.rst +++ b/docs/setup/setup_conf.rst @@ -212,7 +212,8 @@ When it is not ``'none'``, two extra arrays are written out: * ``specific_energy_nu`` -- the specific energy absorbed in each cell as a function of frequency, in erg/s/g. * ``specific_energy_nu_frequencies`` -- the frequencies (in Hz) corresponding to - the leading axis of ``specific_energy_nu``. + the leading axis of ``specific_energy_nu``. These are bin centers, not edges, + so this array has exactly the same length as that axis (one entry per bin). This is the *absorbed* (deposited) energy spectrum, not the mean intensity: each contribution is weighted by the dust absorption opacity, so it is proportional @@ -228,7 +229,8 @@ properties:: import numpy as np m.set_specific_energy_nu_frequencies(np.logspace(11., 16., 100)) -Photons are binned to the nearest frequency in log space. This works for all +Photons are binned to the nearest of these frequencies in log space (so the +supplied values act as bin centers). This works for all grid types, including AMR and Voronoi. ``specific_energy_nu`` can be retrieved like other grid quantities, as an array diff --git a/hyperion/conf/conf_files.py b/hyperion/conf/conf_files.py index 6046828c..12071148 100644 --- a/hyperion/conf/conf_files.py +++ b/hyperion/conf/conf_files.py @@ -420,7 +420,8 @@ def set_specific_energy_nu_frequencies(self, frequencies): This is only relevant if ``conf.output.output_specific_energy_nu`` is set to ``'all'`` or ``'last'``. If this method is not called, the frequency grid of the first dust type is used. Photons are binned to the nearest - frequency in log space. + frequency in log space, so the supplied values are bin centers (not + edges) and ``specific_energy_nu`` has one entry per supplied frequency. Parameters ---------- From 1f38f1d6477429a2b0132df9c4621fbd71569212 Mon Sep 17 00:00:00 2001 From: Thomas Robitaille Date: Wed, 10 Jun 2026 18:52:15 +0000 Subject: [PATCH 36/38] Rename specific_energy_nu to specific_energy_spectrum throughout --- docs/setup/setup_conf.rst | 16 +-- hyperion/conf/conf_files.py | 44 ++++----- hyperion/conf/tests/test_conf_io.py | 10 +- hyperion/grid/cartesian_grid.py | 2 +- hyperion/grid/cylindrical_polar_grid.py | 2 +- hyperion/grid/octree_grid.py | 2 +- hyperion/grid/spherical_polar_grid.py | 2 +- hyperion/grid/voronoi_grid.py | 2 +- hyperion/model/model.py | 18 ++-- ...nu.py => test_specific_energy_spectrum.py} | 98 +++++++++---------- src/grid/grid_generic.f90 | 20 ++-- src/grid/grid_physics_3d.f90 | 70 ++++++------- src/grid/grid_propagate_3d.f90 | 12 +-- src/main/settings.f90 | 8 +- src/main/setup_rt.f90 | 26 ++--- src/mpi/mpi_routines.f90 | 8 +- 16 files changed, 170 insertions(+), 170 deletions(-) rename hyperion/model/tests/{test_specific_energy_nu.py => test_specific_energy_spectrum.py} (71%) diff --git a/docs/setup/setup_conf.rst b/docs/setup/setup_conf.rst index 2be2b9db..a0026479 100644 --- a/docs/setup/setup_conf.rst +++ b/docs/setup/setup_conf.rst @@ -203,23 +203,23 @@ In addition to the total specific energy absorbed in each cell a function of frequency. This is controlled in the same way as the other output quantities:: - m.conf.output.output_specific_energy_nu = 'all' # every iteration + m.conf.output.output_specific_energy_spectrum = 'all' # every iteration # 'last' -> only final iteration # 'none' -> not computed or saved (default) When it is not ``'none'``, two extra arrays are written out: -* ``specific_energy_nu`` -- the specific energy absorbed in each cell as a +* ``specific_energy_spectrum`` -- the specific energy absorbed in each cell as a function of frequency, in erg/s/g. -* ``specific_energy_nu_frequencies`` -- the frequencies (in Hz) corresponding to - the leading axis of ``specific_energy_nu``. These are bin centers, not edges, +* ``specific_energy_spectrum_frequencies`` -- the frequencies (in Hz) corresponding to + the leading axis of ``specific_energy_spectrum``. These are bin centers, not edges, so this array has exactly the same length as that axis (one entry per bin). This is the *absorbed* (deposited) energy spectrum, not the mean intensity: each contribution is weighted by the dust absorption opacity, so it is proportional to :math:`\kappa_\nu J_\nu`. To recover the mean intensity :math:`J_\nu` (the radiation field), divide by the dust absorption opacity at each frequency. -Summed over frequency, ``specific_energy_nu`` recovers the total +Summed over frequency, ``specific_energy_spectrum`` recovers the total ``specific_energy``. By default the binning uses the frequency grid of the first dust type. You can @@ -227,16 +227,16 @@ instead provide your own frequency grid (in Hz), independent of the dust properties:: import numpy as np - m.set_specific_energy_nu_frequencies(np.logspace(11., 16., 100)) + m.set_specific_energy_spectrum_frequencies(np.logspace(11., 16., 100)) Photons are binned to the nearest of these frequencies in log space (so the supplied values act as bin centers). This works for all grid types, including AMR and Voronoi. -``specific_energy_nu`` can be retrieved like other grid quantities, as an array +``specific_energy_spectrum`` can be retrieved like other grid quantities, as an array with an extra leading frequency axis:: - senu = m.get_quantities()['specific_energy_nu'] + senu = m.get_quantities()['specific_energy_spectrum'] Advanced parameters ------------------- diff --git a/hyperion/conf/conf_files.py b/hyperion/conf/conf_files.py index 12071148..20cd1b14 100644 --- a/hyperion/conf/conf_files.py +++ b/hyperion/conf/conf_files.py @@ -18,7 +18,7 @@ def __init__(self): self.output_density = 'none' self.output_density_diff = 'none' self.output_specific_energy = 'last' - self.output_specific_energy_nu = 'none' + self.output_specific_energy_spectrum = 'none' self.output_n_photons = 'none' self._freeze() @@ -28,10 +28,10 @@ def read(cls, group): self.output_density = group.attrs['output_density'].decode('utf-8') self.output_density_diff = group.attrs['output_density_diff'].decode('utf-8') self.output_specific_energy = group.attrs['output_specific_energy'].decode('utf-8') - if 'output_specific_energy_nu' in group.attrs: - self.output_specific_energy_nu = group.attrs['output_specific_energy_nu'].decode('utf-8') + if 'output_specific_energy_spectrum' in group.attrs: + self.output_specific_energy_spectrum = group.attrs['output_specific_energy_spectrum'].decode('utf-8') else: - self.output_specific_energy_nu = 'none' + self.output_specific_energy_spectrum = 'none' self.output_n_photons = group.attrs['output_n_photons'].decode('utf-8') return self @@ -39,7 +39,7 @@ def write(self, group): group.attrs['output_density'] = np.bytes_(self.output_density.encode('utf-8')) group.attrs['output_density_diff'] = np.bytes_(self.output_density_diff.encode('utf-8')) group.attrs['output_specific_energy'] = np.bytes_(self.output_specific_energy.encode('utf-8')) - group.attrs['output_specific_energy_nu'] = np.bytes_(self.output_specific_energy_nu.encode('utf-8')) + group.attrs['output_specific_energy_spectrum'] = np.bytes_(self.output_specific_energy_spectrum.encode('utf-8')) group.attrs['output_n_photons'] = np.bytes_(self.output_n_photons.encode('utf-8')) @@ -58,7 +58,7 @@ def __init__(self): self.set_max_reabsorptions(1000000) self.set_pda(False) self.set_mrw(False) - self.specific_energy_nu_frequencies = None + self.specific_energy_spectrum_frequencies = None self.set_convergence(False) self.set_kill_on_absorb(False) @@ -399,41 +399,41 @@ def _read_pda(self, group): def _write_pda(self, group): group.attrs['pda'] = bool2str(self.pda) - def _read_specific_energy_nu_frequencies(self, group): - if 'specific_energy_nu_frequencies' in group: - self.specific_energy_nu_frequencies = np.array(group['specific_energy_nu_frequencies']['nu']) + def _read_specific_energy_spectrum_frequencies(self, group): + if 'specific_energy_spectrum_frequencies' in group: + self.specific_energy_spectrum_frequencies = np.array(group['specific_energy_spectrum_frequencies']['nu']) else: - self.specific_energy_nu_frequencies = None + self.specific_energy_spectrum_frequencies = None - def _write_specific_energy_nu_frequencies(self, group): - if self.specific_energy_nu_frequencies is not None: - group.create_dataset('specific_energy_nu_frequencies', - data=np.array(list(zip(self.specific_energy_nu_frequencies)), + def _write_specific_energy_spectrum_frequencies(self, group): + if self.specific_energy_spectrum_frequencies is not None: + group.create_dataset('specific_energy_spectrum_frequencies', + data=np.array(list(zip(self.specific_energy_spectrum_frequencies)), dtype=[('nu', float)])) - def set_specific_energy_nu_frequencies(self, frequencies): + def set_specific_energy_spectrum_frequencies(self, frequencies): ''' Set the frequency grid onto which the frequency-resolved specific energy - (``specific_energy_nu``) is binned. + (``specific_energy_spectrum``) is binned. - This is only relevant if ``conf.output.output_specific_energy_nu`` is set + This is only relevant if ``conf.output.output_specific_energy_spectrum`` is set to ``'all'`` or ``'last'``. If this method is not called, the frequency grid of the first dust type is used. Photons are binned to the nearest frequency in log space, so the supplied values are bin centers (not - edges) and ``specific_energy_nu`` has one entry per supplied frequency. + edges) and ``specific_energy_spectrum`` has one entry per supplied frequency. Parameters ---------- frequencies : iterable of float - The frequencies (in Hz) onto which to bin ``specific_energy_nu``. + The frequencies (in Hz) onto which to bin ``specific_energy_spectrum``. ''' frequencies = np.asarray(frequencies, dtype=float) if frequencies.ndim != 1 or frequencies.size < 1: raise ValueError("frequencies should be a 1-d array of at least one value") if np.any(frequencies <= 0.): raise ValueError("frequencies should be positive (in Hz)") - self.specific_energy_nu_frequencies = np.sort(frequencies) + self.specific_energy_spectrum_frequencies = np.sort(frequencies) def set_mrw(self, mrw, gamma=1.0, inter_max=1000, warn=True): @@ -764,7 +764,7 @@ def read_run_conf(self, group): # not a class method because inherited self._read_max_reabsorptions(group) self._read_pda(group) self._read_mrw(group) - self._read_specific_energy_nu_frequencies(group) + self._read_specific_energy_spectrum_frequencies(group) self._read_convergence(group) self._read_kill_on_absorb(group) self._read_kill_on_scatter(group) @@ -793,7 +793,7 @@ def write_run_conf(self, group): self._write_max_reabsorptions(group) self._write_pda(group) self._write_mrw(group) - self._write_specific_energy_nu_frequencies(group) + self._write_specific_energy_spectrum_frequencies(group) self._write_convergence(group) self._write_kill_on_absorb(group) self._write_kill_on_scatter(group) diff --git a/hyperion/conf/tests/test_conf_io.py b/hyperion/conf/tests/test_conf_io.py index fb667e72..9cfe7ed5 100644 --- a/hyperion/conf/tests/test_conf_io.py +++ b/hyperion/conf/tests/test_conf_io.py @@ -14,7 +14,7 @@ @pytest.mark.parametrize(('attribute', 'value'), list(product(['output_density', 'output_density_diff', - 'output_specific_energy', 'output_specific_energy_nu', + 'output_specific_energy', 'output_specific_energy_spectrum', 'output_n_photons'], ['none', 'last', 'all']))) def test_io_output_conf(attribute, value): @@ -251,16 +251,16 @@ def test_io_run_conf_mrw(value): r2.read_run_conf(v) assert r2.mrw == r1.mrw -def test_io_run_conf_specific_energy_nu_frequencies(): +def test_io_run_conf_specific_energy_spectrum_frequencies(): r1 = RunConf() - r1.set_specific_energy_nu_frequencies(np.logspace(11., 16., 10)) + r1.set_specific_energy_spectrum_frequencies(np.logspace(11., 16., 10)) r1.set_n_photons(1, 2) v = virtual_file() r1.write_run_conf(v) r2 = RunConf() r2.read_run_conf(v) - np.testing.assert_allclose(r2.specific_energy_nu_frequencies, - r1.specific_energy_nu_frequencies) + np.testing.assert_allclose(r2.specific_energy_spectrum_frequencies, + r1.specific_energy_spectrum_frequencies) @pytest.mark.parametrize(('value'), [False, True]) def test_io_run_conf_convergence(value): diff --git a/hyperion/grid/cartesian_grid.py b/hyperion/grid/cartesian_grid.py index dc9f71cb..b1f00bde 100644 --- a/hyperion/grid/cartesian_grid.py +++ b/hyperion/grid/cartesian_grid.py @@ -279,7 +279,7 @@ def read_quantities(self, group, quantities='all'): # Read in physical quantities if quantities is not None: for quantity in group: - if quantity == 'specific_energy_nu_frequencies': + if quantity == 'specific_energy_spectrum_frequencies': continue # per-frequency metadata, not a per-cell grid quantity if quantities == 'all' or quantity in quantities: array = np.array(group[quantity]) diff --git a/hyperion/grid/cylindrical_polar_grid.py b/hyperion/grid/cylindrical_polar_grid.py index a6a8b7cf..2dc5b149 100644 --- a/hyperion/grid/cylindrical_polar_grid.py +++ b/hyperion/grid/cylindrical_polar_grid.py @@ -307,7 +307,7 @@ def read_quantities(self, group, quantities='all'): # Read in physical quantities if quantities is not None: for quantity in group: - if quantity == 'specific_energy_nu_frequencies': + if quantity == 'specific_energy_spectrum_frequencies': continue # per-frequency metadata, not a per-cell grid quantity if quantities == 'all' or quantity in quantities: array = np.array(group[quantity]) diff --git a/hyperion/grid/octree_grid.py b/hyperion/grid/octree_grid.py index 31696047..38485b16 100644 --- a/hyperion/grid/octree_grid.py +++ b/hyperion/grid/octree_grid.py @@ -372,7 +372,7 @@ def read_quantities(self, group, quantities='all'): # Read in physical quantities if quantities is not None: for quantity in group: - if quantity == 'specific_energy_nu_frequencies': + if quantity == 'specific_energy_spectrum_frequencies': continue # per-frequency metadata, not a per-cell grid quantity if quantities == 'all' or quantity in quantities: array = np.array(group[quantity]) diff --git a/hyperion/grid/spherical_polar_grid.py b/hyperion/grid/spherical_polar_grid.py index 9fa45bdc..e6d925d4 100644 --- a/hyperion/grid/spherical_polar_grid.py +++ b/hyperion/grid/spherical_polar_grid.py @@ -317,7 +317,7 @@ def read_quantities(self, group, quantities='all'): # Read in physical quantities if quantities is not None: for quantity in group: - if quantity == 'specific_energy_nu_frequencies': + if quantity == 'specific_energy_spectrum_frequencies': continue # per-frequency metadata, not a per-cell grid quantity if quantities == 'all' or quantity in quantities: array = np.array(group[quantity]) diff --git a/hyperion/grid/voronoi_grid.py b/hyperion/grid/voronoi_grid.py index 81b48d2d..d15a70a3 100644 --- a/hyperion/grid/voronoi_grid.py +++ b/hyperion/grid/voronoi_grid.py @@ -396,7 +396,7 @@ def read_quantities(self, group, quantities='all'): # Read in physical quantities if quantities is not None: for quantity in group: - if quantity == 'specific_energy_nu_frequencies': + if quantity == 'specific_energy_spectrum_frequencies': continue # per-frequency metadata, not a per-cell grid quantity if quantities == 'all' or quantity in quantities: array = np.array(group[quantity]) diff --git a/hyperion/model/model.py b/hyperion/model/model.py index a9143c03..7701d0d5 100644 --- a/hyperion/model/model.py +++ b/hyperion/model/model.py @@ -223,7 +223,7 @@ def use_geometry(self, filename): # Close the file f.close() - def use_quantities(self, filename, quantities=['density', 'specific_energy', 'specific_energy_nu'], + def use_quantities(self, filename, quantities=['density', 'specific_energy', 'specific_energy_spectrum'], use_minimum_specific_energy=True, use_dust=True, copy=True, only_initial=False): ''' @@ -297,12 +297,12 @@ def use_quantities(self, filename, quantities=['density', 'specific_energy', 'sp else: quantities_path['specific_energy'] = last_iteration - if 'specific_energy_nu' in quantities: + if 'specific_energy_spectrum' in quantities: if only_initial or last_iteration is None: - if 'specific_energy_nu' in f['/Input/Grid/Quantities']: - quantities_path['specific_energy_nu'] = '/Input/Grid/Quantities' - elif 'specific_energy_nu' in f[last_iteration]: - quantities_path['specific_energy_nu'] = last_iteration + if 'specific_energy_spectrum' in f['/Input/Grid/Quantities']: + quantities_path['specific_energy_spectrum'] = '/Input/Grid/Quantities' + elif 'specific_energy_spectrum' in f[last_iteration]: + quantities_path['specific_energy_spectrum'] = last_iteration # Minimum specific energy if use_minimum_specific_energy: @@ -320,9 +320,9 @@ def use_quantities(self, filename, quantities=['density', 'specific_energy', 'sp if 'specific_energy' in f['/Grid/Quantities']: quantities_path['specific_energy'] = '/Grid/Quantities' - if 'specific_energy_nu' in quantities: - if 'specific_energy_nu' in f['/Grid/Quantities']: - quantities_path['specific_energy_nu'] = '/Grid/Quantities' + if 'specific_energy_spectrum' in quantities: + if 'specific_energy_spectrum' in f['/Grid/Quantities']: + quantities_path['specific_energy_spectrum'] = '/Grid/Quantities' # Minimum specific energy if use_minimum_specific_energy: diff --git a/hyperion/model/tests/test_specific_energy_nu.py b/hyperion/model/tests/test_specific_energy_spectrum.py similarity index 71% rename from hyperion/model/tests/test_specific_energy_nu.py rename to hyperion/model/tests/test_specific_energy_spectrum.py index d7483da0..f609e5e2 100644 --- a/hyperion/model/tests/test_specific_energy_nu.py +++ b/hyperion/model/tests/test_specific_energy_spectrum.py @@ -25,13 +25,13 @@ def _read_dataset(filename, suffix): return np.asarray(f[matches[-1]][()], dtype=float) -def _assert_nu_sums_to_specific_energy(filename): - # Summed over frequency, specific_energy_nu must reproduce specific_energy. +def _assert_spectrum_sums_to_specific_energy(filename): + # 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_nu), + # (which only affects specific_energy, not the unclamped specific_energy_spectrum), # so we compare only cells that were genuinely heated above that floor. se = _read_dataset(filename, '/specific_energy') - se_nu = _read_dataset(filename, '/specific_energy_nu') + se_nu = _read_dataset(filename, '/specific_energy_spectrum') assert se_nu is not None nu_sum = se_nu.sum(axis=0) floor = se.min() @@ -40,7 +40,7 @@ def _assert_nu_sums_to_specific_energy(filename): np.testing.assert_allclose(nu_sum[heated], se[heated], rtol=1.e-6) -def _cartesian_model(output_specific_energy_nu='last', n_cells=2, n_photons=100000, +def _cartesian_model(output_specific_energy_spectrum='last', n_cells=2, n_photons=100000, density=1.e-18, frequencies=None): m = Model() edges = np.linspace(-1., 1., n_cells + 1) @@ -54,20 +54,20 @@ def _cartesian_model(output_specific_energy_nu='last', n_cells=2, n_photons=1000 m.set_n_photons(initial=n_photons, imaging=0) m.set_seed(-12345) m.conf.output.output_specific_energy = 'last' - m.conf.output.output_specific_energy_nu = output_specific_energy_nu + m.conf.output.output_specific_energy_spectrum = output_specific_energy_spectrum if frequencies is not None: - m.set_specific_energy_nu_frequencies(frequencies) + m.set_specific_energy_spectrum_frequencies(frequencies) return m @pytest.mark.requires_hyperion_binaries -def test_specific_energy_nu_does_not_change_specific_energy(tmpdir): +def test_specific_energy_spectrum_does_not_change_specific_energy(tmpdir): # Computing the frequency-resolved specific energy must not perturb the # ordinary specific-energy (temperature) calculation. This guards against the # regression where specific_energy_sum accumulation was gated on the option. se = {} for output in ('none', 'last'): - m = _cartesian_model(output_specific_energy_nu=output) + m = _cartesian_model(output_specific_energy_spectrum=output) m.write(tmpdir.join(random_id()).strpath) out = m.run(tmpdir.join(random_id()).strpath) se[output] = _read_dataset(out.filename, '/specific_energy') @@ -77,29 +77,29 @@ def test_specific_energy_nu_does_not_change_specific_energy(tmpdir): @pytest.mark.requires_hyperion_binaries -def test_specific_energy_nu_sums_to_specific_energy(tmpdir): - # specific_energy_nu, summed over frequency, must reproduce specific_energy. +def test_specific_energy_spectrum_sums_to_specific_energy(tmpdir): + # specific_energy_spectrum, summed over frequency, must reproduce specific_energy. m = _cartesian_model('last') m.write(tmpdir.join(random_id()).strpath) out = m.run(tmpdir.join(random_id()).strpath) - _assert_nu_sums_to_specific_energy(out.filename) + _assert_spectrum_sums_to_specific_energy(out.filename) @pytest.mark.requires_hyperion_binaries -def test_specific_energy_nu_frequencies_written(tmpdir): - # specific_energy_nu_frequencies is a per-frequency (1-D) quantity, not +def test_specific_energy_spectrum_frequencies_written(tmpdir): + # specific_energy_spectrum_frequencies is a per-frequency (1-D) quantity, not # per-cell. This guards against the regression where it was written as a grid # array, which segfaulted whenever n_nu was smaller than n_cells. m = _cartesian_model('last', n_cells=8) # 512 cells, dust has 2 frequencies m.write(tmpdir.join(random_id()).strpath) out = m.run(tmpdir.join(random_id()).strpath) - bins = _read_dataset(out.filename, '/specific_energy_nu_frequencies') + bins = _read_dataset(out.filename, '/specific_energy_spectrum_frequencies') assert bins is not None assert bins.shape == (2,) # == number of dust frequencies, not n_cells @pytest.mark.requires_hyperion_binaries -def test_specific_energy_nu_dustless_model_runs(tmpdir): +def test_specific_energy_spectrum_dustless_model_runs(tmpdir): # A model with no dust must still run (and image) without crashing. This # guards against the regression where the instrumentation dereferenced the # first dust type unconditionally, segfaulting dustless models. @@ -119,8 +119,8 @@ def test_specific_energy_nu_dustless_model_runs(tmpdir): @pytest.mark.requires_hyperion_binaries -def test_specific_energy_nu_amr(tmpdir): - # specific_energy_nu must work on AMR grids: it should run and, summed over +def test_specific_energy_spectrum_amr(tmpdir): + # specific_energy_spectrum must work on AMR grids: it should run and, summed over # frequency, reproduce specific_energy in every cell. m = Model() amr = AMRGrid() @@ -140,14 +140,14 @@ def test_specific_energy_nu_amr(tmpdir): m.set_n_photons(initial=100000, imaging=0) m.set_seed(-9) m.conf.output.output_specific_energy = 'last' - m.conf.output.output_specific_energy_nu = '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_nu_sums_to_specific_energy(out.filename) + _assert_spectrum_sums_to_specific_energy(out.filename) @pytest.mark.requires_hyperion_binaries -def test_specific_energy_nu_custom_frequency_grid(tmpdir): +def test_specific_energy_spectrum_custom_frequency_grid(tmpdir): # A user-specified frequency grid (independent of the dust frequency grid) # should be used for the output bins, and energy must still be conserved # when summing over frequency. @@ -155,40 +155,40 @@ def test_specific_energy_nu_custom_frequency_grid(tmpdir): m = _cartesian_model('last', density=1.e-16, frequencies=frequencies) m.write(tmpdir.join(random_id()).strpath) out = m.run(tmpdir.join(random_id()).strpath) - bins = _read_dataset(out.filename, '/specific_energy_nu_frequencies') + bins = _read_dataset(out.filename, '/specific_energy_spectrum_frequencies') np.testing.assert_allclose(bins, frequencies, rtol=1.e-6) - se_nu = _read_dataset(out.filename, '/specific_energy_nu') + se_nu = _read_dataset(out.filename, '/specific_energy_spectrum') assert se_nu.shape[0] == len(frequencies) - _assert_nu_sums_to_specific_energy(out.filename) + _assert_spectrum_sums_to_specific_energy(out.filename) -def test_specific_energy_nu_frequencies_roundtrip(tmpdir): +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 frequencies = np.logspace(10., 16., 8) conf = RunConf() - conf.set_specific_energy_nu_frequencies(frequencies) + conf.set_specific_energy_spectrum_frequencies(frequencies) with h5py.File(tmpdir.join('conf.h5').strpath, 'w') as f: - conf._write_specific_energy_nu_frequencies(f) + conf._write_specific_energy_spectrum_frequencies(f) conf2 = RunConf() - conf2._read_specific_energy_nu_frequencies(f) - np.testing.assert_allclose(conf2.specific_energy_nu_frequencies, frequencies) + conf2._read_specific_energy_spectrum_frequencies(f) + np.testing.assert_allclose(conf2.specific_energy_spectrum_frequencies, frequencies) -def test_specific_energy_nu_frequencies_validation(): +def test_specific_energy_spectrum_frequencies_validation(): from ...conf.conf_files import RunConf conf = RunConf() with pytest.raises(ValueError): - conf.set_specific_energy_nu_frequencies([[1.e10, 1.e12], [1.e14, 1.e16]]) + conf.set_specific_energy_spectrum_frequencies([[1.e10, 1.e12], [1.e14, 1.e16]]) with pytest.raises(ValueError): - conf.set_specific_energy_nu_frequencies([1.e10, -1.e12]) + conf.set_specific_energy_spectrum_frequencies([1.e10, -1.e12]) # Default: no custom grid set - assert conf.specific_energy_nu_frequencies is None + assert conf.specific_energy_spectrum_frequencies is None @pytest.mark.requires_hyperion_binaries -def test_specific_energy_nu_voronoi(tmpdir): - # specific_energy_nu must also work on (unstructured) Voronoi grids, where +def test_specific_energy_spectrum_voronoi(tmpdir): + # specific_energy_spectrum must also work on (unstructured) Voronoi grids, where # cells are a flat list rather than a structured grid. rng = np.random.RandomState(12345) n = 30 @@ -205,25 +205,25 @@ def test_specific_energy_nu_voronoi(tmpdir): 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_nu = '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_nu_sums_to_specific_energy(out.filename) + _assert_spectrum_sums_to_specific_energy(out.filename) @pytest.mark.requires_hyperion_binaries -def test_specific_energy_nu_get_quantities(tmpdir): - # specific_energy_nu should be retrievable through get_quantities, the +def test_specific_energy_spectrum_get_quantities(tmpdir): + # specific_energy_spectrum should be retrievable through get_quantities, the # per-frequency frequencies metadata should not pollute the grid quantities, # and the retrieved array should still conserve energy. m = _cartesian_model('last', density=1.e-16) m.write(tmpdir.join(random_id()).strpath) out = m.run(tmpdir.join(random_id()).strpath) g = out.get_quantities() - assert 'specific_energy_nu' in g.quantities - assert 'specific_energy_nu_frequencies' not in g.quantities + assert 'specific_energy_spectrum' in g.quantities + assert 'specific_energy_spectrum_frequencies' not in g.quantities se = np.array(g.quantities['specific_energy']) # (dust, nz, ny, nx) - se_nu = np.array(g.quantities['specific_energy_nu']) # (nu, dust, nz, ny, nx) + se_nu = np.array(g.quantities['specific_energy_spectrum']) # (nu, dust, nz, ny, nx) assert se_nu.shape[0] == 2 # number of dust frequencies nu_sum = se_nu.sum(axis=0) heated = se > se.min() * (1. + 1.e-6) @@ -231,8 +231,8 @@ def test_specific_energy_nu_get_quantities(tmpdir): @pytest.mark.requires_hyperion_binaries -def test_specific_energy_nu_yt_export_skips_frequency_resolved(tmpdir): - # The frequency-resolved specific_energy_nu is not a per-cell scalar, so it +def test_specific_energy_spectrum_yt_export_skips_frequency_resolved(tmpdir): + # The frequency-resolved specific_energy_spectrum is not a per-cell scalar, so it # must be skipped (not exported as a malformed field) when converting to yt, # while the ordinary scalar quantities are still exported. pytest.importorskip("yt") @@ -240,15 +240,15 @@ def test_specific_energy_nu_yt_export_skips_frequency_resolved(tmpdir): m.write(tmpdir.join(random_id()).strpath) out = m.run(tmpdir.join(random_id()).strpath) g = out.get_quantities() - assert 'specific_energy_nu' in g.quantities + assert 'specific_energy_spectrum' in g.quantities ds = g.to_yt() fields = [str(f) for f in ds.field_list] - assert not any('specific_energy_nu' in f for f in fields) + assert not any('specific_energy_spectrum' in f for f in fields) assert any(f.endswith("'density')") for f in fields) @pytest.mark.requires_hyperion_binaries -def test_specific_energy_nu_mpi_matches_serial(tmpdir): +def test_specific_energy_spectrum_mpi_matches_serial(tmpdir): # The saved spectrum must not depend on the number of MPI processes. We # compare the total (which is conserved, hence robust to Monte Carlo noise) # between a serial run and a 2-process MPI run. @@ -261,6 +261,6 @@ def test_specific_energy_nu_mpi_matches_serial(tmpdir): out_serial = m.run(tmpdir.join(random_id()).strpath) out_mpi = m.run(tmpdir.join(random_id()).strpath, mpi=True, n_processes=2) - total_serial = np.nansum(_read_dataset(out_serial.filename, '/specific_energy_nu')) - total_mpi = np.nansum(_read_dataset(out_mpi.filename, '/specific_energy_nu')) + 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) diff --git a/src/grid/grid_generic.f90 b/src/grid/grid_generic.f90 index 786296b5..c7ec2641 100644 --- a/src/grid/grid_generic.f90 +++ b/src/grid/grid_generic.f90 @@ -6,8 +6,8 @@ module grid_generic use grid_io, only : write_grid_3d, write_grid_4d, write_grid_5d use grid_geometry, only : geo - use grid_physics, only : n_photons, last_photon_id, specific_energy_sum, specific_energy_sum_nu, specific_energy, specific_energy_nu, nu_bins, density, density_original - use settings, only : output_n_photons, output_specific_energy, output_specific_energy_nu, output_density, output_density_diff, physics_io_type, compute_specific_energy_nu + use grid_physics, only : n_photons, last_photon_id, specific_energy_sum, specific_energy_sum_spectrum, specific_energy, specific_energy_spectrum, nu_bins, density, density_original + use settings, only : output_n_photons, output_specific_energy, output_specific_energy_spectrum, output_density, output_density_diff, physics_io_type, compute_specific_energy_spectrum implicit none save @@ -23,7 +23,7 @@ subroutine grid_reset_energy() if(allocated(n_photons)) n_photons = 0 if(allocated(last_photon_id)) last_photon_id = 0 if(allocated(specific_energy_sum)) specific_energy_sum = 0._dp - if(allocated(specific_energy_sum_nu)) specific_energy_sum_nu = 0._dp + if(allocated(specific_energy_sum_spectrum)) specific_energy_sum_spectrum = 0._dp end subroutine grid_reset_energy subroutine output_grid(group, iter, n_iter) @@ -65,25 +65,25 @@ subroutine output_grid(group, iter, n_iter) ! WRITE THE FREQUENCY-RESOLVED SPECIFIC ENERGY - if (compute_specific_energy_nu) then + if (compute_specific_energy_spectrum) then - if(trim(output_specific_energy_nu)=='all' .or. (trim(output_specific_energy_nu)=='last'.and.iter==n_iter)) then + if(trim(output_specific_energy_spectrum)=='all' .or. (trim(output_specific_energy_spectrum)=='last'.and.iter==n_iter)) then ! The frequencies are a per-frequency quantity, not per-cell, so they ! are written as a plain 1-D dataset rather than a grid array. - call mp_write_array(group, 'specific_energy_nu_frequencies', nu_bins) + call mp_write_array(group, 'specific_energy_spectrum_frequencies', nu_bins) - if(allocated(specific_energy_nu)) then + if(allocated(specific_energy_spectrum)) then select case(physics_io_type) case(sp) - call write_grid_5d(group, 'specific_energy_nu', real(specific_energy_nu, sp), geo) + call write_grid_5d(group, 'specific_energy_spectrum', real(specific_energy_spectrum, sp), geo) case(dp) - call write_grid_5d(group, 'specific_energy_nu', real(specific_energy_nu, dp), geo) + call write_grid_5d(group, 'specific_energy_spectrum', real(specific_energy_spectrum, dp), geo) case default call error("output_grid","unexpected value of physics_io_type (should be sp or dp)") end select else - call warn("output_grid","specific_energy_nu array is not allocated") + call warn("output_grid","specific_energy_spectrum array is not allocated") end if end if diff --git a/src/grid/grid_physics_3d.f90 b/src/grid/grid_physics_3d.f90 index e48a5958..26578fc6 100644 --- a/src/grid/grid_physics_3d.f90 +++ b/src/grid/grid_physics_3d.f90 @@ -36,14 +36,14 @@ module grid_physics integer(idp),allocatable, public :: n_photons(:) integer(idp),allocatable, public :: last_photon_id(:) real(dp),allocatable, public :: specific_energy(:,:) - real(dp),allocatable, public :: specific_energy_nu(:,:,:) + real(dp),allocatable, public :: specific_energy_spectrum(:,:,:) real(dp),allocatable, public :: specific_energy_sum(:,:) - real(dp),allocatable, public :: specific_energy_sum_nu(:,:,:) + real(dp),allocatable, public :: specific_energy_sum_spectrum(:,:,:) real(dp),allocatable, public :: nu_bins(:) real(dp),allocatable, public :: log_nu_bins(:) real(dp),allocatable, public :: specific_energy_additional(:,:) - real(dp),allocatable, public :: specific_energy_additional_nu(:,:,:) + real(dp),allocatable, public :: specific_energy_additional_spectrum(:,:,:) real(dp),allocatable, public :: energy_abs_tot(:) real(dp),allocatable, public :: minimum_specific_energy(:) @@ -99,25 +99,25 @@ integer function select_dust_specific_energy_rho(icell) result(id_select) id_select = sample_pdf(absorption) end function select_dust_specific_energy_rho - subroutine setup_grid_physics(group, use_mrw, use_pda, compute_specific_energy_nu) + subroutine setup_grid_physics(group, use_mrw, use_pda, compute_specific_energy_spectrum) implicit none integer(hid_t),intent(in) :: group - logical,intent(in) :: use_mrw, use_pda, compute_specific_energy_nu + logical,intent(in) :: use_mrw, use_pda, compute_specific_energy_spectrum integer :: n_nu_bins - ! specific_energy_nu is binned onto a user-specified frequency grid if one was given, + ! 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_nu_frequencies)) then - n_nu_bins = size(specific_energy_nu_frequencies) + if (allocated(specific_energy_spectrum_frequencies)) then + n_nu_bins = size(specific_energy_spectrum_frequencies) else n_nu_bins = d(1)%n_nu end if ! Density allocate(density(geo%n_cells, n_dust)) allocate(specific_energy(geo%n_cells, n_dust)) - allocate(specific_energy_nu(geo%n_cells,n_dust,n_nu_bins)) + allocate(specific_energy_spectrum(geo%n_cells,n_dust,n_nu_bins)) if(n_dust > 0) then @@ -175,10 +175,10 @@ subroutine setup_grid_physics(group, use_mrw, use_pda, compute_specific_energy_n specific_energy(:, id) = 0. end where - if (compute_specific_energy_nu) then + if (compute_specific_energy_spectrum) then do idx=1,n_nu_bins where(.not.geo%mask) - specific_energy_nu(:, id, idx) = 0. + specific_energy_spectrum(:, id, idx) = 0. endwhere end do end if @@ -188,21 +188,21 @@ subroutine setup_grid_physics(group, use_mrw, use_pda, compute_specific_energy_n if(trim(specific_energy_type) == 'additional') then allocate(specific_energy_additional(geo%n_cells, n_dust)) - if (compute_specific_energy_nu) then - allocate(specific_energy_additional_nu(geo%n_cells,n_dust,n_nu_bins)) + if (compute_specific_energy_spectrum) then + allocate(specific_energy_additional_spectrum(geo%n_cells,n_dust,n_nu_bins)) end if ! We store a copy of the initial specific energy in a separate ! array, and we set the specific energy to the minimum specific ! 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_nu = specific_energy_nu + specific_energy_additional_spectrum = specific_energy_spectrum do id=1,n_dust specific_energy(:,id) = minimum_specific_energy(id) - if (compute_specific_energy_nu) then + if (compute_specific_energy_spectrum) then do idx=1,n_nu_bins - specific_energy_nu(:,id,idx) = minimum_specific_energy(id) + specific_energy_spectrum(:,id,idx) = minimum_specific_energy(id) end do end if @@ -220,9 +220,9 @@ subroutine setup_grid_physics(group, use_mrw, use_pda, compute_specific_energy_n do id=1,n_dust specific_energy(:,id) = minimum_specific_energy(id) - if (compute_specific_energy_nu) then + if (compute_specific_energy_spectrum) then do idx = 1,n_nu_bins - specific_energy_nu(:,id,idx) = minimum_specific_energy(id) + specific_energy_spectrum(:,id,idx) = minimum_specific_energy(id) end do end if @@ -240,16 +240,16 @@ subroutine setup_grid_physics(group, use_mrw, use_pda, compute_specific_energy_n specific_energy_sum = 0._dp ! Set up basics for the frequency-resolved specific energy - allocate(specific_energy_sum_nu(geo%n_cells, n_dust, n_nu_bins)) - specific_energy_sum_nu = 0._dp + allocate(specific_energy_sum_spectrum(geo%n_cells, n_dust, n_nu_bins)) + specific_energy_sum_spectrum = 0._dp ! 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 ! log-frequency space (see grid_propagate). - if (compute_specific_energy_nu) then + if (compute_specific_energy_spectrum) then allocate(nu_bins(n_nu_bins)) - if (allocated(specific_energy_nu_frequencies)) then - nu_bins = specific_energy_nu_frequencies + if (allocated(specific_energy_spectrum_frequencies)) then + nu_bins = specific_energy_spectrum_frequencies else do idx=1,n_nu_bins nu_bins(idx) = d(1)%nu(idx) @@ -327,7 +327,7 @@ subroutine sublimate_dust() integer :: reset integer :: n_nu_bins - n_nu_bins = size(specific_energy_nu, 3) + n_nu_bins = size(specific_energy_spectrum, 3) reset = 0 @@ -342,9 +342,9 @@ subroutine sublimate_dust() specific_energy(ic, id) = minimum_specific_energy(id) reset = reset + 1 - if (compute_specific_energy_nu) then + if (compute_specific_energy_spectrum) then do idx=1,n_nu_bins - specific_energy_nu(ic,id,idx) = minimum_specific_energy(id) + specific_energy_spectrum(ic,id,idx) = minimum_specific_energy(id) end do end if @@ -365,9 +365,9 @@ subroutine sublimate_dust() reset = reset + 1 - if (compute_specific_energy_nu) then + if (compute_specific_energy_spectrum) then do idx=1,n_nu_bins - specific_energy_nu(ic,id,idx) = minimum_specific_energy(id) + specific_energy_spectrum(ic,id,idx) = minimum_specific_energy(id) end do end if @@ -385,9 +385,9 @@ subroutine sublimate_dust() end if - if (compute_specific_energy_nu) then + if (compute_specific_energy_spectrum) then do idx=1,n_nu_bins - specific_energy_nu(ic,id,idx) = minimum_specific_energy(id) + specific_energy_spectrum(ic,id,idx) = minimum_specific_energy(id) end do end if @@ -414,16 +414,16 @@ subroutine update_energy_abs(scale) integer :: id,idx integer :: n_nu_bins - n_nu_bins = size(specific_energy_nu, 3) + n_nu_bins = size(specific_energy_spectrum, 3) if(main_process()) write(*,'(" [grid_physics] updating energy_abs")') do id=1,n_dust specific_energy(:,id) = specific_energy_sum(:,id) * scale / geo%volume - if (compute_specific_energy_nu) then + if (compute_specific_energy_spectrum) then do idx=1,n_nu_bins - specific_energy_nu(:,id,idx) = specific_energy_sum_nu(:,id,idx) * scale/geo%volume + specific_energy_spectrum(:,id,idx) = specific_energy_sum_spectrum(:,id,idx) * scale/geo%volume end do end if @@ -444,8 +444,8 @@ subroutine update_energy_abs(scale) specific_energy = specific_energy + specific_energy_additional - if (compute_specific_energy_nu) then - specific_energy_nu = specific_energy_nu + specific_energy_additional_nu + if (compute_specific_energy_spectrum) then + specific_energy_spectrum = specific_energy_spectrum + specific_energy_additional_spectrum end if end if diff --git a/src/grid/grid_propagate_3d.f90 b/src/grid/grid_propagate_3d.f90 index 3ccac335..3769704c 100644 --- a/src/grid/grid_propagate_3d.f90 +++ b/src/grid/grid_propagate_3d.f90 @@ -5,10 +5,10 @@ module grid_propagate use type_grid_cell use dust_main, only : n_dust use grid_geometry, only : escaped, find_wall, in_correct_cell, next_cell, opposite_wall - use grid_physics, only : specific_energy_sum, specific_energy_sum_nu, log_nu_bins, density, n_photons, last_photon_id + use grid_physics, only : specific_energy_sum, specific_energy_sum_spectrum, log_nu_bins, density, n_photons, last_photon_id use sources use counters - use settings, only : frac_check => propagation_check_frequency, compute_specific_energy_nu + use settings, only : frac_check => propagation_check_frequency, compute_specific_energy_spectrum implicit none save @@ -133,15 +133,15 @@ subroutine grid_integrate(p,tau_required,tau_achieved) tau_achieved = tau_achieved + tau_cell - if (compute_specific_energy_nu) idx = minloc(abs(log_nu_bins - log10(p%nu)), DIM=1) + 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) = & & specific_energy_sum(p%icell%ic, id) + tmin * p%current_kappa(id) * p%energy - if (compute_specific_energy_nu) then - specific_energy_sum_nu(p%icell%ic, id, idx) = & - & specific_energy_sum_nu(p%icell%ic, id, idx) + tmin * p%current_kappa(id) * p%energy + if (compute_specific_energy_spectrum) then + specific_energy_sum_spectrum(p%icell%ic, id, idx) = & + & specific_energy_sum_spectrum(p%icell%ic, id, idx) + tmin * p%current_kappa(id) * p%energy end if end if end do diff --git a/src/main/settings.f90 b/src/main/settings.f90 index faf0ccd7..5fd3bf65 100644 --- a/src/main/settings.f90 +++ b/src/main/settings.f90 @@ -20,7 +20,7 @@ module settings integer(idp),public :: n_last_photons_dust = 0 integer(idp),public :: n_raytracing_photons_sources = 0 integer(idp),public :: n_raytracing_photons_dust = 0 - logical,public :: use_raytracing, use_mrw, use_pda, compute_specific_energy_nu + logical,public :: use_raytracing, use_mrw, use_pda, compute_specific_energy_spectrum logical, public :: kill_on_absorb, kill_on_scatter real(dp),public :: mrw_gamma logical, public :: forced_first_interaction @@ -30,16 +30,16 @@ module settings real(dp) :: monochromatic_energy_threshold real(dp),allocatable :: frequencies(:) - ! Optional user-specified frequency grid for specific_energy_nu. If not + ! Optional user-specified frequency grid for specific_energy_spectrum. If not ! allocated, the frequency grid of the first dust type is used instead. - real(dp),allocatable :: specific_energy_nu_frequencies(:) + real(dp),allocatable :: specific_energy_spectrum_frequencies(:) integer :: physics_io_type character(len=4) :: output_density character(len=4) :: output_density_diff character(len=4) :: output_specific_energy - character(len=4) :: output_specific_energy_nu + character(len=4) :: output_specific_energy_spectrum character(len=4) :: output_n_photons logical :: check_convergence = .false. diff --git a/src/main/setup_rt.f90 b/src/main/setup_rt.f90 index 73ef5641..9bee928c 100644 --- a/src/main/setup_rt.f90 +++ b/src/main/setup_rt.f90 @@ -75,28 +75,28 @@ subroutine setup_initial(input_handle) call mp_read_keyword(input_handle, '/', 'pda', use_pda) call mp_read_keyword(input_handle, '/', 'mrw', use_mrw) ! The frequency-resolved specific energy is computed whenever it is to be - ! output, i.e. unless output_specific_energy_nu is 'none'. This is read here + ! output, i.e. unless output_specific_energy_spectrum is 'none'. This is read here ! (ahead of the main OUTPUT block below) because it is needed to set up the ! grid physics arrays. g_output = mp_open_group(input_handle, '/Output') - if(mp_exists_keyword(g_output, '.', 'output_specific_energy_nu')) then - call mp_read_keyword(g_output, '.', 'output_specific_energy_nu', output_specific_energy_nu) + if(mp_exists_keyword(g_output, '.', 'output_specific_energy_spectrum')) then + call mp_read_keyword(g_output, '.', 'output_specific_energy_spectrum', output_specific_energy_spectrum) else - output_specific_energy_nu = 'none' + output_specific_energy_spectrum = 'none' end if call mp_close_group(g_output) - if(trim(output_specific_energy_nu).ne.'all' & - & .and.trim(output_specific_energy_nu).ne.'last' & - & .and.trim(output_specific_energy_nu).ne.'none') & - & call error("setup_initial","output_specific_energy_nu should be one of all/last/none") + if(trim(output_specific_energy_spectrum).ne.'all' & + & .and.trim(output_specific_energy_spectrum).ne.'last' & + & .and.trim(output_specific_energy_spectrum).ne.'none') & + & call error("setup_initial","output_specific_energy_spectrum should be one of all/last/none") - compute_specific_energy_nu = trim(output_specific_energy_nu).ne.'none' + compute_specific_energy_spectrum = trim(output_specific_energy_spectrum).ne.'none' ! Optional user-specified frequency grid (otherwise the dust grid is used) - if(compute_specific_energy_nu) then - if(mp_path_exists(input_handle, 'specific_energy_nu_frequencies')) then - call mp_table_read_column_auto(input_handle, 'specific_energy_nu_frequencies', 'nu', specific_energy_nu_frequencies) + if(compute_specific_energy_spectrum) then + if(mp_path_exists(input_handle, 'specific_energy_spectrum_frequencies')) then + call mp_table_read_column_auto(input_handle, 'specific_energy_spectrum_frequencies', 'nu', specific_energy_spectrum_frequencies) end if end if @@ -198,7 +198,7 @@ subroutine setup_initial(input_handle) call mp_close_group(g_geometry) g_physics = mp_open_group(input_handle, '/Grid/Quantities') - call setup_grid_physics(g_physics, use_mrw, use_pda, compute_specific_energy_nu) + call setup_grid_physics(g_physics, use_mrw, use_pda, compute_specific_energy_spectrum) call mp_close_group(g_physics) call mp_read_keyword(input_handle, '/', 'physics_io_bytes', physics_io_bytes) diff --git a/src/mpi/mpi_routines.f90 b/src/mpi/mpi_routines.f90 index 6adbf7e0..bb585432 100644 --- a/src/mpi/mpi_routines.f90 +++ b/src/mpi/mpi_routines.f90 @@ -281,12 +281,12 @@ subroutine mp_collect_physical_arrays() end if if(main_process()) then - allocate(tmp_3d(size(specific_energy_sum_nu,1),size(specific_energy_sum_nu,2),size(specific_energy_sum_nu,3))) - call mpi_reduce(specific_energy_sum_nu, tmp_3d, size(specific_energy_sum_nu), mpi_real8, mpi_sum, rank_main, mpi_comm_world, ierr) - specific_energy_sum_nu = tmp_3d + 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_nu, dummy_dp, size(specific_energy_sum_nu), mpi_real8, mpi_sum, rank_main, mpi_comm_world, ierr) + 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 if(allocated(n_photons)) then From 01c272f7e99e20c39a05418751cf2641a3e2be7e Mon Sep 17 00:00:00 2001 From: Thomas Robitaille Date: Wed, 10 Jun 2026 19:18:41 +0000 Subject: [PATCH 37/38] Warn in the docs that the spectrum end bins collect all out-of-range flux --- docs/setup/setup_conf.rst | 9 +++++++++ 1 file changed, 9 insertions(+) diff --git a/docs/setup/setup_conf.rst b/docs/setup/setup_conf.rst index a0026479..21c7747c 100644 --- a/docs/setup/setup_conf.rst +++ b/docs/setup/setup_conf.rst @@ -233,6 +233,15 @@ Photons are binned to the nearest of these frequencies in log space (so the supplied values act as bin centers). This works for all grid types, including AMR and Voronoi. +.. warning:: The binning has no outer edges: the first and last bins collect + *all* photons below and above the outermost bin centers + respectively, however far away in frequency. This keeps the spectrum + energy-conserving (no photons are dropped), but it means a grid that + does not span the full frequency range will pile up out-of-range + flux in its end bins. When supplying a custom grid, make sure it + brackets the full range of frequencies present in the simulation + (the default dust-based grid already does this). + ``specific_energy_spectrum`` can be retrieved like other grid quantities, as an array with an extra leading frequency axis:: From 7906578da9be9063f777d2667cf1a51f6fd799bd Mon Sep 17 00:00:00 2001 From: Thomas Robitaille Date: Wed, 10 Jun 2026 19:40:18 +0000 Subject: [PATCH 38/38] Skip the MPI spectrum test when the MPI build of the binaries is unavailable --- hyperion/model/tests/test_specific_energy_spectrum.py | 6 +++++- 1 file changed, 5 insertions(+), 1 deletion(-) diff --git a/hyperion/model/tests/test_specific_energy_spectrum.py b/hyperion/model/tests/test_specific_energy_spectrum.py index f609e5e2..78a639b4 100644 --- a/hyperion/model/tests/test_specific_energy_spectrum.py +++ b/hyperion/model/tests/test_specific_energy_spectrum.py @@ -251,9 +251,13 @@ def test_specific_energy_spectrum_yt_export_skips_frequency_resolved(tmpdir): def test_specific_energy_spectrum_mpi_matches_serial(tmpdir): # The saved spectrum must not depend on the number of MPI processes. We # compare the total (which is conserved, hence robust to Monte Carlo noise) - # between a serial run and a 2-process MPI run. + # between a serial run and a 2-process MPI run. This is skipped unless both an + # MPI launcher and the MPI build of the binary are available (e.g. CI only + # builds the serial binaries). if shutil.which('mpirun') is None and shutil.which('mpiexec') is None: pytest.skip("no MPI launcher available") + if shutil.which('hyperion_car_mpi') is None: + pytest.skip("MPI build of the Hyperion binaries not available") m = _cartesian_model('last', n_photons=200000) m.write(tmpdir.join(random_id()).strpath)