Line data    Source code

       1              : // Copyright Amazon.com, Inc. or its affiliates. All Rights Reserved.
       2              : // SPDX-License-Identifier: Apache-2.0
       3              : 
       4              : #ifndef PALACE_MODELS_ROM_OPERATOR_HPP
       5              : #define PALACE_MODELS_ROM_OPERATOR_HPP
       6              : 
       7              : #include <complex>
       8              : #include <memory>
       9              : #include <vector>
      10              : #include <Eigen/Dense>
      11              : #include "linalg/ksp.hpp"
      12              : #include "linalg/operator.hpp"
      13              : #include "linalg/vector.hpp"
      14              : 
      15              : namespace palace
      16              : {
      17              : 
      18              : class IoData;
      19              : class SpaceOperator;
      20              : 
      21              : // Class for handling minimal-rational interpolation of solutions in frequency space. Used
      22              : // as an error indicator and for selecting the next frequency sample points in PROM
      23              : // construction. Each excitation gets a separate MRI, so sample frequencies are not shared.
      24              : class MinimalRationalInterpolation
      25              : {
      26              : private:
      27              :   // (Complex-valued) upper-trianglar matrix R from orthogonalization of the HDM samples.
      28              :   // Minimal rational interpolant (MRI) defined by the vector q of interpolation weights and
      29              :   // support points z is used as an error indicator.
      30              :   std::vector<ComplexVector> Q;
      31              :   std::size_t dim_Q = 0;
      32              :   Eigen::MatrixXcd R;
      33              :   Eigen::VectorXcd q;
      34              :   std::vector<double> z;
      35              : 
      36              : public:
      37              :   MinimalRationalInterpolation(int max_size);
      38              :   void AddSolutionSample(double omega, const ComplexVector &u,
      39              :                          const SpaceOperator &space_op, Orthogonalization orthog_type);
      40              :   std::vector<double> FindMaxError(int N) const;
      41              : 
      42            0 :   const auto &GetSamplePoints() const { return z; }
      43              : };
      44              : 
      45              : //
      46              : // A class handling projection-based reduced order model (PROM) construction and use for
      47              : // adaptive fast frequency sweeps.
      48              : //
      49              : class RomOperator
      50              : {
      51              : private:
      52              :   // Reference to HDM discretization (not owned).
      53              :   SpaceOperator &space_op;
      54              : 
      55              :   // Used for constructing & resuse of RHS1.
      56              :   int excitation_idx_cache = 0;
      57              : 
      58              :   // HDM system matrices and excitation RHS.
      59              :   std::unique_ptr<ComplexOperator> K, M, C, A2;
      60              :   ComplexVector RHS1, RHS2, r;
      61              :   // Defaults: will be toggeled by SetExcitationIndex & SolveHDM.
      62              :   bool has_A2 = true;
      63              :   bool has_RHS1 = true;
      64              :   bool has_RHS2 = true;
      65              : 
      66              :   // HDM linear system solver and preconditioner.
      67              :   std::unique_ptr<ComplexKspSolver> ksp;
      68              : 
      69              :   // PROM matrices and vectors.
      70              :   Eigen::MatrixXcd Kr, Mr, Cr, Ar;
      71              :   Eigen::VectorXcd RHS1r;
      72              :   Eigen::VectorXcd RHSr;
      73              : 
      74              :   // PROM reduced-order basis (real-valued) and active dimension.
      75              :   std::vector<Vector> V;
      76              :   std::size_t dim_V = 0;
      77              :   Orthogonalization orthog_type;
      78              : 
      79              :   // MRIs: one for each excitation index.
      80              :   std::map<int, MinimalRationalInterpolation> mri;
      81              : 
      82              : public:
      83              :   RomOperator(const IoData &iodata, SpaceOperator &space_op, int max_size_per_excitation);
      84              : 
      85              :   // Return the HDM linear solver.
      86              :   const ComplexKspSolver &GetLinearSolver() const { return *ksp; }
      87              : 
      88              :   // Return PROM dimension.
      89            0 :   auto GetReducedDimension() const { return dim_V; }
      90              : 
      91              :   // Return set of sampled parameter points for basis construction.
      92              :   const auto &GetSamplePoints(int excitation_idx) const
      93              :   {
      94            0 :     return mri.at(excitation_idx).GetSamplePoints();
      95              :   }
      96              : 
      97              :   // Set excitation index to build corresponding RHS vector (linear in frequency part).
      98              :   void SetExcitationIndex(int excitation_idx);
      99              : 
     100              :   // Assemble and solve the HDM at the specified frequency.
     101              :   void SolveHDM(int excitation_idx, double omega, ComplexVector &u);
     102              : 
     103              :   // Add field configuration to the reduced-order basis and update the PROM.
     104              :   void UpdatePROM(const ComplexVector &u);
     105              : 
     106              :   // Add solution u to the minimal-rational interpolation for error estimation. MRI are
     107              :   // separated by excitation index.
     108              :   void UpdateMRI(int excitation_idx, double omega, const ComplexVector &u);
     109              : 
     110              :   // Assemble and solve the PROM at the specified frequency, expanding the solution back
     111              :   // into the high-dimensional space.
     112              :   void SolvePROM(int excitation_idx, double omega, ComplexVector &u);
     113              : 
     114              :   // Compute the location(s) of the maximum error in the range of the previously sampled
     115              :   // parameter points.
     116              :   std::vector<double> FindMaxError(int excitation_idx, int N = 1) const
     117              :   {
     118            0 :     return mri.at(excitation_idx).FindMaxError(N);
     119              :   }
     120              : 
     121              :   // Compute eigenvalue estimates for the current PROM system.
     122              :   std::vector<std::complex<double>> ComputeEigenvalueEstimates() const;
     123              : };
     124              : 
     125              : }  // namespace palace
     126              : 
     127              : #endif  // PALACE_MODELS_ROM_OPERATOR_HPP