From b2de54a5a223ace604229b110d2687d4226e93b7 Mon Sep 17 00:00:00 2001 From: "copilot-swe-agent[bot]" <198982749+Copilot@users.noreply.github.com> Date: Fri, 17 Apr 2026 15:03:04 +0000 Subject: [PATCH 1/2] Initial plan From fcf6afb36918dac9fdaf328f09570ddc00f84d8a Mon Sep 17 00:00:00 2001 From: "copilot-swe-agent[bot]" <198982749+Copilot@users.noreply.github.com> Date: Fri, 17 Apr 2026 15:19:48 +0000 Subject: [PATCH 2/2] Add interior penalty elliptic operator and convergence tests Agent-Logs-Url: https://github.com/inducer/grudge/sessions/1e832e02-7aae-4c75-9aa6-12a32585fcd5 Co-authored-by: inducer <352067+inducer@users.noreply.github.com> --- grudge/models/__init__.py | 5 ++ grudge/models/elliptic.py | 140 ++++++++++++++++++++++++++++++++++++++ test/test_elliptic.py | 63 +++++++++++++++++ 3 files changed, 208 insertions(+) create mode 100644 grudge/models/elliptic.py create mode 100644 test/test_elliptic.py diff --git a/grudge/models/__init__.py b/grudge/models/__init__.py index 833b01bc..407078cf 100644 --- a/grudge/models/__init__.py +++ b/grudge/models/__init__.py @@ -84,3 +84,8 @@ def estimate_rk4_timestep( ) return nodal_min(dcoll, "vol", local_timesteps) + + +from grudge.models.elliptic import ( + InteriorPenaltyEllipticOperator as InteriorPenaltyEllipticOperator, +) diff --git a/grudge/models/elliptic.py b/grudge/models/elliptic.py new file mode 100644 index 00000000..007cfa09 --- /dev/null +++ b/grudge/models/elliptic.py @@ -0,0 +1,140 @@ +"""Interior penalty DG discretization of elliptic operators.""" +from __future__ import annotations + + +__copyright__ = """ +Copyright (C) 2026 University of Illinois Board of Trustees +""" + +__license__ = """ +Permission is hereby granted, free of charge, to any person obtaining a copy +of this software and associated documentation files (the "Software"), to deal +in the Software without restriction, including without limitation the rights +to use, copy, modify, merge, publish, distribute, sublicense, and/or sell +copies of the Software, and to permit persons to whom the Software is +furnished to do so, subject to the following conditions: + +The above copyright notice and this permission notice shall be included in +all copies or substantial portions of the Software. + +THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR +IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, +FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE +AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER +LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, +OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN +THE SOFTWARE. +""" + +from collections.abc import Callable +from typing import Any + +import numpy as np + +from meshmode.mesh import BTAG_ALL, BTAG_NONE + +import grudge.geometry as geo +from grudge import op +from grudge.dof_desc import FACE_RESTR_ALL, as_dofdesc +from grudge.dt_utils import characteristic_lengthscales +from grudge.models import Operator + + +class InteriorPenaltyEllipticOperator(Operator): + r"""Discretize :math:`-\nabla\cdot(\kappa \nabla u)` with SIPG. + + The implementation follows the symmetric interior penalty bilinear form + + .. math:: + + a(u, v) + = \int_\Omega \nabla u \cdot \kappa \nabla v \,dx + - \sum_{e \in \mathcal{E}} + \int_e [u]\cdot\{\kappa \nabla v\} + + [v]\cdot\{\kappa \nabla u\} + - \frac{\alpha_e}{h_e}[v]\cdot[u] \,ds. + """ + + def __init__( + self, + dcoll, + *, + kappa: Any = 1.0, + penalty_factor: float = 10.0, + dirichlet_tag=BTAG_ALL, + dirichlet_bc: Any = 0.0, + neumann_tag=BTAG_NONE, + neumann_bc: Any = 0.0, + comm_tag=None): + self.dcoll = dcoll + self.kappa = kappa + self.penalty_factor = penalty_factor + self.dirichlet_tag = dirichlet_tag + self.dirichlet_bc = dirichlet_bc + self.neumann_tag = neumann_tag + self.neumann_bc = neumann_bc + self.comm_tag = comm_tag + + @staticmethod + def _is_active_boundary_tag(tag: Any) -> bool: + return tag is not None and tag != BTAG_NONE + + def _get_boundary_data(self, actx, dd, bc): + if isinstance(bc, Callable): + return bc(actx.thaw(self.dcoll.nodes(dd=dd))) + + if np.isscalar(bc): + return self.dcoll.discr_from_dd(dd).zeros(actx) + bc + + return bc + + def operator(self, u): + dcoll = self.dcoll + actx = u.array_context + assert actx is not None + + all_faces_dd = as_dofdesc(FACE_RESTR_ALL) + h = characteristic_lengthscales(actx, dcoll) + h = h * (dcoll.discr_from_dd(as_dofdesc("vol")).zeros(actx) + 1.0) + + k_grad_u = self.kappa * op.local_grad(dcoll, u) + + # {{{ divergence operator with SIPG penalty flux + q_flux = 0 + + u_tpairs = op.interior_trace_pairs(dcoll, u, comm_tag=self.comm_tag) + h_tpairs = op.interior_trace_pairs(dcoll, h, comm_tag=self.comm_tag) + q_tpairs = op.interior_trace_pairs(dcoll, k_grad_u, comm_tag=self.comm_tag) + + for u_tpair, h_tpair, q_tpair in zip( + u_tpairs, h_tpairs, q_tpairs, strict=True): + normal = geo.normal(actx, dcoll, u_tpair.dd) + penalty = self.penalty_factor / h_tpair.avg + q_flux = q_flux + op.project( + dcoll, + u_tpair.dd, + all_faces_dd, + np.dot(q_tpair.avg, normal) - penalty * (u_tpair.int - u_tpair.ext)) + + if self._is_active_boundary_tag(self.dirichlet_tag): + dir_dd = as_dofdesc(self.dirichlet_tag) + dir_normal = geo.normal(actx, dcoll, dir_dd) + dir_h = op.project(dcoll, "vol", dir_dd, h) + dir_u = op.project(dcoll, "vol", dir_dd, u) + dir_q = op.project(dcoll, "vol", dir_dd, k_grad_u) + dir_bc = self._get_boundary_data(actx, dir_dd, self.dirichlet_bc) + dir_penalty = self.penalty_factor / dir_h + + q_flux = q_flux + op.project( + dcoll, dir_dd, all_faces_dd, + np.dot(dir_q, dir_normal) - dir_penalty * (dir_u - dir_bc)) + + if self._is_active_boundary_tag(self.neumann_tag): + neu_dd = as_dofdesc(self.neumann_tag) + neu_bc = self._get_boundary_data(actx, neu_dd, self.neumann_bc) + q_flux = q_flux + op.project(dcoll, neu_dd, all_faces_dd, neu_bc) + + return op.inverse_mass( + dcoll, + op.weak_local_div(dcoll, k_grad_u) - op.face_mass(dcoll, q_flux)) + # }}} diff --git a/test/test_elliptic.py b/test/test_elliptic.py new file mode 100644 index 00000000..6132c054 --- /dev/null +++ b/test/test_elliptic.py @@ -0,0 +1,63 @@ +from __future__ import annotations + +import pytest + +import meshmode.mesh.generation as mgen +from arraycontext import NumpyArrayContext +from meshmode.mesh import BTAG_ALL, BTAG_NONE +from pytools.convergence import EOCRecorder + +import grudge.geometry as geo +from grudge import op +from grudge.discretization import make_discretization_collection +from grudge.dt_utils import h_max_from_volume +from grudge.models.elliptic import InteriorPenaltyEllipticOperator + + +@pytest.mark.parametrize("bc_kind", ["dirichlet", "neumann"]) +def test_sipg_elliptic_h_convergence_quadratic(bc_kind: str) -> None: + actx = NumpyArrayContext() + eoc = EOCRecorder() + + for n in [2, 4, 8]: + mesh = mgen.generate_regular_rect_mesh( + a=(-1.0, -1.0), b=(1.0, 1.0), nelements_per_axis=(n, n)) + dcoll = make_discretization_collection(actx, mesh, order=2) + + x = actx.thaw(dcoll.nodes()) + u_exact = x[0]**2 + x[0]*x[1] + x[1]**2 + 1.0 + rhs = 0.0 * x[0] - 4.0 + + if bc_kind == "dirichlet": + dir_bc = op.project(dcoll, "vol", BTAG_ALL, u_exact) + + elliptic_op = InteriorPenaltyEllipticOperator( + dcoll, + penalty_factor=10.0, + dirichlet_tag=BTAG_ALL, + dirichlet_bc=dir_bc, + neumann_tag=BTAG_NONE, + ) + elif bc_kind == "neumann": + bdry_dd = BTAG_ALL + bdry_normal = geo.normal(actx, dcoll, bdry_dd) + grad_u = op.local_grad(dcoll, u_exact) + grad_u_bdry = op.project(dcoll, "vol", bdry_dd, grad_u) + neu_bc = grad_u_bdry[0]*bdry_normal[0] + grad_u_bdry[1]*bdry_normal[1] + + elliptic_op = InteriorPenaltyEllipticOperator( + dcoll, + penalty_factor=10.0, + dirichlet_tag=BTAG_NONE, + neumann_tag=BTAG_ALL, + neumann_bc=neu_bc, + ) + else: + raise ValueError(f"invalid boundary condition type: '{bc_kind}'") + + resid = elliptic_op.operator(u_exact) - rhs + err = actx.to_numpy(op.norm(dcoll, resid, 2)) + h_max = actx.to_numpy(h_max_from_volume(dcoll)) + eoc.add_data_point(h_max, err) + + assert eoc.max_error() < 5.0e-9 or eoc.order_estimate() > 1.0