From 5071458be0e27c7228dc92f639d7fe56814a608a Mon Sep 17 00:00:00 2001 From: oawile <oawile@7c7fe9aa-52eb-4d9e-b0a8-ba7d787348e9> Date: Tue, 9 Mar 2010 14:01:26 +0000 Subject: [PATCH] - added back the ppm_template_comp_pp_* files git-svn-id: https://ppm.inf.ethz.ch/svn/ppmnumerics/branches/ngtopo/libppmnumerics@556 7c7fe9aa-52eb-4d9e-b0a8-ba7d787348e9 --- src/ppm_template_comp_pp_cell.f | 469 ++++++++++++++++++++++++++++++ src/ppm_template_comp_pp_ring.f | 218 ++++++++++++++ src/ppm_template_comp_pp_verlet.f | 303 +++++++++++++++++++ 3 files changed, 990 insertions(+) create mode 100644 src/ppm_template_comp_pp_cell.f create mode 100644 src/ppm_template_comp_pp_ring.f create mode 100644 src/ppm_template_comp_pp_verlet.f diff --git a/src/ppm_template_comp_pp_cell.f b/src/ppm_template_comp_pp_cell.f new file mode 100644 index 0000000..f10f11d --- /dev/null +++ b/src/ppm_template_comp_pp_cell.f @@ -0,0 +1,469 @@ + !------------------------------------------------------------------------- + ! Subroutine : ppm_template_comp_pp_cell + !------------------------------------------------------------------------- + ! + ! Purpose : Template for the user subroutine which computes + ! direct particle-particle interactions using cell + ! lists. How to use this template: + ! (1) copy this file to your application code dir + ! (2) rename the file and this routine (make sure + ! to also replace all the occurrences of the + ! subroutine name in all WRITE statements). + ! (3) change address in header and add your CVS + ! Log (if needed). + ! (4) declare and add additional arguments + ! (5) add USE and INCLUDE statements for your + ! modules and header files + ! (6) define the precision of floating point data + ! (maybe provide different versions of this + ! subroutine for different precisions) + ! (7) add particle-particle interactions where + ! needed (search for USER CODE). Make sure to + ! replace all occurrences of dimens with the + ! variable which actually contains the problem + ! dimensionality in your code. + ! (8) Optional: Make argument checking conditional + ! on the debug level of your application + ! (9) Optional: Make this routine use your own + ! error logging, message writing, allocation, + ! subroutine start, subroutine stop and + ! timing routines. + ! (10)Optional: Move the variable declarations for + ! clist and Nm to your client program and + ! add them as arguments to this routine. This + ! makes it possible to use several sets of + ! lists. + ! (11)Optional: Remove the imode=-1 and imode=1 + ! cases from this routine and build/destroy + ! the lists directly in your client program by + ! calling ppm_neighlist_clist/..destroy. You + ! can then remove imode from the argument list + ! of this routine. + ! + ! Input : xp(:,:) (F) particle co-ordinates + ! pdata(:,:) (F) particle data used for interaction + ! (e.g. vorticity, strength, ...) + ! ... (.) USER: add more input arguments as + ! needed + ! Np (I) number of particles on this proc. + ! cutoff (F) cutoff radius for PP interactions + ! lsymm (L) Do symmetric PP interactions? + ! .TRUE. : Yes + ! .FALSE. : No + ! imode (I) Mode of action. Any of the folowing: + ! 0 PP interactions + ! 1 build cell lists and return + ! -1 destroy cell lists and return + ! nsublist (I) number of subdomains on the local + ! processor + ! + ! Input/output : params(:) (F) user defined parameters and/or + ! output from the interaction (i.e. + ! potential energy) + ! + ! Output : dpd(:,:) (F) Change of particle data (pdata) due to + ! interaction. + ! info (I) return status. =0 if no error. + ! + ! Routines : ppm_neighlist_clist + ! ppm_clist_destroy + ! ppm_neighlist_MkNeighIdx + ! + ! Remarks : If particles have moved between calls or lsymm has + ! changed, always call with imode=1. + ! + ! seach for USER in this file and fill in + ! the particle-particle interactions in all 4 places. + ! Use the dpd array to return the result. It is your + ! responsibility to properly allocate and initialize + ! this array BEFORE calling this routine. + ! + ! After finishing all time steps, always call this + ! routine with imode=-1 to free the memory of the + ! cell lists. Failure to do so will result in a + ! memory leak in your program. + ! + ! If the CYCLE command is a problem on your hardware + ! (e.g. prevents vectorization), split the DO-loop in + ! two parts as: + ! DO jpart = istart,ipart-1 + ! ... + ! ENDDO + ! DO jpart = ipart+1,iend + ! ... + ! ENDDO + ! + ! References : + ! + ! Revisions : + !------------------------------------------------------------------------- + ! $Revision: 1.1.1.1 $ + ! Revision 1.4 2004/02/24 11:35:53 ivos + ! Revision 1.3 2004/02/19 14:01:39 gonnetp + ! Revision 1.2 2004/02/04 17:20:48 ivos + ! Revision 1.1 2004/01/26 17:24:35 ivos + !------------------------------------------------------------------------- + ! Parallel Particle Mesh Library (PPM) + ! Institute of Computational Science + ! ETH Zentrum, Hirschengraben 84 + ! CH-8092 Zurich, Switzerland + !------------------------------------------------------------------------- + + SUBROUTINE ppm_template_comp_pp_cell(xp,pdata,Np,cutoff,lsymm, & + & imode,nsublist,dpd,params,info) + + !------------------------------------------------------------------------- + ! Modules + !------------------------------------------------------------------------- + ! USER: USE modules here + USE ppm_module_neighlist + IMPLICIT NONE + !------------------------------------------------------------------------- + ! Includes + !------------------------------------------------------------------------- + INCLUDE 'ppm_param.h' + !------------------------------------------------------------------------- + ! USER: Define precision + !------------------------------------------------------------------------- +! INTEGER, PARAMETER :: MK = ppm_kind_single + INTEGER, PARAMETER :: MK = ppm_kind_double + !------------------------------------------------------------------------- + ! Arguments + !------------------------------------------------------------------------- + REAL(MK), DIMENSION(:,:), INTENT(IN ) :: xp + REAL(MK), DIMENSION(:,:), INTENT(IN ) :: pdata + ! USER: add more input arguments as needed + INTEGER , INTENT(IN ) :: Np + REAL(MK) , INTENT(IN ) :: cutoff + LOGICAL , INTENT(IN ) :: lsymm + INTEGER , INTENT(IN ) :: imode,nsublist + REAL(MK), DIMENSION(:,:), INTENT(INOUT) :: dpd + REAL(MK), DIMENSION(:), INTENT(INOUT) :: params + INTEGER , INTENT( OUT) :: info + !------------------------------------------------------------------------- + ! Local variables + !------------------------------------------------------------------------- + ! counters + INTEGER :: i,idom,ibox,jbox + INTEGER :: ipart,jpart,ip,jp,isize + INTEGER :: cbox,iinter,j,k + ! coordinate differences + REAL(MK) :: dx,dy,dz + ! square of inter particle distance + REAL(MK) :: dij + ! cutoff squared + REAL(MK), SAVE :: cut2 + ! start and end particle in a box + INTEGER :: istart,iend,jstart,jend + ! box size for cell list + REAL(MK), DIMENSION(3) :: bsize + ! cell list + ! USER: if allocation of the following varaible fails due to stack + ! size limitations, try putting it in a module and add a USE + ! statement for it above. In order to use several Cell lists, move + ! these declarations to the client program and add clist and Nm + ! to the argument list of this routine. + TYPE(ppm_type_ptr_to_clist), DIMENSION(:), POINTER, SAVE :: clist + ! number of cells in all directions + INTEGER, DIMENSION(:,:), POINTER, SAVE :: Nm + ! cell neighbor lists + INTEGER, DIMENSION(:,:), POINTER, SAVE :: inp,jnp + ! number of interactions for each cell + INTEGER, SAVE :: nnp + ! cell offsets for box index + INTEGER :: n1,n2,nz + !------------------------------------------------------------------------- + ! Externals + !------------------------------------------------------------------------- + + !------------------------------------------------------------------------- + ! Initialise + !------------------------------------------------------------------------- + info = 0 + + !------------------------------------------------------------------------- + ! Check arguments. + ! USER: you may want to make this conditional on the debug level of + ! your application... + !------------------------------------------------------------------------- + IF (cutoff .LE. 0.0_MK) THEN + WRITE(*,'(2A)') '(ppm_template_comp_pp_cell): ', & + & 'cutoff must be >0' + info = -1 + GOTO 9999 + ENDIF + IF (Np .LE. 0) THEN + WRITE(*,'(2A)') '(ppm_template_comp_pp_cell): ', & + & 'Np must be >0' + info = -1 + GOTO 9999 + ENDIF + + !------------------------------------------------------------------------- + ! If imode = -1, destroy cell lists and return + ! USER: This could also be done by your client directly. The code is + ! just here as a template. + !------------------------------------------------------------------------- + IF (imode .EQ. -1) THEN + CALL ppm_clist_destroy(clist,info) + DEALLOCATE(inp,jnp,Nm,STAT=info) + GOTO 9999 + ENDIF + + !------------------------------------------------------------------------- + ! Build new cell lists if needed + ! USER: This could also be done by your client directly. The code is + ! just here as a template. + !------------------------------------------------------------------------- + IF (imode .EQ. 1) THEN + !--------------------------------------------------------------------- + ! Destroy old cell list (if there already was one) + !--------------------------------------------------------------------- + CALL ppm_clist_destroy(clist,info) + + !--------------------------------------------------------------------- + ! The size of the cells in each spatial direction should at last + ! be the size of the cutoff. + !--------------------------------------------------------------------- + bsize(1:3) = cutoff + cut2 = cutoff*cutoff + + !--------------------------------------------------------------------- + ! Generate cell lists + !--------------------------------------------------------------------- + CALL ppm_neighlist_clist(xp,Np,bsize,lsymm,clist,Nm,info) + IF (info .NE. 0) THEN + WRITE(*,'(2A)') '(ppm_template_comp_pp_cell): ', & + & 'Building cell lists failed ' + info = -1 + GOTO 9999 + ENDIF + + !--------------------------------------------------------------------- + ! Generate cell neighbor lists + !--------------------------------------------------------------------- + CALL ppm_neighlist_MkNeighIdx(lsymm,inp,jnp,nnp,info) + GOTO 9999 + ENDIF + + !------------------------------------------------------------------------- + ! PARTICLE-PARTICLE INTERACTIONS using symmetry + ! This will not vectorize and has therefore been moved into a + ! separate loop !! + !------------------------------------------------------------------------- + IF (lsymm) THEN + DO idom=1,nsublist + n1 = Nm(1,idom) + n2 = Nm(1,idom)*Nm(2,idom) + nz = Nm(3,idom) + IF (dimens .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 cell 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-1 + ip = clist(idom)%lpdx(ipart) + 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 (dimens .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. cut2) CYCLE + !------------------------------------- + ! Particle ip interacts with + ! particle jp here... and + ! vice versa to use symmetry. + !------------------------------------- + ! + ! USER CODE HERE: + ! dpd(:,ip) = f(pdata(:,ip), + ! pdata(:,jp),dij) + ! dpd(:,jp) = -dpd(:,ip) + ENDDO + 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 empty + IF (jend .LT. jstart) CYCLE + ! loop over all particles inside this cell + DO ipart=istart,iend + ip = clist(idom)%lpdx(ipart) + ! 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 (dimens .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. cut2) CYCLE + !------------------------------------- + ! Particle ip interacts with + ! particle jp here... and + ! vice versa to use symmetry. + !------------------------------------- + ! + ! USER CODE HERE: + ! dpd(:,ip) = f(pdata(:,ip), + ! pdata(:,jp),dij) + ! dpd(:,jp) = -dpd(:,ip) + ENDDO + 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 (dimens .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 cell 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) + 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 (dimens .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. cut2) CYCLE + !--------------------------------- + ! Particle ip interacts with + ! particle jp. + !--------------------------------- + ! + ! USER CODE HERE: + ! dpd(:,ip) = f(pdata(:,ip), + ! pdata(:,jp),dij) + ENDDO + 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 this iinter if empty + IF (jend .LT. jstart) CYCLE + ! loop over all particles inside this cell + DO ipart=istart,iend + ip = clist(idom)%lpdx(ipart) + ! 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 (dimens .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. cut2) CYCLE + !--------------------------------- + ! Particle ip interacts with + ! particle jp. + !--------------------------------- + ! + ! USER CODE HERE: + ! dpd(:,ip) = f(pdata(:,ip), + ! pdata(:,jp),dij) + ENDDO + ENDDO + ENDIF ! ibox .EQ. jbox + ENDDO ! iinter + ENDDO ! i + ENDDO ! j + ENDDO ! k + ENDDO ! idom + ENDIF + + !------------------------------------------------------------------------- + ! Return + !------------------------------------------------------------------------- + 9999 CONTINUE + RETURN + END SUBROUTINE ppm_template_comp_pp_cell diff --git a/src/ppm_template_comp_pp_ring.f b/src/ppm_template_comp_pp_ring.f new file mode 100644 index 0000000..e131234 --- /dev/null +++ b/src/ppm_template_comp_pp_ring.f @@ -0,0 +1,218 @@ + !------------------------------------------------------------------------- + ! Subroutine : ppm_template_comp_pp_ring + !------------------------------------------------------------------------- + ! + ! Purpose : Template for the user subroutine which computes + ! direct particle-particle interactions using the + ! ring topology. How to use this template: + ! (1) copy this file to your application code dir + ! (2) rename the file and this routine (make sure + ! to also replace all the occurrences of the + ! subroutine name in all WRITE statements). + ! (3) change address in header and add your CVS + ! Log (if needed). + ! (4) declare and add additional arguments + ! (5) add USE and INCLUDE statements for your + ! modules and header files + ! (6) define the precision of floating point data + ! (maybe provide different versions of this + ! subroutine for different precisions) + ! (7) add particle-particle interactions where + ! needed (search for USER CODE) + ! (8) Optional: Make argument checking conditional + ! on the debug level of your application + ! (9) Optional: Make this routine use your own + ! error logging, message writing, allocation, + ! subroutine start, subroutine stop and + ! timing routines. + ! + ! Input : xp(:,:) (F) : particle co-ordinates [Group 1] + ! vp(:,:) (F) : particle data used for interaction + ! (e.g. vorticity, strength, ...) + ! [Group 1] + ! ... (.) : USER: add more input arguments for + ! group 1 here if needed. + ! needed + ! Npart (I) : number of particles on this proc. + ! [Group 1] + ! xp2(:,:) (F) : particle co-ordinates [Group 2] + ! vp2(:,:) (F) : particle data used for interaction + ! (e.g. vorticity, strength, ...) + ! [Group 2] + ! ... (.) : USER: add more input arguments for + ! group 1 here if needed. + ! Lpart (I) : number of particles [Group 2] + ! lsymm (L) : Whether to use symmetry or not: + ! .FALSE. pp interaction w/o symmetry + ! .TRUE. pp interaction w/ symmetry + ! mode (I) : Whether the two groups are the + ! same or not: + ! 0 not the same group + ! 1 the same group + ! + ! Input/output : params(:) (F) : user defined parameters and/or + ! output from the interaction (i.e. + ! potential energy) + ! fp(:,:) (F) : Change of particle data due + ! to interaction [Group 1] + ! fp2(:,:) (F) : Change of particle data due + ! to interaction [Group 2] + ! ... (.) : USER: add more output arguments for + ! if needed. + ! Output : info (I) : return status. 0 if no error. + ! + ! Routines : + ! + ! Remarks : Search for USER in this file and fill in + ! the particle-particle interactions in all 4 places. + ! Use the fp and fp2 array to return the result. + ! + ! + ! References : + ! + ! Revisions : + !------------------------------------------------------------------------- + ! $Revision: 1.1.1.1 $ + ! Revision 1.2 2004/04/23 17:25:33 oingo + ! Revision 1.1 2004/04/22 08:29:06 oingo + !------------------------------------------------------------------------- + ! Parallel Particle Mesh Library (PPM) + ! Institute of Computational Science + ! ETH Zentrum, Hirschengraben 84 + ! CH-8092 Zurich, Switzerland + !------------------------------------------------------------------------- + + SUBROUTINE ppm_template_comp_pp_ring(xp, vp, fp ,Npart, & + & xp2,vp2,fp2,Lpart, & + & lsymm,params,mode,info) + + !------------------------------------------------------------------------- + ! Modules + !------------------------------------------------------------------------- + ! USER: USE statements here... + IMPLICIT NONE + !------------------------------------------------------------------------- + ! Includes + !------------------------------------------------------------------------- + INCLUDE 'ppm_param.h' + !------------------------------------------------------------------------- + ! USER: Define precision + !------------------------------------------------------------------------- +! INTEGER, PARAMETER :: MK = ppm_kind_single + INTEGER, PARAMETER :: MK = ppm_kind_double + !------------------------------------------------------------------------- + ! Arguments + !------------------------------------------------------------------------- + REAL(MK), DIMENSION(:,:), INTENT(IN ) :: xp + REAL(MK), DIMENSION(:,:), INTENT(IN ) :: vp + ! USER: add more arguments here if needed + REAL(MK), DIMENSION(:,:), INTENT(INOUT) :: fp + INTEGER , INTENT(IN ) :: Npart + REAL(MK), DIMENSION(:,:), INTENT(IN ) :: xp2 + REAL(MK), DIMENSION(:,:), INTENT(IN ) :: vp2 + ! USER: add more arguments here if needed + REAL(MK), DIMENSION(:,:), INTENT(INOUT) :: fp2 + INTEGER , INTENT(IN ) :: Lpart + LOGICAL , INTENT(IN ) :: lsymm + REAL(MK), DIMENSION(:), INTENT(INOUT) :: params + INTEGER , INTENT(IN ) :: mode + INTEGER , INTENT( OUT) :: info + !------------------------------------------------------------------------- + ! Local variables + !------------------------------------------------------------------------- + ! counters + INTEGER :: i,j + ! coordinate differences + REAL(MK) :: dx,dy,dz + ! square of inter particle distance + REAL(MK) :: dij + + !------------------------------------------------------------------------- + ! Externals + !------------------------------------------------------------------------- + + !------------------------------------------------------------------------- + ! Initialise + !------------------------------------------------------------------------- + info = 0 + + !------------------------------------------------------------------------- + ! Check arguments + ! USER: You might want to make this conditional on the debug level + ! of your application... + !------------------------------------------------------------------------- + IF (Npart .LE. 0) THEN + WRITE(*,'(2A)') '(ppm_template_comp_pp_ring): ', & + & 'Npart must be >0' + info = -1 + GOTO 9999 + ENDIF + IF (Lpart .LE. 0) THEN + WRITE(*,'(2A)') '(ppm_template_comp_pp_ring): ', & + & 'Lpart must be >0' + info = -1 + GOTO 9999 + ENDIF + IF ((mode .NE. 0) .AND. (mode .NE. 1)) THEN + WRITE(*,'(2A)') '(ppm_template_comp_pp_ring): ', & + & 'MODE must be either 0 or 1' + info = -1 + GOTO 9999 + ENDIF + + !------------------------------------------------------------------------- + ! Check if we are computing the selfinteraction (i.e. Group1 .EQ. Group2) + !------------------------------------------------------------------------- + IF (mode .EQ. 1) THEN + IF (lsymm) THEN + DO i = 1,Npart + DO j = i+1,Lpart + !----------------------------------------------------------- + ! USER CODE HERE (with symmetry) + ! dx = xp(1,i) - xp2(1,j) + ! dy = xp(2,i) - xp2(2,j) + ! dz = xp(3,i) - xp2(3,j) + ! dij = (dx*dx)+(dy*dy)+(dz*dz) + ! ... + !----------------------------------------------------------- + ENDDO + ENDDO + ELSE + DO i = 1,Npart + DO j = 1,Lpart + IF (i .EQ. j) CYCLE + !----------------------------------------------------------- + ! USER CODE HERE (without symmetry) + !----------------------------------------------------------- + ENDDO + ENDDO + ENDIF + !------------------------------------------------------------------------- + ! Here we compute the interaction between two different groups + !------------------------------------------------------------------------- + ELSE + IF (lsymm) THEN + DO i = 1,Npart + DO j = 1,Lpart + !----------------------------------------------------------- + ! USER CODE HERE (with symmetry) + !----------------------------------------------------------- + ENDDO + ENDDO + ELSE + DO i = 1,Npart + DO j = 1,Lpart + !----------------------------------------------------------- + ! USER CODE HERE (without symmetry) + !----------------------------------------------------------- + ENDDO + ENDDO + ENDIF + ENDIF + + !------------------------------------------------------------------------- + ! Return + !------------------------------------------------------------------------- + 9999 CONTINUE + RETURN + END SUBROUTINE ppm_template_comp_pp_ring diff --git a/src/ppm_template_comp_pp_verlet.f b/src/ppm_template_comp_pp_verlet.f new file mode 100644 index 0000000..5904db3 --- /dev/null +++ b/src/ppm_template_comp_pp_verlet.f @@ -0,0 +1,303 @@ + !------------------------------------------------------------------------- + ! Subroutine : ppm_template_comp_pp_verlet + !------------------------------------------------------------------------- + ! + ! Purpose : Template for the user subroutine which computes + ! direct particle-particle interactions using Verlet + ! lists. How to use this template: + ! (1) copy this file to your application code dir + ! (2) rename the file and this routine (make sure + ! to also replace all the occurrences of the + ! subroutine name in all WRITE statements). + ! (3) change address in header and add your CVS + ! Log (if needed). + ! (4) declare and add additional arguments + ! (5) add USE and INCLUDE statements for your + ! modules and header files + ! (6) define the precision of floating point data + ! (maybe provide different versions of this + ! subroutine for different precisions) + ! (7) add particle-particle interactions where + ! needed (search for USER CODE). Make sure to + ! replace dimens with the variable which + ! actually contains the problem + ! dimensionality in your code. + ! (8) Optional: Make argument checking conditional + ! on the debug level of your application + ! (9) Optional: Make this routine use your own + ! error logging, message writing, allocation, + ! subroutine start, subroutine stop and + ! timing routines. + ! (10)Optional: Move the variable declarations for + ! vlist and nvlist to your client program and + ! add them as arguments to this routine. This + ! makes it possible to use several sets of + ! lists. + ! (11)Optional: Remove the imode=-1 and imode=1 + ! cases from this routine and build/destroy + ! the lists directly in your client program by + ! calling ppm_neighlist_vlist/DEALLOCATE. You + ! can then remove imode from the argument list + ! of this routine. + ! + ! Input : xp(:,:) (F) particle co-ordinates + ! pdata(:,:) (F) particle data used for interaction + ! (e.g. vorticity, strength, ...) + ! ... (.) USER: add more input arguments as + ! needed + ! Np (I) number of particles on this proc. + ! cutoff (F) cutoff radius for PP interactions + ! skin (F) skin thikness for the Verlet list + ! lsymm (L) use symmetry for PP interactions? + ! .TRUE. : Yes + ! .FALSE.: No + ! imode (I) Mode of action. Any of the folowing: + ! 0 PP interactions + ! 1 build Verlet lists and return + ! -1 destroy Verlet lists and return + ! + ! Input/output : + ! + ! Output : dpd(:,:) (F) Change of particle data (pdata) due to + ! interaction. + ! info (I) return status. =0 if no error. + ! + ! Routines : ppm_neighlist_vlist + ! + ! Remarks : If particles have moved by more than the skin + ! thikness between calls or lsymm has + ! changed, always call with imode=1. + ! + ! seach for USER in this file and fill in + ! the particle-particle interactions in both places. + ! Use the dpd array to return the result. It is your + ! responsibility to properly allocate and initialize + ! this array BEFORE calling this routine. + ! + ! After finishing all time steps, always call this + ! routine with imode=-1 to free the memory of the + ! cell lists. Failure to do so will result in a + ! memory leak in your program. + ! + ! If the CYCLE command is a problem on your hardware + ! (e.g. prevents vectorization), split the DO-loop in + ! two parts as: + ! DO jpart = istart,ipart-1 + ! ... + ! ENDDO + ! DO jpart = ipart+1,iend + ! ... + ! ENDDO + ! + ! References : + ! + ! Revisions : + !------------------------------------------------------------------------- + ! $Revision: 1.1.1.1 $ + ! Revision 1.2 2004/02/24 11:35:54 ivos + ! Revision 1.1 2004/01/26 17:24:35 ivos + !------------------------------------------------------------------------- + ! Parallel Particle Mesh Library (PPM) + ! Institute of Computational Science + ! ETH Zentrum, Hirschengraben 84 + ! CH-8092 Zurich, Switzerland + !------------------------------------------------------------------------- + + SUBROUTINE ppm_template_comp_pp_verlet(xp,pdata,Np,cutoff,skin, & + & lsymm,imode,dpd,info) + + !------------------------------------------------------------------------- + ! Modules + !------------------------------------------------------------------------- + ! USER: USE statements here... + USE ppm_module_neighlist + IMPLICIT NONE + !------------------------------------------------------------------------- + ! Includes + !------------------------------------------------------------------------- + INCLUDE 'ppm_param.h' + !------------------------------------------------------------------------- + ! USER: Define precision + !------------------------------------------------------------------------- +! INTEGER, PARAMETER :: MK = ppm_kind_single + INTEGER, PARAMETER :: MK = ppm_kind_double + !------------------------------------------------------------------------- + ! Arguments + !------------------------------------------------------------------------- + REAL(MK), DIMENSION(:,:), INTENT(IN ) :: xp + REAL(MK), DIMENSION(:,:), INTENT(IN ) :: pdata + ! USER: add more input arguments as needed + INTEGER , INTENT(IN ) :: Np + REAL(MK) , INTENT(IN ) :: cutoff + REAL(MK) , INTENT(IN ) :: skin + LOGICAL , INTENT(IN ) :: lsymm + INTEGER , INTENT(IN ) :: imode + REAL(MK), DIMENSION(:,:), INTENT(INOUT) :: dpd + INTEGER , INTENT( OUT) :: info + !------------------------------------------------------------------------- + ! Local variables + !------------------------------------------------------------------------- + ! Verlet lists + ! USER: if allocation of the following two varaibles fails due to stack + ! size limitations, try putting them in a module and add a USE + ! statement for it above. In order to use several Verlet lists, move + ! these declarations to the client program and add vlist and nvlist + ! to the argument list of this routine. + INTEGER, DIMENSION(:,:), POINTER, SAVE :: vlist + INTEGER, DIMENSION( :), POINTER, SAVE :: nvlist + ! counters + INTEGER :: jpart,ip,jp + ! coordinate differences + REAL(MK) :: dx,dy,dz + ! square of inter particle distance + REAL(MK) :: dij + ! cutoff squared + REAL(MK) :: cut2 + !------------------------------------------------------------------------- + ! Externals + !------------------------------------------------------------------------- + + !------------------------------------------------------------------------- + ! Initialise + !------------------------------------------------------------------------- + info = 0 + cut2 = cutoff*cutoff + + !------------------------------------------------------------------------- + ! Check arguments. + ! USER: You might want to make this conditional on the debug level of + ! your application... + !------------------------------------------------------------------------- + IF (cutoff .LE. 0.0_MK) THEN + WRITE(*,'(2A)') '(ppm_template_comp_pp_verlet): ', & + & 'cutoff must be >0' + info = -1 + GOTO 9999 + ENDIF + IF (skin .LT. 0.0_MK) THEN + WRITE(*,'(2A)') '(ppm_template_comp_pp_verlet): ', & + & 'skin must be >0' + info = -1 + GOTO 9999 + ENDIF + IF (Np .LE. 0) THEN + WRITE(*,'(2A)') '(ppm_template_comp_pp_verlet):', & + & 'Np must be >0' + info = -1 + GOTO 9999 + ENDIF + + !------------------------------------------------------------------------- + ! If imode = -1, destroy Verlet lists and return. + ! USER: This could also be done by your client directly. The code is + ! just here as a template. + !------------------------------------------------------------------------- + IF (imode .EQ. -1) THEN + DEALLOCATE(vlist,nvlist,STAT=info) + IF (info .NE. 0) THEN + WRITE(*,'(2A,I)') '(ppm_template_comp_pp_verlet): ', & + & 'DEALLOCATE failed on line ',__LINE__ + info = -1 + ENDIF + GOTO 9999 + ENDIF + + !------------------------------------------------------------------------- + ! Build new Verlet lists if needed. + ! USER: This could also be done by your client program calling + ! ppm_neighlist_vlist directly. The code is just here as a template. + !------------------------------------------------------------------------- + IF (imode .EQ. 1) THEN + !--------------------------------------------------------------------- + ! Generate Verlet lists + !--------------------------------------------------------------------- + CALL ppm_neighlist_vlist(xp,Np,cutoff,skin,lsymm,vlist,nvlist,info) + IF (info .NE. 0) THEN + WRITE(*,'(2A,I)') '(ppm_template_comp_pp_verlet): ',& + & 'Building Verlet lists failed on line ',__LINE__ + info = -1 + ENDIF + GOTO 9999 + ENDIF + + !------------------------------------------------------------------------- + ! PARTICLE-PARTICLE INTERACTIONS using symmetry + !------------------------------------------------------------------------- + IF (lsymm) THEN + DO ip=1,Np + DO jpart=1,nvlist(ip) + jp = vlist(jpart,ip) + !------------------------------------------------------------- + ! Calculate the square of the distance between the two + ! particles. It will always be .LE. (cutoff+skin)**2 by + ! construction of the Verlet list. + ! COMMENT THIS IF YOU DO NOT NEED THE DISTANCE! + !------------------------------------------------------------- + dx = xp(1,ip) - xp(1,jp) + dy = xp(2,ip) - xp(2,jp) + IF (dimens .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 + ! skip this interaction if the particles are further + ! apart than the given cutoff + IF (dij .GT. cut2) CYCLE + + !------------------------------------------------------------- + ! Particle ip interacts with particle jp here... and + ! vice versa to use symmetry. + !------------------------------------------------------------- + ! + ! USER CODE HERE: + ! dpd(:,ip) = f(pdata(:,ip),pdata(:,jp),dx,dy,dz) + ! dpd(:,jp) = -dpd(:,ip) + ENDDO + ENDDO + + !------------------------------------------------------------------------- + ! PARTICLE-PARTICLE INTERACTIONS not using symmetry + ! This will not vectorize and has therefore been moved into a + ! separate loop !! + !------------------------------------------------------------------------- + ELSE + DO ip=1,Np + DO jpart=1,nvlist(ip) + jp = vlist(jpart,ip) + !------------------------------------------------------------- + ! Calculate the square of the distance between the two + ! particles. It will always be .LE. (cutoff+skin)**2 by + ! construction of the Verlet list. + ! COMMENT THIS IF YOU DO NOT NEED THE DISTANCE! + !------------------------------------------------------------- + dx = xp(1,ip) - xp(1,jp) + dy = xp(2,ip) - xp(2,jp) + IF (dimens .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 + ! skip this interaction if the particles are further + ! apart than the given cutoff + IF (dij .GT. cut2) CYCLE + + !------------------------------------------------------------- + ! Particle ip interacts with particle jp. + !------------------------------------------------------------- + ! + ! USER CODE HERE: + ! dpd(:,ip) = f(pdata(:,ip),pdata(:,jp),dx,dy,dz + ENDDO + ENDDO + ENDIF + + !------------------------------------------------------------------------- + ! Return + !------------------------------------------------------------------------- + 9999 CONTINUE + RETURN + END SUBROUTINE ppm_template_comp_pp_verlet -- GitLab