From fd551c5b010c0f8ac6135937883b71be04ad986a Mon Sep 17 00:00:00 2001 From: camillelandri Date: Wed, 1 Jul 2026 11:44:46 +0200 Subject: [PATCH 01/28] include hdf5 librairies when compiling with Krome --- build/Makefile | 1 + 1 file changed, 1 insertion(+) diff --git a/build/Makefile b/build/Makefile index c0c7b49e1..3c0bfcc76 100644 --- a/build/Makefile +++ b/build/Makefile @@ -266,6 +266,7 @@ ifeq ($(SINKTREE), yes) endif ifeq ($(KROME), yes) + HDF5 := yes FPPFLAGS += -DKROME LDFLAGS += -llapack endif From 213dce01faa7195c70c9d1140d66cc3ad8c6cf41 Mon Sep 17 00:00:00 2001 From: camillelandri Date: Thu, 2 Jul 2026 11:30:12 +0200 Subject: [PATCH 02/28] add Mats' new version of the raytracer (was never merged) --- src/main/ptmass_radiation.f90 | 12 ++-- src/main/utils_raytracer.f90 | 111 ++++++++++++++++++---------------- 2 files changed, 68 insertions(+), 55 deletions(-) diff --git a/src/main/ptmass_radiation.f90 b/src/main/ptmass_radiation.f90 index ef7cf2177..c31b09fcb 100644 --- a/src/main/ptmass_radiation.f90 +++ b/src/main/ptmass_radiation.f90 @@ -284,16 +284,20 @@ subroutine get_dust_temperature(npart,xyzh,eos_vars,nptmass,xyzmh_ptmass,dust_te if (iget_tdust == 4) then ! update tau_Lucy if (idust_opacity == 2) then - call get_all_tau(npart, nptmass, xyzmh_ptmass, xyzh, nucleation(:,ikappa), iray_resolution, tau_lucy) + call get_all_tau(npart, nptmass, xyzmh_ptmass, xyzh, nucleation(:,ikappa),& + iray_resolution, .false., tau_lucy) else - call get_all_tau(npart, nptmass, xyzmh_ptmass, xyzh, calc_kappa_bowen(dust_temp(1:npart)), iray_resolution, tau_lucy) + call get_all_tau(npart, nptmass, xyzmh_ptmass, xyzh, calc_kappa_bowen(dust_temp(1:npart)),& + iray_resolution, .false., tau_lucy) endif elseif (itau_alloc == 1) then ! update tau if (idust_opacity == 2) then - call get_all_tau(npart, nptmass, xyzmh_ptmass, xyzh, nucleation(:,ikappa), iray_resolution, tau) + call get_all_tau(npart, nptmass, xyzmh_ptmass, xyzh, nucleation(:,ikappa),& + iray_resolution, .true., tau) else - call get_all_tau(npart, nptmass, xyzmh_ptmass, xyzh, calc_kappa_bowen(dust_temp(1:npart)), iray_resolution, tau) + call get_all_tau(npart, nptmass, xyzmh_ptmass, xyzh, calc_kappa_bowen(dust_temp(1:npart)),& + iray_resolution, .true., tau) endif endif ! diff --git a/src/main/utils_raytracer.f90 b/src/main/utils_raytracer.f90 index 1a14eddf8..3a1595ced 100644 --- a/src/main/utils_raytracer.f90 +++ b/src/main/utils_raytracer.f90 @@ -1,6 +1,6 @@ !--------------------------------------------------------------------------! ! The Phantom Smoothed Particle Hydrodynamics code, by Daniel Price et al. ! -! Copyright (c) 2007-2026 The Authors (see AUTHORS) ! +! Copyright (c) 2007-2025 The Authors (see AUTHORS) ! ! See LICENCE file for usage and distribution conditions ! ! http://phantomsph.github.io/ ! !--------------------------------------------------------------------------! @@ -25,7 +25,7 @@ module raytracer use healpix implicit none - public :: get_all_tau + public :: get_all_tau, get_all_tau_single private @@ -43,24 +43,26 @@ module raytracer ! IN: xyzh: The array containing the particles position+smooting lenght ! IN: kappa_cgs: The array containing the opacities of all SPH particles ! IN: order: The healpix order which is used for the uniform ray sampling + ! IN: outwards: The direction of the ray, true for outwards false for inwards !+ ! OUT: tau: The array of optical depths for each SPH particle !+ !------------------------------------------------------------------------------------ -subroutine get_all_tau(npart, nptmass, xyzmh_ptmass, xyzh, kappa_cgs, order, tau) +subroutine get_all_tau(npart, nptmass, xyzmh_ptmass, xyzh, kappa_cgs, order, outwards, tau) use part, only: iReff - integer, intent(in) :: npart, order, nptmass - real, intent(in) :: kappa_cgs(:), xyzh(:,:), xyzmh_ptmass(:,:) - real, intent(out) :: tau(:) + integer, intent(in) :: npart, order, nptmass + real, intent(in) :: kappa_cgs(:), xyzh(:,:), xyzmh_ptmass(:,:) + logical, intent(in) :: outwards + real, intent(out) :: tau(:) real :: Rinject Rinject = xyzmh_ptmass(iReff,1) if (nptmass == 2 ) then call get_all_tau_companion(npart, xyzmh_ptmass(1:3,1), xyzmh_ptmass(iReff,1), xyzh, kappa_cgs, & - Rinject, xyzmh_ptmass(1:3,2), xyzmh_ptmass(iReff,2), order, tau) + Rinject, xyzmh_ptmass(1:3,2), xyzmh_ptmass(iReff,2), order, outwards, tau) else call get_all_tau_single(npart, xyzmh_ptmass(1:3,1), xyzmh_ptmass(iReff,1), xyzh,& - kappa_cgs, Rinject, order, tau) + kappa_cgs, Rinject, order, outwards, tau) endif end subroutine get_all_tau @@ -78,20 +80,22 @@ end subroutine get_all_tau ! IN: Rstar: The radius of the primary star ! IN: Rinject: The particles injection radius ! IN: order: The healpix order which is used for the uniform ray sampling + ! IN: outwards: The direction of the ray, true for outwards false for inwards !+ ! OUT: taus: The array of optical depths to each SPH particle !+ !--------------------------------------------------------------------------------- -subroutine get_all_tau_single(npart, primary, Rstar, xyzh, kappa, Rinject, order, tau) - use part, only:isdead_or_accreted - integer, intent(in) :: npart,order - real, intent(in) :: primary(3), kappa(:), Rstar, Rinject, xyzh(:,:) - real, intent(out) :: tau(:) +subroutine get_all_tau_single(npart, primary, Rstar, xyzh, kappa, Rinject, order, outwards, tau) + use part, only : isdead_or_accreted + integer, intent(in) :: npart,order + real, intent(in) :: primary(3), Rstar, xyzh(:,:), kappa(:), Rinject + logical, intent(in) :: outwards + real, intent(out) :: tau(:) integer :: i, nrays, nsides real :: ray_dir(3),part_dir(3) - real, allocatable :: rays_dist(:,:), rays_tau(:,:) - integer, allocatable :: rays_dim(:) + real, dimension(:,:), allocatable :: rays_dist, rays_tau + integer, dimension(:), allocatable :: rays_dim integer, parameter :: ndim = 200 ! maximum number of points along the ray where tau is calculated nrays = 12*4**order ! The number of rays traced given the healpix order @@ -108,17 +112,18 @@ subroutine get_all_tau_single(npart, primary, Rstar, xyzh, kappa, Rinject, order !$omp parallel default(none) & !$omp private(ray_dir) & - !$omp shared(nrays,nsides,primary,kappa,xyzh,Rstar,Rinject,rays_dist,rays_tau,rays_dim) + !$omp shared(nrays,nsides,primary,kappa,xyzh,Rstar,Rinject,rays_dist,rays_tau,rays_dim,outwards) !$omp do do i = 1, nrays !returns ray_dir, the unit vector identifying a ray (index i-1 because healpix starts counting from index 0) call pix2vec_nest(nsides, i-1, ray_dir) !calculate the properties along the ray (tau, distance, number of points) - call ray_tracer(primary,ray_dir,xyzh,kappa,Rstar,Rinject,rays_tau(:,i),rays_dist(:,i),rays_dim(i)) + call ray_tracer(primary,ray_dir,xyzh,kappa,Rstar,Rinject,outwards,rays_tau(:,i),rays_dist(:,i),rays_dim(i)) enddo !$omp enddo !$omp end parallel + !_---------------------------------------------- ! DETERMINE the optical depth for each particle ! using the values available on the HEALPix rays @@ -157,21 +162,23 @@ end subroutine get_all_tau_single ! IN: companion: The xyz coordinates of the companion ! IN: Rcomp: The radius of the companion ! IN: order: The healpix order which is used for the uniform ray sampling + ! IN: outwards: The direction of the ray, true for outwards false for inwards !+ ! OUT: tau: The array of optical depths for each SPH particle !+ !-------------------------------------------------------------------------- -subroutine get_all_tau_companion(npart, primary, Rstar, xyzh, kappa, Rinject, companion, Rcomp, order, tau) - use part, only:isdead_or_accreted - integer, intent(in) :: npart, order - real, intent(in) :: primary(3), companion(3), kappa(:), Rstar, Rinject, xyzh(:,:), Rcomp - real, intent(out) :: tau(:) +subroutine get_all_tau_companion(npart, primary, Rstar, xyzh, kappa, Rinject, companion, Rcomp, order, outwards, tau) + use part, only : isdead_or_accreted + integer, intent(in) :: npart, order + real, intent(in) :: primary(3), companion(3), kappa(:), Rstar, Rinject, xyzh(:,:), Rcomp + logical, intent(in) :: outwards + real, intent(out) :: tau(:) integer :: i, nrays, nsides real :: normCompanion,theta0,phi,cosphi,sinphi,theta,sep,root real :: ray_dir(3),part_dir(3),uvecCompanion(3) - real, allocatable :: rays_dist(:,:), rays_tau(:,:) - integer, allocatable :: rays_dim(:) + real, dimension(:,:), allocatable :: rays_dist, rays_tau + integer, dimension(:), allocatable :: rays_dim integer, parameter :: ndim = 200 ! maximum number of points along the ray where tau is calculated nrays = 12*4**order ! The number of rays traced given the healpix order @@ -197,7 +204,7 @@ subroutine get_all_tau_companion(npart, primary, Rstar, xyzh, kappa, Rinject, co !$omp parallel default(none) & !$omp private(ray_dir,theta,root,sep) & !$omp shared(nrays,nsides,primary,kappa,xyzh,Rstar,Rinject,Rcomp,rays_dist,rays_tau,rays_dim) & - !$omp shared(uvecCompanion,normCompanion,cosphi,sinphi,theta0) + !$omp shared(uvecCompanion,normCompanion,cosphi,sinphi,theta0,outwards) !$omp do do i = 1, nrays !returns ray_dir, the unit vector identifying a ray (index i-1 because healpix starts counting from index 0) @@ -210,9 +217,9 @@ subroutine get_all_tau_companion(npart, primary, Rstar, xyzh, kappa, Rinject, co if (theta < theta0) then root = sqrt(Rcomp**2-normCompanion**2*sin(theta)**2) sep = normCompanion*cos(theta)-root - call ray_tracer(primary,ray_dir,xyzh,kappa,Rstar,Rinject,rays_tau(:,i),rays_dist(:,i),rays_dim(i), sep) + call ray_tracer(primary,ray_dir,xyzh,kappa,Rstar,Rinject,outwards,rays_tau(:,i),rays_dist(:,i),rays_dim(i), sep) else - call ray_tracer(primary,ray_dir,xyzh,kappa,Rstar,Rinject,rays_tau(:,i),rays_dist(:,i),rays_dim(i)) + call ray_tracer(primary,ray_dir,xyzh,kappa,Rstar,Rinject,outwards,rays_tau(:,i),rays_dist(:,i),rays_dim(i)) endif enddo !$omp enddo @@ -262,9 +269,9 @@ end subroutine get_all_tau_companion !+ !-------------------------------------------------------------------------- subroutine interpolate_tau(nsides, vec, rays_tau, rays_dist, rays_dim, tau) - integer, intent(in) :: nsides, rays_dim(:) - real, intent(in) :: vec(:), rays_dist(:,:), rays_tau(:,:) - real, intent(out) :: tau + integer, intent(in) :: nsides, rays_dim(:) + real, intent(in) :: vec(:), rays_dist(:,:), rays_tau(:,:) + real, intent(out) :: tau integer :: rayIndex, neighbours(8), nneigh, i, k real :: tautemp, ray(3), vectemp(3), weight, tempdist(8), distRay_sq, vec_norm2 @@ -312,6 +319,7 @@ subroutine interpolate_tau(nsides, vec, rays_tau, rays_dist, rays_dim, tau) tau = tau / weight end subroutine interpolate_tau + !-------------------------------------------------------------------------- !+ ! Interpolation of the optical depth for an arbitrary point on the ray, @@ -327,9 +335,9 @@ end subroutine interpolate_tau !+ !-------------------------------------------------------------------------- subroutine get_tau_on_ray(distance, tau_along_ray, dist_along_ray, len, tau) - real, intent(in) :: distance, tau_along_ray(:), dist_along_ray(:) - integer, intent(in) :: len - real, intent(out) :: tau + real, intent(in) :: distance, tau_along_ray(:), dist_along_ray(:) + integer, intent(in) :: len + real, intent(out) :: tau integer :: L, R, m ! left, right and middle index for binary search @@ -366,6 +374,7 @@ end subroutine get_tau_on_ray ! IN: kappa: The array containing the particles opacity ! IN: Rstar: The radius of the primary star ! IN: Rinject: The particles injection radius + ! IN: outwards: The direction of the ray, true for outwards false for inwards !+ ! OUT: tau_along_ray: The vector of cumulative optical depth along the ray ! OUT: dist_along_ray: The vector of distances from the primary along the ray @@ -374,13 +383,14 @@ end subroutine get_tau_on_ray ! OPT: maxDistance: The maximal distance the ray needs to be traced !+ !-------------------------------------------------------------------------- -subroutine ray_tracer(primary, ray, xyzh, kappa, Rstar, Rinject, tau_along_ray, dist_along_ray, len, maxDistance) +subroutine ray_tracer(primary, ray, xyzh, kappa, Rstar, Rinject, outwards, tau_along_ray, dist_along_ray, len, maxDistance) use units, only:unit_opacity use part, only:itauL_alloc - real, intent(in) :: primary(3), ray(3), Rstar, Rinject, xyzh(:,:), kappa(:) - real, intent(out) :: dist_along_ray(:), tau_along_ray(:) - integer, intent(out) :: len + real, intent(in) :: primary(3), ray(3), Rstar, Rinject, xyzh(:,:), kappa(:) + logical, optional :: outwards real, optional :: maxDistance + real, intent(out) :: dist_along_ray(:), tau_along_ray(:) + integer, intent(out) :: len real, parameter :: tau_max = 99. real :: dr, next_dr, h, dtaudr, previousdtaudr, nextdtaudr, distance @@ -418,9 +428,11 @@ subroutine ray_tracer(primary, ray, xyzh, kappa, Rstar, Rinject, tau_along_ray, endif len = i - if (itauL_alloc > 0) then - !reverse integration start from zero inward + if (.not. outwards) then + !reverse integration start from zero outward and integrate inwards tau_along_ray(1:len) = tau_along_ray(len) - tau_along_ray(1:len) + endif + if (itauL_alloc > 0) then !find the first point where tau_lucy < 2/3 if (tau_along_ray(1) > 2./3.) then L = 1 @@ -442,7 +454,7 @@ end subroutine ray_tracer logical function hasNext(inext, tau, distance, maxDistance) integer, intent(in) :: inext - real, intent(in) :: tau, distance + real, intent(in) :: tau, distance real, optional :: maxDistance real :: tau_max = 99. if (present(maxDistance)) then @@ -474,9 +486,8 @@ end function hasNext !-------------------------------------------------------------------------- subroutine find_next(inpoint, h, ray, xyzh, kappa, dtaudr, distance, inext) use neighkdtree, only:getneigh_pos,leaf_is_active,listneigh - use kernel, only:radkern,cnormk,wkern - use part, only:hfact,rhoh,massoftype,igas - use dim, only:maxpsph + use kernel, only:radkern,cnormk,wkern + use part, only:hfact,rhoh,massoftype,igas real, intent(in) :: xyzh(:,:), kappa(:), inpoint(:), ray(:), h integer, intent(inout) :: inext real, intent(out) :: distance, dtaudr @@ -484,7 +495,7 @@ subroutine find_next(inpoint, h, ray, xyzh, kappa, dtaudr, distance, inext) integer, parameter :: nmaxcache = 0 real :: xyzcache(0,nmaxcache) - integer :: nneigh, i, prev,j + integer :: nneigh, i, prev real :: dmin, vec(3), dr, raydistance, q, norm_sq prev = inext @@ -498,23 +509,21 @@ subroutine find_next(inpoint, h, ray, xyzh, kappa, dtaudr, distance, inext) dmin = huge(0.) !loop over all neighbours do i=1,nneigh - j = listneigh(i) - if (j > maxpsph) cycle - vec = xyzh(1:3,j) - inpoint + vec = xyzh(1:3,listneigh(i)) - inpoint norm_sq = dot_product(vec,vec) - q = sqrt(norm_sq)/xyzh(4,j) + q = sqrt(norm_sq)/xyzh(4,listneigh(i)) !add optical depth contribution from each particle - dtaudr = dtaudr+wkern(q*q,q)*kappa(j)*rhoh(xyzh(4,j), massoftype(igas)) + dtaudr = dtaudr+wkern(q*q,q)*kappa(listneigh(i))*rhoh(xyzh(4,listneigh(i)), massoftype(igas)) ! find the next particle : among the neighbours find the particle located the closest to the ray - if (j /= prev) then + if (listneigh(i) /= prev) then dr = dot_product(vec,ray) !projected distance along the ray if (dr>0.) then !distance perpendicular to the ray direction raydistance = norm_sq - dr**2 if (raydistance < dmin) then dmin = raydistance - inext = j + inext = listneigh(i) distance = dr endif endif From 5ea5fe6b6689b54bf75ef05c71a92c08990a3b83 Mon Sep 17 00:00:00 2001 From: camillelandri Date: Thu, 2 Jul 2026 11:46:03 +0200 Subject: [PATCH 03/28] - write hdf5 output (just one group with abundances for now) - adapt number of threads to 1 if Krome is compiled with gfortran --- build/Makefile | 2 +- src/utils/analysis_krome.f90 | 185 +++++++++++++++++++++++++---------- 2 files changed, 134 insertions(+), 53 deletions(-) diff --git a/build/Makefile b/build/Makefile index 3c0bfcc76..fc9c2e907 100644 --- a/build/Makefile +++ b/build/Makefile @@ -885,7 +885,7 @@ allanalysis: analysis_structurefn.f90) \ analysisfile="utils_timing.f90 io_structurefn.f90 utils_sort.f90 random.f90 struct_part.f90 $$analysisfile"; \ ;; \ - analysis_krome.f90|analysis_mcfost.f90) \ + analysis_krome.F90|analysis_mcfost.f90) \ continue; \ ;; \ esac; \ diff --git a/src/utils/analysis_krome.f90 b/src/utils/analysis_krome.f90 index 4d69e7f0f..89fd85fa6 100644 --- a/src/utils/analysis_krome.f90 +++ b/src/utils/analysis_krome.f90 @@ -20,12 +20,15 @@ module analysis use krome_user, only:krome_nmols use part, only: maxp use raytracer, only: get_all_tau + use hdf5 + use omp_lib + implicit none character(len=20), parameter, public :: analysistype = 'krome' public :: do_analysis logical, allocatable :: mask(:) - real, allocatable :: abundance(:,:), abundance_prev(:,:), one(:) + real(8), allocatable :: abundance(:,:), abundance_prev(:,:), one(:) character(len=16) :: abundance_label(krome_nmols) integer(8), allocatable :: iorig_old(:) integer, allocatable :: iprev(:) @@ -35,11 +38,12 @@ module analysis contains subroutine do_analysis(dumpfile,num,xyzh,vxyzu,particlemass,npart,time,iunit) - use part, only: isdead_or_accreted, iorig, rhoh, nptmass, xyzmh_ptmass, iReff - use neighkdtree, only:build_tree - use units, only: utime,unit_density + use part, only: isdead_or_accreted, iorig, rhoh, nptmass, xyzmh_ptmass, iReff, iboundary, igas, iphase, iamtype + use neighkdtree, only: build_tree + use units, only: utime,unit_density,udist + use physcon, only: atomic_mass_unit use eos, only: get_temperature, ieos, gamma,gmw, init_eos - use io, only: fatal + use io, only: fatal, iverbose use krome_main, only: krome_init, krome use krome_user, only: krome_get_names,krome_set_user_Auv,krome_set_user_xi,& krome_set_user_alb,krome_set_user_AuvAv @@ -49,17 +53,25 @@ subroutine do_analysis(dumpfile,num,xyzh,vxyzu,particlemass,npart,time,iunit) real, intent(in) :: particlemass,time real, save :: tprev = 0. integer, save :: nprev = 0 - real :: dt_cgs, rho_cgs, numberdensity, T_gas, gammai, mui, AUV, xi + real :: dt_cgs, rho_cgs, rholist(maxp), Tlist(maxp), mulist(maxp), Auvlist(maxp), xilist(maxp) + real :: numberdensity, T_gas, gammai, mui, AUV, xi real :: abundance_part(krome_nmols), Y(krome_nmols), column_density(npart), xyzh_copy(4,npart) - integer :: i, j, ierr, completed_iterations, npart_copy = 0 + real :: max_radius, radius + integer :: i, j, ierr, completed_iterations, npart_copy = 0, hdferr, i_radius = 1 + +#ifdef __GFORTRAN__ + print*, "Setting number of threads to 1 (KROME is not thread-safe when compiled with gfortran)" + call omp_set_num_threads(1) +#endif if (.not.done_init) then done_init = .true. + print*, "initialising KROME" call krome_init() print*, "Initialised KROME" + abundance_label(:) = krome_get_names() - maxp = maxp /2 allocate(abundance(krome_nmols,maxp)) abundance = 0. allocate(abundance_prev(krome_nmols,maxp)) @@ -70,6 +82,7 @@ subroutine do_analysis(dumpfile,num,xyzh,vxyzu,particlemass,npart,time,iunit) iorig_old = 0 allocate(iprev(maxp)) iprev = 0 + print*, "setting abundances" !$omp parallel do default(none) & !$omp shared(npart,xyzh,vxyzu,dt_cgs,nprev,iorig,iorig_old,iprev) & @@ -84,58 +97,100 @@ subroutine do_analysis(dumpfile,num,xyzh,vxyzu,particlemass,npart,time,iunit) enddo call init_eos(ieos, ierr) if (ierr /= 0) call fatal(analysistype, "Failed to initialise EOS") - else + + ! initialize HDF5 + call h5open_f(hdferr) + if (hdferr /= 0) then + print*,'ERROR: Failed to initialize HDF5' + return + endif + + else dt_cgs = (time - tprev)*utime completed_iterations = 0 - print*, "not first step data, timestep = ",dt_cgs, "npart = ",npart, "nprev = ",nprev + print*, dumpfile, ": not first step data, timestep = ",dt_cgs, "npart = ",npart, "nprev = ",nprev + xyzmh_ptmass(iReff,1) = 2. npart_copy = npart xyzh_copy = xyzh(:,:npart) call build_tree(npart_copy,npart_copy,xyzh_copy,vxyzu) call get_all_tau(npart, nptmass, xyzmh_ptmass, xyzh, one, 5, .false., column_density) + + max_radius = 0.0 + do i = 1, npart + if (.not.isdead_or_accreted(xyzh(4, i))) then + radius = sqrt(xyzh(1, i)**2 + xyzh(2, i)**2 + xyzh(3, i)**2) + if (radius > max_radius) then + max_radius = radius + i_radius = i + endif + endif + enddo + column_density = column_density + rhoh(xyzh(4,i_radius),particlemass)*unit_density * max_radius * udist + + rholist = 0. + Tlist = 0. + mulist = 0. + Auvlist = 0. + xilist = 0. + !$omp parallel do default(none) & - !$omp shared(npart,xyzh,vxyzu,dt_cgs,nprev,iorig,iorig_old,iprev) & - !$omp shared(abundance,abundance_prev,particlemass,unit_density) & + !$omp shared(npart,xyzh,vxyzu,dt_cgs,nprev,iorig,iorig_old,iprev,iverbose) & + !$omp shared(abundance,abundance_prev,particlemass,unit_density,udist) & !$omp shared(ieos,gamma,gmw,time,completed_iterations,column_density,AuvAv,albedo) & + !$omp shared(rholist,Tlist,mulist,Auvlist,xilist,iphase) & !$omp private(i,j,abundance_part,Y,rho_cgs,numberdensity,T_gas,gammai,mui,AUV,xi) + outer: do i=1,npart if (.not.isdead_or_accreted(xyzh(4,i))) then + inner: do j=1,nprev if (iorig(i) == iorig_old(j)) then iprev(i) = j exit inner endif enddo inner + if (j == iprev(i)) then abundance_part(:) = abundance_prev(:,iprev(i)) else call chem_init(abundance_part) endif - !Thermodynamic quantities - rho_cgs = rhoh(xyzh(4,i),particlemass)*unit_density - gammai = gamma - mui = gmw - numberdensity = rho_cgs / (mui * 1.6605E-24) - T_gas = get_temperature(ieos,xyzh(1:3, i),rhoh(xyzh(4,i),particlemass),vxyzu(:,i),gammai,mui) - T_gas = max(T_gas,20.0d0) - - !Radiation quantities - AUV = 4.65 * column_density(i) / 1.87e21 - xi = get_xi(AUV) - - call krome_set_user_Auv(AUV) - call krome_set_user_xi(xi) - call krome_set_user_alb(ALBEDO) - call krome_set_user_AuvAv(AuvAv) - - Y = abundance_part*numberdensity - call krome(Y,T_gas,dt_cgs) - abundance_part = Y/numberdensity - abundance(:,i) = abundance_part + if (iamtype(iphase(i)) /= iboundary .and. i > 2460) then ! 2460 is the amount of boundary particles + !Thermodynamic quantities + rho_cgs = rhoh(xyzh(4,i),particlemass)*unit_density + gammai = gamma + mui = gmw + numberdensity = rho_cgs / (mui * atomic_mass_unit) + T_gas = get_temperature(ieos,xyzh(1:3, i),rhoh(xyzh(4,i),particlemass),vxyzu(:,i),gammai,mui) + T_gas = max(T_gas,20.0d0) + + !Radiation quantities + AUV = AuvAv * column_density(i) / (mui * atomic_mass_unit) / 1.87e21 + xi = get_xi(AUV) + + call krome_set_user_Auv(AUV) + call krome_set_user_xi(xi) + call krome_set_user_alb(ALBEDO) + call krome_set_user_AuvAv(AuvAv) + + Y = abundance_part*numberdensity + call krome(Y,T_gas,dt_cgs) + abundance_part = Y/numberdensity + abundance(:,i) = abundance_part + + rholist(i) = rho_cgs + Tlist(i) = T_gas + mulist(i) = mui + Auvlist(i) = AUV + xilist(i) = xi + endif + endif + if (iverbose > 1) then + !$omp atomic + completed_iterations = completed_iterations + 1 + print*, 'Completed ', completed_iterations, ' of ', npart endif - !$omp atomic - completed_iterations = completed_iterations + 1 - print*, 'Completed ', completed_iterations, ' of ', npart enddo outer endif @@ -178,25 +233,51 @@ real function get_xi(AUV) end function get_xi subroutine write_chem(npart, dumpfile) - use krome_user, only:krome_idx_He,krome_idx_C,krome_idx_N,krome_idx_O,& - krome_idx_H,krome_idx_S,krome_idx_Fe,krome_idx_Si,krome_idx_Mg,& - krome_idx_Na,krome_idx_P,krome_idx_F,krome_idx_CO,krome_idx_C2H2,& - krome_idx_C2H,krome_idx_H2,krome_idx_SiNC,krome_idx_e + use krome_user, only:krome_idx_He integer, intent(in) :: npart character(len=*), intent(in) :: dumpfile - integer :: i, iu - - open(newunit=iu,file=dumpfile//'.comp',status='replace',action='write') - write(iu, *) '# H, He, C, N, O, S, Fe, Si, Mg, Na, P, F, CO, C2H2, C2H, H2, SiNC, e-' - do i=1, npart - write(iu, *) abundance(krome_idx_H, i), abundance(krome_idx_He, i), abundance(krome_idx_C, i), & - abundance(krome_idx_N, i), abundance(krome_idx_O, i), abundance(krome_idx_S, i), & - abundance(krome_idx_Fe, i), abundance(krome_idx_Si, i), abundance(krome_idx_Mg, i), & - abundance(krome_idx_Na, i), abundance(krome_idx_P, i), abundance(krome_idx_F, i), & - abundance(krome_idx_CO, i), abundance(krome_idx_C2H2, i), abundance(krome_idx_C2H, i), & - abundance(krome_idx_H2, i), abundance(krome_idx_SiNC, i), abundance(krome_idx_e, i) + integer :: hdferr + integer :: i + integer(hid_t) :: file_id, dspace_id, dset_id, group_id + integer(hsize_t) :: dims(1) + + ! create HDF5 file + print "(1x,a)",'Writing to '//trim(dumpfile)//'.h5' + call h5fcreate_f(trim(dumpfile)//'.h5', H5F_ACC_TRUNC_F, file_id, hdferr) + if (hdferr /= 0) then + print*,'ERROR: Failed to create HDF5 file' + return + endif + + ! create group for particle data + call h5gcreate_f(file_id, 'chemistry', group_id, hdferr) + if (hdferr /= 0) then + print*,'ERROR: Failed to create HDF5 group' + return + endif + + ! create dataspace for particle datasets + dims(1) = npart + call h5screate_simple_f(1, dims, dspace_id, hdferr) + if (hdferr /= 0) then + print*,'ERROR: Failed to create HDF5 dataspace' + return + endif + do i=1,krome_nmols + ! create dataset for each abundance + call h5dcreate_f(group_id, trim(adjustl(abundance_label(i))), H5T_NATIVE_DOUBLE, dspace_id, dset_id, hdferr) + ! write abundance data to dataset + call h5dwrite_f(dset_id, H5T_NATIVE_DOUBLE, abundance(i,1:npart), dims, hdferr) + if (hdferr /= 0) then + print*,'ERROR: Failed to create/write dataset for '//trim(adjustl(abundance_label(i))) + return + endif + call h5dclose_f(dset_id, hdferr) enddo - close(iu) + call h5sclose_f(dspace_id, hdferr) + call h5gclose_f(group_id, hdferr) + call h5fclose_f(file_id, hdferr) + end subroutine write_chem From c6c562830d97360047a7b3a692f74d4cce1b5559 Mon Sep 17 00:00:00 2001 From: Camille Landri Date: Thu, 2 Jul 2026 11:47:33 +0200 Subject: [PATCH 04/28] Rename analysis_krome.f90 to analysis_krome.F90 allow for preprocessing of the file --- src/utils/{analysis_krome.f90 => analysis_krome.F90} | 0 1 file changed, 0 insertions(+), 0 deletions(-) rename src/utils/{analysis_krome.f90 => analysis_krome.F90} (100%) diff --git a/src/utils/analysis_krome.f90 b/src/utils/analysis_krome.F90 similarity index 100% rename from src/utils/analysis_krome.f90 rename to src/utils/analysis_krome.F90 From b4da6a9b99feebf2aefdb7aaac1e09c1e1cff29a Mon Sep 17 00:00:00 2001 From: camillelandri Date: Fri, 3 Jul 2026 10:18:16 +0200 Subject: [PATCH 05/28] read chemistry file of first dump if it exists --- src/utils/analysis_krome.F90 | 102 +++++++++++++++++++++++++++-------- 1 file changed, 79 insertions(+), 23 deletions(-) diff --git a/src/utils/analysis_krome.F90 b/src/utils/analysis_krome.F90 index 89fd85fa6..e10c3f86c 100644 --- a/src/utils/analysis_krome.F90 +++ b/src/utils/analysis_krome.F90 @@ -57,7 +57,7 @@ subroutine do_analysis(dumpfile,num,xyzh,vxyzu,particlemass,npart,time,iunit) real :: numberdensity, T_gas, gammai, mui, AUV, xi real :: abundance_part(krome_nmols), Y(krome_nmols), column_density(npart), xyzh_copy(4,npart) real :: max_radius, radius - integer :: i, j, ierr, completed_iterations, npart_copy = 0, hdferr, i_radius = 1 + integer :: i, j, isize=0, ierr, completed_iterations, npart_copy = 0, hdferr, i_radius = 1 #ifdef __GFORTRAN__ print*, "Setting number of threads to 1 (KROME is not thread-safe when compiled with gfortran)" @@ -67,6 +67,13 @@ subroutine do_analysis(dumpfile,num,xyzh,vxyzu,particlemass,npart,time,iunit) if (.not.done_init) then done_init = .true. + ! initialize HDF5 + call h5open_f(hdferr) + if (hdferr /= 0) then + print*,'ERROR: Failed to initialize HDF5' + return + endif + print*, "initialising KROME" call krome_init() print*, "Initialised KROME" @@ -83,28 +90,29 @@ subroutine do_analysis(dumpfile,num,xyzh,vxyzu,particlemass,npart,time,iunit) allocate(iprev(maxp)) iprev = 0 - print*, "setting abundances" - !$omp parallel do default(none) & - !$omp shared(npart,xyzh,vxyzu,dt_cgs,nprev,iorig,iorig_old,iprev) & - !$omp shared(abundance,abundance_prev,particlemass,unit_density) & - !$omp shared(ieos,rho_cgs,T_gas,j) & - !$omp private(i,abundance_part) - do i=1, npart - if (.not.isdead_or_accreted(xyzh(4,i))) then - call chem_init(abundance_part) - abundance(:,i) = abundance_part - endif - enddo + !check if a file with abundances already exists, if so read it in and use it as initial abundances + inquire(file=trim(dumpfile)//'.h5', size=isize) + if (isize > 0) then + print*, "Found existing abundance file, reading in abundances" + call read_chem(npart, dumpfile) + else + print*, "setting abundances" + !$omp parallel do default(none) & + !$omp shared(npart,xyzh,vxyzu,dt_cgs,nprev,iorig,iorig_old,iprev) & + !$omp shared(abundance,abundance_prev,particlemass,unit_density) & + !$omp shared(ieos,rho_cgs,T_gas,j) & + !$omp private(i,abundance_part) + do i=1, npart + if (.not.isdead_or_accreted(xyzh(4,i))) then + call chem_init(abundance_part) + abundance(:,i) = abundance_part + endif + enddo + call write_chem(npart, dumpfile) + endif call init_eos(ieos, ierr) if (ierr /= 0) call fatal(analysistype, "Failed to initialise EOS") - ! initialize HDF5 - call h5open_f(hdferr) - if (hdferr /= 0) then - print*,'ERROR: Failed to initialize HDF5' - return - endif - else dt_cgs = (time - tprev)*utime completed_iterations = 0 @@ -192,9 +200,10 @@ subroutine do_analysis(dumpfile,num,xyzh,vxyzu,particlemass,npart,time,iunit) print*, 'Completed ', completed_iterations, ' of ', npart endif enddo outer + call write_chem(npart, dumpfile) endif - call write_chem(npart, dumpfile) + nprev = npart tprev = time iorig_old(1:npart) = iorig(1:npart) @@ -277,10 +286,57 @@ subroutine write_chem(npart, dumpfile) call h5sclose_f(dspace_id, hdferr) call h5gclose_f(group_id, hdferr) call h5fclose_f(file_id, hdferr) - - end subroutine write_chem +subroutine read_chem(npart, dumpfile) + use krome_user, only:krome_idx_He + integer, intent(in) :: npart + character(len=*), intent(in) :: dumpfile + integer :: hdferr + integer :: i + integer(hid_t) :: file_id, dset_id, group_id, filespace_id + integer(hsize_t) :: max_dims(1), file_dims(1) + + ! open HDF5 file + print "(1x,a)",'Reading from '//trim(dumpfile)//'.h5' + call h5fopen_f(trim(dumpfile)//'.h5', H5F_ACC_RDONLY_F, file_id, hdferr) + if (hdferr /= 0) then + print*,'ERROR: Failed to open HDF5 file' + return + endif + + ! open group for particle data + call h5gopen_f(file_id, 'chemistry', group_id, hdferr) + if (hdferr /= 0) then + print*,'ERROR: Failed to open HDF5 group' + return + endif + + + + do i=1,krome_nmols + ! open dataset for each abundance + call h5dopen_f(group_id, trim(adjustl(abundance_label(i))), dset_id, hdferr) + ! check number of particles in file + call h5dget_space_f(dset_id, filespace_id, hdferr) + call h5sget_simple_extent_dims_f(filespace_id, file_dims, max_dims, hdferr) + if (file_dims(1) /= npart) then + print*,'ERROR: Number of particles in HDF5 file does not match current dump' + return + endif + ! read abundance data from dataset + call h5dread_f(dset_id, H5T_NATIVE_DOUBLE, abundance(i,1:npart), file_dims, hdferr) + if (hdferr /= 0) then + print*,'ERROR: Failed to open/read dataset for '//trim(adjustl(abundance_label(i))) + return + endif + call h5dclose_f(dset_id, hdferr) + call h5sclose_f(filespace_id, hdferr) + enddo + call h5gclose_f(group_id, hdferr) + call h5fclose_f(file_id, hdferr) +end subroutine read_chem + subroutine chem_init(abundance_part) use krome_user, only:krome_idx_H2,krome_idx_He,krome_idx_CO,krome_idx_C2H2,& krome_idx_HCN,krome_idx_N2,krome_idx_SiC2,krome_idx_CS,& From 55eb05fedaecf60f80e1c3c45a840d0b1418e7cd Mon Sep 17 00:00:00 2001 From: camillelandri Date: Thu, 2 Jul 2026 15:38:57 +0200 Subject: [PATCH 06/28] fix compilation with self compiled hdf5 libraries --- build/Makefile | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/build/Makefile b/build/Makefile index fc9c2e907..335123eec 100644 --- a/build/Makefile +++ b/build/Makefile @@ -405,7 +405,7 @@ ifeq (X$(HDF5LIB), X) endif FFLAGS+= -I$(HDF5INCLUDE) CCFLAGS+= -I$(HDF5INCLUDE) - LDFLAGS+= -L$(HDF5_DIR)/lib -lhdf5 -lhdf5_fortran + LDFLAGS+= -L$(HDF5LIB) -Wl,-rpath,$(HDF5LIB) -lhdf5 -lhdf5_fortran endif IDFLAGS=$(FPPFLAGS) From 462d2ad52df8e033963ed58906fdbef42aba0068 Mon Sep 17 00:00:00 2001 From: camillelandri Date: Fri, 3 Jul 2026 10:38:00 +0200 Subject: [PATCH 07/28] add timings --- src/utils/analysis_krome.F90 | 28 +++++++++++++++++----------- 1 file changed, 17 insertions(+), 11 deletions(-) diff --git a/src/utils/analysis_krome.F90 b/src/utils/analysis_krome.F90 index e10c3f86c..54ae45061 100644 --- a/src/utils/analysis_krome.F90 +++ b/src/utils/analysis_krome.F90 @@ -56,12 +56,14 @@ subroutine do_analysis(dumpfile,num,xyzh,vxyzu,particlemass,npart,time,iunit) real :: dt_cgs, rho_cgs, rholist(maxp), Tlist(maxp), mulist(maxp), Auvlist(maxp), xilist(maxp) real :: numberdensity, T_gas, gammai, mui, AUV, xi real :: abundance_part(krome_nmols), Y(krome_nmols), column_density(npart), xyzh_copy(4,npart) - real :: max_radius, radius + real :: max_radius, radius, t_start, t_stop integer :: i, j, isize=0, ierr, completed_iterations, npart_copy = 0, hdferr, i_radius = 1 #ifdef __GFORTRAN__ print*, "Setting number of threads to 1 (KROME is not thread-safe when compiled with gfortran)" call omp_set_num_threads(1) +#else + print*, "running with ", omp_get_max_threads(), " threads" #endif if (.not.done_init) then @@ -116,8 +118,9 @@ subroutine do_analysis(dumpfile,num,xyzh,vxyzu,particlemass,npart,time,iunit) else dt_cgs = (time - tprev)*utime completed_iterations = 0 - print*, dumpfile, ": not first step data, timestep = ",dt_cgs, "npart = ",npart, "nprev = ",nprev - + print*, "not first step data, timestep = ",dt_cgs, "npart = ",npart, "nprev = ",nprev + print*, "Calculating column density" + call cpu_time(t_start) xyzmh_ptmass(iReff,1) = 2. npart_copy = npart xyzh_copy = xyzh(:,:npart) @@ -135,13 +138,11 @@ subroutine do_analysis(dumpfile,num,xyzh,vxyzu,particlemass,npart,time,iunit) endif enddo column_density = column_density + rhoh(xyzh(4,i_radius),particlemass)*unit_density * max_radius * udist + call cpu_time(t_stop) + print*, " - Took ", t_stop - t_start, " seconds" - rholist = 0. - Tlist = 0. - mulist = 0. - Auvlist = 0. - xilist = 0. - + print*, "Running KROME" + call cpu_time(t_start) !$omp parallel do default(none) & !$omp shared(npart,xyzh,vxyzu,dt_cgs,nprev,iorig,iorig_old,iprev,iverbose) & !$omp shared(abundance,abundance_prev,particlemass,unit_density,udist) & @@ -151,7 +152,6 @@ subroutine do_analysis(dumpfile,num,xyzh,vxyzu,particlemass,npart,time,iunit) outer: do i=1,npart if (.not.isdead_or_accreted(xyzh(4,i))) then - inner: do j=1,nprev if (iorig(i) == iorig_old(j)) then iprev(i) = j @@ -200,10 +200,16 @@ subroutine do_analysis(dumpfile,num,xyzh,vxyzu,particlemass,npart,time,iunit) print*, 'Completed ', completed_iterations, ' of ', npart endif enddo outer + call cpu_time(t_stop) + print*, " - Took ", (t_stop - t_start)/omp_get_max_threads(), " seconds" + + call cpu_time(t_start) call write_chem(npart, dumpfile) + call cpu_time(t_stop) + print*, " - Took ", t_stop - t_start, " seconds" endif - + ! keep track of previous abundances for next dump nprev = npart tprev = time iorig_old(1:npart) = iorig(1:npart) From faf63a76582bd6c4913ecbc0fb6f6759df92ae78 Mon Sep 17 00:00:00 2001 From: camillelandri Date: Fri, 3 Jul 2026 12:04:45 +0200 Subject: [PATCH 08/28] change hdf5 layout - one dataset with chem matrix instread of nmols datasets with particle array --- src/utils/analysis_krome.F90 | 198 +++++++++++++++++++++++++---------- 1 file changed, 140 insertions(+), 58 deletions(-) diff --git a/src/utils/analysis_krome.F90 b/src/utils/analysis_krome.F90 index 54ae45061..0c5210192 100644 --- a/src/utils/analysis_krome.F90 +++ b/src/utils/analysis_krome.F90 @@ -53,10 +53,10 @@ subroutine do_analysis(dumpfile,num,xyzh,vxyzu,particlemass,npart,time,iunit) real, intent(in) :: particlemass,time real, save :: tprev = 0. integer, save :: nprev = 0 - real :: dt_cgs, rho_cgs, rholist(maxp), Tlist(maxp), mulist(maxp), Auvlist(maxp), xilist(maxp) + real :: dt_cgs, rho_cgs real :: numberdensity, T_gas, gammai, mui, AUV, xi real :: abundance_part(krome_nmols), Y(krome_nmols), column_density(npart), xyzh_copy(4,npart) - real :: max_radius, radius, t_start, t_stop + real :: max_radius, radius, tstart integer :: i, j, isize=0, ierr, completed_iterations, npart_copy = 0, hdferr, i_radius = 1 #ifdef __GFORTRAN__ @@ -119,14 +119,17 @@ subroutine do_analysis(dumpfile,num,xyzh,vxyzu,particlemass,npart,time,iunit) dt_cgs = (time - tprev)*utime completed_iterations = 0 print*, "not first step data, timestep = ",dt_cgs, "npart = ",npart, "nprev = ",nprev - print*, "Calculating column density" - call cpu_time(t_start) + print*, "Building neighbour tree..." + tstart = omp_get_wtime() xyzmh_ptmass(iReff,1) = 2. npart_copy = npart xyzh_copy = xyzh(:,:npart) call build_tree(npart_copy,npart_copy,xyzh_copy,vxyzu) + print*, " - Took ", omp_get_wtime() - tstart, " seconds" + + print*, "Calculating column density..." + tstart = omp_get_wtime() call get_all_tau(npart, nptmass, xyzmh_ptmass, xyzh, one, 5, .false., column_density) - max_radius = 0.0 do i = 1, npart if (.not.isdead_or_accreted(xyzh(4, i))) then @@ -138,16 +141,14 @@ subroutine do_analysis(dumpfile,num,xyzh,vxyzu,particlemass,npart,time,iunit) endif enddo column_density = column_density + rhoh(xyzh(4,i_radius),particlemass)*unit_density * max_radius * udist - call cpu_time(t_stop) - print*, " - Took ", t_stop - t_start, " seconds" + print*, " - Took ", omp_get_wtime() - tstart, " seconds" print*, "Running KROME" - call cpu_time(t_start) + tstart = omp_get_wtime() !$omp parallel do default(none) & !$omp shared(npart,xyzh,vxyzu,dt_cgs,nprev,iorig,iorig_old,iprev,iverbose) & - !$omp shared(abundance,abundance_prev,particlemass,unit_density,udist) & + !$omp shared(abundance,abundance_prev,particlemass,unit_density,udist,iphase) & !$omp shared(ieos,gamma,gmw,time,completed_iterations,column_density,AuvAv,albedo) & - !$omp shared(rholist,Tlist,mulist,Auvlist,xilist,iphase) & !$omp private(i,j,abundance_part,Y,rho_cgs,numberdensity,T_gas,gammai,mui,AUV,xi) outer: do i=1,npart @@ -186,12 +187,6 @@ subroutine do_analysis(dumpfile,num,xyzh,vxyzu,particlemass,npart,time,iunit) call krome(Y,T_gas,dt_cgs) abundance_part = Y/numberdensity abundance(:,i) = abundance_part - - rholist(i) = rho_cgs - Tlist(i) = T_gas - mulist(i) = mui - Auvlist(i) = AUV - xilist(i) = xi endif endif if (iverbose > 1) then @@ -200,13 +195,11 @@ subroutine do_analysis(dumpfile,num,xyzh,vxyzu,particlemass,npart,time,iunit) print*, 'Completed ', completed_iterations, ' of ', npart endif enddo outer - call cpu_time(t_stop) - print*, " - Took ", (t_stop - t_start)/omp_get_max_threads(), " seconds" - - call cpu_time(t_start) + print*, " - Took ", omp_get_wtime() - tstart, " seconds" + + tstart = omp_get_wtime() call write_chem(npart, dumpfile) - call cpu_time(t_stop) - print*, " - Took ", t_stop - t_start, " seconds" + print*, " - Took ", omp_get_wtime() - tstart, " seconds" endif ! keep track of previous abundances for next dump @@ -248,13 +241,13 @@ real function get_xi(AUV) end function get_xi subroutine write_chem(npart, dumpfile) - use krome_user, only:krome_idx_He integer, intent(in) :: npart character(len=*), intent(in) :: dumpfile integer :: hdferr - integer :: i - integer(hid_t) :: file_id, dspace_id, dset_id, group_id - integer(hsize_t) :: dims(1) + integer(hid_t) :: file_id, dspace_id, dset_id, group_id, type_id + integer(hsize_t) :: dims(2) + integer(hsize_t) :: dims_labels(1) + integer(size_t) :: str_len ! create HDF5 file print "(1x,a)",'Writing to '//trim(dumpfile)//'.h5' @@ -272,36 +265,79 @@ subroutine write_chem(npart, dumpfile) endif ! create dataspace for particle datasets - dims(1) = npart - call h5screate_simple_f(1, dims, dspace_id, hdferr) + dims(1) = krome_nmols + dims(2) = npart + call h5screate_simple_f(2, dims, dspace_id, hdferr) if (hdferr /= 0) then print*,'ERROR: Failed to create HDF5 dataspace' return endif - do i=1,krome_nmols - ! create dataset for each abundance - call h5dcreate_f(group_id, trim(adjustl(abundance_label(i))), H5T_NATIVE_DOUBLE, dspace_id, dset_id, hdferr) - ! write abundance data to dataset - call h5dwrite_f(dset_id, H5T_NATIVE_DOUBLE, abundance(i,1:npart), dims, hdferr) - if (hdferr /= 0) then - print*,'ERROR: Failed to create/write dataset for '//trim(adjustl(abundance_label(i))) - return - endif - call h5dclose_f(dset_id, hdferr) - enddo + + ! create one 2D dataset for all species abundances: + call h5dcreate_f(group_id, 'abundances', H5T_NATIVE_DOUBLE, dspace_id, dset_id, hdferr) + if (hdferr /= 0) then + print*,'ERROR: Failed to create dataset abundances' + return + endif + + call h5dwrite_f(dset_id, H5T_NATIVE_DOUBLE, abundance(:,1:npart), dims, hdferr) + if (hdferr /= 0) then + print*,'ERROR: Failed to write dataset abundances' + return + endif + + call h5dclose_f(dset_id, hdferr) + + ! store species labels so abundances row i is always identifiable + dims_labels(1) = krome_nmols + call h5screate_simple_f(1, dims_labels, dspace_id, hdferr) + if (hdferr /= 0) then + print*,'ERROR: Failed to create HDF5 dataspace for species labels' + return + endif + + call h5tcopy_f(H5T_FORTRAN_S1, type_id, hdferr) + if (hdferr /= 0) then + print*,'ERROR: Failed to create HDF5 string datatype for species labels' + return + endif + + str_len = len(abundance_label(1)) + call h5tset_size_f(type_id, str_len, hdferr) + if (hdferr /= 0) then + print*,'ERROR: Failed to set HDF5 string size for species labels' + return + endif + + call h5dcreate_f(group_id, 'species_labels', type_id, dspace_id, dset_id, hdferr) + if (hdferr /= 0) then + print*,'ERROR: Failed to create dataset species_labels' + return + endif + + call h5dwrite_f(dset_id, type_id, abundance_label, dims_labels, hdferr) + if (hdferr /= 0) then + print*,'ERROR: Failed to write dataset species_labels' + return + endif + + call h5dclose_f(dset_id, hdferr) + call h5tclose_f(type_id, hdferr) call h5sclose_f(dspace_id, hdferr) call h5gclose_f(group_id, hdferr) call h5fclose_f(file_id, hdferr) end subroutine write_chem subroutine read_chem(npart, dumpfile) - use krome_user, only:krome_idx_He integer, intent(in) :: npart character(len=*), intent(in) :: dumpfile integer :: hdferr integer :: i - integer(hid_t) :: file_id, dset_id, group_id, filespace_id - integer(hsize_t) :: max_dims(1), file_dims(1) + integer(hid_t) :: file_id, dset_id, group_id, filespace_id, type_id + integer(hsize_t) :: max_dims(2), file_dims(2) + integer(hsize_t) :: max_dims_labels(1), file_dims_labels(1) + integer(size_t) :: str_len + character(len=16) :: labels_file(krome_nmols) ! open HDF5 file print "(1x,a)",'Reading from '//trim(dumpfile)//'.h5' @@ -317,28 +353,74 @@ subroutine read_chem(npart, dumpfile) print*,'ERROR: Failed to open HDF5 group' return endif + ! open matrix dataset storing all species abundances + call h5dopen_f(group_id, 'abundance_matrix', dset_id, hdferr) + if (hdferr /= 0) then + print*,'ERROR: Failed to open dataset abundance_matrix' + return + endif + call h5dget_space_f(dset_id, filespace_id, hdferr) + call h5sget_simple_extent_dims_f(filespace_id, file_dims, max_dims, hdferr) + if (file_dims(1) /= krome_nmols .or. file_dims(2) /= npart) then + print*,'ERROR: abundance_matrix shape mismatch in HDF5 file' + return + endif + call h5dread_f(dset_id, H5T_NATIVE_DOUBLE, abundance(:,1:npart), file_dims, hdferr) + if (hdferr /= 0) then + print*,'ERROR: Failed to read dataset abundance_matrix' + return + endif - do i=1,krome_nmols - ! open dataset for each abundance - call h5dopen_f(group_id, trim(adjustl(abundance_label(i))), dset_id, hdferr) - ! check number of particles in file - call h5dget_space_f(dset_id, filespace_id, hdferr) - call h5sget_simple_extent_dims_f(filespace_id, file_dims, max_dims, hdferr) - if (file_dims(1) /= npart) then - print*,'ERROR: Number of particles in HDF5 file does not match current dump' - return - endif - ! read abundance data from dataset - call h5dread_f(dset_id, H5T_NATIVE_DOUBLE, abundance(i,1:npart), file_dims, hdferr) - if (hdferr /= 0) then - print*,'ERROR: Failed to open/read dataset for '//trim(adjustl(abundance_label(i))) + call h5dclose_f(dset_id, hdferr) + call h5sclose_f(filespace_id, hdferr) + + ! read and validate species labels to guarantee row-label mapping + call h5dopen_f(group_id, 'species_labels', dset_id, hdferr) + if (hdferr /= 0) then + print*,'ERROR: Failed to open dataset species_labels' + return + endif + + call h5dget_space_f(dset_id, filespace_id, hdferr) + call h5sget_simple_extent_dims_f(filespace_id, file_dims_labels, max_dims_labels, hdferr) + if (file_dims_labels(1) /= krome_nmols) then + print*,'ERROR: species_labels length mismatch in HDF5 file' + return + endif + + call h5tcopy_f(H5T_FORTRAN_S1, type_id, hdferr) + if (hdferr /= 0) then + print*,'ERROR: Failed to create HDF5 string datatype for species_labels read' + return + endif + + str_len = len(labels_file(1)) + call h5tset_size_f(type_id, str_len, hdferr) + if (hdferr /= 0) then + print*,'ERROR: Failed to set HDF5 string size for species_labels read' + return + endif + + call h5dread_f(dset_id, type_id, labels_file, file_dims_labels, hdferr) + if (hdferr /= 0) then + print*,'ERROR: Failed to read dataset species_labels' + return + endif + + do i = 1, krome_nmols + if (trim(adjustl(labels_file(i))) /= trim(adjustl(abundance_label(i)))) then + print*,'ERROR: species_labels mismatch at index ', i + print*,' file: ', trim(adjustl(labels_file(i))) + print*,' run : ', trim(adjustl(abundance_label(i))) return endif - call h5dclose_f(dset_id, hdferr) - call h5sclose_f(filespace_id, hdferr) enddo + + call h5dclose_f(dset_id, hdferr) + call h5sclose_f(filespace_id, hdferr) + call h5tclose_f(type_id, hdferr) call h5gclose_f(group_id, hdferr) call h5fclose_f(file_id, hdferr) end subroutine read_chem From d2eda1ede6778ea186640f3ddea86a2fb87d5e71 Mon Sep 17 00:00:00 2001 From: camillelandri Date: Wed, 8 Jul 2026 10:20:05 +0200 Subject: [PATCH 09/28] include download of HDF5 libs in CI, specify HDF5_DIR after installing it, fix gfortran compilation failure, change rundir to CI workspace --- .github/workflows/krome.yml | 20 ++++++++++++++------ 1 file changed, 14 insertions(+), 6 deletions(-) diff --git a/.github/workflows/krome.yml b/.github/workflows/krome.yml index 0e3c3c715..2331e83c7 100644 --- a/.github/workflows/krome.yml +++ b/.github/workflows/krome.yml @@ -7,7 +7,7 @@ env: PREFIX: /usr/local/ PHANTOM_DIR: ${{ github.workspace }} KROMEPATH: ${{ github.workspace }}/krome - FC: gfortran # default if not overwritten by matrix + PHANTOM_FC: gfortran # default if not overwritten by matrix # A workflow run is made up of one or more jobs that can run sequentially or in parallel jobs: @@ -31,7 +31,15 @@ jobs: - name: "Install LAPACK dependencies (Ubuntu)" if: matrix.os == 'ubuntu-latest' run: sudo apt-get update && sudo apt-get install -y liblapack-dev - + + - name: "Install HDF5 dependencies (Ubuntu)" + if: matrix.os == 'ubuntu-latest' + run: sudo apt-get update && sudo apt-get install -y libhdf5-serial-dev && echo "HDF5_DIR=/usr/lib/x86_64-linux-gnu/hdf5/serial" >> "$GITHUB_ENV" + + - name: "Install HDF5 dependencies (macOS)" + if: matrix.os == 'macos-latest' + run: brew update && brew install hdf5 && echo "HDF5_DIR=$(brew --prefix hdf5)" >> "$GITHUB_ENV" + - name: "Clone phantom" uses: actions/checkout@v4 @@ -42,13 +50,13 @@ jobs: run: echo -e "-n=networks/react_umist\n-iRHS\n-noSinkCheck\n-noRecCheck\n-noTlimits\n-unsafe\n-skipODEthermo\n-skipJacobian" > krome.setup && cat krome.setup - name: "Compile phantom and link with krome" - run: make SYSTEM=${{ env.FC }} SETUP=wind KROME=yes KROMEPATH=${KROMEPATH} PREFIX=${PREFIX} RUNDIR=${{ github.workspace }} + run: make SYSTEM=$PHANTOM_FC SETUP=wind KROME=yes RUNDIR=${{ github.workspace }} - name: "Compile phantomsetup and link with krome" - run: make SYSTEM=${{ env.FC }} SETUP=wind KROME=yes KROMEPATH=${KROMEPATH} PREFIX=${PREFIX} RUNDIR=${PREFIX} setup + run: make SYSTEM=$PHANTOM_FC SETUP=wind KROME=yes RUNDIR=${{ github.workspace }} setup - name: "Compile phantomanalysis and link with krome" - run: make SYSTEM=${{ env.FC }} SETUP=wind KROME=yes KROMEPATH=${KROMEPATH} PREFIX=${PREFIX} RUNDIR=${PREFIX} analysis + run: make SYSTEM=$PHANTOM_FC SETUP=wind KROME=yes RUNDIR=${{ github.workspace }} analysis - name: "Compile phantommoddump and link with krome" - run: make SYSTEM=${{ env.FC }} SETUP=wind KROME=yes KROMEPATH=${KROMEPATH} PREFIX=${PREFIX} RUNDIR=${PREFIX} moddump + run: make SYSTEM=$PHANTOM_FC SETUP=wind KROME=yes RUNDIR=${{ github.workspace }} moddump From 4b2685ec49d2fba04c3bc6d5e4af343c572f8e3f Mon Sep 17 00:00:00 2001 From: Camille Landri Date: Thu, 2 Jul 2026 11:47:33 +0200 Subject: [PATCH 10/28] Rename analysis_krome.f90 to analysis_krome.F90 allow for preprocessing of the file --- src/utils/{analysis_krome.f90 => analysis_krome.F90} | 0 1 file changed, 0 insertions(+), 0 deletions(-) rename src/utils/{analysis_krome.f90 => analysis_krome.F90} (100%) diff --git a/src/utils/analysis_krome.f90 b/src/utils/analysis_krome.F90 similarity index 100% rename from src/utils/analysis_krome.f90 rename to src/utils/analysis_krome.F90 From a93b5864a4d975900893677836cccaeda9a99561 Mon Sep 17 00:00:00 2001 From: Camille Landri Date: Thu, 2 Jul 2026 15:38:57 +0200 Subject: [PATCH 11/28] fix compilation with self compiled hdf5 libraries --- build/Makefile | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/build/Makefile b/build/Makefile index fc9c2e907..335123eec 100644 --- a/build/Makefile +++ b/build/Makefile @@ -405,7 +405,7 @@ ifeq (X$(HDF5LIB), X) endif FFLAGS+= -I$(HDF5INCLUDE) CCFLAGS+= -I$(HDF5INCLUDE) - LDFLAGS+= -L$(HDF5_DIR)/lib -lhdf5 -lhdf5_fortran + LDFLAGS+= -L$(HDF5LIB) -Wl,-rpath,$(HDF5LIB) -lhdf5 -lhdf5_fortran endif IDFLAGS=$(FPPFLAGS) From b9b81a590e33f182638f3ef4332e1a7ec4286e21 Mon Sep 17 00:00:00 2001 From: Camille Landri Date: Thu, 2 Jul 2026 15:38:57 +0200 Subject: [PATCH 12/28] fix compilation with self compiled hdf5 libraries --- build/Makefile | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/build/Makefile b/build/Makefile index fc9c2e907..335123eec 100644 --- a/build/Makefile +++ b/build/Makefile @@ -405,7 +405,7 @@ ifeq (X$(HDF5LIB), X) endif FFLAGS+= -I$(HDF5INCLUDE) CCFLAGS+= -I$(HDF5INCLUDE) - LDFLAGS+= -L$(HDF5_DIR)/lib -lhdf5 -lhdf5_fortran + LDFLAGS+= -L$(HDF5LIB) -Wl,-rpath,$(HDF5LIB) -lhdf5 -lhdf5_fortran endif IDFLAGS=$(FPPFLAGS) From dd01a5e99a273608e02479af7c2b3dbf2c3c0cba Mon Sep 17 00:00:00 2001 From: camillelandri Date: Fri, 3 Jul 2026 10:18:16 +0200 Subject: [PATCH 13/28] read chemistry file of first dump if it exists --- src/utils/analysis_krome.F90 | 102 +++++++++++++++++++++++++++-------- 1 file changed, 79 insertions(+), 23 deletions(-) diff --git a/src/utils/analysis_krome.F90 b/src/utils/analysis_krome.F90 index 89fd85fa6..e10c3f86c 100644 --- a/src/utils/analysis_krome.F90 +++ b/src/utils/analysis_krome.F90 @@ -57,7 +57,7 @@ subroutine do_analysis(dumpfile,num,xyzh,vxyzu,particlemass,npart,time,iunit) real :: numberdensity, T_gas, gammai, mui, AUV, xi real :: abundance_part(krome_nmols), Y(krome_nmols), column_density(npart), xyzh_copy(4,npart) real :: max_radius, radius - integer :: i, j, ierr, completed_iterations, npart_copy = 0, hdferr, i_radius = 1 + integer :: i, j, isize=0, ierr, completed_iterations, npart_copy = 0, hdferr, i_radius = 1 #ifdef __GFORTRAN__ print*, "Setting number of threads to 1 (KROME is not thread-safe when compiled with gfortran)" @@ -67,6 +67,13 @@ subroutine do_analysis(dumpfile,num,xyzh,vxyzu,particlemass,npart,time,iunit) if (.not.done_init) then done_init = .true. + ! initialize HDF5 + call h5open_f(hdferr) + if (hdferr /= 0) then + print*,'ERROR: Failed to initialize HDF5' + return + endif + print*, "initialising KROME" call krome_init() print*, "Initialised KROME" @@ -83,28 +90,29 @@ subroutine do_analysis(dumpfile,num,xyzh,vxyzu,particlemass,npart,time,iunit) allocate(iprev(maxp)) iprev = 0 - print*, "setting abundances" - !$omp parallel do default(none) & - !$omp shared(npart,xyzh,vxyzu,dt_cgs,nprev,iorig,iorig_old,iprev) & - !$omp shared(abundance,abundance_prev,particlemass,unit_density) & - !$omp shared(ieos,rho_cgs,T_gas,j) & - !$omp private(i,abundance_part) - do i=1, npart - if (.not.isdead_or_accreted(xyzh(4,i))) then - call chem_init(abundance_part) - abundance(:,i) = abundance_part - endif - enddo + !check if a file with abundances already exists, if so read it in and use it as initial abundances + inquire(file=trim(dumpfile)//'.h5', size=isize) + if (isize > 0) then + print*, "Found existing abundance file, reading in abundances" + call read_chem(npart, dumpfile) + else + print*, "setting abundances" + !$omp parallel do default(none) & + !$omp shared(npart,xyzh,vxyzu,dt_cgs,nprev,iorig,iorig_old,iprev) & + !$omp shared(abundance,abundance_prev,particlemass,unit_density) & + !$omp shared(ieos,rho_cgs,T_gas,j) & + !$omp private(i,abundance_part) + do i=1, npart + if (.not.isdead_or_accreted(xyzh(4,i))) then + call chem_init(abundance_part) + abundance(:,i) = abundance_part + endif + enddo + call write_chem(npart, dumpfile) + endif call init_eos(ieos, ierr) if (ierr /= 0) call fatal(analysistype, "Failed to initialise EOS") - ! initialize HDF5 - call h5open_f(hdferr) - if (hdferr /= 0) then - print*,'ERROR: Failed to initialize HDF5' - return - endif - else dt_cgs = (time - tprev)*utime completed_iterations = 0 @@ -192,9 +200,10 @@ subroutine do_analysis(dumpfile,num,xyzh,vxyzu,particlemass,npart,time,iunit) print*, 'Completed ', completed_iterations, ' of ', npart endif enddo outer + call write_chem(npart, dumpfile) endif - call write_chem(npart, dumpfile) + nprev = npart tprev = time iorig_old(1:npart) = iorig(1:npart) @@ -277,10 +286,57 @@ subroutine write_chem(npart, dumpfile) call h5sclose_f(dspace_id, hdferr) call h5gclose_f(group_id, hdferr) call h5fclose_f(file_id, hdferr) - - end subroutine write_chem +subroutine read_chem(npart, dumpfile) + use krome_user, only:krome_idx_He + integer, intent(in) :: npart + character(len=*), intent(in) :: dumpfile + integer :: hdferr + integer :: i + integer(hid_t) :: file_id, dset_id, group_id, filespace_id + integer(hsize_t) :: max_dims(1), file_dims(1) + + ! open HDF5 file + print "(1x,a)",'Reading from '//trim(dumpfile)//'.h5' + call h5fopen_f(trim(dumpfile)//'.h5', H5F_ACC_RDONLY_F, file_id, hdferr) + if (hdferr /= 0) then + print*,'ERROR: Failed to open HDF5 file' + return + endif + + ! open group for particle data + call h5gopen_f(file_id, 'chemistry', group_id, hdferr) + if (hdferr /= 0) then + print*,'ERROR: Failed to open HDF5 group' + return + endif + + + + do i=1,krome_nmols + ! open dataset for each abundance + call h5dopen_f(group_id, trim(adjustl(abundance_label(i))), dset_id, hdferr) + ! check number of particles in file + call h5dget_space_f(dset_id, filespace_id, hdferr) + call h5sget_simple_extent_dims_f(filespace_id, file_dims, max_dims, hdferr) + if (file_dims(1) /= npart) then + print*,'ERROR: Number of particles in HDF5 file does not match current dump' + return + endif + ! read abundance data from dataset + call h5dread_f(dset_id, H5T_NATIVE_DOUBLE, abundance(i,1:npart), file_dims, hdferr) + if (hdferr /= 0) then + print*,'ERROR: Failed to open/read dataset for '//trim(adjustl(abundance_label(i))) + return + endif + call h5dclose_f(dset_id, hdferr) + call h5sclose_f(filespace_id, hdferr) + enddo + call h5gclose_f(group_id, hdferr) + call h5fclose_f(file_id, hdferr) +end subroutine read_chem + subroutine chem_init(abundance_part) use krome_user, only:krome_idx_H2,krome_idx_He,krome_idx_CO,krome_idx_C2H2,& krome_idx_HCN,krome_idx_N2,krome_idx_SiC2,krome_idx_CS,& From 37f7a93e626d7deddebca80d7fccca7024431822 Mon Sep 17 00:00:00 2001 From: Camille Landri Date: Fri, 3 Jul 2026 10:38:00 +0200 Subject: [PATCH 14/28] add timings --- src/utils/analysis_krome.F90 | 28 +++++++++++++++++----------- 1 file changed, 17 insertions(+), 11 deletions(-) diff --git a/src/utils/analysis_krome.F90 b/src/utils/analysis_krome.F90 index e10c3f86c..54ae45061 100644 --- a/src/utils/analysis_krome.F90 +++ b/src/utils/analysis_krome.F90 @@ -56,12 +56,14 @@ subroutine do_analysis(dumpfile,num,xyzh,vxyzu,particlemass,npart,time,iunit) real :: dt_cgs, rho_cgs, rholist(maxp), Tlist(maxp), mulist(maxp), Auvlist(maxp), xilist(maxp) real :: numberdensity, T_gas, gammai, mui, AUV, xi real :: abundance_part(krome_nmols), Y(krome_nmols), column_density(npart), xyzh_copy(4,npart) - real :: max_radius, radius + real :: max_radius, radius, t_start, t_stop integer :: i, j, isize=0, ierr, completed_iterations, npart_copy = 0, hdferr, i_radius = 1 #ifdef __GFORTRAN__ print*, "Setting number of threads to 1 (KROME is not thread-safe when compiled with gfortran)" call omp_set_num_threads(1) +#else + print*, "running with ", omp_get_max_threads(), " threads" #endif if (.not.done_init) then @@ -116,8 +118,9 @@ subroutine do_analysis(dumpfile,num,xyzh,vxyzu,particlemass,npart,time,iunit) else dt_cgs = (time - tprev)*utime completed_iterations = 0 - print*, dumpfile, ": not first step data, timestep = ",dt_cgs, "npart = ",npart, "nprev = ",nprev - + print*, "not first step data, timestep = ",dt_cgs, "npart = ",npart, "nprev = ",nprev + print*, "Calculating column density" + call cpu_time(t_start) xyzmh_ptmass(iReff,1) = 2. npart_copy = npart xyzh_copy = xyzh(:,:npart) @@ -135,13 +138,11 @@ subroutine do_analysis(dumpfile,num,xyzh,vxyzu,particlemass,npart,time,iunit) endif enddo column_density = column_density + rhoh(xyzh(4,i_radius),particlemass)*unit_density * max_radius * udist + call cpu_time(t_stop) + print*, " - Took ", t_stop - t_start, " seconds" - rholist = 0. - Tlist = 0. - mulist = 0. - Auvlist = 0. - xilist = 0. - + print*, "Running KROME" + call cpu_time(t_start) !$omp parallel do default(none) & !$omp shared(npart,xyzh,vxyzu,dt_cgs,nprev,iorig,iorig_old,iprev,iverbose) & !$omp shared(abundance,abundance_prev,particlemass,unit_density,udist) & @@ -151,7 +152,6 @@ subroutine do_analysis(dumpfile,num,xyzh,vxyzu,particlemass,npart,time,iunit) outer: do i=1,npart if (.not.isdead_or_accreted(xyzh(4,i))) then - inner: do j=1,nprev if (iorig(i) == iorig_old(j)) then iprev(i) = j @@ -200,10 +200,16 @@ subroutine do_analysis(dumpfile,num,xyzh,vxyzu,particlemass,npart,time,iunit) print*, 'Completed ', completed_iterations, ' of ', npart endif enddo outer + call cpu_time(t_stop) + print*, " - Took ", (t_stop - t_start)/omp_get_max_threads(), " seconds" + + call cpu_time(t_start) call write_chem(npart, dumpfile) + call cpu_time(t_stop) + print*, " - Took ", t_stop - t_start, " seconds" endif - + ! keep track of previous abundances for next dump nprev = npart tprev = time iorig_old(1:npart) = iorig(1:npart) From 30c08b651aa6da640977659a8cc11fc9a1996bbd Mon Sep 17 00:00:00 2001 From: Camille Landri Date: Fri, 3 Jul 2026 12:04:45 +0200 Subject: [PATCH 15/28] change hdf5 layout - one dataset with chem matrix instread of nmols datasets with particle array --- src/utils/analysis_krome.F90 | 198 +++++++++++++++++++++++++---------- 1 file changed, 140 insertions(+), 58 deletions(-) diff --git a/src/utils/analysis_krome.F90 b/src/utils/analysis_krome.F90 index 54ae45061..0c5210192 100644 --- a/src/utils/analysis_krome.F90 +++ b/src/utils/analysis_krome.F90 @@ -53,10 +53,10 @@ subroutine do_analysis(dumpfile,num,xyzh,vxyzu,particlemass,npart,time,iunit) real, intent(in) :: particlemass,time real, save :: tprev = 0. integer, save :: nprev = 0 - real :: dt_cgs, rho_cgs, rholist(maxp), Tlist(maxp), mulist(maxp), Auvlist(maxp), xilist(maxp) + real :: dt_cgs, rho_cgs real :: numberdensity, T_gas, gammai, mui, AUV, xi real :: abundance_part(krome_nmols), Y(krome_nmols), column_density(npart), xyzh_copy(4,npart) - real :: max_radius, radius, t_start, t_stop + real :: max_radius, radius, tstart integer :: i, j, isize=0, ierr, completed_iterations, npart_copy = 0, hdferr, i_radius = 1 #ifdef __GFORTRAN__ @@ -119,14 +119,17 @@ subroutine do_analysis(dumpfile,num,xyzh,vxyzu,particlemass,npart,time,iunit) dt_cgs = (time - tprev)*utime completed_iterations = 0 print*, "not first step data, timestep = ",dt_cgs, "npart = ",npart, "nprev = ",nprev - print*, "Calculating column density" - call cpu_time(t_start) + print*, "Building neighbour tree..." + tstart = omp_get_wtime() xyzmh_ptmass(iReff,1) = 2. npart_copy = npart xyzh_copy = xyzh(:,:npart) call build_tree(npart_copy,npart_copy,xyzh_copy,vxyzu) + print*, " - Took ", omp_get_wtime() - tstart, " seconds" + + print*, "Calculating column density..." + tstart = omp_get_wtime() call get_all_tau(npart, nptmass, xyzmh_ptmass, xyzh, one, 5, .false., column_density) - max_radius = 0.0 do i = 1, npart if (.not.isdead_or_accreted(xyzh(4, i))) then @@ -138,16 +141,14 @@ subroutine do_analysis(dumpfile,num,xyzh,vxyzu,particlemass,npart,time,iunit) endif enddo column_density = column_density + rhoh(xyzh(4,i_radius),particlemass)*unit_density * max_radius * udist - call cpu_time(t_stop) - print*, " - Took ", t_stop - t_start, " seconds" + print*, " - Took ", omp_get_wtime() - tstart, " seconds" print*, "Running KROME" - call cpu_time(t_start) + tstart = omp_get_wtime() !$omp parallel do default(none) & !$omp shared(npart,xyzh,vxyzu,dt_cgs,nprev,iorig,iorig_old,iprev,iverbose) & - !$omp shared(abundance,abundance_prev,particlemass,unit_density,udist) & + !$omp shared(abundance,abundance_prev,particlemass,unit_density,udist,iphase) & !$omp shared(ieos,gamma,gmw,time,completed_iterations,column_density,AuvAv,albedo) & - !$omp shared(rholist,Tlist,mulist,Auvlist,xilist,iphase) & !$omp private(i,j,abundance_part,Y,rho_cgs,numberdensity,T_gas,gammai,mui,AUV,xi) outer: do i=1,npart @@ -186,12 +187,6 @@ subroutine do_analysis(dumpfile,num,xyzh,vxyzu,particlemass,npart,time,iunit) call krome(Y,T_gas,dt_cgs) abundance_part = Y/numberdensity abundance(:,i) = abundance_part - - rholist(i) = rho_cgs - Tlist(i) = T_gas - mulist(i) = mui - Auvlist(i) = AUV - xilist(i) = xi endif endif if (iverbose > 1) then @@ -200,13 +195,11 @@ subroutine do_analysis(dumpfile,num,xyzh,vxyzu,particlemass,npart,time,iunit) print*, 'Completed ', completed_iterations, ' of ', npart endif enddo outer - call cpu_time(t_stop) - print*, " - Took ", (t_stop - t_start)/omp_get_max_threads(), " seconds" - - call cpu_time(t_start) + print*, " - Took ", omp_get_wtime() - tstart, " seconds" + + tstart = omp_get_wtime() call write_chem(npart, dumpfile) - call cpu_time(t_stop) - print*, " - Took ", t_stop - t_start, " seconds" + print*, " - Took ", omp_get_wtime() - tstart, " seconds" endif ! keep track of previous abundances for next dump @@ -248,13 +241,13 @@ real function get_xi(AUV) end function get_xi subroutine write_chem(npart, dumpfile) - use krome_user, only:krome_idx_He integer, intent(in) :: npart character(len=*), intent(in) :: dumpfile integer :: hdferr - integer :: i - integer(hid_t) :: file_id, dspace_id, dset_id, group_id - integer(hsize_t) :: dims(1) + integer(hid_t) :: file_id, dspace_id, dset_id, group_id, type_id + integer(hsize_t) :: dims(2) + integer(hsize_t) :: dims_labels(1) + integer(size_t) :: str_len ! create HDF5 file print "(1x,a)",'Writing to '//trim(dumpfile)//'.h5' @@ -272,36 +265,79 @@ subroutine write_chem(npart, dumpfile) endif ! create dataspace for particle datasets - dims(1) = npart - call h5screate_simple_f(1, dims, dspace_id, hdferr) + dims(1) = krome_nmols + dims(2) = npart + call h5screate_simple_f(2, dims, dspace_id, hdferr) if (hdferr /= 0) then print*,'ERROR: Failed to create HDF5 dataspace' return endif - do i=1,krome_nmols - ! create dataset for each abundance - call h5dcreate_f(group_id, trim(adjustl(abundance_label(i))), H5T_NATIVE_DOUBLE, dspace_id, dset_id, hdferr) - ! write abundance data to dataset - call h5dwrite_f(dset_id, H5T_NATIVE_DOUBLE, abundance(i,1:npart), dims, hdferr) - if (hdferr /= 0) then - print*,'ERROR: Failed to create/write dataset for '//trim(adjustl(abundance_label(i))) - return - endif - call h5dclose_f(dset_id, hdferr) - enddo + + ! create one 2D dataset for all species abundances: + call h5dcreate_f(group_id, 'abundances', H5T_NATIVE_DOUBLE, dspace_id, dset_id, hdferr) + if (hdferr /= 0) then + print*,'ERROR: Failed to create dataset abundances' + return + endif + + call h5dwrite_f(dset_id, H5T_NATIVE_DOUBLE, abundance(:,1:npart), dims, hdferr) + if (hdferr /= 0) then + print*,'ERROR: Failed to write dataset abundances' + return + endif + + call h5dclose_f(dset_id, hdferr) + + ! store species labels so abundances row i is always identifiable + dims_labels(1) = krome_nmols + call h5screate_simple_f(1, dims_labels, dspace_id, hdferr) + if (hdferr /= 0) then + print*,'ERROR: Failed to create HDF5 dataspace for species labels' + return + endif + + call h5tcopy_f(H5T_FORTRAN_S1, type_id, hdferr) + if (hdferr /= 0) then + print*,'ERROR: Failed to create HDF5 string datatype for species labels' + return + endif + + str_len = len(abundance_label(1)) + call h5tset_size_f(type_id, str_len, hdferr) + if (hdferr /= 0) then + print*,'ERROR: Failed to set HDF5 string size for species labels' + return + endif + + call h5dcreate_f(group_id, 'species_labels', type_id, dspace_id, dset_id, hdferr) + if (hdferr /= 0) then + print*,'ERROR: Failed to create dataset species_labels' + return + endif + + call h5dwrite_f(dset_id, type_id, abundance_label, dims_labels, hdferr) + if (hdferr /= 0) then + print*,'ERROR: Failed to write dataset species_labels' + return + endif + + call h5dclose_f(dset_id, hdferr) + call h5tclose_f(type_id, hdferr) call h5sclose_f(dspace_id, hdferr) call h5gclose_f(group_id, hdferr) call h5fclose_f(file_id, hdferr) end subroutine write_chem subroutine read_chem(npart, dumpfile) - use krome_user, only:krome_idx_He integer, intent(in) :: npart character(len=*), intent(in) :: dumpfile integer :: hdferr integer :: i - integer(hid_t) :: file_id, dset_id, group_id, filespace_id - integer(hsize_t) :: max_dims(1), file_dims(1) + integer(hid_t) :: file_id, dset_id, group_id, filespace_id, type_id + integer(hsize_t) :: max_dims(2), file_dims(2) + integer(hsize_t) :: max_dims_labels(1), file_dims_labels(1) + integer(size_t) :: str_len + character(len=16) :: labels_file(krome_nmols) ! open HDF5 file print "(1x,a)",'Reading from '//trim(dumpfile)//'.h5' @@ -317,28 +353,74 @@ subroutine read_chem(npart, dumpfile) print*,'ERROR: Failed to open HDF5 group' return endif + ! open matrix dataset storing all species abundances + call h5dopen_f(group_id, 'abundance_matrix', dset_id, hdferr) + if (hdferr /= 0) then + print*,'ERROR: Failed to open dataset abundance_matrix' + return + endif + call h5dget_space_f(dset_id, filespace_id, hdferr) + call h5sget_simple_extent_dims_f(filespace_id, file_dims, max_dims, hdferr) + if (file_dims(1) /= krome_nmols .or. file_dims(2) /= npart) then + print*,'ERROR: abundance_matrix shape mismatch in HDF5 file' + return + endif + call h5dread_f(dset_id, H5T_NATIVE_DOUBLE, abundance(:,1:npart), file_dims, hdferr) + if (hdferr /= 0) then + print*,'ERROR: Failed to read dataset abundance_matrix' + return + endif - do i=1,krome_nmols - ! open dataset for each abundance - call h5dopen_f(group_id, trim(adjustl(abundance_label(i))), dset_id, hdferr) - ! check number of particles in file - call h5dget_space_f(dset_id, filespace_id, hdferr) - call h5sget_simple_extent_dims_f(filespace_id, file_dims, max_dims, hdferr) - if (file_dims(1) /= npart) then - print*,'ERROR: Number of particles in HDF5 file does not match current dump' - return - endif - ! read abundance data from dataset - call h5dread_f(dset_id, H5T_NATIVE_DOUBLE, abundance(i,1:npart), file_dims, hdferr) - if (hdferr /= 0) then - print*,'ERROR: Failed to open/read dataset for '//trim(adjustl(abundance_label(i))) + call h5dclose_f(dset_id, hdferr) + call h5sclose_f(filespace_id, hdferr) + + ! read and validate species labels to guarantee row-label mapping + call h5dopen_f(group_id, 'species_labels', dset_id, hdferr) + if (hdferr /= 0) then + print*,'ERROR: Failed to open dataset species_labels' + return + endif + + call h5dget_space_f(dset_id, filespace_id, hdferr) + call h5sget_simple_extent_dims_f(filespace_id, file_dims_labels, max_dims_labels, hdferr) + if (file_dims_labels(1) /= krome_nmols) then + print*,'ERROR: species_labels length mismatch in HDF5 file' + return + endif + + call h5tcopy_f(H5T_FORTRAN_S1, type_id, hdferr) + if (hdferr /= 0) then + print*,'ERROR: Failed to create HDF5 string datatype for species_labels read' + return + endif + + str_len = len(labels_file(1)) + call h5tset_size_f(type_id, str_len, hdferr) + if (hdferr /= 0) then + print*,'ERROR: Failed to set HDF5 string size for species_labels read' + return + endif + + call h5dread_f(dset_id, type_id, labels_file, file_dims_labels, hdferr) + if (hdferr /= 0) then + print*,'ERROR: Failed to read dataset species_labels' + return + endif + + do i = 1, krome_nmols + if (trim(adjustl(labels_file(i))) /= trim(adjustl(abundance_label(i)))) then + print*,'ERROR: species_labels mismatch at index ', i + print*,' file: ', trim(adjustl(labels_file(i))) + print*,' run : ', trim(adjustl(abundance_label(i))) return endif - call h5dclose_f(dset_id, hdferr) - call h5sclose_f(filespace_id, hdferr) enddo + + call h5dclose_f(dset_id, hdferr) + call h5sclose_f(filespace_id, hdferr) + call h5tclose_f(type_id, hdferr) call h5gclose_f(group_id, hdferr) call h5fclose_f(file_id, hdferr) end subroutine read_chem From 70f6c5dc9dba4a71c68e2f71db14f613cd3fcdaf Mon Sep 17 00:00:00 2001 From: Camille Landri Date: Wed, 8 Jul 2026 10:20:05 +0200 Subject: [PATCH 16/28] include download of HDF5 libs in CI, specify HDF5_DIR after installing it, fix gfortran compilation failure, change rundir to CI workspace --- .github/workflows/krome.yml | 20 ++++++++++++++------ 1 file changed, 14 insertions(+), 6 deletions(-) diff --git a/.github/workflows/krome.yml b/.github/workflows/krome.yml index 0e3c3c715..2331e83c7 100644 --- a/.github/workflows/krome.yml +++ b/.github/workflows/krome.yml @@ -7,7 +7,7 @@ env: PREFIX: /usr/local/ PHANTOM_DIR: ${{ github.workspace }} KROMEPATH: ${{ github.workspace }}/krome - FC: gfortran # default if not overwritten by matrix + PHANTOM_FC: gfortran # default if not overwritten by matrix # A workflow run is made up of one or more jobs that can run sequentially or in parallel jobs: @@ -31,7 +31,15 @@ jobs: - name: "Install LAPACK dependencies (Ubuntu)" if: matrix.os == 'ubuntu-latest' run: sudo apt-get update && sudo apt-get install -y liblapack-dev - + + - name: "Install HDF5 dependencies (Ubuntu)" + if: matrix.os == 'ubuntu-latest' + run: sudo apt-get update && sudo apt-get install -y libhdf5-serial-dev && echo "HDF5_DIR=/usr/lib/x86_64-linux-gnu/hdf5/serial" >> "$GITHUB_ENV" + + - name: "Install HDF5 dependencies (macOS)" + if: matrix.os == 'macos-latest' + run: brew update && brew install hdf5 && echo "HDF5_DIR=$(brew --prefix hdf5)" >> "$GITHUB_ENV" + - name: "Clone phantom" uses: actions/checkout@v4 @@ -42,13 +50,13 @@ jobs: run: echo -e "-n=networks/react_umist\n-iRHS\n-noSinkCheck\n-noRecCheck\n-noTlimits\n-unsafe\n-skipODEthermo\n-skipJacobian" > krome.setup && cat krome.setup - name: "Compile phantom and link with krome" - run: make SYSTEM=${{ env.FC }} SETUP=wind KROME=yes KROMEPATH=${KROMEPATH} PREFIX=${PREFIX} RUNDIR=${{ github.workspace }} + run: make SYSTEM=$PHANTOM_FC SETUP=wind KROME=yes RUNDIR=${{ github.workspace }} - name: "Compile phantomsetup and link with krome" - run: make SYSTEM=${{ env.FC }} SETUP=wind KROME=yes KROMEPATH=${KROMEPATH} PREFIX=${PREFIX} RUNDIR=${PREFIX} setup + run: make SYSTEM=$PHANTOM_FC SETUP=wind KROME=yes RUNDIR=${{ github.workspace }} setup - name: "Compile phantomanalysis and link with krome" - run: make SYSTEM=${{ env.FC }} SETUP=wind KROME=yes KROMEPATH=${KROMEPATH} PREFIX=${PREFIX} RUNDIR=${PREFIX} analysis + run: make SYSTEM=$PHANTOM_FC SETUP=wind KROME=yes RUNDIR=${{ github.workspace }} analysis - name: "Compile phantommoddump and link with krome" - run: make SYSTEM=${{ env.FC }} SETUP=wind KROME=yes KROMEPATH=${KROMEPATH} PREFIX=${PREFIX} RUNDIR=${PREFIX} moddump + run: make SYSTEM=$PHANTOM_FC SETUP=wind KROME=yes RUNDIR=${{ github.workspace }} moddump From e6ddc35e8b9e251db24a6a217afd90f6c4078072 Mon Sep 17 00:00:00 2001 From: camillelandri Date: Mon, 13 Jul 2026 11:54:18 +0200 Subject: [PATCH 17/28] ran the bots --- AUTHORS | 4 ++-- src/main/eos.f90 | 18 ++++++++--------- src/main/eos_zerotemp.f90 | 38 +++++++++++++++++++----------------- src/main/utils_raytracer.f90 | 10 ++++------ src/setup/setup_grtde.f90 | 6 +++--- src/tests/test_setstar.f90 | 10 ++-------- src/utils/analysis_krome.F90 | 24 +++++++++++------------ 7 files changed, 52 insertions(+), 58 deletions(-) diff --git a/AUTHORS b/AUTHORS index 2b95c98b5..00582b408 100644 --- a/AUTHORS +++ b/AUTHORS @@ -35,7 +35,9 @@ Spencer Magnall Caitlyn Hardiman Amber Tilly Steven Rieder +Ali Pourmand Sergei Biriukov +Camille Landri Giovanni Dipierro Roberto Iaconi Siméon @@ -45,10 +47,8 @@ Stephen Nielson Sahl Rowther Martina Toscani Benedetta Veronesi -Camille Landri Simon Glover Thomas Reichardt -Ali Pourmand Chunliang Mu <86601204+chunliangmu@users.noreply.github.com> Jean-François Gonzalez Christopher Russell diff --git a/src/main/eos.f90 b/src/main/eos.f90 index d75ad799b..01c4da0a0 100644 --- a/src/main/eos.f90 +++ b/src/main/eos.f90 @@ -50,8 +50,9 @@ module eos ! ! :Dependencies: dim, dump_utils, eos_HIIR, eos_barotropic, eos_gasradrec, ! eos_helmholtz, eos_idealplusrad, eos_mesa, eos_piecewise, eos_shen, -! eos_stamatellos, eos_stratified, eos_tillotson, infile_utils, io, -! ionization_mod, mesa_microphysics, part, physcon, units +! eos_stamatellos, eos_stratified, eos_tillotson, eos_zerotemp, +! infile_utils, io, ionization_mod, mesa_microphysics, part, physcon, +! units ! use part, only:ien_etotal,ien_entropy,ien_type use dim, only:gr,do_radiation @@ -529,7 +530,7 @@ subroutine equationofstate(eos_type,ponrhoi,spsoundi,rhoi,xi,yi,zi,tempi,eni,gam ponrhoi = presi/rhoi gammai = 1.d0 + presi/(eni*rhoi) spsoundi = sqrt(gammai*ponrhoi) - case (25) ! zero temperature EOS + case (25) ! zero temperature EOS cgsrhoi = rhoi * unit_density call get_zerotemp_pressure(cgsrhoi,cgspresi) @@ -538,7 +539,7 @@ subroutine equationofstate(eos_type,ponrhoi,spsoundi,rhoi,xi,yi,zi,tempi,eni,gam presi = cgspresi/unit_pressure ponrhoi = presi/rhoi spsoundi = cgsspsoundi / unit_velocity - tempi = 0. + tempi = 0. case default spsoundi = 0. ! avoids compiler warnings ponrhoi = 0. @@ -647,7 +648,6 @@ subroutine init_eos(eos_type,ierr) if (ierr_ra > 0) call warning('init_eos','Failed to read EOS file') call init_coolra - case(25) ! ! zero temperature @@ -989,7 +989,7 @@ subroutine calc_temp_and_ene(eos_type,rho,pres,ene,temp,ierr,guesseint,mu_local, case(24) ! Stamatellos temp = pres /(rho * Rg) * mu call getintenerg_opdep(temp, rho, ene) - case(25) ! zero temp eos + case(25) ! zero temp eos call get_zerotemp_u(rho,ene) temp = 0 case default @@ -1035,7 +1035,7 @@ subroutine calc_rho_from_PT(eos_type,pres,temp,rho,ierr,mu_local,X_local,Z_local rho = pres / (temp * Rg) * mu case(12) ! Ideal gas + radiation call get_idealplusrad_rhofrompresT(pres,temp,mu,rho) - case(25) ! zero temperature eos + case(25) ! zero temperature eos call get_zerotemp_rhofrompres(pres,rho,ierr) case default ierr = 1 @@ -1511,7 +1511,7 @@ logical function eos_allows_shock_and_work(ieos) integer, intent(in) :: ieos select case(ieos) - case(2,5,10,12,15,16,21,22,24,25) + case(2,5,10,12,15,16,21,22,24,25) eos_allows_shock_and_work = .true. case default eos_allows_shock_and_work = .false. @@ -1538,7 +1538,7 @@ end function eos_requires_polyk ! a non-zero pressure even if no thermal energy is set !+ !----------------------------------------------------------------------- -logical function eos_has_pressure_without_u(ieos) +logical function eos_has_pressure_without_u(ieos) integer, intent(in) :: ieos eos_has_pressure_without_u = eos_requires_isothermal(ieos) .or. & diff --git a/src/main/eos_zerotemp.f90 b/src/main/eos_zerotemp.f90 index 4c59d9cda..326ad7b19 100644 --- a/src/main/eos_zerotemp.f90 +++ b/src/main/eos_zerotemp.f90 @@ -6,18 +6,27 @@ !--------------------------------------------------------------------------! module eos_zerotemp ! -! Implements zero temperature equation of state, e.g. for white dwarfs. Meant to +! Implements zero temperature equation of state, e.g. for white dwarfs. Meant to ! ! :References: Kippenhahn & Weigert, Stellar Structure and Evolution, section 15.2 ! -! :Dependencies: physcon +! :Owner: Ali Pourmand +! +! :Runtime parameters: +! - xc : *Carbon mass fraction* +! - xh : *Hydrogen mass fraction* +! - xhe : *Helium mass fraction* +! - xmg : *Magnesium mass fraction* +! - xne : *Neon mass fraction* +! - xo : *Oxygen mass fraction* +! +! :Dependencies: infile_utils, io, physcon, units ! use units, only:unit_density,unit_velocity use physcon, only:pi,atomic_mass_unit,mass_electron_cgs,planckh,c implicit none real :: mu_e ! mean molecular weight per free electron - public :: eos_zerotemp_init,read_options_eos_zerotemp,write_options_eos_zerotemp,& f_chandra,get_zerotemp_pressure,get_zerotemp_u,get_zerotemp_spsoundi,& get_zerotemp_rhofrompres,eos_zerotemp_calc_mu_e,eos_zerotemp_eosinfo @@ -31,8 +40,8 @@ module eos_zerotemp real :: xc = 0.5 real :: xo = 0.5 real :: xne = 0.0 - real :: xmg = 0.0 - + real :: xmg = 0.0 + integer, parameter :: speciesmax = 6 character(len=10) :: speciesname(speciesmax) real :: xmass(speciesmax) ! mass fraction of species @@ -40,7 +49,6 @@ module eos_zerotemp real :: Zion(speciesmax) ! number of protons contains - !---------------------------------------------------------------- !+ ! initialise zero-temperature EOS composition @@ -48,7 +56,7 @@ module eos_zerotemp !---------------------------------------------------------------- subroutine eos_zerotemp_init(ierr) - use io, only: warning, fatal + use io, only:warning, fatal integer, intent(out) :: ierr ierr = 0 @@ -98,7 +106,7 @@ end subroutine eos_zerotemp_init !---------------------------------------------------------------- subroutine read_options_eos_zerotemp(db,nerr) - use infile_utils, only: inopts, read_inopt + use infile_utils, only:inopts, read_inopt type(inopts), intent(inout) :: db(:) integer, intent(inout) :: nerr @@ -118,7 +126,7 @@ end subroutine read_options_eos_zerotemp !---------------------------------------------------------------- subroutine write_options_eos_zerotemp(iunit) - use infile_utils, only: write_inopt + use infile_utils, only:write_inopt integer, intent(in) :: iunit call write_inopt(xh ,'xh' ,'Hydrogen mass fraction',iunit) @@ -171,8 +179,6 @@ real function f_chandra(x) result(fx) fx = x*(2*x**2 - 3)*sqrt(1+x**2) + 3*log(x + sqrt(1+x**2)) end function f_chandra - - ! ---------------------------------------------------------------- !+ ! Calculates the zero temperature pressure for a fully degenerate electron gas, i.e. a @@ -191,7 +197,6 @@ subroutine get_zerotemp_pressure(rhoi,presi) ! This is a good approximation for white dwarfs where the electrons are highly degenerate but the ions are not. ! Note also that this assumes a fully ionised gas, so mu is the mean molecular weight per free electron - ne = rhoi / (mu_e * atomic_mass_unit) x = ( (3.0 * ne * planckh**3) / & @@ -201,7 +206,6 @@ subroutine get_zerotemp_pressure(rhoi,presi) end subroutine get_zerotemp_pressure - !----------------------------------------------------------------------- !+! Inputs and outputs in cgs units !+! get internal energy from density for zero temperature EOS @@ -221,7 +225,6 @@ subroutine get_zerotemp_u(rhoi,u) u = u/rhoi !previous equation is erg/cm^3 so here I convert it end subroutine get_zerotemp_u - !---------------------------------------------------------------- !+ ! Calculates sound speed from density (derivative of f(x), analytically found), input output in cgs @@ -273,7 +276,7 @@ subroutine get_zerotemp_rhofrompres(presi,densi,ierr) do while (f_chandra(xhi) < fx) xhi = 2.0*xhi - end do + enddo ierr = 0 @@ -285,11 +288,11 @@ subroutine get_zerotemp_rhofrompres(presi,densi,ierr) xhi = xmid else xlo = xmid - end if + endif if (abs(xhi-xlo)/(xmid+1.e-300) < tolerance) exit - end do + enddo if (iter == iter_max) ierr = 1 @@ -302,5 +305,4 @@ subroutine get_zerotemp_rhofrompres(presi,densi,ierr) end subroutine get_zerotemp_rhofrompres - end module eos_zerotemp \ No newline at end of file diff --git a/src/main/utils_raytracer.f90 b/src/main/utils_raytracer.f90 index 3a1595ced..52ed8a3f3 100644 --- a/src/main/utils_raytracer.f90 +++ b/src/main/utils_raytracer.f90 @@ -1,6 +1,6 @@ !--------------------------------------------------------------------------! ! The Phantom Smoothed Particle Hydrodynamics code, by Daniel Price et al. ! -! Copyright (c) 2007-2025 The Authors (see AUTHORS) ! +! Copyright (c) 2007-2026 The Authors (see AUTHORS) ! ! See LICENCE file for usage and distribution conditions ! ! http://phantomsph.github.io/ ! !--------------------------------------------------------------------------! @@ -20,7 +20,7 @@ module raytracer ! ! :Runtime parameters: None ! -! :Dependencies: dim, healpix, kernel, neighkdtree, part, units +! :Dependencies: healpix, kernel, neighkdtree, part, units ! use healpix @@ -86,7 +86,7 @@ end subroutine get_all_tau !+ !--------------------------------------------------------------------------------- subroutine get_all_tau_single(npart, primary, Rstar, xyzh, kappa, Rinject, order, outwards, tau) - use part, only : isdead_or_accreted + use part, only:isdead_or_accreted integer, intent(in) :: npart,order real, intent(in) :: primary(3), Rstar, xyzh(:,:), kappa(:), Rinject logical, intent(in) :: outwards @@ -123,7 +123,6 @@ subroutine get_all_tau_single(npart, primary, Rstar, xyzh, kappa, Rinject, order !$omp enddo !$omp end parallel - !_---------------------------------------------- ! DETERMINE the optical depth for each particle ! using the values available on the HEALPix rays @@ -168,7 +167,7 @@ end subroutine get_all_tau_single !+ !-------------------------------------------------------------------------- subroutine get_all_tau_companion(npart, primary, Rstar, xyzh, kappa, Rinject, companion, Rcomp, order, outwards, tau) - use part, only : isdead_or_accreted + use part, only:isdead_or_accreted integer, intent(in) :: npart, order real, intent(in) :: primary(3), companion(3), kappa(:), Rstar, Rinject, xyzh(:,:), Rcomp logical, intent(in) :: outwards @@ -319,7 +318,6 @@ subroutine interpolate_tau(nsides, vec, rays_tau, rays_dist, rays_dim, tau) tau = tau / weight end subroutine interpolate_tau - !-------------------------------------------------------------------------- !+ ! Interpolation of the optical depth for an arbitrary point on the ray, diff --git a/src/setup/setup_grtde.f90 b/src/setup/setup_grtde.f90 index fb5482a02..7074af903 100644 --- a/src/setup/setup_grtde.f90 +++ b/src/setup/setup_grtde.f90 @@ -38,9 +38,9 @@ module setup ! - z2 : *pos z star 2* ! ! :Dependencies: eos, externalforces, gravwaveutils, infile_utils, io, -! kernel, mpidomain, options, orbits, part, physcon, relaxstar, -! setbinary, setorbit, setstar, setunits, setup_params, systemutils, -! timestep, units, vectorutils +! kernel, metric_tools, mpidomain, options, orbits, part, physcon, +! relaxstar, setbinary, setorbit, setstar, setunits, setup_params, +! systemutils, timestep, units, vectorutils ! use setstar, only:star_t diff --git a/src/tests/test_setstar.f90 b/src/tests/test_setstar.f90 index 89fc47f69..d5cbada09 100644 --- a/src/tests/test_setstar.f90 +++ b/src/tests/test_setstar.f90 @@ -52,7 +52,6 @@ subroutine test_setstar(ntests,npass) ! test red supergiant call test_redsupergiant(ntests,npass) - ! test white dwarf call test_whitedwarf(ntests,npass) @@ -364,11 +363,9 @@ subroutine test_redsupergiant(ntests,npass) end subroutine test_redsupergiant - - !----------------------------------------------------------------------- !+ -! test that we can successfully relax a white dwarf with the +! test that we can successfully relax a white dwarf with the ! zero temperature EoS and Helmholtz eos(Ali Pourmand, similar to test_redsupergiant) !+ !----------------------------------------------------------------------- @@ -448,14 +445,13 @@ subroutine test_whitedwarf(ntests,npass) rhozero=rhozero,npart_total=ntot,mask=i_belong,ierr=ierr,& write_files=.false.,density_error=rmserr,energy_error=ekin) - call checkval(ierr,expected_error,0,nfail(1),'set_star runs with expected ierr') call update_test_scores(ntests,nfail,npass) call check_setup(nerror,nwarn,restart=.false.) call checkval(nerror+nwarn,0,0,nfail(1),'no errors or warnings') call update_test_scores(ntests,nfail,npass) - + call checkval(rmserr,0.0,0.05,nfail(1),'error in density profile') call update_test_scores(ntests,nfail,npass) @@ -478,6 +474,4 @@ subroutine test_whitedwarf(ntests,npass) end subroutine test_whitedwarf - - end module testsetstar diff --git a/src/utils/analysis_krome.F90 b/src/utils/analysis_krome.F90 index 0c5210192..c148cf16f 100644 --- a/src/utils/analysis_krome.F90 +++ b/src/utils/analysis_krome.F90 @@ -10,12 +10,12 @@ module analysis ! ! :References: None ! -! :Owner: Mats Esseldeurs +! :Owner: Camille Landri ! ! :Runtime parameters: None ! -! :Dependencies: eos, io, krome_main, krome_user, neighkdtree, part, -! physcon, raytracer, units +! :Dependencies: eos, hdf5, io, krome_main, krome_user, neighkdtree, +! omp_lib, part, physcon, raytracer, units ! use krome_user, only:krome_nmols use part, only: maxp @@ -39,7 +39,7 @@ module analysis subroutine do_analysis(dumpfile,num,xyzh,vxyzu,particlemass,npart,time,iunit) use part, only: isdead_or_accreted, iorig, rhoh, nptmass, xyzmh_ptmass, iReff, iboundary, igas, iphase, iamtype - use neighkdtree, only: build_tree + use neighkdtree, only:build_tree use units, only: utime,unit_density,udist use physcon, only: atomic_mass_unit use eos, only: get_temperature, ieos, gamma,gmw, init_eos @@ -58,8 +58,8 @@ subroutine do_analysis(dumpfile,num,xyzh,vxyzu,particlemass,npart,time,iunit) real :: abundance_part(krome_nmols), Y(krome_nmols), column_density(npart), xyzh_copy(4,npart) real :: max_radius, radius, tstart integer :: i, j, isize=0, ierr, completed_iterations, npart_copy = 0, hdferr, i_radius = 1 - -#ifdef __GFORTRAN__ + +#ifdef __GFORTRAN__ print*, "Setting number of threads to 1 (KROME is not thread-safe when compiled with gfortran)" call omp_set_num_threads(1) #else @@ -126,7 +126,7 @@ subroutine do_analysis(dumpfile,num,xyzh,vxyzu,particlemass,npart,time,iunit) xyzh_copy = xyzh(:,:npart) call build_tree(npart_copy,npart_copy,xyzh_copy,vxyzu) print*, " - Took ", omp_get_wtime() - tstart, " seconds" - + print*, "Calculating column density..." tstart = omp_get_wtime() call get_all_tau(npart, nptmass, xyzmh_ptmass, xyzh, one, 5, .false., column_density) @@ -150,7 +150,7 @@ subroutine do_analysis(dumpfile,num,xyzh,vxyzu,particlemass,npart,time,iunit) !$omp shared(abundance,abundance_prev,particlemass,unit_density,udist,iphase) & !$omp shared(ieos,gamma,gmw,time,completed_iterations,column_density,AuvAv,albedo) & !$omp private(i,j,abundance_part,Y,rho_cgs,numberdensity,T_gas,gammai,mui,AUV,xi) - + outer: do i=1,npart if (.not.isdead_or_accreted(xyzh(4,i))) then inner: do j=1,nprev @@ -173,11 +173,11 @@ subroutine do_analysis(dumpfile,num,xyzh,vxyzu,particlemass,npart,time,iunit) numberdensity = rho_cgs / (mui * atomic_mass_unit) T_gas = get_temperature(ieos,xyzh(1:3, i),rhoh(xyzh(4,i),particlemass),vxyzu(:,i),gammai,mui) T_gas = max(T_gas,20.0d0) - + !Radiation quantities AUV = AuvAv * column_density(i) / (mui * atomic_mass_unit) / 1.87e21 xi = get_xi(AUV) - + call krome_set_user_Auv(AUV) call krome_set_user_xi(xi) call krome_set_user_alb(ALBEDO) @@ -196,12 +196,12 @@ subroutine do_analysis(dumpfile,num,xyzh,vxyzu,particlemass,npart,time,iunit) endif enddo outer print*, " - Took ", omp_get_wtime() - tstart, " seconds" - + tstart = omp_get_wtime() call write_chem(npart, dumpfile) print*, " - Took ", omp_get_wtime() - tstart, " seconds" endif - + ! keep track of previous abundances for next dump nprev = npart tprev = time From 57f83fe5997693dc3a9ccc05f470aecc7e016038 Mon Sep 17 00:00:00 2001 From: camillelandri Date: Mon, 13 Jul 2026 12:17:39 +0200 Subject: [PATCH 18/28] fix name mismatch in hdf5 groups --- src/utils/analysis_krome.F90 | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/src/utils/analysis_krome.F90 b/src/utils/analysis_krome.F90 index c148cf16f..ad1708155 100644 --- a/src/utils/analysis_krome.F90 +++ b/src/utils/analysis_krome.F90 @@ -354,22 +354,22 @@ subroutine read_chem(npart, dumpfile) return endif ! open matrix dataset storing all species abundances - call h5dopen_f(group_id, 'abundance_matrix', dset_id, hdferr) + call h5dopen_f(group_id, 'abundances', dset_id, hdferr) if (hdferr /= 0) then - print*,'ERROR: Failed to open dataset abundance_matrix' + print*,'ERROR: Failed to open dataset abundances' return endif call h5dget_space_f(dset_id, filespace_id, hdferr) call h5sget_simple_extent_dims_f(filespace_id, file_dims, max_dims, hdferr) if (file_dims(1) /= krome_nmols .or. file_dims(2) /= npart) then - print*,'ERROR: abundance_matrix shape mismatch in HDF5 file' + print*,'ERROR: abundances shape mismatch in HDF5 file' return endif call h5dread_f(dset_id, H5T_NATIVE_DOUBLE, abundance(:,1:npart), file_dims, hdferr) if (hdferr /= 0) then - print*,'ERROR: Failed to read dataset abundance_matrix' + print*,'ERROR: Failed to read dataset abundances' return endif From 1c5982e126fe07dae0d7288030b1857323a132db Mon Sep 17 00:00:00 2001 From: camillelandri Date: Mon, 13 Jul 2026 14:50:25 +0200 Subject: [PATCH 19/28] add preprocessor flags when using omp --- src/utils/analysis_krome.F90 | 19 +++++++++++++++---- 1 file changed, 15 insertions(+), 4 deletions(-) diff --git a/src/utils/analysis_krome.F90 b/src/utils/analysis_krome.F90 index ad1708155..eb695f4fa 100644 --- a/src/utils/analysis_krome.F90 +++ b/src/utils/analysis_krome.F90 @@ -21,8 +21,9 @@ module analysis use part, only: maxp use raytracer, only: get_all_tau use hdf5 - use omp_lib - + #ifdef _OPENMP + use omp_lib, only: omp_set_num_threads, omp_get_max_threads, omp_get_wtime + #endif implicit none character(len=20), parameter, public :: analysistype = 'krome' public :: do_analysis @@ -60,8 +61,10 @@ subroutine do_analysis(dumpfile,num,xyzh,vxyzu,particlemass,npart,time,iunit) integer :: i, j, isize=0, ierr, completed_iterations, npart_copy = 0, hdferr, i_radius = 1 #ifdef __GFORTRAN__ +#ifdef _OPENMP print*, "Setting number of threads to 1 (KROME is not thread-safe when compiled with gfortran)" call omp_set_num_threads(1) +#endif #else print*, "running with ", omp_get_max_threads(), " threads" #endif @@ -120,13 +123,16 @@ subroutine do_analysis(dumpfile,num,xyzh,vxyzu,particlemass,npart,time,iunit) completed_iterations = 0 print*, "not first step data, timestep = ",dt_cgs, "npart = ",npart, "nprev = ",nprev print*, "Building neighbour tree..." +#ifdef _OPENMP tstart = omp_get_wtime() +#endif xyzmh_ptmass(iReff,1) = 2. npart_copy = npart xyzh_copy = xyzh(:,:npart) call build_tree(npart_copy,npart_copy,xyzh_copy,vxyzu) +#ifdef _OPENMP print*, " - Took ", omp_get_wtime() - tstart, " seconds" - +#endif print*, "Calculating column density..." tstart = omp_get_wtime() call get_all_tau(npart, nptmass, xyzmh_ptmass, xyzh, one, 5, .false., column_density) @@ -144,7 +150,9 @@ subroutine do_analysis(dumpfile,num,xyzh,vxyzu,particlemass,npart,time,iunit) print*, " - Took ", omp_get_wtime() - tstart, " seconds" print*, "Running KROME" +#ifdef _OPENMP tstart = omp_get_wtime() +#endif !$omp parallel do default(none) & !$omp shared(npart,xyzh,vxyzu,dt_cgs,nprev,iorig,iorig_old,iprev,iverbose) & !$omp shared(abundance,abundance_prev,particlemass,unit_density,udist,iphase) & @@ -195,11 +203,14 @@ subroutine do_analysis(dumpfile,num,xyzh,vxyzu,particlemass,npart,time,iunit) print*, 'Completed ', completed_iterations, ' of ', npart endif enddo outer +#ifdef _OPENMP print*, " - Took ", omp_get_wtime() - tstart, " seconds" - tstart = omp_get_wtime() +#endif call write_chem(npart, dumpfile) +#ifdef _OPENMP print*, " - Took ", omp_get_wtime() - tstart, " seconds" +#endif endif ! keep track of previous abundances for next dump From 383041733c39ba9b1cc0643656493e7757044e52 Mon Sep 17 00:00:00 2001 From: camillelandri Date: Mon, 13 Jul 2026 15:02:08 +0200 Subject: [PATCH 20/28] more efficient loop to arrange particle indices --- src/utils/analysis_krome.F90 | 19 ++++++++++++++----- 1 file changed, 14 insertions(+), 5 deletions(-) diff --git a/src/utils/analysis_krome.F90 b/src/utils/analysis_krome.F90 index eb695f4fa..cb8efa207 100644 --- a/src/utils/analysis_krome.F90 +++ b/src/utils/analysis_krome.F90 @@ -161,12 +161,21 @@ subroutine do_analysis(dumpfile,num,xyzh,vxyzu,particlemass,npart,time,iunit) outer: do i=1,npart if (.not.isdead_or_accreted(xyzh(4,i))) then - inner: do j=1,nprev - if (iorig(i) == iorig_old(j)) then - iprev(i) = j - exit inner + if (i <= nprev) then + if (iorig(i) == iorig_old(i)) then + iprev(i) = i + j = i endif - enddo inner + endif + if (iprev(i) == 0) then + inner: do k=1,nprev + if (iorig(i) == iorig_old(k)) then + iprev(i) = k + j = k + exit inner + endif + enddo inner + endif if (j == iprev(i)) then abundance_part(:) = abundance_prev(:,iprev(i)) From 185260271f94dfebb656f8b7bc5256cda7fed89f Mon Sep 17 00:00:00 2001 From: camillelandri Date: Mon, 13 Jul 2026 17:37:39 +0200 Subject: [PATCH 21/28] remove a lot of error catching meant for debugging --- src/utils/analysis_krome.F90 | 85 ++---------------------------------- 1 file changed, 4 insertions(+), 81 deletions(-) diff --git a/src/utils/analysis_krome.F90 b/src/utils/analysis_krome.F90 index cb8efa207..eb02d8758 100644 --- a/src/utils/analysis_krome.F90 +++ b/src/utils/analysis_krome.F90 @@ -58,7 +58,7 @@ subroutine do_analysis(dumpfile,num,xyzh,vxyzu,particlemass,npart,time,iunit) real :: numberdensity, T_gas, gammai, mui, AUV, xi real :: abundance_part(krome_nmols), Y(krome_nmols), column_density(npart), xyzh_copy(4,npart) real :: max_radius, radius, tstart - integer :: i, j, isize=0, ierr, completed_iterations, npart_copy = 0, hdferr, i_radius = 1 + integer :: i, j, k, isize=0, ierr, completed_iterations, npart_copy = 0, hdferr, i_radius = 1 #ifdef __GFORTRAN__ #ifdef _OPENMP @@ -76,7 +76,7 @@ subroutine do_analysis(dumpfile,num,xyzh,vxyzu,particlemass,npart,time,iunit) call h5open_f(hdferr) if (hdferr /= 0) then print*,'ERROR: Failed to initialize HDF5' - return + fatal endif print*, "initialising KROME" @@ -279,67 +279,27 @@ subroutine write_chem(npart, dumpfile) ! create group for particle data call h5gcreate_f(file_id, 'chemistry', group_id, hdferr) - if (hdferr /= 0) then - print*,'ERROR: Failed to create HDF5 group' - return - endif - ! create dataspace for particle datasets dims(1) = krome_nmols dims(2) = npart call h5screate_simple_f(2, dims, dspace_id, hdferr) - if (hdferr /= 0) then - print*,'ERROR: Failed to create HDF5 dataspace' - return - endif - ! create one 2D dataset for all species abundances: call h5dcreate_f(group_id, 'abundances', H5T_NATIVE_DOUBLE, dspace_id, dset_id, hdferr) - if (hdferr /= 0) then - print*,'ERROR: Failed to create dataset abundances' - return - endif - call h5dwrite_f(dset_id, H5T_NATIVE_DOUBLE, abundance(:,1:npart), dims, hdferr) if (hdferr /= 0) then print*,'ERROR: Failed to write dataset abundances' return endif - call h5dclose_f(dset_id, hdferr) ! store species labels so abundances row i is always identifiable dims_labels(1) = krome_nmols call h5screate_simple_f(1, dims_labels, dspace_id, hdferr) - if (hdferr /= 0) then - print*,'ERROR: Failed to create HDF5 dataspace for species labels' - return - endif - call h5tcopy_f(H5T_FORTRAN_S1, type_id, hdferr) - if (hdferr /= 0) then - print*,'ERROR: Failed to create HDF5 string datatype for species labels' - return - endif - str_len = len(abundance_label(1)) call h5tset_size_f(type_id, str_len, hdferr) - if (hdferr /= 0) then - print*,'ERROR: Failed to set HDF5 string size for species labels' - return - endif - call h5dcreate_f(group_id, 'species_labels', type_id, dspace_id, dset_id, hdferr) - if (hdferr /= 0) then - print*,'ERROR: Failed to create dataset species_labels' - return - endif - call h5dwrite_f(dset_id, type_id, abundance_label, dims_labels, hdferr) - if (hdferr /= 0) then - print*,'ERROR: Failed to write dataset species_labels' - return - endif call h5dclose_f(dset_id, hdferr) call h5tclose_f(type_id, hdferr) @@ -369,16 +329,6 @@ subroutine read_chem(npart, dumpfile) ! open group for particle data call h5gopen_f(file_id, 'chemistry', group_id, hdferr) - if (hdferr /= 0) then - print*,'ERROR: Failed to open HDF5 group' - return - endif - ! open matrix dataset storing all species abundances - call h5dopen_f(group_id, 'abundances', dset_id, hdferr) - if (hdferr /= 0) then - print*,'ERROR: Failed to open dataset abundances' - return - endif call h5dget_space_f(dset_id, filespace_id, hdferr) call h5sget_simple_extent_dims_f(filespace_id, file_dims, max_dims, hdferr) @@ -388,53 +338,26 @@ subroutine read_chem(npart, dumpfile) endif call h5dread_f(dset_id, H5T_NATIVE_DOUBLE, abundance(:,1:npart), file_dims, hdferr) - if (hdferr /= 0) then - print*,'ERROR: Failed to read dataset abundances' - return - endif - call h5dclose_f(dset_id, hdferr) call h5sclose_f(filespace_id, hdferr) ! read and validate species labels to guarantee row-label mapping call h5dopen_f(group_id, 'species_labels', dset_id, hdferr) - if (hdferr /= 0) then - print*,'ERROR: Failed to open dataset species_labels' - return - endif - call h5dget_space_f(dset_id, filespace_id, hdferr) call h5sget_simple_extent_dims_f(filespace_id, file_dims_labels, max_dims_labels, hdferr) if (file_dims_labels(1) /= krome_nmols) then - print*,'ERROR: species_labels length mismatch in HDF5 file' - return + call fatal(analysistype, 'species_labels length mismatch in HDF5 file') endif call h5tcopy_f(H5T_FORTRAN_S1, type_id, hdferr) - if (hdferr /= 0) then - print*,'ERROR: Failed to create HDF5 string datatype for species_labels read' - return - endif - str_len = len(labels_file(1)) call h5tset_size_f(type_id, str_len, hdferr) - if (hdferr /= 0) then - print*,'ERROR: Failed to set HDF5 string size for species_labels read' - return - endif - call h5dread_f(dset_id, type_id, labels_file, file_dims_labels, hdferr) - if (hdferr /= 0) then - print*,'ERROR: Failed to read dataset species_labels' - return - endif do i = 1, krome_nmols if (trim(adjustl(labels_file(i))) /= trim(adjustl(abundance_label(i)))) then print*,'ERROR: species_labels mismatch at index ', i - print*,' file: ', trim(adjustl(labels_file(i))) - print*,' run : ', trim(adjustl(abundance_label(i))) - return + call fatal(analysistype, 'species_labels mismatch in HDF5 file') endif enddo From c0b93b6420185b2d5ded80eb559b30372641dd3d Mon Sep 17 00:00:00 2001 From: camillelandri Date: Mon, 13 Jul 2026 17:43:19 +0200 Subject: [PATCH 22/28] fatal if hdf5 isn't initialised --- src/utils/analysis_krome.F90 | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/src/utils/analysis_krome.F90 b/src/utils/analysis_krome.F90 index eb02d8758..f625fda80 100644 --- a/src/utils/analysis_krome.F90 +++ b/src/utils/analysis_krome.F90 @@ -75,8 +75,7 @@ subroutine do_analysis(dumpfile,num,xyzh,vxyzu,particlemass,npart,time,iunit) ! initialize HDF5 call h5open_f(hdferr) if (hdferr /= 0) then - print*,'ERROR: Failed to initialize HDF5' - fatal + call fatal(analysistype, 'Failed to initialize HDF5') endif print*, "initialising KROME" From 8b591fe0a4daaef1e693c5aefaf0c2bddc316f54 Mon Sep 17 00:00:00 2001 From: camillelandri Date: Mon, 13 Jul 2026 18:01:35 +0200 Subject: [PATCH 23/28] copy point mass array to avoid overwriting sink radius --- src/utils/analysis_krome.F90 | 7 ++++--- 1 file changed, 4 insertions(+), 3 deletions(-) diff --git a/src/utils/analysis_krome.F90 b/src/utils/analysis_krome.F90 index f625fda80..ac709f3c5 100644 --- a/src/utils/analysis_krome.F90 +++ b/src/utils/analysis_krome.F90 @@ -50,7 +50,7 @@ subroutine do_analysis(dumpfile,num,xyzh,vxyzu,particlemass,npart,time,iunit) krome_set_user_alb,krome_set_user_AuvAv character(len=*), intent(in) :: dumpfile integer, intent(in) :: num,npart,iunit - real, intent(in) :: xyzh(:,:),vxyzu(:,:) + real, intent(in) :: xyzh(:,:),vxyzu(:,:), xzymh_ptmass_copy(:,:) real, intent(in) :: particlemass,time real, save :: tprev = 0. integer, save :: nprev = 0 @@ -125,10 +125,11 @@ subroutine do_analysis(dumpfile,num,xyzh,vxyzu,particlemass,npart,time,iunit) #ifdef _OPENMP tstart = omp_get_wtime() #endif - xyzmh_ptmass(iReff,1) = 2. + xyzmh_ptmass_copy(:,:) = xyzmh_ptmass(:,:) !to avoid overwriting the original ptmass array + xyzmh_ptmass_copy(iReff,1) = 2. npart_copy = npart xyzh_copy = xyzh(:,:npart) - call build_tree(npart_copy,npart_copy,xyzh_copy,vxyzu) + call build_tree(npart_copy,npart_copy,xyzh_copy,vxyzu,xyzmh_ptmass_copy) #ifdef _OPENMP print*, " - Took ", omp_get_wtime() - tstart, " seconds" #endif From ac7d6af1285f0c3b153de03366684b96d2a6c63d Mon Sep 17 00:00:00 2001 From: camillelandri Date: Mon, 13 Jul 2026 18:02:29 +0200 Subject: [PATCH 24/28] copy point mass array to avoid overwriting sink radius --- src/utils/analysis_krome.F90 | 9 +++++---- 1 file changed, 5 insertions(+), 4 deletions(-) diff --git a/src/utils/analysis_krome.F90 b/src/utils/analysis_krome.F90 index f625fda80..b28e6f89c 100644 --- a/src/utils/analysis_krome.F90 +++ b/src/utils/analysis_krome.F90 @@ -50,7 +50,7 @@ subroutine do_analysis(dumpfile,num,xyzh,vxyzu,particlemass,npart,time,iunit) krome_set_user_alb,krome_set_user_AuvAv character(len=*), intent(in) :: dumpfile integer, intent(in) :: num,npart,iunit - real, intent(in) :: xyzh(:,:),vxyzu(:,:) + real, intent(in) :: xyzh(:,:),vxyzu(:,:), xzymh_ptmass_copy(:,:) real, intent(in) :: particlemass,time real, save :: tprev = 0. integer, save :: nprev = 0 @@ -125,16 +125,17 @@ subroutine do_analysis(dumpfile,num,xyzh,vxyzu,particlemass,npart,time,iunit) #ifdef _OPENMP tstart = omp_get_wtime() #endif - xyzmh_ptmass(iReff,1) = 2. + xyzmh_ptmass_copy(:,:) = xyzmh_ptmass(:,:) !to avoid overwriting the original ptmass array + xyzmh_ptmass_copy(iReff,1) = 2. npart_copy = npart xyzh_copy = xyzh(:,:npart) - call build_tree(npart_copy,npart_copy,xyzh_copy,vxyzu) + call build_tree(npart_copy,npart_copy,xyzh_copy,vxyzu,xyzmh_ptmass_copy) #ifdef _OPENMP print*, " - Took ", omp_get_wtime() - tstart, " seconds" #endif print*, "Calculating column density..." tstart = omp_get_wtime() - call get_all_tau(npart, nptmass, xyzmh_ptmass, xyzh, one, 5, .false., column_density) + call get_all_tau(npart, nptmass, xyzmh_ptmass_copy, xyzh, one, 5, .false., column_density) max_radius = 0.0 do i = 1, npart if (.not.isdead_or_accreted(xyzh(4, i))) then From 1bbfac2e00c3792683fbeaa8ab776f8f39af6f53 Mon Sep 17 00:00:00 2001 From: Camille Landri Date: Tue, 14 Jul 2026 11:28:20 +0200 Subject: [PATCH 25/28] Documentation for krome --- docs/physics/composition.rst | 200 +++++++++++++++++++++++++++++++---- 1 file changed, 177 insertions(+), 23 deletions(-) diff --git a/docs/physics/composition.rst b/docs/physics/composition.rst index e395dde06..86a550acd 100644 --- a/docs/physics/composition.rst +++ b/docs/physics/composition.rst @@ -1,10 +1,6 @@ -Composition tracking in phantom +Composition tracking in Phantom ===================================== -This file documents how to track chemical abundances and mixing in phantom. - -Composition tracking with fixed abundances --------------------------------------------- Tracking the chemical composition of the gas in phantom with fixed abundances is straightforward, since phantom is a Lagrangian code and the particle identifiers are preserved throughout the simulation. @@ -14,27 +10,185 @@ always written to the dump files in the same order, so the particle id is simply the particle index in the dump file. In the MPI code, the particle id is stored in the 'iorig' array in the dump file. -Hence composition tracking can done as a post-processing step. +Hence composition tracking can be done as a post-processing step. + +Tracking chemical abundances with KROME +===================================== + +Phantom contains routines to post-process simulation outputs with KROME, a chemical kinetics package that solves time-dependent chemical evolution according to a specified network of chemical reactions. +KROME has only been robustly tested in Phantom for post-processing, and not for on-the-fly chemistry. +The publicly available routines of Phantom+KROME are meant to be used to process wind simulations, and are not guaranteed to be accurate for other types of simulations. +Below we provide a guide on how to install KROME, and run phantomanalysis to calculate the evolution of chemical species in a simulation. + +To calculate chemical abundances in Phantom dumps, you first need to install KROME and HDF5. + +Downloading KROME +~~~~~~~~~~~~~~~~~~ +You can download KROME from the publicly available repository (https://bitbucket.org/tgrassi/krome/src/master/). For an introduction on how KROME works, follow the tutorials provided in their documentation. +Make sure to store it in a directory that will be accessible to your Phantom installation. You do not need to compile KROME right now, it will be compiled automatically when you compile Phantom with KROME support. +You should set the environment variable ``KROMEPATH`` to point to the installation directory of KROME: +:: + export KROMEPATH=/path/to/krome + +Installing HDF5 +~~~~~~~~~~~~~~~~~~ +HDF5 is a library and file format for storing large amounts of data. The KROME routines within Phantom make use of the HDF5 library to read and write output files, because chemical abundance files are typically too large +to be stored in ASCII format. +There are several ways to install HDF5, you can get it directly from your package manager (on Linux or macOS), you can download precompiled binaries from the official release on the GitHub page (https://github.com/HDFGroup/hdf5/releases/tag/2.1.1), +or you can compile it from source (also using their GitHub repository). If you work on a cluster, it is likely to already be installed as a module, and you should be able to load the HDF5 module installed on the system (for questions, contact your cluster admin). + +To decide how to install HDF5, **you must check which compiler you intend to use for Phantom**. When linking libraries, it is important that they are compiled with the same compiler and version, otherwise you may encounter linking errors when compiling your code. +Thus, if you intend to compile Phantom with gfortran, you should install HDF5 with gfortran, same goes for ifort or ifx. Pre-compiled binaries available on the HDF5 GitHub releases or your own package manager are compiled with gfortran or ifx, but not ifort. +If you are using ifort, you will need to compile HDF5 from source. +Another constraint to account for is that **KROME is only parallelized with ifort and ifx**, so if you intend to use KROME in parallel, you will need to compile HDF5 with ifort or ifx as well. +You can technically compile KROME with gfortran, but this will strongly limit the performance of the post-processing routines. + +If you already have HDF5 installed, you can check which compiler was used to compile it by running the command ``h5cc -showconfig`` in the terminal. + +Below we detail the three installation methods for gfortran, ifort, and ifx, and more specific instructions for macOS users. + +gfortran +--------- +If you're planning on compiling Phantom with gfortran, you can simply install HDF5 from your package manager, for example on Ubuntu/Debian, with the following commands: +:: + # Ubuntu/Debian + sudo apt-get install libhdf5-serial-dev + +You should then set the environment variable ``HDF5_DIR`` to point to the installation directory of HDF5, which is typically ``/usr/lib/x86_64-linux-gnu/hdf5/serial``: +:: + export HDF5_DIR=/usr/lib/x86_64-linux-gnu/hdf5/serial +You can check where it is installed by running the command ``h5cc -showconfig`` in the terminal, which will show details about the installation of HDF5. + + +ifx +--------- +You can download ifx precompiled binaries from the official release on the GitHub page (https://github.com/HDFGroup/hdf5/releases/tag/2.1.1). +Make sure to download the file containing "intel" in the filename, for example ``hdf5-2.1.1-ubuntu-2404-intel.tar.gz``. +The downloaded archive can be extracted: +:: + tar -xvf hdf5-2.1.1-ubuntu-2404-intel.tar.gz + +which will create a directory containing another archive, which needs to be extracted as well: +:: + tar -xvf HDF5-2.1.1-Linux.tar.gz +The extracted directory contains a series of subdirectories, including a ``HDF5/``, which contains all the libraries and binaries needed to link HDF5 to Phantom. +You can move this directory to a permanent location that is more accessible, and set the environment variable ``HDF5_DIR`` to point to it: +:: + export HDF5_DIR=/path/to/HDF5 + +ifort +----- +Since no precompiled binaries are available for ifort, you will need to compile HDF5 from source. You can download the source code from the official release on the GitHub page (https://github.com/HDFGroup/hdf5/releases/tag/2.1.1). +Extract the downloaded archive to a permanent location and store the location in the environment variable ``HDF5_DIR``: +:: + tar -xvf hdf5-2.1.1.tar.gz + cd /path/to/hdf5-2.1.1 + export HDF5_DIR=/path/to/hdf5-2.1.1 +Then, you can set up the compilation of HDF5 with ifort using the following commands: +:: + cmake --preset ci-StdShar-Intel -B build/ifort \ + -DCMAKE_Fortran_COMPILER=ifort \ + -DCMAKE_INSTALL_PREFIX=$HDF5_DIR \ + -DHDF5_BUILD_FORTRAN=ON \ + -DHDF5_BUILD_JAVA=OFF \ + -DHDF5_ENABLE_JNI=OFF \ + -DHDF5_BUILD_CPP_LIB=OFF \ + -DHDF5_BUILD_TOOLS=OFF \ + -DHDF5_BUILD_EXAMPLES=OFF \ + -DBUILD_TESTING=OFF \ + -DHDF5_ENABLE_PLUGIN_SUPPORT=OFF +You can then compile and install HDF5 with: +:: + cmake --build build/ifort -j + cmake --install build/ifort +Finally, you need to specify two new environment variables to point to the HDF5 libraries and binaries: +:: + export HDF5INCLUDE=/home/camille/software/hdf5/2.1.1-ifort/mod/shared + export HDF5LIB=/home/camille/software/hdf5/2.1.1-ifort/lib + +macOS +----- +On macOS, you can only compile Phantom with gfortran, since Intel compilers are not supported anymore. +This means that you can only compile Phantom and KROME with gfortran, and **you can only use KROME in serial**. +You can install HDF5 with Homebrew: +:: + brew install hdf5 +Using the following command, you can store the path to HDF5 (likely something like ``/opt/homebrew/opt/hdf5``) and automatically store it to the environment variable ``HDF5_DIR``: +:: + export HDF5_DIR=$(brew --prefix hdf5) + + +Compiling Phantom with KROME +~~~~~~~~~~~~~~~~~~ +Once both KROME and HDF5 are installed, you can compile Phantom with KROME libraries. +You can use KROME post-processing on pre-existing models, or you can create one for the occasion. + +Once you have selected the model to process, you need to create a setup file for KROME in the directory where your phantom makefile is stored. +This file, called krome.setup, contains compilation options for KROME. For wind models, it can be created using the following command: +:: + echo -e "-n=networks/react_umist\n-iRHS\n-noSinkCheck\n-noRecCheck\n-noTlimits\n-unsafe\n-skipODEthermo\n-skipJacobian" > krome.setup + +Each flag is explained in the KROME documentation, but the most important one is ``-n``, which specifies the chemical network to use. +The chemical network files are stored in the ``networks/`` directory of the KROME directory, and you may choose any network, +although **there is no guarantee of the chemical network being suitable for your simulation**. +The network we are using here, ``react_umist`` is a detailed network for winds of AGB stars based on observations and detailed 1D chemical simulations, +but you may want to explore other networks depending on your needs. + +Once the setup file is created, you can compile Phantom with KROME by running: +:: + make KROME=yes + +The compilation will only work if the paths to HDF5 and KROME directories are set up properly, so remember to set the environment variables ``KROMEPATH`` and ``HDF5_DIR`` before compiling Phantom. +If you intend to use KROME with Phantom often, you may consider adding these variables to your ``~/.bashrc`` file. + +You can then compile the KROME routines for post-processing by running: +:: + make analysis KROME=yes ANALYSIS=analysis_krome.F90 + + +.. note:: + - You need to have compiled Phantom with KROME support before attempting to compile or use the KROME analysis routines. + - Forgetting to set KROME=yes will result in a regular compilation of Phantom **without** KROME support, and the analysis routines will not compile. + - If you need to compile other phantom utilities (phantomsetup, phantommoddump, etc.) while you are using KROME analysis, you need to compile them **with** KROME support as well, otherwise you will get linking errors when trying to run the different routines. + - Switching between KROME and non-KROME compilations of Phantom is possible, but you need to recompile Phantom with the appropriate flags each time you switch and **you need to clean the build directory** with ``make clean KROME=yes`` to avoid conflicts. + - When using ``make KROME=yes``, the compilation of Phantom takes longer than usual, since it compiles KROME as well, thus you may want to compile Phantom with KROME only when you need to use the KROME routines. + - The compilation of KROME is quiet, so you will not see much output from the compilation process itself, and since it takes a few minutes to complete it may look like it is stuck. + - If the compilation fails with missing routines or variables with "krome" in their name, KROME likely failed to compile properly. Run ``make clean KROME=yes``, check your environment variables then recompile Phantom with KROME. + + +Phantomanalysis with KROME +~~~~~~~~~~~~~~~~~~ +Your simulation is now ready to be post-processed with KROME the same way you would post-process it with any other analysis routine: +:: + ./phantomanalysis dump_????? -Writing a .comp file -~~~~~~~~~~~~~~~~~~~~~ -In practice, the composition of the gas can be tracked by writing a .comp file -either once for the entire simulation, or one per dump file. The .comp file -is a simple ascii file with one row per particle, where the columns are the abundances. -You should also give the first line a header with labels for each element:: +The routines will read the first dump file, if it finds a chemistry file corresponding to the dump (dump_?????.h5), it will read the chemical abundances from that file and use them as initial conditions for the chemical evolution, otherwise it will use default initial abundances. +The default initial abundances are set in ``analysis_krome.F90`` (located in ``src/analysis/``), and correspond to AGB stellar winds abundances. - # h1, he3, he4, c12, n14, o16, ne20, mg24, si28, s32, ar36, ca40, ti44, cr48, fe52, ni56 - 7.1119142075E-01 9.3180507858E-05 2.7341075111E-01 ... - 7.1119141814E-01 9.3180815248E-05 2.7341075311E-01 ... - 7.0179434509E-01 3.6580185766E-06 2.8264795926E-01 ... - ... +In case you need to change the initial abundances, you can edit the subroutine ``set_initial_abundances`` in ``analysis_krome.F90``, recompile the analysis routines, and run ``phantomanalysis`` again. +Please copy the ``analysis_krome.F90`` file to a different name before editing it, so we can keep the original version intact as you make any changes you want to your copy. You then need to change how ``phantomanalysis`` is compiled: +:: + make analysis KROME=yes ANALYSIS=analysis_krome_mycopy.F90 -If the phantom dump is called `foo_0000`, the .comp file should be called `foo_0000.comp` -or simply `foo.comp` if the composition is constant on each particle. +If your analysis stopped before processing all the dumps, you can restart it by using a bit of shell scripting. You can store the following commands in a bash script and run it, or include them in the script you use to run ``phantomanalysis`` on your cluster: +:: + export FILE=dump_name + export DIR=/path/to/dumps -This file is automatically read by splash and used to create extra columns in the visualisation. + start=00010 #the last dump number that was processed + files=$(printf '%s\n' $DIR/"$FILE"_* \ + | awk -F_ -v s="$start" '$NF ~ /^[0-9]{5}$/ && ($NF+0)>=s' \ + | sort) -See also --------- + if [ -z "$files" ]; then + echo "No files found in $DIR matching pattern $FILE_* with number >= $start" + exit 1 + fi + ./phantomanalysis $files +This script will find all the dump files in the directory that match the pattern and have a number greater than or equal to the last processed dump (that you must indicate by hand), and run ``phantomanalysis`` on the list of files. -- :doc:`Setting up stars and tidal disruption events ` +.. WARNING:: + The KROME routines are not guaranteed to work for all types of simulations, and have only been robustly tested for post-processing wind simulations. + More importantly, no one but you can guarantee that the chemical network you are using is suitable for your simulation. + KROME developers provide a large number of chemical networks, but it is up to you to choose the one that is appropriate for your simulation, or build your own. + Please make sure to check the KROME documentation and tutorials to understand how to choose a suitable chemical network, and how to set up initial abundances and physical conditions for your simulation. From 13deae52a718b9206d492aba6406c05b32c3bc8f Mon Sep 17 00:00:00 2001 From: Camille Landri Date: Wed, 15 Jul 2026 10:35:54 +0200 Subject: [PATCH 26/28] adress the bot's comments --- docs/physics/composition.rst | 10 +++++----- src/utils/analysis_krome.F90 | 5 +++-- 2 files changed, 8 insertions(+), 7 deletions(-) diff --git a/docs/physics/composition.rst b/docs/physics/composition.rst index 86a550acd..bf5fbd709 100644 --- a/docs/physics/composition.rst +++ b/docs/physics/composition.rst @@ -43,7 +43,7 @@ If you are using ifort, you will need to compile HDF5 from source. Another constraint to account for is that **KROME is only parallelized with ifort and ifx**, so if you intend to use KROME in parallel, you will need to compile HDF5 with ifort or ifx as well. You can technically compile KROME with gfortran, but this will strongly limit the performance of the post-processing routines. -If you already have HDF5 installed, you can check which compiler was used to compile it by running the command ``h5cc -showconfig`` in the terminal. +If you already have HDF5 installed, you can check which compiler was used to compile it by running the command ``h5fc -showconfig`` in the terminal. Below we detail the three installation methods for gfortran, ifort, and ifx, and more specific instructions for macOS users. @@ -57,7 +57,7 @@ If you're planning on compiling Phantom with gfortran, you can simply install HD You should then set the environment variable ``HDF5_DIR`` to point to the installation directory of HDF5, which is typically ``/usr/lib/x86_64-linux-gnu/hdf5/serial``: :: export HDF5_DIR=/usr/lib/x86_64-linux-gnu/hdf5/serial -You can check where it is installed by running the command ``h5cc -showconfig`` in the terminal, which will show details about the installation of HDF5. +You can check where it is installed by running the command ``h5fc -showconfig`` in the terminal, which will show details about the installation of HDF5. ifx @@ -103,8 +103,8 @@ You can then compile and install HDF5 with: cmake --install build/ifort Finally, you need to specify two new environment variables to point to the HDF5 libraries and binaries: :: - export HDF5INCLUDE=/home/camille/software/hdf5/2.1.1-ifort/mod/shared - export HDF5LIB=/home/camille/software/hdf5/2.1.1-ifort/lib + export HDF5INCLUDE=$HDF5_DIR/mod/shared + export HDF5LIB=$HDF5_DIR/lib macOS ----- @@ -163,7 +163,7 @@ Your simulation is now ready to be post-processed with KROME the same way you wo ./phantomanalysis dump_????? The routines will read the first dump file, if it finds a chemistry file corresponding to the dump (dump_?????.h5), it will read the chemical abundances from that file and use them as initial conditions for the chemical evolution, otherwise it will use default initial abundances. -The default initial abundances are set in ``analysis_krome.F90`` (located in ``src/analysis/``), and correspond to AGB stellar winds abundances. +The default initial abundances are set in ``analysis_krome.F90`` (located in ``src/utils/``), and correspond to AGB stellar winds abundances. In case you need to change the initial abundances, you can edit the subroutine ``set_initial_abundances`` in ``analysis_krome.F90``, recompile the analysis routines, and run ``phantomanalysis`` again. Please copy the ``analysis_krome.F90`` file to a different name before editing it, so we can keep the original version intact as you make any changes you want to your copy. You then need to change how ``phantomanalysis`` is compiled: diff --git a/src/utils/analysis_krome.F90 b/src/utils/analysis_krome.F90 index 5e17da82f..1b48408a1 100644 --- a/src/utils/analysis_krome.F90 +++ b/src/utils/analysis_krome.F90 @@ -60,15 +60,16 @@ subroutine do_analysis(dumpfile,num,xyzh,vxyzu,particlemass,npart,time,iunit) real :: max_radius, radius, tstart integer :: i, j, k, isize=0, ierr, completed_iterations, npart_copy = 0, hdferr, i_radius = 1 -#ifdef __GFORTRAN__ + #ifdef _OPENMP +#ifdef __GFORTRAN__ print*, "Setting number of threads to 1 (KROME is not thread-safe when compiled with gfortran)" call omp_set_num_threads(1) #else print*, "running with ", omp_get_max_threads(), " threads" #endif #else - print*, "running with ", omp_get_max_threads(), " threads" + print*, "running without OpenMP" #endif if (.not.done_init) then From 17b34f90ae49665f612abf63bbc5835e9af99731 Mon Sep 17 00:00:00 2001 From: Camille Landri Date: Wed, 15 Jul 2026 10:53:30 +0200 Subject: [PATCH 27/28] fix zero-size array warning --- src/main/utils_raytracer.f90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/main/utils_raytracer.f90 b/src/main/utils_raytracer.f90 index 52ed8a3f3..0f1e3a45e 100644 --- a/src/main/utils_raytracer.f90 +++ b/src/main/utils_raytracer.f90 @@ -491,7 +491,7 @@ subroutine find_next(inpoint, h, ray, xyzh, kappa, dtaudr, distance, inext) real, intent(out) :: distance, dtaudr integer, parameter :: nmaxcache = 0 - real :: xyzcache(0,nmaxcache) + real :: xyzcache(3, max(1,nmaxcache)) integer :: nneigh, i, prev real :: dmin, vec(3), dr, raydistance, q, norm_sq From d20c1489ac66e36023d2e62abe943751bd64bde7 Mon Sep 17 00:00:00 2001 From: Camille Landri Date: Wed, 15 Jul 2026 14:53:17 +0200 Subject: [PATCH 28/28] fix call of optical depth routines in analysis_raytracer + some warnings --- src/utils/analysis_raytracer.f90 | 22 +++++++++++----------- src/utils/utils_raytracer_all.f90 | 2 +- 2 files changed, 12 insertions(+), 12 deletions(-) diff --git a/src/utils/analysis_raytracer.f90 b/src/utils/analysis_raytracer.f90 index 8b022f824..73fc2e611 100644 --- a/src/utils/analysis_raytracer.f90 +++ b/src/utils/analysis_raytracer.f90 @@ -351,11 +351,11 @@ subroutine do_analysis(dumpfile,num,xyzh,vxyzu,particlemass,npart,time,iunit) print*, 'Start calculating optical depth outward: ', trim(jstring) if (primsec(1,2) == 0. .and. primsec(2,2) == 0. .and. primsec(3,2) == 0.) then call system_clock(start) - call get_all_tau(npart2, 1, xyzmh_ptmass, xyzh2, kappa, j, tau) + call get_all_tau(npart2, 1, xyzmh_ptmass, xyzh2, kappa, j, .true., tau) call system_clock(finish) else call system_clock(start) - call get_all_tau(npart2, 2, xyzmh_ptmass, xyzh2, kappa, j, tau) + call get_all_tau(npart2, 2, xyzmh_ptmass, xyzh2, kappa, j, .true., tau) call system_clock(finish) endif timeTau = (finish-start)/1000. @@ -501,11 +501,11 @@ subroutine do_analysis(dumpfile,num,xyzh,vxyzu,particlemass,npart,time,iunit) call build_tree(npart2,npart2,xyzh2,vxyzu) if (primsec(1,2) == 0. .and. primsec(2,2) == 0. .and. primsec(3,2) == 0.) then call system_clock(start) - call get_all_tau(npart2, 1, xyzmh_ptmass, xyzh2, kappa, order, tau) + call get_all_tau(npart2, 1, xyzmh_ptmass, xyzh2, kappa, order, .true., tau) call system_clock(finish) else call system_clock(start) - call get_all_tau(npart2, 2, xyzmh_ptmass, xyzh2, kappa, order, tau) + call get_all_tau(npart2, 2, xyzmh_ptmass, xyzh2, kappa, order, .true., tau) call system_clock(finish) endif timeTau = (finish-start)/1000. @@ -521,11 +521,11 @@ subroutine do_analysis(dumpfile,num,xyzh,vxyzu,particlemass,npart,time,iunit) print*,'Start doing scaling analysis with order =',order if (primsec(1,2) == 0. .and. primsec(2,2) == 0. .and. primsec(3,2) == 0.) then call system_clock(start) - call get_all_tau(npart2, 1, xyzmh_ptmass, xyzh2, kappa, order, tau) + call get_all_tau(npart2, 1, xyzmh_ptmass, xyzh2, kappa, order, .true., tau) call system_clock(finish) else call system_clock(start) - call get_all_tau(npart2, 2, xyzmh_ptmass, xyzh2, kappa, order, tau) + call get_all_tau(npart2, 2, xyzmh_ptmass, xyzh2, kappa, order, .true., tau) call system_clock(finish) endif timeTau = (finish-start)/1000. @@ -574,11 +574,11 @@ subroutine do_analysis(dumpfile,num,xyzh,vxyzu,particlemass,npart,time,iunit) print*, 'Start calculating optical depth outward: ', trim(jstring) if (primsec(1,2) == 0. .and. primsec(2,2) == 0. .and. primsec(3,2) == 0.) then call system_clock(start) - call get_all_tau(npart2, 1, xyzmh_ptmass, xyzh2, kappa, j, tau) + call get_all_tau(npart2, 1, xyzmh_ptmass, xyzh2, kappa, j, .true., tau) call system_clock(finish) else call system_clock(start) - call get_all_tau(npart2, 2, xyzmh_ptmass, xyzh2, kappa, j, tau) + call get_all_tau(npart2, 2, xyzmh_ptmass, xyzh2, kappa, j, .true., tau) call system_clock(finish) endif timeTau = (finish-start)/1000. @@ -636,11 +636,11 @@ subroutine do_analysis(dumpfile,num,xyzh,vxyzu,particlemass,npart,time,iunit) print*,'Start calculating optical depth' if (primsec(1,2) == 0. .and. primsec(2,2) == 0. .and. primsec(3,2) == 0.) then call system_clock(start) - call get_all_tau(npart2, 1, xyzmh_ptmass, xyzh2, kappa, order, tau) + call get_all_tau(npart2, 1, xyzmh_ptmass, xyzh2, kappa, order, .true., tau) call system_clock(finish) else call system_clock(start) - call get_all_tau(npart2, 2, xyzmh_ptmass, xyzh2, kappa, order, tau) + call get_all_tau(npart2, 2, xyzmh_ptmass, xyzh2, kappa, order, .true., tau) call system_clock(finish) endif timeTau = (finish-start)/1000. @@ -668,7 +668,7 @@ subroutine do_analysis(dumpfile,num,xyzh,vxyzu,particlemass,npart,time,iunit) order = 7 print*, 'Start calculating optical depth outward, order=',order call system_clock(start) - call get_all_tau(npart2, 2, xyzmh_ptmass, xyzh2, kappa, order, tau) + call get_all_tau(npart2, 2, xyzmh_ptmass, xyzh2, kappa, order, .true., tau) call system_clock(finish) timeTau = (finish-start)/1000. print*,'Time = ',timeTau,' seconds.' diff --git a/src/utils/utils_raytracer_all.f90 b/src/utils/utils_raytracer_all.f90 index 966b5f98d..4d3b8f6ef 100644 --- a/src/utils/utils_raytracer_all.f90 +++ b/src/utils/utils_raytracer_all.f90 @@ -1064,7 +1064,7 @@ subroutine get_tau_inwards(point, primary, xyzh, neighbors, kappa, Rstar, tau) integer :: i, next, previous, nneigh integer, parameter :: nmaxcache = 0 - real :: xyzcache(0,nmaxcache) + real :: xyzcache(3, max(1,nmaxcache)) real :: ray(3), nextDist, previousDist, maxDist, dtaudr, previousdtaudr, nextdtaudr ray = primary - xyzh(1:3,point)