From 8413cd55cedd6c22bea42b190afc776c3d12d7d4 Mon Sep 17 00:00:00 2001 From: cflag Date: Thu, 22 Jan 2026 19:41:49 +0100 Subject: [PATCH 1/6] Avoid huge 3D arrays when iibm = 2 --- src/genepsi3d.f90 | 79 ++++++++++++++++++++++++----------------------- src/variables.f90 | 6 ++-- 2 files changed, 43 insertions(+), 42 deletions(-) diff --git a/src/genepsi3d.f90 b/src/genepsi3d.f90 index 47a98d83..1b1cb9b5 100644 --- a/src/genepsi3d.f90 +++ b/src/genepsi3d.f90 @@ -215,9 +215,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) @@ -234,18 +232,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 @@ -277,13 +269,14 @@ 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) + 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) inum=0 - if(xepsi(1,j,k).eq.1.)then + if(xepsi(1,1,1).eq.1.)then inum=1 nobjxraf(j,k)=1 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).eq.zero.and.xepsi(i+1,1,1).eq.one)then inum=inum+1 nobjxraf(j,k)=nobjxraf(j,k)+1 endif @@ -333,13 +326,14 @@ 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) + 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) jnum=0 - if(yepsi(i,1,k) == one)then + if(yepsi(1,1,1) == one)then jnum=1 nobjyraf(i,k)=1 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 endif @@ -389,13 +383,14 @@ 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) + 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) knum=0 - if(zepsi(i,j,1) == one)then + if(zepsi(1,1,1) == one)then knum=1 nobjzraf(i,j)=1 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 endif @@ -417,20 +412,21 @@ subroutine gene_epsi_3D(ep1,nx,ny,nz,dx,dy,dz,xlx,yly,zlz ,& !x-pencil 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) inum=0 - if(xepsi(1,j,k) == one)then + if(xepsi(1,1,1) == 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 + if(xepsi(i,1,1) == zero .and. xepsi(i+1,1,1) == 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 + 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,j,k)==1.)then + if(xepsi(nxraf,1,1)==1.)then xf(inum,j,k)=xlx+dx!2.*xlx endif enddo @@ -439,6 +435,7 @@ subroutine gene_epsi_3D(ep1,nx,ny,nz,dx,dy,dz,xlx,yly,zlz ,& 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 @@ -447,10 +444,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 @@ -485,20 +482,21 @@ subroutine gene_epsi_3D(ep1,nx,ny,nz,dx,dy,dz,xlx,yly,zlz ,& !y-pencil 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) jnum=0 - if(yepsi(i,1,k) == one)then + if(yepsi(1,1,1) == 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 + if(yepsi(1,j,1) == zero .and. yepsi(1,j+1,1) == 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 + 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(i,nyraf,k) == one)then + if(yepsi(1,nyraf,1) == one)then yf(jnum,i,k)=yly+(yp(ny)-yp(ny-1))*half!2.*yly endif enddo @@ -507,6 +505,7 @@ subroutine gene_epsi_3D(ep1,nx,ny,nz,dx,dy,dz,xlx,yly,zlz ,& 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 @@ -515,10 +514,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 @@ -553,20 +552,21 @@ subroutine gene_epsi_3D(ep1,nx,ny,nz,dx,dy,dz,xlx,yly,zlz ,& !z-pencil 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) knum=0 - if(zepsi(i,j,1) == one)then + if(zepsi(1,1,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 + if(zepsi(1,1,k) == zero .and. zepsi(1,1,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 + 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(i,j,nzraf) == one)then + if(zepsi(1,1,nzraf) == one)then zf(knum,i,j)=zlz+dz!2.*zlz endif enddo @@ -576,6 +576,7 @@ subroutine gene_epsi_3D(ep1,nx,ny,nz,dx,dy,dz,xlx,yly,zlz ,& 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 @@ -584,10 +585,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 diff --git a/src/variables.f90 b/src/variables.f90 index 92b4fdab..0e4e40bb 100644 --- a/src/variables.f90 +++ b/src/variables.f90 @@ -656,11 +656,11 @@ subroutine init_variables zi=zero allocate(zf(nobjmax,zsize(1),zsize(2))) zf=zero - allocate(xepsi(nxraf,xsize(2),xsize(3))) + allocate(xepsi(nxraf,1,1)) xepsi=zero - allocate(yepsi(ysize(1),nyraf,ysize(3))) + allocate(yepsi(1,nyraf,1)) yepsi=zero - allocate(zepsi(zsize(1),zsize(2),nzraf)) + allocate(zepsi(1,1,nzraf)) zepsi=zero endif From 5598e835afaab9d2fb7a0afc026d42ef1e841919 Mon Sep 17 00:00:00 2001 From: cflag Date: Sat, 24 Jan 2026 13:07:21 +0100 Subject: [PATCH 2/6] xepsi, yepsi and zepsi are local variables --- src/genepsi3d.f90 | 9 +++++---- src/module_param.f90 | 1 - src/variables.f90 | 6 ------ 3 files changed, 5 insertions(+), 11 deletions(-) diff --git a/src/genepsi3d.f90 b/src/genepsi3d.f90 index 1b1cb9b5..6073c124 100644 --- a/src/genepsi3d.f90 +++ b/src/genepsi3d.f90 @@ -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) @@ -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 @@ -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 diff --git a/src/module_param.f90 b/src/module_param.f90 index 8ca94d3b..92b83ee9 100644 --- a/src/module_param.f90 +++ b/src/module_param.f90 @@ -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 !############################################################################ diff --git a/src/variables.f90 b/src/variables.f90 index 0e4e40bb..853d4408 100644 --- a/src/variables.f90 +++ b/src/variables.f90 @@ -656,12 +656,6 @@ subroutine init_variables zi=zero allocate(zf(nobjmax,zsize(1),zsize(2))) zf=zero - allocate(xepsi(nxraf,1,1)) - xepsi=zero - allocate(yepsi(1,nyraf,1)) - yepsi=zero - allocate(zepsi(1,1,nzraf)) - zepsi=zero endif From 137a1e034730633a60d2501c51f18c9fd4d5e775 Mon Sep 17 00:00:00 2001 From: cflag Date: Wed, 28 Jan 2026 08:12:11 +0100 Subject: [PATCH 3/6] Reduce the number of calls to geomcomplex --- src/genepsi3d.f90 | 126 +++++++++++++++++++++------------------------- 1 file changed, 57 insertions(+), 69 deletions(-) diff --git a/src/genepsi3d.f90 b/src/genepsi3d.f90 index 6073c124..2cc52c76 100644 --- a/src/genepsi3d.f90 +++ b/src/genepsi3d.f90 @@ -270,7 +270,9 @@ 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 and ibug inum=0 if(xepsi(1,1,1).eq.1.)then inum=1 @@ -288,6 +290,23 @@ subroutine gene_epsi_3D(ep1,nx,ny,nz,dx,dy,dz,xlx,yly,zlz ,& if(nobjx(j,k).ne.nobjxraf(j,k))then ibug=ibug+1 endif + ! Update xi and xf + inum=0 + if(xepsi(1,1,1) == one)then + inum=inum+1 + xi(inum,j,k)=-dx!-xlx + endif + do i=1,nxraf-1 + if(xepsi(i,1,1) == zero .and. xepsi(i+1,1,1) == one)then + inum=inum+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)==1.)then + xf(inum,j,k)=xlx+dx!2.*xlx + endif enddo enddo call MPI_REDUCE(nobjxmaxraf,mpi_aux_i,1,MPI_INTEGER,MPI_MAX,0,MPI_COMM_WORLD,code) @@ -327,7 +346,9 @@ 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 and jbug jnum=0 if(yepsi(1,1,1) == one)then jnum=1 @@ -345,6 +366,23 @@ subroutine gene_epsi_3D(ep1,nx,ny,nz,dx,dy,dz,xlx,yly,zlz ,& if(nobjy(i,k).ne.nobjyraf(i,k))then jbug=jbug+1 endif + ! Update yi and yf + jnum=0 + if(yepsi(1,1,1) == one)then + jnum=jnum+1 + yi(jnum,i,k)=-(yp(2)-yp(1))!-yly + endif + do j=1,nyraf-1 + if(yepsi(1,j,1) == zero .and. yepsi(1,j+1,1) == one)then + jnum=jnum+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 enddo enddo call MPI_REDUCE(nobjymaxraf,mpi_aux_i,1,MPI_INTEGER,MPI_MAX,0,MPI_COMM_WORLD,code) @@ -384,7 +422,9 @@ 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 and kbug knum=0 if(zepsi(1,1,1) == one)then knum=1 @@ -402,6 +442,23 @@ subroutine gene_epsi_3D(ep1,nx,ny,nz,dx,dy,dz,xlx,yly,zlz ,& if(nobjz(i,j).ne.nobjzraf(i,j))then kbug=kbug+1 endif + ! Update zi and zf + knum=0 + if(zepsi(1,1,1) == one)then + knum=knum+1 + zi(knum,i,j)=-dz!zlz + endif + do k=1,nzraf-1 + if(zepsi(1,1,k) == zero .and. zepsi(1,1,k+1) == one)then + knum=knum+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 enddo enddo call MPI_REDUCE(nobjzmaxraf,mpi_aux_i,1,MPI_INTEGER,MPI_MAX,0,MPI_COMM_WORLD,code) @@ -410,29 +467,6 @@ 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) - 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) - inum=0 - if(xepsi(1,1,1) == one)then - inum=inum+1 - xi(inum,j,k)=-dx!-xlx - endif - do i=1,nxraf-1 - if(xepsi(i,1,1) == zero .and. xepsi(i+1,1,1) == one)then - inum=inum+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)==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) @@ -480,29 +514,6 @@ 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) - 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) - jnum=0 - if(yepsi(1,1,1) == one)then - jnum=jnum+1 - yi(jnum,i,k)=-(yp(2)-yp(1))!-yly - endif - do j=1,nyraf-1 - if(yepsi(1,j,1) == zero .and. yepsi(1,j+1,1) == one)then - jnum=jnum+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 - enddo - enddo - if(jbug /= 0)then do k=1,ysize(3) do i=1,ysize(1) @@ -550,29 +561,6 @@ 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) - 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) - knum=0 - if(zepsi(1,1,1) == one)then - knum=knum+1 - zi(knum,i,j)=-dz!zlz - endif - do k=1,nzraf-1 - if(zepsi(1,1,k) == zero .and. zepsi(1,1,k+1) == one)then - knum=knum+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 - enddo - enddo - kdebraf=0 if(kbug.ne.0)then do j=1,zsize(2) From 5f0dde1a8db247a8637e3383a705b6fc2d3c6c09 Mon Sep 17 00:00:00 2001 From: cflag Date: Wed, 28 Jan 2026 08:17:48 +0100 Subject: [PATCH 4/6] Fuse loop in x pencil --- src/genepsi3d.f90 | 30 ++++++++++-------------------- 1 file changed, 10 insertions(+), 20 deletions(-) diff --git a/src/genepsi3d.f90 b/src/genepsi3d.f90 index 2cc52c76..ff41f567 100644 --- a/src/genepsi3d.f90 +++ b/src/genepsi3d.f90 @@ -272,41 +272,31 @@ subroutine gene_epsi_3D(ep1,nx,ny,nz,dx,dy,dz,xlx,yly,zlz ,& 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 and ibug + ! Update nobjxraf, nobjxmaxraf, xi, xf and ibug inum=0 - if(xepsi(1,1,1).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,1,1).eq.zero.and.xepsi(i+1,1,1).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 if(nobjx(j,k).ne.nobjxraf(j,k))then ibug=ibug+1 endif - ! Update xi and xf - inum=0 - if(xepsi(1,1,1) == one)then - inum=inum+1 - xi(inum,j,k)=-dx!-xlx - endif - do i=1,nxraf-1 - if(xepsi(i,1,1) == zero .and. xepsi(i+1,1,1) == one)then - inum=inum+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)==1.)then - xf(inum,j,k)=xlx+dx!2.*xlx - endif enddo enddo call MPI_REDUCE(nobjxmaxraf,mpi_aux_i,1,MPI_INTEGER,MPI_MAX,0,MPI_COMM_WORLD,code) From a325ead5a08ad833ee6bfdbab5b89ff4d0994113 Mon Sep 17 00:00:00 2001 From: cflag Date: Wed, 28 Jan 2026 08:20:09 +0100 Subject: [PATCH 5/6] Fuse loop in y pencil --- src/genepsi3d.f90 | 26 ++++++++------------------ 1 file changed, 8 insertions(+), 18 deletions(-) diff --git a/src/genepsi3d.f90 b/src/genepsi3d.f90 index ff41f567..e736ddf8 100644 --- a/src/genepsi3d.f90 +++ b/src/genepsi3d.f90 @@ -338,33 +338,17 @@ subroutine gene_epsi_3D(ep1,nx,ny,nz,dx,dy,dz,xlx,yly,zlz ,& 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 and jbug + ! Update nobjyraf, nobjymaxraf, yi, yf and jbug jnum=0 if(yepsi(1,1,1) == one)then jnum=1 nobjyraf(i,k)=1 - endif - do j=1,nyraf-1 - if(yepsi(1,j,1) == zero .and. yepsi(1,j+1,1) == one)then - jnum=jnum+1 - nobjyraf(i,k)=nobjyraf(i,k)+1 - endif - enddo - if(jnum.gt.nobjymaxraf)then - nobjymaxraf=jnum - endif - if(nobjy(i,k).ne.nobjyraf(i,k))then - jbug=jbug+1 - endif - ! Update yi and yf - jnum=0 - if(yepsi(1,1,1) == one)then - jnum=jnum+1 yi(jnum,i,k)=-(yp(2)-yp(1))!-yly endif do j=1,nyraf-1 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. @@ -373,6 +357,12 @@ subroutine gene_epsi_3D(ep1,nx,ny,nz,dx,dy,dz,xlx,yly,zlz ,& 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 + if(nobjy(i,k).ne.nobjyraf(i,k))then + jbug=jbug+1 + endif enddo enddo call MPI_REDUCE(nobjymaxraf,mpi_aux_i,1,MPI_INTEGER,MPI_MAX,0,MPI_COMM_WORLD,code) From 1ac9e4f0eeac034533cccbb6ed546720557e34c6 Mon Sep 17 00:00:00 2001 From: cflag Date: Wed, 28 Jan 2026 08:21:24 +0100 Subject: [PATCH 6/6] Fuse loop in z pencil --- src/genepsi3d.f90 | 26 ++++++++------------------ 1 file changed, 8 insertions(+), 18 deletions(-) diff --git a/src/genepsi3d.f90 b/src/genepsi3d.f90 index e736ddf8..ee3c9090 100644 --- a/src/genepsi3d.f90 +++ b/src/genepsi3d.f90 @@ -404,33 +404,17 @@ subroutine gene_epsi_3D(ep1,nx,ny,nz,dx,dy,dz,xlx,yly,zlz ,& 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 and kbug + ! Update nobjzraf, nobjzmaxraf, zi, zf and kbug knum=0 if(zepsi(1,1,1) == one)then knum=1 nobjzraf(i,j)=1 - endif - do k=1,nzraf-1 - if(zepsi(1,1,k) == zero .and. zepsi(1,1,k+1) == one)then - knum=knum+1 - nobjzraf(i,j)=nobjzraf(i,j)+1 - endif - enddo - if(knum.gt.nobjzmaxraf)then - nobjzmaxraf=knum - endif - if(nobjz(i,j).ne.nobjzraf(i,j))then - kbug=kbug+1 - endif - ! Update zi and zf - knum=0 - if(zepsi(1,1,1) == one)then - knum=knum+1 zi(knum,i,j)=-dz!zlz endif do k=1,nzraf-1 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 @@ -439,6 +423,12 @@ subroutine gene_epsi_3D(ep1,nx,ny,nz,dx,dy,dz,xlx,yly,zlz ,& if(zepsi(1,1,nzraf) == one)then zf(knum,i,j)=zlz+dz!2.*zlz endif + if(knum.gt.nobjzmaxraf)then + nobjzmaxraf=knum + endif + if(nobjz(i,j).ne.nobjzraf(i,j))then + kbug=kbug+1 + endif enddo enddo call MPI_REDUCE(nobjzmaxraf,mpi_aux_i,1,MPI_INTEGER,MPI_MAX,0,MPI_COMM_WORLD,code)