diff --git a/.basedpyright/baseline.json b/.basedpyright/baseline.json index 53900ca2..f474cceb 100644 --- a/.basedpyright/baseline.json +++ b/.basedpyright/baseline.json @@ -2419,22 +2419,6 @@ "lineCount": 1 } }, - { - "code": "reportUnknownMemberType", - "range": { - "startColumn": 50, - "endColumn": 70, - "lineCount": 1 - } - }, - { - "code": "reportUnknownArgumentType", - "range": { - "startColumn": 50, - "endColumn": 70, - "lineCount": 1 - } - }, { "code": "reportUnknownMemberType", "range": { @@ -2443,22 +2427,6 @@ "lineCount": 1 } }, - { - "code": "reportUnknownMemberType", - "range": { - "startColumn": 50, - "endColumn": 70, - "lineCount": 1 - } - }, - { - "code": "reportUnknownArgumentType", - "range": { - "startColumn": 50, - "endColumn": 70, - "lineCount": 1 - } - }, { "code": "reportUnknownMemberType", "range": { @@ -2467,22 +2435,6 @@ "lineCount": 1 } }, - { - "code": "reportUnknownMemberType", - "range": { - "startColumn": 50, - "endColumn": 70, - "lineCount": 1 - } - }, - { - "code": "reportUnknownArgumentType", - "range": { - "startColumn": 50, - "endColumn": 70, - "lineCount": 1 - } - }, { "code": "reportUnknownMemberType", "range": { @@ -5941,30 +5893,6 @@ "lineCount": 1 } }, - { - "code": "reportUnknownMemberType", - "range": { - "startColumn": 46, - "endColumn": 62, - "lineCount": 1 - } - }, - { - "code": "reportUnknownArgumentType", - "range": { - "startColumn": 46, - "endColumn": 62, - "lineCount": 1 - } - }, - { - "code": "reportUnknownArgumentType", - "range": { - "startColumn": 23, - "endColumn": 26, - "lineCount": 1 - } - }, { "code": "reportUnknownMemberType", "range": { @@ -10727,6 +10655,22 @@ "lineCount": 1 } }, + { + "code": "reportArgumentType", + "range": { + "startColumn": 25, + "endColumn": 48, + "lineCount": 1 + } + }, + { + "code": "reportArgumentType", + "range": { + "startColumn": 25, + "endColumn": 48, + "lineCount": 1 + } + }, { "code": "reportIncompatibleMethodOverride", "range": { @@ -10935,6 +10879,22 @@ "lineCount": 1 } }, + { + "code": "reportUnannotatedClassAttribute", + "range": { + "startColumn": 4, + "endColumn": 26, + "lineCount": 1 + } + }, + { + "code": "reportUnannotatedClassAttribute", + "range": { + "startColumn": 4, + "endColumn": 30, + "lineCount": 1 + } + }, { "code": "reportUnusedParameter", "range": { @@ -10991,6 +10951,22 @@ "lineCount": 1 } }, + { + "code": "reportUnannotatedClassAttribute", + "range": { + "startColumn": 4, + "endColumn": 26, + "lineCount": 1 + } + }, + { + "code": "reportUnannotatedClassAttribute", + "range": { + "startColumn": 4, + "endColumn": 30, + "lineCount": 1 + } + }, { "code": "reportUnannotatedClassAttribute", "range": { @@ -11129,22 +11105,6 @@ "lineCount": 1 } }, - { - "code": "reportUnknownMemberType", - "range": { - "startColumn": 50, - "endColumn": 66, - "lineCount": 1 - } - }, - { - "code": "reportUnknownArgumentType", - "range": { - "startColumn": 50, - "endColumn": 66, - "lineCount": 1 - } - }, { "code": "reportUnknownArgumentType", "range": { @@ -11971,14 +11931,6 @@ "lineCount": 1 } }, - { - "code": "reportUnknownArgumentType", - "range": { - "startColumn": 34, - "endColumn": 37, - "lineCount": 1 - } - }, { "code": "reportUnknownArgumentType", "range": { @@ -12699,14 +12651,6 @@ "lineCount": 1 } }, - { - "code": "reportUnknownArgumentType", - "range": { - "startColumn": 44, - "endColumn": 47, - "lineCount": 1 - } - }, { "code": "reportUnknownMemberType", "range": { @@ -12715,22 +12659,6 @@ "lineCount": 1 } }, - { - "code": "reportUnknownArgumentType", - "range": { - "startColumn": 50, - "endColumn": 53, - "lineCount": 1 - } - }, - { - "code": "reportUnknownArgumentType", - "range": { - "startColumn": 48, - "endColumn": 51, - "lineCount": 1 - } - }, { "code": "reportUnknownMemberType", "range": { @@ -12739,78 +12667,6 @@ "lineCount": 1 } }, - { - "code": "reportUnknownArgumentType", - "range": { - "startColumn": 42, - "endColumn": 45, - "lineCount": 1 - } - }, - { - "code": "reportUnknownArgumentType", - "range": { - "startColumn": 33, - "endColumn": 36, - "lineCount": 1 - } - }, - { - "code": "reportUnknownArgumentType", - "range": { - "startColumn": 33, - "endColumn": 36, - "lineCount": 1 - } - }, - { - "code": "reportUnknownArgumentType", - "range": { - "startColumn": 32, - "endColumn": 35, - "lineCount": 1 - } - }, - { - "code": "reportUnknownArgumentType", - "range": { - "startColumn": 33, - "endColumn": 36, - "lineCount": 1 - } - }, - { - "code": "reportUnknownArgumentType", - "range": { - "startColumn": 38, - "endColumn": 41, - "lineCount": 1 - } - }, - { - "code": "reportUnknownArgumentType", - "range": { - "startColumn": 57, - "endColumn": 60, - "lineCount": 1 - } - }, - { - "code": "reportUnknownArgumentType", - "range": { - "startColumn": 40, - "endColumn": 43, - "lineCount": 1 - } - }, - { - "code": "reportUnknownArgumentType", - "range": { - "startColumn": 40, - "endColumn": 43, - "lineCount": 1 - } - }, { "code": "reportUnknownMemberType", "range": { @@ -12819,14 +12675,6 @@ "lineCount": 1 } }, - { - "code": "reportUnknownArgumentType", - "range": { - "startColumn": 30, - "endColumn": 33, - "lineCount": 1 - } - }, { "code": "reportAny", "range": { @@ -22465,6 +22313,22 @@ "lineCount": 1 } }, + { + "code": "reportImplicitOverride", + "range": { + "startColumn": 8, + "endColumn": 16, + "lineCount": 1 + } + }, + { + "code": "reportUnknownArgumentType", + "range": { + "startColumn": 19, + "endColumn": 30, + "lineCount": 1 + } + }, { "code": "reportUnknownParameterType", "range": { @@ -23577,14 +23441,6 @@ "lineCount": 1 } }, - { - "code": "reportAny", - "range": { - "startColumn": 4, - "endColumn": 12, - "lineCount": 1 - } - }, { "code": "reportUnknownParameterType", "range": { diff --git a/sumpy/kernel.py b/sumpy/kernel.py index df9c0468..236d3b58 100644 --- a/sumpy/kernel.py +++ b/sumpy/kernel.py @@ -106,6 +106,12 @@ .. autoclass:: LineOfCompressionKernel :show-inheritance: :members: mapper_method +.. autoclass:: BrinkmanletKernel + :show-inheritance: + :members: mapper_method +.. autoclass:: BrinkmanStressKernel + :show-inheritance: + :members: mapper_method Derivatives ----------- @@ -724,7 +730,7 @@ def __reduce__(self) -> tuple[object, ...]: @property @override def is_complex_valued(self) -> bool: - # FIXME + # FIXME: 2D uses Hankel functions (complex-valued); 3D uses real exp return True @override @@ -1119,6 +1125,273 @@ def get_pde_as_diff_op(self) -> LinearPDESystemOperator: return laplacian(w) +@dataclass(frozen=True, repr=False) +class BrinkmanletKernel(ExpressionKernel): + r"""A kernel for the Brinkman equations. + + .. math:: + + \begin{cases} + -\mu (\Delta K_{ij}(\mathbf{x}, \mathbf{y}) + - k^2 K_{ij}(\mathbf{x}, \mathbf{y})) + + \nabla_i P_j(\mathbf{x}, \mathbf{y}) + = \delta_{ij}(\mathbf{x} - \mathbf{y}), \\ + \nabla_i K_{ij} = 0, + \end{cases} + + where :math:`P_j` is the pressure kernel. + """ + + mapper_method: ClassVar[str] = "map_brinkmanlet_kernel" + + icomp: int + """Component index for the kernel.""" + jcomp: int + """Component index for the kernel.""" + + viscosity_mu_name: str + """The argument name to use for the dynamic viscosity when generating + functions to evaluate this kernel. + """ + darcy_impermeability_name: str + """The argument name to use for the Darcy impermeability when generating + functions to evaluate this kernel. + """ + + def __init__(self, + dim: int, + icomp: int, + jcomp: int, + viscosity_mu_name: str = "mu", + darcy_impermeability_name: str = "k") -> None: + mu = SpatialConstant(viscosity_mu_name) + k = SpatialConstant(darcy_impermeability_name) + + d = make_sym_vector("d", dim) + r = pymbolic_real_norm_2(d) + R = k * r # noqa: N806 + delta_ij = 1 if icomp == jcomp else 0 + + # NOTE: + # [1] C. Pozrikidis, A Practical Guide to Boundary Element Methods, + # CRC Press, 2002. + # [2] https://dlmf.nist.gov/10.27#E8 + # [3] https://doi.org/10.1080/00036811.2011.614604 + # [4] https://doi.org/10.1002/mana.200710797 + + if dim == 2: + # transforming Bessel functions to Hankel functions using [2] + K0 = var("pi") * var("I") / 2 * var("hankel_1")(0, var("I") * R) # noqa: N806 + K1 = -var("pi") / 2 * var("hankel_1")(1, var("I") * R) # noqa: N806 + + # [1] Equations 7.7.5 and 7.7.6 + # [3] Equations 5.2 and 5.3 (for the scaling we use here) + a = 2 * (K0 + K1 / R - 1 / R**2) + b = 2 * (2 / R**2 - 2 * K1 / R - K0) + expr = a * delta_ij + b * d[icomp] * d[jcomp] / r ** 2 + scaling = -1 / (4 * var("pi") * mu) + elif dim == 3: + # [4] Equations 4.2 and 4.3 + a = 2 * var("exp")(-R) * (1 + 1 / R + 1 / R**2) - 2 / R**2 + b = 6 / R**2 - 2 * var("exp")(-R) * (1 + 3 / R + 3 / R**2) + expr = a * delta_ij / r + b * d[icomp] * d[jcomp] / r ** 3 + scaling = -1 / (8 * var("pi") * mu) + else: + raise NotImplementedError(f"unsupported dimension: '{dim}'") + + from pymbolic.mapper.flattener import flatten + + super().__init__(dim, + expression=flatten(expr), + global_scaling_const=flatten(scaling)) + + object.__setattr__(self, "icomp", icomp) + object.__setattr__(self, "jcomp", jcomp) + object.__setattr__(self, "viscosity_mu_name", viscosity_mu_name) + object.__setattr__(self, "darcy_impermeability_name", darcy_impermeability_name) + + @override + def __reduce__(self) -> tuple[object, ...]: + return ( + self.__class__, + (self.dim, self.icomp, self.jcomp, + self.viscosity_mu_name, self.darcy_impermeability_name), + ) + + @property + @override + def is_complex_valued(self) -> bool: + # FIXME: 2D uses Hankel functions (complex-valued); 3D uses real exp + return self.dim == 2 + + @override + def __str__(self) -> str: + return ( + f"BrinkmanletKnl{self.dim}D_{self.icomp}{self.jcomp}" + f"({self.viscosity_mu_name}, {self.darcy_impermeability_name})") + + @override + def prepare_loopy_kernel(self, loopy_knl: lp.TranslationUnit) -> lp.TranslationUnit: + from sumpy.codegen import register_bessel_callables + return register_bessel_callables(loopy_knl) + + @memoize_method + @override + def get_args(self) -> Sequence[KernelArgument]: + return [ + KernelArgument(loopy_arg=lp.ValueArg(self.viscosity_mu_name, np.float64)), + KernelArgument(loopy_arg=lp.ValueArg(self.darcy_impermeability_name, np.float64)), # noqa: E501 + ] + + @override + def get_pde_as_diff_op(self) -> LinearPDESystemOperator: + from sumpy.expansion.diff_op import laplacian, make_identity_diff_op + + w = make_identity_diff_op(self.dim) + k = sym.Symbol(self.darcy_impermeability_name) + return laplacian(laplacian(w) - k**2 * w) + + +@dataclass(frozen=True, repr=False) +class BrinkmanStressKernel(ExpressionKernel): + r"""A kernel for the Brinkman equations. + + .. math:: + + K_{ijk} = + -p_j \delta_{ik} + + \mu (\partial_k K_{ij} + \partial_i K_{jk}), + + where the two-index :math:`K_{ij}` is the :class:`~sumpy.kernel.BrinkmanletKernel` + and :math:`P_j` is the pressure kernel. + """ + + mapper_method: ClassVar[str] = "map_brinkman_stress_kernel" + + icomp: int + """Component index for the kernel.""" + jcomp: int + """Component index for the kernel.""" + kcomp: int + """Component index for the kernel.""" + + # NOTE: we keep the viscosity_mu_name for symmetry with the BrinkmanKernel, + # even though it isn't actually used in the stress kernel + viscosity_mu_name: str + """The argument name to use for the dynamic viscosity when generating + functions to evaluate this kernel. + """ + darcy_impermeability_name: str + """The argument name to use for the Darcy impermeability when generating + functions to evaluate this kernel. + """ + + def __init__(self, + dim: int, + icomp: int, + jcomp: int, + kcomp: int, + viscosity_mu_name: str = "mu", + darcy_impermeability_name: str = "k") -> None: + k = SpatialConstant(darcy_impermeability_name) + + d = make_sym_vector("d", dim) + r = pymbolic_real_norm_2(d) + R = k * r # noqa: N806 + delta_ij = 1 if icomp == jcomp else 0 + delta_ik = 1 if icomp == kcomp else 0 + delta_kj = 1 if jcomp == kcomp else 0 + + # NOTE: + # [1] C. Pozrikidis, A Practical Guide to Boundary Element Methods, + # CRC Press, 2002. + # [2] https://dlmf.nist.gov/10.27#E8 + # [3] https://doi.org/10.1080/00036811.2011.614604 + # [4] https://doi.org/10.1002/mana.200710797 + + if dim == 2: + # transforming Bessel functions to Hankel functions using [2] + K0 = var("pi") * var("I") / 2 * var("hankel_1")(0, var("I") * R) # noqa: N806 + K1 = -var("pi") / 2 * var("hankel_1")(1, var("I") * R) # noqa: N806 + + # [1] Equations 7.7.7 and 7.7.8 + # [3] Equations 5.4-5.6 (for the scaling we use here) + a = 2 * (2 / R**2 - 2 * K1 / R - K0) + b = 8 / R**2 - 4*K0 - 2*(R + 4 / R)*K1 + c = b + R * K1 + expr = ( + 2 * (a - 1) * d[jcomp] / r**2 * delta_ik + + b * (d[kcomp] * delta_ij + d[icomp] * delta_kj) / r**2 + - 4 * c * d[icomp] * d[jcomp] * d[kcomp] / r**4 + ) + scaling = -1 / (4 * var("pi")) + elif dim == 3: + # [4] Equations 4.4-4.6 + d1 = 2 * var("exp")(-R) * (1 + 3 / R + 3 / R**2) - 6 / R**2 + 1 + d2 = var("exp")(-R) * (R + 3 + 6 / R + 6 / R**2) - 6 / R**2 + d3 = var("exp")(-R) * (-2 * R - 12 - 30 / R - 30 / R**2) + 30 / R**2 + expr = ( + d1 * d[jcomp] / r**3 * delta_ik + + d2 * (d[kcomp] * delta_ij + d[icomp] * delta_kj) / r**3 + + d3 * d[icomp] * d[jcomp] * d[kcomp] / r**5 + ) + scaling = 1/(4*var("pi")) + else: + raise NotImplementedError(f"unsupported dimension: '{dim}'") + + from pymbolic.mapper.flattener import flatten + + super().__init__(dim, + expression=flatten(expr), + global_scaling_const=flatten(scaling)) + + object.__setattr__(self, "icomp", icomp) + object.__setattr__(self, "jcomp", jcomp) + object.__setattr__(self, "kcomp", kcomp) + object.__setattr__(self, "viscosity_mu_name", viscosity_mu_name) + object.__setattr__(self, "darcy_impermeability_name", darcy_impermeability_name) + + @override + def __reduce__(self) -> tuple[object, ...]: + return ( + self.__class__, + (self.dim, self.icomp, self.jcomp, self.kcomp, + self.viscosity_mu_name, self.darcy_impermeability_name), + ) + + @property + @override + def is_complex_valued(self) -> bool: + # 2D uses Hankel functions (complex-valued); 3D uses real exp + return self.dim == 2 + + @override + def __str__(self) -> str: + return ( + f"BrinkmanStressKnl{self.dim}D_{self.icomp}{self.jcomp}{self.kcomp}" + f"({self.viscosity_mu_name}, {self.darcy_impermeability_name})") + + @override + def prepare_loopy_kernel(self, loopy_knl: lp.TranslationUnit) -> lp.TranslationUnit: + from sumpy.codegen import register_bessel_callables + return register_bessel_callables(loopy_knl) + + @memoize_method + @override + def get_args(self) -> Sequence[KernelArgument]: + return [ + KernelArgument(loopy_arg=lp.ValueArg(self.viscosity_mu_name, np.float64)), + KernelArgument(loopy_arg=lp.ValueArg(self.darcy_impermeability_name, np.float64)), # noqa: E501 + ] + + @override + def get_pde_as_diff_op(self) -> LinearPDESystemOperator: + from sumpy.expansion.diff_op import laplacian, make_identity_diff_op + + w = make_identity_diff_op(self.dim) + k = sym.Symbol(self.darcy_impermeability_name) + return laplacian(laplacian(w) - k**2 * w) + # }}} @@ -1590,6 +1863,8 @@ 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_brinkmanlet_kernel = map_expression_kernel + map_brinkman_stress_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 +1937,8 @@ 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_brinkmanlet_kernel = map_expression_kernel + map_brinkman_stress_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 24fa7108..297bfe09 100644 --- a/sumpy/test/test_kernels.py +++ b/sumpy/test/test_kernels.py @@ -62,6 +62,8 @@ from sumpy.kernel import ( AxisTargetDerivative, BiharmonicKernel, + BrinkmanletKernel, + BrinkmanStressKernel, DirectionalSourceDerivative, ElasticityKernel, HelmholtzKernel, @@ -528,6 +530,10 @@ def test_p2e2p( VolumeTaylorMultipoleExpansion, False), (StokesletKernel(2, 0, 0), LinearPDEConformingVolumeTaylorLocalExpansion, LinearPDEConformingVolumeTaylorMultipoleExpansion, False), + (BrinkmanletKernel(2, 0, 0), VolumeTaylorLocalExpansion, + VolumeTaylorMultipoleExpansion, False), + (BrinkmanStressKernel(2, 0, 0, 0), VolumeTaylorLocalExpansion, + VolumeTaylorMultipoleExpansion, False), ]) def test_translations( actx_factory: ArrayContextFactory, @@ -547,8 +553,16 @@ def test_translations( extra_kwargs = {} if isinstance(knl, HelmholtzKernel): extra_kwargs["k"] = 0.05 - if isinstance(knl, StokesletKernel): + elif isinstance(knl, YukawaKernel): + extra_kwargs["lam"] = 0.05 + elif isinstance(knl, (StokesletKernel, StressletKernel)): extra_kwargs["mu"] = 0.05 + elif isinstance(knl, ElasticityKernel): + extra_kwargs["mu"] = 0.05 + extra_kwargs["nu"] = 0.1 + elif isinstance(knl, (BrinkmanletKernel, BrinkmanStressKernel)): + extra_kwargs["mu"] = 0.1 + extra_kwargs["k"] = 0.1 # Just to make sure things also work away from the origin rng = np.random.default_rng(18) @@ -587,8 +601,11 @@ def test_translations( del eval_offset - if knl.dim == 2: # noqa: SIM108 - orders = [2, 3, 4] + if knl.dim == 2: + if isinstance(knl, BrinkmanStressKernel): # noqa: SIM108 + orders = [3, 4, 5] + else: + orders = [2, 3, 4] else: orders = [3, 4, 5] @@ -944,6 +961,14 @@ def get_kernel_name_for_test(knl: Kernel) -> Callable[[str], str]: @pytest.mark.parametrize("knl", [ BiharmonicKernel(2), + BrinkmanletKernel(2, 0, 1), + BrinkmanletKernel(3, 0, 0, + viscosity_mu_name="viscosity", + darcy_impermeability_name="kappa"), + BrinkmanStressKernel(2, 0, 1, 0), + BrinkmanStressKernel(3, 0, 0, 1, + viscosity_mu_name="viscosity", + darcy_impermeability_name="kappa"), ElasticityKernel(2, 0, 0), HelmholtzKernel(3, helmholtz_k_name="kay"), LaplaceKernel(3), diff --git a/sumpy/test/test_misc.py b/sumpy/test/test_misc.py index ce85f30b..df3fc208 100644 --- a/sumpy/test/test_misc.py +++ b/sumpy/test/test_misc.py @@ -57,6 +57,8 @@ ) from sumpy.kernel import ( BiharmonicKernel, + BrinkmanletKernel, + BrinkmanStressKernel, ElasticityKernel, ExpressionKernel, HelmholtzKernel, @@ -91,6 +93,9 @@ def __init__(self, kernel, **kwargs): eq = diff_op.eqs[0] self.eq = eq + def __repr__(self) -> str: + return str(self.kernel) + def pde_func(self, cp, pot): subs_dict = {sym.Symbol(k): v for k, v in self.extra_kwargs.items()} result = 0 @@ -130,7 +135,15 @@ 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(BrinkmanletKernel(2, 0, 1), mu=5, k=3), + KernelInfo(BrinkmanletKernel(2, 1, 1), mu=5, k=3), + KernelInfo(BrinkmanletKernel(3, 0, 1), mu=5, k=3), + KernelInfo(BrinkmanletKernel(3, 1, 1), mu=5, k=3), + KernelInfo(BrinkmanStressKernel(2, 0, 1, 0), mu=5, k=3), + KernelInfo(BrinkmanStressKernel(2, 0, 1, 1), mu=5, k=3), + KernelInfo(BrinkmanStressKernel(3, 0, 1, 0), mu=5, k=3), + KernelInfo(BrinkmanStressKernel(3, 0, 1, 2), mu=5, k=3), + ], ids=repr) def test_pde_check_kernels(actx_factory: ArrayContextFactory, knl_info, order=5): actx = actx_factory()