quadrature comparison

examples/planewave-quad-comparison.py

@@ -0,0 +1,197 @@
+"""
+Tensor Product Quadrature vs. Vioreanu-Rokhlin Quadrature for Plane Wave on Sphere
+==================================================================================
+
+This test compares the absolute error of **Tensor Product Quadrature** and
+**Vioreanu-Rokhlin Quadrature** against the number of discretization nodes
+with matched total polynomial degree exactness.
+
+Comparison of Polynomial Exactness
+----------------------------------
+
+The following table summarizes the total degree of polynomial exactness of both
+quadrature methods based on the order:
+
+Order     VR Exact_to     Tensor Exact_to
+-----------------------------------------
+1         2              3
+2         4              5
+3         5              7
+4         7              9
+5         8              11
+6         10             13
+7         12             15
+8         14             17
+9         15             19
+10        17             21
+11        19             23
+12        20             25
+13        22             27
+14        24             29
+15        25             31
+16        27             33
+17        28             35
+18        30             37
+19        32             39
+
+
+Wave Function and Sphere Integral
+---------------------------------
+
+The normal direction is:
+
+    d = [-5, 4, 1], n = d / <d, d>
+
+The plane wave is defined as:
+
+    f(x) = exp(1j * n · x)
+
+We compute the integral of the plane wave over a sphere of radius 1:
+
+    ∫_sphere f dS ≈ 10.57423625632583807548
+
+This value is obtained using Mathematica with a working precision of 21 digits.
+
+Mathematica Code
+----------------
+
+Below is the Mathematica code used to define the wave function and
+compute the integral numerically:
+
+    n = Normalize[{-5, 4, 1}];
+    r = 1;
+    wave[θ_, φ_] :=
+      Exp[I r (n[[1]] Sin[θ] Cos[φ] +
+               n[[2]] Sin[θ] Sin[φ] +
+               n[[3]] Cos[θ])];
+
+    NIntegrate[
+     wave[θ, φ] * r^2 Sin[θ],
+     {θ, 0, Pi},
+     {φ, 0, 2 Pi},
+     WorkingPrecision -> 21
+    ]
+
+"""
+
+import meshmode.mesh.generation as mgen
+import numpy as np
+import matplotlib.pyplot as plt
+from meshmode.discretization import Discretization
+from meshmode.discretization.poly_element import InterpolatoryQuadratureGroupFactory
+from meshmode.mesh import (
+    SimplexElementGroup,
+    TensorProductElementGroup,
+)
+from pytential.qbx import QBXLayerPotentialSource
+from pytential import GeometryCollection, bind, sym
+from arraycontext import flatten
+from meshmode.array_context import PyOpenCLArrayContext
+import pyopencl as cl
+
+cl_ctx = cl.create_some_context()
+queue = cl.CommandQueue(cl_ctx)
+actx = PyOpenCLArrayContext(queue)
+
+
+def quadrature(level, target_order, qbx_order, group_cls=SimplexElementGroup):
+    mesh = mgen.generate_sphere(1, target_order,
+            uniform_refinement_rounds=level, group_cls=group_cls)
+    pre_density_discr = Discretization(actx, mesh,
+                    InterpolatoryQuadratureGroupFactory(target_order))
+
+    qbx = QBXLayerPotentialSource(
+        pre_density_discr, target_order, qbx_order,
+        fmm_order=False)
+    discr_stage = sym.QBX_SOURCE_STAGE1
+    places = GeometryCollection({"qbx": qbx}, auto_where="qbx")
+    density_discr = places.get_discretization("qbx", discr_stage)
+    ambient_dim = qbx.ambient_dim
+    dofdesc = sym.DOFDescriptor("qbx", discr_stage)
+
+    sources = density_discr.nodes()
+    weights_nodes = bind(places,
+            sym.weights_and_area_elements(ambient_dim=3, dim=2, dofdesc=dofdesc))(actx)
+
+    sources_host = actx.to_numpy(flatten(sources, actx)).reshape(ambient_dim, -1)
+    weights_nodes_host = actx.to_numpy(flatten(weights_nodes, actx))
+
+    return sources_host, weights_nodes_host
+
+
+def wave(x):
+    n = np.array([-5, 4, 1])
+    n = n / np.linalg.norm(n)
+    return np.exp(1j * np.dot(n, x))
+
+
+def run_test(vr_target_orders, tensor_target_orders, refine_levels):
+    ref = 10.57423625632583807548
+
+    for vr_target_order, tensor_target_order in zip(
+        vr_target_orders, tensor_target_orders, strict=False):
+        print(
+            f"{'VR Order'}: {vr_target_order} "
+            f"Tensor Order: {tensor_target_order}")
+        print(
+            f"{'VR Nodes':<15}{'Tensor Nodes':<15}"
+            f"{'VR Error':<20}{'Tensor Error':<20}")
+        print("-" * 70)
+        vr_result = []
+        tensor_result = []
+        vr_nodes = []
+        tensor_nodes = []
+        vr_err = []
+        tensor_err = []
+
+        for level in refine_levels:
+            # VR quadrature
+            qbx_order = vr_target_order
+            sources_h, weights_nodes_h = quadrature(level,
+                            vr_target_order, qbx_order=qbx_order)
+            vr_value = np.dot(wave(sources_h), weights_nodes_h)
+            vr_result.append(vr_value)
+            vr_nodes.append(len(sources_h[0]))
+            vr_err.append(np.abs(vr_value - ref))
+
+            # Tensor quadrature
+            qbx_order = tensor_target_order
+            sources_h, weights_nodes_h = quadrature(level,
+                            tensor_target_order, qbx_order=qbx_order,
+                            group_cls=TensorProductElementGroup)
+            tensor_value = np.dot(wave(sources_h), weights_nodes_h)
+            tensor_result.append(tensor_value)
+            tensor_nodes.append(len(sources_h[0]))
+            tensor_err.append(np.abs(tensor_value - ref))
+
+            print(f"{vr_nodes[-1]:<15}{tensor_nodes[-1]:<15}"
+                  f"{vr_err[-1]:<20.12e}{tensor_err[-1]:<20.12e}")
+
+            if tensor_err[-1] <= 1e-13 or vr_err[-1] <= 1e-13:
+                break
+
+        print("\n")
+
+        plt.figure()
+        plt.semilogy(vr_nodes, vr_err, "o-",
+                label=f"Vioreanu-Rokhlin (Order {vr_target_order})")
+        plt.semilogy(tensor_nodes, tensor_err, "o-",
+                label=f"Tensor (Order {tensor_target_order})")
+        plt.xlabel(r"$\# \mathrm{nodes}$")
+        plt.ylabel(r"$\log_{10}(|\mathrm{abs\ err}|)$")
+        plt.legend()
+        plt.grid(True)
+        plt.title(
+            r"$\log_{10}(|\mathrm{{abs\ err}}|) \ \mathrm{{vs}} \ \# \mathrm{{nodes}}$"
+            "\n"
+            rf"$\mathrm{{VR\ order}} = {vr_target_order},"
+            rf" \mathrm{{Tensor\ order}} = {tensor_target_order}$"
+        )
+        plt.show()
+
+
+if __name__ == "__main__":
+    refine_levels = [0, 1, 2, 3, 4, 5, 6]
+    vr_target_orders = [4, 9, 16]
+    tensor_target_orders = [3, 7, 13]
+    run_test(vr_target_orders, tensor_target_orders, refine_levels)