src/gray_plasma.f90: use mollifier to interpolate density
This methods ensures n_e(ψ) is everywhere positive and C²-continuous. The density boundary is now dependent on the width w of the mollifier, which also controls the smoothness of the curve, specifically: ψ_bnd = ψ_last + w where ψ_last = 1.01 currently is an extra point added to improve the convergence of the mollified density to the original data as w→0.
This commit is contained in:
Michele Guerini Rocco
parent
98a855bdaa
commit
3a2dd3697f
GPG Key ID:
BFBAF4C975F76450
@ -198,7 +198,7 @@ MHD equilibrium parameters
|
||||
- *0*, use sign from COCOS
|
||||
|
||||
**sgni** (default: **0**)
|
||||
: Force the sign of the plasma current. One of: *+1*,
|
||||
: Force the sign of the plasma current. One of:
|
||||
|
||||
- *+1*, counter-clockwise (viewed from above)
|
||||
- *-1*, clockwise
|
||||
@ -226,10 +226,11 @@ MHD equilibrium parameters
|
||||
**filenm**
|
||||
: Filepath of the plasma profiles data (relative to the gray.ini file)
|
||||
|
||||
**sspld** (default: **0.1**)
|
||||
: Tension of the electron density spline
|
||||
**sspld** (default: **0.020**)
|
||||
: Width of the the electron density mollifier (>0)
|
||||
|
||||
Note: 0 means perfect interpolation
|
||||
Large values produce smoother profiles but also increase
|
||||
the density boundary (`psnbnd = 1.01 + sspld`).
|
||||
|
||||
**factte** (default: **1.0**)
|
||||
: Rescaling factor for the electron temperature factte = 1
|
||||
@ -305,7 +306,6 @@ Below is an example of a gray.ini file
|
||||
iprof = PROF_ANALYTIC
|
||||
irho = RHO_TOR
|
||||
filenm = "profiles.txt"
|
||||
sspld = 0.01
|
||||
|
||||
[output]
|
||||
ipec = RHO_POL
|
||||
|
||||
@ -21,7 +21,7 @@ equiln.txt x filename for equilibrium
|
||||
|
||||
1 2 x iprof= 0/1 analytical/numerical profiles, irho=0/1/2 radial coordinate rhot/rhop/psi
|
||||
profilesn.txt x filename for profiles
|
||||
1.05 0.005 : psi plasma boundary, sspld spline coeff density
|
||||
1.05 0.005 : psi plasma boundary, sspld width of density mollifier
|
||||
1 1 1 x factT, factn, iscal=0/1/2 nustar=const/ngreenw=const/no_scal.
|
||||
|
||||
1 501 x ipec=±1/2 profiles rhop/rhot; nnd=no. radial points (no. of radial intervals +1)
|
||||
|
||||
@ -123,8 +123,8 @@ module gray_params
|
||||
|
||||
! Kinetic plasma profiles
|
||||
type profiles_parameters
|
||||
real(wp_) :: psnbnd = 1.01_wp_ ! Plasma boundary (ψ: ne(ψ)=0)
|
||||
real(wp_) :: sspld = 0.1_wp_ ! Density spline tension
|
||||
real(wp_) :: psnbnd = 1.0_wp_ ! Plasma boundary (ψ: ne(ψ)=0)
|
||||
real(wp_) :: sspld = 0.020_wp_ ! Density mollifier width
|
||||
real(wp_) :: factne = 1.0_wp_ ! Rescale factor for n_e
|
||||
real(wp_) :: factte = 1.0_wp_ ! Rescale factor for T_e
|
||||
integer(kind(scaling_enum)) :: iscal = SCALE_OFF ! Rescaling model for n_e,T_e
|
||||
|
||||
@ -8,7 +8,7 @@
|
||||
!
|
||||
module gray_plasma
|
||||
use const_and_precisions, only : wp_, zero
|
||||
use splines, only : spline_1d
|
||||
use splines, only : spline_1d, mollifier_1d
|
||||
|
||||
implicit none
|
||||
|
||||
@ -67,22 +67,11 @@ module gray_plasma
|
||||
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_1d) :: temp_spline
|
||||
type(spline_1d) :: zeff_spline
|
||||
type(density_tail) :: tail
|
||||
type(mollifier_1d) :: dens_spline
|
||||
type(spline_1d) :: temp_spline, zeff_spline
|
||||
contains
|
||||
procedure :: init => init_splines
|
||||
procedure :: density => numeric_density
|
||||
@ -173,59 +162,13 @@ contains
|
||||
! 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
|
||||
! Past ψ=1.01 + the width of the mollifier, both the
|
||||
! density and its derivative are identically zero
|
||||
if (psin >= 1.01_wp_ + self%dens_spline%w) 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 = 0
|
||||
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
|
||||
! Interpolate the data
|
||||
dens = self%dens_spline%eval(psin)
|
||||
ddens = self%dens_spline%eval(psin, order=1)
|
||||
|
||||
end subroutine numeric_density
|
||||
|
||||
@ -410,77 +353,14 @@ contains
|
||||
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)
|
||||
! Mollified linear interpolation of density
|
||||
! Note: we add an extra datum to ensure n_e=0
|
||||
call self%dens_spline%init( &
|
||||
x=[psi, 1.01_wp_], &
|
||||
y=[dens, 0.00_wp_], &
|
||||
width=params%profiles%sspld)
|
||||
|
||||
! 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
|
||||
! Make sure the wave launcher does not fall inside the plasma
|
||||
! Note: if it does, the initial wave conditions become
|
||||
! invalid as they are given assuming a vacuum (N=1)
|
||||
block
|
||||
@ -497,32 +377,34 @@ contains
|
||||
! Note: this returns -1 when the data is not available
|
||||
call equil%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
|
||||
if (psi > 0 .and. psi < 1) 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
|
||||
'launcher: ψ=', psi, ' boundary: ψ=1'
|
||||
call log_error(msg, mod='coreprofiles', proc='init_splines')
|
||||
err = 2
|
||||
return
|
||||
end if
|
||||
|
||||
if (psi > 1 .and. psi < 1.01_wp_ + self%dens_spline%w) then
|
||||
! Shrink the density mollifier width
|
||||
call log_warning('reduced density mollification to not' // &
|
||||
' reach the wave launcher', &
|
||||
mod='coreprofiles', proc='init_splines')
|
||||
write (msg, '(a, g0.3, a, g0.3, " → ", g0.3)') &
|
||||
'launcher: ψ=', psi, ' boundary: ψ=', &
|
||||
1.01_wp_+self%dens_spline%w, (1.01_wp_+psi)/2
|
||||
call log_debug(msg, mod='coreprofiles', proc='init_splines')
|
||||
self%dens_spline%w = (psi - 1.01_wp_)/2
|
||||
end if
|
||||
end block
|
||||
|
||||
! Set the density boundary ψ
|
||||
! Note: this is used to detect entrance in the plasma
|
||||
params%profiles%psnbnd = self%tail%end
|
||||
params%profiles%psnbnd = 1.01_wp_ + self%dens_spline%w
|
||||
|
||||
write (msg, '(a,g0.4)') 'density boundary: ψ=', self%tail%end
|
||||
write (msg, '(a,g0.4)') 'density boundary: ψ=', params%profiles%psnbnd
|
||||
call log_info(msg, mod='coreprofiles', proc='init_splines')
|
||||
end subroutine init_splines
|
||||
|
||||
|
||||
@ -11,7 +11,7 @@ contains
|
||||
|
||||
pure subroutine locate(array, value, index)
|
||||
! Given an `array`, and a `value` returns the `index`
|
||||
! such that `array(index) < value < array(index+1)`
|
||||
! such that `array(index) < value ≤ array(index+1)`
|
||||
!
|
||||
! Notes:
|
||||
! 1. `array` must be monotonic, either increasing or decreasing.
|
||||
|
||||
Loading…
Reference in New Issue
Block a user