Added all changes from struphy-with-cupy

params_LinearMHDDriftkineticCC.py

@@ -0,0 +1,116 @@
+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
+
+# import model, set verbosity
+from struphy.models.hybrid import LinearMHDDriftkineticCC
+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,
+    )

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

@@ -409,7 +409,7 @@ def parse_path(directory):
         for filename in files:
             if re.search(r"__\w+__", root):
                 continue
-            if (filename.endswith(".py") or filename.endswith(".ipynb")) and not re.search(r"__\w+__", filename):
+            if filename.endswith(".py") and not re.search(r"__\w+__", filename):
                 file_path = os.path.join(root, filename)
                 python_files.append(file_path)
     # exit()
@@ -484,9 +484,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

@@ -42,14 +42,14 @@ def struphy_test(
                 str(mpi),
                 "pytest",
                 "-k",
-                "not _models and not _tutorial and not pproc and not _verif_",
+                "not _models and not _tutorial and not pproc",
                 "--with-mpi",
             ]
         else:
             cmd = [
                 "pytest",
                 "-k",
-                "not _models and not _tutorial and not pproc and not _verif_",
+                "not _models and not _tutorial and not pproc",
             ]
 
         if with_desc:

src/struphy/console/tests/test_console.py

@@ -7,6 +7,7 @@
 import pytest
 
 # from psydac.ddm.mpi import mpi as MPI
+import struphy
 import struphy as struphy_lib
 from struphy.console.compile import struphy_compile
 from struphy.console.main import struphy

src/struphy/diagnostics/diagnostics_pic.ipynb

@@ -7,13 +7,11 @@
    "outputs": [],
    "source": [
     "import os\n",
-    "\n",
+    "import struphy\n",
     "import numpy as np\n",
     "from matplotlib import pyplot as plt\n",
     "\n",
-    "import struphy\n",
-    "\n",
-    "path_out = os.path.join(struphy.__path__[0], \"io/out\", \"sim_1\")\n",
+    "path_out = os.path.join(struphy.__path__[0], 'io/out', 'sim_1')\n",
     "\n",
     "print(path_out)\n",
     "os.listdir(path_out)"
@@ -30,7 +28,7 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "data_path = os.path.join(path_out, \"post_processing\")\n",
+    "data_path = os.path.join(path_out, 'post_processing')\n",
     "\n",
     "os.listdir(data_path)"
    ]
@@ -41,7 +39,7 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "t_grid = np.load(os.path.join(data_path, \"t_grid.npy\"))\n",
+    "t_grid = np.load(os.path.join(data_path, 't_grid.npy'))\n",
     "t_grid"
    ]
   },
@@ -51,7 +49,7 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "f_path = os.path.join(data_path, \"kinetic_data\", \"ions\", \"distribution_function\")\n",
+    "f_path = os.path.join(data_path, 'kinetic_data', 'ions', 'distribution_function')\n",
     "\n",
     "print(os.listdir(f_path))"
    ]
@@ -62,7 +60,7 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "path = os.path.join(f_path, \"e1\")\n",
+    "path = os.path.join(f_path, 'e1')\n",
     "print(os.listdir(path))"
    ]
   },
@@ -72,9 +70,9 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "grid = np.load(os.path.join(f_path, \"e1/\", \"grid_e1.npy\"))\n",
-    "f_binned = np.load(os.path.join(f_path, \"e1/\", \"f_binned.npy\"))\n",
-    "delta_f_e1_binned = np.load(os.path.join(f_path, \"e1/\", \"delta_f_binned.npy\"))\n",
+    "grid = np.load(os.path.join(f_path, 'e1/', 'grid_e1.npy'))\n",
+    "f_binned = np.load(os.path.join(f_path, 'e1/', 'f_binned.npy'))\n",
+    "delta_f_e1_binned = np.load(os.path.join(f_path, 'e1/', 'delta_f_binned.npy'))\n",
     "\n",
     "print(grid.shape)\n",
     "print(f_binned.shape)\n",
@@ -89,18 +87,18 @@
    "source": [
     "steps = list(np.arange(10))\n",
     "\n",
-    "plt.figure(figsize=(12, 5 * len(steps)))\n",
+    "plt.figure(figsize=(12, 5*len(steps)))\n",
     "for n, step in enumerate(steps):\n",
-    "    plt.subplot(len(steps), 2, 2 * n + 1)\n",
-    "    plt.plot(grid, f_binned[step], label=f\"time = {t_grid[step]}\")\n",
-    "    plt.xlabel(\"e1\")\n",
-    "    # plt.ylim([.5, 1.5])\n",
-    "    plt.title(\"full-f\")\n",
-    "    plt.subplot(len(steps), 2, 2 * n + 2)\n",
-    "    plt.plot(grid, delta_f_e1_binned[step], label=f\"time = {t_grid[step]}\")\n",
-    "    plt.xlabel(\"e1\")\n",
-    "    # plt.ylim([-3e-3, 3e-3])\n",
-    "    plt.title(r\"$\\delta f$\")\n",
+    "    plt.subplot(len(steps), 2, 2*n + 1)\n",
+    "    plt.plot(grid, f_binned[step], label=f'time = {t_grid[step]}')\n",
+    "    plt.xlabel('e1')\n",
+    "    #plt.ylim([.5, 1.5])\n",
+    "    plt.title('full-f')\n",
+    "    plt.subplot(len(steps), 2, 2*n + 2)\n",
+    "    plt.plot(grid, delta_f_e1_binned[step], label=f'time = {t_grid[step]}')\n",
+    "    plt.xlabel('e1')\n",
+    "    #plt.ylim([-3e-3, 3e-3])\n",
+    "    plt.title(r'$\\delta f$')\n",
     "    plt.legend()"
    ]
   },
@@ -110,7 +108,7 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "path = os.path.join(f_path, \"e1_v1\")\n",
+    "path = os.path.join(f_path, 'e1_v1')\n",
     "print(os.listdir(path))"
    ]
   },
@@ -120,10 +118,10 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "grid_e1 = np.load(os.path.join(f_path, \"e1_v1/\", \"grid_e1.npy\"))\n",
-    "grid_v1 = np.load(os.path.join(f_path, \"e1_v1/\", \"grid_v1.npy\"))\n",
-    "f_binned = np.load(os.path.join(f_path, \"e1_v1/\", \"f_binned.npy\"))\n",
-    "delta_f_binned = np.load(os.path.join(f_path, \"e1_v1/\", \"delta_f_binned.npy\"))\n",
+    "grid_e1 = np.load(os.path.join(f_path, 'e1_v1/', 'grid_e1.npy'))\n",
+    "grid_v1 = np.load(os.path.join(f_path, 'e1_v1/', 'grid_v1.npy'))\n",
+    "f_binned = np.load(os.path.join(f_path, 'e1_v1/', 'f_binned.npy'))\n",
+    "delta_f_binned = np.load(os.path.join(f_path, 'e1_v1/', 'delta_f_binned.npy'))\n",
     "\n",
     "print(grid_e1.shape)\n",
     "print(grid_v1.shape)\n",
@@ -139,20 +137,20 @@
    "source": [
     "steps = list(np.arange(10))\n",
     "\n",
-    "plt.figure(figsize=(12, 5 * len(steps)))\n",
+    "plt.figure(figsize=(12, 5*len(steps)))\n",
     "for n, step in enumerate(steps):\n",
-    "    plt.subplot(len(steps), 2, 2 * n + 1)\n",
-    "    plt.pcolor(grid_e1, grid_v1, f_binned[step].T, label=f\"time = {t_grid[step]}\")\n",
-    "    plt.xlabel(\"$e1$\")\n",
-    "    plt.ylabel(r\"$v_\\parallel$\")\n",
-    "    plt.title(\"full-f\")\n",
+    "    plt.subplot(len(steps), 2, 2*n + 1)\n",
+    "    plt.pcolor(grid_e1, grid_v1, f_binned[step].T, label=f'time = {t_grid[step]}')\n",
+    "    plt.xlabel('$e1$')\n",
+    "    plt.ylabel(r'$v_\\parallel$')\n",
+    "    plt.title('full-f')\n",
     "    plt.legend()\n",
     "    plt.colorbar()\n",
-    "    plt.subplot(len(steps), 2, 2 * n + 2)\n",
-    "    plt.pcolor(grid_e1, grid_v1, delta_f_binned[step].T, label=f\"time = {t_grid[step]}\")\n",
-    "    plt.xlabel(\"$e1$\")\n",
-    "    plt.ylabel(r\"$v_\\parallel$\")\n",
-    "    plt.title(r\"$\\delta f$\")\n",
+    "    plt.subplot(len(steps), 2, 2*n + 2)\n",
+    "    plt.pcolor(grid_e1, grid_v1, delta_f_binned[step].T, label=f'time = {t_grid[step]}')\n",
+    "    plt.xlabel('$e1$')\n",
+    "    plt.ylabel(r'$v_\\parallel$')\n",
+    "    plt.title(r'$\\delta f$')\n",
     "    plt.legend()\n",
     "    plt.colorbar()"
    ]
@@ -163,7 +161,7 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "fields_path = os.path.join(data_path, \"fields_data\")\n",
+    "fields_path = os.path.join(data_path, 'fields_data')\n",
     "\n",
     "print(os.listdir(fields_path))"
    ]
@@ -176,7 +174,7 @@
    "source": [
     "import pickle\n",
     "\n",
-    "with open(os.path.join(fields_path, \"grids_phy.bin\"), \"rb\") as file:\n",
+    "with open(os.path.join(fields_path, 'grids_phy.bin'), 'rb') as file:\n",
     "    x_grid, y_grid, z_grid = pickle.load(file)\n",
     "\n",
     "print(type(x_grid))\n",
@@ -189,7 +187,7 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "with open(os.path.join(fields_path, \"em_fields\", \"phi_phy.bin\"), \"rb\") as file:\n",
+    "with open(os.path.join(fields_path, 'em_fields', 'phi_phy.bin'), 'rb') as file:\n",
     "    phi = pickle.load(file)\n",
     "\n",
     "plt.figure(figsize=(12, 12))\n",
@@ -199,9 +197,9 @@
     "    t = t_grid[step]\n",
     "    print(phi[t][0].shape)\n",
     "    plt.subplot(2, 2, n + 1)\n",
-    "    plt.plot(x_grid[:, 0, 0], phi[t][0][:, 0, 0], label=f\"time = {t}\")\n",
-    "    plt.xlabel(\"x\")\n",
-    "    plt.ylabel(r\"$\\phi$(x)\")\n",
+    "    plt.plot(x_grid[:, 0, 0], phi[t][0][:, 0, 0], label=f'time = {t}')\n",
+    "    plt.xlabel('x')\n",
+    "    plt.ylabel(r'$\\phi$(x)')\n",
     "    plt.legend()"
    ]
   },

src/struphy/feec/basis_projection_ops.py

@@ -51,7 +51,7 @@ def __init__(self, derham, domain, verbose=True, **weights):
 
         self._rank = derham.comm.Get_rank() if derham.comm is not None else 0
 
-        if xp.any([p == 1 and Nel > 1 for p, Nel in zip(derham.p, derham.Nel)]):
+        if xp.any(xp.array([p == 1 and Nel > 1 for p, Nel in zip(derham.p, derham.Nel)])):
             if MPI.COMM_WORLD.Get_rank() == 0:
                 print(
                     f'\nWARNING: Class "BasisProjectionOperators" called with p={derham.p} (interpolation of piece-wise constants should be avoided).',

src/struphy/feec/preconditioner.py

@@ -318,6 +318,11 @@ 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"""
@@ -699,6 +704,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