diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index 6f3dd2f1..df2e7114 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 13bca8b8..51e5a8f9 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 00000000..fca53353 --- /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 00000000..7bc4cc6a --- /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 f3cd75e6..c4749e52 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -149,6 +149,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 34b888de..1065cf8c 100644 --- a/sumpy/e2e.py +++ b/sumpy/e2e.py @@ -324,7 +324,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 209b4947..b04b93f5 100644 --- a/sumpy/expansion/m2l.py +++ b/sumpy/expansion/m2l.py @@ -885,10 +885,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 c7a71f5c..250aa4e7 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 ----------- @@ -1084,6 +1087,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 + + # }}} @@ -1554,6 +1627,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)) @@ -1626,6 +1700,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 2167a56a..0576b5df 100644 --- a/sumpy/test/test_kernels.py +++ b/sumpy/test/test_kernels.py @@ -62,6 +62,7 @@ AxisTargetDerivative, BiharmonicKernel, DirectionalSourceDerivative, + HeatKernel, HelmholtzKernel, Kernel, LaplaceKernel, @@ -258,6 +259,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), @@ -302,6 +312,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") @@ -330,12 +342,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)) @@ -345,22 +368,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)) @@ -458,7 +480,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 @@ -469,6 +491,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 @@ -495,6 +521,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, @@ -541,10 +575,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]) @@ -557,20 +596,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 @@ -579,7 +618,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 41bf03b2..fc62efc3 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(("kernel_cls", "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 032b74f3..4f7648c4 100644 --- a/sumpy/toys.py +++ b/sumpy/toys.py @@ -386,6 +386,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, @@ -490,9 +492,9 @@ def _m2l(actx: ArrayContext, return ret - # }}} + # {{{ potential source classes Number_ish = int | float | complex | np.number