diff --git a/CHANGELOG.md b/CHANGELOG.md index 1a2e7d157..88b1eec9f 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -15,6 +15,11 @@ The format of this changelog is based on - Fixed a small regression bug for boundary postprocessing when specifying `"Side": "LargerRefractiveIndex"`, introduced as part of v0.13.0. + - Added an improvement to numeric wave ports to avoid targetting evanescent modes at + higher operating frequencies. Also finite conductivity boundaries + (`config["Boundaries"]["Conductivity"]`) are automatically marked as PEC for the wave + port mode solve (previously these were marked as PMC unless specified under + `"WavePortPEC"`). ## [0.13.0] - 2024-05-20 @@ -61,7 +66,7 @@ The format of this changelog is based on configuration file keyword changes to for consistency were made to `config["Domains"]["Postprocessing"]["Probe"]` and `config["Model"]["Refinement"]["Boxes"]`. - - Fix a bug in MFEM for nonconformal AMR meshes with internal boundaries affecting + - Fixed a bug in MFEM for nonconformal AMR meshes with internal boundaries affecting non-homogeneous Dirichlet boundary conditions for electrostatic simulations (see [#236](https://github.com/awslabs/palace/pull/236)). diff --git a/docs/src/config/boundaries.md b/docs/src/config/boundaries.md index f595bc86a..a202f2c11 100644 --- a/docs/src/config/boundaries.md +++ b/docs/src/config/boundaries.md @@ -355,7 +355,8 @@ corresponding coordinate system. "Excitation": , "Active": , "Mode": , - "Offset": + "Offset": , + "SolverType": }, ... ] @@ -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 diff --git a/docs/src/config/solver.md b/docs/src/config/solver.md index e079818bf..1d4889eab 100644 --- a/docs/src/config/solver.md +++ b/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]` diff --git a/palace/drivers/eigensolver.cpp b/palace/drivers/eigensolver.cpp index 42320fab3..8f12cba7e 100644 --- a/palace/drivers/eigensolver.cpp +++ b/palace/drivers/eigensolver.cpp @@ -60,10 +60,6 @@ EigenSolver::Solve(const std::vector> &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> &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"); diff --git a/palace/models/waveportoperator.cpp b/palace/models/waveportoperator.cpp index e18ca8772..250c5dd4d 100644 --- a/palace/models/waveportoperator.cpp +++ b/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(std::move(Br), std::move(Bi)); + opB = std::make_unique(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(port_comm, print); #endif } - else // config::EigenSolverData::Type::SLEPC + else // config::WavePortData::EigenSolverType::SLEPC { #if defined(PALACE_WITH_SLEPC) auto slepc = std::make_unique(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 A; + // port mode. The B matrix is operating frequency-independent and has already been + // constructed. + std::unique_ptr 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(std::move(Ar), std::move(Ai)); + opA = std::make_unique(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 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!"); @@ -1092,7 +1097,8 @@ void WavePortOperator::SetUpBoundaryProperties(const IoData &iodata, // are touching and share one or more edges. mfem::Array dbc_bcs; dbc_bcs.Reserve(static_cast(iodata.boundaries.pec.attributes.size() + - iodata.boundaries.auxpec.attributes.size())); + iodata.boundaries.auxpec.attributes.size() + + iodata.boundaries.conductivity.size())); for (auto attr : iodata.boundaries.pec.attributes) { if (attr <= 0 || attr > bdr_attr_max) @@ -1109,6 +1115,17 @@ void WavePortOperator::SetUpBoundaryProperties(const IoData &iodata, } dbc_bcs.Append(attr); } + for (const auto &data : iodata.boundaries.conductivity) + { + for (auto attr : data.attributes) + { + if (attr <= 0 || attr > bdr_attr_max) + { + continue; // Can just ignore if wrong + } + dbc_bcs.Append(attr); + } + } // If user accidentally specifies a surface as both "PEC" and "WavePortPEC", this is fine // so allow for duplicates in the attribute list. dbc_bcs.Sort(); diff --git a/palace/models/waveportoperator.hpp b/palace/models/waveportoperator.hpp index 3b703e97a..e581fcf5c 100644 --- a/palace/models/waveportoperator.hpp +++ b/palace/models/waveportoperator.hpp @@ -66,7 +66,7 @@ class WavePortData // Operator storage for repeated boundary mode eigenvalue problem solves. std::unique_ptr Atnr, Atni, Antr, Anti, Annr, Anni; - std::unique_ptr B0; + std::unique_ptr opB; ComplexVector v0, e0; // Eigenvalue solver for boundary modes. diff --git a/palace/utils/configfile.cpp b/palace/utils/configfile.cpp index 6b3c3155d..6b1fbc32c 100644 --- a/palace/utils/configfile.cpp +++ b/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'; diff --git a/palace/utils/configfile.hpp b/palace/utils/configfile.hpp index 550165978..0014dc8ef 100644 --- a/palace/utils/configfile.hpp +++ b/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); }; diff --git a/palace/utils/iodata.cpp b/palace/utils/iodata.cpp index ef48debd1..4409abf69 100644 --- a/palace/utils/iodata.cpp +++ b/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; diff --git a/scripts/schema/config/boundaries.json b/scripts/schema/config/boundaries.json index 6b86c4ede..f14b08fb1 100644 --- a/scripts/schema/config/boundaries.json +++ b/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" } }