a4a39571ce
When combining multiple independent simulations with `gray -s`, the statistics in sum-summary.txt are computed by summing the profiles pointwise and calling `postproc_profiles` with their (supposedly common) MHD equilibrium. The subroutine assumes the normalisations of dP/dV⋅dV and J_φ⋅dA are P_inside(ρ=1), I_cd_inside(ρ=1), however if the individual profiles have been computed with different versions of Gray, built with different compilers, or on different machines, the values may be very slightly off. Further, due to the use of the relation Δρ² = ⟨ρ²⟩ - ⟨ρ⟩², the profile widths are very sensible to the normalisation, producing wildly incorrect values. For example, due to the refactor of `magsurf_data`, a change of about 0.05% in dV resulted in Δρ² being off by 300%.
405 lines
11 KiB
Fortran
405 lines
11 KiB
Fortran
module pec
|
||
use const_and_precisions, only : wp_
|
||
|
||
implicit none
|
||
|
||
real(wp_), dimension(:), allocatable, save :: rhop_tab,rhot_tab
|
||
real(wp_), dimension(:), allocatable, save :: rtabpsi1
|
||
real(wp_), dimension(:), allocatable, save :: dvol,darea
|
||
real(wp_), dimension(:), allocatable, save :: ratjav,ratjbv,ratjplv
|
||
|
||
contains
|
||
|
||
subroutine pec_init(params, equil, rt_in)
|
||
use gray_params, only : output_parameters, RHO_POL
|
||
use gray_equil, only : abstract_equil
|
||
|
||
! 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
|
||
integer :: it
|
||
real(wp_) :: drt, rt, rt1, rhop1
|
||
real(wp_) :: ratjai, ratjbi, ratjpli
|
||
real(wp_) :: voli0, voli1, areai0, areai1
|
||
|
||
associate(n => params%nrho)
|
||
|
||
! rt_in present: read input grid
|
||
! else: build equidistant grid dimension n
|
||
|
||
! ipec=0 ρ_t grid
|
||
! ipec=1 ρ_p grid
|
||
call dealloc_pec
|
||
allocate(rhop_tab(n), rhot_tab(n), rtabpsi1(0:n))
|
||
allocate(dvol(n), darea(n), ratjav(n), ratjbv(n), ratjplv(n))
|
||
|
||
voli0 = 0
|
||
areai0 = 0
|
||
rtabpsi1(0) = 0
|
||
|
||
do it = 1, n
|
||
if(present(rt_in)) then
|
||
! read radial grid from input
|
||
rt = rt_in(it)
|
||
if (it < n) then
|
||
drt = rt_in(it+1)-rt
|
||
end if
|
||
else
|
||
! build equidistant radial grid
|
||
drt = 1/dble(n - 1)
|
||
rt = (it - 1)*drt
|
||
end if
|
||
! radial coordinate of i-(i+1) interval mid point
|
||
if (it < n) then
|
||
rt1 = rt + drt/2
|
||
else
|
||
rt1 = 1
|
||
end if
|
||
if (params%ipec == RHO_POL) then
|
||
rhop_tab(it) = rt
|
||
rhot_tab(it) = equil%pol2tor(rt)
|
||
rhop1 = rt1
|
||
else
|
||
rhot_tab(it) = rt
|
||
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
|
||
rtabpsi1(it) = rhop1**2
|
||
|
||
call equil%flux_average(rhop1, area=areai1, volume=voli1)
|
||
dvol(it) = abs(voli1 - voli0)
|
||
darea(it) = abs(areai1 - areai0)
|
||
voli0 = voli1
|
||
areai0 = areai1
|
||
|
||
call equil%flux_average(rhop_tab(it), ratio_astra_tor=ratjai, &
|
||
ratio_jintrac_tor=ratjbi, ratio_paral_tor=ratjpli)
|
||
ratjav(it) = ratjai
|
||
ratjbv(it) = ratjbi
|
||
ratjplv(it) = ratjpli
|
||
end do
|
||
|
||
end associate
|
||
end subroutine pec_init
|
||
|
||
|
||
subroutine spec(psjki, ppabs, ccci, iiv, pabs, &
|
||
currt, dpdv, ajphiv, ajcd, pins, currins)
|
||
! subroutine arguments
|
||
real(wp_), dimension(:, :), intent(in) :: psjki, ppabs, ccci
|
||
integer, intent(in) :: iiv(:)
|
||
real(wp_), intent(in) :: pabs, currt
|
||
real(wp_), dimension(:), intent(out) :: dpdv, ajphiv, ajcd, pins, currins
|
||
|
||
! local variables
|
||
integer :: i, ii, jk
|
||
real(wp_) :: spds, sccs, facpds, facjs
|
||
real(wp_), dimension(size(psjki, dim=2)) :: xxi, ypt, yamp
|
||
real(wp_), dimension(size(dpdv)) :: wdpdv, wajphiv
|
||
|
||
! calculation of dP and dI over radial grid
|
||
dpdv = 0
|
||
ajphiv = 0
|
||
do jk = 1, size(iiv)
|
||
ii = iiv(jk)
|
||
if (ii < size(psjki, dim=2)) then
|
||
if (psjki(jk,ii+1) /= 0) ii = ii + 1 !!! CHECK
|
||
end if
|
||
xxi = 0
|
||
ypt = 0
|
||
yamp = 0
|
||
do i = 1, ii
|
||
xxi(i)=abs(psjki(jk,i))
|
||
if(xxi(i) <= 1) then
|
||
ypt(i) = ppabs(jk,i)
|
||
yamp(i) = ccci(jk,i)
|
||
end if
|
||
end do
|
||
call pec_tab(xxi, ypt, yamp, ii, rtabpsi1, wdpdv, wajphiv)
|
||
dpdv = dpdv + wdpdv
|
||
ajphiv = ajphiv + wajphiv
|
||
end do
|
||
|
||
! here dpdv is still dP (not divided yet by dV)
|
||
! here ajphiv is still dI (not divided yet by dA)
|
||
spds = 0
|
||
sccs = 0
|
||
do i = 1, size(dpdv)
|
||
spds = spds + dpdv(i)
|
||
sccs = sccs + ajphiv(i)
|
||
pins(i) = spds
|
||
currins(i) = sccs
|
||
end do
|
||
|
||
facpds = 1
|
||
facjs = 1
|
||
if(spds > 0) facpds = pabs / spds
|
||
if(sccs /= 0) facjs = currt / sccs
|
||
|
||
dpdv = facpds * (dpdv / dvol)
|
||
ajphiv = facjs * (ajphiv / darea)
|
||
ajcd = ajphiv * ratjbv
|
||
pins = facpds * pins
|
||
currins = facjs * currins
|
||
|
||
! now dpdv is dP/dV [MW/m^3]
|
||
! now ajphiv is J_phi=dI/dA [MA/m^2]
|
||
end subroutine spec
|
||
|
||
|
||
subroutine pec_tab(xxi, ypt, yamp, ii, xtab1, wdpdv, wajphiv)
|
||
! Power and current projected on ψ grid - mid points
|
||
use utils, only : locate, linear_interp
|
||
|
||
! subroutine arguments
|
||
real(wp_), dimension(:), intent(in) :: xxi, ypt, yamp
|
||
integer, intent(in) :: ii
|
||
real(wp_), dimension(0:), intent(in) :: xtab1
|
||
real(wp_), dimension(ubound(xtab1, 1)), intent(out) :: wdpdv, wajphiv
|
||
|
||
! local variables
|
||
integer :: isev(21)
|
||
real(wp_) :: ppa1, ppa2, cci1, cci2, dppa, didst, rt1
|
||
integer :: i, is, ise0, idecr, iise0, iise, iis, iis1
|
||
integer :: ind1, ind2, iind, ind, indi, itb1
|
||
|
||
isev = 0
|
||
ise0 = 0
|
||
idecr = -1
|
||
is = 1
|
||
wdpdv = 0
|
||
wajphiv = 0
|
||
do i=1,ii
|
||
if(ise0 == 0) then
|
||
if(xxi(i) < 1) then
|
||
ise0 = i
|
||
isev(is) = max(i - 1, 1)
|
||
is = is + 1
|
||
end if
|
||
else
|
||
if (idecr == -1) then
|
||
if(xxi(i) > xxi(i-1)) then
|
||
isev(is) = i - 1
|
||
is = is + 1
|
||
idecr = 1
|
||
end if
|
||
else
|
||
if(xxi(i) > 1) exit
|
||
if(xxi(i) < xxi(i-1)) then
|
||
isev(is) = i - 1
|
||
is = is + 1
|
||
idecr = -1
|
||
end if
|
||
end if
|
||
end if
|
||
end do
|
||
|
||
isev(is) = i-1
|
||
ppa1 = 0
|
||
cci1 = 0
|
||
|
||
do iis = 1, is-1
|
||
iis1 = iis + 1
|
||
iise0 = isev(iis)
|
||
iise = isev(iis1)
|
||
if (mod(iis,2) /= 0) then
|
||
idecr = -1
|
||
ind1 = size(wdpdv)
|
||
ind2 = 2
|
||
iind = -1
|
||
else
|
||
idecr = 1
|
||
ind1 = 1
|
||
ind2 = size(wdpdv)
|
||
iind = 1
|
||
end if
|
||
do ind = ind1, ind2, iind
|
||
indi = ind
|
||
if (idecr == -1) indi = ind - 1
|
||
rt1 = xtab1(indi)
|
||
call locate(xxi(iise0:iise), rt1, itb1)
|
||
if (itb1 >= 1 .and. itb1 <= iise) then
|
||
itb1 = itb1 + (iise0 - 1)
|
||
if(itb1 == iise) then
|
||
ppa2=ypt(itb1)
|
||
cci2=yamp(itb1)
|
||
else
|
||
ppa2 = linear_interp(xxi(itb1), ypt(itb1), xxi(itb1+1), ypt(itb1+1), rt1)
|
||
cci2 = linear_interp(xxi(itb1), yamp(itb1), xxi(itb1+1), yamp(itb1+1), rt1)
|
||
end if
|
||
dppa = ppa2 - ppa1
|
||
didst = cci2 - cci1
|
||
wdpdv(ind) = wdpdv(ind) + dppa
|
||
wajphiv(ind) = wajphiv(ind) + didst
|
||
ppa1 = ppa2
|
||
cci1 = cci2
|
||
end if
|
||
end do
|
||
end do
|
||
end subroutine pec_tab
|
||
|
||
|
||
subroutine postproc_profiles(equil, 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 gray_equil, only : abstract_equil
|
||
|
||
class(abstract_equil), intent(in) :: equil
|
||
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_) :: P_tot, I_tot, ratjplmx, rhopjava, rhoppav
|
||
real(wp_) :: rhotjav, dvdrhotav, dadrhotava
|
||
|
||
! ⟨ρ_t⟩, Δρ_t using dP/dV⋅dV as a PDF
|
||
! Note: the factor √8 is to match the FWHM of a Gaussian
|
||
P_tot = sum(dpdv*dvol)
|
||
if (P_tot > 0) then
|
||
rhotpav = sum(rhot_tab * dpdv*dvol)/P_tot
|
||
drhotpav = sqrt(8 * sum((rhot_tab - rhotpav)**2 * dpdv*dvol)/P_tot)
|
||
else
|
||
rhotpav = 0
|
||
drhotpav = 0
|
||
end if
|
||
|
||
! ⟨ρ_t⟩ using J_φ⋅dA as a PDF
|
||
I_tot = sum(ajphiv*darea)
|
||
if (abs(I_tot) > 0) then
|
||
rhotjav = sum(rhot_tab*ajphiv*darea)/I_tot
|
||
else
|
||
rhotjav = 0
|
||
end if
|
||
|
||
! ⟨ρ_t⟩, Δρ_t using |J_φ|⋅dA as a PDF
|
||
I_tot = sum(abs(ajphiv)*darea)
|
||
if (I_tot > 0) then
|
||
rhotjava = sum(rhot_tab * abs(ajphiv)*darea)/I_tot
|
||
drhotjava = sqrt(8* sum((rhot_tab - rhotjava)**2 * abs(ajphiv)*darea)/I_tot)
|
||
else
|
||
rhotjava = 0
|
||
drhotjava = 0
|
||
end if
|
||
|
||
rhoppav = equil%tor2pol(rhotpav)
|
||
rhopjava = equil%tor2pol(rhotjava)
|
||
|
||
if (P_tot > 0) then
|
||
call equil%flux_average(rhoppav, dVdrho_t=dvdrhotav)
|
||
dpdvp = P_tot * 2/(sqrt(pi)*drhotpav*dvdrhotav)
|
||
call profwidth(rhot_tab,dpdv,rhotp,dpdvmx,drhotp)
|
||
else
|
||
dpdvp = 0
|
||
rhotp = 0
|
||
dpdvmx = 0
|
||
drhotp = 0
|
||
end if
|
||
|
||
if (I_tot > 0) then
|
||
call equil%flux_average(rhopjava, dAdrho_t=dadrhotava, ratio_astra_tor=ratjamx, &
|
||
ratio_jintrac_tor=ratjbmx, ratio_paral_tor=ratjplmx)
|
||
ajphip = I_tot * 2/(sqrt(pi)*drhotjava*dadrhotava)
|
||
call profwidth(rhot_tab,ajphiv,rhotjfi,ajmxfi,drhotjfi)
|
||
else
|
||
ajphip = 0
|
||
rhotjfi = 0
|
||
ajmxfi = 0
|
||
drhotjfi = 0
|
||
ratjamx = 0
|
||
ratjbmx = 0
|
||
end if
|
||
end subroutine postproc_profiles
|
||
|
||
|
||
subroutine profwidth(xx,yy,xpk,ypk,dxxe)
|
||
use const_and_precisions, only : emn1
|
||
use utils, only : locate, linear_interp, vmaxmin
|
||
|
||
! subroutine arguments
|
||
real(wp_), intent(in) :: xx(:), yy(:)
|
||
real(wp_), intent(out) :: xpk, ypk, dxxe
|
||
|
||
! local variables
|
||
integer :: imn,imx,ipk,ie
|
||
real(wp_) :: xmn,xmx,ymn,ymx,xpkp,xpkm,yye,rte1,rte2
|
||
real(wp_) :: ypkp,ypkm
|
||
|
||
call vmaxmin(yy, ymn, ymx, imn, imx)
|
||
ypk = 0
|
||
xmx = xx(imx)
|
||
xmn = xx(imn)
|
||
if (abs(ymx) > abs(ymn)) then
|
||
ipk = imx
|
||
ypkp = ymx
|
||
xpkp = xmx
|
||
if(abs(ymn/ymx) < 1.0e-2_wp_) ymn = 0.0_wp_
|
||
ypkm = ymn
|
||
xpkm = xmn
|
||
else
|
||
ipk = imn
|
||
ypkp = ymn
|
||
xpkp = xmn
|
||
if(abs(ymx/ymn) < 1.0e-2_wp_) ymx = 0.0_wp_
|
||
ypkm = ymx
|
||
xpkm = xmx
|
||
end if
|
||
if(xpkp > 0) then
|
||
xpk = xpkp
|
||
ypk = ypkp
|
||
yye = ypk*emn1
|
||
call locate(yy(1:ipk), yye, ie)
|
||
if(ie > 0 .and. ie < size(yy)) then
|
||
rte1 = linear_interp(yy(ie), xx(ie), yy(ie+1), xx(ie+1), yye)
|
||
else
|
||
rte1 = 0
|
||
end if
|
||
call locate(yy(ipk:size(yy)), yye, ie)
|
||
if(ie > 0 .and. ie < size(yy)) then
|
||
ie = ie + (ipk - 1)
|
||
rte2 = linear_interp(yy(ie), xx(ie), yy(ie+1), xx(ie+1), yye)
|
||
else
|
||
rte2 = 0
|
||
end if
|
||
else
|
||
ipk=2
|
||
xpk=xx(2)
|
||
ypk=yy(2)
|
||
rte1=0.0_wp_
|
||
yye=ypk*emn1
|
||
call locate(yy, yye, ie)
|
||
if(ie > 0 .and. ie < size(yy)) then
|
||
rte2 = linear_interp(yy(ie), xx(ie), yy(ie+1), xx(ie+1), yye)
|
||
else
|
||
rte2 = 0
|
||
end if
|
||
end if
|
||
dxxe = rte2 - rte1
|
||
if(ymx /= 0 .and. ymn /= 0) dxxe = -dxxe
|
||
end subroutine profwidth
|
||
|
||
|
||
subroutine dealloc_pec
|
||
if (allocated(rhop_tab)) deallocate(rhop_tab)
|
||
if (allocated(rhot_tab)) deallocate(rhot_tab)
|
||
if (allocated(rtabpsi1)) deallocate(rtabpsi1)
|
||
if (allocated(dvol)) deallocate(dvol)
|
||
if (allocated(darea)) deallocate(darea)
|
||
if (allocated(ratjav)) deallocate(ratjav)
|
||
if (allocated(ratjbv)) deallocate(ratjbv)
|
||
if (allocated(ratjplv)) deallocate(ratjplv)
|
||
end subroutine dealloc_pec
|
||
|
||
end module pec
|