istp/gray
2
0
Fork 0
Code Issues Pull Requests Packages Projects Releases 1 Wiki Activity

src/gray_core.f90: print_* do not depend on q_spline

Browse Source 
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.
This commit is contained in:
Michele Guerini Rocco 2023-10-26 18:15:10 +02:00
parent 8d86d70e91
commit 24edfdc43a
Signed by: rnhmjoj
GPG Key ID: BFBAF4C975F76450
2 changed files with 131 additions and 102 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

10
src/equilibrium.f90
View File

@ -47,7 +47,7 @@ module equilibrium
real(wp_), save :: rcen ! Major radius R₀, computed as fpolas/btrcen 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 :: rmnm, rmxm ! R interval of the equilibrium definition
real(wp_), save :: zmnm, zmxm ! z 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 private
public read_eqdsk, read_equil_an ! Reading data files public read_eqdsk, read_equil_an ! Reading data files
@ -59,7 +59,7 @@ module equilibrium
public unset_equil_spline ! Deinitialising internal state public unset_equil_spline ! Deinitialising internal state
! Members exposed for magsurf_data ! 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 ! Members exposed to gray_core and more
public psia, phitedge public psia, phitedge
@ -1099,9 +1099,13 @@ contains
real(wp_), intent(in) :: psi_n real(wp_), intent(in) :: psi_n
real(wp_) :: J_phi real(wp_) :: J_phi
! local constants
real(wp_), parameter :: eps = 1.e-4_wp_
! local variables ! local variables
real(wp_) :: R1, R2 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) J_phi = tor_curr(R2, zmaxis)
end function tor_curr_psi end function tor_curr_psi

217
src/gray_core.f90
View File

@ -2175,50 +2175,34 @@ bb: do
! Prints the (input) plasma kinetic profiles (unit 55) ! Prints the (input) plasma kinetic profiles (unit 55)
use gray_params, only : profiles_parameters 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 coreprofiles, only : density, temp
use units, only : uprfin, unit_active use units, only : uprfin, unit_active
use magsurf_data, only : npsi
implicit none implicit none
! suborutine arguments ! suborutine arguments
type(profiles_parameters), intent(in) :: params type(profiles_parameters), intent(in) :: params
! local constants
real(wp_), parameter :: eps = 1.e-4_wp_
! local variables ! local variables
integer :: i, N_data, N_tail integer :: i
real(wp_) :: rhot, jphi, dens, ddens real(wp_) :: rho_t, J_phi, dens, ddens
real(wp_) :: psin, dpsin, psin_last real(wp_) :: psi_n, dpsi_n
if (.not. unit_active(uprfin)) return if (.not. unit_active(uprfin)) return
! N of profiles data points ! Δψ for the uniform ψ grid
N_data = q_spline%ndata dpsi_n = params%psnbnd / (npsi - 1)
! 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
write (uprfin, *) '#psi rhot ne Te q Jphi' write (uprfin, *) '#psi rhot ne Te q Jphi'
do i = 1, N_data + N_tail do i = 0, npsi - 1
if (i > N_data) then psi_n = dpsi_n * i
psin = psin_last + dpsin*(i - N_data) rho_t = frhotor(sqrt(psi_n))
else J_phi = tor_curr_psi(psi_n)
psin = q_spline%data(i) call density(psi_n, dens, ddens)
end if
rhot = frhotor(sqrt(psin))
call density(psin, dens, ddens)
jphi = tor_curr_psi(max(eps, psin))
write (uprfin, '(12(1x,e12.5))') & 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 do
end subroutine print_prof end subroutine print_prof
@ -2226,8 +2210,9 @@ bb: do
subroutine print_bres(bres) subroutine print_bres(bres)
! Prints the EC resonance surface table (unit 70) ! 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 units, only : ubres, unit_active
use magsurf_data, only : npsi
implicit none implicit none
@ -2242,14 +2227,15 @@ bb: do
integer, dimension(10) :: ncpts integer, dimension(10) :: ncpts
real(wp_) :: dr,dz,btmx,btmn,zzk,rrj,bbphi,bbr,bbz,bbb real(wp_) :: dr,dz,btmx,btmn,zzk,rrj,bbphi,bbr,bbz,bbb
real(wp_), dimension(icmx) :: rrcb,zzcb real(wp_), dimension(icmx) :: rrcb,zzcb
real(wp_) :: rv(q_spline%ndata), zv(q_spline%ndata) real(wp_) :: rv(npsi), zv(npsi)
real(wp_), dimension(q_spline%ndata,q_spline%ndata) :: btotal real(wp_), dimension(npsi,npsi) :: btotal
if (.not. unit_active(ubres)) return if (.not. unit_active(ubres)) return
dr = (rmxm - rmnm)/(q_spline%ndata - 1) ! Build a regular (R, z) grid
dz = (zmxm - zmnm)/(q_spline%ndata - 1) dr = (rmxm - rmnm)/(npsi - 1)
do j=1,q_spline%ndata dz = (zmxm - zmnm)/(npsi - 1)
do j=1,npsi
rv(j) = rmnm + dr*(j - 1) rv(j) = rmnm + dr*(j - 1)
zv(j) = zmnm + dz*(j - 1) zv(j) = zmnm + dz*(j - 1)
end do end do
@ -2257,9 +2243,9 @@ bb: do
! Btotal on psi grid ! Btotal on psi grid
btmx = -1.0e30_wp_ btmx = -1.0e30_wp_
btmn = 1.0e30_wp_ btmn = 1.0e30_wp_
do k = 1, q_spline%ndata do k = 1, npsi
zzk = zv(k) zzk = zv(k)
do j = 1, q_spline%ndata do j = 1, npsi
rrj = rv(j) rrj = rv(j)
call bfield(rrj, zzk, bbr, bbz, bbphi) call bfield(rrj, zzk, bbr, bbz, bbphi)
btotal(j,k) = sqrt(bbr**2 + bbz**2 + bbphi**2) btotal(j,k) = sqrt(bbr**2 + bbz**2 + bbphi**2)
@ -2275,7 +2261,7 @@ bb: do
if (bbb >= btmn .and. bbb <= btmx) then if (bbb >= btmn .and. bbb <= btmx) then
nconts = size(ncpts) nconts = size(ncpts)
nctot = size(rrcb) 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) nconts, ncpts, nctot, rrcb, zzcb)
do j = 1, nctot do j = 1, nctot
write (ubres, '(i6,12(1x,e12.5))') j, bbb, rrcb(j), zzcb(j) write (ubres, '(i6,12(1x,e12.5))') j, bbb, rrcb(j), zzcb(j)
@ -2287,52 +2273,51 @@ bb: do
end subroutine print_bres 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) ! Prints several input quantities on a regular grid (unit 72)
use equilibrium, only : rmnm, rmxm, zmnm, zmxm, & use equilibrium, only : rmnm, rmxm, zmnm, zmxm, pol_flux, bfield
pol_flux, pol_curr, q_spline
use coreprofiles, only : density, temp use coreprofiles, only : density, temp
use units, only : umaps, unit_active use units, only : umaps, unit_active
use magsurf_data, only : npsi
implicit none implicit none
! subroutine arguments ! 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 ! local variables
integer :: j, k integer :: j, k
real(wp_) :: dr, dz, zk, rj, bphi, br, bz, btot, & real(wp_) :: dR, dz, B_R, B_phi, B_z, B
psin, ne, dne, te, xg, yg, anpl real(wp_) :: psi_n, ne, dne, Te, X, Y, Npl
real(wp_), dimension(q_spline%ndata) :: r, z real(wp_), dimension(npsi) :: R, z
if (.not. unit_active(umaps)) return if (.not. unit_active(umaps)) return
dr = (rmxm-rmnm)/(q_spline%ndata - 1) ! Build a regular (R, z) grid
dz = (zmxm-zmnm)/(q_spline%ndata - 1) dR = (rmxm - rmnm)/(npsi - 1)
do j=1,q_spline%ndata dz = (zmxm - zmnm)/(npsi - 1)
r(j) = rmnm + dr*(j - 1) do j = 1, npsi
R(j) = rmnm + dR*(j - 1)
z(j) = zmnm + dz*(j - 1) z(j) = zmnm + dz*(j - 1)
end do end do
write (umaps, *) '#R z psin Br Bphi Bz Btot ne Te X Y Npl' write (umaps, *) '#R z psin Br Bphi Bz Btot ne Te X Y Npl'
do j = 1, q_spline%ndata do j = 1, npsi
rj = r(j) Npl = Npl0 * R0/r(j)
anpl = anpl0 * r0/rj do k = 1, npsi
do k = 1, q_spline%ndata call pol_flux(r(j), z(k), psi_n)
zk = z(k) call bfield(r(j), z(k), B_R, B_z, B_phi)
call pol_flux(rj, zk, psin, dpsidr=bz, dpsidz=br) call density(psi_n, ne, dne)
call pol_curr(psin, fpol=bphi) B = sqrt(B_R**2 + B_phi**2 + B_z**2)
br = -br/rj X = xgcn*ne
bphi = bphi/rj Y = B/B_res
bz = bz/rj Te = temp(psi_n)
btot = sqrt(br**2 + bphi**2 + bz**2)
yg = btot/bres
te = temp(psin)
call density(psin, ne, dne)
xg = xgcn*ne
write (umaps,'(12(x,e12.5))') & 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 end do
write (umaps,*) write (umaps,*)
end do end do
@ -2340,54 +2325,94 @@ bb: do
end subroutine print_maps end subroutine print_maps
subroutine print_surfq(qval) subroutine print_surfq(qvals)
! Print constant ψ surfaces for a given `q` value ! Prints the constant ψ surfaces for a set of values of q
use equilibrium, only : q_spline, fq, frhotor, & use equilibrium, only : fq, frhotor, rmaxis, zmaxis, zbsup, zbinf
rmaxis, zmaxis, zbsup, zbinf
use magsurf_data, only : contours_psi, npoints, print_contour use magsurf_data, only : contours_psi, npoints, print_contour
use utils, only : locate, intlin
use logger, only : log_info use logger, only : log_info
use minpack, only : hybrj1
implicit none implicit none
! subroutine arguments ! subroutine arguments
real(wp_), dimension(:), intent(in) :: qval real(wp_), intent(in) :: qvals(:)
! local variables ! local variables
integer :: i1,i integer :: i
real(wp_) :: rup,zup,rlw,zlw,rhot,psival real(wp_) :: rho_t, rho_p, psi_n
real(wp_), dimension(npoints) :: rcn,zcn
real(wp_), dimension(q_spline%ndata) :: qpsi
character(256) :: msg ! for log messages formatting 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:', & call log_info('constant ψ surfaces for:', &
mod='gray_core', proc='print_surfq') mod='gray_core', proc='print_surfq')
do i=1, size(qval) do i = 1, size(qvals)
! FIXME: check for non monotonous q profile ! Find value of ψ_n for the current q
call locate(abs(qpsi), q_spline%ndata, qval(i), i1) block
if (i1 > 0 .and. i1 < q_spline%ndata) then real(wp_) :: sol(1), fvec(1), fjac(1,1), wa(7)
call intlin(abs(qpsi(i1)), q_spline%data(i1), abs(qpsi(i1+1)), & integer :: info
q_spline%data(i1+1), qval(i), psival)
rup = rmaxis ! Note: since we are mostly interested in q=3/2,2 we start
rlw = rmaxis ! searching near the boundary in case q is not monotonic.
zup = (zbsup + zmaxis)/2.0_wp_ sol = [0.8_wp_] ! first guess
zlw = (zmaxis + zbinf)/2.0_wp_
call contours_psi(psival, rcn, zcn, rup, zup, rlw, zlw) ! Solve fq(ψ_n) = qvals(i) for ψ_n
call print_contour(psival, rcn, zcn) call hybrj1(equation, n=1, x=sol, fvec=fvec, fjac=fjac, &
rhot = frhotor(sqrt(psival)) 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))') & write (msg, '(4(x,a,"=",g0.3))') &
'q', qval(i), 'ψ', psival, 'rhop', sqrt(psival), 'rhot', rhot 'q', qvals(i), 'ψ_n', psi_n, 'ρ_p', rho_p, 'ρ_t', rho_t
call log_info(msg, mod='gray_core', proc='print_surfq') call log_info(msg, mod='gray_core', proc='print_surfq')
end if
end do 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 end subroutine print_surfq