Translate calls to low level functions (SeLaLib/LAPACK)

.github/workflows/testing.yml

@@ -27,51 +27,51 @@ on:
 
 jobs:
 
-  Coverage:
-
-    runs-on: ubuntu-latest
-
-    steps:
-      - uses: actions/checkout@v3
-      - name: Set up Python
-        uses: actions/setup-python@v4
-        with:
-            python-version: '3.10'
-      - name: Install dependencies
-        uses: ./.github/actions/linux_install
-      - name: Install SeLaLib
-        uses: ./.github/actions/selalib_installation
-      - name: Cache pip
-        uses: actions/cache@v3
-        with:
-          # This path is specific to Ubuntu
-          path: ~/.cache/pip
-          # Look to see if there is a cache hit for the corresponding requirements file
-          key: ${{ runner.os }}-pip-${{ hashFiles('requirements.txt') }}
-          restore-keys: |
-            ${{ runner.os }}-pip-
-            ${{ runner.os }}-
-      - name: Install python dependencies
-        uses: ./.github/actions/pip_installation
-      - name: Coverage install
-        uses: ./.github/actions/coverage_install
-      - name: Unit tests
-        run: |
-          coverage run -p --source=pygyro -m pytest pygyro -m "serial and not long" --short
-          make
-          coverage run -p --source=pygyro -m pytest pygyro -m serial
-          mpirun -n 1 coverage run -p mpi_tester.py pygyro -m parallel --mpisub
-          mpirun -n 4 --oversubscribe coverage run -p mpi_tester.py pygyro -m parallel --mpisub
-          mpirun -n 6 --oversubscribe coverage run -p mpi_tester.py pygyro -m parallel --mpisub
-      - name: Collect coverage information
-        continue-on-error: True
-        uses: ./.github/actions/coverage_collection
-      - name: Run codacy-coverage-reporter
-        uses: codacy/codacy-coverage-reporter-action@master
-        continue-on-error: True
-        with:
-          project-token: ${{ secrets.CODACY_PROJECT_TOKEN }}
-          coverage-reports: current_coverage.xml
+    #Coverage:
+
+    #  runs-on: ubuntu-latest
+
+    #  steps:
+    #    - uses: actions/checkout@v3
+    #    - name: Set up Python
+    #      uses: actions/setup-python@v4
+    #      with:
+    #          python-version: '3.10'
+    #    - name: Install dependencies
+    #      uses: ./.github/actions/linux_install
+    #    - name: Install SeLaLib
+    #      uses: ./.github/actions/selalib_installation
+    #    - name: Cache pip
+    #      uses: actions/cache@v3
+    #      with:
+    #        # This path is specific to Ubuntu
+    #        path: ~/.cache/pip
+    #        # Look to see if there is a cache hit for the corresponding requirements file
+    #        key: ${{ runner.os }}-pip-${{ hashFiles('requirements.txt') }}
+    #        restore-keys: |
+    #          ${{ runner.os }}-pip-
+    #          ${{ runner.os }}-
+    #    - name: Install python dependencies
+    #      uses: ./.github/actions/pip_installation
+    #    - name: Coverage install
+    #      uses: ./.github/actions/coverage_install
+    #    - name: Unit tests
+    #      run: |
+    #        coverage run -p --source=pygyro -m pytest pygyro -m "serial and not long" --short
+    #        make
+    #        coverage run -p --source=pygyro -m pytest pygyro -m serial
+    #        mpirun -n 1 coverage run -p mpi_tester.py pygyro -m parallel --mpisub
+    #        mpirun -n 4 --oversubscribe coverage run -p mpi_tester.py pygyro -m parallel --mpisub
+    #        mpirun -n 6 --oversubscribe coverage run -p mpi_tester.py pygyro -m parallel --mpisub
+    #    - name: Collect coverage information
+    #      continue-on-error: True
+    #      uses: ./.github/actions/coverage_collection
+    #    - name: Run codacy-coverage-reporter
+    #      uses: codacy/codacy-coverage-reporter-action@master
+    #      continue-on-error: True
+    #      with:
+    #        project-token: ${{ secrets.CODACY_PROJECT_TOKEN }}
+    #        coverage-reports: current_coverage.xml
 
   Pyccel_Fortran:
 

pygyro/advection/accelerated_advection_steps.py

@@ -1,6 +1,9 @@
+import numpy as np
 from typing import Final
 from pyccel.decorators import pure
 from ..splines.splines import Spline1D, Spline2D
+from ..splines.accelerated_spline_interpolators import solve_system_nonperiodic, solve_system_periodic, solve_2d_system
+from ..splines.sll_m_spline_matrix_periodic_banded import PeriodicBandedMatrix
 from ..initialisation.initialiser_funcs import f_eq
 
 
@@ -37,6 +40,7 @@ def poloidal_advection_step_expl(f: 'float[:,:]',
 
     multFactor = dt / B0
     multFactor_half = 0.5 * multFactor
+
     phi_spline.eval_vector(qPts, rPts, drPhi_0, 0, 1)
     phi_spline.eval_vector(qPts, rPts, dthetaPhi_0, 1, 0)
 
@@ -141,6 +145,24 @@ def v_parallel_advection_eval_step(f: 'float[:]', vPts: 'float[:]',
             f[i] = spl.eval(v)
 
 
+def v_parallel_advection_eval_step_loop(f: 'float[:,:,:]', vPts: 'float[:]',
+                                           rPos: 'float', vMin: 'float', vMax: 'float',
+                                        spl: Spline1D, bmat: 'float[:,:](order=F)', l: np.int32, u: np.int32, ipiv: 'int32[:]',
+                                        parGradVals: 'float[:,:,:]', dt: float, i: int,
+                                           CN0: 'float', kN0: 'float', deltaRN0: 'float', rp: 'float',
+                                           CTi: 'float', kTi: 'float', deltaRTi: 'float', bound: 'int'):
+    n1, n2, _ = f.shape
+    for j in range(n1):  # z
+        for k in range(n2):  # q
+            coeffs = spl.coeffs
+            solve_system_nonperiodic(f[j, k, :], coeffs, bmat, l, u, ipiv)
+
+            vPts -= parGradVals[i, j, k]*dt
+
+            v_parallel_advection_eval_step(f[j, k, :], vPts, rPos, vMin, vMax, spl,
+                                           CN0, kN0, deltaRN0, rp, CTi, kTi, deltaRTi, bound)
+
+
 def get_lagrange_vals(i: 'int', shifts: 'int[:]',
                       vals: 'float[:,:,:]', qVals: 'float[:]',
                       thetaShifts: 'float[:]', spl: Spline1D):
@@ -170,6 +192,28 @@ def flux_advection(nq: 'int', nr: 'int',
                 f[j, i] += coeffs[k]*vals[i, j, k]
 
 
+def flux_advection_loop(f: 'float[:,:,:,:]', thetaSpline: Spline1D, theta_offset: int, theta_splu: PeriodicBandedMatrix,
+                        shifts: 'int[:,:,:]', lagrange_vals: 'float[:,:,:]', qVals: 'float[:]',
+                        thetaShifts: 'float[:,:,:]', lagrange_coeffs: 'float[:,:,:]'):
+    nr, nv, nq, nz = f.shape
+
+    for rIdx in range(nr):  # r
+        for cIdx in range(nv):  # v
+            # find the values of the function at each required point
+            for k in range(nz):
+                solve_system_periodic(
+                    f[rIdx, cIdx, :, k], thetaSpline, theta_offset, theta_splu)
+
+                get_lagrange_vals(k, shifts[rIdx, cIdx],
+                                  lagrange_vals, qVals,
+                                  thetaShifts[rIdx, cIdx],
+                                  thetaSpline)
+
+            flux_advection(nq, nz, f[rIdx, cIdx],
+                           lagrange_coeffs[rIdx, cIdx],
+                           lagrange_vals)
+
+
 def poloidal_advection_step_impl(f: 'float[:,:]', dt: 'float', v: 'float', rPts: 'float[:]', qPts: 'float[:]',
                                  drPhi_0: 'float[:,:]', dthetaPhi_0: 'float[:,:]', drPhi_k: 'float[:,:]', dthetaPhi_k: 'float[:,:]',
                                  endPts_k1_q: 'float[:,:]', endPts_k1_r: 'float[:,:]', endPts_k2_q: 'float[:,:]', endPts_k2_r: 'float[:,:]',
@@ -294,3 +338,37 @@ def poloidal_advection_step_impl(f: 'float[:,:]', dt: 'float', v: 'float', rPts:
                     endPts_k2_q[i, j] = endPts_k2_q[i, j] % (2*pi)
                     f[i, j] = pol_spline.eval(
                         endPts_k2_q[i, j], endPts_k2_r[i, j])
+
+
+def poloidal_advection_loop(f: 'float[:,:,:,:]', phi: 'float[:,:,:]', dt: 'float',
+                            is_explicit: 'bool', rPts: 'float[:]', qPts: 'float[:]',
+                            vPts: 'float[:]', interp_wt: 'float[:,:]',
+                            r_bmat: 'float[:,:](order=F)', r_l: np.int32, r_u: np.int32,
+                            r_ipiv: 'int32[:]', theta_offset: int, theta_splu: PeriodicBandedMatrix,
+                            drPhi_0: 'float[:,:]', dthetaPhi_0: 'float[:,:]', drPhi_k: 'float[:,:]', dthetaPhi_k: 'float[:,:]',
+                            endPts_k1_q: 'float[:,:]', endPts_k1_r: 'float[:,:]', endPts_k2_q: 'float[:,:]', endPts_k2_r: 'float[:,:]',
+                            phi_spline: Spline2D, pol_spline: Spline2D,
+                            CN0: 'float', kN0: 'float', deltaRN0: 'float', rp: 'float', CTi: 'float', kTi: 'float', deltaRTi: 'float',
+                            B0: 'float', tol: 'float', nulBound: 'bool'):
+    _, nz, _, _ = f.shape
+    for j in range(nz):
+        solve_2d_system(phi[j], phi_spline, interp_wt, r_bmat, r_l, r_u,
+                        r_ipiv, theta_offset, theta_splu)
+
+        # Do step
+        for i, v in enumerate(vPts):
+            solve_2d_system(f[i, j], pol_spline, interp_wt, r_bmat, r_l, r_u,
+                            r_ipiv, theta_offset, theta_splu)
+
+            if is_explicit:
+                poloidal_advection_step_expl(f[i, j], float(dt), v, rPts, qPts,
+                                             drPhi_0, dthetaPhi_0, drPhi_k, dthetaPhi_k,
+                                             endPts_k1_q, endPts_k1_r, endPts_k2_q, endPts_k2_r,
+                                             phi_spline, pol_spline, CN0, kN0, deltaRN0, rp, CTi,
+                                             kTi, deltaRTi, B0, nulBound)
+            else:
+                poloidal_advection_step_impl(f[i, j], float(dt), v, rPts, qPts,
+                                             drPhi_0, dthetaPhi_0, drPhi_k, dthetaPhi_k,
+                                             endPts_k1_q, endPts_k1_r, endPts_k2_q, endPts_k2_r,
+                                             phi_spline, pol_spline, CN0, kN0, deltaRN0, rp, CTi,
+                                             kTi, deltaRTi, B0, tol, nulBound)

pygyro/advection/advection.py

@@ -4,12 +4,15 @@
 
 from ..splines.splines import BSplines, Spline1D, Spline2D
 from ..splines.spline_interpolators import SplineInterpolator1D, SplineInterpolator2D
+from ..splines.accelerated_spline_interpolators import solve_system_periodic, solve_system_nonperiodic
 from ..model.layout import Layout
 from ..model.grid import Grid
-from .accelerated_advection_steps import get_lagrange_vals, flux_advection, \
+from .accelerated_advection_steps import get_lagrange_vals, flux_advection, flux_advection_loop, \
+    v_parallel_advection_eval_step_loop, \
     v_parallel_advection_eval_step, \
     poloidal_advection_step_expl, \
-    poloidal_advection_step_impl
+    poloidal_advection_step_impl, \
+    poloidal_advection_loop
 
 
 def fieldline(theta, z_diff, iota, r, R0):
@@ -278,7 +281,9 @@ def step(self, f: np.ndarray, cIdx: int, rIdx: int = 0):
 
         # find the values of the function at each required point
         for i in range(self._nPoints[1]):
-            self._interpolator.compute_interpolant(f[:, i], self._thetaSpline)
+            solve_system_periodic(
+                f[:, i], self._thetaSpline, self._interpolator._offset, self._interpolator._splu)
+            # self._interpolator.compute_interpolant(f[:, i], self._thetaSpline)
 
             get_lagrange_vals(i, self._shifts[rIdx, cIdx],
                               self._LagrangeVals, self._points[0],
@@ -291,9 +296,13 @@ def step(self, f: np.ndarray, cIdx: int, rIdx: int = 0):
 
     def gridStep(self, grid: Grid):
         assert grid.getLayout(grid.currentLayout).dims_order == (0, 3, 1, 2)
-        for i, _ in grid.getCoords(0):  # r
-            for j, _ in grid.getCoords(1):  # v
-                self.step(grid.get2DSlice(i, j), j)
+        flux_advection_loop(grid.getAllData(), self._thetaSpline,
+                            self._interpolator._offset, self._interpolator._splu,
+                            self._shifts, self._LagrangeVals, self._points[0],
+                            self._thetaShifts, self._lagrangeCoeffs)
+        # for i, _ in grid.getCoords(0):  # r
+        #    for j, _ in grid.getCoords(1):  # v
+        #        self.step(grid.get2DSlice(i, j), j)
 
 
 class VParallelAdvection:
@@ -360,7 +369,9 @@ def step(self, f: np.ndarray, dt: float, c: float, r: float):
 
         """
         assert f.shape == self._nPoints
-        self._interpolator.compute_interpolant(f, self._spline)
+        # self._interpolator.compute_interpolant(f, self._spline)
+        solve_system_nonperiodic(f, self._spline.coeffs, self._interpolator._bmat,
+                                 self._interpolator._l, self._interpolator._u, self._interpolator._ipiv)
 
         v_parallel_advection_eval_step(f, self._points-c*dt, r, self._points[0],
                                        self._points[-1], self._spline,
@@ -373,17 +384,33 @@ def gridStep(self, grid: Grid, phi: Grid, parGrad: ParallelGradient, parGradVals
         for i, r in grid.getCoords(0):
             parGrad.parallel_gradient(
                 np.real(phi.get2DSlice(i)), i, parGradVals[i])
-            for j, _ in grid.getCoords(1):  # z
-                for k, _ in grid.getCoords(2):  # q
-                    self.step(grid.get1DSlice(
-                        i, j, k), dt, parGradVals[i, j, k], r)
+            v_parallel_advection_eval_step_loop(grid.get3DSlice(i), self._points, r, self._points[0],
+                                                self._points[-1], self._spline, self._interpolator._bmat,
+                                                self._interpolator._l, self._interpolator._u,
+                                                self._interpolator._ipiv, parGradVals, dt, i,
+                                                self._constants.CN0, self._constants.kN0,
+                                                self._constants.deltaRN0, self._constants.rp,
+                                                self._constants.CTi, self._constants.kTi,
+                                                self._constants.deltaRTi, self._edgeType)
+            # for j, _ in grid.getCoords(1):  # z
+            #    for k, _ in grid.getCoords(2):  # q
+            #        self.step(grid.get1DSlice(
+            #            i, j, k), dt, parGradVals[i, j, k], r)
 
     def gridStepKeepGradient(self, grid: Grid, parGradVals, dt: float):
         for i, r in grid.getCoords(0):
-            for j, _ in grid.getCoords(1):  # z
-                for k, _ in grid.getCoords(2):  # q
-                    self.step(grid.get1DSlice(
-                        i, j, k), dt, parGradVals[i, j, k], r)
+            v_parallel_advection_eval_step_loop(grid.get3DSlice(i), self._points, r, self._points[0],
+                                                self._points[-1], self._spline, self._interpolator._bmat,
+                                                self._interpolator._l, self._interpolator._u,
+                                                self._interpolator._ipiv, parGradVals, dt, i,
+                                                self._constants.CN0, self._constants.kN0,
+                                                self._constants.deltaRN0, self._constants.rp,
+                                                self._constants.CTi, self._constants.kTi,
+                                                self._constants.deltaRTi, self._edgeType)
+            # for j, _ in grid.getCoords(1):  # z
+            #    for k, _ in grid.getCoords(2):  # q
+            #        self.step(grid.get1DSlice(
+            #            i, j, k), dt, parGradVals[i, j, k], r)
 
 
 class PoloidalAdvection:
@@ -508,14 +535,30 @@ def gridStep(self, grid: Grid, phi: Grid, dt: float):
         phiLayout = phi.getLayout(grid.currentLayout)
         assert gridLayout.dims_order[1:] == phiLayout.dims_order
         assert gridLayout.dims_order == (3, 2, 1, 0)
+        vPts = grid.getCoordVals(0)
+        poloidal_advection_loop(grid.getAllData(), np.real(phi.getAllData()), float(dt), self._explicit,
+                                self._points[1], self._points[0], vPts,
+                                self._interpolator._bwork, self._interpolator._interp2._bmat,
+                                self._interpolator._interp2._l, self._interpolator._interp2._u,
+                                self._interpolator._interp2._ipiv, self._interpolator._interp1._offset,
+                                self._interpolator._interp1._splu,
+                                self._drPhi_0, self._dqPhi_0, self._drPhi_k,
+                                self._dqPhi_k, self._endPts_k1_q,
+                                self._endPts_k1_r, self._endPts_k2_q,
+                                self._endPts_k2_r, self._phiSplines[0],
+                                self._spline, self._constants.CN0,
+                                self._constants.kN0, self._constants.deltaRN0,
+                                self._constants.rp, self._constants.CTi,
+                                self._constants.kTi, self._constants.deltaRTi,
+                                self._constants.B0, self._TOL, self._nulEdge)
         # Evaluate splines
-        for j, _ in grid.getCoords(1):  # z
-            self._interpolator.compute_interpolant(
-                np.real(phi.get2DSlice(j)), self._phiSplines[j])
+        # for j, _ in grid.getCoords(1):  # z
+        #    self._interpolator.compute_interpolant(
+        #        np.real(phi.get2DSlice(j)), self._phiSplines[j])
         # Do step
-        for i, v in grid.getCoords(0):
-            for j, _ in grid.getCoords(1):  # z
-                self.step(grid.get2DSlice(i, j), dt, self._phiSplines[j], v)
+        # for i, v in grid.getCoords(0):
+        #    for j, _ in grid.getCoords(1):  # z
+        #        self.step(grid.get2DSlice(i, j), dt, self._phiSplines[j], v)
 
     def gridStep_SplinesUnchanged(self, grid: Grid, dt: float):
         gridLayout = grid.getLayout(grid.currentLayout)