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

replace pec module with an object

Browse Source 
This was the final module with global variables to be rewritten.
The functionaly of pec: `pec_init`, `spec`, `postproc_profiles` has been
replaced by the `ray_projector` object in `gray_project.f90` with the
following methods: `projector%init`, `projector%project` and
`projector%statistics`.

The new code is functionally identically with only breaking change being
in Δρ_J, the full-width at max/e of the current density.
Before this change Δρ_J could be negative to signal the J_φ profile had
at least one positive and one negative peak, after the value is always
positive. Note: in either case Δρ_J was given by the largest peak only.
This commit is contained in:
Michele Guerini Rocco 2024-10-30 17:43:08 +01:00
parent 12f15239df
commit 6f66317541
Signed by: rnhmjoj
GPG Key ID: BFBAF4C975F76450
8 changed files with 565 additions and 330 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

28
src/gray_core.f90
View File

@ -7,13 +7,14 @@ module gray_core
contains
subroutine gray_main(params, equil, plasma, limiter, results, error, rhout)
subroutine gray_main(params, equil, plasma, limiter, results, error)
use const_and_precisions, only : zero, one, comp_tiny
use polarization, only : ellipse_to_field
use types, only : contour, table, wrap
use gray_params, only : gray_parameters, gray_results, EQ_VACUUM, ABSORP_OFF
use gray_equil, only : abstract_equil
use gray_plasma, only : abstract_plasma
use gray_project, only : ray_projector
use gray_tables, only : init_tables, store_ec_profiles, store_ray_data, &
store_beam_shape, find_flux_surfaces, &
flux_surfaces, kinetic_profiles, flux_averages, &
@ -22,8 +23,6 @@ contains
print_err_raytracing, print_err_ecrh_cd
use beams, only : xgygcoeff, launchangles2n
use beamdata, only : pweight
use pec, only : pec_init, spec, postproc_profiles, dealloc_pec, &
rhop_tab, rhot_tab
use utils, only : vmaxmin
use multipass, only : initbeam, initmultipass, turnoffray, &
plasma_in, plasma_out, wall_out
@ -37,9 +36,6 @@ contains
type(gray_results), intent(out) :: results
integer(kind=gray_error), intent(out) :: error
! Predefined grid for the output profiles (optional)
real(wp_), dimension(:), intent(in), optional :: rhout
! local variables
real(wp_), parameter :: taucr = 12._wp_, etaucr = exp(-taucr)
character, dimension(2), parameter :: mode=['O','X']
@ -67,6 +63,8 @@ contains
! buffer for formatting log messages
character(256) :: msg
type(ray_projector) :: projector
real(wp_), dimension(2) :: pabs_pass,icd_pass,cpl
real(wp_), dimension(3) :: xv,anv0,anv,bv,derxg
@ -120,7 +118,7 @@ contains
if (params%equilibrium%iequil /= EQ_VACUUM) then
! Initialise the output profiles
call pec_init(params%output, equil, rhout)
call projector%init(params%output, equil)
end if
! Allocate memory for the results...
@ -561,9 +559,8 @@ contains
if (params%equilibrium%iequil /= EQ_VACUUM) then
! compute power and current density profiles for all rays
call spec(psjki, ppabs, ccci, iiv, pabs_beam, &
icd_beam, dpdv_beam, jphi_beam, jcd_beam, &
pins_beam, currins_beam)
call projector%project(psjki, ppabs, ccci, iiv, dpdv_beam, &
jphi_beam, jcd_beam, pins_beam, currins_beam)
end if
pabs_pass(iox) = pabs_pass(iox) + pabs_beam ! 0D results for current pass, sum on O/X mode beams
@ -611,11 +608,11 @@ contains
if (params%equilibrium%iequil /= EQ_VACUUM) then
call store_ec_profiles( &
results%tables%ec_profiles, rhop_tab, rhot_tab, jphi_beam, &
jcd_beam, dpdv_beam, currins_beam, pins_beam, index_rt)
results%tables%ec_profiles, projector%rho_p, projector%rho_t, &
jphi_beam, jcd_beam, dpdv_beam, currins_beam, pins_beam, index_rt)
call postproc_profiles(equil, rhot_tab, dpdv_beam, &
jphi_beam, rhotpav, drhotpav, rhotjava, drhotjava, dpdvp, jphip, &
call projector%statistics(equil, 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
@ -657,9 +654,6 @@ contains
call log_debug('pass loop end', mod='gray_core', proc='gray_main')
! ============ main loop END ============
! Free memory
call dealloc_pec
end block
end associate
end subroutine gray_main

3
src/gray_jetto1beam.f90
View File

@ -49,6 +49,7 @@ subroutine gray_jetto1beam(ijetto, mr, mz, r, z, psin, psia, rax, zax, &
params%profiles%iprof = 1
params%output%ipec = 1
params%output%nrho = nrho
params%output%grid = sqrt(psrad)
end block init_params
! Set a simple limiter following the boundary of the data grid
@ -111,7 +112,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, equil, plasma, limiter, res, err, rhout=sqrt(psrad))
call gray_main(params, equil, plasma, limiter, res, err)
write_debug_files: block
integer :: i, err

9
src/gray_params.f90
View File

@ -155,10 +155,11 @@ module gray_params
! Output data parameters
type output_parameters
integer :: ipec = RHO_TOR ! Radial coordinate of the EC profiles
integer :: nrho = 501 ! Number of grid steps for EC profiles + 1
integer :: istpr = 5 ! Subsampling factor for beam cross-section (tables 8, 12)
integer :: istpl = 5 ! Subsampling factor for outer rays (table 33)
real(wp_), allocatable :: grid(:) ! Radial grid, defaults to uniform ρ ∈ [0, 1]
integer :: ipec = RHO_TOR ! Radial coordinate of the EC profiles
integer :: nrho = 501 ! Number of grid steps for EC profiles + 1
integer :: istpr = 5 ! Subsampling factor for beam cross-section (tables 8, 12)
integer :: istpl = 5 ! Subsampling factor for outer rays (table 33)
end type
! Other parameters

310
src/gray_project.f90 Normal file
View File

@ -0,0 +1,310 @@
! This module provides the ray_projector object that is used to
! compute the final ECCD power and current density radial profiles
module gray_project
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
public ray_projector
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, full-width at max/e and max of dPdV
real(wp_), intent(out) :: rho_max_J, drho_max_J, J_phi_max ! argmax, full-width at max/e 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 full-width at max/e 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 full-width at max/e 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 full-width at max/e of J_φ
call profwidth(self%rho_t, abs(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, width)
! Computes the argmax, max and full-width at max/e of the curve y(x)
use const_and_precisions, only : emn1
use utils, only : locate_unord, linear_interp
! subroutine arguments
real(wp_), intent(in) :: x(:), y(:)
real(wp_), intent(out) :: x_max, y_max, width
! local variables
integer :: i, i_max, left, right
integer :: nsols, sols(size(y))
real(wp_) :: x_left, x_right
! Find maximum
i_max = maxloc(y, dim=1)
x_max = x(i_max)
y_max = y(i_max)
! Find width
width = 0
! Find solutions to y(x) = y_max⋅e⁻¹
call locate_unord(y, y_max*emn1, sols, size(y), nsols)
if (nsols < 2) return
! Find the two closest solutions to the maximum
left = sols(1)
right = sols(nsols)
do i = 1, nsols
if (sols(i) > left .and. sols(i) < i_max) left = sols(i)
if (sols(i) < right .and. sols(i) > i_max) right = sols(i)
end do
! Interpolate linearly x(y) for better accuracy
x_left = linear_interp(y(left), x(left), y(left+1), x(left+1), y_max*emn1)
x_right = linear_interp(y(right), x(right), y(right+1), x(right+1), y_max*emn1)
width = x_right - x_left
end subroutine profwidth
end module gray_project

2
src/gray_tables.f90
View File

@ -271,7 +271,7 @@ contains
integer :: j, k, n
integer :: n_conts, n_points(10), n_tot
real(wp_) :: dr, dz, B_min, B_max
real(wp_) :: B, B_R, B_z, B_phi, B_tot(n_grid, n_grid)
real(wp_) :: B, B_tot(n_grid, n_grid)
real(wp_), dimension(n_cont) :: R_cont, z_cont
real(wp_), dimension(n_grid) :: R, z

14
src/main.f90
View File

@ -238,8 +238,7 @@ contains
! (of different beams) and recomputes the summary table
use gray_params, only : gray_parameters
use gray_equil, only : abstract_equil
use pec, only : pec_init, postproc_profiles, dealloc_pec, &
rhot_tab
use gray_project, only : ray_projector
! subroutine arguments
type(gray_parameters), intent(inout) :: params
@ -263,7 +262,8 @@ contains
integer, allocatable :: beams(:)
! Initialise the output profiles
call pec_init(params%output, equil)
type(ray_projector) :: projector
call projector%init(params%output, equil)
associate(nrho =>params%output%nrho)
allocate(jphi(nrho), currins(nrho), pins(nrho), rtin(nrho), rpin(nrho))
@ -357,9 +357,9 @@ contains
results%icd = currins(params%output%nrho)
! Recompute the summary values
call postproc_profiles(equil, &
rhot_tab, results%dPdV, jphi, rhotpav, drhotpav, rhotjava, &
drhotjava, dpdvp, jphip, rhotp, drhotp, rhotj, drhotj, &
call projector%statistics(equil, &
results%dPdV, jphi, rhotpav, drhotpav, rhotjava, &
drhotjava, dpdvp, jphip, rhotp, drhotp, rhotj, drhotj, &
dpdvmx, jphimx, ratjamx, ratjbmx)
write(fmtstr, '("(",i0,"(1x,e12.5),i5,7(1x,e12.5))")') nextracols+12
@ -378,8 +378,6 @@ contains
close(beams(i))
end do
! Free memory
call dealloc_pec
end subroutine sum_profiles
end program main

527
src/pec.f90
View File

@ -1,345 +1,288 @@
module pec
! 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
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
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 pec_init(params, equil, rt_in)
use gray_params, only : output_parameters, RHO_POL
use gray_equil, only : abstract_equil
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
type(output_parameters), intent(in) :: params
class(abstract_equil), intent(in) :: equil
real(wp_), intent(in), optional :: rt_in(params%nrho)
class(ray_projector), intent(inout) :: self
type(output_parameters), intent(in) :: params
class(abstract_equil), intent(in) :: equil
! local variables
integer :: it
real(wp_) :: drt, rt, rt1, rhop1
real(wp_) :: ratjai, ratjbi, ratjpli
real(wp_) :: voli0, voli1, areai0, areai1
integer :: i
real(wp_) :: drho, rho, rho_mid
real(wp_) :: V0, V1, A0, A1
associate(n => params%nrho)
associate (n => params%nrho)
! rt_in present: read input grid
! else: build equidistant grid dimension n
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))
! 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))
V0 = 0
A0 = 0
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
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 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
! 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
rhop_tab(it) = rt
rhot_tab(it) = equil%pol2tor(rt)
rhop1 = rt1
self%rho_p(i) = rho
self%rho_t(i) = equil%pol2tor(rho)
rho_mid = rho + drho/2
else
rhot_tab(it) = rt
rhop_tab(it) = equil%tor2pol(rt)
rhop1 = equil%tor2pol(rt1)
self%rho_t(i) = rho
self%rho_p(i) = equil%tor2pol(rho)
rho_mid = equil%tor2pol(rho + drho/2)
end if
! psi grid at mid points, size n+1, for use in pec_tab
rtabpsi1(it) = rhop1**2
! ψ_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
call equil%flux_average(rhop1, area=areai1, volume=voli1)
dvol(it) = abs(voli1 - voli0)
darea(it) = abs(areai1 - areai0)
voli0 = voli1
areai0 = areai1
! 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
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
! 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 pec_init
end subroutine 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 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
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
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 :: 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
integer :: i, ray, bin
integer :: cross(maxval(last_step)+1), ncross
real(wp_) :: I_entry, I_exit, P_entry, P_exit
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
! 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
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
! 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))
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)
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
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)
! 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
dppa = ppa2 - ppa1
didst = cci2 - cci1
wdpdv(ind) = wdpdv(ind) + dppa
wajphiv(ind) = wajphiv(ind) + didst
ppa1 = ppa2
cci1 = cci2
end if
end do
! 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
end subroutine pec_tab
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 postproc_profiles(equil, rhot_tab, dpdv, ajphiv, &
rhotpav,drhotpav,rhotjava,drhotjava,dpdvp, ajphip, &
rhotp, drhotp, rhotjfi, drhotjfi, dpdvmx,ajmxfi, ratjamx, ratjbmx)
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:
! Radial average values over power and current density profile
use const_and_precisions, only : pi, comp_eps
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) :: 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
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*dvol)
P_tot = sum(dPdV*self%dV)
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)
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
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
rho_avg_P = 0
drho_avg_P = 0
end if
! ⟨ρ_t⟩, Δρ_t using |J_φ|⋅dA as a PDF
I_tot = sum(abs(ajphiv)*darea)
I_tot = sum(abs(J_phi)*self%dA)
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)
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
rhotjava = 0
drhotjava = 0
rho_avg_J = 0
drho_avg_J = 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)
! 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
dpdvp = 0
rhotp = 0
dpdvmx = 0
drhotp = 0
dPdV_peak = 0
rho_max_P = 0
dPdV_max = 0
drho_max_P = 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)
! 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
ajphip = 0
rhotjfi = 0
ajmxfi = 0
drhotjfi = 0
ratjamx = 0
ratjbmx = 0
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 postproc_profiles
end subroutine statistics
subroutine profwidth(xx,yy,xpk,ypk,dxxe)
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) :: xx(:), yy(:)
real(wp_), intent(out) :: xpk, ypk, dxxe
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
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)
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
@ -356,49 +299,37 @@ contains
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)
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(yy(ipk:size(yy)), yye, ie)
if(ie > 0 .and. ie < size(yy)) then
call locate(y(ipk:size(y)), yye, ie)
if(ie > 0 .and. ie < size(y)) then
ie = ie + (ipk - 1)
rte2 = linear_interp(yy(ie), xx(ie), yy(ie+1), xx(ie+1), yye)
rte2 = linear_interp(y(ie), x(ie), y(ie+1), x(ie+1), yye)
else
rte2 = 0
end if
else
ipk=2
xpk=xx(2)
ypk=yy(2)
x_max=x(2)
y_max=y(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)
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
dxxe = rte2 - rte1
if(ymx /= 0 .and. ymn /= 0) dxxe = -dxxe
fwhm = rte2 - rte1
if(ymx /= 0 .and. ymn /= 0) fwhm = -fwhm
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
end module gray_ec_profiles

2
src/utils.f90
View File

@ -45,7 +45,7 @@ contains
pure subroutine locate_unord(array, value, locs, n, nlocs)
! Given an `array` of size `n` and a `value`, finds at most
! `n` locations `locs` such that `value` is between
! `array(locs(i))` and `array(locs(i+i))`, in whichever order.
! `array(locs(i))` and `array(locs(i)+1)`, in whichever order.
! subroutine arguments
real(wp_), intent(in) :: array(:)