Skip to content

Replace numpy with cunumpy #44

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Draft
max-models wants to merge 34 commits into
base: devel-tiny
Choose a base branch
from
Draft

Replace numpy with cunumpy #44

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
34 commits
Select commit Hold shift + click to select a range
196cc1d
Added psydac/arrays.py
max-models Oct 10, 2025
e6d1fba
Replaced import numpy as np with from psydac.arrays import xp as np
max-models Oct 10, 2025
32d5797
Bugfix
max-models Oct 10, 2025
0a78b05
Fixes for running psydac with cupy arrays
max-models Oct 10, 2025
721ff30
Merge branch 'devel-tiny' into add-xp-backend
max-models Oct 20, 2025
367a7be
Remplaced numpy imports
max-models Oct 20, 2025
6bb76a7
Set arguments to int in view.shape = (int(self._numrhs), int(self._di…
max-models Oct 20, 2025
95faef0
Convert gridsizes tuple to array
max-models Oct 20, 2025
aa4ff3e
Always call _splu.solve with numpy array
max-models Oct 20, 2025
c262eb3
Convert to tuple of integers
max-models Oct 20, 2025
9831134
Call self._solver_function with numpy arrays
max-models Oct 20, 2025
6898e83
Convert self._perm components to int
max-models Oct 20, 2025
39dd8e8
Merge branch 'devel-tiny' into add-xp-backend
max-models Oct 22, 2025
6c5f003
Replaces arrays.py with cunumpy
max-models Oct 22, 2025
46d30d3
Replaced np. with xp.
max-models Oct 22, 2025
86498e3
Use if array_backend.backend == 'cupy' to check if cupy or numpy is used
max-models Oct 28, 2025
e1a1005
Use if array_backend.backend == 'cupy' to check if cupy or numpy is used
max-models Oct 28, 2025
3bfc615
Merge branch 'devel-tiny' into add-xp-backend
max-models Oct 28, 2025
c322ee3
Merge remote-tracking branch 'origin/devel-tiny' into add-xp-backend
Copilot May 5, 2026
316b080
Update petsc4py install
max-models Jul 23, 2026
e7cecca
Added oversubscribe flag
max-models Jul 23, 2026
420d008
Only run PR - test Struphy in container with fortran
max-models Jul 23, 2026
a24495d
Merge branch '69-error-failed-building-wheel-for-petsc4py' into add-x…
max-models Jul 23, 2026
cc84bb7
np --> xp
max-models Jul 23, 2026
eb92044
Use np in kernels
max-models Jul 23, 2026
e79bde0
Merge branch 'devel-tiny' into add-xp-backend
max-models Jul 24, 2026
4675106
int --> int64
max-models Jul 24, 2026
8e45218
Update nonblocking_data_exchanger.py, cp --> np
max-models Jul 29, 2026
4d33ad3
Update blocking_data_exchanger.py, cp --> np
max-models Jul 29, 2026
c926ed7
Update blocking_data_exchanger.py
max-models Jul 29, 2026
61f90a9
Update nonblocking_data_exchanger.py
max-models Jul 29, 2026
b142a59
Use sim(M) since xp.sum(M) doesn't work for python lists
max-models Jul 30, 2026
9e278b7
Use xp.max for xp array in splines.py
max-models Jul 30, 2026
b8ee897
Handle default case explicitly
max-models Jul 30, 2026
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
60 changes: 30 additions & 30 deletions feectools/api/fem_bilinear_form.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -2,7 +2,7 @@
import os
import importlib

import numpy as np
import cunumpy as xp

from sympy import ImmutableDenseMatrix, Matrix, Symbol, sympify
from sympy.tensor.indexed import Indexed, IndexedBase
Expand Down Expand Up @@ -465,8 +465,8 @@ def allocate_matrices(self, backend=None):

expr = self.kernel_expr.expr
target = self.kernel_expr.target
test_degree = np.array(self.test_basis.space.degree)
trial_degree = np.array(self.trial_basis.space.degree)
test_degree = xp.array(self.test_basis.space.degree)
trial_degree = xp.array(self.trial_basis.space.degree)
test_space = self.spaces[1].coeff_space
trial_space = self.spaces[0].coeff_space
test_fem_space = self.spaces[1]
Expand All @@ -476,20 +476,20 @@ def allocate_matrices(self, backend=None):
is_conformal = True

if isinstance(expr, (ImmutableDenseMatrix, Matrix)):
if not isinstance(test_degree[0],(list, tuple, np.ndarray)):
if not isinstance(test_degree[0],(list, tuple, xp.ndarray)):
test_degree = [test_degree]

if not isinstance(trial_degree[0],(list, tuple, np.ndarray)):
if not isinstance(trial_degree[0],(list, tuple, xp.ndarray)):
trial_degree = [trial_degree]

pads = np.empty((len(test_degree),len(trial_degree),len(test_degree[0])), dtype=int)
pads = xp.empty((len(test_degree),len(trial_degree),len(test_degree[0])), dtype=int)
for i in range(len(test_degree)):
for j in range(len(trial_degree)):
td = test_degree[i]
trd = trial_degree[j]
pads[i,j][:] = np.array([td, trd]).max(axis=0)
pads[i,j][:] = xp.array([td, trd]).max(axis=0)
else:
pads = np.maximum(test_degree, trial_degree)
pads = xp.maximum(test_degree, trial_degree)

if self._matrix is None and (is_broken or isinstance(expr, (ImmutableDenseMatrix, Matrix))):
self._matrix = BlockLinearOperator(trial_space, test_space)
Expand Down Expand Up @@ -714,8 +714,8 @@ def assemble(self, *, reset=True, **kwargs):
bs, d, s, p, mult = construct_test_space_arguments(basis_v)
basis += bs
spans += s
degrees += [np.int64(a) for a in d]
pads += [np.int64(a) for a in p]
degrees += [xp.int64(a) for a in d]
pads += [xp.int64(a) for a in p]
if v.space.is_multipatch or v.space.is_vector_valued:
coeffs += (e._data for e in v.coeffs)
else:
Expand Down Expand Up @@ -1382,12 +1382,12 @@ def make_file(self, temps, ordered_stmts, field_derivatives, max_logical_derivat
I_1 = f'int(floor(i_1/{test_mult[0]})*{trial_mult[0]})' if max(test_mult[0], trial_mult[0]) > 1 else 'i_1'
I_2 = f'int(floor(i_2/{test_mult[1]})*{trial_mult[1]})' if max(test_mult[1], trial_mult[1]) > 1 else 'i_2'
I_3 = f'int(floor(i_3/{test_mult[2]})*{trial_mult[2]})' if max(test_mult[2], trial_mult[2]) > 1 else 'i_3'
#MAX_P1 = max(int( ( MAX_P1 + np.floor(MAX_P1 / test_mult[0]) * trial_mult[0] ) / 2 ), MAX_P1) if max(test_mult[0], trial_mult[0]) > 1 else MAX_P1
#MAX_P2 = max(int( ( MAX_P2 + np.floor(MAX_P2 / test_mult[1]) * trial_mult[1] ) / 2 ), MAX_P2) if max(test_mult[1], trial_mult[1]) > 1 else MAX_P2
#MAX_P3 = max(int( ( MAX_P3 + np.floor(MAX_P3 / test_mult[2]) * trial_mult[2] ) / 2 ), MAX_P3) if max(test_mult[2], trial_mult[2]) > 1 else MAX_P3
n_cols_x1 = max( int(MAX_P1 + 1 + np.floor(MAX_P1 / test_mult[0]) * trial_mult[0]), 2*MAX_P1+1 )
n_cols_x2 = max( int(MAX_P2 + 1 + np.floor(MAX_P2 / test_mult[1]) * trial_mult[1]), 2*MAX_P2+1 )
n_cols_x3 = max( int(MAX_P3 + 1 + np.floor(MAX_P3 / test_mult[2]) * trial_mult[2]), 2*MAX_P3+1 )
#MAX_P1 = max(int( ( MAX_P1 + xp.floor(MAX_P1 / test_mult[0]) * trial_mult[0] ) / 2 ), MAX_P1) if max(test_mult[0], trial_mult[0]) > 1 else MAX_P1
#MAX_P2 = max(int( ( MAX_P2 + xp.floor(MAX_P2 / test_mult[1]) * trial_mult[1] ) / 2 ), MAX_P2) if max(test_mult[1], trial_mult[1]) > 1 else MAX_P2
#MAX_P3 = max(int( ( MAX_P3 + xp.floor(MAX_P3 / test_mult[2]) * trial_mult[2] ) / 2 ), MAX_P3) if max(test_mult[2], trial_mult[2]) > 1 else MAX_P3
n_cols_x1 = max( int(MAX_P1 + 1 + xp.floor(MAX_P1 / test_mult[0]) * trial_mult[0]), 2*MAX_P1+1 )
n_cols_x2 = max( int(MAX_P2 + 1 + xp.floor(MAX_P2 / test_mult[1]) * trial_mult[1]), 2*MAX_P2+1 )
n_cols_x3 = max( int(MAX_P3 + 1 + xp.floor(MAX_P3 / test_mult[2]) * trial_mult[2]), 2*MAX_P3+1 )
MAX_P1 = n_cols_x1 - MAX_P1 - 1
MAX_P2 = n_cols_x2 - MAX_P2 - 1
MAX_P3 = n_cols_x3 - MAX_P3 - 1
Expand Down Expand Up @@ -2024,9 +2024,9 @@ def construct_arguments_generate_assembly_file(self):

# keys_2[(u[0], v[1])][3] = (1,2) means that the fourth sub-expression corresponding to the trial-test-function-component-product
# u[0] * v[1] involves a first derivative in x2 direction of the trial function and a second derivative in x2 direction of the test function
keys_1[block] = np.array([(alpha_1, beta_1) for alpha_1, beta_1 in zip(x1_trial_keys[block], x1_test_keys[block])])
keys_2[block] = np.array([(alpha_2, beta_2) for alpha_2, beta_2 in zip(x2_trial_keys[block], x2_test_keys[block])])
keys_3[block] = np.array([(alpha_3, beta_3) for alpha_3, beta_3 in zip(x3_trial_keys[block], x3_test_keys[block])])
keys_1[block] = xp.array([(alpha_1, beta_1) for alpha_1, beta_1 in zip(x1_trial_keys[block], x1_test_keys[block])])
keys_2[block] = xp.array([(alpha_2, beta_2) for alpha_2, beta_2 in zip(x2_trial_keys[block], x2_test_keys[block])])
keys_3[block] = xp.array([(alpha_3, beta_3) for alpha_3, beta_3 in zip(x3_trial_keys[block], x3_test_keys[block])])

# Those are the function values in each direction of a particular component of the trial/test function
global_basis_u_1, global_basis_u_2, global_basis_u_3 = global_basis_u[u_i]
Expand Down Expand Up @@ -2056,9 +2056,9 @@ def construct_arguments_generate_assembly_file(self):
# of non-zero product.
# Hence, we assign zeros for each element, each quadrature point on the element, each test and trial function combination,
# and each (or even more than required) appearing partial derivative combination of these functions - in each direction
test_trial_1 = np.zeros((n_element_1, k1, test_v_p1 + 1, trial_u_p1 + 1, max_block_trial_x1_derivative+1, max_block_test_x1_derivative+1), dtype='float64')
test_trial_2 = np.zeros((n_element_2, k2, test_v_p2 + 1, trial_u_p2 + 1, max_block_trial_x2_derivative+1, max_block_test_x2_derivative+1), dtype='float64')
test_trial_3 = np.zeros((n_element_3, k3, test_v_p3 + 1, trial_u_p3 + 1, max_block_trial_x3_derivative+1, max_block_test_x3_derivative+1), dtype='float64')
test_trial_1 = xp.zeros((n_element_1, k1, test_v_p1 + 1, trial_u_p1 + 1, max_block_trial_x1_derivative+1, max_block_test_x1_derivative+1), dtype='float64')
test_trial_2 = xp.zeros((n_element_2, k2, test_v_p2 + 1, trial_u_p2 + 1, max_block_trial_x2_derivative+1, max_block_test_x2_derivative+1), dtype='float64')
test_trial_3 = xp.zeros((n_element_3, k3, test_v_p3 + 1, trial_u_p3 + 1, max_block_trial_x3_derivative+1, max_block_test_x3_derivative+1), dtype='float64')

# And that's how we fill the test_trial arrays
if self._pyccelize_test_trial_computation and assembly_backend['name'] == 'pyccel':
Expand All @@ -2068,7 +2068,7 @@ def construct_arguments_generate_assembly_file(self):
[max_block_trial_x1_derivative, max_block_trial_x2_derivative, max_block_trial_x3_derivative], [max_block_test_x1_derivative, max_block_test_x2_derivative, max_block_test_x3_derivative],
[test_trial_1, test_trial_2, test_trial_3]):

args = tuple(np.int64(a) if isinstance(a, int) else a for a in args)
args = tuple(xp.int64(a) if isinstance(a, int) else a for a in args)

test_trial_func(*args)
else:
Expand Down Expand Up @@ -2118,13 +2118,13 @@ def construct_arguments_generate_assembly_file(self):
# to store local information correctly. 2*degree+1 in the simplest case.
n_funs_x2 = n_element_2 + test_v_p2 + (test_mult[1]-1)*(n_element_2-1)
n_funs_x3 = n_element_3 + test_v_p3 + (test_mult[2]-1)*(n_element_3-1)
n_cols_x2 = max( int(max_p_2 + 1 + np.floor(max_p_2 / test_mult[1]) * trial_mult[1]), 2*max_p_2+1 )
n_cols_x3 = max( int(max_p_3 + 1 + np.floor(max_p_3 / test_mult[2]) * trial_mult[2]), 2*max_p_3+1 )
n_cols_x2 = max( int(max_p_2 + 1 + xp.floor(max_p_2 / test_mult[1]) * trial_mult[1]), 2*max_p_2+1 )
n_cols_x3 = max( int(max_p_3 + 1 + xp.floor(max_p_3 / test_mult[2]) * trial_mult[2]), 2*max_p_3+1 )

a3[block] = np.zeros((n_expr, n_funs_x3, n_cols_x3), dtype='float64')
a2[block] = np.zeros((n_expr, n_funs_x2, n_funs_x3, n_cols_x2, n_cols_x3), dtype='float64')
a3[block] = xp.zeros((n_expr, n_funs_x3, n_cols_x3), dtype='float64')
a2[block] = xp.zeros((n_expr, n_funs_x2, n_funs_x3, n_cols_x2, n_cols_x3), dtype='float64')

coupling_terms[block] = np.zeros((n_element_2, k2, n_element_3, k3, n_expr), dtype='float64')
coupling_terms[block] = xp.zeros((n_element_2, k2, n_element_3, k3, n_expr), dtype='float64')

# We gather the socalled new args - all other args are being obtained in a similar way using the old assembly implementation
new_args = (*list(test_trial_1s.values()),
Expand Down Expand Up @@ -2158,8 +2158,8 @@ def construct_arguments_generate_assembly_file(self):

threads_args = ()

args = tuple(np.int64(a) if isinstance(a, int) else a for a in args)
threads_args = tuple(np.int64(a) if isinstance(a, int) else a for a in threads_args)
args = tuple(xp.int64(a) if isinstance(a, int) else a for a in args)
threads_args = tuple(xp.int64(a) if isinstance(a, int) else a for a in threads_args)

#---------- We now generate the assembly file ----------

Expand Down
6 changes: 3 additions & 3 deletions feectools/api/fem_common.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -2,7 +2,7 @@

from sympy import Expr, ImmutableDenseMatrix, Matrix

import numpy as np
import cunumpy as xp

from sympde.expr.basic import BasicForm
from sympde.expr.evaluation import KernelExpression
Expand Down Expand Up @@ -215,8 +215,8 @@ def collect_spaces(space, *args):

#==============================================================================
def compute_diag_len(p, md, mc):
n = ((np.ceil((p+1)/mc)-1)*md).astype('int')
n = n-np.minimum(0, n-p)+p+1
n = ((xp.ceil((p+1)/mc)-1)*md).astype('int')
n = n-xp.minimum(0, n-p)+p+1
return n.astype('int')

#==============================================================================
Expand Down
4 changes: 2 additions & 2 deletions feectools/api/fem_sum_form.py
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1,4 +1,4 @@
import numpy as np
import cunumpy as xp

from sympde.expr.expr import (
BilinearForm as sym_BilinearForm,
Expand Down Expand Up @@ -119,5 +119,5 @@ def assemble(self, *, reset=True, **kwargs):
return self._operator
else:
M = [form.assemble(**kwargs) for form in self.forms]
M = np.sum(M)
M = sum(M)
return M
Comment thread
max-models marked this conversation as resolved.
Loading
Loading