From 8c44d1649bb0f6c662264226def1c0bddf0c7bf7 Mon Sep 17 00:00:00 2001 From: Andreas Kloeckner Date: Thu, 6 Jun 2024 14:01:10 -0500 Subject: [PATCH 01/12] Fix, test node count claims for Gauss quadrature --- modepy/nodes.py | 6 ++++-- modepy/quadrature/jacobi_gauss.py | 13 +++++++------ modepy/test/test_quadrature.py | 2 ++ 3 files changed, 13 insertions(+), 8 deletions(-) diff --git a/modepy/nodes.py b/modepy/nodes.py index 8bffe2d..04b49a7 100644 --- a/modepy/nodes.py +++ b/modepy/nodes.py @@ -395,7 +395,8 @@ def tensor_product_nodes( def legendre_gauss_tensor_product_nodes(dims: int, n: int) -> np.ndarray: """ - :arg n: the number of points in one dimension. + :arg n: the degree of polynomial exactly interpolated by the nodes. + The one-dimensional base quadrature has *n+1* nodes. .. versionadded:: 2024.2 """ @@ -408,7 +409,8 @@ def legendre_gauss_tensor_product_nodes(dims: int, n: int) -> np.ndarray: def legendre_gauss_lobatto_tensor_product_nodes(dims: int, n: int) -> np.ndarray: """ - :arg n: the number of points in one dimension. + :arg n: the degree of polynomial exactly interpolated by the nodes. + The one-dimensional base quadrature has *n+1* nodes. """ from modepy.quadrature.jacobi_gauss import legendre_gauss_lobatto_nodes return tensor_product_nodes(dims, diff --git a/modepy/quadrature/jacobi_gauss.py b/modepy/quadrature/jacobi_gauss.py index 17e440c..b1ca70e 100644 --- a/modepy/quadrature/jacobi_gauss.py +++ b/modepy/quadrature/jacobi_gauss.py @@ -38,7 +38,7 @@ class JacobiGaussQuadrature(Quadrature): I[f] = \int_{-1}^1 f(x) (1 - x)^\alpha (1 + x)^\beta\, \mathrm{d}x, where :math:`\alpha, \beta > -1`. The quadrature rule is exact - up to degree :math:`2N + 1`. + up to degree :math:`2N + 1`. The quadrature has *N+1* nodes. .. automethod:: __init__ """ @@ -167,7 +167,7 @@ def b(n): class LegendreGaussQuadrature(JacobiGaussQuadrature): - r"""A Gauss quadrature rule with weight :math:`1`. + r"""A Gauss quadrature rule with weight :math:`1` and *N+1* nodes. Corresponds to a Gauss-Jacobi quadrature rule with :math:`\alpha = \beta = 0`. @@ -183,7 +183,8 @@ def __init__(self, class ChebyshevGaussQuadrature(JacobiGaussQuadrature): - r"""A Gauss quadrature rule with weight :math:`\sqrt{1-x^2}^{\mp 1}`. + r"""A Gauss quadrature rule with weight :math:`\sqrt{1-x^2}^{\mp 1}` + and *N+1* nodes. The Chebyshev-Gauss quadrature rules of the first kind and second kind correspond to Gauss-Jacobi quadrature rules with @@ -212,7 +213,7 @@ def __init__(self, class GaussGegenbauerQuadrature(JacobiGaussQuadrature): r"""Gauss-Gegenbauer quadrature is a special case of Gauss-Jacobi quadrature - with :math:`\alpha = \beta`. + with :math:`\alpha = \beta` and *N+1* nodes. .. versionadded:: 2019.1 """ @@ -232,7 +233,7 @@ def jacobi_gauss_lobatto_nodes( force_dim_axis: bool = False) -> np.ndarray: """Compute the Gauss-Lobatto quadrature nodes corresponding to the :class:`~modepy.JacobiGaussQuadrature` - with the same parameters. + with the same parameters. There will be *N+1* nodes. Exact to degree :math:`2N - 3`. """ @@ -261,7 +262,7 @@ def legendre_gauss_lobatto_nodes( backend: Optional[str] = None, force_dim_axis: bool = False) -> np.ndarray: """Compute the Legendre-Gauss-Lobatto quadrature nodes. - *N* is the number of points. + *N+1* is the number of nodes. Exact to degree :math:`2N - 1`. """ diff --git a/modepy/test/test_quadrature.py b/modepy/test/test_quadrature.py index 729f9b5..641cf58 100644 --- a/modepy/test/test_quadrature.py +++ b/modepy/test/test_quadrature.py @@ -68,6 +68,7 @@ def test_gauss_quadrature(backend): for s in range(9 + 1): quad = LegendreGaussQuadrature(s, backend, force_dim_axis=True) + 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 @@ -83,6 +84,7 @@ def test_clenshaw_curtis_quadrature(): for s in range(1, 9 + 1): quad = ClenshawCurtisQuadrature(s, force_dim_axis=True) + assert quad.nodes.shape[1] == s+1 for deg in range(quad.exact_to + 1): def f(x): return x**deg # noqa: B023 From 60a7be171a73f09075807982aecae9fd59f3c516 Mon Sep 17 00:00:00 2001 From: Andreas Kloeckner Date: Thu, 6 Jun 2024 14:01:42 -0500 Subject: [PATCH 02/12] Fix test_nodal_mass_matrix_against_quad: Do not overwrite nodes --- modepy/test/test_matrices.py | 6 ++---- 1 file changed, 2 insertions(+), 4 deletions(-) diff --git a/modepy/test/test_matrices.py b/modepy/test/test_matrices.py index 37f4d4e..e355258 100644 --- a/modepy/test/test_matrices.py +++ b/modepy/test/test_matrices.py @@ -46,10 +46,10 @@ def test_nodal_mass_matrix_against_quad( if isinstance(shape, mp.Hypercube) or shape == mp.Simplex(1): if nodes_on_bdry: nodes = mp.legendre_gauss_lobatto_tensor_product_nodes( - shape.dim, order + 1, + shape.dim, order, ) else: - nodes = mp.legendre_gauss_tensor_product_nodes(shape.dim, order + 1) + nodes = mp.legendre_gauss_tensor_product_nodes(shape.dim, order) elif isinstance(shape, mp.Simplex): if nodes_on_bdry: nodes = mp.warp_and_blend_nodes(shape.dim, order) @@ -59,8 +59,6 @@ def test_nodal_mass_matrix_against_quad( else: raise AssertionError() - nodes = mp.edge_clustered_nodes_for_space(space, shape) - basis = mp.orthonormal_basis_for_space(space, shape) quad_mass_mat = mp.nodal_quad_mass_matrix(quad, basis.functions, nodes) From 6fc8ad72883454db69ad4b9664af3863cbd8b9e4 Mon Sep 17 00:00:00 2001 From: Andreas Kloeckner Date: Thu, 6 Jun 2024 14:02:10 -0500 Subject: [PATCH 03/12] test_matrices: Add commandline test run snippet --- modepy/test/test_matrices.py | 14 ++++++++++++++ 1 file changed, 14 insertions(+) diff --git a/modepy/test/test_matrices.py b/modepy/test/test_matrices.py index e355258..28f2cb5 100644 --- a/modepy/test/test_matrices.py +++ b/modepy/test/test_matrices.py @@ -68,3 +68,17 @@ def test_nodal_mass_matrix_against_quad( err = la.norm(quad_mass_mat@nodes_to_quad - vdm_mass_mat)/la.norm(vdm_mass_mat) assert err < 1e-14 + + +# You can test individual routines by typing +# $ python test_modes.py 'test_routine()' + +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 487486e5df3de292ed342520ec43db91acba01a7 Mon Sep 17 00:00:00 2001 From: Andreas Kloeckner Date: Thu, 6 Jun 2024 14:52:39 -0500 Subject: [PATCH 04/12] Move Jacobi-Gauss quadrature docs into source file --- doc/quadrature.rst | 21 --------------------- modepy/quadrature/jacobi_gauss.py | 23 +++++++++++++++++++++++ 2 files changed, 23 insertions(+), 21 deletions(-) diff --git a/doc/quadrature.rst b/doc/quadrature.rst index 83cff9f..0ed5199 100644 --- a/doc/quadrature.rst +++ b/doc/quadrature.rst @@ -11,27 +11,6 @@ Jacobi-Gauss quadrature in one dimension .. automodule:: modepy.quadrature.jacobi_gauss -.. currentmodule:: modepy - -.. autoclass:: JacobiGaussQuadrature - :members: - :show-inheritance: - -.. autoclass:: LegendreGaussQuadrature - :show-inheritance: - -.. autoclass:: ChebyshevGaussQuadrature - :show-inheritance: - -.. autoclass:: GaussGegenbauerQuadrature - :show-inheritance: - -.. currentmodule:: modepy.quadrature.jacobi_gauss - -.. autofunction:: jacobi_gauss_lobatto_nodes - -.. autofunction:: legendre_gauss_lobatto_nodes - Clenshaw-Curtis and Fejér quadrature in one dimension ----------------------------------------------------- diff --git a/modepy/quadrature/jacobi_gauss.py b/modepy/quadrature/jacobi_gauss.py index b1ca70e..c7edae9 100644 --- a/modepy/quadrature/jacobi_gauss.py +++ b/modepy/quadrature/jacobi_gauss.py @@ -1,3 +1,26 @@ +""" +.. currentmodule:: modepy + +.. autoclass:: JacobiGaussQuadrature + :members: + :show-inheritance: + +.. autoclass:: LegendreGaussQuadrature + :show-inheritance: + +.. autoclass:: ChebyshevGaussQuadrature + :show-inheritance: + +.. autoclass:: GaussGegenbauerQuadrature + :show-inheritance: + +.. currentmodule:: modepy.quadrature.jacobi_gauss + +.. autofunction:: jacobi_gauss_lobatto_nodes + +.. autofunction:: legendre_gauss_lobatto_nodes +""" + __copyright__ = "Copyright (C) 2007 Andreas Kloeckner" __license__ = """ From ab73229cbe4d885b419ffc249f175579900314df Mon Sep 17 00:00:00 2001 From: Andreas Kloeckner Date: Fri, 7 Jun 2024 12:35:02 -0500 Subject: [PATCH 05/12] Add, test scaled_jacobi --- modepy/__init__.py | 4 +-- modepy/modes.py | 51 +++++++++++++++++++++++++++++++++++++++ modepy/test/test_modes.py | 21 ++++++++++++++++ 3 files changed, 74 insertions(+), 2 deletions(-) diff --git a/modepy/__init__.py b/modepy/__init__.py index 4943aaa..ade3513 100644 --- a/modepy/__init__.py +++ b/modepy/__init__.py @@ -30,7 +30,7 @@ nodal_quad_mass_matrix_for_face, resampling_matrix, vandermonde) from modepy.modes import ( Basis, BasisNotOrthonormal, TensorProductBasis, basis_for_space, grad_jacobi, - jacobi, monomial_basis_for_space, orthonormal_basis_for_space, + jacobi, monomial_basis_for_space, orthonormal_basis_for_space, scaled_jacobi, symbolicize_function) from modepy.nodes import ( edge_clustered_nodes_for_space, equidistant_nodes, equispaced_nodes_for_space, @@ -68,7 +68,7 @@ "FunctionSpace", "TensorProductSpace", "PN", "QN", "space_for_shape", - "jacobi", "grad_jacobi", + "jacobi", "grad_jacobi", "scaled_jacobi", "symbolicize_function", "Basis", "BasisNotOrthonormal", "TensorProductBasis", "basis_for_space", "orthonormal_basis_for_space", "monomial_basis_for_space", diff --git a/modepy/modes.py b/modepy/modes.py index c4ca59c..31b62f7 100644 --- a/modepy/modes.py +++ b/modepy/modes.py @@ -66,6 +66,7 @@ .. autofunction:: jacobi .. autofunction:: grad_jacobi +.. autofunction:: scaled_jacobi Conversion to Symbolic ---------------------- @@ -209,6 +210,56 @@ def grad_jacobi(alpha: float, beta: float, n: int, x: RealValueT) -> RealValueT: else: return math.sqrt(n*(n+alpha+beta+1)) * jacobi(alpha+1, beta+1, n-1, x) + +def binom(x, y): + # FIXME This may overflow quickly. + # mpmath has clever 'gammaprod' pole handling in case that's ever needed: + # https://github.com/mpmath/mpmath/blob/75a2ed37c4f2c576a9d01d360ee4c94ead57c7ff/mpmath/functions/factorials.py#L61-L65 + from math import gamma + return gamma(x + 1) / (gamma(y + 1) * gamma(x - y + 1)) + + +def scaled_jacobi(alpha: float, beta: float, n: int, x: RealValueT) -> RealValueT: + r""" + Same as :func:`jacobi`, but scaled so that + + .. math:: + + P_n^{(\alpha, \beta)}(1) = \binom{ n+\alpha \\ n}, + + and + + .. math:: + + P_n^{(\alpha, \beta)}(-1) = (-1)^n \binom{ n+\beta \\ n}. + + This is the more conventional definition of the Jacobi polynomials. + """ + from math import factorial, gamma, sqrt + + eps = 100 * np.finfo((np.float32(0) + np.asarray(x)).dtype).eps + if n == 0 and abs(alpha + beta + 1) < eps: + # maxima claims limit(1/(x*gamma(x)),x, 0) == 1 + scaling = ( + # 2**(alpha + beta + 1) \approx 1 + gamma(n + alpha + 1) + * gamma(n + beta + 1) + / factorial(n) + ) + else: + # source: + # https://en.wikipedia.org/w/index.php?title=Jacobi_polynomials&oldid=1219118998 + scaling = ( + 2**(alpha + beta + 1) + / (2*n + alpha + beta + 1) + * gamma(n + alpha + 1) + * gamma(n + beta + 1) + / gamma(n + alpha + beta + 1) + / factorial(n) + ) + + return sqrt(scaling) * jacobi(alpha, beta, n, x) + # }}} diff --git a/modepy/test/test_modes.py b/modepy/test/test_modes.py index 3ed9463..879bab6 100644 --- a/modepy/test/test_modes.py +++ b/modepy/test/test_modes.py @@ -37,6 +37,27 @@ logger = logging.getLogger(__name__) +@pytest.mark.parametrize(("alpha", "beta"), [ + (0, 0), # Gauss-Legendre + (-0.5, -0.5), # Chebyshev-Gauss (first kind) + (0.5, 0.5), # Chebyshev-Gauss (second kind) + (1, 0), + (3, 2), + (0, 2), + (5, 0), + (3, 4) + ]) +def test_scaled_jacobi(alpha, beta): + from modepy.modes import binom, scaled_jacobi as jacobi + + for n in range(10): + err1 = abs(jacobi(alpha, beta, n, 1) - binom(n + alpha, n)) + errm1 = abs(jacobi(alpha, beta, n, -1) - (-1)**n*binom(n + beta, n)) + + assert err1 < 1e-11 + assert errm1 < 1e-11 + + # {{{ test_orthonormality_jacobi_1d @pytest.mark.parametrize(("alpha", "beta", "ebound"), [ From 44b60adf3642a121983420a1908b957bfcf3b5d0 Mon Sep 17 00:00:00 2001 From: Andreas Kloeckner Date: Fri, 7 Jun 2024 12:43:06 -0500 Subject: [PATCH 06/12] Add JacobiGaussLobattoQuadrature, LegendreGaussLobattoQuadrature NB: The actual-Jacobi part of the quadrature is kind of untested right now. --- modepy/__init__.py | 8 ++- modepy/quadrature/jacobi_gauss.py | 94 +++++++++++++++++++++++++++++++ modepy/test/test_quadrature.py | 15 +++-- 3 files changed, 111 insertions(+), 6 deletions(-) diff --git a/modepy/__init__.py b/modepy/__init__.py index ade3513..d01a465 100644 --- a/modepy/__init__.py +++ b/modepy/__init__.py @@ -44,8 +44,9 @@ ClenshawCurtisQuadrature, FejerQuadrature) from modepy.quadrature.grundmann_moeller import GrundmannMoellerSimplexQuadrature from modepy.quadrature.jacobi_gauss import ( - ChebyshevGaussQuadrature, GaussGegenbauerQuadrature, JacobiGaussQuadrature, - LegendreGaussQuadrature) + ChebyshevGaussQuadrature, GaussGegenbauerQuadrature, + JacobiGaussLobattoQuadrature, JacobiGaussQuadrature, + LegendreGaussLobattoQuadrature, LegendreGaussQuadrature) from modepy.quadrature.jaskowiec_sukumar import JaskowiecSukumarQuadrature from modepy.quadrature.vioreanu_rokhlin import VioreanuRokhlinSimplexQuadrature from modepy.quadrature.witherden_vincent import WitherdenVincentQuadrature @@ -98,6 +99,9 @@ "JacobiGaussQuadrature", "LegendreGaussQuadrature", "GaussLegendreQuadrature", "ChebyshevGaussQuadrature", "GaussGegenbauerQuadrature", + "JacobiGaussLobattoQuadrature", + "LegendreGaussLobattoQuadrature", + "XiaoGimbutasSimplexQuadrature", "GrundmannMoellerSimplexQuadrature", "VioreanuRokhlinSimplexQuadrature", "ClenshawCurtisQuadrature", "FejerQuadrature", diff --git a/modepy/quadrature/jacobi_gauss.py b/modepy/quadrature/jacobi_gauss.py index c7edae9..dc3129c 100644 --- a/modepy/quadrature/jacobi_gauss.py +++ b/modepy/quadrature/jacobi_gauss.py @@ -19,6 +19,11 @@ .. autofunction:: jacobi_gauss_lobatto_nodes .. autofunction:: legendre_gauss_lobatto_nodes + +.. currentmodule:: modepy + +.. autoclass:: JacobiGaussLobattoQuadrature +.. autoclass:: LegendreGaussLobattoQuadrature """ __copyright__ = "Copyright (C) 2007 Andreas Kloeckner" @@ -291,3 +296,92 @@ def legendre_gauss_lobatto_nodes( """ return jacobi_gauss_lobatto_nodes(0, 0, N, backend=backend, force_dim_axis=force_dim_axis) + + +class JacobiGaussLobattoQuadrature(Quadrature): + r""" + Compute the Jacobi-Gauss-Lobatto quadrature with respect + to the Jacobi weight function + + .. math:: + + I[f] = \int_{-1}^1 f(x) (1 - x)^\alpha (1 + x)^\beta\, \mathrm{d}x, + + There will be *N+1* nodes. Exact to degree :math:`2N - 3`. + + .. versionadded:: 2024.2 + """ + + def __init__(self, + alpha: float, beta: float, N: int, # noqa: N803 + *, backend: Optional[str] = None, + force_dim_axis: bool = False) -> None: + nodes = jacobi_gauss_lobatto_nodes(alpha, beta, N, backend) + + from math import gamma + + from modepy.modes import binom, scaled_jacobi + + # Formula numbers refer to https://doi.org/10.1023/A:1016689830453 + if N + 1 < 2: + raise ValueError("Lobatto rules must have at least two nodes") + + # Alternate reference via chebfun, using same source paper but + # using the beta function in implementation: + # https://github.com/chebfun/chebfun/blob/db207bc9f48278ca4def15bf90591bfa44d0801d/lobpts.m + + n = N - 1 + + # FIXME: Recurrences might be better than just typing up the formulas. + + # (3.10) + common_factor = ( + 2**(alpha + beta + 1) + * binom(n + alpha + 1, n) + / binom(n + beta + 1, n) + / binom(n + alpha + beta + 2, n) + * gamma(alpha + 2) + / gamma(alpha + beta + 3) + ) + edge_weight = ( + common_factor + * gamma(beta + 1) + ) + + # (4.7) + int_nodes = nodes[1:-1] + frac = ( + 4*(n+alpha+1) * (n+beta+1) + (alpha-beta)**2 + ) / ( + 2*n + alpha + beta + 2 + )**2 + + int_weights = ( + common_factor + * gamma(beta + 2) + + # FIXME: This is part of the paper, but tests only pass without it: + # / (n + 1)**2 + + / (frac - int_nodes**2) + * (1 - int_nodes**2) + / scaled_jacobi(alpha, beta, n + 1, int_nodes)**2 + ) + + weights = np.empty_like(nodes) + weights[0] = edge_weight + weights[-1] = edge_weight + weights[1:-1] = int_weights + + if force_dim_axis: + nodes = nodes.reshape(1, -1) + + super().__init__(nodes, weights, 2*N - 3) + + +class LegendreGaussLobattoQuadrature(JacobiGaussLobattoQuadrature): + def __init__( + self, N, *, backend: Optional[str] = None, # noqa: N803 + force_dim_axis: bool = False + ) -> None: + super().__init__(0, 0, N, backend=backend, force_dim_axis=force_dim_axis) diff --git a/modepy/test/test_quadrature.py b/modepy/test/test_quadrature.py index 641cf58..d05951b 100644 --- a/modepy/test/test_quadrature.py +++ b/modepy/test/test_quadrature.py @@ -63,11 +63,18 @@ def gaussian_density(x, mu, sigma): @pytest.mark.parametrize("backend", BACKENDS) -def test_gauss_quadrature(backend): - from modepy.quadrature.jacobi_gauss import LegendreGaussQuadrature - +@pytest.mark.parametrize("quad_type", [ + mp.LegendreGaussQuadrature, + mp.LegendreGaussLobattoQuadrature, + ]) +def test_gauss_quadrature(backend, quad_type): for s in range(9 + 1): - quad = LegendreGaussQuadrature(s, backend, force_dim_axis=True) + if quad_type == mp.LegendreGaussLobattoQuadrature and s == 0: + # no one-node Lobatto rule + continue + + quad = quad_type(s, backend=backend, force_dim_axis=True) + assert quad.nodes.shape[1] == s+1 for deg in range(quad.exact_to + 1): def f(x): From 553e5258e9b789217558fb0c4dbe08d7d6f861b7 Mon Sep 17 00:00:00 2001 From: Andreas Kloeckner Date: Fri, 7 Jun 2024 12:44:15 -0500 Subject: [PATCH 07/12] Drop a stray period in jacobi_gauss --- modepy/quadrature/jacobi_gauss.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/modepy/quadrature/jacobi_gauss.py b/modepy/quadrature/jacobi_gauss.py index dc3129c..d8caec1 100644 --- a/modepy/quadrature/jacobi_gauss.py +++ b/modepy/quadrature/jacobi_gauss.py @@ -144,7 +144,7 @@ def compute_weights_and_nodes( if abs(alpha + 0.5) < 1.0e-14 and abs(beta + 0.5) < 1.0e-14: # NOTE: these are hardcoded for two reasons: # * the algorithm below doesn't work for alpha = beta = -0.5 - # * and, well, we know them explicitly so why not.. + # * and, well, we know them explicitly so why not. return ( np.cos(np.pi / (2 * (N+1)) * (2*np.arange(N+1, 0, -1) - 1)), np.full(N+1, np.pi / (N + 1)) From fd2d5f6bf165db9c4447524d615a8232e3155b8a Mon Sep 17 00:00:00 2001 From: Andreas Kloeckner Date: Fri, 7 Jun 2024 14:14:15 -0500 Subject: [PATCH 08/12] Add LegendreGaussLobattoTensorProductQuadrature --- modepy/__init__.py | 2 ++ modepy/quadrature/__init__.py | 14 ++++++++++++++ 2 files changed, 16 insertions(+) diff --git a/modepy/__init__.py b/modepy/__init__.py index d01a465..a54c9cd 100644 --- a/modepy/__init__.py +++ b/modepy/__init__.py @@ -38,6 +38,7 @@ node_tuples_for_space, random_nodes_for_shape, tensor_product_nodes, warp_and_blend_nodes) from modepy.quadrature import ( + LegendreGaussLobattoTensorProductQuadrature, LegendreGaussTensorProductQuadrature, Quadrature, QuadratureRuleUnavailable, TensorProductQuadrature, ZeroDimensionalQuadrature, quadrature_for_space) from modepy.quadrature.clenshaw_curtis import ( @@ -94,6 +95,7 @@ "Quadrature", "QuadratureRuleUnavailable", "ZeroDimensionalQuadrature", "TensorProductQuadrature", "LegendreGaussTensorProductQuadrature", + "LegendreGaussLobattoTensorProductQuadrature", "quadrature_for_space", "JacobiGaussQuadrature", "LegendreGaussQuadrature", diff --git a/modepy/quadrature/__init__.py b/modepy/quadrature/__init__.py index 690f8bf..7637c3c 100644 --- a/modepy/quadrature/__init__.py +++ b/modepy/quadrature/__init__.py @@ -182,6 +182,20 @@ def __init__(self, ] * dims) +class LegendreGaussLobattoTensorProductQuadrature(TensorProductQuadrature): + """A tensor product using only :class:`~modepy.LegendreGaussLobattoQuadrature` + one-dimenisonal rules. + """ + + def __init__(self, + N: int, dims: int, # noqa: N803 + backend: Optional[str] = None) -> None: + from modepy.quadrature.jacobi_gauss import LegendreGaussLobattoQuadrature + super().__init__([ + LegendreGaussLobattoQuadrature(N, backend=backend, force_dim_axis=True) + ] * dims) + + # {{{ quadrature @singledispatch From f516b750c24b4f1dd719ae32e0b09c5afde1d871 Mon Sep 17 00:00:00 2001 From: Andreas Kloeckner Date: Fri, 7 Jun 2024 14:19:19 -0500 Subject: [PATCH 09/12] TensorProductQuadrature: expose constituent quadrature rules --- modepy/matrices.py | 3 ++- modepy/quadrature/__init__.py | 12 ++++++++++-- 2 files changed, 12 insertions(+), 3 deletions(-) diff --git a/modepy/matrices.py b/modepy/matrices.py index c9e35b5..74c0ec0 100644 --- a/modepy/matrices.py +++ b/modepy/matrices.py @@ -28,7 +28,7 @@ import numpy.linalg as la from modepy.modes import Basis, BasisNotOrthonormal -from modepy.quadrature import Quadrature +from modepy.quadrature import Quadrature, TensorProductQuadrature from modepy.shapes import Face @@ -62,6 +62,7 @@ .. autofunction:: nodal_mass_matrix_for_face .. autofunction:: modal_quad_mass_matrix_for_face .. autofunction:: nodal_quad_mass_matrix_for_face +.. autofunction:: spectral_diag_nodal_mass_matrix Differentiation is also convenient to express by using :math:`V^{-1}` to obtain modal values and then using a Vandermonde matrix for the derivatives diff --git a/modepy/quadrature/__init__.py b/modepy/quadrature/__init__.py index 7637c3c..d6f6e32 100644 --- a/modepy/quadrature/__init__.py +++ b/modepy/quadrature/__init__.py @@ -35,7 +35,7 @@ """ from functools import singledispatch -from typing import Callable, Iterable, Optional +from typing import Callable, Iterable, Optional, Sequence import numpy as np @@ -143,9 +143,15 @@ def __init__(self, quad: Quadrature, left: float, right: float) -> None: class TensorProductQuadrature(Quadrature): r"""A tensor product quadrature of one-dimensional :class:`Quadrature`\ s. + .. autoattribute:: quadratures .. automethod:: __init__ """ + quadratures: Sequence[Quadrature] + """The lower-dimensional quadratures from which the tensor product quadrature is + composed. + """ + def __init__(self, quads: Iterable[Quadrature]) -> None: """ :arg quad: a iterable of :class:`Quadrature` objects for one-dimensional @@ -154,7 +160,7 @@ def __init__(self, quads: Iterable[Quadrature]) -> None: from modepy.nodes import tensor_product_nodes - quads = list(quads) + quads = tuple(quads) x = tensor_product_nodes([quad.nodes for quad in quads]) w = np.prod(tensor_product_nodes([quad.weights for quad in quads]), axis=0) assert w.size == x.shape[1] @@ -167,6 +173,8 @@ def __init__(self, quads: Iterable[Quadrature]) -> None: super().__init__(x, w, exact_to=exact_to) + self.quadratures = quads + class LegendreGaussTensorProductQuadrature(TensorProductQuadrature): """A tensor product using only :class:`~modepy.LegendreGaussQuadrature` From e86497430738c08146c9a93f4381c0a4991d45b3 Mon Sep 17 00:00:00 2001 From: Andreas Kloeckner Date: Fri, 7 Jun 2024 14:20:42 -0500 Subject: [PATCH 10/12] Add, test spectral_diag_nodal_mass_matrix --- modepy/__init__.py | 4 +++- modepy/matrices.py | 21 +++++++++++++++++ modepy/test/test_matrices.py | 45 ++++++++++++++++++++++++++++++++++++ 3 files changed, 69 insertions(+), 1 deletion(-) diff --git a/modepy/__init__.py b/modepy/__init__.py index a54c9cd..11e4047 100644 --- a/modepy/__init__.py +++ b/modepy/__init__.py @@ -27,7 +27,8 @@ inverse_mass_matrix, mass_matrix, modal_mass_matrix_for_face, modal_quad_mass_matrix, modal_quad_mass_matrix_for_face, multi_vandermonde, nodal_mass_matrix_for_face, nodal_quad_mass_matrix, - nodal_quad_mass_matrix_for_face, resampling_matrix, vandermonde) + nodal_quad_mass_matrix_for_face, resampling_matrix, + spectral_diag_nodal_mass_matrix, vandermonde) from modepy.modes import ( Basis, BasisNotOrthonormal, TensorProductBasis, basis_for_space, grad_jacobi, jacobi, monomial_basis_for_space, orthonormal_basis_for_space, scaled_jacobi, @@ -88,6 +89,7 @@ "diff_matrix_permutation", "inverse_mass_matrix", "mass_matrix", "modal_quad_mass_matrix", "nodal_quad_mass_matrix", + "spectral_diag_nodal_mass_matrix", "modal_mass_matrix_for_face", "nodal_mass_matrix_for_face", "modal_quad_mass_matrix_for_face", "nodal_quad_mass_matrix_for_face", diff --git a/modepy/matrices.py b/modepy/matrices.py index 74c0ec0..c730435 100644 --- a/modepy/matrices.py +++ b/modepy/matrices.py @@ -397,6 +397,27 @@ def nodal_quad_mass_matrix( modal_quad_mass_matrix(quadrature, test_functions)) +def spectral_diag_nodal_mass_matrix( + quadrature: TensorProductQuadrature + ) -> np.ndarray: + """Return the diagonal mass matrix for use in the spectral element method. + This mass matrix is exact for Lagrange polynomials with respect to + Gauss-Legendre (GL) nodes, using GL nodal degrees of freedom. + It is approximate for Lagrange polynomials with respect to the + Gauss-Lobatto-Legendre (GLL) nodes, using GLL nodal degrees of freedom. + + Returns the vector of diagonal entries. + + .. versionadded:: 2024.2 + """ + if not isinstance(quadrature, TensorProductQuadrature): + raise ValueError("only applicable to tensor product discretizations") + if not all(quad.dim == 1 for quad in quadrature.quadratures): + raise ValueError("constituent quadratures of TP quadrature must be 1D") + + return quadrature.weights + + def modal_mass_matrix_for_face( face: Face, face_quad: Quadrature, trial_functions: Sequence[Callable[[np.ndarray], np.ndarray]], diff --git a/modepy/test/test_matrices.py b/modepy/test/test_matrices.py index 28f2cb5..6b7fcc6 100644 --- a/modepy/test/test_matrices.py +++ b/modepy/test/test_matrices.py @@ -21,6 +21,9 @@ """ +from typing import Tuple, cast + +import numpy as np import numpy.linalg as la import pytest @@ -70,6 +73,48 @@ def test_nodal_mass_matrix_against_quad( assert err < 1e-14 +@pytest.mark.parametrize("shape", [ + mp.Hypercube(1), + mp.Hypercube(2), + mp.Hypercube(3), + ]) +def test_tensor_product_diag_mass_matrix(shape: mp.Shape) -> None: + shape = mp.Simplex(1) + + for order in range(16): + space = mp.space_for_shape(shape, order) + basis = mp.orthonormal_basis_for_space(space, shape) + + gl_quad = mp.LegendreGaussTensorProductQuadrature(order, shape.dim) + gl_ref_mass_mat = mp.mass_matrix(basis, gl_quad.nodes) + gl_diag_mass_mat = np.diag(mp.spectral_diag_nodal_mass_matrix(gl_quad)) + + gl_err = ( + la.norm(gl_ref_mass_mat - gl_diag_mass_mat, "fro") + / la.norm(gl_ref_mass_mat, "fro") + ) + + assert gl_err < 1e-14 + + if order == 0: + # no single-node Lobatto quadratures + continue + + gll_quad = mp.LegendreGaussLobattoTensorProductQuadrature(order, shape.dim) + gll_ref_mass_mat = mp.mass_matrix(basis, gll_quad.nodes) + gll_diag_mass_mat = np.diag(mp.spectral_diag_nodal_mass_matrix(gll_quad)) + + # Note that gll_diag_mass_mat is not a good approximation of gll_ref_mass_mat + # in the matrix norm sense! + + for mid, func in zip(basis.mode_ids, basis.functions): + if max(cast(Tuple[int, ...], mid)) < order - 1: + err = np.abs( + gll_ref_mass_mat @ func(gll_quad.nodes) + - gll_diag_mass_mat @ func(gll_quad.nodes)) + assert np.max(err) < 1e-14 + + # You can test individual routines by typing # $ python test_modes.py 'test_routine()' From f90962752f2b6a265896b30191267a430cfc8079 Mon Sep 17 00:00:00 2001 From: Andreas Kloeckner Date: Fri, 7 Jun 2024 13:01:44 -0500 Subject: [PATCH 11/12] Revamp README to better advertise capabilities --- README.rst | 35 ++++++++++++++++++++++++++++++++--- 1 file changed, 32 insertions(+), 3 deletions(-) diff --git a/README.rst b/README.rst index 2e6a72d..c3b4713 100644 --- a/README.rst +++ b/README.rst @@ -1,5 +1,5 @@ -modepy: Basis Functions and Node Sets for Interpolation -======================================================= +modepy: Basis Functions, Node Sets, Quadratures +=============================================== .. image:: https://gitlab.tiker.net/inducer/modepy/badges/main/pipeline.svg :alt: Gitlab Build Status @@ -18,12 +18,41 @@ modepy: Basis Functions and Node Sets for Interpolation simplices (i.e. segments, triangles and tetrahedra) and tensor products of simplices (i.e. squares, cubes, prisms, etc.). These are a key building block for high-order unstructured discretizations, as often used in a finite -element context. It closely follows the approach taken in the book +element context. Features include: + +- Support for simplex and tensor product elements in any dimension. +- Orthogonal bases: + - Jacobi polynomials with derivatives + - Orthogonal polynomials for simplices up to 3D and tensor product elements + and their derivatives. + - All bases permit symbolic evaluation, for code generation. +- Access to numerous quadrature rules: + - Jacobi-Gauss, Jacobi-Gauss-Lobatto in 1D + (includes Legendre, Chebyshev, ultraspherical, Gegenbauer) + - Clenshaw-Curtis and Fejér in 1D + - Grundmann-Möller on the simplex + - Xiao-Gimbutas on the simplex + - Vioreanu-Rokhlin on the simplex + - Jaśkowiec-Sukumar on the tetrahedron + - Witherden-Vincent on the hypercube + - Generic tensor products built on the above, e.g. for prisms and hypercubes +- Matrices for FEM, usable across all element types: + - generalized Vandermonde, + - mass matrices (including lumped diagonal), + - face mass matrices, + - differentiation matrices, and + - resampling matrices. +- Objects to represent 'element shape' and 'function space', + generic node/mode/quadrature retrieval based on them. + +Its roots closely followed the approach taken in the book Hesthaven, Jan S., and Tim Warburton. "Nodal Discontinuous Galerkin Methods: Algorithms, Analysis, and Applications". 1st ed. Springer, 2007. `Book web page `_ +but much has been added beyond that basic functionality. + Resources: * `documentation `_ From 40e3d9b284514f21ebd5f38828507c13d8e3b641 Mon Sep 17 00:00:00 2001 From: Andreas Kloeckner Date: Fri, 7 Jun 2024 14:27:37 -0500 Subject: [PATCH 12/12] Fix implicit use of deprecated force_dim_axis=False --- modepy/nodes.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/modepy/nodes.py b/modepy/nodes.py index 04b49a7..e933a83 100644 --- a/modepy/nodes.py +++ b/modepy/nodes.py @@ -401,7 +401,7 @@ def legendre_gauss_tensor_product_nodes(dims: int, n: int) -> np.ndarray: .. versionadded:: 2024.2 """ from modepy.quadrature.jacobi_gauss import LegendreGaussQuadrature - gl_quad = LegendreGaussQuadrature(n) + gl_quad = LegendreGaussQuadrature(n, force_dim_axis=True) gl_nodes = gl_quad.nodes.reshape(1, -1) return tensor_product_nodes(dims, gl_nodes)