Add op.h_max_from_volume

grudge/op.py

@@ -3,6 +3,7 @@
 
 .. autofunction:: nodes
 .. autofunction:: normal
+.. autofunction:: h_max_from_volume
 
 .. autofunction:: local_grad
 .. autofunction:: local_d_dx
@@ -153,6 +154,32 @@ def normal(dcoll, dd):
     actx = dcoll.discr_from_dd(dd)._setup_actx
     return freeze(normal(actx, dcoll, dd))
 
+
+@memoize_on_first_arg
+def h_max_from_volume(dcoll, dim=None, dd=None):
+    """Returns a characteristic length based on the volume of the elements.
+    This length may not be representative if the elements have very high
+    aspect ratios.
+
+    :arg dcoll: a :class:`grudge.discretization.DiscretizationCollection`.
+    :arg dim: an integer denoting topological dimension. If *None*, the
+        spatial dimension specified by :attr:`dcoll.dim` is used.
+    :arg dd: a :class:`~grudge.dof_desc.DOFDesc`, or a value convertible to one.
+        Defaults to the base volume discretization if not provided.
+    :returns: an integer denoting the characteristic length.
+    """
+    if dd is None:
+        dd = dof_desc.DD_VOLUME
+    dd = dof_desc.as_dofdesc(dd)
+
+    if dim is None:
+        dim = dcoll.dim
+
+    ones_volm = dcoll._volume_discr.zeros(dcoll._setup_actx) + 1.0
+    return nodal_max(
+        dcoll, dd, elementwise_sum(dcoll, mass(dcoll, dd, ones_volm))
+    ) ** (1.0 / dim)
+
 # }}}
 
 

test/test_grudge.py

@@ -249,11 +249,8 @@ def test_mass_surface_area(actx_factory, name):
 
         # }}}
 
-        # {{{ compute h_max using mass weights
-
-        h_max = actx.np.max(flattened_mass_weights) ** (1/dcoll.dim)
-
-        # }}}
+        # compute max element size
+        h_max = op.h_max_from_volume(dcoll)
 
         eoc.add_data_point(h_max, area_error)
 
@@ -315,13 +312,8 @@ def test_surface_mass_operator_inverse(actx_factory, name):
 
         # }}}
 
-        # {{{ compute h_max from mass weights
-
-        ones_volm = volume_discr.zeros(actx) + 1
-        flattened_mass_weights = flatten(op.mass(dcoll, dd, ones_volm))
-        h_max = actx.np.max(flattened_mass_weights) ** (1/dcoll.dim)
-
-        # }}}
+        # compute max element size
+        h_max = op.h_max_from_volume(dcoll)
 
         eoc.add_data_point(h_max, inv_error)
 
@@ -630,9 +622,7 @@ def test_surface_divergence_theorem(actx_factory, mesh_name, visualize=False):
         logger.info("errors: global %.5e local %.5e", err_global, err_local)
 
         # compute max element size
-        ones_volm = volume.zeros(actx) + 1
-        sum_mass_weights = op.elementwise_sum(dcoll, op.mass(dcoll, dd, ones_volm))
-        h_max = op.nodal_max(dcoll, dd, sum_mass_weights) ** (1/dcoll.dim)
+        h_max = op.h_max_from_volume(dcoll)
 
         eoc_global.add_data_point(h_max, err_global)
         eoc_local.add_data_point(h_max, err_local)