From 53862f3102c108c9682f38f85d115bce4b341506 Mon Sep 17 00:00:00 2001 From: dnarayanan Date: Mon, 4 Oct 2021 09:34:48 -0400 Subject: [PATCH 01/14] 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 5aab88c2d53982ce054a2971ec18716eabb2445c Mon Sep 17 00:00:00 2001 From: dnarayanan Date: Wed, 10 Nov 2021 10:27:49 -0500 Subject: [PATCH 02/14] 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 ce7680009f9df7064b85738745850a521eaace8a Mon Sep 17 00:00:00 2001 From: dnarayanan Date: Tue, 23 Nov 2021 09:36:38 -0500 Subject: [PATCH 03/14] 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 6d4f86dc..96d8f701 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.string_(self.output_density.encode('utf-8')) group.attrs['output_density_diff'] = np.string_(self.output_density_diff.encode('utf-8')) group.attrs['output_specific_energy'] = np.string_(self.output_specific_energy.encode('utf-8')) + group.attrs['output_specific_energy_nu'] = np.string_(self.output_specific_energ_nu.encode('utf-8')) group.attrs['output_n_photons'] = np.string_(self.output_n_photons.encode('utf-8')) diff --git a/hyperion/model/model.py b/hyperion/model/model.py index 2a0ce3d4..9fe7b53d 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 8c8145db1b551df06bfc9af3db3206c1c780d95e Mon Sep 17 00:00:00 2001 From: dnarayanan Date: Wed, 24 Nov 2021 15:54:52 -0500 Subject: [PATCH 04/14] 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 96d8f701..d6d1acb6 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.string_(self.output_density.encode('utf-8')) group.attrs['output_density_diff'] = np.string_(self.output_density_diff.encode('utf-8')) group.attrs['output_specific_energy'] = np.string_(self.output_specific_energy.encode('utf-8')) - group.attrs['output_specific_energy_nu'] = np.string_(self.output_specific_energ_nu.encode('utf-8')) + group.attrs['output_specific_energy_nu'] = np.string_(self.output_specific_energy_nu.encode('utf-8')) group.attrs['output_n_photons'] = np.string_(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 d38e048c54d331c0c45bac61da1d92fefd5b5694 Mon Sep 17 00:00:00 2001 From: dnarayanan Date: Thu, 2 Dec 2021 08:29:06 -0500 Subject: [PATCH 05/14] 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 033472e0403578696df58bfcec43fba782f2c583 Mon Sep 17 00:00:00 2001 From: dnarayanan Date: Thu, 9 Dec 2021 15:32:13 -0500 Subject: [PATCH 06/14] 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 74fd6e63b8beadabb759eb24a81009474be5eb0d Mon Sep 17 00:00:00 2001 From: dnarayanan Date: Mon, 13 Dec 2021 13:09:07 -0500 Subject: [PATCH 07/14] 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 b54d60891ed12d47299c20d8c9a8c59365a077ce Mon Sep 17 00:00:00 2001 From: dnarayanan Date: Tue, 14 Dec 2021 11:09:55 -0500 Subject: [PATCH 08/14] 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 d6d1acb6..66743cc9 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 7c588b1f41ae02cc666884ae68ed3f1719c3e038 Mon Sep 17 00:00:00 2001 From: dnarayanan Date: Tue, 14 Dec 2021 11:50:39 -0500 Subject: [PATCH 09/14] 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 66743cc9..16f97d7e 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 97b2dc683dc97760f1fdb30239b78dd727281a6d Mon Sep 17 00:00:00 2001 From: dnarayanan Date: Mon, 9 May 2022 08:51:16 -0400 Subject: [PATCH 10/14] 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 9fe7b53d..788662f5 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 6d0f51de562aefc4ddb9161dde6362b67e1ebbe7 Mon Sep 17 00:00:00 2001 From: dnarayanan Date: Mon, 9 May 2022 12:22:59 -0400 Subject: [PATCH 11/14] 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 788662f5..b5e1442c 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 2ffe74dd53030f80f22582de836cd47b75fad707 Mon Sep 17 00:00:00 2001 From: dnarayanan Date: Tue, 10 May 2022 12:14:36 -0400 Subject: [PATCH 12/14] 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 8717014c8fb61c588d3fdbb682b90ca3edb3df7f Mon Sep 17 00:00:00 2001 From: dnarayanan Date: Tue, 24 May 2022 12:42:56 -0400 Subject: [PATCH 13/14] 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 a0184e5b8f6e8ec6778d8cf1283e7398c7a1d2b1 Mon Sep 17 00:00:00 2001 From: dnarayanan Date: Thu, 11 Aug 2022 11:49:45 -0400 Subject: [PATCH 14/14] 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.