Merge pull request #192 from awslabs/sjg/wave-ports-dev

Fix bug for wave ports

CHANGELOG.md

@@ -20,6 +20,7 @@ The format of this changelog is based on
     domain electric field energy.
   - Added documentation for various timer categories and improved timing breakdown of
     various sections of a simulation.
+  - Fixed bug in implementation of numeric wave ports for driven simulations.
 
 ## [0.12.0] - 2023-12-21
 

palace/models/curlcurloperator.cpp

@@ -36,11 +36,6 @@ CurlCurlOperator::CurlCurlOperator(const IoData &iodata,
     mat_op(iodata, *mesh.back()), surf_j_op(iodata, GetH1Space())
 {
   // Finalize setup.
-  BilinearForm::pa_order_threshold = iodata.solver.pa_order_threshold;
-  fem::DefaultIntegrationOrder::p_trial = iodata.solver.order;
-  fem::DefaultIntegrationOrder::q_order_jac = iodata.solver.q_order_jac;
-  fem::DefaultIntegrationOrder::q_order_extra_pk = iodata.solver.q_order_extra;
-  fem::DefaultIntegrationOrder::q_order_extra_qk = iodata.solver.q_order_extra;
   CheckBoundaryProperties();
 
   // Print essential BC information.

palace/models/laplaceoperator.cpp

@@ -29,13 +29,6 @@ LaplaceOperator::LaplaceOperator(const IoData &iodata,
     nd_fespace(h1_fespaces.GetFinestFESpace(), *mesh.back(), nd_fec.get()),
     mat_op(iodata, *mesh.back()), source_attr_lists(ConstructSources(iodata))
 {
-  // Finalize setup.
-  BilinearForm::pa_order_threshold = iodata.solver.pa_order_threshold;
-  fem::DefaultIntegrationOrder::p_trial = iodata.solver.order;
-  fem::DefaultIntegrationOrder::q_order_jac = iodata.solver.q_order_jac;
-  fem::DefaultIntegrationOrder::q_order_extra_pk = iodata.solver.q_order_extra;
-  fem::DefaultIntegrationOrder::q_order_extra_qk = iodata.solver.q_order_extra;
-
   // Print essential BC information.
   if (dbc_attr.Size())
   {

palace/models/postoperator.cpp

@@ -493,8 +493,8 @@ std::complex<double> PostOperator::GetSParameter(const WavePortOperator &wave_po
   {
     Sij.real(Sij.real() - 1.0);
   }
-  // Port de-embedding: S_demb = S exp(-ikₙᵢ dᵢ) exp(-ikₙⱼ dⱼ) (distance offset is default
-  // 0 unless specified).
+  // Port de-embedding: S_demb = S exp(ikₙᵢ dᵢ) exp(ikₙⱼ dⱼ) (distance offset is default 0
+  // unless specified).
   Sij *= std::exp(1i * src_data.kn0 * src_data.d_offset);
   Sij *= std::exp(1i * data.kn0 * data.d_offset);
   return Sij;

palace/models/spaceoperator.cpp

@@ -45,11 +45,6 @@ SpaceOperator::SpaceOperator(const IoData &iodata,
     surf_j_op(iodata, GetH1Space())
 {
   // Finalize setup.
-  BilinearForm::pa_order_threshold = iodata.solver.pa_order_threshold;
-  fem::DefaultIntegrationOrder::p_trial = iodata.solver.order;
-  fem::DefaultIntegrationOrder::q_order_jac = iodata.solver.q_order_jac;
-  fem::DefaultIntegrationOrder::q_order_extra_pk = iodata.solver.q_order_extra;
-  fem::DefaultIntegrationOrder::q_order_extra_qk = iodata.solver.q_order_extra;
   CheckBoundaryProperties();
 
   // Print essential BC information.

palace/models/waveportoperator.hpp

@@ -50,22 +50,22 @@ class WavePortData
   mfem::Vector port_normal;
 
 private:
+  // List of all boundary attributes making up this port boundary.
+  mfem::Array<int> attr_list;
+
   // SubMesh data structures to define finite element spaces and grid functions on the
   // SubMesh corresponding to this port boundary.
   std::unique_ptr<Mesh> port_mesh;
   std::unique_ptr<mfem::FiniteElementCollection> port_nd_fec, port_h1_fec;
   std::unique_ptr<FiniteElementSpace> port_nd_fespace, port_h1_fespace;
   std::unique_ptr<mfem::ParTransferMap> port_nd_transfer, port_h1_transfer;
   std::unordered_map<int, int> submesh_parent_elems;
-
-  // List of all boundary attributes making up this port boundary.
-  mfem::Array<int> attr_list;
+  mfem::Array<int> port_dbc_tdof_list;
+  double mu_eps_min;
 
   // Operator storage for repeated boundary mode eigenvalue problem solves.
-  double mu_eps_max;
-  std::unique_ptr<mfem::HypreParMatrix> Btt, Btn, BtnT, Dtt;
-  ComplexVector v0, e0, e0t, e0n;
-  mfem::Array<int> port_dbc_tdof_list;
+  std::unique_ptr<mfem::HypreParMatrix> Atnr, Atni, Antr, Anti, Annr, Anni, Br, Bi;
+  ComplexVector v0, e0, e0n;
 
   // Eigenvalue solver for boundary modes.
   MPI_Comm port_comm;

palace/utils/iodata.cpp

@@ -14,6 +14,8 @@
 #include <string_view>
 #include <mfem.hpp>
 #include <nlohmann/json.hpp>
+#include "fem/bilinearform.hpp"
+#include "fem/integrator.hpp"
 #include "utils/communication.hpp"
 #include "utils/constants.hpp"
 #include "utils/geodata.hpp"
@@ -417,6 +419,13 @@ void IoData::CheckConfiguration()
   {
     solver.linear.mg_smooth_order = std::max(2 * solver.order, 4);
   }
+
+  // Configure settings for quadrature rules and partial assembly.
+  BilinearForm::pa_order_threshold = solver.pa_order_threshold;
+  fem::DefaultIntegrationOrder::p_trial = solver.order;
+  fem::DefaultIntegrationOrder::q_order_jac = solver.q_order_jac;
+  fem::DefaultIntegrationOrder::q_order_extra_pk = solver.q_order_extra;
+  fem::DefaultIntegrationOrder::q_order_extra_qk = solver.q_order_extra;
 }
 
 namespace
@@ -437,7 +446,7 @@ constexpr config::SymmetricMatrixData<N> &operator/=(config::SymmetricMatrixData
 
 void IoData::NondimensionalizeInputs(mfem::ParMesh &mesh)
 {
-  // Nondimensionalization of the equations is based on a given length Lc in[m], typically
+  // Nondimensionalization of the equations is based on a given length Lc in [m], typically
   // the largest domain dimension. Configuration file lengths and the mesh coordinates are
   // provided with units of model.L0 x [m].
   MFEM_VERIFY(!init, "NondimensionalizeInputs should only be called once!");