From c55e866ba651c79e9da8e915bd2a606c3fc5b647 Mon Sep 17 00:00:00 2001 From: Michele Guerini Rocco Date: Thu, 13 Feb 2025 09:54:13 +0100 Subject: [PATCH] src/gray_plasma.f90: use mollifier to interpolate density MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit 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. --- doc/man/gray.ini.5.md | 10 +- doc/man/gray_params.data.5.md | 2 +- src/gray_params.f90 | 4 +- src/gray_plasma.f90 | 182 ++++++---------------------------- src/utils.f90 | 2 +- 5 files changed, 41 insertions(+), 159 deletions(-) diff --git a/doc/man/gray.ini.5.md b/doc/man/gray.ini.5.md index cd135cc..50ac1c5 100644 --- a/doc/man/gray.ini.5.md +++ b/doc/man/gray.ini.5.md @@ -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 diff --git a/doc/man/gray_params.data.5.md b/doc/man/gray_params.data.5.md index ce1380a..bf98803 100644 --- a/doc/man/gray_params.data.5.md +++ b/doc/man/gray_params.data.5.md @@ -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) diff --git a/src/gray_params.f90 b/src/gray_params.f90 index d76dd94..faf0cb4 100644 --- a/src/gray_params.f90 +++ b/src/gray_params.f90 @@ -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 diff --git a/src/gray_plasma.f90 b/src/gray_plasma.f90 index 1392f5a..f736b58 100644 --- a/src/gray_plasma.f90 +++ b/src/gray_plasma.f90 @@ -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 diff --git a/src/utils.f90 b/src/utils.f90 index bd1d05c..c57b07e 100644 --- a/src/utils.f90 +++ b/src/utils.f90 @@ -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.