diff --git a/CUDA_STRATEGY.md b/CUDA_STRATEGY.md index e81695fc7..218b346ba 100644 --- a/CUDA_STRATEGY.md +++ b/CUDA_STRATEGY.md @@ -28,6 +28,7 @@ Progress is tracked in struphy-hub/struphy#650. - [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 were host-only at this stage; their current CuPy support is described in [Polar splines on CuPy](#polar-splines-on-cupy-695). See also the [PR 19 implementation notes](#pr-19-implementation-notes). - [x] **`linear_vlasov_ampere` on the GPU** (part of #688) — the first matrix accumulation: CUDA version of `linear_vlasov_ampere` (`EfieldWeightsCoupling`) on cunumpy's `Array6D` views (cunumpy 0.6.1), device versions of `fill_mat`, `fill_mat_vec`, `m_v_fill_b_v1_symm` and `outer`, and the matrix path of `Accumulator` checked on the CuPy backend (device data, ghost regions, transposed blocks, Schur operator) without a GPU through cunumpy's fake CuPy with emulated kernel launches (see the [implementation notes](#linear_vlasov_ampere-implementation-notes)). - [ ] **Diagnostics on the device** (struphy-hub/struphy#697) — scalar quantities, marker binning, saved markers and spline evaluation for output run on the device; the scalars reach the host in one copy per update (`Scalars.to_host`), every other saved dataset in one counted copy per output step (`DataContainer.save_data`), and the time state is host bookkeeping. CUDA version of `eval_spline_mpi_tensor_product_fixed` (`SplineFunction.eval_tp_fixed_loc`), moved into its own folder. Marker diagnostics kernels (`pic/diagnostics/kernels`) and SPH kernel density plots are follow-ups (see [Diagnostics on the device](#diagnostics-on-the-device-697)). +- [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)**: one complete model (`VlasovAmpereOneSpecies`) on the GPU end to end, which needs the [feectools stack](#feectools) merged and released first; the guiding-center kernels (step 3), by hand (decided in #690, see [Hand-written CUDA](#hand-written-cuda-decision-for-690)). - [ ] **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. @@ -66,7 +67,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()`, `write_local_projectors_header()` and `write_spline_header()` | | `src/struphy/kernel_arguments/pusher_args.cuh`, `local_projectors_args.cuh`, `spline_args.cuh` | the C structs `MarkerArgs`, `DerhamArgs`, `DomainArgs`, `LocalProjectorsArgs` and `SplineArgs`, 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`). Sixteen have CUDA versions (`linear_vlasov_ampere` the latest), 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, 4 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`). Eighteen have CUDA versions (the latest `linear_vlasov_ampere`, `vlasov_maxwell`, `bstar_parallel_3form` and `eval_spline_mpi_tensor_product_fixed`), 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 | @@ -865,6 +866,23 @@ What a CuPy run copies to the host per output step, and where: - **Follow-ups**: the marker diagnostics kernels in `pic/diagnostics/kernels/` (5D energies, magnetic moments, canonical momenta, `thermal_energy`; used by the guiding-center and hybrid models) have no CUDA version; SPH kernel density plots need the SPH kernels (step 6). + +## 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. extraction/DOF operators, derivatives and transposes, vector arithmetic, all three projector solvers, transposed and boundary-condition solves, mass matrices/preconditioners, callable and sampled weights, and spline coefficient extraction/restart. Repeated Kronecker and preconditioner applications are checked 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 dea89f68b..c939235e4 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 f1364f368..a22d9e05f 100644 --- a/src/struphy/pic/tests/test_kernel_backends.py +++ b/src/struphy/pic/tests/test_kernel_backends.py @@ -198,7 +198,7 @@ def test_catalog_signatures(): "pic.sph.kernels": 4, "bsplines.kernels": 4, "geometry.kernels": 4, - "feec.kernels": 1, + "feec.kernels": 17, "feec.local_projectors.kernels": 8, }