diff --git a/feectools b/feectools index f091889ea..9abb31730 160000 --- a/feectools +++ b/feectools @@ -1 +1 @@ -Subproject commit f091889ea41bee9b7cc97dbec7b9f95ce1f4384c +Subproject commit 9abb31730b281d482bcfd8061bf4475c6cfb4ade diff --git a/src/struphy/feec/boundary_mass.py b/src/struphy/feec/boundary_mass.py new file mode 100644 index 000000000..77064786b --- /dev/null +++ b/src/struphy/feec/boundary_mass.py @@ -0,0 +1,628 @@ +import logging +from typing import Callable + +import cunumpy as xp +from feectools.api.settings import PSYDAC_BACKEND_GPYCCEL +from feectools.linalg.block import BlockLinearOperator, BlockVector +from feectools.linalg.stencil import StencilMatrix, StencilVector + +from struphy.feec import mass_kernels +from struphy.feec.linear_operators import LinOpWithTransp +from struphy.feec.mass import WeightedMassOperators +from struphy.feec.psydac_derham import Derham, SplineFunction +from struphy.geometry.base import Domain +from struphy.utils.pyccel import Pyccelkernel + +logger = logging.getLogger("struphy") + + +class BoundaryIntegralOperators: + """ + Collection of boundary integral operators and boundary mass operators + for the H1, H(curl) and H(div) spaces. + + Analogous to WeightedMassOperators but for surface integrals. + + Parameters + ---------- + mass_ops : WeightedMassOperators + Mass operators object, contains geometry and derham. + """ + + def __init__( + self, + mass_ops: WeightedMassOperators, + active_faces: list[bool] | None = None, + ): + + self._mass_ops = mass_ops + self._derham = mass_ops.derham + self._domain = mass_ops.domain + + # shared surface setup for all spaces + # active faces based on bcs + if active_faces is not None: + # use provided active faces directly + self._active_faces = active_faces + else: + # default: integrate on free faces based on bcs + self._active_faces = [] + for face_idx in range(6): + normal_dir = face_idx % 3 + bc = self._derham.bcs[normal_dir] + if bc is None: + self._active_faces.append(False) + elif face_idx < 3: + self._active_faces.append(bc[0] == "free") + else: + self._active_faces.append(bc[1] == "free") + + # TODO: shared surface quad grids, geom weights, spans, wts, bases + # for each space (H1, Hcurl, Hdiv) — these differ because the + # quadrature grids are different for each space + + ################################################## + # H1 boundary operators (scalar, normal trace) # + ################################################## + + @property + def S0(self) -> "BoundaryMassOperatorH1": + """ + Boundary mass matrix for H1: + + S0_{ijk,lmn} = int_{partial Omega} Lambda^0_{ijk} Lambda^0_{lmn} sqrt(g) |DF^-T n| dS + """ + if not hasattr(self, "_S0"): + self._S0 = BoundaryMassOperatorH1(self._mass_ops, self._active_faces) + return self._S0 + + ################################################## + # H(curl) boundary operators (tangential trace) # + ################################################## + + @property + def S1(self) -> "BoundaryMassOperatorHCurl": + """ + Boundary mass matrix for H(curl): + + S1_{(mu,ijk),(nu,lmn)} = int_{partial Omega} (Lambda^1_{mu,ijk} x n) . Lambda^1_{nu,lmn} sqrt(g) |DF^-T n| dS + + Encodes the bilinear form for the tangential trace u x n against H(curl) test functions. + """ + if not hasattr(self, "_S1"): + self._S1 = BoundaryMassOperatorHCurl(self._mass_ops, self._active_faces) + return self._S1 + + +class BoundaryMassOperatorH1(LinOpWithTransp): + """ + Assembles the boundary mass matrix for H1 basis functions. + + Computes the six surface integrals + + S_i'_{ijk,lmn} = int_{partial Omega_i'} Lambda^0_{ijk} Lambda^0_{lmn} sqrt(g) |DF^-T n_hat_i| dS + + and adds them together into a single StencilMatrix S such that + + I = psi^T S alpha + + for any discrete test function psi_h and spline function alpha_h in V^0_h. + + Parameters + ---------- + mass_ops : WeightedMassOperators + Mass operators object, contains geometry and derham. + """ + + def __init__(self, mass_ops: WeightedMassOperators, active_faces: list[bool]): + self._mass_ops = mass_ops + self._derham = mass_ops.derham + self._domain_obj = mass_ops.domain + + # H1 space info + self._space = self._derham.fem_spaces["0"] + self._space_key = "0" + + # 3D quadrature grid info for H1 space + self._quad_grid_pts = self._derham.spline_attributes[self._space_key].quad_grid_pts + self._spans_l = self._derham.spline_attributes[self._space_key].quad_grid_spans + self._wts_l = self._derham.spline_attributes[self._space_key].quad_grid_wts + self._bases_l = self._derham.spline_attributes[self._space_key].quad_grid_bases + self._tensor_fem_spaces = self._derham.spline_attributes[self._space_key].tensor_spaces + self._nbasis = self._derham.spline_attributes[self._space_key].nbasis + + # boundary and extraction operators + self._V_extraction_op = self._derham.extraction_ops[self._space_key] + self._W_extraction_op = self._derham.extraction_ops[self._space_key] + self._V_boundary_op = self._derham.boundary_ops[self._space_key] + self._W_boundary_op = self._derham.boundary_ops[self._space_key] + + self._V_extraction_op_T = self._V_extraction_op.T + self._W_extraction_op_T = self._W_extraction_op.T + self._V_boundary_op_T = self._V_boundary_op.T + self._W_boundary_op_T = self._W_boundary_op.T + + # initialize StencilMatrix + fem_space = self._tensor_fem_spaces[0] + self._mat = StencilMatrix( + fem_space.coeff_space, + fem_space.coeff_space, + backend=PSYDAC_BACKEND_GPYCCEL, + precompiled=True, + ) + + # build composite operator B * E * M * E^T * B^T + self._M = self._W_extraction_op @ self._mat @ self._V_extraction_op_T + self._M0 = self._W_boundary_op @ self._M @ self._V_boundary_op_T + + # set domain and codomain + self._domain = self._M0.domain + self._codomain = self._M0.codomain + self._dtype = fem_space.coeff_space.dtype + + # allocate temporaries + self._temp_WB = self._W_boundary_op.domain.zeros() + self._temp_WE = self._W_extraction_op.domain.zeros() + self._temp_VB = self._V_boundary_op.domain.zeros() + self._temp_VE = self._V_extraction_op.domain.zeros() + self._temp_mat = self._mat.domain.zeros() + + # determine which faces to integrate over based on bcs + self._active_faces = active_faces + + # for each active face, extract surface quadrature grid and geometric weights + self._surface_quad_grid_meshes = [] + self._surface_geom_weights = [] + self._surface_spans = [] + self._surface_wts = [] + self._surface_bases = [] + + for face_idx in range(6): + if not self._active_faces[face_idx]: + self._surface_quad_grid_meshes.append(None) + self._surface_geom_weights.append(None) + self._surface_spans.append(None) + self._surface_wts.append(None) + self._surface_bases.append(None) + continue + + normal_dir = face_idx % 3 + surf_dirs = [d for d in range(3) if d != normal_dir] + fixed_val = 0.0 if face_idx < 3 else 1.0 + + surf_pts = [self._quad_grid_pts[0][d].flatten() for d in surf_dirs] + self._surface_quad_grid_meshes.append(xp.meshgrid(*surf_pts, indexing="ij")) + + surf_pts_1d = [self._quad_grid_pts[0][d].flatten() for d in surf_dirs] + e_1d = [None, None, None] + e_1d[surf_dirs[0]] = surf_pts_1d[0] + e_1d[surf_dirs[1]] = surf_pts_1d[1] + e_1d[normal_dir] = xp.array([fixed_val]) + + sqrt_g = xp.abs(self._domain_obj.jacobian_det(*e_1d)) + DFinv = self._domain_obj.jacobian_inv(*e_1d, change_out_order=True) + DFinv_n = DFinv[..., normal_dir, :] + norm_DFinv_n = xp.sqrt(xp.sum(DFinv_n**2, axis=-1)) + + surface_geom_weights = xp.squeeze(sqrt_g * norm_DFinv_n) + self._surface_geom_weights.append(surface_geom_weights) + + self._surface_spans.append([self._spans_l[0][d] for d in surf_dirs]) + self._surface_wts.append([self._wts_l[0][d] for d in surf_dirs]) + self._surface_bases.append([self._bases_l[0][d] for d in surf_dirs]) + + # load assembly kernel + self._assembly_kernel = Pyccelkernel(mass_kernels.surface_kernel_3d_mat_h1) + + self.assemble() + + @property + def domain(self): + return self._domain + + @property + def codomain(self): + return self._codomain + + @property + def dtype(self): + return self._dtype + + def _assemble_face( + self, + face_idx: int, + mat: StencilMatrix, + ): + """ + Assembles the contribution of a single face to the boundary mass matrix. + + Parameters + ---------- + face_idx : int + Index of the face (0 to 5). + + mat : StencilMatrix + Output matrix to accumulate into. + """ + normal_dir = face_idx % 3 + fem_space = self._tensor_fem_spaces[0] + starts = [int(start) for start in fem_space.coeff_space.starts] + ends = [int(end) for end in fem_space.coeff_space.ends] + pads = fem_space.coeff_space.pads + + boundary_index = 0 if face_idx < 3 else self._nbasis[0][normal_dir] - 1 + + logger.debug(f"{normal_dir=}, {face_idx=} {boundary_index=}, {starts=}, {ends=}, {pads=}") + + # only assemble if current rank is a true boundary (not an interior partition boundary) + if starts[normal_dir] == boundary_index or ends[normal_dir] == boundary_index: + logger.debug("Assembling face", face_idx) + self._assembly_kernel( + *self._surface_spans[face_idx], + *fem_space.degree, + *fem_space.degree, + *starts, + *pads, + *self._surface_wts[face_idx], + *self._surface_bases[face_idx], + *self._surface_bases[face_idx], + boundary_index, + normal_dir, + self._surface_geom_weights[face_idx], + mat._data, + ) + + def assemble( + self, + clear: bool = True, + ): + """ + Assembles the boundary mass matrix. + + Parameters + ---------- + clear : bool, optional + Whether to zero the matrix before assembly. + """ + if clear: + self._mat._data[:] = 0.0 + + for face_idx in range(6): + if not self._active_faces[face_idx]: + continue + self._assemble_face(face_idx, self._mat) + + self._mat.exchange_assembly_data() + self._mat.update_ghost_regions() + + def dot(self, v, out=None, apply_bc=True): + """ + Applies the boundary mass matrix to a vector. + + Parameters + ---------- + v : StencilVector + Input vector (spline coefficients of alpha_h). + + out : StencilVector, optional + Output vector. If None, a new zero vector is created. + + apply_bc : bool + Whether to apply boundary operators. + + Returns + ------- + out : StencilVector + The result S * v. + """ + if out is None: + out = self.codomain.zeros() + + if apply_bc: + self._V_boundary_op_T.dot(v, out=self._temp_VB) + self._V_extraction_op_T.dot(self._temp_VB, out=self._temp_mat) + self._mat.dot(self._temp_mat, out=self._temp_WE) + self._W_extraction_op.dot(self._temp_WE, out=self._temp_WB) + self._W_boundary_op.dot(self._temp_WB, out=out) + else: + self._V_extraction_op_T.dot(v, out=self._temp_mat) + self._mat.dot(self._temp_mat, out=self._temp_WE) + self._W_extraction_op.dot(self._temp_WE, out=out) + + return out + + def transpose(self, conjugate=False): + """ + Returns self since the boundary mass matrix is symmetric. + """ + return self + + def __call__( + self, + clear: bool = True, + ): + """ + Assembles the boundary mass matrix. + + Parameters + ---------- + clear : bool, optional + Whether to zero the matrix before assembly. + """ + self.assemble(clear=clear) + return self + + def toarray(self): + return self._M0.toarray() + + def tosparse(self): + return self._M0.tosparse() + + +class BoundaryMassOperatorHCurl(LinOpWithTransp): + """ + Assembles the boundary mass matrix for H(curl) basis functions. + + Computes the surface integrals + + W^{mu,nu}_{ijk,lmn} = int_{partial Omega} hat_Lambda^1_{mu,ijk} hat_R_n^{mu,nu} hat_Lambda^1_{nu,lmn} dS + + where hat_R_n is the pullback of [n]_x to logical coordinates. + + The result is a 3x3 BlockLinearOperator where diagonal blocks are zero + (skew-symmetry of [n]_x) and off-diagonal blocks are assembled via + surface_kernel_3d_mat_h1. + + Parameters + ---------- + mass_ops : WeightedMassOperators + Mass operators object, contains geometry and derham. + active_faces : list[bool] + Which of the six faces to integrate over. + """ + + def __init__( + self, + mass_ops: WeightedMassOperators, + active_faces: list[bool], + ): + self._mass_ops = mass_ops + self._derham = mass_ops.derham + self._domain_obj = mass_ops.domain + self._active_faces = active_faces + + self._space_key = "1" + self._space = self._derham.fem_spaces[self._space_key] + self._quad_grid_pts = self._derham.spline_attributes[self._space_key].quad_grid_pts + self._spans_l = self._derham.spline_attributes[self._space_key].quad_grid_spans + self._wts_l = self._derham.spline_attributes[self._space_key].quad_grid_wts + self._bases_l = self._derham.spline_attributes[self._space_key].quad_grid_bases + self._tensor_fem_spaces = self._derham.spline_attributes[self._space_key].tensor_spaces + self._nbasis = self._derham.spline_attributes[self._space_key].nbasis + + self._V_extraction_op = self._derham.extraction_ops[self._space_key] + self._W_extraction_op = self._derham.extraction_ops[self._space_key] + self._V_boundary_op = self._derham.boundary_ops[self._space_key] + self._W_boundary_op = self._derham.boundary_ops[self._space_key] + + self._V_extraction_op_T = self._V_extraction_op.T + self._W_extraction_op_T = self._W_extraction_op.T + self._V_boundary_op_T = self._V_boundary_op.T + self._W_boundary_op_T = self._W_boundary_op.T + + V = self._space + W = self._space + + blocks = [ + [ + StencilMatrix( + Vs.coeff_space, + Ws.coeff_space, + backend=PSYDAC_BACKEND_GPYCCEL, + precompiled=True, + ) + if i != j + else None + for j, Vs in enumerate(V.spaces) + ] + for i, Ws in enumerate(W.spaces) + ] + + self._mat = BlockLinearOperator( + V.coeff_space, + W.coeff_space, + blocks=blocks, + ) + + self._M = self._W_extraction_op @ self._mat @ self._V_extraction_op_T + self._M0 = self._W_boundary_op @ self._M @ self._V_boundary_op_T + + self._domain = self._M0.domain + self._codomain = self._M0.codomain + self._dtype = self._tensor_fem_spaces[0].coeff_space.dtype + + # allocate temporaries + self._temp_VB = self._V_boundary_op.domain.zeros() + self._temp_VE = self._V_extraction_op.domain.zeros() + self._temp_WB = self._W_boundary_op.domain.zeros() + self._temp_WE = self._W_extraction_op.domain.zeros() + self._temp_mat = self._mat.domain.zeros() + + # for each active face extract surface data + self._surface_R_n = [] + self._surface_spans = [] + self._surface_wts = [] + self._surface_bases = [] + + for face_idx in range(6): + if not self._active_faces[face_idx]: + self._surface_R_n.append(None) + self._surface_spans.append(None) + self._surface_wts.append(None) + self._surface_bases.append(None) + continue + + normal_dir = face_idx % 3 + surf_dirs = [d for d in range(3) if d != normal_dir] + + sign = 1.0 if face_idx < 3 else -1.0 + n_hat = xp.zeros(3) + n_hat[normal_dir] = sign + + # constant skew-symmetric cross-product matrix R_n such that R_n v = n_hat x v + R_n_const = xp.zeros((3, 3)) + R_n_const[0, 1] = -n_hat[2] + R_n_const[0, 2] = n_hat[1] + R_n_const[1, 0] = n_hat[2] + R_n_const[1, 2] = -n_hat[0] + R_n_const[2, 0] = -n_hat[1] + R_n_const[2, 1] = n_hat[0] + + # store R_n per component mu on its own quadrature grid shape + surface_R_n_per_mu = [None, None, None] + for mu in surf_dirs: + nq1 = self._spans_l[mu][surf_dirs[0]].size * self._wts_l[mu][surf_dirs[0]].shape[1] + nq2 = self._spans_l[mu][surf_dirs[1]].size * self._wts_l[mu][surf_dirs[1]].shape[1] + R_n_mu = xp.zeros((nq1, nq2, 3, 3)) + R_n_mu[..., :, :] = R_n_const + surface_R_n_per_mu[mu] = R_n_mu + + self._surface_R_n.append(surface_R_n_per_mu) + surface_spans_per_mu = [None, None, None] + surface_wts_per_mu = [None, None, None] + surface_bases_per_mu = [None, None, None] + + for mu in surf_dirs: + surface_spans_per_mu[mu] = [self._spans_l[mu][d] for d in surf_dirs] + surface_wts_per_mu[mu] = [self._wts_l[mu][d] for d in surf_dirs] + surface_bases_per_mu[mu] = [self._bases_l[mu][d] for d in surf_dirs] + + self._surface_spans.append(surface_spans_per_mu) + self._surface_wts.append(surface_wts_per_mu) + self._surface_bases.append(surface_bases_per_mu) + + self._assembly_kernel = Pyccelkernel(mass_kernels.surface_kernel_3d_mat_h1) + + self.assemble() + + @property + def domain(self): + return self._domain + + @property + def codomain(self): + return self._codomain + + @property + def dtype(self): + return self._dtype + + def _assemble_face( + self, + face_idx: int, + mat: BlockLinearOperator, + ): + normal_dir = face_idx % 3 + surf_dirs = [d for d in range(3) if d != normal_dir] + + mu, nu = surf_dirs[0], surf_dirs[1] + + fem_space_mu = self._tensor_fem_spaces[mu] + fem_space_nu = self._tensor_fem_spaces[nu] + + starts_mu = [int(s) for s in fem_space_mu.coeff_space.starts] + ends_mu = [int(e) for e in fem_space_mu.coeff_space.ends] + pads_mu = fem_space_mu.coeff_space.pads + + starts_nu = [int(s) for s in fem_space_nu.coeff_space.starts] + ends_nu = [int(e) for e in fem_space_nu.coeff_space.ends] + pads_nu = fem_space_nu.coeff_space.pads + + boundary_index_mu = 0 if face_idx < 3 else self._nbasis[mu][normal_dir] - 1 + boundary_index_nu = 0 if face_idx < 3 else self._nbasis[nu][normal_dir] - 1 + + logger.debug(f"{normal_dir=}, {face_idx=} {boundary_index_mu=}, {starts_mu=}, {ends_mu=}, {pads_mu=}") + logger.debug(f"{normal_dir=}, {face_idx=} {boundary_index_nu=}, {starts_nu=}, {ends_nu=}, {pads_nu=}") + + mat_fun_mu_nu = self._surface_R_n[face_idx][mu][..., mu, nu] + mat_fun_nu_mu = self._surface_R_n[face_idx][nu][..., nu, mu] + + if starts_mu[normal_dir] == boundary_index_mu or ends_mu[normal_dir] == boundary_index_mu: + logger.debug(f"Assembling face {face_idx} for block ({mu},{nu})") + self._assembly_kernel( + *self._surface_spans[face_idx][mu], + *fem_space_mu.degree, + *fem_space_nu.degree, + *starts_mu, + *pads_mu, + *self._surface_wts[face_idx][mu], + *self._surface_bases[face_idx][mu], + *self._surface_bases[face_idx][nu], + boundary_index_mu, + normal_dir, + mat_fun_mu_nu, + mat.blocks[mu][nu]._data, + ) + + if starts_nu[normal_dir] == boundary_index_nu or ends_nu[normal_dir] == boundary_index_nu: + logger.debug(f"Assembling face {face_idx} for block ({nu},{mu})") + self._assembly_kernel( + *self._surface_spans[face_idx][nu], + *fem_space_nu.degree, + *fem_space_mu.degree, + *starts_nu, + *pads_nu, + *self._surface_wts[face_idx][nu], + *self._surface_bases[face_idx][nu], + *self._surface_bases[face_idx][mu], + boundary_index_nu, + normal_dir, + mat_fun_nu_mu, + mat.blocks[nu][mu]._data, + ) + + def assemble(self, clear: bool = True): + """Assembles the H(curl) boundary mass matrix.""" + if clear: + for mu in range(3): + for nu in range(3): + if mu != nu: + self._mat.blocks[mu][nu]._data[:] = 0.0 + + for face_idx in range(6): + if not self._active_faces[face_idx]: + continue + self._assemble_face(face_idx, self._mat) + + for mu in range(3): + for nu in range(3): + if mu != nu: + self._mat.blocks[mu][nu].exchange_assembly_data() + self._mat.blocks[mu][nu].update_ghost_regions() + + def dot(self, v, out=None, apply_bc=True): + """Applies the H(curl) boundary mass matrix to a BlockVector.""" + if out is None: + out = self.codomain.zeros() + + if apply_bc: + self._V_boundary_op_T.dot(v, out=self._temp_VB) + self._V_extraction_op_T.dot(self._temp_VB, out=self._temp_mat) + self._mat.dot(self._temp_mat, out=self._temp_WE) + self._W_extraction_op.dot(self._temp_WE, out=self._temp_WB) + self._W_boundary_op.dot(self._temp_WB, out=out) + else: + self._V_extraction_op_T.dot(v, out=self._temp_mat) + self._mat.dot(self._temp_mat, out=self._temp_WE) + self._W_extraction_op.dot(self._temp_WE, out=out) + + return out + + def transpose(self, conjugate=False): + return self + + def toarray(self): + return self._M0.toarray() + + def tosparse(self): + return self._M0.tosparse() diff --git a/src/struphy/feec/mass_kernels.py b/src/struphy/feec/mass_kernels.py index 7b4f09720..8a7244918 100644 --- a/src/struphy/feec/mass_kernels.py +++ b/src/struphy/feec/mass_kernels.py @@ -21,13 +21,39 @@ def kernel_1d_mat( data: "float[:,:]", ): """ - Performs the integration of Lambda_i * mat_fun(eta1) * Lambda_l for the basis functions (i, l) available on the calling process. + 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): @@ -58,9 +84,29 @@ def kernel_1d_vec( data: "float[:]", ): """ - Performs the integration of Lambda_i * mat_fun(eta1) for the basis functions (i) available on the calling process. + 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 @@ -93,9 +139,28 @@ def kernel_1d_eval( values: "float[:]", ): """ - Evaluates sum_i [ coeffs_i * Lambda_i(quad_eta1) ] for all quadrature points on the calling process. + 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 @@ -140,9 +205,35 @@ def kernel_2d_mat( data: "float[:,:,:,:]", ): """ - Performs the integration of Lambda_ij * mat_fun(eta1, eta2) * Lambda_lm for the basis functions (ij, lm) available on the calling process. + 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 @@ -195,9 +286,31 @@ def kernel_2d_vec( data: "float[:,:]", ): """ - Performs the integration of Lambda_ij * mat_fun(eta1, eta2) for the basis functions (ij) available on the calling process. + 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 @@ -244,9 +357,29 @@ def kernel_2d_eval( values: "float[:,:]", ): """ - Evaluates sum_ij [ coeffs_ij * Lambda_ij(quad_eta1, quad_eta2) ] for all quadrature points on the calling process. + 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 @@ -310,9 +443,35 @@ def kernel_3d_mat( data: "float[:,:,:,:,:,:]", ): """ - Performs the integration of Lambda_ijk * mat_fun(eta1, eta2, eta3) * Lambda_lmn for the basis functions (ijk, lmn) available on the calling process. + 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 @@ -421,9 +580,31 @@ def kernel_3d_vec( data: "float[:,:,:]", ): """ - Performs the integration of Lambda_ijk * mat_fun(eta1, eta2, eta3) for the basis functions (ijk) available on the calling process. + 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 @@ -489,9 +670,29 @@ def kernel_3d_eval( values: "float[:,:,:]", ): """ - Evaluates sum_ijk [ coeffs_ijk * Lambda_ijk(quad_eta1, quad_eta2, quad_eta3) ] for all quadrature points on the calling process. + 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 @@ -570,10 +771,45 @@ def kernel_3d_matrixfree( data_in: "float[:,:,:]", ): """ - Performs the integration of Lambda_ijk * mat_fun(eta1, eta2, eta3) * f(eta1, eta2, eta3) for the basis functions (ijk) available on the calling process, + 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 (attention: data is NOT set to zero first, but the results are added to data). + 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 @@ -691,6 +927,30 @@ def kernel_3d_diag( 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 (unused for data, which is a + StencilDiagonalMatrix and therefore has no padding, but kept for a uniform kernel signature). + 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 StencilDiagonalMatrix to store the results. Periodic wrap-around (index -= nb) is applied + when a local index runs beyond the array bounds, since there are no ghost regions on this matrix type. """ ne1 = spans1.size @@ -766,3 +1026,309 @@ def kernel_3d_diag( # No padding on StencilDiagonalMatrix 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_h1( + 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[:,:,:,:,:,:]", +): + """ + Computes a boundary (surface) H1 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] + + starts = [starts0, starts1, starts2] + pads = [pads0, pads1, pads2] + pi = [pi0, pi1, pi2] + pj = [pj0, pj1, pj2] + + surf_dirs = [d for d in range(3) if d != normal_dir] + + pi_s1 = pi[surf_dirs[0]] + pi_s2 = pi[surf_dirs[1]] + + pj_s1 = pj[surf_dirs[0]] + pj_s2 = pj[surf_dirs[1]] + + starts_n = starts[normal_dir] + starts_s1 = starts[surf_dirs[0]] + starts_s2 = starts[surf_dirs[1]] + + pads_n = pads[normal_dir] + pads_s1 = pads[surf_dirs[0]] + pads_s2 = pads[surf_dirs[1]] + + 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 + + +def surface_kernel_3d_mat_hcurl( + spans1: "int[:]", + spans2: "int[:]", + pi0: int, + pi1: int, + pi2: int, + qi0: int, + qi1: int, + qi2: 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, + n_cross_weight: "float[:,:]", + data: "float[:,:,:,:,:,:]", +): + """ + Computes a boundary (surface) H(curl) mass matrix: the integration of (n x Lambda_i) * mat_fun(eta_s1, eta_s2) + * (n x Lambda_j), weighted by n_cross_weight, over the boundary surface at the fixed global index + boundary_index, for the codomain ("i") and domain ("j") tangential-trace basis functions available on the + calling process. Mirrors surface_kernel_3d_mat_h1, but for H(curl) spaces where only the tangential + components of the basis functions couple through the cross product with the surface normal n. + + 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; 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 Cartesian axes 0, 1 and 2. + qi0, qi1, qi2 : int + Degree of the domain basis functions along Cartesian axes 0, 1 and 2. + starts0, starts1, starts2 : int + Starting index on the current rank along Cartesian axes 0, 1 and 2. + pads0, pads1, pads2 : int + Padding (=spline degree) for ghost regions in data, along Cartesian 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 in the normal direction at which the boundary surface is located. + n_cross_weight : "float[:,:]" + Weight from the cross product with the surface normal n (and mat_fun) evaluated at surface quadrature + points, indexed like mat_fun in surface_kernel_3d_mat_h1. + data : "float[:,:,:,:,:,:]" + _data array of StencilMatrix to store the results, following the same storage convention as in + surface_kernel_3d_mat_h1. + """ + pass diff --git a/src/struphy/feec/psydac_derham.py b/src/struphy/feec/psydac_derham.py index 82bc49377..e6ee50e60 100644 --- a/src/struphy/feec/psydac_derham.py +++ b/src/struphy/feec/psydac_derham.py @@ -507,22 +507,30 @@ def proj_loc_grid_wts(self) -> tuple[tuple[xp.ndarray]] | None: @property def quad_grid_pts(self) -> tuple[tuple[xp.ndarray]]: - """Tuple of quadrature grid points in each direction for each component of the vector space.""" + """Tuple of quadrature grid points in each direction for each component of the vector space. + The indexing is [global element, local quadrature point]. + The length of the first dimension is the number of elements (cells).""" return self._quad_grid_pts @property def quad_grid_wts(self) -> tuple[tuple[xp.ndarray]]: - """Tuple of quadrature grid weights in each direction for each component of the vector space.""" + """Tuple of quadrature grid weights in each direction for each component of the vector space. + The indexing is [global element, local quadrature point]. + The length of the first dimension is the number of elements (cells).""" return self._quad_grid_wts @property def quad_grid_spans(self) -> tuple[tuple[xp.ndarray]]: - """Tuple of quadrature grid basis function spans in each direction for each component of the vector space.""" + """Tuple of quadrature grid basis function spans in each direction for each component of the vector space. + 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).""" return self._quad_grid_spans @property def quad_grid_bases(self) -> tuple[tuple[xp.ndarray]]: - """Tuple of quadrature grid basis function values in each direction for each component of the vector space.""" + """Tuple of quadrature grid basis function values in each direction for each component of the vector space. + Indexing is [global element, basis function, derivative, local quadrature point]. + The length of the first dimension is the number of elements (cells).""" return self._quad_grid_bases diff --git a/src/struphy/feec/tests/test_boundary_integrals.py b/src/struphy/feec/tests/test_boundary_integrals.py new file mode 100644 index 000000000..e79294a18 --- /dev/null +++ b/src/struphy/feec/tests/test_boundary_integrals.py @@ -0,0 +1,371 @@ +import logging +from typing import Callable + +import cunumpy as xp +import pytest +from feectools.ddm.mpi import mpi as MPI + +from struphy import domains +from struphy.feec.boundary_mass import BoundaryIntegralOperators +from struphy.feec.mass import L2Projector, WeightedMassOperators +from struphy.feec.psydac_derham import Derham +from struphy.io.options import DerhamOptions +from struphy.topology.grids import TensorProductGrid + +logger = logging.getLogger("struphy") + + +@pytest.mark.parametrize( + "num_elements", + [[8, 9, 10]], +) +@pytest.mark.parametrize("degree", [[1, 2, 3]]) +@pytest.mark.parametrize( + "bcs", + [ + (None, None, None), + (("free", "free"), None, None), + (None, ("free", "free"), None), + (None, None, ("free", "free")), + (("free", "free"), ("free", "free"), None), + (("free", "free"), ("free", "free"), ("free", "free")), + ], +) +def test_boundary_mass_unit_cube_constant(num_elements, degree, bcs): + """ + Tests the boundary mass operator for alpha = 1 on the unit cube. + """ + comm = MPI.COMM_WORLD + + grid = TensorProductGrid(num_elements=num_elements) + derham_opts = DerhamOptions(degree=degree, bcs=bcs) + derham = Derham(grid, derham_opts, comm=comm) + + domain = domains.Cuboid(l1=0.0, r1=1.0, l2=0.0, r2=1.0, l3=0.0, r3=1.0) + mass_ops = WeightedMassOperators(derham, domain) + + face_value = 2.1 + alpha = lambda e1, e2, e3: xp.ones_like(e1) * face_value + + num_faces = 0 + for face_tuple in bcs: + if face_tuple is None: + continue + if face_tuple[0] == "free": + num_faces += 1 + if face_tuple[1] == "free": + num_faces += 1 + + exact = num_faces * face_value + + P = L2Projector("H1", mass_ops) + alpha_h = P(alpha) + + bnd_ops = BoundaryIntegralOperators(mass_ops) + v = bnd_ops.S0.dot(alpha_h) + arr = v.toarray() + + if comm is None: + coeffs = arr + else: + coeffs = xp.zeros_like(arr) + comm.Allreduce(arr, coeffs, op=MPI.SUM) + + numerical = xp.sum(coeffs) + + logger.info(f"numerical = {numerical}, exact = {exact}, error = {xp.abs(numerical - exact)}") + + assert xp.abs(numerical - exact) < 1e-3 + + +@pytest.mark.parametrize("num_elements", [[8, 9, 10]]) +@pytest.mark.parametrize("degree", [[1, 2, 3]]) +@pytest.mark.parametrize( + "bcs", + [ + (("dirichlet", "free"), ("free", "free"), ("free", "free")), + (("free", "dirichlet"), ("free", "free"), ("free", "free")), + (("dirichlet", "dirichlet"), ("free", "free"), ("free", "free")), + ], +) +def test_boundary_mass_unit_cube_nonconstant(num_elements, degree, bcs): + """ + Tests the boundary mass operator for alpha = eta1 + eta2 + eta3 on the unit cube. + """ + comm = MPI.COMM_WORLD + + grid = TensorProductGrid(num_elements=num_elements) + derham_opts = DerhamOptions(degree=degree, bcs=bcs) + derham = Derham(grid, derham_opts, comm=comm) + + domain = domains.Cuboid(l1=0.0, r1=1.0, l2=0.0, r2=1.0, l3=0.0, r3=1.0) + mass_ops = WeightedMassOperators(derham, domain) + + if bcs[0] == ("dirichlet", "free"): + alpha = lambda e1, e2, e3: e1 + 0 * e2 + 0 * e3 + exact = 3.0 + elif bcs[0] == ("free", "dirichlet"): + alpha = lambda e1, e2, e3: 1.0 - e1 + 0 * e2 + 0 * e3 + exact = 3.0 + else: + assert bcs[0] == ("dirichlet", "dirichlet") + alpha = lambda e1, e2, e3: e1 * (1.0 - e1) + 0 * e2 + 0 * e3 + exact = 2.0 / 3.0 + + P = L2Projector("H1", mass_ops) + alpha_h = P(alpha, apply_bc=True) + + bnd_ops = BoundaryIntegralOperators(mass_ops) + v = bnd_ops.S0.dot(alpha_h) + arr = v.toarray() + + if comm is None: + coeffs = arr + else: + coeffs = xp.zeros_like(arr) + comm.Allreduce(arr, coeffs, op=MPI.SUM) + + numerical = xp.sum(coeffs) + + logger.info(f"numerical = {numerical}, exact = {exact}, error = {xp.abs(numerical - exact)}") + + assert xp.abs(numerical - exact) < 2e-2 + + +@pytest.mark.parametrize("num_elements", [[8, 9, 10]]) +@pytest.mark.parametrize("degree", [[1, 2, 3]]) +@pytest.mark.parametrize("bcs", [(("free", "free"), ("free", "free"), ("free", "free"))]) +def test_boundary_mass_cuboid_nontrivial(num_elements, degree, bcs): + """ + Tests the boundary mass operator for alpha = eta1 + eta2 + eta3 + on a non-unit cuboid [-1,1] x [-1,3] x [0,3]. + """ + comm = MPI.COMM_WORLD + + grid = TensorProductGrid(num_elements=num_elements) + derham_opts = DerhamOptions(degree=degree, bcs=bcs) + derham = Derham(grid, derham_opts, comm=comm) + + domain = domains.Cuboid(l1=-1.0, r1=1.0, l2=-1.0, r2=3.0, l3=0.0, r3=3.0) + mass_ops = WeightedMassOperators(derham, domain) + + alpha = lambda e1, e2, e3: e1 + e2 + e3 + exact = 78.0 + + P = L2Projector("H1", mass_ops) + alpha_h = P(alpha) + + bnd_ops = BoundaryIntegralOperators(mass_ops) + v = bnd_ops.S0.dot(alpha_h) + arr = v.toarray() + + if comm is None: + coeffs = arr + else: + coeffs = xp.zeros_like(arr) + comm.Allreduce(arr, coeffs, op=MPI.SUM) + + numerical = xp.sum(coeffs) + + logger.info(f"numerical = {numerical}, exact = {exact}, error = {xp.abs(numerical - exact)}") + + assert xp.abs(numerical - exact) < 1e-3 + + +@pytest.mark.parametrize("num_elements", [[8, 9, 10]]) +@pytest.mark.parametrize("degree", [[1, 2, 3]]) +@pytest.mark.parametrize("bcs", [(("free", "free"), None, ("free", "free"))]) +def test_boundary_mass_hollow_cylinder_nonconstant(num_elements, degree, bcs): + """ + Tests the boundary mass operator for alpha = exp(eta3) on a HollowCylinder. + """ + import math + + comm = MPI.COMM_WORLD + + grid = TensorProductGrid(num_elements=num_elements) + derham_opts = DerhamOptions(degree=degree, bcs=bcs) + derham = Derham(grid, derham_opts, comm=comm) + + a1 = 0.2 + a2 = 1.0 + Lz = 4.0 + + domain = domains.HollowCylinder(a1=a1, a2=a2, Lz=Lz) + mass_ops = WeightedMassOperators(derham, domain) + + alpha = lambda e1, e2, e3: xp.exp(e3) + e = math.e + exact = xp.pi * (2 * a1 * Lz * (e - 1) + 2 * a2 * Lz * (e - 1) + (a2**2 - a1**2) * (1 + e)) + + P = L2Projector("H1", mass_ops) + alpha_h = P(alpha) + + bnd_ops = BoundaryIntegralOperators(mass_ops) + v = bnd_ops.S0.dot(alpha_h) + arr = v.toarray() + + if comm is None: + coeffs = arr + else: + coeffs = xp.zeros_like(arr) + comm.Allreduce(arr, coeffs, op=MPI.SUM) + + numerical = xp.sum(coeffs) + + logger.info(f"numerical = {numerical}, exact = {exact}, error = {xp.abs(numerical - exact)}") + + assert xp.abs(numerical - exact) < 1e-2 + + +@pytest.mark.parametrize("num_elements", [[10, 10, 10]]) +@pytest.mark.parametrize("degree", [[2, 2, 2]]) +@pytest.mark.parametrize("bcs", [(("free", "free"), ("free", "free"), ("free", "free"))]) +@pytest.mark.parametrize( + "active_faces, u_idx, v_idx, exact", + [ + ([True, False, False, False, False, False], 1, 2, 1.0), + ([False, True, False, False, False, False], 2, 0, 1.0), + ([False, False, True, False, False, False], 0, 1, 1.0), + ([False, False, False, True, False, False], 1, 2, -1.0), + ([False, False, False, False, True, False], 2, 0, -1.0), + ([False, False, False, False, False, True], 0, 1, -1.0), + ], +) +def test_boundary_mass_hcurl_per_face(num_elements, degree, bcs, active_faces, u_idx, v_idx, exact): + comm = MPI.COMM_WORLD + + grid = TensorProductGrid(num_elements=num_elements) + derham_opts = DerhamOptions(degree=degree, bcs=bcs) + derham = Derham(grid, derham_opts, comm=comm) + + domain = domains.Cuboid(l1=0.0, r1=1.0, l2=0.0, r2=1.0, l3=0.0, r3=1.0) + mass_ops = WeightedMassOperators(derham, domain) + + u_funs = [ + lambda e1, e2, e3: xp.ones_like(e1) if 0 == u_idx else xp.zeros_like(e1), + lambda e1, e2, e3: xp.ones_like(e1) if 1 == u_idx else xp.zeros_like(e1), + lambda e1, e2, e3: xp.ones_like(e1) if 2 == u_idx else xp.zeros_like(e1), + ] + + v_funs = [ + lambda e1, e2, e3: xp.ones_like(e1) if 0 == v_idx else xp.zeros_like(e1), + lambda e1, e2, e3: xp.ones_like(e1) if 1 == v_idx else xp.zeros_like(e1), + lambda e1, e2, e3: xp.ones_like(e1) if 2 == v_idx else xp.zeros_like(e1), + ] + + P = L2Projector("Hcurl", mass_ops) + u_h = P(u_funs) + v_h = P(v_funs) + + bnd_ops = BoundaryIntegralOperators(mass_ops, active_faces=active_faces) + numerical = bnd_ops.S1.dot_inner(u_h, v_h) + + logger.info(f"numerical = {numerical}, exact = {exact}, error = {xp.abs(numerical - exact)}") + + assert xp.abs(numerical - exact) < 1e-1 + + +@pytest.mark.parametrize("num_elements", [[10, 10, 10]]) +@pytest.mark.parametrize("degree", [[2, 2, 2]]) +@pytest.mark.parametrize("bcs", [(("free", "free"), ("free", "free"), ("free", "free"))]) +@pytest.mark.parametrize( + "active_faces, u_idx, v_idx, exact", + [ + ([True, False, False, False, False, False], 1, 2, 12.0), + ([False, True, False, False, False, False], 2, 0, 6.0), + ([False, False, True, False, False, False], 0, 1, 8.0), + ([False, False, False, True, False, False], 1, 2, -12.0), + ([False, False, False, False, True, False], 2, 0, -6.0), + ([False, False, False, False, False, True], 0, 1, -8.0), + ], +) +def test_boundary_mass_hcurl_cuboid_nontrivial(num_elements, degree, bcs, active_faces, u_idx, v_idx, exact): + """ + Tests the H(curl) boundary mass operator on a non-unit cuboid [-1,1] x [-1,3] x [0,3] + with constant unit vector fields u = e_{u_idx} and v = e_{v_idx}. + """ + comm = MPI.COMM_WORLD + + grid = TensorProductGrid(num_elements=num_elements) + derham_opts = DerhamOptions(degree=degree, bcs=bcs) + derham = Derham(grid, derham_opts, comm=comm) + + domain = domains.Cuboid(l1=-1.0, r1=1.0, l2=-1.0, r2=3.0, l3=0.0, r3=3.0) + mass_ops = WeightedMassOperators(derham, domain) + + def make_pulled(domain, idx): + phys_funs = [ + lambda x, y, z: xp.ones_like(x) if 0 == idx else xp.zeros_like(x), + lambda x, y, z: xp.ones_like(x) if 1 == idx else xp.zeros_like(x), + lambda x, y, z: xp.ones_like(x) if 2 == idx else xp.zeros_like(x), + ] + + def pulled(*etas): + return domain.pull(phys_funs, *etas, kind="1") + + return [ + lambda *etas, p=pulled: p(*etas)[0], + lambda *etas, p=pulled: p(*etas)[1], + lambda *etas, p=pulled: p(*etas)[2], + ] + + P = L2Projector("Hcurl", mass_ops) + u_h = P(make_pulled(domain, u_idx)) + v_h = P(make_pulled(domain, v_idx)) + + bnd_ops = BoundaryIntegralOperators(mass_ops, active_faces=active_faces) + numerical = xp.dot(v_h.toarray(), bnd_ops.S1.dot(u_h).toarray()) + + logger.info(f"numerical = {numerical}, exact = {exact}, error = {xp.abs(numerical - exact)}") + + assert xp.abs(numerical - exact) < 1 + + +if __name__ == "__main__": + from struphy import set_logging_level + + set_logging_level(logging.INFO) + + test_boundary_mass_unit_cube_constant( + [8, 9, 10], + [1, 2, 3], + (("free", "free"), ("free", "free"), ("free", "free")), + ) + + test_boundary_mass_unit_cube_nonconstant( + [8, 9, 10], + [1, 2, 3], + (("dirichlet", "free"), ("free", "free"), ("free", "free")), + ) + + test_boundary_mass_cuboid_nontrivial( + [8, 9, 10], + [1, 2, 3], + (("free", "free"), ("free", "free"), ("free", "free")), + ) + test_boundary_mass_hollow_cylinder_nonconstant( + [8, 9, 10], + [1, 2, 3], + (("free", "free"), None, ("free", "free")), + ) + + test_boundary_mass_hcurl_per_face( + [10, 10, 10], + [2, 2, 2], + (("free", "free"), ("free", "free"), ("free", "free")), + [True, False, False, False, False, False], + 1, + 2, + 1.0, + ) + + test_boundary_mass_hcurl_cuboid_nontrivial( + [10, 10, 10], + [1, 2, 3], + (("free", "free"), ("free", "free"), ("free", "free")), + [True, False, False, False, False, False], + 1, + 2, + 12.0, + )