! File: submodule_sph_particles_handle_positions.f90 ! Authors: Francesco Torsello (FT) !************************************************************************ ! Copyright (C) 2020-2023 Francesco Torsello * ! * ! This file is part of SPHINCS_ID * ! * ! SPHINCS_ID is free software: you can redistribute it and/or modify * ! it under the terms of the GNU General Public License as published by * ! the Free Software Foundation, either version 3 of the License, or * ! (at your option) any later version. * ! * ! SPHINCS_ID is distributed in the hope that it will be useful, * ! but WITHOUT ANY WARRANTY; without even the implied warranty of * ! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * ! GNU General Public License for more details. * ! * ! You should have received a copy of the GNU General Public License * ! along with SPHINCS_ID. If not, see <https://www.gnu.org/licenses/>. * ! The copy of the GNU General Public License should be in the file * ! 'COPYING'. * !************************************************************************ SUBMODULE (sph_particles) handle_positions !************************************************** ! !# This SUBMODULE contains the implementation of ! the PROCEDURES to handle particle positions. ! ! FT 24.03.2022 ! !************************************************** IMPLICIT NONE CONTAINS MODULE PROCEDURE find_particles_above_xy_plane !************************************************************* ! !# Find the particles above the \(xy\) plane ! ! FT 25.03.2022 ! !************************************************************* USE utility, ONLY: zero IMPLICIT NONE INTEGER:: a, npart_diff INTEGER, DIMENSION(npart):: above_xy_plane, below_xy_plane DOUBLE PRECISION, DIMENSION(:), ALLOCATABLE:: & above_xy_plane_a_tmp, below_xy_plane_a_tmp above_xy_plane= zero below_xy_plane= zero npart_above_xy= zero !$OMP PARALLEL DO DEFAULT( NONE ) & !$OMP SHARED( pos, above_xy_plane, below_xy_plane, npart ) & !$OMP PRIVATE( a ) & !$OMP REDUCTION(+:npart_above_xy) DO a= 1, npart, 1 IF( pos(3,a) > zero )THEN npart_above_xy= npart_above_xy + 1 above_xy_plane(a)= a ELSEIF( pos(3,a) < zero )THEN below_xy_plane(a)= a ENDIF ENDDO !$OMP END PARALLEL DO npart_diff= 0 IF( npart/2 /= npart_above_xy )THEN PRINT *, "** WARNING! Mismatch in the number of particles above the xy ",& "plane in SUBROUTINE reflect_particles_xy_plane!" PRINT *, " * npart/2= ", npart/2 PRINT *, " * npart_above_xy= ", npart_above_xy PRINT *, " * If you are inside the APM iteration, you're safe since ", & "this is taken care of." PRINT *, " Otherwise, you may want to double check that ", & "you know what's going on." PRINT * IF( npart/2 > npart_above_xy )THEN npart_diff= npart/2 - npart_above_xy ENDIF npart_above_xy= npart/2 ENDIF ALLOCATE(above_xy_plane_a(npart_above_xy)) above_xy_plane_a_tmp= PACK( above_xy_plane, above_xy_plane /= zero ) below_xy_plane_a_tmp= PACK( below_xy_plane, below_xy_plane /= zero ) IF( npart_diff > 0 )THEN !$OMP PARALLEL DO DEFAULT( NONE ) & !$OMP SHARED( pos, above_xy_plane_a, above_xy_plane_a_tmp, npart ) & !$OMP PRIVATE( a ) DO a= 1, SIZE(above_xy_plane_a_tmp), 1 above_xy_plane_a(a)= above_xy_plane_a_tmp(a) ENDDO !$OMP END PARALLEL DO IF( npart_above_xy /= SIZE(above_xy_plane_a_tmp) + npart_diff )THEN PRINT *, "** ERROR! Mismatch in the number of particles above the xy ",& "plane in SUBROUTINE reflect_particles_xy_plane!" PRINT *, " * SIZE(above_xy_plane_a_tmp)= ", SIZE(above_xy_plane_a_tmp) PRINT *, " * npart_diff= ", npart_diff PRINT *, " * npart_above_xy= ", npart_above_xy PRINT *, " * Stopping..." PRINT * STOP ENDIF !$OMP PARALLEL DO DEFAULT( NONE ) & !$OMP SHARED( pos, above_xy_plane_a, above_xy_plane_a_tmp, & !$OMP npart_above_xy, below_xy_plane_a_tmp ) & !$OMP PRIVATE( a ) DO a= SIZE(above_xy_plane_a_tmp) + 1, npart_above_xy, 1 above_xy_plane_a(a)= below_xy_plane_a_tmp(a) ENDDO !$OMP END PARALLEL DO ELSE !$OMP PARALLEL DO DEFAULT( NONE ) & !$OMP SHARED( pos, above_xy_plane_a, above_xy_plane_a_tmp, npart_above_xy ) & !$OMP PRIVATE( a ) DO a= 1, npart_above_xy, 1 above_xy_plane_a(a)= above_xy_plane_a_tmp(a) ENDDO !$OMP END PARALLEL DO ENDIF END PROCEDURE find_particles_above_xy_plane MODULE PROCEDURE reflect_particles_xy_plane !************************************************************* ! !# Reflect the particle with z>0 with respect to the xy plane ! ! FT 25.03.2022 ! !************************************************************* IMPLICIT NONE INTEGER:: a DOUBLE PRECISION, DIMENSION(3,npart):: pos_tmp DOUBLE PRECISION, DIMENSION(npart) :: nu_tmp !CHARACTER( LEN= : ), ALLOCATABLE:: finalnamefile !LOGICAL:: exist IF( PRESENT(nu) .NEQV. PRESENT(nu_below) )THEN PRINT *, "** ERROR! In SUBROUTINE reflect_particles_xy_plane, the ", & "arguments 'nu' and 'nu_below must be either both present or ", & "both absent." PRINT *, " * Stopping..." PRINT * STOP ENDIF pos_tmp= pos IF( PRESENT(nu) ) nu_tmp = nu !$OMP PARALLEL DO DEFAULT( NONE ) & !$OMP SHARED( pos, pos_tmp, above_xy_plane_a, npart_above_xy, & !$OMP nu, nu_tmp ) & !$OMP PRIVATE( a ) DO a= 1, npart_above_xy, 1 pos( :, a )= pos_tmp( :, above_xy_plane_a(a) ) IF( PRESENT(nu) ) nu( a )= nu_tmp( above_xy_plane_a(a) ) ENDDO !$OMP END PARALLEL DO !$OMP PARALLEL DO DEFAULT( NONE ) & !$OMP SHARED( pos, pos_below, npart_above_xy, nu, nu_below ) & !$OMP PRIVATE( a ) DO a= 1, npart_above_xy, 1 pos_below( 1, a )= pos( 1, a ) pos_below( 2, a )= pos( 2, a ) pos_below( 3, a )= - pos( 3, a ) IF( PRESENT(nu) ) nu_below( a )= nu( a ) ENDDO !$OMP END PARALLEL DO END PROCEDURE reflect_particles_xy_plane MODULE PROCEDURE impose_equatorial_plane_symmetry !************************************************************* ! !# Mirror the particle with z>0 with respect to the xy plane, ! to impose the equatorial-plane symmetry ! ! FT 1.09.2021 ! !************************************************************* USE analyze, ONLY: COM IMPLICIT NONE INTEGER:: a, npart_above_xy DOUBLE PRECISION:: com_x, com_y, com_z, com_d INTEGER, DIMENSION(:), ALLOCATABLE:: above_xy_plane_a DOUBLE PRECISION, DIMENSION(:,:), ALLOCATABLE:: pos_below DOUBLE PRECISION, DIMENSION(:), ALLOCATABLE:: nu_below ! DOUBLE PRECISION, DIMENSION(3,npart+npart_ghost):: pos_tmp ! DOUBLE PRECISION, DIMENSION(npart+npart_ghost) :: nu_tmp IF( MOD(npart,2) /= 0 )THEN PRINT *, "** ERROR! If the equatorial symmetry has to be imposed, ", & "the particle number must be even." PRINT *, " * npart= ", npart PRINT *, " * Stopping..." PRINT * STOP ENDIF ! IF( PRESENT(pos_prev) )THEN ! ! ! If some of the particles crossed the xy plane top-down in the ! ! last step, replace their z coordinate with their previous ! ! z coordinate ! !$OMP PARALLEL DO DEFAULT( NONE ) & ! !$OMP SHARED( pos, pos_prev, npart ) & ! !$OMP PRIVATE( a ) ! DO a= 1, npart, 1 ! ! IF( (pos_prev(3,a) > zero) .AND. (pos(3,a) <= zero) )THEN ! ! pos(3,a)= pos_prev(3,a) ! ! ENDIF ! ! ENDDO ! !$OMP END PARALLEL DO ! ! ENDIF ! pos_tmp= pos ! nu_tmp = nu ! itr= 0 ! DO a= 1, npart, 1 ! ! IF( pos_tmp( 3, a ) > zero & ! .AND. & ! itr <= npart/2 )THEN ! ! itr= itr + 1 ! pos( 1, itr )= pos_tmp( 1, a ) ! pos( 2, itr )= pos_tmp( 2, a ) ! pos( 3, itr )= pos_tmp( 3, a ) ! IF( PRESENT(nu) ) nu( itr )= nu_tmp( a ) ! ! ENDIF ! ! ENDDO ! npart_above_xy= itr CALL find_particles_above_xy_plane( npart, pos, npart_above_xy, & above_xy_plane_a ) ! IF(npart/2 /= npart_above_xy )THEN ! ! ! ! ENDIF ALLOCATE( pos_below(3,npart_above_xy) ) IF( PRESENT(nu) )THEN ALLOCATE( nu_below(npart_above_xy) ) CALL reflect_particles_xy_plane( npart, pos, pos_below, npart_above_xy, & above_xy_plane_a, nu, nu_below ) ELSE CALL reflect_particles_xy_plane( npart, pos, pos_below, npart_above_xy, & above_xy_plane_a ) ENDIF !$OMP PARALLEL DO DEFAULT( NONE ) & !$OMP SHARED( pos, npart_above_xy, nu, & !$OMP pos_below, nu_below ) & !$OMP PRIVATE( a ) DO a= 1, npart_above_xy, 1 pos( :, npart_above_xy + a )= pos_below( :, a ) IF( PRESENT(nu) )THEN nu( npart_above_xy + a )= nu_below( a ) ENDIF ENDDO !$OMP END PARALLEL DO IF( PRESENT(verbose) )THEN IF( verbose .EQV. .TRUE. )THEN CALL COM( npart, pos, nu, & ! input com_x, com_y, com_z, com_d ) ! output PRINT *, "** After mirroring particles:" IF( PRESENT(com_star) ) PRINT *, & " * x coordinate of the center of mass of the star, ", & "from LORENE: com_star= ", com_star, "Msun_geo" PRINT *, " * x coordinate of the center of mass of the particle ", & "distribution: com_x= ", com_x, "Msun_geo" PRINT *, " * y coordinate of the center of mass of the particle ", & "distribution: com_y= ", com_y, "Msun_geo" PRINT *, " * z coordinate of the center of mass of the particle ", & "distribution: com_z= ", com_z, "Msun_geo" PRINT *, " * Distance of the center of mass of the particle ", & "distribution from the origin: com_d= ", com_d IF( PRESENT(com_star) ) PRINT *, " * |com_x-com_star/com_star|=", & ABS( com_x-com_star )/ABS( com_star + 1 ) PRINT * ENDIF ENDIF END PROCEDURE impose_equatorial_plane_symmetry MODULE PROCEDURE check_particle_positions !************************************************* ! !# Check that the particles are not at the same ! positions ! ! FT 1.9.2021 ! !************************************************* USE NR, ONLY: indexx IMPLICIT NONE INTEGER:: itr !! Iterator INTEGER:: itr2 !! Iterator INTEGER:: x_idx !# Index at which a new value of the \(x\) coordinate appears, ! in the array `pos` sorted so that the \(x\) coordinate does not decrease INTEGER, DIMENSION(:), ALLOCATABLE:: x_sort !# Array storing the sorted indices of array `pos`, so that the \(x\) ! coordinate of the particles is in nondecreasing order INTEGER, DIMENSION(:), ALLOCATABLE:: x_number !# Array storing, for each \(x\) coordinate, the number of particles ! having that \(x\) coordinate PRINT *, "** Checking that there are not multiple particles", & " at the same position..." PRINT * ALLOCATE( x_sort( npart ) ) ALLOCATE( x_number( npart ) ) ! Sort x coordinates of the particles CALL indexx( npart, pos( 1, : ), x_sort ) x_number= 1 itr2= 1 ! Find the number of times each x appears ! TODO: parallelize this DO itr= 1, npart - 1, 1 IF( pos( 1, x_sort(itr) ) == & pos( 1, x_sort(itr+1) ) )THEN x_number(itr2)= x_number(itr2) + 1 ELSE itr2= itr2 + 1 ENDIF ENDDO x_number= x_number(1:itr2) IF( SUM( x_number ) /= npart )THEN PRINT *, "** ERROR! The sum of the numbers of particles with the same", & " x is not equal to the particle number." PRINT *, " * SUM( x_number )=", SUM( x_number ), ", ", & "npart=", npart PRINT *, " * Stopping..." PRINT * STOP ENDIF IF( PRESENT(debug) )THEN IF( debug .EQV. .TRUE. )THEN ! These two IF statements are nested because gfortran doesn't like them ! together when debug is not PRESENT !$OMP PARALLEL DO DEFAULT( NONE ) & !$OMP SHARED( pos, x_sort, x_number ) & !$OMP PRIVATE( itr, itr2, x_idx ) DO itr= 1, SIZE(x_number), 1 IF( itr == 1 )THEN x_idx= 1 ELSE x_idx= SUM(x_number(1:itr-1)) + 1 ENDIF DO itr2= x_idx, x_idx + x_number(itr) - 2, 1 ! If they do not have the same x IF( pos( 1, x_sort(itr2) ) /= & pos( 1, x_sort(itr2+1) ) )THEN PRINT *, "** ERROR! ", "The two particles ", x_sort(itr2), & " and", x_sort(itr2+1), & " do not have the same x, but should!" PRINT *, pos( :, x_sort(itr2) ) PRINT *, pos( :, x_sort(itr2+1) ) PRINT *, " * Stopping..." PRINT * STOP ENDIF ENDDO ENDDO !$OMP END PARALLEL DO ENDIF ENDIF !$OMP PARALLEL DO DEFAULT( NONE ) & !$OMP SHARED( pos, x_sort, x_number ) & !$OMP PRIVATE( itr, itr2, x_idx ) DO itr= 1, SIZE(x_number), 1 IF( itr == 1 )THEN x_idx= 1 ELSE x_idx= SUM(x_number(1:itr-1)) + 1 ENDIF DO itr2= x_idx, x_idx + x_number(itr) - 2, 1 ! If they have the same y IF( pos( 2, x_sort(itr2) ) == & pos( 2, x_sort(itr2+1) ) )THEN ! If they have the same z IF( pos( 3, x_sort(itr2) ) == & pos( 3, x_sort(itr2+1) ) )THEN ! They are the same PRINT *, "** ERROR! ", "The two particles ", x_sort(itr2), & " and", x_sort(itr2+1), " have the same coordinates!" PRINT *, pos( :, x_sort(itr2) ) PRINT *, pos( :, x_sort(itr2+1) ) PRINT *, " * Stopping..." PRINT * STOP ENDIF ENDIF ENDDO ENDDO !$OMP END PARALLEL DO DEALLOCATE( x_sort ) DEALLOCATE( x_number ) END PROCEDURE check_particle_positions MODULE PROCEDURE check_particle_position !***************************************************** ! !# Return the number of times that pos_a appears ! in the array pos ! @todo This algorithm scales as O(npart**2) ! if used in a loop over the particles... ! To be documented, after it's fixed ! ! FT 13.10.2021 ! !***************************************************** !USE NR, ONLY: indexx IMPLICIT NONE INTEGER:: itr, itr2, size_x!, cnt INTEGER, DIMENSION(npart):: x_sort, cnts INTEGER, DIMENSION(npart):: x_number INTEGER, DIMENSION(:), ALLOCATABLE:: x_number_filt ! Sort x coordinates of the particles !CALL indexx( npart, pos( 1, : ), x_sort ) x_number= 0 itr2= 0 ! Find the number of times that the x coordinate of pos_a appears in pos !$OMP PARALLEL DO DEFAULT( NONE ) & !$OMP SHARED( pos, pos_a, x_sort, x_number, npart ) & !$OMP PRIVATE( itr ) DO itr= 1, npart, 1 IF( pos( 1, itr ) == pos_a( 1 ) )THEN !itr2= itr2 + 1 x_number(itr)= itr !ELSEIF( pos( 1, x_sort(itr) ) > pos_a( 1 ) )THEN ! ! EXIT ENDIF ENDDO !$OMP END PARALLEL DO x_number_filt= PACK( x_number, x_number /= 0 ) size_x= SIZE(x_number_filt) cnts= 0 !$OMP PARALLEL DO DEFAULT( NONE ) & !$OMP SHARED( pos, pos_a, x_sort, x_number_filt, size_x, cnts )& !$OMP PRIVATE( itr ) DO itr= 1, size_x, 1 ! If they have the same y IF( pos( 2, x_number_filt(itr) ) == pos_a( 2 ) )THEN ! If they have the same z IF( pos( 3, x_number_filt(itr) ) == pos_a( 3 ) )THEN cnts(itr)= cnts(itr) + 1 ENDIF ENDIF ENDDO !$OMP END PARALLEL DO cnt= SUM( cnts ) END PROCEDURE check_particle_position MODULE PROCEDURE correct_center_of_mass !*********************************************************** ! !# Translate the particles so that their center of mass ! coincides with the center of mass of the star, given by ! |id| ! ! FT 1.09.2021 ! !*********************************************************** USE analyze, ONLY: COM USE utility, ONLY: is_finite_number, zero IMPLICIT NONE INTEGER:: a, i DOUBLE PRECISION:: com_x, com_y, com_z, com_d DOUBLE PRECISION, DIMENSION(3):: pos_corr_tmp CALL COM( npart, pos, nu, & ! input com_x, com_y, com_z, com_d ) ! output IF( PRESENT(verbose) )THEN IF( verbose .EQV. .TRUE. )THEN PRINT *, "** Before center of mass correction:" PRINT *, " * x coordinate of the center of mass of the star, ", & "from the ID: com_star= ", com_star(1), "Msun_geo" PRINT *, " * y coordinate of the center of mass of the star, ", & "from the ID: com_star= ", com_star(2), "Msun_geo" PRINT *, " * z coordinate of the center of mass of the star, ", & "from the ID: com_star= ", com_star(3), "Msun_geo" PRINT *, " * x coordinate of the center of mass of the particle ", & "distribution: com_x= ", com_x, "Msun_geo" PRINT *, " * y coordinate of the center of mass of the particle ", & "distribution: com_y= ", com_y, "Msun_geo" PRINT *, " * z coordinate of the center of mass of the particle ", & "distribution: com_z= ", com_z, "Msun_geo" PRINT *, " * Distance of the center of mass of the particle ", & "distribution from the origin: com_d= ", com_d PRINT *, " * |com_x-com_star_x/com_star_x|=", & ABS( com_x-com_star(1) )/ABS( com_star(1) + 1 ) PRINT * ENDIF ENDIF !$OMP PARALLEL DO DEFAULT( NONE ) & !$OMP SHARED( pos, com_star, & !$OMP com_x, com_y, com_z, npart ) & !$OMP PRIVATE( pos_corr_tmp, a ) DO a= 1, npart, 1 pos_corr_tmp(1)= pos(1,a) - ( com_x - com_star(1) ) pos_corr_tmp(2)= pos(2,a) - ( com_y - com_star(2) ) pos_corr_tmp(3)= pos(3,a) - ( com_z - com_star(3) ) DO i= 1, 3, 1 IF( .NOT.is_finite_number(pos_corr_tmp(i)) )THEN PRINT *, pos_corr_tmp STOP ENDIF ENDDO IF( get_density( & pos_corr_tmp(1), pos_corr_tmp(2), pos_corr_tmp(3) ) > zero & .AND. & validate_pos( & pos_corr_tmp(1), pos_corr_tmp(2), pos_corr_tmp(3) ) & )THEN pos(:,a)= pos_corr_tmp ENDIF ENDDO !$OMP END PARALLEL DO CALL COM( npart, pos, nu, & ! input com_x, com_y, com_z, com_d ) ! output IF( PRESENT(verbose) )THEN IF( verbose .EQV. .TRUE. )THEN PRINT *, "** After center of mass correction:" PRINT *, " * x coordinate of the center of mass of the star, ", & "from the ID: com_star= ", com_star(1), "Msun_geo" PRINT *, " * y coordinate of the center of mass of the star, ", & "from the ID: com_star= ", com_star(2), "Msun_geo" PRINT *, " * z coordinate of the center of mass of the star, ", & "from the ID: com_star= ", com_star(3), "Msun_geo" PRINT *, " * x coordinate of the center of mass of the particle ", & "distribution: com_x= ", com_x, "Msun_geo" PRINT *, " * y coordinate of the center of mass of the particle ", & "distribution: com_y= ", com_y, "Msun_geo" PRINT *, " * z coordinate of the center of mass of the particle ", & "distribution: com_z= ", com_z, "Msun_geo" PRINT *, " * Distance of the center of mass of the particle ", & "distribution from the origin: com_d= ", com_d PRINT *, " * |com_x-com_star_x/com_star_x|=", & ABS( com_x-com_star(1) )/ABS( com_star(1) + 1 ) PRINT * ENDIF ENDIF END PROCEDURE correct_center_of_mass MODULE PROCEDURE get_neighbours_bf !************************************************************** ! !# just for test purposes: get neighbours of particle ipart in ! a "brute force" way; ipart is ALSO on the neighbour list; ! SKR 8.2.2010 ! ! Removed ipart from its own neighbors' list ! FT 04.06.2021 ! !************************************************************** IMPLICIT NONE INTEGER:: a DOUBLE PRECISION:: diff(dimensions), d2, r_int2 ! square of interaction radius r_int2= (2.D0*h(ipart))**2 nnei= 0 !$OMP PARALLEL DO DEFAULT(NONE) & !$OMP SHARED(pos,dimensions,ipart,npart,r_int2,nnei,neilist)& !$OMP PRIVATE(a,diff,d2) DO a= 1, npart, 1 IF( a /= ipart )THEN diff= pos(1:dimensions,a)-pos(1:dimensions,ipart) d2= DOT_PRODUCT(diff,diff) ! neighbour? IF(d2 < r_int2)THEN nnei= nnei + 1 neilist(nnei)= a ENDIF ENDIF ENDDO !$OMP END PARALLEL DO END PROCEDURE get_neighbours_bf MODULE PROCEDURE get_object_of_particle !************************************************************** ! !# Returns the number of the matter object asociated with the particle ! number given as input. Example: give number \(n\) as input; this ! particle number corresponds to a particle on matter object \(m\). ! This functions returns \(m\). ! ! FT 05.12.2022 ! !************************************************************** IMPLICIT NONE INTEGER:: i_matter DO i_matter= 1, this% n_matter, 1 IF(a > this% npart_fin(i_matter - 1) & .AND. a <= this% npart_fin(i_matter))THEN get_object_of_particle= i_matter RETURN ENDIF ENDDO END PROCEDURE get_object_of_particle END SUBMODULE handle_positions