diff --git a/doc/conf.py b/doc/conf.py
index 3becabe61a7206322542e48a615be62742af0dae..560efec5951aa5de9a4b732e304a121c0c0c8baa 100644
--- a/doc/conf.py
+++ b/doc/conf.py
@@ -320,11 +320,11 @@ texinfo_documents = [
 # Example configuration for intersphinx: refer to the Python standard library.
 intersphinx_mapping = {
     'http://docs.python.org/': None,
-    'http://documen.tician.de/boxtree/': None,
     'http://docs.scipy.org/doc/numpy/': None,
-    'http://documen.tician.de/modepy/': None,
     'http://documen.tician.de/pyopencl/': None,
+    'http://documen.tician.de/modepy/': None,
     'http://documen.tician.de/pymbolic/': None,
+    'http://documen.tician.de/meshmode/': None,
     #'http://documen.tician.de/loopy/': None,
     }
 autoclass_content = "both"
diff --git a/doc/index.rst b/doc/index.rst
index a20ec5093527240204a0da00b19b4188f08317c2..b996c470ffc998618b118e815b05b1c42ff7e9be 100644
--- a/doc/index.rst
+++ b/doc/index.rst
@@ -1,9 +1,4 @@
-.. grudge documentation master file, created by
-   sphinx-quickstart on Sun Sep 27 13:08:30 2015.
-   You can adapt this file completely to your liking, but it should at least
-   contain the root `toctree` directive.
-
-Welcome to grudge's documentation!
+Welcome to grudge's Documentation!
 ==================================
 
 Contents:
@@ -12,10 +7,10 @@ Contents:
     :maxdepth: 2
 
     misc
+    symbolic
 
 
-
-Indices and tables
+Indices and Tables
 ==================
 
 * :ref:`genindex`
diff --git a/doc/symbolic.rst b/doc/symbolic.rst
new file mode 100644
index 0000000000000000000000000000000000000000..632b1f861beaa80c3f3c26e93b08272ec4a31259
--- /dev/null
+++ b/doc/symbolic.rst
@@ -0,0 +1,14 @@
+Symbolic Operator Representation
+================================
+
+Based on :mod:`pymbolic`.
+
+Basic Objects
+-------------
+
+.. automodule:: grudge.symbolic.primitives
+
+Operators
+---------
+
+.. automodule:: grudge.symbolic.operators
diff --git a/examples/dagrt-fusion.py b/examples/dagrt-fusion.py
index 27e35abd82161e528481344eb28a6f0f81c22839..b9ef892e675ec70b9df309541013ebc89dd36878 100755
--- a/examples/dagrt-fusion.py
+++ b/examples/dagrt-fusion.py
@@ -50,6 +50,7 @@ from pymbolic.mapper import Mapper
 from pymbolic.mapper.evaluator import EvaluationMapper \
         as PymbolicEvaluationMapper
 from pytools import memoize
+from pytools.obj_array import join_fields
 
 from grudge import sym, bind, DGDiscretizationWithBoundaries
 from leap.rk import LSRK4MethodBuilder
@@ -245,7 +246,6 @@ class RK4TimeStepperBase(object):
         self.component_getter = component_getter
 
     def get_initial_context(self, fields, t_start, dt):
-        from pytools.obj_array import join_fields
 
         # Flatten fields.
         flattened_fields = []
@@ -286,7 +286,6 @@ class RK4TimeStepperBase(object):
         assert len(output_states) == num_fields
         assert len(output_states) == len(output_residuals)
 
-        from pytools.obj_array import join_fields
         flattened_results = join_fields(output_t, output_dt, *output_states)
 
         self.bound_op = bind(
@@ -354,8 +353,6 @@ class RK4TimeStepper(RK4TimeStepperBase):
                     + tuple(
                         Variable(field_var_name, dd=sym.DD_VOLUME)[i]
                         for i in range(num_fields)))))
-
-        from pytools.obj_array import join_fields
         sym_rhs = join_fields(*(call[i] for i in range(num_fields)))
 
         self.queue = queue
@@ -377,7 +374,6 @@ class RK4TimeStepper(RK4TimeStepperBase):
         self.component_getter = component_getter
 
     def _bound_op(self, queue, t, *args, profile_data=None):
-        from pytools.obj_array import join_fields
         context = {
                 "t": t,
                 self.field_var_name: join_fields(*args)}
@@ -452,16 +448,15 @@ def dg_flux(c, tpair):
     dims = len(v)
 
     normal = sym.normal(tpair.dd, dims)
-
-    flux_weak = sym.join_fields(
+    flux_weak = join_fields(
             np.dot(v.avg, normal),
             u.avg * normal)
 
-    flux_weak -= (1 if c > 0 else -1)*sym.join_fields(
+    flux_weak -= (1 if c > 0 else -1)*join_fields(
             0.5*(u.int-u.ext),
             0.5*(normal * np.dot(normal, v.int-v.ext)))
 
-    flux_strong = sym.join_fields(
+    flux_strong = join_fields(
             np.dot(v.int, normal),
             u.int * normal) - flux_weak
 
@@ -507,13 +502,13 @@ def get_strong_wave_op_with_discr_direct(cl_ctx, dims=2, order=4):
     rad_u = sym.cse(sym.interp("vol", BTAG_ALL)(u))
     rad_v = sym.cse(sym.interp("vol", BTAG_ALL)(v))
 
-    rad_bc = sym.cse(sym.join_fields(
+    rad_bc = sym.cse(join_fields(
             0.5*(rad_u - sign*np.dot(rad_normal, rad_v)),
             0.5*rad_normal*(np.dot(rad_normal, rad_v) - sign*rad_u)
             ), "rad_bc")
 
     sym_operator = (
-            - sym.join_fields(
+            - join_fields(
                 -c*np.dot(sym.nabla(dims), v) - source_f,
                 -c*(sym.nabla(dims)*u)
                 )
@@ -550,7 +545,6 @@ def test_stepper_equivalence(ctx_factory, order=4):
     elif dims == 3:
         dt = 0.02
 
-    from pytools.obj_array import join_fields
     ic = join_fields(discr.zeros(queue),
             [discr.zeros(queue) for i in range(discr.dim)])
 
@@ -804,7 +798,6 @@ def test_stepper_mem_ops(ctx_factory, use_fusion):
     dt = 0.04
     t_end = 0.2
 
-    from pytools.obj_array import join_fields
     ic = join_fields(discr.zeros(queue),
             [discr.zeros(queue) for i in range(discr.dim)])
 
@@ -975,7 +968,6 @@ def test_stepper_timing(ctx_factory, use_fusion):
     dt = 0.04
     t_end = 0.1
 
-    from pytools.obj_array import join_fields
     ic = join_fields(discr.zeros(queue),
             [discr.zeros(queue) for i in range(discr.dim)])
 
@@ -1040,7 +1032,6 @@ def get_example_stepper(queue, dims=2, order=3, use_fusion=True,
                 exec_mapper_factory=exec_mapper_factory)
 
     if return_ic:
-        from pytools.obj_array import join_fields
         ic = join_fields(discr.zeros(queue),
                 [discr.zeros(queue) for i in range(discr.dim)])
         return stepper, ic
diff --git a/grudge/__init__.py b/grudge/__init__.py
index b854007a48c227d58541a13a972a3a7e25c4e062..17358dd474304bec731bedb58a98cd692ec69935 100644
--- a/grudge/__init__.py
+++ b/grudge/__init__.py
@@ -22,6 +22,11 @@ OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
 THE SOFTWARE.
 """
 
-import grudge.sym  # noqa
-from grudge.execution import bind  # noqa
-from grudge.discretization import DGDiscretizationWithBoundaries  # noqa
+import grudge.symbolic as sym
+from grudge.execution import bind
+
+from grudge.discretization import DGDiscretizationWithBoundaries
+
+__all__ = [
+    "sym", "bind", "DGDiscretizationWithBoundaries"
+]
diff --git a/grudge/models/em.py b/grudge/models/em.py
index 3422c3be14f8ec10d872d46adc29caaa155131dd..225b3b4a411beae96359df0548ff5366f00ea675 100644
--- a/grudge/models/em.py
+++ b/grudge/models/em.py
@@ -437,7 +437,7 @@ def get_rectangular_cavity_mode(E_0, mode_indices):  # noqa: N803
     if dims == 2:
         tfac = sym.ScalarVariable("t") * omega
 
-        result = sym.join_fields(
+        result = join_fields(
             0,
             0,
             sym.sin(kx * x) * sym.sin(ky * y) * sym.cos(tfac),  # ez
@@ -449,7 +449,7 @@ def get_rectangular_cavity_mode(E_0, mode_indices):  # noqa: N803
         tdep = sym.exp(-1j * omega * sym.ScalarVariable("t"))
 
         gamma_squared = ky**2 + kx**2
-        result = sym.join_fields(
+        result = join_fields(
             -kx * kz * E_0*cx*sy*sz*tdep / gamma_squared,  # ex
             -ky * kz * E_0*sx*cy*sz*tdep / gamma_squared,  # ey
             E_0 * sx*sy*cz*tdep,  # ez
diff --git a/grudge/sym.py b/grudge/sym.py
deleted file mode 100644
index bf9c7de615c82c11a5a6d91fa9ed872341d0d32c..0000000000000000000000000000000000000000
--- a/grudge/sym.py
+++ /dev/null
@@ -1,28 +0,0 @@
-from __future__ import division, absolute_import
-
-__copyright__ = "Copyright (C) 2015 Andreas Kloeckner"
-
-__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.symbolic.primitives import *  # noqa
-from grudge.symbolic.operators import *  # noqa
-
-from grudge.symbolic.tools import pretty  # noqa
diff --git a/grudge/symbolic/__init__.py b/grudge/symbolic/__init__.py
index acd34e3532cba7516fab14dc1a9d2af46ddd649d..59fa5b26b5943a24982db5616632cef7dc8c042f 100644
--- a/grudge/symbolic/__init__.py
+++ b/grudge/symbolic/__init__.py
@@ -23,3 +23,9 @@ 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.symbolic.primitives import *        # noqa: F401,F403
+from grudge.symbolic.operators import *         # noqa: F401,F403
+from grudge.symbolic.tools import pretty        # noqa: F401
+
+from pymbolic.primitives import If, Comparison  # noqa: F401
diff --git a/grudge/symbolic/operators.py b/grudge/symbolic/operators.py
index a1899ed21ef6431f7b57a81fcf143680433b7b53..c1360b1627a79ccafe3e52405a099f550012a410 100644
--- a/grudge/symbolic/operators.py
+++ b/grudge/symbolic/operators.py
@@ -1,5 +1,3 @@
-"""Building blocks and mappers for operator expression trees."""
-
 from __future__ import division, absolute_import
 
 __copyright__ = "Copyright (C) 2008-2017 Andreas Kloeckner, Bogdan Enache"
@@ -27,16 +25,63 @@ THE SOFTWARE.
 from six.moves import intern
 
 import numpy as np
-import numpy.linalg as la  # noqa
 import pymbolic.primitives
 
+__doc__ = """
+
+Building blocks and mappers for operator expression trees.
+
+Basic Operators
+^^^^^^^^^^^^^^^
+
+.. autoclass:: Operator
+.. autoclass:: ElementwiseLinearOperator
+.. autoclass:: InterpolationOperator
+
+.. data:: interp
+
+Reductions
+^^^^^^^^^^
+
+.. autoclass:: ElementwiseMaxOperator
+
+.. autoclass:: NodalReductionOperator
+.. autoclass:: NodalSum
+.. autoclass:: NodalMax
+.. autoclass:: NodalMin
+
+Differentiation
+^^^^^^^^^^^^^^^
+
+.. autoclass:: StrongFormDiffOperatorBase
+.. autoclass:: WeakFormDiffOperatorBase
+.. autoclass:: StiffnessOperator
+.. autoclass:: DiffOperator
+.. autoclass:: StiffnessTOperator
+.. autoclass:: MInvSTOperator
+
+.. autoclass:: RefDiffOperator
+.. autoclass:: RefStiffnessTOperator
+
+.. autofunction:: nabla
+.. autofunction:: minv_stiffness_t
+.. autofunction:: stiffness
+.. autofunction:: stiffness_t
+
+Mass Operators
+^^^^^^^^^^^^^^
 
-def _sym():
-    # A hack to make referring to grudge.sym possible without
-    # circular imports and tons of typing.
+.. autoclass:: MassOperatorBase
 
-    from grudge import sym
-    return sym
+.. autoclass:: MassOperator
+.. autoclass:: InverseMassOperator
+.. autoclass:: FaceMassOperator
+
+.. autoclass:: RefMassOperator
+.. autoclass:: RefInverseMassOperator
+.. autoclass:: RefFaceMassOperator
+
+"""
 
 
 # {{{ base classes
@@ -45,18 +90,19 @@ class Operator(pymbolic.primitives.Expression):
     """
     .. attribute:: dd_in
 
-        an instance of :class:`grudge.sym.DOFDesc` describing the
+        an instance of :class:`~grudge.symbolic.primitives.DOFDesc` describing the
         input discretization.
 
     .. attribute:: dd_out
 
-        an instance of :class:`grudge.sym.DOFDesc` describing the
+        an instance of :class:`~grudge.symbolic.primitives.DOFDesc` describing the
         output discretization.
     """
 
     def __init__(self, dd_in, dd_out):
-        self.dd_in = _sym().as_dofdesc(dd_in)
-        self.dd_out = _sym().as_dofdesc(dd_out)
+        import grudge.symbolic.primitives as prim
+        self.dd_in = prim.as_dofdesc(dd_in)
+        self.dd_out = prim.as_dofdesc(dd_out)
 
     def stringifier(self):
         from grudge.symbolic.mappers import StringifyMapper
@@ -100,8 +146,10 @@ class ElementwiseLinearOperator(Operator):
 
 class InterpolationOperator(Operator):
     def __init__(self, dd_in, dd_out):
-        official_dd_in = _sym().as_dofdesc(dd_in)
-        official_dd_out = _sym().as_dofdesc(dd_out)
+        import grudge.symbolic.primitives as prim
+        official_dd_in = prim.as_dofdesc(dd_in)
+        official_dd_out = prim.as_dofdesc(dd_out)
+
         if official_dd_in == official_dd_out:
             raise ValueError("Interpolating from {} to {}"
             " does not do anything.".format(official_dd_in, official_dd_out))
@@ -123,7 +171,8 @@ class ElementwiseMaxOperator(Operator):
 class NodalReductionOperator(Operator):
     def __init__(self, dd_in, dd_out=None):
         if dd_out is None:
-            dd_out = _sym().DD_SCALAR
+            import grudge.symbolic.primitives as prim
+            dd_out = prim.DD_SCALAR
 
         assert dd_out.is_scalar()
 
@@ -150,13 +199,14 @@ class NodalMin(NodalReductionOperator):
 
 class DiffOperatorBase(Operator):
     def __init__(self, xyz_axis, dd_in=None, dd_out=None):
+        import grudge.symbolic.primitives as prim
         if dd_in is None:
-            dd_in = _sym().DD_VOLUME
+            dd_in = prim.DD_VOLUME
 
         if dd_out is None:
-            dd_out = dd_in.with_qtag(_sym().QTAG_NONE)
+            dd_out = dd_in.with_qtag(prim.QTAG_NONE)
         else:
-            dd_out = _sym().as_dofdesc(dd_out)
+            dd_out = prim.as_dofdesc(dd_out)
 
         if dd_out.uses_quadrature():
             raise ValueError("differentiation outputs are not on "
@@ -207,10 +257,13 @@ class MInvSTOperator(WeakFormDiffOperatorBase):
 
 class RefDiffOperatorBase(ElementwiseLinearOperator):
     def __init__(self, rst_axis, dd_in=None, dd_out=None):
+        import grudge.symbolic.primitives as prim
         if dd_in is None:
-            dd_in = _sym().DD_VOLUME
+            dd_in = prim.DD_VOLUME
+
         if dd_out is None:
-            dd_out = dd_in.with_qtag(_sym().QTAG_NONE)
+            dd_out = dd_in.with_qtag(prim.QTAG_NONE)
+
         if dd_out.uses_quadrature():
             raise ValueError("differentiation outputs are not on "
                     "quadrature grids")
@@ -275,8 +328,10 @@ class FilterOperator(ElementwiseLinearOperator):
           (For example an instance of
           :class:`ExponentialFilterResponseFunction`.
         """
+        import grudge.symbolic.primitives as prim
         if dd_in is None:
-            dd_in = _sym().DD_VOLUME
+            dd_in = prim.DD_VOLUME
+
         if dd_out is None:
             dd_out = dd_in
 
@@ -314,8 +369,10 @@ class FilterOperator(ElementwiseLinearOperator):
 
 class AveragingOperator(ElementwiseLinearOperator):
     def __init__(self, dd_in=None, dd_out=None):
+        import grudge.symbolic.primitives as prim
         if dd_in is None:
-            dd_in = _sym().DD_VOLUME
+            dd_in = prim.DD_VOLUME
+
         if dd_out is None:
             dd_out = dd_in
 
@@ -362,8 +419,9 @@ class MassOperatorBase(Operator):
     """
 
     def __init__(self, dd_in=None, dd_out=None):
+        import grudge.symbolic.primitives as prim
         if dd_in is None:
-            dd_in = _sym().DD_VOLUME
+            dd_in = prim.DD_VOLUME
         if dd_out is None:
             dd_out = dd_in
 
@@ -429,15 +487,14 @@ class OppositeInteriorFaceSwap(Operator):
     """
 
     def __init__(self, dd_in=None, dd_out=None, unique_id=None):
-        sym = _sym()
-
+        import grudge.symbolic.primitives as prim
         if dd_in is None:
-            dd_in = sym.DOFDesc(sym.FACE_RESTR_INTERIOR, None)
+            dd_in = prim.DOFDesc(prim.FACE_RESTR_INTERIOR, None)
         if dd_out is None:
             dd_out = dd_in
 
         super(OppositeInteriorFaceSwap, self).__init__(dd_in, dd_out)
-        if self.dd_in.domain_tag is not sym.FACE_RESTR_INTERIOR:
+        if self.dd_in.domain_tag is not prim.FACE_RESTR_INTERIOR:
             raise ValueError("dd_in must be an interior faces domain")
         if self.dd_out != self.dd_in:
             raise ValueError("dd_out and dd_in must be identical")
@@ -462,7 +519,7 @@ class OppositePartitionFaceSwap(Operator):
         MPI tag offset to keep different subexpressions apart in MPI traffic.
     """
     def __init__(self, dd_in=None, dd_out=None, unique_id=None):
-        sym = _sym()
+        import grudge.symbolic.primitives as prim
 
         if dd_in is None and dd_out is None:
             raise ValueError("dd_in or dd_out must be specified")
@@ -472,7 +529,7 @@ class OppositePartitionFaceSwap(Operator):
             dd_out = dd_in
 
         super(OppositePartitionFaceSwap, self).__init__(dd_in, dd_out)
-        if not isinstance(self.dd_in.domain_tag, sym.BTAG_PARTITION):
+        if not isinstance(self.dd_in.domain_tag, prim.BTAG_PARTITION):
             raise ValueError("dd_in must be a partition boundary faces domain")
         if self.dd_out != self.dd_in:
             raise ValueError("dd_out and dd_in must be identical")
@@ -492,20 +549,20 @@ class OppositePartitionFaceSwap(Operator):
 
 class FaceMassOperatorBase(ElementwiseLinearOperator):
     def __init__(self, dd_in=None, dd_out=None):
-        sym = _sym()
-
+        import grudge.symbolic.primitives as prim
         if dd_in is None:
-            dd_in = sym.DOFDesc(sym.FACE_RESTR_ALL, None)
+            dd_in = prim.DOFDesc(prim.FACE_RESTR_ALL, None)
 
         if dd_out is None or dd_out == "vol":
-            dd_out = sym.DOFDesc("vol", sym.QTAG_NONE)
+            dd_out = prim.DOFDesc("vol", prim.QTAG_NONE)
+
         if dd_out.uses_quadrature():
             raise ValueError("face mass operator outputs are not on "
                     "quadrature grids")
 
         if not dd_out.is_volume():
             raise ValueError("dd_out must be a volume domain")
-        if dd_in.domain_tag is not sym.FACE_RESTR_ALL:
+        if dd_in.domain_tag is not prim.FACE_RESTR_ALL:
             raise ValueError("dd_in must be an interior faces domain")
 
         super(FaceMassOperatorBase, self).__init__(dd_in, dd_out)
@@ -575,43 +632,40 @@ def stiffness_t(dim, dd_in=None, dd_out=None):
 
 
 def integral(arg, dd=None):
-    sym = _sym()
+    import grudge.symbolic.primitives as prim
 
     if dd is None:
-        dd = sym.DD_VOLUME
-
-    dd = sym.as_dofdesc(dd)
+        dd = prim.DD_VOLUME
+    dd = prim.as_dofdesc(dd)
 
-    return sym.NodalSum(dd)(
-            arg * sym.cse(
-                sym.MassOperator(dd_in=dd)(sym.Ones(dd)),
+    return NodalSum(dd)(
+            arg * prim.cse(
+                MassOperator(dd_in=dd)(prim.Ones(dd)),
                 "mass_quad_weights",
-                sym.cse_scope.DISCRETIZATION))
+                prim.cse_scope.DISCRETIZATION))
 
 
 def norm(p, arg, dd=None):
     """
     :arg arg: is assumed to be a vector, i.e. have shape ``(n,)``.
     """
-    sym = _sym()
+    import grudge.symbolic.primitives as prim
 
     if dd is None:
-        dd = sym.DD_VOLUME
-
-    dd = sym.as_dofdesc(dd)
+        dd = prim.DD_VOLUME
+    dd = prim.as_dofdesc(dd)
 
     if p == 2:
-        norm_squared = sym.NodalSum(dd_in=dd)(
-                sym.FunctionSymbol("fabs")(
-                    arg * sym.MassOperator()(arg)))
+        norm_squared = NodalSum(dd_in=dd)(
+                prim.fabs(arg * MassOperator()(arg)))
 
         if isinstance(norm_squared, np.ndarray):
             norm_squared = norm_squared.sum()
 
-        return sym.FunctionSymbol("sqrt")(norm_squared)
+        return prim.sqrt(norm_squared)
 
     elif p == np.Inf:
-        result = sym.NodalMax(dd_in=dd)(sym.FunctionSymbol("fabs")(arg))
+        result = NodalMax(dd_in=dd)(prim.fabs(arg))
         from pymbolic.primitives import Max
 
         if isinstance(result, np.ndarray):
diff --git a/grudge/symbolic/primitives.py b/grudge/symbolic/primitives.py
index 0b9fac20332e505b038ee2ed9bda94485048d8d2..ea9e89f6bb6a607ec96944b45c3f7f8c5f8a3193 100644
--- a/grudge/symbolic/primitives.py
+++ b/grudge/symbolic/primitives.py
@@ -27,45 +27,45 @@ THE SOFTWARE.
 from six.moves import range, intern
 
 import numpy as np
-import pymbolic.primitives
-from meshmode.mesh import BTAG_ALL, BTAG_REALLY_ALL, BTAG_NONE, BTAG_PARTITION  # noqa
-from meshmode.discretization.connection import (  # noqa
-        FACE_RESTR_ALL, FACE_RESTR_INTERIOR)
-
-from pymbolic.primitives import (  # noqa
+from pytools.obj_array import make_obj_array
+
+from meshmode.mesh import (
+        BTAG_ALL,
+        BTAG_REALLY_ALL,
+        BTAG_NONE,
+        BTAG_PARTITION)
+from meshmode.discretization.connection import (
+        FACE_RESTR_ALL,
+        FACE_RESTR_INTERIOR)
+
+import pymbolic.primitives as prim
+from pymbolic.primitives import (
+        Variable as VariableBase,
+        CommonSubexpression as CommonSubexpressionBase,
         cse_scope as cse_scope_base,
-        make_common_subexpression as cse, If, Comparison)
+        make_common_subexpression as cse)
 from pymbolic.geometric_algebra import MultiVector
-from pytools.obj_array import join_fields, make_obj_array  # noqa
 
 
-class ExpressionBase(pymbolic.primitives.Expression):
+class ExpressionBase(prim.Expression):
     def make_stringifier(self, originating_stringifier=None):
         from grudge.symbolic.mappers import StringifyMapper
         return StringifyMapper()
 
 
-def _sym():
-    # A hack to make referring to grudge.sym possible without
-    # circular imports and tons of typing.
-
-    from grudge import sym
-    return sym
-
-
 __doc__ = """
 
-.. currentmodule:: grudge.sym
-
-.. autoclass:: If
-
 DOF description
 ^^^^^^^^^^^^^^^
 
 .. autoclass:: DTAG_SCALAR
 .. autoclass:: DTAG_VOLUME_ALL
+.. autoclass:: DTAG_BOUNDARY
 .. autoclass:: QTAG_NONE
+
 .. autoclass:: DOFDesc
+.. autofunction:: as_dofdesc
+
 .. data:: DD_SCALAR
 .. data:: DD_VOLUME
 
@@ -74,10 +74,12 @@ Symbols
 
 .. autoclass:: Variable
 .. autoclass:: ScalarVariable
-.. autoclass:: make_sym_array
-.. autoclass:: make_sym_mv
 .. autoclass:: FunctionSymbol
 
+.. autofunction:: make_sym_array
+.. autofunction:: make_sym_mv
+
+.. function :: fabs(arg)
 .. function :: sqrt(arg)
 .. function :: exp(arg)
 .. function :: sin(arg)
@@ -90,11 +92,11 @@ Helpers
 
 .. autoclass:: OperatorBinding
 .. autoclass:: PrioritizedSubexpression
-.. autoclass:: BoundaryPair
 .. autoclass:: Ones
 
 Geometry data
 ^^^^^^^^^^^^^
+
 .. autoclass:: NodeCoordinateComponent
 .. autofunction:: nodes
 .. autofunction:: mv_nodes
@@ -111,20 +113,20 @@ Geometry data
 
 # {{{ DOF description
 
-class DTAG_SCALAR:  # noqa
+class DTAG_SCALAR:          # noqa: N801
     pass
 
 
-class DTAG_VOLUME_ALL:  # noqa
+class DTAG_VOLUME_ALL:      # noqa: N801
     pass
 
 
-class DTAG_BOUNDARY:  # noqa
+class DTAG_BOUNDARY:        # noqa: N801
     def __init__(self, tag):
         self.tag = tag
 
 
-class QTAG_NONE:  # noqa
+class QTAG_NONE:            # noqa: N801
     pass
 
 
@@ -133,14 +135,18 @@ class DOFDesc(object):
 
     .. attribute:: domain_tag
     .. attribute:: quadrature_tag
+
     .. automethod:: is_scalar
     .. automethod:: is_discretized
     .. automethod:: is_volume
     .. automethod:: is_boundary
     .. automethod:: is_trace
+
     .. automethod:: uses_quadrature
+
     .. automethod:: with_qtag
     .. automethod:: with_dtag
+
     .. automethod:: __eq__
     .. automethod:: __ne__
     .. automethod:: __hash__
@@ -149,46 +155,36 @@ class DOFDesc(object):
     def __init__(self, domain_tag, quadrature_tag=None):
         """
         :arg domain_tag: One of the following:
-            :class:`grudge.sym.DTAG_SCALAR` (or the string ``"scalar"``),
-            :class:`grudge.sym.DTAG_VOLUME_ALL` (or the string ``"vol"``)
+            :class:`DTAG_SCALAR` (or the string ``"scalar"``),
+            :class:`DTAG_VOLUME_ALL` (or the string ``"vol"``)
             for the default volume discretization,
-            :class:`meshmode.discretization.connection.FACE_RESTR_ALL`
-            (or the string ``"all_faces"``),
-            or
-            :class:`meshmode.discretization.connection.FACE_RESTR_INTERIOR`
-            (or the string ``"int_faces"``),
-            or one of
-            :class:`meshmode.discretization.BTAG_ALL`,
-            :class:`meshmode.discretization.BTAG_NONE`,
-            :class:`meshmode.discretization.BTAG_REALLY_ALL`,
-            :class:`meshmode.discretization.PARTITION`,
-            or :class
+            :data:`~meshmode.discretization.connection.FACE_RESTR_ALL`
+            (or the string ``"all_faces"``), or
+            :data:`~meshmode.discretization.connection.FACE_RESTR_INTERIOR`
+            (or the string ``"int_faces"``), or one of
+            :class:`~meshmode.mesh.BTAG_ALL`,
+            :class:`~meshmode.mesh.BTAG_NONE`,
+            :class:`~meshmode.mesh.BTAG_REALLY_ALL`,
+            :class:`~meshmode.mesh.BTAG_PARTITION`,
             or *None* to indicate that the geometry is not yet known.
 
         :arg quadrature_tag:
-            *None* to indicate that the quadrature grid is not known,or
+            *None* to indicate that the quadrature grid is not known, or
             :class:`QTAG_NONE` to indicate the use of the basic discretization
             grid, or a string to indicate the use of the thus-tagged quadratue
             grid.
         """
-        if domain_tag == "scalar":
-            domain_tag = DTAG_SCALAR
-        elif domain_tag is DTAG_SCALAR:
+
+        if domain_tag is None:
+            pass
+        elif domain_tag in [DTAG_SCALAR, "scalar"]:
             domain_tag = DTAG_SCALAR
-        elif domain_tag == "vol":
+        elif domain_tag in [DTAG_VOLUME_ALL, "vol"]:
             domain_tag = DTAG_VOLUME_ALL
-        elif domain_tag is DTAG_VOLUME_ALL:
-            pass
-        elif domain_tag == "all_faces":
+        elif domain_tag in [FACE_RESTR_ALL, "all_faces"]:
             domain_tag = FACE_RESTR_ALL
-        elif domain_tag is FACE_RESTR_ALL:
-            pass
-        elif domain_tag == "int_faces":
+        elif domain_tag in [FACE_RESTR_INTERIOR, "int_faces"]:
             domain_tag = FACE_RESTR_INTERIOR
-        elif domain_tag is FACE_RESTR_INTERIOR:
-            pass
-        elif domain_tag is None:
-            pass
         elif isinstance(domain_tag, BTAG_PARTITION):
             pass
         elif domain_tag in [BTAG_ALL, BTAG_REALLY_ALL, BTAG_NONE]:
@@ -273,8 +269,7 @@ DD_VOLUME = DOFDesc(DTAG_VOLUME_ALL, None)
 def as_dofdesc(dd):
     if isinstance(dd, DOFDesc):
         return dd
-    else:
-        return DOFDesc(dd, None)
+    return DOFDesc(dd, quadrature_tag=None)
 
 # }}}
 
@@ -314,11 +309,11 @@ class HasDOFDesc(object):
 
 # {{{ variables
 
-class cse_scope(cse_scope_base):  # noqa
+class cse_scope(cse_scope_base):        # noqa: N801
     DISCRETIZATION = "grudge_discretization"
 
 
-class Variable(HasDOFDesc, ExpressionBase, pymbolic.primitives.Variable):
+class Variable(HasDOFDesc, ExpressionBase, VariableBase):
     """A user-supplied input variable with a known :class:`DOFDesc`.
     """
     init_arg_names = ("name", "dd")
@@ -347,7 +342,7 @@ def make_sym_array(name, shape, dd=None):
     def var_factory(name):
         return Variable(name, dd)
 
-    return pymbolic.primitives.make_sym_array(name, shape, var_factory)
+    return prim.make_sym_array(name, shape, var_factory)
 
 
 def make_sym_mv(name, dim, var_factory=None):
@@ -355,15 +350,15 @@ def make_sym_mv(name, dim, var_factory=None):
             make_sym_array(name, dim, var_factory))
 
 
-class FunctionSymbol(ExpressionBase, pymbolic.primitives.Variable):
+class FunctionSymbol(ExpressionBase, VariableBase):
     """A symbol to be used as the function argument of
-    :class:`pymbolic.primitives.Call`.
-
+    :class:`~pymbolic.primitives.Call`.
     """
 
     mapper_method = "map_function_symbol"
 
 
+fabs = FunctionSymbol("fabs")
 sqrt = FunctionSymbol("sqrt")
 exp = FunctionSymbol("exp")
 sin = FunctionSymbol("sin")
@@ -399,7 +394,7 @@ class OperatorBinding(ExpressionBase):
         return hash((self.__class__, self.op, obj_array_to_hashable(self.field)))
 
 
-class PrioritizedSubexpression(pymbolic.primitives.CommonSubexpression):
+class PrioritizedSubexpression(CommonSubexpressionBase):
     """When the optemplate-to-code transformation is performed,
     prioritized subexpressions  work like common subexpression in
     that they are assigned their own separate identifier/register
@@ -409,7 +404,7 @@ class PrioritizedSubexpression(pymbolic.primitives.CommonSubexpression):
     """
 
     def __init__(self, child, priority=0):
-        pymbolic.primitives.CommonSubexpression.__init__(self, child)
+        super(PrioritizedSubexpression, self).__init__(child)
         self.priority = priority
 
     def __getinitargs__(self):
@@ -470,19 +465,21 @@ def forward_metric_derivative(xyz_axis, rst_axis, dd=None):
 
     .. math::
 
-        \frac{d x_{\mathtt{xyz\_axis}} }{d r_{\mathtt{rst\_axis}} }
+        \frac{\partial x_{\mathrm{xyz\_axis}} }{\partial r_{\mathrm{rst\_axis}} }
     """
     if dd is None:
         dd = DD_VOLUME
 
     inner_dd = dd.with_qtag(QTAG_NONE)
 
-    diff_op = _sym().RefDiffOperator(rst_axis, inner_dd)
+    from grudge.symbolic.operators import RefDiffOperator
+    diff_op = RefDiffOperator(rst_axis, inner_dd)
 
     result = diff_op(NodeCoordinateComponent(xyz_axis, inner_dd))
 
     if dd.uses_quadrature():
-        result = _sym().interp(inner_dd, dd)(result)
+        from grudge.symbolic.operators import interp
+        result = interp(inner_dd, dd)(result)
 
     return cse(
             result,
@@ -595,7 +592,7 @@ def area_element(ambient_dim, dim=None, dd=None):
 
 
 def mv_normal(dd, ambient_dim, dim=None):
-    """Exterior unit normal as a MultiVector."""
+    """Exterior unit normal as a :class:`~pymbolic.geometric_algebra.MultiVector`."""
 
     dd = as_dofdesc(dd)
 
@@ -676,13 +673,15 @@ class TracePair:
 
 
 def int_tpair(expression, qtag=None):
-    i = _sym().interp("vol", "int_faces")(expression)
-    e = cse(_sym().OppositeInteriorFaceSwap()(i))
+    from grudge.symbolic.operators import interp, OppositeInteriorFaceSwap
+
+    i = interp("vol", "int_faces")(expression)
+    e = cse(OppositeInteriorFaceSwap()(i))
 
-    if qtag is not None and qtag != _sym().QTAG_NONE:
-        q_dd = _sym().DOFDesc("int_faces", qtag)
-        i = cse(_sym().interp("int_faces", q_dd)(i))
-        e = cse(_sym().interp("int_faces", q_dd)(e))
+    if qtag is not None and qtag != QTAG_NONE:
+        q_dd = DOFDesc("int_faces", qtag)
+        i = cse(interp("int_faces", q_dd)(i))
+        e = cse(interp("int_faces", q_dd)(e))
     else:
         q_dd = "int_faces"
 
@@ -710,7 +709,8 @@ def bv_tpair(dd, interior, exterior):
         representing the exterior value to be used
         for the flux.
     """
-    interior = _sym().cse(_sym().interp("vol", dd)(interior))
+    from grudge.symbolic.operators import interp
+    interior = cse(interp("vol", dd)(interior))
     return TracePair(dd, interior, exterior)
 
 # }}}
diff --git a/test/test_grudge.py b/test/test_grudge.py
index e373e1335a86a9f384b19a54bb2a8842130decc8..1d2732a9acfb9b51c0aadf4534ddb85d49d46995 100644
--- a/test/test_grudge.py
+++ b/test/test_grudge.py
@@ -23,11 +23,14 @@ THE SOFTWARE.
 """
 
 
-import numpy as np  # noqa
-import numpy.linalg as la  # noqa
-import pyopencl as cl  # noqa
-import pyopencl.array  # noqa
-import pyopencl.clmath  # noqa
+import numpy as np
+import numpy.linalg as la
+
+import pyopencl as cl
+import pyopencl.array
+import pyopencl.clmath
+
+from pytools.obj_array import join_fields
 
 import pytest  # noqa
 
@@ -185,7 +188,7 @@ def test_2d_gauss_theorem(ctx_factory):
     discr = DGDiscretizationWithBoundaries(cl_ctx, mesh, order=2)
 
     def f(x):
-        return sym.join_fields(
+        return join_fields(
                 sym.sin(3*x[0])+sym.cos(3*x[1]),
                 sym.sin(2*x[0])+sym.cos(x[1]))
 
@@ -408,7 +411,6 @@ def test_improvement_quadrature(ctx_factory, order):
     from meshmode.mesh.generation import generate_regular_rect_mesh
     from grudge.models.advection import VariableCoefficientAdvectionOperator
     from pytools.convergence import EOCRecorder
-    from pytools.obj_array import join_fields
     from meshmode.discretization.poly_element import QuadratureSimplexGroupFactory
 
     cl_ctx = cl.create_some_context()