From 96011ba1cb47e8279a0a65d3e1a9e8f1421f7a2b Mon Sep 17 00:00:00 2001 From: Max Date: Wed, 7 Oct 2026 22:59:19 +0200 Subject: [PATCH] Move FEEC assembly kernels into kernel folders The 16 entry kernels of feec/mass_kernels.py and feec/basis_projection_kernels.py (used by WeightedMassOperator, StencilMatrixFreeMassOperator, L2Projector, BoundaryMassOperator and BasisProjectionOperator) move into one folder each under feec/kernels// with _kernels.py and an __init__.py declaring the cunumpy Kernel. Call sites import the kernels instead of wrapping PyccelKernel(getattr(module, name)). The function bodies are unchanged. The kernels have no CUDA version yet and are declared with missing_cuda="fallback", which keeps the former PyccelKernel behaviour on the CuPy backend. Solves #411. Co-Authored-By: Claude Opus 5.5 --- CUDA_STRATEGY.md | 20 +- src/struphy/feec/basis_projection_kernels.py | 568 -------- src/struphy/feec/basis_projection_ops.py | 19 +- src/struphy/feec/boundary_mass.py | 5 +- src/struphy/feec/kernels/__init__.py | 10 +- .../__init__.py | 11 + ..._dofs_for_weighted_basisfuns_1d_kernels.py | 128 ++ .../__init__.py | 11 + ..._dofs_for_weighted_basisfuns_2d_kernels.py | 186 +++ .../__init__.py | 11 + ..._dofs_for_weighted_basisfuns_3d_kernels.py | 259 ++++ .../feec/kernels/kernel_1d_eval/__init__.py | 11 + .../kernel_1d_eval/kernel_1d_eval_kernels.py | 55 + .../feec/kernels/kernel_1d_mat/__init__.py | 11 + .../kernel_1d_mat/kernel_1d_mat_kernels.py | 68 + .../feec/kernels/kernel_1d_vec/__init__.py | 11 + .../kernel_1d_vec/kernel_1d_vec_kernels.py | 59 + .../feec/kernels/kernel_2d_eval/__init__.py | 11 + .../kernel_2d_eval/kernel_2d_eval_kernels.py | 72 + .../feec/kernels/kernel_2d_mat/__init__.py | 11 + .../kernel_2d_mat/kernel_2d_mat_kernels.py | 88 ++ .../feec/kernels/kernel_2d_vec/__init__.py | 11 + .../kernel_2d_vec/kernel_2d_vec_kernels.py | 76 + .../feec/kernels/kernel_3d_diag/__init__.py | 11 + .../kernel_3d_diag/kernel_3d_diag_kernels.py | 118 ++ .../feec/kernels/kernel_3d_eval/__init__.py | 11 + .../kernel_3d_eval/kernel_3d_eval_kernels.py | 85 ++ .../feec/kernels/kernel_3d_mat/__init__.py | 11 + .../kernel_3d_mat/kernel_3d_mat_kernels.py | 147 ++ .../kernels/kernel_3d_matrixfree/__init__.py | 11 + .../kernel_3d_matrixfree_kernels.py | 173 +++ .../feec/kernels/kernel_3d_vec/__init__.py | 11 + .../kernel_3d_vec/kernel_3d_vec_kernels.py | 96 ++ .../kernels/surface_kernel_3d_mat/__init__.py | 11 + .../surface_kernel_3d_mat_kernels.py | 170 +++ .../kernels/surface_kernel_3d_vec/__init__.py | 11 + .../surface_kernel_3d_vec_kernels.py | 91 ++ src/struphy/feec/mass.py | 47 +- src/struphy/feec/mass_kernels.py | 1270 ----------------- src/struphy/pic/tests/test_kernel_backends.py | 2 +- 40 files changed, 2111 insertions(+), 1877 deletions(-) delete mode 100644 src/struphy/feec/basis_projection_kernels.py create mode 100644 src/struphy/feec/kernels/assemble_dofs_for_weighted_basisfuns_1d/__init__.py create mode 100644 src/struphy/feec/kernels/assemble_dofs_for_weighted_basisfuns_1d/assemble_dofs_for_weighted_basisfuns_1d_kernels.py create mode 100644 src/struphy/feec/kernels/assemble_dofs_for_weighted_basisfuns_2d/__init__.py create mode 100644 src/struphy/feec/kernels/assemble_dofs_for_weighted_basisfuns_2d/assemble_dofs_for_weighted_basisfuns_2d_kernels.py create mode 100644 src/struphy/feec/kernels/assemble_dofs_for_weighted_basisfuns_3d/__init__.py create mode 100644 src/struphy/feec/kernels/assemble_dofs_for_weighted_basisfuns_3d/assemble_dofs_for_weighted_basisfuns_3d_kernels.py create mode 100644 src/struphy/feec/kernels/kernel_1d_eval/__init__.py create mode 100644 src/struphy/feec/kernels/kernel_1d_eval/kernel_1d_eval_kernels.py create mode 100644 src/struphy/feec/kernels/kernel_1d_mat/__init__.py create mode 100644 src/struphy/feec/kernels/kernel_1d_mat/kernel_1d_mat_kernels.py create mode 100644 src/struphy/feec/kernels/kernel_1d_vec/__init__.py create mode 100644 src/struphy/feec/kernels/kernel_1d_vec/kernel_1d_vec_kernels.py create mode 100644 src/struphy/feec/kernels/kernel_2d_eval/__init__.py create mode 100644 src/struphy/feec/kernels/kernel_2d_eval/kernel_2d_eval_kernels.py create mode 100644 src/struphy/feec/kernels/kernel_2d_mat/__init__.py create mode 100644 src/struphy/feec/kernels/kernel_2d_mat/kernel_2d_mat_kernels.py create mode 100644 src/struphy/feec/kernels/kernel_2d_vec/__init__.py create mode 100644 src/struphy/feec/kernels/kernel_2d_vec/kernel_2d_vec_kernels.py create mode 100644 src/struphy/feec/kernels/kernel_3d_diag/__init__.py create mode 100644 src/struphy/feec/kernels/kernel_3d_diag/kernel_3d_diag_kernels.py create mode 100644 src/struphy/feec/kernels/kernel_3d_eval/__init__.py create mode 100644 src/struphy/feec/kernels/kernel_3d_eval/kernel_3d_eval_kernels.py create mode 100644 src/struphy/feec/kernels/kernel_3d_mat/__init__.py create mode 100644 src/struphy/feec/kernels/kernel_3d_mat/kernel_3d_mat_kernels.py create mode 100644 src/struphy/feec/kernels/kernel_3d_matrixfree/__init__.py create mode 100644 src/struphy/feec/kernels/kernel_3d_matrixfree/kernel_3d_matrixfree_kernels.py create mode 100644 src/struphy/feec/kernels/kernel_3d_vec/__init__.py create mode 100644 src/struphy/feec/kernels/kernel_3d_vec/kernel_3d_vec_kernels.py create mode 100644 src/struphy/feec/kernels/surface_kernel_3d_mat/__init__.py create mode 100644 src/struphy/feec/kernels/surface_kernel_3d_mat/surface_kernel_3d_mat_kernels.py create mode 100644 src/struphy/feec/kernels/surface_kernel_3d_vec/__init__.py create mode 100644 src/struphy/feec/kernels/surface_kernel_3d_vec/surface_kernel_3d_vec_kernels.py delete mode 100644 src/struphy/feec/mass_kernels.py diff --git a/CUDA_STRATEGY.md b/CUDA_STRATEGY.md index daf0d25f0..02d81ca0e 100644 --- a/CUDA_STRATEGY.md +++ b/CUDA_STRATEGY.md @@ -26,6 +26,7 @@ Progress is tracked in struphy-hub/struphy#650. - [ ] **PR 17: struphy on the feectools CUDA stack** — the `feectools` submodule points at the top of the feectools CUDA stack (`cuda-4-device-kernels`, [feectools#88](https://github.com/struphy-hub/feectools/pull/88)) instead of `devel-tiny`, so the struphy CUDA PRs run against feectools with device stencil operations, MPI with device buffers and one GPU per rank. Moves along with the stack; before the struphy stack is merged into `devel`, the feectools stack is merged into `devel-tiny` and the submodule points there again (see [feectools](#feectools)). - [x] **PR 18: Geometry evaluation on the GPU for all analytic mappings** (#681, open) — CUDA versions of the four geometry entry kernels (`kernel_evaluate_pic`, `kernel_evaluate`, `kernel_pullpush_pic`, `kernel_pullpush` in `geometry/kernels/`), built on device versions of the whole metric chain (`f`, `df`, `det_df`, `df_inv`, `g`, `g_inv`, `select_metric_coeff`, `pull`/`push`/`tran`) for every analytic mapping (`kind_map` 10–12, 20–22, 30–32). Restores CuPy particle runs (weight initialization evaluates `jacobian_det`), and removes the Cuboid-only checks in `Pusher`, the accumulators and `reflect`. Parity arguments cover every analytic mapping (see [PR 18](#pr-18-geometry-evaluation-for-all-analytic-mappings) and the [PR 18 implementation notes](#pr-18-implementation-notes)). - [x] **PR 19: Spline mappings on the GPU** ([#683](https://github.com/struphy-hub/struphy/pull/683)) — `kind_map` 0–2 (`IGAPolarCylinder`, `IGAPolarTorus`, `Tokamak`, GVEC, DESC): `DomainArgs` gets array views for `t1..3`, `ind1..3` and `cx/cy/cz` (shapes needed on the device), spline-mapped `Domain`s can be created on the CuPy backend (control points fitted on the host), and `spline_3d`, `spline_2d_straight`, `spline_2d_torus` get device versions in the `kind_map` switch, so every mapping runs on CuPy and `check_mapping_on_device` is gone. Parity arguments and device-helper tests add four spline mappings. Polar splines and `EQDSKequilibrium` stay host-only (see [PR 19](#pr-19-spline-mappings) and the [PR 19 implementation notes](#pr-19-implementation-notes)). +- [x] **FEEC assembly kernels in kernel folders** (#411) — the 16 assembly kernels of `feec/mass.py`, `feec/boundary_mass.py` and `feec/basis_projection_ops.py` move from `feec/mass_kernels.py` and `feec/basis_projection_kernels.py` into one folder each under `feec/kernels/`; pure restructuring, ready for CUDA ports (see [FEEC assembly kernels](#feec-assembly-kernels-411)). - [ ] **Next (order to be confirmed)**: 6D array views in cunumpy and the blocked matrix accumulations `vlasov_maxwell`, `linear_vlasov_ampere` (steps 1–2 of the [porting order](#porting-order)); one complete model (`VlasovAmpereOneSpecies`) on the GPU end to end, including when to compile the kernels, which needs the [feectools stack](#feectools) merged and released first; a decision on hand-written CUDA vs. code generation before the guiding-center kernels (step 3). - [ ] **feectools**: the FEEC side (stencil vectors and matrices, MPI exchange, GPU binding) in feectools, stacked PRs [#85](https://github.com/struphy-hub/feectools/pull/85) (merged into `cuda-development`), [#86](https://github.com/struphy-hub/feectools/pull/86), [#87](https://github.com/struphy-hub/feectools/pull/87), [#88](https://github.com/struphy-hub/feectools/pull/88), integrated by [#90](https://github.com/struphy-hub/feectools/pull/90); needed before the end-to-end model run, not for PR 18/19 (see [feectools](#feectools)). - [x] **PR 13: Move the kernel infrastructure to cunumpy** — `Kernel`, `KernelCatalog`, `CudaKernel` and `Argument` are replaced by their cunumpy counterparts; each owner has a single `args_*` object on both backends, and `pusher_args.cuh` is generated (see [Moving to cunumpy](#moving-to-cunumpy-pr-13)). Kernels and device helpers only change their includes. @@ -64,7 +65,7 @@ CUDA kernels can be added one by one. If the code runs on the GPU and needs a ke | `src/struphy/utils/cuda_arguments.py` | `CUDA_STRUCTS`, `CUDA_OPTIONS`, `write_pusher_header()` and `write_local_projectors_header()` | | `src/struphy/kernel_arguments/pusher_args.cuh`, `local_projectors_args.cuh` | the C structs `MarkerArgs`, `DerhamArgs`, `DomainArgs` and `LocalProjectorsArgs`, generated from the CUDA classes, using cunumpy array views | | `src/struphy/geometry/base.py`, `src/struphy/pic/base.py`, `src/struphy/feec/psydac_derham.py` | `Domain.args_domain`, `Particles.args_markers` and `Derham.args_derham` are the pyccel class on NumPy and the CUDA class on CuPy, chosen once at construction; every kernel call goes through a `Kernel` | -| `src/struphy/*/kernels/` (`pic/pushing`, `pic/accumulation`, `pic/diagnostics`, `pic/sph`, `bsplines`, `geometry`, `feec`, `feec/local_projectors`) | every kernel called from Python with argument objects, one folder each: 44 pusher/evaluation (incl. `reflect`), 16 accumulation, 10 marker diagnostics, 4 SPH evaluation, 3 spline evaluation, 4 geometry, 1 FEEC utility and 8 local projector kernels. Each folder's `__init__.py` declares its `Kernel`; the code imports it (`from struphy.geometry.kernels.kernel_evaluate import kernel_evaluate`). Fifteen have CUDA versions, with parity cases in `pic/tests/cuda_parity_cases.py`, among them the four geometry kernels (PR 18) | +| `src/struphy/*/kernels/` (`pic/pushing`, `pic/accumulation`, `pic/diagnostics`, `pic/sph`, `bsplines`, `geometry`, `feec`, `feec/local_projectors`) | every kernel called from Python with argument objects, one folder each: 44 pusher/evaluation (incl. `reflect`), 16 accumulation, 10 marker diagnostics, 4 SPH evaluation, 3 spline evaluation, 4 geometry, 1 FEEC utility, 16 FEEC assembly (#411) and 8 local projector kernels. Each folder's `__init__.py` declares its `Kernel`; the code imports it (`from struphy.geometry.kernels.kernel_evaluate import kernel_evaluate`). Fifteen have CUDA versions, with parity cases in `pic/tests/cuda_parity_cases.py`, among them the four geometry kernels (PR 18) | | `src/struphy/geometry/evaluation_kernels.cuh`, `transform_kernels.cuh`, `spline_mappings_kernels.cuh`, `geometry/domains//_cuda.cuh` | device versions of the mapping helpers for every mapping: `f`/`df` per analytic domain (`kind_map` 10–12, 20–22, 30–32, PR 18) and the spline mappings `spline_3d`, `spline_2d_straight`, `spline_2d_torus` (`kind_map` 0–2, PR 19), the `kind_map` switch, `det_df`, `df_inv`, `g`, `g_inv`, `select_metric_coeff`, `pull`, `push`, `tran` | | `src/struphy/pic/tests/test_cuda_parity.py`, `test_cuda_emulation.py` | parity of every CUDA kernel with pyccel on the cases of `pic/tests/cuda_parity_cases.py` (`PARITY_CASES`): on a GPU (`assert_kernels_agree`), and without one by CPU emulation | | `src/struphy/pic/tests/test_kernel_backends.py` | the generated-header tests, the check that each CUDA class mirrors its pyccel class, signature check of all kernels, owner class selection, and (on a GPU) struct layout and copy/pickle checks | @@ -626,3 +627,20 @@ Every mapping runs on the GPU: the spline mappings (`kind_map` 0–2) join the a The GPU tests have not run on an H100; the CPU emulation and the CPU regression tests are the gate. - **Parity cases.** The spline mappings are added to the geometry kernels' entries of `pic/tests/cuda_parity_cases.py` (the per-folder `_test_args.py` modules were replaced by that module in PR 16). + +## FEEC assembly kernels (#411) + +The assembly kernels of the weighted mass operators, the boundary mass operators and the basis projection +operators are entry kernels like the pushers, so they follow the folder rule: `feec/kernels//` with +`_kernels.py` and an `__init__.py` that declares the `Kernel`. Moved, without numerical change: +`kernel_{1,2,3}d_mat`, `kernel_{1,2,3}d_vec`, `kernel_{1,2,3}d_eval`, `kernel_3d_matrixfree`, `kernel_3d_diag`, +`surface_kernel_3d_mat`, `surface_kernel_3d_vec` (from `feec/mass_kernels.py`) and +`assemble_dofs_for_weighted_basisfuns_{1,2,3}d` (from `feec/basis_projection_kernels.py`); both modules are +gone. The call sites import the kernels and pick the one for the dimension from a small dict +(`_MAT_KERNELS[ldim]`, ...) instead of `PyccelKernel(getattr(module, name))`. + +None of them has a CUDA version yet. They are declared with `missing_cuda="fallback"`, which keeps the behaviour +of the former `PyccelKernel` calls on the CuPy backend (host copies of the arrays at every call, with a warning +once) instead of raising, so FEEC setups that worked on CuPy keep working. **They are ready for CUDA ports**: a +port adds `_cuda.cu` to the folder and its parity cases to `pic/tests/cuda_parity_cases.py`, and drops the +fallback for that kernel. They take only arrays and scalars, no argument objects. diff --git a/src/struphy/feec/basis_projection_kernels.py b/src/struphy/feec/basis_projection_kernels.py deleted file mode 100644 index a7d1ca1a1..000000000 --- a/src/struphy/feec/basis_projection_kernels.py +++ /dev/null @@ -1,568 +0,0 @@ -def assemble_dofs_for_weighted_basisfuns_1d( - mat: "float[:,:]", - starts_in: "int[:]", - ends_in: "int[:]", - pads_in: "int[:]", - starts_out: "int[:]", - ends_out: "int[:]", - pads_out: "int[:]", - fun_q: "float[:]", - wts1: "float[:,:]", - span1: "int[:,:]", - basis1: "float[:,:,:]", - sub1: "int[:]", - dim1_in: int, - dim1_out: int, - p1_out: int, -): - """Kernel for assembling the matrix - - A_(i,j) = DOFS_i(fun*Lambda^in_j) , - - into the _data attribute of a StencilMatrix. - Here, DOFS_i are the degrees-of-freedom of the output space (codomain, must not be a product space), - Lambda^in_j are the basis functions of the input space (domain, must not be a product space), and fun is an arbitrary function. - - Parameters - ---------- - mat : 2d float array - _data attribute of StencilMatrix. - - starts_in : int - Starting index of the input space (domain) of a distributed StencilMatrix. - - ends_in : int - Ending index of the input space (domain) of a distributed StencilMatrix. - - pads_in : int - Paddings of the input space (domain) of a distributed StencilMatrix. - - starts_out : int - Starting indices of the output space (codomain) of a distributed StencilMatrix. - - ends_out : int - Ending indices of the output space (codomain) of a distributed StencilMatrix. - - pads_out : int - Paddings of the output space (codomain) of a distributed StencilMatrix. - - fun_q : 1d float array - The function evaluated at the points (nq*ii + iq), where iq a local quadrature point of interval ii. - - wts1 : 2d float array - Quadrature weights in format (ii, iq). - - span1 : 2d int array - Knot span indices in direction eta1 in format (ii, iq). - - basis1 : 3d float array - Values of p1 + 1 non-zero eta-1 basis functions at quadrature points in format (ii, iq, basis function). - - sub1 : 1d int array - Sub-interval indices in direction 1. - - dim1_in : int - Dimension of the first direction of the input space - - p1_out : int - Spline degree of the first direction of the output space - """ - - # Start/end indices and paddings for distributed stencil matrix of input space - # si1 = starts_in[0} - # ei1 = ends_in[0] - pi1 = pads_in[0] - - # Start/end indices for distributed stencil matrix of output space - so1 = starts_out[0] - # eo1 = ends_out[0] - po1 = pads_out[0] - - # Spline degrees of input space - p1 = basis1.shape[2] - 1 - - # number of quadrature points - nq1 = span1.shape[1] - - # Set output to zero - mat[:] = 0.0 - - # Find column index for _data: - if dim1_out <= dim1_in: - cut1 = p1 - else: - cut1 = p1_out - - # Interval (either element or sub-interval thereof) - # ------------------------------------------------- - cumsub_i = 0 # Cumulative sub-interval index - for ii in range(span1.shape[0]): - cumsub_i += sub1[ii] - i = ii - cumsub_i # local DOF index - - # Quadrature point index in interval - # ---------------------------------- - for iq in range(nq1): - funval = fun_q[nq1 * ii + iq] * wts1[ii, iq] - - # Basis function of input space: - # ------------------------------ - for b1 in range(p1 + 1): - m = span1[ii, iq] - p1 + b1 # global index - # basis value - value = funval * basis1[ii, iq, b1] - - # Diff of global indices, needs to be adjusted for boundary conditions --> col1 - col1_tmp = m - (i + so1) - if col1_tmp > cut1: - m = m - dim1_in - elif col1_tmp < -cut1: - m = m + dim1_in - # add padding - col1 = pi1 + m - (i + so1) - - # Row index: padding + local index. - mat[po1 + i, col1] += value - - -def assemble_dofs_for_weighted_basisfuns_2d( - mat: "float[:,:,:,:]", - starts_in: "int[:]", - ends_in: "int[:]", - pads_in: "int[:]", - starts_out: "int[:]", - ends_out: "int[:]", - pads_out: "int[:]", - fun_q: "float[:,:]", - wts1: "float[:,:]", - wts2: "float[:,:]", - span1: "int[:,:]", - span2: "int[:,:]", - basis1: "float[:,:,:]", - basis2: "float[:,:,:]", - sub1: "int[:]", - sub2: "int[:]", - dim1_in: int, - dim2_in: int, - dim1_out: int, - dim2_out: int, - p1_out: int, - p2_out: int, -): - """Kernel for assembling the matrix - - A_(ij,kl) = DOFS_ij(fun*Lambda^in_kl) , - - into the _data attribute of a StencilMatrix. - Here, DOFS_ij are the degrees-of-freedom of the output space (codomain, must not be a product space), - Lambda^in_kl are the basis functions of the input space (domain, must not be a product space), and fun is an arbitrary function. - - Parameters - ---------- - mat : 4d float array - _data attribute of StencilMatrix. - - starts_in : 1d int array - Starting indices of the input space (domain) of a distributed StencilMatrix. - - ends_in : 1d int array - Ending indices of the input space (domain) of a distributed StencilMatrix. - - pads_in : 1d int array - Paddings of the input space (domain) of a distributed StencilMatrix. - - starts_out : 1d int array - Starting indices of the output space (codomain) of a distributed StencilMatrix. - - ends_out : 1d int array - Ending indices of the output space (codomain) of a distributed StencilMatrix. - - pads_out : 1d int array - Paddings of the output space (codomain) of a distributed StencilMatrix. - - fun_q : 2d float array - The function evaluated at the points (nq_i*ii + iq, nq_j*jj + jq), where iq a local quadrature point of interval ii. - - wts1 : 2d float array - Quadrature weights in direction eta1 in format (ii, iq). - - wts2 : 2d float array - Quadrature weights in direction eta2 in format (jj, jq). - - span1 : 2d int array - Knot span indices in direction eta1 in format (ii, iq). - - span2 : 2d int array - Knot span indices in direction eta2 in format (jj, jq). - - basis1 : 3d float array - Values of p1 + 1 non-zero eta-1 basis functions at quadrature points in format (ii, iq, basis function). - - basis2 : 3d float array - Values of p2 + 1 non-zero eta-2 basis functions at quadrature points in format (jj, jq, basis function). - - sub1 : 1d int array - Sub-interval indices in direction 1. - - sub2 : 1d int array - Sub-interval indices in direction 2. - - dim1_in : int - Dimension of the first direction of the input space - - dim2_in : int - Dimension of the second direction of the input space - - p1_out : int - Spline degree of the first direction of the output space - - p2_out : int - Spline degree of the second direction of the output space - """ - - # Start/end indices and paddings for distributed stencil matrix of input space - # si1 = starts_in[0] - # si2 = starts_in[1] - # ei1 = ends_in[0] - # ei2 = ends_in[1] - pi1 = pads_in[0] - pi2 = pads_in[1] - - # Start/end indices for distributed stencil matrix of output space - so1 = starts_out[0] - so2 = starts_out[1] - # eo1 = ends_out[0] - # eo2 = ends_out[1] - po1 = pads_out[0] - po2 = pads_out[1] - - # Spline degrees of input space - p1 = basis1.shape[2] - 1 - p2 = basis2.shape[2] - 1 - - # number of quadrature points - nq1 = span1.shape[1] - nq2 = span2.shape[1] - - # Set output to zero - mat[:] = 0.0 - - # Find column index for _data: - if dim1_out <= dim1_in: - cut1 = p1 - else: - cut1 = p1_out - - if dim2_out <= dim2_in: - cut2 = p2 - else: - cut2 = p2_out - - # Interval (either element or sub-interval thereof) - # ------------------------------------------------- - cumsub_i = 0 # Cumulative sub-interval index - for ii in range(span1.shape[0]): - cumsub_i += sub1[ii] - i = ii - cumsub_i # local DOF index - - cumsub_j = 0 # Cumulative sub-interval index - for jj in range(span2.shape[0]): - cumsub_j += sub2[jj] - j = jj - cumsub_j # local DOF index - - # Quadrature point index in interval - # ---------------------------------- - for iq in range(nq1): - for jq in range(nq2): - funval = fun_q[nq1 * ii + iq, nq2 * jj + jq] * wts1[ii, iq] * wts2[jj, jq] - - # Basis function of input space: - # ------------------------------ - for b1 in range(p1 + 1): - m = span1[ii, iq] - p1 + b1 # global index - # basis value - val1 = funval * basis1[ii, iq, b1] - - # Diff of global indices, needs to be adjusted for boundary conditions --> col1 - col1_tmp = m - (i + so1) - if col1_tmp > cut1: - m = m - dim1_in - elif col1_tmp < -cut1: - m = m + dim1_in - # add padding - col1 = pi1 + m - (i + so1) - - for b2 in range(p2 + 1): - # global index - n = span2[jj, jq] - p2 + b2 - value = val1 * basis2[jj, jq, b2] - - # Diff of global indices, needs to be adjusted for boundary conditions --> col2 - col2_tmp = n - (j + so2) - if col2_tmp > cut2: - n = n - dim2_in - elif col2_tmp < -cut2: - n = n + dim2_in - # add padding - col2 = pi2 + n - (j + so2) - - # Row index: padding + local index. - mat[po1 + i, po2 + j, col1, col2] += value - - -def assemble_dofs_for_weighted_basisfuns_3d( - mat: "float[:,:,:,:,:,:]", - starts_in: "int[:]", - ends_in: "int[:]", - pads_in: "int[:]", - starts_out: "int[:]", - ends_out: "int[:]", - pads_out: "int[:]", - fun_q: "float[:,:,:]", - wts1: "float[:,:]", - wts2: "float[:,:]", - wts3: "float[:,:]", - span1: "int[:,:]", - span2: "int[:,:]", - span3: "int[:,:]", - basis1: "float[:,:,:]", - basis2: "float[:,:,:]", - basis3: "float[:,:,:]", - sub1: "int[:]", - sub2: "int[:]", - sub3: "int[:]", - dim1_in: int, - dim2_in: int, - dim3_in: int, - dim1_out: int, - dim2_out: int, - dim3_out: int, - p1_out: int, - p2_out: int, - p3_out: int, -): - """Kernel for assembling the matrix - - A_(ijk,mno) = DOFS_ijk(fun*Lambda^in_mno) , - - into the _data attribute of a StencilMatrix. - Here, DOFS_ijk are the degrees-of-freedom of the output space (codomain, must not be a product space), - Lambda^in_mno are the basis functions of the input space (domain, must not be a product space), and fun is an arbitrary function. - - Parameters - ---------- - mat : 6d float array - _data attribute of StencilMatrix. - - starts_in : 1d int array - Starting indices of the input space (domain) of a distributed StencilMatrix. - - ends_in : 1d int array - Ending indices of the input space (domain) of a distributed StencilMatrix. - - pads_in : 1d int array - Paddings of the input space (domain) of a distributed StencilMatrix. - - starts_out : 1d int array - Starting indices of the output space (codomain) of a distributed StencilMatrix. - - ends_out : 1d int array - Ending indices of the output space (codomain) of a distributed StencilMatrix. - - pads_out : 1d int array - Paddings of the output space (codomain) of a distributed StencilMatrix. - - fun_q : 3d float array - The function evaluated at the points (nq_i*ii + iq, nq_j*jj + jq, nq_k*kk + kq), where iq a local quadrature point of interval ii. - - wts1 : 2d float array - Quadrature weights in direction eta1 in format (ii, iq). - - wts2 : 2d float array - Quadrature weights in direction eta2 in format (jj, jq). - - wts3 : 2d float array - Quadrature weights in direction eta3 in format (kk, kq). - - span1 : 2d int array - Knot span indices in direction eta1 in format (ii, iq). - - span2 : 2d int array - Knot span indices in direction eta2 in format (jj, jq). - - span3 : 2d int array - Knot span indices in direction eta3 in format (kk, kq). - - basis1 : 3d float array - Values of p1 + 1 non-zero eta-1 basis functions at quadrature points in format (ii, iq, basis function). - - basis2 : 3d float array - Values of p2 + 1 non-zero eta-2 basis functions at quadrature points in format (jj, jq, basis function). - - basis3 : 3d float array - Values of p3 + 1 non-zero eta-3 basis functions at quadrature points in format (kk, kq, basis function). - - sub1 : 1d int array - Sub-interval indices in direction 1. - - sub2 : 1d int array - Sub-interval indices in direction 2. - - sub3 : 1d int array - Sub-interval indices in direction 3. - - dim1_in : int - Dimension of the first direction of the input space - - dim2_in : int - Dimension of the second direction of the input space - - dim3_in : int - Dimension of the third direction of the input space - - p1_out : int - Spline degree of the first direction of the output space - - p2_out : int - Spline degree of the second direction of the output space - - p3_out : int - Spline degree of the third direction of the output space - """ - - # Start/end indices and paddings for distributed stencil matrix of input space - # si1 = starts_in[0] - # si2 = starts_in[1] - # si3 = starts_in[2] - # ei1 = ends_in[0] - # ei2 = ends_in[1] - # ei3 = ends_in[2] - pi1 = pads_in[0] - pi2 = pads_in[1] - pi3 = pads_in[2] - - # Start/end indices for distributed stencil matrix of output space - so1 = starts_out[0] - so2 = starts_out[1] - so3 = starts_out[2] - # eo1 = ends_out[0] - # eo2 = ends_out[1] - # eo3 = ends_out[2] - po1 = pads_out[0] - po2 = pads_out[1] - po3 = pads_out[2] - - # Spline degrees of input space - p1 = basis1.shape[2] - 1 - p2 = basis2.shape[2] - 1 - p3 = basis3.shape[2] - 1 - - # number of quadrature points - nq1 = span1.shape[1] - nq2 = span2.shape[1] - nq3 = span3.shape[1] - - # Set output to zero - mat[:] = 0.0 - - # Find column index for _data: - if dim1_out <= dim1_in: - cut1 = p1 - else: - cut1 = p1_out - - if dim2_out <= dim2_in: - cut2 = p2 - else: - cut2 = p2_out - - if dim3_out <= dim3_in: - cut3 = p3 - else: - cut3 = p3_out - - # Interval (either element or sub-interval thereof) - # ------------------------------------------------- - cumsub_i = 0 # Cumulative sub-interval index - for ii in range(span1.shape[0]): - cumsub_i += sub1[ii] - i = ii - cumsub_i # local DOF index - - cumsub_j = 0 # Cumulative sub-interval index - for jj in range(span2.shape[0]): - cumsub_j += sub2[jj] - j = jj - cumsub_j # local DOF index - - cumsub_k = 0 # Cumulative sub-interval index - for kk in range(span3.shape[0]): - cumsub_k += sub3[kk] - k = kk - cumsub_k # local DOF index - - # Quadrature point index in interval - # ---------------------------------- - for iq in range(nq1): - for jq in range(nq2): - for kq in range(nq3): - funval = ( - fun_q[ - nq1 * ii + iq, - nq2 * jj + jq, - nq3 * kk + kq, - ] - * wts1[ii, iq] - * wts2[jj, jq] - * wts3[kk, kq] - ) - - # Basis function of input space: - # ------------------------------ - for b1 in range(p1 + 1): - m = span1[ii, iq] - p1 + b1 # global index - # basis value - val1 = funval * basis1[ii, iq, b1] - - # Diff of global indices, needs to be adjusted for boundary conditions --> col1 - col1_tmp = m - (i + so1) - if col1_tmp > cut1: - m = m - dim1_in - elif col1_tmp < -cut1: - m = m + dim1_in - # add padding - col1 = pi1 + m - (i + so1) - - for b2 in range(p2 + 1): - # global index - n = span2[jj, jq] - p2 + b2 - val2 = val1 * basis2[jj, jq, b2] - - # Diff of global indices, needs to be adjusted for boundary conditions --> col2 - col2_tmp = n - (j + so2) - if col2_tmp > cut2: - n = n - dim2_in - elif col2_tmp < -cut2: - n = n + dim2_in - # add padding - col2 = pi2 + n - (j + so2) - - for b3 in range(p3 + 1): - # global index - o = span3[kk, kq] - p3 + b3 - value = val2 * basis3[kk, kq, b3] - - # Diff of global indices, needs to be adjusted for boundary conditions --> col3 - col3_tmp = o - (k + so3) - if col3_tmp > cut3: - o = o - dim3_in - elif col3_tmp < -cut3: - o = o + dim3_in - # add padding - col3 = pi3 + o - (k + so3) - # Row index: padding + local index. - mat[ - po1 + i, - po2 + j, - po3 + k, - col1, - col2, - col3, - ] += value diff --git a/src/struphy/feec/basis_projection_ops.py b/src/struphy/feec/basis_projection_ops.py index ba2c29dec..56ec67fbf 100644 --- a/src/struphy/feec/basis_projection_ops.py +++ b/src/struphy/feec/basis_projection_ops.py @@ -2,7 +2,6 @@ import cunumpy as xp import numpy as np -from cunumpy.kernels import PyccelKernel from feectools.api.settings import PSYDAC_BACKEND_GPYCCEL from feectools.fem.basic import FemSpace from feectools.fem.tensor import TensorFemSpace @@ -11,7 +10,9 @@ from feectools.linalg.stencil import StencilMatrix, StencilVector, StencilVectorSpace from maybempi import MPI -from struphy.feec import basis_projection_kernels +from struphy.feec.kernels.assemble_dofs_for_weighted_basisfuns_1d import assemble_dofs_for_weighted_basisfuns_1d +from struphy.feec.kernels.assemble_dofs_for_weighted_basisfuns_2d import assemble_dofs_for_weighted_basisfuns_2d +from struphy.feec.kernels.assemble_dofs_for_weighted_basisfuns_3d import assemble_dofs_for_weighted_basisfuns_3d from struphy.feec.linear_operators import BoundaryOperator from struphy.feec.local_projectors_kernels import assemble_basis_projection_operator_local from struphy.feec.projectors import CommutingProjector, CommutingProjectorLocal @@ -24,6 +25,13 @@ logger = logging.getLogger("struphy") +# assembly kernels by dimension of the problem (1d, 2d or 3d) +_ASSEMBLY_KERNELS = { + 1: assemble_dofs_for_weighted_basisfuns_1d, + 2: assemble_dofs_for_weighted_basisfuns_2d, + 3: assemble_dofs_for_weighted_basisfuns_3d, +} + class BasisProjectionOperators: r""" @@ -2079,12 +2087,7 @@ def assemble(self, weights=None): ) dofs_mat = self._dof_mat[i, j] - kernel = PyccelKernel( - getattr( - basis_projection_kernels, - "assemble_dofs_for_weighted_basisfuns_" + str(V.ldim) + "d", - ), - ) + kernel = _ASSEMBLY_KERNELS[V.ldim] logger.debug(f"Assemble block {i, j}") kernel( diff --git a/src/struphy/feec/boundary_mass.py b/src/struphy/feec/boundary_mass.py index 9c2953245..7d2d6daf1 100644 --- a/src/struphy/feec/boundary_mass.py +++ b/src/struphy/feec/boundary_mass.py @@ -2,13 +2,12 @@ from typing import Callable import cunumpy as xp -from cunumpy.kernels import PyccelKernel from feectools.api.settings import PSYDAC_BACKEND_GPYCCEL from feectools.linalg.basic import LinearOperator from feectools.linalg.block import BlockLinearOperator, BlockVector from feectools.linalg.stencil import StencilMatrix, StencilVector -from struphy.feec import mass_kernels +from struphy.feec.kernels.surface_kernel_3d_mat import surface_kernel_3d_mat from struphy.feec.mass import WeightedMassOperators from struphy.feec.psydac_derham import Derham, SplineFunction from struphy.geometry.base import Domain @@ -174,7 +173,7 @@ def __init__( self._setup_surface_data() # load assembly kernel - self._assembly_kernel = PyccelKernel(mass_kernels.surface_kernel_3d_mat) + self._assembly_kernel = surface_kernel_3d_mat self.assemble() diff --git a/src/struphy/feec/kernels/__init__.py b/src/struphy/feec/kernels/__init__.py index 31570db0e..c3fa8fa33 100644 --- a/src/struphy/feec/kernels/__init__.py +++ b/src/struphy/feec/kernels/__init__.py @@ -1,4 +1,4 @@ -"""FEEC utility kernels, one folder per kernel (see ``CUDA_STRATEGY.md``). +"""FEEC kernels, one folder per kernel (see ``CUDA_STRATEGY.md``). Each folder ```` holds the pyccel kernel ``_kernels.py`` and, once ported, its CUDA version ``_cuda.cu``. The folder's @@ -6,4 +6,12 @@ backend:: from struphy.feec.kernels.hybrid_weight import hybrid_weight + +Besides the utility kernel ``hybrid_weight``, the folders hold the assembly kernels of +:mod:`struphy.feec.mass` (``kernel_{1,2,3}d_mat``, ``kernel_{1,2,3}d_vec``, ``kernel_{1,2,3}d_eval``, +``kernel_3d_matrixfree``, ``kernel_3d_diag``), of :mod:`struphy.feec.boundary_mass` +(``surface_kernel_3d_mat``, ``surface_kernel_3d_vec``) and of :mod:`struphy.feec.basis_projection_ops` +(``assemble_dofs_for_weighted_basisfuns_{1,2,3}d``). None of the assembly kernels has a CUDA version yet; +they are declared with ``missing_cuda="fallback"``, so on the CuPy backend they run on host copies of the +arrays, as before. """ diff --git a/src/struphy/feec/kernels/assemble_dofs_for_weighted_basisfuns_1d/__init__.py b/src/struphy/feec/kernels/assemble_dofs_for_weighted_basisfuns_1d/__init__.py new file mode 100644 index 000000000..aa277960d --- /dev/null +++ b/src/struphy/feec/kernels/assemble_dofs_for_weighted_basisfuns_1d/__init__.py @@ -0,0 +1,11 @@ +"""The kernel of this folder: pyccel in ``_kernels.py``, CUDA in ``_cuda.cu`` once ported. + +No CUDA version yet: on the CuPy backend the pyccel kernel runs on host copies of the arrays (a warning is +emitted once), as before the move into this folder. +""" + +from cunumpy.kernels import Kernel + +from struphy.utils.cuda_arguments import CUDA_STRUCTS + +assemble_dofs_for_weighted_basisfuns_1d = Kernel.from_folder(__name__, structs=CUDA_STRUCTS, missing_cuda="fallback") diff --git a/src/struphy/feec/kernels/assemble_dofs_for_weighted_basisfuns_1d/assemble_dofs_for_weighted_basisfuns_1d_kernels.py b/src/struphy/feec/kernels/assemble_dofs_for_weighted_basisfuns_1d/assemble_dofs_for_weighted_basisfuns_1d_kernels.py new file mode 100644 index 000000000..abf7c22ee --- /dev/null +++ b/src/struphy/feec/kernels/assemble_dofs_for_weighted_basisfuns_1d/assemble_dofs_for_weighted_basisfuns_1d_kernels.py @@ -0,0 +1,128 @@ +"""Kernel assemble_dofs_for_weighted_basisfuns_1d: assembly kernel for basis projection operators.""" + + +def assemble_dofs_for_weighted_basisfuns_1d( + mat: "float[:,:]", + starts_in: "int[:]", + ends_in: "int[:]", + pads_in: "int[:]", + starts_out: "int[:]", + ends_out: "int[:]", + pads_out: "int[:]", + fun_q: "float[:]", + wts1: "float[:,:]", + span1: "int[:,:]", + basis1: "float[:,:,:]", + sub1: "int[:]", + dim1_in: int, + dim1_out: int, + p1_out: int, +): + """Kernel for assembling the matrix + + A_(i,j) = DOFS_i(fun*Lambda^in_j) , + + into the _data attribute of a StencilMatrix. + Here, DOFS_i are the degrees-of-freedom of the output space (codomain, must not be a product space), + Lambda^in_j are the basis functions of the input space (domain, must not be a product space), and fun is an arbitrary function. + + Parameters + ---------- + mat : 2d float array + _data attribute of StencilMatrix. + + starts_in : int + Starting index of the input space (domain) of a distributed StencilMatrix. + + ends_in : int + Ending index of the input space (domain) of a distributed StencilMatrix. + + pads_in : int + Paddings of the input space (domain) of a distributed StencilMatrix. + + starts_out : int + Starting indices of the output space (codomain) of a distributed StencilMatrix. + + ends_out : int + Ending indices of the output space (codomain) of a distributed StencilMatrix. + + pads_out : int + Paddings of the output space (codomain) of a distributed StencilMatrix. + + fun_q : 1d float array + The function evaluated at the points (nq*ii + iq), where iq a local quadrature point of interval ii. + + wts1 : 2d float array + Quadrature weights in format (ii, iq). + + span1 : 2d int array + Knot span indices in direction eta1 in format (ii, iq). + + basis1 : 3d float array + Values of p1 + 1 non-zero eta-1 basis functions at quadrature points in format (ii, iq, basis function). + + sub1 : 1d int array + Sub-interval indices in direction 1. + + dim1_in : int + Dimension of the first direction of the input space + + p1_out : int + Spline degree of the first direction of the output space + """ + + # Start/end indices and paddings for distributed stencil matrix of input space + # si1 = starts_in[0} + # ei1 = ends_in[0] + pi1 = pads_in[0] + + # Start/end indices for distributed stencil matrix of output space + so1 = starts_out[0] + # eo1 = ends_out[0] + po1 = pads_out[0] + + # Spline degrees of input space + p1 = basis1.shape[2] - 1 + + # number of quadrature points + nq1 = span1.shape[1] + + # Set output to zero + mat[:] = 0.0 + + # Find column index for _data: + if dim1_out <= dim1_in: + cut1 = p1 + else: + cut1 = p1_out + + # Interval (either element or sub-interval thereof) + # ------------------------------------------------- + cumsub_i = 0 # Cumulative sub-interval index + for ii in range(span1.shape[0]): + cumsub_i += sub1[ii] + i = ii - cumsub_i # local DOF index + + # Quadrature point index in interval + # ---------------------------------- + for iq in range(nq1): + funval = fun_q[nq1 * ii + iq] * wts1[ii, iq] + + # Basis function of input space: + # ------------------------------ + for b1 in range(p1 + 1): + m = span1[ii, iq] - p1 + b1 # global index + # basis value + value = funval * basis1[ii, iq, b1] + + # Diff of global indices, needs to be adjusted for boundary conditions --> col1 + col1_tmp = m - (i + so1) + if col1_tmp > cut1: + m = m - dim1_in + elif col1_tmp < -cut1: + m = m + dim1_in + # add padding + col1 = pi1 + m - (i + so1) + + # Row index: padding + local index. + mat[po1 + i, col1] += value diff --git a/src/struphy/feec/kernels/assemble_dofs_for_weighted_basisfuns_2d/__init__.py b/src/struphy/feec/kernels/assemble_dofs_for_weighted_basisfuns_2d/__init__.py new file mode 100644 index 000000000..34ccbda42 --- /dev/null +++ b/src/struphy/feec/kernels/assemble_dofs_for_weighted_basisfuns_2d/__init__.py @@ -0,0 +1,11 @@ +"""The kernel of this folder: pyccel in ``_kernels.py``, CUDA in ``_cuda.cu`` once ported. + +No CUDA version yet: on the CuPy backend the pyccel kernel runs on host copies of the arrays (a warning is +emitted once), as before the move into this folder. +""" + +from cunumpy.kernels import Kernel + +from struphy.utils.cuda_arguments import CUDA_STRUCTS + +assemble_dofs_for_weighted_basisfuns_2d = Kernel.from_folder(__name__, structs=CUDA_STRUCTS, missing_cuda="fallback") diff --git a/src/struphy/feec/kernels/assemble_dofs_for_weighted_basisfuns_2d/assemble_dofs_for_weighted_basisfuns_2d_kernels.py b/src/struphy/feec/kernels/assemble_dofs_for_weighted_basisfuns_2d/assemble_dofs_for_weighted_basisfuns_2d_kernels.py new file mode 100644 index 000000000..61db5cb4a --- /dev/null +++ b/src/struphy/feec/kernels/assemble_dofs_for_weighted_basisfuns_2d/assemble_dofs_for_weighted_basisfuns_2d_kernels.py @@ -0,0 +1,186 @@ +"""Kernel assemble_dofs_for_weighted_basisfuns_2d: assembly kernel for basis projection operators.""" + + +def assemble_dofs_for_weighted_basisfuns_2d( + mat: "float[:,:,:,:]", + starts_in: "int[:]", + ends_in: "int[:]", + pads_in: "int[:]", + starts_out: "int[:]", + ends_out: "int[:]", + pads_out: "int[:]", + fun_q: "float[:,:]", + wts1: "float[:,:]", + wts2: "float[:,:]", + span1: "int[:,:]", + span2: "int[:,:]", + basis1: "float[:,:,:]", + basis2: "float[:,:,:]", + sub1: "int[:]", + sub2: "int[:]", + dim1_in: int, + dim2_in: int, + dim1_out: int, + dim2_out: int, + p1_out: int, + p2_out: int, +): + """Kernel for assembling the matrix + + A_(ij,kl) = DOFS_ij(fun*Lambda^in_kl) , + + into the _data attribute of a StencilMatrix. + Here, DOFS_ij are the degrees-of-freedom of the output space (codomain, must not be a product space), + Lambda^in_kl are the basis functions of the input space (domain, must not be a product space), and fun is an arbitrary function. + + Parameters + ---------- + mat : 4d float array + _data attribute of StencilMatrix. + + starts_in : 1d int array + Starting indices of the input space (domain) of a distributed StencilMatrix. + + ends_in : 1d int array + Ending indices of the input space (domain) of a distributed StencilMatrix. + + pads_in : 1d int array + Paddings of the input space (domain) of a distributed StencilMatrix. + + starts_out : 1d int array + Starting indices of the output space (codomain) of a distributed StencilMatrix. + + ends_out : 1d int array + Ending indices of the output space (codomain) of a distributed StencilMatrix. + + pads_out : 1d int array + Paddings of the output space (codomain) of a distributed StencilMatrix. + + fun_q : 2d float array + The function evaluated at the points (nq_i*ii + iq, nq_j*jj + jq), where iq a local quadrature point of interval ii. + + wts1 : 2d float array + Quadrature weights in direction eta1 in format (ii, iq). + + wts2 : 2d float array + Quadrature weights in direction eta2 in format (jj, jq). + + span1 : 2d int array + Knot span indices in direction eta1 in format (ii, iq). + + span2 : 2d int array + Knot span indices in direction eta2 in format (jj, jq). + + basis1 : 3d float array + Values of p1 + 1 non-zero eta-1 basis functions at quadrature points in format (ii, iq, basis function). + + basis2 : 3d float array + Values of p2 + 1 non-zero eta-2 basis functions at quadrature points in format (jj, jq, basis function). + + sub1 : 1d int array + Sub-interval indices in direction 1. + + sub2 : 1d int array + Sub-interval indices in direction 2. + + dim1_in : int + Dimension of the first direction of the input space + + dim2_in : int + Dimension of the second direction of the input space + + p1_out : int + Spline degree of the first direction of the output space + + p2_out : int + Spline degree of the second direction of the output space + """ + + # Start/end indices and paddings for distributed stencil matrix of input space + # si1 = starts_in[0] + # si2 = starts_in[1] + # ei1 = ends_in[0] + # ei2 = ends_in[1] + pi1 = pads_in[0] + pi2 = pads_in[1] + + # Start/end indices for distributed stencil matrix of output space + so1 = starts_out[0] + so2 = starts_out[1] + # eo1 = ends_out[0] + # eo2 = ends_out[1] + po1 = pads_out[0] + po2 = pads_out[1] + + # Spline degrees of input space + p1 = basis1.shape[2] - 1 + p2 = basis2.shape[2] - 1 + + # number of quadrature points + nq1 = span1.shape[1] + nq2 = span2.shape[1] + + # Set output to zero + mat[:] = 0.0 + + # Find column index for _data: + if dim1_out <= dim1_in: + cut1 = p1 + else: + cut1 = p1_out + + if dim2_out <= dim2_in: + cut2 = p2 + else: + cut2 = p2_out + + # Interval (either element or sub-interval thereof) + # ------------------------------------------------- + cumsub_i = 0 # Cumulative sub-interval index + for ii in range(span1.shape[0]): + cumsub_i += sub1[ii] + i = ii - cumsub_i # local DOF index + + cumsub_j = 0 # Cumulative sub-interval index + for jj in range(span2.shape[0]): + cumsub_j += sub2[jj] + j = jj - cumsub_j # local DOF index + + # Quadrature point index in interval + # ---------------------------------- + for iq in range(nq1): + for jq in range(nq2): + funval = fun_q[nq1 * ii + iq, nq2 * jj + jq] * wts1[ii, iq] * wts2[jj, jq] + + # Basis function of input space: + # ------------------------------ + for b1 in range(p1 + 1): + m = span1[ii, iq] - p1 + b1 # global index + # basis value + val1 = funval * basis1[ii, iq, b1] + + # Diff of global indices, needs to be adjusted for boundary conditions --> col1 + col1_tmp = m - (i + so1) + if col1_tmp > cut1: + m = m - dim1_in + elif col1_tmp < -cut1: + m = m + dim1_in + # add padding + col1 = pi1 + m - (i + so1) + + for b2 in range(p2 + 1): + # global index + n = span2[jj, jq] - p2 + b2 + value = val1 * basis2[jj, jq, b2] + + # Diff of global indices, needs to be adjusted for boundary conditions --> col2 + col2_tmp = n - (j + so2) + if col2_tmp > cut2: + n = n - dim2_in + elif col2_tmp < -cut2: + n = n + dim2_in + # add padding + col2 = pi2 + n - (j + so2) + + # Row index: padding + local index. + mat[po1 + i, po2 + j, col1, col2] += value diff --git a/src/struphy/feec/kernels/assemble_dofs_for_weighted_basisfuns_3d/__init__.py b/src/struphy/feec/kernels/assemble_dofs_for_weighted_basisfuns_3d/__init__.py new file mode 100644 index 000000000..f1c1f6984 --- /dev/null +++ b/src/struphy/feec/kernels/assemble_dofs_for_weighted_basisfuns_3d/__init__.py @@ -0,0 +1,11 @@ +"""The kernel of this folder: pyccel in ``_kernels.py``, CUDA in ``_cuda.cu`` once ported. + +No CUDA version yet: on the CuPy backend the pyccel kernel runs on host copies of the arrays (a warning is +emitted once), as before the move into this folder. +""" + +from cunumpy.kernels import Kernel + +from struphy.utils.cuda_arguments import CUDA_STRUCTS + +assemble_dofs_for_weighted_basisfuns_3d = Kernel.from_folder(__name__, structs=CUDA_STRUCTS, missing_cuda="fallback") diff --git a/src/struphy/feec/kernels/assemble_dofs_for_weighted_basisfuns_3d/assemble_dofs_for_weighted_basisfuns_3d_kernels.py b/src/struphy/feec/kernels/assemble_dofs_for_weighted_basisfuns_3d/assemble_dofs_for_weighted_basisfuns_3d_kernels.py new file mode 100644 index 000000000..fe9a50597 --- /dev/null +++ b/src/struphy/feec/kernels/assemble_dofs_for_weighted_basisfuns_3d/assemble_dofs_for_weighted_basisfuns_3d_kernels.py @@ -0,0 +1,259 @@ +"""Kernel assemble_dofs_for_weighted_basisfuns_3d: assembly kernel for basis projection operators.""" + + +def assemble_dofs_for_weighted_basisfuns_3d( + mat: "float[:,:,:,:,:,:]", + starts_in: "int[:]", + ends_in: "int[:]", + pads_in: "int[:]", + starts_out: "int[:]", + ends_out: "int[:]", + pads_out: "int[:]", + fun_q: "float[:,:,:]", + wts1: "float[:,:]", + wts2: "float[:,:]", + wts3: "float[:,:]", + span1: "int[:,:]", + span2: "int[:,:]", + span3: "int[:,:]", + basis1: "float[:,:,:]", + basis2: "float[:,:,:]", + basis3: "float[:,:,:]", + sub1: "int[:]", + sub2: "int[:]", + sub3: "int[:]", + dim1_in: int, + dim2_in: int, + dim3_in: int, + dim1_out: int, + dim2_out: int, + dim3_out: int, + p1_out: int, + p2_out: int, + p3_out: int, +): + """Kernel for assembling the matrix + + A_(ijk,mno) = DOFS_ijk(fun*Lambda^in_mno) , + + into the _data attribute of a StencilMatrix. + Here, DOFS_ijk are the degrees-of-freedom of the output space (codomain, must not be a product space), + Lambda^in_mno are the basis functions of the input space (domain, must not be a product space), and fun is an arbitrary function. + + Parameters + ---------- + mat : 6d float array + _data attribute of StencilMatrix. + + starts_in : 1d int array + Starting indices of the input space (domain) of a distributed StencilMatrix. + + ends_in : 1d int array + Ending indices of the input space (domain) of a distributed StencilMatrix. + + pads_in : 1d int array + Paddings of the input space (domain) of a distributed StencilMatrix. + + starts_out : 1d int array + Starting indices of the output space (codomain) of a distributed StencilMatrix. + + ends_out : 1d int array + Ending indices of the output space (codomain) of a distributed StencilMatrix. + + pads_out : 1d int array + Paddings of the output space (codomain) of a distributed StencilMatrix. + + fun_q : 3d float array + The function evaluated at the points (nq_i*ii + iq, nq_j*jj + jq, nq_k*kk + kq), where iq a local quadrature point of interval ii. + + wts1 : 2d float array + Quadrature weights in direction eta1 in format (ii, iq). + + wts2 : 2d float array + Quadrature weights in direction eta2 in format (jj, jq). + + wts3 : 2d float array + Quadrature weights in direction eta3 in format (kk, kq). + + span1 : 2d int array + Knot span indices in direction eta1 in format (ii, iq). + + span2 : 2d int array + Knot span indices in direction eta2 in format (jj, jq). + + span3 : 2d int array + Knot span indices in direction eta3 in format (kk, kq). + + basis1 : 3d float array + Values of p1 + 1 non-zero eta-1 basis functions at quadrature points in format (ii, iq, basis function). + + basis2 : 3d float array + Values of p2 + 1 non-zero eta-2 basis functions at quadrature points in format (jj, jq, basis function). + + basis3 : 3d float array + Values of p3 + 1 non-zero eta-3 basis functions at quadrature points in format (kk, kq, basis function). + + sub1 : 1d int array + Sub-interval indices in direction 1. + + sub2 : 1d int array + Sub-interval indices in direction 2. + + sub3 : 1d int array + Sub-interval indices in direction 3. + + dim1_in : int + Dimension of the first direction of the input space + + dim2_in : int + Dimension of the second direction of the input space + + dim3_in : int + Dimension of the third direction of the input space + + p1_out : int + Spline degree of the first direction of the output space + + p2_out : int + Spline degree of the second direction of the output space + + p3_out : int + Spline degree of the third direction of the output space + """ + + # Start/end indices and paddings for distributed stencil matrix of input space + # si1 = starts_in[0] + # si2 = starts_in[1] + # si3 = starts_in[2] + # ei1 = ends_in[0] + # ei2 = ends_in[1] + # ei3 = ends_in[2] + pi1 = pads_in[0] + pi2 = pads_in[1] + pi3 = pads_in[2] + + # Start/end indices for distributed stencil matrix of output space + so1 = starts_out[0] + so2 = starts_out[1] + so3 = starts_out[2] + # eo1 = ends_out[0] + # eo2 = ends_out[1] + # eo3 = ends_out[2] + po1 = pads_out[0] + po2 = pads_out[1] + po3 = pads_out[2] + + # Spline degrees of input space + p1 = basis1.shape[2] - 1 + p2 = basis2.shape[2] - 1 + p3 = basis3.shape[2] - 1 + + # number of quadrature points + nq1 = span1.shape[1] + nq2 = span2.shape[1] + nq3 = span3.shape[1] + + # Set output to zero + mat[:] = 0.0 + + # Find column index for _data: + if dim1_out <= dim1_in: + cut1 = p1 + else: + cut1 = p1_out + + if dim2_out <= dim2_in: + cut2 = p2 + else: + cut2 = p2_out + + if dim3_out <= dim3_in: + cut3 = p3 + else: + cut3 = p3_out + + # Interval (either element or sub-interval thereof) + # ------------------------------------------------- + cumsub_i = 0 # Cumulative sub-interval index + for ii in range(span1.shape[0]): + cumsub_i += sub1[ii] + i = ii - cumsub_i # local DOF index + + cumsub_j = 0 # Cumulative sub-interval index + for jj in range(span2.shape[0]): + cumsub_j += sub2[jj] + j = jj - cumsub_j # local DOF index + + cumsub_k = 0 # Cumulative sub-interval index + for kk in range(span3.shape[0]): + cumsub_k += sub3[kk] + k = kk - cumsub_k # local DOF index + + # Quadrature point index in interval + # ---------------------------------- + for iq in range(nq1): + for jq in range(nq2): + for kq in range(nq3): + funval = ( + fun_q[ + nq1 * ii + iq, + nq2 * jj + jq, + nq3 * kk + kq, + ] + * wts1[ii, iq] + * wts2[jj, jq] + * wts3[kk, kq] + ) + + # Basis function of input space: + # ------------------------------ + for b1 in range(p1 + 1): + m = span1[ii, iq] - p1 + b1 # global index + # basis value + val1 = funval * basis1[ii, iq, b1] + + # Diff of global indices, needs to be adjusted for boundary conditions --> col1 + col1_tmp = m - (i + so1) + if col1_tmp > cut1: + m = m - dim1_in + elif col1_tmp < -cut1: + m = m + dim1_in + # add padding + col1 = pi1 + m - (i + so1) + + for b2 in range(p2 + 1): + # global index + n = span2[jj, jq] - p2 + b2 + val2 = val1 * basis2[jj, jq, b2] + + # Diff of global indices, needs to be adjusted for boundary conditions --> col2 + col2_tmp = n - (j + so2) + if col2_tmp > cut2: + n = n - dim2_in + elif col2_tmp < -cut2: + n = n + dim2_in + # add padding + col2 = pi2 + n - (j + so2) + + for b3 in range(p3 + 1): + # global index + o = span3[kk, kq] - p3 + b3 + value = val2 * basis3[kk, kq, b3] + + # Diff of global indices, needs to be adjusted for boundary conditions --> col3 + col3_tmp = o - (k + so3) + if col3_tmp > cut3: + o = o - dim3_in + elif col3_tmp < -cut3: + o = o + dim3_in + # add padding + col3 = pi3 + o - (k + so3) + # Row index: padding + local index. + mat[ + po1 + i, + po2 + j, + po3 + k, + col1, + col2, + col3, + ] += value diff --git a/src/struphy/feec/kernels/kernel_1d_eval/__init__.py b/src/struphy/feec/kernels/kernel_1d_eval/__init__.py new file mode 100644 index 000000000..200b2335f --- /dev/null +++ b/src/struphy/feec/kernels/kernel_1d_eval/__init__.py @@ -0,0 +1,11 @@ +"""The kernel of this folder: pyccel in ``_kernels.py``, CUDA in ``_cuda.cu`` once ported. + +No CUDA version yet: on the CuPy backend the pyccel kernel runs on host copies of the arrays (a warning is +emitted once), as before the move into this folder. +""" + +from cunumpy.kernels import Kernel + +from struphy.utils.cuda_arguments import CUDA_STRUCTS + +kernel_1d_eval = Kernel.from_folder(__name__, structs=CUDA_STRUCTS, missing_cuda="fallback") diff --git a/src/struphy/feec/kernels/kernel_1d_eval/kernel_1d_eval_kernels.py b/src/struphy/feec/kernels/kernel_1d_eval/kernel_1d_eval_kernels.py new file mode 100644 index 000000000..f5b192122 --- /dev/null +++ b/src/struphy/feec/kernels/kernel_1d_eval/kernel_1d_eval_kernels.py @@ -0,0 +1,55 @@ +"""Kernel kernel_1d_eval: integral kernel for mass matrices and L2-projections.""" + +from numpy import shape + + +def kernel_1d_eval( + spans1: "int[:]", + pi1: int, + starts1: int, + pads1: int, + bi1: "float[:,:,:,:]", + coeffs_data: "float[:]", + values: "float[:]", +): + """ + Evaluates sum_i1 [ coeffs_i1 * Lambda_i1(quad_eta1) ] for all quadrature points on the calling process. + + The results are written into values. + + Parameters + ---------- + spans1 : array[int] + Array of span indices; the span is the index of the last non-vanishing spline on each grid element + (cell). The length of the returned array is the number of elements (cells). + pi1 : int + Degree of the basis functions. + starts1 : int + Starting index on the current rank. + pads1 : int + Padding (=spline degree) for ghost regions in coeffs_data. + bi1 : "float[:,:,:,:]" + Values of basis functions. The indexing is [global element, local basis function, derivative, quadrature point]. + coeffs_data : "float[:]" + _data array of StencilVector holding the spline coefficients of the function to be evaluated. + values : "float[:]" + Output array (flattened over elements and quadrature points) holding the evaluated function values; + it is set to zero at the start of the kernel, i.e. it is overwritten, not added to. + """ + + values[:] = 0.0 + + ne1 = spans1.size + + nq1 = shape(bi1)[3] + + for iel1 in range(ne1): + for il1 in range(pi1 + 1): + # global spline indices + i_global1 = spans1[iel1] - pi1 + il1 + + # local spline indices (- starts --> can be negative, will therefore be written to ghost regions) + i_local1 = i_global1 - starts1 + + for q1 in range(nq1): + values[iel1 * nq1 + q1] += coeffs_data[pads1 + i_local1] * bi1[iel1, il1, 0, q1] diff --git a/src/struphy/feec/kernels/kernel_1d_mat/__init__.py b/src/struphy/feec/kernels/kernel_1d_mat/__init__.py new file mode 100644 index 000000000..0ba760166 --- /dev/null +++ b/src/struphy/feec/kernels/kernel_1d_mat/__init__.py @@ -0,0 +1,11 @@ +"""The kernel of this folder: pyccel in ``_kernels.py``, CUDA in ``_cuda.cu`` once ported. + +No CUDA version yet: on the CuPy backend the pyccel kernel runs on host copies of the arrays (a warning is +emitted once), as before the move into this folder. +""" + +from cunumpy.kernels import Kernel + +from struphy.utils.cuda_arguments import CUDA_STRUCTS + +kernel_1d_mat = Kernel.from_folder(__name__, structs=CUDA_STRUCTS, missing_cuda="fallback") diff --git a/src/struphy/feec/kernels/kernel_1d_mat/kernel_1d_mat_kernels.py b/src/struphy/feec/kernels/kernel_1d_mat/kernel_1d_mat_kernels.py new file mode 100644 index 000000000..8a60d2367 --- /dev/null +++ b/src/struphy/feec/kernels/kernel_1d_mat/kernel_1d_mat_kernels.py @@ -0,0 +1,68 @@ +"""Kernel kernel_1d_mat: integral kernel for mass matrices and L2-projections.""" + +from numpy import shape + + +def kernel_1d_mat( + spans1: "int[:]", + pi1: int, + pj1: int, + starts1: int, + pads1: int, + w1: "float[:,:]", + bi1: "float[:,:,:,:]", + bj1: "float[:,:,:,:]", + mat_fun: "float[:]", + data: "float[:,:]", +): + """ + Performs the integration of Lambda_(i1) * mat_fun(eta1) * Lambda_(j1) for the basis functions (i1, j1) available on the calling process. + + The results are written into data (attention: data is NOT set to zero first, but the results are added to data). + + Parameters + ---------- + spans1 : array[int] + Array of span indices; the span is the index of the last non-vanishing spline on each grid element + (cell). The length of the returned array is the number of elements (cells). + pi1 : int + Degree of the codomain basis functions. + pj1 : int + Degree of the domain basis functions. + starts1 : int + Starting index on the current rank. + pads1 : int + Padding (=spline degree) for ghost regions in data. + w1 : "float[:,:]" + Quadrature weights. The indexing is [global element, quadrature point]. + bi1 : "float[:,:,:,:]" + Values of codomain basis functions. The indexing is [global element, local basis function, derivative, quadrature point]. + bj1 : "float[:,:,:,:]" + Values of domain basis functions. The indexing is [global element, local basis function, derivative, quadrature point]. + mat_fun : "float[:]" + Function under the integral evaluated at quadrature points (flattened). + data : "float[:,:]" + _data array of StencilMatrix to store the results. + """ + + # number of elements + ne1 = spans1.size + + # number of quadrature points in each element + nq1 = shape(w1)[1] + + for iel1 in range(ne1): + for il1 in range(pi1 + 1): + # global spline indices + i_global1 = spans1[iel1] - pi1 + il1 + + # local spline indices (- starts --> can be negative, will therefore be written to ghost regions) + i_local1 = i_global1 - starts1 + + for jl1 in range(pj1 + 1): + value = 0.0 + + for q1 in range(nq1): + value += w1[iel1, q1] * bi1[iel1, il1, 0, q1] * bj1[iel1, jl1, 0, q1] * mat_fun[iel1 * nq1 + q1] + + data[pads1 + i_local1, pads1 + jl1 - il1] += value diff --git a/src/struphy/feec/kernels/kernel_1d_vec/__init__.py b/src/struphy/feec/kernels/kernel_1d_vec/__init__.py new file mode 100644 index 000000000..1229129d5 --- /dev/null +++ b/src/struphy/feec/kernels/kernel_1d_vec/__init__.py @@ -0,0 +1,11 @@ +"""The kernel of this folder: pyccel in ``_kernels.py``, CUDA in ``_cuda.cu`` once ported. + +No CUDA version yet: on the CuPy backend the pyccel kernel runs on host copies of the arrays (a warning is +emitted once), as before the move into this folder. +""" + +from cunumpy.kernels import Kernel + +from struphy.utils.cuda_arguments import CUDA_STRUCTS + +kernel_1d_vec = Kernel.from_folder(__name__, structs=CUDA_STRUCTS, missing_cuda="fallback") diff --git a/src/struphy/feec/kernels/kernel_1d_vec/kernel_1d_vec_kernels.py b/src/struphy/feec/kernels/kernel_1d_vec/kernel_1d_vec_kernels.py new file mode 100644 index 000000000..6d80200eb --- /dev/null +++ b/src/struphy/feec/kernels/kernel_1d_vec/kernel_1d_vec_kernels.py @@ -0,0 +1,59 @@ +"""Kernel kernel_1d_vec: integral kernel for mass matrices and L2-projections.""" + +from numpy import shape + + +def kernel_1d_vec( + spans1: "int[:]", + pi1: int, + starts1: int, + pads1: int, + w1: "float[:,:]", + bi1: "float[:,:,:,:]", + mat_fun: "float[:]", + data: "float[:]", +): + """ + Performs the integration of Lambda_(i1) * mat_fun(eta1) for the basis functions (i1) available on the calling process. + + The results are written into data (attention: data is NOT set to zero first, but the results are added to data). + + Parameters + ---------- + spans1 : array[int] + Array of span indices; the span is the index of the last non-vanishing spline on each grid element + (cell). The length of the returned array is the number of elements (cells). + pi1 : int + Degree of the basis functions. + starts1 : int + Starting index on the current rank. + pads1 : int + Padding (=spline degree) for ghost regions in data. + w1 : "float[:,:]" + Quadrature weights. The indexing is [global element, quadrature point]. + bi1 : "float[:,:,:,:]" + Values of basis functions. The indexing is [global element, local basis function, derivative, quadrature point]. + mat_fun : "float[:]" + Function under the integral evaluated at quadrature points (flattened). + data : "float[:]" + _data array of StencilVector to store the results. + """ + + ne1 = spans1.size + + nq1 = shape(w1)[1] + + for iel1 in range(ne1): + for il1 in range(pi1 + 1): + # global spline indices + i_global1 = spans1[iel1] - pi1 + il1 + + # local spline indices (- starts --> can be negative, will therefore be written to ghost regions) + i_local1 = i_global1 - starts1 + + value = 0.0 + + for q1 in range(nq1): + value += w1[iel1, q1] * bi1[iel1, il1, 0, q1] * mat_fun[iel1 * nq1 + q1] + + data[pads1 + i_local1] += value diff --git a/src/struphy/feec/kernels/kernel_2d_eval/__init__.py b/src/struphy/feec/kernels/kernel_2d_eval/__init__.py new file mode 100644 index 000000000..57e1d745a --- /dev/null +++ b/src/struphy/feec/kernels/kernel_2d_eval/__init__.py @@ -0,0 +1,11 @@ +"""The kernel of this folder: pyccel in ``_kernels.py``, CUDA in ``_cuda.cu`` once ported. + +No CUDA version yet: on the CuPy backend the pyccel kernel runs on host copies of the arrays (a warning is +emitted once), as before the move into this folder. +""" + +from cunumpy.kernels import Kernel + +from struphy.utils.cuda_arguments import CUDA_STRUCTS + +kernel_2d_eval = Kernel.from_folder(__name__, structs=CUDA_STRUCTS, missing_cuda="fallback") diff --git a/src/struphy/feec/kernels/kernel_2d_eval/kernel_2d_eval_kernels.py b/src/struphy/feec/kernels/kernel_2d_eval/kernel_2d_eval_kernels.py new file mode 100644 index 000000000..9a379cc62 --- /dev/null +++ b/src/struphy/feec/kernels/kernel_2d_eval/kernel_2d_eval_kernels.py @@ -0,0 +1,72 @@ +"""Kernel kernel_2d_eval: integral kernel for mass matrices and L2-projections.""" + +from numpy import shape + + +def kernel_2d_eval( + spans1: "int[:]", + spans2: "int[:]", + pi1: int, + pi2: int, + starts1: int, + starts2: int, + pads1: int, + pads2: int, + bi1: "float[:,:,:,:]", + bi2: "float[:,:,:,:]", + coeffs_data: "float[:,:]", + values: "float[:,:]", +): + """ + Evaluates sum_(i1, i2) [ coeffs_{i1,i2} * Lambda_{i1, i2}(quad_eta1, quad_eta2) ] for all quadrature points on the calling process. + + The results are written into values. + + Parameters + ---------- + spans1, spans2 : array[int] + Arrays of span indices in direction 1 and 2; the span is the index of the last non-vanishing spline + on each grid element (cell). The length of each array is the number of elements (cells) in that direction. + pi1, pi2 : int + Degree of the basis functions in direction 1 and 2. + starts1, starts2 : int + Starting index on the current rank, in direction 1 and 2. + pads1, pads2 : int + Padding (=spline degree) for ghost regions in coeffs_data, in direction 1 and 2. + bi1, bi2 : "float[:,:,:,:]" + Values of basis functions in direction 1 and 2. The indexing is + [global element, local basis function, derivative, quadrature point]. + coeffs_data : "float[:,:]" + _data array of StencilVector holding the spline coefficients of the function to be evaluated. + values : "float[:,:]" + Output array (flattened over elements and quadrature points in each direction) holding the evaluated + function values; it is set to zero at the start of the kernel, i.e. it is overwritten, not added to. + """ + + values[:, :] = 0.0 + + ne1 = spans1.size + ne2 = spans2.size + + nq1 = shape(bi1)[3] + nq2 = shape(bi2)[3] + + for iel1 in range(ne1): + for iel2 in range(ne2): + for il1 in range(pi1 + 1): + for il2 in range(pi2 + 1): + # global spline indices + i_global1 = spans1[iel1] - pi1 + il1 + i_global2 = spans2[iel2] - pi2 + il2 + + # local spline indices (- starts --> can be negative, will therefore be written to ghost regions) + i_local1 = i_global1 - starts1 + i_local2 = i_global2 - starts2 + + for q1 in range(nq1): + for q2 in range(nq2): + values[iel1 * nq1 + q1, iel2 * nq2 + q2] += ( + coeffs_data[pads1 + i_local1, pads2 + i_local2] + * bi1[iel1, il1, 0, q1] + * bi2[iel2, il2, 0, q2] + ) diff --git a/src/struphy/feec/kernels/kernel_2d_mat/__init__.py b/src/struphy/feec/kernels/kernel_2d_mat/__init__.py new file mode 100644 index 000000000..006c8b25d --- /dev/null +++ b/src/struphy/feec/kernels/kernel_2d_mat/__init__.py @@ -0,0 +1,11 @@ +"""The kernel of this folder: pyccel in ``_kernels.py``, CUDA in ``_cuda.cu`` once ported. + +No CUDA version yet: on the CuPy backend the pyccel kernel runs on host copies of the arrays (a warning is +emitted once), as before the move into this folder. +""" + +from cunumpy.kernels import Kernel + +from struphy.utils.cuda_arguments import CUDA_STRUCTS + +kernel_2d_mat = Kernel.from_folder(__name__, structs=CUDA_STRUCTS, missing_cuda="fallback") diff --git a/src/struphy/feec/kernels/kernel_2d_mat/kernel_2d_mat_kernels.py b/src/struphy/feec/kernels/kernel_2d_mat/kernel_2d_mat_kernels.py new file mode 100644 index 000000000..559260cba --- /dev/null +++ b/src/struphy/feec/kernels/kernel_2d_mat/kernel_2d_mat_kernels.py @@ -0,0 +1,88 @@ +"""Kernel kernel_2d_mat: integral kernel for mass matrices and L2-projections.""" + +from numpy import shape + + +def kernel_2d_mat( + spans1: "int[:]", + spans2: "int[:]", + pi1: int, + pi2: int, + pj1: int, + pj2: int, + starts1: int, + starts2: int, + pads1: int, + pads2: int, + w1: "float[:,:]", + w2: "float[:,:]", + bi1: "float[:,:,:,:]", + bi2: "float[:,:,:,:]", + bj1: "float[:,:,:,:]", + bj2: "float[:,:,:,:]", + mat_fun: "float[:,:]", + data: "float[:,:,:,:]", +): + """ + Performs the integration of Lambda_(i1, i2) * mat_fun(eta1, eta2) * Lambda_(j1, j2) for the basis functions (i1, i2, j1, j2) available on the calling process. + + The results are written into data (attention: data is NOT set to zero first, but the results are added to data). + + Parameters + ---------- + spans1, spans2 : array[int] + Arrays of span indices in direction 1 and 2; the span is the index of the last non-vanishing spline + on each grid element (cell). The length of each array is the number of elements (cells) in that direction. + pi1, pi2 : int + Degree of the codomain basis functions in direction 1 and 2. + pj1, pj2 : int + Degree of the domain basis functions in direction 1 and 2. + starts1, starts2 : int + Starting index on the current rank, in direction 1 and 2. + pads1, pads2 : int + Padding (=spline degree) for ghost regions in data, in direction 1 and 2. + w1, w2 : "float[:,:]" + Quadrature weights in direction 1 and 2. The indexing is [global element, quadrature point]. + bi1, bi2 : "float[:,:,:,:]" + Values of codomain basis functions in direction 1 and 2. The indexing is + [global element, local basis function, derivative, quadrature point]. + bj1, bj2 : "float[:,:,:,:]" + Values of domain basis functions in direction 1 and 2, same indexing convention as bi1, bi2. + mat_fun : "float[:,:]" + Function under the integral evaluated at quadrature points (flattened in each direction). + The indexing is [flattened quadrature point in direction 1, flattened quadrature point in direction 2]. + data : "float[:,:,:,:]" + _data array of StencilMatrix to store the results. + """ + + ne1 = spans1.size + ne2 = spans2.size + + nq1 = shape(w1)[1] + nq2 = shape(w2)[1] + + for iel1 in range(ne1): + for iel2 in range(ne2): + for il1 in range(pi1 + 1): + for il2 in range(pi2 + 1): + # global spline indices + i_global1 = spans1[iel1] - pi1 + il1 + i_global2 = spans2[iel2] - pi2 + il2 + + # local spline indices (- starts --> can be negative, will therefore be written to ghost regions) + i_local1 = i_global1 - starts1 + i_local2 = i_global2 - starts2 + + for jl1 in range(pj1 + 1): + for jl2 in range(pj2 + 1): + value = 0.0 + + for q1 in range(nq1): + for q2 in range(nq2): + wvol = w1[iel1, q1] * w2[iel2, q2] * mat_fun[iel1 * nq1 + q1, iel2 * nq2 + q2] + bi = bi1[iel1, il1, 0, q1] * bi2[iel2, il2, 0, q2] + bj = bj1[iel1, jl1, 0, q1] * bj2[iel2, jl2, 0, q2] + + value += wvol * bi * bj + + data[pads1 + i_local1, pads2 + i_local2, pads1 + jl1 - il1, pads2 + jl2 - il2] += value diff --git a/src/struphy/feec/kernels/kernel_2d_vec/__init__.py b/src/struphy/feec/kernels/kernel_2d_vec/__init__.py new file mode 100644 index 000000000..e1b6f44c1 --- /dev/null +++ b/src/struphy/feec/kernels/kernel_2d_vec/__init__.py @@ -0,0 +1,11 @@ +"""The kernel of this folder: pyccel in ``_kernels.py``, CUDA in ``_cuda.cu`` once ported. + +No CUDA version yet: on the CuPy backend the pyccel kernel runs on host copies of the arrays (a warning is +emitted once), as before the move into this folder. +""" + +from cunumpy.kernels import Kernel + +from struphy.utils.cuda_arguments import CUDA_STRUCTS + +kernel_2d_vec = Kernel.from_folder(__name__, structs=CUDA_STRUCTS, missing_cuda="fallback") diff --git a/src/struphy/feec/kernels/kernel_2d_vec/kernel_2d_vec_kernels.py b/src/struphy/feec/kernels/kernel_2d_vec/kernel_2d_vec_kernels.py new file mode 100644 index 000000000..48ada9248 --- /dev/null +++ b/src/struphy/feec/kernels/kernel_2d_vec/kernel_2d_vec_kernels.py @@ -0,0 +1,76 @@ +"""Kernel kernel_2d_vec: integral kernel for mass matrices and L2-projections.""" + +from numpy import shape + + +def kernel_2d_vec( + spans1: "int[:]", + spans2: "int[:]", + pi1: int, + pi2: int, + starts1: int, + starts2: int, + pads1: int, + pads2: int, + w1: "float[:,:]", + w2: "float[:,:]", + bi1: "float[:,:,:,:]", + bi2: "float[:,:,:,:]", + mat_fun: "float[:,:]", + data: "float[:,:]", +): + """ + Performs the integration of Lambda_(i1, i2) * mat_fun(eta1, eta2) for the basis functions (i1, i2) available on the calling process. + + The results are written into data (attention: data is NOT set to zero first, but the results are added to data). + + Parameters + ---------- + spans1, spans2 : array[int] + Arrays of span indices in direction 1 and 2; the span is the index of the last non-vanishing spline + on each grid element (cell). The length of each array is the number of elements (cells) in that direction. + pi1, pi2 : int + Degree of the basis functions in direction 1 and 2. + starts1, starts2 : int + Starting index on the current rank, in direction 1 and 2. + pads1, pads2 : int + Padding (=spline degree) for ghost regions in data, in direction 1 and 2. + w1, w2 : "float[:,:]" + Quadrature weights in direction 1 and 2. The indexing is [global element, quadrature point]. + bi1, bi2 : "float[:,:,:,:]" + Values of basis functions in direction 1 and 2. The indexing is + [global element, local basis function, derivative, quadrature point]. + mat_fun : "float[:,:]" + Function under the integral evaluated at quadrature points (flattened in each direction). + The indexing is [flattened quadrature point in direction 1, flattened quadrature point in direction 2]. + data : "float[:,:]" + _data array of StencilVector to store the results. + """ + + ne1 = spans1.size + ne2 = spans2.size + + nq1 = shape(w1)[1] + nq2 = shape(w2)[1] + + for iel1 in range(ne1): + for iel2 in range(ne2): + for il1 in range(pi1 + 1): + for il2 in range(pi2 + 1): + # global spline indices + i_global1 = spans1[iel1] - pi1 + il1 + i_global2 = spans2[iel2] - pi2 + il2 + + # local spline indices (- starts --> can be negative, will therefore be written to ghost regions) + i_local1 = i_global1 - starts1 + i_local2 = i_global2 - starts2 + + value = 0.0 + + for q1 in range(nq1): + for q2 in range(nq2): + wvol = w1[iel1, q1] * w2[iel2, q2] * mat_fun[iel1 * nq1 + q1, iel2 * nq2 + q2] + + value += wvol * bi1[iel1, il1, 0, q1] * bi2[iel2, il2, 0, q2] + + data[pads1 + i_local1, pads2 + i_local2] += value diff --git a/src/struphy/feec/kernels/kernel_3d_diag/__init__.py b/src/struphy/feec/kernels/kernel_3d_diag/__init__.py new file mode 100644 index 000000000..430e38e20 --- /dev/null +++ b/src/struphy/feec/kernels/kernel_3d_diag/__init__.py @@ -0,0 +1,11 @@ +"""The kernel of this folder: pyccel in ``_kernels.py``, CUDA in ``_cuda.cu`` once ported. + +No CUDA version yet: on the CuPy backend the pyccel kernel runs on host copies of the arrays (a warning is +emitted once), as before the move into this folder. +""" + +from cunumpy.kernels import Kernel + +from struphy.utils.cuda_arguments import CUDA_STRUCTS + +kernel_3d_diag = Kernel.from_folder(__name__, structs=CUDA_STRUCTS, missing_cuda="fallback") diff --git a/src/struphy/feec/kernels/kernel_3d_diag/kernel_3d_diag_kernels.py b/src/struphy/feec/kernels/kernel_3d_diag/kernel_3d_diag_kernels.py new file mode 100644 index 000000000..937af9040 --- /dev/null +++ b/src/struphy/feec/kernels/kernel_3d_diag/kernel_3d_diag_kernels.py @@ -0,0 +1,118 @@ +"""Kernel kernel_3d_diag: integral kernel for mass matrices and L2-projections.""" + +import numpy as np +from numpy import shape + + +def kernel_3d_diag( + spans1: "int[:]", + spans2: "int[:]", + spans3: "int[:]", + pi1: int, + pi2: int, + pi3: int, + starts1: int, + starts2: int, + starts3: int, + pads1: int, + pads2: int, + pads3: int, + w1: "float[:,:]", + w2: "float[:,:]", + w3: "float[:,:]", + bi1: "float[:,:,:,:]", + bi2: "float[:,:,:,:]", + bi3: "float[:,:,:,:]", + mat_fun: "float[:,:,:]", + data: "float[:,:,:]", +): + """ + Computes the diagonal of a mass matrix, assuming that the domain and the codomain are the same. + + The results are written into data (attention: data is NOT set to zero first, but the results are added to data). + + Parameters + ---------- + spans1, spans2, spans3 : array[int] + Arrays of span indices in direction 1, 2 and 3; the span is the index of the last non-vanishing spline + on each grid element (cell). The length of each array is the number of elements (cells) in that direction. + pi1, pi2, pi3 : int + Degree of the basis functions in direction 1, 2 and 3. + starts1, starts2, starts3 : int + Starting index on the current rank, in direction 1, 2 and 3. + pads1, pads2, pads3 : int + Padding (=spline degree) for ghost regions, in direction 1, 2 and 3. + w1, w2, w3 : "float[:,:]" + Quadrature weights in direction 1, 2 and 3. The indexing is [global element, quadrature point]. + bi1, bi2, bi3 : "float[:,:,:,:]" + Values of basis functions in direction 1, 2 and 3. The indexing is + [global element, local basis function, derivative, quadrature point]. + mat_fun : "float[:,:,:]" + Function under the integral evaluated at quadrature points (flattened in each direction). + The indexing is [flattened quad. point dir. 1, flattened quad. point dir. 2, flattened quad. point dir. 3]. + data : "float[:,:,:]" + _data array of a StencilVector (with ghost regions) to store the results. Contributions to entries owned + by other processes (or periodic images) are written to the ghost regions and must be added to their + owners afterwards with exchange_assembly_data(). + """ + + ne1 = spans1.size + ne2 = spans2.size + ne3 = spans3.size + + nq1 = shape(w1)[1] + nq2 = shape(w2)[1] + nq3 = shape(w3)[1] + + tmp_bi1 = np.zeros(nq1) + tmp_bi2 = np.zeros(nq2) + tmp_bi3 = np.zeros(nq3) + + tmp_w1 = np.zeros(nq1) + tmp_w2 = np.zeros(nq2) + tmp_w3 = np.zeros(nq3) + + tmp_mat_fun = np.zeros((nq1, nq2, nq3)) + + for iel1 in range(ne1): + for iel2 in range(ne2): + for iel3 in range(ne3): + tmp_mat_fun[:, :, :] = mat_fun[ + iel1 * nq1 : (iel1 + 1) * nq1, + iel2 * nq2 : (iel2 + 1) * nq2, + iel3 * nq3 : (iel3 + 1) * nq3, + ] + + tmp_w1[:] = w1[iel1, :] + tmp_w2[:] = w2[iel2, :] + tmp_w3[:] = w3[iel3, :] + + for il1 in range(pi1 + 1): + for il2 in range(pi2 + 1): + for il3 in range(pi3 + 1): + tmp_bi1[:] = bi1[iel1, il1, 0, :] + tmp_bi2[:] = bi2[iel2, il2, 0, :] + tmp_bi3[:] = bi3[iel3, il3, 0, :] + + # global spline indices + i_global1 = spans1[iel1] - pi1 + il1 + i_global2 = spans2[iel2] - pi2 + il2 + i_global3 = spans3[iel3] - pi3 + il3 + + # local spline indices (- starts --> can be negative, will therefore be written to ghost regions) + i_local1 = i_global1 - starts1 + pads1 + i_local2 = i_global2 - starts2 + pads2 + i_local3 = i_global3 - starts3 + pads3 + + value = 0.0 + + for q1 in range(nq1): + for q2 in range(nq2): + for q3 in range(nq3): + wvol = tmp_w1[q1] * tmp_w2[q2] * tmp_w3[q3] * tmp_mat_fun[q1, q2, q3] + + bi = tmp_bi1[q1] * tmp_bi2[q2] * tmp_bi3[q3] + + value += wvol * bi * bi + + data[i_local1, i_local2, i_local3] += value diff --git a/src/struphy/feec/kernels/kernel_3d_eval/__init__.py b/src/struphy/feec/kernels/kernel_3d_eval/__init__.py new file mode 100644 index 000000000..03c74e294 --- /dev/null +++ b/src/struphy/feec/kernels/kernel_3d_eval/__init__.py @@ -0,0 +1,11 @@ +"""The kernel of this folder: pyccel in ``_kernels.py``, CUDA in ``_cuda.cu`` once ported. + +No CUDA version yet: on the CuPy backend the pyccel kernel runs on host copies of the arrays (a warning is +emitted once), as before the move into this folder. +""" + +from cunumpy.kernels import Kernel + +from struphy.utils.cuda_arguments import CUDA_STRUCTS + +kernel_3d_eval = Kernel.from_folder(__name__, structs=CUDA_STRUCTS, missing_cuda="fallback") diff --git a/src/struphy/feec/kernels/kernel_3d_eval/kernel_3d_eval_kernels.py b/src/struphy/feec/kernels/kernel_3d_eval/kernel_3d_eval_kernels.py new file mode 100644 index 000000000..86afed5af --- /dev/null +++ b/src/struphy/feec/kernels/kernel_3d_eval/kernel_3d_eval_kernels.py @@ -0,0 +1,85 @@ +"""Kernel kernel_3d_eval: integral kernel for mass matrices and L2-projections.""" + +from numpy import shape + + +def kernel_3d_eval( + spans1: "int[:]", + spans2: "int[:]", + spans3: "int[:]", + pi1: int, + pi2: int, + pi3: int, + starts1: int, + starts2: int, + starts3: int, + pads1: int, + pads2: int, + pads3: int, + bi1: "float[:,:,:,:]", + bi2: "float[:,:,:,:]", + bi3: "float[:,:,:,:]", + coeffs_data: "float[:,:,:]", + values: "float[:,:,:]", +): + """ + Evaluates sum_(i1,i2,i3) [ coeffs_{i1,i2,i3} * Lambda_{i1,i2,i3}(quad_eta1, quad_eta2, quad_eta3) ] for all quadrature points on the calling process. + + The results are written into values. + + Parameters + ---------- + spans1, spans2, spans3 : array[int] + Arrays of span indices in direction 1, 2 and 3; the span is the index of the last non-vanishing spline + on each grid element (cell). The length of each array is the number of elements (cells) in that direction. + pi1, pi2, pi3 : int + Degree of the basis functions in direction 1, 2 and 3. + starts1, starts2, starts3 : int + Starting index on the current rank, in direction 1, 2 and 3. + pads1, pads2, pads3 : int + Padding (=spline degree) for ghost regions in coeffs_data, in direction 1, 2 and 3. + bi1, bi2, bi3 : "float[:,:,:,:]" + Values of basis functions in direction 1, 2 and 3. The indexing is + [global element, local basis function, derivative, quadrature point]. + coeffs_data : "float[:,:,:]" + _data array of StencilVector holding the spline coefficients of the function to be evaluated. + values : "float[:,:,:]" + Output array (flattened over elements and quadrature points in each direction) holding the evaluated + function values; it is set to zero at the start of the kernel, i.e. it is overwritten, not added to. + """ + + values[:, :, :] = 0.0 + + ne1 = spans1.size + ne2 = spans2.size + ne3 = spans3.size + + nq1 = shape(bi1)[3] + nq2 = shape(bi2)[3] + nq3 = shape(bi3)[3] + + for iel1 in range(ne1): + for iel2 in range(ne2): + for iel3 in range(ne3): + for il1 in range(pi1 + 1): + for il2 in range(pi2 + 1): + for il3 in range(pi3 + 1): + # global spline indices + i_global1 = spans1[iel1] - pi1 + il1 + i_global2 = spans2[iel2] - pi2 + il2 + i_global3 = spans3[iel3] - pi3 + il3 + + # local spline indices (- starts --> can be negative, will therefore be written to ghost regions) + i_local1 = i_global1 - starts1 + i_local2 = i_global2 - starts2 + i_local3 = i_global3 - starts3 + + for q1 in range(nq1): + for q2 in range(nq2): + for q3 in range(nq3): + values[iel1 * nq1 + q1, iel2 * nq2 + q2, iel3 * nq3 + q3] += ( + coeffs_data[pads1 + i_local1, pads2 + i_local2, pads3 + i_local3] + * bi1[iel1, il1, 0, q1] + * bi2[iel2, il2, 0, q2] + * bi3[iel3, il3, 0, q3] + ) diff --git a/src/struphy/feec/kernels/kernel_3d_mat/__init__.py b/src/struphy/feec/kernels/kernel_3d_mat/__init__.py new file mode 100644 index 000000000..2c87c0a26 --- /dev/null +++ b/src/struphy/feec/kernels/kernel_3d_mat/__init__.py @@ -0,0 +1,11 @@ +"""The kernel of this folder: pyccel in ``_kernels.py``, CUDA in ``_cuda.cu`` once ported. + +No CUDA version yet: on the CuPy backend the pyccel kernel runs on host copies of the arrays (a warning is +emitted once), as before the move into this folder. +""" + +from cunumpy.kernels import Kernel + +from struphy.utils.cuda_arguments import CUDA_STRUCTS + +kernel_3d_mat = Kernel.from_folder(__name__, structs=CUDA_STRUCTS, missing_cuda="fallback") diff --git a/src/struphy/feec/kernels/kernel_3d_mat/kernel_3d_mat_kernels.py b/src/struphy/feec/kernels/kernel_3d_mat/kernel_3d_mat_kernels.py new file mode 100644 index 000000000..26677b187 --- /dev/null +++ b/src/struphy/feec/kernels/kernel_3d_mat/kernel_3d_mat_kernels.py @@ -0,0 +1,147 @@ +"""Kernel kernel_3d_mat: integral kernel for mass matrices and L2-projections.""" + +import numpy as np +from numpy import shape + + +def kernel_3d_mat( + spans1: "int[:]", + spans2: "int[:]", + spans3: "int[:]", + pi1: int, + pi2: int, + pi3: int, + pj1: int, + pj2: int, + pj3: int, + starts1: int, + starts2: int, + starts3: int, + pads1: int, + pads2: int, + pads3: int, + w1: "float[:,:]", + w2: "float[:,:]", + w3: "float[:,:]", + bi1: "float[:,:,:,:]", + bi2: "float[:,:,:,:]", + bi3: "float[:,:,:,:]", + bj1: "float[:,:,:,:]", + bj2: "float[:,:,:,:]", + bj3: "float[:,:,:,:]", + mat_fun: "float[:,:,:]", + data: "float[:,:,:,:,:,:]", +): + """ + Performs the integration of Lambda_(i1,i2,i3) * mat_fun(eta1, eta2, eta3) * Lambda_(j1,j2,j3) for the basis functions (i1,i2,i3, j1,j2,j3) available on the calling process. + + The results are written into data (attention: data is NOT set to zero first, but the results are added to data). + + Parameters + ---------- + spans1, spans2, spans3 : array[int] + Arrays of span indices in direction 1, 2 and 3; the span is the index of the last non-vanishing spline + on each grid element (cell). The length of each array is the number of elements (cells) in that direction. + pi1, pi2, pi3 : int + Degree of the codomain basis functions in direction 1, 2 and 3. + pj1, pj2, pj3 : int + Degree of the domain basis functions in direction 1, 2 and 3. + starts1, starts2, starts3 : int + Starting index on the current rank, in direction 1, 2 and 3. + pads1, pads2, pads3 : int + Padding (=spline degree) for ghost regions in data, in direction 1, 2 and 3. + w1, w2, w3 : "float[:,:]" + Quadrature weights in direction 1, 2 and 3. The indexing is [global element, quadrature point]. + bi1, bi2, bi3 : "float[:,:,:,:]" + Values of codomain basis functions in direction 1, 2 and 3. The indexing is + [global element, local basis function, derivative, quadrature point]. + bj1, bj2, bj3 : "float[:,:,:,:]" + Values of domain basis functions in direction 1, 2 and 3, same indexing convention as bi1, bi2, bi3. + mat_fun : "float[:,:,:]" + Function under the integral evaluated at quadrature points (flattened in each direction). + The indexing is [flattened quad. point dir. 1, flattened quad. point dir. 2, flattened quad. point dir. 3]. + data : "float[:,:,:,:,:,:]" + _data array of StencilMatrix to store the results. + """ + + ne1 = spans1.size + ne2 = spans2.size + ne3 = spans3.size + + nq1 = shape(w1)[1] + nq2 = shape(w2)[1] + nq3 = shape(w3)[1] + + tmp_bi1 = np.zeros(nq1) + tmp_bi2 = np.zeros(nq2) + tmp_bi3 = np.zeros(nq3) + + tmp_bj1 = np.zeros(nq1) + tmp_bj2 = np.zeros(nq2) + tmp_bj3 = np.zeros(nq3) + + tmp_w1 = np.zeros(nq1) + tmp_w2 = np.zeros(nq2) + tmp_w3 = np.zeros(nq3) + + tmp_mat_fun = np.zeros((nq1, nq2, nq3)) + + for iel1 in range(ne1): + for iel2 in range(ne2): + for iel3 in range(ne3): + tmp_mat_fun[:, :, :] = mat_fun[ + iel1 * nq1 : (iel1 + 1) * nq1, + iel2 * nq2 : (iel2 + 1) * nq2, + iel3 * nq3 : (iel3 + 1) * nq3, + ] + + tmp_w1[:] = w1[iel1, :] + tmp_w2[:] = w2[iel2, :] + tmp_w3[:] = w3[iel3, :] + + for il1 in range(pi1 + 1): + for il2 in range(pi2 + 1): + for il3 in range(pi3 + 1): + tmp_bi1[:] = bi1[iel1, il1, 0, :] + tmp_bi2[:] = bi2[iel2, il2, 0, :] + tmp_bi3[:] = bi3[iel3, il3, 0, :] + + # global spline indices + i_global1 = spans1[iel1] - pi1 + il1 + i_global2 = spans2[iel2] - pi2 + il2 + i_global3 = spans3[iel3] - pi3 + il3 + + # local spline indices (- starts --> can be negative, will therefore be written to ghost regions) + i_local1 = i_global1 - starts1 + i_local2 = i_global2 - starts2 + i_local3 = i_global3 - starts3 + + for jl1 in range(pj1 + 1): + for jl2 in range(pj2 + 1): + for jl3 in range(pj3 + 1): + tmp_bj1[:] = bj1[iel1, jl1, 0, :] + tmp_bj2[:] = bj2[iel2, jl2, 0, :] + tmp_bj3[:] = bj3[iel3, jl3, 0, :] + + value = 0.0 + + for q1 in range(nq1): + for q2 in range(nq2): + for q3 in range(nq3): + wvol = ( + tmp_w1[q1] * tmp_w2[q2] * tmp_w3[q3] * tmp_mat_fun[q1, q2, q3] + ) + + bi = tmp_bi1[q1] * tmp_bi2[q2] * tmp_bi3[q3] + bj = tmp_bj1[q1] * tmp_bj2[q2] * tmp_bj3[q3] + + value += wvol * bi * bj + + data[ + pads1 + i_local1, + pads2 + i_local2, + pads3 + i_local3, + pads1 + jl1 - il1, + pads2 + jl2 - il2, + pads3 + jl3 - il3, + ] += value diff --git a/src/struphy/feec/kernels/kernel_3d_matrixfree/__init__.py b/src/struphy/feec/kernels/kernel_3d_matrixfree/__init__.py new file mode 100644 index 000000000..012ea20b8 --- /dev/null +++ b/src/struphy/feec/kernels/kernel_3d_matrixfree/__init__.py @@ -0,0 +1,11 @@ +"""The kernel of this folder: pyccel in ``_kernels.py``, CUDA in ``_cuda.cu`` once ported. + +No CUDA version yet: on the CuPy backend the pyccel kernel runs on host copies of the arrays (a warning is +emitted once), as before the move into this folder. +""" + +from cunumpy.kernels import Kernel + +from struphy.utils.cuda_arguments import CUDA_STRUCTS + +kernel_3d_matrixfree = Kernel.from_folder(__name__, structs=CUDA_STRUCTS, missing_cuda="fallback") diff --git a/src/struphy/feec/kernels/kernel_3d_matrixfree/kernel_3d_matrixfree_kernels.py b/src/struphy/feec/kernels/kernel_3d_matrixfree/kernel_3d_matrixfree_kernels.py new file mode 100644 index 000000000..4a21d26bd --- /dev/null +++ b/src/struphy/feec/kernels/kernel_3d_matrixfree/kernel_3d_matrixfree_kernels.py @@ -0,0 +1,173 @@ +"""Kernel kernel_3d_matrixfree: integral kernel for mass matrices and L2-projections.""" + +import numpy as np +from numpy import shape + + +def kernel_3d_matrixfree( + spansi1: "int[:]", + spansi2: "int[:]", + spansi3: "int[:]", + spansj1: "int[:]", + spansj2: "int[:]", + spansj3: "int[:]", + pi1: int, + pi2: int, + pi3: int, + pj1: int, + pj2: int, + pj3: int, + startsi1: int, + startsi2: int, + startsi3: int, + startsj1: int, + startsj2: int, + startsj3: int, + padsi1: int, + padsi2: int, + padsi3: int, + padsj1: int, + padsj2: int, + padsj3: int, + w1: "float[:,:]", + w2: "float[:,:]", + w3: "float[:,:]", + bi1: "float[:,:,:,:]", + bi2: "float[:,:,:,:]", + bi3: "float[:,:,:,:]", + bj1: "float[:,:,:,:]", + bj2: "float[:,:,:,:]", + bj3: "float[:,:,:,:]", + mat_fun: "float[:,:,:]", + data_out: "float[:,:,:]", + data_in: "float[:,:,:]", +): + """ + Performs the integration of Lambda_(i1, i2, i3) * mat_fun(eta1, eta2, eta3) * f(eta1, eta2, eta3) for the basis functions (i1, i2, i3) available on the calling process, + where f is the spline function represented by the coefficients in data_in. + + The results are written into data_out (attention: data_out is NOT set to zero first, but the results are added to data_out). + This computes the action of the mass matrix on a vector without ever assembling the matrix itself. + + Parameters + ---------- + spansi1, spansi2, spansi3 : array[int] + Arrays of span indices in direction 1, 2 and 3 for the codomain ("i") basis functions; the span is the + index of the last non-vanishing spline on each grid element (cell). + spansj1, spansj2, spansj3 : array[int] + Arrays of span indices in direction 1, 2 and 3 for the domain ("j") basis functions. + pi1, pi2, pi3 : int + Degree of the codomain basis functions in direction 1, 2 and 3. + pj1, pj2, pj3 : int + Degree of the domain basis functions in direction 1, 2 and 3. + startsi1, startsi2, startsi3 : int + Starting index on the current rank for the codomain basis functions, in direction 1, 2 and 3. + startsj1, startsj2, startsj3 : int + Starting index on the current rank for the domain basis functions, in direction 1, 2 and 3. + padsi1, padsi2, padsi3 : int + Padding (=spline degree) for ghost regions in data_out, in direction 1, 2 and 3. + padsj1, padsj2, padsj3 : int + Padding (=spline degree) for ghost regions in data_in, in direction 1, 2 and 3. + w1, w2, w3 : "float[:,:]" + Quadrature weights in direction 1, 2 and 3. The indexing is [global element, quadrature point]. + bi1, bi2, bi3 : "float[:,:,:,:]" + Values of codomain basis functions in direction 1, 2 and 3. The indexing is + [global element, local basis function, derivative, quadrature point]. + bj1, bj2, bj3 : "float[:,:,:,:]" + Values of domain basis functions in direction 1, 2 and 3, same indexing convention as bi1, bi2, bi3. + mat_fun : "float[:,:,:]" + Function under the integral evaluated at quadrature points (flattened in each direction). + The indexing is [flattened quad. point dir. 1, flattened quad. point dir. 2, flattened quad. point dir. 3]. + data_out : "float[:,:,:]" + _data array of StencilVector to store the results of the matrix-vector product. + data_in : "float[:,:,:]" + _data array of StencilVector holding the spline coefficients of the input function f. + """ + + ne1 = spansi1.size + ne2 = spansi2.size + ne3 = spansi3.size + + nq1 = shape(w1)[1] + nq2 = shape(w2)[1] + nq3 = shape(w3)[1] + + tmp_w1 = np.zeros(nq1) + tmp_w2 = np.zeros(nq2) + tmp_w3 = np.zeros(nq3) + + tmp_bi1 = np.zeros(pi1 + 1) + tmp_bi2 = np.zeros(pi2 + 1) + tmp_bi3 = np.zeros(pi3 + 1) + + tmp_bj1 = np.zeros(pj1 + 1) + tmp_bj2 = np.zeros(pj2 + 1) + tmp_bj3 = np.zeros(pj3 + 1) + + tmp_mat_fun = np.zeros((nq1, nq2, nq3)) + + for iel1 in range(ne1): + for iel2 in range(ne2): + for iel3 in range(ne3): + tmp_mat_fun[:, :, :] = mat_fun[ + iel1 * nq1 : (iel1 + 1) * nq1, + iel2 * nq2 : (iel2 + 1) * nq2, + iel3 * nq3 : (iel3 + 1) * nq3, + ] + + tmp_w1[:] = w1[iel1, :] + tmp_w2[:] = w2[iel2, :] + tmp_w3[:] = w3[iel3, :] + + for q1 in range(nq1): + for q2 in range(nq2): + for q3 in range(nq3): + tmp_bi1[:] = bi1[iel1, :, 0, q1] + tmp_bi2[:] = bi2[iel2, :, 0, q2] + tmp_bi3[:] = bi3[iel3, :, 0, q3] + + tmp_bj1[:] = bj1[iel1, :, 0, q1] + tmp_bj2[:] = bj2[iel2, :, 0, q2] + tmp_bj3[:] = bj3[iel3, :, 0, q3] + + bj = 0.0 + for jl1 in range(pj1 + 1): + for jl2 in range(pj2 + 1): + for jl3 in range(pj3 + 1): + # global spline indices + j_global1 = spansj1[iel1] - pj1 + jl1 + j_global2 = spansj2[iel2] - pj2 + jl2 + j_global3 = spansj3[iel3] - pj3 + jl3 + + # local spline indices (- starts --> can be negative, will therefore be written to ghost regions) + j_local1 = j_global1 - startsj1 + padsj1 + j_local2 = j_global2 - startsj2 + padsj2 + j_local3 = j_global3 - startsj3 + padsj3 + + bj += ( + tmp_bj1[jl1] + * tmp_bj2[jl2] + * tmp_bj3[jl3] + * data_in[j_local1, j_local2, j_local3] + ) + + for il1 in range(pi1 + 1): + for il2 in range(pi2 + 1): + for il3 in range(pi3 + 1): + # global spline indices + i_global1 = spansi1[iel1] - pi1 + il1 + i_global2 = spansi2[iel2] - pi2 + il2 + i_global3 = spansi3[iel3] - pi3 + il3 + + # local spline indices (- starts --> can be negative, will therefore be written to ghost regions) + i_local1 = i_global1 - startsi1 + padsi1 + i_local2 = i_global2 - startsi2 + padsi2 + i_local3 = i_global3 - startsi3 + padsi3 + + wvol = tmp_w1[q1] * tmp_w2[q2] * tmp_w3[q3] * tmp_mat_fun[q1, q2, q3] + + bi = tmp_bi1[il1] * tmp_bi2[il2] * tmp_bi3[il3] + + value = wvol * bi * bj + + data_out[i_local1, i_local2, i_local3] += value diff --git a/src/struphy/feec/kernels/kernel_3d_vec/__init__.py b/src/struphy/feec/kernels/kernel_3d_vec/__init__.py new file mode 100644 index 000000000..7cd7274a3 --- /dev/null +++ b/src/struphy/feec/kernels/kernel_3d_vec/__init__.py @@ -0,0 +1,11 @@ +"""The kernel of this folder: pyccel in ``_kernels.py``, CUDA in ``_cuda.cu`` once ported. + +No CUDA version yet: on the CuPy backend the pyccel kernel runs on host copies of the arrays (a warning is +emitted once), as before the move into this folder. +""" + +from cunumpy.kernels import Kernel + +from struphy.utils.cuda_arguments import CUDA_STRUCTS + +kernel_3d_vec = Kernel.from_folder(__name__, structs=CUDA_STRUCTS, missing_cuda="fallback") diff --git a/src/struphy/feec/kernels/kernel_3d_vec/kernel_3d_vec_kernels.py b/src/struphy/feec/kernels/kernel_3d_vec/kernel_3d_vec_kernels.py new file mode 100644 index 000000000..a4f636659 --- /dev/null +++ b/src/struphy/feec/kernels/kernel_3d_vec/kernel_3d_vec_kernels.py @@ -0,0 +1,96 @@ +"""Kernel kernel_3d_vec: integral kernel for mass matrices and L2-projections.""" + +from numpy import shape + + +def kernel_3d_vec( + spans1: "int[:]", + spans2: "int[:]", + spans3: "int[:]", + pi1: int, + pi2: int, + pi3: int, + starts1: int, + starts2: int, + starts3: int, + pads1: int, + pads2: int, + pads3: int, + w1: "float[:,:]", + w2: "float[:,:]", + w3: "float[:,:]", + bi1: "float[:,:,:,:]", + bi2: "float[:,:,:,:]", + bi3: "float[:,:,:,:]", + mat_fun: "float[:,:,:]", + data: "float[:,:,:]", +): + """ + Performs the integration of Lambda_(i1,i2,i3) * mat_fun(eta1, eta2, eta3) for the basis functions (i1,i2,i3) available on the calling process. + + The results are written into data (attention: data is NOT set to zero first, but the results are added to data). + + Parameters + ---------- + spans1, spans2, spans3 : array[int] + Arrays of span indices in direction 1, 2 and 3; the span is the index of the last non-vanishing spline + on each grid element (cell). The length of each array is the number of elements (cells) in that direction. + pi1, pi2, pi3 : int + Degree of the basis functions in direction 1, 2 and 3. + starts1, starts2, starts3 : int + Starting index on the current rank, in direction 1, 2 and 3. + pads1, pads2, pads3 : int + Padding (=spline degree) for ghost regions in data, in direction 1, 2 and 3. + w1, w2, w3 : "float[:,:]" + Quadrature weights in direction 1, 2 and 3. The indexing is [global element, quadrature point]. + bi1, bi2, bi3 : "float[:,:,:,:]" + Values of basis functions in direction 1, 2 and 3. The indexing is + [global element, local basis function, derivative, quadrature point]. + mat_fun : "float[:,:,:]" + Function under the integral evaluated at quadrature points (flattened in each direction). + The indexing is [flattened quad. point dir. 1, flattened quad. point dir. 2, flattened quad. point dir. 3]. + data : "float[:,:,:]" + _data array of StencilVector to store the results. + """ + + ne1 = spans1.size + ne2 = spans2.size + ne3 = spans3.size + + nq1 = shape(w1)[1] + nq2 = shape(w2)[1] + nq3 = shape(w3)[1] + + for iel1 in range(ne1): + for iel2 in range(ne2): + for iel3 in range(ne3): + for il1 in range(pi1 + 1): + for il2 in range(pi2 + 1): + for il3 in range(pi3 + 1): + # global spline indices + i_global1 = spans1[iel1] - pi1 + il1 + i_global2 = spans2[iel2] - pi2 + il2 + i_global3 = spans3[iel3] - pi3 + il3 + + # local spline indices (- starts --> can be negative, will therefore be written to ghost regions) + i_local1 = i_global1 - starts1 + i_local2 = i_global2 - starts2 + i_local3 = i_global3 - starts3 + + value = 0.0 + + for q1 in range(nq1): + for q2 in range(nq2): + for q3 in range(nq3): + wvol = ( + w1[iel1, q1] + * w2[iel2, q2] + * w3[iel3, q3] + * mat_fun[iel1 * nq1 + q1, iel2 * nq2 + q2, iel3 * nq3 + q3] + ) + + value += ( + wvol * bi1[iel1, il1, 0, q1] * bi2[iel2, il2, 0, q2] * bi3[iel3, il3, 0, q3] + ) + + data[pads1 + i_local1, pads2 + i_local2, pads3 + i_local3] += value diff --git a/src/struphy/feec/kernels/surface_kernel_3d_mat/__init__.py b/src/struphy/feec/kernels/surface_kernel_3d_mat/__init__.py new file mode 100644 index 000000000..014e9f236 --- /dev/null +++ b/src/struphy/feec/kernels/surface_kernel_3d_mat/__init__.py @@ -0,0 +1,11 @@ +"""The kernel of this folder: pyccel in ``_kernels.py``, CUDA in ``_cuda.cu`` once ported. + +No CUDA version yet: on the CuPy backend the pyccel kernel runs on host copies of the arrays (a warning is +emitted once), as before the move into this folder. +""" + +from cunumpy.kernels import Kernel + +from struphy.utils.cuda_arguments import CUDA_STRUCTS + +surface_kernel_3d_mat = Kernel.from_folder(__name__, structs=CUDA_STRUCTS, missing_cuda="fallback") diff --git a/src/struphy/feec/kernels/surface_kernel_3d_mat/surface_kernel_3d_mat_kernels.py b/src/struphy/feec/kernels/surface_kernel_3d_mat/surface_kernel_3d_mat_kernels.py new file mode 100644 index 000000000..f45197fab --- /dev/null +++ b/src/struphy/feec/kernels/surface_kernel_3d_mat/surface_kernel_3d_mat_kernels.py @@ -0,0 +1,170 @@ +"""Kernel surface_kernel_3d_mat: integral kernel for mass matrices and L2-projections.""" + +from numpy import shape + + +def surface_kernel_3d_mat( + spans1: "int[:]", + spans2: "int[:]", + pi0: int, + pi1: int, + pi2: int, + pj0: int, + pj1: int, + pj2: int, + starts0: int, + starts1: int, + starts2: int, + pads0: int, + pads1: int, + pads2: int, + w1: "float[:,:]", + w2: "float[:,:]", + bi1: "float[:,:,:,:]", + bi2: "float[:,:,:,:]", + bj1: "float[:,:,:,:]", + bj2: "float[:,:,:,:]", + boundary_index: int, + normal_dir: int, + mat_fun: "float[:,:]", + data: "float[:,:,:,:,:,:]", +): + """ + Assembles a boundary (surface) mass matrix: the integration of Lambda_i * mat_fun(eta_s1, eta_s2) * Lambda_j + over the boundary surface at the fixed global index boundary_index in the normal_dir direction, for the + codomain ("i") and domain ("j") basis functions available on the calling process in the two tangential + directions orthogonal to normal_dir. + + The results are written into data (attention: data is NOT set to zero first, but the results are added to data). + + Parameters + ---------- + spans1, spans2 : array[int] + Arrays of span indices in the two tangential grid directions used for the surface quadrature (as + determined by normal_dir); the span is the index of the last non-vanishing spline on each grid element. + pi0, pi1, pi2 : int + Degree of the codomain basis functions along logical axes 0, 1 and 2. + pj0, pj1, pj2 : int + Degree of the domain basis functions along logical axes 0, 1 and 2. + starts0, starts1, starts2 : int + Starting index on the current rank along logical axes 0, 1 and 2. + pads0, pads1, pads2 : int + Padding (=spline degree) for ghost regions in data, along logical axes 0, 1 and 2. + w1, w2 : "float[:,:]" + Quadrature weights in the two tangential directions. The indexing is [global element, quadrature point]. + bi1, bi2 : "float[:,:,:,:]" + Values of codomain basis functions in the two tangential directions. The indexing is + [global element, local basis function, derivative, quadrature point]. + bj1, bj2 : "float[:,:,:,:]" + Values of domain basis functions in the two tangential directions, same indexing convention as bi1, bi2. + boundary_index : int + Global index along normal_dir at which the boundary surface is located. + normal_dir : int + Logical direction (0, 1 or 2) normal to the surface; the remaining two directions are the tangential + directions used for the surface integration. + mat_fun : "float[:,:]" + Function under the integral evaluated at surface quadrature points (flattened in each tangential direction). + data : "float[:,:,:,:,:,:]" + _data array of StencilMatrix to store the results. + """ + + ne1 = spans1.size + ne2 = spans2.size + + nq1 = shape(w1)[1] + nq2 = shape(w2)[1] + + # Select the normal (n) and the two tangential (s1, s2) directions with scalars only: + # Python lists would make pyccel generate gFTL containers, which we want to avoid as a dependency. + if normal_dir == 0: + pi_s1 = pi1 + pi_s2 = pi2 + pj_s1 = pj1 + pj_s2 = pj2 + starts_n = starts0 + starts_s1 = starts1 + starts_s2 = starts2 + pads_n = pads0 + pads_s1 = pads1 + pads_s2 = pads2 + elif normal_dir == 1: + pi_s1 = pi0 + pi_s2 = pi2 + pj_s1 = pj0 + pj_s2 = pj2 + starts_n = starts1 + starts_s1 = starts0 + starts_s2 = starts2 + pads_n = pads1 + pads_s1 = pads0 + pads_s2 = pads2 + else: + pi_s1 = pi0 + pi_s2 = pi1 + pj_s1 = pj0 + pj_s2 = pj1 + starts_n = starts2 + starts_s1 = starts0 + starts_s2 = starts1 + pads_n = pads2 + pads_s1 = pads0 + pads_s2 = pads1 + + i_local_n = boundary_index - starts_n + + for iel1 in range(ne1): + for iel2 in range(ne2): + for il1 in range(pi_s1 + 1): + for il2 in range(pi_s2 + 1): + i_global1 = spans1[iel1] - pi_s1 + il1 + i_global2 = spans2[iel2] - pi_s2 + il2 + + i_local1 = i_global1 - starts_s1 + i_local2 = i_global2 - starts_s2 + + for jl1 in range(pj_s1 + 1): + for jl2 in range(pj_s2 + 1): + j_local1 = jl1 - il1 + j_local2 = jl2 - il2 + + value = 0.0 + + for q1 in range(nq1): + for q2 in range(nq2): + wvol = w1[iel1, q1] * w2[iel2, q2] * mat_fun[iel1 * nq1 + q1, iel2 * nq2 + q2] + + value += ( + wvol + * bi1[iel1, il1, 0, q1] + * bi2[iel2, il2, 0, q2] + * bj1[iel1, jl1, 0, q1] + * bj2[iel2, jl2, 0, q2] + ) + + if normal_dir == 0: + data[ + pads_n + i_local_n, + pads_s1 + i_local1, + pads_s2 + i_local2, + pads_n, + pads_s1 + j_local1, + pads_s2 + j_local2, + ] += value + elif normal_dir == 1: + data[ + pads_s1 + i_local1, + pads_n + i_local_n, + pads_s2 + i_local2, + pads_s1 + j_local1, + pads_n, + pads_s2 + j_local2, + ] += value + else: + data[ + pads_s1 + i_local1, + pads_s2 + i_local2, + pads_n + i_local_n, + pads_s1 + j_local1, + pads_s2 + j_local2, + pads_n, + ] += value diff --git a/src/struphy/feec/kernels/surface_kernel_3d_vec/__init__.py b/src/struphy/feec/kernels/surface_kernel_3d_vec/__init__.py new file mode 100644 index 000000000..7f6b69b34 --- /dev/null +++ b/src/struphy/feec/kernels/surface_kernel_3d_vec/__init__.py @@ -0,0 +1,11 @@ +"""The kernel of this folder: pyccel in ``_kernels.py``, CUDA in ``_cuda.cu`` once ported. + +No CUDA version yet: on the CuPy backend the pyccel kernel runs on host copies of the arrays (a warning is +emitted once), as before the move into this folder. +""" + +from cunumpy.kernels import Kernel + +from struphy.utils.cuda_arguments import CUDA_STRUCTS + +surface_kernel_3d_vec = Kernel.from_folder(__name__, structs=CUDA_STRUCTS, missing_cuda="fallback") diff --git a/src/struphy/feec/kernels/surface_kernel_3d_vec/surface_kernel_3d_vec_kernels.py b/src/struphy/feec/kernels/surface_kernel_3d_vec/surface_kernel_3d_vec_kernels.py new file mode 100644 index 000000000..0c584f2d0 --- /dev/null +++ b/src/struphy/feec/kernels/surface_kernel_3d_vec/surface_kernel_3d_vec_kernels.py @@ -0,0 +1,91 @@ +"""Kernel surface_kernel_3d_vec: integral kernel for mass matrices and L2-projections.""" + +from numpy import shape + + +def surface_kernel_3d_vec( + spans1: "int[:]", + spans2: "int[:]", + pi0: int, + pi1: int, + pi2: int, + starts0: int, + starts1: int, + starts2: int, + pads0: int, + pads1: int, + pads2: int, + w1: "float[:,:]", + w2: "float[:,:]", + bi1: "float[:,:,:,:]", + bi2: "float[:,:,:,:]", + boundary_index: int, + mat_fun: "float[:,:]", + data: "float[:,:,:]", +): + """ + Performs the integration of Lambda_0ij * mat_fun(eta1, eta2) over the boundary surface at the fixed + (normal-direction) global index boundary_index, for the basis functions (ij) available on the calling + process in the two surface (tangential) directions. + + The results are written into data (attention: data is NOT set to zero first, but the results are added to data). + + Parameters + ---------- + spans1, spans2 : array[int] + Arrays of span indices in the two surface (tangential) directions; the span is the index of the last + non-vanishing spline on each grid element (cell) in that direction. + pi0 : int + Degree of the basis function in the normal direction (kept for a uniform kernel signature; not used + directly since the normal index is fixed to boundary_index). + pi1, pi2 : int + Degree of the basis functions in the two surface directions. + starts0 : int + Starting index on the current rank in the normal direction. + starts1, starts2 : int + Starting index on the current rank in the two surface directions. + pads0 : int + Padding (=spline degree) for ghost regions in data, in the normal direction. + pads1, pads2 : int + Padding (=spline degree) for ghost regions in data, in the two surface directions. + w1, w2 : "float[:,:]" + Quadrature weights in the two surface directions. The indexing is [global element, quadrature point]. + bi1, bi2 : "float[:,:,:,:]" + Values of basis functions in the two surface directions. The indexing is + [global element, local basis function, derivative, quadrature point]. + boundary_index : int + Global index in the normal direction at which the boundary surface is located. + mat_fun : "float[:,:]" + Function under the integral evaluated at surface quadrature points (flattened in each surface direction). + data : "float[:,:,:]" + _data array of StencilVector to store the results; only the slice at the fixed normal index + (pads0 + i_local0) is written. + """ + + ne1 = spans1.size + ne2 = spans2.size + + nq1 = shape(w1)[1] + nq2 = shape(w2)[1] + + i_local0 = boundary_index - starts0 + + for iel1 in range(ne1): + for iel2 in range(ne2): + for il1 in range(pi1 + 1): + for il2 in range(pi2 + 1): + i_global1 = spans1[iel1] - pi1 + il1 + i_global2 = spans2[iel2] - pi2 + il2 + + i_local1 = i_global1 - starts1 + i_local2 = i_global2 - starts2 + + value = 0.0 + + for q1 in range(nq1): + for q2 in range(nq2): + wvol = w1[iel1, q1] * w2[iel2, q2] * mat_fun[iel1 * nq1 + q1, iel2 * nq2 + q2] + + value += wvol * bi1[iel1, il1, 0, q1] * bi2[iel2, il2, 0, q2] + + data[pads0 + i_local0, pads1 + i_local1, pads2 + i_local2] += value diff --git a/src/struphy/feec/mass.py b/src/struphy/feec/mass.py index e56a19af0..ff40e394b 100644 --- a/src/struphy/feec/mass.py +++ b/src/struphy/feec/mass.py @@ -5,7 +5,6 @@ from typing import Callable import cunumpy as xp -from cunumpy.kernels import PyccelKernel from feectools.api.settings import PSYDAC_BACKEND_GPYCCEL from feectools.fem.tensor import FemSpace, TensorFemSpace from feectools.fem.vector import VectorFemSpace @@ -17,7 +16,15 @@ from scope_profiler import ProfileManager from struphy import equils -from struphy.feec import mass_kernels +from struphy.feec.kernels.kernel_1d_eval import kernel_1d_eval +from struphy.feec.kernels.kernel_1d_mat import kernel_1d_mat +from struphy.feec.kernels.kernel_2d_eval import kernel_2d_eval +from struphy.feec.kernels.kernel_2d_mat import kernel_2d_mat +from struphy.feec.kernels.kernel_3d_diag import kernel_3d_diag +from struphy.feec.kernels.kernel_3d_eval import kernel_3d_eval +from struphy.feec.kernels.kernel_3d_mat import kernel_3d_mat +from struphy.feec.kernels.kernel_3d_matrixfree import kernel_3d_matrixfree +from struphy.feec.kernels.kernel_3d_vec import kernel_3d_vec from struphy.feec.linear_operators import BoundaryOperator from struphy.feec.psydac_derham import Derham, SplineFunction from struphy.feec.utilities import LocalProjectionMatrix, LocalRotationMatrix, get_quad_grids @@ -32,6 +39,12 @@ logger = logging.getLogger("struphy") +# assembly kernels by dimension of the problem (1d, 2d or 3d) +_MAT_KERNELS = {1: kernel_1d_mat, 2: kernel_2d_mat, 3: kernel_3d_mat} +_EVAL_KERNELS = {1: kernel_1d_eval, 2: kernel_2d_eval, 3: kernel_3d_eval} +_MATRIXFREE_KERNELS = {3: kernel_3d_matrixfree} +_DIAG_KERNELS = {3: kernel_3d_diag} + # space identifiers with scalar-valued and vector-valued elements _SCALAR_SPACES = ("H1", "L2") _VECTOR_SPACES = ("Hcurl", "Hdiv", "H1vec") @@ -1929,13 +1942,8 @@ def _allocate_dot_temporaries(self): self._temp_mat = self._mat.domain.zeros() def _load_assembly_kernel(self): - """Load the pyccelized assembly kernel for the dimension of the problem (1d, 2d or 3d).""" - self._assembly_kernel = PyccelKernel( - getattr( - mass_kernels, - "kernel_" + str(self._V.ldim) + "d_mat", - ), - ) + """Load the assembly kernel for the dimension of the problem (1d, 2d or 3d).""" + self._assembly_kernel = _MAT_KERNELS[self._V.ldim] @property def derham(self): @@ -2570,7 +2578,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 = _EVAL_KERNELS[W.ldim] # loop over components for a, wspace in enumerate(Wspaces): @@ -2668,19 +2676,8 @@ def __init__(self, derham, V, W, weights=None, nquads=None): self._nquads = nquads self._dtype = V.coeff_space.dtype - self._dot_kernel = PyccelKernel( - getattr( - mass_kernels, - "kernel_" + str(self._V.ldim) + "d_matrixfree", - ), - ) - - self._diag_kernel = PyccelKernel( - getattr( - mass_kernels, - "kernel_" + str(self._V.ldim) + "d_diag", - ), - ) + self._dot_kernel = _MATRIXFREE_KERNELS[self._V.ldim] + self._diag_kernel = _DIAG_KERNELS[self._V.ldim] # temporary with ghost regions for the diagonal (contributions to other processes are exchanged) self._diag_tmp = W.coeff_space.zeros() @@ -3265,7 +3262,7 @@ def get_dofs( pads = fem_space.coeff_space.pads if isinstance(vec, StencilVector): - mass_kernels.kernel_3d_vec( + kernel_3d_vec( *spans, *fem_space.degree, *starts, @@ -3276,7 +3273,7 @@ def get_dofs( vec._data, ) else: - mass_kernels.kernel_3d_vec( + kernel_3d_vec( *spans, *fem_space.degree, *starts, diff --git a/src/struphy/feec/mass_kernels.py b/src/struphy/feec/mass_kernels.py deleted file mode 100644 index 883c58209..000000000 --- a/src/struphy/feec/mass_kernels.py +++ /dev/null @@ -1,1270 +0,0 @@ -""" -Integral kernels for mass matrices and L2-projections. -""" - -import numpy as np -from numpy import shape - -# ================= 1d ================================= - - -def kernel_1d_mat( - spans1: "int[:]", - pi1: int, - pj1: int, - starts1: int, - pads1: int, - w1: "float[:,:]", - bi1: "float[:,:,:,:]", - bj1: "float[:,:,:,:]", - mat_fun: "float[:]", - data: "float[:,:]", -): - """ - Performs the integration of Lambda_(i1) * mat_fun(eta1) * Lambda_(j1) for the basis functions (i1, j1) available on the calling process. - - The results are written into data (attention: data is NOT set to zero first, but the results are added to data). - - Parameters - ---------- - spans1 : array[int] - Array of span indices; the span is the index of the last non-vanishing spline on each grid element - (cell). The length of the returned array is the number of elements (cells). - pi1 : int - Degree of the codomain basis functions. - pj1 : int - Degree of the domain basis functions. - starts1 : int - Starting index on the current rank. - pads1 : int - Padding (=spline degree) for ghost regions in data. - w1 : "float[:,:]" - Quadrature weights. The indexing is [global element, quadrature point]. - bi1 : "float[:,:,:,:]" - Values of codomain basis functions. The indexing is [global element, local basis function, derivative, quadrature point]. - bj1 : "float[:,:,:,:]" - Values of domain basis functions. The indexing is [global element, local basis function, derivative, quadrature point]. - mat_fun : "float[:]" - Function under the integral evaluated at quadrature points (flattened). - data : "float[:,:]" - _data array of StencilMatrix to store the results. - """ - - # number of elements - ne1 = spans1.size - - # number of quadrature points in each element - nq1 = shape(w1)[1] - - for iel1 in range(ne1): - for il1 in range(pi1 + 1): - # global spline indices - i_global1 = spans1[iel1] - pi1 + il1 - - # local spline indices (- starts --> can be negative, will therefore be written to ghost regions) - i_local1 = i_global1 - starts1 - - for jl1 in range(pj1 + 1): - value = 0.0 - - for q1 in range(nq1): - value += w1[iel1, q1] * bi1[iel1, il1, 0, q1] * bj1[iel1, jl1, 0, q1] * mat_fun[iel1 * nq1 + q1] - - data[pads1 + i_local1, pads1 + jl1 - il1] += value - - -def kernel_1d_vec( - spans1: "int[:]", - pi1: int, - starts1: int, - pads1: int, - w1: "float[:,:]", - bi1: "float[:,:,:,:]", - mat_fun: "float[:]", - data: "float[:]", -): - """ - Performs the integration of Lambda_(i1) * mat_fun(eta1) for the basis functions (i1) available on the calling process. - - The results are written into data (attention: data is NOT set to zero first, but the results are added to data). - - Parameters - ---------- - spans1 : array[int] - Array of span indices; the span is the index of the last non-vanishing spline on each grid element - (cell). The length of the returned array is the number of elements (cells). - pi1 : int - Degree of the basis functions. - starts1 : int - Starting index on the current rank. - pads1 : int - Padding (=spline degree) for ghost regions in data. - w1 : "float[:,:]" - Quadrature weights. The indexing is [global element, quadrature point]. - bi1 : "float[:,:,:,:]" - Values of basis functions. The indexing is [global element, local basis function, derivative, quadrature point]. - mat_fun : "float[:]" - Function under the integral evaluated at quadrature points (flattened). - data : "float[:]" - _data array of StencilVector to store the results. - """ - - ne1 = spans1.size - - nq1 = shape(w1)[1] - - for iel1 in range(ne1): - for il1 in range(pi1 + 1): - # global spline indices - i_global1 = spans1[iel1] - pi1 + il1 - - # local spline indices (- starts --> can be negative, will therefore be written to ghost regions) - i_local1 = i_global1 - starts1 - - value = 0.0 - - for q1 in range(nq1): - value += w1[iel1, q1] * bi1[iel1, il1, 0, q1] * mat_fun[iel1 * nq1 + q1] - - data[pads1 + i_local1] += value - - -def kernel_1d_eval( - spans1: "int[:]", - pi1: int, - starts1: int, - pads1: int, - bi1: "float[:,:,:,:]", - coeffs_data: "float[:]", - values: "float[:]", -): - """ - Evaluates sum_i1 [ coeffs_i1 * Lambda_i1(quad_eta1) ] for all quadrature points on the calling process. - - The results are written into values. - - Parameters - ---------- - spans1 : array[int] - Array of span indices; the span is the index of the last non-vanishing spline on each grid element - (cell). The length of the returned array is the number of elements (cells). - pi1 : int - Degree of the basis functions. - starts1 : int - Starting index on the current rank. - pads1 : int - Padding (=spline degree) for ghost regions in coeffs_data. - bi1 : "float[:,:,:,:]" - Values of basis functions. The indexing is [global element, local basis function, derivative, quadrature point]. - coeffs_data : "float[:]" - _data array of StencilVector holding the spline coefficients of the function to be evaluated. - values : "float[:]" - Output array (flattened over elements and quadrature points) holding the evaluated function values; - it is set to zero at the start of the kernel, i.e. it is overwritten, not added to. - """ - - values[:] = 0.0 - - ne1 = spans1.size - - nq1 = shape(bi1)[3] - - for iel1 in range(ne1): - for il1 in range(pi1 + 1): - # global spline indices - i_global1 = spans1[iel1] - pi1 + il1 - - # local spline indices (- starts --> can be negative, will therefore be written to ghost regions) - i_local1 = i_global1 - starts1 - - for q1 in range(nq1): - values[iel1 * nq1 + q1] += coeffs_data[pads1 + i_local1] * bi1[iel1, il1, 0, q1] - - -# ================= 2d ================================= - - -def kernel_2d_mat( - spans1: "int[:]", - spans2: "int[:]", - pi1: int, - pi2: int, - pj1: int, - pj2: int, - starts1: int, - starts2: int, - pads1: int, - pads2: int, - w1: "float[:,:]", - w2: "float[:,:]", - bi1: "float[:,:,:,:]", - bi2: "float[:,:,:,:]", - bj1: "float[:,:,:,:]", - bj2: "float[:,:,:,:]", - mat_fun: "float[:,:]", - data: "float[:,:,:,:]", -): - """ - Performs the integration of Lambda_(i1, i2) * mat_fun(eta1, eta2) * Lambda_(j1, j2) for the basis functions (i1, i2, j1, j2) available on the calling process. - - The results are written into data (attention: data is NOT set to zero first, but the results are added to data). - - Parameters - ---------- - spans1, spans2 : array[int] - Arrays of span indices in direction 1 and 2; the span is the index of the last non-vanishing spline - on each grid element (cell). The length of each array is the number of elements (cells) in that direction. - pi1, pi2 : int - Degree of the codomain basis functions in direction 1 and 2. - pj1, pj2 : int - Degree of the domain basis functions in direction 1 and 2. - starts1, starts2 : int - Starting index on the current rank, in direction 1 and 2. - pads1, pads2 : int - Padding (=spline degree) for ghost regions in data, in direction 1 and 2. - w1, w2 : "float[:,:]" - Quadrature weights in direction 1 and 2. The indexing is [global element, quadrature point]. - bi1, bi2 : "float[:,:,:,:]" - Values of codomain basis functions in direction 1 and 2. The indexing is - [global element, local basis function, derivative, quadrature point]. - bj1, bj2 : "float[:,:,:,:]" - Values of domain basis functions in direction 1 and 2, same indexing convention as bi1, bi2. - mat_fun : "float[:,:]" - Function under the integral evaluated at quadrature points (flattened in each direction). - The indexing is [flattened quadrature point in direction 1, flattened quadrature point in direction 2]. - data : "float[:,:,:,:]" - _data array of StencilMatrix to store the results. - """ - - ne1 = spans1.size - ne2 = spans2.size - - nq1 = shape(w1)[1] - nq2 = shape(w2)[1] - - for iel1 in range(ne1): - for iel2 in range(ne2): - for il1 in range(pi1 + 1): - for il2 in range(pi2 + 1): - # global spline indices - i_global1 = spans1[iel1] - pi1 + il1 - i_global2 = spans2[iel2] - pi2 + il2 - - # local spline indices (- starts --> can be negative, will therefore be written to ghost regions) - i_local1 = i_global1 - starts1 - i_local2 = i_global2 - starts2 - - for jl1 in range(pj1 + 1): - for jl2 in range(pj2 + 1): - value = 0.0 - - for q1 in range(nq1): - for q2 in range(nq2): - wvol = w1[iel1, q1] * w2[iel2, q2] * mat_fun[iel1 * nq1 + q1, iel2 * nq2 + q2] - bi = bi1[iel1, il1, 0, q1] * bi2[iel2, il2, 0, q2] - bj = bj1[iel1, jl1, 0, q1] * bj2[iel2, jl2, 0, q2] - - value += wvol * bi * bj - - data[pads1 + i_local1, pads2 + i_local2, pads1 + jl1 - il1, pads2 + jl2 - il2] += value - - -def kernel_2d_vec( - spans1: "int[:]", - spans2: "int[:]", - pi1: int, - pi2: int, - starts1: int, - starts2: int, - pads1: int, - pads2: int, - w1: "float[:,:]", - w2: "float[:,:]", - bi1: "float[:,:,:,:]", - bi2: "float[:,:,:,:]", - mat_fun: "float[:,:]", - data: "float[:,:]", -): - """ - Performs the integration of Lambda_(i1, i2) * mat_fun(eta1, eta2) for the basis functions (i1, i2) available on the calling process. - - The results are written into data (attention: data is NOT set to zero first, but the results are added to data). - - Parameters - ---------- - spans1, spans2 : array[int] - Arrays of span indices in direction 1 and 2; the span is the index of the last non-vanishing spline - on each grid element (cell). The length of each array is the number of elements (cells) in that direction. - pi1, pi2 : int - Degree of the basis functions in direction 1 and 2. - starts1, starts2 : int - Starting index on the current rank, in direction 1 and 2. - pads1, pads2 : int - Padding (=spline degree) for ghost regions in data, in direction 1 and 2. - w1, w2 : "float[:,:]" - Quadrature weights in direction 1 and 2. The indexing is [global element, quadrature point]. - bi1, bi2 : "float[:,:,:,:]" - Values of basis functions in direction 1 and 2. The indexing is - [global element, local basis function, derivative, quadrature point]. - mat_fun : "float[:,:]" - Function under the integral evaluated at quadrature points (flattened in each direction). - The indexing is [flattened quadrature point in direction 1, flattened quadrature point in direction 2]. - data : "float[:,:]" - _data array of StencilVector to store the results. - """ - - ne1 = spans1.size - ne2 = spans2.size - - nq1 = shape(w1)[1] - nq2 = shape(w2)[1] - - for iel1 in range(ne1): - for iel2 in range(ne2): - for il1 in range(pi1 + 1): - for il2 in range(pi2 + 1): - # global spline indices - i_global1 = spans1[iel1] - pi1 + il1 - i_global2 = spans2[iel2] - pi2 + il2 - - # local spline indices (- starts --> can be negative, will therefore be written to ghost regions) - i_local1 = i_global1 - starts1 - i_local2 = i_global2 - starts2 - - value = 0.0 - - for q1 in range(nq1): - for q2 in range(nq2): - wvol = w1[iel1, q1] * w2[iel2, q2] * mat_fun[iel1 * nq1 + q1, iel2 * nq2 + q2] - - value += wvol * bi1[iel1, il1, 0, q1] * bi2[iel2, il2, 0, q2] - - data[pads1 + i_local1, pads2 + i_local2] += value - - -def kernel_2d_eval( - spans1: "int[:]", - spans2: "int[:]", - pi1: int, - pi2: int, - starts1: int, - starts2: int, - pads1: int, - pads2: int, - bi1: "float[:,:,:,:]", - bi2: "float[:,:,:,:]", - coeffs_data: "float[:,:]", - values: "float[:,:]", -): - """ - Evaluates sum_(i1, i2) [ coeffs_{i1,i2} * Lambda_{i1, i2}(quad_eta1, quad_eta2) ] for all quadrature points on the calling process. - - The results are written into values. - - Parameters - ---------- - spans1, spans2 : array[int] - Arrays of span indices in direction 1 and 2; the span is the index of the last non-vanishing spline - on each grid element (cell). The length of each array is the number of elements (cells) in that direction. - pi1, pi2 : int - Degree of the basis functions in direction 1 and 2. - starts1, starts2 : int - Starting index on the current rank, in direction 1 and 2. - pads1, pads2 : int - Padding (=spline degree) for ghost regions in coeffs_data, in direction 1 and 2. - bi1, bi2 : "float[:,:,:,:]" - Values of basis functions in direction 1 and 2. The indexing is - [global element, local basis function, derivative, quadrature point]. - coeffs_data : "float[:,:]" - _data array of StencilVector holding the spline coefficients of the function to be evaluated. - values : "float[:,:]" - Output array (flattened over elements and quadrature points in each direction) holding the evaluated - function values; it is set to zero at the start of the kernel, i.e. it is overwritten, not added to. - """ - - values[:, :] = 0.0 - - ne1 = spans1.size - ne2 = spans2.size - - nq1 = shape(bi1)[3] - nq2 = shape(bi2)[3] - - for iel1 in range(ne1): - for iel2 in range(ne2): - for il1 in range(pi1 + 1): - for il2 in range(pi2 + 1): - # global spline indices - i_global1 = spans1[iel1] - pi1 + il1 - i_global2 = spans2[iel2] - pi2 + il2 - - # local spline indices (- starts --> can be negative, will therefore be written to ghost regions) - i_local1 = i_global1 - starts1 - i_local2 = i_global2 - starts2 - - for q1 in range(nq1): - for q2 in range(nq2): - values[iel1 * nq1 + q1, iel2 * nq2 + q2] += ( - coeffs_data[pads1 + i_local1, pads2 + i_local2] - * bi1[iel1, il1, 0, q1] - * bi2[iel2, il2, 0, q2] - ) - - -# ================= 3d ================================= - - -def kernel_3d_mat( - spans1: "int[:]", - spans2: "int[:]", - spans3: "int[:]", - pi1: int, - pi2: int, - pi3: int, - pj1: int, - pj2: int, - pj3: int, - starts1: int, - starts2: int, - starts3: int, - pads1: int, - pads2: int, - pads3: int, - w1: "float[:,:]", - w2: "float[:,:]", - w3: "float[:,:]", - bi1: "float[:,:,:,:]", - bi2: "float[:,:,:,:]", - bi3: "float[:,:,:,:]", - bj1: "float[:,:,:,:]", - bj2: "float[:,:,:,:]", - bj3: "float[:,:,:,:]", - mat_fun: "float[:,:,:]", - data: "float[:,:,:,:,:,:]", -): - """ - Performs the integration of Lambda_(i1,i2,i3) * mat_fun(eta1, eta2, eta3) * Lambda_(j1,j2,j3) for the basis functions (i1,i2,i3, j1,j2,j3) available on the calling process. - - The results are written into data (attention: data is NOT set to zero first, but the results are added to data). - - Parameters - ---------- - spans1, spans2, spans3 : array[int] - Arrays of span indices in direction 1, 2 and 3; the span is the index of the last non-vanishing spline - on each grid element (cell). The length of each array is the number of elements (cells) in that direction. - pi1, pi2, pi3 : int - Degree of the codomain basis functions in direction 1, 2 and 3. - pj1, pj2, pj3 : int - Degree of the domain basis functions in direction 1, 2 and 3. - starts1, starts2, starts3 : int - Starting index on the current rank, in direction 1, 2 and 3. - pads1, pads2, pads3 : int - Padding (=spline degree) for ghost regions in data, in direction 1, 2 and 3. - w1, w2, w3 : "float[:,:]" - Quadrature weights in direction 1, 2 and 3. The indexing is [global element, quadrature point]. - bi1, bi2, bi3 : "float[:,:,:,:]" - Values of codomain basis functions in direction 1, 2 and 3. The indexing is - [global element, local basis function, derivative, quadrature point]. - bj1, bj2, bj3 : "float[:,:,:,:]" - Values of domain basis functions in direction 1, 2 and 3, same indexing convention as bi1, bi2, bi3. - mat_fun : "float[:,:,:]" - Function under the integral evaluated at quadrature points (flattened in each direction). - The indexing is [flattened quad. point dir. 1, flattened quad. point dir. 2, flattened quad. point dir. 3]. - data : "float[:,:,:,:,:,:]" - _data array of StencilMatrix to store the results. - """ - - ne1 = spans1.size - ne2 = spans2.size - ne3 = spans3.size - - nq1 = shape(w1)[1] - nq2 = shape(w2)[1] - nq3 = shape(w3)[1] - - tmp_bi1 = np.zeros(nq1) - tmp_bi2 = np.zeros(nq2) - tmp_bi3 = np.zeros(nq3) - - tmp_bj1 = np.zeros(nq1) - tmp_bj2 = np.zeros(nq2) - tmp_bj3 = np.zeros(nq3) - - tmp_w1 = np.zeros(nq1) - tmp_w2 = np.zeros(nq2) - tmp_w3 = np.zeros(nq3) - - tmp_mat_fun = np.zeros((nq1, nq2, nq3)) - - for iel1 in range(ne1): - for iel2 in range(ne2): - for iel3 in range(ne3): - tmp_mat_fun[:, :, :] = mat_fun[ - iel1 * nq1 : (iel1 + 1) * nq1, - iel2 * nq2 : (iel2 + 1) * nq2, - iel3 * nq3 : (iel3 + 1) * nq3, - ] - - tmp_w1[:] = w1[iel1, :] - tmp_w2[:] = w2[iel2, :] - tmp_w3[:] = w3[iel3, :] - - for il1 in range(pi1 + 1): - for il2 in range(pi2 + 1): - for il3 in range(pi3 + 1): - tmp_bi1[:] = bi1[iel1, il1, 0, :] - tmp_bi2[:] = bi2[iel2, il2, 0, :] - tmp_bi3[:] = bi3[iel3, il3, 0, :] - - # global spline indices - i_global1 = spans1[iel1] - pi1 + il1 - i_global2 = spans2[iel2] - pi2 + il2 - i_global3 = spans3[iel3] - pi3 + il3 - - # local spline indices (- starts --> can be negative, will therefore be written to ghost regions) - i_local1 = i_global1 - starts1 - i_local2 = i_global2 - starts2 - i_local3 = i_global3 - starts3 - - for jl1 in range(pj1 + 1): - for jl2 in range(pj2 + 1): - for jl3 in range(pj3 + 1): - tmp_bj1[:] = bj1[iel1, jl1, 0, :] - tmp_bj2[:] = bj2[iel2, jl2, 0, :] - tmp_bj3[:] = bj3[iel3, jl3, 0, :] - - value = 0.0 - - for q1 in range(nq1): - for q2 in range(nq2): - for q3 in range(nq3): - wvol = ( - tmp_w1[q1] * tmp_w2[q2] * tmp_w3[q3] * tmp_mat_fun[q1, q2, q3] - ) - - bi = tmp_bi1[q1] * tmp_bi2[q2] * tmp_bi3[q3] - bj = tmp_bj1[q1] * tmp_bj2[q2] * tmp_bj3[q3] - - value += wvol * bi * bj - - data[ - pads1 + i_local1, - pads2 + i_local2, - pads3 + i_local3, - pads1 + jl1 - il1, - pads2 + jl2 - il2, - pads3 + jl3 - il3, - ] += value - - -def kernel_3d_vec( - spans1: "int[:]", - spans2: "int[:]", - spans3: "int[:]", - pi1: int, - pi2: int, - pi3: int, - starts1: int, - starts2: int, - starts3: int, - pads1: int, - pads2: int, - pads3: int, - w1: "float[:,:]", - w2: "float[:,:]", - w3: "float[:,:]", - bi1: "float[:,:,:,:]", - bi2: "float[:,:,:,:]", - bi3: "float[:,:,:,:]", - mat_fun: "float[:,:,:]", - data: "float[:,:,:]", -): - """ - Performs the integration of Lambda_(i1,i2,i3) * mat_fun(eta1, eta2, eta3) for the basis functions (i1,i2,i3) available on the calling process. - - The results are written into data (attention: data is NOT set to zero first, but the results are added to data). - - Parameters - ---------- - spans1, spans2, spans3 : array[int] - Arrays of span indices in direction 1, 2 and 3; the span is the index of the last non-vanishing spline - on each grid element (cell). The length of each array is the number of elements (cells) in that direction. - pi1, pi2, pi3 : int - Degree of the basis functions in direction 1, 2 and 3. - starts1, starts2, starts3 : int - Starting index on the current rank, in direction 1, 2 and 3. - pads1, pads2, pads3 : int - Padding (=spline degree) for ghost regions in data, in direction 1, 2 and 3. - w1, w2, w3 : "float[:,:]" - Quadrature weights in direction 1, 2 and 3. The indexing is [global element, quadrature point]. - bi1, bi2, bi3 : "float[:,:,:,:]" - Values of basis functions in direction 1, 2 and 3. The indexing is - [global element, local basis function, derivative, quadrature point]. - mat_fun : "float[:,:,:]" - Function under the integral evaluated at quadrature points (flattened in each direction). - The indexing is [flattened quad. point dir. 1, flattened quad. point dir. 2, flattened quad. point dir. 3]. - data : "float[:,:,:]" - _data array of StencilVector to store the results. - """ - - ne1 = spans1.size - ne2 = spans2.size - ne3 = spans3.size - - nq1 = shape(w1)[1] - nq2 = shape(w2)[1] - nq3 = shape(w3)[1] - - for iel1 in range(ne1): - for iel2 in range(ne2): - for iel3 in range(ne3): - for il1 in range(pi1 + 1): - for il2 in range(pi2 + 1): - for il3 in range(pi3 + 1): - # global spline indices - i_global1 = spans1[iel1] - pi1 + il1 - i_global2 = spans2[iel2] - pi2 + il2 - i_global3 = spans3[iel3] - pi3 + il3 - - # local spline indices (- starts --> can be negative, will therefore be written to ghost regions) - i_local1 = i_global1 - starts1 - i_local2 = i_global2 - starts2 - i_local3 = i_global3 - starts3 - - value = 0.0 - - for q1 in range(nq1): - for q2 in range(nq2): - for q3 in range(nq3): - wvol = ( - w1[iel1, q1] - * w2[iel2, q2] - * w3[iel3, q3] - * mat_fun[iel1 * nq1 + q1, iel2 * nq2 + q2, iel3 * nq3 + q3] - ) - - value += ( - wvol * bi1[iel1, il1, 0, q1] * bi2[iel2, il2, 0, q2] * bi3[iel3, il3, 0, q3] - ) - - data[pads1 + i_local1, pads2 + i_local2, pads3 + i_local3] += value - - -def kernel_3d_eval( - spans1: "int[:]", - spans2: "int[:]", - spans3: "int[:]", - pi1: int, - pi2: int, - pi3: int, - starts1: int, - starts2: int, - starts3: int, - pads1: int, - pads2: int, - pads3: int, - bi1: "float[:,:,:,:]", - bi2: "float[:,:,:,:]", - bi3: "float[:,:,:,:]", - coeffs_data: "float[:,:,:]", - values: "float[:,:,:]", -): - """ - Evaluates sum_(i1,i2,i3) [ coeffs_{i1,i2,i3} * Lambda_{i1,i2,i3}(quad_eta1, quad_eta2, quad_eta3) ] for all quadrature points on the calling process. - - The results are written into values. - - Parameters - ---------- - spans1, spans2, spans3 : array[int] - Arrays of span indices in direction 1, 2 and 3; the span is the index of the last non-vanishing spline - on each grid element (cell). The length of each array is the number of elements (cells) in that direction. - pi1, pi2, pi3 : int - Degree of the basis functions in direction 1, 2 and 3. - starts1, starts2, starts3 : int - Starting index on the current rank, in direction 1, 2 and 3. - pads1, pads2, pads3 : int - Padding (=spline degree) for ghost regions in coeffs_data, in direction 1, 2 and 3. - bi1, bi2, bi3 : "float[:,:,:,:]" - Values of basis functions in direction 1, 2 and 3. The indexing is - [global element, local basis function, derivative, quadrature point]. - coeffs_data : "float[:,:,:]" - _data array of StencilVector holding the spline coefficients of the function to be evaluated. - values : "float[:,:,:]" - Output array (flattened over elements and quadrature points in each direction) holding the evaluated - function values; it is set to zero at the start of the kernel, i.e. it is overwritten, not added to. - """ - - values[:, :, :] = 0.0 - - ne1 = spans1.size - ne2 = spans2.size - ne3 = spans3.size - - nq1 = shape(bi1)[3] - nq2 = shape(bi2)[3] - nq3 = shape(bi3)[3] - - for iel1 in range(ne1): - for iel2 in range(ne2): - for iel3 in range(ne3): - for il1 in range(pi1 + 1): - for il2 in range(pi2 + 1): - for il3 in range(pi3 + 1): - # global spline indices - i_global1 = spans1[iel1] - pi1 + il1 - i_global2 = spans2[iel2] - pi2 + il2 - i_global3 = spans3[iel3] - pi3 + il3 - - # local spline indices (- starts --> can be negative, will therefore be written to ghost regions) - i_local1 = i_global1 - starts1 - i_local2 = i_global2 - starts2 - i_local3 = i_global3 - starts3 - - for q1 in range(nq1): - for q2 in range(nq2): - for q3 in range(nq3): - values[iel1 * nq1 + q1, iel2 * nq2 + q2, iel3 * nq3 + q3] += ( - coeffs_data[pads1 + i_local1, pads2 + i_local2, pads3 + i_local3] - * bi1[iel1, il1, 0, q1] - * bi2[iel2, il2, 0, q2] - * bi3[iel3, il3, 0, q3] - ) - - -def kernel_3d_matrixfree( - spansi1: "int[:]", - spansi2: "int[:]", - spansi3: "int[:]", - spansj1: "int[:]", - spansj2: "int[:]", - spansj3: "int[:]", - pi1: int, - pi2: int, - pi3: int, - pj1: int, - pj2: int, - pj3: int, - startsi1: int, - startsi2: int, - startsi3: int, - startsj1: int, - startsj2: int, - startsj3: int, - padsi1: int, - padsi2: int, - padsi3: int, - padsj1: int, - padsj2: int, - padsj3: int, - w1: "float[:,:]", - w2: "float[:,:]", - w3: "float[:,:]", - bi1: "float[:,:,:,:]", - bi2: "float[:,:,:,:]", - bi3: "float[:,:,:,:]", - bj1: "float[:,:,:,:]", - bj2: "float[:,:,:,:]", - bj3: "float[:,:,:,:]", - mat_fun: "float[:,:,:]", - data_out: "float[:,:,:]", - data_in: "float[:,:,:]", -): - """ - Performs the integration of Lambda_(i1, i2, i3) * mat_fun(eta1, eta2, eta3) * f(eta1, eta2, eta3) for the basis functions (i1, i2, i3) available on the calling process, - where f is the spline function represented by the coefficients in data_in. - - The results are written into data_out (attention: data_out is NOT set to zero first, but the results are added to data_out). - This computes the action of the mass matrix on a vector without ever assembling the matrix itself. - - Parameters - ---------- - spansi1, spansi2, spansi3 : array[int] - Arrays of span indices in direction 1, 2 and 3 for the codomain ("i") basis functions; the span is the - index of the last non-vanishing spline on each grid element (cell). - spansj1, spansj2, spansj3 : array[int] - Arrays of span indices in direction 1, 2 and 3 for the domain ("j") basis functions. - pi1, pi2, pi3 : int - Degree of the codomain basis functions in direction 1, 2 and 3. - pj1, pj2, pj3 : int - Degree of the domain basis functions in direction 1, 2 and 3. - startsi1, startsi2, startsi3 : int - Starting index on the current rank for the codomain basis functions, in direction 1, 2 and 3. - startsj1, startsj2, startsj3 : int - Starting index on the current rank for the domain basis functions, in direction 1, 2 and 3. - padsi1, padsi2, padsi3 : int - Padding (=spline degree) for ghost regions in data_out, in direction 1, 2 and 3. - padsj1, padsj2, padsj3 : int - Padding (=spline degree) for ghost regions in data_in, in direction 1, 2 and 3. - w1, w2, w3 : "float[:,:]" - Quadrature weights in direction 1, 2 and 3. The indexing is [global element, quadrature point]. - bi1, bi2, bi3 : "float[:,:,:,:]" - Values of codomain basis functions in direction 1, 2 and 3. The indexing is - [global element, local basis function, derivative, quadrature point]. - bj1, bj2, bj3 : "float[:,:,:,:]" - Values of domain basis functions in direction 1, 2 and 3, same indexing convention as bi1, bi2, bi3. - mat_fun : "float[:,:,:]" - Function under the integral evaluated at quadrature points (flattened in each direction). - The indexing is [flattened quad. point dir. 1, flattened quad. point dir. 2, flattened quad. point dir. 3]. - data_out : "float[:,:,:]" - _data array of StencilVector to store the results of the matrix-vector product. - data_in : "float[:,:,:]" - _data array of StencilVector holding the spline coefficients of the input function f. - """ - - ne1 = spansi1.size - ne2 = spansi2.size - ne3 = spansi3.size - - nq1 = shape(w1)[1] - nq2 = shape(w2)[1] - nq3 = shape(w3)[1] - - tmp_w1 = np.zeros(nq1) - tmp_w2 = np.zeros(nq2) - tmp_w3 = np.zeros(nq3) - - tmp_bi1 = np.zeros(pi1 + 1) - tmp_bi2 = np.zeros(pi2 + 1) - tmp_bi3 = np.zeros(pi3 + 1) - - tmp_bj1 = np.zeros(pj1 + 1) - tmp_bj2 = np.zeros(pj2 + 1) - tmp_bj3 = np.zeros(pj3 + 1) - - tmp_mat_fun = np.zeros((nq1, nq2, nq3)) - - for iel1 in range(ne1): - for iel2 in range(ne2): - for iel3 in range(ne3): - tmp_mat_fun[:, :, :] = mat_fun[ - iel1 * nq1 : (iel1 + 1) * nq1, - iel2 * nq2 : (iel2 + 1) * nq2, - iel3 * nq3 : (iel3 + 1) * nq3, - ] - - tmp_w1[:] = w1[iel1, :] - tmp_w2[:] = w2[iel2, :] - tmp_w3[:] = w3[iel3, :] - - for q1 in range(nq1): - for q2 in range(nq2): - for q3 in range(nq3): - tmp_bi1[:] = bi1[iel1, :, 0, q1] - tmp_bi2[:] = bi2[iel2, :, 0, q2] - tmp_bi3[:] = bi3[iel3, :, 0, q3] - - tmp_bj1[:] = bj1[iel1, :, 0, q1] - tmp_bj2[:] = bj2[iel2, :, 0, q2] - tmp_bj3[:] = bj3[iel3, :, 0, q3] - - bj = 0.0 - for jl1 in range(pj1 + 1): - for jl2 in range(pj2 + 1): - for jl3 in range(pj3 + 1): - # global spline indices - j_global1 = spansj1[iel1] - pj1 + jl1 - j_global2 = spansj2[iel2] - pj2 + jl2 - j_global3 = spansj3[iel3] - pj3 + jl3 - - # local spline indices (- starts --> can be negative, will therefore be written to ghost regions) - j_local1 = j_global1 - startsj1 + padsj1 - j_local2 = j_global2 - startsj2 + padsj2 - j_local3 = j_global3 - startsj3 + padsj3 - - bj += ( - tmp_bj1[jl1] - * tmp_bj2[jl2] - * tmp_bj3[jl3] - * data_in[j_local1, j_local2, j_local3] - ) - - for il1 in range(pi1 + 1): - for il2 in range(pi2 + 1): - for il3 in range(pi3 + 1): - # global spline indices - i_global1 = spansi1[iel1] - pi1 + il1 - i_global2 = spansi2[iel2] - pi2 + il2 - i_global3 = spansi3[iel3] - pi3 + il3 - - # local spline indices (- starts --> can be negative, will therefore be written to ghost regions) - i_local1 = i_global1 - startsi1 + padsi1 - i_local2 = i_global2 - startsi2 + padsi2 - i_local3 = i_global3 - startsi3 + padsi3 - - wvol = tmp_w1[q1] * tmp_w2[q2] * tmp_w3[q3] * tmp_mat_fun[q1, q2, q3] - - bi = tmp_bi1[il1] * tmp_bi2[il2] * tmp_bi3[il3] - - value = wvol * bi * bj - - data_out[i_local1, i_local2, i_local3] += value - - -def kernel_3d_diag( - spans1: "int[:]", - spans2: "int[:]", - spans3: "int[:]", - pi1: int, - pi2: int, - pi3: int, - starts1: int, - starts2: int, - starts3: int, - pads1: int, - pads2: int, - pads3: int, - w1: "float[:,:]", - w2: "float[:,:]", - w3: "float[:,:]", - bi1: "float[:,:,:,:]", - bi2: "float[:,:,:,:]", - bi3: "float[:,:,:,:]", - mat_fun: "float[:,:,:]", - data: "float[:,:,:]", -): - """ - Computes the diagonal of a mass matrix, assuming that the domain and the codomain are the same. - - The results are written into data (attention: data is NOT set to zero first, but the results are added to data). - - Parameters - ---------- - spans1, spans2, spans3 : array[int] - Arrays of span indices in direction 1, 2 and 3; the span is the index of the last non-vanishing spline - on each grid element (cell). The length of each array is the number of elements (cells) in that direction. - pi1, pi2, pi3 : int - Degree of the basis functions in direction 1, 2 and 3. - starts1, starts2, starts3 : int - Starting index on the current rank, in direction 1, 2 and 3. - pads1, pads2, pads3 : int - Padding (=spline degree) for ghost regions, in direction 1, 2 and 3. - w1, w2, w3 : "float[:,:]" - Quadrature weights in direction 1, 2 and 3. The indexing is [global element, quadrature point]. - bi1, bi2, bi3 : "float[:,:,:,:]" - Values of basis functions in direction 1, 2 and 3. The indexing is - [global element, local basis function, derivative, quadrature point]. - mat_fun : "float[:,:,:]" - Function under the integral evaluated at quadrature points (flattened in each direction). - The indexing is [flattened quad. point dir. 1, flattened quad. point dir. 2, flattened quad. point dir. 3]. - data : "float[:,:,:]" - _data array of a StencilVector (with ghost regions) to store the results. Contributions to entries owned - by other processes (or periodic images) are written to the ghost regions and must be added to their - owners afterwards with exchange_assembly_data(). - """ - - ne1 = spans1.size - ne2 = spans2.size - ne3 = spans3.size - - nq1 = shape(w1)[1] - nq2 = shape(w2)[1] - nq3 = shape(w3)[1] - - tmp_bi1 = np.zeros(nq1) - tmp_bi2 = np.zeros(nq2) - tmp_bi3 = np.zeros(nq3) - - tmp_w1 = np.zeros(nq1) - tmp_w2 = np.zeros(nq2) - tmp_w3 = np.zeros(nq3) - - tmp_mat_fun = np.zeros((nq1, nq2, nq3)) - - for iel1 in range(ne1): - for iel2 in range(ne2): - for iel3 in range(ne3): - tmp_mat_fun[:, :, :] = mat_fun[ - iel1 * nq1 : (iel1 + 1) * nq1, - iel2 * nq2 : (iel2 + 1) * nq2, - iel3 * nq3 : (iel3 + 1) * nq3, - ] - - tmp_w1[:] = w1[iel1, :] - tmp_w2[:] = w2[iel2, :] - tmp_w3[:] = w3[iel3, :] - - for il1 in range(pi1 + 1): - for il2 in range(pi2 + 1): - for il3 in range(pi3 + 1): - tmp_bi1[:] = bi1[iel1, il1, 0, :] - tmp_bi2[:] = bi2[iel2, il2, 0, :] - tmp_bi3[:] = bi3[iel3, il3, 0, :] - - # global spline indices - i_global1 = spans1[iel1] - pi1 + il1 - i_global2 = spans2[iel2] - pi2 + il2 - i_global3 = spans3[iel3] - pi3 + il3 - - # local spline indices (- starts --> can be negative, will therefore be written to ghost regions) - i_local1 = i_global1 - starts1 + pads1 - i_local2 = i_global2 - starts2 + pads2 - i_local3 = i_global3 - starts3 + pads3 - - value = 0.0 - - for q1 in range(nq1): - for q2 in range(nq2): - for q3 in range(nq3): - wvol = tmp_w1[q1] * tmp_w2[q2] * tmp_w3[q3] * tmp_mat_fun[q1, q2, q3] - - bi = tmp_bi1[q1] * tmp_bi2[q2] * tmp_bi3[q3] - - value += wvol * bi * bi - - data[i_local1, i_local2, i_local3] += value - - -def surface_kernel_3d_vec( - spans1: "int[:]", - spans2: "int[:]", - pi0: int, - pi1: int, - pi2: int, - starts0: int, - starts1: int, - starts2: int, - pads0: int, - pads1: int, - pads2: int, - w1: "float[:,:]", - w2: "float[:,:]", - bi1: "float[:,:,:,:]", - bi2: "float[:,:,:,:]", - boundary_index: int, - mat_fun: "float[:,:]", - data: "float[:,:,:]", -): - """ - Performs the integration of Lambda_0ij * mat_fun(eta1, eta2) over the boundary surface at the fixed - (normal-direction) global index boundary_index, for the basis functions (ij) available on the calling - process in the two surface (tangential) directions. - - The results are written into data (attention: data is NOT set to zero first, but the results are added to data). - - Parameters - ---------- - spans1, spans2 : array[int] - Arrays of span indices in the two surface (tangential) directions; the span is the index of the last - non-vanishing spline on each grid element (cell) in that direction. - pi0 : int - Degree of the basis function in the normal direction (kept for a uniform kernel signature; not used - directly since the normal index is fixed to boundary_index). - pi1, pi2 : int - Degree of the basis functions in the two surface directions. - starts0 : int - Starting index on the current rank in the normal direction. - starts1, starts2 : int - Starting index on the current rank in the two surface directions. - pads0 : int - Padding (=spline degree) for ghost regions in data, in the normal direction. - pads1, pads2 : int - Padding (=spline degree) for ghost regions in data, in the two surface directions. - w1, w2 : "float[:,:]" - Quadrature weights in the two surface directions. The indexing is [global element, quadrature point]. - bi1, bi2 : "float[:,:,:,:]" - Values of basis functions in the two surface directions. The indexing is - [global element, local basis function, derivative, quadrature point]. - boundary_index : int - Global index in the normal direction at which the boundary surface is located. - mat_fun : "float[:,:]" - Function under the integral evaluated at surface quadrature points (flattened in each surface direction). - data : "float[:,:,:]" - _data array of StencilVector to store the results; only the slice at the fixed normal index - (pads0 + i_local0) is written. - """ - - ne1 = spans1.size - ne2 = spans2.size - - nq1 = shape(w1)[1] - nq2 = shape(w2)[1] - - i_local0 = boundary_index - starts0 - - for iel1 in range(ne1): - for iel2 in range(ne2): - for il1 in range(pi1 + 1): - for il2 in range(pi2 + 1): - i_global1 = spans1[iel1] - pi1 + il1 - i_global2 = spans2[iel2] - pi2 + il2 - - i_local1 = i_global1 - starts1 - i_local2 = i_global2 - starts2 - - value = 0.0 - - for q1 in range(nq1): - for q2 in range(nq2): - wvol = w1[iel1, q1] * w2[iel2, q2] * mat_fun[iel1 * nq1 + q1, iel2 * nq2 + q2] - - value += wvol * bi1[iel1, il1, 0, q1] * bi2[iel2, il2, 0, q2] - - data[pads0 + i_local0, pads1 + i_local1, pads2 + i_local2] += value - - -def surface_kernel_3d_mat( - spans1: "int[:]", - spans2: "int[:]", - pi0: int, - pi1: int, - pi2: int, - pj0: int, - pj1: int, - pj2: int, - starts0: int, - starts1: int, - starts2: int, - pads0: int, - pads1: int, - pads2: int, - w1: "float[:,:]", - w2: "float[:,:]", - bi1: "float[:,:,:,:]", - bi2: "float[:,:,:,:]", - bj1: "float[:,:,:,:]", - bj2: "float[:,:,:,:]", - boundary_index: int, - normal_dir: int, - mat_fun: "float[:,:]", - data: "float[:,:,:,:,:,:]", -): - """ - Assembles a boundary (surface) mass matrix: the integration of Lambda_i * mat_fun(eta_s1, eta_s2) * Lambda_j - over the boundary surface at the fixed global index boundary_index in the normal_dir direction, for the - codomain ("i") and domain ("j") basis functions available on the calling process in the two tangential - directions orthogonal to normal_dir. - - The results are written into data (attention: data is NOT set to zero first, but the results are added to data). - - Parameters - ---------- - spans1, spans2 : array[int] - Arrays of span indices in the two tangential grid directions used for the surface quadrature (as - determined by normal_dir); the span is the index of the last non-vanishing spline on each grid element. - pi0, pi1, pi2 : int - Degree of the codomain basis functions along logical axes 0, 1 and 2. - pj0, pj1, pj2 : int - Degree of the domain basis functions along logical axes 0, 1 and 2. - starts0, starts1, starts2 : int - Starting index on the current rank along logical axes 0, 1 and 2. - pads0, pads1, pads2 : int - Padding (=spline degree) for ghost regions in data, along logical axes 0, 1 and 2. - w1, w2 : "float[:,:]" - Quadrature weights in the two tangential directions. The indexing is [global element, quadrature point]. - bi1, bi2 : "float[:,:,:,:]" - Values of codomain basis functions in the two tangential directions. The indexing is - [global element, local basis function, derivative, quadrature point]. - bj1, bj2 : "float[:,:,:,:]" - Values of domain basis functions in the two tangential directions, same indexing convention as bi1, bi2. - boundary_index : int - Global index along normal_dir at which the boundary surface is located. - normal_dir : int - Logical direction (0, 1 or 2) normal to the surface; the remaining two directions are the tangential - directions used for the surface integration. - mat_fun : "float[:,:]" - Function under the integral evaluated at surface quadrature points (flattened in each tangential direction). - data : "float[:,:,:,:,:,:]" - _data array of StencilMatrix to store the results. - """ - - ne1 = spans1.size - ne2 = spans2.size - - nq1 = shape(w1)[1] - nq2 = shape(w2)[1] - - # Select the normal (n) and the two tangential (s1, s2) directions with scalars only: - # Python lists would make pyccel generate gFTL containers, which we want to avoid as a dependency. - if normal_dir == 0: - pi_s1 = pi1 - pi_s2 = pi2 - pj_s1 = pj1 - pj_s2 = pj2 - starts_n = starts0 - starts_s1 = starts1 - starts_s2 = starts2 - pads_n = pads0 - pads_s1 = pads1 - pads_s2 = pads2 - elif normal_dir == 1: - pi_s1 = pi0 - pi_s2 = pi2 - pj_s1 = pj0 - pj_s2 = pj2 - starts_n = starts1 - starts_s1 = starts0 - starts_s2 = starts2 - pads_n = pads1 - pads_s1 = pads0 - pads_s2 = pads2 - else: - pi_s1 = pi0 - pi_s2 = pi1 - pj_s1 = pj0 - pj_s2 = pj1 - starts_n = starts2 - starts_s1 = starts0 - starts_s2 = starts1 - pads_n = pads2 - pads_s1 = pads0 - pads_s2 = pads1 - - i_local_n = boundary_index - starts_n - - for iel1 in range(ne1): - for iel2 in range(ne2): - for il1 in range(pi_s1 + 1): - for il2 in range(pi_s2 + 1): - i_global1 = spans1[iel1] - pi_s1 + il1 - i_global2 = spans2[iel2] - pi_s2 + il2 - - i_local1 = i_global1 - starts_s1 - i_local2 = i_global2 - starts_s2 - - for jl1 in range(pj_s1 + 1): - for jl2 in range(pj_s2 + 1): - j_local1 = jl1 - il1 - j_local2 = jl2 - il2 - - value = 0.0 - - for q1 in range(nq1): - for q2 in range(nq2): - wvol = w1[iel1, q1] * w2[iel2, q2] * mat_fun[iel1 * nq1 + q1, iel2 * nq2 + q2] - - value += ( - wvol - * bi1[iel1, il1, 0, q1] - * bi2[iel2, il2, 0, q2] - * bj1[iel1, jl1, 0, q1] - * bj2[iel2, jl2, 0, q2] - ) - - if normal_dir == 0: - data[ - pads_n + i_local_n, - pads_s1 + i_local1, - pads_s2 + i_local2, - pads_n, - pads_s1 + j_local1, - pads_s2 + j_local2, - ] += value - elif normal_dir == 1: - data[ - pads_s1 + i_local1, - pads_n + i_local_n, - pads_s2 + i_local2, - pads_s1 + j_local1, - pads_n, - pads_s2 + j_local2, - ] += value - else: - data[ - pads_s1 + i_local1, - pads_s2 + i_local2, - pads_n + i_local_n, - pads_s1 + j_local1, - pads_s2 + j_local2, - pads_n, - ] += value diff --git a/src/struphy/pic/tests/test_kernel_backends.py b/src/struphy/pic/tests/test_kernel_backends.py index 41e64dade..67b6c35e7 100644 --- a/src/struphy/pic/tests/test_kernel_backends.py +++ b/src/struphy/pic/tests/test_kernel_backends.py @@ -182,7 +182,7 @@ def test_catalog_signatures(): "pic.sph.kernels": 4, "bsplines.kernels": 3, "geometry.kernels": 4, - "feec.kernels": 1, + "feec.kernels": 17, "feec.local_projectors.kernels": 8, }