From eea3c671468d2f2fda68c8ccd1223961f42061a0 Mon Sep 17 00:00:00 2001 From: Shawn Lin Date: Mon, 13 Jan 2025 19:09:41 -0600 Subject: [PATCH 1/3] quadrature comparison --- examples/planewave_quad_comparison.py | 177 ++++++++++++++++++++++++++ 1 file changed, 177 insertions(+) create mode 100644 examples/planewave_quad_comparison.py diff --git a/examples/planewave_quad_comparison.py b/examples/planewave_quad_comparison.py new file mode 100644 index 000000000..140162909 --- /dev/null +++ b/examples/planewave_quad_comparison.py @@ -0,0 +1,177 @@ +""" +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 InterpolatoryQuadratureSimplexGroupFactory +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, force_device_scalars=True) + +def quadrature(level, target_order, qbx_order, tensor=False): + if tensor: + from meshmode.mesh import TensorProductElementGroup + mesh = mgen.generate_sphere(1, target_order, uniform_refinement_rounds=level, group_cls=TensorProductElementGroup) + from meshmode.discretization.poly_element import InterpolatoryQuadratureGroupFactory + pre_density_discr = Discretization(actx, mesh, InterpolatoryQuadratureGroupFactory(target_order)) + else: + mesh = mgen.generate_sphere(1, target_order, uniform_refinement_rounds=level) + pre_density_discr = Discretization(actx, mesh, InterpolatoryQuadratureSimplexGroupFactory(target_order)) + + qbx = QBXLayerPotentialSource(pre_density_discr, target_order, qbx_order, fmm_order=False) + dis_stage = sym.QBX_SOURCE_STAGE1 + places = GeometryCollection({"qbx": qbx}, auto_where=('qbx')) + density_discr = places.get_discretization("qbx", dis_stage) + ambient_dim = qbx.ambient_dim + dofdesc = sym.DOFDescriptor("qbx", dis_stage) + + sources = density_discr.nodes() + weights_nodes = bind(places, sym.weights_and_area_elements(ambient_dim=3, dim=2, dofdesc=dofdesc))(actx) + + sources_h = actx.to_numpy(flatten(sources, actx)).reshape(ambient_dim, -1) + weights_nodes_h = actx.to_numpy(flatten(weights_nodes, actx)) + + return sources_h, weights_nodes_h + +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): + print(f"{'VR Order'}: {vr_target_order}, Tensor Order: {tensor_target_order}") + print(f"{'VR Nodes':<15}{'Tensor Nodes':<15}{'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, tensor=True) + 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( + rf"$\log_{{10}}(|\mathrm{{abs\ err}}|) \ \mathrm{{vs}} \ \# \mathrm{{nodes}}$" "\n" + rf"$\mathrm{{VR\ order}} = {vr_target_order}, \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) \ No newline at end of file From e8cd815c55b8721967a45e1ab1bf8348d8e0185e Mon Sep 17 00:00:00 2001 From: Shawn Lin Date: Tue, 14 Jan 2025 10:55:30 -0600 Subject: [PATCH 2/3] fix + name change --- ...{planewave_quad_comparison.py => planewave-quad-comparison.py} | 0 1 file changed, 0 insertions(+), 0 deletions(-) rename examples/{planewave_quad_comparison.py => planewave-quad-comparison.py} (100%) diff --git a/examples/planewave_quad_comparison.py b/examples/planewave-quad-comparison.py similarity index 100% rename from examples/planewave_quad_comparison.py rename to examples/planewave-quad-comparison.py From f2061d70e78858f303206401e3c77d941d0a4b4a Mon Sep 17 00:00:00 2001 From: Shawn Lin Date: Tue, 14 Jan 2025 11:06:47 -0600 Subject: [PATCH 3/3] fix --- examples/planewave-quad-comparison.py | 126 +++++++++++++++----------- 1 file changed, 73 insertions(+), 53 deletions(-) diff --git a/examples/planewave-quad-comparison.py b/examples/planewave-quad-comparison.py index 140162909..65ac19b00 100644 --- a/examples/planewave-quad-comparison.py +++ b/examples/planewave-quad-comparison.py @@ -2,15 +2,15 @@ Tensor Product Quadrature vs. Vioreanu-Rokhlin Quadrature for Plane Wave on Sphere ================================================================================== -This test compares the absolute error of **Tensor Product Quadrature** and +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: +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 ----------------------------------------- @@ -55,19 +55,20 @@ Mathematica Code ---------------- -Below is the Mathematica code used to define the wave function and compute the integral numerically: +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 = 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}, + wave[θ, φ] * r^2 Sin[θ], + {θ, 0, Pi}, + {φ, 0, 2 Pi}, WorkingPrecision -> 21 ] @@ -77,7 +78,11 @@ import numpy as np import matplotlib.pyplot as plt from meshmode.discretization import Discretization -from meshmode.discretization.poly_element import InterpolatoryQuadratureSimplexGroupFactory +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 @@ -86,44 +91,51 @@ cl_ctx = cl.create_some_context() queue = cl.CommandQueue(cl_ctx) -actx = PyOpenCLArrayContext(queue, force_device_scalars=True) - -def quadrature(level, target_order, qbx_order, tensor=False): - if tensor: - from meshmode.mesh import TensorProductElementGroup - mesh = mgen.generate_sphere(1, target_order, uniform_refinement_rounds=level, group_cls=TensorProductElementGroup) - from meshmode.discretization.poly_element import InterpolatoryQuadratureGroupFactory - pre_density_discr = Discretization(actx, mesh, InterpolatoryQuadratureGroupFactory(target_order)) - else: - mesh = mgen.generate_sphere(1, target_order, uniform_refinement_rounds=level) - pre_density_discr = Discretization(actx, mesh, InterpolatoryQuadratureSimplexGroupFactory(target_order)) - - qbx = QBXLayerPotentialSource(pre_density_discr, target_order, qbx_order, fmm_order=False) - dis_stage = sym.QBX_SOURCE_STAGE1 - places = GeometryCollection({"qbx": qbx}, auto_where=('qbx')) - density_discr = places.get_discretization("qbx", dis_stage) +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", dis_stage) + 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) + weights_nodes = bind(places, + sym.weights_and_area_elements(ambient_dim=3, dim=2, dofdesc=dofdesc))(actx) - sources_h = actx.to_numpy(flatten(sources, actx)).reshape(ambient_dim, -1) - weights_nodes_h = actx.to_numpy(flatten(weights_nodes, 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 - return sources_h, weights_nodes_h 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): - print(f"{'VR Order'}: {vr_target_order}, Tensor Order: {tensor_target_order}") - print(f"{'VR Nodes':<15}{'Tensor Nodes':<15}{'VR Error':<20}{'Tensor Error':<20}") + 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 = [] @@ -131,11 +143,12 @@ def run_test(vr_target_orders, tensor_target_orders, refine_levels): 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) + 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])) @@ -143,7 +156,9 @@ def run_test(vr_target_orders, tensor_target_orders, refine_levels): # Tensor quadrature qbx_order = tensor_target_order - sources_h, weights_nodes_h = quadrature(level, tensor_target_order, qbx_order=qbx_order, tensor=True) + 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])) @@ -151,27 +166,32 @@ def run_test(vr_target_orders, tensor_target_orders, refine_levels): 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.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( - rf"$\log_{{10}}(|\mathrm{{abs\ err}}|) \ \mathrm{{vs}} \ \# \mathrm{{nodes}}$" "\n" - rf"$\mathrm{{VR\ order}} = {vr_target_order}, \mathrm{{Tensor\ order}} = {tensor_target_order}$" + 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) \ No newline at end of file + vr_target_orders = [4, 9, 16] + tensor_target_orders = [3, 7, 13] + run_test(vr_target_orders, tensor_target_orders, refine_levels)