From 9e132f40d8b67e9710c98eb7361cc44a3d059caf Mon Sep 17 00:00:00 2001 From: xywei Date: Mon, 23 Feb 2026 20:11:59 -0600 Subject: [PATCH 1/6] feat: add transplanted 1d quadrature --- doc/quadrature.rst | 119 ++++++ examples/plot-qbx-transplanted-vs-gauss.py | 285 ++++++++++++++ modepy/__init__.py | 4 + modepy/quadrature/__init__.py | 5 + modepy/quadrature/transplanted.py | 414 +++++++++++++++++++++ modepy/test/test_quadrature.py | 112 ++++++ 6 files changed, 939 insertions(+) create mode 100644 examples/plot-qbx-transplanted-vs-gauss.py create mode 100644 modepy/quadrature/transplanted.py diff --git a/doc/quadrature.rst b/doc/quadrature.rst index 0ed51996..b0b4501a 100644 --- a/doc/quadrature.rst +++ b/doc/quadrature.rst @@ -25,6 +25,125 @@ Clenshaw-Curtis and Fejér quadrature in one dimension :members: :show-inheritance: +Transplanted quadrature in one dimension +---------------------------------------- + +The transplanted maps implemented here include the Hale-Trefethen +conformal-map family and the Kosloff-Tal-Ezer map. + +Given a base rule :math:`(s_i, w_i)` on :math:`[-1,1]`, transplanted quadrature +uses a map :math:`x=g(s)` to build + +.. math:: + + x_i = g(s_i), \qquad \tilde w_i = w_i g'(s_i), + +so that + +.. math:: + + \int_{-1}^1 f(x)\,dx = \int_{-1}^1 f(g(s))g'(s)\,ds + \approx \sum_i \tilde w_i f(x_i). + +Map functions +~~~~~~~~~~~~~ + +.. currentmodule:: modepy.quadrature.transplanted + +Identity map +^^^^^^^^^^^^ + +Use ``map_name="identity"`` for the unmodified base rule. + +.. autofunction:: map_identity + +Sausage polynomial maps +^^^^^^^^^^^^^^^^^^^^^^^ + +Use ``map_name="sausage_d{odd}"`` (for example ``"sausage_d5"``, +``"sausage_d9"``, ``"sausage_d17"``) for odd-degree normalized polynomial +truncations of :math:`\arcsin`. + +.. autofunction:: map_sausage + +Kosloff-Tal-Ezer map +^^^^^^^^^^^^^^^^^^^^ + +Use ``map_name="kte"`` (or ``"kosloff_tal_ezer"``). + +* ``kte_rho`` (``>1``) sets the default parameterization + :math:`\alpha = 2/(\rho + \rho^{-1})`. +* ``kte_alpha`` explicitly sets :math:`\alpha` (must satisfy ``0 1``. + +.. note:: + + The strip map requires interior nodes (``abs(s)<1``), so endpoint rules + (for example Gauss-Lobatto or Clenshaw-Curtis) are not valid with + ``map_name="strip"``. + +.. autofunction:: map_strip + +Map dispatcher +^^^^^^^^^^^^^^ + +.. autofunction:: map_trefethen_transplant + +Quadrature classes +~~~~~~~~~~~~~~~~~~ + +.. currentmodule:: modepy + +Transplanted1DQuadrature +^^^^^^^^^^^^^^^^^^^^^^^^ + +.. autoclass:: Transplanted1DQuadrature + :show-inheritance: + +TransplantedLegendreGaussQuadrature +^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ + +.. autoclass:: TransplantedLegendreGaussQuadrature + :show-inheritance: + +Example +~~~~~~~ + +.. code-block:: python + + import modepy as mp + + q_kte = mp.TransplantedLegendreGaussQuadrature( + 20, + map_name="kte", + kte_rho=1.4, + force_dim_axis=True, + ) + + q_sausage = mp.TransplantedLegendreGaussQuadrature( + 20, + map_name="sausage_d9", + force_dim_axis=True, + ) + +References +~~~~~~~~~~ + +* N. Hale and L. N. Trefethen, *New Quadrature Formulas from Conformal + Maps*, ``SIAM J. Numer. Anal.`` 46(2), 930-948 (2008), + doi:10.1137/07068607X. +* D. Kosloff and H. Tal-Ezer, *A Modified Chebyshev Pseudospectral Method + with an O(N^{-1}) Time Step Restriction*, + ``Journal of Computational Physics`` 104(2), 457-469 (1993), + doi:10.1006/jcph.1993.1044. + Quadratures on the simplex -------------------------- diff --git a/examples/plot-qbx-transplanted-vs-gauss.py b/examples/plot-qbx-transplanted-vs-gauss.py new file mode 100644 index 00000000..7bbf1c06 --- /dev/null +++ b/examples/plot-qbx-transplanted-vs-gauss.py @@ -0,0 +1,285 @@ +from functools import lru_cache +import warnings + +import numpy as np +import pyopencl as cl +import pyopencl.tools as cl_tools +from arraycontext import flatten +import matplotlib.pyplot as plt + +import meshmode.mesh.generation as mgen +import modepy as mp + +from meshmode.discretization import Discretization +from meshmode.discretization.poly_element import ( + PolynomialGivenNodesElementGroup, +) +from modepy.quadrature import Transformed1DQuadrature + +from pytential import GeometryCollection, bind, sym +from pytential.array_context import PyOpenCLArrayContext +from pytential.qbx import QBXLayerPotentialSource + +from scipy.optimize import root_scalar +from scipy.special import ellipk + +from sumpy.expansion.local import LineTaylorLocalExpansion +from sumpy.kernel import LaplaceKernel +from sumpy.qbx import LayerPotential + + +NPANELS, MODE = 8, 8 +QBX_ORDER = 20 +ASSOC_TOL = 0.05 +NPTS = list(range(4, 30)) +MAPS = [ + ("gauss", None), + ("kte", "kte"), + ("strip", "strip"), + ("sausage_d5", "sausage_d5"), + ("sausage_d9", "sausage_d9"), +] +OUT = "/tmp/qbx-transplanted-vs-gauss-2d.png" + +STRIP_RHO: float | None = None +STRIP_SAFETY = 0.5 + + +@lru_cache(maxsize=16) +def strip_map_parameter_m(rho: float) -> float: + target = 4.0 * np.log(rho) / np.pi + + def f(m: float) -> float: + return float(ellipk(1.0 - m) / ellipk(m) - target) + + upper = 1.0 - 1.0e-8 + while f(upper) > 0.0 and 1.0 - upper > 1.0e-16: + upper = 1.0 - (1.0 - upper) / 10.0 + result = root_scalar(f, bracket=(1.0e-14, upper), method="brentq") + if not result.converged: + raise RuntimeError("failed to solve strip-map parameter m") + return float(result.root) + + +def strip_half_width(rho: float) -> float: + return float(np.pi / (4.0 * np.arctanh(strip_map_parameter_m(rho) ** 0.25))) + + +def strip_rho_for_half_width( + target_half_width: float, rho_min: float = 1.05, rho_max: float = 5.0 +) -> float: + if target_half_width <= strip_half_width(rho_min): + return float(rho_min) + if target_half_width >= strip_half_width(rho_max): + return float(rho_max) + + def f(rho: float) -> float: + return strip_half_width(rho) - target_half_width + + return float(root_scalar(f, bracket=(rho_min, rho_max), method="brentq").root) + + +def kte_alpha_for_rho(rho: float) -> float: + return float(2.0 / (rho + 1.0 / rho)) + + +def make_quad(npts: int, map_name: str | None, strip_rho: float) -> mp.Quadrature: + if map_name is None: + return mp.LegendreGaussQuadrature(npts - 1, force_dim_axis=True) + return mp.TransplantedLegendreGaussQuadrature( + npts - 1, + map_name=map_name, + strip_rho=strip_rho, + kte_rho=strip_rho, + force_dim_axis=True, + ) + + +def make_group(order: int, quad: mp.Quadrature): + class _G(PolynomialGivenNodesElementGroup): + def __init__(self, meg): + super().__init__(meg, order, quad.nodes) + + def quadrature_rule(self): + return quad + + def discretization_key(self): + return ( + type(self), + self.dim, + self.order, + tuple(quad.nodes.ravel()), + tuple(quad.weights.ravel()), + ) + + return _G + + +def make_mesh_and_t(panel_edges: np.ndarray, npts: int, unit_nodes: np.ndarray): + with warnings.catch_warnings(): + warnings.filterwarnings( + "ignore", + message=( + "Unimplemented: Cannot check element orientation for a mesh with " + "mesh.dim != mesh.ambient_dim" + ), + category=UserWarning, + ) + return mgen.make_curve_mesh( + mgen.circle, + panel_edges, + order=npts - 1, + unit_nodes=unit_nodes, + node_vertex_consistency_tolerance=False, + return_parametrization_points=True, + ) + + +def source_ds_weights(quad: mp.Quadrature, panel_edges: np.ndarray) -> np.ndarray: + dtw = np.concatenate([ + Transformed1DQuadrature(quad, a, b).weights + for a, b in zip(panel_edges[:-1], panel_edges[1:], strict=True) + ]) + return (2.0 * np.pi) * dtw + + +def gauss_centers_radii(actx, panel_edges: np.ndarray, npts: int): + qg = make_quad(npts, None, 2.0) + mesh, _ = make_mesh_and_t(panel_edges, npts, qg.nodes) + discr = Discretization(actx, mesh, make_group(npts - 1, qg)) + qbx = QBXLayerPotentialSource( + discr, + fine_order=1, + qbx_order=1, + fmm_order=False, + target_association_tolerance=ASSOC_TOL, + ) + places = GeometryCollection(qbx) + centers = actx.to_numpy( + flatten(bind(places, sym.expansion_centers(2, +1))(actx), actx) + ).reshape(2, -1) + radii = actx.to_numpy(flatten(bind(places, sym.expansion_radii(2))(actx), actx)) + return centers, radii + + +def auto_strip_rho(actx, panel_edges: np.ndarray) -> float: + _, radii = gauss_centers_radii(actx, panel_edges, max(NPTS)) + eta_min = float(np.min(radii) / (np.pi / NPANELS)) + target_half_width = STRIP_SAFETY * eta_min + rho = strip_rho_for_half_width(target_half_width) + print( + f"auto strip rho: eta_min={eta_min:.6f}, target_half_width={target_half_width:.6f}, rho={rho:.4f}" + ) + return rho + + +def eval_rule( + actx, + lpot, + panel_edges: np.ndarray, + npts: int, + map_name: str | None, + strip_rho: float, + targets: np.ndarray, + centers: np.ndarray, + radii: np.ndarray, +) -> np.ndarray: + quad = make_quad(npts, map_name, strip_rho) + mesh, t_src = make_mesh_and_t(panel_edges, npts, quad.nodes) + sources = mesh.groups[0].nodes.reshape(2, -1) + sigma = np.cos(MODE * 2.0 * np.pi * t_src) + strengths = (actx.from_numpy(sigma * source_ds_weights(quad, panel_edges)),) + (result,) = lpot( + actx, + actx.from_numpy(targets), + actx.from_numpy(sources), + actx.from_numpy(centers), + strengths, + expansion_radii=actx.from_numpy(radii), + ) + return actx.to_numpy(result) + + +def main() -> None: + panel_edges = np.linspace(0.0, 1.0, NPANELS + 1) + queue = cl.CommandQueue(cl.create_some_context(interactive=False)) + allocator = cl_tools.ImmediateAllocator(queue) + actx = PyOpenCLArrayContext(queue, allocator=allocator) + strip_rho = ( + STRIP_RHO if STRIP_RHO is not None else auto_strip_rho(actx, panel_edges) + ) + print( + f"using strip rho={strip_rho:.4f} (half-width={strip_half_width(strip_rho):.6f})" + ) + + lknl = LaplaceKernel(2) + lpot = LayerPotential( + expansion=LineTaylorLocalExpansion(lknl, QBX_ORDER), + target_kernels=(lknl,), + source_kernels=(lknl,), + ) + + orders, totals = [], [] + errors = {name: [] for name, _ in MAPS} + names = [n for n, _ in MAPS] + + print("QBX convergence on meshmode circle (frozen Gauss targets+centers)") + print("order total_nodes " + " ".join(f"{n:>10s}" for n in names)) + print("-" * (33 + 12 * len(names))) + + for i, npts in enumerate(NPTS): + with warnings.catch_warnings(): + warnings.filterwarnings("ignore", category=cl.CompilerWarning) + + qg = make_quad(npts, None, strip_rho) + tgt_mesh, t_tgt = make_mesh_and_t(panel_edges, npts, qg.nodes) + targets = tgt_mesh.groups[0].nodes.reshape(2, -1) + centers, radii = gauss_centers_radii(actx, panel_edges, npts) + ref = np.cos(MODE * 2.0 * np.pi * t_tgt) / (2.0 * MODE) + + orders.append(npts - 1) + totals.append(NPANELS * npts) + for name, map_name in MAPS: + values = eval_rule( + actx, + lpot, + panel_edges, + npts, + map_name, + strip_rho, + targets, + centers, + radii, + ) + errors[name].append(float(np.max(np.abs(values - ref)))) + + vals = " ".join(f"{errors[name][i]:10.3e}" for name in names) + print(f"{orders[i]:5d} {totals[i]:11d} {vals}") + + fig, ax = plt.subplots(figsize=(7.8, 4.8), constrained_layout=True) + style = { + "gauss": ("o-", "Gauss-Legendre"), + "kte": ( + "D-", + f"KTE (rho={strip_rho:.3f}, alpha={kte_alpha_for_rho(strip_rho):.3f})", + ), + "strip": ("s-", f"Strip (rho={strip_rho:.3f})"), + "sausage_d5": ("^-", "Sausage d5"), + "sausage_d9": ("v-", "Sausage d9"), + } + for name in names: + marker, label = style[name] + ax.semilogy(orders, errors[name], marker, label=label) + ax.set_xlabel("per-panel quadrature order") + ax.set_ylabel("max abs error vs circle eigenvalue") + ax.set_title( + f"2D direct QBX on circle, frozen centers, mode={MODE}, qbx_order={QBX_ORDER}" + ) + ax.grid(True, which="both", alpha=0.35) + ax.legend(loc="best") + fig.savefig(OUT, dpi=160) + print(f"Saved plot to: {OUT}") + + +if __name__ == "__main__": + main() diff --git a/modepy/__init__.py b/modepy/__init__.py index 6f22b142..598173d9 100644 --- a/modepy/__init__.py +++ b/modepy/__init__.py @@ -71,6 +71,8 @@ Quadrature, QuadratureRuleUnavailable, TensorProductQuadrature, + Transplanted1DQuadrature, + TransplantedLegendreGaussQuadrature, ZeroDimensionalQuadrature, quadrature_for_space, ) @@ -137,6 +139,8 @@ "TensorProductQuadrature", "TensorProductShape", "TensorProductSpace", + "Transplanted1DQuadrature", + "TransplantedLegendreGaussQuadrature", "VioreanuRokhlinSimplexQuadrature", "WitherdenVincentQuadrature", "XiaoGimbutasSimplexQuadrature", diff --git a/modepy/quadrature/__init__.py b/modepy/quadrature/__init__.py index 996f56fc..be1444b8 100644 --- a/modepy/quadrature/__init__.py +++ b/modepy/quadrature/__init__.py @@ -357,4 +357,9 @@ def _quadrature_for_tp( # }}} +from modepy.quadrature.transplanted import ( + Transplanted1DQuadrature, + TransplantedLegendreGaussQuadrature, +) + # vim: foldmethod=marker diff --git a/modepy/quadrature/transplanted.py b/modepy/quadrature/transplanted.py new file mode 100644 index 00000000..a96b5c79 --- /dev/null +++ b/modepy/quadrature/transplanted.py @@ -0,0 +1,414 @@ +from __future__ import annotations + + +r""" +.. currentmodule:: modepy.quadrature.transplanted + +Transplanted quadrature applies a smooth map :math:`x=g(s)` to an existing +one-dimensional rule on :math:`[-1,1]`. + +Given base nodes/weights :math:`(s_i, w_i)`, the transplanted rule is + +.. math:: + + x_i = g(s_i), \qquad \tilde w_i = w_i g'(s_i), + +so that + +.. math:: + + \int_{-1}^1 f(x)\,dx = \int_{-1}^1 f(g(s)) g'(s)\,ds + \approx \sum_i \tilde w_i f(x_i). + +The map dispatcher :func:`map_trefethen_transplant` recognizes these map names: + +* ``"identity"`` +* ``"sausage_d{odd}"`` (for example ``"sausage_d5"``, ``"sausage_d9"``, + ``"sausage_d17"``) +* ``"kte"`` or ``"kosloff_tal_ezer"`` +* ``"strip"`` + +Parameter notes: + +* ``strip_rho`` controls the strip-map conformal parameter, with ``strip_rho > 1``. +* ``kte_rho`` controls the default KTE parameterization through + :math:`\alpha = 2 / (\rho + \rho^{-1})`, with ``kte_rho > 1``. +* ``kte_alpha`` explicitly sets :math:`\alpha` (must satisfy ``0 < kte_alpha < 1``) + and overrides ``kte_rho``. + +.. note:: + + The strip map requires interior nodes (``abs(s) < 1``). Endpoint-including + base rules (for example Gauss-Lobatto or Clenshaw-Curtis) are therefore not + valid with ``map_name="strip"``. + +.. autofunction:: map_identity +.. autofunction:: map_sausage +.. autofunction:: map_kosloff_tal_ezer +.. autofunction:: map_strip +.. autofunction:: map_trefethen_transplant + +.. currentmodule:: modepy + +.. autoclass:: Transplanted1DQuadrature +.. autoclass:: TransplantedLegendreGaussQuadrature +""" + +from functools import lru_cache +from math import isfinite +from typing import TYPE_CHECKING, Any + +import numpy as np + +from modepy.quadrature import Quadrature + + +if TYPE_CHECKING: + from modepy.typing import ArrayF + + +def map_identity(s: ArrayF) -> tuple[ArrayF, ArrayF]: + """Identity transplant map on :math:`[-1, 1]`. + + Returns ``(s, 1)``. + """ + return s.copy(), np.ones_like(s) + + +def _arcsin_taylor_coefficients(max_odd_degree: int) -> ArrayF: + if max_odd_degree < 1 or max_odd_degree % 2 == 0: + raise ValueError(f"sausage degree must be positive and odd: {max_odd_degree}") + + nterms = (max_odd_degree + 1) // 2 + coeffs = np.empty(nterms, dtype=np.float64) + coeffs[0] = 1.0 + + for k in range(1, nterms): + km1 = k - 1 + coeffs[k] = coeffs[km1] * (2.0 * km1 + 1.0) ** 2 / (2.0 * k * (2.0 * km1 + 3.0)) + + return coeffs + + +def map_sausage(s: ArrayF, degree: int) -> tuple[ArrayF, ArrayF]: + r"""Odd-degree polynomial sausage map from Hale-Trefethen (2008). + + This is the normalized odd Taylor truncation of :math:`\arcsin(s)` + through the monomial of degree *degree*. + + :arg degree: positive odd degree in ``{1, 3, 5, ...}``. + """ + coeffs = _arcsin_taylor_coefficients(degree) + denom = float(np.sum(coeffs)) + + g = np.zeros_like(s) + gp = np.zeros_like(s) + + for k, c_k in enumerate(coeffs): + power = 2 * k + 1 + g = g + c_k * s**power + gp = gp + c_k * power * s ** (power - 1) + + return g / denom, gp / denom + + +def _default_kte_alpha(rho: float) -> float: + if rho <= 1.0: + raise ValueError(f"KTE parameter rho must be > 1: {rho}") + + return float(2.0 / (rho + 1.0 / rho)) + + +def _kte_alpha(*, rho: float, alpha: float | None) -> float: + if alpha is None: + alpha = _default_kte_alpha(rho) + + if not (0.0 < alpha < 1.0) or not isfinite(alpha): + raise ValueError(f"KTE parameter alpha must satisfy 0 < alpha < 1: {alpha}") + + return float(alpha) + + +def map_kosloff_tal_ezer( + s: ArrayF, + *, + rho: float = 1.4, + alpha: float | None = None, +) -> tuple[ArrayF, ArrayF]: + r"""Kosloff-Tal-Ezer map. + + The map is + + .. math:: + + g(s) = \frac{\arcsin(\alpha s)}{\arcsin(\alpha)}, + + where :math:`0 < \alpha < 1`. + + If *alpha* is not provided, it is chosen from *rho* using + + .. math:: + + \alpha = \frac{2}{\rho + \rho^{-1}}, + + matching the parameter choice discussed by Hale-Trefethen for a + :math:`\rho`-ellipse analyticity model. + + Reference: + D. Kosloff and H. Tal-Ezer, "A Modified Chebyshev Pseudospectral + Method with an O(N^{-1}) Time Step Restriction," Journal of + Computational Physics 104(2), 457-469 (1993), + doi:10.1006/jcph.1993.1044. + """ + alpha = _kte_alpha(rho=rho, alpha=alpha) + denom = np.arcsin(alpha) + + g = np.arcsin(alpha * s) / denom + gp = alpha / (denom * np.sqrt(1.0 - alpha**2 * s**2)) + + return g, gp + + +def _sausage_degree_from_map_name(map_name: str) -> int | None: + if not map_name.startswith("sausage_d"): + return None + + degree_text = map_name[len("sausage_d") :] + if not degree_text.isdigit(): + raise ValueError( + f"unsupported sausage map '{map_name}'. Expected format: sausage_d{{odd}}" + ) + + degree = int(degree_text) + if degree < 1 or degree % 2 == 0 or not isfinite(degree): + raise ValueError( + f"unsupported sausage degree in '{map_name}'. " + "Expected a positive odd degree, e.g. sausage_d5" + ) + + return degree + + +def _require_scipy_for_strip_map() -> tuple[Any, Any, Any]: + try: + from scipy.optimize import root_scalar + from scipy.special import ellipj, ellipk + except ImportError as exc: + raise RuntimeError( + "The Trefethen strip map requires scipy. " + "Install modepy with scipy support to use map_name='strip'." + ) from exc + + return root_scalar, ellipj, ellipk + + +@lru_cache(maxsize=16) +def _strip_map_parameter_m(rho: float) -> float: + if rho <= 1.0: + raise ValueError(f"strip map parameter rho must be > 1: {rho}") + + root_scalar, _, ellipk = _require_scipy_for_strip_map() + + target = 4.0 * np.log(rho) / np.pi + + def f(m: float) -> float: + return float(ellipk(1.0 - m) / ellipk(m) - target) + + upper = 1.0 - 1.0e-8 + while f(upper) > 0.0 and 1.0 - upper > 1.0e-16: + upper = 1.0 - (1.0 - upper) / 10.0 + + if f(upper) > 0.0: + raise RuntimeError(f"failed to bracket strip-map parameter for rho={rho}") + + result = root_scalar(f, bracket=(1.0e-14, upper), method="brentq") + if not result.converged: + raise RuntimeError("failed to solve strip-map parameter m") + + return float(result.root) + + +def map_strip(s: ArrayF, *, rho: float = 1.4) -> tuple[ArrayF, ArrayF]: + r"""Strip map from Hale-Trefethen transplanted quadrature. + + :arg rho: strip parameter, must satisfy ``rho > 1``. + + .. important:: + + This map requires interior nodes (``abs(s) < 1``), so it is intended + for base rules such as Legendre-Gauss that do not include endpoints. + """ + if np.any(np.abs(s) >= 1.0): + raise ValueError("strip map expects interior nodes, i.e. abs(s) < 1") + + _, ellipj, ellipk = _require_scipy_for_strip_map() + + m = _strip_map_parameter_m(rho) + m4 = m**0.25 + K = float(ellipk(m)) + + omega = 2.0 * K * np.arcsin(s) / np.pi + sn, cn, dn, _ = ellipj(omega, m) + + g = np.arctanh(m4 * sn) / np.arctanh(m4) + gp = ( + 2.0 + * K + * m4 + * cn + * dn + / (np.pi * np.sqrt(1.0 - s**2) * (1.0 - np.sqrt(m) * sn**2) * np.arctanh(m4)) + ) + + return g, gp + + +def map_trefethen_transplant( + s: ArrayF, + map_name: str, + *, + strip_rho: float = 1.4, + kte_rho: float = 1.4, + kte_alpha: float | None = None, +) -> tuple[ArrayF, ArrayF]: + """Map 1D nodes to a Trefethen transplanted quadrature rule. + + :arg s: nodes on :math:`[-1, 1]`. + :arg map_name: one of ``identity``, ``sausage_d{odd}``, ``kte``, + ``kosloff_tal_ezer``, ``strip``. + :arg strip_rho: strip-map parameter for ``map_name="strip"``. + :arg kte_rho: KTE parameter for ``map_name in {"kte", "kosloff_tal_ezer"}`` + when ``kte_alpha`` is not supplied. + :arg kte_alpha: explicit KTE :math:`\alpha` override. + + :returns: ``(mapped_nodes, jacobian)``. + + The supported maps are: + + * ``identity``: :func:`map_identity` + * ``sausage_d{odd}``: :func:`map_sausage` + * ``kte`` / ``kosloff_tal_ezer``: :func:`map_kosloff_tal_ezer` + * ``strip``: :func:`map_strip` + + Reference: + N. Hale and L. N. Trefethen, "New Quadrature Formulas from + Conformal Maps," SIAM J. Numer. Anal. 46(2), 930-948 (2008), + doi:10.1137/07068607X. + + D. Kosloff and H. Tal-Ezer, "A Modified Chebyshev Pseudospectral + Method with an O(N^{-1}) Time Step Restriction," Journal of + Computational Physics 104(2), 457-469 (1993), + doi:10.1006/jcph.1993.1044. + """ + if map_name == "identity": + return map_identity(s) + + sausage_degree = _sausage_degree_from_map_name(map_name) + if sausage_degree is not None: + return map_sausage(s, sausage_degree) + + if map_name == "strip": + return map_strip(s, rho=strip_rho) + + if map_name in {"kte", "kosloff_tal_ezer"}: + return map_kosloff_tal_ezer(s, rho=kte_rho, alpha=kte_alpha) + + raise ValueError( + "unsupported map_name " + f"'{map_name}'. Expected one of: " + "identity, sausage_d{odd}, kte, kosloff_tal_ezer, strip" + ) + + +class Transplanted1DQuadrature(Quadrature): + r"""Map an existing 1D quadrature rule using a Trefethen transplant map. + + The transformed rule approximates + + .. math:: + + \int_{-1}^1 f(x)\,dx = \int_{-1}^1 f(g(s)) g'(s)\,ds, + + by mapping existing nodes :math:`s_i` and scaling existing weights :math:`w_i` + with :math:`g'(s_i)`. + + Reference: + N. Hale and L. N. Trefethen, "New Quadrature Formulas from + Conformal Maps," SIAM J. Numer. Anal. 46(2), 930-948 (2008), + doi:10.1137/07068607X. + + D. Kosloff and H. Tal-Ezer, "A Modified Chebyshev Pseudospectral + Method with an O(N^{-1}) Time Step Restriction," Journal of + Computational Physics 104(2), 457-469 (1993), + doi:10.1006/jcph.1993.1044. + """ + + def __init__( + self, + quadrature: Quadrature, + map_name: str = "sausage_d9", + *, + strip_rho: float = 1.4, + kte_rho: float = 1.4, + kte_alpha: float | None = None, + ) -> None: + base_nodes = quadrature.nodes + if base_nodes.ndim == 1: + nodes_1d = base_nodes + force_dim_axis = False + elif base_nodes.ndim == 2 and base_nodes.shape[0] == 1: + nodes_1d = base_nodes[0] + force_dim_axis = True + else: + raise ValueError( + "Transplanted1DQuadrature requires a one-dimensional base quadrature" + ) + + mapped_nodes, jacobian = map_trefethen_transplant( + nodes_1d, + map_name=map_name, + strip_rho=strip_rho, + kte_rho=kte_rho, + kte_alpha=kte_alpha, + ) + mapped_weights = quadrature.weights * jacobian + + if force_dim_axis: + mapped_nodes = mapped_nodes.reshape(1, -1) + + super().__init__(mapped_nodes, mapped_weights) + + self.base_quadrature = quadrature + self.map_name = map_name + self.strip_rho = strip_rho + self.kte_rho = kte_rho + self.kte_alpha = kte_alpha + + +class TransplantedLegendreGaussQuadrature(Transplanted1DQuadrature): + r"""Legendre-Gauss quadrature transplanted by a Trefethen map.""" + + def __init__( + self, + N: int, # noqa: N803 + map_name: str = "sausage_d9", + *, + strip_rho: float = 1.4, + kte_rho: float = 1.4, + kte_alpha: float | None = None, + backend: str | None = None, + force_dim_axis: bool = False, + ) -> None: + from modepy.quadrature.jacobi_gauss import LegendreGaussQuadrature + + super().__init__( + LegendreGaussQuadrature( + N, + backend=backend, + force_dim_axis=force_dim_axis, + ), + map_name=map_name, + strip_rho=strip_rho, + kte_rho=kte_rho, + kte_alpha=kte_alpha, + ) diff --git a/modepy/test/test_quadrature.py b/modepy/test/test_quadrature.py index 955f6d21..c9de124f 100644 --- a/modepy/test/test_quadrature.py +++ b/modepy/test/test_quadrature.py @@ -106,6 +106,118 @@ def f(x): assert err < 2.0e-15, (s, deg, err, i_f, i_f_true) +@pytest.mark.parametrize("map_name", [ + "identity", + "sausage_d5", + "sausage_d9", + "sausage_d17", + "kte", + "kosloff_tal_ezer", +]) +def test_transplanted_legendre_gauss_quadrature(map_name: str) -> None: + base = mp.LegendreGaussQuadrature(15, force_dim_axis=True) + transplanted = mp.TransplantedLegendreGaussQuadrature( + 15, + map_name=map_name, + force_dim_axis=True, + ) + + assert transplanted.nodes.shape == base.nodes.shape + assert transplanted.weights.shape == base.weights.shape + + from modepy.quadrature.transplanted import map_trefethen_transplant + + mapped_nodes, mapped_jacobian = map_trefethen_transplant( + base.nodes[0], map_name=map_name + ) + + assert la.norm(transplanted.nodes[0] - mapped_nodes, np.inf) < 1.0e-14 + assert ( + la.norm(transplanted.weights - base.weights * mapped_jacobian, np.inf) < 1.0e-14 + ) + + # Sausage maps are polynomial maps, so integrating constants + # should remain exact. + if map_name.startswith("sausage_d"): + err = abs(np.sum(transplanted.weights) - 2.0) + assert err < 1.0e-14 + + +def test_transplanted_strip_map_quadrature() -> None: + pytest.importorskip("scipy") + + transplanted = mp.TransplantedLegendreGaussQuadrature( + 32, + map_name="strip", + strip_rho=1.4, + force_dim_axis=True, + ) + + assert np.all(np.diff(transplanted.nodes[0]) > 0) + assert np.all(transplanted.weights > 0) + assert abs(np.sum(transplanted.weights) - 2.0) < 1.0e-9 + + +def test_transplanted_strip_map_rejects_endpoint_rules() -> None: + pytest.importorskip("scipy") + + with pytest.raises(ValueError, match="interior nodes"): + mp.Transplanted1DQuadrature( + mp.ClenshawCurtisQuadrature(5, force_dim_axis=True), map_name="strip" + ) + + +def test_transplanted_sausage_map_rejects_even_degrees() -> None: + with pytest.raises(ValueError, match="positive odd degree"): + mp.TransplantedLegendreGaussQuadrature( + 8, map_name="sausage_d4", force_dim_axis=True + ) + + +def test_transplanted_sausage_map_name_matches_direct_map() -> None: + from modepy.quadrature.transplanted import map_sausage, map_trefethen_transplant + + s = np.linspace(-0.95, 0.95, 33) + + for degree in [5, 9, 17]: + g, gp = map_sausage(s, degree=degree) + g_ref, gp_ref = map_trefethen_transplant(s, map_name=f"sausage_d{degree}") + assert la.norm(g - g_ref, np.inf) < 1.0e-15 + assert la.norm(gp - gp_ref, np.inf) < 1.0e-15 + + +def test_transplanted_kte_map_name_matches_direct_map() -> None: + from modepy.quadrature.transplanted import ( + map_kosloff_tal_ezer, + map_trefethen_transplant, + ) + + s = np.linspace(-1.0, 1.0, 33) + rho = 1.4 + alpha = 2.0 / (rho + 1.0 / rho) + + g, gp = map_kosloff_tal_ezer(s, rho=rho) + g_ref, gp_ref = map_kosloff_tal_ezer(s, alpha=alpha) + assert la.norm(g - g_ref, np.inf) < 1.0e-15 + assert la.norm(gp - gp_ref, np.inf) < 1.0e-15 + + g_name, gp_name = map_trefethen_transplant(s, map_name="kte", kte_rho=rho) + assert la.norm(g - g_name, np.inf) < 1.0e-15 + assert la.norm(gp - gp_name, np.inf) < 1.0e-15 + + +def test_transplanted_kte_map_rejects_invalid_parameters() -> None: + from modepy.quadrature.transplanted import map_trefethen_transplant + + s = np.array([0.0]) + + with pytest.raises(ValueError, match="rho must be > 1"): + map_trefethen_transplant(s, map_name="kte", kte_rho=1.0) + + with pytest.raises(ValueError, match="0 < alpha < 1"): + map_trefethen_transplant(s, map_name="kte", kte_alpha=1.0) + + @pytest.mark.parametrize("kind", [1, 2]) def test_fejer_quadrature(kind: int) -> None: from modepy.quadrature.clenshaw_curtis import FejerQuadrature From f9884a8a33155360e2ccb807575b3d83587d13d4 Mon Sep 17 00:00:00 2001 From: xywei Date: Tue, 24 Feb 2026 09:58:23 -0600 Subject: [PATCH 2/6] fix: unblock CI --- doc/quadrature.rst | 2 +- examples/plot-qbx-transplanted-vs-gauss.py | 32 ++++++++++++---------- modepy/quadrature/__init__.py | 5 ++-- modepy/quadrature/transplanted.py | 12 ++++---- 4 files changed, 29 insertions(+), 22 deletions(-) diff --git a/doc/quadrature.rst b/doc/quadrature.rst index b0b4501a..98b9bf84 100644 --- a/doc/quadrature.rst +++ b/doc/quadrature.rst @@ -137,7 +137,7 @@ References ~~~~~~~~~~ * N. Hale and L. N. Trefethen, *New Quadrature Formulas from Conformal - Maps*, ``SIAM J. Numer. Anal.`` 46(2), 930-948 (2008), + Maps*, ``SIAM Journal on Numerical Analysis`` 46(2), 930-948 (2008), doi:10.1137/07068607X. * D. Kosloff and H. Tal-Ezer, *A Modified Chebyshev Pseudospectral Method with an O(N^{-1}) Time Step Restriction*, diff --git a/examples/plot-qbx-transplanted-vs-gauss.py b/examples/plot-qbx-transplanted-vs-gauss.py index 7bbf1c06..77f41380 100644 --- a/examples/plot-qbx-transplanted-vs-gauss.py +++ b/examples/plot-qbx-transplanted-vs-gauss.py @@ -1,32 +1,32 @@ -from functools import lru_cache +from __future__ import annotations + +# pyright: reportMissingImports=false import warnings +from functools import lru_cache +from itertools import pairwise +import matplotlib.pyplot as plt +import meshmode.mesh.generation as mgen import numpy as np import pyopencl as cl import pyopencl.tools as cl_tools from arraycontext import flatten -import matplotlib.pyplot as plt - -import meshmode.mesh.generation as mgen -import modepy as mp - from meshmode.discretization import Discretization from meshmode.discretization.poly_element import ( PolynomialGivenNodesElementGroup, ) -from modepy.quadrature import Transformed1DQuadrature - from pytential import GeometryCollection, bind, sym from pytential.array_context import PyOpenCLArrayContext from pytential.qbx import QBXLayerPotentialSource - from scipy.optimize import root_scalar from scipy.special import ellipk - from sumpy.expansion.local import LineTaylorLocalExpansion from sumpy.kernel import LaplaceKernel from sumpy.qbx import LayerPotential +import modepy as mp +from modepy.quadrature import Transformed1DQuadrature + NPANELS, MODE = 8, 8 QBX_ORDER = 20 @@ -137,8 +137,7 @@ def make_mesh_and_t(panel_edges: np.ndarray, npts: int, unit_nodes: np.ndarray): def source_ds_weights(quad: mp.Quadrature, panel_edges: np.ndarray) -> np.ndarray: dtw = np.concatenate([ - Transformed1DQuadrature(quad, a, b).weights - for a, b in zip(panel_edges[:-1], panel_edges[1:], strict=True) + Transformed1DQuadrature(quad, a, b).weights for a, b in pairwise(panel_edges) ]) return (2.0 * np.pi) * dtw @@ -168,7 +167,10 @@ def auto_strip_rho(actx, panel_edges: np.ndarray) -> float: target_half_width = STRIP_SAFETY * eta_min rho = strip_rho_for_half_width(target_half_width) print( - f"auto strip rho: eta_min={eta_min:.6f}, target_half_width={target_half_width:.6f}, rho={rho:.4f}" + "auto strip rho: " + f"eta_min={eta_min:.6f}, " + f"target_half_width={target_half_width:.6f}, " + f"rho={rho:.4f}" ) return rho @@ -209,7 +211,9 @@ def main() -> None: STRIP_RHO if STRIP_RHO is not None else auto_strip_rho(actx, panel_edges) ) print( - f"using strip rho={strip_rho:.4f} (half-width={strip_half_width(strip_rho):.6f})" + "using strip " + f"rho={strip_rho:.4f} " + f"(half-width={strip_half_width(strip_rho):.6f})" ) lknl = LaplaceKernel(2) diff --git a/modepy/quadrature/__init__.py b/modepy/quadrature/__init__.py index be1444b8..2a12503e 100644 --- a/modepy/quadrature/__init__.py +++ b/modepy/quadrature/__init__.py @@ -357,9 +357,10 @@ def _quadrature_for_tp( # }}} + from modepy.quadrature.transplanted import ( - Transplanted1DQuadrature, - TransplantedLegendreGaussQuadrature, + Transplanted1DQuadrature as Transplanted1DQuadrature, + TransplantedLegendreGaussQuadrature as TransplantedLegendreGaussQuadrature, ) # vim: foldmethod=marker diff --git a/modepy/quadrature/transplanted.py b/modepy/quadrature/transplanted.py index a96b5c79..c7f306a8 100644 --- a/modepy/quadrature/transplanted.py +++ b/modepy/quadrature/transplanted.py @@ -245,15 +245,15 @@ def map_strip(s: ArrayF, *, rho: float = 1.4) -> tuple[ArrayF, ArrayF]: m = _strip_map_parameter_m(rho) m4 = m**0.25 - K = float(ellipk(m)) + k = float(ellipk(m)) - omega = 2.0 * K * np.arcsin(s) / np.pi + omega = 2.0 * k * np.arcsin(s) / np.pi sn, cn, dn, _ = ellipj(omega, m) g = np.arctanh(m4 * sn) / np.arctanh(m4) gp = ( 2.0 - * K + * k * m4 * cn * dn @@ -292,7 +292,8 @@ def map_trefethen_transplant( Reference: N. Hale and L. N. Trefethen, "New Quadrature Formulas from - Conformal Maps," SIAM J. Numer. Anal. 46(2), 930-948 (2008), + Conformal Maps," SIAM Journal on Numerical Analysis 46(2), + 930-948 (2008), doi:10.1137/07068607X. D. Kosloff and H. Tal-Ezer, "A Modified Chebyshev Pseudospectral @@ -334,7 +335,8 @@ class Transplanted1DQuadrature(Quadrature): Reference: N. Hale and L. N. Trefethen, "New Quadrature Formulas from - Conformal Maps," SIAM J. Numer. Anal. 46(2), 930-948 (2008), + Conformal Maps," SIAM Journal on Numerical Analysis 46(2), + 930-948 (2008), doi:10.1137/07068607X. D. Kosloff and H. Tal-Ezer, "A Modified Chebyshev Pseudospectral From 5f2d7a96857d345c9b17bdf0889d73ded0669b74 Mon Sep 17 00:00:00 2001 From: xywei Date: Tue, 24 Feb 2026 10:29:26 -0600 Subject: [PATCH 3/6] fix: less strict checking for the qbx example script --- examples/plot-qbx-transplanted-vs-gauss.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/examples/plot-qbx-transplanted-vs-gauss.py b/examples/plot-qbx-transplanted-vs-gauss.py index 77f41380..056d04e8 100644 --- a/examples/plot-qbx-transplanted-vs-gauss.py +++ b/examples/plot-qbx-transplanted-vs-gauss.py @@ -1,6 +1,7 @@ from __future__ import annotations -# pyright: reportMissingImports=false +# Optional QBX dependencies (meshmode/pytential/sumpy) are not installed in CI. +# pyright: basic, reportMissingImports=false import warnings from functools import lru_cache from itertools import pairwise From 3efe064499db878a7af1beb53e4f3ceb3ffa7001 Mon Sep 17 00:00:00 2001 From: xywei Date: Tue, 24 Feb 2026 11:00:11 -0600 Subject: [PATCH 4/6] fix: add typing to transplanted.py --- modepy/quadrature/transplanted.py | 140 ++++++++++++++++++++---------- 1 file changed, 96 insertions(+), 44 deletions(-) diff --git a/modepy/quadrature/transplanted.py b/modepy/quadrature/transplanted.py index c7f306a8..c895facb 100644 --- a/modepy/quadrature/transplanted.py +++ b/modepy/quadrature/transplanted.py @@ -55,8 +55,9 @@ """ from functools import lru_cache -from math import isfinite -from typing import TYPE_CHECKING, Any +from importlib import import_module +from math import asin, isfinite, log, sqrt +from typing import TYPE_CHECKING, Protocol, cast import numpy as np @@ -64,30 +65,75 @@ if TYPE_CHECKING: + from collections.abc import Callable + from modepy.typing import ArrayF +class _RootScalarResult(Protocol): + converged: bool + root: float + + +class _RootScalarFn(Protocol): + def __call__( + self, + f: Callable[[float], float], + *, + bracket: tuple[float, float], + method: str, + ) -> _RootScalarResult: ... + + +class _EllipkFn(Protocol): + def __call__(self, x: float) -> float: ... + + +class _EllipjFn(Protocol): + def __call__( + self, u: ArrayF, m: float + ) -> tuple[ArrayF, ArrayF, ArrayF, ArrayF]: ... + + +def _scipy_attr(module_name: str, attr_name: str) -> object: + try: + module = import_module(module_name) + except ImportError as exc: + raise RuntimeError( + "The Trefethen strip map requires scipy. " + "Install modepy with scipy support to use map_name='strip'." + ) from exc + + try: + return cast("object", getattr(module, attr_name)) + except AttributeError as exc: + raise RuntimeError( + f"scipy module '{module_name}' is missing required attribute '{attr_name}'" + ) from exc + + def map_identity(s: ArrayF) -> tuple[ArrayF, ArrayF]: """Identity transplant map on :math:`[-1, 1]`. Returns ``(s, 1)``. """ - return s.copy(), np.ones_like(s) + return np.array(s, dtype=np.float64, copy=True), np.ones_like(s, dtype=np.float64) -def _arcsin_taylor_coefficients(max_odd_degree: int) -> ArrayF: +def _arcsin_taylor_coefficients(max_odd_degree: int) -> tuple[float, ...]: if max_odd_degree < 1 or max_odd_degree % 2 == 0: raise ValueError(f"sausage degree must be positive and odd: {max_odd_degree}") nterms = (max_odd_degree + 1) // 2 - coeffs = np.empty(nterms, dtype=np.float64) - coeffs[0] = 1.0 + coeffs = [1.0] for k in range(1, nterms): km1 = k - 1 - coeffs[k] = coeffs[km1] * (2.0 * km1 + 1.0) ** 2 / (2.0 * k * (2.0 * km1 + 3.0)) + coeffs.append( + coeffs[km1] * (2.0 * km1 + 1.0) ** 2 / (2.0 * k * (2.0 * km1 + 3.0)) + ) - return coeffs + return tuple(coeffs) def map_sausage(s: ArrayF, degree: int) -> tuple[ArrayF, ArrayF]: @@ -99,10 +145,10 @@ def map_sausage(s: ArrayF, degree: int) -> tuple[ArrayF, ArrayF]: :arg degree: positive odd degree in ``{1, 3, 5, ...}``. """ coeffs = _arcsin_taylor_coefficients(degree) - denom = float(np.sum(coeffs)) + denom = float(sum(coeffs)) - g = np.zeros_like(s) - gp = np.zeros_like(s) + g = np.zeros_like(s, dtype=np.float64) + gp = np.zeros_like(s, dtype=np.float64) for k, c_k in enumerate(coeffs): power = 2 * k + 1 @@ -161,10 +207,10 @@ def map_kosloff_tal_ezer( doi:10.1006/jcph.1993.1044. """ alpha = _kte_alpha(rho=rho, alpha=alpha) - denom = np.arcsin(alpha) + denom = asin(alpha) - g = np.arcsin(alpha * s) / denom - gp = alpha / (denom * np.sqrt(1.0 - alpha**2 * s**2)) + g = np.asarray(np.arcsin(alpha * s) / denom, dtype=np.float64) + gp = np.asarray(alpha / (denom * np.sqrt(1.0 - alpha**2 * s**2)), dtype=np.float64) return g, gp @@ -189,27 +235,15 @@ def _sausage_degree_from_map_name(map_name: str) -> int | None: return degree -def _require_scipy_for_strip_map() -> tuple[Any, Any, Any]: - try: - from scipy.optimize import root_scalar - from scipy.special import ellipj, ellipk - except ImportError as exc: - raise RuntimeError( - "The Trefethen strip map requires scipy. " - "Install modepy with scipy support to use map_name='strip'." - ) from exc - - return root_scalar, ellipj, ellipk - - @lru_cache(maxsize=16) def _strip_map_parameter_m(rho: float) -> float: if rho <= 1.0: raise ValueError(f"strip map parameter rho must be > 1: {rho}") - root_scalar, _, ellipk = _require_scipy_for_strip_map() + root_scalar = cast("_RootScalarFn", _scipy_attr("scipy.optimize", "root_scalar")) + ellipk = cast("_EllipkFn", _scipy_attr("scipy.special", "ellipk")) - target = 4.0 * np.log(rho) / np.pi + target = 4.0 * log(rho) / np.pi def f(m: float) -> float: return float(ellipk(1.0 - m) / ellipk(m) - target) @@ -241,23 +275,35 @@ def map_strip(s: ArrayF, *, rho: float = 1.4) -> tuple[ArrayF, ArrayF]: if np.any(np.abs(s) >= 1.0): raise ValueError("strip map expects interior nodes, i.e. abs(s) < 1") - _, ellipj, ellipk = _require_scipy_for_strip_map() + ellipj = cast("_EllipjFn", _scipy_attr("scipy.special", "ellipj")) + ellipk = cast("_EllipkFn", _scipy_attr("scipy.special", "ellipk")) m = _strip_map_parameter_m(rho) - m4 = m**0.25 + m4 = sqrt(sqrt(m)) k = float(ellipk(m)) omega = 2.0 * k * np.arcsin(s) / np.pi - sn, cn, dn, _ = ellipj(omega, m) - - g = np.arctanh(m4 * sn) / np.arctanh(m4) - gp = ( - 2.0 - * k - * m4 - * cn - * dn - / (np.pi * np.sqrt(1.0 - s**2) * (1.0 - np.sqrt(m) * sn**2) * np.arctanh(m4)) + sn_jacobi, cn_jacobi, dn_jacobi, _ = ellipj(omega, m) + sn = np.asarray(sn_jacobi, dtype=np.float64) + cn = np.asarray(cn_jacobi, dtype=np.float64) + dn = np.asarray(dn_jacobi, dtype=np.float64) + + g = np.asarray(np.arctanh(m4 * sn) / np.arctanh(m4), dtype=np.float64) + gp = np.asarray( + ( + 2.0 + * k + * m4 + * cn + * dn + / ( + np.pi + * np.sqrt(1.0 - s**2) + * (1.0 - np.sqrt(m) * sn**2) + * np.arctanh(m4) + ) + ), + dtype=np.float64, ) return g, gp @@ -345,6 +391,12 @@ class Transplanted1DQuadrature(Quadrature): doi:10.1006/jcph.1993.1044. """ + base_quadrature: Quadrature + map_name: str + strip_rho: float + kte_rho: float + kte_alpha: float | None + def __init__( self, quadrature: Quadrature, @@ -356,10 +408,10 @@ def __init__( ) -> None: base_nodes = quadrature.nodes if base_nodes.ndim == 1: - nodes_1d = base_nodes + nodes_1d = np.asarray(base_nodes, dtype=np.float64) force_dim_axis = False elif base_nodes.ndim == 2 and base_nodes.shape[0] == 1: - nodes_1d = base_nodes[0] + nodes_1d = np.asarray(base_nodes[0], dtype=np.float64) force_dim_axis = True else: raise ValueError( @@ -376,7 +428,7 @@ def __init__( mapped_weights = quadrature.weights * jacobian if force_dim_axis: - mapped_nodes = mapped_nodes.reshape(1, -1) + mapped_nodes = np.reshape(mapped_nodes, (1, mapped_nodes.shape[0])) super().__init__(mapped_nodes, mapped_weights) From 0b90d2fd0a12cb703fadc491ae5490e7b426cc25 Mon Sep 17 00:00:00 2001 From: xywei Date: Tue, 24 Feb 2026 11:11:50 -0600 Subject: [PATCH 5/6] fix: add typing to test_quadrature.py --- modepy/test/test_quadrature.py | 27 ++++++++++++++++++++------- 1 file changed, 20 insertions(+), 7 deletions(-) diff --git a/modepy/test/test_quadrature.py b/modepy/test/test_quadrature.py index c9de124f..6ddfaf4a 100644 --- a/modepy/test/test_quadrature.py +++ b/modepy/test/test_quadrature.py @@ -25,6 +25,7 @@ import logging +from typing import TYPE_CHECKING, cast import numpy as np import numpy.linalg as la @@ -33,6 +34,10 @@ import modepy as mp +if TYPE_CHECKING: + from modepy.typing import ArrayF + + logger = logging.getLogger(__name__) @@ -122,24 +127,29 @@ def test_transplanted_legendre_gauss_quadrature(map_name: str) -> None: force_dim_axis=True, ) + base_nodes = cast("ArrayF", np.asarray(base.nodes[0], dtype=np.float64)) + trans_nodes = cast("ArrayF", np.asarray(transplanted.nodes[0], dtype=np.float64)) + base_weights = cast("ArrayF", np.asarray(base.weights, dtype=np.float64)) + trans_weights = cast("ArrayF", np.asarray(transplanted.weights, dtype=np.float64)) + assert transplanted.nodes.shape == base.nodes.shape assert transplanted.weights.shape == base.weights.shape from modepy.quadrature.transplanted import map_trefethen_transplant mapped_nodes, mapped_jacobian = map_trefethen_transplant( - base.nodes[0], map_name=map_name + base_nodes, map_name=map_name ) - assert la.norm(transplanted.nodes[0] - mapped_nodes, np.inf) < 1.0e-14 + assert la.norm(trans_nodes - mapped_nodes, np.inf) < 1.0e-14 assert ( - la.norm(transplanted.weights - base.weights * mapped_jacobian, np.inf) < 1.0e-14 + la.norm(trans_weights - base_weights * mapped_jacobian, np.inf) < 1.0e-14 ) # Sausage maps are polynomial maps, so integrating constants # should remain exact. if map_name.startswith("sausage_d"): - err = abs(np.sum(transplanted.weights) - 2.0) + err = abs(float(np.sum(trans_weights)) - 2.0) assert err < 1.0e-14 @@ -153,9 +163,12 @@ def test_transplanted_strip_map_quadrature() -> None: force_dim_axis=True, ) - assert np.all(np.diff(transplanted.nodes[0]) > 0) - assert np.all(transplanted.weights > 0) - assert abs(np.sum(transplanted.weights) - 2.0) < 1.0e-9 + strip_nodes = cast("ArrayF", np.asarray(transplanted.nodes[0], dtype=np.float64)) + strip_weights = cast("ArrayF", np.asarray(transplanted.weights, dtype=np.float64)) + + assert np.all(np.diff(strip_nodes) > 0) + assert np.all(strip_weights > 0) + assert abs(float(np.sum(strip_weights)) - 2.0) < 1.0e-9 def test_transplanted_strip_map_rejects_endpoint_rules() -> None: From dc4f215f3daa70a11cf5e32044b2b20fd3dfa69d Mon Sep 17 00:00:00 2001 From: xywei Date: Tue, 24 Feb 2026 11:38:29 -0600 Subject: [PATCH 6/6] fix: address copilot reviews --- examples/plot-qbx-transplanted-vs-gauss.py | 4 +++- modepy/quadrature/transplanted.py | 11 ++++++++++- 2 files changed, 13 insertions(+), 2 deletions(-) diff --git a/examples/plot-qbx-transplanted-vs-gauss.py b/examples/plot-qbx-transplanted-vs-gauss.py index 056d04e8..934be1a1 100644 --- a/examples/plot-qbx-transplanted-vs-gauss.py +++ b/examples/plot-qbx-transplanted-vs-gauss.py @@ -2,9 +2,11 @@ # Optional QBX dependencies (meshmode/pytential/sumpy) are not installed in CI. # pyright: basic, reportMissingImports=false +import tempfile import warnings from functools import lru_cache from itertools import pairwise +from pathlib import Path import matplotlib.pyplot as plt import meshmode.mesh.generation as mgen @@ -40,7 +42,7 @@ ("sausage_d5", "sausage_d5"), ("sausage_d9", "sausage_d9"), ] -OUT = "/tmp/qbx-transplanted-vs-gauss-2d.png" +OUT = Path(tempfile.gettempdir()) / "qbx-transplanted-vs-gauss-2d.png" STRIP_RHO: float | None = None STRIP_SAFETY = 0.5 diff --git a/modepy/quadrature/transplanted.py b/modepy/quadrature/transplanted.py index c895facb..c4ae4170 100644 --- a/modepy/quadrature/transplanted.py +++ b/modepy/quadrature/transplanted.py @@ -235,6 +235,14 @@ def _sausage_degree_from_map_name(map_name: str) -> int | None: return degree +def _map_preserves_exact_to(map_name: str) -> bool: + if map_name == "identity": + return True + + sausage_degree = _sausage_degree_from_map_name(map_name) + return sausage_degree == 1 + + @lru_cache(maxsize=16) def _strip_map_parameter_m(rho: float) -> float: if rho <= 1.0: @@ -430,7 +438,8 @@ def __init__( if force_dim_axis: mapped_nodes = np.reshape(mapped_nodes, (1, mapped_nodes.shape[0])) - super().__init__(mapped_nodes, mapped_weights) + exact_to = quadrature._exact_to if _map_preserves_exact_to(map_name) else None + super().__init__(mapped_nodes, mapped_weights, exact_to=exact_to) self.base_quadrature = quadrature self.map_name = map_name