From 97906348ba3b17721378f1c1f2284ffaf7a1dfff Mon Sep 17 00:00:00 2001
From: Michele Guerini Rocco <m.guerinirocco2@campus.unimib.it>
Date: Thu, 13 Feb 2025 09:49:25 +0100
Subject: [PATCH] src/splines.f90: remove spline_simple
---
src/beams.f90 | 136 +++++++++++++++++---------
src/gray_equil.f90 | 34 ++++---
src/gray_plasma.f90 | 10 +-
src/splines.f90 | 228 ++++++++------------------------------------
4 files changed, 150 insertions(+), 258 deletions(-)
diff --git a/src/beams.f90 b/src/beams.f90
index 3210ee4..efc924d 100644
--- a/src/beams.f90
+++ b/src/beams.f90
@@ -95,7 +95,7 @@ contains
! ellipses in the transverse plane at the launch point (deg)
use gray_params, only : antenna_parameters
- use splines, only : spline_simple
+ use splines, only : spline_1d, linear_1d
use utils, only : locate
use logger, only : log_error
@@ -105,13 +105,11 @@ contains
! local variables
integer :: u, nisteer, i, k, ii
- real(wp_) :: steer, dal
- real(wp_), dimension(:), allocatable :: &
- alphastv, betastv, x00v, y00v, &
+ character(256) :: msg
+ real(wp_) :: steer
+ real(wp_), dimension(:), allocatable :: &
+ alphastv, betastv, x00v, y00v, &
z00v, waist1v, waist2v, rci1v, rci2v, phi1v, phi2v
- type(spline_simple) :: beta, waist1, waist2, &
- rci1, rci2, phi1, phi2, &
- x0, y0, z0
open(newunit=u, file=params%filenm, status='old', action='read', iostat=err)
if (err /= 0) then
@@ -145,47 +143,95 @@ contains
rci1v = 10 * rci1v
rci2v = 10 * rci2v
- call beta%init(alphastv, betastv, nisteer)
- call waist1%init(alphastv, waist1v, nisteer)
- call waist2%init(alphastv, waist2v, nisteer)
- call rci1%init(alphastv, rci1v, nisteer)
- call rci2%init(alphastv, rci2v, nisteer)
- call phi1%init(alphastv, phi1v, nisteer)
- call phi2%init(alphastv, phi2v, nisteer)
- call x0%init(alphastv, x00v, nisteer)
- call y0%init(alphastv, y00v, nisteer)
- call z0%init(alphastv, z00v, nisteer)
+ call locate(alphastv, params%alpha, k)
- if((params%alpha > alphastv(1)) .and. (params%alpha < alphastv(nisteer))) then
- call locate(alphastv, params%alpha, k)
- dal = params%alpha - alphastv(k)
- params%beta = beta%raw_eval(k, dal)
- params%pos(1) = x0%raw_eval(k, dal)
- params%pos(2) = y0%raw_eval(k, dal)
- params%pos(3) = z0%raw_eval(k, dal)
- params%w(1) = waist1%raw_eval(k, dal)
- params%w(2) = waist2%raw_eval(k, dal)
- params%ri(1) = rci1%raw_eval(k, dal)
- params%ri(2) = rci2%raw_eval(k, dal)
- params%phi(1) = phi1%raw_eval(k, dal)
- params%phi(2) = phi2%raw_eval(k, dal)
- else
- ! params%alpha outside table range
- if(params%alpha >= alphastv(nisteer)) ii=nisteer
- if(params%alpha <= alphastv(1)) ii=1
- params%alpha = alphastv(ii)
- params%beta = betastv(ii)
- params%pos(1) = x00v(ii)
- params%pos(2) = y00v(ii)
- params%pos(3) = z00v(ii)
- params%w(1) = waist1v(ii)
- params%w(2) = waist2v(ii)
- params%ri(1) = rci1v(ii)
- params%ri(2) = rci2v(ii)
- params%phi(1) = phi1v(ii)
- params%phi(2) = phi2v(ii)
+ ! No interpolation: 0D table or out of range
+ if (nisteer == 1 .or. k == 0 .or. k == nisteer) then
+ if (k == 0) k = k + 1
+ params%alpha = alphastv(k)
+ params%beta = betastv(k)
+ params%pos(1) = x00v(k)
+ params%pos(2) = y00v(k)
+ params%pos(3) = z00v(k)
+ params%w(1) = waist1v(k)
+ params%w(2) = waist2v(k)
+ params%ri(1) = rci1v(k)
+ params%ri(2) = rci2v(k)
+ params%phi(1) = phi1v(k)
+ params%phi(2) = phi2v(k)
+ return
end if
+ if (nisteer < 4) then
+ ! Linear interpolation for small tables
+ block
+ type(linear_1d) :: beta, waist1, waist2, &
+ rci1, rci2, phi1, phi2, &
+ x0, y0, z0
+
+ call beta%init(alphastv, betastv, nisteer)
+ call waist1%init(alphastv, waist1v, nisteer)
+ call waist2%init(alphastv, waist2v, nisteer)
+ call rci1%init(alphastv, rci1v, nisteer)
+ call rci2%init(alphastv, rci2v, nisteer)
+ call phi1%init(alphastv, phi1v, nisteer)
+ call phi2%init(alphastv, phi2v, nisteer)
+ call x0%init(alphastv, x00v, nisteer)
+ call y0%init(alphastv, y00v, nisteer)
+ call z0%init(alphastv, z00v, nisteer)
+
+ params%beta = beta%raw_eval(k, params%alpha)
+ params%pos(1) = x0%raw_eval(k, params%alpha)
+ params%pos(2) = y0%raw_eval(k, params%alpha)
+ params%pos(3) = z0%raw_eval(k, params%alpha)
+ params%w(1) = waist1%raw_eval(k, params%alpha)
+ params%w(2) = waist2%raw_eval(k, params%alpha)
+ params%ri(1) = rci1%raw_eval(k, params%alpha)
+ params%ri(2) = rci2%raw_eval(k, params%alpha)
+ params%phi(1) = phi1%raw_eval(k, params%alpha)
+ params%phi(2) = phi2%raw_eval(k, params%alpha)
+ return
+ end block
+ end if
+
+ ! Cubic interpolation otherwise
+ block
+ type(spline_1d) :: beta, waist1, waist2, &
+ rci1, rci2, phi1, phi2, &
+ x0, y0, z0
+
+ call beta%init(alphastv, betastv, nisteer, err=err)
+ call waist1%init(alphastv, waist1v, nisteer, err=err)
+ call waist2%init(alphastv, waist2v, nisteer, err=err)
+ call rci1%init(alphastv, rci1v, nisteer, err=err)
+ call rci2%init(alphastv, rci2v, nisteer, err=err)
+ call phi1%init(alphastv, phi1v, nisteer, err=err)
+ call phi2%init(alphastv, phi2v, nisteer, err=err)
+ call x0%init(alphastv, x00v, nisteer, err=err)
+ call y0%init(alphastv, y00v, nisteer, err=err)
+ call z0%init(alphastv, z00v, nisteer, err=err)
+
+ if (err > 0) then
+ write (msg, '(a, a, g0)') 'interpolation failed', &
+ '`curfit` returned ier=', err
+ call log_error(msg, mod='beams', proc='read_beam1')
+ return
+ end if
+ err = 0
+
+ call locate(beta%knots, params%alpha, k)
+ params%beta = beta%raw_eval(k, params%alpha)
+ params%pos(1) = x0%raw_eval(k, params%alpha)
+ params%pos(2) = y0%raw_eval(k, params%alpha)
+ params%pos(3) = z0%raw_eval(k, params%alpha)
+ params%w(1) = waist1%raw_eval(k, params%alpha)
+ params%w(2) = waist2%raw_eval(k, params%alpha)
+ params%ri(1) = rci1%raw_eval(k, params%alpha)
+ params%ri(2) = rci2%raw_eval(k, params%alpha)
+ params%phi(1) = phi1%raw_eval(k, params%alpha)
+ params%phi(2) = phi2%raw_eval(k, params%alpha)
+ end block
+
end subroutine read_beam1
diff --git a/src/gray_equil.f90 b/src/gray_equil.f90
index 38c413c..8e34295 100644
--- a/src/gray_equil.f90
+++ b/src/gray_equil.f90
@@ -4,7 +4,7 @@
!
module gray_equil
use const_and_precisions, only : wp_, comp_huge
- use splines, only : spline_simple, spline_1d, spline_2d, linear_1d
+ use splines, only : spline_1d, spline_2d, linear_1d
use types, only : contour
implicit none
@@ -26,12 +26,12 @@ module gray_equil
real(wp_) :: z_boundary(2) ! z range of the plasma boundary
! Flux average splines (see `flux_average`)
- type(spline_simple) :: spline_area, spline_volume
- type(spline_simple) :: spline_dAdrho_t, spline_dVdrho_t
- type(spline_simple) :: spline_R_J, spline_B_avg, spline_B_min, spline_B_max
- type(spline_simple) :: spline_f_c, spline_q, spline_I_p, spline_J_phi_avg
- type(spline_simple) :: spline_astra_tor, spline_jintrac_tor, spline_paral_tor
- type(spline_2d) :: spline_cd_eff
+ type(spline_1d) :: spline_area, spline_volume
+ type(spline_1d) :: spline_dAdrho_t, spline_dVdrho_t
+ type(spline_1d) :: spline_R_J, spline_B_avg, spline_B_min, spline_B_max
+ type(spline_1d) :: spline_f_c, spline_q, spline_I_p, spline_J_phi_avg
+ type(spline_1d) :: spline_astra_tor, spline_jintrac_tor, spline_paral_tor
+ type(spline_2d) :: spline_cd_eff
contains
procedure(pol_flux_sub), deferred :: pol_flux
@@ -147,11 +147,10 @@ module gray_equil
private
real(wp_) :: fpol_a ! Poloidal current at the edge (r=a), F_a
! Splines
- type(spline_2d) :: psi_spline
- type(contour) :: psi_domain
- type(spline_1d) :: fpol_spline
- type(spline_simple) :: q_spline
- type(linear_1d) :: rhop_spline, rhot_spline
+ type(spline_2d) :: psi_spline
+ type(contour) :: psi_domain
+ type(spline_1d) :: fpol_spline, q_spline
+ type(linear_1d) :: rhop_spline, rhot_spline
contains
procedure :: pol_flux => numeric_pol_flux
procedure :: pol_curr => numeric_pol_curr
@@ -773,6 +772,7 @@ contains
use gray_params, only : EQ_EQDSK_FULL, EQ_EQDSK_PARTIAL
use gray_params, only : X_AT_TOP, X_AT_BOTTOM, X_IS_MISSING
use utils, only : vmaxmin
+ use numint, only : simpson
use logger, only : log_debug, log_info
! subroutine arguments
@@ -1017,13 +1017,11 @@ contains
call self%q_spline%init(data%psi, data%q, npsi)
- ! Toroidal flux as Φ(ψ) = 2π ∫ q(ψ)dψ
- ! Specifically, Φ(ψ_n) = 2π|ψ_a| ∫₀ψ_n q(ψ_n)dψ_n, with ψ_n ∈ [0,1]
+ ! Toroidal flux as Φ(ψ) = 2π ∫₀^ψ q(ψ) dψ
phi(1) = 0
- do k = 1, npsi-1
- dx = self%q_spline%data(k+1) - self%q_spline%data(k)
- phi(k+1) = phi(k) + dx*(self%q_spline%coeffs(k,1) + dx*(self%q_spline%coeffs(k,2)/2 + &
- dx*(self%q_spline%coeffs(k,3)/3 + dx* self%q_spline%coeffs(k,4)/4) ) )
+ do k = 2, npsi
+ call simpson(2, h=data%psi(2)-data%psi(1), fi=data%q(k-1:k), s=phi(k))
+ phi(k) = phi(k) + phi(k-1)
end do
self%phi_a = phi(npsi)
diff --git a/src/gray_plasma.f90 b/src/gray_plasma.f90
index 30a198e..1392f5a 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_simple, spline_1d
+ use splines, only : spline_1d
implicit none
@@ -79,10 +79,10 @@ module gray_plasma
! 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
+ type(spline_1d) :: dens_spline
+ type(spline_1d) :: temp_spline
+ type(spline_1d) :: zeff_spline
+ type(density_tail) :: tail
contains
procedure :: init => init_splines
procedure :: density => numeric_density
diff --git a/src/splines.f90 b/src/splines.f90
index dc2d599..14f25a6 100644
--- a/src/splines.f90
+++ b/src/splines.f90
@@ -1,7 +1,6 @@
! This module provides a high-level interface for creating and evaluating
! several kind of splines:
!
-! `spline_simple` is a simple interpolating cubic spline,
! `spline_1d` and `spline_2d` are wrapper around the DIERCKX cubic splines,
! `linear_1d` is a linear interpolation.
module splines
@@ -9,27 +8,16 @@ module splines
implicit none
- ! A 1D interpolating cubic spline
- type spline_simple
- integer :: ndata ! Number of data points
- real(wp_), allocatable :: data(:) ! Data points (ndata)
- real(wp_), allocatable :: coeffs(:,:) ! Spline coefficients (ndata, 4)
- contains
- procedure :: init => spline_simple_init
- procedure :: eval => spline_simple_eval
- procedure :: raw_eval => spline_simple_raw_eval
- procedure :: deriv => spline_simple_deriv
- end type
-
! A 1D smoothing/interpolating cubic spline
type spline_1d
integer :: nknots ! Number of spline knots
real(wp_), allocatable :: knots(:) ! Knots positions
real(wp_), allocatable :: coeffs(:) ! B-spline coefficients
contains
- procedure :: init => spline_1d_init
- procedure :: eval => spline_1d_eval
- procedure :: deriv => spline_1d_deriv
+ procedure :: init => spline_1d_init
+ procedure :: eval => spline_1d_eval
+ procedure :: raw_eval => spline_1d_raw_eval
+ procedure :: deriv => spline_1d_deriv
end type
! Wrapper to store pointers in an array
@@ -70,178 +58,10 @@ module splines
end type
private
- public spline_simple, spline_1d, spline_2d, linear_1d
+ public spline_1d, spline_2d, linear_1d
contains
- subroutine spline_simple_init(self, x, y, n)
- ! Initialises the spline
-
- ! subroutine arguments
- class(spline_simple), intent(inout) :: self
- integer, intent(in) :: n
- real(wp_), dimension(n), intent(in) :: x, y
-
- self%ndata = n
- allocate(self%data(n))
- allocate(self%coeffs(n, 4))
-
- self%data = x
- call spline_simple_coeffs(x, y, n, self%coeffs)
- end subroutine spline_simple_init
-
-
- pure function spline_simple_eval(self, x) result(y)
- ! Evaluates the spline at x
- use utils, only : locate
-
- ! subroutine arguments
- class(spline_simple), intent(in) :: self
- real(wp_), intent(in) :: x
- real(wp_) :: y
-
- ! local variables
- integer :: i
- real(wp_) :: dx
-
- call locate(self%data, x, i)
- i = min(max(1, i), self%ndata - 1)
- dx = x - self%data(i)
- y = self%raw_eval(i, dx)
- end function spline_simple_eval
-
-
- pure function spline_simple_raw_eval(self, i, dx) result(y)
- ! Evaluates the i-th polynomial of the spline at dx
-
- ! subroutine arguments
- class(spline_simple), intent(in) :: self
- integer, intent(in) :: i
- real(wp_), intent(in) :: dx
- real(wp_) :: y
-
- y = self%coeffs(i,1) + dx*(self%coeffs(i,2) &
- + dx*(self%coeffs(i,3) + dx*self%coeffs(i,4)))
- end function spline_simple_raw_eval
-
-
- pure function spline_simple_deriv(self, x) result(y)
- ! Computes the derivative of the spline at x
- use utils, only : locate
-
-
- ! subroutine arguments
- class(spline_simple), intent(in) :: self
- real(wp_), intent(in) :: x
- real(wp_) :: y
-
- ! local variables
- integer :: i
- real(wp_) :: dx
-
- call locate(self%data, x, i)
- i = min(max(1, i), self%ndata - 1)
- dx = x - self%data(i)
- y = self%coeffs(i,2) + dx*(2*self%coeffs(i,3) + 3*dx*self%coeffs(i,4))
- end function spline_simple_deriv
-
-
- pure subroutine spline_simple_coeffs(x, y, n, c)
- ! Computes the cubic coefficients of all n polynomials
-
- ! subroutine arguments
- integer, intent(in) :: n
- real(wp_), intent(in) :: x(n), y(n)
- real(wp_), intent(inout) :: c(n*4)
-
- ! local variables
- integer :: jmp, i, ii, j, j1, j2, j3, k
- real(wp_) :: xb, xc, ya, yb, h, a, r, dya, dyb, dy2
- jmp = 1
-
- if (n <= 1) return
-
- ! initialisation
- xc = x(1)
- yb = y(1)
- h = 0
- a = 0
- r = 0
- dyb = 0
- dy2 = c(2)
-
- ! form the system of linear equations
- ! and eliminate subsequentially
- j = 1
- do i = 1, n
- j2 = n + i
- j3 = j2 + n
- a = h*(2 - a)
- dya = dyb + h*r
- if (i >= n) then
- ! set derivative dy2 at last point
- dyb = dy2
- h = 0
- dyb = dya
- goto 13
- else
- j = j+jmp
- xb = xc
- xc = x(j)
- h = xc-xb
-
- ! ii=0 - increasing abscissae
- ! ii=1 - decreasing abscissae
- ii = 0
- if (h == 0) return
- if (h < 0) ii = 1
- ya = yb
- yb = y(j)
- dyb = (yb - ya)/h
-
- if (i <= 1) then
- j1 = ii
- goto 13
- end if
- end if
-
- if (j1-ii /= 0) return
- a = 1 / (2*h + a)
- 13 continue
- r = a*(dyb - dya)
- c(j3) = r
- a = h*a
- c(j2) = a
- c(i) = dyb
- end do
-
- ! back substitution of the system of linear equations
- ! and computation of the other coefficients
- a = 1
- j1 = j3+n+ii-ii*n
- i = n
- do k = 1, n
- xb = x(j)
- h = xc - xb
- xc = xb
- a = a+h
- yb = r
- r = c(j3)-r*c(j2)
- ya = 2*r
- c(j3) = ya + r
- c(j2) = c(i) - h*(ya+yb)
- c(j1) = (yb - r)/a
- c(i) = y(j)
- a = 0
- j = j-jmp
- i = i-1
- j2 = j2-1
- j3 = j3-1
- j1 = j3+n+ii
- end do
- end subroutine spline_simple_coeffs
-
-
subroutine spline_1d_init(self, x, y, n, range, weights, tension, err)
! Initialises a spline_1d.
! Takes:
@@ -260,7 +80,7 @@ contains
real(wp_), intent(in) :: x(n)
real(wp_), intent(in) :: y(n)
integer, intent(in) :: n
- real(wp_), intent(in) :: range(2)
+ real(wp_), intent(in), optional :: range(2)
real(wp_), intent(in), optional :: weights(n)
real(wp_), intent(in), optional :: tension
integer, intent(out), optional :: err
@@ -273,10 +93,12 @@ contains
! default values
integer :: err_def
- real(wp_) :: weights_def(n), tension_def
+ real(wp_) :: range_def(2), weights_def(n), tension_def
+ range_def = [x(1), x(n)]
weights_def = 1
tension_def = 0
+ if (present(range)) range_def = range
if (present(weights)) weights_def = weights
if (present(tension)) tension_def = tension
@@ -286,9 +108,10 @@ contains
allocate(self%knots(nknots_est), self%coeffs(nknots_est))
end if
- call curfit(0, n, x, y, weights_def, range(1), range(2), 3, tension_def, &
- nknots_est, self%nknots, self%knots, self%coeffs, residuals, &
- work_real, size(work_real), work_int, err_def)
+ call curfit(0, n, x, y, weights_def, range_def(1), range_def(2), &
+ 3, tension_def, nknots_est, self%nknots, self%knots, &
+ self%coeffs, residuals, work_real, size(work_real), &
+ work_int, err_def)
if (present(err)) err = err_def
end subroutine spline_1d_init
@@ -312,6 +135,31 @@ contains
end function spline_1d_eval
+ pure function spline_1d_raw_eval(self, i, x) result(y)
+ ! Evaluates the spline at: t(i) ≤ x < t(i+1), where t are the knots
+ use dierckx, only : fpbspl
+
+ ! subroutine arguments
+ class(spline_1d), intent(in) :: self
+ integer, intent(in) :: i
+ real(wp_), intent(in) :: x
+ real(wp_) :: y
+
+ ! local variables
+ integer :: err, j
+ real(wp_) :: h(6)
+
+ ! Evaluate the basis splines
+ call fpbspl(self%knots, self%nknots, 3, x, i, h)
+
+ ! Take the linear combination with the coefficients
+ y = 0
+ do j = 1, 4
+ y = y + self%coeffs(i-4+j)*h(j)
+ end do
+ end function spline_1d_raw_eval
+
+
pure function spline_1d_deriv(self, x, order) result(y)
! Evaluates the spline n-th order derivative at x
use dierckx, only : splder