Major changes for optimization

README.md

@@ -27,6 +27,13 @@ tested in both repositories before release.
 
 All of the examples for getting started can be found within the `freegsnke/examples` directory.
 
+## Parallel execution
+
+FreeGS4E contains multiple parallelized components for faster execution, through the use of multithreading. 
+
+The number of threads can be controlled through the environment variable `OMP_NUM_THREADS` (recommended). Note that by default the number of threads is limited to 32: for use in HPC systems, it is recommended to set the environment variable `NUMEXPR_MAX_THREADS` to the number of available cores.
+
+More fine grained parallel control can be achieved by using the environment variable `NUMEXPR_NUM_THREADS` or programatically through the use of the function `set_num_threads()` in `freegs4e.parallel_funcs`. These will change the number of threads used by parallel functions in FreeGS4E without affecting the number of threads used by OMP-based libraries (e.g. numpy). 
 
 ## Contributing
 

freegs4e/__init__.py

@@ -36,7 +36,7 @@
 
 __version__ = metadata.version("freegs4e")
 
-from . import control, jtor, machine, plotting
+from . import control, jtor, machine, multigrid, plotting
 from .dump import OutputFile
 from .equilibrium import Equilibrium
 from .picard import solve

freegs4e/equilibrium.py

@@ -31,13 +31,9 @@
 from scipy.optimize import least_squares
 from scipy.spatial.distance import pdist, squareform
 
-from . import critical, machine, multigrid, polygons  # multigrid solver
+from . import critical, machine, polygons
 from .boundary import fixedBoundary, freeBoundary  # finds free-boundary
-from .gradshafranov import (  # operators which define the G-S equation
-    GSsparse,
-    GSsparse4thOrder,
-    mu0,
-)
+from .gradshafranov import mu0
 from .plotting import plotEquilibrium
 
 
@@ -49,7 +45,7 @@ class Equilibrium:
 
     def __init__(
         self,
-        tokamak=machine.EmptyTokamak(),
+        tokamak=None,
         Rmin=0.1,
         Rmax=2.0,
         Zmin=-1.0,
@@ -59,7 +55,7 @@ def __init__(
         boundary=freeBoundary,
         psi=None,
         current=0.0,
-        order=4,
+        order=None,
     ):
         """
         Initializes a plasma equilibrium.
@@ -90,12 +86,15 @@ def __init__(
         current : float
             Plasma current [A].
         order : int
-            Order of differential operators used in calculations.
-            Must be either 2 or 4.
+            Deprecated. Kept for backwards compatibility.
+
         """
 
         # assign tokamak object
-        self.tokamak = tokamak
+        if tokamak:
+            self.tokamak = tokamak
+        else:
+            self.tokamak = machine.EmptyTokamak()
 
         # assign bounds of computational domain
         if Rmin > Rmax:
@@ -127,29 +126,24 @@ def __init__(
 
         # generate Greens function mappings (used
         # in self.psi() to speed up calculations)
-        self._pgreen = tokamak.createPsiGreens(self.R, self.Z)
-        self._vgreen = tokamak.createPsiGreensVec(self.R, self.Z)
+        self._pgreen = self.tokamak.createPsiGreens(self.R, self.Z)
+        self._vgreen = self.tokamak.createPsiGreensVec(self.R, self.Z)
         # self._updatePlasmaPsi(psi)  # Needs to be after _pgreen
 
         # assign plasma current
         self._current = current
 
-        # deinfe the GS solver
-        if order == 2:
-            generator = GSsparse(Rmin, Rmax, Zmin, Zmax)
-        elif order == 4:
-            generator = GSsparse4thOrder(Rmin, Rmax, Zmin, Zmax)
-        else:
-            raise ValueError(
-                "Invalid choice of order ({}). Valid values are 2 or 4.".format(
-                    order
-                )
+        # order attribute is kept for backwards compatibility only
+        if order is not None:
+            warnings.warn(
+                "Order attribute of Equilibrium objects is deprecated. No solver is created inside Equilibrium.",
+                DeprecationWarning,
             )
-        self.order = order
+        else:
+            order = 4
 
-        self._solver = multigrid.createVcycle(
-            nx, ny, generator, nlevels=1, ncycle=1, niter=2, direct=True
-        )
+        self._order = order
+        self.__solver = None
 
         # assign initial guess for plasma flux (if None)
         if psi is None:
@@ -166,6 +160,38 @@ def __init__(
             )
         self._updatePlasmaPsi(psi)
 
+    @property
+    def order(self):
+        warnings.warn(
+            "Order attribute of Equilibrium objects is deprecated, as no solver is created",
+            DeprecationWarning,
+        )
+        return self._order
+
+    @order.setter
+    def order(self, value):
+        warnings.warn(
+            "Order attribute of Equilibrium objects is deprecated, as no solver is created",
+            DeprecationWarning,
+        )
+        self._order = value
+
+    @property
+    def _solver(self):
+        warnings.warn(
+            "Solver attribute of Equilibrium objects is deprecated. Equilibrium contains no linear solver.",
+            DeprecationWarning,
+        )
+        return self.__solver
+
+    @order.setter
+    def _solver(self, value):
+        warnings.warn(
+            "Solver attribute of Equilibrium objects is deprecated. Equilibrium contains no linear solver.",
+            DeprecationWarning,
+        )
+        self.__solver = value
+
     def create_psi_plasma_default(
         self, adaptive_centre=False, gpars=(0.5, 0.5, 0, 2)
     ):
@@ -174,63 +200,14 @@ def create_psi_plasma_default(
         """
         return PsiGuessGaussian(adaptive_centre, gpars)(self)
 
-    def setSolverVcycle(self, nlevels=1, ncycle=1, niter=1, direct=True):
-        """
-        Sets a new linear solver based on the multigrid scheme.
-
-        This method configures a multigrid V-cycle solver and assigns it to
-        `self._solver`.
-
-        Parameters
-        ----------
-        nlevels : int
-            Number of resolution levels, including the original.
-        ncycle : int
-            Number of V-cycles to use.
-        niter : int
-            Number of linear solver (Jacobi) iterations per level.
-        direct : bool
-            If True, uses a direct solver at the coarsest level.
-
-        Returns
-        -------
-        None
-            This function modifies `self._solver` but does not return a value.
-        """
-
-        # set the solver
-        self._solver = multigrid.createVcycle(
-            nx=self.nx,
-            ny=self.ny,
-            generator=GSsparse(self.Rmin, self.Rmax, self.Zmin, self.Zmax),
-            nlevels=nlevels,
-            ncycle=ncycle,
-            niter=niter,
-            direct=direct,
-        )
+    def setSolverVcycle(
+        self, nlevels=1, ncycle=1, niter=1, direct=True, order=4
+    ):
+        """Deprecated. Equilibrium objects do not contain solvers"""
+        self.setSolver(None)
 
     def setSolver(self, solver):
-        """
-        Sets the linear solver to use. The given object/function must have a __call__ method
-        which takes two inputs:
-
-        solver(x, b)
-
-        where x is the initial guess and b is the right hand side (this should solve Ax = b,
-        returning the result).
-
-
-        Parameters
-        ----------
-        solver : object
-            The solver object.
-
-        Returns
-        -------
-        None
-            This function modifies `self._solver` but does not return a value.
-        """
-
+        """Deprecated. Equilibrium objects do not contain solvers"""
         self._solver = solver
 
     def callSolver(self, psi, rhs):