pulling changes from master

.github/workflows/build.yml

@@ -8,7 +8,6 @@ jobs:
       fail-fast: false
     env:
       SPEC_PATH: ${{ github.workspace }}
-      PYTHONPATH: ${{ github.workspace }}/Utilities/pythontools
       OMPI_ALLOW_RUN_AS_ROOT: 1
       OMPI_ALLOW_RUN_AS_ROOT_CONFIRM: 1
     steps:
@@ -22,6 +21,8 @@ jobs:
         sudo apt-get update
         sudo apt-get install gfortran mpi-default-bin mpi-default-dev libhdf5-103 libhdf5-dev libfftw3-bin libfftw3-dev libopenblas0-openmp libopenblas-dev 
         pip3 install --user numpy f90nml scikit-build scipy h5py matplotlib
+        cd ${{ github.workspace }}/Utilities/pythontools
+        pip3 install -v .
     - name: compile_xspec
       run: |
         cd ${SPEC_PATH}
@@ -33,14 +34,12 @@ jobs:
     - name: run_fast_cartesian
       run: |
         cd ${SPEC_PATH}/ci/G1V03L2Fi
-        echo ${PYTHONPATH}
         export OMP_NUM_THREADS=1
         mpiexec -n 2 --allow-run-as-root ${SPEC_PATH}/xspec G1V03L2Fi.001.sp
         python3 -m py_spec.ci.test compare.h5 G1V03L2Fi.001.sp.h5 
     - name: run_fast_cylinder
       run: |
         cd ${SPEC_PATH}/ci/G2V32L1Fi
-        echo ${PYTHONPATH}
         export OMP_NUM_THREADS=1
         mpiexec -n 2 --allow-run-as-root ${SPEC_PATH}/xspec G2V32L1Fi.001.sp
         python3 -m py_spec.ci.test compare.h5 G2V32L1Fi.001.sp.h5 

.github/workflows/build_cmake.yml

@@ -8,7 +8,6 @@ jobs:
       fail-fast: false
     env:
       SPEC_PATH: ${{ github.workspace }}
-      PYTHONPATH: ${{ github.workspace }}/Utilities/pythontools
       OMPI_ALLOW_RUN_AS_ROOT: 1
       OMPI_ALLOW_RUN_AS_ROOT_CONFIRM: 1
     steps:
@@ -20,6 +19,8 @@ jobs:
         pip3 install --upgrade pip
         pip3 install --user ninja cmake scipy
         pip3 install --user numpy f90nml scikit-build scipy h5py matplotlib
+        cd ${{ github.workspace }}/Utilities/pythontools
+        pip3 install -v .
     - name: Build & Test
       uses: ashutoshvarma/action-cmake-build@master
       with:
@@ -39,14 +40,12 @@ jobs:
     - name: run_fast_cartesian
       run: |
         cd ${SPEC_PATH}/ci/G1V03L2Fi
-        echo ${PYTHONPATH}
         export OMP_NUM_THREADS=1
         mpiexec -n 2 --allow-run-as-root $SPEC_PATH/install/bin/xspec G1V03L2Fi.001.sp
         python3 -m py_spec.ci.test compare.h5 G1V03L2Fi.001.sp.h5 
     - name: run_fast_cylinder
       run: |
         cd ${SPEC_PATH}/ci/G2V32L1Fi
-        echo ${PYTHONPATH}
         export OMP_NUM_THREADS=1
         mpiexec -n 2 --allow-run-as-root $SPEC_PATH/install/bin/xspec G2V32L1Fi.001.sp
         python3 -m py_spec.ci.test compare.h5 G2V32L1Fi.001.sp.h5 

.github/workflows/py_spec.yml

@@ -19,7 +19,6 @@ jobs:
       working-directory: ${{ env.PY_SPEC_DIR }}
       run: |
         pip install --upgrade pip
-        pip3 install -r requirements.txt
         pip3 install setuptools wheel twine
 
     - name: Build py_spec

.github/workflows/python_wrapper.yml

@@ -6,26 +6,25 @@ jobs:
     name: python_wrapper build
 
     steps:
-    - uses: actions/checkout@v2
+    - uses: actions/checkout@v4
 
     # Python3 should be pre-installed on 'ubuntu-latest'
     - name: Python version info
       run: |
-        python3 --version
-        pip3 --version
+        python --version
+        pip --version
 
     - name: Install dependencies
       run: |
         sudo apt-get update
-        sudo apt-get install gfortran mpi-default-bin mpi-default-dev libhdf5-103 libhdf5-dev libfftw3-bin libfftw3-dev libopenblas0-openmp libopenblas-dev cmake ninja-build
-        pip3 install --user numpy f90nml scikit-build scipy
-        pip3 install --user git+https://github.com/zhucaoxiang/f90wrap
+        sudo apt-get install gfortran mpi-default-bin mpi-default-dev libhdf5-dev libfftw3-bin libfftw3-dev  libopenblas-dev cmake ninja-build
+        pip install  numpy f90nml scikit-build scipy meson meson-python
+        pip install f90wrap
 
     - name: Build python_wrapper
       run: |
-        python3 setup.py bdist_wheel
-        pip3 install --user dist/*.whl
+        pip install .
         
     - name: Test if installation is ok
       run: |
-        python3 -c "import spec; print('success')"
+        python -c "import spec; print('success')"

Utilities/pythontools/py_spec/__init__.py

@@ -1,5 +1,10 @@
-# import of all SPEC-related python scripts.
-__version__ = "3.3.3"
+try:
+    from importlib import metadata
+except ImportError:
+    # Running on pre-3.8 Python; use importlib-metadata package
+    import importlib_metadata as metadata
+
+__version__ = metadata.version(__package__ or __name__)
 
 from .ci import test
 from .input.spec_namelist import SPECNamelist

Utilities/pythontools/py_spec/output/_plot_pressure.py

@@ -9,8 +9,8 @@ def plot_pressure(self, normalize=True, ax=None, **kwargs):
     import numpy as np
     import matplotlib.pyplot as plt
 
-    pressure = self.input.physics.pressure * self.input.physics.pscale
-    tflux = self.output.tflux[: len(pressure)]
+    pressure = np.atleast_1d(self.input.physics.pressure) * self.input.physics.pscale
+    tflux = np.atleast_1d(self.output.tflux)[: len(pressure)]
     if not normalize:
         #  remove  mu_0
         pressure /= 4 * np.pi * 1.0e-7

Utilities/pythontools/py_spec/output/_processing.py

@@ -1,5 +1,6 @@
 import numpy as np
 from scipy import integrate
+import typing
 
 def get_RZ_derivatives(
     self,
@@ -504,7 +505,7 @@ def get_average_beta(self, ns=64, nt=64, nz=64):
     """Get beta averaged in plasma volume"""
 
     # Read pressure
-    press = self.input.physics.pressure * self.input.physics.pscale
+    press = np.atleast_1d(self.input.physics.pressure) * self.input.physics.pscale
 
     # Create coordinate grid
     nfp = self.input.physics.Nfp
@@ -606,11 +607,32 @@ def test_derivatives(self, lvol=0, s=0.3, t=0.4, z=0.5, delta=1e-6, tol=1e-6):
     print((g[0,1,0,:,:] - g[0,0,0,:,:])/ds/2-dg[0,0,0,1,:,:])
     print((g[0,0,1,:,:] - g[0,0,0,:,:])/ds/2-dg[0,0,0,2,:,:])
 
-def get_surface_current_density(self, lsurf:int=None, nt:int=64, nz:int=64):
+def _validate_lsurf(lsurf:np.ndarray, mvol:int)->np.ndarray:
+    """Check 
+    
+    Args:
+        - lsurf: Interface number(s), between 1 and Mvol-1. default is np.arange(1, mvol)
+        - mvol: Number of volumes
+    Returns:
+        - lsurf: 1d array of interface numbers
+        
+    Raises:
+        - ValueError: if input is outside of range or wrong type (lsurf)
+    """
+
+    if lsurf is None:
+        lsurf = np.arange(1,mvol)
+    else:
+        lsurf = np.atleast_1d(lsurf)
+    if (lsurf<1).any() or (lsurf>mvol-1).any(): raise ValueError('lsurf should be in [1,mvol-1]')
+
+    return lsurf
+
+def get_surface_current_density(self, lsurf:np.ndarray, nt:int=64, nz:int=64)->typing.Tuple[np.ndarray, np.ndarray, np.ndarray]:
     """Compute j_surf.B on each side of the provided interfaces
     
     Args:
-        - lsurf: Interface number, between 1 and Mvol-1.
+        - lsurf: Interface number(s), between 1 and Mvol-1. default is np.arange(1, mvol)
         - nt: Number of poloidal points
         - nz: Number of toroidal points
         
@@ -628,9 +650,7 @@ def get_surface_current_density(self, lsurf:int=None, nt:int=64, nz:int=64):
     nfp = self.input.physics.Nfp
 
     if mvol==1: raise ValueError('Mvol=1; no interface current!')
-    if not isinstance(lsurf, np.ndarray): raise ValueError('lsurf should be a np.ndarray')
-    if lsurf is None: raise ValueError('Need to provide lsurf')
-    if (lsurf<1).any() or (lsurf>mvol-1).any(): raise ValueError('lsurf should be in [1,mvol-1]')
+    lsurf = _validate_lsurf(lsurf, mvol)
     if nt<1: raise ValueError('nt should greater than zero')
     if nz<1: raise ValueError('nz should greater than zero')
 
@@ -683,11 +703,11 @@ def get_surface_current_density(self, lsurf:int=None, nt:int=64, nz:int=64):
 
     return j_dot_B, tarr, zarr
 
-def get_surface(self, lsurf:int=None, nt:int=64, nz:int=64):
+def get_surface(self, lsurf:np.ndarray=None, nt:int=64, nz:int=64):
     """Compute the surface area of a volume interface
     
     Args:
-        - lsurf: Interface number, between 1 and Mvol-1.
+        - lsurf: Interface number(s), between 1 and Mvol-1. default is np.arange(1, mvol)
         - nt: Number of poloidal points for integration, default is 64
         - nz: Number of toroidal points for integration, default is 64
         
@@ -700,15 +720,10 @@ def get_surface(self, lsurf:int=None, nt:int=64, nz:int=64):
 
     mvol = self.output.Mvol
     nfp = self.input.physics.Nfp
-    if lsurf is None:
-        lsurf = np.arange(1,mvol)
 
     if mvol==1: 
         raise ValueError('Mvol=1; no interface current!')
-    if not isinstance(lsurf, np.ndarray): 
-        raise ValueError('lsurf should be a np.ndarray')
-    if (lsurf<1).any() or (lsurf>mvol-1).any(): 
-        raise ValueError('lsurf should be in [1,mvol-1]')
+    lsurf = _validate_lsurf(lsurf, mvol)
     if nt<1: 
         raise ValueError('nt should greater than zero')
     if nz<1: 
@@ -744,15 +759,19 @@ def get_flux_surface_average( self, lsurf, f, tarr, zarr ):
     inner side of the interface (i.e. lvol=lsurf-1, sarr=1)
     
     Args:
-        - lsurf (1D numpy array): Interface number, between 1 and Mvol-1
+        - lsurf: Interface number(s), between 1 and Mvol-1. default is np.arange(1, mvol)
         - f (2D numpy array): function evaluated on a grid
         - tgrid (2D numpy array): theta grid
         - zgrid (2D numpy array): phi grid
     
     Returns>
         - fsavg (1D numpy array): The flux surface average of f on each surface
           given in lsurf
+    Raises:
+        - ValueError: if input is wrong (invalid lsurf)
     """
+
+    lsurf = _validate_lsurf(lsurf, self.output.Mvol)
 
     # Get jacobian
     output = np.zeros(lsurf.shape)

Utilities/pythontools/py_spec/output/spec.py

@@ -95,15 +95,9 @@ def __init__(self, *args, **kwargs):
         if isinstance(_content, h5py.File):
             _content.close()
 
-            # make sure that Lrad is always an array
-            if np.isscalar(self.input.physics.Lrad):
-                self.input.physics.Lrad = np.array([self.input.physics.Lrad])
-            # make sure that im always an array
-            if np.isscalar(self.output.im):
-                self.output.im = np.array([self.output.im])
-            # make sure that in_ is always an array
-            if np.isscalar(self.output.in_):
-                self.output.in_ = np.array([self.output.in_])
+            self.input.physics.Lrad = np.atleast_1d(self.input.physics.Lrad)
+            self.output.im = np.atleast_1d(self.output.im)
+            self.output.in_ = np.atleast_1d(self.output.in_)
 
             # these define the target dimensions in the radial direction
             Nvol = self.input.physics.Nvol

Utilities/pythontools/pyproject.toml

@@ -0,0 +1,31 @@
+[build-system]
+requires = ["setuptools"]
+build-backend = "setuptools.build_meta"
+
+
+[project]
+name="py_spec"
+version="3.3.5"
+dependencies = ["numpy>=1.21.1", 
+            "f90nml",
+            "h5py", 
+            "matplotlib",
+            "coilpy; python_version<'3.12'",
+            "scipy>=1.7.0"]
+description="SPEC(Stepped-Pressure Equilibrium Code) python utilities"
+readme="README.md"
+authors = [
+    { name = "Christopher Berg Smiet", email = "[email protected]" },
+    { name = "Caoxiang Zhu", email = "[email protected]" },
+    { name = "SPEC developers"}
+]
+maintainers = [
+    { name = "Christopher Berg Smiet", email = "[email protected]" },
+]
+classifiers=[
+    "Development Status :: 3 - Alpha",
+    "License :: OSI Approved :: GNU Lesser General Public License v3 or later (LGPLv3+)",
+    "Programming Language :: Python :: 3",
+    "Topic :: Scientific/Engineering",
+]
+license = {text = "GNU 3.0"}

Utilities/pythontools/requirements.txt

@@ -1,7 +1,6 @@
 h5py
 matplotlib
 f90nml
-numpy
 
 # Version 1.7.0 or higher is required - otherwise scipy.integrate.simpson does not exist
 scipy>=1.7.0

Utilities/pythontools/setup.py

@@ -1,23 +1,6 @@
 import setuptools
-from py_spec import __version__
-
-with open("README.md", "r") as fh:
-    long_description = fh.read()
 
 setuptools.setup(
     name="py_spec",
-    version=__version__,
-    description="SPEC(Stepped-Pressure Equilibrium Code) python utilities",
-    long_description=long_description,
-    long_description_content_type="text/markdown",
-    classifiers=[
-        "Development Status :: 3 - Alpha",
-        "License :: OSI Approved :: GNU Lesser General Public License v3 or later (LGPLv3+)",
-        "Programming Language :: Python :: 3",
-        "Topic :: Scientific/Engineering",
-    ],
-    url="https://princetonuniversity.github.io/SPEC/",
-    author="SPEC developers",
-    license="GNU 3.0",
-    packages=['py_spec', 'py_spec.input', 'py_spec.output', 'py_spec.ci']
+    packages=['py_spec', 'py_spec.input', 'py_spec.output', 'py_spec.ci', 'py_spec.math']
 )