diff --git a/.github/workflows/testing.yml b/.github/workflows/testing.yml index 26fa0463..78c7631d 100644 --- a/.github/workflows/testing.yml +++ b/.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: diff --git a/pygyro/advection/accelerated_advection_steps.py b/pygyro/advection/accelerated_advection_steps.py index 4800543a..10e5b503 100644 --- a/pygyro/advection/accelerated_advection_steps.py +++ b/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) diff --git a/pygyro/advection/advection.py b/pygyro/advection/advection.py index 5a70f893..7a751348 100644 --- a/pygyro/advection/advection.py +++ b/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) diff --git a/pygyro/splines/accelerated_spline_interpolators.py b/pygyro/splines/accelerated_spline_interpolators.py index 42dc80af..b60abf03 100644 --- a/pygyro/splines/accelerated_spline_interpolators.py +++ b/pygyro/splines/accelerated_spline_interpolators.py @@ -1,7 +1,7 @@ from typing import TypeVar, Final import numpy as np from pyccel.stdlib.internal.lapack import dgbtrs, zgbtrs -from .splines import Spline1D, BSplines +from .splines import Spline1D, Spline2D, BSplines from .sll_m_spline_matrix_periodic_banded import PeriodicBandedMatrix T = TypeVar('T', float, complex) @@ -46,5 +46,53 @@ def solve_system_nonperiodic(ug: 'Final[T[:]]', c: 'T[:]', bmat: 'Final[T[:,:](o else: zgbtrs('N', np.int32(bmat.shape[1]), l, u, np.int32(1), bmat, np.int32( bmat.shape[0]), ipiv, c, np.int32(c.shape[0]), sinfo) + assert sinfo == 0 return sinfo + + +def solve_2d_system(ug: 'float[:,:]', spl: Spline2D, wt: 'float[:,:]', + r_bmat: 'float[:,:](order=F)', r_l: np.int32, r_u: np.int32, r_ipiv: 'int32[:]', + theta_offset: int, theta_splu: PeriodicBandedMatrix): + basis1 = spl.basis1 + basis2 = spl.basis2 + n1, n2 = basis1.nbasis, basis2.nbasis + p1, _ = basis1.degree, basis2.degree + assert ug.shape[0] == n1 + assert ug.shape[1] == n2 + + spline1 = Spline1D(basis1) + + w = spl.coeffs + + s1, s2 = w.shape + + # Cycle over x1 position and interpolate f along x2 direction. + # Work on spl.coeffs + sinfo: np.int32 + dgbtrs('N', np.int32(r_bmat.shape[1]), r_l, r_u, np.int32( + n1), r_bmat, np.int32(r_bmat.shape[0]), r_ipiv, ug.T, np.int32(s2), sinfo) + assert sinfo == 0 + + # Transpose coefficients to self._bwork + for i1 in range(n1): + for i2 in range(n2): + wt[i2, i1] = ug[i1, i2] + + # Cycle over x2 position and interpolate w along x1 direction. + # Work on self._bwork + for i2 in range(n2): + solve_system_periodic(wt[i2, :n1], spline1, theta_offset, theta_splu) + # self._interp1.compute_interpolant(wt[i2, :n1], self._spline1) + c = spline1.coeffs + wt[i2, :] = c + + # Transpose coefficients to spl.coeffs + for i1 in range(s1): + for i2 in range(s2): + w[i1, i2] = wt[i2, i1] + + # x1-periodic only: "wrap around" coefficients onto extended array + for i1 in range(p1): + for i2 in range(s2): + w[n1 + i1, i2] = w[i1, i2] diff --git a/pygyro/splines/spline_interpolators.py b/pygyro/splines/spline_interpolators.py index 92364b1a..dc039825 100644 --- a/pygyro/splines/spline_interpolators.py +++ b/pygyro/splines/spline_interpolators.py @@ -10,7 +10,7 @@ from .splines import BSplines, Spline2D, Spline1D, Spline1DComplex from .spline_eval_funcs import nu_find_span, nu_basis_funs from .cubic_uniform_spline_eval_funcs import cu_find_span, cu_basis_funs -from .accelerated_spline_interpolators import solve_system_periodic, solve_system_nonperiodic +from .accelerated_spline_interpolators import solve_system_periodic, solve_system_nonperiodic, solve_2d_system __all__ = ["SplineInterpolator1D", "SplineInterpolator2D"] @@ -234,6 +234,12 @@ def compute_interpolant(self, ug, spl): # assert basis1 is self._basis1 # assert basis2 is self._basis2 + if basis1.periodic and not basis2.periodic and self._interp1._splu: + solve_2d_system(ug, spl, self._bwork, self._interp2._bmat, self._interp2._l, + self._interp2._u, self._interp2._ipiv, self._interp1._offset, + self._interp1._splu) + return + n1, n2 = basis1.nbasis, basis2.nbasis p1, p2 = basis1.degree, basis2.degree assert ug.shape == (n1, n2) diff --git a/pygyro/splines/tests/test_spline_interpolators.py b/pygyro/splines/tests/test_spline_interpolators.py index 3b023324..fd5aab7d 100644 --- a/pygyro/splines/tests/test_spline_interpolators.py +++ b/pygyro/splines/tests/test_spline_interpolators.py @@ -283,3 +283,55 @@ def test_SplineInterpolator2D_cosine(ncells, degree, periodic1, periodic2): breaks1).max(), np.diff(breaks2).max(), deg1, deg2) assert max_norm_err < err_bound + +# =============================================================================== + + +@pytest.mark.parametrize("ncells", [10, 20, 40, 80, 160]) +@pytest.mark.parametrize("degree", range(1, 5)) +def test_SplineInterpolator2D_radial(ncells, degree): + """ + TODO + """ + periodic1 = True + periodic2 = False + + nc1 = nc2 = ncells + deg1 = deg2 = degree + + f = AnalyticalProfile1D_Cos() + poly = AnalyticalProfile1D_Poly(degree) + domain1 = f.domain + domain2 = poly.domain + + # Along x1 + breaks1 = random_grid(domain1, nc1, 0.0) + knots1 = make_knots(breaks1, deg1, periodic1) + basis1 = BSplines(knots1, deg1, periodic1, False) + + # Along x2 + breaks2 = random_grid(domain2, nc2, 0.0) + knots2 = make_knots(breaks2, deg2, periodic2) + basis2 = BSplines(knots2, deg2, periodic2, False) + + # 2D spline and interpolator on tensor-product space + spline = Spline2D(basis1, basis2) + interp = SplineInterpolator2D(basis1, basis2) + + x1g = basis1.greville + x2g = basis2.greville + ug = f.eval(x1g)[:, None] * poly.eval(x2g)[None, :] + + interp.compute_interpolant(ug, spline) + + x1t = np.linspace(*domain1, num=20) + x2t = np.linspace(*domain2, num=20) + vals = np.empty((20, 20)) + spline.eval_vector(x1t, x2t, vals) + err = vals - f.eval(x1t)[:, None] * poly.eval(x2t)[None, :] + + max_norm_err = np.max(abs(err)) + err_bound = spline_1d_error_bound( + f, max(np.diff(breaks1).max(), np.diff(breaks2).max()), degree) + + assert max_norm_err < err_bound