Add both dual and primary vector to PROM

- Add orthogonality tolerance in UpdatePROM
- Always use refined Gram-Schmidt in RomOp with circuit synthesis

Author: phdum-a

palace/models/romoperator.cpp

@@ -328,6 +328,12 @@ RomOperator::RomOperator(const IoData &iodata, SpaceOperator &space_op,
                          std::size_t max_size_per_excitation)
   : space_op(space_op), orthog_type(iodata.solver.linear.gs_orthog)
 {
+  // Always use refined CGS with adaptive synthesis for higher accuracy on print-out.
+  if (iodata.solver.driven.adaptive_circuit_synthesis)
+  {
+    orthog_type = Orthogonalization::CGS2;
+  }
+
   // Construct the system matrices defining the linear operator. PEC boundaries are handled
   // simply by setting diagonal entries of the system matrix for the corresponding dofs.
   // Because the Dirichlet BC is always homogeneous, no special elimination is required on
@@ -433,29 +439,60 @@ void RomOperator::SolveHDM(int excitation_idx, double omega, ComplexVector &u)
 
 void RomOperator::AddLumpedPortModesForSynthesis(const IoData &iodata)
 {
-  auto &lumped_port_op = space_op.GetLumpedPortOp();
-  for (const auto &[port_idx, port_data] : lumped_port_op)
+  // Add modes for lumped port to use them a circuit matrices.
+  //
+  // The excitation vector that we expect to add to the PROM is just the primary vector e
+  // associated with that port. However, when computing the response (scattering) matrix, we
+  // have an overlap g * e, where g is the boundary bilinear form. For Nédélec elements e
+  // and g * e are not proportional to each other. This is true on simple tet meshes, for
+  // simple hex meshes they are equal empirically, but this was not tested generally (e.g.
+  // on non-conforming meshes).
+  //
+  // To preserve orthogonality while also extracting the "full" port profile, we add both to
+  // the PROM — first g * e and then e (orthogonalized to g * e), which accounts for a
+  // "finite element correction". When e and g * e are the same, we should only add a single
+  // vector to the PROM.
+
+  // Workspace vector: Lumped Ports are Real, but interface is complex.
+  ComplexVector vec;
+  vec.SetSize(space_op.GetNDSpace().GetTrueVSize());
+  vec.UseDevice(true);
+  vec = 0.0;
+
+  // Add all dual vectors first: we want this port to be orthogonal to the outside world and
+  // check orthogonality internally.
+  for (const auto &[port_idx, port_data] : space_op.GetLumpedPortOp())
+  {
+    space_op.GetLumpedPortExcitationVectorDual(port_idx, vec, true);
+    UpdatePROM(vec, fmt::format("port_{:d}", port_idx));
+  }
+
+  // Check that the ports don't have any overlap.
+  MFEM_VERIFY(GetRomOrthogonalityMatrix().isDiagonal(),
+              "Lumped port fields on the mesh should have exactly zero overlap. This may "
+              "be non-zero if attributes share edges.");
+
+  // Add primary vector second, since this needs to be orthogonalized to the dual
+  // "outside" vector.
+  for (const auto &[port_idx, port_data] : space_op.GetLumpedPortOp())
   {
-    ComplexVector port_excitation_E;
-    port_excitation_E.UseDevice(true);
-    space_op.GetLumpedPortExcitationVector(port_idx, port_excitation_E, true);
-    UpdatePROM(port_excitation_E, fmt::format("port_{:d}", port_idx));
+    space_op.GetLumpedPortExcitationVectorPrimary(port_idx, vec, true);
+    // Drop if primary vector is same as dual vector.
+    auto orthog_condition = ORTHOG_TOL * linalg::Norml2(space_op.GetComm(), vec.Real());
+    UpdatePROM(vec, fmt::format("port_{:d}_correction", port_idx), orthog_condition);
   }
 
-  // Debug Print ROM Matrices
+  // Debug Print
   if constexpr (true)  // ignore
   {
     fs::path folder_tmp = fs::path(iodata.problem.output) / "prom_port_debug";
     fs::create_directories(folder_tmp);
     PrintPROMMatrices(iodata.units, folder_tmp);
   }
-
-  const auto &orth_R = GetRomOrthogonalityMatrix();
-  MFEM_VERIFY(orth_R.isDiagonal(), "Lumped port modes should have exactly zero overlap.");
 }
 
-
-void RomOperator::UpdatePROM(const ComplexVector &u, std::string_view node_label)
+void RomOperator::UpdatePROM(const ComplexVector &u, std::string_view node_label,
+                             double drop_degenerate_vector_norm_tol)
 {
   // Update PROM basis V. The basis is always real (each complex solution adds two basis
   // vectors, if it has a nonzero real and imaginary parts).
@@ -469,30 +506,44 @@ void RomOperator::UpdatePROM(const ComplexVector &u, std::string_view node_label
   const bool has_imag = (norm_im > norm_tol);
 
   const std::size_t dim_V_old = V.size();
-  const std::size_t dim_V_new =
+  std::size_t dim_V_new =
       V.size() + static_cast<std::size_t>(has_real) + static_cast<std::size_t>(has_imag);
 
   orth_R.conservativeResizeLike(Eigen::MatrixXd::Zero(dim_V_new, dim_V_new));
 
-  auto add_real_vector_to_basis = [this](const Vector &vector)
+  auto add_real_vector_to_basis = [this, drop_degenerate_vector_norm_tol](
+                                      const Vector &vector, std::string_view node_label)
   {
     auto dim_V = V.size();
     auto &v = V.emplace_back(vector);
     OrthogonalizeColumn(orthog_type, space_op.GetComm(), V, v, orth_R.col(dim_V).data(),
                         dim_V);
     orth_R(dim_V, dim_V) = linalg::Norml2(space_op.GetComm(), v);
+    if (orth_R(dim_V, dim_V) < drop_degenerate_vector_norm_tol)
+    {
+      V.pop_back();
+      return false;
+    }
     v *= 1.0 / orth_R(dim_V, dim_V);
+    v_node_label.emplace_back(node_label);
+    return true;
   };
 
   if (has_real)
   {
-    add_real_vector_to_basis(u.Real());
-    v_node_label.emplace_back(fmt::format("{}_re", node_label));
+    add_real_vector_to_basis(u.Real(), fmt::format("{}_re", node_label));
   }
   if (has_imag)
   {
-    add_real_vector_to_basis(u.Imag());
-    v_node_label.emplace_back(fmt::format("{}_im", node_label));
+    add_real_vector_to_basis(u.Imag(), fmt::format("{}_im", node_label));
+  }
+
+  // Vectors might have been dropped due to orthogonality check.
+  dim_V_new = V.size();
+  orth_R.conservativeResize(dim_V_new, dim_V_new);  // Might shrink
+  if (dim_V_new == dim_V_old)
+  {
+    return;
   }
 
   // Update reduced-order operators. Resize preserves the upper dim0 x dim0 block of each

palace/models/romoperator.hpp

@@ -134,7 +134,8 @@ class RomOperator
   // "node_label". This will be printed in the header of the csv files when printing PROM
   // matrices. It is needed to distinguish port and solution field configuration as well as
   // to reconstruct if field configuration are pure real, imaginary or complex.
-  void UpdatePROM(const ComplexVector &u, std::string_view node_label);
+  void UpdatePROM(const ComplexVector &u, std::string_view node_label,
+                  double drop_degenerate_vector_norm_tol = 0.0);
 
   // Add solution u to the minimal-rational interpolation for error estimation. MRI are
   // separated by excitation index.

palace/models/spaceoperator.cpp

@@ -841,43 +841,79 @@ bool SpaceOperator::GetExcitationVector(int excitation_idx, double omega,
   return nnz1 || nnz2;
 }
 
-bool SpaceOperator::GetLumpedPortExcitationVector(int port_idx, ComplexVector &RHS,
-                                                  bool zero_metal)
+bool SpaceOperator::GetLumpedPortExcitationVectorPrimary(int port_idx,
+                                                         ComplexVector &RHS_primary,
+                                                         bool zero_metal)
 {
-  RHS.SetSize(GetNDSpace().GetTrueVSize());
-  RHS.UseDevice(true);
-  RHS = 0.0;
+  const auto &data = lumped_port_op.GetPort(port_idx);
 
-  MFEM_VERIFY(RHS.Size() == GetNDSpace().GetTrueVSize(),
-              "Invalid T-vector size for AddExcitationVector1Internal!");
   SumVectorCoefficient fb(GetMesh().SpaceDimension());
+  mfem::Array<int> attr_list;
+  mfem::Array<int> attr_marker;
+  for (const auto &elem : data.elems)
+  {
+    attr_list.Append(elem->GetAttrList());
+    // TODO: Deal with effective reference normalization
+    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));
+  }
+  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);
+
+  RHS_primary.SetSize(GetNDSpace().GetTrueVSize());
+  RHS_primary.UseDevice(true);
+  RHS_primary = 0.0;
 
+  GridFunction rhs(GetNDSpace());
+  rhs = 0.0;
+  rhs.Real().ProjectBdrCoefficientTangent(fb, attr_marker);
+  GetNDSpace().GetProlongationMatrix()->AddMultTranspose(rhs.Real(), RHS_primary.Real());
+  if (zero_metal)
+  {
+    linalg::SetSubVector(RHS_primary.Real(), nd_dbc_tdof_lists.back(), 0.0);
+  }
+  return true;
+}
+
+bool SpaceOperator::GetLumpedPortExcitationVectorDual(int port_idx, ComplexVector &RHS_dual,
+                                                      bool zero_metal)
+{
   const auto &data = lumped_port_op.GetPort(port_idx);
 
+  SumVectorCoefficient fb(GetMesh().SpaceDimension());
   mfem::Array<int> attr_list;
   mfem::Array<int> attr_marker;
-
   for (const auto &elem : data.elems)
   {
     attr_list.Append(elem->GetAttrList());
-    fb.AddCoefficient(
-        elem->GetModeCoefficient(1.0 / (elem->GetGeometryWidth() * data.elems.size())));
+    // TODO: Deal with effective reference normalization
+    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));
   }
   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);
 
-  mfem::LinearForm rhs1(&GetNDSpace().Get());
-  rhs1.AddBoundaryIntegrator(new VectorFEBoundaryLFIntegrator(fb), attr_marker);
-  rhs1.UseFastAssembly(false);
-  rhs1.UseDevice(false);
-  rhs1.Assemble();
-  rhs1.UseDevice(true);
-  GetNDSpace().GetProlongationMatrix()->AddMultTranspose(rhs1, RHS.Real());
+  RHS_dual.SetSize(GetNDSpace().GetTrueVSize());
+  RHS_dual.UseDevice(true);
+  RHS_dual = 0.0;
+
+  mfem::LinearForm rhs(&GetNDSpace().Get());
+  rhs.AddBoundaryIntegrator(new VectorFEBoundaryLFIntegrator(fb), attr_marker);
+  rhs.UseFastAssembly(false);
+  rhs.UseDevice(false);
+  rhs.Assemble();
+  rhs.UseDevice(true);
+  GetNDSpace().GetProlongationMatrix()->Mult(rhs, RHS_dual.Real());
 
   if (zero_metal)
   {
-    linalg::SetSubVector(RHS.Real(), nd_dbc_tdof_lists.back(), 0.0);
+    linalg::SetSubVector(RHS_dual.Real(), nd_dbc_tdof_lists.back(), 0.0);
   }
   return true;
 }
@@ -934,8 +970,8 @@ bool SpaceOperator::AddExcitationVector1Internal(int excitation_idx, Vector &RHS
 bool SpaceOperator::AddExcitationVector2Internal(int excitation_idx, double omega,
                                                  ComplexVector &RHS2)
 {
-  // Assemble the contribution of wave ports to the frequency domain excitation term at the
-  // specified frequency.
+  // Assemble the contribution of wave ports to the frequency domain excitation term at
+  // the specified frequency.
   MFEM_VERIFY(RHS2.Size() == GetNDSpace().GetTrueVSize(),
               "Invalid T-vector size for AddExcitationVector2Internal!");
   SumVectorCoefficient fbr(GetMesh().SpaceDimension()), fbi(GetMesh().SpaceDimension());

palace/models/spaceoperator.hpp

@@ -215,13 +215,16 @@ class SpaceOperator
   bool GetExcitationVector(int excitation_idx, Vector &RHS);
   bool GetExcitationVector(int excitation_idx, double omega, ComplexVector &RHS);
 
-  bool GetLumpedPortExcitationVector(int port_idx, ComplexVector &RHS, bool zero_metal);
-
   // Separate out RHS vector as RHS = iω RHS1 + RHS2(ω). The return value indicates whether
   // or not the excitation is nonzero (and thus is true most of the time).
   bool GetExcitationVector1(int excitation_idx, ComplexVector &RHS1);
   bool GetExcitationVector2(int excitation_idx, double omega, ComplexVector &RHS2);
 
+  bool GetLumpedPortExcitationVectorPrimary(int port_idx, ComplexVector &RHS_primary,
+                                            bool zero_metal);
+  bool GetLumpedPortExcitationVectorDual(int port_idx, ComplexVector &RHS_dual,
+                                         bool zero_metal);
+
   // Construct a constant or randomly initialized vector which satisfies the PEC essential
   // boundary conditions.
   void GetRandomInitialVector(ComplexVector &v);

test/unit/test-drivensolver.cpp

@@ -63,5 +63,5 @@ TEST_CASE("PortOrthogonalityAdaptivePROM", "[driven_solver][Serial][Parallel]")
                       utils::ConfigureOmp(), ""};
   CHECK_THROWS(solver.SolveEstimateMarkRefine(mesh_),
                Catch::Matchers::ContainsSubstring(
-                   "Lumped port modes should have exactly zero overlap"));
+                   "should have exactly zero overlap"));
 }