add contextmanager to play with fudge factor

pytential/symbolic/primitives.py

@@ -23,6 +23,7 @@
 from sys import intern
 from warnings import warn
 from functools import wraps, partial
+from contextlib import contextmanager
 
 import numpy as np
 from pymbolic.primitives import (  # noqa: F401,N813
@@ -261,36 +262,38 @@ class QBX_SOURCE_STAGE1:   # noqa: N801
     """Symbolic identifier for the Stage 1 discretization of a
     :class:`pytential.qbx.QBXLayerPotentialSource`.
     """
-    pass
 
 
 class QBX_SOURCE_STAGE2:   # noqa: N801
     """Symbolic identifier for the Stage 2 discretization of a
     :class:`pytential.qbx.QBXLayerPotentialSource`.
     """
-    pass
 
 
-class QBX_SOURCE_QUAD_STAGE2:   # noqa: N801
+class QBX_SOURCE_GEOMETRY_STAGE2:       # noqa: N801
     """Symbolic identifier for the upsampled Stage 2 discretization of a
-    :class:`pytential.qbx.QBXLayerPotentialSource`.
+    :class:`pytential.qbx.QBXLayerPotentialSource` that can
+    be used for taking geometry derivatives.
+    """
+
+
+class QBX_SOURCE_QUAD_STAGE2:           # noqa: N801
+    """Symbolic identifier for the upsampled Stage 2 discretization of a
+    :class:`pytential.qbx.QBXLayerPotentialSource` that can be used for
+    quadrature.
     """
-    pass
 
 
 class GRANULARITY_NODE:     # noqa: N801
     """DOFs are per node."""
-    pass
 
 
 class GRANULARITY_CENTER:   # noqa: N801
     """DOFs interleaved per expansion center (two per node, one on each side)."""
-    pass
 
 
 class GRANULARITY_ELEMENT:  # noqa: N801
     """DOFs per discretization element."""
-    pass
 
 
 class DOFDescriptor:
@@ -734,6 +737,7 @@ def sqrt_jac_q_weight(ambient_dim, dim=None, dofdesc=None):
 @_deprecate_kwargs("where", "dofdesc")
 def normal(ambient_dim, dim=None, dofdesc=None):
     """Exterior unit normals."""
+    # return nodes(ambient_dim, dofdesc=dofdesc)
 
     # Don't be tempted to add a sign here. As it is, it produces
     # exterior normals for positively oriented curves and surfaces.
@@ -1008,11 +1012,29 @@ def _scaled_max_curvature(ambient_dim, dim=None, dofdesc=None):
 
 # {{{ qbx-specific geometry
 
+_EXPANSION_RADII_FUDGE_FACTOR = None
+
+
+@contextmanager
+def _expansion_radii_fudge_factor(factor: float):
+    global _EXPANSION_RADII_FUDGE_FACTOR
+    _EXPANSION_RADII_FUDGE_FACTOR = factor
+
+    try:
+        yield None
+    finally:
+        _EXPANSION_RADII_FUDGE_FACTOR = None
+
+
 def _expansion_radii_factor(ambient_dim, dim):
     if dim is None:
         dim = ambient_dim - 1
 
-    dim_fudge_factor = 0.5 if dim == 2 else 1.0
+    if _EXPANSION_RADII_FUDGE_FACTOR is None:
+        dim_fudge_factor = 0.5 if dim == 2 else 1.0
+    else:
+        dim_fudge_factor = _EXPANSION_RADII_FUDGE_FACTOR
+
     return 0.5 * dim_fudge_factor
 
 
@@ -1064,9 +1086,13 @@ def _close_target_tunnel_radii(ambient_dim, dim=None,
 
 @_deprecate_kwargs("where", "dofdesc")
 def expansion_radii(ambient_dim, dim=None, granularity=None, dofdesc=None):
+    dofdesc = as_dofdesc(dofdesc)
+
     factor = _expansion_radii_factor(ambient_dim, dim)
-    return cse(factor * _quad_resolution(ambient_dim, dim=dim,
-        granularity=granularity, dofdesc=dofdesc),
+    element_radius = _quad_resolution(ambient_dim, dim=dim,
+            granularity=granularity, dofdesc=dofdesc)
+
+    return cse(factor * element_radius,
         "expansion_radii",
         cse_scope.DISCRETIZATION)
 
@@ -1108,6 +1134,17 @@ def h_max(ambient_dim, dim=None, dofdesc=None):
             cse_scope.DISCRETIZATION)
 
 
+def h_min(ambient_dim, dim=None, dofdesc=None):
+    """Defines a maximum element size in the discretization."""
+
+    dofdesc = as_dofdesc(dofdesc).copy(granularity=GRANULARITY_ELEMENT)
+    r = _quad_resolution(ambient_dim, dim=dim, dofdesc=dofdesc)
+
+    return cse(NodeMin(r),
+            "h_min",
+            cse_scope.DISCRETIZATION)
+
+
 @_deprecate_kwargs("where", "dofdesc")
 def weights_and_area_elements(ambient_dim, dim=None, dofdesc=None):
     """Combines :func:`area_element` and :class:`QWeight`."""

test/test_stokes.py

@@ -304,12 +304,15 @@ def run_stokes_identity(actx_factory, case, identity, resolution, visualize=Fals
     result = bind(places, identity.apply_operator())(actx)
     ref_result = bind(places, identity.ref_result())(actx)
 
+    h_min = actx.to_numpy(
+            bind(places, sym.h_min(places.ambient_dim))(actx)
+            )
     h_max = actx.to_numpy(
             bind(places, sym.h_max(places.ambient_dim))(actx)
             )
     error = rnorm2(actx, density_discr, result, ref_result)
-    logger.info("resolution %4d h_max %.5e error %.5e",
-            resolution, h_max, error)
+    logger.info("resolution %4d h_min %.5e h_max %.5e error %.5e",
+            resolution, h_min, h_max, error)
 
     # }}}
 
@@ -345,11 +348,13 @@ def run_stokes_identity(actx_factory, case, identity, resolution, visualize=Fals
             plt.close(fig)
         else:
             from meshmode.discretization.visualization import make_visualizer
-            vis = make_visualizer(actx, density_discr, vis_order=case.target_order)
+            vis = make_visualizer(actx, density_discr,
+                    vis_order=case.target_order,
+                    force_equidistant=True)
 
             vis.write_vtk_file(f"{filename}.vtu", [
                 ("result", result), ("ref", ref_result),
-                ], overwrite=True)
+                ], use_high_order=True, overwrite=True)
 
     return h_max, error
 
@@ -383,21 +388,28 @@ def test_stokeslet_identity(actx_factory, cls, visualize=False):
     from pytools.convergence import EOCRecorder
     eoc = EOCRecorder()
 
+    from pytential.symbolic.primitives import _expansion_radii_fudge_factor
+
     source_ovsmp = 4 if cls.func.ambient_dim == 2 else 4
-    case = cls(fmm_backend=None,
+    case = cls(fmm_backend=None, fmm_tol=None,
             target_order=3, qbx_order=3, 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)
+    with _expansion_radii_fudge_factor(1.0):
+        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)
+            eoc.add_data_point(h_max, err)
+
+    print(eoc.pretty_print(
+        abscissa_format="%.8e",
+        error_format="%.8e",
+        eoc_format="%.2f"))
 
-    print(eoc)
     # assert eoc.order_estimate() > case.target_order - 0.5
 
 # }}}