diff --git a/src/genepsi3d.f90 b/src/genepsi3d.f90 index 47a98d83..ee3c9090 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 @@ -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) @@ -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 @@ -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 @@ -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 @@ -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 @@ -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 @@ -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 @@ -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 @@ -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 @@ -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 @@ -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 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 92b4fdab..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,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