[WIP] Power Iteration Method for Numerical Stability Analysis of the Non-Linear System

SU2_CFD/include/drivers/CSinglezoneDriver.hpp

@@ -46,6 +46,52 @@ class CSinglezoneDriver : public CDriver {
      * \return Boolean indicating whether the problem is converged.
      */
   virtual bool GetTimeConvergence() const;
+  vector<passivedouble> v_estimate; 
+
+  /*!
+   * \brief Seed derivatives for all solvers in the zone.
+   * \param[in] derivatives - Vector of derivative values
+   */
+  void SeedAllDerivatives(const vector<passivedouble>& derivatives, CVariable* nodes, CGeometry *geometry) {
+    unsigned long offset = 0;
+    for (auto iSol = 0u; iSol < MAX_SOLS; ++iSol) {
+      CSolver* solver = solver_container[ZONE_0][INST_0][MESH_0][iSol];
+      if (!solver) continue;
+
+      unsigned short nVar = solver->GetnVar();
+      unsigned long nPoint = geometry->GetnPointDomain();
+      unsigned long nTotal = nVar * nPoint;
+
+      vector<passivedouble> solver_derivs(derivatives.begin() + offset, 
+                                              derivatives.begin() + offset + nTotal);
+      solver->SetSolutionDerivatives(solver_derivs, nodes);
+      offset += nTotal;
+    }
+  }
+
+  /*!
+   * \brief Extract derivatives from all solvers in the zone.
+   * \param[out] derivatives - Vector to store derivative values
+   */
+  void GetAllDerivatives(vector<passivedouble>& derivatives, CVariable* nodes, CGeometry *geometry) {
+      unsigned long offset = 0;
+      for (auto iSol = 0u; iSol < MAX_SOLS; ++iSol) {
+        CSolver* solver = solver_container[ZONE_0][INST_0][MESH_0][iSol];
+        if (!solver) continue;
+
+        unsigned short nVar = solver->GetnVar();
+        unsigned long nPoint = geometry->GetnPointDomain();
+        unsigned long nTotal = nVar * nPoint;
+
+        vector<passivedouble> solver_derivs(nTotal);
+        solver->GetSolutionDerivatives(solver_derivs, nodes);
+
+        for (unsigned long i = 0; i < nTotal; i++) {
+          derivatives[offset + i] = solver_derivs[i];
+        }
+        offset += nTotal;
+      }
+  }
 
 public:
 
@@ -110,4 +156,10 @@ class CSinglezoneDriver : public CDriver {
    */
   bool Monitor(unsigned long TimeIter) override;
 
+  /*!
+   * \brief Calculate the spectral radius using power iteration method.
+   */
+  void PreRunSpectralRadius();
+
+  void PostRunSpectralRadius();
 };

SU2_CFD/include/solvers/CSolver.hpp

@@ -56,7 +56,7 @@
 #include "../../../Common/include/graph_coloring_structure.hpp"
 #include "../../../Common/include/toolboxes/MMS/CVerificationSolution.hpp"
 #include "../variables/CVariable.hpp"
-
+// #include "../drivers/CDriver.hpp"
 #ifdef HAVE_LIBROM
 #include "librom.h"
 #endif
@@ -103,6 +103,9 @@
   su2double Total_Custom_ObjFunc = 0.0; /*!< \brief Total custom objective function. */
   su2double Total_ComboObj = 0.0;       /*!< \brief Total 'combo' objective for all monitored boundaries */
 
+  CSysVector<su2double> seed_vector; // Current eigenvector estimate
+  su2double spectral_radius;
+
   /*--- Variables that need to go. ---*/
 
   su2double *Residual,      /*!< \brief Auxiliary nVar vector. */
@@ -4352,6 +4355,34 @@
   AD::EndNoSharedReading();
 }
 
+  /*!
+   * \brief Seed derivatives for all solution variables using a flat vector.
+   * \param[in] derivatives - Flat vector of derivative values (size nVar * nPoint)
+   */
+  void SetSolutionDerivatives(const vector<passivedouble>& derivatives, CVariable* nodes) {
+    unsigned long offset = 0;
+    for (unsigned long iPoint = 0; iPoint < nPoint; iPoint++) {
+      su2double* solution = nodes->GetSolution(iPoint);
+      for (unsigned short iVar = 0; iVar < nVar; iVar++) {
+        SU2_TYPE::SetDerivative(solution[iVar], derivatives[offset++]);
+      }
+    }
+  }
+
+  /*!
+   * \brief Get derivatives from all solution variables into a flat vector.
+   * \param[out] derivatives - Flat vector to store derivative values (size nVar * nPoint)
+   */
+  void GetSolutionDerivatives(vector<passivedouble>& derivatives, CVariable* nodes) const {
+    unsigned long offset = 0;
+    for (unsigned long iPoint = 0; iPoint < nPoint; iPoint++) {
+      su2double* solution = nodes->GetSolution(iPoint);
+      for (unsigned short iVar = 0; iVar < nVar; iVar++) {
+        derivatives[offset++] = SU2_TYPE::GetDerivative(solution[iVar]);
+      }
+    }
+  }
+
 protected:
   /*!
    * \brief Allocate the memory for the verification solution, if necessary.

SU2_CFD/src/drivers/CSinglezoneDriver.cpp

@@ -73,7 +73,7 @@ void CSinglezoneDriver::StartSolver() {
     Preprocess(TimeIter);
 
     /*--- Run a time-step iteration of the single-zone problem. ---*/
-
+    PreRunSpectralRadius();
     Run();
 
     /*--- Perform some postprocessing on the solution before the update ---*/
@@ -83,7 +83,7 @@ void CSinglezoneDriver::StartSolver() {
     /*--- Update the solution for dual time stepping strategy ---*/
 
     Update();
-
+    PostRunSpectralRadius();
     /*--- Monitor the computations after each iteration. ---*/
 
     Monitor(TimeIter);
@@ -312,3 +312,78 @@ bool CSinglezoneDriver::Monitor(unsigned long TimeIter){
 bool CSinglezoneDriver::GetTimeConvergence() const{
   return output_container[ZONE_0]->GetCauchyCorrectedTimeConvergence(config_container[ZONE_0]);
 }
+
+void CSinglezoneDriver::PreRunSpectralRadius() {
+  CGeometry* geometry = nullptr;
+  CVariable* nodes = nullptr;
+
+  // Get total number of variables
+  unsigned long nVarTotal = 0;
+  unsigned long nPoint = geometry->GetnPointDomain();
+  for (auto iSol = 0u; iSol < MAX_SOLS; ++iSol) {
+      CSolver* solver = solver_container[ZONE_0][INST_0][MESH_0][iSol];
+      if (solver) nVarTotal += solver->GetnVar();
+  }
+
+  unsigned long nTotal = nVarTotal * nPoint;
+
+  // Initialize power iteration vector 
+  if (v_estimate.empty()) {
+      v_estimate.resize(nTotal);
+      passivedouble norm = 0.0;
+      for (unsigned long i = 0; i < nTotal; i++) {
+          v_estimate[i] = static_cast<passivedouble>(rand()) / RAND_MAX;
+          norm += v_estimate[i] * v_estimate[i];
+      }
+      norm = sqrt(norm);
+      for (unsigned long i = 0; i < nTotal; i++) v_estimate[i] /= norm;
+  }
+
+  // Seed derivatives for all solvers
+  SeedAllDerivatives(v_estimate, nodes, geometry);
+}
+
+void CSinglezoneDriver::PostRunSpectralRadius() {
+    CGeometry* geometry = nullptr;
+    CVariable* nodes = nullptr;
+
+    // Get total number of variables (var*nodes)
+    unsigned long nVarTotal = 0;
+    unsigned long nPoint = geometry->GetnPointDomain();
+    for (auto iSol = 0u; iSol < MAX_SOLS; ++iSol) {
+        CSolver* solver = solver_container[ZONE_0][INST_0][MESH_0][iSol];
+        if (solver) nVarTotal += solver->GetnVar();
+    }
+    unsigned long nTotal = nVarTotal * nPoint;
+
+    // Extract derivatives from all solvers
+    vector<passivedouble> w(nTotal);
+    GetAllDerivatives(w, nodes, geometry);
+
+    // Compute norm and update eigenvector estimate
+    passivedouble rho_new = 0.0;
+    for (unsigned long i = 0; i < nTotal; i++) rho_new += w[i] * w[i];
+    rho_new = sqrt(rho_new);
+
+    for (unsigned long i = 0; i < nTotal; i++) v_estimate[i] = w[i] / rho_new;
+
+    // Check convergence
+    static passivedouble rho_old = 0.0;
+    static int iter = 0;
+
+    int max_iter = 500; //add config option later
+    passivedouble tol = 1e-7; //add config option later
+
+    if (abs(rho_new - rho_old) < tol || iter >= max_iter) {
+        if (rank == MASTER_NODE) {
+            std::cout << "Spectral Radius Estimate: " << rho_new << " after " << iter << " iterations." << std::endl;
+        }
+
+        rho_old = 0.0;
+        iter = 0;
+        v_estimate.clear();
+    } else {
+        rho_old = rho_new;
+        iter++;
+    }
+}