Stokes biharmonic

doc/index.rst

@@ -23,6 +23,7 @@ Contents
     tools
     misc
     qbx
+    merge
     🚀 Github <https://github.com/inducer/pytential>
     💾 Download Releases <https://pypi.org/project/pytential>
 

doc/merge.rst

@@ -0,0 +1,4 @@
+Utilities to merge and reduce IntG expressions
+==============================================
+
+.. automodule:: pytential.symbolic.pde.system_utils

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
+
+

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

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")
 

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

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)

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)
 
 # }}}