diff --git a/build_python.sh b/build_python.sh new file mode 100755 index 0000000..336a614 --- /dev/null +++ b/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 \ No newline at end of file diff --git a/meson.build b/meson.build new file mode 100644 index 0000000..faeab68 --- /dev/null +++ b/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, +) diff --git a/python/.f2py_f2cmap b/python/.f2py_f2cmap new file mode 100644 index 0000000..4f5336c --- /dev/null +++ b/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"), +) diff --git a/python/environment.yml b/python/environment.yml new file mode 100644 index 0000000..e7990ee --- /dev/null +++ b/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 diff --git a/python/radbeltpy.pyf b/python/radbeltpy.pyf new file mode 100644 index 0000000..6967b34 --- /dev/null +++ b/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 diff --git a/src/radbelt_c_module.f90 b/src/radbelt_c_module.f90 new file mode 100644 index 0000000..5d7d754 --- /dev/null +++ b/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 diff --git a/src/shellig.f90 b/src/shellig.f90 index 20846ba..68c65ce 100644 --- a/src/shellig.f90 +++ b/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 diff --git a/src/trmfun.f90 b/src/trmfun.f90 index 3c784e9..5565c31 100644 --- a/src/trmfun.f90 +++ b/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) diff --git a/test/radbeltpy_test.py b/test/radbeltpy_test.py new file mode 100644 index 0000000..bfc940b --- /dev/null +++ b/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)') +