module nlps ! ! (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. ! ! "NonLinear PseudoSpectral" BD ! Given the fields a, b, ka, and kb this subroutine calculates the ! Poisson bracket ! {a,b} = da/dx db/dy - da/dy db/dx. ! ! The fields arrive in (x,ky,kz) space and must be transformed ! with the appropriate derivatives into (x,y,z) space. ! use gryffin_redist, only: gryffin_redist_type implicit none private public :: nlpsc, init_nlps, maxv, reality integer :: idum real, dimension(:,:), allocatable, save :: kxa2, kya2, kxb2, kyb2 real, dimension(:,:), allocatable :: km real, dimension(:), allocatable :: workx real, dimension(:), allocatable :: worky real, dimension(:), allocatable :: tablex ! mcfft trig array real, dimension(:), allocatable :: tabley ! mcfft trig array integer :: ntablex, ntabley, nworkx, nworky integer, parameter :: nredist = 2 type (gryffin_redist_type), dimension (nredist), save :: redist logical, dimension (nredist) :: redist_initialized = .false. contains subroutine nlpsc(isave, a, b, nl_term, lmax) use itg_data, only: shr use gryffin_layouts use constants, only: ii use gryffin_grid, only: mrr, y0, r0 integer :: iredist integer l, m, n, isave, mr, lmax, idx, idx_low, idx_high integer, save :: lmax_last integer :: nfirst = 1 complex, dimension(:,m_low:,n_low:) :: a, b, nl_term complex, dimension(:,:,:), allocatable :: kxa, kya, kxb, kyb call get_iredist (lmax, lmax_last, iredist) ! Don't recalculate FFT's unnecessarily! ! ! isave determines which fields are new for this call: ! isave=0 -> both fields are new ! isave=1 -> "a" fields are the same as last time nlps was called ! isave=2 -> "b" fields are the same as last time nlps was called ! if(nfirst == 1) then nfirst = 0 call mod_alloc(redist(iredist)) else if(lmax /= lmax_last) then deallocate (kxa2, kya2, kxb2, kyb2, km) call mod_alloc(redist(iredist)) endif endif lmax_last = lmax ! Calculate the d/dy and d/dx terms: if(isave /= 1) then allocate (kxa(lmax, m_low:m_alloc, n_low:n_alloc)) allocate (kya(lmax, m_low:m_alloc, n_low:n_alloc)) do n = n_low, n_high do m = m_low,m_high mr=mrr(m,n) do l = 1, lmax kya(l,m,n) = ii*mr/y0*a(l,m,n) kxa(l,m,n) = ii*max(mr,1)/y0*shr*r0(m,n)*a(l,m,n) enddo enddo enddo endif if(isave /= 2) then allocate (kxb(lmax, m_low:m_alloc, n_low:n_alloc)) allocate (kyb(lmax, m_low:m_alloc, n_low:n_alloc)) do n = n_low, n_high do m = m_low,m_high mr=mrr(m,n) do l = 1, lmax kyb(l,m,n) = ii*mr/y0*b(l,m,n) kxb(l,m,n) = ii*max(mr,1)/y0*shr*r0(m,n)*b(l,m,n) enddo enddo enddo endif ! Transform the vectors i*kx*a, i*kx*b, i*ky*b and i*ky*a into x space: if(isave /= 1) then call real_space(kxa, kxa2, redist(iredist)) deallocate (kxa) call real_space(kya, kya2, redist(iredist)) deallocate (kya) endif if(isave /= 2) then call real_space(kxb, kxb2, redist(iredist)) deallocate (kxb) call real_space(kyb, kyb2, redist(iredist)) deallocate (kyb) endif ! Perform the multiplications: idx_low = redist(iredist)%y_lo%llim_proc idx_high = redist(iredist)%y_lo%ulim_proc do idx = idx_low, idx_high km(:,idx) = kxa2(:,idx) * kyb2(:,idx) - kxb2(:,idx) * kya2(:,idx) enddo ! Transform back. call k_space(nl_term, km, redist(iredist)) end subroutine nlpsc subroutine real_space(ak, a, r) ! ! Transform a variable from k space into real space and dealias. ! use mp, only: iproc use gryffin_grid, only: nddp, nadd, nalias, md, malias, & r0, z0, mrr, nrr, ldb, nd, nmax use gryffin_layouts use redistribute, only: gather use fft_work use itg_data, only: ncycle, nprnt type (gryffin_redist_type), intent (in out) :: r complex, dimension(:,m_low:,n_low:) :: ak complex, dimension(ldb, m_low:m_alloc, n_low:n_alloc) :: aksh ! complex, dimension(size(ak,1),size(ak,2),size(ak,3)) :: aksh real, dimension(:,r%y_lo%llim_proc:) :: a complex, dimension(nalias,r%x_lo%llim_proc:r%x_lo%ulim_alloc) :: ka complex, dimension(malias/2+1,r%y_lo%llim_proc:r%y_lo%ulim_alloc) :: aky integer :: l, m, n, index, imode, isign, idx, np, dn real :: scale ! ! Protect against uninitialized vars on procs with nothing to do. ! if (r%x_lo%llim_proc > r%x_lo%ulim_proc) ka = 0. if (r%y_lo%llim_proc > r%y_lo%ulim_proc) aky = 0. ! ! ExB shift, serial only ! if (n_low.ne.1) write (*,*) n_low,n_high,'parallelization problem ExB' if (n_high.ne.nd) write (*,*) n_low,n_high,'parallelization problem ExB' do n=n_low,n_high do m=m_low,m_high if (mrr(m,n) /= 0) then dn=r0(m,n)*float(nmax)*float(mrr(m,n))/z0-float(nrr(n))+0.5 if (r0(m,n) < z0*float(nrr(n))/float(nmax)/float(mrr(m,n))- & 1.e-10) dn=dn-1 else dn=0 endif np=n+dn if (np > nd) np=np-nd if (np < 1) np=np+nd ! if ((n /= np) .and. (mod(ncycle,nprnt) == 0)) & ! write(*,*) 'n,np: ',n,np, ' (l,m,n): ',m,n,r0(m,n), & ! z0*float(nrr(n))/float(nmax)/float(mrr(m,n)) ! if (n /= np) then ! write(*,*) 'n,np: ',n,np, ' (l,m,n): ',m,n,r0(m,n), & ! z0*float(nrr(n))/float(nmax)/float(mrr(m,n)) ! endif do l=1,ldb aksh(l,m,np)=ak(l,m,n) enddo enddo enddo ! ! Redistribute: ! call gather(r%kx2x, aksh, ka) ! ! de-alias ! do n=nddp+1,nddp+nadd ka(n,:)=0.0 enddo ! ! Transform: ! scale=1. isign=1 do idx = r%x_lo%llim_proc, r%x_lo%ulim_proc call ccfft(isign, nalias, scale, ka(1:nalias,idx), & ka(1:nalias,idx), tablex, workx, 0) enddo ! ! Redistribute: ! call gather (r%x2y, ka, aky) ! ! Normalize the ky /= components: ! aky(2:md,:) = aky(2:md,:)/2.0 ! ! dealias for high ky: ! aky(md+1:malias/2+1,:)=0.0 scale = 1. isign = 1 do imode=r%y_lo%llim_proc, r%y_lo%ulim_proc call csfft(isign, malias, scale, aky(1:malias/2+1,imode), & a(1:malias,imode), tabley, worky, 0) enddo end subroutine real_space subroutine k_space(ak, a, r) use gryffin_grid, only: malias, md, nalias, nddp, nadd, & r0, z0, mrr, nrr, nmax, ldb, nd use mp, only: iproc use gryffin_layouts use redistribute, only: scatter use fft_work use itg_data, only: ncycle, nprnt type (gryffin_redist_type), intent (in out) :: r real, dimension(:,r%y_lo%llim_proc:) :: a complex, dimension(:,m_low:,n_low:) :: ak complex, dimension(ldb, m_low:m_alloc, n_low:n_alloc) :: aksh ! complex, dimension(size(ak,1),size(ak,2),size(ak,3)) :: aksh complex, dimension(nalias, r%x_lo%llim_proc:r%x_lo%ulim_alloc) :: ka complex, dimension(malias/2+1,r%y_lo%llim_proc:r%y_lo%ulim_alloc) :: aky real scale integer isign integer m,idx,l,n, dn, np if (r%x_lo%llim_proc > r%x_lo%ulim_proc) ka = 0. if (r%y_lo%llim_proc > r%y_lo%ulim_proc) aky = 0. scale = 1./float(malias) isign = -1 do m=r%y_lo%llim_proc, r%y_lo%ulim_proc call scfft(isign, malias, scale, a(1:malias,m), & aky(1:malias/2+1,m), tabley, worky, 0) enddo ! ! get ky>0 components (factor of 2 for sin/cos vs. complex conventions). ! aky(2:md,:) = aky(2:md,:)*2.0 ! ! Redistribute ! call scatter (r%x2y, aky, ka) scale=1./float(nalias) isign=-1 do idx = r%x_lo%llim_proc, r%x_lo%ulim_proc call ccfft(isign, nalias, scale, ka(1:nalias,idx), & ka(1:nalias,idx), tablex, workx, 0) enddo ! ! Enforce reality for ky = 0 modes ! call reality(r, ka) ! ! Redistribute ! call scatter (r%kx2x, ka, aksh) ! ! ExB shift, serial only ! do n=n_low,n_high do m=m_low,m_high if (mrr(m,n) /= 0) then dn=r0(m,n)*float(nmax)*float(mrr(m,n))/z0-float(nrr(n))+0.5 if (r0(m,n) < z0*float(nrr(n))/float(nmax)/float(mrr(m,n))- & 1.e-10) dn=dn-1 else dn=0 endif np=n+dn if (np > nd) np=np-nd if (np < 1) np=np+nd ! if ((n /= np) .and. (mod(ncycle,nprnt) == 0)) & ! write(*,*) 'n,np: ',n,np, ' (l,m,n): ',m,n,r0(m,n), & ! z0*float(nrr(n))/float(nmax)/float(mrr(m,n)) ! if (n /= np) then ! write(*,*) 'n,np: ',n,np, ' (l,m,n): ',m,n,r0(m,n), & ! z0*float(nrr(n))/float(nmax)/float(mrr(m,n)) ! endif do l=1,ldb ak(l,m,n)=aksh(l,m,np) enddo enddo enddo end subroutine k_space subroutine init_nlps use gryffin_grid, only: malias, nalias use fft_work complex, dimension(malias/2+1) :: zdum real, dimension(malias) :: xdum idum=0 ntabley = csfft_table0 + csfft_table1*malias nworky = csfft_work0 + csfft_work1*malias allocate (worky(nworky), tabley(ntabley)) call csfft(0,malias,1.,zdum,xdum,tabley,worky,0) ntablex=ccfft_table0 + ccfft_table1*nalias nworkx=ccfft_work0 + ccfft_work1*nalias allocate (tablex(ntablex), workx(nworkx)) call ccfft(0, nalias, 1.0, zdum, zdum, tablex, workx, 0) end subroutine init_nlps subroutine mod_alloc(r) use gryffin_grid, only: malias type (gryffin_redist_type), intent (in) :: r allocate (kxa2(malias, r%y_lo%llim_proc:r%y_lo%ulim_alloc), & kya2(malias, r%y_lo%llim_proc:r%y_lo%ulim_alloc), & kxb2(malias, r%y_lo%llim_proc:r%y_lo%ulim_alloc), & kyb2(malias, r%y_lo%llim_proc:r%y_lo%ulim_alloc), & km(malias, r%y_lo%llim_proc:r%y_lo%ulim_alloc)) end subroutine mod_alloc subroutine get_iredist (lmax, lmax_old, iredist) use gryffin_redist, only: init_gryffin_redist, delete_gryffin_redist use gryffin_grid, only: md, nd, malias, nalias, nddp, nadd implicit none integer, intent (in) :: lmax, lmax_old integer, intent (out) :: iredist integer :: i do i = 1, nredist if (redist_initialized(i) .and. lmax == redist(i)%y_lo%lmax) then iredist = i return end if end do do i = 1, nredist if (.not. redist_initialized(i)) then call init_gryffin_redist (lmax, md, nd, malias, nalias, nddp, nadd, redist(i)) iredist = i redist_initialized(i) = .true. return end if end do do i = 1, nredist if (lmax_old /= redist(i)%y_lo%lmax) then call delete_gryffin_redist (redist(i)) call init_gryffin_redist (lmax, md, nd, malias, nalias, nddp, nadd, redist(i)) iredist = i return end if end do call delete_gryffin_redist (redist(1)) call init_gryffin_redist (lmax, md, nd, malias, nalias, nddp, nadd, redist(1)) iredist = 1 end subroutine get_iredist subroutine maxv(b, max_vel) use itg_data, only: shr, cflx, cfly use gryffin_layouts use mp, only: max_allreduce use gryffin_grid, only: ldb, mrr, malias, y0, r0 use constants, only: ii complex, dimension(:,m_low:,n_low:) :: b complex, dimension(:,:,:), pointer :: kxb, kyb real max_vel integer l, m, n, imode, mr, iredist call get_iredist (ldb, 0, iredist) allocate (kxb(ldb, m_low:m_alloc, n_low:n_alloc)) allocate (kyb(ldb, m_low:m_alloc, n_low:n_alloc)) do n = n_low, n_high do m = m_low, m_high mr=mrr(m,n) do l = 1, ldb kyb(l,m,n)=ii*mr/y0*b(l,m,n) kxb(l,m,n)=ii*max(mr,1)/y0*shr*r0(m,n)*b(l,m,n) enddo enddo enddo deallocate (kxb2, kyb2) allocate (kxb2(malias, redist(iredist)%y_lo%llim_proc:redist(iredist)%y_lo%ulim_alloc), & kyb2(malias, redist(iredist)%y_lo%llim_proc:redist(iredist)%y_lo%ulim_alloc)) call real_space(kxb, kxb2, redist(iredist)) deallocate (kxb) call real_space(kyb, kyb2, redist(iredist)) deallocate (kyb) do imode=redist(iredist)%y_lo%llim_proc, redist(iredist)%y_lo%ulim_proc do m=1,malias max_vel=max(max(max_vel,abs(kyb2(m,imode))*cflx), & abs(kxb2(m,imode))*cfly) enddo enddo call max_allreduce(max_vel) end subroutine maxv subroutine reality (r, ka) use gryffin_grid, only: nddm, nddp, nadd use gryffin_redist, only: gryffin_redist_type, ky_is_zero type (gryffin_redist_type), intent (in) :: r complex, dimension(:,r%x_lo%llim_proc:) :: ka integer :: n, idx do idx = r%x_lo%llim_proc, r%x_lo%ulim_proc if (ky_is_zero(r, idx)) then do n = 1, nddm ka(nddp+nadd+n,idx) = conjg(ka(nddp+1-n,idx)) enddo endif enddo end subroutine reality end module nlps