diff --git a/pytential/symbolic/elasticity.py b/pytential/symbolic/elasticity.py new file mode 100644 index 000000000..64e6e766b --- /dev/null +++ b/pytential/symbolic/elasticity.py @@ -0,0 +1,166 @@ +__copyright__ = "Copyright (C) 2021 Isuru Fernando" + +__license__ = """ +Permission is hereby granted, free of charge, to any person obtaining a copy +of this software and associated documentation files (the "Software"), to deal +in the Software without restriction, including without limitation the rights +to use, copy, modify, merge, publish, distribute, sublicense, and/or sell +copies of the Software, and to permit persons to whom the Software is +furnished to do so, subject to the following conditions: + +The above copyright notice and this permission notice shall be included in +all copies or substantial portions of the Software. + +THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR +IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, +FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE +AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER +LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, +OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN +THE SOFTWARE. +""" + +import numpy as np + +from pytential import sym +from sumpy.kernel import (AxisTargetDerivative, AxisSourceDerivative, + TargetPointMultiplier, LaplaceKernel) +from pytential.symbolic.stokes import (StressletWrapperBase, StokesletWrapperBase, + _MU_SYM_DEFAULT) + + +class StressletWrapperYoshida(StressletWrapperBase): + """Stresslet Wrapper using Yoshida et al's method [1] which uses Laplace + derivatives. + + [1] Yoshida, K. I., Nishimura, N., & Kobayashi, S. (2001). Application of + fast multipole Galerkin boundary integral equation method to elastostatic + crack problems in 3D. + International Journal for Numerical Methods in Engineering, 50(3), 525-547. + """ + + def __init__(self, dim=None, mu_sym=_MU_SYM_DEFAULT, nu_sym=0.5): + self.dim = dim + if dim != 3: + raise ValueError("unsupported dimension given to " + "StressletWrapperYoshida") + self.kernel = LaplaceKernel(dim=3) + self.mu = mu_sym + self.nu = nu_sym + + def apply(self, density_vec_sym, dir_vec_sym, qbx_forced_limit, + extra_deriv_dirs=()): + return self.apply_stokeslet_and_stresslet([0]*self.dim, + density_vec_sym, dir_vec_sym, qbx_forced_limit, 0, 1, + extra_deriv_dirs) + + def apply_stokeslet_and_stresslet(self, stokeslet_density_vec_sym, + stresslet_density_vec_sym, dir_vec_sym, + qbx_forced_limit, stokeslet_weight, stresslet_weight, + extra_deriv_dirs=()): + + mu = self.mu + nu = self.nu + lam = 2*nu*mu/(1-2*nu) + stokeslet_weight *= -1 + + def C(i, j, k, l): # noqa: E741 + res = 0 + if i == j and k == l: + res += lam + if i == k and j == l: + res += mu + if i == l and j == k: + res += mu + return res * stresslet_weight + + def add_extra_deriv_dirs(target_kernel): + for deriv_dir in extra_deriv_dirs: + target_kernel = AxisTargetDerivative(deriv_dir, target_kernel) + return target_kernel + + def P(i, j, int_g): + int_g = int_g.copy(target_kernel=add_extra_deriv_dirs( + int_g.target_kernel)) + res = -int_g.copy(target_kernel=TargetPointMultiplier(j, + AxisTargetDerivative(i, int_g.target_kernel))) + if i == j: + res += (3 - 4*nu)*int_g + return res / (4*mu*(1 - nu)) + + def Q(i, int_g): + res = int_g.copy(target_kernel=add_extra_deriv_dirs( + AxisTargetDerivative(i, int_g.target_kernel))) + return res / (4*mu*(1 - nu)) + + sym_expr = np.zeros((3,), dtype=object) + + kernel = self.kernel + source = [sym.NodeCoordinateComponent(d) for d in range(3)] + normal = dir_vec_sym + sigma = stresslet_density_vec_sym + + source_kernels = [None]*4 + for i in range(3): + source_kernels[i] = AxisSourceDerivative(i, kernel) + source_kernels[3] = kernel + + for i in range(3): + for k in range(3): + densities = [0]*4 + for l in range(3): # noqa: E741 + for j in range(3): + for m in range(3): + densities[l] += C(k, l, m, j)*normal[m]*sigma[j] + densities[3] += stokeslet_weight * stokeslet_density_vec_sym[k] + int_g = sym.IntG(target_kernel=kernel, + source_kernels=tuple(source_kernels), + densities=tuple(densities), + qbx_forced_limit=qbx_forced_limit) + sym_expr[i] += P(i, k, int_g) + + densities = [0]*4 + for k in range(3): + for m in range(3): + for j in range(3): + for l in range(3): # noqa: E741 + densities[l] += \ + C(k, l, m, j)*normal[m]*sigma[j]*source[k] + if k == l: + densities[3] += \ + C(k, l, m, j)*normal[m]*sigma[j] + densities[3] += stokeslet_weight * source[k] \ + * stokeslet_density_vec_sym[k] + int_g = sym.IntG(target_kernel=kernel, + source_kernels=tuple(source_kernels), + densities=tuple(densities), + qbx_forced_limit=qbx_forced_limit) + sym_expr[i] += Q(i, int_g) + + return sym_expr + + +class StokesletWrapperYoshida(StokesletWrapperBase): + """Stokeslet Wrapper using Yoshida et al's method [1] which uses Laplace + derivatives. + + [1] Yoshida, K. I., Nishimura, N., & Kobayashi, S. (2001). Application of + fast multipole Galerkin boundary integral equation method to elastostatic + crack problems in 3D. + International Journal for Numerical Methods in Engineering, 50(3), 525-547. + """ + + def __init__(self, dim=None, mu_sym=_MU_SYM_DEFAULT, nu_sym=0.5): + self.dim = dim + if dim != 3: + raise ValueError("unsupported dimension given to " + "StokesletWrapperYoshida") + self.kernel = LaplaceKernel(dim=3) + self.mu = mu_sym + self.nu = nu_sym + + def apply(self, density_vec_sym, qbx_forced_limit, extra_deriv_dirs=()): + stresslet = StressletWrapperYoshida(3, self.mu, self.nu) + return stresslet.apply_stokeslet_and_stresslet(density_vec_sym, + [0]*self.dim, [0]*self.dim, qbx_forced_limit, 1, 0, + extra_deriv_dirs) diff --git a/pytential/symbolic/pde/system_utils.py b/pytential/symbolic/pde/system_utils.py new file mode 100644 index 000000000..061c3e2e8 --- /dev/null +++ b/pytential/symbolic/pde/system_utils.py @@ -0,0 +1,466 @@ +__copyright__ = "Copyright (C) 2020 Isuru Fernando" + +__license__ = """ +Permission is hereby granted, free of charge, to any person obtaining a copy +of this software and associated documentation files (the "Software"), to deal +in the Software without restriction, including without limitation the rights +to use, copy, modify, merge, publish, distribute, sublicense, and/or sell +copies of the Software, and to permit persons to whom the Software is +furnished to do so, subject to the following conditions: + +The above copyright notice and this permission notice shall be included in +all copies or substantial portions of the Software. + +THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR +IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, +FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE +AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER +LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, +OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN +THE SOFTWARE. +""" + +import numpy as np + +from sumpy.symbolic import make_sym_vector, SympyToPymbolicMapper +import sumpy.symbolic as sym +from sumpy.kernel import (AxisTargetDerivative, AxisSourceDerivative, + ExpressionKernel, KernelWrapper, TargetPointMultiplier) +from pytools import (memoize_on_first_arg, + generate_nonnegative_integer_tuples_summing_to_at_most + as gnitstam) + +from pytential.symbolic.primitives import (NodeCoordinateComponent, + hashable_kernel_args) +from pytential.symbolic.mappers import IdentityMapper +from pytential.utils import chop, lu_solve_with_expand +import pytential + +import logging +logger = logging.getLogger(__name__) + +__all__ = ( + "rewrite_using_base_kernel", + "get_deriv_relation", + ) + +__doc__ = """ +.. autofunction:: rewrite_using_base_kernel +.. autofunction:: get_deriv_relation +""" + + +# {{{ rewrite_using_base_kernel + +_NO_ARG_SENTINEL = object() + + +def rewrite_using_base_kernel(exprs, base_kernel=_NO_ARG_SENTINEL): + """Rewrites an expression with :class:`~pytential.symbolic.primitives.IntG` + objects using *base_kernel*. + + For example, if *base_kernel* is the Biharmonic kernel, and a Laplace kernel + is encountered, this will (forcibly) rewrite the Laplace kernel in terms of + derivatives of the Biharmonic kernel. + + The routine will fail if this process cannot be completed. + """ + if base_kernel is None: + return list(exprs) + mapper = RewriteUsingBaseKernelMapper(base_kernel) + return [mapper(expr) for expr in exprs] + + +class RewriteUsingBaseKernelMapper(IdentityMapper): + """Rewrites IntGs using the base kernel. First this method replaces + IntGs with :class:`sumpy.kernel.AxisTargetDerivative` to IntGs + :class:`sumpy.kernel.AxisSourceDerivative` and then replaces + IntGs with :class:`sumpy.kernel.TargetPointMultiplier` to IntGs + without them using :class:`sumpy.kernel.ExpressionKernel` + and then finally converts them to the base kernel by finding + a relationship between the derivatives. + """ + def __init__(self, base_kernel): + self.base_kernel = base_kernel + + def map_int_g(self, expr): + # First convert IntGs with target derivatives to source derivatives + expr = convert_target_deriv_to_source(expr) + # Convert IntGs with TargetMultiplier to a sum of IntGs without + # TargetMultipliers + new_int_gs = convert_target_multiplier_to_source(expr) + # Convert IntGs with different kernels to expressions containing + # IntGs with base_kernel or its derivatives + return sum(convert_int_g_to_base(new_int_g, + self.base_kernel) for new_int_g in new_int_gs) + + +def convert_target_deriv_to_source(int_g): + """Converts AxisTargetDerivatives to AxisSourceDerivative instances + from an IntG. If there are outer TargetPointMultiplier transformations + they are preserved. + """ + knl = int_g.target_kernel + source_kernels = list(int_g.source_kernels) + coeff = 1 + multipliers = [] + while isinstance(knl, TargetPointMultiplier): + multipliers.append(knl.axis) + knl = knl.inner_kernel + + while isinstance(knl, AxisTargetDerivative): + coeff *= -1 + source_kernels = [AxisSourceDerivative(knl.axis, source_knl) for + source_knl in source_kernels] + knl = knl.inner_kernel + + # TargetPointMultiplier has to be the outermost kernel + # If it is the inner kernel, return early + if isinstance(knl, TargetPointMultiplier): + return int_g + + for axis in reversed(multipliers): + knl = TargetPointMultiplier(axis, knl) + + new_densities = tuple(density*coeff for density in int_g.densities) + return int_g.copy(target_kernel=knl, + densities=new_densities, + source_kernels=tuple(source_kernels)) + + +def _get_kernel_expression(expr, kernel_arguments): + from pymbolic.mapper.substitutor import substitute + from sumpy.symbolic import PymbolicToSympyMapperWithSymbols + + pymbolic_expr = substitute(expr, kernel_arguments) + + res = PymbolicToSympyMapperWithSymbols()(pymbolic_expr) + return res + + +def _monom_to_expr(monom, variables): + prod = 1 + for i, nrepeats in enumerate(monom): + for _ in range(nrepeats): + prod *= variables[i] + return prod + + +def convert_target_multiplier_to_source(int_g): + """Convert an IntG with TargetMultiplier to a sum of IntGs without + TargetMultiplier and only source dependent transformations + """ + import sympy + import sumpy.symbolic as sym + from sumpy.symbolic import SympyToPymbolicMapper + conv = SympyToPymbolicMapper() + + knl = int_g.target_kernel + # we use a symbol for d = (x - y) + ds = sympy.symbols(f"d0:{knl.dim}") + sources = sympy.symbols(f"y0:{knl.dim}") + # instead of just x, we use x = (d + y) + targets = [d + source for d, source in zip(ds, sources)] + orig_expr = sympy.Function("f")(*ds) # pylint: disable=not-callable + expr = orig_expr + found = False + while isinstance(knl, KernelWrapper): + if isinstance(knl, TargetPointMultiplier): + expr = targets[knl.axis] * expr + found = True + elif isinstance(knl, AxisTargetDerivative): + # sympy can't differentiate w.r.t target because + # it's not a symbol, but d/d(x) = d/d(d) + expr = expr.diff(ds[knl.axis]) + else: + return [int_g] + knl = knl.inner_kernel + + if not found: + return [int_g] + + sources_pymbolic = [NodeCoordinateComponent(i) for i in range(knl.dim)] + expr = expr.expand() + # Now the expr is an Add and looks like + # u''(d, s)*d[0] + u(d, s) + assert isinstance(expr, sympy.Add) + result = [] + + for arg in expr.args: + deriv_terms = arg.atoms(sympy.Derivative) + if len(deriv_terms) == 1: + deriv_term = deriv_terms.pop() + rest_terms = sympy.Poly(arg.xreplace({deriv_term: 1}), *ds, *sources) + derivatives = deriv_term.args[1:] + elif len(deriv_terms) == 0: + rest_terms = sympy.Poly(arg.xreplace({orig_expr: 1}), *ds, *sources) + derivatives = [(d, 0) for d in ds] + else: + raise AssertionError("impossible condition") + assert len(rest_terms.terms()) == 1 + monom, coeff = rest_terms.terms()[0] + expr_multiplier = _monom_to_expr(monom[:len(ds)], ds) + density_multiplier = _monom_to_expr(monom[len(ds):], sources_pymbolic) \ + * conv(coeff) + + new_int_gs = _multiply_int_g(int_g, sym.sympify(expr_multiplier), + density_multiplier) + for new_int_g in new_int_gs: + knl = new_int_g.target_kernel + for axis_var, nrepeats in derivatives: + axis = ds.index(axis_var) + for _ in range(nrepeats): + knl = AxisTargetDerivative(axis, knl) + result.append(new_int_g.copy(target_kernel=knl)) + return result + + +def _multiply_int_g(int_g, expr_multiplier, density_multiplier): + """Multiply the exprssion in IntG with the *expr_multiplier* + which is a symbolic expression and multiply the densities + with *density_multiplier* which is a pymbolic expression. + """ + from sumpy.symbolic import SympyToPymbolicMapper + result = [] + + base_kernel = int_g.target_kernel.get_base_kernel() + sym_d = make_sym_vector("d", base_kernel.dim) + base_kernel_expr = _get_kernel_expression(base_kernel.expression, + int_g.kernel_arguments) + conv = SympyToPymbolicMapper() + + for knl, density in zip(int_g.source_kernels, int_g.densities): + if expr_multiplier == 1: + new_knl = knl.get_base_kernel() + else: + new_expr = conv(knl.postprocess_at_source(base_kernel_expr, sym_d) + * expr_multiplier) + new_knl = ExpressionKernel(knl.dim, new_expr, + knl.get_base_kernel().global_scaling_const, + knl.is_complex_valued) + result.append(int_g.copy(target_kernel=new_knl, + densities=(density*density_multiplier,), + source_kernels=(new_knl,))) + return result + + +def convert_int_g_to_base(int_g, base_kernel): + result = 0 + for knl, density in zip(int_g.source_kernels, int_g.densities): + result += _convert_int_g_to_base( + int_g.copy(source_kernels=(knl,), densities=(density,)), + base_kernel) + return result + + +def _convert_int_g_to_base(int_g, base_kernel): + target_kernel = int_g.target_kernel.replace_base_kernel(base_kernel) + dim = target_kernel.dim + + result = 0 + for density, source_kernel in zip(int_g.densities, int_g.source_kernels): + deriv_relation = get_deriv_relation_kernel(source_kernel.get_base_kernel(), + base_kernel, hashable_kernel_arguments=( + hashable_kernel_args(int_g.kernel_arguments))) + + const = deriv_relation[0] + # NOTE: we set a dofdesc here to force the evaluation of this integral + # on the source instead of the target when using automatic tagging + # see :meth:`pytential.symbolic.mappers.LocationTagger._default_dofdesc` + dd = pytential.sym.DOFDescriptor(None, + discr_stage=pytential.sym.QBX_SOURCE_STAGE1) + const *= pytential.sym.integral(dim, dim-1, density, dofdesc=dd) + + if const != 0 and target_kernel != target_kernel.get_base_kernel(): + # There might be some TargetPointMultipliers hanging around. + # FIXME: handle them instead of bailing out + return [int_g] + + if source_kernel != source_kernel.get_base_kernel(): + # We assume that any source transformation is a derivative + # and the constant when applied becomes zero. + const = 0 + + result += const + + new_kernel_args = filter_kernel_arguments([base_kernel], + int_g.kernel_arguments) + + for mi, c in deriv_relation[1]: + knl = source_kernel.replace_base_kernel(base_kernel) + for d, val in enumerate(mi): + for _ in range(val): + knl = AxisSourceDerivative(d, knl) + c *= -1 + result += int_g.copy(source_kernels=(knl,), target_kernel=target_kernel, + densities=(density * c,), kernel_arguments=new_kernel_args) + return result + + +def get_deriv_relation(kernels, base_kernel, tol=1e-10, order=None, + kernel_arguments=None): + res = [] + for knl in kernels: + res.append(get_deriv_relation_kernel(knl, base_kernel, tol, order, + hashable_kernel_arguments=hashable_kernel_args(kernel_arguments))) + return res + + +@memoize_on_first_arg +def get_deriv_relation_kernel(kernel, base_kernel, tol=1e-10, order=None, + hashable_kernel_arguments=None): + kernel_arguments = dict(hashable_kernel_arguments) + (L, U, perm), rand, mis = _get_base_kernel_matrix(base_kernel, order=order) + dim = base_kernel.dim + sym_vec = make_sym_vector("d", dim) + sympy_conv = SympyToPymbolicMapper() + + expr = _get_kernel_expression(kernel.expression, kernel_arguments) + vec = [] + for i in range(len(mis)): + vec.append(evalf(expr.xreplace(dict((k, v) for + k, v in zip(sym_vec, rand[:, i]))))) + vec = sym.Matrix(vec) + result = [] + const = 0 + logger.debug("%s = ", kernel) + + sol = lu_solve_with_expand(L, U, perm, vec) + for i, coeff in enumerate(sol): + coeff = chop(coeff, tol) + if coeff == 0: + continue + if mis[i] != (-1, -1, -1): + coeff *= _get_kernel_expression(kernel.global_scaling_const, + kernel_arguments) + coeff /= _get_kernel_expression(base_kernel.global_scaling_const, + kernel_arguments) + result.append((mis[i], sympy_conv(coeff))) + logger.debug(" + %s.diff(%s)*%s", base_kernel, mis[i], coeff) + else: + const = sympy_conv(coeff * _get_kernel_expression( + kernel.global_scaling_const, kernel_arguments)) + logger.debug(" + %s", const) + return (const, result) + + +@memoize_on_first_arg +def _get_base_kernel_matrix(base_kernel, order=None, retries=3, + kernel_arguments=None): + dim = base_kernel.dim + + pde = base_kernel.get_pde_as_diff_op() + if order is None: + order = pde.order + + if order > pde.order: + raise NotImplementedError(f"order ({order}) cannot be greater than the order" + f"of the PDE ({pde.order}) yet.") + + mis = sorted(gnitstam(order, dim), key=sum) + # (-1, -1, -1) represent a constant + mis.append((-1, -1, -1)) + + if order == pde.order: + pde_mis = [ident.mi for eq in pde.eqs for ident in eq.keys()] + pde_mis = [mi for mi in pde_mis if sum(mi) == order] + logger.debug(f"Removing {pde_mis[-1]} to avoid linear dependent mis") + mis.remove(pde_mis[-1]) + + rand = np.random.randint(1, 10**15, (dim, len(mis))) + rand = rand.astype(object) + for i in range(rand.shape[0]): + for j in range(rand.shape[1]): + rand[i, j] = sym.sympify(rand[i, j])/10**15 + sym_vec = make_sym_vector("d", dim) + + base_expr = _get_kernel_expression(base_kernel.expression, kernel_arguments) + + mat = [] + for rand_vec_idx in range(rand.shape[1]): + row = [] + for mi in mis[:-1]: + expr = base_expr + for var_idx, nderivs in enumerate(mi): + if nderivs == 0: + continue + expr = expr.diff(sym_vec[var_idx], nderivs) + replace_dict = dict( + (k, v) for k, v in zip(sym_vec, rand[:, rand_vec_idx]) + ) + eval_expr = evalf(expr.xreplace(replace_dict)) + row.append(eval_expr) + row.append(1) + mat.append(row) + + mat = sym.Matrix(mat) + failed = False + try: + L, U, perm = mat.LUdecomposition() + except RuntimeError: + # symengine throws an error when rank deficient + # and sympy returns U with last row zero + failed = True + + if not sym.USE_SYMENGINE and all(expr == 0 for expr in U[-1, :]): + failed = True + + if failed: + if retries == 0: + raise RuntimeError("Failed to find a base kernel") + return _get_base_kernel_matrix( + base_kernel=base_kernel, + order=order, + retries=retries-1, + ) + + return (L, U, perm), rand, mis + + +def evalf(expr, prec=100): + """evaluate an expression numerically using ``prec`` + number of bits. + """ + from sumpy.symbolic import USE_SYMENGINE + if USE_SYMENGINE: + return expr.n(prec=prec) + else: + import sympy + dps = int(sympy.log(2**prec, 10)) + return expr.n(n=dps) + + +def filter_kernel_arguments(knls, kernel_arguments): + """From a dictionary of kernel arguments, filter out arguments + that are not needed for the kernels given as a list and return a new + dictionary. + """ + kernel_arg_names = set() + + for kernel in knls: + for karg in (kernel.get_args() + kernel.get_source_args()): + kernel_arg_names.add(karg.loopy_arg.name) + + return {k: v for (k, v) in kernel_arguments.items() if k in kernel_arg_names} + +# }}} + + +if __name__ == "__main__": + logging.basicConfig(level=logging.DEBUG) + from sumpy.kernel import (StokesletKernel, BiharmonicKernel, # noqa:F401 + StressletKernel, ElasticityKernel, LaplaceKernel) + base_kernel = BiharmonicKernel(3) + #base_kernel = LaplaceKernel(3) + kernels = [StokesletKernel(3, 0, 2), StokesletKernel(3, 0, 0)] + kernels += [StressletKernel(3, 0, 0, 0), StressletKernel(3, 0, 0, 1), + StressletKernel(3, 0, 0, 2), StressletKernel(3, 0, 1, 2)] + + sym_d = make_sym_vector("d", base_kernel.dim) + sym_r = sym.sqrt(sum(a**2 for a in sym_d)) + conv = SympyToPymbolicMapper() + expression_knl = ExpressionKernel(3, conv(sym_d[0]*sym_d[1]/sym_r**3), 1, False) + expression_knl2 = ExpressionKernel(3, conv(1/sym_r + sym_d[0]*sym_d[0]/sym_r**3), + 1, False) + kernels = [expression_knl, expression_knl2] + get_deriv_relation(kernels, base_kernel, tol=1e-10, order=4) diff --git a/pytential/symbolic/stokes.py b/pytential/symbolic/stokes.py index cfe8bef3c..c325dd3a0 100644 --- a/pytential/symbolic/stokes.py +++ b/pytential/symbolic/stokes.py @@ -23,7 +23,12 @@ import numpy as np from pytential import sym -from sumpy.kernel import StokesletKernel, StressletKernel, LaplaceKernel +from pytential.symbolic.pde.system_utils import rewrite_using_base_kernel +from sumpy.kernel import (StressletKernel, LaplaceKernel, + ElasticityKernel, BiharmonicKernel, + AxisTargetDerivative, AxisSourceDerivative, TargetPointMultiplier) +from sumpy.symbolic import SpatialConstant +from abc import ABC __doc__ = """ .. autoclass:: StokesletWrapper @@ -35,9 +40,12 @@ """ -# {{{ StokesletWrapper +# {{{ StokesletWrapper/StressletWrapper ABCs -class StokesletWrapper: +_MU_SYM_DEFAULT = SpatialConstant("mu") + + +class StokesletWrapperBase(ABC): """Wrapper class for the :class:`~sumpy.kernel.StokesletKernel` kernel. This class is meant to shield the user from the messiness of writing @@ -59,43 +67,17 @@ class StokesletWrapper: :meth:`apply_stress` (applies symmetric viscous stress tensor in the requested direction). - .. attribute:: kernel_dict - - The dictionary allows us to exploit symmetry -- that - :math:`S_{01}` is identical to :math:`S_{10}` -- and avoid creating - multiple expansions for the same kernel in a different ordering. - - .. automethod:: __init__ .. automethod:: apply .. automethod:: apply_pressure .. automethod:: apply_derivative .. automethod:: apply_stress """ - - def __init__(self, dim=None): + def __init__(self, dim, mu_sym, nu_sym): self.dim = dim + self.mu = mu_sym + self.nu = nu_sym - if dim == 2: - self.kernel_dict = { - (2, 0): StokesletKernel(dim=2, icomp=0, jcomp=0), - (1, 1): StokesletKernel(dim=2, icomp=0, jcomp=1), - (0, 2): StokesletKernel(dim=2, icomp=1, jcomp=1) - } - - elif dim == 3: - self.kernel_dict = { - (2, 0, 0): StokesletKernel(dim=3, icomp=0, jcomp=0), - (1, 1, 0): StokesletKernel(dim=3, icomp=0, jcomp=1), - (1, 0, 1): StokesletKernel(dim=3, icomp=0, jcomp=2), - (0, 2, 0): StokesletKernel(dim=3, icomp=1, jcomp=1), - (0, 1, 1): StokesletKernel(dim=3, icomp=1, jcomp=2), - (0, 0, 2): StokesletKernel(dim=3, icomp=2, jcomp=2) - } - - else: - raise ValueError("unsupported dimension given to StokesletWrapper") - - def apply(self, density_vec_sym, mu_sym, qbx_forced_limit): + def apply(self, density_vec_sym, qbx_forced_limit, extra_deriv_dirs=()): """Symbolic expressions for integrating Stokeslet kernel. Returns an object array of symbolic expressions for the vector @@ -103,58 +85,30 @@ def apply(self, density_vec_sym, mu_sym, qbx_forced_limit): variable *density_vec_sym*. :arg density_vec_sym: a symbolic vector variable for the density vector. - :arg mu_sym: a symbolic variable for the viscosity. :arg qbx_forced_limit: the *qbx_forced_limit* argument to be passed on to :class:`~pytential.symbolic.primitives.IntG`. + :arg extra_deriv_dirs: adds target derivatives to all the integral + objects with the given derivative axis. """ + raise NotImplementedError - sym_expr = np.empty((self.dim,), dtype=object) - - for comp in range(self.dim): - - # Start variable count for kernel with 1 for the requested result - # component - base_count = np.zeros(self.dim, dtype=np.int32) - base_count[comp] += 1 - - for i in range(self.dim): - var_ctr = base_count.copy() - var_ctr[i] += 1 - ctr_key = tuple(var_ctr) - - if i < 1: - sym_expr[comp] = sym.int_g_vec( - self.kernel_dict[ctr_key], density_vec_sym[i], - qbx_forced_limit=qbx_forced_limit, mu=mu_sym) - - else: - sym_expr[comp] = sym_expr[comp] + sym.int_g_vec( - self.kernel_dict[ctr_key], density_vec_sym[i], - qbx_forced_limit=qbx_forced_limit, mu=mu_sym) - - return sym_expr - - def apply_pressure(self, density_vec_sym, mu_sym, qbx_forced_limit): + def apply_pressure(self, density_vec_sym, qbx_forced_limit): """Symbolic expression for pressure field associated with the Stokeslet.""" + # Pressure representation doesn't differ depending on the implementation + # and is implemented in base class here. from pytential.symbolic.mappers import DerivativeTaker kernel = LaplaceKernel(dim=self.dim) + sym_expr = 0 for i in range(self.dim): - - if i < 1: - sym_expr = DerivativeTaker(i).map_int_g( - sym.int_g_vec(kernel, density_vec_sym[i], - qbx_forced_limit=qbx_forced_limit)) - else: - sym_expr = sym_expr + (DerivativeTaker(i).map_int_g( - sym.int_g_vec(kernel, density_vec_sym[i], - qbx_forced_limit=qbx_forced_limit))) + sym_expr += (DerivativeTaker(i).map_int_g( + sym.S(kernel, density_vec_sym[i], + qbx_forced_limit=qbx_forced_limit))) return sym_expr - def apply_derivative(self, deriv_dir, density_vec_sym, - mu_sym, qbx_forced_limit): + def apply_derivative(self, deriv_dir, density_vec_sym, qbx_forced_limit): """Symbolic derivative of velocity from Stokeslet. Returns an object array of symbolic expressions for the vector @@ -163,46 +117,12 @@ def apply_derivative(self, deriv_dir, density_vec_sym, :arg deriv_dir: integer denoting the axis direction for the derivative. :arg density_vec_sym: a symbolic vector variable for the density vector. - :arg mu_sym: a symbolic variable for the viscosity. :arg qbx_forced_limit: the *qbx_forced_limit* argument to be passed on to :class:`~pytential.symbolic.primitives.IntG`. """ + return self.apply(density_vec_sym, qbx_forced_limit, [deriv_dir]) - from pytential.symbolic.mappers import DerivativeTaker - - sym_expr = np.empty((self.dim,), dtype=object) - - for comp in range(self.dim): - - # Start variable count for kernel with 1 for the requested result - # component - base_count = np.zeros(self.dim, dtype=np.int32) - base_count[comp] += 1 - - for i in range(self.dim): - var_ctr = base_count.copy() - var_ctr[i] += 1 - ctr_key = tuple(var_ctr) - - if i < 1: - sym_expr[comp] = DerivativeTaker(deriv_dir).map_int_g( - sym.int_g_vec(self.kernel_dict[ctr_key], - density_vec_sym[i], - qbx_forced_limit=qbx_forced_limit, - mu=mu_sym)) - - else: - sym_expr[comp] = sym_expr[comp] + DerivativeTaker( - deriv_dir).map_int_g( - sym.int_g_vec(self.kernel_dict[ctr_key], - density_vec_sym[i], - qbx_forced_limit=qbx_forced_limit, - mu=mu_sym)) - - return sym_expr - - def apply_stress(self, density_vec_sym, dir_vec_sym, - mu_sym, qbx_forced_limit): + def apply_stress(self, density_vec_sym, dir_vec_sym, qbx_forced_limit): r"""Symbolic expression for viscous stress applied to a direction. Returns a vector of symbolic expressions for the force resulting @@ -222,50 +142,13 @@ def apply_stress(self, density_vec_sym, dir_vec_sym, :arg density_vec_sym: a symbolic vector variable for the density vector. :arg dir_vec_sym: a symbolic vector for the application direction. - :arg mu_sym: a symbolic variable for the viscosity. :arg qbx_forced_limit: the *qbx_forced_limit* argument to be passed on to :class:`~pytential.symbolic.primitives.IntG`. """ + raise NotImplementedError - import itertools - - sym_expr = np.empty((self.dim,), dtype=object) - stresslet_obj = StressletWrapper(dim=self.dim) - - for comp in range(self.dim): - - # Start variable count for kernel with 1 for the requested result - # component - base_count = np.zeros(self.dim, dtype=np.int32) - base_count[comp] += 1 - - for i, j in itertools.product(range(self.dim), range(self.dim)): - var_ctr = base_count.copy() - var_ctr[i] += 1 - var_ctr[j] += 1 - ctr_key = tuple(var_ctr) - - if i + j < 1: - sym_expr[comp] = dir_vec_sym[i] * sym.int_g_vec( - stresslet_obj.kernel_dict[ctr_key], - density_vec_sym[j], - qbx_forced_limit=qbx_forced_limit, mu=mu_sym) - - else: - sym_expr[comp] = sym_expr[comp] + dir_vec_sym[i] * sym.int_g_vec( - stresslet_obj.kernel_dict[ctr_key], - density_vec_sym[j], - qbx_forced_limit=qbx_forced_limit, - mu=mu_sym) - - return sym_expr - -# }}} - - -# {{{ StressletWrapper -class StressletWrapper: +class StressletWrapperBase(ABC): """Wrapper class for the :class:`~sumpy.kernel.StressletKernel` kernel. This class is meant to shield the user from the messiness of writing @@ -286,48 +169,18 @@ class StressletWrapper: :meth:`apply_stress` (applies symmetric viscous stress tensor in the requested direction). - .. attribute:: kernel_dict - - The dictionary allows us to exploit symmetry -- that - :math:`T_{012}` is identical to :math:`T_{120}` -- and avoid creating - multiple expansions for the same kernel in a different ordering. - - .. automethod:: __init__ .. automethod:: apply .. automethod:: apply_pressure .. automethod:: apply_derivative .. automethod:: apply_stress """ - - def __init__(self, dim=None): + def __init__(self, dim, mu_sym, nu_sym): self.dim = dim + self.mu = mu_sym + self.nu = nu_sym - if dim == 2: - self.kernel_dict = { - (3, 0): StressletKernel(dim=2, icomp=0, jcomp=0, kcomp=0), - (2, 1): StressletKernel(dim=2, icomp=0, jcomp=0, kcomp=1), - (1, 2): StressletKernel(dim=2, icomp=0, jcomp=1, kcomp=1), - (0, 3): StressletKernel(dim=2, icomp=1, jcomp=1, kcomp=1) - } - - elif dim == 3: - self.kernel_dict = { - (3, 0, 0): StressletKernel(dim=3, icomp=0, jcomp=0, kcomp=0), - (2, 1, 0): StressletKernel(dim=3, icomp=0, jcomp=0, kcomp=1), - (2, 0, 1): StressletKernel(dim=3, icomp=0, jcomp=0, kcomp=2), - (1, 2, 0): StressletKernel(dim=3, icomp=0, jcomp=1, kcomp=1), - (1, 1, 1): StressletKernel(dim=3, icomp=0, jcomp=1, kcomp=2), - (1, 0, 2): StressletKernel(dim=3, icomp=0, jcomp=2, kcomp=2), - (0, 3, 0): StressletKernel(dim=3, icomp=1, jcomp=1, kcomp=1), - (0, 2, 1): StressletKernel(dim=3, icomp=1, jcomp=1, kcomp=2), - (0, 1, 2): StressletKernel(dim=3, icomp=1, jcomp=2, kcomp=2), - (0, 0, 3): StressletKernel(dim=3, icomp=2, jcomp=2, kcomp=2) - } - - else: - raise ValueError("unsupported dimension given to StressletWrapper") - - def apply(self, density_vec_sym, dir_vec_sym, mu_sym, qbx_forced_limit): + def apply(self, density_vec_sym, dir_vec_sym, qbx_forced_limit, + extra_deriv_dirs=()): """Symbolic expressions for integrating Stresslet kernel. Returns an object array of symbolic expressions for the vector @@ -336,124 +189,55 @@ def apply(self, density_vec_sym, dir_vec_sym, mu_sym, qbx_forced_limit): :arg density_vec_sym: a symbolic vector variable for the density vector. :arg dir_vec_sym: a symbolic vector variable for the direction vector. - :arg mu_sym: a symbolic variable for the viscosity. :arg qbx_forced_limit: the *qbx_forced_limit* argument to be passed on to :class:`~pytential.symbolic.primitives.IntG`. + :arg extra_deriv_dirs: adds target derivatives to all the integral + objects with the given derivative axis. """ + raise NotImplementedError - import itertools - - sym_expr = np.empty((self.dim,), dtype=object) - - for comp in range(self.dim): - - # Start variable count for kernel with 1 for the requested result - # component - base_count = np.zeros(self.dim, dtype=np.int32) - base_count[comp] += 1 - - for i, j in itertools.product(range(self.dim), range(self.dim)): - var_ctr = base_count.copy() - var_ctr[i] += 1 - var_ctr[j] += 1 - ctr_key = tuple(var_ctr) - - if i + j < 1: - sym_expr[comp] = sym.int_g_vec( - self.kernel_dict[ctr_key], - dir_vec_sym[i] * density_vec_sym[j], - qbx_forced_limit=qbx_forced_limit, mu=mu_sym) - - else: - sym_expr[comp] = sym_expr[comp] + sym.int_g_vec( - self.kernel_dict[ctr_key], - dir_vec_sym[i] * density_vec_sym[j], - qbx_forced_limit=qbx_forced_limit, - mu=mu_sym) - - return sym_expr - - def apply_pressure(self, density_vec_sym, dir_vec_sym, mu_sym, qbx_forced_limit): - """Symbolic expression for pressure field associated with the Stresslet.""" + def apply_pressure(self, density_vec_sym, dir_vec_sym, qbx_forced_limit): + """Symbolic expression for pressure field associated with the Stresslet. + """ + # Pressure representation doesn't differ depending on the implementation + # and is implemented in base class here. import itertools from pytential.symbolic.mappers import DerivativeTaker kernel = LaplaceKernel(dim=self.dim) - factor = (2. * mu_sym) + factor = (2. * self.mu) - for i, j in itertools.product(range(self.dim), range(self.dim)): + sym_expr = 0 - if i + j < 1: - sym_expr = factor * DerivativeTaker(i).map_int_g( - DerivativeTaker(j).map_int_g( - sym.int_g_vec(kernel, - density_vec_sym[i] * dir_vec_sym[j], - qbx_forced_limit=qbx_forced_limit))) - else: - sym_expr = sym_expr + ( - factor * DerivativeTaker(i).map_int_g( + for i, j in itertools.product(range(self.dim), range(self.dim)): + sym_expr += factor * DerivativeTaker(i).map_int_g( DerivativeTaker(j).map_int_g( sym.int_g_vec(kernel, density_vec_sym[i] * dir_vec_sym[j], - qbx_forced_limit=qbx_forced_limit)))) + qbx_forced_limit=qbx_forced_limit))) return sym_expr def apply_derivative(self, deriv_dir, density_vec_sym, dir_vec_sym, - mu_sym, qbx_forced_limit): - """Symbolic derivative of velocity from stresslet. + qbx_forced_limit): + """Symbolic derivative of velocity from Stokeslet. Returns an object array of symbolic expressions for the vector resulting from integrating the *deriv_dir* target derivative of the - dyadic Stresslet kernel with variable *density_vec_sym* and source - direction vectors *dir_vec_sym*. + dyadic Stokeslet kernel with variable *density_vec_sym*. :arg deriv_dir: integer denoting the axis direction for the derivative. :arg density_vec_sym: a symbolic vector variable for the density vector. :arg dir_vec_sym: a symbolic vector variable for the normal direction. - :arg mu_sym: a symbolic variable for the viscosity. :arg qbx_forced_limit: the *qbx_forced_limit* argument to be passed on to :class:`~pytential.symbolic.primitives.IntG`. """ - - import itertools - from pytential.symbolic.mappers import DerivativeTaker - - sym_expr = np.empty((self.dim,), dtype=object) - - for comp in range(self.dim): - - # Start variable count for kernel with 1 for the requested result - # component - base_count = np.zeros(self.dim, dtype=np.int32) - base_count[comp] += 1 - - for i, j in itertools.product(range(self.dim), range(self.dim)): - var_ctr = base_count.copy() - var_ctr[i] += 1 - var_ctr[j] += 1 - ctr_key = tuple(var_ctr) - - if i + j < 1: - sym_expr[comp] = DerivativeTaker(deriv_dir).map_int_g( - sym.int_g_vec(self.kernel_dict[ctr_key], - dir_vec_sym[i] * density_vec_sym[j], - qbx_forced_limit=qbx_forced_limit, - mu=mu_sym)) - - else: - sym_expr[comp] = sym_expr[comp] + DerivativeTaker( - deriv_dir).map_int_g( - sym.int_g_vec(self.kernel_dict[ctr_key], - dir_vec_sym[i] * density_vec_sym[j], - qbx_forced_limit=qbx_forced_limit, - mu=mu_sym)) - - return sym_expr + return self.apply(density_vec_sym, dir_vec_sym, qbx_forced_limit, + [deriv_dir]) def apply_stress(self, density_vec_sym, normal_vec_sym, dir_vec_sym, - mu_sym, qbx_forced_limit): + qbx_forced_limit): r"""Symbolic expression for viscous stress applied to a direction. Returns a vector of symbolic expressions for the force resulting @@ -469,10 +253,234 @@ def apply_stress(self, density_vec_sym, normal_vec_sym, dir_vec_sym, :arg normal_vec_sym: a symbolic vector variable for the normal vectors (outward facing normals at source locations). :arg dir_vec_sym: a symbolic vector for the application direction. - :arg mu_sym: a symbolic variable for the viscosity. :arg qbx_forced_limit: the *qbx_forced_limit* argument to be passed on to :class:`~pytential.symbolic.primitives.IntG`. """ + raise NotImplementedError + +# }}} + + +# {{{ StokesletWrapper Naive and Biharmonic impl + +def _create_int_g(knl, deriv_dirs, density, **kwargs): + for deriv_dir in deriv_dirs: + knl = AxisTargetDerivative(deriv_dir, knl) + + kernel_arg_names = set(karg.loopy_arg.name + for karg in (knl.get_args() + knl.get_source_args())) + + # When the kernel is Laplace, mu and nu are not kernel arguments + # Also when nu==0.5, it's not a kernel argument to StokesletKernel + for var_name in ["mu", "nu"]: + if var_name not in kernel_arg_names: + kwargs.pop(var_name) + + res = sym.int_g_vec(knl, density, **kwargs) + return res + + +class _StokesletWrapperNaiveOrBiharmonic(StokesletWrapperBase): + def __init__(self, dim=None, mu_sym=_MU_SYM_DEFAULT, nu_sym=0.5, + method="biharmonic"): + super().__init__(dim, mu_sym, nu_sym) + if not (dim == 3 or dim == 2): + raise ValueError("unsupported dimension given to StokesletWrapper") + + self.method = method + if method == "biharmonic": + self.base_kernel = BiharmonicKernel(dim) + elif method == "naive": + self.base_kernel = None + else: + raise ValueError("method has to be one of biharmonic/naive") + + self.kernel_dict = {} + # The two cases of nu=0.5 and nu!=0.5 differ significantly and + # ElasticityKernel needs to know if nu=0.5 or not at creation time + poisson_ratio = "nu" if nu_sym != 0.5 else 0.5 + + for i in range(dim): + for j in range(i, dim): + self.kernel_dict[(i, j)] = ElasticityKernel(dim=dim, icomp=i, + jcomp=j, poisson_ratio=poisson_ratio) + + # The dictionary allows us to exploit symmetry -- that + # :math:`T_{01}` is identical to :math:`T_{10}` -- and avoid creating + # multiple expansions for the same kernel in a different ordering. + for i in range(dim): + for j in range(i): + self.kernel_dict[(i, j)] = self.kernel_dict[(j, i)] + + def get_int_g(self, idx, density_sym, dir_vec_sym, qbx_forced_limit, + deriv_dirs): + """ + Returns the Integral of the Stokeslet/Stresslet kernel given by `idx` + and its derivatives. + """ + res = _create_int_g(self.kernel_dict[idx], deriv_dirs, + density=density_sym*dir_vec_sym[idx[-1]], + qbx_forced_limit=qbx_forced_limit, mu=self.mu, + nu=self.nu)/(2*(1-self.nu)) + return res + + def apply(self, density_vec_sym, qbx_forced_limit, extra_deriv_dirs=()): + + sym_expr = np.zeros((self.dim,), dtype=object) + + # For stokeslet, there's no direction vector involved + # passing a list of ones instead to remove its usage. + for comp in range(self.dim): + for i in range(self.dim): + sym_expr[comp] += self.get_int_g((comp, i), + density_vec_sym[i], [1]*self.dim, + qbx_forced_limit, deriv_dirs=extra_deriv_dirs) + + return np.array(rewrite_using_base_kernel(sym_expr, + base_kernel=self.base_kernel)) + + def apply_stress(self, density_vec_sym, dir_vec_sym, qbx_forced_limit): + + sym_expr = np.zeros((self.dim,), dtype=object) + stresslet_obj = StressletWrapper(dim=self.dim, + mu_sym=self.mu, nu_sym=self.nu, method=self.method) + + # For stokeslet, there's no direction vector involved + # passing a list of ones instead to remove its usage. + for comp in range(self.dim): + for i in range(self.dim): + for j in range(self.dim): + # pylint does not like __new__ returning new object + # pylint: disable=no-member + sym_expr[comp] += dir_vec_sym[i] * \ + stresslet_obj.get_int_g((comp, i, j), + density_vec_sym[j], [1]*self.dim, + qbx_forced_limit, deriv_dirs=[]) + # pylint: enable=no-member + + return sym_expr + + +class StokesletWrapperNaive(_StokesletWrapperNaiveOrBiharmonic): + def __init__(self, dim=None, mu_sym=_MU_SYM_DEFAULT, nu_sym=0.5): + super().__init__(dim=dim, mu_sym=mu_sym, nu_sym=nu_sym, + method="naive") + + +class StokesletWrapperBiharmonic(_StokesletWrapperNaiveOrBiharmonic): + def __init__(self, dim=None, mu_sym=_MU_SYM_DEFAULT, nu_sym=0.5): + super().__init__(dim=dim, mu_sym=mu_sym, nu_sym=nu_sym, + method="biharmonic") + + +# }}} + + +# {{{ StressletWrapper Naive and Biharmonic impl + +class _StressletWrapperNaiveOrBiharmonic(StressletWrapperBase): + def __init__(self, dim=None, mu_sym=_MU_SYM_DEFAULT, nu_sym=0.5, + method="biharmonic"): + super().__init__(dim, mu_sym, nu_sym) + if not (dim == 3 or dim == 2): + raise ValueError("unsupported dimension given to StokesletWrapper") + + self.method = method + if method == "biharmonic": + self.base_kernel = BiharmonicKernel(dim) + elif method == "naive": + self.base_kernel = None + else: + raise ValueError("method has to be one of biharmonic/naive") + + self.kernel_dict = {} + + for i in range(dim): + for j in range(i, dim): + for k in range(j, dim): + self.kernel_dict[(i, j, k)] = StressletKernel(dim=dim, icomp=i, + jcomp=j, kcomp=k) + + # The dictionary allows us to exploit symmetry -- that + # :math:`T_{012}` is identical to :math:`T_{120}` -- and avoid creating + # multiple expansions for the same kernel in a different ordering. + for i in range(dim): + for j in range(dim): + for k in range(dim): + if (i, j, k) in self.kernel_dict: + continue + s = tuple(sorted([i, j, k])) + self.kernel_dict[(i, j, k)] = self.kernel_dict[s] + + # For elasticity (nu != 0.5), we need the LaplaceKernel + self.kernel_dict["laplace"] = LaplaceKernel(self.dim) + + def get_int_g(self, idx, density_sym, dir_vec_sym, qbx_forced_limit, + deriv_dirs): + """ + Returns the Integral of the Stresslet kernel given by `idx` + and its derivatives. + """ + + nu = self.nu + kernel_indices = [idx] + dir_vec_indices = [idx[-1]] + coeffs = [1] + extra_deriv_dirs_vec = [[]] + + kernel_indices = [idx, "laplace", "laplace", "laplace"] + dir_vec_indices = [idx[-1], idx[1], idx[0], idx[2]] + coeffs = [1, (1 - 2*nu)/self.dim, -(1 - 2*nu)/self.dim, -(1 - 2*nu)] + extra_deriv_dirs_vec = [[], [idx[0]], [idx[1]], [idx[2]]] + if idx[0] != idx[1]: + coeffs[-1] = 0 + + result = 0 + for kernel_idx, dir_vec_idx, coeff, extra_deriv_dirs in \ + zip(kernel_indices, dir_vec_indices, coeffs, + extra_deriv_dirs_vec): + knl = self.kernel_dict[kernel_idx] + result += _create_int_g(knl, tuple(deriv_dirs) + tuple(extra_deriv_dirs), + density=density_sym*dir_vec_sym[dir_vec_idx], + qbx_forced_limit=qbx_forced_limit, mu=self.mu, nu=self.nu) * \ + coeff + return result/(2*(1 - nu)) + + def apply(self, density_vec_sym, dir_vec_sym, qbx_forced_limit, + extra_deriv_dirs=()): + + sym_expr = np.zeros((self.dim,), dtype=object) + + for comp in range(self.dim): + for i in range(self.dim): + for j in range(self.dim): + sym_expr[comp] += self.get_int_g((comp, i, j), + density_vec_sym[i], dir_vec_sym, + qbx_forced_limit, deriv_dirs=extra_deriv_dirs) + + return np.array(rewrite_using_base_kernel(sym_expr, + base_kernel=self.base_kernel)) + + def apply_stokeslet_and_stresslet(self, stokeslet_density_vec_sym, + stresslet_density_vec_sym, dir_vec_sym, + qbx_forced_limit, stokeslet_weight, stresslet_weight, + extra_deriv_dirs=()): + + stokeslet_obj = StokesletWrapper(dim=self.dim, + mu_sym=self.mu, nu_sym=self.nu, method=self.method) + + sym_expr = 0 + if stresslet_weight != 0: + sym_expr += self.apply(stresslet_density_vec_sym, dir_vec_sym, + qbx_forced_limit, extra_deriv_dirs) * stresslet_weight + if stokeslet_weight != 0: + sym_expr += stokeslet_obj.apply(stokeslet_density_vec_sym, + qbx_forced_limit, extra_deriv_dirs) * stokeslet_weight + + return sym_expr + + def apply_stress(self, density_vec_sym, normal_vec_sym, dir_vec_sym, + qbx_forced_limit): sym_expr = np.empty((self.dim,), dtype=object) @@ -480,25 +488,208 @@ def apply_stress(self, density_vec_sym, normal_vec_sym, dir_vec_sym, sym_grad_matrix = np.empty((self.dim, self.dim), dtype=object) for i in range(self.dim): sym_grad_matrix[:, i] = self.apply_derivative(i, density_vec_sym, - normal_vec_sym, mu_sym, qbx_forced_limit) + normal_vec_sym, qbx_forced_limit) for comp in range(self.dim): # First, add the pressure term: sym_expr[comp] = - dir_vec_sym[comp] * self.apply_pressure( density_vec_sym, normal_vec_sym, - mu_sym, qbx_forced_limit) + qbx_forced_limit) # Now add the velocity derivative components for j in range(self.dim): sym_expr[comp] = sym_expr[comp] + ( - dir_vec_sym[j] * mu_sym * ( + dir_vec_sym[j] * self.mu * ( sym_grad_matrix[comp][j] + sym_grad_matrix[j][comp]) ) return sym_expr + +class StressletWrapperNaive(_StressletWrapperNaiveOrBiharmonic): + def __init__(self, dim=None, mu_sym=_MU_SYM_DEFAULT, nu_sym=0.5): + super().__init__(dim=dim, mu_sym=mu_sym, nu_sym=nu_sym, + method="naive") + + +class StressletWrapperBiharmonic(_StressletWrapperNaiveOrBiharmonic): + def __init__(self, dim=None, mu_sym=_MU_SYM_DEFAULT, nu_sym=0.5): + super().__init__(dim=dim, mu_sym=mu_sym, nu_sym=nu_sym, + method="biharmonic") + +# }}} + + +# {{{ Stokeslet/Stresslet using Laplace (Tornberg) + +class StokesletWrapperTornberg(StokesletWrapperBase): + """A Stresslet wrapper using Tornberg and Greengard's method which + uses Laplace derivatives. + + [1] Tornberg, A. K., & Greengard, L. (2008). A fast multipole method for the + three-dimensional Stokes equations. + Journal of Computational Physics, 227(3), 1613-1619. + """ + + def __init__(self, dim=None, mu_sym=_MU_SYM_DEFAULT, nu_sym=0.5): + self.dim = dim + if dim != 3: + raise ValueError("unsupported dimension given to " + "StokesletWrapperTornberg") + if nu_sym != 0.5: + raise ValueError("nu != 0.5 is not supported") + self.kernel = LaplaceKernel(dim=self.dim) + self.mu = mu_sym + self.nu = nu_sym + + def apply(self, density_vec_sym, qbx_forced_limit, extra_deriv_dirs=()): + stresslet = StressletWrapperTornberg(3, self.mu, self.nu) + return stresslet.apply_stokeslet_and_stresslet(density_vec_sym, + [0]*self.dim, [0]*self.dim, qbx_forced_limit, 1, 0, + extra_deriv_dirs) + + +class StressletWrapperTornberg(StressletWrapperBase): + """A Stresslet wrapper using Tornberg and Greengard's method which + uses Laplace derivatives. + + [1] Tornberg, A. K., & Greengard, L. (2008). A fast multipole method for the + three-dimensional Stokes equations. + Journal of Computational Physics, 227(3), 1613-1619. + """ + def __init__(self, dim=None, mu_sym=_MU_SYM_DEFAULT, nu_sym=0.5): + self.dim = dim + if dim != 3: + raise ValueError("unsupported dimension given to " + "StressletWrapperTornberg") + if nu_sym != 0.5: + raise ValueError("nu != 0.5 is not supported") + self.kernel = LaplaceKernel(dim=self.dim) + self.mu = mu_sym + self.nu = nu_sym + + def apply(self, density_vec_sym, dir_vec_sym, qbx_forced_limit, + extra_deriv_dirs=()): + return self.apply_stokeslet_and_stresslet([0]*self.dim, + density_vec_sym, dir_vec_sym, qbx_forced_limit, 0, 1, extra_deriv_dirs) + + def apply_stokeslet_and_stresslet(self, stokeslet_density_vec_sym, + stresslet_density_vec_sym, dir_vec_sym, + qbx_forced_limit, stokeslet_weight, stresslet_weight, + extra_deriv_dirs=()): + + sym_expr = np.zeros((self.dim,), dtype=object) + + source = [sym.NodeCoordinateComponent(d) for d in range(self.dim)] + common_source_kernels = [AxisSourceDerivative(k, self.kernel) for + k in range(self.dim)] + common_source_kernels.append(self.kernel) + + stresslet_weight *= 3.0/6 + stokeslet_weight *= -0.5*self.mu**(-1) + + for i in range(self.dim): + for j in range(self.dim): + densities = [stresslet_weight*( + stresslet_density_vec_sym[k] * dir_vec_sym[j] + + stresslet_density_vec_sym[j] * dir_vec_sym[k]) + for k in range(self.dim)] + densities.append(stokeslet_weight*stokeslet_density_vec_sym[j]) + target_kernel = TargetPointMultiplier(j, + AxisTargetDerivative(i, self.kernel)) + for deriv_dir in extra_deriv_dirs: + target_kernel = AxisTargetDerivative(deriv_dir, target_kernel) + sym_expr[i] -= sym.IntG(target_kernel=target_kernel, + source_kernels=tuple(common_source_kernels), + densities=tuple(densities), + qbx_forced_limit=qbx_forced_limit) + + if i == j: + target_kernel = self.kernel + for deriv_dir in extra_deriv_dirs: + target_kernel = AxisTargetDerivative( + deriv_dir, target_kernel) + + sym_expr[i] += sym.IntG(target_kernel=target_kernel, + source_kernels=common_source_kernels, + densities=densities, + qbx_forced_limit=qbx_forced_limit) + + common_density0 = sum(source[k] * stresslet_density_vec_sym[k] for + k in range(self.dim)) + common_density1 = sum(source[k] * dir_vec_sym[k] for + k in range(self.dim)) + common_density2 = sum(source[k] * stokeslet_density_vec_sym[k] for + k in range(self.dim)) + densities = [stresslet_weight*(common_density0 * dir_vec_sym[k] + + common_density1 * stresslet_density_vec_sym[k]) for + k in range(self.dim)] + densities.append(stokeslet_weight * common_density2) + + target_kernel = AxisTargetDerivative(i, self.kernel) + for deriv_dir in extra_deriv_dirs: + target_kernel = AxisTargetDerivative(deriv_dir, target_kernel) + sym_expr[i] += sym.IntG(target_kernel=target_kernel, + source_kernels=tuple(common_source_kernels), + densities=tuple(densities), + qbx_forced_limit=qbx_forced_limit) + + return sym_expr + +# }}} + + +# {{{ StokesletWrapper dispatch class + +class StokesletWrapper(StokesletWrapperBase): + def __new__(cls, dim=None, mu_sym=_MU_SYM_DEFAULT, nu_sym=0.5, method=None): + if method is None: + import warnings + warnings.warn("method argument not given. falling back to 'naive'" + "method argument will be required in the future.") + method = "naive" + if method == "naive": + return StokesletWrapperNaive(dim=dim, mu_sym=mu_sym, nu_sym=nu_sym) + elif method == "biharmonic": + return StokesletWrapperBiharmonic(dim=dim, mu_sym=mu_sym, nu_sym=nu_sym) + elif method == "laplace": + if nu_sym == 0.5: + return StokesletWrapperTornberg(dim=dim, + mu_sym=mu_sym, nu_sym=nu_sym) + else: + from pytential.symbolic.elasticity import StokesletWrapperYoshida + return StokesletWrapperYoshida(dim=dim, + mu_sym=mu_sym, nu_sym=nu_sym) + else: + raise ValueError(f"invalid method: {method}." + "Needs to be one of naive, laplace, biharmonic") + + +class StressletWrapper(StressletWrapperBase): + def __new__(cls, dim=None, mu_sym=_MU_SYM_DEFAULT, nu_sym=0.5, method=None): + if method is None: + import warnings + warnings.warn("method argument not given. falling back to 'naive'" + "method argument will be required in the future.") + method = "naive" + if method == "naive": + return StressletWrapperNaive(dim=dim, mu_sym=mu_sym, nu_sym=nu_sym) + elif method == "biharmonic": + return StressletWrapperBiharmonic(dim=dim, mu_sym=mu_sym, nu_sym=nu_sym) + elif method == "laplace": + if nu_sym == 0.5: + return StressletWrapperTornberg(dim=dim, + mu_sym=mu_sym, nu_sym=nu_sym) + else: + from pytential.symbolic.elasticity import StressletWrapperYoshida + return StressletWrapperYoshida(dim=dim, + mu_sym=mu_sym, nu_sym=nu_sym) + else: + raise ValueError(f"invalid method: {method}." + "Needs to be one of naive, laplace, biharmonic") + # }}} @@ -518,20 +709,23 @@ class StokesOperator: .. automethod:: pressure """ - def __init__(self, ambient_dim, side): + def __init__(self, ambient_dim, side, method, mu_sym, nu_sym): """ :arg ambient_dim: dimension of the ambient space. :arg side: :math:`+1` for exterior or :math:`-1` for interior. """ - if side not in [+1, -1]: raise ValueError(f"invalid evaluation side: {side}") self.ambient_dim = ambient_dim self.side = side + self.mu = mu_sym + self.nu = nu_sym - self.stresslet = StressletWrapper(dim=self.ambient_dim) - self.stokeslet = StokesletWrapper(dim=self.ambient_dim) + self.stresslet = StressletWrapper(dim=self.ambient_dim, + mu_sym=mu_sym, nu_sym=nu_sym, method=method) + self.stokeslet = StokesletWrapper(dim=self.ambient_dim, + mu_sym=mu_sym, nu_sym=nu_sym, method=method) @property def dim(self): @@ -543,7 +737,7 @@ def get_density_var(self, name="sigma"): """ return sym.make_sym_vector(name, self.ambient_dim) - def prepare_rhs(self, b, *, mu): + def prepare_rhs(self, b): """ :returns: a (potentially) modified right-hand side *b* that matches requirements of the representation. @@ -557,13 +751,13 @@ def operator(self, sigma): """ raise NotImplementedError - def velocity(self, sigma, *, normal, mu, qbx_forced_limit=None): + def velocity(self, sigma, *, normal, qbx_forced_limit=None): """ :returns: a representation of the velocity field in the Stokes flow. """ raise NotImplementedError - def pressure(self, sigma, *, normal, mu, qbx_forced_limit=None): + def pressure(self, sigma, *, normal, qbx_forced_limit=None): """ :returns: a representation of the pressure in the Stokes flow. """ @@ -587,7 +781,8 @@ class HsiaoKressExteriorStokesOperator(StokesOperator): .. automethod:: __init__ """ - def __init__(self, *, omega, alpha=None, eta=None): + def __init__(self, *, omega, alpha=1.0, eta=1.0, method="biharmonic", + mu_sym=_MU_SYM_DEFAULT, nu_sym=0.5): r""" :arg omega: farfield behaviour of the velocity field, as defined by :math:`A` in [HsiaoKress1985]_ Equation 2.3. @@ -595,7 +790,8 @@ def __init__(self, *, omega, alpha=None, eta=None): :arg eta: real parameter :math:`\eta > 0`. Choosing this parameter well can have a non-trivial effect on the conditioning. """ - super().__init__(ambient_dim=2, side=+1) + super().__init__(ambient_dim=2, side=+1, method=method, + mu_sym=mu_sym, nu_sym=nu_sym) # NOTE: in [hsiao-kress], there is an analysis on a circle, which # recommends values in @@ -603,61 +799,60 @@ def __init__(self, *, omega, alpha=None, eta=None): # so we choose alpha = eta = 1, which seems to be in line with some # of the presented numerical results too. - if alpha is None: - alpha = 1.0 - - if eta is None: - eta = 1.0 - self.omega = omega self.alpha = alpha self.eta = eta - def _farfield(self, mu, qbx_forced_limit): - length = sym.integral(self.ambient_dim, self.dim, 1) - return self.stokeslet.apply( - -self.omega / length, - mu, - qbx_forced_limit=qbx_forced_limit) - - def _operator(self, sigma, normal, mu, qbx_forced_limit): + def _farfield(self, qbx_forced_limit): + source_dofdesc = sym.DOFDescriptor(None, discr_stage=sym.QBX_SOURCE_STAGE1) + length = sym.integral(self.ambient_dim, self.dim, 1, dofdesc=source_dofdesc) + # pylint does not like __new__ returning new object + # pylint: disable=no-member + result = self.stresslet.apply_stokeslet_and_stresslet( + -self.omega / length, [0]*self.ambient_dim, [0]*self.ambient_dim, + qbx_forced_limit=qbx_forced_limit, stokeslet_weight=1, + stresslet_weight=0) + # pylint: enable=no-member + return result + + def _operator(self, sigma, normal, qbx_forced_limit): slp_qbx_forced_limit = qbx_forced_limit if slp_qbx_forced_limit == "avg": - slp_qbx_forced_limit = +1 + slp_qbx_forced_limit = "avg" # NOTE: we set a dofdesc here to force the evaluation of this integral # on the source instead of the target when using automatic tagging # see :meth:`pytential.symbolic.mappers.LocationTagger._default_dofdesc` dd = sym.DOFDescriptor(None, discr_stage=sym.QBX_SOURCE_STAGE1) - int_sigma = sym.integral(self.ambient_dim, self.dim, sigma, dofdesc=dd) - meanless_sigma = sym.cse(sigma - sym.mean(self.ambient_dim, self.dim, sigma)) - op_k = self.stresslet.apply(sigma, normal, mu, - qbx_forced_limit=qbx_forced_limit) - op_s = ( - self.alpha / (2.0 * np.pi) * int_sigma - - self.stokeslet.apply(meanless_sigma, mu, - qbx_forced_limit=slp_qbx_forced_limit) - ) + meanless_sigma = sym.cse(sigma - sym.mean(self.ambient_dim, + self.dim, sigma, dofdesc=dd)) - return op_k + self.eta * op_s + result = self.eta * self.alpha / (2.0 * np.pi) * int_sigma + # pylint does not like __new__ returning new object + # pylint: disable=no-member + result += self.stresslet.apply_stokeslet_and_stresslet(meanless_sigma, + sigma, normal, qbx_forced_limit=qbx_forced_limit, + stokeslet_weight=-self.eta, stresslet_weight=1) + # pylint: enable=no-member - def prepare_rhs(self, b, *, mu): - return b + self._farfield(mu, qbx_forced_limit=+1) + return result - def operator(self, sigma, *, normal, mu): + def prepare_rhs(self, b): + return b + self._farfield(qbx_forced_limit=+1) + + def operator(self, sigma, *, normal, qbx_forced_limit="avg"): # NOTE: H. K. 1985 Equation 2.18 - return -0.5 * self.side * sigma - self._operator(sigma, normal, mu, "avg") + return -0.5 * self.side * sigma - self._operator( + sigma, normal, qbx_forced_limit) - def velocity(self, sigma, *, normal, mu, qbx_forced_limit=2): + def velocity(self, sigma, *, normal, qbx_forced_limit=2): # NOTE: H. K. 1985 Equation 2.16 - return ( - -self._farfield(mu, qbx_forced_limit) - - self._operator(sigma, normal, mu, qbx_forced_limit) - ) + return -self._farfield(qbx_forced_limit) \ + - self._operator(sigma, normal, qbx_forced_limit) - def pressure(self, sigma, *, normal, mu, qbx_forced_limit=2): + def pressure(self, sigma, *, normal, qbx_forced_limit=2): # FIXME: H. K. 1985 Equation 2.17 raise NotImplementedError @@ -675,13 +870,15 @@ class HebekerExteriorStokesOperator(StokesOperator): .. automethod:: __init__ """ - def __init__(self, *, eta=None): + def __init__(self, *, eta=None, method="laplace", mu_sym=_MU_SYM_DEFAULT, + nu_sym=0.5): r""" :arg eta: a parameter :math:`\eta > 0`. Choosing this parameter well can have a non-trivial effect on the conditioning of the operator. """ - super().__init__(ambient_dim=3, side=+1) + super().__init__(ambient_dim=3, side=+1, method=method, + mu_sym=mu_sym, nu_sym=nu_sym) # NOTE: eta is chosen here based on H. 1986 Figure 1, which is # based on solving on the unit sphere @@ -689,28 +886,28 @@ def __init__(self, *, eta=None): eta = 0.75 self.eta = eta - - def _operator(self, sigma, normal, mu, qbx_forced_limit): - slp_qbx_forced_limit = qbx_forced_limit - if slp_qbx_forced_limit == "avg": - slp_qbx_forced_limit = self.side - - op_w = self.stresslet.apply(sigma, normal, mu, - qbx_forced_limit=qbx_forced_limit) - op_v = self.stokeslet.apply(sigma, mu, - qbx_forced_limit=slp_qbx_forced_limit) - - return op_w + self.eta * op_v - - def operator(self, sigma, *, normal, mu): + self.laplace_kernel = LaplaceKernel(3) + + def _operator(self, sigma, normal, qbx_forced_limit): + # pylint does not like __new__ returning new object + # pylint: disable=no-member + result = self.stresslet.apply_stokeslet_and_stresslet(sigma, + sigma, normal, qbx_forced_limit=qbx_forced_limit, + stokeslet_weight=self.eta, stresslet_weight=1, + extra_deriv_dirs=()) + # pylint: enable=no-member + return result + + def operator(self, sigma, *, normal, qbx_forced_limit="avg"): # NOTE: H. 1986 Equation 17 - return -0.5 * self.side * sigma - self._operator(sigma, normal, mu, "avg") + return -0.5 * self.side * sigma - self._operator(sigma, + normal, qbx_forced_limit) - def velocity(self, sigma, *, normal, mu, qbx_forced_limit=2): + def velocity(self, sigma, *, normal, qbx_forced_limit=2): # NOTE: H. 1986 Equation 16 - return -self._operator(sigma, normal, mu, qbx_forced_limit) + return -self._operator(sigma, normal, qbx_forced_limit) - def pressure(self, sigma, *, normal, mu, qbx_forced_limit=2): + def pressure(self, sigma, *, normal, qbx_forced_limit=2): # FIXME: not given in H. 1986, but should be easy to derive using the # equivalent single-/double-layer pressure kernels raise NotImplementedError diff --git a/pytential/utils.py b/pytential/utils.py index b64176909..64c369ba4 100644 --- a/pytential/utils.py +++ b/pytential/utils.py @@ -1,5 +1,5 @@ __copyright__ = """ -Copyright (C) 2020 Matt Wala +Copyright (C) 2020 Isuru Fernando """ __license__ = """ @@ -22,6 +22,8 @@ THE SOFTWARE. """ +import sumpy.symbolic as sym + def sort_arrays_together(*arys, key=None): """Sort a sequence of arrays by considering them @@ -32,4 +34,61 @@ def sort_arrays_together(*arys, key=None): """ return zip(*sorted([x for x in zip(*arys)], key=key)) + +def chop(expr, tol): + """Given a symbolic expression, remove all occurences of numbers + with absolute value less than a given tolerance and replace floating + point numbers that are close to an integer up to a given relative + tolerance by the integer. + """ + nums = expr.atoms(sym.Number) + replace_dict = {} + for num in nums: + if float(abs(num)) < tol: + replace_dict[num] = 0 + else: + new_num = float(num) + if abs((int(new_num) - new_num)/new_num) < tol: + new_num = int(new_num) + replace_dict[num] = new_num + return expr.xreplace(replace_dict) + + +def lu_solve_with_expand(L, U, perm, b): + """Given an LU factorization and a vector, solve a linear + system with intermediate results expanded to avoid + an explosion of the expression trees + + :param L: lower triangular matrix + :param U: upper triangular matrix + :param perm: permutation matrix + :param b: column vector to solve for + """ + def forward_substitution(L, b): + n = len(b) + res = sym.Matrix(b) + for i in range(n): + for j in range(i): + res[i] -= L[i, j]*res[j] + res[i] = (res[i] / L[i, i]).expand() + return res + + def backward_substitution(U, b): + n = len(b) + res = sym.Matrix(b) + for i in range(n-1, -1, -1): + for j in range(n - 1, i, -1): + res[i] -= U[i, j]*res[j] + res[i] = (res[i] / U[i, i]).expand() + return res + + def permuteFwd(b, perm): + res = sym.Matrix(b) + for p, q in perm: + res[p], res[q] = res[q], res[p] + return res + + return backward_substitution(U, + forward_substitution(L, permuteFwd(b, perm))) + # vim: foldmethod=marker diff --git a/requirements.txt b/requirements.txt index 88777f180..456c4e61e 100644 --- a/requirements.txt +++ b/requirements.txt @@ -11,5 +11,5 @@ git+https://github.com/inducer/loopy.git#egg=loopy git+https://github.com/inducer/boxtree.git#egg=boxtree git+https://github.com/inducer/arraycontext.git#egg=arraycontext git+https://github.com/inducer/meshmode.git#egg=meshmode -git+https://github.com/inducer/sumpy.git#egg=sumpy +git+https://github.com/isuruf/sumpy.git@spatial_constant#egg=sumpy git+https://github.com/inducer/pyfmmlib.git#egg=pyfmmlib diff --git a/test/test_stokes.py b/test/test_stokes.py index dca8a4dbe..ed68dcb43 100644 --- a/test/test_stokes.py +++ b/test/test_stokes.py @@ -31,17 +31,17 @@ from meshmode.discretization.poly_element import \ InterpolatoryQuadratureGroupFactory from pytools.obj_array import make_obj_array +from sumpy.symbolic import SpatialConstant from meshmode import _acf # noqa: F401 from arraycontext import pytest_generate_tests_for_array_contexts -from meshmode.array_context import PytestPyOpenCLArrayContextFactory import extra_int_eq_data as eid import logging logger = logging.getLogger(__name__) pytest_generate_tests = pytest_generate_tests_for_array_contexts([ - PytestPyOpenCLArrayContextFactory, + "pyopencl-deprecated", ]) @@ -68,11 +68,13 @@ def dof_array_rel_error(actx, x, xref, p=None): def run_exterior_stokes(actx_factory, *, ambient_dim, target_order, qbx_order, resolution, - fmm_order=False, # FIXME: FMM is slower than direct evaluation + fmm_order=None, # FIXME: FMM is slower than direct evaluation source_ovsmp=None, radius=1.5, aspect_ratio=1.0, mu=1.0, + nu=0.4, visualize=False, + method="naive", _target_association_tolerance=0.05, _expansions_in_tree_have_extent=True): @@ -152,27 +154,41 @@ def run_exterior_stokes(actx_factory, *, # {{{ symbolic sym_normal = sym.make_sym_vector("normal", ambient_dim) - sym_mu = sym.var("mu") + sym_mu = SpatialConstant("mu2") + + if nu == 0.5: + sym_nu = 0.5 + else: + sym_nu = SpatialConstant("nu2") if ambient_dim == 2: from pytential.symbolic.stokes import HsiaoKressExteriorStokesOperator sym_omega = sym.make_sym_vector("omega", ambient_dim) - op = HsiaoKressExteriorStokesOperator(omega=sym_omega) + op = HsiaoKressExteriorStokesOperator(omega=sym_omega, method=method, + mu_sym=sym_mu, nu_sym=sym_nu) elif ambient_dim == 3: from pytential.symbolic.stokes import HebekerExteriorStokesOperator - op = HebekerExteriorStokesOperator() + op = HebekerExteriorStokesOperator(method=method, + mu_sym=sym_mu, nu_sym=sym_nu) else: raise AssertionError() sym_sigma = op.get_density_var("sigma") sym_bc = op.get_density_var("bc") - sym_op = op.operator(sym_sigma, normal=sym_normal, mu=sym_mu) - sym_rhs = op.prepare_rhs(sym_bc, mu=mu) + sym_op = op.operator(sym_sigma, normal=sym_normal) + sym_rhs = op.prepare_rhs(sym_bc) - sym_velocity = op.velocity(sym_sigma, normal=sym_normal, mu=sym_mu) + sym_velocity = op.velocity(sym_sigma, normal=sym_normal) - sym_source_pot = op.stokeslet.apply(sym_sigma, sym_mu, qbx_forced_limit=None) + if ambient_dim == 3: + sym_source_pot = op.stokeslet.apply(sym_sigma, qbx_forced_limit=None) + else: + # Use the naive method here as biharmonic requires source derivatives + # of point_source + from pytential.symbolic.stokes import StokesletWrapper + sym_source_pot = StokesletWrapper(ambient_dim, mu_sym=sym_mu, + nu_sym=sym_nu, method="naive").apply(sym_sigma, qbx_forced_limit=None) # }}} @@ -193,20 +209,43 @@ def run_exterior_stokes(actx_factory, *, omega = bind(places, total_charge * sym.Ones())(actx) if ambient_dim == 2: - bc_context = {"mu": mu, "omega": omega} - op_context = {"mu": mu, "omega": omega, "normal": normal} + bc_context = {"mu2": mu, "omega": omega} + op_context = {"mu2": mu, "omega": omega, "normal": normal} else: bc_context = {} - op_context = {"mu": mu, "normal": normal} + op_context = {"mu2": mu, "normal": normal} + direct_context = {"mu2": mu} - bc = bind(places, sym_source_pot, - auto_where=("point_source", "source"))(actx, sigma=charges, mu=mu) + if sym_nu != 0.5: + bc_context["nu2"] = nu + op_context["nu2"] = nu + direct_context["nu2"] = nu + + bc_op = bind(places, sym_source_pot, + auto_where=("point_source", "source")) + bc = bc_op(actx, sigma=charges, **direct_context) rhs = bind(places, sym_rhs)(actx, bc=bc, **bc_context) bound_op = bind(places, sym_op) # }}} + fmm_timing_data = {} + bound_op.eval({"sigma": rhs, **op_context}, array_context=actx, + timing_data=fmm_timing_data) + + def print_timing_data(timings, name): + result = {k: 0 for k in list(timings.values())[0].keys()} + total = 0 + for k, timing in timings.items(): + for k, v in timing.items(): + result[k] += v["wall_elapsed"] + total += v["wall_elapsed"] + result["total"] = total + print(f"{name}={result}") + + # print_timing_data(fmm_timing_data, method) + # {{{ solve from pytential.solve import gmres @@ -219,7 +258,6 @@ def run_exterior_stokes(actx_factory, *, progress=visualize, stall_iterations=0, hard_failure=True) - sigma = result.solution # }}} @@ -229,7 +267,8 @@ def run_exterior_stokes(actx_factory, *, velocity = bind(places, sym_velocity, auto_where=("source", "point_target"))(actx, sigma=sigma, **op_context) ref_velocity = bind(places, sym_source_pot, - auto_where=("point_source", "point_target"))(actx, sigma=charges, mu=mu) + auto_where=("point_source", "point_target"))(actx, sigma=charges, + **direct_context) v_error = [ dof_array_rel_error(actx, u, uref) @@ -265,11 +304,18 @@ def run_exterior_stokes(actx_factory, *, return h_max, v_error -@pytest.mark.parametrize("ambient_dim", [ - 2, - pytest.param(3, marks=pytest.mark.slowtest) +@pytest.mark.parametrize("ambient_dim, method, nu", [ + (2, "naive", 0.5), + (2, "biharmonic", 0.5), + pytest.param(3, "naive", 0.5, marks=pytest.mark.slowtest), + pytest.param(3, "biharmonic", 0.5, marks=pytest.mark.slowtest), + pytest.param(3, "laplace", 0.5, marks=pytest.mark.slowtest), + + (2, "biharmonic", 0.4), + pytest.param(3, "biharmonic", 0.4, marks=pytest.mark.slowtest), + pytest.param(3, "laplace", 0.4, marks=pytest.mark.slowtest), ]) -def test_exterior_stokes(actx_factory, ambient_dim, visualize=False): +def test_exterior_stokes(actx_factory, ambient_dim, method, nu, visualize=False): if visualize: logging.basicConfig(level=logging.INFO) @@ -281,8 +327,10 @@ def test_exterior_stokes(actx_factory, ambient_dim, visualize=False): qbx_order = 3 if ambient_dim == 2: + fmm_order = 10 resolutions = [20, 35, 50] elif ambient_dim == 3: + fmm_order = 6 resolutions = [0, 1, 2] else: raise ValueError(f"unsupported dimension: {ambient_dim}") @@ -291,9 +339,12 @@ def test_exterior_stokes(actx_factory, ambient_dim, visualize=False): h_max, errors = run_exterior_stokes(actx_factory, ambient_dim=ambient_dim, target_order=target_order, + fmm_order=fmm_order, qbx_order=qbx_order, source_ovsmp=source_ovsmp, resolution=resolution, + method=method, + nu=nu, visualize=visualize) for eoc, e in zip(eocs, errors): @@ -305,12 +356,17 @@ def test_exterior_stokes(actx_factory, ambient_dim, visualize=False): error_format="%.8e", eoc_format="%.2f")) + extra_order = 0 + if method == "biharmonic": + extra_order += 1 + elif nu != 0.5: + extra_order += 0.5 for eoc in eocs: # This convergence data is not as clean as it could be. See # https://github.com/inducer/pytential/pull/32 # for some discussion. order = min(target_order, qbx_order) - assert eoc.order_estimate() > order - 0.5 + assert eoc.order_estimate() > order - 0.5 - extra_order # }}} @@ -404,13 +460,13 @@ class StokesletIdentity: def __init__(self, ambient_dim): from pytential.symbolic.stokes import StokesletWrapper self.ambient_dim = ambient_dim - self.stokeslet = StokesletWrapper(self.ambient_dim) + self.stokeslet = StokesletWrapper(self.ambient_dim, mu_sym=1) def apply_operator(self): sym_density = sym.normal(self.ambient_dim).as_vector() return self.stokeslet.apply( sym_density, - mu_sym=1, qbx_forced_limit=+1) + qbx_forced_limit=+1) def ref_result(self): return make_obj_array([1.0e-15 * sym.Ones()] * self.ambient_dim) @@ -463,13 +519,13 @@ class StressletIdentity: def __init__(self, ambient_dim): from pytential.symbolic.stokes import StokesletWrapper self.ambient_dim = ambient_dim - self.stokeslet = StokesletWrapper(self.ambient_dim) + self.stokeslet = StokesletWrapper(self.ambient_dim, mu_sym=1) def apply_operator(self): sym_density = sym.normal(self.ambient_dim).as_vector() return self.stokeslet.apply_stress( sym_density, sym_density, - mu_sym=1, qbx_forced_limit="avg") + qbx_forced_limit="avg") def ref_result(self): return -0.5 * sym.normal(self.ambient_dim).as_vector() @@ -520,6 +576,14 @@ def test_stresslet_identity(actx_factory, cls, visualize=False): if __name__ == "__main__": import sys if len(sys.argv) > 1: + import pyopencl as cl + from arraycontext import PyOpenCLArrayContext + context = cl._csc() + queue = cl.CommandQueue(context) + + def actx_factory(): + return PyOpenCLArrayContext(queue) + exec(sys.argv[1]) else: from pytest import main