Skip to content

Avoid huge 3D arrays when iibm = 2 #402

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

Open
mathrack wants to merge 6 commits into
base: master
Choose a base branch
from
Open

Avoid huge 3D arrays when iibm = 2 #402

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
8413cd5
Avoid huge 3D arrays when iibm = 2
mathrack Jan 22, 2026
5598e83
xepsi, yepsi and zepsi are local variables
mathrack Jan 24, 2026
137a1e0
Reduce the number of calls to geomcomplex
mathrack Jan 28, 2026
5f0dde1
Fuse loop in x pencil
mathrack Jan 28, 2026
a325ead
Fuse loop in y pencil
mathrack Jan 28, 2026
1ac9e4f
Fuse loop in z pencil
mathrack Jan 28, 2026
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
154 changes: 57 additions & 97 deletions src/genepsi3d.f90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -65,6 +65,8 @@ subroutine geomcomplex(epsi, nxi, nxf, ny, nyi, nyf, nzi, nzf, dx, yp, dz, remp)
REAL(mytype),DIMENSION(ny) :: yp
REAL(mytype) :: remp

epsi = 0._mytype

IF (itype.EQ.itype_cyl) THEN

CALL geomcomplex_cyl(epsi, nxi, nxf, ny, nyi, nyf, nzi, nzf, dx, yp, remp)
Expand Down Expand Up @@ -159,7 +161,6 @@ subroutine gene_epsi_3D(ep1,nx,ny,nz,dx,dy,dz,xlx,yly,zlz ,&
use var, only : ta2, ta3
use decomp_2d
use MPI
use complex_geometry, only: xepsi, yepsi, zepsi
implicit none
!
real(mytype),dimension(xsize(1),xsize(2),xsize(3)), intent(inout):: ep1
Expand All @@ -184,9 +185,9 @@ subroutine gene_epsi_3D(ep1,nx,ny,nz,dx,dy,dz,xlx,yly,zlz ,&
integer, dimension(xsize(2),xsize(3)) :: nobjxraf
integer, dimension(ysize(1),ysize(3)) :: nobjyraf
integer, dimension(zsize(1),zsize(2)) :: nobjzraf
!real(mytype),dimension(nxraf,xsize(2),xsize(3)) :: xepsi
!real(mytype),dimension(ysize(1),nyraf,ysize(3)) :: yepsi
!real(mytype),dimension(zsize(1),zsize(2),nzraf) :: zepsi
real(mytype),dimension(nxraf,1,1) :: xepsi
real(mytype),dimension(1,nyraf,1) :: yepsi
real(mytype),dimension(1,1,nzraf) :: zepsi
real(mytype),dimension(nyraf) :: ypraf
real(mytype) :: dxraf,dyraf,dzraf
integer :: i,j,k
Expand Down Expand Up @@ -215,9 +216,7 @@ subroutine gene_epsi_3D(ep1,nx,ny,nz,dx,dy,dz,xlx,yly,zlz ,&
else
dxraf =xlx/real(nxraf-1, mytype)
endif
xepsi=zero
call geomcomplex(xepsi,1,nxraf,ny,xstart(2),xend(2),xstart(3),xend(3),dxraf,yp,dz,one)
! if (nrank==0) print*,' step 2'

!y-pencil
if(ncly)then
dyraf =yly/real(nyraf, mytype)
Expand All @@ -234,18 +233,12 @@ subroutine gene_epsi_3D(ep1,nx,ny,nz,dx,dy,dz,xlx,yly,zlz ,&
call stretching(nyraf, ypraf, opt_write = .false.)
end if

yepsi=zero
call geomcomplex(yepsi,ystart(1),yend(1),nyraf,1,nyraf,ystart(3),yend(3),dx,ypraf,dz,one)
! if (nrank==0) print*,' step 3'
!z-pencil
if(nclz)then
dzraf=zlz/real(nzraf, mytype)
else
dzraf=zlz/real(nzraf-1, mytype)
endif
zepsi=zero
call geomcomplex(zepsi,zstart(1),zend(1),ny,zstart(2),zend(2),1,nzraf,dx,yp,dzraf,one)
! if (nrank==0) print*,' step 4'

!x-pencil
nobjx(:,:)=0
Expand Down Expand Up @@ -277,17 +270,27 @@ subroutine gene_epsi_3D(ep1,nx,ny,nz,dx,dy,dz,xlx,yly,zlz ,&
inum=0
do k=1,xsize(3)
do j=1,xsize(2)
! Generate 1D profile on the refined grid
call geomcomplex(xepsi,1,nxraf,ny,xstart(2)+j-1,xstart(2)+j-1,xstart(3)+k-1,xstart(3)+k-1,dxraf,yp,dz,one)
! Update nobjxraf, nobjxmaxraf, xi, xf and ibug
inum=0
if(xepsi(1,j,k).eq.1.)then
if(xepsi(1,1,1) == one)then
inum=1
nobjxraf(j,k)=1
xi(inum,j,k)=-dx!-xlx
endif
do i=1,nxraf-1
if(xepsi(i,j,k).eq.zero.and.xepsi(i+1,j,k).eq.one)then
if(xepsi(i,1,1) == zero.and.xepsi(i+1,1,1) == one)then
inum=inum+1
nobjxraf(j,k)=nobjxraf(j,k)+1
xi(inum,j,k)=dxraf*(i-1)+dxraf/2.
elseif(xepsi(i,1,1) == one .and. xepsi(i+1,1,1) == zero)then
xf(inum,j,k)=dxraf*(i-1)+dxraf/2.
endif
enddo
if(xepsi(nxraf,1,1) == one)then
xf(inum,j,k)=xlx+dx!2.*xlx
end if
if(inum.gt.nobjxmaxraf)then
nobjxmaxraf=inum
endif
Expand Down Expand Up @@ -333,17 +336,27 @@ subroutine gene_epsi_3D(ep1,nx,ny,nz,dx,dy,dz,xlx,yly,zlz ,&
jnum=0
do k=1,ysize(3)
do i=1,ysize(1)
! Generate 1D profile on the refined grid
call geomcomplex(yepsi,ystart(1)+i-1,ystart(1)+i-1,nyraf,1,nyraf,ystart(3)+k-1,ystart(3)+k-1,dx,ypraf,dz,one)
! Update nobjyraf, nobjymaxraf, yi, yf and jbug
jnum=0
if(yepsi(i,1,k) == one)then
if(yepsi(1,1,1) == one)then
jnum=1
nobjyraf(i,k)=1
yi(jnum,i,k)=-(yp(2)-yp(1))!-yly
endif
do j=1,nyraf-1
if(yepsi(i,j,k) == zero .and. yepsi(i,j+1,k) == one)then
if(yepsi(1,j,1) == zero .and. yepsi(1,j+1,1) == one)then
jnum=jnum+1
nobjyraf(i,k)=nobjyraf(i,k)+1
yi(jnum,i,k)=ypraf(j)+(ypraf(j+1)-ypraf(j))*half!dyraf*(j-1)+dyraf/2.
elseif(yepsi(1,j,1) == one .and. yepsi(1,j+1,1) == zero)then
yf(jnum,i,k)=ypraf(j)+(ypraf(j+1)-ypraf(j))*half!dyraf*(j-1)+dyraf/2.
endif
enddo
if(yepsi(1,nyraf,1) == one)then
yf(jnum,i,k)=yly+(yp(ny)-yp(ny-1))*half!2.*yly
endif
if(jnum.gt.nobjymaxraf)then
nobjymaxraf=jnum
endif
Expand Down Expand Up @@ -389,17 +402,27 @@ subroutine gene_epsi_3D(ep1,nx,ny,nz,dx,dy,dz,xlx,yly,zlz ,&
knum=0
do j=1,zsize(2)
do i=1,zsize(1)
! Generate 1D profile on the refined grid
call geomcomplex(zepsi,zstart(1)+i-1,zstart(1)+i-1,ny,zstart(2)+j-1,zstart(2)+j-1,1,nzraf,dx,yp,dzraf,one)
! Update nobjzraf, nobjzmaxraf, zi, zf and kbug
knum=0
if(zepsi(i,j,1) == one)then
if(zepsi(1,1,1) == one)then
knum=1
nobjzraf(i,j)=1
zi(knum,i,j)=-dz!zlz
endif
do k=1,nzraf-1
if(zepsi(i,j,k) == zero .and. zepsi(i,j,k+1) == one)then
if(zepsi(1,1,k) == zero .and. zepsi(1,1,k+1) == one)then
knum=knum+1
nobjzraf(i,j)=nobjzraf(i,j)+1
zi(knum,i,j)=dzraf*(k-1)+dzraf*half
elseif(zepsi(1,1,k) == one .and. zepsi(1,1,k+1) == zero)then
zf(knum,i,j)=dzraf*(k-1)+dzraf*half
endif
enddo
if(zepsi(1,1,nzraf) == one)then
zf(knum,i,j)=zlz+dz!2.*zlz
endif
if(knum.gt.nobjzmaxraf)then
nobjzmaxraf=knum
endif
Expand All @@ -414,31 +437,10 @@ subroutine gene_epsi_3D(ep1,nx,ny,nz,dx,dy,dz,xlx,yly,zlz ,&
! if (nrank==0) print*,' kbug=',mpi_aux_i
! if (nrank==0) print*,' step 7'

!x-pencil
do k=1,xsize(3)
do j=1,xsize(2)
inum=0
if(xepsi(1,j,k) == one)then
inum=inum+1
xi(inum,j,k)=-dx!-xlx
endif
do i=1,nxraf-1
if(xepsi(i,j,k) == zero .and. xepsi(i+1,j,k) == one)then
inum=inum+1
xi(inum,j,k)=dxraf*(i-1)+dxraf/2.
elseif(xepsi(i,j,k) == one .and. xepsi(i+1,j,k)== zero)then
xf(inum,j,k)=dxraf*(i-1)+dxraf/2.
endif
enddo
if(xepsi(nxraf,j,k)==1.)then
xf(inum,j,k)=xlx+dx!2.*xlx
endif
enddo
enddo

if(ibug /= 0)then
do k=1,xsize(3)
do j=1,xsize(2)
call geomcomplex(xepsi,1,nxraf,ny,xstart(2)+j-1,xstart(2)+j-1,xstart(3)+k-1,xstart(3)+k-1,dxraf,yp,dz,one)
if(nobjx(j,k) /= nobjxraf(j,k))then
iobj=0
if(ep1(1,j,k) == one)iobj=iobj+1
Expand All @@ -447,10 +449,10 @@ subroutine gene_epsi_3D(ep1,nx,ny,nz,dx,dy,dz,xlx,yly,zlz ,&
if(ep1(i,j,k) == zero .and. ep1(i+1,j,k) == zero)iflu=1
if(ep1(i,j,k) == one .and. ep1(i+1,j,k) == one)isol=1
do iraf=1,nraf
if(xepsi(iraf+nraf*(i-1) ,j,k) == zero .and.&
xepsi(iraf+nraf*(i-1)+1,j,k) == one)idebraf=iraf+nraf*(i-1)+1
if(xepsi(iraf+nraf*(i-1) ,j,k) == one .and.&
xepsi(iraf+nraf*(i-1)+1,j,k) == zero)ifinraf=iraf+nraf*(i-1)+1
if(xepsi(iraf+nraf*(i-1) ,1,1) == zero .and.&
xepsi(iraf+nraf*(i-1)+1,1,1) == one)idebraf=iraf+nraf*(i-1)+1
if(xepsi(iraf+nraf*(i-1) ,1,1) == one .and.&
xepsi(iraf+nraf*(i-1)+1,1,1) == zero)ifinraf=iraf+nraf*(i-1)+1
enddo
if(idebraf /= 0 .and. ifinraf /= 0 .and.&
idebraf < ifinraf .and. iflu == 1)then
Expand Down Expand Up @@ -482,31 +484,10 @@ subroutine gene_epsi_3D(ep1,nx,ny,nz,dx,dy,dz,xlx,yly,zlz ,&
endif
!if (nrank==0) write(*,*) ' step 8'

!y-pencil
do k=1,ysize(3)
do i=1,ysize(1)
jnum=0
if(yepsi(i,1,k) == one)then
jnum=jnum+1
yi(jnum,i,k)=-(yp(2)-yp(1))!-yly
endif
do j=1,nyraf-1
if(yepsi(i,j,k) == zero .and. yepsi(i,j+1,k) == one)then
jnum=jnum+1
yi(jnum,i,k)=ypraf(j)+(ypraf(j+1)-ypraf(j))*half!dyraf*(j-1)+dyraf/2.
elseif(yepsi(i,j,k) == one .and. yepsi(i,j+1,k) == zero)then
yf(jnum,i,k)=ypraf(j)+(ypraf(j+1)-ypraf(j))*half!dyraf*(j-1)+dyraf/2.
endif
enddo
if(yepsi(i,nyraf,k) == one)then
yf(jnum,i,k)=yly+(yp(ny)-yp(ny-1))*half!2.*yly
endif
enddo
enddo

if(jbug /= 0)then
do k=1,ysize(3)
do i=1,ysize(1)
call geomcomplex(yepsi,ystart(1)+i-1,ystart(1)+i-1,nyraf,1,nyraf,ystart(3)+k-1,ystart(3)+k-1,dx,ypraf,dz,one)
if(nobjy(i,k) /= nobjyraf(i,k))then
jobj=0
if(ta2(i,1,k) == one)jobj=jobj+1
Expand All @@ -515,10 +496,10 @@ subroutine gene_epsi_3D(ep1,nx,ny,nz,dx,dy,dz,xlx,yly,zlz ,&
if(ta2(i,j,k) == zero .and. ta2(i,j+1,k) == zero)jflu=1
if(ta2(i,j,k) == one .and. ta2(i,j+1,k) == one)jsol=1
do jraf=1,nraf
if(yepsi(i,jraf+nraf*(j-1) ,k) == zero .and.&
yepsi(i,jraf+nraf*(j-1)+1,k) == one)jdebraf=jraf+nraf*(j-1)+1
if(yepsi(i,jraf+nraf*(j-1) ,k) == one .and.&
yepsi(i,jraf+nraf*(j-1)+1,k) == zero)jfinraf=jraf+nraf*(j-1)+1
if(yepsi(1,jraf+nraf*(j-1) ,1) == zero .and.&
yepsi(1,jraf+nraf*(j-1)+1,1) == one)jdebraf=jraf+nraf*(j-1)+1
if(yepsi(1,jraf+nraf*(j-1) ,1) == one .and.&
yepsi(1,jraf+nraf*(j-1)+1,1) == zero)jfinraf=jraf+nraf*(j-1)+1
enddo
if(jdebraf /= 0 .and. jfinraf /= 0 .and.&
jdebraf < jfinraf.and.jflu == 1)then
Expand Down Expand Up @@ -550,32 +531,11 @@ subroutine gene_epsi_3D(ep1,nx,ny,nz,dx,dy,dz,xlx,yly,zlz ,&
endif
!if (nrank==0) write(*,*) ' step 9'

!z-pencil
do j=1,zsize(2)
do i=1,zsize(1)
knum=0
if(zepsi(i,j,1) == one)then
knum=knum+1
zi(knum,i,j)=-dz!zlz
endif
do k=1,nzraf-1
if(zepsi(i,j,k) == zero .and. zepsi(i,j,k+1) == one)then
knum=knum+1
zi(knum,i,j)=dzraf*(k-1)+dzraf*half
elseif(zepsi(i,j,k) == one .and. zepsi(i,j,k+1) == zero)then
zf(knum,i,j)=dzraf*(k-1)+dzraf*half
endif
enddo
if(zepsi(i,j,nzraf) == one)then
zf(knum,i,j)=zlz+dz!2.*zlz
endif
enddo
enddo

kdebraf=0
if(kbug.ne.0)then
do j=1,zsize(2)
do i=1,zsize(1)
call geomcomplex(zepsi,zstart(1)+i-1,zstart(1)+i-1,ny,zstart(2)+j-1,zstart(2)+j-1,1,nzraf,dx,yp,dzraf,one)
if(nobjz(i,j) /= nobjzraf(i,j))then
kobj=0
if(ta3(i,j,1) == one)kobj=kobj+1
Expand All @@ -584,10 +544,10 @@ subroutine gene_epsi_3D(ep1,nx,ny,nz,dx,dy,dz,xlx,yly,zlz ,&
if(ta3(i,j,k) == zero .and. ta3(i,j,k+1) == zero)kflu=1
if(ta3(i,j,k) == one .and. ta3(i,j,k+1) == one)ksol=1
do kraf=1,nraf
if(zepsi(i,j,kraf+nraf*(k-1) ) == zero .and.&
zepsi(i,j,kraf+nraf*(k-1)+1) == one)kdebraf=kraf+nraf*(k-1)+1
if(zepsi(i,j,kraf+nraf*(k-1) ) == one .and.&
zepsi(i,j,kraf+nraf*(k-1)+1) == zero)kfinraf=kraf+nraf*(k-1)+1
if(zepsi(1,1,kraf+nraf*(k-1) ) == zero .and.&
zepsi(1,1,kraf+nraf*(k-1)+1) == one)kdebraf=kraf+nraf*(k-1)+1
if(zepsi(1,1,kraf+nraf*(k-1) ) == one .and.&
zepsi(1,1,kraf+nraf*(k-1)+1) == zero)kfinraf=kraf+nraf*(k-1)+1
enddo
if(kdebraf /= 0 .and. kfinraf /= 0 .and.&
kdebraf < kfinraf .and. kflu == 1)then
Expand Down
1 change: 0 additions & 1 deletion src/module_param.f90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -551,7 +551,6 @@ module complex_geometry
integer ,allocatable,dimension(:,:) :: nobjx,nobjy,nobjz
integer ,allocatable,dimension(:,:,:) :: nxipif,nxfpif,nyipif,nyfpif,nzipif,nzfpif
real(mytype),allocatable,dimension(:,:,:) :: xi,xf,yi,yf,zi,zf
real(mytype),allocatable,dimension(:,:,:) :: xepsi, yepsi, zepsi
integer :: nxraf,nyraf,nzraf,nraf,nobjmax
end module complex_geometry
!############################################################################
Expand Down
6 changes: 0 additions & 6 deletions src/variables.f90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -656,12 +656,6 @@ subroutine init_variables
zi=zero
allocate(zf(nobjmax,zsize(1),zsize(2)))
zf=zero
allocate(xepsi(nxraf,xsize(2),xsize(3)))
xepsi=zero
allocate(yepsi(ysize(1),nyraf,ysize(3)))
yepsi=zero
allocate(zepsi(zsize(1),zsize(2),nzraf))
zepsi=zero

endif

Expand Down