[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/iteration/CFluidIteration.hpp

@@ -133,7 +133,19 @@ class CFluidIteration : public CIteration {
                    CVolumetricMovement*** grid_movement, CFreeFormDefBox*** FFDBox, unsigned short val_iZone,
                    unsigned short val_iInst) override;
 
+  /*--- Spectral radius analysis functions ---*/
+  void SeedAllDerivatives(const su2matrix<passivedouble>& derivatives, CGeometry**** geometry, CSolver***** solver, 
+                         unsigned short val_iZone, unsigned short val_iInst);
+  void GetAllDerivatives(su2matrix<passivedouble>& derivatives, CGeometry**** geometry, CSolver***** solver, 
+                        unsigned short val_iZone, unsigned short val_iInst);
+  void PreRunSpectralRadius(CGeometry**** geometry, CSolver***** solver, CConfig** config, 
+                           unsigned short val_iZone, unsigned short val_iInst);
+  void PostRunSpectralRadius(COutput* output, CIntegration**** integration, CGeometry**** geometry, CSolver***** solver,
+                            CNumerics****** numerics, CConfig** config, CSurfaceMovement** surface_movement,
+                            CVolumetricMovement*** grid_movement, CFreeFormDefBox*** FFDBox, unsigned short val_iZone,
+                            unsigned short val_iInst);
  private:
+  su2matrix<passivedouble> v_estimate;
 
   /*!
    * \brief Imposes a gust via the grid velocities.

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

SU2_CFD/src/iteration/CFluidIteration.cpp

@@ -398,6 +398,11 @@
   Preprocess(output, integration, geometry, solver, numerics, config, surface_movement, grid_movement, FFDBox,
              val_iZone, INST_0);
 
+  /*--- Spectral Radius Analysis ---*/
+  if (config[val_iZone]->GetSpectralRadius_Analysis()) {
+    PreRunSpectralRadius(geometry, solver, config, val_iZone, val_iInst);
+  }
+
   /*--- For steady-state flow simulations, we need to loop over ExtIter for the number of time steps ---*/
   /*--- However, ExtIter is the number of FSI iterations, so nIntIter is used in this case ---*/
 
@@ -408,6 +413,17 @@
     Iterate(output, integration, geometry, solver, numerics, config, surface_movement, grid_movement, FFDBox, val_iZone,
             INST_0);
 
+    /*--- Spectral Radius Analysis ---*/
+    // Spectral radius analysis
+    if (config[val_iZone]->GetSpectralRadius_Analysis()) {
+      PostRunSpectralRadius(output, integration, geometry, solver, numerics, config, 
+                            surface_movement, grid_movement, FFDBox, val_iZone, val_iInst);
+
+      // After spectral radius analysis, we break out of the inner loop
+      // since we only want to run one iteration for the analysis
+      break;
+    }
+
     /*--- Monitor the pseudo-time ---*/
     StopCalc = Monitor(output, integration, geometry, solver, numerics, config, surface_movement, grid_movement, FFDBox,
                        val_iZone, INST_0);
@@ -699,3 +715,238 @@
   }
 
 }
+
+void CFluidIteration::SeedAllDerivatives(const su2matrix<passivedouble>& derivatives, CGeometry**** geometry, 
+                                        CSolver***** solver, unsigned short val_iZone, unsigned short val_iInst) {
+  unsigned short targetSolver = 2;
+  unsigned short targetVar = 0;
+
+  // Only process the target solver
+  CSolver* solver_ptr = solver[val_iZone][val_iInst][MESH_0][targetSolver];
+
+  CVariable* nodes = solver_ptr->GetNodes();
+  unsigned long nPoint = geometry[val_iZone][val_iInst][MESH_0]->GetnPointDomain();
+
+  // Seed only the target variable
+  for (unsigned long iPoint = 0; iPoint < nPoint; iPoint++) {
+    su2double* solution = nodes->GetSolution(iPoint);
+    double derivative_value = derivatives(iPoint, 0); // Only one column now
+    SU2_TYPE::SetDerivative(solution[targetVar], derivative_value);
+  }
+
+  if (rank == MASTER_NODE) { 
+    std::cout << "Finished SeedAllDerivatives for solver " << targetSolver 
+              << ", variable " << targetVar << std::endl;
+  }
+}
+
+void CFluidIteration::GetAllDerivatives(su2matrix<passivedouble>& derivatives, CGeometry**** geometry, 
+                                       CSolver***** solver, unsigned short val_iZone, unsigned short val_iInst) {
+  unsigned short target_Solver = 2;
+  unsigned short target_Var = 0;
+
+  //process the required solver
+  CSolver* solver_ptr = solver[val_iZone][val_iInst][MESH_0][target_Solver];
+  if (!solver_ptr) {
+    if (rank == MASTER_NODE) {std::cout << "Target solver " << target_Solver << " is not allocated." << std::endl;}
+    return;
+  }
+
+  CVariable* nodes = solver_ptr->GetNodes();
+  unsigned long nPoint = geometry[val_iZone][val_iInst][MESH_0]->GetnPointDomain();
+
+  //extract derivatives
+  for (unsigned long iPoint = 0; iPoint < nPoint; iPoint++) {
+    su2double* solution = nodes->GetSolution(iPoint);
+    double derivative_value = SU2_TYPE::GetDerivative(solution[target_Var]);
+    derivatives(iPoint, 0) = derivative_value; // Only one column now
+  }
+
+  if (rank == MASTER_NODE) { 
+    std::cout << "Finished GetAllDerivatives for solver " << target_Solver 
+              << ", variable " << target_Var << std::endl;
+  }
+}
+
+void CFluidIteration::PreRunSpectralRadius(CGeometry**** geometry, CSolver***** solver, CConfig** config, 
+                                          unsigned short val_iZone, unsigned short val_iInst) {
+  if (!config[val_iZone]->GetSpectralRadius_Analysis()) return;
+
+  CGeometry* geom_ptr = geometry[val_iZone][val_iInst][MESH_0];
+
+  //get total number of variables and points
+  unsigned long nPoint = geom_ptr->GetnPointDomain();
+
+  //requested solver and variable
+  unsigned short targetSolver = 2;
+  unsigned short targetVar = 0;
+  unsigned long nVarTotal = 1;
+
+    //initialize power iteration matrix
+  if (v_estimate.rows() != nPoint || v_estimate.cols() != nVarTotal) {
+    v_estimate.resize(nPoint, nVarTotal);
+
+    for (unsigned long iPoint = 0; iPoint < nPoint; iPoint++) {
+      v_estimate(iPoint, 0) = static_cast<passivedouble>(1.0);
+    }
+
+    //calculate norm
+    passivedouble local_norm = 0.0;
+    for (unsigned long iPoint = 0; iPoint < nPoint; iPoint++) {
+      local_norm += v_estimate(iPoint, 0) * v_estimate(iPoint, 0);
+    }
+
+    passivedouble global_norm;
+    SU2_MPI::Allreduce(&local_norm, &global_norm, 1, MPI_DOUBLE, MPI_SUM, SU2_MPI::GetComm());
+    global_norm = sqrt(global_norm);
+
+    //normalise
+    for (unsigned long iPoint = 0; iPoint < nPoint; iPoint++) {
+      for (unsigned long iVar = 0; iVar < nVarTotal; iVar++) {
+        v_estimate(iPoint, iVar) /= global_norm;
+      }
+    }
+  }
+
+  //seed all solvers
+  SeedAllDerivatives(v_estimate, geometry, solver, val_iZone, val_iInst);
+
+  if (rank == MASTER_NODE) { 
+    std::cout << "Seeding done" << std::endl;
+  }
+}
+
+void CFluidIteration::PostRunSpectralRadius(COutput* output, CIntegration**** integration, CGeometry**** geometry, CSolver***** solver,
+                            CNumerics****** numerics, CConfig** config, CSurfaceMovement** surface_movement,
+                            CVolumetricMovement*** grid_movement, CFreeFormDefBox*** FFDBox, unsigned short val_iZone,
+                            unsigned short val_iInst) {
+  if (!config[val_iZone]->GetSpectralRadius_Analysis()) return;
+
+  CGeometry* geom_ptr = geometry[val_iZone][val_iInst][MESH_0];
+
+  unsigned long nPoint = geom_ptr->GetnPointDomain();
+  unsigned long nVarTotal = 0;
+
+  for (auto iSol = 0u; iSol < MAX_SOLS; ++iSol) {
+    CSolver* solver_ptr = solver[val_iZone][val_iInst][MESH_0][iSol];
+    if (solver_ptr) nVarTotal += solver_ptr->GetnVar();
+  }
+
+  if (rank == MASTER_NODE) {
+    std::cout << "Total variables: " << nVarTotal << ", Total points: " << nPoint << std::endl;
+  }
+
+  su2matrix<passivedouble> w_k(nPoint, 1);
+
+  if (rank == MASTER_NODE) {
+    std::cout << "Matrix w initialized with size: " << w_k.rows() << " x " << w_k.cols() << std::endl;
+  }
+
+  //vector to store eigenvalue history
+  std::vector<std::pair<int, passivedouble>> eigenvalue_history;
+
+  passivedouble rho_old = 0.0;
+  int max_iter = config[val_iZone]->GetnPowerMethodIter();
+  su2double tol = config[val_iZone]->GetPowerMethodTolerance();
+  int k = 0;
+
+  //for Cauchy convergence criterion (turned off in current code)
+  const int cauchy_window = 10;
+  std::vector<passivedouble> recent_eigenvalues;
+
+  //write matrix to CSV
+  auto write_matrix = [](const su2matrix<passivedouble>& matrix, const std::string& filename) {
+    std::ofstream outfile(filename);
+    outfile << std::scientific << std::setprecision(15);
+    for (unsigned long iPoint = 0; iPoint < matrix.rows(); ++iPoint) {
+      for (unsigned long iVar = 0; iVar < matrix.cols(); ++iVar) {
+        outfile << matrix(iPoint, iVar);
+        if (iVar < matrix.cols() - 1) outfile << ",";
+      }
+      outfile << "\n";
+    }
+    outfile.close();
+    std::cout << "Matrix written to " << filename << std::endl;
+  };
+
+  //write eigenvalue history to CSV - written when converged or when solver ends
+  auto write_eigenvalue_history = [](const std::vector<std::pair<int, passivedouble>>& history, const std::string& filename) {
+    std::ofstream eigenvalue_file(filename);
+    eigenvalue_file << std::scientific << std::setprecision(15);
+    eigenvalue_file << "Iteration,Eigenvalue\n";
+    for (const auto& entry : history) {
+      eigenvalue_file << entry.first << "," << entry.second << "\n";
+    }
+    eigenvalue_file.close();
+    std::cout << "Eigenvalue history written to " << filename << std::endl;
+  };
+
+  while (k < max_iter) {
+    if (rank == MASTER_NODE) {
+      std::cout << "Power iteration " << k << std::endl;
+    }
+    GetAllDerivatives(w_k, geometry, solver, val_iZone, val_iInst);
+    passivedouble rho_local = 0.0;
+    for (unsigned long iPoint = 0; iPoint < nPoint; iPoint++) {
+      rho_local += w_k(iPoint, 0) * w_k(iPoint, 0); // Only column 0
+    }
+
+    passivedouble rho_global = 0.0;
+    SU2_MPI::Allreduce(&rho_local, &rho_global, 1, MPI_DOUBLE, MPI_SUM, SU2_MPI::GetComm());
+    passivedouble rho_k = sqrt(rho_global);
+
+    //store eigenvalue for this iteration with high precision
+    eigenvalue_history.emplace_back(k, rho_k);
+
+    //store recent eigenvalues for Cauchy convergence criterion
+    recent_eigenvalues.push_back(rho_k);
+    if (recent_eigenvalues.size() > cauchy_window) {
+      recent_eigenvalues.erase(recent_eigenvalues.begin());
+    }
+
+    if (rank == MASTER_NODE) {
+      std::cout << std::scientific << std::setprecision(15);
+      std::cout << "Spectral Radius Estimate for Iteration " << k << " is " << rho_k << std::endl;
+      std::cout << std::defaultfloat; // Reset to default formatting
+    }
+
+    //convergence criterion: Cauchy criterion over the last 10 iterations
+    if (std::abs(rho_k - rho_old) < tol) {
+      if (rank == MASTER_NODE) {
+        std::cout << "Converged Eigenvalue after " << k << " iterations is "
+                  << std::scientific << std::setprecision(15) << rho_k << std::endl;
+
+        write_matrix(w_k, "w_k.csv");
+
+        write_eigenvalue_history(eigenvalue_history, "eigenvalue_history.csv");
+      }
+      break;
+    }
+
+    //normalize w_k to form next seed v_estimate
+    for (unsigned long iPoint = 0; iPoint < nPoint; iPoint++) {
+      v_estimate(iPoint, 0) = w_k(iPoint, 0) / rho_k;
+    }
+
+    //seed the derivatives for this iteration
+    SeedAllDerivatives(v_estimate, geometry, solver, val_iZone, val_iInst);
+    rho_old = rho_k;
+    k++;
+
+    //run a single iteration of the solver to apply the Jacobian
+    Iterate(output, integration, geometry, solver, numerics, config, 
+            surface_movement, grid_movement, FFDBox, val_iZone, val_iInst);
+
+    if (k == max_iter) {
+      if (rank == MASTER_NODE) {
+        write_matrix(w_k, "w_k.csv");
+        write_eigenvalue_history(eigenvalue_history, "eigenvalue_history.csv");
+
+        std::cout << "Reached maximum iterations without convergence. Final eigenvalue: " 
+                  << std::scientific << std::setprecision(15) << rho_k << std::endl;
+      }
+    }
+  }
+
+  w_k.resize(0, 0);
+}