PROM Circuit Extraction

CHANGELOG.md

@@ -36,6 +36,8 @@ The format of this changelog is based on
   - Added an option to drop small entries (below machine epsilon) from the matrix used in the sparse
     direct solver. This can be specified with `config["Solver"]["Linear"]["DropSmallEntries"]`.
     [PR 476](https://github.com/awslabs/palace/pull/476).
+  - Added support for circuit synthesis from the PROM in adaptive driven simulations
+    (`config["Solver"]["Driven"]["AdaptiveCircuitSynthesis"]`). Modifies ROM to add `LumpedPort` fields to make circuit. Prints out the projected matrices of the driven problem as well as the basis orthogonalization matrix.
 
 #### Interface Changes
 

docs/src/config/boundaries.md

@@ -265,7 +265,7 @@ accounts for nonzero metal thickness.
         "Elements":
         [
             {
-                "Attributes": <string>,
+                "Attributes": [<int array>],
                 "Direction": <string> or [<float array>],
                 "CoordinateSystem": <string>
             },

docs/src/config/solver.md

@@ -173,7 +173,8 @@ number of eigenmodes of the problem. The available options are:
     "Restart": <int>,
     "AdaptiveTol": <float>,
     "AdaptiveMaxSamples": <int>,
-    "AdaptiveConvergenceMemory": <int>
+    "AdaptiveConvergenceMemory": <int>,
+    "AdaptiveCircuitSynthesis": <bool>
 }
 ```
 
@@ -225,6 +226,16 @@ sampling algorithm for constructing the reduced-order model for adaptive fast fr
 sweep. For example, a memory of "2" requires two consecutive samples which satisfy the
 error tolerance.
 
+`"AdaptiveCircuitSynthesis" [False]` : Usese adaptive reduced order model to print circuit-like
+matrices (inverse inductance ``L^{-1}``, inverse resistance ``R^{-1}``, capacitance ``C`` and basis
+orthogonalization matrix). These matrices are directly normalized to the conventional voltage for
+the external ports. This adds the port fields as a basis function `LumpedPort` to the reduced order
+model. Requires:
+
+  - Adaptive frequency sweep `AdaptiveTol > 0.0` are turned on,
+  - All `LumpedPort` fields are orthogonal to each other,
+  - Only terms with LRC like frequency dependence are currently supported. This means no `WavePort` or `WavePortPEC`, no `Conductivity`, and no second-order `Farfield` boundary conditions.
+
 ### `solver["Driven"]["Samples"]`
 
 ```json

palace/drivers/drivensolver.cpp

@@ -4,6 +4,12 @@
 #include "drivensolver.hpp"
 
 #include <complex>
+#include <iostream>
+#include <Eigen/Core>
+#include <Eigen/Dense>
+#include <Eigen/src/Core/IO.h>
+#include <fmt/core.h>
+#include <fmt/ostream.h>
 #include <mfem.hpp>
 #include "fem/errorindicator.hpp"
 #include "fem/mesh.hpp"
@@ -41,7 +47,8 @@ DrivenSolver::Solve(const std::vector<std::unique_ptr<Mesh>> &mesh) const
   const auto &omega_sample = iodata.solver.driven.sample_f;
 
   bool adaptive = (iodata.solver.driven.adaptive_tol > 0.0);
-  if (adaptive && omega_sample.size() <= iodata.solver.driven.prom_indices.size())
+  if (adaptive && omega_sample.size() <= iodata.solver.driven.prom_indices.size() &&
+      !iodata.solver.driven.adaptive_circuit_synthesis)
   {
     Mpi::Warning("Adaptive frequency sweep requires > {} total frequency samples!\n"
                  "Reverting to uniform sweep!\n",
@@ -213,14 +220,12 @@ ErrorIndicator DrivenSolver::SweepAdaptive(SpaceOperator &space_op) const
   // Configure PROM parameters if not specified.
   double offline_tol = iodata.solver.driven.adaptive_tol;
   int convergence_memory = iodata.solver.driven.adaptive_memory;
-  int max_size_per_excitation = iodata.solver.driven.adaptive_max_size;
+  auto max_size_per_excitation =
+      static_cast<std::size_t>(iodata.solver.driven.adaptive_max_size);
   int nprom_indices = static_cast<int>(iodata.solver.driven.prom_indices.size());
   MFEM_VERIFY(max_size_per_excitation <= 0 || max_size_per_excitation >= nprom_indices,
               "Adaptive frequency sweep must sample at least " << nprom_indices
                                                                << " frequency points!");
-  // Maximum size — no more than nr steps needed.
-  max_size_per_excitation =
-      std::min(max_size_per_excitation, static_cast<int>(omega_sample.size()));
 
   // Allocate negative curl matrix for postprocessing the B-field and vectors for the
   // high-dimensional field solution.
@@ -262,14 +267,20 @@ ErrorIndicator DrivenSolver::SweepAdaptive(SpaceOperator &space_op) const
   RomOperator prom_op(iodata, space_op, max_size_per_excitation);
   space_op.GetWavePortOp().SetSuppressOutput(true);
 
+  // Add ports to PROM if we do synthesis.
+  if (iodata.solver.driven.adaptive_circuit_synthesis)
+  {
+    prom_op.AddLumpedPortModesForSynthesis(iodata);
+  }
+
   // Initialize the basis with samples from the top and bottom of the frequency
   // range of interest. Each call for an HDM solution adds the frequency sample to P_S and
   // removes it from P \ P_S. Timing for the HDM construction and solve is handled inside
   // of the RomOperator.
-  auto UpdatePROM = [&](int excitation_idx, double omega)
+  auto UpdatePROM = [&](int excitation_idx, double omega, std::size_t sample_idx)
   {
     // Add the HDM solution to the PROM reduced basis.
-    prom_op.UpdatePROM(E);
+    prom_op.UpdatePROM(E, fmt::format("sample_e{:d}_s{:d}", excitation_idx, sample_idx));
     prom_op.UpdateMRI(excitation_idx, omega, E);
 
     // Compute B = -1/(iω) ∇ x E on the true dofs, and set the internal GridFunctions in
@@ -315,7 +326,7 @@ ErrorIndicator DrivenSolver::SweepAdaptive(SpaceOperator &space_op) const
       linalg::AXPY(-1.0, E, Eh);
       max_errors.push_back(linalg::Norml2(space_op.GetComm(), Eh) /
                            linalg::Norml2(space_op.GetComm(), E));
-      UpdatePROM(excitation_idx, omega);
+      UpdatePROM(excitation_idx, omega, i);
     }
     // The estimates associated to the end points are assumed inaccurate.
     max_errors[0] = std::numeric_limits<double>::infinity();
@@ -337,18 +348,16 @@ ErrorIndicator DrivenSolver::SweepAdaptive(SpaceOperator &space_op) const
       prom_op.SolveHDM(excitation_idx, omega_star, E);
       prom_op.SolvePROM(excitation_idx, omega_star, Eh);
       linalg::AXPY(-1.0, E, Eh);
+
       max_errors.push_back(linalg::Norml2(space_op.GetComm(), Eh) /
                            linalg::Norml2(space_op.GetComm(), E));
       memory = max_errors.back() < offline_tol ? memory + 1 : 0;
 
       Mpi::Print("\nGreedy iteration {:d} (n = {:d}): ω* = {:.3e} GHz ({:.3e}), error = "
-                 "{:.3e}{}\n",
+                 "{:.3e}, memory = {:d}/{:d}\n",
                  it - it0 + 1, prom_op.GetReducedDimension(), omega_star * unit_GHz,
-                 omega_star, max_errors.back(),
-                 (memory == 0)
-                     ? ""
-                     : fmt::format(", memory = {:d}/{:d}", memory, convergence_memory));
-      UpdatePROM(excitation_idx, omega_star);
+                 omega_star, max_errors.back(), memory, convergence_memory);
+      UpdatePROM(excitation_idx, omega_star, it);
     }
     Mpi::Print("\nAdaptive sampling{} {:d} frequency samples:\n"
                " n = {:d}, error = {:.3e}, tol = {:.3e}, memory = {:d}/{:d}\n",
@@ -363,6 +372,11 @@ ErrorIndicator DrivenSolver::SweepAdaptive(SpaceOperator &space_op) const
   Mpi::Print(" Total offline phase elapsed time: {:.2e} s\n",
              Timer::Duration(Timer::Now() - t0).count());  // Timing on root
 
+  if (iodata.solver.driven.adaptive_circuit_synthesis)
+  {
+    prom_op.PrintPROMMatrices(iodata.units, iodata.problem.output);
+  }
+
   // XX TODO: Add output of eigenvalue estimates from the PROM system (and nonlinear EVP
   // in the general case with wave ports, etc.?)
 

palace/drivers/drivensolver.hpp

@@ -14,8 +14,6 @@ namespace palace
 
 class ErrorIndicator;
 class Mesh;
-template <ProblemType>
-class PostOperator;
 class SpaceOperator;
 
 //

palace/linalg/iterative.cpp

@@ -5,6 +5,7 @@
 
 #include <algorithm>
 #include <cmath>
+#include <cstddef>
 #include <limits>
 #include <string>
 #include "linalg/orthog.hpp"
@@ -306,7 +307,7 @@ inline void ApplyBA(PreconditionerSide side, const OperType *A, const Solver<Ope
 template <typename VecType, typename ScalarType>
 inline void OrthogonalizeIteration(Orthogonalization type, MPI_Comm comm,
                                    const std::vector<VecType> &V, VecType &w,
-                                   ScalarType *Hj, int j)
+                                   ScalarType *Hj, std::size_t j)
 {
   // Orthogonalize w against the leading j + 1 columns of V.
   switch (type)

palace/linalg/nleps.cpp

@@ -283,7 +283,7 @@ void QuasiNewtonSolver::SetInitialGuess()
 
   // If the number of initial guesses is greater than the number of requested modes
   // de-prioritize the initial guesses that have larger errors.
-  std::vector<size_t> indices(nev_linear);
+  std::vector<std::size_t> indices(nev_linear);
   std::iota(indices.begin(), indices.end(), 0);
   if (nev_linear > nev)
   {
@@ -330,8 +330,8 @@ namespace
 ComplexVector MatVecMult(const std::vector<ComplexVector> &X, const Eigen::VectorXcd &y)
 {
   MFEM_ASSERT(X.size() == y.size(), "Mismatch in dimension of input vectors!");
-  const size_t k = X.size();
-  const size_t n = X[0].Size();
+  const std::size_t k = X.size();
+  const std::size_t n = X[0].Size();
   const bool use_dev = X[0].UseDevice();
   ComplexVector z;
   z.SetSize(n);

palace/linalg/orthog.hpp

@@ -4,7 +4,11 @@
 #ifndef PALACE_LINALG_ORTHOG_HPP
 #define PALACE_LINALG_ORTHOG_HPP
 
+#include <complex>
+#include <cstddef>
+#include <memory>
 #include <vector>
+#include "linalg/operator.hpp"
 #include "linalg/vector.hpp"
 #include "utils/communication.hpp"
 
@@ -15,36 +19,131 @@ namespace palace::linalg
 // Orthogonalization functions for orthogonalizing a vector against a number of basis
 // vectors using modified or classical Gram-Schmidt.
 //
+// Assumes that the input vectors are normalized, but does not normalize the output vectors!
+// If done in a loop, normalization has to be managed by hand! (TODO: Reconsider).
+//
+
+// Concept: InnerProductHelper has function InnerProduct(VecType, VecType) -> ScalarType,
+// acting on local degrees of freedom. Also add MPI reduction.
+
+// Simplest case is canonical inner product on R & C.
+
+class InnerProductStandard
+{
+public:
+  double InnerProduct(const Vector &x, const Vector &y) const { return LocalDot(x, y); }
+
+  double InnerProduct(MPI_Comm comm, const Vector &x, const Vector &y) const
+  {
+    return Dot(comm, x, y);
+  }
+
+  std::complex<double> InnerProduct(const ComplexVector &x, const ComplexVector &y) const
+  {
+    return LocalDot(x, y);
+  }
+
+  std::complex<double> InnerProduct(MPI_Comm comm, const ComplexVector &x,
+                                    const ComplexVector &y) const
+  {
+    return Dot(comm, x, y);
+  }
+};
+
+class InnerProductRealWeight
+{
+  // Choose generic operator, although can improve by refining for specialized type.
+  std::shared_ptr<Operator> weight_op;
+
+  // Don't have inner product wrapper / implemented in Operator, so need to allocate a
+  // vector as a workspace. (TODO: Optimize this away).
+  mutable Vector v_workspace = {};
+
+  void SetWorkspace(const Vector &blueprint) const
+  {
+    v_workspace.SetSize(blueprint.Size());
+    v_workspace.UseDevice(blueprint.UseDevice());
+  }
+
+public:
+  template <typename OpType>
+  explicit InnerProductRealWeight(const std::shared_ptr<OpType> &weight_op_)
+    : weight_op(weight_op_)
+  {
+  }
+  // Follow same conventions as Dot:  yᴴ x or yᵀ x (not y comes second in the arguments).
+  double InnerProduct(const Vector &x, const Vector &y) const
+  {
+    SetWorkspace(x);
+    weight_op->Mult(x, v_workspace);
+    return LocalDot(v_workspace, y);
+  }
 
-template <typename VecType, typename ScalarType>
+  double InnerProduct(MPI_Comm comm, const Vector &x, const Vector &y) const
+  {
+    SetWorkspace(x);
+    weight_op->Mult(x, v_workspace);
+    return Dot(comm, v_workspace, y);
+  }
+
+  std::complex<double> InnerProduct(const ComplexVector &x, const ComplexVector &y) const
+  {
+    using namespace std::complex_literals;
+    SetWorkspace(x.Real());
+
+    // weight_op is real.
+    weight_op->Mult(x.Real(), v_workspace);
+    std::complex<double> dot = {+LocalDot(v_workspace, y.Real()),
+                                -LocalDot(v_workspace, y.Imag())};
+
+    weight_op->Mult(x.Imag(), v_workspace);
+    dot += std::complex<double>{LocalDot(v_workspace, y.Imag()),
+                                LocalDot(v_workspace, y.Real())};
+
+    return dot;
+  }
+
+  std::complex<double> InnerProduct(MPI_Comm comm, const ComplexVector &x,
+                                    const ComplexVector &y) const
+  {
+    auto dot = InnerProduct(x, y);
+    Mpi::GlobalSum(1, &dot, comm);
+    return dot;
+  }
+};
+
+template <typename VecType, typename ScalarType,
+          typename InnerProductW = InnerProductStandard>
 inline void OrthogonalizeColumnMGS(MPI_Comm comm, const std::vector<VecType> &V, VecType &w,
-                                   ScalarType *H, int m)
+                                   ScalarType *H, std::size_t m,
+                                   const InnerProductW &dot_op = {})
 {
-  MFEM_ASSERT(static_cast<std::size_t>(m) <= V.size(),
-              "Out of bounds number of columns for MGS orthogonalization!");
-  for (int j = 0; j < m; j++)
+  MFEM_ASSERT(m <= V.size(), "Out of bounds number of columns for MGS orthogonalization!");
+  for (std::size_t j = 0; j < m; j++)
   {
-    H[j] = linalg::Dot(comm, w, V[j]);  // Global inner product
+    // Global inner product: Note order is important for complex vectors.
+    H[j] = dot_op.InnerProduct(comm, w, V[j]);
     w.Add(-H[j], V[j]);
   }
 }
 
-template <typename VecType, typename ScalarType>
+template <typename VecType, typename ScalarType,
+          typename InnerProductW = InnerProductStandard>
 inline void OrthogonalizeColumnCGS(MPI_Comm comm, const std::vector<VecType> &V, VecType &w,
-                                   ScalarType *H, int m, bool refine = false)
+                                   ScalarType *H, std::size_t m, bool refine = false,
+                                   const InnerProductW &dot_op = {})
 {
-  MFEM_ASSERT(static_cast<std::size_t>(m) <= V.size(),
-              "Out of bounds number of columns for CGS orthogonalization!");
+  MFEM_ASSERT(m <= V.size(), "Out of bounds number of columns for CGS orthogonalization!");
   if (m == 0)
   {
     return;
   }
-  for (int j = 0; j < m; j++)
+  for (std::size_t j = 0; j < m; j++)
   {
-    H[j] = w * V[j];  // Local inner product
+    H[j] = dot_op.InnerProduct(w, V[j]);  // Local inner product
   }
   Mpi::GlobalSum(m, H, comm);
-  for (int j = 0; j < m; j++)
+  for (std::size_t j = 0; j < m; j++)
   {
     w.Add(-H[j], V[j]);
   }
@@ -53,10 +152,10 @@ inline void OrthogonalizeColumnCGS(MPI_Comm comm, const std::vector<VecType> &V,
     std::vector<ScalarType> dH(m);
     for (int j = 0; j < m; j++)
     {
-      dH[j] = w * V[j];  // Local inner product
+      dH[j] = dot_op.InnerProduct(w, V[j]);  // Local inner product
     }
     Mpi::GlobalSum(m, dH.data(), comm);
-    for (int j = 0; j < m; j++)
+    for (std::size_t j = 0; j < m; j++)
     {
       H[j] += dH[j];
       w.Add(-dH[j], V[j]);

palace/linalg/rap.cpp

@@ -958,12 +958,17 @@ BuildParSumOperator(ScalarType (&&coeff_in)[N], const OperType *(&&ops_in)[N],
 }
 
 // Explicit instantiation.
+template std::unique_ptr<ParOperator> BuildParSumOperator(double (&&)[1],
+                                                          const Operator *(&&)[1], bool);
 template std::unique_ptr<ParOperator> BuildParSumOperator(double (&&)[2],
                                                           const Operator *(&&)[2], bool);
 template std::unique_ptr<ParOperator> BuildParSumOperator(double (&&)[3],
                                                           const Operator *(&&)[3], bool);
 template std::unique_ptr<ParOperator> BuildParSumOperator(double (&&)[4],
                                                           const Operator *(&&)[4], bool);
+
+template std::unique_ptr<ComplexParOperator>
+    BuildParSumOperator(std::complex<double> (&&)[1], const ComplexOperator *(&&)[1], bool);
 template std::unique_ptr<ComplexParOperator>
     BuildParSumOperator(std::complex<double> (&&)[2], const ComplexOperator *(&&)[2], bool);
 template std::unique_ptr<ComplexParOperator>

palace/models/lumpedportoperator.hpp

@@ -45,7 +45,7 @@ class LumpedPortData
   int excitation;
   bool active;
 
-private:
+protected:
   // Linear forms for postprocessing integrated quantities on the port.
   mutable std::unique_ptr<mfem::LinearForm> s, v;