Commit 476df09b authored by Roy Fabrice's avatar Roy Fabrice
Browse files

returns arrays of properties instead of slices

+ use of some subroutines
+ test where and do concurrent constructs
+ use of allocation_m module
parent 81eec015
......@@ -27,6 +27,8 @@
!> Subroutine to read one RAMSES part.out file
module read_ramses_part_file_m
use common_shared_variables_m, only : mpi_process, timer
implicit none
private
......@@ -34,39 +36,31 @@ module read_ramses_part_file_m
contains
!> Read a part. file produced by Ramses
subroutine Read_ramses_part_file(fpart, lpart, local_nparts, deltasd, nsd, filename, parameters, ramses_particles)
!> Read one particles file produced by Ramses
!>
subroutine Read_ramses_part_file(particles_number, particles, partial_particles_number_per_cube, file_name, parameters)
use allocation_m, only : Alloc
use arrays_of_particle_properties_m, only : arrays_of_particle_properties_t
use common_shared_variables_m, only : mpi_process, timer
use constants_m, only : IDKIND
use error_handling_m, only : ERR_MSG_LEN, Allocate_error, IO_error
use index_m, only : Coords_to_id
use iso_fortran_env, only : ERROR_UNIT
use pfof_snap_parameters_m, only : pfof_snap_parameters_t
integer(kind=4), dimension(mpi_process%comm%size), intent(inout) :: local_nparts
type(arrays_of_particle_properties_t), intent(inout) :: ramses_particles
integer, intent(in) :: fpart
integer, intent(in) :: lpart
real(kind=4), intent(in) :: deltasd
integer(kind=4), intent(in) :: nsd
character(len=*), intent(in) :: filename
type(pfof_snap_parameters_t), intent(in) :: parameters
use iso_fortran_env, only : ERROR_UNIT, INT32
use p2ch_parameters_m, only : p2ch_parameters_t
integer(INT32), intent(out) :: particles_number
type(arrays_of_particle_properties_t), intent(out) :: particles
integer(INT32), intent(out), allocatable, dimension(:) :: partial_particles_number_per_cube
character(len=*), intent(in) :: file_name
type(p2ch_parameters_t), intent(in) :: parameters
integer(INT32) :: cubes_number
character(ERR_MSG_LEN) :: error_message
integer(kind=4) :: ioerr
integer(kind=4) :: part_unit
integer(kind=4) :: ndim
integer(kind=4) :: npart
integer(kind=4) :: allocstat
integer(kind=4) :: idim
integer(kind=4) :: ipart
integer(kind=4) :: n_i, n_j, n_k
integer(kind=IDKIND) :: ind
integer :: ioerr
integer :: part_unit
integer(INT32) :: ndim
integer(INT32) :: idim
integer(kind=4) :: int4_tmp
real(kind=4) :: real4_tmp
real(kind=8) :: real8_tmp
real(kind=8), allocatable, dimension(:) :: real8_tmps
character(*), parameter :: ROUTINE_NAME = 'Read_ramses_part_file'
......@@ -74,9 +68,9 @@ contains
write(ERROR_UNIT,'(a,i0)') ROUTINE_NAME//' begins on progess ', mpi_process%rank
#endif
open(newunit=part_unit, file=filename, status='old', action='read', form='unformatted', iostat=ioerr, iomsg=error_message)
open(newunit=part_unit, file=file_name, status='old', action='read', form='unformatted', iostat=ioerr, iomsg=error_message)
if(ioerr/=0) then
call IO_error('open file '//trim(filename), ROUTINE_NAME, error_message, ioerr, mpi_process%rank)
call IO_error('open file '//trim(file_name), ROUTINE_NAME, error_message, ioerr, mpi_process%rank)
end if
read(part_unit, iostat=ioerr, iomsg=error_message) int4_tmp
......@@ -88,90 +82,28 @@ contains
call IO_error('read ndim', ROUTINE_NAME, error_message, ioerr, mpi_process%rank)
end if
read(part_unit, iostat=ioerr, iomsg=error_message) npart
read(part_unit, iostat=ioerr, iomsg=error_message) particles_number
if(ioerr/=0) then
call IO_error('read npart', ROUTINE_NAME, error_message, ioerr, mpi_process%rank)
end if
ramsesV2 : if(parameters%code_index.eq.'RA2') then
! Read positions
do idim = 1, ndim
read(part_unit, iostat=ioerr, iomsg=error_message) ramses_particles%positions(idim,fpart:lpart)
if(ioerr/=0) then
call IO_error('read ramses_particles%positions', ROUTINE_NAME, error_message, ioerr, mpi_process%rank)
end if
end do
! Read velocities
do idim = 1, ndim
read(part_unit, iostat=ioerr, iomsg=error_message) ramses_particles%velocities(idim,fpart:lpart)
if(ioerr/=0) then
call IO_error('read ramses_particles%velocities', ROUTINE_NAME, error_message, ioerr, mpi_process%rank)
end if
end do
call Alloc(particles%positions, [3, particles_number], 'particles%positions', ROUTINE_NAME)
call Alloc(particles%velocities, [3, particles_number], 'particles%velocities', ROUTINE_NAME)
call Alloc(particles%mass, particles_number, 'particles%mass', ROUTINE_NAME)
if(parameters%do_skip_mass) then
read(part_unit, iostat=ioerr, iomsg=error_message) real4_tmp
if(ioerr/=0) then
call IO_error('read real4_tmp', ROUTINE_NAME, error_message, ioerr, mpi_process%rank)
end if
else
! Read masses if needed
read(part_unit, iostat=ioerr, iomsg=error_message) ramses_particles%mass(fpart:lpart)
if(ioerr/=0) then
call IO_error('read ramses_particles%mass', ROUTINE_NAME, error_message, ioerr, mpi_process%rank)
end if
end if
ramsesV2 : if(parameters%ramses_code.eq.'RA2') then
call Read_RA2_pvm(particles%positions, particles%velocities, particles%mass, part_unit, ndim, particles_number)
end if ramsesV2
ramsesV3 : if(parameters%code_index == 'RA3') then
allocate(real8_tmps(npart), stat=allocstat, errmsg=error_message)
if(allocstat /= 0) then
call Allocate_error('real8_tmps', ROUTINE_NAME,error_message,allocstat,mpi_process%rank)
end if
read(part_unit, iostat=ioerr)
read(part_unit, iostat=ioerr)
read(part_unit, iostat=ioerr)
read(part_unit, iostat=ioerr)
read(part_unit, iostat=ioerr)
! Read positions
do idim = 1, ndim
read(part_unit, iostat=ioerr, iomsg=error_message) real8_tmps
if(ioerr/=0) then
call IO_error('read ramses_particles%positions', ROUTINE_NAME, error_message, ioerr, mpi_process%rank)
end if
ramses_particles%positions(idim,fpart:lpart) = real(real8_tmps, kind=4)
end do
! Read velocities
do idim = 1, ndim
read(part_unit, iostat=ioerr, iomsg=error_message) real8_tmps
if(ioerr/=0) then
call IO_error('read ramses_particles%velocities', ROUTINE_NAME, error_message, ioerr, mpi_process%rank)
end if
ramses_particles%velocities(idim,fpart:lpart) = real(real8_tmps(:),kind=4)
end do
if(parameters%do_skip_mass) then
read(part_unit, iostat=ioerr, iomsg=error_message) real8_tmp
if(ioerr/=0) then
call IO_error('read ramses_particles%mass', ROUTINE_NAME, error_message, ioerr, mpi_process%rank)
end if
else
! Read masses if needed
read(part_unit, iostat=ioerr, iomsg=error_message) real8_tmps
if(ioerr/=0) then
call IO_error('read ramses_particles%mass', ROUTINE_NAME, error_message, ioerr, mpi_process%rank)
end if
ramses_particles%mass(fpart:lpart) = real(real8_tmps, kind=4)
end if
ramsesV3 : if(parameters%ramses_code == 'RA3') then
call Read_RA3_pvm(particles%positions, particles%velocities, particles%mass, part_unit, ndim, particles_number)
end if ramsesV3
call Alloc(particles%identity_p, particles_number, 'particles%identity_p', ROUTINE_NAME)
! Read particle id
read(part_unit, iostat=ioerr, iomsg=error_message) ramses_particles%identity_p(fpart:lpart)
read(part_unit, iostat=ioerr, iomsg=error_message) particles%identity_p
if(ioerr/=0) then
call IO_error('read ramses_particles%identity_p', ROUTINE_NAME, error_message, ioerr, mpi_process%rank)
call IO_error('read particles%identity_p', ROUTINE_NAME, error_message, ioerr, mpi_process%rank)
end if
! First: skip the level (integer 4)
......@@ -180,76 +112,229 @@ contains
call IO_error('read dummy', ROUTINE_NAME, error_message, ioerr, mpi_process%rank)
end if
! Read potential if potential parameters is .true.
if(parameters%do_read_potential) then
! same kind as position and velocity
read(part_unit, iostat=ioerr, iomsg=error_message) ramses_particles%potential(fpart:lpart)
call Alloc(particles%potential, particles_number, 'particles%potential', ROUTINE_NAME)
! Read potential
! same kind as position and velocity
read(part_unit, iostat=ioerr, iomsg=error_message) particles%potential
if(ioerr/=0) then
call IO_error('read particles%potential', ROUTINE_NAME, error_message, ioerr, mpi_process%rank)
end if
call Alloc(particles%gravitational_field, [3, particles_number], 'particles%gravitational_field', ROUTINE_NAME)
! Read force
! same kind as position and velocity
do idim = 1,ndim
read(part_unit, iostat=ioerr, iomsg=error_message) particles%gravitational_field(idim,:)
if(ioerr/=0) then
call IO_error('read ramses_particles%potential', ROUTINE_NAME, error_message, ioerr, mpi_process%rank)
call IO_error('read particles%gravitational_field', ROUTINE_NAME, error_message, ioerr, mpi_process%rank)
end if
end if
end do
! Read force if force parameters is .true.
if(parameters%do_read_gravitational_field) then
! same kind as position and velocity
do idim = 1,ndim
read(part_unit, iostat=ioerr, iomsg=error_message) ramses_particles%gravitational_field(idim,fpart:lpart)
if(parameters%stars_in_files) then
call Alloc(particles%birth_date, particles_number, 'particles%birth_date', ROUTINE_NAME)
read(part_unit, iostat=ioerr, iomsg=error_message) real8_tmps
if(ioerr/=0) then
call IO_error('read particles%birth_date', ROUTINE_NAME, error_message, ioerr, mpi_process%rank)
end if
particles%birth_date(:) = real(real8_tmps(:),kind=4)
if(parameters%metallicity_in_files) then
call Alloc(particles%metallicity, particles_number, 'particles%metallicity', ROUTINE_NAME)
read(part_unit, iostat=ioerr, iomsg=error_message) real8_tmps
if(ioerr/=0) then
call IO_error('read ramses_particles%gravitational_field', ROUTINE_NAME, error_message, ioerr, mpi_process%rank)
call IO_error('read particles%metallicity', ROUTINE_NAME, error_message, ioerr, mpi_process%rank)
end if
end do
particles%metallicity(:) = real(real8_tmps(:),kind=4)
end if
end if
close(part_unit)
deallocate(real8_tmps)
if(parameters%star) then
if(parameters%do_skip_star) then
read(part_unit, iostat=ioerr, iomsg=error_message) real8_tmp
if(ioerr/=0) then
call IO_error('read ramses_particles%birth_date', ROUTINE_NAME, error_message, ioerr, mpi_process%rank)
end if
else
read(part_unit, iostat=ioerr, iomsg=error_message) real8_tmps
if(ioerr/=0) then
call IO_error('read ramses_particles%birth_date', ROUTINE_NAME, error_message, ioerr, mpi_process%rank)
end if
ramses_particles%birth_date(fpart:lpart) = real(real8_tmps(:),kind=4)
cubes_number = (2**parameters%cube_level)**3
call Alloc(partial_particles_number_per_cube, cubes_number, 'partial_particles_number_per_cube', ROUTINE_NAME)
call Alloc(particles%cube_identity_p, particles_number, 'particles%cube_identity_p', ROUTINE_NAME)
call Count_particles_per_cube(particles%cube_identity_p, partial_particles_number_per_cube, particles%positions, parameters%cube_level, ndim, particles_number)
#ifdef DEBUG
write(ERROR_UNIT,'(a,i0)') ROUTINE_NAME//' ends on progess ', mpi_process%rank
#endif
end subroutine Read_ramses_part_file
!------------------------------------------------------------------------------------------------------------------------------------
subroutine Count_particles_per_cube(cube_identity, particles_numbers_per_cube, positions, cube_level, ndim, particles_number)
use constants_m, only : IDKIND
use index_m, only : Coords_to_id
use iso_fortran_env, only : ERROR_UNIT, INT32, REAL32
integer(IDKIND), dimension(*), intent(inout) :: cube_identity
integer(INT32), dimension(*), intent(inout) :: particles_numbers_per_cube
real(REAL32), dimension(ndim,*), intent(in) :: positions
integer(INT32), intent(in) :: cube_level
integer(INT32), intent(in) :: ndim
integer(INT32), intent(in) :: particles_number
integer(INT32) :: cube_grid_dim
real(REAL32) :: cube_size
integer(INT32) :: index
integer(INT32), dimension(ndim) :: indices
integer(INT32) :: ipart
character(*), parameter :: ROUTINE_NAME = 'Count_particles_per_cube'
#ifdef DEBUG
write(ERROR_UNIT,'(a,i0)') ROUTINE_NAME//' begins on progess ', mpi_process%rank
#endif
cube_grid_dim = 2**cube_level
cube_size = 1.e0 / real(cube_grid_dim,kind=4)
! count the particles in each cube
do ipart = 1, particles_number
indices(:) = int(positions(:,ipart)/cube_size)
index = Coords_to_id( indices,(/cube_grid_dim,cube_grid_dim,cube_grid_dim/) )
particles_numbers_per_cube(index) = particles_numbers_per_cube(index) + 1
cube_identity(ipart) = index
end do
do concurrent ( ipart = 1:particles_number )
indices(:) = int(positions(:,ipart)/cube_size)
index = Coords_to_id( indices,(/cube_grid_dim,cube_grid_dim,cube_grid_dim/) )
particles_numbers_per_cube(index) = particles_numbers_per_cube(index) + 1
cube_identity(ipart) = index
end do
#ifdef DEBUG
write(ERROR_UNIT,'(a,i0)') ROUTINE_NAME//' ends on progess ', mpi_process%rank
#endif
end subroutine Count_particles_per_cube
!------------------------------------------------------------------------------------------------------------------------------------
subroutine Read_RA2_pvm(positions, velocities, mass, file_unit, ndim, npart)
use constants_m, only : EPSILON
use error_handling_m, only : ERR_MSG_LEN, IO_error
use iso_fortran_env, only : ERROR_UNIT, INT32, REAL32
real(REAL32), dimension(3,*), intent(out) :: positions
real(REAL32), dimension(3,*), intent(out) :: velocities
real(REAL32), dimension(*), intent(out) :: mass
integer, intent(in) :: file_unit
integer(INT32), intent(in) :: ndim
integer(INT32), intent(in) :: npart
character(ERR_MSG_LEN) :: error_message
integer(INT32) :: idim
integer :: ioerr
character(*), parameter :: ROUTINE_NAME = 'Read_RA2_pvm'
#ifdef DEBUG
write(ERROR_UNIT,'(a,i0)') ROUTINE_NAME//' begins on progess ', mpi_process%rank
#endif
! Read positions
do idim = 1, ndim
read(file_unit, iostat=ioerr, iomsg=error_message) positions(idim,1:npart)
if(ioerr/=0) then
call IO_error('read positions', ROUTINE_NAME, error_message, ioerr, mpi_process%rank)
end if
if(parameters%metal) then
if(parameters%do_skip_metal) then
read(part_unit, iostat=ioerr, iomsg=error_message) real8_tmp
if(ioerr/=0) then
call IO_error('read ramses_particles%metallicity', ROUTINE_NAME, error_message, ioerr, mpi_process%rank)
end if
else
read(part_unit, iostat=ioerr, iomsg=error_message) real8_tmps
if(ioerr/=0) then
call IO_error('read ramses_particles%metallicity', ROUTINE_NAME, error_message, ioerr, mpi_process%rank)
end if
ramses_particles%metallicity(:) = real(real8_tmps(:),kind=4)
end if
where( abs(positions(idim,1:npart)-1.e0) < EPSILON )
positions(idim,1:npart) = 0.e0
end where
end do
! Read velocities
do idim = 1, ndim
read(file_unit, iostat=ioerr, iomsg=error_message) velocities(idim,1:npart)
if(ioerr/=0) then
call IO_error('read velocities', ROUTINE_NAME, error_message, ioerr, mpi_process%rank)
end if
end do
! Read masses if needed
read(file_unit, iostat=ioerr, iomsg=error_message) mass(1:npart)
if(ioerr/=0) then
call IO_error('read mass', ROUTINE_NAME, error_message, ioerr, mpi_process%rank)
end if
deallocate(real8_tmps)
! count the particles in each cube
do ipart = fpart, lpart
do idim = 1, ndim
if(abs(ramses_particles%positions(idim,ipart)-1.0e0) < 1.e-8) ramses_particles%positions(idim,ipart) = 0.e0
end do
n_i = int(ramses_particles%positions(1,ipart)/deltasd)
n_j = int(ramses_particles%positions(2,ipart)/deltasd)
n_k = int(ramses_particles%positions(3,ipart)/deltasd)
ind = Coords_to_id((/n_i,n_j,n_k/),(/nsd,nsd,nsd/))
local_nparts(ind) = local_nparts(ind)+1
ramses_particles%cube_identity_p(ipart) = ind
#ifdef DEBUG
write(ERROR_UNIT,'(a,i0)') ROUTINE_NAME//' ends on progess ', mpi_process%rank
#endif
end subroutine Read_RA2_pvm
!------------------------------------------------------------------------------------------------------------------------------------
subroutine Read_RA3_pvm(positions, velocities, mass, file_unit, ndim, npart)
use allocation_m, only : Alloc
use constants_m, only : EPSILON
use error_handling_m, only : ERR_MSG_LEN, IO_error
use iso_fortran_env, only : ERROR_UNIT, INT32, REAL32, REAL64
real(REAL32), dimension(3,*), intent(out) :: positions
real(REAL32), dimension(3,*), intent(out) :: velocities
real(REAL32), dimension(*), intent(out) :: mass
integer, intent(in) :: file_unit
integer(INT32), intent(in) :: ndim
integer(INT32), intent(in) :: npart
character(ERR_MSG_LEN) :: error_message
integer(INT32) :: idim
integer :: ioerr
real(REAL64), dimension(:), allocatable :: real8_tmps
character(*), parameter :: ROUTINE_NAME = 'Read_RA3_pvm'
#ifdef DEBUG
write(ERROR_UNIT,'(a,i0)') ROUTINE_NAME//' begins on progess ', mpi_process%rank
#endif
call Alloc(real8_tmps, npart, 'real8_tmps', ROUTINE_NAME)
! skip useless info
read(file_unit, iostat=ioerr)
read(file_unit, iostat=ioerr)
read(file_unit, iostat=ioerr)
read(file_unit, iostat=ioerr)
read(file_unit, iostat=ioerr)
! Read positions
do idim = 1, ndim
read(file_unit, iostat=ioerr, iomsg=error_message) real8_tmps
if(ioerr/=0) then
call IO_error('read positions', ROUTINE_NAME, error_message, ioerr, mpi_process%rank)
end if
positions(idim,1:npart) = real(real8_tmps(:), kind=4)
where( abs(positions(idim,1:npart)-1.e0) < EPSILON )
positions(idim,1:npart) = 0.e0
end where
end do
close(part_unit)
! Read velocities
do idim = 1, ndim
read(file_unit, iostat=ioerr, iomsg=error_message) real8_tmps
if(ioerr/=0) then
call IO_error('read velocities', ROUTINE_NAME, error_message, ioerr, mpi_process%rank)
end if
velocities(idim,1:npart) = real(real8_tmps(:),kind=4)
end do
! Read masses if needed
read(file_unit, iostat=ioerr, iomsg=error_message) real8_tmps
if(ioerr/=0) then
call IO_error('read ramses_particles%mass', ROUTINE_NAME, error_message, ioerr, mpi_process%rank)
end if
mass(1:npart) = real(real8_tmps(:), kind=4)
#ifdef DEBUG
write(ERROR_UNIT,'(a,i0)') ROUTINE_NAME//' ends on progess ', mpi_process%rank
#endif
end subroutine Read_ramses_part_file
end subroutine Read_RA3_pvm
end module read_ramses_part_file_m
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment