From 5b88083a54479137ff992fa6aa7c9f3b5b4df671 Mon Sep 17 00:00:00 2001 From: Frederick Groth Date: Wed, 27 Aug 2025 11:24:49 +0200 Subject: [PATCH 01/16] add option for emissivity weighting --- src/Makefile | 4 ++-- src/particle_data.f | 8 +++++++- src/particles.f | 25 ++++++++++++++++++++++--- src/reader.f | 10 ++++++++-- src/readers/arepo_hdf5.f | 21 ++++++++++++++++++++- src/readers/gadget_unformatted.f | 18 +++++++++++++++++- 6 files changed, 76 insertions(+), 10 deletions(-) diff --git a/src/Makefile b/src/Makefile index 1a92023..abaf845 100644 --- a/src/Makefile +++ b/src/Makefile @@ -8,7 +8,7 @@ $(info FFTW=$(FFTW)) FILTER := $(if $(FILTER),$(FILTER),0) $(info FILTER=$(FILTER)) -# weighting scheme: number-weighted (0), mass-weighted (1), or volume-weighted (2) +# weighting scheme: number-weighted (0), mass-weighted (1), volume-weighted (2), or emissivity-weighted (3) WEIGHT := $(if $(WEIGHT),$(WEIGHT),0) $(info WEIGHT=$(WEIGHT)) @@ -16,7 +16,7 @@ $(info WEIGHT=$(WEIGHT)) KERNEL := $(if $(KERNEL),$(KERNEL),0) $(info KERNEL=$(KERNEL)) -# weighting scheme for the multiscale filter (0: volume-weighted; 1: mass-weighted) +# weighting scheme for the multiscale filter (0: volume-weighted; 1: mass-weighted, 2: emissivity-weighted) WEIGHT_FILTER := $(if $(WEIGHT_FILTER),$(WEIGHT_FILTER),0) $(info WEIGHT_FILTER=$(WEIGHT_FILTER)) diff --git a/src/particle_data.f b/src/particle_data.f index 7f6d1ab..9dc4cdc 100644 --- a/src/particle_data.f +++ b/src/particle_data.f @@ -20,7 +20,13 @@ MODULE PARTICLE_DATA #else ! Dummy variable REAL VOL -#endif +#endif + +#if weight_scheme == 3 + REAL,ALLOCATABLE::EMISSIVITY(:) +#else + REAL EMISSIVITY +#endif INTEGER,ALLOCATABLE::LIHAL(:) INTEGER,ALLOCATABLE::LIHAL_IX(:),LIHAL_JY(:),LIHAL_KZ(:) diff --git a/src/particles.f b/src/particles.f index 60fd640..7f9f2f3 100644 --- a/src/particles.f +++ b/src/particles.f @@ -2484,6 +2484,7 @@ SUBROUTINE INTERPOLATE_VELOCITIES(NX,NY,NZ,NL,NPATCH,PARE, !$OMP PARALLEL DO SHARED(NX,NY,NZ,STEP,I1,I2,J1,J2,K1,K2,RADX,RADY, !$OMP+ RADZ,U2DM,U3DM,U4DM,L0,U2,U3,U4,MASAP,VOL, +!$OMP+ EMISSIVITY, !$OMP+ KNEIGHBOURS,DX,VISC0,ABVC,TREE,XTREE, !$OMP+ YTREE,ZTREE,PI,FLAG_MASS), !$OMP+ PRIVATE(IX,JY,KZ,DIST,NEIGH,CONTA,H_KERN,BAS8, @@ -2532,6 +2533,10 @@ SUBROUTINE INTERPOLATE_VELOCITIES(NX,NY,NZ,NL,NPATCH,PARE, DO I=1,CONTA DIST(I)=DIST(I)*VOL(NEIGH(I)) END DO +#elif weight_scheme == 3 + DO I=1,CONTA + DIST(I)=DIST(I)*EMISSIVITY(NEIGH(I)) + END DO #endif BAS8=0.D0 @@ -2652,7 +2657,7 @@ SUBROUTINE INTERPOLATE_VELOCITIES(NX,NY,NZ,NL,NPATCH,PARE, !$OMP+ JJ2,KK1,KK2,RADX,RADY,RADZ,U2DM,U3DM, !$OMP+ U4DM,L0,U2,U3,U4,KNEIGHBOURS,DX,VISC0,ABVC, !$OMP+ TREE,XTREE,YTREE,ZTREE,FLAG_MASS, -!$OMP+ MASAP,VOL,PI), +!$OMP+ MASAP,VOL,EMISSIVITY,PI), !$OMP+ PRIVATE(IX,JY,KZ,DIST,NEIGH,CONTA,H_KERN,BAS8, !$OMP+ BAS8X,BAS8Y,BAS8Z,BAS8M,I,BASMASS,DA,SEARCH), !$OMP+ SCHEDULE(DYNAMIC)!, DEFAULT(NONE) @@ -2699,6 +2704,10 @@ SUBROUTINE INTERPOLATE_VELOCITIES(NX,NY,NZ,NL,NPATCH,PARE, DO I=1,CONTA DIST(I)=DIST(I)*VOL(NEIGH(I)) END DO +#elif weight_scheme == 3 + DO I=1,CONTA + DIST(I)=DIST(I)*EMISSIVITY(NEIGH(I)) + END DO #endif BAS8=0.D0 @@ -2819,7 +2828,8 @@ SUBROUTINE INTERPOLATE_VELOCITIES(NX,NY,NZ,NL,NPATCH,PARE, !$OMP PARALLEL DO SHARED(NX,NY,NZ,II1,II2,JJ1,JJ2,KK1,KK2,RADX,RADY, !$OMP+ RADZ,U2DM,U3DM,U4DM,L0,U2,U3,U4, !$OMP+ KNEIGHBOURS,DX,CR0AMR,VISC0,ABVC, -!$OMP+ TREE,XTREE,YTREE,ZTREE,FLAG_MASS,PI,MASAP,VOL), +!$OMP+ TREE,XTREE,YTREE,ZTREE,FLAG_MASS,PI,MASAP,VOL, +!$OMP+ EMISSIVITY), !$OMP+ PRIVATE(IX,JY,KZ,DIST,NEIGH,CONTA,H_KERN,BAS8, !$OMP+ BAS8X,BAS8Y,BAS8Z,BAS8M,I,BASMASS,DA,SEARCH), !$OMP+ SCHEDULE(DYNAMIC)!, DEFAULT(NONE) @@ -2861,6 +2871,10 @@ SUBROUTINE INTERPOLATE_VELOCITIES(NX,NY,NZ,NL,NPATCH,PARE, DO I=1,CONTA DIST(I)=DIST(I)*VOL(NEIGH(I)) END DO +#elif weight_scheme == 3 + DO I=1,CONTA + DIST(I)=DIST(I)*EMISSIVITY(NEIGH(I)) + END DO #endif BAS8=0.D0 @@ -2910,7 +2924,8 @@ SUBROUTINE INTERPOLATE_VELOCITIES(NX,NY,NZ,NL,NPATCH,PARE, !$OMP PARALLEL DO SHARED(LOW1,LOW2,PATCHNX,PATCHNY,PATCHNZ,CR0AMR1, !$OMP+ RX,RY,RZ,U2DM,U3DM,U4DM,L1,U12,U13,U14, !$OMP+ KNEIGHBOURS,DXPA,DX,VISC1,ABVC,SOLAP, -!$OMP+ TREE,XTREE,YTREE,ZTREE,FLAG_MASS,PI,MASAP,VOL), +!$OMP+ TREE,XTREE,YTREE,ZTREE,FLAG_MASS,PI,MASAP,VOL, +!$OMP+ EMISSIVITY), !$OMP+ PRIVATE(IPATCH,N1,N2,N3,IX,JY,KZ,DIST,NEIGH,DA, !$OMP+ CONTA,H_KERN,BAS8,BAS8X,BAS8Y,BAS8Z,BAS8M,I, !$OMP+ BASMASS,SEARCH,DO_CELL), @@ -2973,6 +2988,10 @@ SUBROUTINE INTERPOLATE_VELOCITIES(NX,NY,NZ,NL,NPATCH,PARE, DO I=1,CONTA DIST(I)=DIST(I)*VOL(NEIGH(I)) END DO +#elif weight_scheme == 3 + DO I=1,CONTA + DIST(I)=DIST(I)*EMISSIVITY(NEIGH(I)) + END DO #endif BAS8=0.D0 diff --git a/src/reader.f b/src/reader.f index 9ce4556..7fd86b6 100644 --- a/src/reader.f +++ b/src/reader.f @@ -155,6 +155,8 @@ subroutine read_particles(iter,files_per_snap,nx,ny,nz,t,zeta, #if weight_scheme == 2 deallocate(vol) +#elif weight_scheme == 3 + deallocate(emissivity) #endif end if @@ -167,7 +169,9 @@ subroutine read_particles(iter,files_per_snap,nx,ny,nz,t,zeta, #endif #if weight_scheme == 2 - allocate(vol(parti)) + allocate(vol(parti)) +#elif weight_scheme == 3 + allocate(emissivity(parti)) #endif npart(:)=0 @@ -192,7 +196,9 @@ subroutine read_particles(iter,files_per_snap,nx,ny,nz,t,zeta, !$omp parallel do shared(vol, masap, parti), private(i), default(none) do i=1,parti vol(i)=masap(i)/vol(i) - end do + end do +#elif weight_scheme == 3 + ! this is claculated already during the reading process. #endif npart(0)=low2 !retrocompatibility with general reader diff --git a/src/readers/arepo_hdf5.f b/src/readers/arepo_hdf5.f index a39ad11..8b87435 100644 --- a/src/readers/arepo_hdf5.f +++ b/src/readers/arepo_hdf5.f @@ -250,6 +250,25 @@ subroutine read_arepo_hdf5(iter, files_per_snap, vol(low1:low2)=scr4(1:numpart_thisfile(1)) call h5dclose_f(attr_id, status) deallocate(scr4) +#elif weight_scheme == 3 + allocate(scr4(numpart_thisfile(1))) + write(*,*) 'reading density ...' + call h5dopen_f(group_id, "density", attr_id, status) + call h5dget_type_f(attr_id, memtype_id, status) + call h5dread_f(attr_id, memtype_id, scr4, dims1d, status) + emissivity(low1:low2)=scr4(1:numpart_thisfile(1)) + write(*,*) 'reading internal energy ...' + call h5dopen_f(group_id, "InternalEnergy", attr_id, status) + call h5dget_type_f(attr_id, memtype_id, status) + call h5dread_f(attr_id, memtype_id, scr4, dims1d, status) + do i=1,numpart_thisfile(1) + ! calculate the weight as rho^2*sqrt(T) + ! as this is only used as weight, we don't convert values to actual temperatures! + emissivity(low1+i-1) = emissivity(low1+i-1)* + & emissivity(low1+i-1)*sqrt(scr4(i)) + end do + call h5dclose_f(attr_id, status) + deallocate(scr4) #endif call h5gclose_f(group_id, status) @@ -259,4 +278,4 @@ subroutine read_arepo_hdf5(iter, files_per_snap, return end -*********************************************************************** \ No newline at end of file +*********************************************************************** diff --git a/src/readers/gadget_unformatted.f b/src/readers/gadget_unformatted.f index 38a129c..a0e06f2 100644 --- a/src/readers/gadget_unformatted.f +++ b/src/readers/gadget_unformatted.f @@ -63,7 +63,7 @@ subroutine read_gadget_unformatted(iter, files_per_snap, integer npart_gadget(6), nall(6), blocksize real*8 massarr(6) - integer basint1,basint2,basint3,basint4 + integer basint1,basint2,basint3,basint4,i real*8 bas81,bas82,bas83,bas84,bas85,bas86 integer low1,low2 real*4,allocatable::scr4(:) @@ -148,6 +148,22 @@ subroutine read_gadget_unformatted(iter, files_per_snap, write(*,*) ' found for ',(blocksize-8)/4,' particles' vol(low1:low2)=scr4(1:npart_gadget(1)) deallocate(scr4) +#elif weight_scheme == 3 + allocate(scr4(sum(npart_gadget(1:6)))) + write(*,*) 'reading density ...' + call read_float(fil2,'RHO ',scr4,blocksize) + write(*,*) ' found for ',(blocksize-8)/4,' particles' + emissivity(low1:low2)=scr4(1:npart_gadget(1)) + write(*,*) 'reading internal energy ...' + call read_float(fil2,'U ',scr4,blocksize) + write(*,*) ' found for ',(blocksize-8)/4,' particles' + do i=1,sum(npart_gadget(1:6)) + ! calculate the weight as rho^2*sqrt(T) + ! as this is only used as weight, we don't convert values to actual temperatures! + emissivity(low1+i-1) = emissivity(low1+i-1)* + & emissivity(low1+i-1)*sqrt(scr4(i)) + end do + deallocate(scr4) #endif end do !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ From 9408658507a3d9b5adab38c2af1e08907dc779b3 Mon Sep 17 00:00:00 2001 From: Frederick Groth Date: Wed, 27 Aug 2025 13:23:42 +0200 Subject: [PATCH 02/16] add missing h5close --- src/readers/arepo_hdf5.f | 1 + 1 file changed, 1 insertion(+) diff --git a/src/readers/arepo_hdf5.f b/src/readers/arepo_hdf5.f index 8b87435..d0c838f 100644 --- a/src/readers/arepo_hdf5.f +++ b/src/readers/arepo_hdf5.f @@ -257,6 +257,7 @@ subroutine read_arepo_hdf5(iter, files_per_snap, call h5dget_type_f(attr_id, memtype_id, status) call h5dread_f(attr_id, memtype_id, scr4, dims1d, status) emissivity(low1:low2)=scr4(1:numpart_thisfile(1)) + call h5dclose_f(attr_id, status) write(*,*) 'reading internal energy ...' call h5dopen_f(group_id, "InternalEnergy", attr_id, status) call h5dget_type_f(attr_id, memtype_id, status) From a6204b5978b399f2a60d9017515d82bc8c5c0a3e Mon Sep 17 00:00:00 2001 From: Frederick Groth Date: Wed, 27 Aug 2025 13:26:32 +0200 Subject: [PATCH 03/16] add parallelization --- src/readers/arepo_hdf5.f | 2 ++ 1 file changed, 2 insertions(+) diff --git a/src/readers/arepo_hdf5.f b/src/readers/arepo_hdf5.f index d0c838f..fc2cb76 100644 --- a/src/readers/arepo_hdf5.f +++ b/src/readers/arepo_hdf5.f @@ -262,6 +262,8 @@ subroutine read_arepo_hdf5(iter, files_per_snap, call h5dopen_f(group_id, "InternalEnergy", attr_id, status) call h5dget_type_f(attr_id, memtype_id, status) call h5dread_f(attr_id, memtype_id, scr4, dims1d, status) +!$omp parallel do shared(numpart_thisfile, scr4, emissivity, low1), +!$omp+ private(i), default(none) do i=1,numpart_thisfile(1) ! calculate the weight as rho^2*sqrt(T) ! as this is only used as weight, we don't convert values to actual temperatures! From 3ce5b17bf39d266dc637ff65764f25743f793eae Mon Sep 17 00:00:00 2001 From: Frederick Groth Date: Wed, 27 Aug 2025 14:40:39 +0200 Subject: [PATCH 04/16] add parallelization --- src/readers/gadget_unformatted.f | 2 ++ 1 file changed, 2 insertions(+) diff --git a/src/readers/gadget_unformatted.f b/src/readers/gadget_unformatted.f index a0e06f2..c374f40 100644 --- a/src/readers/gadget_unformatted.f +++ b/src/readers/gadget_unformatted.f @@ -157,6 +157,8 @@ subroutine read_gadget_unformatted(iter, files_per_snap, write(*,*) 'reading internal energy ...' call read_float(fil2,'U ',scr4,blocksize) write(*,*) ' found for ',(blocksize-8)/4,' particles' +!$omp parallel do shared(npart_gadget, scr4, emissivity, low1), +!$omp+ private(i), default(none) do i=1,sum(npart_gadget(1:6)) ! calculate the weight as rho^2*sqrt(T) ! as this is only used as weight, we don't convert values to actual temperatures! From c4e49cce27815b8de3c28520296b1b4adc36d0b9 Mon Sep 17 00:00:00 2001 From: Frederick Groth Date: Wed, 27 Aug 2025 14:56:10 +0200 Subject: [PATCH 05/16] also add filter weight option for emission weighting --- src/filter.f | 25 ++++++++++++++ src/readers/fixgrid_hdf5.f | 45 ++++++++++++++++++++++++- src/readers/masclet.f | 69 ++++++++++++++++++++++++++++++++++++-- 3 files changed, 135 insertions(+), 4 deletions(-) diff --git a/src/filter.f b/src/filter.f index 0e3719e..f1f6542 100644 --- a/src/filter.f +++ b/src/filter.f @@ -1013,6 +1013,10 @@ subroutine multiscale_filter(nx,ny,nz,nl,npatch,pare, REAL DENSI_IN0(0:nmax+1,0:nmay+1,0:nmaz+1) REAL DENSI_IN1(NAMRX,NAMRY,NAMRZ,NPALEV) COMMON /DENSI/ DENSI_IN0,DENSI_IN1 +#elif weight_filter == 2 + REAL EMISS_IN0(0:nmax+1,0:nmay+1,0:nmaz+1) + REAL EMISS_IN1(NAMRX,NAMRY,NAMRZ,NPALEV) + COMMON /EMISS/ EMISS_IN0,EMISS_IN1 #endif * Auxiliary variables @@ -1078,6 +1082,27 @@ subroutine multiscale_filter(nx,ny,nz,nl,npatch,pare, dens1(:,:,:,i) = densi_in1(:,:,:,i) * dens1(:,:,:,i) end do end do +#elif weight_filter == 2 +!$omp parallel do shared(nx,ny,nz,emiss_in0,dens0), +!$omp+ private(i,j,k), +!$omp+ default(none) + do k=1,nz + do j=1,ny + do i=1,nx + dens0(i,j,k) = emiss_in0(i,j,k) * dens0(i,j,k) + end do + end do + end do + + do ir=1,nl + low1=sum(npatch(0:ir-1))+1 + low2=sum(npatch(0:ir)) +!$omp parallel do shared(low1,low2,emiss_in1,dens1,ir), +!$omp+ private(i), default(none) + do i=low1,low2 + dens1(:,:,:,i) = emiss_in1(:,:,:,i) * dens1(:,:,:,i) + end do + end do #endif write(*,*) 'in filter' diff --git a/src/readers/fixgrid_hdf5.f b/src/readers/fixgrid_hdf5.f index d98a6d2..5eacff3 100644 --- a/src/readers/fixgrid_hdf5.f +++ b/src/readers/fixgrid_hdf5.f @@ -109,6 +109,10 @@ subroutine read_fixgrid_hdf5_fields(iter, flag_filter, nx, ny, real dens0(0:nmax+1,0:nmay+1,0:nmaz+1) real dens1(namrx,namry,namrz,npalev) common /densi/ dens0,dens1 +#elif weight_filter == 2 + real emis0(0:nmax+1,0:nmay+1,0:nmaz+1) + real emis1(namrx,namry,namrz,npalev) + common /emiss/ emis0,emis1 #endif ! I assume input data is simple precision. @@ -149,6 +153,16 @@ subroutine read_fixgrid_hdf5_fields(iter, flag_filter, nx, ny, end do end do end do +#elif weight_filter == 2 +!$omp parallel do shared(emis0,nx,ny,nz), +!$omp+ private(ix,jy,kz), default(none) + do kz=0,nz+1 + do jy=0,ny+1 + do ix=0,nx+1 + emis0(ix,jy,kz) = 0.0 + end do + end do + end do #endif #endif @@ -217,6 +231,31 @@ subroutine read_fixgrid_hdf5_fields(iter, flag_filter, nx, ny, call CtoF_order_3d_real_parallel(nx,ny,nz,scr3d) dens0(1:nx,1:ny,1:nz) = scr3d(1:nx,1:ny,1:nz) call h5dclose_f(attr_id, status) +#elif weight_filter == 2 +* Read density and internal energy (if necessary: filter + emissivity-weighting) + write(*,*) 'Reading density' + call h5dopen_f(file_id, "density", attr_id, status) + call h5dget_type_f(attr_id, memtype_id, status) + call h5dread_f(attr_id, memtype_id, scr3d, dims3d, status) + call CtoF_order_3d_real_parallel(nx,ny,nz,scr3d) + emis0(1:nx,1:ny,1:nz) = scr3d(1:nx,1:ny,1:nz) + call h5dclose_f(attr_id, status) + write(*,*) 'Reading internal energy' + call h5dopen_f(file_id, "InternalEnergy", attr_id, status) + call h5dget_type_f(attr_id, memtype_id, status) + call h5dread_f(attr_id, memtype_id, scr3d, dims3d, status) + call CtoF_order_3d_real_parallel(nx,ny,nz,scr3d) +!$omp parallel do shared(emis0,scr3d,nx,ny,nz), +!$omp+ private(ix,jy,kz), default(none) + do kz=0,nz+1 + do jy=0,ny+1 + do ix=0,nx+1 + emis0(ix,jy,kz) = emis(ix,iy,iz)*emis(ix,iy,iz)* + & sqrt(scr3d(ix,iy,iz)) + end do + end do + end do + call h5dclose_f(attr_id, status) #endif end if @@ -244,6 +283,10 @@ subroutine read_fixgrid_hdf5_fields(iter, flag_filter, nx, ny, call p_minmax_ir(dens0,dens1,1,0,nx,ny,nz,nl,patchnx,patchny, & patchnz,npatch,0,basx,basy) write(*,*) 'dens min,max',basx,basy +#elif weight_filter == 2 + call p_minmax_ir(emis0,emis1,1,0,nx,ny,nz,nl,patchnx,patchny, + & patchnz,npatch,0,basx,basy) + write(*,*) 'emis min,max',basx,basy #endif write(*,*) end if @@ -337,4 +380,4 @@ subroutine write_trivial_gridsfile(iter,nx,ny,nz,nl,npatch, close(23) end -*********************************************************************** \ No newline at end of file +*********************************************************************** diff --git a/src/readers/masclet.f b/src/readers/masclet.f index 6141419..a4e3c24 100644 --- a/src/readers/masclet.f +++ b/src/readers/masclet.f @@ -148,6 +148,10 @@ subroutine read_masclet_fields(iter, flag_filter, nx, ny, nz, real dens0(0:nmax+1,0:nmay+1,0:nmaz+1) real dens1(namrx,namry,namrz,npalev) common /densi/ dens0,dens1 +#elif weight_filter == 2 + real emissivity0(0:nmax+1,0:nmay+1,0:nmaz+1) + real emissivity1(namrx,namry,namrz,npalev) + common /emiss/ emissivity0,emissivity1 #endif integer is_mascletB @@ -187,7 +191,17 @@ subroutine read_masclet_fields(iter, flag_filter, nx, ny, nz, dens0(ix,jy,kz) = 0.0 end do end do + end do +#elif weight_filter == 2 +!$omp parallel do shared(emissivity0,nx,ny,nz), +!$omp+ private(ix,jy,kz), default(none) + do kz=0,nz+1 + do jy=0,ny+1 + do ix=0,nx+1 + emissivity0(ix,jy,kz) = 0.0 end do + end do + end do #endif #endif @@ -215,6 +229,12 @@ subroutine read_masclet_fields(iter, flag_filter, nx, ny, nz, do ip = low1,low2 dens1(:,:,:,ip) = 0.0 end do +#elif weight_filter == 2 +!$omp parallel do shared(emissivity1,low1,low2) +!$omp+ private(ip), default(none) + do ip = low1,low2 + emissivity1(:,:,:,ip) = 0.0 + end do #endif #endif @@ -231,7 +251,10 @@ subroutine read_masclet_fields(iter, flag_filter, nx, ny, nz, #if weight_filter == 1 read(31) (((scr4(i,j,k),i=1,nx),j=1,ny),k=1,nz) dens0(1:nx,1:ny,1:nz) = 1.0 + scr4(1:nx,1:ny,1:nz) -#elif weight_filter == 0 +#elif weight_filter == 2 + read(31) (((scr4(i,j,k),i=1,nx),j=1,ny),k=1,nz) + emissivity0(1:nx,1:ny,1:nz) = scr4(1:nx,1:ny,1:nz) +#else read(31) #endif @@ -245,7 +268,21 @@ subroutine read_masclet_fields(iter, flag_filter, nx, ny, nz, read(31) ! pressure read(31) ! pot read(31) ! opot +#if weight_filter == 2 + read(31) (((scr4(i,j,k),i=1,nx),j=1,ny),k=1,nz) +!$omp parallel do shared(emissivity0,scr4,nx,ny,nz), +!$omp+ private(ix,jy,kz), default(none) + do kz=0,nz+1 + do jy=0,ny+1 + do ix=0,nx+1 + emissivity0(ix,jy,kz) = emissivity0(ix,jy,kz)* + & emissivity0(ix,jy,kz)*sqrt(scr4(ix,jy,kz) + end do + end do + end do +#else read(31) ! temp +#endif read(31) ! metal read(31) ! cr0 if (is_mascletB.eq.1) then @@ -267,7 +304,10 @@ subroutine read_masclet_fields(iter, flag_filter, nx, ny, nz, #if weight_filter == 1 read(31) (((scr4(i,j,k),i=1,n1),j=1,n2),k=1,n3) dens1(1:n1,1:n2,1:n3,ip) = 1.0 + scr4(1:n1,1:n2,1:n3) -#elif weight_filter == 0 +#elif weight_filter == 2 + read(31) (((scr4(i,j,k),i=1,n1),j=1,n2),k=1,n3) + emissivity1(1:n1,1:n2,1:n3,ip) = scr4(1:n1,1:n2,1:n3) +#else read(31) #endif read(31) (((scr4(i,j,k),i=1,n1),j=1,n2),k=1,n3) @@ -280,6 +320,21 @@ subroutine read_masclet_fields(iter, flag_filter, nx, ny, nz, read(31) ! pressure read(31) ! pot read(31) ! opot +#if weight_filter == 2 + read(31) (((scr4(i,j,k),i=1,nx),j=1,ny),k=1,nz) +!$omp parallel do shared(emissivity1,scr4,n1,n2,n3), +!$omp+ private(ix,jy,kz), default(none) + do kz=0,n1+1 + do jy=0,n2+1 + do ix=0,n3+1 + emissivity1(ix,jy,kz,ip) = emissivity1(ix,jy,kz,ip)* + & emissivity1(ix,jy,kz,ip)*sqrt(scr4(ix,jy,kz) + end do + end do + end do +#else + read(31) ! temp +#endif read(31) ! temp read(31) ! metal read(31) ! cr0 @@ -344,6 +399,10 @@ subroutine read_masclet_fields(iter, flag_filter, nx, ny, nz, call p_minmax_ir(dens0,dens1,1,0,nx,ny,nz,nl,patchnx,patchny, & patchnz,npatch,0,basx,basy) write(*,*) 'dens min,max',basx,basy +#elif weight_filter == 2 + call p_minmax_ir(emis0,emis1,1,0,nx,ny,nz,nl,patchnx,patchny, + & patchnz,npatch,0,basx,basy) + write(*,*) 'emis min,max',basx,basy #endif call p_minmax_ir(u2,u12,1,1,nx,ny,nz,nl,patchnx,patchny,patchnz, & npatch,0,basx,basy) @@ -372,6 +431,10 @@ subroutine read_masclet_fields(iter, flag_filter, nx, ny, nz, call p_minmax_ir(dens0,dens1,1,0,nx,ny,nz,nl,patchnx,patchny, & patchnz,npatch,ir,basx,basy) write(*,*) 'dens min,max',basx,basy +#elif weight_filter == 2 + call p_minmax_ir(emis0,emis1,1,0,nx,ny,nz,nl,patchnx,patchny, + & patchnz,npatch,ir,basx,basy) + write(*,*) 'emis min,max',basx,basy #endif call p_minmax_ir(u2,u12,1,1,nx,ny,nz,nl,patchnx,patchny,patchnz, & npatch,ir,basx,basy) @@ -450,4 +513,4 @@ subroutine read_masclet_fix_grid_level(nl, nl_vortex, npatch, end if return - end \ No newline at end of file + end From 59988dc13fc28f60a6ad56854266bf05a81658f3 Mon Sep 17 00:00:00 2001 From: Frederick Groth Date: Wed, 27 Aug 2025 16:11:29 +0200 Subject: [PATCH 06/16] assign emissivity to grid. So far, this uses the same weighting as all other quantities (even though it probably should be volume weighted) --- src/particles.f | 106 +++++++++++++++++++++++++++++++++++++++++------- 1 file changed, 91 insertions(+), 15 deletions(-) diff --git a/src/particles.f b/src/particles.f index 7f9f2f3..c0136d6 100644 --- a/src/particles.f +++ b/src/particles.f @@ -1225,6 +1225,9 @@ SUBROUTINE ERROR_PARTICLES(NX,NY,NZ,NL,NPATCH,PATCHNX,PATCHNY, * Local variables INTEGER IPATCH,IX,JY,KZ,LOW1,LOW2,IR,IP,II REAL BASX,BASY,BASZ,BASXX,BASYY,BASZZ,BAS,BAS2,DXPA,DYPA,DZPA +#if weight_filter == 2 + REAL BASEMISS +#endif REAL UBAS(3,3,3),RXBAS(3),RYBAS(3),RZBAS(3),AAA,BBB,CCC REAL PATCHLEV(NPALEV) @@ -1295,7 +1298,6 @@ SUBROUTINE ERROR_PARTICLES(NX,NY,NZ,NL,NPATCH,PATCHNX,PATCHNY, UBAS(1:3,1:3,1:3)=U14(IX-1:IX+1,JY-1:JY+1,KZ-1:KZ+1,IPATCH) CALL LININT52D_NEW_REAL(AAA,BBB,CCC,RXBAS,RYBAS,RZBAS,UBAS, & BASZ) - ELSE IF (IX.GT.1.AND.IX.LT.NX.AND. & JY.GT.1.AND.JY.LT.NY.AND. @@ -2315,12 +2317,20 @@ SUBROUTINE INTERPOLATE_VELOCITIES(NX,NY,NZ,NL,NPATCH,PARE, COMMON /DENSI/ L0,L1 #endif +#if weight_filter == 2 + REAL EMISS0(0:NMAX+1,0:NMAY+1,0:NMAZ+1) + REAL EMISS1(NAMRX,NAMRY,NAMRZ,NPALEV) + COMMON /EMISS/ EMISS0,EMISS1 +#else + REAL EMISS0, EMISS1 +#endif + INTEGER IX,JY,KZ,IR,I,J,K,IPATCH,LOW1,LOW2,CONTA,KNEIGHBOURS INTEGER N1,N2,N3,II,JJ,KK,JPATCH,I1,I2,J1,J2,K1,K2,STEP INTEGER NPART_TOT,II1,II2,JJ1,JJ2,KK1,KK2,IIP1,JJP1,KKP1 REAL DXPA,DYPA,DZPA,BASX,BASY,BASZ,BAS,RBAS,BASXX,BASYY,BASZZ REAL MEDIOLADO0,PI,MINX,MAXX,MINY,MAXY,MINZ,MAXZ,FRAC_INT,H_KERN - REAL*8 BAS8,BAS8X,BAS8Y,BAS8Z,BAS8M,BASMASS + REAL*8 BAS8,BAS8X,BAS8Y,BAS8Z,BAS8M,BASMASS,BAS8EMISS REAL,ALLOCATABLE::DIST(:) INTEGER,ALLOCATABLE::NEIGH(:),LB(:) REAL FUIN,U(2,2,2),UW(2,2,2) @@ -2484,11 +2494,12 @@ SUBROUTINE INTERPOLATE_VELOCITIES(NX,NY,NZ,NL,NPATCH,PARE, !$OMP PARALLEL DO SHARED(NX,NY,NZ,STEP,I1,I2,J1,J2,K1,K2,RADX,RADY, !$OMP+ RADZ,U2DM,U3DM,U4DM,L0,U2,U3,U4,MASAP,VOL, -!$OMP+ EMISSIVITY, +!$OMP+ EMISSIVITY,EMISS0, !$OMP+ KNEIGHBOURS,DX,VISC0,ABVC,TREE,XTREE, !$OMP+ YTREE,ZTREE,PI,FLAG_MASS), !$OMP+ PRIVATE(IX,JY,KZ,DIST,NEIGH,CONTA,H_KERN,BAS8, -!$OMP+ BAS8X,BAS8Y,BAS8Z,BAS8M,I,BASMASS,DA,SEARCH), +!$OMP+ BAS8X,BAS8Y,BAS8Z,BAS8M,I,BASMASS, +!$OMP+ BAS8EMISS,DA,SEARCH), !$OMP+ SCHEDULE(DYNAMIC)!, DEFAULT(NONE) DO KZ=1,NZ+1,STEP DO JY=1,NY+1,STEP @@ -2545,6 +2556,7 @@ SUBROUTINE INTERPOLATE_VELOCITIES(NX,NY,NZ,NL,NPATCH,PARE, BAS8Z=0.D0 BAS8M=0.D0 BASMASS=0.D0 + BAS8EMISS=0.D0 DO I=1,CONTA BAS8=BAS8+DIST(I) BAS8X=BAS8X+DIST(I)*U2DM(NEIGH(I)) @@ -2552,6 +2564,11 @@ SUBROUTINE INTERPOLATE_VELOCITIES(NX,NY,NZ,NL,NPATCH,PARE, BAS8Z=BAS8Z+DIST(I)*U4DM(NEIGH(I)) #if use_filter == 1 BAS8M=BAS8M+DIST(I)*ABVC(NEIGH(I)) +#endif +#if filter_weight == 2 + ! emissivity is always volume weighted + ! todo! + BAS8EMISS=BAS8EMISS+DIST(I)*EMISSIVITY(NEIGH(I)) #endif !BAS8M=MAX(BAS8M,ABVC(NEIGH(I))) BASMASS=BASMASS+MASAP(NEIGH(I)) @@ -2563,7 +2580,10 @@ SUBROUTINE INTERPOLATE_VELOCITIES(NX,NY,NZ,NL,NPATCH,PARE, VISC0(IX,JY,KZ)=BAS8M/BAS8 #endif !VISC0(IX,JY,KZ)=BAS8M - +#if filter_weight == 2 + ! todo: different weighting + EMISS0(IX,IY,IZ)=BAS8EMISS/BAS8 +#endif IF (FLAG_MASS.EQ.1) L0(IX,JY,KZ)=BASMASS/(4*PI/3)/H_KERN**3 @@ -2575,7 +2595,8 @@ SUBROUTINE INTERPOLATE_VELOCITIES(NX,NY,NZ,NL,NPATCH,PARE, ! Fill the blanks by interpolation !$OMP PARALLEL DO SHARED(NX,NY,NZ,STEP,I1,I2,J1,J2,K1,K2,RADX,RADY, -!$OMP+ RADZ,U2,U3,U4,L0,DX,DY,DZ,VISC0), +!$OMP+ RADZ,U2,U3,U4,L0,DX,DY,DZ,VISC0, +!$OMP+ EMISS0), !$OMP+ PRIVATE(IX,JY,KZ,II,JJ,KK,IIP1,JJP1,KKP1,BASX,BASY, !$OMP+ BASZ), !$OMP+ DEFAULT(NONE), SCHEDULE(DYNAMIC) @@ -2644,6 +2665,17 @@ SUBROUTINE INTERPOLATE_VELOCITIES(NX,NY,NZ,NL,NPATCH,PARE, & VISC0(IIP1,JJP1,KKP1)* BASX * BASY * BASZ #endif +#if filter_weight == 2 + EMISS0(IX,JY,KZ)=EMISS0(II,JJ,KK) *(1-BASX)*(1-BASY)*(1-BASZ) + + & EMISS0(IIP1,JJ,KK)* BASX *(1-BASY)*(1-BASZ) + + & EMISS0(II,JJP1,KK)*(1-BASX)* BASY *(1-BASZ) + + & EMISS0(IIP1,JJP1,KK)* BASX * BASY *(1-BASZ) + + & EMISS0(II,JJ,KKP1)*(1-BASX)*(1-BASY)* BASZ + + & EMISS0(IIP1,JJ,KKP1)* BASX *(1-BASY)* BASZ + + & EMISS0(II,JJP1,KKP1)*(1-BASX)* BASY * BASZ + + & EMISS0(IIP1,JJP1,KKP1)* BASX * BASY * BASZ +#endif + END DO END DO END DO @@ -2657,9 +2689,10 @@ SUBROUTINE INTERPOLATE_VELOCITIES(NX,NY,NZ,NL,NPATCH,PARE, !$OMP+ JJ2,KK1,KK2,RADX,RADY,RADZ,U2DM,U3DM, !$OMP+ U4DM,L0,U2,U3,U4,KNEIGHBOURS,DX,VISC0,ABVC, !$OMP+ TREE,XTREE,YTREE,ZTREE,FLAG_MASS, -!$OMP+ MASAP,VOL,EMISSIVITY,PI), +!$OMP+ MASAP,VOL,EMISSIVITY,EMISS0,PI), !$OMP+ PRIVATE(IX,JY,KZ,DIST,NEIGH,CONTA,H_KERN,BAS8, -!$OMP+ BAS8X,BAS8Y,BAS8Z,BAS8M,I,BASMASS,DA,SEARCH), +!$OMP+ BAS8X,BAS8Y,BAS8Z,BAS8M,I,BASMASS, +!$OMP+ BAS8EMISS,DA,SEARCH), !$OMP+ SCHEDULE(DYNAMIC)!, DEFAULT(NONE) DO KZ=K1,K2+1,STEP DO JY=J1,J2+1,STEP @@ -2716,6 +2749,7 @@ SUBROUTINE INTERPOLATE_VELOCITIES(NX,NY,NZ,NL,NPATCH,PARE, BAS8Z=0.D0 BAS8M=0.D0 BASMASS=0.D0 + BAS8EMISS=0.D0 DO I=1,CONTA BAS8=BAS8+DIST(I) BAS8X=BAS8X+DIST(I)*U2DM(NEIGH(I)) @@ -2726,6 +2760,9 @@ SUBROUTINE INTERPOLATE_VELOCITIES(NX,NY,NZ,NL,NPATCH,PARE, #endif !BAS8M=MAX(BAS8M,ABVC(NEIGH(I))) BASMASS=BASMASS+MASAP(NEIGH(I)) +#if weight_filter == 2 + BAS8EMISS=BAS8EMISS+DIST(I)*EMISSIVITY(NEIGH(I)) +#endif END DO U2(IX,JY,KZ)=BAS8X/BAS8 @@ -2735,7 +2772,10 @@ SUBROUTINE INTERPOLATE_VELOCITIES(NX,NY,NZ,NL,NPATCH,PARE, VISC0(IX,JY,KZ)=BAS8M/BAS8 #endif !VISC0(IX,JY,KZ)=BAS8M - +#if weight_filter == 2 + EMISS0(IX,JY,KZ)=BAS8EMISS/BAS8 +#endif + IF (FLAG_MASS.EQ.1) L0(IX,JY,KZ)=BASMASS/(4*PI/3)/H_KERN**3 DEALLOCATE(DIST, NEIGH) @@ -2747,7 +2787,7 @@ SUBROUTINE INTERPOLATE_VELOCITIES(NX,NY,NZ,NL,NPATCH,PARE, !$OMP PARALLEL DO SHARED(NX,NY,NZ,STEP,II1,II2,JJ1,JJ2,KK1,KK2,RADX, !$OMP+ RADY,RADZ,U2,U3,U4,L0,DX,DY,DZ,I1,I2,J1,J2,K1, -!$OMP+ K2,VISC0), +!$OMP+ K2,VISC0,EMISS0), !$OMP+ PRIVATE(IX,JY,KZ,II,JJ,KK,IIP1,JJP1,KKP1,BASX,BASY, !$OMP+ BASZ), !$OMP+ DEFAULT(NONE), SCHEDULE(DYNAMIC) @@ -2816,6 +2856,17 @@ SUBROUTINE INTERPOLATE_VELOCITIES(NX,NY,NZ,NL,NPATCH,PARE, & VISC0(IIP1,JJP1,KKP1)* BASX * BASY * BASZ #endif +#if weight_filter == 2 + EMISS0(IX,JY,KZ)=EMISS0(II,JJ,KK) *(1-BASX)*(1-BASY)*(1-BASZ) + + & EMISS0(IIP1,JJ,KK)* BASX *(1-BASY)*(1-BASZ) + + & EMISS0(II,JJP1,KK)*(1-BASX)* BASY *(1-BASZ) + + & EMISS0(IIP1,JJP1,KK)* BASX * BASY *(1-BASZ) + + & EMISS0(II,JJ,KKP1)*(1-BASX)*(1-BASY)* BASZ + + & EMISS0(IIP1,JJ,KKP1)* BASX *(1-BASY)* BASZ + + & EMISS0(II,JJP1,KKP1)*(1-BASX)* BASY * BASZ + + & EMISS0(IIP1,JJP1,KKP1)* BASX * BASY * BASZ +#endif + END DO END DO END DO @@ -2829,9 +2880,10 @@ SUBROUTINE INTERPOLATE_VELOCITIES(NX,NY,NZ,NL,NPATCH,PARE, !$OMP+ RADZ,U2DM,U3DM,U4DM,L0,U2,U3,U4, !$OMP+ KNEIGHBOURS,DX,CR0AMR,VISC0,ABVC, !$OMP+ TREE,XTREE,YTREE,ZTREE,FLAG_MASS,PI,MASAP,VOL, -!$OMP+ EMISSIVITY), +!$OMP+ EMISSIVITY,EMISS0), !$OMP+ PRIVATE(IX,JY,KZ,DIST,NEIGH,CONTA,H_KERN,BAS8, -!$OMP+ BAS8X,BAS8Y,BAS8Z,BAS8M,I,BASMASS,DA,SEARCH), +!$OMP+ BAS8X,BAS8Y,BAS8Z,BAS8M,I,BASMASS, +!$OMP+ BAS8EMISS,DA,SEARCH), !$OMP+ SCHEDULE(DYNAMIC)!, DEFAULT(NONE) DO KZ=KK1,KK2 DO JY=JJ1,JJ2 @@ -2883,6 +2935,7 @@ SUBROUTINE INTERPOLATE_VELOCITIES(NX,NY,NZ,NL,NPATCH,PARE, BAS8Z=0.D0 BAS8M=0.D0 BASMASS=0.D0 + BAS8EMISS=0.D0 DO I=1,CONTA BAS8=BAS8+DIST(I) BAS8X=BAS8X+DIST(I)*U2DM(NEIGH(I)) @@ -2893,6 +2946,9 @@ SUBROUTINE INTERPOLATE_VELOCITIES(NX,NY,NZ,NL,NPATCH,PARE, #endif !BAS8M=MAX(BAS8M,ABVC(NEIGH(I))) BASMASS=BASMASS+MASAP(NEIGH(I)) +#if weight_filter == 2 + BAS8EMISS=BAS8EMISS+DIST(I)*EMISSIVITY(NEIGH(I)) +#endif END DO U2(IX,JY,KZ)=BAS8X/BAS8 @@ -2902,7 +2958,10 @@ SUBROUTINE INTERPOLATE_VELOCITIES(NX,NY,NZ,NL,NPATCH,PARE, VISC0(IX,JY,KZ)=BAS8M/BAS8 #endif !VISC0(IX,JY,KZ)=BAS8M - +#if weight_filter == 2 + EMISS0(IX,JY,KZ)=BAS8EMISS/BAS8 +#endif + IF (FLAG_MASS.EQ.1) L0(IX,JY,KZ)=BASMASS/(4*PI/3)/H_KERN**3 DEALLOCATE(DIST, NEIGH) @@ -2925,10 +2984,10 @@ SUBROUTINE INTERPOLATE_VELOCITIES(NX,NY,NZ,NL,NPATCH,PARE, !$OMP+ RX,RY,RZ,U2DM,U3DM,U4DM,L1,U12,U13,U14, !$OMP+ KNEIGHBOURS,DXPA,DX,VISC1,ABVC,SOLAP, !$OMP+ TREE,XTREE,YTREE,ZTREE,FLAG_MASS,PI,MASAP,VOL, -!$OMP+ EMISSIVITY), +!$OMP+ EMISSIVITY,EMISS1), !$OMP+ PRIVATE(IPATCH,N1,N2,N3,IX,JY,KZ,DIST,NEIGH,DA, !$OMP+ CONTA,H_KERN,BAS8,BAS8X,BAS8Y,BAS8Z,BAS8M,I, -!$OMP+ BASMASS,SEARCH,DO_CELL), +!$OMP+ BASMASS,BAS8EMISS,SEARCH,DO_CELL), !$OMP+ SCHEDULE(DYNAMIC,1)!, DEFAULT(NONE) DO IPATCH=LOW1,LOW2 !write(*,*) ir,ipatch @@ -3000,6 +3059,7 @@ SUBROUTINE INTERPOLATE_VELOCITIES(NX,NY,NZ,NL,NPATCH,PARE, BAS8Z=0.D0 BAS8M=0.D0 BASMASS=0.D0 + BAS8EMISS=0.D0 DO I=1,CONTA BAS8=BAS8+DIST(I) BAS8X=BAS8X+DIST(I)*U2DM(NEIGH(I)) @@ -3010,6 +3070,9 @@ SUBROUTINE INTERPOLATE_VELOCITIES(NX,NY,NZ,NL,NPATCH,PARE, #endif !BAS8M=MAX(BAS8M,ABVC(NEIGH(I))) BASMASS=BASMASS+MASAP(NEIGH(I)) +#if weight_filter == 2 + BAS8EMISS=BAS8EMISS+DIST(I)*EMISSIVITY(NEIGH(I)) +#endif END DO DEALLOCATE(DIST, NEIGH) @@ -3017,6 +3080,9 @@ SUBROUTINE INTERPOLATE_VELOCITIES(NX,NY,NZ,NL,NPATCH,PARE, U12(IX,JY,KZ,IPATCH)=BAS8X/BAS8 U13(IX,JY,KZ,IPATCH)=BAS8Y/BAS8 U14(IX,JY,KZ,IPATCH)=BAS8Z/BAS8 +#if weight_filter == 2 + EMISS1(IX,JY,KZ,IPATCH)=BAS8EMISS/BAS8 +#endif if (ix.eq.0.or.ix.eq.n1+1.or. & jy.eq.0.or.jy.eq.n2+1.or. @@ -3155,6 +3221,11 @@ SUBROUTINE INTERPOLATE_VELOCITIES(NX,NY,NZ,NL,NPATCH,PARE, WRITE(*,*) 'Mach min,max',BASX,BASY END IF #endif +#if weight_filter == 2 + CALL P_MINMAX_IR(EMISS0,EMISS1,1,1,NX,NY,NZ,NL,PATCHNX,PATCHNY, + & PATCHNZ,NPATCH,0,BASX,BASY) + write(*,*) 'emissivity min,max',BASX,BASY +#endif DO IR=1,NL @@ -3181,6 +3252,11 @@ SUBROUTINE INTERPOLATE_VELOCITIES(NX,NY,NZ,NL,NPATCH,PARE, ELSE WRITE(*,*) 'Mach min,max',BASX,BASY END IF +#endif +#if weight_filter == 2 + CALL P_MINMAX_IR(EMISS0,EMISS1,1,1,NX,NY,NZ,NL,PATCHNX,PATCHNY, + & PATCHNZ,NPATCH,IR,BASX,BASY) + write(*,*) 'emissivity min,max',BASX,BASY #endif END DO From 588a6a940d9dfd81e4a61cf975432258625309e9 Mon Sep 17 00:00:00 2001 From: Frederick Groth Date: Fri, 29 Aug 2025 13:30:43 +0200 Subject: [PATCH 07/16] ensure the emissivity used for weight_filter=2 is always assigned to the grid volume weighted --- src/particle_data.f | 2 +- src/particles.f | 112 ++++++++++++++++++++++++++----- src/reader.f | 14 ++-- src/readers/arepo_hdf5.f | 5 +- src/readers/gadget_unformatted.f | 5 +- 5 files changed, 109 insertions(+), 29 deletions(-) diff --git a/src/particle_data.f b/src/particle_data.f index 9dc4cdc..c52bd07 100644 --- a/src/particle_data.f +++ b/src/particle_data.f @@ -15,7 +15,7 @@ MODULE PARTICLE_DATA REAL ABVC #endif -#if weight_scheme == 2 +#if weight_scheme == 2 || weight_filter == 2 REAL,ALLOCATABLE::VOL(:) #else ! Dummy variable diff --git a/src/particles.f b/src/particles.f index c0136d6..618e775 100644 --- a/src/particles.f +++ b/src/particles.f @@ -1225,9 +1225,6 @@ SUBROUTINE ERROR_PARTICLES(NX,NY,NZ,NL,NPATCH,PATCHNX,PATCHNY, * Local variables INTEGER IPATCH,IX,JY,KZ,LOW1,LOW2,IR,IP,II REAL BASX,BASY,BASZ,BASXX,BASYY,BASZZ,BAS,BAS2,DXPA,DYPA,DZPA -#if weight_filter == 2 - REAL BASEMISS -#endif REAL UBAS(3,3,3),RXBAS(3),RYBAS(3),RZBAS(3),AAA,BBB,CCC REAL PATCHLEV(NPALEV) @@ -2330,8 +2327,13 @@ SUBROUTINE INTERPOLATE_VELOCITIES(NX,NY,NZ,NL,NPATCH,PARE, INTEGER NPART_TOT,II1,II2,JJ1,JJ2,KK1,KK2,IIP1,JJP1,KKP1 REAL DXPA,DYPA,DZPA,BASX,BASY,BASZ,BAS,RBAS,BASXX,BASYY,BASZZ REAL MEDIOLADO0,PI,MINX,MAXX,MINY,MAXY,MINZ,MAXZ,FRAC_INT,H_KERN - REAL*8 BAS8,BAS8X,BAS8Y,BAS8Z,BAS8M,BASMASS,BAS8EMISS + REAL*8 BAS8,BAS8X,BAS8Y,BAS8Z,BAS8M,BASMASS,BAS8EMISS,BAS8_VOL REAL,ALLOCATABLE::DIST(:) +#if weight_filter == 2 + REAL,ALLOCATABLE::VOL_DIST(:) +#else + REAL VOL_DIST ! dummy variable +#endif INTEGER,ALLOCATABLE::NEIGH(:),LB(:) REAL FUIN,U(2,2,2),UW(2,2,2) logical do_cell @@ -2499,7 +2501,7 @@ SUBROUTINE INTERPOLATE_VELOCITIES(NX,NY,NZ,NL,NPATCH,PARE, !$OMP+ YTREE,ZTREE,PI,FLAG_MASS), !$OMP+ PRIVATE(IX,JY,KZ,DIST,NEIGH,CONTA,H_KERN,BAS8, !$OMP+ BAS8X,BAS8Y,BAS8Z,BAS8M,I,BASMASS, -!$OMP+ BAS8EMISS,DA,SEARCH), +!$OMP+ BAS8EMISS,BAS8_VOL,DA,SEARCH,VOL_DIST), !$OMP+ SCHEDULE(DYNAMIC)!, DEFAULT(NONE) DO KZ=1,NZ+1,STEP DO JY=1,NY+1,STEP @@ -2512,6 +2514,9 @@ SUBROUTINE INTERPOLATE_VELOCITIES(NX,NY,NZ,NL,NPATCH,PARE, & KZ.LT.INT(K2/STEP)*STEP) CYCLE ALLOCATE(DIST(KNEIGHBOURS), NEIGH(KNEIGHBOURS)) +#if weight_filter == 2 + ALLOCATE(VOL_DIST(KNEIGHBOURS)) +#endif DA = SEARCH%KNEAREST(TREE, XTREE, YTREE, ZTREE, & XQUERY=DBLE(RADX(IX)), YQUERY=DBLE(RADY(JY)), & ZQUERY=DBLE(RADZ(KZ)), K=KNEIGHBOURS) @@ -2521,12 +2526,18 @@ SUBROUTINE INTERPOLATE_VELOCITIES(NX,NY,NZ,NL,NPATCH,PARE, IF (DIST(KNEIGHBOURS).GT.DX) THEN CONTA=KNEIGHBOURS ELSE +#if weight_filter == 2 + DEALLOCATE(VOL_DIST) +#endif DEALLOCATE(DIST,NEIGH) DA = SEARCH%KNEAREST(TREE, XTREE, YTREE, ZTREE, & XQUERY=DBLE(RADX(IX)), YQUERY=DBLE(RADY(JY)), & ZQUERY=DBLE(RADZ(KZ)), RADIUS=DBLE(DX)) CONTA = DA%SIZE() ALLOCATE(DIST(CONTA), NEIGH(CONTA)) +#if weight_filter == 2 + ALLOCATE(VOL_DIST(CONTA)) +#endif DIST = DA%V%VALUES NEIGH = DA%I%VALUES END IF @@ -2536,6 +2547,11 @@ SUBROUTINE INTERPOLATE_VELOCITIES(NX,NY,NZ,NL,NPATCH,PARE, L0(IX,JY,KZ)=H_KERN CALL KERNEL_FUNC(CONTA,CONTA,H_KERN,DIST) +#if weight_filter == 2 + DO I=1,CONTA + VOL_DIST(I)=DIST(I)*VOL(NEIGH(I)) + END DO +#endif #if weight_scheme == 1 DO I=1,CONTA DIST(I)=DIST(I)*MASAP(NEIGH(I)) @@ -2557,6 +2573,7 @@ SUBROUTINE INTERPOLATE_VELOCITIES(NX,NY,NZ,NL,NPATCH,PARE, BAS8M=0.D0 BASMASS=0.D0 BAS8EMISS=0.D0 + BAS8_VOL=0.D0 DO I=1,CONTA BAS8=BAS8+DIST(I) BAS8X=BAS8X+DIST(I)*U2DM(NEIGH(I)) @@ -2567,8 +2584,8 @@ SUBROUTINE INTERPOLATE_VELOCITIES(NX,NY,NZ,NL,NPATCH,PARE, #endif #if filter_weight == 2 ! emissivity is always volume weighted - ! todo! - BAS8EMISS=BAS8EMISS+DIST(I)*EMISSIVITY(NEIGH(I)) + BAS8_VOL=BAS8_VOL+VOL_DIST(I) + BAS8EMISS=BAS8EMISS+VOL_DIST(I)*EMISSIVITY(NEIGH(I)) #endif !BAS8M=MAX(BAS8M,ABVC(NEIGH(I))) BASMASS=BASMASS+MASAP(NEIGH(I)) @@ -2582,11 +2599,14 @@ SUBROUTINE INTERPOLATE_VELOCITIES(NX,NY,NZ,NL,NPATCH,PARE, !VISC0(IX,JY,KZ)=BAS8M #if filter_weight == 2 ! todo: different weighting - EMISS0(IX,IY,IZ)=BAS8EMISS/BAS8 + EMISS0(IX,IY,IZ)=BAS8EMISS/BAS8_VOL #endif IF (FLAG_MASS.EQ.1) L0(IX,JY,KZ)=BASMASS/(4*PI/3)/H_KERN**3 +#if weight_filter == 2 + DEALLOCATE(VOL_DIST) +#endif DEALLOCATE(DIST, NEIGH) END DO END DO @@ -2692,7 +2712,7 @@ SUBROUTINE INTERPOLATE_VELOCITIES(NX,NY,NZ,NL,NPATCH,PARE, !$OMP+ MASAP,VOL,EMISSIVITY,EMISS0,PI), !$OMP+ PRIVATE(IX,JY,KZ,DIST,NEIGH,CONTA,H_KERN,BAS8, !$OMP+ BAS8X,BAS8Y,BAS8Z,BAS8M,I,BASMASS, -!$OMP+ BAS8EMISS,DA,SEARCH), +!$OMP+ BAS8EMISS,BAS8_VOL,DA,SEARCH,VOL_DIST), !$OMP+ SCHEDULE(DYNAMIC)!, DEFAULT(NONE) DO KZ=K1,K2+1,STEP DO JY=J1,J2+1,STEP @@ -2705,6 +2725,9 @@ SUBROUTINE INTERPOLATE_VELOCITIES(NX,NY,NZ,NL,NPATCH,PARE, & KZ.LT.INT(KK2/STEP)*STEP) CYCLE ALLOCATE(DIST(KNEIGHBOURS), NEIGH(KNEIGHBOURS)) +#if weight_filter == 2 + ALLOCATE(VOL_DIST(KNEIGHBOURS)) +#endif DA = SEARCH%KNEAREST(TREE, XTREE, YTREE, ZTREE, & XQUERY=DBLE(RADX(IX)), YQUERY=DBLE(RADY(JY)), & ZQUERY=DBLE(RADZ(KZ)), K=KNEIGHBOURS) @@ -2714,12 +2737,18 @@ SUBROUTINE INTERPOLATE_VELOCITIES(NX,NY,NZ,NL,NPATCH,PARE, IF (DIST(KNEIGHBOURS).GT.DX) THEN CONTA=KNEIGHBOURS ELSE +#if weight_filter == 2 + DEALLOCATE(VOL_DIST) +#endif DEALLOCATE(DIST,NEIGH) DA = SEARCH%KNEAREST(TREE, XTREE, YTREE, ZTREE, & XQUERY=DBLE(RADX(IX)), YQUERY=DBLE(RADY(JY)), & ZQUERY=DBLE(RADZ(KZ)), RADIUS=DBLE(DX)) CONTA = DA%SIZE() ALLOCATE(DIST(CONTA), NEIGH(CONTA)) +#if weight_filter == 2 + ALLOCATE(VOL_DIST(CONTA)) +#endif DIST = DA%V%VALUES NEIGH = DA%I%VALUES END IF @@ -2729,6 +2758,11 @@ SUBROUTINE INTERPOLATE_VELOCITIES(NX,NY,NZ,NL,NPATCH,PARE, L0(IX,JY,KZ)=H_KERN CALL KERNEL_FUNC(CONTA,CONTA,H_KERN,DIST) +#if weight_filter == 2 + DO I=1,CONTA + VOL_DIST(I)=DIST(I)*VOL(NEIGH(I)) + END DO +#endif #if weight_scheme == 1 DO I=1,CONTA DIST(I)=DIST(I)*MASAP(NEIGH(I)) @@ -2750,6 +2784,7 @@ SUBROUTINE INTERPOLATE_VELOCITIES(NX,NY,NZ,NL,NPATCH,PARE, BAS8M=0.D0 BASMASS=0.D0 BAS8EMISS=0.D0 + BAS8_VOL=0.D0 DO I=1,CONTA BAS8=BAS8+DIST(I) BAS8X=BAS8X+DIST(I)*U2DM(NEIGH(I)) @@ -2761,7 +2796,8 @@ SUBROUTINE INTERPOLATE_VELOCITIES(NX,NY,NZ,NL,NPATCH,PARE, !BAS8M=MAX(BAS8M,ABVC(NEIGH(I))) BASMASS=BASMASS+MASAP(NEIGH(I)) #if weight_filter == 2 - BAS8EMISS=BAS8EMISS+DIST(I)*EMISSIVITY(NEIGH(I)) + BAS8_VOL=BAS8_VOL+VOL_DIST(I) + BAS8EMISS=BAS8EMISS+VOL_DIST(I)*EMISSIVITY(NEIGH(I)) #endif END DO @@ -2773,11 +2809,14 @@ SUBROUTINE INTERPOLATE_VELOCITIES(NX,NY,NZ,NL,NPATCH,PARE, #endif !VISC0(IX,JY,KZ)=BAS8M #if weight_filter == 2 - EMISS0(IX,JY,KZ)=BAS8EMISS/BAS8 + EMISS0(IX,JY,KZ)=BAS8EMISS/BAS8_VOL #endif IF (FLAG_MASS.EQ.1) L0(IX,JY,KZ)=BASMASS/(4*PI/3)/H_KERN**3 +#if weight_filter == 2 + DEALLOCATE(VOL_DIST) +#endif DEALLOCATE(DIST, NEIGH) END DO END DO @@ -2883,7 +2922,7 @@ SUBROUTINE INTERPOLATE_VELOCITIES(NX,NY,NZ,NL,NPATCH,PARE, !$OMP+ EMISSIVITY,EMISS0), !$OMP+ PRIVATE(IX,JY,KZ,DIST,NEIGH,CONTA,H_KERN,BAS8, !$OMP+ BAS8X,BAS8Y,BAS8Z,BAS8M,I,BASMASS, -!$OMP+ BAS8EMISS,DA,SEARCH), +!$OMP+ BAS8EMISS,BAS8_VOL,DA,SEARCH,VOL_DIST), !$OMP+ SCHEDULE(DYNAMIC)!, DEFAULT(NONE) DO KZ=KK1,KK2 DO JY=JJ1,JJ2 @@ -2891,6 +2930,9 @@ SUBROUTINE INTERPOLATE_VELOCITIES(NX,NY,NZ,NL,NPATCH,PARE, DO IX=II1,II2 IF (CR0AMR(IX,JY,KZ).EQ.1) THEN ALLOCATE(DIST(KNEIGHBOURS), NEIGH(KNEIGHBOURS)) +#if weight_filter == 2 + ALLOCATE(VOL_DIST(KNEIGHBOURS), NEIGH(KNEIGHBOURS)) +#endif DA = SEARCH%KNEAREST(TREE, XTREE, YTREE, ZTREE, & XQUERY=DBLE(RADX(IX)), YQUERY=DBLE(RADY(JY)), & ZQUERY=DBLE(RADZ(KZ)), K=KNEIGHBOURS) @@ -2900,12 +2942,18 @@ SUBROUTINE INTERPOLATE_VELOCITIES(NX,NY,NZ,NL,NPATCH,PARE, IF (DIST(KNEIGHBOURS).GT.DX) THEN CONTA=KNEIGHBOURS ELSE +#if weight_filter == 2 + DEALLOCATE(VOL_DIST) +#endif DEALLOCATE(DIST,NEIGH) DA = SEARCH%KNEAREST(TREE, XTREE, YTREE, ZTREE, & XQUERY=DBLE(RADX(IX)), YQUERY=DBLE(RADY(JY)), & ZQUERY=DBLE(RADZ(KZ)), RADIUS=DBLE(DX)) CONTA = DA%SIZE() ALLOCATE(DISt(CONTA), NEIGH(CONTA)) +#if weight_filter == 2 + ALLOCATE(VOL_DIST(CONTA)) +#endif DIST = DA%V%VALUES NEIGH = DA%I%VALUES END IF @@ -2915,6 +2963,11 @@ SUBROUTINE INTERPOLATE_VELOCITIES(NX,NY,NZ,NL,NPATCH,PARE, L0(IX,JY,KZ)=H_KERN CALL KERNEL_FUNC(CONTA,CONTA,H_KERN,DIST) +#if weight_filter == 2 + DO I=1,CONTA + VOL_DIST(I)=DIST(I)*VOL(NEIGH(I)) + END DO +#endif #if weight_scheme == 1 DO I=1,CONTA DIST(I)=DIST(I)*MASAP(NEIGH(I)) @@ -2936,6 +2989,7 @@ SUBROUTINE INTERPOLATE_VELOCITIES(NX,NY,NZ,NL,NPATCH,PARE, BAS8M=0.D0 BASMASS=0.D0 BAS8EMISS=0.D0 + BAS8_VOL=0.D0 DO I=1,CONTA BAS8=BAS8+DIST(I) BAS8X=BAS8X+DIST(I)*U2DM(NEIGH(I)) @@ -2947,7 +3001,8 @@ SUBROUTINE INTERPOLATE_VELOCITIES(NX,NY,NZ,NL,NPATCH,PARE, !BAS8M=MAX(BAS8M,ABVC(NEIGH(I))) BASMASS=BASMASS+MASAP(NEIGH(I)) #if weight_filter == 2 - BAS8EMISS=BAS8EMISS+DIST(I)*EMISSIVITY(NEIGH(I)) + BAS8_VOL=BAS8_VOL+VOL_DIST(I) + BAS8EMISS=BAS8EMISS+VOL_DIST(I)*EMISSIVITY(NEIGH(I)) #endif END DO @@ -2959,11 +3014,14 @@ SUBROUTINE INTERPOLATE_VELOCITIES(NX,NY,NZ,NL,NPATCH,PARE, #endif !VISC0(IX,JY,KZ)=BAS8M #if weight_filter == 2 - EMISS0(IX,JY,KZ)=BAS8EMISS/BAS8 + EMISS0(IX,JY,KZ)=BAS8EMISS/BAS8_VOL #endif IF (FLAG_MASS.EQ.1) L0(IX,JY,KZ)=BASMASS/(4*PI/3)/H_KERN**3 +#if weight_filter == 2 + DEALLOCATE(VOL_DIST) +#endif DEALLOCATE(DIST, NEIGH) END IF END DO @@ -2987,7 +3045,8 @@ SUBROUTINE INTERPOLATE_VELOCITIES(NX,NY,NZ,NL,NPATCH,PARE, !$OMP+ EMISSIVITY,EMISS1), !$OMP+ PRIVATE(IPATCH,N1,N2,N3,IX,JY,KZ,DIST,NEIGH,DA, !$OMP+ CONTA,H_KERN,BAS8,BAS8X,BAS8Y,BAS8Z,BAS8M,I, -!$OMP+ BASMASS,BAS8EMISS,SEARCH,DO_CELL), +!$OMP+ BASMASS,BAS8EMISS,BAS8_VOL,SEARCH,DO_CELL, +!$OMP+ VOL_DIST), !$OMP+ SCHEDULE(DYNAMIC,1)!, DEFAULT(NONE) DO IPATCH=LOW1,LOW2 !write(*,*) ir,ipatch @@ -3015,6 +3074,9 @@ SUBROUTINE INTERPOLATE_VELOCITIES(NX,NY,NZ,NL,NPATCH,PARE, !Q(3)=RZ(KZ,IPATCH) ALLOCATE(DIST(KNEIGHBOURS), NEIGH(KNEIGHBOURS)) +#if weight_filter == 2 + ALLOCATE(VOL_DIST(KNEIGHBOURS), NEIGH(KNEIGHBOURS)) +#endif DA = SEARCH%KNEAREST(TREE, XTREE, YTREE, ZTREE, & XQUERY=DBLE(RX(IX,IPATCH)), YQUERY=DBLE(RY(JY,IPATCH)), & ZQUERY=DBLE(RZ(KZ,IPATCH)), K=KNEIGHBOURS) @@ -3024,12 +3086,18 @@ SUBROUTINE INTERPOLATE_VELOCITIES(NX,NY,NZ,NL,NPATCH,PARE, IF (DIST(KNEIGHBOURS).GT.DXPA) THEN CONTA=KNEIGHBOURS ELSE +#if weight_filter == 2 + DEALLOCATE(VOL_DIST) +#endif DEALLOCATE(DIST,NEIGH) DA = SEARCH%KNEAREST(TREE, XTREE, YTREE, ZTREE, & XQUERY=DBLE(RX(IX,IPATCH)), YQUERY=DBLE(RY(JY,IPATCH)), & ZQUERY=DBLE(RZ(KZ,IPATCH)), RADIUS=DBLE(DXPA)) CONTA = DA%SIZE() ALLOCATE(DIST(CONTA), NEIGH(CONTA)) +#if weight_filter == 2 + ALLOCATE(VOL_DIST(CONTA)) +#endif DIST = DA%V%VALUES NEIGH = DA%I%VALUES END IF @@ -3039,6 +3107,11 @@ SUBROUTINE INTERPOLATE_VELOCITIES(NX,NY,NZ,NL,NPATCH,PARE, ! write(*,*) ipatch,ix,jy,kz,conta,h_kern,dxpa ! end if CALL KERNEL_FUNC(CONTA,CONTA,H_KERN,DIST) +#if weight_filter == 2 + DO I=1,CONTA + VOL_DIST(I)=DIST(I)*VOL(NEIGH(I)) + END DO +#endif #if weight_scheme == 1 DO I=1,CONTA DIST(I)=DIST(I)*MASAP(NEIGH(I)) @@ -3060,6 +3133,7 @@ SUBROUTINE INTERPOLATE_VELOCITIES(NX,NY,NZ,NL,NPATCH,PARE, BAS8M=0.D0 BASMASS=0.D0 BAS8EMISS=0.D0 + BAS8_VOL=0.D0 DO I=1,CONTA BAS8=BAS8+DIST(I) BAS8X=BAS8X+DIST(I)*U2DM(NEIGH(I)) @@ -3071,17 +3145,21 @@ SUBROUTINE INTERPOLATE_VELOCITIES(NX,NY,NZ,NL,NPATCH,PARE, !BAS8M=MAX(BAS8M,ABVC(NEIGH(I))) BASMASS=BASMASS+MASAP(NEIGH(I)) #if weight_filter == 2 - BAS8EMISS=BAS8EMISS+DIST(I)*EMISSIVITY(NEIGH(I)) + BAS8_VOL=BAS8_VOL+VOL_DIST(I) + BAS8EMISS=BAS8EMISS+VOL_DIST(I)*EMISSIVITY(NEIGH(I)) #endif END DO +#if weight_filter == 2 + DEALLOCATE(VOL_DIST) +#endif DEALLOCATE(DIST, NEIGH) U12(IX,JY,KZ,IPATCH)=BAS8X/BAS8 U13(IX,JY,KZ,IPATCH)=BAS8Y/BAS8 U14(IX,JY,KZ,IPATCH)=BAS8Z/BAS8 #if weight_filter == 2 - EMISS1(IX,JY,KZ,IPATCH)=BAS8EMISS/BAS8 + EMISS1(IX,JY,KZ,IPATCH)=BAS8EMISS/BAS8_VOL #endif if (ix.eq.0.or.ix.eq.n1+1.or. diff --git a/src/reader.f b/src/reader.f index 7fd86b6..9bbff71 100644 --- a/src/reader.f +++ b/src/reader.f @@ -153,9 +153,10 @@ subroutine read_particles(iter,files_per_snap,nx,ny,nz,t,zeta, deallocate(abvc) #endif -#if weight_scheme == 2 +#if weight_scheme == 2 || weight_filter == 2 deallocate(vol) -#elif weight_scheme == 3 +#endif +#if weight_scheme == 3 deallocate(emissivity) #endif end if @@ -168,9 +169,10 @@ subroutine read_particles(iter,files_per_snap,nx,ny,nz,t,zeta, allocate(abvc(parti)) #endif -#if weight_scheme == 2 +#if weight_scheme == 2 || weight_filter == 2 allocate(vol(parti)) -#elif weight_scheme == 3 +#endif +#if weight_scheme == 3 allocate(emissivity(parti)) #endif @@ -192,13 +194,11 @@ subroutine read_particles(iter,files_per_snap,nx,ny,nz,t,zeta, ************************************************************************ ************************************************************************ -#if weight_scheme == 2 +#if weight_scheme == 2 || weight_filter == 2 !$omp parallel do shared(vol, masap, parti), private(i), default(none) do i=1,parti vol(i)=masap(i)/vol(i) end do -#elif weight_scheme == 3 - ! this is claculated already during the reading process. #endif npart(0)=low2 !retrocompatibility with general reader diff --git a/src/readers/arepo_hdf5.f b/src/readers/arepo_hdf5.f index fc2cb76..4d319c3 100644 --- a/src/readers/arepo_hdf5.f +++ b/src/readers/arepo_hdf5.f @@ -241,7 +241,7 @@ subroutine read_arepo_hdf5(iter, files_per_snap, end if #endif -#if weight_scheme == 2 +#if weight_scheme == 2 || weight_filter == 2 allocate(scr4(numpart_thisfile(1))) write(*,*) 'reading density ...' call h5dopen_f(group_id, "density", attr_id, status) @@ -250,7 +250,8 @@ subroutine read_arepo_hdf5(iter, files_per_snap, vol(low1:low2)=scr4(1:numpart_thisfile(1)) call h5dclose_f(attr_id, status) deallocate(scr4) -#elif weight_scheme == 3 +#endif +#if weight_scheme == 3 allocate(scr4(numpart_thisfile(1))) write(*,*) 'reading density ...' call h5dopen_f(group_id, "density", attr_id, status) diff --git a/src/readers/gadget_unformatted.f b/src/readers/gadget_unformatted.f index c374f40..e79ed9d 100644 --- a/src/readers/gadget_unformatted.f +++ b/src/readers/gadget_unformatted.f @@ -141,14 +141,15 @@ subroutine read_gadget_unformatted(iter, files_per_snap, end if #endif -#if weight_scheme == 2 +#if weight_scheme == 2 || weight_filter == 2 allocate(scr4(sum(npart_gadget(1:6)))) write(*,*) 'reading density ...' call read_float(fil2,'RHO ',scr4,blocksize) write(*,*) ' found for ',(blocksize-8)/4,' particles' vol(low1:low2)=scr4(1:npart_gadget(1)) deallocate(scr4) -#elif weight_scheme == 3 +#endif +#if weight_scheme == 3 allocate(scr4(sum(npart_gadget(1:6)))) write(*,*) 'reading density ...' call read_float(fil2,'RHO ',scr4,blocksize) From f0fb3e5adeacc3e62195930b8e7464968098edc4 Mon Sep 17 00:00:00 2001 From: Frederick Groth Date: Fri, 29 Aug 2025 13:31:28 +0200 Subject: [PATCH 08/16] remove leftover comment --- src/particles.f | 1 - 1 file changed, 1 deletion(-) diff --git a/src/particles.f b/src/particles.f index 618e775..2bf8eb8 100644 --- a/src/particles.f +++ b/src/particles.f @@ -2598,7 +2598,6 @@ SUBROUTINE INTERPOLATE_VELOCITIES(NX,NY,NZ,NL,NPATCH,PARE, #endif !VISC0(IX,JY,KZ)=BAS8M #if filter_weight == 2 - ! todo: different weighting EMISS0(IX,IY,IZ)=BAS8EMISS/BAS8_VOL #endif From de9d7482a186ce61f1c7d15c83143e0ea4aff242 Mon Sep 17 00:00:00 2001 From: Frederick Groth Date: Fri, 29 Aug 2025 15:13:40 +0200 Subject: [PATCH 09/16] remove duplicate allocation --- src/particles.f | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/particles.f b/src/particles.f index 2bf8eb8..dfd9bd7 100644 --- a/src/particles.f +++ b/src/particles.f @@ -2930,7 +2930,7 @@ SUBROUTINE INTERPOLATE_VELOCITIES(NX,NY,NZ,NL,NPATCH,PARE, IF (CR0AMR(IX,JY,KZ).EQ.1) THEN ALLOCATE(DIST(KNEIGHBOURS), NEIGH(KNEIGHBOURS)) #if weight_filter == 2 - ALLOCATE(VOL_DIST(KNEIGHBOURS), NEIGH(KNEIGHBOURS)) + ALLOCATE(VOL_DIST(KNEIGHBOURS)) #endif DA = SEARCH%KNEAREST(TREE, XTREE, YTREE, ZTREE, & XQUERY=DBLE(RADX(IX)), YQUERY=DBLE(RADY(JY)), @@ -3074,7 +3074,7 @@ SUBROUTINE INTERPOLATE_VELOCITIES(NX,NY,NZ,NL,NPATCH,PARE, ALLOCATE(DIST(KNEIGHBOURS), NEIGH(KNEIGHBOURS)) #if weight_filter == 2 - ALLOCATE(VOL_DIST(KNEIGHBOURS), NEIGH(KNEIGHBOURS)) + ALLOCATE(VOL_DIST(KNEIGHBOURS)) #endif DA = SEARCH%KNEAREST(TREE, XTREE, YTREE, ZTREE, & XQUERY=DBLE(RX(IX,IPATCH)), YQUERY=DBLE(RY(JY,IPATCH)), From bcf31be541b7610734c319d7bee34cc3fbcc3a1a Mon Sep 17 00:00:00 2001 From: Frederick Groth Date: Fri, 29 Aug 2025 15:14:20 +0200 Subject: [PATCH 10/16] fix the range of array indices --- src/readers/gadget_unformatted.f | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/readers/gadget_unformatted.f b/src/readers/gadget_unformatted.f index e79ed9d..13dc845 100644 --- a/src/readers/gadget_unformatted.f +++ b/src/readers/gadget_unformatted.f @@ -160,7 +160,7 @@ subroutine read_gadget_unformatted(iter, files_per_snap, write(*,*) ' found for ',(blocksize-8)/4,' particles' !$omp parallel do shared(npart_gadget, scr4, emissivity, low1), !$omp+ private(i), default(none) - do i=1,sum(npart_gadget(1:6)) + do i=1,npart_gadget(1) ! calculate the weight as rho^2*sqrt(T) ! as this is only used as weight, we don't convert values to actual temperatures! emissivity(low1+i-1) = emissivity(low1+i-1)* From 97124ef71d1c8a43f3862d65847144bcb1764e37 Mon Sep 17 00:00:00 2001 From: Frederick Groth Date: Fri, 29 Aug 2025 16:45:25 +0200 Subject: [PATCH 11/16] do not calculate emissivity in ghost zones --- src/particles.f | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/particles.f b/src/particles.f index dfd9bd7..394edb5 100644 --- a/src/particles.f +++ b/src/particles.f @@ -3157,9 +3157,6 @@ SUBROUTINE INTERPOLATE_VELOCITIES(NX,NY,NZ,NL,NPATCH,PARE, U12(IX,JY,KZ,IPATCH)=BAS8X/BAS8 U13(IX,JY,KZ,IPATCH)=BAS8Y/BAS8 U14(IX,JY,KZ,IPATCH)=BAS8Z/BAS8 -#if weight_filter == 2 - EMISS1(IX,JY,KZ,IPATCH)=BAS8EMISS/BAS8_VOL -#endif if (ix.eq.0.or.ix.eq.n1+1.or. & jy.eq.0.or.jy.eq.n2+1.or. @@ -3171,6 +3168,9 @@ SUBROUTINE INTERPOLATE_VELOCITIES(NX,NY,NZ,NL,NPATCH,PARE, VISC1(IX,JY,KZ,IPATCH)=BAS8M/BAS8 #endif !VISC1(IX,JY,KZ,IPATCH)=BAS8M +#if weight_filter == 2 + EMISS1(IX,JY,KZ,IPATCH)=BAS8EMISS/BAS8_VOL +#endif if (flag_mass.eq.0) then L1(IX,JY,KZ,IPATCH)=H_KERN From 3686241035212247f660b37b45984d651b5cfc15 Mon Sep 17 00:00:00 2001 From: Frederick Groth Date: Tue, 6 Jan 2026 17:34:19 +0100 Subject: [PATCH 12/16] fix some typos in pre-compilation statements. Add some missing grid assignment for the emissivity weighting --- src/particles.f | 24 ++++++++++++++++++++---- 1 file changed, 20 insertions(+), 4 deletions(-) diff --git a/src/particles.f b/src/particles.f index 394edb5..ee86da3 100644 --- a/src/particles.f +++ b/src/particles.f @@ -2319,6 +2319,7 @@ SUBROUTINE INTERPOLATE_VELOCITIES(NX,NY,NZ,NL,NPATCH,PARE, REAL EMISS1(NAMRX,NAMRY,NAMRZ,NPALEV) COMMON /EMISS/ EMISS0,EMISS1 #else + ! Dummy variables REAL EMISS0, EMISS1 #endif @@ -2582,7 +2583,7 @@ SUBROUTINE INTERPOLATE_VELOCITIES(NX,NY,NZ,NL,NPATCH,PARE, #if use_filter == 1 BAS8M=BAS8M+DIST(I)*ABVC(NEIGH(I)) #endif -#if filter_weight == 2 +#if weight_filter == 2 ! emissivity is always volume weighted BAS8_VOL=BAS8_VOL+VOL_DIST(I) BAS8EMISS=BAS8EMISS+VOL_DIST(I)*EMISSIVITY(NEIGH(I)) @@ -2597,8 +2598,8 @@ SUBROUTINE INTERPOLATE_VELOCITIES(NX,NY,NZ,NL,NPATCH,PARE, VISC0(IX,JY,KZ)=BAS8M/BAS8 #endif !VISC0(IX,JY,KZ)=BAS8M -#if filter_weight == 2 - EMISS0(IX,IY,IZ)=BAS8EMISS/BAS8_VOL +#if weight_filter == 2 + EMISS0(IX,JY,KZ)=BAS8EMISS/BAS8_VOL #endif IF (FLAG_MASS.EQ.1) L0(IX,JY,KZ)=BASMASS/(4*PI/3)/H_KERN**3 @@ -2684,7 +2685,7 @@ SUBROUTINE INTERPOLATE_VELOCITIES(NX,NY,NZ,NL,NPATCH,PARE, & VISC0(IIP1,JJP1,KKP1)* BASX * BASY * BASZ #endif -#if filter_weight == 2 +#if weight_filter == 2 EMISS0(IX,JY,KZ)=EMISS0(II,JJ,KK) *(1-BASX)*(1-BASY)*(1-BASZ) + & EMISS0(IIP1,JJ,KK)* BASX *(1-BASY)*(1-BASZ) + & EMISS0(II,JJP1,KK)*(1-BASX)* BASY *(1-BASZ) + @@ -3206,6 +3207,11 @@ SUBROUTINE INTERPOLATE_VELOCITIES(NX,NY,NZ,NL,NPATCH,PARE, & PATCHX,PATCHY,PATCHZ,PATCHRX,PATCHRY,PATCHRZ, & VISC1(1:NAMRX,1:NAMRY,1:NAMRZ,:),NL) #endif +#if weight_filter == 2 + CALL SYNC_AMR_FILTER(IR,NPATCH,PARE,PATCHNX,PATCHNY,PATCHNZ, + & PATCHX,PATCHY,PATCHZ,PATCHRX,PATCHRY,PATCHRZ, + & EMISS1(1:NAMRX,1:NAMRY,1:NAMRZ,:),NL) +#endif LOW1=SUM(NPATCH(0:IR-1))+1 LOW2=SUM(NPATCH(0:IR)) @@ -3244,6 +3250,11 @@ SUBROUTINE INTERPOLATE_VELOCITIES(NX,NY,NZ,NL,NPATCH,PARE, u(1:2,1:2,1:2) = VISC1(I:I+1,J:J+1,K:K+1,IPATCH) call finer_to_coarser(u,uw,fuin) VISC1(II,JJ,KK,JPATCH) = FUIN +#endif +#if weight_filter == 2 + u(1:2,1:2,1:2) = EMISS1(I:I+1,J:J+1,K:K+1,IPATCH) + call finer_to_coarser(u,uw,fuin) + EMISS1(II,JJ,KK,JPATCH) = FUIN #endif else uw(1:2,1:2,1:2) = 1. @@ -3268,6 +3279,11 @@ SUBROUTINE INTERPOLATE_VELOCITIES(NX,NY,NZ,NL,NPATCH,PARE, u(1:2,1:2,1:2) = VISC1(I:I+1,J:J+1,K:K+1,IPATCH) call finer_to_coarser(u,uw,fuin) VISC0(II,JJ,KK) = FUIN +#endif +#if weight_filter == 2 + u(1:2,1:2,1:2) = EMISS1(I:I+1,J:J+1,K:K+1,IPATCH) + call finer_to_coarser(u,uw,fuin) + EMISS0(II,JJ,KK) = FUIN #endif end if END DO From 8182f607ff6a6f5e5507ac0c1a83e13f3d283c51 Mon Sep 17 00:00:00 2001 From: Frederick Groth Date: Tue, 6 Jan 2026 17:34:49 +0100 Subject: [PATCH 13/16] add additional output for emissivity weighting --- src/reader.f | 8 ++++++++ 1 file changed, 8 insertions(+) diff --git a/src/reader.f b/src/reader.f index 9bbff71..a442254 100644 --- a/src/reader.f +++ b/src/reader.f @@ -239,6 +239,10 @@ subroutine read_particles(iter,files_per_snap,nx,ny,nz,t,zeta, end if end if #endif +#if weight_scheme == 3 + write(*,*) 'emis=',minval(emissivity(low1:low2)), + & maxval(emissivity(low1:low2)) +#endif if (xmin.lt.ddxl.or.xmax.gt.ddxr.or. & ymin.lt.ddyl.or.ymax.gt.ddyr.or. @@ -377,6 +381,10 @@ subroutine read_particles(iter,files_per_snap,nx,ny,nz,t,zeta, end if end if #endif +#if weight_scheme == 3 + write(*,*) 'emis=',minval(emissivity(low1:low2)), + & maxval(emissivity(low1:low2)) +#endif write(*,*) 'routine create mesh ------------------------------' From 271d31df7231abdc7928e03ca2a316c699c06a90 Mon Sep 17 00:00:00 2001 From: Frederick Groth Date: Thu, 8 Jan 2026 16:20:27 +0100 Subject: [PATCH 14/16] masclet stores delta_rho. Convert to rho by adding 1 (as done for the density weighting) --- src/readers/masclet.f | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/readers/masclet.f b/src/readers/masclet.f index a4e3c24..1c09f1e 100644 --- a/src/readers/masclet.f +++ b/src/readers/masclet.f @@ -253,7 +253,7 @@ subroutine read_masclet_fields(iter, flag_filter, nx, ny, nz, dens0(1:nx,1:ny,1:nz) = 1.0 + scr4(1:nx,1:ny,1:nz) #elif weight_filter == 2 read(31) (((scr4(i,j,k),i=1,nx),j=1,ny),k=1,nz) - emissivity0(1:nx,1:ny,1:nz) = scr4(1:nx,1:ny,1:nz) + emissivity0(1:nx,1:ny,1:nz) = 1.0 + scr4(1:nx,1:ny,1:nz) #else read(31) #endif @@ -306,7 +306,7 @@ subroutine read_masclet_fields(iter, flag_filter, nx, ny, nz, dens1(1:n1,1:n2,1:n3,ip) = 1.0 + scr4(1:n1,1:n2,1:n3) #elif weight_filter == 2 read(31) (((scr4(i,j,k),i=1,n1),j=1,n2),k=1,n3) - emissivity1(1:n1,1:n2,1:n3,ip) = scr4(1:n1,1:n2,1:n3) + emissivity1(1:n1,1:n2,1:n3,ip) = 1.0 + scr4(1:n1,1:n2,1:n3) #else read(31) #endif From f0e8775dcb3955a0dca625cdf60f79d457ebb63b Mon Sep 17 00:00:00 2001 From: Frederick Groth Date: Thu, 8 Jan 2026 16:24:55 +0100 Subject: [PATCH 15/16] adjust comment, as the weighting can now also be by emissivity --- src/filter.f | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/filter.f b/src/filter.f index f1f6542..7058fb5 100644 --- a/src/filter.f +++ b/src/filter.f @@ -1049,7 +1049,7 @@ subroutine multiscale_filter(nx,ny,nz,nl,npatch,pare, lado0 = nx * dx -* DENS0, DENS1 PROPORTIONAL TO CELLS MASSES or VOLUME!!! +* DENS0, DENS1 PROPORTIONAL TO CELLS MASSES, VOLUME, OR EMISSIVITY!!! dens0 = 1.0 do ir=1,nl low1=sum(npatch(0:ir-1))+1 From ec79a8cf3d17f8012e329c4a4803e6fbf1b96005 Mon Sep 17 00:00:00 2001 From: Frederick Groth Date: Thu, 8 Jan 2026 16:48:05 +0100 Subject: [PATCH 16/16] properly read the emissivity in both cases: emissivity-weighted grid assignment and emissivity-weighted filter, to make the two options independent --- src/particle_data.f | 2 +- src/reader.f | 8 ++++---- 2 files changed, 5 insertions(+), 5 deletions(-) diff --git a/src/particle_data.f b/src/particle_data.f index c52bd07..eb484c8 100644 --- a/src/particle_data.f +++ b/src/particle_data.f @@ -22,7 +22,7 @@ MODULE PARTICLE_DATA REAL VOL #endif -#if weight_scheme == 3 +#if weight_scheme == 3 || weight_filter == 2 REAL,ALLOCATABLE::EMISSIVITY(:) #else REAL EMISSIVITY diff --git a/src/reader.f b/src/reader.f index a442254..6967883 100644 --- a/src/reader.f +++ b/src/reader.f @@ -156,7 +156,7 @@ subroutine read_particles(iter,files_per_snap,nx,ny,nz,t,zeta, #if weight_scheme == 2 || weight_filter == 2 deallocate(vol) #endif -#if weight_scheme == 3 +#if weight_scheme == 3 || weight_filter == 2 deallocate(emissivity) #endif end if @@ -172,7 +172,7 @@ subroutine read_particles(iter,files_per_snap,nx,ny,nz,t,zeta, #if weight_scheme == 2 || weight_filter == 2 allocate(vol(parti)) #endif -#if weight_scheme == 3 +#if weight_scheme == 3 || weight_filter == 2 allocate(emissivity(parti)) #endif @@ -239,7 +239,7 @@ subroutine read_particles(iter,files_per_snap,nx,ny,nz,t,zeta, end if end if #endif -#if weight_scheme == 3 +#if weight_scheme == 3 || weight_filter == 2 write(*,*) 'emis=',minval(emissivity(low1:low2)), & maxval(emissivity(low1:low2)) #endif @@ -381,7 +381,7 @@ subroutine read_particles(iter,files_per_snap,nx,ny,nz,t,zeta, end if end if #endif -#if weight_scheme == 3 +#if weight_scheme == 3 || weight_filter == 2 write(*,*) 'emis=',minval(emissivity(low1:low2)), & maxval(emissivity(low1:low2)) #endif