bug fix for particle with periodic condition (#345)

* bug fix for particle with periodic condition

* add lenpe define

* remove iface

* add periodic_bc to decomp_2d_init

* Define periodic_bc in one line

* add a comment and format change

---------

Co-authored-by: Jian Fang <jian.fang@stfc.ac.uk>
Co-authored-by: CFLAG <mathrack@users.noreply.github.com>

Author: 3 people

src/module_param.f90

@@ -538,6 +538,8 @@ module param
   real(mytype),parameter :: twopi=two*acos(-one)
 #endif
 
+  logical :: periodic_bc(3)
+  
 end module param
 !############################################################################
 !############################################################################

src/parameters.f90

@@ -137,6 +137,9 @@ subroutine parameter(input_i3d)
   nclxS1 = nclx1; nclxSn = nclxn
   nclyS1 = ncly1; nclySn = nclyn
   nclzS1 = nclz1; nclzSn = nclzn
+
+  !! set periodic direction
+  periodic_bc = (/nclx1.eq.0, ncly1.eq.0, nclz1.eq.0/)
 
   if (numscalar.ne.0) then
      iscalar = 1
@@ -842,4 +845,6 @@ subroutine parameter_defaults()
   ts_tr_tbl=1.402033_mytype
   x0_tr_tbl=3.505082_mytype
 
+  periodic_bc=.false.
+
 end subroutine parameter_defaults

src/particle.f90

@@ -1990,7 +1990,7 @@ subroutine particle_domain_check(particle,num_active_particles)
 
       if(npart>npexit) exit
 
-      if(pa%x(1)<0) then
+      if(pa%x(1)<0._mytype) then
 
         ! xmin face
         call particle_bc(face=1,bctype=bc_particle(1),particle=pa,particle_deduce=counter)
@@ -2002,7 +2002,7 @@ subroutine particle_domain_check(particle,num_active_particles)
 
       endif
 
-      if(pa%x(2)<0) then
+      if(pa%x(2)<0._mytype) then
 
         ! ymin face
         call particle_bc(face=3,bctype=bc_particle(3),particle=pa,particle_deduce=counter)
@@ -2014,7 +2014,7 @@ subroutine particle_domain_check(particle,num_active_particles)
 
       endif
 
-      if(pa%x(3)<0) then
+      if(pa%x(3)<0._mytype) then
 
         ! zmin face
         call particle_bc(face=5,bctype=bc_particle(5),particle=pa,particle_deduce=counter)
@@ -2087,8 +2087,7 @@ subroutine particle_bc(face,bctype,particle,particle_deduce)
     integer,intent(inout) :: particle_deduce
 
     ! local data
-    real(mytype) :: iface
-    real(mytype),save :: bcord(6)
+    real(mytype), dimension(6), save :: bcord, lenpe
     logical,save :: firstcal=.true.
     integer :: idir
 
@@ -2101,6 +2100,15 @@ subroutine particle_bc(face,bctype,particle,particle_deduce)
         bcord(5)=0.0_mytype
         bcord(6)=zlz
 
+        ! lenpe is to get the particle back to the domain for periodic boundaries.
+        ! defined as the distance (+ or -) between the two paring periodic boundaries.
+        lenpe(1)=xlx
+        lenpe(2)=-xlx
+        lenpe(3)=yly
+        lenpe(4)=-yly
+        lenpe(5)=zlz
+        lenpe(6)=-zlz
+
         firstcal=.false.
     endif
 
@@ -2111,14 +2119,8 @@ subroutine particle_bc(face,bctype,particle,particle_deduce)
         call decomp_2d_abort(1,"idir error @ particle_bc")
     endif
 
-    if(mod(face,2)==0) then
-        iface=-1._mytype
-    else
-        iface= 1._mytype
-    endif
-
     if(bctype=='periodic') then
-        particle%x(idir)=particle%x(idir)+bcord(face)*iface
+        particle%x(idir)=particle%x(idir)+lenpe(face)
     elseif(bctype=='reflective') then
         particle%x(idir)=-particle%x(idir)+2.0_mytype*bcord(face)
     elseif(bctype=='outflow') then

src/xcompact3d.f90

@@ -126,6 +126,7 @@ subroutine init_xcompact3d()
 
   use param, only : ilesmod, jles,itype
   use param, only : irestart, mhd_active
+  use param, only : periodic_bc
 
   use variables, only : nx, ny, nz, nxm, nym, nzm
   use variables, only : p_row, p_col
@@ -187,7 +188,8 @@ subroutine init_xcompact3d()
 
   call parameter(InputFN)
 
-  call decomp_2d_init(nx,ny,nz,p_row,p_col)
+  call decomp_2d_init(nx,ny,nz,p_row,p_col,periodic_bc)
+
   call decomp_2d_io_init()
   call init_coarser_mesh_statS(nstat,nstat,nstat,.true.)    !start from 1 == true
   call init_coarser_mesh_statV(nvisu,nvisu,nvisu,.true.)    !start from 1 == true