From 6f66317541aa94f08dee30411c8a4ef40e3cc8f2 Mon Sep 17 00:00:00 2001 From: Michele Guerini Rocco Date: Wed, 30 Oct 2024 17:43:08 +0100 Subject: [PATCH] replace pec module with an object MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit 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. --- src/gray_core.f90 | 28 +-- src/gray_jetto1beam.f90 | 3 +- src/gray_params.f90 | 9 +- src/gray_project.f90 | 310 +++++++++++++++++++++++ src/gray_tables.f90 | 2 +- src/main.f90 | 14 +- src/pec.f90 | 527 +++++++++++++++++----------------------- src/utils.f90 | 2 +- 8 files changed, 565 insertions(+), 330 deletions(-) create mode 100644 src/gray_project.f90 diff --git a/src/gray_core.f90 b/src/gray_core.f90 index e172a71..13a89ce 100644 --- a/src/gray_core.f90 +++ b/src/gray_core.f90 @@ -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 diff --git a/src/gray_jetto1beam.f90 b/src/gray_jetto1beam.f90 index 35b0c6e..83e06b1 100644 --- a/src/gray_jetto1beam.f90 +++ b/src/gray_jetto1beam.f90 @@ -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 diff --git a/src/gray_params.f90 b/src/gray_params.f90 index 574ab4a..d76dd94 100644 --- a/src/gray_params.f90 +++ b/src/gray_params.f90 @@ -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 diff --git a/src/gray_project.f90 b/src/gray_project.f90 new file mode 100644 index 0000000..f05ce1b --- /dev/null +++ b/src/gray_project.f90 @@ -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 diff --git a/src/gray_tables.f90 b/src/gray_tables.f90 index b7b84a9..34488b7 100644 --- a/src/gray_tables.f90 +++ b/src/gray_tables.f90 @@ -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 diff --git a/src/main.f90 b/src/main.f90 index ef34ad6..de992b7 100644 --- a/src/main.f90 +++ b/src/main.f90 @@ -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 diff --git a/src/pec.f90 b/src/pec.f90 index b191826..24fbcc3 100644 --- a/src/pec.f90 +++ b/src/pec.f90 @@ -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 diff --git a/src/utils.f90 b/src/utils.f90 index 88a7779..bd1d05c 100644 --- a/src/utils.f90 +++ b/src/utils.f90 @@ -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(:)