move to FluidIteration

Author: Bot-Enigma-0

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

@@ -102,9 +102,6 @@ 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 @@ 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,213 +312,3 @@ bool CSinglezoneDriver::Monitor(unsigned long TimeIter){
 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 << 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 << 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);
-      }
-    }
-    // Calculate the global norm
-    passivedouble local_norm = norm;
-    passivedouble global_norm;
-    SU2_MPI::Allreduce(&local_norm, &global_norm, 1, MPI_DOUBLE, MPI_SUM, SU2_MPI::GetComm());
-    global_norm = sqrt(local_norm);
-
-    // normalize 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];
-
-  unsigned long nPoint = geometry->GetnPointDomain();
-  unsigned short nVarTotal = 0;
-  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;
-  }
-
-  su2matrix<passivedouble> w_k(nPoint, nVarTotal);
-  if (rank == MASTER_NODE) {std::cout << "matrix w initialized with size: " << w_k.rows() << " x " << w_k.cols() << endl;}
-
-  passivedouble rho_old = 0.0;
-  int max_iter = config_container[ZONE_0]->GetnPowerMethodIter();
-  su2double tol = config_container[ZONE_0]->GetPowerMethodTolerance();
-  int k = 0;
-
-  while (k < max_iter) {
-    GetAllDerivatives(w_k, geometry); 
-
-    // compute local squared norm
-    passivedouble rho_local = 0.0;
-    for (unsigned long iPoint = 0; iPoint < nPoint; iPoint++) {
-      for (unsigned short iVar = 0; iVar < nVarTotal; iVar++) {
-          rho_local += w_k(iPoint, iVar) * w_k(iPoint, iVar);
-      }
-    }
-
-    // global reduction
-    passivedouble rho_global = 0.0;
-    SU2_MPI::Allreduce(&rho_local, &rho_global, 1, MPI_DOUBLE, MPI_MAX, SU2_MPI::GetComm());
-    passivedouble rho_k = sqrt(rho_local);
-
-    if (rank == MASTER_NODE) {
-      std::cout << "Spectral Radius Estimate for Iteration " << k << " is " << rho_k << endl;
-    }
-
-    // convergence check
-    if (abs(rho_k - rho_old) < tol) {
-      if (rank == MASTER_NODE) {
-        std::cout << "Converged Eigenvalue after " << k << " iterations is " << rho_k << endl;
-      }
-      break; 
-    }
-
-    // normalize w_k to form next seed w_{k+1}
-    for (unsigned long iPoint = 0; iPoint < nPoint; ++iPoint) {
-      for (unsigned short iVar = 0; iVar < nVarTotal; ++iVar) {
-        w_k(iPoint, iVar) /= rho_k;
-      }
-    }
-
-    SeedAllDerivatives(w_k, geometry);
-    rho_old = rho_k;
-    k++;
-
-  } //while
-
-  w_k.resize(0, 0);
-}