diff --git a/.github/workflows/regression.yml b/.github/workflows/regression.yml index cef9d68216ed..73533f6ce341 100644 --- a/.github/workflows/regression.yml +++ b/.github/workflows/regression.yml @@ -23,7 +23,7 @@ jobs: strategy: fail-fast: false matrix: - config_set: [BaseMPI, ReverseMPI, ForwardMPI, BaseNoMPI, ReverseNoMPI, ForwardNoMPI, BaseOMP, ReverseOMP, ForwardOMP] + config_set: [BaseMPI, ReverseMPI, ForwardMPI, BaseNoMPI, ReverseNoMPI, ForwardNoMPI, ReverseTagNoMPI, BaseOMP, ReverseOMP, ForwardOMP] include: - config_set: BaseMPI flags: '-Denable-pywrapper=true -Denable-coolprop=true -Denable-mpp=true -Dinstall-mpp=true -Denable-mlpcpp=true -Denable-tests=true --warnlevel=2' @@ -37,6 +37,8 @@ jobs: flags: '-Denable-autodiff=true -Denable-normal=false -Dwith-mpi=disabled -Denable-pywrapper=true -Denable-tests=true --warnlevel=3 --werror' - config_set: ForwardNoMPI flags: '-Denable-directdiff=true -Denable-normal=false -Dwith-mpi=disabled -Denable-tests=true --warnlevel=3 --werror' + - config_set: ReverseTagNoMPI + flags: '-Denable-autodiff=true -Denable-normal=false -Dwith-mpi=disabled -Denable-pywrapper=true -Denable-tests=true --warnlevel=3 --werror -Dcodi-tape=Tag' - config_set: BaseOMP flags: '-Dwith-omp=true -Denable-mixedprec=true -Denable-pywrapper=true -Denable-tecio=false --warnlevel=3 --werror' - config_set: ReverseOMP @@ -215,6 +217,55 @@ jobs: with: entrypoint: /bin/rm args: -rf install install_bin.tgz src ccache ${{ matrix.config_set }} + discadj_tape_tests: + if: inputs.runner != 'ARM64' + runs-on: ${{ inputs.runner || 'ubuntu-latest' }} + name: Tape Tests + needs: build + strategy: + fail-fast: false + matrix: + testscript: ['serial_regression_AD.py'] + include: + - testscript: 'serial_regression_AD.py' + tag: TagNoMPI + steps: + - name: Pre Cleanup + uses: docker://ghcr.io/su2code/su2/test-su2:240320-1536 + with: + entrypoint: /bin/rm + args: -rf install install_bin.tgz src ccache ${{ matrix.config_set }} + - name: Download All artifacts + uses: actions/download-artifact@v4 + - name: Uncompress and Move Binaries + run: | + BIN_FOLDER="$PWD/install/bin" + mkdir -p $BIN_FOLDER + ls -lah $BIN_FOLDER + for type in Base Reverse Forward; do + TYPE_FOLDER="${type}${{matrix.tag}}" + echo "Processing '$TYPE_FOLDER' ..." + if [ -d $TYPE_FOLDER ]; then + pushd $TYPE_FOLDER + ls -lah + tar -zxvf install_bin.tgz + ls -lah install/bin/ + cp -r install/* $BIN_FOLDER/../ + popd; + fi + done + chmod a+x $BIN_FOLDER/* + ls -lahR $BIN_FOLDER + - name: Run Tests in Container + uses: docker://ghcr.io/su2code/su2/test-su2:240320-1536 + with: + # -t -c + args: -b ${{github.ref}} -t develop -c develop -s ${{matrix.testscript}} -a "--tapetests" + - name: Cleanup + uses: docker://ghcr.io/su2code/su2/test-su2:240320-1536 + with: + entrypoint: /bin/rm + args: -rf install install_bin.tgz src ccache ${{ matrix.config_set }} thread_sanitizer_tests: if: inputs.runner != 'ARM64' runs-on: ${{ inputs.runner || 'ubuntu-latest' }} diff --git a/Common/include/CConfig.hpp b/Common/include/CConfig.hpp index 1d98b92534e7..66ceff88955b 100644 --- a/Common/include/CConfig.hpp +++ b/Common/include/CConfig.hpp @@ -591,14 +591,17 @@ class CConfig { MUSCL_AdjTurb; /*!< \brief MUSCL scheme for the adj turbulence equations.*/ bool MUSCL_Species; /*!< \brief MUSCL scheme for the species equations.*/ bool Use_Accurate_Jacobians; /*!< \brief Use numerically computed Jacobians for AUSM+up(2) and SLAU(2). */ + bool Use_Accurate_Turb_Jacobians; /*!< \brief Use numerically computed Jacobians for standard SA turbulence model. */ bool EulerPersson; /*!< \brief Boolean to determine whether this is an Euler simulation with Persson shock capturing. */ bool FSI_Problem = false,/*!< \brief Boolean to determine whether the simulation is FSI or not. */ Multizone_Problem; /*!< \brief Boolean to determine whether we are solving a multizone problem. */ //bool ContactResistance = false; /*!< \brief Apply contact resistance for conjugate heat transfer. */ unsigned short nID_DV; /*!< \brief ID for the region of FEM when computed using direct differentiation. */ - bool AD_Mode; /*!< \brief Algorithmic Differentiation support. */ - bool AD_Preaccumulation; /*!< \brief Enable or disable preaccumulation in the AD mode. */ + bool AD_Mode; /*!< \brief Algorithmic Differentiation support. */ + bool AD_Preaccumulation; /*!< \brief Enable or disable preaccumulation in the AD mode. */ + CHECK_TAPE_TYPE AD_CheckTapeType; /*!< \brief Type of tape that is checked in a tape debug run. */ + CHECK_TAPE_VARIABLES AD_CheckTapeVariables; /*!< \brief Type of variables that are checked in a tape debug run. */ STRUCT_COMPRESS Kind_Material_Compress; /*!< \brief Determines if the material is compressible or incompressible (structural analysis). */ STRUCT_MODEL Kind_Material; /*!< \brief Determines the material model to be used (structural analysis). */ STRUCT_DEFORMATION Kind_Struct_Solver; /*!< \brief Determines the geometric condition (small or large deformations) for structural analysis. */ @@ -636,6 +639,9 @@ class CConfig { su2double Relaxation_Factor_Adjoint; /*!< \brief Relaxation coefficient for variable updates of adjoint solvers. */ su2double Relaxation_Factor_CHT; /*!< \brief Relaxation coefficient for the update of conjugate heat variables. */ su2double EntropyFix_Coeff; /*!< \brief Entropy fix coefficient. */ + su2double Max_Update_SST; /*!< \brief Cap for the Under-Relaxation Factor for SST Turbulent Variables*/ + su2double Max_Update_SA; /*!< \brief Cap for the Under-Relaxation Factor for SA Turbulent Variables*/ + su2double Max_Update_Flow; /*!< \brief Cap for the Under-Relaxation Factor for Flow Density and Energy Variables*/ unsigned short nLocationStations, /*!< \brief Number of section cuts to make when outputting mesh and cp . */ nWingStations; /*!< \brief Number of section cuts to make when calculating internal volume. */ su2double Kappa_1st_AdjFlow, /*!< \brief Lax 1st order dissipation coefficient for adjoint flow equations (coarse multigrid levels). */ @@ -822,6 +828,7 @@ class CConfig { SurfSens_FileName, /*!< \brief Output file for the sensitivity on the surface (discrete adjoint). */ VolSens_FileName, /*!< \brief Output file for the sensitivity in the volume (discrete adjoint). */ ObjFunc_Hess_FileName; /*!< \brief Hessian approximation obtained by the Sobolev smoothing solver. */ + bool Multizone_Adapt_FileName; /*!< \brief Append zone number to solution and restart file names. */ bool Wrt_Performance, /*!< \brief Write the performance summary at the end of a calculation. */ @@ -1055,7 +1062,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. */ @@ -1897,6 +1905,12 @@ class CConfig { */ su2double GetPressure_FreeStreamND(void) const { return Pressure_FreeStreamND; } + /*! + * \brief Get a reference to the non-dimensionalized freestream pressure (used for AD tracking). + * \return Reference to non-dimensionalized freestream pressure. + */ + su2double& GetPressure_FreeStreamND(void) { return Pressure_FreeStreamND; } + /*! * \brief Get the value of the thermodynamic pressure. * \return Thermodynamic pressure. @@ -1922,6 +1936,12 @@ class CConfig { */ su2double GetTemperature_FreeStreamND(void) const { return Temperature_FreeStreamND; } + /*! + * \brief Get a reference to the non-dimensionalized freestream temperature (used for AD tracking). + * \return Reference to non-dimensionalized freestream temperature. + */ + su2double& GetTemperature_FreeStreamND(void) { return Temperature_FreeStreamND; } + /*! * \brief Get the value of the non-dimensionalized vibrational-electronic freestream temperature. * \return Non-dimensionalized vibrational-electronic freestream temperature. @@ -4208,6 +4228,24 @@ class CConfig { */ su2double GetRelaxation_Factor_CHT(void) const { return Relaxation_Factor_CHT; } + /*! + * \brief Get the under-relaxation for flow variables density and energy. + * \return under-relaxation for flow variables. + */ + su2double GetUnderRelax_Flow(void) const { return Max_Update_Flow; } + + /*! + * \brief Get the under-relaxation for SA variable, nu_tilde. + * \return under-relaxation for SA variables. + */ + su2double GetUnderRelax_SA(void) const { return Max_Update_SA; } + + /*! + * \brief Get the under-relaxation for SST turbulence variables k and omega. + * \return under-relaxation for SST variables. + */ + su2double GetUnderRelax_SST(void) const { return Max_Update_SST; } + /*! * \brief Get the number of samples used in quasi-Newton methods. */ @@ -4500,6 +4538,12 @@ class CConfig { */ bool GetUse_Accurate_Jacobians(void) const { return Use_Accurate_Jacobians; } + /*! + * \brief Get whether to "Use Accurate Jacobians" for Standard SA turbulence model. + * \return yes/no. + */ + bool GetUse_Accurate_Turb_Jacobians(void) const { return Use_Accurate_Turb_Jacobians; } + /*! * \brief Get the kind of integration scheme (explicit or implicit) * for the flow equations. @@ -5384,19 +5428,25 @@ class CConfig { bool GetWrt_Volume_Overwrite(void) const { return Wrt_Volume_Overwrite; } /*! - * \brief Provides the number of varaibles. + * \brief Get whether filenames are appended the zone number automatically (multiphysics solver). + * \return Flag for appending zone numbers to restart and solution filenames. If Flag=true, zone numer is appended. + */ + bool GetMultizone_AdaptFilename(void) const { return Multizone_Adapt_FileName; } + + /*! + * \brief Provides the number of variables. * \return Number of variables. */ unsigned short GetnVar(void); /*! - * \brief Provides the number of varaibles. + * \brief Provides the number of variables. * \return Number of variables. */ unsigned short GetnZone(void) const { return nZone; } /*! - * \brief Provides the number of varaibles. + * \brief Provides the number of variables. * \return Number of variables. */ unsigned short GetiZone(void) const { return iZone; } @@ -8779,6 +8829,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. @@ -9216,6 +9272,16 @@ class CConfig { */ su2double GetConst_DES(void) const { return Const_DES; } + /*! + * \brief Get the type of tape that will be checked in a tape debug run. + */ + CHECK_TAPE_TYPE GetAD_CheckTapeType(void) const { return AD_CheckTapeType; } + + /*! + * \brief Get the type of variables that will be checked for in a tape debug run. + */ + CHECK_TAPE_VARIABLES GetAD_CheckTapeVariables(void) const { return AD_CheckTapeVariables; } + /*! * \brief Get if AD preaccumulation should be performed. */ @@ -9605,6 +9671,11 @@ class CConfig { */ unsigned short GetnVolumeOutputFiles() const { return nVolumeOutputFiles; } + /*! + * \brief GetnVolumeOutputFrequencies + */ + unsigned short GetnVolumeOutputFrequencies() const { return nVolumeOutputFrequencies; } + /*! * \brief GetVolumeOutputFrequency * \param[in] iFile: index of file number for which the writing frequency needs to be returned. diff --git a/Common/include/basic_types/ad_structure.hpp b/Common/include/basic_types/ad_structure.hpp index 0bd979b9663a..5a02b23c2c67 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, @@ -45,6 +48,17 @@ namespace AD { */ inline void StartRecording() {} +/*! + * \brief Pause the recording of the operations and involved variables. + * If called, all operations occuring after the call will not be stored on the computational graph. + */ +inline void PauseRecording() {} + +/*! + * \brief Resume the recording of the operations and variables after the recording had been paused. + */ +inline void ResumeRecording() {} + /*! * \brief Stops the recording of the operations and variables. */ @@ -101,14 +115,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 +279,50 @@ 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 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 Sets the tag of a variable to 0. + * \param[in] v - the variable whose tag is cleared. + */ +inline void ClearTagOnVariable(su2double& v) {} + +/*! + * \brief Struct to store information about errors during a tag debug run. + */ +struct ErrorReport {}; + +/*! + * \brief Set a reference to the output file of an ErrorReport. + * \param[in] report - the ErrorReport whose output file is set. + * \param[in] output_file - pointer to the output file. + */ +inline void SetDebugReportFile(ErrorReport& report, std::ostream* output_file) {} + +/*! + * \brief Set the ErrorReport to which error information from a tag debug recording is written. + * \param[in] report - the ErrorReport to which error information is written. + */ +inline void SetTagErrorCallback(ErrorReport& report) {} + +/*! + * \brief Reset the error counter in an ErrorReport. + * \param[in] report - the ErrorReport whose error counter is resetted. + */ +inline void ResetErrorCounter(ErrorReport& report) {} + +/*! + * \brief Get the error count of an ErrorReport. + * \param[in] report - the ErrorReport whose pointer to its error counter is returned. + * \return Value of the error counter. + */ +inline unsigned long GetErrorCount(const ErrorReport& report) { return 0; } /*! * \brief Pushes back the current tape position to the tape position's vector. @@ -304,6 +367,7 @@ inline void EndNoSharedReading() {} using CheckpointHandler = codi::ExternalFunctionUserData; using Tape = su2double::Tape; +using Identifier = su2double::Identifier; #ifdef HAVE_OPDI using ExtFuncHelper = codi::OpenMPExternalFunctionHelper; @@ -349,6 +413,10 @@ FORCEINLINE void ResetInput(su2double& data) { data = data.getValue(); } FORCEINLINE void StartRecording() { AD::getTape().setActive(); } +FORCEINLINE void PauseRecording() { AD::getTape().setPassive(); } + +FORCEINLINE void ResumeRecording() { AD::getTape().setActive(); } + FORCEINLINE void StopRecording() { AD::getTape().setPassive(); } FORCEINLINE bool TapeActive() { return AD::getTape().isActive(); } @@ -470,15 +538,15 @@ 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. - return; +FORCEINLINE void SetDerivative(Identifier index, const double val) { + // Allow multiple threads to "set the derivative" of passive variables without causing data races. + if (!AD::getTape().isIdentifierActive(index)) return; AD::getTape().setGradient(index, val, codi::AdjointsManagement::Manual); } @@ -488,13 +556,11 @@ 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 bool IsIdentifierActive(su2double const& value) { - return getTape().isIdentifierActive(value.getIdentifier()); -} +FORCEINLINE Identifier GetPassiveIndex() { return AD::getTape().getPassiveIndex(); } /*--- Base case for parameter pack expansion. ---*/ FORCEINLINE void SetPreaccIn() {} @@ -502,7 +568,7 @@ FORCEINLINE void SetPreaccIn() {} template ::value> = 0> FORCEINLINE void SetPreaccIn(const T& data, Ts&&... moreData) { if (!PreaccActive) return; - if (IsIdentifierActive(data)) PreaccHelper.addInput(data); + PreaccHelper.addInput(data); SetPreaccIn(moreData...); } @@ -515,9 +581,7 @@ template FORCEINLINE void SetPreaccIn(const T& data, const int size) { if (PreaccActive) { for (int i = 0; i < size; i++) { - if (IsIdentifierActive(data[i])) { - PreaccHelper.addInput(data[i]); - } + PreaccHelper.addInput(data[i]); } } } @@ -527,9 +591,7 @@ 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])) { - PreaccHelper.addInput(data[i][j]); - } + PreaccHelper.addInput(data[i][j]); } } } @@ -547,7 +609,7 @@ FORCEINLINE void SetPreaccOut() {} template ::value> = 0> FORCEINLINE void SetPreaccOut(T& data, Ts&&... moreData) { if (!PreaccActive) return; - if (IsIdentifierActive(data)) PreaccHelper.addOutput(data); + PreaccHelper.addOutput(data); SetPreaccOut(moreData...); } @@ -555,9 +617,7 @@ template FORCEINLINE void SetPreaccOut(T&& data, const int size) { if (PreaccActive) { for (int i = 0; i < size; i++) { - if (IsIdentifierActive(data[i])) { - PreaccHelper.addOutput(data[i]); - } + PreaccHelper.addOutput(data[i]); } } } @@ -567,9 +627,7 @@ 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])) { - PreaccHelper.addOutput(data[i][j]); - } + PreaccHelper.addOutput(data[i][j]); } } } @@ -676,6 +734,40 @@ FORCEINLINE void ResumePreaccumulation(bool wasActive) { SU2_OMP_SAFE_GLOBAL_ACCESS(PreaccEnabled = true;) } +struct ErrorReport { + unsigned long ErrorCounter = 0; + std::ostream* out = &std::cout; +}; + +FORCEINLINE void ResetErrorCounter(ErrorReport& report) { report.ErrorCounter = 0; } + +FORCEINLINE void SetDebugReportFile(ErrorReport& report, std::ostream* output_file) { report.out = output_file; } + +FORCEINLINE unsigned long GetErrorCount(const ErrorReport& report) { return report.ErrorCounter; } + +#ifdef CODI_TAG_TAPE + +FORCEINLINE void SetTag(int tag) { AD::getTape().setCurTag(tag); } +FORCEINLINE void ClearTagOnVariable(su2double& v) { AD::getTape().clearTagOnVariable(v); } + +static void tagErrorCallback(const int& correctTag, const int& wrongTag, void* userData) { + auto* report = static_cast(userData); + + report->ErrorCounter += 1; + *(report->out) << "Use of variable with bad tag '" << wrongTag << "', should be '" << correctTag << "'." << std::endl; +} + +FORCEINLINE void SetTagErrorCallback(ErrorReport& report) { + AD::getTape().setTagErrorCallback(tagErrorCallback, &report); +} + +#else +FORCEINLINE void SetTag(int tag) {} +FORCEINLINE void ClearTagOnVariable(su2double& v) {} +FORCEINLINE void SetTagErrorCallback(ErrorReport report) {} + +#endif // CODI_TAG_TAPE + #endif // CODI_REVERSE_TYPE void Initialize(); diff --git a/Common/include/basic_types/datatype_structure.hpp b/Common/include/basic_types/datatype_structure.hpp index e16c223e8a69..9d5eb4559b02 100644 --- a/Common/include/basic_types/datatype_structure.hpp +++ b/Common/include/basic_types/datatype_structure.hpp @@ -89,8 +89,15 @@ void SetDerivative(su2double& data, const passivedouble& val); FORCEINLINE void SetValue(su2double& data, const passivedouble& val) { data.setValue(val); } +FORCEINLINE passivedouble GetValue(const double& data) { return data; } + FORCEINLINE passivedouble GetValue(const su2double& data) { return data.getValue(); } +template +FORCEINLINE passivedouble GetValue(const codi::ExpressionInterface& data) { + return data.cast().getValue(); +} + FORCEINLINE void SetSecondary(su2double& data, const passivedouble& val) { data.setGradient(val); } FORCEINLINE void SetDerivative(su2double& data, const passivedouble& val) { data.setGradient(val); } diff --git a/Common/include/code_config.hpp b/Common/include/code_config.hpp index 48625f20bf44..5509dd8d525f 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/Common/include/geometry/dual_grid/CPoint.hpp b/Common/include/geometry/dual_grid/CPoint.hpp index f3912aa2a2de..05d64e25482a 100644 --- a/Common/include/geometry/dual_grid/CPoint.hpp +++ b/Common/include/geometry/dual_grid/CPoint.hpp @@ -113,9 +113,10 @@ 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 before solver iteration. */ + su2matrix + AD_OutputIndex; /*!< \brief Indices of Coord variables in the adjoint vector after solver iteration. */ /*! * \brief Allocate fields required by the minimal constructor. diff --git a/Common/include/linear_algebra/CSysMatrix.hpp b/Common/include/linear_algebra/CSysMatrix.hpp index 7d0e87b92e35..d3212e7a9f5f 100644 --- a/Common/include/linear_algebra/CSysMatrix.hpp +++ b/Common/include/linear_algebra/CSysMatrix.hpp @@ -553,6 +553,22 @@ class CSysMatrix { SetBlock(block_i, block_j, val_block, alpha); } + /*! + * \brief Returns the 4 blocks ii, ij, ji, jj used by "UpdateBlocks". + * \note This method assumes an FVM-type sparse pattern. + * \param[in] edge - Index of edge that connects iPoint and jPoint. + * \param[in] iPoint - Row to which we add the blocks. + * \param[in] jPoint - Row from which we subtract the blocks. + * \param[out] bii, bij, bji, bjj - Blocks of the matrix. + */ + inline void GetBlocks(unsigned long iEdge, unsigned long iPoint, unsigned long jPoint, ScalarType*& bii, + ScalarType*& bij, ScalarType*& bji, ScalarType*& bjj) { + bii = &matrix[dia_ptr[iPoint] * nVar * nEqn]; + bjj = &matrix[dia_ptr[jPoint] * nVar * nEqn]; + bij = &matrix[edge_ptr(iEdge, 0) * nVar * nEqn]; + bji = &matrix[edge_ptr(iEdge, 1) * nVar * nEqn]; + } + /*! * \brief Update 4 blocks ii, ij, ji, jj (add to i* sub from j*). * \note This method assumes an FVM-type sparse pattern. @@ -566,10 +582,8 @@ class CSysMatrix { template inline void UpdateBlocks(unsigned long iEdge, unsigned long iPoint, unsigned long jPoint, const MatrixType& block_i, const MatrixType& block_j, OtherType scale = 1) { - ScalarType* bii = &matrix[dia_ptr[iPoint] * nVar * nEqn]; - ScalarType* bjj = &matrix[dia_ptr[jPoint] * nVar * nEqn]; - ScalarType* bij = &matrix[edge_ptr(iEdge, 0) * nVar * nEqn]; - ScalarType* bji = &matrix[edge_ptr(iEdge, 1) * nVar * nEqn]; + ScalarType *bii, *bij, *bji, *bjj; + GetBlocks(iEdge, iPoint, jPoint, bii, bij, bji, bjj); unsigned long iVar, jVar, offset = 0; diff --git a/Common/include/option_structure.hpp b/Common/include/option_structure.hpp index 223a53da742b..6f118ce29197 100644 --- a/Common/include/option_structure.hpp +++ b/Common/include/option_structure.hpp @@ -2543,6 +2543,29 @@ static const MapType DirectDiff_Var_Map = { MakePair("ELECTRIC_FIELD", D_EFIELD) }; +/*! + * \brief Types of tapes that can be checked in a tape debug run. + */ +enum class CHECK_TAPE_TYPE { + OBJECTIVE_FUNCTION, /*!< \brief Tape that only includes dependencies and objective function calculation. */ + FULL_SOLVER /*!< \brief Tape that includes dependencies and all solvers. */ +}; +static const MapType CheckTapeType_Map = { + MakePair("OBJECTIVE_FUNCTION_TAPE", CHECK_TAPE_TYPE::OBJECTIVE_FUNCTION) + MakePair("FULL_SOLVER_TAPE", CHECK_TAPE_TYPE::FULL_SOLVER) +}; + +/*! + * \brief Types of variables that can be checked for in a tape debug run. + */ +enum class CHECK_TAPE_VARIABLES { + SOLVER_VARIABLES, /*!< \brief A (debug) tag will be assigned to solver/conservative variables. */ + MESH_COORDINATES /*!< \brief A (debug) tag will be assigned to solver/conservative variables and mesh coordinates. */ +}; +static const MapType CheckTapeVariables_Map = { + MakePair("SOLVER_VARIABLES", CHECK_TAPE_VARIABLES::SOLVER_VARIABLES) + MakePair("SOLVER_VARIABLES_AND_MESH_COORDINATES", CHECK_TAPE_VARIABLES::MESH_COORDINATES) +}; enum class RECORDING { CLEAR_INDICES, @@ -2550,6 +2573,10 @@ enum class RECORDING { MESH_COORDS, MESH_DEFORM, SOLUTION_AND_MESH, + TAG_INIT_SOLVER_VARIABLES, + TAG_CHECK_SOLVER_VARIABLES, + TAG_INIT_SOLVER_AND_MESH, + TAG_CHECK_SOLVER_AND_MESH }; /*! diff --git a/Common/include/parallelization/vectorization.hpp b/Common/include/parallelization/vectorization.hpp index f08017a9bee9..1e4b8b77f327 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 * alignof(Scalar) }; static constexpr bool StoreAsRef = true; private: - alignas(Size * sizeof(Scalar)) Scalar x_[N]; + alignas(Size * alignof(Scalar)) Scalar x_[N]; public: #define ARRAY_BOILERPLATE \ diff --git a/Common/src/CConfig.cpp b/Common/src/CConfig.cpp index 9222ed7f8928..596001ad0b09 100644 --- a/Common/src/CConfig.cpp +++ b/Common/src/CConfig.cpp @@ -1108,10 +1108,16 @@ 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) + DiscreteAdjointDebug = true; +# else + DiscreteAdjointDebug = false; +# endif #else const bool discAdjDefault = false; + DiscreteAdjointDebug = false; #endif /*!\brief MATH_PROBLEM \n DESCRIPTION: Mathematical problem \n Options: DIRECT, ADJOINT \ingroup Config*/ addMathProblemOption("MATH_PROBLEM", ContinuousAdjoint, false, DiscreteAdjoint, discAdjDefault, Restart_Flow, discAdjDefault); @@ -1186,6 +1192,8 @@ void CConfig::SetConfig_Options() { addBoolOption("WRT_VOLUME_OVERWRITE", Wrt_Volume_Overwrite, true); /*!\brief SYSTEM_MEASUREMENTS \n DESCRIPTION: System of measurements \n OPTIONS: see \link Measurements_Map \endlink \n DEFAULT: SI \ingroup Config*/ addEnumOption("SYSTEM_MEASUREMENTS", SystemMeasurements, Measurements_Map, SI); + /*!\brief MULTIZONE_ADAPT_FILENAME \n DESCRIPTION: Append zone number to restart and solution filenames. \ingroup Config*/ + addBoolOption("MULTIZONE_ADAPT_FILENAME", Multizone_Adapt_FileName, YES); /*!\par CONFIG_CATEGORY: FluidModel \ingroup Config*/ /*!\brief FLUID_MODEL \n DESCRIPTION: Fluid model \n OPTIONS: See \link FluidModel_Map \endlink \n DEFAULT: STANDARD_AIR \ingroup Config*/ @@ -1889,6 +1897,12 @@ void CConfig::SetConfig_Options() { addEnumOption("DISCADJ_LIN_PREC", Kind_DiscAdj_Linear_Prec, Linear_Solver_Prec_Map, ILU); /* DESCRIPTION: Linear solver for the discete adjoint systems */ + addDoubleOption("MAX_UPDATE_CAP_FLOW", Max_Update_Flow, 0.2); + /* DESCRIPTION: Max value for under-relaxation cap for density and energy variables */ + addDoubleOption("MAX_UPDATE_CAP_SA", Max_Update_SA, 0.99); + /* DESCRIPTION: Max value for under-relaxation cap for SA turbulence variables */ + addDoubleOption("MAX_UPDATE_CAP_SST", Max_Update_SST, 1.0); + /* DESCRIPTION: Max value for under-relaxation cap for SST turbulence variables */ /*!\par CONFIG_CATEGORY: Convergence\ingroup Config*/ /*--- Options related to convergence ---*/ @@ -1998,6 +2012,8 @@ void CConfig::SetConfig_Options() { /*!\brief CONV_NUM_METHOD_TURB * \n DESCRIPTION: Convective numerical method \ingroup Config*/ addConvectOption("CONV_NUM_METHOD_TURB", Kind_ConvNumScheme_Turb, Kind_Centered_Turb, Kind_Upwind_Turb); + /*!\brief USE_ACCURATE_TURB_JACOBIANS \n DESCRIPTION: Use numerically computed Jacobians for Standard SA model \ingroup Config*/ + addBoolOption("USE_ACCURATE_TURB_JACOBIANS", Use_Accurate_Turb_Jacobians, false); /*!\brief MUSCL_ADJTURB \n DESCRIPTION: Check if the MUSCL scheme should be used \ingroup Config*/ addBoolOption("MUSCL_ADJTURB", MUSCL_AdjTurb, false); @@ -2793,6 +2809,12 @@ void CConfig::SetConfig_Options() { /* DESCRIPTION: Preaccumulation in the AD mode. */ addBoolOption("PREACC", AD_Preaccumulation, YES); + /* DESCRIPTION: Specify the tape which is checked in a tape debug run. */ + addEnumOption("CHECK_TAPE_TYPE", AD_CheckTapeType, CheckTapeType_Map, CHECK_TAPE_TYPE::FULL_SOLVER); + + /* DESCRIPTION: Specify the tape which is checked in a tape debug run. */ + addEnumOption("CHECK_TAPE_VARIABLES", AD_CheckTapeVariables, CheckTapeVariables_Map, CHECK_TAPE_VARIABLES::SOLVER_VARIABLES); + /*--- options that are used in the python optimization scripts. These have no effect on the c++ toolsuite ---*/ /*!\par CONFIG_CATEGORY:Python Options\ingroup Config*/ @@ -3413,7 +3435,12 @@ 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] = Time_Domain ? 1 : nOuterIter; + } + else { + VolumeOutputFrequencies[iVolumeFreq] = Time_Domain ? 1 : 250; + } } } else if (nVolumeOutputFrequencies < nVolumeOutputFiles) { /*--- If there are fewer frequencies than files, repeat the last frequency. @@ -8378,7 +8405,7 @@ string CConfig::GetFilename(string filename, const string& ext, int timeIter) co filename = filename + string(ext); /*--- Append the zone number if multizone problems ---*/ - if (Multizone_Problem) + if (Multizone_Problem && Multizone_Adapt_FileName) filename = GetMultizone_FileName(filename, GetiZone(), ext); /*--- Append the zone number if multiple instance problems ---*/ diff --git a/Common/src/geometry/CGeometry.cpp b/Common/src/geometry/CGeometry.cpp index 742bf4be65ae..dd962f3f9701 100644 --- a/Common/src/geometry/CGeometry.cpp +++ b/Common/src/geometry/CGeometry.cpp @@ -792,15 +792,20 @@ void CGeometry::PreprocessPeriodicComms(CGeometry* geometry, CConfig* config) { nPoint_Send_All[0] = 0; int* nPoint_Recv_All = new int[size + 1]; nPoint_Recv_All[0] = 0; - int* nPoint_Flag = new int[size]; - for (iRank = 0; iRank < size; iRank++) { - nPoint_Send_All[iRank] = 0; - nPoint_Recv_All[iRank] = 0; - nPoint_Flag[iRank] = -1; - } - nPoint_Send_All[size] = 0; - nPoint_Recv_All[size] = 0; + /*--- Store a set of unique (point, marker) pairs per destination rank. ---*/ + using PointMarkerPair = std::pair; + + auto pairHash = [](const PointMarkerPair& p) -> std::size_t { + // Use the golden ratio constant and bit mixing to generate pair hash + std::size_t h1 = std::hash{}(p.first); + std::size_t h2 = std::hash{}(p.second); + + return h1 ^ (h2 + 0x9e3779b9 + (h1 << 6) + (h1 >> 2)); + }; + + using PointMarkerSet = std::unordered_set; + std::vector Points_Send_All(size, PointMarkerSet(0, pairHash)); /*--- Loop through all of our periodic markers and track our sends with each rank. ---*/ @@ -821,19 +826,17 @@ void CGeometry::PreprocessPeriodicComms(CGeometry* geometry, CConfig* config) { iRank = static_cast(geometry->vertex[iMarker][iVertex]->GetDonorProcessor()); - /*--- If we have not visited this point last, increment our - number of points that must be sent to a particular proc. ---*/ - - if ((nPoint_Flag[iRank] != static_cast(iPoint))) { - nPoint_Flag[iRank] = static_cast(iPoint); - nPoint_Send_All[iRank + 1] += 1; - } + /*--- Store the (point, marker) pair in the set for the destination rank. ---*/ + Points_Send_All[iRank].insert(std::make_pair(iPoint, static_cast(iMarker))); } } } } - delete[] nPoint_Flag; + for (iRank = 0; iRank < size; iRank++) { + nPoint_Send_All[iRank + 1] = Points_Send_All[iRank].size(); + nPoint_Recv_All[iRank + 1] = 0; + } /*--- Communicate the number of points to be sent/recv'd amongst all processors. After this communication, each proc knows how @@ -943,11 +946,15 @@ void CGeometry::PreprocessPeriodicComms(CGeometry* geometry, CConfig* config) { auto* idSend = new unsigned long[nPoint_PeriodicSend[nPeriodicSend] * nPackets]; for (iSend = 0; iSend < nPoint_PeriodicSend[nPeriodicSend] * nPackets; iSend++) idSend[iSend] = 0; + /*--- Re-use set of unique points per destination rank. ---*/ + for (iRank = 0; iRank < size; iRank++) Points_Send_All[iRank].clear(); + /*--- Build the lists of local index and periodic marker index values. ---*/ ii = 0; jj = 0; for (iSend = 0; iSend < nPeriodicSend; iSend++) { + int destRank = Neighbors_PeriodicSend[iSend]; for (iMarker = 0; iMarker < config->GetnMarker_All(); iMarker++) { if (config->GetMarker_All_KindBC(iMarker) == PERIODIC_BOUNDARY) { iPeriodic = config->GetMarker_All_PerBound(iMarker); @@ -969,14 +976,21 @@ void CGeometry::PreprocessPeriodicComms(CGeometry* geometry, CConfig* config) { index on the matching periodic point and the periodic marker index to be communicated to the recv rank. ---*/ - if (iRank == Neighbors_PeriodicSend[iSend]) { - Local_Point_PeriodicSend[ii] = iPoint; - Local_Marker_PeriodicSend[ii] = static_cast(iMarker); - jj = ii * nPackets; - idSend[jj] = geometry->vertex[iMarker][iVertex]->GetDonorPoint(); - jj++; - idSend[jj] = static_cast(iPeriodic); - ii++; + if (iRank == destRank) { + /*--- Check if we have already added this (point, marker) pair + for this destination rank. Use the result if insert(), which is + a pair whose second element is success. ---*/ + const auto pointMarkerPair = std::make_pair(iPoint, static_cast(iMarker)); + const auto insertResult = Points_Send_All[destRank].insert(pointMarkerPair); + if (insertResult.second) { + Local_Point_PeriodicSend[ii] = iPoint; + Local_Marker_PeriodicSend[ii] = static_cast(iMarker); + jj = ii * nPackets; + idSend[jj] = geometry->vertex[iMarker][iVertex]->GetDonorPoint(); + jj++; + idSend[jj] = static_cast(iPeriodic); + ii++; + } } } } diff --git a/Common/src/geometry/dual_grid/CPoint.cpp b/Common/src/geometry/dual_grid/CPoint.cpp index 675fa232ab84..16da2e70b605 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 dde873162087..8eec26d2e82b 100644 --- a/SU2_CFD/include/drivers/CDiscAdjMultizoneDriver.hpp +++ b/SU2_CFD/include/drivers/CDiscAdjMultizoneDriver.hpp @@ -69,10 +69,10 @@ 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 + FULL_SOLVER_TAPE, /*!< \brief Entire derivative information for a coupled adjoint solution update. */ OBJECTIVE_FUNCTION_TAPE, /*!< \brief Record only the dependence of the objective function w.r.t. solver variables (from all zones). */ @@ -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. */ @@ -141,6 +141,16 @@ class CDiscAdjMultizoneDriver : public CMultizoneDriver { */ void StartSolver() override; + /*! + * \brief [Overload] Launch the tape test mode for the discrete adjoint multizone solver. + */ + void TapeTest (); + + /*! + * \brief [Overload] Get error numbers after a tape test run of the discrete adjoint multizone solver. + */ + int TapeTestGatherErrors(AD::ErrorReport& error_report) const; + /*! * \brief Preprocess the multizone iteration */ diff --git a/SU2_CFD/include/numerics/CNumerics.hpp b/SU2_CFD/include/numerics/CNumerics.hpp index 2ad25b65f1b6..1ee3a2c56b6c 100644 --- a/SU2_CFD/include/numerics/CNumerics.hpp +++ b/SU2_CFD/include/numerics/CNumerics.hpp @@ -190,6 +190,8 @@ class CNumerics { bool bounded_scalar = false; /*!< \brief Flag for bounded scalar problem */ + su2double DiffCoeff_kw[6]; /*!< \brief Storage for diffusion coefficient*/ + public: /*! * \brief Return type used in some "ComputeResidual" overloads to give a @@ -706,6 +708,7 @@ class CNumerics { * \param[in] val_CDkw_i - Value of the cross diffusion at point i. */ virtual void SetCrossDiff(su2double val_CDkw_i) {/* empty */}; + inline su2double GetDiffCoeff_kw(int index) const {return DiffCoeff_kw[index];} /*! * \brief Set the value of the effective intermittency for the LM model. diff --git a/SU2_CFD/include/numerics/turbulent/turb_diffusion.hpp b/SU2_CFD/include/numerics/turbulent/turb_diffusion.hpp index b8edcd9007a7..4a9aff95dacb 100644 --- a/SU2_CFD/include/numerics/turbulent/turb_diffusion.hpp +++ b/SU2_CFD/include/numerics/turbulent/turb_diffusion.hpp @@ -52,6 +52,9 @@ class CAvgGrad_TurbSA final : public CAvgGrad_Scalar { using Base::Jacobian_j; const su2double sigma = 2.0/3.0; + const su2double cb2 = 0.622; + + const bool use_accurate_jacobians; /*! * \brief Adds any extra variables to AD @@ -67,17 +70,39 @@ class CAvgGrad_TurbSA final : public CAvgGrad_Scalar { /*--- Compute mean effective viscosity ---*/ + /*--- First Term. Normal diffusion, and conservative part of the quadratic diffusion. + * ||grad nu_t||^2 = div(nu_t grad nu_t) - nu_t div grad nu_t ---*/ const su2double nu_i = Laminar_Viscosity_i/Density_i; const su2double nu_j = Laminar_Viscosity_j/Density_j; - const su2double nu_e = 0.5*(nu_i+nu_j+ScalarVar_i[0]+ScalarVar_j[0]); + const su2double nu_e = 0.5 * (nu_i + nu_j + (1 + cb2) * (ScalarVar_i[0] + ScalarVar_j[0])); + const su2double term_1 = nu_e; - Flux[0] = nu_e*Proj_Mean_GradScalarVar[0]/sigma; + /* Second Term (quadratic diffusion, non conservative). */ + const su2double nu_tilde_i = ScalarVar_i[0]; + const su2double term_2 = cb2 * nu_tilde_i; - /*--- For Jacobians -> Use of TSL approx. to compute derivatives of the gradients ---*/ + const su2double diffusion_coefficient = term_1 - term_2; + Flux[0] = diffusion_coefficient * Proj_Mean_GradScalarVar[0] / sigma; if (implicit) { - Jacobian_i[0][0] = (0.5*Proj_Mean_GradScalarVar[0]-nu_e*proj_vector_ij)/sigma; - Jacobian_j[0][0] = (0.5*Proj_Mean_GradScalarVar[0]+nu_e*proj_vector_ij)/sigma; + /*--- For Jacobians -> Use of TSL approx. to compute derivatives of the gradients ---*/ + Jacobian_i[0][0] = -diffusion_coefficient * proj_vector_ij / sigma; + Jacobian_j[0][0] = diffusion_coefficient * proj_vector_ij / sigma; + + if (use_accurate_jacobians) { + /*--- The diffusion coefficient is also a function of nu_t. ---*/ + const su2double dTerm1_dnut_i = (1 + cb2) * 0.5; + const su2double dTerm1_dnut_j = (1 + cb2) * 0.5; + + const su2double dTerm2_dnut_i = cb2; + const su2double dTerm2_dnut_j = 0.0; + + const su2double dDC_dnut_i = dTerm1_dnut_i - dTerm2_dnut_i; + const su2double dDC_dnut_j = dTerm1_dnut_j - dTerm2_dnut_j; + + Jacobian_i[0][0] += dDC_dnut_i * Proj_Mean_GradScalarVar[0] / sigma; + Jacobian_j[0][0] += dDC_dnut_j * Proj_Mean_GradScalarVar[0] / sigma; + } } } @@ -91,7 +116,8 @@ class CAvgGrad_TurbSA final : public CAvgGrad_Scalar { */ CAvgGrad_TurbSA(unsigned short val_nDim, unsigned short val_nVar, bool correct_grad, const CConfig* config) - : CAvgGrad_Scalar(val_nDim, val_nVar, correct_grad, config) {} + : CAvgGrad_Scalar(val_nDim, val_nVar, correct_grad, config), + use_accurate_jacobians(config->GetUse_Accurate_Turb_Jacobians()) {} }; /*! @@ -118,6 +144,7 @@ class CAvgGrad_TurbSA_Neg final : public CAvgGrad_Scalar { const su2double sigma = 2.0/3.0; const su2double cn1 = 16.0; + const su2double cb2 = 0.622; /*! * \brief Adds any extra variables to AD @@ -136,27 +163,39 @@ class CAvgGrad_TurbSA_Neg final : public CAvgGrad_Scalar { const su2double nu_i = Laminar_Viscosity_i/Density_i; const su2double nu_j = Laminar_Viscosity_j/Density_j; - const su2double nu_ij = 0.5*(nu_i+nu_j); - const su2double nu_tilde_ij = 0.5*(ScalarVar_i[0] + ScalarVar_j[0]); - - su2double nu_e; - - if (nu_tilde_ij > 0.0) { - nu_e = nu_ij + nu_tilde_ij; - } - else { - const su2double Xi = nu_tilde_ij/nu_ij; - const su2double fn = (cn1 + Xi*Xi*Xi)/(cn1 - Xi*Xi*Xi); - nu_e = nu_ij + fn*nu_tilde_ij; + const su2double nu_ij = 0.5 * (nu_i + nu_j); + const su2double nu_tilde_i = ScalarVar_i[0]; + const su2double nu_tilde_j = ScalarVar_j[0]; + const su2double nu_tilde_ij = 0.5 * (nu_tilde_i + nu_tilde_j); + + /*--- Following Diskin's implementation from 10.2514/1.J064629, they propose a new fn function + * to be evaluated at the cell to maintain positivity in the diffusion coefficient, which is + * used in both terms. The new fn term averaged across the face reverts to the original fn + * function. ---*/ + + /*--- Second Term (LHS) ---*/ + const su2double zeta_i = ((1 + cb2) * nu_tilde_ij - cb2 * nu_tilde_i) / nu_ij; + su2double fn_i = 1.0; + if (zeta_i < 0.0) { + fn_i = (cn1 + pow(zeta_i,3)) / (cn1 - pow(zeta_i,3)); } - Flux[0] = nu_e*Proj_Mean_GradScalarVar[0]/sigma; + const su2double term_1 = (nu_ij + (1 + cb2) * nu_tilde_ij * fn_i); + const su2double term_2 = cb2 * nu_tilde_i * fn_i; + Flux[0] = (term_1 - term_2) * Proj_Mean_GradScalarVar[0] / sigma; - /*--- For Jacobians -> Use of TSL approx. to compute derivatives of the gradients ---*/ + /*--- For Jacobians -> Use of TSL approx. to compute derivatives of the gradients + * Exact Jacobians were tested on multiple cases but resulted in divergence of all + * simulations, hence only frozen diffusion coefficient (approximate) Jacobians are used. ---*/ if (implicit) { - Jacobian_i[0][0] = (0.5*Proj_Mean_GradScalarVar[0]-nu_e*proj_vector_ij)/sigma; - Jacobian_j[0][0] = (0.5*Proj_Mean_GradScalarVar[0]+nu_e*proj_vector_ij)/sigma; + const su2double diffusion_coefficient = (term_1 - term_2); + + const su2double dGrad_dnut_i = -proj_vector_ij; + const su2double dGrad_dnut_j = proj_vector_ij; + + Jacobian_i[0][0] = diffusion_coefficient * dGrad_dnut_i / sigma; + Jacobian_j[0][0] = diffusion_coefficient * dGrad_dnut_j / sigma; } } @@ -201,6 +240,7 @@ class CAvgGrad_TurbSST final : public CAvgGrad_Scalar { const su2double sigma_k2; const su2double sigma_om1; const su2double sigma_om2; + const bool use_accurate_jacobians; su2double F1_i, F1_j; /*!< \brief Menter's first blending function */ @@ -231,20 +271,60 @@ class CAvgGrad_TurbSST final : public CAvgGrad_Scalar { const su2double diff_j_omega = Laminar_Viscosity_j + sigma_omega_j*Eddy_Viscosity_j; const su2double diff_kine = 0.5*(diff_i_kine + diff_j_kine); - const su2double diff_omega = 0.5*(diff_i_omega + diff_j_omega); - + const su2double diff_omega_T1 = 0.5*(diff_i_omega + diff_j_omega); + + /*--- We aim to treat the cross-diffusion as a diffusion treat rather than a source term. + * Re-writing the cross-diffusion contribution as λ/w ∇w ∇k, where λ = (2 (1- F1) ρ σ_ω2) + * and expanding using the product rule for divergence theorem gives: ∇(λ∇k) - w ∇(λ∇k). + * Discretising using FVM, gives: (λ_dw)_ij ∇k - w ∇(λ∇k); where λ_dw is λ_ij/w_ij. ---*/ + + const su2double lambda_i = 2 * (1 - F1_i) * Density_i * sigma_omega_i; + const su2double lambda_j = 2 * (1 - F1_j) * Density_j * sigma_omega_j; + const su2double lambda_ij = 0.5 * (lambda_i + lambda_j); + const su2double w_ij = 0.5 * (ScalarVar_i[1] + ScalarVar_j[1]); + + const su2double diff_omega_T2 = lambda_ij; + + const su2double diff_omega_T3 = -ScalarVar_i[1] * lambda_ij/w_ij; + + const su2double diff_omega = diff_omega_T1 + diff_omega_T2 + diff_omega_T3; + + /* Store some variables for debugging, REMOVE LATER.*/ + CNumerics::DiffCoeff_kw[0] = diff_kine; + CNumerics::DiffCoeff_kw[1] = diff_omega; + CNumerics::DiffCoeff_kw[2] = diff_omega_T1; + CNumerics::DiffCoeff_kw[3] = diff_omega_T2; + CNumerics::DiffCoeff_kw[4] = diff_omega_T3; + CNumerics::DiffCoeff_kw[5] = lambda_ij; + Flux[0] = diff_kine*Proj_Mean_GradScalarVar[0]; - Flux[1] = diff_omega*Proj_Mean_GradScalarVar[1]; + Flux[1] = diff_omega_T1*Proj_Mean_GradScalarVar[1] + (diff_omega_T2 + diff_omega_T3)*Proj_Mean_GradScalarVar[0]; /*--- For Jacobians -> Use of TSL (Thin Shear Layer) approx. to compute derivatives of the gradients ---*/ if (implicit) { const su2double proj_on_rho_i = proj_vector_ij/Density_i; - Jacobian_i[0][0] = -diff_kine*proj_on_rho_i; Jacobian_i[0][1] = 0.0; - Jacobian_i[1][0] = 0.0; Jacobian_i[1][1] = -diff_omega*proj_on_rho_i; - const su2double proj_on_rho_j = proj_vector_ij/Density_j; - Jacobian_j[0][0] = diff_kine*proj_on_rho_j; Jacobian_j[0][1] = 0.0; - Jacobian_j[1][0] = 0.0; Jacobian_j[1][1] = diff_omega*proj_on_rho_j; + Jacobian_i[0][0] = -diff_kine*proj_on_rho_i; + Jacobian_i[0][1] = 0.0; + Jacobian_i[1][0] = (diff_omega_T2+diff_omega_T3)*-proj_on_rho_i; + Jacobian_i[1][1] = -diff_omega_T1*proj_on_rho_i; + + Jacobian_j[0][0] = diff_kine*proj_on_rho_j; + Jacobian_j[0][1] = 0.0; + Jacobian_j[1][0] = (diff_omega_T2+diff_omega_T3)*proj_on_rho_j; + Jacobian_j[1][1] = diff_omega_T1*proj_on_rho_j; + + if (use_accurate_jacobians) { + Jacobian_i[0][0] = -diff_kine*proj_on_rho_i; + Jacobian_i[0][1] = 0.0; + Jacobian_i[1][0] = (diff_omega_T2 + diff_omega_T3)*-proj_on_rho_i; + Jacobian_i[1][1] = -proj_on_rho_i * diff_omega_T1 - 2*lambda_ij*ScalarVar_j[1]/pow(ScalarVar_i[1]+ScalarVar_j[1],2) * Proj_Mean_GradScalarVar[0]; + + Jacobian_j[0][0] = diff_kine*proj_on_rho_j; + Jacobian_j[0][1] = 0.0; + Jacobian_j[1][0] = (diff_omega_T2 + diff_omega_T3)*proj_on_rho_j; + Jacobian_j[1][1] = proj_on_rho_j * diff_omega_T1 + 2*lambda_ij*ScalarVar_i[1]/pow(ScalarVar_i[1]+ScalarVar_j[1],2) * Proj_Mean_GradScalarVar[0]; + } } } @@ -263,7 +343,8 @@ class CAvgGrad_TurbSST final : public CAvgGrad_Scalar { sigma_k1(constants[0]), sigma_k2(constants[1]), sigma_om1(constants[2]), - sigma_om2(constants[3]) { + sigma_om2(constants[3]), + use_accurate_jacobians(config->GetUse_Accurate_Turb_Jacobians()) { } /*! diff --git a/SU2_CFD/include/numerics/turbulent/turb_sources.hpp b/SU2_CFD/include/numerics/turbulent/turb_sources.hpp index 4fa9080e1ad5..11af3f6e7245 100644 --- a/SU2_CFD/include/numerics/turbulent/turb_sources.hpp +++ b/SU2_CFD/include/numerics/turbulent/turb_sources.hpp @@ -86,7 +86,7 @@ class CSourceBase_TurbSA : public CNumerics { */ inline void ResidualAxisymmetricDiffusion(su2double sigma) { if (Coord_i[1] < EPS) return; - + const su2double yinv = 1.0 / Coord_i[1]; const su2double& nue = ScalarVar_i[0]; @@ -243,14 +243,14 @@ class CSourceBase_TurbSA : public CNumerics { var.interDestrFactor = 1.0; } - /*--- Compute production, destruction and cross production and jacobian ---*/ - su2double Production = 0.0, Destruction = 0.0, CrossProduction = 0.0; - SourceTerms::get(ScalarVar_i[0], var, Production, Destruction, CrossProduction, Jacobian_i[0]); + /*--- Compute production, destruction and jacobian ---*/ + su2double Production = 0.0, Destruction = 0.0; + SourceTerms::get(ScalarVar_i[0], var, Production, Destruction, Jacobian_i[0]); - Residual = (Production - Destruction + CrossProduction) * Volume; + Residual = (Production - Destruction) * Volume; if (axisymmetric) ResidualAxisymmetricDiffusion(var.sigma); - + Jacobian_i[0] *= Volume; } @@ -423,13 +423,13 @@ struct Edw { }; /*! - * \brief SA source terms classes: production, destruction and cross-productions term and their derivative. + * \brief SA source terms classes: production and destruction term and their derivative. + * \note Quadratic diffusion is included in the viscous fluxes. * \ingroup SourceDiscr * \param[in] nue: SA variable. * \param[in] var: Common SA variables struct. * \param[out] production: Production term. * \param[out] destruction: Destruction term. - * \param[out] cross_production: CrossProduction term. * \param[out] jacobian: Derivative of the combined source term wrt nue. */ struct SourceTerms { @@ -437,10 +437,9 @@ struct SourceTerms { /*! \brief Baseline (Original SA model). */ struct Bsl { static void get(const su2double& nue, const CSAVariables& var, su2double& production, su2double& destruction, - su2double& cross_production, su2double& jacobian) { + su2double& jacobian) { ComputeProduction(nue, var, production, jacobian); ComputeDestruction(nue, var, destruction, jacobian); - ComputeCrossProduction(nue, var, cross_production, jacobian); } static void ComputeProduction(const su2double& nue, const CSAVariables& var, su2double& production, @@ -458,23 +457,17 @@ struct Bsl { jacobian -= var.interDestrFactor * ((var.cw1 * var.d_fw - cb1_k2 * var.d_ft2) * pow(nue, 2) + factor * 2 * nue) / var.dist_i_2; } - static void ComputeCrossProduction(const su2double& nue, const CSAVariables& var, su2double& cross_production, - su2double&) { - cross_production = var.cb2_sigma * var.norm2_Grad; - /*--- No contribution to the jacobian. ---*/ - } }; /*! \brief Negative. */ struct Neg { static void get(const su2double& nue, const CSAVariables& var, su2double& production, su2double& destruction, - su2double& cross_production, su2double& jacobian) { + su2double& jacobian) { if (nue > 0.0) { - Bsl::get(nue, var, production, destruction, cross_production, jacobian); + Bsl::get(nue, var, production, destruction, jacobian); } else { ComputeProduction(nue, var, production, jacobian); ComputeDestruction(nue, var, destruction, jacobian); - ComputeCrossProduction(nue, var, cross_production, jacobian); } } @@ -493,10 +486,6 @@ struct Neg { jacobian -= 2 * dD_dnu * var.interDestrFactor; } - static void ComputeCrossProduction(const su2double& nue, const CSAVariables& var, su2double& cross_production, - su2double& jacobian) { - Bsl::ComputeCrossProduction(nue, var, cross_production, jacobian); - } }; }; @@ -562,7 +551,7 @@ class CCompressibilityCorrection final : public ParentClass { const su2double v = V_i[idx.Velocity() + 1]; const su2double d_axiCorrection = 2.0 * c5 * nue * pow(v * yinv / sound_speed, 2) * Volume; - const su2double axiCorrection = 0.5 * nue * d_axiCorrection; + const su2double axiCorrection = 0.5 * nue * d_axiCorrection; this->Residual -= axiCorrection; this->Jacobian_i[0] -= d_axiCorrection; @@ -924,7 +913,7 @@ class CSourcePieceWise_TurbSST final : public CNumerics { /*--- Cross diffusion ---*/ - Residual[1] += (1.0 - F1_i) * CDkw_i * Volume; + // Residual[1] += (1.0 - F1_i) * CDkw_i * Volume; /*--- Contribution due to 2D axisymmetric formulation ---*/ diff --git a/SU2_CFD/include/solvers/CDiscAdjSolver.hpp b/SU2_CFD/include/solvers/CDiscAdjSolver.hpp index 7071c31750f9..d196eb1bff4f 100644 --- a/SU2_CFD/include/solvers/CDiscAdjSolver.hpp +++ b/SU2_CFD/include/solvers/CDiscAdjSolver.hpp @@ -57,7 +57,7 @@ class CDiscAdjSolver final : public CSolver { su2double Total_Sens_BPress; /*!< \brief Total sensitivity to outlet pressure. */ su2double Total_Sens_Density; /*!< \brief Total sensitivity to initial density (incompressible). */ su2double Total_Sens_ModVel; /*!< \brief Total sensitivity to inlet velocity (incompressible). */ - su2double Mach, Alpha, Beta, Pressure, Temperature, BPressure, ModVel; + su2double Mach, Alpha, Beta, Temperature, BPressure, ModVel; su2double TemperatureRad, Total_Sens_Temp_Rad; CDiscAdjVariable* nodes = nullptr; /*!< \brief The highest level in the variable hierarchy this solver can safely use. */ diff --git a/SU2_CFD/include/solvers/CFVMFlowSolverBase.inl b/SU2_CFD/include/solvers/CFVMFlowSolverBase.inl index d3673e1b763a..bb9099934282 100644 --- a/SU2_CFD/include/solvers/CFVMFlowSolverBase.inl +++ b/SU2_CFD/include/solvers/CFVMFlowSolverBase.inl @@ -558,7 +558,7 @@ void CFVMFlowSolverBase::ComputeUnderRelaxationFactor(const CConfig* confi /* Loop over the solution update given by relaxing the linear system for this nonlinear iteration. */ - const su2double allowableRatio = 0.2; + const su2double allowableRatio = config->GetUnderRelax_Flow(); SU2_OMP_FOR_STAT(omp_chunk_size) for (unsigned long iPoint = 0; iPoint < nPointDomain; iPoint++) { @@ -671,12 +671,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(); + strainMax = SU2_TYPE::GetValue(max(strainMax, StrainMag(iPoint))); + omegaMax = SU2_TYPE::GetValue(max(omegaMax, GeometryToolbox::Norm(3, Vorticity))); } END_SU2_OMP_FOR diff --git a/SU2_CFD/include/solvers/CScalarSolver.hpp b/SU2_CFD/include/solvers/CScalarSolver.hpp index 0fc61aa7ab25..363226ebd8c3 100644 --- a/SU2_CFD/include/solvers/CScalarSolver.hpp +++ b/SU2_CFD/include/solvers/CScalarSolver.hpp @@ -77,6 +77,7 @@ class CScalarSolver : public CSolver { /*--- Edge fluxes for reducer strategy (see the notes in CEulerSolver.hpp). ---*/ CSysVector EdgeFluxes; /*!< \brief Flux across each edge. */ + CSysVector EdgeFluxesDiff; /*!< \brief Flux difference between ij and ji for non-conservative discretisation. */ /*! * \brief The highest level in the variable hierarchy this solver can safely use. diff --git a/SU2_CFD/include/solvers/CScalarSolver.inl b/SU2_CFD/include/solvers/CScalarSolver.inl index 314354116ff1..bffe6f49dd0d 100644 --- a/SU2_CFD/include/solvers/CScalarSolver.inl +++ b/SU2_CFD/include/solvers/CScalarSolver.inl @@ -341,15 +341,21 @@ void CScalarSolver::Upwind_Residual(CGeometry* geometry, CSolver** template void CScalarSolver::SumEdgeFluxes(CGeometry* geometry) { + const bool nonConservative = EdgeFluxesDiff.GetLocSize() > 0; + SU2_OMP_FOR_STAT(omp_chunk_size) for (unsigned long iPoint = 0; iPoint < nPoint; ++iPoint) { LinSysRes.SetBlock_Zero(iPoint); for (auto iEdge : geometry->nodes->GetEdges(iPoint)) { - if (iPoint == geometry->edges->GetNode(iEdge, 0)) + if (iPoint == geometry->edges->GetNode(iEdge, 0)) { LinSysRes.AddBlock(iPoint, EdgeFluxes.GetBlock(iEdge)); - else + } else { LinSysRes.SubtractBlock(iPoint, EdgeFluxes.GetBlock(iEdge)); + if (nonConservative) { + LinSysRes.SubtractBlock(iPoint, EdgeFluxesDiff.GetBlock(iEdge)); + } + } } } END_SU2_OMP_FOR diff --git a/SU2_CFD/include/solvers/CSolver.hpp b/SU2_CFD/include/solvers/CSolver.hpp index 417e57d08178..1e744e1c4db7 100644 --- a/SU2_CFD/include/solvers/CSolver.hpp +++ b/SU2_CFD/include/solvers/CSolver.hpp @@ -191,12 +191,12 @@ class CSolver { CSysVector LinSysSol; /*!< \brief vector to store iterative solution of implicit linear system. */ CSysVector LinSysRes; /*!< \brief vector to store iterative residual of implicit linear system. */ #ifndef CODI_FORWARD_TYPE - CSysMatrix Jacobian; /*!< \brief Complete sparse Jacobian structure for implicit computations. */ - CSysSolve System; /*!< \brief Linear solver/smoother. */ + using JacobianScalarType = su2mixedfloat; #else - CSysMatrix Jacobian; - CSysSolve System; + using JacobianScalarType = su2double; #endif + CSysMatrix Jacobian; /*!< \brief Complete sparse Jacobian structure for implicit computations. */ + CSysSolve System; /*!< \brief Linear solver/smoother. */ CSysVector OutputVariables; /*!< \brief vector to store the extra variables to be written. */ string* OutputHeadingNames; /*!< \brief vector of strings to store the headings for the exra variables */ diff --git a/SU2_CFD/include/solvers/CTurbSSTSolver.hpp b/SU2_CFD/include/solvers/CTurbSSTSolver.hpp index 5bdc3549aed3..352270f6011a 100644 --- a/SU2_CFD/include/solvers/CTurbSSTSolver.hpp +++ b/SU2_CFD/include/solvers/CTurbSSTSolver.hpp @@ -53,6 +53,14 @@ class CTurbSSTSolver final : public CTurbSolver { CSolver **solver_container, const CConfig *config, unsigned short val_marker); + + /*! + * \brief Compute a suitable under-relaxation parameter to limit the change in the solution variables over + * a nonlinear iteration for stability. + * \param[in] config - Definition of the particular problem. + */ + void ComputeUnderRelaxationFactor(const CConfig *config); + public: /*! * \brief Constructor. diff --git a/SU2_CFD/include/variables/CTurbVariable.hpp b/SU2_CFD/include/variables/CTurbVariable.hpp index 2e9de32fbb24..8dc37e9910c4 100644 --- a/SU2_CFD/include/variables/CTurbVariable.hpp +++ b/SU2_CFD/include/variables/CTurbVariable.hpp @@ -43,6 +43,28 @@ class CTurbVariable : public CScalarVariable { static constexpr size_t MAXNVAR = 2; VectorType turb_index; VectorType intermittency; /*!< \brief Value of the intermittency for the trans. model. */ + VectorType DC_TKE; /*!< \brief TKE Diffusion Coefficient. */ + VectorType DC_Omega; /*!< \brief Omega Diffusion Coefficient. */ + VectorType DC_OmegaT1; + VectorType DC_OmegaT2; + VectorType DC_OmegaT3; + VectorType DC_Omega_kc; + VectorType DC_Omega_so2i; + VectorType DC_Omega_Coi; + VectorType DC_Omega_Coj; + VectorType DC_Omega_F1i; + VectorType DC_Omega_F1j; + VectorType DC_Omega_so2j; + VectorType DC_Omega_mui; + VectorType DC_Omega_muj; + VectorType DC_Omega_muti; + VectorType DC_Omega_mutj; + VectorType DC_Omega_rhoi; + VectorType DC_Omega_rhoj; + VectorType DC_Omega_ki; + VectorType DC_Omega_kj; + VectorType DC_Omega_wi; + VectorType DC_Omega_wj; /*! * \brief Constructor of the class. @@ -100,6 +122,68 @@ class CTurbVariable : public CScalarVariable { */ inline void SetIntermittency(unsigned long iPoint, su2double val_intermittency) final { intermittency(iPoint) = val_intermittency; } + /*! + * \brief Set the Diffusion Coefficients of TKE and omega equations. + * \param[in] iPoint - Point index. + * \param[in] val_DC_kw - diffusion coefficient value + */ + inline void SetDiffCoeff_kw(unsigned long iPoint, su2double* val_DC_kw) { + DC_TKE(iPoint) = val_DC_kw[0]; + DC_Omega(iPoint) = val_DC_kw[1]; + DC_OmegaT1(iPoint) = val_DC_kw[2]; + DC_OmegaT2(iPoint) = val_DC_kw[3]; + DC_OmegaT3(iPoint) = val_DC_kw[4]; + DC_Omega_kc(iPoint) = val_DC_kw[5]; + } + + /*! + * \brief Get the diffusion coefficient value of TKE. + * \param[in] iPoint - Point index. + */ + inline su2double GetTKE_DC(unsigned long iPoint) const final { + return DC_TKE(iPoint); + } + + /*! + * \brief Get the diffusion coefficient value of omega. + * \param[in] iPoint - Point index. + */ + inline su2double GetOmega_DC(unsigned long iPoint) const final { + return DC_Omega(iPoint); + } + + /*! + * \brief Get omega diffusion coefficient Term 1. + * \param[in] iPoint - Point index. + */ + inline su2double GetOmega_DCT1(unsigned long iPoint) const final { + return DC_OmegaT1(iPoint); + } + + /*! + * \brief Get omega diffusion coefficient Term 2. + * \param[in] iPoint - Point index. + */ + inline su2double GetOmega_DCT2(unsigned long iPoint) const final { + return DC_OmegaT2(iPoint); + } + + /*! + * \brief Get omega diffusion coefficient Term 3. + * \param[in] iPoint - Point index. + */ + inline su2double GetOmega_DCT3(unsigned long iPoint) const final { + return DC_OmegaT3(iPoint); + } + + /*! + * \brief Get the clip value kc used in omega diffusion coeffient. + * \param[in] iPoint - Point index. + */ + inline su2double GetOmega_DC_kc(unsigned long iPoint) const final { + return DC_Omega_kc(iPoint); + } + /*! * \brief Register eddy viscosity (muT) as Input or Output of an AD recording. * \param[in] input - Boolean whether In- or Output should be registered. diff --git a/SU2_CFD/include/variables/CVariable.hpp b/SU2_CFD/include/variables/CVariable.hpp index 7088d2ea5a6e..554f632f7a8e 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 before solver iteration. */ + su2matrix AD_OutputIndex; /*!< \brief Indices of Solution variables in the adjoint vector after solver iteration. */ 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,11 +133,11 @@ 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); } - void RegisterContainer(bool input, su2activevector& variable, su2vector* ad_index = nullptr) { + void RegisterContainer(bool input, su2activevector& variable, su2vector* 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) { @@ -1664,6 +1664,19 @@ class CVariable { */ inline virtual su2double GetCrossDiff(unsigned long iPoint) const { return 0.0; } + /*! + * \brief Get the value of the diffusion coefficient of tke and omega equations. + * \param[in] ipoint + * \param[in] val_diff_coeff + */ + inline virtual void SetDiffCoeff_kw(unsigned long iPoint, su2double* val_DiffCoeff) const { } + inline virtual su2double GetTKE_DC(unsigned long iPoint) const {return 0.0; } + inline virtual su2double GetOmega_DC(unsigned long iPoint) const {return 0.0; } + inline virtual su2double GetOmega_DCT1(unsigned long iPoint) const {return 0.0; } + inline virtual su2double GetOmega_DCT2(unsigned long iPoint) const {return 0.0; } + inline virtual su2double GetOmega_DCT3(unsigned long iPoint) const {return 0.0; } + inline virtual su2double GetOmega_DC_kc(unsigned long iPoint) const {return 0.0; } + /*! * \brief Get the value of the eddy viscosity. * \return the value of the eddy viscosity. @@ -2187,7 +2200,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); @@ -2195,7 +2208,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 c8dd0b1e98c3..42d2de2c1c59 100644 --- a/SU2_CFD/src/drivers/CDiscAdjMultizoneDriver.cpp +++ b/SU2_CFD/src/drivers/CDiscAdjMultizoneDriver.cpp @@ -159,6 +159,18 @@ 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()) { + + if (rank == MASTER_NODE) { + cout << "\nSU2_CFD_AD is compiled for debug mode recording. To resume the discrete adjoint solver, adjust -Dcodi-tape (-Dcodi-tape=JacobianLinear by default) and recompile." << endl; + } + Preprocess(0); + TapeTest(); + 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 +229,96 @@ void CDiscAdjMultizoneDriver::StartSolver() { } +void CDiscAdjMultizoneDriver::TapeTest() { + + if (rank == MASTER_NODE) { + cout <<"\n---------------------------- Start Debug Run ----------------------------" << endl; + } + + int total_errors = 0; + AD::ErrorReport error_report; + AD::SetTagErrorCallback(error_report); + std::ofstream out1("run1_process" + to_string(rank) + ".out"); + std::ofstream out2("run2_process" + to_string(rank) + ".out"); + + AD::ResetErrorCounter(error_report); + AD::SetDebugReportFile(error_report, &out1); + + /*--- This recording will assign the initial (same) tag to each registered variable. + * During the recording, each dependent variable will be assigned the same tag. ---*/ + + if(driver_config->GetAD_CheckTapeType() == CHECK_TAPE_TYPE::OBJECTIVE_FUNCTION) { + if(driver_config->GetAD_CheckTapeVariables() == CHECK_TAPE_VARIABLES::MESH_COORDINATES) { + if (rank == MASTER_NODE) cout << "\nChecking OBJECTIVE_FUNCTION_TAPE for SOLVER_VARIABLES_AND_MESH_COORDINATES." << endl; + SetRecording(RECORDING::TAG_INIT_SOLVER_AND_MESH, Kind_Tape::OBJECTIVE_FUNCTION_TAPE, ZONE_0); + } + else { + if (rank == MASTER_NODE) cout << "\nChecking OBJECTIVE_FUNCTION_TAPE for SOLVER_VARIABLES." << endl; + SetRecording(RECORDING::TAG_INIT_SOLVER_VARIABLES, Kind_Tape::OBJECTIVE_FUNCTION_TAPE, ZONE_0); + } + } + else { + if(driver_config->GetAD_CheckTapeVariables() == CHECK_TAPE_VARIABLES::MESH_COORDINATES) { + if (rank == MASTER_NODE) cout << "\nChecking FULL_SOLVER_TAPE for SOLVER_VARIABLES_AND_MESH_COORDINATES." << endl; + SetRecording(RECORDING::TAG_INIT_SOLVER_AND_MESH, Kind_Tape::FULL_SOLVER_TAPE, ZONE_0); + } + else { + if (rank == MASTER_NODE) cout << "\nChecking FULL_SOLVER_TAPE for SOLVER_VARIABLES." << endl; + SetRecording(RECORDING::TAG_INIT_SOLVER_VARIABLES, Kind_Tape::FULL_SOLVER_TAPE, ZONE_0); + } + } + total_errors = TapeTestGatherErrors(error_report); + + AD::ResetErrorCounter(error_report); + AD::SetDebugReportFile(error_report, &out2); + + /*--- 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. ---*/ + + if(driver_config->GetAD_CheckTapeType() == CHECK_TAPE_TYPE::OBJECTIVE_FUNCTION) { + if(driver_config->GetAD_CheckTapeVariables() == CHECK_TAPE_VARIABLES::MESH_COORDINATES) + SetRecording(RECORDING::TAG_CHECK_SOLVER_AND_MESH, Kind_Tape::OBJECTIVE_FUNCTION_TAPE, ZONE_0); + else + SetRecording(RECORDING::TAG_CHECK_SOLVER_VARIABLES, Kind_Tape::OBJECTIVE_FUNCTION_TAPE, ZONE_0); + } + else { + if(driver_config->GetAD_CheckTapeVariables() == CHECK_TAPE_VARIABLES::MESH_COORDINATES) + SetRecording(RECORDING::TAG_CHECK_SOLVER_AND_MESH, Kind_Tape::FULL_SOLVER_TAPE, ZONE_0); + else + SetRecording(RECORDING::TAG_CHECK_SOLVER_VARIABLES, Kind_Tape::FULL_SOLVER_TAPE, ZONE_0); + } + total_errors += TapeTestGatherErrors(error_report); + + if (rank == MASTER_NODE) { + cout << "\n------------------------- Tape Test Run Summary -------------------------" << endl; + cout << "\nTotal number of tape inconsistencies: " << total_errors << endl; + cout << "\n--------------------------- End Tape Test Run ---------------------------" << endl; + } +} + +int CDiscAdjMultizoneDriver::TapeTestGatherErrors(AD::ErrorReport& error_report) const { + + int num_errors = AD::GetErrorCount(error_report); + int total_errors = 0; + std::vector process_error(size); + SU2_MPI::Allreduce(&num_errors, &total_errors, 1, MPI_INT, MPI_SUM, SU2_MPI::GetComm()); + SU2_MPI::Gather(&num_errors, 1, MPI_INT, process_error.data(), 1, MPI_INT, MASTER_NODE, SU2_MPI::GetComm()); + + if (rank == MASTER_NODE) { + std::cout << "\nTotal number of detected tape inconsistencies: " << total_errors << std::endl; + if(total_errors > 0 && size > 1) { + std::cout << "\n"; + for (int irank = 0; irank < size; irank++) { + std::cout << "Number of inconsistencies from process " << irank << ": " << process_error[irank] << std::endl; + } + } + } + return total_errors; +} + bool CDiscAdjMultizoneDriver::Iterate(unsigned short iZone, unsigned long iInnerIter, bool KrylovMode) { config_container[iZone]->SetInnerIter(iInnerIter); @@ -310,6 +412,21 @@ void CDiscAdjMultizoneDriver::Run() { const bool time_domain = driver_config->GetTime_Domain(); driver_config->Set_StartTime(SU2_MPI::Wtime()); + /*--- Temporary warning because we need to test writing intermediate output to file (requires re-recording). ---*/ + for(iZone = 0; iZone < nZone; iZone++) { + for (auto iVolumeFreq = 0; iVolumeFreq < config_container[iZone]->GetnVolumeOutputFrequencies(); iVolumeFreq++){ + if (config_container[iZone]->GetVolumeOutputFrequency(iVolumeFreq) < nOuterIter) { + if (rank == MASTER_NODE) { + cout << "\nWARNING (iZone = " << iZone + << "): " + "Writing out restart files during solver iterations is not tested for the discrete adjoint multizone solver.\n" + "It is recommended to remove OUTPUT_WRT_FREQ from the config file, output files will be written when solver finalizes." << std::endl; + } + break; + } + } + } + /*--- If the gradient of the objective function is 0 so are the adjoint variables. * Unless in unsteady problems where there are other contributions to the RHS. ---*/ @@ -352,8 +469,8 @@ void CDiscAdjMultizoneDriver::Run() { * here. Otherwise, the whole tape of a coupled run will be created. ---*/ if (RecordingState != RECORDING::SOLUTION_VARIABLES) { - SetRecording(RECORDING::CLEAR_INDICES, Kind_Tape::FULL_TAPE, ZONE_0); - SetRecording(RECORDING::SOLUTION_VARIABLES, Kind_Tape::FULL_TAPE, ZONE_0); + SetRecording(RECORDING::CLEAR_INDICES, Kind_Tape::FULL_SOLVER_TAPE, ZONE_0); + SetRecording(RECORDING::SOLUTION_VARIABLES, Kind_Tape::FULL_SOLVER_TAPE, ZONE_0); } /*-- Start loop over zones. ---*/ @@ -499,11 +616,11 @@ void CDiscAdjMultizoneDriver::EvaluateSensitivities(unsigned long Iter, bool for /*--- SetRecording stores the computational graph on one iteration of the direct problem. Calling it with NONE * as argument ensures that all information from a previous recording is removed. ---*/ - SetRecording(RECORDING::CLEAR_INDICES, Kind_Tape::FULL_TAPE, ZONE_0); + SetRecording(RECORDING::CLEAR_INDICES, Kind_Tape::FULL_SOLVER_TAPE, ZONE_0); /*--- Store the computational graph of one direct iteration with the mesh coordinates as input. ---*/ - SetRecording(RECORDING::MESH_COORDS, Kind_Tape::FULL_TAPE, ZONE_0); + SetRecording(RECORDING::MESH_COORDS, Kind_Tape::FULL_SOLVER_TAPE, ZONE_0); /*--- Initialize the adjoint of the output variables of the iteration with the adjoint solution * of the current iteration. The values are passed to the AD tool. ---*/ @@ -589,6 +706,10 @@ void CDiscAdjMultizoneDriver::SetRecording(RECORDING kind_recording, Kind_Tape t case RECORDING::CLEAR_INDICES: cout << "Clearing the computational graph." << endl; break; case RECORDING::MESH_COORDS: cout << "Storing computational graph wrt MESH COORDINATES." << endl; break; case RECORDING::SOLUTION_VARIABLES: cout << "Storing computational graph wrt CONSERVATIVE VARIABLES." << endl; break; + case RECORDING::TAG_INIT_SOLVER_VARIABLES: cout << "Simulating recording with tag 1 on conservative variables." << endl; AD::SetTag(1); break; + case RECORDING::TAG_CHECK_SOLVER_VARIABLES: cout << "Checking first recording with tag 2 on conservative variables." << endl; AD::SetTag(2); break; + case RECORDING::TAG_INIT_SOLVER_AND_MESH: cout << "Simulating recording with tag 1 on conservative variables and mesh coordinates." << endl; AD::SetTag(1); break; + case RECORDING::TAG_CHECK_SOLVER_AND_MESH: cout << "Checking first recording with tag 2 on conservative variables and mesh coordinates." << endl; AD::SetTag(2); break; default: break; } } @@ -741,12 +862,12 @@ void CDiscAdjMultizoneDriver::SetObjFunction(RECORDING kind_recording) { if (rank == MASTER_NODE) { AD::RegisterOutput(ObjFunc); AD::SetIndex(ObjFunc_Index, ObjFunc); - if (kind_recording == RECORDING::SOLUTION_VARIABLES) { - cout << " Objective function : " << ObjFunc; - if (driver_config->GetWrt_AD_Statistics()){ - cout << " (" << ObjFunc_Index << ")\n"; - } - cout << endl; + if (kind_recording == RECORDING::SOLUTION_VARIABLES || + kind_recording == RECORDING::TAG_INIT_SOLVER_VARIABLES || + kind_recording == RECORDING::TAG_CHECK_SOLVER_VARIABLES || + kind_recording == RECORDING::TAG_INIT_SOLVER_AND_MESH || + kind_recording == RECORDING::TAG_CHECK_SOLVER_AND_MESH) { + cout << " Objective function : " << ObjFunc << endl; } } } diff --git a/SU2_CFD/src/iteration/CDiscAdjFluidIteration.cpp b/SU2_CFD/src/iteration/CDiscAdjFluidIteration.cpp index 7e51ff89c4b7..e0134adc0124 100644 --- a/SU2_CFD/src/iteration/CDiscAdjFluidIteration.cpp +++ b/SU2_CFD/src/iteration/CDiscAdjFluidIteration.cpp @@ -404,7 +404,13 @@ 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_SOLVER_VARIABLES || + kind_recording == RECORDING::TAG_CHECK_SOLVER_VARIABLES || + kind_recording == RECORDING::TAG_INIT_SOLVER_AND_MESH || + kind_recording == RECORDING::TAG_CHECK_SOLVER_AND_MESH || + kind_recording == RECORDING::SOLUTION_AND_MESH) { + /*--- Register flow and turbulent variables as input ---*/ if (config[iZone]->GetFluidProblem()) { @@ -429,9 +435,12 @@ void CDiscAdjFluidIteration::RegisterInput(CSolver***** solver, CGeometry**** ge } } - if (kind_recording == RECORDING::MESH_COORDS || kind_recording == RECORDING::SOLUTION_AND_MESH) { - /*--- Register node coordinates as input ---*/ + if (kind_recording == RECORDING::MESH_COORDS || + kind_recording == RECORDING::SOLUTION_AND_MESH || + kind_recording == RECORDING::TAG_INIT_SOLVER_AND_MESH || + kind_recording == RECORDING::TAG_CHECK_SOLVER_AND_MESH) { + /*--- Register node coordinates as input ---*/ geometry0->RegisterCoordinates(); } diff --git a/SU2_CFD/src/output/CFlowOutput.cpp b/SU2_CFD/src/output/CFlowOutput.cpp index e9f70fee2768..7a7c8aba263d 100644 --- a/SU2_CFD/src/output/CFlowOutput.cpp +++ b/SU2_CFD/src/output/CFlowOutput.cpp @@ -1256,6 +1256,14 @@ void CFlowOutput::SetVolumeOutputFieldsScalarSolution(const CConfig* config){ case TURB_FAMILY::KW: AddVolumeOutput("TKE", "Turb_Kin_Energy", "SOLUTION", "Turbulent kinetic energy"); AddVolumeOutput("DISSIPATION", "Omega", "SOLUTION", "Rate of dissipation"); + AddVolumeOutput("DIFF_COEFF_TKE", "TKE_Diff_Coeff", "SOLUTION", "Diffusion Coefficient for TKE equation"); + AddVolumeOutput("DIFF_COEFF_Omega", "Omega_Diff_Coeff", "SOLUTION", "Toal Diffusion Coefficient for Omega equation"); + AddVolumeOutput("DIFF_COEFF_OmegaT1", "w_DiffCoeff_T1", "SOLUTION", "Term 1 in Diffusion Coefficient for Omega equation"); + AddVolumeOutput("DIFF_COEFF_OmegaT2", "w_DiffCoeff_T2", "SOLUTION", "Term 2 in Diffusion Coefficient for Omega equation"); + AddVolumeOutput("DIFF_COEFF_OmegaT3", "w_DiffCoeff_T3", "SOLUTION", "Term 3 in Diffusion Coefficient for Omega equation"); + AddVolumeOutput("DIFF_COEFF_Omegakc", "w_DiffCoeff_kc", "SOLUTION", "kc value in Diffusion Coefficient for Omega equation"); + AddVolumeOutput("CDkw", "CDkw", "SOLUTION", "Cross-Diffusion term"); + AddVolumeOutput("F1", "F1", "SOLUTION", "F1 blending function"); break; case TURB_FAMILY::NONE: @@ -1536,6 +1544,14 @@ void CFlowOutput::LoadVolumeDataScalar(const CConfig* config, const CSolver* con case TURB_FAMILY::KW: SetVolumeOutputValue("TKE", iPoint, Node_Turb->GetSolution(iPoint, 0)); SetVolumeOutputValue("DISSIPATION", iPoint, Node_Turb->GetSolution(iPoint, 1)); + SetVolumeOutputValue("DIFF_COEFF_TKE", iPoint, Node_Turb->GetTKE_DC(iPoint)); + SetVolumeOutputValue("DIFF_COEFF_Omega", iPoint, Node_Turb->GetOmega_DC(iPoint)); + SetVolumeOutputValue("CDkw", iPoint, Node_Turb->GetCrossDiff(iPoint)); + SetVolumeOutputValue("F1", iPoint, Node_Turb->GetF1blending(iPoint)); + SetVolumeOutputValue("DIFF_COEFF_OmegaT1", iPoint, Node_Turb->GetOmega_DCT1(iPoint)); + SetVolumeOutputValue("DIFF_COEFF_OmegaT2", iPoint, Node_Turb->GetOmega_DCT2(iPoint)); + SetVolumeOutputValue("DIFF_COEFF_OmegaT3", iPoint, Node_Turb->GetOmega_DCT3(iPoint)); + SetVolumeOutputValue("DIFF_COEFF_Omegakc", iPoint, Node_Turb->GetOmega_DC_kc(iPoint)); SetVolumeOutputValue("RES_TKE", iPoint, turb_solver->LinSysRes(iPoint, 0)); SetVolumeOutputValue("RES_DISSIPATION", iPoint, turb_solver->LinSysRes(iPoint, 1)); if (limiter) { diff --git a/SU2_CFD/src/solvers/CDiscAdjSolver.cpp b/SU2_CFD/src/solvers/CDiscAdjSolver.cpp index 7a7e172a84d1..07caefe91c3c 100644 --- a/SU2_CFD/src/solvers/CDiscAdjSolver.cpp +++ b/SU2_CFD/src/solvers/CDiscAdjSolver.cpp @@ -178,17 +178,21 @@ void CDiscAdjSolver::RegisterVariables(CGeometry *geometry, CConfig *config, boo /*--- Register farfield values as input ---*/ - if((config->GetKind_Regime() == ENUM_REGIME::COMPRESSIBLE) && (KindDirect_Solver == RUNTIME_FLOW_SYS && !config->GetBoolTurbomachinery())) { + if((config->GetKind_Regime() == ENUM_REGIME::COMPRESSIBLE) && (KindDirect_Solver == RUNTIME_FLOW_SYS && !config->GetBoolTurbomachinery())){ su2double Velocity_Ref = config->GetVelocity_Ref(); Alpha = config->GetAoA()*PI_NUMBER/180.0; Beta = config->GetAoS()*PI_NUMBER/180.0; Mach = config->GetMach(); - Pressure = config->GetPressure_FreeStreamND(); - Temperature = config->GetTemperature_FreeStreamND(); + /*--- Pressure and Temperature can be registered directly via their config file value + * (no further treatment required here). ---*/ + su2double& Pressure = config->GetPressure_FreeStreamND(); + su2double& Temperature = config->GetTemperature_FreeStreamND(); su2double SoundSpeed = 0.0; + // Treat Velocity_FreeStreamND config value as non-dependent (in debug mode) + AD::ClearTagOnVariable(config->GetVelocity_FreeStreamND()[0]); if (nDim == 2) { SoundSpeed = config->GetVelocity_FreeStreamND()[0]*Velocity_Ref/(cos(Alpha)*Mach); } if (nDim == 3) { SoundSpeed = config->GetVelocity_FreeStreamND()[0]*Velocity_Ref/(cos(Alpha)*cos(Beta)*Mach); } @@ -211,11 +215,8 @@ void CDiscAdjSolver::RegisterVariables(CGeometry *geometry, CConfig *config, boo config->GetVelocity_FreeStreamND()[2] = sin(Alpha)*cos(Beta)*Mach*SoundSpeed/Velocity_Ref; } - config->SetTemperature_FreeStreamND(Temperature); direct_solver->SetTemperature_Inf(Temperature); - config->SetPressure_FreeStreamND(Pressure); direct_solver->SetPressure_Inf(Pressure); - } if ((config->GetKind_Regime() == ENUM_REGIME::COMPRESSIBLE) && (KindDirect_Solver == RUNTIME_FLOW_SYS) && config->GetBoolTurbomachinery()){ @@ -399,6 +400,9 @@ void CDiscAdjSolver::ExtractAdjoint_Variables(CGeometry *geometry, CConfig *conf if ((config->GetKind_Regime() == ENUM_REGIME::COMPRESSIBLE) && (KindDirect_Solver == RUNTIME_FLOW_SYS) && !config->GetBoolTurbomachinery()) { su2double Local_Sens_Press, Local_Sens_Temp, Local_Sens_AoA, Local_Sens_Mach; + su2double& Pressure = config->GetPressure_FreeStreamND(); + su2double& Temperature = config->GetTemperature_FreeStreamND(); + Local_Sens_Mach = SU2_TYPE::GetDerivative(Mach); Local_Sens_AoA = SU2_TYPE::GetDerivative(Alpha); Local_Sens_Temp = SU2_TYPE::GetDerivative(Temperature); diff --git a/SU2_CFD/src/solvers/CSolver.cpp b/SU2_CFD/src/solvers/CSolver.cpp index f2e6e4b6d43f..796e009f181b 100644 --- a/SU2_CFD/src/solvers/CSolver.cpp +++ b/SU2_CFD/src/solvers/CSolver.cpp @@ -4108,7 +4108,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/solvers/CTurbSASolver.cpp b/SU2_CFD/src/solvers/CTurbSASolver.cpp index 4a0d54a95176..3cece057cdc9 100644 --- a/SU2_CFD/src/solvers/CTurbSASolver.cpp +++ b/SU2_CFD/src/solvers/CTurbSASolver.cpp @@ -73,8 +73,10 @@ CTurbSASolver::CTurbSASolver(CGeometry *geometry, CConfig *config, unsigned shor LinSysRes.Initialize(nPoint, nPointDomain, nVar, 0.0); System.SetxIsZero(true); - if (ReducerStrategy) + if (ReducerStrategy) { EdgeFluxes.Initialize(geometry->GetnEdge(), geometry->GetnEdge(), nVar, nullptr); + EdgeFluxesDiff.Initialize(geometry->GetnEdge(), geometry->GetnEdge(), nVar, nullptr); + } if (config->GetExtraOutput()) { if (nDim == 2) { nOutputVariables = 13; } @@ -294,15 +296,77 @@ void CTurbSASolver::Postprocessing(CGeometry *geometry, CSolver **solver_contain void CTurbSASolver::Viscous_Residual(const unsigned long iEdge, const CGeometry* geometry, CSolver** solver_container, CNumerics* numerics, const CConfig* config) { - /*--- Define an object to set solver specific numerics contribution. ---*/ - auto SolverSpecificNumerics = [&](unsigned long iPoint, unsigned long jPoint) { + const bool implicit = (config->GetKind_TimeIntScheme() == EULER_IMPLICIT); + CFlowVariable* flowNodes = solver_container[FLOW_SOL] ? + su2staticcast_p(solver_container[FLOW_SOL]->GetNodes()) : nullptr; + + /*--- Points in edge. ---*/ + const auto iPoint = geometry->edges->GetNode(iEdge, 0); + const auto jPoint = geometry->edges->GetNode(iEdge, 1); + + /*--- Helper function to compute the flux ---*/ + auto ComputeFlux = [&](unsigned long iPoint, unsigned long jPoint, const su2double* normal) { + numerics->SetCoord(geometry->nodes->GetCoord(iPoint), geometry->nodes->GetCoord(jPoint)); + numerics->SetNormal(normal); + + if (flowNodes) { + numerics->SetPrimitive(flowNodes->GetPrimitive(iPoint), flowNodes->GetPrimitive(jPoint)); + } /*--- Roughness heights. ---*/ numerics->SetRoughness(geometry->nodes->GetRoughnessHeight(iPoint), geometry->nodes->GetRoughnessHeight(jPoint)); + + numerics->SetScalarVar(nodes->GetSolution(iPoint), nodes->GetSolution(jPoint)); + numerics->SetScalarVarGradient(nodes->GetGradient(iPoint), nodes->GetGradient(jPoint)); + + return numerics->ComputeResidual(config); }; - /*--- Now instantiate the generic implementation with the functor above. ---*/ + /*--- Compute fluxes and jacobians i->j ---*/ + const su2double* normal = geometry->edges->GetNormal(iEdge); + auto residual_ij = ComputeFlux(iPoint, jPoint, normal); - Viscous_Residual_impl(SolverSpecificNumerics, iEdge, geometry, solver_container, numerics, config); + JacobianScalarType *Block_ii = nullptr, *Block_ij = nullptr, *Block_ji = nullptr, *Block_jj = nullptr; + if (implicit) { + Jacobian.GetBlocks(iEdge, iPoint, jPoint, Block_ii, Block_ij, Block_ji, Block_jj); + } + if (ReducerStrategy) { + EdgeFluxes.SubtractBlock(iEdge, residual_ij); + EdgeFluxesDiff.SetBlock(iEdge, residual_ij); + if (implicit) { + /*--- For the reducer strategy the Jacobians are averaged for simplicity. ---*/ + assert(nVar == 1); + Block_ij[0] -= 0.5 * SU2_TYPE::GetValue(residual_ij.jacobian_j[0][0]); + Block_ji[0] += 0.5 * SU2_TYPE::GetValue(residual_ij.jacobian_i[0][0]); + } + } else { + LinSysRes.SubtractBlock(iPoint, residual_ij); + if (implicit) { + assert(nVar == 1); + Block_ii[0] -= SU2_TYPE::GetValue(residual_ij.jacobian_i[0][0]); + Block_ij[0] -= SU2_TYPE::GetValue(residual_ij.jacobian_j[0][0]); + } + } + + /*--- Compute fluxes and jacobians j->i ---*/ + su2double flipped_normal[MAXNDIM]; + for (auto iDim = 0u; iDim < nDim; iDim++) flipped_normal[iDim] = -normal[iDim]; + + auto residual_ji = ComputeFlux(jPoint, iPoint, flipped_normal); + if (ReducerStrategy) { + EdgeFluxesDiff.AddBlock(iEdge, residual_ji); + if (implicit) { + Block_ij[0] += 0.5 * SU2_TYPE::GetValue(residual_ji.jacobian_i[0][0]); + Block_ji[0] -= 0.5 * SU2_TYPE::GetValue(residual_ji.jacobian_j[0][0]); + } + } else { + LinSysRes.SubtractBlock(jPoint, residual_ji); + if (implicit) { + /*--- The order of arguments were flipped in the evaluation of residual_ji, the Jacobian + * associated with point i is stored in jacobian_j and point j in jacobian_i. ---*/ + Block_ji[0] -= SU2_TYPE::GetValue(residual_ji.jacobian_j[0][0]); + Block_jj[0] -= SU2_TYPE::GetValue(residual_ji.jacobian_i[0][0]); + } + } } void CTurbSASolver::Source_Residual(CGeometry *geometry, CSolver **solver_container, @@ -1560,7 +1624,7 @@ void CTurbSASolver::ComputeUnderRelaxationFactor(const CConfig *config) { system for this nonlinear iteration. */ su2double localUnderRelaxation = 1.00; - const su2double allowableRatio = 0.99; + const su2double allowableRatio = config->GetUnderRelax_SA(); SU2_OMP_FOR_STAT(omp_chunk_size) for (unsigned long iPoint = 0; iPoint < nPointDomain; iPoint++) { diff --git a/SU2_CFD/src/solvers/CTurbSSTSolver.cpp b/SU2_CFD/src/solvers/CTurbSSTSolver.cpp index 9fee88d18760..48b449e6e12d 100644 --- a/SU2_CFD/src/solvers/CTurbSSTSolver.cpp +++ b/SU2_CFD/src/solvers/CTurbSSTSolver.cpp @@ -75,8 +75,10 @@ CTurbSSTSolver::CTurbSSTSolver(CGeometry *geometry, CConfig *config, unsigned sh LinSysRes.Initialize(nPoint, nPointDomain, nVar, 0.0); System.SetxIsZero(true); - if (ReducerStrategy) + if (ReducerStrategy) { EdgeFluxes.Initialize(geometry->GetnEdge(), geometry->GetnEdge(), nVar, nullptr); + EdgeFluxesDiff.Initialize(geometry->GetnEdge(), geometry->GetnEdge(), nVar, nullptr); + } /*--- Initialize the BGS residuals in multizone problems. ---*/ if (multizone){ @@ -295,8 +297,100 @@ void CTurbSSTSolver::Viscous_Residual(const unsigned long iEdge, const CGeometry }; /*--- Now instantiate the generic implementation with the functor above. ---*/ + const bool implicit = (config->GetKind_TimeIntScheme() == EULER_IMPLICIT); + CFlowVariable* flowNodes = solver_container[FLOW_SOL] ? + su2staticcast_p(solver_container[FLOW_SOL]->GetNodes()) : nullptr; + + /*--- Points in edge ---*/ + auto iPoint = geometry->edges->GetNode(iEdge, 0); + auto jPoint = geometry->edges->GetNode(iEdge, 1); + + /*--- Helper function to compute the flux ---*/ + auto ComputeFlux = [&](unsigned long point_i, unsigned long point_j, const su2double* normal) { + numerics->SetCoord(geometry->nodes->GetCoord(point_i),geometry->nodes->GetCoord(point_j)); + numerics->SetNormal(normal); + + if (flowNodes) { + numerics->SetPrimitive(flowNodes->GetPrimitive(point_i), flowNodes->GetPrimitive(point_j)); + } + + numerics->SetScalarVar(nodes->GetSolution(point_i), nodes->GetSolution(point_j)); + numerics->SetScalarVarGradient(nodes->GetGradient(point_i), nodes->GetGradient(point_j)); + + return numerics->ComputeResidual(config); + }; + + SolverSpecificNumerics(iPoint, jPoint); + + /*--- Compute fluxes and jacobians i->j ---*/ + const su2double* normal = geometry->edges->GetNormal(iEdge); + auto residual_ij = ComputeFlux(iPoint, jPoint, normal); + + JacobianScalarType *Block_ii = nullptr, *Block_ij = nullptr, *Block_ji = nullptr, *Block_jj = nullptr; + if (implicit) { + Jacobian.GetBlocks(iEdge, iPoint, jPoint, Block_ii, Block_ij, Block_ji, Block_jj); + } + if (ReducerStrategy) { + EdgeFluxes.SubtractBlock(iEdge, residual_ij); + EdgeFluxesDiff.SetBlock(iEdge, residual_ij); + if (implicit) { + /*--- For the reducer strategy the Jacobians are averaged for simplicity. ---*/ + assert(nVar == 2); + for (int iVar=0; iVari ---*/ + su2double flipped_normal[MAXNDIM]; + for (auto iDim = 0u; iDim < nDim; iDim++) flipped_normal[iDim] = -normal[iDim]; - Viscous_Residual_impl(SolverSpecificNumerics, iEdge, geometry, solver_container, numerics, config); + auto residual_ji = ComputeFlux(jPoint, iPoint, flipped_normal); + if (ReducerStrategy) { + EdgeFluxesDiff.AddBlock(iEdge, residual_ji); + if (implicit) { + assert(nVar == 2); + for (int iVar=0; iVarGetDiffCoeff_kw(0); + DC_kw[1] = numerics->GetDiffCoeff_kw(1); + DC_kw[2] = numerics->GetDiffCoeff_kw(2); + DC_kw[3] = numerics->GetDiffCoeff_kw(3); + DC_kw[4] = numerics->GetDiffCoeff_kw(4); + DC_kw[5] = numerics->GetDiffCoeff_kw(5); + + nodes->SetDiffCoeff_kw(iPoint, DC_kw); } void CTurbSSTSolver::Source_Residual(CGeometry *geometry, CSolver **solver_container, @@ -1019,3 +1113,38 @@ void CTurbSSTSolver::SetUniformInlet(const CConfig* config, unsigned short iMark } } + +void CTurbSSTSolver::ComputeUnderRelaxationFactor(const CConfig *config) { + + const su2double allowableRatio = config->GetUnderRelax_SST(); + + SU2_OMP_FOR_STAT(omp_chunk_size) + /* Loop over the solution update given by relaxing the linear + system for this nonlinear iteration. */ + for (unsigned long iPoint = 0; iPoint < nPointDomain; iPoint++) { + su2double localUnderRelaxation = 1.0; + + for (unsigned short iVar =0; iVar < nVar; iVar++) { + const unsigned long index = iPoint * nVar + iVar; + su2double ratio = fabs(LinSysSol[index])/(fabs(nodes->GetSolution(iPoint, iVar)) + EPS); + /* We impose a limit on the maximum percentage that the + turbulence variables can change over a nonlinear iteration. */ + if (ratio > allowableRatio) { + localUnderRelaxation = min(allowableRatio / ratio, localUnderRelaxation); + } + + } + + /* Threshold the relaxation factor in the event that there is + a very small value. This helps avoid catastrophic crashes due + to non-realizable states by canceling the update. */ + + if (localUnderRelaxation < 1e-10) localUnderRelaxation = 0.0; + + /* Store the under-relaxation factor for this point. */ + + nodes->SetUnderRelaxation(iPoint, localUnderRelaxation); + } + + END_SU2_OMP_FOR +} \ No newline at end of file diff --git a/SU2_CFD/src/variables/CTurbVariable.cpp b/SU2_CFD/src/variables/CTurbVariable.cpp index 25a6a6801ef6..f0216ccb7ea5 100644 --- a/SU2_CFD/src/variables/CTurbVariable.cpp +++ b/SU2_CFD/src/variables/CTurbVariable.cpp @@ -35,6 +35,14 @@ CTurbVariable::CTurbVariable(unsigned long npoint, unsigned long ndim, unsigned turb_index.resize(nPoint) = su2double(1.0); intermittency.resize(nPoint) = su2double(1.0); + if (TurbModelFamily(config->GetKind_Turb_Model()) == TURB_FAMILY::KW) { + DC_TKE.resize(nPoint) = su2double(1.0); + DC_Omega.resize(nPoint) = su2double(1.0); + DC_OmegaT1.resize(nPoint) = su2double(1.0); + DC_OmegaT2.resize(nPoint) = su2double(1.0); + DC_OmegaT3.resize(nPoint) = su2double(1.0); + DC_Omega_kc.resize(nPoint) = su2double(1.0); + } } void CTurbVariable::RegisterEddyViscosity(bool input) { diff --git a/SU2_CFD/src/variables/CVariable.cpp b/SU2_CFD/src/variables/CVariable.cpp index 75bf88822826..3370fbb7a104 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/TestCases/TestCase.py b/TestCases/TestCase.py index 03b431ba3a7d..2da91f390462 100644 --- a/TestCases/TestCase.py +++ b/TestCases/TestCase.py @@ -46,6 +46,7 @@ def parse_args(description: str): parser = argparse.ArgumentParser(description=description) parser.add_argument('--tsan', action='store_true', help='Run thread sanitizer tests. Requires a tsan-enabled SU2 build.') parser.add_argument('--asan', action='store_true', help='Run address sanitizer tests. Requires an asan-enabled SU2 build.') + parser.add_argument('--tapetests', action='store_true', help='Run discrete adjoint tests in tape debug mode. Requires a SU2_CFD_AD build with Tag reverse type.') return parser.parse_args() class TestCase: @@ -107,10 +108,14 @@ def __init__(self,tag_in): # multizone problem self.multizone = False + # Indicate tapetest mode + self.enabled_with_tapetests = False + # The test condition. These must be set after initialization self.test_iter = 1 self.ntest_vals = 4 self.test_vals = [] + self.tapetest_vals = [] self.test_vals_aarch64 = [] self.cpu_arch = platform.machine().casefold() self.enabled_on_cpu_arch = ["x86_64","amd64","aarch64","arm64"] @@ -119,6 +124,7 @@ def __init__(self,tag_in): self.command = self.Command() self.timeout = 0 self.tol = 0.0 + self.tapetest_tol = 0 self.tol_file_percent = 0.0 self.comp_threshold = 0.0 @@ -127,9 +133,10 @@ def __init__(self,tag_in): self.reference_file_aarch64 = "" self.test_file = "of_grad.dat" - def run_test(self, with_tsan=False, with_asan=False): + def run_test(self, with_tsan=False, with_asan=False, with_tapetests=False): - if not self.is_enabled(with_tsan, with_asan): + # Check whether this test is valid and can be continued + if not self.is_enabled(with_tsan, with_asan, with_tapetests): return True print('==================== Start Test: %s ===================='%self.tag) @@ -138,6 +145,7 @@ def run_test(self, with_tsan=False, with_asan=False): timed_out = False iter_missing = True start_solver = True + tapetest_out = True # if root, add flag to mpirun self.command.allow_mpi_as_root() @@ -188,10 +196,10 @@ def run_test(self, with_tsan=False, with_asan=False): delta_vals = [] sim_vals = [] - if not with_tsan and not with_asan: # sanitizer findings result in non-zero return code, no need to examine the output + if not with_tsan and not with_asan and not with_tapetests: # Sanitizer findings result in non-zero return code, no need to examine the output. Tapetest output is examined separately. # Examine the output - f = open(logfilename,'r') - output = f.readlines() + with open(logfilename,'r') as f: + output = f.readlines() if not timed_out and len(self.test_vals) != 0: start_solver = False for line in output: @@ -240,6 +248,32 @@ def run_test(self, with_tsan=False, with_asan=False): #for j in output: # print(j) + if with_tapetests and self.enabled_with_tapetests: # examine the tapetest output + with open(logfilename,'r') as f: + output = f.readlines() + if not timed_out and len(self.tapetest_vals) != 0: + tapetest_out = False + for line in output: + if not tapetest_out: # Don't bother parsing anything before "Total number of tape inconsistencies" + if line.find('Total number of tape inconsistencies:') > -1: + tapetest_out = True + raw_data = line.split(':') # Split line into description and the string representing the number of errors + data = raw_data[1].strip() # Clear the string representing the number of errors (for now, expecting a single integer, but zone-wise error numbers are planned) + + if not len(self.tapetest_vals)==len(data): # something went wrong... probably bad input + print("Error in tapetest_vals!") + passed = False + break + for j in range(len(data)): + sim_vals.append( int(data[j]) ) + delta_vals.append( abs(int(data[j])-self.tapetest_vals[j]) ) + if delta_vals[j] > self.tapetest_tol: + exceed_tol = True + passed = False + + if not tapetest_out: + passed = False + process.communicate() if process.returncode != 0: @@ -259,15 +293,18 @@ def run_test(self, with_tsan=False, with_asan=False): if exceed_tol: print('ERROR: Difference between computed input and test_vals exceeded tolerance. TOL=%f'%self.tol) - if not start_solver: + if not start_solver and not with_tapetests: print('ERROR: The code was not able to get to the "Begin solver" section.') - if not with_tsan and not with_asan and iter_missing: + if not tapetest_out and with_tapetests: + print('ERROR: The code was not able to get to the "Total number of tape inconsistencies" line.') + + if not with_tsan and not with_asan and not with_tapetests and iter_missing: print('ERROR: The iteration number %d could not be found.'%self.test_iter) print('CPU architecture=%s' % self.cpu_arch) - if len(self.test_vals) != 0: + if len(self.test_vals) != 0 and not with_tapetests: print('test_iter=%d' % self.test_iter) print_vals(self.test_vals, name="test_vals (stored)") @@ -283,9 +320,9 @@ def run_test(self, with_tsan=False, with_asan=False): os.chdir(workdir) return passed - def run_filediff(self, with_tsan=False, with_asan=False): + def run_filediff(self, with_tsan=False, with_asan=False, with_tapetests=False): - if not self.is_enabled(with_tsan, with_asan): + if not self.is_enabled(with_tsan, with_asan, with_tapetests): return True print('==================== Start Test: %s ===================='%self.tag) @@ -480,9 +517,9 @@ def run_filediff(self, with_tsan=False, with_asan=False): os.chdir(workdir) return passed - def run_opt(self, with_tsan=False, with_asan=False): + def run_opt(self, with_tsan=False, with_asan=False, with_tapetests=False): - if not self.is_enabled(with_tsan, with_asan): + if not self.is_enabled(with_tsan, with_asan, with_tapetests): return True print('==================== Start Test: %s ===================='%self.tag) @@ -605,9 +642,9 @@ def run_opt(self, with_tsan=False, with_asan=False): os.chdir(workdir) return passed - def run_geo(self, with_tsan=False, with_asan=False): + def run_geo(self, with_tsan=False, with_asan=False, with_tapetests=False): - if not self.is_enabled(with_tsan, with_asan): + if not self.is_enabled(with_tsan, with_asan, with_tapetests): return True print('==================== Start Test: %s ===================='%self.tag) @@ -696,7 +733,7 @@ def run_geo(self, with_tsan=False, with_asan=False): delta_vals.append( abs(float(data[j])-self.test_vals[j]) ) if delta_vals[j] > self.tol: exceed_tol = True - passed = False + passed = False else: iter_missing = True @@ -745,9 +782,9 @@ def run_geo(self, with_tsan=False, with_asan=False): os.chdir(workdir) return passed - def run_def(self, with_tsan=False, with_asan=False): + def run_def(self, with_tsan=False, with_asan=False, with_tapetests=False): - if not self.is_enabled(with_tsan, with_asan): + if not self.is_enabled(with_tsan, with_asan, with_tapetests): return True print('==================== Start Test: %s ===================='%self.tag) @@ -970,19 +1007,24 @@ def disable_restart(self): return - def is_enabled(self, with_tsan=False, with_asan=False): + def is_enabled(self, with_tsan=False, with_asan=False, with_tapetests=False): is_enabled_on_arch = self.cpu_arch in self.enabled_on_cpu_arch if not is_enabled_on_arch: print('Ignoring test "%s" because it is not enabled for the current CPU architecture: %s' % (self.tag, self.cpu_arch)) + # A test case is valid to be continued if neither of the special modes (tsan, asan, tapetests) is active, or if so, the corresponding test case is enabled, too. tsan_compatible = not with_tsan or self.enabled_with_tsan asan_compatible = not with_asan or self.enabled_with_asan + tapetests_compatible = not with_tapetests or self.enabled_with_tapetests if not tsan_compatible: print('Ignoring test "%s" because it is not enabled to run with the thread sanitizer.' % self.tag) - return is_enabled_on_arch and tsan_compatible and asan_compatible + if not tapetests_compatible: + print('Ignoring test "%s" because it is not enabled to run a test of the tape.' % self.tag) + + return is_enabled_on_arch and tsan_compatible and asan_compatible and tapetests_compatible and tapetests_compatible def adjust_test_data(self): diff --git a/TestCases/cont_adj_euler/naca0012/inv_NACA0012_discadj_multizone.cfg b/TestCases/cont_adj_euler/naca0012/inv_NACA0012_discadj_multizone.cfg new file mode 100644 index 000000000000..85d1d03c5a24 --- /dev/null +++ b/TestCases/cont_adj_euler/naca0012/inv_NACA0012_discadj_multizone.cfg @@ -0,0 +1,112 @@ +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% % +% SU2 configuration file % +% Case description: Adjoint transonic inviscid flow around a NACA0012 airfoil % +% Author: Thomas D. Economon % +% Institution: Stanford University % +% Date: 2011.11.02 % +% File Version 8.2.0 "Harrier" % +% % +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +% ------------- DIRECT, ADJOINT, AND LINEARIZED PROBLEM DEFINITION ------------% +% +MULTIZONE= YES +MULTIZONE_ADAPT_FILENAME= NO +SOLVER= EULER +MATH_PROBLEM= DISCRETE_ADJOINT +RESTART_SOL= NO +READ_BINARY_RESTART= NO + +% ----------- COMPRESSIBLE AND INCOMPRESSIBLE FREE-STREAM DEFINITION ----------% +% +MACH_NUMBER= 0.8 +AOA= 1.25 +FREESTREAM_PRESSURE= 101325.0 +FREESTREAM_TEMPERATURE= 288.15 + +% ---------------------- REFERENCE VALUE DEFINITION ---------------------------% +% +REF_ORIGIN_MOMENT_X = 0.25 +REF_ORIGIN_MOMENT_Y = 0.00 +REF_ORIGIN_MOMENT_Z = 0.00 +REF_LENGTH= 1.0 +REF_AREA= 1.0 + +% ----------------------- BOUNDARY CONDITION DEFINITION -----------------------% +% +MARKER_EULER= ( airfoil ) +MARKER_FAR= ( farfield ) +MARKER_PLOTTING= ( airfoil ) +MARKER_MONITORING= ( airfoil ) + +% ------------- COMMON PARAMETERS TO DEFINE THE NUMERICAL METHOD --------------% +% +NUM_METHOD_GRAD= GREEN_GAUSS +OBJECTIVE_FUNCTION= DRAG +CFL_NUMBER= 5.0 +CFL_ADAPT= NO +CFL_ADAPT_PARAM= ( 1.5, 0.5, 1.0, 100.0 ) +RK_ALPHA_COEFF= ( 0.66667, 0.66667, 1.000000 ) +OUTER_ITER= 150 +INNER_ITER= 1 + +% -------------------- FLOW NUMERICAL METHOD DEFINITION -----------------------% +% +CONV_NUM_METHOD_FLOW= JST +JST_SENSOR_COEFF= ( 0.5, 0.02 ) +TIME_DISCRE_FLOW= EULER_IMPLICIT + +% ---------------- ADJOINT-FLOW NUMERICAL METHOD DEFINITION -------------------% +% +CONV_NUM_METHOD_ADJFLOW= ROE +MUSCL_ADJFLOW= YES +SLOPE_LIMITER_ADJFLOW= NONE +ADJ_SHARP_LIMITER_COEFF= 3.0 +ADJ_JST_SENSOR_COEFF= ( 0.5, 0.02 ) +CFL_REDUCTION_ADJFLOW= 0.5 +TIME_DISCRE_ADJFLOW= EULER_IMPLICIT + +% ----------------------- DESIGN VARIABLE PARAMETERS --------------------------% +% +DV_KIND= HICKS_HENNE +DV_MARKER= ( airfoil ) +DV_PARAM= ( 1, 0.5 ) +DV_VALUE= 0.01 + +% ------------------------ GRID DEFORMATION PARAMETERS ------------------------% +% +DEFORM_STIFFNESS_TYPE= WALL_DISTANCE + +% --------------------------- CONVERGENCE PARAMETERS --------------------------% +% +CONV_RESIDUAL_MINVAL= -12 +CONV_STARTITER= 10 +CONV_CAUCHY_ELEMS= 100 +CONV_CAUCHY_EPS= 1E-6 + +% ------------------------- INPUT/OUTPUT INFORMATION --------------------------% +% +MESH_FILENAME= mesh_NACA0012_inv.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= 250 +SCREEN_OUTPUT= (OUTER_ITER, BGS_RES[0]) +OUTPUT_FILES=(RESTART_ASCII) + +% --------------------- OPTIMAL SHAPE DESIGN DEFINITION -----------------------% +% +OPT_OBJECTIVE= DRAG * 0.001 +OPT_CONSTRAINT= ( LIFT > 0.327 ) * 0.001; ( MOMENT_Z > 0.0 ) * 0.001; ( AIRFOIL_THICKNESS > 0.12 ) * 0.001 +DEFINITION_DV= ( 30, 1.0 | airfoil | 0, 0.05 ); ( 30, 1.0 | airfoil | 0, 0.10 ); ( 30, 1.0 | airfoil | 0, 0.15 ); ( 30, 1.0 | airfoil | 0, 0.20 ); ( 30, 1.0 | airfoil | 0, 0.25 ); ( 30, 1.0 | airfoil | 0, 0.30 ); ( 30, 1.0 | airfoil | 0, 0.35 ); ( 30, 1.0 | airfoil | 0, 0.40 ); ( 30, 1.0 | airfoil | 0, 0.45 ); ( 30, 1.0 | airfoil | 0, 0.50 ); ( 30, 1.0 | airfoil | 0, 0.55 ); ( 30, 1.0 | airfoil | 0, 0.60 ); ( 30, 1.0 | airfoil | 0, 0.65 ); ( 30, 1.0 | airfoil | 0, 0.70 ); ( 30, 1.0 | airfoil | 0, 0.75 ); ( 30, 1.0 | airfoil | 0, 0.80 ); ( 30, 1.0 | airfoil | 0, 0.85 ); ( 30, 1.0 | airfoil | 0, 0.90 ); ( 30, 1.0 | airfoil | 0, 0.95 ); ( 30, 1.0 | airfoil | 1, 0.05 ); ( 30, 1.0 | airfoil | 1, 0.10 ); ( 30, 1.0 | airfoil | 1, 0.15 ); ( 30, 1.0 | airfoil | 1, 0.20 ); ( 30, 1.0 | airfoil | 1, 0.25 ); ( 30, 1.0 | airfoil | 1, 0.30 ); ( 30, 1.0 | airfoil | 1, 0.35 ); ( 30, 1.0 | airfoil | 1, 0.40 ); ( 30, 1.0 | airfoil | 1, 0.45 ); ( 30, 1.0 | airfoil | 1, 0.50 ); ( 30, 1.0 | airfoil | 1, 0.55 ); ( 30, 1.0 | airfoil | 1, 0.60 ); ( 30, 1.0 | airfoil | 1, 0.65 ); ( 30, 1.0 | airfoil | 1, 0.70 ); ( 30, 1.0 | airfoil | 1, 0.75 ); ( 30, 1.0 | airfoil | 1, 0.80 ); ( 30, 1.0 | airfoil | 1, 0.85 ); ( 30, 1.0 | airfoil | 1, 0.90 ); ( 30, 1.0 | airfoil | 1, 0.95 ) diff --git a/TestCases/hybrid_regression.py b/TestCases/hybrid_regression.py index 366589a7994f..e4a747f2c913 100644 --- a/TestCases/hybrid_regression.py +++ b/TestCases/hybrid_regression.py @@ -181,7 +181,7 @@ def main(): turb_flatplate.cfg_dir = "rans/flatplate" turb_flatplate.cfg_file = "turb_SA_flatplate.cfg" turb_flatplate.test_iter = 20 - turb_flatplate.test_vals = [-4.316134, -6.737979, -0.187461, 0.057468] + turb_flatplate.test_vals = [-4.316127, -6.738720, -0.187461, 0.057469] test_list.append(turb_flatplate) # ONERA M6 Wing @@ -189,7 +189,7 @@ def main(): turb_oneram6.cfg_dir = "rans/oneram6" turb_oneram6.cfg_file = "turb_ONERAM6.cfg" turb_oneram6.test_iter = 10 - turb_oneram6.test_vals = [-2.408675, -6.662904, 0.238578, 0.158968, 0.000000] + turb_oneram6.test_vals = [-2.408655, -6.628338, 0.238580, 0.158951, 0.000000] test_list.append(turb_oneram6) # NACA0012 (SA, FUN3D finest grid results: CL=1.0983, CD=0.01242) @@ -197,8 +197,8 @@ def main(): turb_naca0012_sa.cfg_dir = "rans/naca0012" turb_naca0012_sa.cfg_file = "turb_NACA0012_sa.cfg" turb_naca0012_sa.test_iter = 5 - turb_naca0012_sa.test_vals = [-12.098325, -14.149988, 1.057665, 0.022971, 20.000000, -2.292707, 0.000000, -12.068169, 0] - turb_naca0012_sa.test_vals_aarch64 = [-12.098325, -14.149988, 1.057665, 0.022971, 20.000000, -2.292707, 0.000000, -12.068169, 0] + turb_naca0012_sa.test_vals = [-12.050637, -16.149098, 1.058588, 0.022984, 20.000000, -2.832819, 0.000000, -14.067279, 0] + turb_naca0012_sa.test_vals_aarch64 = [-12.050637, -16.149098, 1.058588, 0.022984, 20.000000, -2.832819, 0.000000, -14.067279, 0] test_list.append(turb_naca0012_sa) # NACA0012 (SST, FUN3D finest grid results: CL=1.0840, CD=0.01253) @@ -238,7 +238,7 @@ def main(): propeller.cfg_dir = "rans/propeller" propeller.cfg_file = "propeller.cfg" propeller.test_iter = 10 - propeller.test_vals = [-3.389724, -8.409502, 0.000048, 0.056344] + propeller.test_vals = [-3.389724, -8.410479, 0.000048, 0.056344] test_list.append(propeller) ####################################### @@ -332,7 +332,7 @@ def main(): hb_rans_preconditioning.cfg_dir = "harmonic_balance/hb_rans_preconditioning" hb_rans_preconditioning.cfg_file = "davis.cfg" hb_rans_preconditioning.test_iter = 25 - hb_rans_preconditioning.test_vals = [-1.902111, 0.484080, 0.601469, 3.608991, -5.949369] + hb_rans_preconditioning.test_vals = [-1.902098, 0.484244, 0.601482, 3.609005, -5.943887] test_list.append(hb_rans_preconditioning) ############################# @@ -403,7 +403,7 @@ def main(): inc_turb_naca0012.cfg_dir = "incomp_rans/naca0012" inc_turb_naca0012.cfg_file = "naca0012.cfg" inc_turb_naca0012.test_iter = 20 - inc_turb_naca0012.test_vals = [-4.788405, -11.040560, 0.000008, 0.309505] + inc_turb_naca0012.test_vals = [-4.788405, -11.040877, 0.000008, 0.309505] test_list.append(inc_turb_naca0012) # NACA0012, SST_SUST @@ -419,7 +419,7 @@ def main(): inc_weakly_coupled.cfg_dir = "disc_adj_heat" inc_weakly_coupled.cfg_file = "primal.cfg" inc_weakly_coupled.test_iter = 10 - inc_weakly_coupled.test_vals = [-18.894811, -17.879327, -18.412938, -17.855948, -18.343462, -15.659612, 5.545700] + inc_weakly_coupled.test_vals = [-18.095922, -16.331787, -16.514014, -13.703679, -18.203865, -14.053738, 5.545900] test_list.append(inc_weakly_coupled) ###################################### @@ -503,7 +503,7 @@ def main(): ddes_flatplate.cfg_dir = "ddes/flatplate" ddes_flatplate.cfg_file = "ddes_flatplate.cfg" ddes_flatplate.test_iter = 10 - ddes_flatplate.test_vals = [-2.714713, -5.788302, -0.214960, 0.023758, 0.000000] + ddes_flatplate.test_vals = [-2.714713, -5.763293, -0.214960, 0.023758, 0.000000] ddes_flatplate.unsteady = True test_list.append(ddes_flatplate) @@ -582,8 +582,8 @@ def main(): 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] + multi_interface.test_vals = [-8.632229, -8.894737, -9.348730] + multi_interface.test_vals_aarch64 = [-8.632229, -8.894737, -9.348730] test_list.append(multi_interface) ###################################### diff --git a/TestCases/hybrid_regression_AD.py b/TestCases/hybrid_regression_AD.py index 0727b5888735..9b4bdb23c52f 100644 --- a/TestCases/hybrid_regression_AD.py +++ b/TestCases/hybrid_regression_AD.py @@ -78,7 +78,7 @@ def main(): discadj_rans_naca0012_sa.cfg_dir = "disc_adj_rans/naca0012" discadj_rans_naca0012_sa.cfg_file = "turb_NACA0012_sa.cfg" discadj_rans_naca0012_sa.test_iter = 10 - discadj_rans_naca0012_sa.test_vals = [-2.997064, -0.196172, 0.000003, -0.000000, 5.000000, -2.919675, 5.000000, -7.323218] + discadj_rans_naca0012_sa.test_vals = [-2.997050, -0.199287, 0.000003, -0.000000, 5.000000, -2.919668, 5.000000, -7.320138] test_list.append(discadj_rans_naca0012_sa) # Adjoint turbulent NACA0012 SST @@ -123,7 +123,7 @@ def main(): discadj_incomp_turb_NACA0012_sa.cfg_dir = "disc_adj_incomp_rans/naca0012" discadj_incomp_turb_NACA0012_sa.cfg_file = "turb_naca0012_sa.cfg" discadj_incomp_turb_NACA0012_sa.test_iter = 10 - discadj_incomp_turb_NACA0012_sa.test_vals = [10.000000, -3.845995, -1.031096, 0.000000] + discadj_incomp_turb_NACA0012_sa.test_vals = [10.000000, -3.845995, -1.023534, 0.000000] test_list.append(discadj_incomp_turb_NACA0012_sa) # Adjoint Incompressible Turbulent NACA 0012 SST diff --git a/TestCases/mms/fvm_euler/inv_mms_vortex.cfg b/TestCases/mms/fvm_euler/inv_mms_vortex.cfg new file mode 100644 index 000000000000..c069ba281de7 --- /dev/null +++ b/TestCases/mms/fvm_euler/inv_mms_vortex.cfg @@ -0,0 +1,89 @@ +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% % +% SU2 configuration file % +% Case description: Isentropic vortex (FVM) % +% Author: Brian Munguía % +% Date: 2025.25.06 % +% File Version 8.2.0 "Harrier" % +% % +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +% ------------- DIRECT, ADJOINT, AND LINEARIZED PROBLEM DEFINITION ------------% +% +KIND_VERIFICATION_SOLUTION= INVISCID_VORTEX +SOLVER= EULER +MATH_PROBLEM= DIRECT +RESTART_SOL= NO + +% ----------- COMPRESSIBLE AND INCOMPRESSIBLE FREE-STREAM DEFINITION ----------% +% +MACH_NUMBER= 0.5 +AOA= 45.0 +FREESTREAM_DENSITY= 1.0 +FREESTREAM_PRESSURE= 1.0 +FREESTREAM_TEMPERATURE= 1.0 + +% -------------- COMPRESSIBLE AND INCOMPRESSIBLE FLUID CONSTANTS --------------% +% +FLUID_MODEL= IDEAL_GAS +GAMMA_VALUE= 1.4 +GAS_CONSTANT= 1.0 + +% ---------------------- REFERENCE VALUE DEFINITION ---------------------------% +% +REF_ORIGIN_MOMENT_X= 0.00 +REF_ORIGIN_MOMENT_Y= 0.00 +REF_ORIGIN_MOMENT_Z= 0.00 +REF_AREA= 1.0 +REF_DIMENSIONALIZATION= DIMENSIONAL + +% ------------------------- UNSTEADY SIMULATION -------------------------------% +% +TIME_DOMAIN= YES +TIME_MARCHING= TIME_STEPPING +TIME_STEP= 2.0e-3 +MAX_TIME= 50.0 + +% -------------------- BOUNDARY CONDITION DEFINITION --------------------------% +% +MARKER_PERIODIC= ( PeriodicBottom, PeriodicTop, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 2.0, 0.0, \ + PeriodicLeft, PeriodicRight, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 2.0, 0.0, 0.0 ) + +% ------------- COMMON PARAMETERS TO DEFINE THE NUMERICAL METHOD --------------% +% +CFL_NUMBER= 0.2 +TIME_ITER= 50 +RK_ALPHA_COEFF= ( 0.666667, 0.666667, 1.0 ) + +% ------------------ FLOW NUMERICAL METHOD DEFINITION ----------------------% +% +CONV_NUM_METHOD_FLOW= ROE +MUSCL_FLOW= YES +SLOPE_LIMITER_FLOW= NONE +TIME_DISCRE_FLOW= EULER_IMPLICIT + +% --------------------------- CONVERGENCE PARAMETERS --------------------------% +% +CONV_RESIDUAL_MINVAL= -15 +CONV_STARTITER= 10 +CONV_CAUCHY_ELEMS= 100 +CONV_CAUCHY_EPS= 1E-6 + +% ------------------------- INPUT/OUTPUT INFORMATION --------------------------% +% +MESH_FILENAME= TriAdapt.su2 +MESH_FORMAT= SU2 +MESH_OUT_FILENAME= mesh_out.su2 +SOLUTION_FILENAME= solution_flow.dat +SOLUTION_ADJ_FILENAME= solution_adj.dat +TABULAR_FORMAT= CSV +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= 50 +SCREEN_OUTPUT= (TIME_ITER, RMS_DENSITY, RMS_ENERGY, LIFT, DRAG) diff --git a/TestCases/parallel_regression.py b/TestCases/parallel_regression.py index a6c8cf296637..fe18c70ed446 100644 --- a/TestCases/parallel_regression.py +++ b/TestCases/parallel_regression.py @@ -357,7 +357,7 @@ def main(): rae2822_sa.cfg_dir = "rans/rae2822" rae2822_sa.cfg_file = "turb_SA_RAE2822.cfg" rae2822_sa.test_iter = 20 - rae2822_sa.test_vals = [-2.004689, -5.265797, 0.809463, 0.062016, 0] + rae2822_sa.test_vals = [-2.004689, -5.265730, 0.809462, 0.062011, 0.000000] test_list.append(rae2822_sa) # RAE2822 SST @@ -381,7 +381,7 @@ def main(): turb_flatplate.cfg_dir = "rans/flatplate" turb_flatplate.cfg_file = "turb_SA_flatplate.cfg" turb_flatplate.test_iter = 20 - turb_flatplate.test_vals = [-4.297198, -6.730442, -0.187632, 0.057700] + turb_flatplate.test_vals = [-4.297192, -6.731227, -0.187632, 0.057700] test_list.append(turb_flatplate) # Flat plate (compressible) with species inlet @@ -389,7 +389,7 @@ def main(): turb_flatplate_species.cfg_dir = "rans/flatplate" turb_flatplate_species.cfg_file = "turb_SA_flatplate_species.cfg" turb_flatplate_species.test_iter = 20 - turb_flatplate_species.test_vals = [-4.249479, -0.634915, -1.716291, 1.223196, -3.307930, 9.000000, -6.634088, 5.000000, -6.985954, 10.000000, -6.255640, 0.999903, 0.999903] + turb_flatplate_species.test_vals = [-4.249474, -0.634908, -1.716288, 1.223201, -3.307930, 9.000000, -6.634095, 5.000000, -6.986784, 10.000000, -6.255641, 0.999903, 0.999903] test_list.append(turb_flatplate_species) # Flat plate SST compressibility correction Wilcox @@ -413,7 +413,7 @@ def main(): turb_oneram6.cfg_dir = "rans/oneram6" turb_oneram6.cfg_file = "turb_ONERAM6.cfg" turb_oneram6.test_iter = 10 - turb_oneram6.test_vals = [-2.408685, -6.662907, 0.238579, 0.158968, 0.000000] + turb_oneram6.test_vals = [-2.408664, -6.628340, 0.238581, 0.158952, 0.000000] turb_oneram6.timeout = 3200 test_list.append(turb_oneram6) @@ -422,7 +422,7 @@ def main(): turb_oneram6_vc.cfg_dir = "rans/oneram6" turb_oneram6_vc.cfg_file = "turb_ONERAM6_vc.cfg" turb_oneram6_vc.test_iter = 15 - turb_oneram6_vc.test_vals = [-2.282318, -6.614780, 0.234330, 0.143024, 0.000000] + turb_oneram6_vc.test_vals = [-2.282278, -6.568458, 0.234350, 0.142989, 0.000000] turb_oneram6_vc.timeout = 3200 test_list.append(turb_oneram6_vc) @@ -431,7 +431,7 @@ def main(): turb_oneram6_nk.cfg_dir = "rans/oneram6" turb_oneram6_nk.cfg_file = "turb_ONERAM6_nk.cfg" turb_oneram6_nk.test_iter = 20 - turb_oneram6_nk.test_vals = [-4.828066, -4.426250, -11.417591, 0.224679, 0.044309, 1.000000, -0.642977, 31.384000] + turb_oneram6_nk.test_vals = [-4.827571, -4.425650, -11.379658, 0.224787, 0.044208, 1.000000, -0.642711, 31.384000] turb_oneram6_nk.timeout = 600 turb_oneram6_nk.tol = 0.0001 test_list.append(turb_oneram6_nk) @@ -441,7 +441,7 @@ def main(): turb_naca0012_sa.cfg_dir = "rans/naca0012" turb_naca0012_sa.cfg_file = "turb_NACA0012_sa.cfg" turb_naca0012_sa.test_iter = 5 - turb_naca0012_sa.test_vals = [-12.094695, -14.685268, 1.057665, 0.022971, 20.000000, -1.692967, 20.000000, -4.037673, 0] + turb_naca0012_sa.test_vals = [-12.050143, -16.185204, 1.058588, 0.022984, 20.000000, -1.562786, 20.000000, -3.900124, 0.000000] turb_naca0012_sa.timeout = 3200 test_list.append(turb_naca0012_sa) @@ -507,7 +507,7 @@ def main(): propeller.cfg_dir = "rans/propeller" propeller.cfg_file = "propeller.cfg" propeller.test_iter = 10 - propeller.test_vals = [-3.389724, -8.409502, 0.000048, 0.056344] + propeller.test_vals = [-3.389724, -8.410479, 0.000048, 0.056344] propeller.timeout = 3200 test_list.append(propeller) @@ -646,7 +646,7 @@ def main(): inc_turb_naca0012.cfg_dir = "incomp_rans/naca0012" inc_turb_naca0012.cfg_file = "naca0012.cfg" inc_turb_naca0012.test_iter = 20 - inc_turb_naca0012.test_vals = [-4.788595, -11.040625, -0.000002, 0.309519] + inc_turb_naca0012.test_vals = [-4.788595, -11.040942, -0.000002, 0.309519] test_list.append(inc_turb_naca0012) # NACA0012, SST_SUST @@ -732,7 +732,7 @@ def main(): turbmod_sa_bsl_rae2822.cfg_dir = "turbulence_models/sa/rae2822" turbmod_sa_bsl_rae2822.cfg_file = "turb_SA_BSL_RAE2822.cfg" turbmod_sa_bsl_rae2822.test_iter = 20 - turbmod_sa_bsl_rae2822.test_vals = [-2.004689, 0.742306, 0.497308, -5.265797, 0.809463, 0.062016] + turbmod_sa_bsl_rae2822.test_vals = [-2.004689, 0.742307, 0.497308, -5.265730, 0.809462, 0.062011] test_list.append(turbmod_sa_bsl_rae2822) # SA Negative @@ -740,7 +740,7 @@ def main(): turbmod_sa_neg_rae2822.cfg_dir = "turbulence_models/sa/rae2822" turbmod_sa_neg_rae2822.cfg_file = "turb_SA_NEG_RAE2822.cfg" turbmod_sa_neg_rae2822.test_iter = 10 - turbmod_sa_neg_rae2822.test_vals = [-1.204800, 1.611685, 1.349330, 1.489602, 1.263603, 0.466487, 0] + turbmod_sa_neg_rae2822.test_vals = [-1.345491, 1.448526, 1.208572, -0.846542, 1.263765, 0.494751, 0.000000] turbmod_sa_neg_rae2822.test_vals_aarch64 = [-1.359612, 1.493629, 1.218367, -1.441703, 1.248499, 0.457987, 0] test_list.append(turbmod_sa_neg_rae2822) @@ -749,7 +749,7 @@ def main(): turbmod_sa_comp_rae2822.cfg_dir = "turbulence_models/sa/rae2822" turbmod_sa_comp_rae2822.cfg_file = "turb_SA_COMP_RAE2822.cfg" turbmod_sa_comp_rae2822.test_iter = 20 - turbmod_sa_comp_rae2822.test_vals = [-2.004687, 0.742304, 0.497309, -5.266084, 0.809467, 0.062029] + turbmod_sa_comp_rae2822.test_vals = [-2.004688, 0.742305, 0.497309, -5.266066, 0.809466, 0.062026] test_list.append(turbmod_sa_comp_rae2822) # SA Edwards @@ -757,7 +757,7 @@ def main(): turbmod_sa_edw_rae2822.cfg_dir = "turbulence_models/sa/rae2822" turbmod_sa_edw_rae2822.cfg_file = "turb_SA_EDW_RAE2822.cfg" turbmod_sa_edw_rae2822.test_iter = 20 - turbmod_sa_edw_rae2822.test_vals = [-2.004687, 0.742306, 0.497310, -5.290769, 0.809485, 0.062036] + turbmod_sa_edw_rae2822.test_vals = [-2.004687, 0.742306, 0.497310, -5.290787, 0.809485, 0.062034] test_list.append(turbmod_sa_edw_rae2822) # SA Compressibility and Edwards @@ -765,7 +765,7 @@ def main(): turbmod_sa_comp_edw_rae2822.cfg_dir = "turbulence_models/sa/rae2822" turbmod_sa_comp_edw_rae2822.cfg_file = "turb_SA_COMP_EDW_RAE2822.cfg" turbmod_sa_comp_edw_rae2822.test_iter = 20 - turbmod_sa_comp_edw_rae2822.test_vals = [-2.004685, 0.742307, 0.497311, -5.290750, 0.809487, 0.062045] + turbmod_sa_comp_edw_rae2822.test_vals = [-2.004686, 0.742307, 0.497311, -5.290776, 0.809487, 0.062044] test_list.append(turbmod_sa_comp_edw_rae2822) # SA QCR @@ -773,7 +773,7 @@ def main(): turbmod_sa_qcr_rae2822.cfg_dir = "turbulence_models/sa/rae2822" turbmod_sa_qcr_rae2822.cfg_file = "turb_SA_QCR_RAE2822.cfg" turbmod_sa_qcr_rae2822.test_iter = 20 - turbmod_sa_qcr_rae2822.test_vals = [-2.004793, 0.742353, 0.497315, -5.265977, 0.807841, 0.062027] + turbmod_sa_qcr_rae2822.test_vals = [-2.004794, 0.742354, 0.497315, -5.265910, 0.807835, 0.062022] test_list.append(turbmod_sa_qcr_rae2822) ############################ @@ -785,7 +785,7 @@ def main(): schubauer_klebanoff_transition.cfg_dir = "transition/Schubauer_Klebanoff" schubauer_klebanoff_transition.cfg_file = "transitional_BC_model_ConfigFile.cfg" schubauer_klebanoff_transition.test_iter = 10 - schubauer_klebanoff_transition.test_vals = [-8.215651, -13.240283, 0.000048, 0.007983] + schubauer_klebanoff_transition.test_vals = [-8.215651, -13.239431, 0.000048, 0.007983] test_list.append(schubauer_klebanoff_transition) ##################################### @@ -942,7 +942,7 @@ def main(): hb_rans_preconditioning.cfg_file = "davis.cfg" hb_rans_preconditioning.test_iter = 25 hb_rans_preconditioning.tol = 0.00001 - hb_rans_preconditioning.test_vals = [-1.902098, 0.484070, 0.601481, 3.609002, -5.949356] + hb_rans_preconditioning.test_vals = [-1.902085, 0.484234, 0.601494, 3.609016, -5.943874] test_list.append(hb_rans_preconditioning) ###################################### @@ -1009,7 +1009,7 @@ def main(): ddes_flatplate.cfg_dir = "ddes/flatplate" ddes_flatplate.cfg_file = "ddes_flatplate.cfg" ddes_flatplate.test_iter = 10 - ddes_flatplate.test_vals = [-2.714713, -5.788290, -0.214960, 0.023758, 0.000000] + ddes_flatplate.test_vals = [-2.714713, -5.763284, -0.214960, 0.023758, 0.000000] ddes_flatplate.unsteady = True test_list.append(ddes_flatplate) @@ -1086,7 +1086,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 = [-7.701448, -8.512241, -6.014939, -6.468738, -5.801759, -4.607179, -5.551037, -5.300771, -3.804188, -5.256055, -5.765160, -3.609605, -2.229276, -2.883962, -0.563469] + Aachen_3D_restart.test_vals = [-7.701447, -8.512235, -6.014939, -6.468414, -5.801735, -4.607179, -5.550688, -5.300764, -3.804188, -5.256009, -5.764984, -3.609605, -2.229276, -2.883960, -0.563469] test_list.append(Aachen_3D_restart) # Jones APU Turbocharger restart @@ -1119,8 +1119,8 @@ def main(): 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] + multi_interface.test_vals = [-8.632229, -8.894737, -9.348730] + multi_interface.test_vals_aarch64 = [-8.632229, -8.894737, -9.348730] test_list.append(multi_interface) ###################################### @@ -1479,6 +1479,15 @@ def main(): ### Method of Manufactured Solutions (MMS) ### ############################################## + # FVM, compressible, euler (periodic isentropic vortex) + mms_fvm_vortex = TestCase('mms_fvm_vortex') + mms_fvm_vortex.cfg_dir = "mms/fvm_euler" + mms_fvm_vortex.cfg_file = "inv_mms_vortex.cfg" + mms_fvm_vortex.test_iter = 10 + mms_fvm_vortex.test_vals = [-5.704300, -4.848072, 0.000000, 0.000000] + mms_fvm_vortex.unsteady = True + test_list.append(mms_fvm_vortex) + # FVM, compressible, laminar N-S mms_fvm_ns = TestCase('mms_fvm_ns') mms_fvm_ns.cfg_dir = "mms/fvm_navierstokes" @@ -1558,7 +1567,7 @@ def main(): species2_primitiveVenturi_mixingmodel_viscosity.cfg_dir = "species_transport/venturi_primitive_3species" species2_primitiveVenturi_mixingmodel_viscosity.cfg_file = "species2_primitiveVenturi_mixingmodel_viscosity.cfg" species2_primitiveVenturi_mixingmodel_viscosity.test_iter = 50 - species2_primitiveVenturi_mixingmodel_viscosity.test_vals = [-5.232339, -3.617118, -3.857221, -7.533847, -5.126646, 5.000000, -1.682962, 5.000000, -3.474078, 5.000000, -2.086859, 2.495548, 0.985490, 0.600234, 0.909824] + species2_primitiveVenturi_mixingmodel_viscosity.test_vals = [-5.236181, -3.624224, -3.871405, -7.533947, -5.128207, 5.000000, -1.678081, 5.000000, -3.489246, 5.000000, -2.086184, 2.495408, 0.985469, 0.600247, 0.909691] test_list.append(species2_primitiveVenturi_mixingmodel_viscosity) # 2 species (1 eq) primitive venturi mixing using mixing model including heat capacity and mass diffusivity diff --git a/TestCases/parallel_regression_AD.py b/TestCases/parallel_regression_AD.py index 21ec615f9ff4..ce5ed630cc01 100644 --- a/TestCases/parallel_regression_AD.py +++ b/TestCases/parallel_regression_AD.py @@ -83,7 +83,7 @@ def main(): discadj_rans_naca0012_sa.cfg_dir = "disc_adj_rans/naca0012" discadj_rans_naca0012_sa.cfg_file = "turb_NACA0012_sa.cfg" discadj_rans_naca0012_sa.test_iter = 10 - discadj_rans_naca0012_sa.test_vals = [-2.996963, -0.196020, 0.000004, -0.000000, 5.000000, -3.430615, 5.000000, -7.411381] + discadj_rans_naca0012_sa.test_vals = [-2.996948, -0.199135, 0.000004, -0.000000, 5.000000, -3.430662, 5.000000, -7.402623] test_list.append(discadj_rans_naca0012_sa) # Adjoint turbulent NACA0012 SST @@ -128,8 +128,8 @@ def main(): discadj_incomp_turb_NACA0012_sa.cfg_dir = "disc_adj_incomp_rans/naca0012" discadj_incomp_turb_NACA0012_sa.cfg_file = "turb_naca0012_sa.cfg" discadj_incomp_turb_NACA0012_sa.test_iter = 10 - discadj_incomp_turb_NACA0012_sa.test_vals = [10.000000, -3.846020, -1.031079, 0.000000] - discadj_incomp_turb_NACA0012_sa.test_vals_aarch64 = [10.000000, -3.846020, -1.031078, 0.000000] + discadj_incomp_turb_NACA0012_sa.test_vals = [10.000000, -3.846020, -1.023517, 0.000000] + discadj_incomp_turb_NACA0012_sa.test_vals_aarch64 = [10.000000, -3.846020, -1.023517, 0.000000] test_list.append(discadj_incomp_turb_NACA0012_sa) # Adjoint Incompressible Turbulent NACA 0012 SST @@ -256,7 +256,7 @@ def main(): discadj_heat.cfg_dir = "disc_adj_heat" discadj_heat.cfg_file = "disc_adj_heat.cfg" discadj_heat.test_iter = 10 - discadj_heat.test_vals = [-1.999668, 0.670563, 0.000000, 0.006210] + discadj_heat.test_vals = [-1.999669, 0.657451, 0.000000, 0.006210] discadj_heat.test_vals_aarch64 = [-2.226539, 0.605868, 0.000000, -6.256400] test_list.append(discadj_heat) diff --git a/TestCases/serial_regression.py b/TestCases/serial_regression.py index 8cb540e0519f..aba223dcefde 100644 --- a/TestCases/serial_regression.py +++ b/TestCases/serial_regression.py @@ -218,7 +218,7 @@ def main(): rae2822_sa.cfg_dir = "rans/rae2822" rae2822_sa.cfg_file = "turb_SA_RAE2822.cfg" rae2822_sa.test_iter = 20 - rae2822_sa.test_vals = [-2.020123, -5.269339, 0.807147, 0.060499, 0] + rae2822_sa.test_vals = [-2.020123, -5.269264, 0.807147, 0.060494, 0.000000] test_list.append(rae2822_sa) # RAE2822 SST @@ -242,7 +242,7 @@ def main(): turb_flatplate.cfg_dir = "rans/flatplate" turb_flatplate.cfg_file = "turb_SA_flatplate.cfg" turb_flatplate.test_iter = 20 - turb_flatplate.test_vals = [-4.316135, -6.737979, -0.187461, 0.057468] + turb_flatplate.test_vals = [-4.316128, -6.738720, -0.187461, 0.057469] test_list.append(turb_flatplate) # FLAT PLATE, WALL FUNCTIONS, COMPRESSIBLE SST @@ -258,7 +258,7 @@ def main(): turb_wallfunction_flatplate_sa.cfg_dir = "wallfunctions/flatplate/compressible_SA" turb_wallfunction_flatplate_sa.cfg_file = "turb_SA_flatplate.cfg" turb_wallfunction_flatplate_sa.test_iter = 10 - turb_wallfunction_flatplate_sa.test_vals = [-4.487562, -2.016144, -2.169439, 0.834724, -5.382532, 10.000000, -1.508186, 0.034484, 0.002639] + turb_wallfunction_flatplate_sa.test_vals = [-4.456430, -1.959113, -2.105871, 0.883035, -5.248648, 10.000000, -1.537437, 0.034538, 0.002637] test_list.append(turb_wallfunction_flatplate_sa) # ONERA M6 Wing @@ -266,7 +266,7 @@ def main(): turb_oneram6.cfg_dir = "rans/oneram6" turb_oneram6.cfg_file = "turb_ONERAM6.cfg" turb_oneram6.test_iter = 10 - turb_oneram6.test_vals = [-2.408685, -6.662908, 0.238580, 0.158968, 0.000000] + turb_oneram6.test_vals = [-2.408665, -6.628342, 0.238581, 0.158951, 0.000000] turb_oneram6.timeout = 3200 test_list.append(turb_oneram6) @@ -275,7 +275,7 @@ def main(): turb_naca0012_sa.cfg_dir = "rans/naca0012" turb_naca0012_sa.cfg_file = "turb_NACA0012_sa.cfg" turb_naca0012_sa.test_iter = 5 - turb_naca0012_sa.test_vals = [-12.091778, -14.685322, 1.057665, 0.022971, 20.000000, -2.686203, 20.000000, -4.459916, 0] + turb_naca0012_sa.test_vals = [-12.049550, -16.193702, 1.058588, 0.022984, 20.000000, -3.437275, 20.000000, -4.705227, 0.000000] turb_naca0012_sa.timeout = 3200 test_list.append(turb_naca0012_sa) @@ -322,7 +322,7 @@ def main(): propeller.cfg_dir = "rans/propeller" propeller.cfg_file = "propeller.cfg" propeller.test_iter = 10 - propeller.test_vals = [-3.389724, -8.409502, 0.000048, 0.056344] + propeller.test_vals = [-3.389724, -8.410479, 0.000048, 0.056344] propeller.timeout = 3200 test_list.append(propeller) @@ -441,7 +441,7 @@ def main(): inc_turb_naca0012.cfg_dir = "incomp_rans/naca0012" inc_turb_naca0012.cfg_file = "naca0012.cfg" inc_turb_naca0012.test_iter = 20 - inc_turb_naca0012.test_vals = [-4.788495, -11.040578, 0.000023, 0.309503] + inc_turb_naca0012.test_vals = [-4.788495, -11.040895, 0.000023, 0.309502] test_list.append(inc_turb_naca0012) # NACA0012, SST_SUST @@ -465,7 +465,7 @@ def main(): inc_turb_wallfunction_flatplate_sa.cfg_dir = "wallfunctions/flatplate/incompressible_SA" inc_turb_wallfunction_flatplate_sa.cfg_file = "turb_SA_flatplate.cfg" inc_turb_wallfunction_flatplate_sa.test_iter = 10 - inc_turb_wallfunction_flatplate_sa.test_vals = [-6.894429, -5.717193, -6.743475, -4.242550, -9.587276, 10.000000, -2.879802, 0.001021, 0.003759] + inc_turb_wallfunction_flatplate_sa.test_vals = [-6.894206, -5.715970, -6.743740, -4.242551, -9.550272, 10.000000, -2.879370, 0.001021, 0.003759] test_list.append(inc_turb_wallfunction_flatplate_sa) #################### @@ -558,7 +558,7 @@ def main(): schubauer_klebanoff_transition.cfg_dir = "transition/Schubauer_Klebanoff" schubauer_klebanoff_transition.cfg_file = "transitional_BC_model_ConfigFile.cfg" schubauer_klebanoff_transition.test_iter = 10 - schubauer_klebanoff_transition.test_vals = [-8.284308, -13.240273, 0.000057, 0.007982] + schubauer_klebanoff_transition.test_vals = [-8.284308, -13.239421, 0.000057, 0.007982] test_list.append(schubauer_klebanoff_transition) ##################################### @@ -729,7 +729,7 @@ def main(): hb_rans_preconditioning.cfg_dir = "harmonic_balance/hb_rans_preconditioning" hb_rans_preconditioning.cfg_file = "davis.cfg" hb_rans_preconditioning.test_iter = 25 - hb_rans_preconditioning.test_vals = [-1.902097, 0.484069, 0.601483, 3.609005, -5.949355] + hb_rans_preconditioning.test_vals = [-1.902084, 0.484233, 0.601496, 3.609018, -5.943873] test_list.append(hb_rans_preconditioning) ###################################### @@ -796,7 +796,7 @@ def main(): ddes_flatplate.cfg_dir = "ddes/flatplate" ddes_flatplate.cfg_file = "ddes_flatplate.cfg" ddes_flatplate.test_iter = 10 - ddes_flatplate.test_vals = [-2.714713, -5.788301, -0.214960, 0.023758, 0.000000] + ddes_flatplate.test_vals = [-2.714713, -5.763293, -0.214960, 0.023758, 0.000000] ddes_flatplate.unsteady = True test_list.append(ddes_flatplate) @@ -856,7 +856,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 = [-7.701448, -8.512359, -6.014939, -6.468741, -5.801762, -4.607173, -5.551041, -5.300777, -3.804188, -5.256055, -5.765225, -3.609601, -2.229277, -2.883896, -0.563470] + Aachen_3D_restart.test_vals = [-7.701448, -8.512353, -6.014939, -6.468417, -5.801739, -4.607173, -5.550692, -5.300771, -3.804187, -5.256008, -5.765048, -3.609601, -2.229277, -2.883894, -0.563470] test_list.append(Aachen_3D_restart) # Jones APU Turbocharger restart @@ -889,8 +889,8 @@ def main(): 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] + multi_interface.test_vals = [-8.632229, -8.894737, -9.348730] + multi_interface.test_vals_aarch64 = [-8.632229, -8.894737, -9.348730] test_list.append(multi_interface) @@ -1129,6 +1129,15 @@ def main(): ### Method of Manufactured Solutions (MMS) ### ############################################## + # FVM, compressible, euler (periodic isentropic vortex) + mms_fvm_vortex = TestCase('mms_fvm_vortex') + mms_fvm_vortex.cfg_dir = "mms/fvm_euler" + mms_fvm_vortex.cfg_file = "inv_mms_vortex.cfg" + mms_fvm_vortex.test_iter = 10 + mms_fvm_vortex.test_vals = [-5.704300, -4.848072, 0.000000, 0.000000] + mms_fvm_vortex.unsteady = True + test_list.append(mms_fvm_vortex) + # FVM, compressible, laminar N-S mms_fvm_ns = TestCase('mms_fvm_ns') mms_fvm_ns.cfg_dir = "mms/fvm_navierstokes" diff --git a/TestCases/serial_regression_AD.py b/TestCases/serial_regression_AD.py index 89f64d54c22b..9cfb3177caa1 100644 --- a/TestCases/serial_regression_AD.py +++ b/TestCases/serial_regression_AD.py @@ -53,6 +53,16 @@ def main(): discadj_naca0012.test_vals = [-3.562611, -8.932639, -0.000000, 0.005608] test_list.append(discadj_naca0012) + # Inviscid NACA0012 (via discadj multizone driver) + discadj_naca0012_via_mz = TestCase('discadj_naca0012_via_mz') + discadj_naca0012_via_mz.cfg_dir = "cont_adj_euler/naca0012" + discadj_naca0012_via_mz.cfg_file = "inv_NACA0012_discadj_multizone.cfg" + discadj_naca0012_via_mz.test_iter = 100 + discadj_naca0012_via_mz.test_vals = [-3.563784, -5.975640, -6.326231, -8.929567] + discadj_naca0012_via_mz.enabled_with_tapetests = True + discadj_naca0012_via_mz.tapetest_vals = [0] + test_list.append(discadj_naca0012_via_mz) + # Inviscid Cylinder 3D (multiple markers) discadj_cylinder3D = TestCase('discadj_cylinder3D') discadj_cylinder3D.cfg_dir = "disc_adj_euler/cylinder3D" @@ -78,7 +88,7 @@ def main(): discadj_rans_naca0012_sa.cfg_dir = "disc_adj_rans/naca0012" discadj_rans_naca0012_sa.cfg_file = "turb_NACA0012_sa.cfg" discadj_rans_naca0012_sa.test_iter = 10 - discadj_rans_naca0012_sa.test_vals = [-2.996976, -0.196055, 0.000004, -0.000000, 5.000000, -3.971736, 5.000000, -10.464319] + discadj_rans_naca0012_sa.test_vals = [-2.996961, -0.199169, 0.000004, -0.000000, 5.000000, -3.971752, 5.000000, -10.457920] test_list.append(discadj_rans_naca0012_sa) # Adjoint turbulent NACA0012 SST @@ -183,7 +193,7 @@ def main(): discadj_heat.cfg_dir = "disc_adj_heat" discadj_heat.cfg_file = "disc_adj_heat.cfg" discadj_heat.test_iter = 10 - discadj_heat.test_vals = [-2.174678, 0.591525, 0.000000, 0.008748] + discadj_heat.test_vals = [-2.174674, 0.581857, 0.000000, 0.008748] test_list.append(discadj_heat) ################################### @@ -226,7 +236,7 @@ def main(): if test.tol == 0.0: test.tol = 0.00001 - pass_list = [ test.run_test(args.tsan, args.asan) for test in test_list ] + pass_list = [ test.run_test(args.tsan, args.asan, args.tapetests) for test in test_list ] ################################### ### Coupled RHT-CFD Adjoint ### @@ -243,7 +253,7 @@ def main(): discadj_rht.reference_file_aarch64 = "of_grad_cd_aarch64.csv.ref" discadj_rht.test_file = "of_grad_cd.csv" discadj_rht.enabled_with_asan = False - pass_list.append(discadj_rht.run_filediff(args.tsan, args.asan)) + pass_list.append(discadj_rht.run_filediff(args.tsan, args.asan, args.tapetests)) test_list.append(discadj_rht) ###################################### @@ -261,7 +271,7 @@ def main(): discadj_euler_py.reference_file_aarch64 = "of_grad_cd_disc_aarch64.dat.ref" discadj_euler_py.test_file = "of_grad_cd.dat" discadj_euler_py.enabled_with_asan = False - pass_list.append(discadj_euler_py.run_filediff(args.tsan, args.asan)) + pass_list.append(discadj_euler_py.run_filediff(args.tsan, args.asan, args.tapetests)) test_list.append(discadj_euler_py) # test discrete_adjoint with multiple ffd boxes @@ -275,7 +285,7 @@ def main(): discadj_multiple_ffd_py.reference_file_aarch64 = "of_grad_cd_aarch64.dat.ref" discadj_multiple_ffd_py.test_file = "of_grad_cd.dat" discadj_multiple_ffd_py.enabled_with_asan = False - pass_list.append(discadj_multiple_ffd_py.run_filediff(args.tsan, args.asan)) + pass_list.append(discadj_multiple_ffd_py.run_filediff(args.tsan, args.asan, args.tapetests)) test_list.append(discadj_multiple_ffd_py) # test direct_differentiation.py @@ -289,7 +299,7 @@ def main(): directdiff_euler_py.reference_file_aarch64 = "of_grad_directdiff_aarch64.dat.ref" directdiff_euler_py.test_file = "DIRECTDIFF/of_grad_directdiff.dat" directdiff_euler_py.enabled_with_asan = False - pass_list.append(directdiff_euler_py.run_filediff(args.tsan, args.asan)) + pass_list.append(directdiff_euler_py.run_filediff(args.tsan, args.asan, args.tapetests)) test_list.append(directdiff_euler_py) # test direct_differentiation.py with multiple ffd boxes @@ -303,7 +313,7 @@ def main(): directdiff_multiple_ffd_py.reference_file_aarch64 = "of_grad_directdiff_aarch64.dat.ref" directdiff_multiple_ffd_py.test_file = "DIRECTDIFF/of_grad_directdiff.dat" directdiff_multiple_ffd_py.enabled_with_asan = False - pass_list.append(directdiff_multiple_ffd_py.run_filediff(args.tsan, args.asan)) + pass_list.append(directdiff_multiple_ffd_py.run_filediff(args.tsan, args.asan, args.tapetests)) test_list.append(directdiff_multiple_ffd_py) # test continuous_adjoint.py, with multiple objectives @@ -330,7 +340,7 @@ def main(): pywrapper_FEA_AD_FlowLoad.new_output = False pywrapper_FEA_AD_FlowLoad.enabled_with_asan = False test_list.append(pywrapper_FEA_AD_FlowLoad) - pass_list.append(pywrapper_FEA_AD_FlowLoad.run_test(args.tsan, args.asan)) + pass_list.append(pywrapper_FEA_AD_FlowLoad.run_test(args.tsan, args.asan, args.tapetests)) # Flow AD Mesh Displacement Sensitivity pywrapper_CFD_AD_MeshDisp = TestCase('pywrapper_CFD_AD_MeshDisp') @@ -344,7 +354,7 @@ def main(): pywrapper_CFD_AD_MeshDisp.new_output = False pywrapper_CFD_AD_MeshDisp.enabled_with_asan = False test_list.append(pywrapper_CFD_AD_MeshDisp) - pass_list.append(pywrapper_CFD_AD_MeshDisp.run_test(args.tsan, args.asan)) + pass_list.append(pywrapper_CFD_AD_MeshDisp.run_test(args.tsan, args.asan, args.tapetests)) ################################### @@ -361,7 +371,7 @@ def main(): grad_smooth_naca0012.reference_file_aarch64 = "of_hess_aarch64.dat.ref" grad_smooth_naca0012.test_file = "of_hess.dat" grad_smooth_naca0012.enabled_with_asan = False - pass_list.append(grad_smooth_naca0012.run_filediff(args.tsan, args.asan)) + pass_list.append(grad_smooth_naca0012.run_filediff(args.tsan, args.asan, args.tapetests)) test_list.append(grad_smooth_naca0012) # Tests summary diff --git a/TestCases/tutorials.py b/TestCases/tutorials.py index 457d5697a0e0..4fcc604b40ac 100644 --- a/TestCases/tutorials.py +++ b/TestCases/tutorials.py @@ -291,7 +291,7 @@ def main(): tutorial_nicfd_nozzle_pinn.cfg_dir = "../Tutorials/compressible_flow/NICFD_nozzle/PhysicsInformed" tutorial_nicfd_nozzle_pinn.cfg_file = "config_NICFD_PINN.cfg" tutorial_nicfd_nozzle_pinn.test_iter = 20 - tutorial_nicfd_nozzle_pinn.test_vals = [-3.181747, -1.638856, -1.277037, 2.445964, -11.769570] + tutorial_nicfd_nozzle_pinn.test_vals = [-3.181747, -1.638856, -1.277037, 2.445964, -11.759632] tutorial_nicfd_nozzle_pinn.no_restart = True test_list.append(tutorial_nicfd_nozzle_pinn) @@ -301,7 +301,7 @@ def main(): tutorial_unst_naca0012.cfg_dir = "../Tutorials/compressible_flow/Unsteady_NACA0012" tutorial_unst_naca0012.cfg_file = "unsteady_naca0012.cfg" tutorial_unst_naca0012.test_iter = 520 - tutorial_unst_naca0012.test_vals = [520.000000, 0.000000, -5.290694, 0.000000, 0.317272, 0.820972, 0.002144, 0.012805] + tutorial_unst_naca0012.test_vals = [520.000000, 0.000000, -5.292632, 0.000000, 0.300303, 0.770888, 0.002399, 0.014070] tutorial_unst_naca0012.test_vals_aarch64 = [520.000000, 0.000000, -5.298777, 0.000000, 0.288956, 0.736706, 0.002419, 0.007134] tutorial_unst_naca0012.unsteady = True test_list.append(tutorial_unst_naca0012) @@ -311,7 +311,7 @@ def main(): propeller_var_load.cfg_dir = "../Tutorials/compressible_flow/ActuatorDisk_VariableLoad" propeller_var_load.cfg_file = "propeller_variable_load.cfg" propeller_var_load.test_iter = 20 - propeller_var_load.test_vals = [-1.830257, -4.535041, -0.000323, 0.171647] + propeller_var_load.test_vals = [-1.830257, -4.534990, -0.000323, 0.171646] propeller_var_load.timeout = 3200 test_list.append(propeller_var_load) diff --git a/TestCases/vandv.py b/TestCases/vandv.py index d563cc653de0..8c942e27b631 100644 --- a/TestCases/vandv.py +++ b/TestCases/vandv.py @@ -44,9 +44,9 @@ def main(): p30n30 = TestCase('30P30N') p30n30.cfg_dir = "vandv/rans/30p30n" p30n30.cfg_file = "config.cfg" - p30n30.test_iter = 20 - p30n30.test_vals = [-10.582183, -10.106601, -10.474910, -10.182549, -12.679336, 0.052181, 2.829820, 1.318613, -0.221374] - p30n30.test_vals_aarch64 = [-10.582183, -10.106601, -10.474910, -10.182549, -12.679336, 0.052181, 2.829820, 1.318613, -0.221374] + p30n30.test_iter = 5 + p30n30.test_vals = [-11.502310, -11.511459, -11.981995, -11.704990, -14.235571, 0.052235, 2.830394, 1.318894, -0.291726] + p30n30.test_vals_aarch64 = [-11.502310, -11.511459, -11.981995, -11.704990, -14.235571, 0.052235, 2.830394, 1.318894, -0.291726] test_list.append(p30n30) # flat plate - sst-v1994m @@ -72,8 +72,8 @@ def main(): swbli_sa.cfg_dir = "vandv/rans/swbli" swbli_sa.cfg_file = "config_sa.cfg" swbli_sa.test_iter = 5 - swbli_sa.test_vals = [-11.511182, -10.750503, -11.853919, -10.320019, -14.316261, 0.002238, -1.585259, 1.276300] - swbli_sa.test_vals_aarch64 = [-11.511278, -10.750583, -11.854073, -10.320108, -14.316261, 0.002238, -1.585354, 1.276300] + swbli_sa.test_vals = [-11.504424, -10.941741, -12.049925, -10.586263, -16.090385, 0.002242, -1.614365, 1.340100] + swbli_sa.test_vals_aarch64 = [-11.504424, -10.941741, -12.049925, -10.586263, -16.090385, 0.002242, -1.614365, 1.340100] test_list.append(swbli_sa) @@ -90,7 +90,7 @@ def main(): dsma661_sa.cfg_dir = "vandv/rans/dsma661" dsma661_sa.cfg_file = "dsma661_sa_config.cfg" dsma661_sa.test_iter = 5 - dsma661_sa.test_vals = [-11.013046, -8.140606, -8.989695, -5.978550, -10.593381, 0.155689, 0.024173] + dsma661_sa.test_vals = [-11.270155, -8.240208, -9.000574, -5.954878, -10.737828, 0.155687, 0.024232] test_list.append(dsma661_sa) # DSMA661 - SST-V2003m diff --git a/config_template.cfg b/config_template.cfg index 32e6369e6e33..07bf2474bc87 100644 --- a/config_template.cfg +++ b/config_template.cfg @@ -1852,6 +1852,10 @@ CFL_REDUCTION_TURB= 1.0 LOWER_LIMIT_K_FACTOR= 1e-15 LOWER_LIMIT_OMEGA_FACTOR= 1e-5 +% Use numerically computed exact Jacobians for turbulence models +% Slower per iteration but potentialy more stable and capable of higher CFL +USE_ACCURATE_TURB_JACOBIANS= NO +% % --------------------- HEAT NUMERICAL METHOD DEFINITION ----------------------% % % Value of the thermal diffusivity @@ -2041,13 +2045,20 @@ TOTAL_DV_PENALTY= 0.0 % % Parameters for the corresponding OF (allowed stress and KS multiplier). STRESS_PENALTY_PARAM= (1.0, 10.0) -% % ---------------- AUTOMATIC DIFFERENTIATION -------------------% % % Preaccumulation in the AD mode. PREACC= YES +% ---------------- AUTOMATIC DIFFERENTIATION (TAPE TEST DEBUG MODE) -------------------% +% +% Specify the kind of tape that is checked (OBJECTIVE_FUNCTION_TAPE, FULL_SOLVER_TAPE) +CHECK_TAPE_TYPE = FULL_SOLVER_TAPE +% +% Specify the variables for which the tape is checked (SOLVER_VARIABLES, SOLVER_VARIABLES_AND_MESH_COORDINATES) +CHECK_TAPE_VARIABLES = SOLVER_VARIABLES + % ---------------- MESH DEFORMATION PARAMETERS (NEW SOLVER) -------------------% % % Use the reformatted pseudo-elastic solver for grid deformation @@ -2376,6 +2387,9 @@ TABULAR_FORMAT= CSV % Set .precision(value) to specified value for SU2_DOT and HISTORY output. Useful for exact gradient validation. OUTPUT_PRECISION= 10 % +% For multizone problems, extend solution and restart filenames automatically by zone number +MULTIZONE_ADAPT_FILENAME= YES +% % Files to output % Possible formats : (TECPLOT_ASCII, TECPLOT, SURFACE_TECPLOT_ASCII, % SURFACE_TECPLOT, CSV, SURFACE_CSV, PARAVIEW_ASCII, PARAVIEW_LEGACY, SURFACE_PARAVIEW_ASCII, diff --git a/meson.build b/meson.build index 443cb1f4e543..2807fcc8f1b2 100644 --- a/meson.build +++ b/meson.build @@ -6,7 +6,7 @@ project('SU2', 'c', 'cpp', default_options: ['buildtype=release', 'warning_level=0', 'c_std=c99', - 'cpp_std=c++11']) + 'cpp_std=c++17']) fsmod = import('fs') if not fsmod.exists('su2preconfig.timestamp') @@ -101,6 +101,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 830b5f6574ee..58305769f163 100644 --- a/meson_options.txt +++ b/meson_options.txt @@ -24,7 +24,7 @@ option('enable-coolprop', type : 'boolean', value : false, description: 'enable option('enable-mlpcpp', type : 'boolean', value : false, description: 'enable MLPCpp support') option('enable-gprof', type : 'boolean', value : false, description: 'enable profiling through gprof') 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')