src/pec.f90: remove leftover file
This commit is contained in:
Michele Guerini Rocco
parent
1a8011f2b7
commit
a8a35bc783
GPG Key ID:
BFBAF4C975F76450
335
src/pec.f90
335
src/pec.f90
@ -1,335 +0,0 @@
|
||||
! This module provides the ray_projector object that is used to
|
||||
! compute the final ECCD power and current density radial profiles
|
||||
module gray_ec_profiles
|
||||
|
||||
use const_and_precisions, only : wp_
|
||||
|
||||
implicit none
|
||||
|
||||
type ray_projector
|
||||
! Projects the volumetric rays data onto a radial grid to
|
||||
! obtain the ECCD power and current density profiles
|
||||
real(wp_), public, allocatable :: rho_p(:) ! radial bin edges (ρ_p)
|
||||
real(wp_), public, allocatable :: rho_t(:) ! radial bin edges (ρ_t)
|
||||
real(wp_), private, allocatable :: psi_mid(:) ! radial bin midpoints (ψ_n)
|
||||
real(wp_), private, allocatable :: dA(:), dV(:) ! area, volume elements
|
||||
real(wp_), private, allocatable :: ratio_J_cd_phi(:) ! J_cd_jintrac / J_φ
|
||||
contains
|
||||
procedure :: init ! Initialises the grids etc.
|
||||
procedure :: project ! Projects rays to radial profile
|
||||
procedure :: statistics ! Compute statistics of a profile
|
||||
end type
|
||||
|
||||
private
|
||||
|
||||
contains
|
||||
|
||||
subroutine init(self, params, equil)
|
||||
! Initialises the radial grid and other quantities needed to
|
||||
! compute the radial profiles (dP/dV, J_φ, J_cd)
|
||||
use const_and_precisions, only : one
|
||||
use gray_params, only : output_parameters, RHO_POL
|
||||
use gray_equil, only : abstract_equil
|
||||
|
||||
! subroutine arguments
|
||||
class(ray_projector), intent(inout) :: self
|
||||
type(output_parameters), intent(in) :: params
|
||||
class(abstract_equil), intent(in) :: equil
|
||||
|
||||
! local variables
|
||||
integer :: i
|
||||
real(wp_) :: drho, rho, rho_mid
|
||||
real(wp_) :: V0, V1, A0, A1
|
||||
|
||||
associate (n => params%nrho)
|
||||
|
||||
allocate(self%rho_p(n), self%rho_t(n), self%psi_mid(n))
|
||||
allocate(self%dV(n), self%dA(n), self%ratio_J_cd_phi(n))
|
||||
|
||||
V0 = 0
|
||||
A0 = 0
|
||||
|
||||
do i = 1, n
|
||||
if(allocated(params%grid)) then
|
||||
! Read points from the input grid
|
||||
rho = params%grid(i)
|
||||
if (i < n) drho = params%grid(i + 1) - rho
|
||||
else
|
||||
! Build a uniform grid
|
||||
drho = one/(n - 1)
|
||||
rho = (i - 1)*drho
|
||||
end if
|
||||
|
||||
! Switch between uniform grid in ρ_p or ρ_t
|
||||
if (params%ipec == RHO_POL) then
|
||||
self%rho_p(i) = rho
|
||||
self%rho_t(i) = equil%pol2tor(rho)
|
||||
rho_mid = rho + drho/2
|
||||
else
|
||||
self%rho_t(i) = rho
|
||||
self%rho_p(i) = equil%tor2pol(rho)
|
||||
rho_mid = equil%tor2pol(rho + drho/2)
|
||||
end if
|
||||
|
||||
! ψ_n at the midpoint of the bin
|
||||
!
|
||||
! Note: each bin is associated with a volume element ΔV₁=V₁-V₀,
|
||||
! with V₁,₀ computed at the bin midpoint. For consistency the ΔP
|
||||
! is also computed at the midpoint: ΔP₁ = P_in(ψ₁) - P_in(ψ₀),
|
||||
! where ψ₁,₀ are shifted by Δρ/2 w.r.t ρ₁.
|
||||
self%psi_mid(i) = rho_mid**2
|
||||
|
||||
! Compute ΔV, ΔA for each bin
|
||||
call equil%flux_average(rho_mid, area=A1, volume=V1)
|
||||
self%dA(i) = abs(A1 - A0)
|
||||
self%dV(i) = abs(V1 - V0)
|
||||
V0 = V1
|
||||
A0 = A1
|
||||
|
||||
! Compute conversion factors between J_φ and J_cd
|
||||
call equil%flux_average(self%rho_p(i), ratio_jintrac_tor=self%ratio_J_cd_phi(i))
|
||||
end do
|
||||
|
||||
end associate
|
||||
end subroutine init
|
||||
|
||||
|
||||
subroutine project(self, psi_n, P_abs, I_cd, last_step, &
|
||||
dPdV, J_phi, J_cd, P_in, I_in)
|
||||
! Projects the volumetric rays data ψ_n, P_abs, I_cd onto a radial grid
|
||||
! to obtain the density profiles dP/dV, J_φ, J_cd and the cumulatives
|
||||
! P_in and I_in
|
||||
use const_and_precisions, only : zero
|
||||
use utils, only : locate_unord, linear_interp
|
||||
|
||||
! subroutine arguments
|
||||
class(ray_projector), intent(in) :: self
|
||||
real(wp_), intent(in) :: psi_n(:, :) ! normalised flux (ray,step)
|
||||
real(wp_), intent(in) :: P_abs(:, :) ! absorbed power (ray,step)
|
||||
real(wp_), intent(in) :: I_cd(:, :) ! driven current (ray,step)
|
||||
integer, intent(in) :: last_step(:) ! last step absorption took place
|
||||
real(wp_), intent(out) :: dPdV(:), J_phi(:) ! power, toroidal current density
|
||||
real(wp_), intent(out) :: J_cd(:) ! current drive density, JINTRAC def.
|
||||
real(wp_), intent(out) :: P_in(:), I_in(:) ! power, current within ρ
|
||||
|
||||
! local variables
|
||||
integer :: i, ray, bin
|
||||
integer :: cross(maxval(last_step)+1), ncross
|
||||
real(wp_) :: I_entry, I_exit, P_entry, P_exit
|
||||
|
||||
! We start by computing the cumulatives first, and obtain the densities
|
||||
! as finite differences later. This should avoid precision losses.
|
||||
P_in = 0
|
||||
I_in = 0
|
||||
|
||||
! For each bin in the radial grid
|
||||
bins: do concurrent (bin = 1 : size(self%rho_t))
|
||||
! For each ray
|
||||
rays: do ray = 1, size(psi_n, dim=1)
|
||||
! Skip rays with no absorption
|
||||
if (last_step(ray) == 1) cycle
|
||||
|
||||
! Find where the ray crosses the flux surface of the current bin
|
||||
call locate_unord(psi_n(ray, :), self%psi_mid(bin), cross, last_step(ray), ncross)
|
||||
|
||||
if (psi_n(ray, 1) < 0) then
|
||||
! Skip unphysical crossing from -1 (ψ undefined)
|
||||
cross = cross(2:ncross)
|
||||
ncross = ncross - 1
|
||||
else if (psi_n(ray, 1) < self%psi_mid(bin)) then
|
||||
! Add a crossing when starting within ψ (reflections!)
|
||||
cross = [1, cross(1:ncross)]
|
||||
end if
|
||||
|
||||
! For each crossing
|
||||
crossings: do i = 1, ncross, 2
|
||||
! Interpolate P_abs(ψ) and I_cd(ψ) at the entry
|
||||
P_entry = interp(psi_n(ray,:), P_abs(ray,:), cross(i), self%psi_mid(bin))
|
||||
I_entry = interp(psi_n(ray,:), I_cd(ray,:), cross(i), self%psi_mid(bin))
|
||||
|
||||
if (i+1 <= ncross) then
|
||||
! Interpolate P_abs(ψ) and I_cd(ψ) at the exit
|
||||
P_exit = interp(psi_n(ray,:), P_abs(ray,:), cross(i+1), self%psi_mid(bin))
|
||||
I_exit = interp(psi_n(ray,:), I_cd(ray,:), cross(i+1), self%psi_mid(bin))
|
||||
else
|
||||
! The ray never left, use last_step as exit point
|
||||
P_exit = P_abs(ray, last_step(ray))
|
||||
I_exit = I_cd(ray, last_step(ray))
|
||||
end if
|
||||
|
||||
! Accumulate the deposition from this crossing
|
||||
P_in(bin) = P_in(bin) + (P_exit - P_entry)
|
||||
I_in(bin) = I_in(bin) + (I_exit - I_entry)
|
||||
end do crossings
|
||||
end do rays
|
||||
end do bins
|
||||
|
||||
! Compute densities
|
||||
dPdV(1) = P_in(1) / self%dV(1)
|
||||
J_phi(1) = I_in(1) / self%dA(1)
|
||||
|
||||
do concurrent (bin = 2 : size(self%rho_t))
|
||||
dPdV(bin) = max(P_in(bin) - P_in(bin-1), zero) / self%dV(bin)
|
||||
J_phi(bin) = (I_in(bin) - I_in(bin-1)) / self%dA(bin)
|
||||
end do
|
||||
|
||||
J_cd = J_phi * self%ratio_J_cd_phi
|
||||
|
||||
contains
|
||||
|
||||
pure real(wp_) function interp(x, y, i, x1)
|
||||
! Interpolate linearly y(x) at x=x(i)
|
||||
real(wp_), intent(in) :: x(:), y(:), x1
|
||||
integer, intent(in) :: i
|
||||
interp = linear_interp(x(i), y(i), x(i+1), y(i+1), x1)
|
||||
end function interp
|
||||
|
||||
end subroutine project
|
||||
|
||||
|
||||
subroutine statistics(self, equil, dPdV, J_phi, rho_avg_P, drho_avg_P, &
|
||||
rho_avg_J, drho_avg_J, dPdV_peak, J_phi_peak, &
|
||||
rho_max_P, drho_max_P, rho_max_J, drho_max_J, &
|
||||
dPdV_max, J_phi_max, ratio_Ja_max, ratio_Jb_max)
|
||||
! Given the radial profiles dP/dV and J_φ computes the following statics:
|
||||
|
||||
use const_and_precisions, only : pi
|
||||
use gray_equil, only : abstract_equil
|
||||
|
||||
class(ray_projector), intent(in) :: self
|
||||
class(abstract_equil), intent(in) :: equil
|
||||
|
||||
real(wp_), intent(in) :: dPdV(:), J_phi(:) ! dP/dV(ρ), J_φ(ρ) profiles
|
||||
real(wp_), intent(out) :: rho_avg_P, drho_avg_P, dPdV_peak ! ⟨ρ_t⟩, Δρ_t, max dP/dV using dP/dV⋅dV as a PDF
|
||||
real(wp_), intent(out) :: rho_avg_J, drho_avg_J, J_phi_peak ! ⟨ρ_t⟩, Δρ_t, max |J_φ| using |J_φ|⋅dA as a PDF
|
||||
real(wp_), intent(out) :: rho_max_P, drho_max_P, dPdV_max ! argmax, FWHM and max of dPdV
|
||||
real(wp_), intent(out) :: rho_max_J, drho_max_J, J_phi_max ! argmax, FWHM and max of |J_φ|
|
||||
real(wp_), intent(out) :: ratio_Ja_max, ratio_Jb_max ! max of J_cd_astra/J_φ, J_cd_jintrac/J_φ
|
||||
|
||||
! local variables
|
||||
real(wp_) :: P_tot, I_tot
|
||||
real(wp_) :: dVdrho_t, dAdrho_t
|
||||
|
||||
! ⟨ρ_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*self%dV)
|
||||
if (P_tot > 0) then
|
||||
rho_avg_P = sum(self%rho_t * dPdV*self%dV)/P_tot
|
||||
drho_avg_P = sqrt(8 * sum((self%rho_t - rho_avg_P)**2 * dPdV*self%dV)/P_tot)
|
||||
else
|
||||
rho_avg_P = 0
|
||||
drho_avg_P = 0
|
||||
end if
|
||||
|
||||
! ⟨ρ_t⟩, Δρ_t using |J_φ|⋅dA as a PDF
|
||||
I_tot = sum(abs(J_phi)*self%dA)
|
||||
if (I_tot > 0) then
|
||||
rho_avg_J = sum(self%rho_t * abs(J_phi)*self%dA)/I_tot
|
||||
drho_avg_J = sqrt(8* sum((self%rho_t - rho_avg_J)**2 * abs(J_phi)*self%dA)/I_tot)
|
||||
else
|
||||
rho_avg_J = 0
|
||||
drho_avg_J = 0
|
||||
end if
|
||||
|
||||
if (P_tot > 0) then
|
||||
! Compute max of dP/dV as if it were a Gaussian
|
||||
call equil%flux_average(equil%tor2pol(rho_avg_P), dVdrho_t=dVdrho_t)
|
||||
dPdV_peak = P_tot * 2/(sqrt(pi) * drho_avg_P * dVdrho_t)
|
||||
|
||||
! Compute numerically argmax, max and FWHM of dP/dV
|
||||
call profwidth(self%rho_t, dPdV, rho_max_P, dPdV_max, drho_max_P)
|
||||
else
|
||||
dPdV_peak = 0
|
||||
rho_max_P = 0
|
||||
dPdV_max = 0
|
||||
drho_max_P = 0
|
||||
end if
|
||||
|
||||
if (I_tot > 0) then
|
||||
! Compute max of J_φ as if it were a Gaussian
|
||||
call equil%flux_average(equil%tor2pol(rho_avg_J), dAdrho_t=dAdrho_t, &
|
||||
ratio_astra_tor=ratio_Ja_max, ratio_jintrac_tor=ratio_Jb_max)
|
||||
J_phi_peak = I_tot * 2/(sqrt(pi)*drho_avg_J*dAdrho_t)
|
||||
|
||||
! Compute numerically argmax, max and FWHM of dP/dV
|
||||
call profwidth(self%rho_t, J_phi, rho_max_J, J_phi_max, drho_max_J)
|
||||
else
|
||||
J_phi_peak = 0
|
||||
rho_max_J = 0
|
||||
J_phi_max = 0
|
||||
drho_max_J = 0
|
||||
ratio_Ja_max = 0
|
||||
ratio_Jb_max = 0
|
||||
end if
|
||||
end subroutine statistics
|
||||
|
||||
|
||||
subroutine profwidth(x, y, x_max, y_max, fwhm)
|
||||
! Computes the argmax, max, and FWHM of the curve y(x)
|
||||
|
||||
use const_and_precisions, only : emn1
|
||||
use utils, only : locate, linear_interp, vmaxmin
|
||||
|
||||
! subroutine arguments
|
||||
real(wp_), intent(in) :: x(:), y(:)
|
||||
real(wp_), intent(out) :: x_max, y_max, fwhm
|
||||
|
||||
! local variables
|
||||
integer :: imn, imx, ipk, ie
|
||||
real(wp_) :: xmn, xmx, ymn, ymx, xpkp, xpkm, yye, rte1, rte2
|
||||
real(wp_) :: ypkp, ypkm
|
||||
|
||||
call vmaxmin(y, ymn, ymx, imn, imx)
|
||||
y_max = 0
|
||||
xmx = x(imx)
|
||||
xmn = x(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
|
||||
x_max = xpkp
|
||||
y_max = ypkp
|
||||
yye = y_max*emn1
|
||||
call locate(y(1:ipk), yye, ie)
|
||||
if(ie > 0 .and. ie < size(y)) then
|
||||
rte1 = linear_interp(y(ie), x(ie), y(ie+1), x(ie+1), yye)
|
||||
else
|
||||
rte1 = 0
|
||||
end if
|
||||
call locate(y(ipk:size(y)), yye, ie)
|
||||
if(ie > 0 .and. ie < size(y)) then
|
||||
ie = ie + (ipk - 1)
|
||||
rte2 = linear_interp(y(ie), x(ie), y(ie+1), x(ie+1), yye)
|
||||
else
|
||||
rte2 = 0
|
||||
end if
|
||||
else
|
||||
ipk=2
|
||||
x_max=x(2)
|
||||
y_max=y(2)
|
||||
rte1=0.0_wp_
|
||||
yye=y_max*emn1
|
||||
call locate(y, yye, ie)
|
||||
if(ie > 0 .and. ie < size(y)) then
|
||||
rte2 = linear_interp(y(ie), x(ie), y(ie+1), x(ie+1), yye)
|
||||
else
|
||||
rte2 = 0
|
||||
end if
|
||||
end if
|
||||
fwhm = rte2 - rte1
|
||||
if(ymx /= 0 .and. ymn /= 0) fwhm = -fwhm
|
||||
end subroutine profwidth
|
||||
|
||||
end module gray_ec_profiles
|
||||
Loading…
Reference in New Issue
Block a user