Skip to content
Snippets Groups Projects
Commit dd7b0757 authored by odemirel's avatar odemirel
Browse files

- move comp_pp_* to libppmnumerics (ticket #27) 

git-svn-id: https://ppm.inf.ethz.ch/svn/ppmnumerics/branches/ngtopo/libppmnumerics@521 7c7fe9aa-52eb-4d9e-b0a8-ba7d787348e9
parent 517b63b1
No related branches found
No related tags found
No related merge requests found
Expand all Hide whitespace changes
Inline Side-by-side
Showing
with 3601 additions and 0 deletions
src/ppm_comp_pp_cell.f 0 → 100644
+ 799
− 0
View file @ dd7b0757
!-------------------------------------------------------------------------
! Subroutine : ppm_comp_pp_cell
!-------------------------------------------------------------------------
! Parallel Particle Mesh Library (PPM)
! ETH Zurich
! CH-8092 Zurich, Switzerland
!-------------------------------------------------------------------------
#if __KERNEL == __INTERNAL
#if __KIND == __SINGLE_PRECISION
SUBROUTINE ppm_comp_pp_cell_si(xp,Np,pdata,lda,lsymm,kernel,kpar, &
& nsublist,Nm,clist,cutoff2,dpd,info)
#elif __KIND == __DOUBLE_PRECISION
SUBROUTINE ppm_comp_pp_cell_di(xp,Np,pdata,lda,lsymm,kernel,kpar, &
& nsublist,Nm,clist,cutoff2,dpd,info)
#elif __KIND == __SINGLE_PRECISION_COMPLEX
SUBROUTINE ppm_comp_pp_cell_sci(xp,Np,pdata,lda,lsymm,kernel,kpar, &
& nsublist,Nm,clist,cutoff2,dpd,info)
#elif __KIND == __DOUBLE_PRECISION_COMPLEX
SUBROUTINE ppm_comp_pp_cell_dci(xp,Np,pdata,lda,lsymm,kernel,kpar, &
& nsublist,Nm,clist,cutoff2,dpd,info)
#endif
#elif __KERNEL == __USER_FUNCTION
#if __KIND == __SINGLE_PRECISION
SUBROUTINE ppm_comp_pp_cell_su(xp,Np,pdata,lda,lsymm,kernel,kpar, &
& nsublist,Nm,clist,cutoff2,dpd,info)
#elif __KIND == __DOUBLE_PRECISION
SUBROUTINE ppm_comp_pp_cell_du(xp,Np,pdata,lda,lsymm,kernel,kpar, &
& nsublist,Nm,clist,cutoff2,dpd,info)
#elif __KIND == __SINGLE_PRECISION_COMPLEX
SUBROUTINE ppm_comp_pp_cell_scu(xp,Np,pdata,lda,lsymm,kernel,kpar, &
& nsublist,Nm,clist,cutoff2,dpd,info)
#elif __KIND == __DOUBLE_PRECISION_COMPLEX
SUBROUTINE ppm_comp_pp_cell_dcu(xp,Np,pdata,lda,lsymm,kernel,kpar, &
& nsublist,Nm,clist,cutoff2,dpd,info)
#endif
#elif __KERNEL == __LOOKUP_TABLE
#if __KIND == __SINGLE_PRECISION
SUBROUTINE ppm_comp_pp_cell_st(xp,Np,pdata,lda,lsymm,kernel,kpar, &
& nsublist,Nm,clist,cutoff2,dpd,info)
#elif __KIND == __DOUBLE_PRECISION
SUBROUTINE ppm_comp_pp_cell_dt(xp,Np,pdata,lda,lsymm,kernel,kpar, &
& nsublist,Nm,clist,cutoff2,dpd,info)
#elif __KIND == __SINGLE_PRECISION_COMPLEX
SUBROUTINE ppm_comp_pp_cell_sct(xp,Np,pdata,lda,lsymm,kernel,kpar, &
& nsublist,Nm,clist,cutoff2,dpd,info)
#elif __KIND == __DOUBLE_PRECISION_COMPLEX
SUBROUTINE ppm_comp_pp_cell_dct(xp,Np,pdata,lda,lsymm,kernel,kpar, &
& nsublist,Nm,clist,cutoff2,dpd,info)
#endif
#endif
!!! This routine computes kernel interactions by direct
!!! particle-particle interactions using cell lists.
!!!
!!! [NOTE]
!!! Loops over lda have been explicitly unrolled for the cases of
!!! `lda.EQ.1` and `lda.EQ.2` to allow vectorization in these cases.
!!!
!!! [WARNING]
!!! clist has to be created by the user program before
!!! calling this routine. ppm_neighlist_clist should be
!!! used for this. The reason is that this allows the
!!! user to have multiple cell lists and call this
!!! routine for the appropriate one. Do not forget to
!!! call ppm_neighlist_clist_destroy and to DELLOCATE
!!! Nm after the cell lists are no longer needed.
!!!
!!! [NOTE]
!!! dpd needs to be allocated to proper size before
!!! calling this routine. Also, this routine is not
!!! resetting dpd to zero before doing the PP
!!! interactions. This allows contributions from
!!! different kernels to be accumulated. If needed,
!!! set it to zero before calling this routine the
!!! first time.
!-------------------------------------------------------------------------
! Modules
!-------------------------------------------------------------------------
USE ppm_module_data
USE ppm_module_data_neighlist
USE ppm_module_substart
USE ppm_module_substop
USE ppm_module_error
USE ppm_module_neighlist_MkNeighIdx
IMPLICIT NONE
#if __KIND == __SINGLE_PRECISION | __KIND == __SINGLE_PRECISION_COMPLEX
INTEGER, PARAMETER :: MK = ppm_kind_single
#elif __KIND == __DOUBLE_PRECISION | __KIND == __DOUBLE_PRECISION_COMPLEX
INTEGER, PARAMETER :: MK = ppm_kind_double
#endif
!-------------------------------------------------------------------------
! Includes
!-------------------------------------------------------------------------
#include "ppm_define.h"
!-------------------------------------------------------------------------
! Arguments
!-------------------------------------------------------------------------
REAL(MK) , DIMENSION(:,:), INTENT(IN ) :: xp
!!! particle co-ordinates
#if __KIND == __SINGLE_PRECISION | __KIND == __DOUBLE_PRECISION
REAL(MK) , DIMENSION(:,:), INTENT(IN ) :: pdata
!!! particle strengths used for interaction.
!!! Overloaded type: single,double
REAL(MK) , DIMENSION(:,:), INTENT( OUT) :: dpd
!!! Change of particle data (pdata) due to interaction.
!!! This is not initialized by this routine!
!!! Overloaded types: single,double.
#elif __KIND == __SINGLE_PRECISION_COMPLEX | __KIND == __DOUBLE_PRECISION_COMPLEX
COMPLEX(MK), DIMENSION(:,:), INTENT(IN ) :: pdata
!!! particle strengths used for interaction.
!!! Overloaded types: single complex, double complex.
COMPLEX(MK), DIMENSION(:,:), INTENT( OUT) :: dpd
!!! Change of particle data (pdata) due to interaction.
!!! This is not initialized by this routine!
!!! Overloaded types: single complex, double complex.
#endif
INTEGER , INTENT(IN ) :: Np
!!! number of particles on this proc.
INTEGER , INTENT(IN ) :: lda
!!! lading dimension of pdata.
INTEGER , INTENT(IN ) :: nsublist
!!! number of subdomains on the local processor
#if __KERNEL == __INTERNAL
INTEGER , INTENT(IN ) :: kernel
!!! kernel to be used for PP interactions. To use ppm-internal
!!! kernels, specify one of:
!!!
!!! ---------------------------------------
!!! ppm_param_kerel_laplace2d_2p
!!! (2nd order Laplacian, polynomial in 2D)
!!! ppm_param_kerel_laplace3d_2p
!!! (2nd order Laplacian, polynomial in 3D)
!!! ---------------------------------------
!!!
!!! To use your own kernel function, pass the function pointer here.
!!! Your function should take one argument and return one value.
!!! The third possibility is to pass a lookup table with tabulated
!!! kernel values. Such a table can be created using
!!! ppm_comp_pp_mk_table.
#if __KIND == __SINGLE_PRECISION | __KIND == __DOUBLE_PRECISION
REAL(MK) , DIMENSION(:) , INTENT(IN ) :: kpar
!!! Kernel parameters. See documentation or ppm_comp_pp_kernels.inc
!!! for description. Type can be single, double
#elif __KIND == __SINGLE_PRECISION_COMPLEX | __KIND == __DOUBLE_PRECISION_COMPLEX
COMPLEX(MK), DIMENSION(:) , INTENT(IN ) :: kpar
!!! Kernel parameters. See documentation or ppm_comp_pp_kernels.inc
!!! for description. Type can be single complex or double complex.
#endif
#elif __KERNEL == __LOOKUP_TABLE
REAL(MK) , INTENT(IN ) :: kpar
!!! Kernel parameters. Lookup table version.
!!! Pass dxtableinv (the inverse of the table spacing) as a scalar here.
#if __KIND == __SINGLE_PRECISION | __KIND == __DOUBLE_PRECISION
REAL(MK) , DIMENSION(: ), INTENT(IN ) :: kernel
#elif __KIND == __SINGLE_PRECISION_COMPLEX | __KIND == __DOUBLE_PRECISION_COMPLEX
COMPLEX(MK), DIMENSION(: ), INTENT(IN ) :: kernel
#endif
!!! lookup table with tabulated kernel values. Such a table can
!!! be created using <<ppm_comp_pp_mk_table>>.
#endif
LOGICAL , INTENT(IN ) :: lsymm
!!! Do symmetric PP interactions?
TYPE(ppm_type_ptr_to_clist), DIMENSION(:), INTENT(IN) :: clist
!!! Cell list as a list of ptr_to_clist.
!!!
!!! particle index list of isub: `clist(isub)%lpdx(:)` +
!!! pointer to first particle in ibox of isub: ` clist(isub)%lhbx(ibox)`
INTEGER , DIMENSION(:,:), INTENT(IN ) :: Nm
!!! number of cells in x,y,(z) directions (including the ghosts cells)
!!! in each subdomain.
!!!
!!! 1st index: direction. second index: subid.
REAL(MK) , INTENT(IN ) :: cutoff2
!!! Squared PP interaction cutoff. Should be .LE. cell size.
INTEGER , INTENT( OUT) :: info
!!! Returns status, 0 upon success
!-------------------------------------------------------------------------
! Local variables
!-------------------------------------------------------------------------
! counters
INTEGER :: i,idom,ibox,jbox,idx
INTEGER :: ipart,jpart,ip,jp
INTEGER :: cbox,iinter,j,k,ispec
! coordinate differences
REAL(MK) :: dx,dy,dz
REAL(MK) :: factor,factor2
#if __KIND == __SINGLE_PRECISION | __KIND == __DOUBLE_PRECISION
REAL(MK) :: summ,summ2
REAL(MK) :: eta,dm
#elif __KIND == __SINGLE_PRECISION_COMPLEX | __KIND == __DOUBLE_PRECISION_COMPLEX
COMPLEX(MK) :: summ,summ2
COMPLEX(MK) :: eta,dm
#endif
! square of inter particle distance
REAL(MK) :: dij,dij2,dij4,dij5
! start and end particle in a box
INTEGER :: istart,iend,jstart,jend
! cell neighbor lists
INTEGER, DIMENSION(:,:), POINTER :: inp,jnp
! number of interactions for each cell
INTEGER, SAVE :: nnp
! cell offsets for box index
INTEGER :: n1,n2,nz
REAL(MK) :: t0
!-------------------------------------------------------------------------
! Externals
!-------------------------------------------------------------------------
#if __KERNEL == __USER_FUNCTION
#if __KIND == __SINGLE_PRECISION | __KIND == __DOUBLE_PRECISION
REAL(MK) , DIMENSION(:) , INTENT(IN ) :: kpar
INTERFACE
FUNCTION kernel(x,kpar)
IMPLICIT NONE
#if __KIND == __SINGLE_PRECISION | __KIND == __SINGLE_PRECISION_COMPLEX
INTEGER, PARAMETER :: MK = KIND(1.0E0)
#elif __KIND == __DOUBLE_PRECISION | __KIND == __DOUBLE_PRECISION_COMPLEX
INTEGER, PARAMETER :: MK = KIND(1.0D0)
#endif
REAL(MK), INTENT(IN) :: x
REAL(MK), DIMENSION(:), INTENT(IN) :: kpar
REAL(MK) :: kernel
END FUNCTION kernel
END INTERFACE
#elif __KIND == __SINGLE_PRECISION_COMPLEX | __KIND == __DOUBLE_PRECISION_COMPLEX
COMPLEX(MK), DIMENSION(:) , INTENT(IN ) :: kpar
INTERFACE
FUNCTION kernel(x,kpar)
IMPLICIT NONE
#if __KIND == __SINGLE_PRECISION | __KIND == __SINGLE_PRECISION_COMPLEX
INTEGER, PARAMETER :: MK = KIND(1.0E0)
#elif __KIND == __DOUBLE_PRECISION | __KIND == __DOUBLE_PRECISION_COMPLEX
INTEGER, PARAMETER :: MK = KIND(1.0D0)
#endif
REAL(MK), INTENT(IN) :: x
COMPLEX(MK), DIMENSION(:), INTENT(IN) :: kpar
COMPLEX(MK) :: kernel
END FUNCTION kernel
END INTERFACE
#endif
#endif
!-------------------------------------------------------------------------
! Initialise
!-------------------------------------------------------------------------
CALL substart('ppm_comp_pp_cell',t0,info)
!-------------------------------------------------------------------------
! Check arguments.
!-------------------------------------------------------------------------
IF (ppm_debug .GT. 0) THEN
IF (.NOT. ppm_initialized) THEN
info = ppm_error_error
CALL ppm_error(ppm_err_ppm_noinit,'ppm_comp_pp_cell', &
& 'Please call ppm_init first!',__LINE__,info)
GOTO 9999
ENDIF
IF (Np .LT. 0) THEN
info = ppm_error_error
CALL ppm_error(ppm_err_argument,'ppm_comp_pp_cell', &
& 'Np must be >=0',__LINE__,info)
GOTO 9999
ENDIF
IF (lda .LT. 1) THEN
info = ppm_error_error
CALL ppm_error(ppm_err_argument,'ppm_comp_pp_cell', &
& 'lda must be >0',__LINE__,info)
GOTO 9999
ENDIF
IF (nsublist .LT. 0) THEN
info = ppm_error_error
CALL ppm_error(ppm_err_argument,'ppm_comp_pp_cell', &
& 'nsublist must be >=0',__LINE__,info)
GOTO 9999
ENDIF
IF (cutoff2 .LT. 0.0_MK) THEN
info = ppm_error_error
CALL ppm_error(ppm_err_argument,'ppm_comp_pp_cell', &
& 'cutoff2 must be >=0',__LINE__,info)
GOTO 9999
ENDIF
#if __KERNEL == __LOOKUP_TABLE
IF (kpar .LT. 0.0_MK) THEN
info = ppm_error_error
CALL ppm_error(ppm_err_argument,'ppm_comp_pp_cell', &
& 'kpar (dxtableinv) must be >=0',__LINE__,info)
GOTO 9999
ENDIF
#endif
ENDIF
!-------------------------------------------------------------------------
! Build interaction index lists
!-------------------------------------------------------------------------
CALL ppm_neighlist_MkNeighIdx(lsymm,inp,jnp,nnp,info)
!-------------------------------------------------------------------------
! PARTICLE-PARTICLE INTERACTIONS using symmetry
!-------------------------------------------------------------------------
IF (lsymm) THEN
DO idom=1,nsublist
n1 = Nm(1,idom)
n2 = Nm(1,idom)*Nm(2,idom)
nz = Nm(3,idom)
IF (ppm_dim .EQ. 2) THEN
n2 = 0
nz = 2
ENDIF
! loop over all REAL cells (the -2 in the end does this)
DO k=0,nz-2
DO j=0,Nm(2,idom)-2
DO i=0,Nm(1,idom)-2
! index of the center box
cbox = i + 1 + n1*j + n2*k
! loop over all box-box interactions
DO iinter=1,nnp
! determine box indices for this interaction
ibox = cbox+(inp(1,iinter)+n1*inp(2,iinter)+ &
& n2*inp(3,iinter))
jbox = cbox+(jnp(1,iinter)+n1*jnp(2,iinter)+ &
& n2*jnp(3,iinter))
!-------------------------------------------------
! Read indices and check if box is empty
!-------------------------------------------------
istart = clist(idom)%lhbx(ibox)
iend = clist(idom)%lhbx(ibox+1)-1
IF (iend .LT. istart) CYCLE
!-------------------------------------------------
! Within the box itself use symmetry and avoid
! adding the particle itself to its own list
!-------------------------------------------------
IF (ibox .EQ. jbox) THEN
DO ipart=istart,iend
ip = clist(idom)%lpdx(ipart)
IF (lda .EQ. 1) THEN
summ = 0.0_MK
#ifdef __SXF90
!CDIR NODEP
#endif
DO jpart=(ipart+1),iend
jp = clist(idom)%lpdx(jpart)
dx = xp(1,ip) - xp(1,jp)
dy = xp(2,ip) - xp(2,jp)
IF (ppm_dim .GT. 2) THEN
dz = xp(3,ip) - xp(3,jp)
dij = (dx*dx)+(dy*dy)+(dz*dz)
ELSE
dz = 0.0_MK
dij = (dx*dx)+(dy*dy)
ENDIF
IF (dij .GT. cutoff2) CYCLE
!---------------------------------
! Particle ip interacts with
! particle jp here... and
! vice versa to use symmetry.
!---------------------------------
#include "ppm_comp_pp_kernels.inc"
dm = eta*(pdata(1,jp) - &
& pdata(1,ip))
summ = summ + dm
dpd(1,jp) = dpd(1,jp) - dm
ENDDO
dpd(1,ip) = dpd(1,ip) + summ
ELSEIF (lda .EQ. 2) THEN
summ = 0.0_MK
summ2 = 0.0_MK
#ifdef __SXF90
!CDIR NODEP
#endif
DO jpart=(ipart+1),iend
jp = clist(idom)%lpdx(jpart)
dx = xp(1,ip) - xp(1,jp)
dy = xp(2,ip) - xp(2,jp)
IF (ppm_dim .GT. 2) THEN
dz = xp(3,ip) - xp(3,jp)
dij = (dx*dx)+(dy*dy)+(dz*dz)
ELSE
dz = 0.0_MK
dij = (dx*dx)+(dy*dy)
ENDIF
IF (dij .GT. cutoff2) CYCLE
!---------------------------------
! Particle ip interacts with
! particle jp here... and
! vice versa to use symmetry.
!---------------------------------
#include "ppm_comp_pp_kernels.inc"
dm = eta*(pdata(1,jp) - &
& pdata(1,ip))
summ = summ + dm
dpd(1,jp) = dpd(1,jp) - dm
dm = eta*(pdata(2,jp) - &
& pdata(2,ip))
summ2 = summ2 + dm
dpd(2,jp) = dpd(2,jp) - dm
ENDDO
dpd(1,ip) = dpd(1,ip) + summ
dpd(2,ip) = dpd(2,ip) + summ2
ELSE
DO jpart=(ipart+1),iend
jp = clist(idom)%lpdx(jpart)
dx = xp(1,ip) - xp(1,jp)
dy = xp(2,ip) - xp(2,jp)
IF (ppm_dim .GT. 2) THEN
dz = xp(3,ip) - xp(3,jp)
dij = (dx*dx)+(dy*dy)+(dz*dz)
ELSE
dz = 0.0_MK
dij = (dx*dx)+(dy*dy)
ENDIF
IF (dij .GT. cutoff2) CYCLE
!---------------------------------
! Particle ip interacts with
! particle jp here... and
! vice versa to use symmetry.
!---------------------------------
#include "ppm_comp_pp_kernels.inc"
DO ispec=1,lda
dm = eta*(pdata(ispec,jp)- &
& pdata(ispec,ip))
dpd(ispec,ip)=dpd(ispec,ip)+dm
dpd(ispec,jp)=dpd(ispec,jp)-dm
ENDDO
ENDDO
ENDIF
ENDDO
!-------------------------------------------------
! For the other boxes check all particles
!-------------------------------------------------
ELSE
! get pointers to first and last particle
jstart = clist(idom)%lhbx(jbox)
jend = clist(idom)%lhbx(jbox+1)-1
! skip this iinter if other box is empty
IF (jend .LT. jstart) CYCLE
! loop over all particles inside this cell
DO ipart=istart,iend
ip = clist(idom)%lpdx(ipart)
IF (lda .EQ. 1) THEN
summ = 0.0_MK
! check against all particles
! in the other cell
#ifdef __SXF90
!CDIR NODEP
#endif
DO jpart=jstart,jend
jp = clist(idom)%lpdx(jpart)
dx = xp(1,ip) - xp(1,jp)
dy = xp(2,ip) - xp(2,jp)
IF (ppm_dim .GT. 2) THEN
dz = xp(3,ip) - xp(3,jp)
dij = (dx*dx)+(dy*dy)+(dz*dz)
ELSE
dz = 0.0_MK
dij = (dx*dx)+(dy*dy)
ENDIF
IF (dij .GT. cutoff2) CYCLE
!---------------------------------
! Particle ip interacts with
! particle jp here... and
! vice versa to use symmetry.
!---------------------------------
#include "ppm_comp_pp_kernels.inc"
dm = eta*(pdata(1,jp) - &
& pdata(1,ip))
summ = summ +dm
dpd(1,jp) = dpd(1,jp) - dm
ENDDO
dpd(1,ip) = dpd(1,ip) + summ
ELSEIF (lda .EQ. 2) THEN
summ = 0.0_MK
summ2 = 0.0_MK
! check against all particles
! in the other cell
#ifdef __SXF90
!CDIR NODEP
#endif
DO jpart=jstart,jend
jp = clist(idom)%lpdx(jpart)
dx = xp(1,ip) - xp(1,jp)
dy = xp(2,ip) - xp(2,jp)
IF (ppm_dim .GT. 2) THEN
dz = xp(3,ip) - xp(3,jp)
dij = (dx*dx)+(dy*dy)+(dz*dz)
ELSE
dz = 0.0_MK
dij = (dx*dx)+(dy*dy)
ENDIF
IF (dij .GT. cutoff2) CYCLE
!---------------------------------
! Particle ip interacts with
! particle jp here... and
! vice versa to use symmetry.
!---------------------------------
#include "ppm_comp_pp_kernels.inc"
dm = eta*(pdata(1,jp) - &
& pdata(1,ip))
summ = summ + dm
dpd(1,jp) = dpd(1,jp) - dm
dm = eta*(pdata(2,jp) - &
& pdata(2,ip))
summ2 = summ2 + dm
dpd(2,jp) = dpd(2,jp) - dm
ENDDO
dpd(1,ip) = dpd(1,ip) + summ
dpd(2,ip) = dpd(2,ip) + summ2
ELSE
! check against all particles
! in the other cell
DO jpart=jstart,jend
jp = clist(idom)%lpdx(jpart)
dx = xp(1,ip) - xp(1,jp)
dy = xp(2,ip) - xp(2,jp)
IF (ppm_dim .GT. 2) THEN
dz = xp(3,ip) - xp(3,jp)
dij = (dx*dx)+(dy*dy)+(dz*dz)
ELSE
dz = 0.0_MK
dij = (dx*dx)+(dy*dy)
ENDIF
IF (dij .GT. cutoff2) CYCLE
!---------------------------------
! Particle ip interacts with
! particle jp here... and
! vice versa to use symmetry.
!---------------------------------
#include "ppm_comp_pp_kernels.inc"
DO ispec=1,lda
dm = eta*(pdata(ispec,jp)- &
& pdata(ispec,ip))
dpd(ispec,ip)=dpd(ispec,ip)+dm
dpd(ispec,jp)=dpd(ispec,jp)-dm
ENDDO
ENDDO
ENDIF
ENDDO
ENDIF ! ibox .EQ. jbox
ENDDO ! iinter
ENDDO ! i
ENDDO ! j
ENDDO ! k
ENDDO ! idom
!-------------------------------------------------------------------------
! PARTICLE-PARTICLE INTERACTIONS not using symmetry
!-------------------------------------------------------------------------
ELSE
DO idom=1,nsublist
n1 = Nm(1,idom)
n2 = Nm(1,idom)*Nm(2,idom)
nz = Nm(3,idom)
IF (ppm_dim .EQ. 2) THEN
n2 = 0
nz = 2
ENDIF
! loop over all REAL cells (the -2 in the end does this)
DO k=1,nz-2
DO j=1,Nm(2,idom)-2
DO i=1,Nm(1,idom)-2
! index of the center box
cbox = i + 1 + n1*j + n2*k
! loop over all box-box interactions
DO iinter=1,nnp
! determine box indices for this interaction
ibox = cbox+(inp(1,iinter)+n1*inp(2,iinter)+ &
& n2*inp(3,iinter))
jbox = cbox+(jnp(1,iinter)+n1*jnp(2,iinter)+ &
& n2*jnp(3,iinter))
!-------------------------------------------------
! Read indices and check if box is empty
!-------------------------------------------------
istart = clist(idom)%lhbx(ibox)
iend = clist(idom)%lhbx(ibox+1)-1
IF (iend .LT. istart) CYCLE
!-------------------------------------------------
! Do all interactions within the box itself
!-------------------------------------------------
IF (ibox .EQ. jbox) THEN
DO ipart=istart,iend
ip = clist(idom)%lpdx(ipart)
IF (lda .EQ. 1) THEN
DO jpart=istart,iend
jp = clist(idom)%lpdx(jpart)
! No particle interacts with
! itself
IF (ip .EQ. jp) CYCLE
dx = xp(1,ip) - xp(1,jp)
dy = xp(2,ip) - xp(2,jp)
IF (ppm_dim .GT. 2) THEN
dz = xp(3,ip) - xp(3,jp)
dij = (dx*dx)+(dy*dy)+(dz*dz)
ELSE
dz = 0.0_MK
dij = (dx*dx)+(dy*dy)
ENDIF
IF (dij .GT. cutoff2) CYCLE
!---------------------------------
! Particle ip interacts with
! particle jp.
!---------------------------------
#include "ppm_comp_pp_kernels.inc"
dm = eta*(pdata(1,jp) - &
& pdata(1,ip))
dpd(1,ip) = dpd(1,ip) + dm
ENDDO
ELSEIF (lda .EQ. 2) THEN
DO jpart=istart,iend
jp = clist(idom)%lpdx(jpart)
! No particle interacts with
! itself
IF (ip .EQ. jp) CYCLE
dx = xp(1,ip) - xp(1,jp)
dy = xp(2,ip) - xp(2,jp)
IF (ppm_dim .GT. 2) THEN
dz = xp(3,ip) - xp(3,jp)
dij = (dx*dx)+(dy*dy)+(dz*dz)
ELSE
dz = 0.0_MK
dij = (dx*dx)+(dy*dy)
ENDIF
IF (dij .GT. cutoff2) CYCLE
!---------------------------------
! Particle ip interacts with
! particle jp.
!---------------------------------
#include "ppm_comp_pp_kernels.inc"
dm = eta*(pdata(1,jp) - &
& pdata(1,ip))
dpd(1,ip) = dpd(1,ip) + dm
dm = eta*(pdata(2,jp) - &
& pdata(2,ip))
dpd(2,ip) = dpd(2,ip) + dm
ENDDO
ELSE
DO jpart=istart,iend
jp = clist(idom)%lpdx(jpart)
! No particle interacts with
! itself
IF (ip .EQ. jp) CYCLE
dx = xp(1,ip) - xp(1,jp)
dy = xp(2,ip) - xp(2,jp)
IF (ppm_dim .GT. 2) THEN
dz = xp(3,ip) - xp(3,jp)
dij = (dx*dx)+(dy*dy)+(dz*dz)
ELSE
dz = 0.0_MK
dij = (dx*dx)+(dy*dy)
ENDIF
IF (dij .GT. cutoff2) CYCLE
!---------------------------------
! Particle ip interacts with
! particle jp.
!---------------------------------
#include "ppm_comp_pp_kernels.inc"
DO ispec=1,lda
dm = eta*(pdata(ispec,jp)- &
& pdata(ispec,ip))
dpd(ispec,ip)=dpd(ispec,ip)+dm
ENDDO
ENDDO
ENDIF
ENDDO
!-------------------------------------------------
! Do interactions with all neighboring boxes
!-------------------------------------------------
ELSE
! get pointers to first and last particle
jstart = clist(idom)%lhbx(jbox)
jend = clist(idom)%lhbx(jbox+1)-1
! skip if empty
IF (jend .LT. jstart) CYCLE
! loop over all particles inside this cell
DO ipart=istart,iend
ip = clist(idom)%lpdx(ipart)
IF (lda .EQ. 1) THEN
! check against all particles
! in the other cell
DO jpart=jstart,jend
jp = clist(idom)%lpdx(jpart)
dx = xp(1,ip) - xp(1,jp)
dy = xp(2,ip) - xp(2,jp)
IF (ppm_dim .GT. 2) THEN
dz = xp(3,ip) - xp(3,jp)
dij = (dx*dx)+(dy*dy)+(dz*dz)
ELSE
dz = 0.0_MK
dij = (dx*dx)+(dy*dy)
ENDIF
IF (dij .GT. cutoff2) CYCLE
!---------------------------------
! Particle ip interacts with
! particle jp.
!---------------------------------
#include "ppm_comp_pp_kernels.inc"
dm = eta*(pdata(1,jp) - &
& pdata(1,ip))
dpd(1,ip) = dpd(1,ip) + dm
ENDDO
ELSEIF (lda .EQ. 2) THEN
! check against all particles
! in the other cell
DO jpart=jstart,jend
jp = clist(idom)%lpdx(jpart)
dx = xp(1,ip) - xp(1,jp)
dy = xp(2,ip) - xp(2,jp)
IF (ppm_dim .GT. 2) THEN
dz = xp(3,ip) - xp(3,jp)
dij = (dx*dx)+(dy*dy)+(dz*dz)
ELSE
dz = 0.0_MK
dij = (dx*dx)+(dy*dy)
ENDIF
IF (dij .GT. cutoff2) CYCLE
!---------------------------------
! Particle ip interacts with
! particle jp.
!---------------------------------
#include "ppm_comp_pp_kernels.inc"
dm = eta*(pdata(1,jp) - &
& pdata(1,ip))
dpd(1,ip) = dpd(1,ip) + dm
dm = eta*(pdata(2,jp) - &
& pdata(2,ip))
dpd(2,ip) = dpd(2,ip) + dm
ENDDO
ELSE
! check against all particles
! in the other cell
DO jpart=jstart,jend
jp = clist(idom)%lpdx(jpart)
dx = xp(1,ip) - xp(1,jp)
dy = xp(2,ip) - xp(2,jp)
IF (ppm_dim .GT. 2) THEN
dz = xp(3,ip) - xp(3,jp)
dij = (dx*dx)+(dy*dy)+(dz*dz)
ELSE
dz = 0.0_MK
dij = (dx*dx)+(dy*dy)
ENDIF
IF (dij .GT. cutoff2) CYCLE
!---------------------------------
! Particle ip interacts with
! particle jp.
!---------------------------------
#include "ppm_comp_pp_kernels.inc"
DO ispec=1,lda
dm = eta*(pdata(ispec,jp)- &
& pdata(ispec,ip))
dpd(ispec,ip)=dpd(ispec,ip)+dm
ENDDO
ENDDO
ENDIF
ENDDO
ENDIF ! ibox .EQ. jbox
ENDDO ! iinter
ENDDO ! i
ENDDO ! j
ENDDO ! k
ENDDO ! idom
ENDIF
!-------------------------------------------------------------------------
! Return
!-------------------------------------------------------------------------
9999 CONTINUE
CALL substop('ppm_comp_pp_cell',t0,info)
RETURN
#if __KERNEL == __INTERNAL
#if __KIND == __SINGLE_PRECISION
END SUBROUTINE ppm_comp_pp_cell_si
#elif __KIND == __DOUBLE_PRECISION
END SUBROUTINE ppm_comp_pp_cell_di
#elif __KIND == __SINGLE_PRECISION_COMPLEX
END SUBROUTINE ppm_comp_pp_cell_sci
#elif __KIND == __DOUBLE_PRECISION_COMPLEX
END SUBROUTINE ppm_comp_pp_cell_dci
#endif
#elif __KERNEL == __USER_FUNCTION
#if __KIND == __SINGLE_PRECISION
END SUBROUTINE ppm_comp_pp_cell_su
#elif __KIND == __DOUBLE_PRECISION
END SUBROUTINE ppm_comp_pp_cell_du
#elif __KIND == __SINGLE_PRECISION_COMPLEX
END SUBROUTINE ppm_comp_pp_cell_scu
#elif __KIND == __DOUBLE_PRECISION_COMPLEX
END SUBROUTINE ppm_comp_pp_cell_dcu
#endif
#elif __KERNEL == __LOOKUP_TABLE
#if __KIND == __SINGLE_PRECISION
END SUBROUTINE ppm_comp_pp_cell_st
#elif __KIND == __DOUBLE_PRECISION
END SUBROUTINE ppm_comp_pp_cell_dt
#elif __KIND == __SINGLE_PRECISION_COMPLEX
END SUBROUTINE ppm_comp_pp_cell_sct
#elif __KIND == __DOUBLE_PRECISION_COMPLEX
END SUBROUTINE ppm_comp_pp_cell_dct
#endif
#endif
src/ppm_comp_pp_correct.f 0 → 100644
+ 343
− 0
View file @ dd7b0757
!-------------------------------------------------------------------------
! Subroutine : ppm_comp_pp_correct
!-------------------------------------------------------------------------
! Parallel Particle Mesh Library (PPM)
! ETH Zurich
! CH-8092 Zurich, Switzerland
!-------------------------------------------------------------------------
#if __KERNEL == __INTERNAL
#if __KIND == __SINGLE_PRECISION
SUBROUTINE ppm_comp_pp_correct_si(dh,cutoff,kernel,eps,kpar,targt, &
& kappa,info)
#elif __KIND == __DOUBLE_PRECISION
SUBROUTINE ppm_comp_pp_correct_di(dh,cutoff,kernel,eps,kpar,targt, &
& kappa,info)
#elif __KIND == __SINGLE_PRECISION_COMPLEX
SUBROUTINE ppm_comp_pp_correct_sci(dh,cutoff,kernel,eps,kpar,targt, &
& kappa,info)
#elif __KIND == __DOUBLE_PRECISION_COMPLEX
SUBROUTINE ppm_comp_pp_correct_dci(dh,cutoff,kernel,eps,kpar,targt, &
& kappa,info)
#endif
#elif __KERNEL == __USER_FUNCTION
#if __KIND == __SINGLE_PRECISION
SUBROUTINE ppm_comp_pp_correct_su(dh,cutoff,kernel,eps,kpar,targt, &
& kappa,info)
#elif __KIND == __DOUBLE_PRECISION
SUBROUTINE ppm_comp_pp_correct_du(dh,cutoff,kernel,eps,kpar,targt, &
& kappa,info)
#elif __KIND == __SINGLE_PRECISION_COMPLEX
SUBROUTINE ppm_comp_pp_correct_scu(dh,cutoff,kernel,eps,kpar,targt, &
& kappa,info)
#elif __KIND == __DOUBLE_PRECISION_COMPLEX
SUBROUTINE ppm_comp_pp_correct_dcu(dh,cutoff,kernel,eps,kpar,targt, &
& kappa,info)
#endif
#elif __KERNEL == __LOOKUP_TABLE
#if __KIND == __SINGLE_PRECISION
SUBROUTINE ppm_comp_pp_correct_st(dh,cutoff,kernel,eps,kpar,targt, &
& kappa,info)
#elif __KIND == __DOUBLE_PRECISION
SUBROUTINE ppm_comp_pp_correct_dt(dh,cutoff,kernel,eps,kpar,targt, &
& kappa,info)
#elif __KIND == __SINGLE_PRECISION_COMPLEX
SUBROUTINE ppm_comp_pp_correct_sct(dh,cutoff,kernel,eps,kpar,targt, &
& kappa,info)
#elif __KIND == __DOUBLE_PRECISION_COMPLEX
SUBROUTINE ppm_comp_pp_correct_dct(dh,cutoff,kernel,eps,kpar,targt, &
& kappa,info)
#endif
#endif
!!! Computes the discretisation correction factor
!!! needed to ensure that the discrete second-order
!!! moment of the kernel function is equal to some
!!! value. The kernel is evaluated on a regular mesh
!!! here, so this factor can only be used if the
!!! particles are on a regular grid (e.g. after
!!! remeshing).
!-------------------------------------------------------------------------
! Modules
!-------------------------------------------------------------------------
USE ppm_module_data
USE ppm_module_substart
USE ppm_module_substop
USE ppm_module_error
USE ppm_module_write
IMPLICIT NONE
#if __KIND == __SINGLE_PRECISION | __KIND == __SINGLE_PRECISION_COMPLEX
INTEGER, PARAMETER :: MK = ppm_kind_single
#elif __KIND == __DOUBLE_PRECISION | __KIND == __DOUBLE_PRECISION_COMPLEX
INTEGER, PARAMETER :: MK = ppm_kind_double
#endif
!-------------------------------------------------------------------------
! Includes
!-------------------------------------------------------------------------
#include "ppm_define.h"
!-------------------------------------------------------------------------
! Arguments
!-------------------------------------------------------------------------
REAL(MK), DIMENSION(:), INTENT(IN ) :: dh
!!! grid spacings in all directions 1..ppm_dim
REAL(MK) , INTENT(IN ) :: cutoff
!!! PP interaction cutoff used. If no cutoff is used, pass a negative
!!! value here.
REAL(MK) , INTENT(IN ) :: eps
!!! Kernel size (width)
#if __KIND == __SINGLE_PRECISION | __KIND == __DOUBLE_PRECISION
REAL(MK) , INTENT(IN ) :: targt
!!! target value for the second order kernel moment. Overloaded types:
!!! single,double
REAL(MK) , INTENT( OUT) :: kappa
!!! computed discretisation correction factor. Every kernel evaluation has
!!! to be multiplied with this in order to correct for discretisation
!!! effects. Overloaded types: single,double
#elif __KIND == __SINGLE_PRECISION_COMPLEX | __KIND == __DOUBLE_PRECISION_COMPLEX
COMPLEX(MK) , INTENT(IN ) :: targt
!!! target value for the second order kernel moment. Overloaded types:
!!! single complex, double complex
COMPLEX(MK) , INTENT( OUT) :: kappa
!!! computed discretisation correction factor. Every kernel evaluation has
!!! to be multiplied with this in order to correct for discretisation
!!! effects. Overloaded types: single complex, double complex
#endif
#if __KERNEL == __INTERNAL
INTEGER , INTENT(IN ) :: kernel
!!! kernel for which to compute the correction. To use ppm-internal
!!! kernels, specify one of:
!!!
!!! ---------------------------------------
!!! ppm_param_kerel_laplace2d_2p
!!! (2nd order Laplacian,
!!! polynomial in 2D)
!!! ppm_param_kerel_laplace3d_2p
!!! (2nd order Laplacian,
!!! polynomial in 3D)
!!! ---------------------------------------
!!!
!!! To use your own kernel function, pass the function pointer here. Your
!!! function should take one argument and return one value.
#if __KIND == __SINGLE_PRECISION | __KIND == __DOUBLE_PRECISION
REAL(MK) , DIMENSION(:), INTENT(IN ) :: kpar
!!! Kernel parameters. See documentation or ppm_comp_pp_kernels.inc for
!!! description. Type can be single or complex.
#elif __KIND == __SINGLE_PRECISION_COMPLEX | __KIND == __DOUBLE_PRECISION_COMPLEX
COMPLEX(MK), DIMENSION(:), INTENT(IN ) :: kpar
!!! Kernel parameters. See documentation or ppm_comp_pp_kernels.inc for
!!! description. Type can be single complex or double complex.
#endif
#elif __KERNEL == __LOOKUP_TABLE
REAL(MK) , INTENT(IN ) :: kpar
!!! Kernel parameters. See documentation or ppm_comp_pp_kernels.inc for
!!! description. Pass dxtableinv (the inverse of the table spacing) as
!!! a scalar here.
#if __KIND == __SINGLE_PRECISION | __KIND == __DOUBLE_PRECISION
REAL(MK) , DIMENSION(: ), INTENT(IN ) :: kernel
#elif __KIND == __SINGLE_PRECISION_COMPLEX | __KIND == __DOUBLE_PRECISION_COMPLEX
COMPLEX(MK), DIMENSION(: ), INTENT(IN ) :: kernel
#endif
!!! Lookup table with tabulated kernel values.
!!! Such a table can be created using <<ppm_comp_pp_mk_table>>.
#endif
INTEGER , INTENT( OUT) :: info
!!! Returns 0 upon success
!-------------------------------------------------------------------------
! Local variables
!-------------------------------------------------------------------------
REAL(MK) :: t0,dij2,dij4,dij5
REAL(MK) :: hi,dx,dy,dz,dij,factor,factor2
#if __KIND == __SINGLE_PRECISION | __KIND == __DOUBLE_PRECISION
REAL(MK) :: eta,summ
#elif __KIND == __SINGLE_PRECISION_COMPLEX | __KIND == __DOUBLE_PRECISION_COMPLEX
COMPLEX(MK) :: eta,summ
#endif
INTEGER :: lhx,lhy,lhz,i,j,k,idx
CHARACTER(LEN=ppm_char) :: cbuf
!-------------------------------------------------------------------------
! Externals
!-------------------------------------------------------------------------
#if __KERNEL == __USER_FUNCTION
#if __KIND == __SINGLE_PRECISION | __KIND == __DOUBLE_PRECISION
REAL(MK) , DIMENSION(:), INTENT(IN ) :: kpar
INTERFACE
FUNCTION kernel(x,kpar)
IMPLICIT NONE
#if __KIND == __SINGLE_PRECISION | __KIND == __SINGLE_PRECISION_COMPLEX
INTEGER, PARAMETER :: MK = KIND(1.0E0)
#elif __KIND == __DOUBLE_PRECISION | __KIND == __DOUBLE_PRECISION_COMPLEX
INTEGER, PARAMETER :: MK = KIND(1.0D0)
#endif
REAL(MK), INTENT(IN) :: x
REAL(MK), DIMENSION(:), INTENT(IN) :: kpar
REAL(MK) :: kernel
END FUNCTION kernel
END INTERFACE
#elif __KIND == __SINGLE_PRECISION_COMPLEX | __KIND == __DOUBLE_PRECISION_COMPLEX
COMPLEX(MK), DIMENSION(:), INTENT(IN ) :: kpar
INTERFACE
FUNCTION kernel(x,kpar)
IMPLICIT NONE
#if __KIND == __SINGLE_PRECISION | __KIND == __SINGLE_PRECISION_COMPLEX
INTEGER, PARAMETER :: MK = KIND(1.0E0)
#elif __KIND == __DOUBLE_PRECISION | __KIND == __DOUBLE_PRECISION_COMPLEX
INTEGER, PARAMETER :: MK = KIND(1.0D0)
#endif
REAL(MK), INTENT(IN) :: x
COMPLEX(MK), DIMENSION(:), INTENT(IN) :: kpar
COMPLEX(MK) :: kernel
END FUNCTION kernel
END INTERFACE
#endif
#endif
!-------------------------------------------------------------------------
! Initialise
!-------------------------------------------------------------------------
CALL substart('ppm_comp_pp_correct',t0,info)
kappa = 1.0_MK
!-------------------------------------------------------------------------
! Check arguments.
!-------------------------------------------------------------------------
IF (ppm_debug .GT. 0) THEN
IF (.NOT. ppm_initialized) THEN
info = ppm_error_error
CALL ppm_error(ppm_err_ppm_noinit,'ppm_comp_pp_correct', &
& 'Please call ppm_init first!',__LINE__,info)
GOTO 9999
ENDIF
IF (eps .LT. 0.0_MK) THEN
info = ppm_error_error
CALL ppm_error(ppm_err_argument,'ppm_comp_pp_correct', &
& 'eps must be > 0',__LINE__,info)
GOTO 9999
ENDIF
IF (SIZE(dh,1) .LT. ppm_dim) THEN
info = ppm_error_error
CALL ppm_error(ppm_err_argument,'ppm_comp_pp_correct', &
& 'dh must be at least of length ppm_dim',__LINE__,info)
GOTO 9999
ENDIF
DO i=1,ppm_dim
IF (dh(i) .LE. 0.0_MK) THEN
info = ppm_error_error
CALL ppm_error(ppm_err_argument,'ppm_comp_pp_correct', &
& 'dh must be >0 in all directions',__LINE__,info)
GOTO 9999
ENDIF
ENDDO
#if __KERNEL == __LOOKUP_TABLE
IF (kpar .LT. 0.0_MK) THEN
info = ppm_error_error
CALL ppm_error(ppm_err_argument,'ppm_comp_pp_correct', &
& 'kpar (dxtableinv) must be >=0',__LINE__,info)
GOTO 9999
ENDIF
#endif
ENDIF
!-------------------------------------------------------------------------
! Integration bounds
!-------------------------------------------------------------------------
IF (cutoff .GE. 0.0_MK) THEN
hi = cutoff
ELSE
hi = 50.0_MK
ENDIF
!-------------------------------------------------------------------------
! Summation bounds
!-------------------------------------------------------------------------
lhx = INT(hi/dh(1))
lhy = INT(hi/dh(2))
IF (ppm_dim .GT. 2) lhz = INT(hi/dh(3))
!-------------------------------------------------------------------------
! Discrete sum
!-------------------------------------------------------------------------
summ = 0.0_MK
IF (ppm_dim .GT. 2) THEN
DO i=-lhx,lhx
dx = REAL(i,MK)*dh(1)
DO j=-lhy,lhy
dy = REAL(j,MK)*dh(2)
DO k=-lhz,lhz
dz = REAL(k,MK)*dh(3)
dij = dx*dx + dy*dy + dz*dz
#include "ppm_comp_pp_kernels.inc"
! second order moment
eta = dij*eta
summ = summ + eta
ENDDO
ENDDO
ENDDO
summ = summ*dh(1)*dh(2)*dh(3)
ELSE
DO i=-lhx,lhx
dx = REAL(i,MK)*dh(1)
DO j=-lhy,lhy
dy = REAL(j,MK)*dh(2)
dij = dx*dx + dy*dy
#include "ppm_comp_pp_kernels.inc"
! second order moment
eta = dij*eta
summ = summ + eta
ENDDO
ENDDO
summ = summ*dh(1)*dh(2)
ENDIF
!-------------------------------------------------------------------------
! Correction factor
!-------------------------------------------------------------------------
#if __KERNEL != __LOOKUP_TABLE
kappa = (REAL(ppm_dim,MK)*targt*eps)/summ
#endif
!-------------------------------------------------------------------------
! Debug output
!-------------------------------------------------------------------------
IF (ppm_debug .GT. 1) THEN
WRITE(cbuf,'(A,F16.10)') 'Discretization correction factor: ',kappa
CALL ppm_write(ppm_rank,'ppm_comp_pp_correct',cbuf,info)
ENDIF
!-------------------------------------------------------------------------
! Return
!-------------------------------------------------------------------------
9999 CONTINUE
CALL substop('ppm_comp_pp_correct',t0,info)
RETURN
#if __KERNEL == __INTERNAL
#if __KIND == __SINGLE_PRECISION
END SUBROUTINE ppm_comp_pp_correct_si
#elif __KIND == __DOUBLE_PRECISION
END SUBROUTINE ppm_comp_pp_correct_di
#elif __KIND == __SINGLE_PRECISION_COMPLEX
END SUBROUTINE ppm_comp_pp_correct_sci
#elif __KIND == __DOUBLE_PRECISION_COMPLEX
END SUBROUTINE ppm_comp_pp_correct_dci
#endif
#elif __KERNEL == __USER_FUNCTION
#if __KIND == __SINGLE_PRECISION
END SUBROUTINE ppm_comp_pp_correct_su
#elif __KIND == __DOUBLE_PRECISION
END SUBROUTINE ppm_comp_pp_correct_du
#elif __KIND == __SINGLE_PRECISION_COMPLEX
END SUBROUTINE ppm_comp_pp_correct_scu
#elif __KIND == __DOUBLE_PRECISION_COMPLEX
END SUBROUTINE ppm_comp_pp_correct_dcu
#endif
#elif __KERNEL == __LOOKUP_TABLE
#if __KIND == __SINGLE_PRECISION
END SUBROUTINE ppm_comp_pp_correct_st
#elif __KIND == __DOUBLE_PRECISION
END SUBROUTINE ppm_comp_pp_correct_dt
#elif __KIND == __SINGLE_PRECISION_COMPLEX
END SUBROUTINE ppm_comp_pp_correct_sct
#elif __KIND == __DOUBLE_PRECISION_COMPLEX
END SUBROUTINE ppm_comp_pp_correct_dct
#endif
#endif
src/ppm_comp_pp_doring.f 0 → 100644
+ 611
− 0
View file @ dd7b0757
This diff is collapsed.
src/ppm_comp_pp_mk_table.f 0 → 100644
+ 247
− 0
View file @ dd7b0757
!-------------------------------------------------------------------------
! Subroutine : ppm_comp_pp_mk_table
!-------------------------------------------------------------------------
! Parallel Particle Mesh Library (PPM)
! ETH Zurich
! CH-8092 Zurich, Switzerland
!-------------------------------------------------------------------------
#if __KERNEL == __INTERNAL
#if __KIND == __SINGLE_PRECISION
SUBROUTINE ppm_comp_pp_mk_table_si(kernel,kpar,ntable,cutoff2, &
& ktab,dxtableinv,info)
#elif __KIND == __DOUBLE_PRECISION
SUBROUTINE ppm_comp_pp_mk_table_di(kernel,kpar,ntable,cutoff2, &
& ktab,dxtableinv,info)
#elif __KIND == __SINGLE_PRECISION_COMPLEX
SUBROUTINE ppm_comp_pp_mk_table_sci(kernel,kpar,ntable,cutoff2, &
& ktab,dxtableinv,info)
#elif __KIND == __DOUBLE_PRECISION_COMPLEX
SUBROUTINE ppm_comp_pp_mk_table_dci(kernel,kpar,ntable,cutoff2, &
& ktab,dxtableinv,info)
#endif
#elif __KERNEL == __USER_FUNCTION
#if __KIND == __SINGLE_PRECISION
SUBROUTINE ppm_comp_pp_mk_table_su(kernel,kpar,ntable,cutoff2, &
& ktab,dxtableinv,info)
#elif __KIND == __DOUBLE_PRECISION
SUBROUTINE ppm_comp_pp_mk_table_du(kernel,kpar,ntable,cutoff2, &
& ktab,dxtableinv,info)
#elif __KIND == __SINGLE_PRECISION_COMPLEX
SUBROUTINE ppm_comp_pp_mk_table_scu(kernel,kpar,ntable,cutoff2, &
& ktab,dxtableinv,info)
#elif __KIND == __DOUBLE_PRECISION_COMPLEX
SUBROUTINE ppm_comp_pp_mk_table_dcu(kernel,kpar,ntable,cutoff2, &
& ktab,dxtableinv,info)
#endif
#endif
!!! Creates a lookup table for a PP interaction kernel.
!!! The result is returned in the variables ktab with
!!! indices running from 0 to ntable-1. The table
!!! contains kernel values eta as a function of the
!!! _squared_ distance x**2.
!-------------------------------------------------------------------------
! Modules
!-------------------------------------------------------------------------
USE ppm_module_data
USE ppm_module_substart
USE ppm_module_substop
USE ppm_module_error
USE ppm_module_alloc
IMPLICIT NONE
#if __KIND == __SINGLE_PRECISION | __KIND == __SINGLE_PRECISION_COMPLEX
INTEGER, PARAMETER :: MK = ppm_kind_single
#elif __KIND == __DOUBLE_PRECISION | __KIND == __DOUBLE_PRECISION_COMPLEX
INTEGER, PARAMETER :: MK = ppm_kind_double
#endif
!-------------------------------------------------------------------------
! Includes
!-------------------------------------------------------------------------
#include "ppm_define.h"
!-------------------------------------------------------------------------
! Arguments
!-------------------------------------------------------------------------
INTEGER , INTENT(IN ) :: ntable
!!! Number of table entries to be created.
REAL(MK) , INTENT(IN ) :: cutoff2
!!! Square of the cutoff used for PP interactions. The kernel will be
!!! tabulated up to this value.
#if __KIND == __SINGLE_PRECISION | __KIND == __DOUBLE_PRECISION
REAL(MK) , DIMENSION(:), POINTER :: ktab
!!! Tabulated kernel values 0...ntable-1.
!!! Overloaded types: single,double.
#elif __KIND == __SINGLE_PRECISION_COMPLEX | __KIND == __DOUBLE_PRECISION_COMPLEX
COMPLEX(MK) , DIMENSION(:), POINTER :: ktab
!!! Tabulated kernel values 0...ntable-1.
!!! Overloaded types: single complex,double complex.
#endif
#if __KERNEL == __INTERNAL
INTEGER , INTENT(IN ) :: kernel
!!! kernel for which to compute the correction. To use ppm-internal
!!! kernels, specify one of:
!!!
!!! ---------------------------------------
!!! ppm_param_kerel_laplace2d_2p
!!! (2nd order Laplacian,
!!! polynomial in 2D)
!!! ppm_param_kerel_laplace3d_2p
!!! (2nd order Laplacian,
!!! polynomial in 3D)
!!! ---------------------------------------
!!!
!!! To use your own kernel function, pass the function pointer here. Your
!!! function should take one argument and return one value.
#if __KIND == __SINGLE_PRECISION | __KIND == __DOUBLE_PRECISION
REAL(MK) , DIMENSION(:) , INTENT(IN ) :: kpar
#elif __KIND == __SINGLE_PRECISION_COMPLEX | __KIND == __DOUBLE_PRECISION_COMPLEX
COMPLEX(MK), DIMENSION(:) , INTENT(IN ) :: kpar
#endif
!!! Kernel parameters. See documentation or ppm_comp_pp_kernels.inc for
!!! description. Type can be single, double, single complex or double
!!! complex. Omit this argument when using a lookup table.
#endif
REAL(MK) , INTENT( OUT) :: dxtableinv
!!! inverse of the dx (kernel evaluation locations) used to create
!!! the table. To use the table, look up the value
!!! `ktab(INT(dxtableinv*x))`.
INTEGER , INTENT( OUT) :: info
!!! Returns status, 0 upon success
!-------------------------------------------------------------------------
! Local variables
!-------------------------------------------------------------------------
INTEGER :: i,atable,iopt,idx
INTEGER, DIMENSION(1) :: ldl,ldu
REAL(MK), PARAMETER :: Rmin = 0.0_MK
REAL(MK) :: dxtable,Rmax,factor,factor2
REAL(MK) :: t0,dij,dij2,dij4,dij5,dx,dy,dz
#if __KIND == __SINGLE_PRECISION | __KIND == __DOUBLE_PRECISION
REAL(MK) :: eta
#elif __KIND == __SINGLE_PRECISION_COMPLEX | __KIND == __DOUBLE_PRECISION_COMPLEX
COMPLEX(MK) :: eta
#endif
!-------------------------------------------------------------------------
! Externals
!-------------------------------------------------------------------------
#if __KERNEL == __USER_FUNCTION
#if __KIND == __SINGLE_PRECISION | __KIND == __DOUBLE_PRECISION
REAL(MK) , DIMENSION(:) , INTENT(IN ) :: kpar
INTERFACE
FUNCTION kernel(x,kpar)
IMPLICIT NONE
#if __KIND == __SINGLE_PRECISION | __KIND == __SINGLE_PRECISION_COMPLEX
INTEGER, PARAMETER :: MK = KIND(1.0E0)
#elif __KIND == __DOUBLE_PRECISION | __KIND == __DOUBLE_PRECISION_COMPLEX
INTEGER, PARAMETER :: MK = KIND(1.0D0)
#endif
REAL(MK), INTENT(IN) :: x
REAL(MK), DIMENSION(:), INTENT(IN) :: kpar
REAL(MK) :: kernel
END FUNCTION kernel
END INTERFACE
#elif __KIND == __SINGLE_PRECISION_COMPLEX | __KIND == __DOUBLE_PRECISION_COMPLEX
COMPLEX(MK), DIMENSION(:) , INTENT(IN ) :: kpar
INTERFACE
FUNCTION kernel(x,kpar)
IMPLICIT NONE
#if __KIND == __SINGLE_PRECISION | __KIND == __SINGLE_PRECISION_COMPLEX
INTEGER, PARAMETER :: MK = KIND(1.0E0)
#elif __KIND == __DOUBLE_PRECISION | __KIND == __DOUBLE_PRECISION_COMPLEX
INTEGER, PARAMETER :: MK = KIND(1.0D0)
#endif
REAL(MK), INTENT(IN) :: x
COMPLEX(MK), DIMENSION(:), INTENT(IN) :: kpar
COMPLEX(MK) :: kernel
END FUNCTION kernel
END INTERFACE
#endif
#endif
!-------------------------------------------------------------------------
! Initialise
!-------------------------------------------------------------------------
CALL substart('ppm_comp_pp_mk_table',t0,info)
!-------------------------------------------------------------------------
! Check arguments.
!-------------------------------------------------------------------------
IF (ppm_debug .GT. 0) THEN
IF (ntable .LT. 1) THEN
info = ppm_error_error
CALL ppm_error(ppm_err_argument,'ppm_comp_pp_mk_table', &
& 'ntable must be >0',__LINE__,info)
GOTO 9999
ENDIF
IF (cutoff2 .LT. 0.0_MK) THEN
info = ppm_error_error
CALL ppm_error(ppm_err_argument,'ppm_comp_pp_mk_table', &
& 'cutoff2 must be >=0',__LINE__,info)
GOTO 9999
ENDIF
ENDIF
!-------------------------------------------------------------------------
! Compute Rmax
!-------------------------------------------------------------------------
Rmax = cutoff2
!-------------------------------------------------------------------------
! Allocate memory for table lookup
!-------------------------------------------------------------------------
atable = ntable-1
iopt = ppm_param_alloc_fit
ldl(1) = 0
ldu(1) = atable
CALL ppm_alloc(ktab,ldl,ldu,iopt,info)
IF (info .NE. ppm_param_success) THEN
info = ppm_error_fatal
CALL ppm_error(ppm_err_alloc,'ppm_comp_pp_mk_table', &
& 'kernel lookup table KTAB',__LINE__,info)
GOTO 9999
ENDIF
!-------------------------------------------------------------------------
! Compute values
!-------------------------------------------------------------------------
dxtable = (Rmax - Rmin)/REAL(atable,MK)
dxtableinv = 1.0_MK/dxtable
! This hack is needed since some kernels use dx,dy,dz explicitly. DO
! NOT MAKE TABLES FOR THOSE KERNELS!
dx = 0.0_MK
dy = 0.0_MK
dz = 0.0_MK
DO i=0,atable
dij = Rmin + REAL(i,MK)*dxtable
#include "ppm_comp_pp_kernels.inc"
ktab(i) = eta
ENDDO
!-------------------------------------------------------------------------
! Return
!-------------------------------------------------------------------------
9999 CONTINUE
CALL substop('ppm_comp_pp_mk_table',t0,info)
RETURN
#if __KERNEL == __INTERNAL
#if __KIND == __SINGLE_PRECISION
END SUBROUTINE ppm_comp_pp_mk_table_si
#elif __KIND == __DOUBLE_PRECISION
END SUBROUTINE ppm_comp_pp_mk_table_di
#elif __KIND == __SINGLE_PRECISION_COMPLEX
END SUBROUTINE ppm_comp_pp_mk_table_sci
#elif __KIND == __DOUBLE_PRECISION_COMPLEX
END SUBROUTINE ppm_comp_pp_mk_table_dci
#endif
#elif __KERNEL == __USER_FUNCTION
#if __KIND == __SINGLE_PRECISION
END SUBROUTINE ppm_comp_pp_mk_table_su
#elif __KIND == __DOUBLE_PRECISION
END SUBROUTINE ppm_comp_pp_mk_table_du
#elif __KIND == __SINGLE_PRECISION_COMPLEX
END SUBROUTINE ppm_comp_pp_mk_table_scu
#elif __KIND == __DOUBLE_PRECISION_COMPLEX
END SUBROUTINE ppm_comp_pp_mk_table_dcu
#endif
#endif
src/ppm_comp_pp_ring.f 0 → 100644
+ 503
− 0
View file @ dd7b0757
This diff is collapsed.
src/ppm_comp_pp_verlet.f 0 → 100644
+ 493
− 0
View file @ dd7b0757
This diff is collapsed.
src/ppm_module_comp_part.f 0 → 100644
+ 21
− 0
View file @ dd7b0757
!-------------------------------------------------------------------------
! Module : ppm_module_comp_part
!-------------------------------------------------------------------------
! Parallel Particle Mesh Library (PPM)
! ETH Zurich
! CH-8092 Zurich, Switzerland
!-------------------------------------------------------------------------
MODULE ppm_module_comp_part
!!! This module contains all user-callable routines to
!!! compute kernel-PP interactions on a set of particles.
!----------------------------------------------------------------------
! PPM modules
!----------------------------------------------------------------------
USE ppm_module_comp_pp_verlet
USE ppm_module_comp_pp_cell
USE ppm_module_comp_pp_ring
USE ppm_module_comp_pp_correct
USE ppm_module_comp_pp_mk_table
END MODULE ppm_module_comp_part
src/ppm_module_comp_pp_cell.f 0 → 100644
+ 101
− 0
View file @ dd7b0757
!-------------------------------------------------------------------------
! Module : ppm_module_comp_pp_cell
!-------------------------------------------------------------------------
! Parallel Particle Mesh Library (PPM)
! ETH Zurich
! CH-8092 Zurich, Switzerland
!-------------------------------------------------------------------------
!-------------------------------------------------------------------------
! Define types
!-------------------------------------------------------------------------
#define __SINGLE_PRECISION 1
#define __DOUBLE_PRECISION 2
#define __SINGLE_PRECISION_COMPLEX 3
#define __DOUBLE_PRECISION_COMPLEX 4
#define __INTERNAL 5
#define __USER_FUNCTION 6
#define __LOOKUP_TABLE 7
MODULE ppm_module_comp_pp_cell
!!! This module provides the routines concerned
!!! with computing kernel interactions within a set of particles.
!----------------------------------------------------------------------
! Define interfaces to the cell list versions
!----------------------------------------------------------------------
INTERFACE ppm_comp_pp_cell
MODULE PROCEDURE ppm_comp_pp_cell_si
MODULE PROCEDURE ppm_comp_pp_cell_di
MODULE PROCEDURE ppm_comp_pp_cell_sci
MODULE PROCEDURE ppm_comp_pp_cell_dci
MODULE PROCEDURE ppm_comp_pp_cell_su
MODULE PROCEDURE ppm_comp_pp_cell_du
MODULE PROCEDURE ppm_comp_pp_cell_scu
MODULE PROCEDURE ppm_comp_pp_cell_dcu
MODULE PROCEDURE ppm_comp_pp_cell_st
MODULE PROCEDURE ppm_comp_pp_cell_dt
MODULE PROCEDURE ppm_comp_pp_cell_sct
MODULE PROCEDURE ppm_comp_pp_cell_dct
END INTERFACE
!----------------------------------------------------------------------
! include the source
!----------------------------------------------------------------------
CONTAINS
#define __KERNEL __INTERNAL
#define __KIND __SINGLE_PRECISION
#include "ppm_comp_pp_cell.f"
#undef __KIND
#define __KIND __DOUBLE_PRECISION
#include "ppm_comp_pp_cell.f"
#undef __KIND
#define __KIND __SINGLE_PRECISION_COMPLEX
#include "ppm_comp_pp_cell.f"
#undef __KIND
#define __KIND __DOUBLE_PRECISION_COMPLEX
#include "ppm_comp_pp_cell.f"
#undef __KIND
#undef __KERNEL
#define __KERNEL __USER_FUNCTION
#define __KIND __SINGLE_PRECISION
#include "ppm_comp_pp_cell.f"
#undef __KIND
#define __KIND __DOUBLE_PRECISION
#include "ppm_comp_pp_cell.f"
#undef __KIND
#define __KIND __SINGLE_PRECISION_COMPLEX
#include "ppm_comp_pp_cell.f"
#undef __KIND
#define __KIND __DOUBLE_PRECISION_COMPLEX
#include "ppm_comp_pp_cell.f"
#undef __KIND
#undef __KERNEL
#define __KERNEL __LOOKUP_TABLE
#define __KIND __SINGLE_PRECISION
#include "ppm_comp_pp_cell.f"
#undef __KIND
#define __KIND __DOUBLE_PRECISION
#include "ppm_comp_pp_cell.f"
#undef __KIND
#define __KIND __SINGLE_PRECISION_COMPLEX
#include "ppm_comp_pp_cell.f"
#undef __KIND
#define __KIND __DOUBLE_PRECISION_COMPLEX
#include "ppm_comp_pp_cell.f"
#undef __KIND
#undef __KERNEL
END MODULE ppm_module_comp_pp_cell
src/ppm_module_comp_pp_correct.f 0 → 100644
+ 101
− 0
View file @ dd7b0757
!-------------------------------------------------------------------------
! Module : ppm_module_comp_pp_correct
!-------------------------------------------------------------------------
! Parallel Particle Mesh Library (PPM)
! ETH Zurich
! CH-8092 Zurich, Switzerland
!-------------------------------------------------------------------------
!-------------------------------------------------------------------------
! Define types
!-------------------------------------------------------------------------
#define __SINGLE_PRECISION 1
#define __DOUBLE_PRECISION 2
#define __SINGLE_PRECISION_COMPLEX 3
#define __DOUBLE_PRECISION_COMPLEX 4
#define __INTERNAL 5
#define __USER_FUNCTION 6
#define __LOOKUP_TABLE 7
MODULE ppm_module_comp_pp_correct
!!! This module provides the routines concerned
!!! with computing kernel interactions within a set of particles.
!----------------------------------------------------------------------
! Define interfaces to the discretisation correction routine
!----------------------------------------------------------------------
INTERFACE ppm_comp_pp_correct
MODULE PROCEDURE ppm_comp_pp_correct_si
MODULE PROCEDURE ppm_comp_pp_correct_di
MODULE PROCEDURE ppm_comp_pp_correct_sci
MODULE PROCEDURE ppm_comp_pp_correct_dci
MODULE PROCEDURE ppm_comp_pp_correct_su
MODULE PROCEDURE ppm_comp_pp_correct_du
MODULE PROCEDURE ppm_comp_pp_correct_scu
MODULE PROCEDURE ppm_comp_pp_correct_dcu
MODULE PROCEDURE ppm_comp_pp_correct_st
MODULE PROCEDURE ppm_comp_pp_correct_dt
MODULE PROCEDURE ppm_comp_pp_correct_sct
MODULE PROCEDURE ppm_comp_pp_correct_dct
END INTERFACE
!----------------------------------------------------------------------
! include the source
!----------------------------------------------------------------------
CONTAINS
#define __KERNEL __INTERNAL
#define __KIND __SINGLE_PRECISION
#include "ppm_comp_pp_correct.f"
#undef __KIND
#define __KIND __DOUBLE_PRECISION
#include "ppm_comp_pp_correct.f"
#undef __KIND
#define __KIND __SINGLE_PRECISION_COMPLEX
#include "ppm_comp_pp_correct.f"
#undef __KIND
#define __KIND __DOUBLE_PRECISION_COMPLEX
#include "ppm_comp_pp_correct.f"
#undef __KIND
#undef __KERNEL
#define __KERNEL __USER_FUNCTION
#define __KIND __SINGLE_PRECISION
#include "ppm_comp_pp_correct.f"
#undef __KIND
#define __KIND __DOUBLE_PRECISION
#include "ppm_comp_pp_correct.f"
#undef __KIND
#define __KIND __SINGLE_PRECISION_COMPLEX
#include "ppm_comp_pp_correct.f"
#undef __KIND
#define __KIND __DOUBLE_PRECISION_COMPLEX
#include "ppm_comp_pp_correct.f"
#undef __KIND
#undef __KERNEL
#define __KERNEL __LOOKUP_TABLE
#define __KIND __SINGLE_PRECISION
#include "ppm_comp_pp_correct.f"
#undef __KIND
#define __KIND __DOUBLE_PRECISION
#include "ppm_comp_pp_correct.f"
#undef __KIND
#define __KIND __SINGLE_PRECISION_COMPLEX
#include "ppm_comp_pp_correct.f"
#undef __KIND
#define __KIND __DOUBLE_PRECISION_COMPLEX
#include "ppm_comp_pp_correct.f"
#undef __KIND
#undef __KERNEL
END MODULE ppm_module_comp_pp_correct
src/ppm_module_comp_pp_doring.f 0 → 100644
+ 101
− 0
View file @ dd7b0757
!-------------------------------------------------------------------------
! Module : ppm_module_comp_pp_doring
!-------------------------------------------------------------------------
! Parallel Particle Mesh Library (PPM)
! ETH Zurich
! CH-8092 Zurich, Switzerland
!-------------------------------------------------------------------------
!-------------------------------------------------------------------------
! Define types
!-------------------------------------------------------------------------
#define __SINGLE_PRECISION 1
#define __DOUBLE_PRECISION 2
#define __SINGLE_PRECISION_COMPLEX 3
#define __DOUBLE_PRECISION_COMPLEX 4
#define __INTERNAL 5
#define __USER_FUNCTION 6
#define __LOOKUP_TABLE 7
MODULE ppm_module_comp_pp_doring
!!! This module provides the routines concerned
!!! with computing kernel interactions within a set of particles.
!----------------------------------------------------------------------
! Define interfaces to the ring interaction subroutine
!----------------------------------------------------------------------
INTERFACE ppm_comp_pp_doring
MODULE PROCEDURE ppm_comp_pp_doring_si
MODULE PROCEDURE ppm_comp_pp_doring_di
MODULE PROCEDURE ppm_comp_pp_doring_sci
MODULE PROCEDURE ppm_comp_pp_doring_dci
MODULE PROCEDURE ppm_comp_pp_doring_su
MODULE PROCEDURE ppm_comp_pp_doring_du
MODULE PROCEDURE ppm_comp_pp_doring_scu
MODULE PROCEDURE ppm_comp_pp_doring_dcu
MODULE PROCEDURE ppm_comp_pp_doring_st
MODULE PROCEDURE ppm_comp_pp_doring_dt
MODULE PROCEDURE ppm_comp_pp_doring_sct
MODULE PROCEDURE ppm_comp_pp_doring_dct
END INTERFACE
!----------------------------------------------------------------------
! include the source
!----------------------------------------------------------------------
CONTAINS
#define __KERNEL __INTERNAL
#define __KIND __SINGLE_PRECISION
#include "ppm_comp_pp_doring.f"
#undef __KIND
#define __KIND __DOUBLE_PRECISION
#include "ppm_comp_pp_doring.f"
#undef __KIND
#define __KIND __SINGLE_PRECISION_COMPLEX
#include "ppm_comp_pp_doring.f"
#undef __KIND
#define __KIND __DOUBLE_PRECISION_COMPLEX
#include "ppm_comp_pp_doring.f"
#undef __KIND
#undef __KERNEL
#define __KERNEL __USER_FUNCTION
#define __KIND __SINGLE_PRECISION
#include "ppm_comp_pp_doring.f"
#undef __KIND
#define __KIND __DOUBLE_PRECISION
#include "ppm_comp_pp_doring.f"
#undef __KIND
#define __KIND __SINGLE_PRECISION_COMPLEX
#include "ppm_comp_pp_doring.f"
#undef __KIND
#define __KIND __DOUBLE_PRECISION_COMPLEX
#include "ppm_comp_pp_doring.f"
#undef __KIND
#undef __KERNEL
#define __KERNEL __LOOKUP_TABLE
#define __KIND __SINGLE_PRECISION
#include "ppm_comp_pp_doring.f"
#undef __KIND
#define __KIND __DOUBLE_PRECISION
#include "ppm_comp_pp_doring.f"
#undef __KIND
#define __KIND __SINGLE_PRECISION_COMPLEX
#include "ppm_comp_pp_doring.f"
#undef __KIND
#define __KIND __DOUBLE_PRECISION_COMPLEX
#include "ppm_comp_pp_doring.f"
#undef __KIND
#undef __KERNEL
END MODULE ppm_module_comp_pp_doring
src/ppm_module_comp_pp_mk_table.f 0 → 100644
+ 79
− 0
View file @ dd7b0757
!-------------------------------------------------------------------------
! Module : ppm_module_comp_pp_mk_table
!-------------------------------------------------------------------------
! Parallel Particle Mesh Library (PPM)
! ETH Zurich
! CH-8092 Zurich, Switzerland
!-------------------------------------------------------------------------
!-------------------------------------------------------------------------
! Define types
!-------------------------------------------------------------------------
#define __SINGLE_PRECISION 1
#define __DOUBLE_PRECISION 2
#define __SINGLE_PRECISION_COMPLEX 3
#define __DOUBLE_PRECISION_COMPLEX 4
#define __INTERNAL 5
#define __USER_FUNCTION 6
MODULE ppm_module_comp_pp_mk_table
!!! [[ppm_comp_pp_mk_table]]
!!! This module provides the routines concerned
!!! with computing kernel interactions within a set of particles.
!----------------------------------------------------------------------
! Define interfaces to the lookup table creation routine
!----------------------------------------------------------------------
INTERFACE ppm_comp_pp_mk_table
MODULE PROCEDURE ppm_comp_pp_mk_table_si
MODULE PROCEDURE ppm_comp_pp_mk_table_di
MODULE PROCEDURE ppm_comp_pp_mk_table_sci
MODULE PROCEDURE ppm_comp_pp_mk_table_dci
MODULE PROCEDURE ppm_comp_pp_mk_table_su
MODULE PROCEDURE ppm_comp_pp_mk_table_du
MODULE PROCEDURE ppm_comp_pp_mk_table_scu
MODULE PROCEDURE ppm_comp_pp_mk_table_dcu
END INTERFACE
!----------------------------------------------------------------------
! include the source
!----------------------------------------------------------------------
CONTAINS
#define __KERNEL __INTERNAL
#define __KIND __SINGLE_PRECISION
#include "ppm_comp_pp_mk_table.f"
#undef __KIND
#define __KIND __DOUBLE_PRECISION
#include "ppm_comp_pp_mk_table.f"
#undef __KIND
#define __KIND __SINGLE_PRECISION_COMPLEX
#include "ppm_comp_pp_mk_table.f"
#undef __KIND
#define __KIND __DOUBLE_PRECISION_COMPLEX
#include "ppm_comp_pp_mk_table.f"
#undef __KIND
#undef __KERNEL
#define __KERNEL __USER_FUNCTION
#define __KIND __SINGLE_PRECISION
#include "ppm_comp_pp_mk_table.f"
#undef __KIND
#define __KIND __DOUBLE_PRECISION
#include "ppm_comp_pp_mk_table.f"
#undef __KIND
#define __KIND __SINGLE_PRECISION_COMPLEX
#include "ppm_comp_pp_mk_table.f"
#undef __KIND
#define __KIND __DOUBLE_PRECISION_COMPLEX
#include "ppm_comp_pp_mk_table.f"
#undef __KIND
#undef __KERNEL
END MODULE ppm_module_comp_pp_mk_table
src/ppm_module_comp_pp_ring.f 0 → 100644
+ 101
− 0
View file @ dd7b0757
!-------------------------------------------------------------------------
! Module : ppm_module_comp_pp_ring
!-------------------------------------------------------------------------
! Parallel Particle Mesh Library (PPM)
! ETH Zurich
! CH-8092 Zurich, Switzerland
!-------------------------------------------------------------------------
!-------------------------------------------------------------------------
! Define types
!-------------------------------------------------------------------------
#define __SINGLE_PRECISION 1
#define __DOUBLE_PRECISION 2
#define __SINGLE_PRECISION_COMPLEX 3
#define __DOUBLE_PRECISION_COMPLEX 4
#define __INTERNAL 5
#define __USER_FUNCTION 6
#define __LOOKUP_TABLE 7
MODULE ppm_module_comp_pp_ring
!!! This module provides the routines concerned
!!! with computing kernel interactions within a set of particles.
!----------------------------------------------------------------------
! Define interfaces to the ring versions
!----------------------------------------------------------------------
INTERFACE ppm_comp_pp_ring
MODULE PROCEDURE ppm_comp_pp_ring_si
MODULE PROCEDURE ppm_comp_pp_ring_di
MODULE PROCEDURE ppm_comp_pp_ring_sci
MODULE PROCEDURE ppm_comp_pp_ring_dci
MODULE PROCEDURE ppm_comp_pp_ring_su
MODULE PROCEDURE ppm_comp_pp_ring_du
MODULE PROCEDURE ppm_comp_pp_ring_scu
MODULE PROCEDURE ppm_comp_pp_ring_dcu
MODULE PROCEDURE ppm_comp_pp_ring_st
MODULE PROCEDURE ppm_comp_pp_ring_dt
MODULE PROCEDURE ppm_comp_pp_ring_sct
MODULE PROCEDURE ppm_comp_pp_ring_dct
END INTERFACE
!----------------------------------------------------------------------
! include the source
!----------------------------------------------------------------------
CONTAINS
#define __KERNEL __INTERNAL
#define __KIND __SINGLE_PRECISION
#include "ppm_comp_pp_ring.f"
#undef __KIND
#define __KIND __DOUBLE_PRECISION
#include "ppm_comp_pp_ring.f"
#undef __KIND
#define __KIND __SINGLE_PRECISION_COMPLEX
#include "ppm_comp_pp_ring.f"
#undef __KIND
#define __KIND __DOUBLE_PRECISION_COMPLEX
#include "ppm_comp_pp_ring.f"
#undef __KIND
#undef __KERNEL
#define __KERNEL __USER_FUNCTION
#define __KIND __SINGLE_PRECISION
#include "ppm_comp_pp_ring.f"
#undef __KIND
#define __KIND __DOUBLE_PRECISION
#include "ppm_comp_pp_ring.f"
#undef __KIND
#define __KIND __SINGLE_PRECISION_COMPLEX
#include "ppm_comp_pp_ring.f"
#undef __KIND
#define __KIND __DOUBLE_PRECISION_COMPLEX
#include "ppm_comp_pp_ring.f"
#undef __KIND
#undef __KERNEL
#define __KERNEL __LOOKUP_TABLE
#define __KIND __SINGLE_PRECISION
#include "ppm_comp_pp_ring.f"
#undef __KIND
#define __KIND __DOUBLE_PRECISION
#include "ppm_comp_pp_ring.f"
#undef __KIND
#define __KIND __SINGLE_PRECISION_COMPLEX
#include "ppm_comp_pp_ring.f"
#undef __KIND
#define __KIND __DOUBLE_PRECISION_COMPLEX
#include "ppm_comp_pp_ring.f"
#undef __KIND
#undef __KERNEL
END MODULE ppm_module_comp_pp_ring
src/ppm_module_comp_pp_verlet.f 0 → 100644
+ 101
− 0
View file @ dd7b0757
!-------------------------------------------------------------------------
! Module : ppm_module_comp_pp_verlet
!-------------------------------------------------------------------------
! Parallel Particle Mesh Library (PPM)
! ETH Zurich
! CH-8092 Zurich, Switzerland
!-------------------------------------------------------------------------
!-------------------------------------------------------------------------
! Define types
!-------------------------------------------------------------------------
#define __SINGLE_PRECISION 1
#define __DOUBLE_PRECISION 2
#define __SINGLE_PRECISION_COMPLEX 3
#define __DOUBLE_PRECISION_COMPLEX 4
#define __INTERNAL 5
#define __USER_FUNCTION 6
#define __LOOKUP_TABLE 7
MODULE ppm_module_comp_pp_verlet
!!! This module provides the routines concerned
!!! with computing kernel interactions within a set of particles.
!----------------------------------------------------------------------
! Define interfaces to the Verlet list versions
!----------------------------------------------------------------------
INTERFACE ppm_comp_pp_verlet
MODULE PROCEDURE ppm_comp_pp_verlet_si
MODULE PROCEDURE ppm_comp_pp_verlet_di
MODULE PROCEDURE ppm_comp_pp_verlet_sci
MODULE PROCEDURE ppm_comp_pp_verlet_dci
MODULE PROCEDURE ppm_comp_pp_verlet_su
MODULE PROCEDURE ppm_comp_pp_verlet_du
MODULE PROCEDURE ppm_comp_pp_verlet_scu
MODULE PROCEDURE ppm_comp_pp_verlet_dcu
MODULE PROCEDURE ppm_comp_pp_verlet_st
MODULE PROCEDURE ppm_comp_pp_verlet_dt
MODULE PROCEDURE ppm_comp_pp_verlet_sct
MODULE PROCEDURE ppm_comp_pp_verlet_dct
END INTERFACE
!----------------------------------------------------------------------
! include the source
!----------------------------------------------------------------------
CONTAINS
#define __KERNEL __INTERNAL
#define __KIND __SINGLE_PRECISION
#include "ppm_comp_pp_verlet.f"
#undef __KIND
#define __KIND __DOUBLE_PRECISION
#include "ppm_comp_pp_verlet.f"
#undef __KIND
#define __KIND __SINGLE_PRECISION_COMPLEX
#include "ppm_comp_pp_verlet.f"
#undef __KIND
#define __KIND __DOUBLE_PRECISION_COMPLEX
#include "ppm_comp_pp_verlet.f"
#undef __KIND
#undef __KERNEL
#define __KERNEL __USER_FUNCTION
#define __KIND __SINGLE_PRECISION
#include "ppm_comp_pp_verlet.f"
#undef __KIND
#define __KIND __DOUBLE_PRECISION
#include "ppm_comp_pp_verlet.f"
#undef __KIND
#define __KIND __SINGLE_PRECISION_COMPLEX
#include "ppm_comp_pp_verlet.f"
#undef __KIND
#define __KIND __DOUBLE_PRECISION_COMPLEX
#include "ppm_comp_pp_verlet.f"
#undef __KIND
#undef __KERNEL
#define __KERNEL __LOOKUP_TABLE
#define __KIND __SINGLE_PRECISION
#include "ppm_comp_pp_verlet.f"
#undef __KIND
#define __KIND __DOUBLE_PRECISION
#include "ppm_comp_pp_verlet.f"
#undef __KIND
#define __KIND __SINGLE_PRECISION_COMPLEX
#include "ppm_comp_pp_verlet.f"
#undef __KIND
#define __KIND __DOUBLE_PRECISION_COMPLEX
#include "ppm_comp_pp_verlet.f"
#undef __KIND
#undef __KERNEL
END MODULE ppm_module_comp_pp_verlet
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