replace coreprofiles module with an object
This change replaces the `coreprofiles` module with a new `gray_plasma`
module providing the same functionality without using global variables.
- `read_profiles`, `read_profiles_an` are replaced by a single `load_plasma`
routines that handles both profiles kind (numerical, analytical).
- `scale_profiles` is merged into `load_plasma`, which besides reading
the profiles from file peforms the rescaling and interpolation based
on the `gray_parameters` settings.
- `set_profiles_spline`, `set_profiles_an`, `unset_profiles_spline`
are completely removed as the module no longer has any internal state.
- `density`, `ftemp`, `fzeff` are replaced by the `abstract_plasma`
type which provides the `dens`, `temp` and `zeff` methods for
either `numeric_plasma` or `analytic_plasma` subtypes.
This commit is contained in:
Michele Guerini Rocco
parent
a4ab741341
commit
2c441668bb
GPG Key ID:
BFBAF4C975F76450
@ -1,519 +0,0 @@
|
||||
! This module handles the loading, interpolation and evaluation of the
|
||||
! plasma profiles (density, temperature, effective charge)
|
||||
!
|
||||
! Two kinds of profiles are supported: analytical (suffix `_an` in the
|
||||
! subroutine names) or numerical (suffix `_spline`). For the latter, the
|
||||
! the data is interpolated using splines.
|
||||
module coreprofiles
|
||||
use const_and_precisions, only : wp_, zero, one
|
||||
use splines, only : spline_simple, spline_1d
|
||||
|
||||
implicit none
|
||||
|
||||
! Parameters of the C² polynomial tail of the density spline
|
||||
type density_tail
|
||||
real(wp_) :: start ! ψ₀, start of the tail
|
||||
real(wp_) :: end ! ψ₁, end of the end
|
||||
real(wp_) :: value ! s(ψ₀), value at the start
|
||||
real(wp_) :: deriv1 ! s'(ψ₀), first derivative at the start
|
||||
real(wp_) :: deriv2 ! s"(ψ₀), second derivative at the start
|
||||
end type
|
||||
|
||||
! Parameters of the analytical profiles model
|
||||
type analytic_model
|
||||
real(wp_) :: dens0 ! Density scaling factor
|
||||
real(wp_) :: n1, n2 ! Density exponents
|
||||
real(wp_) :: te0, te1 ! Temperature at ψ=0, ψ=1
|
||||
real(wp_) :: t1, t2 ! Temperature exponents
|
||||
real(wp_) :: zeff ! Effective charge
|
||||
end type
|
||||
|
||||
! Global variables storing the state of the module
|
||||
type(spline_1d), save :: dens_spline
|
||||
type(spline_simple), save :: temp_spline
|
||||
type(spline_simple), save :: zeff_spline
|
||||
type(density_tail), save :: tail
|
||||
type(analytic_model), save :: model
|
||||
|
||||
integer, save :: iprof
|
||||
|
||||
private
|
||||
public read_profiles, read_profiles_an ! Reading data files
|
||||
public scale_profiles ! Applying rescaling
|
||||
public density, temp, fzeff ! Accessing interpolated data
|
||||
public set_profiles_spline, set_profiles_an ! Initialising internal state
|
||||
public unset_profiles_spline ! Deinitialising internal state
|
||||
|
||||
contains
|
||||
|
||||
subroutine density(psin, dens, ddens)
|
||||
! Computes the density its first derivative as a function of
|
||||
! normalised poloidal flux.
|
||||
!
|
||||
! Note: density has units of 10¹⁹ m⁻³.
|
||||
use logger, only : log_error
|
||||
|
||||
! subroutine arguments
|
||||
real(wp_), intent(in) :: psin ! normalised poloidal flux
|
||||
real(wp_), intent(out) :: dens, ddens ! density and first derivative
|
||||
|
||||
! local variables
|
||||
character(256) :: msg ! for log messages formatting
|
||||
|
||||
! Initialise both to zero
|
||||
dens = 0
|
||||
ddens = 0
|
||||
|
||||
! Bail out when ψ is not available
|
||||
if (psin < 0) return
|
||||
|
||||
if (iprof == 0) then
|
||||
! Use the analytical model
|
||||
!
|
||||
! n(ψ) = dens0⋅(1 - ψ^aln1)^aln2
|
||||
!
|
||||
if (psin > 1) return
|
||||
dens = model%dens0 * (1 - psin**model%n1)**model%n2
|
||||
ddens = -model%dens0 * model%n1*model%n2 * psin**(model%n1 - 1) &
|
||||
* (1 - psin**model%n1)**(model%n2 - 1)
|
||||
else
|
||||
! Use the numerical data
|
||||
|
||||
! Past the tail end both density and its
|
||||
! derivatives are identically zero
|
||||
if (psin >= tail%end) return
|
||||
|
||||
if (psin < tail%start) then
|
||||
! Use the interpolating spline when in range
|
||||
|
||||
! Evaluate the spline
|
||||
dens = dens_spline%eval(psin)
|
||||
ddens = dens_spline%deriv(psin)
|
||||
|
||||
! Evaluate the spline 1st derivative
|
||||
if (abs(dens) < 1.0e-10_wp_) dens = zero
|
||||
else
|
||||
! Use a C² polynomial extension outside (ψ > ψ₀)
|
||||
|
||||
! The tail consists of the product p(ψ)⋅h(t), where:
|
||||
!
|
||||
! - p(ψ) is the 2nd order Taylor polynomial of the spline,
|
||||
! centered at ψ₀. See set_profiles_spline for details.
|
||||
!
|
||||
! - h(t) is a "smoothing" polynomial in the variable
|
||||
! t = (ψ - ψ₀)/(ψ₁ - ψ₀), defined as:
|
||||
!
|
||||
! h(t) = (1 - t)³⋅(1 + 3t + 6t²)
|
||||
!
|
||||
! with the following properties:
|
||||
!
|
||||
! h(0) = 1 h'(0)=0 h"(0)=0
|
||||
! h(1) = 0 h'(1)=0 h"(1)=0
|
||||
block
|
||||
real(wp_) :: dpsi, t, p, dp, h, dh
|
||||
dpsi = psin - tail%start ! Δψ = (ψ - ψ₀)
|
||||
|
||||
! Taylor polynomial p(ψ) and its derivative
|
||||
p = tail%value + dpsi*tail%deriv1 + dpsi**2*tail%deriv2/2
|
||||
dp = tail%deriv1 + dpsi*tail%deriv2
|
||||
|
||||
! Smoothing polynomial h(t) and its derivative
|
||||
t = dpsi/(tail%end - tail%start)
|
||||
h = (1 - t)**3 * (1 + 3*t + 6*t**2)
|
||||
dh = -30*(1 - t)**2 * t**2 / (tail%end - tail%start)
|
||||
|
||||
dens = p*h
|
||||
ddens = dp*h + p*dh
|
||||
end block
|
||||
end if
|
||||
|
||||
if (dens < 0) then
|
||||
write (msg, '("negative density:", 2(x,a,"=",g0.3))') &
|
||||
'ne', dens, 'ψ', psin
|
||||
call log_error(msg, mod='coreprofiles', proc='density')
|
||||
end if
|
||||
|
||||
end if
|
||||
end subroutine density
|
||||
|
||||
|
||||
function temp(psin)
|
||||
! Computes the temperature as a function of the
|
||||
! normalised poloidal flux.
|
||||
!
|
||||
! Note: temperature has units of keV.
|
||||
|
||||
! subroutine arguments
|
||||
real(wp_), intent(in) :: psin
|
||||
real(wp_) :: temp
|
||||
|
||||
! local variables
|
||||
real(wp_) :: proft
|
||||
|
||||
temp = 0
|
||||
if (psin >= 1 .or. psin < 0) return
|
||||
if (iprof == 0) then
|
||||
! Use the analytical model
|
||||
!
|
||||
! T(ψ) = (te0 - te1)⋅(1 - ψ^t1)^t2 + te1
|
||||
!
|
||||
proft = (1 - psin**model%t1)**model%t2
|
||||
temp = (model%te0 - model%te1)*proft + model%te1
|
||||
else
|
||||
! Use the interpolated numerical data
|
||||
temp = temp_spline%eval(psin)
|
||||
endif
|
||||
end function temp
|
||||
|
||||
|
||||
function fzeff(psin)
|
||||
! Computes the effective charge Zeff as a
|
||||
! function of the normalised poloidal flux.
|
||||
|
||||
! subroutine arguments
|
||||
real(wp_), intent(in) :: psin
|
||||
real(wp_) :: fzeff
|
||||
|
||||
fzeff = one
|
||||
if (psin >= 1 .or. psin < 0) return
|
||||
if (iprof == 0) then
|
||||
! Use the analytical model (just a constant)
|
||||
fzeff = model%zeff
|
||||
else
|
||||
! Use the interpolated numerical data
|
||||
fzeff = zeff_spline%eval(psin)
|
||||
endif
|
||||
end function fzeff
|
||||
|
||||
|
||||
subroutine read_profiles(filenm, data, err, unit)
|
||||
! Reads the radial plasma profiles from `file` and store them
|
||||
! into `data`. If given, the file is opened in the `unit` number.
|
||||
! Format notes:
|
||||
! 1. The file is formatted as a table with the following columns:
|
||||
! radial coordinate, temperature, density, effective charge.
|
||||
! 2. The first line is a header specifying the number of rows.
|
||||
use utils, only : get_free_unit
|
||||
use gray_params, only : profiles_data, PROF_NUMERIC
|
||||
use logger, only : log_error
|
||||
|
||||
! subroutine arguments
|
||||
character(len=*), intent(in) :: filenm
|
||||
type(profiles_data), intent(out) :: data
|
||||
integer, intent(out) :: err
|
||||
integer, optional, intent(in) :: unit
|
||||
|
||||
! local variables
|
||||
integer :: u, i, nrows
|
||||
|
||||
iprof = PROF_NUMERIC
|
||||
|
||||
! Free the arrays when already allocated
|
||||
if (allocated(data%psrad)) deallocate(data%psrad)
|
||||
if (allocated(data%terad)) deallocate(data%terad)
|
||||
if (allocated(data%derad)) deallocate(data%derad)
|
||||
if (allocated(data%zfc)) deallocate(data%zfc)
|
||||
|
||||
u = get_free_unit(unit)
|
||||
|
||||
! Read number of rows and allocate the arrays
|
||||
open(file=filenm, status='old', action='read', unit=u, iostat=err)
|
||||
if (err /= 0) then
|
||||
call log_error('opening profiles file ('//trim(filenm)//') failed!', &
|
||||
mod='coreprofiles', proc="read_profiles")
|
||||
err = 1
|
||||
return
|
||||
end if
|
||||
|
||||
read(u, *) nrows
|
||||
allocate(data%psrad(nrows), data%terad(nrows), &
|
||||
data%derad(nrows), data%zfc(nrows))
|
||||
|
||||
! Read the table rows
|
||||
do i=1,nrows
|
||||
read(u, *) data%psrad(i), data%terad(i), data%derad(i), data%zfc(i)
|
||||
end do
|
||||
close(u)
|
||||
|
||||
! ??
|
||||
data%psrad(1) = max(data%psrad(1), zero)
|
||||
|
||||
end subroutine read_profiles
|
||||
|
||||
|
||||
subroutine read_profiles_an(filenm, data, err, unit)
|
||||
! Reads the plasma profiles `data` in the analytical format
|
||||
! from params%filenm.
|
||||
! If given, the file is opened in the `unit` number.
|
||||
!
|
||||
! The file should be formatted as follows:
|
||||
!
|
||||
! 1 dens0 n1 n2
|
||||
! 2 te0 te1 t1 t2
|
||||
! 3 zeff
|
||||
!
|
||||
! See `density`, `temp`, `fzeff` subroutines for the meaning
|
||||
! of the parameters (i.e. the formulae for n,T,Zeff).
|
||||
use gray_params, only : profiles_data, PROF_ANALYTIC
|
||||
use utils, only : get_free_unit
|
||||
use logger, only : log_error
|
||||
|
||||
! subroutine arguments
|
||||
character(len=*), intent(in) :: filenm
|
||||
type(profiles_data), intent(out) :: data
|
||||
integer, intent(out) :: err
|
||||
integer, optional, intent(in) :: unit
|
||||
|
||||
! local variables
|
||||
integer :: u
|
||||
|
||||
u = get_free_unit(unit)
|
||||
|
||||
iprof = PROF_ANALYTIC
|
||||
|
||||
if (allocated(data%terad)) deallocate(data%terad)
|
||||
if (allocated(data%derad)) deallocate(data%derad)
|
||||
if (allocated(data%zfc)) deallocate(data%zfc)
|
||||
allocate(data%terad(4), data%derad(3), data%zfc(1))
|
||||
|
||||
open(file=filenm, status='old', action='read', unit=u, iostat=err)
|
||||
if (err /= 0) then
|
||||
call log_error('opening profiles file ('//trim(filenm)//') failed!', &
|
||||
mod='coreprofiles', proc='read_profiles_an')
|
||||
err = 1
|
||||
return
|
||||
end if
|
||||
|
||||
read (u,*) data%derad(1:3) ! dens0, n1, n2
|
||||
read (u,*) data%terad(1:4) ! te0, te1, t1, t2
|
||||
read (u,*) data%zfc(1) ! zeff
|
||||
close(u)
|
||||
end subroutine read_profiles_an
|
||||
|
||||
|
||||
subroutine scale_profiles(params, factb, data)
|
||||
! Rescales the temperature and density profiles by `params%factte`
|
||||
! and `params%factne`, respectively, according to the model
|
||||
! specified by `params%iscal`.
|
||||
! See the GRAY user manual for the explanation.
|
||||
use gray_params, only : profiles_parameters, profiles_data
|
||||
|
||||
! subroutine arguments
|
||||
type(profiles_parameters), intent(in) :: params
|
||||
real(wp_), intent(in) :: factb
|
||||
type(profiles_data), intent(inout) :: data
|
||||
|
||||
! local variables
|
||||
real(wp_) :: aat, aan, ffact
|
||||
integer :: last_te, last_ne
|
||||
|
||||
if (params%iscal==0) then
|
||||
aat = 2.0_wp_/3.0_wp_
|
||||
aan = 4.0_wp_/3.0_wp_
|
||||
else
|
||||
aat = one
|
||||
aan = one
|
||||
end if
|
||||
|
||||
if (params%iscal==2) then
|
||||
ffact = one
|
||||
else
|
||||
ffact = factb
|
||||
end if
|
||||
|
||||
if (params%iprof==0) then
|
||||
last_te = 2
|
||||
last_ne = 1
|
||||
else
|
||||
last_te = size(data%terad)
|
||||
last_ne = size(data%derad)
|
||||
end if
|
||||
|
||||
data%terad(1:last_te) = data%terad(1:last_te) * ffact**aat * params%factte
|
||||
data%derad(1:last_ne) = data%derad(1:last_ne) * ffact**aan * params%factne
|
||||
end subroutine scale_profiles
|
||||
|
||||
|
||||
subroutine set_profiles_spline(params, data, err, launch_pos)
|
||||
! Computes splines for the plasma profiles data and stores them
|
||||
! in their respective global variables, see the top of this file.
|
||||
!
|
||||
! When `launch_pos` (cartesian launch coordinates in cm) is present,
|
||||
! the subroutine will also check that the wave launcher is strictly
|
||||
! outside the reconstructed plasma density boundary.
|
||||
!
|
||||
! `err` is 1 if I/O errors occured, 2 if other initialisation failed.
|
||||
use gray_params, only : profiles_parameters, profiles_data
|
||||
use logger, only : log_debug, log_info, log_warning, log_error
|
||||
|
||||
! subroutine arguments
|
||||
type(profiles_parameters), intent(inout) :: params
|
||||
type(profiles_data), intent(inout) :: data
|
||||
integer, intent(out) :: err
|
||||
real(wp_), optional, intent(in) :: launch_pos(3)
|
||||
|
||||
! local variables
|
||||
integer :: n
|
||||
|
||||
! for log messages formatting
|
||||
character(256) :: msg
|
||||
|
||||
n = size(data%psrad)
|
||||
|
||||
! Spline interpolation of temperature and effective charge
|
||||
call temp_spline%init(data%psrad, data%terad, n)
|
||||
call zeff_spline%init(data%psrad, data%zfc, n)
|
||||
|
||||
! Spline interpolation of density (smooth spline)
|
||||
call dens_spline%init(data%psrad, data%derad, n, range=[zero, data%psrad(n)], &
|
||||
tension=params%sspld, err=err)
|
||||
|
||||
! Tension is very small, re-fit with an interpolating spline
|
||||
if (err == -1) then
|
||||
call log_warning('re-fitting with density curve with zero tension', &
|
||||
mod='coreprofiles', proc='density')
|
||||
call dens_spline%init(data%psrad, data%derad, n, &
|
||||
range=[zero, data%psrad(n)], tension=zero)
|
||||
else if (err > 0) then
|
||||
write (msg, '(a, g0)') 'density fit failed with error ', err
|
||||
call log_error(msg, mod='coreprofiles', proc='density')
|
||||
err = 2
|
||||
return
|
||||
else
|
||||
err = 0
|
||||
end if
|
||||
|
||||
! Computation of the polynomial tail parameters
|
||||
!
|
||||
! Note: The density is the only quantity that needs to be evaluated
|
||||
! at the edge. The spline thus has to be extended to transition
|
||||
! smoothly from the last profile point to 0 outside the plasma.
|
||||
block
|
||||
real(wp_) :: s0, s1, s2 ! spline, 1st, 2nd derivative
|
||||
real(wp_) :: delta4 ! discriminant Δ/4 of q(x)
|
||||
real(wp_) :: x0, x1 ! vertex of q(x), solution
|
||||
|
||||
! Compute the coefficients of a 2nd order Taylor polinomial to
|
||||
! extend the spline beyond the last point:
|
||||
!
|
||||
! p(ψ) = s(ψ₀) + (ψ - ψ₀)s'(ψ₀) + ½(ψ - ψ₀)²s"(ψ₀)
|
||||
!
|
||||
! where s(ψ) is the spline and ψ₀ the last point.
|
||||
!
|
||||
s0 = dens_spline%eval(data%psrad(n))
|
||||
s1 = dens_spline%deriv(data%psrad(n), order=1)
|
||||
s2 = dens_spline%deriv(data%psrad(n), order=2)
|
||||
|
||||
! Determine where to end the tail (to ensure the density remains
|
||||
! positive) from the zeros of the Taylor polynomial p(ψ)
|
||||
!
|
||||
! Define x=(ψ - ψ₀), then p(ψ)=0 is rewritten as
|
||||
!
|
||||
! q(x) = x² + 2s'/s" x + 2s/s" = 0
|
||||
!
|
||||
! The discriminant is Δ/4 = (s'/s")² - 2(s/s") and
|
||||
! the solutions are x = x₀ ± √(Δ/4), with x₀ = -s'/s".
|
||||
!
|
||||
x0 = -s1 / s2 ! vertex of parabola y=q(x)
|
||||
delta4 = (s1 / s2)**2 - 2*s0/s2 ! Δ/4 of q(x)
|
||||
|
||||
if (delta4 > 0) then
|
||||
! Pick the smallest positive zero (implying >ψ₀)
|
||||
x1 = x0 + sign(sqrt(delta4), sqrt(delta4) - x0)
|
||||
else
|
||||
! There are no zeros, use the parabola vertex
|
||||
x1 = x0
|
||||
call log_debug('spline extension has no zeros', &
|
||||
mod='coreprofiles', proc='set_profiles_spline')
|
||||
end if
|
||||
|
||||
! Store the tail parameters
|
||||
tail%start = data%psrad(n)
|
||||
tail%end = tail%start + x1
|
||||
tail%value = s0
|
||||
tail%deriv1 = s1
|
||||
tail%deriv2 = s2
|
||||
end block
|
||||
|
||||
! Make sure the wave launcher does not fall inside the tail
|
||||
! Note: if it does, the initial wave conditions become
|
||||
! invalid as they are given assuming a vacuum (N=1)
|
||||
if (present(launch_pos)) then
|
||||
block
|
||||
use equilibrium, only : pol_flux
|
||||
use const_and_precisions, only : cm
|
||||
real(wp_) :: R, Z, psi
|
||||
|
||||
! Convert from cartesian to cylindrical coordinates
|
||||
R = hypot(launch_pos(1), launch_pos(2)) * cm ! R = √(x²+y²)
|
||||
z = launch_pos(3) * cm
|
||||
|
||||
! Get the poloidal flux at the launcher
|
||||
! Note: this returns -1 when the data is not available
|
||||
call pol_flux(R, z, psi)
|
||||
|
||||
if (psi > tail%start .and. psi < tail%end) then
|
||||
! Fall back to the midpoint of ψ₀ and the launcher ψ
|
||||
call log_warning('downscaled tail to not reach the wave launcher', &
|
||||
mod='coreprofiles', proc='set_profiles_spline')
|
||||
write (msg, '(a, g0.3, a, g0.3, " → ", g0.3)') &
|
||||
'launcher: ψ=', psi, ' boundary: ψ=', tail%end, (tail%start + psi)/2
|
||||
call log_debug(msg, mod='coreprofiles', proc='set_profiles_spline')
|
||||
tail%end = (tail%start + psi)/2
|
||||
end if
|
||||
|
||||
if (psi > 0 .and. psi < tail%start) then
|
||||
! This must be a user error, stop here
|
||||
write (msg, '(a, a, g0.3, a, g0.3)') &
|
||||
'wave launcher is inside the plasma! ', &
|
||||
'launcher: ψ=', psi, ' boundary: ψ=', tail%end
|
||||
call log_error(msg, mod='coreprofiles', proc='set_profiles_spline')
|
||||
err = 2
|
||||
return
|
||||
end if
|
||||
end block
|
||||
end if
|
||||
|
||||
! Set the density boundary ψ
|
||||
! Note: this is used to detect entrance in the plasma
|
||||
params%psnbnd = tail%end
|
||||
|
||||
write (msg, '(a,g0.4)') 'density boundary: ψ=', tail%end
|
||||
call log_info(msg, mod='coreprofiles', proc='set_profiles_spline')
|
||||
end subroutine set_profiles_spline
|
||||
|
||||
|
||||
subroutine unset_profiles_spline
|
||||
! Unsets the splines global variables, see the top of this file.
|
||||
|
||||
call dens_spline%deinit
|
||||
call temp_spline%deinit
|
||||
call zeff_spline%deinit
|
||||
end subroutine unset_profiles_spline
|
||||
|
||||
|
||||
subroutine set_profiles_an(params, data)
|
||||
! Stores the analytical profiles data in their respective
|
||||
! global variables, see the top of this file.
|
||||
use gray_params, only : profiles_parameters, profiles_data
|
||||
|
||||
! subroutine arguments
|
||||
type(profiles_parameters), intent(inout) :: params
|
||||
type(profiles_data), intent(in) :: data
|
||||
|
||||
model%te0 = data%terad(1)
|
||||
model%te1 = data%terad(2)
|
||||
model%t1 = data%terad(3)
|
||||
model%t2 = data%terad(4)
|
||||
model%dens0 = data%derad(1)
|
||||
model%n1 = data%derad(2)
|
||||
model%n2 = data%derad(3)
|
||||
model%zeff = data%zfc(1)
|
||||
|
||||
! Define the plasma boundary to be exactly ψ=1
|
||||
! Note: this is used to detect entrance in the plasma
|
||||
params%psnbnd = one
|
||||
|
||||
end subroutine set_profiles_an
|
||||
|
||||
end module coreprofiles
|
||||
@ -7,12 +7,12 @@ module gray_core
|
||||
|
||||
contains
|
||||
|
||||
subroutine gray_main(params, data, results, error, rhout)
|
||||
subroutine gray_main(params, data, plasma, results, error, rhout)
|
||||
use const_and_precisions, only : zero, one, comp_tiny
|
||||
use coreprofiles, only : temp, fzeff
|
||||
use polarization, only : ellipse_to_field
|
||||
use types, only : table, wrap
|
||||
use gray_params, only : gray_parameters, gray_data, gray_results, EQ_VACUUM
|
||||
use gray_plasma, only : abstract_plasma
|
||||
use gray_tables, only : init_tables, store_ec_profiles, store_ray_data, &
|
||||
store_beam_shape, find_flux_surfaces, &
|
||||
kinetic_profiles, ec_resonance, inputs_maps
|
||||
@ -30,6 +30,7 @@ contains
|
||||
! subroutine arguments
|
||||
type(gray_parameters), intent(inout) :: params
|
||||
type(gray_data), intent(in) :: data
|
||||
class(abstract_plasma), intent(in) :: plasma
|
||||
type(gray_results), intent(out) :: results
|
||||
|
||||
! Predefined grid for the output profiles (optional)
|
||||
@ -132,9 +133,9 @@ contains
|
||||
! ========= set environment END =========
|
||||
|
||||
! Pre-determinted tables
|
||||
results%tables%kinetic_profiles = kinetic_profiles(params)
|
||||
results%tables%kinetic_profiles = kinetic_profiles(params, plasma)
|
||||
results%tables%ec_resonance = ec_resonance(params, bres)
|
||||
results%tables%inputs_maps = inputs_maps(params, bres, xgcn)
|
||||
results%tables%inputs_maps = inputs_maps(params, plasma, bres, xgcn)
|
||||
|
||||
! print ψ surface for q=3/2 and q=2/1
|
||||
call find_flux_surfaces( &
|
||||
@ -257,8 +258,8 @@ contains
|
||||
do jk=1,params%raytracing%nray
|
||||
if(iwait(jk)) cycle ! jk ray is waiting for next pass
|
||||
stv(jk) = stv(jk) + params%raytracing%dst ! current ray step
|
||||
call rkstep(params, sox, bres, xgcn, yw(:,jk), ypw(:,jk), &
|
||||
gri(:,jk), ggri(:,:,jk), igrad_b)
|
||||
call rkstep(params, plasma, sox, bres, xgcn, yw(:,jk), &
|
||||
ypw(:,jk), gri(:,jk), ggri(:,:,jk), igrad_b)
|
||||
end do
|
||||
! update position and grad
|
||||
if(igrad_b == 1) call gradi_upd(params%raytracing, yw, ak0, xc, du1, gri, ggri)
|
||||
@ -273,9 +274,10 @@ contains
|
||||
! compute derivatives with updated gradient and local plasma values
|
||||
xv = yw(1:3,jk)
|
||||
anv = yw(4:6,jk)
|
||||
call ywppla_upd(params, xv, anv, gri(:,jk), ggri(:,:,jk), sox, bres, &
|
||||
xgcn,ypw(:,jk), psinv, dens, btot, bv, xg, yg, derxg, &
|
||||
anpl, anpr, ddr, ddi, dersdst, derdnm, ierrn, igrad_b)
|
||||
call ywppla_upd(params, plasma, xv, anv, gri(:,jk), ggri(:,:,jk), &
|
||||
sox, bres, xgcn,ypw(:,jk), psinv, dens, btot, bv, &
|
||||
xg, yg, derxg, anpl, anpr, ddr, ddi, dersdst, &
|
||||
derdnm, ierrn, igrad_b)
|
||||
! update global error code and print message
|
||||
if(ierrn/=0) then
|
||||
error = ior(error,ierrn)
|
||||
@ -372,9 +374,9 @@ contains
|
||||
yw(:,jk) = [xv, anv] ! * updated coordinates (reflected)
|
||||
igrad_b = 0 ! * switch to ray-tracing
|
||||
|
||||
call ywppla_upd(params, xv, anv, gri(:,jk), ggri(:,:,jk), sox, &
|
||||
bres, xgcn, ypw(:,jk), psinv, dens, btot, bv, &
|
||||
xg, yg, derxg, anpl, anpr, ddr, ddi, dersdst, &
|
||||
call ywppla_upd(params, plasma, xv, anv, gri(:,jk), ggri(:,:,jk), &
|
||||
sox, bres, xgcn, ypw(:,jk), psinv, dens, btot, &
|
||||
bv, xg, yg, derxg, anpl, anpr, ddr, ddi, dersdst, &
|
||||
derdnm, ierrn, igrad_b) ! * update derivatives after reflection
|
||||
if(ierrn/=0) then ! * update global error code and print message
|
||||
error = ior(error,ierrn)
|
||||
@ -429,12 +431,13 @@ contains
|
||||
|
||||
if(ierrn==0 .and. params%ecrh_cd%iwarm>0 .and. ins_pl .and. &
|
||||
(tau1(jk)+tau0(jk)+lgcpl1(jk))<=taucr .and. .not.iwait(jk)) then ! H&CD computation check
|
||||
tekev=temp(psinv)
|
||||
tekev = plasma%temp(psinv)
|
||||
block
|
||||
complex(wp_) :: Npr_warm
|
||||
call alpha_effj(params%ecrh_cd, psinv, xg, yg, dens, tekev, &
|
||||
ak0, bres, derdnm, anpl, anpr, sox, Npr_warm, &
|
||||
alpha, didp, nharm, nhf, iokhawa, ierrhcd)
|
||||
call alpha_effj(params%ecrh_cd, plasma, psinv, xg, yg, dens, &
|
||||
tekev, ak0, bres, derdnm, anpl, anpr, sox, &
|
||||
Npr_warm, alpha, didp, nharm, nhf, iokhawa, &
|
||||
ierrhcd)
|
||||
anprre = Npr_warm%re
|
||||
anprim = Npr_warm%im
|
||||
if(ierrhcd /= 0) then
|
||||
@ -995,12 +998,15 @@ contains
|
||||
|
||||
|
||||
|
||||
subroutine rkstep(params, sox, bres, xgcn, y, yp, dgr, ddgr, igrad)
|
||||
subroutine rkstep(params, plasma, sox, bres, xgcn, y, yp, dgr, ddgr, igrad)
|
||||
! Runge-Kutta integrator
|
||||
use gray_params, only : gray_parameters
|
||||
use gray_plasma, only : abstract_plasma
|
||||
|
||||
! subroutine arguments
|
||||
type(gray_parameters), intent(in) :: params
|
||||
class(abstract_plasma), intent(in) :: plasma
|
||||
|
||||
real(wp_), intent(in) :: bres, xgcn
|
||||
real(wp_), intent(inout) :: y(6)
|
||||
real(wp_), intent(in) :: yp(6)
|
||||
@ -1024,18 +1030,20 @@ contains
|
||||
function f(y)
|
||||
real(wp_), intent(in) :: y(6)
|
||||
real(wp_) :: f(6)
|
||||
call rhs(params, sox, bres, xgcn, y, dgr, ddgr, f, igrad)
|
||||
call rhs(params, plasma, sox, bres, xgcn, y, dgr, ddgr, f, igrad)
|
||||
end function
|
||||
end subroutine rkstep
|
||||
|
||||
|
||||
subroutine rhs(params, sox, bres, xgcn, y, dgr, ddgr, dery, igrad)
|
||||
subroutine rhs(params, plasma, sox, bres, xgcn, y, dgr, ddgr, dery, igrad)
|
||||
! Compute right-hand side terms of the ray equations (dery)
|
||||
! used in R-K integrator
|
||||
use gray_params, only : gray_parameters
|
||||
use gray_plasma, only : abstract_plasma
|
||||
|
||||
! subroutine arguments
|
||||
type(gray_parameters), intent(in) :: params
|
||||
class(abstract_plasma), intent(in) :: plasma
|
||||
real(wp_), intent(in) :: y(6)
|
||||
real(wp_), intent(in) :: bres, xgcn
|
||||
real(wp_), intent(in) :: dgr(3)
|
||||
@ -1049,7 +1057,7 @@ contains
|
||||
real(wp_), dimension(3,3) :: derbv
|
||||
|
||||
xv = y(1:3)
|
||||
call plas_deriv(params, xv, bres, xgcn, dens, btot, &
|
||||
call plas_deriv(params, plasma, xv, bres, xgcn, dens, btot, &
|
||||
bv, derbv, xg, yg, derxg, deryg)
|
||||
anv = y(4:6)
|
||||
call disp_deriv(params, anv, sox, xg, yg, derxg, deryg, &
|
||||
@ -1057,16 +1065,18 @@ contains
|
||||
end subroutine rhs
|
||||
|
||||
|
||||
subroutine ywppla_upd(params, xv, anv, dgr, ddgr, sox, bres, xgcn, dery, &
|
||||
subroutine ywppla_upd(params, plasma, xv, anv, dgr, ddgr, sox, bres, xgcn, dery, &
|
||||
psinv, dens, btot, bv, xg, yg, derxg, anpl, anpr, &
|
||||
ddr, ddi, dersdst, derdnm, error, igrad)
|
||||
! Compute right-hand side terms of the ray equations (dery)
|
||||
! used after full R-K step and grad(S_I) update
|
||||
use gray_errors, only : raise_error, large_npl
|
||||
use gray_params, only : gray_parameters
|
||||
use gray_plasma, only : abstract_plasma
|
||||
|
||||
! subroutine rguments
|
||||
type(gray_parameters), intent(in) :: params
|
||||
class(abstract_plasma), intent(in) :: plasma
|
||||
|
||||
real(wp_), intent(in) :: xv(3), anv(3)
|
||||
real(wp_), intent(in) :: dgr(3)
|
||||
@ -1086,7 +1096,7 @@ contains
|
||||
real(wp_), dimension(3,3) :: derbv
|
||||
real(wp_), parameter :: anplth1 = 0.99_wp_, anplth2 = 1.05_wp_
|
||||
|
||||
call plas_deriv(params, xv,bres,xgcn,dens,btot,bv,derbv,xg,yg,derxg,deryg,psinv)
|
||||
call plas_deriv(params, plasma, xv,bres,xgcn,dens,btot,bv,derbv,xg,yg,derxg,deryg,psinv)
|
||||
call disp_deriv(params, anv,sox,xg,yg,derxg,deryg,bv,derbv,dgr,ddgr,igrad, &
|
||||
dery,anpl,anpr,ddr,ddi,dersdst,derdnm)
|
||||
|
||||
@ -1305,15 +1315,16 @@ contains
|
||||
end subroutine solg3
|
||||
|
||||
|
||||
subroutine plas_deriv(params, xv, bres, xgcn, dens, btot, bv, &
|
||||
subroutine plas_deriv(params, plasma, xv, bres, xgcn, dens, btot, bv, &
|
||||
derbv, xg, yg, derxg, deryg, psinv_opt)
|
||||
use const_and_precisions, only : zero, cm
|
||||
use gray_params, only : gray_parameters, EQ_VACUUM
|
||||
use gray_plasma, only : abstract_plasma
|
||||
use equilibrium, only : psia, pol_flux, pol_curr, sgnbphi
|
||||
use coreprofiles, only : density
|
||||
|
||||
! subroutine arguments
|
||||
type(gray_parameters), intent(in) :: params
|
||||
class(abstract_plasma), intent(in) :: plasma
|
||||
real(wp_), dimension(3), intent(in) :: xv
|
||||
real(wp_), intent(in) :: xgcn, bres
|
||||
real(wp_), intent(out) :: dens, btot, xg, yg
|
||||
@ -1385,7 +1396,7 @@ contains
|
||||
end if
|
||||
|
||||
! compute xg and derivative
|
||||
call density(psinv, dens, ddensdpsi)
|
||||
call plasma%density(psinv, dens, ddensdpsi)
|
||||
xg = xgcn*dens
|
||||
dxgdpsi = xgcn*ddensdpsi
|
||||
|
||||
@ -1735,14 +1746,14 @@ contains
|
||||
|
||||
|
||||
|
||||
subroutine alpha_effj(params, psinv, X, Y, density, temperature, &
|
||||
k0, Bres, derdnm, Npl, Npr_cold, sox, &
|
||||
Npr, alpha, dIdp, nhmin, nhmax, iokhawa, error)
|
||||
subroutine alpha_effj(params, plasma, psinv, X, Y, density, temperature, &
|
||||
k0, Bres, derdnm, Npl, Npr_cold, sox, Npr, &
|
||||
alpha, dIdp, nhmin, nhmax, iokhawa, error)
|
||||
! Computes the absorption coefficient and effective current
|
||||
|
||||
use const_and_precisions, only : pi, mc2=>mc2_
|
||||
use gray_params, only : ecrh_cd_parameters
|
||||
use coreprofiles, only : fzeff
|
||||
use gray_plasma, only : abstract_plasma
|
||||
use equilibrium, only : sgnbphi
|
||||
use dispersion, only : harmnumber, warmdisp
|
||||
use eccd, only : setcdcoeff, eccdeff, fjch0, fjch, fjncl
|
||||
@ -1755,7 +1766,9 @@ contains
|
||||
! Inputs
|
||||
|
||||
! ECRH & CD parameters
|
||||
type(ecrh_cd_parameters) :: params
|
||||
type(ecrh_cd_parameters), intent(in) :: params
|
||||
! plasma object
|
||||
class(abstract_plasma), intent(in) :: plasma
|
||||
! poloidal flux ψ
|
||||
real(wp_), intent(in) :: psinv
|
||||
! CMA diagram variables: X=(ω_pe/ω)², Y=ω_ce/ω
|
||||
@ -1855,7 +1868,7 @@ contains
|
||||
! Current drive computation
|
||||
if (params%ieccd <= 0) return
|
||||
|
||||
zeff = fzeff(psinv)
|
||||
zeff = plasma%zeff(psinv)
|
||||
ithn = 1
|
||||
if (nlarmor > nhmin) ithn = 2
|
||||
rhop = sqrt(psinv)
|
||||
|
||||
@ -4,7 +4,9 @@ subroutine gray_jetto1beam(ijetto, mr, mz, r, z, psin, psia, rax, zax, &
|
||||
use const_and_precisions, only: wp_
|
||||
use gray_params, only: gray_parameters, gray_data, gray_results
|
||||
use gray_core, only: gray_main
|
||||
use gray_plasma, only: numeric_plasma
|
||||
|
||||
implicit none
|
||||
|
||||
! subroutine arguments
|
||||
integer, intent(in) :: ijetto, mr, mz, nbnd, nrho, ibeam
|
||||
@ -21,6 +23,7 @@ subroutine gray_jetto1beam(ijetto, mr, mz, r, z, psin, psia, rax, zax, &
|
||||
! local variables
|
||||
type(gray_parameters) :: params
|
||||
type(gray_data) :: data
|
||||
type(numeric_plasma) :: plasma
|
||||
type(gray_results) :: res
|
||||
logical, save :: firstcall = .true.
|
||||
|
||||
@ -86,19 +89,9 @@ subroutine gray_jetto1beam(ijetto, mr, mz, r, z, psin, psia, rax, zax, &
|
||||
if (err /= 0) return
|
||||
end block init_equilibrium
|
||||
|
||||
! Set plasma kinetic profiles
|
||||
init_profiles: block
|
||||
use coreprofiles, only : set_profiles_spline
|
||||
|
||||
! Copy argument arrays
|
||||
data%profiles%derad = dne
|
||||
data%profiles%terad = te
|
||||
data%profiles%zfc = zeff
|
||||
|
||||
! Compute splines
|
||||
call set_profiles_spline(params%profiles, data%profiles, err)
|
||||
! Compute splines of kinetic profiles
|
||||
call plasma%init(params, psrad, te, dne, zeff, err)
|
||||
if (err /= 0) return
|
||||
end block init_profiles
|
||||
|
||||
! Set wave launcher parameters
|
||||
init_antenna: block
|
||||
@ -115,7 +108,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, data, res, err, rhout=sqrt(psrad))
|
||||
call gray_main(params, data, plasma, res, err, rhout=sqrt(psrad))
|
||||
|
||||
write_debug_files: block
|
||||
integer :: i, err
|
||||
@ -130,13 +123,9 @@ subroutine gray_jetto1beam(ijetto, mr, mz, r, z, psin, psia, rax, zax, &
|
||||
! Free memory
|
||||
free_memory: block
|
||||
use equilibrium, only : unset_equil_spline
|
||||
use coreprofiles, only : unset_profiles_spline
|
||||
|
||||
! Unset global variables of the `equilibrium` module
|
||||
call unset_equil_spline
|
||||
|
||||
! Unset global variables of the `coreprofiles` module
|
||||
call unset_profiles_spline
|
||||
end block free_memory
|
||||
|
||||
! Copy over the results
|
||||
|
||||
536
src/gray_plasma.f90
Normal file
536
src/gray_plasma.f90
Normal file
@ -0,0 +1,536 @@
|
||||
! This module handles the loading, interpolation and evaluation of the
|
||||
! plasma profiles (density, temperature, effective charge)
|
||||
!
|
||||
! The type `abstract_plasma` provides an abstract interface for
|
||||
! accessing the profiles data; these implementations are available:
|
||||
! - `numeric_plasma`, numerical data interpolated using splines
|
||||
! - `analytic_plasma`, simple analytic model
|
||||
!
|
||||
module gray_plasma
|
||||
use const_and_precisions, only : wp_, zero
|
||||
use splines, only : spline_simple, spline_1d
|
||||
|
||||
implicit none
|
||||
|
||||
! Generic plasma interface
|
||||
type, abstract :: abstract_plasma
|
||||
contains
|
||||
procedure(density_sub), deferred :: density
|
||||
procedure(temp_fun), deferred :: temp
|
||||
procedure(zeff_fun), deferred :: zeff
|
||||
end type
|
||||
|
||||
! Analytical plasma description
|
||||
type, extends(abstract_plasma) :: analytic_plasma
|
||||
private
|
||||
real(wp_) :: dens0 ! Density scaling factor
|
||||
real(wp_) :: n1, n2 ! Density exponents
|
||||
real(wp_) :: te0, te1 ! Temperature at ψ=0, ψ=1
|
||||
real(wp_) :: t1, t2 ! Temperature exponents
|
||||
real(wp_) :: z ! Effective charge
|
||||
contains
|
||||
procedure :: density => analytic_density
|
||||
procedure :: temp => analytic_temp
|
||||
procedure :: zeff => analytic_zeff
|
||||
end type
|
||||
|
||||
! Parameters of the C² polynomial tail of the density spline
|
||||
type density_tail
|
||||
real(wp_) :: start ! ψ₀, start of the tail
|
||||
real(wp_) :: end ! ψ₁, end of the end
|
||||
real(wp_) :: value ! s(ψ₀), value at the start
|
||||
real(wp_) :: deriv1 ! s'(ψ₀), first derivative at the start
|
||||
real(wp_) :: deriv2 ! s"(ψ₀), second derivative at the start
|
||||
end type
|
||||
|
||||
! Numerical plasma description
|
||||
type, extends(abstract_plasma) :: numeric_plasma
|
||||
private
|
||||
type(spline_1d) :: dens_spline
|
||||
type(spline_simple) :: temp_spline
|
||||
type(spline_simple) :: zeff_spline
|
||||
type(density_tail) :: tail
|
||||
contains
|
||||
procedure :: init => init_splines
|
||||
procedure :: density => numeric_density
|
||||
procedure :: temp => numeric_temp
|
||||
procedure :: zeff => numeric_zeff
|
||||
end type
|
||||
|
||||
abstract interface
|
||||
subroutine density_sub(self, psin, dens, ddens)
|
||||
! Computes the density its first derivative as a function of
|
||||
! normalised poloidal flux.
|
||||
!
|
||||
! Note: density has units of 10¹⁹ m⁻³.
|
||||
import :: abstract_plasma, wp_
|
||||
class(abstract_plasma), intent(in) :: self
|
||||
real(wp_), intent(in) :: psin ! normalised poloidal flux
|
||||
real(wp_), intent(out) :: dens, ddens ! density and first derivative
|
||||
end subroutine density_sub
|
||||
|
||||
function temp_fun(self, psin) result(temp)
|
||||
! Computes the temperature as a function of the
|
||||
! normalised poloidal flux.
|
||||
!
|
||||
! Note: temperature has units of keV.
|
||||
import :: abstract_plasma, wp_
|
||||
class(abstract_plasma), intent(in) :: self
|
||||
real(wp_), intent(in) :: psin
|
||||
real(wp_) :: temp
|
||||
end function temp_fun
|
||||
|
||||
function zeff_fun(self, psin) result(zeff)
|
||||
! Computes the effective charge Z_eff as a
|
||||
! function of the normalised poloidal flux.
|
||||
import :: abstract_plasma, wp_
|
||||
class(abstract_plasma), intent(in) :: self
|
||||
real(wp_), intent(in) :: psin
|
||||
real(wp_) :: zeff
|
||||
end function zeff_fun
|
||||
end interface
|
||||
|
||||
! Parameters of the analytical profiles model
|
||||
type analytic_model
|
||||
real(wp_) :: dens0 ! Density scaling factor
|
||||
real(wp_) :: n1, n2 ! Density exponents
|
||||
real(wp_) :: te0, te1 ! Temperature at ψ=0, ψ=1
|
||||
real(wp_) :: t1, t2 ! Temperature exponents
|
||||
real(wp_) :: zeff ! Effective charge
|
||||
end type
|
||||
|
||||
private
|
||||
public abstract_plasma ! The abstract plasma object
|
||||
public analytic_plasma, numeric_plasma ! Implementations
|
||||
public load_plasma ! To load plasma from file
|
||||
|
||||
contains
|
||||
|
||||
subroutine analytic_density(self, psin, dens, ddens)
|
||||
! subroutine arguments
|
||||
class(analytic_plasma), intent(in) :: self
|
||||
real(wp_), intent(in) :: psin ! normalised poloidal flux
|
||||
real(wp_), intent(out) :: dens, ddens ! density and first derivative
|
||||
|
||||
if (psin > 1) then
|
||||
! Outside: the density and every derivative is identically zero
|
||||
dens = 0
|
||||
ddens = 0
|
||||
else
|
||||
! Inside: the density is given by n(ψ) = dens0⋅(1 - ψ^n1)^n2
|
||||
dens = self%dens0 * (1 - psin**self%n1)**self%n2
|
||||
ddens = -self%dens0 * self%n1*self%n2 * psin**(self%n1 - 1) &
|
||||
* (1 - psin**self%n1)**(self%n2 - 1)
|
||||
end if
|
||||
end subroutine analytic_density
|
||||
|
||||
|
||||
function analytic_temp(self, psin) result(temp)
|
||||
! function arguments
|
||||
class(analytic_plasma), intent(in) :: self
|
||||
real(wp_), intent(in) :: psin
|
||||
real(wp_) :: temp
|
||||
|
||||
! local variables
|
||||
real(wp_) :: proft
|
||||
|
||||
if (psin >= 1 .or. psin < 0) then
|
||||
! Outside: temperature is undefined, return 0
|
||||
temp = 0
|
||||
else
|
||||
! Inside: the temperature is given by
|
||||
!
|
||||
! T(ψ) = (te0 - te1)⋅(1 - ψ^t1)^t2 + te1
|
||||
!
|
||||
proft = (1 - psin**self%t1)**self%t2
|
||||
temp = (self%te0 - self%te1)*proft + self%te1
|
||||
end if
|
||||
end function analytic_temp
|
||||
|
||||
|
||||
function analytic_zeff(self, psin) result(zeff)
|
||||
! function arguments
|
||||
class(analytic_plasma), intent(in) :: self
|
||||
real(wp_), intent(in) :: psin
|
||||
real(wp_) :: zeff
|
||||
|
||||
if (psin >= 1 .or. psin < 0) then
|
||||
! Outside: Zeff is undefined, return 0
|
||||
zeff = 0
|
||||
else
|
||||
! Inside: just a constant
|
||||
zeff = self%z
|
||||
end if
|
||||
end function analytic_zeff
|
||||
|
||||
|
||||
subroutine numeric_density(self, psin, dens, ddens)
|
||||
use logger, only : log_error
|
||||
|
||||
! subroutine arguments
|
||||
class(numeric_plasma), intent(in) :: self
|
||||
real(wp_), intent(in) :: psin
|
||||
real(wp_), intent(out) :: dens, ddens
|
||||
|
||||
! local variables
|
||||
character(256) :: msg ! for log messages formatting
|
||||
|
||||
! Initialise both to zero
|
||||
dens = 0
|
||||
ddens = 0
|
||||
|
||||
! Bail out when ψ is not available
|
||||
if (psin < 0) return
|
||||
|
||||
! Past the tail end both density and its
|
||||
! derivatives are identically zero
|
||||
if (psin >= self%tail%end) return
|
||||
|
||||
if (psin < self%tail%start) then
|
||||
! Use the interpolating spline when in range
|
||||
|
||||
! Evaluate the spline
|
||||
dens = self%dens_spline%eval(psin)
|
||||
ddens = self%dens_spline%deriv(psin)
|
||||
|
||||
! Evaluate the spline 1st derivative
|
||||
if (abs(dens) < 1.0e-10_wp_) dens = zero
|
||||
else
|
||||
! Use a C² polynomial extension outside (ψ > ψ₀)
|
||||
|
||||
! The tail consists of the product p(ψ)⋅h(t), where:
|
||||
!
|
||||
! - p(ψ) is the 2nd order Taylor polynomial of the spline,
|
||||
! centered at ψ₀. See set_profiles_spline for details.
|
||||
!
|
||||
! - h(t) is a "smoothing" polynomial in the variable
|
||||
! t = (ψ - ψ₀)/(ψ₁ - ψ₀), defined as:
|
||||
!
|
||||
! h(t) = (1 - t)³⋅(1 + 3t + 6t²)
|
||||
!
|
||||
! with the following properties:
|
||||
!
|
||||
! h(0) = 1 h'(0)=0 h"(0)=0
|
||||
! h(1) = 0 h'(1)=0 h"(1)=0
|
||||
block
|
||||
real(wp_) :: dpsi, t, p, dp, h, dh
|
||||
dpsi = psin - self%tail%start ! Δψ = (ψ - ψ₀)
|
||||
|
||||
! Taylor polynomial p(ψ) and its derivative
|
||||
p = self%tail%value + dpsi*self%tail%deriv1 + dpsi**2*self%tail%deriv2/2
|
||||
dp = self%tail%deriv1 + dpsi*self%tail%deriv2
|
||||
|
||||
! Smoothing polynomial h(t) and its derivative
|
||||
t = dpsi/(self%tail%end - self%tail%start)
|
||||
h = (1 - t)**3 * (1 + 3*t + 6*t**2)
|
||||
dh = -30*(1 - t)**2 * t**2 / (self%tail%end - self%tail%start)
|
||||
|
||||
dens = p*h
|
||||
ddens = dp*h + p*dh
|
||||
end block
|
||||
end if
|
||||
|
||||
if (dens < 0) then
|
||||
write (msg, '("negative density:", 2(x,a,"=",g0.3))') &
|
||||
'ne', dens, 'ψ', psin
|
||||
call log_error(msg, mod='coreprofiles', proc='density')
|
||||
end if
|
||||
|
||||
end subroutine numeric_density
|
||||
|
||||
|
||||
function numeric_temp(self, psin) result(temp)
|
||||
! function arguments
|
||||
class(numeric_plasma), intent(in) :: self
|
||||
real(wp_), intent(in) :: psin
|
||||
real(wp_) :: temp
|
||||
|
||||
temp = 0
|
||||
|
||||
! Outside: temperature is undefined, return 0
|
||||
if (psin >= 1 .or. psin < 0) return
|
||||
|
||||
! Inside: use the interpolated numerical data
|
||||
temp = self%temp_spline%eval(psin)
|
||||
end function numeric_temp
|
||||
|
||||
|
||||
function numeric_zeff(self, psin) result(zeff)
|
||||
! function arguments
|
||||
class(numeric_plasma), intent(in) :: self
|
||||
real(wp_), intent(in) :: psin
|
||||
real(wp_) :: zeff
|
||||
|
||||
zeff = 1
|
||||
|
||||
! Outside: Zeff is undefined, return 1
|
||||
if (psin >= 1 .or. psin < 0) return
|
||||
|
||||
! Inside: use the interpolated numerical data
|
||||
zeff = self%zeff_spline%eval(psin)
|
||||
end function numeric_zeff
|
||||
|
||||
|
||||
subroutine load_plasma(params, plasma, err)
|
||||
use gray_params, only : gray_parameters
|
||||
use gray_params, only : PROF_ANALYTIC, PROF_NUMERIC
|
||||
use gray_params, only : RHO_TOR, RHO_POL, RHO_PSI
|
||||
use gray_params, only : SCALE_COLLISION, SCALE_GREENWALD, SCALE_OFF
|
||||
|
||||
use logger, only : log_error, log_debug
|
||||
use equilibrium, only : frhopol
|
||||
|
||||
! subroutine arguments
|
||||
type(gray_parameters), intent(inout) :: params
|
||||
class(abstract_plasma), allocatable, intent(out) :: plasma
|
||||
integer, intent(out) :: err
|
||||
|
||||
! local variables
|
||||
integer :: u ! unit number
|
||||
real(wp_) :: at, an, fact ! rescaling exponents and factor
|
||||
type(numeric_plasma) :: np
|
||||
type(analytic_plasma) :: ap
|
||||
|
||||
open(file=params%profiles%filenm, status='old', &
|
||||
action='read', newunit=u, iostat=err)
|
||||
if (err /= 0) then
|
||||
call log_error( &
|
||||
'opening profiles file ('//trim(params%profiles%filenm)//') failed!', &
|
||||
mod='coreprofiles', proc='load_plasma')
|
||||
err = 1
|
||||
return
|
||||
end if
|
||||
|
||||
! Set appropriate rescaling factors
|
||||
select case(params%profiles%iscal)
|
||||
case (SCALE_COLLISION)
|
||||
call log_debug('rescaling preserving ν*', mod='coreprofiles', proc='load_plasma')
|
||||
at = 2/3.0_wp_
|
||||
an = 4/3.0_wp_
|
||||
fact = params%equilibrium%factb
|
||||
case (SCALE_GREENWALD)
|
||||
call log_debug('rescaling preserving n_G', mod='coreprofiles', proc='load_plasma')
|
||||
at = 1
|
||||
an = 1
|
||||
fact = params%equilibrium%factb
|
||||
case (SCALE_OFF)
|
||||
at = 1
|
||||
an = 1
|
||||
fact = 1
|
||||
end select
|
||||
|
||||
select case(params%profiles%iprof)
|
||||
case (PROF_ANALYTIC)
|
||||
! The analytical file format is:
|
||||
!
|
||||
! 1 dens0 n1 n2
|
||||
! 2 te0 te1 t1 t2
|
||||
! 3 zeff
|
||||
call log_debug('loading analytical file', mod='coreprofiles', proc='load_plasma')
|
||||
|
||||
! Read model parameters
|
||||
read (u, *) ap%dens0, ap%n1, ap%n2
|
||||
read (u, *) ap%te0, ap%te1, ap%t1, ap%t2
|
||||
read (u, *) ap%z
|
||||
|
||||
! Apply rescaling factors
|
||||
ap%dens0 = ap%dens0 * fact**an * params%profiles%factne
|
||||
ap%te0 = ap%te0 * fact**at * params%profiles%factte
|
||||
ap%te1 = ap%te1 * fact**at * params%profiles%factte
|
||||
|
||||
allocate(plasma, source=ap)
|
||||
|
||||
case (PROF_NUMERIC)
|
||||
! The numerical file format is:
|
||||
!
|
||||
! 1 nrows
|
||||
! 2 ψ_n T_e n_e Z_eff
|
||||
! 3 ...
|
||||
call log_debug('loading numerical file', mod='coreprofiles', proc='load_plasma')
|
||||
|
||||
block
|
||||
integer :: i, nrows
|
||||
real(wp_), allocatable, dimension(:) :: psi, temp, dens, zeff
|
||||
|
||||
! Read header and allocate arrays
|
||||
read(u, *) nrows
|
||||
allocate(psi(nrows), temp(nrows), dens(nrows), zeff(nrows))
|
||||
|
||||
! Read the table rows
|
||||
do i=1, nrows
|
||||
read(u, *) psi(i), temp(i), dens(i), zeff(i)
|
||||
end do
|
||||
|
||||
! ??
|
||||
psi(1) = max(psi(1), zero)
|
||||
|
||||
! Apply rescaling factors
|
||||
dens = dens * fact**an * params%profiles%factne
|
||||
temp = temp * fact**at * params%profiles%factte
|
||||
|
||||
! Convert psi to ψ_n (normalised poloidal flux)
|
||||
select case (params%profiles%irho)
|
||||
case (RHO_TOR) ! psi is ρ_t = √Φ_n (toroidal flux)
|
||||
psi = [(frhopol(psi(i))**2, i=1,nrows)]
|
||||
case (RHO_POL) ! psi is ρ_p = √ψ_n (poloidal flux)
|
||||
psi = psi**2
|
||||
case (RHO_PSI) ! psi is already ψ_n
|
||||
end select
|
||||
|
||||
! Interpolate
|
||||
call np%init(params, psi, temp, dens, zeff, err)
|
||||
end block
|
||||
|
||||
allocate(plasma, source=np)
|
||||
|
||||
end select
|
||||
|
||||
close(u)
|
||||
|
||||
end subroutine load_plasma
|
||||
|
||||
|
||||
subroutine init_splines(self, params, psi, temp, dens, zeff, err)
|
||||
! Computes splines for the plasma profiles data and stores them
|
||||
! in their respective global variables, see the top of this file.
|
||||
!
|
||||
! `err` is 1 if I/O errors occured, 2 if other initialisation failed.
|
||||
use gray_params, only : gray_parameters
|
||||
use logger, only : log_debug, log_info, log_warning, log_error
|
||||
|
||||
! subroutine arguments
|
||||
class(numeric_plasma), intent(out) :: self
|
||||
type(gray_parameters), intent(inout) :: params
|
||||
real(wp_), dimension(:), intent(in) :: psi, temp, dens, zeff
|
||||
integer, intent(out) :: err
|
||||
|
||||
! local variables
|
||||
integer :: n
|
||||
|
||||
! for log messages formatting
|
||||
character(256) :: msg
|
||||
|
||||
n = size(psi)
|
||||
|
||||
! Spline interpolation of temperature and effective charge
|
||||
call self%temp_spline%init(psi, temp, n)
|
||||
call self%zeff_spline%init(psi, zeff, n)
|
||||
|
||||
! Spline interpolation of density (smooth spline)
|
||||
call self%dens_spline%init(psi, dens, n, range=[zero, psi(n)], &
|
||||
tension=params%profiles%sspld, err=err)
|
||||
|
||||
! Tension is very small, re-fit with an interpolating spline
|
||||
if (err == -1) then
|
||||
call log_warning('re-fitting with density curve with zero tension', &
|
||||
mod='coreprofiles', proc='init_splines')
|
||||
call self%dens_spline%init(psi, dens, n, range=[zero, psi(n)], tension=zero)
|
||||
else if (err > 0) then
|
||||
write (msg, '(a, g0)') 'density fit failed with error ', err
|
||||
call log_error(msg, mod='coreprofiles', proc='init_splines')
|
||||
err = 2
|
||||
return
|
||||
else
|
||||
err = 0
|
||||
end if
|
||||
|
||||
! Computation of the polynomial tail parameters
|
||||
!
|
||||
! Note: The density is the only quantity that needs to be evaluated
|
||||
! at the edge. The spline thus has to be extended to transition
|
||||
! smoothly from the last profile point to 0 outside the plasma.
|
||||
block
|
||||
real(wp_) :: s0, s1, s2 ! spline, 1st, 2nd derivative
|
||||
real(wp_) :: delta4 ! discriminant Δ/4 of q(x)
|
||||
real(wp_) :: x0, x1 ! vertex of q(x), solution
|
||||
|
||||
! Compute the coefficients of a 2nd order Taylor polinomial to
|
||||
! extend the spline beyond the last point:
|
||||
!
|
||||
! p(ψ) = s(ψ₀) + (ψ - ψ₀)s'(ψ₀) + ½(ψ - ψ₀)²s"(ψ₀)
|
||||
!
|
||||
! where s(ψ) is the spline and ψ₀ the last point.
|
||||
!
|
||||
s0 = self%dens_spline%eval(psi(n))
|
||||
s1 = self%dens_spline%deriv(psi(n), order=1)
|
||||
s2 = self%dens_spline%deriv(psi(n), order=2)
|
||||
|
||||
! Determine where to end the tail (to ensure the density remains
|
||||
! positive) from the zeros of the Taylor polynomial p(ψ)
|
||||
!
|
||||
! Define x=(ψ - ψ₀), then p(ψ)=0 is rewritten as
|
||||
!
|
||||
! q(x) = x² + 2s'/s" x + 2s/s" = 0
|
||||
!
|
||||
! The discriminant is Δ/4 = (s'/s")² - 2(s/s") and
|
||||
! the solutions are x = x₀ ± √(Δ/4), with x₀ = -s'/s".
|
||||
!
|
||||
x0 = -s1 / s2 ! vertex of parabola y=q(x)
|
||||
delta4 = (s1 / s2)**2 - 2*s0/s2 ! Δ/4 of q(x)
|
||||
|
||||
if (delta4 > 0) then
|
||||
! Pick the smallest positive zero (implying >ψ₀)
|
||||
x1 = x0 + sign(sqrt(delta4), sqrt(delta4) - x0)
|
||||
else
|
||||
! There are no zeros, use the parabola vertex
|
||||
x1 = x0
|
||||
call log_debug('spline extension has no zeros', &
|
||||
mod='coreprofiles', proc='init_splines')
|
||||
end if
|
||||
|
||||
! Store the tail parameters
|
||||
self%tail%start = psi(n)
|
||||
self%tail%end = self%tail%start + x1
|
||||
self%tail%value = s0
|
||||
self%tail%deriv1 = s1
|
||||
self%tail%deriv2 = s2
|
||||
end block
|
||||
|
||||
! Make sure the wave launcher does not fall inside the tail
|
||||
! Note: if it does, the initial wave conditions become
|
||||
! invalid as they are given assuming a vacuum (N=1)
|
||||
block
|
||||
use equilibrium, only : pol_flux
|
||||
use const_and_precisions, only : cm
|
||||
real(wp_) :: R, Z, psi
|
||||
|
||||
! Convert from cartesian to cylindrical coordinates
|
||||
associate (launch_pos => params%antenna%pos)
|
||||
R = hypot(launch_pos(1), launch_pos(2)) * cm ! R = √(x²+y²)
|
||||
z = launch_pos(3) * cm
|
||||
end associate
|
||||
|
||||
! Get the poloidal flux at the launcher
|
||||
! Note: this returns -1 when the data is not available
|
||||
call pol_flux(R, z, psi)
|
||||
|
||||
if (psi > self%tail%start .and. psi < self%tail%end) then
|
||||
! Fall back to the midpoint of ψ₀ and the launcher ψ
|
||||
call log_warning('downscaled tail to not reach the wave launcher', &
|
||||
mod='coreprofiles', proc='init_splines')
|
||||
write (msg, '(a, g0.3, a, g0.3, " → ", g0.3)') &
|
||||
'launcher: ψ=', psi, ' boundary: ψ=', self%tail%end, (self%tail%start + psi)/2
|
||||
call log_debug(msg, mod='coreprofiles', proc='init_splines')
|
||||
self%tail%end = (self%tail%start + psi)/2
|
||||
end if
|
||||
|
||||
if (psi > 0 .and. psi < self%tail%start) then
|
||||
! This must be a user error, stop here
|
||||
write (msg, '(a, a, g0.3, a, g0.3)') &
|
||||
'wave launcher is inside the plasma! ', &
|
||||
'launcher: ψ=', psi, ' boundary: ψ=', self%tail%end
|
||||
call log_error(msg, mod='coreprofiles', proc='init_splines')
|
||||
err = 2
|
||||
return
|
||||
end if
|
||||
end block
|
||||
|
||||
! Set the density boundary ψ
|
||||
! Note: this is used to detect entrance in the plasma
|
||||
params%profiles%psnbnd = self%tail%end
|
||||
|
||||
write (msg, '(a,g0.4)') 'density boundary: ψ=', self%tail%end
|
||||
call log_info(msg, mod='coreprofiles', proc='init_splines')
|
||||
end subroutine init_splines
|
||||
|
||||
end module gray_plasma
|
||||
@ -110,16 +110,17 @@ contains
|
||||
end subroutine init_tables
|
||||
|
||||
|
||||
function kinetic_profiles(params) result(tbl)
|
||||
function kinetic_profiles(params, plasma) result(tbl)
|
||||
! Generates the plasma kinetic profiles
|
||||
|
||||
use gray_params, only : gray_parameters, EQ_VACUUM
|
||||
use gray_plasma, only : abstract_plasma
|
||||
use equilibrium, only : fq, frhotor
|
||||
use coreprofiles, only : density, temp
|
||||
use magsurf_data, only : npsi, vajphiav
|
||||
|
||||
! function arguments
|
||||
type(gray_parameters), intent(in) :: params
|
||||
class(abstract_plasma), intent(in) :: plasma
|
||||
type(table) :: tbl
|
||||
|
||||
! local variables
|
||||
@ -151,8 +152,8 @@ contains
|
||||
else
|
||||
J_phi = 0
|
||||
end if
|
||||
call density(psi_n, dens, ddens)
|
||||
call tbl%append([psi_n, rho_t, dens, temp(psi_n), fq(psi_n), J_phi])
|
||||
call plasma%density(psi_n, dens, ddens)
|
||||
call tbl%append([psi_n, rho_t, dens, plasma%temp(psi_n), fq(psi_n), J_phi])
|
||||
end do
|
||||
|
||||
end function kinetic_profiles
|
||||
@ -226,17 +227,18 @@ contains
|
||||
end function ec_resonance
|
||||
|
||||
|
||||
function inputs_maps(params, B_res, X_norm) result(tbl)
|
||||
function inputs_maps(params, plasma, B_res, X_norm) result(tbl)
|
||||
! Generates 2D maps of several input quantities
|
||||
|
||||
use const_and_precisions, only : comp_eps, cm, degree
|
||||
use gray_params, only : gray_parameters, EQ_VACUUM
|
||||
use gray_plasma, only : abstract_plasma
|
||||
use equilibrium, only : rmnm, rmxm, zmnm, zmxm, pol_flux, bfield
|
||||
use coreprofiles, only : density, temp
|
||||
use magsurf_data, only : npsi
|
||||
|
||||
! function arguments
|
||||
type(gray_parameters), intent(in) :: params
|
||||
class(abstract_plasma), intent(in) :: plasma
|
||||
real(wp_), intent(in) :: B_res ! resonant magnetic field, e/m_eω
|
||||
real(wp_), intent(in) :: X_norm ! X normalisation, e²/ε₀m_eω²
|
||||
type(table) :: tbl
|
||||
@ -271,11 +273,11 @@ contains
|
||||
do k = 1, npsi
|
||||
call pol_flux(r(j), z(k), psi_n)
|
||||
call bfield(r(j), z(k), B_R, B_z, B_phi)
|
||||
call density(psi_n, ne, dne)
|
||||
call plasma%density(psi_n, ne, dne)
|
||||
B = sqrt(B_R**2 + B_phi**2 + B_z**2)
|
||||
X = X_norm*ne
|
||||
Y = B/B_res
|
||||
Te = temp(psi_n)
|
||||
Te = plasma%temp(psi_n)
|
||||
call tbl%append([R(j), z(k), psi_n, B_R, B_phi, &
|
||||
B_z, B, ne, Te, X, Y, Npl])
|
||||
end do
|
||||
|
||||
73
src/main.f90
73
src/main.f90
@ -5,6 +5,7 @@ program main
|
||||
use gray_cli, only : cli_options, parse_cli_options, &
|
||||
parse_param_overrides
|
||||
use gray_core, only : gray_main
|
||||
use gray_plasma, only : abstract_plasma, load_plasma
|
||||
use gray_params, only : gray_parameters, gray_data, gray_results, &
|
||||
read_gray_params, read_gray_config, &
|
||||
make_gray_header
|
||||
@ -18,6 +19,7 @@ program main
|
||||
! gray_main subroutine arguments
|
||||
type(gray_parameters) :: params ! Inputs
|
||||
type(gray_data) :: data !
|
||||
class(abstract_plasma), allocatable :: plasma !
|
||||
type(gray_results) :: results ! Outputs
|
||||
integer :: err ! Exit code
|
||||
|
||||
@ -81,8 +83,8 @@ program main
|
||||
call init_equilibrium(params, data, err)
|
||||
if (err /= 0) call exit(err)
|
||||
|
||||
call init_profiles(params%profiles, params%equilibrium%factb, &
|
||||
params%antenna%pos, data%profiles, err)
|
||||
! Read and initialise the plasma state object
|
||||
call load_plasma(params, plasma, err)
|
||||
if (err /= 0) call exit(err)
|
||||
|
||||
call init_misc(params, data)
|
||||
@ -98,7 +100,7 @@ program main
|
||||
params%misc%active_tables = opts%tables
|
||||
|
||||
! Finally run the simulation
|
||||
call gray_main(params, data, results, err)
|
||||
call gray_main(params, data, plasma, results, err)
|
||||
end if
|
||||
|
||||
print_res: block
|
||||
@ -136,10 +138,7 @@ program main
|
||||
|
||||
! Free memory
|
||||
free_mem: block
|
||||
use coreprofiles, only : unset_profiles_spline
|
||||
use equilibrium, only : unset_equil_spline
|
||||
|
||||
call unset_profiles_spline
|
||||
call unset_equil_spline
|
||||
end block free_mem
|
||||
|
||||
@ -199,68 +198,6 @@ contains
|
||||
end subroutine init_equilibrium
|
||||
|
||||
|
||||
subroutine init_profiles(params, factb, launch_pos, data, err)
|
||||
! Reads the plasma kinetic profiles file (containing the elecron
|
||||
! temperature, density and plasma effective charge) and initialises
|
||||
! the respective GRAY data structure.
|
||||
use gray_params, only : profiles_parameters, profiles_data, &
|
||||
RHO_TOR, RHO_POL, RHO_PSI, &
|
||||
PROF_ANALYTIC, PROF_NUMERIC
|
||||
use equilibrium, only : frhopol
|
||||
use coreprofiles, only : read_profiles_an, read_profiles, &
|
||||
scale_profiles, set_profiles_an, &
|
||||
set_profiles_spline
|
||||
use logger, only : log_debug
|
||||
|
||||
! subroutine arguments
|
||||
type(profiles_parameters), intent(inout) :: params
|
||||
real(wp_), intent(in) :: factb
|
||||
real(wp_), intent(in) :: launch_pos(3)
|
||||
type(profiles_data), intent(out) :: data
|
||||
integer, intent(out) :: err
|
||||
|
||||
! local variables
|
||||
integer :: i
|
||||
|
||||
select case (params%iprof)
|
||||
case (PROF_ANALYTIC)
|
||||
call log_debug('loading analytical file', &
|
||||
mod='main', proc='init_profiles')
|
||||
call read_profiles_an(params%filenm, data, err)
|
||||
if (err /= 0) return
|
||||
|
||||
case (PROF_NUMERIC)
|
||||
call log_debug('loading numerical file', &
|
||||
mod='main', proc='init_profiles')
|
||||
call read_profiles(params%filenm, data, err)
|
||||
if (err /= 0) return
|
||||
|
||||
! Convert psrad to ψ
|
||||
select case (params%irho)
|
||||
case (RHO_TOR) ! psrad is ρ_t = √Φ (toroidal flux)
|
||||
data%psrad = [(frhopol(data%psrad(i))**2, i=1,size(data%psrad))]
|
||||
case (RHO_POL) ! psrad is ρ_p = √ψ (poloidal flux)
|
||||
data%psrad = data%psrad**2
|
||||
case (RHO_PSI) ! psrad is already ψ
|
||||
end select
|
||||
end select
|
||||
|
||||
! Rescale input data
|
||||
call scale_profiles(params, factb, data)
|
||||
|
||||
! Set global variables
|
||||
select case (params%iprof)
|
||||
case (PROF_ANALYTIC)
|
||||
call set_profiles_an(params, data)
|
||||
case (PROF_NUMERIC)
|
||||
call log_debug('computing splines...', mod='main', proc='init_profiles')
|
||||
call set_profiles_spline(params, data, err, launch_pos)
|
||||
call log_debug('splines computed', mod='main', proc='init_profiles')
|
||||
end select
|
||||
|
||||
end subroutine init_profiles
|
||||
|
||||
|
||||
subroutine init_antenna(params, err)
|
||||
! Reads the wave launcher file (containing the wave frequency, launcher
|
||||
! position, direction and beam description) and initialises the respective
|
||||
|
||||
Loading…
Reference in New Issue
Block a user