Skip to content

bug fix for particle with periodic condition #345

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
mathrack merged 6 commits into from
Mar 20, 2025
Merged

bug fix for particle with periodic condition #345

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
6 commits
Select commit Hold shift + click to select a range
56b8272
bug fix for particle with periodic condition
fangjian1119 Mar 11, 2025
3c807dc
add lenpe define
fangjian1119 Mar 11, 2025
6c123ad
remove iface
fangjian19 Mar 12, 2025
a247b03
add periodic_bc to decomp_2d_init
fangjian19 Mar 12, 2025
9cb09e2
Define periodic_bc in one line
mathrack Mar 12, 2025
6cb0672
add a comment and format change
fangjian19 Mar 20, 2025
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
2 changes: 2 additions & 0 deletions src/module_param.f90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -538,6 +538,8 @@ module param
real(mytype),parameter :: twopi=two*acos(-one)
#endif

logical :: periodic_bc(3)

end module param
!############################################################################
!############################################################################
Expand Down
5 changes: 5 additions & 0 deletions src/parameters.f90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -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
26 changes: 14 additions & 12 deletions src/particle.f90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Expand All @@ -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)
Expand All @@ -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)
Expand Down Expand Up @@ -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

Expand All @@ -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
Comment on lines +2105 to +2110
Copy link
Member

@pbartholomew08 pbartholomew08 Mar 20, 2025

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Why is the length being defined as 2xL in each direction?

Copy link
Member

@pbartholomew08 pbartholomew08 Mar 20, 2025

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

perhaps add a comment to explain this


firstcal=.false.
endif

Expand All @@ -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
Expand Down
4 changes: 3 additions & 1 deletion src/xcompact3d.f90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -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
Expand Down