diff --git a/src/equilibrium.f90 b/src/equilibrium.f90 deleted file mode 100644 index 646cc09..0000000 --- a/src/equilibrium.f90 +++ /dev/null @@ -1,1343 +0,0 @@ -! This module handles the loading, interpolation and evaluation of the -! MHD equilibrium data (poloidal current function, safety factor, -! poloidal flux, magnetic field, etc.) -! -! Two kinds of data are supported: analytical (suffix `_an` in the -! subroutine names) or numerical (suffix `_spline`). For the latter, the -! the data is interpolated using splines. -module equilibrium - use const_and_precisions, only : wp_ - use types, only : contour - use splines, only : spline_simple, spline_1d, spline_2d, linear_1d - use gray_params, only : EQ_VACUUM, EQ_ANALYTICAL, & - EQ_EQDSK_FULL, EQ_EQDSK_PARTIAL - - implicit none - - ! Parameters of the analytical equilibrium model - type analytic_model - real(wp_) :: q0 ! Safety factor at the magnetic axis - real(wp_) :: q1 ! Safety factor at the edge - real(wp_) :: alpha ! Exponent for the q(ρ_p) power law - real(wp_) :: R0 ! R of the magnetic axis (m) - real(wp_) :: z0 ! z of the magnetic axis (m) - real(wp_) :: a ! Minor radius (m) - real(wp_) :: B0 ! Magnetic field at the magnetic axis (T) - end type - - ! Order of the splines - integer, parameter :: kspl=3, ksplp=kspl + 1 - - ! Global variable storing the state of the module - - ! Splines - type(spline_1d), save :: fpol_spline ! Poloidal current function F(ψ) - type(spline_2d), save :: psi_spline ! Normalised poloidal flux ψ(R,z) - type(spline_simple), save :: q_spline ! Safey factor q(ψ) - type(linear_1d), save :: rhop_spline, rhot_spline ! Normalised radii ρ_p(ρ_t), ρ_t(ρ_p) - - ! Analytic model - type(analytic_model), save :: model - - ! Contour of the region where the ψ(R,z) spline is monotonically - ! increasing. The mapping from ψ to other plasma parameters is - ! physically meaningful only inside this domain. - type(contour), save :: psi_domain - - ! More parameters - real(wp_), save :: fpolas ! Poloidal current at the edge (r=a), F_a - real(wp_), save :: psia ! Poloidal flux at the edge (r=a), ψ_a - real(wp_), save :: phitedge ! Toroidal flux at the edge - real(wp_), save :: btaxis ! B₀: B_φ = B₀R₀/R - real(wp_), save :: rmaxis, zmaxis ! R,z of the magnetic axis (O point) - real(wp_), save :: sgnbphi ! Sign of B_φ (>0 means CCW) - real(wp_), save :: btrcen ! used only for jcd_astra def. - real(wp_), save :: rcen ! Major radius R₀, computed as fpolas/btrcen - real(wp_), save :: rmnm, rmxm ! R interval of the equilibrium grid - real(wp_), save :: zmnm, zmxm ! z interval of the equilibrium grid - real(wp_), save :: zbinf, zbsup ! z interval of the plasma boundary - - integer, save :: iequil - - private - public read_eqdsk, read_equil_an ! Reading data files - public scale_equil, change_cocos ! Transforming data - public fq, bfield, tor_curr ! Accessing local quantities - public pol_flux, pol_curr ! - public frhotor, frhopol ! Converting between poloidal/toroidal flux - public set_equil_spline, set_equil_an ! Initialising internal state - public unset_equil_spline ! Deinitialising internal state - - ! Members exposed for magsurf_data - public kspl, psi_spline, find_htg_point - public model, btaxis, btrcen - - ! Members exposed to gray_core and more - public psia, phitedge - public rmaxis, zmaxis, sgnbphi - public rmnm, rmxm, zmnm, zmxm - public zbinf, zbsup - -contains - - subroutine read_eqdsk(params, data, err, unit) - ! Reads the MHD equilibrium `data` from a G-EQDSK file (params%filenm). - ! If given, the file is opened in the `unit` number. - ! For a description of the G-EQDSK, see the GRAY user manual. - use const_and_precisions, only : one - use gray_params, only : equilibrium_parameters, equilibrium_data - use utils, only : get_free_unit - use logger, only : log_error - - ! subroutine arguments - type(equilibrium_parameters), intent(in) :: params - type(equilibrium_data), intent(out) :: data - integer, intent(out) :: err - integer, optional, intent(in) :: unit - - ! local variables - integer :: u, idum, i, j, nr, nz, nbnd, nlim - character(len=48) :: string - real(wp_) :: dr, dz, dps, rleft, zmid, zleft, psiedge, psiaxis - real(wp_) :: xdum ! dummy variable, used to discard data - real(wp_), allocatable :: R(:), z(:) - - u = get_free_unit(unit) - - ! Open the G-EQDSK file - open(file=params%filenm, status='old', action='read', unit=u, iostat=err) - if (err /= 0) then - call log_error('opening eqdsk file ('//trim(params%filenm)//') failed!', & - mod='equilibrium', proc='read_eqdsk') - err = 1 - return - end if - - ! get size of main arrays and allocate them - if (params%idesc) then - read (u,'(a48,3i4)') string,idum,nr,nz - else - read (u,*) nr, nz - end if - if (allocated(data%rv)) deallocate(data%rv) - if (allocated(data%zv)) deallocate(data%zv) - if (allocated(data%psin)) deallocate(data%psin) - if (allocated(data%psinr)) deallocate(data%psinr) - if (allocated(data%fpol)) deallocate(data%fpol) - if (allocated(data%qpsi)) deallocate(data%qpsi) - allocate(data%rv(nr), data%zv(nz), & - data%psin(nr, nz), & - data%psinr(nr), & - data%fpol(nr), & - data%qpsi(nr)) - - ! Store 0D data and main arrays - if (params%ifreefmt) then - read (u, *) dr, dz, data%rvac, rleft, zmid - read (u, *) data%rax, data%zax, psiaxis, psiedge, xdum - read (u, *) xdum, xdum, xdum, xdum, xdum - read (u, *) xdum, xdum, xdum, xdum, xdum - read (u, *) (data%fpol(i), i=1,nr) - read (u, *) (xdum,i=1, nr) - read (u, *) (xdum,i=1, nr) - read (u, *) (xdum,i=1, nr) - read (u, *) ((data%psin(i,j), i=1,nr), j=1,nz) - read (u, *) (data%qpsi(i), i=1,nr) - else - read (u, '(5e16.9)') dr,dz,data%rvac,rleft,zmid - read (u, '(5e16.9)') data%rax,data%zax,psiaxis,psiedge,xdum - read (u, '(5e16.9)') xdum,xdum,xdum,xdum,xdum - read (u, '(5e16.9)') xdum,xdum,xdum,xdum,xdum - read (u, '(5e16.9)') (data%fpol(i),i=1,nr) - read (u, '(5e16.9)') (xdum,i=1,nr) - read (u, '(5e16.9)') (xdum,i=1,nr) - read (u, '(5e16.9)') (xdum,i=1,nr) - read (u, '(5e16.9)') ((data%psin(i,j),i=1,nr),j=1,nz) - read (u, '(5e16.9)') (data%qpsi(i),i=1,nr) - end if - - ! Get size of boundary and limiter arrays and allocate them - read (u, *) nbnd, nlim - - ! Load plasma boundary data - if(nbnd > 0) then - allocate(R(nbnd), z(nbnd)) - if (params%ifreefmt) then - read(u, *) (R(i), z(i), i=1,nbnd) - else - read(u, '(5e16.9)') (R(i), z(i), i=1,nbnd) - end if - data%boundary = contour(R, z) - deallocate(R, z) - end if - - ! Load limiter data - if(nlim > 0) then - allocate(R(nlim), z(nlim)) - if (params%ifreefmt) then - read(u, *) (R(i), z(i), i=1,nlim) - else - read(u, '(5e16.9)') (R(i), z(i), i=1,nlim) - end if - data%limiter = contour(R, z) - deallocate(R, z) - end if - - ! End of G-EQDSK file - close(u) - - ! Build rv,zv,psinr arrays - zleft = zmid-0.5_wp_*dz - dr = dr/(nr-1) - dz = dz/(nz-1) - dps = one/(nr-1) - do i=1,nr - data%psinr(i) = (i-1)*dps - data%rv(i) = rleft + (i-1)*dr - end do - do i=1,nz - data%zv(i) = zleft + (i-1)*dz - end do - - ! Normalize psin - data%psia = psiedge - psiaxis - if(.not. params%ipsinorm) data%psin = (data%psin - psiaxis)/data%psia - - end subroutine read_eqdsk - - - subroutine read_equil_an(filenm, ipass, data, err, unit) - ! Reads the MHD equilibrium `data` in the analytical format - ! from params%filenm. - ! If given, the file is opened in the `unit` number. - ! - ! TODO: add format description - use gray_params, only : equilibrium_data - use utils, only : get_free_unit - use logger, only : log_error - - ! subroutine arguments - character(len=*), intent(in) :: filenm - integer, intent(in) :: ipass - type(equilibrium_data), intent(out) :: data - integer, intent(out) :: err - integer, optional, intent(in) :: unit - - ! local variables - integer :: i, u, nlim - real(wp_), allocatable :: R(:), z(:) - - u = get_free_unit(unit) - - open(file=filenm, status='old', action='read', unit=u, iostat=err) - if (err /= 0) then - call log_error('opening equilibrium file ('//trim(filenm)//') failed!', & - mod='equilibrium', proc='read_equil_an') - err = 1 - return - end if - - read(u, *) model%R0, model%z0, model%a - read(u, *) model%B0 - read(u, *) model%q0, model%q1, model%alpha - - if(ipass >= 2) then - ! get size of boundary and limiter arrays and allocate them - read (u,*) nlim - - ! store boundary and limiter data - if (nlim > 0) then - allocate(R(nlim), z(nlim)) - read(u,*) (R(i), z(i), i = 1, nlim) - data%limiter = contour(R, z) - end if - end if - close(u) - end subroutine read_equil_an - - - subroutine change_cocos(data, cocosin, cocosout, error) - ! Convert the MHD equilibrium data from one coordinate convention - ! (COCOS) to another. These are specified by `cocosin` and - ! `cocosout`, respectively. - ! - ! For more information, see: https://doi.org/10.1016/j.cpc.2012.09.010 - use const_and_precisions, only : zero, one, pi - use gray_params, only : equilibrium_data - - ! subroutine arguments - type(equilibrium_data), intent(inout) :: data - integer, intent(in) :: cocosin, cocosout - integer, intent(out), optional :: error - - ! local variables - real(wp_) :: isign, bsign - integer :: exp2pi, exp2piout - logical :: phiccw, psiincr, qpos, phiccwout, psiincrout, qposout - - call decode_cocos(cocosin, exp2pi, phiccw, psiincr, qpos) - call decode_cocos(cocosout, exp2piout, phiccwout, psiincrout, qposout) - - ! Check sign consistency - isign = sign(one, data%psia) - if (.not.psiincr) isign = -isign - bsign = sign(one, data%fpol(size(data%fpol))) - if (qpos .neqv. isign * bsign * data%qpsi(size(data%qpsi)) > zero) then - ! Warning: sign inconsistency found among q, Ipla and Bref - data%qpsi = -data%qpsi - if (present(error)) error = 1 - else - if (present(error)) error = 0 - end if - - ! Convert cocosin to cocosout - - ! Opposite direction of toroidal angle phi in cocosin and cocosout - if (phiccw .neqv. phiccwout) data%fpol = -data%fpol - - ! q has opposite sign for given sign of Bphi*Ip - if (qpos .neqv. qposout) data%qpsi = -data%qpsi - - ! psi and Ip signs don't change accordingly - if ((phiccw .eqv. phiccwout) .neqv. (psiincr .eqv. psiincrout)) & - data%psia = -data%psia - - ! Convert Wb to Wb/rad or viceversa - data%psia = data%psia * (2.0_wp_*pi)**(exp2piout - exp2pi) - end subroutine change_cocos - - - subroutine decode_cocos(cocos, exp2pi, phiccw, psiincr, qpos) - ! Extracts the sign and units conventions from a COCOS index - - ! subroutine arguments - integer, intent(in) :: cocos - integer, intent(out) :: exp2pi - logical, intent(out) :: phiccw, psiincr, qpos - - ! local variables - integer :: cmod10, cmod4 - - cmod10 = mod(cocos, 10) - cmod4 = mod(cmod10, 4) - - ! cocos>10 ψ in Wb, cocos<10 ψ in Wb/rad - exp2pi = (cocos - cmod10)/10 - - ! cocos mod 10 = 1,3,5,7: toroidal angle φ increasing CCW - phiccw = (mod(cmod10, 2)== 1) - - ! cocos mod 10 = 1,2,5,6: ψ increasing with positive Ip - psiincr = (cmod4==1 .or. cmod4==2) - - ! cocos mod 10 = 1,2,7,8: q positive for positive Bφ*Ip - qpos = (cmod10<3 .or. cmod10>6) - end subroutine decode_cocos - - - subroutine scale_equil(params, data) - ! Rescale the magnetic field (B) and the plasma current (I_p) - ! and/or force their signs. - ! - ! Notes: - ! 1. signi and signb are ignored on input if equal to 0. - ! They are used to assign the direction of B_φ and I_p BEFORE scaling. - ! 2. cocos=3 assumed: positive toroidal direction is CCW from above - ! 3. B_φ and I_p scaled by the same factor factb to keep q unchanged - ! 4. factb<0 reverses the directions of Bphi and Ipla - - use const_and_precisions, only : one - use gray_params, only : equilibrium_parameters, equilibrium_data - - ! subroutine arguments - type(equilibrium_parameters), intent(inout) :: params - type(equilibrium_data), intent(inout) :: data - - ! Notes on cocos=3 - ! 1. With both I_p,B_φ > 0 we have ∂ψ/∂r<0 and ∂Φ/∂r>0. - ! 2. ψ_a = ψ_edge - ψ_axis < 0. - ! 3. q = 1/2π ∂Φ/∂ψ ~ ∂Φ/∂r⋅∂r/∂ψ < 0. - ! 4. In general, sgn(q) = -sgn(I_p)⋅sgn(B_φ). - - iequil = params%iequil - - select case(iequil) - case (EQ_ANALYTICAL) ! Analytical model - ! Apply signs - if (params%sgni /= 0) then - model%q0 = sign(model%q0, -params%sgni*one) - model%q1 = sign(model%q1, -params%sgni*one) - end if - if (params%sgnb /= 0) then - model%B0 = sign(model%B0, +params%sgnb*one) - end if - ! Rescale - model%B0 = model%B0 * params%factb - - case (EQ_EQDSK_FULL, EQ_EQDSK_PARTIAL) ! Numeric data - ! Apply signs - if (params%sgni /= 0) & - data%psia = sign(data%psia, -params%sgni*one) - if (params%sgnb /= 0) & - data%fpol = sign(data%fpol, +params%sgnb*one) - ! Rescale - data%psia = data%psia * params%factb - data%fpol = data%fpol * params%factb - - ! Compute the signs to be shown in the outputs header - ! Note: this is needed if sgni, sgnb = 0 in gray.ini - params%sgni = int(sign(one, -data%psia)) - params%sgnb = int(sign(one, +data%fpol(size(data%fpol)))) - end select - end subroutine scale_equil - - - subroutine set_equil_spline(params, data, err) - ! Computes splines for the MHD equilibrium data and stores them - ! in their respective global variables, see the top of this file. - use const_and_precisions, only : zero, one - use gray_params, only : equilibrium_parameters, equilibrium_data - use gray_params, only : X_AT_TOP, X_AT_BOTTOM, X_IS_MISSING - use utils, only : vmaxmin, vmaxmini, inside - use logger, only : log_info - - ! subroutine arguments - type(equilibrium_parameters), intent(in) :: params - type(equilibrium_data), intent(in) :: data - integer, intent(out) :: err - - ! local variables - integer :: nr, nz, nrest, nzest, npsest, nrz, npsi, nbnd, ibinf, ibsup - real(wp_) :: psinoptmp, psinxptmp - real(wp_) :: rbinf, rbsup, r1, z1 - real(wp_) :: psiant, psinop - real(wp_) :: rhotn(size(data%psinr)) - real(wp_), allocatable :: rv1d(:), zv1d(:), fvpsi(:) - integer :: i, j, ij - character(256) :: msg ! for log messages formatting - - ! compute array sizes - nr = size(data%rv) - nz = size(data%zv) - nrz = nr*nz - npsi = size(data%psinr) - nrest = nr + ksplp - nzest = nz + ksplp - npsest = npsi + ksplp - - ! length in m !!! - rmnm = data%rv(1) - rmxm = data%rv(nr) - zmnm = data%zv(1) - zmxm = data%zv(nz) - - ! Spline interpolation of ψ(R, z) - - select case (iequil) - - case (EQ_EQDSK_PARTIAL) - ! Data valid only inside boundary (data%psin=0 outside), - ! presence of boundary anticipated here to filter invalid data - nbnd = size(data%boundary%R) - - ! allocate knots and spline coefficients arrays - if (allocated(psi_spline%knots_x)) deallocate(psi_spline%knots_x) - if (allocated(psi_spline%knots_y)) deallocate(psi_spline%knots_y) - if (allocated(psi_spline%coeffs)) deallocate(psi_spline%coeffs) - allocate(psi_spline%knots_x(nrest), psi_spline%knots_y(nzest)) - allocate(psi_spline%coeffs(nrz)) - - ! determine number of valid grid points - nrz=0 - do j=1,nz - do i=1,nr - if (nbnd.gt.0) then - if (.not. data%boundary%contains(data%rv(i), data%zv(j))) cycle - else - if(data%psin(i,j).le.0.0d0) cycle - end if - nrz=nrz+1 - end do - end do - - ! store valid data - allocate(rv1d(nrz),zv1d(nrz),fvpsi(nrz)) - ij=0 - do j=1,nz - do i=1,nr - if (nbnd.gt.0) then - if (.not. data%boundary%contains(data%rv(i), data%zv(j))) cycle - else - if(data%psin(i,j).le.0.0d0) cycle - end if - ij=ij+1 - rv1d(ij)=data%rv(i) - zv1d(ij)=data%zv(j) - fvpsi(ij)=data%psin(i,j) - end do - end do - - ! Fit as a scattered set of points - ! use reduced number of knots to limit memory comsumption ? - psi_spline%nknots_x=nr/4+4 - psi_spline%nknots_y=nz/4+4 - call psi_spline%init_nonreg(rv1d, zv1d, fvpsi, tension=params%ssplps, & - range=[rmnm, rmxm, zmnm, zmxm], err=err) - - ! if failed, re-fit with an interpolating spline (zero tension) - if(err == -1) then - err = 0 - psi_spline%nknots_x=nr/4+4 - psi_spline%nknots_y=nz/4+4 - call psi_spline%init_nonreg(rv1d, zv1d, fvpsi, tension=zero, & - range=[rmnm, rmxm, zmnm, zmxm], err=err) - end if - deallocate(rv1d, zv1d, fvpsi) - ! reset nrz to the total number of grid points for next allocations - nrz = nr*nz - - case (EQ_EQDSK_FULL) - ! Data are valid on the full R,z grid - - ! reshape 2D ψ array to 1D (transposed) - allocate(fvpsi(nrz)) - fvpsi = reshape(transpose(data%psin), [nrz]) - - ! compute spline coefficients - call psi_spline%init(data%rv, data%zv, fvpsi, nr, nz, & - range=[rmnm, rmxm, zmnm, zmxm], & - tension=params%ssplps, err=err) - - ! if failed, re-fit with an interpolating spline (zero tension) - if(err == -1) then - call psi_spline%init(data%rv, data%zv, fvpsi, nr, nz, & - range=[rmnm, rmxm, zmnm, zmxm], & - tension=zero) - err = 0 - end if - deallocate(fvpsi) - end select - - if (err /= 0) then - err = 2 - return - end if - - ! compute spline coefficients for ψ(R,z) partial derivatives - call psi_spline%init_deriv(nr, nz, 1, 0) ! ∂ψ/∂R - call psi_spline%init_deriv(nr, nz, 0, 1) ! ∂ψ/∂z - call psi_spline%init_deriv(nr, nz, 1, 1) ! ∂²ψ/∂R∂z - call psi_spline%init_deriv(nr, nz, 2, 0) ! ∂²ψ/∂R² - call psi_spline%init_deriv(nr, nz, 0, 2) ! ∂²ψ/∂z² - - ! set the initial ψ(R,z) domain to the grid boundary - ! - ! Note: this is required to bootstrap `flux_pol` calls - ! within this very subroutine. - psi_domain%R = [rmnm, rmnm, rmxm, rmxm] - psi_domain%z = [zmnm, zmxm, zmxm, zmnm] - - ! Spline interpolation of F(ψ) - - ! give a small weight to the last point - block - real(wp_) :: w(npsi) - w(1:npsi-1) = 1 - w(npsi) = 1.0e2_wp_ - call fpol_spline%init(data%psinr, data%fpol, npsi, range=[zero, one], & - weights=w, tension=params%ssplf) - end block - - ! set vacuum value used outside 0 ≤ ψ ≤ 1 range - fpolas = fpol_spline%eval(data%psinr(npsi)) - sgnbphi = sign(one ,fpolas) - - ! Re-normalize ψ_n after the spline computation - ! Note: this ensures 0 ≤ ψ_n'(R,z) < 1 inside the plasma - - ! Start with un-corrected ψ_n - psia = data%psia - psinop = 0 ! ψ_n(O point) - psiant = 1 ! ψ_n(X point) - ψ_n(O point) - - ! Use provided boundary to set an initial guess - ! for the search of O/X points - nbnd = size(data%boundary%R) - if (nbnd > 0) then - call vmaxmini(data%boundary%z, nbnd, zbinf, zbsup, ibinf, ibsup) - rbinf = data%boundary%R(ibinf) - rbsup = data%boundary%R(ibsup) - else - zbinf=data%zv(2) - zbsup=data%zv(nz-1) - rbinf=data%rv((nr+1)/2) - rbsup=rbinf - end if - - ! Search for exact location of the magnetic axis - call find_ox_point(data%rax,data%zax,rmaxis,zmaxis,psinoptmp) - - write (msg, '("O-point found:", 3(x,a,"=",g0.3))') & - 'r', rmaxis, 'z', zmaxis, 'ψ', psinoptmp - call log_info(msg, mod='equilibrium', proc='set_equil_spline') - - ! Search for X-point - select case (params%ixp) - case (X_AT_BOTTOM) - call find_ox_point(rbinf,zbinf,r1,z1,psinxptmp) - if(psinxptmp/=-1.0_wp_) then - write (msg, '("X-point found:", 3(x,a,"=",g0.3))') & - 'r', r1, 'z', z1, 'ψ', psinxptmp - call log_info(msg, mod='equilibrium', proc='set_equil_spline') - - zbinf=z1 - psinop=psinoptmp - psiant=psinxptmp-psinop - call find_htg_point(rmaxis,0.5_wp_*(zmaxis+zbsup),r1,z1,one) - zbsup=z1 - end if - - case (X_AT_TOP) - call find_ox_point(rbsup,zbsup,r1,z1,psinxptmp) - if(psinxptmp.ne.-1.0_wp_) then - write (msg, '("X-point found:", 3(x,a,"=",g0.3))') & - 'r', r1, 'z', z1, 'ψ', psinxptmp - call log_info(msg, mod='equilibrium', proc='set_equil_spline') - - zbsup=z1 - psinop=psinoptmp - psiant=psinxptmp-psinop - call find_htg_point(rmaxis,0.5_wp_*(zmaxis+zbinf),r1,z1,one) - zbinf=z1 - end if - - case (X_IS_MISSING) - psinop=psinoptmp - psiant=one-psinop - ! Find upper horizontal tangent point - call find_htg_point(rmaxis,0.5_wp_*(zmaxis+zbsup),r1,z1,one) - zbsup=z1 - rbsup=r1 - ! Find lower horizontal tangent point - call find_htg_point(rmaxis,0.5_wp_*(zmaxis+zbinf),r1,z1,one) - zbinf=z1 - rbinf=r1 - write (msg, '("X-point not found in", 2(x,a,"∈[",g0.3,",",g0.3,"]"))') & - 'r', rbinf, rbsup, 'z', zbinf, zbsup - call log_info(msg, mod='equilibrium', proc='set_equil_spline') - end select - - ! Correct the spline coefficients: ψ_n → (ψ_n - psinop)/psiant - call psi_spline%transform(1/psiant, -psinop/psiant) - - ! Do the opposite scaling to preserve un-normalised values - ! Note: this is only used for the poloidal magnetic field - psia = psia * psiant - - ! Save Bt value on axis (required in flux_average and used in Jcd def) - ! and vacuum value B0 at ref. radius data%rvac (used in Jcd_astra def) - - call pol_curr(zero, btaxis) - btaxis = btaxis/rmaxis - btrcen = fpolas/data%rvac - rcen = data%rvac - - write (msg, '(2(a,g0.3))') 'Bt_center=', btrcen, ' Bt_axis=', btaxis - call log_info(msg, mod='equilibrium', proc='set_equil_spline') - - ! Compute ρ_p/ρ_t mapping based on the input q profile - call setqphi_num(data%psinr, abs(data%qpsi), abs(psia), rhotn) - call set_rho_spline(sqrt(data%psinr), rhotn) - - ! Compute the domain of the ψ mapping - psi_domain = data%boundary - call rescale_boundary(psi_domain, O=[rmaxis, zmaxis], t0=1.5_wp_) - - end subroutine set_equil_spline - - - subroutine rescale_boundary(cont, O, t0) - ! Given the plasma boundary contour `cont`, the position of the - ! magnetic axis `O`, and a scaling factor `t0`; this subroutine - ! rescales the contour by `t0` about `O` while ensuring the - ! psi_spline stays monotonic within the new boundary. - - ! subroutine arguments - type(contour), intent(inout) :: cont ! (R,z) contour - real(wp_), intent(in) :: O(2) ! center point - real(wp_), intent(in) :: t0 ! scaling factor - - ! subroutine variables - integer :: i - real(wp_) :: t - real(wp_), parameter :: dt = 0.05 - real(wp_) :: P(2), N(2) - - do i = 1, size(cont%R) - ! For each point on the contour compute: - P = [cont%R(i), cont%z(i)] ! point on the contour - N = P - O ! direction of the line from O to P - - ! Find the max t: s(t) = ψ(O + tN) is monotonic in [1, t] - t = 1 - do while (t < t0) - if (s(t + dt) < s(t)) exit - t = t + dt - end do - - ! The final point is thus O + tN - P = O + t * N - cont%R(i) = P(1) - cont%z(i) = P(2) - end do - - contains - - function s(t) - ! Rescriction of ψ(R, z) on the line Q(t) = O + tN - - real(wp_), intent(in) :: t - real(wp_) :: s, Q(2) - - Q = O + t * N - s = psi_spline%eval(Q(1), Q(2)) - end function - - end subroutine rescale_boundary - - - subroutine setqphi_num(psinq,q,psia,rhotn) - ! Computes the spline of the safety factor q(ψ) - use const_and_precisions, only : pi - - ! subroutine arguments - real(wp_), dimension(:), intent(in) :: psinq,q - real(wp_), intent(in) :: psia - real(wp_), dimension(:), intent(out), optional :: rhotn - - ! local variables - real(wp_), dimension(size(q)) :: phit - real(wp_) :: dx - integer :: k - - call q_spline%init(psinq, q, size(q)) - - ! Toroidal flux as Φ(ψ) = 2π ∫q(ψ)dψ - phit(1)=0 - do k=1,q_spline%ndata-1 - dx=q_spline%data(k+1)-q_spline%data(k) - phit(k+1)=phit(k) + dx*(q_spline%coeffs(k,1) + dx*(q_spline%coeffs(k,2)/2 + & - dx*(q_spline%coeffs(k,3)/3 + dx* q_spline%coeffs(k,4)/4) ) ) - end do - phitedge=phit(q_spline%ndata) - if(present(rhotn)) rhotn(1:q_spline%ndata)=sqrt(phit/phitedge) - phitedge=2*pi*psia*phitedge - end subroutine setqphi_num - - - subroutine set_equil_an - ! Computes the analytical equilibrium data and stores them - ! in their respective global variables, see the top of this file. - use const_and_precisions, only : pi, one - - real(wp_) :: dq, gamma - - btaxis = model%B0 - rmaxis = model%R0 - zmaxis = model%z0 - btrcen = model%B0 - rcen = model%R0 - zbinf = model%z0 - model%a - zbsup = model%z0 + model%a - sgnbphi = sign(one, model%B0) - - rmxm = model%R0 + model%a - rmnm = model%R0 - model%a - zmxm = zbsup - zmnm = zbinf - - ! Toroidal flux at r=a: - ! - ! Φ(a) = B₀πa² 2γ/(γ + 1) - ! - ! where γ=1/√(1-ε²), - ! ε=a/R₀ is the tokamak aspect ratio - gamma = 1/sqrt(1 - (model%a/model%R0)**2) - phitedge = model%B0 * pi * model%a**2 * 2*gamma/(gamma + 1) - - ! In cocos=3 the safety factor is - ! - ! q(ψ) = 1/2π ∂Φ/∂ψ. - ! - ! Given the power law of the model - ! - ! q(ψ) = q₀ + (q₁-q₀) (ψ/ψa)^(α/2), - ! - ! we can find ψ_a = ψ(r=a) by integrating: - ! - ! ∫ q(ψ)dψ = 1/2π ∫ dΦ - ! ∫₀^ψ_a q(ψ)dψ = 1/2π Φ_a - ! ψa [q₀ + (q₁-q₀)/(α/2+1)] = Φa/2π - ! - ! ⇒ ψ_a = Φ_a 1/2π 1/(q₀ + Δq) - ! - ! where Δq = (q₁ - q₀)/(α/2 + 1) - dq = (model%q1 - model%q0) / (model%alpha/2 + 1) - psia = 1/(2*pi) * phitedge / (model%q0 + dq) - end subroutine set_equil_an - - - subroutine set_rho_spline(rhop, rhot) - ! Computes the splines for converting between the poloidal (ρ_p) - ! and toroidal (ρ_t) normalised radii - - ! subroutine arguments - real(wp_), dimension(:), intent(in) :: rhop, rhot - - call rhop_spline%init(rhot, rhop, size(rhop)) - call rhot_spline%init(rhop, rhot, size(rhot)) - end subroutine set_rho_spline - - - subroutine pol_flux(R, z, psi_n, dpsidr, dpsidz, & - ddpsidrr, ddpsidzz, ddpsidrz) - ! Computes the normalised poloidal flux ψ_n and its - ! derivatives wrt (R, z) up to the second order. - ! - ! Note: all output arguments are optional. - use utils, only : inside - use const_and_precisions, only : pi - - ! subroutine arguments - real(wp_), intent(in) :: R, z - real(wp_), intent(out), optional :: psi_n, dpsidr, dpsidz, & - ddpsidrr, ddpsidzz, ddpsidrz - - ! local variables - real(wp_) :: r_g, rho_t, rho_p ! geometric radius, √Φ_n, √ψ_n - real(wp_) :: gamma ! γ = 1/√(1 - r²/R₀²) - real(wp_) :: dpsidphi ! (∂ψ_n/∂Φ_n) - real(wp_) :: ddpsidphidr, ddpsidphidz ! ∇(∂ψ_n/∂Φ_n) - real(wp_) :: phi_n ! Φ_n - real(wp_) :: dphidr, dphidz ! ∇Φ_n - real(wp_) :: ddphidrdr, ddphidzdz ! ∇∇Φ_n - real(wp_) :: ddphidrdz ! - real(wp_) :: q, dq ! q(ρ_p), Δq=(q₁-q₀)/(α/2 + 1) - real(wp_) :: dqdr, dqdz ! ∇q - real(wp_) :: dphidr2, ddphidr2dr2 ! dΦ_n/d(r²), d²Φ_n/d(r²)² - - ! Values for vacuum/outside the domain - if (present(psi_n)) psi_n = -1 - if (present(dpsidr)) dpsidr = 0 - if (present(dpsidz)) dpsidz = 0 - if (present(ddpsidrr)) ddpsidrr = 0 - if (present(ddpsidzz)) ddpsidzz = 0 - if (present(ddpsidrz)) ddpsidrz = 0 - - select case (iequil) - case (EQ_ANALYTICAL) - ! Analytical model - ! - ! The normalised poloidal flux ψ_n(R, z) is computed as follows: - ! 1. ψ_n = ρ_p² - ! 2. ρ_p = ρ_p(ρ_t), using `frhopol`, which in turns uses q(ψ) - ! 3. ρ_t = √Φ_n - ! 4. Φ_n = Φ(r)/Φ(a), where Φ(r) is the flux of B_φ=B₀R₀/R - ! through a circular surface - ! 5. r = √[(R-R₀)²+(z-z₀)²] is the geometric minor radius - r_g = hypot(R - model%R0, z - model%z0) - - ! The exact flux of the toroidal field B_φ = B₀R₀/R is: - ! - ! Φ(r) = B₀πr² 2γ/(γ + 1) where γ=1/√(1 - r²/R₀²). - ! - ! Notes: - ! 1. the function Φ(r) is defined for r≤R₀ only. - ! 2. as r → 0, γ → 1, so Φ ~ B₀πr². - ! 3. as r → 1⁻, Φ → 2B₀πr² but dΦ/dr → -∞. - ! 4. |B_R|, |B_z| → +-∞. - ! - if (r_g > model%R0) then - if (present(psi_n)) psi_n = -1 - if (present(dpsidr)) dpsidr = 0 - if (present(dpsidz)) dpsidz = 0 - if (present(ddpsidrr)) ddpsidrr = 0 - if (present(ddpsidzz)) ddpsidzz = 0 - if (present(ddpsidrz)) ddpsidrz = 0 - return - end if - - gamma = 1 / sqrt(1 - (r_g/model%R0)**2) - phi_n = model%B0 * pi*r_g**2 * 2*gamma/(gamma + 1) / phitedge - rho_t = sqrt(phi_n) - rho_p = frhopol(rho_t) - - ! For ∇Φ_n and ∇∇Φ_n we also need: - ! - ! ∂Φ∂(r²) = B₀π γ(r) - ! ∂²Φ∂(r²)² = B₀π γ³(r) / (2 R₀²) - ! - dphidr2 = model%B0 * pi * gamma / phitedge - ddphidr2dr2 = model%B0 * pi * gamma**3/(2 * model%R0**2) / phitedge - - ! ∇Φ_n = ∂Φ_n/∂(r²) ∇(r²) - ! where ∇(r²) = 2[(R-R₀), (z-z₀)] - dphidr = dphidr2 * 2*(R - model%R0) - dphidz = dphidr2 * 2*(z - model%z0) - - ! ∇∇Φ_n = ∇[∂Φ_n/∂(r²)] ∇(r²) + ∂Φ_n/∂(r²) ∇∇(r²) - ! = ∂²Φ_n/∂(r²)² ∇(r²)∇(r²) + ∂Φ_n/∂(r²) ∇∇(r²) - ! where ∇∇(r²) = 2I - ddphidrdr = ddphidr2dr2 * 4*(R - model%R0)*(R - model%R0) + dphidr2*2 - ddphidzdz = ddphidr2dr2 * 4*(z - model%z0)*(z - model%z0) + dphidr2*2 - ddphidrdz = ddphidr2dr2 * 4*(R - model%R0)*(z - model%z0) - - ! ψ_n = ρ_p(ρ_t)² - if (present(psi_n)) psi_n = rho_p**2 - - ! Using the definitions in `frhotor`: - ! - ! ∇ψ_n = ∂ψ_n/∂Φ_n ∇Φ_n - ! - ! ∂ψ_n/∂Φ_n = Φ_a/ψ_a ∂ψ/∂Φ - ! = Φ_a/ψ_a 1/2πq - ! - ! Using ψ_a = 1/2π Φ_a / (q₀ + Δq), then: - ! - ! ∂ψ_n/∂Φ_n = (q₀ + Δq)/q - ! - q = model%q0 + (model%q1 - model%q0) * rho_p**model%alpha - dq = (model%q1 - model%q0) / (model%alpha/2 + 1) - dpsidphi = (model%q0 + dq) / q - - ! Using the above, ∇ψ_n = ∂ψ_n/∂Φ_n ∇Φ_n - if (present(dpsidr)) dpsidr = dpsidphi * dphidr - if (present(dpsidz)) dpsidz = dpsidphi * dphidz - - ! For the second derivatives: - ! - ! ∇∇ψ_n = ∇(∂ψ_n/∂Φ_n) ∇Φ_n + (∂ψ_n/∂Φ_n) ∇∇Φ_n - ! - ! ∇(∂ψ_n/∂Φ_n) = - (∂ψ_n/∂Φ_n) ∇q/q - ! - ! From q(ψ) = q₀ + (q₁-q₀) ψ_n^α/2, we have: - ! - ! ∇q = α/2 (q-q₀) ∇ψ_n/ψ_n - ! = α/2 (q-q₀)/ψ_n (∂ψ_n/∂Φ_n) ∇Φ_n. - ! - dqdr = model%alpha/2 * (model%q1 - model%q0)*rho_p**(model%alpha-2) * dpsidphi * dphidr - dqdz = model%alpha/2 * (model%q1 - model%q0)*rho_p**(model%alpha-2) * dpsidphi * dphidz - ddpsidphidr = - dpsidphi * dqdr/q - ddpsidphidz = - dpsidphi * dqdz/q - - ! Combining all of the above: - ! - ! ∇∇ψ_n = ∇(∂ψ_n/∂Φ_n) ∇Φ_n + (∂ψ_n/∂Φ_n) ∇∇Φ_n - ! - if (present(ddpsidrr)) ddpsidrr = ddpsidphidr * dphidr + dpsidphi * ddphidrdr - if (present(ddpsidzz)) ddpsidzz = ddpsidphidz * dphidz + dpsidphi * ddphidzdz - if (present(ddpsidrz)) ddpsidrz = ddpsidphidr * dphidz + dpsidphi * ddphidrdz - - case (EQ_EQDSK_FULL, EQ_EQDSK_PARTIAL) - ! Numerical data - if (psi_domain%contains(R, z)) then - ! Within the interpolation range - if (present(psi_n)) psi_n = psi_spline%eval(R, z) - if (present(dpsidr)) dpsidr = psi_spline%deriv(R, z, 1, 0) - if (present(dpsidz)) dpsidz = psi_spline%deriv(R, z, 0, 1) - if (present(ddpsidrr)) ddpsidrr = psi_spline%deriv(R, z, 2, 0) - if (present(ddpsidzz)) ddpsidzz = psi_spline%deriv(R, z, 0, 2) - if (present(ddpsidrz)) ddpsidrz = psi_spline%deriv(R, z, 1, 1) - end if - - end select - end subroutine pol_flux - - - subroutine pol_curr(psi_n, fpol, dfpol) - ! Computes the poloidal current function F(ψ_n) - ! and (optionally) its derivative dF/dψ_n given ψ_n. - - ! function arguments - real(wp_), intent(in) :: psi_n ! normalised poloidal flux - real(wp_), intent(out) :: fpol ! poloidal current - real(wp_), intent(out), optional :: dfpol ! derivative - - select case (iequil) - case (EQ_VACUUM) - ! Vacuum, no plasma - fpol = 0 - if (present(dfpol)) dfpol = 0 - - case (EQ_ANALYTICAL) - ! Analytical model - ! F(ψ) = B₀⋅R₀, a constant - fpol = model%B0 * model%R0 - if (present(dfpol)) dfpol = 0 - - case (EQ_EQDSK_FULL, EQ_EQDSK_PARTIAL) - ! Numerical data - if(psi_n <= 1 .and. psi_n >= 0) then - fpol = fpol_spline%eval(psi_n) - if (present(dfpol)) dfpol = fpol_spline%deriv(psi_n) - else - fpol = fpolas - if (present(dfpol)) dfpol = 0 - end if - end select - end subroutine pol_curr - - - function frhotor(rho_p) - ! Converts from poloidal (ρ_p) to toroidal (ρ_t) normalised radius - - ! function arguments - real(wp_), intent(in) :: rho_p - real(wp_) :: frhotor - - select case (iequil) - - case (EQ_ANALYTICAL) - ! Analytical model - block - ! The change of variable is obtained by integrating - ! - ! q(ψ) = 1/2π ∂Φ/∂ψ - ! - ! and defining ψ = ψ_a ρ_p², Φ = Φ_a ρ_t². - ! The result is: - ! - ! - ψ_a = 1/2π Φ_a / [q₀ + Δq] - ! - ! - ρ_t = ρ_p √[(q₀ + Δq ρ_p^α)/(q₀ + Δq)] - ! - ! where Δq = (q₁ - q₀)/(α/2 + 1) - real(wp_) :: dq - dq = (model%q1 - model%q0) / (model%alpha/2 + 1) - frhotor = rho_p * sqrt((model%q0 + dq*rho_p**model%alpha) & - / (model%q0 + dq)) - end block - - case (EQ_EQDSK_FULL, EQ_EQDSK_PARTIAL) - ! Numerical data - frhotor = rhot_spline%eval(rho_p) - - end select - end function frhotor - - - function frhopol(rho_t) - ! Converts from toroidal (ρ_t) to poloidal (ρ_p) normalised radius - use const_and_precisions, only : comp_eps - - ! function arguments - real(wp_), intent(in) :: rho_t - real(wp_) :: frhopol - - select case (iequil) - case (EQ_VACUUM) - ! Vacuum, no plasma - frhopol = 0 - - case (EQ_ANALYTICAL) - ! Analytical model - block - ! In general there is no closed form for ρ_p(ρ_t) in the - ! analytical model, we thus solve numerically the equation - ! ρ_t(ρ_p) = ρ_t₀ for ρ_p. - use minpack, only : hybrj1 - - real(wp_) :: rho_p(1), fvec(1), fjac(1,1), wa(7) - integer :: info - - rho_p = [rho_t] ! first guess, ρ_p ≈ ρ_t - call hybrj1(equation, n=1, x=rho_p, fvec=fvec, fjac=fjac, & - ldfjac=1, tol=comp_eps, info=info, wa=wa, lwa=7) - frhopol = rho_p(1) - end block - - case (EQ_EQDSK_FULL, EQ_EQDSK_PARTIAL) - ! Numerical data - frhopol = rhop_spline%eval(rho_t) - end select - - contains - - subroutine equation(n, x, f, df, ldf, flag) - ! The equation to solve: f(x) = ρ_t(x) - ρ_t₀ = 0 - - ! optimal step size - real(wp_), parameter :: e = comp_eps**(1/3.0_wp_) - - ! subroutine arguments - integer, intent(in) :: n, ldf, flag - real(wp_), intent(in) :: x(n) - real(wp_), intent(inout) :: f(n), df(ldf,n) - - if (flag == 1) then - ! return f(x) - f(1) = frhotor(x(1)) - rho_t - else - ! return f'(x), computed numerically - if (x(1) - e > 0) then - df(1,1) = (frhotor(x(1) + e) - frhotor(x(1) - e)) / (2*e) - else - ! single-sided when close to ρ=0 - df(1,1) = (frhotor(x(1) + e) - frhotor(x(1))) / e - end if - end if - end subroutine - - end function frhopol - - - function fq(psin) - ! Computes the safety factor q as a function of the normalised - ! poloidal flux ψ. - ! - ! Note: this returns the absolute value of q. - - ! function arguments - real(wp_), intent(in) :: psin - real(wp_) :: fq - - select case(iequil) - case (EQ_VACUUM) - ! Vacuum, q is undefined - fq = 0 - - case (EQ_ANALYTICAL) - ! Analytical model - ! The safety factor is a power law in ρ_p: - ! q(ρ_p) = q₀ + (q₁-q₀) ρ_p^α - block - real(wp_) :: rho_p - rho_p = sqrt(psin) - fq = abs(model%q0 + (model%q1 - model%q0) * rho_p**model%alpha) - end block - - case (EQ_EQDSK_FULL, EQ_EQDSK_PARTIAL) - ! Numerical data - if (psin < 1) then - fq = q_spline%eval(psin) - else - ! q is undefined or infinite - fq = 0 - end if - - end select - end function fq - - - subroutine bfield(R, z, B_R, B_z, B_phi) - ! Computes the magnetic field as a function of - ! (R, z) in cylindrical coordinates - ! - ! Note: all output arguments are optional. - - ! subroutine arguments - real(wp_), intent(in) :: R, z - real(wp_), intent(out), optional :: B_R, B_z, B_phi - - ! local variables - real(wp_) :: psi_n, fpol, dpsidr, dpsidz - - if (iequil == EQ_VACUUM) then - ! Vacuum, no plasma nor field - if (present(B_R)) B_R = 0 - if (present(B_z)) B_z = 0 - if (present(B_phi)) B_phi = 0 - return - end if - - call pol_flux(R, z, psi_n, dpsidr, dpsidz) - call pol_curr(psi_n, fpol) - - ! The field in cocos=3 is given by - ! - ! B = F(ψ)∇φ + ∇ψ×∇φ. - ! - ! Writing the gradient of ψ=ψ(R,z) as - ! - ! ∇ψ = ∂ψ/∂R ∇R + ∂ψ/∂z ∇z, - ! - ! and carrying out the cross products: - ! - ! B = F(ψ)∇φ - ∂ψ/∂z ∇R/R + ∂ψ/∂R ∇z/R - ! - if (present(B_R)) B_R = - 1/R * dpsidz * psia - if (present(B_z)) B_z = + 1/R * dpsidr * psia - if (present(B_phi)) B_phi = fpol / R - end subroutine bfield - - - function tor_curr(R, z) result(J_phi) - ! Computes the toroidal current J_φ as a function of (R, z) - use const_and_precisions, only : mu0_ - - ! function arguments - real(wp_), intent(in) :: R, z - real(wp_) :: J_phi - - ! local variables - real(wp_) :: dB_Rdz, dB_zdR ! derivatives of B_R, B_z - real(wp_) :: dpsidr, ddpsidrr, ddpsidzz ! derivatives of ψ_n - - call pol_flux(R, z, dpsidr=dpsidr, ddpsidrr=ddpsidrr, ddpsidzz=ddpsidzz) - - ! In the usual MHD limit we have ∇×B = μ₀J. Using the - ! curl in cylindrical coords the toroidal current is - ! - ! J_φ = 1/μ₀ (∇×B)_φ = 1/μ₀ [∂B_R/∂z - ∂B_z/∂R]. - ! - ! Finally, from B = F(ψ)∇φ + ∇ψ×∇φ we find: - ! - ! B_R = - 1/R ∂ψ/∂z, - ! B_z = + 1/R ∂ψ/∂R, - ! - ! from which: - ! - ! ∂B_R/∂z = - 1/R ∂²ψ/∂z² - ! ∂B_z/∂R = + 1/R ∂²ψ/∂R² - 1/R² ∂ψ/∂R. - ! - dB_Rdz = - 1/R * ddpsidzz * psia - dB_zdR = + 1/R * (ddpsidrr - 1/R * dpsidr) * psia - J_phi = 1/mu0_ * (dB_Rdz - dB_zdR) - end function tor_curr - - - subroutine find_ox_point(R0, z0, R1, z1, psi1) - ! Given the point (R₀,z₀) as an initial guess, finds - ! the exact location (R₁,z₁) where ∇ψ(R₁,z₁) = 0. - ! It also returns ψ₁=ψ(R₁,z₁). - ! - ! Note: this is used to find the magnetic X and O point, - ! because both are stationary points for ψ(R,z). - use const_and_precisions, only : comp_eps - use minpack, only : hybrj1 - use logger, only : log_error, log_debug - - ! subroutine arguments - real(wp_), intent(in) :: R0, z0 - real(wp_), intent(out) :: R1, z1, psi1 - - ! local variables - integer :: info - real(wp_) :: sol(2), f(2), df(2,2), wa(15) - character(256) :: msg - - sol = [R0, z0] ! first guess - - call hybrj1(equation, n=2, x=sol, fvec=f, fjac=df, ldfjac=2, & - tol=sqrt(comp_eps), info=info, wa=wa, lwa=15) - if (info /= 1) then - write (msg, '("guess:", 2(x,", ",g0.3))') R0, z0 - call log_debug(msg, mod='equilibrium', proc='find_ox_point') - write (msg, '("solution:", 2(x,", ",g0.3))') sol - call log_debug(msg, mod='equilibrium', proc='find_ox_point') - write (msg, '("hybrj1 failed with error ",g0)') info - call log_error(msg, mod='equilibrium', proc='find_ox_point') - end if - - R1 = sol(1) ! solution - z1 = sol(2) - call pol_flux(R1, z1, psi1) - - contains - - subroutine equation(n, x, f, df, ldf, flag) - ! The equation to solve: f(R,z) = ∇ψ(R,z) = 0 - - ! subroutine arguments - integer, intent(in) :: n, flag, ldf - real(wp_), intent(in) :: x(n) - real(wp_), intent(inout) :: f(n), df(ldf,n) - - if (flag == 1) then - ! return f(R,z) = ∇ψ(R,z) - call pol_flux(R=x(1), z=x(2), dpsidr=f(1), dpsidz=f(2)) - else - ! return ∇f(R,z) = ∇∇ψ(R,z) - call pol_flux(R=x(1), z=x(2), ddpsidrr=df(1,1), & - ddpsidzz=df(2,2), ddpsidrz=df(1,2)) - df(2,1) = df(1,2) - end if - end subroutine equation - - end subroutine find_ox_point - - - subroutine find_htg_point(R0, z0, R1, z1, psi0) - ! Given the point (R₀,z₀) as an initial guess, finds - ! the exact location (R₁,z₁) where: - ! { ψ(R₁,z₁) = ψ₀ - ! { ∂ψ/∂R(R₁,z₁) = 0 . - ! - ! Note: this is used to find the horizontal tangent - ! point of the contour ψ(R,z)=ψ₀. - use const_and_precisions, only : comp_eps - use minpack, only : hybrj1 - use logger, only : log_error, log_debug - - ! subroutine arguments - real(wp_), intent(in) :: R0, z0, psi0 - real(wp_), intent(out) :: R1, z1 - - ! local variables - integer :: info - real(wp_) :: sol(2), f(2), df(2,2), wa(15) - character(256) :: msg - - sol = [R0, z0] ! first guess - - call hybrj1(equation, n=2, x=sol, fvec=f, fjac=df, ldfjac=2, & - tol=sqrt(comp_eps), info=info, wa=wa, lwa=15) - if (info /= 1) then - write (msg, '("guess:", 2(x,", ",g0.3))') R0, z0 - call log_debug(msg, mod='equilibrium', proc='find_htg_point') - write (msg, '("solution:", 2(x,", ",g0.3))') sol - call log_debug(msg, mod='equilibrium', proc='find_htg_point') - write (msg, '("hybrj1 failed with error ",g0)') info - call log_error(msg, mod='equilibrium', proc='find_htg_point') - end if - - R1 = sol(1) ! solution - z1 = sol(2) - contains - - subroutine equation(n, x, f, df, ldf, flag) - ! The equation to solve: f(R,z) = [ψ(R,z)-ψ₀, ∂ψ/∂R] = 0 - - ! subroutine arguments - integer, intent(in) :: n, ldf, flag - real(wp_), intent(in) :: x(n) - real(wp_), intent(inout) :: f(n), df(ldf, n) - - if (flag == 1) then - ! return f(R,z) = [ψ(R,z)-ψ₀, ∂ψ/∂R] - call pol_flux(R=x(1), z=x(2), psi_n=f(1), dpsidr=f(2)) - f(1) = f(1) - psi0 - else - ! return ∇f(R,z) = [[∂ψ/∂R, ∂ψ/∂z], [∂²ψ/∂R², ∂²ψ/∂R∂z]] - call pol_flux(R=x(1), z=x(2), dpsidr=df(1,1), dpsidz=df(1,2), & - ddpsidrr=df(2,1), ddpsidrz=df(2,2)) - end if - end subroutine equation - - end subroutine find_htg_point - - - subroutine unset_equil_spline - ! Unsets the splines global variables, see the top of this file. - - call fpol_spline%deinit - call psi_spline%deinit - call q_spline%deinit - call rhop_spline%deinit - call rhot_spline%deinit - - if (allocated(psi_domain%R)) deallocate(psi_domain%R, psi_domain%z) - - end subroutine unset_equil_spline - -end module equilibrium diff --git a/src/gray_core.f90 b/src/gray_core.f90 index dd244d8..b12d988 100644 --- a/src/gray_core.f90 +++ b/src/gray_core.f90 @@ -7,11 +7,12 @@ module gray_core contains - subroutine gray_main(params, data, plasma, results, error, rhout) + subroutine gray_main(params, equil, plasma, limiter, results, error, rhout) use const_and_precisions, only : zero, one, comp_tiny use polarization, only : ellipse_to_field - use types, only : table, wrap - use gray_params, only : gray_parameters, gray_data, gray_results, EQ_VACUUM + use types, only : contour, table, wrap + use gray_params, only : gray_parameters, gray_results, EQ_VACUUM + use gray_equil, only : abstract_equil use gray_plasma, only : abstract_plasma use gray_tables, only : init_tables, store_ec_profiles, store_ray_data, & store_beam_shape, find_flux_surfaces, & @@ -29,8 +30,9 @@ contains ! subroutine arguments type(gray_parameters), intent(inout) :: params - type(gray_data), intent(in) :: data + class(abstract_equil), intent(in) :: equil class(abstract_plasma), intent(in) :: plasma + type(contour), intent(in) :: limiter type(gray_results), intent(out) :: results ! Predefined grid for the output profiles (optional) @@ -115,10 +117,10 @@ contains if (params%equilibrium%iequil /= EQ_VACUUM) then ! Initialise the magsurf_data module - call compute_flux_averages(params, results%tables) + call compute_flux_averages(params, equil, results%tables) ! Initialise the output profiles - call pec_init(params%output, rhout) + call pec_init(params%output, equil, rhout) end if ! Allocate memory for the results... @@ -133,14 +135,14 @@ contains ! ========= set environment END ========= ! Pre-determinted tables - results%tables%kinetic_profiles = kinetic_profiles(params, plasma) - results%tables%ec_resonance = ec_resonance(params, bres) - results%tables%inputs_maps = inputs_maps(params, plasma, bres, xgcn) + results%tables%kinetic_profiles = kinetic_profiles(params, equil, plasma) + results%tables%ec_resonance = ec_resonance(params, equil, bres) + results%tables%inputs_maps = inputs_maps(params, equil, plasma, bres, xgcn) ! print ψ surface for q=3/2 and q=2/1 - call find_flux_surfaces( & - qvals=[1.5_wp_, 2.0_wp_], params=params, & - tbl=results%tables%flux_surfaces) + call find_flux_surfaces(qvals=[1.5_wp_, 2.0_wp_], & + params=params, equil=equil, & + tbl=results%tables%flux_surfaces) ! print initial position write (msg, '("initial position:",3(x,g0.3))') params%antenna%pos @@ -218,14 +220,14 @@ contains lgcpl1 = zero p0ray = p0jk ! * initial beam power - call ic_gb(params, anv0, ak0, yw, ypw, stv, xc, du1, & - gri, ggri, index_rt, results%tables) ! * initial conditions + call ic_gb(params, equil, anv0, ak0, yw, ypw, stv, xc, & + du1, gri, ggri, index_rt, results%tables) ! * initial conditions do jk=1,params%raytracing%nray zzm = yw(3,jk)*0.01_wp_ rrm = sqrt(yw(1,jk)*yw(1,jk)+yw(2,jk)*yw(2,jk))*0.01_wp_ - if (data%equilibrium%limiter%contains(rrm, zzm)) then ! * start propagation in/outside vessel? + if (limiter%contains(rrm, zzm)) then ! * start propagation in/outside vessel? iow(jk) = 1 ! + inside else iow(jk) = 0 ! + outside @@ -258,7 +260,7 @@ contains do jk=1,params%raytracing%nray if(iwait(jk)) cycle ! jk ray is waiting for next pass stv(jk) = stv(jk) + params%raytracing%dst ! current ray step - call rkstep(params, plasma, sox, bres, xgcn, yw(:,jk), & + call rkstep(params, equil, plasma, sox, bres, xgcn, yw(:,jk), & ypw(:,jk), gri(:,jk), ggri(:,:,jk), igrad_b) end do ! update position and grad @@ -274,10 +276,11 @@ contains ! compute derivatives with updated gradient and local plasma values xv = yw(1:3,jk) anv = yw(4:6,jk) - call ywppla_upd(params, plasma, xv, anv, gri(:,jk), ggri(:,:,jk), & - sox, bres, xgcn,ypw(:,jk), psinv, dens, btot, bv, & - xg, yg, derxg, anpl, anpr, ddr, ddi, dersdst, & - derdnm, ierrn, igrad_b) + call ywppla_upd(params, equil, plasma, & + xv, anv, gri(:,jk), ggri(:,:,jk), sox, bres, & + xgcn,ypw(:,jk), psinv, dens, btot, bv, xg, yg, & + derxg, anpl, anpr, ddr, ddi, dersdst, derdnm, & + ierrn, igrad_b) ! update global error code and print message if(ierrn/=0) then error = ior(error,ierrn) @@ -289,7 +292,7 @@ contains rrm = sqrt(xv(1)*xv(1)+xv(2)*xv(2))*0.01_wp_ ins_pl = (psinv>=zero .and. psinvin else if(ext_wl) then ! ray exit vessel - call wall_out(jk, data%equilibrium%limiter, ins_pl, xv, anv, & + call wall_out(jk, equil, limiter, ins_pl, xv, anv, & params%raytracing%dst, bres, sox, psipol, chipol, & iow, iop, ext, eyt) yw(:,jk) = [xv, anv] ! * updated coordinates (reflected) igrad_b = 0 ! * switch to ray-tracing - call ywppla_upd(params, plasma, xv, anv, gri(:,jk), ggri(:,:,jk), & - sox, bres, xgcn, ypw(:,jk), psinv, dens, btot, & - bv, xg, yg, derxg, anpl, anpr, ddr, ddi, dersdst, & + call ywppla_upd(params, equil, plasma, & + xv, anv, gri(:,jk), ggri(:,:,jk), sox, bres, & + xgcn, ypw(:,jk), psinv, dens, btot, bv, xg, & + yg, derxg, anpl, anpr, ddr, ddi, dersdst, & derdnm, ierrn, igrad_b) ! * update derivatives after reflection if(ierrn/=0) then ! * update global error code and print message error = ior(error,ierrn) @@ -396,8 +401,8 @@ contains yypnext(:,jk,index_rt) = ypw(:,jk) ! for next pass = reflection point stnext(jk,index_rt) = stv(jk) ! . starting step for next pass = step after reflection - call plasma_in(jk, xv, anv, bres, sox, cpl, psipol, chipol, & ! . ray re-enters plasma after reflection - iop, ext, eyt, perfect=.false.) + call plasma_in(jk, equil, xv, anv, bres, sox, cpl, & ! . ray re-enters plasma after reflection + psipol, chipol, iop, ext, eyt, perfect=.false.) if(cpl(1) < etaucr) then ! . low coupled power for O-mode? => de-activate derived rays call turnoffray(jk,ip+1,npass,2*ib-1,iroff) @@ -434,10 +439,10 @@ contains tekev = plasma%temp(psinv) block complex(wp_) :: Npr_warm - call alpha_effj(params%ecrh_cd, plasma, psinv, xg, yg, dens, & - tekev, ak0, bres, derdnm, anpl, anpr, sox, & - Npr_warm, alpha, didp, nharm, nhf, iokhawa, & - ierrhcd) + call alpha_effj(params%ecrh_cd, equil, plasma, & + psinv, xg, yg, dens, tekev, ak0, bres, & + derdnm, anpl, anpr, sox, Npr_warm, alpha, & + didp, nharm, nhf, iokhawa, ierrhcd) anprre = Npr_warm%re anprim = Npr_warm%im if(ierrhcd /= 0) then @@ -481,7 +486,7 @@ contains ccci(jk,i:params%raytracing%nstep) = ccci(jk,i-1) psjki(jk,i:params%raytracing%nstep) = psjki(jk,i-1) else - call store_ray_data(params, results%tables, & + call store_ray_data(params, equil, results%tables, & i, jk, stv(jk), p0ray(jk), xv, psinv, btot, bv, ak0, & anpl, anpr, anv, anprim, dens, tekev, alpha, tau, dids, & nharm, nhf, iokhawa, index_rt, ddr, ddi, xg, yg, derxg) ! p0ray/etau1 [dids normalization] = fraction of p0 coupled to this ray (not including absorption from previous passes) @@ -583,9 +588,9 @@ contains results%tables%ec_profiles, rhop_tab, rhot_tab, jphi_beam, & jcd_beam, dpdv_beam, currins_beam, pins_beam, ip) - call postproc_profiles(pabs_beam,icd_beam,rhot_tab,dpdv_beam, & - jphi_beam, rhotpav,drhotpav,rhotjava,drhotjava,dpdvp,jphip, & - rhotp,drhotp,rhotj,drhotj,dpdvmx,jphimx,ratjamx,ratjbmx) ! *compute profiles width for current beam + call postproc_profiles(equil, pabs_beam, icd_beam, rhot_tab, dpdv_beam, & + jphi_beam, rhotpav, drhotpav, rhotjava, drhotjava, dpdvp, jphip, & + rhotp, drhotp, rhotj, drhotj, dpdvmx, jphimx, ratjamx, ratjbmx) ! *compute profiles width for current beam if (results%tables%summary%active) then call results%tables%summary%append([ & @@ -665,18 +670,20 @@ contains end subroutine vectinit - subroutine ic_gb(params, anv0c, ak0, ywrk0, ypwrk0, & + subroutine ic_gb(params, equil, anv0c, ak0, ywrk0, ypwrk0, & stv, xc0, du10, gri, ggri, index_rt, & tables) ! beam tracing initial conditions igrad=1 ! !!!!!! check ray tracing initial conditions igrad=0 !!!!!! use const_and_precisions, only : zero,one,pi,half,two,degree,ui=>im use gray_params, only : gray_parameters, gray_tables + use gray_equil, only : abstract_equil use types, only : table use gray_tables, only : store_ray_data ! subroutine arguments type(gray_parameters), intent(in) :: params + class(abstract_equil), intent(in) :: equil real(wp_), dimension(3), intent(in) :: anv0c real(wp_), intent(in) :: ak0 real(wp_), dimension(6, params%raytracing%nray), intent(out) :: ywrk0, ypwrk0 @@ -984,7 +991,7 @@ contains ! save step "zero" data if (present(tables)) & - call store_ray_data(params, tables, & + call store_ray_data(params, equil, tables, & i=0, jk=jk, s=stv(jk), P0=one, pos=xc0(:,k,j), & psi_n=-one, B=zero, b_n=[zero,zero,zero], k0=ak0, & N_pl=zero, N_pr=zero, N=ywrk0(:, jk), N_pr_im=zero, & @@ -998,21 +1005,23 @@ contains - subroutine rkstep(params, plasma, sox, bres, xgcn, y, yp, dgr, ddgr, igrad) + subroutine rkstep(params, equil, plasma, & + sox, bres, xgcn, y, yp, dgr, ddgr, igrad) ! Runge-Kutta integrator use gray_params, only : gray_parameters + use gray_equil, only : abstract_equil use gray_plasma, only : abstract_plasma ! subroutine arguments - type(gray_parameters), intent(in) :: params - class(abstract_plasma), intent(in) :: plasma - - real(wp_), intent(in) :: bres, xgcn - real(wp_), intent(inout) :: y(6) - real(wp_), intent(in) :: yp(6) - real(wp_), intent(in) :: dgr(3) - real(wp_), intent(in) :: ddgr(3,3) - integer, intent(in) :: igrad,sox + type(gray_parameters), intent(in) :: params + class(abstract_equil), intent(in) :: equil + class(abstract_plasma), intent(in) :: plasma + real(wp_), intent(in) :: bres, xgcn + real(wp_), intent(inout) :: y(6) + real(wp_), intent(in) :: yp(6) + real(wp_), intent(in) :: dgr(3) + real(wp_), intent(in) :: ddgr(3,3) + integer, intent(in) :: igrad,sox ! local variables real(wp_), dimension(6) :: k1, k2, k3, k4 @@ -1030,19 +1039,22 @@ contains function f(y) real(wp_), intent(in) :: y(6) real(wp_) :: f(6) - call rhs(params, plasma, sox, bres, xgcn, y, dgr, ddgr, f, igrad) + call rhs(params, equil, plasma, sox, bres, xgcn, y, dgr, ddgr, f, igrad) end function end subroutine rkstep - subroutine rhs(params, plasma, sox, bres, xgcn, y, dgr, ddgr, dery, igrad) + subroutine rhs(params, equil, plasma, & + sox, bres, xgcn, y, dgr, ddgr, dery, igrad) ! Compute right-hand side terms of the ray equations (dery) ! used in R-K integrator use gray_params, only : gray_parameters + use gray_equil, only : abstract_equil use gray_plasma, only : abstract_plasma ! subroutine arguments type(gray_parameters), intent(in) :: params + class(abstract_equil), intent(in) :: equil class(abstract_plasma), intent(in) :: plasma real(wp_), intent(in) :: y(6) real(wp_), intent(in) :: bres, xgcn @@ -1057,7 +1069,7 @@ contains real(wp_), dimension(3,3) :: derbv xv = y(1:3) - call plas_deriv(params, plasma, xv, bres, xgcn, dens, btot, & + call plas_deriv(equil, plasma, xv, bres, xgcn, dens, btot, & bv, derbv, xg, yg, derxg, deryg) anv = y(4:6) call disp_deriv(params, anv, sox, xg, yg, derxg, deryg, & @@ -1065,18 +1077,21 @@ contains end subroutine rhs - subroutine ywppla_upd(params, plasma, xv, anv, dgr, ddgr, sox, bres, xgcn, dery, & - psinv, dens, btot, bv, xg, yg, derxg, anpl, anpr, & - ddr, ddi, dersdst, derdnm, error, igrad) + subroutine ywppla_upd(params, equil, plasma, & + xv, anv, dgr, ddgr, sox, bres, xgcn, dery, & + psinv, dens, btot, bv, xg, yg, derxg, anpl, & + anpr, ddr, ddi, dersdst, derdnm, error, igrad) ! Compute right-hand side terms of the ray equations (dery) ! used after full R-K step and grad(S_I) update use gray_errors, only : raise_error, large_npl use gray_params, only : gray_parameters + use gray_equil, only : abstract_equil use gray_plasma, only : abstract_plasma ! subroutine rguments - type(gray_parameters), intent(in) :: params - class(abstract_plasma), intent(in) :: plasma + type(gray_parameters), intent(in) :: params + class(abstract_equil), intent(in) :: equil + class(abstract_plasma), intent(in) :: plasma real(wp_), intent(in) :: xv(3), anv(3) real(wp_), intent(in) :: dgr(3) @@ -1096,9 +1111,10 @@ contains real(wp_), dimension(3,3) :: derbv real(wp_), parameter :: anplth1 = 0.99_wp_, anplth2 = 1.05_wp_ - call plas_deriv(params, plasma, xv,bres,xgcn,dens,btot,bv,derbv,xg,yg,derxg,deryg,psinv) - call disp_deriv(params, anv,sox,xg,yg,derxg,deryg,bv,derbv,dgr,ddgr,igrad, & - dery,anpl,anpr,ddr,ddi,dersdst,derdnm) + call plas_deriv(equil, plasma, xv, bres, xgcn, dens, btot, & + bv, derbv, xg, yg, derxg, deryg, psinv) + call disp_deriv(params, anv, sox, xg, yg, derxg, deryg, bv, derbv, dgr, & + ddgr, igrad, dery, anpl, anpr, ddr, ddi, dersdst, derdnm) if (abs(anpl) > anplth2) then error = raise_error(error, large_npl, 1) @@ -1315,15 +1331,14 @@ contains end subroutine solg3 - subroutine plas_deriv(params, plasma, xv, bres, xgcn, dens, btot, bv, & - derbv, xg, yg, derxg, deryg, psinv_opt) + subroutine plas_deriv(equil, plasma, xv, bres, xgcn, dens, btot, & + bv, derbv, xg, yg, derxg, deryg, psinv_opt) use const_and_precisions, only : zero, cm - use gray_params, only : gray_parameters, EQ_VACUUM + use gray_equil, only : abstract_equil, vacuum use gray_plasma, only : abstract_plasma - use equilibrium, only : psia, pol_flux, pol_curr, sgnbphi ! subroutine arguments - type(gray_parameters), intent(in) :: params + class(abstract_equil), intent(in) :: equil class(abstract_plasma), intent(in) :: plasma real(wp_), dimension(3), intent(in) :: xv real(wp_), intent(in) :: xgcn, bres @@ -1352,11 +1367,12 @@ contains bv = zero derbv = zero - if (params%equilibrium%iequil == EQ_VACUUM) then - ! copy optional output - if (present(psinv_opt)) psinv_opt = psinv - return - end if + select type (equil) + type is (vacuum) + ! copy optional output + if (present(psinv_opt)) psinv_opt = psinv + return + end select dbtot = zero dbv = zero @@ -1374,15 +1390,16 @@ contains csphi = xx/rr snphi = yy/rr - bv(1) = -snphi*sgnbphi - bv(2) = csphi*sgnbphi + bv(1) = -snphi*equil%sgn_bphi + bv(2) = csphi*equil%sgn_bphi ! convert from cm to meters zzm = 1.0e-2_wp_*zz rrm = 1.0e-2_wp_*rr - call pol_flux(rrm, zzm, psinv, dpsidr, dpsidz, ddpsidrr, ddpsidzz, ddpsidrz) - call pol_curr(psinv, fpolv, dfpv) + call equil%pol_flux(rrm, zzm, psinv, dpsidr, dpsidz, & + ddpsidrr, ddpsidzz, ddpsidrz) + call equil%pol_curr(psinv, fpolv, dfpv) ! copy optional output if (present(psinv_opt)) psinv_opt = psinv @@ -1402,8 +1419,8 @@ contains ! B = f(psi)/R e_phi+ grad psi x e_phi/R bphi = fpolv/rrm - brr = -dpsidz*psia/rrm - bzz = +dpsidr*psia/rrm + brr = -dpsidz*equil%psi_a/rrm + bzz = +dpsidr*equil%psi_a/rrm ! bvc(i) = B_i in cylindrical coordinates bvc = [brr, bphi, bzz] @@ -1414,12 +1431,12 @@ contains bv(3)=bvc(3) ! dbvcdc(iv,jv) = d Bcil(iv) / dxvcil(jv) - dbvcdc(1,1) = -ddpsidrz*psia/rrm - brr/rrm - dbvcdc(1,3) = -ddpsidzz*psia/rrm + dbvcdc(1,1) = -ddpsidrz*equil%psi_a/rrm - brr/rrm + dbvcdc(1,3) = -ddpsidzz*equil%psi_a/rrm dbvcdc(2,1) = dfpv*dpsidr/rrm - bphi/rrm dbvcdc(2,3) = dfpv*dpsidz/rrm - dbvcdc(3,1) = +ddpsidrr*psia/rrm - bzz/rrm - dbvcdc(3,3) = +ddpsidrz*psia/rrm + dbvcdc(3,1) = +ddpsidrr*equil%psi_a/rrm - bzz/rrm + dbvcdc(3,3) = +ddpsidrz*equil%psi_a/rrm ! dbvdc(iv,jv) = d Bcart(iv) / dxvcil(jv) dbvdc(1,1) = dbvcdc(1,1)*csphi - dbvcdc(2,1)*snphi @@ -1746,15 +1763,15 @@ contains - subroutine alpha_effj(params, plasma, psinv, X, Y, density, temperature, & - k0, Bres, derdnm, Npl, Npr_cold, sox, Npr, & - alpha, dIdp, nhmin, nhmax, iokhawa, error) + subroutine alpha_effj(params, equil, plasma, psinv, X, Y, density, & + temperature, k0, Bres, derdnm, Npl, Npr_cold, & + sox, Npr, alpha, dIdp, nhmin, nhmax, iokhawa, error) ! Computes the absorption coefficient and effective current use const_and_precisions, only : pi, mc2=>mc2_ use gray_params, only : ecrh_cd_parameters + use gray_equil, only : abstract_equil use gray_plasma, only : abstract_plasma - use equilibrium, only : sgnbphi use dispersion, only : harmnumber, warmdisp use eccd, only : setcdcoeff, eccdeff, fjch0, fjch, fjncl use gray_errors, only : negative_absorption, raise_error @@ -1767,8 +1784,9 @@ contains ! ECRH & CD parameters type(ecrh_cd_parameters), intent(in) :: params - ! plasma object - class(abstract_plasma), intent(in) :: plasma + ! MHD equilibrium, plasma object + class(abstract_equil), intent(in) :: equil + class(abstract_plasma), intent(in) :: plasma ! poloidal flux ψ real(wp_), intent(in) :: psinv ! CMA diagram variables: X=(ω_pe/ω)², Y=ω_ce/ω @@ -1899,7 +1917,7 @@ contains ! current drive efficiency [A⋅m/W] effjcdav = rbavi*effjcd - dIdp = sgnbphi * effjcdav / (2*pi*rrii) + dIdp = equil%sgn_bphi * effjcdav / (2*pi*rrii) end subroutine alpha_effj diff --git a/src/gray_equil.f90 b/src/gray_equil.f90 new file mode 100644 index 0000000..c1f1c5a --- /dev/null +++ b/src/gray_equil.f90 @@ -0,0 +1,1620 @@ +! This module handles the loading, interpolation and evaluation of the +! MHD equilibrium data (poloidal current function, safety factor, +! poloidal flux, magnetic field, etc.) +! +module gray_equil + use const_and_precisions, only : wp_, comp_huge + use splines, only : spline_simple, spline_1d, spline_2d, linear_1d + use types, only : contour + + implicit none + + ! macro for suppressing unused variable warnings +# define unused(x) associate(x => x); end associate + + type, abstract :: abstract_equil + ! Generic equilibrium interface + real(wp_) :: psi_a = 0 ! Poloidal flux at the edge minus flux on axis, ψ_a + real(wp_) :: phi_a = 0 ! Toroidal flux at the edge (r=a), Φ_a + real(wp_) :: b_axis = 0 ! Value of B_φ at the magnetic axis (used in J_cd def) + real(wp_) :: b_centre = 0 ! Value of B_φ at R_centre (used in Jcd_astra def) + real(wp_) :: r_centre = 1 ! Alternative reference radius for B_φ + real(wp_) :: sgn_bphi = 0 ! Sign of B_φ (>0 means counter-clockwise) + real(wp_) :: axis(2) = [0, 0] ! Magnetic axis position (R₀, z₀) + real(wp_) :: r_range(2) = [-comp_huge, comp_huge] ! R range of the equilibrium domain + real(wp_) :: z_range(2) = [-comp_huge, comp_huge] ! z range of the equilibrium domain + real(wp_) :: z_boundary(2) = [0, 0] ! z range of the plasma boundary + contains + procedure(pol_flux_sub), deferred :: pol_flux + procedure(pol_curr_sub), deferred :: pol_curr + procedure(safety_fun), deferred :: safety + procedure(rho_conv_fun), deferred :: pol2tor, tor2pol + procedure(flux_contour_sub), deferred :: flux_contour + procedure :: b_field + procedure :: tor_curr + end type + + abstract interface + subroutine pol_flux_sub(self, R, z, psi_n, dpsidr, dpsidz, & + ddpsidrr, ddpsidzz, ddpsidrz) + ! Computes the normalised poloidal flux ψ_n and its + ! derivatives wrt (R, z) up to the second order. + ! + ! Note: all output arguments are optional. + import :: abstract_equil, wp_ + class(abstract_equil), intent(in) :: self + real(wp_), intent(in) :: R, z + real(wp_), intent(out), optional :: & + psi_n, dpsidr, dpsidz, ddpsidrr, ddpsidzz, ddpsidrz + end subroutine pol_flux_sub + + subroutine pol_curr_sub(self, psi_n, fpol, dfpol) + ! Computes the poloidal current function F(ψ_n) + ! and (optionally) its derivative dF/dψ_n given ψ_n. + import :: abstract_equil, wp_ + class(abstract_equil), intent(in) :: self + real(wp_), intent(in) :: psi_n ! normalised poloidal flux + real(wp_), intent(out) :: fpol ! poloidal current + real(wp_), intent(out), optional :: dfpol ! derivative + end subroutine pol_curr_sub + + function safety_fun(self, psi_n) result(q) + ! Computes the safety factor q as a function of the + ! normalised poloidal flux ψ_n. + ! + ! Note: this returns the absolute value of q. + import :: abstract_equil, wp_ + class(abstract_equil), intent(in) :: self + real(wp_), intent(in) :: psi_n + real(wp_) :: q + end function safety_fun + + function rho_conv_fun(self, rho_in) result(rho_out) + ! Converts between poloidal (ρ_p) and toroidal (ρ_t) normalised radius + import :: abstract_equil, wp_ + class(abstract_equil), intent(in) :: self + real(wp_), intent(in) :: rho_in + real(wp_) :: rho_out + end function rho_conv_fun + + subroutine flux_contour_sub(self, psi0, R_min, R, z, & + R_hi, z_hi, R_lo, z_lo) + ! Computes a contour of the ψ(R,z)=ψ₀ flux surface. + ! Notes: + ! - R,z are the contour points + ! - R_min is a value such that R>R_min for any contour point + ! - (R,z)_hi and (R,z)_lo are a guess for the higher and lower + ! horizontal-tangent points of the contour. These variables + ! will be updated with the exact value on success. + import :: abstract_equil, wp_ + class(abstract_equil), intent(in) :: self + real(wp_), intent(in) :: psi0 + real(wp_), intent(in) :: R_min + real(wp_), intent(out) :: R(:), z(:) + real(wp_), intent(inout) :: R_hi, z_hi, R_lo, z_lo + end subroutine flux_contour_sub + end interface + + + type, extends(abstract_equil) :: analytic_equil + ! Analytical equilibrium + private + real(wp_) :: q0 ! Safety factor at the magnetic axis + real(wp_) :: q1 ! Safety factor at the edge + real(wp_) :: alpha ! Exponent for the q(ρ_p) power law + real(wp_) :: R0 ! R of the magnetic axis (m) + real(wp_) :: z0 ! z of the magnetic axis (m) + real(wp_) :: a ! Minor radius (m) + real(wp_) :: B0 ! Magnetic field at the magnetic axis (T) + contains + procedure :: pol_flux => analytic_pol_flux + procedure :: pol_curr => analytic_pol_curr + procedure :: safety => analytic_safety + procedure :: tor2pol => analytic_tor2pol + procedure :: pol2tor => analytic_pol2tor + procedure :: flux_contour => analytic_flux_contour + end type + + + type, extends(abstract_equil) :: vacuum + ! Vacuum + contains + procedure :: pol_flux => vacuum_pol_flux + procedure :: pol_curr => vacuum_pol_curr + procedure :: safety => vacuum_safety + procedure :: tor2pol => vacuum_conv + procedure :: pol2tor => vacuum_conv + procedure :: flux_contour => vacuum_flux_contour + end type + + + type, extends(abstract_equil) :: numeric_equil + ! Numerical equilibrium + private + real(wp_) :: fpol_a ! Poloidal current at the edge (r=a), F_a + ! Splines + type(spline_2d) :: psi_spline + type(contour) :: psi_domain + type(spline_1d) :: fpol_spline + type(spline_simple) :: q_spline + type(linear_1d) :: rhop_spline, rhot_spline + contains + procedure :: pol_flux => numeric_pol_flux + procedure :: pol_curr => numeric_pol_curr + procedure :: safety => numeric_safety + procedure :: tor2pol => numeric_tor2pol + procedure :: pol2tor => numeric_pol2tor + procedure :: flux_contour => numeric_flux_contour + procedure :: init => numeric_init + procedure :: find_ox_point + procedure :: find_htg_point + end type + + + type eqdsk_data + ! MHD equilibrium data from G-EQDSK file format + real(wp_), allocatable :: grid_r(:) ! R values of the uniform grid + real(wp_), allocatable :: grid_z(:) ! z values of the uniform grid + real(wp_), allocatable :: fpol(:) ! Poloidal current function, F(ψ_n) + real(wp_), allocatable :: q(:) ! Safety factor, q(ψ_n) + real(wp_), allocatable :: psi(:) ! Normalised poloidal flux, ψ_n(R) + real(wp_), allocatable :: psi_map(:,:) ! Normalised poloidal flux 2D map, ψ(R, z) + real(wp_) :: psi_a ! Poloidal flux at the edge minus flux on axis, ψ_a + real(wp_) :: r_ref ! Reference R₀ (B = B₀R₀/R without the plasma) + real(wp_) :: axis(2) ! Magnetic axis position (R₀, z₀) + type(contour) :: limiter ! limiter contour (wall) + type(contour) :: boundary ! boundary contour (plasma) + end type + + + private + public abstract_equil ! The abstract equilibrium object + public analytic_equil, numeric_equil, vacuum ! Implementations + public load_equil ! To load equilibrium from file + public eqdsk_data ! G-EQDSK data structure + public contour ! re-export contours eqdsk_data + +contains + + subroutine b_field(self, R, z, B_R, B_z, B_phi) + ! Computes the magnetic field as a function of + ! (R, z) in cylindrical coordinates + ! + ! Note: all output arguments are optional. + + ! subroutine arguments + class(abstract_equil), intent(in) :: self + real(wp_), intent(in) :: R, z + real(wp_), intent(out), optional :: B_R, B_z, B_phi + + ! local variables + real(wp_) :: psi_n, fpol, dpsidr, dpsidz + + call self%pol_flux(R, z, psi_n, dpsidr, dpsidz) + call self%pol_curr(psi_n, fpol) + + ! The field in cocos=3 is given by + ! + ! B = F(ψ)∇φ + ∇ψ×∇φ. + ! + ! Writing the gradient of ψ=ψ(R,z) as + ! + ! ∇ψ = ∂ψ/∂R ∇R + ∂ψ/∂z ∇z, + ! + ! and carrying out the cross products: + ! + ! B = F(ψ)∇φ - ∂ψ/∂z ∇R/R + ∂ψ/∂R ∇z/R + ! + if (present(B_R)) B_R = - 1/R * dpsidz * self%psi_a + if (present(B_z)) B_z = + 1/R * dpsidr * self%psi_a + if (present(B_phi)) B_phi = fpol / R + end subroutine b_field + + + function tor_curr(self, R, z) result(J_phi) + ! Computes the toroidal current J_φ as a function of (R, z) + use const_and_precisions, only : mu0_ + + ! function arguments + class(abstract_equil), intent(in) :: self + real(wp_), intent(in) :: R, z + real(wp_) :: J_phi + + ! local variables + real(wp_) :: dB_Rdz, dB_zdR ! derivatives of B_R, B_z + real(wp_) :: dpsidr, ddpsidrr, ddpsidzz ! derivatives of ψ_n + + call self%pol_flux(R, z, dpsidr=dpsidr, ddpsidrr=ddpsidrr, ddpsidzz=ddpsidzz) + + ! In the usual MHD limit we have ∇×B = μ₀J. Using the + ! curl in cylindrical coords the toroidal current is + ! + ! J_φ = 1/μ₀ (∇×B)_φ = 1/μ₀ [∂B_R/∂z - ∂B_z/∂R]. + ! + ! Finally, from B = F(ψ)∇φ + ∇ψ×∇φ we find: + ! + ! B_R = - 1/R ∂ψ/∂z, + ! B_z = + 1/R ∂ψ/∂R, + ! + ! from which: + ! + ! ∂B_R/∂z = - 1/R ∂²ψ/∂z² + ! ∂B_z/∂R = + 1/R ∂²ψ/∂R² - 1/R² ∂ψ/∂R. + ! + dB_Rdz = - 1/R * ddpsidzz * self%psi_a + dB_zdR = + 1/R * (ddpsidrr - 1/R * dpsidr) * self%psi_a + J_phi = 1/mu0_ * (dB_Rdz - dB_zdR) + end function tor_curr + + + ! + ! Analytical model + ! + + subroutine analytic_pol_flux(self, R, z, psi_n, dpsidr, dpsidz, & + ddpsidrr, ddpsidzz, ddpsidrz) + use const_and_precisions, only : pi + + ! function arguments + class(analytic_equil), intent(in) :: self + real(wp_), intent(in) :: R, z + real(wp_), intent(out), optional :: & + psi_n, dpsidr, dpsidz, ddpsidrr, ddpsidzz, ddpsidrz + + ! local variables + real(wp_) :: r_g, rho_t, rho_p ! geometric radius, √Φ_n, √ψ_n + real(wp_) :: gamma ! γ = 1/√(1 - r²/R₀²) + real(wp_) :: dpsidphi ! (∂ψ_n/∂Φ_n) + real(wp_) :: ddpsidphidr, ddpsidphidz ! ∇(∂ψ_n/∂Φ_n) + real(wp_) :: phi_n ! Φ_n + real(wp_) :: dphidr, dphidz ! ∇Φ_n + real(wp_) :: ddphidrdr, ddphidzdz ! ∇∇Φ_n + real(wp_) :: ddphidrdz ! ∂²Φ_n/∂R∂z + real(wp_) :: q, dq ! q(ρ_p), Δq=(q₁-q₀)/(α/2 + 1) + real(wp_) :: dqdr, dqdz ! ∇q + real(wp_) :: dphidr2, ddphidr2dr2 ! dΦ_n/d(r²), d²Φ_n/d(r²)² + + ! The normalised poloidal flux ψ_n(R, z) is computed as follows: + ! 1. ψ_n = ρ_p² + ! 2. ρ_p = ρ_p(ρ_t), using `self%tor2pol`, which in turns uses q(ψ) + ! 3. ρ_t = √Φ_n + ! 4. Φ_n = Φ(r_g)/Φ(a), where Φ(r) is the flux of B_φ=B₀R₀/R + ! through a circular surface + ! 5. r_g = √[(R-R₀)²+(z-z₀)²] is the geometric minor radius + r_g = hypot(R - self%R0, z - self%z0) + + ! The exact flux of the toroidal field B_φ = B₀R₀/R is: + ! + ! Φ(r) = B₀πr² 2γ/(γ + 1) where γ=1/√(1 - r²/R₀²). + ! + ! Notes: + ! 1. the function Φ(r) is defined for r≤R₀ only. + ! 2. as r → 0, γ → 1, so Φ ~ B₀πr². + ! 3. as r → 1⁻, Φ → 2B₀πr² but dΦ/dr → -∞. + ! 4. |B_R|, |B_z| → +-∞. + ! + if (r_g > self%R0) then + if (present(psi_n)) psi_n = -1 + if (present(dpsidr)) dpsidr = 0 + if (present(dpsidz)) dpsidz = 0 + if (present(ddpsidrr)) ddpsidrr = 0 + if (present(ddpsidzz)) ddpsidzz = 0 + if (present(ddpsidrz)) ddpsidrz = 0 + return + end if + + gamma = 1 / sqrt(1 - (r_g/self%R0)**2) + phi_n = self%B0 * pi*r_g**2 * 2*gamma/(gamma + 1) / self%phi_a + rho_t = sqrt(phi_n) + rho_p = self%tor2pol(rho_t) + + ! For ∇Φ_n and ∇∇Φ_n we also need: + ! + ! ∂Φ∂(r²) = B₀π γ(r) + ! ∂²Φ∂(r²)² = B₀π γ³(r) / (2 R₀²) + ! + dphidr2 = self%B0 * pi * gamma / self%phi_a + ddphidr2dr2 = self%B0 * pi * gamma**3/(2 * self%R0**2) / self%phi_a + + ! ∇Φ_n = ∂Φ_n/∂(r²) ∇(r²) + ! where ∇(r²) = 2[(R-R₀), (z-z₀)] + dphidr = dphidr2 * 2*(R - self%R0) + dphidz = dphidr2 * 2*(z - self%z0) + + ! ∇∇Φ_n = ∇[∂Φ_n/∂(r²)] ∇(r²) + ∂Φ_n/∂(r²) ∇∇(r²) + ! = ∂²Φ_n/∂(r²)² ∇(r²)∇(r²) + ∂Φ_n/∂(r²) ∇∇(r²) + ! where ∇∇(r²) = 2I + ddphidrdr = ddphidr2dr2 * 4*(R - self%R0)*(R - self%R0) + dphidr2*2 + ddphidzdz = ddphidr2dr2 * 4*(z - self%z0)*(z - self%z0) + dphidr2*2 + ddphidrdz = ddphidr2dr2 * 4*(R - self%R0)*(z - self%z0) + + ! ψ_n = ρ_p(ρ_t)² + if (present(psi_n)) psi_n = rho_p**2 + + ! Using the definitions in `frhotor`: + ! + ! ∇ψ_n = ∂ψ_n/∂Φ_n ∇Φ_n + ! + ! ∂ψ_n/∂Φ_n = Φ_a/ψ_a ∂ψ/∂Φ + ! = Φ_a/ψ_a 1/2πq + ! + ! Using ψ_a = 1/2π Φ_a / (q₀ + Δq), then: + ! + ! ∂ψ_n/∂Φ_n = (q₀ + Δq)/q + ! + q = self%q0 + (self%q1 - self%q0) * rho_p**self%alpha + dq = (self%q1 - self%q0) / (self%alpha/2 + 1) + dpsidphi = (self%q0 + dq) / q + + ! Using the above, ∇ψ_n = ∂ψ_n/∂Φ_n ∇Φ_n + if (present(dpsidr)) dpsidr = dpsidphi * dphidr + if (present(dpsidz)) dpsidz = dpsidphi * dphidz + + ! For the second derivatives: + ! + ! ∇∇ψ_n = ∇(∂ψ_n/∂Φ_n) ∇Φ_n + (∂ψ_n/∂Φ_n) ∇∇Φ_n + ! + ! ∇(∂ψ_n/∂Φ_n) = - (∂ψ_n/∂Φ_n) ∇q/q + ! + ! From q(ψ) = q₀ + (q₁-q₀) ψ_n^α/2, we have: + ! + ! ∇q = α/2 (q-q₀) ∇ψ_n/ψ_n + ! = α/2 (q-q₀)/ψ_n (∂ψ_n/∂Φ_n) ∇Φ_n. + ! + dqdr = self%alpha/2 * (self%q1 - self%q0)*rho_p**(self%alpha-2) * dpsidphi * dphidr + dqdz = self%alpha/2 * (self%q1 - self%q0)*rho_p**(self%alpha-2) * dpsidphi * dphidz + ddpsidphidr = - dpsidphi * dqdr/q + ddpsidphidz = - dpsidphi * dqdz/q + + ! Combining all of the above: + ! + ! ∇∇ψ_n = ∇(∂ψ_n/∂Φ_n) ∇Φ_n + (∂ψ_n/∂Φ_n) ∇∇Φ_n + ! + if (present(ddpsidrr)) ddpsidrr = ddpsidphidr * dphidr + dpsidphi * ddphidrdr + if (present(ddpsidzz)) ddpsidzz = ddpsidphidz * dphidz + dpsidphi * ddphidzdz + if (present(ddpsidrz)) ddpsidrz = ddpsidphidr * dphidz + dpsidphi * ddphidrdz + end subroutine analytic_pol_flux + + + subroutine analytic_pol_curr(self, psi_n, fpol, dfpol) + ! subroutine arguments + class(analytic_equil), intent(in) :: self + real(wp_), intent(in) :: psi_n + real(wp_), intent(out) :: fpol ! poloidal current + real(wp_), intent(out), optional :: dfpol ! its derivative + + unused(psi_n) + + ! The poloidal current function F(ψ) is just a constant: + ! + ! B_φ = B₀R₀/R φ^ ≡ F(ψ)∇φ ⇒ F(ψ)=B₀R₀ + ! + fpol = self%B0 * self%R0 + if (present(dfpol)) dfpol = 0 + end subroutine analytic_pol_curr + + + function analytic_safety(self, psi_n) result(q) + ! function arguments + class(analytic_equil), intent(in) :: self + real(wp_), intent(in) :: psi_n + real(wp_) :: q + + ! local variables + real(wp_) :: rho_p + + ! The safety factor is a power law in ρ_p: + ! + ! q(ρ_p) = q₀ + (q₁-q₀) ρ_p^α + ! + rho_p = sqrt(psi_n) + q = abs(self%q0 + (self%q1 - self%q0) * rho_p**self%alpha) + end function analytic_safety + + + function analytic_pol2tor(self, rho_in) result(rho_out) + ! function arguments + class(analytic_equil), intent(in) :: self + real(wp_), intent(in) :: rho_in + real(wp_) :: rho_out + + ! local variables + real(wp_) :: dq + + ! The change of variable is obtained by integrating + ! + ! q(ψ) = 1/2π ∂Φ/∂ψ + ! + ! and defining ψ = ψ_a ρ_p², Φ = Φ_a ρ_t². + ! The result is: + ! + ! - ψ_a = 1/2π Φ_a / [q₀ + Δq] + ! + ! - ρ_t = ρ_p √[(q₀ + Δq ρ_p^α)/(q₀ + Δq)] + ! + ! where Δq = (q₁ - q₀)/(α/2 + 1) + dq = (self%q1 - self%q0) / (self%alpha/2 + 1) + rho_out = rho_in * sqrt((self%q0 + dq*rho_in**self%alpha) / (self%q0 + dq)) + end function analytic_pol2tor + + + function analytic_tor2pol(self, rho_in) result(rho_out) + ! Converts from toroidal (ρ_t) to poloidal (ρ_p) normalised radius + use const_and_precisions, only : comp_eps + use minpack, only : hybrj1 + + ! function arguments + class(analytic_equil), intent(in) :: self + real(wp_), intent(in) :: rho_in + real(wp_) :: rho_out + + ! local variables + real(wp_) :: rho_p(1), fvec(1), fjac(1,1), wa(7) + integer :: info + + ! In general there is no closed form for ρ_p(ρ_t) in the + ! analytical model, we thus solve numerically the equation + ! ρ_t(ρ_p) = ρ_t₀ for ρ_p. + rho_p = [rho_in] ! first guess, ρ_p ≈ ρ_t + call hybrj1(equation, n=1, x=rho_p, fvec=fvec, fjac=fjac, & + ldfjac=1, tol=comp_eps, info=info, wa=wa, lwa=7) + rho_out = rho_p(1) + + contains + + subroutine equation(n, x, f, df, ldf, flag) + ! The equation to solve: f(x) = ρ_t(x) - ρ_t₀ = 0 + + ! optimal step size + real(wp_), parameter :: e = comp_eps**(1/3.0_wp_) + + ! subroutine arguments + integer, intent(in) :: n, ldf, flag + real(wp_), intent(in) :: x(n) + real(wp_), intent(inout) :: f(n), df(ldf,n) + + if (flag == 1) then + ! return f(x) + f(1) = self%pol2tor(x(1)) - rho_in + else + ! return f'(x), computed numerically + if (x(1) - e > 0) then + df(1,1) = (self%pol2tor(x(1) + e) - self%pol2tor(x(1) - e)) / (2*e) + else + ! single-sided when close to ρ=0 + df(1,1) = (self%pol2tor(x(1) + e) - self%pol2tor(x(1))) / e + end if + end if + end subroutine + + end function analytic_tor2pol + + + subroutine analytic_flux_contour(self, psi0, R_min, R, z, & + R_hi, z_hi, R_lo, z_lo) + use const_and_precisions, only : pi + + ! subroutine arguments + class(analytic_equil), intent(in) :: self + real(wp_), intent(in) :: psi0 + real(wp_), intent(in) :: R_min + real(wp_), intent(out) :: R(:), z(:) + real(wp_), intent(inout) :: R_hi, z_hi, R_lo, z_lo + + ! local variables + integer :: n, i + real(wp_) :: r_g ! geometric minor radius + real(wp_) :: theta ! geometric poloidal angle + + unused(R_min) + unused(R_hi); unused(z_hi); unused(R_lo); unused(z_lo) + + n = size(R) + r_g = sqrt(psi0) * self%a + theta = 2*pi / (n - 1) + + do concurrent (i=1:n) + R(i) = self%R0 + r_g * cos(theta*(i-1)) + z(i) = self%z0 + r_g * sin(theta*(i-1)) + end do + + end subroutine analytic_flux_contour + + + ! + ! Numerical equilibrium + ! + + subroutine numeric_pol_flux(self, R, z, psi_n, dpsidr, dpsidz, & + ddpsidrr, ddpsidzz, ddpsidrz) + ! function arguments + class(numeric_equil), intent(in) :: self + real(wp_), intent(in) :: R, z + real(wp_), intent(out), optional :: & + psi_n, dpsidr, dpsidz, ddpsidrr, ddpsidzz, ddpsidrz + + if (self%psi_domain%contains(R, z)) then + ! Within the safety domain + if (present(psi_n)) psi_n = self%psi_spline%eval(R, z) + if (present(dpsidr)) dpsidr = self%psi_spline%deriv(R, z, 1, 0) + if (present(dpsidz)) dpsidz = self%psi_spline%deriv(R, z, 0, 1) + if (present(ddpsidrr)) ddpsidrr = self%psi_spline%deriv(R, z, 2, 0) + if (present(ddpsidzz)) ddpsidzz = self%psi_spline%deriv(R, z, 0, 2) + if (present(ddpsidrz)) ddpsidrz = self%psi_spline%deriv(R, z, 1, 1) + else + ! Values for outside the domain + if (present(psi_n)) psi_n = -1 + if (present(dpsidr)) dpsidr = 0 + if (present(dpsidz)) dpsidz = 0 + if (present(ddpsidrr)) ddpsidrr = 0 + if (present(ddpsidzz)) ddpsidzz = 0 + if (present(ddpsidrz)) ddpsidrz = 0 + end if + + end subroutine numeric_pol_flux + + + subroutine numeric_pol_curr(self, psi_n, fpol, dfpol) + ! subroutine arguments + class(numeric_equil), intent(in) :: self + real(wp_), intent(in) :: psi_n + real(wp_), intent(out) :: fpol ! poloidal current + real(wp_), intent(out), optional :: dfpol ! its derivative + + if (psi_n <= 1 .and. psi_n >= 0) then + ! Inside plasma + fpol = self%fpol_spline%eval(psi_n) + if (present(dfpol)) dfpol = self%fpol_spline%deriv(psi_n) + else + ! Outside plasma + fpol = self%fpol_a + if (present(dfpol)) dfpol = 0 + end if + end subroutine numeric_pol_curr + + + function numeric_safety(self, psi_n) result(q) + ! function arguments + class(numeric_equil), intent(in) :: self + real(wp_), intent(in) :: psi_n + real(wp_) :: q + + if (psi_n < 1) then + ! Inside plasma + q = self%q_spline%eval(psi_n) + else + ! Outside plasma, q is undefined + q = 0 + end if + end function numeric_safety + + + function numeric_pol2tor(self, rho_in) result(rho_out) + ! function arguments + class(numeric_equil), intent(in) :: self + real(wp_), intent(in) :: rho_in + real(wp_) :: rho_out + + rho_out = self%rhot_spline%eval(rho_in) + end function numeric_pol2tor + + + function numeric_tor2pol(self, rho_in) result(rho_out) + ! function arguments + class(numeric_equil), intent(in) :: self + real(wp_), intent(in) :: rho_in + real(wp_) :: rho_out + + rho_out = self%rhop_spline%eval(rho_in) + end function numeric_tor2pol + + + subroutine numeric_flux_contour(self, psi0, R_min, R, z, & + R_hi, z_hi, R_lo, z_lo) + use const_and_precisions, only : pi + use logger, only : log_warning + use dierckx, only : profil, sproota + + ! subroutine arguments + class(numeric_equil), intent(in) :: self + real(wp_), intent(in) :: psi0 + real(wp_), intent(in) :: R_min + real(wp_), intent(out) :: R(:), z(:) + real(wp_), intent(inout) :: R_hi, z_hi, R_lo, z_lo + + ! local variables + integer :: n, np, i + integer :: ier, iopt, m + real(wp_) :: theta + real(wp_) :: R_hi1, z_hi1, R_lo1, z_lo1, zc + real(wp_) :: czc(self%psi_spline%nknots_x), zeroc(4) + character(256) :: msg + + n = size(R) + np = (n - 1)/2 + theta = pi / np + + call self%find_htg_point(R_hi, z_hi, R_hi1, z_hi1, psi0) + call self%find_htg_point(R_lo, z_lo, R_lo1, z_lo1, psi0) + + R(1) = R_lo1 + z(1) = z_lo1 + R(n) = R_lo1 + z(n) = z_lo1 + R(np+1) = R_hi1 + z(np+1) = z_hi1 + + do i = 2, np + zc = z_lo1 + (z_hi1 - z_lo1) * (1 - cos(theta*(i-1))) / 2 + iopt = 1 + associate (s => self%psi_spline) + call profil(iopt, s%knots_x, s%nknots_x, s%knots_y, s%nknots_y, & + s%coeffs, 3, 3, zc, s%nknots_x, czc, ier) + if (ier > 0) then + write(msg, '(a, a, g0)') & + 'when computing ψ(R,z) contour `profil` returned ier=', ier + call log_warning(msg, mod='gray_equil', proc='numeric_flux_contour') + end if + call sproota(psi0, s%knots_x, s%nknots_x, czc, zeroc, 4, m, ier) + end associate + + if (zeroc(1) > R_min) then + R(i) = zeroc(1) + R(n+1-i) = zeroc(2) + else + R(i) = zeroc(2) + R(n+1-i) = zeroc(3) + end if + z(i) = zc + z(n+1-i) = zc + end do + + ! Replace the initial guess with the exact values + R_hi = R_hi1 + z_hi = z_hi1 + R_lo = R_lo1 + z_lo = z_lo1 + end subroutine numeric_flux_contour + + + subroutine numeric_init(self, params, data, err) + ! Initialises a numeric equilibrium + use const_and_precisions, only : zero, one + use gray_params, only : equilibrium_parameters + use gray_params, only : EQ_EQDSK_FULL, EQ_EQDSK_PARTIAL + use gray_params, only : X_AT_TOP, X_AT_BOTTOM, X_IS_MISSING + use utils, only : vmaxmin, vmaxmini + use logger, only : log_debug, log_info + + ! subroutine arguments + class(numeric_equil), intent(out) :: self + type(equilibrium_parameters), intent(in) :: params + type(eqdsk_data), intent(in) :: data + integer, intent(out) :: err + + ! local variables + integer :: nr, nz, nrz, npsi, nbnd, ibinf, ibsup + real(wp_) :: psinoptmp, psinxptmp + real(wp_) :: rbinf, rbsup, zbinf, zbsup, R1, z1 + real(wp_) :: psiant, psinop + real(wp_), dimension(:), allocatable :: rv1d, zv1d, fvpsi + integer :: i, j, ij + character(256) :: msg ! for log messages formatting + + ! compute array sizes + nr = size(data%grid_r) + nz = size(data%grid_z) + nrz = nr*nz + npsi = size(data%psi) + + ! length in m !!! + self%r_range = [data%grid_r(1), data%grid_r(nr)] + self%z_range = [data%grid_z(1), data%grid_z(nz)] + + call log_debug('computing splines...', mod='gray_equil', proc='numeric_init') + + ! Spline interpolation of ψ(R, z) + + select case (params%iequil) + + case (EQ_EQDSK_PARTIAL) + ! Data valid only inside boundary (data%psi=0 outside), + ! presence of boundary anticipated here to filter invalid data + nbnd = size(data%boundary%R) + + ! allocate knots and spline coefficients arrays + allocate(self%psi_spline%knots_x(nr + 4), self%psi_spline%knots_y(nz + 4)) + allocate(self%psi_spline%coeffs(nrz)) + + ! determine number of valid grid points + nrz=0 + do j=1,nz + do i=1,nr + if (nbnd > 0) then + if (.not. data%boundary%contains(data%grid_r(i), data%grid_z(j))) cycle + else + if (data%psi_map(i,j) <= 0) cycle + end if + nrz=nrz+1 + end do + end do + + ! store valid data + allocate(rv1d(nrz), zv1d(nrz), fvpsi(nrz)) + ij=0 + do j=1,nz + do i=1,nr + if (nbnd > 0) then + if (.not. data%boundary%contains(data%grid_r(i), data%grid_z(j))) cycle + else + if (data%psi_map(i,j) <= 0) cycle + end if + ij = ij + 1 + rv1d(ij) = data%grid_r(i) + zv1d(ij) = data%grid_z(j) + fvpsi(ij) = data%psi_map(i,j) + end do + end do + + ! Fit as a scattered set of points + ! use reduced number of knots to limit memory comsumption ? + self%psi_spline%nknots_x=nr/4+4 + self%psi_spline%nknots_y=nz/4+4 + call self%psi_spline%init_nonreg( & + rv1d, zv1d, fvpsi, tension=params%ssplps, & + range=[self%r_range, self%z_range], err=err) + + ! if failed, re-fit with an interpolating spline (zero tension) + if (err == -1) then + err = 0 + self%psi_spline%nknots_x=nr/4+4 + self%psi_spline%nknots_y=nz/4+4 + call self%psi_spline%init_nonreg( & + rv1d, zv1d, fvpsi, tension=zero, & + range=[self%r_range, self%z_range], err=err) + end if + deallocate(rv1d, zv1d, fvpsi) + ! reset nrz to the total number of grid points for next allocations + nrz = nr*nz + + case (EQ_EQDSK_FULL) + ! Data are valid on the full R,z grid + + ! reshape 2D ψ array to 1D (transposed) + allocate(fvpsi(nrz)) + fvpsi = reshape(transpose(data%psi_map), [nrz]) + + ! compute spline coefficients + call self%psi_spline%init(data%grid_r, data%grid_z, fvpsi, nr, nz, & + range=[self%r_range, self%z_range], & + tension=params%ssplps, err=err) + + ! if failed, re-fit with an interpolating spline (zero tension) + if (err == -1) then + call self%psi_spline%init(data%grid_r, data%grid_z, fvpsi, nr, nz, & + range=[self%r_range, self%z_range], tension=zero) + err = 0 + end if + deallocate(fvpsi) + end select + + if (err /= 0) then + err = 2 + return + end if + + ! compute spline coefficients for ψ(R,z) partial derivatives + call self%psi_spline%init_deriv(nr, nz, 1, 0) ! ∂ψ/∂R + call self%psi_spline%init_deriv(nr, nz, 0, 1) ! ∂ψ/∂z + call self%psi_spline%init_deriv(nr, nz, 1, 1) ! ∂²ψ/∂R∂z + call self%psi_spline%init_deriv(nr, nz, 2, 0) ! ∂²ψ/∂R² + call self%psi_spline%init_deriv(nr, nz, 0, 2) ! ∂²ψ/∂z² + + ! set the initial ψ(R,z) domain to the grid boundary + ! + ! Note: this is required to bootstrap `flux_pol` calls + ! within this very subroutine. + self%psi_domain = contour(Rmin=self%r_range(1), Rmax=self%r_range(2), & + zmin=self%z_range(1), zmax=self%z_range(2)) + + ! Spline interpolation of F(ψ) + + ! give a small weight to the last point + block + real(wp_) :: w(npsi) + w(1:npsi-1) = 1 + w(npsi) = 1.0e2_wp_ + call self%fpol_spline%init(data%psi, data%fpol, npsi, weights=w, & + range=[zero, one], tension=params%ssplf) + end block + + ! set vacuum value used outside 0 ≤ ψ ≤ 1 range + self%fpol_a = self%fpol_spline%eval(data%psi(npsi)) + self%sgn_bphi = sign(one, self%fpol_a) + + ! Re-normalize ψ_n after the spline computation + ! Note: this ensures 0 ≤ ψ_n'(R,z) < 1 inside the plasma + + ! Start with un-corrected ψ_n + self%psi_a = data%psi_a + psinop = 0 ! ψ_n(O point) + psiant = 1 ! ψ_n(X point) - ψ_n(O point) + + ! Use provided boundary to set an initial guess + ! for the search of O/X points + nbnd = size(data%boundary%R) + if (nbnd > 0) then + call vmaxmini(data%boundary%z, nbnd, zbinf, zbsup, ibinf, ibsup) + rbinf = data%boundary%R(ibinf) + rbsup = data%boundary%R(ibsup) + else + zbinf = data%grid_z(2) + zbsup = data%grid_z(nz-1) + rbinf = data%grid_r((nr+1)/2) + rbsup = rbinf + end if + + ! Search for exact location of the magnetic axis + call self%find_ox_point(R0=data%axis(1), z0=data%axis(2), & + R1=self%axis(1), z1=self%axis(2), psi1=psinoptmp) + + write (msg, '("O-point found:", 3(x,a,"=",g0.3))') & + 'r', self%axis(1), 'z', self%axis(2), 'ψ', psinoptmp + call log_info(msg, mod='gray_equil', proc='numeric_init') + + ! Search for X-point + select case (params%ixp) + case (X_AT_BOTTOM) + call self%find_ox_point(R0=rbinf, z0=zbinf, R1=R1, z1=z1, psi1=psinxptmp) + rbinf = z1 + if (psinxptmp /= -1) then + write (msg, '("X-point found:", 3(x,a,"=",g0.3))') & + 'R', R1, 'z', z1, 'ψ', psinxptmp + call log_info(msg, mod='gray_equil', proc='numeric_init') + psinop = psinoptmp + psiant = psinxptmp-psinop + call self%find_htg_point(R0=self%axis(1), z0=(self%axis(2)+zbsup)/2, & + R1=r1, z1=zbsup, psi0=one) + end if + + case (X_AT_TOP) + call self%find_ox_point(R0=rbsup, z0=zbsup, R1=R1, z1=z1, psi1=psinxptmp) + zbsup = z1 + if (psinxptmp /= -1) then + write (msg, '("X-point found:", 3(x,a,"=",g0.3))') & + 'R', r1, 'z', zbsup, 'ψ', psinxptmp + call log_info(msg, mod='gray_equil', proc='numeric_init') + psinop = psinoptmp + psiant = psinxptmp - psinop + call self%find_htg_point(R0=self%axis(1), z0=(self%axis(2)+zbinf)/2, & + R1=r1, z1=zbinf, psi0=one) + end if + + case (X_IS_MISSING) + psinop = psinoptmp + psiant = 1 - psinop + ! Find upper horizontal tangent point + call self%find_htg_point(R0=self%axis(1), z0=(self%axis(2)+zbsup)/2, & + R1=rbsup, z1=zbsup, psi0=one) + ! Find lower horizontal tangent point + call self%find_htg_point(R0=self%axis(1), z0=(self%axis(2)+zbinf)/2, & + R1=rbinf, z1=zbinf, psi0=one) + write (msg, '("X-point not found in", 2(x,a,"∈[",g0.3,",",g0.3,"]"))') & + 'R', rbinf, rbsup, 'z', zbinf, zbsup + call log_info(msg, mod='gray_equil', proc='numeric_init') + end select + + ! Correct the spline coefficients: ψ_n → (ψ_n - psinop)/psiant + call self%psi_spline%transform(1/psiant, -psinop/psiant) + + ! Do the opposite scaling to preserve un-normalised values + ! Note: this is only used for the poloidal magnetic field + self%psi_a = self%psi_a * psiant + + ! Compute other constants + call self%pol_curr(zero, self%b_axis) + self%b_axis = self%b_axis / self%axis(1) + self%b_centre = self%fpol_a / data%r_ref + self%r_centre = data%r_ref + self%z_boundary = [zbinf, zbsup] + + write (msg, '(2(a,g0.3))') 'B_centre=', self%b_centre, ' B_axis=', self%b_axis + call log_info(msg, mod='gray_equil', proc='numeric_init') + + ! Compute ρ_p/ρ_t mapping based on the input q profile + block + use const_and_precisions, only : pi + real(wp_), dimension(npsi) :: phi, rho_p, rho_t + real(wp_) :: dx + integer :: k + + call self%q_spline%init(data%psi, abs(data%q), npsi) + + ! Toroidal flux as Φ(ψ) = 2π ∫q(ψ)dψ + phi(1) = 0 + do k = 1, npsi-1 + dx = self%q_spline%data(k+1) - self%q_spline%data(k) + phi(k+1) = phi(k) + dx*(self%q_spline%coeffs(k,1) + dx*(self%q_spline%coeffs(k,2)/2 + & + dx*(self%q_spline%coeffs(k,3)/3 + dx* self%q_spline%coeffs(k,4)/4) ) ) + end do + + self%phi_a = phi(npsi) + rho_p = sqrt(data%psi) + rho_t = sqrt(phi/self%phi_a) + self%phi_a = 2*pi * abs(self%psi_a) * self%phi_a + + call self%rhop_spline%init(rho_t, rho_p, size(rho_p)) + call self%rhot_spline%init(rho_p, rho_t, size(rho_t)) + end block + + call log_debug('splines computed', mod='gray_equil', proc='numeric_init') + + ! Compute the domain of the ψ mapping + self%psi_domain = data%boundary + call rescale_boundary(self%psi_domain, self%psi_spline, O=self%axis, t0=1.5_wp_) + + end subroutine numeric_init + + + subroutine rescale_boundary(cont, psi, O, t0) + ! Given the plasma boundary contour `cont`, the position of the + ! magnetic axis `O`, and a scaling factor `t0`; this subroutine + ! rescales the contour by `t0` about `O` while ensuring the + ! psi_spline stays monotonic within the new boundary. + + ! subroutine arguments + type(contour), intent(inout) :: cont ! (R,z) contour + type(spline_2d), intent(inout) :: psi ! ψ(R,z) spline + real(wp_), intent(in) :: O(2) ! center point + real(wp_), intent(in) :: t0 ! scaling factor + + ! subroutine variables + integer :: i + real(wp_) :: t + real(wp_), parameter :: dt = 0.05 + real(wp_) :: P(2), N(2) + + do i = 1, size(cont%R) + ! For each point on the contour compute: + P = [cont%R(i), cont%z(i)] ! point on the contour + N = P - O ! direction of the line from O to P + + ! Find the max t: s(t) = ψ(O + tN) is monotonic in [1, t] + t = 1 + do while (t < t0) + if (s(t + dt) < s(t)) exit + t = t + dt + end do + + ! The final point is thus O + tN + P = O + t * N + cont%R(i) = P(1) + cont%z(i) = P(2) + end do + + contains + + function s(t) + ! Rescriction of ψ(R, z) on the line Q(t) = O + tN + + real(wp_), intent(in) :: t + real(wp_) :: s, Q(2) + + Q = O + t * N + s = psi%eval(Q(1), Q(2)) + end function + + end subroutine rescale_boundary + + + subroutine find_ox_point(self, R0, z0, R1, z1, psi1) + ! Given the point (R₀,z₀) as an initial guess, finds + ! the exact location (R₁,z₁) where ∇ψ(R₁,z₁) = 0. + ! It also returns ψ₁=ψ(R₁,z₁). + ! + ! Note: this is used to find the magnetic X and O point, + ! because both are stationary points for ψ(R,z). + use const_and_precisions, only : comp_eps + use minpack, only : hybrj1 + use logger, only : log_error, log_debug + + ! subroutine arguments + class(numeric_equil) :: self + real(wp_), intent(in) :: R0, z0 + real(wp_), intent(out) :: R1, z1, psi1 + + ! local variables + integer :: info + real(wp_) :: sol(2), f(2), df(2,2), wa(15) + character(256) :: msg + + sol = [R0, z0] ! first guess + + call hybrj1(equation, n=2, x=sol, fvec=f, fjac=df, ldfjac=2, & + tol=sqrt(comp_eps), info=info, wa=wa, lwa=15) + if (info /= 1) then + write (msg, '("guess: ", g0.3, ", ", g0.3)') R0, z0 + call log_debug(msg, mod='equilibrium', proc='find_ox_point') + write (msg, '("solution: ", g0.3, ", ", g0.3)') sol + call log_debug(msg, mod='equilibrium', proc='find_ox_point') + write (msg, '("hybrj1 failed with error ",g0)') info + call log_error(msg, mod='equilibrium', proc='find_ox_point') + end if + + R1 = sol(1) ! solution + z1 = sol(2) + call self%pol_flux(R1, z1, psi1) + + contains + + subroutine equation(n, x, f, df, ldf, flag) + ! The equation to solve: f(R,z) = ∇ψ(R,z) = 0 + + ! subroutine arguments + integer, intent(in) :: n, flag, ldf + real(wp_), intent(in) :: x(n) + real(wp_), intent(inout) :: f(n), df(ldf,n) + + if (flag == 1) then + ! return f(R,z) = ∇ψ(R,z) + call self%pol_flux(R=x(1), z=x(2), dpsidr=f(1), dpsidz=f(2)) + else + ! return ∇f(R,z) = ∇∇ψ(R,z) + call self%pol_flux(R=x(1), z=x(2), ddpsidrr=df(1,1), & + ddpsidzz=df(2,2), ddpsidrz=df(1,2)) + df(2,1) = df(1,2) + end if + end subroutine equation + + end subroutine find_ox_point + + + subroutine find_htg_point(self, R0, z0, R1, z1, psi0) + ! Given the point (R₀,z₀) as an initial guess, finds + ! the exact location (R₁,z₁) where: + ! { ψ(R₁,z₁) = ψ₀ + ! { ∂ψ/∂R(R₁,z₁) = 0 . + ! + ! Note: this is used to find the horizontal tangent + ! point of the contour ψ(R,z)=ψ₀. + use const_and_precisions, only : comp_eps + use minpack, only : hybrj1 + use logger, only : log_error, log_debug + + ! subroutine arguments + class(numeric_equil) :: self + real(wp_), intent(in) :: R0, z0, psi0 + real(wp_), intent(out) :: R1, z1 + + ! local variables + integer :: info + real(wp_) :: sol(2), f(2), df(2,2), wa(15) + character(256) :: msg + + sol = [R0, z0] ! first guess + + call hybrj1(equation, n=2, x=sol, fvec=f, fjac=df, ldfjac=2, & + tol=sqrt(comp_eps), info=info, wa=wa, lwa=15) + if (info /= 1) then + write (msg, '("guess: ", g0.3, ", ", g0.3)') R0, z0 + call log_debug(msg, mod='equilibrium', proc='find_ox_point') + write (msg, '("solution: ", g0.3, ", ", g0.3)') sol + call log_debug(msg, mod='equilibrium', proc='find_htg_point') + write (msg, '("hybrj1 failed with error ",g0)') info + call log_error(msg, mod='equilibrium', proc='find_htg_point') + end if + + R1 = sol(1) ! solution + z1 = sol(2) + contains + + subroutine equation(n, x, f, df, ldf, flag) + ! The equation to solve: f(R,z) = [ψ(R,z)-ψ₀, ∂ψ/∂R] = 0 + + ! subroutine arguments + integer, intent(in) :: n, ldf, flag + real(wp_), intent(in) :: x(n) + real(wp_), intent(inout) :: f(n), df(ldf, n) + + if (flag == 1) then + ! return f(R,z) = [ψ(R,z)-ψ₀, ∂ψ/∂R] + call self%pol_flux(R=x(1), z=x(2), psi_n=f(1), dpsidr=f(2)) + f(1) = f(1) - psi0 + else + ! return ∇f(R,z) = [[∂ψ/∂R, ∂ψ/∂z], [∂²ψ/∂R², ∂²ψ/∂R∂z]] + call self%pol_flux(R=x(1), z=x(2), dpsidr=df(1,1), dpsidz=df(1,2), & + ddpsidrr=df(2,1), ddpsidrz=df(2,2)) + end if + end subroutine equation + + end subroutine find_htg_point + + + ! + ! Vacuum + ! + + subroutine vacuum_pol_flux(self, R, z, psi_n, dpsidr, dpsidz, & + ddpsidrr, ddpsidzz, ddpsidrz) + ! function arguments + class(vacuum), intent(in) :: self + real(wp_), intent(in) :: R, z + real(wp_), intent(out), optional :: & + psi_n, dpsidr, dpsidz, ddpsidrr, ddpsidzz, ddpsidrz + + unused(self); unused(R); unused(z) + + ! ψ(R,z) is undefined everywhere, return -1 + if (present(psi_n)) psi_n = -1 + if (present(dpsidr)) dpsidr = 0 + if (present(dpsidz)) dpsidz = 0 + if (present(ddpsidrr)) ddpsidrr = 0 + if (present(ddpsidzz)) ddpsidzz = 0 + if (present(ddpsidrz)) ddpsidrz = 0 + end subroutine vacuum_pol_flux + + + subroutine vacuum_pol_curr(self, psi_n, fpol, dfpol) + ! subroutine arguments + class(vacuum), intent(in) :: self + real(wp_), intent(in) :: psi_n + real(wp_), intent(out) :: fpol ! poloidal current + real(wp_), intent(out), optional :: dfpol ! its derivative + + unused(self); unused(psi_n) + + ! There is no current, F(ψ)=0 + fpol = 0 + if (present(dfpol)) dfpol = 0 + end subroutine vacuum_pol_curr + + + function vacuum_safety(self, psi_n) result(q) + ! function arguments + class(vacuum), intent(in) :: self + real(wp_), intent(in) :: psi_n + real(wp_) :: q + + unused(self); unused(psi_n) + + ! q(ψ) is undefined, return 0 + q = 0 + end function vacuum_safety + + + function vacuum_conv(self, rho_in) result(rho_out) + ! function arguments + class(vacuum), intent(in) :: self + real(wp_), intent(in) :: rho_in + real(wp_) :: rho_out + + unused(self) + + ! Neither ρ_p nor ρ_t are defined, do nothing + rho_out = rho_in + end function vacuum_conv + + + subroutine vacuum_flux_contour(self, psi0, R_min, R, z, & + R_hi, z_hi, R_lo, z_lo) + ! subroutine arguments + class(vacuum), intent(in) :: self + real(wp_), intent(in) :: psi0 + real(wp_), intent(in) :: R_min + real(wp_), intent(out) :: R(:), z(:) + real(wp_), intent(inout) :: R_hi, z_hi, R_lo, z_lo + + unused(self); unused(psi0); unused(R_min) + unused(R_hi); unused(z_hi); unused(R_lo); unused(z_lo) + + ! flux surfaces are undefined + R = 0 + z = 0 + end subroutine vacuum_flux_contour + + + ! + ! Helpers + ! + + subroutine load_equil(params, equil, limiter, err) + ! Loads a generic MHD equilibrium and limiter + ! contour from file (params%filenm) + use gray_params, only : equilibrium_parameters, EQ_VACUUM + use gray_params, only : EQ_ANALYTICAL, EQ_EQDSK_FULL, EQ_EQDSK_PARTIAL + use types, only : contour + use logger, only : log_warning, log_debug + + ! subroutine arguments + type(equilibrium_parameters), intent(inout) :: params + class(abstract_equil), allocatable, intent(out) :: equil + type(contour), intent(out) :: limiter + integer, intent(out) :: err + + ! local variables + type(numeric_equil) :: ne + type(analytic_equil) :: ae + type(eqdsk_data) :: data + + select case (params%iequil) + case (EQ_VACUUM) + call log_debug('vacuum, doing nothing', mod='gray_equil', proc='load_equil') + allocate(equil, source=vacuum()) + return + + case (EQ_ANALYTICAL) + call log_debug('loading analytical file', mod='gray_equil', proc='load_equil') + call load_analytical(params, ae, limiter, err) + if (err /= 0) return + allocate(equil, source=ae) + + case (EQ_EQDSK_FULL, EQ_EQDSK_PARTIAL) + call log_debug('loading G-EQDSK file', mod='gray_equil', proc='load_equil') + call load_eqdsk(params, data, err) + if (err /= 0) return + + call change_cocos(data, params%icocos, 3) + call scale_eqdsk(params, data) + + allocate(equil, mold=ne) + select type (equil) + type is (numeric_equil) + call equil%init(params, data, err) + if (err /= 0) return + end select + + limiter = data%limiter + + end select + + end subroutine load_equil + + + subroutine load_analytical(params, equil, limiter, err) + ! Loads the analytical equilibrium and limiter contour from file + use const_and_precisions, only : pi, one + use gray_params, only : equilibrium_parameters + use logger, only : log_error + + ! subroutine arguments + type(equilibrium_parameters), intent(in) :: params + type(analytic_equil), intent(out) :: equil + type(contour), intent(out) :: limiter + integer, intent(out) :: err + + ! local variables + integer :: i, u, nlim + real(wp_), allocatable :: R(:), z(:) + + open(file=params%filenm, status='old', action='read', newunit=u, iostat=err) + if (err /= 0) then + call log_error('opening equilibrium ('//trim(params%filenm)//') failed!', & + mod='gray_equil', proc='load_analytical') + err = 1 + return + end if + + ! The analytical file format is: + ! + ! 1 R0 z0 a + ! 2 B0 + ! 3 q0 q1 alpha + ! 4 nlim + ! 5 R z + ! ... + + ! Read model parameters + read (u, *) equil%R0, equil%z0, equil%a + read (u, *) equil%B0 + read (u, *) equil%q0, equil%q1, equil%alpha + read (u, *) nlim + + ! Read limiter points + if (nlim > 0) then + allocate(R(nlim), z(nlim)) + read (u, *) (R(i), z(i), i = 1, nlim) + limiter = contour(R, z) + end if + close(u) + + ! Notes on cocos=3 + ! 1. With both I_p,B_φ > 0 we have ∂ψ/∂r<0 and ∂Φ/∂r>0. + ! 2. ψ_a = ψ_edge - ψ_axis < 0. + ! 3. q = 1/2π ∂Φ/∂ψ ~ ∂Φ/∂r⋅∂r/∂ψ < 0. + ! 4. In general, sgn(q) = -sgn(I_p)⋅sgn(B_φ). + + ! Apply forced signs + if (params%sgni /= 0) then + equil%q0 = sign(equil%q0, -params%sgni*one) + equil%q1 = sign(equil%q1, -params%sgni*one) + end if + if (params%sgnb /= 0) then + equil%B0 = sign(equil%B0, +params%sgnb*one) + end if + + ! Apply rescaling factors + equil%B0 = equil%B0 * params%factb + + ! Toroidal flux at r=a: + ! + ! Φ(a) = B₀πa² 2γ/(γ + 1) + ! + ! where γ=1/√(1-ε²), + ! ε=a/R₀ is the tokamak aspect ratio + block + real(wp_) :: gamma + gamma = 1/sqrt(1 - (equil%a/equil%R0)**2) + equil%phi_a = equil%B0 * pi * equil%a**2 * 2*gamma/(gamma + 1) + end block + + ! In cocos=3 the safety factor is + ! + ! q(ψ) = 1/2π ∂Φ/∂ψ. + ! + ! Given the power law of the model + ! + ! q(ψ) = q₀ + (q₁-q₀) (ψ/ψa)^(α/2), + ! + ! we can find ψ_a = ψ(r=a) by integrating: + ! + ! ∫ q(ψ)dψ = 1/2π ∫ dΦ + ! ∫₀^ψ_a q(ψ)dψ = 1/2π Φ_a + ! ψa [q₀ + (q₁-q₀)/(α/2+1)] = Φa/2π + ! + ! ⇒ ψ_a = Φ_a 1/2π 1/(q₀ + Δq) + ! + ! where Δq = (q₁ - q₀)/(α/2 + 1) + block + real(wp_) :: dq + dq = (equil%q1 - equil%q0) / (equil%alpha/2 + 1) + equil%psi_a = 1/(2*pi) * equil%phi_a / (equil%q0 + dq) + end block + + ! Compute other constants (see abstract_equil) + equil%b_axis = equil%B0 + equil%b_centre = equil%B0 + equil%r_centre = equil%R0 + equil%sgn_bphi = sign(one, equil%B0) + equil%axis = [equil%R0, equil%z0] + equil%r_range = equil%R0 + [-equil%a, equil%a] + equil%z_range = equil%z0 + [-equil%a, equil%a] + equil%z_boundary = equil%z0 + [-equil%a, equil%a] + end subroutine load_analytical + + + subroutine load_eqdsk(params, data, err) + ! Loads the equilibrium `data` from a G-EQDSK file (params%filenm). + ! For a description of the G-EQDSK format, see the GRAY user manual. + use const_and_precisions, only : one + use gray_params, only : equilibrium_parameters + use utils, only : get_free_unit + use logger, only : log_error + + ! subroutine arguments + type(equilibrium_parameters), intent(in) :: params + type(eqdsk_data), intent(out) :: data + integer, intent(out) :: err + + ! local variables + integer :: u, i, j, nr, nz, nlim, nbnd + character(len=48) :: string + real(wp_) :: r_dim, z_dim, r_left, z_mid, psi_edge, psi_axis + real(wp_) :: skip_r ! dummy variables, used to discard data + integer :: skip_i ! + + ! Open the G-EQDSK file + open(file=params%filenm, status='old', action='read', newunit=u, iostat=err) + if (err /= 0) then + call log_error('opening eqdsk file ('//trim(params%filenm)//') failed!', & + mod='gray_equil', proc='load_eqdsk') + err = 1 + return + end if + + ! Get size of main arrays and allocate them + if (params%idesc) then + read (u,'(a48,3i4)') string, skip_i, nr, nz + else + read (u,*) nr, nz + end if + allocate(data%grid_r(nr), data%grid_z(nz), data%psi_map(nr, nz)) + allocate(data%psi(nr), data%fpol(nr), data%q(nr)) + + ! Store 0D data and main arrays + if (params%ifreefmt) then + read (u, *) r_dim, z_dim, data%r_ref, r_left, z_mid + read (u, *) data%axis, psi_axis, psi_edge, skip_r + read (u, *) skip_r, skip_r, skip_r, skip_r, skip_r + read (u, *) skip_r, skip_r, skip_r, skip_r, skip_r + read (u, *) (data%fpol(i), i=1,nr) + read (u, *) (skip_r,i=1, nr) + read (u, *) (skip_r,i=1, nr) + read (u, *) (skip_r,i=1, nr) + read (u, *) ((data%psi_map(i,j), i=1,nr), j=1,nz) + read (u, *) (data%q(i), i=1,nr) + else + read (u, '(5e16.9)') r_dim, z_dim, data%r_ref, r_left, z_mid + read (u, '(5e16.9)') data%axis, psi_axis, psi_edge, skip_r + read (u, '(5e16.9)') skip_r, skip_r, skip_r, skip_r, skip_r + read (u, '(5e16.9)') skip_r, skip_r, skip_r, skip_r, skip_r + read (u, '(5e16.9)') (data%fpol(i), i=1,nr) + read (u, '(5e16.9)') (skip_r, i=1,nr) + read (u, '(5e16.9)') (skip_r, i=1,nr) + read (u, '(5e16.9)') (skip_r, i=1,nr) + read (u, '(5e16.9)') ((data%psi_map(i,j), i=1,nr), j=1,nz) + read (u, '(5e16.9)') (data%q(i), i=1,nr) + end if + + ! Get size of boundary and limiter arrays and allocate them + read (u, *) nbnd, nlim + + ! Load plasma boundary data + if (nbnd > 0) then + block + real(wp_) :: R(nbnd), z(nbnd) + if (params%ifreefmt) then + read(u, *) (R(i), z(i), i=1,nbnd) + else + read(u, '(5e16.9)') (R(i), z(i), i=1,nbnd) + end if + data%boundary = contour(R, z) + end block + end if + + ! Load limiter data + if (nlim > 0) then + block + real(wp_) :: R(nlim), z(nlim) + if (params%ifreefmt) then + read(u, *) (R(i), z(i), i=1,nlim) + else + read(u, '(5e16.9)') (R(i), z(i), i=1,nlim) + end if + data%limiter = contour(R, z) + end block + end if + + ! End of G-EQDSK file + close(u) + + ! Build the grid arrays + block + real(wp_) :: dpsi, dr, dz + dr = r_dim/(nr - 1) + dz = z_dim/(nz - 1) + dpsi = one/(nr - 1) + do i = 1, nr + data%psi(i) = (i-1)*dpsi + data%grid_r(i) = r_left + (i-1)*dr + end do + do i = 1, nz + data%grid_z(i) = z_mid - z_dim/2 + (i-1)*dz + end do + end block + + ! Normalize the poloidal flux + data%psi_a = psi_edge - psi_axis + if (.not. params%ipsinorm) data%psi_map = (data%psi_map - psi_axis)/data%psi_a + + end subroutine load_eqdsk + + + subroutine scale_eqdsk(params, data) + ! Rescale the magnetic field (B) and the plasma current (I_p) + ! and/or force their signs. + ! + ! Notes: + ! 1. signi and signb are ignored on input if equal to 0. + ! They are used to assign the direction of B_φ and I_p BEFORE scaling. + ! 2. cocos=3 assumed: positive toroidal direction is CCW from above + ! 3. B_φ and I_p scaled by the same factor factb to keep q unchanged + ! 4. factb<0 reverses the directions of Bphi and Ipla + use const_and_precisions, only : one + use gray_params, only : equilibrium_parameters + + ! subroutine arguments + type(equilibrium_parameters), intent(inout) :: params + type(eqdsk_data), intent(inout) :: data + + ! Notes on cocos=3 + ! 1. With both I_p,B_φ > 0 we have ∂ψ/∂r<0 and ∂Φ/∂r>0. + ! 2. ψ_a = ψ_edge - ψ_axis < 0. + ! 3. q = 1/2π ∂Φ/∂ψ ~ ∂Φ/∂r⋅∂r/∂ψ < 0. + ! 4. In general, sgn(q) = -sgn(I_p)⋅sgn(B_φ). + + ! Apply signs + if (params%sgni /= 0) data%psi_a = sign(data%psi_a, -params%sgni*one) + if (params%sgnb /= 0) data%fpol = sign(data%fpol, +params%sgnb*one) + + ! Rescale + data%psi_a = data%psi_a * params%factb + data%fpol = data%fpol * params%factb + + ! Compute the signs to be shown in the outputs header + ! Note: this is needed if sgni, sgnb = 0 in gray.ini + params%sgni = int(sign(one, -data%psi_a)) + params%sgnb = int(sign(one, +data%fpol(size(data%fpol)))) + end subroutine scale_eqdsk + + + subroutine change_cocos(data, cocosin, cocosout, err) + ! Convert the MHD equilibrium data from one coordinate convention + ! (COCOS) to another. These are specified by `cocosin` and + ! `cocosout`, respectively. + ! + ! For more information, see: https://doi.org/10.1016/j.cpc.2012.09.010 + use const_and_precisions, only : zero, one, pi + use logger, only : log_warning + + ! subroutine arguments + type(eqdsk_data), intent(inout) :: data + integer, intent(in) :: cocosin, cocosout + integer, intent(out), optional :: err + + ! local variables + real(wp_) :: isign, bsign + integer :: exp2pi, exp2piout + logical :: phiccw, psiincr, qpos, phiccwout, psiincrout, qposout + + call decode_cocos(cocosin, exp2pi, phiccw, psiincr, qpos) + call decode_cocos(cocosout, exp2piout, phiccwout, psiincrout, qposout) + + ! Check sign consistency + isign = sign(one, data%psi_a) + if (.not. psiincr) isign = -isign + bsign = sign(one, data%fpol(size(data%fpol))) + if (qpos .neqv. isign * bsign * data%q(size(data%q)) > zero) then + ! Warning: sign inconsistency found among q, I_p and B_ref + data%q = -data%q + call log_warning('data not consistent with cocosin', & + mod='gray_equil', proc='change_cocos') + if (present(err)) err = 1 + else + if (present(err)) err = 0 + end if + + ! Convert cocosin to cocosout + + ! Opposite direction of toroidal angle phi in cocosin and cocosout + if (phiccw .neqv. phiccwout) data%fpol = -data%fpol + + ! q has opposite sign for given sign of B_phi⋅I_p + if (qpos .neqv. qposout) data%q = -data%q + + ! psi and Ip signs don't change accordingly + if ((phiccw .eqv. phiccwout) .neqv. (psiincr .eqv. psiincrout)) & + data%psi_a = -data%psi_a + + ! Convert Wb to Wb/rad or viceversa + data%psi_a = data%psi_a * (2*pi)**(exp2piout - exp2pi) + end subroutine change_cocos + + + subroutine decode_cocos(cocos, exp2pi, phiccw, psiincr, qpos) + ! Extracts the sign and units conventions from a COCOS index + + ! subroutine arguments + integer, intent(in) :: cocos + integer, intent(out) :: exp2pi + logical, intent(out) :: phiccw, psiincr, qpos + + ! local variables + integer :: cmod10, cmod4 + + cmod10 = mod(cocos, 10) + cmod4 = mod(cmod10, 4) + + ! cocos>10 ψ in Wb, cocos<10 ψ in Wb/rad + exp2pi = (cocos - cmod10)/10 + + ! cocos mod 10 = 1,3,5,7: toroidal angle φ increasing CCW + phiccw = (mod(cmod10, 2)== 1) + + ! cocos mod 10 = 1,2,5,6: ψ increasing with positive Ip + psiincr = (cmod4==1 .or. cmod4==2) + + ! cocos mod 10 = 1,2,7,8: q positive for positive Bφ*Ip + qpos = (cmod10<3 .or. cmod10>6) + end subroutine decode_cocos + +end module diff --git a/src/gray_jetto1beam.f90 b/src/gray_jetto1beam.f90 index 8b94cee..d69dae3 100644 --- a/src/gray_jetto1beam.f90 +++ b/src/gray_jetto1beam.f90 @@ -1,10 +1,11 @@ subroutine gray_jetto1beam(ijetto, mr, mz, r, z, psin, psia, rax, zax, & nbnd, rbnd, zbnd, nrho, psrad, fpol, te, dne, zeff, qpsi, ibeam, & p0mw, alphain, betain, dpdv, jcd, pabs, icd, err) - use const_and_precisions, only: wp_ - use gray_params, only: gray_parameters, gray_data, gray_results - use gray_core, only: gray_main - use gray_plasma, only: numeric_plasma + use const_and_precisions, only : wp_ + use gray_params, only : gray_parameters, gray_results + use gray_core, only : gray_main + use gray_equil, only : numeric_equil, contour + use gray_plasma, only : numeric_plasma implicit none @@ -21,11 +22,12 @@ subroutine gray_jetto1beam(ijetto, mr, mz, r, z, psin, psia, rax, zax, & integer, intent(out) :: err ! local variables - type(gray_parameters) :: params - type(gray_data) :: data - type(numeric_plasma) :: plasma - type(gray_results) :: res - logical, save :: firstcall = .true. + type(gray_parameters), save :: params + type(numeric_equil) :: equil + type(numeric_plasma) :: plasma + type(contour), save :: limiter + type(gray_results) :: res + logical, save :: firstcall = .true. ! Initialisation tasks for the first call only first_call: if (firstcall) then @@ -52,7 +54,6 @@ subroutine gray_jetto1beam(ijetto, mr, mz, r, z, psin, psia, rax, zax, & ! Set a simple limiter following the boundary of the data grid simple_limiter: block use const_and_precisions, only : cm - use types, only : contour ! Avoid clipping out the launcher real(wp_) :: R_launch, R_max @@ -60,37 +61,37 @@ subroutine gray_jetto1beam(ijetto, mr, mz, r, z, psin, psia, rax, zax, & R_max = max(R_launch, R(mr)) ! Convert to a closed contour - data%equilibrium%limiter = contour( & - Rmin=params%misc%rwall, Rmax=R_max, zmin=z(1), zmax=z(mz)) + limiter = contour(Rmin=params%misc%rwall, Rmax=R_max, & + zmin=z(1), zmax=z(mz)) end block simple_limiter end if first_call ! Set MHD equilibrium data init_equilibrium: block - use equilibrium, only : set_equil_spline - use types, only : contour + use gray_equil, only : eqdsk_data - ! Copy argument arrays - data%equilibrium%rv = r - data%equilibrium%zv = z - data%equilibrium%rax = rax - data%equilibrium%rvac = rax - data%equilibrium%zax = zax - data%equilibrium%psinr = psrad - data%equilibrium%fpol = fpol - data%equilibrium%psia = psia - data%equilibrium%psin = psin - data%equilibrium%qpsi = qpsi - data%equilibrium%boundary = contour(rbnd, zbnd) + type(eqdsk_data) :: data + + ! Build EQDSK structure from in-memory data + data%grid_r = r + data%grid_z = z + data%axis = [rax, zax] + data%r_ref = rax + data%psi = psrad + data%fpol = fpol + data%psi_a = psia + data%psi_map = psin + data%q = qpsi + data%boundary = contour(rbnd, zbnd) ! Compute splines - call set_equil_spline(params%equilibrium, data%equilibrium, err) + call equil%init(params%equilibrium, data, err) if (err /= 0) return end block init_equilibrium ! Compute splines of kinetic profiles - call plasma%init(params, psrad, te, dne, zeff, err) + call plasma%init(params, equil, psrad, te, dne, zeff, err) if (err /= 0) return ! Set wave launcher parameters @@ -108,7 +109,7 @@ subroutine gray_jetto1beam(ijetto, mr, mz, r, z, psin, psia, rax, zax, & end block init_antenna ! Call main subroutine for the ibeam-th beam - call gray_main(params, data, plasma, res, err, rhout=sqrt(psrad)) + call gray_main(params, equil, plasma, limiter, res, err, rhout=sqrt(psrad)) write_debug_files: block integer :: i, err @@ -120,14 +121,6 @@ subroutine gray_jetto1beam(ijetto, mr, mz, r, z, psin, psia, rax, zax, & end do end block write_debug_files - ! Free memory - free_memory: block - use equilibrium, only : unset_equil_spline - - ! Unset global variables of the `equilibrium` module - call unset_equil_spline - end block free_memory - ! Copy over the results pabs = res%pabs icd = res%icd diff --git a/src/gray_params.f90 b/src/gray_params.f90 index 97210d7..e08543a 100644 --- a/src/gray_params.f90 +++ b/src/gray_params.f90 @@ -167,30 +167,6 @@ module gray_params integer, allocatable :: active_tables(:) ! IDs of output tables to fill in end type - ! MHD equilibrium data - type equilibrium_data - real(wp_), allocatable :: rv(:) ! R of the uniform grid - real(wp_), allocatable :: zv(:) ! Z of the uniform grid - type(contour) :: limiter ! limiter contour (wall) - type(contour) :: boundary ! boundary contour (plasma) - real(wp_), allocatable :: fpol(:) ! Poloidal current function - real(wp_), allocatable :: qpsi(:) ! Safety factor on the flux grid - real(wp_), allocatable :: psin(:,:) ! Poloidal flux on a uniform grid - real(wp_), allocatable :: psinr(:) ! Poloidal flux - real(wp_) :: psia ! Poloidal flux at edge - flux at magnetic axis - real(wp_) :: rvac ! Reference R₀ (B = B₀R₀/R without the plasma) - real(wp_) :: rax ! R of the magnetic axis - real(wp_) :: zax ! Z of the magnetic axis - end type - - ! Kinetic plasma profiles data - type profiles_data - real(wp_), allocatable :: psrad(:) ! Radial coordinate - real(wp_), allocatable :: terad(:) ! Electron temperature profile - real(wp_), allocatable :: derad(:) ! Electron density profile - real(wp_), allocatable :: zfc(:) ! Effective charge profile - end type - ! All GRAY parameters type gray_parameters type(antenna_parameters) :: antenna @@ -221,12 +197,6 @@ module gray_params "misc.rwall" & ] - ! All GRAY input data - type gray_data - type(equilibrium_data) :: equilibrium - type(profiles_data) :: profiles - end type - ! Wrapper type for array of pointers type table_ptr type(table), pointer :: ptr => null() diff --git a/src/gray_plasma.f90 b/src/gray_plasma.f90 index c44cb19..e5e799f 100644 --- a/src/gray_plasma.f90 +++ b/src/gray_plasma.f90 @@ -262,17 +262,18 @@ contains end function numeric_zeff - subroutine load_plasma(params, plasma, err) + subroutine load_plasma(params, equil, plasma, err) + ! Loads a generic plasma description from file (params%filenm) use gray_params, only : gray_parameters use gray_params, only : PROF_ANALYTIC, PROF_NUMERIC use gray_params, only : RHO_TOR, RHO_POL, RHO_PSI use gray_params, only : SCALE_COLLISION, SCALE_GREENWALD, SCALE_OFF - + use gray_equil, only : abstract_equil use logger, only : log_error, log_debug - use equilibrium, only : frhopol ! subroutine arguments type(gray_parameters), intent(inout) :: params + class(abstract_equil), intent(in) :: equil class(abstract_plasma), allocatable, intent(out) :: plasma integer, intent(out) :: err @@ -362,14 +363,14 @@ contains ! Convert psi to ψ_n (normalised poloidal flux) select case (params%profiles%irho) case (RHO_TOR) ! psi is ρ_t = √Φ_n (toroidal flux) - psi = [(frhopol(psi(i))**2, i=1,nrows)] + psi = [(equil%tor2pol(psi(i))**2, i=1,nrows)] case (RHO_POL) ! psi is ρ_p = √ψ_n (poloidal flux) psi = psi**2 case (RHO_PSI) ! psi is already ψ_n end select ! Interpolate - call np%init(params, psi, temp, dens, zeff, err) + call np%init(params, equil, psi, temp, dens, zeff, err) end block allocate(plasma, source=np) @@ -381,17 +382,19 @@ contains end subroutine load_plasma - subroutine init_splines(self, params, psi, temp, dens, zeff, err) + subroutine init_splines(self, params, equil, psi, temp, dens, zeff, err) ! Computes splines for the plasma profiles data and stores them ! in their respective global variables, see the top of this file. ! ! `err` is 1 if I/O errors occured, 2 if other initialisation failed. use gray_params, only : gray_parameters + use gray_equil, only : abstract_equil use logger, only : log_debug, log_info, log_warning, log_error ! subroutine arguments class(numeric_plasma), intent(out) :: self type(gray_parameters), intent(inout) :: params + class(abstract_equil), intent(in) :: equil real(wp_), dimension(:), intent(in) :: psi, temp, dens, zeff integer, intent(out) :: err @@ -481,7 +484,6 @@ contains ! Note: if it does, the initial wave conditions become ! invalid as they are given assuming a vacuum (N=1) block - use equilibrium, only : pol_flux use const_and_precisions, only : cm real(wp_) :: R, Z, psi @@ -493,7 +495,7 @@ contains ! Get the poloidal flux at the launcher ! Note: this returns -1 when the data is not available - call pol_flux(R, z, psi) + call equil%pol_flux(R, z, psi) if (psi > self%tail%start .and. psi < self%tail%end) then ! Fall back to the midpoint of ψ₀ and the launcher ψ diff --git a/src/gray_tables.f90 b/src/gray_tables.f90 index 22d6f56..600d876 100644 --- a/src/gray_tables.f90 +++ b/src/gray_tables.f90 @@ -110,16 +110,17 @@ contains end subroutine init_tables - function kinetic_profiles(params, plasma) result(tbl) + function kinetic_profiles(params, equil, plasma) result(tbl) ! Generates the plasma kinetic profiles use gray_params, only : gray_parameters, EQ_VACUUM + use gray_equil, only : abstract_equil use gray_plasma, only : abstract_plasma - use equilibrium, only : fq, frhotor use magsurf_data, only : npsi, vajphiav ! function arguments type(gray_parameters), intent(in) :: params + class(abstract_equil), intent(in) :: equil class(abstract_plasma), intent(in) :: plasma type(table) :: tbl @@ -145,7 +146,7 @@ contains do i = 0, npsi + ntail rho_p = i * drho_p - rho_t = frhotor(rho_p) + rho_t = equil%pol2tor(rho_p) psi_n = rho_p**2 if (psi_n < 1) then J_phi = vajphiav(i+1) * 1.e-6_wp_ @@ -153,24 +154,26 @@ contains J_phi = 0 end if call plasma%density(psi_n, dens, ddens) - call tbl%append([psi_n, rho_t, dens, plasma%temp(psi_n), fq(psi_n), J_phi]) + call tbl%append([psi_n, rho_t, dens, plasma%temp(psi_n), & + equil%safety(psi_n), J_phi]) end do end function kinetic_profiles - function ec_resonance(params, B_res) result(tbl) + function ec_resonance(params, equil, B_res) result(tbl) ! Generates the EC resonance surface table use const_and_precisions, only : comp_eps use gray_params, only : gray_parameters, EQ_VACUUM - use equilibrium, only : rmnm, rmxm, zmnm, zmxm, bfield + use gray_equil, only : abstract_equil use magsurf_data, only : npsi use types, only : item use cniteq, only : cniteq_f ! function arguments type(gray_parameters), intent(in) :: params + class(abstract_equil), intent(in) :: equil real(wp_), intent(in) :: B_res type(table) :: tbl @@ -190,11 +193,11 @@ contains real(wp_), dimension(npsi) :: R, z ! Build a regular (R, z) grid - dr = (rmxm - rmnm - comp_eps)/(npsi - 1) - dz = (zmxm - zmnm)/(npsi - 1) + dr = (equil%r_range(2) - equil%r_range(1) - comp_eps)/(npsi - 1) + dz = (equil%z_range(2) - equil%z_range(1))/(npsi - 1) do j=1,npsi - R(j) = comp_eps + rmnm + dr*(j - 1) - z(j) = zmnm + dz*(j - 1) + R(j) = comp_eps + equil%r_range(1) + dr*(j - 1) + z(j) = equil%z_range(1) + dz*(j - 1) end do ! B_tot on psi grid @@ -202,7 +205,7 @@ contains B_min = +1.0e30_wp_ do k = 1, npsi do j = 1, npsi - call bfield(R(j), z(k), B_R, B_z, B_phi) + call equil%b_field(R(j), z(k), B_R, B_z, B_phi) B_tot(j,k) = sqrt(B_R**2 + B_z**2 + B_phi**2) if(B_tot(j,k) >= B_max) B_max = B_tot(j,k) if(B_tot(j,k) <= B_min) B_min = B_tot(j,k) @@ -227,17 +230,18 @@ contains end function ec_resonance - function inputs_maps(params, plasma, B_res, X_norm) result(tbl) + function inputs_maps(params, equil, plasma, B_res, X_norm) result(tbl) ! Generates 2D maps of several input quantities use const_and_precisions, only : comp_eps, cm, degree use gray_params, only : gray_parameters, EQ_VACUUM + use gray_equil, only : abstract_equil use gray_plasma, only : abstract_plasma - use equilibrium, only : rmnm, rmxm, zmnm, zmxm, pol_flux, bfield use magsurf_data, only : npsi ! function arguments type(gray_parameters), intent(in) :: params + class(abstract_equil), intent(in) :: equil class(abstract_plasma), intent(in) :: plasma real(wp_), intent(in) :: B_res ! resonant magnetic field, e/m_eω real(wp_), intent(in) :: X_norm ! X normalisation, e²/ε₀m_eω² @@ -261,18 +265,18 @@ contains Npl0 = sin(params%antenna%beta*degree) ! initial value of N∥ ! Build a regular (R, z) grid - dR = (rmxm - rmnm - comp_eps)/(npsi - 1) - dz = (zmxm - zmnm)/(npsi - 1) + dR = (equil%r_range(2) - equil%r_range(1) - comp_eps)/(npsi - 1) + dz = (equil%z_range(2) - equil%z_range(1))/(npsi - 1) do j = 1, npsi - R(j) = comp_eps + rmnm + dR*(j - 1) - z(j) = zmnm + dz*(j - 1) + R(j) = comp_eps + equil%r_range(1) + dR*(j - 1) + z(j) = equil%z_range(1) + dz*(j - 1) end do do j = 1, npsi Npl = Npl0 * R0/r(j) do k = 1, npsi - call pol_flux(r(j), z(k), psi_n) - call bfield(r(j), z(k), B_R, B_z, B_phi) + call equil%pol_flux(r(j), z(k), psi_n) + call equil%b_field(r(j), z(k), B_R, B_z, B_phi) call plasma%density(psi_n, ne, dne) B = sqrt(B_R**2 + B_phi**2 + B_z**2) X = X_norm*ne @@ -286,13 +290,13 @@ contains end function inputs_maps - subroutine find_flux_surfaces(qvals, params, tbl) + subroutine find_flux_surfaces(qvals, params, equil, tbl) ! Finds the ψ for a set of values of q and stores the ! associated surfaces to the flux surface table use gray_params, only : gray_parameters - use equilibrium, only : fq, frhotor, rmaxis, zmaxis, zbsup, zbinf - use magsurf_data, only : contours_psi, npoints + use gray_equil, only : abstract_equil + use magsurf_data, only : npoints use logger, only : log_info use minpack, only : hybrj1 use types, only : table @@ -300,6 +304,7 @@ contains ! subroutine arguments real(wp_), intent(in) :: qvals(:) type(gray_parameters), intent(in) :: params + class(abstract_equil), intent(in) :: equil type(table), intent(inout) :: tbl ! local variables @@ -319,7 +324,7 @@ contains ! searching near the boundary in case q is not monotonic. sol = [0.8_wp_] ! first guess - ! Solve fq(ψ_n) = qvals(i) for ψ_n + ! Solve q(ψ_n) = qvals(i) for ψ_n call hybrj1(equation, n=1, x=sol, fvec=fvec, fjac=fjac, & ldfjac=1, tol=1e-3_wp_, info=info, wa=wa, lwa=7) ! Solution @@ -334,19 +339,20 @@ contains real(wp_), dimension(npoints) :: R_cont, z_cont real(wp_) :: R_hi, R_lo, z_hi, z_lo - ! Guesses for the high,low horzizontal-tangent points - R_hi = rmaxis; - z_hi = (zbsup + zmaxis)/2 - R_lo = rmaxis - z_lo = (zbinf + zmaxis)/2 + ! Guesses for the high,low horizontal-tangent points + R_hi = equil%axis(1) + z_hi = (equil%z_boundary(2) + equil%axis(2))/2 + R_lo = equil%axis(1) + z_lo = (equil%z_boundary(1) + equil%axis(2))/2 - call contours_psi(params, psi_n, R_cont, z_cont, R_hi, z_hi, R_lo, z_lo) + call equil%flux_contour(psi_n, params%misc%rwall, & + R_cont, z_cont, R_hi, z_hi, R_lo, z_lo) call store_flux_surface(tbl, psi_n, R_cont, z_cont) end block ! Log the values found for ψ_n, ρ_p, ρ_t rho_p = sqrt(psi_n) - rho_t = frhotor(rho_p) + rho_t = equil%pol2tor(rho_p) write (msg, '(4(x,a,"=",g0.3))') & 'q', qvals(i), 'ψ_n', psi_n, 'ρ_p', rho_p, 'ρ_t', rho_t call log_info(msg, mod='gray_tables', proc='find_flux_surfaces') @@ -368,10 +374,10 @@ contains if (flag == 1) then ! return f(x) - f(1) = fq(x(1)) - qvals(i) + f(1) = equil%safety(x(1)) - qvals(i) else ! return f'(x), computed numerically - df(1,1) = (fq(x(1) + e) - fq(x(1) - e)) / (2*e) + df(1,1) = (equil%safety(x(1) + e) - equil%safety(x(1) - e)) / (2*e) end if end subroutine @@ -402,7 +408,7 @@ contains end subroutine store_flux_surface - subroutine store_ray_data(params, tables, & + subroutine store_ray_data(params, equil, tables, & i, jk, s, P0, pos, psi_n, B, b_n, k0, & N_pl, N_pr, N, N_pr_im, n_e, T_e, & alpha, tau, dIds, nhm, nhf, iokhawa, & @@ -410,15 +416,16 @@ contains ! Stores some ray variables and local quantities use const_and_precisions, only : degree, cm - use equilibrium, only : frhotor + use gray_equil, only : abstract_equil use gray_params, only : gray_parameters, gray_tables use beamdata, only : unfold_index ! subroutine arguments ! tables where to store the data - type(gray_parameters), intent(in) :: params - type(gray_tables), intent(inout) :: tables + type(gray_parameters), intent(in) :: params + class(abstract_equil), intent(in) :: equil + type(gray_tables), intent(inout) :: tables integer, intent(in) :: i, jk ! step, ray index real(wp_), intent(in) :: s ! arclength @@ -453,7 +460,7 @@ contains if (jk == 1 .and. tables%central_ray%active) then phi_deg = atan2(pos_m(2), pos_m(1)) / degree if(psi_n >= 0 .and. psi_n <= 1) then - rho_t = frhotor(sqrt(psi_n)) + rho_t = equil%pol2tor(sqrt(psi_n)) else rho_t = 1.0_wp_ end if diff --git a/src/magsurf_data.f90 b/src/magsurf_data.f90 index 58e8a87..e38d1d2 100644 --- a/src/magsurf_data.f90 +++ b/src/magsurf_data.f90 @@ -95,16 +95,16 @@ contains end subroutine dealloc_surfvec - subroutine compute_flux_averages(params, tables) + subroutine compute_flux_averages(params, equil, tables) use const_and_precisions, only : wp_,zero,one,pi,ccj=>mu0inv - use equilibrium, only : btrcen,btaxis,rmaxis,zmaxis,phitedge,zbsup,zbinf, & - bfield,frhotor,fq,tor_curr,psia,pol_flux - use dierckx, only : regrid,coeff_parder + use dierckx, only : regrid, coeff_parder use types, only : table, wrap use gray_params, only : gray_parameters, gray_tables + use gray_equil, only : abstract_equil ! subroutine arguments type(gray_parameters), intent(in) :: params + class(abstract_equil), intent(in) :: equil type(gray_tables), intent(inout), optional :: tables ! local constants @@ -130,8 +130,7 @@ contains real(wp_), dimension(nnintp*nlam) :: ffhlam,dffhlam real(wp_), dimension(lwrk) :: wrk real(wp_), dimension(:), allocatable :: rctemp,zctemp -! common/external functions/variables - real(wp_) :: fpolv,ddpsidrr,ddpsidzz + real(wp_) :: fpolv, ddpsidrr, ddpsidzz npsi=nnintp npoints = 2*ncnt+1 @@ -162,30 +161,31 @@ contains dffhlam(nlam)=-99999.0_wp_ jp=1 - anorm=2.0_wp_*pi*rmaxis/abs(btaxis) + anorm=2.0_wp_*pi*equil%axis(1)/abs(equil%b_axis) dvdpsi=2.0_wp_*pi*anorm - dadpsi=2.0_wp_*pi/abs(btaxis) - b2av=btaxis**2 - ratio_cdator=abs(btaxis/btrcen) + dadpsi=2.0_wp_*pi/abs(equil%b_axis) + b2av=equil%b_axis**2 + ratio_cdator=abs(equil%b_axis/equil%b_centre) ratio_cdbtor=1.0_wp_ ratio_pltor=1.0_wp_ fc=1.0_wp_ - call pol_flux(rmaxis, zmaxis, ddpsidrr=ddpsidrr, ddpsidzz=ddpsidzz) - qq=abs(btaxis)/sqrt(ddpsidrr*ddpsidzz)/psia - ajphiav=-ccj*(ddpsidrr+ddpsidzz)*psia/rmaxis + call equil%pol_flux(equil%axis(1), equil%axis(2), & + ddpsidrr=ddpsidrr, ddpsidzz=ddpsidzz) + qq=abs(equil%b_axis)/sqrt(ddpsidrr*ddpsidzz)/equil%psi_a + ajphiav=-ccj*(ddpsidrr+ddpsidzz)*equil%psi_a/equil%axis(1) psicon(1)=0.0_wp_ - rcon(:,1)=rmaxis - zcon(:,1)=zmaxis + rcon(:,1)=equil%axis(1) + zcon(:,1)=equil%axis(2) pstab(1)=0.0_wp_ rpstab(1)=0.0_wp_ vcurrp(1)=0.0_wp_ vajphiav(1)=ajphiav - bmxpsi(1)=abs(btaxis) - bmnpsi(1)=abs(btaxis) - bav(1)=abs(btaxis) + bmxpsi(1)=abs(equil%b_axis) + bmnpsi(1)=abs(equil%b_axis) + bav(1)=abs(equil%b_axis) rbav(1)=1.0_wp_ - rri(1)=rmaxis + rri(1)=equil%axis(1) varea(1)=0.0_wp_ vvol(1)=0.0_wp_ vratjpl(1)=ratio_pltor @@ -196,21 +196,22 @@ contains dadrhotv(1)=0.0_wp_ dvdrhotv(1)=0.0_wp_ - rup=rmaxis - rlw=rmaxis - zup=zmaxis+(zbsup-zmaxis)/10.0_wp_ - zlw=zmaxis-(zmaxis-zbinf)/10.0_wp_ + rup=equil%axis(1) + rlw=equil%axis(1) + zup=equil%axis(2)+(equil%z_boundary(2)-equil%axis(2))/10.0_wp_ + zlw=equil%axis(2)-(equil%axis(2)-equil%z_boundary(1))/10.0_wp_ do jp=2,npsi height=dble(jp-1)/dble(npsi-1) if(jp.eq.npsi) height=0.9999_wp_ rhopjp=height psinjp=height*height - rhotjp=frhotor(rhopjp) + rhotjp=equil%pol2tor(rhopjp) if (rhotjp /= rhotjp) print *, 'Nan!!!!!!!!!!' psicon(jp)=height - call contours_psi(params, psinjp, rctemp, zctemp, rup, zup, rlw, zlw) + call equil%flux_contour(psinjp, params%misc%rwall, & + rctemp, zctemp, rup, zup, rlw, zlw) rcon(:,jp) = rctemp zcon(:,jp) = zctemp @@ -226,8 +227,8 @@ contains bmmx=-1.0e+30_wp_ bmmn=1.0e+30_wp_ - ajphi0 = tor_curr(rctemp(1),zctemp(1)) - call bfield(rctemp(1), zctemp(1), brr, bzz, bphi) + ajphi0 = equil%tor_curr(rctemp(1), zctemp(1)) + call equil%b_field(rctemp(1), zctemp(1), brr, bzz, bphi) fpolv=bphi*rctemp(1) btot0=sqrt(bphi**2+brr**2+bzz**2) bpoloid0=sqrt(brr**2+bzz**2) @@ -237,8 +238,8 @@ contains do inc=1,npoints-1 inc1=inc+1 - dla=sqrt((rctemp(inc)-rmaxis)**2+(zctemp(inc)-zmaxis)**2) - dlb=sqrt((rctemp(inc1)-rmaxis)**2+(zctemp(inc1)-zmaxis)**2) + dla=sqrt((rctemp(inc)-equil%axis(1))**2+(zctemp(inc)-equil%axis(2))**2) + dlb=sqrt((rctemp(inc1)-equil%axis(1))**2+(zctemp(inc1)-equil%axis(2))**2) dlp=sqrt((rctemp(inc1)-rctemp(inc))**2+(zctemp(inc1)-zctemp(inc))**2) drc=(rctemp(inc1)-rctemp(inc)) @@ -253,8 +254,8 @@ contains ! compute line integrals on the contour psi=psinjp=height^2 rpsim=rctemp(inc1) zpsim=zctemp(inc1) - call bfield(rpsim, zpsim, brr, bzz) - ajphi = tor_curr(rpsim,zpsim) + call equil%b_field(rpsim, zpsim, brr, bzz) + ajphi = equil%tor_curr(rpsim, zpsim) bphi=fpolv/rpsim btot=sqrt(bphi**2+brr**2+bzz**2) bpoloid=sqrt(brr**2+bzz**2) @@ -310,7 +311,7 @@ contains bmxpsi(jp)=bmmx bmnpsi(jp)=bmmn rri(jp)=bav(jp)/abs(fpolv*r2iav) - ratio_cdator=abs(b2av*riav/(fpolv*r2iav*btrcen)) + ratio_cdator=abs(b2av*riav/(fpolv*r2iav*equil%b_centre)) ratio_cdbtor=abs(b2av*riav/(fpolv*r2iav*bbav)) ratio_pltor=abs(bbav*riav/(fpolv*r2iav)) vratjpl(jp)=ratio_pltor @@ -318,8 +319,8 @@ contains vratjb(jp)=ratio_cdbtor qq=abs(dvdpsi*fpolv*r2iav/(4.0_wp_*pi*pi)) qqv(jp)=qq - dadrhotv(jp)=phitedge*rhotjp/fq(psinjp)*dadpsi/pi - dvdrhotv(jp)=phitedge*rhotjp/fq(psinjp)*dvdpsi/pi + dadrhotv(jp)=equil%phi_a*rhotjp/equil%safety(psinjp)*dadpsi/pi + dvdrhotv(jp)=equil%phi_a*rhotjp/equil%safety(psinjp)*dvdpsi/pi ! computation of fraction of circulating/trapped fraction fc, ft ! and of function H(lambda,rhop) @@ -397,8 +398,8 @@ contains do jp=1,npsi if (.not. tables%flux_averages%active) exit call tables%flux_averages%append([ & - rpstab(jp), frhotor(rpstab(jp)), bav(jp), bmxpsi(jp), & - bmnpsi(jp), varea(jp), vvol(jp), vcurrp(jp), vajphiav(jp), & + rpstab(jp), equil%pol2tor(rpstab(jp)), bav(jp), bmxpsi(jp), & + bmnpsi(jp), varea(jp), vvol(jp), vcurrp(jp), vajphiav(jp), & ffc(jp), vratja(jp), vratjb(jp), qqv(jp)]) end do @@ -444,88 +445,4 @@ contains end subroutine fluxval - - - subroutine contours_psi(params, h, rcn, zcn, rup, zup, rlw, zlw) - use const_and_precisions, only : wp_, pi - use gray_params, only : gray_parameters - use logger, only : log_warning - use dierckx, only : profil, sproota - use equilibrium, only : model, frhotor, psi_spline, & - kspl, find_htg_point - - ! local constants - integer, parameter :: mest=4 - - ! subroutine arguments - type(gray_parameters), intent(in) :: params - real(wp_), intent(in) :: h - real(wp_), intent(out) :: rcn(:), zcn(:) - real(wp_), intent(inout) :: rup, zup, rlw, zlw - - ! local variables - integer :: npoints,np,ic,ier,iopt,m - real(wp_) :: ra,rb,za,zb,th,zc - real(wp_), dimension(mest) :: zeroc - real(wp_), dimension(psi_spline%nknots_x) :: czc - - npoints=size(rcn) - np=(npoints-1)/2 - th=pi/dble(np) - - if (params%equilibrium%iequil<2) then - block - real(wp_) :: r_p ! poloidal radius - r_p = sqrt(h) * model%a - do ic=1,npoints - rcn(ic) = model%R0 + r_p * cos(th*(ic-1)) - zcn(ic) = model%z0 + r_p * sin(th*(ic-1)) - end do - end block - else - - ra=rup - rb=rlw - za=zup - zb=zlw - call find_htg_point(ra,za,rup,zup,h) - call find_htg_point(rb,zb,rlw,zlw,h) - - rcn(1)=rlw - zcn(1)=zlw - rcn(npoints)=rlw - zcn(npoints)=zlw - rcn(np+1)=rup - zcn(np+1)=zup - do ic=2,np - zc=zlw+(zup-zlw)*(1.0_wp_-cos(th*(ic-1)))/2.0_wp_ - iopt=1 - call profil(iopt, psi_spline%knots_x, psi_spline%nknots_x, & - psi_spline%knots_y, psi_spline%nknots_y, & - psi_spline%coeffs, kspl, kspl, zc, & - psi_spline%nknots_x, czc, ier) - if (ier > 0) then - block - character(256) :: msg - write(msg, '(a, a, g0)') & - 'when computing ψ(R,z) contour `profil` returned ier=', ier - call log_warning(msg, mod='magsurf_data', proc='contours_psi') - end block - end if - call sproota(h, psi_spline%knots_x, psi_spline%nknots_x, & - czc, zeroc, mest, m, ier) - if (zeroc(1) > params%misc%rwall) then - rcn(ic)=zeroc(1) - rcn(npoints+1-ic)=zeroc(2) - else - rcn(ic)=zeroc(2) - rcn(npoints+1-ic)=zeroc(3) - end if - zcn(ic)=zc - zcn(npoints+1-ic)=zc - end do - - end if - end subroutine contours_psi - end module magsurf_data diff --git a/src/main.f90 b/src/main.f90 index 897fcd2..14e7554 100644 --- a/src/main.f90 +++ b/src/main.f90 @@ -2,12 +2,14 @@ program main use const_and_precisions, only : wp_ use logger, only : INFO, ERROR, WARNING, set_log_level, log_message use utils, only : dirname + use types, only : contour use gray_cli, only : cli_options, parse_cli_options, & parse_param_overrides use gray_core, only : gray_main + use gray_equil, only : abstract_equil, load_equil use gray_plasma, only : abstract_plasma, load_plasma - use gray_params, only : gray_parameters, gray_data, gray_results, & - read_gray_params, read_gray_config, & + use gray_params, only : gray_parameters, gray_results, & + read_gray_params, read_gray_config, & make_gray_header implicit none @@ -18,8 +20,9 @@ program main ! gray_main subroutine arguments type(gray_parameters) :: params ! Inputs - type(gray_data) :: data ! + class(abstract_equil), allocatable :: equil ! class(abstract_plasma), allocatable :: plasma ! + type(contour) :: limiter ! type(gray_results) :: results ! Outputs integer :: err ! Exit code @@ -68,7 +71,7 @@ program main associate (p => params%raytracing) if (p%nrayr < 5) then p%igrad = 0 - call log_message(level=WARNING, mod='main', proc='main', & + call log_message(level=WARNING, mod='main', & msg='nrayr < 5 ⇒ optical case only') end if if (p%nrayr == 1) p%nrayth = 1 @@ -78,29 +81,33 @@ program main end if end associate - ! Read the input data and set the global variables - ! of the respective module. Note: order matters. - call init_equilibrium(params, data, err) + ! Read and initialise the equilibrium and limiter objects + call load_equil(params%equilibrium, equil, limiter, err) if (err /= 0) call exit(err) ! Read and initialise the plasma state object - call load_plasma(params, plasma, err) + call load_plasma(params, equil, plasma, err) if (err /= 0) call exit(err) - call init_misc(params, data) + ! Create a simple limiter, if necessary + if (.not. allocated(limiter%R) .or. params%raytracing%ipass < 0) then + call log_message('using fallback limiter', level=INFO, mod='main') + params%raytracing%ipass = abs(params%raytracing%ipass) + limiter = make_fallback_limiter(params, equil) + end if ! Get back to the initial directory err = chdir(cwd) if (allocated(opts%sum_filelist)) then call log_message(level=INFO, mod='main', msg='summing profiles') - call sum_profiles(params, opts%sum_filelist, opts%output_dir, results) + call sum_profiles(params, equil, opts%sum_filelist, opts%output_dir, results) else ! Activate the given output tables params%misc%active_tables = opts%tables ! Finally run the simulation - call gray_main(params, data, plasma, results, err) + call gray_main(params, equil, plasma, limiter, results, err) end if print_res: block @@ -136,68 +143,10 @@ program main end do end block write_res - ! Free memory - free_mem: block - use equilibrium, only : unset_equil_spline - call unset_equil_spline - end block free_mem - end block no_save contains - subroutine init_equilibrium(params, data, err) - ! Reads the MHD equilibrium file (either in the G-EQDSK format - ! or an analytical description) and initialises the respective - ! GRAY parameters and data. - use gray_params, only : EQ_VACUUM, EQ_ANALYTICAL, EQ_EQDSK_FULL, EQ_EQDSK_PARTIAL - use equilibrium, only : read_equil_an, read_eqdsk, change_cocos, & - set_equil_an, set_equil_spline, scale_equil - use logger, only : log_debug - - ! subroutine arguments - type(gray_parameters), intent(inout) :: params - type(gray_data), intent(out) :: data - integer, intent(out) :: err - - select case (params%equilibrium%iequil) - case (EQ_VACUUM) - call log_debug('vacuum, no MHD equilibrium', & - mod='main', proc='init_equilibrium') - - case (EQ_ANALYTICAL) - call log_debug('loading analytical file', & - mod='main', proc='init_equilibrium') - call read_equil_an(params%equilibrium%filenm, & - params%raytracing%ipass, & - data%equilibrium, err) - if (err /= 0) return - - case (EQ_EQDSK_FULL, EQ_EQDSK_PARTIAL) - call log_debug('loading G-EQDK file', & - mod='main', proc='init_equilibrium') - call read_eqdsk(params%equilibrium, data%equilibrium, err) - if (err /= 0) return - call change_cocos(data%equilibrium, params%equilibrium%icocos, 3) - end select - - ! Rescale B, I and/or force their signs - call scale_equil(params%equilibrium, data%equilibrium) - - ! Set global variables - select case (params%equilibrium%iequil) - case (EQ_ANALYTICAL) - call set_equil_an - - case (EQ_EQDSK_FULL, EQ_EQDSK_PARTIAL) - call log_debug('computing splines...', mod='main', proc='init_equilibrium') - call set_equil_spline(params%equilibrium, data%equilibrium, err) - if (err /= 0) return - call log_debug('splines computed', mod='main', proc='init_equilibrium') - end select - end subroutine init_equilibrium - - subroutine init_antenna(params, err) ! Reads the wave launcher file (containing the wave frequency, launcher ! position, direction and beam description) and initialises the respective @@ -227,82 +176,61 @@ contains end subroutine init_antenna - subroutine init_misc(params, data) - ! Performs miscellanous initial tasks, before the gray_main subroutine. - use gray_params, only : EQ_VACUUM, EQ_ANALYTICAL, & - EQ_EQDSK_FULL, EQ_EQDSK_PARTIAL - use types, only : contour - use const_and_precisions, only : cm, comp_huge + function make_fallback_limiter(params, equil) result(limiter) + ! Creates a simple rectangular limiter: + ! + ! (R_wall, z_launch+Δz)╔═════╗(R_max, z_launch+Δz) + ! ║4 3║ + ! ║ ║ + ! ║1 2║ + ! (R_wall, z_launch-Δz)╚═════╝(R_max, z_launch-Δz) + ! + ! Note: R_max = { +∞ for vacuum, + ! { R₀+a for analytical model, + ! { grid_r[-1] for numerical data. + ! + use logger, only : log_info + use const_and_precisions, only : cm - ! subroutine arguments - type(gray_parameters), intent(inout) :: params - type(gray_data), intent(inout) :: data + ! function arguments + type(gray_parameters), intent(in) :: params + class(abstract_equil), intent(in) :: equil + type(contour) :: limiter - ! Build a basic limiter if one is not provided by equilibrium.txt - if (.not. allocated(data%equilibrium%limiter%R) & - .or. params%raytracing%ipass < 0) then + ! local variables + real(wp_) :: R_launch, z_launch, R_max, delta_z - ! Restore sign of ipass - params%raytracing%ipass = abs(params%raytracing%ipass) + ! Launcher coordinates + R_launch = norm2(params%antenna%pos(1:2)) * cm + z_launch = params%antenna%pos(3) * cm - block - real(wp_) :: R_launch, z_launch, R_max, delta_z + ! Max vertical distance a ray can travel + delta_z = abs(params%raytracing%ipass) * & + params%raytracing%dst * params%raytracing%nstep * cm - ! Launcher coordinates - R_launch = norm2(params%antenna%pos(1:2)) * cm - z_launch = params%antenna%pos(3) * cm + ! Avoid clipping out the launcher + R_max = max(R_launch, equil%r_range(2)) - ! Max vertical distance a ray can travel - delta_z = abs(params%raytracing%ipass) * & - params%raytracing%dst * params%raytracing%nstep * cm + ! Build a rectangle + limiter = contour( & + Rmin=params%misc%rwall, Rmax=R_max, & + zmin=z_launch - delta_z, zmax=z_launch + delta_z) - ! Max radius, either due to the plasma extent or equilibrium grid - select case (params%equilibrium%iequil) - case (EQ_VACUUM) - ! Use a very large R, ~ unbounded - R_max = comp_huge - - case (EQ_ANALYTICAL) - ! use R₀+a - block - use equilibrium, only : model - R_max = model%R0 + model%a - end block - - case (EQ_EQDSK_FULL, EQ_EQDSK_PARTIAL) - ! use max R of the grid - R_max = data%equilibrium%rv(size(data%equilibrium%rv)) - end select - - ! Avoid clipping out the launcher - R_max = max(R_launch, R_max) - - ! Convert to a list of R,z: - ! - ! (R_wall, z_launch+Δz)╔═════╗(R_max, z_launch+Δz) - ! ║4 3║ - ! ║ ║ - ! ║1 2║ - ! (R_wall, z_launch-Δz)╚═════╝(R_max, z_launch-Δz) - ! - data%equilibrium%limiter = contour( & - Rmin=params%misc%rwall, Rmax=R_max, & - zmin=z_launch - delta_z, zmax=z_launch + delta_z) - end block - end if - end subroutine init_misc + end function make_fallback_limiter - subroutine sum_profiles(params, filelist, output_dir, results) + subroutine sum_profiles(params, equil, filelist, output_dir, results) ! Combines the EC profiles obtained from multiple gray runs ! (of different beams) and recomputes the summary table use gray_params, only : gray_parameters + use gray_equil, only : abstract_equil use magsurf_data, only : compute_flux_averages, dealloc_surfvec use pec, only : pec_init, postproc_profiles, dealloc_pec, & rhot_tab ! subroutine arguments type(gray_parameters), intent(inout) :: params + class(abstract_equil), intent(in) :: equil character(len=*), intent(in) :: filelist, output_dir type(gray_results), intent(inout) :: results @@ -322,10 +250,10 @@ contains integer, allocatable :: beams(:) ! Initialise the magsurf_data module - call compute_flux_averages(params) + call compute_flux_averages(params, equil) ! Initialise the output profiles - call pec_init(params%output) + call pec_init(params%output, equil) associate(nrho =>params%output%nrho) allocate(jphi(nrho), currins(nrho), pins(nrho), rtin(nrho), rpin(nrho)) @@ -419,7 +347,7 @@ contains results%icd = currins(params%output%nrho) ! Recompute the summary values - call postproc_profiles( & + call postproc_profiles(equil, & results%Pabs, results%Icd, rhot_tab, results%dPdV, jphi, & rhotpav, drhotpav, rhotjava, drhotjava, dpdvp, jphip, & rhotp, drhotp, rhotj, drhotj, dpdvmx, jphimx, ratjamx, ratjbmx) diff --git a/src/multipass.f90 b/src/multipass.f90 index ef1e123..22264da 100644 --- a/src/multipass.f90 +++ b/src/multipass.f90 @@ -2,26 +2,28 @@ module multipass use const_and_precisions, only : wp_ use polarization, only : pol_limit, field_to_ellipse use reflections, only : wall_refl - use equilibrium, only : bfield + use gray_equil, only : abstract_equil implicit none contains - subroutine plasma_in(i, x, N, Bres, sox, cpl, psi, chi, iop, ext, eyt, perfect) + subroutine plasma_in(i, equil, x, N, Bres, sox, cpl, & + psi, chi, iop, ext, eyt, perfect) ! Computes the ray polarisation and power couplings when it enteres the plasma use const_and_precisions, only : cm ! subroutine arguments - integer, intent(in) :: i ! ray index - real(wp_), intent(in) :: x(3), N(3) ! position, refactive index - real(wp_), intent(in) :: Bres ! resonant B field - integer, intent(in) :: sox ! sign of polarisation mode: -1 ⇒ O, +1 ⇒ X - real(wp_), intent(out) :: cpl(2) ! power coupling vector (O, X) - real(wp_), intent(out) :: psi, chi ! polarisation ellipse angles - integer, intent(inout) :: iop(:) ! inside/outside plasma flag - complex(wp_), intent(inout) :: ext(:), eyt(:) ! ray polarisation vector (e_x, e_y) - logical, intent(in) :: perfect ! whether to assume perfect coupling + integer, intent(in) :: i ! ray index + class(abstract_equil), intent(in) :: equil ! MHD equilibrium + real(wp_), intent(in) :: x(3), N(3) ! position, refactive index + real(wp_), intent(in) :: Bres ! resonant B field + integer, intent(in) :: sox ! sign of polarisation mode: -1 ⇒ O, +1 ⇒ X + real(wp_), intent(out) :: cpl(2) ! power coupling vector (O, X) + real(wp_), intent(out) :: psi, chi ! polarisation ellipse angles + integer, intent(inout) :: iop(:) ! inside/outside plasma flag + complex(wp_), intent(inout) :: ext(:), eyt(:) ! ray polarisation vector (e_x, e_y) + logical, intent(in) :: perfect ! whether to assume perfect coupling ! local variables real(wp_) :: R, z, cosphi, sinphi, B_phi, B_R, B_z @@ -37,7 +39,7 @@ contains z = x(3) * cm cosphi = x(1)/R * cm sinphi = x(2)/R * cm - call bfield(R, z, B_R, B_z, B_phi) + call equil%b_field(R, z, B_R, B_z, B_phi) B(1) = B_R*cosphi - B_phi*sinphi B(2) = B_R*sinphi + B_phi*cosphi B(3) = B_z @@ -70,16 +72,17 @@ contains end subroutine plasma_in - subroutine plasma_out(i, xv, anv, bres, sox, iop, ext, eyt) + subroutine plasma_out(i, equil, xv, anv, bres, sox, iop, ext, eyt) ! Ray exits plasma ! subroutine arguments - integer, intent(in) :: i ! ray index - real(wp_), intent(in) :: xv(3), anv(3) - real(wp_), intent(in) :: bres - integer, intent(in) :: sox - integer, intent(inout) :: iop(:) ! in/out plasma flag - complex(wp_), intent(out) :: ext(:), eyt(:) + integer, intent(in) :: i ! ray index + class(abstract_equil), intent(in) :: equil ! MHD equilibrium + real(wp_), intent(in) :: xv(3), anv(3) + real(wp_), intent(in) :: bres + integer, intent(in) :: sox + integer, intent(inout) :: iop(:) ! in/out plasma flag + complex(wp_), intent(out) :: ext(:), eyt(:) ! local variables real(wp_) :: rm, csphi, snphi, bphi, br, bz @@ -91,7 +94,7 @@ contains rm=sqrt(xmv(1)**2+xmv(2)**2) csphi=xmv(1)/rm snphi=xmv(2)/rm - call bfield(rm,xmv(3),br,bz,bphi) + call equil%b_field(rm,xmv(3),br,bz,bphi) bv(1)=br*csphi-bphi*snphi bv(2)=br*snphi+bphi*csphi bv(3)=bz @@ -99,34 +102,35 @@ contains end subroutine plasma_out - subroutine wall_out(i, wall, ins, xv, anv, dst, bres, sox, & - psipol1, chipol1, iow, iop, ext, eyt) + subroutine wall_out(i, equil, wall, ins, xv, anv, dst, bres, & + sox, psipol1, chipol1, iow, iop, ext, eyt) ! Ray exits vessel use types, only : contour ! subroutine arguments - integer, intent(in) :: i ! ray index - type(contour), intent(in) :: wall ! wall contour - logical, intent(in) :: ins ! inside plasma? (ins=1 plasma/wall overlap) - real(wp_), intent(inout) :: xv(3), anv(3) - real(wp_), intent(in) :: dst ! step size - real(wp_), intent(in) :: bres - integer, intent(in) :: sox - real(wp_), intent(out) :: psipol1, chipol1 - integer, intent(inout) :: iow(:) ! in/out vessel and plasma flags - integer, intent(inout) :: iop(:) ! in/out vessel and plasma flags - complex(wp_), intent(inout) :: ext(:), eyt(:) + integer, intent(in) :: i ! ray index + class(abstract_equil), intent(in) :: equil ! MHD equilibrium + type(contour), intent(in) :: wall ! wall contour + logical, intent(in) :: ins ! inside plasma? (ins=1 plasma/wall overlap) + real(wp_), intent(inout) :: xv(3), anv(3) + real(wp_), intent(in) :: dst ! step size + real(wp_), intent(in) :: bres + integer, intent(in) :: sox + real(wp_), intent(out) :: psipol1, chipol1 + integer, intent(inout) :: iow(:) ! in/out vessel and plasma flags + integer, intent(inout) :: iop(:) ! in/out vessel and plasma flags + complex(wp_), intent(inout) :: ext(:), eyt(:) ! local variables integer :: irfl real(wp_), dimension(3) :: xvrfl,anvrfl,walln complex(wp_) :: ext1,eyt1 - iow(i)=iow(i)+1 ! out->in - if(ins) call plasma_out(i,xv,anv,bres,sox,iop,ext,eyt) ! plasma-wall overlapping + iow(i) = iow(i) + 1 ! out->in + if(ins) call plasma_out(i, equil, xv, anv, bres, sox, iop, ext, eyt) ! plasma-wall overlapping call wall_refl(wall, xv-dst*anv, anv, ext(i), eyt(i), & - xvrfl, anvrfl, ext1, eyt1, walln, irfl) ! ray reflects at wall - ext(i) = ext1 ! save parameters at wall reflection + xvrfl, anvrfl, ext1, eyt1, walln, irfl) ! ray reflects at wall + ext(i) = ext1 ! save parameters at wall reflection eyt(i) = eyt1 xv = xvrfl anv = anvrfl diff --git a/src/pec.f90 b/src/pec.f90 index eb41908..8bfc2d4 100644 --- a/src/pec.f90 +++ b/src/pec.f90 @@ -10,13 +10,14 @@ module pec contains - subroutine pec_init(params, rt_in) + subroutine pec_init(params, equil, rt_in) use gray_params, only : output_parameters, RHO_POL - use equilibrium, only : frhotor, frhopol + use gray_equil, only : abstract_equil use magsurf_data, only : fluxval ! subroutine arguments type(output_parameters), intent(in) :: params + class(abstract_equil), intent(in) :: equil real(wp_), intent(in), optional :: rt_in(params%nrho) ! local variables @@ -60,12 +61,12 @@ contains end if if (params%ipec == RHO_POL) then rhop_tab(it) = rt - rhot_tab(it) = frhotor(rt) + rhot_tab(it) = equil%pol2tor(rt) rhop1 = rt1 else rhot_tab(it) = rt - rhop_tab(it) = frhopol(rt) - rhop1 = frhopol(rt1) + rhop_tab(it) = equil%tor2pol(rt) + rhop1 = equil%tor2pol(rt1) end if ! psi grid at mid points, size n+1, for use in pec_tab @@ -242,23 +243,24 @@ contains end subroutine pec_tab - subroutine postproc_profiles(pabs,currt,rhot_tab,dpdv,ajphiv, & + subroutine postproc_profiles(equil, pabs, currt, rhot_tab, dpdv, ajphiv, & rhotpav,drhotpav,rhotjava,drhotjava,dpdvp, ajphip, & rhotp, drhotp, rhotjfi, drhotjfi, dpdvmx,ajmxfi, ratjamx, ratjbmx) ! Radial average values over power and current density profile use const_and_precisions, only : pi, comp_eps - use equilibrium, only : frhopol - use magsurf_data, only : fluxval + use gray_equil, only : abstract_equil + use magsurf_data, only : fluxval - real(wp_), intent(in) :: pabs,currt - real(wp_), intent(in) :: rhot_tab(:) - real(wp_), intent(in) :: dpdv(:), ajphiv(:) - real(wp_), intent(out) :: rhotpav, drhotpav, dpdvp - real(wp_), intent(out) :: rhotjava, drhotjava, ajphip - real(wp_), intent(out) :: rhotp, drhotp, dpdvmx - real(wp_), intent(out) :: rhotjfi, drhotjfi, ajmxfi - real(wp_), intent(out) :: ratjamx, ratjbmx + class(abstract_equil), intent(in) :: equil + real(wp_), intent(in) :: pabs,currt + real(wp_), intent(in) :: rhot_tab(:) + real(wp_), intent(in) :: dpdv(:), ajphiv(:) + real(wp_), intent(out) :: rhotpav, drhotpav, dpdvp + real(wp_), intent(out) :: rhotjava, drhotjava, ajphip + real(wp_), intent(out) :: rhotp, drhotp, dpdvmx + real(wp_), intent(out) :: rhotjfi, drhotjfi, ajmxfi + real(wp_), intent(out) :: ratjamx, ratjbmx real(wp_) :: sccsa, ratjplmx, rhopjava, rhoppav real(wp_) :: rhotjav, rhot2pav, rhot2java, dvdrhotav, dadrhotava @@ -287,8 +289,8 @@ contains drhotpav = sqrt(8._wp_*(max(rhot2pav -rhotpav**2,comp_eps))) drhotjava = sqrt(8._wp_*(max(rhot2java-rhotjava**2,comp_eps))) - rhoppav = frhopol(rhotpav) - rhopjava = frhopol(rhotjava) + rhoppav = equil%tor2pol(rhotpav) + rhopjava = equil%tor2pol(rhotjava) if (pabs > 0) then call fluxval(rhoppav,dvdrhot=dvdrhotav) diff --git a/src/splines.f90 b/src/splines.f90 index e1175d5..a3a83c6 100644 --- a/src/splines.f90 +++ b/src/splines.f90 @@ -58,6 +58,7 @@ module splines procedure :: init_deriv => spline_2d_init_deriv procedure :: deriv => spline_2d_deriv procedure :: transform => spline_2d_transform + final :: spline_2d_final end type ! A simple 1D linear interpolation @@ -507,6 +508,23 @@ contains end subroutine spline_2d_deinit + subroutine spline_2d_final(self) + ! Deallocates pointer components in a spline_2d + + type(spline_2d), intent(inout) :: self + + if (allocated(self%partial)) then + deallocate(self%partial(1, 0)%ptr) + deallocate(self%partial(0, 1)%ptr) + deallocate(self%partial(1, 1)%ptr) + deallocate(self%partial(2, 0)%ptr) + deallocate(self%partial(0, 2)%ptr) + deallocate(self%partial) + end if + + end subroutine spline_2d_final + + function spline_2d_eval(self, x, y) result(z) ! Evaluates the spline at (x, y) use dierckx, only : fpbisp