From f2061d70e78858f303206401e3c77d941d0a4b4a Mon Sep 17 00:00:00 2001 From: Shawn Lin Date: Tue, 14 Jan 2025 11:06:47 -0600 Subject: [PATCH] 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)