GHOST (the Geophysical High-Order Suite for Turbulence) is an accurate and highly scalable pseudospectral code that solves a variety of PDEs often encountered in studies of turbulent flows.
I found a bug when using the VDB interface with hybrid (distributed+shared memory) parallelisation.
It is due to how the threads update the particle's id. My solution was to define a new variable this%idj which maps the local id of each thread to the global id. Doing that I found the same with MPI and Hybrid. What I did is simply as follows.
SUBROUTINE GPart_GlobatoLocalIndex(this, gvdb, ngvdb)
!-----------------------------------------------------------------
!-----------------------------------------------------------------
! METHOD : GPart_GlobatoLocalIndex
! DESCRIPTION: Determines a global to local index,
! Call before synchronisation and Copywork
!
! ARGUMENTS :
! this : 'this' class instance (IN)
! gvdb : global VDB containing part. position records. Location
! gives particle id.
! ngvdb : no. records in global VDB
!-----------------------------------------------------------------
USE fprecision
USE commtypes
IMPLICIT NONE
CLASS(GPart), INTENT(INOUT) :: this
INTEGER :: nl
INTEGER, INTENT(IN) :: ngvdb
INTEGER :: i, j
REAL(KIND=GP), INTENT(IN), DIMENSION(3, ngvdb) :: gvdb
nl = 0
this%idj_ = GPNULL
!$omp parallel do
DO j = 1, ngvdb
IF (gvdb(3, j) .GE. this%lxbnds_(3, 1) .AND. gvdb(3, j) .LT. this%lxbnds_(3, 2)) THEN
!$omp critical
nl = nl + 1
this%idj_(j) = nl
!$omp end critical
END IF
END DO
END SUBROUTINE GPart_GlobatoLocalIndex
!-----------------------------------------------------------------
!-----------------------------------------------------------------
SUBROUTINE GPart_GetLocalWrk_aux(this, id, lx, ly, lz, tx, ty, tz, nl, gvdb, gtmp, ngvdb)
!-----------------------------------------------------------------
!-----------------------------------------------------------------
! METHOD : GPart_GetLocalWrk_aux
! DESCRIPTION: Removes from PDB NULL particles, concatenates list,
! and sets new number of particles. This auxiliary
! subroutines also updates arrays used during the
! intermediate steps of the RK solver, and is needed
! if local work is recomputed in the midde of a RK
! iteration.
! ARGUMENTS :
! this : 'this' class instance (IN)
! id : part ids
! lx,ly,lz: local part. d.b. vectors
! tx,ty,tz: local initial part. d.b. vectors
! nl : no. parts. in local pdb
! gvdb : global VDB containing part. position records. Location
! gives particle id.
! gtmp : global VDB containing part. position records at the
! beginning of the RK loop. Location gives particle id.
! ngvdb : no. records in global VDB
!-----------------------------------------------------------------
USE fprecision
USE commtypes
IMPLICIT NONE
CLASS(GPart), INTENT(INOUT) :: this
INTEGER, INTENT(INOUT) :: nl
INTEGER, INTENT(INOUT), DIMENSION(this%maxparts_) :: id
INTEGER, INTENT(IN) :: ngvdb
INTEGER :: i, j
REAL(KIND=GP), INTENT(INOUT), DIMENSION(this%maxparts_) :: lx, ly, lz
REAL(KIND=GP), INTENT(INOUT), DIMENSION(this%maxparts_) :: tx, ty, tz
REAL(KIND=GP), INTENT(IN), DIMENSION(3, ngvdb) :: gvdb
REAL(KIND=GP), INTENT(IN), DIMENSION(3, ngvdb) :: gtmp
nl = 0
id = GPNULL
!$omp parallel do
DO j = 1, ngvdb
IF (gvdb(3, j) .GE. this%lxbnds_(3, 1) .AND. gvdb(3, j) .LT. this%lxbnds_(3, 2)) THEN
!$omp critical
nl = nl + 1
id(this%idj_(j)) = j - 1
lx(this%idj_(j)) = gvdb(1, j)
ly(this%idj_(j)) = gvdb(2, j)
lz(this%idj_(j)) = gvdb(3, j)
tx(this%idj_(j)) = gtmp(1, j)
ty(this%idj_(j)) = gtmp(2, j)
tz(this%idj_(j)) = gtmp(3, j)
!$omp end critical
END IF
END DO
END SUBROUTINE GPart_GetLocalWrk_aux
!-----------------------------------------------------------------
!-----------------------------------------------------------------
SUBROUTINE GPart_CopyLocalWrk(this, lx, ly, lz, gvdb, vgvdb, ngvdb)
!-----------------------------------------------------------------
!-----------------------------------------------------------------
! METHOD : GPart_CopyLocalWrk
! DESCRIPTION: Updates records of the VDB.
! ARGUMENTS :
! this : 'this' class instance (IN)
! lx,ly,lz: local part. d.b. vectors
! gvdb : global VDB containing part. position records. Location
! gives particle id.
! vgvdb : global VDB containing part. property records
! (can be velocity or anything associated to the particle).
! ngvdb : no. records in global VDB
!-----------------------------------------------------------------
USE fprecision
USE commtypes
IMPLICIT NONE
CLASS(GPart), INTENT(INOUT) :: this
INTEGER, INTENT(IN) :: ngvdb
INTEGER :: i, j, nll
REAL(KIND=GP), INTENT(INOUT), DIMENSION(this%maxparts_) :: lx, ly, lz
REAL(KIND=GP), INTENT(IN), DIMENSION(3, ngvdb) :: gvdb
REAL(KIND=GP), INTENT(IN), DIMENSION(3, ngvdb) :: vgvdb
nll = 0
!$omp parallel do
DO j = 1, ngvdb
IF (gvdb(3, j) .GE. this%lxbnds_(3, 1) .AND. gvdb(3, j) .LT. this%lxbnds_(3, 2)) THEN
!$omp critical
nll = nll + 1
lx(this%idj_(j)) = vgvdb(1, j)
ly(this%idj_(j)) = vgvdb(2, j)
lz(this%idj_(j)) = vgvdb(3, j)
!$omp end critical
END IF
END DO
END SUBROUTINE GPart_CopyLocalWrk
!-----------------------------------------------------------------
!-----------------------------------------------------------------
Dear developers,
I found a bug when using the VDB interface with hybrid (distributed+shared memory) parallelisation. It is due to how the threads update the particle's id. My solution was to define a new variable
this%idjwhich maps the local id of each thread to the global id. Doing that I found the same with MPI and Hybrid. What I did is simply as follows.and of course allocation of the variable
ALLOCATE (this%idj_(this%maxparts_))and de-allocation
IF (ALLOCATED(this%idj_)) DEALLOCATE (this%idj_)Thanks!