diff --git a/examples/planewave-quad-comparison.py b/examples/planewave-quad-comparison.py new file mode 100644 index 000000000..65ac19b00 --- /dev/null +++ b/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 / + +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)