!> FIXME : Add documentation 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, geoType 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, 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