experimenting with a python interface

work in progress

build_python.sh

@@ -0,0 +1,15 @@
+#!/bin/bash
+
+#
+# Build the python extension module using f2py
+#
+
+# this is to eliminate the "ld: warning: could not create compact unwind" on MacOS
+export LDFLAGS=-Wl,-no_compact_unwind
+
+f2py -c --backend meson ./python/radbeltpy.pyf \
+    ./src/radbelt_kinds_module.F90 \
+    ./src/trmfun.f90 \
+    ./src/shellig.f90 \
+    ./src/radbelt_module.f90 \
+    ./src/radbelt_c_module.f90

meson.build

@@ -0,0 +1,26 @@
+project('radbelt', ['c', 'fortran'],
+        version: '0.0.1',
+        default_options : [
+          'warning_level=everything',
+          'buildtype=release',]
+        )
+
+add_languages('fortran')
+
+_args = [] # Extra arguments
+_deps = [] # Dependencies
+# _deps += dependency('lapack')
+
+cc = meson.get_compiler('c')
+
+radbeltlib = library('RadBelt',
+              sources: [
+                'src/radbelt_kinds_module.F90' ,
+                'src/trmfun.f90' ,
+                'src/shellig.f90' ,
+                'src/radbelt_module.f90' ,
+                'src/radbelt_c_module.f90' ,
+              ],
+              dependencies: _deps,
+              install: true,
+)

python/.f2py_f2cmap

@@ -0,0 +1,13 @@
+dict(
+    real=dict(
+        real64="double",
+        real32="float",
+        c_double="double",
+        c_float="float"),
+    integer=dict(
+        int64="long long",
+        int32="int",
+        c_int="int",
+        c_long="long",
+        c_long_long="long long"),
+)

python/environment.yml

@@ -0,0 +1,20 @@
+# To use:
+#
+#   $ conda env create --prefix ./env -f ./python/environment.yml
+#   $ conda activate ./env
+
+channels:
+  - nodefaults
+  - conda-forge
+
+dependencies:
+  #- compilers
+  #- gfortran
+  - fortran-compiler
+  - meson
+  - pkg-config
+  - fpm
+  - ford
+  - graphviz
+  - fprettify
+  - numpy

python/radbeltpy.pyf

@@ -0,0 +1,27 @@
+!    -*- f90 -*-
+! Note: the context of this file is case sensitive.
+
+python module radbelt ! in
+    interface  ! in :radbelt
+        module radbelt_c_module ! in :radbelt:radbelt_c_module.f90
+            use iso_c_binding, only: c_double,c_int,c_char
+            subroutine get_flux_g_c(lon,lat,height,year,e,imname,flux)  ! in :radbelt:radbelt_c_module.f90:radbelt_c_module
+                real(c_double),intent(in) :: lon
+                real(c_double),intent(in) :: lat
+                real(c_double),intent(in) :: height
+                real(c_double),intent(in) :: year
+                real(c_double),intent(in) :: e
+                integer(c_int),intent(in) :: imname
+                real(c_double),intent(out) :: flux
+            end subroutine get_flux_g_c
+            subroutine set_data_files_paths_c(aep8_dir, igrf_dir)  ! in :radbelt:radbelt_c_module.f90:radbelt_c_module
+                character(kind=c_char,len=256),intent(in) :: aep8_dir
+                character(kind=c_char,len=256),intent(in) :: igrf_dir
+                !character(kind=c_char,len=1),dimension(*),intent(in) :: aep8_dir ! doesn't work
+                !character(kind=c_char,len=1),dimension(*),intent(in) :: igrf_dir
+            end subroutine set_data_files_paths_c
+        end module radbelt_c_module
+    end interface
+end python module radbelt
+
+! This file was manually created

src/radbelt_c_module.f90

@@ -0,0 +1,93 @@
+!*****************************************************************************************
+!>
+!  Experimental C interface to the radbelt module.
+!
+!@note Currently, this makes use of a global variable to store the state,
+!      so should not be considered threadsafe.
+
+module radbelt_c_module
+
+use iso_c_binding, only: c_double, c_int, c_char, c_null_char
+use radbelt_module, only: radbelt_type
+
+implicit none
+
+type(radbelt_type) :: radbelt !! global variable
+
+integer,parameter :: STR_LEN = 256 !! string length for paths
+
+contains
+
+! function c2f_str(cstr) result(fstr)
+!     character(kind=c_char,len=1),dimension(*),intent(in) :: cstr
+!     character(len=:),allocatable :: fstr
+!     integer :: i
+!     fstr = ''
+!     i = 0
+!     do
+!         i = i + 1
+!         if (cstr(i)==c_null_char) exit
+!         fstr = fstr // cstr(i)
+!     end do
+! end function c2f_str
+
+!*****************************************************************************************
+!>
+!  C interface for setting the data file paths
+
+subroutine set_data_files_paths_c(aep8_dir, igrf_dir) bind(C, name="set_data_files_paths_c")
+
+    character(kind=c_char,len=1),dimension(STR_LEN),intent(in) :: aep8_dir
+    character(kind=c_char,len=1),dimension(STR_LEN),intent(in) :: igrf_dir
+
+    character(len=:),allocatable :: aep8_dir_, igrf_dir_
+
+    c2f : block
+        ! convert to c strings to fortran
+        integer :: i
+        aep8_dir_ = ''
+        igrf_dir_ = ''
+        do i = 1, STR_LEN
+            aep8_dir_ = aep8_dir_//aep8_dir(i)
+            igrf_dir_ = igrf_dir_//igrf_dir(i)
+        end do
+        aep8_dir_ = trim(aep8_dir_)
+        igrf_dir_ = trim(igrf_dir_)
+    end block c2f
+
+    ! !... doesn't work this way... ??
+    ! aep8_dir_ = c2f_str(aep8_dir)
+    ! igrf_dir_ = c2f_str(igrf_dir)
+    ! write(*,*) trim(aep8_dir_)
+    ! write(*,*) trim(igrf_dir_)
+
+    call radbelt%set_data_files_paths(aep8_dir_, igrf_dir_)
+
+ end subroutine set_data_files_paths_c
+!*****************************************************************************************
+
+!*****************************************************************************************
+!>
+!  C interface to [[get_flux_g]].
+
+subroutine get_flux_g_c(lon,lat,height,year,e,imname,flux) bind(C, name="get_flux_g_c")
+
+    real(c_double),intent(in) :: lon !! geodetic longitude in degrees (east)
+    real(c_double),intent(in) :: lat !! geodetic latitude in degrees (north)
+    real(c_double),intent(in) :: height !! altitude in km above sea level
+    real(c_double),intent(in) :: year !! decimal year for which geomagnetic field is to
+                                      !! be calculated (e.g.:1995.5 for day 185 of 1995)
+    real(c_double),intent(in) :: e !! minimum energy
+    integer(c_int),intent(in) :: imname !! which method to use:
+                                 !!
+                                 !! * 1 -- particle species: electrons, solar activity: min
+                                 !! * 2 -- particle species: electrons, solar activity: max
+                                 !! * 3 -- particle species: protons, solar activity: min
+                                 !! * 4 -- particle species: protons, solar activity: max
+    real(c_double),intent(out) :: flux !! The flux of particles above the given energy, in units of cm^-2 s^-1.
+
+    flux = radbelt%get_flux(lon,lat,height,year,e,imname)
+
+end subroutine get_flux_g_c
+
+end module radbelt_c_module

src/shellig.f90

@@ -214,6 +214,7 @@ subroutine findb0(me,stps,bdel,value,bequ,rr0)
 
       step=stps
       irun=0
+      rold = 0.0_wp ! to avoid -Wmaybe-uninitialized warnings
 
       main : do
         irun=irun+1

src/trmfun.f90

@@ -152,6 +152,11 @@ subroutine trara1(me,descr,map,fl,bb0,e,f,n)
    integer :: i0 , i1 , i2 , i3 , ie , l3 , nb , nl
    logical :: s0 , s1 , s2
 
+   e0 = 0.0_wp  ! to avoid -Wmaybe-uninitialized warnings
+   f0 = 0.0_wp  ! to avoid -Wmaybe-uninitialized warnings
+   i0 = 0       ! to avoid -Wmaybe-uninitialized warnings
+   s0 = .false. ! to avoid -Wmaybe-uninitialized warnings  -- but not sure what default value here should be !  -JW
+
    bb0_ = bb0
    me%fistep = descr(7)/descr(2)
    escale = descr(4)

test/radbeltpy_test.py

@@ -0,0 +1,94 @@
+#
+# Test of the f2py-created Python interface.
+#
+
+from pathlib import Path
+
+dir = Path(__file__).resolve().parents[1] # root directory
+
+import sys
+sys.path.insert(0,'.')   # assuming the radbelt lib is in the main directory
+
+import radbelt
+import time
+import os
+# import numpy as np
+
+#########################################################################################################
+def set_data_files_paths(aep8_dir : str, igrf_dir : str) -> None:
+    """Python function to set the file paths"""
+    radbelt.radbelt_c_module.set_data_files_paths_c(aep8_dir, igrf_dir)
+
+#########################################################################################################
+#@np.vectorize   <-- makes it run slower for scalar inputs??
+def get_flux(lon : float, lat : float , height : float, year : float, e : float, imname :int) -> float:
+    """Python function to get the flux"""
+    return radbelt.radbelt_c_module.get_flux_g_c(lon,lat,height,year,e,imname)
+
+
+# set location of the data files:
+# dir = Path(__file__).resolve().parents[1] # root directory
+aep8_dir = str(dir / 'data' / 'aep8')
+igrf_dir = str(dir / 'data' / 'igrf')
+set_data_files_paths(aep8_dir, igrf_dir)
+
+EPS = sys.float_info.epsilon # machine precision for error checking
+lon = -45.0
+lat = -30.0
+height = 500.0
+year = 2021.1616438356164  # decimal year
+imname = 4 # 'p', 'max'
+e = 20.0
+
+flux = get_flux(lon,lat,height,year,e,imname)
+
+print(f'flux = {flux}')
+
+error = flux - 2642.50370051985726336559603128948869  #difference from real128 version
+relerror = abs(error/flux)
+
+print(f'Flux      = {flux}')
+print(f'Error     = {error}')
+print(f'Rel Error = {relerror}')
+if relerror>10*EPS:
+    raise Exception('error')
+
+print('')
+
+n_cases = 0
+tstart = time.time()
+for lat in range(-90, 91, 5):
+    for lon in range(-180, 181, 45):
+        for alt in range(500, 1001, 100):
+            n_cases = n_cases + 1
+            f = get_flux(lon,lat,height,year, e, 4)
+tend = time.time()
+print(f'Python version runtime: {tend-tstart} sec. {int(n_cases/(tend-tstart))} (cases/sec)')
+
+# print('... VECTORIZE ...')
+# n_cases = 0
+# tstart = time.time()
+
+# lon_vec = []
+# lat_vec = []
+# height_vec = []
+# year_vec = []
+# e_vec = []
+# ifile_vec = []
+
+# for lat in range(-90, 91, 5):
+#     for lon in range(-180, 181, 45):
+#         for alt in range(500, 1001, 100):
+#             n_cases = n_cases + 1
+#             lon_vec.append(lon)
+#             lat_vec.append(lat)
+#             height_vec.append(height)
+#             year_vec.append(year)
+#             e_vec.append(e)
+#             ifile_vec.append(4)
+#             #f = get_flux(lon,lat,height,year, e, 4)
+
+# f_vec = get_flux(lon_vec,lat_vec,height_vec,year_vec,e_vec,ifile_vec)
+# tend = time.time()
+# print(f'Python version runtime: {tend-tstart} sec. {int(n_cases/(tend-tstart))} (cases/sec)')
+