diff --git a/doc/index.rst b/doc/index.rst index 598cd4641..cdfa3cc26 100644 --- a/doc/index.rst +++ b/doc/index.rst @@ -23,6 +23,7 @@ Contents tools misc qbx + merge 🚀 Github 💾 Download Releases diff --git a/doc/merge.rst b/doc/merge.rst new file mode 100644 index 000000000..2c4b15ea2 --- /dev/null +++ b/doc/merge.rst @@ -0,0 +1,4 @@ +Utilities to merge and reduce IntG expressions +============================================== + +.. automodule:: pytential.symbolic.pde.system_utils diff --git a/doc/private.rst b/doc/private.rst new file mode 100644 index 000000000..c84fae8ad --- /dev/null +++ b/doc/private.rst @@ -0,0 +1,14 @@ +:orphan: + +Internal APIs +============= + +These are private API that are documented for developers +and should not be used by end users. + +Reduce Number of FMMs +--------------------- + +.. automodule:: pytential.symbolic.pde.reduce_fmms + + diff --git a/experiments/stokes-2d-interior.py b/experiments/stokes-2d-interior.py index bcb15fde8..57727d273 100644 --- a/experiments/stokes-2d-interior.py +++ b/experiments/stokes-2d-interior.py @@ -21,6 +21,8 @@ qbx_order = 2 fmm_order = 7 mu = 3 +# method has to be one of biharmonic/naive for 2D +method = "biharmonic" # Test solution type -- either 'fundamental' or 'couette' (default is couette) soln_type = 'couette' @@ -87,10 +89,10 @@ def get_obj_array(obj_array): loc_sign = -1 # Create stresslet object - stresslet_obj = StressletWrapper(dim=2) + stresslet_obj = StressletWrapper(dim=2, mu_sym=mu_sym, method=method) # Describe boundary operator - bdry_op_sym = loc_sign * 0.5 * sigma_sym + sqrt_w * stresslet_obj.apply(inv_sqrt_w_sigma, nvec_sym, mu_sym, qbx_forced_limit='avg') + bdry_op_sym = loc_sign * 0.5 * sigma_sym + sqrt_w * stresslet_obj.apply(inv_sqrt_w_sigma, nvec_sym, qbx_forced_limit='avg') # Bind to the qbx discretization bound_op = bind(qbx, bdry_op_sym) @@ -140,7 +142,7 @@ def couette_soln(x, y, dp, h): sigma = gmres_result.solution # Describe representation of solution for evaluation in domain - representation_sym = stresslet_obj.apply(inv_sqrt_w_sigma, nvec_sym, mu_sym, qbx_forced_limit=-2) + representation_sym = stresslet_obj.apply(inv_sqrt_w_sigma, nvec_sym, qbx_forced_limit=-2) from sumpy.visualization import FieldPlotter nsamp = 10 @@ -194,11 +196,11 @@ def stride_hack(arr): print("exact velocity at max error points: x -> ", err[0][max_error_loc[0]], ", y -> ", err[1][max_error_loc[1]]) from pytential.symbolic.mappers import DerivativeTaker - rep_pressure = stresslet_obj.apply_pressure(inv_sqrt_w_sigma, nvec_sym, mu_sym, qbx_forced_limit=-2) + rep_pressure = stresslet_obj.apply_pressure(inv_sqrt_w_sigma, nvec_sym, qbx_forced_limit=-2) pressure = bind((qbx, PointsTarget(eval_points_dev)), rep_pressure)(queue, sigma=sigma, mu=mu, normal=normal) pressure = pressure.get() - print "pressure = ", pressure + print(f"pressure = {pressure}") x_dir_vecs = np.zeros((2,len(eval_points[0]))) x_dir_vecs[0,:] = 1.0 @@ -207,7 +209,7 @@ def stride_hack(arr): x_dir_vecs = cl.array.to_device(queue, x_dir_vecs) y_dir_vecs = cl.array.to_device(queue, y_dir_vecs) dir_vec_sym = sym.make_sym_vector("force_direction", dim) - rep_stress = stresslet_obj.apply_stress(inv_sqrt_w_sigma, nvec_sym, dir_vec_sym, mu_sym, qbx_forced_limit=-2) + rep_stress = stresslet_obj.apply_stress(inv_sqrt_w_sigma, nvec_sym, dir_vec_sym, qbx_forced_limit=-2) applied_stress_x = bind((qbx, PointsTarget(eval_points_dev)), rep_stress)(queue, sigma=sigma, normal=normal, force_direction=x_dir_vecs, mu=mu) @@ -216,8 +218,8 @@ def stride_hack(arr): rep_stress)(queue, sigma=sigma, normal=normal, force_direction=y_dir_vecs, mu=mu) applied_stress_y = get_obj_array(applied_stress_y) - print "stress applied to x direction: ", applied_stress_x - print "stress applied to y direction: ", applied_stress_y + print(f"stress applied to x direction: {applied_stress_x}") + print(f"stress applied to y direction: {applied_stress_y}") import matplotlib.pyplot as plt diff --git a/pytential/qbx/fmm.py b/pytential/qbx/fmm.py index 91cce641a..d507eca16 100644 --- a/pytential/qbx/fmm.py +++ b/pytential/qbx/fmm.py @@ -112,7 +112,7 @@ class QBXExpansionWrangler(SumpyExpansionWrangler): def __init__(self, tree_indep, geo_data, dtype, qbx_order, fmm_level_to_order, source_extra_kwargs, kernel_extra_kwargs, - _use_target_specific_qbx=None): + *, translation_classes_data=None, _use_target_specific_qbx=None): if _use_target_specific_qbx: raise ValueError("TSQBX is not implemented in sumpy") diff --git a/pytential/source.py b/pytential/source.py index abf64d6bc..5ee30efbe 100644 --- a/pytential/source.py +++ b/pytential/source.py @@ -143,6 +143,7 @@ def ambient_dim(self): return self._nodes.shape[0] def op_group_features(self, expr): + from sumpy.kernel import TargetTransformationRemover from pytential.utils import sort_arrays_together # since IntGs with the same source kernels and densities calculations # for P2E and E2E are the same and only differs in E2P depending on the @@ -152,7 +153,7 @@ def op_group_features(self, expr): result = ( expr.source, *sort_arrays_together(expr.source_kernels, expr.densities, key=str), - expr.target_kernel, + TargetTransformationRemover()(expr.target_kernel), ) return result 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/execution.py b/pytential/symbolic/execution.py index c58de1762..c4816ccea 100644 --- a/pytential/symbolic/execution.py +++ b/pytential/symbolic/execution.py @@ -843,10 +843,12 @@ def __call__(self, *args, **kwargs): raise TypeError("More than one positional argument supplied. " "None or an PyOpenCLArrayContext expected.") - return self.eval(kwargs, array_context=array_context) + timing_data = kwargs.pop("timing_data", None) + return self.eval(kwargs, array_context=array_context, + timing_data=timing_data) -def bind(places, expr, auto_where=None): +def bind(places, expr, auto_where=None, _merge_exprs=True): """ :arg places: a :class:`pytential.collection.GeometryCollection`. Alternatively, any list or mapping that is a valid argument for its @@ -865,6 +867,18 @@ def bind(places, expr, auto_where=None): auto_where = places.auto_where expr = _prepare_expr(places, expr, auto_where=auto_where) + if _merge_exprs: + from pytential.symbolic.pde.system_utils import merge_int_g_exprs + from pymbolic.primitives import Expression + from pytential.qbx import QBXLayerPotentialSource + fmmlib = any(value.fmm_backend == "fmmlib" for value + in places.places.values() if isinstance(value, QBXLayerPotentialSource)) + if not fmmlib: + if isinstance(expr, (np.ndarray, list, tuple)): + expr = np.array(merge_int_g_exprs(list(expr)), dtype=object) + elif isinstance(expr, Expression): + expr = merge_int_g_exprs([expr])[0] + return BoundExpression(places, expr) # }}} diff --git a/pytential/symbolic/pde/reduce_fmms.py b/pytential/symbolic/pde/reduce_fmms.py new file mode 100644 index 000000000..6bfab4f3e --- /dev/null +++ b/pytential/symbolic/pde/reduce_fmms.py @@ -0,0 +1,497 @@ +__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. +""" + +from sumpy.kernel import (AxisTargetDerivative, AxisSourceDerivative, + KernelWrapper) + +from pymbolic.interop.sympy import PymbolicToSympyMapper, SympyToPymbolicMapper +from pymbolic.mapper import Mapper +from pymbolic.geometric_algebra.mapper import WalkMapper +from pymbolic.primitives import Product +import sympy +import functools +from collections import defaultdict + +import logging +logger = logging.getLogger(__name__) + + +__all__ = ( + "reduce_number_of_fmms", + ) + +__doc__ = """ +.. autofunction:: reduce_number_of_fmms +""" + + +# {{{ Reduce number of FMMs - main routine + +def reduce_number_of_fmms(int_gs, source_dependent_variables): + """ + Reduce the number of FMMs needed for a system of + :class:`~pytential.symbolic.primitives.IntG` objects. + + This is done by converting the ``IntG`` object to a matrix of polynomials + with d variables corresponding to d dimensions, where each variable represents + a (target) derivative operator along one of the axes. All the properties of + derivative operator that we want are reflected in the properties of the + polynomial including addition, multiplication and exact polynomial division. + + This matrix is factored into two matrices, where the left hand side matrix + represents a transformation at the target, and the right hand side matrix + represents a transformation at the source. + + If the expressions given are not linear, then the input expressions are + returned as is. + + :arg int_gs: list of ``IntG`` objects. + + :arg source_dependent_variables: list of :class:`pymbolic.primitives.Expression` + objects. When reducing FMMs, consider only these variables as dependent + on source. For eg: densities, source derivative vectors. + + Note: there is no argument for target-dependent variables as the algorithm + assumes that there are no target-dependent variables passed to this function. + (where a "source/target-dependent variable" is a symbolic variable that evaluates + to a vector discretized on the sources/targets) + """ + + dim = int_gs[0].target_kernel.dim + axis_vars = sympy.symbols(f"_x0:{dim}") + + # A high level driver for this function should send int_gs that share all the + # properties except for the densities, source transformations and target + # transformations. + assert _check_int_gs_common(int_gs) + + if source_dependent_variables is None: + source_dependent_variables = get_all_source_dependent_variables(int_gs) + + try: + mat, source_exprs = _create_matrix(int_gs, source_dependent_variables, + axis_vars) + except ValueError: + logger.debug("could not create matrix from %s", int_gs) + return int_gs + + mat = sympy.Matrix(mat) + + # Factor the matrix into two + try: + left_factor = _factor_left(mat, axis_vars) + right_factor = _factor_right(mat, left_factor) + except ValueError: + logger.debug("could not find a factorization for %s", mat) + return int_gs + + # If there are n inputs and m outputs, + # + # - matrix: R^{m x n}, + # - LHS: R^{m x k}, + # - RHS: R^{k x n}. + # + # If k is greater than or equal to n we are gaining nothing. + # Return as is. + if right_factor.shape[0] >= mat.shape[0]: + return int_gs + + base_kernel = int_gs[0].source_kernels[0].get_base_kernel() + base_int_g = int_gs[0].copy(target_kernel=base_kernel, + source_kernels=(base_kernel,), densities=(1,)) + + # Convert polynomials back to IntGs with source derivatives + source_int_gs = [[_convert_source_poly_to_int_g_derivs( + expr.as_poly(*axis_vars, domain=sympy.EX), base_int_g, + axis_vars) for expr in row] for row in right_factor.tolist()] + + # For each row in the right factor, merge the IntGs to one IntG + # to get a total of k IntGs. + source_int_gs_merged = [] + for i in range(right_factor.shape[0]): + source_kernels = [] + densities = [] + for j in range(right_factor.shape[1]): + new_densities = [density * source_exprs[j] for density in + source_int_gs[i][j].densities] + source_kernels.extend(source_int_gs[i][j].source_kernels) + densities.extend(new_densities) + source_int_gs_merged.append(source_int_gs[i][0].copy( + source_kernels=tuple(source_kernels), densities=tuple(densities))) + + # Now that we have the IntG expressions depending on the source + # we now have to attach the target dependent derivatives. + res = [0]*left_factor.shape[0] + for i in range(left_factor.shape[0]): + for j in range(left_factor.shape[1]): + res[i] += _convert_target_poly_to_int_g_derivs( + left_factor[i, j].as_poly(*axis_vars, domain=sympy.EX), + int_gs[i], source_int_gs_merged[j]) + + return res + + +class GatherAllSourceDependentVariables(WalkMapper): + def __init__(self): + self.vars = {} + + def map_variable(self, expr): + from sumpy.symbolic import SpatialConstant + if not isinstance(expr, SpatialConstant): + self.vars[expr] = True + + def map_list(self, exprs): + for expr in exprs: + self.rec(expr) + + def map_int_g(self, expr): + self.vars[expr] = True + for density in expr.densities: + self.rec(density) + + def map_node_coordinate_component(self, expr): + self.vars[expr] = True + + map_num_reference_derivative = map_node_coordinate_component + map_interpolation = map_node_coordinate_component + + +def get_all_source_dependent_variables(int_gs): + mapper = GatherAllSourceDependentVariables() + for int_g in int_gs: + for density in int_g.densities: + mapper(density) + return list(mapper.vars.keys()) + +# }}} + + +# {{{ convert IntG expressions to a matrix + +def _check_int_gs_common(int_gs): + """Checks that the :class:`~pytential.symbolic.primtive.IntG` objects + have the same base kernel and other properties that would allow + merging them. + """ + from pytential.symbolic.pde.system_utils import merge_kernel_arguments + + kernel_arguments = {} + base_kernel = int_gs[0].source_kernels[0].get_base_kernel() + common_int_g = int_gs[0].copy(target_kernel=base_kernel, + source_kernels=(base_kernel,), densities=(1,)) + + for int_g in int_gs: + for source_kernel in int_g.source_kernels: + if source_kernel.get_base_kernel() != base_kernel: + return False + + if common_int_g.qbx_forced_limit != int_g.qbx_forced_limit: + return False + + if common_int_g.source != int_g.source: + return False + + try: + kernel_arguments = merge_kernel_arguments(kernel_arguments, + int_g.kernel_arguments) + except ValueError: + return False + return True + + +def _create_matrix(int_gs, source_dependent_variables, axis_vars): + """Create a matrix from a list of :class:`~pytential.symbolic.primitives.IntG` + objects and returns the matrix and the expressions corresponding to each column. + Each expression is an expression containing ``source_dependent_variables``. + Each element in the matrix is a multi-variate polynomial and the variables + in the polynomial are from ``axis_vars`` input. Each polynomial represents + a derivative operator. + + Number of rows of the returned matrix is equal to the number of ``int_gs`` and + the number of columns is equal to the number of input source dependent + expressions. + """ + source_exprs = [] + coefficient_collector = CoefficientCollector(source_dependent_variables) + to_sympy = PymbolicToSympyMapper() + matrix = [] + + for int_g in int_gs: + row = [0]*len(source_exprs) + for density, source_kernel in zip(int_g.densities, int_g.source_kernels): + d = coefficient_collector(density) + for source_expr, coeff in d.items(): + if source_expr not in source_exprs: + source_exprs.append(source_expr) + row += [0] + poly = _kernel_source_derivs_as_poly(source_kernel, axis_vars) + row[source_exprs.index(source_expr)] += poly * to_sympy(coeff) + matrix.append(row) + + # At the beginning, we didn't know the number of columns of the matrix. + # Therefore we used a list for rows and they kept expanding. + # Here we are adding zero padding to make the result look like a matrix. + for row in matrix: + row += [0]*(len(source_exprs) - len(row)) + + return matrix, source_exprs + + +class CoefficientCollector(Mapper): + """From a density expression, extracts expressions that need to be + evaluated for each source and coefficients for each expression. + + For eg: when this mapper is given as ``s*(s + 2) + 3`` input, + it returns {s**2: 1, s: 2, 1: 3}. + + This is more general than + :class:`pymbolic.mapper.coefficient.CoefficientCollector` as that deals + only with linear expressions, but this collector works for polynomial + expressions too. + """ + def __init__(self, source_dependent_variables): + self.source_dependent_variables = source_dependent_variables + + def __call__(self, expr): + if expr in self.source_dependent_variables: + return {expr: 1} + return super().__call__(expr) + + def map_sum(self, expr): + stride_dicts = [self.rec(ch) for ch in expr.children] + + result = defaultdict(lambda: 0) + for stride_dict in stride_dicts: + for var, stride in stride_dict.items(): + result[var] += stride + return dict(result) + + def map_algebraic_leaf(self, expr): + if expr in self.source_dependent_variables: + return {expr: 1} + else: + return {1: expr} + + def map_node_coordinate_component(self, expr): + return {expr: 1} + + map_num_reference_derivative = map_node_coordinate_component + + def map_common_subexpression(self, expr): + return {expr: 1} + + def map_subscript(self, expr): + if expr in self.source_dependent_variables or \ + expr.aggregate in self.source_dependent_variables: + return {expr: 1} + else: + return {1: expr} + + def map_constant(self, expr): + return {1: expr} + + def map_product(self, expr): + if len(expr.children) > 2: + # rewrite products of more than two children as a nested + # product and recurse to make it easier to handle. + left = expr.children[0] + right = Product(tuple(expr.children[1:])) + new_prod = Product((left, right)) + return self.rec(new_prod) + elif len(expr.children) == 1: + return self.rec(expr.children[0]) + elif len(expr.children) == 0: + return {1: 1} + left, right = expr.children + d_left = self.rec(left) + d_right = self.rec(right) + d = defaultdict(lambda: 0) + for var_left, coeff_left in d_left.items(): + for var_right, coeff_right in d_right.items(): + d[var_left*var_right] += coeff_left*coeff_right + return dict(d) + + def map_quotient(self, expr): + d_num = self.rec(expr.numerator) + d_den = self.rec(expr.denominator) + if len(d_den) > 1: + raise ValueError + den_var, den_coeff = list(d_den.items())[0] + + return {num_var/den_var: num_coeff/den_coeff for + num_var, num_coeff in d_num.items()} + + def map_power(self, expr): + d_base = self.rec(expr.base) + d_exponent = self.rec(expr.exponent) + # d_exponent should look like {1: k} + if len(d_exponent) > 1 or 1 not in d_exponent: + raise RuntimeError("nonlinear expression") + exp, = d_exponent.values() + if exp == 1: + return d_base + if len(d_base) > 1: + raise NotImplementedError("powers are not implemented") + (var, coeff), = d_base.items() + return {var**exp: coeff**exp} + + rec = __call__ + + +def _kernel_source_derivs_as_poly(kernel, axis_vars): + """Converts a :class:`sumpy.kernel.Kernel` object to a polynomial. + A :class:`sumpy.kernel.Kernel` represents a derivative operator + and the derivative operator is converted to a polynomial with + variables given by `axis_vars`. + + For eg: for source x the derivative operator, + d/dx_1 dx_2 + d/dx_1 is converted to x_2 * x_1 + x_1. + """ + if isinstance(kernel, AxisSourceDerivative): + poly = _kernel_source_derivs_as_poly(kernel.inner_kernel, axis_vars) + return -axis_vars[kernel.axis]*poly + if isinstance(kernel, KernelWrapper): + raise ValueError + return 1 + +# }}} + + +# {{{ factor the matrix + +def _syzygy_module_groebner_basis_mat(m, generators): + """Takes as input a module of polynomials with domain :class:`sympy.EX` + represented as a matrix and returns the syzygy module as a matrix of polynomials + in the same domain. The syzygy module *S* that is returned as a matrix + satisfies S m = 0 and is a left nullspace of the matrix. + Using :class:`sympy.EX` because that represents the domain with any symbolic + element. Usually we need an Integer or Rational domain, but since there can be + unrelated symbols like *mu* in the expression, we need to use a symbolic domain. + """ + from sympy.polys.orderings import grevlex + + def _convert_to_matrix(module, *generators): + result = [] + for syzygy in module: + row = [] + for dmp in syzygy.data: + row.append(sympy.Poly(dmp.to_dict(), *generators, + domain=sympy.EX).as_expr()) + result.append(row) + return sympy.Matrix(result) + + ring = sympy.EX.old_poly_ring(*generators, order=grevlex) + column_ideals = [ring.free_module(1).submodule(*m[:, i].tolist(), order=grevlex) + for i in range(m.shape[1])] + column_syzygy_modules = [ideal.syzygy_module() for ideal in column_ideals] + + intersection = functools.reduce(lambda x, y: x.intersect(y), + column_syzygy_modules) + + # _groebner_vec returns a groebner basis of the syzygy module as a list + groebner_vec = intersection._groebner_vec() + return _convert_to_matrix(groebner_vec, *generators) + + +def _factor_left(mat, axis_vars): + """Return the left hand side of the factorisation of the matrix + For a matrix M, we want to find a factorisation such that M = L R + with minimum number of columns of L. The polynomials represent + derivative operators and therefore division is not well defined. + To avoid divisions, we work in a polynomial ring which doesn't + have division either. + + To get a good factorisation, what we do is first find a matrix + such that S M = 0 where S is the syzygy module converted to a matrix. + It can also be referred to as the left nullspace of the matrix. + Then, M.T S.T = 0 which implies that M.T is in the space spanned by + the syzygy module of S.T and to get M we get the transpose of that. + """ + syzygy_of_mat = _syzygy_module_groebner_basis_mat(mat, axis_vars) + if len(syzygy_of_mat) == 0: + raise ValueError("could not find a factorization") + return _syzygy_module_groebner_basis_mat(syzygy_of_mat.T, axis_vars).T + + +def _factor_right(mat, factor_left): + """Return the right hand side of the factorisation of the matrix + Note that from the construction of *factor_left*, we know that + the factor on the right has elements in the same polynomial ring + as the input matrix *mat*. Therefore, doing divisions are fine + as they should result in exact polynomial divisions and will have + no remainders. + """ + return factor_left.LUsolve(sympy.Matrix(mat), + iszerofunc=lambda x: x.simplify() == 0) + +# }}} + + +# {{{ convert factors back into IntGs + +def _convert_source_poly_to_int_g_derivs(poly, orig_int_g, axis_vars): + """This does the opposite of :func:`_kernel_source_derivs_as_poly` + and converts a polynomial back to a source derivative + operator. First it is converted to a :class:`sumpy.kernel.Kernel` + and then to a :class:`~pytential.symbolic.primitives.IntG`. + """ + from pytential.symbolic.pde.system_utils import simplify_densities + to_pymbolic = SympyToPymbolicMapper() + + orig_kernel = orig_int_g.source_kernels[0] + source_kernels = [] + densities = [] + for monom, coeff in poly.terms(): + kernel = orig_kernel + for idim, rep in enumerate(monom): + for _ in range(rep): + kernel = AxisSourceDerivative(idim, kernel) + source_kernels.append(kernel) + # (-1) below is because d/dx f(c - x) = - f'(c - x) + densities.append(to_pymbolic(coeff) * (-1)**sum(monom)) + return orig_int_g.copy(source_kernels=tuple(source_kernels), + densities=tuple(simplify_densities(densities))) + + +def _convert_target_poly_to_int_g_derivs(poly, orig_int_g, rhs_int_g): + """This does the opposite of :func:`_kernel_source_derivs_as_poly` + and converts a polynomial back to a target derivative + operator. It is applied to a :class:`~pytential.symbolic.primitives.IntG` + object and returns a new instance. + """ + to_pymbolic = SympyToPymbolicMapper() + + result = 0 + for monom, coeff in poly.terms(): + kernel = orig_int_g.target_kernel + for idim, rep in enumerate(monom): + for _ in range(rep): + kernel = AxisTargetDerivative(idim, kernel) + result += orig_int_g.copy(target_kernel=kernel, + source_kernels=rhs_int_g.source_kernels, + densities=rhs_int_g.densities) * to_pymbolic(coeff) + + return result + +# }}} + +# vim: fdm=marker diff --git a/pytential/symbolic/pde/system_utils.py b/pytential/symbolic/pde/system_utils.py new file mode 100644 index 000000000..8891db22d --- /dev/null +++ b/pytential/symbolic/pde/system_utils.py @@ -0,0 +1,950 @@ +__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, + DirectionalSourceDerivative, ExpressionKernel, + KernelWrapper, TargetPointMultiplier, DirectionalDerivative) +from pytools import (memoize_on_first_arg, + generate_nonnegative_integer_tuples_summing_to_at_most + as gnitstam) +from collections import defaultdict + +from pymbolic.geometric_algebra.mapper import WalkMapper +from pymbolic.mapper import CombineMapper +from pymbolic.mapper.coefficient import CoefficientCollector +from pytential.symbolic.primitives import (IntG, NodeCoordinateComponent, + hashable_kernel_args, hashable_kernel_arg_value) +from pytential.symbolic.mappers import IdentityMapper +from pytential.utils import chop, lu_solve_with_expand +import pytential + +from pytential.symbolic.pde.reduce_fmms import reduce_number_of_fmms + +import logging +logger = logging.getLogger(__name__) + +__all__ = ( + "merge_int_g_exprs", + "get_deriv_relation", + ) + +__doc__ = """ +.. autofunction:: rewrite_using_base_kernel +.. autofunction:: merge_int_g_exprs +.. 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)] + f = sympy.Function("f") + orig_expr = f(*ds) + 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} + +# }}} + + +def merge_int_g_exprs(exprs, source_dependent_variables=None): + """ + Merge expressions involving :class:`~pytential.symbolic.primitives.IntG` + objects. + + Several techniques are used for merging and reducing number of FMMs + + * :class:`sumpy.kernel.AxisTargetDerivative` instances are converted + to :class:`sumpy.kernel.AxisSourceDerivative` instances. + (by flipping signs, assuming translation-invariance). + Target derivatives will be brought back by the syzygy module + construction below if beneficial. + (For example, `D + d/dx(S)` can be re-written as `D - d/dy(S)` which can be + done in one FMM) + + * If there is a sum of two *IntG* s with same target derivative and different + source derivatives of the same kernel, they are merged into one FMM. + + * Reduce the number of FMMs by converting the *IntG* expression to + a matrix and factoring the matrix where the left operand matrix represents + a transformation at target and the right matrix represents a transformation + at source. For this to work, we need to know which variables depend on + source so that they do not end up in the left operand. User needs to supply + this as the argument *source_dependent_variable*. This is done by the + call to :func:`pytential.symbolic.pde.reduce_fmms.reduce_number_of_fmms`. + + :arg base_kernel: A :class:`sumpy.kernel.Kernel` object if given will be used + for converting a :class:`~pytential.symbolic.primitives.IntG` to a linear + expression of same type with the kernel replaced by base_kernel and its + derivatives + + :arg source_dependent_variable: When merging expressions, consider only these + variables as dependent on source. This is important when reducing the + number of FMMs needed for the output. + """ + # Using a dictionary instead of a set because sets are unordered + all_source_group_identifiers = {} + + result = np.array([0 for _ in exprs], dtype=object) + + int_g_cc = IntGCoefficientCollector() + int_gs_by_source_group = defaultdict(list) + + def add_int_gs_in_expr(expr): + for int_g in get_int_g_s([expr]): + source_group_identifier = get_int_g_source_group_identifier(int_g) + int_gs_by_source_group[source_group_identifier].append(int_g) + for density in int_g.densities: + add_int_gs_in_expr(density) + + for i, expr in enumerate(exprs): + int_gs_by_group = {} + try: + int_g_coeff_map = int_g_cc(expr) + except (RuntimeError, AssertionError): + # Don't touch this expression, because it's not linear. + # FIXME: if there's ever any use case, then we can extract + # some IntGs from them. + logger.debug("%s is not linear", expr) + result[i] += expr + add_int_gs_in_expr(expr) + continue + for int_g, coeff in int_g_coeff_map.items(): + if int_g == 1: + # coeff map may have some constant terms, add them to + result[i] += coeff + continue + + # convert DirectionalSourceDerivative to AxisSourceDerivative + # as kernel arguments need to be the same for merging + int_g = convert_directional_source_to_axis_source(int_g) + # convert TargetDerivative to source before checking the group + # as the target kernel has to be the same for merging + int_g = convert_target_deriv_to_source(int_g) + if not is_expr_target_dependent(coeff): + # move the coefficient inside + int_g = int_g.copy(densities=[density*coeff for density in + int_g.densities]) + coeff = 1 + + source_group_identifier = get_int_g_source_group_identifier(int_g) + target_group_identifier = get_int_g_target_group_identifier(int_g) + group = (source_group_identifier, target_group_identifier, coeff) + + all_source_group_identifiers[source_group_identifier] = 1 + + if group not in int_gs_by_group: + new_int_g = int_g + else: + prev_int_g = int_gs_by_group[group] + # Let's merge IntGs with the same group + new_int_g = merge_two_int_gs(int_g, prev_int_g) + int_gs_by_group[group] = new_int_g + + # Do some simplifications after merging. Not stricty necessary + for (_, _, coeff), int_g in int_gs_by_group.items(): + # replace an IntG with d axis source derivatives to an IntG + # with one directional source derivative + # TODO: reenable this later + # result_int_g = convert_axis_source_to_directional_source(int_g) + # simplify the densities as they may become large due to pymbolic + # not doing automatic simplifications unlike sympy/symengine + result_int_g = int_g.copy( + densities=simplify_densities(int_g.densities)) + result[i] += result_int_g * coeff + add_int_gs_in_expr(result_int_g) + + # No IntGs found + if all(not int_gs for int_gs in int_gs_by_source_group): + return exprs + + # Do the calculation for each source_group_identifier separately + # and assemble them + replacements = {} + for int_gs in int_gs_by_source_group.values(): + # For each output, we now have a sum of int_gs with + # different target attributes. + # for eg: {+}S + {-}D (where {x} is the QBX limit). + # We can't merge them together, because merging implies + # that everything happens at the source level and therefore + # require same target attributes. + # + # To handle this case, we can treat them separately as in + # different source base kernels, but that would imply more + # FMMs than necessary. + # + # Take the following example, + # + # [ {+}(S + D), {-}S + {avg}D, {avg}S + {-}D] + # + # If we treated the target attributes separately, then we + # will be reducing [{+}(S + D), 0, 0], [0, {-}S, {-}D], + # [0, {avg}D, {avg}S] separately which results in + # [{+}(S + D)], [{-}S, {-}D], [{avg}S, {avg}D] as + # the reduced FMMs and pytential will calculate + # [S + D, S, D] as three separate FMMs and then assemble + # the three outputs by applying target attributes. + # + # Instead, we can do S, D as two separate FMMs and get the + # result for all three outputs. To do that, we will first + # get all five expressions in the example + # [ {+}(S + D), {-}S, {avg}D, {avg}S, {-}D] + # and then remove the target attributes to get, + # [S + D, S, D]. We will reduce these and restore the target + # attributes at the end + + targetless_int_g_mapping = defaultdict(list) + for int_g in int_gs: + common_int_g = remove_target_attributes(int_g) + targetless_int_g_mapping[common_int_g].append(int_g) + + insns_to_reduce = list(targetless_int_g_mapping.keys()) + reduced_insns = reduce_number_of_fmms(insns_to_reduce, + source_dependent_variables) + + for insn, reduced_insn in zip(insns_to_reduce, reduced_insns): + for int_g in targetless_int_g_mapping[insn]: + replacements[int_g] = restore_target_attributes(reduced_insn, int_g) + + mapper = IntGSubstitutor(replacements) + result = [mapper(expr) for expr in result] + + orig_count = get_number_of_fmms(exprs) + new_count = get_number_of_fmms(result) + if orig_count < new_count: + raise RuntimeError("merge_int_g_exprs failed. " + "Please open an issue in pytential bug tracker.") + + return result + + +def get_number_of_fmms(exprs): + fmms = set() + for int_g in get_int_g_s(exprs): + fmms.add(remove_target_attributes(int_g)) + return len(fmms) + + +class IntGCoefficientCollector(CoefficientCollector): + def __init__(self): + super().__init__({}) + + def map_int_g(self, expr): + return {expr: 1} + + def map_algebraic_leaf(self, expr, *args, **kwargs): + return {1: expr} + + handle_unsupported_expression = map_algebraic_leaf + + +def is_expr_target_dependent(expr): + mapper = IsExprTargetDependent() + return mapper(expr) + + +class IsExprTargetDependent(CombineMapper): + def combine(self, values): + import operator + from functools import reduce + return reduce(operator.or_, values, False) + + def map_constant(self, expr): + return False + + map_variable = map_constant + map_wildcard = map_constant + map_function_symbol = map_constant + + def map_common_subexpression(self, expr): + return self.rec(expr.child) + + def map_coordinate_component(self, expr): + return True + + def map_num_reference_derivative(self, expr): + return True + + def map_q_weight(self, expr): + return True + + +def get_hashable_kernel_argument(arg): + if hasattr(arg, "__iter__"): + try: + return tuple(arg) + except TypeError: + pass + return arg + + +def get_int_g_source_group_identifier(int_g): + """Return a identifier for a group for the *int_g* so that all elements in that + group have the same source attributes. + """ + target_arg_names = get_normal_vector_names(int_g.target_kernel) + args = dict((k, v) for k, v in sorted( + int_g.kernel_arguments.items()) if k not in target_arg_names) + return (int_g.source, hashable_kernel_args(args), + int_g.target_kernel.get_base_kernel()) + + +def get_int_g_target_group_identifier(int_g): + """Return a identifier for a group for the *int_g* so that all elements in that + group have the same target attributes. + """ + target_arg_names = get_normal_vector_names(int_g.target_kernel) + args = dict((k, v) for k, v in sorted( + int_g.kernel_arguments.items()) if k in target_arg_names) + return (int_g.target, int_g.qbx_forced_limit, int_g.target_kernel, + hashable_kernel_args(args)) + + +def get_normal_vector_names(kernel): + """Return the normal vector names in a kernel + """ + normal_vectors = set() + while isinstance(kernel, KernelWrapper): + if isinstance(kernel, DirectionalDerivative): + normal_vectors.add(kernel.dir_vec_name) + kernel = kernel.inner_kernel + return normal_vectors + + +class IntGSubstitutor(IdentityMapper): + """Replaces IntGs with pymbolic expression given by the + replacements dictionary + """ + def __init__(self, replacements): + self.replacements = replacements + + def map_int_g(self, expr): + if expr in self.replacements: + new_expr = self.replacements[expr] + if new_expr != expr: + return self.rec(new_expr) + else: + expr = new_expr + + densities = [self.rec(density) for density in expr.densities] + return expr.copy(densities=tuple(densities)) + + +def remove_target_attributes(int_g): + """Remove target attributes from *int_g* and return an expression + that is common to all expression in the same source group. + """ + normals = get_normal_vector_names(int_g.target_kernel) + kernel_arguments = {k: v for k, v in int_g.kernel_arguments.items() if + k not in normals} + return int_g.copy(target=None, qbx_forced_limit=None, + target_kernel=int_g.target_kernel.get_base_kernel(), + kernel_arguments=kernel_arguments) + + +def restore_target_attributes(expr, orig_int_g): + """Restore target attributes from *orig_int_g* to all the + :class:`~pytential.symbolic.primitives.IntG` objects in the + input *expr*. + """ + int_gs = get_int_g_s([expr]) + + replacements = { + int_g: int_g.copy(target=orig_int_g.target, + qbx_forced_limit=orig_int_g.qbx_forced_limit, + target_kernel=orig_int_g.target_kernel.replace_base_kernel( + int_g.target_kernel), + kernel_arguments=orig_int_g.kernel_arguments) + for int_g in int_gs} + + substitutor = IntGSubstitutor(replacements) + return substitutor(expr) + + +def merge_two_int_gs(int_g_1, int_g_2): + kernel_arguments = merge_kernel_arguments(int_g_1.kernel_arguments, + int_g_2.kernel_arguments) + source_kernels = int_g_1.source_kernels + int_g_2.source_kernels + densities = int_g_1.densities + int_g_2.densities + + return int_g_1.copy( + source_kernels=tuple(source_kernels), + densities=tuple(densities), + kernel_arguments=kernel_arguments, + ) + + +class GetIntGs(WalkMapper): + """A Mapper that walks expressions and collects + :class:`~pytential.symbolic.primitives.IntG` objects + """ + def __init__(self): + self.int_g_s = set() + + def map_int_g(self, expr): + self.int_g_s.add(expr) + + def map_constant(self, expr): + pass + + map_variable = map_constant + handle_unsupported_expression = map_constant + + +def have_int_g_s(expr): + """Checks if a :mod:`pymbolic` expression has + :class:`~pytential.symbolic.primitives.IntG` objects + """ + mapper = GetIntGs() + mapper(expr) + return bool(mapper.int_g_s) + + +def get_int_g_s(exprs): + """Returns all :class:`~pytential.symbolic.primitives.IntG` objects + in a list of :mod:`pymbolic` expressions. + """ + get_int_g_mapper = GetIntGs() + [get_int_g_mapper(expr) for expr in exprs] + return get_int_g_mapper.int_g_s + + +def convert_directional_source_to_axis_source(int_g): + """Convert an IntG with a DirectionalSourceDerivative instance + to an IntG with d AxisSourceDerivative instances. + """ + source_kernels = [] + densities = [] + for source_kernel, density in zip(int_g.source_kernels, int_g.densities): + knl_result = _convert_directional_source_knl_to_axis_source(source_kernel, + int_g.kernel_arguments) + for knl, coeff in knl_result: + source_kernels.append(knl) + densities.append(coeff * density) + + kernel_arguments = filter_kernel_arguments( + list(source_kernels) + [int_g.target_kernel], int_g.kernel_arguments) + return int_g.copy(source_kernels=tuple(source_kernels), + densities=tuple(densities), kernel_arguments=kernel_arguments) + + +def _convert_directional_source_knl_to_axis_source(knl, knl_arguments): + if isinstance(knl, DirectionalSourceDerivative): + dim = knl.dim + dir_vec = knl_arguments[knl.dir_vec_name] + + res = [] + inner_result = _convert_directional_source_knl_to_axis_source( + knl.inner_kernel, knl_arguments) + for inner_knl, coeff in inner_result: + for d in range(dim): + res.append((AxisSourceDerivative(d, inner_knl), coeff*dir_vec[d])) + return res + elif isinstance(knl, KernelWrapper): + inner_result = _convert_directional_source_knl_to_axis_source( + knl.inner_kernel, knl_arguments) + return [(knl.replace_inner_kernel(inner_knl), coeff) for + inner_knl, coeff in inner_result] + else: + return [(knl, 1)] + + +def convert_axis_source_to_directional_source(int_g): + """Convert an IntG with d AxisSourceDerivative instances to + an IntG with one DirectionalSourceDerivative instance. + """ + from pytential.symbolic.primitives import _DIR_VEC_NAME + if not isinstance(int_g, IntG): + return int_g + knls = list(int_g.source_kernels) + dim = knls[0].dim + if len(knls) != dim: + return int_g + if not any(isinstance(knl, AxisSourceDerivative) for knl in knls): + return int_g + # TODO: sort + axes = [knl.axis for knl in knls] + if axes != list(range(dim)): + return int_g + base_knls = set(knl.inner_kernel for knl in knls) + if len(base_knls) > 1: + return int_g + base_knl = base_knls.pop() + kernel_arguments = int_g.kernel_arguments.copy() + + name = _DIR_VEC_NAME + count = 0 + while name in kernel_arguments: + name = _DIR_VEC_NAME + f"_{count}" + count += 1 + + kernel_arguments[name] = \ + np.array(int_g.densities, dtype=np.object) + res = int_g.copy( + source_kernels=( + DirectionalSourceDerivative(base_knl, dir_vec_name=name),), + densities=(1,), + kernel_arguments=kernel_arguments) + return res + + +def merge_kernel_arguments(x, y): + """merge two kernel argument dictionaries and raise a ValueError if + the two dictionaries do not agree for duplicate keys. + """ + res = x.copy() + for k, v in y.items(): + if k in res: + if hashable_kernel_arg_value(res[k]) \ + != hashable_kernel_arg_value(v): + raise ValueError(f"Error merging values for {k}." + f"values were {res[k]} and {v}") + else: + res[k] = v + return res + + +def simplify_densities(densities): + """Simplify densities by converting to sympy and converting back + to trigger sympy's automatic simplification routines. + """ + from sumpy.symbolic import (SympyToPymbolicMapper, PymbolicToSympyMapper) + from pymbolic.mapper import UnsupportedExpressionError + to_sympy = PymbolicToSympyMapper() + to_pymbolic = SympyToPymbolicMapper() + result = [] + for density in densities: + try: + result.append(to_pymbolic(to_sympy(density))) + except (ValueError, NotImplementedError, UnsupportedExpressionError): + logger.debug("%s cannot be simplified", density) + result.append(density) + return tuple(result) + + +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/primitives.py b/pytential/symbolic/primitives.py index a4a652c71..8c071faa5 100644 --- a/pytential/symbolic/primitives.py +++ b/pytential/symbolic/primitives.py @@ -23,6 +23,7 @@ from sys import intern from warnings import warn from functools import partial +from collections import OrderedDict import numpy as np @@ -1260,14 +1261,15 @@ def laplace(ambient_dim, operand): # {{{ potentials -def hashable_kernel_args(kernel_arguments): - hashable_args = [] - for key, val in sorted(kernel_arguments.items()): - if isinstance(val, np.ndarray): - val = tuple(val) - hashable_args.append((key, val)) +def hashable_kernel_arg_value(val): + if isinstance(val, np.ndarray): + val = tuple(val) + return val + - return tuple(hashable_args) +def hashable_kernel_args(kernel_arguments): + return tuple([(key, hashable_kernel_arg_value(val)) for key, val in + sorted(kernel_arguments.items())]) class IntG(Expression): @@ -1352,8 +1354,18 @@ def __init__(self, target_kernel, source_kernels, densities, raise ValueError("invalid value (%s) of qbx_forced_limit" % qbx_forced_limit) - source_kernels = tuple(source_kernels) - densities = tuple(densities) + # Fold duplicates in source_kernels + knl_density_dict = OrderedDict() + for density, source_kernel in zip(densities, source_kernels): + if source_kernel in knl_density_dict: + knl_density_dict[source_kernel] += density + else: + knl_density_dict[source_kernel] = density + knl_density_dict = OrderedDict( + [(k, v) for k, v in knl_density_dict.items() if v]) + densities = tuple(knl_density_dict.values()) + source_kernels = tuple(knl_density_dict.keys()) + kernel_arg_names = set() for kernel in source_kernels + (target_kernel,): @@ -1378,6 +1390,7 @@ def __init__(self, target_kernel, source_kernels, densities, provided_arg_names = set(kernel_arguments.keys()) missing_args = kernel_arg_names - provided_arg_names + if missing_args: raise TypeError("kernel argument(s) '%s' not supplied" % ", ".join(missing_args)) 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/unregularized.py b/pytential/unregularized.py index 7c46697bb..9c6b53516 100644 --- a/pytential/unregularized.py +++ b/pytential/unregularized.py @@ -119,10 +119,11 @@ def evaluate_wrapper(expr): def op_group_features(self, expr): from pytential.utils import sort_arrays_together + from sumpy.kernel import TargetTransformationRemover result = ( expr.source, *sort_arrays_together(expr.source_kernels, expr.densities, key=str), - expr.target_kernel, + TargetTransformationRemover()(expr.target_kernel), ) return result 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..ae460c7d1 100644 --- a/requirements.txt +++ b/requirements.txt @@ -3,7 +3,7 @@ numpy != 1.22.0 git+https://github.com/inducer/pytools.git#egg=pytools git+https://github.com/inducer/pymbolic.git#egg=pymbolic -sympy +sympy>=1.10 git+https://github.com/inducer/modepy.git#egg=modepy git+https://github.com/inducer/pyopencl.git#egg=pyopencl git+https://github.com/inducer/islpy.git#egg=islpy @@ -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_pde_system_utils.py b/test/test_pde_system_utils.py new file mode 100644 index 000000000..d78176f00 --- /dev/null +++ b/test/test_pde_system_utils.py @@ -0,0 +1,320 @@ +__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. +""" + +from sumpy.kernel import (LaplaceKernel, AxisSourceDerivative, + AxisTargetDerivative, TargetPointMultiplier, BiharmonicKernel, HelmholtzKernel, + DirectionalTargetDerivative) +from pytential.symbolic.primitives import (int_g_vec, D, IntG, + NodeCoordinateComponent) +from pytential.symbolic.pde.system_utils import (merge_int_g_exprs, + rewrite_using_base_kernel) +from pymbolic.primitives import make_sym_vector, Variable + + +def test_reduce_number_of_fmms(): + dim = 3 + knl = LaplaceKernel(dim) + densities = make_sym_vector("sigma", 2) + mu = Variable("mu") + + int_g1 = \ + mu * int_g_vec(AxisSourceDerivative(1, AxisSourceDerivative(0, knl)), + densities[0], qbx_forced_limit=1) + \ + int_g_vec(AxisSourceDerivative(1, AxisSourceDerivative(1, knl)), + densities[1] * mu, qbx_forced_limit=1) + + int_g2 = \ + int_g_vec(AxisSourceDerivative(2, AxisSourceDerivative(0, knl)), + densities[0], qbx_forced_limit=1) + \ + int_g_vec(AxisSourceDerivative(2, AxisSourceDerivative(1, knl)), + densities[1], qbx_forced_limit=1) + + # Merging reduces 4 FMMs to 2 FMMs and then further reduced to 1 FMM + result = merge_int_g_exprs([int_g1, int_g2], source_dependent_variables=[]) + + int_g3 = \ + IntG(target_kernel=AxisTargetDerivative(1, knl), + source_kernels=[AxisSourceDerivative(0, knl), + AxisSourceDerivative(1, knl)], + densities=[-mu * densities[0], -mu * densities[1]], + qbx_forced_limit=1) + + int_g4 = \ + IntG(target_kernel=AxisTargetDerivative(2, knl), + source_kernels=[AxisSourceDerivative(0, knl), + AxisSourceDerivative(1, knl)], + densities=[-mu * densities[0], -mu * densities[1]], + qbx_forced_limit=1) + + assert result[0] == int_g3 + assert result[1] == int_g4 * mu**(-1) + + +def test_source_dependent_variable(): + # Same example as test_reduce_number_of_fmms, but with + # mu marked as a source dependent variable + dim = 3 + knl = LaplaceKernel(dim) + densities = make_sym_vector("sigma", 2) + mu = Variable("mu") + nu = Variable("nu") + + int_g1 = \ + int_g_vec(AxisSourceDerivative(1, AxisSourceDerivative(0, knl)), + mu * nu * densities[0], qbx_forced_limit=1) + \ + int_g_vec(AxisSourceDerivative(1, AxisSourceDerivative(1, knl)), + mu * nu * densities[1], qbx_forced_limit=1) + + int_g2 = \ + int_g_vec(AxisSourceDerivative(2, AxisSourceDerivative(0, knl)), + densities[0], qbx_forced_limit=1) + \ + int_g_vec(AxisSourceDerivative(2, AxisSourceDerivative(1, knl)), + densities[1], qbx_forced_limit=1) + + result = merge_int_g_exprs([int_g1, int_g2], + source_dependent_variables=[mu, nu]) + + # Merging reduces 4 FMMs to 2 FMMs. No further reduction of FMMs. + int_g3 = \ + IntG(target_kernel=knl, + source_kernels=[AxisSourceDerivative(1, AxisSourceDerivative(1, knl)), + AxisSourceDerivative(1, AxisSourceDerivative(0, knl))], + densities=[mu * nu * densities[1], mu * nu * densities[0]], + qbx_forced_limit=1) + + int_g4 = \ + IntG(target_kernel=knl, + source_kernels=[AxisSourceDerivative(2, AxisSourceDerivative(1, knl)), + AxisSourceDerivative(2, AxisSourceDerivative(0, knl))], + densities=[densities[1], densities[0]], + qbx_forced_limit=1) + + assert result[0] == int_g3 + assert result[1] == int_g4 + + +def test_base_kernel_merge(): + # Same example as test_reduce_number_of_fmms, but with + # mu marked as a source dependent variable + dim = 3 + knl = LaplaceKernel(dim) + biharm_knl = BiharmonicKernel(dim) + density = make_sym_vector("sigma", 1)[0] + + int_g1 = \ + int_g_vec(TargetPointMultiplier(0, knl), + density, qbx_forced_limit=1) + + int_g2 = \ + int_g_vec(TargetPointMultiplier(1, knl), + density, qbx_forced_limit=1) + + exprs_rewritten = rewrite_using_base_kernel([int_g1, int_g2], + base_kernel=biharm_knl) + result = merge_int_g_exprs(exprs_rewritten, source_dependent_variables=[]) + + sources = [NodeCoordinateComponent(i) for i in range(dim)] + + source_kernels = list(reversed([ + AxisSourceDerivative(i, AxisSourceDerivative(i, biharm_knl)) + for i in range(dim)])) + + int_g3 = IntG(target_kernel=biharm_knl, + source_kernels=source_kernels + [AxisSourceDerivative(0, biharm_knl)], + densities=[density*sources[0]*(-1.0) for _ in range(dim)] + [2*density], + qbx_forced_limit=1) + int_g4 = IntG(target_kernel=biharm_knl, + source_kernels=source_kernels + [AxisSourceDerivative(1, biharm_knl)], + densities=[density*sources[1]*(-1.0) for _ in range(dim)] + [2*density], + qbx_forced_limit=1) + + assert result[0] == int_g3 + assert result[1] == int_g4 + + +def test_merge_different_kernels(): + # Test different kernels Laplace, Helmholtz(k=1), Helmholtz(k=2) + dim = 3 + laplace_knl = LaplaceKernel(dim) + helmholtz_knl = HelmholtzKernel(dim) + density = make_sym_vector("sigma", 1)[0] + + int_g1 = int_g_vec(laplace_knl, density, qbx_forced_limit=1) \ + + int_g_vec(helmholtz_knl, density, qbx_forced_limit=1, k=1) \ + + int_g_vec(AxisTargetDerivative(0, helmholtz_knl), + density, qbx_forced_limit=1, k=1) \ + + int_g_vec(helmholtz_knl, density, qbx_forced_limit=1, k=2) + + int_g2 = int_g_vec(AxisTargetDerivative(0, laplace_knl), + density, qbx_forced_limit=1) + + result = merge_int_g_exprs([int_g1, int_g2], + source_dependent_variables=[]) + + int_g3 = int_g_vec(laplace_knl, density, qbx_forced_limit=1) \ + + IntG(target_kernel=helmholtz_knl, + source_kernels=[AxisSourceDerivative(0, helmholtz_knl), helmholtz_knl], + densities=[-density, density], + qbx_forced_limit=1, k=1) \ + + int_g_vec(helmholtz_knl, density, qbx_forced_limit=1, k=2) + + assert result[0] == int_g3 + assert result[1] == int_g2 + + +def test_merge_different_qbx_forced_limit(): + dim = 3 + laplace_knl = LaplaceKernel(dim) + density = make_sym_vector("sigma", 1)[0] + + int_g1 = int_g_vec(laplace_knl, density, qbx_forced_limit=1) + int_g2 = int_g1.copy(target_kernel=AxisTargetDerivative(0, laplace_knl)) + + int_g3, = merge_int_g_exprs([int_g2 + int_g1]) + int_g4 = int_g1.copy(qbx_forced_limit=2) + int_g2.copy(qbx_forced_limit=-2) + int_g5 = int_g1.copy(qbx_forced_limit=-2) + int_g2.copy(qbx_forced_limit=2) + + result = merge_int_g_exprs([int_g3, int_g4, int_g5], + source_dependent_variables=[]) + + int_g6 = int_g_vec(laplace_knl, -density, qbx_forced_limit=1) + int_g7 = int_g6.copy(target_kernel=AxisTargetDerivative(0, laplace_knl)) + int_g8 = int_g7 * (-1) + int_g6 * (-1) + int_g9 = int_g6.copy(qbx_forced_limit=2) * (-1) \ + + int_g7.copy(qbx_forced_limit=-2) * (-1) + int_g10 = int_g6.copy(qbx_forced_limit=-2) * (-1) \ + + int_g7.copy(qbx_forced_limit=2) * (-1) + + assert result[0] == int_g8 + assert result[1] == int_g9 + assert result[2] == int_g10 + + +def test_merge_directional_source(): + from pymbolic.primitives import Variable + from pytential.symbolic.primitives import cse + + dim = 3 + laplace_knl = LaplaceKernel(dim) + density = Variable("density") + + int_g1 = int_g_vec(laplace_knl, density, qbx_forced_limit=1) + int_g2 = D(laplace_knl, density, qbx_forced_limit=1) + + source_kernels = [AxisSourceDerivative(d, laplace_knl) + for d in range(dim)] + [laplace_knl] + dsource = int_g2.kernel_arguments["dsource_vec"] + densities = [dsource[d]*cse(density) for d in range(dim)] + [density] + int_g3 = int_g2.copy(source_kernels=source_kernels, densities=densities, + kernel_arguments={}) + + result = merge_int_g_exprs([int_g1 + int_g2], + source_dependent_variables=[density]) + assert result[0] == int_g3 + + result = merge_int_g_exprs([int_g1 + int_g2]) + assert result[0] == int_g3 + + +def test_merge_directional_target(): + from pymbolic.primitives import Variable + + dim = 3 + knl = DirectionalTargetDerivative(LaplaceKernel(dim), "target_dir") + density = Variable("density") + target_dir1 = make_sym_vector("target_dir1", dim) + target_dir2 = make_sym_vector("target_dir2", dim) + + int_g1 = int_g_vec(knl, density, qbx_forced_limit=1, target_dir=target_dir1) + int_g2 = int_g_vec(knl, density, qbx_forced_limit=1, target_dir=target_dir2) + + result = merge_int_g_exprs([int_g1 + int_g2]) + assert result[0] == int_g1 + int_g2 + + +def test_restoring_target_attributes(): + from pymbolic.primitives import Variable + dim = 3 + laplace_knl = LaplaceKernel(dim) + density = Variable("density") + + int_g1 = int_g_vec(TargetPointMultiplier(0, AxisTargetDerivative(0, + laplace_knl)), density, qbx_forced_limit=1) + int_g2 = int_g_vec(AxisTargetDerivative(1, laplace_knl), + density, qbx_forced_limit=1) + + result = merge_int_g_exprs([int_g1, int_g2], + source_dependent_variables=[]) + + assert result[0] == int_g1 + assert result[1] == int_g2 + + +def test_int_gs_in_densities(): + from pymbolic.primitives import Variable, Quotient + dim = 3 + laplace_knl = LaplaceKernel(dim) + density = Variable("density") + + int_g1 = \ + int_g_vec(laplace_knl, + int_g_vec(AxisSourceDerivative(2, laplace_knl), density, + qbx_forced_limit=1), qbx_forced_limit=1) + \ + int_g_vec(AxisTargetDerivative(0, laplace_knl), + int_g_vec(AxisSourceDerivative(1, laplace_knl), 2*density, + qbx_forced_limit=1), qbx_forced_limit=1) + + # In the above example the two inner source derivatives should + # be converted to target derivatives and the two outermost + # IntGs should be merged into one by converting the target + # derivative in the last term to a source derivative + result = merge_int_g_exprs([int_g1]) + + source_kernels = [AxisSourceDerivative(0, laplace_knl), laplace_knl] + densities = [ + (-1)*int_g_vec(AxisTargetDerivative(1, laplace_knl), + (-2)*density, qbx_forced_limit=1), + int_g_vec(AxisTargetDerivative(2, laplace_knl), + (-2)*density, qbx_forced_limit=1) * Quotient(1, 2) + ] + int_g3 = IntG(target_kernel=laplace_knl, + source_kernels=tuple(source_kernels), + densities=tuple(densities), + qbx_forced_limit=1) + + assert result[0] == int_g3 + + +# You can test individual routines by typing +# $ python test_pde_system_tools.py 'test_reduce_number_of_fmms()' + +if __name__ == "__main__": + import logging + logging.basicConfig(level=logging.DEBUG) + import sys + if len(sys.argv) > 1: + exec(sys.argv[1]) + else: + from pytest import main + main([__file__]) + +# vim: fdm=marker diff --git a/test/test_stokes.py b/test/test_stokes.py index dca8a4dbe..5d342c5aa 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), + (3, "biharmonic", 0.5), + (3, "laplace", 0.5), + + (2, "biharmonic", 0.4), + (3, "biharmonic", 0.4), + (3, "laplace", 0.4), ]) -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 += 2 + 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