module bounds ! ! (c) Copyright 1991 to 1998 by Michael A. Beer, William D. Dorland, ! P. B. Snyder, Q. P. Liu, and Gregory W. Hammett. ALL RIGHTS RESERVED. ! implicit none save integer ndomz ! # of different FFT lengths. integer nconz ! # of connected FFT domains integer, dimension(:), allocatable :: ndom integer, dimension(:,:), allocatable :: conpos, mcon, icon integer, dimension(:,:), allocatable :: ncon, scon integer, dimension(:,:,:), allocatable :: nncon, mmcon contains subroutine bcon(a) ! ! Enforce boundary conditions in the x (l) direction ! use constants, only: ii use itg_data, only: iperiod, n0, qsf use gryffin_layouts use gryffin_grid, only: ld, md, nd, mrr_all, nrr, l_left, l_right, xp, z0 use gryffin_redist, only: gryffin_fill_type, init_fill_ends use redistribute, only: fill complex, dimension(:,m_low:,n_low:) :: a real :: b1, znshift, wgt1, wgt2, wwidth integer :: l, m, n, nshift, n_max, mr, nr integer, dimension(-nd:nd) :: ninv logical :: first = .true. type (gryffin_fill_type), save :: ends if(iperiod == 0) then a(1,:,:)=0. a(ld,:,:)=0. endif if (iperiod == 1) then a(ld,:,:)=a(1,:,:) endif ! Twisted periodicity if (iperiod == 2) then n_max=(nd-1)/2 ninv = 0 do n=1,nd do nr=-n_max,n_max if (nrr(n) == nr) ninv(nr)=n enddo enddo znshift=2.*n_max*(xp+.0001)/z0 ! handle the special case of a single n mode (nd=1, which ! might be used for linear runs) properly: if(n_max == 0) znshift=1 b1 = 2.*xp*n0*qsf if (first) then call init_fill_ends(ends, mrr_all, nrr, ninv, b1, znshift, & l_left, l_right, ld, n_max, nd) first = .false. endif a(1:l_left-1,:,:) = 0. a(l_right:ld,:,:) = 0. call fill (ends%fill, a, a) ! Window to filter near the edges of the box. ! ! The window weight function is 1 in most of the box, but near the end ! of the extension goes like 2x**2/(1+x**4), where x is the distance ! from the edge of the box, and is normalized to match smoothly to the ! weight=1 region. ! ! Another option might be to use sine transforms... ! ! wwidth is the fraction of the extension region over which the window ! is dropping to zero: wwidth=0.5 ! (Careful, skip m=0 mode, which is supposed to be periodic with itself...). if(md >= 2) then do m=m_low, m_high if(m == 1) cycle do n=n_low, n_high do l=1,l_left-1 ! wgt2 is a linear function, 0 at the end, and reaches a maximum of 1: wgt2=(l-1)/float(l_left-1)/wwidth ! wgt2=1.0 ! to revert to no window if(wgt2 >= 1.0) wgt2=1.0 wgt1=2*wgt2**2/(1+wgt2**4) a(l,m,n)=a(l,m,n)*wgt1 enddo ! a(1,m,n)=0.0 ! a(2,m,n)=a(2,m,n)*0.5 do l=l_right,ld wgt2=(ld-l)/float(ld-l_right)/wwidth ! wgt2=1.0 ! to revert to no window if(wgt2 >= 1.0) wgt2=1.0 wgt1=2*wgt2**2/(1+wgt2**4) a(l,m,n)=a(l,m,n)*wgt1 enddo ! a(ld,m,n)=0.0 ! a(ld-1,m,n)=a(ld-1,m,n)*0.5 enddo enddo endif endif ! ! for connected FFTs, zero leftmost point of leftmost mode ! and rightmost point of rightmost mode, except for m=1 (ky=0) ! if(iperiod == 3) then do m=m_low, m_high if(m == 1) cycle do n=n_low, n_high if (conpos(m,n) == 1) a(1,m,n)=0. if (conpos(m,n) == icon(m,n)) a(ld,m,n)=0. enddo enddo m = 1 do n = n_low, n_high if(idx_local(m_lo, m, n)) a(ld,m,n)=a(1,m,n) enddo endif end subroutine bcon subroutine connect ! ! (c) Copyright 1991 to 1998 by Michael A. Beer, William D. Dorland, ! P. B. Snyder, Q. P. Liu, and Gregory W. Hammett. ALL RIGHTS RESERVED. ! ! ! loads arrays to do ffts on connected domain, ! instead of using extensions ! use itg_data, only: lin, linlay use mp, only: proc0 use gryffin_grid, only: md, nd, nrr, mrr_all, x0, xp, z0 use constants, only: pi use gryffin_layouts integer :: jmax,kmax,jshift,jshift0 integer, allocatable, dimension(:,:) :: l integer :: t, m, n, count, jj, i, j, k, ncon_tmp, m0 logical :: skip_code allocate (l(-nd:nd,0:md-1)) l = 0 allocate (conpos(md, nd), mcon(md, nd), icon(md, nd), ndom(md)) conpos = 0; mcon = 0; icon = 0; ndom = 0 jmax=(nd-1)/2 kmax=md-1 jshift0=(xp+.0001)/z0*(nd-1) if (jshift0 == 0) jshift0=nd**2 if((x0 /= xp) .or. (abs(int((x0+.0001)/pi)*pi-x0) > 0.001) & .or. (abs(int((xp+.0001)/pi)*pi-xp) > 0.001)) then if(proc0) write(*,*) 'failure in connect.f:' if(proc0) write(*,*) ' for iconnect=3, must have x0=xp=pi * integer' call aborter(6,' ') endif ! write(*,*) 'jmax=',jmax,' kmax=',kmax,' jshift(k=1)=',jshift0 ! ! Calculate how many domains are connected by the BC ! do j=-jmax,jmax l(j,0)=1 enddo do k=1,kmax jshift=jshift0*k if (lin == 1 .and. linlay) jshift=jshift0 do j=-jmax,jmax l(j,k)=(jmax-j)/jshift+(j+jmax)/jshift+1 if (j+jmax >= jshift) l(j,k)=0 ! if ((j < 0).and.(l(j,k) > 1)) l(j,k)=0 enddo enddo ! write(*,*) 'jshift(k=1): ',jshift0 ! write(*,100) 0,0,1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16 100 format(i5,17(i3)) ! do j=-jmax,jmax ! write(*,100) j,(l(j,k),k=0,kmax) ! enddo nconz = 1 ndomz = 1 ndom=1 do k=0,kmax do m=1,jmax*2+1 ncon_tmp=0 do j=-jmax,jmax if (l(j,k) == m) then ncon_tmp=ncon_tmp+1 nconz = max(nconz, m) endif enddo if (ncon_tmp > 0) ndom(k+1)=ndom(k+1)+1 enddo ndomz = max(ndomz, ndom(k+1)) enddo allocate (nncon(md, ndomz, md*nd), mmcon(md, ndomz, md*nd)) nncon = 0; mmcon = 0 allocate (ncon(md, ndomz), scon(md, ndomz)) ncon = 0; scon = 0 t=0 do m0=1,md ndom(m0)=1 do m=1,jmax*2+1 ncon(m0,ndom(m0))=0 k=m0-1 do j=-jmax,jmax if (l(j,k) == m) then ncon(m0,ndom(m0))=ncon(m0,ndom(m0))+1 scon(m0,ndom(m0))=m endif enddo t=t+ncon(m0,ndom(m0))*m if (ncon(m0,ndom(m0)) > 0) ndom(m0)=ndom(m0)+1 enddo ndom(m0)=ndom(m0)-1 enddo if(proc0) then do m0=1,md write(*,*) 'total # of domains: ',ndom(m0),' for m=',m0 enddo write(*,102) t,(2*jmax+1)*(kmax+1) endif 102 format('Total modes accounted for: ',i4,'=(2*jmax+1)*(kmax+1)=',i4) do n=1,nd do m=1,md j=nrr(n) k=mrr_all(m,n) mcon(m,n)=0 jshift=jshift0*k if (lin == 1 .and. linlay) jshift=jshift0 if (k /= 0) then icon(m,n)=(jmax-j)/jshift+(j+jmax)/jshift+1 conpos(m,n)=1+(jmax-j)/jshift else icon(m,n)=1 conpos(m,n)=1 endif enddo enddo do m0=1,md do i=1,ndom(m0) count=1 k = m0-1 do j=-jmax,jmax if (l(j,k) == scon(m0,i)) then do jj=1,l(j,k) jshift=jshift0*k if (lin == 1 .and. linlay) jshift=jshift0 n=1+j+(jj-1)*jshift if (n < 1) n=n+nd mcon(k+1,n)=count enddo count=count+1 endif enddo enddo enddo skip_code=.true. ! to skip some code that needs upgrading some day. if(skip_code) goto 7 if(proc0) write(*,*) 'conpos(m,n)' if(proc0) write(*,100) 0,0,1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16 ! 100 format(i5,17(i3)) do j=-jmax,jmax n=1+j if (n < 1) n=n+nd if(proc0) write(*,100) j,(conpos(m,n),m=1,md) enddo if(proc0) write(*,*) 'icon(m,n)' if(proc0) write(*,100) 0,0,1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16 ! 100 format(i5,17(i3)) do j=-jmax,jmax n=1+j if (n < 1) n=n+nd if(proc0) write(*,100) j,(icon(m,n),m=1,md) enddo if(proc0) write(*,*) 'mcon(m,n)' if(proc0) write(*,100) 0,0,1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16 ! 100 format(i5,17(i3)) do j=-jmax,jmax n=1+j if (n < 1) n=n+nd if(proc0) write(*,100) j,(mcon(m,n),m=1,md) enddo if(proc0) then do m=1,md write(*,*) 'ndom(m)=',ndom(m) do j=1,ndom(m) write(*,*) m,j,' ncon(m,j)=',ncon(m,j),' scon(m,j)=',scon(m,j) enddo enddo endif 7 continue do m0=1,md do j=1,ndom(m0) do i=1,ncon(m0,j) do k=1,scon(m0,j) m=m0 do n=1,nd if ((i == mcon(m,n)).and.(k == conpos(m,n)).and. & (scon(m0,j) == icon(m,n))) then mmcon(m0,j,k+(i-1)*scon(m0,j))=m nncon(m0,j,k+(i-1)*scon(m0,j))=n endif enddo enddo enddo enddo enddo ! do j=1,ndom ! do k=1,ncon(j)*scon(j) ! write(*,103) j,k,mmcon(j,k),nncon(j,k) !103 format(i4,i4,' mmcon: ',i4,' nncon: ',i4) ! enddo ! enddo deallocate (l) end subroutine connect subroutine init_theta_fft use itg_data, only: iperiod use gryffin_grid, only: md, nd ! local variables: integer m, n ! ! Prepare for theta FFT's ! if (iperiod == 3) then call connect else allocate (mcon(md, nd), ncon(md, 1), scon(md, 1)) allocate (conpos(md, nd), icon(md, nd), ndom(md)) ndomz=1 nconz=1 ndom=1 ncon=md*nd scon=1 conpos=1 icon=1 do m=1,md do n=1,nd mcon(m,n)=n+(m-1)*nd enddo enddo endif end subroutine init_theta_fft end module bounds