diff --git a/palace/drivers/basesolver.cpp b/palace/drivers/basesolver.cpp
index d171d99c0..6f039129c 100644
--- a/palace/drivers/basesolver.cpp
+++ b/palace/drivers/basesolver.cpp
@@ -568,7 +568,7 @@ void BaseSolver::PostprocessProbes(const PostOperator &postop, const std::string
   {
     return;
   }
-  const bool has_imaginary = postop.HasImaginary();
+  const bool has_imaginary = postop.HasImag();
   for (int f = 0; f < 2; f++)
   {
     // Probe data is ordered as [Fx1, Fy1, Fz1, Fx2, Fy2, Fz2, ...].
diff --git a/palace/fem/CMakeLists.txt b/palace/fem/CMakeLists.txt
index 9312e08f7..4bb308811 100644
--- a/palace/fem/CMakeLists.txt
+++ b/palace/fem/CMakeLists.txt
@@ -11,6 +11,7 @@ target_sources(${LIB_TARGET_NAME}
   ${CMAKE_CURRENT_SOURCE_DIR}/coefficient.cpp
   ${CMAKE_CURRENT_SOURCE_DIR}/errorindicator.cpp
   ${CMAKE_CURRENT_SOURCE_DIR}/fespace.cpp
+  ${CMAKE_CURRENT_SOURCE_DIR}/gridfunction.cpp
   ${CMAKE_CURRENT_SOURCE_DIR}/integrator.cpp
   ${CMAKE_CURRENT_SOURCE_DIR}/interpolator.cpp
   ${CMAKE_CURRENT_SOURCE_DIR}/lumpedelement.cpp
diff --git a/palace/fem/bilinearform.cpp b/palace/fem/bilinearform.cpp
index 1661cf9c9..b97f33d93 100644
--- a/palace/fem/bilinearform.cpp
+++ b/palace/fem/bilinearform.cpp
@@ -260,14 +260,20 @@ std::unique_ptr<ceed::Operator> DiscreteLinearOperator::PartialAssemble() const
   Vector test_multiplicity(test_fespace.GetVSize());
   test_multiplicity = 0.0;
   auto *h_mult = test_multiplicity.HostReadWrite();
-  mfem::Array<int> dofs;
-  for (int i = 0; i < test_fespace.GetMesh().GetNE(); i++)
+  PalacePragmaOmp(parallel)
   {
-    test_fespace.Get().GetElementVDofs(i, dofs);
-    for (int j = 0; j < dofs.Size(); j++)
+    mfem::Array<int> dofs;
+    mfem::DofTransformation dof_trans;
+    PalacePragmaOmp(for schedule(static))
+    for (int i = 0; i < test_fespace.GetMesh().GetNE(); i++)
     {
-      const int k = dofs[j];
-      h_mult[(k >= 0) ? k : -1 - k] += 1.0;
+      test_fespace.Get().GetElementVDofs(i, dofs, dof_trans);
+      for (int j = 0; j < dofs.Size(); j++)
+      {
+        const int k = dofs[j];
+        PalacePragmaOmp(atomic update)
+        h_mult[(k >= 0) ? k : -1 - k] += 1.0;
+      }
     }
   }
   test_multiplicity.UseDevice(true);
diff --git a/palace/fem/coefficient.hpp b/palace/fem/coefficient.hpp
index aebacb3ec..0d57a256a 100644
--- a/palace/fem/coefficient.hpp
+++ b/palace/fem/coefficient.hpp
@@ -9,6 +9,7 @@
 #include <utility>
 #include <vector>
 #include <mfem.hpp>
+#include "fem/gridfunction.hpp"
 #include "models/materialoperator.hpp"
 
 namespace palace
@@ -350,11 +351,11 @@ enum class EnergyDensityType
 // Returns the local energy density evaluated as 1/2 Dᴴ E or 1/2 Bᴴ H for real-valued
 // material coefficients. For internal boundary elements, the solution is taken on the side
 // of the element with the larger-valued material property (permittivity or permeability).
-template <EnergyDensityType Type, typename GridFunctionType>
+template <EnergyDensityType Type>
 class EnergyDensityCoefficient : public mfem::Coefficient, public BdrGridFunctionCoefficient
 {
 private:
-  const GridFunctionType &U;
+  const GridFunction &U;
   const MaterialOperator &mat_op;
   mfem::Vector V;
 
@@ -362,7 +363,7 @@ class EnergyDensityCoefficient : public mfem::Coefficient, public BdrGridFunctio
                                const mfem::IntegrationPoint &ip, int attr);
 
 public:
-  EnergyDensityCoefficient(const GridFunctionType &gf, const MaterialOperator &mat_op)
+  EnergyDensityCoefficient(const GridFunction &gf, const MaterialOperator &mat_op)
     : mfem::Coefficient(), BdrGridFunctionCoefficient(*gf.ParFESpace()->GetParMesh()),
       U(gf), mat_op(mat_op), V(mat_op.SpaceDimension())
   {
@@ -399,49 +400,33 @@ class EnergyDensityCoefficient : public mfem::Coefficient, public BdrGridFunctio
 };
 
 template <>
-inline double
-EnergyDensityCoefficient<EnergyDensityType::ELECTRIC, mfem::ParComplexGridFunction>::
-    GetLocalEnergyDensity(mfem::ElementTransformation &T, const mfem::IntegrationPoint &ip,
-                          int attr)
+inline double EnergyDensityCoefficient<EnergyDensityType::ELECTRIC>::GetLocalEnergyDensity(
+    mfem::ElementTransformation &T, const mfem::IntegrationPoint &ip, int attr)
 {
   // Only the real part of the permittivity contributes to the energy (imaginary part
   // cancels out in the inner product due to symmetry).
-  U.real().GetVectorValue(T, ip, V);
-  double res = mat_op.GetPermittivityReal(attr).InnerProduct(V, V);
-  U.imag().GetVectorValue(T, ip, V);
-  res += mat_op.GetPermittivityReal(attr).InnerProduct(V, V);
-  return 0.5 * res;
-}
-
-template <>
-inline double EnergyDensityCoefficient<EnergyDensityType::ELECTRIC, mfem::ParGridFunction>::
-    GetLocalEnergyDensity(mfem::ElementTransformation &T, const mfem::IntegrationPoint &ip,
-                          int attr)
-{
-  U.GetVectorValue(T, ip, V);
-  return 0.5 * mat_op.GetPermittivityReal(attr).InnerProduct(V, V);
-}
-
-template <>
-inline double
-EnergyDensityCoefficient<EnergyDensityType::MAGNETIC, mfem::ParComplexGridFunction>::
-    GetLocalEnergyDensity(mfem::ElementTransformation &T, const mfem::IntegrationPoint &ip,
-                          int attr)
-{
-  U.real().GetVectorValue(T, ip, V);
-  double res = mat_op.GetInvPermeability(attr).InnerProduct(V, V);
-  U.imag().GetVectorValue(T, ip, V);
-  res += mat_op.GetInvPermeability(attr).InnerProduct(V, V);
-  return 0.5 * res;
+  U.Real().GetVectorValue(T, ip, V);
+  double dot = mat_op.GetPermittivityReal(attr).InnerProduct(V, V);
+  if (U.HasImag())
+  {
+    U.Imag().GetVectorValue(T, ip, V);
+    dot += mat_op.GetPermittivityReal(attr).InnerProduct(V, V);
+  }
+  return 0.5 * dot;
 }
 
 template <>
-inline double EnergyDensityCoefficient<EnergyDensityType::MAGNETIC, mfem::ParGridFunction>::
-    GetLocalEnergyDensity(mfem::ElementTransformation &T, const mfem::IntegrationPoint &ip,
-                          int attr)
+inline double EnergyDensityCoefficient<EnergyDensityType::MAGNETIC>::GetLocalEnergyDensity(
+    mfem::ElementTransformation &T, const mfem::IntegrationPoint &ip, int attr)
 {
-  U.GetVectorValue(T, ip, V);
-  return 0.5 * mat_op.GetInvPermeability(attr).InnerProduct(V, V);
+  U.Real().GetVectorValue(T, ip, V);
+  double dot = mat_op.GetInvPermeability(attr).InnerProduct(V, V);
+  if (U.HasImag())
+  {
+    U.Imag().GetVectorValue(T, ip, V);
+    dot += mat_op.GetInvPermeability(attr).InnerProduct(V, V);
+  }
+  return 0.5 * dot;
 }
 
 // Returns the local field evaluated on a boundary element. For internal boundary elements,
diff --git a/palace/fem/gridfunction.cpp b/palace/fem/gridfunction.cpp
new file mode 100644
index 000000000..3bc89b4f4
--- /dev/null
+++ b/palace/fem/gridfunction.cpp
@@ -0,0 +1,60 @@
+// Copyright Amazon.com, Inc. or its affiliates. All Rights Reserved.
+// SPDX-License-Identifier: Apache-2.0
+
+#include "gridfunction.hpp"
+
+#include "fem/fespace.hpp"
+
+namespace palace
+{
+
+GridFunction::GridFunction(mfem::ParFiniteElementSpace &fespace, bool complex)
+  : gfr(&fespace)
+{
+  if (complex)
+  {
+    gfi.SetSpace(&fespace);
+  }
+}
+
+GridFunction::GridFunction(FiniteElementSpace &fespace, bool complex)
+  : GridFunction(fespace.Get(), complex)
+{
+}
+
+GridFunction &GridFunction::operator=(std::complex<double> s)
+{
+  Real() = s.real();
+  if (HasImag())
+  {
+    Imag() = s.imag();
+  }
+  else
+  {
+    MFEM_ASSERT(
+        s.imag() == 0.0,
+        "Cannot assign complex scalar to a non-complex-valued GridFunction object!");
+  }
+  return *this;
+}
+
+GridFunction &GridFunction::operator*=(double s)
+{
+  Real() *= s;
+  if (HasImag())
+  {
+    Imag() *= s;
+  }
+  return *this;
+}
+
+void GridFunction::Update()
+{
+  Real().Update();
+  if (HasImag())
+  {
+    Imag().Update();
+  }
+}
+
+}  // namespace palace
diff --git a/palace/fem/gridfunction.hpp b/palace/fem/gridfunction.hpp
new file mode 100644
index 000000000..01db3aa23
--- /dev/null
+++ b/palace/fem/gridfunction.hpp
@@ -0,0 +1,73 @@
+// Copyright Amazon.com, Inc. or its affiliates. All Rights Reserved.
+// SPDX-License-Identifier: Apache-2.0
+
+#ifndef PALACE_FEM_GRIDFUNCTION_HPP
+#define PALACE_FEM_GRIDFUNCTION_HPP
+
+#include <mfem.hpp>
+
+namespace palace
+{
+
+class FiniteElementSpace;
+
+//
+// A real- or complex-valued grid function represented as two real grid functions, one for
+// each component. This unifies mfem::ParGridFunction and mfem::ParComplexGridFunction, and
+// replaces the latter due to some issues observed with memory aliasing on GPUs.
+//
+class GridFunction
+{
+private:
+  mfem::ParGridFunction gfr, gfi;
+
+public:
+  GridFunction(mfem::ParFiniteElementSpace &fespace, bool complex = false);
+  GridFunction(FiniteElementSpace &fespace, bool complex = false);
+
+  // Get access to the real and imaginary grid function parts.
+  const mfem::ParGridFunction &Real() const { return gfr; }
+  mfem::ParGridFunction &Real() { return gfr; }
+  const mfem::ParGridFunction &Imag() const
+  {
+    MFEM_ASSERT(HasImag(),
+                "Invalid access of imaginary part of a real-valued GridFunction object!");
+    return gfi;
+  }
+  mfem::ParGridFunction &Imag()
+  {
+    MFEM_ASSERT(HasImag(),
+                "Invalid access of imaginary part of a real-valued GridFunction object!");
+    return gfi;
+  }
+
+  // Check if the grid function is suited for storing complex-valued fields.
+  bool HasImag() const { return (gfi.ParFESpace() != nullptr); }
+
+  // Get access to the underlying finite element space (match MFEM interface).
+  mfem::FiniteElementSpace *FESpace() { return gfr.FESpace(); }
+  const mfem::FiniteElementSpace *FESpace() const { return gfr.FESpace(); }
+  mfem::ParFiniteElementSpace *ParFESpace() { return gfr.ParFESpace(); }
+  const mfem::ParFiniteElementSpace *ParFESpace() const { return gfr.ParFESpace(); }
+
+  // Set all entries equal to s.
+  GridFunction &operator=(std::complex<double> s);
+  GridFunction &operator=(double s)
+  {
+    *this = std::complex<double>(s, 0.0);
+    return *this;
+  }
+
+  // Scale all entries by s.
+  GridFunction &operator*=(double s);
+
+  // Transform for space update (for example on mesh change).
+  void Update();
+
+  // Get the associated MPI communicator.
+  MPI_Comm GetComm() const { return ParFESpace()->GetComm(); }
+};
+
+}  // namespace palace
+
+#endif  // PALACE_FEM_GRIDFUNCTION_HPP
diff --git a/palace/fem/interpolator.cpp b/palace/fem/interpolator.cpp
index 5458ba730..9efd7aa85 100644
--- a/palace/fem/interpolator.cpp
+++ b/palace/fem/interpolator.cpp
@@ -4,6 +4,7 @@
 #include "interpolator.hpp"
 
 #include <algorithm>
+#include "fem/gridfunction.hpp"
 #include "utils/communication.hpp"
 #include "utils/iodata.hpp"
 
@@ -88,13 +89,12 @@ std::vector<double> InterpolationOperator::ProbeField(const mfem::ParGridFunctio
 #endif
 }
 
-std::vector<std::complex<double>>
-InterpolationOperator::ProbeField(const mfem::ParComplexGridFunction &U, bool has_imaginary)
+std::vector<std::complex<double>> InterpolationOperator::ProbeField(const GridFunction &U)
 {
-  std::vector<double> vr = ProbeField(U.real());
-  if (has_imaginary)
+  std::vector<double> vr = ProbeField(U.Real());
+  if (U.HasImag())
   {
-    std::vector<double> vi = ProbeField(U.imag());
+    std::vector<double> vi = ProbeField(U.Imag());
     std::vector<std::complex<double>> vals(vr.size());
     std::transform(vr.begin(), vr.end(), vi.begin(), vals.begin(),
                    [](double xr, double xi) { return std::complex<double>(xr, xi); });
diff --git a/palace/fem/interpolator.hpp b/palace/fem/interpolator.hpp
index 9b5ec9685..460caa2e9 100644
--- a/palace/fem/interpolator.hpp
+++ b/palace/fem/interpolator.hpp
@@ -11,6 +11,7 @@
 namespace palace
 {
 
+class GridFunction;
 class IoData;
 
 //
@@ -24,15 +25,14 @@ class InterpolationOperator
 #endif
   std::vector<int> op_idx;
 
+  std::vector<double> ProbeField(const mfem::ParGridFunction &U);
+
 public:
   InterpolationOperator(const IoData &iodata, mfem::ParMesh &mesh);
 
   const auto &GetProbes() const { return op_idx; }
 
-  std::vector<double> ProbeField(const mfem::ParGridFunction &U);
-
-  std::vector<std::complex<double>> ProbeField(const mfem::ParComplexGridFunction &U,
-                                               bool has_imaginary);
+  std::vector<std::complex<double>> ProbeField(const GridFunction &U);
 };
 
 }  // namespace palace
diff --git a/palace/models/domainpostoperator.cpp b/palace/models/domainpostoperator.cpp
index 11674277d..4a72c7efd 100644
--- a/palace/models/domainpostoperator.cpp
+++ b/palace/models/domainpostoperator.cpp
@@ -6,6 +6,7 @@
 #include <mfem.hpp>
 #include "fem/bilinearform.hpp"
 #include "fem/fespace.hpp"
+#include "fem/gridfunction.hpp"
 #include "fem/integrator.hpp"
 #include "models/materialoperator.hpp"
 #include "utils/communication.hpp"
@@ -74,70 +75,46 @@ DomainPostOperator::DomainPostOperator(const IoData &iodata, const MaterialOpera
   }
 }
 
-double
-DomainPostOperator::GetElectricFieldEnergy(const mfem::ParComplexGridFunction &E) const
+double DomainPostOperator::GetElectricFieldEnergy(const GridFunction &E) const
 {
   if (M_ND)
   {
-    M_ND->Mult(E.real(), D);
-    double res = linalg::LocalDot(E.real(), D);
-    M_ND->Mult(E.imag(), D);
-    res += linalg::LocalDot(E.imag(), D);
-    Mpi::GlobalSum(1, &res, E.ParFESpace()->GetComm());
-    return 0.5 * res;
-  }
-  MFEM_ABORT(
-      "Domain postprocessing is not configured for electric field energy calculation!");
-  return 0.0;
-}
-
-double DomainPostOperator::GetElectricFieldEnergy(const mfem::ParGridFunction &E) const
-{
-  if (M_ND)
-  {
-    M_ND->Mult(E, D);
-    double res = linalg::LocalDot(E, D);
-    Mpi::GlobalSum(1, &res, E.ParFESpace()->GetComm());
-    return 0.5 * res;
+    M_ND->Mult(E.Real(), D);
+    double dot = linalg::LocalDot(E.Real(), D);
+    if (E.HasImag())
+    {
+      M_ND->Mult(E.Imag(), D);
+      dot += linalg::LocalDot(E.Imag(), D);
+    }
+    Mpi::GlobalSum(1, &dot, E.GetComm());
+    return 0.5 * dot;
   }
   MFEM_ABORT(
       "Domain postprocessing is not configured for electric field energy calculation!");
   return 0.0;
 }
 
-double
-DomainPostOperator::GetMagneticFieldEnergy(const mfem::ParComplexGridFunction &B) const
+double DomainPostOperator::GetMagneticFieldEnergy(const GridFunction &B) const
 {
   if (M_RT)
   {
-    M_RT->Mult(B.real(), H);
-    double res = linalg::LocalDot(B.real(), H);
-    M_RT->Mult(B.imag(), H);
-    res += linalg::LocalDot(B.imag(), H);
-    Mpi::GlobalSum(1, &res, B.ParFESpace()->GetComm());
-    return 0.5 * res;
-  }
-  MFEM_ABORT(
-      "Domain postprocessing is not configured for magnetic field energy calculation!");
-  return 0.0;
-}
-
-double DomainPostOperator::GetMagneticFieldEnergy(const mfem::ParGridFunction &B) const
-{
-  if (M_RT)
-  {
-    M_RT->Mult(B, H);
-    double res = linalg::LocalDot(B, H);
-    Mpi::GlobalSum(1, &res, B.ParFESpace()->GetComm());
-    return 0.5 * res;
+    M_RT->Mult(B.Real(), H);
+    double dot = linalg::LocalDot(B.Real(), H);
+    if (B.HasImag())
+    {
+      M_RT->Mult(B.Imag(), H);
+      dot += linalg::LocalDot(B.Imag(), H);
+    }
+    Mpi::GlobalSum(1, &dot, B.GetComm());
+    return 0.5 * dot;
   }
   MFEM_ABORT(
       "Domain postprocessing is not configured for magnetic field energy calculation!");
   return 0.0;
 }
 
-double DomainPostOperator::GetDomainElectricFieldEnergy(
-    int idx, const mfem::ParComplexGridFunction &E) const
+double DomainPostOperator::GetDomainElectricFieldEnergy(int idx,
+                                                        const GridFunction &E) const
 {
   // Compute the electric field energy integral for only a portion of the domain.
   auto it = M_i.find(idx);
@@ -147,33 +124,19 @@ double DomainPostOperator::GetDomainElectricFieldEnergy(
   {
     return 0.0;
   }
-  it->second.first->Mult(E.real(), D);
-  double res = linalg::LocalDot(E.real(), D);
-  it->second.first->Mult(E.imag(), D);
-  res += linalg::LocalDot(E.imag(), D);
-  Mpi::GlobalSum(1, &res, E.ParFESpace()->GetComm());
-  return 0.5 * res;
-}
-
-double
-DomainPostOperator::GetDomainElectricFieldEnergy(int idx,
-                                                 const mfem::ParGridFunction &E) const
-{
-  auto it = M_i.find(idx);
-  MFEM_VERIFY(it != M_i.end() && it->second.first,
-              "Invalid domain index when postprocessing domain electric field energy!");
-  if (!it->second.first)
+  it->second.first->Mult(E.Real(), D);
+  double dot = linalg::LocalDot(E.Real(), D);
+  if (E.HasImag())
   {
-    return 0.0;
+    it->second.first->Mult(E.Imag(), D);
+    dot += linalg::LocalDot(E.Imag(), D);
   }
-  it->second.first->Mult(E, D);
-  double res = linalg::LocalDot(E, D);
-  Mpi::GlobalSum(1, &res, E.ParFESpace()->GetComm());
-  return 0.5 * res;
+  Mpi::GlobalSum(1, &dot, E.GetComm());
+  return 0.5 * dot;
 }
 
-double DomainPostOperator::GetDomainMagneticFieldEnergy(
-    int idx, const mfem::ParComplexGridFunction &B) const
+double DomainPostOperator::GetDomainMagneticFieldEnergy(int idx,
+                                                        const GridFunction &B) const
 {
   // Compute the magnetic field energy integral for only a portion of the domain.
   auto it = M_i.find(idx);
@@ -183,29 +146,15 @@ double DomainPostOperator::GetDomainMagneticFieldEnergy(
   {
     return 0.0;
   }
-  it->second.second->Mult(B.real(), H);
-  double res = linalg::LocalDot(B.real(), H);
-  it->second.second->Mult(B.imag(), H);
-  res += linalg::LocalDot(B.imag(), H);
-  Mpi::GlobalSum(1, &res, B.ParFESpace()->GetComm());
-  return 0.5 * res;
-}
-
-double
-DomainPostOperator::GetDomainMagneticFieldEnergy(int idx,
-                                                 const mfem::ParGridFunction &B) const
-{
-  auto it = M_i.find(idx);
-  MFEM_VERIFY(it != M_i.end() && it->second.second,
-              "Invalid domain index when postprocessing domain magnetic field energy!");
-  if (!it->second.second)
+  it->second.second->Mult(B.Real(), H);
+  double dot = linalg::LocalDot(B.Real(), H);
+  if (B.HasImag())
   {
-    return 0.0;
+    it->second.second->Mult(B.Imag(), H);
+    dot += linalg::LocalDot(B.Imag(), H);
   }
-  it->second.second->Mult(B, H);
-  double res = linalg::LocalDot(B, H);
-  Mpi::GlobalSum(1, &res, B.ParFESpace()->GetComm());
-  return 0.5 * res;
+  Mpi::GlobalSum(1, &dot, B.GetComm());
+  return 0.5 * dot;
 }
 
 }  // namespace palace
diff --git a/palace/models/domainpostoperator.hpp b/palace/models/domainpostoperator.hpp
index dc80c7ed9..980700ce9 100644
--- a/palace/models/domainpostoperator.hpp
+++ b/palace/models/domainpostoperator.hpp
@@ -10,17 +10,10 @@
 #include "linalg/operator.hpp"
 #include "linalg/vector.hpp"
 
-namespace mfem
-{
-
-class ParComplexGridFunction;
-class ParGridFunction;
-
-}  // namespace mfem
-
 namespace palace
 {
 
+class GridFunction;
 class FiniteElementSpace;
 class IoData;
 class MaterialOperator;
@@ -48,17 +41,13 @@ class DomainPostOperator
 
   // Get volume integrals computing the electric or magnetic field energy in the entire
   // domain.
-  double GetElectricFieldEnergy(const mfem::ParComplexGridFunction &E) const;
-  double GetElectricFieldEnergy(const mfem::ParGridFunction &E) const;
-  double GetMagneticFieldEnergy(const mfem::ParComplexGridFunction &B) const;
-  double GetMagneticFieldEnergy(const mfem::ParGridFunction &B) const;
+  double GetElectricFieldEnergy(const GridFunction &E) const;
+  double GetMagneticFieldEnergy(const GridFunction &B) const;
 
   // Get volume integrals for the electric or magnetic field energy in a portion of the
   // domain.
-  double GetDomainElectricFieldEnergy(int idx, const mfem::ParComplexGridFunction &E) const;
-  double GetDomainElectricFieldEnergy(int idx, const mfem::ParGridFunction &E) const;
-  double GetDomainMagneticFieldEnergy(int idx, const mfem::ParComplexGridFunction &E) const;
-  double GetDomainMagneticFieldEnergy(int idx, const mfem::ParGridFunction &E) const;
+  double GetDomainElectricFieldEnergy(int idx, const GridFunction &E) const;
+  double GetDomainMagneticFieldEnergy(int idx, const GridFunction &E) const;
 };
 
 }  // namespace palace
diff --git a/palace/models/lumpedportoperator.cpp b/palace/models/lumpedportoperator.cpp
index 78b3bf18d..59917b50e 100644
--- a/palace/models/lumpedportoperator.cpp
+++ b/palace/models/lumpedportoperator.cpp
@@ -4,6 +4,7 @@
 #include "lumpedportoperator.hpp"
 
 #include "fem/coefficient.hpp"
+#include "fem/gridfunction.hpp"
 #include "fem/integrator.hpp"
 #include "models/materialoperator.hpp"
 #include "utils/communication.hpp"
@@ -215,101 +216,85 @@ void LumpedPortData::InitializeLinearForms(mfem::ParFiniteElementSpace &nd_fespa
   }
 }
 
-std::complex<double> LumpedPortData::GetSParameter(mfem::ParComplexGridFunction &E) const
-{
-  // Compute port S-parameter, or the projection of the field onto the port mode.
-  InitializeLinearForms(*E.ParFESpace());
-  std::complex<double> dot((*s) * E.real(), (*s) * E.imag());
-  Mpi::GlobalSum(1, &dot, E.ParFESpace()->GetComm());
-  return dot;
-}
-
-double LumpedPortData::GetPower(mfem::ParGridFunction &E, mfem::ParGridFunction &B) const
+std::complex<double> LumpedPortData::GetPower(GridFunction &E, GridFunction &B) const
 {
   // Compute port power, (E x H) ⋅ n = E ⋅ (-n x H), integrated over the port surface
   // using the computed E and H = μ⁻¹ B fields. The linear form is reconstructed from
   // scratch each time due to changing H. The BdrCurrentVectorCoefficient computes -n x H,
   // where n is an outward normal.
+  MFEM_VERIFY((E.HasImag() && B.HasImag()) || (!E.HasImag() && !B.HasImag()),
+              "Mismatch between real- and complex-valued E and B fields in port power "
+              "calculation!");
+  const bool has_imag = E.HasImag();
   auto &nd_fespace = *E.ParFESpace();
   const auto &mesh = *nd_fespace.GetParMesh();
-  SumVectorCoefficient fb(mesh.SpaceDimension());
-  mfem::Array<int> attr_list;
-  for (const auto &elem : elems)
-  {
-    fb.AddCoefficient(
-        std::make_unique<RestrictedVectorCoefficient<BdrCurrentVectorCoefficient>>(
-            elem->GetAttrList(), B, mat_op));
-    attr_list.Append(elem->GetAttrList());
-  }
-  int bdr_attr_max = mesh.bdr_attributes.Size() ? mesh.bdr_attributes.Max() : 0;
-  mfem::Array<int> attr_marker = mesh::AttrToMarker(bdr_attr_max, attr_list);
-  mfem::LinearForm p(&nd_fespace);
-  p.AddBoundaryIntegrator(new VectorFEBoundaryLFIntegrator(fb), attr_marker);
-  p.UseFastAssembly(false);
-  p.UseDevice(false);
-  p.Assemble();
-  p.UseDevice(true);
-  double dot = p * E;
-  Mpi::GlobalSum(1, &dot, E.ParFESpace()->GetComm());
-  return dot;
-}
-
-std::complex<double> LumpedPortData::GetPower(mfem::ParComplexGridFunction &E,
-                                              mfem::ParComplexGridFunction &B) const
-{
-  // Compute port power, (E x H⋆) ⋅ n = E ⋅ (-n x H⋆), integrated over the port surface
-  // using the computed E and H = μ⁻¹ B fields. The linear form is reconstructed from
-  // scratch each time due to changing H. The BdrCurrentVectorCoefficient computes -n x H,
-  // where n is an outward normal.
-  auto &nd_fespace = *E.ParFESpace();
-  const auto &mesh = *nd_fespace.GetParMesh();
-  SumVectorCoefficient fbr(mesh.SpaceDimension());
-  SumVectorCoefficient fbi(mesh.SpaceDimension());
+  SumVectorCoefficient fbr(mesh.SpaceDimension()), fbi(mesh.SpaceDimension());
   mfem::Array<int> attr_list;
   for (const auto &elem : elems)
   {
     fbr.AddCoefficient(
         std::make_unique<RestrictedVectorCoefficient<BdrCurrentVectorCoefficient>>(
-            elem->GetAttrList(), B.real(), mat_op));
-    fbi.AddCoefficient(
-        std::make_unique<RestrictedVectorCoefficient<BdrCurrentVectorCoefficient>>(
-            elem->GetAttrList(), B.imag(), mat_op));
+            elem->GetAttrList(), B.Real(), mat_op));
+    if (has_imag)
+    {
+      fbi.AddCoefficient(
+          std::make_unique<RestrictedVectorCoefficient<BdrCurrentVectorCoefficient>>(
+              elem->GetAttrList(), B.Imag(), mat_op));
+    }
     attr_list.Append(elem->GetAttrList());
   }
   int bdr_attr_max = mesh.bdr_attributes.Size() ? mesh.bdr_attributes.Max() : 0;
   mfem::Array<int> attr_marker = mesh::AttrToMarker(bdr_attr_max, attr_list);
-  mfem::LinearForm pr(&nd_fespace), pi(&nd_fespace);
+  mfem::LinearForm pr(&nd_fespace);
   pr.AddBoundaryIntegrator(new VectorFEBoundaryLFIntegrator(fbr), attr_marker);
-  pi.AddBoundaryIntegrator(new VectorFEBoundaryLFIntegrator(fbi), attr_marker);
   pr.UseFastAssembly(false);
-  pi.UseFastAssembly(false);
   pr.UseDevice(false);
-  pi.UseDevice(false);
   pr.Assemble();
-  pi.Assemble();
   pr.UseDevice(true);
-  pi.UseDevice(true);
-  std::complex<double> dot((pr * E.real()) + (pi * E.imag()),
-                           (pr * E.imag()) - (pi * E.real()));
-  Mpi::GlobalSum(1, &dot, E.ParFESpace()->GetComm());
-  return dot;
+  if (has_imag)
+  {
+    mfem::LinearForm pi(&nd_fespace);
+    pi.AddBoundaryIntegrator(new VectorFEBoundaryLFIntegrator(fbi), attr_marker);
+    pi.UseFastAssembly(false);
+    pi.UseDevice(false);
+    pi.Assemble();
+    pi.UseDevice(true);
+    std::complex<double> dot((pr * E.Real()) + (pi * E.Imag()),
+                             (pr * E.Imag()) - (pi * E.Real()));
+    Mpi::GlobalSum(1, &dot, E.ParFESpace()->GetComm());
+    return dot;
+  }
+  else
+  {
+    double dot = pr * E.Real();
+    Mpi::GlobalSum(1, &dot, E.ParFESpace()->GetComm());
+    return dot;
+  }
 }
 
-double LumpedPortData::GetVoltage(mfem::ParGridFunction &E) const
+std::complex<double> LumpedPortData::GetSParameter(GridFunction &E) const
 {
-  // Compute the average voltage across the port.
+  // Compute port S-parameter, or the projection of the field onto the port mode.
   InitializeLinearForms(*E.ParFESpace());
-  double dot = (*v) * E;
-  Mpi::GlobalSum(1, &dot, E.ParFESpace()->GetComm());
+  std::complex<double> dot((*s) * E.Real(), 0.0);
+  if (E.HasImag())
+  {
+    dot.imag((*s) * E.Imag());
+  }
+  Mpi::GlobalSum(1, &dot, E.GetComm());
   return dot;
 }
 
-std::complex<double> LumpedPortData::GetVoltage(mfem::ParComplexGridFunction &E) const
+std::complex<double> LumpedPortData::GetVoltage(GridFunction &E) const
 {
   // Compute the average voltage across the port.
   InitializeLinearForms(*E.ParFESpace());
-  std::complex<double> dot((*v) * E.real(), (*v) * E.imag());
-  Mpi::GlobalSum(1, &dot, E.ParFESpace()->GetComm());
+  std::complex<double> dot((*v) * E.Real(), 0.0);
+  if (E.HasImag())
+  {
+    dot.imag((*v) * E.Imag());
+  }
+  Mpi::GlobalSum(1, &dot, E.GetComm());
   return dot;
 }
 
diff --git a/palace/models/lumpedportoperator.hpp b/palace/models/lumpedportoperator.hpp
index a230d3f10..c9567b3c0 100644
--- a/palace/models/lumpedportoperator.hpp
+++ b/palace/models/lumpedportoperator.hpp
@@ -14,6 +14,7 @@
 namespace palace
 {
 
+class GridFunction;
 class IoData;
 class MaterialOperator;
 class MaterialPropertyCoefficient;
@@ -63,12 +64,9 @@ class LumpedPortData
   double GetExcitationPower() const;
   double GetExcitationVoltage() const;
 
-  std::complex<double> GetSParameter(mfem::ParComplexGridFunction &E) const;
-  std::complex<double> GetPower(mfem::ParComplexGridFunction &E,
-                                mfem::ParComplexGridFunction &B) const;
-  double GetPower(mfem::ParGridFunction &E, mfem::ParGridFunction &B) const;
-  std::complex<double> GetVoltage(mfem::ParComplexGridFunction &E) const;
-  double GetVoltage(mfem::ParGridFunction &E) const;
+  std::complex<double> GetPower(GridFunction &E, GridFunction &B) const;
+  std::complex<double> GetSParameter(GridFunction &E) const;
+  std::complex<double> GetVoltage(GridFunction &E) const;
 };
 
 //
diff --git a/palace/models/postoperator.cpp b/palace/models/postoperator.cpp
index a1d3ce918..745d5051d 100644
--- a/palace/models/postoperator.cpp
+++ b/palace/models/postoperator.cpp
@@ -43,39 +43,37 @@ PostOperator::PostOperator(const IoData &iodata, SpaceOperator &spaceop,
     surf_post_op(iodata, spaceop.GetMaterialOp(), spaceop.GetH1Space()),
     dom_post_op(iodata, spaceop.GetMaterialOp(), &spaceop.GetNDSpace(),
                 &spaceop.GetRTSpace()),
-    has_imaginary(iodata.problem.type != config::ProblemData::Type::TRANSIENT),
-    E(&spaceop.GetNDSpace().Get()), B(&spaceop.GetRTSpace().Get()), V(std::nullopt),
-    A(std::nullopt), lumped_port_init(false), wave_port_init(false),
+    E(std::in_place, spaceop.GetNDSpace(),
+      iodata.problem.type != config::ProblemData::Type::TRANSIENT),
+    B(std::in_place, spaceop.GetRTSpace(),
+      iodata.problem.type != config::ProblemData::Type::TRANSIENT),
+    V(std::nullopt), A(std::nullopt), lumped_port_init(false), wave_port_init(false),
     paraview(CreateParaviewPath(iodata, name), &spaceop.GetNDSpace().GetParMesh()),
     paraview_bdr(CreateParaviewPath(iodata, name) + "_boundary",
                  &spaceop.GetNDSpace().GetParMesh()),
     interp_op(iodata, spaceop.GetNDSpace().GetParMesh())
 {
-  Esr = std::make_unique<BdrFieldVectorCoefficient>(E->real(), mat_op);
-  Bsr = std::make_unique<BdrFieldVectorCoefficient>(B->real(), mat_op);
-  Jsr = std::make_unique<BdrCurrentVectorCoefficient>(B->real(), mat_op);
-  Qsr = std::make_unique<BdrChargeCoefficient>(E->real(), mat_op);
-  if (has_imaginary)
+  Esr = std::make_unique<BdrFieldVectorCoefficient>(E->Real(), mat_op);
+  Bsr = std::make_unique<BdrFieldVectorCoefficient>(B->Real(), mat_op);
+  Jsr = std::make_unique<BdrCurrentVectorCoefficient>(B->Real(), mat_op);
+  Qsr = std::make_unique<BdrChargeCoefficient>(E->Real(), mat_op);
+  if (HasImag())
   {
-    Esi = std::make_unique<BdrFieldVectorCoefficient>(E->imag(), mat_op);
-    Bsi = std::make_unique<BdrFieldVectorCoefficient>(B->imag(), mat_op);
-    Jsi = std::make_unique<BdrCurrentVectorCoefficient>(B->imag(), mat_op);
-    Qsi = std::make_unique<BdrChargeCoefficient>(E->imag(), mat_op);
-    Ue = std::make_unique<EnergyDensityCoefficient<EnergyDensityType::ELECTRIC,
-                                                   mfem::ParComplexGridFunction>>(*E,
-                                                                                  mat_op);
-    Um = std::make_unique<EnergyDensityCoefficient<EnergyDensityType::MAGNETIC,
-                                                   mfem::ParComplexGridFunction>>(*B,
-                                                                                  mat_op);
+    Esi = std::make_unique<BdrFieldVectorCoefficient>(E->Imag(), mat_op);
+    Bsi = std::make_unique<BdrFieldVectorCoefficient>(B->Imag(), mat_op);
+    Jsi = std::make_unique<BdrCurrentVectorCoefficient>(B->Imag(), mat_op);
+    Qsi = std::make_unique<BdrChargeCoefficient>(E->Imag(), mat_op);
+    Ue =
+        std::make_unique<EnergyDensityCoefficient<EnergyDensityType::ELECTRIC>>(*E, mat_op);
+    Um =
+        std::make_unique<EnergyDensityCoefficient<EnergyDensityType::MAGNETIC>>(*B, mat_op);
   }
   else
   {
-    Ue = std::make_unique<
-        EnergyDensityCoefficient<EnergyDensityType::ELECTRIC, mfem::ParGridFunction>>(
-        E->real(), mat_op);
-    Um = std::make_unique<
-        EnergyDensityCoefficient<EnergyDensityType::MAGNETIC, mfem::ParGridFunction>>(
-        B->real(), mat_op);
+    Ue =
+        std::make_unique<EnergyDensityCoefficient<EnergyDensityType::ELECTRIC>>(*E, mat_op);
+    Um =
+        std::make_unique<EnergyDensityCoefficient<EnergyDensityType::MAGNETIC>>(*B, mat_op);
   }
 
   // Add wave port boundary mode postprocessing when available.
@@ -95,9 +93,8 @@ PostOperator::PostOperator(const IoData &iodata, LaplaceOperator &laplaceop,
   : mat_op(laplaceop.GetMaterialOp()),
     surf_post_op(iodata, laplaceop.GetMaterialOp(), laplaceop.GetH1Space()),
     dom_post_op(iodata, laplaceop.GetMaterialOp(), &laplaceop.GetNDSpace(), nullptr),
-    has_imaginary(false), E(&laplaceop.GetNDSpace().Get()), B(std::nullopt),
-    V(&laplaceop.GetH1Space().Get()), A(std::nullopt), lumped_port_init(false),
-    wave_port_init(false),
+    E(std::in_place, laplaceop.GetNDSpace()), B(std::nullopt), V(laplaceop.GetH1Space()),
+    A(std::nullopt), lumped_port_init(false), wave_port_init(false),
     paraview(CreateParaviewPath(iodata, name), &laplaceop.GetNDSpace().GetParMesh()),
     paraview_bdr(CreateParaviewPath(iodata, name) + "_boundary",
                  &laplaceop.GetNDSpace().GetParMesh()),
@@ -106,12 +103,10 @@ PostOperator::PostOperator(const IoData &iodata, LaplaceOperator &laplaceop,
   // Note: When using this constructor, you should not use any of the magnetic field related
   // postprocessing functions (magnetic field energy, inductor energy, surface currents,
   // etc.), since only V and E fields are supplied.
-  Esr = std::make_unique<BdrFieldVectorCoefficient>(E->real(), mat_op);
-  Vs = std::make_unique<BdrFieldCoefficient>(*V, mat_op);
-  Ue = std::make_unique<
-      EnergyDensityCoefficient<EnergyDensityType::ELECTRIC, mfem::ParGridFunction>>(
-      E->real(), mat_op);
-  Qsr = std::make_unique<BdrChargeCoefficient>(E->real(), mat_op);
+  Esr = std::make_unique<BdrFieldVectorCoefficient>(E->Real(), mat_op);
+  Vs = std::make_unique<BdrFieldCoefficient>(V->Real(), mat_op);
+  Ue = std::make_unique<EnergyDensityCoefficient<EnergyDensityType::ELECTRIC>>(*E, mat_op);
+  Qsr = std::make_unique<BdrChargeCoefficient>(E->Real(), mat_op);
 
   // Initialize data collection objects.
   InitializeDataCollection(iodata);
@@ -122,9 +117,8 @@ PostOperator::PostOperator(const IoData &iodata, CurlCurlOperator &curlcurlop,
   : mat_op(curlcurlop.GetMaterialOp()),
     surf_post_op(iodata, curlcurlop.GetMaterialOp(), curlcurlop.GetH1Space()),
     dom_post_op(iodata, curlcurlop.GetMaterialOp(), nullptr, &curlcurlop.GetRTSpace()),
-    has_imaginary(false), E(std::nullopt), B(&curlcurlop.GetRTSpace().Get()),
-    V(std::nullopt), A(&curlcurlop.GetNDSpace().Get()), lumped_port_init(false),
-    wave_port_init(false),
+    E(std::nullopt), B(std::in_place, curlcurlop.GetRTSpace()), V(std::nullopt),
+    A(curlcurlop.GetNDSpace()), lumped_port_init(false), wave_port_init(false),
     paraview(CreateParaviewPath(iodata, name), &curlcurlop.GetNDSpace().GetParMesh()),
     paraview_bdr(CreateParaviewPath(iodata, name) + "_boundary",
                  &curlcurlop.GetNDSpace().GetParMesh()),
@@ -133,12 +127,10 @@ PostOperator::PostOperator(const IoData &iodata, CurlCurlOperator &curlcurlop,
   // Note: When using this constructor, you should not use any of the electric field related
   // postprocessing functions (electric field energy, capacitor energy, surface charge,
   // etc.), since only the B field is supplied.
-  Bsr = std::make_unique<BdrFieldVectorCoefficient>(B->real(), mat_op);
-  As = std::make_unique<BdrFieldVectorCoefficient>(*A, mat_op);
-  Um = std::make_unique<
-      EnergyDensityCoefficient<EnergyDensityType::MAGNETIC, mfem::ParGridFunction>>(
-      B->real(), mat_op);
-  Jsr = std::make_unique<BdrCurrentVectorCoefficient>(B->real(), mat_op);
+  Bsr = std::make_unique<BdrFieldVectorCoefficient>(B->Real(), mat_op);
+  As = std::make_unique<BdrFieldVectorCoefficient>(A->Real(), mat_op);
+  Um = std::make_unique<EnergyDensityCoefficient<EnergyDensityType::MAGNETIC>>(*B, mat_op);
+  Jsr = std::make_unique<BdrCurrentVectorCoefficient>(B->Real(), mat_op);
 
   // Initialize data collection objects.
   InitializeDataCollection(iodata);
@@ -155,8 +147,8 @@ void PostOperator::InitializeDataCollection(const IoData &iodata)
   const int compress = 0;
 #endif
   const bool use_ho = true;
-  const int refine_ho =
-      (E) ? E->ParFESpace()->GetMaxElementOrder() : B->ParFESpace()->GetMaxElementOrder();
+  const int refine_ho = HasE() ? E->ParFESpace()->GetMaxElementOrder()
+                               : B->ParFESpace()->GetMaxElementOrder();
   mesh_Lc0 = iodata.GetLengthScale();
 
   // Output mesh coordinate units same as input.
@@ -180,42 +172,42 @@ void PostOperator::InitializeDataCollection(const IoData &iodata)
   // permeability.
   if (E)
   {
-    if (has_imaginary)
+    if (HasImag())
     {
-      paraview.RegisterField("E_real", &E->real());
-      paraview.RegisterField("E_imag", &E->imag());
+      paraview.RegisterField("E_real", &E->Real());
+      paraview.RegisterField("E_imag", &E->Imag());
       paraview_bdr.RegisterVCoeffField("E_real", Esr.get());
       paraview_bdr.RegisterVCoeffField("E_imag", Esi.get());
     }
     else
     {
-      paraview.RegisterField("E", &E->real());
+      paraview.RegisterField("E", &E->Real());
       paraview_bdr.RegisterVCoeffField("E", Esr.get());
     }
   }
   if (B)
   {
-    if (has_imaginary)
+    if (HasImag())
     {
-      paraview.RegisterField("B_real", &B->real());
-      paraview.RegisterField("B_imag", &B->imag());
+      paraview.RegisterField("B_real", &B->Real());
+      paraview.RegisterField("B_imag", &B->Imag());
       paraview_bdr.RegisterVCoeffField("B_real", Bsr.get());
       paraview_bdr.RegisterVCoeffField("B_imag", Bsi.get());
     }
     else
     {
-      paraview.RegisterField("B", &B->real());
+      paraview.RegisterField("B", &B->Real());
       paraview_bdr.RegisterVCoeffField("B", Bsr.get());
     }
   }
   if (V)
   {
-    paraview.RegisterField("V", &*V);
+    paraview.RegisterField("V", &V->Real());
     paraview_bdr.RegisterCoeffField("V", Vs.get());
   }
   if (A)
   {
-    paraview.RegisterField("A", &*A);
+    paraview.RegisterField("A", &A->Real());
     paraview_bdr.RegisterVCoeffField("A", As.get());
   }
 
@@ -224,7 +216,7 @@ void PostOperator::InitializeDataCollection(const IoData &iodata)
   // currents are single-valued at internal boundaries.
   if (Qsr)
   {
-    if (has_imaginary)
+    if (HasImag())
     {
       paraview_bdr.RegisterCoeffField("Qs_real", Qsr.get());
       paraview_bdr.RegisterCoeffField("Qs_imag", Qsi.get());
@@ -236,7 +228,7 @@ void PostOperator::InitializeDataCollection(const IoData &iodata)
   }
   if (Jsr)
   {
-    if (has_imaginary)
+    if (HasImag())
     {
       paraview_bdr.RegisterVCoeffField("Js_real", Jsr.get());
       paraview_bdr.RegisterVCoeffField("Js_imag", Jsi.get());
@@ -270,66 +262,64 @@ void PostOperator::InitializeDataCollection(const IoData &iodata)
 
 void PostOperator::SetEGridFunction(const ComplexVector &e)
 {
-  MFEM_VERIFY(
-      has_imaginary,
-      "SetEGridFunction for complex-valued output called when has_imaginary == false!");
+  MFEM_VERIFY(HasImag(),
+              "SetEGridFunction for complex-valued output called when HasImag() == false!");
   MFEM_VERIFY(E, "Incorrect usage of PostOperator::SetEGridFunction!");
-  E->real().SetFromTrueDofs(e.Real());  // Parallel distribute
-  E->imag().SetFromTrueDofs(e.Imag());
-  E->real().ExchangeFaceNbrData();  // Ready for parallel comm on shared faces
-  E->imag().ExchangeFaceNbrData();
+  E->Real().SetFromTrueDofs(e.Real());  // Parallel distribute
+  E->Imag().SetFromTrueDofs(e.Imag());
+  E->Real().ExchangeFaceNbrData();  // Ready for parallel comm on shared faces
+  E->Imag().ExchangeFaceNbrData();
   lumped_port_init = wave_port_init = false;
 }
 
 void PostOperator::SetBGridFunction(const ComplexVector &b)
 {
-  MFEM_VERIFY(
-      has_imaginary,
-      "SetBGridFunction for complex-valued output called when has_imaginary == false!");
+  MFEM_VERIFY(HasImag(),
+              "SetBGridFunction for complex-valued output called when HasImag() == false!");
   MFEM_VERIFY(B, "Incorrect usage of PostOperator::SetBGridFunction!");
-  B->real().SetFromTrueDofs(b.Real());  // Parallel distribute
-  B->imag().SetFromTrueDofs(b.Imag());
-  B->real().ExchangeFaceNbrData();  // Ready for parallel comm on shared faces
-  B->imag().ExchangeFaceNbrData();
+  B->Real().SetFromTrueDofs(b.Real());  // Parallel distribute
+  B->Imag().SetFromTrueDofs(b.Imag());
+  B->Real().ExchangeFaceNbrData();  // Ready for parallel comm on shared faces
+  B->Imag().ExchangeFaceNbrData();
   lumped_port_init = wave_port_init = false;
 }
 
 void PostOperator::SetEGridFunction(const Vector &e)
 {
-  MFEM_VERIFY(!has_imaginary,
-              "SetEGridFunction for real-valued output called when has_imaginary == true!");
+  MFEM_VERIFY(!HasImag(),
+              "SetEGridFunction for real-valued output called when HasImag() == true!");
   MFEM_VERIFY(E, "Incorrect usage of PostOperator::SetEGridFunction!");
-  E->real().SetFromTrueDofs(e);
-  E->real().ExchangeFaceNbrData();
+  E->Real().SetFromTrueDofs(e);
+  E->Real().ExchangeFaceNbrData();
   lumped_port_init = wave_port_init = false;
 }
 
 void PostOperator::SetBGridFunction(const Vector &b)
 {
-  MFEM_VERIFY(!has_imaginary,
-              "SetBGridFunction for real-valued output called when has_imaginary == true!");
+  MFEM_VERIFY(!HasImag(),
+              "SetBGridFunction for real-valued output called when HasImag() == true!");
   MFEM_VERIFY(B, "Incorrect usage of PostOperator::SetBGridFunction!");
-  B->real().SetFromTrueDofs(b);
-  B->real().ExchangeFaceNbrData();
+  B->Real().SetFromTrueDofs(b);
+  B->Real().ExchangeFaceNbrData();
   lumped_port_init = wave_port_init = false;
 }
 
 void PostOperator::SetVGridFunction(const Vector &v)
 {
-  MFEM_VERIFY(!has_imaginary,
-              "SetVGridFunction for real-valued output called when has_imaginary == true!");
+  MFEM_VERIFY(!HasImag(),
+              "SetVGridFunction for real-valued output called when HasImag() == true!");
   MFEM_VERIFY(V, "Incorrect usage of PostOperator::SetVGridFunction!");
-  V->SetFromTrueDofs(v);
-  V->ExchangeFaceNbrData();
+  V->Real().SetFromTrueDofs(v);
+  V->Real().ExchangeFaceNbrData();
 }
 
 void PostOperator::SetAGridFunction(const Vector &a)
 {
-  MFEM_VERIFY(!has_imaginary,
-              "SetAGridFunction for real-valued output called when has_imaginary == true!");
+  MFEM_VERIFY(!HasImag(),
+              "SetAGridFunction for real-valued output called when HasImag() == true!");
   MFEM_VERIFY(A, "Incorrect usage of PostOperator::SetAGridFunction!");
-  A->SetFromTrueDofs(a);
-  A->ExchangeFaceNbrData();
+  A->Real().SetFromTrueDofs(a);
+  A->Real().ExchangeFaceNbrData();
 }
 
 void PostOperator::UpdatePorts(const LumpedPortOperator &lumped_port_op, double omega)
@@ -342,20 +332,18 @@ void PostOperator::UpdatePorts(const LumpedPortOperator &lumped_port_op, double
   for (const auto &[idx, data] : lumped_port_op)
   {
     auto &vi = lumped_port_vi[idx];
-    if (has_imaginary)
+    vi.P = data.GetPower(*E, *B);
+    vi.V = data.GetVoltage(*E);
+    if (HasImag())
     {
       MFEM_VERIFY(
           omega > 0.0,
           "Frequency domain lumped port postprocessing requires nonzero frequency!");
       vi.S = data.GetSParameter(*E);
-      vi.P = data.GetPower(*E, *B);
-      vi.V = data.GetVoltage(*E);
       vi.Z = data.GetCharacteristicImpedance(omega);
     }
     else
     {
-      vi.P = data.GetPower(E->real(), B->real());
-      vi.V = data.GetVoltage(E->real());
       vi.S = vi.Z = 0.0;
     }
   }
@@ -364,7 +352,7 @@ void PostOperator::UpdatePorts(const LumpedPortOperator &lumped_port_op, double
 
 void PostOperator::UpdatePorts(const WavePortOperator &wave_port_op, double omega)
 {
-  MFEM_VERIFY(has_imaginary && E && B, "Incorrect usage of PostOperator::UpdatePorts!");
+  MFEM_VERIFY(HasImag() && E && B, "Incorrect usage of PostOperator::UpdatePorts!");
   if (wave_port_init)
   {
     return;
@@ -374,8 +362,8 @@ void PostOperator::UpdatePorts(const WavePortOperator &wave_port_op, double omeg
     MFEM_VERIFY(omega > 0.0,
                 "Frequency domain wave port postprocessing requires nonzero frequency!");
     auto &vi = wave_port_vi[idx];
-    vi.S = data.GetSParameter(*E);
     vi.P = data.GetPower(*E, *B);
+    vi.S = data.GetSParameter(*E);
     vi.V = vi.Z = 0.0;  // Not yet implemented (Z = V² / P, I = V / Z)
   }
   wave_port_init = true;
@@ -388,8 +376,7 @@ double PostOperator::GetEFieldEnergy() const
   // voltages/currents which are 2x the time-averaged values. This correctly yields an EPR
   // of 1 in cases where expected.
   MFEM_VERIFY(E, "PostOperator is not configured for electric field energy calculation!");
-  return has_imaginary ? dom_post_op.GetElectricFieldEnergy(*E)
-                       : dom_post_op.GetElectricFieldEnergy(E->real());
+  return dom_post_op.GetElectricFieldEnergy(*E);
 }
 
 double PostOperator::GetHFieldEnergy() const
@@ -399,22 +386,19 @@ double PostOperator::GetHFieldEnergy() const
   // voltages/currents which are 2x the time-averaged values. This correctly yields an EPR
   // of 1 in cases where expected.
   MFEM_VERIFY(B, "PostOperator is not configured for magnetic field energy calculation!");
-  return has_imaginary ? dom_post_op.GetMagneticFieldEnergy(*B)
-                       : dom_post_op.GetMagneticFieldEnergy(B->real());
+  return dom_post_op.GetMagneticFieldEnergy(*B);
 }
 
 double PostOperator::GetEFieldEnergy(int idx) const
 {
   MFEM_VERIFY(E, "PostOperator is not configured for electric field energy calculation!");
-  return has_imaginary ? dom_post_op.GetDomainElectricFieldEnergy(idx, *E)
-                       : dom_post_op.GetDomainElectricFieldEnergy(idx, E->real());
+  return dom_post_op.GetDomainElectricFieldEnergy(idx, *E);
 }
 
 double PostOperator::GetHFieldEnergy(int idx) const
 {
   MFEM_VERIFY(B, "PostOperator is not configured for magnetic field energy calculation!");
-  return has_imaginary ? dom_post_op.GetDomainMagneticFieldEnergy(idx, *B)
-                       : dom_post_op.GetDomainMagneticFieldEnergy(idx, B->real());
+  return dom_post_op.GetDomainMagneticFieldEnergy(idx, *B);
 }
 
 double PostOperator::GetLumpedInductorEnergy(const LumpedPortOperator &lumped_port_op) const
@@ -594,9 +578,7 @@ double PostOperator::GetInterfaceParticipation(int idx, double Em) const
   // with:
   //          p_mj = 1/2 t_j Re{∫_{Γ_j} (ε_j E_m)ᴴ E_m dS} /(E_elec + E_cap).
   MFEM_VERIFY(E, "Surface Q not defined, no electric field solution found!");
-  double Esurf = has_imaginary
-                     ? surf_post_op.GetInterfaceElectricFieldEnergy(idx, *E)
-                     : surf_post_op.GetInterfaceElectricFieldEnergy(idx, E->real());
+  double Esurf = surf_post_op.GetInterfaceElectricFieldEnergy(idx, *E);
   return Esurf / Em;
 }
 
@@ -607,8 +589,7 @@ double PostOperator::GetSurfaceCharge(int idx) const
   // for both sides of the surface. This then yields the capacitive coupling to the
   // excitation as C_jk = Q_j / V_k where V_k is the excitation voltage.
   MFEM_VERIFY(E, "Surface capacitance not defined, no electric field solution found!");
-  double Q = has_imaginary ? surf_post_op.GetSurfaceElectricCharge(idx, *E)
-                           : surf_post_op.GetSurfaceElectricCharge(idx, E->real());
+  double Q = surf_post_op.GetSurfaceElectricCharge(idx, *E);
   return Q;
 }
 
@@ -620,8 +601,7 @@ double PostOperator::GetSurfaceFlux(int idx) const
   // which are discontinuous at interior boundary elements.
   MFEM_VERIFY(B,
               "Surface inductance not defined, no magnetic flux density solution found!");
-  double Phi = has_imaginary ? surf_post_op.GetSurfaceMagneticFlux(idx, *B)
-                             : surf_post_op.GetSurfaceMagneticFlux(idx, B->real());
+  double Phi = surf_post_op.GetSurfaceMagneticFlux(idx, *B);
   return Phi;
 }
 
@@ -631,7 +611,7 @@ void PostOperator::WriteFields(int step, double time, const ErrorIndicator *indi
   // visualization. Write the mesh coordinates in the same units as originally input.
   bool first_save = (paraview.GetCycle() < 0);
   mfem::ParMesh &mesh =
-      (E) ? *E->ParFESpace()->GetParMesh() : *B->ParFESpace()->GetParMesh();
+      HasE() ? *E->ParFESpace()->GetParMesh() : *B->ParFESpace()->GetParMesh();
   auto ScaleGridFunctions = [&mesh, this](double L)
   {
     // For fields on H(curl) and H(div) spaces, we "undo" the effect of redimensionalizing
@@ -645,31 +625,31 @@ void PostOperator::WriteFields(int step, double time, const ErrorIndicator *indi
     if (E)
     {
       // Piola transform: J^-T
-      E->real() *= L;
-      E->real().FaceNbrData() *= L;
-      if (has_imaginary)
+      E->Real() *= L;
+      E->Real().FaceNbrData() *= L;
+      if (HasImag())
       {
-        E->imag() *= L;
-        E->imag().FaceNbrData() *= L;
+        E->Imag() *= L;
+        E->Imag().FaceNbrData() *= L;
       }
     }
     if (B)
     {
       // Piola transform: J / |J|
       const auto Ld = std::pow(L, mesh.Dimension() - 1);
-      B->real() *= Ld;
-      B->real().FaceNbrData() *= Ld;
-      if (has_imaginary)
+      B->Real() *= Ld;
+      B->Real().FaceNbrData() *= Ld;
+      if (HasImag())
       {
-        B->imag() *= Ld;
-        B->imag().FaceNbrData() *= Ld;
+        B->Imag() *= Ld;
+        B->Imag().FaceNbrData() *= Ld;
       }
     }
     if (A)
     {
       // Piola transform: J^-T
-      *A *= L;
-      A->FaceNbrData() *= L;
+      A->Real() *= L;
+      A->Real().FaceNbrData() *= L;
     }
   };
   mesh::DimensionalizeMesh(mesh, mesh_Lc0);
@@ -724,13 +704,13 @@ void PostOperator::WriteFields(int step, double time, const ErrorIndicator *indi
 std::vector<std::complex<double>> PostOperator::ProbeEField() const
 {
   MFEM_VERIFY(E, "PostOperator is not configured for electric field probes!");
-  return interp_op.ProbeField(*E, has_imaginary);
+  return interp_op.ProbeField(*E);
 }
 
 std::vector<std::complex<double>> PostOperator::ProbeBField() const
 {
   MFEM_VERIFY(B, "PostOperator is not configured for magnetic flux density probes!");
-  return interp_op.ProbeField(*B, has_imaginary);
+  return interp_op.ProbeField(*B);
 }
 
 }  // namespace palace
diff --git a/palace/models/postoperator.hpp b/palace/models/postoperator.hpp
index 10de7c8c5..7d90cd702 100644
--- a/palace/models/postoperator.hpp
+++ b/palace/models/postoperator.hpp
@@ -11,6 +11,7 @@
 #include <string>
 #include <vector>
 #include <mfem.hpp>
+#include "fem/gridfunction.hpp"
 #include "fem/interpolator.hpp"
 #include "linalg/operator.hpp"
 #include "linalg/vector.hpp"
@@ -44,9 +45,7 @@ class PostOperator
   const DomainPostOperator dom_post_op;
 
   // Objects for grid function postprocessing from the FE solution.
-  const bool has_imaginary;
-  mutable std::optional<mfem::ParComplexGridFunction> E, B;
-  mutable std::optional<mfem::ParGridFunction> V, A;
+  mutable std::optional<GridFunction> E, B, V, A;
   std::unique_ptr<mfem::VectorCoefficient> Esr, Esi, Bsr, Bsi, As, Jsr, Jsi;
   std::unique_ptr<mfem::Coefficient> Vs, Ue, Um, Qsr, Qsi;
 
@@ -82,9 +81,9 @@ class PostOperator
   const auto &GetDomainPostOp() const { return dom_post_op; }
 
   // Return options for postprocessing configuration.
-  bool HasImaginary() const { return has_imaginary; }
   bool HasE() const { return E.has_value(); }
   bool HasB() const { return B.has_value(); }
+  bool HasImag() const { return HasE() ? E->HasImag() : B->HasImag(); }
 
   // Populate the grid function solutions for the E- and B-field using the solution vectors
   // on the true dofs. For the real-valued overload, the electric scalar potential can be
@@ -130,8 +129,8 @@ class PostOperator
                                      int source_idx) const;
 
   // Postprocess the circuit voltage and current across lumped port index using the electric
-  // field solution. When has_imaginary is false, the returned voltage has only a nonzero
-  // real part.
+  // field solution. When the internal grid functions are real-valued, the returned voltage
+  // has only a nonzero real part.
   std::complex<double> GetPortPower(const LumpedPortOperator &lumped_port_op,
                                     int idx) const;
   std::complex<double> GetPortPower(const WavePortOperator &wave_port_op, int idx) const;
@@ -174,8 +173,9 @@ class PostOperator
 
   // Probe the E- and B-fields for their vector-values at speceified locations in space.
   // Locations of probes are set up in constructor from configuration file data. If
-  // has_imaginary is false, the returned fields have only nonzero real parts. Output
-  // vectors are ordered by vector dimension, that is [v1x, v1y, v1z, v2x, v2y, v2z, ...].
+  // the internal grid functions are real-valued, the returned fields have only nonzero real
+  // parts. Output vectors are ordered by vector dimension, that is [v1x, v1y, v1z, v2x,
+  // v2y, v2z, ...].
   const auto &GetProbes() const { return interp_op.GetProbes(); }
   std::vector<std::complex<double>> ProbeEField() const;
   std::vector<std::complex<double>> ProbeBField() const;
diff --git a/palace/models/surfacepostoperator.cpp b/palace/models/surfacepostoperator.cpp
index 6fea74d98..94f3e8b9f 100644
--- a/palace/models/surfacepostoperator.cpp
+++ b/palace/models/surfacepostoperator.cpp
@@ -4,6 +4,7 @@
 #include "surfacepostoperator.hpp"
 
 #include <complex>
+#include "fem/gridfunction.hpp"
 #include "fem/integrator.hpp"
 #include "linalg/vector.hpp"
 #include "models/materialoperator.hpp"
@@ -171,78 +172,49 @@ double SurfacePostOperator::GetInterfaceLossTangent(int idx) const
   return it->second.tandelta;
 }
 
-double SurfacePostOperator::GetInterfaceElectricFieldEnergy(
-    int idx, const mfem::ParComplexGridFunction &E) const
+double SurfacePostOperator::GetInterfaceElectricFieldEnergy(int idx,
+                                                            const GridFunction &E) const
 {
   auto it = eps_surfs.find(idx);
   MFEM_VERIFY(it != eps_surfs.end(),
               "Unknown dielectric loss postprocessing surface index requested!");
-  double dot = GetLocalSurfaceIntegral(it->second, E.real()) +
-               GetLocalSurfaceIntegral(it->second, E.imag());
-  Mpi::GlobalSum(1, &dot, E.ParFESpace()->GetComm());
-  return dot;
-}
-
-double
-SurfacePostOperator::GetInterfaceElectricFieldEnergy(int idx,
-                                                     const mfem::ParGridFunction &E) const
-{
-  auto it = eps_surfs.find(idx);
-  MFEM_VERIFY(it != eps_surfs.end(),
-              "Unknown dielectric loss postprocessing surface index requested!");
-  double dot = GetLocalSurfaceIntegral(it->second, E);
-  Mpi::GlobalSum(1, &dot, E.ParFESpace()->GetComm());
+  double dot = GetLocalSurfaceIntegral(it->second, E.Real());
+  if (E.HasImag())
+  {
+    dot += GetLocalSurfaceIntegral(it->second, E.Imag());
+  }
+  Mpi::GlobalSum(1, &dot, E.GetComm());
   return dot;
 }
 
-double
-SurfacePostOperator::GetSurfaceElectricCharge(int idx,
-                                              const mfem::ParComplexGridFunction &E) const
+double SurfacePostOperator::GetSurfaceElectricCharge(int idx, const GridFunction &E) const
 {
   auto it = charge_surfs.find(idx);
   MFEM_VERIFY(it != charge_surfs.end(),
               "Unknown capacitance postprocessing surface index requested!");
-  std::complex<double> dot(GetLocalSurfaceIntegral(it->second, E.real()),
-                           GetLocalSurfaceIntegral(it->second, E.imag()));
-  Mpi::GlobalSum(1, &dot, E.ParFESpace()->GetComm());
+  std::complex<double> dot(GetLocalSurfaceIntegral(it->second, E.Real()), 0.0);
+  if (E.HasImag())
+  {
+    dot.imag(GetLocalSurfaceIntegral(it->second, E.Imag()));
+  }
+  Mpi::GlobalSum(1, &dot, E.GetComm());
   return std::copysign(std::abs(dot), dot.real());
 }
 
-double SurfacePostOperator::GetSurfaceElectricCharge(int idx,
-                                                     const mfem::ParGridFunction &E) const
-{
-  auto it = charge_surfs.find(idx);
-  MFEM_VERIFY(it != charge_surfs.end(),
-              "Unknown capacitance postprocessing surface index requested!");
-  double dot = GetLocalSurfaceIntegral(it->second, E);
-  Mpi::GlobalSum(1, &dot, E.ParFESpace()->GetComm());
-  return dot;
-}
-
-double
-SurfacePostOperator::GetSurfaceMagneticFlux(int idx,
-                                            const mfem::ParComplexGridFunction &B) const
+double SurfacePostOperator::GetSurfaceMagneticFlux(int idx, const GridFunction &B) const
 {
   auto it = flux_surfs.find(idx);
   MFEM_VERIFY(it != flux_surfs.end(),
               "Unknown inductance postprocessing surface index requested!");
-  std::complex<double> dot(GetLocalSurfaceIntegral(it->second, B.real()),
-                           GetLocalSurfaceIntegral(it->second, B.imag()));
-  Mpi::GlobalSum(1, &dot, B.ParFESpace()->GetComm());
+  std::complex<double> dot(GetLocalSurfaceIntegral(it->second, B.Real()), 0.0);
+  if (B.HasImag())
+  {
+    dot.imag(GetLocalSurfaceIntegral(it->second, B.Imag()));
+  }
+  Mpi::GlobalSum(1, &dot, B.GetComm());
   return std::copysign(std::abs(dot), dot.real());
 }
 
-double SurfacePostOperator::GetSurfaceMagneticFlux(int idx,
-                                                   const mfem::ParGridFunction &B) const
-{
-  auto it = flux_surfs.find(idx);
-  MFEM_VERIFY(it != flux_surfs.end(),
-              "Unknown inductance postprocessing surface index requested!");
-  double dot = GetLocalSurfaceIntegral(it->second, B);
-  Mpi::GlobalSum(1, &dot, B.ParFESpace()->GetComm());
-  return dot;
-}
-
 double SurfacePostOperator::GetLocalSurfaceIntegral(const SurfaceData &data,
                                                     const mfem::ParGridFunction &U) const
 {
diff --git a/palace/models/surfacepostoperator.hpp b/palace/models/surfacepostoperator.hpp
index ab4db12b0..f7e936f80 100644
--- a/palace/models/surfacepostoperator.hpp
+++ b/palace/models/surfacepostoperator.hpp
@@ -13,6 +13,7 @@
 namespace palace
 {
 
+class GridFunction;
 class IoData;
 class MaterialOperator;
 
@@ -100,13 +101,9 @@ class SurfacePostOperator
   // Get surface integrals computing dielectric interface energy, surface charge, or
   // surface magnetic flux.
   double GetInterfaceLossTangent(int idx) const;
-  double GetInterfaceElectricFieldEnergy(int idx,
-                                         const mfem::ParComplexGridFunction &E) const;
-  double GetInterfaceElectricFieldEnergy(int idx, const mfem::ParGridFunction &E) const;
-  double GetSurfaceElectricCharge(int idx, const mfem::ParComplexGridFunction &E) const;
-  double GetSurfaceElectricCharge(int idx, const mfem::ParGridFunction &E) const;
-  double GetSurfaceMagneticFlux(int idx, const mfem::ParComplexGridFunction &B) const;
-  double GetSurfaceMagneticFlux(int idx, const mfem::ParGridFunction &B) const;
+  double GetInterfaceElectricFieldEnergy(int idx, const GridFunction &E) const;
+  double GetSurfaceElectricCharge(int idx, const GridFunction &E) const;
+  double GetSurfaceMagneticFlux(int idx, const GridFunction &B) const;
 };
 
 }  // namespace palace
diff --git a/palace/models/waveportoperator.cpp b/palace/models/waveportoperator.cpp
index bbc2b7bf9..7b9360bde 100644
--- a/palace/models/waveportoperator.cpp
+++ b/palace/models/waveportoperator.cpp
@@ -26,9 +26,6 @@ namespace palace
 
 using namespace std::complex_literals;
 
-// XX TODO: All these tiny solves should happen only on the CPU (if we can configure Hypre
-//          device at runtime).
-
 namespace
 {
 
@@ -284,9 +281,8 @@ ComplexHypreParMatrix GetSystemMatrixB(mfem::HypreParMatrix *Bttr,
   return {std::move(Br), std::move(Bi)};
 }
 
-void Normalize(const mfem::ParGridFunction &S0t, mfem::ParComplexGridFunction &E0t,
-               mfem::ParComplexGridFunction &E0n, mfem::LinearForm &sr,
-               mfem::LinearForm &si)
+void Normalize(const GridFunction &S0t, GridFunction &E0t, GridFunction &E0n,
+               mfem::LinearForm &sr, mfem::LinearForm &si)
 {
   // Normalize grid functions to a chosen polarization direction and unit power, |E x H⋆| ⋅
   // n, integrated over the port surface (+n is the direction of propagation). The n x H
@@ -297,20 +293,20 @@ void Normalize(const mfem::ParGridFunction &S0t, mfem::ParComplexGridFunction &E
   // |E x H⋆| ⋅ n = |E ⋅ (-n x H⋆)|. This also updates the n x H coefficients depending on
   // Et, En. Update linear forms for postprocessing too.
   std::complex<double> dot[2] = {
-      {sr * S0t, si * S0t},
-      {-(sr * E0t.real()) - (si * E0t.imag()), -(sr * E0t.imag()) + (si * E0t.real())}};
+      {sr * S0t.Real(), si * S0t.Real()},
+      {-(sr * E0t.Real()) - (si * E0t.Imag()), -(sr * E0t.Imag()) + (si * E0t.Real())}};
   Mpi::GlobalSum(2, dot, S0t.ParFESpace()->GetComm());
   auto scale = std::abs(dot[0]) / (dot[0] * std::sqrt(std::abs(dot[1])));
-  ComplexVector::AXPBY(scale, E0t.real(), E0t.imag(), 0.0, E0t.real(), E0t.imag());
-  ComplexVector::AXPBY(scale, E0n.real(), E0n.imag(), 0.0, E0n.real(), E0n.imag());
+  ComplexVector::AXPBY(scale, E0t.Real(), E0t.Imag(), 0.0, E0t.Real(), E0t.Imag());
+  ComplexVector::AXPBY(scale, E0n.Real(), E0n.Imag(), 0.0, E0n.Real(), E0n.Imag());
   ComplexVector::AXPBY(scale, sr, si, 0.0, sr, si);
 
   // This parallel communication is not required since wave port boundaries are true one-
   // sided boundaries.
-  // E0t.real().ExchangeFaceNbrData();  // Ready for parallel comm on shared faces for n x H
-  // E0t.imag().ExchangeFaceNbrData();  // coefficients evaluation
-  // E0n.real().ExchangeFaceNbrData();
-  // E0n.imag().ExchangeFaceNbrData();
+  // E0t.Real().ExchangeFaceNbrData();  // Ready for parallel comm on shared faces for n x H
+  // E0t.Imag().ExchangeFaceNbrData();  // coefficients evaluation
+  // E0n.Real().ExchangeFaceNbrData();
+  // E0n.Imag().ExchangeFaceNbrData();
 }
 
 // Helper for BdrSubmeshEVectorCoefficient and BdrSubmeshHVectorCoefficient.
@@ -325,17 +321,16 @@ template <ValueType Type>
 class BdrSubmeshEVectorCoefficient : public mfem::VectorCoefficient
 {
 private:
-  const mfem::ParComplexGridFunction &Et, &En;
+  const GridFunction &Et, &En;
   const mfem::ParSubMesh &submesh;
   const std::unordered_map<int, int> &submesh_parent_elems;
   mfem::IsoparametricTransformation T_loc;
 
 public:
-  BdrSubmeshEVectorCoefficient(const mfem::ParComplexGridFunction &Et,
-                               const mfem::ParComplexGridFunction &En,
+  BdrSubmeshEVectorCoefficient(const GridFunction &Et, const GridFunction &En,
                                const mfem::ParSubMesh &submesh,
                                const std::unordered_map<int, int> &submesh_parent_elems)
-    : mfem::VectorCoefficient(Et.real().VectorDim()), Et(Et), En(En), submesh(submesh),
+    : mfem::VectorCoefficient(Et.Real().VectorDim()), Et(Et), En(En), submesh(submesh),
       submesh_parent_elems(submesh_parent_elems)
   {
   }
@@ -383,14 +378,14 @@ class BdrSubmeshEVectorCoefficient : public mfem::VectorCoefficient
     BdrGridFunctionCoefficient::GetNormal(*T_submesh, nor);
     if constexpr (Type == ValueType::REAL)
     {
-      Et.real().GetVectorValue(*T_submesh, ip, V);
-      auto Vn = En.real().GetValue(*T_submesh, ip);
+      Et.Real().GetVectorValue(*T_submesh, ip, V);
+      auto Vn = En.Real().GetValue(*T_submesh, ip);
       V.Add(-Vn, nor);
     }
     else
     {
-      Et.imag().GetVectorValue(*T_submesh, ip, V);
-      auto Vn = En.imag().GetValue(*T_submesh, ip);
+      Et.Imag().GetVectorValue(*T_submesh, ip, V);
+      auto Vn = En.Imag().GetValue(*T_submesh, ip);
       V.Add(-Vn, nor);
     }
   }
@@ -403,7 +398,7 @@ template <ValueType Type>
 class BdrSubmeshHVectorCoefficient : public mfem::VectorCoefficient
 {
 private:
-  const mfem::ParComplexGridFunction &Et, &En;
+  const GridFunction &Et, &En;
   const MaterialOperator &mat_op;
   const mfem::ParSubMesh &submesh;
   const std::unordered_map<int, int> &submesh_parent_elems;
@@ -412,13 +407,12 @@ class BdrSubmeshHVectorCoefficient : public mfem::VectorCoefficient
   double omega;
 
 public:
-  BdrSubmeshHVectorCoefficient(const mfem::ParComplexGridFunction &Et,
-                               const mfem::ParComplexGridFunction &En,
+  BdrSubmeshHVectorCoefficient(const GridFunction &Et, const GridFunction &En,
                                const MaterialOperator &mat_op,
                                const mfem::ParSubMesh &submesh,
                                const std::unordered_map<int, int> &submesh_parent_elems,
                                std::complex<double> kn, double omega)
-    : mfem::VectorCoefficient(Et.real().VectorDim()), Et(Et), En(En), mat_op(mat_op),
+    : mfem::VectorCoefficient(Et.Real().VectorDim()), Et(Et), En(En), mat_op(mat_op),
       submesh(submesh), submesh_parent_elems(submesh_parent_elems), kn(kn), omega(omega)
   {
   }
@@ -492,20 +486,20 @@ class BdrSubmeshHVectorCoefficient : public mfem::VectorCoefficient
     mfem::Vector U;
     if constexpr (Type == ValueType::REAL)
     {
-      Et.real().GetVectorValue(*T_submesh, ip, U);
+      Et.Real().GetVectorValue(*T_submesh, ip, U);
       U *= -kn.real();
 
       mfem::Vector dU;
-      En.imag().GetGradient(*T_submesh, dU);
+      En.Imag().GetGradient(*T_submesh, dU);
       U -= dU;
     }
     else
     {
-      Et.imag().GetVectorValue(*T_submesh, ip, U);
+      Et.Imag().GetVectorValue(*T_submesh, ip, U);
       U *= -kn.real();
 
       mfem::Vector dU;
-      En.real().GetGradient(*T_submesh, dU);
+      En.Real().GetGradient(*T_submesh, dU);
       U += dU;
     }
 
@@ -547,15 +541,15 @@ WavePortData::WavePortData(const config::WavePortData &data,
   port_nd_fespace = std::make_unique<FiniteElementSpace>(*port_mesh, port_nd_fec.get());
   port_h1_fespace = std::make_unique<FiniteElementSpace>(*port_mesh, port_h1_fec.get());
 
-  mfem::ParGridFunction E0t(&nd_fespace);
-  mfem::ParGridFunction E0n(&h1_fespace);
-  port_E0t = std::make_unique<mfem::ParComplexGridFunction>(&port_nd_fespace->Get());
-  port_E0n = std::make_unique<mfem::ParComplexGridFunction>(&port_h1_fespace->Get());
+  GridFunction E0t(nd_fespace), E0n(h1_fespace);
+  port_E0t = std::make_unique<GridFunction>(*port_nd_fespace, true);
+  port_E0n = std::make_unique<GridFunction>(*port_h1_fespace, true);
+  port_E = std::make_unique<GridFunction>(*port_nd_fespace, true);
 
   port_nd_transfer = std::make_unique<mfem::ParTransferMap>(
-      mfem::ParSubMesh::CreateTransferMap(E0t, port_E0t->real()));
+      mfem::ParSubMesh::CreateTransferMap(E0t.Real(), port_E0t->Real()));
   port_h1_transfer = std::make_unique<mfem::ParTransferMap>(
-      mfem::ParSubMesh::CreateTransferMap(E0n, port_E0n->real()));
+      mfem::ParSubMesh::CreateTransferMap(E0n.Real(), port_E0n->Real()));
 
   // Construct mapping from parent (boundary) element indices to submesh (domain)
   // elements.
@@ -571,9 +565,9 @@ WavePortData::WavePortData(const config::WavePortData &data,
   // Extract Dirichlet BC true dofs for the port FE spaces.
   {
     mfem::Array<int> port_nd_dbc_tdof_list, port_h1_dbc_tdof_list;
-    GetEssentialTrueDofs(E0t, E0n, port_E0t->real(), port_E0n->real(), *port_nd_transfer,
-                         *port_h1_transfer, dbc_attr, port_nd_dbc_tdof_list,
-                         port_h1_dbc_tdof_list);
+    GetEssentialTrueDofs(E0t.Real(), E0n.Real(), port_E0t->Real(), port_E0n->Real(),
+                         *port_nd_transfer, *port_h1_transfer, dbc_attr,
+                         port_nd_dbc_tdof_list, port_h1_dbc_tdof_list);
     int nd_tdof_offset = port_nd_fespace->GetTrueVSize();
     port_dbc_tdof_list.Reserve(port_nd_dbc_tdof_list.Size() + port_h1_dbc_tdof_list.Size());
     for (auto tdof : port_nd_dbc_tdof_list)
@@ -623,7 +617,7 @@ WavePortData::WavePortData(const config::WavePortData &data,
         port_h1_fespace->Get().GetTrueDofOffsets(), &diag);
     auto [Bttr, Btti] = GetBtt(mat_op, *port_nd_fespace);
     auto [Br, Bi] = GetSystemMatrixB(Bttr.get(), Btti.get(), Dnn.get(), port_dbc_tdof_list);
-    B = std::make_unique<ComplexWrapperOperator>(std::move(Br), std::move(Bi));
+    B0 = std::make_unique<ComplexWrapperOperator>(std::move(Br), std::move(Bi));
   }
 
   // Configure a communicator for the processes which have elements for this port.
@@ -822,8 +816,8 @@ WavePortData::WavePortData(const config::WavePortData &data,
       }
     };
     mfem::VectorFunctionCoefficient tfunc(dim, TDirection);
-    port_S0t = std::make_unique<mfem::ParGridFunction>(&port_nd_fespace->Get());
-    port_S0t->ProjectCoefficient(tfunc);
+    port_S0t = std::make_unique<GridFunction>(*port_nd_fespace);
+    port_S0t->Real().ProjectCoefficient(tfunc);
   }
 }
 
@@ -866,7 +860,7 @@ void WavePortData::Initialize(double omega)
   {
     ComplexWrapperOperator P(A->Real(), nullptr);  // Non-owning constructor
     ksp->SetOperators(*A, P);
-    eigen->SetOperators(*B, *A, EigenvalueSolver::ScaleType::NONE);
+    eigen->SetOperators(*B0, *A, EigenvalueSolver::ScaleType::NONE);
     eigen->SetInitialSpace(v0);
     int num_conv = eigen->Solve();
     MFEM_VERIFY(num_conv >= mode_idx, "Wave port eigensolver did not converge!");
@@ -907,10 +901,10 @@ void WavePortData::Initialize(double omega)
     e0ti.UseDevice(true);
     e0ni.UseDevice(true);
     ComplexVector::AXPBY(1.0 / (1i * kn0), e0nr, e0ni, 0.0, e0nr, e0ni);
-    port_E0t->real().SetFromTrueDofs(e0tr);  // Parallel distribute
-    port_E0t->imag().SetFromTrueDofs(e0ti);
-    port_E0n->real().SetFromTrueDofs(e0nr);
-    port_E0n->imag().SetFromTrueDofs(e0ni);
+    port_E0t->Real().SetFromTrueDofs(e0tr);  // Parallel distribute
+    port_E0t->Imag().SetFromTrueDofs(e0ti);
+    port_E0n->Real().SetFromTrueDofs(e0nr);
+    port_E0n->Imag().SetFromTrueDofs(e0ni);
   }
 
   // Configure the linear forms for computing S-parameters (projection of the field onto the
@@ -981,30 +975,19 @@ double WavePortData::GetExcitationPower() const
   return excitation ? 1.0 : 0.0;
 }
 
-std::complex<double> WavePortData::GetSParameter(mfem::ParComplexGridFunction &E) const
-{
-  // Compute port S-parameter, or the projection of the field onto the port mode:
-  // (E x H_inc⋆) ⋅ n = E ⋅ (-n x H_inc⋆), integrated over the port surface.
-  mfem::ParComplexGridFunction port_E(&port_nd_fespace->Get());
-  port_nd_transfer->Transfer(E.real(), port_E.real());
-  port_nd_transfer->Transfer(E.imag(), port_E.imag());
-  std::complex<double> dot(-((*port_sr) * port_E.real()) - ((*port_si) * port_E.imag()),
-                           -((*port_sr) * port_E.imag()) + ((*port_si) * port_E.real()));
-  Mpi::GlobalSum(1, &dot, port_nd_fespace->GetComm());
-  return dot;
-}
-
-std::complex<double> WavePortData::GetPower(mfem::ParComplexGridFunction &E,
-                                            mfem::ParComplexGridFunction &B) const
+std::complex<double> WavePortData::GetPower(GridFunction &E, GridFunction &B) const
 {
   // Compute port power, (E x H) ⋅ n = E ⋅ (-n x H), integrated over the port surface
   // using the computed E and H = μ⁻¹ B fields. The linear form is reconstructed from
   // scratch each time due to changing H. The BdrCurrentVectorCoefficient computes -n x H,
   // where n is an outward normal.
+  MFEM_VERIFY(E.HasImag() && B.HasImag(),
+              "Wave ports expect complex-valued E and B fields in port power "
+              "calculation!");
   auto &nd_fespace = *E.ParFESpace();
   const auto &mesh = *nd_fespace.GetParMesh();
-  BdrCurrentVectorCoefficient nxHr_func(B.real(), mat_op);
-  BdrCurrentVectorCoefficient nxHi_func(B.imag(), mat_op);
+  BdrCurrentVectorCoefficient nxHr_func(B.Real(), mat_op);
+  BdrCurrentVectorCoefficient nxHi_func(B.Imag(), mat_op);
   int bdr_attr_max = mesh.bdr_attributes.Size() ? mesh.bdr_attributes.Max() : 0;
   mfem::Array<int> attr_marker = mesh::AttrToMarker(bdr_attr_max, attr_list);
   mfem::LinearForm pr(&nd_fespace), pi(&nd_fespace);
@@ -1018,12 +1001,27 @@ std::complex<double> WavePortData::GetPower(mfem::ParComplexGridFunction &E,
   pi.Assemble();
   pr.UseDevice(true);
   pi.UseDevice(true);
-  std::complex<double> dot(-(pr * E.real()) - (pi * E.imag()),
-                           -(pr * E.imag()) + (pi * E.real()));
+  std::complex<double> dot(-(pr * E.Real()) - (pi * E.Imag()),
+                           -(pr * E.Imag()) + (pi * E.Real()));
   Mpi::GlobalSum(1, &dot, nd_fespace.GetComm());
   return dot;
 }
 
+std::complex<double> WavePortData::GetSParameter(GridFunction &E) const
+{
+  // Compute port S-parameter, or the projection of the field onto the port mode:
+  // (E x H_inc⋆) ⋅ n = E ⋅ (-n x H_inc⋆), integrated over the port surface.
+  MFEM_VERIFY(E.HasImag(),
+              "Wave ports expect complex-valued E and B fields in port S-parameter "
+              "calculation!");
+  port_nd_transfer->Transfer(E.Real(), port_E->Real());
+  port_nd_transfer->Transfer(E.Imag(), port_E->Imag());
+  std::complex<double> dot(-((*port_sr) * port_E->Real()) - ((*port_si) * port_E->Imag()),
+                           -((*port_sr) * port_E->Imag()) + ((*port_si) * port_E->Real()));
+  Mpi::GlobalSum(1, &dot, port_nd_fespace->GetComm());
+  return dot;
+}
+
 WavePortOperator::WavePortOperator(const IoData &iodata, const MaterialOperator &mat_op,
                                    mfem::ParFiniteElementSpace &nd_fespace,
                                    mfem::ParFiniteElementSpace &h1_fespace)
diff --git a/palace/models/waveportoperator.hpp b/palace/models/waveportoperator.hpp
index dc28df234..3b703e97a 100644
--- a/palace/models/waveportoperator.hpp
+++ b/palace/models/waveportoperator.hpp
@@ -10,6 +10,7 @@
 #include <unordered_map>
 #include <mfem.hpp>
 #include "fem/fespace.hpp"
+#include "fem/gridfunction.hpp"
 #include "fem/mesh.hpp"
 #include "linalg/eps.hpp"
 #include "linalg/ksp.hpp"
@@ -65,7 +66,7 @@ class WavePortData
 
   // Operator storage for repeated boundary mode eigenvalue problem solves.
   std::unique_ptr<mfem::HypreParMatrix> Atnr, Atni, Antr, Anti, Annr, Anni;
-  std::unique_ptr<ComplexOperator> B;
+  std::unique_ptr<ComplexOperator> B0;
   ComplexVector v0, e0;
 
   // Eigenvalue solver for boundary modes.
@@ -74,12 +75,10 @@ class WavePortData
   std::unique_ptr<EigenvalueSolver> eigen;
   std::unique_ptr<ComplexKspSolver> ksp;
 
-  // Grid functions storing the last computed electric field mode on the port.
-  std::unique_ptr<mfem::ParComplexGridFunction> port_E0t, port_E0n;
-
-  // Stored objects for computing functions of the port modes for use as an excitation or
-  // in postprocessing.
-  std::unique_ptr<mfem::ParGridFunction> port_S0t;
+  // Grid functions storing the last computed electric field mode on the port, and stored
+  // objects for computing functions of the port modes for use as an excitation or in
+  // postprocessing.
+  std::unique_ptr<GridFunction> port_E0t, port_E0n, port_S0t, port_E;
   std::unique_ptr<mfem::LinearForm> port_sr, port_si;
 
 public:
@@ -114,10 +113,9 @@ class WavePortData
     return 0.0;
   }
 
-  std::complex<double> GetSParameter(mfem::ParComplexGridFunction &E) const;
-  std::complex<double> GetPower(mfem::ParComplexGridFunction &E,
-                                mfem::ParComplexGridFunction &B) const;
-  std::complex<double> GetVoltage(mfem::ParComplexGridFunction &E) const
+  std::complex<double> GetPower(GridFunction &E, GridFunction &B) const;
+  std::complex<double> GetSParameter(GridFunction &E) const;
+  std::complex<double> GetVoltage(GridFunction &E) const
   {
     MFEM_ABORT("GetVoltage is not yet implemented for wave port boundaries!");
     return 0.0;