diff --git a/.basedpyright/baseline.json b/.basedpyright/baseline.json index 42f0e464..bd6079ca 100644 --- a/.basedpyright/baseline.json +++ b/.basedpyright/baseline.json @@ -9241,14 +9241,6 @@ "lineCount": 1 } }, - { - "code": "reportArgumentType", - "range": { - "startColumn": 25, - "endColumn": 48, - "lineCount": 1 - } - }, { "code": "reportReturnType", "range": { @@ -9313,6 +9305,86 @@ "lineCount": 1 } }, + { + "code": "reportUnknownArgumentType", + "range": { + "startColumn": 30, + "endColumn": 73, + "lineCount": 1 + } + }, + { + "code": "reportUnknownLambdaType", + "range": { + "startColumn": 37, + "endColumn": 38, + "lineCount": 1 + } + }, + { + "code": "reportUnknownLambdaType", + "range": { + "startColumn": 40, + "endColumn": 41, + "lineCount": 1 + } + }, + { + "code": "reportUnknownLambdaType", + "range": { + "startColumn": 43, + "endColumn": 44, + "lineCount": 1 + } + }, + { + "code": "reportUnknownLambdaType", + "range": { + "startColumn": 46, + "endColumn": 73, + "lineCount": 1 + } + }, + { + "code": "reportUnknownArgumentType", + "range": { + "startColumn": 52, + "endColumn": 53, + "lineCount": 1 + } + }, + { + "code": "reportUnknownArgumentType", + "range": { + "startColumn": 55, + "endColumn": 56, + "lineCount": 1 + } + }, + { + "code": "reportUnknownArgumentType", + "range": { + "startColumn": 63, + "endColumn": 64, + "lineCount": 1 + } + }, + { + "code": "reportUnknownArgumentType", + "range": { + "startColumn": 66, + "endColumn": 67, + "lineCount": 1 + } + }, + { + "code": "reportReturnType", + "range": { + "startColumn": 15, + "endColumn": 78, + "lineCount": 1 + } + }, { "code": "reportUnknownArgumentType", "range": { @@ -22517,6 +22589,54 @@ "lineCount": 1 } }, + { + "code": "reportAttributeAccessIssue", + "range": { + "startColumn": 58, + "endColumn": 63, + "lineCount": 1 + } + }, + { + "code": "reportUnknownMemberType", + "range": { + "startColumn": 20, + "endColumn": 33, + "lineCount": 1 + } + }, + { + "code": "reportAttributeAccessIssue", + "range": { + "startColumn": 28, + "endColumn": 33, + "lineCount": 1 + } + }, + { + "code": "reportUnknownMemberType", + "range": { + "startColumn": 35, + "endColumn": 48, + "lineCount": 1 + } + }, + { + "code": "reportAttributeAccessIssue", + "range": { + "startColumn": 43, + "endColumn": 48, + "lineCount": 1 + } + }, + { + "code": "reportUnknownMemberType", + "range": { + "startColumn": 50, + "endColumn": 63, + "lineCount": 1 + } + }, { "code": "reportAttributeAccessIssue", "range": { diff --git a/sumpy/kernel.py b/sumpy/kernel.py index 55183331..5342c8d3 100644 --- a/sumpy/kernel.py +++ b/sumpy/kernel.py @@ -104,6 +104,8 @@ :show-inheritance: :members: mapper_method +.. autoclass:: StokesComponentKernelBase + :show-inheritance: .. autoclass:: StokesletComponentKernel :show-inheritance: :members: mapper_method @@ -111,13 +113,20 @@ :show-inheritance: :members: mapper_method +.. autoclass:: ElasticityComponentKernelBase + :show-inheritance: .. autoclass:: ElasticityComponentKernel :show-inheritance: :members: mapper_method +.. autoclass:: ElasticityStressComponentKernel + :show-inheritance: + :members: mapper_method .. autoclass:: LineOfCompressionKernel :show-inheritance: :members: mapper_method +.. autoclass:: BrinkmanComponentKernelBase + :show-inheritance: .. autoclass:: BrinkmanletComponentKernel :show-inheritance: :members: mapper_method @@ -134,6 +143,9 @@ .. autoclass:: ElasticitySystemKernel :show-inheritance: :members: mapper_method +.. autoclass:: ElasticityStressSystemKernel + :show-inheritance: + :members: mapper_method .. autoclass:: StokesletSystemKernel :show-inheritance: @@ -929,8 +941,47 @@ def get_pde_as_diff_op(self) -> LinearPDESystemOperator: @dataclass(frozen=True, repr=False) -class ElasticityComponentKernel(ExpressionKernel): - r"""A kernel for the linear elasticity (Navier-Cauchy) equations +class ElasticityComponentKernelBase(ExpressionKernel): + r"""Base kernel class for the linear elasticity (Navier-Cauchy) equations + (see e.g. Section 2.2 in [Hsiao2008]_). + + .. autoattribute:: viscosity_mu_name + .. autoattribute:: poisson_ratio_name + """ + + viscosity_mu_name: str + r"""The argument name to use for the dynamic viscosity :math:`\mu` when + generating functions to evaluate this kernel. + """ + poisson_ratio_name: str + r"""The argument name to use for Poisson's ratio :math:`\nu` when generating + functions to evaluate this kernel. + """ + + @property + @override + def is_complex_valued(self) -> bool: + return False + + @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.poisson_ratio_name, np.float64)), + ] + + @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) + return laplacian(laplacian(w)) + + +@dataclass(frozen=True, repr=False) +class ElasticityComponentKernel(ElasticityComponentKernelBase): + r"""The displacement kernel for the linear elasticity (Navier-Cauchy) equations (see e.g. Section 2.2 in [Hsiao2008]_). .. math:: @@ -941,8 +992,6 @@ class ElasticityComponentKernel(ExpressionKernel): .. autoattribute:: icomp .. autoattribute:: jcomp - .. autoattribute:: viscosity_mu_name - .. autoattribute:: poisson_ratio_name """ mapper_method: ClassVar[str] = "map_elasticity_kernel" @@ -952,15 +1001,6 @@ class ElasticityComponentKernel(ExpressionKernel): jcomp: int """Component index for the kernel.""" - viscosity_mu_name: str - r"""The argument name to use for the dynamic viscosity :math:`\mu` when - generating functions to evaluate this kernel. - """ - poisson_ratio_name: str - r"""The argument name to use for Poisson's ratio :math:`\nu` when generating - functions to evaluate this kernel. - """ - def __init__(self, dim: int, icomp: int, @@ -995,12 +1035,12 @@ def __init__(self, else: raise NotImplementedError(f"unsupported dimension: '{dim}'") - super().__init__(dim, expression=expr, global_scaling_const=scaling) + super().__init__(dim, expression=expr, global_scaling_const=scaling, + viscosity_mu_name=viscosity_mu_name, + poisson_ratio_name=poisson_ratio_name) object.__setattr__(self, "icomp", icomp) object.__setattr__(self, "jcomp", jcomp) - object.__setattr__(self, "viscosity_mu_name", viscosity_mu_name) - object.__setattr__(self, "poisson_ratio_name", poisson_ratio_name) @override def __str__(self) -> str: @@ -1016,6 +1056,126 @@ def __reduce__(self): self.viscosity_mu_name, self.poisson_ratio_name)) + @override + @memoize_method + def get_pde_system_kernel(self) -> tuple[SystemKernel, tuple[int, ...]]: + return ElasticitySystemKernel( + self.dim, + viscosity_mu_name=self.viscosity_mu_name, + poisson_ratio_name=self.poisson_ratio_name + ), (self.icomp, self.jcomp) + + +@dataclass(frozen=True, repr=False) +class ElasticityStressComponentKernel(ElasticityComponentKernelBase): + r"""The stress kernel for the linear elasticity (Navier-Cauchy) equations + (see e.g. Section 2.2 in [Hsiao2008]_). + + .. math:: + + K_{ijk}(\mathbf{x}, \mathbf{y}) = + \lambda \partial_l K_{kl}(\mathbf{x}, \mathbf{y}) \delta_{ij} + + \mu (\partial_j K_{ik}(\mathbf{x}, \mathbf{y}) + + \partial_i K_{jk}(\mathbf{x}, \mathbf{y})), + + where the two-index :math:`K_{ij}` is the + :class:`~sumpy.kernel.ElasticityComponentKernel`. + + .. autoattribute:: icomp + .. autoattribute:: jcomp + .. autoattribute:: kcomp + """ + + mapper_method: ClassVar[str] = "map_elasticity_stress_kernel" + + icomp: int + """Component index for the kernel.""" + jcomp: int + """Component index for the kernel.""" + kcomp: int + """Component index for the kernel.""" + + def __init__(self, + dim: int, + icomp: int, + jcomp: int, + kcomp: int, + viscosity_mu_name: str = "mu", + poisson_ratio_name: str = "nu") -> None: + nu = sym.SpatialConstant(poisson_ratio_name) + + d = make_sym_vector("d", dim) + r = sym.pymbolic_real_norm_2(d) + delta_ij = 1 if icomp == jcomp else 0 + delta_ik = 1 if icomp == kcomp else 0 + delta_jk = 1 if jcomp == kcomp else 0 + + if dim == 2: + expr = ( + (1 - 2*nu) * ( + d[icomp] / r**2 * delta_jk + + d[jcomp] / r**2 * delta_ik + - d[kcomp] / r**2 * delta_ij) + + 3 * d[icomp] * d[jcomp] * d[kcomp] / r**4 + ) + scaling = -1/(4*var("pi")*(1 - nu)) + elif dim == 3: + expr = ( + (1 - 2*nu) * ( + d[icomp] / r**3 * delta_jk + + d[jcomp] / r**3 * delta_ik + - d[kcomp] / r**3 * delta_ij) + + 3 * d[icomp] * d[jcomp] * d[kcomp] / r**5 + ) + scaling = -1/(8*var("pi")*(1 - nu)) + else: + raise NotImplementedError(f"unsupported dimension: '{dim}'") + + super().__init__(dim, expression=expr, global_scaling_const=scaling, + viscosity_mu_name=viscosity_mu_name, + poisson_ratio_name=poisson_ratio_name) + + object.__setattr__(self, "icomp", icomp) + object.__setattr__(self, "jcomp", jcomp) + object.__setattr__(self, "kcomp", kcomp) + + @override + def __str__(self) -> str: + return ( + f"ElasticityStressKnl{self.dim}D_{self.icomp}{self.jcomp}{self.kcomp}" + f"({self.viscosity_mu_name}, {self.poisson_ratio_name})") + + @override + def __reduce__(self): + return ( + type(self), + (self.dim, self.icomp, self.jcomp, self.kcomp, + self.viscosity_mu_name, + self.poisson_ratio_name)) + + @override + @memoize_method + def get_pde_system_kernel(self) -> tuple[SystemKernel, tuple[int, ...]]: + return ElasticityStressSystemKernel( + self.dim, + viscosity_mu_name=self.viscosity_mu_name, + poisson_ratio_name=self.poisson_ratio_name + ), (self.icomp, self.jcomp, self.kcomp) + + +@dataclass(frozen=True, repr=False) +class StokesComponentKernelBase(ExpressionKernel): + """Base class for kernels of the Stokes equations + (see e.g. Chapter 2 in [Pozrikidis1992]_). + + .. autoattribute:: viscosity_mu_name + """ + + viscosity_mu_name: str + r"""The argument name to use for the dynamic viscosity :math:`\mu` when + generating functions to evaluate this kernel. + """ + @property @override def is_complex_valued(self) -> bool: @@ -1026,18 +1186,8 @@ def is_complex_valued(self) -> bool: def get_args(self) -> Sequence[KernelArgument]: return [ KernelArgument(loopy_arg=lp.ValueArg(self.viscosity_mu_name, np.float64)), - KernelArgument(loopy_arg=lp.ValueArg(self.poisson_ratio_name, np.float64)), ] - @override - @memoize_method - def get_pde_system_kernel(self) -> tuple[SystemKernel, tuple[int, ...]]: - return ElasticitySystemKernel( - self.dim, - viscosity_mu_name=self.viscosity_mu_name, - poisson_ratio_name=self.poisson_ratio_name - ), (self.icomp, self.jcomp) - @override def get_pde_as_diff_op(self) -> LinearPDESystemOperator: from sumpy.expansion.diff_op import laplacian, make_identity_diff_op @@ -1047,8 +1197,9 @@ def get_pde_as_diff_op(self) -> LinearPDESystemOperator: @dataclass(frozen=True, repr=False) -class StokesletComponentKernel(ExpressionKernel): - r"""A kernel for the Stokes equations (see e.g. Chapter 2 in [Pozrikidis1992]_). +class StokesletComponentKernel(StokesComponentKernelBase): + r"""The velocity kernel for the Stokes equations (see e.g. Chapter 2 in + [Pozrikidis1992]_). .. math:: @@ -1065,7 +1216,6 @@ class StokesletComponentKernel(ExpressionKernel): .. autoattribute:: icomp .. autoattribute:: jcomp - .. autoattribute:: viscosity_mu_name """ mapper_method: ClassVar[str] = "map_stokeslet_kernel" @@ -1075,11 +1225,6 @@ class StokesletComponentKernel(ExpressionKernel): jcomp: int """Component index for the kernel.""" - viscosity_mu_name: str - r"""The argument name to use for the dynamic viscosity :math:`\mu` when - generating functions to evaluate this kernel. - """ - def __init__(self, dim: int, icomp: int, @@ -1104,10 +1249,10 @@ def __init__(self, else: raise NotImplementedError(f"unsupported dimension: '{dim}'") - super().__init__(dim, expression=expr, global_scaling_const=scaling) + super().__init__(dim, expression=expr, global_scaling_const=scaling, + viscosity_mu_name=viscosity_mu_name) object.__setattr__(self, "icomp", icomp) object.__setattr__(self, "jcomp", jcomp) - object.__setattr__(self, "viscosity_mu_name", viscosity_mu_name) @override def __str__(self) -> str: @@ -1121,20 +1266,6 @@ def __reduce__(self): type(self), (self.dim, self.icomp, self.jcomp, self.viscosity_mu_name)) - @property - @override - def is_complex_valued(self) -> bool: - return False - - @memoize_method - @override - def get_args(self) -> Sequence[KernelArgument]: - return [ - KernelArgument( - loopy_arg=lp.ValueArg(self.viscosity_mu_name, np.float64), - ) - ] - @override @memoize_method def get_pde_system_kernel(self) -> tuple[SystemKernel, tuple[int, ...]]: @@ -1143,17 +1274,11 @@ def get_pde_system_kernel(self) -> tuple[SystemKernel, tuple[int, ...]]: viscosity_mu_name=self.viscosity_mu_name, ), (self.icomp, self.jcomp) - @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) - return laplacian(laplacian(w)) - @dataclass(frozen=True, repr=False) -class StressletComponentKernel(ExpressionKernel): - r"""A kernel for the Stokes equations (see e.g. Chapter 2 in [Pozrikidis1992]_). +class StressletComponentKernel(StokesComponentKernelBase): + r"""The stress kernel for the Stokes equations (see e.g. Chapter 2 in + [Pozrikidis1992]_). .. math:: @@ -1168,7 +1293,6 @@ class StressletComponentKernel(ExpressionKernel): .. autoattribute:: icomp .. autoattribute:: jcomp .. autoattribute:: kcomp - .. autoattribute:: viscosity_mu_name """ mapper_method: ClassVar[str] = "map_stresslet_kernel" @@ -1179,11 +1303,6 @@ class StressletComponentKernel(ExpressionKernel): kcomp: int """Component index for the kernel.""" - viscosity_mu_name: str - r"""The argument name to use for the dynamic viscosity :math:`\mu` when - generating functions to evaluate this kernel. - """ - def __init__(self, dim: int, icomp: int, @@ -1208,12 +1327,12 @@ def __init__(self, else: raise NotImplementedError(f"unsupported dimension: '{dim}'") - super().__init__(dim, expression=expr, global_scaling_const=scaling) + super().__init__(dim, expression=expr, global_scaling_const=scaling, + viscosity_mu_name=viscosity_mu_name) object.__setattr__(self, "icomp", icomp) object.__setattr__(self, "jcomp", jcomp) object.__setattr__(self, "kcomp", kcomp) - object.__setattr__(self, "viscosity_mu_name", viscosity_mu_name) @override def __reduce__(self) -> tuple[object, ...]: @@ -1222,19 +1341,6 @@ def __reduce__(self) -> tuple[object, ...]: (self.dim, self.icomp, self.jcomp, self.kcomp, self.viscosity_mu_name)) - @property - @override - def is_complex_valued(self) -> bool: - return False - - @memoize_method - @override - def get_args(self) -> Sequence[KernelArgument]: - # NOTE: this is not used, kept for symmetry with the StokesletComponentKernel - return [ - KernelArgument(loopy_arg=lp.ValueArg(self.viscosity_mu_name, np.float64)) - ] - @override def __str__(self) -> str: return ( @@ -1249,13 +1355,6 @@ def get_pde_system_kernel(self) -> tuple[SystemKernel, tuple[int, ...]]: viscosity_mu_name=self.viscosity_mu_name, ), (self.icomp, self.jcomp, self.kcomp) - @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) - return laplacian(laplacian(w)) - @dataclass(frozen=True, repr=False) class LineOfCompressionKernel(ExpressionKernel): @@ -1359,8 +1458,54 @@ def get_pde_as_diff_op(self) -> LinearPDESystemOperator: @dataclass(frozen=True, repr=False) -class BrinkmanletComponentKernel(ExpressionKernel): - r"""A kernel for the Brinkman equations. +class BrinkmanComponentKernelBase(ExpressionKernel): + """Base class for the Brinkman equations. + + .. autoattribute:: viscosity_mu_name + .. autoattribute:: darcy_impermeability_name + """ + + 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. + """ + + @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 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 BrinkmanletComponentKernel(BrinkmanComponentKernelBase): + r"""The velocity kernel for the Brinkman equations. .. math:: @@ -1373,6 +1518,9 @@ class BrinkmanletComponentKernel(ExpressionKernel): \end{cases} where :math:`P_j` is the pressure kernel. + + .. autoattribute:: icomp + .. autoattribute:: jcomp """ mapper_method: ClassVar[str] = "map_brinkmanlet_kernel" @@ -1382,15 +1530,6 @@ class BrinkmanletComponentKernel(ExpressionKernel): 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, @@ -1432,12 +1571,12 @@ def __init__(self, else: raise NotImplementedError(f"unsupported dimension: '{dim}'") - super().__init__(dim, expression=expr, global_scaling_const=scaling) + super().__init__(dim, expression=expr, global_scaling_const=scaling, + viscosity_mu_name=viscosity_mu_name, + darcy_impermeability_name=darcy_impermeability_name) 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, ...]: @@ -1447,31 +1586,12 @@ def __reduce__(self) -> tuple[object, ...]: 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 @memoize_method def get_pde_system_kernel(self) -> tuple[SystemKernel, tuple[int, ...]]: @@ -1481,18 +1601,9 @@ def get_pde_system_kernel(self) -> tuple[SystemKernel, tuple[int, ...]]: darcy_impermeability_name=self.darcy_impermeability_name, ), (self.icomp, self.jcomp) - @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 BrinkmanStressComponentKernel(ExpressionKernel): +class BrinkmanStressComponentKernel(BrinkmanComponentKernelBase): r"""A kernel for the Brinkman equations. .. math:: @@ -1504,6 +1615,10 @@ class BrinkmanStressComponentKernel(ExpressionKernel): where the two-index :math:`K_{ij}` is the :class:`~sumpy.kernel.BrinkmanletComponentKernel` and :math:`P_j` is the pressure kernel. + + .. autoattribute:: icomp + .. autoattribute:: jcomp + .. autoattribute:: kcomp """ mapper_method: ClassVar[str] = "map_brinkman_stress_kernel" @@ -1515,17 +1630,6 @@ class BrinkmanStressComponentKernel(ExpressionKernel): 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, @@ -1579,13 +1683,13 @@ def __init__(self, else: raise NotImplementedError(f"unsupported dimension: '{dim}'") - super().__init__(dim, expression=expr, global_scaling_const=scaling) + super().__init__(dim, expression=expr, global_scaling_const=scaling, + viscosity_mu_name=viscosity_mu_name, + darcy_impermeability_name=darcy_impermeability_name) 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, ...]: @@ -1595,31 +1699,12 @@ def __reduce__(self) -> tuple[object, ...]: 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 @memoize_method def get_pde_system_kernel(self) -> tuple[SystemKernel, tuple[int, ...]]: @@ -1629,14 +1714,6 @@ def get_pde_system_kernel(self) -> tuple[SystemKernel, tuple[int, ...]]: darcy_impermeability_name=self.darcy_impermeability_name, ), (self.icomp, self.jcomp, self.kcomp) - @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 HeatKernel(ExpressionKernel): @@ -1710,7 +1787,7 @@ def get_pde_as_diff_op(self) -> LinearPDESystemOperator: @dataclass(frozen=True, repr=False) class ElasticitySystemKernel(SystemKernel): - r"""A kernel for the linear elasticity (Navier-Cauchy) equations + r"""The displacement kernel for the linear elasticity (Navier-Cauchy) equations (see e.g. Section 2.2 in [Hsiao2008]_). This kernel uses :class:`ElasticityComponentKernel` for its components. @@ -1775,6 +1852,73 @@ def get_pde_as_diff_op(self) -> LinearPDESystemOperator: return mu * laplacian(u) + mu / (1 - 2 * nu) * gradient(divergence(u)) +@dataclass(frozen=True, repr=False) +class ElasticityStressSystemKernel(SystemKernel): + r"""The stress kernel for the linear elasticity (Navier-Cauchy) equations + (see e.g. Section 2.2 in [Hsiao2008]_). + + This kernel uses :class:`ElasticityStressComponentKernel` for its components. + + .. autoattribute:: viscosity_mu_name + .. autoattribute:: poisson_ratio_name + """ + + mapper_method: ClassVar[str] = "map_elasticity_stress_system_kernel" + + viscosity_mu_name: str = "mu" + r"""The argument name to use for the dynamic viscosity :math:`\mu` when + generating functions to evaluate this kernel. + """ + poisson_ratio_name: str = "nu" + r"""The argument name to use for Poisson's ratio :math:`\nu` when generating + functions to evaluate this kernel. + """ + + @override + def __str__(self) -> str: + return ( + f"ElasticityStressKnl{self.dim}D" + f"({self.viscosity_mu_name}, {self.poisson_ratio_name})") + + @override + def __reduce__(self): + return ( + type(self), + (self.dim, self.viscosity_mu_name, self.poisson_ratio_name)) + + @property + @override + def shape(self) -> tuple[int, ...]: + return (self.dim, self.dim, self.dim) + + @override + @keyed_memoize_method(key=lambda i, j, k: ((min(i, j), max(i, j), k))) + def get_scalar_component(self, i: int, j: int, k: int, /) -> ScalarKernel: + if i > j: + i, j = j, i + + return ElasticityStressComponentKernel( + self.dim, i, j, k, + viscosity_mu_name=self.viscosity_mu_name, + poisson_ratio_name=self.poisson_ratio_name, + ) + + @override + def get_pde_as_diff_op(self) -> LinearPDESystemOperator: + from sumpy.expansion.diff_op import ( + divergence, + gradient, + laplacian, + make_identity_diff_op, + ) + + mu = sym.Symbol(self.viscosity_mu_name) + nu = sym.Symbol(self.poisson_ratio_name) + u = make_identity_diff_op(self.dim, self.dim) + + return mu * laplacian(u) + mu / (1 - 2 * nu) * gradient(divergence(u)) + + @dataclass(frozen=True, repr=False) class StokesletSystemKernel(SystemKernel): r"""A kernel for the Stokes equations (see e.g. Chapter 2 in [Pozrikidis1992]_). @@ -2504,6 +2648,7 @@ def map_expression_kernel(self, kernel: ExpressionKernel) -> ScalarKernel: map_helmholtz_kernel: Callable[[Self, HelmholtzKernel], ScalarKernel] = map_expression_kernel # noqa: E501 map_yukawa_kernel: Callable[[Self, YukawaKernel], ScalarKernel] = map_expression_kernel # noqa: E501 map_elasticity_kernel: Callable[[Self, ElasticityComponentKernel], ScalarKernel] = map_expression_kernel # noqa: E501 + map_elasticity_stress_kernel: Callable[[Self, ElasticityStressComponentKernel], ScalarKernel] = map_expression_kernel # noqa: E501 map_line_of_compression_kernel: Callable[[Self, LineOfCompressionKernel], ScalarKernel] = map_expression_kernel # noqa: E501 map_stokeslet_kernel: Callable[[Self, StokesletComponentKernel], ScalarKernel] = map_expression_kernel # noqa: E501 map_stresslet_kernel: Callable[[Self, StressletComponentKernel], ScalarKernel] = map_expression_kernel # noqa: E501 @@ -2587,6 +2732,8 @@ def map_expression_kernel(self, kernel: ExpressionKernel) -> int: Callable[[Self, YukawaKernel], int] = map_expression_kernel map_elasticity_kernel: \ Callable[[Self, ElasticityComponentKernel], int] = map_expression_kernel + map_elasticity_stress_kernel: \ + Callable[[Self, ElasticityStressComponentKernel], int] = map_expression_kernel map_line_of_compression_kernel: \ Callable[[Self, LineOfCompressionKernel], int] = map_expression_kernel map_stokeslet_kernel: \ diff --git a/sumpy/test/test_misc.py b/sumpy/test/test_misc.py index c61e0149..27570571 100644 --- a/sumpy/test/test_misc.py +++ b/sumpy/test/test_misc.py @@ -63,6 +63,8 @@ BrinkmanStressComponentKernel, BrinkmanStressSystemKernel, ElasticityComponentKernel, + ElasticityStressComponentKernel, + ElasticityStressSystemKernel, ElasticitySystemKernel, ExpressionKernel, HeatKernel, @@ -143,6 +145,10 @@ def nderivs(self): KernelInfo(ElasticityComponentKernel(2, 0, 0), mu=5, nu=0.2), KernelInfo(ElasticityComponentKernel(3, 0, 1), mu=5, nu=0.2), KernelInfo(ElasticityComponentKernel(3, 0, 0), mu=5, nu=0.2), + KernelInfo(ElasticityStressComponentKernel(2, 0, 1, 0), mu=5, nu=0.2), + KernelInfo(ElasticityStressComponentKernel(2, 0, 0, 1), mu=5, nu=0.2), + KernelInfo(ElasticityStressComponentKernel(3, 0, 1, 0), mu=5, nu=0.2), + KernelInfo(ElasticityStressComponentKernel(3, 0, 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(BrinkmanletComponentKernel(2, 0, 1), mu=5, k=3), @@ -801,6 +807,8 @@ def test_get_storage_index(order, knl, compressed): def _get_scalar_cls(knl: SystemKernel) -> type[ScalarKernel]: if isinstance(knl, ElasticitySystemKernel): return ElasticityComponentKernel + if isinstance(knl, ElasticityStressSystemKernel): + return ElasticityStressComponentKernel elif isinstance(knl, StokesletSystemKernel): return StokesletComponentKernel elif isinstance(knl, StressletSystemKernel): @@ -815,6 +823,7 @@ def _get_scalar_cls(knl: SystemKernel) -> type[ScalarKernel]: @pytest.mark.parametrize("dim", [2, 3]) @pytest.mark.parametrize("cls", [ + ElasticityStressSystemKernel, ElasticitySystemKernel, StokesletSystemKernel, StressletSystemKernel, @@ -861,24 +870,31 @@ def test_system_kernel_components(dim: int, cls: type[SystemKernel]) -> None: knl_ij = knl[i, j] expected_ij = tuple(sorted([i, j])) - assert knl_ij == knl[j, i] + assert knl_ij is knl[j, i] assert (knl_ij.icomp, knl_ij.jcomp) == expected_ij elif isinstance(knl, (StressletSystemKernel,)): for i, j, k in product(range(dim), repeat=3): knl_ijk = knl[i, j, k] expected_ij = tuple(sorted([i, j, k])) - assert knl_ijk == knl[expected_ij] + assert knl_ijk is knl[expected_ij] assert (knl_ijk.icomp, knl_ijk.jcomp, knl_ijk.kcomp) == expected_ij elif isinstance(knl, (BrinkmanStressSystemKernel,)): for i, j, k in product(range(dim), repeat=3): knl_ijk = knl[i, j, k] expected_ij = min(i, k), j, max(i, k) - assert knl_ijk == knl[expected_ij] + assert knl_ijk is knl[expected_ij] + assert (knl_ijk.icomp, knl_ijk.jcomp, knl_ijk.kcomp) == expected_ij + elif isinstance(knl, (ElasticityStressSystemKernel,)): + for i, j, k in product(range(dim), repeat=3): + knl_ijk = knl[i, j, k] + expected_ij = min(i, j), max(i, j), k + + assert knl_ijk is knl[expected_ij] assert (knl_ijk.icomp, knl_ijk.jcomp, knl_ijk.kcomp) == expected_ij else: - raise AssertionError + raise AssertionError(knl) # }}}