Grudge Refactoring

doc/operators.rst

@@ -1,5 +1,16 @@
-Discontinuous Galerkin Operators
+Discontinuous Galerkin operators
 ================================
 
 .. automodule:: grudge.op
 .. automodule:: grudge.trace_pair
+
+
+Transfering data between discretizations
+========================================
+
+.. automodule:: grudge.projection
+
+Reductions
+==========
+
+.. automodule:: grudge.reductions

examples/advection/surface.py

@@ -172,7 +172,7 @@ def main(ctx_factory, dim=2, order=4, use_quad=False, visualize=False):
     # {{{ Surface advection operator
 
     # velocity field
-    x = thaw(op.nodes(dcoll), actx)
+    x = thaw(dcoll.nodes(), actx)
     c = make_obj_array([-x[1], x[0], 0.0])[:dim]
 
     def f_initial_condition(x):
@@ -234,7 +234,7 @@ def rhs(t, u):
 
         df = dof_desc.DOFDesc(FACE_RESTR_INTERIOR)
         face_discr = dcoll.discr_from_dd(df)
-        face_normal = thaw(op.normal(dcoll, dd=df), actx)
+        face_normal = thaw(dcoll.normal(dd=df), actx)
 
         from meshmode.discretization.visualization import make_visualizer
         vis = make_visualizer(actx, face_discr)

examples/advection/var-velocity.py

@@ -169,7 +169,7 @@ def zero_inflow_bc(dtag, t=0):
 
     from grudge.models.advection import VariableCoefficientAdvectionOperator
 
-    x = thaw(op.nodes(dcoll), actx)
+    x = thaw(dcoll.nodes(), actx)
 
     # velocity field
     if dim == 1:

examples/advection/weak.py

@@ -152,13 +152,13 @@ def u_analytic(x, t=0):
         dcoll,
         c,
         inflow_u=lambda t: u_analytic(
-            thaw(op.nodes(dcoll, dd=BTAG_ALL), actx),
+            thaw(dcoll.nodes(dd=BTAG_ALL), actx),
             t=t
         ),
         flux_type=flux_type
     )
 
-    nodes = thaw(op.nodes(dcoll), actx)
+    nodes = thaw(dcoll.nodes(), actx)
     u = u_analytic(nodes, t=0)
 
     def rhs(t, u):

examples/geometry.py

@@ -32,8 +32,6 @@
 
 from arraycontext import PyOpenCLArrayContext, thaw
 
-import grudge.op as op
-
 from grudge import DiscretizationCollection, shortcuts
 
 
@@ -51,9 +49,9 @@ def main(write_output=True):
 
     dcoll = DiscretizationCollection(actx, mesh, order=4)
 
-    nodes = thaw(op.nodes(dcoll), actx)
-    bdry_nodes = thaw(op.nodes(dcoll, dd=BTAG_ALL), actx)
-    bdry_normals = thaw(op.normal(dcoll, dd=BTAG_ALL), actx)
+    nodes = thaw(dcoll.nodes(), actx)
+    bdry_nodes = thaw(dcoll.nodes(dd=BTAG_ALL), actx)
+    bdry_normals = thaw(dcoll.normal(dd=BTAG_ALL), actx)
 
     if write_output:
         vis = shortcuts.make_visualizer(dcoll)

examples/maxwell/cavities.py

@@ -84,7 +84,7 @@ def cavity_mode(x, t=0):
         else:
             return get_rectangular_cavity_mode(actx, x, t, 1, (2, 3))
 
-    fields = cavity_mode(thaw(op.nodes(dcoll), actx), t=0)
+    fields = cavity_mode(thaw(dcoll.nodes(), actx), t=0)
 
     maxwell_operator.check_bc_coverage(mesh)
 

examples/wave/var-propagation-speed.py

@@ -63,7 +63,7 @@ def source_f(actx, dcoll, t=0):
         source_center = np.array([0.1, 0.22, 0.33])[:dcoll.dim]
         source_width = 0.05
         source_omega = 3
-        nodes = thaw(op.nodes(dcoll), actx)
+        nodes = thaw(dcoll.nodes(), actx)
         source_center_dist = flat_obj_array(
             [nodes[i] - source_center[i] for i in range(dcoll.dim)]
         )
@@ -75,7 +75,7 @@ def source_f(actx, dcoll, t=0):
             )
         )
 
-    x = thaw(op.nodes(dcoll), actx)
+    x = thaw(dcoll.nodes(), actx)
     ones = dcoll.zeros(actx) + 1
     c = actx.np.where(np.dot(x, x) < 0.15, 0.1 * ones, 0.2 * ones)
 

examples/wave/wave-min-mpi.py

@@ -84,7 +84,7 @@ def source_f(actx, dcoll, t=0):
         source_center = np.array([0.1, 0.22, 0.33])[:dcoll.dim]
         source_width = 0.05
         source_omega = 3
-        nodes = thaw(op.nodes(dcoll), actx)
+        nodes = thaw(dcoll.nodes(), actx)
         source_center_dist = flat_obj_array(
             [nodes[i] - source_center[i] for i in range(dcoll.dim)]
         )

examples/wave/wave-min.py

@@ -65,7 +65,7 @@ def source_f(actx, dcoll, t=0):
         source_center = np.array([0.1, 0.22, 0.33])[:dcoll.dim]
         source_width = 0.05
         source_omega = 3
-        nodes = thaw(op.nodes(dcoll), actx)
+        nodes = thaw(dcoll.nodes(), actx)
         source_center_dist = flat_obj_array(
             [nodes[i] - source_center[i] for i in range(dcoll.dim)]
         )

examples/wave/wave-op-mpi.py

@@ -54,7 +54,7 @@ def wave_flux(dcoll, c, w_tpair):
     u = w_tpair[0]
     v = w_tpair[1:]
 
-    normal = thaw(op.normal(dcoll, w_tpair.dd), u.int.array_context)
+    normal = thaw(dcoll.normal(w_tpair.dd), u.int.array_context)
 
     flux_weak = flat_obj_array(
             np.dot(v.avg, normal),
@@ -136,7 +136,7 @@ def bump(actx, dcoll, t=0):
     source_width = 0.05
     source_omega = 3
 
-    nodes = thaw(op.nodes(dcoll), actx)
+    nodes = thaw(dcoll.nodes(), actx)
     center_dist = flat_obj_array([
         nodes[i] - source_center[i]
         for i in range(dcoll.dim)

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

@@ -56,7 +56,7 @@ def wave_flux(dcoll, c, w_tpair):
     u = w_tpair[0]
     v = w_tpair[1:]
 
-    normal = thaw(op.normal(dcoll, dd), u.int.array_context)
+    normal = thaw(dcoll.normal(dd), u.int.array_context)
 
     flux_weak = flat_obj_array(
             np.dot(v.avg, normal),
@@ -146,7 +146,7 @@ def bump(actx, dcoll, t=0, width=0.05, center=None):
     center = center[:dcoll.dim]
     source_omega = 3
 
-    nodes = thaw(op.nodes(dcoll), actx)
+    nodes = thaw(dcoll.nodes(), actx)
     center_dist = flat_obj_array([
         nodes[i] - center[i]
         for i in range(dcoll.dim)

examples/wave/wave-op.py

@@ -52,7 +52,7 @@ def wave_flux(dcoll, c, w_tpair):
     u = w_tpair[0]
     v = w_tpair[1:]
 
-    normal = thaw(op.normal(dcoll, w_tpair.dd), u.int.array_context)
+    normal = thaw(dcoll.normal(w_tpair.dd), u.int.array_context)
 
     flux_weak = flat_obj_array(
             np.dot(v.avg, normal),
@@ -127,7 +127,7 @@ def bump(actx, dcoll, t=0):
     source_width = 0.05
     source_omega = 3
 
-    nodes = thaw(op.nodes(dcoll), actx)
+    nodes = thaw(dcoll.nodes(), actx)
     center_dist = flat_obj_array([
         nodes[i] - source_center[i]
         for i in range(dcoll.dim)

grudge/dt_utils.py

@@ -2,6 +2,8 @@
 
 .. autofunction:: dt_non_geometric_factor
 .. autofunction:: dt_geometric_factors
+.. autofunction:: h_max_from_volume
+.. autofunction:: h_min_from_volume
 """
 
 __copyright__ = """
@@ -33,7 +35,7 @@
 
 from arraycontext import FirstAxisIsElementsTag
 
-from grudge.dof_desc import DD_VOLUME, DOFDesc
+from grudge.dof_desc import DD_VOLUME, DOFDesc, as_dofdesc
 from grudge.discretization import DiscretizationCollection
 
 import grudge.op as op
@@ -90,6 +92,70 @@ def dt_non_geometric_factor(dcoll: DiscretizationCollection, dd=None) -> float:
     return min(min_delta_rs)
 
 
+# {{{ Mesh size functions
+
+@memoize_on_first_arg
+def h_max_from_volume(
+        dcoll: DiscretizationCollection, dim=None, dd=None) -> float:
+    """Returns a (maximum) characteristic length based on the volume of the
+    elements. This length may not be representative if the elements have very
+    high aspect ratios.
+
+    :arg dim: an integer denoting topological dimension. If *None*, the
+        spatial dimension specified by
+        :attr:`grudge.DiscretizationCollection.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: a scalar denoting the maximum characteristic length.
+    """
+    from grudge.reductions import nodal_max, elementwise_sum
+
+    if dd is None:
+        dd = DD_VOLUME
+    dd = as_dofdesc(dd)
+
+    if dim is None:
+        dim = dcoll.dim
+
+    ones = dcoll.discr_from_dd(dd).zeros(dcoll._setup_actx) + 1.0
+    return nodal_max(
+        dcoll,
+        dd,
+        elementwise_sum(dcoll, op.mass(dcoll, dd, ones))
+    ) ** (1.0 / dim)
+
+
+@memoize_on_first_arg
+def h_min_from_volume(
+        dcoll: DiscretizationCollection, dim=None, dd=None) -> float:
+    """Returns a (minimum) characteristic length based on the volume of the
+    elements. This length may not be representative if the elements have very
+    high aspect ratios.
+
+    :arg dim: an integer denoting topological dimension. If *None*, the
+        spatial dimension specified by
+        :attr:`grudge.DiscretizationCollection.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: a scalar denoting the minimum characteristic length.
+    """
+    from grudge.reductions import nodal_min, elementwise_sum
+
+    if dd is None:
+        dd = DD_VOLUME
+    dd = as_dofdesc(dd)
+
+    if dim is None:
+        dim = dcoll.dim
+
+    ones = dcoll.discr_from_dd(dd).zeros(dcoll._setup_actx) + 1.0
+    return nodal_min(
+        dcoll,
+        dd,
+        elementwise_sum(dcoll, op.mass(dcoll, dd, ones))
+    ) ** (1.0 / dim)
+
+
 @memoize_on_first_arg
 def dt_geometric_factors(
         dcoll: DiscretizationCollection, dd=None) -> DOFArray:
@@ -175,3 +241,8 @@ def dt_geometric_factors(
             for cv_i, sae_i in zip(cell_vols, surface_areas)
         )
     )
+
+# }}}
+
+
+# vim: foldmethod=marker

grudge/interpolation.py

@@ -0,0 +1,48 @@
+"""
+.. currentmodule:: grudge.op
+
+Interpolation
+-------------
+
+.. autofunction:: interp
+"""
+
+__copyright__ = """
+Copyright (C) 2021 University of Illinois Board of Trustees
+"""
+
+__license__ = """
+Permission is hereby granted, free of charge, to any person obtaining a copy
+of this software and associated documentation files (the "Software"), to deal
+in the Software without restriction, including without limitation the rights
+to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
+copies of the Software, and to permit persons to whom the Software is
+furnished to do so, subject to the following conditions:
+
+The above copyright notice and this permission notice shall be included in
+all copies or substantial portions of the Software.
+
+THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
+OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
+THE SOFTWARE.
+"""
+
+
+from grudge.discretization import DiscretizationCollection
+
+
+# FIXME: Should revamp interp and make clear distinctions
+# between projection and interpolations.
+# Related issue: https://github.com/inducer/grudge/issues/38
+def interp(dcoll: DiscretizationCollection, src, tgt, vec):
+    from warnings import warn
+    warn("'interp' currently calls to 'project'",
+         UserWarning, stacklevel=2)
+
+    from grudge.projection import project
+
+    return project(dcoll, src, tgt, vec)

grudge/models/advection.py

@@ -43,7 +43,7 @@ def advection_weak_flux(dcoll, flux_type, u_tpair, velocity):
     """
     actx = u_tpair.int.array_context
     dd = u_tpair.dd
-    normal = thaw(op.normal(dcoll, dd), actx)
+    normal = thaw(dcoll.normal(dd), actx)
     v_dot_n = np.dot(velocity, normal)
 
     flux_type = flux_type.lower()
@@ -92,7 +92,7 @@ class StrongAdvectionOperator(AdvectionOperatorBase):
     def flux(self, u_tpair):
         actx = u_tpair.int.array_context
         dd = u_tpair.dd
-        normal = thaw(op.normal(self.dcoll, dd), actx)
+        normal = thaw(self.dcoll.normal(dd), actx)
         v_dot_normal = np.dot(self.v, normal)
 
         return u_tpair.int * v_dot_normal - self.weak_flux(u_tpair)
@@ -285,7 +285,7 @@ def v_dot_n_tpair(actx, dcoll, velocity, trace_dd):
     from grudge.trace_pair import TracePair
     from meshmode.discretization.connection import FACE_RESTR_INTERIOR
 
-    normal = thaw(op.normal(dcoll, trace_dd.with_discr_tag(None)), actx)
+    normal = thaw(dcoll.normal(trace_dd.with_discr_tag(None)), actx)
     v_dot_n = velocity.dot(normal)
     i = op.project(dcoll, trace_dd.with_discr_tag(None), trace_dd, v_dot_n)
 

grudge/models/em.py

@@ -120,7 +120,7 @@ def flux(self, wtpair):
         As per Hesthaven and Warburton page 433.
         """
 
-        normal = thaw(op.normal(self.dcoll, wtpair.dd), self.dcoll._setup_actx)
+        normal = thaw(self.dcoll.normal(wtpair.dd), self.dcoll._setup_actx)
 
         if self.fixed_material:
             e, h = self.split_eh(wtpair)
@@ -220,7 +220,7 @@ def absorbing_bc(self, w):
         absorbing boundary conditions.
         """
 
-        absorb_normal = thaw(op.normal(self.dcoll, dd=self.absorb_tag),
+        absorb_normal = thaw(self.dcoll.normal(dd=self.absorb_tag),
                              self.dcoll._setup_actx)
 
         e, h = self.split_eh(w)