Update to wave ports to exclude evanescent modes from boundary mode c…

…alculation

Also removes old FEAST code from docs, simplifies ability to specify eigenmode solver to wave ports.

docs/src/config/boundaries.md

@@ -355,7 +355,8 @@ corresponding coordinate system.
         "Excitation": <bool>,
         "Active": <bool>,
         "Mode": <int>,
-        "Offset": <float>
+        "Offset": <float>,
+        "SolverType": <string>
     },
     ...
 ]
@@ -380,6 +381,10 @@ boundary. Ranked in order of decreasing wave number.
 `"Offset" [0.0]` :  Offset distance used for scattering parameter de-embedding for this wave
 port boundary, specified in mesh length units.
 
+`"SolverType" ["Default"]` :  Specifies the eigenvalue solver to be used in computing
+the boundary mode for this wave port. See
+[`config["Solver"]["Eigenmode"]["Type"]`](solver.md#solver%5B%22Eigenmode%22%5D).
+
 ## `boundaries["WavePortPEC"]`
 
 ```json

docs/src/config/solver.md

@@ -138,21 +138,9 @@ number of eigenmodes of the problem. The available options are:
 
   - `"SLEPc"`
   - `"ARPACK"`
-  - `"FEAST"`
   - `"Default"` :  Use the default eigensolver. Currently, this is the Krylov-Schur
     eigenvalue solver from `"SLEPc"`.
 
-`"ContourTargetUpper" [None]` :  Specifies the upper frequency target of the contour used
-for the FEAST eigenvalue solver, GHz. This option is relevant only for `"Type": "FEAST"`.
-
-`"ContourAspectRatio" [None]` :  Specifies the aspect ratio of the contour used for the
-FEAST eigenvalue solver. This should be greater than zero, where the aspect ratio is the
-ratio of the contour width to the frequency range(`"ContourTargetUpper"` - `"Target"`).
-This option is relevant only for `"Type": "FEAST"`.
-
-`"ContourNPoints" [4]` :  Number of contour integration points used for the FEAST eigenvalue
-solver. This option is relevant only for `"Type": "FEAST"`.
-
 ### Advanced eigenmode solver options
 
   - `"PEPLinear" [true]`

palace/drivers/eigensolver.cpp

@@ -60,10 +60,6 @@ EigenSolver::Solve(const std::vector<std::unique_ptr<Mesh>> &mesh) const
   {
     Mpi::Warning("SLEPc eigensolver not available, using ARPACK!\n");
   }
-  else if (iodata.solver.eigenmode.type == config::EigenSolverData::Type::FEAST)
-  {
-    Mpi::Warning("FEAST eigensolver requires SLEPc, using ARPACK!\n");
-  }
   type = config::EigenSolverData::Type::ARPACK;
 #elif defined(PALACE_WITH_SLEPC)
   if (iodata.solver.eigenmode.type == config::EigenSolverData::Type::ARPACK)
@@ -74,11 +70,7 @@ EigenSolver::Solve(const std::vector<std::unique_ptr<Mesh>> &mesh) const
 #else
 #error "Eigenmode solver requires building with ARPACK or SLEPc!"
 #endif
-  if (type == config::EigenSolverData::Type::FEAST)
-  {
-    MFEM_ABORT("FEAST eigenvalue solver is currently not supported!");
-  }
-  else if (type == config::EigenSolverData::Type::ARPACK)
+  if (type == config::EigenSolverData::Type::ARPACK)
   {
 #if defined(PALACE_WITH_ARPACK)
     Mpi::Print("\nConfiguring ARPACK eigenvalue solver:\n");

palace/models/waveportoperator.cpp

@@ -599,12 +599,13 @@ WavePortData::WavePortData(const config::WavePortData &data,
   //            given frequency, Math. Comput. (2003).
   // See also: Halla and Monk, On the analysis of waveguide modes in an electromagnetic
   //           transmission line, arXiv:2302.11994 (2023).
-  mu_eps_min = 0.0;  // Use standard inverse transformation to avoid conditioning issues
-                     // associated with shift
-  // const double c_max = mat_op.GetLightSpeedMax().Max();
-  // MFEM_VERIFY(c_max > 0.0 && c_max < mfem::infinity(),
-  //             "Invalid material speed of light detected in WavePortOperator!");
-  // mu_eps_min = 1.0 / (c_max * c_max);
+  const double c_max = mat_op.GetLightSpeedMax().Max();
+  MFEM_VERIFY(c_max > 0.0 && c_max < mfem::infinity(),
+              "Invalid material speed of light detected in WavePortOperator!");
+  mu_eps_min = 1.0 / (c_max * c_max) * 0.5;  // Add a safety factor for minimum propagation
+                                             // constant possible
+  // mu_eps_min = 0.0;  // Use standard inverse transformation to avoid conditioning issues
+  //                    // associated with shift
   std::tie(Atnr, Atni) = GetAtn(mat_op, *port_nd_fespace, *port_h1_fespace);
   std::tie(Antr, Anti) = GetAnt(mat_op, *port_h1_fespace, *port_nd_fespace);
   std::tie(Annr, Anni) = GetAnn(mat_op, *port_h1_fespace, port_normal);
@@ -621,7 +622,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);
-    B0 = std::make_unique<ComplexWrapperOperator>(std::move(Br), std::move(Bi));
+    opB = std::make_unique<ComplexWrapperOperator>(std::move(Br), std::move(Bi));
   }
 
   // Configure a communicator for the processes which have elements for this port.
@@ -729,15 +730,15 @@ WavePortData::WavePortData(const config::WavePortData &data,
 
     // Define the eigenvalue solver.
     constexpr int print = 0;
-    config::EigenSolverData::Type type = solver.eigenmode.type;
-    if (type == config::EigenSolverData::Type::SLEPC)
+    config::WavePortData::EigenSolverType type = data.eigen_type;
+    if (type == config::WavePortData::EigenSolverType::SLEPC)
     {
 #if !defined(PALACE_WITH_SLEPC)
       MFEM_ABORT("Solver was not built with SLEPc support, please choose a "
                  "different solver!");
 #endif
     }
-    else if (type == config::EigenSolverData::Type::ARPACK)
+    else if (type == config::WavePortData::EigenSolverType::ARPACK)
     {
 #if !defined(PALACE_WITH_ARPACK)
       MFEM_ABORT("Solver was not built with ARPACK support, please choose a "
@@ -747,20 +748,20 @@ WavePortData::WavePortData(const config::WavePortData &data,
     else  // Default choice
     {
 #if defined(PALACE_WITH_SLEPC)
-      type = config::EigenSolverData::Type::SLEPC;
+      type = config::WavePortData::EigenSolverType::SLEPC;
 #elif defined(PALACE_WITH_ARPACK)
-      type = config::EigenSolverData::Type::ARPACK;
+      type = config::WavePortData::EigenSolverType::ARPACK;
 #else
 #error "Wave port solver requires building with ARPACK or SLEPc!"
 #endif
     }
-    if (type == config::EigenSolverData::Type::ARPACK)
+    if (type == config::WavePortData::EigenSolverType::ARPACK)
     {
 #if defined(PALACE_WITH_ARPACK)
       eigen = std::make_unique<arpack::ArpackEPSSolver>(port_comm, print);
 #endif
     }
-    else  // config::EigenSolverData::Type::SLEPC
+    else  // config::WavePortData::EigenSolverType::SLEPC
     {
 #if defined(PALACE_WITH_SLEPC)
       auto slepc = std::make_unique<slepc::SlepcEPSSolver>(port_comm, print);
@@ -772,8 +773,12 @@ WavePortData::WavePortData(const config::WavePortData &data,
     constexpr double tol = 1.0e-6;
     eigen->SetNumModes(mode_idx, std::max(2 * mode_idx + 1, 5));
     eigen->SetTol(tol);
-    eigen->SetWhichEigenpairs(EigenvalueSolver::WhichType::LARGEST_MAGNITUDE);
     eigen->SetLinearSolver(*ksp);
+
+    // We want to ignore evanescent modes (kₙ with large imaginary component). The
+    // eigenvalue 1 / (-kₙ² - σ) of the shifted problem will be a large-magnitude negative
+    // real number for an eigenvalue kₙ² with real part close to but not below the cutoff σ.
+    eigen->SetWhichEigenpairs(EigenvalueSolver::WhichType::SMALLEST_REAL);
   }
 
   // Configure port mode sign convention: 1ᵀ Re{-n x H} >= 0 on the "upper-right quadrant"
@@ -844,16 +849,16 @@ void WavePortData::Initialize(double omega)
   }
 
   // Construct matrices and solve the generalized eigenvalue problem for the desired wave
-  // port mode. B uses the non-owning constructor since the matrices Br, Bi are not
-  // functions of frequency (constructed once for all).
-  std::unique_ptr<ComplexOperator> A;
+  // port mode. The B matrix is operating frequency-independent and has already been
+  // constructed.
+  std::unique_ptr<ComplexOperator> opA;
   const double sigma = -omega * omega * mu_eps_min;
   {
     auto [Attr, Atti] = GetAtt(mat_op, *port_nd_fespace, port_normal, omega, sigma);
     auto [Ar, Ai] =
         GetSystemMatrixA(Attr.get(), Atti.get(), Atnr.get(), Atni.get(), Antr.get(),
                          Anti.get(), Annr.get(), Anni.get(), port_dbc_tdof_list);
-    A = std::make_unique<ComplexWrapperOperator>(std::move(Ar), std::move(Ai));
+    opA = std::make_unique<ComplexWrapperOperator>(std::move(Ar), std::move(Ai));
   }
 
   // Configure and solve the (inverse) eigenvalue problem for the desired boundary mode.
@@ -862,9 +867,9 @@ void WavePortData::Initialize(double omega)
   std::complex<double> lambda;
   if (port_comm != MPI_COMM_NULL)
   {
-    ComplexWrapperOperator P(A->Real(), nullptr);  // Non-owning constructor
-    ksp->SetOperators(*A, P);
-    eigen->SetOperators(*B0, *A, EigenvalueSolver::ScaleType::NONE);
+    ComplexWrapperOperator opP(opA->Real(), nullptr);  // Non-owning constructor
+    ksp->SetOperators(*opA, opP);
+    eigen->SetOperators(*opB, *opA, EigenvalueSolver::ScaleType::NONE);
     eigen->SetInitialSpace(v0);
     int num_conv = eigen->Solve();
     MFEM_VERIFY(num_conv >= mode_idx, "Wave port eigensolver did not converge!");

palace/models/waveportoperator.hpp

@@ -66,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> B0;
+  std::unique_ptr<ComplexOperator> opB;
   ComplexVector v0, e0;
 
   // Eigenvalue solver for boundary modes.

palace/utils/configfile.cpp

@@ -1008,6 +1008,12 @@ void LumpedPortBoundaryData::SetUp(json &boundaries)
   }
 }
 
+// Helper for converting string keys to enum for WavePortData::EigenSolverType.
+PALACE_JSON_SERIALIZE_ENUM(WavePortData::EigenSolverType,
+                           {{WavePortData::EigenSolverType::DEFAULT, "Default"},
+                            {WavePortData::EigenSolverType::SLEPC, "SLEPc"},
+                            {WavePortData::EigenSolverType::ARPACK, "ARPACK"}})
+
 void WavePortBoundaryData::SetUp(json &boundaries)
 {
   auto port = boundaries.find("WavePort");
@@ -1034,6 +1040,7 @@ void WavePortBoundaryData::SetUp(json &boundaries)
     MFEM_VERIFY(data.mode_idx > 0,
                 "\"WavePort\" boundary \"Mode\" must be positive (1-based)!");
     data.d_offset = it->value("Offset", data.d_offset);
+    data.eigen_type = it->value("SolverType", data.eigen_type);
     data.excitation = it->value("Excitation", data.excitation);
     data.active = it->value("Active", data.active);
 
@@ -1042,6 +1049,7 @@ void WavePortBoundaryData::SetUp(json &boundaries)
     it->erase("Attributes");
     it->erase("Mode");
     it->erase("Offset");
+    it->erase("SolverType");
     it->erase("Excitation");
     it->erase("Active");
     MFEM_VERIFY(it->empty(),
@@ -1055,6 +1063,7 @@ void WavePortBoundaryData::SetUp(json &boundaries)
       std::cout << "Attributes: " << data.attributes << '\n';
       std::cout << "Mode: " << data.mode_idx << '\n';
       std::cout << "Offset: " << data.d_offset << '\n';
+      std::cout << "SolverType: " << data.eigen_type << '\n';
       std::cout << "Excitation: " << data.excitation << '\n';
       std::cout << "Active: " << data.active << '\n';
     }
@@ -1426,13 +1435,6 @@ void DrivenSolverData::SetUp(json &solver)
   }
 }
 
-// Helper for converting string keys to enum for EigenSolverData::Type.
-PALACE_JSON_SERIALIZE_ENUM(EigenSolverData::Type,
-                           {{EigenSolverData::Type::DEFAULT, "Default"},
-                            {EigenSolverData::Type::SLEPC, "SLEPc"},
-                            {EigenSolverData::Type::ARPACK, "ARPACK"},
-                            {EigenSolverData::Type::FEAST, "FEAST"}})
-
 void EigenSolverData::SetUp(json &solver)
 {
   auto eigenmode = solver.find("Eigenmode");
@@ -1451,17 +1453,6 @@ void EigenSolverData::SetUp(json &solver)
   n_post = eigenmode->value("Save", n_post);
   type = eigenmode->value("Type", type);
   pep_linear = eigenmode->value("PEPLinear", pep_linear);
-  feast_contour_np = eigenmode->value("ContourNPoints", feast_contour_np);
-  if (type == EigenSolverData::Type::FEAST && feast_contour_np > 1)
-  {
-    MFEM_VERIFY(eigenmode->find("ContourTargetUpper") != eigenmode->end() &&
-                    eigenmode->find("ContourAspectRatio") != eigenmode->end(),
-                "Missing \"Eigenmode\" solver \"ContourTargetUpper\" or "
-                "\"ContourAspectRatio\" for FEAST solver in the configuration file!");
-  }
-  feast_contour_ub = eigenmode->value("ContourTargetUpper", feast_contour_ub);
-  feast_contour_ar = eigenmode->value("ContourAspectRatio", feast_contour_ar);
-  feast_moments = eigenmode->value("ContourMoments", feast_moments);
   scale = eigenmode->value("Scaling", scale);
   init_v0 = eigenmode->value("StartVector", init_v0);
   init_v0_const = eigenmode->value("StartVectorConstant", init_v0_const);
@@ -1476,10 +1467,6 @@ void EigenSolverData::SetUp(json &solver)
   eigenmode->erase("Save");
   eigenmode->erase("Type");
   eigenmode->erase("PEPLinear");
-  eigenmode->erase("ContourNPoints");
-  eigenmode->erase("ContourTargetUpper");
-  eigenmode->erase("ContourAspectRatio");
-  eigenmode->erase("ContourMoments");
   eigenmode->erase("Scaling");
   eigenmode->erase("StartVector");
   eigenmode->erase("StartVectorConstant");
@@ -1499,10 +1486,6 @@ void EigenSolverData::SetUp(json &solver)
     std::cout << "Save: " << n_post << '\n';
     std::cout << "Type: " << type << '\n';
     std::cout << "PEPLinear: " << pep_linear << '\n';
-    std::cout << "ContourNPoints: " << feast_contour_np << '\n';
-    std::cout << "ContourTargetUpper: " << feast_contour_ub << '\n';
-    std::cout << "ContourAspectRatio: " << feast_contour_ar << '\n';
-    std::cout << "ContourMoments: " << feast_moments << '\n';
     std::cout << "Scaling: " << scale << '\n';
     std::cout << "StartVector: " << init_v0 << '\n';
     std::cout << "StartVectorConstant: " << init_v0_const << '\n';

palace/utils/configfile.hpp

@@ -443,6 +443,15 @@ struct WavePortData
   // Port offset for de-embedding [m].
   double d_offset = 0.0;
 
+  // Eigenvalue solver type for boundary mode calculation.
+  enum class EigenSolverType
+  {
+    DEFAULT,
+    SLEPC,
+    ARPACK
+  };
+  EigenSolverType eigen_type = EigenSolverType::DEFAULT;
+
   // Flag for source term in driven and transient simulations.
   bool excitation = false;
 
@@ -647,29 +656,13 @@ struct EigenSolverData
   bool mass_orthog = false;
 
   // Eigenvalue solver type.
-  enum class Type
-  {
-    DEFAULT,
-    SLEPC,
-    ARPACK,
-    FEAST
-  };
+  using Type = WavePortData::EigenSolverType;
   Type type = Type::DEFAULT;
 
   // For SLEPc eigenvalue solver, use linearized formulation for quadratic eigenvalue
   // problems.
   bool pep_linear = true;
 
-  // Number of integration points used for the FEAST eigenvalue solver contour.
-  int feast_contour_np = 4;
-
-  // Parameters for the FEAST eigenvalue solver contour.
-  double feast_contour_ub = 0.0;
-  double feast_contour_ar = 1.0;
-
-  // Use more than just the standard single moment for FEAST subspace construction.
-  int feast_moments = 1;
-
   void SetUp(json &solver);
 };
 

palace/utils/iodata.cpp

@@ -549,7 +549,6 @@ void IoData::NondimensionalizeInputs(mfem::ParMesh &mesh)
 
   // For eigenmode simulations:
   solver.eigenmode.target *= 2.0 * M_PI * tc;
-  solver.eigenmode.feast_contour_ub *= 2.0 * M_PI * tc;
 
   // For driven simulations:
   solver.driven.min_f *= 2.0 * M_PI * tc;

scripts/schema/config/boundaries.json

@@ -132,6 +132,7 @@
           "Attributes": { "$ref": "#/$defs/Attributes" },
           "Mode": { "type": "integer", "exclusiveMinimum": 0 },
           "Offset": { "type": "number", "minimum": 0.0 },
+          "SolverType": { "type": "string" },
           "Excitation": { "type": "boolean" },
           "Active": { "type": "boolean" }
         }