FIXME : Add documentation
Type | Intent | Optional | Attributes | Name | ||
---|---|---|---|---|---|---|
integer, | intent(inout) | :: | ntheta_input |
Sets the number of theta grid points. If odd will be updated to be ntheta_input - 1 |
||
integer, | intent(in) | :: | nperiod | |||
type(eikcoefs_output_type), | intent(out) | :: | outputs |
A type containing all of the outputs of eikcoefs |
||
integer, | intent(in), | optional | :: | job_id_in |
job_id is the job id in a Trinity or list mode run. Only for debug |
subroutine eikcoefs (ntheta_input, nperiod, outputs, job_id_in)
use constants, only: twopi
use file_utils, only: open_output_file, close_output_file, finish_file_utils
use file_utils, only: futils_initialized => initialized, init_file_utils
use optionals, only: get_option_with_default
implicit none
!> Sets the number of theta grid points. If odd will be updated to be ntheta_input - 1
integer, intent(in out) :: ntheta_input
integer, intent(in) :: nperiod
!> A type containing all of the outputs of eikcoefs
type(eikcoefs_output_type), intent(out) :: outputs
!> job_id is the job id in a Trinity or list mode run. Only for debug
integer, intent(in), optional :: job_id_in
logical :: start_file_utils, dummy_flag
integer :: job_id, i, iBmax, fort11_unit
real, allocatable, dimension(:) :: seik, rgrid, invR, trip, dsdthet, rmajor, dummy
real, allocatable, dimension(:,:) :: thgrad, rpgrad, crpgrad, grads
real, allocatable, dimension(:,:) :: gradrptot, gradstot, gradztot, gradthtot, gradbtot
real :: shat, dbetadrho, dqdrp, qval, dum, bi
real :: rp, pbar, drhodrp, drhodrhod, drhodpsi, dpsidrp, dpsidrho
! Used to translate definition of the poloidal angle such that the location
! of the maximium magnetic field is at pi (see section 3.2.3 of Ball MIT
! Masters thesis or section 2.4 of "GS2 analytic geometry specification")
real :: thetaShift
job_id = get_option_with_default(job_id_in, -1)
start_file_utils = .not. futils_initialized
if (start_file_utils) then
! Initialize file utils with a run name we provide
! and by setting trin_run prevent it from searching
! the command line arguments.
write(*,*) 'WARNING: Automatically initializing file_utils&
& in eikcoefs. This is a hack to &
& support older programs like eiktest. You should&
& initialize file utils before calling eikcoefs'
call init_file_utils(dummy_flag, name='geometry_module', trin_run=.true., input=.false.)
end if
!Create module level output file and get unit
call open_output_file(fort11_unit,".fort.11")
!Compute the initial constants
debug = (verb > 2)
if (debug) writelots = .true.
if (debug) write(6,*) "eikcoefs: local_eq=",local_eq, 'jid=',job_id
if (local_eq) then
if (irho /= 2) write(*,*) 'Forcing irho = 2 for local_eq'
irho = 2
if (bishop == 0) then
write(*,*) 'Bishop = 0 not compatible with local_eq. Forcing bishop = 1'
bishop = 1 !Not clear that this is a good choice for Miller
end if
end if
if (debug) write(6,*) "eikcoefs: call check ",'jid=',job_id
call check(gen_eq, efit_eq, ppl_eq, local_eq, dfit_eq, idfit_eq, chs_eq, gs2d_eq, transp_eq, bishop)
call init_uniform_theta_grid(ntheta_input, outputs%theta, nperiod, ntgrid, nth)
surf%nt = ntheta_input
if (debug) write(6,*) 'eikcoefs: alloc_local_arrays'
if (.not. allocated(seik)) call alloc_local_arrays(ntgrid)
thetaShift = 0.0
allocate(outputs%bmag(-ntgrid:ntgrid), outputs%bpol(-ntgrid:ntgrid))
if (.not. allocated(geom)) call allocate_geom(geom)
call initialise_geometry_and_get_bmag(geom, surf%r, thetaShift, outputs%theta, outputs%bmag, outputs%bpol)
! Start theta shift calc. This is used to modify the definition of
! theta such that the maxiumum magnetic field is at pi. This is
! necessary to correctly calculate bounce points. Can only shift with local_eq
if (local_eq) then
! Locate the index of the maximum magnetic field in the range [-pi,pi)
iBmax = maxloc(outputs%bmag(-nth:nth-1), dim = 1) - (nth + 1)
! If max is not at -pi and theta then shift theta.
if (iBmax /= -nth) then
thetaShift = outputs%theta(iBmax) - outputs%theta(-nth)
call geom%finalise !Clean up geometry to allow reinitialisation
write(*,*) "Bmax is not located at +/-pi, translating theta by",thetaShift," so that it is..."
call initialise_geometry_and_get_bmag(geom, surf%r, thetaShift, outputs%theta, outputs%bmag, outputs%bpol)
end if
end if
! compute |grad rho|:
allocate(outputs%grho(-ntgrid:ntgrid))
outputs%grho(-nth:nth) = hypot(rpgrad(-nth:nth,1), rpgrad(-nth:nth,2)) * abs(drhodrp)
call periodic_copy(outputs%grho, 0., nperiod)
if (force_sym) call sym(outputs%grho, 0, ntgrid)
if (writelots) write(fort11_unit,*) 'dpsidrp= ',dpsidrp
drhodpsi = drhodrp / dpsidrp ; dpsidrho = 1.0 / drhodpsi
! compute derivatives of eikonal, grad B and grad Z
if (bishop /= 0) then
call bishop_grads(surf%r, surf%rmaj, outputs%bmag, outputs%bpol)
else
call bishop_large_aspect_ratio_approximation(surf%r, surf%dr, ntheta_input)
end if
if (writelots) then
write(fort11_unit,*) 'i,theta,grads1,2= '
do i = -ntgrid, ntgrid
write(fort11_unit,*) i, outputs%theta(i), grads(i, 1), grads(i, 2)
end do
end if
! compute total grad S including toroidal part
call gradstottgrid(invR, grads, gradstot)
gradrptot(-nth:nth, 1:2) = rpgrad(-nth:nth, :)
call periodic_copy(gradrptot(:, 1), 0., nperiod)
call periodic_copy(gradrptot(:, 2), 0., nperiod)
gradrptot(:, 3) = 0.
! gradthtot is grad theta since ballooning, theta goes from
! -ntheta*nperiod -> ntheta*nperiod, but grad theta is periodic in theta
gradthtot(-nth:nth, 1:2) = thgrad(-nth:nth, :)
call periodic_copy(gradthtot(:, 1), 0., nperiod)
call periodic_copy(gradthtot(:, 2), 0., nperiod)
gradthtot(:, 3) = 0.
! compute the magnitude squared of grad S
allocate(outputs%gds2(-ntgrid:ntgrid), outputs%gds21(-ntgrid:ntgrid), &
outputs%gds22(-ntgrid:ntgrid), outputs%gds23(-ntgrid:ntgrid), &
outputs%gds24(-ntgrid:ntgrid), outputs%gds24_noq(-ntgrid:ntgrid))
call calculate_metric_terms(gradstot, gradrptot, gradthtot, dpsidrho, dqdrp, &
qval, surf%r, outputs%bmag, force_sym, outputs%gds2, outputs%gds21, outputs%gds22, &
outputs%gds23, outputs%gds24, outputs%gds24_noq)
! compute curvature, grad b, and coriolis drift terms
allocate(outputs%gbdrift(-ntgrid:ntgrid), outputs%gbdrift0(-ntgrid:ntgrid))
allocate(outputs%cvdrift(-ntgrid:ntgrid), outputs%cvdrift0(-ntgrid:ntgrid))
allocate(outputs%cdrift(-ntgrid:ntgrid), outputs%cdrift0(-ntgrid:ntgrid))
call drift(rgrid, rp, gradstot, gradrptot, gradztot, gradbtot, dqdrp, dpsidrho, dpsidrp,&
outputs%bmag, invR, outputs%theta, dbetadrho, bi, force_sym, nperiod, &
outputs%gbdrift, outputs%gbdrift0, outputs%cvdrift, outputs%cvdrift0, &
outputs%cdrift, outputs%cdrift0)
!Recompute the magnetic field along theta if use_large_aspect - should this come earlier?
if (use_large_aspect) then
outputs%bmag(-nth:nth)= 1 / (1 + rgrid(-nth:nth) * cos(outputs%theta(-nth:nth)) / surf%rmaj)
call periodic_copy(outputs%bmag, 0., nperiod)
end if
if (force_sym) call sym(outputs%bmag, 0, ntgrid)
! compute the grad-parallel coefficient
allocate(outputs%gradpar(-ntgrid:ntgrid))
if (dfit_eq) then
outputs%gradpar = trip / outputs%bmag
else
outputs%gradpar(-nth:nth) = -bi / (rmajor(-nth:nth)**2 * outputs%bmag(-nth:nth) * dsdthet(-nth:nth))
call periodic_copy(outputs%gradpar, 0., nperiod)
end if
if (force_sym) call sym(outputs%gradpar, 0, ntgrid)
! Note : We should call this before we potentially redefine the theta grid
! due to equal_arc below as some specific equilibrium types rely on theta
! having a consistent definition with that used to calculate rgrid initially.
! In particular this is important for Efit and GS2D equilibria (in eeq). Other
! equilibrium types use interpolation so shouldn't depend on the specific theta grid
! (although it's possible there could still be some interaction, given that one of the
! outputs of the below 'aplot' is given by `seik` which is evaluated on the original
! grid rather than the "new" one defined below).
! @note Check if restriction on cartesian geo types still holds.
allocate(outputs%rplot(-ntgrid:ntgrid), outputs%rprime(-ntgrid:ntgrid))
allocate(outputs%zplot(-ntgrid:ntgrid), outputs%zprime(-ntgrid:ntgrid))
allocate(outputs%aplot(-ntgrid:ntgrid), outputs%aprime(-ntgrid:ntgrid))
call plotdata (rgrid, seik, grads, dpsidrho, surf%dr, outputs%theta, &
outputs%rplot, outputs%zplot, outputs%aplot, &
outputs%rprime, outputs%zprime, outputs%aprime)
! In the lines below we redefine theta to be theta', and gradpar
! to be gradpar'. This means that the values of theta in the gs2
! output file are not values of the poloidal angle but values of theta_prime
if (equal_arc) then
allocate(dummy(-ntgrid:ntgrid))
call arclength (ntheta_input, nperiod, outputs%theta, outputs%gradpar, dummy)
outputs%theta = dummy
end if
call calculate_surface_area(outputs%theta, outputs%gradpar, outputs%bmag, outputs%grho, &
dpsidrho, surf%r, nth, outputs%surfarea, outputs%dvdrhon)
if (writelots) write (fort11_unit,*) 'surfarea=', outputs%surfarea
!Close output file
call close_output_file(fort11_unit)
outputs%kxfac = abs(qval * dpsidrho / surf%r)
if (qval == 0. .and. (dfit_eq .or. idfit_eq)) then
shat = max(1.e-7, shat)
!RN> I think kxfac should include the sign. However, from GS2 convention,
!RN> eventually kxfac = abs(kxfac). See get_grids in theta_grid. Even if the sign is
!RN> effective, probably it has no effect due to symmetry. I just leave the sign here.
outputs%kxfac = - dpsidrho
outputs%gds22 = (outputs%grho * dpsidrho)**2
end if
if (start_file_utils) call finish_file_utils
! Pack data into output type
outputs%ntheta = ntheta_input ; outputs%drhodpsi = drhodpsi ; outputs%rhoc = surf%r
outputs%qsf = qval ; outputs%shat = shat ; outputs%dbetadrho = dbetadrho
contains
!> Initialises the specific geometry backend
subroutine initialise_geometry_and_get_bmag(geom, rhoc, thetaShift, theta, bmag, bpol)
use geo_utils, only: geo_input_type
class(abstract_geo_type), intent(in out) :: geom
real, intent(in) :: rhoc, thetaShift
real, dimension(-ntgrid:), intent(in) :: theta
real, dimension(-ntgrid:), intent(out) :: bmag, bpol
type(geo_input_type) :: geom_inputs
real :: avgrmid, rbar, B_T0
integer :: fort78_unit, fort99_unit
! Initialise the specific geometry types (mostly reading data
! from file and calculating derived quantities).
geom_inputs%surf = surf ; geom_inputs%thShift = thetaShift ; geom_inputs%big = big
geom_inputs%eqfile = eqfile ; geom_inputs%eqnormfile = eqnormfile ; geom_inputs%theta = theta
call geom%initialise(geom_inputs, rpmin, rpmax, surf%rmaj, B_T0, avgrmid)
if (debug) write(6,*) "eikcoefs: rpmin,rpmax=", rpmin, rpmax
if (.not. local_eq) surf%r_geo = geom%rcenter(rpmax) !Only used for writelots/checks
if (surf%r_geo > surf%rmaj + 1.e-5) then
write(6,*)
write(6,*) 'Warning: the center of the LCFS is outboard of the '
write(6,*) 'center of the reference surface. This is probably wrong.'
write(6,*)
end if
if (writelots) then
write(fort11_unit,*) 'All lengths normalized to the avg midplane minor radius'
if (irho /= 2) then
write(fort11_unit,*) 'except rho, which always goes from 0 to 1.'
write(fort11_unit,*) 'The rho normalization is determined by irho.'
end if
write(fort11_unit,*) 'avgrmid = ',avgrmid,' meters'
if (.not. local_eq) write(fort11_unit,*) 'R_mag = ',surf%rmaj
if (.not. local_eq) write(fort11_unit,*) 'B_T0 = ',B_T0
end if
if (debug) write(6,*) "eikcoefs: find rp ", 'jid=',job_id
rp = rpofrho(rhoc) ; drhodrp = drho_drp(rp, surf%dr)
if (debug) write(6,*) "eikcoefs: rp=",rp, 'jid=',job_id
if (debug) write(6,*) "eikcoefs: find rgrid"
call rtgrid(rgrid, rp, theta)
call periodic_copy(rgrid, 0.0, nperiod)
pbar = pbarfun(rgrid(0), theta(0))
bi = geom%btori(pbar) ! Get the toroidal field function (i.e. fp = R B_T)
if (writelots) then
write(fort11_unit,*) 'rp of rho = ',rp
write(fort11_unit,*) 'Rcenter(rho) = ',geom%rcenter(rp)
write(fort11_unit,*) 'R_geo = ',surf%r_geo,' !or ',bi
call open_output_file(fort78_unit,".fort.78")
call open_output_file(fort99_unit,".fort.99")
write(fort78_unit,*) '# theta_pol, R, Z'
do i = -nth, nth
write(fort78_unit,*) theta(i), geom%Rpos(rgrid(i), theta(i)), &
geom%Zpos(rgrid(i), theta(i))
write(fort99_unit,*) 'i=', i, ' thet= ', theta(i), ' rgrid= ', rgrid(i)
end do
call close_output_file(fort78_unit) ; call close_output_file(fort99_unit)
rbar = 0.
do i = -nth, nth - 1
rbar = rbar + rgrid(i) * (theta(i + 1) - theta(i))
end do
rbar = rbar / twopi
write (*,*) 'r_bar = ', rbar
write (*,*)
write(fort11_unit,*) 'drhodrp=', drhodrp
if (irho == 1) then
drhodrhod = drho_drhod(rp, drhodrp, surf%dr)
write(fort11_unit,*) 'drho_drhod=', drhodrhod
end if
end if
if (debug) write(fort11_unit,*) 'eikcoefs: various logs', gen_eq, ppl_eq, chs_eq, efit_eq, dfit_eq, idfit_eq, gs2d_eq, transp_eq
if (debug) write(6,*) "eikcoefs: call rmajortgrid"
call rmajortgrid(rgrid, theta, rmajor) ; call periodic_copy(rmajor, 0., nperiod)
if (use_large_aspect) then
invR = 1.0 / surf%rmaj
else
call invRtgrid(rgrid, theta, invR) ; call periodic_copy(invR, 0., nperiod)
end if
! compute gradient of rp in R, Z, phi (cartesian - crpgrad)
! Note that rpgrad is given in (R,Z,phi) coordinates if bishop == 0
! and (rho,l,phi) coordinates otherwise
! Also, rpgrad is grad psi if not local_eq and grad rho if local_eq -- MAB
! When bishop > 0
! rpgrad(:,1) is |grad rp| ; rpgrad(:,2) is 0
! crpgrad(:,1) is drp/dR ; crpgrad(:,2) is drp/dZ
!
! bishop = 0 is high aspect ratio coefficients for debugging -- EGH
call geom%gradient(rgrid, theta, crpgrad, 'P', dum, nth, ntgrid, .false.)
call periodic_copy(crpgrad(:,1), 0., nperiod)
call periodic_copy(crpgrad(:,2), 0., nperiod)
! compute coordinate gradients
call geom%gradient(rgrid, theta, rpgrad, 'P', dum, nth, ntgrid, bishop /= 0)
call geom%gradient(rgrid, theta, thgrad, 'T', dum, nth, ntgrid, bishop /= 0)
call tripprod2dtgrid(rpgrad, thgrad, invR, trip)
call periodic_copy(trip, 0., nperiod)
! compute eikonal S
if (local_eq) qval = geom%qfun(pbar) !leq:qfun ignores input pbar
if (debug) write(6,*) "eikcoefs: call eikonal"
call eikonal(theta, invR, trip, qval, bi, seik, dsdthet, dpsidrp)
if (debug) write(6,*) "eikcoefs: done eikonal"
if (writelots) write(fort11_unit,*) 'q= ', qval
if (local_eq) trip = trip * dpsidrp
if (bishop > 0) then !Following bits _could_ be calculated for all bishop
call B_pol(invR, crpgrad, dpsidrp, bpol, nth) !<Should this use rpgrad for local?
call periodic_copy(bpol, 0., nperiod)
call B_mod(invR, bpol, bmag, bi, nth)
call periodic_copy(bmag, 0., nperiod)
end if !Note bmag not set for bishop == 0
end subroutine initialise_geometry_and_get_bmag
subroutine bishop_grads(rhoc, rmaj, bmag, bpol)
real, intent(in) :: rhoc, rmaj
real, dimension(-ntgrid:), intent(in) :: bmag, bpol
logical, parameter :: save_bishop = .true.
integer :: bishop_unit
real :: a_b, b_b, c_b, dp, dp_new, di, di_new, s_hat, s_hat_new, pressure, tmp
real, allocatable, dimension(:) :: a_bish, b_bish, c_bish, da_bish, db_bish, dc_bish, &
zoftheta, dzdl, dsdl, ltheta, th_bish, invRpol
allocate(a_bish(-ntgrid:ntgrid), b_bish(-ntgrid:ntgrid), c_bish(-ntgrid:ntgrid), &
da_bish(-ntgrid:ntgrid), db_bish(-ntgrid:ntgrid), dc_bish(-ntgrid:ntgrid), &
zoftheta(-ntgrid:ntgrid), dzdl(-ntgrid:ntgrid), dsdl(-ntgrid:ntgrid), &
ltheta(-ntgrid:ntgrid), th_bish(-ntgrid:ntgrid), invRpol(-ntgrid:ntgrid))
call th_bishop(crpgrad, th_bish, nth)
call periodic_copy(th_bish, 0., nperiod)
call loftheta(rgrid, outputs%theta, ltheta)
call periodic_copy(ltheta, ltheta(nth) - ltheta(-nth), nperiod)
call inv_Rpol(outputs%theta, th_bish, ltheta, invRpol, nth)
call periodic_copy(invRpol, 0., nperiod)
if (debug) write(6,*) "eikcoefs: gradients"
if (debug) write(6,*) 'drhodpsi', drhodpsi
if (debug) write(6,*) 'bi', bi
if (verb > 5) write(6,*) 'rmajor', rmajor(:)
if (verb > 5) write(6,*) 'bpol', bpol(:)
if (verb > 5) write(6,*) rpgrad(-nth:nth,1)
if (verb > 5) write(6,*) rpgrad(-nth:nth,2)
if (verb > 5) write(6,*) thgrad(-nth:nth,1)
if (verb > 5) write(6,*) thgrad(-nth:nth,2)
if (verb > 5) write(6,*) 1.0 / trip(-nth:nth)
if (verb > 5) write(6,*) -1.0 / (bpol(-nth:nth) * invRpol(-nth:nth))
if (debug) write(6,*) "eikcoefs: end gradients"
da_bish = (1 / rmajor**2) + (bi / (bpol * rmajor**2))**2
db_bish = bi / (bpol * rmajor)**2
dc_bish = bi * (sin(th_bish) + rmajor * invRpol) / (-bpol * rmajor**4)
da_bish = da_bish / trip
db_bish = db_bish / trip
dc_bish = dc_bish / trip
call integrate(da_bish, outputs%theta, a_bish, ntgrid)
call integrate(db_bish, outputs%theta, b_bish, ntgrid)
call integrate(dc_bish, outputs%theta, c_bish, ntgrid)
a_b = a_bish(nth) - a_bish(-nth)
b_b = b_bish(nth) - b_bish(-nth)
c_b = c_bish(nth) - c_bish(-nth)
a_bish = - a_bish * bpol * rmajor
b_bish = - b_bish * bpol * rmajor
c_bish = - c_bish * bpol * rmajor
! find shat, pressure' and report them
dp = geom%dpfun(pbar) ; di = geom%dbtori(pbar)
if (abs(qval) > epsilon(0.)) then
tmp = rhoc / (qval * twopi * drhodpsi)
s_hat = tmp * (a_b * di + b_b * dp + c_b * 2) !< Calculate the shear from equib.
if (local_eq) s_hat = surf%shat !Can't determine this by calculation in local_eq
if (debug) write(6,*) "a_b (f'), b_b (p'), c_b terms =", a_b*di, b_b*dp, c_b*2
if (debug) write(6,*) 'qval, rho, drhodpsi, s_hat =', qval, rhoc, drhodpsi, s_hat
else
tmp = 1. ; s_hat = 0.
end if
write(fort11_unit,*)
write(fort11_unit,*) 'Quantity, equilibrium value, value used'
! Decide how we want to set s_hat_new in the default case
s_hat_new = s_hat
select case (bishop)
case (2, 3, 4, 5, 8, 9) !1, 6, 7 keep the initial value
s_hat_new = surf%shat
end select
! Decide how we want to set dp_new
dp_new = dp
select case (bishop)
case (2)
dp_new = p_prime_input
case (3)
pressure = geom%pfun(pbar)
dp_new = -invLp_input * pressure * drhodpsi
case (4, 6, 9)
dp_new = 0.5 * beta_prime_input * drhodpsi
case (5)
dp_new = -alpha_input * drhodpsi / (qval**2 * rmaj)
case (7)
dp_new = dp * dp_mult
case (8)
dp_new = 0.5 * beta_prime_input * drhodpsi * dp_mult
end select
!Report on updated dp
select case(bishop)
case (2)
write(fort11_unit,*) 'p_prime = ',dp,', ',dp_new
if (writelots) write(*,*) 'p_prime = ',dp,', ',dp_new
case (3)
write(fort11_unit,*) '1/L_p = ',-dp / (drhodpsi * pressure), ', ', &
-dp_new / (drhodpsi *pressure)
if (writelots) write(*,*) '1/L_p = ',-dp / (drhodpsi * pressure), ', ', &
-dp_new / (drhodpsi * pressure)
case (5)
write(fort11_unit,*) 'alpha = ', -qval**2 * rmaj * dp / drhodpsi,', ', &
-qval**2 * rmaj * dp_new / drhodpsi
if (writelots) write(*,*) 'alpha = ', -qval**2 * rmaj * dp / drhodpsi,', ', &
-qval**2 * rmaj * dp_new / drhodpsi
case default
write(fort11_unit,*) 'd beta/d rho = ',2 * dp / drhodpsi,', ',2 * dp_new / drhodpsi
if (writelots) write(*,*) 'd beta/d rho = ',2 * dp / drhodpsi,', ',2 * dp_new / drhodpsi
end select
write(fort11_unit,*) 's_hat ',s_hat,', ',s_hat_new
if (writelots) write(*,*) 's_hat ',s_hat,', ',s_hat_new
write(fort11_unit,*)
if ((.not. local_eq) .and. (bishop == 1 .or. bishop > 9)) then
di_new = di
else
di_new = (s_hat_new / tmp - 2 * c_b - dp_new * b_b) / a_b
end if
! Calculate grad(S)
! here the coordinates are rho and l
call gradl(ltheta, seik, dSdl, twopi * qval, nth)
call periodic_copy(dSdl, 0., nperiod)
grads(:,1) = a_bish * di_new + b_bish * dp_new + c_bish * 2
grads(:,2) = dSdl
! Calculate grad(Z)
! obtain Z(theta) and hence dZ/dl
call Ztgrid(rgrid, outputs%theta, Zoftheta)
call periodic_copy(Zoftheta, 0., nperiod)
call gradl(ltheta, Zoftheta, dZdl, 0., nth)
call periodic_copy(dZdl, 0., nperiod)
! gradztot(:,1) is grad Z . rhohat
gradztot(:,1) = crpgrad(:,2) / outputs%grho
! gradztot(:,2) is grad Z . lhat
gradztot(:,2) = dZdl
! gradztot(:,3) is grad Z . phihat
gradztot(:,3) = 0.0
shat = s_hat_new ! Output
dqdrp = s_hat_new * qval * drhodrp / rhoc ! Output
dbetadrho = 2 * dp_new * dpsidrho ! Output
call bishop_gradB(bmag, bpol, invRpol, th_bish, ltheta, rmajor, dp_new, bi, gradbtot)
if (save_bishop) then
call open_output_file(bishop_unit, '.bishop')
write(bishop_unit, fmt="('#di, di_new = ',F12.7,' ',F12.7)") di, di_new
write(bishop_unit, fmt="('#dp, dp_new = ',F12.7,' ',F12.7)") dp, dp_new
write(bishop_unit, fmt="('#bi =', F12.7)") bi
write(bishop_unit, fmt="('#drdpsi =', F12.7)") drhodpsi
write(bishop_unit, fmt="('#a_b (fp)=', F20.10, ' ,b_b (pp)=', F20.10,' ,c_b =', F20.10,' ,s_hat =', F20.10)") a_b*di, b_b*dp, c_b, s_hat
write(bishop_unit, fmt ="('# theta, a_bish, b_bish, c_bish, rmajor, Z, Bpol, invRpol, th_bish')")
do i = -ntgrid, ntgrid
write(bishop_unit, '( 10(2X, F25.17) )') outputs%theta(i), a_bish(i), b_bish(i), &
c_bish(i), rmajor(i), Zoftheta(i), bpol(i), invRpol(i), th_bish(i)
end do
call close_output_file(bishop_unit)
end if
end subroutine bishop_grads
!> @Note local_eq _must_ be false if we end in this routine. We make use of this.
subroutine bishop_large_aspect_ratio_approximation(rhoc, dr, ntheta)
real, intent(in) :: rhoc, dr
integer, intent(in) :: ntheta
real :: delrp, rp1, rp2, qval1, qval2, bi1, bi2, dqdrho, pbar_1, pbar_2
real, dimension(:), allocatable :: seik1, seik2, dsdrp
integer :: itot, k
! rp derivative
delrp = dr / drhodrp
rp1 = rp - delrp ; rp2 = rp + delrp
pbar_1 = pbarfun(rp1, outputs%theta(0)) ; pbar_2 = pbarfun(rp2, outputs%theta(0))
! Calculate modified btori -- note we use the fact that rgrid == rp in
! seikon to avoid having to form rgrid1/rgrid2 here.
bi1 = geom%btori(pbar_1) ; bi2 = geom%btori(pbar_2)
allocate(seik1(-ntgrid:ntgrid), seik2(-ntgrid:ntgrid), dsdrp(-ntgrid:ntgrid))
call seikon(rp1, outputs%theta, invR, qval1, bi1, nperiod, seik1)
call seikon(rp2, outputs%theta, invR, qval2, bi2, nperiod, seik2)
dqdrho = (qval2 - qval1) / (2 * dr)
shat = dqdrho * rhoc / qval ! Output
dqdrp = dqdrho * drhodrp ! Output
dbetadrho = 2 * geom%dpfun(pbar) / drhodrp ! Output
if (writelots) then
write (fort11_unit,*) 'dqdrho=',dqdrho,' qval1,2= ',qval1, qval2
write (fort11_unit,*) 's_hat= ',shat
write (fort11_unit,*) 'i,rgrid,rgrid1,rgrid2:'
do i = -nth, nth
write (fort11_unit,*) i, rgrid(i), rp1, rp2
end do
end if
do i = -nth, nth
dsdrp(i) = (seik2(i) - seik1(i)) / (2 * delrp)
end do
call periodic_copy(dsdrp, -twopi * dqdrp, nperiod)
! compute gradient of S
! this is really gradient of alpha
! note that coordinates are R and Z (not rho and l, as in bishop case)
do k = -nperiod+1, nperiod-1
do i = -nth, nth
itot = i + k * ntheta
grads(itot, :) = rpgrad(i, :) * dsdrp(itot) + thgrad(i, :) * dsdthet(i) ! Output
end do
end do
! gradztot(:,1) is grad Z . Rhat ! Output
gradztot(:,1) = 0.
! gradztot(:,2) is grad Z . Zhat ! Output
gradztot(:,2) = 1.
! gradztot(:,3) is grad Z . phihat ! Output
gradztot(:,3) = 0.
call geom%gradient(rgrid, outputs%theta, gradbtot(:, 1:2), 'B', rp, nth, ntgrid, bishop /= 0)
gradbtot(:, 3) = 0.0
end subroutine bishop_large_aspect_ratio_approximation
!> FIXME : Add documentation
subroutine alloc_local_arrays(n)
implicit none
integer, intent(in) :: n
allocate(seik (-n:n), &
dsdthet (-n:n), &
invR (-n:n), &
trip (-n:n), &
rgrid (-n:n), &
rmajor (-n:n))
allocate(thgrad (-n:n,2), &
rpgrad (-n:n,2), &
crpgrad (-n:n,2), &
grads (-n:n,2))
allocate(gradrptot (-n:n,3), &
gradstot (-n:n,3), &
gradztot (-n:n,3), &
gradthtot (-n:n,3), & ! MAB
gradbtot (-n:n,3))
end subroutine alloc_local_arrays
end subroutine eikcoefs