diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index a22d16eb..094ad829 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/sumpy/e2e.py b/sumpy/e2e.py index 0e5bdfd7..d4eb0b02 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/kernel.py b/sumpy/kernel.py index df9c0468..a3d1af5a 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_heat_translations.py b/sumpy/test/test_heat_translations.py new file mode 100644 index 00000000..62bfcbbe --- /dev/null +++ b/sumpy/test/test_heat_translations.py @@ -0,0 +1,291 @@ +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) # noqa: N806 + + 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) # noqa: N806 + + 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) # noqa: N806 + + 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__]) diff --git a/sumpy/test/test_kernels.py b/sumpy/test/test_kernels.py index 24fa7108..01431aaf 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 @@ -316,17 +316,6 @@ def test_p2e2p( 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() @@ -339,16 +328,11 @@ def test_p2e2p( rng = np.random.default_rng(19) center = np.array([2, 1, 0][:knl.dim], np.float64) - sources = actx.from_numpy( + sources = ( 0.7 * (-0.5 + rng.random((knl.dim, nsources), dtype=np.float64)) + center[:, np.newaxis]) - strengths = actx.from_numpy(np.ones(nsources, dtype=np.float64) / nsources) - - 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): @@ -356,103 +340,68 @@ def test_p2e2p( 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): 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 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 = 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) @@ -469,7 +418,6 @@ def test_p2e2p( grad_slack = 0.5 else: tgt_order_grad = tgt_order + 1 - slack = 0.5 grad_slack = 1 @@ -565,20 +513,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 +535,7 @@ def test_translations( del eval_offset - if knl.dim == 2: # noqa: SIM108 + if knl.dim == 2: orders = [2, 3, 4] else: orders = [3, 4, 5] diff --git a/sumpy/test/test_misc.py b/sumpy/test/test_misc.py index ce85f30b..80e9c35b 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 @@ -139,9 +148,12 @@ def test_pde_check_kernels(actx_factory: ArrayContextFactory, knl_info, order=5) 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 @@ -150,7 +162,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, 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) pde = knl_info.pde_func(cp, pot) diff --git a/sumpy/toys.py b/sumpy/toys.py index 3ccda9c8..6d77405f 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):