diff --git a/doc/operators.rst b/doc/operators.rst
index a51d0f2a0e0e70c0989f5ac36c70d22323e0bc40..550a78023e955d47a2a7e13aef747c6b1a0c9fbe 100644
--- a/doc/operators.rst
+++ b/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
diff --git a/examples/advection/surface.py b/examples/advection/surface.py
index 9dcdc93ec550acdb0e50d210ce470c897bb04d24..c81e23386c21f4284451f9225a76f9dc5114cf18 100644
--- a/examples/advection/surface.py
+++ b/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 main(ctx_factory, dim=2, order=4, use_quad=False, visualize=False):
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)
diff --git a/examples/advection/var-velocity.py b/examples/advection/var-velocity.py
index b050c60b8e9648a6ae51159346ae27eee5dfa50b..fca97d27276f768b781bc844623752aca0104de0 100644
--- a/examples/advection/var-velocity.py
+++ b/examples/advection/var-velocity.py
@@ -169,7 +169,7 @@ def main(ctx_factory, dim=2, order=4, use_quad=False, visualize=False):
from grudge.models.advection import VariableCoefficientAdvectionOperator
- x = thaw(op.nodes(dcoll), actx)
+ x = thaw(dcoll.nodes(), actx)
# velocity field
if dim == 1:
diff --git a/examples/advection/weak.py b/examples/advection/weak.py
index 9300ab01f4e9a7aa7288bf847429f49c05d5c914..0d0abb6af1ef3f62d984033bd6798aa07da4fa24 100644
--- a/examples/advection/weak.py
+++ b/examples/advection/weak.py
@@ -152,13 +152,13 @@ def main(ctx_factory, dim=2, order=4, visualize=False):
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):
diff --git a/examples/geometry.py b/examples/geometry.py
index faf0fc42c1553ea764dbc1621142e21c9f3f4599..a40a000a63333be7771fc6d3fd1cd3ad76c4a957 100644
--- a/examples/geometry.py
+++ b/examples/geometry.py
@@ -32,8 +32,6 @@ import pyopencl.tools as cl_tools
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)
diff --git a/examples/maxwell/cavities.py b/examples/maxwell/cavities.py
index 3d68e513fdb51b5b66801d01e90bb2dfa39860e6..5791a3d8e98866ad0cd3bae7c6ca11f7efc6a09d 100644
--- a/examples/maxwell/cavities.py
+++ b/examples/maxwell/cavities.py
@@ -84,7 +84,7 @@ def main(ctx_factory, dim=3, order=4, visualize=False):
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)
diff --git a/examples/wave/var-propagation-speed.py b/examples/wave/var-propagation-speed.py
index 0b8d51b9bdfa0388612366c1d95cedcc931eedb2..04605b7e6d6ff92f8607be09fcf3ac58fa755acd 100644
--- a/examples/wave/var-propagation-speed.py
+++ b/examples/wave/var-propagation-speed.py
@@ -63,7 +63,7 @@ def main(ctx_factory, dim=2, order=4, visualize=False):
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 main(ctx_factory, dim=2, order=4, visualize=False):
)
)
- 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)
diff --git a/examples/wave/wave-min-mpi.py b/examples/wave/wave-min-mpi.py
index a3e75c7b525628bcaa6bd861a1009ec9cd79d985..c803a897a8eb32acc7aa40c6f7eef176d5d148d1 100644
--- a/examples/wave/wave-min-mpi.py
+++ b/examples/wave/wave-min-mpi.py
@@ -84,7 +84,7 @@ def main(ctx_factory, dim=2, order=4, visualize=False):
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)]
)
diff --git a/examples/wave/wave-min.py b/examples/wave/wave-min.py
index 7b18d3dfaec68ab38af936c45f395bd10ad29a75..ab983681f41ba9b043aa7f18e4365d70cd0a93a2 100644
--- a/examples/wave/wave-min.py
+++ b/examples/wave/wave-min.py
@@ -65,7 +65,7 @@ def main(ctx_factory, dim=2, order=4, visualize=False):
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)]
)
diff --git a/examples/wave/wave-op-mpi.py b/examples/wave/wave-op-mpi.py
index 3fb49d7d1b65eb8d9caa5ac6a5d68c79625c5a94..cc8578c9983ccd24c57278f9e40dd658b08ef307 100644
--- a/examples/wave/wave-op-mpi.py
+++ b/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)
diff --git a/examples/wave/wave-op-var-velocity.py b/examples/wave/wave-op-var-velocity.py
index e48c8d7141285b40d0a0935c032449e7bbe5a2a4..07b40b8d6a5422a8d4c94962e0441308ee456953 100644
--- a/examples/wave/wave-op-var-velocity.py
+++ b/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)
diff --git a/examples/wave/wave-op.py b/examples/wave/wave-op.py
index c0494d9f0196da5cb7d951a6950f82250a2e9aa6..56a901716efd60f30f813aacc29539f2b40bf2d2 100644
--- a/examples/wave/wave-op.py
+++ b/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)
diff --git a/grudge/dt_utils.py b/grudge/dt_utils.py
index 14c955595fc5ca41ddcd707066e0238865b273e5..bed073cf100a917fd130950d36641544c1d01024 100644
--- a/grudge/dt_utils.py
+++ b/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 @@ import numpy as np
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
diff --git a/grudge/interpolation.py b/grudge/interpolation.py
new file mode 100644
index 0000000000000000000000000000000000000000..61bdf1a133b984a6c26ebb26992e15368099b611
--- /dev/null
+++ b/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)
diff --git a/grudge/models/advection.py b/grudge/models/advection.py
index 2469498aef61fff823d794123992f586cfa1cb39..0ba4acb3515a543e4cda9431167c55d546d00500 100644
--- a/grudge/models/advection.py
+++ b/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)
diff --git a/grudge/models/em.py b/grudge/models/em.py
index 38b73d3d2a598acef851586cd74c982761efc102..8960e0d3dbb9c3902f962848b5f0b960a60530b6 100644
--- a/grudge/models/em.py
+++ b/grudge/models/em.py
@@ -120,7 +120,7 @@ class MaxwellOperator(HyperbolicOperator):
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 @@ class MaxwellOperator(HyperbolicOperator):
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)
diff --git a/grudge/models/wave.py b/grudge/models/wave.py
index e101fd0443918c3cbb3b5f530a002ff2437578fd..fa42879670f24e5393b8fc0a4002bdf2d5f9227c 100644
--- a/grudge/models/wave.py
+++ b/grudge/models/wave.py
@@ -90,7 +90,7 @@ class WeakWaveOperator(HyperbolicOperator):
u = wtpair[0]
v = wtpair[1:]
actx = u.int.array_context
- normal = thaw(op.normal(self.dcoll, wtpair.dd), actx)
+ normal = thaw(self.dcoll.normal(wtpair.dd), actx)
central_flux_weak = -self.c*flat_obj_array(
np.dot(v.avg, normal),
@@ -133,7 +133,7 @@ class WeakWaveOperator(HyperbolicOperator):
neu_bc = flat_obj_array(neu_u, -neu_v)
# radiation BCs -------------------------------------------------------
- rad_normal = thaw(op.normal(dcoll, dd=self.radiation_tag), actx)
+ rad_normal = thaw(dcoll.normal(dd=self.radiation_tag), actx)
rad_u = op.project(dcoll, "vol", self.radiation_tag, u)
rad_v = op.project(dcoll, "vol", self.radiation_tag, v)
@@ -233,7 +233,7 @@ class VariableCoefficientWeakWaveOperator(HyperbolicOperator):
u = wtpair[1]
v = wtpair[2:]
actx = u.int.array_context
- normal = thaw(op.normal(self.dcoll, wtpair.dd), actx)
+ normal = thaw(self.dcoll.normal(wtpair.dd), actx)
flux_central_weak = -0.5 * flat_obj_array(
np.dot(v.int*c.int + v.ext*c.ext, normal),
@@ -282,7 +282,7 @@ class VariableCoefficientWeakWaveOperator(HyperbolicOperator):
neu_bc = flat_obj_array(neu_c, neu_u, -neu_v)
# radiation BCs -------------------------------------------------------
- rad_normal = thaw(op.normal(dcoll, dd=self.radiation_tag), actx)
+ rad_normal = thaw(dcoll.normal(dd=self.radiation_tag), actx)
rad_c = op.project(dcoll, "vol", self.radiation_tag, self.c)
rad_u = op.project(dcoll, "vol", self.radiation_tag, u)
diff --git a/grudge/op.py b/grudge/op.py
index 200c89324a6ba82d2c33f55e0c7bd79a80b48368..1ab764e08dfc28dd0a3e4a07330fe1fba5bf6200 100644
--- a/grudge/op.py
+++ b/grudge/op.py
@@ -1,20 +1,4 @@
"""
-Data transfer and geometry
-^^^^^^^^^^^^^^^^^^^^^^^^^^
-
-Projection and interpolation
-----------------------------
-
-.. autofunction:: project
-
-Geometric quantities
---------------------
-
-.. autofunction:: nodes
-.. autofunction:: normal
-.. autofunction:: h_max_from_volume
-.. autofunction:: h_min_from_volume
-
Core DG routines
^^^^^^^^^^^^^^^^
@@ -38,24 +22,6 @@ Mass, inverse mass, and face mass operators
.. autofunction:: mass
.. autofunction:: inverse_mass
.. autofunction:: face_mass
-
-Support functions
-^^^^^^^^^^^^^^^^^
-
-Nodal reductions
-----------------
-
-.. autofunction:: norm
-.. autofunction:: nodal_sum
-.. autofunction:: nodal_min
-.. autofunction:: nodal_max
-.. autofunction:: integral
-
-Elementwise reductions
-----------------------
-
-.. autofunction:: elementwise_sum
-.. autofunction:: elementwise_integral
"""
__copyright__ = """
@@ -84,9 +50,6 @@ THE SOFTWARE.
"""
-from numbers import Number
-from functools import reduce
-
from arraycontext import (
ArrayContext,
FirstAxisIsElementsTag,
@@ -95,19 +58,29 @@ from arraycontext import (
from grudge.discretization import DiscretizationCollection
-from pytools import (
- memoize_in,
- memoize_on_first_arg,
- keyed_memoize_in
-)
+from pytools import memoize_in, keyed_memoize_in
from pytools.obj_array import obj_array_vectorize, make_obj_array
from meshmode.dof_array import DOFArray
import numpy as np
+
import grudge.dof_desc as dof_desc
-from grudge.trace_pair import ( # noqa
+from grudge.interpolation import interp # noqa: F401
+from grudge.projection import project # noqa: F401
+
+from grudge.reductions import ( # noqa: F401
+ norm,
+ nodal_sum,
+ nodal_min,
+ nodal_max,
+ integral,
+ elementwise_sum,
+ elementwise_integral,
+)
+
+from grudge.trace_pair import ( # noqa: F401
interior_trace_pair,
interior_trace_pairs,
connected_ranks,
@@ -130,129 +103,6 @@ class HasElementwiseMatvecTag(FirstAxisIsElementsTag):
# }}}
-# {{{ Interpolation and projection
-
-# FIXME: Should reintroduce 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("using 'interp' is deprecated, use 'project' instead.",
-# DeprecationWarning, stacklevel=2)
-#
-# return project(dcoll, src, tgt, vec)
-
-
-def project(dcoll: DiscretizationCollection, src, tgt, vec):
- """Project from one discretization to another, e.g. from the
- volume to the boundary, or from the base to the an overintegrated
- quadrature discretization.
-
- :arg src: a :class:`~grudge.dof_desc.DOFDesc`, or a value convertible to one.
- :arg tgt: a :class:`~grudge.dof_desc.DOFDesc`, or a value convertible to one.
- :arg vec: a :class:`~meshmode.dof_array.DOFArray` or a
- :class:`~arraycontext.ArrayContainer`.
- """
- src = dof_desc.as_dofdesc(src)
- tgt = dof_desc.as_dofdesc(tgt)
- if src == tgt:
- return vec
-
- if isinstance(vec, np.ndarray):
- return obj_array_vectorize(
- lambda el: project(dcoll, src, tgt, el), vec)
-
- if isinstance(vec, Number):
- return vec
-
- return dcoll.connection_from_dds(src, tgt)(vec)
-
-# }}}
-
-
-# {{{ geometric properties
-
-def nodes(dcoll: DiscretizationCollection, dd=None) -> np.ndarray:
- r"""Return the nodes of a discretization.
-
- :arg dd: a :class:`~grudge.dof_desc.DOFDesc`, or a value convertible to one.
- Defaults to the base volume discretization.
- :returns: an object array of :class:`~meshmode.dof_array.DOFArray`\ s.
- """
- if dd is None:
- dd = dof_desc.DD_VOLUME
- dd = dof_desc.as_dofdesc(dd)
-
- return dcoll.nodes(dd)
-
-
-def normal(dcoll: DiscretizationCollection, dd) -> np.ndarray:
- r"""Get the unit normal to the specified surface discretization, *dd*.
-
- :arg dd: a :class:`~grudge.dof_desc.DOFDesc` as the surface discretization.
- :returns: an object array of :class:`~meshmode.dof_array.DOFArray`\ s.
- """
- return dcoll.normal(dd)
-
-
-@memoize_on_first_arg
-def h_max_from_volume(dcoll: DiscretizationCollection, dim=None, dd=None):
- """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.
- """
- if dd is None:
- dd = dof_desc.DD_VOLUME
- dd = dof_desc.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, mass(dcoll, dd, ones))
- ) ** (1.0 / dim)
-
-
-@memoize_on_first_arg
-def h_min_from_volume(dcoll: DiscretizationCollection, dim=None, dd=None):
- """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.
- """
- if dd is None:
- dd = dof_desc.DD_VOLUME
- dd = dof_desc.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, mass(dcoll, dd, ones))
- ) ** (1.0 / dim)
-
-# }}}
-
-
# {{{ Derivative operators
def reference_derivative_matrices(actx: ArrayContext, element_group):
@@ -1019,210 +869,4 @@ def face_mass(dcoll: DiscretizationCollection, *args):
# }}}
-# {{{ Nodal reductions
-
-def _norm(dcoll: DiscretizationCollection, vec, p, dd):
- if isinstance(vec, Number):
- return np.fabs(vec)
- if p == 2:
- return np.real_if_close(np.sqrt(
- nodal_sum(
- dcoll,
- dd,
- vec.conj() * _apply_mass_operator(dcoll, dd, dd, vec)
- )
- ))
- elif p == np.inf:
- return nodal_max(dcoll, dd, abs(vec))
- else:
- raise NotImplementedError("Unsupported value of p")
-
-
-def norm(dcoll: DiscretizationCollection, vec, p, dd=None):
- r"""Return the vector p-norm of a function represented
- by its vector of degrees of freedom *vec*.
-
- :arg vec: a :class:`~meshmode.dof_array.DOFArray` or an object array of
- a :class:`~meshmode.dof_array.DOFArray`\ s,
- where the last axis of the array must have length
- matching the volume dimension.
- :arg p: an integer denoting the order of the integral norm. Currently,
- only values of 2 or `numpy.inf` are supported.
- :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 nonegative scalar denoting the norm.
- """
- # FIXME: We should make the nodal reductions respect MPI.
- # See: https://github.com/inducer/grudge/issues/112 and
- # https://github.com/inducer/grudge/issues/113
-
- if dd is None:
- dd = dof_desc.DD_VOLUME
-
- dd = dof_desc.as_dofdesc(dd)
-
- if isinstance(vec, np.ndarray):
- if p == 2:
- return sum(
- norm(dcoll, vec[idx], p, dd=dd)**2
- for idx in np.ndindex(vec.shape))**0.5
- elif p == np.inf:
- return max(
- norm(dcoll, vec[idx], np.inf, dd=dd)
- for idx in np.ndindex(vec.shape))
- else:
- raise ValueError("unsupported norm order")
-
- return _norm(dcoll, vec, p, dd)
-
-
-def nodal_sum(dcoll: DiscretizationCollection, dd, vec):
- r"""Return the nodal sum of a vector of degrees of freedom *vec*.
-
- :arg dd: a :class:`~grudge.dof_desc.DOFDesc`, or a value
- convertible to one.
- :arg vec: a :class:`~meshmode.dof_array.DOFArray`.
- :returns: a scalar denoting the nodal sum.
- """
- # FIXME: We should make the nodal reductions respect MPI.
- # See: https://github.com/inducer/grudge/issues/112 and
- # https://github.com/inducer/grudge/issues/113
- actx = vec.array_context
- return sum([actx.np.sum(grp_ary) for grp_ary in vec])
-
-
-def nodal_min(dcoll: DiscretizationCollection, dd, vec):
- r"""Return the nodal minimum of a vector of degrees of freedom *vec*.
-
- :arg dd: a :class:`~grudge.dof_desc.DOFDesc`, or a value
- convertible to one.
- :arg vec: a :class:`~meshmode.dof_array.DOFArray`.
- :returns: a scalar denoting the nodal minimum.
- """
- # FIXME: We should make the nodal reductions respect MPI.
- # See: https://github.com/inducer/grudge/issues/112 and
- # https://github.com/inducer/grudge/issues/113
- actx = vec.array_context
- return reduce(lambda acc, grp_ary: actx.np.minimum(acc, actx.np.min(grp_ary)),
- vec, np.inf)
-
-
-def nodal_max(dcoll: DiscretizationCollection, dd, vec):
- r"""Return the nodal maximum of a vector of degrees of freedom *vec*.
-
- :arg dd: a :class:`~grudge.dof_desc.DOFDesc`, or a value
- convertible to one.
- :arg vec: a :class:`~meshmode.dof_array.DOFArray`.
- :returns: a scalar denoting the nodal maximum.
- """
- # FIXME: We should make the nodal reductions respect MPI.
- # See: https://github.com/inducer/grudge/issues/112 and
- # https://github.com/inducer/grudge/issues/113
- actx = vec.array_context
- return reduce(lambda acc, grp_ary: actx.np.maximum(acc, actx.np.max(grp_ary)),
- vec, -np.inf)
-
-
-def integral(dcoll: DiscretizationCollection, dd, vec):
- """Numerically integrates a function represented by a
- :class:`~meshmode.dof_array.DOFArray` of degrees of freedom.
-
- :arg dd: a :class:`~grudge.dof_desc.DOFDesc`, or a value convertible to one.
- :arg vec: a :class:`~meshmode.dof_array.DOFArray`
- :returns: a scalar denoting the evaluated integral.
- """
-
- dd = dof_desc.as_dofdesc(dd)
-
- ones = dcoll.discr_from_dd(dd).zeros(vec.array_context) + 1.0
- return nodal_sum(
- dcoll, dd, vec * _apply_mass_operator(dcoll, dd, dd, ones)
- )
-
-# }}}
-
-
-# {{{ Elementwise reductions
-
-def _map_elementwise_reduction(actx: ArrayContext, op_name):
- @memoize_in(actx, (_map_elementwise_reduction,
- "elementwise_%s_prg" % op_name))
- def prg():
- return make_loopy_program(
- [
- "{[iel]: 0 <= iel < nelements}",
- "{[idof, jdof]: 0 <= idof, jdof < ndofs}"
- ],
- """
- result[iel, idof] = %s(jdof, operand[iel, jdof])
- """ % op_name,
- name="grudge_elementwise_%s_knl" % op_name
- )
- return prg()
-
-
-def elementwise_sum(dcoll: DiscretizationCollection, *args):
- r"""Returns a vector of DOFs with all entries on each element set
- to the sum of DOFs on that element.
-
- May be called with ``(dcoll, vec)`` or ``(dcoll, dd, vec)``.
-
- :arg dd: a :class:`~grudge.dof_desc.DOFDesc`, or a value convertible to one.
- Defaults to the base volume discretization if not provided.
- :arg vec: a :class:`~meshmode.dof_array.DOFArray`
- :returns: a :class:`~meshmode.dof_array.DOFArray` whose entries
- denote the element-wise sum of *vec*.
- """
-
- if len(args) == 1:
- vec, = args
- dd = dof_desc.DOFDesc("vol", dof_desc.DISCR_TAG_BASE)
- elif len(args) == 2:
- dd, vec = args
- else:
- raise TypeError("invalid number of arguments")
-
- dd = dof_desc.as_dofdesc(dd)
-
- if isinstance(vec, np.ndarray):
- return obj_array_vectorize(
- lambda vi: elementwise_sum(dcoll, dd, vi), vec
- )
-
- actx = vec.array_context
-
- return DOFArray(
- actx,
- data=tuple(
- actx.call_loopy(
- _map_elementwise_reduction(actx, "sum"),
- operand=vec_i
- )["result"]
- for vec_i in vec
- )
- )
-
-
-def elementwise_integral(dcoll: DiscretizationCollection, dd, vec):
- """Numerically integrates a function represented by a
- :class:`~meshmode.dof_array.DOFArray` of degrees of freedom in
- each element of a discretization, given by *dd*.
-
- :arg dd: a :class:`~grudge.dof_desc.DOFDesc`, or a value convertible to one.
- :arg vec: a :class:`~meshmode.dof_array.DOFArray`
- :returns: a :class:`~meshmode.dof_array.DOFArray` containing the
- elementwise integral of *vec*, represented as a constant function on each
- cell.
- """
-
- dd = dof_desc.as_dofdesc(dd)
-
- ones = dcoll.discr_from_dd(dd).zeros(vec.array_context) + 1.0
- return elementwise_sum(
- dcoll, dd, vec * _apply_mass_operator(dcoll, dd, dd, ones)
- )
-
-# }}}
-
-
# vim: foldmethod=marker
diff --git a/grudge/projection.py b/grudge/projection.py
new file mode 100644
index 0000000000000000000000000000000000000000..3129253ca0b9c355334a452da4a60627dd5ea2cc
--- /dev/null
+++ b/grudge/projection.py
@@ -0,0 +1,68 @@
+"""
+.. currentmodule:: grudge.op
+
+Projections
+-----------
+
+.. autofunction:: project
+"""
+
+__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.
+"""
+
+
+import numpy as np
+
+from grudge.discretization import DiscretizationCollection
+from grudge.dof_desc import as_dofdesc
+
+from numbers import Number
+
+from pytools.obj_array import obj_array_vectorize
+
+
+def project(dcoll: DiscretizationCollection, src, tgt, vec):
+ """Project from one discretization to another, e.g. from the
+ volume to the boundary, or from the base to the an overintegrated
+ quadrature discretization.
+
+ :arg src: a :class:`~grudge.dof_desc.DOFDesc`, or a value convertible to one.
+ :arg tgt: a :class:`~grudge.dof_desc.DOFDesc`, or a value convertible to one.
+ :arg vec: a :class:`~meshmode.dof_array.DOFArray` or a
+ :class:`~arraycontext.ArrayContainer`.
+ """
+ src = as_dofdesc(src)
+ tgt = as_dofdesc(tgt)
+
+ if src == tgt:
+ return vec
+
+ if isinstance(vec, np.ndarray):
+ return obj_array_vectorize(
+ lambda el: project(dcoll, src, tgt, el), vec)
+
+ if isinstance(vec, Number):
+ return vec
+
+ return dcoll.connection_from_dds(src, tgt)(vec)
diff --git a/grudge/reductions.py b/grudge/reductions.py
new file mode 100644
index 0000000000000000000000000000000000000000..b8a8a868bf2549041205bb4b603d4033ab38f9f4
--- /dev/null
+++ b/grudge/reductions.py
@@ -0,0 +1,272 @@
+"""
+.. currentmodule:: grudge.op
+
+Nodal reductions
+----------------
+
+.. autofunction:: norm
+.. autofunction:: nodal_sum
+.. autofunction:: nodal_min
+.. autofunction:: nodal_max
+.. autofunction:: integral
+
+Elementwise reductions
+----------------------
+
+.. autofunction:: elementwise_sum
+.. autofunction:: elementwise_integral
+"""
+
+__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 numbers import Number
+from functools import reduce
+
+from arraycontext import (
+ ArrayContext,
+ make_loopy_program
+)
+
+from grudge.discretization import DiscretizationCollection
+
+from pytools import memoize_in
+from pytools.obj_array import obj_array_vectorize
+
+from meshmode.dof_array import DOFArray
+
+import numpy as np
+import grudge.dof_desc as dof_desc
+
+
+# {{{ Nodal reductions
+
+def _norm(dcoll: DiscretizationCollection, vec, p, dd):
+ if isinstance(vec, Number):
+ return np.fabs(vec)
+ if p == 2:
+ from grudge.op import _apply_mass_operator
+ return np.real_if_close(np.sqrt(
+ nodal_sum(
+ dcoll,
+ dd,
+ vec.conj() * _apply_mass_operator(dcoll, dd, dd, vec)
+ )
+ ))
+ elif p == np.inf:
+ return nodal_max(dcoll, dd, abs(vec))
+ else:
+ raise NotImplementedError("Unsupported value of p")
+
+
+def norm(dcoll: DiscretizationCollection, vec, p, dd=None) -> float:
+ r"""Return the vector p-norm of a function represented
+ by its vector of degrees of freedom *vec*.
+
+ :arg vec: a :class:`~meshmode.dof_array.DOFArray` or an object array of
+ a :class:`~meshmode.dof_array.DOFArray`\ s,
+ where the last axis of the array must have length
+ matching the volume dimension.
+ :arg p: an integer denoting the order of the integral norm. Currently,
+ only values of 2 or `numpy.inf` are supported.
+ :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 nonegative scalar denoting the norm.
+ """
+ # FIXME: We should make the nodal reductions respect MPI.
+ # See: https://github.com/inducer/grudge/issues/112 and
+ # https://github.com/inducer/grudge/issues/113
+
+ if dd is None:
+ dd = dof_desc.DD_VOLUME
+
+ dd = dof_desc.as_dofdesc(dd)
+
+ if isinstance(vec, np.ndarray):
+ if p == 2:
+ return sum(
+ norm(dcoll, vec[idx], p, dd=dd)**2
+ for idx in np.ndindex(vec.shape))**0.5
+ elif p == np.inf:
+ return max(
+ norm(dcoll, vec[idx], np.inf, dd=dd)
+ for idx in np.ndindex(vec.shape))
+ else:
+ raise ValueError("unsupported norm order")
+
+ return _norm(dcoll, vec, p, dd)
+
+
+def nodal_sum(dcoll: DiscretizationCollection, dd, vec) -> float:
+ r"""Return the nodal sum of a vector of degrees of freedom *vec*.
+
+ :arg dd: a :class:`~grudge.dof_desc.DOFDesc`, or a value
+ convertible to one.
+ :arg vec: a :class:`~meshmode.dof_array.DOFArray`.
+ :returns: a scalar denoting the nodal sum.
+ """
+ # FIXME: We should make the nodal reductions respect MPI.
+ # See: https://github.com/inducer/grudge/issues/112 and
+ # https://github.com/inducer/grudge/issues/113
+ actx = vec.array_context
+ return sum([actx.np.sum(grp_ary) for grp_ary in vec])
+
+
+def nodal_min(dcoll: DiscretizationCollection, dd, vec) -> float:
+ r"""Return the nodal minimum of a vector of degrees of freedom *vec*.
+
+ :arg dd: a :class:`~grudge.dof_desc.DOFDesc`, or a value
+ convertible to one.
+ :arg vec: a :class:`~meshmode.dof_array.DOFArray`.
+ :returns: a scalar denoting the nodal minimum.
+ """
+ # FIXME: We should make the nodal reductions respect MPI.
+ # See: https://github.com/inducer/grudge/issues/112 and
+ # https://github.com/inducer/grudge/issues/113
+ actx = vec.array_context
+ return reduce(lambda acc, grp_ary: actx.np.minimum(acc, actx.np.min(grp_ary)),
+ vec, np.inf)
+
+
+def nodal_max(dcoll: DiscretizationCollection, dd, vec) -> float:
+ r"""Return the nodal maximum of a vector of degrees of freedom *vec*.
+
+ :arg dd: a :class:`~grudge.dof_desc.DOFDesc`, or a value
+ convertible to one.
+ :arg vec: a :class:`~meshmode.dof_array.DOFArray`.
+ :returns: a scalar denoting the nodal maximum.
+ """
+ # FIXME: We should make the nodal reductions respect MPI.
+ # See: https://github.com/inducer/grudge/issues/112 and
+ # https://github.com/inducer/grudge/issues/113
+ actx = vec.array_context
+ return reduce(lambda acc, grp_ary: actx.np.maximum(acc, actx.np.max(grp_ary)),
+ vec, -np.inf)
+
+
+def integral(dcoll: DiscretizationCollection, dd, vec) -> float:
+ """Numerically integrates a function represented by a
+ :class:`~meshmode.dof_array.DOFArray` of degrees of freedom.
+
+ :arg dd: a :class:`~grudge.dof_desc.DOFDesc`, or a value convertible to one.
+ :arg vec: a :class:`~meshmode.dof_array.DOFArray`
+ :returns: a scalar denoting the evaluated integral.
+ """
+ from grudge.op import _apply_mass_operator
+
+ dd = dof_desc.as_dofdesc(dd)
+
+ ones = dcoll.discr_from_dd(dd).zeros(vec.array_context) + 1.0
+ return nodal_sum(
+ dcoll, dd, vec * _apply_mass_operator(dcoll, dd, dd, ones)
+ )
+
+# }}}
+
+
+# {{{ Elementwise reductions
+
+def _map_elementwise_reduction(actx: ArrayContext, op_name):
+ @memoize_in(actx, (_map_elementwise_reduction,
+ "elementwise_%s_prg" % op_name))
+ def prg():
+ return make_loopy_program(
+ [
+ "{[iel]: 0 <= iel < nelements}",
+ "{[idof, jdof]: 0 <= idof, jdof < ndofs}"
+ ],
+ """
+ result[iel, idof] = %s(jdof, operand[iel, jdof])
+ """ % op_name,
+ name="grudge_elementwise_%s_knl" % op_name
+ )
+ return prg()
+
+
+def elementwise_sum(dcoll: DiscretizationCollection, *args) -> DOFArray:
+ r"""Returns a vector of DOFs with all entries on each element set
+ to the sum of DOFs on that element.
+
+ May be called with ``(dcoll, vec)`` or ``(dcoll, dd, vec)``.
+
+ :arg dd: a :class:`~grudge.dof_desc.DOFDesc`, or a value convertible to one.
+ Defaults to the base volume discretization if not provided.
+ :arg vec: a :class:`~meshmode.dof_array.DOFArray`
+ :returns: a :class:`~meshmode.dof_array.DOFArray` whose entries
+ denote the element-wise sum of *vec*.
+ """
+
+ if len(args) == 1:
+ vec, = args
+ dd = dof_desc.DOFDesc("vol", dof_desc.DISCR_TAG_BASE)
+ elif len(args) == 2:
+ dd, vec = args
+ else:
+ raise TypeError("invalid number of arguments")
+
+ dd = dof_desc.as_dofdesc(dd)
+
+ if isinstance(vec, np.ndarray):
+ return obj_array_vectorize(
+ lambda vi: elementwise_sum(dcoll, dd, vi), vec
+ )
+
+ actx = vec.array_context
+
+ return DOFArray(
+ actx,
+ data=tuple(
+ actx.call_loopy(
+ _map_elementwise_reduction(actx, "sum"),
+ operand=vec_i
+ )["result"]
+ for vec_i in vec
+ )
+ )
+
+
+def elementwise_integral(dcoll: DiscretizationCollection, dd, vec) -> DOFArray:
+ """Numerically integrates a function represented by a
+ :class:`~meshmode.dof_array.DOFArray` of degrees of freedom in
+ each element of a discretization, given by *dd*.
+
+ :arg dd: a :class:`~grudge.dof_desc.DOFDesc`, or a value convertible to one.
+ :arg vec: a :class:`~meshmode.dof_array.DOFArray`
+ :returns: a :class:`~meshmode.dof_array.DOFArray` containing the
+ elementwise integral if *vec*.
+ """
+ from grudge.op import _apply_mass_operator
+
+ dd = dof_desc.as_dofdesc(dd)
+
+ ones = dcoll.discr_from_dd(dd).zeros(vec.array_context) + 1.0
+ return elementwise_sum(
+ dcoll, dd, vec * _apply_mass_operator(dcoll, dd, dd, ones)
+ )
+
+# }}}
+
+
+# vim: foldmethod=marker
diff --git a/grudge/trace_pair.py b/grudge/trace_pair.py
index ad2e9c815f6b38d48f7dabf95314bdfceda553fc..a0667d056d3023baeecdec6a3240095489830eb4 100644
--- a/grudge/trace_pair.py
+++ b/grudge/trace_pair.py
@@ -1,3 +1,25 @@
+"""
+Trace Pairs
+^^^^^^^^^^^
+
+Container class
+---------------
+
+.. autoclass:: TracePair
+
+Boundary trace functions
+------------------------
+
+.. autofunction:: bdry_trace_pair
+.. autofunction:: bv_trace_pair
+
+Interior and cross-rank trace functions
+---------------------------------------
+
+.. autofunction:: interior_trace_pairs
+.. autofunction:: cross_rank_trace_pairs
+"""
+
__copyright__ = """
Copyright (C) 2021 University of Illinois Board of Trustees
"""
@@ -45,29 +67,6 @@ import numpy as np
import grudge.dof_desc as dof_desc
-__doc__ = """
-Trace Pairs
-^^^^^^^^^^^
-
-Container class
----------------
-
-.. autoclass:: TracePair
-
-Boundary trace functions
-------------------------
-
-.. autofunction:: bdry_trace_pair
-.. autofunction:: bv_trace_pair
-
-Interior and cross-rank trace functions
----------------------------------------
-
-.. autofunction:: interior_trace_pairs
-.. autofunction:: cross_rank_trace_pairs
-"""
-
-
# {{{ Trace pair container class
@with_container_arithmetic(
diff --git a/test/test_grudge.py b/test/test_grudge.py
index 4963642ecf46cb300f719686f87b7e63049afd73..522028fe3f0d449e5f112e4dd29ef9c302e8351b 100644
--- a/test/test_grudge.py
+++ b/test/test_grudge.py
@@ -250,7 +250,9 @@ def test_mass_surface_area(actx_factory, name):
# }}}
# compute max element size
- h_max = op.h_max_from_volume(dcoll)
+ from grudge.dt_utils import h_max_from_volume
+
+ h_max = h_max_from_volume(dcoll)
eoc.add_data_point(h_max, area_error)
@@ -321,7 +323,9 @@ def test_mass_operator_inverse(actx_factory, name):
# }}}
# compute max element size
- h_max = op.h_max_from_volume(dcoll)
+ from grudge.dt_utils import h_max_from_volume
+
+ h_max = h_max_from_volume(dcoll)
eoc.add_data_point(h_max, inv_error)
@@ -378,7 +382,7 @@ def test_face_normal_surface(actx_factory, mesh_name):
)
surf_normal = surf_normal / actx.np.sqrt(sum(surf_normal**2))
- face_normal_i = thaw(op.normal(dcoll, df), actx)
+ face_normal_i = thaw(dcoll.normal(df), actx)
face_normal_e = dcoll.opposite_face_connection()(face_normal_i)
if mesh.ambient_dim == 3:
@@ -498,7 +502,7 @@ def test_2d_gauss_theorem(actx_factory):
int_1 = op.integral(dcoll, "vol", op.local_div(dcoll, f_volm))
prj_f = op.project(dcoll, "vol", BTAG_ALL, f_volm)
- normal = thaw(op.normal(dcoll, BTAG_ALL), actx)
+ normal = thaw(dcoll.normal(BTAG_ALL), actx)
int_2 = op.integral(dcoll, BTAG_ALL, prj_f.dot(normal))
assert abs(int_1 - int_2) < 1e-13
@@ -600,14 +604,14 @@ def test_surface_divergence_theorem(actx_factory, mesh_name, visualize=False):
ambient_dim = dcoll.ambient_dim
# variables
- f_num = f(thaw(op.nodes(dcoll, dd=dd), actx))
- f_quad_num = f(thaw(op.nodes(dcoll, dd=dq), actx))
+ f_num = f(thaw(dcoll.nodes(dd=dd), actx))
+ f_quad_num = f(thaw(dcoll.nodes(dd=dq), actx))
from grudge.geometry import normal, summed_curvature
kappa = summed_curvature(actx, dcoll, dd=dq)
normal = normal(actx, dcoll, dd=dq)
- face_normal = thaw(op.normal(dcoll, df), actx)
+ face_normal = thaw(dcoll.normal(df), actx)
face_f = op.project(dcoll, dd, df, f_num)
# operators
@@ -627,7 +631,9 @@ 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
- h_max = op.h_max_from_volume(dcoll)
+ from grudge.dt_utils import h_max_from_volume
+
+ h_max = h_max_from_volume(dcoll)
eoc_global.add_data_point(h_max, err_global)
eoc_local.add_data_point(h_max, err_local)
@@ -747,12 +753,12 @@ def test_convergence_advec(actx_factory, mesh_name, mesh_pars, op_type, flux_typ
"weak": WeakAdvectionOperator}[op_type]
adv_operator = op_class(dcoll, v,
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):
@@ -763,7 +769,9 @@ def test_convergence_advec(actx_factory, mesh_name, mesh_pars, op_type, flux_typ
else:
final_time = 0.2
- h_max = op.h_max_from_volume(dcoll, dim=dcoll.ambient_dim)
+ from grudge.dt_utils import h_max_from_volume
+
+ h_max = h_max_from_volume(dcoll, dim=dcoll.ambient_dim)
dt = dt_factor * h_max/order**2
nsteps = (final_time // dt) + 1
dt = final_time/nsteps + 1e-15
@@ -842,7 +850,7 @@ def test_convergence_maxwell(actx_factory, order):
def analytic_sol(x, t=0):
return get_rectangular_cavity_mode(actx, x, t, 1, (1, 2, 2))
- nodes = thaw(op.nodes(dcoll), actx)
+ nodes = thaw(dcoll.nodes(), actx)
fields = analytic_sol(nodes, t=0)
from grudge.models.em import MaxwellOperator
@@ -939,7 +947,7 @@ def test_improvement_quadrature(actx_factory, order):
discr_tag_to_group_factory=discr_tag_to_group_factory
)
- nodes = thaw(op.nodes(dcoll), actx)
+ nodes = thaw(dcoll.nodes(), actx)
def zero_inflow(dtag, t=0):
dd = dof_desc.DOFDesc(dtag, qtag)
@@ -988,7 +996,7 @@ def test_bessel(actx_factory):
dcoll = DiscretizationCollection(actx, mesh, order=3)
- nodes = thaw(op.nodes(dcoll), actx)
+ nodes = thaw(dcoll.nodes(), actx)
r = actx.np.sqrt(nodes[0]**2 + nodes[1]**2)
# FIXME: Bessel functions need to brought out of the symbolic
@@ -1079,50 +1087,13 @@ def test_empty_boundary(actx_factory):
a=(-0.5,)*dim, b=(0.5,)*dim,
nelements_per_axis=(8,)*dim, order=4)
dcoll = DiscretizationCollection(actx, mesh, order=4)
- normal = op.normal(dcoll, BTAG_NONE)
+ normal = dcoll.normal(BTAG_NONE)
from meshmode.dof_array import DOFArray
for component in normal:
assert isinstance(component, DOFArray)
assert len(component) == len(dcoll.discr_from_dd(BTAG_NONE).groups)
-# {{{ DiscretizationCollection testing
-
-def test_dcoll_nodes_and_normals(actx_factory):
- actx = actx_factory()
-
- from meshmode.mesh import BTAG_ALL
- from meshmode.mesh.generation import generate_warped_rect_mesh
-
- mesh = generate_warped_rect_mesh(dim=2, order=4, nelements_side=6)
-
- dcoll = DiscretizationCollection(actx, mesh, order=3)
-
- # Nodes should be *indentical*
- nodes = thaw(op.nodes(dcoll), actx)
- dcoll_nodes = thaw(dcoll.nodes(), actx)
-
- assert op.norm(dcoll, nodes - dcoll_nodes, np.inf) == 0.0
-
- nodes_btag = thaw(op.nodes(dcoll, BTAG_ALL), actx)
- dcoll_nodes_btag = thaw(dcoll.nodes(BTAG_ALL), actx)
-
- assert op.norm(dcoll, nodes_btag - dcoll_nodes_btag, np.inf) == 0.0
-
- # Normals should be *indentical*
- normals = thaw(op.normal(dcoll, "int_faces"), actx)
- dcoll_normals = thaw(dcoll.normal("int_faces"), actx)
-
- assert op.norm(dcoll, normals - dcoll_normals, np.inf) == 0.0
-
- normals_btag = thaw(op.normal(dcoll, BTAG_ALL), actx)
- dcoll_normals_btag = thaw(dcoll.normal(BTAG_ALL), actx)
-
- assert op.norm(dcoll, normals_btag - dcoll_normals_btag, np.inf) == 0.0
-
-# }}}
-
-
# You can test individual routines by typing
# $ python test_grudge.py 'test_routine()'
diff --git a/test/test_mpi_communication.py b/test/test_mpi_communication.py
index 879d3bf87dc3d4ad8dfeb1c3512f1423ec1589f7..4e3484e9c1bf9c1335761001e872c97eda3907d9 100644
--- a/test/test_mpi_communication.py
+++ b/test/test_mpi_communication.py
@@ -76,7 +76,7 @@ def simple_mpi_communication_entrypoint():
dcoll = DiscretizationCollection(actx, local_mesh, order=5,
mpi_communicator=comm)
- x = thaw(op.nodes(dcoll), actx)
+ x = thaw(dcoll.nodes(), actx)
myfunc = actx.np.sin(np.dot(x, [2, 3]))
from grudge.dof_desc import as_dofdesc
@@ -147,7 +147,7 @@ def mpi_communication_entrypoint():
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)]
)
diff --git a/test/test_op.py b/test/test_op.py
index 2d12409ee3a20c8cfce8c02a539dbd7e31f9e2ed..6aa3ba7915ec4767c9477bde6fb8355aa48ff4ec 100644
--- a/test/test_op.py
+++ b/test/test_op.py
@@ -88,7 +88,7 @@ def test_gradient(actx_factory, form, dim, order, vectorize, nested,
result[dim-1] = result[dim-1] * (-np.pi/2*actx.np.sin(np.pi/2*x[dim-1]))
return result
- x = thaw(op.nodes(dcoll), actx)
+ x = thaw(dcoll.nodes(), actx)
if vectorize:
u = make_obj_array([(i+1)*f(x) for i in range(dim)])
@@ -98,7 +98,7 @@ def test_gradient(actx_factory, form, dim, order, vectorize, nested,
def get_flux(u_tpair):
dd = u_tpair.dd
dd_allfaces = dd.with_dtag("all_faces")
- normal = thaw(op.normal(dcoll, dd), actx)
+ normal = thaw(dcoll.normal(dd), actx)
u_avg = u_tpair.avg
if vectorize:
if nested:
@@ -212,7 +212,7 @@ def test_divergence(actx_factory, form, dim, order, vectorize, nested,
result = result + deriv
return result
- x = thaw(op.nodes(dcoll), actx)
+ x = thaw(dcoll.nodes(), actx)
if vectorize:
u = make_obj_array([(i+1)*f(x) for i in range(dim)])
@@ -224,7 +224,7 @@ def test_divergence(actx_factory, form, dim, order, vectorize, nested,
def get_flux(u_tpair):
dd = u_tpair.dd
dd_allfaces = dd.with_dtag("all_faces")
- normal = thaw(op.normal(dcoll, dd), actx)
+ normal = thaw(dcoll.normal(dd), actx)
flux = u_tpair.avg @ normal
return op.project(dcoll, dd, dd_allfaces, flux)