diff --git a/doc/conf.py b/doc/conf.py
index dc8f0e244b6a4212b817e4751f54e348bbe9661b..1b1cb89eb3abdf59de9093124c4a8bdc5ca4e67e 100644
--- a/doc/conf.py
+++ b/doc/conf.py
@@ -1,6 +1,3 @@
-from __future__ import absolute_import
-# -*- coding: utf-8 -*-
-#
 # grudge documentation build configuration file, created by
 # sphinx-quickstart on Sun Sep 27 13:08:30 2015.
 #
@@ -20,44 +17,44 @@ import shlex  # noqa
 # If extensions (or modules to document with autodoc) are in another directory,
 # add these directories to sys.path here. If the directory is relative to the
 # documentation root, use os.path.abspath to make it absolute, like shown here.
-#sys.path.insert(0, os.path.abspath('.'))
+#sys.path.insert(0, os.path.abspath("."))
 
 # -- General configuration ------------------------------------------------
 
 # If your documentation needs a minimal Sphinx version, state it here.
-#needs_sphinx = '1.0'
+#needs_sphinx = "1.0"
 
 # Add any Sphinx extension module names here, as strings. They can be
-# extensions coming with Sphinx (named 'sphinx.ext.*') or your custom
+# extensions coming with Sphinx (named "sphinx.ext.*") or your custom
 # ones.
 extensions = [
-    'sphinx.ext.autodoc',
-    'sphinx.ext.doctest',
-    'sphinx.ext.intersphinx',
-    'sphinx.ext.mathjax',
-    'sphinx.ext.viewcode',
+    "sphinx.ext.autodoc",
+    "sphinx.ext.doctest",
+    "sphinx.ext.intersphinx",
+    "sphinx.ext.mathjax",
+    "sphinx.ext.viewcode",
 ]
 
 # Add any paths that contain templates here, relative to this directory.
-templates_path = ['_templates']
+templates_path = ["_templates"]
 
 # The suffix(es) of source filenames.
 # You can specify multiple suffix as a list of string:
-# source_suffix = ['.rst', '.md']
-source_suffix = '.rst'
+# source_suffix = [".rst", ".md"]
+source_suffix = ".rst"
 
 # The encoding of source files.
-#source_encoding = 'utf-8-sig'
+#source_encoding = "utf-8-sig"
 
 # The master toctree document.
-master_doc = 'index'
+master_doc = "index"
 
 # General information about the project.
-project = u'grudge'
-copyright = u'2015, Andreas Kloeckner'
-author = u'Andreas Kloeckner'
+project = "grudge"
+copyright = "2015, Andreas Kloeckner"
+author = "Andreas Kloeckner"
 
-# The version info for the project you're documenting, acts as replacement for
+# The version info for the project you"re documenting, acts as replacement for
 # |version| and |release|, also used in various other places throughout the
 # built documents.
 #
@@ -68,10 +65,11 @@ def get_version():
     conf = {}
     src = "../grudge/version.py"
     exec(
-            compile(open(src).read(), src, 'exec'),
+            compile(open(src).read(), src, "exec"),
             conf)
     return conf["VERSION_TEXT"]
 
+
 version = get_version()
 
 # The full version, including alpha/beta/rc tags.
@@ -86,19 +84,19 @@ language = None
 
 # There are two options for replacing |today|: either, you set today to some
 # non-false value, then it is used:
-#today = ''
+#today = ""
 # Else, today_fmt is used as the format for a strftime call.
-#today_fmt = '%B %d, %Y'
+#today_fmt = "%B %d, %Y"
 
 # List of patterns, relative to source directory, that match files and
 # directories to ignore when looking for source files.
-exclude_patterns = ['_build']
+exclude_patterns = ["_build"]
 
 # The reST default role (used for this markup: `text`) to use for all
 # documents.
 #default_role = None
 
-# If true, '()' will be appended to :func: etc. cross-reference text.
+# If true, "()" will be appended to :func: etc. cross-reference text.
 #add_function_parentheses = True
 
 # If true, the current module name will be prepended to all description
@@ -110,7 +108,7 @@ exclude_patterns = ['_build']
 #show_authors = False
 
 # The name of the Pygments (syntax highlighting) style to use.
-pygments_style = 'sphinx'
+pygments_style = "sphinx"
 
 # A list of ignored prefixes for module index sorting.
 #modindex_common_prefix = []
@@ -134,11 +132,11 @@ html_theme_options = {
         }
 
 html_sidebars = {
-    '**': [
-        'about.html',
-        'navigation.html',
-        'relations.html',
-        'searchbox.html',
+    "**": [
+        "about.html",
+        "navigation.html",
+        "relations.html",
+        "searchbox.html",
     ]
 }
 
@@ -170,16 +168,16 @@ html_sidebars = {
 # Add any paths that contain custom static files (such as style sheets) here,
 # relative to this directory. They are copied after the builtin static files,
 # so a file named "default.css" will overwrite the builtin "default.css".
-#html_static_path = ['_static']
+#html_static_path = ["_static"]
 
 # Add any extra paths that contain custom files (such as robots.txt or
 # .htaccess) here, relative to this directory. These files are copied
 # directly to the root of the documentation.
 #html_extra_path = []
 
-# If not '', a 'Last updated on:' timestamp is inserted at every page bottom,
+# If not "", a "Last updated on:" timestamp is inserted at every page bottom,
 # using the given strftime format.
-#html_last_updated_fmt = '%b %d, %Y'
+#html_last_updated_fmt = "%b %d, %Y"
 
 # If true, SmartyPants will be used to convert quotes and dashes to
 # typographically correct entities.
@@ -213,50 +211,50 @@ html_sidebars = {
 # If true, an OpenSearch description file will be output, and all pages will
 # contain a <link> tag referring to it.  The value of this option must be the
 # base URL from which the finished HTML is served.
-#html_use_opensearch = ''
+#html_use_opensearch = ""
 
 # This is the file name suffix for HTML files (e.g. ".xhtml").
 #html_file_suffix = None
 
 # Language to be used for generating the HTML full-text search index.
 # Sphinx supports the following languages:
-#   'da', 'de', 'en', 'es', 'fi', 'fr', 'hu', 'it', 'ja'
-#   'nl', 'no', 'pt', 'ro', 'ru', 'sv', 'tr'
-#html_search_language = 'en'
+#   "da", "de", "en", "es", "fi", "fr", "hu", "it", "ja"
+#   "nl", "no", "pt", "ro", "ru", "sv", "tr"
+#html_search_language = "en"
 
 # A dictionary with options for the search language support, empty by default.
-# Now only 'ja' uses this config value
-#html_search_options = {'type': 'default'}
+# Now only "ja" uses this config value
+#html_search_options = {"type": "default"}
 
 # The name of a javascript file (relative to the configuration directory) that
 # implements a search results scorer. If empty, the default will be used.
-#html_search_scorer = 'scorer.js'
+#html_search_scorer = "scorer.js"
 
 # Output file base name for HTML help builder.
-htmlhelp_basename = 'grudgedoc'
+htmlhelp_basename = "grudgedoc"
 
 # -- Options for LaTeX output ---------------------------------------------
 
 latex_elements = {
-    # The paper size ('letterpaper' or 'a4paper').
-    #'papersize': 'letterpaper',
+    # The paper size ("letterpaper" or "a4paper").
+    #"papersize": "letterpaper",
 
-    # The font size ('10pt', '11pt' or '12pt').
-    #'pointsize': '10pt',
+    # The font size ("10pt", "11pt" or "12pt").
+    #"pointsize": "10pt",
 
     # Additional stuff for the LaTeX preamble.
-    #'preamble': '',
+    #"preamble": "",
 
     # Latex figure (float) alignment
-    #'figure_align': 'htbp',
+    #"figure_align": "htbp",
     }
 
 # Grouping the document tree into LaTeX files. List of tuples
 # (source start file, target name, title,
 #  author, documentclass [howto, manual, or own class]).
 latex_documents = [
-    (master_doc, 'grudge.tex', u'grudge Documentation',
-     u'Andreas Kloeckner', 'manual'),
+    (master_doc, "grudge.tex", "grudge Documentation",
+     "Andreas Kloeckner", "manual"),
 ]
 
 # The name of an image file (relative to this directory) to place at the top of
@@ -285,7 +283,7 @@ latex_documents = [
 # One entry per manual page. List of tuples
 # (source start file, name, description, authors, manual section).
 man_pages = [
-    (master_doc, 'grudge', u'grudge Documentation',
+    (master_doc, "grudge", "grudge Documentation",
      [author], 1)
 ]
 
@@ -299,9 +297,9 @@ man_pages = [
 # (source start file, target name, title, author,
 #  dir menu entry, description, category)
 texinfo_documents = [
-    (master_doc, 'grudge', u'grudge Documentation',
-     author, 'grudge', 'One line description of project.',
-     'Miscellaneous'),
+    (master_doc, "grudge", "grudge Documentation",
+     author, "grudge", "One line description of project.",
+     "Miscellaneous"),
 ]
 
 # Documents to append as an appendix to all manuals.
@@ -310,21 +308,21 @@ texinfo_documents = [
 # If false, no module index is generated.
 #texinfo_domain_indices = True
 
-# How to display URL addresses: 'footnote', 'no', or 'inline'.
-#texinfo_show_urls = 'footnote'
+# How to display URL addresses: "footnote", "no", or "inline".
+#texinfo_show_urls = "footnote"
 
-# If true, do not generate a @detailmenu in the "Top" node's menu.
+# If true, do not generate a @detailmenu in the "Top" node"s menu.
 #texinfo_no_detailmenu = False
 
 
 # Example configuration for intersphinx: refer to the Python standard library.
 intersphinx_mapping = {
-    'https://docs.python.org/': None,
-    'https://docs.scipy.org/doc/numpy/': None,
-    'https://documen.tician.de/pyopencl/': None,
-    'https://documen.tician.de/modepy/': None,
-    'https://documen.tician.de/pymbolic/': None,
-    'https://documen.tician.de/meshmode/': None,
-    'http://documen.tician.de/loopy/': None,
+    "https://docs.python.org/": None,
+    "https://docs.scipy.org/doc/numpy/": None,
+    "https://documen.tician.de/pyopencl/": None,
+    "https://documen.tician.de/modepy/": None,
+    "https://documen.tician.de/pymbolic/": None,
+    "https://documen.tician.de/meshmode/": None,
+    "https://documen.tician.de/loopy/": None,
     }
 autoclass_content = "class"
diff --git a/examples/advection/surface.py b/examples/advection/surface.py
index c0bb772cd4c9b3c811530b2835466c079ef97184..2ab36d516416796c6760023000638ead0004fbda 100644
--- a/examples/advection/surface.py
+++ b/examples/advection/surface.py
@@ -73,11 +73,11 @@ class Plotter:
 
             ax = self.fig.gca()
             ax.grid()
-            ax.plot(self.x, u, '-')
-            ax.plot(self.x, u, 'k.')
+            ax.plot(self.x, u, "-")
+            ax.plot(self.x, u, "k.")
             ax.set_xlabel(r"$\theta$")
             ax.set_ylabel("$u$")
-            ax.set_title("t = {:.2f}".format(evt.t))
+            ax.set_title(f"t = {evt.t:.2f}")
 
             self.fig.savefig(filename)
             self.fig.clf()
diff --git a/examples/advection/var-velocity.py b/examples/advection/var-velocity.py
index cc2e62f06125d7b4f16b4c291911d9184f6fafd0..9aa4d5b5faf5535fd0b0b1dbd4c958a690fbfec0 100644
--- a/examples/advection/var-velocity.py
+++ b/examples/advection/var-velocity.py
@@ -1,5 +1,3 @@
-from __future__ import division, absolute_import
-
 __copyright__ = "Copyright (C) 2017 Bogdan Enache"
 
 __license__ = """
@@ -72,14 +70,14 @@ class Plotter:
                 raise FileExistsError("output file '%s' already exists" % filename)
 
             ax = self.fig.gca()
-            ax.plot(self.x, u, '-')
-            ax.plot(self.x, u, 'k.')
+            ax.plot(self.x, u, "-")
+            ax.plot(self.x, u, "k.")
             if self.ylim is not None:
                 ax.set_ylim(self.ylim)
 
             ax.set_xlabel("$x$")
             ax.set_ylabel("$u$")
-            ax.set_title("t = {:.2f}".format(evt.t))
+            ax.set_title(f"t = {evt.t:.2f}")
             self.fig.savefig(filename)
             self.fig.clf()
         else:
diff --git a/examples/advection/weak.py b/examples/advection/weak.py
index 77ce17bc54ed03e84da4f5761d0c251abeb52824..6e4edd0170175f4fb705d235317c6d80f997c752 100644
--- a/examples/advection/weak.py
+++ b/examples/advection/weak.py
@@ -1,5 +1,3 @@
-from __future__ import division, absolute_import
-
 __copyright__ = "Copyright (C) 2007 Andreas Kloeckner"
 
 __license__ = """
@@ -72,14 +70,14 @@ class Plotter:
                 raise FileExistsError("output file '%s' already exists" % filename)
 
             ax = self.fig.gca()
-            ax.plot(self.x, u, '-')
-            ax.plot(self.x, u, 'k.')
+            ax.plot(self.x, u, "-")
+            ax.plot(self.x, u, "k.")
             if self.ylim is not None:
                 ax.set_ylim(self.ylim)
 
             ax.set_xlabel("$x$")
             ax.set_ylabel("$u$")
-            ax.set_title("t = {:.2f}".format(evt.t))
+            ax.set_title(f"t = {evt.t:.2f}")
 
             self.fig.savefig(filename)
             self.fig.clf()
diff --git a/examples/dagrt-fusion.py b/examples/dagrt-fusion.py
index bab82d545df9a7ba9bf4f6b41937dc10f45191df..2a285ad533b0d1ffb545da65f69efc8665072a10 100755
--- a/examples/dagrt-fusion.py
+++ b/examples/dagrt-fusion.py
@@ -2,8 +2,6 @@
 """Study of operator fusion (inlining) for time integration operators in Grudge.
 """
 
-from __future__ import division, print_function
-
 __copyright__ = """
 Copyright (C) 2015 Andreas Kloeckner
 Copyright (C) 2019 Matt Wala
@@ -270,7 +268,7 @@ def transcribe_phase(dag, field_var_name, field_components, phase_name,
 
 # {{{ time integrator implementations
 
-class RK4TimeStepperBase(object):
+class RK4TimeStepperBase:
 
     def __init__(self, component_getter):
         self.component_getter = component_getter
@@ -829,7 +827,7 @@ def test_stepper_mem_ops(ctx_factory, use_fusion):
 SECONDS_PER_NANOSECOND = 10**9
 
 
-class TimingFuture(object):
+class TimingFuture:
 
     def __init__(self, start_event, stop_event):
         self.start_event = start_event
diff --git a/examples/geometry.py b/examples/geometry.py
index df81298d4865b0dce384952633ee4fd00717c72b..0affdc8ecf5ba887a409758152c58420609f74db 100644
--- a/examples/geometry.py
+++ b/examples/geometry.py
@@ -1,7 +1,5 @@
 """Minimal example of a grudge driver."""
 
-from __future__ import division, print_function
-
 __copyright__ = "Copyright (C) 2015 Andreas Kloeckner"
 
 __license__ = """
@@ -24,6 +22,7 @@ OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
 THE SOFTWARE.
 """
 
+
 import numpy as np  # noqa
 import pyopencl as cl
 from grudge import sym, bind, DGDiscretizationWithBoundaries, shortcuts
diff --git a/examples/maxwell/cavities.py b/examples/maxwell/cavities.py
index eb2ee28dae6b2ed74c2760d2359d9542a9df2f59..d31e08219f2fb944e980d93d5f63b9aee754d6d2 100644
--- a/examples/maxwell/cavities.py
+++ b/examples/maxwell/cavities.py
@@ -1,7 +1,5 @@
 """Minimal example of a grudge driver."""
 
-from __future__ import division, print_function
-
 __copyright__ = "Copyright (C) 2015 Andreas Kloeckner"
 
 __license__ = """
diff --git a/examples/wave/var-propagation-speed.py b/examples/wave/var-propagation-speed.py
index 858e72549413b8120962a8052ebd7174622fff18..6a075650184f94e7ee820806c7706fb337e7278c 100644
--- a/examples/wave/var-propagation-speed.py
+++ b/examples/wave/var-propagation-speed.py
@@ -1,7 +1,5 @@
 """Minimal example of a grudge driver."""
 
-from __future__ import division, print_function
-
 __copyright__ = "Copyright (C) 2015 Andreas Kloeckner"
 
 __license__ = """
diff --git a/examples/wave/wave-eager-mpi.py b/examples/wave/wave-eager-mpi.py
index a6ca4381c9c84345d5d7eec64e470c43e05a20e1..5971b2ab9317935316ae33406c09e07a903ac06b 100644
--- a/examples/wave/wave-eager-mpi.py
+++ b/examples/wave/wave-eager-mpi.py
@@ -1,5 +1,3 @@
-from __future__ import division, print_function
-
 __copyright__ = "Copyright (C) 2020 Andreas Kloeckner"
 
 __license__ = """
diff --git a/examples/wave/wave-eager-var-velocity.py b/examples/wave/wave-eager-var-velocity.py
index 5164562cd03453f4a544fd2feb724bc44f8ce1d6..1ea155cb7422a49b22d62c486c628f417d23cf42 100644
--- a/examples/wave/wave-eager-var-velocity.py
+++ b/examples/wave/wave-eager-var-velocity.py
@@ -1,5 +1,3 @@
-from __future__ import division, print_function
-
 __copyright__ = "Copyright (C) 2020 Andreas Kloeckner"
 
 __license__ = """
diff --git a/examples/wave/wave-eager.py b/examples/wave/wave-eager.py
index 21c6d73a79d32adaae116e2183016ad88e0a82b5..4ac1a646810aab460c4969005dcc4db302624b8e 100644
--- a/examples/wave/wave-eager.py
+++ b/examples/wave/wave-eager.py
@@ -1,5 +1,3 @@
-from __future__ import division, print_function
-
 __copyright__ = "Copyright (C) 2020 Andreas Kloeckner"
 
 __license__ = """
diff --git a/examples/wave/wave-min-mpi.py b/examples/wave/wave-min-mpi.py
index 3e0063d866942bcc921f583c9464ee6d78f33bd7..a65582b7dd1227868a3ae81f67ca08b1c37b4237 100644
--- a/examples/wave/wave-min-mpi.py
+++ b/examples/wave/wave-min-mpi.py
@@ -1,7 +1,5 @@
 """Minimal example of a grudge driver."""
 
-from __future__ import division, print_function
-
 __copyright__ = "Copyright (C) 2015 Andreas Kloeckner"
 
 __license__ = """
diff --git a/examples/wave/wave-min.py b/examples/wave/wave-min.py
index e567251aa9631b1fe1da7e648208d91802c42ad7..c5d0e6017274f60996e9db838d0b4707e00ca4e2 100644
--- a/examples/wave/wave-min.py
+++ b/examples/wave/wave-min.py
@@ -1,6 +1,5 @@
 """Minimal example of a grudge driver."""
 
-from __future__ import division, print_function
 
 __copyright__ = "Copyright (C) 2015 Andreas Kloeckner"
 
diff --git a/grudge/__init__.py b/grudge/__init__.py
index 17358dd474304bec731bedb58a98cd692ec69935..cbc7de30a8847e908ff013038d69d3f022e5061d 100644
--- a/grudge/__init__.py
+++ b/grudge/__init__.py
@@ -1,5 +1,3 @@
-from __future__ import division, absolute_import
-
 __copyright__ = "Copyright (C) 2015 Andreas Kloeckner"
 
 __license__ = """
diff --git a/grudge/discretization.py b/grudge/discretization.py
index a5452a24b9e690ec14cd0b5bf455f711f3b7f2e7..ce130aef8dbad0cb165653774960ae96b722ea8c 100644
--- a/grudge/discretization.py
+++ b/grudge/discretization.py
@@ -1,5 +1,3 @@
-from __future__ import division, absolute_import
-
 __copyright__ = "Copyright (C) 2015-2017 Andreas Kloeckner, Bogdan Enache"
 
 __license__ = """
@@ -22,8 +20,6 @@ OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
 THE SOFTWARE.
 """
 
-
-import six
 from pytools import memoize_method
 from grudge import sym
 import numpy as np  # noqa: F401
@@ -35,7 +31,7 @@ __doc__ = """
 """
 
 
-class DGDiscretizationWithBoundaries(object):
+class DGDiscretizationWithBoundaries:
     """
     .. automethod :: discr_from_dd
     .. automethod :: connection_from_dds
@@ -141,7 +137,7 @@ class DGDiscretizationWithBoundaries(object):
                     i_remote_part, grp_factory)
             setup_helper.post_sends()
 
-        for i_remote_part, setup_helper in six.iteritems(setup_helpers):
+        for i_remote_part, setup_helper in setup_helpers.items():
             boundary_connections[i_remote_part] = setup_helper.complete_setup()
 
         return boundary_connections
diff --git a/grudge/dt_finding.py b/grudge/dt_finding.py
index 86e8f9a4be9210e825c25388e2c62cb3925de58d..fd5a3d8e2ddcc2d13170182656eb59d40a23b798 100644
--- a/grudge/dt_finding.py
+++ b/grudge/dt_finding.py
@@ -1,7 +1,5 @@
 """Helpers for estimating a stable time step."""
 
-from __future__ import division, print_function
-
 __copyright__ = "Copyright (C) 2015 Andreas Kloeckner"
 
 __license__ = """
@@ -29,7 +27,7 @@ from meshmode.discretization.poly_element import PolynomialWarpAndBlendElementGr
 import numpy.linalg as la
 
 
-class WarpAndBlendTimestepInfo(object):
+class WarpAndBlendTimestepInfo:
     @staticmethod
     def dt_non_geometric_factor(discr, grp):
         if grp.dim == 1:
diff --git a/grudge/eager.py b/grudge/eager.py
index 38efbc02739560344b1c8084761dce5cd9bbc6f1..efc60549e339f2b80ef634e73cb5613465efd8a8 100644
--- a/grudge/eager.py
+++ b/grudge/eager.py
@@ -1,5 +1,3 @@
-from __future__ import division, print_function
-
 __copyright__ = "Copyright (C) 2020 Andreas Kloeckner"
 
 __license__ = """
diff --git a/grudge/execution.py b/grudge/execution.py
index 59276ec05aabfcdc0eaee55f93749be07318086f..54f93a63602f8c34f47f238d4c3a14d2f2063767 100644
--- a/grudge/execution.py
+++ b/grudge/execution.py
@@ -1,5 +1,3 @@
-from __future__ import division, absolute_import
-
 __copyright__ = "Copyright (C) 2015-2017 Andreas Kloeckner, Bogdan Enache"
 
 __license__ = """
@@ -58,7 +56,7 @@ class ExecutionMapper(mappers.Evaluator,
         mappers.BoundOpMapperMixin,
         mappers.LocalOpReducerMixin):
     def __init__(self, array_context, context, bound_op):
-        super(ExecutionMapper, self).__init__(context)
+        super().__init__(context)
         self.discrwb = bound_op.discrwb
         self.bound_op = bound_op
         self.function_registry = bound_op.function_registry
@@ -95,7 +93,7 @@ class ExecutionMapper(mappers.Evaluator,
         return value
 
     def map_subscript(self, expr):
-        value = super(ExecutionMapper, self).map_subscript(expr)
+        value = super().map_subscript(expr)
 
         if isinstance(expr.aggregate, sym.Variable):
             dd = expr.aggregate.dd
@@ -556,7 +554,7 @@ class ExecutionMapper(mappers.Evaluator,
 
 # {{{ futures
 
-class MPIRecvFuture(object):
+class MPIRecvFuture:
     def __init__(self, array_context, bdry_discr, recv_req, insn_name,
             remote_data_host):
         self.array_context = array_context
@@ -576,7 +574,7 @@ class MPIRecvFuture(object):
         return [(self.insn_name, remote_data)], []
 
 
-class MPISendFuture(object):
+class MPISendFuture:
     def __init__(self, send_request):
         self.send_request = send_request
 
@@ -592,7 +590,7 @@ class MPISendFuture(object):
 
 # {{{ bound operator
 
-class BoundOperator(object):
+class BoundOperator:
     def __init__(self, discrwb, discr_code, eval_code, debug_flags,
             function_registry, exec_mapper_factory):
         self.discrwb = discrwb
diff --git a/grudge/function_registry.py b/grudge/function_registry.py
index 102932261bffbe3c52218611f872f2c912c41e65..89629deffcf450ab41c6869d37e6b5e7aa666a0e 100644
--- a/grudge/function_registry.py
+++ b/grudge/function_registry.py
@@ -1,5 +1,3 @@
-from __future__ import division, with_statement
-
 __copyright__ = """
 Copyright (C) 2013 Andreas Kloeckner
 Copyright (C) 2019 Matt Wala
@@ -65,7 +63,7 @@ class Function(RecordWithoutPickling):
     """
 
     def __init__(self, identifier, **kwargs):
-        super(Function, self).__init__(identifier=identifier, **kwargs)
+        super().__init__(identifier=identifier, **kwargs)
 
     def __call__(self, queue, *args, **kwargs):
         """Call the function implementation, if available."""
@@ -120,7 +118,7 @@ class FixedDOFDescExternalFunction(Function):
     supports_codegen = False
 
     def __init__(self, identifier, implementation, dd):
-        super(FixedDOFDescExternalFunction, self).__init__(
+        super().__init__(
                 identifier,
                 implementation=implementation,
                 dd=dd)
@@ -141,8 +139,7 @@ class FunctionRegistry(RecordWithoutPickling):
         if id_to_function is None:
             id_to_function = {}
 
-        super(FunctionRegistry, self).__init__(
-                id_to_function=id_to_function)
+        super().__init__(id_to_function=id_to_function)
 
     def register(self, function):
         """Return a copy of *self* with *function* registered."""
diff --git a/grudge/models/__init__.py b/grudge/models/__init__.py
index 9ccd1ab95a4ace69555c1cb1d030023cdaf85b04..6b2d4e59e8c503d123a1790112afc9ce7c1eeef1 100644
--- a/grudge/models/__init__.py
+++ b/grudge/models/__init__.py
@@ -1,8 +1,5 @@
 """Base classes for operators."""
 
-from __future__ import division
-from __future__ import absolute_import
-
 __copyright__ = "Copyright (C) 2007 Andreas Kloeckner"
 
 __license__ = """
@@ -26,7 +23,7 @@ THE SOFTWARE.
 """
 
 
-class Operator(object):
+class Operator:
     """A base class for Discontinuous Galerkin operators.
 
     You may derive your own operators from this class, but, at present
@@ -52,7 +49,7 @@ class HyperbolicOperator(Operator):
     """A base class for hyperbolic Discontinuous Galerkin operators."""
 
     def estimate_rk4_timestep(self, discr, t=None, fields=None):
-        u"""Estimate the largest stable timestep for an RK4 method.
+        """Estimate the largest stable timestep for an RK4 method.
         """
 
         from grudge.dt_finding import (
diff --git a/grudge/models/advection.py b/grudge/models/advection.py
index b752a0dccea4a9dc0c6036a6432786f16dd32ec1..55840b5a110a80825c87e6bd13963942c0ac1a13 100644
--- a/grudge/models/advection.py
+++ b/grudge/models/advection.py
@@ -1,8 +1,5 @@
-# -*- coding: utf8 -*-
 """Operators modeling advective phenomena."""
 
-from __future__ import division, absolute_import
-
 __copyright__ = "Copyright (C) 2009-2017 Andreas Kloeckner, Bogdan Enache"
 
 __license__ = """
@@ -53,7 +50,7 @@ def advection_weak_flux(flux_type, u, velocity):
                 )
         return u_upwind * v_dot_n
     else:
-        raise ValueError("flux `{}` is not implemented".format(flux_type))
+        raise ValueError(f"flux '{flux_type}' is not implemented")
 
 # }}}
 
@@ -74,7 +71,7 @@ class AdvectionOperatorBase(HyperbolicOperator):
         self.flux_type = flux_type
 
         if flux_type not in self.flux_types:
-            raise ValueError("unknown flux type: {}".format(flux_type))
+            raise ValueError(f"unknown flux type: '{flux_type}'")
 
     def weak_flux(self, u):
         return advection_weak_flux(self.flux_type, u, self.v)
@@ -143,7 +140,7 @@ class WeakAdvectionOperator(AdvectionOperatorBase):
 
 class VariableCoefficientAdvectionOperator(AdvectionOperatorBase):
     def __init__(self, v, inflow_u, flux_type="central", quad_tag="product"):
-        super(VariableCoefficientAdvectionOperator, self).__init__(
+        super().__init__(
                 v, inflow_u, flux_type=flux_type)
 
         self.quad_tag = quad_tag
@@ -209,12 +206,12 @@ def surface_advection_weak_flux(flux_type, u, velocity):
     elif flux_type == "lf":
         return u.avg * v_dot_n + 0.5 * sym.fabs(v_dot_n) * (u.int - u.ext)
     else:
-        raise ValueError("flux `{}` is not implemented".format(flux_type))
+        raise ValueError(f"flux '{flux_type}' is not implemented")
 
 
 class SurfaceAdvectionOperator(AdvectionOperatorBase):
     def __init__(self, v, flux_type="central", quad_tag=None):
-        super(SurfaceAdvectionOperator, self).__init__(
+        super().__init__(
                 v, inflow_u=None, flux_type=flux_type)
         self.quad_tag = quad_tag
 
diff --git a/grudge/models/burgers.py b/grudge/models/burgers.py
index ade3858c38a58339a813555bf1354480a7c91a6e..0be0b6ffffe5027fd360d21fa999cceb99ca64ed 100644
--- a/grudge/models/burgers.py
+++ b/grudge/models/burgers.py
@@ -1,9 +1,5 @@
-# -*- coding: utf8 -*-
 """Burgers operator."""
 
-from __future__ import division
-from __future__ import absolute_import
-
 __copyright__ = "Copyright (C) 2009 Andreas Kloeckner"
 
 __license__ = """
@@ -26,18 +22,13 @@ OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
 THE SOFTWARE.
 """
 
-
-
-
 from grudge.models import HyperbolicOperator
 import numpy
 from grudge.second_order import CentralSecondDerivative
 
 
-
-
 class BurgersOperator(HyperbolicOperator):
-    def __init__(self, dimensions, viscosity=None, 
+    def __init__(self, dimensions, viscosity=None,
             viscosity_scheme=CentralSecondDerivative()):
         # yes, you read that right--no BCs, 1D only.
         # (well--you can run the operator on a 2D grid. If you must.)
@@ -95,8 +86,8 @@ class BurgersOperator(HyperbolicOperator):
         num_flux = make_lax_friedrichs_flux(
                 #u0,
                 to_if_quad(emax_u),
-                make_obj_array([to_if_quad(u)]), 
-                [make_obj_array([flux(to_if_quad(u))])], 
+                make_obj_array([to_if_quad(u)]),
+                [make_obj_array([flux(to_if_quad(u))])],
                 [], strong=False)[0]
 
         from grudge.second_order import SecondDerivativeTarget
diff --git a/grudge/models/diffusion.py b/grudge/models/diffusion.py
index 04c5cf57c5368db4b88dfd4008b8962490f4d317..4c2db682f92ce985121e140568732e166ca4c9e3 100644
--- a/grudge/models/diffusion.py
+++ b/grudge/models/diffusion.py
@@ -1,9 +1,5 @@
-# -*- coding: utf8 -*-
 """Operators modeling diffusive phenomena."""
 
-from __future__ import division
-from __future__ import absolute_import
-
 __copyright__ = "Copyright (C) 2009 Andreas Kloeckner"
 
 __license__ = """
@@ -26,8 +22,6 @@ OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
 THE SOFTWARE.
 """
 
-
-
 import numpy
 
 import grudge.data
@@ -36,8 +30,6 @@ from grudge.models.poisson import LaplacianOperatorBase
 from grudge.second_order import CentralSecondDerivative
 
 
-
-
 class DiffusionOperator(TimeDependentOperator, LaplacianOperatorBase):
     def __init__(self, dimensions, diffusion_tensor=None,
             dirichlet_bc=grudge.data.make_tdep_constant(0), dirichlet_tag="dirichlet",
@@ -66,10 +58,10 @@ class DiffusionOperator(TimeDependentOperator, LaplacianOperatorBase):
 
         return BoundDiffusionOperator(self, discr)
 
-    def estimate_timestep(self, discr, 
+    def estimate_timestep(self, discr,
             stepper=None, stepper_class=None, stepper_args=None,
             t=None, fields=None):
-        u"""Estimate the largest stable timestep, given a time stepper
+        """Estimate the largest stable timestep, given a time stepper
         `stepper_class`. If none is given, RK4 is assumed.
         """
 
diff --git a/grudge/models/em.py b/grudge/models/em.py
index e21ad4994544f93a6cd135988f963ddec17c1d47..d8fbdb4ee0d967b8531530c9c6ff8c1752151f03 100644
--- a/grudge/models/em.py
+++ b/grudge/models/em.py
@@ -1,8 +1,5 @@
-# -*- coding: utf8 -*-
 """grudge operators modelling electromagnetic phenomena."""
 
-from __future__ import division, absolute_import
-
 __copyright__ = """
 Copyright (C) 2007-2017 Andreas Kloeckner
 Copyright (C) 2010 David Powell
@@ -186,7 +183,8 @@ class MaxwellOperator(HyperbolicOperator):
                 )
 
     def pec_bc(self, w):
-        "Construct part of the flux operator template for PEC boundary conditions"
+        """Construct part of the flux operator template for PEC boundary conditions
+        """
         e, h = self.split_eh(w)
 
         pec_e = sym.cse(sym.project("vol", self.pec_tag)(e))
@@ -195,7 +193,8 @@ class MaxwellOperator(HyperbolicOperator):
         return flat_obj_array(-pec_e, pec_h)
 
     def pmc_bc(self, w):
-        "Construct part of the flux operator template for PMC boundary conditions"
+        """Construct part of the flux operator template for PMC boundary conditions
+        """
         e, h = self.split_eh(w)
 
         pmc_e = sym.cse(sym.project("vol", self.pmc_tag)(e))
@@ -234,7 +233,7 @@ class MaxwellOperator(HyperbolicOperator):
         return bc
 
     def incident_bc(self, w):
-        "Flux terms for incident boundary conditions"
+        """Flux terms for incident boundary conditions"""
         # NOTE: Untested for inhomogeneous materials, but would usually be
         # physically meaningless anyway (are there exceptions to this?)
 
@@ -302,7 +301,7 @@ class MaxwellOperator(HyperbolicOperator):
             [e_subset, h_subset]))
 
     def split_eh(self, w):
-        "Splits an array into E and H components"
+        """Splits an array into E and H components"""
         e_idx, h_idx = self.partial_to_eh_subsets()
         e, h = w[e_idx], w[h_idx]
 
diff --git a/grudge/models/gas_dynamics/__init__.py b/grudge/models/gas_dynamics/__init__.py
index 4df408d811a826a1f924b855f0ff8a9fce8a29f2..8c86c8e0aed0c4cae043cec274f171122668f4cc 100644
--- a/grudge/models/gas_dynamics/__init__.py
+++ b/grudge/models/gas_dynamics/__init__.py
@@ -1,10 +1,5 @@
 """Operator for compressible Navier-Stokes and Euler equations."""
 
-from __future__ import division
-from __future__ import absolute_import
-import six
-from six.moves import range
-
 __copyright__ = "Copyright (C) 2007 Hendrik Riedmann, Andreas Kloeckner"
 
 __license__ = """
@@ -28,8 +23,6 @@ THE SOFTWARE.
 """
 
 
-
-
 import numpy
 import grudge.tools
 import grudge.mesh
@@ -65,7 +58,7 @@ to_bdry_quad = QuadratureGridUpsampler("gasdyn_face")
 
 
 # {{{ equations of state
-class EquationOfState(object):
+class EquationOfState:
     def q_to_p(self, op, q):
         raise NotImplementedError
 
@@ -264,9 +257,9 @@ class GasDynamicsOperator(TimeDependentOperator):
         return cse(self.temperature(q), "temperature")
 
     def get_mu(self, q, to_quad_op):
-        """
+        r"""
         :param to_quad_op: If not *None*, represents an operator which transforms
-          nodal values onto a quadrature grid on which the returned :math:`\mu` 
+          nodal values onto a quadrature grid on which the returned :math:`\mu`
           needs to be represented. In that case, *q* is assumed to already be on the
           same quadrature grid.
         """
@@ -280,7 +273,7 @@ class GasDynamicsOperator(TimeDependentOperator):
             t_s = 110.4
             mu_inf = 1.735e-5
             result = cse(
-                    mu_inf * self.cse_temperature(q) ** 1.5 * (1 + t_s) 
+                    mu_inf * self.cse_temperature(q) ** 1.5 * (1 + t_s)
                     / (self.cse_temperature(q) + t_s),
                     "sutherland_mu")
         else:
@@ -471,7 +464,7 @@ class GasDynamicsOperator(TimeDependentOperator):
 
                     # flux rho_u
                     make_obj_array([
-                        self.rho_u(q)[i]*self.cse_u(q)[j] 
+                        self.rho_u(q)[i]*self.cse_u(q)[j]
                         + delta(i,j) * self.cse_p(q)
                         for j in range(self.dimensions)
                         ])
@@ -508,11 +501,11 @@ class GasDynamicsOperator(TimeDependentOperator):
             rho0=rho0, p0=p0, u0=u0, c0=c0,
 
             # notation: suffix "m" for "minus", i.e. "interior"
-            drhom=cse(self.rho(cse(to_bdry_quad(bdrize_op(state)))) 
+            drhom=cse(self.rho(cse(to_bdry_quad(bdrize_op(state))))
                 - rho0, "drhom"),
-            dumvec=cse(self.cse_u(cse(to_bdry_quad(bdrize_op(state)))) 
+            dumvec=cse(self.cse_u(cse(to_bdry_quad(bdrize_op(state))))
                 - u0, "dumvec"),
-            dpm=cse(self.cse_p(cse(to_bdry_quad(bdrize_op(state)))) 
+            dpm=cse(self.cse_p(cse(to_bdry_quad(bdrize_op(state))))
                 - p0, "dpm"))
 
     def outflow_state(self, state):
@@ -602,7 +595,7 @@ class GasDynamicsOperator(TimeDependentOperator):
                 make_sym_vector("bc_q_supersonic_in", self.dimensions+2),
                 self.supersonic_outflow_tag:
                 RestrictToBoundary(self.supersonic_outflow_tag)(
-                            (state)),
+                            state),
                 self.wall_tag: self.wall_state(state),
                 }
 
@@ -659,7 +652,7 @@ class GasDynamicsOperator(TimeDependentOperator):
             from pymbolic.primitives import Subscript, Variable
 
             if isinstance(expr, Subscript):
-                assert (isinstance(expr.aggregate, Variable) 
+                assert (isinstance(expr.aggregate, Variable)
                         and expr.aggregate.name == "q")
 
                 return cbstate[expr.index]
@@ -697,7 +690,7 @@ class GasDynamicsOperator(TimeDependentOperator):
         faceq_tau_mat = self.tau(to_int_face_quad, state)
 
         return join_fields(
-                0, 
+                0,
                 self.div(
                     numpy.sum(volq_tau_mat*self.cse_u(volq_state), axis=1)
                     + self.heat_conduction_grad(to_vol_quad)
@@ -734,7 +727,7 @@ class GasDynamicsOperator(TimeDependentOperator):
             return make_obj_array([
                 self.div(
                     to_vol_quad(self.sensor())*to_vol_quad(dq[i]),
-                    to_int_face_quad(self.sensor())*to_int_face_quad(dq[i])) 
+                    to_int_face_quad(self.sensor())*to_int_face_quad(dq[i]))
                 for i in range(dq.shape[0])])
         # }}}
 
@@ -758,10 +751,10 @@ class GasDynamicsOperator(TimeDependentOperator):
 
         from grudge.symbolic.tools import make_stiffness_t
 
-        primitive_bcs_as_quad_conservative = dict(
-                (tag, self.primitive_to_conservative(to_bdry_quad(bc)))
-                for tag, bc in 
-                six.iteritems(self.get_primitive_boundary_conditions()))
+        primitive_bcs_as_quad_conservative = {
+                tag: self.primitive_to_conservative(to_bdry_quad(bc))
+                for tag, bc in
+                self.get_primitive_boundary_conditions().items()}
 
         def get_bc_tuple(tag):
             state = self.state()
@@ -783,7 +776,7 @@ class GasDynamicsOperator(TimeDependentOperator):
             first_order_part = 0*first_order_part
 
         result = join_fields(
-                first_order_part 
+                first_order_part
                 + self.make_second_order_part()
                 + make_artificial_diffusion()
                 + self.make_extra_terms(),
@@ -803,7 +796,7 @@ class GasDynamicsOperator(TimeDependentOperator):
 
     # {{{ operator binding ----------------------------------------------------
     def bind(self, discr, sensor=None, sensor_scaling=None, viscosity_only=False):
-        if (sensor is None and 
+        if (sensor is None and
                 self.artificial_viscosity_mode is not None):
             raise ValueError("must specify a sensor if using "
                     "artificial viscosity")
@@ -843,7 +836,7 @@ class GasDynamicsOperator(TimeDependentOperator):
                     bc_q_noslip=self.bc_noslip.boundary_interpolant(
                         t, discr, self.noslip_tag),
                     bc_q_supersonic_in=self.bc_supersonic_inflow
-                    .boundary_interpolant(t, discr, 
+                    .boundary_interpolant(t, discr,
                         self.supersonic_inflow_tag),
                     **extra_kwargs
                     )
@@ -858,10 +851,10 @@ class GasDynamicsOperator(TimeDependentOperator):
 
     # {{{ timestep estimation -------------------------------------------------
 
-    def estimate_timestep(self, discr, 
+    def estimate_timestep(self, discr,
             stepper=None, stepper_class=None, stepper_args=None,
             t=None, max_eigenvalue=None):
-        u"""Estimate the largest stable timestep, given a time stepper
+        """Estimate the largest stable timestep, given a time stepper
         `stepper_class`. If none is given, RK4 is assumed.
         """
 
diff --git a/grudge/models/gas_dynamics/lbm.py b/grudge/models/gas_dynamics/lbm.py
index 4eb65a61757b4e2fca38b05e5e118c5ab4fe8ee5..548e38aaeb0422442362c8f910b705e6c842e01d 100644
--- a/grudge/models/gas_dynamics/lbm.py
+++ b/grudge/models/gas_dynamics/lbm.py
@@ -1,11 +1,5 @@
-# -*- coding: utf8 -*-
 """Lattice-Boltzmann operator."""
 
-from __future__ import division
-from __future__ import absolute_import
-from six.moves import range
-from six.moves import zip
-
 __copyright__ = "Copyright (C) 2011 Andreas Kloeckner"
 
 __license__ = """
@@ -34,10 +28,7 @@ from grudge.models import HyperbolicOperator
 from pytools.obj_array import make_obj_array
 
 
-
-
-
-class LBMMethodBase(object):
+class LBMMethodBase:
     def __len__(self):
         return len(self.direction_vectors)
 
@@ -51,7 +42,7 @@ class LBMMethodBase(object):
             found = False
             for alpha_2 in range(alpha, len(self)):
                 if la.norm(
-                        self.direction_vectors[alpha] 
+                        self.direction_vectors[alpha]
                         + self.direction_vectors[alpha_2]) < 1e-12:
                     self.opposites[alpha] = alpha_2
                     self.opposites[alpha_2] = alpha
@@ -59,7 +50,7 @@ class LBMMethodBase(object):
 
             if not found:
                 raise RuntimeError(
-                        "direction %s had no opposite" 
+                        "direction %s had no opposite"
                         % self.direction_vectors[alpha])
 
 
@@ -197,9 +188,9 @@ class LatticeBoltzmannOperator(HyperbolicOperator):
         from grudge.symbolic.mappers.type_inference import (
                 type_info, NodalRepresentation)
 
-        type_hints = dict(
-                (f_bar_i, type_info.VolumeVector(NodalRepresentation()))
-                for f_bar_i in f_bar_sym)
+        type_hints = {
+                f_bar_i: type_info.VolumeVector(NodalRepresentation())
+                for f_bar_i in f_bar_sym}
 
         compiled_sym_operator = discr.compile(what(f_bar_sym), type_hints=type_hints)
 
diff --git a/grudge/models/nd_calculus.py b/grudge/models/nd_calculus.py
index 96be72343e02ea3e1496f42e3cfcbbfa75de2db1..bcfbd57d758723412dc9b4c73292759bab0ea2ae 100644
--- a/grudge/models/nd_calculus.py
+++ b/grudge/models/nd_calculus.py
@@ -1,9 +1,5 @@
-# -*- coding: utf8 -*-
 """Canned operators for multivariable calculus."""
 
-from __future__ import division
-from __future__ import absolute_import
-
 __copyright__ = "Copyright (C) 2009 Andreas Kloeckner"
 
 __license__ = """
@@ -26,14 +22,9 @@ OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
 THE SOFTWARE.
 """
 
-
-
-
 from grudge.models import Operator
 
 
-
-
 class GradientOperator(Operator):
     def __init__(self, dimensions):
         self.dimensions = dimensions
diff --git a/grudge/models/pml.py b/grudge/models/pml.py
index 8799a2c3bd93a1a9d54bb63be1187e27d22edd29..bd8bfd0468ebc96ea23cbd691de0eaa3dee8344e 100644
--- a/grudge/models/pml.py
+++ b/grudge/models/pml.py
@@ -1,11 +1,5 @@
-# -*- coding: utf8 -*-
 """Models describing absorbing boundary layers."""
 
-from __future__ import division
-from __future__ import absolute_import
-from six.moves import range
-from six.moves import zip
-
 __copyright__ = "Copyright (C) 2007 Andreas Kloeckner"
 
 __license__ = """
@@ -28,9 +22,6 @@ OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
 THE SOFTWARE.
 """
 
-
-
-
 import numpy
 
 from pytools import memoize_method, Record
@@ -40,8 +31,6 @@ from grudge.models.em import \
         TEMaxwellOperator
 
 
-
-
 class AbarbanelGottliebPMLMaxwellOperator(MaxwellOperator):
     """Implements a PML as in
 
@@ -68,8 +57,8 @@ class AbarbanelGottliebPMLMaxwellOperator(MaxwellOperator):
 
         def map(self, f):
             return self.__class__(
-                    **dict((name, f(getattr(self, name)))
-                        for name in self.fields))
+                    **{name: f(getattr(self, name))
+                        for name in self.fields})
 
     def __init__(self, *args, **kwargs):
         self.add_decay = kwargs.pop("add_decay", True)
diff --git a/grudge/models/poisson.py b/grudge/models/poisson.py
index a9f7145744c2c1a5e6e927680f341c3dbbf3e3e7..6636b94fe4d7ea90894bf92a86d30189fe93348e 100644
--- a/grudge/models/poisson.py
+++ b/grudge/models/poisson.py
@@ -1,9 +1,5 @@
-# -*- coding: utf8 -*-
 """Operators for Poisson problems."""
 
-from __future__ import division
-from __future__ import absolute_import
-
 __copyright__ = "Copyright (C) 2007 Andreas Kloeckner"
 
 __license__ = """
@@ -26,7 +22,6 @@ OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
 THE SOFTWARE.
 """
 
-
 import numpy as np
 
 from grudge.models import Operator
@@ -35,7 +30,7 @@ import grudge.data
 import grudge.iterative
 
 
-class LaplacianOperatorBase(object):
+class LaplacianOperatorBase:
     def sym_operator(self, apply_minv, u=None, dir_bc=None, neu_bc=None):
         """
         :param apply_minv: :class:`bool` specifying whether to compute a complete
diff --git a/grudge/models/second_order.py b/grudge/models/second_order.py
index 61aa72a097d914c2213c095d0556bf2738454ffb..8f7251d329ea1cbbb960e368fe3c272f481ca36a 100644
--- a/grudge/models/second_order.py
+++ b/grudge/models/second_order.py
@@ -1,8 +1,5 @@
-# -*- coding: utf8 -*-
 """Schemes for second-order derivatives."""
 
-from __future__ import division
-
 __copyright__ = "Copyright (C) 2009 Andreas Kloeckner"
 
 __license__ = """
@@ -25,7 +22,6 @@ OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
 THE SOFTWARE.
 """
 
-
 import numpy as np
 import grudge.symbolic.mappers as mappers
 from grudge import sym
@@ -35,10 +31,10 @@ from grudge import sym
 
 class StabilizationTermGenerator(mappers.IdentityMapper):
     def __init__(self, flux_args):
-        super(StabilizationTermGenerator, self).__init__()
+        super().__init__()
         self.flux_args = flux_args
-        self.flux_arg_lookup = dict(
-                (flux_arg, i) for i, flux_arg in enumerate(flux_args))
+        self.flux_arg_lookup = {
+                flux_arg: i for i, flux_arg in enumerate(flux_args)}
 
     def get_flux_arg_idx(self, expr, quad_above):
         from grudge.symbolic.mappers import QuadratureDetector
@@ -138,7 +134,7 @@ class StabilizationTermGenerator(mappers.IdentityMapper):
 
 class NeumannBCGenerator(mappers.IdentityMapper):
     def __init__(self, tag, bc):
-        super(NeumannBCGenerator, self).__init__()
+        super().__init__()
         self.tag = tag
         self.bc = bc
 
@@ -194,7 +190,7 @@ class IPDGDerivativeGenerator(mappers.IdentityMapper):
             return self.rec(expr.field)
         elif isinstance(expr.op,
                 sym.QuadratureInteriorFacesGridUpsampler):
-            return super(IPDGDerivativeGenerator, self).map_operator_binding(expr)
+            return super().map_operator_binding(expr)
         else:
             from grudge.symbolic.tools import pretty
             raise ValueError("IPDG derivative generator doesn't know "
@@ -203,7 +199,7 @@ class IPDGDerivativeGenerator(mappers.IdentityMapper):
 
 # {{{ second derivative target
 
-class SecondDerivativeTarget(object):
+class SecondDerivativeTarget:
     def __init__(self, dimensions, strong_form, operand,
             int_flux_operand=None,
             bdry_flux_int_operand=None):
@@ -319,7 +315,7 @@ class SecondDerivativeTarget(object):
 
 # {{{ second derivative schemes
 
-class SecondDerivativeBase(object):
+class SecondDerivativeBase:
     def grad(self, tgt, bc_getter, dirichlet_tags, neumann_tags):
         """
         :param bc_getter: a function (tag, volume_expr) -> boundary expr.
diff --git a/grudge/models/wave.py b/grudge/models/wave.py
index 863b0112052f71853fd746a80e29385b5ebf187f..a06295bb4743c671c9ce33c5bd00258733f96923 100644
--- a/grudge/models/wave.py
+++ b/grudge/models/wave.py
@@ -1,8 +1,5 @@
-# -*- coding: utf8 -*-
 """Wave equation operators."""
 
-from __future__ import division, absolute_import
-
 __copyright__ = "Copyright (C) 2009 Andreas Kloeckner"
 
 __license__ = """
@@ -25,7 +22,6 @@ OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
 THE SOFTWARE.
 """
 
-
 import numpy as np
 from grudge.models import HyperbolicOperator
 from meshmode.mesh import BTAG_ALL, BTAG_NONE
diff --git a/grudge/shortcuts.py b/grudge/shortcuts.py
index 94b8fb2ac6674667f8753aed501fbcc1806a1ffc..a4de4bb7a13d6dc7462f17698a23f858f7af6a2b 100644
--- a/grudge/shortcuts.py
+++ b/grudge/shortcuts.py
@@ -1,7 +1,5 @@
 """Minimal example of a grudge driver."""
 
-from __future__ import division, print_function
-
 __copyright__ = "Copyright (C) 2009 Andreas Kloeckner"
 
 __license__ = """
diff --git a/grudge/symbolic/__init__.py b/grudge/symbolic/__init__.py
index 59fa5b26b5943a24982db5616632cef7dc8c042f..c075674070b0138888b92c03a9a782f8344f9e90 100644
--- a/grudge/symbolic/__init__.py
+++ b/grudge/symbolic/__init__.py
@@ -1,7 +1,5 @@
 """Building blocks and mappers for operator expression trees."""
 
-from __future__ import division, absolute_import
-
 __copyright__ = "Copyright (C) 2008 Andreas Kloeckner"
 
 __license__ = """
diff --git a/grudge/symbolic/compiler.py b/grudge/symbolic/compiler.py
index e5b7f7e872541aacd7e352eb6e247e99915077c4..075c0b07f30375f3a010651d305005fd920c61e2 100644
--- a/grudge/symbolic/compiler.py
+++ b/grudge/symbolic/compiler.py
@@ -1,7 +1,5 @@
 """Compiler to turn operator expression tree into (imperative) bytecode."""
 
-from __future__ import division, absolute_import, print_function
-
 __copyright__ = "Copyright (C) 2008-15 Andreas Kloeckner"
 
 __license__ = """
@@ -26,16 +24,14 @@ THE SOFTWARE.
 
 import numpy as np
 
-import six  # noqa
-from six.moves import zip, reduce
-
 from pytools import Record, memoize_method, memoize
 from pytools.obj_array import obj_array_vectorize
 
 from grudge import sym
 import grudge.symbolic.mappers as mappers
 from pymbolic.primitives import Variable, Subscript
-from six.moves import intern
+from sys import intern
+from functools import reduce
 
 from loopy.version import LOOPY_USE_LANGUAGE_VERSION_2018_2  # noqa: F401
 
@@ -75,7 +71,7 @@ def _make_dep_mapper(include_subscripts):
 
 # {{{ loopy kernel instruction
 
-class LoopyKernelDescriptor(object):
+class LoopyKernelDescriptor:
     def __init__(self, loopy_kernel, input_mappings, output_mappings,
             fixed_arguments, governing_dd):
         self.loopy_kernel = loopy_kernel
@@ -97,14 +93,12 @@ class LoopyKernelInstruction(Instruction):
     scope_indicator = ""
 
     def __init__(self, kernel_descriptor):
-        super(LoopyKernelInstruction, self).__init__()
+        super().__init__()
         self.kernel_descriptor = kernel_descriptor
 
     @memoize_method
     def get_assignees(self):
-        return set(
-                k
-                for k in self.kernel_descriptor.output_mappings.keys())
+        return {k for k in self.kernel_descriptor.output_mappings.keys()}
 
     @memoize_method
     def get_dependencies(self):
@@ -119,7 +113,7 @@ class LoopyKernelInstruction(Instruction):
 
     def __str__(self):
         knl_str = "\n".join(
-                "%s = %s" % (insn.assignee, insn.expression)
+                f"{insn.assignee} = {insn.expression}"
                 for insn in self.kernel_descriptor.loopy_kernel.instructions)
 
         knl_str = knl_str.replace("grdg_", "")
@@ -141,7 +135,7 @@ class AssignBase(Instruction):
             if comment:
                 comment = "/* %s */ " % comment
 
-            return "%s <-%s %s%s" % (
+            return "{} <-{} {}{}".format(
                     self.names[0], self.scope_indicator, comment,
                     self.exprs[0])
         else:
@@ -156,7 +150,7 @@ class AssignBase(Instruction):
                 else:
                     dnr_indicator = ""
 
-                lines.append("  %s <%s-%s %s" % (
+                lines.append("  {} <{}-{} {}".format(
                     n, dnr_indicator, self.scope_indicator, e))
             lines.append("}")
             return "\n".join(lines)
@@ -197,7 +191,7 @@ class Assign(AssignBase):
                 for expr in self.exprs))
 
         from pymbolic.primitives import Variable
-        deps -= set(Variable(name) for name in self.names)
+        deps -= {Variable(name) for name in self.names}
 
         if not each_vector:
             self._dependencies = deps
@@ -228,11 +222,11 @@ class RankDataSwapAssign(Instruction):
         self.dd_out = op.dd_out
         self.send_tag = self.MPI_TAG_GRUDGE_DATA_BASE + op.unique_id
         self.recv_tag = self.MPI_TAG_GRUDGE_DATA_BASE + op.unique_id
-        self.comment = "Swap data with rank %02d" % self.i_remote_rank
+        self.comment = f"Swap data with rank {self.i_remote_rank:02d}"
 
     @memoize_method
     def get_assignees(self):
-        return set([self.name])
+        return {self.name}
 
     @memoize_method
     def get_dependencies(self):
@@ -240,10 +234,10 @@ class RankDataSwapAssign(Instruction):
 
     def __str__(self):
         return ("{\n"
-              + "   /* %s */\n" % self.comment
-              + "   send_tag = %s\n" % self.send_tag
-              + "   recv_tag = %s\n" % self.recv_tag
-              + "   %s <- %s\n" % (self.name, self.field)
+              + f"   /* {self.comment} */\n"
+              + f"   send_tag = {self.send_tag}\n"
+              + f"   recv_tag = {self.recv_tag}\n"
+              + f"   {self.name} <- {self.field}\n"
               + "}")
 
     mapper_method = intern("map_insn_rank_data_swap")
@@ -260,7 +254,7 @@ class FromDiscretizationScopedAssign(AssignBase):
     neglect_for_dofdesc_inference = True
 
     def __init__(self, name, **kwargs):
-        super(FromDiscretizationScopedAssign, self).__init__(name=name, **kwargs)
+        super().__init__(name=name, **kwargs)
 
     @memoize_method
     def flop_count(self):
@@ -305,11 +299,11 @@ class DiffBatchAssign(Instruction):
         if len(self.names) > 1:
             lines.append("{")
             for n, d in zip(self.names, self.operators):
-                lines.append("  %s <- %s(%s)" % (n, d, self.field))
+                lines.append(f"  {n} <- {d}({self.field})")
             lines.append("}")
         else:
             for n, d in zip(self.names, self.operators):
-                lines.append("%s <- %s(%s)" % (n, d, self.field))
+                lines.append(f"{n} <- {d}({self.field})")
 
         return "\n".join(lines)
 
@@ -326,8 +320,8 @@ def dot_dataflow_graph(code, max_node_label_length=30,
     node_names = {}
 
     result = [
-            "initial [label=\"initial\"]"
-            "result [label=\"result\"]"]
+            'initial [label="initial"]'
+            'result [label="result"]']
 
     for num, insn in enumerate(code.instructions):
         node_name = "node%d" % num
@@ -345,8 +339,8 @@ def dot_dataflow_graph(code, max_node_label_length=30,
 
         node_label = node_label.replace("\n", "\\l") + "\\l"
 
-        result.append("%s [ label=\"p%d: %s\" shape=box ];" % (
-            node_name, insn.priority, node_label))
+        result.append(f"{node_name} [ "
+                f'label="p{insn.priority}: {node_label}" shape=box ];')
 
         for assignee in insn.get_assignees():
             origins[assignee] = node_name
@@ -359,8 +353,8 @@ def dot_dataflow_graph(code, max_node_label_length=30,
             return "initial"
 
     def gen_expr_arrow(expr, target_node):
-        result.append("%s -> %s [label=\"%s\"];"
-                % (get_orig_node(expr), target_node, expr))
+        orig_node = get_orig_node(expr)
+        result.append(f'{orig_node} -> {target_node} [label="{expr}"];')
 
     for insn in code.instructions:
         for dep in insn.get_dependencies():
@@ -379,7 +373,7 @@ def dot_dataflow_graph(code, max_node_label_length=30,
 
 # {{{ code representation
 
-class Code(object):
+class Code:
     def __init__(self, instructions, result):
         self.instructions = instructions
         self.result = result
@@ -399,10 +393,10 @@ class Code(object):
             outf.write(dot_dataflow_graph(self, max_node_label_length=None))
 
     def __str__(self):
-        var_to_writer = dict(
-                (var_name, insn)
+        var_to_writer = {
+                var_name: insn
                 for insn in self.instructions
-                for var_name in insn.get_assignees())
+                for var_name in insn.get_assignees()}
 
         # {{{ topological sort
 
@@ -494,10 +488,10 @@ class Code(object):
             from time import time
             start_time = time()
             if profile_data == {}:
-                profile_data['insn_eval_time'] = 0
-                profile_data['future_eval_time'] = 0
-                profile_data['busy_wait_time'] = 0
-                profile_data['total_time'] = 0
+                profile_data["insn_eval_time"] = 0
+                profile_data["future_eval_time"] = 0
+                profile_data["busy_wait_time"] = 0
+                profile_data["total_time"] = 0
         if log_quantities is not None:
             exec_sub_timer = log_quantities["exec_timer"].start_sub_timer()
         context = exec_mapper.context
@@ -539,7 +533,7 @@ class Code(object):
 
                 futures.extend(new_futures)
                 if profile_data is not None:
-                    profile_data['insn_eval_time'] += time() - insn_start_time
+                    profile_data["insn_eval_time"] += time() - insn_start_time
                 if log_quantities is not None:
                     insn_sub_timer.stop().submit()
             except self.NoInstructionAvailable:
@@ -559,7 +553,7 @@ class Code(object):
                     for i in range(len(futures)):
                         if futures[i].is_ready():
                             if profile_data is not None:
-                                profile_data['busy_wait_time'] +=\
+                                profile_data["busy_wait_time"] +=\
                                         time() - busy_wait_start_time
                                 future_start_time = time()
                             if log_quantities is not None:
@@ -580,7 +574,7 @@ class Code(object):
                             did_eval_future = True
 
                             if profile_data is not None:
-                                profile_data['future_eval_time'] +=\
+                                profile_data["future_eval_time"] +=\
                                         time() - future_start_time
                             if log_quantities is not None:
                                 future_sub_timer.stop().submit()
@@ -593,7 +587,7 @@ class Code(object):
         if log_quantities is not None:
             exec_sub_timer.stop().submit()
         if profile_data is not None:
-            profile_data['total_time'] = time() - start_time
+            profile_data["total_time"] = time() - start_time
             return (obj_array_vectorize(exec_mapper, self.result),
                     profile_data)
         return obj_array_vectorize(exec_mapper, self.result)
@@ -642,8 +636,7 @@ def aggregate_assignments(inf_mapper, instructions, result,
         try:
             return var_assignees_cache[insn]
         except KeyError:
-            result = set(Variable(assignee)
-                    for assignee in insn.get_assignees())
+            result = {Variable(assignee) for assignee in insn.get_assignees()}
             var_assignees_cache[insn] = result
             return result
 
@@ -652,7 +645,7 @@ def aggregate_assignments(inf_mapper, instructions, result,
 
         from pymbolic.primitives import Variable
         deps = (ass_1.get_dependencies() | ass_2.get_dependencies()) \
-                - set(Variable(name) for name in names)
+                - {Variable(name) for name in names}
 
         return Assign(
                 names=names, exprs=ass_1.exprs + ass_2.exprs,
@@ -665,10 +658,10 @@ def aggregate_assignments(inf_mapper, instructions, result,
 
     insn_to_origins_cache = {}
 
-    origins_map = dict(
-                (assignee, insn)
+    origins_map = {
+                assignee: insn
                 for insn in instructions
-                for assignee in insn.get_assignees())
+                for assignee in insn.get_assignees()}
 
     from pytools import partition
     from grudge.symbolic.primitives import DTAG_SCALAR
@@ -760,26 +753,26 @@ def aggregate_assignments(inf_mapper, instructions, result,
         if not did_work:
             processed_assigns.append(my_assign)
 
-    externally_used_names = set(
+    externally_used_names = {
             expr
             for insn in processed_assigns + other_insns
-            for expr in insn.get_dependencies())
+            for expr in insn.get_dependencies()}
 
     if isinstance(result, np.ndarray) and result.dtype.char == "O":
-        externally_used_names |= set(expr for expr in result)
+        externally_used_names |= {expr for expr in result}
     else:
-        externally_used_names |= set([result])
+        externally_used_names |= {result}
 
     def schedule_and_finalize_assignment(ass):
         dep_mapper = _make_dep_mapper(include_subscripts=False)
 
         names_exprs = list(zip(ass.names, ass.exprs))
 
-        my_assignees = set(name for name, expr in names_exprs)
+        my_assignees = {name for name, expr in names_exprs}
         names_exprs_deps = [
                 (name, expr,
-                    set(dep.name for dep in dep_mapper(expr) if
-                        isinstance(dep, Variable)) & my_assignees)
+                    {dep.name for dep in dep_mapper(expr) if
+                        isinstance(dep, Variable)} & my_assignees)
                 for name, expr in names_exprs]
 
         ordered_names_exprs = []
@@ -855,7 +848,7 @@ class ToLoopyExpressionMapper(mappers.IdentityMapper):
     def map_name(self, name):
         dot_idx = name.find(".")
         if dot_idx != -1:
-            return "grdg_sub_%s_%s" % (name[:dot_idx], name[dot_idx+1:])
+            return "grdg_sub_{}_{}".format(name[:dot_idx], name[dot_idx+1:])
         else:
             return name
 
@@ -871,7 +864,7 @@ class ToLoopyExpressionMapper(mappers.IdentityMapper):
 
             suffix_nr = 0
             while name in self.used_names:
-                name = "%s_%s" % (name_prefix, suffix_nr)
+                name = f"{name_prefix}_{suffix_nr}"
                 suffix_nr += 1
             self.used_names.add(name)
 
@@ -982,7 +975,7 @@ def bessel_function_mangler(kernel, name, arg_dtypes):
 # }}}
 
 
-class ToLoopyInstructionMapper(object):
+class ToLoopyInstructionMapper:
     def __init__(self, dd_inference_mapper):
         self.dd_inference_mapper = dd_inference_mapper
         self.function_registry = dd_inference_mapper.function_registry
@@ -1065,7 +1058,7 @@ class ToLoopyInstructionMapper(object):
         from grudge.symbolic.mappers import DependencyMapper
         dep_mapper = DependencyMapper(composite_leaves=False)
 
-        for expr, name in six.iteritems(expr_mapper.expr_to_name):
+        for expr, name in expr_mapper.expr_to_name.items():
             deps = dep_mapper(expr)
             assert len(deps) <= 1
             if not deps:
@@ -1130,7 +1123,7 @@ class CodeGenerationState(Record):
 class OperatorCompiler(mappers.IdentityMapper):
     def __init__(self, discr, function_registry,
             prefix="_expr", max_vectors_in_batch_expr=None):
-        super(OperatorCompiler, self).__init__()
+        super().__init__()
         self.prefix = prefix
 
         self.max_vectors_in_batch_expr = max_vectors_in_batch_expr
@@ -1171,7 +1164,7 @@ class OperatorCompiler(mappers.IdentityMapper):
 
         codegen_state = CodeGenerationState(generating_discr_code=False)
         # Finally, walk the expression and build the code.
-        result = super(OperatorCompiler, self).__call__(expr, codegen_state)
+        result = super().__call__(expr, codegen_state)
 
         eval_code = self.eval_code
         del self.eval_code
@@ -1301,7 +1294,7 @@ class OperatorCompiler(mappers.IdentityMapper):
 
     def map_call(self, expr, codegen_state):
         if is_function_loopyable(expr.function, self.function_registry):
-            return super(OperatorCompiler, self).map_call(expr, codegen_state)
+            return super().map_call(expr, codegen_state)
         else:
             # If it's not a C-level function, it shouldn't get muddled up into
             # a vector math expression.
diff --git a/grudge/symbolic/dofdesc_inference.py b/grudge/symbolic/dofdesc_inference.py
index f5f11f0a22827299f340c5bed53d10d014e220f3..36f34fb3d30c78bb8c6e76669d83e919149c401e 100644
--- a/grudge/symbolic/dofdesc_inference.py
+++ b/grudge/symbolic/dofdesc_inference.py
@@ -1,5 +1,3 @@
-from __future__ import division, absolute_import
-
 __copyright__ = "Copyright (C) 2017 Andreas Kloeckner"
 
 __license__ = """
@@ -57,7 +55,7 @@ def unify_dofdescs(dd_a, dd_b, expr=None):
     return dd_a
 
 
-class InferrableMultiAssignment(object):
+class InferrableMultiAssignment:
     """An assignemnt 'instruction' which may be used as part of type
     inference.
 
@@ -86,11 +84,11 @@ class DOFDescInferenceMapper(RecursiveMapper, CSECachingMapperMixin):
 
         self.check = check
 
-        self.name_to_assignment = dict(
-                (name, a)
+        self.name_to_assignment = {
+                name: a
                 for a in assignments
                 if not a.neglect_for_dofdesc_inference
-                for name in a.get_assignees())
+                for name in a.get_assignees()}
 
         if name_to_dofdesc is None:
             name_to_dofdesc = {}
diff --git a/grudge/symbolic/mappers/__init__.py b/grudge/symbolic/mappers/__init__.py
index a6914e1fd56af9cb43210177a05e2006aed45d34..555b5ffe8358233bc1d10331a0d2f40f66c55c92 100644
--- a/grudge/symbolic/mappers/__init__.py
+++ b/grudge/symbolic/mappers/__init__.py
@@ -1,7 +1,5 @@
 """Mappers to transform symbolic operators."""
 
-from __future__ import division
-
 __copyright__ = "Copyright (C) 2008 Andreas Kloeckner"
 
 __license__ = """
@@ -24,9 +22,6 @@ OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
 THE SOFTWARE.
 """
 
-
-import six
-
 import numpy as np
 import pymbolic.primitives
 import pymbolic.mapper.stringifier
@@ -45,7 +40,7 @@ from grudge.tools import OrderedSet
 
 # {{{ mixins
 
-class LocalOpReducerMixin(object):
+class LocalOpReducerMixin:
     """Reduces calls to mapper methods for all local differentiation
     operators to a single mapper method, and likewise for mass
     operators.
@@ -107,7 +102,7 @@ class LocalOpReducerMixin(object):
     # }}}
 
 
-class FluxOpReducerMixin(object):
+class FluxOpReducerMixin:
     """Reduces calls to mapper methods for all flux
     operators to a smaller number of mapper methods.
     """
@@ -159,7 +154,7 @@ class OperatorReducerMixin(LocalOpReducerMixin, FluxOpReducerMixin):
     map_ref_face_mass_operator = _map_op_base
 
 
-class CombineMapperMixin(object):
+class CombineMapperMixin:
     def map_operator_binding(self, expr):
         return self.combine([self.rec(expr.op), self.rec(expr.field)])
 
@@ -227,7 +222,7 @@ class IdentityMapperMixin(LocalOpReducerMixin, FluxOpReducerMixin):
     # }}}
 
 
-class BoundOpMapperMixin(object):
+class BoundOpMapperMixin:
     def map_operator_binding(self, expr, *args, **kwargs):
         return getattr(self, expr.op.mapper_method)(
                 expr.op, expr.field, *args, **kwargs)
@@ -261,7 +256,7 @@ class DependencyMapper(
 
     def map_operator_binding(self, expr):
         if self.include_operator_bindings:
-            return set([expr])
+            return {expr}
         else:
             return CombineMapperMixin.map_operator_binding(self, expr)
 
@@ -269,7 +264,7 @@ class DependencyMapper(
         return set()
 
     def map_grudge_variable(self, expr):
-        return set([expr])
+        return {expr}
 
     def _map_leaf(self, expr):
         return set()
@@ -383,7 +378,7 @@ class OppositeInteriorFaceSwapUniqueIDAssigner(
     map_common_subexpression_uncached = IdentityMapper.map_common_subexpression
 
     def __init__(self):
-        super(OppositeInteriorFaceSwapUniqueIDAssigner, self).__init__()
+        super().__init__()
         self._next_id = 0
         self.seen_ids = set()
 
@@ -726,7 +721,9 @@ class StringifyMapper(pymbolic.mapper.stringifier.StringifyMapper):
         return result
 
     def _format_op_dd(self, op):
-        return ":%s->%s" % (self._format_dd(op.dd_in), self._format_dd(op.dd_out))
+        return ":{}->{}".format(
+                self._format_dd(op.dd_in),
+                self._format_dd(op.dd_out))
 
     # {{{ elementwise ops
 
@@ -801,7 +798,7 @@ class StringifyMapper(pymbolic.mapper.stringifier.StringifyMapper):
     # }}}
 
     def map_elementwise_linear(self, expr, enclosing_prec):
-        return "ElWLin:%s%s" % (
+        return "ElWLin:{}{}".format(
                 expr.__class__.__name__,
                 self._format_op_dd(expr))
 
@@ -836,12 +833,12 @@ class StringifyMapper(pymbolic.mapper.stringifier.StringifyMapper):
 
     def map_operator_binding(self, expr, enclosing_prec):
         from pymbolic.mapper.stringifier import PREC_NONE
-        return "<%s>(%s)" % (
+        return "<{}>({})".format(
                 self.rec(expr.op, PREC_NONE),
                 self.rec(expr.field, PREC_NONE))
 
     def map_grudge_variable(self, expr, enclosing_prec):
-        return "%s:%s" % (expr.name, self._format_dd(expr.dd))
+        return "{}:{}".format(expr.name, self._format_dd(expr.dd))
 
     def map_function_symbol(self, expr, enclosing_prec):
         return expr.name
@@ -1022,7 +1019,7 @@ class _InnerDerivativeJoiner(pymbolic.mapper.RecursiveMapper):
                 else:
                     return self.rec(expr, derivatives)
 
-            for operator, operands in six.iteritems(sub_derivatives):
+            for operator, operands in sub_derivatives.items():
                 for operand in operands:
                     derivatives.setdefault(operator, []).append(
                             factor*operand)
@@ -1072,7 +1069,7 @@ class DerivativeJoiner(CSECachingMapperMixin, IdentityMapper):
             if not sub_derivatives:
                 return expr
             else:
-                for operator, operands in six.iteritems(sub_derivatives):
+                for operator, operands in sub_derivatives.items():
                     derivatives.setdefault(operator, []).extend(operands)
 
                 return result
@@ -1081,7 +1078,7 @@ class DerivativeJoiner(CSECachingMapperMixin, IdentityMapper):
         new_children = [invoke_idj(child)
                 for child in expr.children]
 
-        for operator, operands in six.iteritems(derivatives):
+        for operator, operands in derivatives.items():
             new_children.insert(0, operator(
                 sum(self.rec(operand) for operand in operands)))
 
@@ -1295,9 +1292,9 @@ class SymbolicEvaluator(pymbolic.mapper.evaluator.EvaluationMapper):
                 expr.function,
                 tuple(self.rec(child, *args, **kwargs)
                     for child in expr.parameters),
-                dict(
-                    (key, self.rec(val, *args, **kwargs))
-                    for key, val in six.iteritems(expr.kw_parameters))
+                {
+                    key: self.rec(val, *args, **kwargs)
+                    for key, val in expr.kw_parameters.items()}
                     )
 
     def map_common_subexpression(self, expr):
diff --git a/grudge/symbolic/operators.py b/grudge/symbolic/operators.py
index 6a5a982ce29cf3a1168b1c3851f1c9a7195f59ec..5155cbb6ad213691f56054001b697e505db40878 100644
--- a/grudge/symbolic/operators.py
+++ b/grudge/symbolic/operators.py
@@ -1,5 +1,3 @@
-from __future__ import division, absolute_import
-
 __copyright__ = "Copyright (C) 2008-2017 Andreas Kloeckner, Bogdan Enache"
 
 __license__ = """
@@ -22,7 +20,7 @@ OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
 THE SOFTWARE.
 """
 
-from six.moves import intern
+from sys import intern
 
 import numpy as np
 import pymbolic.primitives
@@ -150,7 +148,7 @@ class ElementwiseLinearOperator(Operator):
 
 class ProjectionOperator(Operator):
     def __init__(self, dd_in, dd_out):
-        super(ProjectionOperator, self).__init__(dd_in, dd_out)
+        super().__init__(dd_in, dd_out)
 
     def __call__(self, expr):
         from pytools.obj_array import obj_array_vectorize
@@ -181,7 +179,7 @@ class InterpolationOperator(ProjectionOperator):
                 "use 'ProjectionOperator' instead.",
                 DeprecationWarning, stacklevel=2)
 
-        super(InterpolationOperator, self).__init__(dd_in, dd_out)
+        super().__init__(dd_in, dd_out)
 
 
 def interp(dd_in, dd_out):
@@ -196,7 +194,7 @@ def interp(dd_in, dd_out):
 
 class ElementwiseReductionOperator(Operator):
     def __init__(self, dd):
-        super(ElementwiseReductionOperator, self).__init__(dd_in=dd, dd_out=dd)
+        super().__init__(dd_in=dd, dd_out=dd)
 
 
 class ElementwiseSumOperator(ElementwiseReductionOperator):
@@ -235,7 +233,7 @@ class NodalReductionOperator(Operator):
 
         assert dd_out.is_scalar()
 
-        super(NodalReductionOperator, self).__init__(dd_out=dd_out, dd_in=dd_in)
+        super().__init__(dd_out=dd_out, dd_in=dd_in)
 
 
 class NodalSum(NodalReductionOperator):
@@ -271,7 +269,7 @@ class DiffOperatorBase(Operator):
             raise ValueError("differentiation outputs are not on "
                     "quadrature grids")
 
-        super(DiffOperatorBase, self).__init__(dd_in, dd_out)
+        super().__init__(dd_in, dd_out)
 
         self.xyz_axis = xyz_axis
 
@@ -327,7 +325,7 @@ class RefDiffOperatorBase(ElementwiseLinearOperator):
             raise ValueError("differentiation outputs are not on "
                     "quadrature grids")
 
-        super(RefDiffOperatorBase, self).__init__(dd_in, dd_out)
+        super().__init__(dd_in, dd_out)
 
         self.rst_axis = rst_axis
 
@@ -371,7 +369,7 @@ class RefStiffnessTOperator(RefDiffOperatorBase):
             grad_vand = (grad_vand,)
 
         weights = in_elem_grp.weights
-        return np.einsum('c,bz,acz->abc', weights, vand_inv_t, grad_vand)
+        return np.einsum("c,bz,acz->abc", weights, vand_inv_t, grad_vand)
 
 
 # }}}
@@ -402,7 +400,7 @@ class FilterOperator(ElementwiseLinearOperator):
         if dd_in != dd_out:
             raise ValueError("dd_in and dd_out must be identical")
 
-        super(FilterOperator, self).__init__(dd_in, dd_out)
+        super().__init__(dd_in, dd_out)
 
         self.mode_response_func = mode_response_func
 
@@ -490,7 +488,7 @@ class MassOperatorBase(Operator):
         if dd_out is None:
             dd_out = dd_in
 
-        super(MassOperatorBase, self).__init__(dd_in, dd_out)
+        super().__init__(dd_in, dd_out)
 
 
 class MassOperator(MassOperatorBase):
@@ -517,7 +515,7 @@ class RefMassOperator(RefMassOperatorBase):
         vand_inv_t = np.linalg.inv(vand).T
 
         weights = in_element_group.weights
-        return np.einsum('j,ik,jk->ij', weights, vand_inv_t, o_vand)
+        return np.einsum("j,ik,jk->ij", weights, vand_inv_t, o_vand)
 
     mapper_method = intern("map_ref_mass")
 
@@ -557,7 +555,7 @@ class OppositeInteriorFaceSwap(Operator):
         if dd_out is None:
             dd_out = dd_in
 
-        super(OppositeInteriorFaceSwap, self).__init__(dd_in, dd_out)
+        super().__init__(dd_in, dd_out)
         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:
@@ -592,7 +590,7 @@ class OppositePartitionFaceSwap(Operator):
         elif dd_out is None:
             dd_out = dd_in
 
-        super(OppositePartitionFaceSwap, self).__init__(dd_in, dd_out)
+        super().__init__(dd_in, dd_out)
         if not (isinstance(self.dd_in.domain_tag, prim.DTAG_BOUNDARY)
                 and isinstance(self.dd_in.domain_tag.tag, prim.BTAG_PARTITION)):
             raise ValueError(
@@ -632,7 +630,7 @@ class FaceMassOperatorBase(ElementwiseLinearOperator):
         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)
+        super().__init__(dd_in, dd_out)
 
 
 class FaceMassOperator(FaceMassOperatorBase):
diff --git a/grudge/symbolic/primitives.py b/grudge/symbolic/primitives.py
index 5dab6139b5e360589f10034548b23d5537d8bd0d..b106de1cb2b67eceb916e4031d3d3f7462b0f081 100644
--- a/grudge/symbolic/primitives.py
+++ b/grudge/symbolic/primitives.py
@@ -1,7 +1,5 @@
 """Operator template language: primitives."""
 
-from __future__ import division, absolute_import
-
 __copyright__ = "Copyright (C) 2008-2017 Andreas Kloeckner, Bogdan Enache"
 
 __license__ = """
@@ -24,7 +22,7 @@ OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
 THE SOFTWARE.
 """
 
-from six.moves import range, intern
+from sys import intern
 
 import numpy as np
 from pytools.obj_array import make_obj_array
@@ -145,14 +143,14 @@ class DTAG_BOUNDARY:        # noqa: N801
         return hash(type(self)) ^ hash(self.tag)
 
     def __repr__(self):
-        return "<%s(%s)>" % (type(self).__name__, repr(self.tag))
+        return "<{}({})>".format(type(self).__name__, repr(self.tag))
 
 
 class QTAG_NONE:            # noqa: N801
     pass
 
 
-class DOFDesc(object):
+class DOFDesc:
     """Describes the meaning of degrees of freedom.
 
     .. attribute:: domain_tag
@@ -275,7 +273,9 @@ class DOFDesc(object):
             else:
                 return repr(s)
 
-        return "DOFDesc(%s, %s)" % (fmt(self.domain_tag), fmt(self.quadrature_tag))
+        return "DOFDesc({}, {})".format(
+                fmt(self.domain_tag),
+                fmt(self.quadrature_tag))
 
 
 DD_SCALAR = DOFDesc(DTAG_SCALAR, None)
@@ -293,7 +293,7 @@ def as_dofdesc(dd):
 
 # {{{ has-dof-desc mix-in
 
-class HasDOFDesc(object):
+class HasDOFDesc:
     """
     .. attribute:: dd
 
@@ -310,7 +310,7 @@ class HasDOFDesc(object):
             dd = args[-1]
             args = args[:-1]
 
-        super(HasDOFDesc, self).__init__(*args, **kwargs)
+        super().__init__(*args, **kwargs)
         self.dd = dd
 
     def __getinitargs__(self):
@@ -339,7 +339,7 @@ class Variable(HasDOFDesc, ExpressionBase, VariableBase):
         if dd is None:
             dd = DD_VOLUME
 
-        super(Variable, self).__init__(name, dd)
+        super().__init__(name, dd)
 
     def __getinitargs__(self):
         return (self.name, self.dd,)
@@ -352,7 +352,7 @@ var = Variable
 
 class ScalarVariable(Variable):
     def __init__(self, name):
-        super(ScalarVariable, self).__init__(name, DD_SCALAR)
+        super().__init__(name, DD_SCALAR)
 
 
 def make_sym_array(name, shape, dd=None):
@@ -375,7 +375,7 @@ class FunctionSymbol(ExpressionBase, VariableBase):
     def __call__(self, *exprs):
         from pytools.obj_array import obj_array_vectorize_n_args
         return obj_array_vectorize_n_args(
-                super(FunctionSymbol, self).__call__, *exprs)
+                super().__call__, *exprs)
 
     mapper_method = "map_function_symbol"
 
@@ -426,7 +426,7 @@ class PrioritizedSubexpression(CommonSubexpressionBase):
     """
 
     def __init__(self, child, priority=0):
-        super(PrioritizedSubexpression, self).__init__(child)
+        super().__init__(child)
         self.priority = priority
 
     def __getinitargs__(self):
@@ -459,7 +459,7 @@ class _SignedFaceOnes(HasDOFDesc, ExpressionBase):
     def __init__(self, dd):
         dd = as_dofdesc(dd)
         assert dd.is_trace()
-        super(_SignedFaceOnes, self).__init__(dd)
+        super().__init__(dd)
 
     mapper_method = intern("map_signed_face_ones")
 
@@ -476,7 +476,7 @@ class NodeCoordinateComponent(DiscretizationProperty):
             raise ValueError("dd must be a discretization for "
                     "NodeCoordinateComponent")
 
-        super(NodeCoordinateComponent, self).__init__(dd)
+        super().__init__(dd)
         self.axis = axis
 
         assert dd.domain_tag is not None
diff --git a/grudge/symbolic/tools.py b/grudge/symbolic/tools.py
index 55a654d619e9f908f851482388191732d9774bb7..12d8c47fe66a473e2d380c83e3ea09c95b0dd1aa 100644
--- a/grudge/symbolic/tools.py
+++ b/grudge/symbolic/tools.py
@@ -1,7 +1,5 @@
 """Operator templates: extra bits of functionality."""
 
-from __future__ import division, absolute_import
-
 __copyright__ = "Copyright (C) 2008 Andreas Kloeckner"
 
 __license__ = """
@@ -24,8 +22,6 @@ OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
 THE SOFTWARE.
 """
 
-
-from six.moves import range
 import numpy as np
 
 
diff --git a/grudge/tools.py b/grudge/tools.py
index 67107e48048b4e78727a97e5d828f9a53cffdcad..4872962e521c5b2a512111589c699b7a11bb28f0 100644
--- a/grudge/tools.py
+++ b/grudge/tools.py
@@ -1,7 +1,5 @@
 """Miscellaneous helper facilities."""
 
-from __future__ import division
-
 __copyright__ = "Copyright (C) 2007 Andreas Kloeckner"
 
 __license__ = """
@@ -29,7 +27,7 @@ from pytools import levi_civita
 
 
 def is_zero(x):
-    # DO NOT try to replace this with an attempted '== 0' comparison.
+    # DO NOT try to replace this with an attempted "== 0" comparison.
     # This will become an elementwise numpy operation and not do what
     # you want.
 
@@ -185,15 +183,15 @@ class OrderedSet(MutableSet):
 
     def pop(self, last=True):
         if not self:
-            raise KeyError('set is empty')
+            raise KeyError("set is empty")
         key = self.end[1][0] if last else self.end[2][0]
         self.discard(key)
         return key
 
     def __repr__(self):
         if not self:
-            return '%s()' % (self.__class__.__name__,)
-        return '%s(%r)' % (self.__class__.__name__, list(self))
+            return f"{self.__class__.__name__}()"
+        return "{}({!r})".format(self.__class__.__name__, list(self))
 
     def __eq__(self, other):
         if isinstance(other, OrderedSet):
diff --git a/setup.cfg b/setup.cfg
index b2d98b73fe5e4ba7299a6e5aa4ce24b5e9591691..e8791cc48f86a67e726046a2f7510f32a39a98a5 100644
--- a/setup.cfg
+++ b/setup.cfg
@@ -9,6 +9,11 @@ exclude=
   grudge/models/nd_calculus.py,
   grudge/dt_finding.py
 
+
+inline-quotes = "
+docstring-quotes = """
+multiline-quotes = """
+
 [mypy]
 ignore_missing_imports = True
 
diff --git a/setup.py b/setup.py
index dd85b0dbc1f6b20bdae108442ba1bea58a7f2f1f..22aaa32c258a5e73d2bff6cee64490fa4403478d 100644
--- a/setup.py
+++ b/setup.py
@@ -54,8 +54,6 @@ def main():
               "cgen>=2013.1.2",
               "leap>=2019.1",
               "dagrt>=2019.1",
-
-              "six>=1.8",
               ])
 
 
diff --git a/test/mesh_data.py b/test/mesh_data.py
index 2ea2cc4782e8157485ed60991e118ab7363d2172..14cee46d54aa1fbb823479d51ede6c5c1e419fe2 100644
--- a/test/mesh_data.py
+++ b/test/mesh_data.py
@@ -1,4 +1,3 @@
-import six
 import numpy as np
 
 
@@ -7,7 +6,7 @@ class MeshBuilder:
     mesh_order = None
 
     def __init__(self, **kwargs):
-        for k, v in six.iteritems(kwargs):
+        for k, v in kwargs.items():
             setattr(self, k, v)
 
         if self.mesh_order is None:
diff --git a/test/test_grudge.py b/test/test_grudge.py
index 3ce0f8ae89f00cd6d60fa660f79b76be0d895809..9fd6ca3d1cc7935ef25aa4db0b180ef16fe3920b 100644
--- a/test/test_grudge.py
+++ b/test/test_grudge.py
@@ -1,5 +1,3 @@
-from __future__ import division, absolute_import, print_function
-
 __copyright__ = "Copyright (C) 2015 Andreas Kloeckner"
 
 __license__ = """
diff --git a/test/test_mpi_communication.py b/test/test_mpi_communication.py
index 91656bda4d1e123cc6a69eafe0c2f1a9b4eb1fa4..4ed9fdb2e0d58c7bfc3df069abea34ad097174b6 100644
--- a/test/test_mpi_communication.py
+++ b/test/test_mpi_communication.py
@@ -1,5 +1,3 @@
-from __future__ import division, absolute_import, print_function
-
 __copyright__ = """
 Copyright (C) 2017 Ellis Hoag
 Copyright (C) 2017 Andreas Kloeckner
@@ -170,7 +168,7 @@ def mpi_communication_entrypoint():
             IntervalTimer, EventCounter
     log_filename = None
     # NOTE: LogManager hangs when using a file on a shared directory.
-    # log_filename = 'grudge_log.dat'
+    # log_filename = "grudge_log.dat"
     logmgr = LogManager(log_filename, "w", comm)
     add_run_info(logmgr)
     add_general_quantities(logmgr)
@@ -243,10 +241,10 @@ def mpi_communication_entrypoint():
             \tBusy Wait: %g\n
             \tTotal: %g seconds""",
             i_local_rank,
-            data['insn_eval_time'] / data['total_time'] * 100,
-            data['future_eval_time'] / data['total_time'] * 100,
-            data['busy_wait_time'] / data['total_time'] * 100,
-            data['total_time'])
+            data["insn_eval_time"] / data["total_time"] * 100,
+            data["future_eval_time"] / data["total_time"] * 100,
+            data["busy_wait_time"] / data["total_time"] * 100,
+            data["total_time"])
 
     print_profile_data(rhs.profile_data)
     logmgr.close()