From 3801e4db3ed558647c6227180e177771fc0202b3 Mon Sep 17 00:00:00 2001 From: Ole Burghardt Date: Wed, 14 Aug 2024 18:17:14 +0200 Subject: [PATCH 01/16] 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/16] 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/16] 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/16] 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/16] 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/16] 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/16] 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 6e65dd17bae693cb21c2b893f2123fcd9d25e722 Mon Sep 17 00:00:00 2001 From: Ole Burghardt Date: Tue, 24 Sep 2024 17:26:02 +0200 Subject: [PATCH 08/16] 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 09/16] 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 10/16] 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 11/16] 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 12/16] 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 13/16] 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 14/16] 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 15/16] 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 16/16] 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