diff --git a/src/faceGradientIntegration.cpp b/src/faceGradientIntegration.cpp
index 6d67cce56..2b19e2f67 100644
--- a/src/faceGradientIntegration.cpp
+++ b/src/faceGradientIntegration.cpp
@@ -60,18 +60,17 @@ void GradFaceIntegrator::AssembleFaceVector(const FiniteElement &el1, const Fini
   elvect.SetSize((dof1 + dof2) * num_equation * dim);
   elvect = 0.0;
 
-// #ifdef _GPU_
-//   DenseMatrix elfun1_mat(elfun.GetData(), dof1, num_equation);
-//   DenseMatrix elfun2_mat(elfun.GetData()+dof1*num_equation,
-//                          dof2, num_equation);
-#ifndef _GPU_
-  DenseMatrix elfun1_mat(elfun.GetData(), dof1, num_equation * dim);
-  DenseMatrix elfun2_mat(elfun.GetData() + dof1 * num_equation * dim, dof2, num_equation * dim);
+  // NB: Incoming Vector &elfun is in gradient space.  The first
+  // component is the state, and the rest is zero.  See
+  // GradNonLinearForm::Mult in gradNonLinearForm.cpp for details.
+  DenseMatrix elfun1_mat(elfun.GetData(), dof1, num_equation);
+  DenseMatrix elfun2_mat(elfun.GetData() + dof1 * num_equation * dim, dof2, num_equation);
+
+
   DenseMatrix elvect1_mat(elvect.GetData(), dof1, num_equation * dim);
   DenseMatrix elvect2_mat(elvect.GetData() + dof1 * num_equation * dim, dof2, num_equation * dim);
   elvect1_mat = 0.;
   elvect2_mat = 0.;
-#endif
 
   // Integration order calculation from DGTraceIntegrator
   int intorder;
@@ -104,20 +103,12 @@ void GradFaceIntegrator::AssembleFaceVector(const FiniteElement &el1, const Fini
     for (int eq = 0; eq < num_equation; eq++) {
       double sum = 0.;
       for (int k = 0; k < dof1; k++) {
-#ifdef _GPU_
-        sum += shape1(k) * elfun(k + eq * dof1);
-#else
         sum += shape1[k] * elfun1_mat(k, eq);
-#endif
       }
       iUp1(eq) = sum;
       sum = 0.;
       for (int k = 0; k < dof2; k++) {
-#ifdef _GPU_
-        sum += shape2(k) * elfun(k + eq * dof2 + dof1 * num_equation);
-#else
         sum += shape2[k] * elfun2_mat(k, eq);
-#endif
       }
       iUp2(eq) = sum;
       mean(eq) = (iUp1(eq) + iUp2(eq)) / 2.;
@@ -133,18 +124,10 @@ void GradFaceIntegrator::AssembleFaceVector(const FiniteElement &el1, const Fini
         const double du1n = (mean(eq) - iUp1(eq)) * nor(d);
         const double du2n = (mean(eq) - iUp2(eq)) * nor(d);
         for (int k = 0; k < dof1; k++) {
-#ifdef _GPU_
-          elvect(k + eq + d * num_equation) += du1n * shape1(k);
-#else
           elvect1_mat(k, eq + d * num_equation) += du1n * shape1(k);
-#endif
         }
         for (int k = 0; k < dof2; k++) {
-#ifdef _GPU_
-          elvect(k + eq + d * num_equation + dof1 * num_equation * dim) -= du2n * shape2(k);
-#else
           elvect2_mat(k, eq + d * num_equation) -= du2n * shape2(k);
-#endif
         }
       }
     }
diff --git a/src/gradNonLinearForm.cpp b/src/gradNonLinearForm.cpp
index 87ad42ada..243a645a1 100644
--- a/src/gradNonLinearForm.cpp
+++ b/src/gradNonLinearForm.cpp
@@ -39,82 +39,14 @@ GradNonLinearForm::GradNonLinearForm(ParFiniteElementSpace *_vfes, IntegrationRu
 
 void GradNonLinearForm::Mult(const ParGridFunction *Up, Vector &y) {
   Vector x;
-
-#ifdef _GPU_
-  //   TestNonLinearForm::setToZero_gpu(y,y.Size()); already done in gradient.cpd
-  x.UseDevice(true);
-  if (fnfi.Size()) {
-    x.NewMemoryAndSize(Up->GetMemory(), vfes->GetNDofs() * num_equation, false);
-    Mult_gpu(x, y);
-  }
-#else
   const double *data = Up->GetData();
+
+  // NB: Setting x here accounts for the fact that the space Up lives
+  // in is different from the space of the gradient.
   x.SetSize(vfes->GetNDofs() * num_equation * dim);
   x = 0.;
   for (int i = 0; i < vfes->GetNDofs() * num_equation; i++) x(i) = data[i];
-  ParNonlinearForm::Mult(x, y);
-#endif
-}
-
-#ifdef _GPU_
-void GradNonLinearForm::Mult_gpu(const Vector &x, Vector &y) {
-  // INTEGRATION SHARED FACES
-  const Vector &px = Prolongate(x);
-
-  const int dofs = vfes->GetNDofs();
-
-  // Terms over shared interior faces in parallel.
-  ParFiniteElementSpace *pfes = ParFESpace();
-  ParMesh *pmesh = pfes->GetParMesh();
-  FaceElementTransformations *tr;
-  const FiniteElement *fe1, *fe2;
-  Array<int> vdofs1, vdofs2;
-  Vector el_x, el_y;
-  el_x.UseDevice(true);
-  el_y.UseDevice(true);
-
-  // if (bfnfi.Size()==0) y.HostReadWrite();
-  X.MakeRef(aux1, 0);  // aux1 contains P.x
-  X.ExchangeFaceNbrData();
-  const int n_shared_faces = pmesh->GetNSharedFaces();
-  for (int i = 0; i < n_shared_faces; i++) {
-    tr = pmesh->GetSharedFaceTransformations(i, true);
-    int Elem2NbrNo = tr->Elem2No - pmesh->GetNE();
 
-    fe1 = pfes->GetFE(tr->Elem1No);
-    fe2 = pfes->GetFaceNbrFE(Elem2NbrNo);
-
-    pfes->GetElementVDofs(tr->Elem1No, vdofs1);
-    pfes->GetFaceNbrElementVDofs(Elem2NbrNo, vdofs2);
-
-    el_x.SetSize(vdofs1.Size() + vdofs2.Size());
-    X.GetSubVector(vdofs1, el_x.GetData());
-    X.FaceNbrData().GetSubVector(vdofs2, el_x.GetData() + vdofs1.Size());
-
-    const int elDof1 = fe1->GetDof();
-    const int elDof2 = fe2->GetDof();
-    for (int k = 0; k < fnfi.Size(); k++) {
-      fnfi[k]->AssembleFaceVector(*fe1, *fe2, *tr, el_x, el_y);
-      // y.AddElementVector(vdofs1, el_y.GetData());
-      GradNonLinearForm::IndexedAddToGlobalMemory(vdofs1, y, el_y, dofs, elDof1, num_equation, dim);
-    }
-  }
-}
-
-void GradNonLinearForm::IndexedAddToGlobalMemory(const Array<int> &vdofs, Vector &y, const Vector &el_y, const int &dof,
-                                                 const int &elDof, const int &num_equation, const int &dim) {
-  double *dy = y.ReadWrite();
-  const double *d_ely = el_y.Read();
-  auto d_vdofs = vdofs.Read();
-
-  MFEM_FORALL(i, elDof, {
-    const int index = d_vdofs[i];
-    for (int eq = 0; eq < num_equation; eq++) {
-      for (int d = 0; d < dim; d++) {
-        dy[index + eq * dof + d * num_equation * dof] += d_ely[i + eq * elDof + d * num_equation * elDof];
-      }
-    }
-  });
+  // After setting x as above, base class Mult does the right thing
+  ParNonlinearForm::Mult(x, y);
 }
-
-#endif  // _GPU_
diff --git a/src/gradNonLinearForm.hpp b/src/gradNonLinearForm.hpp
index 6be197070..7ec050d3a 100644
--- a/src/gradNonLinearForm.hpp
+++ b/src/gradNonLinearForm.hpp
@@ -51,12 +51,6 @@ class GradNonLinearForm : public ParNonlinearForm {
   GradNonLinearForm(ParFiniteElementSpace *f, IntegrationRules *intRules, const int dim, const int num_equation);
 
   void Mult(const ParGridFunction *Up, Vector &y);
-
-#ifdef _GPU_
-  void Mult_gpu(const Vector &x, Vector &y);
-  static void IndexedAddToGlobalMemory(const Array<int> &vdofs, Vector &y, const Vector &el_y, const int &dof,
-                                       const int &elDof, const int &num_equation, const int &dim);
-#endif
 };
 
 #endif  // GRADNONLINEARFORM_HPP_