From 45ef9c5eae24c475c30770d3bdb2a7740e174eba Mon Sep 17 00:00:00 2001 From: Michele Guerini Rocco Date: Sat, 21 May 2022 22:56:57 +0200 Subject: [PATCH] src/coreprofiles: make psnbnd fully automatic 1. Fix the mismatch between the psnbnd in coreprofiles and gray_core. This happens whenever gray overrides the externally provided one (i.e. the density tail would become negative before psnbnd and is so rescaled to end exactly on the zero). 2. Make psnbnd no longer required by always computing it as in 1. It hasn't been removed, because gray_params.data is sacrosant, but it no longer has any effect. 3. Cleanup: mark public functions, restructure the global variables into three categories; add comments explaining the analytical profiles format, formulae and how the polynomial tail is computed. --- input/gray.ini | 6 - src/coreprofiles.f90 | 546 +++++++++++++++++++++++++++++-------------- src/main.f90 | 16 +- 3 files changed, 375 insertions(+), 193 deletions(-) diff --git a/input/gray.ini b/input/gray.ini index 1a2dd97..60a0aee 100644 --- a/input/gray.ini +++ b/input/gray.ini @@ -145,12 +145,6 @@ irho = 0 ; Filepath of the equilibrium (relative to this file) filenm = "profiles.txt" -; Value of ψ at the plasma boundary [multipass module] -; Notes: -; 1. boundary means zero density; -; 2. determines when a ray is considered inside the plasma -psnbnd = 1.007 - ; Tension of the density spline ; Note: 0 means perfect interpolation sspld = 0.1 diff --git a/src/coreprofiles.f90 b/src/coreprofiles.f90 index 615cc11..888b132 100644 --- a/src/coreprofiles.f90 +++ b/src/coreprofiles.f90 @@ -1,136 +1,228 @@ +! This modules 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 + implicit none - integer, save :: npp, nsfd - real(wp_), save :: psdbnd, psnpp, denpp, ddenpp, d2denpp - real(wp_), dimension(:), allocatable, save :: tfn, cfn, psrad - real(wp_), dimension(:, :), allocatable, save :: ct, cz - real(wp_), save :: dens0, aln1, aln2, te0, dte0, alt1, alt2, zeffan + ! Parameters of the plasma profiles splines + type spline_parameters + integer :: ndata ! Number of data points + integer :: nknots ! Number of spline knots + ! Density spline (ψ, knots, B-spline coefficients) + real(wp_), dimension(:), allocatable :: knots, coeffs + ! Temperature and effective charge arrays (ψ, T(ψ), Zeff(ψ)) + real(wp_), dimension(:), allocatable :: psi + real(wp_), dimension(:, :), allocatable :: temp, 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 + + ! 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 variable storing the state of the module + type(spline_parameters), save :: spline + type(density_tail), save :: tail + type(analytic_model), save :: model + + 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) + 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 gray_params, only : iprof - use dierckx, only : splev,splder + use dierckx, only : splev, splder use logger, only : log_error implicit none ! subroutine arguments - real(wp_), intent(in) :: psin - real(wp_), intent(out) :: dens,ddens + real(wp_), intent(in) :: psin ! normalised poloidal flux + real(wp_), intent(out) :: dens, ddens ! density and first derivative ! local variables - integer :: ier,nu - real(wp_) :: profd,dprofd,dpsib,tt,fp,dfp,fh,dfh - real(wp_), dimension(1) :: xxs,ffs - real(wp_), dimension(npp+4) :: wrkfd - character(256) :: msg + integer :: ier ! dierck error code + real(wp_) :: f(1) ! dierck output (must be an array) + real(wp_) :: wrkfd(spline%ndata+4) ! dierck working space array + character(256) :: msg ! for log messages formatting - ! - ! Computation of density [10¹⁹ m⁻³] and derivative wrt ψ - ! - dens=zero - ddens=zero - if((psin >= psdbnd) .or. (psin < zero)) return + ! Initialise both to zero + dens = zero + ddens = zero - if(iprof == 0) then - if(psin > one) return - profd=(one-psin**aln1)**aln2 - dens=dens0*profd - dprofd=-aln1*aln2*psin**(aln1-one) & - *(one-psin**aln1)**(aln2-one) - ddens=dens0*dprofd - else - if(psin > psnpp) then + ! Outside the tail end both density and its + ! derivatives are identically zero + if (psin >= tail%end .or. psin < 0) return - ! Smooth interpolation for psnpp < psi < psdbnd - ! dens = fp * fh - ! fp: parabola matched at psi=psnpp with given profile density - ! fh=(1-t)^3(1+3t+6t^2) is a smoothing function: - ! fh(0)=1, fh(1)=0 and zero first and second deriv at t=0,1 + if (iprof == 0) then + ! Use the analytical model ! - dpsib=psin-psnpp - fp=denpp+dpsib*ddenpp+0.5_wp_*dpsib**2*d2denpp - dfp=ddenpp+dpsib*d2denpp - tt=dpsib/(psdbnd-psnpp) - fh=(one-tt)**3*(one+3.0_wp_*tt+6.0_wp_*tt**2) - dfh=-30.0_wp_*(one-tt)**2*tt**2/(psdbnd-psnpp) - dens=fp*fh - ddens=dfp*fh+fp*dfh + ! 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 + + if (psin < tail%start) then + ! Use the interpolating spline when in range + + ! Evaluate the spline + ier = 0 + call splev(spline%knots, spline%nknots, spline%coeffs, & + 3, [psin], f, 1, ier) + dens = f(1) + + ! Evaluate the spline 1st derivative + ier = 0 + call splder(spline%knots, spline%nknots, spline%coeffs, & + 3, 1, [psin], f, 1, wrkfd, ier) + ddens = f(1) + if (abs(dens) < 1.0e-10_wp_) dens = zero else - xxs(1)=psin - ier=0 - call splev(tfn,nsfd,cfn,3,xxs,ffs,1,ier) - dens=ffs(1) - nu=1 - ier=0 - call splder(tfn,nsfd,cfn,3,nu,xxs,ffs,1,wrkfd,ier) - ddens=ffs(1) - if(abs(dens) < 1.0e-10_wp_) dens=zero + ! 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 < zero) then + + 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) - use const_and_precisions, only : wp_,zero,one - use gray_params, only : iprof - use utils, only : locate - use simplespline, only :spli + ! Computes the temperature as a function of the + ! normalised poloidal flux. + ! + ! Note: temperature has units of keV. + use gray_params, only : iprof + use utils, only : locate + use simplespline, only : spli + implicit none -! arguments + + ! subroutine arguments real(wp_), intent(in) :: psin real(wp_) :: temp -! local variables - integer :: k - real(wp_) :: proft,dps - temp=zero - if((psin >= one).or.(psin < zero)) return - if(iprof == 0) then - proft=(1.0_wp_-psin**alt1)**alt2 - temp=(te0-dte0)*proft+dte0 + ! local variables + integer :: k + real(wp_) :: proft, dps + + temp = zero + 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 - call locate(psrad,npp,psin,k) - k=max(1,min(k,npp-1)) - dps=psin-psrad(k) - temp=spli(ct,npp,k,dps) + ! Use the interpolated numerical data + call locate(spline%psi, spline%ndata, psin, k) + k = max(1, min(k, spline%ndata - 1)) + dps = psin - spline%psi(k) + temp = spli(spline%temp, spline%ndata, k, dps) endif end function temp + function fzeff(psin) - use const_and_precisions, only : wp_,zero,one - use gray_params, only : iprof - use utils, only : locate - use simplespline, only :spli + ! Computes the effective charge Zeff as a + ! function of the normalised poloidal flux. + use gray_params, only : iprof + use utils, only : locate + use simplespline, only : spli + implicit none -! arguments + + ! subroutine arguments real(wp_), intent(in) :: psin real(wp_) :: fzeff -! local variables + + ! local variables integer :: k real(wp_) :: dps - fzeff=one - if((psin >= one).or.(psin < zero)) return - if(iprof == 0) then - fzeff=zeffan + fzeff = one + if (psin >= 1 .or. psin < 0) return + if (iprof == 0) then + ! Use the analytical model (just a constant) + fzeff = model%zeff else - call locate(psrad,npp,psin,k) - k=max(1,min(k,npp-1)) - dps=psin-psrad(k) - fzeff=spli(cz,npp,k,dps) + ! Use the interpolated numerical data + call locate(spline%psi, spline%ndata, psin, k) + k = max(1, min(k, spline%ndata - 1)) + dps = psin - spline%psi(k) + fzeff = spli(spline%zeff, spline%ndata, k, dps) endif end function fzeff - + subroutine read_profiles(filenm, data, unit) ! Reads the radial plasma profiles from `file` and store them ! into `data`. If given, the file is opened in the `unit` number. @@ -154,10 +246,10 @@ contains integer :: err ! 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) + 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) @@ -190,7 +282,14 @@ contains ! from params%filenm. ! If given, the file is opened in the `unit` number. ! - ! TODO: add format description + ! 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 use utils, only : get_free_unit use logger, only : log_error @@ -205,7 +304,7 @@ contains ! local variables integer :: u integer :: err - + u = get_free_unit(unit) if (allocated(data%terad)) deallocate(data%terad) @@ -220,9 +319,9 @@ contains call exit(1) end if - read (u,*) data%derad(1:3) ! dens0, aln1, aln2 - read (u,*) data%terad(1:4) ! te0, dte0, alt1, alt2 - read (u,*) data%zfc(1) ! zeffan + 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 @@ -253,13 +352,13 @@ contains aan = one end if - if(params%iscal==2) then + if (params%iscal==2) then ffact = one else ffact = factb end if - if(params%iprof==0) then + if (params%iprof==0) then last_te = 2 last_ne = 1 else @@ -272,140 +371,227 @@ contains end subroutine scale_profiles - subroutine set_profiles_spline(params, data) + subroutine set_profiles_spline(params, data, 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. use simplespline, only : difcs use dierckx, only : curfit, splev, splder use gray_params, only : profiles_parameters, profiles_data - use logger, only : log_info, log_warning + use logger, only : log_debug, log_info, log_warning, log_error implicit none ! subroutine arguments - type(profiles_parameters), intent(in) :: params + type(profiles_parameters), intent(inout) :: params type(profiles_data), intent(inout) :: data + real(wp_), optional, intent(in) :: launch_pos(3) + + ! curfit parameters + integer, parameter :: iopt = 0 ! smoothing spline mode + integer, parameter :: kspl = 3 ! order of spline (cubic) ! local variables - integer, parameter :: iopt=0, kspl=3 - integer :: n, npest, lwrkf, ier - real(wp_) :: xb, xe, fp, xnv, xxp,xxm,delta2,ssplne_loc - real(wp_), dimension(:), allocatable :: wf, wrkf - integer, dimension(:), allocatable :: iwrkf - real(wp_), dimension(1) :: dedge,ddedge,d2dedge - character(256) :: msg ! for log messages formatting + integer :: n, npest, ier + real(wp_) :: xb, xe, fp, ssplne_loc - n=size(data%psrad) - npest=n+4 - lwrkf=n*4+npest*16 - allocate(wrkf(lwrkf),iwrkf(npest),wf(n)) + ! working space arrays for the dierckx functions + integer :: lwrkf + real(wp_), dimension(:), allocatable :: wf, wrkf + integer, dimension(:), allocatable :: iwrkf + + ! for log messages formatting + character(256) :: msg + + n = size(data%psrad) + npest = n + 4 + lwrkf = n*4 + npest*16 + allocate(wrkf(lwrkf), iwrkf(npest), wf(n)) ssplne_loc=params%sspld - ! if necessary, reallocate spline arrays - if(.not.allocated(psrad)) then - allocate(psrad(n),ct(n,4),cz(n,4)) + ! If necessary, reallocate the spline arrays + if (.not. allocated(spline%psi)) then + allocate(spline%psi(n), spline%temp(n, 4), spline%zeff(n, 4)) else - if(size(psrad) psnpp) psdbnd=min(psdbnd,xnv) - else - xxm=xnv-sqrt(delta2) - xxp=xnv+sqrt(delta2) - if(xxm > psnpp) then - psdbnd=min(psdbnd,xxm) - else if (xxp > psnpp) then - psdbnd=min(psdbnd,xxp) + ! 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(1) > 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 - write (msg, '(a,g0.3)') 'density boundary: ψ=', psdbnd - call log_info(msg, mod="coreprofiles", proc="set_profiles_spline") + + ! Store the tail parameters + tail%start = data%psrad(n) + tail%end = tail%start + x1(1) + tail%value = s0(1) + tail%deriv1 = s1(1) + tail%deriv2 = s2(1) + 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 : equinum_psi + real(wp_) :: R, Z, psi + real(wp_), parameter :: cm = 1.0e-2_wp_ + + ! 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 equinum_psi(R, z, psi) + + if (psi > tail%start .and. psi < tail%end) then + ! Fall back to the midpoint of ψ₀ and the launcher ψ + tail%end = (tail%start + psi)/2 + call log_warning('downscaled tail to not reach the wave launcher', & + mod='coreprofiles', proc='set_profiles_spline') + 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') + call exit(2) + end if + end block end if - deallocate(iwrkf,wrkf,wf) + ! 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') + + deallocate(iwrkf, wrkf, wf) end subroutine set_profiles_spline subroutine unset_profiles_spline + ! Unsets the splines global variables, see the top of this file. + implicit none - if (allocated(psrad)) deallocate(psrad) - if (allocated(ct)) deallocate(ct) - if (allocated(cz)) deallocate(cz) - if (allocated(tfn)) deallocate(tfn) - if (allocated(cfn)) deallocate(cfn) + if (allocated(spline%psi)) deallocate(spline%psi) + if (allocated(spline%temp)) deallocate(spline%temp) + if (allocated(spline%zeff)) deallocate(spline%zeff) + if (allocated(spline%knots)) deallocate(spline%knots) + if (allocated(spline%coeffs)) deallocate(spline%coeffs) end subroutine unset_profiles_spline - subroutine set_profiles_an(data) + 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_data + use gray_params, only : profiles_parameters, profiles_data implicit none ! subroutine arguments - type(profiles_data), intent(in) :: data + 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 - te0 = data%terad(1) - dte0 = data%terad(2) - alt1 = data%terad(3) - alt2 = data%terad(4) - dens0 = data%derad(1) - aln1 = data%derad(2) - aln2 = data%derad(3) - zeffan = data%zfc(1) - psdbnd = one end subroutine set_profiles_an - + end module coreprofiles diff --git a/src/main.f90 b/src/main.f90 index e46e96d..275b586 100644 --- a/src/main.f90 +++ b/src/main.f90 @@ -63,8 +63,9 @@ program main ! Read the input data and set the global variables ! of the respective module. Note: order matters. call init_equilibrium(params, data) - call init_profiles(params%profiles, params%equilibrium%factb, data%profiles) call init_antenna(params%antenna) + call init_profiles(params%profiles, params%equilibrium%factb, & + params%antenna%pos, data%profiles) call init_misc(params, data) ! Change the current directory to output files there @@ -241,7 +242,7 @@ contains end subroutine deinit_equilibrium - subroutine init_profiles(params, factb, data) + subroutine init_profiles(params, factb, launch_pos, data) ! Reads the plasma kinetic profiles file (containing the elecron ! temperature, density and plasma effective charge) and initialises ! the respective GRAY data structure. @@ -253,9 +254,10 @@ contains implicit none ! subroutine arguments - type(profiles_parameters), intent(in) :: params - real(wp_), intent(in) :: factb - type(profiles_data), intent(out) :: data + type(profiles_parameters), intent(inout) :: params + real(wp_), intent(in) :: factb + real(wp_), intent(in) :: launch_pos(3) + type(profiles_data), intent(out) :: data if (params%iprof == 0) then ! Analytical profiles @@ -280,10 +282,10 @@ contains ! Set global variables if (params%iprof == 0) then ! Analytical profiles - call set_profiles_an(data) + call set_profiles_an(params, data) else ! Numerical profiles - call set_profiles_spline(params, data) + call set_profiles_spline(params, data, launch_pos) end if end subroutine init_profiles