Skip to content

Commit

Permalink
fix
Browse files Browse the repository at this point in the history
  • Loading branch information
@ShawnL00
ShawnL00 committed Jan 14, 2025
1 parent e8cd815 commit f2061d7
Showing 1 changed file with 73 additions and 53 deletions.
126 changes: 73 additions & 53 deletions examples/planewave-quad-comparison.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
-----------------------------------------
Expand Down Expand Up @@ -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
]
Expand All @@ -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
Expand All @@ -86,92 +91,107 @@

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 = []
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)
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)
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.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)
vr_target_orders = [4, 9, 16]
tensor_target_orders = [3, 7, 13]
run_test(vr_target_orders, tensor_target_orders, refine_levels)

0 comments on commit f2061d7

Please sign in to comment.