From 3801e4db3ed558647c6227180e177771fc0202b3 Mon Sep 17 00:00:00 2001 From: Ole Burghardt Date: Wed, 14 Aug 2024 18:17:14 +0200 Subject: [PATCH 01/22] Throw error message in case the multizone discrete adjoint solver will write restart files during runtime (there is a bug, potentially when clearing the tape). --- Common/src/CConfig.cpp | 11 ++++++++++- 1 file changed, 10 insertions(+), 1 deletion(-) diff --git a/Common/src/CConfig.cpp b/Common/src/CConfig.cpp index ec4a555d6af..7aa9d7af84a 100644 --- a/Common/src/CConfig.cpp +++ b/Common/src/CConfig.cpp @@ -3391,8 +3391,17 @@ void CConfig::SetPostprocessing(SU2_COMPONENT val_software, unsigned short val_i /*--- Using default frequency of 250 for all files when steady, and 1 for unsteady. ---*/ for (auto iVolumeFreq = 0; iVolumeFreq < nVolumeOutputFrequencies; iVolumeFreq++){ - VolumeOutputFrequencies[iVolumeFreq] = Time_Domain ? 1 : 250; + if (Multizone_Problem && DiscreteAdjoint) { + VolumeOutputFrequencies[iVolumeFreq] = nOuterIter; + } + else { + VolumeOutputFrequencies[iVolumeFreq] = Time_Domain ? 1 : 250; + } } + } else if (Multizone_Problem && DiscreteAdjoint) { + SU2_MPI::Error(string("OUTPUT_WRT_FREQ cannot be specified for this solver " + "(writing of restart and sensitivity files not possible for multizone discrete adjoint during runtime yet).\n" + "Please remove this option from the config file, output files will be written when solver finalizes.\n"), CURRENT_FUNCTION); } else if (nVolumeOutputFrequencies < nVolumeOutputFiles) { /*--- If there are fewer frequencies than files, repeat the last frequency. * This is useful to define 1 frequency for the restart file and 1 frequency for all the visualization files. ---*/ From 0e5f95b75d780a317fe3aef65c11ff17883009e8 Mon Sep 17 00:00:00 2001 From: Ole Burghardt Date: Fri, 16 Aug 2024 00:22:57 +0200 Subject: [PATCH 02/22] Softcode indentifier/index. --- Common/include/basic_types/ad_structure.hpp | 26 ++++++++++++++----- Common/include/geometry/dual_grid/CPoint.hpp | 5 ++-- Common/src/geometry/dual_grid/CPoint.cpp | 4 +-- .../drivers/CDiscAdjMultizoneDriver.hpp | 2 +- SU2_CFD/include/variables/CVariable.hpp | 12 ++++----- .../src/drivers/CDiscAdjMultizoneDriver.cpp | 2 +- SU2_CFD/src/solvers/CSolver.cpp | 2 +- SU2_CFD/src/variables/CVariable.cpp | 4 +-- externals/codi | 2 +- 9 files changed, 35 insertions(+), 24 deletions(-) diff --git a/Common/include/basic_types/ad_structure.hpp b/Common/include/basic_types/ad_structure.hpp index c44374f7eb7..ecd2c93b1c1 100644 --- a/Common/include/basic_types/ad_structure.hpp +++ b/Common/include/basic_types/ad_structure.hpp @@ -38,6 +38,9 @@ */ namespace AD { #ifndef CODI_REVERSE_TYPE + +using Identifier = int; + /*! * \brief Start the recording of the operations and involved variables. * If called, the computational graph of all operations occuring after the call will be stored, @@ -101,14 +104,20 @@ inline void EndUseAdjoints() {} * \param[in] index - Position in the adjoint vector. * \param[in] val - adjoint value to be set. */ -inline void SetDerivative(int index, const double val) {} +inline void SetDerivative(Identifier index, const double val) {} /*! * \brief Extracts the adjoint value at index * \param[in] index - position in the adjoint vector where the derivative will be extracted. * \return Derivative value. */ -inline double GetDerivative(int index) { return 0.0; } +inline double GetDerivative(Identifier index) { return 0.0; } + +/*! + * \brief Returns the identifier that represents an inactive variable. + * \return Passive index. + */ +inline Identifier GetPassiveIndex() { return 0; } /*! * \brief Clears the currently stored adjoints but keeps the computational graph. @@ -259,7 +268,7 @@ inline void SetExtFuncOut(T&& data, const int size_x, const int size_y) {} * \param[in] data - variable whose gradient information will be extracted. * \param[in] index - where obtained gradient information will be stored. */ -inline void SetIndex(int& index, const su2double& data) {} +inline void SetIndex(Identifier& index, const su2double& data) {} /*! * \brief Pushes back the current tape position to the tape position's vector. @@ -304,6 +313,7 @@ inline void EndNoSharedReading() {} using CheckpointHandler = codi::ExternalFunctionUserData; using Tape = su2double::Tape; +using Identifier = su2double::Identifier; #ifdef HAVE_OPDI using ExtFuncHelper = codi::OpenMPExternalFunctionHelper; @@ -470,14 +480,14 @@ FORCEINLINE void BeginUseAdjoints() { AD::getTape().beginUseAdjointVector(); } FORCEINLINE void EndUseAdjoints() { AD::getTape().endUseAdjointVector(); } -FORCEINLINE void SetIndex(int& index, const su2double& data) { index = data.getIdentifier(); } +FORCEINLINE void SetIndex(Identifier& index, const su2double& data) { index = data.getIdentifier(); } // 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. // It should be safeguarded by calls to AD::BeginUseAdjoints() and AD::EndUseAdjoints(). -FORCEINLINE void SetDerivative(int index, const double val) { - if (index == 0) // Allow multiple threads to "set the derivative" of passive variables without causing data races. +FORCEINLINE void SetDerivative(Identifier index, const double val) { + if (!AD::getTape().isIdentifierActive(index)) // Allow multiple threads to "set the derivative" of passive variables without causing data races. return; AD::getTape().setGradient(index, val, codi::AdjointsManagement::Manual); @@ -488,10 +498,12 @@ FORCEINLINE void SetDerivative(int index, const double val) { // Otherwise, please ensure sufficient adjoint vector size by a call to AD::ResizeAdjoints(). // This method does not perform locking either. // It should be safeguarded by calls to AD::BeginUseAdjoints() and AD::EndUseAdjoints(). -FORCEINLINE double GetDerivative(int index) { +FORCEINLINE double GetDerivative(Identifier index) { return AD::getTape().getGradient(index, codi::AdjointsManagement::Manual); } +FORCEINLINE Identifier GetPassiveIndex() { return AD::getTape().getPassiveIndex(); } + FORCEINLINE bool IsIdentifierActive(su2double const& value) { return getTape().isIdentifierActive(value.getIdentifier()); } diff --git a/Common/include/geometry/dual_grid/CPoint.hpp b/Common/include/geometry/dual_grid/CPoint.hpp index c297b1b5df9..acae3cd172d 100644 --- a/Common/include/geometry/dual_grid/CPoint.hpp +++ b/Common/include/geometry/dual_grid/CPoint.hpp @@ -113,9 +113,8 @@ class CPoint { su2activevector MaxLength; /*!< \brief The maximum cell-center to cell-center length. */ su2activevector RoughnessHeight; /*!< \brief Roughness of the nearest wall. */ - su2matrix AD_InputIndex; /*!< \brief Indices of Coord variables in the adjoint vector. */ - su2matrix - AD_OutputIndex; /*!< \brief Indices of Coord variables in the adjoint vector after having been updated. */ + su2matrix AD_InputIndex; /*!< \brief Indices of Coord variables in the adjoint vector. */ + su2matrix AD_OutputIndex; /*!< \brief Indices of Coord variables in the adjoint vector after having been updated. */ /*! * \brief Allocate fields required by the minimal constructor. diff --git a/Common/src/geometry/dual_grid/CPoint.cpp b/Common/src/geometry/dual_grid/CPoint.cpp index f11b7dd8945..9753d71f645 100644 --- a/Common/src/geometry/dual_grid/CPoint.cpp +++ b/Common/src/geometry/dual_grid/CPoint.cpp @@ -67,8 +67,8 @@ void CPoint::FullAllocation(unsigned short imesh, const CConfig* config) { } if (config->GetDiscrete_Adjoint()) { - AD_InputIndex.resize(npoint, nDim) = 0; - AD_OutputIndex.resize(npoint, nDim) = 0; + AD_InputIndex.resize(npoint, nDim) = AD::GetPassiveIndex(); + AD_OutputIndex.resize(npoint, nDim) = AD::GetPassiveIndex(); } /*--- Multigrid structures. ---*/ diff --git a/SU2_CFD/include/drivers/CDiscAdjMultizoneDriver.hpp b/SU2_CFD/include/drivers/CDiscAdjMultizoneDriver.hpp index 13da05669bd..f7c8929c805 100644 --- a/SU2_CFD/include/drivers/CDiscAdjMultizoneDriver.hpp +++ b/SU2_CFD/include/drivers/CDiscAdjMultizoneDriver.hpp @@ -96,7 +96,7 @@ class CDiscAdjMultizoneDriver : public CMultizoneDriver { bool eval_transfer = false; /*!< \brief Evaluate the transfer section of the tape. */ su2double ObjFunc; /*!< \brief Value of the objective function. */ - int ObjFunc_Index; /*!< \brief Index of the value of the objective function. */ + AD::Identifier ObjFunc_Index; /*!< \brief Index of the value of the objective function. */ CIteration*** direct_iteration; /*!< \brief Array of pointers to the direct iterations. */ COutput** direct_output; /*!< \brief Array of pointers to the direct outputs. */ diff --git a/SU2_CFD/include/variables/CVariable.hpp b/SU2_CFD/include/variables/CVariable.hpp index 920f5680915..67e298afdd6 100644 --- a/SU2_CFD/include/variables/CVariable.hpp +++ b/SU2_CFD/include/variables/CVariable.hpp @@ -92,8 +92,8 @@ class CVariable { MatrixType Solution_BGS_k; /*!< \brief Old solution container for BGS iterations. */ - su2matrix AD_InputIndex; /*!< \brief Indices of Solution variables in the adjoint vector. */ - su2matrix AD_OutputIndex; /*!< \brief Indices of Solution variables in the adjoint vector after having been updated. */ + su2matrix AD_InputIndex; /*!< \brief Indices of Solution variables in the adjoint vector. */ + su2matrix AD_OutputIndex; /*!< \brief Indices of Solution variables in the adjoint vector after having been updated. */ VectorType SolutionExtra; /*!< \brief Stores adjoint solution for extra solution variables. Currently only streamwise periodic pressure-drop for massflow prescribed flows. */ @@ -118,7 +118,7 @@ class CVariable { assert(false && "A base method of CVariable was used, but it should have been overridden by the derived class."); } - void RegisterContainer(bool input, su2activematrix& variable, su2matrix* ad_index = nullptr) { + void RegisterContainer(bool input, su2activematrix& variable, su2matrix* ad_index = nullptr) { const auto nPoint = variable.rows(); SU2_OMP_FOR_STAT(roundUpDiv(nPoint,omp_get_num_threads())) for (unsigned long iPoint = 0; iPoint < nPoint; ++iPoint) { @@ -133,7 +133,7 @@ class CVariable { END_SU2_OMP_FOR } - void RegisterContainer(bool input, su2activematrix& variable, su2matrix& ad_index) { + void RegisterContainer(bool input, su2activematrix& variable, su2matrix& ad_index) { RegisterContainer(input, variable, &ad_index); } @@ -2180,7 +2180,7 @@ class CVariable { } inline void GetAdjointSolution_time_n(unsigned long iPoint, su2double *adj_sol) const { - int index = 0; + AD::Identifier index = AD::GetPassiveIndex(); for (unsigned long iVar = 0; iVar < Solution_time_n.cols(); iVar++) { AD::SetIndex(index, Solution_time_n(iPoint, iVar)); adj_sol[iVar] = AD::GetDerivative(index); @@ -2188,7 +2188,7 @@ class CVariable { } inline void GetAdjointSolution_time_n1(unsigned long iPoint, su2double *adj_sol) const { - int index = 0; + AD::Identifier index = AD::GetPassiveIndex(); for (unsigned long iVar = 0; iVar < Solution_time_n1.cols(); iVar++) { AD::SetIndex(index, Solution_time_n1(iPoint, iVar)); adj_sol[iVar] = AD::GetDerivative(index); diff --git a/SU2_CFD/src/drivers/CDiscAdjMultizoneDriver.cpp b/SU2_CFD/src/drivers/CDiscAdjMultizoneDriver.cpp index 14d84721fd7..6a8c51fc0fa 100644 --- a/SU2_CFD/src/drivers/CDiscAdjMultizoneDriver.cpp +++ b/SU2_CFD/src/drivers/CDiscAdjMultizoneDriver.cpp @@ -744,7 +744,7 @@ void CDiscAdjMultizoneDriver::SetObjFunction(RECORDING kind_recording) { if (kind_recording == RECORDING::SOLUTION_VARIABLES) { cout << " Objective function : " << ObjFunc; if (driver_config->GetWrt_AD_Statistics()){ - cout << " (" << ObjFunc_Index << ")\n"; + // cout << " (" << ObjFunc_Index << ")\n"; } cout << endl; } diff --git a/SU2_CFD/src/solvers/CSolver.cpp b/SU2_CFD/src/solvers/CSolver.cpp index bfbd4edad58..e2f5927190d 100644 --- a/SU2_CFD/src/solvers/CSolver.cpp +++ b/SU2_CFD/src/solvers/CSolver.cpp @@ -4071,7 +4071,7 @@ void CSolver::SetVertexTractionsAdjoint(CGeometry *geometry, const CConfig *conf unsigned short iMarker, iDim; unsigned long iVertex, iPoint; - int index; + AD::Identifier index; /*--- Loop over all the markers ---*/ for (iMarker = 0; iMarker < config->GetnMarker_All(); iMarker++) { diff --git a/SU2_CFD/src/variables/CVariable.cpp b/SU2_CFD/src/variables/CVariable.cpp index 9f307fd49f2..a7ad2114c75 100644 --- a/SU2_CFD/src/variables/CVariable.cpp +++ b/SU2_CFD/src/variables/CVariable.cpp @@ -69,8 +69,8 @@ CVariable::CVariable(unsigned long npoint, unsigned long ndim, unsigned long nva External.resize(nPoint,nVar) = su2double(0.0); if (!adjoint) { - AD_InputIndex.resize(nPoint,nVar) = -1; - AD_OutputIndex.resize(nPoint,nVar) = -1; + AD_InputIndex.resize(nPoint,nVar) = AD::GetPassiveIndex(); + AD_OutputIndex.resize(nPoint,nVar) = AD::GetPassiveIndex(); } } diff --git a/externals/codi b/externals/codi index c6b039e5c9e..b4c6bcb9d5f 160000 --- a/externals/codi +++ b/externals/codi @@ -1 +1 @@ -Subproject commit c6b039e5c9edb7675f90ffc725f9dd8e66571264 +Subproject commit b4c6bcb9d5f7dda87303f822267251430d41a370 From dc38cc33eef1304e22b9642663d55cbc9ea6b928 Mon Sep 17 00:00:00 2001 From: Ole Burghardt Date: Fri, 16 Aug 2024 00:48:49 +0200 Subject: [PATCH 03/22] Add tag type option to build system. --- meson.build | 2 ++ meson_options.txt | 2 +- 2 files changed, 3 insertions(+), 1 deletion(-) diff --git a/meson.build b/meson.build index 66dde2dfc9d..491b1f8ab5b 100644 --- a/meson.build +++ b/meson.build @@ -95,6 +95,8 @@ if get_option('enable-autodiff') and not omp codi_rev_args += '-DCODI_PRIMAL_REUSE_TAPE' elif get_option('codi-tape') == 'PrimalMultiUse' codi_rev_args += '-DCODI_PRIMAL_MULTIUSE_TAPE' + elif get_option('codi-tape') == 'Tag' + codi_rev_args += '-DCODI_TAG_TAPE' else error('Invalid CoDiPack tape choice @0@'.format(get_option('codi-tape'))) endif diff --git a/meson_options.txt b/meson_options.txt index 08fbac80669..6bbc43a928a 100644 --- a/meson_options.txt +++ b/meson_options.txt @@ -23,7 +23,7 @@ option('enable-coolprop', type : 'boolean', value : false, description: 'enable option('enable-mlpcpp', type : 'boolean', value : false, description: 'enable profiling through gprof') option('enable-gprof', type : 'boolean', value : false, description: 'enable MLPCpp support') option('opdi-backend', type : 'combo', choices : ['auto', 'macro', 'ompt'], value : 'auto', description: 'OpDiLib backend choice') -option('codi-tape', type : 'combo', choices : ['JacobianLinear', 'JacobianReuse', 'JacobianMultiUse', 'PrimalLinear', 'PrimalReuse', 'PrimalMultiUse'], value : 'JacobianLinear', description: 'CoDiPack tape choice') +option('codi-tape', type : 'combo', choices : ['JacobianLinear', 'JacobianReuse', 'JacobianMultiUse', 'PrimalLinear', 'PrimalReuse', 'PrimalMultiUse', 'Tag'], value : 'JacobianLinear', description: 'CoDiPack tape choice') option('opdi-shared-read-opt', type : 'boolean', value : true, description : 'OpDiLib shared reading optimization') option('librom_root', type : 'string', value : '', description: 'libROM base directory') option('enable-librom', type : 'boolean', value : false, description: 'enable LLNL libROM support') From 19a1a70a520d50d4194bf609ceb3da97a5fec1b4 Mon Sep 17 00:00:00 2001 From: Ole Burghardt Date: Mon, 26 Aug 2024 00:07:43 +0200 Subject: [PATCH 04/22] Add option to use RealReverseTag as AD type. --- Common/include/code_config.hpp | 2 ++ externals/codi | 2 +- 2 files changed, 3 insertions(+), 1 deletion(-) diff --git a/Common/include/code_config.hpp b/Common/include/code_config.hpp index 53e18639950..321c1630245 100644 --- a/Common/include/code_config.hpp +++ b/Common/include/code_config.hpp @@ -110,6 +110,8 @@ using su2double = codi::RealReversePrimal; using su2double = codi::RealReversePrimalIndexGen >; #elif defined(CODI_PRIMAL_MULTIUSE_TAPE) using su2double = codi::RealReversePrimalIndex; +#elif defined(CODI_TAG_TAPE) +using su2double = codi::RealReverseTag; #else #error "Please define a CoDiPack tape." #endif diff --git a/externals/codi b/externals/codi index b4c6bcb9d5f..762ba7698e3 160000 --- a/externals/codi +++ b/externals/codi @@ -1 +1 @@ -Subproject commit b4c6bcb9d5f7dda87303f822267251430d41a370 +Subproject commit 762ba7698e3ceaa1b17aa299421ddd418f00b823 From 36c922ef983b47b56d25a08e4be07454fa4b1001 Mon Sep 17 00:00:00 2001 From: Ole Burghardt Date: Mon, 26 Aug 2024 00:28:44 +0200 Subject: [PATCH 05/22] Temporarily adjust/hardcode type size for memory alignment in vectorization.hpp (to cover RealReverseTag). --- Common/include/parallelization/vectorization.hpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/Common/include/parallelization/vectorization.hpp b/Common/include/parallelization/vectorization.hpp index 9f6db32a912..5a5c17ac23b 100644 --- a/Common/include/parallelization/vectorization.hpp +++ b/Common/include/parallelization/vectorization.hpp @@ -90,11 +90,11 @@ class Array : public CVecExpr, Scalar_t> { public: using Scalar = Scalar_t; enum : size_t { Size = N }; - enum : size_t { Align = Size * sizeof(Scalar) }; + enum : size_t { Align = Size * 32 }; static constexpr bool StoreAsRef = true; private: - alignas(Size * sizeof(Scalar)) Scalar x_[N]; + alignas(Size * 32) Scalar x_[N]; public: #define ARRAY_BOILERPLATE \ From 84e97f42e78eb29cf0d3c458eacde96fe854f315 Mon Sep 17 00:00:00 2001 From: Josh Kelly Date: Thu, 29 Aug 2024 16:37:21 +0200 Subject: [PATCH 06/22] added turbo obj funcs --- Common/include/CConfig.hpp | 47 ++++++++++++++++++- Common/include/option_structure.hpp | 6 +++ Common/src/CConfig.cpp | 13 +++++ SU2_CFD/include/output/CTurboOutput.hpp | 10 +++- .../include/solvers/CFVMFlowSolverBase.inl | 6 +++ SU2_CFD/src/iteration/CFluidIteration.cpp | 20 ++++++-- SU2_CFD/src/output/CFlowCompOutput.cpp | 2 +- SU2_CFD/src/output/CTurboOutput.cpp | 19 ++++++-- 8 files changed, 113 insertions(+), 10 deletions(-) diff --git a/Common/include/CConfig.hpp b/Common/include/CConfig.hpp index df9a0dba3b7..7bfb9cb3173 100644 --- a/Common/include/CConfig.hpp +++ b/Common/include/CConfig.hpp @@ -441,6 +441,11 @@ class CConfig { TURBO_PERF_KIND *Kind_TurboPerf; /*!< \brief Kind of turbomachynery architecture.*/ TURBOMACHINERY_TYPE *Kind_TurboMachinery; + /* Turbomachinery objective functions */ + su2double *EntropyGeneration; + su2double *TotalPressureLoss; + su2double *KineticEnergyLoss; + /* Gradient smoothing options */ su2double SmoothingEps1; /*!< \brief Parameter for the identity part in gradient smoothing. */ su2double SmoothingEps2; /*!< \brief Parameter for the Laplace part in gradient smoothing. */ @@ -7930,7 +7935,28 @@ class CConfig { void SetSurface_Species_Variance(unsigned short val_marker, su2double val_surface_species_variance) { Surface_Species_Variance[val_marker] = val_surface_species_variance; } /*! - * \brief Get the back pressure (static) at an outlet boundary. + * \brief Set entropy generation for a turbomachinery zone + * \param[in] val_iZone - zone index + * \param[in] val_entropy_generation - value of entropy generation + */ + void SetEntropyGeneration(unsigned short val_iZone, su2double val_entropy_generation) { EntropyGeneration[val_iZone] = val_entropy_generation; } + + /*! + * \brief Get total pressure loss for a turbomachinery zone + * \param[in] val_iZone - zone index + * \param[in] val_total_pressure_loss - value of total pressure loss + */ + void SetTotalPressureLoss(unsigned short val_iZone, su2double val_total_pressure_loss) { TotalPressureLoss[val_iZone] = val_total_pressure_loss; } + + /*! + * \brief Get kinetic energy loss for a turbomachinery zone + * \param[in] val_iZone - zone index + * \param[in] val_kinetic_energy_loss - value of kinetic energy loss + */ + void SetKineticEnergyLoss(unsigned short val_iZone, su2double val_kinetic_energy_loss) { KineticEnergyLoss[val_iZone] = val_kinetic_energy_loss; } + + /*! + * \brief Get the back pressure (static) at an outlet boundary.ß * \param[in] val_index - Index corresponding to the outlet boundary. * \return The outlet pressure. */ @@ -8209,6 +8235,25 @@ class CConfig { */ su2double GetSurface_Species_Variance(unsigned short val_marker) const { return Surface_Species_Variance[val_marker]; } + /*! + * \brief Get entropy generation for a turbomachine at a boundary + * \param[in] val_iZone - zone index + */ + su2double GetEntropyGeneration(unsigned short val_iZone) const { return EntropyGeneration[val_iZone]; } + + /*! + * \brief Get total pressure loss for a turbomachinery zone + * \param[in] val_iZone - zone index + */ + su2double GetTotalPressureLoss(unsigned short val_iZone) const { return TotalPressureLoss[val_iZone]; } + + /*! + * \brief Get kinetic energy loss for a turbomachinery zone + * \param[in] val_iZone - zone index + */ + su2double GetKineticEnergyLoss(unsigned short val_iZone) const { return KineticEnergyLoss[val_iZone]; } + + /*! * \brief Get the back pressure (static) at an outlet boundary. * \param[in] val_index - Index corresponding to the outlet boundary. diff --git a/Common/include/option_structure.hpp b/Common/include/option_structure.hpp index 64c781a34f8..d81c49448aa 100644 --- a/Common/include/option_structure.hpp +++ b/Common/include/option_structure.hpp @@ -1993,6 +1993,9 @@ enum ENUM_OBJECTIVE { TOPOL_DISCRETENESS = 63, /*!< \brief Measure of the discreteness of the current topology. */ TOPOL_COMPLIANCE = 64, /*!< \brief Measure of the discreteness of the current topology. */ STRESS_PENALTY = 65, /*!< \brief Penalty function of VM stresses above a maximum value. */ + ENTROPY_GENERATION = 80, /*!< \brief Entropy generation turbomachinery objective function. */ + TOTAL_PRESSURE_LOSS = 81, /*!< \brief Total pressure loss turbomachinery objective function. */ + KINETIC_ENERGY_LOSS = 82 /*!< \breif Kinetic energy loss coefficient turbomachinery objective function. */ }; static const MapType Objective_Map = { MakePair("DRAG", DRAG_COEFFICIENT) @@ -2035,6 +2038,9 @@ static const MapType Objective_Map = { MakePair("TOPOL_DISCRETENESS", TOPOL_DISCRETENESS) MakePair("TOPOL_COMPLIANCE", TOPOL_COMPLIANCE) MakePair("STRESS_PENALTY", STRESS_PENALTY) + MakePair("ENTROPY_GENERATION", ENTROPY_GENERATION) + MakePair("TOTAL_PRESSURE_LOSS", TOTAL_PRESSURE_LOSS) + MakePair("KINETIC_ENERGY_LOSS", KINETIC_ENERGY_LOSS) }; /*! diff --git a/Common/src/CConfig.cpp b/Common/src/CConfig.cpp index 7aa9d7af84a..285dbd8465e 100644 --- a/Common/src/CConfig.cpp +++ b/Common/src/CConfig.cpp @@ -1043,6 +1043,11 @@ void CConfig::SetPointersNull() { nBlades = nullptr; FreeStreamTurboNormal = nullptr; + /*--- Turbomachinery Objective Functions ---*/ + EntropyGeneration = nullptr; + TotalPressureLoss = nullptr; + KineticEnergyLoss = nullptr; + top_optim_kernels = nullptr; top_optim_kernel_params = nullptr; top_optim_filter_radius = nullptr; @@ -5993,6 +5998,11 @@ void CConfig::SetMarkers(SU2_COMPONENT val_software) { Marker_CfgFile_ZoneInterface[iMarker_CfgFile] = YES; } + /*--- Allocate memory for turbomachinery objective functions ---*/ + EntropyGeneration = new su2double[nZone] (); + TotalPressureLoss = new su2double[nZone] (); + KineticEnergyLoss = new su2double[nZone] (); + /*--- Identification of Turbomachinery markers and flag them---*/ for (iMarker_CfgFile = 0; iMarker_CfgFile < nMarker_CfgFile; iMarker_CfgFile++) { @@ -8426,6 +8436,9 @@ string CConfig::GetObjFunc_Extension(string val_filename) const { case TOPOL_DISCRETENESS: AdjExt = "_topdisc"; break; case TOPOL_COMPLIANCE: AdjExt = "_topcomp"; break; case STRESS_PENALTY: AdjExt = "_stress"; break; + case ENTROPY_GENERATION: AdjExt = "_entg"; break; + case TOTAL_PRESSURE_LOSS: AdjExt = "_tot_press_loss"; break; + case KINETIC_ENERGY_LOSS: AdjExt = "_kin_en_loss"; break; } } else{ diff --git a/SU2_CFD/include/output/CTurboOutput.hpp b/SU2_CFD/include/output/CTurboOutput.hpp index cc7e2356f07..5a1741530cd 100644 --- a/SU2_CFD/include/output/CTurboOutput.hpp +++ b/SU2_CFD/include/output/CTurboOutput.hpp @@ -87,7 +87,7 @@ class CTurbomachineryCombinedPrimitiveStates { class CTurbomachineryState { private: su2double Density, Pressure, Entropy, Enthalpy, Temperature, TotalTemperature, TotalPressure, TotalEnthalpy; - su2double AbsFlowAngle, FlowAngle, MassFlow, Rothalpy, TotalRelPressure; + su2double AbsFlowAngle, FlowAngle, MassFlow, Rothalpy, TangVelocity, TotalRelPressure; vector Velocity, RelVelocity, Mach, RelMach; su2double Area, Radius; @@ -124,6 +124,8 @@ class CTurbomachineryState { const su2double& GetRothalpy() const { return Rothalpy; } + const su2double& GetTangVelocity() const { return TangVelocity; } + const vector& GetVelocity() const { return Velocity; } const vector& GetMach() const { return Mach; } @@ -207,7 +209,7 @@ class CPropellorBladePerformance : public CTurbomachineryBladePerformance { */ class CTurbomachineryStagePerformance { protected: - su2double TotalStaticEfficiency, TotalTotalEfficiency, NormEntropyGen, TotalStaticPressureRatio, TotalTotalPressureRatio, EulerianWork; + su2double TotalStaticEfficiency, TotalTotalEfficiency, NormEntropyGen, TotalStaticPressureRatio, TotalTotalPressureRatio, EulerianWork, TotalPressureLoss, KineticEnergyLoss; CFluidModel& fluidModel; public: @@ -232,6 +234,10 @@ class CTurbomachineryStagePerformance { su2double GetTotalStaticPressureRatio() const { return TotalStaticPressureRatio; } su2double GetTotalTotalPressureRatio() const { return TotalTotalPressureRatio; } + + su2double GetTotalPressureLoss() const { return TotalPressureLoss; } + + su2double GetKineticEnergyLoss() const { return KineticEnergyLoss; } }; /*! diff --git a/SU2_CFD/include/solvers/CFVMFlowSolverBase.inl b/SU2_CFD/include/solvers/CFVMFlowSolverBase.inl index f60db438e18..3ca3e58834f 100644 --- a/SU2_CFD/include/solvers/CFVMFlowSolverBase.inl +++ b/SU2_CFD/include/solvers/CFVMFlowSolverBase.inl @@ -2908,6 +2908,12 @@ su2double CFVMFlowSolverBase::EvaluateCommonObjFunc(const CConfig& config) case SURFACE_SPECIES_VARIANCE: objFun += weight * config.GetSurface_Species_Variance(0); break; + case ENTROPY_GENERATION: + objFun += weight * config.GetEntropyGeneration(config.GetnMarker_Turbomachinery()); + case TOTAL_PRESSURE_LOSS: + objFun += weight * config.GetTotalPressureLoss(config.GetnMarker_Turbomachinery()); + case KINETIC_ENERGY_LOSS: + objFun += weight * config.GetKineticEnergyLoss(config.GetnMarker_Turbomachinery()); case CUSTOM_OBJFUNC: objFun += weight * Total_Custom_ObjFunc; break; diff --git a/SU2_CFD/src/iteration/CFluidIteration.cpp b/SU2_CFD/src/iteration/CFluidIteration.cpp index bd60c38a961..180e40c9b44 100644 --- a/SU2_CFD/src/iteration/CFluidIteration.cpp +++ b/SU2_CFD/src/iteration/CFluidIteration.cpp @@ -325,16 +325,15 @@ void CFluidIteration::TurboMonitor(CGeometry**** geometry_container, CConfig** c void CFluidIteration::ComputeTurboPerformance(CSolver***** solver, CGeometry**** geometry_container, CConfig** config_container, unsigned long ExtIter) { unsigned short nDim = geometry_container[ZONE_0][INST_0][MESH_0]->GetnDim(); unsigned short nBladesRow = config_container[ZONE_0]->GetnMarker_Turbomachinery(); - unsigned short iBlade=0, iSpan; vector TurboPrimitiveIn, TurboPrimitiveOut; std::vector> bladesPrimitives; if (rank == MASTER_NODE) { - for (iBlade = 0; iBlade < nBladesRow; iBlade++){ + for (auto iBlade = 0u; iBlade < nBladesRow; iBlade++){ /* Blade Primitive initialized per blade */ std::vector bladePrimitives; auto nSpan = config_container[iBlade]->GetnSpanWiseSections(); - for (iSpan = 0; iSpan < nSpan + 1; iSpan++) { + for (auto iSpan = 0u; iSpan < nSpan + 1; iSpan++) { TurboPrimitiveIn= solver[iBlade][INST_0][MESH_0][FLOW_SOL]->GetTurboPrimitive(iBlade, iSpan, true); TurboPrimitiveOut= solver[iBlade][INST_0][MESH_0][FLOW_SOL]->GetTurboPrimitive(iBlade, iSpan, false); auto spanInletPrimitive = CTurbomachineryPrimitiveState(TurboPrimitiveIn, nDim, geometry_container[iBlade][INST_0][MESH_0]->GetTangGridVelIn(iBlade, iSpan)); @@ -352,6 +351,21 @@ void CFluidIteration::ComputeTurboPerformance(CSolver***** solver, CGeometry**** auto OutState = TurbomachineryPerformance->GetBladesPerformances().at(nZone-1).at(nSpan)->GetOutletState(); TurbomachineryStagePerformance->ComputePerformanceStage(InState, OutState, config_container[nZone-1]); + + /*--- Set turbomachinery objective function value in each zone ---*/ + for (auto iBlade = 0u; iBlade < nBladesRow; iBlade++) { + // Should we set in ZONE_0 or nZone-1? + auto iBladePerf = TurbomachineryPerformance->GetBladesPerformances().at(iBlade).at(nSpan); + InState = iBladePerf->GetInletState(); + OutState = iBladePerf->GetOutletState(); + config_container[nZone-1]->SetEntropyGeneration(iBlade, (OutState.GetEntropy() - InState.GetEntropy())/InState.GetEntropy() * 100); + config_container[nZone-1]->SetTotalPressureLoss(iBlade, iBladePerf->GetTotalPressureLoss()); + config_container[nZone-1]->SetKineticEnergyLoss(iBlade, iBladePerf->GetKineticEnergyLoss()); + } + /*--- Set global turbomachinery objective function ---*/ + config_container[nZone-1]->SetEntropyGeneration(nBladesRow, TurbomachineryStagePerformance->GetNormEntropyGen()*100); + config_container[nZone-1]->SetTotalPressureLoss(nBladesRow, TurbomachineryStagePerformance->GetTotalPressureLoss()); + config_container[nZone-1]->SetKineticEnergyLoss(nBladesRow, TurbomachineryStagePerformance->GetKineticEnergyLoss()); } } diff --git a/SU2_CFD/src/output/CFlowCompOutput.cpp b/SU2_CFD/src/output/CFlowCompOutput.cpp index 1666191e07f..5bb5577fe76 100644 --- a/SU2_CFD/src/output/CFlowCompOutput.cpp +++ b/SU2_CFD/src/output/CFlowCompOutput.cpp @@ -577,7 +577,7 @@ void CFlowCompOutput::SetTurboMultiZonePerformance_Output(std::shared_ptr TurboStagePerf, std::shared_ptr TurboPerf, CConfig *config) { auto BladePerformance = TurboPerf->GetBladesPerformances(); - for (unsigned short iZone = 0; iZone <= config->GetnZone()-1; iZone++) { + for (auto iZone = 0u; iZone < config->GetnZone(); iZone++) { auto nSpan = config->GetnSpan_iZones(iZone); const auto& BladePerf = BladePerformance.at(iZone).at(nSpan); diff --git a/SU2_CFD/src/output/CTurboOutput.cpp b/SU2_CFD/src/output/CTurboOutput.cpp index b0e61f10e25..6551d0062c1 100644 --- a/SU2_CFD/src/output/CTurboOutput.cpp +++ b/SU2_CFD/src/output/CTurboOutput.cpp @@ -60,7 +60,7 @@ void CTurbomachineryState::ComputeState(CFluidModel& fluidModel, const CTurbomac Pressure = primitiveState.GetPressure(); std::vector velocity = primitiveState.GetVelocity(); Velocity.assign(velocity.begin(), velocity.end()); - su2double tangVel = primitiveState.GetTangVelocity(); + TangVelocity = primitiveState.GetTangVelocity(); /*--- Compute static TD quantities ---*/ fluidModel.SetTDState_Prho(Pressure, Density); @@ -81,9 +81,9 @@ void CTurbomachineryState::ComputeState(CFluidModel& fluidModel, const CTurbomac std::for_each(Mach.begin(), Mach.end(), [&](su2double& el) { el /= soundSpeed; }); /*--- Compute relative kinematic quantities ---*/ - su2double tangVel2 = tangVel * tangVel; + su2double tangVel2 = TangVelocity * TangVelocity; RelVelocity.assign(Velocity.begin(), Velocity.end()); - RelVelocity[1] -= tangVel; + RelVelocity[1] -= TangVelocity; su2double relVel2 = GetRelVelocityValue(); FlowAngle = atan(RelVelocity[1] / RelVelocity[0]); RelMach.assign(RelVelocity.begin(), RelVelocity.end()); @@ -245,6 +245,9 @@ void CTurbomachineryStagePerformance::ComputeTurbineStagePerformance(const CTurb su2double enthalpyOutIs = fluidModel.GetStaticEnergy() + OutState.GetPressure() / fluidModel.GetDensity(); su2double totEnthalpyOutIs = enthalpyOutIs + 0.5 * OutState.GetVelocityValue() * OutState.GetVelocityValue(); + su2double tangVel = OutState.GetTangVelocity(); + su2double relVelOutIs2 = 2 * (OutState.GetRothalpy() - enthalpyOutIs) + tangVel * tangVel; + /*--- Compute turbine stage performance ---*/ NormEntropyGen = (OutState.GetEntropy() - InState.GetEntropy()) / InState.GetEntropy(); EulerianWork = InState.GetTotalEnthalpy() - OutState.GetTotalEnthalpy(); @@ -252,6 +255,10 @@ void CTurbomachineryStagePerformance::ComputeTurbineStagePerformance(const CTurb TotalTotalEfficiency = EulerianWork / (InState.GetTotalEnthalpy() - totEnthalpyOutIs); TotalStaticPressureRatio = InState.GetTotalPressure() / OutState.GetPressure(); TotalTotalPressureRatio = InState.GetTotalPressure() / OutState.GetTotalPressure(); + + TotalPressureLoss = (InState.GetTotalRelPressure() - OutState.GetTotalRelPressure()) / + (OutState.GetTotalRelPressure() - OutState.GetPressure()); + KineticEnergyLoss = 2 * (OutState.GetEnthalpy() - enthalpyOutIs) / relVelOutIs2; } void CTurbomachineryStagePerformance::ComputeCompressorStagePerformance(const CTurbomachineryState& InState, @@ -260,6 +267,8 @@ void CTurbomachineryStagePerformance::ComputeCompressorStagePerformance(const CT fluidModel.SetTDState_Ps(OutState.GetPressure(), InState.GetEntropy()); su2double enthalpyOutIs = fluidModel.GetStaticEnergy() + OutState.GetPressure() / fluidModel.GetDensity(); su2double totEnthalpyOutIs = enthalpyOutIs + 0.5 * OutState.GetVelocityValue() * OutState.GetVelocityValue(); + su2double tangVel = OutState.GetTangVelocity(); + su2double relVelOutIs2 = 2 * (OutState.GetRothalpy() - enthalpyOutIs) + tangVel * tangVel; /*--- Compute compressor stage performance ---*/ NormEntropyGen = (OutState.GetEntropy() - InState.GetEntropy()) / InState.GetEntropy(); @@ -268,4 +277,8 @@ void CTurbomachineryStagePerformance::ComputeCompressorStagePerformance(const CT TotalTotalEfficiency = (totEnthalpyOutIs - InState.GetTotalEnthalpy()) / EulerianWork; TotalStaticPressureRatio = OutState.GetPressure() / InState.GetTotalPressure(); TotalTotalPressureRatio = OutState.GetTotalPressure() / InState.GetTotalPressure(); + + TotalPressureLoss = (InState.GetTotalRelPressure() - OutState.GetTotalRelPressure()) / + (InState.GetTotalRelPressure() - InState.GetPressure()); + KineticEnergyLoss = 2 * (OutState.GetEnthalpy() - enthalpyOutIs) / relVelOutIs2; } \ No newline at end of file From bb0831df61bfa34eb2aa852e06de560df58e3fff Mon Sep 17 00:00:00 2001 From: Josh Kelly Date: Thu, 29 Aug 2024 18:25:50 +0200 Subject: [PATCH 07/22] thems the breaks --- SU2_CFD/include/solvers/CFVMFlowSolverBase.inl | 3 +++ 1 file changed, 3 insertions(+) diff --git a/SU2_CFD/include/solvers/CFVMFlowSolverBase.inl b/SU2_CFD/include/solvers/CFVMFlowSolverBase.inl index 3ca3e58834f..4bd470933f5 100644 --- a/SU2_CFD/include/solvers/CFVMFlowSolverBase.inl +++ b/SU2_CFD/include/solvers/CFVMFlowSolverBase.inl @@ -2910,10 +2910,13 @@ su2double CFVMFlowSolverBase::EvaluateCommonObjFunc(const CConfig& config) break; case ENTROPY_GENERATION: objFun += weight * config.GetEntropyGeneration(config.GetnMarker_Turbomachinery()); + break; case TOTAL_PRESSURE_LOSS: objFun += weight * config.GetTotalPressureLoss(config.GetnMarker_Turbomachinery()); + break; case KINETIC_ENERGY_LOSS: objFun += weight * config.GetKineticEnergyLoss(config.GetnMarker_Turbomachinery()); + break; case CUSTOM_OBJFUNC: objFun += weight * Total_Custom_ObjFunc; break; From 905f1e6f1971700073ecd9f0bb31728fc4fa8c15 Mon Sep 17 00:00:00 2001 From: Josh Kelly Date: Fri, 30 Aug 2024 11:29:30 +0200 Subject: [PATCH 08/22] Resource releasing --- Common/src/CConfig.cpp | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/Common/src/CConfig.cpp b/Common/src/CConfig.cpp index c40a80d40f9..d0e50ebc3d1 100644 --- a/Common/src/CConfig.cpp +++ b/Common/src/CConfig.cpp @@ -8266,6 +8266,10 @@ CConfig::~CConfig() { delete [] nBlades; delete [] FreeStreamTurboNormal; + + delete [] EntropyGeneration; + delete [] TotalPressureLoss; + delete [] KineticEnergyLoss; } string CConfig::GetFilename(string filename, const string& ext, int timeIter) const { From c002911246cce5d2debcec36a77d0cea3e73a7a9 Mon Sep 17 00:00:00 2001 From: Josh Kelly Date: Tue, 3 Sep 2024 15:26:40 +0200 Subject: [PATCH 09/22] Revert "Merge remote-tracking branch 'upstream/feature_tag_debug_tape' into feature_mz_adjoint_for_turbo" This reverts commit 9f1b07ca8fb85f7fb3e1d723087a81abe5c4e792, reversing changes made to 4327942aad1126c5b7c3d2f19153b5cf293d9cfe. --- Common/include/basic_types/ad_structure.hpp | 26 +++++-------------- Common/include/code_config.hpp | 2 -- Common/include/geometry/dual_grid/CPoint.hpp | 5 ++-- .../include/parallelization/vectorization.hpp | 4 +-- Common/src/geometry/dual_grid/CPoint.cpp | 4 +-- .../drivers/CDiscAdjMultizoneDriver.hpp | 2 +- SU2_CFD/include/variables/CVariable.hpp | 12 ++++----- .../src/drivers/CDiscAdjMultizoneDriver.cpp | 2 +- SU2_CFD/src/solvers/CSolver.cpp | 2 +- SU2_CFD/src/variables/CVariable.cpp | 4 +-- externals/codi | 2 +- meson.build | 2 -- meson_options.txt | 2 +- 13 files changed, 27 insertions(+), 42 deletions(-) diff --git a/Common/include/basic_types/ad_structure.hpp b/Common/include/basic_types/ad_structure.hpp index ecd2c93b1c1..c44374f7eb7 100644 --- a/Common/include/basic_types/ad_structure.hpp +++ b/Common/include/basic_types/ad_structure.hpp @@ -38,9 +38,6 @@ */ namespace AD { #ifndef CODI_REVERSE_TYPE - -using Identifier = int; - /*! * \brief Start the recording of the operations and involved variables. * If called, the computational graph of all operations occuring after the call will be stored, @@ -104,20 +101,14 @@ inline void EndUseAdjoints() {} * \param[in] index - Position in the adjoint vector. * \param[in] val - adjoint value to be set. */ -inline void SetDerivative(Identifier index, const double val) {} +inline void SetDerivative(int index, const double val) {} /*! * \brief Extracts the adjoint value at index * \param[in] index - position in the adjoint vector where the derivative will be extracted. * \return Derivative value. */ -inline double GetDerivative(Identifier index) { return 0.0; } - -/*! - * \brief Returns the identifier that represents an inactive variable. - * \return Passive index. - */ -inline Identifier GetPassiveIndex() { return 0; } +inline double GetDerivative(int index) { return 0.0; } /*! * \brief Clears the currently stored adjoints but keeps the computational graph. @@ -268,7 +259,7 @@ inline void SetExtFuncOut(T&& data, const int size_x, const int size_y) {} * \param[in] data - variable whose gradient information will be extracted. * \param[in] index - where obtained gradient information will be stored. */ -inline void SetIndex(Identifier& index, const su2double& data) {} +inline void SetIndex(int& index, const su2double& data) {} /*! * \brief Pushes back the current tape position to the tape position's vector. @@ -313,7 +304,6 @@ inline void EndNoSharedReading() {} using CheckpointHandler = codi::ExternalFunctionUserData; using Tape = su2double::Tape; -using Identifier = su2double::Identifier; #ifdef HAVE_OPDI using ExtFuncHelper = codi::OpenMPExternalFunctionHelper; @@ -480,14 +470,14 @@ FORCEINLINE void BeginUseAdjoints() { AD::getTape().beginUseAdjointVector(); } FORCEINLINE void EndUseAdjoints() { AD::getTape().endUseAdjointVector(); } -FORCEINLINE void SetIndex(Identifier& index, const su2double& data) { index = data.getIdentifier(); } +FORCEINLINE void SetIndex(int& index, const su2double& data) { index = data.getIdentifier(); } // 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. // It should be safeguarded by calls to AD::BeginUseAdjoints() and AD::EndUseAdjoints(). -FORCEINLINE void SetDerivative(Identifier index, const double val) { - if (!AD::getTape().isIdentifierActive(index)) // Allow multiple threads to "set the derivative" of passive variables without causing data races. +FORCEINLINE void SetDerivative(int index, const double val) { + if (index == 0) // Allow multiple threads to "set the derivative" of passive variables without causing data races. return; AD::getTape().setGradient(index, val, codi::AdjointsManagement::Manual); @@ -498,12 +488,10 @@ FORCEINLINE void SetDerivative(Identifier index, const double val) { // Otherwise, please ensure sufficient adjoint vector size by a call to AD::ResizeAdjoints(). // This method does not perform locking either. // It should be safeguarded by calls to AD::BeginUseAdjoints() and AD::EndUseAdjoints(). -FORCEINLINE double GetDerivative(Identifier index) { +FORCEINLINE double GetDerivative(int index) { return AD::getTape().getGradient(index, codi::AdjointsManagement::Manual); } -FORCEINLINE Identifier GetPassiveIndex() { return AD::getTape().getPassiveIndex(); } - FORCEINLINE bool IsIdentifierActive(su2double const& value) { return getTape().isIdentifierActive(value.getIdentifier()); } diff --git a/Common/include/code_config.hpp b/Common/include/code_config.hpp index 321c1630245..53e18639950 100644 --- a/Common/include/code_config.hpp +++ b/Common/include/code_config.hpp @@ -110,8 +110,6 @@ using su2double = codi::RealReversePrimal; using su2double = codi::RealReversePrimalIndexGen >; #elif defined(CODI_PRIMAL_MULTIUSE_TAPE) using su2double = codi::RealReversePrimalIndex; -#elif defined(CODI_TAG_TAPE) -using su2double = codi::RealReverseTag; #else #error "Please define a CoDiPack tape." #endif diff --git a/Common/include/geometry/dual_grid/CPoint.hpp b/Common/include/geometry/dual_grid/CPoint.hpp index acae3cd172d..c297b1b5df9 100644 --- a/Common/include/geometry/dual_grid/CPoint.hpp +++ b/Common/include/geometry/dual_grid/CPoint.hpp @@ -113,8 +113,9 @@ class CPoint { su2activevector MaxLength; /*!< \brief The maximum cell-center to cell-center length. */ su2activevector RoughnessHeight; /*!< \brief Roughness of the nearest wall. */ - su2matrix AD_InputIndex; /*!< \brief Indices of Coord variables in the adjoint vector. */ - su2matrix AD_OutputIndex; /*!< \brief Indices of Coord variables in the adjoint vector after having been updated. */ + su2matrix AD_InputIndex; /*!< \brief Indices of Coord variables in the adjoint vector. */ + su2matrix + AD_OutputIndex; /*!< \brief Indices of Coord variables in the adjoint vector after having been updated. */ /*! * \brief Allocate fields required by the minimal constructor. diff --git a/Common/include/parallelization/vectorization.hpp b/Common/include/parallelization/vectorization.hpp index 5a5c17ac23b..9f6db32a912 100644 --- a/Common/include/parallelization/vectorization.hpp +++ b/Common/include/parallelization/vectorization.hpp @@ -90,11 +90,11 @@ class Array : public CVecExpr, Scalar_t> { public: using Scalar = Scalar_t; enum : size_t { Size = N }; - enum : size_t { Align = Size * 32 }; + enum : size_t { Align = Size * sizeof(Scalar) }; static constexpr bool StoreAsRef = true; private: - alignas(Size * 32) Scalar x_[N]; + alignas(Size * sizeof(Scalar)) Scalar x_[N]; public: #define ARRAY_BOILERPLATE \ diff --git a/Common/src/geometry/dual_grid/CPoint.cpp b/Common/src/geometry/dual_grid/CPoint.cpp index 9753d71f645..f11b7dd8945 100644 --- a/Common/src/geometry/dual_grid/CPoint.cpp +++ b/Common/src/geometry/dual_grid/CPoint.cpp @@ -67,8 +67,8 @@ void CPoint::FullAllocation(unsigned short imesh, const CConfig* config) { } if (config->GetDiscrete_Adjoint()) { - AD_InputIndex.resize(npoint, nDim) = AD::GetPassiveIndex(); - AD_OutputIndex.resize(npoint, nDim) = AD::GetPassiveIndex(); + AD_InputIndex.resize(npoint, nDim) = 0; + AD_OutputIndex.resize(npoint, nDim) = 0; } /*--- Multigrid structures. ---*/ diff --git a/SU2_CFD/include/drivers/CDiscAdjMultizoneDriver.hpp b/SU2_CFD/include/drivers/CDiscAdjMultizoneDriver.hpp index f7c8929c805..13da05669bd 100644 --- a/SU2_CFD/include/drivers/CDiscAdjMultizoneDriver.hpp +++ b/SU2_CFD/include/drivers/CDiscAdjMultizoneDriver.hpp @@ -96,7 +96,7 @@ class CDiscAdjMultizoneDriver : public CMultizoneDriver { bool eval_transfer = false; /*!< \brief Evaluate the transfer section of the tape. */ su2double ObjFunc; /*!< \brief Value of the objective function. */ - AD::Identifier ObjFunc_Index; /*!< \brief Index of the value of the objective function. */ + int ObjFunc_Index; /*!< \brief Index of the value of the objective function. */ CIteration*** direct_iteration; /*!< \brief Array of pointers to the direct iterations. */ COutput** direct_output; /*!< \brief Array of pointers to the direct outputs. */ diff --git a/SU2_CFD/include/variables/CVariable.hpp b/SU2_CFD/include/variables/CVariable.hpp index 67e298afdd6..920f5680915 100644 --- a/SU2_CFD/include/variables/CVariable.hpp +++ b/SU2_CFD/include/variables/CVariable.hpp @@ -92,8 +92,8 @@ class CVariable { MatrixType Solution_BGS_k; /*!< \brief Old solution container for BGS iterations. */ - su2matrix AD_InputIndex; /*!< \brief Indices of Solution variables in the adjoint vector. */ - su2matrix AD_OutputIndex; /*!< \brief Indices of Solution variables in the adjoint vector after having been updated. */ + su2matrix AD_InputIndex; /*!< \brief Indices of Solution variables in the adjoint vector. */ + su2matrix AD_OutputIndex; /*!< \brief Indices of Solution variables in the adjoint vector after having been updated. */ VectorType SolutionExtra; /*!< \brief Stores adjoint solution for extra solution variables. Currently only streamwise periodic pressure-drop for massflow prescribed flows. */ @@ -118,7 +118,7 @@ class CVariable { assert(false && "A base method of CVariable was used, but it should have been overridden by the derived class."); } - void RegisterContainer(bool input, su2activematrix& variable, su2matrix* ad_index = nullptr) { + void RegisterContainer(bool input, su2activematrix& variable, su2matrix* ad_index = nullptr) { const auto nPoint = variable.rows(); SU2_OMP_FOR_STAT(roundUpDiv(nPoint,omp_get_num_threads())) for (unsigned long iPoint = 0; iPoint < nPoint; ++iPoint) { @@ -133,7 +133,7 @@ class CVariable { END_SU2_OMP_FOR } - void RegisterContainer(bool input, su2activematrix& variable, su2matrix& ad_index) { + void RegisterContainer(bool input, su2activematrix& variable, su2matrix& ad_index) { RegisterContainer(input, variable, &ad_index); } @@ -2180,7 +2180,7 @@ class CVariable { } inline void GetAdjointSolution_time_n(unsigned long iPoint, su2double *adj_sol) const { - AD::Identifier index = AD::GetPassiveIndex(); + int index = 0; for (unsigned long iVar = 0; iVar < Solution_time_n.cols(); iVar++) { AD::SetIndex(index, Solution_time_n(iPoint, iVar)); adj_sol[iVar] = AD::GetDerivative(index); @@ -2188,7 +2188,7 @@ class CVariable { } inline void GetAdjointSolution_time_n1(unsigned long iPoint, su2double *adj_sol) const { - AD::Identifier index = AD::GetPassiveIndex(); + int index = 0; for (unsigned long iVar = 0; iVar < Solution_time_n1.cols(); iVar++) { AD::SetIndex(index, Solution_time_n1(iPoint, iVar)); adj_sol[iVar] = AD::GetDerivative(index); diff --git a/SU2_CFD/src/drivers/CDiscAdjMultizoneDriver.cpp b/SU2_CFD/src/drivers/CDiscAdjMultizoneDriver.cpp index 6a8c51fc0fa..14d84721fd7 100644 --- a/SU2_CFD/src/drivers/CDiscAdjMultizoneDriver.cpp +++ b/SU2_CFD/src/drivers/CDiscAdjMultizoneDriver.cpp @@ -744,7 +744,7 @@ void CDiscAdjMultizoneDriver::SetObjFunction(RECORDING kind_recording) { if (kind_recording == RECORDING::SOLUTION_VARIABLES) { cout << " Objective function : " << ObjFunc; if (driver_config->GetWrt_AD_Statistics()){ - // cout << " (" << ObjFunc_Index << ")\n"; + cout << " (" << ObjFunc_Index << ")\n"; } cout << endl; } diff --git a/SU2_CFD/src/solvers/CSolver.cpp b/SU2_CFD/src/solvers/CSolver.cpp index e2f5927190d..bfbd4edad58 100644 --- a/SU2_CFD/src/solvers/CSolver.cpp +++ b/SU2_CFD/src/solvers/CSolver.cpp @@ -4071,7 +4071,7 @@ void CSolver::SetVertexTractionsAdjoint(CGeometry *geometry, const CConfig *conf unsigned short iMarker, iDim; unsigned long iVertex, iPoint; - AD::Identifier index; + int index; /*--- Loop over all the markers ---*/ for (iMarker = 0; iMarker < config->GetnMarker_All(); iMarker++) { diff --git a/SU2_CFD/src/variables/CVariable.cpp b/SU2_CFD/src/variables/CVariable.cpp index a7ad2114c75..9f307fd49f2 100644 --- a/SU2_CFD/src/variables/CVariable.cpp +++ b/SU2_CFD/src/variables/CVariable.cpp @@ -69,8 +69,8 @@ CVariable::CVariable(unsigned long npoint, unsigned long ndim, unsigned long nva External.resize(nPoint,nVar) = su2double(0.0); if (!adjoint) { - AD_InputIndex.resize(nPoint,nVar) = AD::GetPassiveIndex(); - AD_OutputIndex.resize(nPoint,nVar) = AD::GetPassiveIndex(); + AD_InputIndex.resize(nPoint,nVar) = -1; + AD_OutputIndex.resize(nPoint,nVar) = -1; } } diff --git a/externals/codi b/externals/codi index 762ba7698e3..c6b039e5c9e 160000 --- a/externals/codi +++ b/externals/codi @@ -1 +1 @@ -Subproject commit 762ba7698e3ceaa1b17aa299421ddd418f00b823 +Subproject commit c6b039e5c9edb7675f90ffc725f9dd8e66571264 diff --git a/meson.build b/meson.build index 491b1f8ab5b..66dde2dfc9d 100644 --- a/meson.build +++ b/meson.build @@ -95,8 +95,6 @@ if get_option('enable-autodiff') and not omp codi_rev_args += '-DCODI_PRIMAL_REUSE_TAPE' elif get_option('codi-tape') == 'PrimalMultiUse' codi_rev_args += '-DCODI_PRIMAL_MULTIUSE_TAPE' - elif get_option('codi-tape') == 'Tag' - codi_rev_args += '-DCODI_TAG_TAPE' else error('Invalid CoDiPack tape choice @0@'.format(get_option('codi-tape'))) endif diff --git a/meson_options.txt b/meson_options.txt index 6bbc43a928a..08fbac80669 100644 --- a/meson_options.txt +++ b/meson_options.txt @@ -23,7 +23,7 @@ option('enable-coolprop', type : 'boolean', value : false, description: 'enable option('enable-mlpcpp', type : 'boolean', value : false, description: 'enable profiling through gprof') option('enable-gprof', type : 'boolean', value : false, description: 'enable MLPCpp support') option('opdi-backend', type : 'combo', choices : ['auto', 'macro', 'ompt'], value : 'auto', description: 'OpDiLib backend choice') -option('codi-tape', type : 'combo', choices : ['JacobianLinear', 'JacobianReuse', 'JacobianMultiUse', 'PrimalLinear', 'PrimalReuse', 'PrimalMultiUse', 'Tag'], value : 'JacobianLinear', description: 'CoDiPack tape choice') +option('codi-tape', type : 'combo', choices : ['JacobianLinear', 'JacobianReuse', 'JacobianMultiUse', 'PrimalLinear', 'PrimalReuse', 'PrimalMultiUse'], value : 'JacobianLinear', description: 'CoDiPack tape choice') option('opdi-shared-read-opt', type : 'boolean', value : true, description : 'OpDiLib shared reading optimization') option('librom_root', type : 'string', value : '', description: 'libROM base directory') option('enable-librom', type : 'boolean', value : false, description: 'enable LLNL libROM support') From 5a9773c4d501bb3d026353c8b542ae9fa6798078 Mon Sep 17 00:00:00 2001 From: Josh Kelly Date: Wed, 4 Sep 2024 10:19:16 +0200 Subject: [PATCH 10/22] formatting --- SU2_CFD/src/integration/CIntegration.cpp | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/SU2_CFD/src/integration/CIntegration.cpp b/SU2_CFD/src/integration/CIntegration.cpp index 65406ae4e16..2abb1bb197b 100644 --- a/SU2_CFD/src/integration/CIntegration.cpp +++ b/SU2_CFD/src/integration/CIntegration.cpp @@ -85,14 +85,13 @@ void CIntegration::Space_Integration(CGeometry *geometry, if (config->GetBoolGiles() && config->GetSpatialFourier()){ solver_container[MainSolver]->PreprocessBC_Giles(geometry, config, conv_bound_numerics, INFLOW); - solver_container[MainSolver]->PreprocessBC_Giles(geometry, config, conv_bound_numerics, OUTFLOW); } BEGIN_SU2_OMP_SAFE_GLOBAL_ACCESS { if (config->GetBoolTurbomachinery()){ /*--- Average quantities at the inflow and outflow boundaries ---*/ - solver_container[MainSolver]->TurboAverageProcess(solver_container, geometry,config,INFLOW); + solver_container[MainSolver]->TurboAverageProcess(solver_container, geometry,config, INFLOW); solver_container[MainSolver]->TurboAverageProcess(solver_container, geometry, config, OUTFLOW); } } From 6e65dd17bae693cb21c2b893f2123fcd9d25e722 Mon Sep 17 00:00:00 2001 From: Ole Burghardt Date: Tue, 24 Sep 2024 17:26:02 +0200 Subject: [PATCH 11/22] Add config option and indicator for a debug discrete adjoint run. --- Common/include/CConfig.hpp | 9 ++++++++- Common/src/CConfig.cpp | 10 +++++++++- SU2_CFD/src/SU2_CFD.cpp | 1 + 3 files changed, 18 insertions(+), 2 deletions(-) diff --git a/Common/include/CConfig.hpp b/Common/include/CConfig.hpp index df9a0dba3b7..aabfc9e344d 100644 --- a/Common/include/CConfig.hpp +++ b/Common/include/CConfig.hpp @@ -1043,7 +1043,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. */ @@ -8856,6 +8857,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/src/CConfig.cpp b/Common/src/CConfig.cpp index ec4a555d6af..f6a96e5e626 100644 --- a/Common/src/CConfig.cpp +++ b/Common/src/CConfig.cpp @@ -1104,11 +1104,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*/ 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); From 82f69e919e5fff3ff2fe79159977134d11340b50 Mon Sep 17 00:00:00 2001 From: Ole Burghardt Date: Tue, 24 Sep 2024 17:29:45 +0200 Subject: [PATCH 12/22] Add routines for setting tags to AD structure. --- Common/include/basic_types/ad_structure.hpp | 8 ++++++++ 1 file changed, 8 insertions(+) diff --git a/Common/include/basic_types/ad_structure.hpp b/Common/include/basic_types/ad_structure.hpp index ecd2c93b1c1..d3d6cd7c945 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. From 9de8e93216358188e785b6bf7e8033b97929556a Mon Sep 17 00:00:00 2001 From: Ole Burghardt Date: Tue, 24 Sep 2024 17:36:58 +0200 Subject: [PATCH 13/22] Add kinds of recording for tag tapes. --- Common/include/option_structure.hpp | 2 ++ .../include/drivers/CDiscAdjMultizoneDriver.hpp | 2 +- SU2_CFD/src/drivers/CDiscAdjMultizoneDriver.cpp | 17 +++++++++++++++-- .../src/iteration/CDiscAdjFluidIteration.cpp | 5 ++++- 4 files changed, 22 insertions(+), 4 deletions(-) diff --git a/Common/include/option_structure.hpp b/Common/include/option_structure.hpp index 64c781a34f8..659107bdc00 100644 --- a/Common/include/option_structure.hpp +++ b/Common/include/option_structure.hpp @@ -2476,6 +2476,8 @@ enum class RECORDING { MESH_COORDS, MESH_DEFORM, SOLUTION_AND_MESH, + TAG_INIT_SOLUTION_VARIABLES, + TAG_CHECK_SOLUTION_VARIABLES }; /*! diff --git a/SU2_CFD/include/drivers/CDiscAdjMultizoneDriver.hpp b/SU2_CFD/include/drivers/CDiscAdjMultizoneDriver.hpp index f7c8929c805..67d476dddb1 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 diff --git a/SU2_CFD/src/drivers/CDiscAdjMultizoneDriver.cpp b/SU2_CFD/src/drivers/CDiscAdjMultizoneDriver.cpp index 6a8c51fc0fa..9885df474cb 100644 --- a/SU2_CFD/src/drivers/CDiscAdjMultizoneDriver.cpp +++ b/SU2_CFD/src/drivers/CDiscAdjMultizoneDriver.cpp @@ -597,6 +597,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,13 +747,17 @@ 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"; + cout << " (" << ObjFunc_Index << ")\n"; } cout << endl; } diff --git a/SU2_CFD/src/iteration/CDiscAdjFluidIteration.cpp b/SU2_CFD/src/iteration/CDiscAdjFluidIteration.cpp index 27173f3357e..40b75819b81 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()) { From 57fbcea154174efb8b17bc65debefcd19ddb4ca7 Mon Sep 17 00:00:00 2001 From: Ole Burghardt Date: Tue, 24 Sep 2024 17:41:30 +0200 Subject: [PATCH 14/22] Comment out checks in AD structure whether an identifier is active. This is a temporary change to avoid false positives when a tag tape is used. They might be unnecessary anyway though. --- Common/include/basic_types/ad_structure.hpp | 22 +++++++++++---------- 1 file changed, 12 insertions(+), 10 deletions(-) diff --git a/Common/include/basic_types/ad_structure.hpp b/Common/include/basic_types/ad_structure.hpp index d3d6cd7c945..5861faccd1c 100644 --- a/Common/include/basic_types/ad_structure.hpp +++ b/Common/include/basic_types/ad_structure.hpp @@ -522,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...); } @@ -535,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]); - } + // } } } } @@ -547,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]); - } + // } } } } @@ -567,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...); } @@ -575,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]); - } + // } } } } @@ -587,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]); - } + // } } } } From 3da17d81ff64a3be01bf9372c1f1b2538cba9d0d Mon Sep 17 00:00:00 2001 From: Ole Burghardt Date: Tue, 24 Sep 2024 17:54:47 +0200 Subject: [PATCH 15/22] Add first WIP version of the tag debug mode to the (multizone) discrete adjoint solver. --- .../drivers/CDiscAdjMultizoneDriver.hpp | 5 +++ SU2_CFD/include/drivers/CDriver.hpp | 5 +++ .../src/drivers/CDiscAdjMultizoneDriver.cpp | 32 +++++++++++++++++++ 3 files changed, 42 insertions(+) diff --git a/SU2_CFD/include/drivers/CDiscAdjMultizoneDriver.hpp b/SU2_CFD/include/drivers/CDiscAdjMultizoneDriver.hpp index 67d476dddb1..1f5ad0f21f5 100644 --- a/SU2_CFD/include/drivers/CDiscAdjMultizoneDriver.hpp +++ b/SU2_CFD/include/drivers/CDiscAdjMultizoneDriver.hpp @@ -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/src/drivers/CDiscAdjMultizoneDriver.cpp b/SU2_CFD/src/drivers/CDiscAdjMultizoneDriver.cpp index 9885df474cb..63fdf9acfd3 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); From b0e7c86c5d0f27ac6b0bca2315d922503418a1be Mon Sep 17 00:00:00 2001 From: Ole Burghardt Date: Tue, 24 Sep 2024 17:56:47 +0200 Subject: [PATCH 16/22] Fix preaccumulation error in ComputeVorticityAndStrainMag (first debug mode finding :-)) --- SU2_CFD/include/solvers/CFVMFlowSolverBase.inl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/SU2_CFD/include/solvers/CFVMFlowSolverBase.inl b/SU2_CFD/include/solvers/CFVMFlowSolverBase.inl index f60db438e18..10ba11296f2 100644 --- a/SU2_CFD/include/solvers/CFVMFlowSolverBase.inl +++ b/SU2_CFD/include/solvers/CFVMFlowSolverBase.inl @@ -672,11 +672,11 @@ void 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 From 8556c2ae8fcfa57faa971650fee11311cdc58255 Mon Sep 17 00:00:00 2001 From: Ole Burghardt Date: Tue, 24 Sep 2024 18:15:39 +0200 Subject: [PATCH 17/22] Revert small change (will be added to CoDi soon). --- SU2_CFD/src/drivers/CDiscAdjMultizoneDriver.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/SU2_CFD/src/drivers/CDiscAdjMultizoneDriver.cpp b/SU2_CFD/src/drivers/CDiscAdjMultizoneDriver.cpp index 63fdf9acfd3..d70619136d2 100644 --- a/SU2_CFD/src/drivers/CDiscAdjMultizoneDriver.cpp +++ b/SU2_CFD/src/drivers/CDiscAdjMultizoneDriver.cpp @@ -789,7 +789,7 @@ void CDiscAdjMultizoneDriver::SetObjFunction(RECORDING kind_recording) { kind_recording == RECORDING::TAG_CHECK_SOLUTION_VARIABLES) { cout << " Objective function : " << ObjFunc; if (driver_config->GetWrt_AD_Statistics()){ - cout << " (" << ObjFunc_Index << ")\n"; + // cout << " (" << ObjFunc_Index << ")\n"; } cout << endl; } From 683de7777cd765301038b2adc9ecd12f910b5e91 Mon Sep 17 00:00:00 2001 From: Josh Kelly Date: Thu, 26 Sep 2024 17:48:52 +0200 Subject: [PATCH 18/22] Resolves tag errors in turbo mode --- SU2_CFD/src/drivers/CDiscAdjMultizoneDriver.cpp | 2 +- SU2_CFD/src/iteration/CDiscAdjFluidIteration.cpp | 3 +++ SU2_CFD/src/solvers/CDiscAdjSolver.cpp | 6 +++--- externals/codi | 2 +- 4 files changed, 8 insertions(+), 5 deletions(-) diff --git a/SU2_CFD/src/drivers/CDiscAdjMultizoneDriver.cpp b/SU2_CFD/src/drivers/CDiscAdjMultizoneDriver.cpp index d70619136d2..75e54ab30f1 100644 --- a/SU2_CFD/src/drivers/CDiscAdjMultizoneDriver.cpp +++ b/SU2_CFD/src/drivers/CDiscAdjMultizoneDriver.cpp @@ -789,7 +789,7 @@ void CDiscAdjMultizoneDriver::SetObjFunction(RECORDING kind_recording) { kind_recording == RECORDING::TAG_CHECK_SOLUTION_VARIABLES) { cout << " Objective function : " << ObjFunc; if (driver_config->GetWrt_AD_Statistics()){ - // cout << " (" << ObjFunc_Index << ")\n"; +// cout << " (" << ObjFunc_Index << ")\n"; } cout << endl; } diff --git a/SU2_CFD/src/iteration/CDiscAdjFluidIteration.cpp b/SU2_CFD/src/iteration/CDiscAdjFluidIteration.cpp index e7e6feb8e78..a5ae846a686 100644 --- a/SU2_CFD/src/iteration/CDiscAdjFluidIteration.cpp +++ b/SU2_CFD/src/iteration/CDiscAdjFluidIteration.cpp @@ -475,8 +475,11 @@ void CDiscAdjFluidIteration::SetDependencies(CSolver***** solver, CGeometry**** solvers0[FLOW_SOL]->CompleteComms(geometry0, config[iZone], MPI_QUANTITIES::SOLUTION); if (config[iZone]->GetBoolTurbomachinery()) { + solvers0[FLOW_SOL]->PreprocessAverage(solvers0, geometry0, config[iZone], INFLOW); + solvers0[FLOW_SOL]->PreprocessAverage(solvers0, geometry0, config[iZone], OUTFLOW); solvers0[FLOW_SOL]->TurboAverageProcess(solvers0, geometry0, config[iZone], INFLOW); solvers0[FLOW_SOL]->TurboAverageProcess(solvers0, geometry0, config[iZone], OUTFLOW); + solvers0[FLOW_SOL]->GatherInOutAverageValues(config[iZone], geometry0); } if (turbulent && !config[iZone]->GetFrozen_Visc_Disc()) { solvers0[TURB_SOL]->Postprocessing(geometry0, solvers0, diff --git a/SU2_CFD/src/solvers/CDiscAdjSolver.cpp b/SU2_CFD/src/solvers/CDiscAdjSolver.cpp index 1b1989fed7d..39b8fbc2bfb 100644 --- a/SU2_CFD/src/solvers/CDiscAdjSolver.cpp +++ b/SU2_CFD/src/solvers/CDiscAdjSolver.cpp @@ -219,14 +219,14 @@ void CDiscAdjSolver::RegisterVariables(CGeometry *geometry, CConfig *config, boo if ((config->GetKind_Regime() == ENUM_REGIME::COMPRESSIBLE) && (KindDirect_Solver == RUNTIME_FLOW_SYS) && config->GetBoolTurbomachinery()){ - BPressure = config->GetPressureOut_BC(); - Temperature = config->GetTotalTemperatureIn_BC(); - if (!reset){ AD::RegisterInput(BPressure); AD::RegisterInput(Temperature); } + BPressure = config->GetPressureOut_BC(); + Temperature = config->GetTotalTemperatureIn_BC(); + config->SetPressureOut_BC(BPressure); config->SetTotalTemperatureIn_BC(Temperature); } diff --git a/externals/codi b/externals/codi index 762ba7698e3..209e13df41e 160000 --- a/externals/codi +++ b/externals/codi @@ -1 +1 @@ -Subproject commit 762ba7698e3ceaa1b17aa299421ddd418f00b823 +Subproject commit 209e13df41e28fa7e24555e62ea9f044f6a8c9bb From 13ca3f5f3849fa8800596e00be9731314bd176ed Mon Sep 17 00:00:00 2001 From: Josh Kelly Date: Thu, 26 Sep 2024 17:52:32 +0200 Subject: [PATCH 19/22] Revert "Resolves tag errors in turbo mode" This reverts commit 683de7777cd765301038b2adc9ecd12f910b5e91. Working on the wrong branch. --- SU2_CFD/src/drivers/CDiscAdjMultizoneDriver.cpp | 2 +- SU2_CFD/src/iteration/CDiscAdjFluidIteration.cpp | 3 --- SU2_CFD/src/solvers/CDiscAdjSolver.cpp | 6 +++--- externals/codi | 2 +- 4 files changed, 5 insertions(+), 8 deletions(-) diff --git a/SU2_CFD/src/drivers/CDiscAdjMultizoneDriver.cpp b/SU2_CFD/src/drivers/CDiscAdjMultizoneDriver.cpp index 75e54ab30f1..d70619136d2 100644 --- a/SU2_CFD/src/drivers/CDiscAdjMultizoneDriver.cpp +++ b/SU2_CFD/src/drivers/CDiscAdjMultizoneDriver.cpp @@ -789,7 +789,7 @@ void CDiscAdjMultizoneDriver::SetObjFunction(RECORDING kind_recording) { kind_recording == RECORDING::TAG_CHECK_SOLUTION_VARIABLES) { cout << " Objective function : " << ObjFunc; if (driver_config->GetWrt_AD_Statistics()){ -// cout << " (" << ObjFunc_Index << ")\n"; + // cout << " (" << ObjFunc_Index << ")\n"; } cout << endl; } diff --git a/SU2_CFD/src/iteration/CDiscAdjFluidIteration.cpp b/SU2_CFD/src/iteration/CDiscAdjFluidIteration.cpp index a5ae846a686..e7e6feb8e78 100644 --- a/SU2_CFD/src/iteration/CDiscAdjFluidIteration.cpp +++ b/SU2_CFD/src/iteration/CDiscAdjFluidIteration.cpp @@ -475,11 +475,8 @@ void CDiscAdjFluidIteration::SetDependencies(CSolver***** solver, CGeometry**** solvers0[FLOW_SOL]->CompleteComms(geometry0, config[iZone], MPI_QUANTITIES::SOLUTION); if (config[iZone]->GetBoolTurbomachinery()) { - solvers0[FLOW_SOL]->PreprocessAverage(solvers0, geometry0, config[iZone], INFLOW); - solvers0[FLOW_SOL]->PreprocessAverage(solvers0, geometry0, config[iZone], OUTFLOW); solvers0[FLOW_SOL]->TurboAverageProcess(solvers0, geometry0, config[iZone], INFLOW); solvers0[FLOW_SOL]->TurboAverageProcess(solvers0, geometry0, config[iZone], OUTFLOW); - solvers0[FLOW_SOL]->GatherInOutAverageValues(config[iZone], geometry0); } if (turbulent && !config[iZone]->GetFrozen_Visc_Disc()) { solvers0[TURB_SOL]->Postprocessing(geometry0, solvers0, diff --git a/SU2_CFD/src/solvers/CDiscAdjSolver.cpp b/SU2_CFD/src/solvers/CDiscAdjSolver.cpp index 39b8fbc2bfb..1b1989fed7d 100644 --- a/SU2_CFD/src/solvers/CDiscAdjSolver.cpp +++ b/SU2_CFD/src/solvers/CDiscAdjSolver.cpp @@ -219,14 +219,14 @@ void CDiscAdjSolver::RegisterVariables(CGeometry *geometry, CConfig *config, boo if ((config->GetKind_Regime() == ENUM_REGIME::COMPRESSIBLE) && (KindDirect_Solver == RUNTIME_FLOW_SYS) && config->GetBoolTurbomachinery()){ + BPressure = config->GetPressureOut_BC(); + Temperature = config->GetTotalTemperatureIn_BC(); + if (!reset){ AD::RegisterInput(BPressure); AD::RegisterInput(Temperature); } - BPressure = config->GetPressureOut_BC(); - Temperature = config->GetTotalTemperatureIn_BC(); - config->SetPressureOut_BC(BPressure); config->SetTotalTemperatureIn_BC(Temperature); } diff --git a/externals/codi b/externals/codi index 209e13df41e..762ba7698e3 160000 --- a/externals/codi +++ b/externals/codi @@ -1 +1 @@ -Subproject commit 209e13df41e28fa7e24555e62ea9f044f6a8c9bb +Subproject commit 762ba7698e3ceaa1b17aa299421ddd418f00b823 From 887e9e6eb986d55bf6d44471ce5eeea002a75d93 Mon Sep 17 00:00:00 2001 From: Josh Kelly Date: Thu, 26 Sep 2024 18:14:08 +0200 Subject: [PATCH 20/22] Resolving tag errors --- SU2_CFD/src/iteration/CDiscAdjFluidIteration.cpp | 3 +++ SU2_CFD/src/solvers/CDiscAdjSolver.cpp | 6 +++--- externals/codi | 2 +- 3 files changed, 7 insertions(+), 4 deletions(-) diff --git a/SU2_CFD/src/iteration/CDiscAdjFluidIteration.cpp b/SU2_CFD/src/iteration/CDiscAdjFluidIteration.cpp index 24fd1602e6e..78b227c3ead 100644 --- a/SU2_CFD/src/iteration/CDiscAdjFluidIteration.cpp +++ b/SU2_CFD/src/iteration/CDiscAdjFluidIteration.cpp @@ -472,8 +472,11 @@ void CDiscAdjFluidIteration::SetDependencies(CSolver***** solver, CGeometry**** solvers0[FLOW_SOL]->CompleteComms(geometry0, config[iZone], MPI_QUANTITIES::SOLUTION); if (config[iZone]->GetBoolTurbomachinery()) { + solvers0[FLOW_SOL]->PreprocessAverage(solvers0, geometry0, config[iZone], INFLOW); + solvers0[FLOW_SOL]->PreprocessAverage(solvers0, geometry0, config[iZone], OUTFLOW); solvers0[FLOW_SOL]->TurboAverageProcess(solvers0, geometry0, config[iZone], INFLOW); solvers0[FLOW_SOL]->TurboAverageProcess(solvers0, geometry0, config[iZone], OUTFLOW); + solvers0[FLOW_SOL]->GatherInOutAverageValues(config[iZone], geometry0); } if (turbulent && !config[iZone]->GetFrozen_Visc_Disc()) { solvers0[TURB_SOL]->Postprocessing(geometry0, solvers0, diff --git a/SU2_CFD/src/solvers/CDiscAdjSolver.cpp b/SU2_CFD/src/solvers/CDiscAdjSolver.cpp index 1b1989fed7d..39b8fbc2bfb 100644 --- a/SU2_CFD/src/solvers/CDiscAdjSolver.cpp +++ b/SU2_CFD/src/solvers/CDiscAdjSolver.cpp @@ -219,14 +219,14 @@ void CDiscAdjSolver::RegisterVariables(CGeometry *geometry, CConfig *config, boo if ((config->GetKind_Regime() == ENUM_REGIME::COMPRESSIBLE) && (KindDirect_Solver == RUNTIME_FLOW_SYS) && config->GetBoolTurbomachinery()){ - BPressure = config->GetPressureOut_BC(); - Temperature = config->GetTotalTemperatureIn_BC(); - if (!reset){ AD::RegisterInput(BPressure); AD::RegisterInput(Temperature); } + BPressure = config->GetPressureOut_BC(); + Temperature = config->GetTotalTemperatureIn_BC(); + config->SetPressureOut_BC(BPressure); config->SetTotalTemperatureIn_BC(Temperature); } diff --git a/externals/codi b/externals/codi index c6b039e5c9e..209e13df41e 160000 --- a/externals/codi +++ b/externals/codi @@ -1 +1 @@ -Subproject commit c6b039e5c9edb7675f90ffc725f9dd8e66571264 +Subproject commit 209e13df41e28fa7e24555e62ea9f044f6a8c9bb From e40f94319be4fde929246467f95ad43bc1937759 Mon Sep 17 00:00:00 2001 From: Josh Kelly Date: Fri, 27 Sep 2024 15:42:16 +0200 Subject: [PATCH 21/22] Revert "Merge remote-tracking branch 'origin/feature_tag_debug_tape' into feature_mz_adjoint_for_turbo" This reverts commit 89c95ba447726665a3ab109d15c4a7a24c409ac1, reversing changes made to c002911246cce5d2debcec36a77d0cea3e73a7a9. --- .github/workflows/regression.yml | 2 +- Common/include/CConfig.hpp | 49 +- Common/include/basic_types/ad_structure.hpp | 30 +- Common/include/option_structure.hpp | 20 +- Common/include/option_structure.inl | 16 +- Common/src/CConfig.cpp | 78 +-- Common/src/geometry/CPhysicalGeometry.cpp | 3 - .../drivers/CDiscAdjMultizoneDriver.hpp | 7 +- SU2_CFD/include/drivers/CDriver.hpp | 5 - SU2_CFD/include/fluid/CFluidScalar.hpp | 3 +- SU2_CFD/include/interfaces/CInterface.hpp | 9 + .../numerics/turbulent/turb_sources.hpp | 56 +- SU2_CFD/include/solvers/CEulerSolver.hpp | 35 +- .../include/solvers/CFVMFlowSolverBase.inl | 4 +- SU2_CFD/include/solvers/CSolver.hpp | 39 +- SU2_CFD/src/SU2_CFD.cpp | 1 - .../src/drivers/CDiscAdjMultizoneDriver.cpp | 47 +- SU2_CFD/src/drivers/CDriver.cpp | 25 +- SU2_CFD/src/drivers/CMultizoneDriver.cpp | 1 - SU2_CFD/src/drivers/CSinglezoneDriver.cpp | 4 +- SU2_CFD/src/fluid/CFluidScalar.cpp | 5 +- SU2_CFD/src/interfaces/CInterface.cpp | 10 +- .../interfaces/cfd/CMixingPlaneInterface.cpp | 2 +- .../src/iteration/CDiscAdjFluidIteration.cpp | 5 +- SU2_CFD/src/output/CFlowCompOutput.cpp | 8 - SU2_CFD/src/output/CFlowOutput.cpp | 18 +- SU2_CFD/src/output/COutput.cpp | 24 +- SU2_CFD/src/output/CTurboOutput.cpp | 2 + SU2_CFD/src/python_wrapper_structure.cpp | 2 +- SU2_CFD/src/solvers/CEulerSolver.cpp | 555 +++++++----------- SU2_CFD/src/solvers/CIncEulerSolver.cpp | 4 +- SU2_CFD/src/solvers/CSolver.cpp | 16 +- TestCases/hybrid_regression.py | 11 +- TestCases/parallel_regression.py | 13 +- TestCases/serial_regression.py | 13 +- .../multi_interface/multi_interface_rst.cfg | 151 ----- .../turbomachinery/multi_interface/zone_1.cfg | 6 - .../turbomachinery/multi_interface/zone_2.cfg | 6 - .../turbomachinery/multi_interface/zone_3.cfg | 6 - config_template.cfg | 3 - 40 files changed, 372 insertions(+), 922 deletions(-) delete mode 100644 TestCases/turbomachinery/multi_interface/multi_interface_rst.cfg delete mode 100644 TestCases/turbomachinery/multi_interface/zone_1.cfg delete mode 100644 TestCases/turbomachinery/multi_interface/zone_2.cfg delete mode 100644 TestCases/turbomachinery/multi_interface/zone_3.cfg diff --git a/.github/workflows/regression.yml b/.github/workflows/regression.yml index 2502204b64c..cd916ffcfe5 100644 --- a/.github/workflows/regression.yml +++ b/.github/workflows/regression.yml @@ -209,7 +209,7 @@ jobs: uses: docker://ghcr.io/su2code/su2/test-su2:240320-1536 with: # -t -c - args: -b ${{github.ref}} -t develop -c develop -s ${{matrix.testscript}} + args: -b ${{github.ref}} -t feature_solid-solid_cht -c develop -s ${{matrix.testscript}} - name: Cleanup uses: docker://ghcr.io/su2code/su2/test-su2:240320-1536 with: diff --git a/Common/include/CConfig.hpp b/Common/include/CConfig.hpp index 71878276052..af0a9916d98 100644 --- a/Common/include/CConfig.hpp +++ b/Common/include/CConfig.hpp @@ -138,7 +138,6 @@ 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). */ @@ -237,7 +236,6 @@ 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. */ @@ -444,7 +442,6 @@ class CConfig { TURBO_PERF_KIND *Kind_TurboPerf; /*!< \brief Kind of turbomachynery architecture.*/ TURBOMACHINERY_TYPE *Kind_TurboMachinery; - su2vector Kind_TurboInterface; /* Turbomachinery objective functions */ su2double *EntropyGeneration; @@ -473,7 +470,6 @@ 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. */ @@ -754,7 +750,6 @@ 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. */ @@ -772,7 +767,6 @@ 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. */ @@ -1057,8 +1051,7 @@ 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. */ - DiscreteAdjointDebug; /*!< \brief Discrete adjoint debug mode using tags. */ + bool DiscreteAdjoint; /*!< \brief AD-based discrete adjoint mode. */ 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. */ @@ -1386,7 +1379,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_Turbomachinery); + string* & Marker_TurboBoundIn, string* & Marker_TurboBoundOut); void addActDiskOption(const string & name, unsigned short & nMarker_ActDiskInlet, unsigned short & nMarker_ActDiskOutlet, string* & Marker_ActDiskInlet, string* & Marker_ActDiskOutlet, @@ -3524,13 +3517,6 @@ 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). @@ -3677,13 +3663,6 @@ 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. @@ -5360,12 +5339,6 @@ 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. @@ -6429,12 +6402,6 @@ 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. @@ -6553,12 +6520,6 @@ 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. @@ -8950,12 +8911,6 @@ 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 ff2125dcfe3..c44374f7eb7 100644 --- a/Common/include/basic_types/ad_structure.hpp +++ b/Common/include/basic_types/ad_structure.hpp @@ -261,12 +261,6 @@ inline void SetExtFuncOut(T&& data, const int size_x, const int size_y) {} */ inline void SetIndex(int& 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. */ @@ -478,8 +472,6 @@ FORCEINLINE void EndUseAdjoints() { AD::getTape().endUseAdjointVector(); } FORCEINLINE void SetIndex(int& 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. @@ -510,8 +502,7 @@ 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...); } @@ -524,9 +515,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]); - // } + } } } } @@ -536,9 +527,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]); - // } + } } } } @@ -556,8 +547,7 @@ 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...); } @@ -565,9 +555,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]); - // } + } } } } @@ -577,9 +567,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 d456a587460..692dad07331 100644 --- a/Common/include/option_structure.hpp +++ b/Common/include/option_structure.hpp @@ -1622,7 +1622,6 @@ 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. */ @@ -1753,8 +1752,7 @@ enum RIEMANN_TYPE { TOTAL_CONDITIONS_PT_1D = 11, STATIC_PRESSURE_1D = 12, MIXING_IN_1D = 13, - MIXING_OUT_1D = 14, - MASS_FLOW_OUTLET = 15 + MIXING_OUT_1D =14 }; static const MapType Riemann_Map = { MakePair("TOTAL_CONDITIONS_PT", TOTAL_CONDITIONS_PT) @@ -1771,7 +1769,6 @@ 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 = { @@ -1789,7 +1786,6 @@ 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) }; /*! @@ -1866,18 +1862,6 @@ 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. */ @@ -2499,8 +2483,6 @@ 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 01ed6de9bf9..955d1531482 100644 --- a/Common/include/option_structure.inl +++ b/Common/include/option_structure.inl @@ -1606,15 +1606,11 @@ 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, string*& Marker_Turbomachinery) - : size(nMarker_TurboPerf), - marker_turboIn(Marker_TurboBoundIn), - marker_turboOut(Marker_TurboBoundOut), - markers(Marker_Turbomachinery) { + string*& Marker_TurboBoundIn, string*& Marker_TurboBoundOut) + : size(nMarker_TurboPerf), marker_turboIn(Marker_TurboBoundIn), marker_turboOut(Marker_TurboBoundOut) { this->name = option_field_name; } @@ -1628,7 +1624,6 @@ class COptionTurboPerformance : public COptionBase { this->size = 0; this->marker_turboIn = nullptr; this->marker_turboOut = nullptr; - this->markers = nullptr; return ""; } @@ -1639,16 +1634,10 @@ 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]; @@ -1665,7 +1654,6 @@ 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 01b4094d80b..7da2c3f0de6 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_Turbomachinery) { + string* & Marker_TurboBoundIn, string* & Marker_TurboBoundOut) { 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, Marker_Turbomachinery); + COptionBase* val = new COptionTurboPerformance(name, nMarker_TurboPerf, Marker_TurboBoundIn, Marker_TurboBoundOut); option_map.insert(pair(name, val)); } @@ -1037,7 +1037,6 @@ void CConfig::SetPointersNull() { Marker_MixingPlaneInterface = nullptr; Marker_TurboBoundIn = nullptr; Marker_TurboBoundOut = nullptr; - Marker_Turbomachinery = nullptr; Marker_Giles = nullptr; Marker_Shroud = nullptr; @@ -1110,19 +1109,11 @@ 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); -#if defined (CODI_REVERSE_TYPE) +#ifdef 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*/ @@ -1647,8 +1638,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 boundaries for which the turbomachinery settings are applied. \ingroup Config*/ - addTurboPerfOption("MARKER_TURBOMACHINERY", nMarker_Turbomachinery, Marker_TurboBoundIn, Marker_TurboBoundOut, Marker_Turbomachinery); + /*!\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 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. @@ -1750,8 +1741,6 @@ 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 */ @@ -5686,7 +5675,6 @@ 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++) { @@ -5712,7 +5700,6 @@ 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] (); @@ -5833,7 +5820,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] = BC_TYPE::FLUID_INTERFACE; + Marker_CfgFile_KindBC[iMarker_CfgFile] = FLUID_INTERFACE; iMarker_CfgFile++; } @@ -6064,17 +6051,6 @@ 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++) { @@ -6086,37 +6062,6 @@ 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++) @@ -8010,13 +7955,6 @@ 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++) @@ -8204,9 +8142,6 @@ 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; @@ -8347,7 +8282,6 @@ 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 b9b0a0f36cb..a069424c666 100644 --- a/Common/src/geometry/CPhysicalGeometry.cpp +++ b/Common/src/geometry/CPhysicalGeometry.cpp @@ -3370,7 +3370,6 @@ 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)); @@ -3395,7 +3394,6 @@ 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); @@ -3766,7 +3764,6 @@ 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 244d70302ad..13da05669bd 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. + * \brief Kinds of recordings (three different ones). */ enum class Kind_Tape { FULL_TAPE, /*!< \brief Entire derivative information for a coupled adjoint @@ -141,11 +141,6 @@ 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 855f5955f2d..c5e78bd9037 100644 --- a/SU2_CFD/include/drivers/CDriver.hpp +++ b/SU2_CFD/include/drivers/CDriver.hpp @@ -425,11 +425,6 @@ 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 e3295d99b6d..9f5c9108c51 100644 --- a/SU2_CFD/include/fluid/CFluidScalar.hpp +++ b/SU2_CFD/include/fluid/CFluidScalar.hpp @@ -41,6 +41,7 @@ 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.*/ @@ -105,7 +106,7 @@ class CFluidScalar final : public CFluidModel { /*! * \brief Constructor of the class. */ - CFluidScalar(su2double val_operating_pressure, const CConfig* config); + CFluidScalar(su2double val_Cp, su2double val_gas_constant, 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 3f85e232ce9..9e40eb769d0 100644 --- a/SU2_CFD/include/interfaces/CInterface.hpp +++ b/SU2_CFD/include/interfaces/CInterface.hpp @@ -159,6 +159,15 @@ 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. diff --git a/SU2_CFD/include/numerics/turbulent/turb_sources.hpp b/SU2_CFD/include/numerics/turbulent/turb_sources.hpp index f59da30719b..ce75b38f115 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 = 2.0; + const su2double CRot = 1.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, Shat, d_Shat, fv1, d_fv1, fv2, d_fv2, Prod; + 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; /*--- List of helpers ---*/ su2double Omega, dist_i_2, inv_k2_d2, inv_Shat, g_6, norm2_Grad; @@ -151,7 +151,20 @@ 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; @@ -175,24 +188,12 @@ 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); - /*--- Evaluate Omega with a rotational correction term. ---*/ - - Omega::get(Vorticity_i, nDim, PrimVar_Grad_i + idx.Velocity(), var); + /*--- Compute ft2 term ---*/ + ft2::get(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); @@ -233,6 +234,7 @@ 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); @@ -345,12 +347,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.Omega) { - var.Shat = var.Omega + Sbar; + if (Sbar >= - c2 * var.S) { + var.Shat = var.S + Sbar; } else { - 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; + 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; } if (var.Shat <= 1e-10) { var.Shat = 1e-10; @@ -364,12 +366,12 @@ struct Bsl { /*! \brief Edward. */ struct Edw { static void get(const su2double& nue, const su2double& nu, CSAVariables& var) { - var.Shat = max(var.Omega * ((1.0 / max(var.Ji, 1.0e-16)) + var.fv1), 1.0e-16); + var.Shat = max(var.S * ((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.Omega * pow(var.Ji, -2) / nu + var.Omega * var.d_fv1; + var.d_Shat = -var.S * pow(var.Ji, -2) / nu + var.S * var.d_fv1; } } }; @@ -381,7 +383,7 @@ struct Neg { // Baseline solution Bsl::get(nue, nu, var); } else { - var.Shat = var.Omega; + var.Shat = 1.0e-10; var.d_Shat = 0.0; } /*--- Don't check whether Sbar <>= -cv2*S. @@ -445,8 +447,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.Prod * nue; - jacobian += factor * (-var.Prod * nue * var.d_ft2 + (1.0 - var.ft2) * (nue * var.d_Shat + var.Prod)); + 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)); } static void ComputeDestruction(const su2double& nue, const CSAVariables& var, su2double& destruction, @@ -479,7 +481,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.Prod; + const su2double dP_dnu = var.intermittency * var.cb1 * (1.0 - var.ct3) * var.S; production = dP_dnu * nue; jacobian += dP_dnu; } diff --git a/SU2_CFD/include/solvers/CEulerSolver.hpp b/SU2_CFD/include/solvers/CEulerSolver.hpp index d498b84c225..342448c6253 100644 --- a/SU2_CFD/include/solvers/CEulerSolver.hpp +++ b/SU2_CFD/include/solvers/CEulerSolver.hpp @@ -138,7 +138,6 @@ 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 6c9194f31ae..5eb43dc2035 100644 --- a/SU2_CFD/include/solvers/CSolver.hpp +++ b/SU2_CFD/include/solvers/CSolver.hpp @@ -1089,6 +1089,38 @@ 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. @@ -3799,13 +3831,6 @@ 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 3bf7ba05a11..353a744c891 100644 --- a/SU2_CFD/src/SU2_CFD.cpp +++ b/SU2_CFD/src/SU2_CFD.cpp @@ -102,7 +102,6 @@ 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 63fdf9acfd3..14d84721fd7 100644 --- a/SU2_CFD/src/drivers/CDiscAdjMultizoneDriver.cpp +++ b/SU2_CFD/src/drivers/CDiscAdjMultizoneDriver.cpp @@ -159,16 +159,6 @@ 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. ---*/ @@ -227,28 +217,6 @@ 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); @@ -629,15 +597,6 @@ 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 @@ -779,14 +738,10 @@ void CDiscAdjMultizoneDriver::SetObjFunction(RECORDING kind_recording) { } } - - if (rank == MASTER_NODE) { AD::RegisterOutput(ObjFunc); AD::SetIndex(ObjFunc_Index, ObjFunc); - if (kind_recording == RECORDING::SOLUTION_VARIABLES || - kind_recording == RECORDING::TAG_INIT_SOLUTION_VARIABLES || - kind_recording == RECORDING::TAG_CHECK_SOLUTION_VARIABLES) { + if (kind_recording == RECORDING::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 0c3c4a10c4e..65a25b6c9a2 100644 --- a/SU2_CFD/src/drivers/CDriver.cpp +++ b/SU2_CFD/src/drivers/CDriver.cpp @@ -2476,23 +2476,14 @@ 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) { - /*--- 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; - } + /*--- 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; } } else{ diff --git a/SU2_CFD/src/drivers/CMultizoneDriver.cpp b/SU2_CFD/src/drivers/CMultizoneDriver.cpp index 368a598e514..d69c59e30ae 100644 --- a/SU2_CFD/src/drivers/CMultizoneDriver.cpp +++ b/SU2_CFD/src/drivers/CMultizoneDriver.cpp @@ -590,7 +590,6 @@ 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 627ea9c41f7..6b1bd3eaaaa 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. ---*/ - /*--- or until a SIGTERM signal stops the loop. We catch SIGTERM and exit gracefully ---*/ - while ( TimeIter < config_container[ZONE_0]->GetnTime_Iter()) { + 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 5e3a2e2f998..991e29bd562 100644 --- a/SU2_CFD/src/fluid/CFluidScalar.cpp +++ b/SU2_CFD/src/fluid/CFluidScalar.cpp @@ -41,9 +41,12 @@ #include "../../include/fluid/CPolynomialViscosity.hpp" #include "../../include/fluid/CSutherland.hpp" -CFluidScalar::CFluidScalar(su2double value_pressure_operating, const CConfig* config) +CFluidScalar::CFluidScalar(su2double val_Cp, su2double val_gas_constant, 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 b9617c4c16e..7dc8d1e7f0c 100644 --- a/SU2_CFD/src/interfaces/CInterface.cpp +++ b/SU2_CFD/src/interfaces/CInterface.cpp @@ -253,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; @@ -278,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; - Target_Flag = -1; + Marker_Target = -1; + } if (Marker_Target != -1 && Marker_Donor != -1){ @@ -300,7 +300,7 @@ void CInterface::PreprocessAverage(CGeometry *donor_geometry, CGeometry *target_ } if(test2 < dist2){ dist2 = test2; - tSpan = jSpan; + tSpan =jSpan; } } @@ -638,4 +638,6 @@ 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 d1bbcf2ff0b..f35dc776f9e 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/iteration/CDiscAdjFluidIteration.cpp b/SU2_CFD/src/iteration/CDiscAdjFluidIteration.cpp index a5ae846a686..78b227c3ead 100644 --- a/SU2_CFD/src/iteration/CDiscAdjFluidIteration.cpp +++ b/SU2_CFD/src/iteration/CDiscAdjFluidIteration.cpp @@ -404,10 +404,7 @@ void CDiscAdjFluidIteration::RegisterInput(CSolver***** solver, CGeometry**** ge SU2_OMP_PARALLEL_(if(solvers0[ADJFLOW_SOL]->GetHasHybridParallel())) { - 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) { + if (kind_recording == RECORDING::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 e693bbbae47..5bb5577fe76 100644 --- a/SU2_CFD/src/output/CFlowCompOutput.cpp +++ b/SU2_CFD/src/output/CFlowCompOutput.cpp @@ -592,10 +592,6 @@ 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]); @@ -608,8 +604,6 @@ 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()); @@ -617,8 +611,6 @@ 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 895419f264d..5ffee0d23f5 100644 --- a/SU2_CFD/src/output/CFlowOutput.cpp +++ b/SU2_CFD/src/output/CFlowOutput.cpp @@ -4066,10 +4066,6 @@ 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); @@ -4082,16 +4078,12 @@ 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 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); + 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); } diff --git a/SU2_CFD/src/output/COutput.cpp b/SU2_CFD/src/output/COutput.cpp index 499d0d69fd5..d8aa5c16af9 100644 --- a/SU2_CFD/src/output/COutput.cpp +++ b/SU2_CFD/src/output/COutput.cpp @@ -25,9 +25,6 @@ * License along with SU2. If not, see . */ -#include -#include - #include "../../../Common/include/geometry/CGeometry.hpp" #include "../../include/solvers/CSolver.hpp" @@ -50,15 +47,6 @@ #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()), @@ -181,11 +169,7 @@ COutput::COutput(const CConfig *config, unsigned short ndim, bool fem_output): volumeDataSorter = nullptr; surfaceDataSorter = nullptr; - headerNeeded = false; - - /*--- Setup a signal handler for SIGTERM. ---*/ - - signal(SIGTERM, signalHandler); + headerNeeded = false; } COutput::~COutput() { @@ -249,7 +233,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); @@ -962,13 +946,9 @@ 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 008f63d902d..6551d0062c1 100644 --- a/SU2_CFD/src/output/CTurboOutput.cpp +++ b/SU2_CFD/src/output/CTurboOutput.cpp @@ -244,6 +244,7 @@ void CTurbomachineryStagePerformance::ComputeTurbineStagePerformance(const CTurb fluidModel.SetTDState_Ps(OutState.GetPressure(), InState.GetEntropy()); su2double enthalpyOutIs = fluidModel.GetStaticEnergy() + OutState.GetPressure() / fluidModel.GetDensity(); su2double totEnthalpyOutIs = enthalpyOutIs + 0.5 * OutState.GetVelocityValue() * OutState.GetVelocityValue(); + su2double tangVel = OutState.GetTangVelocity(); su2double relVelOutIs2 = 2 * (OutState.GetRothalpy() - enthalpyOutIs) + tangVel * tangVel; @@ -254,6 +255,7 @@ void CTurbomachineryStagePerformance::ComputeTurbineStagePerformance(const CTurb TotalTotalEfficiency = EulerianWork / (InState.GetTotalEnthalpy() - totEnthalpyOutIs); TotalStaticPressureRatio = InState.GetTotalPressure() / OutState.GetPressure(); TotalTotalPressureRatio = InState.GetTotalPressure() / OutState.GetTotalPressure(); + TotalPressureLoss = (InState.GetTotalRelPressure() - OutState.GetTotalRelPressure()) / (OutState.GetTotalRelPressure() - OutState.GetPressure()); KineticEnergyLoss = 2 * (OutState.GetEnthalpy() - enthalpyOutIs) / relVelOutIs2; diff --git a/SU2_CFD/src/python_wrapper_structure.cpp b/SU2_CFD/src/python_wrapper_structure.cpp index 86933d11341..a0018371443 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 (auto iMesh = 0u; iMesh <= main_config->GetnMGLevels(); iMesh++) { + for (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 286a84e21c7..6267635a30e 100644 --- a/SU2_CFD/src/solvers/CEulerSolver.cpp +++ b/SU2_CFD/src/solvers/CEulerSolver.cpp @@ -232,9 +232,6 @@ 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(); @@ -282,8 +279,6 @@ 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. ---*/ @@ -2546,7 +2541,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->GetBc_Eval_Freq())) == 0)) || (config->GetInnerIter() == 1) || (config->GetDiscrete_Adjoint())); + bool Evaluate_BC = ((((config->GetInnerIter() % (config->GetScreen_Wrt_Freq(2)*40)) == 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; @@ -2935,7 +2930,6 @@ 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]); @@ -6121,81 +6115,60 @@ 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) { - 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++) + /*--- 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++) { R_Matrix[iVar] = new su2double[nVar]; c_avg[iVar] = 0.0; dcjs[iVar] = 0.0; } - for (iVar = 0; iVar < nVar-1; iVar++) + for (auto iVar = 0; iVar < nVar-1; iVar++) { R_c[iVar] = new su2double[nVar-1]; R_c_inv[iVar] = new su2double[nVar-1]; } - - complex I, c2ks, c2js, c3ks, c3js, cOutks, cOutjs, Beta_inf; - I = complex(0.0,1.0); + 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; + if (nDim == 3){ extrarelfacAvg = config->GetExtraRelFacGiles(0); - spanPercent = config->GetExtraRelFacGiles(1); - Turbo_Flag = config->GetMarker_All_TurbomachineryFlag(val_marker); + auto const spanPercent = config->GetExtraRelFacGiles(1); + auto const Turbo_Flag = config->GetMarker_All_TurbomachineryFlag(val_marker); SpanWiseValues = geometry->GetSpanWiseValue(Turbo_Flag); deltaSpan = SpanWiseValues[nSpanWiseSections-1]*spanPercent; coeffrelfacAvg = (relfacAvgCfg - extrarelfacAvg)/deltaSpan; } - for (iSpan= 0; iSpan < nSpanWiseSections ; iSpan++){ + for (unsigned short 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){ @@ -6208,18 +6181,16 @@ 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]); - AverageSoundSpeed = GetFluidModel()->GetSoundSpeed(); + auto AverageSoundSpeed = GetFluidModel()->GetSoundSpeed(); AverageTurboMach[0] = AverageTurboVelocity[val_marker][iSpan][0]/AverageSoundSpeed; AverageTurboMach[1] = AverageTurboVelocity[val_marker][iSpan][1]/AverageSoundSpeed; @@ -6227,167 +6198,32 @@ void CEulerSolver::BC_Giles(CGeometry *geometry, CSolver **solver_container, CNu AverageTurboMach[1] -= geometry->GetAverageTangGridVel(val_marker,iSpan)/AverageSoundSpeed; } - AvgMach = AverageTurboMach[0]*AverageTurboMach[0] + AverageTurboMach[1]*AverageTurboMach[1]; - - kend = geometry->GetnFreqSpan(val_marker, iSpan); - kend_max = geometry->GetnFreqSpanMax(config->GetMarker_All_TurbomachineryFlag(val_marker)); + 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); + su2double Pressure_e; //Useful to pass by reference + switch(config->GetKind_Data_Giles(Marker_Tag)){ case TOTAL_CONDITIONS_PT: - /*--- 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]; - } - } + BC_Giles_Total_Inlet(conv_numerics, config, val_marker, c_avg, R, R_c_inv, R_c, iSpan, iSpan); break; case TOTAL_CONDITIONS_PT_1D: - /*--- 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]; - } - } + BC_Giles_Total_Inlet(conv_numerics, config, val_marker, c_avg, R, R_c_inv, R_c, iSpan, nSpanWiseSections); break; case MIXING_IN: case MIXING_OUT: - /* --- 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); + BC_Giles_Mixing(conv_numerics, val_marker, deltaprim, c_avg, iSpan); break; case MIXING_IN_1D: case MIXING_OUT_1D: - /* --- 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); + BC_Giles_Mixing(conv_numerics, val_marker, deltaprim, c_avg, nSpanWiseSections); break; @@ -6396,13 +6232,8 @@ void CEulerSolver::BC_Giles(CGeometry *geometry, CSolver **solver_container, CNu Pressure_e /= config->GetPressure_Ref(); /* --- Compute avg characteristic jump --- */ - 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); - } + c_avg[nDim+1] = -2.0*(AveragePressure[val_marker][iSpan]-Pressure_e); + break; case STATIC_PRESSURE_1D: @@ -6410,13 +6241,8 @@ void CEulerSolver::BC_Giles(CGeometry *geometry, CSolver **solver_container, CNu Pressure_e /= config->GetPressure_Ref(); /* --- Compute avg characteristic jump --- */ - 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); - } + c_avg[nDim+1] = -2.0*(AveragePressure[val_marker][nSpanWiseSections]-Pressure_e); + break; case RADIAL_EQUILIBRIUM: @@ -6427,79 +6253,68 @@ 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); - - /*--- Compute avg characteristic jump ---*/ - c_avg[nDim + 1] = -2.0*(AveragePressure[val_marker][iSpan]-Pressure_e); - break; - } + auto const iZone = config->GetiZone(); + /*--- Loop over all the vertices on this boundary marker ---*/ SU2_OMP_FOR_DYN(OMP_MIN_SIZE) - for (iVertex = 0; iVertex < geometry->GetnVertexSpan(val_marker,iSpan); iVertex++) { + for (unsigned long iVertex = 0; iVertex < geometry->GetnVertexSpan(val_marker,iSpan); iVertex++) { /*--- using the other vertex information for retrieving some information ---*/ - oldVertex = geometry->turbovertex[val_marker][iSpan][iVertex]->GetOldVertex(); - V_boundary= GetCharacPrimVar(val_marker, oldVertex); + auto oldVertex = geometry->turbovertex[val_marker][iSpan][iVertex]->GetOldVertex(); + auto V_boundary= GetCharacPrimVar(val_marker, oldVertex); /*--- Index of the closest interior node ---*/ - Point_Normal = geometry->vertex[val_marker][oldVertex]->GetNormal_Neighbor(); + auto 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 (iDim = 0; iDim < nDim; iDim++) Normal[iDim] = -Normal[iDim]; + for (unsigned short iDim = 0; iDim < nDim; iDim++) Normal[iDim] = -Normal[iDim]; conv_numerics->SetNormal(Normal); /*--- find the node related to the vertex ---*/ - iPoint = geometry->turbovertex[val_marker][iSpan][iVertex]->GetNode(); + auto 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 ---*/ - V_domain = nodes->GetPrimitive(iPoint); + auto V_domain = nodes->GetPrimitive(iPoint); /*--- Retrieve domain Secondary variables ---*/ - S_domain = nodes->GetSecondary(iPoint); + auto S_domain = nodes->GetSecondary(iPoint); /*--- Compute the internal state u_i ---*/ - Velocity2_i = 0; - for (iDim = 0; iDim < nDim; iDim++) + su2double Velocity2_i = 0; + for (unsigned short iDim = 0; iDim < nDim; iDim++) { Velocity_i[iDim] = nodes->GetVelocity(iPoint,iDim); Velocity2_i += Velocity_i[iDim]*Velocity_i[iDim]; } - Density_i = nodes->GetDensity(iPoint); + auto Density_i = nodes->GetDensity(iPoint); - Energy_i = nodes->GetEnergy(iPoint); - StaticEnergy_i = Energy_i - 0.5*Velocity2_i; + auto Energy_i = nodes->GetEnergy(iPoint); + auto StaticEnergy_i = Energy_i - 0.5*Velocity2_i; GetFluidModel()->SetTDState_rhoe(Density_i, StaticEnergy_i); - Pressure_i = GetFluidModel()->GetPressure(); + auto 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]; } @@ -6509,8 +6324,13 @@ void CEulerSolver::BC_Giles(CGeometry *geometry, CSolver **solver_container, CNu conv_numerics->GetCharJump(AverageSoundSpeed, AverageDensity[val_marker][iSpan], deltaprim, cj); - pitch = geometry->GetMaxAngularCoord(val_marker, iSpan) - geometry->GetMinAngularCoord(val_marker,iSpan); - theta = geometry->turbovertex[val_marker][iSpan][iVertex]->GetRelAngularCoord(); + 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; switch(config->GetKind_Data_Giles(Marker_Tag)) { @@ -6520,11 +6340,18 @@ 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){ - 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++){ + for(unsigned long 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; @@ -6543,16 +6370,10 @@ void CEulerSolver::BC_Giles(CGeometry *geometry, CSolver **solver_container, CNu c2js_Re = c2js.real(); c3js_Re = c3js.real(); - 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]; - } + 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 }else{ if (AverageTurboVelocity[val_marker][iSpan][1] >= 0.0){ @@ -6560,40 +6381,21 @@ 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]); - 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; + dcjs[0] = 0.0 - cj[0]; + dcjs[1] = c2js_Re - cj[1]; + dcjs[2] = 0.0 - cj[2]; + dcjs[nDim] = c3js_Re - cj[nDim]; } + }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]; @@ -6608,10 +6410,11 @@ 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 MASS_FLOW_OUTLET: + case STATIC_PRESSURE:case STATIC_PRESSURE_1D:case MIXING_OUT:case RADIAL_EQUILIBRIUM:case MIXING_OUT_1D: /* --- 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){ @@ -6635,9 +6438,9 @@ void CEulerSolver::BC_Giles(CGeometry *geometry, CSolver **solver_container, CNu }else{ /* --- subsonic Giles implementation ---*/ - Beta_inf= I*complex(sqrt(1.0 - AvgMach)); - cOutjs = complex(0.0,0.0); - for(k=0; k < 2*kend_max+1; k++){ + 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++){ freq = k - kend_max; if(freq >= (long)(-kend) && freq <= (long)(kend) && AverageTurboMach[0] > config->GetAverageMachLimit()){ TwoPiThetaFreq_Pitch = 2*PI_NUMBER*freq*theta/pitch; @@ -6652,41 +6455,25 @@ void CEulerSolver::BC_Giles(CGeometry *geometry, CSolver **solver_container, CNu } cOutjs_Re = cOutjs.real(); - if (nDim == 2){ - dcjs[3] = cOutjs_Re - cj[3]; - } - else{ - dcjs[4] = cOutjs_Re - cj[4]; - } + dcjs[nDim+1] = cOutjs_Re - cj[nDim+1]; } } else{ - if (nDim == 2){ - dcjs[3] = 0.0; - } - else{ - dcjs[4] = 0.0; - } + dcjs[nDim+1] = 0.0; } /* --- Impose Outlet BC Non-Reflecting --- */ delta_c[0] = cj[0]; delta_c[1] = cj[1]; delta_c[2] = cj[2]; - 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]; - } - + delta_c[3] = cj[3]; + delta_c[nDim+1] = relfacAvg*c_avg[nDim+1] + relfacFou*dcjs[nDim+1]; /*--- 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[2] = 0.0; + delta_c[3] = 0.0; if (nDim == 3)delta_c[4] = 0.0; } @@ -6698,31 +6485,25 @@ void CEulerSolver::BC_Giles(CGeometry *geometry, CSolver **solver_container, CNu } /*--- Compute primitive jump from characteristic variables ---*/ - for (iVar = 0; iVar < nVar; iVar++) + for (unsigned short iVar = 0; iVar < nVar; iVar++) { deltaprim[iVar]=0.0; - for (jVar = 0; jVar < nVar; jVar++) + for (unsigned short jVar = 0; jVar < nVar; jVar++) { deltaprim[iVar] += R_Matrix[iVar][jVar]*delta_c[jVar]; } } /*--- retrieve boundary variables ---*/ - Density_b = AverageDensity[val_marker][iSpan] + deltaprim[0]; + su2double 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]; - 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]; - } - + turboVelocity[nDim-1] = AverageTurboVelocity[val_marker][iSpan][nDim-1] + deltaprim[nDim]; + su2double Pressure_b = AveragePressure[val_marker][iSpan] + deltaprim[nDim+1]; ComputeBackVelocity(turboVelocity, turboNormal, Velocity_b, config->GetMarker_All_TurbomachineryFlag(val_marker), config->GetKind_TurboMachinery(iZone)); - Velocity2_b = 0.0; - for (iDim = 0; iDim < nDim; iDim++) { + su2double Velocity2_b = 0.0; + for (unsigned short iDim = 0; iDim < nDim; iDim++) { Velocity2_b+= Velocity_b[iDim]*Velocity_b[iDim]; } @@ -6730,20 +6511,20 @@ void CEulerSolver::BC_Giles(CGeometry *geometry, CSolver **solver_container, CNu Pressure_b = Pressure_i; Density_b = Density_i; Velocity2_b = 0.0; - for (iDim = 0; iDim < nDim; iDim++) { + for (unsigned short 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); - Energy_b = GetFluidModel()->GetStaticEnergy() + 0.5*Velocity2_b; - Temperature_b= GetFluidModel()->GetTemperature(); - Enthalpy_b = Energy_b + 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; /*--- Primitive variables, using the derived quantities ---*/ V_boundary[0] = Temperature_b; - for (iDim = 0; iDim < nDim; iDim++) + for (unsigned short iDim = 0; iDim < nDim; iDim++) V_boundary[iDim+1] = Velocity_b[iDim]; V_boundary[nDim+1] = Pressure_b; V_boundary[nDim+2] = Density_b; @@ -6771,12 +6552,12 @@ void CEulerSolver::BC_Giles(CGeometry *geometry, CSolver **solver_container, CNu LinSysRes.AddBlock(iPoint, residual); /*--- Jacobian contribution for implicit integration ---*/ - if (implicit) + if (config->GetKind_TimeIntScheme() == EULER_IMPLICIT) Jacobian.AddBlock2Diag(iPoint, residual.jacobian_i); /*--- Viscous contribution ---*/ - if (viscous) { + if (config->GetViscous()) { /*--- Set laminar and eddy viscosity at the infinity ---*/ @@ -6830,7 +6611,7 @@ void CEulerSolver::BC_Giles(CGeometry *geometry, CSolver **solver_container, CNu /*--- Jacobian contribution for implicit integration ---*/ - if (implicit) + if (config->GetKind_TimeIntScheme() == EULER_IMPLICIT) Jacobian.SubtractBlock2Diag(iPoint, residual.jacobian_i); } @@ -6852,11 +6633,11 @@ void CEulerSolver::BC_Giles(CGeometry *geometry, CSolver **solver_container, CNu delete [] delta_c; delete [] deltaprim; delete [] cj; - for (iVar = 0; iVar < nVar; iVar++) + for (unsigned short iVar = 0; iVar < nVar; iVar++) { delete [] R_Matrix[iVar]; } - for (iVar = 0; iVar < nVar-1; iVar++) + for (unsigned short iVar = 0; iVar < nVar-1; iVar++) { delete [] R_c_inv[iVar]; delete [] R_c[iVar]; @@ -6874,6 +6655,73 @@ 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) { @@ -9099,7 +8947,8 @@ 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->GetMarker_All_Giles(iMarker) || config->GetBoolRiemann()){ // May have to implement something similar for Riemann? + if(config->GetBoolGiles() || config->GetBoolRiemann()){ 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 215445d2c5e..5b81394bd9a 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(Pressure_Thermodynamic, config); + auxFluidModel = new CFluidScalar(config->GetSpecific_Heat_Cp(), config->GetGas_Constant(), 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(Pressure_ThermodynamicND, config); + fluidModel = new CFluidScalar(Specific_Heat_CpND, Gas_ConstantND, 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 f8327c868bd..bfbd4edad58 100644 --- a/SU2_CFD/src/solvers/CSolver.cpp +++ b/SU2_CFD/src/solvers/CSolver.cpp @@ -1730,7 +1730,6 @@ 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. */ @@ -1745,12 +1744,10 @@ 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/species transport. */ - const su2double linRes = max(solverFlow->GetResLinSolver(), max(linResTurb, linResSpecies)); + /* Max linear residual between flow and turbulence. */ + const su2double linRes = max(solverFlow->GetResLinSolver(), linResTurb); /* Tolerance limited to an acceptable value. */ const su2double linTol = max(acceptableLinTol, config->GetLinear_Solver_Error()); @@ -1782,11 +1779,6 @@ 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 @@ -1830,7 +1822,6 @@ 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) { @@ -1897,9 +1888,6 @@ 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 af200ef7603..c893649c1dd 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 = [-10.467026, -2.871699, -19.214627, -13.508254, -11.582396, -6.306163, 73273, 73273, 0.019884, 82.491] + 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] test_list.append(Jones_tc_restart) # 2D axial stage @@ -582,15 +582,6 @@ 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 4047b6b39ba..0137b08806c 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.829186, -8.875103, -9.609509, -8.075211, -7.759491, -4.360714] + Aachen_3D_restart.test_vals = [-9.853215, -8.891480, -9.610418, -8.028579, -7.792957, -4.384186] 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 = [-10.467612, -2.871708, -19.345651, -13.625871, -11.582397, -6.306168, 73273, 73273, 0.019884, 82.491] + Jones_tc_restart.test_vals = [-6.614623, -3.001323, -14.336147, -8.776081, -11.382919, -5.852327, 73273, 73273, 0.019884, 82.491] test_list.append(Jones_tc_restart) # 2D axial stage @@ -1112,15 +1112,6 @@ 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 ### ###################################### diff --git a/TestCases/serial_regression.py b/TestCases/serial_regression.py index 08b45e3b3ff..080b41053b7 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.829185, -8.875103, -9.609505, -8.075194, -7.759490, -4.360713] + Aachen_3D_restart.test_vals = [-9.853207, -8.891479, -9.610411, -8.028566, -7.792947, -4.384177] 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 = [-10.466644, -2.871703, -19.193464, -13.487820, -11.582397, -6.306167, 73273, 73273, 0.019884, 82.491] + 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] test_list.append(Jones_tc_restart) # 2D axial stage @@ -889,15 +889,6 @@ 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 deleted file mode 100644 index a52301e7153..00000000000 --- a/TestCases/turbomachinery/multi_interface/multi_interface_rst.cfg +++ /dev/null @@ -1,151 +0,0 @@ -%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% -% % -% 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 deleted file mode 100644 index 3d7a019754a..00000000000 --- a/TestCases/turbomachinery/multi_interface/zone_1.cfg +++ /dev/null @@ -1,6 +0,0 @@ -% ----------------------- 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 deleted file mode 100644 index 3d7a019754a..00000000000 --- a/TestCases/turbomachinery/multi_interface/zone_2.cfg +++ /dev/null @@ -1,6 +0,0 @@ -% ----------------------- 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 deleted file mode 100644 index 3d7a019754a..00000000000 --- a/TestCases/turbomachinery/multi_interface/zone_3.cfg +++ /dev/null @@ -1,6 +0,0 @@ -% ----------------------- DYNAMIC MESH DEFINITION -----------------------------% -% -% -% Type of dynamic mesh (NONE, ROTATING_FRAME) -GRID_MOVEMENT= NONE -% diff --git a/config_template.cfg b/config_template.cfg index 4044a611957..00334732b30 100644 --- a/config_template.cfg +++ b/config_template.cfg @@ -672,9 +672,6 @@ 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 % From a284b6b85aa28409d3ed23bf5fb83ca72c0620f8 Mon Sep 17 00:00:00 2001 From: Josh Kelly Date: Thu, 9 Jan 2025 14:50:04 +0000 Subject: [PATCH 22/22] Resolved tag tape debugging errors --- SU2_CFD/src/iteration/CDiscAdjFluidIteration.cpp | 6 +----- 1 file changed, 1 insertion(+), 5 deletions(-) diff --git a/SU2_CFD/src/iteration/CDiscAdjFluidIteration.cpp b/SU2_CFD/src/iteration/CDiscAdjFluidIteration.cpp index 78b227c3ead..7bf63d4a735 100644 --- a/SU2_CFD/src/iteration/CDiscAdjFluidIteration.cpp +++ b/SU2_CFD/src/iteration/CDiscAdjFluidIteration.cpp @@ -472,11 +472,7 @@ void CDiscAdjFluidIteration::SetDependencies(CSolver***** solver, CGeometry**** solvers0[FLOW_SOL]->CompleteComms(geometry0, config[iZone], MPI_QUANTITIES::SOLUTION); if (config[iZone]->GetBoolTurbomachinery()) { - solvers0[FLOW_SOL]->PreprocessAverage(solvers0, geometry0, config[iZone], INFLOW); - solvers0[FLOW_SOL]->PreprocessAverage(solvers0, geometry0, config[iZone], OUTFLOW); - solvers0[FLOW_SOL]->TurboAverageProcess(solvers0, geometry0, config[iZone], INFLOW); - solvers0[FLOW_SOL]->TurboAverageProcess(solvers0, geometry0, config[iZone], OUTFLOW); - solvers0[FLOW_SOL]->GatherInOutAverageValues(config[iZone], geometry0); + solvers0[FLOW_SOL]->InitTurboContainers(geometry0, config[iZone]); } if (turbulent && !config[iZone]->GetFrozen_Visc_Disc()) { solvers0[TURB_SOL]->Postprocessing(geometry0, solvers0,