Merge branch 'main' into gbozzola/petsc_metis

Sbozzolo · web-flow · commit f9adb70cab63 · 2026-04-10T14:26:40.000-07:00
diff --git a/CHANGELOG.md b/CHANGELOG.md
@@ -33,6 +33,12 @@ The format of this changelog is based on
   - Fixed linear solver stalling in the nonlinear eigenvalue solver by adding an absolute
     tolerance and adapting the relative tolerance to the outer residual
     [PR 694](https://github.com/awslabs/palace/pull/694).
+  - Replaced boundary mass matrix CG solve workaround for lumped port excitation
+    projection with direct interpolation via MFEM's
+    `ProjectBdrCoefficientTangent`. An MFEM patch fixes the underlying face DOF
+    orientation bug for Nedelec elements at order >= 2 in parallel, removing the
+    need for the expensive mass solve
+    [PR 684](https://github.com/awslabs/palace/pull/684).
 
 ## [0.16.0] - 2026-03-05
 
diff --git a/cmake/ExternalMFEM.cmake b/cmake/ExternalMFEM.cmake
@@ -398,6 +398,7 @@ message(STATUS "MFEM_OPTIONS: ${MFEM_OPTIONS_PRINT}")
 set(MFEM_PATCH_FILES
   "${CMAKE_SOURCE_DIR}/extern/patch/mfem/patch_par_tet_mesh_fix_dev.diff"
   "${CMAKE_SOURCE_DIR}/extern/patch/mfem/patch_gmsh_parser_performance.diff"
+  "${CMAKE_SOURCE_DIR}/extern/patch/mfem/mfem_pr5280.diff"
 )
 
 include(ExternalProject)
diff --git a/extern/patch/mfem/mfem_pr5280.diff b/extern/patch/mfem/mfem_pr5280.diff
@@ -0,0 +1,171 @@
+diff --git a/fem/pfespace.cpp b/fem/pfespace.cpp
+index 7bf4ea162..0b733a630 100644
+--- a/fem/pfespace.cpp
++++ b/fem/pfespace.cpp
+@@ -835,6 +835,38 @@ void ParFiniteElementSpace::CheckNDSTriaDofs()
+    nd_strias = glb_nd_strias > 0;
+ }
+ 
++void ParFiniteElementSpace::GetSharedTriFaceDofOrientations(
++   Array<int> &ltori, Array<int> &ldsize) const
++{
++   const int n = GetNDofs();
++   ltori.SetSize(n);  ltori = 0;
++   ldsize.SetSize(n); ldsize = 0;
++   if (fec->GetOrder() <= 1 || pmesh->Dimension() != 3) { return; }
++
++   // Ensure face orientations have been communicated.
++   pmesh->ExchangeFaceNbrData();
++
++   const int ngrps = pmesh->GetNGroups();
++   const int nedofs = fec->DofForGeometry(Geometry::SEGMENT);
++   Array<int> sdofs;
++   for (int g = 1; g < ngrps; g++)
++   {
++      if (pmesh->gtopo.IAmMaster(g)) { continue; }
++      for (int fi = 0; fi < pmesh->GroupNTriangles(g); fi++)
++      {
++         int face, ori, info1, info2;
++         pmesh->GroupTriangle(g, fi, face, ori);
++         pmesh->GetFaceInfos(face, &info1, &info2);
++         GetSharedTriangleDofs(g, fi, sdofs);
++         for (int j = 3*nedofs; j < sdofs.Size(); j++)
++         {
++            ldsize[sdofs[j]] = 2;
++            ltori[sdofs[j]]  = info2 % 64;
++         }
++      }
++   }
++}
++
+ void ParFiniteElementSpace::Build_Dof_TrueDof_Matrix() const // matrix P
+ {
+    MFEM_ASSERT(Conforming(), "wrong code path");
+@@ -909,9 +941,15 @@ void ParFiniteElementSpace::Build_Dof_TrueDof_Matrix() const // matrix P
+       pmesh->ExchangeFaceNbrData();
+ 
+       // Locate and count non-zeros in off-diagonal portion of P
+-      int nnz_offd = 0;
+       Array<int> ldsize(ldof); ldsize = 0;
+-      Array<int> ltori(ldof);  ltori = 0; // Local triangle orientations
++      Array<int> ltori; // Local triangle orientations
++      GetSharedTriFaceDofOrientations(ltori, ldsize);
++      // Count nnz_offd from triangle face DOFs already marked by the above call.
++      int nnz_offd = 0;
++      for (int ii = 0; ii < ldof; ii++)
++      {
++         if (ldsize[ii] == 2) { nnz_offd += 2; }
++      }
+       {
+          int ngrps = pmesh->GetNGroups();
+          int nedofs = fec->DofForGeometry(Geometry::SEGMENT);
+@@ -934,9 +972,6 @@ void ParFiniteElementSpace::Build_Dof_TrueDof_Matrix() const // matrix P
+             }
+             for (int fi=0; fi<pmesh->GroupNTriangles(g); fi++)
+             {
+-               int face, ori, info1, info2;
+-               pmesh->GroupTriangle(g, fi, face, ori);
+-               pmesh->GetFaceInfos(face, &info1, &info2);
+                this->GetSharedTriangleDofs(g, fi, sdofs);
+                for (int i=0; i<3*nedofs; i++)
+                {
+@@ -944,12 +979,7 @@ void ParFiniteElementSpace::Build_Dof_TrueDof_Matrix() const // matrix P
+                   if (ldsize[ind] == 0) { nnz_offd++; }
+                   ldsize[ind] = 1;
+                }
+-               for (int i=3*nedofs; i<sdofs.Size(); i++)
+-               {
+-                  if (ldsize[sdofs[i]] == 0) { nnz_offd += 2; }
+-                  ldsize[sdofs[i]] = 2;
+-                  ltori[sdofs[i]]  = info2 % 64;
+-               }
++               // Face DOFs already counted by GetSharedTriFaceDofOrientations.
+             }
+             for (int fi=0; fi<pmesh->GroupNQuadrilaterals(g); fi++)
+             {
+diff --git a/fem/pfespace.hpp b/fem/pfespace.hpp
+index 4dc1fec40..d989c26c9 100644
+--- a/fem/pfespace.hpp
++++ b/fem/pfespace.hpp
+@@ -340,6 +340,13 @@ public:
+ 
+    inline ParMesh *GetParMesh() const { return pmesh; }
+ 
++   /** @brief Compute shared triangle face DOF orientations for ND spaces.
++       For each L-DOF that is a face DOF on a shared triangular face (non-master
++       side), sets ldsize[dof] = 2 and ltori[dof] = face orientation relative to
++       the master rank. Used for orientation correction around Reduce/Bcast. */
++   void GetSharedTriFaceDofOrientations(Array<int> &ltori,
++                                        Array<int> &ldsize) const;
++
+    int GetDofSign(int i)
+    { return NURBSext || Nonconforming() ? 1 : ldof_sign[VDofToDof(i)]; }
+    HYPRE_BigInt *GetDofOffsets()     const { return dof_offsets; }
+diff --git a/fem/pgridfunc.cpp b/fem/pgridfunc.cpp
+index 6bf10dcba..8a03e8744 100644
+--- a/fem/pgridfunc.cpp
++++ b/fem/pgridfunc.cpp
+@@ -870,12 +870,46 @@ void ParGridFunction::ProjectBdrCoefficientTangent(VectorCoefficient &vcoeff,
+       values(i) = values_counter[i] ? (*this)(i) : 0.0;
+    }
+ 
++   // Compute shared face DOF orientations for ND face DOF correction.
++   // For order > 1, shared triangular faces at partition boundaries have face
++   // DOFs that need orientation correction around Reduce/Bcast. This is
++   // analogous to ldof_sign for edge DOFs, but uses 2x2 matrices for face DOF
++   // pairs. T[ltori] maps master -> local; T_inv[ltori] maps local -> master.
++   const int ldof = Size();
++   Array<int> ltori, ldsize;
++   pfes->GetSharedTriFaceDofOrientations(ltori, ldsize);
++
++   auto apply_face_dof_transform = [&](real_t *data, bool inverse)
++   {
++      for (int i = 0; i < ldof; i++)
++      {
++         if (ldsize[i] == 2 && ltori[i] != 0)
++         {
++            const DenseMatrix &M = inverse
++                                   ? ND_DofTransformation::GetFaceInverseTransform(ltori[i])
++                                   : ND_DofTransformation::GetFaceTransform(ltori[i]);
++            MFEM_ASSERT(i+1 < ldof && ldsize[i+1] == 2,
++                        "face DOF pair not contiguous");
++            const real_t v0 = data[i], v1 = data[i+1];
++            data[i]   = M(0,0)*v0 + M(0,1)*v1;
++            data[i+1] = M(1,0)*v0 + M(1,1)*v1;
++            i++;
++         }
++      }
++   };
++
++   // Convert non-master shared face DOFs to master orientation before Reduce.
++   apply_face_dof_transform(values.HostReadWrite(), true);
++
+    // Count the values globally.
+    GroupCommunicator &gcomm = pfes->GroupComm();
+    gcomm.Reduce<int>(values_counter.HostReadWrite(), GroupCommunicator::Sum);
+    // Accumulate the values globally.
+    gcomm.Reduce<real_t>(values.HostReadWrite(), GroupCommunicator::Sum);
+ 
++   // Convert back to local orientation after Reduce.
++   apply_face_dof_transform(values.HostReadWrite(), false);
++
+    for (int i = 0; i < values.Size(); i++)
+    {
+       if (values_counter[i])
+@@ -883,9 +917,16 @@ void ParGridFunction::ProjectBdrCoefficientTangent(VectorCoefficient &vcoeff,
+          (*this)(i) = values(i)/values_counter[i];
+       }
+    }
++
++   // Convert to master orientation before Bcast.
++   apply_face_dof_transform((*this).HostReadWrite(), true);
++
+    // Broadcast values to other processors to have a consistent GridFunction
+    gcomm.Bcast<real_t>((*this).HostReadWrite());
+ 
++   // Convert back to local orientation after Bcast.
++   apply_face_dof_transform((*this).HostReadWrite(), false);
++
+ #ifdef MFEM_DEBUG
+    Array<int> ess_vdofs_marker;
+    pfes->GetEssentialVDofs(bdr_attr, ess_vdofs_marker);
diff --git a/extern/patch/mfem/patch_openmp_early_return.diff b/extern/patch/mfem/patch_openmp_early_return.diff
diff --git a/palace/models/spaceoperator.cpp b/palace/models/spaceoperator.cpp
@@ -3,7 +3,6 @@
 
 #include "spaceoperator.hpp"
 
-#include <limits>
 #include <set>
 #include <type_traits>
 #include "fem/bilinearform.hpp"
@@ -12,9 +11,6 @@
 #include "fem/mesh.hpp"
 #include "fem/multigrid.hpp"
 #include "linalg/hypre.hpp"
-#include "linalg/iterative.hpp"
-#include "linalg/jacobi.hpp"
-#include "linalg/ksp.hpp"
 #include "linalg/rap.hpp"
 #include "utils/communication.hpp"
 #include "utils/geodata.hpp"
@@ -549,93 +545,6 @@ auto BuildLevelParOperator<ComplexOperator>(std::unique_ptr<Operator> &&br,
   return std::make_unique<ComplexParOperator>(std::move(br), std::move(bi), fespace);
 }
 
-// Project a boundary coefficient into a Vector via a boundary mass matrix solve.
-//
-// This should be done by ParGridFunction::ProjectBdrCoefficientTangent, including parallel
-// reduction. See also https://github.com/mfem/mfem/pull/606. However there seems to be a
-// bug in MFEM that breaks this for Nédélec elements, perhaps due to orientation signs.
-// TODO(future): Investigate and fix this bug.
-//
-// Here project via CG solve of boundary mass matrix system M_bdr * e = f.
-void ProjectBdrCoefficientViaMassSolve(SumVectorCoefficient &fb, const LumpedPortData &data,
-                                       const MaterialOperator &mat_op,
-                                       FiniteElementSpace &nd_fespace, MPI_Comm comm,
-                                       Vector &result)
-{
-  // Assemble the boundary linear form f = ∫ φ_i · coeff dS (parallel-correct via P^T).
-  mfem::LinearForm rhs(&nd_fespace.Get());
-  rhs.AddBoundaryIntegrator(new VectorFEBoundaryLFIntegrator(fb));
-  rhs.UseFastAssembly(false);
-  rhs.UseDevice(false);
-  rhs.Assemble();
-  rhs.UseDevice(true);
-  nd_fespace.GetProlongationMatrix()->MultTranspose(rhs, result);
-
-  // Assemble boundary mass matrix M_bdr = ∫ φ_i · φ_j dS on port surface.
-  MaterialPropertyCoefficient fb_mass(mat_op.MaxCeedBdrAttribute());
-  for (const auto &elem : data.elems)
-  {
-    fb_mass.AddMaterialProperty(mat_op.GetCeedBdrAttributes(elem->GetAttrList()), 1.0);
-  }
-  BilinearForm m_bdr(nd_fespace);
-  if (!fb_mass.empty())
-  {
-    m_bdr.AddBoundaryIntegrator<VectorFEMassIntegrator>(fb_mass);
-  }
-  auto M_bdr = std::make_unique<ParOperator>(m_bdr.Assemble(false), nd_fespace);
-
-  // The boundary mass matrix is zero for all DOFs not on the port surface, making it
-  // singular on the full space. Set non-port DOFs as essential with DIAG_ONE so that M_bdr
-  // acts as identity there, giving a full-rank SPD system. non_port_tdof_list must outlive
-  // M_bdr because SetEssentialTrueDofs uses pointer not copy.
-  mfem::Array<int> attr_list;
-  for (const auto &elem : data.elems)
-  {
-    attr_list.Append(elem->GetAttrList());
-  }
-  const auto &mesh = nd_fespace.GetParMesh();
-  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, attr_marker);
-
-  mfem::Array<int> port_tdof_list;
-  nd_fespace.Get().GetEssentialTrueDofs(attr_marker, port_tdof_list);
-
-  mfem::Array<int> non_port_tdof_list;
-  {
-    std::vector<bool> is_port(nd_fespace.GetTrueVSize(), false);
-    for (int i = 0; i < port_tdof_list.Size(); i++)
-    {
-      is_port[port_tdof_list[i]] = true;
-    }
-    for (int i = 0; i < nd_fespace.GetTrueVSize(); i++)
-    {
-      if (!is_port[i])
-      {
-        non_port_tdof_list.Append(i);
-      }
-    }
-  }
-  M_bdr->SetEssentialTrueDofs(non_port_tdof_list, Operator::DIAG_ONE);
-
-  // CG solve M_bdr * e = f entirely in T-vector space.
-  // TODO: Make solver parameters configurable from IoData, or inherit other settings.
-  auto pcg = std::make_unique<CgSolver<Operator>>(comm, 0);
-  pcg->SetInitialGuess(false);
-  pcg->SetRelTol(1.0e-14);
-  pcg->SetAbsTol(std::numeric_limits<double>::epsilon());
-  pcg->SetMaxIter(200);
-  auto jac = std::make_unique<JacobiSmoother<Operator>>(comm);
-  auto ksp = std::make_unique<BaseKspSolver<Operator>>(std::move(pcg), std::move(jac));
-  ksp->SetOperators(*M_bdr, *M_bdr);
-
-  Vector sol(nd_fespace.GetTrueVSize());
-  sol.UseDevice(true);
-  sol = 0.0;
-  ksp->Mult(result, sol);
-  result = sol;
-}
-
 }  // namespace
 
 void SpaceOperator::AssemblePreconditioner(
@@ -953,33 +862,29 @@ void SpaceOperator::GetLumpedPortExcitationVectorPrimaryEt(int port_idx,
   const auto &data = GetLumpedPortOp().GetPort(port_idx);
 
   SumVectorCoefficient fb(GetMesh().SpaceDimension());
+  mfem::Array<int> attr_list;
   for (const auto &elem : data.elems)
   {
     const double Rs = 1.0 * data.GetToSquare(*elem);
     const double Einc = std::sqrt(
         Rs / (elem->GetGeometryWidth() * elem->GetGeometryLength() * data.elems.size()));
     fb.AddCoefficient(elem->GetModeCoefficient(Einc));
+    attr_list.Append(elem->GetAttrList());
   }
 
   Et_primary.SetSize(GetNDSpace().GetTrueVSize());
   Et_primary.UseDevice(true);
   Et_primary = 0.0;
 
-  // Broken code that should work using ParGridFunction::ProjectBdrCoefficientTangent.
-  // See ProjectBdrCoefficientViaMassSolve comment above.
-
-  //  mfem::Array<int> attr_marker;
-  //
-  // const auto &mesh = GetNDSpace().GetParMesh();
-  // int bdr_attr_max = mesh.bdr_attributes.Size() ? mesh.bdr_attributes.Max() : 0;
-  // mesh::AttrToMarker(bdr_attr_max, attr_list, attr_marker);
-  //
-  // GridFunction rhs(GetNDSpace());
-  // rhs = 0.0;
-  // rhs.Real().ProjectBdrCoefficientTangent(fb, attr_marker);
+  const auto &mesh = GetNDSpace().GetParMesh();
+  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, attr_marker);
 
-  ProjectBdrCoefficientViaMassSolve(fb, data, mat_op, GetNDSpace(), GetComm(),
-                                    Et_primary.Real());
+  GridFunction rhs(GetNDSpace());
+  rhs = 0.0;
+  rhs.Real().ProjectBdrCoefficientTangent(fb, attr_marker);
+  GetNDSpace().GetRestrictionMatrix()->Mult(rhs.Real(), Et_primary.Real());
 
   if (zero_metal)
   {
@@ -994,21 +899,29 @@ void SpaceOperator::GetLumpedPortExcitationVectorPrimaryHtcn(int port_idx,
   const auto &data = lumped_port_op.GetPort(port_idx);
 
   SumVectorCoefficient fb(GetMesh().SpaceDimension());
+  mfem::Array<int> attr_list;
   for (const auto &elem : data.elems)
   {
     const double Rs = 1.0 * data.GetToSquare(*elem);
     const double Hinc = 1.0 / std::sqrt(Rs * elem->GetGeometryWidth() *
                                         elem->GetGeometryLength() * data.elems.size());
     fb.AddCoefficient(elem->GetModeCoefficient(Hinc));
+    attr_list.Append(elem->GetAttrList());
   }
 
   Htcn_primary.SetSize(GetNDSpace().GetTrueVSize());
   Htcn_primary.UseDevice(true);
   Htcn_primary = 0.0;
 
-  // See ParGridFunction::ProjectBdrCoefficientTangent issue above.
-  ProjectBdrCoefficientViaMassSolve(fb, data, mat_op, GetNDSpace(), GetComm(),
-                                    Htcn_primary.Real());
+  const auto &mesh = GetNDSpace().GetParMesh();
+  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, attr_marker);
+
+  GridFunction rhs(GetNDSpace());
+  rhs = 0.0;
+  rhs.Real().ProjectBdrCoefficientTangent(fb, attr_marker);
+  GetNDSpace().GetRestrictionMatrix()->Mult(rhs.Real(), Htcn_primary.Real());
 
   if (zero_metal)
   {
diff --git a/spack_repo/local/packages/palace/mfem_pr5280.diff b/spack_repo/local/packages/palace/mfem_pr5280.diff
@@ -0,0 +1 @@
+../../../../extern/patch/mfem/mfem_pr5280.diff
diff --git a/spack_repo/local/packages/palace/package.py b/spack_repo/local/packages/palace/package.py
@@ -175,7 +175,7 @@ class Palace(CMakePackage, CudaPackage, ROCmPackage):
         # +lapack means: use external lapack
         depends_on(
             "mfem+mpi+metis+lapack@4.9:",
-            patches=["patch_par_tet_mesh_fix_dev.diff", "patch_gmsh_parser_performance.diff"],
+            patches=["patch_par_tet_mesh_fix_dev.diff", "patch_gmsh_parser_performance.diff", "mfem_pr5280.diff"],
         )
         depends_on("mfem+shared", when="+shared")
         depends_on("mfem~shared", when="~shared")

Original file line number	Diff line number	Diff line change
`@@ -398,6 +398,7 @@ message(STATUS "MFEM_OPTIONS: ${MFEM_OPTIONS_PRINT}")`
`398`	`398`	`set(MFEM_PATCH_FILES`
`399`	`399`	`"${CMAKE_SOURCE_DIR}/extern/patch/mfem/patch_par_tet_mesh_fix_dev.diff"`
`400`	`400`	`"${CMAKE_SOURCE_DIR}/extern/patch/mfem/patch_gmsh_parser_performance.diff"`
	`401`	`+ "${CMAKE_SOURCE_DIR}/extern/patch/mfem/mfem_pr5280.diff"`
`401`	`402`	`)`
`402`	`403`
`403`	`404`	`include(ExternalProject)`
Original file line number	Diff line number	Diff line change
`@@ -0,0 +1 @@`
	`1`	`+../../../../extern/patch/mfem/mfem_pr5280.diff`
Original file line number	Diff line number	Diff line change
`@@ -175,7 +175,7 @@ class Palace(CMakePackage, CudaPackage, ROCmPackage):`
`175`	`175`	`# +lapack means: use external lapack`
`176`	`176`	`depends_on(`
`177`	`177`	`"[email protected]:",`
`178`		`- patches=["patch_par_tet_mesh_fix_dev.diff", "patch_gmsh_parser_performance.diff"],`
	`178`	`+ patches=["patch_par_tet_mesh_fix_dev.diff", "patch_gmsh_parser_performance.diff", "mfem_pr5280.diff"],`
`179`	`179`	`)`
`180`	`180`	`depends_on("mfem+shared", when="+shared")`
`181`	`181`	`depends_on("mfem~shared", when="~shared")`