module bounce_avg ! ! (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 real, dimension(:), allocatable :: dkap ! dkappa for bounce averaging real, dimension(:), allocatable :: taub ! bounce time, trapped and passing ! real, dimension(:,:), allocatable :: pmat ! matrix for poisson inversion ! ! For electron collisions: ! ! bratio: bounce avg. of (B-Bmin)/B ! jac: electron Jacobian ! diff: parameter for bounce averaging operator real, dimension(:), allocatable :: bratio, jac, diff real, dimension(:,:), allocatable :: wgtba1 ! weight for bounce avg. real, dimension(:,:), allocatable :: wgtba2 ! weight for kappa avg. contains subroutine kapav(phi, den) ! ! Average n_e and bounce averaged phi over pitch angle kappa ! i.e. given phi_ba(kappa) and den(kappa), get ! phi_bk(theta) and denk(theta). ! use gryffin_layouts use gryffin_grid, only: ldb, kd, kdpass use fields, only: e_denk, phi_bk complex, dimension(:,m_low:,n_low:) :: phi, den integer l,j,m,n,l0 do m = m_low, m_high if(m == 1) cycle do n = n_low, n_high do j=1,ldb l0=1+abs(j-ldb/2-1) e_denk(j,m,n)=0. phi_bk(j,m,n)=0. do l=l0,kd e_denk(j,m,n)=e_denk(j,m,n)+dkap(l)*den(l,m,n)*wgtba2(l,j) phi_bk(j,m,n)=phi_bk(j,m,n)+dkap(l)*phi(l,m,n)*wgtba2(l,j) enddo enddo enddo enddo m=1 do n = n_low, n_high if(idx_local(m_lo, m, n)) then do j=1,ldb l0=1+abs(j-ldb/2-1) e_denk(j,m,n)=0. phi_bk(j,m,n)=0. do l=l0,kd-1 e_denk(j,m,n)=e_denk(j,m,n)+dkap(l)*(den(l,m,n)+den(l+1,m,n))/2*wgtba2(l,j) phi_bk(j,m,n)=phi_bk(j,m,n)+dkap(l)*(phi(l,m,n)+phi(l+1,m,n))/2*wgtba2(l,j) enddo do l=max(kd,l0),kd+1 e_denk(j,m,n)=e_denk(j,m,n)+dkap(l)*den(l,m,n)*wgtba2(l,j) phi_bk(j,m,n)=phi_bk(j,m,n)+dkap(l)*phi(l,m,n)*wgtba2(l,j) enddo do l=kd+2,kdpass e_denk(j,m,n)=e_denk(j,m,n)+dkap(l)*(den(l,m,n)+den(l-1,m,n))/2*wgtba2(l,j) phi_bk(j,m,n)=phi_bk(j,m,n)+dkap(l)*(phi(l,m,n)+phi(l-1,m,n))/2*wgtba2(l,j) enddo enddo endif enddo end subroutine kapav subroutine baphi(phi) ! ! (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. ! ! ! Get bounce averaged potential, ! i.e. load phi_ba(kappa) from phi(theta) ! use gryffin_layouts use gryffin_grid, only: kd, ldb, kdpass use fields, only: phi_ba complex, dimension(:,m_low:,n_low:) :: phi integer :: j, l, m, n, ntheta0, ntheta1 ! for k_theta>0 modes, we only need to calculate the bounce-averaged ! phi for trapped particles: phi_ba=0. do n = n_low, n_high do m = m_low, m_high do l=1,kd ntheta0=ldb/2+l ntheta1=ldb/2+2-l do j=ntheta1,ntheta0 phi_ba(l,m,n)=phi_ba(l,m,n)+phi(j,m,n)*wgtba1(l,j) enddo phi_ba(l,m,n)=phi_ba(l,m,n)/taub(l) enddo enddo enddo ! for k_theta=0 modes, we need the bounce-averaged phi for passing ! particles also: m=1 do n = n_low, n_high if(idx_local(m_lo, 1, n)) then do l=kd+1,kdpass do j=1,ldb phi_ba(l,m,n)=phi_ba(l,m,n)+phi(j,m,n)*wgtba1(l,j) enddo phi_ba(l,m,n)=phi_ba(l,m,n)/taub(l) enddo endif enddo end subroutine baphi subroutine init_bounce_avg use constants, only: pi use itg_data, only: igeo, omegade, omegade0, omegade1, gpar, qsf use itg_data, only: b_mag, omega_gb, omega_kap use itg_data, only: omega_gb0, omega_kap0, omega_gb1, omega_kap1 use itg_data, only: rho, vt, epsn, epse, shr, alpha use gryffin_grid, only: kdpass, ld, ldb, kd, r, kap, md, & mrr, x0, meq, neq, rky use gryffin_layouts use mp, only: proc0 real, dimension(:), allocatable :: ellk real :: s21bcE, s21bbE, elliptic, kshift real :: zeta1, zeta2, dzeta real :: theta, dtheta, dx, arg, epsb, x real :: a1, b1, c1, d1, kapav real :: kapp, kap2p, taubp, bratiop real :: bmin, b_max, zkapmax integer :: ntheta0, ntheta1, ifail integer :: l, m, n, i, j, l0 external s21bcE, s21bbE ! Local Variables: integer, parameter :: nw = 500 ! # of steps in bounce-averages (must be even) call alloc_bounce_avg allocate (ellk(kdpass)) allocate (omegade(kdpass, m_low:m_lo%mu_alloc, & n_low:m_lo%nu_alloc)) allocate (omegade0(kdpass, m_low:m_lo%mu_alloc, & n_low:m_lo%nu_alloc)) allocate (omegade1(kdpass, m_low:m_lo%mu_alloc, & n_low:m_lo%nu_alloc)) ellk = 0.; omegade = 0.; omegade0 = 0.; omegade1 = 0. ! ! kap(l) = kappa pitch-angle grid for bounce-averaged electrons ! if (igeo == 0) then ! ! circular concentric limit ! do l=1,kd j=ldb/2+l kap(l)=sin(r(j)/2.) kap(l)=max(0.0,kap(l)) enddo ! kap(1)=kap(2)/4 ! put the first grid point slightly off kap=0 do l=1,kd ellk(l)=s21bbE(0.,1.-kap(l)**2,1.,ifail) enddo epsb=epse/(1.+epse) else ! ! general geometry ! bmin=b_mag(1) b_max=b_mag(1) do l=1,ld if (b_mag(l) > b_max) b_max=b_mag(l) if (b_mag(l) < bmin) bmin=b_mag(l) enddo if(proc0) write(6,*) 'bmin,bmax: ',bmin,b_max epsb=(b_max-bmin)/2/b_max do l=1,kd j=ldb/2+l kap(l)=sqrt((1-bmin/b_mag(j))/2/epsb) kap(l)=max(0.0,kap(l)) enddo endif ! BD: 1.5.99 ! This if statement is redundant; has been this way for a long time. ! if(epse > 0.0) then zkapmax=1./sqrt(2.*epsb) else zkapmax=2.0 endif ! eventually would like the kap grid to be more closely space ! around the kap=1 point for kap>1 as well as for kap<1. For now ! just use a uniform grid, but at least put the first point at ! kap>1 equally spaced from the boundary as the kap<1 point, so ! that kap=1 is exactly between two grid points: do l=kd+1,kdpass kap(l)=1.0+(1.0-kap(kd))+(l-kd-1)*(zkapmax-1.0)/(kdpass-kd) enddo kshift=(2-kap(kd)-zkapmax*(1-kap(kd)))/kap(kd) do l=kd+1,kdpass kap(l)=(1-kap(kd-l+1+kd))*(zkapmax-kshift)+kshift ellk(l)=s21bbE(0.,1.-1./kap(l)**2,1.,ifail) enddo ! define a mirror point on the other side of the maximum kappa ! to make differencing formulas easier: kap(kdpass+1)=zkapmax+(zkapmax-kap(kdpass)) ! ! bounce averaged precession frequency: ! !fpp$ skip R do n = n_low, n_high do m = m_low, m_high do l=1,kd if (kap(l) /= 1.0) then elliptic=1.-kap(l)**2/3.*s21bcE(0.,1.-kap(l)**2,1.,ifail) & /s21bbE(0.,1.-kap(l)**2,1.,ifail) else elliptic=0. endif omegade(l,m,n)=2.*epsn*rky(m,n)*(elliptic-0.5 & +2.*shr*(elliptic+kap(l)**2-1.) & -alpha*2./3.*(1-kap(l)**2-(1-2*kap(l)**2)*elliptic)) ! ! special case for k_theta=0, k_r=const. ! if ((m == meq).or.((md == 1).and.mrr(m,n) == 0)) omegade(l,m,n)=0. ! ! local if qsf<0 ! ! If igeo=1, qsf can be negative, so in the new section, ! don't keep this option: ! if (qsf < 0 .and. igeo == 0) omegade(l,m,n)=2.*epsn*rky(m,n) enddo enddo enddo !fpp$ noskip R ! Define dkap as the width of the zone about each grid point: ! dkap(1)=kap(2)/2 ! do l=2,kdpass ! dkap(l)=(kap(l+1)-kap(l-1))/2 ! enddo do l=1,kd-1 dkap(l)=kap(l+1)-kap(l) enddo dkap(kd)=1.-kap(kd) dkap(kd+1)=kap(kd+1)-1. do l=kd+2,kdpass dkap(l)=kap(l)-kap(l-1) enddo ! ! weights for bounce averages ! ! The algorithm for the integration below could probably be more ! effficient and converge faster (so that nw doesn't need to be ! so large), perhaps using unequally spaced points (like sccf.f ! in my FPP code) or using a canned routine. ! But it works for now... ! allocate (wgtba1(kdpass, ld), wgtba2(kdpass, ld)) if (igeo == 0) then ! ! circular concentric limit ! wgtba1 = 0. dtheta=2.*x0/float(ldb) ! ! trapped ! do l=2,kd ntheta0=ldb/2+l ntheta1=ldb/2+2-l j=ntheta1 ! use variable transformation: sin(theta/2)=kappa sin(x) if (sin(r(j)/2)/kap(l) < -1.0) then zeta1=-pi/2 else zeta1=asin(sin(r(j)/2)/kap(l)) endif zeta2=asin(sin((r(j)+dtheta/2)/2)/kap(l)) ! write(*,101) kap(l),r(j),r(j)+dtheta/2, & ! sin(r(j)/2)/kap(l),sin((r(j)+dtheta/2)/2)/kap(l),zeta1,zeta2 101 format(7f15.9) dzeta=zeta2-zeta1 do i=1,nw x=(float(i)-0.5)/float(nw)*dzeta+zeta1 dx=1./float(nw)*dzeta arg=1-kap(l)**2*sin(x)**2 wgtba1(l,j)=wgtba1(l,j)+2.*dx*sqrt((1-epse+2*epse*arg)/arg/(1+epse)) enddo do j=ntheta1+1,ntheta0-1 zeta1=asin(sin((r(j)-dtheta/2)/2)/kap(l)) zeta2=asin(sin((r(j)+dtheta/2)/2)/kap(l)) ! write(*,101) kap(l),r(j)-dtheta/2,r(j)+dtheta/2, & ! sin((r(j)-dtheta/2)/2)/kap(l),sin((r(j)+dtheta/2)/2)/kap(l),zeta1,zeta2 dzeta=zeta2-zeta1 do i=1,nw x=(float(i-1)+0.5)/float(nw)*dzeta+zeta1 dx=1./float(nw)*dzeta arg=1-kap(l)**2*sin(x)**2 wgtba1(l,j)=wgtba1(l,j)+2.*dx*sqrt((1-epse+2*epse*arg)/arg/(1+epse)) enddo enddo j=ntheta0 zeta1=asin(sin((r(j)-dtheta/2)/2)/kap(l)) if (sin(r(j)/2)/kap(l) > 1.0) then zeta2=pi/2 else zeta2=asin(sin(r(j)/2)/kap(l)) endif dzeta=zeta2-zeta1 ! write(*,101) kap(l),r(j)-dtheta/2,r(j), & ! sin((r(j)-dtheta/2)/2)/kap(l),sin(r(j)/2)/kap(l),zeta1,zeta2 do i=1,nw x=(float(i)-0.5)/float(nw)*dzeta+zeta1 dx=1./float(nw)*dzeta arg=1-kap(l)**2*sin(x)**2 wgtba1(l,j)=wgtba1(l,j)+2.*dx*sqrt((1-epse+2*epse*arg)/arg/(1+epse)) enddo enddo ! ! passing ! do l=kd+1,kdpass do j=1,ldb wgtba1(l,j)=0. do i=1,nw x=(float(i-1)/float(nw)+0.5/float(nw))*dtheta-dtheta/2 dx=1./float(nw) theta=r(j)+x arg=2*epse/(1+epse*cos(theta))*(kap(l)**2-sin(theta/2)**2) wgtba1(l,j)=wgtba1(l,j)+dx/sqrt(arg/2./epsb) enddo wgtba1(l,j)=wgtba1(l,j)*dtheta enddo enddo ! ! ntheta0,ntheta1 are chosen so that: ! For trapped particles, the bounce-average symmetrically involves ! phi points above and below the midplane. ! For passing particles, the bounce-average involves all phi points ! from l=1 to ldb (the ld point is identical to the l=1 point). ! ! kap=0 limiting value: wgtba1(1,ldb/2+1)=2*pi ! limiting value... ! calculate bounce-time numerically (for passing particles also...): do l=1,kd ntheta0=ldb/2+l ntheta1=ldb/2+2-l do j=ntheta1,ntheta0 taub(l)=taub(l)+wgtba1(l,j) enddo enddo do l=kd+1,kdpass taub(l)=0. do j=1,ldb taub(l)=taub(l)+wgtba1(l,j) enddo enddo ! debug: do l=1,kdpass if(idx_local(m_lo, 2, 1)) write(6,102) l,kap(l),taub(l), & 4*ellk(l),omegade(l,2,1)/rky(2,1)/epsn enddo 102 format( 'l,kap(l),taub(l),4*ellk(l) =',i3,4f14.7) ellk = taub/4 ! ! ?? I think I've now generalized all of this subroutine to use ! a grid-centered approach, with bounce-averages symmetrically ! involving phi above and below the midplane... ! wgtba2 = 0. do j=ldb/2+1,ld l0=j-ldb/2 do l=l0,kdpass if (l <= kd) then wgtba2(l,j)=sqrt(2*epse/(1+epse*cos(r(j))))* & (sqrt( (kap(l)+dkap(l))**2 - sin(r(j)/2)**2 ) & -sqrt( kap(l)**2 - sin(r(j)/2)**2 )) else wgtba2(l,j)=sqrt(2*epse/(1+epse*cos(r(j))))* & (sqrt( kap(l)**2 - sin(r(j)/2)**2 ) & -sqrt( (kap(l)-dkap(l))**2 - sin(r(j)/2)**2 )) endif wgtba2(l,j)=wgtba2(l,j)/dkap(l) ! B(theta)=1/(1+eps*cos(theta)) included to get particle conservation. ! symmetry, for theta<0: ! ?? GWH: also need to put a 1/B into the flux-surface averages... wgtba2(l,ldb-j+2)=wgtba2(l,j) enddo enddo ! ?? GWH: The bounce-averages (wgtba1), kappa-averages (wgtba2), ! and the form of the coefficients of the pitch-angle-scattering ! operator, are related in various ways. Some of these relations ! are important for particle conservation, and so it may be possible ! to derive some of these coefficients from each other in a way ! that preserves conservation properties exactly (or at least to ! within roundoff error). For now, just stick to the above methods, ! as they should preserve the conservation properties to second ! order in the grid spacings... ! BD: This is not used, so I commented it out. ! load matrix for poisson inversion ! do j=1,ldb ! l1=1+abs(j-ldb/2-1) ! do l=l1,kdpass ! if ((l == l1).and.(l1 < kd)) then ! wt2(l,j)=dkap(l)*wgtba2(l,j)/2 ! elseif ((l > l1).and.(l < kd)) then ! wt2(l,j)=(dkap(l)*wgtba2(l,j)+dkap(l-1)*wgtba2(l-1,j))/2 ! elseif ((l == kd).and.(l1 /= kd)) then ! wt2(l,j)=dkap(l)*wgtba2(l,j)+dkap(l-1)*wgtba2(l-1,j)/2 ! elseif ((l == kd).and.(l1 == kd)) then ! wt2(l,j)=dkap(l)*wgtba2(l,j) ! elseif (l == kd+1) then ! wt2(l,j)=dkap(l)*wgtba2(l,j)+dkap(l+1)*wgtba2(l+1,j)/2 ! elseif ((l > kd+1).and.(l < kdpass)) then ! wt2(l,j)=(dkap(l)*wgtba2(l,j)+dkap(l+1)*wgtba2(l+1,j))/2 ! elseif (l == kdpass) then ! wt2(l,j)=dkap(l)*wgtba2(l,j)/2 ! endif ! enddo ! do j1=1,ldb ! pmat(j,j1)=0. ! l0=1+max(abs(j-ldb/2-1),abs(j1-ldb/2-1)) ! do l=l0,kdpass ! pmat(j,j1)=pmat(j,j1)+wt2(l,j)*wgtba1(l,j1)/taub(l) ! enddo ! temp(j,j1)=pmat(j,j1) ! enddo ! enddo ! ! invert ! ifail=0 ! call f01aaE(temp,ld,ldb/2+1,pmat,ld,wrkinv,ifail) ! not used else ! ! general geometry bounce average weights ! wgtba1 = 0. dtheta=2.*x0/float(ldb) ! ! trapped (assume B and J vary linearly over theta grid pt. interval) ! ! special case for l=2 ! assume up-down symmetric for now... ! ! c1<0 ! l=2 ntheta0=ldb/2+l ntheta1=ldb/2+2-l ! gpar is not a function of theta any more... ! b1=1/gpar(ldb/2+1) ! a1=(1/gpar(ldb/2)-b1)/dtheta**2 b1=1/gpar a1=0. d1=b_mag(ldb/2+1) c1=(b_mag(ldb/2)-d1)/dtheta**2 c1=-c1/bmin*(1-2*epsb*kap(l)**2) d1=1-d1/bmin*(1-2*epsb*kap(l)**2) j=ntheta1 wgtba1(l,j)=b1/sqrt(-c1)*(asin(sqrt(-c1/d1)*(r(j)+dtheta/2))+pi/2) j=ntheta1+1 wgtba1(l,j)=b1/sqrt(-c1)*asin(sqrt(-c1/d1)*(r(j)+dtheta/2)) & -b1/sqrt(-c1)*asin(sqrt(-c1/d1)*(r(j)-dtheta/2)) j=ntheta0 wgtba1(l,j)=b1/sqrt(-c1)*(pi/2-asin(sqrt(-c1/d1)*(r(j)-dtheta/2))) do l=3,kd ntheta0=ldb/2+l ntheta1=ldb/2+2-l j=ntheta1 a1=0. b1=1/gpar c1=(b_mag(j+1)-b_mag(j))/dtheta d1=b_mag(j)-r(j)/dtheta*(b_mag(j+1)-b_mag(j)) c1=-c1/bmin*(1-2*epsb*kap(l)**2) d1=1-d1/bmin*(1-2*epsb*kap(l)**2) wgtba1(l,j)=2/c1*b1*sqrt(c1*(r(j)+dtheta/2)+d1) !& ! -2/c1*b1*sqrt(c1*r(j)+d1) do j=ntheta1+1,ldb/2-1 a1=0. b1=1/gpar c1=(b_mag(j)-b_mag(j-1))/dtheta d1=b_mag(j-1)-r(j-1)/dtheta*(b_mag(j)-b_mag(j-1)) c1=-c1/bmin*(1-2*epsb*kap(l)**2) d1=1-d1/bmin*(1-2*epsb*kap(l)**2) wgtba1(l,j)=2/c1*b1*sqrt(c1* r(j)+d1)-2/c1*b1*sqrt(c1*(r(j)-dtheta/2)+d1) a1=0. b1=1/gpar c1=(b_mag(j+1)-b_mag(j))/dtheta d1=b_mag(j)-r(j)/dtheta*(b_mag(j+1)-b_mag(j)) c1=-c1/bmin*(1-2*epsb*kap(l)**2) d1=1-d1/bmin*(1-2*epsb*kap(l)**2) wgtba1(l,j)=wgtba1(l,j)+2/c1*b1*sqrt(c1*(r(j)+dtheta/2)+d1)& -2/c1*b1*sqrt(c1* r(j)+d1) enddo ! ! near Bmin, use quadratic interpolation ! j=ldb/2 a1=0. b1=1/gpar c1=(b_mag(j)-b_mag(j-1))/dtheta d1=b_mag(j-1)-r(j-1)/dtheta*(b_mag(j)-b_mag(j-1)) c1=-c1/bmin*(1-2*epsb*kap(l)**2) d1=1-d1/bmin*(1-2*epsb*kap(l)**2) wgtba1(l,j)=2/c1*b1*sqrt(c1* r(j)+d1)-2/c1*b1*sqrt(c1*(r(j)-dtheta/2)+d1) b1=1/gpar a1=0. d1=b_mag(ldb/2+1) c1=(b_mag(ldb/2)-d1)/dtheta**2 c1=-c1/bmin*(1-2*epsb*kap(l)**2) d1=1-d1/bmin*(1-2*epsb*kap(l)**2) wgtba1(l,j)=wgtba1(l,j) & +b1/sqrt(-c1)*asin(sqrt(-c1/d1)*(r(j)+dtheta/2)) & -b1/sqrt(-c1)*asin(sqrt(-c1/d1)* r(j)) j=ldb/2+1 b1=1/gpar a1=0. d1=b_mag(ldb/2+1) c1=(b_mag(ldb/2)-d1)/dtheta**2 c1=-c1/bmin*(1-2*epsb*kap(l)**2) d1=1-d1/bmin*(1-2*epsb*kap(l)**2) wgtba1(l,j)= & b1/sqrt(-c1)*asin(sqrt(-c1/d1)*(r(j)+dtheta/2)) & -b1/sqrt(-c1)*asin(sqrt(-c1/d1)*(r(j)-dtheta/2)) j=ldb/2+2 b1=1/gpar a1=0. d1=b_mag(ldb/2+1) c1=(b_mag(ldb/2)-d1)/dtheta**2 c1=-c1/bmin*(1-2*epsb*kap(l)**2) d1=1-d1/bmin*(1-2*epsb*kap(l)**2) wgtba1(l,j)= & b1/sqrt(-c1)*asin(sqrt(-c1/d1)* r(j) ) & -b1/sqrt(-c1)*asin(sqrt(-c1/d1)*(r(j)-dtheta/2)) a1=0. b1=1/gpar c1=(b_mag(j+1)-b_mag(j))/dtheta d1=b_mag(j)-r(j)/dtheta*(b_mag(j+1)-b_mag(j)) c1=-c1/bmin*(1-2*epsb*kap(l)**2) d1=1-d1/bmin*(1-2*epsb*kap(l)**2) wgtba1(l,j)=wgtba1(l,j) & +2/c1*b1*sqrt(c1*(r(j)+dtheta/2)+d1) & -2/c1*b1*sqrt(c1* r(j) +d1) do j=ldb/2+3,ntheta0-1 a1=0. b1=1/gpar c1=(b_mag(j)-b_mag(j-1))/dtheta d1=b_mag(j-1)-r(j-1)/dtheta*(b_mag(j)-b_mag(j-1)) c1=-c1/bmin*(1-2*epsb*kap(l)**2) d1=1-d1/bmin*(1-2*epsb*kap(l)**2) wgtba1(l,j)=2/c1*b1*sqrt(c1* r(j) +d1) & -2/c1*b1*sqrt(c1*(r(j)-dtheta/2)+d1) a1=0. b1=1/gpar c1=(b_mag(j+1)-b_mag(j))/dtheta d1=b_mag(j)-r(j)/dtheta*(b_mag(j+1)-b_mag(j)) c1=-c1/bmin*(1-2*epsb*kap(l)**2) d1=1-d1/bmin*(1-2*epsb*kap(l)**2) wgtba1(l,j)=wgtba1(l,j)+2/c1*b1*sqrt(c1*(r(j)+dtheta/2)+d1) & -2/c1*b1*sqrt(c1* r(j) +d1) enddo j=ntheta0 a1=0. b1=1/gpar c1=(b_mag(j)-b_mag(j-1))/dtheta d1=b_mag(j-1)-r(j-1)/dtheta*(b_mag(j)-b_mag(j-1)) c1=-c1/bmin*(1-2*epsb*kap(l)**2) d1=1-d1/bmin*(1-2*epsb*kap(l)**2) wgtba1(l,j)=-2/c1*b1*sqrt(c1*(r(j)-dtheta/2)+d1) enddo ! ! passing ! do l=kd+1,kdpass-1 j=1 a1=0. b1=1/gpar c1=(b_mag(j)-b_mag(ldb))/dtheta d1=b_mag(ldb)-r(ldb)/dtheta*(b_mag(j)-b_mag(ldb)) c1=-c1/bmin*(1-2*epsb*kap(l)**2) d1=1-d1/bmin*(1-2*epsb*kap(l)**2) wgtba1(l,j)=2/c1*b1*sqrt(c1* r(j) +d1) & -2/c1*b1*sqrt(c1*(r(j)-dtheta/2)+d1) a1=0. b1=1/gpar c1=(b_mag(j+1)-b_mag(j))/dtheta d1=b_mag(j)-r(j)/dtheta*(b_mag(j+1)-b_mag(j)) c1=-c1/bmin*(1-2*epsb*kap(l)**2) d1=1-d1/bmin*(1-2*epsb*kap(l)**2) wgtba1(l,j)=wgtba1(l,j) & +2/c1*b1*sqrt(c1*(r(j)+dtheta/2)+d1) & -2/c1*b1*sqrt(c1* r(j) +d1) do j=2,ldb-1 a1=0. b1=1/gpar c1=(b_mag(j)-b_mag(j-1))/dtheta d1=b_mag(j-1)-r(j-1)/dtheta*(b_mag(j)-b_mag(j-1)) c1=-c1/bmin*(1-2*epsb*kap(l)**2) d1=1-d1/bmin*(1-2*epsb*kap(l)**2) wgtba1(l,j)=2/c1*b1*sqrt(c1* r(j) +d1) & -2/c1*b1*sqrt(c1*(r(j)-dtheta/2)+d1) a1=0. b1=1/gpar c1=(b_mag(j+1)-b_mag(j))/dtheta d1=b_mag(j)-r(j)/dtheta*(b_mag(j+1)-b_mag(j)) c1=-c1/bmin*(1-2*epsb*kap(l)**2) d1=1-d1/bmin*(1-2*epsb*kap(l)**2) wgtba1(l,j)=wgtba1(l,j) & +2/c1*b1*sqrt(c1*(r(j)+dtheta/2)+d1) & -2/c1*b1*sqrt(c1* r(j) +d1) enddo j=ldb a1=0. b1=1/gpar c1=(b_mag(j)-b_mag(j-1))/dtheta d1=b_mag(j-1)-r(j-1)/dtheta*(b_mag(j)-b_mag(j-1)) c1=-c1/bmin*(1-2*epsb*kap(l)**2) d1=1-d1/bmin*(1-2*epsb*kap(l)**2) wgtba1(l,j)=2/c1*b1*sqrt(c1* r(j) +d1) & -2/c1*b1*sqrt(c1*(r(j)-dtheta/2)+d1) a1=0. b1=1/gpar c1=(b_mag(1)-b_mag(j))/dtheta d1=b_mag(j)-r(j)/dtheta*(b_mag(1)-b_mag(j)) c1=-c1/bmin*(1-2*epsb*kap(l)**2) d1=1-d1/bmin*(1-2*epsb*kap(l)**2) wgtba1(l,j)=wgtba1(l,j) & +2/c1*b1*sqrt(c1*(r(j)+dtheta/2)+d1) & -2/c1*b1*sqrt(c1* r(j) +d1) enddo ! ! need special case for l=kdpass, where denominator is unity ! l=kdpass j=1 a1=0. b1=1/gpar wgtba1(l,j)=b1*dtheta/2 a1=0. b1=1/gpar wgtba1(l,j)=wgtba1(l,j)+b1*dtheta/2 do j=2,ldb-1 a1=0. b1=1/gpar wgtba1(l,j)=b1*dtheta/2 a1=0. b1=1/gpar wgtba1(l,j)=wgtba1(l,j)+b1*dtheta/2 enddo j=ldb a1=0. b1=1/gpar wgtba1(l,j)=b1*dtheta/2 a1=0. b1=1/gpar wgtba1(l,j)=wgtba1(l,j)+b1*dtheta/2 ! kap=0 limiting value: wgtba1(1,ldb/2+1)=pi/gpar*sqrt(bmin*dtheta**2/(b_mag(ldb/2)-bmin)) ! ! normalize to epsn/sqrt(2*epsb) to match igeo=0 results ! taub is normalized to q R_0/(v sqrt(2*epsb)) ! do l=1,kd ntheta0=ldb/2+l ntheta1=ldb/2+2-l do j=ntheta1,ntheta0 wgtba1(l,j)=wgtba1(l,j)*epsn*sqrt(2*epsb) enddo enddo do l=kd+1,kdpass do j=1,ldb wgtba1(l,j)=wgtba1(l,j)*epsn*sqrt(2*epsb) enddo enddo ! ! ntheta0,ntheta1 are chosen so that: ! For trapped particles, the bounce-average symmetrically involves ! phi points above and below the midplane. ! For passing particles, the bounce-average involves all phi points ! from l=1 to ldb (the ld point is identical to the l=1 point). ! ! calculate bounce-time numerically (for passing particles also...): do l=1,kd ntheta0=ldb/2+l ntheta1=ldb/2+2-l do j=ntheta1,ntheta0 taub(l)=taub(l)+wgtba1(l,j) enddo enddo do l=kd+1,kdpass taub(l)=0. do j=1,ldb taub(l)=taub(l)+wgtba1(l,j) enddo enddo ! calculate omegade numerically (for passing particles also...): ! This compiler directive is probably ignored by f90?? !ffp$ skip R omegade=0. do m = m_low, m_high do n = n_low, n_high do l=1,kd ntheta0=ldb/2+l ntheta1=ldb/2+2-l do j=ntheta1,ntheta0 omegade(l,m,n)=omegade(l,m,n)+wgtba1(l,j) & *(omega_gb(j,m,n,1)+omega_kap(j,m,n,1))/rho(1)/vt(1) omegade0(l,m,n)=omegade0(l,m,n)+wgtba1(l,j) & *(omega_gb0(j,m,n,1)+omega_kap0(j,m,n,1))/rho(1)/vt(1) omegade1(l,m,n)=omegade1(l,m,n)+wgtba1(l,j) & *(omega_gb1(j,m,n,1)+omega_kap1(j,m,n,1))/rho(1)/vt(1) enddo if ((m == meq).or.((md == 1).and.(mrr(m,n) == 0))) omegade(l,m,n)=0. omegade(l,m,n)=omegade(l,m,n)/taub(l) omegade0(l,m,n)=omegade0(l,m,n)/taub(l) omegade1(l,m,n)=omegade1(l,m,n)/taub(l) enddo enddo enddo !ffp$ skip R ! debug: if(idx_local(m_lo, 2, 1)) then do l=1,kdpass write(6,112) l,kap(l),taub(l),omegade(l,2,1)/epsn/rky(2,1) enddo endif 112 format( 'l,kap(l),taub(l),omegade(l,2,1)/rky(l,2,1)/epsn=',i3,3f14.7) ellk = taub/4 ! ! ?? I think I've now generalized all of this subroutine to use ! a grid-centered approach, with bounce-averages symmetrically ! involving phi above and below the midplane... ! wgtba2 = 0. do j=ldb/2+1,ld l0=j-ldb/2 do l=l0,kdpass if (l <= kd) then a1=1-b_mag(j)/bmin*(1-2*epsb*(kap(l)+dkap(l))**2) b1=1-b_mag(j)/bmin*(1-2*epsb*kap(l)**2) if (a1 < 0) then if(proc0) write(6,*) 'error in wgtba2 j,l,a1: ',j,l,a1 a1=0 endif if (b1 < 0) then if(proc0) write(6,*) 'error in wgtba2 j,l,b1: ',j,l,b1 b1=0 endif wgtba2(l,j)=sqrt(a1)-sqrt(b1) else a1=1-b_mag(j)/bmin*(1-2*epsb*kap(l)**2) b1=1-b_mag(j)/bmin*(1-2*epsb*(kap(l)-dkap(l))**2) if (a1 < 0) then if(proc0) write(6,*) 'error in wgtba2 j,l,a1: ',j,l,a1 a1=0 endif if (b1 < 0) then if(proc0) write(6,*) 'error in wgtba2 j,l,b1: ',j,l,b1 b1=0 endif wgtba2(l,j)=sqrt(a1)-sqrt(b1) endif wgtba2(l,j)=wgtba2(l,j)/dkap(l) ! B(theta)=1/(1+eps*cos(theta)) included to get particle conservation. ! symmetry, for theta<0: ! ?? GWH: also need to put a 1/B into the flux-surface averages... wgtba2(l,ldb-j+2)=wgtba2(l,j) enddo enddo ! ?? GWH: The bounce-averages (wgtba1), kappa-averages (wgtba2), ! and the form of the coefficients of the pitch-angle-scattering ! operator, are related in various ways. Some of these relations ! are important for particle conservation, and so it may be possible ! to derive some of these coefficients from each other in a way ! that preserves conservation properties exactly (or at least to ! within roundoff error). For now, just stick to the above methods, ! as they should preserve the conservation properties to second ! order in the grid spacings... ! BD: Since this is not used, I commented it out. ! load matrix for poisson inversion ! do j=1,ldb ! l1=1+abs(j-ldb/2-1) ! do l=l1,kdpass ! if ((l == l1).and.(l1 < kd)) then ! wt2(l,j)=dkap(l)*wgtba2(l,j)/2 ! elseif ((l > l1).and.(l < kd)) then ! wt2(l,j)=(dkap(l)*wgtba2(l,j)+dkap(l-1)*wgtba2(l-1,j))/2 ! elseif ((l == kd).and.(l1 /= kd)) then ! wt2(l,j)=dkap(l)*wgtba2(l,j)+dkap(l-1)*wgtba2(l-1,j)/2 ! elseif ((l == kd).and.(l1 == kd)) then ! wt2(l,j)=dkap(l)*wgtba2(l,j) ! elseif (l == kd+1) then ! wt2(l,j)=dkap(l)*wgtba2(l,j)+dkap(l+1)*wgtba2(l+1,j)/2 ! elseif ((l > kd+1).and.(l < kdpass)) then ! wt2(l,j)=(dkap(l)*wgtba2(l,j)+dkap(l+1)*wgtba2(l+1,j))/2 ! elseif (l == kdpass) then ! wt2(l,j)=dkap(l)*wgtba2(l,j)/2 ! endif ! enddo ! do j1=1,ldb ! pmat(j,j1)=0. ! l0=1+max(abs(j-ldb/2-1),abs(j1-ldb/2-1)) ! do l=l0,kdpass ! pmat(j,j1)=pmat(j,j1)+wt2(l,j)*wgtba1(l,j1)/taub(l) ! enddo ! temp(j,j1)=pmat(j,j1) ! enddo ! enddo ! c ! invert ! ifail=0 ! call f01aaE(temp,ld,ldb/2+1,pmat,ld,wrkinv,ifail) ! not used endif !(igeo if) ! initialization for electron collisions (ecollis) ! need to calculate a few bounce-averages: ! ! bratio = <(B-Bmin)/B> approx.= ! taub = bounce-time*v (only need the relative kappa-dependence) ! if (igeo == 0) then do l=1,kd bratio(l)=0.0 ntheta0=ldb/2+l ntheta1=ldb/2+2-l do j=ntheta1,ntheta0 bratio(l)=bratio(l)+cos(r(j))*wgtba1(l,j) enddo bratio(l)=bratio(l)/taub(l) enddo bratio(1)=0.0 do l=kd+1,kdpass bratio(l)=0.0 do j=1,ldb bratio(l)=bratio(l)+wgtba1(l,j)*cos(r(j)) enddo bratio(l)=bratio(l)/taub(l) enddo do l=2,kdpass bratio(l)=epsb*(1.-bratio(l)) enddo else do l=1,kd bratio(l)=0.0 ntheta0=ldb/2+l ntheta1=ldb/2+2-l do j=ntheta1,ntheta0 bratio(l)=bratio(l)+(b_mag(j)-bmin)/b_mag(j)*wgtba1(l,j) enddo bratio(l)=bratio(l)/taub(l) enddo bratio(1)=0.0 do l=kd+1,kdpass bratio(l)=0.0 do j=1,ldb bratio(l)=bratio(l)+(b_mag(j)-bmin)/b_mag(j)*wgtba1(l,j) enddo bratio(l)=bratio(l)/taub(l) enddo endif ! debug: if(proc0) then do l=1,kdpass write(6,103) l,bratio(l),taub(l) 103 format('l, bratio(l), taub(l), =',i3,2f14.7) enddo endif ! ! The generic diffusion operator can be written in the form: ! ! df/dt = J(x) d/dx ( D(x) df/dx) ! ! Where J(x)=1/V'(x) is a Jacobian which represents conservation properties. ! In our case, x is the the pitch-angle kappa. The bounce-averaged ! collision operator we use was calculated by J.G. Cordey, Nucl. Fus. 16, ! 499 (976). The notation I use is described on p.69-70 of my thesis, ! "Fast Ion Studies of Ion Cyclotron Heating in the PLT Tokamak", ! G.W. Hammett, Princeton University Ph.D. Disseration, 1986. ! ! Use conservative differencing formulas, including the volume elements ! in the coefficients: ! do l=2,kdpass-1 jac(l)=1.0/kap(l)/taub(l)/ ((kap(l+1)-kap(l-1))/2) enddo jac(kdpass)=1.0/kap(l)/taub(l)/ ((kap(kdpass)-kap(kdpass-1))) ! deal with kap(1)=0.0, by using kap evaluated half-a-grid away: kapav=kap(2)/4.0 jac(1)=1./kapav/taub(1)/ (kap(2)/2) do l=1,kdpass-1 kapp=(kap(l)+kap(l+1))/2 ! kap2p=(kap(l)**2+kap(l+1)**2)/2 kap2p=kapp**2 taubp=(taub(l)+taub(l+1))/2 ! bratiop=(bratio(l)+bratio(l+1))/2 bratiop=(bratio(l)**0.5+bratio(l+1)**0.5)/2 bratiop=bratiop**2 ! the reason we interpolate sqrt(bratio) instead of bratio, is that bratio ! goes like kap**2 near kap=0. There are near-cancellations between ! the kap**2 and bratio terms, and one must be careful to interpolate ! both of them in the same way or the diffusion coefficient could go ! negative near kap=0... diff(l)=(1-2*epsb*kap2p)*taubp/kapp*(-bratiop+2*epsb*kap2p)/(kap(l+1)-kap(l)) enddo diff(0)=0.0 ! guard cell value ! evaluate at kap=0 point (taub logarithmically --> infinity at kap=1): ! diff(kd)=(1-2*epse*kap(kd)**2)*taub(kd)/kap(kd)* & ! (-bratio(kd)+2*epse*kap(kd)**2)/(1.0-kap(kd)) ! no diffusion of particles beyond maximum kappa grid point: diff(kdpass)=0.0 ! debug: if(proc0) then write(6,*) '0, diff(0+0.5),jac =',0,diff(0) do l=1,kdpass write(6,*) 'l, diff(l+0.5),jac =',l,diff(l),jac(l) enddo endif deallocate (ellk) end subroutine init_bounce_avg subroutine alloc_bounce_avg use itg_data, only: epse use gryffin_grid, only: kdpass, kap if (epse > 0.0) then allocate (kap(kdpass+1), dkap(kdpass), taub(kdpass), bratio(kdpass)) allocate (jac(kdpass), diff(0:kdpass)) kap = 0.; dkap = 0.; taub = 0.; bratio = 0.; jac = 0.; diff = 0. endif end subroutine alloc_bounce_avg end module bounce_avg