diff --git a/SU2_CFD/include/iteration/CDiscAdjFEAIteration.hpp b/SU2_CFD/include/iteration/CDiscAdjFEAIteration.hpp index f2246f96cef..12b4403f77d 100644 --- a/SU2_CFD/include/iteration/CDiscAdjFEAIteration.hpp +++ b/SU2_CFD/include/iteration/CDiscAdjFEAIteration.hpp @@ -42,6 +42,9 @@ class CDiscAdjFEAIteration final : public CIteration { private: unsigned short CurrentRecording; /*!< \brief Stores the current status of the recording. */ + /*!< \brief Used for thermo elastic problems to avoid duplicating code. */ + CIteration* DiscAdjHeatIteration = nullptr; + /*! * \brief load solution for dynamic problems * \param[in] geometry - Geometrical definition of the problem. @@ -64,7 +67,7 @@ class CDiscAdjFEAIteration final : public CIteration { /*! * \brief Destructor of the class. */ - ~CDiscAdjFEAIteration(void) override; + ~CDiscAdjFEAIteration() override; /*! * \brief Preprocessing to prepare for an iteration of the physics. diff --git a/SU2_CFD/include/output/CAdjElasticityOutput.hpp b/SU2_CFD/include/output/CAdjElasticityOutput.hpp index 8a1d0559ca3..c6872d5b3e2 100644 --- a/SU2_CFD/include/output/CAdjElasticityOutput.hpp +++ b/SU2_CFD/include/output/CAdjElasticityOutput.hpp @@ -36,6 +36,7 @@ */ class CAdjElasticityOutput final : public COutput { private: + bool coupled_heat; //!< Boolean indicating a thermoelastic analysis unsigned short nVar_FEM; //!< Number of FEM variables public: diff --git a/SU2_CFD/include/output/CAdjHeatOutput.hpp b/SU2_CFD/include/output/CAdjHeatOutput.hpp index 72f8cffc931..7c7a5a0bf86 100644 --- a/SU2_CFD/include/output/CAdjHeatOutput.hpp +++ b/SU2_CFD/include/output/CAdjHeatOutput.hpp @@ -46,7 +46,12 @@ class CAdjHeatOutput final: public COutput { /*! * \brief Destructor of the class. */ - ~CAdjHeatOutput(void) override; + ~CAdjHeatOutput() override; + + /*! + * \brief Set the history output field values in another output instance. + */ + static void LoadHistoryDataImpl(CConfig *config, CGeometry *geometry, CSolver **solver, COutput* output); /*! * \brief Load the history output field values @@ -54,6 +59,11 @@ class CAdjHeatOutput final: public COutput { */ void LoadHistoryData(CConfig *config, CGeometry *geometry, CSolver **solver) override; + /*! + * \brief Set the available history output fields in another output instance. + */ + static void SetHistoryOutputFieldsImpl(CConfig *config, COutput* output); + /*! * \brief Set the available history output fields * \param[in] config - Definition of the particular problem. diff --git a/SU2_CFD/include/output/CElasticityOutput.hpp b/SU2_CFD/include/output/CElasticityOutput.hpp index c3c189512b4..ecf054a631a 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: - unsigned short nVar_FEM; //!< Number of FEM variables bool linear_analysis, //!< Boolean indicating a linear analysis nonlinear_analysis, //!< Boolean indicating a nonlinear analysis coupled_heat, //!< Boolean indicating a thermoelastic analysis diff --git a/SU2_CFD/include/output/COutput.hpp b/SU2_CFD/include/output/COutput.hpp index ce4e29ed2c1..555890b29ff 100644 --- a/SU2_CFD/include/output/COutput.hpp +++ b/SU2_CFD/include/output/COutput.hpp @@ -56,6 +56,7 @@ class CFileWriter; class CParallelDataSorter; class CConfig; class CHeatOutput; +class CAdjHeatOutput; using namespace std; @@ -67,6 +68,7 @@ using namespace std; class COutput { protected: friend class CHeatOutput; + friend class CAdjHeatOutput; /*----------------------------- General ----------------------------*/ diff --git a/SU2_CFD/include/solvers/CFEASolver.hpp b/SU2_CFD/include/solvers/CFEASolver.hpp index ce472e58575..631f729b04e 100644 --- a/SU2_CFD/include/solvers/CFEASolver.hpp +++ b/SU2_CFD/include/solvers/CFEASolver.hpp @@ -550,7 +550,7 @@ class CFEASolver : public CFEASolverBase { * \param[in] config - Definition of the problem. * \param[in] solver - Container vector with all the solutions. */ - void Evaluate_ObjFunc(const CConfig *config, CSolver**) final { + void Evaluate_ObjFunc(const CConfig *config, CSolver** solver) final { Total_ComboObj = 0.0; switch (config->GetKind_ObjFunc()) { case REFERENCE_GEOMETRY: @@ -575,6 +575,10 @@ class CFEASolver : public CFEASolverBase { Total_ComboObj = Total_Custom_ObjFunc; break; } + if (config->GetWeakly_Coupled_Heat()) { + solver[HEAT_SOL]->Evaluate_ObjFunc(config, solver); + Total_ComboObj += solver[HEAT_SOL]->GetTotal_ComboObj(); + } } /*! diff --git a/SU2_CFD/src/drivers/CDriver.cpp b/SU2_CFD/src/drivers/CDriver.cpp index 961ce04832f..f491caa39b3 100644 --- a/SU2_CFD/src/drivers/CDriver.cpp +++ b/SU2_CFD/src/drivers/CDriver.cpp @@ -2769,9 +2769,10 @@ void CDriver::PrintDirectResidual(RECORDING kind_recording) { /*--- Helper lambda func to return lenghty [iVar][iZone] string. ---*/ auto iVar_iZone2string = [&](unsigned short ivar, unsigned short izone) { - if (multizone) + if (multizone) { return "[" + std::to_string(ivar) + "][" + std::to_string(izone) + "]"; - return "[" + std::to_string(ivar) + "]"; + } + return "[" + std::to_string(ivar) + "]"; }; /*--- Print residuals in the first iteration ---*/ @@ -2827,37 +2828,33 @@ void CDriver::PrintDirectResidual(RECORDING kind_recording) { if (!addVals) RMSTable.AddColumn("rms_Rad" + iVar_iZone2string(0, iZone), fieldWidth); else RMSTable << log10(solvers[RAD_SOL]->GetRes_RMS(0)); } - - } - else if (configs->GetStructuralProblem()) { - + } else if (configs->GetStructuralProblem()) { if (configs->GetGeometricConditions() == STRUCT_DEFORMATION::LARGE){ if (!addVals) { RMSTable.AddColumn("UTOL-A", fieldWidth); RMSTable.AddColumn("RTOL-A", fieldWidth); RMSTable.AddColumn("ETOL-A", fieldWidth); - } - else { + } else { RMSTable << log10(solvers[FEA_SOL]->GetRes_FEM(0)) << log10(solvers[FEA_SOL]->GetRes_FEM(1)) << log10(solvers[FEA_SOL]->GetRes_FEM(2)); } - } - else{ + } else { if (!addVals) { RMSTable.AddColumn("log10[RMS Ux]", fieldWidth); RMSTable.AddColumn("log10[RMS Uy]", fieldWidth); if (nDim == 3) RMSTable.AddColumn("log10[RMS Uz]", fieldWidth); - } - else { + } else { RMSTable << log10(solvers[FEA_SOL]->GetRes_FEM(0)) << log10(solvers[FEA_SOL]->GetRes_FEM(1)); if (nDim == 3) RMSTable << log10(solvers[FEA_SOL]->GetRes_FEM(2)); } } - - } - else if (configs->GetHeatProblem()) { + if (configs->GetWeakly_Coupled_Heat()){ + if (!addVals) RMSTable.AddColumn("rms_Heat", fieldWidth); + else RMSTable << log10(solvers[HEAT_SOL]->GetRes_RMS(0)); + } + } else if (configs->GetHeatProblem()) { if (!addVals) RMSTable.AddColumn("rms_Heat" + iVar_iZone2string(0, iZone), fieldWidth); else RMSTable << log10(solvers[HEAT_SOL]->GetRes_RMS(0)); diff --git a/SU2_CFD/src/iteration/CDiscAdjFEAIteration.cpp b/SU2_CFD/src/iteration/CDiscAdjFEAIteration.cpp index a6c3bfaa014..a1425be3b7f 100644 --- a/SU2_CFD/src/iteration/CDiscAdjFEAIteration.cpp +++ b/SU2_CFD/src/iteration/CDiscAdjFEAIteration.cpp @@ -26,13 +26,20 @@ */ #include "../../include/iteration/CDiscAdjFEAIteration.hpp" +#include "../../include/iteration/CDiscAdjHeatIteration.hpp" #include "../../include/iteration/CFEAIteration.hpp" #include "../../include/solvers/CFEASolver.hpp" #include "../../include/output/COutput.hpp" -CDiscAdjFEAIteration::CDiscAdjFEAIteration(const CConfig *config) : CIteration(config), CurrentRecording(NONE) {} +CDiscAdjFEAIteration::CDiscAdjFEAIteration(const CConfig *config) : CIteration(config), CurrentRecording(NONE) { + if (config->GetWeakly_Coupled_Heat()) { + DiscAdjHeatIteration = new CDiscAdjHeatIteration(config); + } +} -CDiscAdjFEAIteration::~CDiscAdjFEAIteration() {} +CDiscAdjFEAIteration::~CDiscAdjFEAIteration() { + delete DiscAdjHeatIteration; +} void CDiscAdjFEAIteration::Preprocess(COutput* output, CIntegration**** integration, CGeometry**** geometry, CSolver***** solver, CNumerics****** numerics, CConfig** config, @@ -90,6 +97,10 @@ void CDiscAdjFEAIteration::Preprocess(COutput* output, CIntegration**** integrat solvers0[ADJFEA_SOL]->Preprocessing(geometry0, solvers0, config[iZone], MESH_0, 0, RUNTIME_ADJFEA_SYS, false); + if (DiscAdjHeatIteration) { + DiscAdjHeatIteration->Preprocess(output, integration, geometry, solver, numerics, config, surface_movement, + grid_movement, FFDBox, iZone, iInst); + } } void CDiscAdjFEAIteration::LoadDynamic_Solution(CGeometry**** geometry, CSolver***** solver, CConfig** config, @@ -115,23 +126,23 @@ void CDiscAdjFEAIteration::IterateDiscAdj(CGeometry**** geometry, CSolver***** s unsigned short iZone, unsigned short iInst, bool CrossTerm) { /*--- Extract the adjoints of the conservative input variables and store them for the next iteration ---*/ - - solver[iZone][iInst][MESH_0][ADJFEA_SOL]->ExtractAdjoint_Solution(geometry[iZone][iInst][MESH_0], config[iZone], - CrossTerm); - - solver[iZone][iInst][MESH_0][ADJFEA_SOL]->ExtractAdjoint_Variables(geometry[iZone][iInst][MESH_0], config[iZone]); + for (const auto iSol : {ADJFEA_SOL, ADJHEAT_SOL}) { + if (auto* sol = solver[iZone][iInst][MESH_0][iSol]; sol != nullptr) { + sol->ExtractAdjoint_Solution(geometry[iZone][iInst][MESH_0], config[iZone], CrossTerm); + sol->ExtractAdjoint_Variables(geometry[iZone][iInst][MESH_0], config[iZone]); + } + } } void CDiscAdjFEAIteration::RegisterInput(CSolver***** solver, CGeometry**** geometry, CConfig** config, unsigned short iZone, unsigned short iInst, RECORDING kind_recording) { if (kind_recording != RECORDING::MESH_COORDS) { - /*--- Register structural displacements as input ---*/ - - solver[iZone][iInst][MESH_0][ADJFEA_SOL]->RegisterSolution(geometry[iZone][iInst][MESH_0], config[iZone]); - - /*--- Register variables as input ---*/ - - solver[iZone][iInst][MESH_0][ADJFEA_SOL]->RegisterVariables(geometry[iZone][iInst][MESH_0], config[iZone]); + for (const auto iSol : {ADJFEA_SOL, ADJHEAT_SOL}) { + if (auto* sol = solver[iZone][iInst][MESH_0][iSol]; sol != nullptr) { + sol->RegisterSolution(geometry[iZone][iInst][MESH_0], config[iZone]); + sol->RegisterVariables(geometry[iZone][iInst][MESH_0], config[iZone]); + } + } } else { /*--- Register topology optimization densities (note direct solver) ---*/ @@ -146,6 +157,10 @@ void CDiscAdjFEAIteration::RegisterInput(CSolver***** solver, CGeometry**** geom void CDiscAdjFEAIteration::SetDependencies(CSolver***** solver, CGeometry**** geometry, CNumerics****** numerics, CConfig** config, unsigned short iZone, unsigned short iInst, RECORDING kind_recording) { + if (DiscAdjHeatIteration) { + DiscAdjHeatIteration->SetDependencies(solver, geometry, numerics, config, iZone, iInst, kind_recording); + } + auto dir_solver = solver[iZone][iInst][MESH_0][FEA_SOL]; auto adj_solver = solver[iZone][iInst][MESH_0][ADJFEA_SOL]; auto structural_geometry = geometry[iZone][iInst][MESH_0]; @@ -263,18 +278,25 @@ void CDiscAdjFEAIteration::SetDependencies(CSolver***** solver, CGeometry**** ge void CDiscAdjFEAIteration::RegisterOutput(CSolver***** solver, CGeometry**** geometry, CConfig** config, unsigned short iZone, unsigned short iInst) { - /*--- Register conservative variables as output of the iteration ---*/ - - solver[iZone][iInst][MESH_0][ADJFEA_SOL]->RegisterOutput(geometry[iZone][iInst][MESH_0], config[iZone]); + /*--- Register solution variables as output of the iteration. ---*/ + for (const auto iSol : {ADJFEA_SOL, ADJHEAT_SOL}) { + if (auto* sol = solver[iZone][iInst][MESH_0][iSol]; sol != nullptr) { + sol->RegisterOutput(geometry[iZone][iInst][MESH_0], config[iZone]); + } + } } void CDiscAdjFEAIteration::InitializeAdjoint(CSolver***** solver, CGeometry**** geometry, CConfig** config, unsigned short iZone, unsigned short iInst) { - /*--- Initialize the adjoints the conservative variables ---*/ + /*--- Initialize the adjoints of the solution variables. ---*/ AD::ResizeAdjoints(); AD::BeginUseAdjoints(); - solver[iZone][iInst][MESH_0][ADJFEA_SOL]->SetAdjoint_Output(geometry[iZone][iInst][MESH_0], config[iZone]); + for (const auto iSol : {ADJFEA_SOL, ADJHEAT_SOL}) { + if (auto* sol = solver[iZone][iInst][MESH_0][iSol]; sol != nullptr) { + sol->SetAdjoint_Output(geometry[iZone][iInst][MESH_0], config[iZone]); + } + } AD::EndUseAdjoints(); } @@ -285,8 +307,8 @@ bool CDiscAdjFEAIteration::Monitor(COutput* output, CIntegration**** integration /*--- Write the convergence history (only screen output) ---*/ output->SetHistoryOutput(geometry[iZone][INST_0][MESH_0], solver[iZone][INST_0][MESH_0], config[iZone], - config[iZone]->GetTimeIter(), config[iZone]->GetOuterIter(), - config[iZone]->GetInnerIter()); + config[iZone]->GetTimeIter(), config[iZone]->GetOuterIter(), + config[iZone]->GetInnerIter()); return output->GetConvergence(); } diff --git a/SU2_CFD/src/output/CAdjElasticityOutput.cpp b/SU2_CFD/src/output/CAdjElasticityOutput.cpp index 2320cabf5a5..6fc59bf7384 100644 --- a/SU2_CFD/src/output/CAdjElasticityOutput.cpp +++ b/SU2_CFD/src/output/CAdjElasticityOutput.cpp @@ -27,6 +27,7 @@ #include "../../include/output/CAdjElasticityOutput.hpp" +#include "../../include/output/CAdjHeatOutput.hpp" #include #include "../../../Common/include/geometry/CGeometry.hpp" @@ -34,6 +35,8 @@ CAdjElasticityOutput::CAdjElasticityOutput(CConfig *config, unsigned short nDim) : COutput(config, nDim, false) { + coupled_heat = config->GetWeakly_Coupled_Heat(); + /*--- Initialize number of variables ---*/ nVar_FEM = nDim; @@ -51,8 +54,10 @@ CAdjElasticityOutput::CAdjElasticityOutput(CConfig *config, unsigned short nDim) requestedScreenFields.emplace_back("INNER_ITER"); requestedScreenFields.emplace_back("ADJOINT_DISP_X"); requestedScreenFields.emplace_back("ADJOINT_DISP_Y"); + if (coupled_heat) requestedScreenFields.emplace_back("RMS_ADJ_TEMPERATURE"); requestedScreenFields.emplace_back("SENS_E_0"); requestedScreenFields.emplace_back("SENS_NU_0"); + if (coupled_heat) requestedScreenFields.emplace_back("SENS_GEO"); nRequestedScreenFields = requestedScreenFields.size(); } @@ -134,6 +139,12 @@ void CAdjElasticityOutput::SetHistoryOutputFields(CConfig *config){ AddHistoryOutput("BGS_ADJ_DISP_Z", "bgs[A_Uz]", ScreenOutputFormat::FIXED, "BGS_RES", "BGS residual of the adjoint Z displacement.", HistoryFieldType::RESIDUAL); } } + + if (coupled_heat) { + CAdjHeatOutput::SetHistoryOutputFieldsImpl(config, this); + AddHistoryOutput("LINSOL_ITER_HEAT", "LinSolIterHeat", ScreenOutputFormat::INTEGER, "LINSOL", "Number of iterations of the linear solver."); + AddHistoryOutput("LINSOL_RESIDUAL_HEAT", "LinSolResHeat", ScreenOutputFormat::FIXED, "LINSOL", "Residual of the linear solver."); + } } inline void CAdjElasticityOutput::LoadHistoryData(CConfig *config, CGeometry *geometry, CSolver **solver) { @@ -173,6 +184,13 @@ inline void CAdjElasticityOutput::LoadHistoryData(CConfig *config, CGeometry *ge } } + /*--- Add heat solver data if available. ---*/ + if (coupled_heat) { + CAdjHeatOutput::LoadHistoryDataImpl(config, geometry, solver, this); + SetHistoryOutputValue("LINSOL_ITER_HEAT", solver[ADJHEAT_SOL]->GetIterLinSolver()); + SetHistoryOutputValue("LINSOL_RESIDUAL_HEAT", log10(solver[ADJHEAT_SOL]->GetResLinSolver())); + } + ComputeSimpleCustomOutputs(config); } @@ -191,6 +209,11 @@ void CAdjElasticityOutput::LoadVolumeData(CConfig *config, CGeometry *geometry, if (nVar_FEM == 3) SetVolumeOutputValue("ADJOINT-Z", iPoint, Node_Struc->GetSolution(iPoint, 2)); + if (coupled_heat) { + const auto* Node_Heat = solver[ADJHEAT_SOL]->GetNodes(); + SetVolumeOutputValue("ADJ_TEMPERATURE", iPoint, Node_Heat->GetSolution(iPoint, 0)); + } + SetVolumeOutputValue("SENSITIVITY-X", iPoint, Node_Struc->GetSensitivity(iPoint, 0)); SetVolumeOutputValue("SENSITIVITY-Y", iPoint, Node_Struc->GetSensitivity(iPoint, 1)); if (nDim == 3) @@ -212,6 +235,14 @@ void CAdjElasticityOutput::LoadVolumeData(CConfig *config, CGeometry *geometry, SetVolumeOutputValue("SENS_ACCEL-Y", iPoint, Node_Struc->GetSolution_time_n(iPoint, 2 * nDim + 1)); if (nDim == 3) SetVolumeOutputValue("SENS_ACCEL-Z", iPoint, Node_Struc->GetSolution_time_n(iPoint, 8)); + + if (coupled_heat) { + const auto* Node_Heat = solver[ADJHEAT_SOL]->GetNodes(); + SetVolumeOutputValue("SENS_TEMP_N", iPoint, Node_Heat->GetSolution_time_n(iPoint, 0)); + if (config->GetTime_Marching() == TIME_MARCHING::DT_STEPPING_2ND) { + SetVolumeOutputValue("SENS_TEMP_N1", iPoint, Node_Heat->GetSolution_time_n1(iPoint, 0)); + } + } } void CAdjElasticityOutput::SetVolumeOutputFields(CConfig *config){ @@ -232,6 +263,10 @@ void CAdjElasticityOutput::SetVolumeOutputFields(CConfig *config){ AddVolumeOutput("ADJOINT-Z", "Adjoint_z", "SOLUTION", "adjoint of displacement in the z direction"); /// END_GROUP + if (coupled_heat) { + AddVolumeOutput("ADJ_TEMPERATURE", "Adjoint_Temperature", "SOLUTION", "Adjoint Temperature"); + } + /// BEGIN_GROUP: SENSITIVITY, DESCRIPTION: Geometrical sensitivities of the current objective function. /// DESCRIPTION: Sensitivity x-component. AddVolumeOutput("SENSITIVITY-X", "Sensitivity_x", "SENSITIVITY", "geometric sensitivity in the x direction"); @@ -261,4 +296,10 @@ void CAdjElasticityOutput::SetVolumeOutputFields(CConfig *config){ if (nDim == 3) AddVolumeOutput("SENS_ACCEL-Z", "SensitivityAccelN_z", "SENSITIVITY_N", "sensitivity to the previous z acceleration"); + if (coupled_heat) { + AddVolumeOutput("SENS_TEMP_N", "SensitivityTempN", "SENSITIVITY_N", "sensitivity to the previous temperature"); + if (config->GetTime_Marching() == TIME_MARCHING::DT_STEPPING_2ND) { + AddVolumeOutput("SENS_TEMP_N1", "SensitivityTempN1", "SENSITIVITY_N", "sensitivity to the previous-1 temperature"); + } + } } diff --git a/SU2_CFD/src/output/CAdjHeatOutput.cpp b/SU2_CFD/src/output/CAdjHeatOutput.cpp index 35875abe825..66cf5a336f7 100644 --- a/SU2_CFD/src/output/CAdjHeatOutput.cpp +++ b/SU2_CFD/src/output/CAdjHeatOutput.cpp @@ -89,25 +89,30 @@ CAdjHeatOutput::CAdjHeatOutput(CConfig *config, unsigned short nDim) : COutput(c CAdjHeatOutput::~CAdjHeatOutput() = default; -void CAdjHeatOutput::SetHistoryOutputFields(CConfig *config){ +void CAdjHeatOutput::SetHistoryOutputFieldsImpl(CConfig *config, COutput* output) { /// BEGIN_GROUP: RMS_RES, DESCRIPTION: The root-mean-square residuals of the conservative variables. /// DESCRIPTION: Root-mean square residual of the adjoint temperature. - AddHistoryOutput("RMS_ADJ_TEMPERATURE", "rms[A_T]", ScreenOutputFormat::FIXED, "RMS_RES", "Root-mean square residual of the adjoint temperature.", HistoryFieldType::RESIDUAL); + output->AddHistoryOutput("RMS_ADJ_TEMPERATURE", "rms[A_T]", ScreenOutputFormat::FIXED, "RMS_RES", "Root-mean square residual of the adjoint temperature.", HistoryFieldType::RESIDUAL); /// END_GROUP /// BEGIN_GROUP: MAX_RES, DESCRIPTION: The maximum residuals of the conservative variables. /// DESCRIPTION: Maximum residual of the adjoint temperature. - AddHistoryOutput("MAX_ADJ_TEMPERATURE", "max[A_T]", ScreenOutputFormat::FIXED, "MAX_RES", "Maximum residual of the adjoint temperature.", HistoryFieldType::RESIDUAL); + output->AddHistoryOutput("MAX_ADJ_TEMPERATURE", "max[A_T]", ScreenOutputFormat::FIXED, "MAX_RES", "Maximum residual of the adjoint temperature.", HistoryFieldType::RESIDUAL); /// BEGIN_GROUP: MAX_RES, DESCRIPTION: The root-mean-square residuals of the conservative variables. /// DESCRIPTION: Root-mean-square residual of the adjoint temperature. - AddHistoryOutput("BGS_ADJ_TEMPERATURE", "bgs[A_T]", ScreenOutputFormat::FIXED, "BGS_RES", "BGS residual of the adjoint temperature.", HistoryFieldType::RESIDUAL); + output->AddHistoryOutput("BGS_ADJ_TEMPERATURE", "bgs[A_T]", ScreenOutputFormat::FIXED, "BGS_RES", "BGS residual of the adjoint temperature.", HistoryFieldType::RESIDUAL); /// BEGIN_GROUP: SENSITIVITY, DESCRIPTION: Sensitivities of different geometrical or boundary values. /// DESCRIPTION: Sum of the geometrical sensitivities on all markers set in MARKER_MONITORING. - AddHistoryOutput("SENS_GEO", "Sens_Geo", ScreenOutputFormat::SCIENTIFIC, "SENSITIVITY", "Sum of the geometrical sensitivities on all markers set in MARKER_MONITORING.", HistoryFieldType::COEFFICIENT); + output->AddHistoryOutput("SENS_GEO", "Sens_Geo", ScreenOutputFormat::SCIENTIFIC, "SENSITIVITY", "Sum of the geometrical sensitivities on all markers set in MARKER_MONITORING.", HistoryFieldType::COEFFICIENT); /// END_GROUP +} + +void CAdjHeatOutput::SetHistoryOutputFields(CConfig *config){ + + SetHistoryOutputFieldsImpl(config, this); AddHistoryOutput("LINSOL_ITER", "LinSolIter", ScreenOutputFormat::INTEGER, "LINSOL", "Number of iterations of the linear solver."); AddHistoryOutput("LINSOL_RESIDUAL", "LinSolRes", ScreenOutputFormat::FIXED, "LINSOL", "Residual of the linear solver."); @@ -119,19 +124,23 @@ void CAdjHeatOutput::SetHistoryOutputFields(CConfig *config){ } -void CAdjHeatOutput::LoadHistoryData(CConfig *config, CGeometry *geometry, CSolver **solver) { +void CAdjHeatOutput::LoadHistoryDataImpl(CConfig *config, CGeometry *geometry, CSolver **solver, COutput* output) { CSolver* adjheat_solver = solver[ADJHEAT_SOL]; - SetHistoryOutputValue("RMS_ADJ_TEMPERATURE", log10(adjheat_solver->GetRes_RMS(0))); + output->SetHistoryOutputValue("RMS_ADJ_TEMPERATURE", log10(adjheat_solver->GetRes_RMS(0))); + output->SetHistoryOutputValue("MAX_ADJ_TEMPERATURE", log10(adjheat_solver->GetRes_Max(0))); + if (config->GetMultizone_Problem()) { + output->SetHistoryOutputValue("BGS_ADJ_TEMPERATURE", log10(adjheat_solver->GetRes_BGS(0))); + } + output->SetHistoryOutputValue("SENS_GEO", adjheat_solver->GetTotal_Sens_Geo()); +} - SetHistoryOutputValue("MAX_ADJ_TEMPERATURE", log10(adjheat_solver->GetRes_Max(0))); +void CAdjHeatOutput::LoadHistoryData(CConfig *config, CGeometry *geometry, CSolver **solver) { - if (multiZone) { - SetHistoryOutputValue("BGS_ADJ_TEMPERATURE", log10(adjheat_solver->GetRes_BGS(0))); - } + CSolver* adjheat_solver = solver[ADJHEAT_SOL]; - SetHistoryOutputValue("SENS_GEO", adjheat_solver->GetTotal_Sens_Geo()); + LoadHistoryDataImpl(config, geometry, solver, this); SetHistoryOutputValue("LINSOL_ITER", adjheat_solver->GetIterLinSolver()); SetHistoryOutputValue("LINSOL_RESIDUAL", log10(adjheat_solver->GetResLinSolver())); diff --git a/SU2_CFD/src/output/CElasticityOutput.cpp b/SU2_CFD/src/output/CElasticityOutput.cpp index cebd509656b..5bf9bfb4d6e 100644 --- a/SU2_CFD/src/output/CElasticityOutput.cpp +++ b/SU2_CFD/src/output/CElasticityOutput.cpp @@ -37,10 +37,6 @@ CElasticityOutput::CElasticityOutput(CConfig *config, unsigned short nDim) : COu coupled_heat = config->GetWeakly_Coupled_Heat(); dynamic = config->GetTime_Domain(); - /*--- Initialize number of variables ---*/ - if (linear_analysis) nVar_FEM = nDim; - if (nonlinear_analysis) nVar_FEM = 3; - /*--- Default fields for screen output ---*/ if (nRequestedHistoryFields == 0){ RequestCommonHistory(dynamic); diff --git a/SU2_CFD/src/solvers/CFEASolver.cpp b/SU2_CFD/src/solvers/CFEASolver.cpp index bce77518731..93184cd1407 100644 --- a/SU2_CFD/src/solvers/CFEASolver.cpp +++ b/SU2_CFD/src/solvers/CFEASolver.cpp @@ -3099,7 +3099,7 @@ void CFEASolver::Stiffness_Penalty(CGeometry *geometry, CNumerics **numerics, CC void CFEASolver::LoadRestart(CGeometry **geometry, CSolver ***solver, CConfig *config, int val_iter, bool val_update_geo) { - const bool dynamic = (config->GetTime_Domain()); + const bool dynamic = config->GetTime_Domain(); const bool fluid_structure = config->GetFSI_Simulation(); const bool discrete_adjoint = config->GetDiscrete_Adjoint(); diff --git a/SU2_CFD/src/solvers/CHeatSolver.cpp b/SU2_CFD/src/solvers/CHeatSolver.cpp index 24b2eba6d6a..f0eaf356271 100644 --- a/SU2_CFD/src/solvers/CHeatSolver.cpp +++ b/SU2_CFD/src/solvers/CHeatSolver.cpp @@ -201,6 +201,10 @@ void CHeatSolver::LoadRestart(CGeometry **geometry, CSolver ***solver, CConfig * // P, vx, vy (,vz) skipVars += 1 + nDim; } + if (config->GetStructuralProblem() && config->GetWeakly_Coupled_Heat()) { + skipVars += nDim; + if (config->GetTime_Domain()) skipVars += 2 * nDim; + } /*--- Read the restart data from either an ASCII or binary SU2 file. ---*/ diff --git a/SU2_CFD/src/solvers/CSolverFactory.cpp b/SU2_CFD/src/solvers/CSolverFactory.cpp index f6ceda06d84..31a31495367 100644 --- a/SU2_CFD/src/solvers/CSolverFactory.cpp +++ b/SU2_CFD/src/solvers/CSolverFactory.cpp @@ -178,8 +178,12 @@ CSolver** CSolverFactory::CreateSolverContainer(MAIN_SOLVER kindMainSolver, CCon } break; case MAIN_SOLVER::DISC_ADJ_FEM: - solver[FEA_SOL] = CreateSubSolver(SUB_SOLVER_TYPE::FEA, solver, geometry, config, iMGLevel); + solver[FEA_SOL] = CreateSubSolver(SUB_SOLVER_TYPE::FEA, solver, geometry, config, iMGLevel); solver[ADJFEA_SOL] = CreateSubSolver(SUB_SOLVER_TYPE::DISC_ADJ_FEA, solver, geometry, config, iMGLevel); + if (config->GetWeakly_Coupled_Heat()) { + solver[HEAT_SOL] = CreateSubSolver(SUB_SOLVER_TYPE::HEAT, solver, geometry, config, iMGLevel); + solver[ADJHEAT_SOL] = CreateSubSolver(SUB_SOLVER_TYPE::DISC_ADJ_HEAT, solver, geometry, config, iMGLevel); + } break; case MAIN_SOLVER::FEM_EULER: solver[FLOW_SOL] = CreateSubSolver(SUB_SOLVER_TYPE::DG_EULER, solver, geometry, config, iMGLevel); diff --git a/TestCases/fea_fsi/ThermalBeam_3d/configBeamNonlinear_3d.cfg b/TestCases/fea_fsi/ThermalBeam_3d/configBeamNonlinear_3d.cfg index cb04c9fba24..d8978b0a2b1 100644 --- a/TestCases/fea_fsi/ThermalBeam_3d/configBeamNonlinear_3d.cfg +++ b/TestCases/fea_fsi/ThermalBeam_3d/configBeamNonlinear_3d.cfg @@ -5,7 +5,6 @@ %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% SOLVER= ELASTICITY -MATH_PROBLEM= DIRECT GEOMETRIC_CONDITIONS= LARGE_DEFORMATIONS MATERIAL_MODEL= NEO_HOOKEAN MESH_FILENAME= ../StatBeam_3d/meshBeam_3d.su2 @@ -29,7 +28,6 @@ CONV_FILENAME= history_beam VOLUME_FILENAME= beam RESTART_FILENAME= restart_beam SOLUTION_FILENAME= restart_beam -OUTPUT_WRT_FREQ= 1 CONV_FIELD= REL_RMS_RTOL, REL_RMS_TEMPERATURE CONV_STARTITER= 0 diff --git a/TestCases/fea_fsi/ThermalBeam_3d/configBeamNonlinear_3d_ad.cfg b/TestCases/fea_fsi/ThermalBeam_3d/configBeamNonlinear_3d_ad.cfg new file mode 100644 index 00000000000..740729bd840 --- /dev/null +++ b/TestCases/fea_fsi/ThermalBeam_3d/configBeamNonlinear_3d_ad.cfg @@ -0,0 +1,55 @@ +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% SU2 configuration file % +% Case description: 3D beam with thermal expansion % +% File Version 8.4.0 "Harrier" % +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +SOLVER= ELASTICITY +GEOMETRIC_CONDITIONS= LARGE_DEFORMATIONS +MATERIAL_MODEL= NEO_HOOKEAN +MESH_FILENAME= ../StatBeam_3d/meshBeam_3d.su2 +ELASTICITY_MODULUS= 3E7 +POISSON_RATIO= 0.3 +MATERIAL_THERMAL_EXPANSION_COEFF= 2e-5 +MATERIAL_REFERENCE_TEMPERATURE= 288.15 +MATERIAL_DENSITY= 7854 +MARKER_CLAMPED= ( left, right ) +MARKER_PRESSURE= ( lower, 0, symleft, 0, symright, 0 ) +MARKER_LOAD= ( upper, 1, 1000, 0, -1, 0 ) +DISCADJ_LIN_SOLVER= FGCRODR +DISCADJ_LIN_PREC= ILU +LINEAR_SOLVER_ERROR= 1E-4 +LINEAR_SOLVER_ITER= 1000 +LINEAR_SOLVER_RESTART_FREQUENCY= 50 +LINEAR_SOLVER_RESTART_DEFLATION= 10 +MESH_FORMAT= SU2 +TABULAR_FORMAT= CSV +CONV_FILENAME= history_beam +VOLUME_FILENAME= beam +RESTART_FILENAME= restart_beam +SOLUTION_FILENAME= restart_beam + +CONV_FIELD= REL_ADJOINT_DISP_X, REL_RMS_ADJ_TEMPERATURE +CONV_STARTITER= 0 +CONV_RESIDUAL_MINVAL= -5 +INNER_ITER= 15 + +% Coupling with heat solver. +WEAKLY_COUPLED_HEAT_EQUATION= YES +FREESTREAM_TEMPERATURE= 300 +SPECIFIC_HEAT_CP= 460 +THERMAL_CONDUCTIVITY_CONSTANT= 45 + +% Copy markers to allow specifying boundary conditions for both solvers on the +% same markers. +MARKER_CREATE_COPY= ( left, left_heat, right, right_heat ) +MARKER_ISOTHERMAL= ( left_heat, 400, right_heat, 300 ) + +NUM_METHOD_GRAD= GREEN_GAUSS +TIME_DISCRE_HEAT= EULER_IMPLICIT +CFL_NUMBER= 1e8 + +MARKER_MONITORING= ( left_heat ) +SCREEN_OUTPUT= INNER_ITER, RMS_RES, LINSOL, SENSITIVITY + +OBJECTIVE_FUNCTION = TOTAL_HEATFLUX diff --git a/TestCases/fea_fsi/ThermalBeam_3d/configBeam_3d.cfg b/TestCases/fea_fsi/ThermalBeam_3d/configBeam_3d.cfg index 60afa90b837..c2f263d2108 100644 --- a/TestCases/fea_fsi/ThermalBeam_3d/configBeam_3d.cfg +++ b/TestCases/fea_fsi/ThermalBeam_3d/configBeam_3d.cfg @@ -5,7 +5,6 @@ %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% SOLVER= ELASTICITY -MATH_PROBLEM= DIRECT GEOMETRIC_CONDITIONS= SMALL_DEFORMATIONS MATERIAL_MODEL= LINEAR_ELASTIC MESH_FILENAME= ../StatBeam_3d/meshBeam_3d.su2 diff --git a/TestCases/parallel_regression.py b/TestCases/parallel_regression.py index f5bea2d55b8..6f367c68bdb 100755 --- a/TestCases/parallel_regression.py +++ b/TestCases/parallel_regression.py @@ -1269,7 +1269,6 @@ def main(): statbeam3d.test_iter = 0 statbeam3d.test_vals = [-6.058758, -5.750933, -5.892188, 110190] statbeam3d.test_vals_aarch64 = [-6.062693, -5.769132, -5.891190, 110190] - statbeam3d.command = TestCase.Command(exec = "parallel_computation_fsi.py", param = "-f") test_list.append(statbeam3d) # Static beam, 3d with coupled temperature @@ -1278,7 +1277,6 @@ def main(): thermal_beam_3d.cfg_file = "configBeam_3d.cfg" thermal_beam_3d.test_iter = 4 thermal_beam_3d.test_vals = [-8.137623, -7.840833, -7.903285, -13.978110, 217, -4.093241, 39, -4.072614, 1.3676e+05, 75.0] - thermal_beam_3d.command = TestCase.Command(exec = "parallel_computation_fsi.py", param = "-f") test_list.append(thermal_beam_3d) # Static beam, 3d with coupled temperature, nonlinear elasticity @@ -1287,7 +1285,6 @@ def main(): thermal_beam_nl_3d.cfg_file = "configBeamNonlinear_3d.cfg" thermal_beam_nl_3d.test_iter = 8 thermal_beam_nl_3d.test_vals = [-7.564309, -2.992893, -12.242503, -14.068322, 57, -4.017665, 24, -4.204804, 138710, 75.233] - thermal_beam_nl_3d.command = TestCase.Command(exec = "parallel_computation_fsi.py", param = "-f") test_list.append(thermal_beam_nl_3d) # Rotating cylinder, 3d diff --git a/TestCases/parallel_regression_AD.py b/TestCases/parallel_regression_AD.py index bf1547db62e..7b931c872a8 100644 --- a/TestCases/parallel_regression_AD.py +++ b/TestCases/parallel_regression_AD.py @@ -247,6 +247,15 @@ def main(): discadj_fea.test_vals_aarch64 = [-2.849646, -3.238577, -0.000364, -8.708700] #last 4 columns test_list.append(discadj_fea) + # Thermoelastic problem + # Derivative of heat flux wrt Poisson's ratio (due to increase in cross section) verified via finite differences. + discadj_thermoelastic = TestCase('discadj_thermoelastic') + discadj_thermoelastic.cfg_dir = "fea_fsi/ThermalBeam_3d" + discadj_thermoelastic.cfg_file = "configBeamNonlinear_3d_ad.cfg" + discadj_thermoelastic.test_iter = 10 + discadj_thermoelastic.test_vals = [-5.355510, -5.293378, -6.164317, -6.433862, 43, -4.049556, 27, -4.164192, 0, 0.192640, 0] + test_list.append(discadj_thermoelastic) + ################################### ### Disc. adj. heat ### ###################################