Merge pull request #20 from LLNL/task/bgarcia/freeqdsk_support

.github/workflows/ingrid_ci.yml

@@ -16,7 +16,7 @@ jobs:
     strategy:
       fail-fast: false
       matrix:
-        python-version: ["3.7", "3.8", "3.9", "3.10", "3.11"]
+        python-version: ["3.8", "3.9", "3.10", "3.11"]
 
     steps:
     - uses: actions/checkout@v3

INGRID/gui/ingrid_gui.py

@@ -1,12 +1,7 @@
 import matplotlib
-import os
 
 try:
-    if os.environ.get('DISPLAY','') == '':
-        print('no display found. Using non-interactive Agg backend')
-        matplotlib.use('Agg')
-    else:
-        matplotlib.use("TkAgg")
+    matplotlib.use("TkAgg")
 except:
    active_backend = matplotlib.get_backend()
    msg = 'Warning: Could not set matplotlib backend to "TkAgg". '

INGRID/ingrid.py

@@ -17,17 +17,14 @@
     pass
 import matplotlib.pyplot as plt
 import pathlib
-import inspect
 import functools
 from scipy.optimize import root, minimize
 
 import yaml as yml
 import os
 from pathlib import Path
 from time import time, perf_counter
-from collections import OrderedDict
 
-from INGRID.OMFITgeqdsk import OMFITgeqdsk
 from INGRID.interpol import EfitData
 from INGRID.utils import IngridUtils
 from INGRID.topologies.snl import SNL
@@ -38,7 +35,6 @@
 from INGRID.topologies.sf135 import SF135
 from INGRID.topologies.sf165 import SF165
 from INGRID.topologies.udn import UDN
-from INGRID.line_tracing import LineTracing
 from INGRID.geometry import Point, Line
 
 
@@ -128,10 +124,6 @@ def wrapper_timer(*args, **kwargs):
             return value
         return wrapper_timer
 
-    def LoadEFIT(self, fpath: str) -> None:
-        self.settings['eqdsk'] = fpath
-        self.OMFIT_read_psi()
-
     def StartGUI(self, test_initialization: bool = False) -> None:
         """
         Start GUI for Ingrid.
@@ -1291,7 +1283,7 @@ def StartSetup(self, **kwargs) -> None:
                 + 'Must be <= 2).'
             raise ValueError(v_error_str)
 
-        self.LoadEFIT(self.settings['eqdsk'])
+        self.LoadGEQDSK(self.settings['eqdsk'])
         self.AutoRefineMagAxis()
         self.AutoRefineXPoint()
         if topology == 'DNL':

INGRID/interpol.py

@@ -1,10 +1,7 @@
-#!/usr/bin/env pythonu
-# -*- coding: utf-8 -*-
 """
 Module containing EfitData class for handling all interpolation
 related computations.
 """
-from __future__ import division, print_function, absolute_import
 import numpy as np
 import matplotlib.pyplot as plt
 from scipy.ndimage import zoom
@@ -42,9 +39,22 @@ class EfitData:
         Specify the title of the figure the data will be plotted on.
     """
 
-    def __init__(self, rmin=0.0, rmax=1.0, nr=10, zmin=0.0, zmax=2.0, nz=20,
-                 rcenter=1.6955000, bcenter=-2.1094041, rlimiter=None, zlimiter=None,
-                 rmagx=0.0, zmagx=0.0, name='unnamed', parent=None):
+    def __init__(self, 
+                 rmin=0.0, 
+                 rmax=1.0, 
+                 nr=10, 
+                 zmin=0.0, 
+                 zmax=2.0, 
+                 nz=20,
+                 rcenter=1.6955000, 
+                 bcenter=-2.1094041, 
+                 rlimiter=None, 
+                 zlimiter=None,
+                 rmagx=0.0, 
+                 zmagx=0.0, 
+                 name='unnamed', 
+                 parent=None):
+
         r, dr = np.linspace(rmin, rmax, nr, retstep=True)
         z, dz = np.linspace(zmin, zmax, nz, retstep=True)
         rgrid, zgrid = np.meshgrid(r, z, indexing='ij')

INGRID/utils.py

@@ -15,20 +15,15 @@
 except:
     pass
 import matplotlib.pyplot as plt
-import pathlib
 import inspect
-from scipy.optimize import root, minimize
+from scipy.optimize import root
 
-import yaml as yml
-import os
 from pathlib import Path
-from time import time
-from collections import OrderedDict
+from freeqdsk import geqdsk
 
-from INGRID.OMFITgeqdsk import OMFITgeqdsk
 from INGRID.interpol import EfitData
 from INGRID.line_tracing import LineTracing
-from INGRID.geometry import Point, Line, Patch, segment_intersect, orientation_between
+from INGRID.geometry import Point, Line, Patch, orientation_between
 
 
 class IngridUtils():
@@ -99,9 +94,6 @@ class IngridUtils():
         Dictionary containing Line objects that correspond to target plates
         loaded by the user.
 
-    OMFIT_psi : OMFITgeqdsk
-        Instance of class used to interface G files generated by EFIT.
-
     PsiUNorm : EfitData
         `EfitData` class object containing unormalized
         psi data from provided neqdsk file in settings.
@@ -463,46 +455,59 @@ def ProcessPaths(self) -> None:
                 self.settings['patch_data']['preferences']['new_fname'] = str((path_obj / self.settings['patch_data']['preferences']['new_fname']).absolute())
                 continue
 
-    def OMFIT_read_psi(self) -> None:
+    def LoadGEQDSK(self, geqdsk_path: str) -> None:
         """
         Python class to read the psi data in from an ascii file.
 
         Saves the boundary information and generates an EfitData instance.
         """
 
-        g = OMFITgeqdsk(self.settings['eqdsk'])
-
-        nxefit = g['NW']
-        nyefit = g['NH']
-        rdim = g['RDIM']
-        zdim = g['ZDIM']
-        zmid = g['ZMID']
-        rgrid1 = g['RLEFT']
-
-        rcenter = g['RCENTR']
-        bcenter = g['BCENTR']
-
-        rmagx = g['RMAXIS']
-        zmagx = g['ZMAXIS']
-
-        rlimiter = g['RLIM']
-        zlimiter = g['ZLIM']
-
-        psi = g['PSIRZ'].T
-
-        # calc array for r and z
-        rmin = rgrid1
+        with open(geqdsk_path, 'r') as f:
+            geqdsk_data = geqdsk.read(f)
+
+        #
+        # Extract quantities needed to initialize EfitData class
+        #
+        nx       = geqdsk_data['nx']
+        ny       = geqdsk_data['ny']
+        rdim     = geqdsk_data['rdim']
+        zdim     = geqdsk_data['zdim']
+        zmid     = geqdsk_data['zmid']
+        rleft    = geqdsk_data['rleft']
+        rcenter  = geqdsk_data['rcentr']
+        bcenter  = geqdsk_data['bcentr']
+        rmagx    = geqdsk_data['rmagx']
+        zmagx    = geqdsk_data['zmagx']
+        rlimiter = geqdsk_data['rlim']
+        zlimiter = geqdsk_data['zlim']
+        psi      = geqdsk_data['psi']
+
+        #
+        # Derived values
+        #
+        rmin = rleft
         rmax = rmin + rdim
         zmin = (zmid - 0.5 * zdim)
         zmax = zmin + zdim
 
-        # reproduce efit grid
-        self.PsiUNorm = EfitData(rmin, rmax, nxefit,
-                                 zmin, zmax, nyefit,
-                                 rcenter, bcenter,
-                                 rlimiter, zlimiter,
-                                 rmagx, zmagx,
-                                 name='Efit Data', parent=self)
+        #
+        # Reproduce efit grid
+        #
+        self.PsiUNorm = EfitData(
+                            rmin=rmin, 
+                            rmax=rmax, 
+                            nr=nx,
+                            zmin=zmin, 
+                            zmax=zmax, 
+                            nz=ny,
+                            rcenter=rcenter, 
+                            bcenter=bcenter,
+                            rlimiter=rlimiter, 
+                            zlimiter=zlimiter,
+                            rmagx=rmagx, 
+                            zmagx=zmagx,
+                            name='Efit Data', 
+                            parent=self)
 
         self.PsiUNorm.init_bivariate_spline(self.PsiUNorm.r[:, 0], 
                                             self.PsiUNorm.z[0, :], 
@@ -511,7 +516,15 @@ def OMFIT_read_psi(self) -> None:
         if self.settings['grid_settings']['rmagx'] is None or self.settings['grid_settings']['zmagx'] is None:
             self.settings['grid_settings']['rmagx'], self.settings['grid_settings']['zmagx'] = (self.PsiUNorm.rmagx, self.PsiUNorm.zmagx)
 
-        self.OMFIT_psi = g
+        #
+        # Save the stored data in the object
+        #
+        self.geqdsk_data = geqdsk_data
+        #
+        # Update the eqdsk file referenced in settings to that of the loaded data
+        #
+        self.settings['eqdsk'] = geqdsk_path
+
 
     def ParseTxtCoordinates(self, fpath: str, rshift: float = 0.0, zshift: float = 0.0) -> tuple:
         """
@@ -628,43 +641,70 @@ def SetLimiter(self, fpath: str = '', coordinates: list = [], rshift: float = 0.
 
         use_efit_bounds = self.settings['limiter']['use_efit_bounds']
 
+        #
+        # TODO: Refactor the codebase to leverage numpy arrays....
+        #
+        coordinates = np.array(coordinates).T
+
+        #
+        # If we have no coordinates or fpath provided, we need to
+        # leverage efit boundaries. Overrides any user settings.
+        #
+        if coordinates.shape[0] == 0 and fpath == '':
+            use_efit_bounds = True
+
         if fpath not in ['', '.']:
             try:
                 print('# Processing file for limiter data : {}'.format(fpath))
-                self.OMFIT_psi['RLIM'], self.OMFIT_psi['ZLIM'] = self.ParseTxtCoordinates(fpath)
+                self.geqdsk_data['rlim'], self.geqdsk_data['rlim'] = self.ParseTxtCoordinates(fpath)
             except:
                 raise ValueError(f"# Error in method 'SetLimiter' with fpath={fpath}")
 
-        elif coordinates == []:
-            g = OMFITgeqdsk(self.settings['eqdsk'])
-            RLIM, ZLIM = g['RLIM'], g['ZLIM']
-            self.OMFIT_psi['RLIM'], self.OMFIT_psi['ZLIM'] = RLIM + rshift, ZLIM + zshift
+        #
+        # List of coordinates provided for initialization
+        #
+        elif len(coordinates) > 0:
+            self.geqdsk_data['rlim'] = coordinates[:, 0] + rshift
+            self.geqdsk_data['zlim'] = coordinates[:, 1] + zshift
 
-        elif coordinates != []:
-            RLIM, ZLIM = coordinates
-            self.OMFIT_psi['RLIM'], self.OMFIT_psi['ZLIM'] = [r + rshift for r in RLIM], [z + zshift for z in ZLIM]
-
-        if len(RLIM) == 0 or len(ZLIM) == 0:
-            use_efit_bounds = True
+        #
+        # Empty list of coordinates falls back on using eqdsk limiter settings
+        #
+        else:
+            self.LoadGEQDSK(geqdsk_path=self.settings['eqdsk'])
+            self.geqdsk_data['rlim'] += rshift
+            self.geqdsk_data['zlim'] += zshift
 
         if use_efit_bounds:
+
+            #
+            # The bounding box starts with using the computational domain bounds.
+            # Buffer values shrink the bounding box to allow us to avoid stepping
+            # outside the domain during integration/line-tracing
+            #
             rbuff = self.settings['limiter']['efit_buffer_r']
             zbuff = self.settings['limiter']['efit_buffer_z']
-            coordinates = [(self.PsiUNorm.rmin + rbuff, self.PsiUNorm.zmin + zbuff),
-                           (self.PsiUNorm.rmax - rbuff, self.PsiUNorm.zmin + zbuff),
-                           (self.PsiUNorm.rmax - rbuff, self.PsiUNorm.zmax - zbuff),
-                           (self.PsiUNorm.rmin + rbuff, self.PsiUNorm.zmax - zbuff),
-                           (self.PsiUNorm.rmin + rbuff, self.PsiUNorm.zmin + zbuff)]
-
-            RLIM, ZLIM = [], []
-            for c in coordinates:
-                r, z = c
-                RLIM.append(r)
-                ZLIM.append(z)
-
-            self.OMFIT_psi['RLIM'], self.OMFIT_psi['ZLIM'] = np.array(RLIM) + rshift, np.array(ZLIM) + zshift
 
-        self.LimiterData = Line([Point(p) for p in zip(self.OMFIT_psi['RLIM'], self.OMFIT_psi['ZLIM'])])
+            #
+            # Define the bounding box vertices as:
+            # - LL : Lower left
+            # - LR : Lower right
+            # - UL : Upper left
+            # - UR : Upper right
+            #
+            LL = np.array([self.PsiUNorm.rmin + rbuff, self.PsiUNorm.zmin + zbuff])
+            LR = np.array([self.PsiUNorm.rmax - rbuff, self.PsiUNorm.zmin + zbuff])
+            UL = np.array([self.PsiUNorm.rmin + rbuff, self.PsiUNorm.zmax - zbuff])
+            UR = np.array([self.PsiUNorm.rmax - rbuff, self.PsiUNorm.zmax - zbuff])
+
+            #
+            # Define the simple limiter path and translate by shift values
+            #
+            limiter_path = np.vstack([LL, LR, UR, UL, LL]) 
+            self.geqdsk_data['rlim'] = limiter_path[:, 0] + rshift
+            self.geqdsk_data['zlim'] = limiter_path[:, 1] + zshift
+
+        self.LimiterData = Line([Point(p) for p in zip(self.geqdsk_data['rlim'], self.geqdsk_data['zlim'])])
 
     def SetTargetPlate(self, settings: dict, rshift: float = 0.0, zshift: float = 0.0) -> None:
         """
@@ -1278,7 +1318,7 @@ def ReconstructPatches(self, raw_patch_list: list) -> dict:
 
             patches[patch.patch_name] = patch
 
-        patches = OrderedDict([(k, v) for k, v in patches.items()])
+        patches = {k: v for k, v in patches.items()}
 
         return patches
 

docs/source/requirements.txt

@@ -12,7 +12,6 @@ mkl-random==1.1.0
 mkl-service==2.3.0
 mpmath==1.3.0
 numpy==1.22.0
-pycrypto==2.6.1
 PyInstaller==3.6
 pyparsing==2.4.2
 python-dateutil==2.8.0

setup.py

@@ -28,30 +28,27 @@
 
     # This should be a valid email address corresponding to the author listed
     # above.
-    author_email='[email protected], [email protected]',
+    author_email='[email protected], [email protected]',
     classifiers=[
 	'Development Status :: 5 - Production/Stable',
         #'Intended Audience :: Any UEDGE users',
         'Topic :: Software Development',
         'License :: OSI Approved :: MIT License',
-        'Programming Language :: Python :: 3',
-        'Programming Language :: Python :: 3.5',
-        'Programming Language :: Python :: 3.6',
-        'Programming Language :: Python :: 3.7',
         'Programming Language :: Python :: 3.8',
-	'Programming Language :: Python :: 3.9',
+	'Programming Language :: Python :: 3.9'
     ],
     keywords='grid, mesh, generator, tokamak, edge, plasma, efit, uedge',
     packages=find_packages(),
-    python_requires='>=3.5, <4',
+    python_requires='>=3.8, <4',
     install_requires=[
         'm2r2',
         'scipy >= 1.3.1',
         'numpy',
         'pyyaml',
         'matplotlib',
         'tk',
-        'sympy'
+        'sympy',
+        'freeqdsk'
     ]
 )