program ball
use geometry, only: eikcoefs, eikcoefs_output_type
!> Inputs
use geometry, only: eqfile, equal_arc, use_large_aspect, surf
use geometry, only: bishop, irho, force_sym, efit_eq, writelots, big
use geometry, only: gen_eq, ppl_eq, local_eq, chs_eq, gs2d_eq, transp_eq
use geometry, only: alpha_input, dp_mult, p_prime_input, invLp_input, beta_prime_input
implicit none
real, dimension(:), pointer :: psi
real, dimension(:), allocatable :: theta
real :: dbetadrho
integer :: nbeta, ntgrid, ntheta, iunstable, ibeta, ishat
integer :: nshat, order, isym, iflux, itor, nperiod
real :: beta_p1, beta_p2, diffscheme, psiend, beta_prime, cflmax, cflmin
real :: beta_prime_times, beta_prime_over
real :: bprev, pprev, stabmhd, bcrit, s_hat_max, shat_0
real :: akappa, tri, akappri, tripri, rmin, rmax
real :: rhoc, rmaj, r_geo, qinp, s_hat_input, geoType, shift, delrho
logical :: diagram
type(eikcoefs_output_type) :: eikcoefs_results
namelist/stuff/ntheta,nperiod,geoType,rmaj,akappri,akappa,shift,equal_arc, &
rhoc,rmin,rmax,itor,qinp,iflux,delrho,tri,bishop, irho, isym, &
tripri,writelots,R_geo, gen_eq,efit_eq,local_eq,eqfile, &
chs_eq, s_hat_input,p_prime_input,invLp_input,ppl_eq,&
diffscheme,nbeta,beta_p1,beta_p2,beta_prime_input, &
alpha_input, beta_prime_times, beta_prime_over, big, gs2d_eq, &
transp_eq, diagram, nshat, s_hat_max, force_sym
! set defaults
ntheta=32; nperiod=1
equal_arc = .false. ; bishop = 0
eqfile='dskeq.cdf'
rhoc=0.5
rmin = -1.0 ; rmax = -1.0
diagram = .false.
nshat = 1 ; s_hat_max = 2.0
geoType = 0
rmaj = 3 ; akappa = 1 ; tri = 0 ; qinp = 2
R_geo = 3 ; akappri = 0 ; tripri=0 ; shift = 0.
s_hat_input = 1
p_prime_input = 1
invLp_input = 5
beta_prime_times = -1. ; beta_prime_over = -1.
irho = 2 ; isym = -1 ; itor = -1; iflux=-1
delrho = 0.01
diffscheme=0.33; nbeta = 1; beta_p1 = 0.; beta_p2 = 0.
! read in data
open(10,file='eik.in')
read(10,stuff)
close(10)
if (local_eq .and. (geoType /= 0)) print*,'Warning : geoType /= 0 not supported'
if (rmin /= -1) write(*,*) "Option deprecated : rmin set but not used"
if (rmax /= -1) write(*,*) "Option deprecated : rmax set but not used"
if (isym /= -1) write(*,*) "Option deprecated : isym set but not used (see force_sym)"
if (iflux /= -1) write(*,*) "Option deprecated : iflux set but not used"
if (itor /= -1) write(*,*) "Option deprecated : itor set but not used (see use_large_aspect)"
writelots = .false.
surf%geoType = geoType
surf%mMode = 2 ; surf%nMode = 3
surf%rmaj = rmaj ; surf%R_geo = R_geo ; surf%r = rhoc ; surf%dr = delrho
surf%sHorz = shift ; surf%delm = akappa
surf%deln = tri ; surf%delmp = akappri
surf%delnp = tripri ; surf%thm = 0
surf%thn = 0 ; surf%q = qinp ; surf%shat = s_hat_input
if(.not.local_eq) then
if(beta_prime_times >= 0) then
if(beta_prime_over == -1) beta_prime_over = beta_prime_times
endif
endif
equal_arc = .false. ! this should not be required
! Note that if .not. local_eq then always choose use_large_aspect = F
if ((.not. local_eq) .and. use_large_aspect) then
use_large_aspect = .false.
write(*,*) 'Forcing use_large_aspect to false'
end if
! Catch if no q profile
if (local_eq .and. irho == 1) then
irho = 2
write(*,*) 'Forcing irho=2'
end if
if (local_eq) then
if(gen_eq) write(*,*) 'Forcing gen_eq to be false because local_eq is set'
gen_eq = .false.
end if
! compute the theta grid
if( (.not. gen_eq) .and. (.not. ppl_eq) &
.and. (.not. transp_eq) .and. (.not. chs_eq)) then
call eikcoefs(ntheta, nperiod, eikcoefs_results)
theta = eikcoefs_results%theta
dbetadrho = eikcoefs_results%dbetadrho
ntgrid = (2*nperiod - 1)*(ntheta/2)
allocate(psi(-ntgrid:ntgrid))
else
allocate(psi(-1:1))
endif
open(11,file='ball.out')
open(29,file='psi.out')
if (.not. diagram) then
order = 1
call beta_loop (bcrit, order, ntheta)
stop
endif
! Attempt to sketch out an s-alpha type diagram
if (bishop /= 8) stop ! just do one case for now
shat_0 = s_hat_input
do order = 1, 2
do ishat = 1, nshat
if (nshat > 1) &
s_hat_input = shat_0 + (ishat-1)*(s_hat_max-shat_0)/(nshat-1)
bcrit = 0.
call beta_loop (bcrit, order, ntheta)
if (bcrit /= 0.) then
write (*,*) 'shat= ',s_hat_input,' alpha= ', -Rmaj*qinp**2*bcrit
end if
end do
write (*,*)
end do
contains
subroutine beta_loop (bcrit, order, ntheta)
integer, intent(inout) :: ntheta
integer :: order
real :: bcrit
logical :: done
integer :: ibmin, ibmax, ibstep
done = .false.
if (order == 1) then
ibmin = 1
ibmax = nbeta
ibstep = 1
else
ibmin = nbeta
ibmax = 1
ibstep = -1
end if
do ibeta=ibmin, ibmax, ibstep
if(bishop >= 7) then
if(nbeta == 1) then
dp_mult = beta_prime_times
else
dp_mult = 1./beta_prime_over+(ibeta-1) &
*(beta_prime_times-1./beta_prime_over)/(nbeta-1)
endif
else
if(nbeta == 1) then
beta_prime=beta_p1
else
beta_prime=beta_p1+(ibeta-1)*(beta_p2-beta_p1)/(nbeta-1)
endif
beta_prime_input = beta_prime
endif
call eikcoefs(ntheta, nperiod, eikcoefs_results)
theta = eikcoefs_results%theta
dbetadrho = eikcoefs_results%dbetadrho
if(.not.local_eq) then
ntgrid = (size(theta)-1)/2
ntheta = (size(theta)-1)/(2*nperiod-1)
deallocate(psi)
allocate(psi(-ntgrid:ntgrid))
endif
if(bishop>= 7) beta_prime = dbetadrho
call mhdballn(ntgrid, beta_prime, diffscheme, iunstable, psiend, &
cflmax, cflmin, psi, eikcoefs_results%bmag, eikcoefs_results%cvdrift, &
eikcoefs_results%gradpar, eikcoefs_results%gds2, theta)
if(ibeta.eq.1) then
stabmhd=0.
else
stabmhd=0.
if(psiend /= pprev) then
stabmhd=beta_prime-psiend*(beta_prime-bprev)/(psiend-pprev)
end if
if (diagram) then
if (iunstable == 1 .and. .not. done) then
bcrit = stabmhd
done = .true.
end if
end if
endif
bprev=beta_prime
pprev=psiend
if (done) cycle
if (.not. diagram) then
write(6,998) beta_prime,iunstable,stabmhd,max(abs(cflmax),abs(cflmin)),dp_mult
end if
enddo
998 format(1x,g13.6,2x,i2,3x,g13.6,2x,2(1pe11.4,1x))
end subroutine beta_loop
subroutine mhdballn(ntgrid, beta_prime, diffscheme, iunstable, psiend, cflmax, &
cflmin, psi, bmag, cvdrift, gradpar, gds2, theta)
implicit none
integer, intent(in) :: ntgrid
real, intent(in) :: beta_prime, diffscheme
integer, intent(out) :: iunstable
real, intent(out) :: psiend, cflmax, cflmin
real, dimension(-ntgrid:ntgrid), intent(out) :: psi
real, dimension(-ntgrid:ntgrid), intent(in) :: bmag, cvdrift, gradpar, gds2, theta
integer :: i
real, dimension(-ntgrid:ntgrid) :: psi_p, c, g, c1, c2, ch, g1, g2, gh, delx
real :: alpha, b_prime
! include 'ball.inc'
! input arguments:
! npt=total number of gridpoints with theta>0 = nperiod*ntheta+1
! beta_prime=d beta/ d rho
! alpha=numerical knob (0 to .333 are good if cflmax,min<1)
!
! output arguments:
! iunstable=number of times the solution crossed the axis. If
! this is greater than 0 it is unstable.
!
! psi is the solution of the Euler Lagrange equation for delta W.
! At the first theta point, psi is zero with slope positive 1.
!
! psiend=psi at the last gridpoint
! psi is psi at every grid point
! cflmax,min are indicators of whether resolution is adequate.
! diffscheme is a numerical difference scheme knob.
! For greatest accuracy diffscheme=0. or .333333 .
! If the absolute value of cflmax,min approaches one, inaccuracy and/or
! numerical instability are possible. If cflmax,min are >1-3, increase
! diffscheme to .666 or 1.0 for stability, but probably resolution is
! inadequate.
!
b_prime = -beta_prime
iunstable=0
alpha=1.-diffscheme
g = gds2*gradpar/bmag
c = b_prime*cvdrift/(2.*bmag*gradpar)
i=-ntgrid
cflmax=(theta(i+1)-theta(i))**2*(c(i+1)+c(i))/(g(i+1)+g(i))
cflmin=cflmax
do i=-ntgrid+1, ntgrid
delx(i)=theta(i)-theta(i-1)
ch(i)=.5*(c(i)+c(i-1))
gh(i)=.5*(g(i)+g(i-1))
cflmax=max(cflmax,delx(i)**2*ch(i)/gh(i))
cflmin=min(cflmin,delx(i)**2*ch(i)/gh(i))
enddo
do i=-ntgrid+1, ntgrid-1
c1(i)=-delx(i+1)*(alpha*c(i)+.5*diffscheme*ch(i+1)) &
-delx(i )*(alpha*c(i)+.5*diffscheme*ch(i )) &
-delx(i)*.5*diffscheme*ch(i)
c1(i)=.5*c1(i)
enddo
do i=-ntgrid+1, ntgrid
c2(i)=-.25*diffscheme*ch(i)*delx(i)
g1(i)=gh(i)/delx(i)
g2(i)=1./(gh(i)/delx(i)+.25*diffscheme*ch(i)*delx(i))
enddo
i=-ntgrid
psi(i)=0.
i=i+1
psi(i)=delx(i)
psi_p(i)=psi(i)/g2(i)-g1(i)*psi(i-1)
do i=-ntgrid+1,ntgrid-1
psi_p(i+1)=psi_p(i)+c1(i)*psi(i)+c2(i)*psi(i-1)
psi(i+1)=(g1(i+1)*psi(i)+psi_p(i+1))*g2(i+1)
enddo
psiend=psi(ntgrid)
do i=-ntgrid+1,ntgrid-1
if(psi(i)*psi(i+1) <= 0.) iunstable=iunstable+1
enddo
end subroutine mhdballn
end program ball