From aaf0a4da617b19290a26d45d60b928fa6b0dcfa1 Mon Sep 17 00:00:00 2001 From: Emily Bourne Date: Sun, 17 Aug 2025 19:08:26 +0200 Subject: [PATCH 01/30] Stop using TypeVar --- pygyro/poisson/poisson_tools.py | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/pygyro/poisson/poisson_tools.py b/pygyro/poisson/poisson_tools.py index 89cc6e32..d5309cf0 100644 --- a/pygyro/poisson/poisson_tools.py +++ b/pygyro/poisson/poisson_tools.py @@ -1,8 +1,9 @@ -from pyccel.decorators import template +from typing import TypeVar +T = TypeVar('T', 'complex128[:,:,:]', 'float[:,:,:]') -@template('T', ('complex128[:,:,:]', 'float[:,:,:]')) -def get_perturbed_rho(rho: 'T', feq: 'float[:,:]', grid: 'float[:,:,:,:]', + +def get_perturbed_rho(rho: T, feq: 'float[:,:]', grid: 'float[:,:,:,:]', quad_coeffs: 'float[:]'): """ Calculate: @@ -33,8 +34,7 @@ def get_perturbed_rho(rho: 'T', feq: 'float[:,:]', grid: 'float[:,:,:,:]', (grid[i, j, k, l] - feq[i, l]) -@template('T', ('complex128[:,:,:]', 'float[:,:,:]')) -def get_rho(rho: 'T', grid: 'float[:,:,:,:]', quad_coeffs: 'float[:]'): +def get_rho(rho: T, grid: 'float[:,:,:,:]', quad_coeffs: 'float[:]'): """ Calculate: rho = \\int f dv From 2ddb5bd22c001d2ca6d5fce3308fe3219a18838a Mon Sep 17 00:00:00 2001 From: Emily Bourne Date: Sun, 17 Aug 2025 19:22:55 +0200 Subject: [PATCH 02/30] Update flag --- Makefile | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Makefile b/Makefile index 2adcca06..50f7384c 100644 --- a/Makefile +++ b/Makefile @@ -42,7 +42,7 @@ SO_EXT := $(shell $(PYTHON) -c "import sysconfig; print(sysconfig.get_config_var ifeq ($(ACC), pycc) TOOL := pyccel - TOOL_FLAGS := --compiler=$(COMP) --flags ' $(FC_FLAGS)' --language=$(LANGUAGE) + TOOL_FLAGS := --compiler-family=$(COMP) --flags ' $(FC_FLAGS)' --language=$(LANGUAGE) NAME_PREFIX := else ifeq ($(ACC), numba) From 39528a4eade9b646740b027e622076cf2d4249aa Mon Sep 17 00:00:00 2001 From: Emily Bourne Date: Sun, 17 Aug 2025 18:48:13 +0200 Subject: [PATCH 03/30] Add final annotations --- .../cubic_uniform_spline_eval_funcs.py | 13 +++++++------ pygyro/splines/spline_eval_funcs.py | 19 ++++++++++--------- 2 files changed, 17 insertions(+), 15 deletions(-) diff --git a/pygyro/splines/cubic_uniform_spline_eval_funcs.py b/pygyro/splines/cubic_uniform_spline_eval_funcs.py index b7c4e9ca..6a1bb9c9 100644 --- a/pygyro/splines/cubic_uniform_spline_eval_funcs.py +++ b/pygyro/splines/cubic_uniform_spline_eval_funcs.py @@ -1,3 +1,4 @@ +from typing import Final from pyccel.decorators import pure, stack_array from numpy import empty @@ -131,7 +132,7 @@ def cu_basis_funs_1st_der(span: 'int', offset: 'float', dx: 'float', ders: 'floa @pure @stack_array('basis') -def cu_eval_spline_1d_scalar(x: 'float', knots: 'float[:]', degree: 'int', coeffs: 'float[:]', der: 'int') -> 'float': +def cu_eval_spline_1d_scalar(x: 'float', knots: 'Final[float[:]]', degree: 'int', coeffs: 'Final[float[:]]', der: 'int') -> 'float': """ TODO """ @@ -153,7 +154,7 @@ def cu_eval_spline_1d_scalar(x: 'float', knots: 'float[:]', degree: 'int', coeff @pure @stack_array('basis') -def cu_eval_spline_1d_vector(x: 'float[:]', knots: 'float[:]', degree: 'int', coeffs: 'float[:]', y: 'float[:]', der: 'int' = 0): +def cu_eval_spline_1d_vector(x: 'Final[float[:]]', knots: 'Final[float[:]]', degree: 'int', coeffs: 'Final[float[:]]', y: 'float[:]', der: 'int' = 0): """ TODO """ @@ -182,8 +183,8 @@ def cu_eval_spline_1d_vector(x: 'float[:]', knots: 'float[:]', degree: 'int', co @pure @stack_array('basis1', 'basis2', 'theCoeffs') -def cu_eval_spline_2d_scalar(x: 'float', y: 'float', kts1: 'float[:]', deg1: 'int', kts2: 'float[:]', deg2: 'int', - coeffs: 'float[:,:]', der1: 'int' = 0, der2: 'int' = 0) -> 'float': +def cu_eval_spline_2d_scalar(x: 'float', y: 'float', kts1: 'Final[float[:]]', deg1: 'int', kts2: 'Final[float[:]]', deg2: 'int', + coeffs: 'Final[float[:,:]]', der1: 'int' = 0, der2: 'int' = 0) -> 'float': """ TODO """ @@ -221,8 +222,8 @@ def cu_eval_spline_2d_scalar(x: 'float', y: 'float', kts1: 'float[:]', deg1: 'in @pure @stack_array('basis1', 'basis2', 'theCoeffs') -def cu_eval_spline_2d_cross(X: 'float[:]', Y: 'float[:]', kts1: 'float[:]', deg1: 'int', kts2: 'float[:]', deg2: 'int', - coeffs: 'float[:,:]', z: 'float[:,:]', der1: 'int' = 0, der2: 'int' = 0): +def cu_eval_spline_2d_cross(X: 'Final[float[:]]', Y: 'Final[float[:]]', kts1: 'Final[float[:]]', deg1: 'int', kts2: 'Final[float[:]]', deg2: 'int', + coeffs: 'Final[float[:,:]]', z: 'float[:,:]', der1: 'int' = 0, der2: 'int' = 0): """ TODO """ diff --git a/pygyro/splines/spline_eval_funcs.py b/pygyro/splines/spline_eval_funcs.py index 2b986714..11653cdb 100644 --- a/pygyro/splines/spline_eval_funcs.py +++ b/pygyro/splines/spline_eval_funcs.py @@ -1,9 +1,10 @@ +from typing import Final from pyccel.decorators import pure, stack_array from numpy import empty @pure -def nu_find_span(knots: 'float[:]', degree: 'int', x: 'float') -> int: +def nu_find_span(knots: 'Final[float[:]]', degree: 'int', x: 'float') -> int: """ Determine the knot span index at location x, given the B-Splines' knot sequence and polynomial degree. See @@ -60,7 +61,7 @@ def nu_find_span(knots: 'float[:]', degree: 'int', x: 'float') -> int: @pure @stack_array('left', 'right') -def nu_basis_funs(knots: 'float[:]', degree: 'int', x: 'float', span: 'int', values: 'float[:]'): +def nu_basis_funs(knots: 'Final[float[:]]', degree: 'int', x: 'float', span: 'int', values: 'float[:]'): """ Compute the non-vanishing B-splines at location x, given the knot sequence, polynomial degree and knot @@ -112,7 +113,7 @@ def nu_basis_funs(knots: 'float[:]', degree: 'int', x: 'float', span: 'int', val @pure @stack_array('values') -def nu_basis_funs_1st_der(knots: 'float[:]', degree: 'int', x: 'float', span: 'int', ders: 'float[:]'): +def nu_basis_funs_1st_der(knots: 'Final[float[:]]', degree: 'int', x: 'float', span: 'int', ders: 'float[:]'): """ Compute the first derivative of the non-vanishing B-splines at location x, given the knot sequence, polynomial degree @@ -166,7 +167,7 @@ def nu_basis_funs_1st_der(knots: 'float[:]', degree: 'int', x: 'float', span: 'i @pure @stack_array('basis') -def nu_eval_spline_1d_scalar(x: 'float', knots: 'float[:]', degree: 'int', coeffs: 'float[:]', der: 'int') -> 'float': +def nu_eval_spline_1d_scalar(x: 'float', knots: 'Final[float[:]]', degree: 'int', coeffs: 'Final[float[:]]', der: 'int') -> 'float': """ TODO """ @@ -186,7 +187,7 @@ def nu_eval_spline_1d_scalar(x: 'float', knots: 'float[:]', degree: 'int', coeff @pure @stack_array('basis') -def nu_eval_spline_1d_vector(x: 'float[:]', knots: 'float[:]', degree: 'int', coeffs: 'float[:]', y: 'float[:]', der: 'int' = 0): +def nu_eval_spline_1d_vector(x: 'Final[float[:]]', knots: 'Final[float[:]]', degree: 'int', coeffs: 'Final[float[:]]', y: 'float[:]', der: 'int' = 0): """ TODO """ @@ -213,8 +214,8 @@ def nu_eval_spline_1d_vector(x: 'float[:]', knots: 'float[:]', degree: 'int', co @pure @stack_array('basis1', 'basis2', 'theCoeffs') -def nu_eval_spline_2d_scalar(x: 'float', y: 'float', kts1: 'float[:]', deg1: 'int', kts2: 'float[:]', deg2: 'int', - coeffs: 'float[:,:]', der1: 'int' = 0, der2: 'int' = 0) -> 'float': +def nu_eval_spline_2d_scalar(x: 'float', y: 'float', kts1: 'Final[float[:]]', deg1: 'int', kts2: 'Final[float[:]]', deg2: 'int', + coeffs: 'Final[float[:,:]]', der1: 'int' = 0, der2: 'int' = 0) -> 'float': """ TODO """ @@ -247,8 +248,8 @@ def nu_eval_spline_2d_scalar(x: 'float', y: 'float', kts1: 'float[:]', deg1: 'in @pure @stack_array('basis1', 'basis2', 'theCoeffs') -def nu_eval_spline_2d_cross(X: 'float[:]', Y: 'float[:]', kts1: 'float[:]', deg1: 'int', kts2: 'float[:]', deg2: 'int', - coeffs: 'float[:,:]', z: 'float[:,:]', der1: 'int' = 0, der2: 'int' = 0): +def nu_eval_spline_2d_cross(X: 'Final[float[:]]', Y: 'Final[float[:]]', kts1: 'Final[float[:]]', deg1: 'int', kts2: 'Final[float[:]]', deg2: 'int', + coeffs: 'Final[float[:,:]]', z: 'float[:,:]', der1: 'int' = 0, der2: 'int' = 0): """ TODO """ From 4fbc9c3360a01a957f2754dd07c2ce14b38946d4 Mon Sep 17 00:00:00 2001 From: Emily Bourne Date: Sun, 17 Aug 2025 19:27:15 +0200 Subject: [PATCH 04/30] Use typing.Final --- pygyro/advection/accelerated_advection_steps.py | 15 ++++++++------- 1 file changed, 8 insertions(+), 7 deletions(-) diff --git a/pygyro/advection/accelerated_advection_steps.py b/pygyro/advection/accelerated_advection_steps.py index b31aea2e..d100b2dd 100644 --- a/pygyro/advection/accelerated_advection_steps.py +++ b/pygyro/advection/accelerated_advection_steps.py @@ -1,3 +1,4 @@ +from typing import Final from pyccel.decorators import pure from ..splines.spline_eval_funcs import nu_eval_spline_1d_scalar, nu_eval_spline_1d_vector from ..splines.spline_eval_funcs import nu_eval_spline_2d_cross, nu_eval_spline_2d_scalar @@ -16,8 +17,8 @@ def general_poloidal_advection_step_expl(f: 'float[:,:]', kts1Pol: 'float[:]', kts2Pol: 'float[:]', coeffsPol: 'float[:,:]', deg1Pol: 'int', deg2Pol: 'int', CN0: 'float', kN0: 'float', deltaRN0: 'float', rp: 'float', CTi: 'float', kTi: 'float', deltaRTi: 'float', B0: 'float', nulBound: 'bool', - eval_spline_2d_cross: '()(const float[:], const float[:], const float[:], int, const float[:], int, const float[:,:], float[:,:], int, int)', - eval_spline_2d_scalar: '(float)(float, float, const float[:], int, const float[:], int, const float[:,:], int, int)'): + eval_spline_2d_cross: '()(Final[float[:]], Final[float[:]], Final[float[:]], int, Final[float[:]], int, Final[float[:,:]], float[:,:], int, int)', + eval_spline_2d_scalar: '(float)(float, float, Final[float[:]], int, Final[float[:]], int, Final[float[:,:]], int, int)'): """ Carry out an advection step for the poloidal advection @@ -147,7 +148,7 @@ def general_v_parallel_advection_eval_step(f: 'float[:]', vPts: 'float[:]', coeffs: 'float[:]', CN0: 'float', kN0: 'float', deltaRN0: 'float', rp: 'float', CTi: 'float', kTi: 'float', deltaRTi: 'float', bound: 'int', - eval_spline_1d_scalar: '(float)(float, const float[:], int, const float[:], int)'): + eval_spline_1d_scalar: '(float)(float, Final[float[:]], int, Final[float[:]], int)'): """ TODO """ @@ -193,8 +194,8 @@ def general_get_lagrange_vals(i: 'int', shifts: 'int[:]', vals: 'float[:,:,:]', qVals: 'float[:]', thetaShifts: 'float[:]', kts: 'float[:]', deg: 'int', coeffs: 'float[:]', - eval_spline_1d_vector: '()(const float[:], const float[:], int, const float[:], float[:], int)', - eval_spline_1d_scalar: '(float)(float, const float[:], int, const float[:], int)'): + eval_spline_1d_vector: '()(Final[float[:]], Final[float[:]], int, Final[float[:]], float[:], int)', + eval_spline_1d_scalar: '(float)(float, Final[float[:]], int, Final[float[:]], int)'): """ TODO """ @@ -242,8 +243,8 @@ def general_poloidal_advection_step_impl(f: 'float[:,:]', dt: 'float', v: 'float kts1Pol: 'float[:]', kts2Pol: 'float[:]', coeffsPol: 'float[:,:]', deg1Pol: 'int', deg2Pol: 'int', CN0: 'float', kN0: 'float', deltaRN0: 'float', rp: 'float', CTi: 'float', kTi: 'float', deltaRTi: 'float', B0: 'float', tol: 'float', nulBound: 'bool', - eval_spline_2d_cross: '()(const float[:], const float[:], const float[:], int, const float[:], int, const float[:,:], float[:,:], int, int)', - eval_spline_2d_scalar: '(float)(float, float, const float[:], int, const float[:], int, const float[:,:], int, int)'): + eval_spline_2d_cross: '()(Final[float[:]], Final[float[:]], Final[float[:]], int, Final[float[:]], int, Final[float[:,:]], float[:,:], int, int)', + eval_spline_2d_scalar: '(float)(float, float, Final[float[:]], int, Final[float[:]], int, Final[float[:,:]], int, int)'): """ Carry out an advection step for the poloidal advection From 5b4fb434110cd2d4971d8030441b8d7a241a5cc1 Mon Sep 17 00:00:00 2001 From: Emily Bourne Date: Mon, 18 Aug 2025 10:19:07 +0200 Subject: [PATCH 05/30] Update requirements --- requirements.txt | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/requirements.txt b/requirements.txt index 3e75544f..801c56a9 100644 --- a/requirements.txt +++ b/requirements.txt @@ -4,7 +4,7 @@ scipy>=1.1.0 pytest>=2.8.7 mpi4py>=3.0.0 h5py>=2.8.0 -pyccel>=1.2.3 +pyccel>=2.0.0 # h5py must be built from source using MPI compiler # and linked to parallel HDF5 library. To do so set From 789d9c0fc211c7e1414bb12de180392f726d6d3a Mon Sep 17 00:00:00 2001 From: Emily Bourne Date: Mon, 18 Aug 2025 10:36:01 +0200 Subject: [PATCH 06/30] Try to fix install --- .github/actions/pip_installation/action.yml | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/.github/actions/pip_installation/action.yml b/.github/actions/pip_installation/action.yml index 729f15f2..7e793b7c 100644 --- a/.github/actions/pip_installation/action.yml +++ b/.github/actions/pip_installation/action.yml @@ -15,8 +15,7 @@ runs: shell: bash run: | export CC="mpicc" HDF5_MPI="ON" - python -m pip install Cython numpy>=1.15.1 pkgconfig mpi4py>=3.0.0 # h5py requirements - python -m pip install --no-build-isolation -r requirements.txt + python -m pip install -r requirements.txt python -m pip install numba python -m pip install pythran From 1a3b3403d74fc93824eb3e666936ccaf6b04b232 Mon Sep 17 00:00:00 2001 From: Emily Bourne Date: Mon, 18 Aug 2025 10:51:09 +0200 Subject: [PATCH 07/30] Scipy's trapz function has been renamed --- pygyro/advection/convergence_test_advection.py | 16 ++++++++-------- pygyro/advection/test_advection.py | 8 ++++---- .../poisson/convergence_test_poisson_solver.py | 8 ++++---- pygyro/poisson/test_poisson_solver.py | 4 ++-- 4 files changed, 18 insertions(+), 18 deletions(-) diff --git a/pygyro/advection/convergence_test_advection.py b/pygyro/advection/convergence_test_advection.py index 041b09b2..78239998 100644 --- a/pygyro/advection/convergence_test_advection.py +++ b/pygyro/advection/convergence_test_advection.py @@ -1,6 +1,6 @@ import pytest import numpy as np -from scipy.integrate import trapz +from scipy.integrate import trapezoid from numpy import pi from .advection import FluxSurfaceAdvection, PoloidalAdvection, VParallelAdvection, ParallelGradient, Layout @@ -60,7 +60,7 @@ def test_vParallelAdvection(): else: fEnd[i] = gaussLike(x[i]-c*dt*N) - l2[j] = np.sqrt(trapz((f-fEnd).flatten()**2, x)) + l2[j] = np.sqrt(trapezoid((f-fEnd).flatten()**2, x)) linf[j] = np.linalg.norm((f-fEnd).flatten(), np.inf) print("l2:", l2) @@ -210,7 +210,7 @@ def test_poloidalAdvection_constantAdv(initConditions, xc, yc): polAdv.exact_step(f_vals[:, :], endPts, v) linf[i] = np.linalg.norm((f_vals-final_f_vals).flatten(), np.inf) - l2[i] = np.sqrt(trapz(trapz((f_vals-final_f_vals)**2, + l2[i] = np.sqrt(trapezoid(trapezoid((f_vals-final_f_vals)**2, eta_grids[1], axis=0)*eta_grids[0], eta_grids[0])) print("l2:", l2[i]) @@ -329,7 +329,7 @@ def test_poloidalAdvection_constantAdv_dt(): polAdv.step(f_vals[:, :], dt, phi, v) linf[i] = np.linalg.norm((f_vals-final_f_vals).flatten(), np.inf) - l2[i] = np.sqrt(trapz(trapz((f_vals-final_f_vals)**2, + l2[i] = np.sqrt(trapezoid(trapezoid((f_vals-final_f_vals)**2, eta_grids[1], axis=0)*eta_grids[0], eta_grids[0])) print(N, "l2:", l2[i]) @@ -448,7 +448,7 @@ def test_fluxAdvection(): linf[i] = np.linalg.norm((f_vals-final_f_vals).flatten(), np.inf) l2[i] = np.sqrt( - trapz(trapz((f_vals-final_f_vals)**2, eta_vals[1], axis=0), eta_vals[2])) + trapezoid(trapezoid((f_vals-final_f_vals)**2, eta_vals[1], axis=0), eta_vals[2])) print(N, "l2:", l2[i]) print(N, "linf:", linf[i]) @@ -561,7 +561,7 @@ def test_fluxAdvectionAligned(): linf[i] = np.linalg.norm((f_vals-final_f_vals).flatten(), np.inf) l2[i] = np.sqrt( - trapz(trapz((f_vals-final_f_vals)**2, eta_vals[1], axis=0), eta_vals[2])) + trapezoid(trapezoid((f_vals-final_f_vals)**2, eta_vals[1], axis=0), eta_vals[2])) print(N, "l2:", l2[i]) print(N, "linf:", linf[i]) @@ -661,7 +661,7 @@ def test_Phi_deriv_dtheta(): err = approxGrad-exactGrad - l2[i] = np.sqrt(np.trapz(np.trapz(err**2, dx=dz), dx=dtheta)) + l2[i] = np.sqrt(np.trapezoid(np.trapezoid(err**2, dx=dz), dx=dtheta)) linf[i] = np.linalg.norm(err.flatten(), np.inf) npts[1] *= 2 @@ -744,7 +744,7 @@ def test_Phi_deriv_dz(): err = approxGrad-exactGrad - l2[i] = np.sqrt(np.trapz(np.trapz(err**2, dx=dz), dx=dtheta)) + l2[i] = np.sqrt(np.trapezoid(np.trapezoid(err**2, dx=dz), dx=dtheta)) linf[i] = np.linalg.norm(err.flatten(), np.inf) npts[2] *= 2 diff --git a/pygyro/advection/test_advection.py b/pygyro/advection/test_advection.py index 596dfe3b..57a52990 100644 --- a/pygyro/advection/test_advection.py +++ b/pygyro/advection/test_advection.py @@ -1,6 +1,6 @@ from mpi4py import MPI import pytest -from scipy.integrate import trapz +from scipy.integrate import trapezoid import numpy as np from ..initialisation.setups import setupCylindricalGrid @@ -257,7 +257,7 @@ def test_poloidalAdvection(dt, v, xc, yc): finalPts[1][:] = np.sqrt(x * x + y * y) final_f_vals[:, :] = initConds(finalPts[1], finalPts[0])+fEdge - l2 = np.sqrt(trapz(trapz((f_vals-final_f_vals)**2, + l2 = np.sqrt(trapezoid(trapezoid((f_vals-final_f_vals)**2, eta_grids[1], axis=0)*eta_grids[0], eta_grids[0])) assert l2 < 0.2 @@ -320,7 +320,7 @@ def test_poloidalAdvectionImplicit(dt, v, xc, yc): finalPts[1][:] = np.sqrt(x * x + y * y) final_f_vals[:, :] = initConds(finalPts[1], finalPts[0]) - l2 = np.sqrt(trapz(trapz((f_vals-final_f_vals)**2, + l2 = np.sqrt(trapezoid(trapezoid((f_vals-final_f_vals)**2, eta_grids[1], axis=0)*eta_grids[0], eta_grids[0])) assert l2 < 0.2 @@ -624,7 +624,7 @@ def test_Phi_deriv_dz(phiOrder, zOrder): err = approxGrad-exactGrad - l2[i] = np.sqrt(np.trapz(np.trapz(err**2, dx=dz), dx=dtheta)) + l2[i] = np.sqrt(np.trapezoid(np.trapezoid(err**2, dx=dz), dx=dtheta)) linf[i] = np.linalg.norm(err.flatten(), np.inf) npts[1] *= 2 diff --git a/pygyro/poisson/convergence_test_poisson_solver.py b/pygyro/poisson/convergence_test_poisson_solver.py index 8e2eeb63..2811bb6a 100644 --- a/pygyro/poisson/convergence_test_poisson_solver.py +++ b/pygyro/poisson/convergence_test_poisson_solver.py @@ -2,7 +2,7 @@ import numpy as np import pytest from math import pi -from scipy.integrate import trapz +from scipy.integrate import trapezoid from numpy.polynomial.legendre import leggauss import matplotlib.pyplot as plt @@ -753,7 +753,7 @@ def test_ddTheta(deg): # ~ phi._f=phi._f*2 err = (phi._f-phi_exact._f)[0, 0] - l2[c] = np.sqrt(trapz(np.real(err)**2, q)) + l2[c] = np.sqrt(trapezoid(np.real(err)**2, q)) lInf[c] = np.max(np.abs(np.real(phi._f-phi_exact._f))) npts[1] *= 2 @@ -912,7 +912,7 @@ def test_QuasiNeutralityEquation_pointConverge(): # ~ l2Q[i]=np.sum((approxSpline.eval(evalPts)-np.cos(rArg)**4)**2 \ # ~ * multFactor*weights) - l2[c] = np.sqrt(trapz(l2Q, eta_grid[1])) + l2[c] = np.sqrt(trapezoid(l2Q, eta_grid[1])) lInf[c] = lI npts[0] *= 2 @@ -1096,7 +1096,7 @@ def test_QuasiNeutralityEquation_degreeConverge(): l2Q[i] = np.sum((approxSpline.eval(evalPts)-np.cos(rArg)**4)**2 * multFactor*weights) - l2[c] = np.sqrt(trapz(l2Q, eta_grid[1])) + l2[c] = np.sqrt(trapezoid(l2Q, eta_grid[1])) lInf[c] = lI degree[0] *= 2 diff --git a/pygyro/poisson/test_poisson_solver.py b/pygyro/poisson/test_poisson_solver.py index 7b6b184d..6daf99f6 100644 --- a/pygyro/poisson/test_poisson_solver.py +++ b/pygyro/poisson/test_poisson_solver.py @@ -2,7 +2,7 @@ import numpy as np import pytest from math import pi -from scipy.integrate import trapz +from scipy.integrate import trapezoid from ..model.process_grid import compute_2d_process_grid from ..model.layout import LayoutSwapper, getLayoutHandler @@ -1258,7 +1258,7 @@ def test_BasicPoissonEquation_exact(deg): x = eta_grid[0] err = (phi._f-phi_exact._f)[0, 0] - l2 = np.sqrt(trapz(np.real(err*err.conj()), x)) + l2 = np.sqrt(trapezoid(np.real(err*err.conj()), x)) lInf = np.max(np.abs(phi._f-phi_exact._f)) assert l2 < 1e-10 From 022be0d10737f0e39f6aad38ce479989f2ed93d3 Mon Sep 17 00:00:00 2001 From: Emily Bourne Date: Tue, 19 Aug 2025 10:44:41 +0200 Subject: [PATCH 08/30] Ensure Final is used somewhere recognisable --- pygyro/advection/accelerated_advection_steps.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/pygyro/advection/accelerated_advection_steps.py b/pygyro/advection/accelerated_advection_steps.py index d100b2dd..e1939581 100644 --- a/pygyro/advection/accelerated_advection_steps.py +++ b/pygyro/advection/accelerated_advection_steps.py @@ -9,7 +9,7 @@ def general_poloidal_advection_step_expl(f: 'float[:,:]', dt: 'float', v: 'float', - rPts: 'float[:]', qPts: 'float[:]', + rPts: 'Final[float[:]]', qPts: 'Final[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[:,:]', From e80c7ebf22f7cb68e69b35ce9ff7cb75503a76f1 Mon Sep 17 00:00:00 2001 From: Emily Bourne Date: Wed, 20 Aug 2025 09:20:16 +0200 Subject: [PATCH 09/30] Prefer assert to raise --- pygyro/splines/splines.py | 10 ++-------- 1 file changed, 2 insertions(+), 8 deletions(-) diff --git a/pygyro/splines/splines.py b/pygyro/splines/splines.py index 980046ed..d2ce2609 100644 --- a/pygyro/splines/splines.py +++ b/pygyro/splines/splines.py @@ -369,14 +369,8 @@ def __init__(self, basis1, basis2): self._basis2 = basis2 self._coeffs = np.zeros(shape) - if basis1.degree > 5: - raise NotImplementedError( - "scipy.interpolate.bisplev needs p1 <= 5") - - if basis2.degree > 5: - raise NotImplementedError( - "scipy.interpolate.bisplev needs p2 <= 5") - + assert basis1.degree <= 5 + assert basis2.degree <= 5 assert basis1.cubic_uniform == basis2.cubic_uniform @property From c095a4a729c9f691c153e6eedaca50eaa4854663 Mon Sep 17 00:00:00 2001 From: Emily Bourne Date: Wed, 20 Aug 2025 09:49:45 +0200 Subject: [PATCH 10/30] Remove unused method find_cell --- pygyro/splines/splines.py | 9 --------- 1 file changed, 9 deletions(-) diff --git a/pygyro/splines/splines.py b/pygyro/splines/splines.py index d2ce2609..bf999c9a 100644 --- a/pygyro/splines/splines.py +++ b/pygyro/splines/splines.py @@ -226,15 +226,6 @@ def __getitem__(self, i): spl.coeffs[n:n+p] = spl.coeffs[0:p] return spl - # ... - def find_cell(self, x): - """ Index i of cell $C_{i} := [x_{i},x_{i+1})$ that contains point x. - Last cell includes right endpoint. - """ - a, b = self.domain - assert a <= x <= b - return int(np.searchsorted(self.breaks, x, side='right') - 1) - def _build_integrals(self): n = self.nbasis d = self.degree From 6751160acf4107a8be968ac881954ec841716193 Mon Sep 17 00:00:00 2001 From: Emily Bourne Date: Wed, 20 Aug 2025 11:10:55 +0200 Subject: [PATCH 11/30] Make Spline1DComplex class to avoid Interface constructor --- pygyro/poisson/poisson_solver.py | 2 +- pygyro/splines/spline_interpolators.py | 15 ++++-- pygyro/splines/splines.py | 67 ++++++++++++++++++++++++-- 3 files changed, 76 insertions(+), 8 deletions(-) diff --git a/pygyro/poisson/poisson_solver.py b/pygyro/poisson/poisson_solver.py index 64ded8ef..396d6b96 100644 --- a/pygyro/poisson/poisson_solver.py +++ b/pygyro/poisson/poisson_solver.py @@ -299,7 +299,7 @@ def __init__(self, degree: int, rspline: BSplines, nr: int, nTheta: int, # Create the tools required for the interpolation self._interpolator = SplineInterpolator1D(self._rspline, dtype=complex) - self._spline = Spline1D(self._rspline, np.complex128) + self._spline = Spline1DComplex(self._rspline) self._real_spline = Spline1D(self._rspline) self._realMem = np.empty(nr) diff --git a/pygyro/splines/spline_interpolators.py b/pygyro/splines/spline_interpolators.py index 17d5ef2e..f67233fb 100644 --- a/pygyro/splines/spline_interpolators.py +++ b/pygyro/splines/spline_interpolators.py @@ -6,7 +6,7 @@ from scipy.sparse import csr_matrix, csc_matrix, dia_matrix from scipy.sparse.linalg import splu -from .splines import BSplines, Spline1D, Spline2D +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 @@ -64,7 +64,7 @@ def compute_interpolant(self, ug, spl): The spline in which the coefficients will be saved """ - assert isinstance(spl, Spline1D) + assert isinstance(spl, (Spline1D, Spline1DComplex)) assert spl.basis is self._basis assert len(ug) == self._basis.nbasis @@ -192,8 +192,15 @@ def __init__(self, basis1, basis2, dtype=float): self._basis1 = basis1 self._basis2 = basis2 - self._spline1 = Spline1D(basis1, dtype) - self._spline2 = Spline1D(basis2, dtype) + + if dtype is float: + self._spline1 = Spline1D(basis1) + self._spline2 = Spline1D(basis2) + else: + assert dtype is np.complex128 + self._spline1 = Spline1DComplex(basis1) + self._spline2 = Spline1DComplex(basis2) + self._interp1 = SplineInterpolator1D(basis1, dtype) self._interp2 = SplineInterpolator1D(basis2, dtype) diff --git a/pygyro/splines/splines.py b/pygyro/splines/splines.py index bf999c9a..75470ace 100644 --- a/pygyro/splines/splines.py +++ b/pygyro/splines/splines.py @@ -9,7 +9,7 @@ from .cubic_uniform_spline_eval_funcs import cu_eval_spline_1d_scalar, cu_eval_spline_1d_vector from .cubic_uniform_spline_eval_funcs import cu_eval_spline_2d_cross, cu_eval_spline_2d_scalar -__all__ = ['make_knots', 'BSplines', 'Spline1D', 'Spline2D'] +__all__ = ['make_knots', 'BSplines', 'Spline1D', 'Spline1DComplex', 'Spline2D'] # =============================================================================== @@ -288,10 +288,71 @@ class Spline1D(): TODO """ - def __init__(self, basis, dtype=float): + def __init__(self, basis): assert isinstance(basis, BSplines) self._basis = basis - self._coeffs = np.zeros(basis.ncells + basis.degree, dtype=dtype) + self._coeffs = np.zeros(basis.ncells + basis.degree, dtype=float) + + @property + def basis(self): + """ + TODO + """ + return self._basis + + @property + def coeffs(self): + """ + TODO + """ + return self._coeffs + + def eval(self, x, der=0): + """ + TODO + """ + if (hasattr(x, '__len__')): + result = np.empty_like(x) + if self._basis.cubic_uniform: + cu_eval_spline_1d_vector(x, self._basis.knots, + self._basis.degree, self._coeffs, result, der) + else: + nu_eval_spline_1d_vector(x, self._basis.knots, + self._basis.degree, self._coeffs, result, der) + else: + if self._basis.cubic_uniform: + result = cu_eval_spline_1d_scalar( + x, self._basis.knots, self._basis.degree, self._coeffs, der) + else: + result = nu_eval_spline_1d_scalar( + x, self._basis.knots, self._basis.degree, self._coeffs, der) + return result + + """ + tck = (self._basis.knots, self._coeffs, self._basis.degree) + return splev( x, tck, der ) + """ + + def eval_vector(self, x, y, der=0): + """ + TODO + """ + if self._basis.cubic_uniform: + cu_eval_spline_1d_vector(x, self._basis.knots, + self._basis.degree, self._coeffs, y, der) + else: + nu_eval_spline_1d_vector(x, self._basis.knots, + self._basis.degree, self._coeffs, y, der) + +class Spline1DComplex(): + """ + TODO + """ + + def __init__(self, basis): + assert isinstance(basis, BSplines) + self._basis = basis + self._coeffs = np.zeros(basis.ncells + basis.degree, dtype=np.complex128) @property def basis(self): From fddba4c4523b00c9b1e21e24a6d8c66ece47be75 Mon Sep 17 00:00:00 2001 From: Emily Bourne Date: Wed, 20 Aug 2025 22:10:14 +0200 Subject: [PATCH 12/30] Comment assertion test --- pygyro/splines/tests/test_splines.py | 18 +++++++++--------- 1 file changed, 9 insertions(+), 9 deletions(-) diff --git a/pygyro/splines/tests/test_splines.py b/pygyro/splines/tests/test_splines.py index 313754f7..979ca75f 100644 --- a/pygyro/splines/tests/test_splines.py +++ b/pygyro/splines/tests/test_splines.py @@ -31,15 +31,15 @@ def test_make_knots_periodic(ncells, degree): # =============================================================================== -@pytest.mark.serial -@pytest.mark.parametrize("ncells", [1, 5, 10, 23]) -def test_make_knots_periodic_should_fail(ncells): - """ - TODO - """ - breaks = np.arange(ncells+1, dtype=float) - with pytest.raises(AssertionError): - _ = make_knots(breaks, degree=ncells+1, periodic=True) +#@pytest.mark.serial +#@pytest.mark.parametrize("ncells", [1, 5, 10, 23]) +#def test_make_knots_periodic_should_fail(ncells): +# """ +# TODO +# """ +# breaks = np.arange(ncells+1, dtype=float) +# with pytest.raises(AssertionError): +# _ = make_knots(breaks, degree=ncells+1, periodic=True) # =============================================================================== From f54b9263ded13a9b163469400c197063882d1c10 Mon Sep 17 00:00:00 2001 From: Emily Bourne Date: Wed, 20 Aug 2025 22:10:38 +0200 Subject: [PATCH 13/30] Add type annotations. Allow for complex splines --- .../splines/cubic_uniform_spline_eval_funcs.py | 17 +++++++++-------- pygyro/splines/spline_eval_funcs.py | 17 +++++++++-------- 2 files changed, 18 insertions(+), 16 deletions(-) diff --git a/pygyro/splines/cubic_uniform_spline_eval_funcs.py b/pygyro/splines/cubic_uniform_spline_eval_funcs.py index 6a1bb9c9..66d6eb26 100644 --- a/pygyro/splines/cubic_uniform_spline_eval_funcs.py +++ b/pygyro/splines/cubic_uniform_spline_eval_funcs.py @@ -1,7 +1,8 @@ -from typing import Final +from typing import Final, TypeVar from pyccel.decorators import pure, stack_array from numpy import empty +CoeffType = TypeVar('CoeffType', float, complex) @pure def cu_find_span(xmin: 'float', xmax: 'float', dx: 'float', x: 'float', ncells: 'int'): @@ -132,7 +133,7 @@ def cu_basis_funs_1st_der(span: 'int', offset: 'float', dx: 'float', ders: 'floa @pure @stack_array('basis') -def cu_eval_spline_1d_scalar(x: 'float', knots: 'Final[float[:]]', degree: 'int', coeffs: 'Final[float[:]]', der: 'int') -> 'float': +def cu_eval_spline_1d_scalar(x: 'float', knots: 'Final[float[:]]', degree: 'int', coeffs: 'Final[CoeffType[:]]', der: 'int') -> 'CoeffType': """ TODO """ @@ -146,7 +147,7 @@ def cu_eval_spline_1d_scalar(x: 'float', knots: 'Final[float[:]]', degree: 'int' elif (der == 1): cu_basis_funs_1st_der(span, offset, dx, basis) - y = 0.0 + y = 0.0*coeffs[0] for j in range(4): y += coeffs[span-3+j]*basis[j] return y @@ -154,7 +155,7 @@ def cu_eval_spline_1d_scalar(x: 'float', knots: 'Final[float[:]]', degree: 'int' @pure @stack_array('basis') -def cu_eval_spline_1d_vector(x: 'Final[float[:]]', knots: 'Final[float[:]]', degree: 'int', coeffs: 'Final[float[:]]', y: 'float[:]', der: 'int' = 0): +def cu_eval_spline_1d_vector(x: 'Final[float[:]]', knots: 'Final[float[:]]', degree: 'int', coeffs: 'Final[CoeffType[:]]', y: 'CoeffType[:]', der: 'int' = 0): """ TODO """ @@ -184,7 +185,7 @@ def cu_eval_spline_1d_vector(x: 'Final[float[:]]', knots: 'Final[float[:]]', deg @pure @stack_array('basis1', 'basis2', 'theCoeffs') def cu_eval_spline_2d_scalar(x: 'float', y: 'float', kts1: 'Final[float[:]]', deg1: 'int', kts2: 'Final[float[:]]', deg2: 'int', - coeffs: 'Final[float[:,:]]', der1: 'int' = 0, der2: 'int' = 0) -> 'float': + coeffs: 'Final[CoeffType[:,:]]', der1: 'int' = 0, der2: 'int' = 0) -> 'CoeffType': """ TODO """ @@ -211,7 +212,7 @@ def cu_eval_spline_2d_scalar(x: 'float', y: 'float', kts1: 'Final[float[:]]', de theCoeffs = empty((4, 4)) theCoeffs[:, :] = coeffs[span1-deg1:span1+1, span2-deg2:span2+1] - z = 0.0 + z = 0.0*coeffs[0,0] for i in range(4): theCoeffs[i, 0] = theCoeffs[i, 0]*basis2[0] for j in range(1, 4): @@ -223,7 +224,7 @@ def cu_eval_spline_2d_scalar(x: 'float', y: 'float', kts1: 'Final[float[:]]', de @pure @stack_array('basis1', 'basis2', 'theCoeffs') def cu_eval_spline_2d_cross(X: 'Final[float[:]]', Y: 'Final[float[:]]', kts1: 'Final[float[:]]', deg1: 'int', kts2: 'Final[float[:]]', deg2: 'int', - coeffs: 'Final[float[:,:]]', z: 'float[:,:]', der1: 'int' = 0, der2: 'int' = 0): + coeffs: 'Final[CoeffType[:,:]]', z: 'CoeffType[:,:]', der1: 'int' = 0, der2: 'int' = 0): """ TODO """ @@ -315,7 +316,7 @@ def cu_eval_spline_2d_cross(X: 'Final[float[:]]', Y: 'Final[float[:]]', kts1: 'F @pure @stack_array('basis1', 'basis2', 'theCoeffs') def cu_eval_spline_2d_vector(x: 'float[:]', y: 'float[:]', kts1: 'float[:]', deg1: 'int', kts2: 'float[:]', deg2: 'int', - coeffs: 'float[:,:]', z: 'float[:]', der1: 'int' = 0, der2: 'int' = 0): + coeffs: 'CoeffType[:,:]', z: 'CoeffType[:]', der1: 'int' = 0, der2: 'int' = 0): """ TODO """ diff --git a/pygyro/splines/spline_eval_funcs.py b/pygyro/splines/spline_eval_funcs.py index 11653cdb..a4864b9d 100644 --- a/pygyro/splines/spline_eval_funcs.py +++ b/pygyro/splines/spline_eval_funcs.py @@ -1,7 +1,8 @@ -from typing import Final +from typing import Final, TypeVar from pyccel.decorators import pure, stack_array from numpy import empty +CoeffType = TypeVar('CoeffType', float, complex) @pure def nu_find_span(knots: 'Final[float[:]]', degree: 'int', x: 'float') -> int: @@ -167,7 +168,7 @@ def nu_basis_funs_1st_der(knots: 'Final[float[:]]', degree: 'int', x: 'float', s @pure @stack_array('basis') -def nu_eval_spline_1d_scalar(x: 'float', knots: 'Final[float[:]]', degree: 'int', coeffs: 'Final[float[:]]', der: 'int') -> 'float': +def nu_eval_spline_1d_scalar(x: 'float', knots: 'Final[float[:]]', degree: 'int', coeffs: 'Final[CoeffType[:]]', der: 'int') -> 'CoeffType': """ TODO """ @@ -179,7 +180,7 @@ def nu_eval_spline_1d_scalar(x: 'float', knots: 'Final[float[:]]', degree: 'int' elif (der == 1): nu_basis_funs_1st_der(knots, degree, x, span, basis) - y = 0.0 + y = 0.0*coeffs[0] for j in range(degree+1): y += coeffs[span-degree+j]*basis[j] return y @@ -187,7 +188,7 @@ def nu_eval_spline_1d_scalar(x: 'float', knots: 'Final[float[:]]', degree: 'int' @pure @stack_array('basis') -def nu_eval_spline_1d_vector(x: 'Final[float[:]]', knots: 'Final[float[:]]', degree: 'int', coeffs: 'Final[float[:]]', y: 'float[:]', der: 'int' = 0): +def nu_eval_spline_1d_vector(x: 'Final[float[:]]', knots: 'Final[float[:]]', degree: 'int', coeffs: 'Final[CoeffType[:]]', y: 'CoeffType[:]', der: 'int' = 0): """ TODO """ @@ -215,7 +216,7 @@ def nu_eval_spline_1d_vector(x: 'Final[float[:]]', knots: 'Final[float[:]]', deg @pure @stack_array('basis1', 'basis2', 'theCoeffs') def nu_eval_spline_2d_scalar(x: 'float', y: 'float', kts1: 'Final[float[:]]', deg1: 'int', kts2: 'Final[float[:]]', deg2: 'int', - coeffs: 'Final[float[:,:]]', der1: 'int' = 0, der2: 'int' = 0) -> 'float': + coeffs: 'Final[CoeffType[:,:]]', der1: 'int' = 0, der2: 'int' = 0) -> 'CoeffType': """ TODO """ @@ -237,7 +238,7 @@ def nu_eval_spline_2d_scalar(x: 'float', y: 'float', kts1: 'Final[float[:]]', de theCoeffs = empty((deg1+1, deg2+1)) theCoeffs[:, :] = coeffs[span1-deg1:span1+1, span2-deg2:span2+1] - z = 0.0 + z = 0.0*coeffs[0,0] for i in range(deg1+1): theCoeffs[i, 0] = theCoeffs[i, 0]*basis2[0] for j in range(1, deg2+1): @@ -249,7 +250,7 @@ def nu_eval_spline_2d_scalar(x: 'float', y: 'float', kts1: 'Final[float[:]]', de @pure @stack_array('basis1', 'basis2', 'theCoeffs') def nu_eval_spline_2d_cross(X: 'Final[float[:]]', Y: 'Final[float[:]]', kts1: 'Final[float[:]]', deg1: 'int', kts2: 'Final[float[:]]', deg2: 'int', - coeffs: 'Final[float[:,:]]', z: 'float[:,:]', der1: 'int' = 0, der2: 'int' = 0): + coeffs: 'Final[CoeffType[:,:]]', z: 'CoeffType[:,:]', der1: 'int' = 0, der2: 'int' = 0): """ TODO """ @@ -337,7 +338,7 @@ def nu_eval_spline_2d_cross(X: 'Final[float[:]]', Y: 'Final[float[:]]', kts1: 'F @pure @stack_array('basis1', 'basis2', 'theCoeffs') def nu_eval_spline_2d_vector(x: 'float[:]', y: 'float[:]', kts1: 'float[:]', deg1: 'int', kts2: 'float[:]', deg2: 'int', - coeffs: 'float[:,:]', z: 'float[:]', der1: 'int' = 0, der2: 'int' = 0): + coeffs: 'CoeffType[:,:]', z: 'CoeffType[:]', der1: 'int' = 0, der2: 'int' = 0): """ TODO """ From d1da492c0657073c3933be76d9125be063fa5793 Mon Sep 17 00:00:00 2001 From: Emily Bourne Date: Wed, 20 Aug 2025 22:10:59 +0200 Subject: [PATCH 14/30] Can't check basis --- pygyro/splines/spline_interpolators.py | 9 +++++---- 1 file changed, 5 insertions(+), 4 deletions(-) diff --git a/pygyro/splines/spline_interpolators.py b/pygyro/splines/spline_interpolators.py index f67233fb..effe1137 100644 --- a/pygyro/splines/spline_interpolators.py +++ b/pygyro/splines/spline_interpolators.py @@ -65,7 +65,7 @@ def compute_interpolant(self, ug, spl): """ assert isinstance(spl, (Spline1D, Spline1DComplex)) - assert spl.basis is self._basis + #assert spl.basis is self._basis assert len(ug) == self._basis.nbasis if self._basis.periodic: @@ -221,9 +221,10 @@ def compute_interpolant(self, ug, spl): """ assert isinstance(spl, Spline2D) - basis1, basis2 = spl.basis - assert basis1 is self._basis1 - assert basis2 is self._basis2 + basis1 = spl.basis1 + basis2 = spl.basis2 + #assert basis1 is self._basis1 + #assert basis2 is self._basis2 n1, n2 = basis1.nbasis, basis2.nbasis p1, p2 = basis1.degree, basis2.degree From 28244bd2394844f7821ee44e54943b0e7a31551d Mon Sep 17 00:00:00 2001 From: Emily Bourne Date: Wed, 20 Aug 2025 22:12:19 +0200 Subject: [PATCH 15/30] Annotate methods. Reorder isinstance check to be Pyccel compatible. Avoid cspan limits. Specify allow_negative_indices. Unpack explicitly --- pygyro/splines/splines.py | 164 +++++++++++++++++++++++--------------- 1 file changed, 99 insertions(+), 65 deletions(-) diff --git a/pygyro/splines/splines.py b/pygyro/splines/splines.py index 75470ace..f78cad0b 100644 --- a/pygyro/splines/splines.py +++ b/pygyro/splines/splines.py @@ -1,5 +1,7 @@ # coding: utf-8 # Copyright 2018 Yaman Güçlü +from pyccel.decorators import allow_negative_index +from typing import TypeVar, Final import numpy as np # from scipy.interpolate import splev, bisplev @@ -9,12 +11,14 @@ from .cubic_uniform_spline_eval_funcs import cu_eval_spline_1d_scalar, cu_eval_spline_1d_vector from .cubic_uniform_spline_eval_funcs import cu_eval_spline_2d_cross, cu_eval_spline_2d_scalar +ScalOrArr = TypeVar('ScalOrArr', float, 'float[:]') + __all__ = ['make_knots', 'BSplines', 'Spline1D', 'Spline1DComplex', 'Spline2D'] # =============================================================================== - -def make_knots(breaks, degree, periodic): +@allow_negative_index('breaks', 'T') +def make_knots(breaks : 'Final[float[:]]', degree : Final[int], periodic : Final[bool]): """ Create spline knots from breakpoints, with appropriate boundary conditions. Let p be spline degree. If domain is periodic, knot sequence is extended @@ -45,19 +49,20 @@ def make_knots(breaks, degree, periodic): # Consistency checks assert len(breaks) > 1 - assert all(np.diff(breaks) > 0) + #assert all(np.diff(breaks) > 0) assert degree > 0 if periodic: assert len(breaks) > degree p = degree - T = np.zeros(len(breaks)+2*p) + n = len(breaks) + T = np.zeros(n+2*p) T[p:-p] = breaks if periodic: period = breaks[-1]-breaks[0] - T[0:p] = [xi-period for xi in breaks[-p-1:-1]] - T[-p:] = [xi+period for xi in breaks[1:p+1]] + T[0:p] = [breaks[i-p-1]-period for i in range(p)] + T[-p:] = [breaks[i]+period for i in range(1, p+1)] else: T[0:p] = breaks[0] T[-p:] = breaks[-1] @@ -92,7 +97,8 @@ class BSplines(): """ - def __init__(self, knots, degree, periodic, uniform): + @allow_negative_index('knots') + def __init__(self, knots : 'float[:]', degree : int, periodic : bool, uniform : bool): xmin = knots[degree] xmax = knots[-degree-1] dx = knots[degree+1]-knots[degree] @@ -104,13 +110,13 @@ def __init__(self, knots, degree, periodic, uniform): self._ncells = len(knots)-2*degree-1 self._nbasis = self._ncells if periodic else self._ncells+degree self._offset = degree//2 if periodic else 0 - self._integrals = None + self._integrals = np.empty(0) if self._cubic_uniform_splines: - self._knots = np.array([xmin, xmax, dx, self._ncells]) + self._knots = np.array([xmin, xmax, dx, float(self._ncells)]) assert (int(self._knots[3]) == self._ncells) else: - self._knots = knots + self._knots = np.array(knots) self._build_integrals() @@ -119,11 +125,19 @@ def __init__(self, knots, degree, periodic, uniform): self._interp_pts = np.linspace( xmin, xmax, self._ncells, endpoint=False) else: - self._interp_pts = np.array([xmin, - xmin+dx/3, - *np.linspace(xmin+dx, xmax-dx, self._nbasis-4), - xmax-dx/3, - xmax]) + self._interp_pts = np.empty(self._nbasis) + self._interp_pts[0] = xmin + self._interp_pts[1] = xmin+dx/3 + self._interp_pts[2:-2] = np.linspace(xmin+dx, xmax-dx, self._nbasis-4) + self._interp_pts[-2] = xmax-dx/3 + self._interp_pts[-1] = xmax + #self._interp_pts = np.array([xmin, + # xmin+dx/3, + # *np.linspace(xmin+dx, xmax-dx, self._nbasis-4), + # xmax-dx/3, + # xmax]) + else: + self._interp_pts = np.empty(0) @property def degree(self): @@ -163,8 +177,9 @@ def breaks(self): xmin, xmax, _, _ = self._knots return np.linspace(xmin, xmax, self._ncells+1) else: + n = len(self._knots) p = self._degree - return self._knots[p:-p] + return np.array(self._knots[p:n-p]) @property def domain(self): @@ -182,7 +197,7 @@ def greville(self): """ Coordinates of all Greville points. """ if self._cubic_uniform_splines: - return self._interp_pts + return np.array(self._interp_pts) else: p = self._degree n = self._nbasis @@ -192,17 +207,18 @@ def greville(self): if self._periodic: a, b = self.domain - x = np.around(x, decimals=15) - x = (x-a) % (b-a) + a + #x = np.around(x, decimals=15) + x[:] = (x-a) % (b-a) + a - return np.around(x, decimals=15) + #return np.around(x, decimals=15) + return x @property def integrals(self): return self._integrals # ... - def __getitem__(self, i): + def __getitem__(self, i : int): """ Get the i-th basis function as a 1D spline. @@ -226,6 +242,7 @@ def __getitem__(self, i): spl.coeffs[n:n+p] = spl.coeffs[0:p] return spl + @allow_negative_index('values') def _build_integrals(self): n = self.nbasis d = self.degree @@ -239,7 +256,7 @@ def _build_integrals(self): self._integrals[:] = dx self._integrals[n:] = 0 else: - self._integrals[d:-d] = dx + self._integrals[d:self.ncells] = dx values = np.empty(d+2) knots = np.linspace(xmin, xmin+dx*11, 12) test_pt = xmin + 4*dx @@ -247,11 +264,15 @@ def _build_integrals(self): nu_basis_funs(knots, 4, test_pt, span, values) for i in range(3): - step = dx*(1 - sum(values[:3-i])) + step = dx*(1 - np.sum(values[:3-i])) self._integrals[i] = step - self._integrals[-i-1] = step + self._integrals[self.ncells+d-i-1] = step else: - knots = np.array([self.knots[0], *self.knots, self.knots[-1]]) + knots = np.empty(len(self.knots)+2) + knots[0] = self.knots[0] + knots[1:-1] = self.knots + knots[-1] = self.knots[-1] + #knots = np.array([self.knots[0], *self.knots, self.knots[-1]]) values = np.empty(d+2) for i in range(n): @@ -265,13 +286,19 @@ def _build_integrals(self): first_available = span_l - integ_deg first_wanted = i+1 min_idx = first_wanted-first_available - l = np.sum(values[min_idx:]) + if min_idx >= len(values): + l = 0.0 + else: + l = np.sum(values[min_idx:]) nu_basis_funs(knots, integ_deg, ubound, span_u, values) first_available = span_u - integ_deg first_wanted = i+1 min_idx = first_wanted-first_available - u = np.sum(values[min_idx:]) + if min_idx >= len(values): + u = 0.0 + else: + u = np.sum(values[min_idx:]) self._integrals[i] = ( knots[d+2+i] - knots[i+1])*inv_deg*(u - l) @@ -288,7 +315,7 @@ class Spline1D(): TODO """ - def __init__(self, basis): + def __init__(self, basis : BSplines): assert isinstance(basis, BSplines) self._basis = basis self._coeffs = np.zeros(basis.ncells + basis.degree, dtype=float) @@ -307,11 +334,18 @@ def coeffs(self): """ return self._coeffs - def eval(self, x, der=0): + def eval(self, x : 'float|float[:]', der : int = 0): """ TODO """ - if (hasattr(x, '__len__')): + if isinstance(x, float): + if self._basis.cubic_uniform: + result = cu_eval_spline_1d_scalar( + x, self._basis.knots, self._basis.degree, self._coeffs, der) + else: + result = nu_eval_spline_1d_scalar( + x, self._basis.knots, self._basis.degree, self._coeffs, der) + else: result = np.empty_like(x) if self._basis.cubic_uniform: cu_eval_spline_1d_vector(x, self._basis.knots, @@ -319,13 +353,6 @@ def eval(self, x, der=0): else: nu_eval_spline_1d_vector(x, self._basis.knots, self._basis.degree, self._coeffs, result, der) - else: - if self._basis.cubic_uniform: - result = cu_eval_spline_1d_scalar( - x, self._basis.knots, self._basis.degree, self._coeffs, der) - else: - result = nu_eval_spline_1d_scalar( - x, self._basis.knots, self._basis.degree, self._coeffs, der) return result """ @@ -333,7 +360,7 @@ def eval(self, x, der=0): return splev( x, tck, der ) """ - def eval_vector(self, x, y, der=0): + def eval_vector(self, x : 'float[:]', y : 'float[:]', der : int=0): """ TODO """ @@ -349,7 +376,7 @@ class Spline1DComplex(): TODO """ - def __init__(self, basis): + def __init__(self, basis : BSplines): assert isinstance(basis, BSplines) self._basis = basis self._coeffs = np.zeros(basis.ncells + basis.degree, dtype=np.complex128) @@ -368,25 +395,25 @@ def coeffs(self): """ return self._coeffs - def eval(self, x, der=0): + def eval(self, x : 'float|float[:]', der : int = 0): """ TODO """ - if (hasattr(x, '__len__')): - result = np.empty_like(x) - if self._basis.cubic_uniform: - cu_eval_spline_1d_vector(x, self._basis.knots, - self._basis.degree, self._coeffs, result, der) - else: - nu_eval_spline_1d_vector(x, self._basis.knots, - self._basis.degree, self._coeffs, result, der) - else: + if isinstance(x, float): if self._basis.cubic_uniform: result = cu_eval_spline_1d_scalar( x, self._basis.knots, self._basis.degree, self._coeffs, der) else: result = nu_eval_spline_1d_scalar( x, self._basis.knots, self._basis.degree, self._coeffs, der) + else: + result = np.empty_like(x, dtype=np.complex128) + if self._basis.cubic_uniform: + cu_eval_spline_1d_vector(x, self._basis.knots, + self._basis.degree, self._coeffs, result, der) + else: + nu_eval_spline_1d_vector(x, self._basis.knots, + self._basis.degree, self._coeffs, result, der) return result """ @@ -394,7 +421,7 @@ def eval(self, x, der=0): return splev( x, tck, der ) """ - def eval_vector(self, x, y, der=0): + def eval_vector(self, x : 'float[:]', y : 'complex[:]', der : int=0): """ TODO """ @@ -413,7 +440,7 @@ class Spline2D(): TODO """ - def __init__(self, basis1, basis2): + def __init__(self, basis1 : BSplines, basis2 : BSplines): assert isinstance(basis1, BSplines) assert isinstance(basis2, BSplines) shape = (basis1.ncells + basis1.degree, basis2.ncells + basis2.degree) @@ -426,11 +453,18 @@ def __init__(self, basis1, basis2): assert basis1.cubic_uniform == basis2.cubic_uniform @property - def basis(self): + def basis1(self): + """ + TODO + """ + return self._basis1 + + @property + def basis2(self): """ TODO """ - return self._basis1, self._basis2 + return self._basis2 @property def coeffs(self): @@ -439,11 +473,20 @@ def coeffs(self): """ return self._coeffs - def eval(self, x1, x2, der1=0, der2=0): + def eval(self, x1 : ScalOrArr, x2 : ScalOrArr, der1 : int=0, der2 : int=0): """ TODO """ - if (hasattr(x1, '__len__')): + if isinstance(x1, float): + if self._basis1.cubic_uniform: + result = cu_eval_spline_2d_scalar(x1, x2, self._basis1.knots, self._basis1.degree, + self._basis2.knots, self._basis2.degree, + self._coeffs, der1, der2) + else: + result = nu_eval_spline_2d_scalar(x1, x2, self._basis1.knots, self._basis1.degree, + self._basis2.knots, self._basis2.degree, + self._coeffs, der1, der2) + else: result = np.empty((len(x1), len(x2))) if self._basis1.cubic_uniform: cu_eval_spline_2d_cross(x1, x2, self._basis1.knots, self._basis1.degree, @@ -453,15 +496,6 @@ def eval(self, x1, x2, der1=0, der2=0): nu_eval_spline_2d_cross(x1, x2, self._basis1.knots, self._basis1.degree, self._basis2.knots, self._basis2.degree, self._coeffs, result, der1, der2) - else: - if self._basis1.cubic_uniform: - result = cu_eval_spline_2d_scalar(x1, x2, self._basis1.knots, self._basis1.degree, - self._basis2.knots, self._basis2.degree, - self._coeffs, der1, der2) - else: - result = nu_eval_spline_2d_scalar(x1, x2, self._basis1.knots, self._basis1.degree, - self._basis2.knots, self._basis2.degree, - self._coeffs, der1, der2) return result """ @@ -475,7 +509,7 @@ def eval(self, x1, x2, der1=0, der2=0): return bisplev( x1, x2, tck, der1, der2 ) """ - def eval_vector(self, x1, x2, y, der1=0, der2=0): + def eval_vector(self, x1 : 'float[:]', x2 : 'float[:]', y : 'float[:,:]', der1 : int=0, der2 : int=0): """ TODO """ From 63cb6933052c98d6d17d4604ae5522c5f6b85c4f Mon Sep 17 00:00:00 2001 From: Emily Bourne Date: Wed, 20 Aug 2025 22:44:01 +0200 Subject: [PATCH 16/30] Update Makefile. C is needed due to #1339 --- Makefile | 2 +- pygyro/advection/Makefile | 2 +- pygyro/splines/Makefile | 9 ++++++++- 3 files changed, 10 insertions(+), 3 deletions(-) diff --git a/Makefile b/Makefile index 50f7384c..24d1e1fa 100644 --- a/Makefile +++ b/Makefile @@ -9,7 +9,7 @@ ACC := pycc COMP := GNU # Target language -LANGUAGE := fortran +LANGUAGE := c #PYTHRAN_FLAGS := -DUSE_XSIMD -fopenmp -march=native PYTHRAN_FLAGS := diff --git a/pygyro/advection/Makefile b/pygyro/advection/Makefile index 901db019..a1685e73 100644 --- a/pygyro/advection/Makefile +++ b/pygyro/advection/Makefile @@ -18,7 +18,7 @@ all: accelerated_advection_steps$(SO_EXT) @rm -f .ACC.* @touch $@ -DEPS := ../initialisation/$(NAME_PREFIX)initialiser_funcs.py ../splines/$(NAME_PREFIX)spline_eval_funcs.py ../splines/$(NAME_PREFIX)cubic_uniform_spline_eval_funcs.py +DEPS := ../initialisation/$(NAME_PREFIX)initialiser_funcs.py ../splines/$(NAME_PREFIX)splines.py ifneq ($(ACC), numba) ifneq ($(ACC), pythran) diff --git a/pygyro/splines/Makefile b/pygyro/splines/Makefile index 38f9d7dc..c2e7a40f 100644 --- a/pygyro/splines/Makefile +++ b/pygyro/splines/Makefile @@ -8,7 +8,7 @@ # Main targets #---------------------------------------------------------- -all: spline_eval_funcs$(SO_EXT) cubic_uniform_spline_eval_funcs$(SO_EXT) +all: spline_eval_funcs$(SO_EXT) cubic_uniform_spline_eval_funcs$(SO_EXT) splines$(SO_EXT) .ACC.$(ACC): @rm -f .ACC.* @@ -28,6 +28,13 @@ else $(TOOL) $< $(TOOL_FLAGS) -o $@ endif +splines$(SO_EXT): $(NAME_PREFIX)splines.py .ACC.$(ACC) spline_eval_funcs$(SO_EXT) cubic_uniform_spline_eval_funcs$(SO_EXT) +ifneq ($(ACC), pythran) + $(TOOL) $< $(TOOL_FLAGS) +else + $(TOOL) $< $(TOOL_FLAGS) -o $@ +endif + clean: rm -f *.o *.so *.mod .ACC.pycc .lock_acquisition.lock From 4737da442384c3dc6484b57e3ce55e9b83f875b3 Mon Sep 17 00:00:00 2001 From: Emily Bourne Date: Wed, 20 Aug 2025 22:46:42 +0200 Subject: [PATCH 17/30] Update imports --- pygyro/poisson/poisson_solver.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/pygyro/poisson/poisson_solver.py b/pygyro/poisson/poisson_solver.py index 396d6b96..5b8114d1 100644 --- a/pygyro/poisson/poisson_solver.py +++ b/pygyro/poisson/poisson_solver.py @@ -5,7 +5,7 @@ from numpy.polynomial.legendre import leggauss from ..model.grid import Grid -from ..splines.splines import BSplines, Spline1D, make_knots +from ..splines.splines import BSplines, Spline1D, Spline1DComplex, make_knots from ..splines.spline_interpolators import SplineInterpolator1D from ..initialisation import initialiser_funcs as init from .poisson_tools import get_perturbed_rho, get_rho From cd6257825247f3842b015dca5d7dc81595a4e673 Mon Sep 17 00:00:00 2001 From: Emily Bourne Date: Wed, 20 Aug 2025 22:46:54 +0200 Subject: [PATCH 18/30] Use splines directly --- .../advection/accelerated_advection_steps.py | 159 ++++-------------- pygyro/advection/advection.py | 33 +--- 2 files changed, 39 insertions(+), 153 deletions(-) diff --git a/pygyro/advection/accelerated_advection_steps.py b/pygyro/advection/accelerated_advection_steps.py index e1939581..5cd5dd3f 100644 --- a/pygyro/advection/accelerated_advection_steps.py +++ b/pygyro/advection/accelerated_advection_steps.py @@ -1,24 +1,18 @@ from typing import Final from pyccel.decorators import pure -from ..splines.spline_eval_funcs import nu_eval_spline_1d_scalar, nu_eval_spline_1d_vector -from ..splines.spline_eval_funcs import nu_eval_spline_2d_cross, nu_eval_spline_2d_scalar -from ..splines.cubic_uniform_spline_eval_funcs import cu_eval_spline_1d_scalar, cu_eval_spline_1d_vector -from ..splines.cubic_uniform_spline_eval_funcs import cu_eval_spline_2d_cross, cu_eval_spline_2d_scalar +from ..splines.splines import Spline1D, Spline2D from ..initialisation.initialiser_funcs import f_eq -def general_poloidal_advection_step_expl(f: 'float[:,:]', +def poloidal_advection_step_expl(f: 'float[:,:]', dt: 'float', v: 'float', rPts: 'Final[float[:]]', qPts: 'Final[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[:,:]', - kts1Phi: 'float[:]', kts2Phi: 'float[:]', coeffsPhi: 'float[:,:]', deg1Phi: 'int', deg2Phi: 'int', - kts1Pol: 'float[:]', kts2Pol: 'float[:]', coeffsPol: 'float[:,:]', deg1Pol: 'int', deg2Pol: 'int', + phi_spline : Spline2D, pol_spline : Spline2D, CN0: 'float', kN0: 'float', deltaRN0: 'float', rp: 'float', CTi: 'float', - kTi: 'float', deltaRTi: 'float', B0: 'float', nulBound: 'bool', - eval_spline_2d_cross: '()(Final[float[:]], Final[float[:]], Final[float[:]], int, Final[float[:]], int, Final[float[:,:]], float[:,:], int, int)', - eval_spline_2d_scalar: '(float)(float, float, Final[float[:]], int, Final[float[:]], int, Final[float[:,:]], int, int)'): + kTi: 'float', deltaRTi: 'float', B0: 'float', nulBound: 'bool'): """ Carry out an advection step for the poloidal advection @@ -43,10 +37,8 @@ def general_poloidal_advection_step_expl(f: 'float[:,:]', multFactor = dt / B0 multFactor_half = 0.5 * multFactor - eval_spline_2d_cross(qPts, rPts, kts1Phi, deg1Phi, - kts2Phi, deg2Phi, coeffsPhi, drPhi_0, 0, 1) - eval_spline_2d_cross(qPts, rPts, kts1Phi, deg1Phi, - kts2Phi, deg2Phi, coeffsPhi, dthetaPhi_0, 1, 0) + phi_spline.eval_vector(qPts, rPts, drPhi_0, 0, 1) + phi_spline.eval_vector(qPts, rPts, dthetaPhi_0, 1, 0) nPts_r = rPts.shape[0] nPts_q = qPts.shape[0] @@ -71,14 +63,12 @@ def general_poloidal_advection_step_expl(f: 'float[:,:]', # Add the new value of phi to the derivatives # x^{n+1} = x^n + 0.5( f(x^n) + f(x^n + f(x^n)) ) # ^^^^^^^^^^^^^^^ - drPhi_k[i, j] = eval_spline_2d_scalar(endPts_k1_q[i, j], endPts_k1_r[i, j], - kts1Phi, deg1Phi, kts2Phi, deg2Phi, - coeffsPhi, 0, 1) + drPhi_k[i, j] = phi_spline.eval(endPts_k1_q[i, j], endPts_k1_r[i, j], + 0, 1) drPhi_k[i, j] /= endPts_k1_r[i, j] - dthetaPhi_k[i, j] = eval_spline_2d_scalar(endPts_k1_q[i, j], endPts_k1_r[i, j], - kts1Phi, deg1Phi, kts2Phi, deg2Phi, - coeffsPhi, 1, 0) + dthetaPhi_k[i, j] = phi_spline.eval(endPts_k1_q[i, j], endPts_k1_r[i, j], + 1, 0) dthetaPhi_k[i, j] /= endPts_k1_r[i, j] else: drPhi_k[i, j] = 0.0 @@ -102,9 +92,7 @@ def general_poloidal_advection_step_expl(f: 'float[:,:]', f[i, j] = 0.0 else: endPts_k2_q[i, j] = endPts_k2_q[i, j] % (2*pi) - f[i, j] = eval_spline_2d_scalar(endPts_k2_q[i, j], endPts_k2_r[i, j], - kts1Pol, deg1Pol, kts2Pol, deg2Pol, - coeffsPol, 0, 0) + f[i, j] = pol_spline.eval(endPts_k2_q[i, j], endPts_k2_r[i, j]) else: for i in range(nPts_q): # theta for j in range(nPts_r): # r @@ -116,39 +104,14 @@ def general_poloidal_advection_step_expl(f: 'float[:,:]', deltaRN0, rp, CTi, kTi, deltaRTi) else: endPts_k2_q[i, j] = endPts_k2_q[i, j] % (2*pi) - f[i, j] = eval_spline_2d_scalar(endPts_k2_q[i, j], endPts_k2_r[i, j], - kts1Pol, deg1Pol, kts2Pol, deg2Pol, coeffsPol, 0, 0) - - -def poloidal_advection_step_expl(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[:,:]', - kts1Phi: 'float[:]', kts2Phi: 'float[:]', coeffsPhi: 'float[:,:]', deg1Phi: 'int', deg2Phi: 'int', - kts1Pol: 'float[:]', kts2Pol: 'float[:]', coeffsPol: 'float[:,:]', deg1Pol: 'int', deg2Pol: 'int', - CN0: 'float', kN0: 'float', deltaRN0: 'float', rp: 'float', CTi: 'float', - kTi: 'float', deltaRTi: 'float', B0: 'float', cubic_uniform_splines: 'bool', nulBound: 'bool' = False): - if cubic_uniform_splines: - general_poloidal_advection_step_expl(f, dt, v, rPts, qPts, drPhi_0, dthetaPhi_0, drPhi_k, dthetaPhi_k, endPts_k1_q, - endPts_k1_r, endPts_k2_q, endPts_k2_r, kts1Phi, kts2Phi, coeffsPhi, - deg1Phi, deg2Phi, kts1Pol, kts2Pol, coeffsPol, deg1Pol, deg2Pol, CN0, kN0, deltaRN0, - rp, CTi, kTi, deltaRTi, B0, nulBound, cu_eval_spline_2d_cross, cu_eval_spline_2d_scalar) - else: - general_poloidal_advection_step_expl(f, dt, v, rPts, qPts, drPhi_0, dthetaPhi_0, drPhi_k, dthetaPhi_k, endPts_k1_q, - endPts_k1_r, endPts_k2_q, endPts_k2_r, kts1Phi, kts2Phi, coeffsPhi, - deg1Phi, deg2Phi, kts1Pol, kts2Pol, coeffsPol, deg1Pol, deg2Pol, CN0, kN0, deltaRN0, - rp, CTi, kTi, deltaRTi, B0, nulBound, nu_eval_spline_2d_cross, nu_eval_spline_2d_scalar) + f[i, j] = pol_spline.eval(endPts_k2_q[i, j], endPts_k2_r[i, j]) -def general_v_parallel_advection_eval_step(f: 'float[:]', vPts: 'float[:]', +def v_parallel_advection_eval_step(f: 'float[:]', vPts: 'float[:]', rPos: 'float', vMin: 'float', vMax: 'float', - kts: 'float[:]', deg: 'int', - coeffs: 'float[:]', + spl : Spline1D, CN0: 'float', kN0: 'float', deltaRN0: 'float', rp: 'float', - CTi: 'float', kTi: 'float', deltaRTi: 'float', bound: 'int', - eval_spline_1d_scalar: '(float)(float, Final[float[:]], int, Final[float[:]], int)'): + CTi: 'float', kTi: 'float', deltaRTi: 'float', bound: 'int'): """ TODO """ @@ -159,13 +122,13 @@ def general_v_parallel_advection_eval_step(f: 'float[:]', vPts: 'float[:]', f[i] = f_eq(rPos, v, CN0, kN0, deltaRN0, rp, CTi, kTi, deltaRTi) else: - f[i] = eval_spline_1d_scalar(v, kts, deg, coeffs, 0) + f[i] = spl.eval(v) elif (bound == 1): for i, v in enumerate(vPts): if (v < vMin or v > vMax): f[i] = 0.0 else: - f[i] = eval_spline_1d_scalar(v, kts, deg, coeffs, 0) + f[i] = spl.eval(v) elif (bound == 2): vDiff = vMax - vMin for i, v in enumerate(vPts): @@ -173,29 +136,12 @@ def general_v_parallel_advection_eval_step(f: 'float[:]', vPts: 'float[:]', v += vDiff while (v > vMax): v -= vDiff - f[i] = eval_spline_1d_scalar(v, kts, deg, coeffs, 0) - - -def v_parallel_advection_eval_step(f: 'float[:]', vPts: 'float[:]', - rPos: 'float', vMin: 'float', vMax: 'float', - kts: 'float[:]', deg: 'int', - coeffs: 'float[:]', - CN0: 'float', kN0: 'float', deltaRN0: 'float', rp: 'float', - CTi: 'float', kTi: 'float', deltaRTi: 'float', bound: 'int', cubic_uniform_splines: 'bool'): - if cubic_uniform_splines: - general_v_parallel_advection_eval_step(f, vPts, rPos, vMin, vMax, kts, deg, coeffs, - CN0, kN0, deltaRN0, rp, CTi, kTi, deltaRTi, bound, cu_eval_spline_1d_scalar) - else: - general_v_parallel_advection_eval_step(f, vPts, rPos, vMin, vMax, kts, deg, coeffs, - CN0, kN0, deltaRN0, rp, CTi, kTi, deltaRTi, bound, nu_eval_spline_1d_scalar) + f[i] = spl.eval(v) -def general_get_lagrange_vals(i: 'int', shifts: 'int[:]', +def get_lagrange_vals(i: 'int', shifts: 'int[:]', vals: 'float[:,:,:]', qVals: 'float[:]', - thetaShifts: 'float[:]', kts: 'float[:]', - deg: 'int', coeffs: 'float[:]', - eval_spline_1d_vector: '()(Final[float[:]], Final[float[:]], int, Final[float[:]], float[:], int)', - eval_spline_1d_scalar: '(float)(float, Final[float[:]], int, Final[float[:]], int)'): + thetaShifts: 'float[:]', spl : Spline1D): """ TODO """ @@ -206,21 +152,7 @@ def general_get_lagrange_vals(i: 'int', shifts: 'int[:]', for j, s in enumerate(shifts): idx = (i - s) % nz new_q[:] = (qVals + thetaShifts[j]) % (2*pi) - # eval_spline_1d_vector(new_q, kts, deg, coeffs, vals[idx, :, j], 0) - for k, q in enumerate(new_q): - vals[idx, k, j] = eval_spline_1d_scalar(q, kts, deg, coeffs, 0) - - -def get_lagrange_vals(i: 'int', shifts: 'int[:]', - vals: 'float[:,:,:]', qVals: 'float[:]', - thetaShifts: 'float[:]', kts: 'float[:]', - deg: 'int', coeffs: 'float[:]', cubic_uniform_splines: 'bool'): - if cubic_uniform_splines: - general_get_lagrange_vals( - i, shifts, vals, qVals, thetaShifts, kts, deg, coeffs, cu_eval_spline_1d_vector, cu_eval_spline_1d_scalar) - else: - general_get_lagrange_vals( - i, shifts, vals, qVals, thetaShifts, kts, deg, coeffs, nu_eval_spline_1d_vector, nu_eval_spline_1d_scalar) + spl.eval_vector(new_q, vals[idx, :, j]) @pure @@ -236,15 +168,12 @@ def flux_advection(nq: 'int', nr: 'int', f[j, i] += coeffs[k]*vals[i, j, k] -def general_poloidal_advection_step_impl(f: 'float[:,:]', dt: 'float', v: 'float', rPts: 'float[:]', qPts: 'float[:]', +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[:,:]', - kts1Phi: 'float[:]', kts2Phi: 'float[:]', coeffsPhi: 'float[:,:]', deg1Phi: 'int', deg2Phi: 'int', - kts1Pol: 'float[:]', kts2Pol: 'float[:]', coeffsPol: 'float[:,:]', deg1Pol: 'int', deg2Pol: 'int', + 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', - eval_spline_2d_cross: '()(Final[float[:]], Final[float[:]], Final[float[:]], int, Final[float[:]], int, Final[float[:,:]], float[:,:], int, int)', - eval_spline_2d_scalar: '(float)(float, float, Final[float[:]], int, Final[float[:]], int, Final[float[:,:]], int, int)'): + B0: 'float', tol: 'float', nulBound: 'bool'): """ Carry out an advection step for the poloidal advection @@ -268,10 +197,8 @@ def general_poloidal_advection_step_impl(f: 'float[:,:]', dt: 'float', v: 'float multFactor = dt/B0 - eval_spline_2d_cross(qPts, rPts, kts1Phi, deg1Phi, - kts2Phi, deg2Phi, coeffsPhi, drPhi_0, 0, 1) - eval_spline_2d_cross(qPts, rPts, kts1Phi, deg1Phi, - kts2Phi, deg2Phi, coeffsPhi, dthetaPhi_0, 1, 0) + phi_spline.eval_vector(qPts, rPts, drPhi_0, 0, 1) + phi_spline.eval_vector(qPts, rPts, dthetaPhi_0, 1, 0) nPts_r = rPts.shape[0] nPts_q = qPts.shape[0] @@ -303,13 +230,9 @@ def general_poloidal_advection_step_impl(f: 'float[:,:]', dt: 'float', v: 'float # Add the new value of phi to the derivatives # x^{n+1} = x^n + 0.5( f(x^n) + f(x^n + f(x^n)) ) # ^^^^^^^^^^^^^^^ - drPhi_k[i, j] = eval_spline_2d_scalar(endPts_k1_q[i, j], endPts_k1_r[i, j], - kts1Phi, deg1Phi, kts2Phi, deg2Phi, - coeffsPhi, 0, 1) + drPhi_k[i, j] = phi_spline.eval(endPts_k1_q[i, j], endPts_k1_r[i, j], 0, 1) drPhi_k[i, j] /= endPts_k1_r[i, j] - dthetaPhi_k[i, j] = eval_spline_2d_scalar(endPts_k1_q[i, j], endPts_k1_r[i, j], - kts1Phi, deg1Phi, kts2Phi, deg2Phi, - coeffsPhi, 1, 0) + dthetaPhi_k[i, j] = phi_spline.eval(endPts_k1_q[i, j], endPts_k1_r[i, j], 1, 0) dthetaPhi_k[i, j] /= endPts_k1_r[i, j] else: drPhi_k[i, j] = 0.0 @@ -352,9 +275,7 @@ def general_poloidal_advection_step_impl(f: 'float[:,:]', dt: 'float', v: 'float f[i, j] = 0.0 else: endPts_k2_q[i, j] = endPts_k2_q[i, j] % (2*pi) - f[i, j] = eval_spline_2d_scalar(endPts_k2_q[i, j], endPts_k2_r[i, j], - kts1Pol, deg1Pol, kts2Pol, deg2Pol, - coeffsPol, 0, 0) + f[i, j] = pol_spline.eval(endPts_k2_q[i, j], endPts_k2_r[i, j]) else: for i in range(nPts_q): for j in range(nPts_r): @@ -366,24 +287,4 @@ def general_poloidal_advection_step_impl(f: 'float[:,:]', dt: 'float', v: 'float deltaRN0, rp, CTi, kTi, deltaRTi) else: endPts_k2_q[i, j] = endPts_k2_q[i, j] % (2*pi) - f[i, j] = eval_spline_2d_scalar(endPts_k2_q[i, j], endPts_k2_r[i, j], - kts1Pol, deg1Pol, kts2Pol, deg2Pol, coeffsPol, 0, 0) - - -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[:,:]', - kts1Phi: 'float[:]', kts2Phi: 'float[:]', coeffsPhi: 'float[:,:]', deg1Phi: 'int', deg2Phi: 'int', - kts1Pol: 'float[:]', kts2Pol: 'float[:]', coeffsPol: 'float[:,:]', deg1Pol: 'int', deg2Pol: 'int', - CN0: 'float', kN0: 'float', deltaRN0: 'float', rp: 'float', CTi: 'float', kTi: 'float', deltaRTi: 'float', - B0: 'float', tol: 'float', cubic_uniform_splines: 'bool', nulBound: 'bool' = False): - if cubic_uniform_splines: - general_poloidal_advection_step_impl(f, dt, v, rPts, qPts, drPhi_0, dthetaPhi_0, drPhi_k, dthetaPhi_k, endPts_k1_q, - endPts_k1_r, endPts_k2_q, endPts_k2_r, kts1Phi, kts2Phi, coeffsPhi, - deg1Phi, deg2Phi, kts1Pol, kts2Pol, coeffsPol, deg1Pol, deg2Pol, CN0, kN0, deltaRN0, - rp, CTi, kTi, deltaRTi, B0, tol, nulBound, cu_eval_spline_2d_cross, cu_eval_spline_2d_scalar) - else: - general_poloidal_advection_step_impl(f, dt, v, rPts, qPts, drPhi_0, dthetaPhi_0, drPhi_k, dthetaPhi_k, endPts_k1_q, - endPts_k1_r, endPts_k2_q, endPts_k2_r, kts1Phi, kts2Phi, coeffsPhi, - deg1Phi, deg2Phi, kts1Pol, kts2Pol, coeffsPol, deg1Pol, deg2Pol, CN0, kN0, deltaRN0, - rp, CTi, kTi, deltaRTi, B0, tol, nulBound, nu_eval_spline_2d_cross, nu_eval_spline_2d_scalar) + f[i, j] = pol_spline.eval(endPts_k2_q[i, j], endPts_k2_r[i, j]) diff --git a/pygyro/advection/advection.py b/pygyro/advection/advection.py index 6bfb5045..79a0bef6 100644 --- a/pygyro/advection/advection.py +++ b/pygyro/advection/advection.py @@ -283,10 +283,7 @@ def step(self, f: np.ndarray, cIdx: int, rIdx: int = 0): get_lagrange_vals(i, self._shifts[rIdx, cIdx], self._LagrangeVals, self._points[0], self._thetaShifts[rIdx, cIdx], - self._thetaSpline.basis.knots, - self._thetaSpline.basis.degree, - self._thetaSpline.coeffs, - self._thetaSpline.basis.cubic_uniform) + self._thetaSpline) flux_advection(*self._nPoints, f, self._lagrangeCoeffs[rIdx, cIdx], @@ -366,12 +363,11 @@ def step(self, f: np.ndarray, dt: float, c: float, r: float): self._interpolator.compute_interpolant(f, self._spline) v_parallel_advection_eval_step(f, self._points-c*dt, r, self._points[0], - self._points[-1], self._spline.basis.knots, - self._spline.basis.degree, self._spline.coeffs, + self._points[-1], self._spline, self._constants.CN0, self._constants.kN0, self._constants.deltaRN0, self._constants.rp, self._constants.CTi, self._constants.kTi, - self._constants.deltaRTi, self._edgeType, self._spline.basis.cubic_uniform) + self._constants.deltaRTi, self._edgeType) def gridStep(self, grid: Grid, phi: Grid, parGrad: ParallelGradient, parGradVals: np.array, dt: float): for i, r in grid.getCoords(0): @@ -474,41 +470,30 @@ def step(self, f: np.ndarray, dt: float, phi: Spline2D, v: float): assert f.shape == self._nPoints self._interpolator.compute_interpolant(f, self._spline) - phiBases = phi.basis - polBases = self._spline.basis - if (self._explicit): poloidal_advection_step_expl(f, float(dt), v, self._points[1], self._points[0], 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, phiBases[0].knots, - phiBases[1].knots, phi.coeffs, - phiBases[0].degree, phiBases[1].degree, - polBases[0].knots, polBases[1].knots, - self._spline.coeffs, polBases[0].degree, - polBases[1].degree, self._constants.CN0, + self._endPts_k2_r, phi, + 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, phiBases[0].cubic_uniform, self._nulEdge) + self._constants.B0, self._nulEdge) else: poloidal_advection_step_impl(f, float(dt), v, self._points[1], self._points[0], 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, phiBases[0].knots, - phiBases[1].knots, phi.coeffs, - phiBases[0].degree, phiBases[1].degree, - polBases[0].knots, polBases[1].knots, - self._spline.coeffs, polBases[0].degree, - polBases[1].degree, self._constants.CN0, + self._endPts_k2_r, phi.coeffs, + 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, phiBases[0].cubic_uniform, self._nulEdge) + self._constants.B0, self._TOL, self._nulEdge) def exact_step(self, f, endPts, v): assert f.shape == self._nPoints From 0d37c280b63554b3fc232e521cb8d55d5c661c4b Mon Sep 17 00:00:00 2001 From: Emily Bourne Date: Thu, 21 Aug 2025 12:39:18 +0200 Subject: [PATCH 19/30] Split eval and eval_vector cleanly to avoid #1339 --- pygyro/advection/advection.py | 2 +- pygyro/poisson/poisson_solver.py | 24 ++++--- pygyro/splines/splines.py | 71 ++++++------------- .../tests/test_spline_interpolators.py | 33 ++++++--- pygyro/splines/tests/test_splines.py | 9 ++- 5 files changed, 68 insertions(+), 71 deletions(-) diff --git a/pygyro/advection/advection.py b/pygyro/advection/advection.py index 79a0bef6..5a70f893 100644 --- a/pygyro/advection/advection.py +++ b/pygyro/advection/advection.py @@ -488,7 +488,7 @@ def step(self, f: np.ndarray, dt: float, phi: Spline2D, v: float): 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, phi.coeffs, + self._endPts_k2_r, phi, self._spline, self._constants.CN0, self._constants.kN0, self._constants.deltaRN0, self._constants.rp, self._constants.CTi, diff --git a/pygyro/poisson/poisson_solver.py b/pygyro/poisson/poisson_solver.py index 5b8114d1..dd5c75f1 100644 --- a/pygyro/poisson/poisson_solver.py +++ b/pygyro/poisson/poisson_solver.py @@ -254,24 +254,32 @@ def __init__(self, degree: int, rspline: BSplines, nr: int, nTheta: int, # Find the integral of the multiplication of these splines # and their coefficients and save the value in the # appropriate place for the matrix + rEval = np.empty_like(evalPts) + rDeriv = np.empty_like(evalPts) + splEval = np.empty_like(evalPts) + splDeriv = np.empty_like(evalPts) + self._rspline[s_j].eval_vector(evalPts, rEval) + spline.eval_vector(evalPts, splEval) + self._rspline[s_j].eval_vector(evalPts, rDeriv, der=1) + spline.eval_vector(evalPts, splDeriv, der=1) massCoeffs[j][i] = np.sum(np.tile(self._weights, end-start) * multFactor * - rhoFactor(evalPts) * self._rspline[s_j].eval(evalPts) * spline.eval(evalPts) * evalPts) + rhoFactor(evalPts) * rEval * splEval * evalPts) k2PhiPsiCoeffs[j][i] = np.sum(np.tile(self._weights, end-start) * multFactor * - ddThetaFactor(evalPts) * self._rspline[s_j].eval(evalPts) * spline.eval(evalPts) * evalPts) + ddThetaFactor(evalPts) * rEval * splEval * evalPts) PhiPsiCoeffs[j][i] = np.sum(np.tile(self._weights, end-start) * multFactor * - rFactor(evalPts) * self._rspline[s_j].eval(evalPts) * spline.eval(evalPts) * evalPts) + rFactor(evalPts) * rEval * splEval * evalPts) dPhidPsi = np.sum(np.tile(self._weights, end-start) * multFactor * - -ddrFactor(evalPts) * self._rspline[s_j].eval(evalPts, 1) * spline.eval(evalPts, 1) * evalPts) + -ddrFactor(evalPts) * rDeriv * splDeriv * evalPts) dPhidPsiCoeffs[j][i] = dPhidPsi + \ np.sum(np.tile(self._weights, end-start) * multFactor * - -ddrFactor(evalPts) * self._rspline[s_j].eval(evalPts, 1) * spline.eval(evalPts)) + -ddrFactor(evalPts) * rDeriv * splEval) dPhidPsiCoeffs[self._rspline.degree*2-j][i] = dPhidPsi + \ np.sum(np.tile(self._weights, end-start) * multFactor * - -ddrFactor(evalPts) * self._rspline[s_j].eval(evalPts) * spline.eval(evalPts, 1)) + -ddrFactor(evalPts) * rEval * splDeriv) dPhiPsiCoeffs[j][i] = np.sum(np.tile(self._weights, end-start) * multFactor * - drFactor(evalPts) * self._rspline[s_j].eval(evalPts, 1) * spline.eval(evalPts) * evalPts) + drFactor(evalPts) * rDeriv * splEval * evalPts) dPhiPsiCoeffs[self._rspline.degree*2-j][i] = np.sum(np.tile(self._weights, end-start) * multFactor * - drFactor(evalPts) * self._rspline[s_j].eval(evalPts) * spline.eval(evalPts, 1) * evalPts) + drFactor(evalPts) * rEval * splDeriv * evalPts) # Create the diagonal matrices # Diagonal matrices contain many 0 valued points so sparse diff --git a/pygyro/splines/splines.py b/pygyro/splines/splines.py index f78cad0b..9ba0e3eb 100644 --- a/pygyro/splines/splines.py +++ b/pygyro/splines/splines.py @@ -11,8 +11,6 @@ from .cubic_uniform_spline_eval_funcs import cu_eval_spline_1d_scalar, cu_eval_spline_1d_vector from .cubic_uniform_spline_eval_funcs import cu_eval_spline_2d_cross, cu_eval_spline_2d_scalar -ScalOrArr = TypeVar('ScalOrArr', float, 'float[:]') - __all__ = ['make_knots', 'BSplines', 'Spline1D', 'Spline1DComplex', 'Spline2D'] # =============================================================================== @@ -334,25 +332,16 @@ def coeffs(self): """ return self._coeffs - def eval(self, x : 'float|float[:]', der : int = 0): + def eval(self, x : float, der : int = 0): """ TODO """ - if isinstance(x, float): - if self._basis.cubic_uniform: - result = cu_eval_spline_1d_scalar( - x, self._basis.knots, self._basis.degree, self._coeffs, der) - else: - result = nu_eval_spline_1d_scalar( - x, self._basis.knots, self._basis.degree, self._coeffs, der) + if self._basis.cubic_uniform: + result = cu_eval_spline_1d_scalar( + x, self._basis.knots, self._basis.degree, self._coeffs, der) else: - result = np.empty_like(x) - if self._basis.cubic_uniform: - cu_eval_spline_1d_vector(x, self._basis.knots, - self._basis.degree, self._coeffs, result, der) - else: - nu_eval_spline_1d_vector(x, self._basis.knots, - self._basis.degree, self._coeffs, result, der) + result = nu_eval_spline_1d_scalar( + x, self._basis.knots, self._basis.degree, self._coeffs, der) return result """ @@ -395,25 +384,16 @@ def coeffs(self): """ return self._coeffs - def eval(self, x : 'float|float[:]', der : int = 0): + def eval(self, x : float, der : int = 0): """ TODO """ - if isinstance(x, float): - if self._basis.cubic_uniform: - result = cu_eval_spline_1d_scalar( - x, self._basis.knots, self._basis.degree, self._coeffs, der) - else: - result = nu_eval_spline_1d_scalar( - x, self._basis.knots, self._basis.degree, self._coeffs, der) + if self._basis.cubic_uniform: + result = cu_eval_spline_1d_scalar( + x, self._basis.knots, self._basis.degree, self._coeffs, der) else: - result = np.empty_like(x, dtype=np.complex128) - if self._basis.cubic_uniform: - cu_eval_spline_1d_vector(x, self._basis.knots, - self._basis.degree, self._coeffs, result, der) - else: - nu_eval_spline_1d_vector(x, self._basis.knots, - self._basis.degree, self._coeffs, result, der) + result = nu_eval_spline_1d_scalar( + x, self._basis.knots, self._basis.degree, self._coeffs, der) return result """ @@ -473,29 +453,18 @@ def coeffs(self): """ return self._coeffs - def eval(self, x1 : ScalOrArr, x2 : ScalOrArr, der1 : int=0, der2 : int=0): + def eval(self, x1 : float, x2 : float, der1 : int=0, der2 : int=0): """ TODO """ - if isinstance(x1, float): - if self._basis1.cubic_uniform: - result = cu_eval_spline_2d_scalar(x1, x2, self._basis1.knots, self._basis1.degree, - self._basis2.knots, self._basis2.degree, - self._coeffs, der1, der2) - else: - result = nu_eval_spline_2d_scalar(x1, x2, self._basis1.knots, self._basis1.degree, - self._basis2.knots, self._basis2.degree, - self._coeffs, der1, der2) + if self._basis1.cubic_uniform: + result = cu_eval_spline_2d_scalar(x1, x2, self._basis1.knots, self._basis1.degree, + self._basis2.knots, self._basis2.degree, + self._coeffs, der1, der2) else: - result = np.empty((len(x1), len(x2))) - if self._basis1.cubic_uniform: - cu_eval_spline_2d_cross(x1, x2, self._basis1.knots, self._basis1.degree, - self._basis2.knots, self._basis2.degree, - self._coeffs, result, der1, der2) - else: - nu_eval_spline_2d_cross(x1, x2, self._basis1.knots, self._basis1.degree, - self._basis2.knots, self._basis2.degree, - self._coeffs, result, der1, der2) + result = nu_eval_spline_2d_scalar(x1, x2, self._basis1.knots, self._basis1.degree, + self._basis2.knots, self._basis2.degree, + self._coeffs, der1, der2) return result """ diff --git a/pygyro/splines/tests/test_spline_interpolators.py b/pygyro/splines/tests/test_spline_interpolators.py index dfdff888..19a76494 100644 --- a/pygyro/splines/tests/test_spline_interpolators.py +++ b/pygyro/splines/tests/test_spline_interpolators.py @@ -40,8 +40,12 @@ def test_SplineInterpolator1D_exact(ncells, degree): interp.compute_interpolant(ug, spline) xt = np.linspace(*domain, num=100) - err = spline.eval(xt) - poly.eval(xt) - derr = spline.eval(xt, der=1) - poly.eval(xt, diff=1) + vals = np.empty_like(xt) + ders = np.empty_like(xt) + spline.eval_vector(xt, vals) + err = vals - poly.eval(xt) + spline.eval_vector(xt, ders, der=1) + derr = ders - poly.eval(xt, diff=1) max_norm_err = np.max(abs(err)) max_norm_derr = np.max(abs(derr)) @@ -74,8 +78,13 @@ def test_SplineInterpolator1D_exact_uniform(ncells, degree): interp.compute_interpolant(ug, spline) xt = np.linspace(*domain, num=100) - err = spline.eval(xt) - poly.eval(xt) - derr = spline.eval(xt, der=1) - poly.eval(xt, diff=1) + vals = np.empty_like(xt) + spline.eval_vector(xt, vals) + err = vals - poly.eval(xt) + + ders = np.empty_like(xt) + spline.eval_vector(xt, ders, der=1) + derr = ders - poly.eval(xt, diff=1) max_norm_err = np.max(abs(err)) max_norm_derr = np.max(abs(derr)) @@ -118,7 +127,9 @@ def test_SplineInterpolator1D_cosine(ncells, degree, periodic): interp.compute_interpolant(ug, spline) xt = np.linspace(*f.domain, num=100) - err = spline.eval(xt) - f.eval(xt) + vals = np.empty_like(xt) + spline.eval_vector(xt, vals) + err = vals - f.eval(xt) max_norm_err = np.max(abs(err)) err_bound = spline_1d_error_bound(f, np.diff(breaks).max(), degree) @@ -160,7 +171,9 @@ def test_SplineInterpolator1D_uniform_cosine(ncells, periodic): interp.compute_interpolant(ug, spline) xt = np.linspace(*f.domain, num=100) - err = spline.eval(xt) - f.eval(xt) + vals = np.empty_like(xt) + spline.eval_vector(xt, vals) + err = vals - f.eval(xt) max_norm_err = np.max(abs(err)) err_bound = spline_1d_error_bound(f, np.diff(breaks).max(), degree) @@ -213,7 +226,9 @@ def f(x1, x2): return poly.eval(x1-0.5*x2) x1t = np.linspace(*domain1, num=100) x2t = np.linspace(*domain2, num=100) - err = spline.eval(x1t, x2t) - f(*np.meshgrid(x1t, x2t, indexing='ij')) + vals = np.empty((100,100)) + spline.eval_vector(x1t, x2t, vals) + err = vals - f(*np.meshgrid(x1t, x2t, indexing='ij')) max_norm_err = np.max(abs(err)) assert max_norm_err < 2.0e-14 @@ -259,7 +274,9 @@ def test_SplineInterpolator2D_cosine(ncells, degree, periodic1, periodic2): x1t = np.linspace(*domain1, num=20) x2t = np.linspace(*domain2, num=20) - err = spline.eval(x1t, x2t) - f.eval(np.meshgrid(x1t, x2t, indexing='ij')) + vals = np.empty((20,20)) + spline.eval_vector(x1t, x2t, vals) + err = vals - f.eval(np.meshgrid(x1t, x2t, indexing='ij')) max_norm_err = np.max(abs(err)) err_bound = spline_2d_error_bound(f, np.diff( diff --git a/pygyro/splines/tests/test_splines.py b/pygyro/splines/tests/test_splines.py index 979ca75f..77419e9b 100644 --- a/pygyro/splines/tests/test_splines.py +++ b/pygyro/splines/tests/test_splines.py @@ -87,7 +87,8 @@ def test_BSplines(ncells, degree, periodic, npts=50, tol=1e-15): f = np.zeros(npts) # Accumulated values of all basis functions for i in range(ncells+degree): - fi = basis[i].eval(x) # Evaluate basis function at all test points + fi = np.empty_like(x) + basis[i].eval_vector(x, fi) # Evaluate basis function at all test points f += fi # Sum contributions from all basis functions assert all(fi >= 0.0) # Check positivity of each basis function assert all(abs(1.0-f) < tol) # Check partition of unity @@ -109,7 +110,8 @@ def test_Spline1D_unit(ncells, degree, periodic, npts=50, tol=1e-15): spline.coeffs.fill(1.0) x = np.linspace(breaks[0], breaks[-1], npts) # Test points - f = spline.eval(x) + f = np.empty_like(x) + spline.eval_vector(x, f) assert all(abs(1.0-f) < tol) @@ -161,6 +163,7 @@ def test_Spline2D_unit(ncells, degree, periodic, npts=10, tol=1e-15): x1 = np.linspace(breaks1[0], breaks1[-1], npts) # Test points x2 = np.linspace(breaks2[0], breaks2[-1], npts) # Test points - f = spline.eval(x1, x2) + f = np.empty((npts, npts)) + spline.eval_vector(x1, x2, f) assert np.all(abs(1.0-f) < tol) From 45a39612b047404740be234a17e876d5a94a4ec3 Mon Sep 17 00:00:00 2001 From: Emily Bourne Date: Thu, 21 Aug 2025 12:47:39 +0200 Subject: [PATCH 20/30] Correct type --- pygyro/splines/cubic_uniform_spline_eval_funcs.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/pygyro/splines/cubic_uniform_spline_eval_funcs.py b/pygyro/splines/cubic_uniform_spline_eval_funcs.py index 66d6eb26..e1bbf07a 100644 --- a/pygyro/splines/cubic_uniform_spline_eval_funcs.py +++ b/pygyro/splines/cubic_uniform_spline_eval_funcs.py @@ -327,7 +327,7 @@ def cu_eval_spline_2d_vector(x: 'float[:]', y: 'float[:]', kts1: 'float[:]', deg basis1 = empty(4) basis2 = empty(4) - theCoeffs = empty((4, 4)) + theCoeffs = empty((4, 4), dtype = type(z[0])) if (der1 == 0): if (der2 == 0): From 44a49eefc62b3d5d4d8833161c90b7a6a57bef8e Mon Sep 17 00:00:00 2001 From: Emily Bourne Date: Thu, 21 Aug 2025 12:47:54 +0200 Subject: [PATCH 21/30] Back to Fortran --- Makefile | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Makefile b/Makefile index 24d1e1fa..50f7384c 100644 --- a/Makefile +++ b/Makefile @@ -9,7 +9,7 @@ ACC := pycc COMP := GNU # Target language -LANGUAGE := c +LANGUAGE := fortran #PYTHRAN_FLAGS := -DUSE_XSIMD -fopenmp -march=native PYTHRAN_FLAGS := From 62b8ad2f2c992d4db729271112890542e5ff532e Mon Sep 17 00:00:00 2001 From: Emily Bourne Date: Thu, 21 Aug 2025 12:50:33 +0200 Subject: [PATCH 22/30] Correct type --- pygyro/splines/cubic_uniform_spline_eval_funcs.py | 4 ++-- pygyro/splines/spline_eval_funcs.py | 6 +++--- 2 files changed, 5 insertions(+), 5 deletions(-) diff --git a/pygyro/splines/cubic_uniform_spline_eval_funcs.py b/pygyro/splines/cubic_uniform_spline_eval_funcs.py index e1bbf07a..44b4248b 100644 --- a/pygyro/splines/cubic_uniform_spline_eval_funcs.py +++ b/pygyro/splines/cubic_uniform_spline_eval_funcs.py @@ -209,7 +209,7 @@ def cu_eval_spline_2d_scalar(x: 'float', y: 'float', kts1: 'Final[float[:]]', de elif (der2 == 1): cu_basis_funs_1st_der(span2, offset2, dy, basis2) - theCoeffs = empty((4, 4)) + theCoeffs = empty((4, 4), dtype=type(coeffs[0,0])) theCoeffs[:, :] = coeffs[span1-deg1:span1+1, span2-deg2:span2+1] z = 0.0*coeffs[0,0] @@ -235,7 +235,7 @@ def cu_eval_spline_2d_cross(X: 'Final[float[:]]', Y: 'Final[float[:]]', kts1: 'F basis1 = empty(4) basis2 = empty(4) - theCoeffs = empty((4, 4)) + theCoeffs = empty((4, 4), dtype=type(coeffs[0,0])) if (der1 == 0 and der2 == 0): for i, x in enumerate(X): diff --git a/pygyro/splines/spline_eval_funcs.py b/pygyro/splines/spline_eval_funcs.py index a4864b9d..74798805 100644 --- a/pygyro/splines/spline_eval_funcs.py +++ b/pygyro/splines/spline_eval_funcs.py @@ -235,7 +235,7 @@ def nu_eval_spline_2d_scalar(x: 'float', y: 'float', kts1: 'Final[float[:]]', de elif (der2 == 1): nu_basis_funs_1st_der(kts2, deg2, y, span2, basis2) - theCoeffs = empty((deg1+1, deg2+1)) + theCoeffs = empty((deg1+1, deg2+1), dtype=type(coeffs[0,0])) theCoeffs[:, :] = coeffs[span1-deg1:span1+1, span2-deg2:span2+1] z = 0.0*coeffs[0,0] @@ -256,7 +256,7 @@ def nu_eval_spline_2d_cross(X: 'Final[float[:]]', Y: 'Final[float[:]]', kts1: 'F """ basis1 = empty(deg1+1) basis2 = empty(deg2+1) - theCoeffs = empty((deg1+1, deg2+1)) + theCoeffs = empty((deg1+1, deg2+1), dtype=type(z[0,0])) if (der1 == 0 and der2 == 0): for i, x in enumerate(X): @@ -344,7 +344,7 @@ def nu_eval_spline_2d_vector(x: 'float[:]', y: 'float[:]', kts1: 'float[:]', deg """ basis1 = empty(deg1+1) basis2 = empty(deg2+1) - theCoeffs = empty((deg1+1, deg2+1)) + theCoeffs = empty((deg1+1, deg2+1), dtype=type(z[0])) if (der1 == 0): if (der2 == 0): From be7790529e76c7899f24c08e2243c8408a054843 Mon Sep 17 00:00:00 2001 From: Emily Bourne Date: Thu, 21 Aug 2025 13:11:29 +0200 Subject: [PATCH 23/30] Ensure Fortran compilation --- pygyro/splines/splines.py | 23 ++++++++++++----------- 1 file changed, 12 insertions(+), 11 deletions(-) diff --git a/pygyro/splines/splines.py b/pygyro/splines/splines.py index 9ba0e3eb..7b966293 100644 --- a/pygyro/splines/splines.py +++ b/pygyro/splines/splines.py @@ -126,7 +126,7 @@ def __init__(self, knots : 'float[:]', degree : int, periodic : bool, uniform : self._interp_pts = np.empty(self._nbasis) self._interp_pts[0] = xmin self._interp_pts[1] = xmin+dx/3 - self._interp_pts[2:-2] = np.linspace(xmin+dx, xmax-dx, self._nbasis-4) + self._interp_pts[2:-2] = [xmin + dx*i for i in range(1, self._nbasis-3)] self._interp_pts[-2] = xmax-dx/3 self._interp_pts[-1] = xmax #self._interp_pts = np.array([xmin, @@ -179,12 +179,12 @@ def breaks(self): p = self._degree return np.array(self._knots[p:n-p]) - @property - def domain(self): - """ Domain boundaries [a,b]. - """ - breaks = self.breaks - return breaks[0], breaks[-1] + #@property + #def domain(self): + # """ Domain boundaries [a,b]. + # """ + # breaks = self.breaks + # return breaks[0], breaks[-1] @property def cubic_uniform(self): @@ -204,7 +204,8 @@ def greville(self): x = np.array([np.sum(T[i:i+p])/p for i in range(s, s+n)]) if self._periodic: - a, b = self.domain + a = self.breaks[0] + b = self.breaks[-1] #x = np.around(x, decimals=15) x[:] = (x-a) % (b-a) + a @@ -349,7 +350,7 @@ def eval(self, x : float, der : int = 0): return splev( x, tck, der ) """ - def eval_vector(self, x : 'float[:]', y : 'float[:]', der : int=0): + def eval_vector(self, x : 'Final[float[:]]', y : 'float[:]', der : int=0): """ TODO """ @@ -401,7 +402,7 @@ def eval(self, x : float, der : int = 0): return splev( x, tck, der ) """ - def eval_vector(self, x : 'float[:]', y : 'complex[:]', der : int=0): + def eval_vector(self, x : 'Final[float[:]]', y : 'complex[:]', der : int=0): """ TODO """ @@ -478,7 +479,7 @@ def eval(self, x1 : float, x2 : float, der1 : int=0, der2 : int=0): return bisplev( x1, x2, tck, der1, der2 ) """ - def eval_vector(self, x1 : 'float[:]', x2 : 'float[:]', y : 'float[:,:]', der1 : int=0, der2 : int=0): + def eval_vector(self, x1 : 'Final[float[:]]', x2 : 'Final[float[:]]', y : 'float[:,:]', der1 : int=0, der2 : int=0): """ TODO """ From 23b4ac0a0c60073c6ed859fbe4bf847638fa25a7 Mon Sep 17 00:00:00 2001 From: Emily Bourne Date: Sat, 20 Sep 2025 20:04:27 +0200 Subject: [PATCH 24/30] Autopep --- .../advection/accelerated_advection_steps.py | 60 ++++++++++--------- .../cubic_uniform_spline_eval_funcs.py | 9 +-- pygyro/splines/spline_eval_funcs.py | 7 ++- pygyro/splines/spline_interpolators.py | 6 +- pygyro/splines/splines.py | 46 +++++++------- .../tests/test_spline_interpolators.py | 4 +- pygyro/splines/tests/test_splines.py | 9 +-- 7 files changed, 77 insertions(+), 64 deletions(-) diff --git a/pygyro/advection/accelerated_advection_steps.py b/pygyro/advection/accelerated_advection_steps.py index 5cd5dd3f..4800543a 100644 --- a/pygyro/advection/accelerated_advection_steps.py +++ b/pygyro/advection/accelerated_advection_steps.py @@ -5,14 +5,14 @@ def poloidal_advection_step_expl(f: 'float[:,:]', - dt: 'float', v: 'float', - rPts: 'Final[float[:]]', qPts: 'Final[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[:,:]', - phi_spline : Spline2D, pol_spline : Spline2D, - CN0: 'float', kN0: 'float', deltaRN0: 'float', rp: 'float', CTi: 'float', - kTi: 'float', deltaRTi: 'float', B0: 'float', nulBound: 'bool'): + dt: 'float', v: 'float', + rPts: 'Final[float[:]]', qPts: 'Final[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[:,:]', + phi_spline: Spline2D, pol_spline: Spline2D, + CN0: 'float', kN0: 'float', deltaRN0: 'float', rp: 'float', CTi: 'float', + kTi: 'float', deltaRTi: 'float', B0: 'float', nulBound: 'bool'): """ Carry out an advection step for the poloidal advection @@ -64,11 +64,11 @@ def poloidal_advection_step_expl(f: 'float[:,:]', # x^{n+1} = x^n + 0.5( f(x^n) + f(x^n + f(x^n)) ) # ^^^^^^^^^^^^^^^ drPhi_k[i, j] = phi_spline.eval(endPts_k1_q[i, j], endPts_k1_r[i, j], - 0, 1) + 0, 1) drPhi_k[i, j] /= endPts_k1_r[i, j] dthetaPhi_k[i, j] = phi_spline.eval(endPts_k1_q[i, j], endPts_k1_r[i, j], - 1, 0) + 1, 0) dthetaPhi_k[i, j] /= endPts_k1_r[i, j] else: drPhi_k[i, j] = 0.0 @@ -92,7 +92,8 @@ def poloidal_advection_step_expl(f: 'float[:,:]', f[i, j] = 0.0 else: 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]) + f[i, j] = pol_spline.eval( + endPts_k2_q[i, j], endPts_k2_r[i, j]) else: for i in range(nPts_q): # theta for j in range(nPts_r): # r @@ -104,14 +105,15 @@ def poloidal_advection_step_expl(f: 'float[:,:]', deltaRN0, rp, CTi, kTi, deltaRTi) else: 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]) + f[i, j] = pol_spline.eval( + endPts_k2_q[i, j], endPts_k2_r[i, j]) def v_parallel_advection_eval_step(f: 'float[:]', vPts: 'float[:]', - rPos: 'float', vMin: 'float', vMax: 'float', - spl : Spline1D, - CN0: 'float', kN0: 'float', deltaRN0: 'float', rp: 'float', - CTi: 'float', kTi: 'float', deltaRTi: 'float', bound: 'int'): + rPos: 'float', vMin: 'float', vMax: 'float', + spl: Spline1D, + CN0: 'float', kN0: 'float', deltaRN0: 'float', rp: 'float', + CTi: 'float', kTi: 'float', deltaRTi: 'float', bound: 'int'): """ TODO """ @@ -140,8 +142,8 @@ def v_parallel_advection_eval_step(f: 'float[:]', vPts: 'float[:]', def get_lagrange_vals(i: 'int', shifts: 'int[:]', - vals: 'float[:,:,:]', qVals: 'float[:]', - thetaShifts: 'float[:]', spl : Spline1D): + vals: 'float[:,:,:]', qVals: 'float[:]', + thetaShifts: 'float[:]', spl: Spline1D): """ TODO """ @@ -169,11 +171,11 @@ def flux_advection(nq: 'int', nr: 'int', 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[:,:]', - 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'): + 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'): """ Carry out an advection step for the poloidal advection @@ -230,9 +232,11 @@ def poloidal_advection_step_impl(f: 'float[:,:]', dt: 'float', v: 'float', rPts: # Add the new value of phi to the derivatives # x^{n+1} = x^n + 0.5( f(x^n) + f(x^n + f(x^n)) ) # ^^^^^^^^^^^^^^^ - drPhi_k[i, j] = phi_spline.eval(endPts_k1_q[i, j], endPts_k1_r[i, j], 0, 1) + drPhi_k[i, j] = phi_spline.eval( + endPts_k1_q[i, j], endPts_k1_r[i, j], 0, 1) drPhi_k[i, j] /= endPts_k1_r[i, j] - dthetaPhi_k[i, j] = phi_spline.eval(endPts_k1_q[i, j], endPts_k1_r[i, j], 1, 0) + dthetaPhi_k[i, j] = phi_spline.eval( + endPts_k1_q[i, j], endPts_k1_r[i, j], 1, 0) dthetaPhi_k[i, j] /= endPts_k1_r[i, j] else: drPhi_k[i, j] = 0.0 @@ -275,7 +279,8 @@ def poloidal_advection_step_impl(f: 'float[:,:]', dt: 'float', v: 'float', rPts: f[i, j] = 0.0 else: 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]) + f[i, j] = pol_spline.eval( + endPts_k2_q[i, j], endPts_k2_r[i, j]) else: for i in range(nPts_q): for j in range(nPts_r): @@ -287,4 +292,5 @@ def poloidal_advection_step_impl(f: 'float[:,:]', dt: 'float', v: 'float', rPts: deltaRN0, rp, CTi, kTi, deltaRTi) else: 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]) + f[i, j] = pol_spline.eval( + endPts_k2_q[i, j], endPts_k2_r[i, j]) diff --git a/pygyro/splines/cubic_uniform_spline_eval_funcs.py b/pygyro/splines/cubic_uniform_spline_eval_funcs.py index 44b4248b..70b622ba 100644 --- a/pygyro/splines/cubic_uniform_spline_eval_funcs.py +++ b/pygyro/splines/cubic_uniform_spline_eval_funcs.py @@ -4,6 +4,7 @@ CoeffType = TypeVar('CoeffType', float, complex) + @pure def cu_find_span(xmin: 'float', xmax: 'float', dx: 'float', x: 'float', ncells: 'int'): """ @@ -209,10 +210,10 @@ def cu_eval_spline_2d_scalar(x: 'float', y: 'float', kts1: 'Final[float[:]]', de elif (der2 == 1): cu_basis_funs_1st_der(span2, offset2, dy, basis2) - theCoeffs = empty((4, 4), dtype=type(coeffs[0,0])) + theCoeffs = empty((4, 4), dtype=type(coeffs[0, 0])) theCoeffs[:, :] = coeffs[span1-deg1:span1+1, span2-deg2:span2+1] - z = 0.0*coeffs[0,0] + z = 0.0*coeffs[0, 0] for i in range(4): theCoeffs[i, 0] = theCoeffs[i, 0]*basis2[0] for j in range(1, 4): @@ -235,7 +236,7 @@ def cu_eval_spline_2d_cross(X: 'Final[float[:]]', Y: 'Final[float[:]]', kts1: 'F basis1 = empty(4) basis2 = empty(4) - theCoeffs = empty((4, 4), dtype=type(coeffs[0,0])) + theCoeffs = empty((4, 4), dtype=type(coeffs[0, 0])) if (der1 == 0 and der2 == 0): for i, x in enumerate(X): @@ -327,7 +328,7 @@ def cu_eval_spline_2d_vector(x: 'float[:]', y: 'float[:]', kts1: 'float[:]', deg basis1 = empty(4) basis2 = empty(4) - theCoeffs = empty((4, 4), dtype = type(z[0])) + theCoeffs = empty((4, 4), dtype=type(z[0])) if (der1 == 0): if (der2 == 0): diff --git a/pygyro/splines/spline_eval_funcs.py b/pygyro/splines/spline_eval_funcs.py index 74798805..e3533f59 100644 --- a/pygyro/splines/spline_eval_funcs.py +++ b/pygyro/splines/spline_eval_funcs.py @@ -4,6 +4,7 @@ CoeffType = TypeVar('CoeffType', float, complex) + @pure def nu_find_span(knots: 'Final[float[:]]', degree: 'int', x: 'float') -> int: """ @@ -235,10 +236,10 @@ def nu_eval_spline_2d_scalar(x: 'float', y: 'float', kts1: 'Final[float[:]]', de elif (der2 == 1): nu_basis_funs_1st_der(kts2, deg2, y, span2, basis2) - theCoeffs = empty((deg1+1, deg2+1), dtype=type(coeffs[0,0])) + theCoeffs = empty((deg1+1, deg2+1), dtype=type(coeffs[0, 0])) theCoeffs[:, :] = coeffs[span1-deg1:span1+1, span2-deg2:span2+1] - z = 0.0*coeffs[0,0] + z = 0.0*coeffs[0, 0] for i in range(deg1+1): theCoeffs[i, 0] = theCoeffs[i, 0]*basis2[0] for j in range(1, deg2+1): @@ -256,7 +257,7 @@ def nu_eval_spline_2d_cross(X: 'Final[float[:]]', Y: 'Final[float[:]]', kts1: 'F """ basis1 = empty(deg1+1) basis2 = empty(deg2+1) - theCoeffs = empty((deg1+1, deg2+1), dtype=type(z[0,0])) + theCoeffs = empty((deg1+1, deg2+1), dtype=type(z[0, 0])) if (der1 == 0 and der2 == 0): for i, x in enumerate(X): diff --git a/pygyro/splines/spline_interpolators.py b/pygyro/splines/spline_interpolators.py index effe1137..d9e278c7 100644 --- a/pygyro/splines/spline_interpolators.py +++ b/pygyro/splines/spline_interpolators.py @@ -65,7 +65,7 @@ def compute_interpolant(self, ug, spl): """ assert isinstance(spl, (Spline1D, Spline1DComplex)) - #assert spl.basis is self._basis + # assert spl.basis is self._basis assert len(ug) == self._basis.nbasis if self._basis.periodic: @@ -223,8 +223,8 @@ def compute_interpolant(self, ug, spl): assert isinstance(spl, Spline2D) basis1 = spl.basis1 basis2 = spl.basis2 - #assert basis1 is self._basis1 - #assert basis2 is self._basis2 + # assert basis1 is self._basis1 + # assert basis2 is self._basis2 n1, n2 = basis1.nbasis, basis2.nbasis p1, p2 = basis1.degree, basis2.degree diff --git a/pygyro/splines/splines.py b/pygyro/splines/splines.py index 7b966293..c35c191d 100644 --- a/pygyro/splines/splines.py +++ b/pygyro/splines/splines.py @@ -15,8 +15,9 @@ # =============================================================================== + @allow_negative_index('breaks', 'T') -def make_knots(breaks : 'Final[float[:]]', degree : Final[int], periodic : Final[bool]): +def make_knots(breaks: 'Final[float[:]]', degree: Final[int], periodic: Final[bool]): """ Create spline knots from breakpoints, with appropriate boundary conditions. Let p be spline degree. If domain is periodic, knot sequence is extended @@ -47,7 +48,7 @@ def make_knots(breaks : 'Final[float[:]]', degree : Final[int], periodic : Final # Consistency checks assert len(breaks) > 1 - #assert all(np.diff(breaks) > 0) + # assert all(np.diff(breaks) > 0) assert degree > 0 if periodic: assert len(breaks) > degree @@ -96,7 +97,7 @@ class BSplines(): """ @allow_negative_index('knots') - def __init__(self, knots : 'float[:]', degree : int, periodic : bool, uniform : bool): + def __init__(self, knots: 'float[:]', degree: int, periodic: bool, uniform: bool): xmin = knots[degree] xmax = knots[-degree-1] dx = knots[degree+1]-knots[degree] @@ -126,10 +127,11 @@ def __init__(self, knots : 'float[:]', degree : int, periodic : bool, uniform : self._interp_pts = np.empty(self._nbasis) self._interp_pts[0] = xmin self._interp_pts[1] = xmin+dx/3 - self._interp_pts[2:-2] = [xmin + dx*i for i in range(1, self._nbasis-3)] + self._interp_pts[2:-2] = [xmin + dx * + i for i in range(1, self._nbasis-3)] self._interp_pts[-2] = xmax-dx/3 self._interp_pts[-1] = xmax - #self._interp_pts = np.array([xmin, + # self._interp_pts = np.array([xmin, # xmin+dx/3, # *np.linspace(xmin+dx, xmax-dx, self._nbasis-4), # xmax-dx/3, @@ -179,8 +181,8 @@ def breaks(self): p = self._degree return np.array(self._knots[p:n-p]) - #@property - #def domain(self): + # @property + # def domain(self): # """ Domain boundaries [a,b]. # """ # breaks = self.breaks @@ -206,10 +208,10 @@ def greville(self): if self._periodic: a = self.breaks[0] b = self.breaks[-1] - #x = np.around(x, decimals=15) + # x = np.around(x, decimals=15) x[:] = (x-a) % (b-a) + a - #return np.around(x, decimals=15) + # return np.around(x, decimals=15) return x @property @@ -217,7 +219,7 @@ def integrals(self): return self._integrals # ... - def __getitem__(self, i : int): + def __getitem__(self, i: int): """ Get the i-th basis function as a 1D spline. @@ -271,7 +273,7 @@ def _build_integrals(self): knots[0] = self.knots[0] knots[1:-1] = self.knots knots[-1] = self.knots[-1] - #knots = np.array([self.knots[0], *self.knots, self.knots[-1]]) + # knots = np.array([self.knots[0], *self.knots, self.knots[-1]]) values = np.empty(d+2) for i in range(n): @@ -314,7 +316,7 @@ class Spline1D(): TODO """ - def __init__(self, basis : BSplines): + def __init__(self, basis: BSplines): assert isinstance(basis, BSplines) self._basis = basis self._coeffs = np.zeros(basis.ncells + basis.degree, dtype=float) @@ -333,7 +335,7 @@ def coeffs(self): """ return self._coeffs - def eval(self, x : float, der : int = 0): + def eval(self, x: float, der: int = 0): """ TODO """ @@ -350,7 +352,7 @@ def eval(self, x : float, der : int = 0): return splev( x, tck, der ) """ - def eval_vector(self, x : 'Final[float[:]]', y : 'float[:]', der : int=0): + def eval_vector(self, x: 'Final[float[:]]', y: 'float[:]', der: int = 0): """ TODO """ @@ -361,15 +363,17 @@ def eval_vector(self, x : 'Final[float[:]]', y : 'float[:]', der : int=0): nu_eval_spline_1d_vector(x, self._basis.knots, self._basis.degree, self._coeffs, y, der) + class Spline1DComplex(): """ TODO """ - def __init__(self, basis : BSplines): + def __init__(self, basis: BSplines): assert isinstance(basis, BSplines) self._basis = basis - self._coeffs = np.zeros(basis.ncells + basis.degree, dtype=np.complex128) + self._coeffs = np.zeros( + basis.ncells + basis.degree, dtype=np.complex128) @property def basis(self): @@ -385,7 +389,7 @@ def coeffs(self): """ return self._coeffs - def eval(self, x : float, der : int = 0): + def eval(self, x: float, der: int = 0): """ TODO """ @@ -402,7 +406,7 @@ def eval(self, x : float, der : int = 0): return splev( x, tck, der ) """ - def eval_vector(self, x : 'Final[float[:]]', y : 'complex[:]', der : int=0): + def eval_vector(self, x: 'Final[float[:]]', y: 'complex[:]', der: int = 0): """ TODO """ @@ -421,7 +425,7 @@ class Spline2D(): TODO """ - def __init__(self, basis1 : BSplines, basis2 : BSplines): + def __init__(self, basis1: BSplines, basis2: BSplines): assert isinstance(basis1, BSplines) assert isinstance(basis2, BSplines) shape = (basis1.ncells + basis1.degree, basis2.ncells + basis2.degree) @@ -454,7 +458,7 @@ def coeffs(self): """ return self._coeffs - def eval(self, x1 : float, x2 : float, der1 : int=0, der2 : int=0): + def eval(self, x1: float, x2: float, der1: int = 0, der2: int = 0): """ TODO """ @@ -479,7 +483,7 @@ def eval(self, x1 : float, x2 : float, der1 : int=0, der2 : int=0): return bisplev( x1, x2, tck, der1, der2 ) """ - def eval_vector(self, x1 : 'Final[float[:]]', x2 : 'Final[float[:]]', y : 'float[:,:]', der1 : int=0, der2 : int=0): + def eval_vector(self, x1: 'Final[float[:]]', x2: 'Final[float[:]]', y: 'float[:,:]', der1: int = 0, der2: int = 0): """ TODO """ diff --git a/pygyro/splines/tests/test_spline_interpolators.py b/pygyro/splines/tests/test_spline_interpolators.py index 19a76494..3b023324 100644 --- a/pygyro/splines/tests/test_spline_interpolators.py +++ b/pygyro/splines/tests/test_spline_interpolators.py @@ -226,7 +226,7 @@ def f(x1, x2): return poly.eval(x1-0.5*x2) x1t = np.linspace(*domain1, num=100) x2t = np.linspace(*domain2, num=100) - vals = np.empty((100,100)) + vals = np.empty((100, 100)) spline.eval_vector(x1t, x2t, vals) err = vals - f(*np.meshgrid(x1t, x2t, indexing='ij')) @@ -274,7 +274,7 @@ def test_SplineInterpolator2D_cosine(ncells, degree, periodic1, periodic2): x1t = np.linspace(*domain1, num=20) x2t = np.linspace(*domain2, num=20) - vals = np.empty((20,20)) + vals = np.empty((20, 20)) spline.eval_vector(x1t, x2t, vals) err = vals - f.eval(np.meshgrid(x1t, x2t, indexing='ij')) diff --git a/pygyro/splines/tests/test_splines.py b/pygyro/splines/tests/test_splines.py index 77419e9b..ce97ff15 100644 --- a/pygyro/splines/tests/test_splines.py +++ b/pygyro/splines/tests/test_splines.py @@ -31,9 +31,9 @@ def test_make_knots_periodic(ncells, degree): # =============================================================================== -#@pytest.mark.serial -#@pytest.mark.parametrize("ncells", [1, 5, 10, 23]) -#def test_make_knots_periodic_should_fail(ncells): +# @pytest.mark.serial +# @pytest.mark.parametrize("ncells", [1, 5, 10, 23]) +# def test_make_knots_periodic_should_fail(ncells): # """ # TODO # """ @@ -88,7 +88,8 @@ def test_BSplines(ncells, degree, periodic, npts=50, tol=1e-15): for i in range(ncells+degree): fi = np.empty_like(x) - basis[i].eval_vector(x, fi) # Evaluate basis function at all test points + # Evaluate basis function at all test points + basis[i].eval_vector(x, fi) f += fi # Sum contributions from all basis functions assert all(fi >= 0.0) # Check positivity of each basis function assert all(abs(1.0-f) < tol) # Check partition of unity From 8569b5fe4eb7a0b2e37b1c8b3cb97bc1c18cdfdf Mon Sep 17 00:00:00 2001 From: Emily Bourne Date: Sat, 20 Sep 2025 20:12:12 +0200 Subject: [PATCH 25/30] Use devel branch --- requirements.txt | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/requirements.txt b/requirements.txt index 801c56a9..6db0acc3 100644 --- a/requirements.txt +++ b/requirements.txt @@ -4,7 +4,8 @@ scipy>=1.1.0 pytest>=2.8.7 mpi4py>=3.0.0 h5py>=2.8.0 -pyccel>=2.0.0 +git+https://github.com/pyccel/pyccel.git@devel +#pyccel>=2.0.0 # h5py must be built from source using MPI compiler # and linked to parallel HDF5 library. To do so set From 638d1518196f8253c526a7df2d4572e85f2f6d7b Mon Sep 17 00:00:00 2001 From: Emily Bourne Date: Sat, 20 Sep 2025 20:15:11 +0200 Subject: [PATCH 26/30] Unused variable --- pygyro/splines/splines.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/pygyro/splines/splines.py b/pygyro/splines/splines.py index c35c191d..e9f678d4 100644 --- a/pygyro/splines/splines.py +++ b/pygyro/splines/splines.py @@ -1,7 +1,7 @@ # coding: utf-8 # Copyright 2018 Yaman Güçlü from pyccel.decorators import allow_negative_index -from typing import TypeVar, Final +from typing import Final import numpy as np # from scipy.interpolate import splev, bisplev From 6cafcd2638f6602b64400707f0652658035d7a8e Mon Sep 17 00:00:00 2001 From: Emily Bourne Date: Sat, 20 Sep 2025 20:22:17 +0200 Subject: [PATCH 27/30] Test with supercomputing branch --- requirements.txt | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/requirements.txt b/requirements.txt index 6db0acc3..b0dbd305 100644 --- a/requirements.txt +++ b/requirements.txt @@ -4,7 +4,7 @@ scipy>=1.1.0 pytest>=2.8.7 mpi4py>=3.0.0 h5py>=2.8.0 -git+https://github.com/pyccel/pyccel.git@devel +git+https://github.com/pyccel/pyccel.git@supercomputing #pyccel>=2.0.0 # h5py must be built from source using MPI compiler From 809e8a9e60093b604a4cc971f7b839492658f451 Mon Sep 17 00:00:00 2001 From: Emily Bourne Date: Wed, 10 Sep 2025 17:28:18 +0200 Subject: [PATCH 28/30] Inline properties and make eval methods pure --- pygyro/splines/splines.py | 25 +++++++++++++++++++++---- 1 file changed, 21 insertions(+), 4 deletions(-) diff --git a/pygyro/splines/splines.py b/pygyro/splines/splines.py index e9f678d4..c36d42ba 100644 --- a/pygyro/splines/splines.py +++ b/pygyro/splines/splines.py @@ -1,6 +1,6 @@ # coding: utf-8 # Copyright 2018 Yaman Güçlü -from pyccel.decorators import allow_negative_index +from pyccel.decorators import allow_negative_index, inline, pure from typing import Final import numpy as np @@ -139,30 +139,35 @@ def __init__(self, knots: 'float[:]', degree: int, periodic: bool, uniform: bool else: self._interp_pts = np.empty(0) + @inline @property def degree(self): """ Degree of B-splines. """ return self._degree + @inline @property def ncells(self): """ Number of cells in domain. """ return self._ncells + @inline @property def nbasis(self): """ Number of basis functions, taking into account periodicity. """ return self._nbasis + @inline @property def periodic(self): """ True if domain is periodic, False otherwise. """ return self._periodic + @inline @property def knots(self): """ Knot sequence. @@ -188,6 +193,7 @@ def breaks(self): # breaks = self.breaks # return breaks[0], breaks[-1] + @inline @property def cubic_uniform(self): return self._cubic_uniform_splines @@ -214,6 +220,7 @@ def greville(self): # return np.around(x, decimals=15) return x + @inline @property def integrals(self): return self._integrals @@ -321,6 +328,7 @@ def __init__(self, basis: BSplines): self._basis = basis self._coeffs = np.zeros(basis.ncells + basis.degree, dtype=float) + @inline @property def basis(self): """ @@ -328,6 +336,7 @@ def basis(self): """ return self._basis + @inline @property def coeffs(self): """ @@ -335,7 +344,8 @@ def coeffs(self): """ return self._coeffs - def eval(self, x: float, der: int = 0): + @pure + def eval(self : 'Final[Spline1D]', x : float, der : int = 0): """ TODO """ @@ -375,6 +385,7 @@ def __init__(self, basis: BSplines): self._coeffs = np.zeros( basis.ncells + basis.degree, dtype=np.complex128) + @inline @property def basis(self): """ @@ -382,6 +393,7 @@ def basis(self): """ return self._basis + @inline @property def coeffs(self): """ @@ -389,7 +401,8 @@ def coeffs(self): """ return self._coeffs - def eval(self, x: float, der: int = 0): + @pure + def eval(self : 'Final[Spline1DComplex]', x : float, der : int = 0): """ TODO """ @@ -437,6 +450,7 @@ def __init__(self, basis1: BSplines, basis2: BSplines): assert basis2.degree <= 5 assert basis1.cubic_uniform == basis2.cubic_uniform + @inline @property def basis1(self): """ @@ -444,6 +458,7 @@ def basis1(self): """ return self._basis1 + @inline @property def basis2(self): """ @@ -451,6 +466,7 @@ def basis2(self): """ return self._basis2 + @inline @property def coeffs(self): """ @@ -458,7 +474,8 @@ def coeffs(self): """ return self._coeffs - def eval(self, x1: float, x2: float, der1: int = 0, der2: int = 0): + @pure + def eval(self : 'Final[Spline2D]', x1 : float, x2 : float, der1 : int=0, der2 : int=0): """ TODO """ From a0f8b10ae484bb56490429f7bf2770664502cf5d Mon Sep 17 00:00:00 2001 From: Emily Bourne Date: Wed, 24 Sep 2025 19:26:46 +0200 Subject: [PATCH 29/30] Pep --- pygyro/splines/splines.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/pygyro/splines/splines.py b/pygyro/splines/splines.py index c36d42ba..87491531 100644 --- a/pygyro/splines/splines.py +++ b/pygyro/splines/splines.py @@ -345,7 +345,7 @@ def coeffs(self): return self._coeffs @pure - def eval(self : 'Final[Spline1D]', x : float, der : int = 0): + def eval(self: 'Final[Spline1D]', x: float, der: int = 0): """ TODO """ @@ -402,7 +402,7 @@ def coeffs(self): return self._coeffs @pure - def eval(self : 'Final[Spline1DComplex]', x : float, der : int = 0): + def eval(self: 'Final[Spline1DComplex]', x: float, der: int = 0): """ TODO """ @@ -475,7 +475,7 @@ def coeffs(self): return self._coeffs @pure - def eval(self : 'Final[Spline2D]', x1 : float, x2 : float, der1 : int=0, der2 : int=0): + def eval(self: 'Final[Spline2D]', x1: float, x2: float, der1: int = 0, der2: int = 0): """ TODO """ From 2a6becfcaae29f1254e551ee53c005a237c7d232 Mon Sep 17 00:00:00 2001 From: Emily Bourne Date: Fri, 14 Nov 2025 18:15:11 +0100 Subject: [PATCH 30/30] Test with devel branch --- requirements.txt | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/requirements.txt b/requirements.txt index b0dbd305..260fe245 100644 --- a/requirements.txt +++ b/requirements.txt @@ -4,8 +4,8 @@ scipy>=1.1.0 pytest>=2.8.7 mpi4py>=3.0.0 h5py>=2.8.0 -git+https://github.com/pyccel/pyccel.git@supercomputing -#pyccel>=2.0.0 +git+https://github.com/pyccel/pyccel.git@devel +#pyccel>=2.1.0 # h5py must be built from source using MPI compiler # and linked to parallel HDF5 library. To do so set