diff --git a/grudge/symbolic/mappers/__init__.py b/grudge/symbolic/mappers/__init__.py
index ca11cbf0d92c96b6c9d00a6e035b7b840c2a85cd..e533fcd7f0182e04b61d63ccf34cf73f0b159480 100644
--- a/grudge/symbolic/mappers/__init__.py
+++ b/grudge/symbolic/mappers/__init__.py
@@ -515,19 +515,33 @@ class GlobalToReferenceMapper(CSECachingMapperMixin, IdentityMapper):
             if with_jacobian:
                 rec_field = jac_tag * rec_field
             return sum(
+                ref_class(rst_axis, dd_in=dd_in)(
                     sym.inverse_metric_derivative(
                         rst_axis, expr.op.xyz_axis,
-                        ambient_dim=self.ambient_dim, dim=self.dim) *
-                    ref_class(rst_axis, dd_in=dd_in)(rec_field)
-                    for rst_axis in range(self.dim))
+                        ambient_dim=self.ambient_dim, dim=self.dim, dd=dd_in)
+                    * rec_field)
+                for rst_axis in range(self.dim))
 
         if isinstance(expr.op, op.MassOperator):
             return op.RefMassOperator(expr.op.dd_in, expr.op.dd_out)(
                     jac_in * self.rec(expr.field))
 
         elif isinstance(expr.op, op.InverseMassOperator):
-            return op.RefInverseMassOperator(expr.op.dd_in, expr.op.dd_out)(
-                1/jac_in * self.rec(expr.field))
+            # follows
+            #
+            # Chan, Hewett, and Warburton. "Weight-Adjusted Discontinuous
+            # Galerkin Methods: Curvilinear Meshes."
+            # SIAM Journal on Scientific Computing
+            #
+            # https://doi.org/10.1137/16M1089198
+
+            in_inv_op = op.RefInverseMassOperator(expr.op.dd_in, expr.op.dd_in)
+            trans_mass_op = op.RefMassOperator(expr.op.dd_in, expr.op.dd_out)
+            out_inv_op = op.RefInverseMassOperator(expr.op.dd_out, expr.op.dd_out)
+            return out_inv_op(
+                trans_mass_op(in_inv_op(self.rec(expr.field))
+                /
+                jac_in))
 
         elif isinstance(expr.op, op.FaceMassOperator):
             jac_in_surf = sym.area_element(self.ambient_dim, self.dim - 1,
diff --git a/test/test_grudge.py b/test/test_grudge.py
index 34532bfc80371e29bde28ea9674c48e36ab889aa..3976f246c6f1714a10730695df4cc52850657ed6 100644
--- a/test/test_grudge.py
+++ b/test/test_grudge.py
@@ -333,6 +333,80 @@ def test_convergence_advec(ctx_factory, mesh_name, mesh_pars, op_type, flux_type
     assert eoc_rec.order_estimate() > order
 
 
+@pytest.mark.parametrize("order", [3, 4, 5])
+def test_convergence_curvi_advection(ctx_factory,  order):
+    """Test whether a curvilinear mesh converges"""
+
+    cl_ctx = cl.create_some_context()
+    queue = cl.CommandQueue(cl_ctx)
+
+    from pytools.convergence import EOCRecorder
+    eoc_rec = EOCRecorder()
+
+    dims = 2
+    ns = [4, 6, 8]
+    for n in ns:
+
+        from meshmode.mesh.generation import generate_warped_rect_mesh
+        mesh = generate_warped_rect_mesh(dims, order, n)
+
+        discr = DGDiscretizationWithBoundaries(cl_ctx, mesh, order=order)
+
+        advec_v = np.array([0.1, 0.1])
+        norm_v = la.norm(advec_v)
+
+        def f(x):
+            return sym.sin(3*x)
+
+        def u_analytic(x):
+            return f(-np.dot(advec_v, x)/norm_v+sym.var("t", sym.DD_SCALAR)*norm_v)
+
+        from grudge.models.advection import WeakAdvectionOperator
+        analytic_sol = bind(discr, u_analytic(sym.nodes(dims)))
+
+        analytic_init = bind(discr, u_analytic(sym.nodes(dims)))
+
+        fields = analytic_init(queue, t=0)
+
+        op = WeakAdvectionOperator(advec_v,
+            inflow_u=u_analytic(sym.nodes(dims, sym.BTAG_ALL)),
+            flux_type="central")
+
+        bound_op = bind(discr, op.sym_operator())
+
+        def rhs(t, u):
+            return bound_op(queue, t=t, u=u)
+
+        final_t = 0.1
+
+        h = 1.0/n
+        dt_factor = 4
+        dt = dt_factor * h/order**2
+        nsteps = (final_t // dt) + 1
+        dt = final_t/nsteps + 1e-15
+
+        from grudge.shortcuts import set_up_rk4
+        dt_stepper = set_up_rk4("u", dt, fields, rhs)
+
+        norm = bind(discr, sym.norm(2, sym.var("u")))
+
+        step = 0
+        for event in dt_stepper.run(t_end=final_t):
+            if isinstance(event, dt_stepper.StateComputed):
+                assert event.component_id == "u"
+                esc = event.state_component
+                step += 1
+
+        sol = analytic_sol(queue, t=step * dt)
+        total_error = norm(queue, u=(esc - sol)) / norm(queue, u=sol)
+        eoc_rec.add_data_point(1.0/n, total_error)
+
+    print(eoc_rec.pretty_print(abscissa_label="h",
+            error_label="L2 Error"))
+
+    assert eoc_rec.order_estimate() > order
+
+
 @pytest.mark.parametrize("order", [3, 4, 5])
 def test_convergence_maxwell(ctx_factory,  order):
     """Test whether 3D maxwells actually converges"""