GPU stack 8/8: profiling cases Poisson cube and Vlasov slab (CPU node vs GPU, GPU scaling)

examples/Poisson/gaussian_cube/params_gaussian_cube.py

@@ -37,10 +37,41 @@
     action="store_true",
     help="Small resolution (num_elements = [16, 16, 16]), for the tests (test_examples.py).",
 )
+parser.add_argument(
+    "--steps",
+    type=int,
+    default=1,
+    help="Number of time steps, i.e. Poisson solves with the same matrices (assembled once in the setup).",
+)
+parser.add_argument(
+    "--solver",
+    default="pcg",
+    choices=["pcg", "pcg_sr"],
+    help="CG variant: 'pcg_sr' does one global reduction per iteration instead of two (for many ranks/GPUs).",
+)
+parser.add_argument(
+    "--gpu-timing",
+    action="store_true",
+    help="Setup of the CPU vs GPU cases: one thread per rank, CUDA-event device times per region on the CuPy "
+    "backend instead of the line profiler, and a run label naming the backend and the rank count.",
+)
+parser.add_argument(
+    "--no-output",
+    action="store_true",
+    help="Write only the profiling data: no field or restart data (the pproc script then does nothing).",
+)
 args, _ = parser.parse_known_args()
 if args.test_run:
     args.num_elements = [16, 16, 16]
 
+if args.gpu_timing:
+    # One thread per MPI rank on both backends: OMP_NUM_THREADS=1, and the same for OpenBLAS/MKL, which would
+    # otherwise start one thread per core in every rank. Set before NumPy is imported (OpenBLAS reads it when loaded).
+    import os
+
+    for _var in ("OMP_NUM_THREADS", "OPENBLAS_NUM_THREADS", "MKL_NUM_THREADS"):
+        os.environ[_var] = "1"
+
 precond = None if args.precond == "None" else args.precond
 
 # preconditioner -> (name suffix, linear solver) for the name and the description
@@ -123,8 +154,8 @@
 model = Poisson(base_units=base_units)
 
 # List all variables and decide whether to save their data
-model.em_fields.phi.save_data = True
-model.em_fields.source.save_data = True
+model.em_fields.phi.save_data = not args.no_output
+model.em_fields.source.save_data = not args.no_output
 
 # --------------------------
 # Instance of the simulation
@@ -135,10 +166,14 @@
     out_folders=args.out_folder,
     sim_folder=args.sim_folder,
     restart=False,
+    save_step=args.steps,  # save the first and the last snapshot only
+    save_restart=not args.no_output,
 )
 
-# Time stepping
-time_opts = Time()
+# Time stepping: with --steps > 1 every step solves the Poisson problem again (from x0 = 0, recycle=False) with the
+# matrices assembled once in the setup, as in a time-dependent simulation; the solution is the same at every step.
+dt = 0.01
+time_opts = Time() if args.steps == 1 else Time(dt=dt, Tend=round(args.steps * dt, 14))
 
 # Geometry
 Lx = 2.0
@@ -155,10 +190,33 @@
 # Derham options
 derham_opts = DerhamOptions(degree=(1, 2, 3), bcs=(("dirichlet", "dirichlet"), None, None))
 
-# Profilinig options
-profiling_opts = ProfilingOptions(
-    use_line_profiler=True,
-)
+# Profiling options
+if args.gpu_timing:
+    # CUDA-event device times per region on the CuPy backend (not on cunumpy's fake CuPy, whose CUDA events have no
+    # elapsed time). No line profiler: it would distort the CPU/GPU comparison.
+    import cunumpy as xp
+    from cunumpy.kernel_testing import fake_cupy_active
+    from maybempi import MPI, SerialMPI
+
+    gpu_timing = xp.get_backend() == "cupy" and not fake_cupy_active()
+
+    # scope-profiler run label, the name of the run in the reports and comparison plots, e.g.
+    # "CPU: 28 MPI ranks" and "GPU: 4 MPI ranks, 4 GPUs"
+    num_ranks = 1 if isinstance(MPI, SerialMPI) else MPI.COMM_WORLD.Get_size()
+    profiling_label = (
+        f"{'GPU' if xp.get_backend() == 'cupy' else 'CPU'}: "
+        f"{num_ranks} MPI rank{'s' if num_ranks != 1 else ''}"
+        + (f", {num_ranks} GPU{'s' if num_ranks != 1 else ''}" if xp.get_backend() == "cupy" else "")
+    )
+    profiling_opts = ProfilingOptions(
+        label=profiling_label,
+        use_gpu_timing=gpu_timing,
+        gpu_timing_backend="cupy",
+    )
+else:
+    profiling_opts = ProfilingOptions(
+        use_line_profiler=True,
+    )
 
 # Simulation object
 sim = Simulation(
@@ -184,7 +242,7 @@
 solver_params = SolverParameters(tol=1e-8, maxiter=3000, info=True, recycle=False)
 model.propagators.poisson.options = model.propagators.poisson.Options(
     stab_eps=0.0,
-    solver="pcg",
+    solver=args.solver,
     precond=precond,
     solver_params=solver_params,
 )
@@ -236,4 +294,4 @@ def rhs_fun(e1, e2, e3):
 
 if __name__ == "__main__":
     # post-processing and plots: pproc_gaussian_cube.py
-    sim.run(profiling_activated=True, one_time_step=True)
+    sim.run(profiling_activated=True, one_time_step=args.steps == 1)

examples/Poisson/gaussian_cube/pproc_gaussian_cube.py

@@ -31,6 +31,10 @@
 
 
 def main(path_out=DEFAULT_OUTPUT, save_figs=True, angle=np.pi):
+    if args.no_output:
+        print("The run wrote no field data (--no-output): nothing to post-process.")
+        return
+
     out = Output(path_out).pproc(create_vtk=True, celldivide=CELLDIVIDE)
 
     # Raw FEEC fields are evaluated directly from the saved spline coefficients at the

examples/Poisson/gaussian_cube/submit_gaussian_cube_cpu_vs_gpu.py

@@ -0,0 +1,108 @@
+"""Poisson on a 3D cube: 28 CPU MPI ranks vs 1 GPU, on the same kind of Pitagora booster node.
+
+This file defines the profiling case (the `ProfilingCase`) and submits its two runs of the
+problem of `submit_gaussian_cube.py` (manufactured Gaussian solution on a Cuboid,
+degree (1, 2, 3), unpreconditioned CG), at 128 x 128 x 128 cells and with 5 time steps, i.e. 5
+Poisson solves with the mass matrices assembled once in the setup, as in a time-dependent simulation:
+
+- CPU: 28 MPI ranks with one thread each (OMP_NUM_THREADS=1), NumPy backend;
+- GPU: 1 MPI rank with one thread on one H100, CuPy backend.
+
+Every rank gets one CPU core (`cpus_per_task=1` of the preset), and the params file sets
+OMP_NUM_THREADS=1 (and OPENBLAS/MKL_NUM_THREADS=1) before NumPy is imported (`--gpu-timing`).
+
+Both runs use the `pitagora_booster` preset, so they run on identical hardware; the CPU run does not
+touch the GPUs. Without a batch system (or with STRUPHY_PROFILING_LOCAL=1) the runs execute directly
+on this machine, where the 28-rank run is skipped unless there are that many cores.
+`finalize_run` then packages and uploads each run as soon as its own job finishes.
+Each generated script runs the simulation itself by invoking `params_gaussian_cube.py`
+directly (its `__main__` block is the worker), and then post-processes it with `pproc_gaussian_cube.py`
+(errors of the RHS and the solution, figures).
+
+Why this size: the GPU's advantage grows with the grid (the CG solve is bound by memory bandwidth, the
+CPU's by its 28 cores); at 128^3 one solve takes 3.3 s on the GPU vs ~66 s on 28 login-node ranks,
+and the setup (mass-matrix assembly in CUDA) 1.4 s vs ~15 s.
+
+The runs are profiled with scope-profiler, the CuPy run with CUDA-event timing per region; the
+runs are labelled "CPU: 28 MPI ranks" and "GPU: 1 MPI rank, 1 GPU".
+"""
+
+import argparse
+import sys
+from pathlib import Path
+
+# the profiling driver (profiling_job.py, clusters.py, ...) lives in profiling/ at the repo root, three levels up
+sys.path.insert(0, str(Path(__file__).resolve().parents[3] / "profiling"))
+
+from profiling_job import ProfilingCase  # noqa: E402
+
+
+def main() -> None:
+
+    # Parse arguments, do not remove --upload
+    parser = argparse.ArgumentParser(
+        description=("Submit profiling jobs to a SLURM cluster and package the results for upload."),
+    )
+    parser.add_argument(
+        "--upload",
+        action="store_true",
+        help="Upload the packaged profiling results to the profiling-data repo.",
+    )
+    parser.add_argument(
+        "--language",
+        default="fortran",
+        choices=["fortran", "c"],
+        help="Pyccel language to compile the Struphy kernels with.",
+    )
+    parser.add_argument(
+        "--compiler",
+        default="GNU",
+        help='Pyccel compiler family ("GNU", "intel", "PGI", "nvidia" or "LLVM").',
+    )
+    args = parser.parse_args()
+
+    # Paths relative to this script's location, so it can be run from anywhere.
+    script_dir = Path(__file__).resolve().parent
+
+    profiling_case = ProfilingCase(
+        label="poisson_cube_cpu_node_vs_gpu",
+        name="Poisson on 3D cube: 28 CPU ranks vs 1 GPU",
+        description="The Poisson model with manufactured solution on a 3D cube (Cuboid, 128^3 cells, degree (1, 2, 3)), 5 solves with unpreconditioned CG, on 28 MPI ranks x 1 thread of a booster node (NumPy) vs 1 MPI rank x 1 thread on 1 H100 GPU (CuPy).",
+        physics_problem="Occurs in many plasma applications.",
+        struphy_model_used="Poisson",
+        params_source=script_dir / "params_gaussian_cube.py",
+        pproc_source=script_dir / "pproc_gaussian_cube.py",
+        language=args.language,
+        compiler=args.compiler,
+        upload=args.upload,
+    )
+
+    # 128^3 cells, 5 solves (see the module docstring); the CPU run takes a few minutes
+    param_flags = ["--num-elements", "128", "128", "128", "--steps", "5", "--gpu-timing"]
+    slurm_options = {"time": "00:30:00"}
+
+    # CPU: 28 MPI ranks x 1 thread on a booster node, NumPy backend
+    profiling_case.launch(
+        28,
+        num_nodes=1,
+        param_flags=param_flags,
+        backend="numpy",
+        slurm_preset="pitagora_booster",
+        slurm_options=slurm_options,
+    )
+    # GPU: 1 MPI rank x 1 thread on one H100 (the rank binds to GPU 0 of the node), CuPy backend
+    profiling_case.launch(
+        1,
+        num_nodes=1,
+        param_flags=param_flags,
+        backend="cupy",
+        slurm_preset="pitagora_booster",
+        slurm_options=slurm_options,
+    )
+
+    # Package and push each run as its own job finishes.
+    profiling_case.finalize_run()
+
+
+if __name__ == "__main__":
+    main()

examples/Poisson/gaussian_cube/submit_gaussian_cube_gpu_scaling.py

@@ -0,0 +1,117 @@
+"""Poisson on a 3D cube: strong scaling on 1, 2 and 4 GPUs of a Pitagora booster node.
+
+This file defines the profiling case (the `ProfilingCase`) and submits its three runs of the same
+problem as `submit_gaussian_cube_cpu_vs_gpu.py` (manufactured Gaussian solution on a Cuboid,
+degree (1, 2, 3), unpreconditioned CG, 5 solves), by default on 192 x 192 x 192 cells (`--cells`), on the
+CuPy backend with 1, 2 and 4 MPI ranks (`--gpus`), one GPU per rank (feectools binds each rank to the GPU of its
+node-local rank), one thread per rank. `--solver pcg_sr` selects the CG variant with one global reduction
+per iteration. The runs write only their profiling data (`--no-output` of the params file), so there is no
+post-processing. Each run uses the `pitagora_booster` preset (one node, its 4 GPUs).
+Without a batch system (or with STRUPHY_PROFILING_LOCAL=1) the runs execute directly on this machine.
+`finalize_run` then packages and uploads each run as soon as its own job finishes.
+Each generated script runs the simulation itself by invoking `params_gaussian_cube.py`
+directly (its `__main__` block is the worker).
+
+The runs are profiled with scope-profiler (CUDA-event timing per region) and labelled
+"GPU: N MPI ranks, N GPUs". The case label contains the size (and the solver if not pcg), e.g.
+`poisson_cube_gpu_scaling_192`.
+
+Recommended runs: the default (192^3 on 1, 2 and 4 GPUs) for the strong scaling from one GPU, and
+`--cells 256 --gpus 2 4` for a larger problem (256^3 does not fit on one H100: the setup assembles M0 and M1 as
+3D stencil matrices and evaluates the metric at all quadrature points at once). The fixed cost per CG
+iteration (ghost exchanges, reductions) weighs less on the larger problem: 4 vs 2 GPUs gave 1.30x at 192^3 and
+1.61x at 256^3 (Oct 2026).
+"""
+
+import argparse
+import sys
+from pathlib import Path
+
+# the profiling driver (profiling_job.py, clusters.py, ...) lives in profiling/ at the repo root, three levels up
+sys.path.insert(0, str(Path(__file__).resolve().parents[3] / "profiling"))
+
+from profiling_job import ProfilingCase  # noqa: E402
+
+
+def main() -> None:
+
+    # Parse arguments, do not remove --upload
+    parser = argparse.ArgumentParser(
+        description=("Submit profiling jobs to a SLURM cluster and package the results for upload."),
+    )
+    parser.add_argument(
+        "--upload",
+        action="store_true",
+        help="Upload the packaged profiling results to the profiling-data repo.",
+    )
+    parser.add_argument(
+        "--language",
+        default="fortran",
+        choices=["fortran", "c"],
+        help="Pyccel language to compile the Struphy kernels with.",
+    )
+    parser.add_argument(
+        "--compiler",
+        default="GNU",
+        help='Pyccel compiler family ("GNU", "intel", "PGI", "nvidia" or "LLVM").',
+    )
+    parser.add_argument(
+        "--cells",
+        type=int,
+        default=192,
+        help="Number of cells in each direction. 192^3 fits on one H100; 256^3 needs at least 2 (--gpus 2 4).",
+    )
+    parser.add_argument(
+        "--gpus",
+        type=int,
+        nargs="+",
+        default=[1, 2, 4],
+        choices=[1, 2, 3, 4],
+        help="GPU counts (MPI ranks, one GPU each) to run, on one booster node.",
+    )
+    parser.add_argument(
+        "--solver",
+        default="pcg",
+        choices=["pcg", "pcg_sr"],
+        help="CG variant: 'pcg_sr' does one global reduction per iteration instead of two.",
+    )
+    args = parser.parse_args()
+
+    # Paths relative to this script's location, so it can be run from anywhere.
+    script_dir = Path(__file__).resolve().parent
+    gpus = sorted(set(args.gpus))
+    gpus_text = ", ".join(str(n) for n in gpus[:-1]) + (" and " if len(gpus) > 1 else "") + str(gpus[-1])
+
+    profiling_case = ProfilingCase(
+        # one label per size and solver: runs of different sizes are separate cases in the profiling-data repo
+        label=f"poisson_cube_gpu_scaling_{args.cells}" + ("" if args.solver == "pcg" else f"_{args.solver}"),
+        name=f"Poisson on 3D cube, {args.cells}^3 cells: strong scaling on {gpus_text} GPUs",
+        description=f"The Poisson model with manufactured solution on a 3D cube (Cuboid, {args.cells}^3 cells, degree (1, 2, 3)), 5 solves with unpreconditioned CG ('{args.solver}'), on {gpus_text} H100 GPUs of a booster node (CuPy, one MPI rank per GPU).",
+        physics_problem="Occurs in many plasma applications.",
+        struphy_model_used="Poisson",
+        params_source=script_dir / "params_gaussian_cube.py",
+        language=args.language,
+        compiler=args.compiler,
+        upload=args.upload,
+    )
+
+    # 5 solves, as in submit_gaussian_cube_cpu_vs_gpu.py; only the profiling data is written
+    cells = [str(args.cells)] * 3
+    param_flags = ["--num-elements", *cells, "--steps", "5", "--solver", args.solver, "--gpu-timing", "--no-output"]
+
+    # one MPI rank per GPU, on one booster node
+    for num_tasks in gpus:
+        profiling_case.launch(
+            num_tasks,
+            num_nodes=1,
+            param_flags=param_flags,
+            backend="cupy",
+            slurm_preset="pitagora_booster",
+        )
+
+    # Package and push each run as its own job finishes.
+    profiling_case.finalize_run()
+
+
+if __name__ == "__main__":
+    main()