Skip to content

Replace GSL Bessel I_n with fortnum (1/5) #29

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Draft
krystophny wants to merge 4 commits into
base: main
Choose a base branch
from
Draft

Replace GSL Bessel I_n with fortnum (1/5) #29

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
4 commits
Select commit Hold shift + click to select a range
dae6879
Replace GSL Bessel I_n with fortnum bessel_in_array
krystophny Jun 14, 2026
8f3e113
Fetch fortnum over https so CI can clone the public repo
krystophny Jun 14, 2026
54232cd
Pin fortnum to the v0.1.0 release tag
krystophny Jun 14, 2026
6d760e6
Re-pin fortnum to current main after upstream history rewrite
krystophny Jun 15, 2026
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
14 changes: 13 additions & 1 deletion CMakeLists.txt
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -100,6 +100,17 @@ endif()
include(Util)
find_or_fetch(libneo)

### fortnum: numerical core replacing the former GSL routines. libneo may pull
### fortnum in transitively, so guard against declaring the target twice.
include(FetchContent)
if(NOT TARGET fortnum)
FetchContent_Declare(fortnum
GIT_REPOSITORY https://github.com/lazy-fortran/fortnum.git
GIT_TAG 92de6e949a772cfffc73bb5295fe5e2b056b9c18
)
FetchContent_MakeAvailable(fortnum)
endif()

### Project source files
# SuiteSparse
set(SUITESPARSE_SRC_FILES
Expand Down Expand Up @@ -152,7 +163,7 @@ if(WITH_MFEM)
)
endif()
set_source_files_properties(src/mephit_run.c ${MEPHIT_C_SRC_FILES} ${MEPHIT_CPP_SRC_FILES}
PROPERTIES COMPILE_FLAGS "${CMAKE_C_FLAGS} -DREAL=double -I${TRIANGLE_INCLUDE_DIR} -I${SUITESPARSE_INCLUDE_DIRS} -L${TRIANGLE_LIB_DIR}")
PROPERTIES COMPILE_FLAGS "${CMAKE_C_FLAGS} -DREAL=double -I${TRIANGLE_INCLUDE_DIR} -I${SUITESPARSE_INCLUDE_DIRS} -I${fortnum_SOURCE_DIR}/include -L${TRIANGLE_LIB_DIR}")


### Define library
Expand All @@ -176,6 +187,7 @@ set(MEPHIT_LIBS
hdf5::hdf5_fortran
hdf5::hdf5_hl_fortran
GSL::gsl
fortnum
BLAS::BLAS
LAPACK::LAPACK
${FFTW_LIBRARIES}
Expand Down
28 changes: 8 additions & 20 deletions src/mephit_pert.f90
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1,7 +1,6 @@
module mephit_pert

use iso_fortran_env, only: dp => real64
use iso_c_binding, only: c_int, c_double

implicit none

Expand Down Expand Up @@ -37,17 +36,6 @@ function vector_element_projection(ktri, weight, R, Z, n, f)
end function vector_element_projection
end interface

interface
function gsl_sf_bessel_icn_array(nmin, nmax, x, result_array) &
bind(C, name='gsl_sf_bessel_In_array')
import :: c_int, c_double
integer(c_int), intent(in), value :: nmin, nmax
real(c_double), intent(in), value :: x
real(c_double), intent(inout), dimension(nmax - nmin + 1) :: result_array
integer(c_int) :: gsl_sf_bessel_icn_array
end function gsl_sf_bessel_icn_array
end interface

type :: L1_t
!> Number of points on which the L1 are defined
integer :: npoint
Expand Down Expand Up @@ -1427,21 +1415,21 @@ end subroutine compute_kilca_vac_coeff
!> field
subroutine kilca_vacuum_fourier(tor_mode, pol_mode, R_0, r, vac_coeff, &
B_rad, B_pol, B_tor)
use mephit_conf, only: logger
use mephit_util, only: imun
use fortnum_special, only: bessel_in_array
integer, intent(in) :: tor_mode, pol_mode
real(dp), intent(in) :: R_0, r
complex(dp), intent(in) :: vac_coeff
complex(dp), intent(out) :: B_rad, B_pol, B_tor
real(c_double) :: I_m(-1:1), k_z_r
integer(c_int) :: status
real(dp) :: I_n(0:abs(pol_mode) + 1), I_m(-1:1), k_z_r

k_z_r = tor_mode / R_0 * r
status = gsl_sf_bessel_icn_array(abs(pol_mode)-1, abs(pol_mode)+1, k_z_r, I_m)
if (status /= 0 .and. logger%err) then
write (logger%msg, '("gsl_sf_bessel_In_array returned error ", i0)') status
call logger%write_msg
end if
! fortnum fills I_n(0:nmax) with I_0..I_nmax; recover the order |pol_mode|-1
! (which is -1 for pol_mode = 0) via the symmetry I_{-n}(x) = I_n(x).
call bessel_in_array(abs(pol_mode) + 1, k_z_r, I_n)
I_m(-1) = I_n(abs(abs(pol_mode) - 1))
I_m(0) = I_n(abs(pol_mode))
I_m(1) = I_n(abs(pol_mode) + 1)
B_rad = 0.5d0 * (I_m(-1) + I_m(1)) * vac_coeff
B_pol = imun * pol_mode / k_z_r * I_m(0) * vac_coeff
B_tor = imun * I_m(0) * vac_coeff
Expand Down
Loading