From 3a0b257fe1e420059d23dc78ac6fb4370e98bd3a Mon Sep 17 00:00:00 2001 From: Vaishnavi Wani Date: Thu, 26 Dec 2024 01:42:52 +0100 Subject: [PATCH 1/6] Instantiation of solver in CSolverFactory and add a call to the SingleGrid_Iteration in CFEAIteration --- SU2_CFD/src/iteration/CFEAIteration.cpp | 7 +++++++ SU2_CFD/src/solvers/CSolverFactory.cpp | 3 +++ 2 files changed, 10 insertions(+) diff --git a/SU2_CFD/src/iteration/CFEAIteration.cpp b/SU2_CFD/src/iteration/CFEAIteration.cpp index 2aa19d27de9..1b38dc580bb 100644 --- a/SU2_CFD/src/iteration/CFEAIteration.cpp +++ b/SU2_CFD/src/iteration/CFEAIteration.cpp @@ -48,6 +48,13 @@ void CFEAIteration::Iterate(COutput* output, CIntegration**** integration, CGeom CIntegration* feaIntegration = integration[val_iZone][val_iInst][FEA_SOL]; CSolver* feaSolver = solver[val_iZone][val_iInst][MESH_0][FEA_SOL]; + /*--- Add heat solver integration step ---*/ + if (config[val_iZone]->GetWeakly_Coupled_Heat()) { + config[val_iZone]->SetGlobalParam(MAIN_SOLVER::HEAT_EQUATION, RUNTIME_HEAT_SYS); + integration[val_iZone][val_iInst][HEAT_SOL]->SingleGrid_Iteration( + geometry, solver, numerics, config, RUNTIME_HEAT_SYS, val_iZone, val_iInst); + } + /*--- FEA equations ---*/ config[val_iZone]->SetGlobalParam(MAIN_SOLVER::FEM_ELASTICITY, RUNTIME_FEA_SYS); diff --git a/SU2_CFD/src/solvers/CSolverFactory.cpp b/SU2_CFD/src/solvers/CSolverFactory.cpp index 743775ad139..7c055241624 100644 --- a/SU2_CFD/src/solvers/CSolverFactory.cpp +++ b/SU2_CFD/src/solvers/CSolverFactory.cpp @@ -173,6 +173,9 @@ CSolver** CSolverFactory::CreateSolverContainer(MAIN_SOLVER kindMainSolver, CCon break; case MAIN_SOLVER::FEM_ELASTICITY: solver[FEA_SOL] = CreateSubSolver(SUB_SOLVER_TYPE::FEA, solver, geometry, config, iMGLevel); + if (config->GetWeakly_Coupled_Heat()) { + solver[HEAT_SOL] = CreateSubSolver(SUB_SOLVER_TYPE::HEAT, solver, geometry, config, iMGLevel); + } break; case MAIN_SOLVER::DISC_ADJ_FEM: solver[FEA_SOL] = CreateSubSolver(SUB_SOLVER_TYPE::FEA, solver, geometry, config, iMGLevel); From 323c27ba3f7ead35603694d0a1685bf6103c4582 Mon Sep 17 00:00:00 2001 From: Vaish-W Date: Thu, 2 Jan 2025 20:37:38 +0100 Subject: [PATCH 2/6] Add nodal temperature integration in Compute_StiffMatrix --- SU2_CFD/include/solvers/CFEASolver.hpp | 17 +++++++++++++++++ SU2_CFD/src/solvers/CFEASolver.cpp | 9 +++++++++ Tutorials | 1 + 3 files changed, 27 insertions(+) create mode 160000 Tutorials diff --git a/SU2_CFD/include/solvers/CFEASolver.hpp b/SU2_CFD/include/solvers/CFEASolver.hpp index 4d8bf74be8e..16019ff5403 100644 --- a/SU2_CFD/include/solvers/CFEASolver.hpp +++ b/SU2_CFD/include/solvers/CFEASolver.hpp @@ -28,6 +28,7 @@ #pragma once #include "CFEASolverBase.hpp" +#include "CHeatSolver.hpp" /*! * \class CFEASolver @@ -99,6 +100,22 @@ class CFEASolver : public CFEASolverBase { bool topol_filter_applied; /*!< \brief True if density filtering has been performed. */ bool initial_calc = true; /*!< \brief Becomes false after first call to Preprocessing. */ + /*! + * \brief Pointer to the heat solver for coupled simulations. + * + * This member stores a pointer to the heat solver, which handles + * the solution of the heat equation in weakly coupled simulations. + * It is initialized during the preprocessing step if the configuration + * enables the weak coupling of heat and elasticity solvers. This solver + * provides temperature information to the finite element elasticity solver + * and contributes to the coupled residuals. + * + * The `heat_solver` pointer remains `nullptr` when the heat solver is not enabled + * in the configuration. Memory for the heat solver is dynamically allocated + * during initialization and released in the destructor to avoid memory leaks. + */ + CSolver* heat_solver = nullptr; + /*! * \brief The highest level in the variable hierarchy this solver can safely use, * CVariable is the common denominator between the FEA and Mesh deformation variables. diff --git a/SU2_CFD/src/solvers/CFEASolver.cpp b/SU2_CFD/src/solvers/CFEASolver.cpp index 44a19e1421a..267c0e43b64 100644 --- a/SU2_CFD/src/solvers/CFEASolver.cpp +++ b/SU2_CFD/src/solvers/CFEASolver.cpp @@ -30,6 +30,7 @@ #include "../../include/numerics/elasticity/CFEAElasticity.hpp" #include "../../../Common/include/toolboxes/printing_toolbox.hpp" #include "../../../Common/include/toolboxes/geometry_toolbox.hpp" +#include "../../include/solvers/CHeatSolver.hpp" #include using namespace GeometryToolbox; @@ -570,6 +571,10 @@ void CFEASolver::Preprocessing(CGeometry *geometry, CSolver **solver_container, const bool body_forces = config->GetDeadLoad(); const bool topology_mode = config->GetTopology_Optimization(); + if (config->GetWeakly_Coupled_Heat()) { + heat_solver = solver_container[HEAT_SOL]; + } + /* * For topology optimization we apply a filter on the design density field to avoid * numerical issues (checkerboards), ensure mesh independence, and impose a length scale. @@ -694,6 +699,10 @@ void CFEASolver::Compute_StiffMatrix(CGeometry *geometry, CNumerics **numerics, su2double val_Sol = nodes->GetSolution(indexNode[iNode],iDim) + val_Coord; element->SetRef_Coord(iNode, iDim, val_Coord); element->SetCurr_Coord(iNode, iDim, val_Sol); + } + if (heat_solver) { + auto nodal_temperatures = dynamic_cast(heat_solver)->GetNodalTemperature(); + element->SetTemperature(iNode, nodal_temperatures[indexNode[iNode]]); } } diff --git a/Tutorials b/Tutorials new file mode 160000 index 00000000000..e8d5cdbb960 --- /dev/null +++ b/Tutorials @@ -0,0 +1 @@ +Subproject commit e8d5cdbb96093e659f268b15cfe4a32572335256 From 60d7b6c94688db07f4a4116ef9d6b392641f7a48 Mon Sep 17 00:00:00 2001 From: Vaish-W Date: Fri, 3 Jan 2025 11:55:23 +0100 Subject: [PATCH 3/6] Handle temperature from the heat_solver --- SU2_CFD/include/solvers/CFEASolver.hpp | 1 - SU2_CFD/src/solvers/CFEASolver.cpp | 20 ++++++++++++++++++-- 2 files changed, 18 insertions(+), 3 deletions(-) diff --git a/SU2_CFD/include/solvers/CFEASolver.hpp b/SU2_CFD/include/solvers/CFEASolver.hpp index c003373f639..a7cc6901418 100644 --- a/SU2_CFD/include/solvers/CFEASolver.hpp +++ b/SU2_CFD/include/solvers/CFEASolver.hpp @@ -28,7 +28,6 @@ #pragma once #include "CFEASolverBase.hpp" -#include "CHeatSolver.hpp" /*! * \class CFEASolver diff --git a/SU2_CFD/src/solvers/CFEASolver.cpp b/SU2_CFD/src/solvers/CFEASolver.cpp index 520fb1e21c2..886da26df69 100644 --- a/SU2_CFD/src/solvers/CFEASolver.cpp +++ b/SU2_CFD/src/solvers/CFEASolver.cpp @@ -694,8 +694,8 @@ void CFEASolver::Compute_StiffMatrix(CGeometry *geometry, CNumerics **numerics, element->SetCurr_Coord(iNode, iDim, val_Sol); } if (heat_solver) { - auto nodal_temperatures = dynamic_cast(heat_solver)->GetNodalTemperature(); - element->SetTemperature(iNode, nodal_temperatures[indexNode[iNode]]); + const su2double temperature = heat_solver->GetNodes()->GetSolution(indexNode[iNode], 0); + element->SetTemperature(iNode, temperature); } } @@ -804,6 +804,14 @@ void CFEASolver::Compute_StiffMatrix_NodalStressRes(CGeometry *geometry, CNumeri de_elem->SetRef_Coord(iNode, iDim, val_Coord); } } + if (heat_solver) { + const su2double temperature = heat_solver->GetNodes()->GetSolution(indexNode[iNode], 0); + fea_elem->SetTemperature(iNode, temperature); + + if (de_effects) { + de_elem->SetTemperature(iNode, temperature); + } + } } /*--- In topology mode determine the penalty to apply to the stiffness. ---*/ @@ -1102,6 +1110,10 @@ void CFEASolver::Compute_NodalStressRes(CGeometry *geometry, CNumerics **numeric element->SetCurr_Coord(iNode, iDim, val_Sol); element->SetRef_Coord(iNode, iDim, val_Coord); } + if (heat_solver) { + const su2double temperature = heat_solver->GetNodes()->GetSolution(indexNode[iNode], 0); + element->SetTemperature(iNode, temperature); + } } /*--- In topology mode determine the penalty to apply to the stiffness ---*/ @@ -1218,6 +1230,10 @@ void CFEASolver::Compute_NodalStress(CGeometry *geometry, CNumerics **numerics, element->SetCurr_Coord(iNode, iDim, val_Sol); element->SetRef_Coord(iNode, iDim, val_Coord); } + if (heat_solver) { + const su2double temperature = heat_solver->GetNodes()->GetSolution(indexNode[iNode], 0); + element->SetTemperature(iNode, temperature); + } } /*--- In topology mode determine the penalty to apply to the stiffness ---*/ From b1499dda5c377b794bea0f62cae985019a471d7b Mon Sep 17 00:00:00 2001 From: Vaish-W Date: Fri, 3 Jan 2025 13:27:36 +0100 Subject: [PATCH 4/6] Output solution for the heat solver in CElasticityOutput --- SU2_CFD/include/output/CElasticityOutput.hpp | 9 +++++++ SU2_CFD/src/output/CElasticityOutput.cpp | 26 ++++++++++++++++++++ 2 files changed, 35 insertions(+) diff --git a/SU2_CFD/include/output/CElasticityOutput.hpp b/SU2_CFD/include/output/CElasticityOutput.hpp index ca65cab7a9c..199982ff77b 100644 --- a/SU2_CFD/include/output/CElasticityOutput.hpp +++ b/SU2_CFD/include/output/CElasticityOutput.hpp @@ -37,6 +37,7 @@ class CElasticityOutput final: public COutput { protected: + const CSolver *heat_solver = nullptr; //!< Pointer to the heat solver unsigned short nVar_FEM; //!< Number of FEM variables bool linear_analysis, //!< Boolean indicating a linear analysis nonlinear_analysis, //!< Boolean indicating a nonlinear analysis @@ -84,4 +85,12 @@ class CElasticityOutput final: public COutput { */ bool SetInitResiduals(const CConfig *config) override ; + /*! + * \brief Set the heat solver pointer for output integration. + * \param[in] solver - Pointer to the heat solver. + */ + void SetHeatSolver(const CSolver *solver) { + heat_solver = solver; + } + }; diff --git a/SU2_CFD/src/output/CElasticityOutput.cpp b/SU2_CFD/src/output/CElasticityOutput.cpp index a2abc4e6d12..79ecd22f305 100644 --- a/SU2_CFD/src/output/CElasticityOutput.cpp +++ b/SU2_CFD/src/output/CElasticityOutput.cpp @@ -152,6 +152,15 @@ void CElasticityOutput::LoadHistoryData(CConfig *config, CGeometry *geometry, CS /*--- Keep this as last, since it uses the history values that were set. ---*/ SetCustomAndComboObjectives(FEA_SOL, config, solver); + /*--- Add heat solver data if available ---*/ + if (heat_solver) { + SetHistoryOutputValue("RMS_TEMPERATURE", log10(heat_solver->GetRes_RMS(0))); + SetHistoryOutputValue("MAX_TEMPERATURE", log10(heat_solver->GetRes_Max(0))); + SetHistoryOutputValue("AVG_TEMPERATURE", heat_solver->GetTotal_AvgTemperature()); + SetHistoryOutputValue("TOTAL_HEATFLUX", heat_solver->GetTotal_HeatFlux()); + SetHistoryOutputValue("MAXIMUM_HEATFLUX", heat_solver->GetTotal_MaxHeatFlux()); + } + } void CElasticityOutput::SetHistoryOutputFields(CConfig *config) { @@ -189,6 +198,12 @@ void CElasticityOutput::SetHistoryOutputFields(CConfig *config) { } AddHistoryOutput("COMBO", "ComboObj", ScreenOutputFormat::SCIENTIFIC, "COMBO", "Combined obj. function value.", HistoryFieldType::COEFFICIENT); + AddHistoryOutput("RMS_TEMPERATURE", "rms[T]", ScreenOutputFormat::FIXED, "RMS_RES", "Root mean square residual of the temperature", HistoryFieldType::RESIDUAL); + AddHistoryOutput("MAX_TEMPERATURE", "max[T]", ScreenOutputFormat::FIXED, "MAX_RES", "Maximum residual of the temperature", HistoryFieldType::RESIDUAL); + AddHistoryOutput("AVG_TEMPERATURE", "avg[T]", ScreenOutputFormat::SCIENTIFIC, "HEAT", "Average temperature on all surfaces", HistoryFieldType::COEFFICIENT); + AddHistoryOutput("TOTAL_HEATFLUX", "HF", ScreenOutputFormat::SCIENTIFIC, "HEAT", "Total heat flux on all surfaces", HistoryFieldType::COEFFICIENT); + AddHistoryOutput("MAXIMUM_HEATFLUX", "MaxHF", ScreenOutputFormat::SCIENTIFIC, "HEAT", "Maximum heat flux on all surfaces", HistoryFieldType::COEFFICIENT); + } void CElasticityOutput::LoadVolumeData(CConfig *config, CGeometry *geometry, CSolver **solver, unsigned long iPoint){ @@ -228,6 +243,13 @@ void CElasticityOutput::LoadVolumeData(CConfig *config, CGeometry *geometry, CSo if (config->GetTopology_Optimization()) { SetVolumeOutputValue("TOPOL_DENSITY", iPoint, Node_Struc->GetAuxVar(iPoint)); } + + if (heat_solver) { + const auto Node_Heat = heat_solver->GetNodes(); + SetVolumeOutputValue("TEMPERATURE", iPoint, Node_Heat->GetSolution(iPoint, 0)); + SetVolumeOutputValue("RES_TEMPERATURE", iPoint, heat_solver->LinSysRes(iPoint, 0)); + } + } void CElasticityOutput::SetVolumeOutputFields(CConfig *config){ @@ -267,6 +289,10 @@ void CElasticityOutput::SetVolumeOutputFields(CConfig *config){ if (config->GetTopology_Optimization()) { AddVolumeOutput("TOPOL_DENSITY", "Topology_Density", "TOPOLOGY", "filtered topology density"); } + + AddVolumeOutput("TEMPERATURE", "Temperature", "SOLUTION", "Temperature"); + AddVolumeOutput("RES_TEMPERATURE", "Residual_Temperature", "RESIDUAL", "Residual of the temperature"); + } bool CElasticityOutput::SetInitResiduals(const CConfig *config){ From 0298a31b2681a71104b40422b5a8bb61133a64b0 Mon Sep 17 00:00:00 2001 From: Vaish-W Date: Fri, 3 Jan 2025 18:25:26 +0100 Subject: [PATCH 5/6] Remove Tutorials folder from the repository --- Tutorials | 1 - 1 file changed, 1 deletion(-) delete mode 160000 Tutorials diff --git a/Tutorials b/Tutorials deleted file mode 160000 index e8d5cdbb960..00000000000 --- a/Tutorials +++ /dev/null @@ -1 +0,0 @@ -Subproject commit e8d5cdbb96093e659f268b15cfe4a32572335256 From 3c83bbef12b0be43fca0a9d03b3e8902395cf76b Mon Sep 17 00:00:00 2001 From: Vaish-W Date: Fri, 3 Jan 2025 18:30:12 +0100 Subject: [PATCH 6/6] Add a check for the heat solver --- SU2_CFD/include/output/CElasticityOutput.hpp | 9 --------- SU2_CFD/src/output/CElasticityOutput.cpp | 20 +++++++++++++------- SU2_CFD/src/solvers/CFEASolver.cpp | 4 ---- 3 files changed, 13 insertions(+), 20 deletions(-) diff --git a/SU2_CFD/include/output/CElasticityOutput.hpp b/SU2_CFD/include/output/CElasticityOutput.hpp index 199982ff77b..ca65cab7a9c 100644 --- a/SU2_CFD/include/output/CElasticityOutput.hpp +++ b/SU2_CFD/include/output/CElasticityOutput.hpp @@ -37,7 +37,6 @@ class CElasticityOutput final: public COutput { protected: - const CSolver *heat_solver = nullptr; //!< Pointer to the heat solver unsigned short nVar_FEM; //!< Number of FEM variables bool linear_analysis, //!< Boolean indicating a linear analysis nonlinear_analysis, //!< Boolean indicating a nonlinear analysis @@ -85,12 +84,4 @@ class CElasticityOutput final: public COutput { */ bool SetInitResiduals(const CConfig *config) override ; - /*! - * \brief Set the heat solver pointer for output integration. - * \param[in] solver - Pointer to the heat solver. - */ - void SetHeatSolver(const CSolver *solver) { - heat_solver = solver; - } - }; diff --git a/SU2_CFD/src/output/CElasticityOutput.cpp b/SU2_CFD/src/output/CElasticityOutput.cpp index 79ecd22f305..62815f4fcb0 100644 --- a/SU2_CFD/src/output/CElasticityOutput.cpp +++ b/SU2_CFD/src/output/CElasticityOutput.cpp @@ -106,6 +106,7 @@ CElasticityOutput::CElasticityOutput(CConfig *config, unsigned short nDim) : COu void CElasticityOutput::LoadHistoryData(CConfig *config, CGeometry *geometry, CSolver **solver) { CSolver* fea_solver = solver[FEA_SOL]; + CSolver* heat_solver = solver[HEAT_SOL]; /*--- Residuals: ---*/ /*--- Linear analysis: RMS of the displacements in the nDim coordinates ---*/ @@ -198,11 +199,13 @@ void CElasticityOutput::SetHistoryOutputFields(CConfig *config) { } AddHistoryOutput("COMBO", "ComboObj", ScreenOutputFormat::SCIENTIFIC, "COMBO", "Combined obj. function value.", HistoryFieldType::COEFFICIENT); - AddHistoryOutput("RMS_TEMPERATURE", "rms[T]", ScreenOutputFormat::FIXED, "RMS_RES", "Root mean square residual of the temperature", HistoryFieldType::RESIDUAL); - AddHistoryOutput("MAX_TEMPERATURE", "max[T]", ScreenOutputFormat::FIXED, "MAX_RES", "Maximum residual of the temperature", HistoryFieldType::RESIDUAL); - AddHistoryOutput("AVG_TEMPERATURE", "avg[T]", ScreenOutputFormat::SCIENTIFIC, "HEAT", "Average temperature on all surfaces", HistoryFieldType::COEFFICIENT); - AddHistoryOutput("TOTAL_HEATFLUX", "HF", ScreenOutputFormat::SCIENTIFIC, "HEAT", "Total heat flux on all surfaces", HistoryFieldType::COEFFICIENT); - AddHistoryOutput("MAXIMUM_HEATFLUX", "MaxHF", ScreenOutputFormat::SCIENTIFIC, "HEAT", "Maximum heat flux on all surfaces", HistoryFieldType::COEFFICIENT); + if (config->GetWeakly_Coupled_Heat()) { + AddHistoryOutput("RMS_TEMPERATURE", "rms[T]", ScreenOutputFormat::FIXED, "RMS_RES", "Root mean square residual of the temperature", HistoryFieldType::RESIDUAL); + AddHistoryOutput("MAX_TEMPERATURE", "max[T]", ScreenOutputFormat::FIXED, "MAX_RES", "Maximum residual of the temperature", HistoryFieldType::RESIDUAL); + AddHistoryOutput("AVG_TEMPERATURE", "avg[T]", ScreenOutputFormat::SCIENTIFIC, "HEAT", "Average temperature on all surfaces", HistoryFieldType::COEFFICIENT); + AddHistoryOutput("TOTAL_HEATFLUX", "HF", ScreenOutputFormat::SCIENTIFIC, "HEAT", "Total heat flux on all surfaces", HistoryFieldType::COEFFICIENT); + AddHistoryOutput("MAXIMUM_HEATFLUX", "MaxHF", ScreenOutputFormat::SCIENTIFIC, "HEAT", "Maximum heat flux on all surfaces", HistoryFieldType::COEFFICIENT); + } } @@ -244,6 +247,7 @@ void CElasticityOutput::LoadVolumeData(CConfig *config, CGeometry *geometry, CSo SetVolumeOutputValue("TOPOL_DENSITY", iPoint, Node_Struc->GetAuxVar(iPoint)); } + CSolver* heat_solver = solver[HEAT_SOL]; if (heat_solver) { const auto Node_Heat = heat_solver->GetNodes(); SetVolumeOutputValue("TEMPERATURE", iPoint, Node_Heat->GetSolution(iPoint, 0)); @@ -290,8 +294,10 @@ void CElasticityOutput::SetVolumeOutputFields(CConfig *config){ AddVolumeOutput("TOPOL_DENSITY", "Topology_Density", "TOPOLOGY", "filtered topology density"); } - AddVolumeOutput("TEMPERATURE", "Temperature", "SOLUTION", "Temperature"); - AddVolumeOutput("RES_TEMPERATURE", "Residual_Temperature", "RESIDUAL", "Residual of the temperature"); + if (config->GetWeakly_Coupled_Heat()) { + AddVolumeOutput("TEMPERATURE", "Temperature", "SOLUTION", "Temperature"); + AddVolumeOutput("RES_TEMPERATURE", "Residual_Temperature", "RESIDUAL", "Residual of the temperature"); + } } diff --git a/SU2_CFD/src/solvers/CFEASolver.cpp b/SU2_CFD/src/solvers/CFEASolver.cpp index 886da26df69..03f2307e5fd 100644 --- a/SU2_CFD/src/solvers/CFEASolver.cpp +++ b/SU2_CFD/src/solvers/CFEASolver.cpp @@ -807,10 +807,6 @@ void CFEASolver::Compute_StiffMatrix_NodalStressRes(CGeometry *geometry, CNumeri if (heat_solver) { const su2double temperature = heat_solver->GetNodes()->GetSolution(indexNode[iNode], 0); fea_elem->SetTemperature(iNode, temperature); - - if (de_effects) { - de_elem->SetTemperature(iNode, temperature); - } } }