From 9e132f40d8b67e9710c98eb7361cc44a3d059caf Mon Sep 17 00:00:00 2001 From: xywei Date: Mon, 23 Feb 2026 20:11:59 -0600 Subject: [PATCH 01/11] 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 02/11] 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 03/11] 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 04/11] 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 05/11] 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 06/11] 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 From 567e9a98f821a15d99184255db75aaed711fef07 Mon Sep 17 00:00:00 2001 From: xywei Date: Wed, 25 Feb 2026 09:03:56 -0600 Subject: [PATCH 07/11] fix: misc fixes in docs and interface --- doc/quadrature.rst | 45 ++-- examples/plot-qbx-transplanted-vs-gauss.py | 30 ++- modepy/__init__.py | 10 +- modepy/quadrature/__init__.py | 6 - modepy/quadrature/transplanted.py | 248 ++++++++------------- 5 files changed, 144 insertions(+), 195 deletions(-) diff --git a/doc/quadrature.rst b/doc/quadrature.rst index 98b9bf84..5691cfee 100644 --- a/doc/quadrature.rst +++ b/doc/quadrature.rst @@ -25,24 +25,26 @@ Clenshaw-Curtis and Fejér quadrature in one dimension :members: :show-inheritance: +.. _quadrature-transplanted-1d: + 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 +Given a base rule :math:`(s_i, w_i^{(s)})` 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), + x_i = g(s_i), \qquad \tilde w_i = w_i^{(s)} g'(s_i), so that .. math:: - \int_{-1}^1 f(x)\,dx = \int_{-1}^1 f(g(s))g'(s)\,ds + \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 @@ -60,8 +62,8 @@ Use ``map_name="identity"`` for the unmodified base rule. Sausage polynomial maps ^^^^^^^^^^^^^^^^^^^^^^^ -Use ``map_name="sausage_d{odd}"`` (for example ``"sausage_d5"``, -``"sausage_d9"``, ``"sausage_d17"``) for odd-degree normalized polynomial +Use ``map_name="sausage"`` with odd ``sausage_degree`` (for example +``sausage_degree=5``, ``9``, ``17``) for odd-degree normalized polynomial truncations of :math:`\arcsin`. .. autofunction:: map_sausage @@ -96,22 +98,14 @@ Map dispatcher .. autofunction:: map_trefethen_transplant -Quadrature classes -~~~~~~~~~~~~~~~~~~ +Quadrature wrappers +~~~~~~~~~~~~~~~~~~~ .. currentmodule:: modepy -Transplanted1DQuadrature -^^^^^^^^^^^^^^^^^^^^^^^^ - -.. autoclass:: Transplanted1DQuadrature - :show-inheritance: +.. autofunction:: transplanted_1d_quadrature -TransplantedLegendreGaussQuadrature -^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ - -.. autoclass:: TransplantedLegendreGaussQuadrature - :show-inheritance: +.. autofunction:: transplanted_legendre_gauss_quadrature Example ~~~~~~~ @@ -120,16 +114,17 @@ Example import modepy as mp - q_kte = mp.TransplantedLegendreGaussQuadrature( + q_kte = mp.transplanted_legendre_gauss_quadrature( 20, map_name="kte", kte_rho=1.4, force_dim_axis=True, ) - q_sausage = mp.TransplantedLegendreGaussQuadrature( + q_sausage = mp.transplanted_legendre_gauss_quadrature( 20, - map_name="sausage_d9", + map_name="sausage", + sausage_degree=9, force_dim_axis=True, ) @@ -137,12 +132,12 @@ References ~~~~~~~~~~ * N. Hale and L. N. Trefethen, *New Quadrature Formulas from Conformal - Maps*, ``SIAM Journal on Numerical Analysis`` 46(2), 930-948 (2008), - doi:10.1137/07068607X. + 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*, - ``Journal of Computational Physics`` 104(2), 457-469 (1993), - doi:10.1006/jcph.1993.1044. + with an :math:`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 index 934be1a1..b7c45764 100644 --- a/examples/plot-qbx-transplanted-vs-gauss.py +++ b/examples/plot-qbx-transplanted-vs-gauss.py @@ -36,11 +36,11 @@ ASSOC_TOL = 0.05 NPTS = list(range(4, 30)) MAPS = [ - ("gauss", None), - ("kte", "kte"), - ("strip", "strip"), - ("sausage_d5", "sausage_d5"), - ("sausage_d9", "sausage_d9"), + ("gauss", None, None), + ("kte", "kte", None), + ("strip", "strip", None), + ("sausage_d5", "sausage", 5), + ("sausage_d9", "sausage", 9), ] OUT = Path(tempfile.gettempdir()) / "qbx-transplanted-vs-gauss-2d.png" @@ -86,12 +86,18 @@ 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: +def make_quad( + npts: int, + map_name: str | None, + strip_rho: float, + sausage_degree: int | None = None, +) -> mp.Quadrature: if map_name is None: return mp.LegendreGaussQuadrature(npts - 1, force_dim_axis=True) - return mp.TransplantedLegendreGaussQuadrature( + return mp.transplanted_legendre_gauss_quadrature( npts - 1, map_name=map_name, + sausage_degree=9 if sausage_degree is None else sausage_degree, strip_rho=strip_rho, kte_rho=strip_rho, force_dim_axis=True, @@ -184,12 +190,13 @@ def eval_rule( panel_edges: np.ndarray, npts: int, map_name: str | None, + sausage_degree: int | None, strip_rho: float, targets: np.ndarray, centers: np.ndarray, radii: np.ndarray, ) -> np.ndarray: - quad = make_quad(npts, map_name, strip_rho) + quad = make_quad(npts, map_name, strip_rho, sausage_degree) 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) @@ -227,8 +234,8 @@ def main() -> None: ) orders, totals = [], [] - errors = {name: [] for name, _ in MAPS} - names = [n for n, _ in MAPS] + 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)) @@ -246,13 +253,14 @@ def main() -> None: orders.append(npts - 1) totals.append(NPANELS * npts) - for name, map_name in MAPS: + for name, map_name, sausage_degree in MAPS: values = eval_rule( actx, lpot, panel_edges, npts, map_name, + sausage_degree, strip_rho, targets, centers, diff --git a/modepy/__init__.py b/modepy/__init__.py index 598173d9..40194029 100644 --- a/modepy/__init__.py +++ b/modepy/__init__.py @@ -71,8 +71,6 @@ Quadrature, QuadratureRuleUnavailable, TensorProductQuadrature, - Transplanted1DQuadrature, - TransplantedLegendreGaussQuadrature, ZeroDimensionalQuadrature, quadrature_for_space, ) @@ -90,6 +88,10 @@ LegendreGaussQuadrature, ) from modepy.quadrature.jaskowiec_sukumar import JaskowiecSukumarQuadrature +from modepy.quadrature.transplanted import ( + transplanted_1d_quadrature, + transplanted_legendre_gauss_quadrature, +) from modepy.quadrature.vioreanu_rokhlin import VioreanuRokhlinSimplexQuadrature from modepy.quadrature.witherden_vincent import WitherdenVincentQuadrature from modepy.quadrature.xiao_gimbutas import XiaoGimbutasSimplexQuadrature @@ -109,6 +111,8 @@ GaussLegendreQuadrature = LegendreGaussQuadrature +Transplanted1DQuadrature = transplanted_1d_quadrature +TransplantedLegendreGaussQuadrature = transplanted_legendre_gauss_quadrature __all__ = [ "PN", @@ -181,6 +185,8 @@ "submesh_for_shape", "symbolicize_function", "tensor_product_nodes", + "transplanted_1d_quadrature", + "transplanted_legendre_gauss_quadrature", "unit_vertices_for_shape", "vandermonde", "warp_and_blend_nodes", diff --git a/modepy/quadrature/__init__.py b/modepy/quadrature/__init__.py index 2a12503e..996f56fc 100644 --- a/modepy/quadrature/__init__.py +++ b/modepy/quadrature/__init__.py @@ -357,10 +357,4 @@ def _quadrature_for_tp( # }}} - -from modepy.quadrature.transplanted import ( - Transplanted1DQuadrature as Transplanted1DQuadrature, - TransplantedLegendreGaussQuadrature as TransplantedLegendreGaussQuadrature, -) - # vim: foldmethod=marker diff --git a/modepy/quadrature/transplanted.py b/modepy/quadrature/transplanted.py index c4ae4170..201fcf63 100644 --- a/modepy/quadrature/transplanted.py +++ b/modepy/quadrature/transplanted.py @@ -7,11 +7,11 @@ 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 +Given base nodes/weights :math:`(s_i, w_i^{(s)})`, the transplanted rule is .. math:: - x_i = g(s_i), \qquad \tilde w_i = w_i g'(s_i), + x_i = g(s_i), \qquad \tilde w_i = w_i^{(s)} g'(s_i), so that @@ -20,38 +20,8 @@ \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 +For map names, parameters, examples, and references, see +:ref:`quadrature-transplanted-1d`. """ from functools import lru_cache @@ -117,7 +87,7 @@ def map_identity(s: ArrayF) -> tuple[ArrayF, ArrayF]: Returns ``(s, 1)``. """ - return np.array(s, dtype=np.float64, copy=True), np.ones_like(s, dtype=np.float64) + return np.array(s, copy=True), np.ones_like(s) def _arcsin_taylor_coefficients(max_odd_degree: int) -> tuple[float, ...]: @@ -145,10 +115,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(sum(coeffs)) + denom = np.asarray(sum(coeffs), dtype=s.dtype) - g = np.zeros_like(s, dtype=np.float64) - gp = np.zeros_like(s, dtype=np.float64) + g = np.zeros_like(s) + gp = np.zeros_like(s) for k, c_k in enumerate(coeffs): power = 2 * k + 1 @@ -210,11 +180,25 @@ def map_kosloff_tal_ezer( denom = asin(alpha) 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) + gp = np.asarray( + alpha / (denom * np.sqrt(1.0 - alpha**2 * s**2)), + dtype=np.float64, + ) return g, gp +def _map_preserves_exact_to(map_name: str, *, sausage_degree: int) -> bool: + if map_name == "identity": + return True + + legacy_sausage_degree = _sausage_degree_from_map_name(map_name) + if legacy_sausage_degree is not None: + return legacy_sausage_degree == 1 + + return map_name == "sausage" and sausage_degree == 1 + + def _sausage_degree_from_map_name(map_name: str) -> int | None: if not map_name.startswith("sausage_d"): return None @@ -225,22 +209,7 @@ def _sausage_degree_from_map_name(map_name: str) -> int | None: 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 _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 + return int(degree_text) @lru_cache(maxsize=16) @@ -321,6 +290,7 @@ def map_trefethen_transplant( s: ArrayF, map_name: str, *, + sausage_degree: int = 9, strip_rho: float = 1.4, kte_rho: float = 1.4, kte_alpha: float | None = None, @@ -328,8 +298,9 @@ def map_trefethen_transplant( """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``, + :arg map_name: one of ``identity``, ``sausage``, ``kte``, ``kosloff_tal_ezer``, ``strip``. + :arg sausage_degree: odd polynomial degree for ``map_name="sausage"``. :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. @@ -340,28 +311,22 @@ def map_trefethen_transplant( The supported maps are: * ``identity``: :func:`map_identity` - * ``sausage_d{odd}``: :func:`map_sausage` + * ``sausage``: :func:`map_sausage` + * ``sausage_d{odd}`` (legacy alias): :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 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," 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: + if map_name == "sausage": return map_sausage(s, sausage_degree) + legacy_sausage_degree = _sausage_degree_from_map_name(map_name) + if legacy_sausage_degree is not None: + return map_sausage(s, legacy_sausage_degree) + if map_name == "strip": return map_strip(s, rho=strip_rho) @@ -371,11 +336,19 @@ def map_trefethen_transplant( raise ValueError( "unsupported map_name " f"'{map_name}'. Expected one of: " - "identity, sausage_d{odd}, kte, kosloff_tal_ezer, strip" + "identity, sausage, sausage_d{odd}, kte, kosloff_tal_ezer, strip" ) -class Transplanted1DQuadrature(Quadrature): +def transplanted_1d_quadrature( + quadrature: Quadrature, + map_name: str = "sausage", + *, + sausage_degree: int = 9, + strip_rho: float = 1.4, + kte_rho: float = 1.4, + kte_alpha: float | None = None, +) -> Quadrature: r"""Map an existing 1D quadrature rule using a Trefethen transplant map. The transformed rule approximates @@ -386,92 +359,65 @@ class Transplanted1DQuadrature(Quadrature): 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 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," Journal of - Computational Physics 104(2), 457-469 (1993), - 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, - 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 = np.asarray(base_nodes, dtype=np.float64) - force_dim_axis = False - elif base_nodes.ndim == 2 and base_nodes.shape[0] == 1: - nodes_1d = np.asarray(base_nodes[0], dtype=np.float64) - 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, + base_nodes = quadrature.nodes + if base_nodes.ndim == 1: + nodes_1d = np.asarray(base_nodes) + force_dim_axis = False + elif base_nodes.ndim == 2 and base_nodes.shape[0] == 1: + nodes_1d = cast("ArrayF", base_nodes[0]) + force_dim_axis = True + else: + raise ValueError( + "transplanted_1d_quadrature requires a one-dimensional base quadrature" ) - mapped_weights = quadrature.weights * jacobian - if force_dim_axis: - mapped_nodes = np.reshape(mapped_nodes, (1, mapped_nodes.shape[0])) + mapped_nodes, jacobian = map_trefethen_transplant( + nodes_1d, + map_name=map_name, + sausage_degree=sausage_degree, + strip_rho=strip_rho, + kte_rho=kte_rho, + kte_alpha=kte_alpha, + ) + mapped_weights = quadrature.weights * jacobian + + if force_dim_axis: + mapped_nodes = np.reshape(mapped_nodes, (1, mapped_nodes.shape[0])) - exact_to = quadrature._exact_to if _map_preserves_exact_to(map_name) else None - super().__init__(mapped_nodes, mapped_weights, exact_to=exact_to) + exact_to = None + if _map_preserves_exact_to(map_name, sausage_degree=sausage_degree): + try: + exact_to = quadrature.exact_to + except AttributeError: + exact_to = None - self.base_quadrature = quadrature - self.map_name = map_name - self.strip_rho = strip_rho - self.kte_rho = kte_rho - self.kte_alpha = kte_alpha + return Quadrature(mapped_nodes, mapped_weights, exact_to=exact_to) -class TransplantedLegendreGaussQuadrature(Transplanted1DQuadrature): +def transplanted_legendre_gauss_quadrature( + n: int, + map_name: str = "sausage", + *, + sausage_degree: int = 9, + strip_rho: float = 1.4, + kte_rho: float = 1.4, + kte_alpha: float | None = None, + backend: str | None = None, + force_dim_axis: bool = False, +) -> Quadrature: r"""Legendre-Gauss quadrature transplanted by a Trefethen map.""" + from modepy.quadrature.jacobi_gauss import LegendreGaussQuadrature - 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, - ) + return transplanted_1d_quadrature( + LegendreGaussQuadrature( + n, + backend=backend, + force_dim_axis=force_dim_axis, + ), + map_name=map_name, + sausage_degree=sausage_degree, + strip_rho=strip_rho, + kte_rho=kte_rho, + kte_alpha=kte_alpha, + ) From 74b8d51d7c41bfb67f6d01a52b06dbe75a3e1897 Mon Sep 17 00:00:00 2001 From: xywei Date: Wed, 25 Feb 2026 09:40:14 -0600 Subject: [PATCH 08/11] fix: test regex, simply kte --- modepy/__init__.py | 4 ---- modepy/quadrature/transplanted.py | 32 +++++++++++++------------------ 2 files changed, 13 insertions(+), 23 deletions(-) diff --git a/modepy/__init__.py b/modepy/__init__.py index 40194029..de1ee5b2 100644 --- a/modepy/__init__.py +++ b/modepy/__init__.py @@ -111,8 +111,6 @@ GaussLegendreQuadrature = LegendreGaussQuadrature -Transplanted1DQuadrature = transplanted_1d_quadrature -TransplantedLegendreGaussQuadrature = transplanted_legendre_gauss_quadrature __all__ = [ "PN", @@ -143,8 +141,6 @@ "TensorProductQuadrature", "TensorProductShape", "TensorProductSpace", - "Transplanted1DQuadrature", - "TransplantedLegendreGaussQuadrature", "VioreanuRokhlinSimplexQuadrature", "WitherdenVincentQuadrature", "XiaoGimbutasSimplexQuadrature", diff --git a/modepy/quadrature/transplanted.py b/modepy/quadrature/transplanted.py index 201fcf63..6189e7e2 100644 --- a/modepy/quadrature/transplanted.py +++ b/modepy/quadrature/transplanted.py @@ -92,7 +92,9 @@ def map_identity(s: ArrayF) -> tuple[ArrayF, 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}") + raise ValueError( + f"sausage must use positive odd degree, got: {max_odd_degree}" + ) nterms = (max_odd_degree + 1) // 2 coeffs = [1.0] @@ -128,23 +130,6 @@ def map_sausage(s: ArrayF, degree: int) -> tuple[ArrayF, ArrayF]: 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, *, @@ -176,7 +161,16 @@ def map_kosloff_tal_ezer( Computational Physics 104(2), 457-469 (1993), doi:10.1006/jcph.1993.1044. """ - alpha = _kte_alpha(rho=rho, alpha=alpha) + if alpha is None: + if rho <= 1.0: + raise ValueError(f"KTE parameter rho must be > 1: {rho}") + + alpha = float(2.0 / (rho + 1.0 / rho)) + + if not (0.0 < alpha < 1.0) or not isfinite(alpha): + raise ValueError(f"KTE parameter alpha must satisfy 0 < alpha < 1: {alpha}") + + alpha = float(alpha) denom = asin(alpha) g = np.asarray(np.arcsin(alpha * s) / denom, dtype=np.float64) From 7aa08e39a051252e471314ef0aeff677b0d29ce7 Mon Sep 17 00:00:00 2001 From: xywei Date: Wed, 25 Feb 2026 09:43:45 -0600 Subject: [PATCH 09/11] fix: tests --- modepy/test/test_quadrature.py | 134 +++++++++++++++++++++------------ 1 file changed, 85 insertions(+), 49 deletions(-) diff --git a/modepy/test/test_quadrature.py b/modepy/test/test_quadrature.py index 6ddfaf4a..82245800 100644 --- a/modepy/test/test_quadrature.py +++ b/modepy/test/test_quadrature.py @@ -46,8 +46,9 @@ def test_transformed_quadrature(): def gaussian_density(x, mu, sigma): return ( - 1 / (sigma * np.sqrt(2*np.pi)) - * np.exp(-np.sum((x-mu)**2, axis=0) / (2 * sigma**2)) + 1 + / (sigma * np.sqrt(2 * np.pi)) + * np.exp(-np.sum((x - mu) ** 2, axis=0) / (2 * sigma**2)) ) from modepy.quadrature import Transformed1DQuadrature @@ -56,8 +57,10 @@ def gaussian_density(x, mu, sigma): mu = 17 sigma = 12 tq = Transformed1DQuadrature( - LegendreGaussQuadrature(20, force_dim_axis=True), - left=mu - 6*sigma, right=mu + 6*sigma) + LegendreGaussQuadrature(20, force_dim_axis=True), + left=mu - 6 * sigma, + right=mu + 6 * sigma, + ) result = tq(lambda x: gaussian_density(x, mu, sigma)) assert abs(result - 1) < 1.0e-9 @@ -72,10 +75,13 @@ def gaussian_density(x, mu, sigma): @pytest.mark.parametrize("backend", BACKENDS) -@pytest.mark.parametrize("quad_type", [ - mp.LegendreGaussQuadrature, - mp.LegendreGaussLobattoQuadrature, - ]) +@pytest.mark.parametrize( + "quad_type", + [ + mp.LegendreGaussQuadrature, + mp.LegendreGaussLobattoQuadrature, + ], +) def test_gauss_quadrature(backend, quad_type): for s in range(9 + 1): if quad_type == mp.LegendreGaussLobattoQuadrature and s == 0: @@ -84,13 +90,14 @@ def test_gauss_quadrature(backend, quad_type): quad = quad_type(s, backend=backend, force_dim_axis=True) - assert quad.nodes.shape[1] == s+1 + assert quad.nodes.shape[1] == s + 1 for deg in range(quad.exact_to + 1): + def f(x): return np.sum(x**deg, axis=0) # noqa: B023 i_f = quad(f) - i_f_true = 1 / (deg+1) * (1 - (-1)**(deg + 1)) + i_f_true = 1 / (deg + 1) * (1 - (-1) ** (deg + 1)) err = abs(i_f - i_f_true) assert err < 3.0e-15, (s, deg, err, i_f, i_f_true) @@ -100,30 +107,38 @@ def test_clenshaw_curtis_quadrature() -> None: for s in range(1, 9 + 1): quad = ClenshawCurtisQuadrature(s, force_dim_axis=True) - assert quad.nodes.shape[1] == s+1 + assert quad.nodes.shape[1] == s + 1 for deg in range(quad.exact_to + 1): + def f(x): return x**deg # noqa: B023 i_f = quad(f) - i_f_true = 1 / (deg+1) * (1 - (-1)**(deg + 1)) + i_f_true = 1 / (deg + 1) * (1 - (-1) ** (deg + 1)) err = abs(i_f - i_f_true) 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: +@pytest.mark.parametrize( + ("map_name", "sausage_degree"), + [ + ("identity", 9), + ("sausage", 5), + ("sausage", 9), + ("sausage", 17), + ("kte", 9), + ("kosloff_tal_ezer", 9), + ], +) +def test_transplanted_legendre_gauss_quadrature( + map_name: str, + sausage_degree: int, +) -> None: base = mp.LegendreGaussQuadrature(15, force_dim_axis=True) - transplanted = mp.TransplantedLegendreGaussQuadrature( + transplanted = mp.transplanted_legendre_gauss_quadrature( 15, map_name=map_name, + sausage_degree=sausage_degree, force_dim_axis=True, ) @@ -138,17 +153,17 @@ def test_transplanted_legendre_gauss_quadrature(map_name: str) -> None: from modepy.quadrature.transplanted import map_trefethen_transplant mapped_nodes, mapped_jacobian = map_trefethen_transplant( - base_nodes, map_name=map_name + base_nodes, + map_name=map_name, + sausage_degree=sausage_degree, ) assert la.norm(trans_nodes - mapped_nodes, np.inf) < 1.0e-14 - assert ( - la.norm(trans_weights - base_weights * mapped_jacobian, np.inf) < 1.0e-14 - ) + assert 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"): + if map_name == "sausage": err = abs(float(np.sum(trans_weights)) - 2.0) assert err < 1.0e-14 @@ -156,7 +171,7 @@ def test_transplanted_legendre_gauss_quadrature(map_name: str) -> None: def test_transplanted_strip_map_quadrature() -> None: pytest.importorskip("scipy") - transplanted = mp.TransplantedLegendreGaussQuadrature( + transplanted = mp.transplanted_legendre_gauss_quadrature( 32, map_name="strip", strip_rho=1.4, @@ -175,15 +190,18 @@ def test_transplanted_strip_map_rejects_endpoint_rules() -> None: pytest.importorskip("scipy") with pytest.raises(ValueError, match="interior nodes"): - mp.Transplanted1DQuadrature( + mp.transplanted_1d_quadrature( 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 + mp.transplanted_legendre_gauss_quadrature( + 8, + map_name="sausage", + sausage_degree=4, + force_dim_axis=True, ) @@ -194,7 +212,11 @@ def test_transplanted_sausage_map_name_matches_direct_map() -> None: 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}") + g_ref, gp_ref = map_trefethen_transplant( + s, + map_name="sausage", + sausage_degree=degree, + ) assert la.norm(g - g_ref, np.inf) < 1.0e-15 assert la.norm(gp - gp_ref, np.inf) < 1.0e-15 @@ -243,17 +265,20 @@ def f(x): return x**deg # noqa: B023 i_f = quad(f) - i_f_true = 1 / (deg+1) * (1 - (-1)**(deg + 1)) + i_f_true = 1 / (deg + 1) * (1 - (-1) ** (deg + 1)) err = abs(i_f - i_f_true) assert err < 2.0e-15, (s, deg, err, i_f, i_f_true) -@pytest.mark.parametrize(("quad_class", "highest_order"), [ - (mp.XiaoGimbutasSimplexQuadrature, None), - (mp.JaskowiecSukumarQuadrature, None), - (mp.VioreanuRokhlinSimplexQuadrature, None), - (mp.GrundmannMoellerSimplexQuadrature, 3), - ]) +@pytest.mark.parametrize( + ("quad_class", "highest_order"), + [ + (mp.XiaoGimbutasSimplexQuadrature, None), + (mp.JaskowiecSukumarQuadrature, None), + (mp.VioreanuRokhlinSimplexQuadrature, None), + (mp.GrundmannMoellerSimplexQuadrature, 3), + ], +) @pytest.mark.parametrize("dim", [2, 3]) def test_simplex_quadrature(quad_class, highest_order, dim) -> None: """Check that quadratures on simplices works as advertised""" @@ -279,6 +304,7 @@ def test_simplex_quadrature(quad_class, highest_order, dim) -> None: if 0: import matplotlib.pyplot as pt + pt.plot(quad.nodes[0], quad.nodes[1]) pt.show() @@ -297,10 +323,13 @@ def test_simplex_quadrature(quad_class, highest_order, dim) -> None: @pytest.mark.parametrize("dim", [2, 3]) -@pytest.mark.parametrize(("quad_class", "max_order"), [ - (mp.WitherdenVincentQuadrature, np.inf), - (mp.LegendreGaussTensorProductQuadrature, 6), - ]) +@pytest.mark.parametrize( + ("quad_class", "max_order"), + [ + (mp.WitherdenVincentQuadrature, np.inf), + (mp.LegendreGaussTensorProductQuadrature, 6), + ], +) def test_hypercube_quadrature(dim, quad_class, max_order): from pytools import ( generate_nonnegative_integer_tuples_summing_to_at_most as gnitstam, @@ -314,8 +343,9 @@ def _check_monomial(quad, comb): int_exact = 2**dim * f.hypercube_integral() error = abs(int_approx - int_exact) / abs(int_exact) - logger.info("%s: %.5e %.5e / rel error %.5e", - comb, int_approx, int_exact, error) + logger.info( + "%s: %.5e %.5e / rel error %.5e", comb, int_approx, int_exact, error + ) return error @@ -329,8 +359,9 @@ def _check_monomial(quad, comb): assert np.all(quad.weights > 0) - logger.info("quadrature: %s %d %d", - quad_class.__name__.lower(), order, quad.exact_to) + logger.info( + "quadrature: %s %d %d", quad_class.__name__.lower(), order, quad.exact_to + ) for comb in gnitstam(quad.exact_to, dim): assert _check_monomial(quad, comb) < 5.0e-15 @@ -354,6 +385,7 @@ def test_mass_quadrature_is_newton_cotes(shape: mp.Shape, order: int) -> None: assert basis.orthonormality_weight() == 1 from math import factorial + shape_volume = 2**shape.dim / factorial(shape.dim) # integrals are orthogonal to the constant @@ -363,8 +395,10 @@ def test_mass_quadrature_is_newton_cotes(shape: mp.Shape, order: int) -> None: newton_cotes_weights = la.solve(vdm.T, integrals) - assert (la.norm(newton_cotes_weights - mass_weights, np.inf) - / la.norm(newton_cotes_weights, np.inf)) < 1e-13 + assert ( + la.norm(newton_cotes_weights - mass_weights, np.inf) + / la.norm(newton_cotes_weights, np.inf) + ) < 1e-13 # You can test individual routines by typing @@ -372,10 +406,12 @@ def test_mass_quadrature_is_newton_cotes(shape: mp.Shape, order: int) -> None: if __name__ == "__main__": import sys + if len(sys.argv) > 1: exec(sys.argv[1]) else: from pytest import main + main([__file__]) # vim: fdm=marker From 200d437f9ecdf8e40a681df813fe37e978fc7713 Mon Sep 17 00:00:00 2001 From: xywei Date: Wed, 25 Feb 2026 09:45:42 -0600 Subject: [PATCH 10/11] fix: remove heavy example --- examples/plot-qbx-transplanted-vs-gauss.py | 300 --------------------- 1 file changed, 300 deletions(-) delete mode 100644 examples/plot-qbx-transplanted-vs-gauss.py diff --git a/examples/plot-qbx-transplanted-vs-gauss.py b/examples/plot-qbx-transplanted-vs-gauss.py deleted file mode 100644 index b7c45764..00000000 --- a/examples/plot-qbx-transplanted-vs-gauss.py +++ /dev/null @@ -1,300 +0,0 @@ -from __future__ import annotations - -# 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 -import numpy as np -import pyopencl as cl -import pyopencl.tools as cl_tools -from arraycontext import flatten -from meshmode.discretization import Discretization -from meshmode.discretization.poly_element import ( - PolynomialGivenNodesElementGroup, -) -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 -ASSOC_TOL = 0.05 -NPTS = list(range(4, 30)) -MAPS = [ - ("gauss", None, None), - ("kte", "kte", None), - ("strip", "strip", None), - ("sausage_d5", "sausage", 5), - ("sausage_d9", "sausage", 9), -] -OUT = Path(tempfile.gettempdir()) / "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, - sausage_degree: int | None = None, -) -> mp.Quadrature: - if map_name is None: - return mp.LegendreGaussQuadrature(npts - 1, force_dim_axis=True) - return mp.transplanted_legendre_gauss_quadrature( - npts - 1, - map_name=map_name, - sausage_degree=9 if sausage_degree is None else sausage_degree, - 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 pairwise(panel_edges) - ]) - 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( - "auto strip rho: " - f"eta_min={eta_min:.6f}, " - f"target_half_width={target_half_width:.6f}, " - f"rho={rho:.4f}" - ) - return rho - - -def eval_rule( - actx, - lpot, - panel_edges: np.ndarray, - npts: int, - map_name: str | None, - sausage_degree: int | None, - strip_rho: float, - targets: np.ndarray, - centers: np.ndarray, - radii: np.ndarray, -) -> np.ndarray: - quad = make_quad(npts, map_name, strip_rho, sausage_degree) - 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( - "using strip " - f"rho={strip_rho:.4f} " - f"(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, sausage_degree in MAPS: - values = eval_rule( - actx, - lpot, - panel_edges, - npts, - map_name, - sausage_degree, - 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() From 89886062d7fa971791f7e2764b621c79312224c0 Mon Sep 17 00:00:00 2001 From: xywei Date: Wed, 25 Feb 2026 10:36:00 -0600 Subject: [PATCH 11/11] fix: test typing --- modepy/test/test_quadrature.py | 32 +++++++++++++++++++------------- 1 file changed, 19 insertions(+), 13 deletions(-) diff --git a/modepy/test/test_quadrature.py b/modepy/test/test_quadrature.py index 82245800..ec7907f9 100644 --- a/modepy/test/test_quadrature.py +++ b/modepy/test/test_quadrature.py @@ -44,11 +44,11 @@ def test_transformed_quadrature(): """Test 1D quadrature on arbitrary intervals""" - def gaussian_density(x, mu, sigma): - return ( - 1 - / (sigma * np.sqrt(2 * np.pi)) - * np.exp(-np.sum((x - mu) ** 2, axis=0) / (2 * sigma**2)) + def gaussian_density(x: ArrayF, mu: float, sigma: float) -> ArrayF: + sq_dist = cast("ArrayF", np.sum((x - mu) ** 2, axis=0)) + return cast( + "ArrayF", + 1 / (sigma * np.sqrt(2 * np.pi)) * np.exp(-sq_dist / (2 * sigma**2)), ) from modepy.quadrature import Transformed1DQuadrature @@ -337,10 +337,14 @@ def test_hypercube_quadrature(dim, quad_class, max_order): from modepy.tools import Monomial - def _check_monomial(quad, comb): + def _check_monomial(quad: mp.Quadrature, comb: tuple[int, ...]) -> float: f = Monomial(comb) - int_approx = quad(f) - int_exact = 2**dim * f.hypercube_integral() + int_approx_obj = quad(f) + if isinstance(int_approx_obj, np.ndarray): + int_approx = float(int_approx_obj.item()) + else: + int_approx = float(int_approx_obj) + int_exact = cast("float", 2**dim * f.hypercube_integral()) error = abs(int_approx - int_exact) / abs(int_exact) logger.info( @@ -357,15 +361,15 @@ def _check_monomial(quad, comb): logger.info("UNAVAILABLE at order %d", order) break + quad_exact_to = cast("int", quad.exact_to) + assert np.all(quad.weights > 0) - logger.info( - "quadrature: %s %d %d", quad_class.__name__.lower(), order, quad.exact_to - ) - for comb in gnitstam(quad.exact_to, dim): + logger.info("quadrature: order %d %d", order, quad_exact_to) + for comb in gnitstam(quad_exact_to, dim): assert _check_monomial(quad, comb) < 5.0e-15 - comb = (0,) * (dim - 1) + (quad.exact_to + 1,) + comb = (0,) * (dim - 1) + (quad_exact_to + 1,) assert _check_monomial(quad, comb) > 5.0e-15 order += 2 @@ -394,6 +398,8 @@ def test_mass_quadrature_is_newton_cotes(shape: mp.Shape, order: int) -> None: integrals[0] = shape_volume * basis.functions[0](np.zeros(shape.dim)) newton_cotes_weights = la.solve(vdm.T, integrals) + mass_weights = cast("ArrayF", mass_weights) + newton_cotes_weights = cast("ArrayF", newton_cotes_weights) assert ( la.norm(newton_cotes_weights - mass_weights, np.inf)