From e4520983dc32907c23490321b8c57b8aa2c14fa7 Mon Sep 17 00:00:00 2001 From: Isuru Fernando Date: Mon, 20 Apr 2026 09:39:05 -0500 Subject: [PATCH 1/8] Add heat kernel Squashed from https://github.com/inducer/sumpy/pull/185 --- .github/workflows/ci.yml | 2 + .gitlab-ci.yml | 4 ++ examples/heat-local.py | 85 +++++++++++++++++++++++++++++++++++++ examples/heat-mpole.py | 87 ++++++++++++++++++++++++++++++++++++++ pyproject.toml | 1 + sumpy/e2e.py | 7 ++- sumpy/expansion/m2l.py | 4 -- sumpy/kernel.py | 75 ++++++++++++++++++++++++++++++++ sumpy/test/test_kernels.py | 73 ++++++++++++++++++++++++-------- sumpy/test/test_misc.py | 36 ++++++++++++---- sumpy/toys.py | 4 +- 11 files changed, 347 insertions(+), 31 deletions(-) create mode 100644 examples/heat-local.py create mode 100644 examples/heat-mpole.py diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index a22d16eb0..094ad8299 100644 --- a/.github/workflows/ci.yml +++ b/.github/workflows/ci.yml @@ -68,6 +68,7 @@ jobs: run: | grep -v symengine .test-conda-env-py3.yml > .test-conda-env.yml CONDA_ENVIRONMENT=.test-conda-env.yml + export PYTEST_ADDOPTS=${PYTEST_ADDOPTS:-"-m 'not slowtest'"} curl -L -O https://gitlab.tiker.net/inducer/ci-support/raw/main/build-and-test-py-project-within-miniconda.sh . ./build-and-test-py-project-within-miniconda.sh @@ -89,6 +90,7 @@ jobs: - uses: actions/checkout@v6 - name: "Main Script" run: | + export PYTEST_ADDOPTS=${PYTEST_ADDOPTS:-"-m 'not slowtest'"} curl -L -O https://gitlab.tiker.net/inducer/ci-support/raw/main/build-and-test-py-project-within-miniconda.sh . ./build-and-test-py-project-within-miniconda.sh diff --git a/.gitlab-ci.yml b/.gitlab-ci.yml index 13bca8b85..51e5a8f93 100644 --- a/.gitlab-ci.yml +++ b/.gitlab-ci.yml @@ -22,6 +22,7 @@ Pytest POCL: - export PY_EXE=python3 - export PYOPENCL_TEST=portable:pthread - export EXTRA_INSTALL="pybind11 numpy mako mpi4py" + - export PYTEST_ADDOPTS=${PYTEST_ADDOPTS:-"-m 'not slowtest'"} - curl -L -O https://gitlab.tiker.net/inducer/ci-support/raw/main/build-and-test-py-project.sh - ". ./build-and-test-py-project.sh" tags: @@ -38,6 +39,7 @@ Pytest Titan V: - py_version=3 - export PYOPENCL_TEST=nvi:titan - EXTRA_INSTALL="pybind11 numpy mako mpi4py" + - export PYTEST_ADDOPTS=${PYTEST_ADDOPTS:-"-m 'not slowtest'"} - curl -L -O https://gitlab.tiker.net/inducer/ci-support/raw/main/build-and-test-py-project.sh - ". ./build-and-test-py-project.sh" tags: @@ -56,6 +58,7 @@ Pytest Conda: - export SUMPY_NO_CACHE=1 - export SUMPY_FORCE_SYMBOLIC_BACKEND=symengine - export PYOPENCL_TEST=portable:pthread + - export PYTEST_ADDOPTS=${PYTEST_ADDOPTS:-"-m 'not slowtest'"} - curl -L -O https://gitlab.tiker.net/inducer/ci-support/raw/main/build-and-test-py-project-within-miniconda.sh - ". ./build-and-test-py-project-within-miniconda.sh" tags: @@ -72,6 +75,7 @@ Pytest POCL Titan V: - export SUMPY_NO_CACHE=1 - export SUMPY_FORCE_SYMBOLIC_BACKEND=symengine - export PYOPENCL_TEST=portable:titan + - export PYTEST_ADDOPTS=${PYTEST_ADDOPTS:-"-m 'not slowtest'"} - curl -L -O https://gitlab.tiker.net/inducer/ci-support/raw/main/build-and-test-py-project-within-miniconda.sh - ". ./build-and-test-py-project-within-miniconda.sh" tags: diff --git a/examples/heat-local.py b/examples/heat-local.py new file mode 100644 index 000000000..fca53353e --- /dev/null +++ b/examples/heat-local.py @@ -0,0 +1,85 @@ +import numpy as np + +import pyopencl as cl + +from sumpy import toys +from sumpy.visualization import FieldPlotter +from sumpy.kernel import ( # noqa: F401 + YukawaKernel, + HeatKernel, + HelmholtzKernel, + LaplaceKernel) + +try: + import matplotlib.pyplot as plt + USE_MATPLOTLIB = True +except ImportError: + USE_MATPLOTLIB = False + + +def main(): + from sumpy.array_context import PyOpenCLArrayContext + ctx = cl.create_some_context() + queue = cl.CommandQueue(ctx) + actx = PyOpenCLArrayContext(queue, force_device_scalars=True) + + tctx = toys.ToyContext( + actx.context, + # LaplaceKernel(2), + #YukawaKernel(2), extra_kernel_kwargs={"lam": 5}, + # HelmholtzKernel(2), extra_kernel_kwargs={"k": 0.3}, + HeatKernel(1), extra_kernel_kwargs={"alpha": 1}, + ) + + src_size = 0.1 + pt_src = toys.PointSources( + tctx, + np.array([ + src_size*(np.random.rand(50) - 0.5), + np.zeros(50)]), + np.random.randn(50)) + + fp = FieldPlotter([0, 0.5], extent=np.array([8, 1])) + + if 0 and USE_MATPLOTLIB: + toys.logplot(fp, pt_src, cmap="jet", aspect=8) + plt.colorbar() + plt.show() + + + p = 5 + center = np.array([0, 1], dtype=np.float64) + mexp = toys.local_expand(pt_src, center, p) + diff = mexp - pt_src + + dist = fp.points - center[:, None] + r = np.sqrt(dist[0]**2 + dist[1]**2) + + error_model = r**(p+1) + + def logplot_fp(fp: FieldPlotter, values, **kwargs) -> None: + fp.show_scalar_in_matplotlib( + np.log10(np.abs(values + 1e-15)), **kwargs) + if USE_MATPLOTLIB: + plt.subplot(131) + logplot_fp(fp, error_model, cmap="jet", vmin=-5, vmax=0, aspect=8) + plt.subplot(132) + logplot_fp(fp, diff.eval(fp.points), cmap="jet", vmin=-5, vmax=0, aspect=8) + plt.subplot(133) + logplot_fp(fp, diff.eval(fp.points)/error_model, cmap="jet", vmin=-5, vmax=0, aspect=8) + plt.colorbar() + plt.show() + 1/0 + mexp2 = toys.multipole_expand(mexp, [0, 0.25]) # noqa: F841 + lexp = toys.local_expand(mexp, [3, 0]) + lexp2 = toys.local_expand(lexp, [3, 1], 3) + + # diff = mexp - pt_src + # diff = mexp2 - pt_src + diff = lexp2 - pt_src + + print(toys.l_inf(diff, 1.2, center=lexp2.center)) + + +if __name__ == "__main__": + main() diff --git a/examples/heat-mpole.py b/examples/heat-mpole.py new file mode 100644 index 000000000..7bc4cc6ac --- /dev/null +++ b/examples/heat-mpole.py @@ -0,0 +1,87 @@ +import numpy as np + +import pyopencl as cl + +from sumpy import toys +from sumpy.visualization import FieldPlotter +from sumpy.kernel import ( # noqa: F401 + YukawaKernel, + HeatKernel, + HelmholtzKernel, + LaplaceKernel) + +try: + import matplotlib.pyplot as plt + USE_MATPLOTLIB = True +except ImportError: + USE_MATPLOTLIB = False + + +def main(): + from sumpy.array_context import PyOpenCLArrayContext + ctx = cl.create_some_context() + queue = cl.CommandQueue(ctx) + actx = PyOpenCLArrayContext(queue, force_device_scalars=True) + + tctx = toys.ToyContext( + actx.context, + # LaplaceKernel(2), + #YukawaKernel(2), extra_kernel_kwargs={"lam": 5}, + # HelmholtzKernel(2), extra_kernel_kwargs={"k": 0.3}, + HeatKernel(1), extra_kernel_kwargs={"alpha": 1}, + ) + + src_size = 0.1 + pt_src = toys.PointSources( + tctx, + np.array([ + src_size*(np.random.rand(50) - 0.5), + np.zeros(50)]), + np.random.randn(50)) + + fp = FieldPlotter([0, 0.5], extent=np.array([8, 1])) + + if 0 and USE_MATPLOTLIB: + toys.logplot(fp, pt_src, cmap="jet", aspect=8) + plt.colorbar() + plt.show() + + + p = 5 + mexp = toys.multipole_expand(pt_src, [0, 0], p) + diff = mexp - pt_src + + x, t = fp.points + + delta = 4 * t + conv_factor = src_size / np.sqrt(2*delta) + + error_model = conv_factor**(p+1)*np.exp(-(x/np.sqrt(delta))**2/2) + #error_model = conv_factor**(p+1)/(1-conv_factor)*np.exp(-(x/np.sqrt(delta))**2) + + def logplot_fp(fp: FieldPlotter, values, **kwargs) -> None: + fp.show_scalar_in_matplotlib( + np.log10(np.abs(values + 1e-15)), **kwargs) + if USE_MATPLOTLIB: + plt.subplot(131) + logplot_fp(fp, error_model, cmap="jet", vmin=-5, vmax=0, aspect=8) + plt.subplot(132) + logplot_fp(fp, diff.eval(fp.points), cmap="jet", vmin=-5, vmax=0, aspect=8) + plt.subplot(133) + logplot_fp(fp, diff.eval(fp.points)/error_model, cmap="jet", vmin=-5, vmax=0, aspect=8) + plt.colorbar() + plt.show() + 1/0 + mexp2 = toys.multipole_expand(mexp, [0, 0.25]) # noqa: F841 + lexp = toys.local_expand(mexp, [3, 0]) + lexp2 = toys.local_expand(lexp, [3, 1], 3) + + # diff = mexp - pt_src + # diff = mexp2 - pt_src + diff = lexp2 - pt_src + + print(toys.l_inf(diff, 1.2, center=lexp2.center)) + + +if __name__ == "__main__": + main() diff --git a/pyproject.toml b/pyproject.toml index 3233fc0e5..7042587ff 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -154,6 +154,7 @@ extend-exclude = [ [tool.pytest.ini_options] markers = [ "mpi: tests distributed FMM", + "slowtest: mark a test as slow", ] [tool.basedpyright] diff --git a/sumpy/e2e.py b/sumpy/e2e.py index 0e5bdfd76..d4eb0b023 100644 --- a/sumpy/e2e.py +++ b/sumpy/e2e.py @@ -320,7 +320,12 @@ def get_translation_loopy_insns(self, result_dtype): [sym.Symbol(f"data{i}") for i in range(m2l_translation_classes_dependent_ndata)] - ncoeff_src = len(self.src_expansion) + m2l_translation = self.tgt_expansion.m2l_translation + if m2l_translation.use_preprocessing: + ncoeff_src = m2l_translation.preprocess_multipole_nexprs( + self.tgt_expansion, self.src_expansion) + else: + ncoeff_src = len(self.src_expansion) src_coeff_exprs = [sym.Symbol(f"src_coeffs{i}") for i in range(ncoeff_src)] diff --git a/sumpy/expansion/m2l.py b/sumpy/expansion/m2l.py index 5a77df51f..b3782c73c 100644 --- a/sumpy/expansion/m2l.py +++ b/sumpy/expansion/m2l.py @@ -879,10 +879,6 @@ def loopy_translate(self, vector = translation_classes_dependent_data[icoeff_src] expr = toeplitz_first_row * vector - domains.append(f"{{[icoeff_src]: icoeff_tgt<=icoeff_src<{ncoeff_src} }}") - - expr = src_coeffs[icoeff_tgt] * translation_classes_dependent_data[icoeff_tgt] - insns = [ lp.Assignment( assignee=tgt_coeffs[icoeff_tgt], diff --git a/sumpy/kernel.py b/sumpy/kernel.py index df9c04686..a3d1af5ad 100644 --- a/sumpy/kernel.py +++ b/sumpy/kernel.py @@ -106,6 +106,9 @@ .. autoclass:: LineOfCompressionKernel :show-inheritance: :members: mapper_method +.. autoclass:: HeatKernel + :show-inheritance: + :members: mapper_method Derivatives ----------- @@ -1119,6 +1122,76 @@ def get_pde_as_diff_op(self) -> LinearPDESystemOperator: return laplacian(w) +class HeatKernel(ExpressionKernel): + r"""The Green's function for the heat equation given by + :math:`e^{-r^2/{4 \alpha t}}/\sqrt{(4 \pi \alpha t)^d}` + where :math:`d` is the number of spatial dimensions. + + .. note:: + + This kernel cannot be used in an FMM yet and can only + be used in expansions and evaluations that occur forward + in the time dimension. + """ + heat_alpha_name: str + + mapper_method: ClassVar[str] = "map_heat_kernel" + + def __init__(self, spatial_dims: int, heat_alpha_name: str = "alpha"): + dim = spatial_dims + 1 + d = make_sym_vector("d", dim) + t = d[-1] + r = pymbolic_real_norm_2(d[:-1]) + alpha = SpatialConstant(heat_alpha_name) + expr = var("exp")(-r**2/(4 * alpha * t)) / var("sqrt")(t**(dim - 1)) + scaling = 1/var("sqrt")((4*var("pi")*alpha)**(dim - 1)) + + super().__init__( + dim, + expression=expr, + global_scaling_const=scaling, + ) + + self.heat_alpha_name = heat_alpha_name + + @property + @override + def is_complex_valued(self) -> bool: + return False + + def update_persistent_hash(self, key_hash, key_builder): + key_hash.update(type(self).__name__.encode("utf8")) + key_builder.rec(key_hash, (self.dim, self.heat_alpha_name)) + + @override + def __repr__(self): + return f"HeatKnl{self.dim - 1}D" + + @override + def get_args(self): + return [ + KernelArgument( + loopy_arg=lp.ValueArg(self.heat_alpha_name, np.float64), + )] + + def get_derivative_taker(self, dvec, rscale, sac): + """Return a :class:`sumpy.derivative_taker.ExprDerivativeTaker` instance + that supports taking derivatives of the base kernel with respect to dvec. + """ + from sumpy.derivative_taker import ExprDerivativeTaker + return ExprDerivativeTaker(self.get_expression(dvec), dvec, rscale, + sac) + + @override + def get_pde_as_diff_op(self): + from sumpy.expansion.diff_op import diff, laplacian, make_identity_diff_op + alpha = sym.Symbol(self.heat_alpha_name) + w = make_identity_diff_op(self.dim - 1, time_dependent=True) + time_diff = [0]*self.dim + time_diff[-1] = 1 + return diff(w, tuple(time_diff)) - laplacian(w) * alpha + + # }}} @@ -1590,6 +1663,7 @@ def map_expression_kernel(self, kernel: ExpressionKernel) -> Kernel: map_elasticity_kernel = map_expression_kernel map_line_of_compression_kernel = map_expression_kernel map_stresslet_kernel = map_expression_kernel + map_heat_kernel = map_expression_kernel def map_axis_target_derivative(self, kernel: AxisTargetDerivative) -> Kernel: return type(kernel)(kernel.axis, self.rec(kernel.inner_kernel)) @@ -1662,6 +1736,7 @@ def map_expression_kernel(self, kernel: ExpressionKernel) -> int: map_yukawa_kernel = map_expression_kernel map_line_of_compression_kernel = map_expression_kernel map_stresslet_kernel = map_expression_kernel + map_heat_kernel = map_expression_kernel @override def map_axis_target_derivative(self, kernel: AxisTargetDerivative) -> int: diff --git a/sumpy/test/test_kernels.py b/sumpy/test/test_kernels.py index 24fa71087..996644829 100644 --- a/sumpy/test/test_kernels.py +++ b/sumpy/test/test_kernels.py @@ -64,6 +64,7 @@ BiharmonicKernel, DirectionalSourceDerivative, ElasticityKernel, + HeatKernel, HelmholtzKernel, Kernel, LaplaceKernel, @@ -266,6 +267,15 @@ def test_p2e_multiple( (LaplaceKernel(2), LinearPDEConformingVolumeTaylorLocalExpansion), (LaplaceKernel(2), LinearPDEConformingVolumeTaylorMultipoleExpansion), + (HeatKernel(1), LinearPDEConformingVolumeTaylorLocalExpansion), + (HeatKernel(1), LinearPDEConformingVolumeTaylorMultipoleExpansion), + (HeatKernel(2), LinearPDEConformingVolumeTaylorLocalExpansion), + (HeatKernel(2), LinearPDEConformingVolumeTaylorMultipoleExpansion), + pytest.param(HeatKernel(3), LinearPDEConformingVolumeTaylorLocalExpansion, + marks=pytest.mark.slowtest), + pytest.param(HeatKernel(3), LinearPDEConformingVolumeTaylorMultipoleExpansion, + marks=pytest.mark.slowtest), + (HelmholtzKernel(2), VolumeTaylorMultipoleExpansion), (HelmholtzKernel(2), VolumeTaylorLocalExpansion), (HelmholtzKernel(2), LinearPDEConformingVolumeTaylorLocalExpansion), @@ -310,6 +320,8 @@ def test_p2e2p( extra_kwargs["k"] = 0.2 if isinstance(base_knl, StokesletKernel): extra_kwargs["mu"] = 0.2 + if isinstance(base_knl, HeatKernel): + extra_kwargs["alpha"] = 1.0 if with_source_derivative: knl = DirectionalSourceDerivative(base_knl, "dir_vec") @@ -338,12 +350,23 @@ def test_p2e2p( h_values = [1/2, 1/3, 1/5] rng = np.random.default_rng(19) - center = np.array([2, 1, 0][:knl.dim], np.float64) - sources = actx.from_numpy( - 0.7 * (-0.5 + rng.random((knl.dim, nsources), dtype=np.float64)) - + center[:, np.newaxis]) + + if isinstance(base_knl, HeatKernel): + # Setup sources so that there are no negative time intervals + center = np.array([0, 0, 0, 0.5][-knl.dim:], np.float64) + sources = (-0.5 + rng.random((knl.dim, nsources), dtype=np.float64) + + center[:, np.newaxis]) + loc_center = np.array([0.0, 0.0, 0.0, 6.0][-knl.dim:]) + center + varying_axis = -1 + else: + center = np.array([2, 1, 0][:knl.dim], np.float64) + sources = 0.7 * (-0.5 + rng.random((knl.dim, nsources), dtype=np.float64) + + center[:, np.newaxis]) + loc_center = np.array([5.5, 0.0, 0.0][:knl.dim]) + center + varying_axis = 0 strengths = actx.from_numpy(np.ones(nsources, dtype=np.float64) / nsources) + sources = actx.from_numpy(sources) source_boxes = actx.from_numpy(np.array([0], dtype=np.int32)) box_source_starts = actx.from_numpy(np.array([0], dtype=np.int32)) @@ -353,22 +376,21 @@ def test_p2e2p( extra_source_kwargs = extra_kwargs.copy() if isinstance(knl, DirectionalSourceDerivative): alpha = np.linspace(0, 2*np.pi, nsources, np.float64) - dir_vec = np.vstack([np.cos(alpha), np.sin(alpha)]) + dir_vec = np.vstack( + [np.cos(alpha), np.sin(alpha), alpha, alpha**2][:knl.dim]) extra_source_kwargs["dir_vec"] = actx.from_numpy(dir_vec) from sumpy.visualization import FieldPlotter for h in h_values: if issubclass(expn_class, LocalExpansionBase): - loc_center = np.array([5.5, 0.0, 0.0][:knl.dim]) + center centers = np.array(loc_center, dtype=np.float64).reshape(knl.dim, 1) fp = FieldPlotter(loc_center, extent=h, npoints=res) else: - eval_center = np.array([1/h, 0.0, 0.0][:knl.dim]) + center + eval_center = center.copy() + eval_center[varying_axis] = 1/h fp = FieldPlotter(eval_center, extent=0.1, npoints=res) - centers = (np.array([0.0, 0.0, 0.0][:knl.dim], - dtype=np.float64).reshape(knl.dim, 1) - + center[:, np.newaxis]) + centers = center[:, np.newaxis] centers = actx.from_numpy(centers) targets = actx.from_numpy(obj_array.new_1d(fp.points)) @@ -466,7 +488,7 @@ def test_p2e2p( if issubclass(expn_class, LocalExpansionBase): tgt_order_grad = tgt_order - 1 slack = 0.7 - grad_slack = 0.5 + grad_slack = 0.8 if isinstance(base_knl, HeatKernel) else 0.5 else: tgt_order_grad = tgt_order + 1 @@ -477,6 +499,10 @@ def test_p2e2p( slack += 1 grad_slack += 1 + if isinstance(base_knl, HeatKernel): + slack += 1.0 + grad_slack += 2.5 + if isinstance(knl, DirectionalSourceDerivative): slack += 1 grad_slack += 2 @@ -503,6 +529,14 @@ def test_p2e2p( # {{{ test_translations @pytest.mark.parametrize("knl, local_expn_class, mpole_expn_class, use_fft", [ + (HeatKernel(1), LinearPDEConformingVolumeTaylorLocalExpansion, + LinearPDEConformingVolumeTaylorMultipoleExpansion, False), + (HeatKernel(2), LinearPDEConformingVolumeTaylorLocalExpansion, + LinearPDEConformingVolumeTaylorMultipoleExpansion, False), + (HeatKernel(2), LinearPDEConformingVolumeTaylorLocalExpansion, + LinearPDEConformingVolumeTaylorMultipoleExpansion, True), # not in original PR + (HeatKernel(3), LinearPDEConformingVolumeTaylorLocalExpansion, + LinearPDEConformingVolumeTaylorMultipoleExpansion, False), (LaplaceKernel(2), VolumeTaylorLocalExpansion, VolumeTaylorMultipoleExpansion, False), (LaplaceKernel(2), LinearPDEConformingVolumeTaylorLocalExpansion, @@ -549,10 +583,15 @@ def test_translations( extra_kwargs["k"] = 0.05 if isinstance(knl, StokesletKernel): extra_kwargs["mu"] = 0.05 + if isinstance(knl, HeatKernel): + extra_kwargs["alpha"] = 0.1 # Just to make sure things also work away from the origin rng = np.random.default_rng(18) - origin = np.array([2, 1, 0][:knl.dim], np.float64) + if isinstance(knl, HeatKernel): + origin = np.array([0, 0, 1, 2][-knl.dim:], np.float64) + else: + origin = np.array([2, 1, 0][:knl.dim], np.float64) sources = ( 0.7 * (-0.5 + rng.random((knl.dim, nsources), dtype=np.float64)) + origin[:, np.newaxis]) @@ -565,20 +604,20 @@ def test_translations( from sumpy.visualization import FieldPlotter - eval_offset = np.array([5.5, 0.0, 0][:knl.dim]) + eval_offset = np.array([0.0, 0.0, 0.0, 5.5][-knl.dim:]) fp = FieldPlotter(eval_offset + origin, extent=0.3, npoints=res) targets = fp.points centers = (np.array( [ # box 0: particles, first mpole here - [0, 0, 0][:knl.dim], + [0, 0, 0, 0][-knl.dim:], # box 1: second mpole here - np.array([-0.2, 0.1, 0][:knl.dim], np.float64), + np.array([0, 0, 0.1, -0.2][-knl.dim:], np.float64), # box 2: first local here - eval_offset + np.array([0.3, -0.2, 0][:knl.dim], np.float64), + eval_offset + np.array([0, 0, -0.2, 0.3][-knl.dim:], np.float64), # box 3: second local and eval here eval_offset @@ -587,7 +626,7 @@ def test_translations( del eval_offset - if knl.dim == 2: # noqa: SIM108 + if knl.dim == 2 and not isinstance(knl, HeatKernel): orders = [2, 3, 4] else: orders = [3, 4, 5] diff --git a/sumpy/test/test_misc.py b/sumpy/test/test_misc.py index ce85f30b5..990e3e3b1 100644 --- a/sumpy/test/test_misc.py +++ b/sumpy/test/test_misc.py @@ -59,6 +59,7 @@ BiharmonicKernel, ElasticityKernel, ExpressionKernel, + HeatKernel, HelmholtzKernel, Kernel, LaplaceKernel, @@ -83,7 +84,9 @@ # {{{ pde check for kernels class KernelInfo: - def __init__(self, kernel, **kwargs): + kernel: Kernel + + def __init__(self, kernel: Kernel, **kwargs): self.kernel = kernel self.extra_kwargs = kwargs diff_op = self.kernel.get_pde_as_diff_op() @@ -130,8 +133,14 @@ def nderivs(self): KernelInfo(ElasticityKernel(3, 0, 0), mu=5, nu=0.2), KernelInfo(LineOfCompressionKernel(3, 0), mu=5, nu=0.2), KernelInfo(LineOfCompressionKernel(3, 1), mu=5, nu=0.2), + KernelInfo(HeatKernel(1), alpha=0.1), + KernelInfo(HeatKernel(2), alpha=0.1), + KernelInfo(HeatKernel(3), alpha=0.1), ]) -def test_pde_check_kernels(actx_factory: ArrayContextFactory, knl_info, order=5): +def test_pde_check_kernels( + actx_factory: ArrayContextFactory, + knl_info: KernelInfo, + order: int = 5): actx = actx_factory() dim = knl_info.kernel.dim @@ -150,7 +159,10 @@ def test_pde_check_kernels(actx_factory: ArrayContextFactory, knl_info, order=5) eoc_rec = EOCRecorder() for h in [0.1, 0.05, 0.025]: - cp = CalculusPatch(np.array([1, 0, 0])[:dim], h=h, order=order) + if isinstance(knl_info.kernel, HeatKernel): + cp = CalculusPatch(np.array([0, 0, 0, 1])[-dim:], h=h, order=order) + else: + cp = CalculusPatch(np.array([1, 0, 0])[:dim], h=h, order=order) pot = pt_src.eval(actx, cp.points) pde = knl_info.pde_func(cp, pot) @@ -312,7 +324,14 @@ def dim(self): @pytest.mark.parametrize("case", P2E2E2P_TEST_CASES) -def test_toy_p2e2e2p(actx_factory: ArrayContextFactory, case): +@pytest.mark.parametrize(("make_kernel", "extra_kernel_kwargs"), [ + (LaplaceKernel, {}), + (lambda dim: HeatKernel(dim - 1), {"alpha": 0.1})]) +def test_toy_p2e2e2p( + actx_factory: ArrayContextFactory, + case, + make_kernel: Callable[[int], Kernel], + extra_kernel_kwargs: dict[str, object]): dim = case.dim src = case.source.reshape(dim, -1) @@ -331,13 +350,14 @@ def test_toy_p2e2e2p(actx_factory: ArrayContextFactory, case): raise ValueError( f"convergence factor not in valid range: {case_conv_factor}") - from sumpy.expansion import VolumeTaylorExpansionFactory + from sumpy.expansion import DefaultExpansionFactory actx = actx_factory() ctx = t.ToyContext( - LaplaceKernel(dim), - expansion_factory=VolumeTaylorExpansionFactory(), - m2l_use_fft=case.m2l_use_fft) + make_kernel(dim), + expansion_factory=DefaultExpansionFactory(), + m2l_use_fft=case.m2l_use_fft, + extra_kernel_kwargs=extra_kernel_kwargs) errors = [] diff --git a/sumpy/toys.py b/sumpy/toys.py index 3ccda9c82..b366951ec 100644 --- a/sumpy/toys.py +++ b/sumpy/toys.py @@ -385,6 +385,8 @@ def _e2e(actx: ArrayContext, actx.to_numpy(to_coeffs[1]), derived_from=psource, **expn_kwargs) +# }}} + def _m2l(actx: ArrayContext, psource, to_center, to_rscale, to_order, e2e, expn_class, expn_kwargs, @@ -489,9 +491,9 @@ def _m2l(actx: ArrayContext, return ret - # }}} + # {{{ potential source classes Number_ish = int | float | complex | np.number From acbb53e0cb4c52bba4ddeaef8697fda1c26db429 Mon Sep 17 00:00:00 2001 From: Shawn Lin Date: Thu, 23 Apr 2026 11:01:37 -0500 Subject: [PATCH 2/8] heat kernel: update test_pde_check & remove test_toy_p2e2e2p --- examples/heat-local.py | 85 ---------------------------------------- examples/heat-mpole.py | 87 ----------------------------------------- sumpy/test/test_misc.py | 25 +++++------- 3 files changed, 10 insertions(+), 187 deletions(-) delete mode 100644 examples/heat-local.py delete mode 100644 examples/heat-mpole.py diff --git a/examples/heat-local.py b/examples/heat-local.py deleted file mode 100644 index fca53353e..000000000 --- a/examples/heat-local.py +++ /dev/null @@ -1,85 +0,0 @@ -import numpy as np - -import pyopencl as cl - -from sumpy import toys -from sumpy.visualization import FieldPlotter -from sumpy.kernel import ( # noqa: F401 - YukawaKernel, - HeatKernel, - HelmholtzKernel, - LaplaceKernel) - -try: - import matplotlib.pyplot as plt - USE_MATPLOTLIB = True -except ImportError: - USE_MATPLOTLIB = False - - -def main(): - from sumpy.array_context import PyOpenCLArrayContext - ctx = cl.create_some_context() - queue = cl.CommandQueue(ctx) - actx = PyOpenCLArrayContext(queue, force_device_scalars=True) - - tctx = toys.ToyContext( - actx.context, - # LaplaceKernel(2), - #YukawaKernel(2), extra_kernel_kwargs={"lam": 5}, - # HelmholtzKernel(2), extra_kernel_kwargs={"k": 0.3}, - HeatKernel(1), extra_kernel_kwargs={"alpha": 1}, - ) - - src_size = 0.1 - pt_src = toys.PointSources( - tctx, - np.array([ - src_size*(np.random.rand(50) - 0.5), - np.zeros(50)]), - np.random.randn(50)) - - fp = FieldPlotter([0, 0.5], extent=np.array([8, 1])) - - if 0 and USE_MATPLOTLIB: - toys.logplot(fp, pt_src, cmap="jet", aspect=8) - plt.colorbar() - plt.show() - - - p = 5 - center = np.array([0, 1], dtype=np.float64) - mexp = toys.local_expand(pt_src, center, p) - diff = mexp - pt_src - - dist = fp.points - center[:, None] - r = np.sqrt(dist[0]**2 + dist[1]**2) - - error_model = r**(p+1) - - def logplot_fp(fp: FieldPlotter, values, **kwargs) -> None: - fp.show_scalar_in_matplotlib( - np.log10(np.abs(values + 1e-15)), **kwargs) - if USE_MATPLOTLIB: - plt.subplot(131) - logplot_fp(fp, error_model, cmap="jet", vmin=-5, vmax=0, aspect=8) - plt.subplot(132) - logplot_fp(fp, diff.eval(fp.points), cmap="jet", vmin=-5, vmax=0, aspect=8) - plt.subplot(133) - logplot_fp(fp, diff.eval(fp.points)/error_model, cmap="jet", vmin=-5, vmax=0, aspect=8) - plt.colorbar() - plt.show() - 1/0 - mexp2 = toys.multipole_expand(mexp, [0, 0.25]) # noqa: F841 - lexp = toys.local_expand(mexp, [3, 0]) - lexp2 = toys.local_expand(lexp, [3, 1], 3) - - # diff = mexp - pt_src - # diff = mexp2 - pt_src - diff = lexp2 - pt_src - - print(toys.l_inf(diff, 1.2, center=lexp2.center)) - - -if __name__ == "__main__": - main() diff --git a/examples/heat-mpole.py b/examples/heat-mpole.py deleted file mode 100644 index 7bc4cc6ac..000000000 --- a/examples/heat-mpole.py +++ /dev/null @@ -1,87 +0,0 @@ -import numpy as np - -import pyopencl as cl - -from sumpy import toys -from sumpy.visualization import FieldPlotter -from sumpy.kernel import ( # noqa: F401 - YukawaKernel, - HeatKernel, - HelmholtzKernel, - LaplaceKernel) - -try: - import matplotlib.pyplot as plt - USE_MATPLOTLIB = True -except ImportError: - USE_MATPLOTLIB = False - - -def main(): - from sumpy.array_context import PyOpenCLArrayContext - ctx = cl.create_some_context() - queue = cl.CommandQueue(ctx) - actx = PyOpenCLArrayContext(queue, force_device_scalars=True) - - tctx = toys.ToyContext( - actx.context, - # LaplaceKernel(2), - #YukawaKernel(2), extra_kernel_kwargs={"lam": 5}, - # HelmholtzKernel(2), extra_kernel_kwargs={"k": 0.3}, - HeatKernel(1), extra_kernel_kwargs={"alpha": 1}, - ) - - src_size = 0.1 - pt_src = toys.PointSources( - tctx, - np.array([ - src_size*(np.random.rand(50) - 0.5), - np.zeros(50)]), - np.random.randn(50)) - - fp = FieldPlotter([0, 0.5], extent=np.array([8, 1])) - - if 0 and USE_MATPLOTLIB: - toys.logplot(fp, pt_src, cmap="jet", aspect=8) - plt.colorbar() - plt.show() - - - p = 5 - mexp = toys.multipole_expand(pt_src, [0, 0], p) - diff = mexp - pt_src - - x, t = fp.points - - delta = 4 * t - conv_factor = src_size / np.sqrt(2*delta) - - error_model = conv_factor**(p+1)*np.exp(-(x/np.sqrt(delta))**2/2) - #error_model = conv_factor**(p+1)/(1-conv_factor)*np.exp(-(x/np.sqrt(delta))**2) - - def logplot_fp(fp: FieldPlotter, values, **kwargs) -> None: - fp.show_scalar_in_matplotlib( - np.log10(np.abs(values + 1e-15)), **kwargs) - if USE_MATPLOTLIB: - plt.subplot(131) - logplot_fp(fp, error_model, cmap="jet", vmin=-5, vmax=0, aspect=8) - plt.subplot(132) - logplot_fp(fp, diff.eval(fp.points), cmap="jet", vmin=-5, vmax=0, aspect=8) - plt.subplot(133) - logplot_fp(fp, diff.eval(fp.points)/error_model, cmap="jet", vmin=-5, vmax=0, aspect=8) - plt.colorbar() - plt.show() - 1/0 - mexp2 = toys.multipole_expand(mexp, [0, 0.25]) # noqa: F841 - lexp = toys.local_expand(mexp, [3, 0]) - lexp2 = toys.local_expand(lexp, [3, 1], 3) - - # diff = mexp - pt_src - # diff = mexp2 - pt_src - diff = lexp2 - pt_src - - print(toys.l_inf(diff, 1.2, center=lexp2.center)) - - -if __name__ == "__main__": - main() diff --git a/sumpy/test/test_misc.py b/sumpy/test/test_misc.py index 990e3e3b1..80e9c35b4 100644 --- a/sumpy/test/test_misc.py +++ b/sumpy/test/test_misc.py @@ -148,9 +148,12 @@ def test_pde_check_kernels( extra_source_kwargs=knl_info.extra_kwargs) rng = np.random.default_rng(42) + source_points = rng.random(size=(dim, 50)) - 0.5 + if isinstance(knl_info.kernel, HeatKernel): + source_points[-1] += 0.5 pt_src = t.PointSources( tctx, - rng.random(size=(dim, 50)) - 0.5, + source_points, np.ones(50)) from pytools.convergence import EOCRecorder @@ -160,7 +163,7 @@ def test_pde_check_kernels( for h in [0.1, 0.05, 0.025]: if isinstance(knl_info.kernel, HeatKernel): - cp = CalculusPatch(np.array([0, 0, 0, 1])[-dim:], h=h, order=order) + cp = CalculusPatch(np.array([0, 0, 0, 2])[-dim:], h=h, order=order) else: cp = CalculusPatch(np.array([1, 0, 0])[:dim], h=h, order=order) pot = pt_src.eval(actx, cp.points) @@ -324,14 +327,7 @@ def dim(self): @pytest.mark.parametrize("case", P2E2E2P_TEST_CASES) -@pytest.mark.parametrize(("make_kernel", "extra_kernel_kwargs"), [ - (LaplaceKernel, {}), - (lambda dim: HeatKernel(dim - 1), {"alpha": 0.1})]) -def test_toy_p2e2e2p( - actx_factory: ArrayContextFactory, - case, - make_kernel: Callable[[int], Kernel], - extra_kernel_kwargs: dict[str, object]): +def test_toy_p2e2e2p(actx_factory: ArrayContextFactory, case): dim = case.dim src = case.source.reshape(dim, -1) @@ -350,14 +346,13 @@ def test_toy_p2e2e2p( raise ValueError( f"convergence factor not in valid range: {case_conv_factor}") - from sumpy.expansion import DefaultExpansionFactory + from sumpy.expansion import VolumeTaylorExpansionFactory actx = actx_factory() ctx = t.ToyContext( - make_kernel(dim), - expansion_factory=DefaultExpansionFactory(), - m2l_use_fft=case.m2l_use_fft, - extra_kernel_kwargs=extra_kernel_kwargs) + LaplaceKernel(dim), + expansion_factory=VolumeTaylorExpansionFactory(), + m2l_use_fft=case.m2l_use_fft) errors = [] From 3839d41d5a368e00d6a43ea2559376c129ccbd40 Mon Sep 17 00:00:00 2001 From: Shawn Lin Date: Thu, 23 Apr 2026 15:53:25 -0500 Subject: [PATCH 3/8] reset the m2l --- sumpy/expansion/m2l.py | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/sumpy/expansion/m2l.py b/sumpy/expansion/m2l.py index b3782c73c..5a77df51f 100644 --- a/sumpy/expansion/m2l.py +++ b/sumpy/expansion/m2l.py @@ -879,6 +879,10 @@ def loopy_translate(self, vector = translation_classes_dependent_data[icoeff_src] expr = toeplitz_first_row * vector + domains.append(f"{{[icoeff_src]: icoeff_tgt<=icoeff_src<{ncoeff_src} }}") + + expr = src_coeffs[icoeff_tgt] * translation_classes_dependent_data[icoeff_tgt] + insns = [ lp.Assignment( assignee=tgt_coeffs[icoeff_tgt], From 4739fd506ad0e847b8699090b059c6928c8e3381 Mon Sep 17 00:00:00 2001 From: Shawn Lin Date: Thu, 23 Apr 2026 16:01:33 -0500 Subject: [PATCH 4/8] reset toys.py --- sumpy/toys.py | 4 +--- 1 file changed, 1 insertion(+), 3 deletions(-) diff --git a/sumpy/toys.py b/sumpy/toys.py index b366951ec..3ccda9c82 100644 --- a/sumpy/toys.py +++ b/sumpy/toys.py @@ -385,8 +385,6 @@ def _e2e(actx: ArrayContext, actx.to_numpy(to_coeffs[1]), derived_from=psource, **expn_kwargs) -# }}} - def _m2l(actx: ArrayContext, psource, to_center, to_rscale, to_order, e2e, expn_class, expn_kwargs, @@ -491,8 +489,8 @@ def _m2l(actx: ArrayContext, return ret -# }}} +# }}} # {{{ potential source classes From ecc476a3d87fa4fb78ae53152aa8a2a630a53188 Mon Sep 17 00:00:00 2001 From: Shawn Lin Date: Tue, 12 May 2026 16:58:22 -0500 Subject: [PATCH 5/8] revert to sumpy --- .gitlab-ci.yml | 4 ---- pyproject.toml | 3 +-- 2 files changed, 1 insertion(+), 6 deletions(-) diff --git a/.gitlab-ci.yml b/.gitlab-ci.yml index 51e5a8f93..13bca8b85 100644 --- a/.gitlab-ci.yml +++ b/.gitlab-ci.yml @@ -22,7 +22,6 @@ Pytest POCL: - export PY_EXE=python3 - export PYOPENCL_TEST=portable:pthread - export EXTRA_INSTALL="pybind11 numpy mako mpi4py" - - export PYTEST_ADDOPTS=${PYTEST_ADDOPTS:-"-m 'not slowtest'"} - curl -L -O https://gitlab.tiker.net/inducer/ci-support/raw/main/build-and-test-py-project.sh - ". ./build-and-test-py-project.sh" tags: @@ -39,7 +38,6 @@ Pytest Titan V: - py_version=3 - export PYOPENCL_TEST=nvi:titan - EXTRA_INSTALL="pybind11 numpy mako mpi4py" - - export PYTEST_ADDOPTS=${PYTEST_ADDOPTS:-"-m 'not slowtest'"} - curl -L -O https://gitlab.tiker.net/inducer/ci-support/raw/main/build-and-test-py-project.sh - ". ./build-and-test-py-project.sh" tags: @@ -58,7 +56,6 @@ Pytest Conda: - export SUMPY_NO_CACHE=1 - export SUMPY_FORCE_SYMBOLIC_BACKEND=symengine - export PYOPENCL_TEST=portable:pthread - - export PYTEST_ADDOPTS=${PYTEST_ADDOPTS:-"-m 'not slowtest'"} - curl -L -O https://gitlab.tiker.net/inducer/ci-support/raw/main/build-and-test-py-project-within-miniconda.sh - ". ./build-and-test-py-project-within-miniconda.sh" tags: @@ -75,7 +72,6 @@ Pytest POCL Titan V: - export SUMPY_NO_CACHE=1 - export SUMPY_FORCE_SYMBOLIC_BACKEND=symengine - export PYOPENCL_TEST=portable:titan - - export PYTEST_ADDOPTS=${PYTEST_ADDOPTS:-"-m 'not slowtest'"} - curl -L -O https://gitlab.tiker.net/inducer/ci-support/raw/main/build-and-test-py-project-within-miniconda.sh - ". ./build-and-test-py-project-within-miniconda.sh" tags: diff --git a/pyproject.toml b/pyproject.toml index 7042587ff..af402d9ad 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -154,7 +154,6 @@ extend-exclude = [ [tool.pytest.ini_options] markers = [ "mpi: tests distributed FMM", - "slowtest: mark a test as slow", ] [tool.basedpyright] @@ -187,6 +186,7 @@ exclude = [ "benchmarks", "contrib", ".conda-root", + ".venv", ] [[tool.basedpyright.executionEnvironments]] @@ -224,4 +224,3 @@ reportUnknownMemberType = "hint" reportMissingImports = "none" reportArgumentType = "hint" reportOperatorIssue = "hint" - From a15d9492dcd310a93fb3df25e6ad1a871090ac76 Mon Sep 17 00:00:00 2001 From: Shawn Lin Date: Tue, 12 May 2026 16:59:07 -0500 Subject: [PATCH 6/8] test_kernels update + toy p2m, p2l --- sumpy/test/test_kernels.py | 213 +++++++++++-------------------------- sumpy/toys.py | 4 +- 2 files changed, 63 insertions(+), 154 deletions(-) diff --git a/sumpy/test/test_kernels.py b/sumpy/test/test_kernels.py index 996644829..01431aafa 100644 --- a/sumpy/test/test_kernels.py +++ b/sumpy/test/test_kernels.py @@ -37,7 +37,7 @@ PyOpenCLArrayContext, pytest_generate_tests_for_array_contexts, ) -from pytools import memoize_on_first_arg, obj_array +from pytools import memoize_on_first_arg from pytools.convergence import PConvergenceVerifier import sumpy.symbolic as sym @@ -64,7 +64,6 @@ BiharmonicKernel, DirectionalSourceDerivative, ElasticityKernel, - HeatKernel, HelmholtzKernel, Kernel, LaplaceKernel, @@ -267,15 +266,6 @@ def test_p2e_multiple( (LaplaceKernel(2), LinearPDEConformingVolumeTaylorLocalExpansion), (LaplaceKernel(2), LinearPDEConformingVolumeTaylorMultipoleExpansion), - (HeatKernel(1), LinearPDEConformingVolumeTaylorLocalExpansion), - (HeatKernel(1), LinearPDEConformingVolumeTaylorMultipoleExpansion), - (HeatKernel(2), LinearPDEConformingVolumeTaylorLocalExpansion), - (HeatKernel(2), LinearPDEConformingVolumeTaylorMultipoleExpansion), - pytest.param(HeatKernel(3), LinearPDEConformingVolumeTaylorLocalExpansion, - marks=pytest.mark.slowtest), - pytest.param(HeatKernel(3), LinearPDEConformingVolumeTaylorMultipoleExpansion, - marks=pytest.mark.slowtest), - (HelmholtzKernel(2), VolumeTaylorMultipoleExpansion), (HelmholtzKernel(2), VolumeTaylorLocalExpansion), (HelmholtzKernel(2), LinearPDEConformingVolumeTaylorLocalExpansion), @@ -320,25 +310,12 @@ def test_p2e2p( extra_kwargs["k"] = 0.2 if isinstance(base_knl, StokesletKernel): extra_kwargs["mu"] = 0.2 - if isinstance(base_knl, HeatKernel): - extra_kwargs["alpha"] = 1.0 if with_source_derivative: knl = DirectionalSourceDerivative(base_knl, "dir_vec") else: knl = base_knl - target_kernels = [ - knl, - AxisTargetDerivative(0, knl), - ] - expn = expn_class(knl, order=order) - - from sumpy import P2P, E2PFromSingleBox, P2EFromSingleBox - p2e = P2EFromSingleBox(expn, kernels=[knl]) - e2p = E2PFromSingleBox(expn, kernels=target_kernels) - p2p = P2P(target_kernels, exclude_self=False) - from pytools.convergence import EOCRecorder eoc_rec_pot = EOCRecorder() eoc_rec_grad_x = EOCRecorder() @@ -350,131 +327,81 @@ def test_p2e2p( h_values = [1/2, 1/3, 1/5] rng = np.random.default_rng(19) + center = np.array([2, 1, 0][:knl.dim], np.float64) + sources = ( + 0.7 * (-0.5 + rng.random((knl.dim, nsources), dtype=np.float64)) + + center[:, np.newaxis]) - if isinstance(base_knl, HeatKernel): - # Setup sources so that there are no negative time intervals - center = np.array([0, 0, 0, 0.5][-knl.dim:], np.float64) - sources = (-0.5 + rng.random((knl.dim, nsources), dtype=np.float64) - + center[:, np.newaxis]) - loc_center = np.array([0.0, 0.0, 0.0, 6.0][-knl.dim:]) + center - varying_axis = -1 - else: - center = np.array([2, 1, 0][:knl.dim], np.float64) - sources = 0.7 * (-0.5 + rng.random((knl.dim, nsources), dtype=np.float64) - + center[:, np.newaxis]) - loc_center = np.array([5.5, 0.0, 0.0][:knl.dim]) + center - varying_axis = 0 - - strengths = actx.from_numpy(np.ones(nsources, dtype=np.float64) / nsources) - sources = actx.from_numpy(sources) - - source_boxes = actx.from_numpy(np.array([0], dtype=np.int32)) - box_source_starts = actx.from_numpy(np.array([0], dtype=np.int32)) - box_source_counts_nonchild = ( - actx.from_numpy(np.array([nsources], dtype=np.int32))) + strengths = np.ones(nsources, dtype=np.float64) / nsources extra_source_kwargs = extra_kwargs.copy() if isinstance(knl, DirectionalSourceDerivative): alpha = np.linspace(0, 2*np.pi, nsources, np.float64) - dir_vec = np.vstack( - [np.cos(alpha), np.sin(alpha), alpha, alpha**2][:knl.dim]) + dir_vec = np.vstack([np.cos(alpha), np.sin(alpha)]) extra_source_kwargs["dir_vec"] = actx.from_numpy(dir_vec) + if issubclass(expn_class, LocalExpansionBase): + toy_ctx = t.ToyContext( + knl, + local_expn_class=expn_class, + extra_source_kwargs=extra_source_kwargs, + extra_kernel_kwargs=extra_kwargs, + ) + toy_ctx_grad_x = t.ToyContext( + AxisTargetDerivative(0, knl), + local_expn_class=expn_class, + extra_source_kwargs=extra_source_kwargs, + extra_kernel_kwargs=extra_kwargs, + ) + else: + toy_ctx = t.ToyContext( + knl, + mpole_expn_class=expn_class, + extra_source_kwargs=extra_source_kwargs, + extra_kernel_kwargs=extra_kwargs, + ) + toy_ctx_grad_x = t.ToyContext( + AxisTargetDerivative(0, knl), + mpole_expn_class=expn_class, + extra_source_kwargs=extra_source_kwargs, + extra_kernel_kwargs=extra_kwargs, + ) + + point_sources = t.PointSources(toy_ctx, sources, strengths) + point_sources_grad_x = t.PointSources(toy_ctx_grad_x, sources, strengths) + from sumpy.visualization import FieldPlotter for h in h_values: if issubclass(expn_class, LocalExpansionBase): - centers = np.array(loc_center, dtype=np.float64).reshape(knl.dim, 1) + loc_center = np.array([5.5, 0.0, 0.0][:knl.dim]) + center fp = FieldPlotter(loc_center, extent=h, npoints=res) else: - eval_center = center.copy() - eval_center[varying_axis] = 1/h + eval_center = np.array([1/h, 0.0, 0.0][:knl.dim]) + center fp = FieldPlotter(eval_center, extent=0.1, npoints=res) - centers = center[:, np.newaxis] - centers = actx.from_numpy(centers) - targets = actx.from_numpy(obj_array.new_1d(fp.points)) + targets = fp.points rscale = 0.5 # pick something non-1 - # {{{ apply p2e - - mpoles = p2e(actx, - source_boxes=source_boxes, - box_source_starts=box_source_starts, - box_source_counts_nonchild=box_source_counts_nonchild, - centers=centers, - sources=sources, - strengths=(strengths,), - nboxes=1, - tgt_base_ibox=0, - rscale=rscale, - **extra_source_kwargs) - - # }}} - - # {{{ apply e2p - - ntargets = fp.points.shape[-1] - - box_target_starts = actx.from_numpy(np.array([0], dtype=np.int32)) - box_target_counts_nonchild = ( - actx.from_numpy(np.array([ntargets], dtype=np.int32))) - - pot, grad_x = e2p( - actx, - src_expansions=mpoles, - src_base_ibox=0, - target_boxes=source_boxes, - box_target_starts=box_target_starts, - box_target_counts_nonchild=box_target_counts_nonchild, - centers=centers, - targets=targets, - rscale=rscale, - **extra_kwargs) - pot = actx.to_numpy(pot) - grad_x = actx.to_numpy(grad_x) - - # }}} - - # {{{ compute (direct) reference solution - - pot_direct, grad_x_direct = p2p( - actx, - targets, sources, (strengths,), - **extra_source_kwargs) - pot_direct = actx.to_numpy(pot_direct) - grad_x_direct = actx.to_numpy(grad_x_direct) - - err_pot = la.norm((pot - pot_direct)/res**2) - err_grad_x = la.norm((grad_x - grad_x_direct)/res**2) - - if 1: - err_pot = err_pot / la.norm((pot_direct)/res**2) - err_grad_x = err_grad_x / la.norm((grad_x_direct)/res**2) - - if 0: - import matplotlib.pyplot as pt - from matplotlib.colors import Normalize - - pt.subplot(131) - im = fp.show_scalar_in_matplotlib(pot.real) - im.set_norm(Normalize(vmin=-0.1, vmax=0.1)) - - pt.subplot(132) - im = fp.show_scalar_in_matplotlib(pot_direct.real) - im.set_norm(Normalize(vmin=-0.1, vmax=0.1)) - pt.colorbar() - - pt.subplot(133) - im = fp.show_scalar_in_matplotlib(np.log10(1e-15+np.abs(pot-pot_direct))) - im.set_norm(Normalize(vmin=-6, vmax=1)) - - pt.colorbar() - pt.show() - - # }}} - + if issubclass(expn_class, LocalExpansionBase): + expn = t.local_expand(actx, point_sources, loc_center, + order=order, rscale=rscale) + expn_grad_x = t.local_expand(actx, point_sources_grad_x, loc_center, + order=order, rscale=rscale) + else: + expn = t.multipole_expand(actx, point_sources, center, + order=order, rscale=rscale) + expn_grad_x = t.multipole_expand(actx, point_sources_grad_x, center, + order=order, rscale=rscale) + + pot = expn.eval(actx, targets) + pot_direct = point_sources.eval(actx, targets) + grad_x = expn_grad_x.eval(actx, targets) + grad_x_direct = point_sources_grad_x.eval(actx, targets) + + err_pot = la.norm(pot - pot_direct) / la.norm(pot_direct) + err_grad_x = la.norm(grad_x - grad_x_direct) / la.norm(grad_x_direct) eoc_rec_pot.add_data_point(h, err_pot) eoc_rec_grad_x.add_data_point(h, err_grad_x) @@ -488,10 +415,9 @@ def test_p2e2p( if issubclass(expn_class, LocalExpansionBase): tgt_order_grad = tgt_order - 1 slack = 0.7 - grad_slack = 0.8 if isinstance(base_knl, HeatKernel) else 0.5 + grad_slack = 0.5 else: tgt_order_grad = tgt_order + 1 - slack = 0.5 grad_slack = 1 @@ -499,10 +425,6 @@ def test_p2e2p( slack += 1 grad_slack += 1 - if isinstance(base_knl, HeatKernel): - slack += 1.0 - grad_slack += 2.5 - if isinstance(knl, DirectionalSourceDerivative): slack += 1 grad_slack += 2 @@ -529,14 +451,6 @@ def test_p2e2p( # {{{ test_translations @pytest.mark.parametrize("knl, local_expn_class, mpole_expn_class, use_fft", [ - (HeatKernel(1), LinearPDEConformingVolumeTaylorLocalExpansion, - LinearPDEConformingVolumeTaylorMultipoleExpansion, False), - (HeatKernel(2), LinearPDEConformingVolumeTaylorLocalExpansion, - LinearPDEConformingVolumeTaylorMultipoleExpansion, False), - (HeatKernel(2), LinearPDEConformingVolumeTaylorLocalExpansion, - LinearPDEConformingVolumeTaylorMultipoleExpansion, True), # not in original PR - (HeatKernel(3), LinearPDEConformingVolumeTaylorLocalExpansion, - LinearPDEConformingVolumeTaylorMultipoleExpansion, False), (LaplaceKernel(2), VolumeTaylorLocalExpansion, VolumeTaylorMultipoleExpansion, False), (LaplaceKernel(2), LinearPDEConformingVolumeTaylorLocalExpansion, @@ -583,15 +497,10 @@ def test_translations( extra_kwargs["k"] = 0.05 if isinstance(knl, StokesletKernel): extra_kwargs["mu"] = 0.05 - if isinstance(knl, HeatKernel): - extra_kwargs["alpha"] = 0.1 # Just to make sure things also work away from the origin rng = np.random.default_rng(18) - if isinstance(knl, HeatKernel): - origin = np.array([0, 0, 1, 2][-knl.dim:], np.float64) - else: - origin = np.array([2, 1, 0][:knl.dim], np.float64) + origin = np.array([2, 1, 0][:knl.dim], np.float64) sources = ( 0.7 * (-0.5 + rng.random((knl.dim, nsources), dtype=np.float64)) + origin[:, np.newaxis]) @@ -626,7 +535,7 @@ def test_translations( del eval_offset - if knl.dim == 2 and not isinstance(knl, HeatKernel): + if knl.dim == 2: orders = [2, 3, 4] else: orders = [3, 4, 5] diff --git a/sumpy/toys.py b/sumpy/toys.py index 3ccda9c82..6d77405fa 100644 --- a/sumpy/toys.py +++ b/sumpy/toys.py @@ -183,14 +183,14 @@ def get_p2m(self, order: int): from sumpy import P2EFromSingleBox return P2EFromSingleBox( self.mpole_expn_class(self.no_target_deriv_kernel, order), - kernels=(self.kernel,)) + kernels=(self.no_target_deriv_kernel,)) @memoize_method def get_p2l(self, order: int): from sumpy import P2EFromSingleBox return P2EFromSingleBox( self.local_expn_class(self.no_target_deriv_kernel, order), - kernels=(self.kernel,)) + kernels=(self.no_target_deriv_kernel,)) @memoize_method def get_m2p(self, order: int): From b6cbf844c67c1091cd9874845d142ad6b7d1aeb1 Mon Sep 17 00:00:00 2001 From: Shawn Lin Date: Tue, 12 May 2026 17:42:45 -0500 Subject: [PATCH 7/8] update heat kernel --- sumpy/test/test_heat_translations.py | 290 +++++++++++++++++++++++++++ 1 file changed, 290 insertions(+) create mode 100644 sumpy/test/test_heat_translations.py diff --git a/sumpy/test/test_heat_translations.py b/sumpy/test/test_heat_translations.py new file mode 100644 index 000000000..44369b329 --- /dev/null +++ b/sumpy/test/test_heat_translations.py @@ -0,0 +1,290 @@ +from __future__ import annotations + +import logging +import sys +from functools import partial + +import numpy as np +import numpy.linalg as la +import pytest + +from arraycontext import ArrayContextFactory, pytest_generate_tests_for_array_contexts +from pytools.convergence import EOCRecorder + +import sumpy.toys as t +from sumpy.array_context import PytestPyOpenCLArrayContextFactory, _acf # noqa: F401 +from sumpy.expansion.local import ( + LinearPDEConformingVolumeTaylorLocalExpansion, + VolumeTaylorLocalExpansion, +) +from sumpy.expansion.m2l import NonFFTM2LTranslationClassFactory +from sumpy.expansion.multipole import ( + LinearPDEConformingVolumeTaylorMultipoleExpansion, + VolumeTaylorMultipoleExpansion, +) +from sumpy.kernel import HeatKernel + + +logger = logging.getLogger(__name__) + +pytest_generate_tests = pytest_generate_tests_for_array_contexts([ + PytestPyOpenCLArrayContextFactory, + ]) + + +# {{{ test_heat_m2m + +@pytest.mark.parametrize("order", [4]) +@pytest.mark.parametrize("alpha", [1.0]) +@pytest.mark.parametrize(("knl", "mpole_expn_class"), [ + (HeatKernel(1), LinearPDEConformingVolumeTaylorMultipoleExpansion), + (HeatKernel(1), VolumeTaylorMultipoleExpansion), + ]) +def test_heat_m2m( + actx_factory: ArrayContextFactory, + knl, + mpole_expn_class, + alpha, + order): + # h-convergence of the heat-kernel M2M translation. + # src_center = (1, 0.5), shrinks with h (half-sizes h_x, h_t). + # tgt_center = (0, 2.5), fixed (half-sizes 1, 0.5). + # M2M from src_center to a shifted center c5 = src_center + (-h_x, h_t). + actx = actx_factory() + + extra_kwargs = {"alpha": alpha} + + t_sep = 1.0 + L = np.sqrt(4 * alpha * t_sep) + + src_center = np.array([1.0, 0.5]) + tgt_center = np.array([0.0, 2.5]) + + grid = np.linspace(-1.0, 1.0, 5) + xs, ts = np.meshgrid(grid, grid, indexing="xy") + targets = np.vstack([ + tgt_center[0] + 1.0 * xs.ravel(), + tgt_center[1] + 0.5 * ts.ravel(), + ]) + + toy_ctx = t.ToyContext( + kernel=knl, + mpole_expn_class=mpole_expn_class, + extra_kernel_kwargs=extra_kwargs, + ) + + h_values = [1/4, 1/8, 1/16, 1/32] + rscale = 1 / order + + eoc_rec = EOCRecorder() + for h in h_values: + h_x = h * L + h_t = h * t_sep + + src_grid = np.linspace(-1.0, 1.0, 11) + sxs, sts = np.meshgrid(src_grid, src_grid, indexing="xy") + sources = np.vstack([ + src_center[0] + h_x * sxs.ravel(), + src_center[1] + h_t * sts.ravel(), + ]) + strengths = np.ones(sources.shape[1]) + + pt_src = t.PointSources(toy_ctx, sources, weights=strengths) + + c5 = src_center + np.array([-h_x, h_t]) + + p2m = t.multipole_expand(actx, pt_src, src_center, + order=order, rscale=rscale) + m2m = t.multipole_expand(actx, p2m, c5, order=order, rscale=rscale) + p2m_direct = t.multipole_expand(actx, pt_src, c5, + order=order, rscale=rscale) + + m2m_vals = np.asarray(m2m.eval(actx, targets)).ravel() + p2m_direct_vals = np.asarray(p2m_direct.eval(actx, targets)).ravel() + err = la.norm(m2m_vals - p2m_direct_vals) / la.norm(p2m_direct_vals) + eoc_rec.add_data_point(h, err) + + logger.info("knl %s order %d", knl, order) + logger.info("M2M:\n%s", eoc_rec) + + tgt_order = order + 1 + slack = 0.5 + assert eoc_rec.order_estimate() > tgt_order - slack + +# }}} + + +# {{{ test_heat_l2l + +@pytest.mark.parametrize("order", [4]) +@pytest.mark.parametrize("alpha", [1.0]) +@pytest.mark.parametrize(("knl", "local_expn_class"), [ + (HeatKernel(1), LinearPDEConformingVolumeTaylorLocalExpansion), + (HeatKernel(1), VolumeTaylorLocalExpansion), + ]) +def test_heat_l2l( + actx_factory: ArrayContextFactory, + knl, + local_expn_class, + alpha, + order): + # h-convergence of the heat-kernel L2L translation. + # src_center = (0, 0.5), fixed (half-widths 1, 0.5). + # tgt_center = (0, 2.5), shrinks with h (half-sizes h_x, h_t). + # L2L from tgt_center to a shifted center c2 = tgt_center - (h_x, h_t). + actx = actx_factory() + + extra_kwargs = {"alpha": alpha} + + src_center = np.array([0.0, 0.5]) + src_hx = 1.0 + src_ht = 0.5 + tgt_center = np.array([0.0, 2.5]) + + t_sep = 2 * src_ht + L = np.sqrt(4 * alpha * t_sep) + + rng = np.random.default_rng(0) + src_grid = np.linspace(-1.0, 1.0, 11) + sxs, sts = np.meshgrid(src_grid, src_grid, indexing="xy") + sources = np.vstack([ + src_center[0] + src_hx * sxs.ravel(), + src_center[1] + src_ht * sts.ravel(), + ]) + strengths = rng.random(sources.shape[1]) + + toy_ctx = t.ToyContext( + kernel=knl, + local_expn_class=local_expn_class, + extra_kernel_kwargs=extra_kwargs, + ) + pt_src = t.PointSources(toy_ctx, sources, weights=strengths) + + rscale = 1 / order + p2l = t.local_expand(actx, pt_src, tgt_center, + order=order, rscale=rscale) + + h_values = [1/4, 1/8, 1/16, 1/32] + + eoc_rec = EOCRecorder() + for h in h_values: + h_x = h * L + h_t = h * t_sep + c2 = tgt_center - np.array([h_x, h_t]) + + l2l = t.local_expand(actx, p2l, c2, order=order, rscale=rscale) + p2l_direct = t.local_expand(actx, pt_src, c2, + order=order, rscale=rscale) + + tgt_grid = np.linspace(-1.0, 1.0, 5) + txs, tts = np.meshgrid(tgt_grid, tgt_grid, indexing="xy") + targets = np.vstack([ + tgt_center[0] + h_x * txs.ravel(), + tgt_center[1] + h_t * tts.ravel(), + ]) + l2l_vals = np.asarray(l2l.eval(actx, targets)).ravel() + p2l_direct_vals = np.asarray(p2l_direct.eval(actx, targets)).ravel() + err = la.norm(l2l_vals - p2l_direct_vals) / la.norm(p2l_direct_vals) + eoc_rec.add_data_point(h, err) + + logger.info("knl %s order %d", knl, order) + logger.info("L2L:\n%s", eoc_rec) + + tgt_order = order + 1 + slack = 0.5 + assert eoc_rec.order_estimate() > tgt_order - slack + +# }}} + + +# {{{ test_heat_m2l + +@pytest.mark.parametrize("order", [4]) +@pytest.mark.parametrize("alpha", [1.0]) +@pytest.mark.parametrize(("knl", "local_expn_class", "mpole_expn_class"), [ + (HeatKernel(1), + LinearPDEConformingVolumeTaylorLocalExpansion, + LinearPDEConformingVolumeTaylorMultipoleExpansion), + ]) +def test_heat_m2l( + actx_factory: ArrayContextFactory, + knl, + local_expn_class, + mpole_expn_class, + alpha, + order): + # h-convergence of the heat-kernel M2L translation. + # src_center = (0, 0.5), fixed (half-widths 1, 0.5). + # tgt_center = (0, 2.5), shrinks with h (half-sizes h_x, h_t). + # Local expansion formed at c2 = tgt_center - (h_x, h_t). + actx = actx_factory() + + extra_kwargs = {"alpha": alpha} + + src_center = np.array([0.0, 0.5]) + src_hx = 1.0 + src_ht = 0.5 + tgt_center = np.array([0.0, 2.5]) + + t_sep = 2 * src_ht + L = np.sqrt(4 * alpha * t_sep) + + rng = np.random.default_rng(0) + src_grid = np.linspace(-1.0, 1.0, 11) + sxs, sts = np.meshgrid(src_grid, src_grid, indexing="xy") + sources = np.vstack([ + src_center[0] + src_hx * sxs.ravel(), + src_center[1] + src_ht * sts.ravel(), + ]) + strengths = rng.random(sources.shape[1]) + + m2l_factory = NonFFTM2LTranslationClassFactory() + m2l_translation = m2l_factory.get_m2l_translation_class(knl, local_expn_class)() + toy_ctx = t.ToyContext( + kernel=knl, + local_expn_class=partial(local_expn_class, + m2l_translation_override=m2l_translation), + mpole_expn_class=mpole_expn_class, + extra_kernel_kwargs=extra_kwargs, + ) + pt_src = t.PointSources(toy_ctx, sources, weights=strengths) + + rscale = 1 / order + p2m = t.multipole_expand(actx, pt_src, src_center, + order=order, rscale=rscale) + + h_values = [1/4, 1/8, 1/16, 1/32] + + eoc_rec = EOCRecorder() + for h in h_values: + h_x = h * L + h_t = h * t_sep + c2 = tgt_center - np.array([h_x, h_t]) + + m2l = t.local_expand(actx, p2m, c2, order=order, rscale=rscale) + + tgt_grid = np.linspace(-1.0, 1.0, 5) + txs, tts = np.meshgrid(tgt_grid, tgt_grid, indexing="xy") + targets = np.vstack([ + tgt_center[0] + h_x * txs.ravel(), + tgt_center[1] + h_t * tts.ravel(), + ]) + m2l_vals = np.asarray(m2l.eval(actx, targets)).ravel() + p2m_vals = np.asarray(p2m.eval(actx, targets)).ravel() + err = la.norm(m2l_vals - p2m_vals) / la.norm(p2m_vals) + eoc_rec.add_data_point(h, err) + + logger.info("knl %s order %d", knl, order) + logger.info("M2L:\n%s", eoc_rec) + + tgt_order = order + 1 + slack = 0.5 + assert eoc_rec.order_estimate() > tgt_order - slack + +# }}} + +if __name__ == "__main__": + if len(sys.argv) > 1: + exec(sys.argv[1]) + else: + pytest.main([__file__]) From 7e3ca38679360efedf60b0847aef6d9f76c71489 Mon Sep 17 00:00:00 2001 From: Shawn Lin Date: Tue, 12 May 2026 17:46:34 -0500 Subject: [PATCH 8/8] update heat kernel --- sumpy/test/test_heat_translations.py | 7 ++++--- 1 file changed, 4 insertions(+), 3 deletions(-) diff --git a/sumpy/test/test_heat_translations.py b/sumpy/test/test_heat_translations.py index 44369b329..62bfcbbe2 100644 --- a/sumpy/test/test_heat_translations.py +++ b/sumpy/test/test_heat_translations.py @@ -55,7 +55,7 @@ def test_heat_m2m( extra_kwargs = {"alpha": alpha} t_sep = 1.0 - L = np.sqrt(4 * alpha * t_sep) + L = np.sqrt(4 * alpha * t_sep) # noqa: N806 src_center = np.array([1.0, 0.5]) tgt_center = np.array([0.0, 2.5]) @@ -142,7 +142,7 @@ def test_heat_l2l( tgt_center = np.array([0.0, 2.5]) t_sep = 2 * src_ht - L = np.sqrt(4 * alpha * t_sep) + L = np.sqrt(4 * alpha * t_sep) # noqa: N806 rng = np.random.default_rng(0) src_grid = np.linspace(-1.0, 1.0, 11) @@ -227,7 +227,7 @@ def test_heat_m2l( tgt_center = np.array([0.0, 2.5]) t_sep = 2 * src_ht - L = np.sqrt(4 * alpha * t_sep) + L = np.sqrt(4 * alpha * t_sep) # noqa: N806 rng = np.random.default_rng(0) src_grid = np.linspace(-1.0, 1.0, 11) @@ -283,6 +283,7 @@ def test_heat_m2l( # }}} + if __name__ == "__main__": if len(sys.argv) > 1: exec(sys.argv[1])