From 4256c9840017a8c2ac1aa2aa7532b6188ce58302 Mon Sep 17 00:00:00 2001 From: Bot-Enigma-0 Date: Sat, 13 Sep 2025 21:11:33 +0100 Subject: [PATCH 1/8] power iteration method --- SU2_CFD/include/drivers/CSinglezoneDriver.hpp | 52 ++++++++++++ SU2_CFD/include/solvers/CSolver.hpp | 33 +++++++- SU2_CFD/src/drivers/CSinglezoneDriver.cpp | 79 ++++++++++++++++++- 3 files changed, 161 insertions(+), 3 deletions(-) diff --git a/SU2_CFD/include/drivers/CSinglezoneDriver.hpp b/SU2_CFD/include/drivers/CSinglezoneDriver.hpp index 355369556ea5..f3ed343242de 100644 --- a/SU2_CFD/include/drivers/CSinglezoneDriver.hpp +++ b/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 v_estimate; + + /*! + * \brief Seed derivatives for all solvers in the zone. + * \param[in] derivatives - Vector of derivative values + */ + void SeedAllDerivatives(const vector& 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 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& 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 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(); }; diff --git a/SU2_CFD/include/solvers/CSolver.hpp b/SU2_CFD/include/solvers/CSolver.hpp index 7c38e054d4a1..0f7628b58352 100644 --- a/SU2_CFD/include/solvers/CSolver.hpp +++ b/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 @@ 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 seed_vector; // Current eigenvector estimate + su2double spectral_radius; + /*--- Variables that need to go. ---*/ su2double *Residual, /*!< \brief Auxiliary nVar vector. */ @@ -4352,6 +4355,34 @@ inline void CustomSourceResidual(CGeometry *geometry, CSolver **solver_container 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& 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& 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. diff --git a/SU2_CFD/src/drivers/CSinglezoneDriver.cpp b/SU2_CFD/src/drivers/CSinglezoneDriver.cpp index 1270eef237ab..20890bc63376 100644 --- a/SU2_CFD/src/drivers/CSinglezoneDriver.cpp +++ b/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(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 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++; + } +} \ No newline at end of file From 4d460791072351ba94095730e4865b938cbb410f Mon Sep 17 00:00:00 2001 From: Bot-Enigma-0 Date: Sun, 14 Sep 2025 16:32:01 +0100 Subject: [PATCH 2/8] fix geometry pointer --- SU2_CFD/include/solvers/CSolver.hpp | 1 - SU2_CFD/src/drivers/CSinglezoneDriver.cpp | 95 ++++++++++++----------- 2 files changed, 50 insertions(+), 46 deletions(-) diff --git a/SU2_CFD/include/solvers/CSolver.hpp b/SU2_CFD/include/solvers/CSolver.hpp index 0f7628b58352..dc37d6aad2d6 100644 --- a/SU2_CFD/include/solvers/CSolver.hpp +++ b/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" -// #include "../drivers/CDriver.hpp" #ifdef HAVE_LIBROM #include "librom.h" #endif diff --git a/SU2_CFD/src/drivers/CSinglezoneDriver.cpp b/SU2_CFD/src/drivers/CSinglezoneDriver.cpp index 20890bc63376..1ec8fdd4b3ba 100644 --- a/SU2_CFD/src/drivers/CSinglezoneDriver.cpp +++ b/SU2_CFD/src/drivers/CSinglezoneDriver.cpp @@ -314,7 +314,7 @@ bool CSinglezoneDriver::GetTimeConvergence() const{ } void CSinglezoneDriver::PreRunSpectralRadius() { - CGeometry* geometry = nullptr; + CGeometry* geometry = geometry_container[ZONE_0][INST_0][MESH_0]; CVariable* nodes = nullptr; // Get total number of variables @@ -327,7 +327,7 @@ void CSinglezoneDriver::PreRunSpectralRadius() { unsigned long nTotal = nVarTotal * nPoint; - // Initialize power iteration vector + // Initialize power iteration vector if (v_estimate.empty()) { v_estimate.resize(nTotal); passivedouble norm = 0.0; @@ -341,49 +341,54 @@ void CSinglezoneDriver::PreRunSpectralRadius() { // Seed derivatives for all solvers SeedAllDerivatives(v_estimate, nodes, geometry); + if (rank == MASTER_NODE) {std::cout << "seeding done" << endl;}; //see where code fails } 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 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++; - } -} \ No newline at end of file + if (rank == MASTER_NODE) {std::cout << "running post spectral function" << endl;}; //see where code fails + + CGeometry* geometry = geometry_container[ZONE_0][INST_0][MESH_0]; + 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; + if (rank == MASTER_NODE) {std::cout <<"extracting tangents" << endl;}; //see where code fails + + // Extract derivatives from all solvers + vector w(nTotal); + GetAllDerivatives(w, nodes, geometry); + if (rank == MASTER_NODE) {std::cout << "computing eigen" << endl;}; //see where code fails + + // Compute norm and update eigen 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; + passivedouble tol = 1e-7; + + 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; + } + // Reset for next analysis if needed + rho_old = 0.0; + iter = 0; + v_estimate.clear(); + } else { + rho_old = rho_new; + iter++; + } +} From 39a2375764225157e2f0831472dd003bab4431a2 Mon Sep 17 00:00:00 2001 From: Bot-Enigma-0 Date: Mon, 15 Sep 2025 21:04:53 +0100 Subject: [PATCH 3/8] add options --- Common/include/CConfig.hpp | 16 +++++++++++++++- Common/src/CConfig.cpp | 7 +++++++ 2 files changed, 22 insertions(+), 1 deletion(-) diff --git a/Common/include/CConfig.hpp b/Common/include/CConfig.hpp index 4439410efff4..51396c8c4af8 100644 --- a/Common/include/CConfig.hpp +++ b/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; } }; diff --git a/Common/src/CConfig.cpp b/Common/src/CConfig.cpp index 3cadb5be5ef1..c006a55abd6a 100644 --- a/Common/src/CConfig.cpp +++ b/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*/ From b40e595f69bb3093a5003d852da0f78d945e8bcc Mon Sep 17 00:00:00 2001 From: Bot-Enigma-0 Date: Mon, 15 Sep 2025 21:05:21 +0100 Subject: [PATCH 4/8] apply suggestions --- SU2_CFD/include/drivers/CSinglezoneDriver.hpp | 51 +--- SU2_CFD/include/solvers/CSolver.hpp | 28 -- SU2_CFD/src/drivers/CSinglezoneDriver.cpp | 243 ++++++++++++++---- 3 files changed, 205 insertions(+), 117 deletions(-) diff --git a/SU2_CFD/include/drivers/CSinglezoneDriver.hpp b/SU2_CFD/include/drivers/CSinglezoneDriver.hpp index f3ed343242de..c697ff67915c 100644 --- a/SU2_CFD/include/drivers/CSinglezoneDriver.hpp +++ b/SU2_CFD/include/drivers/CSinglezoneDriver.hpp @@ -46,52 +46,19 @@ class CSinglezoneDriver : public CDriver { * \return Boolean indicating whether the problem is converged. */ virtual bool GetTimeConvergence() const; - vector v_estimate; + su2matrix v_estimate; /*! - * \brief Seed derivatives for all solvers in the zone. - * \param[in] derivatives - Vector of derivative values - */ - void SeedAllDerivatives(const vector& 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 solver_derivs(derivatives.begin() + offset, - derivatives.begin() + offset + nTotal); - solver->SetSolutionDerivatives(solver_derivs, nodes); - offset += nTotal; - } - } + * \brief Seed derivatives for all solvers in the zone. + * \param[in] derivatives - Matrix of derivative values + */ + void SeedAllDerivatives(const su2matrix& derivatives, CGeometry *geometry); /*! - * \brief Extract derivatives from all solvers in the zone. - * \param[out] derivatives - Vector to store derivative values - */ - void GetAllDerivatives(vector& 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 solver_derivs(nTotal); - solver->GetSolutionDerivatives(solver_derivs, nodes); - - for (unsigned long i = 0; i < nTotal; i++) { - derivatives[offset + i] = solver_derivs[i]; - } - offset += nTotal; - } - } + * \brief Extract derivatives from all solvers in the zone. + * \param[out] derivatives - Matrix to store derivative values + */ + void GetAllDerivatives(su2matrix& derivatives, CGeometry *geometry); public: diff --git a/SU2_CFD/include/solvers/CSolver.hpp b/SU2_CFD/include/solvers/CSolver.hpp index dc37d6aad2d6..811846d11479 100644 --- a/SU2_CFD/include/solvers/CSolver.hpp +++ b/SU2_CFD/include/solvers/CSolver.hpp @@ -4354,34 +4354,6 @@ inline void CustomSourceResidual(CGeometry *geometry, CSolver **solver_container 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& 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& 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. diff --git a/SU2_CFD/src/drivers/CSinglezoneDriver.cpp b/SU2_CFD/src/drivers/CSinglezoneDriver.cpp index 1ec8fdd4b3ba..dca05778a735 100644 --- a/SU2_CFD/src/drivers/CSinglezoneDriver.cpp +++ b/SU2_CFD/src/drivers/CSinglezoneDriver.cpp @@ -313,82 +313,231 @@ bool CSinglezoneDriver::GetTimeConvergence() const{ return output_container[ZONE_0]->GetCauchyCorrectedTimeConvergence(config_container[ZONE_0]); } +void CSinglezoneDriver::SeedAllDerivatives(const su2matrix& 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& 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]; - CVariable* nodes = nullptr; - // Get total number of variables + // 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(); + 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(rand()) / RAND_MAX; - norm += v_estimate[i] * v_estimate[i]; + // 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(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; } - 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); - if (rank == MASTER_NODE) {std::cout << "seeding done" << endl;}; //see where code fails + // seed all solvers + SeedAllDerivatives(v_estimate, geometry); + if (rank == MASTER_NODE) { std::cout << "seeding done" << endl; } } void CSinglezoneDriver::PostRunSpectralRadius() { - if (rank == MASTER_NODE) {std::cout << "running post spectral function" << endl;}; //see where code fails + 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]; - CVariable* nodes = nullptr; - // Get total number of variables + // 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(); + 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; } - unsigned long nTotal = nVarTotal * nPoint; - if (rank == MASTER_NODE) {std::cout <<"extracting tangents" << endl;}; //see where code fails - // Extract derivatives from all solvers - vector w(nTotal); - GetAllDerivatives(w, nodes, geometry); - if (rank == MASTER_NODE) {std::cout << "computing eigen" << endl;}; //see where code fails + // extract derivatives from all solvers + su2matrix 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 and update eigen estimate + // compute norm 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 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; + } - for (unsigned long i = 0; i < nTotal; i++) v_estimate[i] = w[i] / rho_new; + // 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 + // check convergence static passivedouble rho_old = 0.0; static int iter = 0; - int max_iter = 500; - passivedouble tol = 1e-7; + int max_iter = config_container[ZONE_0]->GetnPowerMethodIter(); + su2double tol = config_container[ZONE_0]->GetPowerMethodTolerance(); - 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; - } - // Reset for next analysis if needed - rho_old = 0.0; - iter = 0; - v_estimate.clear(); + 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 = rho_new; - iter++; + rho_old = global_rho_new; + iter++; + + if (rank == MASTER_NODE) { + std::cout << "Continuing to iteration " << iter << " with rho_old = " << rho_old << endl; + } } -} +} \ No newline at end of file From 812331b7ec695aa7d0d6071cef0bee7b517b359f Mon Sep 17 00:00:00 2001 From: Bot-Enigma-0 Date: Tue, 16 Sep 2025 20:27:04 +0100 Subject: [PATCH 5/8] correct power method loop --- SU2_CFD/src/drivers/CSinglezoneDriver.cpp | 131 +++++++++------------- 1 file changed, 56 insertions(+), 75 deletions(-) diff --git a/SU2_CFD/src/drivers/CSinglezoneDriver.cpp b/SU2_CFD/src/drivers/CSinglezoneDriver.cpp index dca05778a735..8ff35fe38c1e 100644 --- a/SU2_CFD/src/drivers/CSinglezoneDriver.cpp +++ b/SU2_CFD/src/drivers/CSinglezoneDriver.cpp @@ -351,7 +351,7 @@ void CSinglezoneDriver::SeedAllDerivatives(const su2matrix& deriv // 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; + << ", Got=" << extracted_derivative << ", Value=" << original_value << endl; } } } @@ -395,7 +395,7 @@ void CSinglezoneDriver::GetAllDerivatives(su2matrix& derivatives, // 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; + std::cout << "extracted derivative for solver " << iSol << ", point 0, var 0: " << derivative_value << endl; } } } @@ -433,13 +433,13 @@ void CSinglezoneDriver::PreRunSpectralRadius() { norm += v_estimate(iPoint, iVar) * v_estimate(iPoint, iVar); } } - - // norm calculation + // Calculate the global norm + passivedouble local_norm = norm; 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 + 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; @@ -454,90 +454,71 @@ void CSinglezoneDriver::PreRunSpectralRadius() { void CSinglezoneDriver::PostRunSpectralRadius() { if (!config_container[ZONE_0]->GetSpectralRadius_Analysis()) return; - - if (rank == MASTER_NODE) { + 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(); + 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) { + + if (rank == MASTER_NODE) { std::cout << "total var: " << nVarTotal << ", total points: " << nPoint << endl; } - // extract derivatives from all solvers - su2matrix 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; - } + su2matrix w_k(nPoint, nVarTotal); + if (rank == MASTER_NODE) {std::cout << "matrix w initialized with size: " << w_k.rows() << " x " << w_k.cols() << 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; - } + passivedouble rho_old = 0.0; + int max_iter = config_container[ZONE_0]->GetnPowerMethodIter(); + su2double tol = config_container[ZONE_0]->GetPowerMethodTolerance(); + int k = 0; - // 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; + 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); + } } - } - // 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) { + // 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: " << global_rho_new << " after " << iter << " iterations." << std::endl; + std::cout << "Spectral Radius Estimate for Iteration " << k << " is " << rho_k << 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; + + // 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; } - } -} \ No newline at end of file + + // 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); +} From 748dd9f5ccb382d555c16f339616ebd5585bfbb4 Mon Sep 17 00:00:00 2001 From: Bot-Enigma-0 Date: Sun, 21 Sep 2025 10:03:38 +0100 Subject: [PATCH 6/8] move to FluidIteration --- SU2_CFD/include/iteration/CFluidIteration.hpp | 12 + SU2_CFD/include/solvers/CSolver.hpp | 3 - SU2_CFD/src/drivers/CSinglezoneDriver.cpp | 214 +----------- SU2_CFD/src/iteration/CFluidIteration.cpp | 309 ++++++++++++++++++ 4 files changed, 323 insertions(+), 215 deletions(-) diff --git a/SU2_CFD/include/iteration/CFluidIteration.hpp b/SU2_CFD/include/iteration/CFluidIteration.hpp index 883540f65ac0..fec794f5a40a 100644 --- a/SU2_CFD/include/iteration/CFluidIteration.hpp +++ b/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& derivatives, CGeometry**** geometry, CSolver***** solver, + unsigned short val_iZone, unsigned short val_iInst); + void GetAllDerivatives(su2matrix& 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 v_estimate; /*! * \brief Imposes a gust via the grid velocities. diff --git a/SU2_CFD/include/solvers/CSolver.hpp b/SU2_CFD/include/solvers/CSolver.hpp index 811846d11479..55e462899324 100644 --- a/SU2_CFD/include/solvers/CSolver.hpp +++ b/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 seed_vector; // Current eigenvector estimate - su2double spectral_radius; - /*--- Variables that need to go. ---*/ su2double *Residual, /*!< \brief Auxiliary nVar vector. */ diff --git a/SU2_CFD/src/drivers/CSinglezoneDriver.cpp b/SU2_CFD/src/drivers/CSinglezoneDriver.cpp index 8ff35fe38c1e..1270eef237ab 100644 --- a/SU2_CFD/src/drivers/CSinglezoneDriver.cpp +++ b/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& 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& 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(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 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); -} diff --git a/SU2_CFD/src/iteration/CFluidIteration.cpp b/SU2_CFD/src/iteration/CFluidIteration.cpp index 529b3e12453c..ab5afabf2b7f 100644 --- a/SU2_CFD/src/iteration/CFluidIteration.cpp +++ b/SU2_CFD/src/iteration/CFluidIteration.cpp @@ -398,6 +398,11 @@ void CFluidIteration::Solve(COutput* output, CIntegration**** integration, CGeom 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 @@ void CFluidIteration::Solve(COutput* output, CIntegration**** integration, CGeom 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,296 @@ void CFluidIteration::SetDualTime_Aeroelastic(CConfig* config) const { } } + +void CFluidIteration::SeedAllDerivatives(const su2matrix& derivatives, CGeometry**** geometry, + CSolver***** solver, unsigned short val_iZone, unsigned short val_iInst) { + unsigned long varOffset = 0; + + // Iterate through all solvers + for (auto iSol = 0u; iSol < MAX_SOLS; ++iSol) { + CSolver* solver_ptr = solver[val_iZone][val_iInst][MESH_0][iSol]; + if (!solver_ptr) { + if (rank == MASTER_NODE) { + std::cout << "Solver " << iSol << " is not allocated, skipping." << std::endl; + } + continue; + } + + CVariable* nodes = solver_ptr->GetNodes(); + + unsigned short nVar = solver_ptr->GetnVar(); + unsigned long nPoint = geometry[val_iZone][val_iInst][MESH_0]->GetnPointDomain(); + + if (rank == MASTER_NODE) { + std::cout << "Solver " << iSol << " has nVar=" << nVar << ", nPoint=" << nPoint << std::endl; + } + + // Check if we have enough columns in the derivatives matrix + if (varOffset + nVar > derivatives.cols()) { + if (rank == MASTER_NODE) { + std::cout << "ERROR: Derivatives matrix doesn't have enough columns for solver " << iSol + << ". Needed: " << varOffset + nVar << ", Available: " << derivatives.cols() << std::endl; + } + varOffset += nVar; + continue; + } + + // Seed derivatives + 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(iPoint, iVar + varOffset); + SU2_TYPE::SetDerivative(solution[iVar], derivative_value); + + } + } + + //variable offset for the next solver + varOffset += nVar; + } + + if (rank == MASTER_NODE) { + std::cout << "Finished SeedAllDerivatives" << std::endl; + } +} + +void CFluidIteration::GetAllDerivatives(su2matrix& derivatives, CGeometry**** geometry, + CSolver***** solver, unsigned short val_iZone, unsigned short val_iInst) { + unsigned long varOffset = 0; + + + //iterate through all solvers + for (auto iSol = 0u; iSol < MAX_SOLS; ++iSol) { + CSolver* solver_ptr = solver[val_iZone][val_iInst][MESH_0][iSol]; + if (!solver_ptr) { + if (rank == MASTER_NODE) { + std::cout << "Solver " << iSol << " is not allocated, skipping." << std::endl; + } + continue; + } + + CVariable* nodes = solver_ptr->GetNodes(); + + unsigned short nVar = solver_ptr->GetnVar(); + unsigned long nPoint = geometry[val_iZone][val_iInst][MESH_0]->GetnPointDomain(); + + //extract derivatives for this 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(iPoint, iVar + varOffset) = derivative_value; + } + } + + //variable offset for the next solver + varOffset += nVar; + } + + if (rank == MASTER_NODE) { + std::cout << "Finished GetAllDerivatives" << 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 nVarTotal = 0; + unsigned long nPoint = geom_ptr->GetnPointDomain(); + + 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(); + } + + //initialise power iteration matrix if needed + if (v_estimate.rows() != nPoint || v_estimate.cols() != nVarTotal) { + v_estimate.resize(nPoint, nVarTotal); + + //assign all ones + for (unsigned long iPoint = 0; iPoint < nPoint; iPoint++) { + for (unsigned long iVar = 0; iVar < nVarTotal; iVar++) { + v_estimate(iPoint, iVar) = static_cast(1.0); + } + } + + //calculate norm + passivedouble local_norm = 0.0; + for (unsigned long iPoint = 0; iPoint < nPoint; iPoint++) { + for (unsigned long iVar = 0; iVar < nVarTotal; iVar++) { + local_norm += v_estimate(iPoint, iVar) * v_estimate(iPoint, iVar); + } + } + + 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; + + if (rank == MASTER_NODE) { + std::cout << "Running PostRunSpectralRadius" << std::endl; + } + + 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 w_k(nPoint, nVarTotal); + + 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> 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 recent_eigenvalues; + + //write matrix to CSV + auto write_matrix = [](const su2matrix& 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>& 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++) { + for (unsigned long iVar = 0; iVar < nVarTotal; iVar++) { + rho_local += w_k(iPoint, iVar) * w_k(iPoint, iVar); + } + } + + 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) { + for (unsigned long iVar = 0; iVar < nVarTotal; ++iVar) { + v_estimate(iPoint, iVar) = w_k(iPoint, iVar) / 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); +} \ No newline at end of file From 0c035c453e0e9b8dda4e60fb458ca6ff486fd3d5 Mon Sep 17 00:00:00 2001 From: Bot-Enigma-0 Date: Sun, 21 Sep 2025 10:06:17 +0100 Subject: [PATCH 7/8] revert singlezone.hpp --- SU2_CFD/include/drivers/CSinglezoneDriver.hpp | 19 ------------------- 1 file changed, 19 deletions(-) diff --git a/SU2_CFD/include/drivers/CSinglezoneDriver.hpp b/SU2_CFD/include/drivers/CSinglezoneDriver.hpp index c697ff67915c..355369556ea5 100644 --- a/SU2_CFD/include/drivers/CSinglezoneDriver.hpp +++ b/SU2_CFD/include/drivers/CSinglezoneDriver.hpp @@ -46,19 +46,6 @@ class CSinglezoneDriver : public CDriver { * \return Boolean indicating whether the problem is converged. */ virtual bool GetTimeConvergence() const; - su2matrix v_estimate; - - /*! - * \brief Seed derivatives for all solvers in the zone. - * \param[in] derivatives - Matrix of derivative values - */ - void SeedAllDerivatives(const su2matrix& derivatives, CGeometry *geometry); - - /*! - * \brief Extract derivatives from all solvers in the zone. - * \param[out] derivatives - Matrix to store derivative values - */ - void GetAllDerivatives(su2matrix& derivatives, CGeometry *geometry); public: @@ -123,10 +110,4 @@ class CSinglezoneDriver : public CDriver { */ bool Monitor(unsigned long TimeIter) override; - /*! - * \brief Calculate the spectral radius using power iteration method. - */ - void PreRunSpectralRadius(); - - void PostRunSpectralRadius(); }; From 60f62ebd957a133b3d061bee92933c47af381b05 Mon Sep 17 00:00:00 2001 From: Bot-Enigma-0 Date: Sun, 21 Sep 2025 20:17:29 +0100 Subject: [PATCH 8/8] seed single variable --- SU2_CFD/src/iteration/CFluidIteration.cpp | 152 +++++++--------------- 1 file changed, 47 insertions(+), 105 deletions(-) diff --git a/SU2_CFD/src/iteration/CFluidIteration.cpp b/SU2_CFD/src/iteration/CFluidIteration.cpp index ab5afabf2b7f..e6ccf83605c4 100644 --- a/SU2_CFD/src/iteration/CFluidIteration.cpp +++ b/SU2_CFD/src/iteration/CFluidIteration.cpp @@ -718,96 +718,53 @@ void CFluidIteration::SetDualTime_Aeroelastic(CConfig* config) const { void CFluidIteration::SeedAllDerivatives(const su2matrix& derivatives, CGeometry**** geometry, CSolver***** solver, unsigned short val_iZone, unsigned short val_iInst) { - unsigned long varOffset = 0; + unsigned short targetSolver = 2; + unsigned short targetVar = 0; - // Iterate through all solvers - for (auto iSol = 0u; iSol < MAX_SOLS; ++iSol) { - CSolver* solver_ptr = solver[val_iZone][val_iInst][MESH_0][iSol]; - if (!solver_ptr) { - if (rank == MASTER_NODE) { - std::cout << "Solver " << iSol << " is not allocated, skipping." << std::endl; - } - continue; - } - - CVariable* nodes = solver_ptr->GetNodes(); - - unsigned short nVar = solver_ptr->GetnVar(); - unsigned long nPoint = geometry[val_iZone][val_iInst][MESH_0]->GetnPointDomain(); - - if (rank == MASTER_NODE) { - std::cout << "Solver " << iSol << " has nVar=" << nVar << ", nPoint=" << nPoint << std::endl; - } - - // Check if we have enough columns in the derivatives matrix - if (varOffset + nVar > derivatives.cols()) { - if (rank == MASTER_NODE) { - std::cout << "ERROR: Derivatives matrix doesn't have enough columns for solver " << iSol - << ". Needed: " << varOffset + nVar << ", Available: " << derivatives.cols() << std::endl; - } - varOffset += nVar; - continue; - } - - // Seed derivatives - 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(iPoint, iVar + varOffset); - SU2_TYPE::SetDerivative(solution[iVar], derivative_value); - - } - } - - //variable offset for the next solver - varOffset += nVar; + // 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" << std::endl; + std::cout << "Finished SeedAllDerivatives for solver " << targetSolver + << ", variable " << targetVar << std::endl; } } void CFluidIteration::GetAllDerivatives(su2matrix& derivatives, CGeometry**** geometry, CSolver***** solver, unsigned short val_iZone, unsigned short val_iInst) { - unsigned long varOffset = 0; + 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; + } - //iterate through all solvers - for (auto iSol = 0u; iSol < MAX_SOLS; ++iSol) { - CSolver* solver_ptr = solver[val_iZone][val_iInst][MESH_0][iSol]; - if (!solver_ptr) { - if (rank == MASTER_NODE) { - std::cout << "Solver " << iSol << " is not allocated, skipping." << std::endl; - } - continue; - } - - CVariable* nodes = solver_ptr->GetNodes(); - - unsigned short nVar = solver_ptr->GetnVar(); - unsigned long nPoint = geometry[val_iZone][val_iInst][MESH_0]->GetnPointDomain(); - - //extract derivatives for this 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(iPoint, iVar + varOffset) = derivative_value; - } - } - - //variable offset for the next solver - varOffset += nVar; + 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" << std::endl; + std::cout << "Finished GetAllDerivatives for solver " << target_Solver + << ", variable " << target_Var << std::endl; } } @@ -818,31 +775,25 @@ void CFluidIteration::PreRunSpectralRadius(CGeometry**** geometry, CSolver***** CGeometry* geom_ptr = geometry[val_iZone][val_iInst][MESH_0]; //get total number of variables and points - unsigned long nVarTotal = 0; unsigned long nPoint = geom_ptr->GetnPointDomain(); - 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(); - } - - //initialise power iteration matrix if needed + //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); - //assign all ones for (unsigned long iPoint = 0; iPoint < nPoint; iPoint++) { - for (unsigned long iVar = 0; iVar < nVarTotal; iVar++) { - v_estimate(iPoint, iVar) = static_cast(1.0); - } + v_estimate(iPoint, 0) = static_cast(1.0); } - + //calculate norm passivedouble local_norm = 0.0; for (unsigned long iPoint = 0; iPoint < nPoint; iPoint++) { - for (unsigned long iVar = 0; iVar < nVarTotal; iVar++) { - local_norm += v_estimate(iPoint, iVar) * v_estimate(iPoint, iVar); - } + local_norm += v_estimate(iPoint, 0) * v_estimate(iPoint, 0); } passivedouble global_norm; @@ -857,7 +808,6 @@ void CFluidIteration::PreRunSpectralRadius(CGeometry**** geometry, CSolver***** } } - //seed all solvers SeedAllDerivatives(v_estimate, geometry, solver, val_iZone, val_iInst); @@ -872,10 +822,6 @@ void CFluidIteration::PostRunSpectralRadius(COutput* output, CIntegration**** in unsigned short val_iInst) { if (!config[val_iZone]->GetSpectralRadius_Analysis()) return; - if (rank == MASTER_NODE) { - std::cout << "Running PostRunSpectralRadius" << std::endl; - } - CGeometry* geom_ptr = geometry[val_iZone][val_iInst][MESH_0]; unsigned long nPoint = geom_ptr->GetnPointDomain(); @@ -890,7 +836,7 @@ void CFluidIteration::PostRunSpectralRadius(COutput* output, CIntegration**** in std::cout << "Total variables: " << nVarTotal << ", Total points: " << nPoint << std::endl; } - su2matrix w_k(nPoint, nVarTotal); + su2matrix w_k(nPoint, 1); if (rank == MASTER_NODE) { std::cout << "Matrix w initialized with size: " << w_k.rows() << " x " << w_k.cols() << std::endl; @@ -942,9 +888,7 @@ void CFluidIteration::PostRunSpectralRadius(COutput* output, CIntegration**** in GetAllDerivatives(w_k, geometry, solver, val_iZone, val_iInst); passivedouble rho_local = 0.0; for (unsigned long iPoint = 0; iPoint < nPoint; iPoint++) { - for (unsigned long iVar = 0; iVar < nVarTotal; iVar++) { - rho_local += w_k(iPoint, iVar) * w_k(iPoint, iVar); - } + rho_local += w_k(iPoint, 0) * w_k(iPoint, 0); // Only column 0 } passivedouble rho_global = 0.0; @@ -980,10 +924,8 @@ void CFluidIteration::PostRunSpectralRadius(COutput* output, CIntegration**** in } //normalize w_k to form next seed v_estimate - for (unsigned long iPoint = 0; iPoint < nPoint; ++iPoint) { - for (unsigned long iVar = 0; iVar < nVarTotal; ++iVar) { - v_estimate(iPoint, iVar) = w_k(iPoint, iVar) / rho_k; - } + for (unsigned long iPoint = 0; iPoint < nPoint; iPoint++) { + v_estimate(iPoint, 0) = w_k(iPoint, 0) / rho_k; } //seed the derivatives for this iteration @@ -1007,4 +949,4 @@ void CFluidIteration::PostRunSpectralRadius(COutput* output, CIntegration**** in } w_k.resize(0, 0); -} \ No newline at end of file +}