Merge branch 'main' into dcoll-refactor

Author: thomasgibson

examples/advection/surface.py

@@ -146,11 +146,11 @@ def main(ctx_factory, dim=2, order=4, use_quad=False, visualize=False):
         qtag = None
 
     from meshmode.discretization.poly_element import \
-            PolynomialWarpAndBlendGroupFactory, \
+            default_simplex_group_factory, \
             QuadratureSimplexGroupFactory
 
     discr_tag_to_group_factory[dof_desc.DISCR_TAG_BASE] = \
-        PolynomialWarpAndBlendGroupFactory(order)
+        default_simplex_group_factory(base_dim=dim-1, order=order)
 
     if use_quad:
         discr_tag_to_group_factory[qtag] = \
@@ -203,7 +203,8 @@ def rhs(t, u):
 
     # {{{ time stepping
 
-    dt = adv_operator.estimate_rk4_timestep(dcoll, fields=u0)
+    # FIXME: dt estimate is not necessarily valid for surfaces
+    dt = 0.45 * adv_operator.estimate_rk4_timestep(actx, dcoll, fields=u0)
     nsteps = int(final_time // dt) + 1
 
     logger.info("dt:        %.5e", dt)

examples/advection/var-velocity.py

@@ -195,7 +195,7 @@ def zero_inflow_bc(dtag, t=0):
     def rhs(t, u):
         return adv_operator.operator(t, u)
 
-    dt = adv_operator.estimate_rk4_timestep(dcoll, fields=u)
+    dt = adv_operator.estimate_rk4_timestep(actx, dcoll, fields=u)
 
     logger.info("Timestep size: %g", dt)
 

examples/advection/weak.py

@@ -164,7 +164,7 @@ def u_analytic(x, t=0):
     def rhs(t, u):
         return adv_operator.operator(t, u)
 
-    dt = adv_operator.estimate_rk4_timestep(dcoll, fields=u)
+    dt = adv_operator.estimate_rk4_timestep(actx, dcoll, fields=u)
 
     logger.info("Timestep size: %g", dt)
 

examples/maxwell/cavities.py

@@ -91,7 +91,7 @@ def cavity_mode(x, t=0):
     def rhs(t, w):
         return maxwell_operator.operator(t, w)
 
-    dt = maxwell_operator.estimate_rk4_timestep(dcoll, fields=fields)
+    dt = maxwell_operator.estimate_rk4_timestep(actx, dcoll, fields=fields)
 
     dt_stepper = set_up_rk4("w", dt, fields, rhs)
 

examples/wave/var-propagation-speed.py

@@ -102,9 +102,7 @@ def source_f(actx, dcoll, t=0):
     def rhs(t, w):
         return wave_op.operator(t, w)
 
-    dt_scaling_const = 2/3
-    dt = dt_scaling_const * wave_op.estimate_rk4_timestep(dcoll, fields=fields)
-
+    dt = 2/3 * wave_op.estimate_rk4_timestep(actx, dcoll, fields=fields)
     dt_stepper = set_up_rk4("w", dt, fields, rhs)
 
     final_t = 1

examples/wave/wave-min-mpi.py

@@ -116,8 +116,7 @@ def source_f(actx, dcoll, t=0):
         [dcoll.zeros(actx) for i in range(dcoll.dim)]
     )
 
-    dt_scaling_const = 2/3
-    dt = dt_scaling_const * wave_op.estimate_rk4_timestep(dcoll, fields=fields)
+    dt = 2/3 * wave_op.estimate_rk4_timestep(actx, dcoll, fields=fields)
 
     wave_op.check_bc_coverage(local_mesh)
 

examples/wave/wave-op-mpi.py

@@ -114,14 +114,12 @@ def rk4_step(y, t, h, f):
     return y + h/6*(k1 + 2*k2 + 2*k3 + k4)
 
 
-def estimate_rk4_timestep(dcoll, c):
-    from grudge.dt_utils import (dt_non_geometric_factor,
-                                 dt_geometric_factors)
+def estimate_rk4_timestep(actx, dcoll, c):
+    from grudge.dt_utils import characteristic_lengthscales
 
-    dt_factor = (dt_non_geometric_factor(dcoll)
-                 * op.nodal_min(dcoll, "vol", dt_geometric_factors(dcoll)))
+    local_dts = characteristic_lengthscales(actx, dcoll) / c
 
-    return dt_factor * (1 / c)
+    return op.nodal_min(dcoll, "vol", local_dts)
 
 
 def bump(actx, dcoll, t=0):
@@ -187,8 +185,7 @@ def main(ctx_factory, dim=2, order=3, visualize=False):
             )
 
     c = 1
-    dt_scaling_const = 0.45
-    dt = dt_scaling_const * estimate_rk4_timestep(dcoll, c)
+    dt = 0.45 * estimate_rk4_timestep(actx, dcoll, c)
 
     vis = make_visualizer(dcoll)
 

examples/wave/wave-op-var-velocity.py

@@ -129,14 +129,12 @@ def rk4_step(y, t, h, f):
     return y + h/6*(k1 + 2*k2 + 2*k3 + k4)
 
 
-def estimate_rk4_timestep(dcoll, c):
-    from grudge.dt_utils import (dt_non_geometric_factor,
-                                 dt_geometric_factors)
+def estimate_rk4_timestep(actx, dcoll, c):
+    from grudge.dt_utils import characteristic_lengthscales
 
-    dt_factor = (dt_non_geometric_factor(dcoll)
-                 * op.nodal_min(dcoll, "vol", dt_geometric_factors(dcoll)))
+    local_dts = characteristic_lengthscales(actx, dcoll) / c
 
-    return dt_factor * (1 / c)
+    return op.nodal_min(dcoll, "vol", local_dts)
 
 
 def bump(actx, dcoll, t=0, width=0.05, center=None):
@@ -182,15 +180,14 @@ def main(ctx_factory, dim=2, order=3, visualize=False):
     dcoll = make_discretization_collection(
         actx, mesh,
         discr_tag_to_group_factory={
-            DISCR_TAG_BASE: PolynomialWarpAndBlendGroupFactory(order),
+            DISCR_TAG_BASE: default_simplex_group_factory(base_dim=dim, order=order),
             DISCR_TAG_QUAD: QuadratureSimplexGroupFactory(3*order),
         }
     )
 
     # bounded above by 1
     c = 0.2 + 0.8*bump(actx, dcoll, center=np.zeros(3), width=0.5)
-    dt_scaling_const = 0.5
-    dt = dt_scaling_const * estimate_rk4_timestep(dcoll, c=1)
+    dt = 0.5 * estimate_rk4_timestep(actx, dcoll, c=1)
 
     fields = flat_obj_array(
             bump(actx, dcoll, ),

grudge/discretization.py

@@ -581,7 +581,7 @@ def make_discretization_collection(
     :returns: A :class:`DiscretizationCollection`.
     """
     from meshmode.discretization.poly_element import \
-        PolynomialWarpAndBlendGroupFactory
+        default_simplex_group_factory
 
     if discr_tag_to_group_factory is None:
         if order is None:
@@ -591,7 +591,8 @@ def make_discretization_collection(
 
         # Default choice: warp and blend simplex element group
         discr_tag_to_group_factory = {
-            DISCR_TAG_BASE: PolynomialWarpAndBlendGroupFactory(order=order)
+            DISCR_TAG_BASE: default_simplex_group_factory(base_dim=mesh.dim,
+                                                          order=order)
         }
     else:
         if order is not None:
@@ -603,7 +604,7 @@ def make_discretization_collection(
                 )
 
             discr_tag_to_group_factory[DISCR_TAG_BASE] = \
-                PolynomialWarpAndBlendGroupFactory(order=order)
+                default_simplex_group_factory(base_dim=mesh.dim, order=order)
 
     # Modal discr should always comes from the base discretization
     discr_tag_to_group_factory[DISCR_TAG_MODAL] = \

grudge/dt_utils.py

@@ -1,6 +1,18 @@
 """Helper functions for estimating stable time steps for RKDG methods.
 
-.. autofunction:: dt_non_geometric_factor
+Characteristic lengthscales
+---------------------------
+
+.. autofunction:: characteristic_lengthscales
+
+Non-geometric quantities
+------------------------
+
+.. autofunction:: dt_non_geometric_factors
+
+Mesh size utilities
+-------------------
+
 .. autofunction:: dt_geometric_factors
 .. autofunction:: h_max_from_volume
 .. autofunction:: h_min_from_volume
@@ -33,7 +45,7 @@
 
 import numpy as np
 
-from arraycontext import FirstAxisIsElementsTag
+from arraycontext import ArrayContext, FirstAxisIsElementsTag, thaw, freeze
 
 from grudge.dof_desc import DD_VOLUME, DOFDesc, as_dofdesc
 from grudge.discretization import DiscretizationCollection
@@ -42,25 +54,75 @@
 
 from meshmode.dof_array import DOFArray
 
-from pytools import memoize_on_first_arg
+from pytools import memoize_on_first_arg, memoize_in
+
+
+def characteristic_lengthscales(
+        actx: ArrayContext, dcoll: DiscretizationCollection) -> DOFArray:
+    r"""Computes the characteristic length scale :math:`h_{\text{loc}}` at
+    each node. The characteristic length scale is mainly useful for estimating
+    the stable time step size. E.g. for a hyperbolic system, an estimate of the
+    stable time step can be estimated as :math:`h_{\text{loc}} / c`, where
+    :math:`c` is the characteristic wave speed. The estimate is obtained using
+    the following formula:
+
+    .. math::
+
+        h_{\text{loc}} = \operatorname{min}\left(\Delta r_i\right) r_D
+
+    where :math:`\operatorname{min}\left(\Delta r_i\right)` is the minimum
+    node distance on the reference cell (see :func:`dt_non_geometric_factors`),
+    and :math:`r_D` is the inradius of the cell (see :func:`dt_geometric_factors`).
+
+    :returns: a frozen :class:`~meshmode.dof_array.DOFArray` containing a
+        characteristic lengthscale for each element, at each nodal location.
+
+    .. note::
+
+        While a prediction of stability is only meaningful in relation to a given
+        time integrator with a known stability region, the lengthscale provided here
+        is not intended to be specific to any one time integrator, though the
+        stability region of standard four-stage, fourth-order Runge-Kutta
+        methods has been used as a guide. Any concrete time integrator will
+        likely require scaling of the values returned by this routine.
+    """
+    @memoize_in(dcoll, (characteristic_lengthscales,
+                        "compute_characteristic_lengthscales"))
+    def _compute_characteristic_lengthscales():
+        return freeze(
+            DOFArray(
+                actx,
+                data=tuple(
+                    # Scale each group array of geometric factors by the
+                    # corresponding group non-geometric factor
+                    cng * geo_facts
+                    for cng, geo_facts in zip(
+                        dt_non_geometric_factors(dcoll),
+                        thaw(dt_geometric_factors(dcoll), actx)
+                    )
+                )
+            )
+        )
+    return thaw(_compute_characteristic_lengthscales(), actx)
 
 
 @memoize_on_first_arg
-def dt_non_geometric_factor(dcoll: DiscretizationCollection, dd=None) -> float:
-    r"""Computes the non-geometric scale factor following [Hesthaven_2008]_,
-    section 6.4:
+def dt_non_geometric_factors(
+        dcoll: DiscretizationCollection, dd=None) -> list:
+    r"""Computes the non-geometric scale factors following [Hesthaven_2008]_,
+    section 6.4, for each element group in the *dd* discretization:
 
     .. math::
 
         c_{ng} = \operatorname{min}\left( \Delta r_i \right),
 
     where :math:`\Delta r_i` denotes the distance between two distinct
-    nodes on the reference element.
+    nodal points on the reference element.
 
     :arg dd: a :class:`~grudge.dof_desc.DOFDesc`, or a value convertible to one.
         Defaults to the base volume discretization if not provided.
-    :returns: a :class:`float` denoting the minimum node distance on the
-        reference element.
+    :returns: a :class:`list` of :class:`float` values denoting the minimum
+        node distance on the reference element for each group.
     """
     if dd is None:
         dd = DD_VOLUME
@@ -88,8 +150,7 @@ def dt_non_geometric_factor(dcoll: DiscretizationCollection, dd=None) -> float:
                 )
             )
 
-    # Return minimum over all element groups in the discretization
-    return min(min_delta_rs)
+    return min_delta_rs
 
 
 # {{{ Mesh size functions
@@ -156,7 +217,6 @@ def h_min_from_volume(
     ) ** (1.0 / dim)
 
 
-@memoize_on_first_arg
 def dt_geometric_factors(
         dcoll: DiscretizationCollection, dd=None) -> DOFArray:
     r"""Computes a geometric scaling factor for each cell following [Hesthaven_2008]_,
@@ -175,8 +235,8 @@ def dt_geometric_factors(
 
     :arg dd: a :class:`~grudge.dof_desc.DOFDesc`, or a value convertible to one.
         Defaults to the base volume discretization if not provided.
-    :returns: a :class:`~meshmode.dof_array.DOFArray` containing the geometric
-        factors for each cell and at each nodal location.
+    :returns: a frozen :class:`~meshmode.dof_array.DOFArray` containing the
+        geometric factors for each cell and at each nodal location.
     """
     from meshmode.discretization.poly_element import SimplexElementGroupBase
 
@@ -192,14 +252,20 @@ def dt_geometric_factors(
             "Geometric factors are only implemented for simplex element groups"
         )
 
+    if volm_discr.dim != volm_discr.ambient_dim:
+        from warnings import warn
+        warn("The geometric factor for the characteristic length scale in "
+                "time step estimation is not necessarily valid for non-volume-"
+                "filling discretizations. Continuing anyway.", stacklevel=3)
+
     cell_vols = abs(
         op.elementwise_integral(
             dcoll, dd, volm_discr.zeros(actx) + 1.0
         )
     )
 
     if dcoll.dim == 1:
-        return cell_vols
+        return freeze(cell_vols)
 
     dd_face = DOFDesc("all_faces", dd.discretization_tag)
     face_discr = dcoll.discr_from_dd(dd_face)
@@ -230,7 +296,7 @@ def dt_geometric_factors(
         )
     )
 
-    return DOFArray(
+    return freeze(DOFArray(
         actx,
         data=tuple(
             actx.einsum("e,ei->ei",
@@ -240,7 +306,7 @@ def dt_geometric_factors(
 
             for cv_i, sae_i in zip(cell_vols, surface_areas)
         )
-    )
+    ))
 
 # }}}