module field_eq implicit none contains subroutine poisson(phi, den_e, den, tp) ! ! (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. ! use constants, only: ii use itg_data use mp, only: proc0 use gryffin_layouts use bounds, only: bcon use linear use gryffin_grid, only: ld, l_left, l_right, ldb, nd, x0, rkx use fields, only: e_denk, n_e, phiav_old, phi_bk ! define the arguments: complex, dimension(:,m_low:,n_low:,:) :: den, tp complex, dimension(:,m_low:,n_low:) :: phi, den_e ! declare local workspace complex, dimension( m_low:m_alloc, n_low:n_alloc) :: phiav complex, dimension(ld, m_low:m_alloc, n_low:n_alloc) :: nbar complex tmp,wgt_avg real, dimension(l_left:l_right) :: dz real znorm ! declare other local variables: integer i,l,m,n ! ! start executable code: ! ! iphi00=9 = set phi=0 for all Fourier components, not just (m=0,n=0). ! (no electric field, for neoclassical comparisons): if(iphi00 == 9) return ! ! load stuff for field line average ! do l=l_left,l_right ! BD: backwards compatibility issue: if((igeo == 0) .and. (iphi00 == 0)) then dz(l)=2.*x0/ldb else dz(l)=2.*x0/ldb*jacobian(l) endif enddo ! ! Calculate the real-space ion density (summed over species) ! nbar = 0. phiav = 0. do i=1,nspecies do l=1,ld-1 nbar(l,:,:)=nbar(l,:,:)+(den(l,:,:,i)*nwgt(l,:,:,i) & +tp(l,:,:,i)*twgt(l,:,:,i))*g0wgtp(l,:,:,i) enddo enddo ! GWH: ! Parallel filtering to model finite-size particles. The perp filter ! was included in nwgt and twgt in flrinit.f. Formally, the parallel ! filtering should go in n_i_bar (as the particle charges are assigned ! to the grid), and not just as a filtering on phi. In general geometry, ! with the calculation, there may be some subtle differences... ! However, I'm not sure which order the particle people actually do it ! in. (If needed, could change this to be a filtering of phi at the end ! of this subroutine...) ! ! if(s_par /= 0.0) then ! call bcon(nbar) ! apply boundary conditions ! call kparfft(nbar, 13, 0.) ! endif ! ! If iphi00=2, require field line average of phi = 0 to ensure no ! radial electron transport: ! ! includes variation of k_x^2+k_y^2 along field line ! ! ! Substantial simplification of this calculation can be had for ! nonlinear runs, where we always know ahead of time that ! rky2(m,n)=0 iff m=1 and ! rkx2(1,1,n) /= 0 for n>1. ! (This is true as long as shr /= 0.) ! However, I am leaving the old logic for linear runs, in which the ! layout in (m,n) is not known ahead of time. Dorland. ! !fpp$ skip R ! ! Calculate for electrostatic runs, or for ! em runs with trapped electrons: if (beta_e <= 0. .or. (beta_e > 0. .and. epse > 0.)) then if (iphi00 == 2 .and. (lin == 1 .or. shr == 0.)) then if(lin == 0 .and. shr == 0.) then ! The definition of phiav is not clear in the code in this case. write(*,*) "Error in poisson?" write(*,*) "should be checked" stop endif ! Since we are dropping the possibility of doing runs with a band of ! k_theta's, this simplifies to: do n = n_low, n_high m = 1 if(idx_local(m_lo, m, n)) then wgt_avg=0. if (n /= 1) then tmp=0. do l=l_left,l_right-1 tmp=tmp+(nbar(l,m,n)-e_denk(l,m,n))*pfilter(l,m,n)*jacobian(l) enddo wgt_avg=tmp/denom0(m,n)*(1.-stochf) endif phiav(m,n)=wgt_avg endif enddo endif ! ! Nonlinear replacement loop: ! m = 1 if (iphi00 == 2 .and. lin == 0) then do n = n_low, n_high if (n == 1) cycle if(idx_local(m_lo, m, n)) then do l=l_left,l_right-1 phiav(m,n)=phiav(m,n)+(nbar(l,m,n)-e_denk(l,m,n))*pfilter(l,m,n) & /denom0(m,n)*(1.-stochf)*jacobian(l) ! ! In adiabatic electron code, phiav was (T_i/T_e), now it is ! enddo endif enddo endif endif !fpp$ noskip R ! Electromagnetic version ! (may need generalization for multi-species case) PS 4/4/97 ! if (beta_e > 0.) then do n = n_low, n_high do m = m_low, m_high do l=1,ld-1 if (pfilter2(l,m,n) == 0.) then ! Need to think about physics involoved here... phi(l,m,n)=0. n_e(l,m,n)=nbar(l,m,n) else ! testing trapped electron version ! used to be =(nbar(l,m,n)-den_e(l,m,n))/(pfilter2(l,m,n)) phi(l,m,n)=(nbar(l,m,n) - den_e(l,m,n) & ! -e_denk(l,m,n) & ! +tiovte*(phiav(m,n)+phi_bk(l,m,n)-phiav_old(m,n)) & )/(pfilter2(l,m,n)) endif enddo enddo enddo if (epse > 0.) phiav_old = phiav else ! ! Solve Poisson's equation ! ! do i=1,3 ! for iteration of Poisson's Eq. ! do l=1,ld-1 phi(l,:,:)=(nbar(l,:,:)-e_denk(l,:,:)+tiovte*( & +phiav(:,:)+phi_bk(l,:,:)-phiav_old(:,:)))*pfilter(l,:,:) ! don't use phiav, just phi_bk from previous timestep (doesn't work) ! phi(l,:,:)=(nbar(l,:,:)-e_denk(l,:,:)+tiovte*( ! . +phi_bk(l,:,:)))*pfilter(l,:,:) ! enddo endif ! Add external phi for Rosenbluth-Hinton test. ! (GWH 4/25/97): m = 1 do n = n_low, n_high if(idx_local(m_lo, m, n)) then if(pmag == 0.0) phi(:,m,n) = phi(:,m,n) -ii/rkx(:,m,n) endif enddo ! call baphi(phi) ! call kapav(phi_ba,den) if (epse > 0.0) then phiav_old = phiav else phiav_old = 0.0 endif ! ! GWH ?? 9/19/94: The use of phi_bk (and phiav_old) from the previous ! time step appears to be a kind of 1-iteration approximation. We ! should eventually test whether additional iterations on each call ! to poisson.f are required to get a converged solution to the ! integral equations implied by the bounce and kappa averages... ! ! I am confused about why the Poisson Eq. is so complicated now. ! ! Also, more thought needs to go into a general geometry formulation. ! This works for now, where pfilter(l,1,n) is independent of l in ! the large-aspect ratio limit... ! ! BD 6/22/95: ! Should be right for general geometry now. pfilter doesn't pick up any ! new theta dependence. ! ! ! set Phi(ky=0, kz=0)=0 if iphi00=0 m=1, n /= 1 modes ! m = 1 if (iphi00 == 0 .and. nd /= 1) then do n = n_low, n_high if (idx_local(m_lo, m, n)) then if (n == 1) cycle wgt_avg=0.0 znorm=0. do l=l_left,l_right-1 znorm=znorm+dz(l) wgt_avg=wgt_avg+phi(l,m,n)*dz(l) enddo do l=l_left,l_right-1 phi(l,m,n)=phi(l,m,n)-wgt_avg/znorm enddo endif enddo endif ! ! Do boundary conditions ! if(s_par /= 0.0) then if(proc0) write(*,*) 'Error in poisson related to s_par. Check it out!' ! BD ! This needs to be fixed. It was a bug, but probably didn't show up for nspecies = 1 ! call kparfft(phi, 13, 0.) endif call bcon(phi) end subroutine poisson end module field_eq