diff --git a/Common/include/ad_structure.hpp b/Common/include/ad_structure.hpp index 2826b3f41b0d..822d68b0022e 100644 --- a/Common/include/ad_structure.hpp +++ b/Common/include/ad_structure.hpp @@ -65,17 +65,42 @@ namespace AD{ bool TapeActive(); /*! - * \brief Registers the variable as an input. I.e. as a leaf of the computational graph. + * \brief Prints out tape statistics. + */ + void PrintStatistics(); + + /*! + * \brief Registers the variable as an input and saves internal data (indices). I.e. as a leaf of the computational graph. * \param[in] data - The variable to be registered as input. */ void RegisterInput(su2double &data); + /*! + * \brief Registers the variable as an input. I.e. as a leaf of the computational graph. + * \param[in] data - The variable to be registered as input. + */ + void RegisterInput_intIndexBased(su2double &data); + /*! * \brief Registers the variable as an output. I.e. as the root of the computational graph. * \param[in] data - The variable to be registered as output. */ void RegisterOutput(su2double &data); + /*! + * \brief Sets the adjoint value at index to val + * \param[in] index - Position in the adjoint vector. + * \param[in] val - adjoint value to be set. + */ + void SetDerivative(int index, const double val); + + /*! + * \brief Extracts the adjoint value at index + * \param[in] index - position in the adjoint vector where the derivative will be extracted. + * \return Derivative value. + */ + double GetDerivative(int index); + /*! * \brief Clears the currently stored adjoints but keeps the computational graph. */ @@ -84,7 +109,14 @@ namespace AD{ /*! * \brief Computes the adjoints, i.e. the derivatives of the output with respect to the input variables. */ - void ComputeAdjoint(); + void ComputeAdjoint(); + + /*! + * \brief Computes the adjoints, i.e. the derivatives of the output with respect to the input variables. + * \param[in] enter - Position where we start evaluating the tape. + * \param[in] leave - Position where we stop evaluating the tape. + */ + void ComputeAdjoint(unsigned short enter, unsigned short leave); /*! * \brief Reset the tape structure to be ready for a new recording. @@ -209,6 +241,17 @@ namespace AD{ */ void EndExtFunc(); + /*! + * \brief Evaluates and saves gradient data from a variable. + * \param[in] data - variable whose gradient information will be extracted. + * \param[in] index - where obtained gradient information will be stored. + */ + void SetAdjIndex(int &index, const su2double &data); + + /*! + * \brief Pushes back the current tape position to the tape position's vector. + */ + void Push_TapePosition(); } /*--- Macro to begin and end sections with a passive tape ---*/ diff --git a/Common/include/ad_structure.inl b/Common/include/ad_structure.inl index 648e896a301b..ad8184f4c3cf 100644 --- a/Common/include/ad_structure.inl +++ b/Common/include/ad_structure.inl @@ -68,6 +68,8 @@ namespace AD{ extern su2double::TapeType::Position StartPosition, EndPosition; + extern std::vector TapePositions; + extern std::vector localInputValues; extern std::vector localOutputValues; @@ -77,6 +79,8 @@ namespace AD{ inline void RegisterInput(su2double &data) {AD::globalTape.registerInput(data); inputValues.push_back(data.getGradientData());} + inline void RegisterInput_intIndexBased(su2double &data) {AD::globalTape.registerInput(data);} + inline void RegisterOutput(su2double& data) {AD::globalTape.registerOutput(data);} inline void ResetInput(su2double &data) {data.getGradientData() = su2double::GradientData();} @@ -85,19 +89,43 @@ namespace AD{ inline void StopRecording() {AD::globalTape.setPassive();} - inline bool TapeActive() {return AD::globalTape.isActive();} + inline bool TapeActive() { return AD::globalTape.isActive(); } + + inline void PrintStatistics() {AD::globalTape.printStatistics();} inline void ClearAdjoints() {AD::globalTape.clearAdjoints(); } inline void ComputeAdjoint() {AD::globalTape.evaluate(); adjointVectorPosition = 0;} + inline void ComputeAdjoint(unsigned short enter, unsigned short leave) { + AD::globalTape.evaluate(TapePositions[enter], TapePositions[leave]); + if (leave == 0) { + adjointVectorPosition = 0; + } + } + inline void Reset() { + globalTape.reset(); if (inputValues.size() != 0) { - globalTape.reset(); adjointVectorPosition = 0; inputValues.clear(); } + if (TapePositions.size() != 0) { + TapePositions.clear(); + } + } + + inline void SetAdjIndex(int &index, const su2double &data) { + index = data.getGradientData(); + } + + inline void SetDerivative(int index, const double val) { + AD::globalTape.setGradient(index, val); + } + + inline double GetDerivative(int index) { + return AD::globalTape.getGradient(index); } inline void SetPreaccIn(const su2double &data) { @@ -167,6 +195,10 @@ namespace AD{ } } + inline void Push_TapePosition() { + TapePositions.push_back(AD::globalTape.getPosition()); + } + inline void EndPreacc(){ if (PreaccActive) { PreaccHelper.finish(false); @@ -239,18 +271,30 @@ namespace AD{ inline void RegisterInput(su2double &data) {} + inline void RegisterInput_intIndexBased(su2double &data) {} + inline void RegisterOutput(su2double& data) {} inline void StartRecording() {} inline void StopRecording() {} - inline bool TapeActive() {return false;} + inline bool TapeActive() { return false; } + + inline void PrintStatistics() {} inline void ClearAdjoints() {} inline void ComputeAdjoint() {} + inline void ComputeAdjoint(unsigned short enter, unsigned short leave) {} + + inline void SetAdjIndex(int &index, const su2double &data) {} + + inline void SetDerivative(int index, const double val) {} + + inline double GetDerivative(int position) { return 0.0; } + inline void Reset() {} inline void ResetInput(su2double &data) {} @@ -270,7 +314,9 @@ namespace AD{ inline void StartPreacc() {} inline void EndPreacc() {} - + + inline void Push_TapePosition() {} + inline void StartExtFunc(bool storePrimalInput, bool storePrimalOutput){} inline void SetExtFuncIn(const su2double &data) {} diff --git a/Common/include/config_structure.hpp b/Common/include/config_structure.hpp index 95c7c13cd9d7..e221482d17a8 100644 --- a/Common/include/config_structure.hpp +++ b/Common/include/config_structure.hpp @@ -417,7 +417,6 @@ class CConfig { CFLRedCoeff_AdjFlow, /*!< \brief CFL reduction coefficient for the adjoint problem. */ CFLRedCoeff_AdjTurb, /*!< \brief CFL reduction coefficient for the adjoint problem. */ CFLFineGrid, /*!< \brief CFL of the finest grid. */ - CFLSolid, /*!< \brief CFL in (heat) solid solvers. */ Max_DeltaTime, /*!< \brief Max delta time. */ Unst_CFL; /*!< \brief Unsteady CFL number. */ bool ReorientElements; /*!< \brief Flag for enabling element reorientation. */ @@ -758,6 +757,7 @@ class CConfig { Wrt_SharpEdges, /*!< \brief Write residuals to solution file */ Wrt_Halo, /*!< \brief Write rind layers in solution files */ Wrt_Performance, /*!< \brief Write the performance summary at the end of a calculation. */ + Wrt_AD_Statistics, /*!< \brief Write the tape statistics (discrete adjoint). */ Wrt_MeshQuality, /*!< \brief Write the mesh quality statistics to the visualization files. */ Wrt_InletFile, /*!< \brief Write a template inlet profile file */ Wrt_Slice, /*!< \brief Write 1D slice of a 2D cartesian solution */ @@ -874,6 +874,7 @@ class CConfig { bool Sine_Load; /*!< \brief option for sine load */ su2double *SineLoad_Coeff; /*!< \brief Stores the load coefficient */ su2double Thermal_Diffusivity; /*!< \brief Thermal diffusivity used in the heat solver. */ + bool CHT_Robin; /*!< \brief Option for boundary condition method at CHT interfaces. */ su2double Cyclic_Pitch, /*!< \brief Cyclic pitch for rotorcraft simulations. */ Collective_Pitch; /*!< \brief Collective pitch for rotorcraft simulations. */ su2double Mach_Motion; /*!< \brief Mach number based on mesh velocity and freestream quantities. */ @@ -975,7 +976,8 @@ class CConfig { Nonphys_Reconstr; /*!< \brief Current number of non-physical reconstructions for 2nd-order upwinding. */ bool ParMETIS; /*!< \brief Boolean for activating ParMETIS mode (while testing). */ unsigned short DirectDiff; /*!< \brief Direct Differentation mode. */ - bool DiscreteAdjoint; /*!< \brief AD-based discrete adjoint mode. */ + bool DiscreteAdjoint, /*!< \brief AD-based discrete adjoint mode. */ + FullTape; /*!< \brief Full tape mode for coupled discrete adjoints. */ unsigned long Wrt_Surf_Freq_DualTime; /*!< \brief Writing surface solution frequency for Dual Time. */ su2double Const_DES; /*!< \brief Detached Eddy Simulation Constant. */ unsigned short Kind_HybridRANSLES; /*!< \brief Kind of Hybrid RANS/LES. */ @@ -1404,7 +1406,7 @@ class CConfig { /*! * \brief Constructor of the class which reads the input file. */ - CConfig(char case_filename[MAX_STRING_SIZE], unsigned short val_software, unsigned short val_nZone, bool verb_high); + CConfig(char case_filename[MAX_STRING_SIZE], unsigned short val_software, bool verb_high); /*! * \brief Constructor of the class which reads the input file and uses default options from another config. @@ -2757,13 +2759,6 @@ class CConfig { * \return CFL number for each grid. */ su2double GetCFL(unsigned short val_mesh); - - /*! - * \brief Get the Courant Friedrich Levi number for solid solvers. - * \param[in] val_mesh - Index of the mesh were the CFL is applied. - * \return CFL number for each grid. - */ - su2double GetCFL_Solid(void); /*! * \brief Get the Courant Friedrich Levi number for each grid. @@ -3252,6 +3247,12 @@ class CConfig { * \return TRUE means that the performance summary will be written at the end of a calculation. */ bool GetWrt_Performance(void); + + /*! + * \brief Get information about the computational graph (e.g. memory usage) when using AD in reverse mode. + * \return TRUE means that the tape statistics will be written after each recording. + */ + bool GetWrt_AD_Statistics(void); /*! * \brief Get information about writing the mesh quality metrics to the visualization files. @@ -8434,6 +8435,12 @@ class CConfig { * \return the discrete adjoint indicator. */ bool GetDiscrete_Adjoint(void); + + /*! + * \brief Get the indicator whether we want to use full (coupled) tapes. + * \return the full tape indicator. + */ + bool GetFull_Tape(void); /*! * \brief Get the indicator whether we want to benchmark the MPI performance of FSI problems @@ -8940,6 +8947,12 @@ class CConfig { */ bool GetWeakly_Coupled_Heat(void); + /*! + * \brief Get the boundary condition method for CHT. + * \return YES if Robin BC is used. + */ + bool GetCHT_Robin(void); + /*! * \brief Check if values passed to the BC_HeatFlux-Routine are already integrated. * \return YES if the passed values is the integrated heat flux over the marker's surface. diff --git a/Common/include/config_structure.inl b/Common/include/config_structure.inl index 1060f486668a..f2f601536d3e 100644 --- a/Common/include/config_structure.inl +++ b/Common/include/config_structure.inl @@ -753,8 +753,6 @@ inline unsigned short CConfig::GetGeometryMode(void) { return GeometryMode; } inline su2double CConfig::GetCFL(unsigned short val_mesh) { return CFL[val_mesh]; } -inline su2double CConfig::GetCFL_Solid(void) { return CFLSolid; } - inline void CConfig::SetCFL(unsigned short val_mesh, su2double val_cfl) { CFL[val_mesh] = val_cfl; } inline su2double CConfig::GetUnst_CFL(void) { return Unst_CFL; } @@ -1669,6 +1667,8 @@ inline bool CConfig::GetWrt_Halo(void) { return Wrt_Halo; } inline bool CConfig::GetWrt_Performance(void) { return Wrt_Performance; } +inline bool CConfig::GetWrt_AD_Statistics(void) { return Wrt_AD_Statistics; } + inline bool CConfig::GetWrt_MeshQuality(void) { return Wrt_MeshQuality; } inline bool CConfig::GetWrt_InletFile(void) { return Wrt_InletFile; } @@ -1905,6 +1905,8 @@ inline unsigned short CConfig::GetDirectDiff() { return DirectDiff;} inline bool CConfig::GetDiscrete_Adjoint() { return DiscreteAdjoint;} +inline bool CConfig::GetFull_Tape() { return FullTape; } + inline unsigned short CConfig::GetRiemann_Solver_FEM(void) {return Riemann_Solver_FEM;} inline su2double CConfig::GetQuadrature_Factor_Straight(void) {return Quadrature_Factor_Straight;} @@ -1925,6 +1927,8 @@ inline bool CConfig::GetJacobian_Spatial_Discretization_Only(void) {return Jacob inline bool CConfig::GetWeakly_Coupled_Heat(void) { return Weakly_Coupled_Heat; } +inline bool CConfig::GetCHT_Robin(void) { return CHT_Robin; } + inline bool CConfig::GetIntegrated_HeatFlux(void) { return Integrated_HeatFlux; } inline bool CConfig::GetAD_Mode(void) { return AD_Mode;} diff --git a/Common/include/dual_grid_structure.hpp b/Common/include/dual_grid_structure.hpp index d749a0a55702..94bd81a44f37 100644 --- a/Common/include/dual_grid_structure.hpp +++ b/Common/include/dual_grid_structure.hpp @@ -176,6 +176,8 @@ class CPoint : public CDualGrid { unsigned short nNeighbor; /*!< \brief Number of neighbors. */ bool Flip_Orientation; /*!< \brief Flip the orientation of the normal. */ su2double MaxLength; /*!< \brief The maximum cell-center to cell-center length. */ + int *Input_AdjIndices, /*!< \brief Indices of Coord variables in the adjoint vector. */ + *Output_AdjIndices; /*!< \brief Indices of Coord variables in the adjoint vector after having been updated. */ public: @@ -289,7 +291,25 @@ class CPoint : public CDualGrid { * \param[in] val_dim - Position to store the coordinate. * \param[in] val_coord - Coordinate for val_dim. */ - void SetCoord(unsigned short val_dim, su2double val_coord); + void SetCoord(unsigned short val_dim, su2double val_coord); + + /*! + * \brief Set the adjoint vector indices of Coord vector. + * \param[in] input - Save them to the input or output indices vector. + */ + void SetAdjIndices(bool input); + + /*! + * \brief Set the adjoint values of the (geometric) coordinates. + * \param[in] adj_sol - Adjoint values of the Coord variables. + */ + void SetAdjointSolution(const su2double *adj_sol); + + /*! + * \brief Get the adjoint values of the (geometric) coordinates. + * \param[in] adj_sol - Adjoint values of the Coord variables. + */ + su2double GetAdjointSolution(unsigned short iDim); /*! * \brief Get the coordinates of the control volume. diff --git a/Common/src/ad_structure.cpp b/Common/src/ad_structure.cpp index 8ed974bd02ef..7511fa41d34c 100644 --- a/Common/src/ad_structure.cpp +++ b/Common/src/ad_structure.cpp @@ -49,6 +49,7 @@ namespace AD { su2double::TapeType& globalTape = su2double::getGlobalTape(); su2double::TapeType::Position StartPosition, EndPosition; + std::vector TapePositions; bool Status = false; bool PreaccActive = false; diff --git a/Common/src/config_structure.cpp b/Common/src/config_structure.cpp index 97f2d017d172..d6b52b87bc32 100644 --- a/Common/src/config_structure.cpp +++ b/Common/src/config_structure.cpp @@ -63,7 +63,7 @@ vector GEMM_Profile_MaxTime; /*!< \brief Maximum time spent for thi #include "../include/ad_structure.hpp" #include "../include/toolboxes/printing_toolbox.hpp" -CConfig::CConfig(char case_filename[MAX_STRING_SIZE], unsigned short val_software, unsigned short val_nZone, bool verb_high) { +CConfig::CConfig(char case_filename[MAX_STRING_SIZE], unsigned short val_software, bool verb_high) { base_config = true; @@ -72,8 +72,8 @@ CConfig::CConfig(char case_filename[MAX_STRING_SIZE], unsigned short val_softwar rank = SU2_MPI::GetRank(); size = SU2_MPI::GetSize(); - iZone = val_nZone; - nZone = val_nZone; + iZone = 0; + nZone = 1; /*--- Initialize pointers to Null---*/ @@ -90,6 +90,10 @@ CConfig::CConfig(char case_filename[MAX_STRING_SIZE], unsigned short val_softwar /*--- Set the default values for all of the options that weren't set ---*/ SetDefault(); + + /*--- Set number of zone ---*/ + + SetnZone(); /*--- Configuration file postprocessing ---*/ @@ -157,6 +161,8 @@ CConfig::CConfig(CConfig* config, char case_filename[MAX_STRING_SIZE], unsigned if ((rank == MASTER_NODE) && verb_high) SetOutput(val_software, val_iZone); + Multizone_Problem = config->GetMultizone_Problem(); + } CConfig::CConfig(char case_filename[MAX_STRING_SIZE], unsigned short val_software) { @@ -874,6 +880,8 @@ void CConfig::SetConfig_Options() { addEnumOption("MULTIZONE_SOLVER", Kind_MZSolver, Multizone_Map, MZ_BLOCK_GAUSS_SEIDEL); /*!\brief MATH_PROBLEM \n DESCRIPTION: Mathematical problem \n Options: DIRECT, ADJOINT \ingroup Config*/ addMathProblemOption("MATH_PROBLEM", ContinuousAdjoint, false, DiscreteAdjoint, false, Restart_Flow, false); + /*!\brief FULL_TAPE \n DESCRIPTION: Use full (coupled) tapes for multiphysics discrete adjoint. \ingroup Config*/ + addBoolOption("FULL_TAPE", FullTape, YES); /*!\brief KIND_TURB_MODEL \n DESCRIPTION: Specify turbulence model \n Options: see \link Turb_Model_Map \endlink \n DEFAULT: NO_TURB_MODEL \ingroup Config*/ addEnumOption("KIND_TURB_MODEL", Kind_Turb_Model, Turb_Model_Map, NO_TURB_MODEL); /*!\brief KIND_TRANS_MODEL \n DESCRIPTION: Specify transition model OPTIONS: see \link Trans_Model_Map \endlink \n DEFAULT: NO_TRANS_MODEL \ingroup Config*/ @@ -1369,8 +1377,6 @@ void CConfig::SetConfig_Options() { addEnumOption("TIME_MARCHING", TimeMarching, TimeMarching_Map, STEADY); /* DESCRIPTION: Courant-Friedrichs-Lewy condition of the finest grid */ addDoubleOption("CFL_NUMBER", CFLFineGrid, 1.25); - /* DESCRIPTION: Courant-Friedrichs-Lewy condition of the finest grid in (heat fvm) solid solvers */ - addDoubleOption("CFL_NUMBER_SOLID", CFLSolid, 1.25); /* DESCRIPTION: Max time step in local time stepping simulations */ addDoubleOption("MAX_DELTA_TIME", Max_DeltaTime, 1000000); /* DESCRIPTION: Activate The adaptive CFL number. */ @@ -1791,6 +1797,8 @@ void CConfig::SetConfig_Options() { addBoolOption("WRT_HALO", Wrt_Halo, false); /* DESCRIPTION: Output the performance summary to the console at the end of SU2_CFD \ingroup Config*/ addBoolOption("WRT_PERFORMANCE", Wrt_Performance, false); + /* DESCRIPTION: Output the tape statistics (discrete adjoint) \ingroup Config*/ + addBoolOption("WRT_AD_STATISTICS", Wrt_AD_Statistics, false); /* DESCRIPTION: Write the mesh quality metrics to the visualization files. \ingroup Config*/ addBoolOption("WRT_MESH_QUALITY", Wrt_MeshQuality, false); /* DESCRIPTION: Output a 1D slice of a 2D cartesian solution \ingroup Config*/ @@ -2270,6 +2278,10 @@ void CConfig::SetConfig_Options() { /*!\par CONFIG_CATEGORY: Heat solver \ingroup Config*/ /*--- options related to the heat solver ---*/ + /* DESCRIPTION: Use Robin (default) or Neumann BC at CHT interface. */ + /* Options: NO, YES \ingroup Config */ + addBoolOption("CHT_ROBIN", CHT_Robin, true); + /* DESCRIPTION: Thermal diffusivity constant */ addDoubleOption("THERMAL_DIFFUSIVITY", Thermal_Diffusivity, 1.172E-5); @@ -4553,6 +4565,9 @@ void CConfig::SetPostprocessing(unsigned short val_software, unsigned short val_ case FEM_ELASTICITY: Kind_Solver = DISC_ADJ_FEM; break; + case HEAT_EQUATION_FVM: + Kind_Solver = DISC_ADJ_HEAT; + break; default: break; } @@ -6237,7 +6252,6 @@ void CConfig::SetOutput(unsigned short val_software, unsigned short val_izone) { cout << endl <<"------------------ Convergence Criteria ( Zone " << iZone << " ) ---------------------" << endl; - cout << "Maximum number of solver subiterations: " << nInnerIter <<"."<< endl; cout << "Maximum number of physical time-steps: " << nTimeIter <<"."<< endl; @@ -7486,7 +7500,7 @@ string CConfig::GetFilename(string filename, string ext, unsigned long Iter){ filename = filename + string(ext); /*--- Append the zone number if multizone problems ---*/ - if (GetnZone() > 1) + if (Multizone_Problem) filename = GetMultizone_FileName(filename, GetiZone(), ext); /*--- Append the zone number if multiple instance problems ---*/ @@ -7539,7 +7553,7 @@ string CConfig::GetMultizone_FileName(string val_filename, int val_iZone, string unsigned short lastindex = multizone_filename.find_last_of("."); multizone_filename = multizone_filename.substr(0, lastindex); - if (GetnZone() > 1 ) { + if (Multizone_Problem) { SPRINTF (buffer, "_%d", SU2_TYPE::Int(val_iZone)); multizone_filename.append(string(buffer)); } @@ -7554,7 +7568,7 @@ string CConfig::GetMultizone_HistoryFileName(string val_filename, int val_iZone, char buffer[50]; unsigned short lastindex = multizone_filename.find_last_of("."); multizone_filename = multizone_filename.substr(0, lastindex); - if (GetnZone() > 1 ) { + if (Multizone_Problem) { SPRINTF (buffer, "_%d", SU2_TYPE::Int(val_iZone)); multizone_filename.append(string(buffer)); } diff --git a/Common/src/dual_grid_structure.cpp b/Common/src/dual_grid_structure.cpp index 233c8c06dc7e..4ed88ce88e0a 100644 --- a/Common/src/dual_grid_structure.cpp +++ b/Common/src/dual_grid_structure.cpp @@ -52,10 +52,11 @@ CPoint::CPoint(unsigned short val_nDim, unsigned long val_globalindex, CConfig * Point.clear(); nPoint = 0; Edge.clear(); - Volume = NULL; Vertex = NULL; - Coord = NULL; Coord_Old = NULL; Coord_Sum = NULL; - Coord_n = NULL; Coord_n1 = NULL; Coord_p1 = NULL; - GridVel = NULL; GridVel_Grad = NULL; + Volume = NULL; Vertex = NULL; + Coord = NULL; Coord_Old = NULL; Coord_Sum = NULL; + Coord_n = NULL; Coord_n1 = NULL; Coord_p1 = NULL; + GridVel = NULL; GridVel_Grad = NULL; + Input_AdjIndices = NULL; Output_AdjIndices = NULL; /*--- Volume (0 -> Vol_nP1, 1-> Vol_n, 2 -> Vol_nM1 ) and coordinates of the control volume ---*/ @@ -72,6 +73,11 @@ CPoint::CPoint(unsigned short val_nDim, unsigned long val_globalindex, CConfig * Coord = new su2double[nDim]; + if(config->GetAD_Mode() && config->GetMultizone_Problem()) { + Input_AdjIndices = new int[nDim]; + Output_AdjIndices = new int[nDim]; + } + /*--- Indicator if the control volume has been agglomerated ---*/ Parent_CV = 0; Agglomerate = false; @@ -161,10 +167,11 @@ CPoint::CPoint(su2double val_coord_0, su2double val_coord_1, unsigned long val_g Point.clear(); nPoint = 0; Edge.clear(); - Volume = NULL; Vertex = NULL; - Coord = NULL; Coord_Old = NULL; Coord_Sum = NULL; - Coord_n = NULL; Coord_n1 = NULL; Coord_p1 = NULL; - GridVel = NULL; GridVel_Grad = NULL; + Volume = NULL; Vertex = NULL; + Coord = NULL; Coord_Old = NULL; Coord_Sum = NULL; + Coord_n = NULL; Coord_n1 = NULL; Coord_p1 = NULL; + GridVel = NULL; GridVel_Grad = NULL; + Input_AdjIndices = NULL; Output_AdjIndices = NULL; /*--- Volume (0 -> Vol_nP1, 1-> Vol_n, 2 -> Vol_nM1 ) and coordinates of the control volume ---*/ @@ -183,6 +190,11 @@ CPoint::CPoint(su2double val_coord_0, su2double val_coord_1, unsigned long val_g Coord[0] = val_coord_0; Coord[1] = val_coord_1; + if(config->GetAD_Mode() && config->GetMultizone_Problem()) { + Input_AdjIndices = new int[nDim]; + Output_AdjIndices = new int[nDim]; + } + /*--- Indicator if the control volume has been agglomerated ---*/ Parent_CV = 0; Agglomerate = false; @@ -271,10 +283,11 @@ CPoint::CPoint(su2double val_coord_0, su2double val_coord_1, su2double val_coord Point.clear(); nPoint = 0; Edge.clear(); - Volume = NULL; Vertex = NULL; - Coord = NULL; Coord_Old = NULL; Coord_Sum = NULL; - Coord_n = NULL; Coord_n1 = NULL; Coord_p1 = NULL; - GridVel = NULL; GridVel_Grad = NULL; + Volume = NULL; Vertex = NULL; + Coord = NULL; Coord_Old = NULL; Coord_Sum = NULL; + Coord_n = NULL; Coord_n1 = NULL; Coord_p1 = NULL; + GridVel = NULL; GridVel_Grad = NULL; + Input_AdjIndices = NULL; Output_AdjIndices = NULL; /*--- Volume (0 -> Vol_nP1, 1-> Vol_n, 2 -> Vol_nM1 ) and coordinates of the control volume ---*/ if ( config->GetTime_Marching() == NO ) { @@ -293,6 +306,11 @@ CPoint::CPoint(su2double val_coord_0, su2double val_coord_1, su2double val_coord Coord[1] = val_coord_1; Coord[2] = val_coord_2; + if(config->GetAD_Mode() && config->GetMultizone_Problem()) { + Input_AdjIndices = new int[nDim]; + Output_AdjIndices = new int[nDim]; + } + /*--- Indicator if the control volume has been agglomerated ---*/ Parent_CV = 0; Agglomerate = false; @@ -389,8 +407,9 @@ CPoint::~CPoint() { delete [] GridVel_Grad[iDim]; delete [] GridVel_Grad; } - -} + if (Input_AdjIndices != NULL) delete[] Input_AdjIndices; + if (Output_AdjIndices != NULL) delete[] Output_AdjIndices; + } void CPoint::SetPoint(unsigned long val_point) { @@ -430,6 +449,27 @@ void CPoint::SetBoundary(unsigned short val_nmarker) { } +void CPoint::SetAdjIndices(bool input) { + for (unsigned short iDim = 0; iDim < nDim; iDim++) { + if(input) { + AD::SetAdjIndex(Input_AdjIndices[iDim], Coord[iDim]); + } + else { + AD::SetAdjIndex(Output_AdjIndices[iDim], Coord[iDim]); + } + } +} + +void CPoint::SetAdjointSolution(const su2double *adj_sol) { + for (unsigned short iDim = 0; iDim < nDim; iDim++) { + AD::SetDerivative(Output_AdjIndices[iDim], SU2_TYPE::GetValue(adj_sol[iDim])); + } +} + +su2double CPoint::GetAdjointSolution(unsigned short iDim) { + return AD::GetDerivative(Input_AdjIndices[iDim]); +} + CEdge::CEdge(unsigned long val_iPoint, unsigned long val_jPoint, unsigned short val_nDim) : CDualGrid(val_nDim) { unsigned short iDim; diff --git a/Common/src/geometry_structure.cpp b/Common/src/geometry_structure.cpp index 0ffdbce65b4e..ab3b550e5a5e 100644 --- a/Common/src/geometry_structure.cpp +++ b/Common/src/geometry_structure.cpp @@ -2682,10 +2682,20 @@ void CGeometry::ComputeAirfoil_Section(su2double *Plane_P0, su2double *Plane_Nor void CGeometry::RegisterCoordinates(CConfig *config) { unsigned short iDim; unsigned long iPoint; + + bool input = true; for (iPoint = 0; iPoint < nPoint; iPoint++) { - for (iDim = 0; iDim < nDim; iDim++) { - AD::RegisterInput(node[iPoint]->GetCoord()[iDim]); + if(config->GetMultizone_Problem()) { + for (iDim = 0; iDim < nDim; iDim++) { + AD::RegisterInput_intIndexBased(node[iPoint]->GetCoord()[iDim]); + node[iPoint]->SetAdjIndices(input); + } + } + else { + for (iDim = 0; iDim < nDim; iDim++) { + AD::RegisterInput(node[iPoint]->GetCoord()[iDim]); + } } } } @@ -2695,8 +2705,18 @@ void CGeometry::RegisterOutput_Coordinates(CConfig *config){ unsigned long iPoint; for (iPoint = 0; iPoint < nPoint; iPoint++){ - for (iDim = 0; iDim < nDim; iDim++){ - AD::RegisterOutput(node[iPoint]->GetCoord()[iDim]); + if(config->GetMultizone_Problem()) { + for (iDim = 0; iDim < nDim; iDim++) { +// In case we are dealing with mesh adjoints the following two lines might be needed +// AD::RegisterOutput(node[iPoint]->GetCoord()[iDim]); +// node[iPoint]->SetAdjIndices(false); + AD::RegisterOutput(node[iPoint]->GetCoord()[iDim]); + } + } + else { + for (iDim = 0; iDim < nDim; iDim++) { + AD::RegisterOutput(node[iPoint]->GetCoord()[iDim]); + } } } } diff --git a/Common/src/grid_movement_structure.cpp b/Common/src/grid_movement_structure.cpp index 067e19a385be..1333f9f19383 100644 --- a/Common/src/grid_movement_structure.cpp +++ b/Common/src/grid_movement_structure.cpp @@ -138,7 +138,6 @@ void CVolumetricMovement::SetVolume_Deformation(CGeometry *geometry, CConfig *co su2double MinVolume, MaxVolume, NumError, Residual = 0.0, Residual_Init = 0.0; bool Screen_Output; - /*--- Retrieve number or iterations, tol, output, etc. from config ---*/ Smoothing_Iter = config->GetGridDef_Linear_Iter(); @@ -270,7 +269,6 @@ void CVolumetricMovement::SetVolume_Deformation(CGeometry *geometry, CConfig *co MaxIter = Smoothing_Iter - IterLinSol; IterLinSol = System.FGMRES_LinSolver(LinSysRes, LinSysSol, *mat_vec, *precond, NumError, MaxIter, &Residual, false, config); - Tot_Iter += IterLinSol; if ((rank == MASTER_NODE) && Screen_Output) { cout << " " << Tot_Iter << " " << Residual/Residual_Init << endl; } diff --git a/Common/src/linear_algebra/CSysSolve.cpp b/Common/src/linear_algebra/CSysSolve.cpp index e99f62a46cdf..1eb0e4fa3186 100644 --- a/Common/src/linear_algebra/CSysSolve.cpp +++ b/Common/src/linear_algebra/CSysSolve.cpp @@ -242,7 +242,7 @@ unsigned long CSysSolve::CG_LinSolver(const CSysVector & norm_r = r.norm(); norm0 = b.norm(); - if ( (norm_r < tol*norm0) || (norm_r < eps) ) { + if ((norm_r < tol*norm0) || (norm_r < eps)) { if (rank == MASTER_NODE) cout << "CSysSolve::ConjugateGradient(): system solved by initial guess." << endl; return 0; } @@ -391,7 +391,7 @@ unsigned long CSysSolve::FGMRES_LinSolver(const CSysVector::BCGSTAB_LinSolver(const CSysVector + * with selected contributions from the open-source community. + * + * The main research teams contributing to the current release are: + * - Prof. Juan J. Alonso's group at Stanford University. + * - Prof. Piero Colonna's group at Delft University of Technology. + * - Prof. Nicolas R. Gauger's group at Kaiserslautern University of Technology. + * - Prof. Alberto Guardone's group at Polytechnic University of Milan. + * - Prof. Rafael Palacios' group at Imperial College London. + * - Prof. Vincent Terrapon's group at the University of Liege. + * - Prof. Edwin van der Weide's group at the University of Twente. + * - Lab. of New Concepts in Aeronautics at Tech. Institute of Aeronautics. + * + * Copyright 2012-2019, Francisco D. Palacios, Thomas D. Economon, + * Tim Albring, and the SU2 contributors. + * + * SU2 is free software; you can redistribute it and/or + * modify it under the terms of the GNU Lesser General Public + * License as published by the Free Software Foundation; either + * version 2.1 of the License, or (at your option) any later version. + * + * SU2 is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + * Lesser General Public License for more details. + * + * You should have received a copy of the GNU Lesser General Public + * License along with SU2. If not, see . + */ + +#pragma once +#include "CMultizoneDriver.hpp" + +class CDiscAdjMultizoneDriver : public CMultizoneDriver { + +protected: + + /*! + * \brief Kinds of recordings (three different ones). + */ + enum class Kind_Tape { + FULL_TAPE, /*!< \brief Entire derivative information for a coupled adjoint + solution update. */ + OBJECTIVE_FUNCTION_TAPE, /*!< \brief Record only the dependence of the objective function + w.r.t. solver variables (from all zones). */ + ZONE_SPECIFIC_TAPE /*!< \brief Record only the dependence of the solution update in a + specified zone w.r.t. solver variables (from all zones). */ + }; + + /*! + * \brief Position markers within a tape. + */ + enum Tape_Positions { + START = 0, /*!< \brief Beginning of the tape. */ + REGISTERED = 1, /*!< \brief Solver variables are registered on the tape. */ + DEPENDENCIES = 2, /*!< \brief Derived values (e.g. gradients) are set. */ + OBJECTIVE_FUNCTION = 3, /*!< \brief Objective function is set. */ + TRANSFER = 4, /*!< \brief Solution data is transferred between coupled solvers of + different physical zones (giving cross contributions). */ + ITERATION_READY = 5 /*!< \brief Until here, all derivative information is gathered so + that it can be connected to a solver update evaluation. */ + }; + + int RecordingState; /*!< \brief The kind of recording that the tape currently holds. */ + bool retape; /*!< \brief Boolean whether a full tape can be kept in memory. */ + + su2double ObjFunc; /*!< \brief Value of the objective function. */ + int ObjFunc_Index; /*!< \brief Index of the value of the objective function. */ + + CIteration*** direct_iteration; /*!< \brief Array of pointers to the direct iterations. */ + COutput** direct_output; /*!< \brief Array of pointers to the direct outputs. */ + unsigned short* direct_nInst; /*!< \brief Total number of instances in the direct problem. */ + unsigned short* nInnerIter; /*!< \brief Number of inner iterations for each zone. */ + + +public: + + /*! + * \brief Constructor of the class. + * \param[in] confFile - Configuration file name. + * \param[in] val_nZone - Total number of zones. + * \param[in] MPICommunicator - MPI communicator for SU2. + */ + CDiscAdjMultizoneDriver(char* confFile, + unsigned short val_nZone, + SU2_Comm MPICommunicator); + + /*! + * \brief Destructor of the class. + */ + ~CDiscAdjMultizoneDriver(void); + + /*! + * \brief [Overload] Launch the computation for discrete adjoint multizone problems. + */ + void StartSolver() override; + + /*! + * \brief [Overload] Run an discrete adjoint update of all solvers within multiple zones. + */ + void Run() override; + + /*! + * \brief Record one iteration of a flow iteration in within multiple zones. + * \param[in] kind_recording - Kind of variables with regard to which we are recording. + * \param[in] tape_type - indicator which part of a solution update will be recorded + * \param[in] record_zone - zone where solution update will be recorded + */ + void SetRecording(unsigned short kind_recording, Kind_Tape tape_type, unsigned short record_zone); + + /*! + * \brief Run one direct iteration in a zone. + * \param[in] iZone - Zone in which we run an iteration. + * \param[in] kind_recording - Kind of variables with regard to which we are recording. + */ + void DirectIteration(unsigned short iZone, unsigned short kind_recording); + + /*! + * \brief Set the objective function. + * \param[in] kind_recording - Kind of variables with regard to which we are recording. + */ + void SetObjFunction(unsigned short kind_recording); + + /*! + * \brief Initialize the adjoint value of the objective function. + */ + void SetAdj_ObjFunction(); + + /*! + * \brief Summary of all routines to evaluate the adjoints in iZone. + * \param[in] iZone - Zone in which adjoints are evaluated depending on their (preceding) seeding. + */ + void ComputeAdjoints(unsigned short iZone); + + /*! + * \brief Add External_Old vector to Solution. + * \param[in] iZone - Zone where data between solvers is transferred. + */ + void Add_ExternalOld_To_Solution(unsigned short iZone); + + /*! + * \brief Sets External to zero. + */ + void SetExternal_Zero(void); + + /*! + * \brief Set External_Old to External. + */ + void Set_OldExternal(void); + + /*! + * \brief Add Solution vector to External. + * \param[in] iZone - Zone where data between solvers is transferred. + */ + void Add_Solution_To_External(unsigned short iZone); + + /*! + * \brief Add Solution vector to External_Old. + * \param[in] iZone - Zone where data between solvers is transferred. + */ + void Add_Solution_To_ExternalOld(unsigned short iZone); + + /*! + * \brief Saves the current (adjoint) Solution vector to Solution_BGS_k. + * \param[in] iZone - Zone where data between solvers is transferred. + */ + void Set_BGSSolution(unsigned short iZone); + + /*! + * \brief Compute BGS residuals. + * \param[in] iZone - Zone where solver residuals are computed. + */ + void SetResidual_BGS(unsigned short iZone); +}; diff --git a/SU2_CFD/include/output/CHeatOutput.hpp b/SU2_CFD/include/output/CHeatOutput.hpp index 8a41c49f45b8..81cad15dcd51 100644 --- a/SU2_CFD/include/output/CHeatOutput.hpp +++ b/SU2_CFD/include/output/CHeatOutput.hpp @@ -59,6 +59,12 @@ class CHeatOutput final: public COutput { */ ~CHeatOutput(void) override; + /*! + * \brief Set the available history output fields + * \param[in] config - Definition of the particular problem. + */ + void SetHistoryOutputFields(CConfig *config) override; + /*! * \brief Load the history output field values * \param[in] config - Definition of the particular problem. @@ -78,11 +84,17 @@ class CHeatOutput final: public COutput { * \param[in] solver - The container holding all solution data. * \param[in] iPoint - Index of the point. */ - void LoadVolumeData(CConfig *config, CGeometry *geometry, CSolver **solver, unsigned long iPoint) override; - + void LoadVolumeData(CConfig *config, CGeometry *geometry, CSolver **solver, unsigned long iPoint) override; + /*! - * \brief Set the available history output fields + * \brief LoadSurfaceData * \param[in] config - Definition of the particular problem. + * \param[in] geometry - Geometrical definition of the problem. + * \param[in] solver - The container holding all solution data. + * \param[in] iPoint - Index of the point. + * \param[in] iMarker - Index of the surface marker. + * \param[in] iVertex - Index of the vertex on the marker. */ - void SetHistoryOutputFields(CConfig *config) override; + void LoadSurfaceData(CConfig *config, CGeometry *geometry, CSolver **solver, unsigned long iPoint, + unsigned short iMarker, unsigned long iVertex); }; diff --git a/SU2_CFD/include/output/COutput.hpp b/SU2_CFD/include/output/COutput.hpp index d9b3cb6be8b3..949b1afe546d 100644 --- a/SU2_CFD/include/output/COutput.hpp +++ b/SU2_CFD/include/output/COutput.hpp @@ -378,6 +378,10 @@ class COutput { su2double GetHistoryFieldValue(string field){ return historyOutput_Map[field].value; } + + su2double GetHistoryFieldValuePerSurface(string field, unsigned short iMarker){ + return historyOutputPerSurface_Map[field][iMarker].value; + } /*! * \brief Get a vector with all output fields in a particular group @@ -742,8 +746,8 @@ class COutput { * \param[in] config - Definition of the particular problem. */ inline virtual void SetVolumeOutputFields(CConfig *config){} - - + + /*! * \brief Load the history output field values * \param[in] config - Definition of the particular problem. diff --git a/SU2_CFD/include/solver_structure.hpp b/SU2_CFD/include/solver_structure.hpp index 89d5061bcc5c..553a8debb405 100644 --- a/SU2_CFD/include/solver_structure.hpp +++ b/SU2_CFD/include/solver_structure.hpp @@ -83,6 +83,7 @@ class CSolver { protected: int rank, /*!< \brief MPI Rank. */ size; /*!< \brief MPI Size. */ + bool adjoint; /*!< \brief Boolean to determine whether solver is initialized as a direct or an adjoint solver. */ unsigned short MGLevel; /*!< \brief Multigrid level of this solver object. */ unsigned short IterLinSolver; /*!< \brief Linear solver iterations. */ unsigned short nVar, /*!< \brief Number of variables of the problem. */ @@ -278,6 +279,12 @@ class CSolver { */ virtual void SetNondimensionalization(CConfig *config, unsigned short iMesh); + /*! + * \brief Get information whether the initialization is an adjoint solver or not. + * \return TRUE means that it is an adjoint solver. + */ + bool GetAdjoint(void); + /*! * \brief Compute the pressure at the infinity. * \return Value of the pressure at the infinity. @@ -487,7 +494,14 @@ class CSolver { * \return Pointer to the location (x, y, z) of the biggest residual for the variable val_var. */ su2double* GetPoint_Max_Coord_BGS(unsigned short val_var); - + + /*! + * \brief Set the value of the RMS residual respective solution. + * \param[in] geometry - Geometrical definition of the problem. + * \param[in] config - Definition of the particular problem. + */ + void SetResidual_Solution(CGeometry *geometry, CConfig *config); + /*! * \brief Set Value of the residual due to the Geometric Conservation Law (GCL) for steady rotating frame problems. * \param[in] geometry - Geometrical definition of the problem. @@ -514,7 +528,49 @@ class CSolver { * \param[in] config - Definition of the particular problem. */ void SetAuxVar_Surface_Gradient(CGeometry *geometry, CConfig *config); + + /*! + * \brief Set the solution vector to solution in Solution_Old. + * \param[in] geometry - The geometrical definition of the problem. + */ + void SetSolution_Old(CGeometry *geometry); + /*! + * \brief Add External_Old to Solution vector. + * \param[in] geometry - The geometrical definition of the problem. + */ + void Add_ExternalOld_To_Solution(CGeometry *geometry); + + /*! + * \brief Set the Solution vector to zero. + * \param[in] geometry - The geometrical definition of the problem. + */ + void SetSolution_Zero(CGeometry *geometry); + + /*! + * \brief Set the External vector to zero. + * \param[in] geometry - The geometrical definition of the problem. + */ + void SetExternal_Zero(CGeometry *geometry); + + /*! + * \brief Add the current Solution vector to External. + * \param[in] geometry - The geometrical definition of the problem. + */ + void Add_Solution_To_External(CGeometry *geometry); + + /*! + * \brief Add the current Solution vector to External_Old. + * \param[in] geometry - The geometrical definition of the problem. + */ + void Add_Solution_To_ExternalOld(CGeometry *geometry); + + /*! + * \brief Set External_Old to External. + * \param[in] geometry - The geometrical definition of the problem. + */ + void Set_OldExternal(CGeometry *geometry); + /*! * \brief Compute the Green-Gauss gradient of the solution. * \param[in] geometry - Geometrical definition of the problem. @@ -11325,8 +11381,8 @@ class CWaveSolver : public CSolver { class CHeatSolverFVM : public CSolver { protected: unsigned short nVarFlow, nMarker, CurrentMesh; - su2double **HeatFlux, *Surface_HF, Total_HeatFlux, AllBound_HeatFlux, - *AvgTemperature, Total_AvgTemperature, AllBound_AvgTemperature, + su2double **HeatFlux, *HeatFlux_per_Marker, *Surface_HF, Total_HeatFlux, AllBound_HeatFlux, + *AverageT_per_Marker, Total_AverageT, AllBound_AverageT, *Primitive, *Primitive_Flow_i, *Primitive_Flow_j, *Surface_Areas, Total_HeatFlux_Areas, Total_HeatFlux_Areas_Monitor; su2double ***ConjugateVar, ***InterfaceVar; @@ -12784,7 +12840,7 @@ class CDiscAdjSolver : public CSolver { * \param[in] val_update_geo - Flag for updating coords and grid velocity. */ void LoadRestart(CGeometry **geometry, CSolver ***solver, CConfig *config, int val_iter, bool val_update_geo); - + void ComputeResidual_Multizone(CGeometry *geometry, CConfig *config); }; diff --git a/SU2_CFD/include/solver_structure.inl b/SU2_CFD/include/solver_structure.inl index 39adf5cc107d..2f0ce508955d 100644 --- a/SU2_CFD/include/solver_structure.inl +++ b/SU2_CFD/include/solver_structure.inl @@ -41,6 +41,8 @@ inline void CSolver::SetIterLinSolver(unsigned short val_iterlinsolver) { IterLi inline void CSolver::SetNondimensionalization(CConfig *config, unsigned short iMesh) { } +inline bool CSolver::GetAdjoint(void) { return adjoint; } + inline unsigned short CSolver::GetIterLinSolver(void) { return IterLinSolver; } inline su2double CSolver::GetCSensitivity(unsigned short val_marker, unsigned long val_vertex) { return 0; } @@ -2263,7 +2265,7 @@ inline su2double CHeatSolverFVM::GetTotal_HeatFlux() { return Total_HeatFlux; } inline su2double CHeatSolverFVM::GetHeatFlux(unsigned short val_marker, unsigned long val_vertex) { return HeatFlux[val_marker][val_vertex]; } -inline su2double CHeatSolverFVM::GetTotal_AvgTemperature() { return Total_AvgTemperature; } +inline su2double CHeatSolverFVM::GetTotal_AvgTemperature() { return Total_AverageT; } inline su2double CHeatSolverFVM::GetConjugateHeatVariable(unsigned short val_marker, unsigned long val_vertex, unsigned short pos_var) { return ConjugateVar[val_marker][val_vertex][pos_var]; } diff --git a/SU2_CFD/include/variables/CDiscAdjVariable.hpp b/SU2_CFD/include/variables/CDiscAdjVariable.hpp index 8fe1a4a04821..da57d66d41de 100644 --- a/SU2_CFD/include/variables/CDiscAdjVariable.hpp +++ b/SU2_CFD/include/variables/CDiscAdjVariable.hpp @@ -157,6 +157,49 @@ class CDiscAdjVariable : public CVariable { Solution_Geometry[val_var] = val_solution_geometry; } + /*! + * \brief Get the external contribution to a solution. + * \return Pointer to the External array. + */ + inline su2double* Get_External(void) const { return External; } + + /*! + * \brief Set external contributions to zero. + */ + inline void SetExternalZero(void) { + for (unsigned short iVar = 0; iVar < nVar; iVar++) { + External[iVar] = 0.0; + } + } + + /*! + * \brief Add a value to the External vector. + * \param[in] val_sol - vector that has to be added component-wise + */ + inline void Add_External(const su2double* val_sol) { + for (unsigned short iVar = 0; iVar < nVar; iVar++) { + External[iVar] += val_sol[iVar]; + } + } + + /*! + * \brief Add a value to the old External vector. + * \param[in] val_solution - Value that we want to add to the solution. + */ + inline void Add_ExternalOld(const su2double *val_sol) { + for (unsigned short iVar = 0; iVar < nVar; iVar++) { + External_Old[iVar] = External_Old[iVar] + val_sol[iVar]; + } + } + + /*! + * \brief Set old External to the value of the current variables. + */ + inline void Set_OldExternal() { + for (unsigned short iVar = 0; iVar < nVar; iVar++) + External_Old[iVar] = External[iVar]; + } + /*! * \brief A virtual member. Get the geometry solution. * \param[in] val_var - Index of the variable. diff --git a/SU2_CFD/include/variables/CVariable.hpp b/SU2_CFD/include/variables/CVariable.hpp index b131948923cb..010d91443aa1 100644 --- a/SU2_CFD/include/variables/CVariable.hpp +++ b/SU2_CFD/include/variables/CVariable.hpp @@ -60,7 +60,9 @@ class CVariable { protected: su2double *Solution, /*!< \brief Solution of the problem. */ - *Solution_Old; /*!< \brief Old solution of the problem R-K. */ + *Solution_Old, /*!< \brief Old solution of the problem R-K. */ + *External, /*!< \brief External (outer) contribution in discrete adjoint multizone problems. */ + *External_Old; /*!< \brief Old external (outer) contribution in discrete adjoint multizone problems. */ bool Non_Physical; /*!< \brief Non-physical points in the solution (force first order). */ su2double *Solution_time_n, /*!< \brief Solution of the problem at time n for dual-time stepping technique. */ *Solution_time_n1; /*!< \brief Solution of the problem at time n-1 for dual-time stepping technique. */ @@ -92,6 +94,11 @@ class CVariable { note that this variable cannnot be static, it is possible to have different number of nVar in the same problem. */ su2double *Solution_Adj_Old; /*!< \brief Solution of the problem in the previous AD-BGS iteration. */ + + int *Input_AdjIndices, /*!< \brief Indices of Solution variables in the adjoint vector. */ + *Output_AdjIndices; /*!< \brief Indices of Solution variables in the adjoint vector after having been updated. */ + + /*--- Old solution container for BGS iterations ---*/ su2double *Solution_BGS_k; @@ -317,6 +324,16 @@ class CVariable { Solution[val_var] = Solution_Old[val_var] + val_solution; } + /*! + * \brief Add a value to the solution. + * \param[in] val_solution - Value that we want to add to the solution. + */ + inline void AddSolution(su2double *val_solution) { + for (unsigned short iVar = 0; iVar < nVar; iVar++) { + Solution[iVar] = Solution[iVar] + val_solution[iVar]; + } + } + /*! * \brief A virtual member. * \param[in] val_var - Index of the variable. @@ -361,6 +378,11 @@ class CVariable { */ inline virtual void SetSolution_New(void) {} + /*! + * \brief A virtual member. + */ + inline virtual void SetExternalZero(void) {} + /*! * \brief A virtual member. * \param[in] val_var - Number of the variable. @@ -368,6 +390,21 @@ class CVariable { */ inline virtual void AddSolution_New(unsigned short val_var, su2double val_solution) {} + /*! + * \brief A virtual member. + */ + inline virtual void Add_External(const su2double* val_sol) {} + + /*! + * \brief A virtual member. + */ + inline virtual void Add_ExternalOld(const su2double* val_sol) {} + + /*! + * \brief A virtual member. + */ + inline virtual void Set_OldExternal(void) {} + /*! * \brief Add a value to the solution, clipping the values. * \param[in] val_var - Index of the variable. @@ -411,6 +448,12 @@ class CVariable { */ inline su2double *GetSolution_Old(void) {return Solution_Old; } + /*! + * \brief Get the old external contributions of the problem. + * \return Pointer to the External_Old vector. + */ + inline su2double *GetSolution_ExternalOld(void) {return External_Old; } + /*! * \brief Get the solution at time n. * \return Pointer to the solution (at time n) vector. @@ -745,6 +788,11 @@ class CVariable { */ inline su2double GetSolution_Min(unsigned short val_var) {return Solution_Min[val_var]; } + /*! + * \brief A virtual member. + */ + inline virtual su2double *Get_External(void) const { return NULL; } + /*! * \brief Get the value of the preconditioner Beta. * \return Value of the low Mach preconditioner variable Beta @@ -2515,6 +2563,36 @@ class CVariable { AD::RegisterOutput(Solution[iVar]);} } + /*! + * \brief Register the variables in the solution array as input/output variable. + * \param[in] input - input or output variables. + */ + inline void RegisterSolution_intIndexBased(bool input) { + if (input) { + for (unsigned short iVar = 0; iVar < nVar; iVar++) + AD::RegisterInput_intIndexBased(Solution[iVar]); + } + else { + for (unsigned short iVar = 0; iVar < nVar; iVar++) + AD::RegisterOutput(Solution[iVar]); + } + } + + /*! + * \brief Saving the adjoint vector position with respect to the solution variables. + * \param[in] input - input or output variables. + */ + inline void SetAdjIndices(bool input) { + for (unsigned short iVar = 0; iVar < nVar; iVar++) { + if (input) { + AD::SetAdjIndex(Input_AdjIndices[iVar], Solution[iVar]); + } + else { + AD::SetAdjIndex(Output_AdjIndices[iVar], Solution[iVar]); + } + } + } + /*! * \brief Register the variables in the solution_time_n array as input/output variable. */ @@ -2540,6 +2618,16 @@ class CVariable { SU2_TYPE::SetDerivative(Solution[iVar], SU2_TYPE::GetValue(adj_sol[iVar])); } + /*! + * \brief Set the adjoint values of the solution. + * \param[in] adj_sol - The adjoint values of the solution. + */ + inline void SetAdjointSolution_intIndexBased(su2double *adj_sol) { + for (unsigned short iVar = 0; iVar < nVar; iVar++) { + AD::SetDerivative(Output_AdjIndices[iVar], SU2_TYPE::GetValue(adj_sol[iVar])); + } + } + /*! * \brief Get the adjoint values of the solution. * \param[out] adj_sol - The adjoint values of the solution. @@ -2549,6 +2637,16 @@ class CVariable { adj_sol[iVar] = SU2_TYPE::GetDerivative(Solution[iVar]); } + /*! + * \brief Get the adjoint values of the solution. + * \param[in] adj_sol - The adjoint values of the solution. + */ + inline void GetAdjointSolution_intIndexBased(su2double *adj_sol) { + for (unsigned short iVar = 0; iVar < nVar; iVar++) { + adj_sol[iVar] = AD::GetDerivative(Input_AdjIndices[iVar]); + } + } + /*! * \brief Set the adjoint values of the solution at time n. * \param[in] adj_sol - The adjoint values of the solution. diff --git a/SU2_CFD/obj/Makefile.am b/SU2_CFD/obj/Makefile.am index 1860eacd7947..ada7e4097442 100644 --- a/SU2_CFD/obj/Makefile.am +++ b/SU2_CFD/obj/Makefile.am @@ -69,6 +69,7 @@ libSU2Core_sources = ../src/definition_structure.cpp \ ../src/drivers/CMultizoneDriver.cpp \ ../src/drivers/CSinglezoneDriver.cpp \ ../src/drivers/CDiscAdjSinglezoneDriver.cpp \ + ../src/drivers/CDiscAdjMultizoneDriver.cpp \ ../src/drivers/CDriver.cpp \ ../src/drivers/CDummyDriver.cpp \ ../src/iteration_structure.cpp \ diff --git a/SU2_CFD/src/SU2_CFD.cpp b/SU2_CFD/src/SU2_CFD.cpp index 7f94c3c0125f..17d5b5c7870f 100644 --- a/SU2_CFD/src/SU2_CFD.cpp +++ b/SU2_CFD/src/SU2_CFD.cpp @@ -109,26 +109,37 @@ int main(int argc, char *argv[]) { if (!dry_run){ - if (((config->GetSinglezone_Driver() || (nZone == 1 && config->GetDiscrete_Adjoint())) - && config->GetTime_Marching() != HARMONIC_BALANCE && (!turbo)) || (turbo && config->GetDiscrete_Adjoint())) { - + if ((!config->GetMultizone_Problem()) + && (config->GetTime_Marching() != HARMONIC_BALANCE && (!turbo)) + || (turbo && config->GetDiscrete_Adjoint())) { /*--- Single zone problem: instantiate the single zone driver class. ---*/ if (nZone > 1 ) { SU2_MPI::Error("The required solver doesn't support multizone simulations", CURRENT_FUNCTION); } - if (config->GetDiscrete_Adjoint()) + if (config->GetDiscrete_Adjoint()) { + driver = new CDiscAdjSinglezoneDriver(config_file_name, nZone, MPICommunicator); + } else driver = new CSinglezoneDriver(config_file_name, nZone, MPICommunicator); } else if (config->GetMultizone_Problem() && !turbo && !fsi) { - /*--- Multizone Driver. ---*/ + /*--- Multizone Drivers. ---*/ + + if (config->GetDiscrete_Adjoint()) { + + driver = new CDiscAdjMultizoneDriver(config_file_name, nZone, MPICommunicator); + + } + else { driver = new CMultizoneDriver(config_file_name, nZone, MPICommunicator); + + } } else if (config->GetTime_Marching() == HARMONIC_BALANCE) { diff --git a/SU2_CFD/src/drivers/CDiscAdjMultizoneDriver.cpp b/SU2_CFD/src/drivers/CDiscAdjMultizoneDriver.cpp new file mode 100644 index 000000000000..472e2d8eddd6 --- /dev/null +++ b/SU2_CFD/src/drivers/CDiscAdjMultizoneDriver.cpp @@ -0,0 +1,866 @@ +/*! + * \file CDiscAdjMultizoneDriver.cpp + * \brief The main subroutines for driving adjoint multi-zone problems + * \author O. Burghardt, T. Albring, R. Sanchez + * \version 6.2.0 "Falcon" + * + * The current SU2 release has been coordinated by the + * SU2 International Developers Society + * with selected contributions from the open-source community. + * + * The main research teams contributing to the current release are: + * - Prof. Juan J. Alonso's group at Stanford University. + * - Prof. Piero Colonna's group at Delft University of Technology. + * - Prof. Nicolas R. Gauger's group at Kaiserslautern University of Technology. + * - Prof. Alberto Guardone's group at Polytechnic University of Milan. + * - Prof. Rafael Palacios' group at Imperial College London. + * - Prof. Vincent Terrapon's group at the University of Liege. + * - Prof. Edwin van der Weide's group at the University of Twente. + * - Lab. of New Concepts in Aeronautics at Tech. Institute of Aeronautics. + * + * Copyright 2012-2018, Francisco D. Palacios, Thomas D. Economon, + * Tim Albring, and the SU2 contributors. + * + * SU2 is free software; you can redistribute it and/or + * modify it under the terms of the GNU Lesser General Public + * License as published by the Free Software Foundation; either + * version 2.1 of the License, or (at your option) any later version. + * + * SU2 is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + * Lesser General Public License for more details. + * + * You should have received a copy of the GNU Lesser General Public + * License along with SU2. If not, see . + */ + +#include "../include/drivers/CDiscAdjMultizoneDriver.hpp" + +CDiscAdjMultizoneDriver::CDiscAdjMultizoneDriver(char* confFile, + unsigned short val_nZone, + SU2_Comm MPICommunicator) + + : CMultizoneDriver(confFile, val_nZone, MPICommunicator) { + + retape = !config_container[ZONE_0]->GetFull_Tape(); + + RecordingState = NONE; + + direct_nInst = new unsigned short[nZone]; + nInnerIter = new unsigned short[nZone]; + + + for (iZone = 0; iZone < nZone; iZone++) { + + direct_nInst[iZone] = 1; + nInnerIter[iZone] = config_container[iZone]->GetnInner_Iter(); + } + + direct_iteration = new CIteration**[nZone]; + direct_output = new COutput*[nZone]; + + for (iZone = 0; iZone < nZone; iZone++) { + + direct_iteration[iZone] = new CIteration*[direct_nInst[iZone]]; + + for(iInst = 0; iInst < direct_nInst[iZone]; iInst++) { + + switch (config_container[iZone]->GetKind_Solver()) { + + case DISC_ADJ_EULER: case DISC_ADJ_NAVIER_STOKES: case DISC_ADJ_RANS: + direct_iteration[iZone][iInst] = new CFluidIteration(config_container[iZone]); + direct_output[iZone] = new CFlowCompOutput(config_container[iZone], nDim); + break; + case DISC_ADJ_INC_EULER: case DISC_ADJ_INC_NAVIER_STOKES: case DISC_ADJ_INC_RANS: + direct_iteration[iZone][iInst] = new CFluidIteration(config_container[iZone]); + direct_output[iZone] = new CFlowIncOutput(config_container[iZone], nDim); + break; + case DISC_ADJ_HEAT: + direct_iteration[iZone][iInst] = new CHeatIteration(config_container[iZone]); + direct_output[iZone] = new CHeatOutput(config_container[iZone], nDim); + break; + case DISC_ADJ_FEM: + SU2_MPI::Error("There is no multizone discrete adjoint functionality for problems including elasticity yet.", + CURRENT_FUNCTION); + direct_iteration[iZone][iInst] = new CFEAIteration(config_container[iZone]); + direct_output[iZone] = new CElasticityOutput(config_container[iZone], nDim); + break; + + default: + SU2_MPI::Error("There is no discrete adjoint functionality for one of the specified solvers yet.", + CURRENT_FUNCTION); + } + } + + direct_output[iZone]->PreprocessHistoryOutput(config_container[iZone], false); + } +} + +CDiscAdjMultizoneDriver::~CDiscAdjMultizoneDriver(){ + + for (iZone = 0; iZone < nZone; iZone++){ + for (iInst = 0; iInst < direct_nInst[iZone]; iInst++){ + delete direct_iteration[iZone][iInst]; + } + delete [] direct_iteration[iZone]; + } + + delete[] direct_iteration; + delete[] direct_nInst; + +} + +void CDiscAdjMultizoneDriver::StartSolver() { + + /*--- Main external loop of the solver. Runs for the number of time steps required. ---*/ + + if (rank == MASTER_NODE) + cout << endl <<"------------------------------ Begin Solver -----------------------------" << endl; + + if (rank == MASTER_NODE){ + cout << endl << "Simulation Run using the Discrete Adjoint Multizone Driver" << endl; + + if (driver_config->GetTime_Domain()) + SU2_MPI::Error("The discrete adjoint multizone driver is not ready for unsteady computations yet.", + CURRENT_FUNCTION); + } + + for (iZone = 0; iZone < nZone; iZone++){ + + /*--- Set the value of the external iteration to TimeIter. -------------------------------------*/ + /*--- TODO: This should be generalised for an homogeneous criteria throughout the code. --------*/ + config_container[iZone]->SetTimeIter(0); + + } + + /*--- We directly start the (steady-state) discrete adjoint computation. ---*/ + + Run(); + + /*--- Output the solution in files. ---*/ + + Output(TimeIter); + +} + +void CDiscAdjMultizoneDriver::Run() { + + bool checkSensitivity = false; + unsigned short jZone = 0, + wrt_sol_freq = config_container[ZONE_0]->GetVolume_Wrt_Freq(); + unsigned long nIter = 0, + iOuter_Iter = 0, + iInner_Iter = 0; + + nIter = driver_config->GetnOuter_Iter(); + + for (iZone = 0; iZone < nZone; iZone++) { + + iteration_container[iZone][INST_0]->Preprocess(output_container[iZone], integration_container, geometry_container, + solver_container, numerics_container, config_container, surface_movement, + grid_movement, FFDBox, iZone, INST_0); + + /*--- Set BGS_Solution_k to Solution. ---*/ + + Set_BGSSolution(iZone); + } + + /*--- Loop over the number of outer iterations. ---*/ + + for (iOuter_Iter = 0; iOuter_Iter < nIter; iOuter_Iter++) { + + for (iZone = 0; iZone < nZone; iZone++) { + config_container[iZone]->SetOuterIter(iOuter_Iter); + driver_config->SetOuterIter(iOuter_Iter); + } + + + /*--- For the adjoint iteration we need the derivatives of the iteration function with + * respect to the state (and possibly the mesh coordinate) variables. + * Since these derivatives do not change in the steady state case we only have to record + * if the current recording is different from them. + * + * To set the tape appropriately, the following recording methods are provided: + * (1) NONE: All information from a previous recording is removed. + * (2) FLOW_CONS_VARS: State variables of all solvers in a zone as input. + * (3) MESH_COORDS: Mesh coordinates as input. + * (4) COMBINED: Mesh coordinates and state variables as input. + * + * By default, all (state and mesh coordinate variables) will be declared as output, + * since it does not change the computational effort. ---*/ + + + /*--- If we want to set up zone-specific tapes later on, + * we just record the objective function section here. + * If not, the whole tape of a coupled run will be created. ---*/ + + if (retape) { + SetRecording(NONE, Kind_Tape::FULL_TAPE, ZONE_0); + SetRecording(FLOW_CONS_VARS, Kind_Tape::OBJECTIVE_FUNCTION_TAPE, ZONE_0); + } + else if (RecordingState != FLOW_CONS_VARS) { + + SetRecording(NONE, Kind_Tape::FULL_TAPE, ZONE_0); + SetRecording(FLOW_CONS_VARS, Kind_Tape::FULL_TAPE, ZONE_0); + } + + + SetExternal_Zero(); + + /*-- Start loop over zones. ---*/ + + for (iZone = 0; iZone < nZone; iZone++) { + + if (retape) { + SetRecording(NONE, Kind_Tape::FULL_TAPE, ZONE_0); + SetRecording(FLOW_CONS_VARS, Kind_Tape::ZONE_SPECIFIC_TAPE, iZone); + } + + /*--- Evaluate the objective function gradient w.r.t. to solution contributions from iZone. + * (We always initialize from Solution_BGS_k and extract to Solution.) ---*/ + + AD::ClearAdjoints(); + + SetAdj_ObjFunction(); + + AD::ComputeAdjoint(OBJECTIVE_FUNCTION, START); + + iteration_container[iZone][INST_0]->Iterate(output_container[iZone], integration_container, geometry_container, + solver_container, numerics_container, config_container, + surface_movement, grid_movement, FFDBox, iZone, INST_0); + + /*--- Add the objective function contribution to the external contributions for solvers in iZone. ---*/ + + Add_Solution_To_ExternalOld(iZone); + + /*--- Inner loop to allow for multiple adjoint updates with respect to solvers in iZone. ---*/ + + for (unsigned short iInnerIter = 0; iInnerIter < nInnerIter[iZone]; iInnerIter++) { + + /*--- Evaluate the tape section belonging to solvers in iZone. ---*/ + + ComputeAdjoints(iZone); + + /*--- Extracting adjoints for solvers in iZone w.r.t. to outputs in iZone (diagonal part). ---*/ + + iteration_container[iZone][INST_0]->Iterate(output_container[iZone], integration_container, geometry_container, + solver_container, numerics_container, config_container, + surface_movement, grid_movement, FFDBox, iZone, INST_0); + + /*--- Add off-diagonal contribution (including the OF gradient) to Solution for next inner evaluation. ---*/ + + Add_ExternalOld_To_Solution(iZone); + } + + + /*--- Off-diagonal (coupling term) update. ---*/ + + for (jZone = 0; jZone < nZone; jZone++) { + + if (jZone != iZone) { + + /*--- Extracting adjoints for solvers in jZone w.r.t. to the output of all solvers in iZone, + * that is, for the cases iZone != jZone we are evaluating cross derivatives between zones. ---*/ + + iteration_container[jZone][INST_0]->Iterate(output_container[jZone], integration_container, geometry_container, + solver_container, numerics_container, config_container, + surface_movement, grid_movement, FFDBox, jZone, INST_0); + + /*--- Add the cross derivatives from iZone<-jZone dependencies to the External vector. ---*/ + + Add_Solution_To_External(jZone); + + } + } + + // TODO: Add an option to _already_ update off-diagonals here (i.e., in the zone-loop) + + /*--- Compute residual from Solution and Solution_BGS_k. ---*/ + + SetResidual_BGS(iZone); + + /*--- Save Solution to Solution_BGS_k for a next outer iteration. + * (Solution might be overwritten when entering another zone because of cross derivatives.) ---*/ + + Set_BGSSolution(iZone); + } + + /*--- Now all coupling terms are summed up, set External_Old to External for next outer iteration. ---*/ + + Set_OldExternal(); + + /*--- Print out the convergence data to screen and history file ---*/ + + for (iZone = 0; iZone < nZone; iZone++) { + output_container[iZone]->SetHistory_Output(geometry_container[iZone][INST_0][MESH_0], solver_container[iZone][INST_0][MESH_0], + config_container[iZone], config_container[iZone]->GetTimeIter(), iOuter_Iter, iInner_Iter); + } + + driver_output->SetMultizoneHistory_Output(output_container, config_container, driver_config, + driver_config->GetTimeIter(), driver_config->GetOuterIter()); + + /*--- Check for convergence. ---*/ + + StopCalc = driver_output->GetConvergence(); + + /*--- Clear the stored adjoint information to be ready for a new evaluation. ---*/ + + AD::ClearAdjoints(); + + /*--- Compute the geometrical sensitivities and write them to file. ---*/ + + checkSensitivity = ((iOuter_Iter+1 >= nIter) || + (iOuter_Iter % wrt_sol_freq == 0)); + + if (checkSensitivity || StopCalc){ + + /*--- SetRecording stores the computational graph on one iteration of the direct problem. Calling it with NONE + * as argument ensures that all information from a previous recording is removed. ---*/ + + SetRecording(NONE, Kind_Tape::FULL_TAPE, ZONE_0); + + /*--- Store the computational graph of one direct iteration with the mesh coordinates as input. ---*/ + + SetRecording(MESH_COORDS, Kind_Tape::FULL_TAPE, ZONE_0); + + /*--- Initialize the adjoint of the output variables of the iteration with the adjoint solution + * of the current iteration. The values are passed to the AD tool. ---*/ + + for (iZone = 0; iZone < nZone; iZone++) { + + iteration_container[iZone][INST_0]->InitializeAdjoint(solver_container, geometry_container, + config_container, iZone, INST_0); + } + + /*--- Initialize the adjoint of the objective function with 1.0. ---*/ + + SetAdj_ObjFunction(); + + /*--- Interpret the stored information by calling the corresponding routine of the AD tool. ---*/ + + AD::ComputeAdjoint(); + + /*--- Extract the computed sensitivity values. ---*/ + + for (iZone = 0; iZone < nZone; iZone++) { + + switch (config_container[iZone]->GetKind_Solver()) { + + case DISC_ADJ_EULER: case DISC_ADJ_NAVIER_STOKES: case DISC_ADJ_RANS: + case DISC_ADJ_INC_EULER: case DISC_ADJ_INC_NAVIER_STOKES: case DISC_ADJ_INC_RANS: + solver_container[iZone][INST_0][MESH_0][ADJFLOW_SOL]->SetSensitivity(geometry_container[iZone][INST_0][MESH_0], + NULL, config_container[iZone]); + break; + case DISC_ADJ_HEAT: + solver_container[iZone][INST_0][MESH_0][ADJHEAT_SOL]->SetSensitivity(geometry_container[iZone][INST_0][MESH_0], + NULL, config_container[iZone]); + break; + case DISC_ADJ_FEM: + solver_container[iZone][INST_0][MESH_0][ADJFEA_SOL]->SetSensitivity(geometry_container[iZone][INST_0][MESH_0], + NULL, config_container[iZone]); + break; + + default: + cout << "WARNING: Sensitivities not set for one of the specified discrete adjoint solvers!" << endl; + break; + } + } + + /*--- Clear the stored adjoint information to be ready for a new evaluation. ---*/ + + AD::ClearAdjoints(); + + for (iZone = 0; iZone < nZone; iZone++) { + + output_container[iZone]->SetResult_Files(geometry_container[iZone][INST_0][MESH_0], + config_container[iZone], + solver_container[iZone][INST_0][MESH_0], iOuter_Iter, StopCalc); + } + } + + if (StopCalc) break; + } +} + +void CDiscAdjMultizoneDriver::SetRecording(unsigned short kind_recording, Kind_Tape tape_type, unsigned short record_zone) { + + unsigned short iZone, jZone, iSol, UpdateMesh; + unsigned long ExtIter = 0; + bool DeformMesh = false; + + AD::Reset(); + + /*--- Prepare for recording by resetting the flow solution to the initial converged solution---*/ + + for(iZone = 0; iZone < nZone; iZone++) { + + for (iSol=0; iSol < MAX_SOLS; iSol++){ + if (solver_container[iZone][INST_0][MESH_0][iSol] != NULL) { + if (solver_container[iZone][INST_0][MESH_0][iSol]->GetAdjoint()) { + solver_container[iZone][INST_0][MESH_0][iSol]->SetRecording(geometry_container[iZone][INST_0][MESH_0], + config_container[iZone]); + } + } + } + } + + /*---Enable recording and register input of the flow iteration (conservative variables or node coordinates) --- */ + + if(kind_recording != NONE) { + + if (rank == MASTER_NODE && kind_recording == FLOW_CONS_VARS) { + cout << endl << "-------------------------------------------------------------------------" << endl; + cout << "Storing computational graph." << endl; + } + + AD::StartRecording(); + + AD::Push_TapePosition(); + + for (iZone = 0; iZone < nZone; iZone++) { + + iteration_container[iZone][INST_0]->RegisterInput(solver_container, geometry_container, + config_container, iZone, INST_0, kind_recording); + } + } + + AD::Push_TapePosition(); + + for (iZone = 0; iZone < nZone; iZone++) { + + iteration_container[iZone][INST_0]->SetDependencies(solver_container, geometry_container, numerics_container, + config_container, iZone, INST_0, kind_recording); + } + + AD::Push_TapePosition(); + + /*--- Extract the objective function and store it --- */ + + SetObjFunction(kind_recording); + + AD::Push_TapePosition(); + + if (tape_type != Kind_Tape::OBJECTIVE_FUNCTION_TAPE) { + + /*--- We do the communication here to not derive wrt updated boundary data. ---*/ + + for(iZone = 0; iZone < nZone; iZone++) { + + /*--- In principle, the mesh does not need to be updated ---*/ + UpdateMesh = 0; + + /*--- Transfer from all the remaining zones ---*/ + for (jZone = 0; jZone < nZone; jZone++){ + /*--- The target zone is iZone ---*/ + if (jZone != iZone && interface_container[iZone][jZone] != NULL){ + DeformMesh = Transfer_Data(jZone, iZone); + if (DeformMesh) UpdateMesh+=1; + } + } + /*--- If a mesh update is required due to the transfer of data ---*/ + if (UpdateMesh > 0) DynamicMeshUpdate(iZone, ExtIter); + } + + AD::Push_TapePosition(); + + for(iZone = 0; iZone < nZone; iZone++) { + + AD::Push_TapePosition(); + + if (tape_type == Kind_Tape::ZONE_SPECIFIC_TAPE) { + if (iZone == record_zone) { + DirectIteration(iZone, kind_recording); + } + } + else { + DirectIteration(iZone, kind_recording); + } + + iteration_container[iZone][INST_0]->RegisterOutput(solver_container, geometry_container, + config_container, output_container[iZone], iZone, INST_0); + + AD::Push_TapePosition(); + } + } + + if (rank == MASTER_NODE && kind_recording == FLOW_CONS_VARS) { + + if(config_container[record_zone]->GetWrt_AD_Statistics()) { + + AD::PrintStatistics(); + } + + cout << "-------------------------------------------------------------------------" << endl << endl; + } + + AD::StopRecording(); + + RecordingState = kind_recording; +} + +void CDiscAdjMultizoneDriver::DirectIteration(unsigned short iZone, unsigned short kind_recording) { + + /*--- Do one iteration of the direct solver ---*/ + direct_iteration[iZone][INST_0]->Preprocess(output_container[iZone], integration_container, geometry_container, solver_container, + numerics_container, config_container, surface_movement, grid_movement, FFDBox, iZone, INST_0); + + /*--- Iterate the zone as a block a single time ---*/ + direct_iteration[iZone][INST_0]->Iterate(output_container[iZone], integration_container, geometry_container, solver_container, + numerics_container, config_container, surface_movement, grid_movement, FFDBox, iZone, INST_0); + + Corrector(iZone); + + /*--- Print residuals in the first iteration ---*/ + + if (rank == MASTER_NODE && kind_recording == FLOW_CONS_VARS) { + + switch (config_container[iZone]->GetKind_Solver()) { + + case DISC_ADJ_EULER: case DISC_ADJ_NAVIER_STOKES: + case DISC_ADJ_INC_EULER: case DISC_ADJ_INC_NAVIER_STOKES: + cout << " Zone " << iZone << " (flow) - log10[RMS Solution_0]: " + << log10(solver_container[iZone][INST_0][MESH_0][FLOW_SOL]->GetRes_RMS(0)) << endl; + break; + case DISC_ADJ_RANS: case DISC_ADJ_INC_RANS: + cout << " Zone " << iZone << " (flow) - log10[RMS Solution_0]: " + << log10(solver_container[iZone][INST_0][MESH_0][FLOW_SOL]->GetRes_RMS(0)) << endl; + + if (!config_container[iZone]->GetFrozen_Visc_Disc()) { + + cout << " Zone " << iZone << " (turb) - log10[RMS k] : " + << log10(solver_container[iZone][INST_0][MESH_0][TURB_SOL]->GetRes_RMS(0)) << endl; + } + case DISC_ADJ_HEAT: + cout << " Zone " << iZone << " (heat) - log10[RMS Solution_0]: " + << log10(solver_container[iZone][INST_0][MESH_0][HEAT_SOL]->GetRes_RMS(0)) << endl; + break; + } + } +} + +void CDiscAdjMultizoneDriver::SetObjFunction(unsigned short kind_recording) { + + ObjFunc = 0.0; + su2double Weight_ObjFunc; + + unsigned short iZone, iMarker_Analyze, nMarker_Analyze; + + + /*--- Call objective function calculations. ---*/ + + for (iZone = 0; iZone < nZone; iZone++){ + + switch (config_container[iZone]->GetKind_Solver()) { + + case DISC_ADJ_EULER: case DISC_ADJ_NAVIER_STOKES: case DISC_ADJ_RANS: + case DISC_ADJ_INC_EULER: case DISC_ADJ_INC_NAVIER_STOKES: case DISC_ADJ_INC_RANS: + solver_container[iZone][INST_0][MESH_0][FLOW_SOL]->Pressure_Forces(geometry_container[iZone][INST_0][MESH_0], + config_container[iZone]); + solver_container[iZone][INST_0][MESH_0][FLOW_SOL]->Momentum_Forces(geometry_container[iZone][INST_0][MESH_0], + config_container[iZone]); + solver_container[iZone][INST_0][MESH_0][FLOW_SOL]->Friction_Forces(geometry_container[iZone][INST_0][MESH_0], + config_container[iZone]); + + if(config_container[iZone]->GetWeakly_Coupled_Heat()) { + + solver_container[iZone][INST_0][MESH_0][HEAT_SOL]->Heat_Fluxes(geometry_container[iZone][INST_0][MESH_0], + solver_container[iZone][INST_0][MESH_0], config_container[iZone]); + } + solver_container[iZone][INST_0][MESH_0][FLOW_SOL]->Evaluate_ObjFunc(config_container[iZone]); + break; + case DISC_ADJ_HEAT: + solver_container[iZone][INST_0][MESH_0][HEAT_SOL]->Heat_Fluxes(geometry_container[iZone][INST_0][MESH_0], + solver_container[iZone][INST_0][MESH_0], config_container[iZone]); + break; + } + + direct_output[iZone]->SetHistory_Output(geometry_container[iZone][INST_0][MESH_0], + solver_container[iZone][INST_0][MESH_0], config_container[iZone]); + } + + /*--- Extract objective function values. ---*/ + + for (iZone = 0; iZone < nZone; iZone++){ + + nMarker_Analyze = config_container[iZone]->GetnMarker_Analyze(); + + for (iMarker_Analyze = 0; iMarker_Analyze < nMarker_Analyze; iMarker_Analyze++) { + + Weight_ObjFunc = config_container[iZone]->GetWeight_ObjFunc(iMarker_Analyze); + + switch (config_container[iZone]->GetKind_Solver()) { + + case DISC_ADJ_EULER: case DISC_ADJ_NAVIER_STOKES: case DISC_ADJ_RANS: + // per-surface output to be added soon + break; + case HEAT_EQUATION_FVM: case DISC_ADJ_HEAT: + // per-surface output to be added soon + break; + default: + break; + } + } + + /*--- Not-per-surface objective functions (shall not be included above) ---*/ + + Weight_ObjFunc = config_container[iZone]->GetWeight_ObjFunc(0); + + switch (config_container[iZone]->GetKind_Solver()) { + + case DISC_ADJ_EULER: case DISC_ADJ_NAVIER_STOKES: case DISC_ADJ_RANS: + case DISC_ADJ_INC_EULER: case DISC_ADJ_INC_NAVIER_STOKES: case DISC_ADJ_INC_RANS: + + switch(config_container[iZone]->GetKind_ObjFunc()) { + + // Aerodynamic coefficients + + case DRAG_COEFFICIENT: + ObjFunc+=direct_output[iZone]->GetHistoryFieldValue("DRAG")*Weight_ObjFunc; + break; + case LIFT_COEFFICIENT: + ObjFunc+=direct_output[iZone]->GetHistoryFieldValue("LIFT")*Weight_ObjFunc; + break; + case SIDEFORCE_COEFFICIENT: + ObjFunc+=direct_output[iZone]->GetHistoryFieldValue("SIDEFORCE")*Weight_ObjFunc; + break; + case EFFICIENCY: + ObjFunc+=direct_output[iZone]->GetHistoryFieldValue("EFFICIENCY")*Weight_ObjFunc; + break; + case MOMENT_X_COEFFICIENT: + ObjFunc+=direct_output[iZone]->GetHistoryFieldValue("MOMENT-X")*Weight_ObjFunc; + break; + case MOMENT_Y_COEFFICIENT: + ObjFunc+=direct_output[iZone]->GetHistoryFieldValue("MOMENT-Y")*Weight_ObjFunc; + break; + case MOMENT_Z_COEFFICIENT: + ObjFunc+=direct_output[iZone]->GetHistoryFieldValue("MOMENT-Z")*Weight_ObjFunc; + break; + case FORCE_X_COEFFICIENT: + ObjFunc+=direct_output[iZone]->GetHistoryFieldValue("FORCE-X")*Weight_ObjFunc; + break; + case FORCE_Y_COEFFICIENT: + ObjFunc+=direct_output[iZone]->GetHistoryFieldValue("FORCE-Y")*Weight_ObjFunc; + break; + case FORCE_Z_COEFFICIENT: + ObjFunc+=direct_output[iZone]->GetHistoryFieldValue("FORCE-Z")*Weight_ObjFunc; + break; + + // Other surface-related output values + + case SURFACE_MASSFLOW: + ObjFunc+=direct_output[iZone]->GetHistoryFieldValue("AVG_MASSFLOW")*Weight_ObjFunc; + break; + case SURFACE_MACH: + ObjFunc+=direct_output[iZone]->GetHistoryFieldValue("AVG_MACH")*Weight_ObjFunc; + break; + case SURFACE_UNIFORMITY: + ObjFunc+=direct_output[iZone]->GetHistoryFieldValue("UNIFORMITY")*Weight_ObjFunc; + break; + case SURFACE_SECONDARY: + ObjFunc+=direct_output[iZone]->GetHistoryFieldValue("SECONDARY_STRENGTH")*Weight_ObjFunc; + break; + case SURFACE_MOM_DISTORTION: + ObjFunc+=direct_output[iZone]->GetHistoryFieldValue("MOMENTUM_DISTORTION")*Weight_ObjFunc; + break; + case SURFACE_SECOND_OVER_UNIFORM: + ObjFunc+=direct_output[iZone]->GetHistoryFieldValue("SECONDARY_OVER_UNIFORMITY")*Weight_ObjFunc; + break; + case TOTAL_AVG_TEMPERATURE: + ObjFunc+=direct_output[iZone]->GetHistoryFieldValue("AVG_TOTALTEMP")*Weight_ObjFunc; + break; + case SURFACE_TOTAL_PRESSURE: + ObjFunc+=direct_output[iZone]->GetHistoryFieldValue("AVG_TOTALPRESS")*Weight_ObjFunc; + break; + + // Not yet covered by new output structure. Be careful these use MARKER_MONITORING. + + case SURFACE_PRESSURE_DROP: + ObjFunc+=config_container[iZone]->GetSurface_PressureDrop(0)*Weight_ObjFunc; + break; + case SURFACE_STATIC_PRESSURE: + ObjFunc+=config_container[iZone]->GetSurface_Pressure(0)*Weight_ObjFunc; + break; + case TOTAL_HEATFLUX: + ObjFunc += solver_container[iZone][INST_0][MESH_0][FLOW_SOL]->GetTotal_HeatFlux()*Weight_ObjFunc; + break; + + default: + cout << "Objective function not covered for discrete adjoint multiphysics." << endl; + break; + } + break; + + case DISC_ADJ_HEAT: + + switch(config_container[iZone]->GetKind_ObjFunc()) { + + // Not yet covered by new output structure. Be careful these use MARKER_MONITORING. + + case TOTAL_HEATFLUX: + ObjFunc += solver_container[iZone][INST_0][MESH_0][HEAT_SOL]->GetTotal_HeatFlux()*Weight_ObjFunc; + break; + case TOTAL_AVG_TEMPERATURE: + ObjFunc += solver_container[iZone][INST_0][MESH_0][HEAT_SOL]->GetTotal_AvgTemperature()*Weight_ObjFunc; + break; + + default: + cout << "Objective function not covered for discrete adjoint multiphysics." << endl; + break; + } + break; + + default: + break; + } + } + + if (rank == MASTER_NODE) { + AD::RegisterOutput(ObjFunc); + AD::SetAdjIndex(ObjFunc_Index, ObjFunc); + if (rank == MASTER_NODE && kind_recording == FLOW_CONS_VARS) { + + cout << " Objective function : " << ObjFunc << " (" << ObjFunc_Index << ")" << endl; + } + } +} + +void CDiscAdjMultizoneDriver::SetAdj_ObjFunction() { + + bool TimeDomain = config_container[ZONE_0]->GetTime_Marching() != STEADY; + unsigned long IterAvg_Obj = config_container[ZONE_0]->GetIter_Avg_Objective(); + + su2double seeding = 1.0; + + if (TimeDomain){ + if (TimeIter < IterAvg_Obj){ + seeding = 1.0/((su2double)IterAvg_Obj); + } + else{ + seeding = 0.0; + } + } + + if (rank == MASTER_NODE) { + AD::SetDerivative(ObjFunc_Index, SU2_TYPE::GetValue(seeding)); + } +} + +void CDiscAdjMultizoneDriver::ComputeAdjoints(unsigned short iZone) { + + unsigned short enter_izone = iZone*2+1 + ITERATION_READY; + unsigned short leave_izone = iZone*2 + ITERATION_READY; + + AD::ClearAdjoints(); + + /*--- Initialize the adjoints in iZone ---*/ + + iteration_container[iZone][INST_0]->InitializeAdjoint(solver_container, geometry_container, + config_container, iZone, INST_0); + + /*--- Interpret the stored information by calling the corresponding routine of the AD tool. ---*/ + + AD::ComputeAdjoint(enter_izone, leave_izone); + + AD::ComputeAdjoint(TRANSFER, OBJECTIVE_FUNCTION); + AD::ComputeAdjoint(DEPENDENCIES, START); +} + +void CDiscAdjMultizoneDriver::Add_ExternalOld_To_Solution(unsigned short iZone) { + + unsigned short iSol; + + for (iSol=0; iSol < MAX_SOLS; iSol++){ + if (solver_container[iZone][INST_0][MESH_0][iSol] != NULL) { + if (solver_container[iZone][INST_0][MESH_0][iSol]->GetAdjoint()) { + solver_container[iZone][INST_0][MESH_0][iSol]->Add_ExternalOld_To_Solution(geometry_container[iZone][INST_0][MESH_0]); + } + } + } +} + +void CDiscAdjMultizoneDriver::SetExternal_Zero(void) { + + unsigned short iZone, iSol; + + for(iZone=0; iZone < nZone; iZone++) { + + for (iSol=0; iSol < MAX_SOLS; iSol++){ + if (solver_container[iZone][INST_0][MESH_0][iSol] != NULL) { + if (solver_container[iZone][INST_0][MESH_0][iSol]->GetAdjoint()) { + solver_container[iZone][INST_0][MESH_0][iSol]->SetExternal_Zero(geometry_container[iZone][INST_0][MESH_0]); + } + } + } + } +} + +void CDiscAdjMultizoneDriver::Set_OldExternal(void) { + + unsigned short iZone, iSol; + + for(iZone=0; iZone < nZone; iZone++) { + + for (iSol=0; iSol < MAX_SOLS; iSol++){ + if (solver_container[iZone][INST_0][MESH_0][iSol] != NULL) { + if (solver_container[iZone][INST_0][MESH_0][iSol]->GetAdjoint()) { + solver_container[iZone][INST_0][MESH_0][iSol]->Set_OldExternal(geometry_container[iZone][INST_0][MESH_0]); + } + } + } + } +} + +void CDiscAdjMultizoneDriver::Add_Solution_To_External(unsigned short iZone) { + + unsigned short iSol; + + for (iSol=0; iSol < MAX_SOLS; iSol++){ + if (solver_container[iZone][INST_0][MESH_0][iSol] != NULL) { + if (solver_container[iZone][INST_0][MESH_0][iSol]->GetAdjoint()) { + solver_container[iZone][INST_0][MESH_0][iSol]->Add_Solution_To_External(geometry_container[iZone][INST_0][MESH_0]); + } + } + } +} + +void CDiscAdjMultizoneDriver::Add_Solution_To_ExternalOld(unsigned short iZone) { + + unsigned short iSol; + + for (iSol=0; iSol < MAX_SOLS; iSol++){ + if (solver_container[iZone][INST_0][MESH_0][iSol] != NULL) { + if (solver_container[iZone][INST_0][MESH_0][iSol]->GetAdjoint()) { + solver_container[iZone][INST_0][MESH_0][iSol]->Add_Solution_To_ExternalOld(geometry_container[iZone][INST_0][MESH_0]); + } + } + } +} + +void CDiscAdjMultizoneDriver::Set_BGSSolution(unsigned short iZone) { + + unsigned short iSol; + + for (iSol=0; iSol < MAX_SOLS; iSol++){ + if (solver_container[iZone][INST_0][MESH_0][iSol] != NULL) { + if (solver_container[iZone][INST_0][MESH_0][iSol]->GetAdjoint()) { + solver_container[iZone][INST_0][MESH_0][iSol]->UpdateSolution_BGS(geometry_container[iZone][INST_0][MESH_0], + config_container[iZone]); + } + } + } +} + +void CDiscAdjMultizoneDriver::SetResidual_BGS(unsigned short iZone) { + + unsigned short iSol; + + for (iSol=0; iSol < MAX_SOLS; iSol++){ + if (solver_container[iZone][INST_0][MESH_0][iSol] != NULL) { + if (solver_container[iZone][INST_0][MESH_0][iSol]->GetAdjoint()) { + solver_container[iZone][INST_0][MESH_0][iSol]->ComputeResidual_Multizone(geometry_container[iZone][INST_0][MESH_0], + config_container[iZone]); + } + } + } +} diff --git a/SU2_CFD/src/drivers/CDriver.cpp b/SU2_CFD/src/drivers/CDriver.cpp index 24118983c5be..5db45c5b3436 100644 --- a/SU2_CFD/src/drivers/CDriver.cpp +++ b/SU2_CFD/src/drivers/CDriver.cpp @@ -570,7 +570,7 @@ void CDriver::Input_Preprocessing(CConfig **&config, CConfig *&driver_config) { /*--- Initialize the configuration of the driver ---*/ - driver_config = new CConfig(config_file_name, SU2_CFD, nZone, false); + driver_config = new CConfig(config_file_name, SU2_CFD, false); for (iZone = 0; iZone < nZone; iZone++) { @@ -1926,7 +1926,7 @@ void CDriver::Numerics_Preprocessing(CConfig *config, CGeometry **geometry, CSol case FEM_NAVIER_STOKES: case DISC_ADJ_FEM_NS : compressible =true; fem_ns = true; break; case FEM_RANS : case DISC_ADJ_FEM_RANS : compressible =true; fem_ns = true; fem_turbulent = true; break; case FEM_LES : compressible =true; fem_ns = true; break; - case HEAT_EQUATION_FVM: heat_fvm = true; break; + case HEAT_EQUATION_FVM: case DISC_ADJ_HEAT: heat_fvm = true; break; case FEM_ELASTICITY: case DISC_ADJ_FEM: fem = true; break; case ADJ_EULER : compressible =true; euler = true; adj_euler = true; break; case ADJ_NAVIER_STOKES : compressible =true; ns = true; turbulent = (config->GetKind_Turb_Model() != NONE); adj_ns = true; break; diff --git a/SU2_CFD/src/interfaces/cht/CConjugateHeatInterface.cpp b/SU2_CFD/src/interfaces/cht/CConjugateHeatInterface.cpp index bddd32cb5756..2f0a0e16a88e 100644 --- a/SU2_CFD/src/interfaces/cht/CConjugateHeatInterface.cpp +++ b/SU2_CFD/src/interfaces/cht/CConjugateHeatInterface.cpp @@ -38,39 +38,31 @@ #include "../../../include/interfaces/cht/CConjugateHeatInterface.hpp" -CConjugateHeatInterface::CConjugateHeatInterface(void) : CInterface() { - -} +CConjugateHeatInterface::CConjugateHeatInterface(void) : CInterface() { } CConjugateHeatInterface::CConjugateHeatInterface(unsigned short val_nVar, unsigned short val_nConst, - CConfig *config) : CInterface(val_nVar, val_nConst, config) { - -} - -CConjugateHeatInterface::~CConjugateHeatInterface(void) { + CConfig *config) : CInterface(val_nVar, val_nConst, config) { } -} +CConjugateHeatInterface::~CConjugateHeatInterface(void) { } void CConjugateHeatInterface::GetDonor_Variable(CSolver *donor_solution, CGeometry *donor_geometry, CConfig *donor_config, unsigned long Marker_Donor, unsigned long Vertex_Donor, unsigned long Point_Donor) { - unsigned long iPoint; - unsigned long PointNormal; unsigned short nDim, iDim; + unsigned long iPoint, PointNormal; + + su2double *Coord, *Coord_Normal, *Normal, *Edge_Vector, dist, dist2, Area, + Twall, Tnormal, dTdn, rho_cp_solid, Prandtl_Lam, laminar_viscosity, + thermal_diffusivity, thermal_conductivity, thermal_conductivityND, + heat_flux_density, conductivity_over_dist; nDim = donor_geometry->GetnDim(); - su2double *Coord, *Coord_Normal, *Normal, *Edge_Vector, dist, dist2, Area, Twall = 0.0, Tnormal = 0.0, - dTdn, cp_fluid, rho_cp_solid, Prandtl_Lam, laminar_viscosity, - thermal_diffusivity, thermal_conductivity, thermal_conductivityND, - heat_flux_density, conductivity_over_dist, Temperature_Ref; - su2double Gamma = donor_config->GetGamma(); - su2double Gas_Constant = donor_config->GetGas_ConstantND(); - su2double Cp = (Gamma / (Gamma - 1.0)) * Gas_Constant; Edge_Vector = new su2double[nDim]; /*--- Check whether the current zone is a solid zone or a fluid zone ---*/ + bool flow = ((donor_config->GetKind_Solver() == NAVIER_STOKES) || (donor_config->GetKind_Solver() == RANS) || (donor_config->GetKind_Solver() == DISC_ADJ_NAVIER_STOKES) @@ -79,23 +71,20 @@ void CConjugateHeatInterface::GetDonor_Variable(CSolver *donor_solution, CGeomet || (donor_config->GetKind_Solver() == INC_RANS) || (donor_config->GetKind_Solver() == DISC_ADJ_INC_NAVIER_STOKES) || (donor_config->GetKind_Solver() == DISC_ADJ_INC_RANS)); - bool compressible_flow = (donor_config->GetKind_Regime() == COMPRESSIBLE) && flow; - bool incompressible_flow = (donor_config->GetEnergy_Equation()) && flow; - bool heat_equation = donor_config->GetKind_Solver() == HEAT_EQUATION_FVM; - Temperature_Ref = donor_config->GetTemperature_Ref(); - Prandtl_Lam = donor_config->GetPrandtl_Lam(); - laminar_viscosity = donor_config->GetMu_ConstantND(); // TDE check for consistency - cp_fluid = donor_config->GetSpecific_Heat_Cp(); - rho_cp_solid = donor_config->GetSpecific_Heat_Cp_Solid()*donor_config->GetDensity_Solid(); + bool compressible_flow = (donor_config->GetKind_Regime() == COMPRESSIBLE) && flow; + bool incompressible_flow = (donor_config->GetEnergy_Equation()) && flow; + bool heat_equation = (donor_config->GetKind_Solver() == HEAT_EQUATION_FVM + || donor_config->GetKind_Solver() == DISC_ADJ_HEAT); - PointNormal = donor_geometry->vertex[Marker_Donor][Vertex_Donor]->GetNormal_Neighbor(); Coord = donor_geometry->node[Point_Donor]->GetCoord(); - Coord_Normal = donor_geometry->node[PointNormal]->GetCoord(); + Normal = donor_geometry->vertex[Marker_Donor][Vertex_Donor]->GetNormal(); + PointNormal = donor_geometry->vertex[Marker_Donor][Vertex_Donor]->GetNormal_Neighbor(); + Coord_Normal = donor_geometry->node[PointNormal]->GetCoord(); + + Twall = 0.0; Tnormal = 0.0; dTdn = 0.0; dist2 = 0.0; Area = 0.0; - dist2 = 0.0; - Area = 0.0; for (iDim = 0; iDim < nDim; iDim++) { Edge_Vector[iDim] = Coord_Normal[iDim] - Coord[iDim]; dist2 += Edge_Vector[iDim]*Edge_Vector[iDim]; @@ -104,28 +93,28 @@ void CConjugateHeatInterface::GetDonor_Variable(CSolver *donor_solution, CGeomet dist = sqrt(dist2); Area = sqrt(Area); - /*--- Retrieve temperature solution (later set is as the first donor variable) and its gradient ---*/ - - dTdn = 0.0; + /*--- Retrieve temperature solution and its gradient ---*/ if (compressible_flow) { - Twall = donor_solution->node[Point_Donor]->GetPrimitive(0)*Temperature_Ref; - Tnormal = donor_solution->node[PointNormal]->GetPrimitive(0)*Temperature_Ref; + Twall = donor_solution->node[Point_Donor]->GetPrimitive(0); + Tnormal = donor_solution->node[PointNormal]->GetPrimitive(0); dTdn = (Twall - Tnormal)/dist; } else if (incompressible_flow) { - Twall = donor_solution->node[Point_Donor]->GetTemperature()*Temperature_Ref; - Tnormal = donor_solution->node[PointNormal]->GetTemperature()*Temperature_Ref; + Twall = donor_solution->node[Point_Donor]->GetTemperature(); + Tnormal = donor_solution->node[PointNormal]->GetTemperature(); dTdn = (Twall - Tnormal)/dist; } - else if (flow || heat_equation) { - Twall = donor_solution->node[Point_Donor]->GetSolution(0)*Temperature_Ref; - Tnormal = donor_solution->node[PointNormal]->GetSolution(0)*Temperature_Ref; + else if (heat_equation) { + + Twall = donor_solution->node[Point_Donor]->GetSolution(0); + Tnormal = donor_solution->node[PointNormal]->GetSolution(0); + // TODO: Check if these improve accuracy, if needed at all // for (iDim = 0; iDim < nDim; iDim++) { // dTdn += (Twall - Tnormal)/dist * (Edge_Vector[iDim]/dist) * (Normal[iDim]/Area); // } @@ -133,66 +122,83 @@ void CConjugateHeatInterface::GetDonor_Variable(CSolver *donor_solution, CGeomet dTdn = (Twall - Tnormal)/dist; } else { - cout << "WARNING: Transfer of conjugate heat variables is called with non-supported donor solver!" << endl; + + SU2_MPI::Error("Transfer of conjugate heat variables failed (non-supported donor solver).", CURRENT_FUNCTION); } - /*--- Calculate the heat flux density (temperature gradient times thermal conductivity) - *--- and set it as second donor variable ---*/ + /*--- Calculate the heat flux density (temperature gradient times thermal conductivity) ---*/ + if (compressible_flow) { - iPoint = donor_geometry->vertex[Marker_Donor][Vertex_Donor]->GetNode(); + su2double Gamma = donor_config->GetGamma(); + su2double Gas_Constant = donor_config->GetGas_ConstantND(); + su2double Cp = (Gamma / (Gamma - 1.0)) * Gas_Constant; + + Prandtl_Lam = donor_config->GetPrandtl_Lam(); + laminar_viscosity = donor_solution->node[Point_Donor]->GetLaminarViscosity(); // TDE check for consistency + Cp = (Gamma / (Gamma - 1.0)) * Gas_Constant; thermal_conductivityND = Cp*(laminar_viscosity/Prandtl_Lam); - thermal_conductivity = thermal_conductivityND*donor_config->GetViscosity_Ref(); + heat_flux_density = thermal_conductivityND*dTdn; + + if (donor_config->GetCHT_Robin()) { - heat_flux_density = thermal_conductivity*dTdn; - conductivity_over_dist = thermal_conductivity/dist; + thermal_conductivity = thermal_conductivityND*donor_config->GetViscosity_Ref(); + conductivity_over_dist = thermal_conductivity/dist; + } } else if (incompressible_flow) { - iPoint = donor_geometry->vertex[Marker_Donor][Vertex_Donor]->GetNode(); - thermal_conductivityND = donor_solution->node[iPoint]->GetThermalConductivity(); - thermal_conductivity = thermal_conductivityND*donor_config->GetConductivity_Ref(); + heat_flux_density = thermal_conductivityND*dTdn; - switch (donor_config->GetKind_ConductivityModel()) { + if (donor_config->GetCHT_Robin()) { - case CONSTANT_CONDUCTIVITY: - thermal_conductivity = thermal_conductivityND*donor_config->GetConductivity_Ref(); - break; + switch (donor_config->GetKind_ConductivityModel()) { - case CONSTANT_PRANDTL: - thermal_conductivity = thermal_conductivityND*donor_config->GetGas_Constant_Ref() - *donor_config->GetViscosity_Ref(); - break; - } + case CONSTANT_CONDUCTIVITY: + thermal_conductivity = thermal_conductivityND*donor_config->GetConductivity_Ref(); + break; - heat_flux_density = thermal_conductivity*dTdn; - conductivity_over_dist = thermal_conductivity/dist; + case CONSTANT_PRANDTL: + thermal_conductivity = thermal_conductivityND*donor_config->GetGas_Constant_Ref() + *donor_config->GetViscosity_Ref(); + break; + } + + conductivity_over_dist = thermal_conductivity/dist; + } } - else if (flow) { + else if (heat_equation) { - iPoint = donor_geometry->vertex[Marker_Donor][Vertex_Donor]->GetNode(); + /*--- Heat solver stand-alone case ---*/ - thermal_conductivityND = laminar_viscosity/Prandtl_Lam; - thermal_conductivity = thermal_conductivityND*donor_config->GetViscosity_Ref()*cp_fluid; + thermal_diffusivity = donor_config->GetThermalDiffusivity_Solid(); + heat_flux_density = thermal_diffusivity*dTdn; - heat_flux_density = thermal_conductivity*dTdn; - conductivity_over_dist = thermal_conductivity/dist; + if (donor_config->GetCHT_Robin()) { + rho_cp_solid = donor_config->GetSpecific_Heat_Cp_Solid()*donor_config->GetDensity_Solid(); + conductivity_over_dist = thermal_diffusivity*rho_cp_solid/dist; + } } - else { - thermal_diffusivity = donor_config->GetThermalDiffusivity_Solid(); - heat_flux_density = (thermal_diffusivity*dTdn)*rho_cp_solid; - conductivity_over_dist = thermal_diffusivity*rho_cp_solid/dist; + /*--- Set the conjugate heat variables that are transferred by default ---*/ + + Donor_Variable[0] = Twall*donor_config->GetTemperature_Ref(); + Donor_Variable[1] = heat_flux_density*donor_config->GetHeat_Flux_Ref(); + + /*--- We only need these for the Robin BC option ---*/ + + if (donor_config->GetCHT_Robin()) { + + Donor_Variable[2] = conductivity_over_dist; + Donor_Variable[3] = Tnormal*donor_config->GetTemperature_Ref(); } - Donor_Variable[0] = Twall; - Donor_Variable[1] = heat_flux_density; - Donor_Variable[2] = conductivity_over_dist; - Donor_Variable[3] = Tnormal; + if (Edge_Vector != NULL) { - delete [] Edge_Vector; + delete [] Edge_Vector; + } } void CConjugateHeatInterface::SetTarget_Variable(CSolver *target_solution, CGeometry *target_geometry, @@ -203,8 +209,12 @@ void CConjugateHeatInterface::SetTarget_Variable(CSolver *target_solution, CGeom target_config->GetRelaxation_Factor_CHT(), Target_Variable[0]); target_solution->SetConjugateHeatVariable(Marker_Target, Vertex_Target, 1, target_config->GetRelaxation_Factor_CHT(), Target_Variable[1]); - target_solution->SetConjugateHeatVariable(Marker_Target, Vertex_Target, 2, - target_config->GetRelaxation_Factor_CHT(), Target_Variable[2]); - target_solution->SetConjugateHeatVariable(Marker_Target, Vertex_Target, 3, - target_config->GetRelaxation_Factor_CHT(), Target_Variable[3]); + + if (target_config->GetCHT_Robin()) { + + target_solution->SetConjugateHeatVariable(Marker_Target, Vertex_Target, 2, + target_config->GetRelaxation_Factor_CHT(), Target_Variable[2]); + target_solution->SetConjugateHeatVariable(Marker_Target, Vertex_Target, 3, + target_config->GetRelaxation_Factor_CHT(), Target_Variable[3]); + } } diff --git a/SU2_CFD/src/iteration_structure.cpp b/SU2_CFD/src/iteration_structure.cpp index 691f21dfc6f7..9ec7e4092ba1 100644 --- a/SU2_CFD/src/iteration_structure.cpp +++ b/SU2_CFD/src/iteration_structure.cpp @@ -2303,10 +2303,6 @@ void CDiscAdjFluidIteration::InitializeAdjoint(CSolver *****solver, CGeometry ** bool heat = config[iZone]->GetWeakly_Coupled_Heat(); bool interface_boundary = (config[iZone]->GetnMarker_Fluid_Load() > 0); - /*--- Initialize the adjoint of the objective function (typically with 1.0) ---*/ - - solver[iZone][iInst][MESH_0][ADJFLOW_SOL]->SetAdj_ObjFunc(geometry[iZone][iInst][MESH_0], config[iZone]); - /*--- Initialize the adjoints the conservative variables ---*/ if ((Kind_Solver == DISC_ADJ_NAVIER_STOKES) || (Kind_Solver == DISC_ADJ_RANS) || (Kind_Solver == DISC_ADJ_EULER) || @@ -2443,7 +2439,7 @@ void CDiscAdjFluidIteration::SetDependencies(CSolver *****solver, } /*--- Compute coupling between flow and turbulent equations ---*/ - + solver[iZone][iInst][MESH_0][FLOW_SOL]->Preprocessing(geometry[iZone][iInst][MESH_0], solver[iZone][iInst][MESH_0], config[iZone], MESH_0, NO_RK_ITER, RUNTIME_FLOW_SYS, true); solver[iZone][iInst][MESH_0][FLOW_SOL]->InitiateComms(geometry[iZone][iInst][MESH_0], config[iZone], SOLUTION); solver[iZone][iInst][MESH_0][FLOW_SOL]->CompleteComms(geometry[iZone][iInst][MESH_0], config[iZone], SOLUTION); @@ -2476,15 +2472,15 @@ void CDiscAdjFluidIteration::RegisterOutput(CSolver *****solver, CGeometry ****g /*--- Register conservative variables as output of the iteration ---*/ - solver[iZone][iInst][MESH_0][FLOW_SOL]->RegisterOutput(geometry[iZone][iInst][MESH_0],config[iZone]); + solver[iZone][iInst][MESH_0][ADJFLOW_SOL]->RegisterOutput(geometry[iZone][iInst][MESH_0],config[iZone]); } if (turbulent && !frozen_visc){ - solver[iZone][iInst][MESH_0][TURB_SOL]->RegisterOutput(geometry[iZone][iInst][MESH_0], + solver[iZone][iInst][MESH_0][ADJTURB_SOL]->RegisterOutput(geometry[iZone][iInst][MESH_0], config[iZone]); } if (heat){ - solver[iZone][iInst][MESH_0][HEAT_SOL]->RegisterOutput(geometry[iZone][iInst][MESH_0], + solver[iZone][iInst][MESH_0][ADJHEAT_SOL]->RegisterOutput(geometry[iZone][iInst][MESH_0], config[iZone]); } if (interface_boundary){ diff --git a/SU2_CFD/src/meson.build b/SU2_CFD/src/meson.build index 3195c085eff5..58e2c341bde1 100644 --- a/SU2_CFD/src/meson.build +++ b/SU2_CFD/src/meson.build @@ -112,6 +112,7 @@ su2_cfd_src += files(['interfaces/CInterface.cpp', su2_cfd_src += files(['drivers/CDriver.cpp', 'drivers/CMultizoneDriver.cpp', 'drivers/CSinglezoneDriver.cpp', + 'drivers/CDiscAdjMultizoneDriver.cpp', 'drivers/CDiscAdjSinglezoneDriver.cpp', 'drivers/CDummyDriver.cpp']) diff --git a/SU2_CFD/src/output/CAdjFlowIncOutput.cpp b/SU2_CFD/src/output/CAdjFlowIncOutput.cpp index 0dd430963119..395d011acc21 100644 --- a/SU2_CFD/src/output/CAdjFlowIncOutput.cpp +++ b/SU2_CFD/src/output/CAdjFlowIncOutput.cpp @@ -115,7 +115,7 @@ void CAdjFlowIncOutput::SetHistoryOutputFields(CConfig *config){ /// DESCRIPTION: Root-mean square residual of the adjoint Velocity z-component. AddHistoryOutput("RMS_ADJ_VELOCITY-Z", "rms[A_W]", ScreenOutputFormat::FIXED, "RMS_RES", "Root-mean square residual of the adjoint Velocity z-component.", HistoryFieldType::RESIDUAL); /// DESCRIPTION: Maximum residual of the temperature. - AddHistoryOutput("RMS_ADJ_HEAT", "rms[A_T]", ScreenOutputFormat::FIXED, "RMS_RES", " Root-mean square residual of the adjoint temperature.", HistoryFieldType::RESIDUAL); + AddHistoryOutput("RMS_ADJ_TEMPERATURE", "rms[A_T]", ScreenOutputFormat::FIXED, "RMS_RES", " Root-mean square residual of the adjoint temperature.", HistoryFieldType::RESIDUAL); if (!config->GetFrozen_Visc_Disc() || !config->GetFrozen_Visc_Cont()){ switch(turb_model){ case SA: case SA_NEG: case SA_E: case SA_COMP: case SA_E_COMP: @@ -143,7 +143,7 @@ void CAdjFlowIncOutput::SetHistoryOutputFields(CConfig *config){ /// DESCRIPTION: Maximum residual of the adjoint Velocity z-component AddHistoryOutput("MAX_ADJ_VELOCITY-Z", "max[A_RhoW]", ScreenOutputFormat::FIXED, "MAX_RES", "Maximum residual of the adjoint Velocity z-component", HistoryFieldType::RESIDUAL); /// DESCRIPTION: Maximum residual of the temperature. - AddHistoryOutput("MAX_ADJ_HEAT", "max[A_T]", ScreenOutputFormat::FIXED, "MAX_RES", "Maximum residual of the temperature.", HistoryFieldType::RESIDUAL); + AddHistoryOutput("MAX_ADJ_TEMPERATURE", "max[A_T]", ScreenOutputFormat::FIXED, "MAX_RES", "Maximum residual of the temperature.", HistoryFieldType::RESIDUAL); if (!config->GetFrozen_Visc_Disc() || !config->GetFrozen_Visc_Cont()){ switch(turb_model){ case SA: case SA_NEG: case SA_E: case SA_COMP: case SA_E_COMP: @@ -171,7 +171,7 @@ void CAdjFlowIncOutput::SetHistoryOutputFields(CConfig *config){ /// DESCRIPTION: BGS residual of the adjoint Velocity z-component AddHistoryOutput("BGS_ADJ_VELOCITY-Z", "bgs[A_RhoW]", ScreenOutputFormat::FIXED, "BGS_RES", "BGS residual of the adjoint Velocity z-component", HistoryFieldType::RESIDUAL); /// DESCRIPTION: BGS residual of the temperature. - AddHistoryOutput("BGS_ADJ_HEAT", "bgs[A_T]", ScreenOutputFormat::FIXED, "BGS_RES", "BGS residual of the adjoint temperature.", HistoryFieldType::RESIDUAL); + AddHistoryOutput("BGS_ADJ_TEMPERATURE", "bgs[A_T]", ScreenOutputFormat::FIXED, "BGS_RES", "BGS residual of the adjoint temperature.", HistoryFieldType::RESIDUAL); if (!config->GetFrozen_Visc_Disc() || !config->GetFrozen_Visc_Cont()){ switch(turb_model){ case SA: case SA_NEG: case SA_E: case SA_COMP: case SA_E_COMP: @@ -217,11 +217,11 @@ void CAdjFlowIncOutput::LoadHistoryData(CConfig *config, CGeometry *geometry, CS SetHistoryOutputValue("RMS_ADJ_VELOCITY-Z", log10(adjflow_solver->GetRes_RMS(3))); } if (weakly_coupled_heat){ - SetHistoryOutputValue("RMS_ADJ_HEAT", log10(adjheat_solver->GetRes_RMS(0))); + SetHistoryOutputValue("RMS_ADJ_TEMPERATURE", log10(adjheat_solver->GetRes_RMS(0))); } if (heat){ - if (nDim == 3) SetHistoryOutputValue("RMS_ADJ_HEAT", log10(adjflow_solver->GetRes_RMS(4))); - else SetHistoryOutputValue("RMS_ADJ_HEAT", log10(adjflow_solver->GetRes_RMS(3))); + if (nDim == 3) SetHistoryOutputValue("RMS_ADJ_TEMPERATURE", log10(adjflow_solver->GetRes_RMS(4))); + else SetHistoryOutputValue("RMS_ADJ_TEMPERATURE", log10(adjflow_solver->GetRes_RMS(3))); } if (!config->GetFrozen_Visc_Disc() || !config->GetFrozen_Visc_Cont()){ switch(turb_model){ @@ -242,11 +242,11 @@ void CAdjFlowIncOutput::LoadHistoryData(CConfig *config, CGeometry *geometry, CS SetHistoryOutputValue("MAX_ADJ_VELOCITY-Z", log10(adjflow_solver->GetRes_Max(3))); } if (weakly_coupled_heat){ - SetHistoryOutputValue("MAX_ADJ_HEAT", log10(adjheat_solver->GetRes_Max(0))); + SetHistoryOutputValue("MAX_ADJ_TEMPERATURE", log10(adjheat_solver->GetRes_Max(0))); } if (heat){ - if (nDim == 3) SetHistoryOutputValue("MAX_ADJ_HEAT", log10(adjflow_solver->GetRes_Max(4))); - else SetHistoryOutputValue("MAX_ADJ_HEAT", log10(adjflow_solver->GetRes_Max(3))); + if (nDim == 3) SetHistoryOutputValue("MAX_ADJ_TEMPERATURE", log10(adjflow_solver->GetRes_Max(4))); + else SetHistoryOutputValue("MAX_ADJ_TEMPERATURE", log10(adjflow_solver->GetRes_Max(3))); } if (!config->GetFrozen_Visc_Disc() || !config->GetFrozen_Visc_Cont()){ switch(turb_model){ @@ -269,11 +269,11 @@ void CAdjFlowIncOutput::LoadHistoryData(CConfig *config, CGeometry *geometry, CS SetHistoryOutputValue("BGS_ADJ_VELOCITY-Z", log10(adjflow_solver->GetRes_BGS(3))); } if (weakly_coupled_heat){ - SetHistoryOutputValue("BGS_ADJ_HEAT", log10(adjheat_solver->GetRes_BGS(0))); + SetHistoryOutputValue("BGS_ADJ_TEMPERATURE", log10(adjheat_solver->GetRes_BGS(0))); } if (heat){ - if (nDim == 3) SetHistoryOutputValue("BGS_ADJ_HEAT", log10(adjflow_solver->GetRes_BGS(4))); - else SetHistoryOutputValue("BGS_ADJ_HEAT", log10(adjflow_solver->GetRes_BGS(3))); + if (nDim == 3) SetHistoryOutputValue("BGS_ADJ_TEMPERATURE", log10(adjflow_solver->GetRes_BGS(4))); + else SetHistoryOutputValue("BGS_ADJ_TEMPERATURE", log10(adjflow_solver->GetRes_BGS(3))); } if (!config->GetFrozen_Visc_Disc() || !config->GetFrozen_Visc_Cont()){ switch(turb_model){ @@ -316,10 +316,9 @@ void CAdjFlowIncOutput::SetVolumeOutputFields(CConfig *config){ if (nDim == 3) /// DESCRIPTION: Adjoint Velocity z-component. AddVolumeOutput("ADJ_VELOCITY-Z", "Adjoint_Velocity_z", "SOLUTION", "z-component of the adjoint velocity vector"); - - if (weakly_coupled_heat || heat){ - AddVolumeOutput("ADJ_HEAT", "Adjoint_Heat", "SOLUTION", "Adjoint heat"); - } + + AddVolumeOutput("ADJ_TEMPERATURE", "Adjoint_Temperature", "SOLUTION", "Adjoint temperature"); + if (!config->GetFrozen_Visc_Disc()){ switch(turb_model){ @@ -355,9 +354,9 @@ void CAdjFlowIncOutput::SetVolumeOutputFields(CConfig *config){ AddVolumeOutput("RES_ADJ_VELOCITY-Y", "Residual_Adjoint_Velocity_y", "RESIDUAL", "Residual of the adjoint y-velocity"); if (nDim == 3) /// DESCRIPTION: Residual of the adjoint Velocity z-component. - AddVolumeOutput("RES_ADJ_Velocity-Z", "Residual_Adjoint_Velocity_z", "RESIDUAL", "Residual of the adjoint z-velocity"); + AddVolumeOutput("RES_ADJ_VELOCITY-Z", "Residual_Adjoint_Velocity_z", "RESIDUAL", "Residual of the adjoint z-velocity"); /// DESCRIPTION: Residual of the adjoint energy. - AddVolumeOutput("RES_ADJ_HEAT", "Residual_Adjoint_Heat", "RESIDUAL", "Residual of the adjoint heat"); + AddVolumeOutput("RES_ADJ_TEMPERATURE", "Residual_Adjoint_Heat", "RESIDUAL", "Residual of the adjoint temperature"); if (!config->GetFrozen_Visc_Disc()){ switch(turb_model){ case SA: case SA_NEG: case SA_E: case SA_COMP: case SA_E_COMP: @@ -416,11 +415,11 @@ void CAdjFlowIncOutput::LoadVolumeData(CConfig *config, CGeometry *geometry, CSo } if (weakly_coupled_heat){ - SetVolumeOutputValue("ADJ_HEAT", iPoint, Node_AdjHeat->GetSolution(0)); + SetVolumeOutputValue("ADJ_TEMPERATURE", iPoint, Node_AdjHeat->GetSolution(0)); } - if (heat){ - if (nDim == 3) SetVolumeOutputValue("ADJ_HEAT", iPoint, Node_AdjFlow->GetSolution(4)); - else SetVolumeOutputValue("ADJ_HEAT", iPoint, Node_AdjFlow->GetSolution(3)); + else { + if (nDim == 3) SetVolumeOutputValue("ADJ_TEMPERATURE", iPoint, Node_AdjFlow->GetSolution(4)); + else SetVolumeOutputValue("ADJ_TEMPERATURE", iPoint, Node_AdjFlow->GetSolution(3)); } // Turbulent if (!config->GetFrozen_Visc_Disc()){ @@ -444,6 +443,9 @@ void CAdjFlowIncOutput::LoadVolumeData(CConfig *config, CGeometry *geometry, CSo SetVolumeOutputValue("RES_ADJ_VELOCITY-Y", iPoint, Node_AdjFlow->GetSolution(2) - Node_AdjFlow->GetSolution_Old(2)); if (nDim == 3){ SetVolumeOutputValue("RES_ADJ_VELOCITY-Z", iPoint, Node_AdjFlow->GetSolution(3) - Node_AdjFlow->GetSolution_Old(3)); + SetVolumeOutputValue("RES_ADJ_TEMPERATURE", iPoint, Node_AdjFlow->GetSolution(4) - Node_AdjFlow->GetSolution_Old(4)); + } else { + SetVolumeOutputValue("RES_ADJ_TEMPERATURE", iPoint, Node_AdjFlow->GetSolution(3) - Node_AdjFlow->GetSolution_Old(3)); } if (!config->GetFrozen_Visc_Disc()){ switch(config->GetKind_Turb_Model()){ diff --git a/SU2_CFD/src/output/CAdjHeatOutput.cpp b/SU2_CFD/src/output/CAdjHeatOutput.cpp index feabf5c5adbd..0c59dff5e80a 100644 --- a/SU2_CFD/src/output/CAdjHeatOutput.cpp +++ b/SU2_CFD/src/output/CAdjHeatOutput.cpp @@ -107,7 +107,7 @@ void CAdjHeatOutput::SetHistoryOutputFields(CConfig *config){ /// BEGIN_GROUP: MAX_RES, DESCRIPTION: The maximum residuals of the conservative variables. /// DESCRIPTION: Maximum residual of the adjoint Pressure. - AddHistoryOutput("MAX_ADJ_TEMPERATURE", "max[A_T]", ScreenOutputFormat::FIXED, "BGS_RES", "BGS residual of the adjoint temperature.", HistoryFieldType::RESIDUAL); + 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. @@ -124,8 +124,9 @@ void CAdjHeatOutput::LoadHistoryData(CConfig *config, CGeometry *geometry, CSolv SetHistoryOutputValue("MAX_ADJ_TEMPERATURE", log10(adjheat_solver->GetRes_Max(0))); - if (multiZone) - SetHistoryOutputValue("MAX_ADJ_TEMPERATURE", log10(adjheat_solver->GetRes_BGS(0))); + if (multiZone) { + SetHistoryOutputValue("BGS_ADJ_TEMPERATURE", log10(adjheat_solver->GetRes_BGS(0))); + } SetHistoryOutputValue("SENS_GEO", adjheat_solver->GetTotal_Sens_Geo()); diff --git a/SU2_CFD/src/output/CFlowIncOutput.cpp b/SU2_CFD/src/output/CFlowIncOutput.cpp index 87633de48fae..21803e8ac2eb 100644 --- a/SU2_CFD/src/output/CFlowIncOutput.cpp +++ b/SU2_CFD/src/output/CFlowIncOutput.cpp @@ -278,7 +278,7 @@ void CFlowIncOutput::LoadHistoryData(CConfig *config, CGeometry *geometry, CSolv SetHistoryOutputValue("TEMPERATURE", heat_solver->GetTotal_AvgTemperature()); SetHistoryOutputValue("RMS_TEMPERATURE", log10(heat_solver->GetRes_RMS(0))); SetHistoryOutputValue("MAX_TEMPERATURE", log10(heat_solver->GetRes_Max(0))); - if (multiZone) SetHistoryOutputValue("BGS_HEAT", log10(heat_solver->GetRes_BGS(0))); + if (multiZone) SetHistoryOutputValue("BGS_TEMPERATURE", log10(heat_solver->GetRes_BGS(0))); } if (heat){ SetHistoryOutputValue("HEATFLUX", flow_solver->GetTotal_HeatFlux()); diff --git a/SU2_CFD/src/output/CHeatOutput.cpp b/SU2_CFD/src/output/CHeatOutput.cpp index 2e51d7387403..cb3dcbf15d9d 100644 --- a/SU2_CFD/src/output/CHeatOutput.cpp +++ b/SU2_CFD/src/output/CHeatOutput.cpp @@ -132,8 +132,11 @@ void CHeatOutput::SetVolumeOutputFields(CConfig *config){ // SOLUTION AddVolumeOutput("TEMPERATURE", "Temperature", "SOLUTION", "Temperature"); + // Primitives + AddVolumeOutput("HEAT_FLUX", "Heat_Flux", "PRIMITIVE", "Heatflux"); + // Residuals - AddVolumeOutput("RES_TEMPERATURE", "Residual_Temperature", "RESIDUAL", "RMS residual of the temperature"); + AddVolumeOutput("RES_TEMPERATURE", "Residual_Temperature", "RESIDUAL", "Residual of the temperature"); } @@ -158,3 +161,10 @@ void CHeatOutput::LoadVolumeData(CConfig *config, CGeometry *geometry, CSolver * } +void CHeatOutput::LoadSurfaceData(CConfig *config, CGeometry *geometry, CSolver **solver, unsigned long iPoint, unsigned short iMarker, unsigned long iVertex){ + + /* Heat flux value at each surface grid node. */ + SetVolumeOutputValue("HEAT_FLUX", iPoint, solver[HEAT_SOL]->GetHeatFlux(iMarker, iVertex)); + +} + diff --git a/SU2_CFD/src/output/COutput.cpp b/SU2_CFD/src/output/COutput.cpp index df756497af54..8614c0999d11 100644 --- a/SU2_CFD/src/output/COutput.cpp +++ b/SU2_CFD/src/output/COutput.cpp @@ -1550,6 +1550,18 @@ void COutput::Postprocess_HistoryData(CConfig *config){ Average[currentField.outputGroup].second++; } + + for (unsigned short iField = 0; iField < historyOutputPerSurface_List.size(); iField++){ + for (unsigned short iMarker = 0; iMarker < historyOutputPerSurface_Map[historyOutputPerSurface_List[iField]].size(); iMarker++){ + HistoryOutputField &Field = historyOutputPerSurface_Map[historyOutputPerSurface_List[iField]][iMarker]; + if (Field.fieldType == HistoryFieldType::COEFFICIENT){ + if (config->GetDirectDiff() != NO_DERIVATIVE){ + SetHistoryOutputValue("D_" + historyOutputPerSurface_List[iField][iMarker], SU2_TYPE::GetDerivative(Field.value)); + } + } + } + } + if (currentField.fieldType == HistoryFieldType::COEFFICIENT){ if(SetUpdate_Averages(config)){ if (config->GetTime_Domain()){ @@ -1636,6 +1648,15 @@ void COutput::Postprocess_HistoryFields(CConfig *config){ } } } + + for (unsigned short iField = 0; iField < historyOutputPerSurface_List.size(); iField++){ + for (unsigned short iMarker = 0; iMarker < historyOutputPerSurface_Map[historyOutputPerSurface_List[iField]].size(); iMarker++){ + HistoryOutputField &Field = historyOutputPerSurface_Map[historyOutputPerSurface_List[iField]][iMarker]; + if (Field.fieldType == HistoryFieldType::COEFFICIENT){ + AddHistoryOutput("D_" + historyOutputPerSurface_List[iField][iMarker], "d[" + Field.fieldName + "]", Field.screenFormat, "D_" + Field.outputGroup, "Derivative values for per-surface output.", HistoryFieldType::AUTO_COEFFICIENT); + } + } + } for (unsigned short iFieldConv = 0; iFieldConv < convFields.size(); iFieldConv++){ const string &convField = convFields[iFieldConv]; diff --git a/SU2_CFD/src/solver_adjoint_discrete.cpp b/SU2_CFD/src/solver_adjoint_discrete.cpp index 4cc42e00e070..ed6f368278ce 100644 --- a/SU2_CFD/src/solver_adjoint_discrete.cpp +++ b/SU2_CFD/src/solver_adjoint_discrete.cpp @@ -54,6 +54,10 @@ CDiscAdjSolver::CDiscAdjSolver(CGeometry *geometry, CConfig *config, CSolver *di ifstream restart_file; string filename, AdjExt; + adjoint = true; + + bool fsi = config->GetFSI_Simulation(); + nVar = direct_solver->GetnVar(); nDim = geometry->GetnDim(); @@ -136,6 +140,20 @@ CDiscAdjSolver::CDiscAdjSolver(CGeometry *geometry, CConfig *config, CSolver *di for (iPoint = 0; iPoint < nPoint; iPoint++) node[iPoint] = new CDiscAdjVariable(Solution, nDim, nVar, config); + switch(KindDirect_Solver){ + case RUNTIME_FLOW_SYS: + SolverName = "ADJ.FLOW"; + break; + case RUNTIME_HEAT_SYS: + SolverName = "ADJ.HEAT"; + break; + case RUNTIME_TURB_SYS: + SolverName = "ADJ.TURB"; + break; + default: + SolverName = "ADJ.SOL"; + break; + } } CDiscAdjSolver::~CDiscAdjSolver(void) { @@ -256,7 +274,13 @@ void CDiscAdjSolver::RegisterSolution(CGeometry *geometry, CConfig *config) { /*--- Register solution at all necessary time instances and other variables on the tape ---*/ for (iPoint = 0; iPoint < nPoint; iPoint++) { - direct_solver->node[iPoint]->RegisterSolution(input); + if(config->GetMultizone_Problem()) { + direct_solver->node[iPoint]->RegisterSolution_intIndexBased(input); + direct_solver->node[iPoint]->SetAdjIndices(input); + } + else { + direct_solver->node[iPoint]->RegisterSolution(input); + } } if (time_n_needed) { for (iPoint = 0; iPoint < nPoint; iPoint++) { @@ -376,7 +400,13 @@ void CDiscAdjSolver::RegisterOutput(CGeometry *geometry, CConfig *config) { /*--- Register output variables on the tape ---*/ for (iPoint = 0; iPoint < nPoint; iPoint++) { - direct_solver->node[iPoint]->RegisterSolution(input); + if(config->GetMultizone_Problem()) { + direct_solver->node[iPoint]->RegisterSolution_intIndexBased(input); + direct_solver->node[iPoint]->SetAdjIndices(input); + } + else { + direct_solver->node[iPoint]->RegisterSolution(input); + } } } @@ -480,11 +510,18 @@ void CDiscAdjSolver::ExtractAdjoint_Solution(CGeometry *geometry, CConfig *confi /*--- Set the old solution ---*/ - node[iPoint]->Set_OldSolution(); + if(!config->GetMultizone_Problem()) { + node[iPoint]->Set_OldSolution(); + } /*--- Extract the adjoint solution ---*/ - direct_solver->node[iPoint]->GetAdjointSolution(Solution); + if(config->GetMultizone_Problem()) { + direct_solver->node[iPoint]->GetAdjointSolution_intIndexBased(Solution); + } + else { + direct_solver->node[iPoint]->GetAdjointSolution(Solution); + } /*--- Store the adjoint solution ---*/ @@ -528,6 +565,7 @@ void CDiscAdjSolver::ExtractAdjoint_Solution(CGeometry *geometry, CConfig *confi } SetResidual_RMS(geometry, config); + } void CDiscAdjSolver::ExtractAdjoint_Variables(CGeometry *geometry, CConfig *config) { @@ -736,7 +774,12 @@ void CDiscAdjSolver::SetAdjoint_Output(CGeometry *geometry, CConfig *config) { for (iPoint = 0; iPoint < nPoint; iPoint++) { for (iVar = 0; iVar < nVar; iVar++) { - Solution[iVar] = node[iPoint]->GetSolution(iVar); + if(config->GetMultizone_Problem()) { + Solution[iVar] = node[iPoint]->Get_BGSSolution_k(iVar); + } + else { + Solution[iVar] = node[iPoint]->GetSolution(iVar); + } } if (fsi) { for (iVar = 0; iVar < nVar; iVar++) { @@ -748,9 +791,13 @@ void CDiscAdjSolver::SetAdjoint_Output(CGeometry *geometry, CConfig *config) { Solution[iVar] += node[iPoint]->GetDual_Time_Derivative(iVar); } } - direct_solver->node[iPoint]->SetAdjointSolution(Solution); + if(config->GetMultizone_Problem()) { + direct_solver->node[iPoint]->SetAdjointSolution_intIndexBased(Solution); + } + else { + direct_solver->node[iPoint]->SetAdjointSolution(Solution); + } } - } void CDiscAdjSolver::SetAdjoint_OutputMesh(CGeometry *geometry, CConfig *config){ @@ -801,11 +848,16 @@ void CDiscAdjSolver::SetSensitivity(CGeometry *geometry, CSolver **solver, CConf for (iDim = 0; iDim < nDim; iDim++) { - Sensitivity = SU2_TYPE::GetDerivative(Coord[iDim]); + if(config->GetMultizone_Problem()) { + Sensitivity = geometry->node[iPoint]->GetAdjointSolution(iDim); + } + else { + Sensitivity = SU2_TYPE::GetDerivative(Coord[iDim]); + } /*--- Set the index manually to zero. ---*/ - AD::ResetInput(Coord[iDim]); + AD::ResetInput(Coord[iDim]); /*--- If sharp edge, set the sensitivity to 0 on that region ---*/ @@ -842,7 +894,8 @@ void CDiscAdjSolver::SetSurface_Sensitivity(CGeometry *geometry, CConfig *config if(config->GetMarker_All_KindBC(iMarker) == EULER_WALL || config->GetMarker_All_KindBC(iMarker) == HEAT_FLUX - || config->GetMarker_All_KindBC(iMarker) == ISOTHERMAL) { + || config->GetMarker_All_KindBC(iMarker) == ISOTHERMAL + || config->GetMarker_All_KindBC(iMarker) == CHT_WALL_INTERFACE) { Sens = 0.0; @@ -1057,7 +1110,12 @@ void CDiscAdjSolver::ComputeResidual_Multizone(CGeometry *geometry, CConfig *con /*--- Compute the BGS solution (adding the cross term) ---*/ for (iPoint = 0; iPoint < nPointDomain; iPoint++){ for (iVar = 0; iVar < nVar; iVar++){ - bgs_sol = node[iPoint]->GetSolution(iVar) + node[iPoint]->GetCross_Term_Derivative(iVar); + if(config->GetMultizone_Problem()) { + bgs_sol = node[iPoint]->GetSolution(iVar); + } + else { + bgs_sol = node[iPoint]->GetSolution(iVar) + node[iPoint]->GetCross_Term_Derivative(iVar); + } node[iPoint]->Set_BGSSolution(iVar, bgs_sol); } } @@ -1076,4 +1134,3 @@ void CDiscAdjSolver::ComputeResidual_Multizone(CGeometry *geometry, CConfig *con SetResidual_BGS(geometry, config); } - diff --git a/SU2_CFD/src/solver_adjoint_elasticity.cpp b/SU2_CFD/src/solver_adjoint_elasticity.cpp index bbd221d050a4..4e12d66d8441 100644 --- a/SU2_CFD/src/solver_adjoint_elasticity.cpp +++ b/SU2_CFD/src/solver_adjoint_elasticity.cpp @@ -79,6 +79,8 @@ CDiscAdjFEASolver::CDiscAdjFEASolver(CGeometry *geometry, CConfig *config) : CS CDiscAdjFEASolver::CDiscAdjFEASolver(CGeometry *geometry, CConfig *config, CSolver *direct_solver, unsigned short Kind_Solver, unsigned short iMesh) : CSolver(){ + adjoint = true; + unsigned short iVar, iMarker, iDim; unsigned long iVertex, iPoint; diff --git a/SU2_CFD/src/solver_adjoint_mean.cpp b/SU2_CFD/src/solver_adjoint_mean.cpp index 7ae32d1296c4..3942c5e5377d 100644 --- a/SU2_CFD/src/solver_adjoint_mean.cpp +++ b/SU2_CFD/src/solver_adjoint_mean.cpp @@ -68,6 +68,8 @@ CAdjEulerSolver::CAdjEulerSolver(CGeometry *geometry, CConfig *config, unsigned string filename, AdjExt; su2double myArea_Monitored, Area, *Normal; + adjoint = true; + bool restart = config->GetRestart(); bool axisymmetric = config->GetAxisymmetric(); diff --git a/SU2_CFD/src/solver_adjoint_turbulent.cpp b/SU2_CFD/src/solver_adjoint_turbulent.cpp index 0e2cf8e5eaff..30357a6844d0 100644 --- a/SU2_CFD/src/solver_adjoint_turbulent.cpp +++ b/SU2_CFD/src/solver_adjoint_turbulent.cpp @@ -44,6 +44,8 @@ CAdjTurbSolver::CAdjTurbSolver(CGeometry *geometry, CConfig *config, unsigned sh unsigned long iPoint; unsigned short iDim, iVar, nLineLets; + adjoint = true; + nDim = geometry->GetnDim(); Gamma = config->GetGamma(); Gamma_Minus_One = Gamma - 1.0; diff --git a/SU2_CFD/src/solver_direct_heat.cpp b/SU2_CFD/src/solver_direct_heat.cpp index 4a43ca0321aa..e86bed9204fc 100644 --- a/SU2_CFD/src/solver_direct_heat.cpp +++ b/SU2_CFD/src/solver_direct_heat.cpp @@ -41,6 +41,7 @@ CHeatSolverFVM::CHeatSolverFVM(void) : CSolver() { ConjugateVar = NULL; + HeatFlux = NULL; } CHeatSolverFVM::CHeatSolverFVM(CGeometry *geometry, CConfig *config, unsigned short iMesh) : CSolver() { @@ -158,13 +159,13 @@ CHeatSolverFVM::CHeatSolverFVM(CGeometry *geometry, CConfig *config, unsigned sh Smatrix[iDim] = new su2double [nDim]; } - HeatFlux = new su2double*[nMarker]; - AvgTemperature = new su2double[nMarker]; - Surface_Areas = new su2double[config->GetnMarker_HeatFlux()]; + HeatFlux_per_Marker = new su2double[nMarker]; + AverageT_per_Marker = new su2double[nMarker]; + Surface_Areas = new su2double[config->GetnMarker_HeatFlux()]; for(iMarker = 0; iMarker < nMarker; iMarker++) { - HeatFlux[iMarker] = new su2double[geometry->GetnVertex(iMarker)]; - AvgTemperature[iMarker] = 0.0; + HeatFlux_per_Marker[iMarker] = 0.0; + AverageT_per_Marker[iMarker] = 0.0; } for(iMarker = 0; iMarker < config->GetnMarker_HeatFlux(); iMarker++) { Surface_Areas[iMarker] = 0.0; @@ -172,7 +173,7 @@ CHeatSolverFVM::CHeatSolverFVM(CGeometry *geometry, CConfig *config, unsigned sh Set_Heatflux_Areas(geometry, config); - /*--- Non-dimensionalization of heat equation */ + /*--- Set the reference values for temperature ---*/ su2double Temperature_FreeStream = config->GetInc_Temperature_Init(); config->SetTemperature_FreeStream(Temperature_FreeStream); @@ -190,17 +191,25 @@ CHeatSolverFVM::CHeatSolverFVM(CGeometry *geometry, CConfig *config, unsigned sh config->SetTemperature_Ref(Temperature_Ref); config->SetTemperature_FreeStreamND(config->GetTemperature_FreeStream()/config->GetTemperature_Ref()); - if (rank == MASTER_NODE) { - cout << "Weakly coupled heat solver's freestream temperature: " << config->GetTemperature_FreeStreamND() << endl; + + /*--- Set the reference values for heat fluxes. If the heat solver runs stand-alone, + * thermal conductivity is read directly from config file ---*/ + + if (heat_equation) { + + su2double rho_cp = config->GetDensity_Solid()*config->GetSpecific_Heat_Cp_Solid(); + + config->SetThermalDiffusivity_Solid(config->GetThermalConductivity_Solid() / rho_cp); + config->SetHeat_Flux_Ref(rho_cp*Temperature_Ref); } + else { - su2double Temperature_Solid_Freestream_ND = config->GetTemperature_Freestream_Solid()/config->GetTemperature_Ref(); - if (heat_equation && (rank == MASTER_NODE)) { - cout << "Heat solver freestream temperature in case for solids: " << Temperature_Solid_Freestream_ND << endl; + config->SetHeat_Flux_Ref(config->GetViscosity_Ref()*config->GetSpecific_Heat_Cp()); } /*--- Store the value of the temperature and the heat flux density at the boundaries, - used for IO with a donor cell ---*/ + used for communications with donor cells ---*/ + unsigned short nConjVariables = 4; ConjugateVar = new su2double** [nMarker]; @@ -216,11 +225,14 @@ CHeatSolverFVM::CHeatSolverFVM(CGeometry *geometry, CConfig *config, unsigned sh } } - /*--- If the heat solver runs stand-alone, we have to set the reference values ---*/ - if(heat_equation) { - su2double rho_cp = config->GetDensity_Solid()*config->GetSpecific_Heat_Cp_Solid(); - su2double thermal_diffusivity_solid = config->GetThermalConductivity_Solid() / rho_cp; - config->SetThermalDiffusivity_Solid(thermal_diffusivity_solid); + /*--- Heat flux in all the markers ---*/ + + HeatFlux = new su2double* [nMarker]; + for (iMarker = 0; iMarker < nMarker; iMarker++) { + HeatFlux[iMarker] = new su2double [geometry->nVertex[iMarker]]; + for (iVertex = 0; iVertex < geometry->nVertex[iMarker]; iVertex++) { + HeatFlux[iMarker][iVertex] = 0.0; + } } if (multizone){ @@ -238,11 +250,18 @@ CHeatSolverFVM::CHeatSolverFVM(CGeometry *geometry, CConfig *config, unsigned sh } } - for (iPoint = 0; iPoint < nPoint; iPoint++) - if (flow) - node[iPoint] = new CHeatFVMVariable(config->GetTemperature_FreeStreamND(), nDim, nVar, config); - else - node[iPoint] = new CHeatFVMVariable(Temperature_Solid_Freestream_ND, nDim, nVar, config); + /*--- Initialize solution array ---*/ + + for (iPoint = 0; iPoint < nPoint; iPoint++) { + + if (flow) { + node[iPoint] = new CHeatFVMVariable(config->GetTemperature_FreeStreamND(), nDim, nVar, config); + } + else { + su2double Temperature_Solid_Freestream_ND = config->GetTemperature_Freestream_Solid()/config->GetTemperature_Ref(); + node[iPoint] = new CHeatFVMVariable(Temperature_Solid_Freestream_ND, nDim, nVar, config); + } + } /*--- MPI solution ---*/ @@ -250,10 +269,22 @@ CHeatSolverFVM::CHeatSolverFVM(CGeometry *geometry, CConfig *config, unsigned sh CompleteComms(geometry, config, SOLUTION); /*--- Add the solver name (max 8 characters) ---*/ + SolverName = "HEAT"; } -CHeatSolverFVM::~CHeatSolverFVM(void) { } +CHeatSolverFVM::~CHeatSolverFVM(void) { + + unsigned short iMarker; + + if (HeatFlux != NULL) { + for (iMarker = 0; iMarker < nMarker; iMarker++) { + delete [] HeatFlux[iMarker]; + } + delete [] HeatFlux; + } + +} void CHeatSolverFVM::Preprocessing(CGeometry *geometry, CSolver **solver_container, CConfig *config, unsigned short iMesh, unsigned short iRKStep, unsigned short RunTime_EqSystem, bool Output) { @@ -268,10 +299,12 @@ void CHeatSolverFVM::Preprocessing(CGeometry *geometry, CSolver **solver_contain for (iPoint = 0; iPoint < nPoint; iPoint ++) { /*--- Initialize the residual vector ---*/ + LinSysRes.SetBlock_Zero(iPoint); } /*--- Initialize the Jacobian matrices ---*/ + Jacobian.SetValZero(); if (config->GetKind_Gradient_Method() == GREEN_GAUSS) SetSolution_Gradient_GG(geometry, config); @@ -283,6 +316,7 @@ void CHeatSolverFVM::Postprocessing(CGeometry *geometry, CSolver **solver_contai void CHeatSolverFVM::LoadRestart(CGeometry **geometry, CSolver ***solver, CConfig *config, int val_iter, bool val_update_geo) { /*--- Restart the solution from file information ---*/ + unsigned short iDim, iVar, iMesh; unsigned long iPoint, index, iChildren, Point_Fine; @@ -1089,8 +1123,8 @@ void CHeatSolverFVM::BC_ConjugateHeat_Interface(CGeometry *geometry, CSolver **s unsigned long iVertex, iPoint, total_index; unsigned short iDim, iVar, iMarker; - su2double Area, rho_cp_solid, - Temperature_Ref, Tinterface, T_Conjugate, Tnormal_Conjugate, Conductance, HeatFluxDensity, HeatFluxValue; + su2double thermal_diffusivity, rho_cp_solid, Temperature_Ref, T_Conjugate, Tinterface, + Tnormal_Conjugate, HeatFluxDensity, HeatFlux, Area; bool implicit = (config->GetKind_TimeIntScheme_Flow() == EULER_IMPLICIT); bool flow = ((config->GetKind_Solver() == INC_NAVIER_STOKES) @@ -1152,20 +1186,28 @@ void CHeatSolverFVM::BC_ConjugateHeat_Interface(CGeometry *geometry, CSolver **s for (iDim = 0; iDim < nDim; iDim++) Area += Normal[iDim]*Normal[iDim]; Area = sqrt (Area); - Tinterface = node[iPoint]->GetSolution(0); - Tnormal_Conjugate = GetConjugateHeatVariable(iMarker, iVertex, 3)/Temperature_Ref; - Conductance = GetConjugateHeatVariable(iMarker, iVertex, 2)/rho_cp_solid; + thermal_diffusivity = GetConjugateHeatVariable(iMarker, iVertex, 2)/rho_cp_solid; + + if (config->GetCHT_Robin()) { - HeatFluxDensity = Conductance*(Tinterface - Tnormal_Conjugate); + Tinterface = node[iPoint]->GetSolution(0); + Tnormal_Conjugate = GetConjugateHeatVariable(iMarker, iVertex, 3)/Temperature_Ref; - HeatFluxValue = HeatFluxDensity * Area; + HeatFluxDensity = thermal_diffusivity*(Tinterface - Tnormal_Conjugate); + HeatFlux = HeatFluxDensity * Area; + } + else { + + HeatFluxDensity = GetConjugateHeatVariable(iMarker, iVertex, 1)/config->GetHeat_Flux_Ref(); + HeatFlux = HeatFluxDensity*Area; + } - Res_Visc[0] = -HeatFluxValue; + Res_Visc[0] = -HeatFlux; LinSysRes.SubtractBlock(iPoint, Res_Visc); if (implicit) { - Jacobian_i[0][0] = Conductance*Area; + Jacobian_i[0][0] = thermal_diffusivity*Area; Jacobian.SubtractBlock(iPoint, iPoint, Jacobian_i); } } @@ -1179,26 +1221,25 @@ void CHeatSolverFVM::Heat_Fluxes(CGeometry *geometry, CSolver **solver_container unsigned long iVertex, iPoint, iPointNormal; unsigned short Boundary, Monitoring, iMarker, iDim; - su2double *Coord, *Coord_Normal, *Normal, Area, dist, Twall, dTdn, cp_fluid, rho_cp_solid, - thermal_conductivity, thermal_diffusivity; + su2double *Coord, *Coord_Normal, *Normal, Area, dist, Twall, dTdn, thermal_diffusivity; string Marker_Tag, HeatFlux_Tag; + bool flow = ((config->GetKind_Solver() == INC_NAVIER_STOKES) || (config->GetKind_Solver() == INC_RANS) || (config->GetKind_Solver() == DISC_ADJ_INC_NAVIER_STOKES) || (config->GetKind_Solver() == DISC_ADJ_INC_RANS)); + #ifdef HAVE_MPI - su2double MyAllBound_HeatFlux, MyAllBound_AvgTemperature; + su2double MyAllBound_HeatFlux, MyAllBound_AverageT; #endif - cp_fluid = config->GetSpecific_Heat_Cp(); - rho_cp_solid = config->GetSpecific_Heat_Cp_Solid()*config->GetDensity_Solid(); - AllBound_HeatFlux = 0.0; - AllBound_AvgTemperature = 0.0; + AllBound_AverageT = 0.0; for ( iMarker = 0; iMarker < nMarker; iMarker++ ) { - AvgTemperature[iMarker] = 0.0; + AverageT_per_Marker[iMarker] = 0.0; + HeatFlux_per_Marker[iMarker] = 0.0; Boundary = config->GetMarker_All_KindBC(iMarker); Marker_Tag = config->GetMarker_All_TagBound(iMarker); @@ -1232,17 +1273,15 @@ void CHeatSolverFVM::Heat_Fluxes(CGeometry *geometry, CSolver **solver_container if(flow) { thermal_diffusivity = config->GetViscosity_FreeStreamND()/config->GetPrandtl_Lam(); - thermal_conductivity = thermal_diffusivity*config->GetViscosity_Ref()*cp_fluid; } else { - thermal_conductivity = config->GetThermalDiffusivity_Solid()*rho_cp_solid; + thermal_diffusivity = config->GetThermalDiffusivity_Solid(); } - HeatFlux[iMarker][iVertex] = thermal_conductivity*dTdn*config->GetTemperature_Ref()*Area; + HeatFlux[iMarker][iVertex] = thermal_diffusivity*dTdn*config->GetHeat_Flux_Ref(); + + HeatFlux_per_Marker[iMarker] += HeatFlux[iMarker][iVertex]*Area; - if (Monitoring){ - AllBound_HeatFlux += HeatFlux[iMarker][iVertex]; - } } } } @@ -1274,47 +1313,46 @@ void CHeatSolverFVM::Heat_Fluxes(CGeometry *geometry, CSolver **solver_container if(flow) { thermal_diffusivity = config->GetViscosity_FreeStreamND()/config->GetPrandtl_Lam(); - thermal_conductivity = thermal_diffusivity*config->GetViscosity_Ref()*cp_fluid; } else { - thermal_conductivity = config->GetThermalDiffusivity_Solid()*rho_cp_solid; + thermal_diffusivity = config->GetThermalDiffusivity_Solid(); } - HeatFlux[iMarker][iVertex] = thermal_conductivity*dTdn*config->GetTemperature_Ref()*Area; - - if (Monitoring){ - AllBound_HeatFlux += HeatFlux[iMarker][iVertex]; - } + HeatFlux[iMarker][iVertex] = thermal_diffusivity*dTdn*config->GetHeat_Flux_Ref(); + HeatFlux_per_Marker[iMarker] += HeatFlux[iMarker][iVertex]*Area; + /*--- We do only aim to compute averaged temperatures on the (interesting) heat flux walls ---*/ + if ( Boundary == HEAT_FLUX ) { - AvgTemperature[iMarker] += Twall*config->GetTemperature_Ref()*Area; + AverageT_per_Marker[iMarker] += Twall*config->GetTemperature_Ref()*Area; } } } } if (Monitoring == YES) { - AllBound_AvgTemperature += AvgTemperature[iMarker]; + + AllBound_HeatFlux += HeatFlux_per_Marker[iMarker]; + AllBound_AverageT += AverageT_per_Marker[iMarker]; } } #ifdef HAVE_MPI MyAllBound_HeatFlux = AllBound_HeatFlux; - MyAllBound_AvgTemperature = AllBound_AvgTemperature; + MyAllBound_AverageT = AllBound_AverageT; SU2_MPI::Allreduce(&MyAllBound_HeatFlux, &AllBound_HeatFlux, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD); - SU2_MPI::Allreduce(&MyAllBound_AvgTemperature, &AllBound_AvgTemperature, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD); + SU2_MPI::Allreduce(&MyAllBound_AverageT, &AllBound_AverageT, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD); #endif if (Total_HeatFlux_Areas_Monitor != 0.0) { - Total_AvgTemperature = AllBound_AvgTemperature/Total_HeatFlux_Areas_Monitor; + Total_AverageT = AllBound_AverageT/Total_HeatFlux_Areas_Monitor; } else { - Total_AvgTemperature = 0.0; + Total_AverageT = 0.0; } - Total_HeatFlux = AllBound_HeatFlux; } diff --git a/SU2_CFD/src/solver_direct_mean_inc.cpp b/SU2_CFD/src/solver_direct_mean_inc.cpp index 735535228893..5c5c87fbc815 100644 --- a/SU2_CFD/src/solver_direct_mean_inc.cpp +++ b/SU2_CFD/src/solver_direct_mean_inc.cpp @@ -8485,12 +8485,9 @@ void CIncNSSolver::BC_Isothermal_Wall(CGeometry *geometry, CSolver **solver_cont void CIncNSSolver::BC_ConjugateHeat_Interface(CGeometry *geometry, CSolver **solver_container, CNumerics *conv_numerics, CConfig *config, unsigned short val_marker) { unsigned short iVar, jVar, iDim, Wall_Function; - unsigned long iVertex, iPoint, Point_Normal, total_index; + unsigned long iVertex, iPoint, total_index; - su2double *GridVel; - su2double *Normal, *Coord_i, *Coord_j, Area, dist_ij; - su2double Tconjugate, dTdn; - su2double thermal_conductivity; + su2double *GridVel, Tconjugate; su2double Temperature_Ref = config->GetTemperature_Ref(); bool implicit = (config->GetKind_TimeIntScheme_Flow() == EULER_IMPLICIT); @@ -8550,63 +8547,11 @@ void CIncNSSolver::BC_ConjugateHeat_Interface(CGeometry *geometry, CSolver **sol Tconjugate = GetConjugateHeatVariable(val_marker, iVertex, 0)/Temperature_Ref; -// node[iPoint]->SetSolution_Old(nDim+1, Tconjugate); -// node[iPoint]->SetEnergy_ResTruncError_Zero(); - - Normal = geometry->vertex[val_marker][iVertex]->GetNormal(); - - Area = 0.0; - for (iDim = 0; iDim < nDim; iDim++) - Area += Normal[iDim]*Normal[iDim]; - Area = sqrt (Area); - - /*--- Compute closest normal neighbor ---*/ - - Point_Normal = geometry->vertex[val_marker][iVertex]->GetNormal_Neighbor(); - - /*--- Get coordinates of i & nearest normal and compute distance ---*/ - - Coord_i = geometry->node[iPoint]->GetCoord(); - Coord_j = geometry->node[Point_Normal]->GetCoord(); - dist_ij = 0; - for (iDim = 0; iDim < nDim; iDim++) - dist_ij += (Coord_j[iDim]-Coord_i[iDim])*(Coord_j[iDim]-Coord_i[iDim]); - dist_ij = sqrt(dist_ij); - - /*--- Compute the normal gradient in temperature using Twall ---*/ - - dTdn = -(node[Point_Normal]->GetTemperature() - Tconjugate)/dist_ij; - - /*--- Get thermal conductivity ---*/ - - thermal_conductivity = node[iPoint]->GetThermalConductivity(); - - /*--- Apply a weak boundary condition for the energy equation. - Compute the residual due to the prescribed heat flux. ---*/ - - Res_Visc[nDim+1] = thermal_conductivity*dTdn*Area; - - /*--- Jacobian contribution for temperature equation. ---*/ - - if (implicit) { - su2double Edge_Vector[3]; - su2double dist_ij_2 = 0, proj_vector_ij = 0; - for (iDim = 0; iDim < nDim; iDim++) { - Edge_Vector[iDim] = Coord_j[iDim]-Coord_i[iDim]; - dist_ij_2 += Edge_Vector[iDim]*Edge_Vector[iDim]; - proj_vector_ij += Edge_Vector[iDim]*Normal[iDim]; - } - if (dist_ij_2 == 0.0) proj_vector_ij = 0.0; - else proj_vector_ij = proj_vector_ij/dist_ij_2; - - Jacobian_i[nDim+1][nDim+1] = -thermal_conductivity*proj_vector_ij; - - Jacobian.SubtractBlock(iPoint, iPoint, Jacobian_i); - } - - /*--- Viscous contribution to the residual at the wall ---*/ - - LinSysRes.SubtractBlock(iPoint, Res_Visc); + /*--- Strong imposition of the temperature on the fluid zone. ---*/ + + LinSysRes.SetBlock_Zero(iPoint, nDim+1); + node[iPoint]->SetSolution_Old(nDim+1, Tconjugate); + node[iPoint]->SetEnergy_ResTruncError_Zero(); } @@ -8618,10 +8563,10 @@ void CIncNSSolver::BC_ConjugateHeat_Interface(CGeometry *geometry, CSolver **sol total_index = iPoint*nVar+iVar; Jacobian.DeleteValsRowi(total_index); } -// if(energy) { -// total_index = iPoint*nVar+nDim+1; -// Jacobian.DeleteValsRowi(total_index); -// } + if(energy) { + total_index = iPoint*nVar+nDim+1; + Jacobian.DeleteValsRowi(total_index); + } } } diff --git a/SU2_CFD/src/solver_structure.cpp b/SU2_CFD/src/solver_structure.cpp index 1ecc24e93b58..681bad8cc2d0 100644 --- a/SU2_CFD/src/solver_structure.cpp +++ b/SU2_CFD/src/solver_structure.cpp @@ -56,6 +56,8 @@ CSolver::CSolver(bool mesh_deform_mode) : System(mesh_deform_mode) { rank = SU2_MPI::GetRank(); size = SU2_MPI::GetSize(); + + adjoint = false; /*--- Set the multigrid level to the finest grid. This can be overwritten in the constructors of the derived classes. ---*/ @@ -2987,6 +2989,48 @@ void CSolver::SetSolution_Gradient_LS(CGeometry *geometry, CConfig *config) { } +void CSolver::SetSolution_Zero(CGeometry *geometry) { + for (unsigned long iPoint = 0; iPoint < geometry->GetnPoint(); iPoint++) { + node[iPoint]->SetSolutionZero(); + } +} + +void CSolver::SetSolution_Old(CGeometry *geometry) { + for (unsigned long iPoint = 0; iPoint < geometry->GetnPoint(); iPoint++) { + node[iPoint]->SetSolution(node[iPoint]->GetSolution_Old()); + } +} + +void CSolver::Add_ExternalOld_To_Solution(CGeometry *geometry) { + for (unsigned long iPoint = 0; iPoint < geometry->GetnPoint(); iPoint++) { + node[iPoint]->AddSolution(node[iPoint]->GetSolution_ExternalOld()); + } +} + +void CSolver::SetExternal_Zero(CGeometry *geometry) { + for (unsigned long iPoint = 0; iPoint < geometry->GetnPoint(); iPoint++) { + node[iPoint]->SetExternalZero(); + } +} + +void CSolver::Add_Solution_To_External(CGeometry *geometry) { + for (unsigned long iPoint = 0; iPoint < geometry->GetnPoint(); iPoint++) { + node[iPoint]->Add_External(node[iPoint]->GetSolution()); + } +} + +void CSolver::Add_Solution_To_ExternalOld(CGeometry *geometry) { + for (unsigned long iPoint = 0; iPoint < geometry->GetnPoint(); iPoint++) { + node[iPoint]->Add_ExternalOld(node[iPoint]->GetSolution()); + } +} + +void CSolver::Set_OldExternal(CGeometry *geometry) { + for (unsigned long iPoint = 0; iPoint < geometry->GetnPoint(); iPoint++) { + node[iPoint]->Set_OldExternal(); + } +} + void CSolver::SetGridVel_Gradient(CGeometry *geometry, CConfig *config) { unsigned short iDim, jDim, iVar, iNeigh; unsigned long iPoint, jPoint; diff --git a/SU2_CFD/src/variables/CDiscAdjVariable.cpp b/SU2_CFD/src/variables/CDiscAdjVariable.cpp index c84c43b8ad4d..def9f45c835f 100644 --- a/SU2_CFD/src/variables/CDiscAdjVariable.cpp +++ b/SU2_CFD/src/variables/CDiscAdjVariable.cpp @@ -42,6 +42,8 @@ CDiscAdjVariable::CDiscAdjVariable() : CVariable() { /*--- Initialize arrays to NULL ---*/ Solution_Direct = NULL; + External = NULL; + Solution_Old = NULL; Sensitivity = NULL; DualTime_Derivative = NULL; @@ -72,6 +74,8 @@ CDiscAdjVariable::CDiscAdjVariable(su2double* val_solution, unsigned short val_n /*--- Initialize arrays to NULL ---*/ Solution_Direct = NULL; + External = NULL; + External_Old = NULL; Sensitivity = NULL; DualTime_Derivative = NULL; @@ -94,6 +98,8 @@ CDiscAdjVariable::CDiscAdjVariable(su2double* val_solution, unsigned short val_n } Solution_Direct = new su2double[nVar]; + External = new su2double[nVar]; + External_Old = new su2double[nVar]; Sensitivity = new su2double[nDim]; @@ -104,10 +110,10 @@ CDiscAdjVariable::CDiscAdjVariable(su2double* val_solution, unsigned short val_n } for (iVar = 0; iVar < nVar; iVar++) { - Solution[iVar] = val_solution[iVar]; + Solution[iVar] = val_solution[iVar]; + External[iVar] = val_solution[iVar]; } - if (dual_time) { for (iVar = 0; iVar < nVar; iVar++) { Solution_time_n[iVar] = 0.0; @@ -140,9 +146,15 @@ CDiscAdjVariable::CDiscAdjVariable(su2double* val_solution, unsigned short val_n } } - if (config->GetMultizone_Problem()) - Set_BGSSolution_k(); + if (config->GetMultizone_Problem()) { + Solution_BGS = new su2double[nVar]; + for (iVar = 0; iVar < nVar; iVar++) { + Solution_BGS[iVar] = 0.0; + } + + Set_BGSSolution_k(); + } } CDiscAdjVariable::~CDiscAdjVariable() { @@ -157,6 +169,8 @@ CDiscAdjVariable::~CDiscAdjVariable() { if (Solution_Geometry_BGS_k != NULL) delete [] Solution_Geometry_BGS_k; if (Solution_Direct != NULL) delete [] Solution_Direct; + if (External != NULL) delete [] External; + if (External_Old != NULL) delete [] External_Old; if (Sensitivity != NULL) delete [] Sensitivity; if (DualTime_Derivative != NULL) delete [] DualTime_Derivative; diff --git a/SU2_CFD/src/variables/CVariable.cpp b/SU2_CFD/src/variables/CVariable.cpp index 1a36edec7b3f..224717cc76fe 100644 --- a/SU2_CFD/src/variables/CVariable.cpp +++ b/SU2_CFD/src/variables/CVariable.cpp @@ -57,7 +57,8 @@ CVariable::CVariable(void) { Residual_Sum = NULL; Solution_Adj_Old = NULL; Solution_BGS_k = NULL; - + Input_AdjIndices = NULL; + Output_AdjIndices = NULL; } CVariable::CVariable(unsigned short val_nvar, CConfig *config) { @@ -79,6 +80,8 @@ CVariable::CVariable(unsigned short val_nvar, CConfig *config) { Residual_Sum = NULL; Solution_Adj_Old = NULL; Solution_BGS_k = NULL; + Input_AdjIndices = NULL; + Output_AdjIndices = NULL; /*--- Initialize the number of solution variables. This version of the constructor will be used primarily for converting the @@ -118,6 +121,8 @@ CVariable::CVariable(unsigned short val_nDim, unsigned short val_nvar, CConfig * Residual_Sum = NULL; Solution_Adj_Old = NULL; Solution_BGS_k = NULL; + Input_AdjIndices = NULL; + Output_AdjIndices = NULL; /*--- Initializate the number of dimension and number of variables ---*/ nDim = val_nDim; @@ -163,6 +168,16 @@ CVariable::CVariable(unsigned short val_nDim, unsigned short val_nvar, CConfig * } } + if(config->GetMultizone_Problem() && config->GetAD_Mode()) { + Input_AdjIndices = new int[nVar]; + Output_AdjIndices = new int[nVar]; + + for (iVar = 0; iVar < nVar; iVar++) { + Input_AdjIndices[iVar] = -1; + Output_AdjIndices[iVar] = -1; + } + } + if (config->GetMultizone_Problem()){ Solution_BGS_k = new su2double[nVar](); } @@ -188,6 +203,8 @@ CVariable::~CVariable(void) { if (Residual_Sum != NULL) delete [] Residual_Sum; if (Solution_Adj_Old != NULL) delete [] Solution_Adj_Old; if (Solution_BGS_k != NULL) delete [] Solution_BGS_k; + if (Input_AdjIndices != NULL) delete [] Input_AdjIndices; + if (Output_AdjIndices != NULL) delete [] Output_AdjIndices; if (Gradient != NULL) { for (iVar = 0; iVar < nVar; iVar++) diff --git a/SU2_DEF/src/SU2_DEF.cpp b/SU2_DEF/src/SU2_DEF.cpp index ad035c7f17e7..a29ad2b050ae 100644 --- a/SU2_DEF/src/SU2_DEF.cpp +++ b/SU2_DEF/src/SU2_DEF.cpp @@ -102,7 +102,7 @@ int main(int argc, char *argv[]) { } /*--- Initialize the configuration of the driver ---*/ - driver_config = new CConfig(config_file_name, SU2_DEF, nZone, false); + driver_config = new CConfig(config_file_name, SU2_DEF, false); /*--- Initialize a char to store the zone filename ---*/ char zone_file_name[MAX_STRING_SIZE]; diff --git a/SU2_DOT/src/SU2_DOT.cpp b/SU2_DOT/src/SU2_DOT.cpp index 0bfca1658ddf..7bd436e4f5af 100644 --- a/SU2_DOT/src/SU2_DOT.cpp +++ b/SU2_DOT/src/SU2_DOT.cpp @@ -106,7 +106,7 @@ int main(int argc, char *argv[]) { } /*--- Initialize the configuration of the driver ---*/ - driver_config = new CConfig(config_file_name, SU2_DOT, nZone, false); + driver_config = new CConfig(config_file_name, SU2_DOT, false); /*--- Initialize a char to store the zone filename ---*/ char zone_file_name[MAX_STRING_SIZE]; diff --git a/SU2_GEO/src/SU2_GEO.cpp b/SU2_GEO/src/SU2_GEO.cpp index 9059ed22de1a..75a5b24eb478 100644 --- a/SU2_GEO/src/SU2_GEO.cpp +++ b/SU2_GEO/src/SU2_GEO.cpp @@ -119,7 +119,7 @@ int main(int argc, char *argv[]) { constructor, the input configuration file is parsed and all options are read and stored. ---*/ - config_container[iZone] = new CConfig(config_file_name, SU2_GEO, nZone, true); + config_container[iZone] = new CConfig(config_file_name, SU2_GEO, true); config_container[iZone]->SetMPICommunicator(MPICommunicator); /*--- Definition of the geometry class to store the primal grid in the partitioning process. ---*/ diff --git a/SU2_MSH/src/SU2_MSH.cpp b/SU2_MSH/src/SU2_MSH.cpp index eb57a65c6ee2..390f84e6459b 100644 --- a/SU2_MSH/src/SU2_MSH.cpp +++ b/SU2_MSH/src/SU2_MSH.cpp @@ -101,7 +101,7 @@ int main(int argc, char *argv[]) { constructor, the input configuration file is parsed and all options are read and stored. ---*/ - config_container[iZone] = new CConfig(config_file_name, SU2_MSH, nZone, true); + config_container[iZone] = new CConfig(config_file_name, SU2_MSH, true); config_container[iZone]->SetMPICommunicator(MPICommunicator); /*--- Definition of the geometry class to store the primal grid in the partitioning process. ---*/ diff --git a/SU2_SOL/src/SU2_SOL.cpp b/SU2_SOL/src/SU2_SOL.cpp index 5ceeb85940cf..6c90b37f54df 100644 --- a/SU2_SOL/src/SU2_SOL.cpp +++ b/SU2_SOL/src/SU2_SOL.cpp @@ -100,7 +100,7 @@ int main(int argc, char *argv[]) { } /*--- Initialize the configuration of the driver ---*/ - driver_config = new CConfig(config_file_name, SU2_SOL, nZone, false); + driver_config = new CConfig(config_file_name, SU2_SOL, false); /*--- Initialize a char to store the zone filename ---*/ char zone_file_name[MAX_STRING_SIZE]; diff --git a/TestCases/coupled_cht/disc_adj_incomp_2d/cht_2d_3cylinders.cfg b/TestCases/coupled_cht/disc_adj_incomp_2d/cht_2d_3cylinders.cfg new file mode 100644 index 000000000000..51f91e815391 --- /dev/null +++ b/TestCases/coupled_cht/disc_adj_incomp_2d/cht_2d_3cylinders.cfg @@ -0,0 +1,79 @@ +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% % +% SU2 configuration file % +% Case description: 2D cylinder array with CHT couplings % +% Author: O. Burghardt, T. Economon % +% Institution: Chair for Scientific Computing, TU Kaiserslautern % +% Date: August 8, 2019 % +% File Version 6.0.1 "Falcon" % +% % +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +% +% Physical governing equations (EULER, NAVIER_STOKES, +% WAVE_EQUATION, HEAT_EQUATION, FEM_ELASTICITY, +% POISSON_EQUATION) +SOLVER= MULTIPHYSICS +% +% Mathematical problem (DIRECT, CONTINUOUS_ADJOINT, DISCRETE_ADJOINT) +MATH_PROBLEM= DISCRETE_ADJOINT +% +% +CONFIG_LIST = (flow_cylinder.cfg, solid_cylinder1.cfg, solid_cylinder2.cfg, solid_cylinder3.cfg) +% +% +MARKER_ZONE_INTERFACE= (cylinder_outer1, cylinder_inner1, cylinder_outer2, cylinder_inner2, cylinder_outer3, cylinder_inner3) +% +% +MARKER_CHT_INTERFACE= (cylinder_outer1, cylinder_inner1, cylinder_outer2, cylinder_inner2, cylinder_outer3, cylinder_inner3) +% +% +TIME_DOMAIN = NO +% +% Number of total iterations +OUTER_ITER = 11 +% +% Mesh input file +MESH_FILENAME= mesh_cht_3cyl.su2 +% +% Mesh input file format (SU2, CGNS, NETCDF_ASCII) +MESH_FORMAT= SU2 + +% These are just default parameters so that we can run SU2_DOT_AD, they have no physical meaning for this test case. + +% ----------------------- DESIGN VARIABLE PARAMETERS --------------------------% +% +% Kind of deformation (NO_DEFORMATION, TRANSLATION, ROTATION, SCALE, +% FFD_SETTING, FFD_NACELLE +% FFD_CONTROL_POINT, FFD_CAMBER, FFD_THICKNESS, FFD_TWIST +% FFD_CONTROL_POINT_2D, FFD_CAMBER_2D, FFD_THICKNESS_2D, FFD_TWIST_2D, +% HICKS_HENNE, SURFACE_BUMP) +DV_KIND= HICKS_HENNE +% +% Marker of the surface in which we are going apply the shape deformation +DV_MARKER= ( cylinder_outer1, cylinder_inner1, cylinder_outer2, cylinder_inner2, cylinder_outer3, cylinder_inner3 ) +% +% Parameters of the shape deformation +% - NO_DEFORMATION ( 1.0 ) +% - TRANSLATION ( x_Disp, y_Disp, z_Disp ), as a unit vector +% - ROTATION ( x_Orig, y_Orig, z_Orig, x_End, y_End, z_End ) +% - SCALE ( 1.0 ) +% - ANGLE_OF_ATTACK ( 1.0 ) +% - FFD_SETTING ( 1.0 ) +% - FFD_CONTROL_POINT ( FFD_BoxTag, i_Ind, j_Ind, k_Ind, x_Disp, y_Disp, z_Disp ) +% - FFD_NACELLE ( FFD_BoxTag, rho_Ind, theta_Ind, phi_Ind, rho_Disp, phi_Disp ) +% - FFD_GULL ( FFD_BoxTag, j_Ind ) +% - FFD_ANGLE_OF_ATTACK ( FFD_BoxTag, 1.0 ) +% - FFD_CAMBER ( FFD_BoxTag, i_Ind, j_Ind ) +% - FFD_THICKNESS ( FFD_BoxTag, i_Ind, j_Ind ) +% - FFD_TWIST ( FFD_BoxTag, j_Ind, x_Orig, y_Orig, z_Orig, x_End, y_End, z_End ) +% - FFD_CONTROL_POINT_2D ( FFD_BoxTag, i_Ind, j_Ind, x_Disp, y_Disp ) +% - FFD_CAMBER_2D ( FFD_BoxTag, i_Ind ) +% - FFD_THICKNESS_2D ( FFD_BoxTag, i_Ind ) +% - FFD_TWIST_2D ( FFD_BoxTag, x_Orig, y_Orig ) +% - HICKS_HENNE ( Lower Surface (0)/Upper Surface (1)/Only one Surface (2), x_Loc ) +% - SURFACE_BUMP ( x_Start, x_End, x_Loc ) +DV_PARAM= (0.0, 0.5) +% +% Value of the shape deformation +DV_VALUE= 0.1 diff --git a/TestCases/coupled_cht/disc_adj_incomp_2d/flow_cylinder.cfg b/TestCases/coupled_cht/disc_adj_incomp_2d/flow_cylinder.cfg new file mode 100644 index 000000000000..651ab1abfe8d --- /dev/null +++ b/TestCases/coupled_cht/disc_adj_incomp_2d/flow_cylinder.cfg @@ -0,0 +1,258 @@ +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% % +% SU2 configuration file % +% Case description: Steady incompressible laminar flow around heated cylinders % +% % +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +% ------------- DIRECT, ADJOINT, AND LINEARIZED PROBLEM DEFINITION ------------% +% +% Physical governing equations (EULER, NAVIER_STOKES, +% WAVE_EQUATION, HEAT_EQUATION, FEM_ELASTICITY, +% POISSON_EQUATION) +SOLVER= INC_NAVIER_STOKES +% +% If Navier-Stokes, kind of turbulent model (NONE, SA) +KIND_TURB_MODEL= NONE +% +% Restart solution (NO, YES) +RESTART_SOL= NO +% +% Objective function in gradient evaluation (DRAG, LIFT, SIDEFORCE, MOMENT_X, +% MOMENT_Y, MOMENT_Z, EFFICIENCY, +% EQUIVALENT_AREA, NEARFIELD_PRESSURE, +% FORCE_X, FORCE_Y, FORCE_Z, THRUST, +% TORQUE, TOTAL_HEATFLUX, +% MAXIMUM_HEATFLUX, INVERSE_DESIGN_PRESSURE, +% INVERSE_DESIGN_HEATFLUX, SURFACE_TOTAL_PRESSURE, +% SURFACE_MASSFLOW, SURFACE_STATIC_PRESSURE, SURFACE_MACH) +% For a weighted sum of objectives: separate by commas, add OBJECTIVE_WEIGHT and MARKER_MONITORING in matching order. +OBJECTIVE_FUNCTION= TOTAL_HEATFLUX +% +% List of weighting values when using more than one OBJECTIVE_FUNCTION. Separate by commas and match with MARKER_MONITORING. +OBJECTIVE_WEIGHT = 1.0 +% +% Read binary restart files (YES, NO) +READ_BINARY_RESTART = YES +% +% Data written to history file +HISTORY_OUTPUT=(ITER, RMS_RES, HEAT ) + +% -------------------- BOUNDARY CONDITION DEFINITION --------------------------% +% +% Farfield boundary marker(s) (NONE = no marker) +MARKER_FAR= ( farfield ) +% +% Marker(s) of the surface to be plotted or designed +MARKER_PLOTTING= (cylinder_outer1, cylinder_outer2, cylinder_outer3) +% +% Marker(s) of the surface where the functional (Cd, Cl, etc.) will be evaluated +MARKER_MONITORING= (cylinder_outer1, cylinder_outer2, cylinder_outer3) + +% ---------------- INCOMPRESSIBLE FLOW CONDITION DEFINITION -------------------% +% +% Density model within the incompressible flow solver. +% Options are CONSTANT (default), BOUSSINESQ, or VARIABLE. If VARIABLE, +% an appropriate fluid model must be selected. +INC_DENSITY_MODEL= VARIABLE +% +% Solve the energy equation in the incompressible flow solver +INC_ENERGY_EQUATION = YES +% +% Initial density for incompressible flows (1.2886 kg/m^3 by default) +INC_DENSITY_INIT= 0.000210322 +% +% Initial velocity for incompressible flows (1.0,0,0 m/s by default) +INC_VELOCITY_INIT= ( 3.40297, 0.0, 0.0 ) +% +% Initial temperature for incompressible flows that include the +% energy equation (288.15 K by default). Value is ignored if +% INC_ENERGY_EQUATION is false. +INC_TEMPERATURE_INIT= 288.15 +% +% Non-dimensionalization scheme for incompressible flows. Options are +% INITIAL_VALUES (default), REFERENCE_VALUES, or DIMENSIONAL. +% INC_*_REF values are ignored unless REFERENCE_VALUES is chosen. +INC_NONDIM= DIMENSIONAL + +% ---- IDEAL GAS, POLYTROPIC, VAN DER WAALS AND PENG ROBINSON CONSTANTS -------% +% +% Fluid model (STANDARD_AIR, IDEAL_GAS, VW_GAS, PR_GAS, +% CONSTANT_DENSITY, INC_IDEAL_GAS) +FLUID_MODEL= INC_IDEAL_GAS +% +% Specific heat at constant pressure, Cp (1004.703 J/kg*K (air)). +% Incompressible fluids with energy eqn. only (CONSTANT_DENSITY, INC_IDEAL_GAS). +SPECIFIC_HEAT_CP= 1004.703 +% +% Molecular weight for an incompressible ideal gas (28.96 g/mol (air) default) +% Incompressible fluids with energy eqn. only (CONSTANT_DENSITY, INC_IDEAL_GAS). +MOLECULAR_WEIGHT= 28.96 + +% --------------------------- VISCOSITY MODEL ---------------------------------% +% +% Viscosity model (SUTHERLAND, CONSTANT_VISCOSITY). +VISCOSITY_MODEL= CONSTANT_VISCOSITY +% +% Molecular Viscosity that would be constant (1.716E-5 by default) +MU_CONSTANT= 1.7893e-05 +% +% Sutherland Viscosity Ref (1.716E-5 default value for AIR SI) +MU_REF= 1.716E-5 +% +% Sutherland Temperature Ref (273.15 K default value for AIR SI) +MU_T_REF= 273.15 +% +% Sutherland constant (110.4 default value for AIR SI) +SUTHERLAND_CONSTANT= 110.4 + +% --------------------------- THERMAL CONDUCTIVITY MODEL ----------------------% +% +% Conductivity model (CONSTANT_CONDUCTIVITY, CONSTANT_PRANDTL). +CONDUCTIVITY_MODEL= CONSTANT_PRANDTL +% +% Molecular Thermal Conductivity that would be constant (0.0257 by default) +KT_CONSTANT= 0.0257 +% +% Laminar Prandtl number (0.72 (air), only for CONSTANT_PRANDTL) +PRANDTL_LAM= 0.72 +% +% Turbulent Prandtl number (0.9 (air), only for CONSTANT_PRANDTL) +PRANDTL_TURB= 0.90 + +% ------------- COMMON PARAMETERS DEFINING THE NUMERICAL METHOD ---------------% +% +% Numerical method for spatial gradients (GREEN_GAUSS, WEIGHTED_LEAST_SQUARES) +NUM_METHOD_GRAD= GREEN_GAUSS +% +% Courant-Friedrichs-Lewy condition of the finest grid +CFL_NUMBER= 10.0 +% +% Adaptive CFL number (NO, YES) +CFL_ADAPT= NO +% +% Parameters of the adaptive CFL number (factor down, factor up, CFL min value, +% CFL max value ) +CFL_ADAPT_PARAM= ( 1.5, 0.5, 10.0, 10000.0 ) +% +% Runge-Kutta alpha coefficients +RK_ALPHA_COEFF= ( 0.66667, 0.66667, 1.000000 ) + +% ------------------------ LINEAR SOLVER DEFINITION ---------------------------% +% +% Linear solver or smoother for implicit formulations (BCGSTAB, FGMRES, SMOOTHER_JACOBI, +% SMOOTHER_ILU, SMOOTHER_LUSGS, +% SMOOTHER_LINELET) +LINEAR_SOLVER= FGMRES +% +% Preconditioner of the Krylov linear solver (ILU, LU_SGS, LINELET, JACOBI) +LINEAR_SOLVER_PREC= ILU +% +% Linael solver ILU preconditioner fill-in level (0 by default) +LINEAR_SOLVER_ILU_FILL_IN= 0 +% +% Minimum error of the linear solver for implicit formulations +LINEAR_SOLVER_ERROR= 1E-15 +% +% Max number of iterations of the linear solver for the implicit formulation +LINEAR_SOLVER_ITER= 5 + +% -------------------- FLOW NUMERICAL METHOD DEFINITION -----------------------% +% +% Convective numerical method (JST, LAX-FRIEDRICH, CUSP, ROE, AUSM, HLLC, +% TURKEL_PREC, MSW) +CONV_NUM_METHOD_FLOW= FDS +% +% Monotonic Upwind Scheme for Conservation Laws (TVD) in the flow equations. +% Required for 2nd order upwind schemes (NO, YES) +MUSCL_FLOW= YES +% +% Slope limiter (NONE, VENKATAKRISHNAN, VENKATAKRISHNAN_WANG, +% BARTH_JESPERSEN, VAN_ALBADA_EDGE) +SLOPE_LIMITER_FLOW= NONE +% +% Time discretization (RUNGE-KUTTA_EXPLICIT, EULER_IMPLICIT, EULER_EXPLICIT) +TIME_DISCRE_FLOW= EULER_IMPLICIT + +% --------------------------- CONVERGENCE PARAMETERS --------------------------% +% +% Min value of the residual (log10 of the residual) +CONV_RESIDUAL_MINVAL= -19 +% +% Start convergence criteria at iteration number +CONV_STARTITER= 10 +% +% Number of elements to apply the criteria +CONV_CAUCHY_ELEMS= 100 +% +% Epsilon to control the series convergence +CONV_CAUCHY_EPS= 1E-6 + +% ------------------------- INPUT/OUTPUT INFORMATION --------------------------% +% +% Restart flow input file +SOLUTION_FILENAME= solution_flow.dat +% +% Restart adjoint input file +SOLUTION_ADJ_FILENAME= solution_adj.dat +% +% Output file format (TECPLOT, TECPLOT_BINARY, PARAVIEW, +% FIELDVIEW, FIELDVIEW_BINARY) +TABULAR_FORMAT= TECPLOT +% +% Output file convergence history (w/o extension) +CONV_FILENAME= history +% +% Output file with the forces breakdown +BREAKDOWN_FILENAME= forces_breakdown.dat +% +% Output file restart flow +RESTART_FILENAME= restart_flow.dat +% +% Output file restart adjoint +RESTART_ADJ_FILENAME= restart_adj.dat +% +% Output file flow (w/o extension) variables +VOLUME_FILENAME= flow +% +% Output file adjoint (w/o extension) variables +VOLUME_ADJ_FILENAME= adjoint +% +% Output Objective function +VALUE_OBJFUNC_FILENAME= of_eval.dat +% +% Output objective function gradient (using continuous adjoint) +GRAD_OBJFUNC_FILENAME= of_grad.dat +% +% Output file surface flow coefficient (w/o extension) +SURFACE_FILENAME= surface_flow +% +% Output file surface adjoint coefficient (w/o extension) +SURFACE_ADJ_FILENAME= surface_adjoint +% +% Writing solution file frequency +WRT_SOL_FREQ= 250 +% +% Writing convergence history frequency +WRT_CON_FREQ= 1 + +% ------------------------ GRID DEFORMATION PARAMETERS ------------------------% +% +% Linear solver or smoother for implicit formulations (FGMRES, RESTARTED_FGMRES, BCGSTAB) +DEFORM_LINEAR_SOLVER= FGMRES +% +% Number of smoothing iterations for mesh deformation +DEFORM_LINEAR_ITER= 200 +% +% Number of nonlinear deformation iterations (surface deformation increments) +DEFORM_NONLINEAR_ITER= 1 +% +% Print the residuals during mesh deformation to the console (YES, NO) +DEFORM_CONSOLE_OUTPUT= YES +% +% Type of element stiffness imposed for FEA mesh deformation (INVERSE_VOLUME, +% WALL_DISTANCE, CONSTANT_STIFFNESS) +DEFORM_STIFFNESS_TYPE= INVERSE_VOLUME +% +% Visualize the deformation (NO, YES) +VISUALIZE_VOLUME_DEF= YES diff --git a/TestCases/coupled_cht/disc_adj_incomp_2d/solid_cylinder1.cfg b/TestCases/coupled_cht/disc_adj_incomp_2d/solid_cylinder1.cfg new file mode 100644 index 000000000000..61e97b996331 --- /dev/null +++ b/TestCases/coupled_cht/disc_adj_incomp_2d/solid_cylinder1.cfg @@ -0,0 +1,186 @@ +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% % +% SU2 configuration file % +% Case description: Steady incompressible laminar flow around heated cylinders % +% % +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +% ------------- DIRECT, ADJOINT, AND LINEARIZED PROBLEM DEFINITION ------------% +% +% Physical governing equations (EULER, NAVIER_STOKES, +% WAVE_EQUATION, HEAT_EQUATION, FEM_ELASTICITY, +% POISSON_EQUATION) +SOLVER= HEAT_EQUATION_FVM +% +% Restart solution (NO, YES) +RESTART_SOL= NO +% +% Objective function in gradient evaluation (DRAG, LIFT, SIDEFORCE, MOMENT_X, +% MOMENT_Y, MOMENT_Z, EFFICIENCY, +% EQUIVALENT_AREA, NEARFIELD_PRESSURE, +% FORCE_X, FORCE_Y, FORCE_Z, THRUST, +% TORQUE, TOTAL_HEATFLUX, +% MAXIMUM_HEATFLUX, INVERSE_DESIGN_PRESSURE, +% INVERSE_DESIGN_HEATFLUX, SURFACE_TOTAL_PRESSURE, +% SURFACE_MASSFLOW, SURFACE_STATIC_PRESSURE, SURFACE_MACH) +% For a weighted sum of objectives: separate by commas, add OBJECTIVE_WEIGHT and MARKER_MONITORING in matching order. +OBJECTIVE_FUNCTION= TOTAL_HEATFLUX +% +% List of weighting values when using more than one OBJECTIVE_FUNCTION. Separate by commas and match with MARKER_MONITORING. +OBJECTIVE_WEIGHT = 1.0 +% +% Read binary restart files (YES, NO) +READ_BINARY_RESTART = YES + +% -------------------- BOUNDARY CONDITION DEFINITION --------------------------% +% +MARKER_ISOTHERMAL= ( core1, 350.0 ) +% +% Marker(s) of the surface to be plotted or designed +MARKER_PLOTTING= (cylinder_inner1 ) +% +% Marker(s) of the surface where the functional (Cd, Cl, etc.) will be evaluated +MARKER_MONITORING= ( NONE ) + +% ---------------- (SOLIDS) CONDUCTION CONDITION DEFINITION -------------------% +% +% We should keep the dimensionalization of the coupled flow solver +INC_NONDIM= DIMENSIONAL +% +% Temperature initialization value +SOLID_TEMPERATURE_INIT= 350.0 +% +% Nettis case: hollow cylinder (air w/ 4x the conductivity) +% +% Solid density (kg/m^3) +SOLID_DENSITY= 0.000210322 +% +% Solid specific heat (J/kg*K) +SPECIFIC_HEAT_CP_SOLID = 1004.703 +% +% Solid thermal conductivity (W/m*K) +THERMAL_CONDUCTIVITY_SOLID= 0.1028 + +% ------------- COMMON PARAMETERS DEFINING THE NUMERICAL METHOD ---------------% +% +% Numerical method for spatial gradients (GREEN_GAUSS, WEIGHTED_LEAST_SQUARES) +NUM_METHOD_GRAD= GREEN_GAUSS +% +% Courant-Friedrichs-Lewy condition of the finest grid +CFL_NUMBER= 10.0 +% +% Adaptive CFL number (NO, YES) +CFL_ADAPT= NO +% +% Parameters of the adaptive CFL number (factor down, factor up, CFL min value, +% CFL max value ) +CFL_ADAPT_PARAM= ( 1.5, 0.5, 10.0, 10000.0 ) +% +% Runge-Kutta alpha coefficients +RK_ALPHA_COEFF= ( 0.66667, 0.66667, 1.000000 ) + +% ------------------------ LINEAR SOLVER DEFINITION ---------------------------% +% +% Linear solver or smoother for implicit formulations (BCGSTAB, FGMRES, SMOOTHER_JACOBI, +% SMOOTHER_ILU, SMOOTHER_LUSGS, +% SMOOTHER_LINELET) +LINEAR_SOLVER= FGMRES +% +% Preconditioner of the Krylov linear solver (ILU, LU_SGS, LINELET, JACOBI) +LINEAR_SOLVER_PREC= ILU +% +% Linael solver ILU preconditioner fill-in level (0 by default) +LINEAR_SOLVER_ILU_FILL_IN= 0 +% +% Minimum error of the linear solver for implicit formulations +LINEAR_SOLVER_ERROR= 1E-15 +% +% Max number of iterations of the linear solver for the implicit formulation +LINEAR_SOLVER_ITER= 5 + +% -------------------- HEAT NUMERICAL METHOD DEFINITION -----------------------% +% +TIME_DISCRE_HEAT= EULER_IMPLICIT + +% --------------------------- CONVERGENCE PARAMETERS --------------------------% +% +% Min value of the residual (log10 of the residual) +CONV_RESIDUAL_MINVAL= -19 +% +% Start convergence criteria at iteration number +CONV_STARTITER= 10 +% +% Number of elements to apply the criteria +CONV_CAUCHY_ELEMS= 100 +% +% Epsilon to control the series convergence +CONV_CAUCHY_EPS= 1E-6 + +% ------------------------- INPUT/OUTPUT INFORMATION --------------------------% +% +% Restart flow input file +SOLUTION_FILENAME= solution_flow.dat +% +% Restart adjoint input file +SOLUTION_ADJ_FILENAME= solution_adj.dat +% +% Output file format (TECPLOT, TECPLOT_BINARY, PARAVIEW, +% FIELDVIEW, FIELDVIEW_BINARY) +TABULAR_FORMAT= TECPLOT +% +% Output file convergence history (w/o extension) +CONV_FILENAME= history +% +% Output file with the forces breakdown +BREAKDOWN_FILENAME= forces_breakdown.dat +% +% Output file restart flow +RESTART_FILENAME= restart_flow.dat +% +% Output file restart adjoint +RESTART_ADJ_FILENAME= restart_adj.dat +% +% Output file flow (w/o extension) variables +VOLUME_FILENAME= flow +% +% Output file adjoint (w/o extension) variables +VOLUME_ADJ_FILENAME= adjoint +% +% Output Objective function +VALUE_OBJFUNC_FILENAME= of_eval.dat +% +% Output objective function gradient (using continuous adjoint) +GRAD_OBJFUNC_FILENAME= of_grad.dat +% +% Output file surface flow coefficient (w/o extension) +SURFACE_FILENAME= surface_flow +% +% Output file surface adjoint coefficient (w/o extension) +SURFACE_ADJ_FILENAME= surface_adjoint +% +% Writing solution file frequency +WRT_SOL_FREQ= 250 +% +% Writing convergence history frequency +WRT_CON_FREQ= 1 + +% ------------------------ GRID DEFORMATION PARAMETERS ------------------------% +% +% Linear solver or smoother for implicit formulations (FGMRES, RESTARTED_FGMRES, BCGSTAB) +DEFORM_LINEAR_SOLVER= FGMRES +% +% Number of smoothing iterations for mesh deformation +DEFORM_LINEAR_ITER= 200 +% +% Number of nonlinear deformation iterations (surface deformation increments) +DEFORM_NONLINEAR_ITER= 1 +% +% Print the residuals during mesh deformation to the console (YES, NO) +DEFORM_CONSOLE_OUTPUT= YES +% +% Type of element stiffness imposed for FEA mesh deformation (INVERSE_VOLUME, +% WALL_DISTANCE, CONSTANT_STIFFNESS) +DEFORM_STIFFNESS_TYPE= INVERSE_VOLUME +% +% Visualize the deformation (NO, YES) +VISUALIZE_VOLUME_DEF= YES diff --git a/TestCases/coupled_cht/disc_adj_incomp_2d/solid_cylinder2.cfg b/TestCases/coupled_cht/disc_adj_incomp_2d/solid_cylinder2.cfg new file mode 100644 index 000000000000..1282a952a7dc --- /dev/null +++ b/TestCases/coupled_cht/disc_adj_incomp_2d/solid_cylinder2.cfg @@ -0,0 +1,196 @@ +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% % +% SU2 configuration file % +% Case description: Steady incompressible laminar flow around heated cylinders % +% % +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +% ------------- DIRECT, ADJOINT, AND LINEARIZED PROBLEM DEFINITION ------------% +% +% Physical governing equations (EULER, NAVIER_STOKES, +% WAVE_EQUATION, HEAT_EQUATION, FEM_ELASTICITY, +% POISSON_EQUATION) +SOLVER= HEAT_EQUATION_FVM +% +% Restart solution (NO, YES) +RESTART_SOL= NO +% +% Objective function in gradient evaluation (DRAG, LIFT, SIDEFORCE, MOMENT_X, +% MOMENT_Y, MOMENT_Z, EFFICIENCY, +% EQUIVALENT_AREA, NEARFIELD_PRESSURE, +% FORCE_X, FORCE_Y, FORCE_Z, THRUST, +% TORQUE, TOTAL_HEATFLUX, +% MAXIMUM_HEATFLUX, INVERSE_DESIGN_PRESSURE, +% INVERSE_DESIGN_HEATFLUX, SURFACE_TOTAL_PRESSURE, +% SURFACE_MASSFLOW, SURFACE_STATIC_PRESSURE, SURFACE_MACH) +% For a weighted sum of objectives: separate by commas, add OBJECTIVE_WEIGHT and MARKER_MONITORING in matching order. +OBJECTIVE_FUNCTION= TOTAL_HEATFLUX +% +% List of weighting values when using more than one OBJECTIVE_FUNCTION. Separate by commas and match with MARKER_MONITORING. +OBJECTIVE_WEIGHT = 1.0 +% +% Read binary restart files (YES, NO) +READ_BINARY_RESTART = YES + +% -------------------- BOUNDARY CONDITION DEFINITION --------------------------% +% +MARKER_ISOTHERMAL= ( core2, 350.0 ) +% +% Marker(s) of the surface to be plotted or designed +MARKER_PLOTTING= (cylinder_inner2) +% +% Marker(s) of the surface where the functional (Cd, Cl, etc.) will be evaluated +MARKER_MONITORING= ( NONE ) + +% ---------------- (SOLIDS) CONDUCTION CONDITION DEFINITION -------------------% +% +% We should keep the dimensionalization of the coupled flow solver +INC_NONDIM= DIMENSIONAL +% +% Temperature initialization value +SOLID_TEMPERATURE_INIT= 350.0 +% +% Nettis case: hollow cylinder (air w/ 4x the conductivity) +% +% Solid density (kg/m^3) +SOLID_DENSITY= 0.000210322 +% +% Solid specific heat (J/kg*K) +SPECIFIC_HEAT_CP_SOLID = 1004.703 +% +% Solid thermal conductivity (W/m*K) +THERMAL_CONDUCTIVITY_SOLID= 0.1028 + +% ------------- COMMON PARAMETERS DEFINING THE NUMERICAL METHOD ---------------% +% +% Numerical method for spatial gradients (GREEN_GAUSS, WEIGHTED_LEAST_SQUARES) +NUM_METHOD_GRAD= GREEN_GAUSS +% +% Courant-Friedrichs-Lewy condition of the finest grid +CFL_NUMBER= 10.0 +% +% Adaptive CFL number (NO, YES) +CFL_ADAPT= NO +% +% Parameters of the adaptive CFL number (factor down, factor up, CFL min value, +% CFL max value ) +CFL_ADAPT_PARAM= ( 1.5, 0.5, 10.0, 10000.0 ) +% +% Runge-Kutta alpha coefficients +RK_ALPHA_COEFF= ( 0.66667, 0.66667, 1.000000 ) +% +% Objective function in gradient evaluation (DRAG, LIFT, SIDEFORCE, MOMENT_X, +% MOMENT_Y, MOMENT_Z, EFFICIENCY, +% EQUIVALENT_AREA, NEARFIELD_PRESSURE, +% FORCE_X, FORCE_Y, FORCE_Z, THRUST, +% TORQUE, TOTAL_HEATFLUX, +% MAXIMUM_HEATFLUX, INVERSE_DESIGN_PRESSURE, +% INVERSE_DESIGN_HEATFLUX, SURFACE_TOTAL_PRESSURE, +% SURFACE_MASSFLOW, SURFACE_STATIC_PRESSURE, SURFACE_MACH) + +% ------------------------ LINEAR SOLVER DEFINITION ---------------------------% +% +% Linear solver or smoother for implicit formulations (BCGSTAB, FGMRES, SMOOTHER_JACOBI, +% SMOOTHER_ILU, SMOOTHER_LUSGS, +% SMOOTHER_LINELET) +LINEAR_SOLVER= FGMRES +% +% Preconditioner of the Krylov linear solver (ILU, LU_SGS, LINELET, JACOBI) +LINEAR_SOLVER_PREC= ILU +% +% Linael solver ILU preconditioner fill-in level (0 by default) +LINEAR_SOLVER_ILU_FILL_IN= 0 +% +% Minimum error of the linear solver for implicit formulations +LINEAR_SOLVER_ERROR= 1E-15 +% +% Max number of iterations of the linear solver for the implicit formulation +LINEAR_SOLVER_ITER= 5 + +% -------------------- HEAT NUMERICAL METHOD DEFINITION -----------------------% +% +% Time discretization (EULER_IMPLICIT, EULER_EXPLICIT) +TIME_DISCRE_HEAT= EULER_IMPLICIT + +% --------------------------- CONVERGENCE PARAMETERS --------------------------% +% +% Min value of the residual (log10 of the residual) +CONV_RESIDUAL_MINVAL= -19 +% +% Start convergence criteria at iteration number +CONV_STARTITER= 10 +% +% Number of elements to apply the criteria +CONV_CAUCHY_ELEMS= 100 +% +% Epsilon to control the series convergence +CONV_CAUCHY_EPS= 1E-6 + +% ------------------------- INPUT/OUTPUT INFORMATION --------------------------% +% +% Restart flow input file +SOLUTION_FILENAME= solution_flow.dat +% +% Restart adjoint input file +SOLUTION_ADJ_FILENAME= solution_adj.dat +% +% Output file format (TECPLOT, TECPLOT_BINARY, PARAVIEW, +% FIELDVIEW, FIELDVIEW_BINARY) +TABULAR_FORMAT= TECPLOT +% +% Output file convergence history (w/o extension) +CONV_FILENAME= history +% +% Output file with the forces breakdown +BREAKDOWN_FILENAME= forces_breakdown.dat +% +% Output file restart flow +RESTART_FILENAME= restart_flow.dat +% +% Output file restart adjoint +RESTART_ADJ_FILENAME= restart_adj.dat +% +% Output file flow (w/o extension) variables +VOLUME_FILENAME= flow +% +% Output file adjoint (w/o extension) variables +VOLUME_ADJ_FILENAME= adjoint +% +% Output Objective function +VALUE_OBJFUNC_FILENAME= of_eval.dat +% +% Output objective function gradient (using continuous adjoint) +GRAD_OBJFUNC_FILENAME= of_grad.dat +% +% Output file surface flow coefficient (w/o extension) +SURFACE_FILENAME= surface_flow +% +% Output file surface adjoint coefficient (w/o extension) +SURFACE_ADJ_FILENAME= surface_adjoint +% +% Writing solution file frequency +WRT_SOL_FREQ= 250 +% +% Writing convergence history frequency +WRT_CON_FREQ= 1 + +% ------------------------ GRID DEFORMATION PARAMETERS ------------------------% +% +% Linear solver or smoother for implicit formulations (FGMRES, RESTARTED_FGMRES, BCGSTAB) +DEFORM_LINEAR_SOLVER= FGMRES +% +% Number of smoothing iterations for mesh deformation +DEFORM_LINEAR_ITER= 200 +% +% Number of nonlinear deformation iterations (surface deformation increments) +DEFORM_NONLINEAR_ITER= 1 +% +% Print the residuals during mesh deformation to the console (YES, NO) +DEFORM_CONSOLE_OUTPUT= YES +% +% Type of element stiffness imposed for FEA mesh deformation (INVERSE_VOLUME, +% WALL_DISTANCE, CONSTANT_STIFFNESS) +DEFORM_STIFFNESS_TYPE= INVERSE_VOLUME +% +% Visualize the deformation (NO, YES) +VISUALIZE_VOLUME_DEF= YES diff --git a/TestCases/coupled_cht/disc_adj_incomp_2d/solid_cylinder3.cfg b/TestCases/coupled_cht/disc_adj_incomp_2d/solid_cylinder3.cfg new file mode 100644 index 000000000000..89c0c5e7a9cc --- /dev/null +++ b/TestCases/coupled_cht/disc_adj_incomp_2d/solid_cylinder3.cfg @@ -0,0 +1,196 @@ +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% % +% SU2 configuration file % +% Case description: Steady incompressible laminar flow around heated cylinders % +% % +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +% ------------- DIRECT, ADJOINT, AND LINEARIZED PROBLEM DEFINITION ------------% +% +% Physical governing equations (EULER, NAVIER_STOKES, +% WAVE_EQUATION, HEAT_EQUATION, FEM_ELASTICITY, +% POISSON_EQUATION) +SOLVER= HEAT_EQUATION_FVM +% +% Restart solution (NO, YES) +RESTART_SOL= NO +% +% Objective function in gradient evaluation (DRAG, LIFT, SIDEFORCE, MOMENT_X, +% MOMENT_Y, MOMENT_Z, EFFICIENCY, +% EQUIVALENT_AREA, NEARFIELD_PRESSURE, +% FORCE_X, FORCE_Y, FORCE_Z, THRUST, +% TORQUE, TOTAL_HEATFLUX, +% MAXIMUM_HEATFLUX, INVERSE_DESIGN_PRESSURE, +% INVERSE_DESIGN_HEATFLUX, SURFACE_TOTAL_PRESSURE, +% SURFACE_MASSFLOW, SURFACE_STATIC_PRESSURE, SURFACE_MACH) +% For a weighted sum of objectives: separate by commas, add OBJECTIVE_WEIGHT and MARKER_MONITORING in matching order. +OBJECTIVE_FUNCTION= TOTAL_HEATFLUX +% +% List of weighting values when using more than one OBJECTIVE_FUNCTION. Separate by commas and match with MARKER_MONITORING. +OBJECTIVE_WEIGHT = 1.0 +% +% Read binary restart files (YES, NO) +READ_BINARY_RESTART = YES + +% -------------------- BOUNDARY CONDITION DEFINITION --------------------------% +% +MARKER_ISOTHERMAL= ( core3, 350.0 ) +% +% Marker(s) of the surface to be plotted or designed +MARKER_PLOTTING= ( cylinder_inner3 ) +% +% Marker(s) of the surface where the functional (Cd, Cl, etc.) will be evaluated +MARKER_MONITORING= ( NONE ) + +% ---------------- (SOLIDS) CONDUCTION CONDITION DEFINITION -------------------% +% +% We should keep the dimensionalization of the coupled flow solver +INC_NONDIM= DIMENSIONAL +% +% Temperature initialization value +SOLID_TEMPERATURE_INIT= 350.0 +% +% Nettis case: hollow cylinder (air w/ 4x the conductivity) +% +% Solid density (kg/m^3) +SOLID_DENSITY= 0.000210322 +% +% Solid specific heat (J/kg*K) +SPECIFIC_HEAT_CP_SOLID = 1004.703 +% +% Solid thermal conductivity (W/m*K) +THERMAL_CONDUCTIVITY_SOLID= 0.1028 + +% ------------- COMMON PARAMETERS DEFINING THE NUMERICAL METHOD ---------------% +% +% Numerical method for spatial gradients (GREEN_GAUSS, WEIGHTED_LEAST_SQUARES) +NUM_METHOD_GRAD= GREEN_GAUSS +% +% Courant-Friedrichs-Lewy condition of the finest grid +CFL_NUMBER= 10.0 +% +% Adaptive CFL number (NO, YES) +CFL_ADAPT= NO +% +% Parameters of the adaptive CFL number (factor down, factor up, CFL min value, +% CFL max value ) +CFL_ADAPT_PARAM= ( 1.5, 0.5, 10.0, 10000.0 ) +% +% Runge-Kutta alpha coefficients +RK_ALPHA_COEFF= ( 0.66667, 0.66667, 1.000000 ) +% +% Objective function in gradient evaluation (DRAG, LIFT, SIDEFORCE, MOMENT_X, +% MOMENT_Y, MOMENT_Z, EFFICIENCY, +% EQUIVALENT_AREA, NEARFIELD_PRESSURE, +% FORCE_X, FORCE_Y, FORCE_Z, THRUST, +% TORQUE, TOTAL_HEATFLUX, +% MAXIMUM_HEATFLUX, INVERSE_DESIGN_PRESSURE, +% INVERSE_DESIGN_HEATFLUX, SURFACE_TOTAL_PRESSURE, +% SURFACE_MASSFLOW, SURFACE_STATIC_PRESSURE, SURFACE_MACH) + +% ------------------------ LINEAR SOLVER DEFINITION ---------------------------% +% +% Linear solver or smoother for implicit formulations (BCGSTAB, FGMRES, SMOOTHER_JACOBI, +% SMOOTHER_ILU, SMOOTHER_LUSGS, +% SMOOTHER_LINELET) +LINEAR_SOLVER= FGMRES +% +% Preconditioner of the Krylov linear solver (ILU, LU_SGS, LINELET, JACOBI) +LINEAR_SOLVER_PREC= ILU +% +% Linael solver ILU preconditioner fill-in level (0 by default) +LINEAR_SOLVER_ILU_FILL_IN= 0 +% +% Minimum error of the linear solver for implicit formulations +LINEAR_SOLVER_ERROR= 1E-15 +% +% Max number of iterations of the linear solver for the implicit formulation +LINEAR_SOLVER_ITER= 5 + +% -------------------- HEAT NUMERICAL METHOD DEFINITION -----------------------% +% +% Time discretization (EULER_IMPLICIT, EULER_EXPLICIT) +TIME_DISCRE_HEAT= EULER_IMPLICIT + +% --------------------------- CONVERGENCE PARAMETERS --------------------------% +% +% Min value of the residual (log10 of the residual) +CONV_RESIDUAL_MINVAL= -19 +% +% Start convergence criteria at iteration number +CONV_STARTITER= 10 +% +% Number of elements to apply the criteria +CONV_CAUCHY_ELEMS= 100 +% +% Epsilon to control the series convergence +CONV_CAUCHY_EPS= 1E-6 + +% ------------------------- INPUT/OUTPUT INFORMATION --------------------------% +% +% Restart flow input file +SOLUTION_FILENAME= solution_flow.dat +% +% Restart adjoint input file +SOLUTION_ADJ_FILENAME= solution_adj.dat +% +% Output file format (TECPLOT, TECPLOT_BINARY, PARAVIEW, +% FIELDVIEW, FIELDVIEW_BINARY) +TABULAR_FORMAT= TECPLOT +% +% Output file convergence history (w/o extension) +CONV_FILENAME= history +% +% Output file with the forces breakdown +BREAKDOWN_FILENAME= forces_breakdown.dat +% +% Output file restart flow +RESTART_FILENAME= restart_flow.dat +% +% Output file restart adjoint +RESTART_ADJ_FILENAME= restart_adj.dat +% +% Output file flow (w/o extension) variables +VOLUME_FILENAME= flow +% +% Output file adjoint (w/o extension) variables +VOLUME_ADJ_FILENAME= adjoint +% +% Output Objective function +VALUE_OBJFUNC_FILENAME= of_eval.dat +% +% Output objective function gradient (using continuous adjoint) +GRAD_OBJFUNC_FILENAME= of_grad.dat +% +% Output file surface flow coefficient (w/o extension) +SURFACE_FILENAME= surface_flow +% +% Output file surface adjoint coefficient (w/o extension) +SURFACE_ADJ_FILENAME= surface_adjoint +% +% Writing solution file frequency +WRT_SOL_FREQ= 250 +% +% Writing convergence history frequency +WRT_CON_FREQ= 1 + +% ------------------------ GRID DEFORMATION PARAMETERS ------------------------% +% +% Linear solver or smoother for implicit formulations (FGMRES, RESTARTED_FGMRES, BCGSTAB) +DEFORM_LINEAR_SOLVER= FGMRES +% +% Number of smoothing iterations for mesh deformation +DEFORM_LINEAR_ITER= 200 +% +% Number of nonlinear deformation iterations (surface deformation increments) +DEFORM_NONLINEAR_ITER= 1 +% +% Print the residuals during mesh deformation to the console (YES, NO) +DEFORM_CONSOLE_OUTPUT= YES +% +% Type of element stiffness imposed for FEA mesh deformation (INVERSE_VOLUME, +% WALL_DISTANCE, CONSTANT_STIFFNESS) +DEFORM_STIFFNESS_TYPE= INVERSE_VOLUME +% +% Visualize the deformation (NO, YES) +VISUALIZE_VOLUME_DEF= YES diff --git a/TestCases/coupled_cht/incomp_2d/cht_2d_3cylinders.cfg b/TestCases/coupled_cht/incomp_2d/cht_2d_3cylinders.cfg new file mode 100644 index 000000000000..88f0428519c2 --- /dev/null +++ b/TestCases/coupled_cht/incomp_2d/cht_2d_3cylinders.cfg @@ -0,0 +1,79 @@ +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% % +% SU2 configuration file % +% Case description: 2D cylinder array with CHT couplings % +% Author: O. Burghardt, T. Economon % +% Institution: Chair for Scientific Computing, TU Kaiserslautern % +% Date: August 8, 2019 % +% File Version 6.0.1 "Falcon" % +% % +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +% +% Physical governing equations (EULER, NAVIER_STOKES, +% WAVE_EQUATION, HEAT_EQUATION, FEM_ELASTICITY, +% POISSON_EQUATION) +SOLVER= MULTIPHYSICS +% +% Mathematical problem (DIRECT, CONTINUOUS_ADJOINT, DISCRETE_ADJOINT) +MATH_PROBLEM= DIRECT +% +% +CONFIG_LIST = (flow_cylinder.cfg, solid_cylinder1.cfg, solid_cylinder2.cfg, solid_cylinder3.cfg) +% +% +MARKER_ZONE_INTERFACE= (cylinder_outer1, cylinder_inner1, cylinder_outer2, cylinder_inner2, cylinder_outer3, cylinder_inner3) +% +% +MARKER_CHT_INTERFACE= (cylinder_outer1, cylinder_inner1, cylinder_outer2, cylinder_inner2, cylinder_outer3, cylinder_inner3) +% +% +TIME_DOMAIN = NO +% +% Number of total iterations +OUTER_ITER = 11 +% +% Mesh input file +MESH_FILENAME= mesh_cht_3cyl.su2 +% +% Mesh input file format (SU2, CGNS, NETCDF_ASCII) +MESH_FORMAT= SU2 + +% These are just default parameters so that we can run SU2_DOT_AD, they have no physical meaning for this test case. + +% ----------------------- DESIGN VARIABLE PARAMETERS --------------------------% +% +% Kind of deformation (NO_DEFORMATION, TRANSLATION, ROTATION, SCALE, +% FFD_SETTING, FFD_NACELLE +% FFD_CONTROL_POINT, FFD_CAMBER, FFD_THICKNESS, FFD_TWIST +% FFD_CONTROL_POINT_2D, FFD_CAMBER_2D, FFD_THICKNESS_2D, FFD_TWIST_2D, +% HICKS_HENNE, SURFACE_BUMP) +DV_KIND= HICKS_HENNE +% +% Marker of the surface in which we are going apply the shape deformation +DV_MARKER= ( cylinder_outer1, cylinder_inner1, cylinder_outer2, cylinder_inner2, cylinder_outer3, cylinder_inner3 ) +% +% Parameters of the shape deformation +% - NO_DEFORMATION ( 1.0 ) +% - TRANSLATION ( x_Disp, y_Disp, z_Disp ), as a unit vector +% - ROTATION ( x_Orig, y_Orig, z_Orig, x_End, y_End, z_End ) +% - SCALE ( 1.0 ) +% - ANGLE_OF_ATTACK ( 1.0 ) +% - FFD_SETTING ( 1.0 ) +% - FFD_CONTROL_POINT ( FFD_BoxTag, i_Ind, j_Ind, k_Ind, x_Disp, y_Disp, z_Disp ) +% - FFD_NACELLE ( FFD_BoxTag, rho_Ind, theta_Ind, phi_Ind, rho_Disp, phi_Disp ) +% - FFD_GULL ( FFD_BoxTag, j_Ind ) +% - FFD_ANGLE_OF_ATTACK ( FFD_BoxTag, 1.0 ) +% - FFD_CAMBER ( FFD_BoxTag, i_Ind, j_Ind ) +% - FFD_THICKNESS ( FFD_BoxTag, i_Ind, j_Ind ) +% - FFD_TWIST ( FFD_BoxTag, j_Ind, x_Orig, y_Orig, z_Orig, x_End, y_End, z_End ) +% - FFD_CONTROL_POINT_2D ( FFD_BoxTag, i_Ind, j_Ind, x_Disp, y_Disp ) +% - FFD_CAMBER_2D ( FFD_BoxTag, i_Ind ) +% - FFD_THICKNESS_2D ( FFD_BoxTag, i_Ind ) +% - FFD_TWIST_2D ( FFD_BoxTag, x_Orig, y_Orig ) +% - HICKS_HENNE ( Lower Surface (0)/Upper Surface (1)/Only one Surface (2), x_Loc ) +% - SURFACE_BUMP ( x_Start, x_End, x_Loc ) +DV_PARAM= (0.0, 0.5) +% +% Value of the shape deformation +DV_VALUE= 0.1 diff --git a/TestCases/coupled_cht/incomp_2d/flow_cylinder.cfg b/TestCases/coupled_cht/incomp_2d/flow_cylinder.cfg new file mode 100644 index 000000000000..651ab1abfe8d --- /dev/null +++ b/TestCases/coupled_cht/incomp_2d/flow_cylinder.cfg @@ -0,0 +1,258 @@ +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% % +% SU2 configuration file % +% Case description: Steady incompressible laminar flow around heated cylinders % +% % +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +% ------------- DIRECT, ADJOINT, AND LINEARIZED PROBLEM DEFINITION ------------% +% +% Physical governing equations (EULER, NAVIER_STOKES, +% WAVE_EQUATION, HEAT_EQUATION, FEM_ELASTICITY, +% POISSON_EQUATION) +SOLVER= INC_NAVIER_STOKES +% +% If Navier-Stokes, kind of turbulent model (NONE, SA) +KIND_TURB_MODEL= NONE +% +% Restart solution (NO, YES) +RESTART_SOL= NO +% +% Objective function in gradient evaluation (DRAG, LIFT, SIDEFORCE, MOMENT_X, +% MOMENT_Y, MOMENT_Z, EFFICIENCY, +% EQUIVALENT_AREA, NEARFIELD_PRESSURE, +% FORCE_X, FORCE_Y, FORCE_Z, THRUST, +% TORQUE, TOTAL_HEATFLUX, +% MAXIMUM_HEATFLUX, INVERSE_DESIGN_PRESSURE, +% INVERSE_DESIGN_HEATFLUX, SURFACE_TOTAL_PRESSURE, +% SURFACE_MASSFLOW, SURFACE_STATIC_PRESSURE, SURFACE_MACH) +% For a weighted sum of objectives: separate by commas, add OBJECTIVE_WEIGHT and MARKER_MONITORING in matching order. +OBJECTIVE_FUNCTION= TOTAL_HEATFLUX +% +% List of weighting values when using more than one OBJECTIVE_FUNCTION. Separate by commas and match with MARKER_MONITORING. +OBJECTIVE_WEIGHT = 1.0 +% +% Read binary restart files (YES, NO) +READ_BINARY_RESTART = YES +% +% Data written to history file +HISTORY_OUTPUT=(ITER, RMS_RES, HEAT ) + +% -------------------- BOUNDARY CONDITION DEFINITION --------------------------% +% +% Farfield boundary marker(s) (NONE = no marker) +MARKER_FAR= ( farfield ) +% +% Marker(s) of the surface to be plotted or designed +MARKER_PLOTTING= (cylinder_outer1, cylinder_outer2, cylinder_outer3) +% +% Marker(s) of the surface where the functional (Cd, Cl, etc.) will be evaluated +MARKER_MONITORING= (cylinder_outer1, cylinder_outer2, cylinder_outer3) + +% ---------------- INCOMPRESSIBLE FLOW CONDITION DEFINITION -------------------% +% +% Density model within the incompressible flow solver. +% Options are CONSTANT (default), BOUSSINESQ, or VARIABLE. If VARIABLE, +% an appropriate fluid model must be selected. +INC_DENSITY_MODEL= VARIABLE +% +% Solve the energy equation in the incompressible flow solver +INC_ENERGY_EQUATION = YES +% +% Initial density for incompressible flows (1.2886 kg/m^3 by default) +INC_DENSITY_INIT= 0.000210322 +% +% Initial velocity for incompressible flows (1.0,0,0 m/s by default) +INC_VELOCITY_INIT= ( 3.40297, 0.0, 0.0 ) +% +% Initial temperature for incompressible flows that include the +% energy equation (288.15 K by default). Value is ignored if +% INC_ENERGY_EQUATION is false. +INC_TEMPERATURE_INIT= 288.15 +% +% Non-dimensionalization scheme for incompressible flows. Options are +% INITIAL_VALUES (default), REFERENCE_VALUES, or DIMENSIONAL. +% INC_*_REF values are ignored unless REFERENCE_VALUES is chosen. +INC_NONDIM= DIMENSIONAL + +% ---- IDEAL GAS, POLYTROPIC, VAN DER WAALS AND PENG ROBINSON CONSTANTS -------% +% +% Fluid model (STANDARD_AIR, IDEAL_GAS, VW_GAS, PR_GAS, +% CONSTANT_DENSITY, INC_IDEAL_GAS) +FLUID_MODEL= INC_IDEAL_GAS +% +% Specific heat at constant pressure, Cp (1004.703 J/kg*K (air)). +% Incompressible fluids with energy eqn. only (CONSTANT_DENSITY, INC_IDEAL_GAS). +SPECIFIC_HEAT_CP= 1004.703 +% +% Molecular weight for an incompressible ideal gas (28.96 g/mol (air) default) +% Incompressible fluids with energy eqn. only (CONSTANT_DENSITY, INC_IDEAL_GAS). +MOLECULAR_WEIGHT= 28.96 + +% --------------------------- VISCOSITY MODEL ---------------------------------% +% +% Viscosity model (SUTHERLAND, CONSTANT_VISCOSITY). +VISCOSITY_MODEL= CONSTANT_VISCOSITY +% +% Molecular Viscosity that would be constant (1.716E-5 by default) +MU_CONSTANT= 1.7893e-05 +% +% Sutherland Viscosity Ref (1.716E-5 default value for AIR SI) +MU_REF= 1.716E-5 +% +% Sutherland Temperature Ref (273.15 K default value for AIR SI) +MU_T_REF= 273.15 +% +% Sutherland constant (110.4 default value for AIR SI) +SUTHERLAND_CONSTANT= 110.4 + +% --------------------------- THERMAL CONDUCTIVITY MODEL ----------------------% +% +% Conductivity model (CONSTANT_CONDUCTIVITY, CONSTANT_PRANDTL). +CONDUCTIVITY_MODEL= CONSTANT_PRANDTL +% +% Molecular Thermal Conductivity that would be constant (0.0257 by default) +KT_CONSTANT= 0.0257 +% +% Laminar Prandtl number (0.72 (air), only for CONSTANT_PRANDTL) +PRANDTL_LAM= 0.72 +% +% Turbulent Prandtl number (0.9 (air), only for CONSTANT_PRANDTL) +PRANDTL_TURB= 0.90 + +% ------------- COMMON PARAMETERS DEFINING THE NUMERICAL METHOD ---------------% +% +% Numerical method for spatial gradients (GREEN_GAUSS, WEIGHTED_LEAST_SQUARES) +NUM_METHOD_GRAD= GREEN_GAUSS +% +% Courant-Friedrichs-Lewy condition of the finest grid +CFL_NUMBER= 10.0 +% +% Adaptive CFL number (NO, YES) +CFL_ADAPT= NO +% +% Parameters of the adaptive CFL number (factor down, factor up, CFL min value, +% CFL max value ) +CFL_ADAPT_PARAM= ( 1.5, 0.5, 10.0, 10000.0 ) +% +% Runge-Kutta alpha coefficients +RK_ALPHA_COEFF= ( 0.66667, 0.66667, 1.000000 ) + +% ------------------------ LINEAR SOLVER DEFINITION ---------------------------% +% +% Linear solver or smoother for implicit formulations (BCGSTAB, FGMRES, SMOOTHER_JACOBI, +% SMOOTHER_ILU, SMOOTHER_LUSGS, +% SMOOTHER_LINELET) +LINEAR_SOLVER= FGMRES +% +% Preconditioner of the Krylov linear solver (ILU, LU_SGS, LINELET, JACOBI) +LINEAR_SOLVER_PREC= ILU +% +% Linael solver ILU preconditioner fill-in level (0 by default) +LINEAR_SOLVER_ILU_FILL_IN= 0 +% +% Minimum error of the linear solver for implicit formulations +LINEAR_SOLVER_ERROR= 1E-15 +% +% Max number of iterations of the linear solver for the implicit formulation +LINEAR_SOLVER_ITER= 5 + +% -------------------- FLOW NUMERICAL METHOD DEFINITION -----------------------% +% +% Convective numerical method (JST, LAX-FRIEDRICH, CUSP, ROE, AUSM, HLLC, +% TURKEL_PREC, MSW) +CONV_NUM_METHOD_FLOW= FDS +% +% Monotonic Upwind Scheme for Conservation Laws (TVD) in the flow equations. +% Required for 2nd order upwind schemes (NO, YES) +MUSCL_FLOW= YES +% +% Slope limiter (NONE, VENKATAKRISHNAN, VENKATAKRISHNAN_WANG, +% BARTH_JESPERSEN, VAN_ALBADA_EDGE) +SLOPE_LIMITER_FLOW= NONE +% +% Time discretization (RUNGE-KUTTA_EXPLICIT, EULER_IMPLICIT, EULER_EXPLICIT) +TIME_DISCRE_FLOW= EULER_IMPLICIT + +% --------------------------- CONVERGENCE PARAMETERS --------------------------% +% +% Min value of the residual (log10 of the residual) +CONV_RESIDUAL_MINVAL= -19 +% +% Start convergence criteria at iteration number +CONV_STARTITER= 10 +% +% Number of elements to apply the criteria +CONV_CAUCHY_ELEMS= 100 +% +% Epsilon to control the series convergence +CONV_CAUCHY_EPS= 1E-6 + +% ------------------------- INPUT/OUTPUT INFORMATION --------------------------% +% +% Restart flow input file +SOLUTION_FILENAME= solution_flow.dat +% +% Restart adjoint input file +SOLUTION_ADJ_FILENAME= solution_adj.dat +% +% Output file format (TECPLOT, TECPLOT_BINARY, PARAVIEW, +% FIELDVIEW, FIELDVIEW_BINARY) +TABULAR_FORMAT= TECPLOT +% +% Output file convergence history (w/o extension) +CONV_FILENAME= history +% +% Output file with the forces breakdown +BREAKDOWN_FILENAME= forces_breakdown.dat +% +% Output file restart flow +RESTART_FILENAME= restart_flow.dat +% +% Output file restart adjoint +RESTART_ADJ_FILENAME= restart_adj.dat +% +% Output file flow (w/o extension) variables +VOLUME_FILENAME= flow +% +% Output file adjoint (w/o extension) variables +VOLUME_ADJ_FILENAME= adjoint +% +% Output Objective function +VALUE_OBJFUNC_FILENAME= of_eval.dat +% +% Output objective function gradient (using continuous adjoint) +GRAD_OBJFUNC_FILENAME= of_grad.dat +% +% Output file surface flow coefficient (w/o extension) +SURFACE_FILENAME= surface_flow +% +% Output file surface adjoint coefficient (w/o extension) +SURFACE_ADJ_FILENAME= surface_adjoint +% +% Writing solution file frequency +WRT_SOL_FREQ= 250 +% +% Writing convergence history frequency +WRT_CON_FREQ= 1 + +% ------------------------ GRID DEFORMATION PARAMETERS ------------------------% +% +% Linear solver or smoother for implicit formulations (FGMRES, RESTARTED_FGMRES, BCGSTAB) +DEFORM_LINEAR_SOLVER= FGMRES +% +% Number of smoothing iterations for mesh deformation +DEFORM_LINEAR_ITER= 200 +% +% Number of nonlinear deformation iterations (surface deformation increments) +DEFORM_NONLINEAR_ITER= 1 +% +% Print the residuals during mesh deformation to the console (YES, NO) +DEFORM_CONSOLE_OUTPUT= YES +% +% Type of element stiffness imposed for FEA mesh deformation (INVERSE_VOLUME, +% WALL_DISTANCE, CONSTANT_STIFFNESS) +DEFORM_STIFFNESS_TYPE= INVERSE_VOLUME +% +% Visualize the deformation (NO, YES) +VISUALIZE_VOLUME_DEF= YES diff --git a/TestCases/coupled_cht/incomp_2d/solid_cylinder1.cfg b/TestCases/coupled_cht/incomp_2d/solid_cylinder1.cfg new file mode 100644 index 000000000000..61e97b996331 --- /dev/null +++ b/TestCases/coupled_cht/incomp_2d/solid_cylinder1.cfg @@ -0,0 +1,186 @@ +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% % +% SU2 configuration file % +% Case description: Steady incompressible laminar flow around heated cylinders % +% % +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +% ------------- DIRECT, ADJOINT, AND LINEARIZED PROBLEM DEFINITION ------------% +% +% Physical governing equations (EULER, NAVIER_STOKES, +% WAVE_EQUATION, HEAT_EQUATION, FEM_ELASTICITY, +% POISSON_EQUATION) +SOLVER= HEAT_EQUATION_FVM +% +% Restart solution (NO, YES) +RESTART_SOL= NO +% +% Objective function in gradient evaluation (DRAG, LIFT, SIDEFORCE, MOMENT_X, +% MOMENT_Y, MOMENT_Z, EFFICIENCY, +% EQUIVALENT_AREA, NEARFIELD_PRESSURE, +% FORCE_X, FORCE_Y, FORCE_Z, THRUST, +% TORQUE, TOTAL_HEATFLUX, +% MAXIMUM_HEATFLUX, INVERSE_DESIGN_PRESSURE, +% INVERSE_DESIGN_HEATFLUX, SURFACE_TOTAL_PRESSURE, +% SURFACE_MASSFLOW, SURFACE_STATIC_PRESSURE, SURFACE_MACH) +% For a weighted sum of objectives: separate by commas, add OBJECTIVE_WEIGHT and MARKER_MONITORING in matching order. +OBJECTIVE_FUNCTION= TOTAL_HEATFLUX +% +% List of weighting values when using more than one OBJECTIVE_FUNCTION. Separate by commas and match with MARKER_MONITORING. +OBJECTIVE_WEIGHT = 1.0 +% +% Read binary restart files (YES, NO) +READ_BINARY_RESTART = YES + +% -------------------- BOUNDARY CONDITION DEFINITION --------------------------% +% +MARKER_ISOTHERMAL= ( core1, 350.0 ) +% +% Marker(s) of the surface to be plotted or designed +MARKER_PLOTTING= (cylinder_inner1 ) +% +% Marker(s) of the surface where the functional (Cd, Cl, etc.) will be evaluated +MARKER_MONITORING= ( NONE ) + +% ---------------- (SOLIDS) CONDUCTION CONDITION DEFINITION -------------------% +% +% We should keep the dimensionalization of the coupled flow solver +INC_NONDIM= DIMENSIONAL +% +% Temperature initialization value +SOLID_TEMPERATURE_INIT= 350.0 +% +% Nettis case: hollow cylinder (air w/ 4x the conductivity) +% +% Solid density (kg/m^3) +SOLID_DENSITY= 0.000210322 +% +% Solid specific heat (J/kg*K) +SPECIFIC_HEAT_CP_SOLID = 1004.703 +% +% Solid thermal conductivity (W/m*K) +THERMAL_CONDUCTIVITY_SOLID= 0.1028 + +% ------------- COMMON PARAMETERS DEFINING THE NUMERICAL METHOD ---------------% +% +% Numerical method for spatial gradients (GREEN_GAUSS, WEIGHTED_LEAST_SQUARES) +NUM_METHOD_GRAD= GREEN_GAUSS +% +% Courant-Friedrichs-Lewy condition of the finest grid +CFL_NUMBER= 10.0 +% +% Adaptive CFL number (NO, YES) +CFL_ADAPT= NO +% +% Parameters of the adaptive CFL number (factor down, factor up, CFL min value, +% CFL max value ) +CFL_ADAPT_PARAM= ( 1.5, 0.5, 10.0, 10000.0 ) +% +% Runge-Kutta alpha coefficients +RK_ALPHA_COEFF= ( 0.66667, 0.66667, 1.000000 ) + +% ------------------------ LINEAR SOLVER DEFINITION ---------------------------% +% +% Linear solver or smoother for implicit formulations (BCGSTAB, FGMRES, SMOOTHER_JACOBI, +% SMOOTHER_ILU, SMOOTHER_LUSGS, +% SMOOTHER_LINELET) +LINEAR_SOLVER= FGMRES +% +% Preconditioner of the Krylov linear solver (ILU, LU_SGS, LINELET, JACOBI) +LINEAR_SOLVER_PREC= ILU +% +% Linael solver ILU preconditioner fill-in level (0 by default) +LINEAR_SOLVER_ILU_FILL_IN= 0 +% +% Minimum error of the linear solver for implicit formulations +LINEAR_SOLVER_ERROR= 1E-15 +% +% Max number of iterations of the linear solver for the implicit formulation +LINEAR_SOLVER_ITER= 5 + +% -------------------- HEAT NUMERICAL METHOD DEFINITION -----------------------% +% +TIME_DISCRE_HEAT= EULER_IMPLICIT + +% --------------------------- CONVERGENCE PARAMETERS --------------------------% +% +% Min value of the residual (log10 of the residual) +CONV_RESIDUAL_MINVAL= -19 +% +% Start convergence criteria at iteration number +CONV_STARTITER= 10 +% +% Number of elements to apply the criteria +CONV_CAUCHY_ELEMS= 100 +% +% Epsilon to control the series convergence +CONV_CAUCHY_EPS= 1E-6 + +% ------------------------- INPUT/OUTPUT INFORMATION --------------------------% +% +% Restart flow input file +SOLUTION_FILENAME= solution_flow.dat +% +% Restart adjoint input file +SOLUTION_ADJ_FILENAME= solution_adj.dat +% +% Output file format (TECPLOT, TECPLOT_BINARY, PARAVIEW, +% FIELDVIEW, FIELDVIEW_BINARY) +TABULAR_FORMAT= TECPLOT +% +% Output file convergence history (w/o extension) +CONV_FILENAME= history +% +% Output file with the forces breakdown +BREAKDOWN_FILENAME= forces_breakdown.dat +% +% Output file restart flow +RESTART_FILENAME= restart_flow.dat +% +% Output file restart adjoint +RESTART_ADJ_FILENAME= restart_adj.dat +% +% Output file flow (w/o extension) variables +VOLUME_FILENAME= flow +% +% Output file adjoint (w/o extension) variables +VOLUME_ADJ_FILENAME= adjoint +% +% Output Objective function +VALUE_OBJFUNC_FILENAME= of_eval.dat +% +% Output objective function gradient (using continuous adjoint) +GRAD_OBJFUNC_FILENAME= of_grad.dat +% +% Output file surface flow coefficient (w/o extension) +SURFACE_FILENAME= surface_flow +% +% Output file surface adjoint coefficient (w/o extension) +SURFACE_ADJ_FILENAME= surface_adjoint +% +% Writing solution file frequency +WRT_SOL_FREQ= 250 +% +% Writing convergence history frequency +WRT_CON_FREQ= 1 + +% ------------------------ GRID DEFORMATION PARAMETERS ------------------------% +% +% Linear solver or smoother for implicit formulations (FGMRES, RESTARTED_FGMRES, BCGSTAB) +DEFORM_LINEAR_SOLVER= FGMRES +% +% Number of smoothing iterations for mesh deformation +DEFORM_LINEAR_ITER= 200 +% +% Number of nonlinear deformation iterations (surface deformation increments) +DEFORM_NONLINEAR_ITER= 1 +% +% Print the residuals during mesh deformation to the console (YES, NO) +DEFORM_CONSOLE_OUTPUT= YES +% +% Type of element stiffness imposed for FEA mesh deformation (INVERSE_VOLUME, +% WALL_DISTANCE, CONSTANT_STIFFNESS) +DEFORM_STIFFNESS_TYPE= INVERSE_VOLUME +% +% Visualize the deformation (NO, YES) +VISUALIZE_VOLUME_DEF= YES diff --git a/TestCases/coupled_cht/incomp_2d/solid_cylinder2.cfg b/TestCases/coupled_cht/incomp_2d/solid_cylinder2.cfg new file mode 100644 index 000000000000..1282a952a7dc --- /dev/null +++ b/TestCases/coupled_cht/incomp_2d/solid_cylinder2.cfg @@ -0,0 +1,196 @@ +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% % +% SU2 configuration file % +% Case description: Steady incompressible laminar flow around heated cylinders % +% % +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +% ------------- DIRECT, ADJOINT, AND LINEARIZED PROBLEM DEFINITION ------------% +% +% Physical governing equations (EULER, NAVIER_STOKES, +% WAVE_EQUATION, HEAT_EQUATION, FEM_ELASTICITY, +% POISSON_EQUATION) +SOLVER= HEAT_EQUATION_FVM +% +% Restart solution (NO, YES) +RESTART_SOL= NO +% +% Objective function in gradient evaluation (DRAG, LIFT, SIDEFORCE, MOMENT_X, +% MOMENT_Y, MOMENT_Z, EFFICIENCY, +% EQUIVALENT_AREA, NEARFIELD_PRESSURE, +% FORCE_X, FORCE_Y, FORCE_Z, THRUST, +% TORQUE, TOTAL_HEATFLUX, +% MAXIMUM_HEATFLUX, INVERSE_DESIGN_PRESSURE, +% INVERSE_DESIGN_HEATFLUX, SURFACE_TOTAL_PRESSURE, +% SURFACE_MASSFLOW, SURFACE_STATIC_PRESSURE, SURFACE_MACH) +% For a weighted sum of objectives: separate by commas, add OBJECTIVE_WEIGHT and MARKER_MONITORING in matching order. +OBJECTIVE_FUNCTION= TOTAL_HEATFLUX +% +% List of weighting values when using more than one OBJECTIVE_FUNCTION. Separate by commas and match with MARKER_MONITORING. +OBJECTIVE_WEIGHT = 1.0 +% +% Read binary restart files (YES, NO) +READ_BINARY_RESTART = YES + +% -------------------- BOUNDARY CONDITION DEFINITION --------------------------% +% +MARKER_ISOTHERMAL= ( core2, 350.0 ) +% +% Marker(s) of the surface to be plotted or designed +MARKER_PLOTTING= (cylinder_inner2) +% +% Marker(s) of the surface where the functional (Cd, Cl, etc.) will be evaluated +MARKER_MONITORING= ( NONE ) + +% ---------------- (SOLIDS) CONDUCTION CONDITION DEFINITION -------------------% +% +% We should keep the dimensionalization of the coupled flow solver +INC_NONDIM= DIMENSIONAL +% +% Temperature initialization value +SOLID_TEMPERATURE_INIT= 350.0 +% +% Nettis case: hollow cylinder (air w/ 4x the conductivity) +% +% Solid density (kg/m^3) +SOLID_DENSITY= 0.000210322 +% +% Solid specific heat (J/kg*K) +SPECIFIC_HEAT_CP_SOLID = 1004.703 +% +% Solid thermal conductivity (W/m*K) +THERMAL_CONDUCTIVITY_SOLID= 0.1028 + +% ------------- COMMON PARAMETERS DEFINING THE NUMERICAL METHOD ---------------% +% +% Numerical method for spatial gradients (GREEN_GAUSS, WEIGHTED_LEAST_SQUARES) +NUM_METHOD_GRAD= GREEN_GAUSS +% +% Courant-Friedrichs-Lewy condition of the finest grid +CFL_NUMBER= 10.0 +% +% Adaptive CFL number (NO, YES) +CFL_ADAPT= NO +% +% Parameters of the adaptive CFL number (factor down, factor up, CFL min value, +% CFL max value ) +CFL_ADAPT_PARAM= ( 1.5, 0.5, 10.0, 10000.0 ) +% +% Runge-Kutta alpha coefficients +RK_ALPHA_COEFF= ( 0.66667, 0.66667, 1.000000 ) +% +% Objective function in gradient evaluation (DRAG, LIFT, SIDEFORCE, MOMENT_X, +% MOMENT_Y, MOMENT_Z, EFFICIENCY, +% EQUIVALENT_AREA, NEARFIELD_PRESSURE, +% FORCE_X, FORCE_Y, FORCE_Z, THRUST, +% TORQUE, TOTAL_HEATFLUX, +% MAXIMUM_HEATFLUX, INVERSE_DESIGN_PRESSURE, +% INVERSE_DESIGN_HEATFLUX, SURFACE_TOTAL_PRESSURE, +% SURFACE_MASSFLOW, SURFACE_STATIC_PRESSURE, SURFACE_MACH) + +% ------------------------ LINEAR SOLVER DEFINITION ---------------------------% +% +% Linear solver or smoother for implicit formulations (BCGSTAB, FGMRES, SMOOTHER_JACOBI, +% SMOOTHER_ILU, SMOOTHER_LUSGS, +% SMOOTHER_LINELET) +LINEAR_SOLVER= FGMRES +% +% Preconditioner of the Krylov linear solver (ILU, LU_SGS, LINELET, JACOBI) +LINEAR_SOLVER_PREC= ILU +% +% Linael solver ILU preconditioner fill-in level (0 by default) +LINEAR_SOLVER_ILU_FILL_IN= 0 +% +% Minimum error of the linear solver for implicit formulations +LINEAR_SOLVER_ERROR= 1E-15 +% +% Max number of iterations of the linear solver for the implicit formulation +LINEAR_SOLVER_ITER= 5 + +% -------------------- HEAT NUMERICAL METHOD DEFINITION -----------------------% +% +% Time discretization (EULER_IMPLICIT, EULER_EXPLICIT) +TIME_DISCRE_HEAT= EULER_IMPLICIT + +% --------------------------- CONVERGENCE PARAMETERS --------------------------% +% +% Min value of the residual (log10 of the residual) +CONV_RESIDUAL_MINVAL= -19 +% +% Start convergence criteria at iteration number +CONV_STARTITER= 10 +% +% Number of elements to apply the criteria +CONV_CAUCHY_ELEMS= 100 +% +% Epsilon to control the series convergence +CONV_CAUCHY_EPS= 1E-6 + +% ------------------------- INPUT/OUTPUT INFORMATION --------------------------% +% +% Restart flow input file +SOLUTION_FILENAME= solution_flow.dat +% +% Restart adjoint input file +SOLUTION_ADJ_FILENAME= solution_adj.dat +% +% Output file format (TECPLOT, TECPLOT_BINARY, PARAVIEW, +% FIELDVIEW, FIELDVIEW_BINARY) +TABULAR_FORMAT= TECPLOT +% +% Output file convergence history (w/o extension) +CONV_FILENAME= history +% +% Output file with the forces breakdown +BREAKDOWN_FILENAME= forces_breakdown.dat +% +% Output file restart flow +RESTART_FILENAME= restart_flow.dat +% +% Output file restart adjoint +RESTART_ADJ_FILENAME= restart_adj.dat +% +% Output file flow (w/o extension) variables +VOLUME_FILENAME= flow +% +% Output file adjoint (w/o extension) variables +VOLUME_ADJ_FILENAME= adjoint +% +% Output Objective function +VALUE_OBJFUNC_FILENAME= of_eval.dat +% +% Output objective function gradient (using continuous adjoint) +GRAD_OBJFUNC_FILENAME= of_grad.dat +% +% Output file surface flow coefficient (w/o extension) +SURFACE_FILENAME= surface_flow +% +% Output file surface adjoint coefficient (w/o extension) +SURFACE_ADJ_FILENAME= surface_adjoint +% +% Writing solution file frequency +WRT_SOL_FREQ= 250 +% +% Writing convergence history frequency +WRT_CON_FREQ= 1 + +% ------------------------ GRID DEFORMATION PARAMETERS ------------------------% +% +% Linear solver or smoother for implicit formulations (FGMRES, RESTARTED_FGMRES, BCGSTAB) +DEFORM_LINEAR_SOLVER= FGMRES +% +% Number of smoothing iterations for mesh deformation +DEFORM_LINEAR_ITER= 200 +% +% Number of nonlinear deformation iterations (surface deformation increments) +DEFORM_NONLINEAR_ITER= 1 +% +% Print the residuals during mesh deformation to the console (YES, NO) +DEFORM_CONSOLE_OUTPUT= YES +% +% Type of element stiffness imposed for FEA mesh deformation (INVERSE_VOLUME, +% WALL_DISTANCE, CONSTANT_STIFFNESS) +DEFORM_STIFFNESS_TYPE= INVERSE_VOLUME +% +% Visualize the deformation (NO, YES) +VISUALIZE_VOLUME_DEF= YES diff --git a/TestCases/coupled_cht/incomp_2d/solid_cylinder3.cfg b/TestCases/coupled_cht/incomp_2d/solid_cylinder3.cfg new file mode 100644 index 000000000000..89c0c5e7a9cc --- /dev/null +++ b/TestCases/coupled_cht/incomp_2d/solid_cylinder3.cfg @@ -0,0 +1,196 @@ +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% % +% SU2 configuration file % +% Case description: Steady incompressible laminar flow around heated cylinders % +% % +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +% ------------- DIRECT, ADJOINT, AND LINEARIZED PROBLEM DEFINITION ------------% +% +% Physical governing equations (EULER, NAVIER_STOKES, +% WAVE_EQUATION, HEAT_EQUATION, FEM_ELASTICITY, +% POISSON_EQUATION) +SOLVER= HEAT_EQUATION_FVM +% +% Restart solution (NO, YES) +RESTART_SOL= NO +% +% Objective function in gradient evaluation (DRAG, LIFT, SIDEFORCE, MOMENT_X, +% MOMENT_Y, MOMENT_Z, EFFICIENCY, +% EQUIVALENT_AREA, NEARFIELD_PRESSURE, +% FORCE_X, FORCE_Y, FORCE_Z, THRUST, +% TORQUE, TOTAL_HEATFLUX, +% MAXIMUM_HEATFLUX, INVERSE_DESIGN_PRESSURE, +% INVERSE_DESIGN_HEATFLUX, SURFACE_TOTAL_PRESSURE, +% SURFACE_MASSFLOW, SURFACE_STATIC_PRESSURE, SURFACE_MACH) +% For a weighted sum of objectives: separate by commas, add OBJECTIVE_WEIGHT and MARKER_MONITORING in matching order. +OBJECTIVE_FUNCTION= TOTAL_HEATFLUX +% +% List of weighting values when using more than one OBJECTIVE_FUNCTION. Separate by commas and match with MARKER_MONITORING. +OBJECTIVE_WEIGHT = 1.0 +% +% Read binary restart files (YES, NO) +READ_BINARY_RESTART = YES + +% -------------------- BOUNDARY CONDITION DEFINITION --------------------------% +% +MARKER_ISOTHERMAL= ( core3, 350.0 ) +% +% Marker(s) of the surface to be plotted or designed +MARKER_PLOTTING= ( cylinder_inner3 ) +% +% Marker(s) of the surface where the functional (Cd, Cl, etc.) will be evaluated +MARKER_MONITORING= ( NONE ) + +% ---------------- (SOLIDS) CONDUCTION CONDITION DEFINITION -------------------% +% +% We should keep the dimensionalization of the coupled flow solver +INC_NONDIM= DIMENSIONAL +% +% Temperature initialization value +SOLID_TEMPERATURE_INIT= 350.0 +% +% Nettis case: hollow cylinder (air w/ 4x the conductivity) +% +% Solid density (kg/m^3) +SOLID_DENSITY= 0.000210322 +% +% Solid specific heat (J/kg*K) +SPECIFIC_HEAT_CP_SOLID = 1004.703 +% +% Solid thermal conductivity (W/m*K) +THERMAL_CONDUCTIVITY_SOLID= 0.1028 + +% ------------- COMMON PARAMETERS DEFINING THE NUMERICAL METHOD ---------------% +% +% Numerical method for spatial gradients (GREEN_GAUSS, WEIGHTED_LEAST_SQUARES) +NUM_METHOD_GRAD= GREEN_GAUSS +% +% Courant-Friedrichs-Lewy condition of the finest grid +CFL_NUMBER= 10.0 +% +% Adaptive CFL number (NO, YES) +CFL_ADAPT= NO +% +% Parameters of the adaptive CFL number (factor down, factor up, CFL min value, +% CFL max value ) +CFL_ADAPT_PARAM= ( 1.5, 0.5, 10.0, 10000.0 ) +% +% Runge-Kutta alpha coefficients +RK_ALPHA_COEFF= ( 0.66667, 0.66667, 1.000000 ) +% +% Objective function in gradient evaluation (DRAG, LIFT, SIDEFORCE, MOMENT_X, +% MOMENT_Y, MOMENT_Z, EFFICIENCY, +% EQUIVALENT_AREA, NEARFIELD_PRESSURE, +% FORCE_X, FORCE_Y, FORCE_Z, THRUST, +% TORQUE, TOTAL_HEATFLUX, +% MAXIMUM_HEATFLUX, INVERSE_DESIGN_PRESSURE, +% INVERSE_DESIGN_HEATFLUX, SURFACE_TOTAL_PRESSURE, +% SURFACE_MASSFLOW, SURFACE_STATIC_PRESSURE, SURFACE_MACH) + +% ------------------------ LINEAR SOLVER DEFINITION ---------------------------% +% +% Linear solver or smoother for implicit formulations (BCGSTAB, FGMRES, SMOOTHER_JACOBI, +% SMOOTHER_ILU, SMOOTHER_LUSGS, +% SMOOTHER_LINELET) +LINEAR_SOLVER= FGMRES +% +% Preconditioner of the Krylov linear solver (ILU, LU_SGS, LINELET, JACOBI) +LINEAR_SOLVER_PREC= ILU +% +% Linael solver ILU preconditioner fill-in level (0 by default) +LINEAR_SOLVER_ILU_FILL_IN= 0 +% +% Minimum error of the linear solver for implicit formulations +LINEAR_SOLVER_ERROR= 1E-15 +% +% Max number of iterations of the linear solver for the implicit formulation +LINEAR_SOLVER_ITER= 5 + +% -------------------- HEAT NUMERICAL METHOD DEFINITION -----------------------% +% +% Time discretization (EULER_IMPLICIT, EULER_EXPLICIT) +TIME_DISCRE_HEAT= EULER_IMPLICIT + +% --------------------------- CONVERGENCE PARAMETERS --------------------------% +% +% Min value of the residual (log10 of the residual) +CONV_RESIDUAL_MINVAL= -19 +% +% Start convergence criteria at iteration number +CONV_STARTITER= 10 +% +% Number of elements to apply the criteria +CONV_CAUCHY_ELEMS= 100 +% +% Epsilon to control the series convergence +CONV_CAUCHY_EPS= 1E-6 + +% ------------------------- INPUT/OUTPUT INFORMATION --------------------------% +% +% Restart flow input file +SOLUTION_FILENAME= solution_flow.dat +% +% Restart adjoint input file +SOLUTION_ADJ_FILENAME= solution_adj.dat +% +% Output file format (TECPLOT, TECPLOT_BINARY, PARAVIEW, +% FIELDVIEW, FIELDVIEW_BINARY) +TABULAR_FORMAT= TECPLOT +% +% Output file convergence history (w/o extension) +CONV_FILENAME= history +% +% Output file with the forces breakdown +BREAKDOWN_FILENAME= forces_breakdown.dat +% +% Output file restart flow +RESTART_FILENAME= restart_flow.dat +% +% Output file restart adjoint +RESTART_ADJ_FILENAME= restart_adj.dat +% +% Output file flow (w/o extension) variables +VOLUME_FILENAME= flow +% +% Output file adjoint (w/o extension) variables +VOLUME_ADJ_FILENAME= adjoint +% +% Output Objective function +VALUE_OBJFUNC_FILENAME= of_eval.dat +% +% Output objective function gradient (using continuous adjoint) +GRAD_OBJFUNC_FILENAME= of_grad.dat +% +% Output file surface flow coefficient (w/o extension) +SURFACE_FILENAME= surface_flow +% +% Output file surface adjoint coefficient (w/o extension) +SURFACE_ADJ_FILENAME= surface_adjoint +% +% Writing solution file frequency +WRT_SOL_FREQ= 250 +% +% Writing convergence history frequency +WRT_CON_FREQ= 1 + +% ------------------------ GRID DEFORMATION PARAMETERS ------------------------% +% +% Linear solver or smoother for implicit formulations (FGMRES, RESTARTED_FGMRES, BCGSTAB) +DEFORM_LINEAR_SOLVER= FGMRES +% +% Number of smoothing iterations for mesh deformation +DEFORM_LINEAR_ITER= 200 +% +% Number of nonlinear deformation iterations (surface deformation increments) +DEFORM_NONLINEAR_ITER= 1 +% +% Print the residuals during mesh deformation to the console (YES, NO) +DEFORM_CONSOLE_OUTPUT= YES +% +% Type of element stiffness imposed for FEA mesh deformation (INVERSE_VOLUME, +% WALL_DISTANCE, CONSTANT_STIFFNESS) +DEFORM_STIFFNESS_TYPE= INVERSE_VOLUME +% +% Visualize the deformation (NO, YES) +VISUALIZE_VOLUME_DEF= YES diff --git a/TestCases/coupled_cht/incompressible/config.cfg b/TestCases/coupled_cht/incompressible/config.cfg deleted file mode 100644 index 0aa8e2473331..000000000000 --- a/TestCases/coupled_cht/incompressible/config.cfg +++ /dev/null @@ -1,25 +0,0 @@ -%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% -% % -% SU2 configuration file % -% Case description: 2D Cylinder test case for CHT coupling % -% Author: Ole Burghardt % -% Institution: Chair for Scientific Computing, TU Kaiserslautern % -% Date: March 12th, 2018 % -% File Version 6.0.1 "Falcon" % -% % -%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% -SOLVER= MULTIPHYSICS - -CONFIG_LIST = (configFlow.cfg, configSolid.cfg) - -MARKER_ZONE_INTERFACE= (PIN, PINSD) - -MESH_FILENAME= coupled_cht_cylinder2d.su2 - -TIME_DOMAIN = NO - -OUTER_ITER = 11 - -% Only required by the python scripts -MATH_PROBLEM = DIRECT -WRT_SOL_FREQ = 100 diff --git a/TestCases/disc_adj_incomp_navierstokes/cylinder/heated_cylinder.cfg b/TestCases/disc_adj_incomp_navierstokes/cylinder/heated_cylinder.cfg index 2a3746d8f612..7608733a0ab1 100644 --- a/TestCases/disc_adj_incomp_navierstokes/cylinder/heated_cylinder.cfg +++ b/TestCases/disc_adj_incomp_navierstokes/cylinder/heated_cylinder.cfg @@ -312,7 +312,7 @@ WRT_SOL_FREQ= 250 WRT_CON_FREQ= 1 % % Screen output -SCREEN_OUTPUT= (INNER_ITER, RMS_ADJ_PRESSURE, RMS_ADJ_HEAT, SENS_PRESS, SENS_AOA) +SCREEN_OUTPUT= (INNER_ITER, RMS_ADJ_PRESSURE, RMS_ADJ_TEMPERATURE, SENS_PRESS, SENS_AOA) % ----------------------- DESIGN VARIABLE PARAMETERS --------------------------% % % Kind of deformation (NO_DEFORMATION, TRANSLATION, ROTATION, SCALE, diff --git a/TestCases/parallel_regression.py b/TestCases/parallel_regression.py index 520867292dd0..f0477b8a9577 100644 --- a/TestCases/parallel_regression.py +++ b/TestCases/parallel_regression.py @@ -1111,11 +1111,11 @@ def main(): # CHT incompressible cht_incompressible = TestCase('cht_incompressible') - cht_incompressible.cfg_dir = "coupled_cht/incompressible" - cht_incompressible.cfg_file = "config.cfg" + cht_incompressible.cfg_dir = "coupled_cht/incomp_2d" + cht_incompressible.cfg_file = "cht_2d_3cylinders.cfg" cht_incompressible.test_iter = 10 - cht_incompressible.test_vals = [10.000000, -2.607260, -2.708775] #last 4 columns - cht_incompressible.su2_exec = "mpirun -n 2 SU2_CFD" + cht_incompressible.test_vals = [-2.132187, -0.579649, -0.579649, -0.579649] #last 4 columns + cht_incompressible.su2_exec = "SU2_CFD" cht_incompressible.timeout = 1600 cht_incompressible.multizone = True cht_incompressible.tol = 0.00001 diff --git a/TestCases/parallel_regression_AD.py b/TestCases/parallel_regression_AD.py index c7debc389a0c..bd1995839977 100644 --- a/TestCases/parallel_regression_AD.py +++ b/TestCases/parallel_regression_AD.py @@ -254,7 +254,7 @@ def main(): discadj_heat.cfg_dir = "disc_adj_heat" discadj_heat.cfg_file = "disc_adj_heat.cfg" discadj_heat.test_iter = 10 - discadj_heat.test_vals = [3.183906, 0.923840, -223.200000, -2059.800000] #last 4 columns + discadj_heat.test_vals = [-2.281765, 0.706785, -0.743990, -6.866000] #last 4 columns discadj_heat.su2_exec = "parallel_computation.py -f" discadj_heat.timeout = 1600 discadj_heat.tol = 0.00001 @@ -275,6 +275,21 @@ def main(): # discadj_fsi.tol = 0.00001 # test_list.append(discadj_fsi) + ################################### + ### Coupled CHT Adjoint ### + ################################### + + # Coupled discrete adjoint for heatflux in heated cylinder array + discadj_cht = TestCase('discadj_cht') + discadj_cht.cfg_dir = "coupled_cht/disc_adj_incomp_2d" + discadj_cht.cfg_file = "cht_2d_3cylinders.cfg" + discadj_cht.test_iter = 10 + discadj_cht.test_vals = [-2.403180, -3.097866, -3.097837, -3.095571] #last 4 columns + discadj_cht.su2_exec = "parallel_computation.py -f" + discadj_cht.timeout = 1600 + discadj_cht.tol = 0.00001 + test_list.append(discadj_cht) + ###################################### ### RUN TESTS ### ###################################### diff --git a/TestCases/serial_regression.py b/TestCases/serial_regression.py index c2bb796f75ff..ac65ead2a7fc 100644 --- a/TestCases/serial_regression.py +++ b/TestCases/serial_regression.py @@ -1207,10 +1207,10 @@ def main(): # CHT incompressible cht_incompressible = TestCase('cht_incompressible') - cht_incompressible.cfg_dir = "coupled_cht/incompressible" - cht_incompressible.cfg_file = "config.cfg" + cht_incompressible.cfg_dir = "coupled_cht/incomp_2d" + cht_incompressible.cfg_file = "cht_2d_3cylinders.cfg" cht_incompressible.test_iter = 10 - cht_incompressible.test_vals = [10.000000, -2.601006, -3.342894] #last 4 columns + cht_incompressible.test_vals = [-2.132187, -0.579649, -0.579649, -0.579649] #last 4 columns cht_incompressible.su2_exec = "SU2_CFD" cht_incompressible.timeout = 1600 cht_incompressible.multizone = True diff --git a/TestCases/serial_regression_AD.py b/TestCases/serial_regression_AD.py index bbe94ac14f5d..d0a9be91b4ff 100644 --- a/TestCases/serial_regression_AD.py +++ b/TestCases/serial_regression_AD.py @@ -239,7 +239,7 @@ def main(): discadj_heat.cfg_dir = "disc_adj_heat" discadj_heat.cfg_file = "disc_adj_heat.cfg" discadj_heat.test_iter = 10 - discadj_heat.test_vals = [3.139355, 1.144919, -1040.600000, -2464.900000] #last 4 columns + discadj_heat.test_vals = [-2.271573, 0.671194, -3.172000, -8.231600] #last 4 columns discadj_heat.su2_exec = "SU2_CFD_AD" discadj_heat.timeout = 1600 discadj_heat.tol = 0.00001 @@ -258,7 +258,22 @@ def main(): # discadj_fsi.su2_exec = "SU2_CFD_AD" # discadj_fsi.timeout = 1600 # discadj_fsi.tol = 0.00001 -# test_list.append(discadj_fsi) +# test_list.append(discadj_fsi) + + ################################### + ### Coupled CHT Adjoint ### + ################################### + + # Coupled discrete adjoint for heatflux in heated cylinder array + discadj_cht = TestCase('discadj_cht') + discadj_cht.cfg_dir = "coupled_cht/disc_adj_incomp_2d" + discadj_cht.cfg_file = "cht_2d_3cylinders.cfg" + discadj_cht.test_iter = 10 + discadj_cht.test_vals = [-2.403785, -3.097865, -3.097836, -3.097833] #last 4 columns + discadj_cht.su2_exec = "parallel_computation.py -f" + discadj_cht.timeout = 1600 + discadj_cht.tol = 0.00001 + test_list.append(discadj_cht) ###################################### ### RUN TESTS ###