add tests for stokes identities

test/extra_int_eq_data.py

@@ -291,6 +291,16 @@ class CircleTestCase(EllipseTestCase):
     aspect_ratio = 1.0
     radius = 1.0
 
+
+class StarfishTestCase(CurveTestCase):
+    name = "starfish"
+    narms = 5
+    amplitude = 0.25
+
+    def _curve_fn(self, t):
+        from meshmode.mesh.generation import NArmedStarfish
+        return NArmedStarfish(self.narms, self.amplitude)(t)
+
 # }}}
 
 

test/test_stokes.py

@@ -20,6 +20,7 @@
 THE SOFTWARE.
 """
 
+from functools import partial
 import pytest
 
 import numpy as np
@@ -34,6 +35,7 @@
 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__)
 
@@ -44,6 +46,16 @@
 
 # {{{ test_exterior_stokes
 
+def rnorm2(actx, x, y):
+    from meshmode.dof_array import flat_norm
+    y_norm = flat_norm(y.dot(y), ord=2)
+    if y_norm < 1.0e-14:
+        y_norm = 1.0
+
+    d = x - y
+    return actx.to_numpy(flat_norm(d.dot(d), ord=2) / y_norm)
+
+
 def run_exterior_stokes(actx_factory, *,
         ambient_dim, target_order, qbx_order, resolution,
         fmm_order=False,    # FIXME: FMM is slower than direct evaluation
@@ -198,20 +210,12 @@ def run_exterior_stokes(actx_factory, *,
 
     # {{{ check velocity at "point_target"
 
-    def rnorm2(x, y):
-        y_norm = actx.np.linalg.norm(y.dot(y), ord=2)
-        if y_norm < 1.0e-14:
-            y_norm = 1.0
-
-        d = x - y
-        return actx.to_numpy(actx.np.linalg.norm(d.dot(d), ord=2) / y_norm)
-
     ps_velocity = bind(places, sym_velocity,
             auto_where=("source", "point_target"))(actx, sigma=sigma, **op_context)
     ex_velocity = bind(places, sym_source_pot,
             auto_where=("point_source", "point_target"))(actx, sigma=charges, mu=mu)
 
-    v_error = rnorm2(ps_velocity, ex_velocity)
+    v_error = rnorm2(actx, ps_velocity, ex_velocity)
     h_max = actx.to_numpy(
             bind(places, sym.h_max(ambient_dim))(actx)
             )
@@ -257,7 +261,6 @@ def test_exterior_stokes(actx_factory, ambient_dim, visualize=False):
     target_order = 3
     qbx_order = 3
 
-    print(ambient_dim)
     if ambient_dim == 2:
         resolutions = [20, 35, 50]
     elif ambient_dim == 3:
@@ -285,6 +288,168 @@ def test_exterior_stokes(actx_factory, ambient_dim, visualize=False):
 # }}}
 
 
+# {{{ test Stokeslet identity
+
+def run_stokes_identity(actx_factory, case, identity, resolution, visualize=False):
+    actx = actx_factory()
+
+    qbx = case.get_layer_potential(actx, resolution, case.target_order)
+    places = GeometryCollection(qbx, auto_where=case.name)
+
+    density_discr = places.get_discretization(case.name)
+    logger.info("ndofs:     %d", density_discr.ndofs)
+    logger.info("nelements: %d", density_discr.mesh.nelements)
+
+    # {{{ evaluate
+
+    result = bind(places, identity.apply_operator())(actx)
+    ref_result = bind(places, identity.ref_result())(actx)
+
+    h_max = actx.to_numpy(
+            bind(places, sym.h_max(places.ambient_dim))(actx)
+            )
+    error = rnorm2(actx, result, ref_result)
+    logger.info("resolution %4d h_max %.5e error %.5e",
+            resolution, h_max, error)
+
+    # }}}
+
+    if visualize:
+        filename = "stokes_{}_{}d_resolution_{}".format(
+                type(identity).__name__.lower(), places.ambient_dim, resolution)
+
+        if places.ambient_dim == 2:
+            from meshmode.dof_array import flatten_to_numpy
+            result = flatten_to_numpy(actx, result)
+            ref_result = flatten_to_numpy(actx, ref_result)
+
+            import matplotlib.pyplot as plt
+            fig = plt.figure()
+            ax = fig.gca()
+
+            ax.plot(result[0], "o-")
+            ax.plot(ref_result[0], "k--")
+            ax.grid()
+            fig.savefig(f"{filename}_x")
+            fig.clf()
+
+            ax = fig.gca()
+            ax.plot(result[1], "o-")
+            ax.plot(ref_result[1], "k--")
+            ax.grid()
+            fig.savefig(f"{filename}_y")
+            plt.close(fig)
+        else:
+            from meshmode.discretization.visualization import make_visualizer
+            vis = make_visualizer(actx, density_discr)
+
+            vis.write_vtk_file(f"{filename}.vtu", [
+                ("result", result), ("ref", ref_result),
+                ], overwrite=True)
+
+    return h_max, error
+
+
+class StokesletIdentity:
+    """[Pozrikidis1992] Problem 3.1.1"""
+
+    def __init__(self, ambient_dim):
+        from pytential.symbolic.stokes import StokesletWrapper
+        self.ambient_dim = ambient_dim
+        self.stokeslet = StokesletWrapper(self.ambient_dim)
+
+    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)
+
+    def ref_result(self):
+        return make_obj_array([0] * self.ambient_dim)
+
+
+@pytest.mark.parametrize("cls", [
+    partial(eid.StarfishTestCase, resolutions=[16, 32, 64, 96, 128]),
+    partial(eid.SphereTestCase, resolutions=[0, 1, 2, 3]),
+    ])
+def test_stokeslet_identity(actx_factory, cls, visualize=False):
+    if visualize:
+        logging.basicConfig(level=logging.INFO)
+
+    from pytools.convergence import EOCRecorder
+    eoc = EOCRecorder()
+
+    source_ovsmp = 4 if cls.ambient_dim == 2 else 8
+    case = cls(fmm_backend=None,
+            target_order=5, qbx_order=5, source_ovsmp=source_ovsmp)
+    identity = StokesletIdentity(case.ambient_dim)
+    logger.info("\n%s", str(case))
+
+    for resolution in case.resolutions:
+        h_max, err = run_stokes_identity(
+                actx_factory, case, identity,
+                resolution=resolution,
+                visualize=visualize)
+
+        eoc.add_data_point(h_max, err)
+
+    print(eoc)
+    assert eoc.order_estimate() > case.target_order - 0.5
+
+# }}}
+
+
+# {{{ test Stresslet identity
+
+class StressletIdentity:
+    """[Pozrikidis1992] Equation 3.2.7"""
+
+    def __init__(self, ambient_dim):
+        from pytential.symbolic.stokes import StokesletWrapper
+        self.ambient_dim = ambient_dim
+        self.stokeslet = StokesletWrapper(self.ambient_dim)
+
+    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")
+
+    def ref_result(self):
+        return -0.5 * sym.normal(self.ambient_dim).as_vector()
+
+
+@pytest.mark.parametrize("cls", [
+    partial(eid.StarfishTestCase, resolutions=[16, 32, 64, 96, 128]),
+    partial(eid.SphereTestCase, resolutions=[0, 1, 2, 3]),
+    ])
+def test_stresslet_identity(actx_factory, cls, visualize=False):
+    if visualize:
+        logging.basicConfig(level=logging.INFO)
+
+    from pytools.convergence import EOCRecorder
+    eoc = EOCRecorder()
+
+    source_ovsmp = 4 if cls.ambient_dim == 2 else 8
+    case = cls(fmm_backend=None,
+            target_order=5, qbx_order=5, source_ovsmp=source_ovsmp)
+    identity = StressletIdentity(case.ambient_dim)
+    logger.info("\n%s", str(case))
+
+    for resolution in case.resolutions:
+        h_max, err = run_stokes_identity(
+                actx_factory, case, identity,
+                resolution=resolution,
+                visualize=visualize)
+
+        eoc.add_data_point(h_max, err)
+
+    print(eoc)
+    assert eoc.order_estimate() > case.target_order - 0.5
+
+# }}}
+
+
 # You can test individual routines by typing
 # $ python test_stokes.py 'test_routine()'