Added all changes from struphy-with-cupy

params_LinearMHDDriftkineticCC.py

@@ -0,0 +1,116 @@
+# import model, set verbosity
+from struphy.models.hybrid import LinearMHDDriftkineticCC
+
+from struphy import main
+from struphy.fields_background import equils
+from struphy.geometry import domains
+from struphy.initial import perturbations
+from struphy.io.options import BaseUnits, DerhamOptions, EnvironmentOptions, FieldsBackground, Time
+from struphy.kinetic_background import maxwellians
+from struphy.pic.utilities import (
+    BinningPlot,
+    BoundaryParameters,
+    KernelDensityPlot,
+    LoadingParameters,
+    WeightsParameters,
+)
+from struphy.topology import grids
+
+# environment options
+env = EnvironmentOptions()
+
+# units
+base_units = BaseUnits()
+
+# time stepping
+time_opts = Time()
+
+# geometry
+domain = domains.Cuboid()
+
+# fluid equilibrium (can be used as part of initial conditions)
+equil = equils.HomogenSlab()
+
+# grid
+grid = grids.TensorProductGrid(Nel=(16, 16, 16))
+
+# derham options
+derham_opts = DerhamOptions()
+
+# light-weight model instance
+model = LinearMHDDriftkineticCC()
+
+# species parameters
+model.mhd.set_phys_params()
+model.energetic_ions.set_phys_params()
+
+loading_params = LoadingParameters(ppc=1000)
+weights_params = WeightsParameters()
+boundary_params = BoundaryParameters()
+model.energetic_ions.set_markers(
+    loading_params=loading_params,
+    weights_params=weights_params,
+    boundary_params=boundary_params,
+)
+model.energetic_ions.set_sorting_boxes()
+model.energetic_ions.set_save_data()
+
+# propagator options
+model.propagators.push_bxe.options = model.propagators.push_bxe.Options(
+    b_tilde=model.em_fields.b_field,
+)
+model.propagators.push_parallel.options = model.propagators.push_parallel.Options(
+    b_tilde=model.em_fields.b_field,
+)
+model.propagators.shearalfen_cc5d.options = model.propagators.shearalfen_cc5d.Options(
+    energetic_ions=model.energetic_ions.var,
+)
+model.propagators.magnetosonic.options = model.propagators.magnetosonic.Options(
+    b_field=model.em_fields.b_field,
+)
+model.propagators.cc5d_density.options = model.propagators.cc5d_density.Options(
+    energetic_ions=model.energetic_ions.var,
+    b_tilde=model.em_fields.b_field,
+)
+model.propagators.cc5d_gradb.options = model.propagators.cc5d_gradb.Options(
+    b_tilde=model.em_fields.b_field,
+)
+model.propagators.cc5d_curlb.options = model.propagators.cc5d_curlb.Options(
+    b_tilde=model.em_fields.b_field,
+)
+
+# background, perturbations and initial conditions
+model.mhd.velocity.add_background(FieldsBackground())
+model.mhd.velocity.add_perturbation(perturbations.TorusModesCos(given_in_basis="v", comp=0))
+model.mhd.velocity.add_perturbation(perturbations.TorusModesCos(given_in_basis="v", comp=1))
+model.mhd.velocity.add_perturbation(perturbations.TorusModesCos(given_in_basis="v", comp=2))
+maxwellian_1 = maxwellians.GyroMaxwellian2D(n=(1.0, None), equil=equil)
+maxwellian_2 = maxwellians.GyroMaxwellian2D(n=(0.1, None), equil=equil)
+background = maxwellian_1 + maxwellian_2
+model.energetic_ions.var.add_background(background)
+
+# if .add_initial_condition is not called, the background is the kinetic initial condition
+perturbation = perturbations.TorusModesCos()
+maxwellian_1pt = maxwellians.GyroMaxwellian2D(n=(1.0, perturbation), equil=equil)
+init = maxwellian_1pt + maxwellian_2
+model.energetic_ions.var.add_initial_condition(init)
+
+# optional: exclude variables from saving
+# model.energetic_ions.var.save_data = False
+
+if __name__ == "__main__":
+    # start run
+    verbose = True
+
+    main.run(
+        model,
+        params_path=__file__,
+        env=env,
+        base_units=base_units,
+        time_opts=time_opts,
+        domain=domain,
+        equil=equil,
+        grid=grid,
+        derham_opts=derham_opts,
+        verbose=verbose,
+    )

pyproject.toml

@@ -24,9 +24,9 @@ classifiers = [
 ]
 dependencies = [
     "numpy<=2.5.0",
-    "cunumpy<=0.1.1",
+    "cunumpy>=0.1.4, <=0.1.5",
     "pyccel>=2.2.0, <=2.2.3",
-    "feectools<=0.1.3",
+    "feectools<=0.1.8",
     "scipy<=1.18.0",
     "h5py<=3.16.0",
     "matplotlib<=3.11.0",

src/struphy/bsplines/bsplines.py

@@ -611,8 +611,8 @@ def make_knots(breaks, degree, periodic):
 
     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] = xp.asarray([xi - period for xi in breaks[-p - 1 : -1]])
+        T[-p:] = xp.asarray([xi + period for xi in breaks[1 : p + 1]])
     else:
         T[0:p] = breaks[0]
         T[-p:] = breaks[-1]

src/struphy/console/format.py

@@ -492,9 +492,7 @@ def get_python_files(input_type, path=None):
 
         # python_files = [f for f in files if f.endswith(".py") and os.path.isfile(f)]
         python_files = [
-            os.path.join(repopath, f)
-            for f in files
-            if (f.endswith(".py") or f.endswith(".ipynb")) and os.path.isfile(os.path.join(repopath, f))
+            os.path.join(repopath, f) for f in files if f.endswith(".py") and os.path.isfile(os.path.join(repopath, f))
         ]
 
         if not python_files:

src/struphy/console/test.py

@@ -73,6 +73,7 @@ def struphy_test(
             cmd = [
                 "pytest",
                 "--testmon",
+                "-xvs",
             ] + list_of_tests
 
         if with_desc:
@@ -156,7 +157,6 @@ def struphy_test(
             "pytest",
             "-m",
             "single",
-            "-xvs",
             # "--with-mpi",
             "--model-name",
             group,

src/struphy/feec/basis_projection_ops.py

@@ -1,6 +1,8 @@
 import logging
 
 import cunumpy as xp
+import numpy as np
+from cunumpy import PyccelKernel
 from feectools.api.settings import PSYDAC_BACKEND_GPYCCEL
 from feectools.ddm.mpi import mpi as MPI
 from feectools.fem.basic import FemSpace
@@ -19,7 +21,6 @@
 from struphy.polar.basic import PolarDerhamSpace, PolarVector
 from struphy.polar.linear_operators import PolarExtractionOperator
 from struphy.utils.docstring_converter import auto_convert_docstring
-from struphy.utils.pyccel import Pyccelkernel
 
 logger = logging.getLogger("struphy")
 
@@ -53,9 +54,7 @@ def __init__(self, derham, domain, **weights):
 
         self._rank = derham.comm.Get_rank() if derham.comm is not None else 0
 
-        if xp.any(
-            [degree == 1 and num_elements > 1 for degree, num_elements in zip(derham.degree, derham.num_elements)]
-        ):
+        if any([degree == 1 and num_elements > 1 for degree, num_elements in zip(derham.degree, derham.num_elements)]):
             logger.warning(
                 f'WARNING: Class "BasisProjectionOperators" called with degree={derham.degree} (interpolation of piece-wise constants should be avoided).',
             )
@@ -1057,11 +1056,11 @@ def __init__(
         if isinstance(V, TensorFemSpace):
             self._Vspaces = [V.coeff_space]
             self._V1ds = [V.spaces]
-            self._VNbasis = xp.array([self._V1ds[0][0].nbasis, self._V1ds[0][1].nbasis, self._V1ds[0][2].nbasis])
+            self._VNbasis = np.array([self._V1ds[0][0].nbasis, self._V1ds[0][1].nbasis, self._V1ds[0][2].nbasis])
         else:
             self._Vspaces = V.coeff_space
             self._V1ds = [comp.spaces for comp in V.spaces]
-            self._VNbasis = xp.array(
+            self._VNbasis = np.array(
                 [
                     [self._V1ds[0][0].nbasis, self._V1ds[0][1].nbasis, self._V1ds[0][2].nbasis],
                     [
@@ -1949,13 +1948,13 @@ def assemble(self, weights=None):
 
             # input vector space (domain), column of block
             for j, (Vspace, V1d, loc_weight) in enumerate(zip(_Vspaces, _V1ds, weight_line)):
-                _starts_in = xp.array(Vspace.starts)
-                _ends_in = xp.array(Vspace.ends)
-                _pads_in = xp.array(Vspace.pads)
+                _starts_in = np.array(Vspace.starts)
+                _ends_in = np.array(Vspace.ends)
+                _pads_in = np.array(Vspace.pads)
 
-                _starts_out = xp.array(Wspace.starts)
-                _ends_out = xp.array(Wspace.ends)
-                _pads_out = xp.array(Wspace.pads)
+                _starts_out = np.array(Wspace.starts)
+                _ends_out = np.array(Wspace.ends)
+                _pads_out = np.array(Wspace.pads)
 
                 # use cached information if asked
                 if self._use_cache:
@@ -2043,7 +2042,7 @@ def assemble(self, weights=None):
                         )
                         dofs_mat = self._dof_mat[i, j]
 
-                    kernel = Pyccelkernel(
+                    kernel = PyccelKernel(
                         getattr(
                             basis_projection_kernels,
                             "assemble_dofs_for_weighted_basisfuns_" + str(V.ldim) + "d",

src/struphy/feec/boundary_mass.py

@@ -2,6 +2,7 @@
 from typing import Callable
 
 import cunumpy as xp
+from cunumpy import PyccelKernel
 from feectools.api.settings import PSYDAC_BACKEND_GPYCCEL
 from feectools.linalg.block import BlockLinearOperator, BlockVector
 from feectools.linalg.stencil import StencilMatrix, StencilVector
@@ -11,7 +12,6 @@
 from struphy.feec.mass import WeightedMassOperators
 from struphy.feec.psydac_derham import Derham, SplineFunction
 from struphy.geometry.base import Domain
-from struphy.utils.pyccel import Pyccelkernel
 
 logger = logging.getLogger("struphy")
 
@@ -174,7 +174,7 @@ def __init__(
         self._setup_surface_data()
 
         # load assembly kernel
-        self._assembly_kernel = Pyccelkernel(mass_kernels.surface_kernel_3d_mat)
+        self._assembly_kernel = PyccelKernel(mass_kernels.surface_kernel_3d_mat)
 
         self.assemble()
 

src/struphy/feec/mass.py

@@ -4,6 +4,7 @@
 from typing import Callable
 
 import cunumpy as xp
+from cunumpy import PyccelKernel
 from feectools.api.settings import PSYDAC_BACKEND_GPYCCEL
 from feectools.ddm.mpi import MockComm
 from feectools.ddm.mpi import mpi as MPI
@@ -26,7 +27,6 @@
 from struphy.polar.basic import PolarVector
 from struphy.polar.linear_operators import PolarExtractionOperator
 from struphy.utils.docstring_converter import auto_convert_docstring, info
-from struphy.utils.pyccel import Pyccelkernel
 from struphy.utils.utils import __class_with_params_repr_no_defaults__
 
 logger = logging.getLogger("struphy")
@@ -1856,7 +1856,7 @@ def __init__(
 
         # load assembly kernel
         if not self._matrix_free:
-            self._assembly_kernel = Pyccelkernel(
+            self._assembly_kernel = PyccelKernel(
                 getattr(
                     mass_kernels,
                     "kernel_" + str(self._V.ldim) + "d_mat",
@@ -2457,7 +2457,7 @@ def eval_quad(self, W, coeffs, out=None):
                 assert isinstance(out, (list, tuple))
 
         # load assembly kernel
-        kernel = Pyccelkernel(getattr(mass_kernels, "kernel_" + str(W.ldim) + "d_eval"))
+        kernel = PyccelKernel(getattr(mass_kernels, "kernel_" + str(W.ldim) + "d_eval"))
 
         # loop over components
         for a, wspace in enumerate(Wspaces):
@@ -2555,14 +2555,14 @@ def __init__(self, derham, V, W, weights=None, nquads=None):
         self._nquads = nquads
 
         self._dtype = V.coeff_space.dtype
-        self._dot_kernel = Pyccelkernel(
+        self._dot_kernel = PyccelKernel(
             getattr(
                 mass_kernels,
                 "kernel_" + str(self._V.ldim) + "d_matrixfree",
             ),
         )
 
-        self._diag_kernel = Pyccelkernel(
+        self._diag_kernel = PyccelKernel(
             getattr(
                 mass_kernels,
                 "kernel_" + str(self._V.ldim) + "d_diag",

src/struphy/feec/preconditioner.py

@@ -371,16 +371,16 @@ def solver(self):
         """KroneckerLinearSolver or BlockDiagonalSolver for exactly inverting the approximate mass matrix self.matrix."""
         return self._solver
 
+    @property
+    def domain(self):
+        """The domain of the linear operator - an element of Vectorspace"""
+        return self._space
+
     @property
     def codomain(self):
         """The codomain of the linear operator - an element of Vectorspace"""
         return self._codomain
 
-    @property
-    def domain(self):
-        """The domain of the linear operator - an element of Vectorspace"""
-        return self._domain
-
     @property
     def dtype(self):
         return self._dtype
@@ -750,6 +750,9 @@ def matrix(self):
     def solver(self):
         """KroneckerLinearSolver or BlockDiagonalSolver for exactly inverting the approximate mass matrix self.matrix."""
         return self._solver
+
+    @property
+    def domain(self):
         """The domain of the linear operator - an element of Vectorspace"""
         return self._space
 
@@ -758,11 +761,6 @@ def codomain(self):
         """The codomain of the linear operator - an element of Vectorspace"""
         return self._codomain
 
-    @property
-    def domain(self):
-        """The domain of the linear operator - an element of Vectorspace"""
-        return self._domain
-
     @property
     def dtype(self):
         return self._dtype