From 24edfdc43ac5eff7a082d691b306163d2ea479d9 Mon Sep 17 00:00:00 2001 From: Michele Guerini Rocco Date: Thu, 26 Oct 2023 18:15:10 +0200 Subject: [PATCH] src/gray_core.f90: print_* do not depend on q_spline MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit This changes all the print_* subroutines to not depend on q_spline, either for the ψ data points or just their number. In fact, in the new analytical model q_spline doesn't exist anymore. Pro: This allows to mark q_spline as private. --- src/equilibrium.f90 | 10 +- src/gray_core.f90 | 223 ++++++++++++++++++++++++-------------------- 2 files changed, 131 insertions(+), 102 deletions(-) diff --git a/src/equilibrium.f90 b/src/equilibrium.f90 index 0c9000e..9d6de83 100644 --- a/src/equilibrium.f90 +++ b/src/equilibrium.f90 @@ -47,7 +47,7 @@ module equilibrium real(wp_), save :: rcen ! Major radius R₀, computed as fpolas/btrcen real(wp_), save :: rmnm, rmxm ! R interval of the equilibrium definition real(wp_), save :: zmnm, zmxm ! z interval of the equilibrium definition - real(wp_), save :: zbinf, zbsup + real(wp_), save :: zbinf, zbsup ! z interval of the plasma boundary private public read_eqdsk, read_equil_an ! Reading data files @@ -59,7 +59,7 @@ module equilibrium public unset_equil_spline ! Deinitialising internal state ! Members exposed for magsurf_data - public kspl, psi_spline, q_spline, points_tgo, model + public kspl, psi_spline, points_tgo, model ! Members exposed to gray_core and more public psia, phitedge @@ -1099,9 +1099,13 @@ contains real(wp_), intent(in) :: psi_n real(wp_) :: J_phi + ! local constants + real(wp_), parameter :: eps = 1.e-4_wp_ + ! local variables real(wp_) :: R1, R2 - call psi_raxis(psi_n, R1, R2) + + call psi_raxis(max(eps, psi_n), R1, R2) J_phi = tor_curr(R2, zmaxis) end function tor_curr_psi diff --git a/src/gray_core.f90 b/src/gray_core.f90 index 792dcb5..8da0ded 100644 --- a/src/gray_core.f90 +++ b/src/gray_core.f90 @@ -2175,50 +2175,34 @@ bb: do ! Prints the (input) plasma kinetic profiles (unit 55) use gray_params, only : profiles_parameters - use equilibrium, only : q_spline, fq, frhotor, tor_curr_psi + use equilibrium, only : fq, frhotor, tor_curr_psi use coreprofiles, only : density, temp use units, only : uprfin, unit_active + use magsurf_data, only : npsi implicit none ! suborutine arguments type(profiles_parameters), intent(in) :: params - ! local constants - real(wp_), parameter :: eps = 1.e-4_wp_ - ! local variables - integer :: i, N_data, N_tail - real(wp_) :: rhot, jphi, dens, ddens - real(wp_) :: psin, dpsin, psin_last + integer :: i + real(wp_) :: rho_t, J_phi, dens, ddens + real(wp_) :: psi_n, dpsi_n if (.not. unit_active(uprfin)) return - ! N of profiles data points - N_data = q_spline%ndata - - ! parameters for the density tail (numerical profiles only) - if (params%iprof == 1) then - N_tail = N_data / 10 ! N of density tail points - psin_last = q_spline%data(N_data) ! last data point - dpsin = (params%psnbnd - psin_last) / N_tail ! Δψ for uniform grid - else - N_tail = 0 - end if + ! Δψ for the uniform ψ grid + dpsi_n = params%psnbnd / (npsi - 1) write (uprfin, *) '#psi rhot ne Te q Jphi' - do i = 1, N_data + N_tail - if (i > N_data) then - psin = psin_last + dpsin*(i - N_data) - else - psin = q_spline%data(i) - end if - rhot = frhotor(sqrt(psin)) - call density(psin, dens, ddens) - jphi = tor_curr_psi(max(eps, psin)) - + do i = 0, npsi - 1 + psi_n = dpsi_n * i + rho_t = frhotor(sqrt(psi_n)) + J_phi = tor_curr_psi(psi_n) + call density(psi_n, dens, ddens) write (uprfin, '(12(1x,e12.5))') & - psin, rhot, dens, temp(psin), fq(psin), jphi*1.e-6_wp_ + psi_n, rho_t, dens, temp(psi_n), fq(psi_n), J_phi*1.e-6_wp_ end do end subroutine print_prof @@ -2226,8 +2210,9 @@ bb: do subroutine print_bres(bres) ! Prints the EC resonance surface table (unit 70) - use equilibrium, only : rmnm, rmxm, zmnm, zmxm, bfield, q_spline + use equilibrium, only : rmnm, rmxm, zmnm, zmxm, bfield use units, only : ubres, unit_active + use magsurf_data, only : npsi implicit none @@ -2242,14 +2227,15 @@ bb: do integer, dimension(10) :: ncpts real(wp_) :: dr,dz,btmx,btmn,zzk,rrj,bbphi,bbr,bbz,bbb real(wp_), dimension(icmx) :: rrcb,zzcb - real(wp_) :: rv(q_spline%ndata), zv(q_spline%ndata) - real(wp_), dimension(q_spline%ndata,q_spline%ndata) :: btotal + real(wp_) :: rv(npsi), zv(npsi) + real(wp_), dimension(npsi,npsi) :: btotal if (.not. unit_active(ubres)) return - dr = (rmxm - rmnm)/(q_spline%ndata - 1) - dz = (zmxm - zmnm)/(q_spline%ndata - 1) - do j=1,q_spline%ndata + ! Build a regular (R, z) grid + dr = (rmxm - rmnm)/(npsi - 1) + dz = (zmxm - zmnm)/(npsi - 1) + do j=1,npsi rv(j) = rmnm + dr*(j - 1) zv(j) = zmnm + dz*(j - 1) end do @@ -2257,9 +2243,9 @@ bb: do ! Btotal on psi grid btmx = -1.0e30_wp_ btmn = 1.0e30_wp_ - do k = 1, q_spline%ndata + do k = 1, npsi zzk = zv(k) - do j = 1, q_spline%ndata + do j = 1, npsi rrj = rv(j) call bfield(rrj, zzk, bbr, bbz, bbphi) btotal(j,k) = sqrt(bbr**2 + bbz**2 + bbphi**2) @@ -2275,7 +2261,7 @@ bb: do if (bbb >= btmn .and. bbb <= btmx) then nconts = size(ncpts) nctot = size(rrcb) - call cniteq(rv, zv, btotal, q_spline%ndata, q_spline%ndata, bbb, & + call cniteq(rv, zv, btotal, npsi, npsi, bbb, & nconts, ncpts, nctot, rrcb, zzcb) do j = 1, nctot write (ubres, '(i6,12(1x,e12.5))') j, bbb, rrcb(j), zzcb(j) @@ -2287,52 +2273,51 @@ bb: do end subroutine print_bres - subroutine print_maps(bres, xgcn, r0, anpl0) + subroutine print_maps(B_res, xgcn, R0, Npl0) ! Prints several input quantities on a regular grid (unit 72) - use equilibrium, only : rmnm, rmxm, zmnm, zmxm, & - pol_flux, pol_curr, q_spline + use equilibrium, only : rmnm, rmxm, zmnm, zmxm, pol_flux, bfield use coreprofiles, only : density, temp use units, only : umaps, unit_active + use magsurf_data, only : npsi implicit none ! subroutine arguments - real(wp_), intent(in) :: bres,xgcn,r0,anpl0 + real(wp_), intent(in) :: B_res ! resonant magnetic field, e/m_eω + real(wp_), intent(in) :: xgcn ! X normalisation, e²/ε₀m_eω² + real(wp_), intent(in) :: R0 ! initial value of R + real(wp_), intent(in) :: Npl0 ! initial value of N∥ ! local variables - integer :: j,k - real(wp_) :: dr, dz, zk, rj, bphi, br, bz, btot, & - psin, ne, dne, te, xg, yg, anpl - real(wp_), dimension(q_spline%ndata) :: r, z + integer :: j, k + real(wp_) :: dR, dz, B_R, B_phi, B_z, B + real(wp_) :: psi_n, ne, dne, Te, X, Y, Npl + real(wp_), dimension(npsi) :: R, z if (.not. unit_active(umaps)) return - dr = (rmxm-rmnm)/(q_spline%ndata - 1) - dz = (zmxm-zmnm)/(q_spline%ndata - 1) - do j=1,q_spline%ndata - r(j) = rmnm + dr*(j - 1) + ! Build a regular (R, z) grid + dR = (rmxm - rmnm)/(npsi - 1) + dz = (zmxm - zmnm)/(npsi - 1) + do j = 1, npsi + R(j) = rmnm + dR*(j - 1) z(j) = zmnm + dz*(j - 1) end do write (umaps, *) '#R z psin Br Bphi Bz Btot ne Te X Y Npl' - do j = 1, q_spline%ndata - rj = r(j) - anpl = anpl0 * r0/rj - do k = 1, q_spline%ndata - zk = z(k) - call pol_flux(rj, zk, psin, dpsidr=bz, dpsidz=br) - call pol_curr(psin, fpol=bphi) - br = -br/rj - bphi = bphi/rj - bz = bz/rj - btot = sqrt(br**2 + bphi**2 + bz**2) - yg = btot/bres - te = temp(psin) - call density(psin, ne, dne) - xg = xgcn*ne + 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 density(psi_n, ne, dne) + B = sqrt(B_R**2 + B_phi**2 + B_z**2) + X = xgcn*ne + Y = B/B_res + Te = temp(psi_n) write (umaps,'(12(x,e12.5))') & - rj, zk, psin, br, bphi, bz, btot, ne, te, xg, yg, anpl + R(j), z(k), psi_n, B_R, B_phi, B_z, B, ne, Te, X, Y, Npl end do write (umaps,*) end do @@ -2340,54 +2325,94 @@ bb: do end subroutine print_maps - subroutine print_surfq(qval) - ! Print constant ψ surfaces for a given `q` value + subroutine print_surfq(qvals) + ! Prints the constant ψ surfaces for a set of values of q - use equilibrium, only : q_spline, fq, frhotor, & - rmaxis, zmaxis, zbsup, zbinf + use equilibrium, only : fq, frhotor, rmaxis, zmaxis, zbsup, zbinf use magsurf_data, only : contours_psi, npoints, print_contour - use utils, only : locate, intlin use logger, only : log_info + use minpack, only : hybrj1 implicit none ! subroutine arguments - real(wp_), dimension(:), intent(in) :: qval + real(wp_), intent(in) :: qvals(:) ! local variables - integer :: i1,i - real(wp_) :: rup,zup,rlw,zlw,rhot,psival - real(wp_), dimension(npoints) :: rcn,zcn - real(wp_), dimension(q_spline%ndata) :: qpsi + integer :: i + real(wp_) :: rho_t, rho_p, psi_n character(256) :: msg ! for log messages formatting - ! build the q profile on the ψ grid - do i = 1, q_spline%ndata - qpsi(i) = fq(q_spline%data(i)) - end do - - ! locate ψ surface for q=qval call log_info('constant ψ surfaces for:', & mod='gray_core', proc='print_surfq') - do i=1, size(qval) - ! FIXME: check for non monotonous q profile - call locate(abs(qpsi), q_spline%ndata, qval(i), i1) - if (i1 > 0 .and. i1 < q_spline%ndata) then - call intlin(abs(qpsi(i1)), q_spline%data(i1), abs(qpsi(i1+1)), & - q_spline%data(i1+1), qval(i), psival) - rup = rmaxis - rlw = rmaxis - zup = (zbsup + zmaxis)/2.0_wp_ - zlw = (zmaxis + zbinf)/2.0_wp_ - call contours_psi(psival, rcn, zcn, rup, zup, rlw, zlw) - call print_contour(psival, rcn, zcn) - rhot = frhotor(sqrt(psival)) - write (msg, '(4(x,a,"=",g0.3))') & - 'q', qval(i), 'ψ', psival, 'rhop', sqrt(psival), 'rhot', rhot - call log_info(msg, mod='gray_core', proc='print_surfq') - end if + do i = 1, size(qvals) + ! Find value of ψ_n for the current q + block + real(wp_) :: sol(1), fvec(1), fjac(1,1), wa(7) + integer :: info + + ! Note: since we are mostly interested in q=3/2,2 we start + ! searching near the boundary in case q is not monotonic. + sol = [0.8_wp_] ! first guess + + ! Solve fq(ψ_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 + psi_n = sol(1) + + ! Handle spurious or no solutions + if (info /= 1 .or. psi_n < 0 .or. psi_n > 1) cycle + end block + + ! Compute and print the ψ_n(R,z) = ψ_n₀ contour + block + 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 + + call contours_psi(psi_n, R_cont, z_cont, R_hi, z_hi, R_lo, z_lo) + call print_contour(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) + 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_core', proc='print_surfq') end do + contains + + subroutine equation(n, x, f, df, ldf, flag) + ! The equation to solve: f(x) = q(x) - q₀ = 0 + use const_and_precisions, only : comp_eps + + implicit none + + ! subroutine arguments + integer, intent(in) :: n, ldf, flag + real(wp_), intent(in) :: x(n) + real(wp_), intent(inout) :: f(n), df(ldf,n) + + ! optimal step size + real(wp_), parameter :: e = comp_eps**(1/3.0_wp_) + + if (flag == 1) then + ! return f(x) + f(1) = fq(x(1)) - qvals(i) + else + ! return f'(x), computed numerically + df(1,1) = (fq(x(1) + e) - fq(x(1) - e)) / (2*e) + end if + end subroutine + end subroutine print_surfq