Skip to content

Parametric recurrence ported to actx #219

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Draft
wants to merge 1 commit into
base: parametric_recurrence
Choose a base branch
from
Draft

Parametric recurrence ported to actx #219

Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
33 changes: 21 additions & 12 deletions test/test_recurrenceqbx.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -6,8 +6,13 @@

import numpy as np
import sympy as sp

from sympy import hankel1

from arraycontext import (
ArrayContext,
pytest_generate_tests_for_array_contexts,
)

from sumpy.array_context import _acf
from sumpy.expansion.diff_op import (
laplacian,
Expand All @@ -19,10 +24,11 @@
from sumpy.recurrenceqbx import _make_sympy_vec, recurrence_qbx_lp


actx_factory = _acf
pytest_generate_tests = pytest_generate_tests_for_array_contexts([
"pytato:pyopencl",
])
ExpnClass = LineTaylorLocalExpansion

actx = actx_factory()
lknl = LaplaceKernel(2)
hlknl = HelmholtzKernel(2, "k")

Expand Down Expand Up @@ -74,29 +80,32 @@ def _qbx_lp_laplace_general(sources, targets, centers, radius, strengths, order)
return result_qbx


def _create_ellipse(n_p):
def _create_ellipse(actx: ArrayContext, n_p):
h = 9.688 / n_p
radius = 7*h
t = np.linspace(0, 2 * np.pi, n_p, endpoint=False)
t = actx.np.linspace(0, 2 * np.pi, n_p, endpoint=False)

unit_circle_param = np.exp(1j * t)
unit_circle = np.array([2 * unit_circle_param.real, unit_circle_param.imag])
unit_circle_param = actx.np.exp(1j * t)
unit_circle = actx.np.stack([2 * unit_circle_param.real, unit_circle_param.imag],
axis=0)

sources = unit_circle
normals = np.array([unit_circle_param.real, 2*unit_circle_param.imag])
normals = normals / np.linalg.norm(normals, axis=0)
normals = actx.np.stack([unit_circle_param.real, 2*unit_circle_param.imag], axis=0)
# normals = normals / actx.np.linalg.norm(normals, axis=0)
normals = normals / np.sum(actx.np.sum(normals**2, axis=0))
centers = sources - normals * radius

mode_nr = 25
density = np.cos(mode_nr * t)
density = actx.np.cos(mode_nr * t)

return sources, centers, normals, density, h, radius


def test_recurrence_laplace_2d_ellipse():
def test_recurrence_laplace_2d_ellipse(actx_factory):
r"""
Tests recurrence + qbx code.
"""
actx = actx_factory()

# ------------- 1. Define PDE, Green's Function
w = make_identity_diff_op(2)
Expand All @@ -110,7 +119,7 @@ def test_recurrence_laplace_2d_ellipse():
p = 4
err = []
for n_p in range(200, 1001, 200):
sources, centers, normals, density, h, radius = _create_ellipse(n_p)
sources, centers, normals, density, h, radius = _create_ellipse(actx, n_p)
strengths = h * density
exp_res = recurrence_qbx_lp(sources, centers, normals,
strengths, radius, laplace2d,
Expand Down
Loading