From a4f734d7ca1da4bd15091a80a831cd9738cf4af8 Mon Sep 17 00:00:00 2001 From: Emily Bourne Date: Wed, 24 Sep 2025 22:00:00 +0200 Subject: [PATCH 01/15] Add OpenMP parallelism to Poisson --- pygyro/poisson/poisson_tools.py | 2 ++ 1 file changed, 2 insertions(+) diff --git a/pygyro/poisson/poisson_tools.py b/pygyro/poisson/poisson_tools.py index d5309cf0..6634a8c1 100644 --- a/pygyro/poisson/poisson_tools.py +++ b/pygyro/poisson/poisson_tools.py @@ -25,6 +25,7 @@ def get_perturbed_rho(rho: T, feq: 'float[:,:]', grid: 'float[:,:,:,:]', nc, = quad_coeffs.shape + #$omp parallel for collapse(3) for i in range(n): for j in range(m): for k in range(p): @@ -53,6 +54,7 @@ def get_rho(rho: T, grid: 'float[:,:,:,:]', quad_coeffs: 'float[:]'): nc, = quad_coeffs.shape + #$omp parallel for collapse(3) for i in range(n): for j in range(m): for k in range(p): From 05997cfe0bec588aae44c994ffca37bd1df470d5 Mon Sep 17 00:00:00 2001 From: Emily Bourne Date: Wed, 24 Sep 2025 22:42:27 +0200 Subject: [PATCH 02/15] Autopep --- pygyro/poisson/poisson_tools.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/pygyro/poisson/poisson_tools.py b/pygyro/poisson/poisson_tools.py index 6634a8c1..480ff689 100644 --- a/pygyro/poisson/poisson_tools.py +++ b/pygyro/poisson/poisson_tools.py @@ -25,7 +25,7 @@ def get_perturbed_rho(rho: T, feq: 'float[:,:]', grid: 'float[:,:,:,:]', nc, = quad_coeffs.shape - #$omp parallel for collapse(3) + # $omp parallel for collapse(3) for i in range(n): for j in range(m): for k in range(p): @@ -54,7 +54,7 @@ def get_rho(rho: T, grid: 'float[:,:,:,:]', quad_coeffs: 'float[:]'): nc, = quad_coeffs.shape - #$omp parallel for collapse(3) + # $omp parallel for collapse(3) for i in range(n): for j in range(m): for k in range(p): From a9db0103e3e9c732e4368de2954125f71344acf9 Mon Sep 17 00:00:00 2001 From: Emily Bourne Date: Wed, 24 Sep 2025 23:34:59 +0200 Subject: [PATCH 03/15] Compile with OpenMP --- Makefile | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Makefile b/Makefile index 50f7384c..e0d5fa02 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-family=$(COMP) --flags ' $(FC_FLAGS)' --language=$(LANGUAGE) + TOOL_FLAGS := --compiler-family=$(COMP) --flags ' $(FC_FLAGS)' --language=$(LANGUAGE) --openmp NAME_PREFIX := else ifeq ($(ACC), numba) From 9fda0dc69d38c645d8917f21aa976831dab8a69f Mon Sep 17 00:00:00 2001 From: Emily Bourne Date: Wed, 24 Sep 2025 23:37:49 +0200 Subject: [PATCH 04/15] Ignore comment space for OpenMP --- .github/workflows/testing.yml | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/.github/workflows/testing.yml b/.github/workflows/testing.yml index 78c7631d..28b7752f 100644 --- a/.github/workflows/testing.yml +++ b/.github/workflows/testing.yml @@ -198,15 +198,15 @@ jobs: python -m pip install autopep8 - name: Check for pep changes run: | - if [[ $(autopep8 -r --diff pygyro/) ]] + if [[ $(autopep8 --ignore=E265 -r --diff pygyro/) ]] then - echo "Please run 'autopep8 -ir pygyro/' to clean up formatting" - autopep8 -r --diff pygyro/ + echo "Please run 'autopep8 --ignore=E265 -ir pygyro/' to clean up formatting" + autopep8 --ignore=E265 -r --diff pygyro/ exit 1 fi - if [[ $(autopep8 -r --diff plotting/) ]] + if [[ $(autopep8 --ignore=E265 -r --diff plotting/) ]] then - echo "Please run 'autopep8 -ir plotting/' to clean up formatting" - autopep8 -r --diff plotting/ + echo "Please run 'autopep8 --ignore=E265 -ir plotting/' to clean up formatting" + autopep8 --ignore=E265 -r --diff plotting/ exit 1 fi From 36aa42f4e462fe02ab9620beb39cde898fa9eccf Mon Sep 17 00:00:00 2001 From: Emily Bourne Date: Wed, 24 Sep 2025 23:43:44 +0200 Subject: [PATCH 05/15] Bad spaces --- pygyro/poisson/poisson_tools.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/pygyro/poisson/poisson_tools.py b/pygyro/poisson/poisson_tools.py index 480ff689..6634a8c1 100644 --- a/pygyro/poisson/poisson_tools.py +++ b/pygyro/poisson/poisson_tools.py @@ -25,7 +25,7 @@ def get_perturbed_rho(rho: T, feq: 'float[:,:]', grid: 'float[:,:,:,:]', nc, = quad_coeffs.shape - # $omp parallel for collapse(3) + #$omp parallel for collapse(3) for i in range(n): for j in range(m): for k in range(p): @@ -54,7 +54,7 @@ def get_rho(rho: T, grid: 'float[:,:,:,:]', quad_coeffs: 'float[:]'): nc, = quad_coeffs.shape - # $omp parallel for collapse(3) + #$omp parallel for collapse(3) for i in range(n): for j in range(m): for k in range(p): From b542dbfe014c35ba4aa8a66db0d84ae59625f299 Mon Sep 17 00:00:00 2001 From: Emily Bourne Date: Wed, 24 Sep 2025 23:57:46 +0200 Subject: [PATCH 06/15] OpenMP on advection --- .../advection/accelerated_advection_steps.py | 20 ++++++++++++++++--- 1 file changed, 17 insertions(+), 3 deletions(-) diff --git a/pygyro/advection/accelerated_advection_steps.py b/pygyro/advection/accelerated_advection_steps.py index 10e5b503..c540afe0 100644 --- a/pygyro/advection/accelerated_advection_steps.py +++ b/pygyro/advection/accelerated_advection_steps.py @@ -41,15 +41,20 @@ def poloidal_advection_step_expl(f: 'float[:,:]', multFactor = dt / B0 multFactor_half = 0.5 * multFactor - phi_spline.eval_vector(qPts, rPts, drPhi_0, 0, 1) - phi_spline.eval_vector(qPts, rPts, dthetaPhi_0, 1, 0) - nPts_r = rPts.shape[0] nPts_q = qPts.shape[0] + #$omp parallel for collapse(2) + for i in range(nPts_r): + for j, q in enumerate(qPts): + r = rPts[i] + drPhi_0[i, j] = phi_spline.eval(q, r, 0, 1) + dthetaPhi_0[i, j] = phi_spline.eval(q, r, 1, 0) + idx = nPts_r-1 rMax = rPts[idx] + #$omp parallel for collapse(2) for i in range(nPts_q): for j in range(nPts_r): # Step one of Heun method @@ -88,6 +93,7 @@ def poloidal_advection_step_expl(f: 'float[:,:]', # Find value at the determined point if (nulBound): + #$omp parallel for collapse(2) for i in range(nPts_q): # theta for j in range(nPts_r): # r if (endPts_k2_r[i, j] < rPts[0]): @@ -99,6 +105,7 @@ def poloidal_advection_step_expl(f: 'float[:,:]', f[i, j] = pol_spline.eval( endPts_k2_q[i, j], endPts_k2_r[i, j]) else: + #$omp parallel for collapse(2) for i in range(nPts_q): # theta for j in range(nPts_r): # r if (endPts_k2_r[i, j] < rPts[0]): @@ -152,6 +159,7 @@ def v_parallel_advection_eval_step_loop(f: 'float[:,:,:]', vPts: 'float[:]', CN0: 'float', kN0: 'float', deltaRN0: 'float', rp: 'float', CTi: 'float', kTi: 'float', deltaRTi: 'float', bound: 'int'): n1, n2, _ = f.shape + #$omp parallel for collapse(2) firstprivate(spl, vPts) private(coeffs) for j in range(n1): # z for k in range(n2): # q coeffs = spl.coeffs @@ -197,6 +205,7 @@ def flux_advection_loop(f: 'float[:,:,:,:]', thetaSpline: Spline1D, theta_offset thetaShifts: 'float[:,:,:]', lagrange_coeffs: 'float[:,:,:]'): nr, nv, nq, nz = f.shape + #$omp parallel for collapse(2) firstprivate(thetaSpline) for rIdx in range(nr): # r for cIdx in range(nv): # v # find the values of the function at each required point @@ -252,6 +261,7 @@ def poloidal_advection_step_impl(f: 'float[:,:]', dt: 'float', v: 'float', rPts: idx = nPts_r-1 rMax = rPts[idx] + #$omp parallel for collapse(2) for i in range(nPts_q): for j in range(nPts_r): # Step one of Heun method @@ -266,6 +276,7 @@ def poloidal_advection_step_impl(f: 'float[:,:]', dt: 'float', v: 'float', rPts: norm = tol+1 while (norm > tol): norm = 0.0 + #$omp parallel for collapse(2) for i in range(nPts_q): for j in range(nPts_r): # Handle theta boundary conditions @@ -315,6 +326,7 @@ def poloidal_advection_step_impl(f: 'float[:,:]', dt: 'float', v: 'float', rPts: # Find value at the determined point if (nulBound): + #$omp parallel for collapse(2) for i in range(nPts_q): for j in range(nPts_r): if (endPts_k2_r[i, j] < rPts[0]): @@ -326,6 +338,7 @@ def poloidal_advection_step_impl(f: 'float[:,:]', dt: 'float', v: 'float', rPts: f[i, j] = pol_spline.eval( endPts_k2_q[i, j], endPts_k2_r[i, j]) else: + #$omp parallel for collapse(2) for i in range(nPts_q): for j in range(nPts_r): if (endPts_k2_r[i, j] < rPts[0]): @@ -351,6 +364,7 @@ def poloidal_advection_loop(f: 'float[:,:,:,:]', phi: 'float[:,:,:]', dt: 'float CN0: 'float', kN0: 'float', deltaRN0: 'float', rp: 'float', CTi: 'float', kTi: 'float', deltaRTi: 'float', B0: 'float', tol: 'float', nulBound: 'bool'): _, nz, _, _ = f.shape + #$omp parallel for firstprivate(phi_spline, pol_spline) private(interp_wt, drPhi_0, dthetaPhi_0, drPhi_k, dthetaPhi_k, endPts_k1_q, endPts_k1_r, endPts_k2_q, endPts_k2_r) for j in range(nz): solve_2d_system(phi[j], phi_spline, interp_wt, r_bmat, r_l, r_u, r_ipiv, theta_offset, theta_splu) From 6b10610975f18b140be89f41c7b8c878d1c5a962 Mon Sep 17 00:00:00 2001 From: Emily Bourne Date: Thu, 25 Sep 2025 08:15:09 +0200 Subject: [PATCH 07/15] Keep ignoring default ignores --- .github/workflows/testing.yml | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/.github/workflows/testing.yml b/.github/workflows/testing.yml index 28b7752f..ea9fd59d 100644 --- a/.github/workflows/testing.yml +++ b/.github/workflows/testing.yml @@ -198,15 +198,15 @@ jobs: python -m pip install autopep8 - name: Check for pep changes run: | - if [[ $(autopep8 --ignore=E265 -r --diff pygyro/) ]] + if [[ $(autopep8 --ignore=E265,E226,E24,W50,W690 -r --diff pygyro/) ]] then - echo "Please run 'autopep8 --ignore=E265 -ir pygyro/' to clean up formatting" - autopep8 --ignore=E265 -r --diff pygyro/ + echo "Please run 'autopep8 --ignore=E265,E226,E24,W50,W690 -ir pygyro/' to clean up formatting" + autopep8 --ignore=E265,E226,E24,W50,W690 -r --diff pygyro/ exit 1 fi - if [[ $(autopep8 --ignore=E265 -r --diff plotting/) ]] + if [[ $(autopep8 --ignore=E265,E226,E24,W50,W690 -r --diff plotting/) ]] then - echo "Please run 'autopep8 --ignore=E265 -ir plotting/' to clean up formatting" - autopep8 --ignore=E265 -r --diff plotting/ + echo "Please run 'autopep8 --ignore=E265,E226,E24,W50,W690 -ir plotting/' to clean up formatting" + autopep8 --ignore=E265,E226,E24,W50,W690 -r --diff plotting/ exit 1 fi From ce7047eda972f6436cb0753e112e1308659507be Mon Sep 17 00:00:00 2001 From: Emily Bourne Date: Thu, 25 Sep 2025 08:21:35 +0200 Subject: [PATCH 08/15] Remove bad OpenMP --- pygyro/advection/accelerated_advection_steps.py | 1 - 1 file changed, 1 deletion(-) diff --git a/pygyro/advection/accelerated_advection_steps.py b/pygyro/advection/accelerated_advection_steps.py index c540afe0..b0184848 100644 --- a/pygyro/advection/accelerated_advection_steps.py +++ b/pygyro/advection/accelerated_advection_steps.py @@ -364,7 +364,6 @@ def poloidal_advection_loop(f: 'float[:,:,:,:]', phi: 'float[:,:,:]', dt: 'float CN0: 'float', kN0: 'float', deltaRN0: 'float', rp: 'float', CTi: 'float', kTi: 'float', deltaRTi: 'float', B0: 'float', tol: 'float', nulBound: 'bool'): _, nz, _, _ = f.shape - #$omp parallel for firstprivate(phi_spline, pol_spline) private(interp_wt, drPhi_0, dthetaPhi_0, drPhi_k, dthetaPhi_k, endPts_k1_q, endPts_k1_r, endPts_k2_q, endPts_k2_r) for j in range(nz): solve_2d_system(phi[j], phi_spline, interp_wt, r_bmat, r_l, r_u, r_ipiv, theta_offset, theta_splu) From 6e85d5a17e5001c98ca13bed060b8fd3f55882d8 Mon Sep 17 00:00:00 2001 From: Emily Bourne Date: Thu, 25 Sep 2025 08:42:22 +0200 Subject: [PATCH 09/15] OpenMP for spline interpolation --- pygyro/splines/accelerated_spline_interpolators.py | 7 +++++-- 1 file changed, 5 insertions(+), 2 deletions(-) diff --git a/pygyro/splines/accelerated_spline_interpolators.py b/pygyro/splines/accelerated_spline_interpolators.py index b60abf03..1259abc1 100644 --- a/pygyro/splines/accelerated_spline_interpolators.py +++ b/pygyro/splines/accelerated_spline_interpolators.py @@ -75,24 +75,27 @@ def solve_2d_system(ug: 'float[:,:]', spl: Spline2D, wt: 'float[:,:]', assert sinfo == 0 # Transpose coefficients to self._bwork + #$omp parallel for collapse(2) for i1 in range(n1): for i2 in range(n2): wt[i2, i1] = ug[i1, i2] # Cycle over x2 position and interpolate w along x1 direction. # Work on self._bwork + #$ omp parallel for firstprivate(spline1) for i2 in range(n2): solve_system_periodic(wt[i2, :n1], spline1, theta_offset, theta_splu) # self._interp1.compute_interpolant(wt[i2, :n1], self._spline1) - c = spline1.coeffs - wt[i2, :] = c + wt[i2, :] = spline1.coeffs # Transpose coefficients to spl.coeffs + #$omp parallel for collapse(2) for i1 in range(s1): for i2 in range(s2): w[i1, i2] = wt[i2, i1] # x1-periodic only: "wrap around" coefficients onto extended array + #$omp parallel for collapse(2) for i1 in range(p1): for i2 in range(s2): w[n1 + i1, i2] = w[i1, i2] From 4ebe36fc9a440aabab38ccf4735b69bab27d5e2a Mon Sep 17 00:00:00 2001 From: Emily Bourne Date: Thu, 25 Sep 2025 08:51:03 +0200 Subject: [PATCH 10/15] Clean up --- pygyro/splines/accelerated_spline_interpolators.py | 6 ++---- 1 file changed, 2 insertions(+), 4 deletions(-) diff --git a/pygyro/splines/accelerated_spline_interpolators.py b/pygyro/splines/accelerated_spline_interpolators.py index 1259abc1..d258aa63 100644 --- a/pygyro/splines/accelerated_spline_interpolators.py +++ b/pygyro/splines/accelerated_spline_interpolators.py @@ -6,7 +6,6 @@ T = TypeVar('T', float, complex) - def solve_system_periodic(ug: 'Final[float[:]]', spl: Spline1D, offset: int, splu: Final[PeriodicBandedMatrix]): """ Compute the coefficients c of the spline which interpolates the points ug @@ -53,7 +52,7 @@ def solve_system_nonperiodic(ug: 'Final[T[:]]', c: 'T[:]', bmat: 'Final[T[:,:](o def solve_2d_system(ug: 'float[:,:]', spl: Spline2D, wt: 'float[:,:]', r_bmat: 'float[:,:](order=F)', r_l: np.int32, r_u: np.int32, r_ipiv: 'int32[:]', - theta_offset: int, theta_splu: PeriodicBandedMatrix): + theta_offset: int, theta_splu: Final[PeriodicBandedMatrix]): basis1 = spl.basis1 basis2 = spl.basis2 n1, n2 = basis1.nbasis, basis2.nbasis @@ -82,10 +81,9 @@ def solve_2d_system(ug: 'float[:,:]', spl: Spline2D, wt: 'float[:,:]', # Cycle over x2 position and interpolate w along x1 direction. # Work on self._bwork - #$ omp parallel for firstprivate(spline1) + ##$ omp parallel for firstprivate(spline1) for i2 in range(n2): solve_system_periodic(wt[i2, :n1], spline1, theta_offset, theta_splu) - # self._interp1.compute_interpolant(wt[i2, :n1], self._spline1) wt[i2, :] = spline1.coeffs # Transpose coefficients to spl.coeffs From 04a6ecdaa57e0ffa065a9dc8cb673fc188222324 Mon Sep 17 00:00:00 2001 From: Emily Bourne Date: Thu, 25 Sep 2025 12:14:51 +0200 Subject: [PATCH 11/15] Add comment about why no parallel --- pygyro/splines/accelerated_spline_interpolators.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/pygyro/splines/accelerated_spline_interpolators.py b/pygyro/splines/accelerated_spline_interpolators.py index d258aa63..dc71cc52 100644 --- a/pygyro/splines/accelerated_spline_interpolators.py +++ b/pygyro/splines/accelerated_spline_interpolators.py @@ -81,7 +81,8 @@ def solve_2d_system(ug: 'float[:,:]', spl: Spline2D, wt: 'float[:,:]', # Cycle over x2 position and interpolate w along x1 direction. # Work on self._bwork - ##$ omp parallel for firstprivate(spline1) + ##$ omp parallel for default(none) shared(wt, n1, n2, theta_offset, theta_splu, basis1) private(spline1) + # theta_splu contains schur which can be modified by thread for i2 in range(n2): solve_system_periodic(wt[i2, :n1], spline1, theta_offset, theta_splu) wt[i2, :] = spline1.coeffs From 2f70a562499222008b52d261590d0b251d257cc8 Mon Sep 17 00:00:00 2001 From: Emily Bourne Date: Thu, 25 Sep 2025 12:30:20 +0200 Subject: [PATCH 12/15] Correct order --- pygyro/advection/accelerated_advection_steps.py | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/pygyro/advection/accelerated_advection_steps.py b/pygyro/advection/accelerated_advection_steps.py index b0184848..c82aae4e 100644 --- a/pygyro/advection/accelerated_advection_steps.py +++ b/pygyro/advection/accelerated_advection_steps.py @@ -45,11 +45,11 @@ def poloidal_advection_step_expl(f: 'float[:,:]', nPts_q = qPts.shape[0] #$omp parallel for collapse(2) - for i in range(nPts_r): - for j, q in enumerate(qPts): - r = rPts[i] - drPhi_0[i, j] = phi_spline.eval(q, r, 0, 1) - dthetaPhi_0[i, j] = phi_spline.eval(q, r, 1, 0) + for j in range(nPts_q): + for i, r in enumerate(rPts): + q = qPts[j] + drPhi_0[j, i] = phi_spline.eval(q, r, 0, 1) + dthetaPhi_0[j, i] = phi_spline.eval(q, r, 1, 0) idx = nPts_r-1 rMax = rPts[idx] From 7907b9cf4d3bf860060467ab808100b58d23767d Mon Sep 17 00:00:00 2001 From: Emily Bourne Date: Thu, 25 Sep 2025 12:34:51 +0200 Subject: [PATCH 13/15] Fix OpenMP call --- pygyro/splines/accelerated_spline_interpolators.py | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/pygyro/splines/accelerated_spline_interpolators.py b/pygyro/splines/accelerated_spline_interpolators.py index dc71cc52..1fa23983 100644 --- a/pygyro/splines/accelerated_spline_interpolators.py +++ b/pygyro/splines/accelerated_spline_interpolators.py @@ -81,8 +81,7 @@ def solve_2d_system(ug: 'float[:,:]', spl: Spline2D, wt: 'float[:,:]', # Cycle over x2 position and interpolate w along x1 direction. # Work on self._bwork - ##$ omp parallel for default(none) shared(wt, n1, n2, theta_offset, theta_splu, basis1) private(spline1) - # theta_splu contains schur which can be modified by thread + #$ omp parallel for default(none) shared(wt, n1, n2, theta_offset) firstprivate(spline1, theta_splu) for i2 in range(n2): solve_system_periodic(wt[i2, :n1], spline1, theta_offset, theta_splu) wt[i2, :] = spline1.coeffs From a558c6c064ddd6e26b224af403db60f939e4b0d8 Mon Sep 17 00:00:00 2001 From: Emily Bourne Date: Thu, 25 Sep 2025 12:36:28 +0200 Subject: [PATCH 14/15] Pep --- pygyro/splines/accelerated_spline_interpolators.py | 1 + 1 file changed, 1 insertion(+) diff --git a/pygyro/splines/accelerated_spline_interpolators.py b/pygyro/splines/accelerated_spline_interpolators.py index 1fa23983..807cc510 100644 --- a/pygyro/splines/accelerated_spline_interpolators.py +++ b/pygyro/splines/accelerated_spline_interpolators.py @@ -6,6 +6,7 @@ T = TypeVar('T', float, complex) + def solve_system_periodic(ug: 'Final[float[:]]', spl: Spline1D, offset: int, splu: Final[PeriodicBandedMatrix]): """ Compute the coefficients c of the spline which interpolates the points ug From 2024cdf44e491237facff2d9e1167c8c54c90142 Mon Sep 17 00:00:00 2001 From: Emily Bourne Date: Thu, 25 Sep 2025 14:57:34 +0200 Subject: [PATCH 15/15] Use static schedule --- pygyro/splines/accelerated_spline_interpolators.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/pygyro/splines/accelerated_spline_interpolators.py b/pygyro/splines/accelerated_spline_interpolators.py index 807cc510..467cecd6 100644 --- a/pygyro/splines/accelerated_spline_interpolators.py +++ b/pygyro/splines/accelerated_spline_interpolators.py @@ -82,7 +82,7 @@ def solve_2d_system(ug: 'float[:,:]', spl: Spline2D, wt: 'float[:,:]', # Cycle over x2 position and interpolate w along x1 direction. # Work on self._bwork - #$ omp parallel for default(none) shared(wt, n1, n2, theta_offset) firstprivate(spline1, theta_splu) + #$ omp parallel for default(none) shared(wt, n1, n2, theta_offset) firstprivate(spline1, theta_splu) schedule(static) for i2 in range(n2): solve_system_periodic(wt[i2, :n1], spline1, theta_offset, theta_splu) wt[i2, :] = spline1.coeffs