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

Common/include/CConfig.hpp

@@ -1078,7 +1078,9 @@ class CConfig {
   WINDOW_FUNCTION Kind_WindowFct;      /*!< \brief Type of window (weight) function for objective functional. */
   unsigned short Kind_HybridRANSLES;   /*!< \brief Kind of Hybrid RANS/LES. */
   unsigned short Kind_RoeLowDiss;      /*!< \brief Kind of Roe scheme with low dissipation for unsteady flows. */
-
+  int PowerIterations;              /*!< \brief Number of iterations for the power method>*/
+  su2double PowerMethodTol;         /*!< \brief Convergence criteria for the dominant eigenvalue using power method>*/
+  bool SpectralRadiusMethod;        /*!< \brief Option to perform spectral radius analysis>*/
   unsigned short nSpanWiseSections; /*!< \brief number of span-wise sections */
   unsigned short nSpanMaxAllZones;  /*!< \brief number of maximum span-wise sections for all zones */
   unsigned short *nSpan_iZones;     /*!< \brief number of span-wise sections for each zones */
@@ -10089,5 +10091,17 @@ class CConfig {
    * \return option data structure for the flamelet fluid model.
    */
   const FluidFlamelet_ParsedOptions& GetFlameletParsedOptions() const { return flamelet_ParsedOptions; }
+
+  /*!
+   * \brief Get requested number of iterations for the Power Method
+   * \return Number of iterations for the Power Method.
+   */
+  int GetnPowerMethodIter(void) const { return PowerIterations; }
+  /*!
+   * \brief Get requested tolerance for the Power Method
+   * \return Value of the power method tolerance.
+   */
+  su2double GetPowerMethodTolerance(void) const { return PowerMethodTol; }
+  bool GetSpectralRadius_Analysis(void) const { return SpectralRadiusMethod; }
 
 };

Common/src/CConfig.cpp

@@ -2829,6 +2829,13 @@ void CConfig::SetConfig_Options() {
   /* DESCRIPTION: Specify the tape which is checked in a tape debug run. */
   addEnumOption("CHECK_TAPE_VARIABLES", AD_CheckTapeVariables, CheckTapeVariables_Map, CHECK_TAPE_VARIABLES::SOLVER_VARIABLES);
 
+  /*!\par CONFIG_CATEGORY: Multi-grid \ingroup Config*/
+  /*!\brief POWER_ITERATION_ITERATION\n DESCRIPTION: Number of power iterations to perform when evaluating the dominant eigenvalue \n DEFAULT: 100 \ingroup Config*/
+  addIntegerOption("POWER_ITERATION_ITERATION", PowerIterations, 100);
+  /*!\brief POWER_ITERATION_TOLERANCE\n DESCRIPTION: Min value of the residual (log10 of the residual)\n DEFAULT: -14.0 \ingroup Config*/
+  addDoubleOption("POWER_ITERATION_TOLERANCE", PowerMethodTol, 1e-7);
+  addBoolOption("SPECTRAL_RADIUS", SpectralRadiusMethod, NO);
+
   /*--- options that are used in the python optimization scripts. These have no effect on the c++ toolsuite ---*/
   /*!\par CONFIG_CATEGORY:Python Options\ingroup Config*/
 

SU2_CFD/include/drivers/CSinglezoneDriver.hpp

@@ -46,6 +46,19 @@ class CSinglezoneDriver : public CDriver {
      * \return Boolean indicating whether the problem is converged.
      */
   virtual bool GetTimeConvergence() const;
+  su2matrix<passivedouble> v_estimate;
+
+  /*!
+  * \brief Seed derivatives for all solvers in the zone.
+  * \param[in] derivatives - Matrix of derivative values
+  */
+  void SeedAllDerivatives(const su2matrix<passivedouble>& derivatives, CGeometry *geometry);
+
+  /*!
+  * \brief Extract derivatives from all solvers in the zone.
+  * \param[out] derivatives - Matrix to store derivative values
+  */
+  void GetAllDerivatives(su2matrix<passivedouble>& derivatives, CGeometry *geometry);
 
 public:
 
@@ -110,4 +123,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,6 @@
 #include "../../../Common/include/graph_coloring_structure.hpp"
 #include "../../../Common/include/toolboxes/MMS/CVerificationSolution.hpp"
 #include "../variables/CVariable.hpp"
-
 #ifdef HAVE_LIBROM
 #include "librom.h"
 #endif
@@ -103,6 +102,9 @@ class CSolver {
   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. */

SU2_CFD/src/drivers/CSinglezoneDriver.cpp

@@ -73,7 +73,7 @@
     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 @@
     /*--- Update the solution for dual time stepping strategy ---*/
 
     Update();
-
+    PostRunSpectralRadius();
     /*--- Monitor the computations after each iteration. ---*/
 
     Monitor(TimeIter);
@@ -312,3 +312,232 @@
 bool CSinglezoneDriver::GetTimeConvergence() const{
   return output_container[ZONE_0]->GetCauchyCorrectedTimeConvergence(config_container[ZONE_0]);
 }
+
+void CSinglezoneDriver::SeedAllDerivatives(const su2matrix<passivedouble>& derivatives, CGeometry *geometry) {
+  unsigned long pointOffset = 0;
+  unsigned long varOffset = 0;
+
+  if (rank == MASTER_NODE) { 
+    std::cout << "runnin SeedAllDerivatives function" << endl;
+    std::cout << "matrix size: " << derivatives.rows() << " x " << derivatives.cols() << endl;
+  }
+
+  for (auto iSol = 0u; iSol < MAX_SOLS; ++iSol) {
+    CSolver* solver = solver_container[ZONE_0][INST_0][MESH_0][iSol];
+    if (!solver) continue;
+
+    CVariable* nodes = solver->GetNodes();
+    if (!nodes) continue;
+
+    unsigned short nVar = solver->GetnVar();
+    unsigned long nPoint = geometry->GetnPointDomain();
+
+    if (rank == MASTER_NODE) { 
+      std::cout << "solver " << iSol << " has nVar=" << nVar << ", nPoint=" << nPoint << endl;
+    }
+
+    // seed derivatives (removed function from csolver to make debug easier)
+    for (unsigned long iPoint = 0; iPoint < nPoint; iPoint++) {
+      su2double* solution = nodes->GetSolution(iPoint);
+      for (unsigned short iVar = 0; iVar < nVar; iVar++) {
+        double derivative_value = derivatives(pointOffset + iPoint, varOffset + iVar);
+
+        // debugging value and derivative 
+        double original_value = SU2_TYPE::GetValue(solution[iVar]);
+        SU2_TYPE::SetDerivative(solution[iVar], derivative_value);
+        double extracted_derivative = SU2_TYPE::GetDerivative(solution[iVar]);
+
+        // see if forward mode works- output for the first few points - 
+        // WORKS WHEN RUNNING SU2_CFD_DIRECTDIFF, when SU2_CFD is run gives "Got= 0"
+        if (iPoint < 5 && iVar == 0 && rank == MASTER_NODE) {
+          std::cout << "Point " << iPoint << ", Var " << iVar << ": Set=" << derivative_value 
+                    << ", Got=" << extracted_derivative << ", Value=" << original_value << std::endl;
+        }
+      }
+    }
+    pointOffset += nPoint;
+    varOffset += nVar;
+  }
+
+  if (rank == MASTER_NODE) { 
+    std::cout << "Finished SeedAllDerivatives" << endl;
+  }
+}
+
+void CSinglezoneDriver::GetAllDerivatives(su2matrix<passivedouble>& derivatives, CGeometry *geometry) {
+  unsigned long pointOffset = 0;
+  unsigned long varOffset = 0;
+
+  if (rank == MASTER_NODE) { 
+    std::cout << "runnin GetAllDerivatives function" << endl;
+  }
+
+  for (auto iSol = 0u; iSol < MAX_SOLS; ++iSol) {
+    CSolver* solver = solver_container[ZONE_0][INST_0][MESH_0][iSol];
+    if (!solver) continue;
+
+    CVariable* nodes = solver->GetNodes();
+    if (!nodes) continue;
+
+    unsigned short nVar = solver->GetnVar();
+    unsigned long nPoint = geometry->GetnPointDomain();
+
+    if (rank == MASTER_NODE) { 
+      std::cout << "solver " << iSol << " has nVar=" << nVar << ", nPoint=" << nPoint << endl;
+    }
+
+    // extract derivatives for solver
+    for (unsigned long iPoint = 0; iPoint < nPoint; iPoint++) {
+      su2double* solution = nodes->GetSolution(iPoint);
+      for (unsigned short iVar = 0; iVar < nVar; iVar++) {
+        double derivative_value = SU2_TYPE::GetDerivative(solution[iVar]);
+        derivatives(pointOffset + iPoint, varOffset + iVar) = derivative_value;
+
+        // Debug: check if derivatives are being extracted correctly
+        if (iPoint == 0 && iVar == 0 && rank == MASTER_NODE) {
+          std::cout << "extracted derivative for solver " << iSol << ", point 0, var 0: " << derivative_value << std::endl;
+        }
+      }
+    }
+    pointOffset += nPoint;
+    varOffset += nVar;
+  }
+
+  if (rank == MASTER_NODE) { 
+    std::cout << "Finished GetAllDerivatives" << endl;
+  }
+}
+
+void CSinglezoneDriver::PreRunSpectralRadius() {
+  if (!config_container[ZONE_0]->GetSpectralRadius_Analysis()) return;
+
+  CGeometry* geometry = geometry_container[ZONE_0][INST_0][MESH_0];
+
+  // get total number of variables and points
+  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();
+  }
+
+  // initialize power iteration matrix if needed
+  if (v_estimate.rows() != nPoint || v_estimate.cols() != nVarTotal) {
+    v_estimate.resize(nPoint, nVarTotal);
+
+    passivedouble norm = 0.0;
+    // initialize with random values
+    for (unsigned long iPoint = 0; iPoint < nPoint; iPoint++) {
+      for (unsigned short iVar = 0; iVar < nVarTotal; iVar++) {
+        v_estimate(iPoint, iVar) = static_cast<passivedouble>(rand()) / RAND_MAX;
+        norm += v_estimate(iPoint, iVar) * v_estimate(iPoint, iVar);
+      }
+    }
+
+    // norm calculation
+    passivedouble global_norm;
+    SU2_MPI::Allreduce(&norm, &global_norm, 1, MPI_DOUBLE, MPI_SUM, SU2_MPI::GetComm());
+    global_norm = sqrt(global_norm);
+
+    // normalize the matrix
+    for (unsigned long iPoint = 0; iPoint < nPoint; iPoint++) {
+      for (unsigned short iVar = 0; iVar < nVarTotal; iVar++) {
+        v_estimate(iPoint, iVar) /= global_norm;
+      }
+    }
+  }
+
+  // seed all solvers
+  SeedAllDerivatives(v_estimate, geometry);
+  if (rank == MASTER_NODE) { std::cout << "seeding done" << endl; }
+}
+
+void CSinglezoneDriver::PostRunSpectralRadius() {
+  if (!config_container[ZONE_0]->GetSpectralRadius_Analysis()) return;
+
+  if (rank == MASTER_NODE) { 
+    std::cout << "Running PostRunSpectralRadius" << endl;
+  }
+
+  CGeometry* geometry = geometry_container[ZONE_0][INST_0][MESH_0];
+
+  // get total number of variables and points
+  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();
+  }
+
+  if (rank == MASTER_NODE) { 
+    std::cout << "total var: " << nVarTotal << ", total points: " << nPoint << endl;
+  }
+
+  // extract derivatives from all solvers
+  su2matrix<passivedouble> w(nPoint, nVarTotal);
+
+  if (rank == MASTER_NODE) { 
+    std::cout << "matrix w initialized with size: " << w.rows() << " x " << w.cols() << endl;
+  }
+
+  GetAllDerivatives(w, geometry);
+
+  if (rank == MASTER_NODE) { 
+    std::cout << "derivatives extracted" << endl;
+  }
+
+  // compute norm 
+  passivedouble rho_new = 0.0;
+  for (unsigned long iPoint = 0; iPoint < nPoint; iPoint++) {
+    for (unsigned short iVar = 0; iVar < nVarTotal; iVar++) {
+      rho_new += w(iPoint, iVar) * w(iPoint, iVar);
+    }
+  }
+
+  passivedouble global_rho_new;
+  SU2_MPI::Allreduce(&rho_new, &global_rho_new, 1, MPI_DOUBLE, MPI_SUM, SU2_MPI::GetComm());
+  global_rho_new = sqrt(global_rho_new);
+
+  if (rank == MASTER_NODE) { 
+    std::cout << "Global rho_new: " << global_rho_new << endl;
+  }
+
+  // update eigenvector estimate
+  for (unsigned long iPoint = 0; iPoint < nPoint; iPoint++) {
+    for (unsigned short iVar = 0; iVar < nVarTotal; iVar++) {
+      v_estimate(iPoint, iVar) = w(iPoint, iVar) / global_rho_new;
+    }
+  }
+
+  // check convergence
+  static passivedouble rho_old = 0.0;
+  static int iter = 0;
+
+  int max_iter = config_container[ZONE_0]->GetnPowerMethodIter();
+  su2double tol = config_container[ZONE_0]->GetPowerMethodTolerance();
+
+  if (rank == MASTER_NODE) { 
+    std::cout << "Power Iterations requested = " << max_iter << endl;
+    std::cout << "Power tolerance requested = " << tol << endl;
+    std::cout << "Current iteration = " << iter << endl;
+    std::cout << "Current rho_old = " << rho_old << endl;
+    std::cout << "Current global_rho_new = " << global_rho_new << endl;
+  }
+
+  if (abs(global_rho_new - rho_old) < tol || iter >= max_iter) {
+    if (rank == MASTER_NODE) {
+      std::cout << "Spectral Radius Estimate: " << global_rho_new << " after " << iter << " iterations." << std::endl;
+    }
+    // Reset
+    rho_old = 0.0;
+    iter = 0;
+    v_estimate.resize(0, 0); 
+  } else {
+    rho_old = global_rho_new;
+    iter++;
+
+    if (rank == MASTER_NODE) { 
+      std::cout << "Continuing to iteration " << iter << " with rho_old = " << rho_old << endl;
+    }
+  }
+}