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

replace coreprofiles module with an object

Browse Source 
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 2024-08-14 10:47:39 +02:00
parent b9ae6681b6
commit eb64803913
Signed by: rnhmjoj
GPG Key ID: BFBAF4C975F76450
6 changed files with 622 additions and 664 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

519
src/coreprofiles.f90
View File

@ -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

75
src/gray_core.f90
View File

@ -7,12 +7,12 @@ module gray_core
contains 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 const_and_precisions, only : zero, one, comp_tiny
use coreprofiles, only : temp, fzeff
use polarization, only : ellipse_to_field use polarization, only : ellipse_to_field
use types, only : table, wrap use types, only : table, wrap
use gray_params, only : gray_parameters, gray_data, gray_results, EQ_VACUUM 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, & use gray_tables, only : init_tables, store_ec_profiles, store_ray_data, &
store_beam_shape, find_flux_surfaces, & store_beam_shape, find_flux_surfaces, &
kinetic_profiles, ec_resonance, inputs_maps kinetic_profiles, ec_resonance, inputs_maps
@ -30,6 +30,7 @@ contains
! subroutine arguments ! subroutine arguments
type(gray_parameters), intent(inout) :: params type(gray_parameters), intent(inout) :: params
type(gray_data), intent(in) :: data type(gray_data), intent(in) :: data
class(abstract_plasma), intent(in) :: plasma
type(gray_results), intent(out) :: results type(gray_results), intent(out) :: results
! Predefined grid for the output profiles (optional) ! Predefined grid for the output profiles (optional)
@ -132,9 +133,9 @@ contains
! ========= set environment END ========= ! ========= set environment END =========
! Pre-determinted tables ! 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%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 ! print ψ surface for q=3/2 and q=2/1
call find_flux_surfaces( & call find_flux_surfaces( &
@ -257,8 +258,8 @@ contains
do jk=1,params%raytracing%nray do jk=1,params%raytracing%nray
if(iwait(jk)) cycle ! jk ray is waiting for next pass if(iwait(jk)) cycle ! jk ray is waiting for next pass
stv(jk) = stv(jk) + params%raytracing%dst ! current ray step stv(jk) = stv(jk) + params%raytracing%dst ! current ray step
call rkstep(params, sox, bres, xgcn, yw(:,jk), ypw(:,jk), & call rkstep(params, plasma, sox, bres, xgcn, yw(:,jk), &
gri(:,jk), ggri(:,:,jk), igrad_b) ypw(:,jk), gri(:,jk), ggri(:,:,jk), igrad_b)
end do end do
! update position and grad ! update position and grad
if(igrad_b == 1) call gradi_upd(params%raytracing, yw, ak0, xc, du1, gri, ggri) 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 ! compute derivatives with updated gradient and local plasma values
xv = yw(1:3,jk) xv = yw(1:3,jk)
anv = yw(4:6,jk) anv = yw(4:6,jk)
call ywppla_upd(params, xv, anv, gri(:,jk), ggri(:,:,jk), sox, bres, & call ywppla_upd(params, plasma, xv, anv, gri(:,jk), ggri(:,:,jk), &
xgcn,ypw(:,jk), psinv, dens, btot, bv, xg, yg, derxg, & sox, bres, xgcn,ypw(:,jk), psinv, dens, btot, bv, &
anpl, anpr, ddr, ddi, dersdst, derdnm, ierrn, igrad_b) xg, yg, derxg, anpl, anpr, ddr, ddi, dersdst, &
derdnm, ierrn, igrad_b)
! update global error code and print message ! update global error code and print message
if(ierrn/=0) then if(ierrn/=0) then
error = ior(error,ierrn) error = ior(error,ierrn)
@ -372,9 +374,9 @@ contains
yw(:,jk) = [xv, anv] ! * updated coordinates (reflected) yw(:,jk) = [xv, anv] ! * updated coordinates (reflected)
igrad_b = 0 ! * switch to ray-tracing igrad_b = 0 ! * switch to ray-tracing
call ywppla_upd(params, xv, anv, gri(:,jk), ggri(:,:,jk), sox, & call ywppla_upd(params, plasma, xv, anv, gri(:,jk), ggri(:,:,jk), &
bres, xgcn, ypw(:,jk), psinv, dens, btot, bv, & sox, bres, xgcn, ypw(:,jk), psinv, dens, btot, &
xg, yg, derxg, anpl, anpr, ddr, ddi, dersdst, & bv, xg, yg, derxg, anpl, anpr, ddr, ddi, dersdst, &
derdnm, ierrn, igrad_b) ! * update derivatives after reflection derdnm, ierrn, igrad_b) ! * update derivatives after reflection
if(ierrn/=0) then ! * update global error code and print message if(ierrn/=0) then ! * update global error code and print message
error = ior(error,ierrn) error = ior(error,ierrn)
@ -429,12 +431,13 @@ contains
if(ierrn==0 .and. params%ecrh_cd%iwarm>0 .and. ins_pl .and. & 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 (tau1(jk)+tau0(jk)+lgcpl1(jk))<=taucr .and. .not.iwait(jk)) then ! H&CD computation check
tekev=temp(psinv) tekev = plasma%temp(psinv)
block block
complex(wp_) :: Npr_warm complex(wp_) :: Npr_warm
call alpha_effj(params%ecrh_cd, psinv, xg, yg, dens, tekev, & call alpha_effj(params%ecrh_cd, plasma, psinv, xg, yg, dens, &
ak0, bres, derdnm, anpl, anpr, sox, Npr_warm, & tekev, ak0, bres, derdnm, anpl, anpr, sox, &
alpha, didp, nharm, nhf, iokhawa, ierrhcd) Npr_warm, alpha, didp, nharm, nhf, iokhawa, &
ierrhcd)
anprre = Npr_warm%re anprre = Npr_warm%re
anprim = Npr_warm%im anprim = Npr_warm%im
if(ierrhcd /= 0) then 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 ! Runge-Kutta integrator
use gray_params, only : gray_parameters use gray_params, only : gray_parameters
use gray_plasma, only : abstract_plasma
! subroutine arguments ! subroutine arguments
type(gray_parameters), intent(in) :: params type(gray_parameters), intent(in) :: params
class(abstract_plasma), intent(in) :: plasma
real(wp_), intent(in) :: bres, xgcn real(wp_), intent(in) :: bres, xgcn
real(wp_), intent(inout) :: y(6) real(wp_), intent(inout) :: y(6)
real(wp_), intent(in) :: yp(6) real(wp_), intent(in) :: yp(6)
@ -1024,18 +1030,20 @@ contains
function f(y) function f(y)
real(wp_), intent(in) :: y(6) real(wp_), intent(in) :: y(6)
real(wp_) :: f(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 function
end subroutine rkstep 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) ! Compute right-hand side terms of the ray equations (dery)
! used in R-K integrator ! used in R-K integrator
use gray_params, only : gray_parameters use gray_params, only : gray_parameters
use gray_plasma, only : abstract_plasma
! subroutine arguments ! subroutine arguments
type(gray_parameters), intent(in) :: params type(gray_parameters), intent(in) :: params
class(abstract_plasma), intent(in) :: plasma
real(wp_), intent(in) :: y(6) real(wp_), intent(in) :: y(6)
real(wp_), intent(in) :: bres, xgcn real(wp_), intent(in) :: bres, xgcn
real(wp_), intent(in) :: dgr(3) real(wp_), intent(in) :: dgr(3)
@ -1049,7 +1057,7 @@ contains
real(wp_), dimension(3,3) :: derbv real(wp_), dimension(3,3) :: derbv
xv = y(1:3) 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) bv, derbv, xg, yg, derxg, deryg)
anv = y(4:6) anv = y(4:6)
call disp_deriv(params, anv, sox, xg, yg, derxg, deryg, & call disp_deriv(params, anv, sox, xg, yg, derxg, deryg, &
@ -1057,16 +1065,18 @@ contains
end subroutine rhs 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, & psinv, dens, btot, bv, xg, yg, derxg, anpl, anpr, &
ddr, ddi, dersdst, derdnm, error, igrad) ddr, ddi, dersdst, derdnm, error, igrad)
! Compute right-hand side terms of the ray equations (dery) ! Compute right-hand side terms of the ray equations (dery)
! used after full R-K step and grad(S_I) update ! used after full R-K step and grad(S_I) update
use gray_errors, only : raise_error, large_npl use gray_errors, only : raise_error, large_npl
use gray_params, only : gray_parameters use gray_params, only : gray_parameters
use gray_plasma, only : abstract_plasma
! subroutine rguments ! subroutine rguments
type(gray_parameters), intent(in) :: params type(gray_parameters), intent(in) :: params
class(abstract_plasma), intent(in) :: plasma
real(wp_), intent(in) :: xv(3), anv(3) real(wp_), intent(in) :: xv(3), anv(3)
real(wp_), intent(in) :: dgr(3) real(wp_), intent(in) :: dgr(3)
@ -1086,7 +1096,7 @@ contains
real(wp_), dimension(3,3) :: derbv real(wp_), dimension(3,3) :: derbv
real(wp_), parameter :: anplth1 = 0.99_wp_, anplth2 = 1.05_wp_ 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, & call disp_deriv(params, anv,sox,xg,yg,derxg,deryg,bv,derbv,dgr,ddgr,igrad, &
dery,anpl,anpr,ddr,ddi,dersdst,derdnm) dery,anpl,anpr,ddr,ddi,dersdst,derdnm)
@ -1305,15 +1315,16 @@ contains
end subroutine solg3 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) derbv, xg, yg, derxg, deryg, psinv_opt)
use const_and_precisions, only : zero, cm use const_and_precisions, only : zero, cm
use gray_params, only : gray_parameters, EQ_VACUUM use gray_params, only : gray_parameters, EQ_VACUUM
use gray_plasma, only : abstract_plasma
use equilibrium, only : psia, pol_flux, pol_curr, sgnbphi use equilibrium, only : psia, pol_flux, pol_curr, sgnbphi
use coreprofiles, only : density
! subroutine arguments ! subroutine arguments
type(gray_parameters), intent(in) :: params type(gray_parameters), intent(in) :: params
class(abstract_plasma), intent(in) :: plasma
real(wp_), dimension(3), intent(in) :: xv real(wp_), dimension(3), intent(in) :: xv
real(wp_), intent(in) :: xgcn, bres real(wp_), intent(in) :: xgcn, bres
real(wp_), intent(out) :: dens, btot, xg, yg real(wp_), intent(out) :: dens, btot, xg, yg
@ -1385,7 +1396,7 @@ contains
end if end if
! compute xg and derivative ! compute xg and derivative
call density(psinv, dens, ddensdpsi) call plasma%density(psinv, dens, ddensdpsi)
xg = xgcn*dens xg = xgcn*dens
dxgdpsi = xgcn*ddensdpsi dxgdpsi = xgcn*ddensdpsi
@ -1735,14 +1746,14 @@ contains
subroutine alpha_effj(params, psinv, X, Y, density, temperature, & subroutine alpha_effj(params, plasma, psinv, X, Y, density, temperature, &
k0, Bres, derdnm, Npl, Npr_cold, sox, & k0, Bres, derdnm, Npl, Npr_cold, sox, Npr, &
Npr, alpha, dIdp, nhmin, nhmax, iokhawa, error) alpha, dIdp, nhmin, nhmax, iokhawa, error)
! Computes the absorption coefficient and effective current ! Computes the absorption coefficient and effective current
use const_and_precisions, only : pi, mc2=>mc2_ use const_and_precisions, only : pi, mc2=>mc2_
use gray_params, only : ecrh_cd_parameters use gray_params, only : ecrh_cd_parameters
use coreprofiles, only : fzeff use gray_plasma, only : abstract_plasma
use equilibrium, only : sgnbphi use equilibrium, only : sgnbphi
use dispersion, only : harmnumber, warmdisp use dispersion, only : harmnumber, warmdisp
use eccd, only : setcdcoeff, eccdeff, fjch0, fjch, fjncl use eccd, only : setcdcoeff, eccdeff, fjch0, fjch, fjncl
@ -1755,7 +1766,9 @@ contains
! Inputs ! Inputs
! ECRH & CD parameters ! 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 ψ ! poloidal flux ψ
real(wp_), intent(in) :: psinv real(wp_), intent(in) :: psinv
! CMA diagram variables: X=(ω_pe/ω)², Y=ω_ce/ω ! CMA diagram variables: X=(ω_pe/ω)², Y=ω_ce/ω
@ -1855,7 +1868,7 @@ contains
! Current drive computation ! Current drive computation
if (params%ieccd <= 0) return if (params%ieccd <= 0) return
zeff = fzeff(psinv) zeff = plasma%zeff(psinv)
ithn = 1 ithn = 1
if (nlarmor > nhmin) ithn = 2 if (nlarmor > nhmin) ithn = 2
rhop = sqrt(psinv) rhop = sqrt(psinv)

23
src/gray_jetto1beam.f90
View File

@ -4,7 +4,9 @@ subroutine gray_jetto1beam(ijetto, mr, mz, r, z, psin, psia, rax, zax, &
use const_and_precisions, only: wp_ use const_and_precisions, only: wp_
use gray_params, only: gray_parameters, gray_data, gray_results use gray_params, only: gray_parameters, gray_data, gray_results
use gray_core, only: gray_main use gray_core, only: gray_main
use gray_plasma, only: numeric_plasma
implicit none
! subroutine arguments ! subroutine arguments
integer, intent(in) :: ijetto, mr, mz, nbnd, nrho, ibeam 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 ! local variables
type(gray_parameters) :: params type(gray_parameters) :: params
type(gray_data) :: data type(gray_data) :: data
type(numeric_plasma) :: plasma
type(gray_results) :: res type(gray_results) :: res
logical, save :: firstcall = .true. logical, save :: firstcall = .true.
@ -86,19 +89,9 @@ subroutine gray_jetto1beam(ijetto, mr, mz, r, z, psin, psia, rax, zax, &
if (err /= 0) return if (err /= 0) return
end block init_equilibrium end block init_equilibrium
! Set plasma kinetic profiles ! Compute splines of kinetic profiles
init_profiles: block call plasma%init(params, psrad, te, dne, zeff, err)
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)
if (err /= 0) return if (err /= 0) return
end block init_profiles
! Set wave launcher parameters ! Set wave launcher parameters
init_antenna: block init_antenna: block
@ -115,7 +108,7 @@ subroutine gray_jetto1beam(ijetto, mr, mz, r, z, psin, psia, rax, zax, &
end block init_antenna end block init_antenna
! Call main subroutine for the ibeam-th beam ! 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 write_debug_files: block
integer :: i, err integer :: i, err
@ -130,13 +123,9 @@ subroutine gray_jetto1beam(ijetto, mr, mz, r, z, psin, psia, rax, zax, &
! Free memory ! Free memory
free_memory: block free_memory: block
use equilibrium, only : unset_equil_spline use equilibrium, only : unset_equil_spline
use coreprofiles, only : unset_profiles_spline
! Unset global variables of the `equilibrium` module ! Unset global variables of the `equilibrium` module
call unset_equil_spline call unset_equil_spline
! Unset global variables of the `coreprofiles` module
call unset_profiles_spline
end block free_memory end block free_memory
! Copy over the results ! Copy over the results

536
src/gray_plasma.f90 Normal file
View 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

18
src/gray_tables.f90
View File

@ -110,16 +110,17 @@ contains
end subroutine init_tables end subroutine init_tables
function kinetic_profiles(params) result(tbl) function kinetic_profiles(params, plasma) result(tbl)
! Generates the plasma kinetic profiles ! Generates the plasma kinetic profiles
use gray_params, only : gray_parameters, EQ_VACUUM use gray_params, only : gray_parameters, EQ_VACUUM
use gray_plasma, only : abstract_plasma
use equilibrium, only : fq, frhotor use equilibrium, only : fq, frhotor
use coreprofiles, only : density, temp
use magsurf_data, only : npsi, vajphiav use magsurf_data, only : npsi, vajphiav
! function arguments ! function arguments
type(gray_parameters), intent(in) :: params type(gray_parameters), intent(in) :: params
class(abstract_plasma), intent(in) :: plasma
type(table) :: tbl type(table) :: tbl
! local variables ! local variables
@ -151,8 +152,8 @@ contains
else else
J_phi = 0 J_phi = 0
end if end if
call density(psi_n, dens, ddens) call plasma%density(psi_n, dens, ddens)
call tbl%append([psi_n, rho_t, dens, temp(psi_n), fq(psi_n), J_phi]) call tbl%append([psi_n, rho_t, dens, plasma%temp(psi_n), fq(psi_n), J_phi])
end do end do
end function kinetic_profiles end function kinetic_profiles
@ -226,17 +227,18 @@ contains
end function ec_resonance 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 ! Generates 2D maps of several input quantities
use const_and_precisions, only : comp_eps, cm, degree use const_and_precisions, only : comp_eps, cm, degree
use gray_params, only : gray_parameters, EQ_VACUUM 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 equilibrium, only : rmnm, rmxm, zmnm, zmxm, pol_flux, bfield
use coreprofiles, only : density, temp
use magsurf_data, only : npsi use magsurf_data, only : npsi
! function arguments ! function arguments
type(gray_parameters), intent(in) :: params 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) :: B_res ! resonant magnetic field, e/m_eω
real(wp_), intent(in) :: X_norm ! X normalisation, e²/ε₀m_eω² real(wp_), intent(in) :: X_norm ! X normalisation, e²/ε₀m_eω²
type(table) :: tbl type(table) :: tbl
@ -271,11 +273,11 @@ contains
do k = 1, npsi do k = 1, npsi
call pol_flux(r(j), z(k), psi_n) call pol_flux(r(j), z(k), psi_n)
call bfield(r(j), z(k), B_R, B_z, B_phi) 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) B = sqrt(B_R**2 + B_phi**2 + B_z**2)
X = X_norm*ne X = X_norm*ne
Y = B/B_res 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, & call tbl%append([R(j), z(k), psi_n, B_R, B_phi, &
B_z, B, ne, Te, X, Y, Npl]) B_z, B, ne, Te, X, Y, Npl])
end do end do

73
src/main.f90
View File

@ -5,6 +5,7 @@ program main
use gray_cli, only : cli_options, parse_cli_options, & use gray_cli, only : cli_options, parse_cli_options, &
parse_param_overrides parse_param_overrides
use gray_core, only : gray_main use gray_core, only : gray_main
use gray_plasma, only : abstract_plasma, load_plasma
use gray_params, only : gray_parameters, gray_data, gray_results, & use gray_params, only : gray_parameters, gray_data, gray_results, &
read_gray_params, read_gray_config, & read_gray_params, read_gray_config, &
make_gray_header make_gray_header
@ -18,6 +19,7 @@ program main
! gray_main subroutine arguments ! gray_main subroutine arguments
type(gray_parameters) :: params ! Inputs type(gray_parameters) :: params ! Inputs
type(gray_data) :: data ! type(gray_data) :: data !
class(abstract_plasma), allocatable :: plasma !
type(gray_results) :: results ! Outputs type(gray_results) :: results ! Outputs
integer :: err ! Exit code integer :: err ! Exit code
@ -81,8 +83,8 @@ program main
call init_equilibrium(params, data, err) call init_equilibrium(params, data, err)
if (err /= 0) call exit(err) if (err /= 0) call exit(err)
call init_profiles(params%profiles, params%equilibrium%factb, & ! Read and initialise the plasma state object
params%antenna%pos, data%profiles, err) call load_plasma(params, plasma, err)
if (err /= 0) call exit(err) if (err /= 0) call exit(err)
call init_misc(params, data) call init_misc(params, data)
@ -98,7 +100,7 @@ program main
params%misc%active_tables = opts%tables params%misc%active_tables = opts%tables
! Finally run the simulation ! Finally run the simulation
call gray_main(params, data, results, err) call gray_main(params, data, plasma, results, err)
end if end if
print_res: block print_res: block
@ -136,10 +138,7 @@ program main
! Free memory ! Free memory
free_mem: block free_mem: block
use coreprofiles, only : unset_profiles_spline
use equilibrium, only : unset_equil_spline use equilibrium, only : unset_equil_spline
call unset_profiles_spline
call unset_equil_spline call unset_equil_spline
end block free_mem end block free_mem
@ -199,68 +198,6 @@ contains
end subroutine init_equilibrium 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) subroutine init_antenna(params, err)
! Reads the wave launcher file (containing the wave frequency, launcher ! Reads the wave launcher file (containing the wave frequency, launcher
! position, direction and beam description) and initialises the respective ! position, direction and beam description) and initialises the respective