diff --git a/palace/fem/coefficient.hpp b/palace/fem/coefficient.hpp
index 3ff76e661..97bf711c4 100644
--- a/palace/fem/coefficient.hpp
+++ b/palace/fem/coefficient.hpp
@@ -508,14 +508,14 @@ class GradFluxCoefficient : public mfem::Coefficient
 private:
   const mfem::ParGridFunction &gf;
   const MaterialOperator &mat_op;
-  mfem::Vector grad, V;
+  mfem::Vector grad, tmp;
   int component;
 
 public:
   GradFluxCoefficient(const mfem::ParGridFunction &gf, const MaterialOperator &mat_op)
     : mfem::Coefficient(), gf(gf), mat_op(mat_op),
       grad(gf.ParFESpace()->GetParMesh()->SpaceDimension()),
-      V(gf.ParFESpace()->GetParMesh()->SpaceDimension()), component(-1)
+      tmp(gf.ParFESpace()->GetParMesh()->SpaceDimension()), component(-1)
   {
   }
 
@@ -532,9 +532,16 @@ class GradFluxCoefficient : public mfem::Coefficient
     MFEM_ASSERT(component >= 0 &&
                     component < gf.ParFESpace()->GetParMesh()->SpaceDimension(),
                 "Invalid component index, try calling SetComponent(int)!");
+    Eval(tmp, T, ip);
+    return tmp(component);
+  }
+
+  // VectorCoefficient style Eval method for use in elemental integral evaluation.
+  inline void Eval(mfem::Vector &V, mfem::ElementTransformation &T,
+                   const mfem::IntegrationPoint &ip)
+  {
     gf.GetGradient(T, grad);
     mat_op.GetPermittivityReal(T.Attribute).Mult(grad, V);
-    return V(component);
   }
 };
 
diff --git a/palace/linalg/errorestimator.cpp b/palace/linalg/errorestimator.cpp
index 34d414821..57e07370f 100644
--- a/palace/linalg/errorestimator.cpp
+++ b/palace/linalg/errorestimator.cpp
@@ -278,7 +278,8 @@ GradFluxErrorEstimator::GradFluxErrorEstimator(
     smooth_flux(h1_fespaces.GetFinestFESpace().GetTrueVSize()),
     flux_rhs(h1_fespaces.GetFinestFESpace().GetTrueVSize()),
     field_gf(&h1_fespaces.GetFinestFESpace()),
-    smooth_flux_gf(&h1_fespaces.GetFinestFESpace())
+    smooth_flux_gf{&h1_fespaces.GetFinestFESpace(), &h1_fespaces.GetFinestFESpace(),
+                   &h1_fespaces.GetFinestFESpace()}
 {
 }
 
@@ -310,39 +311,35 @@ ErrorIndicator GradFluxErrorEstimator::ComputeIndicators(const Vector &v) const
       h1_fespace.GetProlongationMatrix()->MultTranspose(rhs, flux_rhs);
       smooth_projector.Mult(flux_rhs, smooth_flux);
     }
-    smooth_flux_gf.SetFromTrueDofs(smooth_flux);
-    smooth_flux_gf.ExchangeFaceNbrData();
-    for (int e = 0; e < h1_fespace.GetNE(); e++)
+    smooth_flux_gf[c].SetFromTrueDofs(smooth_flux);
+    smooth_flux_gf[c].ExchangeFaceNbrData();
+  }
+
+  Vector coef_eval(3);
+  for (int e = 0; e < h1_fespace.GetNE(); e++)
+  {
+    auto &T = *h1_fespace.GetElementTransformation(e);
+    // integration order 2p + q
+    const auto &ir = mfem::IntRules.Get(T.GetGeometryType(),
+                                        2 * h1_fespace.GetFE(e)->GetOrder() + T.Order());
+    for (const auto &ip : ir)
     {
-      auto &T = *h1_fespace.GetElementTransformation(e);
-      // integration order 2p + q
-      const auto &ir = mfem::IntRules.Get(T.GetGeometryType(),
-                                          2 * h1_fespace.GetFE(e)->GetOrder() + T.Order());
-      for (const auto &ip : ir)
+      T.SetIntPoint(&ip);
+      coef.Eval(coef_eval, T, ip);
+      for (int c = 0; c < sdim; c++)
       {
-        T.SetIntPoint(&ip);
-        double smooth_val = smooth_flux_gf.GetValue(e, ip);
-        double coarse_val = coef.Eval(T, ip);
+        coef.SetComponent(c);
+        double smooth_val = smooth_flux_gf[c].GetValue(e, ip);
         const double w_i = ip.weight * T.Weight();
-        estimates[e] += w_i * std::pow(smooth_val - coarse_val, 2.0);
-        normalization += w_i * std::pow(coarse_val, 2.0);
+        estimates[e] += w_i * std::pow(smooth_val - coef_eval(c), 2.0);
+        normalization += w_i * std::pow(coef_eval(c), 2.0);
       }
     }
-    if constexpr (false)
-    {
-      // Debugging branch generates some intermediate fields for paraview.
-      mfem::ParaViewDataCollection paraview("debug_coeff" + std::to_string(c),
-                                            h1_fespace.GetParMesh());
-      paraview.RegisterCoeffField("Flux", &coef);
-      paraview.RegisterField("SmoothFlux", &smooth_flux_gf);
-      paraview.Save();
-    }
+    estimates[e] = std::sqrt(estimates[e]);
   }
-  linalg::Sqrt(estimates);
 
   Mpi::GlobalSum(1, &normalization, h1_fespace.GetComm());
   normalization = std::sqrt(normalization);
-
   if (normalization > 0)
   {
     estimates /= normalization;
diff --git a/palace/linalg/errorestimator.hpp b/palace/linalg/errorestimator.hpp
index 872ec9f58..c4f1b7c52 100644
--- a/palace/linalg/errorestimator.hpp
+++ b/palace/linalg/errorestimator.hpp
@@ -85,7 +85,8 @@ class GradFluxErrorEstimator
 
   FluxProjector<mfem::H1_FECollection> smooth_projector;
   mutable Vector smooth_flux, flux_rhs;
-  mutable mfem::ParGridFunction field_gf, smooth_flux_gf;
+  mutable mfem::ParGridFunction field_gf;
+  mutable std::array<mfem::ParGridFunction, 3> smooth_flux_gf;
 
 public:
   // Constructor for using geometric and p multigrid.