diff --git a/Common/include/CConfig.hpp b/Common/include/CConfig.hpp index 7bfb9cb3173..71878276052 100644 --- a/Common/include/CConfig.hpp +++ b/Common/include/CConfig.hpp @@ -138,6 +138,7 @@ class CConfig { su2double Buffet_lambda; /*!< \brief Offset parameter for buffet sensor.*/ su2double Damp_Engine_Inflow; /*!< \brief Damping factor for the engine inlet. */ su2double Damp_Engine_Exhaust; /*!< \brief Damping factor for the engine exhaust. */ + unsigned long Bc_Eval_Freq; /*!< \brief Evaluation frequency for Engine and Actuator disk markers. */ su2double Damp_Res_Restric, /*!< \brief Damping factor for the residual restriction. */ Damp_Correc_Prolong; /*!< \brief Damping factor for the correction prolongation. */ su2double Position_Plane; /*!< \brief Position of the Near-Field (y coordinate 2D, and z coordinate 3D). */ @@ -192,6 +193,7 @@ class CConfig { nMarker_Fluid_Load, /*!< \brief Number of markers in which the flow load is computed/employed. */ nMarker_Fluid_InterfaceBound, /*!< \brief Number of fluid interface markers. */ nMarker_CHTInterface, /*!< \brief Number of conjugate heat transfer interface markers. */ + nMarker_ContactResistance, /*!< \brief Number of CHT interfaces with contact resistance. */ nMarker_Inlet, /*!< \brief Number of inlet flow markers. */ nMarker_Inlet_Species, /*!< \brief Number of inlet species markers. */ nSpecies_per_Inlet, /*!< \brief Number of species defined per inlet markers. */ @@ -235,6 +237,7 @@ class CConfig { *Marker_MixingPlaneInterface, /*!< \brief MixingPlane interface boundary markers. */ *Marker_TurboBoundIn, /*!< \brief Turbomachinery performance boundary markers. */ *Marker_TurboBoundOut, /*!< \brief Turbomachinery performance boundary donor markers. */ + *Marker_Turbomachinery, /*!< \breif Turbomachinery markers */ *Marker_NearFieldBound, /*!< \brief Near Field boundaries markers. */ *Marker_Deform_Mesh, /*!< \brief Deformable markers at the boundary. */ *Marker_Deform_Mesh_Sym_Plane, /*!< \brief Marker with symmetric deformation. */ @@ -396,6 +399,7 @@ class CConfig { su2double **Periodic_RotCenter; /*!< \brief Rotational center for each periodic boundary. */ su2double **Periodic_RotAngles; /*!< \brief Rotation angles for each periodic boundary. */ su2double **Periodic_Translation; /*!< \brief Translation vector for each periodic boundary. */ + su2double *CHT_ContactResistance; /*!< \brief Contact resistance values for each solid-solid CHT interface. */ string *Marker_CfgFile_TagBound; /*!< \brief Global index for markers using config file. */ unsigned short *Marker_All_KindBC, /*!< \brief Global index for boundaries using grid information. */ *Marker_CfgFile_KindBC; /*!< \brief Global index for boundaries using config file. */ @@ -440,6 +444,7 @@ class CConfig { TURBO_PERF_KIND *Kind_TurboPerf; /*!< \brief Kind of turbomachynery architecture.*/ TURBOMACHINERY_TYPE *Kind_TurboMachinery; + su2vector Kind_TurboInterface; /* Turbomachinery objective functions */ su2double *EntropyGeneration; @@ -468,6 +473,7 @@ class CConfig { unsigned short* nDV_Value; /*!< \brief Number of values for each design variable (might be different than 1 if we allow arbitrary movement). */ unsigned short nFFDBox; /*!< \brief Number of ffd boxes. */ unsigned short nTurboMachineryKind; /*!< \brief Number turbomachinery types specified. */ + unsigned short nTurboInterfaces; /*!< \brief Number of turbomachiery interfaces */ unsigned short nParamDV; /*!< \brief Number of parameters of the design variable. */ string DV_Filename; /*!< \brief Filename for providing surface positions from an external parameterization. */ string DV_Unordered_Sens_Filename; /*!< \brief Filename of volume sensitivities in an unordered ASCII format. */ @@ -594,6 +600,7 @@ class CConfig { bool EulerPersson; /*!< \brief Boolean to determine whether this is an Euler simulation with Persson shock capturing. */ bool FSI_Problem = false,/*!< \brief Boolean to determine whether the simulation is FSI or not. */ Multizone_Problem; /*!< \brief Boolean to determine whether we are solving a multizone problem. */ + //bool ContactResistance = false; /*!< \brief Apply contact resistance for conjugate heat transfer. */ unsigned short nID_DV; /*!< \brief ID for the region of FEM when computed using direct differentiation. */ bool AD_Mode; /*!< \brief Algorithmic Differentiation support. */ @@ -747,6 +754,7 @@ class CConfig { *Marker_All_Turbomachinery, /*!< \brief Global index for Turbomachinery markers using the grid information. */ *Marker_All_TurbomachineryFlag, /*!< \brief Global index for Turbomachinery markers flag using the grid information. */ *Marker_All_MixingPlaneInterface, /*!< \brief Global index for MixingPlane interface markers using the grid information. */ + *Marker_All_Giles, /*!< \brief Global index for Giles markers using the grid information. */ *Marker_All_DV, /*!< \brief Global index for design variable markers using the grid information. */ *Marker_All_Moving, /*!< \brief Global index for moving surfaces using the grid information. */ *Marker_All_Deform_Mesh, /*!< \brief Global index for deformable markers at the boundary. */ @@ -764,6 +772,7 @@ class CConfig { *Marker_CfgFile_Turbomachinery, /*!< \brief Global index for Turbomachinery using the config information. */ *Marker_CfgFile_TurbomachineryFlag, /*!< \brief Global index for Turbomachinery flag using the config information. */ *Marker_CfgFile_MixingPlaneInterface, /*!< \brief Global index for MixingPlane interface using the config information. */ + *Marker_CfgFile_Giles, /*!< \brief Global index for Giles markers flag using the config information. */ *Marker_CfgFile_Moving, /*!< \brief Global index for moving surfaces using the config information. */ *Marker_CfgFile_Deform_Mesh, /*!< \brief Global index for deformable markers at the boundary. */ *Marker_CfgFile_Deform_Mesh_Sym_Plane, /*!< \brief Global index for markers with symmetric deformations. */ @@ -1048,7 +1057,8 @@ class CConfig { long ParMETIS_pointWgt; /*!< \brief Load balancing weight given to points. */ long ParMETIS_edgeWgt; /*!< \brief Load balancing weight given to edges. */ unsigned short DirectDiff; /*!< \brief Direct Differentation mode. */ - bool DiscreteAdjoint; /*!< \brief AD-based discrete adjoint mode. */ + bool DiscreteAdjoint, /*!< \brief AD-based discrete adjoint mode. */ + DiscreteAdjointDebug; /*!< \brief Discrete adjoint debug mode using tags. */ su2double Const_DES; /*!< \brief Detached Eddy Simulation Constant. */ WINDOW_FUNCTION Kind_WindowFct; /*!< \brief Type of window (weight) function for objective functional. */ unsigned short Kind_HybridRANSLES; /*!< \brief Kind of Hybrid RANS/LES. */ @@ -1376,7 +1386,7 @@ class CConfig { su2double** & RotCenter, su2double** & RotAngles, su2double** & Translation); void addTurboPerfOption(const string & name, unsigned short & nMarker_TurboPerf, - string* & Marker_TurboBoundIn, string* & Marker_TurboBoundOut); + string* & Marker_TurboBoundIn, string* & Marker_TurboBoundOut, string* & Marker_Turbomachinery); void addActDiskOption(const string & name, unsigned short & nMarker_ActDiskInlet, unsigned short & nMarker_ActDiskOutlet, string* & Marker_ActDiskInlet, string* & Marker_ActDiskOutlet, @@ -3514,6 +3524,13 @@ class CConfig { */ void SetMarker_All_MixingPlaneInterface(unsigned short val_marker, unsigned short val_mixpla_interface) { Marker_All_MixingPlaneInterface[val_marker] = val_mixpla_interface; } + /*! + * \brief Set if a marker val_marker is part of the Giles boundary (read from the config file). + * \param[in] val_marker - Index of the marker in which we are interested. + * \param[in] val_giles - 0 if not part of the Giles boundary or greater than 1 if it is part. + */ + void SetMarker_All_Giles(unsigned short val_marker, unsigned short val_giles) { Marker_All_Giles[val_marker] = val_giles; } + /*! * \brief Set if a marker val_marker is going to be affected by design variables val_moving * (read from the config file). @@ -3660,6 +3677,13 @@ class CConfig { */ unsigned short GetMarker_All_TurbomachineryFlag(unsigned short val_marker) const { return Marker_All_TurbomachineryFlag[val_marker]; } +/*! + * \brief Get the Giles boundary information for a marker val_marker. + * \param[in] val_marker value of the marker on the grid. + * \return 0 if is not part of the MixingPlane Interface and greater than 1 if it is part. + */ + unsigned short GetMarker_All_Giles(unsigned short val_marker) const { return Marker_All_Giles[val_marker]; } + /*! * \brief Get the number of FSI interface markers val_marker. * \param[in] void. @@ -3667,6 +3691,13 @@ class CConfig { */ unsigned short GetMarker_n_ZoneInterface(void) const { return nMarker_ZoneInterface; } + /*! + * \brief Get the contact resistance value of a specified interface. + * \param[in] val_interface interface index. + * \return Contact resistance value (zero by default). + */ + su2double GetContactResistance(unsigned short val_interface) const { return (nMarker_ContactResistance > 0) ? CHT_ContactResistance[val_interface] : 0.0; } + /*! * \brief Get the DV information for a marker val_marker. * \param[in] val_marker - 0 or 1 depending if the the marker is going to be affected by design variables. @@ -5329,6 +5360,12 @@ class CConfig { */ TURBO_PERF_KIND GetKind_TurboPerf(unsigned short val_iZone) const { return Kind_TurboPerf[val_iZone]; }; + /*! + * \brief gets interface kind for an interface marker in turbomachinery problem + * \return interface kind + */ + TURBO_INTERFACE_KIND GetKind_TurboInterface(unsigned short interfaceIndex) const { return Kind_TurboInterface[interfaceIndex]; } + /*! * \brief get outlet bounds name for Turbomachinery performance calculation. * \return name of the bound. @@ -6392,6 +6429,12 @@ class CConfig { */ unsigned short GetMarker_CfgFile_MixingPlaneInterface(const string& val_marker) const; + /*! + * \brief Get the Giles boundary information from the config definition for the marker val_marker. + * \return Plotting information of the boundary in the config information for the marker val_marker. + */ + unsigned short GetMarker_CfgFile_Giles(const string& val_marker) const; + /*! * \brief Get the DV information from the config definition for the marker val_marker. * \return DV information of the boundary in the config information for the marker val_marker. @@ -6510,6 +6553,12 @@ class CConfig { */ su2double GetMinLogResidual(void) const { return MinLogResidual; } + /*! + * \brief Evaluation frequency for Engine and Actuator disk markers. + * \return Value Evaluation frequency . + */ + unsigned long GetBc_Eval_Freq(void) const { return Bc_Eval_Freq; } + /*! * \brief Value of the damping factor for the engine inlet bc. * \return Value of the damping factor. @@ -8901,6 +8950,12 @@ class CConfig { */ bool GetDiscrete_Adjoint(void) const { return DiscreteAdjoint; } + /*! + * \brief Get the indicator whether a debug run for the discrete adjoint solver will be started. + * \return the discrete adjoint debug indicator. + */ + bool GetDiscrete_Adjoint_Debug(void) const { return DiscreteAdjointDebug; } + /*! * \brief Get the number of subiterations while a ramp is applied. * \return Number of FSI subiters. diff --git a/Common/include/basic_types/ad_structure.hpp b/Common/include/basic_types/ad_structure.hpp index ecd2c93b1c1..5861faccd1c 100644 --- a/Common/include/basic_types/ad_structure.hpp +++ b/Common/include/basic_types/ad_structure.hpp @@ -270,6 +270,12 @@ inline void SetExtFuncOut(T&& data, const int size_x, const int size_y) {} */ inline void SetIndex(Identifier& index, const su2double& data) {} +/*! + * \brief Sets the tag tape to a specific tag. + * \param[in] tag - the number to which the tag is set. + */ +inline void SetTag(int tag) {} + /*! * \brief Pushes back the current tape position to the tape position's vector. */ @@ -482,6 +488,8 @@ FORCEINLINE void EndUseAdjoints() { AD::getTape().endUseAdjointVector(); } FORCEINLINE void SetIndex(Identifier& index, const su2double& data) { index = data.getIdentifier(); } +FORCEINLINE void SetTag(int tag) { AD::getTape().setCurTag(tag); } + // WARNING: For performance reasons, this method does not perform bounds checking. // When using it, please ensure sufficient adjoint vector size by a call to AD::ResizeAdjoints(). // This method does not perform locking either. @@ -514,7 +522,8 @@ FORCEINLINE void SetPreaccIn() {} template ::value> = 0> FORCEINLINE void SetPreaccIn(const T& data, Ts&&... moreData) { if (!PreaccActive) return; - if (IsIdentifierActive(data)) PreaccHelper.addInput(data); + // if (IsIdentifierActive(data)) + PreaccHelper.addInput(data); SetPreaccIn(moreData...); } @@ -527,9 +536,9 @@ template FORCEINLINE void SetPreaccIn(const T& data, const int size) { if (PreaccActive) { for (int i = 0; i < size; i++) { - if (IsIdentifierActive(data[i])) { + // if (IsIdentifierActive(data[i])) { PreaccHelper.addInput(data[i]); - } + // } } } } @@ -539,9 +548,9 @@ FORCEINLINE void SetPreaccIn(const T& data, const int size_x, const int size_y) if (!PreaccActive) return; for (int i = 0; i < size_x; i++) { for (int j = 0; j < size_y; j++) { - if (IsIdentifierActive(data[i][j])) { + // if (IsIdentifierActive(data[i][j])) { PreaccHelper.addInput(data[i][j]); - } + // } } } } @@ -559,7 +568,8 @@ FORCEINLINE void SetPreaccOut() {} template ::value> = 0> FORCEINLINE void SetPreaccOut(T& data, Ts&&... moreData) { if (!PreaccActive) return; - if (IsIdentifierActive(data)) PreaccHelper.addOutput(data); + // if (IsIdentifierActive(data)) + PreaccHelper.addOutput(data); SetPreaccOut(moreData...); } @@ -567,9 +577,9 @@ template FORCEINLINE void SetPreaccOut(T&& data, const int size) { if (PreaccActive) { for (int i = 0; i < size; i++) { - if (IsIdentifierActive(data[i])) { + // if (IsIdentifierActive(data[i])) { PreaccHelper.addOutput(data[i]); - } + // } } } } @@ -579,9 +589,9 @@ FORCEINLINE void SetPreaccOut(T&& data, const int size_x, const int size_y) { if (!PreaccActive) return; for (int i = 0; i < size_x; i++) { for (int j = 0; j < size_y; j++) { - if (IsIdentifierActive(data[i][j])) { + // if (IsIdentifierActive(data[i][j])) { PreaccHelper.addOutput(data[i][j]); - } + // } } } } diff --git a/Common/include/option_structure.hpp b/Common/include/option_structure.hpp index d81c49448aa..d456a587460 100644 --- a/Common/include/option_structure.hpp +++ b/Common/include/option_structure.hpp @@ -409,6 +409,7 @@ enum ENUM_TRANSFER { CONJUGATE_HEAT_WEAKLY_FS = 17, /*!< \brief Conjugate heat transfer (between incompressible fluids and solids). */ CONJUGATE_HEAT_SF = 18, /*!< \brief Conjugate heat transfer (between solids and compressible fluids). */ CONJUGATE_HEAT_WEAKLY_SF = 19, /*!< \brief Conjugate heat transfer (between solids and incompressible fluids). */ + CONJUGATE_HEAT_SS = 20, /*!< \brief Conjugate heat transfer (between two solids). */ }; /*! @@ -1621,6 +1622,7 @@ enum BC_TYPE { FLUID_INTERFACE = 39, /*!< \brief Domain interface definition. */ DISP_DIR_BOUNDARY = 40, /*!< \brief Boundary displacement definition. */ DAMPER_BOUNDARY = 41, /*!< \brief Damper. */ + MIXING_PLANE_INTERFACE = 42, /*< \breif Mxing plane */ CHT_WALL_INTERFACE = 50, /*!< \brief Domain interface definition. */ SMOLUCHOWSKI_MAXWELL = 55, /*!< \brief Smoluchoski/Maxwell wall boundary condition. */ SEND_RECEIVE = 99, /*!< \brief Boundary send-receive definition. */ @@ -1751,7 +1753,8 @@ enum RIEMANN_TYPE { TOTAL_CONDITIONS_PT_1D = 11, STATIC_PRESSURE_1D = 12, MIXING_IN_1D = 13, - MIXING_OUT_1D =14 + MIXING_OUT_1D = 14, + MASS_FLOW_OUTLET = 15 }; static const MapType Riemann_Map = { MakePair("TOTAL_CONDITIONS_PT", TOTAL_CONDITIONS_PT) @@ -1768,6 +1771,7 @@ static const MapType Riemann_Map = { MakePair("RADIAL_EQUILIBRIUM", RADIAL_EQUILIBRIUM) MakePair("TOTAL_CONDITIONS_PT_1D", TOTAL_CONDITIONS_PT_1D) MakePair("STATIC_PRESSURE_1D", STATIC_PRESSURE_1D) + MakePair("MASS_FLOW_OUTLET", MASS_FLOW_OUTLET) }; static const MapType Giles_Map = { @@ -1785,6 +1789,7 @@ static const MapType Giles_Map = { MakePair("RADIAL_EQUILIBRIUM", RADIAL_EQUILIBRIUM) MakePair("TOTAL_CONDITIONS_PT_1D", TOTAL_CONDITIONS_PT_1D) MakePair("STATIC_PRESSURE_1D", STATIC_PRESSURE_1D) + MakePair("MASS_FLOW_OUTLET", MASS_FLOW_OUTLET) }; /*! @@ -1861,6 +1866,18 @@ static const MapType TurboPerfKind_Map = { MakePair("PROPELLOR", TURBO_PERF_KIND::PROPELLOR) }; +/*! + * \brief Types of Turbomachinery interfaces. + */ +enum class TURBO_INTERFACE_KIND{ + MIXING_PLANE = ENUM_TRANSFER::MIXING_PLANE, + FROZEN_ROTOR = ENUM_TRANSFER::SLIDING_INTERFACE, +}; +static const MapType TurboInterfaceKind_Map = { + MakePair("MIXING_PLANE", TURBO_INTERFACE_KIND::MIXING_PLANE) + MakePair("FROZEN_ROTOR", TURBO_INTERFACE_KIND::FROZEN_ROTOR) +}; + /*! * \brief Types of Turbomachinery performance flag. */ @@ -2482,6 +2499,8 @@ enum class RECORDING { MESH_COORDS, MESH_DEFORM, SOLUTION_AND_MESH, + TAG_INIT_SOLUTION_VARIABLES, + TAG_CHECK_SOLUTION_VARIABLES }; /*! diff --git a/Common/include/option_structure.inl b/Common/include/option_structure.inl index 955d1531482..01ed6de9bf9 100644 --- a/Common/include/option_structure.inl +++ b/Common/include/option_structure.inl @@ -1606,11 +1606,15 @@ class COptionTurboPerformance : public COptionBase { unsigned short& size; string*& marker_turboIn; string*& marker_turboOut; + string*& markers; public: COptionTurboPerformance(const string option_field_name, unsigned short& nMarker_TurboPerf, - string*& Marker_TurboBoundIn, string*& Marker_TurboBoundOut) - : size(nMarker_TurboPerf), marker_turboIn(Marker_TurboBoundIn), marker_turboOut(Marker_TurboBoundOut) { + string*& Marker_TurboBoundIn, string*& Marker_TurboBoundOut, string*& Marker_Turbomachinery) + : size(nMarker_TurboPerf), + marker_turboIn(Marker_TurboBoundIn), + marker_turboOut(Marker_TurboBoundOut), + markers(Marker_Turbomachinery) { this->name = option_field_name; } @@ -1624,6 +1628,7 @@ class COptionTurboPerformance : public COptionBase { this->size = 0; this->marker_turboIn = nullptr; this->marker_turboOut = nullptr; + this->markers = nullptr; return ""; } @@ -1634,10 +1639,16 @@ class COptionTurboPerformance : public COptionBase { this->size = 0; this->marker_turboIn = nullptr; this->marker_turboOut = nullptr; + this->markers = nullptr; ; return newstring; } + this->markers = new string[totalVals]; + for (unsigned long i = 0; i < totalVals; i++) { + this->markers[i].assign(option_value[i]); + } + unsigned long nVals = totalVals / mod_num; this->size = nVals; this->marker_turboIn = new string[nVals]; @@ -1654,6 +1665,7 @@ class COptionTurboPerformance : public COptionBase { this->size = 0; this->marker_turboIn = nullptr; this->marker_turboOut = nullptr; + this->markers = nullptr; } }; diff --git a/Common/src/CConfig.cpp b/Common/src/CConfig.cpp index c40a80d40f9..be45861f59a 100644 --- a/Common/src/CConfig.cpp +++ b/Common/src/CConfig.cpp @@ -539,10 +539,10 @@ void CConfig::addPeriodicOption(const string & name, unsigned short & nMarker_Pe } void CConfig::addTurboPerfOption(const string & name, unsigned short & nMarker_TurboPerf, - string* & Marker_TurboBoundIn, string* & Marker_TurboBoundOut) { + string* & Marker_TurboBoundIn, string* & Marker_TurboBoundOut, string* & Marker_Turbomachinery) { assert(option_map.find(name) == option_map.end()); all_options.insert(pair(name, true)); - COptionBase* val = new COptionTurboPerformance(name, nMarker_TurboPerf, Marker_TurboBoundIn, Marker_TurboBoundOut); + COptionBase* val = new COptionTurboPerformance(name, nMarker_TurboPerf, Marker_TurboBoundIn, Marker_TurboBoundOut, Marker_Turbomachinery); option_map.insert(pair(name, val)); } @@ -1037,6 +1037,7 @@ void CConfig::SetPointersNull() { Marker_MixingPlaneInterface = nullptr; Marker_TurboBoundIn = nullptr; Marker_TurboBoundOut = nullptr; + Marker_Turbomachinery = nullptr; Marker_Giles = nullptr; Marker_Shroud = nullptr; @@ -1109,11 +1110,19 @@ void CConfig::SetConfig_Options() { addBoolOption("MULTIZONE", Multizone_Problem, NO); /*!\brief PHYSICAL_PROBLEM \n DESCRIPTION: Physical governing equations \n Options: see \link Solver_Map \endlink \n DEFAULT: NONE \ingroup Config*/ addEnumOption("MULTIZONE_SOLVER", Kind_MZSolver, Multizone_Map, ENUM_MULTIZONE::MZ_BLOCK_GAUSS_SEIDEL); -#ifdef CODI_REVERSE_TYPE +#if defined (CODI_REVERSE_TYPE) const bool discAdjDefault = true; +# if defined (CODI_TAG_TAPE) + const bool discAdjDebugDefault = true; +# else + const bool discAdjDebugDefault = false; +# endif #else const bool discAdjDefault = false; + const bool discAdjDebugDefault = false; #endif + // TODO Set the indicator through MATH_PROBLEM + DiscreteAdjointDebug = discAdjDebugDefault; /*!\brief MATH_PROBLEM \n DESCRIPTION: Mathematical problem \n Options: DIRECT, ADJOINT \ingroup Config*/ addMathProblemOption("MATH_PROBLEM", ContinuousAdjoint, false, DiscreteAdjoint, discAdjDefault, Restart_Flow, discAdjDefault); /*!\brief KIND_TURB_MODEL \n DESCRIPTION: Specify turbulence model \n Options: see \link Turb_Model_Map \endlink \n DEFAULT: NONE \ingroup Config*/ @@ -1522,6 +1531,8 @@ void CConfig::SetConfig_Options() { addStringListOption("MARKER_ZONE_INTERFACE", nMarker_ZoneInterface, Marker_ZoneInterface); /*!\brief MARKER_CHT_INTERFACE \n DESCRIPTION: CHT interface boundary marker(s) \ingroup Config*/ addStringListOption("MARKER_CHT_INTERFACE", nMarker_CHTInterface, Marker_CHTInterface); + /*!\brief CHT_INTERFACE_CONTACT_RESISTANCE: Thermal contact resistance values for each CHT inerface. \ingroup Config*/ + addDoubleListOption("CHT_INTERFACE_CONTACT_RESISTANCE", nMarker_ContactResistance, CHT_ContactResistance); /* DESCRIPTION: Internal boundary marker(s) */ addStringListOption("MARKER_INTERNAL", nMarker_Internal, Marker_Internal); /* DESCRIPTION: Custom boundary marker(s) */ @@ -1636,8 +1647,8 @@ void CConfig::SetConfig_Options() { addStringListOption("MARKER_MIXINGPLANE_INTERFACE", nMarker_MixingPlaneInterface, Marker_MixingPlaneInterface); /*!\brief TURBULENT_MIXINGPLANE \n DESCRIPTION: Activate mixing plane also for turbulent quantities \ingroup Config*/ addBoolOption("TURBULENT_MIXINGPLANE", turbMixingPlane, false); - /*!\brief MARKER_TURBOMACHINERY \n DESCRIPTION: Identify the inflow and outflow boundaries in which the turbomachinery settings are applied. \ingroup Config*/ - addTurboPerfOption("MARKER_TURBOMACHINERY", nMarker_Turbomachinery, Marker_TurboBoundIn, Marker_TurboBoundOut); + /*!\brief MARKER_TURBOMACHINERY \n DESCRIPTION: Identify the boundaries for which the turbomachinery settings are applied. \ingroup Config*/ + addTurboPerfOption("MARKER_TURBOMACHINERY", nMarker_Turbomachinery, Marker_TurboBoundIn, Marker_TurboBoundOut, Marker_Turbomachinery); /*!\brief NUM_SPANWISE_SECTIONS \n DESCRIPTION: Integer number of spanwise sections to compute 3D turbo BC and Performance for turbomachinery */ addUnsignedShortOption("NUM_SPANWISE_SECTIONS", nSpanWiseSections_User, 1); /*!\brief SPANWISE_KIND \n DESCRIPTION: type of algorithm to identify the span-wise sections at the turbo boundaries. @@ -1739,6 +1750,8 @@ void CConfig::SetConfig_Options() { /*!\brief RAMP_AND_RELEASE\n DESCRIPTION: release the load after applying the ramp*/ addBoolOption("RAMP_AND_RELEASE_LOAD", RampAndRelease, false); + /* DESCRIPTION: Evaluation frequency for Engine and Actuator disk markers. */ + addUnsignedLongOption("BC_EVAL_FREQ", Bc_Eval_Freq, 40); /* DESCRIPTION: Damping factor for engine inlet condition */ addDoubleOption("DAMP_ENGINE_INFLOW", Damp_Engine_Inflow, 0.95); /* DESCRIPTION: Damping factor for engine exhaust condition */ @@ -3574,6 +3587,25 @@ void CConfig::SetPostprocessing(SU2_COMPONENT val_software, unsigned short val_i SU2_MPI::Error(string("You probably want to set INC_ENERGY_EQUATION= YES for the fluid solver. \n"), CURRENT_FUNCTION); } + /*--- Check correctness and consistency of contact resistance options. ---*/ + if (nMarker_ContactResistance > 0) { + + /*--- Set constant contact resistance across CHT interfaces if a single value is provided. ---*/ + if (nMarker_ContactResistance == 1) { + auto val_CHTInterface = CHT_ContactResistance[0]; + delete [] CHT_ContactResistance; + CHT_ContactResistance = new su2double[nMarker_CHTInterface]; + for (auto iCHTMarker=0u; iCHTMarker < nMarker_CHTInterface; iCHTMarker++) + CHT_ContactResistance[iCHTMarker] = val_CHTInterface; + }else if((nMarker_CHTInterface/2) != nMarker_ContactResistance){ + SU2_MPI::Error("Number of CHT interfaces does not match number of contact resistances.", CURRENT_FUNCTION); + } + for (auto iCHTMarker=0u; iCHTMarker < nMarker_ContactResistance; iCHTMarker++){ + if (CHT_ContactResistance[iCHTMarker] < 0) + SU2_MPI::Error("Contact resistance value should be positive.", CURRENT_FUNCTION); + } + } + /*--- By default, in 2D we should use TWOD_AIRFOIL (independenly from the input file) ---*/ if (val_nDim == 2) Geo_Description = TWOD_AIRFOIL; @@ -5654,6 +5686,7 @@ void CConfig::SetMarkers(SU2_COMPONENT val_software) { Marker_All_Turbomachinery = new unsigned short[nMarker_All] (); // Store whether the boundary is in needed for Turbomachinery computations. Marker_All_TurbomachineryFlag = new unsigned short[nMarker_All] (); // Store whether the boundary has a flag for Turbomachinery computations. Marker_All_MixingPlaneInterface = new unsigned short[nMarker_All] (); // Store whether the boundary has a in the MixingPlane interface. + Marker_All_Giles = new unsigned short[nMarker_All] (); // Store whether the boundary has is a Giles boundary. Marker_All_SobolevBC = new unsigned short[nMarker_All] (); // Store wether the boundary should apply to the gradient smoothing. for (iMarker_All = 0; iMarker_All < nMarker_All; iMarker_All++) { @@ -5679,6 +5712,7 @@ void CConfig::SetMarkers(SU2_COMPONENT val_software) { Marker_CfgFile_Turbomachinery = new unsigned short[nMarker_CfgFile] (); Marker_CfgFile_TurbomachineryFlag = new unsigned short[nMarker_CfgFile] (); Marker_CfgFile_MixingPlaneInterface = new unsigned short[nMarker_CfgFile] (); + Marker_CfgFile_Giles = new unsigned short[nMarker_CfgFile] (); Marker_CfgFile_PyCustom = new unsigned short[nMarker_CfgFile] (); Marker_CfgFile_SobolevBC = new unsigned short[nMarker_CfgFile] (); @@ -5799,7 +5833,7 @@ void CConfig::SetMarkers(SU2_COMPONENT val_software) { for (iMarker_Fluid_InterfaceBound = 0; iMarker_Fluid_InterfaceBound < nMarker_Fluid_InterfaceBound; iMarker_Fluid_InterfaceBound++) { Marker_CfgFile_TagBound[iMarker_CfgFile] = Marker_Fluid_InterfaceBound[iMarker_Fluid_InterfaceBound]; - Marker_CfgFile_KindBC[iMarker_CfgFile] = FLUID_INTERFACE; + Marker_CfgFile_KindBC[iMarker_CfgFile] = BC_TYPE::FLUID_INTERFACE; iMarker_CfgFile++; } @@ -6030,6 +6064,17 @@ void CConfig::SetMarkers(SU2_COMPONENT val_software) { } } + /*--- Idenftification fo Giles Markers ---*/ + // This is seperate from MP and Turbomachinery Markers as all mixing plane markers are Giles, + // but not all Giles markers are mixing plane + for (iMarker_CfgFile = 0; iMarker_CfgFile < nMarker_CfgFile; iMarker_CfgFile++) { + Marker_CfgFile_Giles[iMarker_CfgFile] = NO; + for (iMarker_Giles = 0; iMarker_Giles < nMarker_Giles; iMarker_Giles++) { + if (Marker_CfgFile_TagBound[iMarker_CfgFile] == Marker_Giles[iMarker_Giles]) + Marker_CfgFile_Giles[iMarker_CfgFile] = YES; + } + } + /*--- Identification of MixingPlane interface markers ---*/ for (iMarker_CfgFile = 0; iMarker_CfgFile < nMarker_CfgFile; iMarker_CfgFile++) { @@ -6041,6 +6086,37 @@ void CConfig::SetMarkers(SU2_COMPONENT val_software) { Marker_CfgFile_MixingPlaneInterface[iMarker_CfgFile] = indexMarker; } + /*--- Once we have identified the MixingPlane and Turbomachinery markers + * we next need to determine what type of interface between zones is + * used. It is convenient to do this here as it tidies up the interface + * preproccesing in CDriver ---*/ + if (nMarker_Turbomachinery != 0) { + nTurboInterfaces = (nMarker_Turbomachinery - 1)*2; //Two markers per zone minus inlet & outlet + Kind_TurboInterface.resize(nTurboInterfaces); + /*--- Loop over all markers ---*/ + for (iMarker_CfgFile = 0; iMarker_CfgFile < nMarker_CfgFile; iMarker_CfgFile++) { + /*--- Identify mixing plane markers ---*/ + if (Marker_MixingPlaneInterface != nullptr){ // Necessary in cases where no mixing plane interfaces are defined + if (Marker_CfgFile_MixingPlaneInterface[iMarker_CfgFile] != 0) { //Is a mixing plane + /*--- Find which list position this marker is in turbomachinery markers ---*/ + const auto* target = std::find(Marker_Turbomachinery, &Marker_Turbomachinery[nMarker_Turbomachinery*2-1], Marker_CfgFile_TagBound[iMarker_CfgFile]); + const auto target_index = target - Marker_Turbomachinery; + /*--- Assert that we find the marker within the turbomachienry markers ---*/ + assert(target != &Marker_Turbomachinery[nMarker_Turbomachinery*2-1]); + /*--- Assign the correct interface ---*/ + Kind_TurboInterface[target_index-1] = TURBO_INTERFACE_KIND::MIXING_PLANE; // Need to subtract 1 from index as to not consider the inlet an interface + } + } + if (Marker_Fluid_InterfaceBound != nullptr){ // No fluid interfaces are defined in the config file (nullptr if no interfaces defined) + if (Marker_CfgFile_KindBC[iMarker_CfgFile] == BC_TYPE::FLUID_INTERFACE) { // iMarker_CfgFile is a fluid interface + const auto* target = std::find(Marker_Turbomachinery, &Marker_Turbomachinery[nMarker_Turbomachinery*2-1], Marker_CfgFile_TagBound[iMarker_CfgFile]); + const auto target_index = target - Marker_Turbomachinery; + Kind_TurboInterface[target_index-1] = TURBO_INTERFACE_KIND::FROZEN_ROTOR; + } + } + } + } + for (iMarker_CfgFile = 0; iMarker_CfgFile < nMarker_CfgFile; iMarker_CfgFile++) { Marker_CfgFile_DV[iMarker_CfgFile] = NO; for (iMarker_DV = 0; iMarker_DV < nMarker_DV; iMarker_DV++) @@ -7934,6 +8010,13 @@ unsigned short CConfig::GetMarker_CfgFile_MixingPlaneInterface(const string& val return Marker_CfgFile_MixingPlaneInterface[iMarker_CfgFile]; } +unsigned short CConfig::GetMarker_CfgFile_Giles(const string& val_marker) const { + unsigned short iMarker_CfgFile; + for (iMarker_CfgFile = 0; iMarker_CfgFile < nMarker_CfgFile; iMarker_CfgFile++) + if (Marker_CfgFile_TagBound[iMarker_CfgFile] == val_marker) break; + return Marker_CfgFile_Giles[iMarker_CfgFile]; +} + unsigned short CConfig::GetMarker_CfgFile_DV(const string& val_marker) const { unsigned short iMarker_CfgFile; for (iMarker_CfgFile = 0; iMarker_CfgFile < nMarker_CfgFile; iMarker_CfgFile++) @@ -8121,6 +8204,9 @@ CConfig::~CConfig() { delete[] Marker_CfgFile_TurbomachineryFlag; delete[] Marker_All_TurbomachineryFlag; + delete[] Marker_CfgFile_Giles; + delete[] Marker_All_Giles; + delete[] Marker_CfgFile_MixingPlaneInterface; delete[] Marker_All_MixingPlaneInterface; @@ -8261,6 +8347,7 @@ CConfig::~CConfig() { delete [] Marker_TurboBoundIn; delete [] Marker_TurboBoundOut; + delete [] Marker_Turbomachinery; delete [] Marker_Riemann; delete [] Marker_Giles; diff --git a/Common/src/geometry/CPhysicalGeometry.cpp b/Common/src/geometry/CPhysicalGeometry.cpp index a069424c666..b9b0a0f36cb 100644 --- a/Common/src/geometry/CPhysicalGeometry.cpp +++ b/Common/src/geometry/CPhysicalGeometry.cpp @@ -3370,6 +3370,7 @@ void CPhysicalGeometry::SetBoundaries(CConfig* config) { config->SetMarker_All_PerBound(iMarker, config->GetMarker_CfgFile_PerBound(Marker_Tag)); config->SetMarker_All_Turbomachinery(iMarker, config->GetMarker_CfgFile_Turbomachinery(Marker_Tag)); config->SetMarker_All_TurbomachineryFlag(iMarker, config->GetMarker_CfgFile_TurbomachineryFlag(Marker_Tag)); + config->SetMarker_All_Giles(iMarker, config->GetMarker_CfgFile_Giles(Marker_Tag)); config->SetMarker_All_MixingPlaneInterface(iMarker, config->GetMarker_CfgFile_MixingPlaneInterface(Marker_Tag)); config->SetMarker_All_SobolevBC(iMarker, config->GetMarker_CfgFile_SobolevBC(Marker_Tag)); @@ -3394,6 +3395,7 @@ void CPhysicalGeometry::SetBoundaries(CConfig* config) { config->SetMarker_All_PerBound(iMarker, NO); config->SetMarker_All_Turbomachinery(iMarker, NO); config->SetMarker_All_TurbomachineryFlag(iMarker, NO); + config->SetMarker_All_Giles(iMarker, NO); config->SetMarker_All_MixingPlaneInterface(iMarker, NO); config->SetMarker_All_SobolevBC(iMarker, NO); @@ -3764,6 +3766,7 @@ void CPhysicalGeometry::LoadUnpartitionedSurfaceElements(CConfig* config, CMeshR config->SetMarker_All_SendRecv(iMarker, NONE); config->SetMarker_All_Turbomachinery(iMarker, config->GetMarker_CfgFile_Turbomachinery(Marker_Tag)); config->SetMarker_All_TurbomachineryFlag(iMarker, config->GetMarker_CfgFile_TurbomachineryFlag(Marker_Tag)); + config->SetMarker_All_Giles(iMarker, config->GetMarker_CfgFile_Giles(Marker_Tag)); config->SetMarker_All_MixingPlaneInterface(iMarker, config->GetMarker_CfgFile_MixingPlaneInterface(Marker_Tag)); config->SetMarker_All_SobolevBC(iMarker, config->GetMarker_CfgFile_SobolevBC(Marker_Tag)); } diff --git a/SU2_CFD/include/drivers/CDiscAdjMultizoneDriver.hpp b/SU2_CFD/include/drivers/CDiscAdjMultizoneDriver.hpp index f7c8929c805..1f5ad0f21f5 100644 --- a/SU2_CFD/include/drivers/CDiscAdjMultizoneDriver.hpp +++ b/SU2_CFD/include/drivers/CDiscAdjMultizoneDriver.hpp @@ -69,7 +69,7 @@ class CDiscAdjMultizoneDriver : public CMultizoneDriver { }; /*! - * \brief Kinds of recordings (three different ones). + * \brief Kinds of recordings. */ enum class Kind_Tape { FULL_TAPE, /*!< \brief Entire derivative information for a coupled adjoint @@ -141,6 +141,11 @@ class CDiscAdjMultizoneDriver : public CMultizoneDriver { */ void StartSolver() override; + /*! + * \brief [Overload] Launch the debug mode for the discrete adjoint multizone solver. + */ + void DebugRun() override; + /*! * \brief Preprocess the multizone iteration */ diff --git a/SU2_CFD/include/drivers/CDriver.hpp b/SU2_CFD/include/drivers/CDriver.hpp index c5e78bd9037..855f5955f2d 100644 --- a/SU2_CFD/include/drivers/CDriver.hpp +++ b/SU2_CFD/include/drivers/CDriver.hpp @@ -425,6 +425,11 @@ class CDriver : public CDriverBase { */ virtual void StartSolver() {} + /*! + * \brief Launch a debug run of the solver. + */ + virtual void DebugRun() {} + /*! * \brief Deallocation routine */ diff --git a/SU2_CFD/include/fluid/CFluidScalar.hpp b/SU2_CFD/include/fluid/CFluidScalar.hpp index 9f5c9108c51..e3295d99b6d 100644 --- a/SU2_CFD/include/fluid/CFluidScalar.hpp +++ b/SU2_CFD/include/fluid/CFluidScalar.hpp @@ -41,7 +41,6 @@ class CFluidScalar final : public CFluidModel { private: const int n_species_mixture; /*!< \brief Number of species in mixture. */ su2double Gas_Constant; /*!< \brief Specific gas constant. */ - const su2double Gamma; /*!< \brief Ratio of specific heats of the gas. */ const su2double Pressure_Thermodynamic; /*!< \brief Constant pressure thermodynamic. */ const su2double GasConstant_Ref; /*!< \brief Gas constant reference needed for Nondimensional problems. */ const su2double Prandtl_Number; /*!< \brief Prandlt number.*/ @@ -106,7 +105,7 @@ class CFluidScalar final : public CFluidModel { /*! * \brief Constructor of the class. */ - CFluidScalar(su2double val_Cp, su2double val_gas_constant, su2double val_operating_pressure, const CConfig* config); + CFluidScalar(su2double val_operating_pressure, const CConfig* config); /*! * \brief Set viscosity model. diff --git a/SU2_CFD/include/interfaces/CInterface.hpp b/SU2_CFD/include/interfaces/CInterface.hpp index fcc4600b5aa..a3ffc31c0fb 100644 --- a/SU2_CFD/include/interfaces/CInterface.hpp +++ b/SU2_CFD/include/interfaces/CInterface.hpp @@ -159,15 +159,6 @@ class CInterface { const CConfig *target_config, unsigned long Marker_Target, unsigned long Vertex_Target, unsigned long Point_Target) = 0; - // /*! - // * \brief A virtual member. - // * \param[in] target_solution - Solution from the target mesh. - // * \param[in] target_solution - Solution from the target mesh. - // * \param[in] donor_zone - Index of the donorZone. - // */ - // inline virtual void SetAverageValues(CSolver *donor_solution, CSolver *target_solution, - // unsigned short donorZone) { } - /*! * \brief A virtual member. * \param[in] donor_geometry - Geometry of the target mesh. @@ -216,5 +207,21 @@ class CInterface { void AllgatherAverage(CSolver *donor_solution, CSolver *target_solution, CGeometry *donor_geometry, CGeometry *target_geometry, const CConfig *donor_config, const CConfig *target_config, unsigned short iMarkerInt); - -}; \ No newline at end of file + + /*! + * \brief Interpolate data and scatter it into different processors, for matching meshes. + * \param[in] donor_solution - Solution from the donor mesh. + * \param[in] target_solution - Solution from the target mesh. + * \param[in] donor_geometry - Geometry of the donor mesh. + * \param[in] target_geometry - Geometry of the target mesh. + * \param[in] donor_config - Definition of the problem at the donor mesh. + * \param[in] target_config - Definition of the problem at the target mesh. + */ + void GatherAverageValues(CSolver *donor_solution, CSolver *target_solution, unsigned short donorZone); + + /*! + * \brief Set the contact resistance value for the solid-to-solid heat transfer interface. + * \param[in] val_contact_resistance - Contact resistance value in m^2/W + */ + inline virtual void SetContactResistance(su2double val_contact_resistance) {}; +}; diff --git a/SU2_CFD/include/interfaces/cht/CConjugateHeatInterface.hpp b/SU2_CFD/include/interfaces/cht/CConjugateHeatInterface.hpp index 919d4e8f7f1..982733d5e36 100644 --- a/SU2_CFD/include/interfaces/cht/CConjugateHeatInterface.hpp +++ b/SU2_CFD/include/interfaces/cht/CConjugateHeatInterface.hpp @@ -35,6 +35,7 @@ * \ingroup Interfaces */ class CConjugateHeatInterface : public CInterface { + su2double ContactResistance = 0; /*!<\brief Contact resistance value of the current inerface. */ public: /*! * \brief Constructor of the class. @@ -70,4 +71,10 @@ class CConjugateHeatInterface : public CInterface { */ void SetTarget_Variable(CSolver *target_solution, CGeometry *target_geometry, const CConfig *target_config, unsigned long Marker_Target, unsigned long Vertex_Target, unsigned long Point_Target) override; + + /*! + * \brief Set the contact resistance value for the solid-to-solid heat transfer interface. + * \param[in] val_contact_resistance - Contact resistance value in m^2/W + */ + void SetContactResistance(su2double val_contact_resistance) override { ContactResistance = val_contact_resistance; } }; diff --git a/SU2_CFD/include/numerics/turbulent/turb_sources.hpp b/SU2_CFD/include/numerics/turbulent/turb_sources.hpp index ce75b38f115..f59da30719b 100644 --- a/SU2_CFD/include/numerics/turbulent/turb_sources.hpp +++ b/SU2_CFD/include/numerics/turbulent/turb_sources.hpp @@ -48,11 +48,11 @@ struct CSAVariables { const su2double cb2_sigma = cb2 / sigma; const su2double cw1 = cb1 / k2 + (1 + cb2) / sigma; const su2double cr1 = 0.5; - const su2double CRot = 1.0; + const su2double CRot = 2.0; const su2double c2 = 0.7, c3 = 0.9; /*--- List of auxiliary functions ---*/ - su2double ft2, d_ft2, r, d_r, g, d_g, glim, fw, d_fw, Ji, d_Ji, S, Shat, d_Shat, fv1, d_fv1, fv2, d_fv2; + su2double ft2, d_ft2, r, d_r, g, d_g, glim, fw, d_fw, Ji, d_Ji, Shat, d_Shat, fv1, d_fv1, fv2, d_fv2, Prod; /*--- List of helpers ---*/ su2double Omega, dist_i_2, inv_k2_d2, inv_Shat, g_6, norm2_Grad; @@ -151,20 +151,7 @@ class CSourceBase_TurbSA : public CNumerics { Residual = 0.0; Jacobian_i[0] = 0.0; - /*--- Evaluate Omega with a rotational correction term. ---*/ - - Omega::get(Vorticity_i, nDim, PrimVar_Grad_i + idx.Velocity(), var); - - /*--- Dacles-Mariani et. al. rotation correction ("-R"). ---*/ - if (options.rot) { - var.Omega += var.CRot * min(0.0, StrainMag_i - var.Omega); - /*--- Do not allow negative production for SA-neg. ---*/ - if (ScalarVar_i[0] < 0) var.Omega = abs(var.Omega); - } - if (dist_i > 1e-10) { - /*--- Vorticity ---*/ - var.S = var.Omega; var.dist_i_2 = pow(dist_i, 2); const su2double nu = laminar_viscosity / density; @@ -188,12 +175,24 @@ class CSourceBase_TurbSA : public CNumerics { var.fv2 = 1 - ScalarVar_i[0] / (nu + ScalarVar_i[0] * var.fv1); var.d_fv2 = -(1 / nu - Ji_2 * var.d_fv1) / pow(1 + var.Ji * var.fv1, 2); - /*--- Compute ft2 term ---*/ - ft2::get(var); + /*--- Evaluate Omega with a rotational correction term. ---*/ + + Omega::get(Vorticity_i, nDim, PrimVar_Grad_i + idx.Velocity(), var); /*--- Compute modified vorticity ---*/ ModVort::get(ScalarVar_i[0], nu, var); var.inv_Shat = 1.0 / var.Shat; + var.Prod = var.Shat; + + /*--- Dacles-Mariani et. al. rotation correction ("-R"). ---*/ + if (options.rot) { + var.Prod += var.CRot * min(0.0, StrainMag_i - var.Omega); + /*--- Do not allow negative production for SA-neg. ---*/ + if (ScalarVar_i[0] < 0) var.Prod = abs(var.Prod); + } + + /*--- Compute ft2 term ---*/ + ft2::get(var); /*--- Compute auxiliary function r ---*/ rFunc::get(ScalarVar_i[0], var); @@ -234,7 +233,6 @@ class CSourceBase_TurbSA : public CNumerics { } else if (transition_LM){ var.intermittency = intermittency_eff_i; - //var.intermittency = 1.0; // Is wrong the reference from NASA? // Original max(min(gamma, 0.5), 1.0) always gives 1 as result. var.interDestrFactor = min(max(intermittency_i, 0.5), 1.0); @@ -347,12 +345,12 @@ struct Bsl { /*--- Limiting of \hat{S} based on "Modifications and Clarifications for the Implementation of the Spalart-Allmaras Turbulence Model" * Note 1 option c in https://turbmodels.larc.nasa.gov/spalart.html ---*/ - if (Sbar >= - c2 * var.S) { - var.Shat = var.S + Sbar; + if (Sbar >= - c2 * var.Omega) { + var.Shat = var.Omega + Sbar; } else { - const su2double Num = var.S * (c2 * c2 * var.S + c3 * Sbar); - const su2double Den = (c3 - 2 * c2) * var.S - Sbar; - var.Shat = var.S + Num / Den; + const su2double Num = var.Omega * (c2 * c2 * var.Omega + c3 * Sbar); + const su2double Den = (c3 - 2 * c2) * var.Omega - Sbar; + var.Shat = var.Omega + Num / Den; } if (var.Shat <= 1e-10) { var.Shat = 1e-10; @@ -366,12 +364,12 @@ struct Bsl { /*! \brief Edward. */ struct Edw { static void get(const su2double& nue, const su2double& nu, CSAVariables& var) { - var.Shat = max(var.S * ((1.0 / max(var.Ji, 1.0e-16)) + var.fv1), 1.0e-16); + var.Shat = max(var.Omega * ((1.0 / max(var.Ji, 1.0e-16)) + var.fv1), 1.0e-16); var.Shat = max(var.Shat, 1.0e-10); if (var.Shat <= 1.0e-10) { var.d_Shat = 0.0; } else { - var.d_Shat = -var.S * pow(var.Ji, -2) / nu + var.S * var.d_fv1; + var.d_Shat = -var.Omega * pow(var.Ji, -2) / nu + var.Omega * var.d_fv1; } } }; @@ -383,7 +381,7 @@ struct Neg { // Baseline solution Bsl::get(nue, nu, var); } else { - var.Shat = 1.0e-10; + var.Shat = var.Omega; var.d_Shat = 0.0; } /*--- Don't check whether Sbar <>= -cv2*S. @@ -447,8 +445,8 @@ struct Bsl { static void ComputeProduction(const su2double& nue, const CSAVariables& var, su2double& production, su2double& jacobian) { const su2double factor = var.intermittency * var.cb1; - production = factor * (1.0 - var.ft2) * var.Shat * nue; - jacobian += factor * (-var.Shat * nue * var.d_ft2 + (1.0 - var.ft2) * (nue * var.d_Shat + var.Shat)); + production = factor * (1.0 - var.ft2) * var.Prod * nue; + jacobian += factor * (-var.Prod * nue * var.d_ft2 + (1.0 - var.ft2) * (nue * var.d_Shat + var.Prod)); } static void ComputeDestruction(const su2double& nue, const CSAVariables& var, su2double& destruction, @@ -481,7 +479,7 @@ struct Neg { static void ComputeProduction(const su2double& nue, const CSAVariables& var, su2double& production, su2double& jacobian) { - const su2double dP_dnu = var.intermittency * var.cb1 * (1.0 - var.ct3) * var.S; + const su2double dP_dnu = var.intermittency * var.cb1 * (1.0 - var.ct3) * var.Prod; production = dP_dnu * nue; jacobian += dP_dnu; } diff --git a/SU2_CFD/include/solvers/CEulerSolver.hpp b/SU2_CFD/include/solvers/CEulerSolver.hpp index 342448c6253..d498b84c225 100644 --- a/SU2_CFD/include/solvers/CEulerSolver.hpp +++ b/SU2_CFD/include/solvers/CEulerSolver.hpp @@ -138,6 +138,7 @@ class CEulerSolver : public CFVMFlowSolverBase::ComputeVorticityAndStrainMag(const CConfig& confi StrainMag(iPoint) = sqrt(2.0*StrainMag(iPoint)); AD::SetPreaccOut(StrainMag(iPoint)); + AD::EndPreacc(); + /*--- Max is not differentiable, so we not register them for preacc. ---*/ strainMax = max(strainMax, StrainMag(iPoint)); omegaMax = max(omegaMax, GeometryToolbox::Norm(3, Vorticity)); - - AD::EndPreacc(); } END_SU2_OMP_FOR diff --git a/SU2_CFD/include/solvers/CSolver.hpp b/SU2_CFD/include/solvers/CSolver.hpp index 5eb43dc2035..6c9194f31ae 100644 --- a/SU2_CFD/include/solvers/CSolver.hpp +++ b/SU2_CFD/include/solvers/CSolver.hpp @@ -1089,38 +1089,6 @@ class CSolver { CConfig *config, unsigned short val_marker) { } - /*! - * - * \brief Generalized handling of calculation of total inlet boundary condition inputs - * \param[in] conv_numerics - Description of the numerical method. - * \param[in] config - Definition of the particular problem. - * \param[in] val_marker - Surface marker where the boundary condition is applied. - * \param[in] iSpan - Current spanwise position - * \param[in] SpanwisePosition - Spanwise position where flow quantaties are evaluated - */ - inline virtual void BC_Giles_Total_Inlet(CNumerics *conv_numerics, - CConfig *config, - unsigned short val_marker, - su2double *&c_avg, - su2double *&R, - su2double **&R_c_inv, - su2double **&R_c, - unsigned short iSpan, - unsigned short SpanwisePosition) { } - - /*! - * - * \brief Generalized handling of calculation of mixing plane boundary condition inputs - * \param[in] conv_numerics - Description of the numerical method. - * \param[in] val_marker - Surface marker where the boundary condition is applied. - * \param[in] SpanwisePosition - Spanwise position where flow quantaties are evaluated - */ - inline virtual void BC_Giles_Mixing(CNumerics *conv_numerics, - unsigned short val_marker, - su2double *&deltaprim, - su2double *&c_avg, - unsigned short SpanwisePosition) { } - /*! * \brief A virtual member. * \param[in] geometry - Geometrical definition of the problem. @@ -3831,6 +3799,13 @@ class CSolver { */ inline virtual su2double GetAverageDensity(unsigned short valMarker, unsigned short valSpan) const { return 0.0; } + /*! + * \brief virtual member + * \param[in] val_marker - boundary marker + * \return Value of the mass flow rate on the surface val_marker + */ + inline virtual su2double GetAverageMassFlowRate(unsigned short valMarker) const {return 0.0; } + /*! * \brief A virtual member. * \param[in] val_marker - bound marker. diff --git a/SU2_CFD/src/SU2_CFD.cpp b/SU2_CFD/src/SU2_CFD.cpp index 353a744c891..3bf7ba05a11 100644 --- a/SU2_CFD/src/SU2_CFD.cpp +++ b/SU2_CFD/src/SU2_CFD.cpp @@ -102,6 +102,7 @@ int main(int argc, char *argv[]) { and perform all the preprocessing. ---*/ const bool disc_adj = config.GetDiscrete_Adjoint(); + const bool disc_adj_debug = config.GetDiscrete_Adjoint_Debug(); const bool multizone = config.GetMultizone_Problem(); const bool harmonic_balance = (config.GetTime_Marching() == TIME_MARCHING::HARMONIC_BALANCE); diff --git a/SU2_CFD/src/drivers/CDiscAdjMultizoneDriver.cpp b/SU2_CFD/src/drivers/CDiscAdjMultizoneDriver.cpp index 6a8c51fc0fa..d70619136d2 100644 --- a/SU2_CFD/src/drivers/CDiscAdjMultizoneDriver.cpp +++ b/SU2_CFD/src/drivers/CDiscAdjMultizoneDriver.cpp @@ -159,6 +159,16 @@ void CDiscAdjMultizoneDriver::Preprocess(unsigned long TimeIter) { void CDiscAdjMultizoneDriver::StartSolver() { + /*--- Start the debug recording mode for the discrete adjoint solver. ---*/ + + if (driver_config->GetDiscrete_Adjoint_Debug()) { + + Preprocess(0); + + DebugRun(); + return; + } + const bool time_domain = driver_config->GetTime_Domain(); /*--- Main external loop of the solver. Runs for the number of time steps required. ---*/ @@ -217,6 +227,28 @@ void CDiscAdjMultizoneDriver::StartSolver() { } +void CDiscAdjMultizoneDriver::DebugRun() { + + cout <<"\n------------------------------ Start Debug Run -----------------------------" << endl; + + cout <<"\n------------------------------ Check Objective Function Tape ---------------" << endl; + + cout << "Initial recording ..." << endl; + /*--- This recording will assign the initial (same) tag to each registered variable. + * During the recording, each dependent variable will be assigned the same tag. ---*/ + SetRecording(RECORDING::TAG_INIT_SOLUTION_VARIABLES, Kind_Tape::OBJECTIVE_FUNCTION_TAPE, ZONE_0); + + cout << "Second recording to check first recording ..." << endl; + /*--- This recording repeats the initial recording with a different tag. + * If a variable was used before it became dependent on the inputs, this variable will still carry the tag + * from the initial recording and a mismatch with the "check" recording tag will throw an error. + * In such a case, a possible reason could be that such a variable is set by a post-processing routine while + * for a mathematically correct recording this dependency must be included earlier. ---*/ + SetRecording(RECORDING::TAG_CHECK_SOLUTION_VARIABLES, Kind_Tape::OBJECTIVE_FUNCTION_TAPE, ZONE_0); + + cout <<"\n------------------------------ End Debug Run -----------------------------" << endl; +} + bool CDiscAdjMultizoneDriver::Iterate(unsigned short iZone, unsigned long iInnerIter, bool KrylovMode) { config_container[iZone]->SetInnerIter(iInnerIter); @@ -597,6 +629,15 @@ void CDiscAdjMultizoneDriver::SetRecording(RECORDING kind_recording, Kind_Tape t if(kind_recording != RECORDING::CLEAR_INDICES) { + if (kind_recording == RECORDING::TAG_INIT_SOLUTION_VARIABLES) { + cout << "Set Tag 1" << endl; + AD::SetTag(1); + } + else if (kind_recording == RECORDING::TAG_CHECK_SOLUTION_VARIABLES) { + cout << "Set Tag 2" << endl; + AD::SetTag(2); + } + AD::StartRecording(); AD::Push_TapePosition(); /// START @@ -738,10 +779,14 @@ void CDiscAdjMultizoneDriver::SetObjFunction(RECORDING kind_recording) { } } + + if (rank == MASTER_NODE) { AD::RegisterOutput(ObjFunc); AD::SetIndex(ObjFunc_Index, ObjFunc); - if (kind_recording == RECORDING::SOLUTION_VARIABLES) { + if (kind_recording == RECORDING::SOLUTION_VARIABLES || + kind_recording == RECORDING::TAG_INIT_SOLUTION_VARIABLES || + kind_recording == RECORDING::TAG_CHECK_SOLUTION_VARIABLES) { cout << " Objective function : " << ObjFunc; if (driver_config->GetWrt_AD_Statistics()){ // cout << " (" << ObjFunc_Index << ")\n"; diff --git a/SU2_CFD/src/drivers/CDriver.cpp b/SU2_CFD/src/drivers/CDriver.cpp index 73063e2ed7d..0c3c4a10c4e 100644 --- a/SU2_CFD/src/drivers/CDriver.cpp +++ b/SU2_CFD/src/drivers/CDriver.cpp @@ -2476,14 +2476,23 @@ void CDriver::InitializeInterface(CConfig **config, CSolver***** solver, CGeomet if (rank == MASTER_NODE) cout << "boundary displacements from the structural solver." << endl; } else if (fluid_donor && fluid_target) { - /*--- Mixing plane for turbo machinery applications. ---*/ - if (config[donor]->GetBoolMixingPlaneInterface()) { - interface_type = MIXING_PLANE; - auto nVar = solver[donor][INST_0][MESH_0][FLOW_SOL]->GetnVar(); - interface[donor][target] = new CMixingPlaneInterface(nVar, 0); - if (rank == MASTER_NODE) { - cout << "Set mixing-plane interface from donor zone " - << donor << " to target zone " << target << "." << endl; + /*--- Interface handling for turbomachinery applications. ---*/ + if (config[donor]->GetBoolTurbomachinery()) { + auto interfaceIndex = donor+target; // Here we assume that the interfaces at each side are the same kind + switch (config[donor]->GetKind_TurboInterface(interfaceIndex)) { + case TURBO_INTERFACE_KIND::MIXING_PLANE: { + interface_type = MIXING_PLANE; + auto nVar = solver[donor][INST_0][MESH_0][FLOW_SOL]->GetnVar(); + interface[donor][target] = new CMixingPlaneInterface(nVar, 0); + if (rank == MASTER_NODE) cout << "using a mixing-plane interface from donor zone " << donor << " to target zone " << target << "." << endl; + break; + } + case TURBO_INTERFACE_KIND::FROZEN_ROTOR: { + auto nVar = solver[donor][INST_0][MESH_0][FLOW_SOL]->GetnPrimVar(); + interface_type = SLIDING_INTERFACE; + interface[donor][target] = new CSlidingInterface(nVar, 0); + if (rank == MASTER_NODE) cout << "using a fluid interface interface from donor zone " << donor << " to target zone " << target << "." << endl; + } } } else{ @@ -2494,19 +2503,21 @@ void CDriver::InitializeInterface(CConfig **config, CSolver***** solver, CGeomet } } else if (heat_donor || heat_target) { - if (heat_donor && heat_target) - SU2_MPI::Error("Conjugate heat transfer between solids is not implemented.", CURRENT_FUNCTION); + if (heat_donor && heat_target){ + interface_type = CONJUGATE_HEAT_SS; - const auto fluidZone = heat_target? donor : target; - - if (config[fluidZone]->GetEnergy_Equation() || (config[fluidZone]->GetKind_Regime() == ENUM_REGIME::COMPRESSIBLE) - || (config[fluidZone]->GetKind_FluidModel() == ENUM_FLUIDMODEL::FLUID_FLAMELET)) - interface_type = heat_target? CONJUGATE_HEAT_FS : CONJUGATE_HEAT_SF; - else if (config[fluidZone]->GetWeakly_Coupled_Heat()) - interface_type = heat_target? CONJUGATE_HEAT_WEAKLY_FS : CONJUGATE_HEAT_WEAKLY_SF; - else - interface_type = NO_TRANSFER; + } else { + const auto fluidZone = heat_target? donor : target; + if (config[fluidZone]->GetEnergy_Equation() || (config[fluidZone]->GetKind_Regime() == ENUM_REGIME::COMPRESSIBLE) + || (config[fluidZone]->GetKind_FluidModel() == ENUM_FLUIDMODEL::FLUID_FLAMELET)) + interface_type = heat_target? CONJUGATE_HEAT_FS : CONJUGATE_HEAT_SF; + else if (config[fluidZone]->GetWeakly_Coupled_Heat()) + interface_type = heat_target? CONJUGATE_HEAT_WEAKLY_FS : CONJUGATE_HEAT_WEAKLY_SF; + else + interface_type = NO_TRANSFER; + } + if (interface_type != NO_TRANSFER) { auto nVar = 4; interface[donor][target] = new CConjugateHeatInterface(nVar, 0); diff --git a/SU2_CFD/src/drivers/CMultizoneDriver.cpp b/SU2_CFD/src/drivers/CMultizoneDriver.cpp index 607e1d989ee..368a598e514 100644 --- a/SU2_CFD/src/drivers/CMultizoneDriver.cpp +++ b/SU2_CFD/src/drivers/CMultizoneDriver.cpp @@ -563,6 +563,9 @@ bool CMultizoneDriver::TransferData(unsigned short donorZone, unsigned short tar BroadcastData(SPECIES_SOL, SPECIES_SOL); } break; + case CONJUGATE_HEAT_SS: + BroadcastData(HEAT_SOL, HEAT_SOL); + break; case CONJUGATE_HEAT_FS: BroadcastData(FLOW_SOL, HEAT_SOL); break; @@ -587,6 +590,7 @@ bool CMultizoneDriver::TransferData(unsigned short donorZone, unsigned short tar const auto nMarkerInt = config_container[donorZone]->GetnMarker_MixingPlaneInterface() / 2; /*--- Transfer the average value from the donorZone to the targetZone ---*/ + /*--- Loops over the mixing planes defined in the config file to find the correct mixing plane for the donor-target combination ---*/ for (auto iMarkerInt = 1; iMarkerInt <= nMarkerInt; iMarkerInt++) { interface_container[donorZone][targetZone]->AllgatherAverage(solver_container[donorZone][INST_0][MESH_0][FLOW_SOL],solver_container[targetZone][INST_0][MESH_0][FLOW_SOL], geometry_container[donorZone][INST_0][MESH_0],geometry_container[targetZone][INST_0][MESH_0], diff --git a/SU2_CFD/src/drivers/CSinglezoneDriver.cpp b/SU2_CFD/src/drivers/CSinglezoneDriver.cpp index 6b1bd3eaaaa..627ea9c41f7 100644 --- a/SU2_CFD/src/drivers/CSinglezoneDriver.cpp +++ b/SU2_CFD/src/drivers/CSinglezoneDriver.cpp @@ -66,8 +66,8 @@ void CSinglezoneDriver::StartSolver() { TimeIter = config_container[ZONE_0]->GetRestart_Iter(); /*--- Run the problem until the number of time iterations required is reached. ---*/ - while ( TimeIter < config_container[ZONE_0]->GetnTime_Iter() ) { - + /*--- or until a SIGTERM signal stops the loop. We catch SIGTERM and exit gracefully ---*/ + while ( TimeIter < config_container[ZONE_0]->GetnTime_Iter()) { /*--- Perform some preprocessing before starting the time-step simulation. ---*/ Preprocess(TimeIter); diff --git a/SU2_CFD/src/fluid/CFluidScalar.cpp b/SU2_CFD/src/fluid/CFluidScalar.cpp index 991e29bd562..5e3a2e2f998 100644 --- a/SU2_CFD/src/fluid/CFluidScalar.cpp +++ b/SU2_CFD/src/fluid/CFluidScalar.cpp @@ -41,12 +41,9 @@ #include "../../include/fluid/CPolynomialViscosity.hpp" #include "../../include/fluid/CSutherland.hpp" -CFluidScalar::CFluidScalar(su2double val_Cp, su2double val_gas_constant, su2double value_pressure_operating, - const CConfig* config) +CFluidScalar::CFluidScalar(su2double value_pressure_operating, const CConfig* config) : CFluidModel(), n_species_mixture(config->GetnSpecies() + 1), - Gas_Constant(val_gas_constant), - Gamma(config->GetGamma()), Pressure_Thermodynamic(value_pressure_operating), GasConstant_Ref(config->GetGas_Constant_Ref()), Prandtl_Number(config->GetPrandtl_Turb()), diff --git a/SU2_CFD/src/interfaces/CInterface.cpp b/SU2_CFD/src/interfaces/CInterface.cpp index d79c2d316d3..b9617c4c16e 100644 --- a/SU2_CFD/src/interfaces/CInterface.cpp +++ b/SU2_CFD/src/interfaces/CInterface.cpp @@ -108,6 +108,11 @@ void CInterface::BroadcastData(const CInterpolator& interpolator, su2activematrix sendDonorVar(nLocalVertexDonor, nVar); if (markDonor >= 0) { + + /*--- Apply contact resistance if specified. ---*/ + + SetContactResistance(donor_config->GetContactResistance(iMarkerInt)); + for (auto iVertex = 0ul, iSend = 0ul; iVertex < donor_geometry->GetnVertex(markDonor); iVertex++) { const auto iPoint = donor_geometry->vertex[markDonor][iVertex]->GetNode(); @@ -248,7 +253,7 @@ void CInterface::PreprocessAverage(CGeometry *donor_geometry, CGeometry *target_ Donor_Flag= -1; for (int iSize=0; iSize 0.0){ + if(BuffMarkerDonor[iSize] >= 0.0){ Marker_Donor = BuffMarkerDonor[iSize]; Donor_Flag = BuffDonorFlag[iSize]; break; @@ -273,8 +278,8 @@ void CInterface::PreprocessAverage(CGeometry *donor_geometry, CGeometry *target_ break; } /*--- If the tag hasn't matched any tag within the Flow markers ---*/ - Marker_Target = -1; - + Marker_Target = -1; + Target_Flag = -1; } if (Marker_Target != -1 && Marker_Donor != -1){ @@ -295,7 +300,7 @@ void CInterface::PreprocessAverage(CGeometry *donor_geometry, CGeometry *target_ } if(test2 < dist2){ dist2 = test2; - tSpan =jSpan; + tSpan = jSpan; } } @@ -633,6 +638,4 @@ void CInterface::AllgatherAverage(CSolver *donor_solution, CSolver *target_solut delete [] avgNuTarget; delete [] avgKineTarget; delete [] avgOmegaTarget; - - } diff --git a/SU2_CFD/src/interfaces/cfd/CMixingPlaneInterface.cpp b/SU2_CFD/src/interfaces/cfd/CMixingPlaneInterface.cpp index f35dc776f9e..d1bbcf2ff0b 100644 --- a/SU2_CFD/src/interfaces/cfd/CMixingPlaneInterface.cpp +++ b/SU2_CFD/src/interfaces/cfd/CMixingPlaneInterface.cpp @@ -118,7 +118,7 @@ void CMixingPlaneInterface::SetAverageValues(CSolver *donor_solution, CSolver *t unsigned short iSpan; for(iSpan = 0; iSpanSetDensityIn(donor_solution->GetDensityIn(donorZone, iSpan), donorZone, iSpan); target_solution->SetPressureIn(donor_solution->GetPressureIn(donorZone, iSpan), donorZone, iSpan); target_solution->SetTurboVelocityIn(donor_solution->GetTurboVelocityIn(donorZone, iSpan), donorZone, iSpan); diff --git a/SU2_CFD/src/interfaces/cht/CConjugateHeatInterface.cpp b/SU2_CFD/src/interfaces/cht/CConjugateHeatInterface.cpp index 54ad481b488..908ea27ae7f 100644 --- a/SU2_CFD/src/interfaces/cht/CConjugateHeatInterface.cpp +++ b/SU2_CFD/src/interfaces/cht/CConjugateHeatInterface.cpp @@ -132,8 +132,10 @@ void CConjugateHeatInterface::GetDonor_Variable(CSolver *donor_solution, CGeomet if ((donor_config->GetKind_CHT_Coupling() == CHT_COUPLING::DIRECT_TEMPERATURE_ROBIN_HEATFLUX) || (donor_config->GetKind_CHT_Coupling() == CHT_COUPLING::AVERAGED_TEMPERATURE_ROBIN_HEATFLUX)) { + /*--- Apply contact resistance to solid-to-solid heat transfer boundary ---*/ const su2double rho_cp_solid = donor_config->GetSpecific_Heat_Cp()*donor_config->GetMaterialDensity(0); - conductivity_over_dist = thermal_diffusivity*rho_cp_solid/dist; + thermal_conductivity = thermal_diffusivity * rho_cp_solid; + conductivity_over_dist = thermal_conductivity/(dist + thermal_conductivity * ContactResistance); } } diff --git a/SU2_CFD/src/iteration/CDiscAdjFluidIteration.cpp b/SU2_CFD/src/iteration/CDiscAdjFluidIteration.cpp index 24fd1602e6e..e7e6feb8e78 100644 --- a/SU2_CFD/src/iteration/CDiscAdjFluidIteration.cpp +++ b/SU2_CFD/src/iteration/CDiscAdjFluidIteration.cpp @@ -404,7 +404,10 @@ void CDiscAdjFluidIteration::RegisterInput(CSolver***** solver, CGeometry**** ge SU2_OMP_PARALLEL_(if(solvers0[ADJFLOW_SOL]->GetHasHybridParallel())) { - if (kind_recording == RECORDING::SOLUTION_VARIABLES || kind_recording == RECORDING::SOLUTION_AND_MESH) { + if (kind_recording == RECORDING::SOLUTION_VARIABLES || + kind_recording == RECORDING::TAG_INIT_SOLUTION_VARIABLES || + kind_recording == RECORDING::TAG_CHECK_SOLUTION_VARIABLES || + kind_recording == RECORDING::SOLUTION_AND_MESH) { /*--- Register flow and turbulent variables as input ---*/ if (config[iZone]->GetFluidProblem()) { diff --git a/SU2_CFD/src/output/CFlowCompOutput.cpp b/SU2_CFD/src/output/CFlowCompOutput.cpp index 5bb5577fe76..e693bbbae47 100644 --- a/SU2_CFD/src/output/CFlowCompOutput.cpp +++ b/SU2_CFD/src/output/CFlowCompOutput.cpp @@ -592,6 +592,10 @@ void CFlowCompOutput::LoadTurboHistoryData(std::shared_ptrGetOutletState().GetTotalPressure()); SetHistoryOutputValue("PressureIn_" + tag.str(), BladePerf->GetInletState().GetPressure()); SetHistoryOutputValue("PressureOut_" + tag.str(), BladePerf->GetOutletState().GetPressure()); + SetHistoryOutputValue("TotalTemperatureIn_" + tag.str(), BladePerf->GetInletState().GetTotalTemperature()); + SetHistoryOutputValue("TotalTemperatureOut_" + tag.str(), BladePerf->GetOutletState().GetTotalTemperature()); + SetHistoryOutputValue("TemperatureIn_" + tag.str(), BladePerf->GetInletState().GetTemperature()); + SetHistoryOutputValue("TemperatureOut_" + tag.str(), BladePerf->GetOutletState().GetTemperature()); SetHistoryOutputValue("DensityIn_" + tag.str(), BladePerf->GetInletState().GetDensity()); SetHistoryOutputValue("DensityOut_" + tag.str(), BladePerf->GetOutletState().GetDensity()); SetHistoryOutputValue("NormalVelocityIn_" + tag.str(), BladePerf->GetInletState().GetVelocity()[0]); @@ -604,6 +608,8 @@ void CFlowCompOutput::LoadTurboHistoryData(std::shared_ptrGetOutletState().GetMachValue()); SetHistoryOutputValue("AbsFlowAngleIn_" + tag.str(), BladePerf->GetInletState().GetAbsFlowAngle()*180/PI_NUMBER); SetHistoryOutputValue("AbsFlowAngleOut_" + tag.str(), BladePerf->GetOutletState().GetAbsFlowAngle()*180/PI_NUMBER); + SetHistoryOutputValue("KineticEnergyLoss_" + tag.str(), BladePerf->GetKineticEnergyLoss()); + SetHistoryOutputValue("TotPressureLoss_" + tag.str(), BladePerf->GetTotalPressureLoss()); } SetHistoryOutputValue("EntropyGeneration", TurboStagePerf->GetNormEntropyGen()*100); SetHistoryOutputValue("EulerianWork", TurboStagePerf->GetEulerianWork()); @@ -611,6 +617,8 @@ void CFlowCompOutput::LoadTurboHistoryData(std::shared_ptrGetTotalTotalEfficiency()*100); SetHistoryOutputValue("PressureRatioTS", TurboStagePerf->GetTotalStaticPressureRatio()); SetHistoryOutputValue("PressureRatioTT", TurboStagePerf->GetTotalTotalPressureRatio()); + SetHistoryOutputValue("KineticEnergyLoss_Stage", TurboStagePerf->GetKineticEnergyLoss()); + SetHistoryOutputValue("TotPressureLoss_Stage", TurboStagePerf->GetTotalPressureLoss()); } void CFlowCompOutput::WriteTurboSpanwisePerformance(std::shared_ptr TurboPerf, CGeometry *geometry, CConfig **config, unsigned short val_iZone) { diff --git a/SU2_CFD/src/output/CFlowOutput.cpp b/SU2_CFD/src/output/CFlowOutput.cpp index 5ffee0d23f5..895419f264d 100644 --- a/SU2_CFD/src/output/CFlowOutput.cpp +++ b/SU2_CFD/src/output/CFlowOutput.cpp @@ -4066,6 +4066,10 @@ void CFlowOutput::AddTurboOutput(unsigned short nZone){ AddHistoryOutput("TotalPressureOut_" + tag, "TotPressureOut_" + tag, ScreenOutputFormat::SCIENTIFIC, "TURBO_PERF", "Flow angle in " + tag, HistoryFieldType::DEFAULT); AddHistoryOutput("PressureIn_" + tag, "PressureIn_" + tag, ScreenOutputFormat::SCIENTIFIC, "TURBO_PERF", "Pressure ratio " + tag, HistoryFieldType::DEFAULT); AddHistoryOutput("PressureOut_" + tag, "PressureOut_" + tag, ScreenOutputFormat::SCIENTIFIC, "TURBO_PERF", "Flow angle in " + tag, HistoryFieldType::DEFAULT); + AddHistoryOutput("TotalTemperatureIn_" + tag, "TotTemperatureIn_" + tag, ScreenOutputFormat::SCIENTIFIC, "TURBO_PERF", "Temperature ratio " + tag, HistoryFieldType::DEFAULT); + AddHistoryOutput("TotalTemperatureOut_" + tag, "TotTemperatureOut_" + tag, ScreenOutputFormat::SCIENTIFIC, "TURBO_PERF", "Flow angle in " + tag, HistoryFieldType::DEFAULT); + AddHistoryOutput("TemperatureIn_" + tag, "TemperatureIn_" + tag, ScreenOutputFormat::SCIENTIFIC, "TURBO_PERF", "Temperature ratio " + tag, HistoryFieldType::DEFAULT); + AddHistoryOutput("TemperatureOut_" + tag, "TemperatureOut_" + tag, ScreenOutputFormat::SCIENTIFIC, "TURBO_PERF", "Flow angle in " + tag, HistoryFieldType::DEFAULT); AddHistoryOutput("DensityIn_" + tag, "DensityIn_" + tag, ScreenOutputFormat::SCIENTIFIC, "TURBO_PERF", "Flow angle out " + tag, HistoryFieldType::DEFAULT); AddHistoryOutput("DensityOut_" + tag, "DensityOut_" + tag, ScreenOutputFormat::SCIENTIFIC, "TURBO_PERF", "Absolute flow angle in " + tag, HistoryFieldType::DEFAULT); AddHistoryOutput("NormalVelocityIn_" + tag, "NormalVelocityIn_" + tag, ScreenOutputFormat::SCIENTIFIC, "TURBO_PERF", "Absolute flow angle out " + tag, HistoryFieldType::DEFAULT); @@ -4078,12 +4082,16 @@ void CFlowOutput::AddTurboOutput(unsigned short nZone){ AddHistoryOutput("MachOut_" + tag, "MachOut_" + tag, ScreenOutputFormat::SCIENTIFIC, "TURBO_PERF", "Total-to-Static efficiency " + tag, HistoryFieldType::DEFAULT); AddHistoryOutput("AbsFlowAngleIn_" + tag, "AbsFlowAngleIn_" + tag, ScreenOutputFormat::SCIENTIFIC, "TURBO_PERF", "Absolute flow angle in " + tag, HistoryFieldType::DEFAULT); AddHistoryOutput("AbsFlowAngleOut_" + tag, "AbsFlowAngleOut_" + tag, ScreenOutputFormat::SCIENTIFIC, "TURBO_PERF", "Absolute flow angle out " + tag, HistoryFieldType::DEFAULT); + AddHistoryOutput("KineticEnergyLoss_" + tag, "KELC_" + tag, ScreenOutputFormat::SCIENTIFIC, "TURBO_PERF", "Blade Kinetic Energy Loss Coefficient", HistoryFieldType::DEFAULT); + AddHistoryOutput("TotPressureLoss_" + tag, "TPLC_" + tag, ScreenOutputFormat::SCIENTIFIC, "TURBO_PERF", "Blade Pressure Loss Coefficient", HistoryFieldType::DEFAULT); } //Adds turbomachinery machine performance variables AddHistoryOutput("EntropyGeneration", "EntropyGen", ScreenOutputFormat::SCIENTIFIC, "TURBO_PERF", "Machine entropy generation", HistoryFieldType::DEFAULT); - AddHistoryOutput("EulerianWork", "EulerianWork", ScreenOutputFormat::SCIENTIFIC, "TURBO_PERF", "Machine entropy generation", HistoryFieldType::DEFAULT); - AddHistoryOutput("TotalStaticEfficiency", "TotStaticEff", ScreenOutputFormat::SCIENTIFIC, "TURBO_PERF", "Machine entropy generation", HistoryFieldType::DEFAULT); - AddHistoryOutput("TotalTotalEfficiency", "TotTotEff", ScreenOutputFormat::SCIENTIFIC, "TURBO_PERF", "Machine entropy generation", HistoryFieldType::DEFAULT); - AddHistoryOutput("PressureRatioTS", "PRTS", ScreenOutputFormat::SCIENTIFIC, "TURBO_PERF", "Machine entropy generation", HistoryFieldType::DEFAULT); - AddHistoryOutput("PressureRatioTT", "PRTT", ScreenOutputFormat::SCIENTIFIC, "TURBO_PERF", "Machine entropy generation", HistoryFieldType::DEFAULT); + AddHistoryOutput("EulerianWork", "EulerianWork", ScreenOutputFormat::SCIENTIFIC, "TURBO_PERF", "Machine Eulerian work", HistoryFieldType::DEFAULT); + AddHistoryOutput("TotalStaticEfficiency", "TotStaticEff", ScreenOutputFormat::SCIENTIFIC, "TURBO_PERF", "Machine total-to-static efficiency", HistoryFieldType::DEFAULT); + AddHistoryOutput("TotalTotalEfficiency", "TotTotEff", ScreenOutputFormat::SCIENTIFIC, "TURBO_PERF", "Machine total-to-total efficiency", HistoryFieldType::DEFAULT); + AddHistoryOutput("PressureRatioTS", "PRTS", ScreenOutputFormat::SCIENTIFIC, "TURBO_PERF", "Machine total-to-static pressure ratio", HistoryFieldType::DEFAULT); + AddHistoryOutput("PressureRatioTT", "PRTT", ScreenOutputFormat::SCIENTIFIC, "TURBO_PERF", "Machine total-to-toal pressure ratio", HistoryFieldType::DEFAULT); + AddHistoryOutput("KineticEnergyLoss_Stage", "KELC_all", ScreenOutputFormat::SCIENTIFIC, "TURBO_PERF", "Machine Kinetic Energy Loss Coefficient", HistoryFieldType::DEFAULT); + AddHistoryOutput("TotPressureLoss_Stage", "TPLC_all", ScreenOutputFormat::SCIENTIFIC, "TURBO_PERF", "Machine Pressure Loss Coefficient", HistoryFieldType::DEFAULT); } diff --git a/SU2_CFD/src/output/COutput.cpp b/SU2_CFD/src/output/COutput.cpp index d8aa5c16af9..499d0d69fd5 100644 --- a/SU2_CFD/src/output/COutput.cpp +++ b/SU2_CFD/src/output/COutput.cpp @@ -25,6 +25,9 @@ * License along with SU2. If not, see . */ +#include +#include + #include "../../../Common/include/geometry/CGeometry.hpp" #include "../../include/solvers/CSolver.hpp" @@ -47,6 +50,15 @@ #include "../../include/output/filewriter/CSU2BinaryFileWriter.hpp" #include "../../include/output/filewriter/CSU2MeshFileWriter.hpp" +namespace { +volatile sig_atomic_t STOP; + +void signalHandler(int signum) { + std::cout << "Interrupt signal (" << signum << ") received, saving files and exiting.\n"; + STOP = 1; +} +} + COutput::COutput(const CConfig *config, unsigned short ndim, bool fem_output): rank(SU2_MPI::GetRank()), size(SU2_MPI::GetSize()), @@ -169,7 +181,11 @@ COutput::COutput(const CConfig *config, unsigned short ndim, bool fem_output): volumeDataSorter = nullptr; surfaceDataSorter = nullptr; - headerNeeded = false; + headerNeeded = false; + + /*--- Setup a signal handler for SIGTERM. ---*/ + + signal(SIGTERM, signalHandler); } COutput::~COutput() { @@ -233,7 +249,7 @@ void COutput::SetHistoryOutput(CGeometry ****geometry, CSolver *****solver, CCon if (config[ZONE_0]->GetMultizone_Problem()) Iter = OuterIter; - + /*--- Turbomachinery Performance Screen summary output---*/ if (Iter%100 == 0 && rank == MASTER_NODE) { SetTurboPerformance_Output(TurboPerf, config[val_iZone], TimeIter, OuterIter, InnerIter); @@ -946,9 +962,13 @@ bool COutput::ConvergenceMonitoring(CConfig *config, unsigned long Iteration) { if (convFields.empty() || Iteration < config->GetStartConv_Iter()) convergence = false; + /*--- If a SIGTERM signal is sent to one of the processes, we set convergence to true. ---*/ + if (STOP) convergence = true; + /*--- Apply the same convergence criteria to all processors. ---*/ unsigned short local = convergence, global = 0; + SU2_MPI::Allreduce(&local, &global, 1, MPI_UNSIGNED_SHORT, MPI_MAX, SU2_MPI::GetComm()); convergence = global > 0; diff --git a/SU2_CFD/src/output/CTurboOutput.cpp b/SU2_CFD/src/output/CTurboOutput.cpp index 6551d0062c1..86ed50cc804 100644 --- a/SU2_CFD/src/output/CTurboOutput.cpp +++ b/SU2_CFD/src/output/CTurboOutput.cpp @@ -61,6 +61,7 @@ void CTurbomachineryState::ComputeState(CFluidModel& fluidModel, const CTurbomac std::vector velocity = primitiveState.GetVelocity(); Velocity.assign(velocity.begin(), velocity.end()); TangVelocity = primitiveState.GetTangVelocity(); + TangVelocity = primitiveState.GetTangVelocity(); /*--- Compute static TD quantities ---*/ fluidModel.SetTDState_Prho(Pressure, Density); @@ -82,8 +83,10 @@ void CTurbomachineryState::ComputeState(CFluidModel& fluidModel, const CTurbomac /*--- Compute relative kinematic quantities ---*/ su2double tangVel2 = TangVelocity * TangVelocity; + su2double tangVel2 = TangVelocity * TangVelocity; RelVelocity.assign(Velocity.begin(), Velocity.end()); RelVelocity[1] -= TangVelocity; + RelVelocity[1] -= TangVelocity; su2double relVel2 = GetRelVelocityValue(); FlowAngle = atan(RelVelocity[1] / RelVelocity[0]); RelMach.assign(RelVelocity.begin(), RelVelocity.end()); diff --git a/SU2_CFD/src/python_wrapper_structure.cpp b/SU2_CFD/src/python_wrapper_structure.cpp index a0018371443..86933d11341 100644 --- a/SU2_CFD/src/python_wrapper_structure.cpp +++ b/SU2_CFD/src/python_wrapper_structure.cpp @@ -111,7 +111,7 @@ void CSinglezoneDriver::SetInitialMesh() { DynamicMeshUpdate(0); SU2_OMP_PARALLEL { - for (iMesh = 0u; iMesh <= main_config->GetnMGLevels(); iMesh++) { + for (auto iMesh = 0u; iMesh <= main_config->GetnMGLevels(); iMesh++) { SU2_OMP_FOR_STAT(roundUpDiv(geometry_container[selected_zone][INST_0][iMesh]->GetnPoint(), omp_get_max_threads())) for (auto iPoint = 0ul; iPoint < geometry_container[selected_zone][INST_0][iMesh]->GetnPoint(); iPoint++) { /*--- Overwrite fictitious velocities. ---*/ diff --git a/SU2_CFD/src/solvers/CEulerSolver.cpp b/SU2_CFD/src/solvers/CEulerSolver.cpp index 6267635a30e..286a84e21c7 100644 --- a/SU2_CFD/src/solvers/CEulerSolver.cpp +++ b/SU2_CFD/src/solvers/CEulerSolver.cpp @@ -232,6 +232,9 @@ CEulerSolver::CEulerSolver(CGeometry *geometry, CConfig *config, Exhaust_Pressure.resize(nMarker); Exhaust_Area.resize(nMarker); + /*--- Turbomachinery simulation ---*/ + AverageMassFlowRate.resize(nMarker); + /*--- Read farfield conditions from config ---*/ Temperature_Inf = config->GetTemperature_FreeStreamND(); @@ -279,6 +282,8 @@ CEulerSolver::CEulerSolver(CGeometry *geometry, CConfig *config, Exhaust_Temperature[iMarker] = Temperature_Inf; Exhaust_Pressure[iMarker] = Pressure_Inf; Exhaust_Area[iMarker] = 0.0; + + AverageMassFlowRate[iMarker] = 0.0; } /*--- Initialize the solution to the far-field state everywhere. ---*/ @@ -2541,7 +2546,7 @@ void CEulerSolver::GetPower_Properties(CGeometry *geometry, CConfig *config, uns su2double Alpha = config->GetAoA()*PI_NUMBER/180.0; su2double Beta = config->GetAoS()*PI_NUMBER/180.0; bool write_heads = ((((config->GetInnerIter() % (config->GetScreen_Wrt_Freq(2)*40)) == 0) && (config->GetInnerIter()!= 0)) || (config->GetInnerIter() == 1)); - bool Evaluate_BC = ((((config->GetInnerIter() % (config->GetScreen_Wrt_Freq(2)*40)) == 0)) || (config->GetInnerIter() == 1) || (config->GetDiscrete_Adjoint())); + bool Evaluate_BC = ((((config->GetInnerIter() % (config->GetBc_Eval_Freq())) == 0)) || (config->GetInnerIter() == 1) || (config->GetDiscrete_Adjoint())); if ((config->GetnMarker_EngineInflow() != 0) || (config->GetnMarker_EngineExhaust() != 0)) Engine = true; if ((config->GetnMarker_ActDiskInlet() != 0) || (config->GetnMarker_ActDiskOutlet() != 0)) Engine = false; @@ -2930,6 +2935,7 @@ void CEulerSolver::GetPower_Properties(CGeometry *geometry, CConfig *config, uns config->SetInflow_RamDrag(iMarker_Inlet, Inlet_RamDrag_Total[iMarker_Inlet]); config->SetInflow_Force(iMarker_Inlet, Inlet_Force_Total[iMarker_Inlet]); config->SetInflow_Power(iMarker_Inlet, Inlet_Power_Total[iMarker_Inlet]); + config->SetInflow_Mach(iMarker_Inlet, Inlet_Mach_Total[iMarker_Inlet]); } else { config->SetActDiskInlet_MassFlow(iMarker_Inlet, Inlet_MassFlow_Total[iMarker_Inlet]); @@ -6115,60 +6121,81 @@ void CEulerSolver::PreprocessBC_Giles(CGeometry *geometry, CConfig *config, CNum void CEulerSolver::BC_Giles(CGeometry *geometry, CSolver **solver_container, CNumerics *conv_numerics, CNumerics *visc_numerics, CConfig *config, unsigned short val_marker) { - /*--- Initialisation of data structures ---*/ - - su2double *Normal = new su2double[nDim], - *turboNormal = new su2double[nDim], - *UnitNormal = new su2double[nDim], - *turboVelocity = new su2double[nDim], - *Velocity_i = new su2double[nDim], - *Velocity_b = new su2double[nDim], - *AverageTurboMach = new su2double[nDim], - *S_boundary = new su2double[8], - *delta_c = new su2double[nVar], - *cj = new su2double[nVar], - **R_Matrix = new su2double*[nVar], - *dcjs = new su2double[nVar], - *c_avg = new su2double[nVar], - **R_c = new su2double*[nVar-1], - **R_c_inv = new su2double*[nVar-1], - *R = new su2double[nVar-1], - *deltaprim = new su2double[nVar]; - - for (auto iVar = 0; iVar < nVar; iVar++) + unsigned short iDim, iVar, jVar, iSpan; + unsigned long iPoint, Point_Normal, oldVertex, k, kend, kend_max, iVertex; + su2double *UnitNormal, *turboVelocity, *turboNormal; + + su2double *Velocity_b, Velocity2_b, Enthalpy_b, Energy_b, Density_b, Pressure_b, Temperature_b; + su2double *Velocity_i, Velocity2_i, Energy_i, StaticEnergy_i, Density_i, Pressure_i; + su2double Pressure_e; + su2double *V_boundary, *V_domain, *S_boundary, *S_domain; + unsigned short iZone = config->GetiZone(); + bool implicit = (config->GetKind_TimeIntScheme() == EULER_IMPLICIT); + string Marker_Tag = config->GetMarker_All_TagBound(val_marker); + bool viscous = config->GetViscous(); + unsigned short nSpanWiseSections = geometry->GetnSpanWiseSections(config->GetMarker_All_TurbomachineryFlag(val_marker)); + su2double relfacAvgCfg = config->GetGiles_RelaxFactorAverage(Marker_Tag); + su2double relfacFouCfg = config->GetGiles_RelaxFactorFourier(Marker_Tag); + su2double *Normal; + su2double TwoPiThetaFreq_Pitch, pitch,theta; + const su2double *SpanWiseValues = nullptr, *FlowDir; + su2double spanPercent, extrarelfacAvg = 0.0, deltaSpan = 0.0, relfacAvg, relfacFou, coeffrelfacAvg = 0.0; + unsigned short Turbo_Flag; + + Normal = new su2double[nDim]; + turboNormal = new su2double[nDim]; + UnitNormal = new su2double[nDim]; + turboVelocity = new su2double[nDim]; + Velocity_i = new su2double[nDim]; + Velocity_b = new su2double[nDim]; + + + su2double AverageSoundSpeed, *AverageTurboMach, AverageEntropy, AverageEnthalpy; + AverageTurboMach = new su2double[nDim]; + S_boundary = new su2double[8]; + + su2double AvgMach , *cj, GilesBeta, *delta_c, **R_Matrix, *deltaprim, **R_c_inv,**R_c, alphaIn_BC, gammaIn_BC = 0, + P_Total, T_Total, Enthalpy_BC, Entropy_BC, *R, *c_avg,*dcjs, Beta_inf2, c2js_Re, c3js_Re, cOutjs_Re, avgVel2 =0.0; + + long freq; + + delta_c = new su2double[nVar]; + deltaprim = new su2double[nVar]; + cj = new su2double[nVar]; + R_Matrix = new su2double*[nVar]; + R_c = new su2double*[nVar-1]; + R_c_inv = new su2double*[nVar-1]; + R = new su2double[nVar-1]; + c_avg = new su2double[nVar]; + dcjs = new su2double[nVar]; + + for (iVar = 0; iVar < nVar; iVar++) { R_Matrix[iVar] = new su2double[nVar]; c_avg[iVar] = 0.0; dcjs[iVar] = 0.0; } - for (auto iVar = 0; iVar < nVar-1; iVar++) + for (iVar = 0; iVar < nVar-1; iVar++) { R_c[iVar] = new su2double[nVar-1]; R_c_inv[iVar] = new su2double[nVar-1]; } - auto const I = complex(0.0,1.0); - - /*--- Compute coeff for under relaxation of Avg and Fourier Coefficient for hub and shroud---*/ - auto const Marker_Tag = config->GetMarker_All_TagBound(val_marker); - su2double relfacAvgCfg = config->GetGiles_RelaxFactorAverage(Marker_Tag); - - unsigned short nSpanWiseSections = geometry->GetnSpanWiseSections(config->GetMarker_All_TurbomachineryFlag(val_marker)); - - const su2double *SpanWiseValues = nullptr; - su2double extrarelfacAvg = 0.0, deltaSpan = 0.0, relfacAvg, relfacFou = config->GetGiles_RelaxFactorFourier(Marker_Tag), coeffrelfacAvg = 0.0; + complex I, c2ks, c2js, c3ks, c3js, cOutks, cOutjs, Beta_inf; + I = complex(0.0,1.0); + /*--- Compute coeff for under relaxation of Avg and Fourier Coefficient for hub and shroud---*/ if (nDim == 3){ extrarelfacAvg = config->GetExtraRelFacGiles(0); - auto const spanPercent = config->GetExtraRelFacGiles(1); - auto const Turbo_Flag = config->GetMarker_All_TurbomachineryFlag(val_marker); + spanPercent = config->GetExtraRelFacGiles(1); + Turbo_Flag = config->GetMarker_All_TurbomachineryFlag(val_marker); SpanWiseValues = geometry->GetSpanWiseValue(Turbo_Flag); deltaSpan = SpanWiseValues[nSpanWiseSections-1]*spanPercent; coeffrelfacAvg = (relfacAvgCfg - extrarelfacAvg)/deltaSpan; } - for (unsigned short iSpan= 0; iSpan < nSpanWiseSections ; iSpan++){ + for (iSpan= 0; iSpan < nSpanWiseSections ; iSpan++){ /*--- Compute under relaxation for the Hub and Shroud Avg and Fourier Coefficient---*/ if(nDim == 3){ if(SpanWiseValues[iSpan] <= SpanWiseValues[0] + deltaSpan){ @@ -6181,16 +6208,18 @@ void CEulerSolver::BC_Giles(CGeometry *geometry, CSolver **solver_container, CNu } else{ relfacAvg = relfacAvgCfg; + relfacFou = relfacFouCfg; } } else{ { relfacAvg = relfacAvgCfg; + relfacFou = relfacFouCfg; } } GetFluidModel()->SetTDState_Prho(AveragePressure[val_marker][iSpan], AverageDensity[val_marker][iSpan]); - auto AverageSoundSpeed = GetFluidModel()->GetSoundSpeed(); + AverageSoundSpeed = GetFluidModel()->GetSoundSpeed(); AverageTurboMach[0] = AverageTurboVelocity[val_marker][iSpan][0]/AverageSoundSpeed; AverageTurboMach[1] = AverageTurboVelocity[val_marker][iSpan][1]/AverageSoundSpeed; @@ -6198,32 +6227,167 @@ void CEulerSolver::BC_Giles(CGeometry *geometry, CSolver **solver_container, CNu AverageTurboMach[1] -= geometry->GetAverageTangGridVel(val_marker,iSpan)/AverageSoundSpeed; } - auto const kend = geometry->GetnFreqSpan(val_marker, iSpan); - auto const kend_max = geometry->GetnFreqSpanMax(config->GetMarker_All_TurbomachineryFlag(val_marker)); - conv_numerics->GetRMatrix(AverageSoundSpeed, AverageDensity[val_marker][iSpan], R_Matrix); + AvgMach = AverageTurboMach[0]*AverageTurboMach[0] + AverageTurboMach[1]*AverageTurboMach[1]; - su2double Pressure_e; //Useful to pass by reference + kend = geometry->GetnFreqSpan(val_marker, iSpan); + kend_max = geometry->GetnFreqSpanMax(config->GetMarker_All_TurbomachineryFlag(val_marker)); + conv_numerics->GetRMatrix(AverageSoundSpeed, AverageDensity[val_marker][iSpan], R_Matrix); switch(config->GetKind_Data_Giles(Marker_Tag)){ case TOTAL_CONDITIONS_PT: - BC_Giles_Total_Inlet(conv_numerics, config, val_marker, c_avg, R, R_c_inv, R_c, iSpan, iSpan); + /*--- Retrieve the specified total conditions for this inlet. ---*/ + P_Total = config->GetGiles_Var1(Marker_Tag); + T_Total = config->GetGiles_Var2(Marker_Tag); + FlowDir = config->GetGiles_FlowDir(Marker_Tag); + alphaIn_BC = atan(FlowDir[1]/FlowDir[0]); + + gammaIn_BC = 0; + if (nDim == 3){ + gammaIn_BC = FlowDir[2]; //atan(FlowDir[2]/FlowDir[0]); + } + + /*--- Non-dim. the inputs---*/ + P_Total /= config->GetPressure_Ref(); + T_Total /= config->GetTemperature_Ref(); + + /* --- Computes the total state --- */ + GetFluidModel()->SetTDState_PT(P_Total, T_Total); + Enthalpy_BC = GetFluidModel()->GetStaticEnergy()+ GetFluidModel()->GetPressure()/GetFluidModel()->GetDensity(); + Entropy_BC = GetFluidModel()->GetEntropy(); + + + /* --- Computes the inverse matrix R_c --- */ + conv_numerics->ComputeResJacobianGiles(GetFluidModel(), AveragePressure[val_marker][iSpan], AverageDensity[val_marker][iSpan], + AverageTurboVelocity[val_marker][iSpan], alphaIn_BC, gammaIn_BC, R_c, R_c_inv); + + GetFluidModel()->SetTDState_Prho(AveragePressure[val_marker][iSpan], AverageDensity[val_marker][iSpan]); + AverageEnthalpy = GetFluidModel()->GetStaticEnergy() + AveragePressure[val_marker][iSpan]/AverageDensity[val_marker][iSpan]; + AverageEntropy = GetFluidModel()->GetEntropy(); + + avgVel2 = 0.0; + for (iDim = 0; iDim < nDim; iDim++) avgVel2 += AverageVelocity[val_marker][iSpan][iDim]*AverageVelocity[val_marker][iSpan][iDim]; + if (nDim == 2){ + R[0] = -(AverageEntropy - Entropy_BC); + R[1] = -(AverageTurboVelocity[val_marker][iSpan][1] - tan(alphaIn_BC)*AverageTurboVelocity[val_marker][iSpan][0]); + R[2] = -(AverageEnthalpy + 0.5*avgVel2 - Enthalpy_BC); + } + + else{ + R[0] = -(AverageEntropy - Entropy_BC); + R[1] = -(AverageTurboVelocity[val_marker][iSpan][1] - tan(alphaIn_BC)*AverageTurboVelocity[val_marker][iSpan][0]); + R[2] = -(AverageTurboVelocity[val_marker][iSpan][2] - tan(gammaIn_BC)*AverageTurboVelocity[val_marker][iSpan][0]); + R[3] = -(AverageEnthalpy + 0.5*avgVel2 - Enthalpy_BC); + + } + /* --- Compute the avg component c_avg = R_c^-1 * R --- */ + for (iVar = 0; iVar < nVar-1; iVar++){ + c_avg[iVar] = 0.0; + for (jVar = 0; jVar < nVar-1; jVar++){ + c_avg[iVar] += R_c_inv[iVar][jVar]*R[jVar]; + } + } break; case TOTAL_CONDITIONS_PT_1D: - BC_Giles_Total_Inlet(conv_numerics, config, val_marker, c_avg, R, R_c_inv, R_c, iSpan, nSpanWiseSections); + /*--- Retrieve the specified total conditions for this inlet. ---*/ + P_Total = config->GetGiles_Var1(Marker_Tag); + T_Total = config->GetGiles_Var2(Marker_Tag); + FlowDir = config->GetGiles_FlowDir(Marker_Tag); + alphaIn_BC = atan(FlowDir[1]/FlowDir[0]); + + gammaIn_BC = 0; + if (nDim == 3){ + // Review definition of angle + gammaIn_BC = FlowDir[2]; //atan(FlowDir[2]/FlowDir[0]); + } + + /*--- Non-dim. the inputs---*/ + P_Total /= config->GetPressure_Ref(); + T_Total /= config->GetTemperature_Ref(); + + /* --- Computes the total state --- */ + GetFluidModel()->SetTDState_PT(P_Total, T_Total); + Enthalpy_BC = GetFluidModel()->GetStaticEnergy()+ GetFluidModel()->GetPressure()/GetFluidModel()->GetDensity(); + Entropy_BC = GetFluidModel()->GetEntropy(); + + + /* --- Computes the inverse matrix R_c --- */ + conv_numerics->ComputeResJacobianGiles(GetFluidModel(), AveragePressure[val_marker][iSpan], AverageDensity[val_marker][iSpan], + AverageTurboVelocity[val_marker][iSpan], alphaIn_BC, gammaIn_BC, R_c, R_c_inv); + + GetFluidModel()->SetTDState_Prho(AveragePressure[val_marker][nSpanWiseSections], AverageDensity[val_marker][nSpanWiseSections]); + AverageEnthalpy = GetFluidModel()->GetStaticEnergy() + AveragePressure[val_marker][nSpanWiseSections]/AverageDensity[val_marker][nSpanWiseSections]; + AverageEntropy = GetFluidModel()->GetEntropy(); + + + avgVel2 = 0.0; + for (iDim = 0; iDim < nDim; iDim++) avgVel2 += AverageVelocity[val_marker][iSpan][iDim]*AverageVelocity[val_marker][iSpan][iDim]; + if (nDim == 2){ + R[0] = -(AverageEntropy - Entropy_BC); + R[1] = -(AverageTurboVelocity[val_marker][nSpanWiseSections][1] - tan(alphaIn_BC)*AverageTurboVelocity[val_marker][nSpanWiseSections][0]); + R[2] = -(AverageEnthalpy + 0.5*avgVel2 - Enthalpy_BC); + } + + else{ + R[0] = -(AverageEntropy - Entropy_BC); + R[1] = -(AverageTurboVelocity[val_marker][nSpanWiseSections][1] - tan(alphaIn_BC)*AverageTurboVelocity[val_marker][nSpanWiseSections][0]); + R[2] = -(AverageTurboVelocity[val_marker][nSpanWiseSections][2] - tan(gammaIn_BC)*AverageTurboVelocity[val_marker][nSpanWiseSections][0]); + R[3] = -(AverageEnthalpy + 0.5*avgVel2 - Enthalpy_BC); + + } + /* --- Compute the avg component c_avg = R_c^-1 * R --- */ + for (iVar = 0; iVar < nVar-1; iVar++){ + c_avg[iVar] = 0.0; + for (jVar = 0; jVar < nVar-1; jVar++){ + c_avg[iVar] += R_c_inv[iVar][jVar]*R[jVar]; + } + } break; case MIXING_IN: case MIXING_OUT: - BC_Giles_Mixing(conv_numerics, val_marker, deltaprim, c_avg, iSpan); + /* --- Compute average jump of primitive at the mixing-plane interface--- */ + deltaprim[0] = ExtAverageDensity[val_marker][iSpan] - AverageDensity[val_marker][iSpan]; + deltaprim[1] = ExtAverageTurboVelocity[val_marker][iSpan][0] - AverageTurboVelocity[val_marker][iSpan][0]; + deltaprim[2] = ExtAverageTurboVelocity[val_marker][iSpan][1] - AverageTurboVelocity[val_marker][iSpan][1]; + if (nDim == 2){ + deltaprim[3] = ExtAveragePressure[val_marker][iSpan] - AveragePressure[val_marker][iSpan]; + } + else + { + deltaprim[3] = ExtAverageTurboVelocity[val_marker][iSpan][2] - AverageTurboVelocity[val_marker][iSpan][2]; + deltaprim[4] = ExtAveragePressure[val_marker][iSpan] - AveragePressure[val_marker][iSpan]; + } + + + /* --- Compute average jump of charachteristic variable at the mixing-plane interface--- */ + GetFluidModel()->SetTDState_Prho(AveragePressure[val_marker][iSpan], AverageDensity[val_marker][iSpan]); + AverageSoundSpeed = GetFluidModel()->GetSoundSpeed(); + conv_numerics->GetCharJump(AverageSoundSpeed, AverageDensity[val_marker][iSpan], deltaprim, c_avg); break; case MIXING_IN_1D: case MIXING_OUT_1D: - BC_Giles_Mixing(conv_numerics, val_marker, deltaprim, c_avg, nSpanWiseSections); + /* --- Compute average jump of primitive at the mixing-plane interface--- */ + deltaprim[0] = ExtAverageDensity[val_marker][nSpanWiseSections] - AverageDensity[val_marker][nSpanWiseSections]; + deltaprim[1] = ExtAverageTurboVelocity[val_marker][nSpanWiseSections][0] - AverageTurboVelocity[val_marker][nSpanWiseSections][0]; + deltaprim[2] = ExtAverageTurboVelocity[val_marker][nSpanWiseSections][1] - AverageTurboVelocity[val_marker][nSpanWiseSections][1]; + if (nDim == 2){ + deltaprim[3] = ExtAveragePressure[val_marker][nSpanWiseSections] - AveragePressure[val_marker][nSpanWiseSections]; + } + else + { + deltaprim[3] = ExtAverageTurboVelocity[val_marker][nSpanWiseSections][2] - AverageTurboVelocity[val_marker][nSpanWiseSections][2]; + deltaprim[4] = ExtAveragePressure[val_marker][nSpanWiseSections] - AveragePressure[val_marker][nSpanWiseSections]; + } + + /* --- Compute average jump of charachteristic variable at the mixing-plane interface--- */ + GetFluidModel()->SetTDState_Prho(AveragePressure[val_marker][nSpanWiseSections], AverageDensity[val_marker][nSpanWiseSections]); + AverageSoundSpeed = GetFluidModel()->GetSoundSpeed(); + conv_numerics->GetCharJump(AverageSoundSpeed, AverageDensity[val_marker][nSpanWiseSections], deltaprim, c_avg); break; @@ -6232,8 +6396,13 @@ void CEulerSolver::BC_Giles(CGeometry *geometry, CSolver **solver_container, CNu Pressure_e /= config->GetPressure_Ref(); /* --- Compute avg characteristic jump --- */ - c_avg[nDim+1] = -2.0*(AveragePressure[val_marker][iSpan]-Pressure_e); - + if (nDim == 2){ + c_avg[3] = -2.0*(AveragePressure[val_marker][iSpan]-Pressure_e); + } + else + { + c_avg[4] = -2.0*(AveragePressure[val_marker][iSpan]-Pressure_e); + } break; case STATIC_PRESSURE_1D: @@ -6241,8 +6410,13 @@ void CEulerSolver::BC_Giles(CGeometry *geometry, CSolver **solver_container, CNu Pressure_e /= config->GetPressure_Ref(); /* --- Compute avg characteristic jump --- */ - c_avg[nDim+1] = -2.0*(AveragePressure[val_marker][nSpanWiseSections]-Pressure_e); - + if (nDim == 2){ + c_avg[3] = -2.0*(AveragePressure[val_marker][nSpanWiseSections]-Pressure_e); + } + else + { + c_avg[4] = -2.0*(AveragePressure[val_marker][nSpanWiseSections]-Pressure_e); + } break; case RADIAL_EQUILIBRIUM: @@ -6253,68 +6427,79 @@ void CEulerSolver::BC_Giles(CGeometry *geometry, CSolver **solver_container, CNu break; - } + case MASS_FLOW_OUTLET: + auto const MassFlowRate_e = config->GetGiles_Var1(Marker_Tag); + auto const relFacMassFlowRate = config->GetGiles_Var2(Marker_Tag); + + Pressure_e = AveragePressure[val_marker][nSpanWiseSections]+relFacMassFlowRate*GetFluidModel()->GetdPdrho_e()*(AverageMassFlowRate[val_marker]-MassFlowRate_e); - auto const iZone = config->GetiZone(); + /*--- Compute avg characteristic jump ---*/ + c_avg[nDim + 1] = -2.0*(AveragePressure[val_marker][iSpan]-Pressure_e); + break; + + } /*--- Loop over all the vertices on this boundary marker ---*/ SU2_OMP_FOR_DYN(OMP_MIN_SIZE) - for (unsigned long iVertex = 0; iVertex < geometry->GetnVertexSpan(val_marker,iSpan); iVertex++) { + for (iVertex = 0; iVertex < geometry->GetnVertexSpan(val_marker,iSpan); iVertex++) { /*--- using the other vertex information for retrieving some information ---*/ - auto oldVertex = geometry->turbovertex[val_marker][iSpan][iVertex]->GetOldVertex(); - auto V_boundary= GetCharacPrimVar(val_marker, oldVertex); + oldVertex = geometry->turbovertex[val_marker][iSpan][iVertex]->GetOldVertex(); + V_boundary= GetCharacPrimVar(val_marker, oldVertex); /*--- Index of the closest interior node ---*/ - auto Point_Normal = geometry->vertex[val_marker][oldVertex]->GetNormal_Neighbor(); + Point_Normal = geometry->vertex[val_marker][oldVertex]->GetNormal_Neighbor(); /*--- Normal vector for this vertex (negate for outward convention), * this normal is scaled with the area of the face of the element ---*/ geometry->vertex[val_marker][oldVertex]->GetNormal(Normal); - for (unsigned short iDim = 0; iDim < nDim; iDim++) Normal[iDim] = -Normal[iDim]; + for (iDim = 0; iDim < nDim; iDim++) Normal[iDim] = -Normal[iDim]; conv_numerics->SetNormal(Normal); /*--- find the node related to the vertex ---*/ - auto iPoint = geometry->turbovertex[val_marker][iSpan][iVertex]->GetNode(); + iPoint = geometry->turbovertex[val_marker][iSpan][iVertex]->GetNode(); /*--- Normalize Normal vector for this vertex (already for outward convention) ---*/ geometry->turbovertex[val_marker][iSpan][iVertex]->GetNormal(UnitNormal); geometry->turbovertex[val_marker][iSpan][iVertex]->GetTurboNormal(turboNormal); /*--- Retrieve solution at this boundary node ---*/ - auto V_domain = nodes->GetPrimitive(iPoint); + V_domain = nodes->GetPrimitive(iPoint); /*--- Retrieve domain Secondary variables ---*/ - auto S_domain = nodes->GetSecondary(iPoint); + S_domain = nodes->GetSecondary(iPoint); /*--- Compute the internal state u_i ---*/ - su2double Velocity2_i = 0; - for (unsigned short iDim = 0; iDim < nDim; iDim++) + Velocity2_i = 0; + for (iDim = 0; iDim < nDim; iDim++) { Velocity_i[iDim] = nodes->GetVelocity(iPoint,iDim); Velocity2_i += Velocity_i[iDim]*Velocity_i[iDim]; } - auto Density_i = nodes->GetDensity(iPoint); + Density_i = nodes->GetDensity(iPoint); - auto Energy_i = nodes->GetEnergy(iPoint); - auto StaticEnergy_i = Energy_i - 0.5*Velocity2_i; + Energy_i = nodes->GetEnergy(iPoint); + StaticEnergy_i = Energy_i - 0.5*Velocity2_i; GetFluidModel()->SetTDState_rhoe(Density_i, StaticEnergy_i); - auto Pressure_i = GetFluidModel()->GetPressure(); + Pressure_i = GetFluidModel()->GetPressure(); ComputeTurboVelocity(Velocity_i, turboNormal, turboVelocity, config->GetMarker_All_TurbomachineryFlag(val_marker),config->GetKind_TurboMachinery(iZone)); - deltaprim[0] = Density_i - AverageDensity[val_marker][iSpan]; - deltaprim[1] = turboVelocity[0] - AverageTurboVelocity[val_marker][iSpan][0]; - deltaprim[2] = turboVelocity[1] - AverageTurboVelocity[val_marker][iSpan][1]; if (nDim == 2){ + deltaprim[0] = Density_i - AverageDensity[val_marker][iSpan]; + deltaprim[1] = turboVelocity[0] - AverageTurboVelocity[val_marker][iSpan][0]; + deltaprim[2] = turboVelocity[1] - AverageTurboVelocity[val_marker][iSpan][1]; deltaprim[3] = Pressure_i - AveragePressure[val_marker][iSpan]; } else{ + deltaprim[0] = Density_i - AverageDensity[val_marker][iSpan]; + deltaprim[1] = turboVelocity[0] - AverageTurboVelocity[val_marker][iSpan][0]; + deltaprim[2] = turboVelocity[1] - AverageTurboVelocity[val_marker][iSpan][1]; deltaprim[3] = turboVelocity[2] - AverageTurboVelocity[val_marker][iSpan][2]; deltaprim[4] = Pressure_i - AveragePressure[val_marker][iSpan]; } @@ -6324,13 +6509,8 @@ void CEulerSolver::BC_Giles(CGeometry *geometry, CSolver **solver_container, CNu conv_numerics->GetCharJump(AverageSoundSpeed, AverageDensity[val_marker][iSpan], deltaprim, cj); - auto pitch = geometry->GetMaxAngularCoord(val_marker, iSpan) - geometry->GetMinAngularCoord(val_marker,iSpan); - auto theta = geometry->turbovertex[val_marker][iSpan][iVertex]->GetRelAngularCoord(); - - auto const AvgMach = AverageTurboMach[0]*AverageTurboMach[0] + AverageTurboMach[1]*AverageTurboMach[1]; - auto Beta_inf= I*complex(sqrt(1.0 - AvgMach)); - su2double TwoPiThetaFreq_Pitch; - long freq; + pitch = geometry->GetMaxAngularCoord(val_marker, iSpan) - geometry->GetMinAngularCoord(val_marker,iSpan); + theta = geometry->turbovertex[val_marker][iSpan][iVertex]->GetRelAngularCoord(); switch(config->GetKind_Data_Giles(Marker_Tag)) { @@ -6340,18 +6520,11 @@ void CEulerSolver::BC_Giles(CGeometry *geometry, CSolver **solver_container, CNu case TOTAL_CONDITIONS_PT: case MIXING_IN:case TOTAL_CONDITIONS_PT_1D: case MIXING_IN_1D: if(config->GetSpatialFourier()){ - - /* --- Initial definition of variables ---*/ - auto c2js = complex(0.0,0.0); - auto c3js = complex(0.0,0.0); - auto c2ks = complex(0.0,0.0); - auto c3ks = complex(0.0,0.0); - auto c2js_Re = c2js.real(); - auto c3js_Re = c3js.real(); - su2double Beta_inf2; - if (AvgMach <= 1.0){ - for(unsigned long k=0; k < 2*kend_max+1; k++){ + Beta_inf= I*complex(sqrt(1.0 - AvgMach)); + c2js = complex(0.0,0.0); + c3js = complex(0.0,0.0); + for(k=0; k < 2*kend_max+1; k++){ freq = k - kend_max; if(freq >= (long)(-kend) && freq <= (long)(kend) && AverageTurboMach[0] > config->GetAverageMachLimit()){ TwoPiThetaFreq_Pitch = 2*PI_NUMBER*freq*theta/pitch; @@ -6370,10 +6543,16 @@ void CEulerSolver::BC_Giles(CGeometry *geometry, CSolver **solver_container, CNu c2js_Re = c2js.real(); c3js_Re = c3js.real(); - dcjs[0] = 0.0 - cj[0]; - dcjs[1] = c2js_Re - cj[1]; - dcjs[2] = 0.0 - cj[2]; - dcjs[nDim] = c3js_Re - cj[nDim]; // Overwrites previous value in 2D case + if (nDim == 2){ + dcjs[0] = 0.0 - cj[0]; + dcjs[1] = c2js_Re - cj[1]; + dcjs[2] = c3js_Re - cj[2]; + }else{ + dcjs[0] = 0.0 - cj[0]; + dcjs[1] = c2js_Re - cj[1]; + dcjs[2] = 0.0 - cj[2]; + dcjs[3] = c3js_Re - cj[3]; + } }else{ if (AverageTurboVelocity[val_marker][iSpan][1] >= 0.0){ @@ -6381,21 +6560,40 @@ void CEulerSolver::BC_Giles(CGeometry *geometry, CSolver **solver_container, CNu }else{ Beta_inf2= sqrt(AvgMach-1.0); } + if (nDim == 2){ + c2js_Re = -cj[3]*(Beta_inf2 + AverageTurboMach[1])/( 1.0 + AverageTurboMach[0]); + c3js_Re = cj[3]*(Beta_inf2 + AverageTurboMach[1])/( 1.0 + AverageTurboMach[0]); + c3js_Re *= (Beta_inf2 + AverageTurboMach[1])/( 1.0 + AverageTurboMach[0]); + }else{ + c2js_Re = -cj[4]*(Beta_inf2 + AverageTurboMach[1])/( 1.0 + AverageTurboMach[0]); + c3js_Re = cj[4]*(Beta_inf2 + AverageTurboMach[1])/( 1.0 + AverageTurboMach[0]); + c3js_Re *= (Beta_inf2 + AverageTurboMach[1])/( 1.0 + AverageTurboMach[0]); + } - c2js_Re = -cj[nDim+1]*(Beta_inf2 + AverageTurboMach[1])/( 1.0 + AverageTurboMach[0]); - c3js_Re = cj[nDim+1]*(Beta_inf2 + AverageTurboMach[1])/( 1.0 + AverageTurboMach[0]); - c3js_Re *= (Beta_inf2 + AverageTurboMach[1])/( 1.0 + AverageTurboMach[0]); - dcjs[0] = 0.0 - cj[0]; - dcjs[1] = c2js_Re - cj[1]; - dcjs[2] = 0.0 - cj[2]; - dcjs[nDim] = c3js_Re - cj[nDim]; + if (nDim == 2){ + dcjs[0] = 0.0 - cj[0]; + dcjs[1] = c2js_Re - cj[1]; + dcjs[2] = c3js_Re - cj[2]; + }else{ + dcjs[0] = 0.0 - cj[0]; + dcjs[1] = c2js_Re - cj[1]; + dcjs[2] = 0.0 - cj[2]; + dcjs[3] = c3js_Re - cj[3]; + } + } + } + else{ + if (nDim == 2){ + dcjs[0] = 0.0; + dcjs[1] = 0.0; + dcjs[2] = 0.0; + }else{ + dcjs[0] = 0.0; + dcjs[1] = 0.0; + dcjs[2] = 0.0; + dcjs[3] = 0.0; } - }else{ - dcjs[0] = 0.0; - dcjs[1] = 0.0; - dcjs[2] = 0.0; - dcjs[nDim] = 0.0; } /* --- Impose Inlet BC Reflecting--- */ delta_c[0] = relfacAvg*c_avg[0] + relfacFou*dcjs[0]; @@ -6410,11 +6608,10 @@ void CEulerSolver::BC_Giles(CGeometry *geometry, CSolver **solver_container, CNu break; - case STATIC_PRESSURE:case STATIC_PRESSURE_1D:case MIXING_OUT:case RADIAL_EQUILIBRIUM:case MIXING_OUT_1D: + case STATIC_PRESSURE:case STATIC_PRESSURE_1D:case MIXING_OUT:case RADIAL_EQUILIBRIUM:case MIXING_OUT_1D: case MASS_FLOW_OUTLET: /* --- implementation of Giles BC---*/ if(config->GetSpatialFourier()){ - su2double GilesBeta, cOutjs_Re; if (AvgMach > 1.0){ /* --- supersonic Giles implementation ---*/ if (AverageTurboVelocity[val_marker][iSpan][1] >= 0.0){ @@ -6438,9 +6635,9 @@ void CEulerSolver::BC_Giles(CGeometry *geometry, CSolver **solver_container, CNu }else{ /* --- subsonic Giles implementation ---*/ - auto cOutjs = complex(0.0,0.0); - auto cOutks = complex(0.0,0.0); - for(unsigned long k=0; k < 2*kend_max+1; k++){ + Beta_inf= I*complex(sqrt(1.0 - AvgMach)); + cOutjs = complex(0.0,0.0); + for(k=0; k < 2*kend_max+1; k++){ freq = k - kend_max; if(freq >= (long)(-kend) && freq <= (long)(kend) && AverageTurboMach[0] > config->GetAverageMachLimit()){ TwoPiThetaFreq_Pitch = 2*PI_NUMBER*freq*theta/pitch; @@ -6455,25 +6652,41 @@ void CEulerSolver::BC_Giles(CGeometry *geometry, CSolver **solver_container, CNu } cOutjs_Re = cOutjs.real(); - dcjs[nDim+1] = cOutjs_Re - cj[nDim+1]; + if (nDim == 2){ + dcjs[3] = cOutjs_Re - cj[3]; + } + else{ + dcjs[4] = cOutjs_Re - cj[4]; + } } } else{ - dcjs[nDim+1] = 0.0; + if (nDim == 2){ + dcjs[3] = 0.0; + } + else{ + dcjs[4] = 0.0; + } } /* --- Impose Outlet BC Non-Reflecting --- */ delta_c[0] = cj[0]; delta_c[1] = cj[1]; delta_c[2] = cj[2]; - delta_c[3] = cj[3]; - delta_c[nDim+1] = relfacAvg*c_avg[nDim+1] + relfacFou*dcjs[nDim+1]; + if (nDim == 2){ + delta_c[3] = relfacAvg*c_avg[3] + relfacFou*dcjs[3]; + } + else{ + delta_c[3] = cj[3]; + delta_c[4] = relfacAvg*c_avg[4] + relfacFou*dcjs[4]; + } + /*--- Automatically impose supersonic autoflow ---*/ if (abs(AverageTurboMach[0]) > 1.0000){ delta_c[0] = 0.0; delta_c[1] = 0.0; delta_c[2] = 0.0; - delta_c[3] = 0.0; + delta_c[2] = 0.0; if (nDim == 3)delta_c[4] = 0.0; } @@ -6485,25 +6698,31 @@ void CEulerSolver::BC_Giles(CGeometry *geometry, CSolver **solver_container, CNu } /*--- Compute primitive jump from characteristic variables ---*/ - for (unsigned short iVar = 0; iVar < nVar; iVar++) + for (iVar = 0; iVar < nVar; iVar++) { deltaprim[iVar]=0.0; - for (unsigned short jVar = 0; jVar < nVar; jVar++) + for (jVar = 0; jVar < nVar; jVar++) { deltaprim[iVar] += R_Matrix[iVar][jVar]*delta_c[jVar]; } } /*--- retrieve boundary variables ---*/ - su2double Density_b = AverageDensity[val_marker][iSpan] + deltaprim[0]; + Density_b = AverageDensity[val_marker][iSpan] + deltaprim[0]; turboVelocity[0] = AverageTurboVelocity[val_marker][iSpan][0] + deltaprim[1]; turboVelocity[1] = AverageTurboVelocity[val_marker][iSpan][1] + deltaprim[2]; - turboVelocity[nDim-1] = AverageTurboVelocity[val_marker][iSpan][nDim-1] + deltaprim[nDim]; - su2double Pressure_b = AveragePressure[val_marker][iSpan] + deltaprim[nDim+1]; + if(nDim == 2){ + Pressure_b = AveragePressure[val_marker][iSpan] + deltaprim[3]; + } + else{ + turboVelocity[2] = AverageTurboVelocity[val_marker][iSpan][2] + deltaprim[3]; + Pressure_b = AveragePressure[val_marker][iSpan] + deltaprim[4]; + } + ComputeBackVelocity(turboVelocity, turboNormal, Velocity_b, config->GetMarker_All_TurbomachineryFlag(val_marker), config->GetKind_TurboMachinery(iZone)); - su2double Velocity2_b = 0.0; - for (unsigned short iDim = 0; iDim < nDim; iDim++) { + Velocity2_b = 0.0; + for (iDim = 0; iDim < nDim; iDim++) { Velocity2_b+= Velocity_b[iDim]*Velocity_b[iDim]; } @@ -6511,20 +6730,20 @@ void CEulerSolver::BC_Giles(CGeometry *geometry, CSolver **solver_container, CNu Pressure_b = Pressure_i; Density_b = Density_i; Velocity2_b = 0.0; - for (unsigned short iDim = 0; iDim < nDim; iDim++) { + for (iDim = 0; iDim < nDim; iDim++) { Velocity_b[iDim] = Velocity_i[iDim]; Velocity2_b+= Velocity_b[iDim]*Velocity_b[iDim]; } } GetFluidModel()->SetTDState_Prho(Pressure_b, Density_b); - su2double Energy_b = GetFluidModel()->GetStaticEnergy() + 0.5*Velocity2_b; - su2double Temperature_b= GetFluidModel()->GetTemperature(); - su2double Enthalpy_b = Energy_b + Pressure_b/Density_b; + Energy_b = GetFluidModel()->GetStaticEnergy() + 0.5*Velocity2_b; + Temperature_b= GetFluidModel()->GetTemperature(); + Enthalpy_b = Energy_b + Pressure_b/Density_b; /*--- Primitive variables, using the derived quantities ---*/ V_boundary[0] = Temperature_b; - for (unsigned short iDim = 0; iDim < nDim; iDim++) + for (iDim = 0; iDim < nDim; iDim++) V_boundary[iDim+1] = Velocity_b[iDim]; V_boundary[nDim+1] = Pressure_b; V_boundary[nDim+2] = Density_b; @@ -6552,12 +6771,12 @@ void CEulerSolver::BC_Giles(CGeometry *geometry, CSolver **solver_container, CNu LinSysRes.AddBlock(iPoint, residual); /*--- Jacobian contribution for implicit integration ---*/ - if (config->GetKind_TimeIntScheme() == EULER_IMPLICIT) + if (implicit) Jacobian.AddBlock2Diag(iPoint, residual.jacobian_i); /*--- Viscous contribution ---*/ - if (config->GetViscous()) { + if (viscous) { /*--- Set laminar and eddy viscosity at the infinity ---*/ @@ -6611,7 +6830,7 @@ void CEulerSolver::BC_Giles(CGeometry *geometry, CSolver **solver_container, CNu /*--- Jacobian contribution for implicit integration ---*/ - if (config->GetKind_TimeIntScheme() == EULER_IMPLICIT) + if (implicit) Jacobian.SubtractBlock2Diag(iPoint, residual.jacobian_i); } @@ -6633,11 +6852,11 @@ void CEulerSolver::BC_Giles(CGeometry *geometry, CSolver **solver_container, CNu delete [] delta_c; delete [] deltaprim; delete [] cj; - for (unsigned short iVar = 0; iVar < nVar; iVar++) + for (iVar = 0; iVar < nVar; iVar++) { delete [] R_Matrix[iVar]; } - for (unsigned short iVar = 0; iVar < nVar-1; iVar++) + for (iVar = 0; iVar < nVar-1; iVar++) { delete [] R_c_inv[iVar]; delete [] R_c[iVar]; @@ -6655,73 +6874,6 @@ void CEulerSolver::BC_Giles(CGeometry *geometry, CSolver **solver_container, CNu delete [] turboVelocity; } -void CEulerSolver::BC_Giles_Total_Inlet( CNumerics *conv_numerics, CConfig *config, unsigned short val_marker, - su2double *&c_avg, su2double *&R, su2double **&R_c_inv, su2double **&R_c, unsigned short iSpan, unsigned short SpanwisePosition) { - - /*--- Calculates boundry condition for quasi-3D and 1D total boundary conditions ---*/ - - /*--- Retrieve the specified boundary conditions. ---*/ - auto const Marker_Tag = config->GetMarker_All_TagBound(val_marker); - auto P_Total = config->GetGiles_Var1(Marker_Tag); - auto T_Total = config->GetGiles_Var2(Marker_Tag); - auto const FlowDir = config->GetGiles_FlowDir(Marker_Tag); - auto const alphaIn_BC = atan(FlowDir[1]/FlowDir[0]); - - su2double gammaIn_BC = 0; - if (nDim == 3){ - gammaIn_BC = FlowDir[2]; //atan(FlowDir[2]/FlowDir[0]); - } - - /*--- Non-dim. the inputs---*/ - P_Total /= config->GetPressure_Ref(); - T_Total /= config->GetTemperature_Ref(); - - /* --- Computes the total state --- */ - GetFluidModel()->SetTDState_PT(P_Total, T_Total); - auto const Enthalpy_BC = GetFluidModel()->GetStaticEnergy()+ GetFluidModel()->GetPressure()/GetFluidModel()->GetDensity(); - auto const Entropy_BC = GetFluidModel()->GetEntropy(); - - /* --- Computes the inverse matrix R_c --- */ - conv_numerics->ComputeResJacobianGiles(GetFluidModel(), AveragePressure[val_marker][iSpan], AverageDensity[val_marker][iSpan], - AverageTurboVelocity[val_marker][iSpan], alphaIn_BC, gammaIn_BC, R_c, R_c_inv); - - GetFluidModel()->SetTDState_Prho(AveragePressure[val_marker][SpanwisePosition], AverageDensity[val_marker][SpanwisePosition]); - auto AverageEnthalpy = GetFluidModel()->GetStaticEnergy() + AveragePressure[val_marker][SpanwisePosition]/AverageDensity[val_marker][SpanwisePosition]; - auto AverageEntropy = GetFluidModel()->GetEntropy(); - - - su2double avgVel2 = 0.0; - for (unsigned short iDim = 0; iDim < nDim; iDim++) avgVel2 += AverageVelocity[val_marker][iSpan][iDim]*AverageVelocity[val_marker][iSpan][iDim]; - - R[0] = -(AverageEntropy - Entropy_BC); - R[1] = -(AverageTurboVelocity[val_marker][SpanwisePosition][1] - tan(alphaIn_BC)*AverageTurboVelocity[val_marker][SpanwisePosition][0]); - R[2] = -(AverageTurboVelocity[val_marker][SpanwisePosition][2] - tan(gammaIn_BC)*AverageTurboVelocity[val_marker][SpanwisePosition][0]); - R[nDim] = -(AverageEnthalpy + 0.5*avgVel2 - Enthalpy_BC); - - /* --- Compute the avg component c_avg = R_c^-1 * R --- */ - for (unsigned short iVar = 0; iVar < nVar-1; iVar++){ - c_avg[iVar] = 0.0; - for (unsigned short jVar = 0; jVar < nVar-1; jVar++){ - c_avg[iVar] += R_c_inv[iVar][jVar]*R[jVar]; - } - } -} - -void CEulerSolver::BC_Giles_Mixing(CNumerics *conv_numerics, unsigned short val_marker, su2double *&deltaprim, su2double *&c_avg, unsigned short SpanwisePosition) { - - /* --- Compute average jump of primitive at the mixing-plane interface--- */ - deltaprim[0] = ExtAverageDensity[val_marker][SpanwisePosition] - AverageDensity[val_marker][SpanwisePosition]; - deltaprim[1] = ExtAverageTurboVelocity[val_marker][SpanwisePosition][0] - AverageTurboVelocity[val_marker][SpanwisePosition][0]; - deltaprim[2] = ExtAverageTurboVelocity[val_marker][SpanwisePosition][1] - AverageTurboVelocity[val_marker][SpanwisePosition][1]; - deltaprim[3] = ExtAverageTurboVelocity[val_marker][SpanwisePosition][2] - AverageTurboVelocity[val_marker][SpanwisePosition][2]; - deltaprim[nDim+1] = ExtAveragePressure[val_marker][SpanwisePosition] - AveragePressure[val_marker][SpanwisePosition]; - - /* --- Compute average jump of charachteristic variable at the mixing-plane interface--- */ - GetFluidModel()->SetTDState_Prho(AveragePressure[val_marker][SpanwisePosition], AverageDensity[val_marker][SpanwisePosition]); - auto const AverageSoundSpeed = GetFluidModel()->GetSoundSpeed(); - conv_numerics->GetCharJump(AverageSoundSpeed, AverageDensity[val_marker][SpanwisePosition], deltaprim, c_avg); -} - void CEulerSolver::BC_Inlet(CGeometry *geometry, CSolver **solver_container, CNumerics *conv_numerics, CNumerics *visc_numerics, CConfig *config, unsigned short val_marker) { @@ -8947,8 +9099,7 @@ void CEulerSolver::TurboAverageProcess(CSolver **solver, CGeometry *geometry, CC /*--- Compute the averaged value for the boundary of interest for the span of interest ---*/ const bool belowMachLimit = (abs(MachTest)< config->GetAverageMachLimit()); - su2double avgDensity{0}, avgPressure{0}, avgKine{0}, avgOmega{0}, avgNu{0}, - avgVelocity[MAXNDIM] = {0}; + su2double avgDensity{0}, avgPressure{0}, avgKine{0}, avgOmega{0}, avgNu{0}, avgVelocity[MAXNDIM] = {0}; for (auto iVar = 0u; iVarGetMarker_All_Turbomachinery(iMarker) == iMarkerTP){ if (config->GetMarker_All_TurbomachineryFlag(iMarker) == marker_flag){ auto Marker_Tag = config->GetMarker_All_TagBound(iMarker); - if(config->GetBoolGiles() || config->GetBoolRiemann()){ + if(config->GetMarker_All_Giles(iMarker) || config->GetBoolRiemann()){ // May have to implement something similar for Riemann? if(config->GetBoolRiemann()){ if(config->GetKind_Data_Riemann(Marker_Tag) == RADIAL_EQUILIBRIUM){ RadialEquilibriumPressure[iMarker][nSpanWiseSections/2] = config->GetRiemann_Var1(Marker_Tag)/config->GetPressure_Ref(); diff --git a/SU2_CFD/src/solvers/CIncEulerSolver.cpp b/SU2_CFD/src/solvers/CIncEulerSolver.cpp index 5b81394bd9a..215445d2c5e 100644 --- a/SU2_CFD/src/solvers/CIncEulerSolver.cpp +++ b/SU2_CFD/src/solvers/CIncEulerSolver.cpp @@ -325,7 +325,7 @@ void CIncEulerSolver::SetNondimensionalization(CConfig *config, unsigned short i config->SetGas_Constant(UNIVERSAL_GAS_CONSTANT / (config->GetMolecular_Weight() / 1000.0)); Pressure_Thermodynamic = config->GetPressure_Thermodynamic(); - auxFluidModel = new CFluidScalar(config->GetSpecific_Heat_Cp(), config->GetGas_Constant(), Pressure_Thermodynamic, config); + auxFluidModel = new CFluidScalar(Pressure_Thermodynamic, config); auxFluidModel->SetTDState_T(Temperature_FreeStream, config->GetSpecies_Init()); break; @@ -482,7 +482,7 @@ void CIncEulerSolver::SetNondimensionalization(CConfig *config, unsigned short i break; case FLUID_MIXTURE: - fluidModel = new CFluidScalar(Specific_Heat_CpND, Gas_ConstantND, Pressure_ThermodynamicND, config); + fluidModel = new CFluidScalar(Pressure_ThermodynamicND, config); fluidModel->SetTDState_T(Temperature_FreeStreamND, config->GetSpecies_Init()); break; diff --git a/SU2_CFD/src/solvers/CSolver.cpp b/SU2_CFD/src/solvers/CSolver.cpp index e2f5927190d..302aff35560 100644 --- a/SU2_CFD/src/solvers/CSolver.cpp +++ b/SU2_CFD/src/solvers/CSolver.cpp @@ -1730,6 +1730,7 @@ void CSolver::AdaptCFLNumber(CGeometry **geometry, CSolver *solverFlow = solver_container[iMesh][FLOW_SOL]; CSolver *solverTurb = solver_container[iMesh][TURB_SOL]; + CSolver *solverSpecies = solver_container[iMesh][SPECIES_SOL]; /* Compute the reduction factor for CFLs on the coarse levels. */ @@ -1744,10 +1745,12 @@ void CSolver::AdaptCFLNumber(CGeometry **geometry, solver residual within the specified number of linear iterations. */ su2double linResTurb = 0.0; + su2double linResSpecies = 0.0; if ((iMesh == MESH_0) && solverTurb) linResTurb = solverTurb->GetResLinSolver(); + if ((iMesh == MESH_0) && solverSpecies) linResSpecies = solverSpecies->GetResLinSolver(); - /* Max linear residual between flow and turbulence. */ - const su2double linRes = max(solverFlow->GetResLinSolver(), linResTurb); + /* Max linear residual between flow and turbulence/species transport. */ + const su2double linRes = max(solverFlow->GetResLinSolver(), max(linResTurb, linResSpecies)); /* Tolerance limited to an acceptable value. */ const su2double linTol = max(acceptableLinTol, config->GetLinear_Solver_Error()); @@ -1779,6 +1782,11 @@ void CSolver::AdaptCFLNumber(CGeometry **geometry, New_Func += log10(solverTurb->GetRes_RMS(iVar)); } } + if ((iMesh == MESH_0) && solverSpecies) { + for (unsigned short iVar = 0; iVar < solverSpecies->GetnVar(); iVar++) { + New_Func += log10(solverSpecies->GetRes_RMS(iVar)); + } + } /* Compute the difference in the nonlinear residuals between the current and previous iterations, taking care with very low initial @@ -1822,6 +1830,7 @@ void CSolver::AdaptCFLNumber(CGeometry **geometry, su2double myCFLMin = 1e30, myCFLMax = 0.0, myCFLSum = 0.0; const su2double CFLTurbReduction = config->GetCFLRedCoeff_Turb(); + const su2double CFLSpeciesReduction = config->GetCFLRedCoeff_Species(); SU2_OMP_MASTER if ((iMesh == MESH_0) && fullComms) { @@ -1888,6 +1897,9 @@ void CSolver::AdaptCFLNumber(CGeometry **geometry, if ((iMesh == MESH_0) && solverTurb) { solverTurb->GetNodes()->SetLocalCFL(iPoint, CFL * CFLTurbReduction); } + if ((iMesh == MESH_0) && solverSpecies) { + solverSpecies->GetNodes()->SetLocalCFL(iPoint, CFL * CFLSpeciesReduction); + } /* Store min and max CFL for reporting on the fine grid. */ diff --git a/TestCases/hybrid_regression.py b/TestCases/hybrid_regression.py index c893649c1dd..af200ef7603 100644 --- a/TestCases/hybrid_regression.py +++ b/TestCases/hybrid_regression.py @@ -561,7 +561,7 @@ def main(): Jones_tc_restart.cfg_dir = "turbomachinery/APU_turbocharger" Jones_tc_restart.cfg_file = "Jones_restart.cfg" Jones_tc_restart.test_iter = 5 - Jones_tc_restart.test_vals = [-6.614627, -3.001324, -14.336143, -8.776079, -11.382919, -5.852328, 73273.000000, 73273.000000, 0.019884, 82.491000] + Jones_tc_restart.test_vals = [-10.467026, -2.871699, -19.214627, -13.508254, -11.582396, -6.306163, 73273, 73273, 0.019884, 82.491] test_list.append(Jones_tc_restart) # 2D axial stage @@ -582,6 +582,15 @@ def main(): transonic_stator_restart.test_vals_aarch64 = [-5.007735, -3.099310, -2.751696, 1.091966, -3.542819, 2.163237, -471630, 94.866, -0.035738] test_list.append(transonic_stator_restart) + # Multiple turbomachinery interface restart + multi_interface = TestCase('multi_interface') + multi_interface.cfg_dir = "turbomachinery/multi_interface" + multi_interface.cfg_file = "multi_interface_rst.cfg" + multi_interface.test_iter = 5 + multi_interface.test_vals = [-8.632374, -8.895124, -9.350417] + multi_interface.test_vals_aarch64 = [-8.632374, -8.895124, -9.350417] + test_list.append(multi_interface) + ###################################### ### Sliding Mesh ### ###################################### diff --git a/TestCases/parallel_regression.py b/TestCases/parallel_regression.py index 03e6fbaef47..4047b6b39ba 100644 --- a/TestCases/parallel_regression.py +++ b/TestCases/parallel_regression.py @@ -1083,7 +1083,7 @@ def main(): Aachen_3D_restart.cfg_file = "aachen_3D_MP_restart.cfg" Aachen_3D_restart.test_iter = 5 Aachen_3D_restart.enabled_with_asan = False - Aachen_3D_restart.test_vals = [-9.853215, -8.891480, -9.610418, -8.028579, -7.792957, -4.384186] + Aachen_3D_restart.test_vals = [-9.829186, -8.875103, -9.609509, -8.075211, -7.759491, -4.360714] test_list.append(Aachen_3D_restart) # Jones APU Turbocharger restart @@ -1091,7 +1091,7 @@ def main(): Jones_tc_restart.cfg_dir = "turbomachinery/APU_turbocharger" Jones_tc_restart.cfg_file = "Jones_restart.cfg" Jones_tc_restart.test_iter = 5 - Jones_tc_restart.test_vals = [-6.614623, -3.001323, -14.336147, -8.776081, -11.382919, -5.852327, 73273, 73273, 0.019884, 82.491] + Jones_tc_restart.test_vals = [-10.467612, -2.871708, -19.345651, -13.625871, -11.582397, -6.306168, 73273, 73273, 0.019884, 82.491] test_list.append(Jones_tc_restart) # 2D axial stage @@ -1112,6 +1112,15 @@ def main(): transonic_stator_restart.test_vals_aarch64 = [-5.011834, -3.091110, -2.757795, 1.087934, -3.544707, 2.166101, -471630, 94.868, -0.035888] test_list.append(transonic_stator_restart) + # Multiple turbomachinery interface restart + multi_interface = TestCase('multi_interface') + multi_interface.cfg_dir = "turbomachinery/multi_interface" + multi_interface.cfg_file = "multi_interface_rst.cfg" + multi_interface.test_iter = 5 + multi_interface.test_vals = [-8.632374, -8.895124, -9.350417] + multi_interface.test_vals_aarch64 = [-8.632374, -8.895124, -9.350417] + test_list.append(multi_interface) + ###################################### ### Sliding Mesh ### ###################################### @@ -1339,6 +1348,7 @@ def main(): sp_pinArray_3d_cht_mf_hf_tp.multizone = True test_list.append(sp_pinArray_3d_cht_mf_hf_tp) + ########################## ### Python wrapper ### ########################## diff --git a/TestCases/serial_regression.py b/TestCases/serial_regression.py index 080b41053b7..08b45e3b3ff 100644 --- a/TestCases/serial_regression.py +++ b/TestCases/serial_regression.py @@ -860,7 +860,7 @@ def main(): Aachen_3D_restart.cfg_dir = "turbomachinery/Aachen_turbine" Aachen_3D_restart.cfg_file = "aachen_3D_MP_restart.cfg" Aachen_3D_restart.test_iter = 5 - Aachen_3D_restart.test_vals = [-9.853207, -8.891479, -9.610411, -8.028566, -7.792947, -4.384177] + Aachen_3D_restart.test_vals = [-9.829185, -8.875103, -9.609505, -8.075194, -7.759490, -4.360713] test_list.append(Aachen_3D_restart) # Jones APU Turbocharger restart @@ -868,7 +868,7 @@ def main(): Jones_tc_restart.cfg_dir = "turbomachinery/APU_turbocharger" Jones_tc_restart.cfg_file = "Jones_restart.cfg" Jones_tc_restart.test_iter = 5 - Jones_tc_restart.test_vals = [-6.614623, -3.001323, -14.336147, -8.776081, -11.382919, -5.852327, 73273.000000, 73273.000000, 0.019884, 82.491000] + Jones_tc_restart.test_vals = [-10.466644, -2.871703, -19.193464, -13.487820, -11.582397, -6.306167, 73273, 73273, 0.019884, 82.491] test_list.append(Jones_tc_restart) # 2D axial stage @@ -889,6 +889,15 @@ def main(): transonic_stator_restart.test_vals_aarch64 = [-5.008547, -3.102420, -2.752033, 1.091152, -3.543849, 2.169844, -471630, 94.866, -0.035806] test_list.append(transonic_stator_restart) + # Multiple turbomachinery interface restart + multi_interface = TestCase('multi_interface') + multi_interface.cfg_dir = "turbomachinery/multi_interface" + multi_interface.cfg_file = "multi_interface_rst.cfg" + multi_interface.test_iter = 5 + multi_interface.test_vals = [-8.632374, -8.895124, -9.350417] + multi_interface.test_vals_aarch64 = [-8.632374, -8.895124, -9.350417] + test_list.append(multi_interface) + ###################################### ### Sliding Mesh ### diff --git a/TestCases/turbomachinery/multi_interface/multi_interface_rst.cfg b/TestCases/turbomachinery/multi_interface/multi_interface_rst.cfg new file mode 100644 index 00000000000..a52301e7153 --- /dev/null +++ b/TestCases/turbomachinery/multi_interface/multi_interface_rst.cfg @@ -0,0 +1,151 @@ +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% % +% SU2 configuration file % +% Case description: Quarter Concentirc Cylinder Channel w/ % +% Different Interface Methods % +% Author: J. Kelly % +% Institution: University of Liverpool % +% Date: Sep 5th, 2024 % +% File Version 8.0.1 "Harrier" % +% % +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% +MULTIZONE= YES +CONFIG_LIST=(zone_1.cfg,zone_2.cfg,zone_3.cfg) +NZONES=3 +% +% ------------- DIRECT, ADJOINT, AND LINEARIZED PROBLEM DEFINITION ------------% +% +SOLVER= RANS +KIND_TURB_MODEL= SA +MATH_PROBLEM= DIRECT +RESTART_SOL= YES +% +% -------------------- COMPRESSIBLE FREE-STREAM DEFINITION --------------------% +% +MACH_NUMBER= 0.35 +AOA= 0.0 +FREESTREAM_PRESSURE= 100.0E+03 +FREESTREAM_TEMPERATURE= 300.0 +FREESTREAM_DENSITY= 1.205 +FREESTREAM_OPTION= TEMPERATURE_FS +FREESTREAM_TURBULENCEINTENSITY = 0.05 +FREESTREAM_TURB2LAMVISCRATIO = 100.0 +INIT_OPTION= TD_CONDITIONS +% +% ---------------------- REFERENCE VALUE DEFINITION ---------------------------% +% +REF_ORIGIN_MOMENT_X = 0.00 +REF_ORIGIN_MOMENT_Y = 0.00 +REF_ORIGIN_MOMENT_Z = 0.00 +REF_LENGTH= 1.0 +REF_AREA= 1.0 +REF_DIMENSIONALIZATION= DIMENSIONAL +% +% ------------------------------ EQUATION OF STATE ----------------------------% +% +FLUID_MODEL= STANDARD_AIR +% +% --------------------------- VISCOSITY MODEL ---------------------------------% +% +VISCOSITY_MODEL= SUTHERLAND +MU_REF= 1.716E-5 +MU_T_REF= 273.15 +SUTHERLAND_CONSTANT= 110.4 +% +% --------------------------- THERMAL CONDUCTIVITY MODEL ----------------------% +% +CONDUCTIVITY_MODEL= CONSTANT_PRANDTL +% +% -------------------- BOUNDARY CONDITION DEFINITION --------------------------% +% +MARKER_HEATFLUX= (HUB1, 0.0, SHROUD1, 0.0, HUB2, 0.0, SHROUD2, 0.0, HUB3, 0.0, SHROUD3, 0.0) +MARKER_PERIODIC= (PER1, PER2, 0.0, 0.0, 0.0, 0.0, 0.0, 90.0, 0.0, 0.0, 0.0, PER4, PER3, 0.0, 0.0, 0.0, 0.0, 0.0, 90.0, 0.0, 0.0, 0.0, PER6, PER5, 0.0, 0.0, 0.0, 0.0, 0.0, 90.0, 0.0, 0.0, 0.0) +% +%-------- INFLOW/OUTFLOW BOUNDARY CONDITION SPECIFIC FOR TURBOMACHINERY --------% +% +MARKER_TURBOMACHINERY= (INFLOW, OUTMIX_0, INMIX_1, OUTMIX_1, INMIX_2, OUTFLOW) +MARKER_ZONE_INTERFACE= (OUTMIX_0, INMIX_1, OUTMIX_1, INMIX_2) +MARKER_FLUID_INTERFACE= (OUTMIX_0, INMIX_1) +MARKER_MIXINGPLANE_INTERFACE= (OUTMIX_1, INMIX_2) +MARKER_GILES= (INFLOW, TOTAL_CONDITIONS_PT, 104E+03, 300, 1.0, 0.0, 0.0, 1.0, 0.0, OUTMIX_1, MIXING_OUT, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, INMIX_2, MIXING_IN, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, OUTFLOW, STATIC_PRESSURE_1D, 90.0E+03, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0) +GILES_EXTRA_RELAXFACTOR= (0.05, 0.05) +SPATIAL_FOURIER= NO +% +%---------------------------- TURBOMACHINERY SIMULATION -----------------------------% +% +TURBOMACHINERY_KIND= AXIAL, AXIAL, AXIAL +TURBO_PERF_KIND= (TURBINE, TURBINE, TURBINE) +MIXINGPLANE_INTERFACE_KIND= LINEAR_INTERPOLATION +TURBULENT_MIXINGPLANE= YES +RAMP_OUTLET_PRESSURE= NO +AVERAGE_PROCESS_KIND= MIXEDOUT +PERFORMANCE_AVERAGE_PROCESS_KIND= MIXEDOUT +MIXEDOUT_COEFF= (1.0, 1.0E-05, 100) +AVERAGE_MACH_LIMIT= 0.05 +% +% ------------- COMMON PARAMETERS DEFINING THE NUMERICAL METHOD ---------------% +% +NUM_METHOD_GRAD= WEIGHTED_LEAST_SQUARES +CFL_NUMBER= 2.0 +CFL_ADAPT= YES +CFL_ADAPT_PARAM= ( 0.99, 1.01, 1.0, 20, 1E-4) +% +% ------------------------ LINEAR SOLVER DEFINITION ---------------------------% +% +LINEAR_SOLVER= FGMRES +LINEAR_SOLVER_PREC= LU_SGS +LINEAR_SOLVER_ERROR= 1E-4 +LINEAR_SOLVER_ITER= 5 +% +% ----------------------- SLOPE LIMITER DEFINITION ----------------------------% +% +VENKAT_LIMITER_COEFF= 0.01 +LIMITER_ITER= 999999 +% +% -------------------- FLOW NUMERICAL METHOD DEFINITION -----------------------% +% +CONV_NUM_METHOD_FLOW= ROE +MUSCL_FLOW= NO +SLOPE_LIMITER_FLOW= VENKATAKRISHNAN_WANG +TIME_DISCRE_FLOW= EULER_IMPLICIT +% +% -------------------- TURBULENT NUMERICAL METHOD DEFINITION ------------------% +% +CONV_NUM_METHOD_TURB= SCALAR_UPWIND +MUSCL_TURB= NO +SLOPE_LIMITER_TURB= VENKATAKRISHNAN +TIME_DISCRE_TURB= EULER_IMPLICIT +CFL_REDUCTION_TURB= 1.0 +% +% --------------------------- CONVERGENCE PARAMETERS --------------------------% +% +OUTER_ITER= 10 +CONV_FIELD= RMS_DENSITY +CONV_RESIDUAL_MINVAL= -10 +CONV_STARTITER= 0 +CONV_CAUCHY_ELEMS= 999 +CONV_CAUCHY_EPS= 1E-6 +SCREEN_OUTPUT= (OUTER_ITER, RMS_DENSITY[0], RMS_DENSITY[1], RMS_DENSITY[2]) +OUTPUT_FILES= RESTART +% +% ------------------------- INPUT/OUTPUT INFORMATION --------------------------% +% +MESH_FILENAME= three_zone_quarter_cyl.su2 +MESH_FORMAT= SU2 +MESH_OUT_FILENAME= mesh_out.su2 +SOLUTION_FILENAME= solution_flow.dat +SOLUTION_ADJ_FILENAME= solution_adj.dat +TABULAR_FORMAT= TECPLOT +CONV_FILENAME= history +RESTART_FILENAME= restart_flow.dat +RESTART_ADJ_FILENAME= restart_adj.dat +VOLUME_FILENAME= flow +VOLUME_ADJ_FILENAME= adjoint +GRAD_OBJFUNC_FILENAME= of_grad.dat +SURFACE_FILENAME= surface_flow +SURFACE_ADJ_FILENAME= surface_adjoint +OUTPUT_WRT_FREQ= 100 +HISTORY_WRT_FREQ_OUTER= 1 +WRT_ZONE_CONV= NO +WRT_ZONE_HIST= YES diff --git a/TestCases/turbomachinery/multi_interface/zone_1.cfg b/TestCases/turbomachinery/multi_interface/zone_1.cfg new file mode 100644 index 00000000000..3d7a019754a --- /dev/null +++ b/TestCases/turbomachinery/multi_interface/zone_1.cfg @@ -0,0 +1,6 @@ +% ----------------------- DYNAMIC MESH DEFINITION -----------------------------% +% +% +% Type of dynamic mesh (NONE, ROTATING_FRAME) +GRID_MOVEMENT= NONE +% diff --git a/TestCases/turbomachinery/multi_interface/zone_2.cfg b/TestCases/turbomachinery/multi_interface/zone_2.cfg new file mode 100644 index 00000000000..3d7a019754a --- /dev/null +++ b/TestCases/turbomachinery/multi_interface/zone_2.cfg @@ -0,0 +1,6 @@ +% ----------------------- DYNAMIC MESH DEFINITION -----------------------------% +% +% +% Type of dynamic mesh (NONE, ROTATING_FRAME) +GRID_MOVEMENT= NONE +% diff --git a/TestCases/turbomachinery/multi_interface/zone_3.cfg b/TestCases/turbomachinery/multi_interface/zone_3.cfg new file mode 100644 index 00000000000..3d7a019754a --- /dev/null +++ b/TestCases/turbomachinery/multi_interface/zone_3.cfg @@ -0,0 +1,6 @@ +% ----------------------- DYNAMIC MESH DEFINITION -----------------------------% +% +% +% Type of dynamic mesh (NONE, ROTATING_FRAME) +GRID_MOVEMENT= NONE +% diff --git a/TestCases/tutorials.py b/TestCases/tutorials.py index 5d987f8cd9a..b4c8534d8a1 100644 --- a/TestCases/tutorials.py +++ b/TestCases/tutorials.py @@ -64,6 +64,15 @@ def main(): cht_incompressible.multizone = True test_list.append(cht_incompressible) + # Solid-to-solid and solid-to-fluid CHT with contact resistance + cht_CR = TestCase('cht_solid_solid') + cht_CR.cfg_dir = "../Tutorials/multiphysics/contact_resistance_cht" + cht_CR.cfg_file = "master.cfg" + cht_CR.test_iter = 80 + cht_CR.test_vals = [ -8.857438, -9.377593, -10.097769, -2.122358] + cht_CR.multizone = True + test_list.append(cht_CR) + ### Incompressible Flow # 2D pin case massflow periodic with heatflux BC and prescribed extracted outlet heat diff --git a/config_template.cfg b/config_template.cfg index 00334732b30..4044a611957 100644 --- a/config_template.cfg +++ b/config_template.cfg @@ -672,6 +672,9 @@ FAN_POLY_EFF= 1.0 % Only half engine is in the computational grid (NO, YES) ENGINE_HALF_MODEL= NO % +% Evaluation Frequency for Engine and Actuator Disk Markers +BC_EVAL_FREQ = 40 +% % Damping factor for the engine inflow. DAMP_ENGINE_INFLOW= 0.95 %