Stokes biharmonic

pytential/symbolic/stokes.py

@@ -23,7 +23,8 @@
 import numpy as np
 
 from pytential import sym
-from sumpy.kernel import StokesletKernel, StressletKernel, LaplaceKernel
+from sumpy.kernel import (StokesletKernel, StressletKernel, LaplaceKernel,
+    AxisTargetDerivative, BiharmonicKernel)
 
 
 class StokesletWrapper:
@@ -55,28 +56,79 @@ class StokesletWrapper:
        for the same kernel in a different ordering.
     """
 
-    def __init__(self, dim=None):
-
+    def __init__(self, dim=None, use_biharmonic=True):
+        self.use_biharmonic = use_biharmonic
         self.dim = dim
-        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)
-                               }
+        if not (dim == 3 or dim == 2):
+            raise ValueError("unsupported dimension given to StokesletWrapper")
 
-        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)
-                               }
+        self.kernel_dict = {}
 
-        else:
-            raise ValueError("unsupported dimension given to StokesletWrapper")
+        self.base_kernel = BiharmonicKernel(dim=dim)
+
+        if not use_biharmonic:
+            for i in range(dim):
+                for j in range(i, dim):
+                    self.kernel_dict[(i, j)] = StokesletKernel(dim=dim, icomp=i,
+                                                               jcomp=j)
+
+            for i in range(dim):
+                for j in range(i):
+                    self.kernel_dict[(i, j)] = self.kernel_dict[(j, i)]
+
+    def _get_stokeslet_int_g(self, icomp, jcomp, density, mu_sym, qbx_forced_limit,
+            deriv_dirs):
+
+        """
+        Returns the Integral of the Stokeslet kernel given by `icomp` and `jcomp`
+        and its derivatives.
+        Note that,
+
+            Stokeslet(icomp, jcomp)
+                = -1/mu_sym d/(d x_icomp) d/(d x_jcomp) BiharmonicKernel
+
+        for icomp != jcomp. Also note that for 3D,
+
+            Stokeslet(i, i)
+                = -1/mu_sym sum(-d^2/(dx_j^2 BiharmonicKernel for j != i).
+
+        For 2D,
+
+            Stokeslet(i, i)
+                = -1/mu_sym (sum(d^2/dx_j^2 BiharmonicKernel for j != i)
+                             -3/(8*pi)).
+        """
+
+        if not self.use_biharmonic:
+            knl = self.kernel_dict[(icomp, jcomp)]
+            for deriv_dir in deriv_dirs:
+                knl = AxisTargetDerivative(deriv_dir, knl)
+            return sym.IntG(knl, density,
+                    qbx_forced_limit=qbx_forced_limit, mu=mu_sym)
+
+
+        def func(knl, deriv_dirs):
+            for deriv_dir in deriv_dirs:
+                knl = AxisTargetDerivative(deriv_dir, knl)
+
+            res = sym.IntG(knl, density,
+                    qbx_forced_limit=qbx_forced_limit)
+            return res
+
+        mult = -1 / mu_sym
+
+        if icomp != jcomp:
+            return mult * func(self.base_kernel,
+                deriv_dirs=(deriv_dirs + [icomp, jcomp]))
+
+        const = 0
+        if self.dim == 2 and not deriv_dirs:
+            from math import pi
+            const = -3/(8*pi)
+
+        return -mult * (sum(func(self.base_kernel,
+                deriv_dirs=(deriv_dirs + [i, i])) for i in range(self.dim)
+                if i != icomp) + const * density)
 
     def apply(self, density_vec_sym, mu_sym, qbx_forced_limit):
         """ Symbolic expressions for integrating Stokeslet kernel
@@ -93,29 +145,12 @@ def apply(self, density_vec_sym, mu_sym, qbx_forced_limit):
             for the average of the two one-sided boundary limits.
         """
 
-        sym_expr = np.empty((self.dim,), dtype=object)
+        sym_expr = np.zeros((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.int)
-            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.IntG(
-                                     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.IntG(
-                                     self.kernel_dict[ctr_key], density_vec_sym[i],
-                                     qbx_forced_limit=qbx_forced_limit, mu=mu_sym)
+                sym_expr[comp] += self._get_stokeslet_int_g(comp, i,
+                        density_vec_sym[i], mu_sym, qbx_forced_limit, deriv_dirs=[])
 
         return sym_expr
 
@@ -124,17 +159,12 @@ def apply_pressure(self, density_vec_sym, mu_sym, qbx_forced_limit):
 
         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.S(kernel, density_vec_sym[i],
-                                qbx_forced_limit=qbx_forced_limit))
-            else:
-                sym_expr = sym_expr + (DerivativeTaker(i).map_int_g(
-                                sym.S(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
 
@@ -156,36 +186,13 @@ def apply_derivative(self, deriv_dir, density_vec_sym,
             for the average of the two one-sided boundary limits.
         """
 
-        from pytential.symbolic.mappers import DerivativeTaker
-
-        sym_expr = np.empty((self.dim,), dtype=object)
+        sym_expr = self.apply(density_vec_sym, mu_sym, qbx_forced_limit)
 
         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.int)
-            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.IntG(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.IntG(self.kernel_dict[ctr_key],
-                                             density_vec_sym[i],
-                                             qbx_forced_limit=qbx_forced_limit,
-                                             mu=mu_sym))
+                sym_expr[comp] += self._get_stokeslet_int_g(comp, i,
+                        density_vec_sym[i], mu_sym, qbx_forced_limit,
+                        deriv_dirs=[deriv_dir])
 
         return sym_expr
 
@@ -220,7 +227,7 @@ def apply_stress(self, density_vec_sym, dir_vec_sym,
 
         import itertools
 
-        sym_expr = np.empty((self.dim,), dtype=object)
+        sym_expr = np.zeros((self.dim,), dtype=object)
         stresslet_obj = StressletWrapper(dim=self.dim)
 
         for comp in range(self.dim):
@@ -236,18 +243,11 @@ def apply_stress(self, density_vec_sym, dir_vec_sym,
                 var_ctr[j] += 1
                 ctr_key = tuple(var_ctr)
 
-                if i + j < 1:
-                    sym_expr[comp] = dir_vec_sym[i] * sym.IntG(
-                                     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.IntG(
-                                                stresslet_obj.kernel_dict[ctr_key],
-                                                density_vec_sym[j],
-                                                qbx_forced_limit=qbx_forced_limit,
-                                                mu=mu_sym)
+                sym_expr[comp] += dir_vec_sym[i] * sym.IntG(
+                                  stresslet_obj.kernel_dict[ctr_key],
+                                  density_vec_sym[j],
+                                  qbx_forced_limit=qbx_forced_limit,
+                                  mu=mu_sym)
 
         return sym_expr
 

test/test_stokes.py

@@ -44,7 +44,7 @@
 
 def run_exterior_stokes_2d(ctx_factory, nelements,
         mesh_order=4, target_order=4, qbx_order=4,
-        fmm_order=False,  # FIXME: FMM is slower than direct eval
+        fmm_order=7,  # FIXME: FMM is slower than direct eval
         mu=1, circle_rad=1.5, visualize=False):
 
     # This program tests an exterior Stokes flow in 2D using the
@@ -299,7 +299,6 @@ def get_obj_array(obj_array):
     return h_max, l2_err
 
 
-@pytest.mark.slowtest
 def test_exterior_stokes_2d(ctx_factory, qbx_order=3):
     from pytools.convergence import EOCRecorder
     eoc_rec = EOCRecorder()