diff --git a/.github/workflows/regression.yml b/.github/workflows/regression.yml index c001fedf53f..af6f6c60875 100644 --- a/.github/workflows/regression.yml +++ b/.github/workflows/regression.yml @@ -52,12 +52,12 @@ jobs: key: ${{ matrix.config_set }}-${{ github.sha }} restore-keys: ${{ matrix.config_set }} - name: Pre Cleanup - uses: docker://ghcr.io/su2code/su2/build-su2:230813-0103 + uses: docker://ghcr.io/su2code/su2/build-su2:240320-1536 with: entrypoint: /bin/rm args: -rf install install_bin.tgz src ccache ${{ matrix.config_set }} - name: Build - uses: docker://ghcr.io/su2code/su2/build-su2:230813-0103 + uses: docker://ghcr.io/su2code/su2/build-su2:240320-1536 with: args: -b ${{github.ref}} -f "${{matrix.flags}}" - name: Compress binaries @@ -68,7 +68,7 @@ jobs: name: ${{ matrix.config_set }} path: install_bin.tgz - name: Post Cleanup - uses: docker://ghcr.io/su2code/su2/build-su2:230813-0103 + uses: docker://ghcr.io/su2code/su2/build-su2:240320-1536 with: entrypoint: /bin/rm args: -rf install install_bin.tgz src ccache ${{ matrix.config_set }} @@ -97,12 +97,12 @@ jobs: key: ${{ matrix.config_set }}-${{ github.sha }} restore-keys: ${{ matrix.config_set }} - name: Pre Cleanup - uses: docker://ghcr.io/su2code/su2/build-su2-tsan:230813-0103 + uses: docker://ghcr.io/su2code/su2/build-su2-tsan:240320-1536 with: entrypoint: /bin/rm args: -rf install install_bin.tgz src ccache ${{ matrix.config_set }} - name: Build - uses: docker://ghcr.io/su2code/su2/build-su2-tsan:230813-0103 + uses: docker://ghcr.io/su2code/su2/build-su2-tsan:240320-1536 with: args: -b ${{github.ref}} -f "${{matrix.flags}}" - name: Compress binaries @@ -113,7 +113,47 @@ jobs: name: ${{ matrix.config_set }} path: install_bin.tgz - name: Post Cleanup - uses: docker://ghcr.io/su2code/su2/build-su2-tsan:230813-0103 + uses: docker://ghcr.io/su2code/su2/build-su2-tsan:240320-1536 + with: + entrypoint: /bin/rm + args: -rf install install_bin.tgz src ccache ${{ matrix.config_set }} + build_asan: + name: Build SU2 (asan) + strategy: + fail-fast: false + matrix: + config_set: [BaseNoMPI-asan, ReverseNoMPI-asan] + include: + - config_set: BaseNoMPI-asan + flags: '--buildtype=debugoptimized -Denable-openblas=true -Dwith-mpi=disabled -Denable-mlpcpp=true --warnlevel=3 --werror' + - config_set: ReverseNoMPI-asan + flags: '--buildtype=debugoptimized -Denable-autodiff=true -Denable-normal=false -Dwith-mpi=disabled --warnlevel=3 --werror' + runs-on: ${{ inputs.runner || 'ubuntu-latest' }} + steps: + - name: Reduce ASLR entropy for asan + run: sudo sysctl -w vm.mmap_rnd_bits=28 + - name: Cache Object Files + uses: actions/cache@v4 + with: + path: ccache + key: ${{ matrix.config_set }}-${{ github.sha }} + restore-keys: ${{ matrix.config_set }} + - name: Pre Cleanup + uses: docker://ghcr.io/su2code/su2/build-su2-asan:240320-1536 + with: + entrypoint: /bin/rm + args: -rf install install_bin.tgz src ccache ${{ matrix.config_set }} + - name: Build + run: docker run --rm --cap-add SYS_PTRACE -v $(pwd):${{ github.workspace }} -w ${{ github.workspace }} ghcr.io/su2code/su2/build-su2-asan:240320-1536 -b ${{github.ref}} -f "${{matrix.flags}}" + - name: Compress binaries + run: tar -zcvf install_bin.tgz install/* + - name: Upload Binaries + uses: actions/upload-artifact@v4 + with: + name: ${{ matrix.config_set }} + path: install_bin.tgz + - name: Post Cleanup + uses: docker://ghcr.io/su2code/su2/build-su2-asan:240320-1536 with: entrypoint: /bin/rm args: -rf install install_bin.tgz src ccache ${{ matrix.config_set }} @@ -144,7 +184,7 @@ jobs: tag: OMP steps: - name: Pre Cleanup - uses: docker://ghcr.io/su2code/su2/test-su2:230813-0103 + 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 }} @@ -170,12 +210,12 @@ jobs: chmod a+x $BIN_FOLDER/* ls -lahR $BIN_FOLDER - name: Run Tests in Container - uses: docker://ghcr.io/su2code/su2/test-su2:230813-0103 + uses: docker://ghcr.io/su2code/su2/test-su2:240320-1536 with: # -t -c args: -b ${{github.ref}} -t develop -c develop -s ${{matrix.testscript}} - name: Cleanup - uses: docker://ghcr.io/su2code/su2/test-su2:230813-0103 + 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 }} @@ -192,7 +232,7 @@ jobs: - name: Reduce ASLR entropy for tsan run: sudo sysctl -w vm.mmap_rnd_bits=28 - name: Pre Cleanup - uses: docker://ghcr.io/su2code/su2/test-su2-tsan:230813-0103 + uses: docker://ghcr.io/su2code/su2/test-su2-tsan:240320-1536 with: entrypoint: /bin/rm args: -rf install install_bin.tgz src ccache ${{ matrix.config_set }} @@ -218,12 +258,59 @@ jobs: chmod a+x $BIN_FOLDER/* ls -lahR $BIN_FOLDER - name: Run Tests in Container - uses: docker://ghcr.io/su2code/su2/test-su2-tsan:230813-0103 + uses: docker://ghcr.io/su2code/su2/test-su2-tsan:240320-1536 with: # -t -c args: -b ${{github.ref}} -t develop -c develop -s ${{matrix.testscript}} -a "--tsan" - name: Cleanup - uses: docker://ghcr.io/su2code/su2/test-su2-tsan:230813-0103 + uses: docker://ghcr.io/su2code/su2/test-su2-tsan:240320-1536 + with: + entrypoint: /bin/rm + args: -rf install install_bin.tgz src ccache ${{ matrix.config_set }} + address_sanitizer_tests: + runs-on: ${{ inputs.runner || 'ubuntu-latest' }} + name: Address Sanitizer Tests + needs: build_asan + strategy: + fail-fast: false + matrix: + testscript: ['serial_regression.py', 'serial_regression_AD.py'] + steps: + - name: Reduce ASLR entropy for asan + run: sudo sysctl -w vm.mmap_rnd_bits=28 + - name: Pre Cleanup + uses: docker://ghcr.io/su2code/su2/test-su2-asan: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}NoMPI-asan" + 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-asan:240320-1536 + with: + # -t -c + args: -b ${{github.ref}} -t develop -c develop -s ${{matrix.testscript}} -a "--asan" + - name: Cleanup + uses: docker://ghcr.io/su2code/su2/test-su2-asan:240320-1536 with: entrypoint: /bin/rm args: -rf install install_bin.tgz src ccache ${{ matrix.config_set }} @@ -244,7 +331,7 @@ jobs: tag: MPI steps: - name: Pre Cleanup - uses: docker://ghcr.io/su2code/su2/test-su2:230813-0103 + 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 }} @@ -305,11 +392,11 @@ jobs: echo $PWD ls -lahR - name: Run Unit Tests - uses: docker://ghcr.io/su2code/su2/test-su2:230813-0103 + uses: docker://ghcr.io/su2code/su2/test-su2:240320-1536 with: entrypoint: install/bin/${{matrix.testdriver}} - name: Post Cleanup - uses: docker://ghcr.io/su2code/su2/test-su2:230813-0103 + 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 }} diff --git a/Common/include/CConfig.hpp b/Common/include/CConfig.hpp index dffd3e85fdd..bc32b5407b1 100644 --- a/Common/include/CConfig.hpp +++ b/Common/include/CConfig.hpp @@ -1224,7 +1224,13 @@ class CConfig { unsigned short nSpecies_Init; /*!< \brief Number of entries of SPECIES_INIT */ /*--- Additional flamelet solver options ---*/ - su2double flame_init[8]; /*!< \brief Initial solution parameters for flamelet solver.*/ + ///TODO: Add python wrapper initialization option + FLAMELET_INIT_TYPE flame_init_type = FLAMELET_INIT_TYPE::NONE; /*!< \brief Method for solution ignition for flamelet problems. */ + std::array flame_init; /*!< \brief Flame front initialization parameters. */ + std::array spark_init; /*!< \brief Spark ignition initialization parameters. */ + su2double* spark_reaction_rates; /*!< \brief Source terms for flamelet spark ignition option. */ + unsigned short nspark; /*!< \brief Number of source terms for spark initialization. */ + bool preferential_diffusion = false; /*!< \brief Preferential diffusion physics for flamelet solver.*/ /*--- lookup table ---*/ unsigned short n_scalars = 0; /*!< \brief Number of transported scalars for flamelet LUT approach. */ @@ -2138,13 +2144,44 @@ class CConfig { /*! * \brief Get the flame initialization. - * (x1,x2,x3) = flame offset. - * (x4,x5,x6) = flame normal, separating unburnt from burnt. + * (x1,x2,x3) = flame offset/spark center location. + * (x4,x5,x6) = flame normal, separating unburnt from burnt or + * spark radius, spark start iteration, spark duration. * (x7) = flame thickness, the length from unburnt to burnt conditions. * (x8) = flame burnt thickness, the length to stay at burnt conditions. - * \return Flame initialization for the flamelet model. + * \return Ignition initialization parameters for the flamelet model. + */ + const su2double* GetFlameInit() const { + switch (flame_init_type) + { + case FLAMELET_INIT_TYPE::FLAME_FRONT: + return flame_init.data(); + break; + case FLAMELET_INIT_TYPE::SPARK: + return spark_init.data(); + break; + default: + return nullptr; + break; + } + } + + /*! + * \brief Get species net reaction rates applied during spark ignition. + */ + const su2double* GetSpark() const { + return spark_reaction_rates; + } + + /*! + * \brief Preferential diffusion combustion problem. + */ + bool GetPreferentialDiffusion() const { return preferential_diffusion; } + + /*! + * \brief Define preferential diffusion combustion problem. */ - const su2double* GetFlameInit() const { return flame_init; } + inline void SetPreferentialDiffusion(bool input) { preferential_diffusion = input; } /*! * \brief Get the number of control variables for flamelet model. @@ -2190,6 +2227,11 @@ class CConfig { if (n_user_sources > 0) return user_source_names[i_user_source]; else return none; } + /*! + * \brief Get the ignition method used for combustion problems. + */ + FLAMELET_INIT_TYPE GetFlameletInitType() const { return flame_init_type; } + /*! * \brief Get the number of transported scalars for combustion. */ diff --git a/Common/include/containers/CLookUpTable.hpp b/Common/include/containers/CLookUpTable.hpp index 1c79535c5aa..230b009fc92 100644 --- a/Common/include/containers/CLookUpTable.hpp +++ b/Common/include/containers/CLookUpTable.hpp @@ -389,7 +389,6 @@ class CLookUpTable { /*! * \brief Find the table levels with constant z-values directly above and below query val_z. * \param[in] val_CV3 - Value of controlling variable 3. - * \param[in] within_limits - Whether query point lies within table bounds. * \returns Pair of inclusion level indices (first = lower level index, second = upper level index). */ std::pair FindInclusionLevels(const su2double val_CV3); @@ -410,6 +409,11 @@ class CLookUpTable { return limits_table_x[i_level]; } + /*! + * \brief Check whether requested set of variables are included in the table. + */ + bool CheckForVariables(const std::vector& vars_to_check) const; + /*! * \brief Returns the index to the variable in the lookup table. * \param[in] nameVar - Variable name for which to retrieve the column index. diff --git a/Common/include/geometry/CGeometry.hpp b/Common/include/geometry/CGeometry.hpp index 15a726cabec..688f71c7dce 100644 --- a/Common/include/geometry/CGeometry.hpp +++ b/Common/include/geometry/CGeometry.hpp @@ -426,7 +426,7 @@ class CGeometry { * \param[out] COUNT_PER_POINT - Number of communicated variables per point. * \param[out] MPI_TYPE - Enumerated type for the datatype of the quantity to be communicated. */ - void GetCommCountAndType(const CConfig* config, unsigned short commType, unsigned short& COUNT_PER_POINT, + void GetCommCountAndType(const CConfig* config, MPI_QUANTITIES commType, unsigned short& COUNT_PER_POINT, unsigned short& MPI_TYPE) const; /*! @@ -436,14 +436,14 @@ class CGeometry { * \param[in] config - Definition of the particular problem. * \param[in] commType - Enumerated type for the quantity to be communicated. */ - void InitiateComms(CGeometry* geometry, const CConfig* config, unsigned short commType) const; + void InitiateComms(CGeometry* geometry, const CConfig* config, MPI_QUANTITIES commType) const; /*! * \brief Routine to complete the set of non-blocking communications launched by InitiateComms() and unpacking of the * data into the geometry class. \param[in] geometry - Geometrical definition of the problem. \param[in] config - * Definition of the particular problem. \param[in] commType - Enumerated type for the quantity to be unpacked. */ - void CompleteComms(CGeometry* geometry, const CConfig* config, unsigned short commType); + void CompleteComms(CGeometry* geometry, const CConfig* config, MPI_QUANTITIES commType); /*! * \brief Get number of coordinates. diff --git a/Common/include/linear_algebra/CSysMatrix.hpp b/Common/include/linear_algebra/CSysMatrix.hpp index 573e6217d93..02747f8be07 100644 --- a/Common/include/linear_algebra/CSysMatrix.hpp +++ b/Common/include/linear_algebra/CSysMatrix.hpp @@ -92,7 +92,7 @@ struct CSysMatrixComms { */ template static void Initiate(const CSysVector& x, CGeometry* geometry, const CConfig* config, - unsigned short commType = SOLUTION_MATRIX); + MPI_QUANTITIES commType = MPI_QUANTITIES::SOLUTION_MATRIX); /*! * \brief Routine to complete the set of non-blocking communications launched by @@ -104,7 +104,7 @@ struct CSysMatrixComms { */ template static void Complete(CSysVector& x, CGeometry* geometry, const CConfig* config, - unsigned short commType = SOLUTION_MATRIX); + MPI_QUANTITIES commType = MPI_QUANTITIES::SOLUTION_MATRIX); }; /*! diff --git a/Common/include/option_structure.hpp b/Common/include/option_structure.hpp index 7199a310193..9bbe613c4bc 100644 --- a/Common/include/option_structure.hpp +++ b/Common/include/option_structure.hpp @@ -1333,11 +1333,38 @@ enum FLAMELET_SCALAR_SOURCES { * \brief Look-up operations for the flamelet scalar solver. */ enum FLAMELET_LOOKUP_OPS { - TD, /*!< \brief Thermochemical properties (temperature, density, diffusivity, etc.). */ + THERMO, /*!< \brief Thermochemical properties (temperature, density, diffusivity, etc.). */ + PREFDIF, /*!< \brief Preferential diffusion scalars. */ SOURCES, /*!< \brief Scalar source terms (controlling variables, passive species).*/ LOOKUP, /*!< \brief Passive look-up variables specified in config. */ }; +/*! + * \brief The preferential diffusion scalar indices for the preferential diffusion model. + */ +enum FLAMELET_PREF_DIFF_SCALARS { + I_BETA_PROGVAR, /*!< \brief Preferential diffusion scalar for the progress variable. */ + I_BETA_ENTH_THERMAL, /*!< \brief Preferential diffusion scalar for temperature. */ + I_BETA_ENTH, /*!< \brief Preferential diffusion scalar for total enthalpy. */ + I_BETA_MIXFRAC, /*!< \brief Preferential diffusion scalar for mixture fraction. */ + N_BETA_TERMS, /*!< \brief Total number of preferential diffusion scalars. */ +}; + +/*! + * \brief Flame initialization options for the flamelet solver. + */ +enum class FLAMELET_INIT_TYPE { + FLAME_FRONT, /*!< \brief Straight flame front. */ + SPARK, /*!< \brief Species reaction rate in a set location. */ + NONE, /*!< \brief No ignition, cold flow only. */ +}; + +static const MapType Flamelet_Init_Map = { + MakePair("NONE", FLAMELET_INIT_TYPE::NONE) + MakePair("FLAME_FRONT", FLAMELET_INIT_TYPE::FLAME_FRONT) + MakePair("SPARK", FLAMELET_INIT_TYPE::SPARK) +}; + /*! * \brief Types of subgrid scale models */ @@ -2471,7 +2498,7 @@ enum PERIODIC_QUANTITIES { /*! * \brief Vertex-based quantities exchanged in MPI point-to-point communications. */ -enum MPI_QUANTITIES { +enum class MPI_QUANTITIES { SOLUTION , /*!< \brief Conservative solution communication. */ SOLUTION_OLD , /*!< \brief Conservative solution old communication. */ SOLUTION_GRADIENT , /*!< \brief Conservative solution gradient communication. */ diff --git a/Common/src/CConfig.cpp b/Common/src/CConfig.cpp index 18da76556c4..63810ab3e03 100644 --- a/Common/src/CConfig.cpp +++ b/Common/src/CConfig.cpp @@ -987,6 +987,7 @@ void CConfig::SetPointersNull() { Species_Init = nullptr; Species_Clipping_Min = nullptr; Species_Clipping_Max = nullptr; + spark_reaction_rates = nullptr; /*--- Moving mesh pointers ---*/ @@ -1376,14 +1377,23 @@ void CConfig::SetConfig_Options() { /*!\brief SPECIES_CLIPPING_MIN \n DESCRIPTION: Minimum values for scalar clipping \ingroup Config*/ addDoubleListOption("SPECIES_CLIPPING_MIN", nSpecies_Clipping_Min, Species_Clipping_Min); - /*!\brief FLAME_INIT \n DESCRIPTION: flame initialization using the flamelet model \ingroup Config*/ + /*!\brief FLAME_INIT_METHOD \n DESCRIPTION: Ignition method for flamelet solver \n DEFAULT: no ignition; cold flow only. */ + addEnumOption("FLAME_INIT_METHOD", flame_init_type, Flamelet_Init_Map, FLAMELET_INIT_TYPE::NONE); + /*!\brief FLAME_INIT \n DESCRIPTION: flame front initialization using the flamelet model \ingroup Config*/ /*--- flame offset (x,y,z) ---*/ flame_init[0] = 0.0; flame_init[1] = 0.0; flame_init[2] = 0.0; /*--- flame normal (nx, ny, nz) ---*/ flame_init[3] = 1.0; flame_init[4] = 0.0; flame_init[5] = 0.0; /*--- flame thickness (x) and flame burnt thickness (after this thickness, we have unburnt conditions again) ---*/ flame_init[6] = 0.5e-3; flame_init[7] = 1.0; - addDoubleArrayOption("FLAME_INIT", 8,flame_init); + addDoubleArrayOption("FLAME_INIT", 8,flame_init.begin()); + + /*!\brief SPARK_INIT \n DESCRIPTION: spark initialization using the flamelet model \ingroup Config*/ + for (auto iSpark=0u; iSpark<6; ++iSpark) spark_init[iSpark]=0; + addDoubleArrayOption("SPARK_INIT", 6, spark_init.begin()); + + /*!\brief SPARK_REACTION_RATES \n DESCRIPTION: Net source term values applied to species within spark area during spark ignition. \ingroup Config*/ + addDoubleListOption("SPARK_REACTION_RATES", nspark, spark_reaction_rates); /*--- Options related to mass diffusivity and thereby the species solver. ---*/ @@ -2136,6 +2146,9 @@ void CConfig::SetConfig_Options() { /* DESCRIPTION: Names of the user scalar source terms. */ addStringListOption("USER_SOURCE_NAMES", n_user_sources, user_source_names); + /* DESCRIPTION: Enable preferential diffusion for FGM simulations. \n DEFAULT: false */ + addBoolOption("PREFERENTIAL_DIFFUSION", preferential_diffusion, false); + /*!\brief CONV_FILENAME \n DESCRIPTION: Output file convergence history (w/o extension) \n DEFAULT: history \ingroup Config*/ addStringOption("CONV_FILENAME", Conv_FileName, string("history")); /*!\brief BREAKDOWN_FILENAME \n DESCRIPTION: Output file forces breakdown \ingroup Config*/ diff --git a/Common/src/containers/CLookUpTable.cpp b/Common/src/containers/CLookUpTable.cpp index 6efa79bcef9..16c042efba7 100644 --- a/Common/src/containers/CLookUpTable.cpp +++ b/Common/src/containers/CLookUpTable.cpp @@ -802,3 +802,13 @@ bool CLookUpTable::IsInTriangle(su2double val_CV1, su2double val_CV2, unsigned l return (abs(area_tri - (area_0 + area_1 + area_2)) < area_tri * 1e-10); } + +bool CLookUpTable::CheckForVariables(const std::vector& vars_to_check) const { + for (const string& var_to_check : vars_to_check) { + if (!std::any_of(names_var.begin(), names_var.end(), + [var_to_check](const std::string& n) { return var_to_check == n; })) { + return false; + }; + } + return true; +} diff --git a/Common/src/geometry/CGeometry.cpp b/Common/src/geometry/CGeometry.cpp index e1f199999f9..5292fb69a13 100644 --- a/Common/src/geometry/CGeometry.cpp +++ b/Common/src/geometry/CGeometry.cpp @@ -543,29 +543,29 @@ void CGeometry::PostP2PSends(CGeometry* geometry, const CConfig* config, unsigne END_SU2_OMP_MASTER } -void CGeometry::GetCommCountAndType(const CConfig* config, unsigned short commType, unsigned short& COUNT_PER_POINT, +void CGeometry::GetCommCountAndType(const CConfig* config, MPI_QUANTITIES commType, unsigned short& COUNT_PER_POINT, unsigned short& MPI_TYPE) const { switch (commType) { - case COORDINATES: + case MPI_QUANTITIES::COORDINATES: COUNT_PER_POINT = nDim; MPI_TYPE = COMM_TYPE_DOUBLE; break; - case GRID_VELOCITY: + case MPI_QUANTITIES::GRID_VELOCITY: COUNT_PER_POINT = nDim; MPI_TYPE = COMM_TYPE_DOUBLE; break; - case COORDINATES_OLD: + case MPI_QUANTITIES::COORDINATES_OLD: if (config->GetTime_Marching() == TIME_MARCHING::DT_STEPPING_2ND) COUNT_PER_POINT = nDim * 2; else COUNT_PER_POINT = nDim; MPI_TYPE = COMM_TYPE_DOUBLE; break; - case MAX_LENGTH: + case MPI_QUANTITIES::MAX_LENGTH: COUNT_PER_POINT = 1; MPI_TYPE = COMM_TYPE_DOUBLE; break; - case NEIGHBORS: + case MPI_QUANTITIES::NEIGHBORS: COUNT_PER_POINT = 1; MPI_TYPE = COMM_TYPE_UNSIGNED_SHORT; break; @@ -575,7 +575,7 @@ void CGeometry::GetCommCountAndType(const CConfig* config, unsigned short commTy } } -void CGeometry::InitiateComms(CGeometry* geometry, const CConfig* config, unsigned short commType) const { +void CGeometry::InitiateComms(CGeometry* geometry, const CConfig* config, MPI_QUANTITIES commType) const { if (nP2PSend == 0) return; /*--- Local variables ---*/ @@ -633,15 +633,15 @@ void CGeometry::InitiateComms(CGeometry* geometry, const CConfig* config, unsign buf_offset = (msg_offset + iSend) * COUNT_PER_POINT; switch (commType) { - case COORDINATES: + case MPI_QUANTITIES::COORDINATES: vector = nodes->GetCoord(iPoint); for (iDim = 0; iDim < nDim; iDim++) bufDSend[buf_offset + iDim] = vector[iDim]; break; - case GRID_VELOCITY: + case MPI_QUANTITIES::GRID_VELOCITY: vector = nodes->GetGridVel(iPoint); for (iDim = 0; iDim < nDim; iDim++) bufDSend[buf_offset + iDim] = vector[iDim]; break; - case COORDINATES_OLD: + case MPI_QUANTITIES::COORDINATES_OLD: vector = nodes->GetCoord_n(iPoint); for (iDim = 0; iDim < nDim; iDim++) { bufDSend[buf_offset + iDim] = vector[iDim]; @@ -653,10 +653,10 @@ void CGeometry::InitiateComms(CGeometry* geometry, const CConfig* config, unsign } } break; - case MAX_LENGTH: + case MPI_QUANTITIES::MAX_LENGTH: bufDSend[buf_offset] = nodes->GetMaxLength(iPoint); break; - case NEIGHBORS: + case MPI_QUANTITIES::NEIGHBORS: bufSSend[buf_offset] = geometry->nodes->GetnNeighbor(iPoint); break; default: @@ -672,7 +672,7 @@ void CGeometry::InitiateComms(CGeometry* geometry, const CConfig* config, unsign } } -void CGeometry::CompleteComms(CGeometry* geometry, const CConfig* config, unsigned short commType) { +void CGeometry::CompleteComms(CGeometry* geometry, const CConfig* config, MPI_QUANTITIES commType) { if (nP2PRecv == 0) return; /*--- Local variables ---*/ @@ -734,21 +734,21 @@ void CGeometry::CompleteComms(CGeometry* geometry, const CConfig* config, unsign /*--- Store the data correctly depending on the quantity. ---*/ switch (commType) { - case COORDINATES: + case MPI_QUANTITIES::COORDINATES: for (iDim = 0; iDim < nDim; iDim++) nodes->SetCoord(iPoint, iDim, bufDRecv[buf_offset + iDim]); break; - case GRID_VELOCITY: + case MPI_QUANTITIES::GRID_VELOCITY: for (iDim = 0; iDim < nDim; iDim++) nodes->SetGridVel(iPoint, iDim, bufDRecv[buf_offset + iDim]); break; - case COORDINATES_OLD: + case MPI_QUANTITIES::COORDINATES_OLD: nodes->SetCoord_n(iPoint, &bufDRecv[buf_offset]); if (config->GetTime_Marching() == TIME_MARCHING::DT_STEPPING_2ND) nodes->SetCoord_n1(iPoint, &bufDRecv[buf_offset + nDim]); break; - case MAX_LENGTH: + case MPI_QUANTITIES::MAX_LENGTH: nodes->SetMaxLength(iPoint, bufDRecv[buf_offset]); break; - case NEIGHBORS: + case MPI_QUANTITIES::NEIGHBORS: nodes->SetnNeighbor(iPoint, bufSRecv[buf_offset]); break; default: @@ -2327,8 +2327,8 @@ void CGeometry::RegisterCoordinates() const { } void CGeometry::UpdateGeometry(CGeometry** geometry_container, CConfig* config) { - geometry_container[MESH_0]->InitiateComms(geometry_container[MESH_0], config, COORDINATES); - geometry_container[MESH_0]->CompleteComms(geometry_container[MESH_0], config, COORDINATES); + geometry_container[MESH_0]->InitiateComms(geometry_container[MESH_0], config, MPI_QUANTITIES::COORDINATES); + geometry_container[MESH_0]->CompleteComms(geometry_container[MESH_0], config, MPI_QUANTITIES::COORDINATES); geometry_container[MESH_0]->SetControlVolume(config, UPDATE); geometry_container[MESH_0]->SetBoundControlVolume(config, UPDATE); diff --git a/Common/src/geometry/CPhysicalGeometry.cpp b/Common/src/geometry/CPhysicalGeometry.cpp index 8d245ce1206..29f10a492b3 100644 --- a/Common/src/geometry/CPhysicalGeometry.cpp +++ b/Common/src/geometry/CPhysicalGeometry.cpp @@ -6364,8 +6364,8 @@ void CPhysicalGeometry::SetMaxLength(CConfig* config) { } END_SU2_OMP_FOR - InitiateComms(this, config, MAX_LENGTH); - CompleteComms(this, config, MAX_LENGTH); + InitiateComms(this, config, MPI_QUANTITIES::MAX_LENGTH); + CompleteComms(this, config, MPI_QUANTITIES::MAX_LENGTH); } void CPhysicalGeometry::MatchActuator_Disk(const CConfig* config) { diff --git a/Common/src/grid_movement/CFreeFormDefBox.cpp b/Common/src/grid_movement/CFreeFormDefBox.cpp index 78e4b384c9f..6cfd9cc8c85 100644 --- a/Common/src/grid_movement/CFreeFormDefBox.cpp +++ b/Common/src/grid_movement/CFreeFormDefBox.cpp @@ -144,7 +144,8 @@ CFreeFormDefBox::~CFreeFormDefBox() { for (iCornerPoints = 0; iCornerPoints < nCornerPoints; iCornerPoints++) delete[] Coord_Corner_Points[iCornerPoints]; delete[] Coord_Corner_Points; - for (iDim = 0; iDim < nDim; iDim++) { + /*--- nDim is 3 at allocation but might change later, so we used fixed 3 as upper bound for deallocation ---*/ + for (iDim = 0; iDim < 3; iDim++) { delete BlendingFunction[iDim]; } delete[] BlendingFunction; diff --git a/Common/src/grid_movement/CSurfaceMovement.cpp b/Common/src/grid_movement/CSurfaceMovement.cpp index 59b0cf7bfe8..5c29534fb39 100644 --- a/Common/src/grid_movement/CSurfaceMovement.cpp +++ b/Common/src/grid_movement/CSurfaceMovement.cpp @@ -33,12 +33,21 @@ CSurfaceMovement::CSurfaceMovement() : CGridMovement() { size = SU2_MPI::GetSize(); rank = SU2_MPI::GetRank(); + FFDBox = nullptr; nFFDBox = 0; nLevel = 0; FFDBoxDefinition = false; } -CSurfaceMovement::~CSurfaceMovement() = default; +CSurfaceMovement::~CSurfaceMovement() { + if (FFDBox != nullptr) { + for (unsigned int iFFDBox = 0; iFFDBox < MAX_NUMBER_FFD; ++iFFDBox) { + if (FFDBox[iFFDBox] != nullptr) delete FFDBox[iFFDBox]; + } + delete[] FFDBox; + FFDBox = nullptr; + } +} vector > CSurfaceMovement::SetSurface_Deformation(CGeometry* geometry, CConfig* config) { unsigned short iFFDBox, iDV, iLevel, iChild, iParent, jFFDBox, iMarker; @@ -60,7 +69,16 @@ vector > CSurfaceMovement::SetSurface_Deformation(CGeometry* g if (config->GetDesign_Variable(0) == FFD_SETTING) { /*--- Definition of the FFD deformation class ---*/ - FFDBox = new CFreeFormDefBox*[MAX_NUMBER_FFD]; + /*--- As this method might be called multiple times, properly delete old objects before allocating new ones. ---*/ + if (FFDBox != nullptr) { + for (iFFDBox = 0; iFFDBox < MAX_NUMBER_FFD; ++iFFDBox) { + if (FFDBox[iFFDBox] != nullptr) delete FFDBox[iFFDBox]; + } + delete[] FFDBox; + FFDBox = nullptr; + } + + FFDBox = new CFreeFormDefBox*[MAX_NUMBER_FFD](); /*--- Read the FFD information from the config file ---*/ @@ -167,7 +185,16 @@ vector > CSurfaceMovement::SetSurface_Deformation(CGeometry* g (config->GetDesign_Variable(0) == FFD_THICKNESS) || (config->GetDesign_Variable(0) == FFD_ANGLE_OF_ATTACK)) { /*--- Definition of the FFD deformation class ---*/ - FFDBox = new CFreeFormDefBox*[MAX_NUMBER_FFD]; + /*--- As this method might be called multiple times, properly delete old objects before allocating new ones. ---*/ + if (FFDBox != nullptr) { + for (iFFDBox = 0; iFFDBox < MAX_NUMBER_FFD; ++iFFDBox) { + if (FFDBox[iFFDBox] != nullptr) delete FFDBox[iFFDBox]; + } + delete[] FFDBox; + FFDBox = nullptr; + } + + FFDBox = new CFreeFormDefBox*[MAX_NUMBER_FFD](); /*--- Read the FFD information from the grid file ---*/ diff --git a/Common/src/grid_movement/CVolumetricMovement.cpp b/Common/src/grid_movement/CVolumetricMovement.cpp index 5b09a9cdd71..a3f91fa36c1 100644 --- a/Common/src/grid_movement/CVolumetricMovement.cpp +++ b/Common/src/grid_movement/CVolumetricMovement.cpp @@ -76,8 +76,8 @@ void CVolumetricMovement::UpdateGridCoord(CGeometry* geometry, CConfig* config) * Hence we still need a communication of the transformed coordinates, otherwise periodicity * is not maintained. ---*/ - geometry->InitiateComms(geometry, config, COORDINATES); - geometry->CompleteComms(geometry, config, COORDINATES); + geometry->InitiateComms(geometry, config, MPI_QUANTITIES::COORDINATES); + geometry->CompleteComms(geometry, config, MPI_QUANTITIES::COORDINATES); } void CVolumetricMovement::UpdateDualGrid(CGeometry* geometry, CConfig* config) { diff --git a/Common/src/linear_algebra/CSysMatrix.cpp b/Common/src/linear_algebra/CSysMatrix.cpp index 10747fadbcb..280d05b6b7f 100644 --- a/Common/src/linear_algebra/CSysMatrix.cpp +++ b/Common/src/linear_algebra/CSysMatrix.cpp @@ -226,7 +226,7 @@ void CSysMatrix::Initialize(unsigned long npoint, unsigned long npoi template void CSysMatrixComms::Initiate(const CSysVector& x, CGeometry* geometry, const CConfig* config, - unsigned short commType) { + MPI_QUANTITIES commType) { if (geometry->nP2PSend == 0) return; /*--- Local variables ---*/ @@ -236,13 +236,13 @@ void CSysMatrixComms::Initiate(const CSysVector& x, CGeometry* geometry, cons /*--- Create a boolean for reversing the order of comms. ---*/ - const bool reverse = (commType == SOLUTION_MATRIXTRANS); + const bool reverse = (commType == MPI_QUANTITIES::SOLUTION_MATRIXTRANS); /*--- Set the size of the data packet and type depending on quantity. ---*/ switch (commType) { - case SOLUTION_MATRIX: - case SOLUTION_MATRIXTRANS: + case MPI_QUANTITIES::SOLUTION_MATRIX: + case MPI_QUANTITIES::SOLUTION_MATRIXTRANS: break; default: SU2_MPI::Error("Unrecognized quantity for point-to-point MPI comms.", CURRENT_FUNCTION); @@ -265,7 +265,7 @@ void CSysMatrixComms::Initiate(const CSysVector& x, CGeometry* geometry, cons for (auto iMessage = 0; iMessage < geometry->nP2PSend; iMessage++) { switch (commType) { - case SOLUTION_MATRIX: { + case MPI_QUANTITIES::SOLUTION_MATRIX: { su2double* bufDSend = geometry->bufD_P2PSend; /*--- Get the offset for the start of this message. ---*/ @@ -294,7 +294,7 @@ void CSysMatrixComms::Initiate(const CSysVector& x, CGeometry* geometry, cons break; } - case SOLUTION_MATRIXTRANS: { + case MPI_QUANTITIES::SOLUTION_MATRIXTRANS: { /*--- We are going to communicate in reverse, so we use the recv buffer for the send instead. Also, all of the offsets and counts are derived from the recv data structures. ---*/ @@ -341,7 +341,7 @@ void CSysMatrixComms::Initiate(const CSysVector& x, CGeometry* geometry, cons } template -void CSysMatrixComms::Complete(CSysVector& x, CGeometry* geometry, const CConfig* config, unsigned short commType) { +void CSysMatrixComms::Complete(CSysVector& x, CGeometry* geometry, const CConfig* config, MPI_QUANTITIES commType) { if (geometry->nP2PRecv == 0) return; /*--- Local variables ---*/ @@ -366,7 +366,7 @@ void CSysMatrixComms::Complete(CSysVector& x, CGeometry* geometry, const CCon const auto source = status.MPI_SOURCE; switch (commType) { - case SOLUTION_MATRIX: { + case MPI_QUANTITIES::SOLUTION_MATRIX: { const su2double* bufDRecv = geometry->bufD_P2PRecv; /*--- We know the offsets based on the source rank. ---*/ @@ -400,7 +400,7 @@ void CSysMatrixComms::Complete(CSysVector& x, CGeometry* geometry, const CCon break; } - case SOLUTION_MATRIXTRANS: { + case MPI_QUANTITIES::SOLUTION_MATRIXTRANS: { /*--- We are going to communicate in reverse, so we use the send buffer for the recv instead. Also, all of the offsets and counts are derived from the send data structures. ---*/ @@ -1197,8 +1197,8 @@ void CSysMatrix::ComputePastixPreconditioner(const CSysVector(const CSysVector&, CGeometry*, const CConfig*, unsigned short); \ - template void CSysMatrixComms::Complete(CSysVector&, CGeometry*, const CConfig*, unsigned short); + template void CSysMatrixComms::Initiate(const CSysVector&, CGeometry*, const CConfig*, MPI_QUANTITIES); \ + template void CSysMatrixComms::Complete(CSysVector&, CGeometry*, const CConfig*, MPI_QUANTITIES); #define INSTANTIATE_MATRIX(TYPE) \ template class CSysMatrix; \ diff --git a/Docs/docmain.hpp b/Docs/docmain.hpp index 96743d0e371..5fefea50075 100644 --- a/Docs/docmain.hpp +++ b/Docs/docmain.hpp @@ -225,9 +225,3 @@ * \brief Classes for explicit (done by the programmer) vectorization (SIMD) of computations. * \ingroup Toolboxes */ - -/*! - * \defgroup Multi-Layer Perceptrons (MLP) - * \brief Data look up and interpolation via dense, feed-forward multi-layer perceptrons. - * \ingroup Toolboxes - */ diff --git a/SU2_CFD/include/fluid/CFluidFlamelet.hpp b/SU2_CFD/include/fluid/CFluidFlamelet.hpp index 12f370ad5eb..49273b2a9e6 100644 --- a/SU2_CFD/include/fluid/CFluidFlamelet.hpp +++ b/SU2_CFD/include/fluid/CFluidFlamelet.hpp @@ -42,9 +42,11 @@ class CFluidFlamelet final : public CFluidModel { enum LOOKUP_TD { TEMPERATURE, HEATCAPACITY, VISCOSITY, CONDUCTIVITY, DIFFUSIONCOEFFICIENT, MOLARWEIGHT, SIZE }; + su2double beta_progvar, beta_enth_thermal, beta_enth, beta_mixfrac; int rank; - bool include_mixture_fraction = false; /*!< \brief include mixture fraction in controlling variables. */ + bool include_mixture_fraction = false, /*!< \brief include mixture fraction in controlling variables. */ + preferential_diffusion = false; /*!< \brief use preferential diffusion physics. */ unsigned short n_scalars, n_lookups, n_user_scalars, /*!< \brief number of passive reactant species. */ n_control_vars; /*!< \brief number of controlling variables. */ @@ -62,12 +64,15 @@ class CFluidFlamelet final : public CFluidModel { vector LUT_idx_TD, LUT_idx_Sources, - LUT_idx_LookUp; + LUT_idx_LookUp, + LUT_idx_PD; /*--- Class variables for the multi-layer perceptron method ---*/ #ifdef USE_MLPCPP + size_t n_betas; MLPToolbox::CLookUp_ANN* lookup_mlp; /*!< \brief Multi-layer perceptron collection. */ MLPToolbox::CIOMap* iomap_TD; /*!< \brief Input-output map for thermochemical properties. */ + MLPToolbox::CIOMap* iomap_PD; /*!< \brief Input-output map for the preferential diffusion scalars. */ MLPToolbox::CIOMap* iomap_Sources; /*!< \brief Input-output map for species source terms. */ MLPToolbox::CIOMap* iomap_LookUp; /*!< \brief Input-output map for passive look-up terms. */ MLPToolbox::CIOMap* iomap_Current; @@ -76,10 +81,14 @@ class CFluidFlamelet final : public CFluidModel { vector scalars_vector; vector varnames_TD, /*!< \brief Lookup names for thermodynamic state variables. */ - varnames_Sources, varnames_LookUp; + varnames_Sources, /*!< \brief Lookup names for source terms. */ + varnames_LookUp, /*!< \brief Lookup names for passive look-up terms. */ + varnames_PD; /*!< \brief Lookup names for preferential diffusion scalars. */ vector val_vars_TD, /*!< \brief References to thermodynamic state variables. */ - val_vars_Sources, val_vars_LookUp; + val_vars_Sources, /*!< \brief References to source terms. */ + val_vars_LookUp, /*!< \brief References passive look-up terms. */ + val_vars_PD; /*!< \brief References to preferential diffusion scalars. */ void PreprocessLookUp(CConfig* config); @@ -107,14 +116,14 @@ class CFluidFlamelet final : public CFluidModel { void SetTDState_T(su2double val_temperature, const su2double* val_scalars = nullptr) override; /*! - * \brief Evaluate the flamelet manifold. - * \param[in] input_scalar - scalar solution. - * \param[in] input_varnames - names of variables to evaluate. - * \param[in] output_refs - output data. - * \param[out] Extrapolation - scalar solution is within bounds (0) or out of bounds (1). + * \brief Evaluate data-set for flamelet simulations. + * \param[in] input_scalar - controlling variables used to interpolate manifold. + * \param[in] lookup_type - look-up operation to be performed (FLAMELET_LOOKUP_OPS) + * \param[in] output_refs - output variables where interpolated results are stored. + * \param[out] Extrapolation - query data is within manifold bounds (0) or out of bounds (1). */ - inline unsigned long EvaluateDataSet(const vector& input_scalar, unsigned short lookup_type, - vector& output_refs) override; + unsigned long EvaluateDataSet(const vector& input_scalar, unsigned short lookup_type, + vector& output_refs); /*! * \brief Check for out-of-bounds condition for data set interpolation. @@ -142,4 +151,10 @@ class CFluidFlamelet final : public CFluidModel { * \param[out] Mu - value of the laminar viscosity */ inline su2double GetLaminarViscosity() override { return Mu; } + + /*! + * \brief Preferential diffusion as relevant phenomenon in flamelet simulations. + * \return Inclusion of preferential diffusion model. + */ + inline bool GetPreferentialDiffusion() const override { return preferential_diffusion; } }; diff --git a/SU2_CFD/include/fluid/CFluidModel.hpp b/SU2_CFD/include/fluid/CFluidModel.hpp index 7fc767b17d2..29f0becdb0d 100644 --- a/SU2_CFD/include/fluid/CFluidModel.hpp +++ b/SU2_CFD/include/fluid/CFluidModel.hpp @@ -46,34 +46,34 @@ class CLookUpTable; */ class CFluidModel { protected: - su2double StaticEnergy{0.0}; /*!< \brief Internal Energy. */ - su2double Entropy{0.0}; /*!< \brief Entropy. */ - su2double Density{0.0}; /*!< \brief Density. */ - su2double Pressure{0.0}; /*!< \brief Pressure. */ - su2double SoundSpeed2{0.0}; /*!< \brief SpeedSound. */ - su2double Temperature{0.0}; /*!< \brief Temperature. */ - su2double dPdrho_e{0.0}; /*!< \brief DpDd_e. */ - su2double dPde_rho{0.0}; /*!< \brief DpDe_d. */ - su2double dTdrho_e{0.0}; /*!< \brief DTDd_e. */ - su2double dTde_rho{0.0}; /*!< \brief DTDe_d. */ - su2double dhdrho_P{0.0}; /*!< \brief DhDrho_p. */ - su2double dhdP_rho{0.0}; /*!< \brief DhDp_rho. */ - su2double dsdrho_P{0.0}; /*!< \brief DsDrho_p. */ - su2double dsdP_rho{0.0}; /*!< \brief DsDp_rho. */ - su2double Cp{0.0}; /*!< \brief Specific Heat Capacity at constant pressure. */ - su2double Cv{0.0}; /*!< \brief Specific Heat Capacity at constant volume. */ - su2double Mu{0.0}; /*!< \brief Laminar viscosity. */ - su2double Mu_Turb{0.0}; /*!< \brief Eddy viscosity provided by a turbulence model. */ - su2double dmudrho_T{0.0}; /*!< \brief Partial derivative of viscosity w.r.t. density. */ - su2double dmudT_rho{0.0}; /*!< \brief Partial derivative of viscosity w.r.t. temperature. */ - su2double Kt{0.0}; /*!< \brief Thermal conductivity. */ - su2double dktdrho_T{0.0}; /*!< \brief Partial derivative of conductivity w.r.t. density. */ - su2double dktdT_rho{0.0}; /*!< \brief Partial derivative of conductivity w.r.t. temperature. */ - su2double mass_diffusivity{0.0}; /*!< \brief Mass Diffusivity */ + su2double StaticEnergy{0.0}; /*!< \brief Internal Energy. */ + su2double Entropy{0.0}; /*!< \brief Entropy. */ + su2double Density{0.0}; /*!< \brief Density. */ + su2double Pressure{0.0}; /*!< \brief Pressure. */ + su2double SoundSpeed2{0.0}; /*!< \brief SpeedSound. */ + su2double Temperature{0.0}; /*!< \brief Temperature. */ + su2double dPdrho_e{0.0}; /*!< \brief DpDd_e. */ + su2double dPde_rho{0.0}; /*!< \brief DpDe_d. */ + su2double dTdrho_e{0.0}; /*!< \brief DTDd_e. */ + su2double dTde_rho{0.0}; /*!< \brief DTDe_d. */ + su2double dhdrho_P{0.0}; /*!< \brief DhDrho_p. */ + su2double dhdP_rho{0.0}; /*!< \brief DhDp_rho. */ + su2double dsdrho_P{0.0}; /*!< \brief DsDrho_p. */ + su2double dsdP_rho{0.0}; /*!< \brief DsDp_rho. */ + su2double Cp{0.0}; /*!< \brief Specific Heat Capacity at constant pressure. */ + su2double Cv{0.0}; /*!< \brief Specific Heat Capacity at constant volume. */ + su2double Mu{0.0}; /*!< \brief Laminar viscosity. */ + su2double Mu_Turb{0.0}; /*!< \brief Eddy viscosity provided by a turbulence model. */ + su2double dmudrho_T{0.0}; /*!< \brief Partial derivative of viscosity w.r.t. density. */ + su2double dmudT_rho{0.0}; /*!< \brief Partial derivative of viscosity w.r.t. temperature. */ + su2double Kt{0.0}; /*!< \brief Thermal conductivity. */ + su2double dktdrho_T{0.0}; /*!< \brief Partial derivative of conductivity w.r.t. density. */ + su2double dktdT_rho{0.0}; /*!< \brief Partial derivative of conductivity w.r.t. temperature. */ + su2double mass_diffusivity{0.0}; /*!< \brief Mass Diffusivity */ unique_ptr LaminarViscosity; /*!< \brief Laminar Viscosity Model */ unique_ptr ThermalConductivity; /*!< \brief Thermal Conductivity Model */ - unique_ptr MassDiffusivity; /*!< \brief Mass Diffusivity Model */ + unique_ptr MassDiffusivity; /*!< \brief Mass Diffusivity Model */ /*! * \brief Instantiate the right type of viscosity model based on config. @@ -144,10 +144,16 @@ class CFluidModel { virtual inline unsigned short GetNScalars() const { return 0; } /*! - * \brief Evaluate data manifold for flamelet or data-driven fluid problems. - * \param[in] input - input data for manifold regression. + * \brief Evaluate data-set for flamelet or data-driven fluid simulations. + * \param[in] input_scalar - data manifold query data. + * \param[in] lookup_type - look-up operation to be performed. + * \param[in] output_refs - output variables where interpolated results are stored. + * \param[out] Extrapolation - query data is within manifold bounds (0) or out of bounds (1). */ - virtual unsigned long EvaluateDataSet(const vector &input_scalar, unsigned short lookup_type, vector &output_refs) { return 0; } + virtual unsigned long EvaluateDataSet(const vector& input_scalar, unsigned short lookup_type, + vector& output_refs) { + return 0; + } /*! * \brief Get fluid dynamic viscosity. @@ -331,7 +337,7 @@ class CFluidModel { * \brief Virtual member. * \param[in] T - Temperature value at the point. */ - virtual void SetTDState_T(su2double val_Temperature, const su2double* val_scalars = nullptr) { } + virtual void SetTDState_T(su2double val_Temperature, const su2double* val_scalars = nullptr) {} /*! * \brief Set fluid eddy viscosity provided by a turbulence model needed for computing effective thermal conductivity. @@ -344,6 +350,12 @@ class CFluidModel { */ virtual unsigned long GetExtrapolation() const { return 0; } + /*! + * \brief Get the state of the Preferential diffusion model for flamelet simulations. + * \return True if preferential diffusion model is active, false otherwise. + */ + virtual bool GetPreferentialDiffusion() const { return false; } + /*! * \brief Get number of Newton solver iterations. * \return Newton solver iteration count at termination. diff --git a/SU2_CFD/include/solvers/CFVMFlowSolverBase.hpp b/SU2_CFD/include/solvers/CFVMFlowSolverBase.hpp index 5e0e34b7a3e..26a69293004 100644 --- a/SU2_CFD/include/solvers/CFVMFlowSolverBase.hpp +++ b/SU2_CFD/include/solvers/CFVMFlowSolverBase.hpp @@ -710,8 +710,8 @@ class CFVMFlowSolverBase : public CSolver { /*--- MPI parallelization ---*/ - InitiateComms(geometry, config, MAX_EIGENVALUE); - CompleteComms(geometry, config, MAX_EIGENVALUE); + InitiateComms(geometry, config, MPI_QUANTITIES::MAX_EIGENVALUE); + CompleteComms(geometry, config, MPI_QUANTITIES::MAX_EIGENVALUE); } /*! @@ -780,8 +780,8 @@ class CFVMFlowSolverBase : public CSolver { /*--- MPI parallelization ---*/ - InitiateComms(geometry, config, SENSOR); - CompleteComms(geometry, config, SENSOR); + InitiateComms(geometry, config, MPI_QUANTITIES::SENSOR); + CompleteComms(geometry, config, MPI_QUANTITIES::SENSOR); } @@ -871,8 +871,8 @@ class CFVMFlowSolverBase : public CSolver { /*--- MPI solution ---*/ - InitiateComms(geometry, config, SOLUTION); - CompleteComms(geometry, config, SOLUTION); + InitiateComms(geometry, config, MPI_QUANTITIES::SOLUTION); + CompleteComms(geometry, config, MPI_QUANTITIES::SOLUTION); if (!adjoint) { /*--- For verification cases, compute the global error metrics. ---*/ @@ -996,8 +996,8 @@ class CFVMFlowSolverBase : public CSolver { CompletePeriodicComms(geometry, config, iPeriodic, PERIODIC_IMPLICIT); } - InitiateComms(geometry, config, SOLUTION); - CompleteComms(geometry, config, SOLUTION); + InitiateComms(geometry, config, MPI_QUANTITIES::SOLUTION); + CompleteComms(geometry, config, MPI_QUANTITIES::SOLUTION); /*--- For verification cases, compute the global error metrics. ---*/ ComputeVerificationError(geometry, config); diff --git a/SU2_CFD/include/solvers/CFVMFlowSolverBase.inl b/SU2_CFD/include/solvers/CFVMFlowSolverBase.inl index c3051311c30..b8efa1e695d 100644 --- a/SU2_CFD/include/solvers/CFVMFlowSolverBase.inl +++ b/SU2_CFD/include/solvers/CFVMFlowSolverBase.inl @@ -266,8 +266,8 @@ void CFVMFlowSolverBase::CommunicateInitialState(CGeometry* geometry, cons /*--- Perform the MPI communication of the solution ---*/ - InitiateComms(geometry, config, SOLUTION); - CompleteComms(geometry, config, SOLUTION); + InitiateComms(geometry, config, MPI_QUANTITIES::SOLUTION); + CompleteComms(geometry, config, MPI_QUANTITIES::SOLUTION); /*--- Store the initial CFL number for all grid points. ---*/ @@ -383,7 +383,7 @@ void CFVMFlowSolverBase::SetPrimitive_Gradient_GG(CGeometry* geometry, con bool reconstruction) { const auto& primitives = nodes->GetPrimitive(); auto& gradient = reconstruction ? nodes->GetGradient_Reconstruction() : nodes->GetGradient_Primitive(); - const auto comm = reconstruction? PRIMITIVE_GRAD_REC : PRIMITIVE_GRADIENT; + const auto comm = reconstruction? MPI_QUANTITIES::PRIMITIVE_GRAD_REC : MPI_QUANTITIES::PRIMITIVE_GRADIENT; const auto commPer = reconstruction? PERIODIC_PRIM_GG_R : PERIODIC_PRIM_GG; computeGradientsGreenGauss(this, comm, commPer, *geometry, *config, primitives, 0, nPrimVarGrad, gradient); @@ -408,7 +408,7 @@ void CFVMFlowSolverBase::SetPrimitive_Gradient_LS(CGeometry* geometry, con const auto& primitives = nodes->GetPrimitive(); auto& rmatrix = nodes->GetRmatrix(); auto& gradient = reconstruction ? nodes->GetGradient_Reconstruction() : nodes->GetGradient_Primitive(); - const auto comm = reconstruction? PRIMITIVE_GRAD_REC : PRIMITIVE_GRADIENT; + const auto comm = reconstruction? MPI_QUANTITIES::PRIMITIVE_GRAD_REC : MPI_QUANTITIES::PRIMITIVE_GRADIENT; computeGradientsLeastSquares(this, comm, commPer, *geometry, *config, weighted, primitives, 0, nPrimVarGrad, gradient, rmatrix); @@ -423,7 +423,7 @@ void CFVMFlowSolverBase::SetPrimitive_Limiter(CGeometry* geometry, const C auto& primMax = nodes->GetSolution_Max(); auto& limiter = nodes->GetLimiter_Primitive(); - computeLimiters(kindLimiter, this, PRIMITIVE_LIMITER, PERIODIC_LIM_PRIM_1, PERIODIC_LIM_PRIM_2, *geometry, *config, 0, + computeLimiters(kindLimiter, this, MPI_QUANTITIES::PRIMITIVE_LIMITER, PERIODIC_LIM_PRIM_1, PERIODIC_LIM_PRIM_2, *geometry, *config, 0, nPrimVarGrad, primitives, gradient, primMin, primMax, limiter); } @@ -929,8 +929,8 @@ void CFVMFlowSolverBase::LoadRestart_impl(CGeometry **geometry, CSolver ** /*--- Compute the grid velocities on the coarser levels. ---*/ if (iMesh) geometry[iMesh]->SetRestricted_GridVelocity(geometry[iMesh - 1]); else { - geometry[MESH_0]->InitiateComms(geometry[MESH_0], config, GRID_VELOCITY); - geometry[MESH_0]->CompleteComms(geometry[MESH_0], config, GRID_VELOCITY); + geometry[MESH_0]->InitiateComms(geometry[MESH_0], config, MPI_QUANTITIES::GRID_VELOCITY); + geometry[MESH_0]->CompleteComms(geometry[MESH_0], config, MPI_QUANTITIES::GRID_VELOCITY); } } } @@ -941,8 +941,8 @@ void CFVMFlowSolverBase::LoadRestart_impl(CGeometry **geometry, CSolver ** on the fine level in order to have all necessary quantities updated, especially if this is a turbulent simulation (eddy viscosity). ---*/ - solver[MESH_0][FLOW_SOL]->InitiateComms(geometry[MESH_0], config, SOLUTION); - solver[MESH_0][FLOW_SOL]->CompleteComms(geometry[MESH_0], config, SOLUTION); + solver[MESH_0][FLOW_SOL]->InitiateComms(geometry[MESH_0], config, MPI_QUANTITIES::SOLUTION); + solver[MESH_0][FLOW_SOL]->CompleteComms(geometry[MESH_0], config, MPI_QUANTITIES::SOLUTION); /*--- For turbulent/species simulations the flow preprocessing is done by the turbulence/species solver * after it loads its variables (they are needed to compute flow primitives). In case turbulence and species, the @@ -957,8 +957,8 @@ void CFVMFlowSolverBase::LoadRestart_impl(CGeometry **geometry, CSolver ** for (auto iMesh = 1u; iMesh <= config->GetnMGLevels(); iMesh++) { MultigridRestriction(*geometry[iMesh - 1], solver[iMesh - 1][FLOW_SOL]->GetNodes()->GetSolution(), *geometry[iMesh], solver[iMesh][FLOW_SOL]->GetNodes()->GetSolution()); - solver[iMesh][FLOW_SOL]->InitiateComms(geometry[iMesh], config, SOLUTION); - solver[iMesh][FLOW_SOL]->CompleteComms(geometry[iMesh], config, SOLUTION); + solver[iMesh][FLOW_SOL]->InitiateComms(geometry[iMesh], config, MPI_QUANTITIES::SOLUTION); + solver[iMesh][FLOW_SOL]->CompleteComms(geometry[iMesh], config, MPI_QUANTITIES::SOLUTION); if (config->GetKind_Turb_Model() == TURB_MODEL::NONE && config->GetKind_Species_Model() == SPECIES_MODEL::NONE) { diff --git a/SU2_CFD/include/solvers/CHeatSolver.hpp b/SU2_CFD/include/solvers/CHeatSolver.hpp index e82c96881d7..627983ec096 100644 --- a/SU2_CFD/include/solvers/CHeatSolver.hpp +++ b/SU2_CFD/include/solvers/CHeatSolver.hpp @@ -103,8 +103,8 @@ class CHeatSolver final : public CScalarSolver { * \param[in] config - Definition of the particular problem. * \note Calls a generic implementation after defining a SolverSpecificNumerics object. */ - inline void Viscous_Residual(unsigned long iEdge, CGeometry* geometry, CSolver** solver_container, - CNumerics* numerics, CConfig* config) override { + inline void Viscous_Residual(const unsigned long iEdge, const CGeometry* geometry, CSolver** solver_container, + CNumerics* numerics, const CConfig* config) override { const CVariable* flow_nodes = flow ? solver_container[FLOW_SOL]->GetNodes() : nullptr; const su2double const_diffusivity = config->GetThermalDiffusivity(); diff --git a/SU2_CFD/include/solvers/CScalarSolver.hpp b/SU2_CFD/include/solvers/CScalarSolver.hpp index ca84ea7c644..ca4ce3e2483 100644 --- a/SU2_CFD/include/solvers/CScalarSolver.hpp +++ b/SU2_CFD/include/solvers/CScalarSolver.hpp @@ -99,9 +99,9 @@ class CScalarSolver : public CSolver { * \param[in] config - Definition of the particular problem. */ template - FORCEINLINE void Viscous_Residual_impl(const SolverSpecificNumericsFunc& SolverSpecificNumerics, unsigned long iEdge, - CGeometry* geometry, CSolver** solver_container, CNumerics* numerics, - CConfig* config) { + FORCEINLINE void Viscous_Residual_impl(const SolverSpecificNumericsFunc& SolverSpecificNumerics, const unsigned long iEdge, + const CGeometry* geometry, CSolver** solver_container, CNumerics* numerics, + const CConfig* config) { const bool implicit = (config->GetKind_TimeIntScheme() == EULER_IMPLICIT); CFlowVariable* flowNodes = solver_container[FLOW_SOL] ? su2staticcast_p(solver_container[FLOW_SOL]->GetNodes()) : nullptr; @@ -307,8 +307,8 @@ class CScalarSolver : public CSolver { * \param[in] config - Definition of the particular problem. * \note Calls a generic implementation after defining a SolverSpecificNumerics object. */ - inline virtual void Viscous_Residual(unsigned long iEdge, CGeometry* geometry, CSolver** solver_container, - CNumerics* numerics, CConfig* config) { + inline virtual void Viscous_Residual(const unsigned long iEdge, const CGeometry* geometry, CSolver** solver_container, + CNumerics* numerics, const CConfig* config) { /*--- Define an empty object for solver specific numerics contribution. In case there are none, this default *--- implementation will be called ---*/ auto SolverSpecificNumerics = [&](unsigned long iPoint, unsigned long jPoint) {}; diff --git a/SU2_CFD/include/solvers/CScalarSolver.inl b/SU2_CFD/include/solvers/CScalarSolver.inl index 6f9aeafac4a..1777fa2e289 100644 --- a/SU2_CFD/include/solvers/CScalarSolver.inl +++ b/SU2_CFD/include/solvers/CScalarSolver.inl @@ -531,8 +531,8 @@ void CScalarSolver::CompleteImplicitIteration(CGeometry* geometry, CompletePeriodicComms(geometry, config, iPeriodic, PERIODIC_IMPLICIT); } - InitiateComms(geometry, config, SOLUTION); - CompleteComms(geometry, config, SOLUTION); + InitiateComms(geometry, config, MPI_QUANTITIES::SOLUTION); + CompleteComms(geometry, config, MPI_QUANTITIES::SOLUTION); } template diff --git a/SU2_CFD/include/solvers/CSolver.hpp b/SU2_CFD/include/solvers/CSolver.hpp index 2b114a7714b..5a52599cfe6 100644 --- a/SU2_CFD/include/solvers/CSolver.hpp +++ b/SU2_CFD/include/solvers/CSolver.hpp @@ -240,7 +240,7 @@ class CSolver { * \param[out] MPI_TYPE - Enumerated type for the datatype of the quantity to be communicated. */ void GetCommCountAndType(const CConfig* config, - unsigned short commType, + MPI_QUANTITIES commType, unsigned short &COUNT_PER_POINT, unsigned short &MPI_TYPE) const; @@ -252,7 +252,7 @@ class CSolver { */ void InitiateComms(CGeometry *geometry, const CConfig *config, - unsigned short commType); + MPI_QUANTITIES commType); /*! * \brief Routine to complete the set of non-blocking communications launched by InitiateComms() and unpacking of the data in the solver class. @@ -262,7 +262,7 @@ class CSolver { */ void CompleteComms(CGeometry *geometry, const CConfig *config, - unsigned short commType); + MPI_QUANTITIES commType); /*! * \brief Helper function to define the type and number of variables per point for each communication type. diff --git a/SU2_CFD/include/solvers/CSpeciesFlameletSolver.hpp b/SU2_CFD/include/solvers/CSpeciesFlameletSolver.hpp index e821e1d6a9e..1d1bd498230 100644 --- a/SU2_CFD/include/solvers/CSpeciesFlameletSolver.hpp +++ b/SU2_CFD/include/solvers/CSpeciesFlameletSolver.hpp @@ -37,7 +37,6 @@ */ class CSpeciesFlameletSolver final : public CSpeciesSolver { private: - vector conjugate_var; /*!< \brief CHT variables for each boundary and vertex. */ bool include_mixture_fraction = false; /*!< \brief include mixture fraction as a controlling variable. */ /*! * \brief Compute the preconditioner for low-Mach flows. @@ -63,8 +62,8 @@ class CSpeciesFlameletSolver final : public CSpeciesSolver { * \param[in] val_enth_out - pointer to output enthalpy variable. * \param[out] Converged - 0 if Newton solver converged, 1 if not. */ - unsigned long GetEnthFromTemp(CFluidModel* fluid_model, su2double const val_temp, - const su2double* scalar_solution, su2double* val_enth_out); + unsigned long GetEnthFromTemp(CFluidModel* fluid_model, su2double const val_temp, const su2double* scalar_solution, + su2double* val_enth_out); /*! * \brief Find maximum progress variable value within the manifold for the current solution. @@ -97,13 +96,23 @@ class CSpeciesFlameletSolver final : public CSpeciesSolver { unsigned long SetScalarLookUps(const CConfig* config, CFluidModel* fluid_model_local, unsigned long iPoint, const vector& scalars); + /*! + * \brief Retrieve the preferential diffusion scalar values from manifold. + * \param[in] config - definition of particular problem. + * \param[in] fluid_model_local - pointer to flamelet fluid model. + * \param[in] iPoint - node ID. + * \param[in] scalars - local scalar solution. + * \return - within manifold bounds (0) or outside manifold bounds (1). + */ + unsigned long SetPreferentialDiffusionScalars(const CConfig* config, CFluidModel* fluid_model_local, + unsigned long iPoint, const vector& scalars); + public: /*! - * \brief Constructor. + * \brief Define a Flamelet Generated Manifold species solver. * \param[in] geometry - Geometrical definition of the problem. * \param[in] config - Definition of the particular problem. * \param[in] iMesh - Index of the mesh in multigrid computations. - * \param[in] FluidModel */ CSpeciesFlameletSolver(CGeometry* geometry, CConfig* config, unsigned short iMesh); @@ -175,4 +184,16 @@ class CSpeciesFlameletSolver final : public CSpeciesSolver { */ void BC_ConjugateHeat_Interface(CGeometry* geometry, CSolver** solver_container, CNumerics* numerics, CConfig* config, unsigned short val_marker) override; + + /*! + * \brief Compute the fluxes due to viscous and preferential diffusion effects of the flamelet species at a particular edge. + * \param[in] iEdge - Edge for which we want to compute the flux + * \param[in] geometry - Geometrical definition of the problem. + * \param[in] solver_container - Container vector with all the solutions. + * \param[in] numerics - Description of the numerical method. + * \param[in] config - Definition of the particular problem. + * \note Calls a generic implementation after defining a SolverSpecificNumerics object. + */ + void Viscous_Residual(const unsigned long iEdge, const CGeometry* geometry, CSolver** solver_container, CNumerics* numerics, + const CConfig* config) final; }; diff --git a/SU2_CFD/include/solvers/CSpeciesSolver.hpp b/SU2_CFD/include/solvers/CSpeciesSolver.hpp index 1ee78d90286..4fbb9419117 100644 --- a/SU2_CFD/include/solvers/CSpeciesSolver.hpp +++ b/SU2_CFD/include/solvers/CSpeciesSolver.hpp @@ -89,8 +89,8 @@ class CSpeciesSolver : public CScalarSolver { * \param[in] config - Definition of the particular problem. * \note Calls a generic implementation after defining a SolverSpecificNumerics object. */ - void Viscous_Residual(unsigned long iEdge, CGeometry* geometry, CSolver** solver_container, CNumerics* numerics, - CConfig* config) final; + void Viscous_Residual(const unsigned long iEdge, const CGeometry* geometry, CSolver** solver_container, CNumerics* numerics, + const CConfig* config) override; /*! * \brief Impose the inlet boundary condition. diff --git a/SU2_CFD/include/solvers/CTransLMSolver.hpp b/SU2_CFD/include/solvers/CTransLMSolver.hpp index 220f2b01ed7..98595abb266 100644 --- a/SU2_CFD/include/solvers/CTransLMSolver.hpp +++ b/SU2_CFD/include/solvers/CTransLMSolver.hpp @@ -98,8 +98,8 @@ class CTransLMSolver final : public CTurbSolver { * \param[in] config - Definition of the particular problem. * \note Calls a generic implementation after defining a SolverSpecificNumerics object. */ - void Viscous_Residual(unsigned long iEdge, CGeometry* geometry, CSolver** solver_container, - CNumerics* numerics, CConfig* config) override; + void Viscous_Residual(const unsigned long iEdge, const CGeometry* geometry, CSolver** solver_container, + CNumerics* numerics, const CConfig* config) override; /*! * \brief Source term computation. diff --git a/SU2_CFD/include/solvers/CTurbSASolver.hpp b/SU2_CFD/include/solvers/CTurbSASolver.hpp index cb0307817e7..c79649799cc 100644 --- a/SU2_CFD/include/solvers/CTurbSASolver.hpp +++ b/SU2_CFD/include/solvers/CTurbSASolver.hpp @@ -125,8 +125,8 @@ class CTurbSASolver final : public CTurbSolver { * \param[in] config - Definition of the particular problem. * \note Calls a generic implementation after defining a SolverSpecificNumerics object. */ - void Viscous_Residual(unsigned long iEdge, CGeometry* geometry, CSolver** solver_container, - CNumerics* numerics, CConfig* config) override; + void Viscous_Residual(const unsigned long iEdge, const CGeometry* geometry, CSolver** solver_container, + CNumerics* numerics, const CConfig* config) override; /*! * \brief Source term computation. diff --git a/SU2_CFD/include/solvers/CTurbSSTSolver.hpp b/SU2_CFD/include/solvers/CTurbSSTSolver.hpp index f4cf89d38d0..bdb4227a191 100644 --- a/SU2_CFD/include/solvers/CTurbSSTSolver.hpp +++ b/SU2_CFD/include/solvers/CTurbSSTSolver.hpp @@ -107,8 +107,8 @@ class CTurbSSTSolver final : public CTurbSolver { * \param[in] config - Definition of the particular problem. * \note Calls a generic implementation after defining a SolverSpecificNumerics object. */ - void Viscous_Residual(unsigned long iEdge, CGeometry* geometry, CSolver** solver_container, - CNumerics* numerics, CConfig* config) override; + void Viscous_Residual(const unsigned long iEdge, const CGeometry* geometry, CSolver** solver_container, + CNumerics* numerics, const CConfig* config) override; /*! * \brief Source term computation. diff --git a/SU2_CFD/src/drivers/CDriver.cpp b/SU2_CFD/src/drivers/CDriver.cpp index 7160cc93e6f..14c478c0f47 100644 --- a/SU2_CFD/src/drivers/CDriver.cpp +++ b/SU2_CFD/src/drivers/CDriver.cpp @@ -912,8 +912,8 @@ void CDriver::InitializeGeometryFVM(CConfig *config, CGeometry **&geometry) { if ((rank == MASTER_NODE) && (size > SINGLE_NODE) && (iMGlevel == MESH_0)) cout << "Communicating number of neighbors." << endl; - geometry[iMGlevel]->InitiateComms(geometry[iMGlevel], config, NEIGHBORS); - geometry[iMGlevel]->CompleteComms(geometry[iMGlevel], config, NEIGHBORS); + geometry[iMGlevel]->InitiateComms(geometry[iMGlevel], config, MPI_QUANTITIES::NEIGHBORS); + geometry[iMGlevel]->CompleteComms(geometry[iMGlevel], config, MPI_QUANTITIES::NEIGHBORS); } } @@ -2389,7 +2389,9 @@ void CDriver::PreprocessDynamicMesh(CConfig *config, CGeometry **geometry, CSolv cout << "Setting dynamic mesh structure for zone "<< iZone + 1<<"." << endl; grid_movement = new CVolumetricMovement(geometry[MESH_0], config); - surface_movement = new CSurfaceMovement(); + if (surface_movement == nullptr) + surface_movement = new CSurfaceMovement(); + surface_movement->CopyBoundary(geometry[MESH_0], config); if (config->GetTime_Marching() == TIME_MARCHING::HARMONIC_BALANCE){ if (rank == MASTER_NODE) cout << endl << "Instance "<< iInst + 1 <<":" << endl; @@ -2511,7 +2513,7 @@ void CDriver::InitializeInterface(CConfig **config, CSolver***** solver, CGeomet const auto fluidZone = heat_target? donor : target; - if (config[fluidZone]->GetEnergy_Equation() || (config[fluidZone]->GetKind_Regime() == ENUM_REGIME::COMPRESSIBLE) + if (config[fluidZone]->GetEnergy_Equation() || (config[fluidZone]->GetKind_Regime() == ENUM_REGIME::COMPRESSIBLE) || (config[fluidZone]->GetKind_FluidModel() == ENUM_FLUIDMODEL::FLUID_FLAMELET)) interface_type = heat_target? CONJUGATE_HEAT_FS : CONJUGATE_HEAT_SF; else if (config[fluidZone]->GetWeakly_Coupled_Heat()) @@ -2691,7 +2693,7 @@ void CDriver::PreprocessTurbomachinery(CConfig** config, CGeometry**** geometry, } } - for (iZone = 0; iZone < nZone-1; iZone++) { + for (iZone = 0; iZone < nZone-1; iZone++) { geometry[nZone-1][INST_0][MESH_0]->SetAvgTurboGeoValues(config[iZone],geometry[iZone][INST_0][MESH_0], iZone); } diff --git a/SU2_CFD/src/drivers/CDriverBase.cpp b/SU2_CFD/src/drivers/CDriverBase.cpp index c964677e623..602e89e3ace 100644 --- a/SU2_CFD/src/drivers/CDriverBase.cpp +++ b/SU2_CFD/src/drivers/CDriverBase.cpp @@ -391,8 +391,8 @@ vector CDriverBase::GetMarkerVertexNormals(unsigned short iMarker } void CDriverBase::CommunicateMeshDisplacements() { - solver_container[selected_zone][INST_0][MESH_0][MESH_SOL]->InitiateComms(main_geometry, main_config, MESH_DISPLACEMENTS); - solver_container[selected_zone][INST_0][MESH_0][MESH_SOL]->CompleteComms(main_geometry, main_config, MESH_DISPLACEMENTS); + solver_container[selected_zone][INST_0][MESH_0][MESH_SOL]->InitiateComms(main_geometry, main_config, MPI_QUANTITIES::MESH_DISPLACEMENTS); + solver_container[selected_zone][INST_0][MESH_0][MESH_SOL]->CompleteComms(main_geometry, main_config, MPI_QUANTITIES::MESH_DISPLACEMENTS); } map CDriverBase::GetSolverIndices() const { diff --git a/SU2_CFD/src/fluid/CFluidFlamelet.cpp b/SU2_CFD/src/fluid/CFluidFlamelet.cpp index 3e7c67c2c72..18dbcd2a43a 100644 --- a/SU2_CFD/src/fluid/CFluidFlamelet.cpp +++ b/SU2_CFD/src/fluid/CFluidFlamelet.cpp @@ -25,6 +25,8 @@ * License along with SU2. If not, see . */ +#include +#include #include "../include/fluid/CFluidFlamelet.hpp" #include "../../../Common/include/containers/CLookUpTable.hpp" #if defined(HAVE_MLPCPP) @@ -47,7 +49,6 @@ CFluidFlamelet::CFluidFlamelet(CConfig* config, su2double value_pressure_operati cout << "Number of user scalars: " << n_user_scalars << endl; cout << "Number of control variables: " << n_control_vars << endl; } - scalars_vector.resize(n_scalars); table_scalar_names.resize(n_scalars); @@ -93,31 +94,31 @@ CFluidFlamelet::CFluidFlamelet(CConfig* config, su2double value_pressure_operati Pressure = value_pressure_operating; PreprocessLookUp(config); + + if (rank == MASTER_NODE) { + cout << "Preferential diffusion: " << (preferential_diffusion ? "Enabled" : "Disabled") << endl; + } } CFluidFlamelet::~CFluidFlamelet() { - switch (Kind_DataDriven_Method) { - case ENUM_DATADRIVEN_METHOD::LUT: - delete look_up_table; - break; - case ENUM_DATADRIVEN_METHOD::MLP: + if (Kind_DataDriven_Method == ENUM_DATADRIVEN_METHOD::LUT) + delete look_up_table; #ifdef USE_MLPCPP - delete iomap_TD; - delete iomap_Sources; - delete iomap_LookUp; - delete lookup_mlp; + if (Kind_DataDriven_Method == ENUM_DATADRIVEN_METHOD::MLP) { + delete iomap_TD; + delete iomap_Sources; + delete iomap_LookUp; + delete lookup_mlp; + if (preferential_diffusion) delete iomap_PD; + } #endif - break; - default: - break; - } } void CFluidFlamelet::SetTDState_T(su2double val_temperature, const su2double* val_scalars) { for (auto iVar = 0u; iVar < n_scalars; iVar++) scalars_vector[iVar] = val_scalars[iVar]; /*--- Add all quantities and their names to the look up vectors. ---*/ - EvaluateDataSet(scalars_vector, FLAMELET_LOOKUP_OPS::TD, val_vars_TD); + EvaluateDataSet(scalars_vector, FLAMELET_LOOKUP_OPS::THERMO, val_vars_TD); Temperature = val_vars_TD[LOOKUP_TD::TEMPERATURE]; Cp = val_vars_TD[LOOKUP_TD::HEATCAPACITY]; @@ -183,9 +184,51 @@ void CFluidFlamelet::PreprocessLookUp(CConfig* config) { /*--- Passive look-up terms ---*/ size_t n_lookups = config->GetNLookups(); - varnames_LookUp.resize(n_lookups); - val_vars_LookUp.resize(n_lookups); - for (auto iLookup = 0u; iLookup < n_lookups; iLookup++) varnames_LookUp[iLookup] = config->GetLookupName(iLookup); + if (n_lookups == 0) { + varnames_LookUp.resize(1); + val_vars_LookUp.resize(1); + varnames_LookUp[0] = "NULL"; + } else { + varnames_LookUp.resize(n_lookups); + val_vars_LookUp.resize(n_lookups); + for (auto iLookup = 0u; iLookup < n_lookups; iLookup++) varnames_LookUp[iLookup] = config->GetLookupName(iLookup); + } + + /*--- Preferential diffusion scalars ---*/ + varnames_PD.resize(FLAMELET_PREF_DIFF_SCALARS::N_BETA_TERMS); + val_vars_PD.resize(FLAMELET_PREF_DIFF_SCALARS::N_BETA_TERMS); + + varnames_PD[FLAMELET_PREF_DIFF_SCALARS::I_BETA_PROGVAR] = "Beta_ProgVar"; + varnames_PD[FLAMELET_PREF_DIFF_SCALARS::I_BETA_ENTH_THERMAL] = "Beta_Enth_Thermal"; + varnames_PD[FLAMELET_PREF_DIFF_SCALARS::I_BETA_ENTH] = "Beta_Enth"; + varnames_PD[FLAMELET_PREF_DIFF_SCALARS::I_BETA_MIXFRAC] = "Beta_MixFrac"; + + val_vars_PD[FLAMELET_PREF_DIFF_SCALARS::I_BETA_PROGVAR] = beta_progvar; + val_vars_PD[FLAMELET_PREF_DIFF_SCALARS::I_BETA_ENTH_THERMAL] = beta_enth_thermal; + val_vars_PD[FLAMELET_PREF_DIFF_SCALARS::I_BETA_ENTH] = beta_enth; + val_vars_PD[FLAMELET_PREF_DIFF_SCALARS::I_BETA_MIXFRAC] = beta_mixfrac; + + preferential_diffusion = config->GetPreferentialDiffusion(); + switch (Kind_DataDriven_Method) { + case ENUM_DATADRIVEN_METHOD::LUT: + preferential_diffusion = look_up_table->CheckForVariables(varnames_PD); + break; + case ENUM_DATADRIVEN_METHOD::MLP: +#ifdef USE_MLPCPP + n_betas = 0; + for (auto iMLP = 0u; iMLP < config->GetNDataDriven_Files(); iMLP++) { + auto outputMap = lookup_mlp->FindVariableIndices(iMLP, varnames_PD, false); + n_betas += outputMap.size(); + } + preferential_diffusion = (n_betas == varnames_PD.size()); +#endif + break; + default: + break; + } + + if (!preferential_diffusion && config->GetPreferentialDiffusion()) + SU2_MPI::Error("Preferential diffusion scalars not included in flamelet manifold.", CURRENT_FUNCTION); if (Kind_DataDriven_Method == ENUM_DATADRIVEN_METHOD::MLP) { #ifdef USE_MLPCPP @@ -195,6 +238,10 @@ void CFluidFlamelet::PreprocessLookUp(CConfig* config) { lookup_mlp->PairVariableswithMLPs(*iomap_TD); lookup_mlp->PairVariableswithMLPs(*iomap_Sources); lookup_mlp->PairVariableswithMLPs(*iomap_LookUp); + if (preferential_diffusion) { + iomap_PD = new MLPToolbox::CIOMap(controlling_variable_names, varnames_PD); + lookup_mlp->PairVariableswithMLPs(*iomap_PD); + } #endif } else { for (auto iVar=0u; iVar < varnames_TD.size(); iVar++) { @@ -210,37 +257,52 @@ void CFluidFlamelet::PreprocessLookUp(CConfig* config) { LUT_idx_Sources.push_back(LUT_idx); } for (auto iVar=0u; iVar < varnames_LookUp.size(); iVar++) { - LUT_idx_LookUp.push_back(look_up_table->GetIndexOfVar(varnames_LookUp[iVar])); + unsigned long LUT_idx; + if (noSource(varnames_LookUp[iVar])) + LUT_idx = look_up_table->GetNullIndex(); + else + LUT_idx = look_up_table->GetIndexOfVar(varnames_LookUp[iVar]); + LUT_idx_LookUp.push_back(LUT_idx); + } + if (preferential_diffusion) { + for (auto iVar=0u; iVar < varnames_PD.size(); iVar++) { + LUT_idx_PD.push_back(look_up_table->GetIndexOfVar(varnames_PD[iVar])); + } } } } unsigned long CFluidFlamelet::EvaluateDataSet(const vector& input_scalar, unsigned short lookup_type, vector& output_refs) { + AD::StartPreacc(); + for (auto iVar = 0u; iVar < input_scalar.size(); iVar++) AD::SetPreaccIn(input_scalar[iVar]); + su2double val_enth = input_scalar[I_ENTH]; su2double val_prog = input_scalar[I_PROGVAR]; su2double val_mixfrac = include_mixture_fraction ? input_scalar[I_MIXFRAC] : 0.0; - vector varnames; vector val_vars; vector refs_vars; vector LUT_idx; switch (lookup_type) { - case FLAMELET_LOOKUP_OPS::TD: - varnames = varnames_TD; + case FLAMELET_LOOKUP_OPS::THERMO: LUT_idx = LUT_idx_TD; #ifdef USE_MLPCPP iomap_Current = iomap_TD; +#endif + break; + case FLAMELET_LOOKUP_OPS::PREFDIF: + LUT_idx = LUT_idx_PD; +#ifdef USE_MLPCPP + iomap_Current = iomap_PD; #endif break; case FLAMELET_LOOKUP_OPS::SOURCES: - varnames = varnames_Sources; LUT_idx = LUT_idx_Sources; #ifdef USE_MLPCPP iomap_Current = iomap_Sources; #endif break; case FLAMELET_LOOKUP_OPS::LOOKUP: - varnames = varnames_LookUp; LUT_idx = LUT_idx_LookUp; #ifdef USE_MLPCPP iomap_Current = iomap_LookUp; @@ -249,13 +311,14 @@ unsigned long CFluidFlamelet::EvaluateDataSet(const vector& input_sca default: break; } - if (output_refs.size() != varnames.size()) - SU2_MPI::Error(string("Output vector size incompatible with manifold lookup operation."), CURRENT_FUNCTION); + /*--- Add all quantities and their names to the look up vectors. ---*/ bool inside; switch (Kind_DataDriven_Method) { case ENUM_DATADRIVEN_METHOD::LUT: + if (output_refs.size() != LUT_idx.size()) + SU2_MPI::Error(string("Output vector size incompatible with manifold lookup operation."), CURRENT_FUNCTION); if (include_mixture_fraction) { inside = look_up_table->LookUp_XYZ(LUT_idx, output_refs, val_prog, val_enth, val_mixfrac); } else { @@ -274,6 +337,7 @@ unsigned long CFluidFlamelet::EvaluateDataSet(const vector& input_sca default: break; } - + for (auto iVar = 0u; iVar < output_refs.size(); iVar++) AD::SetPreaccOut(output_refs[iVar]); + AD::EndPreacc(); return extrapolation; } diff --git a/SU2_CFD/src/integration/CIntegration.cpp b/SU2_CFD/src/integration/CIntegration.cpp index 87e34cd1c8a..f53d50639cc 100644 --- a/SU2_CFD/src/integration/CIntegration.cpp +++ b/SU2_CFD/src/integration/CIntegration.cpp @@ -171,14 +171,13 @@ void CIntegration::Space_Integration(CGeometry *geometry, break; case CHT_WALL_INTERFACE: - if ( MainSolver == FLOW_SOL && (config->GetKind_FluidModel() == FLUID_FLAMELET) ){ - solver_container[MainSolver]->BC_Isothermal_Wall(geometry, solver_container, conv_bound_numerics, visc_bound_numerics, config, iMarker); - break; - } - if ((MainSolver == SPECIES_SOL) || (MainSolver == HEAT_SOL) || ((MainSolver == FLOW_SOL) && ((config->GetKind_Regime() == ENUM_REGIME::COMPRESSIBLE) || config->GetEnergy_Equation()))) { - solver_container[MainSolver]->BC_ConjugateHeat_Interface(geometry, solver_container, conv_bound_numerics, config, iMarker); - } - else { + if ((MainSolver == FLOW_SOL && (config->GetKind_FluidModel() == FLUID_FLAMELET)) || + (MainSolver == SPECIES_SOL) || (MainSolver == HEAT_SOL) || + ((MainSolver == FLOW_SOL) && + ((config->GetKind_Regime() == ENUM_REGIME::COMPRESSIBLE) || config->GetEnergy_Equation()))) { + solver_container[MainSolver]->BC_ConjugateHeat_Interface(geometry, solver_container, conv_bound_numerics, + config, iMarker); + } else { solver_container[MainSolver]->BC_HeatFlux_Wall(geometry, solver_container, conv_bound_numerics, visc_bound_numerics, config, iMarker); } break; diff --git a/SU2_CFD/src/integration/CMultiGridIntegration.cpp b/SU2_CFD/src/integration/CMultiGridIntegration.cpp index a6750435d97..f8a7c1a9645 100644 --- a/SU2_CFD/src/integration/CMultiGridIntegration.cpp +++ b/SU2_CFD/src/integration/CMultiGridIntegration.cpp @@ -366,8 +366,8 @@ void CMultiGridIntegration::GetProlongated_Correction(unsigned short RunTime_EqS /*--- MPI the set solution old ---*/ - sol_coarse->InitiateComms(geo_coarse, config, SOLUTION_OLD); - sol_coarse->CompleteComms(geo_coarse, config, SOLUTION_OLD); + sol_coarse->InitiateComms(geo_coarse, config, MPI_QUANTITIES::SOLUTION_OLD); + sol_coarse->CompleteComms(geo_coarse, config, MPI_QUANTITIES::SOLUTION_OLD); SU2_OMP_FOR_STAT(roundUpDiv(geo_coarse->GetnPointDomain(), omp_get_num_threads())) for (Point_Coarse = 0; Point_Coarse < geo_coarse->GetnPointDomain(); Point_Coarse++) { @@ -479,8 +479,8 @@ void CMultiGridIntegration::SetProlongated_Correction(CSolver *sol_fine, CGeomet /*--- MPI the new interpolated solution ---*/ - sol_fine->InitiateComms(geo_fine, config, SOLUTION); - sol_fine->CompleteComms(geo_fine, config, SOLUTION); + sol_fine->InitiateComms(geo_fine, config, MPI_QUANTITIES::SOLUTION); + sol_fine->CompleteComms(geo_fine, config, MPI_QUANTITIES::SOLUTION); } @@ -608,8 +608,8 @@ void CMultiGridIntegration::SetRestricted_Solution(unsigned short RunTime_EqSyst /*--- MPI the new interpolated solution ---*/ - sol_coarse->InitiateComms(geo_coarse, config, SOLUTION); - sol_coarse->CompleteComms(geo_coarse, config, SOLUTION); + sol_coarse->InitiateComms(geo_coarse, config, MPI_QUANTITIES::SOLUTION); + sol_coarse->CompleteComms(geo_coarse, config, MPI_QUANTITIES::SOLUTION); } diff --git a/SU2_CFD/src/integration/CSingleGridIntegration.cpp b/SU2_CFD/src/integration/CSingleGridIntegration.cpp index 56b74282cde..de70accb876 100644 --- a/SU2_CFD/src/integration/CSingleGridIntegration.cpp +++ b/SU2_CFD/src/integration/CSingleGridIntegration.cpp @@ -111,8 +111,8 @@ void CSingleGridIntegration::SetRestricted_Solution(unsigned short RunTime_EqSys CGeometry *geo_fine, CGeometry *geo_coarse, CConfig *config) { CSolver::MultigridRestriction(*geo_fine, sol_fine->GetNodes()->GetSolution(), *geo_coarse, sol_coarse->GetNodes()->GetSolution()); - sol_coarse->InitiateComms(geo_coarse, config, SOLUTION); - sol_coarse->CompleteComms(geo_coarse, config, SOLUTION); + sol_coarse->InitiateComms(geo_coarse, config, MPI_QUANTITIES::SOLUTION); + sol_coarse->CompleteComms(geo_coarse, config, MPI_QUANTITIES::SOLUTION); } void CSingleGridIntegration::SetRestricted_EddyVisc(unsigned short RunTime_EqSystem, CSolver *sol_fine, CSolver *sol_coarse, @@ -160,7 +160,7 @@ void CSingleGridIntegration::SetRestricted_EddyVisc(unsigned short RunTime_EqSys /*--- MPI the new interpolated solution (this also includes the eddy viscosity) ---*/ - sol_coarse->InitiateComms(geo_coarse, config, SOLUTION_EDDY); - sol_coarse->CompleteComms(geo_coarse, config, SOLUTION_EDDY); + sol_coarse->InitiateComms(geo_coarse, config, MPI_QUANTITIES::SOLUTION_EDDY); + sol_coarse->CompleteComms(geo_coarse, config, MPI_QUANTITIES::SOLUTION_EDDY); } diff --git a/SU2_CFD/src/iteration/CDiscAdjFEAIteration.cpp b/SU2_CFD/src/iteration/CDiscAdjFEAIteration.cpp index d7848499f62..53f1bff6945 100644 --- a/SU2_CFD/src/iteration/CDiscAdjFEAIteration.cpp +++ b/SU2_CFD/src/iteration/CDiscAdjFEAIteration.cpp @@ -233,11 +233,11 @@ void CDiscAdjFEAIteration::SetDependencies(CSolver***** solver, CGeometry**** ge /*--- MPI dependencies. ---*/ - dir_solver->InitiateComms(structural_geometry, config[iZone], SOLUTION_FEA); - dir_solver->CompleteComms(structural_geometry, config[iZone], SOLUTION_FEA); + dir_solver->InitiateComms(structural_geometry, config[iZone], MPI_QUANTITIES::SOLUTION_FEA); + dir_solver->CompleteComms(structural_geometry, config[iZone], MPI_QUANTITIES::SOLUTION_FEA); - structural_geometry->InitiateComms(structural_geometry, config[iZone], COORDINATES); - structural_geometry->CompleteComms(structural_geometry, config[iZone], COORDINATES); + structural_geometry->InitiateComms(structural_geometry, config[iZone], MPI_QUANTITIES::COORDINATES); + structural_geometry->CompleteComms(structural_geometry, config[iZone], MPI_QUANTITIES::COORDINATES); } END_SU2_OMP_PARALLEL diff --git a/SU2_CFD/src/iteration/CDiscAdjFluidIteration.cpp b/SU2_CFD/src/iteration/CDiscAdjFluidIteration.cpp index 6bde7b99418..27173f3357e 100644 --- a/SU2_CFD/src/iteration/CDiscAdjFluidIteration.cpp +++ b/SU2_CFD/src/iteration/CDiscAdjFluidIteration.cpp @@ -468,26 +468,26 @@ void CDiscAdjFluidIteration::SetDependencies(CSolver***** solver, CGeometry**** /*--- Compute coupling between flow, turbulent and species equations ---*/ solvers0[FLOW_SOL]->Preprocessing(geometry0, solvers0, config[iZone], MESH_0, NO_RK_ITER, RUNTIME_FLOW_SYS, true); - solvers0[FLOW_SOL]->InitiateComms(geometry0, config[iZone], SOLUTION); - solvers0[FLOW_SOL]->CompleteComms(geometry0, config[iZone], SOLUTION); + solvers0[FLOW_SOL]->InitiateComms(geometry0, config[iZone], MPI_QUANTITIES::SOLUTION); + solvers0[FLOW_SOL]->CompleteComms(geometry0, config[iZone], MPI_QUANTITIES::SOLUTION); if (turbulent && !config[iZone]->GetFrozen_Visc_Disc()) { solvers0[TURB_SOL]->Postprocessing(geometry0, solvers0, config[iZone], MESH_0); - solvers0[TURB_SOL]->InitiateComms(geometry0, config[iZone], SOLUTION); - solvers0[TURB_SOL]->CompleteComms(geometry0, config[iZone], SOLUTION); + solvers0[TURB_SOL]->InitiateComms(geometry0, config[iZone], MPI_QUANTITIES::SOLUTION); + solvers0[TURB_SOL]->CompleteComms(geometry0, config[iZone], MPI_QUANTITIES::SOLUTION); } if (config[iZone]->GetKind_Species_Model() != SPECIES_MODEL::NONE) { solvers0[SPECIES_SOL]->Preprocessing(geometry0, solvers0, config[iZone], MESH_0, NO_RK_ITER, RUNTIME_FLOW_SYS, true); - solvers0[SPECIES_SOL]->InitiateComms(geometry0, config[iZone], SOLUTION); - solvers0[SPECIES_SOL]->CompleteComms(geometry0, config[iZone], SOLUTION); + solvers0[SPECIES_SOL]->InitiateComms(geometry0, config[iZone], MPI_QUANTITIES::SOLUTION); + solvers0[SPECIES_SOL]->CompleteComms(geometry0, config[iZone], MPI_QUANTITIES::SOLUTION); } if (config[iZone]->GetWeakly_Coupled_Heat()) { solvers0[HEAT_SOL]->Set_Heatflux_Areas(geometry0, config[iZone]); solvers0[HEAT_SOL]->Preprocessing(geometry0, solvers0, config[iZone], MESH_0, NO_RK_ITER, RUNTIME_HEAT_SYS, true); solvers0[HEAT_SOL]->Postprocessing(geometry0, solvers0, config[iZone], MESH_0); - solvers0[HEAT_SOL]->InitiateComms(geometry0, config[iZone], SOLUTION); - solvers0[HEAT_SOL]->CompleteComms(geometry0, config[iZone], SOLUTION); + solvers0[HEAT_SOL]->InitiateComms(geometry0, config[iZone], MPI_QUANTITIES::SOLUTION); + solvers0[HEAT_SOL]->CompleteComms(geometry0, config[iZone], MPI_QUANTITIES::SOLUTION); } } @@ -495,8 +495,8 @@ void CDiscAdjFluidIteration::SetDependencies(CSolver***** solver, CGeometry**** if (config[iZone]->AddRadiation()) { solvers0[RAD_SOL]->Postprocessing(geometry0, solvers0, config[iZone], MESH_0); - solvers0[RAD_SOL]->InitiateComms(geometry0, config[iZone], SOLUTION); - solvers0[RAD_SOL]->CompleteComms(geometry0, config[iZone], SOLUTION); + solvers0[RAD_SOL]->InitiateComms(geometry0, config[iZone], MPI_QUANTITIES::SOLUTION); + solvers0[RAD_SOL]->CompleteComms(geometry0, config[iZone], MPI_QUANTITIES::SOLUTION); } } diff --git a/SU2_CFD/src/iteration/CDiscAdjHeatIteration.cpp b/SU2_CFD/src/iteration/CDiscAdjHeatIteration.cpp index 1c3666d5216..f0c14d4168a 100644 --- a/SU2_CFD/src/iteration/CDiscAdjHeatIteration.cpp +++ b/SU2_CFD/src/iteration/CDiscAdjHeatIteration.cpp @@ -213,8 +213,8 @@ void CDiscAdjHeatIteration::SetDependencies(CSolver***** solver, CGeometry**** g solvers0[HEAT_SOL]->Preprocessing(geometries[MESH_0], solvers0, config[iZone], MESH_0, NO_RK_ITER, RUNTIME_HEAT_SYS, true); solvers0[HEAT_SOL]->Postprocessing(geometries[MESH_0], solvers0, config[iZone], MESH_0); - solvers0[HEAT_SOL]->InitiateComms(geometries[MESH_0], config[iZone], SOLUTION); - solvers0[HEAT_SOL]->CompleteComms(geometries[MESH_0], config[iZone], SOLUTION); + solvers0[HEAT_SOL]->InitiateComms(geometries[MESH_0], config[iZone], MPI_QUANTITIES::SOLUTION); + solvers0[HEAT_SOL]->CompleteComms(geometries[MESH_0], config[iZone], MPI_QUANTITIES::SOLUTION); } void CDiscAdjHeatIteration::RegisterOutput(CSolver***** solver, CGeometry**** geometry, CConfig** config, diff --git a/SU2_CFD/src/solvers/CAdjEulerSolver.cpp b/SU2_CFD/src/solvers/CAdjEulerSolver.cpp index e43bd91c871..17128c7fec8 100644 --- a/SU2_CFD/src/solvers/CAdjEulerSolver.cpp +++ b/SU2_CFD/src/solvers/CAdjEulerSolver.cpp @@ -328,8 +328,8 @@ CAdjEulerSolver::CAdjEulerSolver(CGeometry *geometry, CConfig *config, unsigned /*--- MPI solution ---*/ - InitiateComms(geometry, config, SOLUTION); - CompleteComms(geometry, config, SOLUTION); + InitiateComms(geometry, config, MPI_QUANTITIES::SOLUTION); + CompleteComms(geometry, config, MPI_QUANTITIES::SOLUTION); SolverName = "ADJ.FLOW"; } @@ -888,8 +888,8 @@ void CAdjEulerSolver::SetInitialCondition(CGeometry **geometry, CSolver ***solve for (auto iMesh = 1ul; iMesh <= config->GetnMGLevels(); iMesh++) { MultigridRestriction(*geometry[iMesh - 1], solver_container[iMesh - 1][ADJFLOW_SOL]->GetNodes()->GetSolution(), *geometry[iMesh], solver_container[iMesh][ADJFLOW_SOL]->GetNodes()->GetSolution()); - solver_container[iMesh][ADJFLOW_SOL]->InitiateComms(geometry[iMesh], config, SOLUTION); - solver_container[iMesh][ADJFLOW_SOL]->CompleteComms(geometry[iMesh], config, SOLUTION); + solver_container[iMesh][ADJFLOW_SOL]->InitiateComms(geometry[iMesh], config, MPI_QUANTITIES::SOLUTION); + solver_container[iMesh][ADJFLOW_SOL]->CompleteComms(geometry[iMesh], config, MPI_QUANTITIES::SOLUTION); } } @@ -1343,8 +1343,8 @@ void CAdjEulerSolver::SetCentered_Dissipation_Sensor(CGeometry *geometry, CConfi /*--- MPI parallelization ---*/ - InitiateComms(geometry, config, SENSOR); - CompleteComms(geometry, config, SENSOR); + InitiateComms(geometry, config, MPI_QUANTITIES::SENSOR); + CompleteComms(geometry, config, MPI_QUANTITIES::SENSOR); } @@ -1377,8 +1377,8 @@ void CAdjEulerSolver::ExplicitRK_Iteration(CGeometry *geometry, CSolver **solver /*--- MPI solution ---*/ - InitiateComms(geometry, config, SOLUTION); - CompleteComms(geometry, config, SOLUTION); + InitiateComms(geometry, config, MPI_QUANTITIES::SOLUTION); + CompleteComms(geometry, config, MPI_QUANTITIES::SOLUTION); /*--- Compute the root mean square residual ---*/ SetResidual_RMS(geometry, config); @@ -1411,8 +1411,8 @@ void CAdjEulerSolver::ExplicitEuler_Iteration(CGeometry *geometry, CSolver **sol /*--- MPI solution ---*/ - InitiateComms(geometry, config, SOLUTION); - CompleteComms(geometry, config, SOLUTION); + InitiateComms(geometry, config, MPI_QUANTITIES::SOLUTION); + CompleteComms(geometry, config, MPI_QUANTITIES::SOLUTION); /*--- Compute the root mean square residual ---*/ SetResidual_RMS(geometry, config); @@ -1491,8 +1491,8 @@ void CAdjEulerSolver::ImplicitEuler_Iteration(CGeometry *geometry, CSolver **sol /*--- MPI solution ---*/ - InitiateComms(geometry, config, SOLUTION); - CompleteComms(geometry, config, SOLUTION); + InitiateComms(geometry, config, MPI_QUANTITIES::SOLUTION); + CompleteComms(geometry, config, MPI_QUANTITIES::SOLUTION); /*--- Compute the root mean square residual ---*/ @@ -3880,8 +3880,8 @@ void CAdjEulerSolver::LoadRestart(CGeometry **geometry, CSolver ***solver, CConf it down to the coarse levels. We also call the preprocessing routine on the fine level in order to have all necessary quantities updated. ---*/ - solver[MESH_0][ADJFLOW_SOL]->InitiateComms(geometry[MESH_0], config, SOLUTION); - solver[MESH_0][ADJFLOW_SOL]->CompleteComms(geometry[MESH_0], config, SOLUTION); + solver[MESH_0][ADJFLOW_SOL]->InitiateComms(geometry[MESH_0], config, MPI_QUANTITIES::SOLUTION); + solver[MESH_0][ADJFLOW_SOL]->CompleteComms(geometry[MESH_0], config, MPI_QUANTITIES::SOLUTION); solver[MESH_0][ADJFLOW_SOL]->Preprocessing(geometry[MESH_0], solver[MESH_0], config, MESH_0, NO_RK_ITER, RUNTIME_FLOW_SYS, false); /*--- Interpolate the solution down to the coarse multigrid levels ---*/ @@ -3889,8 +3889,8 @@ void CAdjEulerSolver::LoadRestart(CGeometry **geometry, CSolver ***solver, CConf for (iMesh = 1; iMesh <= config->GetnMGLevels(); iMesh++) { MultigridRestriction(*geometry[iMesh - 1], solver[iMesh - 1][ADJFLOW_SOL]->GetNodes()->GetSolution(), *geometry[iMesh], solver[iMesh][ADJFLOW_SOL]->GetNodes()->GetSolution()); - solver[iMesh][ADJFLOW_SOL]->InitiateComms(geometry[iMesh], config, SOLUTION); - solver[iMesh][ADJFLOW_SOL]->CompleteComms(geometry[iMesh], config, SOLUTION); + solver[iMesh][ADJFLOW_SOL]->InitiateComms(geometry[iMesh], config, MPI_QUANTITIES::SOLUTION); + solver[iMesh][ADJFLOW_SOL]->CompleteComms(geometry[iMesh], config, MPI_QUANTITIES::SOLUTION); solver[iMesh][ADJFLOW_SOL]->Preprocessing(geometry[iMesh], solver[iMesh], config, iMesh, NO_RK_ITER, RUNTIME_FLOW_SYS, false); } diff --git a/SU2_CFD/src/solvers/CAdjNSSolver.cpp b/SU2_CFD/src/solvers/CAdjNSSolver.cpp index b9da5ad2b63..c07425a2bdc 100644 --- a/SU2_CFD/src/solvers/CAdjNSSolver.cpp +++ b/SU2_CFD/src/solvers/CAdjNSSolver.cpp @@ -281,8 +281,8 @@ CAdjNSSolver::CAdjNSSolver(CGeometry *geometry, CConfig *config, unsigned short /*--- MPI solution ---*/ - InitiateComms(geometry, config, SOLUTION); - CompleteComms(geometry, config, SOLUTION); + InitiateComms(geometry, config, MPI_QUANTITIES::SOLUTION); + CompleteComms(geometry, config, MPI_QUANTITIES::SOLUTION); } diff --git a/SU2_CFD/src/solvers/CAdjTurbSolver.cpp b/SU2_CFD/src/solvers/CAdjTurbSolver.cpp index 58403f541a4..3d5b6eb7024 100644 --- a/SU2_CFD/src/solvers/CAdjTurbSolver.cpp +++ b/SU2_CFD/src/solvers/CAdjTurbSolver.cpp @@ -158,8 +158,8 @@ CAdjTurbSolver::CAdjTurbSolver(CGeometry *geometry, CConfig *config, unsigned sh /*--- MPI solution ---*/ - InitiateComms(geometry, config, SOLUTION); - CompleteComms(geometry, config, SOLUTION); + InitiateComms(geometry, config, MPI_QUANTITIES::SOLUTION); + CompleteComms(geometry, config, MPI_QUANTITIES::SOLUTION); } @@ -511,8 +511,8 @@ void CAdjTurbSolver::ImplicitEuler_Iteration(CGeometry *geometry, CSolver **solv /*--- MPI solution ---*/ - InitiateComms(geometry, config, SOLUTION); - CompleteComms(geometry, config, SOLUTION); + InitiateComms(geometry, config, MPI_QUANTITIES::SOLUTION); + CompleteComms(geometry, config, MPI_QUANTITIES::SOLUTION); /*--- Compute the root mean square residual ---*/ diff --git a/SU2_CFD/src/solvers/CBaselineSolver.cpp b/SU2_CFD/src/solvers/CBaselineSolver.cpp index bc8b9c1ebc5..411cce58397 100644 --- a/SU2_CFD/src/solvers/CBaselineSolver.cpp +++ b/SU2_CFD/src/solvers/CBaselineSolver.cpp @@ -459,8 +459,8 @@ void CBaselineSolver::LoadRestart(CGeometry **geometry, CSolver ***solver, CConf /*--- MPI solution ---*/ - InitiateComms(geometry[iInst], config, SOLUTION); - CompleteComms(geometry[iInst], config, SOLUTION); + InitiateComms(geometry[iInst], config, MPI_QUANTITIES::SOLUTION); + CompleteComms(geometry[iInst], config, MPI_QUANTITIES::SOLUTION); /*--- Update the geometry for flows on dynamic meshes ---*/ @@ -468,11 +468,11 @@ void CBaselineSolver::LoadRestart(CGeometry **geometry, CSolver ***solver, CConf /*--- Communicate the new coordinates and grid velocities at the halos ---*/ - geometry[iInst]->InitiateComms(geometry[iInst], config, COORDINATES); - geometry[iInst]->CompleteComms(geometry[iInst], config, COORDINATES); + geometry[iInst]->InitiateComms(geometry[iInst], config, MPI_QUANTITIES::COORDINATES); + geometry[iInst]->CompleteComms(geometry[iInst], config, MPI_QUANTITIES::COORDINATES); - geometry[iInst]->InitiateComms(geometry[iInst], config, GRID_VELOCITY); - geometry[iInst]->CompleteComms(geometry[iInst], config, GRID_VELOCITY); + geometry[iInst]->InitiateComms(geometry[iInst], config, MPI_QUANTITIES::GRID_VELOCITY); + geometry[iInst]->CompleteComms(geometry[iInst], config, MPI_QUANTITIES::GRID_VELOCITY); } diff --git a/SU2_CFD/src/solvers/CEulerSolver.cpp b/SU2_CFD/src/solvers/CEulerSolver.cpp index 16cdfe59a4b..bdc4db5a473 100644 --- a/SU2_CFD/src/solvers/CEulerSolver.cpp +++ b/SU2_CFD/src/solvers/CEulerSolver.cpp @@ -2351,8 +2351,8 @@ void CEulerSolver::SetUndivided_Laplacian(CGeometry *geometry, const CConfig *co /*--- MPI parallelization ---*/ - InitiateComms(geometry, config, UNDIVIDED_LAPLACIAN); - CompleteComms(geometry, config, UNDIVIDED_LAPLACIAN); + InitiateComms(geometry, config, MPI_QUANTITIES::UNDIVIDED_LAPLACIAN); + CompleteComms(geometry, config, MPI_QUANTITIES::UNDIVIDED_LAPLACIAN); } @@ -2415,8 +2415,8 @@ void CEulerSolver::SetUpwind_Ducros_Sensor(CGeometry *geometry, CConfig *config) } END_SU2_OMP_FOR - InitiateComms(geometry, config, SENSOR); - CompleteComms(geometry, config, SENSOR); + InitiateComms(geometry, config, MPI_QUANTITIES::SENSOR); + CompleteComms(geometry, config, MPI_QUANTITIES::SENSOR); } @@ -5067,8 +5067,8 @@ void CEulerSolver::BC_Riemann(CGeometry *geometry, CSolver **solver_container, su2double Velocity_i[MAXNDIM]={0}; for (auto iDim=0u; iDim < nDim; iDim++) Velocity_i[iDim] = nodes->GetVelocity(iPoint,iDim); - - const auto Velocity2_i = GeometryToolbox::SquaredNorm(nDim, Velocity_i); + + const auto Velocity2_i = GeometryToolbox::SquaredNorm(nDim, Velocity_i); const auto Density_i = nodes->GetDensity(iPoint); @@ -5090,7 +5090,7 @@ void CEulerSolver::BC_Riemann(CGeometry *geometry, CSolver **solver_container, T_Total{0}, P_Total{0}, Density_e{0}, StaticEnthalpy_e{0}, StaticEnergy_e{0}; su2double Velocity2_e{0}, NormalVelocity{0}, TangVelocity{0}, VelMag_e{0}; - su2double Velocity_e[MAXNDIM] = {0}; + su2double Velocity_e[MAXNDIM] = {0}; const su2double * Flow_Dir, * Mach; switch(config->GetKind_Data_Riemann(Marker_Tag)) { @@ -5146,10 +5146,10 @@ void CEulerSolver::BC_Riemann(CGeometry *geometry, CSolver **solver_container, GetFluidModel()->SetTDState_PT(P_static, T_static); /* --- Compute the boundary state u_e --- */ - for (auto iDim = 0u; iDim < nDim; iDim++) + for (auto iDim = 0u; iDim < nDim; iDim++) Velocity_e[iDim] = Mach[iDim]*GetFluidModel()->GetSoundSpeed(); - Velocity2_e = GeometryToolbox::SquaredNorm(nDim, Velocity_e); + Velocity2_e = GeometryToolbox::SquaredNorm(nDim, Velocity_e); Density_e = GetFluidModel()->GetDensity(); StaticEnergy_e = GetFluidModel()->GetStaticEnergy(); Energy_e = StaticEnergy_e + 0.5 * Velocity2_e; @@ -5176,7 +5176,7 @@ void CEulerSolver::BC_Riemann(CGeometry *geometry, CSolver **solver_container, for (auto iDim = 0u; iDim < nDim; iDim++) Velocity_e[iDim] = Mach[iDim]*GetFluidModel()->GetSoundSpeed(); - Velocity2_e = GeometryToolbox::SquaredNorm(nDim, Velocity_e); + Velocity2_e = GeometryToolbox::SquaredNorm(nDim, Velocity_e); Density_e = GetFluidModel()->GetDensity(); StaticEnergy_e = GetFluidModel()->GetStaticEnergy(); Energy_e = StaticEnergy_e + 0.5 * Velocity2_e; @@ -5208,10 +5208,10 @@ void CEulerSolver::BC_Riemann(CGeometry *geometry, CSolver **solver_container, /* --- Compute the boundary state u_e --- */ GetFluidModel()->SetTDState_Prho(Pressure_e, Density_e); - for (auto iDim = 0u; iDim < nDim; iDim++) + for (auto iDim = 0u; iDim < nDim; iDim++) Velocity_e[iDim] = Velocity_i[iDim]; - Velocity2_e = GeometryToolbox::SquaredNorm(nDim, Velocity_e); + Velocity2_e = GeometryToolbox::SquaredNorm(nDim, Velocity_e); Energy_e = GetFluidModel()->GetStaticEnergy() + 0.5*Velocity2_e; break; @@ -5235,7 +5235,7 @@ void CEulerSolver::BC_Riemann(CGeometry *geometry, CSolver **solver_container, } /*--- Flow eigenvalues, boundary state u_e and u_i ---*/ - su2double Lambda_i[MAXNVAR] = {0}, u_e[MAXNVAR] = {0}, u_i[MAXNVAR]={0}, u_b[MAXNVAR]={0}, + su2double Lambda_i[MAXNVAR] = {0}, u_e[MAXNVAR] = {0}, u_i[MAXNVAR]={0}, u_b[MAXNVAR]={0}, dw[MAXNVAR]={0}; u_e[0] = Density_e; u_i[0] = Density_i; @@ -5302,7 +5302,7 @@ void CEulerSolver::BC_Riemann(CGeometry *geometry, CSolver **solver_container, LinSysRes.AddBlock(iPoint, Residual); if (implicit) { - su2double **Jacobian_b = new su2double*[nVar], + su2double **Jacobian_b = new su2double*[nVar], **Jacobian_i = new su2double*[nVar]; su2double **DubDu = new su2double*[nVar]; for (auto iVar = 0u; iVar < nVar; iVar++){ @@ -5333,7 +5333,7 @@ void CEulerSolver::BC_Riemann(CGeometry *geometry, CSolver **solver_container, if (dynamic_grid){ const auto gridVel = geometry->nodes->GetGridVel(iPoint); const auto projVelocity = GeometryToolbox::DotProduct(nDim, gridVel, Normal); - for (auto iVar = 0u; iVar < nVar; iVar++) + for (auto iVar = 0u; iVar < nVar; iVar++) Jacobian_b[iVar][iVar] -= projVelocity; } @@ -5394,7 +5394,7 @@ void CEulerSolver::BC_Riemann(CGeometry *geometry, CSolver **solver_container, visc_numerics->SetPrimVarGradient(nodes->GetGradient_Primitive(iPoint), nodes->GetGradient_Primitive(iPoint)); /*--- Secondary variables ---*/ - + auto S_domain = nodes->GetSecondary(iPoint); /*--- Compute secondary thermodynamic properties (partial derivatives...) ---*/ @@ -8795,9 +8795,9 @@ void CEulerSolver::PreprocessAverage(CSolver **solver, CGeometry *geometry, CCon const auto nSpanWiseSections = config->GetnSpanWiseSections(); const auto iZone = config->GetiZone(); - + for (auto iSpan= 0u; iSpan < nSpanWiseSections; iSpan++){ - su2double TotalAreaVelocity[MAXNDIM]={0.0}, + su2double TotalAreaVelocity[MAXNDIM]={0.0}, TotalAreaPressure{0}, TotalAreaDensity{0}; for (auto iMarker = 0u; iMarker < config->GetnMarker_All(); iMarker++){ @@ -8987,7 +8987,7 @@ void CEulerSolver::TurboAverageProcess(CSolver **solver, CGeometry *geometry, CC TotalAreaVelocity[iDim] += Area*Velocity[iDim]; TotalMassVelocity[iDim] += Area*(Density*TurboVelocity[0] )*Velocity[iDim]; } - + TotalFluxes[0] += Area*(Density*TurboVelocity[0]); TotalFluxes[1] += Area*(Density*TurboVelocity[0]*TurboVelocity[0] + Pressure); for (auto iDim = 2; iDim < nDim+1; iDim++) @@ -9036,7 +9036,7 @@ void CEulerSolver::TurboAverageProcess(CSolver **solver, CGeometry *geometry, CC UpdateTotalQuantities(iMarker, jSpan, iVertex); } } - } + } } // marker_flag match } // iMarkerTP match @@ -9251,7 +9251,7 @@ void CEulerSolver::TurboAverageProcess(CSolver **solver, CGeometry *geometry, CC OmegaOut[iMarkerTP - 1][iSpan] = AverageOmega[iMarker][iSpan]; NuOut[iMarkerTP - 1][iSpan] = AverageNu[iMarker][iSpan]; } - + auto TurboVel = (marker_flag == INFLOW) ? TurboVelocityIn[iMarkerTP - 1][iSpan] : TurboVelocityOut[iMarkerTP - 1][iSpan]; if (performance_average_process == MIXEDOUT) { @@ -9281,7 +9281,7 @@ void CEulerSolver::TurboAverageProcess(CSolver **solver, CGeometry *geometry, CC if(config->GetKind_Data_Giles(Marker_Tag) == RADIAL_EQUILIBRIUM){ RadialEquilibriumPressure[iMarker][nSpanWiseSections/2] = config->GetGiles_Var1(Marker_Tag)/config->GetPressure_Ref(); } - } + } for (auto iSpan= nSpanWiseSections/2; iSpan < nSpanWiseSections-1; iSpan++){ const auto Radius2 = geometry->GetTurboRadius(iMarker,iSpan+1); const auto Radius1 = geometry->GetTurboRadius(iMarker,iSpan); @@ -9325,9 +9325,9 @@ void CEulerSolver::MixedOut_Average(CConfig *config, su2double val_init_pressure su2double vel[MAXNDIM] = {0}; vel[0] = (val_Averaged_Flux[1] - pressure_mix) / val_Averaged_Flux[0]; - for (auto iDim = 1u; iDim < nDim; iDim++) + for (auto iDim = 1u; iDim < nDim; iDim++) vel[iDim] = val_Averaged_Flux[iDim+1] / val_Averaged_Flux[0]; - + const su2double velsq = GeometryToolbox::DotProduct(nDim, vel, vel); diff --git a/SU2_CFD/src/solvers/CFEASolver.cpp b/SU2_CFD/src/solvers/CFEASolver.cpp index f67a0199d01..fd142988bce 100644 --- a/SU2_CFD/src/solvers/CFEASolver.cpp +++ b/SU2_CFD/src/solvers/CFEASolver.cpp @@ -3116,8 +3116,8 @@ void CFEASolver::LoadRestart(CGeometry **geometry, CSolver ***solver, CConfig *c /*--- MPI. If dynamic, we also need to communicate the old solution. ---*/ - InitiateComms(geometry[MESH_0], config, SOLUTION_FEA); - CompleteComms(geometry[MESH_0], config, SOLUTION_FEA); + InitiateComms(geometry[MESH_0], config, MPI_QUANTITIES::SOLUTION_FEA); + CompleteComms(geometry[MESH_0], config, MPI_QUANTITIES::SOLUTION_FEA); /*--- It's important to not push back the solution when this function is used to load solutions for * unsteady discrete adjoints, otherwise we overwrite one of the two solutions needed. ---*/ diff --git a/SU2_CFD/src/solvers/CGradientSmoothingSolver.cpp b/SU2_CFD/src/solvers/CGradientSmoothingSolver.cpp index 1f3a425aae6..47a624faada 100644 --- a/SU2_CFD/src/solvers/CGradientSmoothingSolver.cpp +++ b/SU2_CFD/src/solvers/CGradientSmoothingSolver.cpp @@ -287,8 +287,8 @@ void CGradientSmoothingSolver::ApplyGradientSmoothingDV(CGeometry* geometry, CNu /*--- Matrix vector product with the Laplace-Beltrami stiffness matrix. ---*/ if (config->GetSmoothOnSurface()) { - CSysMatrixComms::Initiate(helperVecIn, geometry, config, SOLUTION_MATRIX); - CSysMatrixComms::Complete(helperVecIn, geometry, config, SOLUTION_MATRIX); + CSysMatrixComms::Initiate(helperVecIn, geometry, config, MPI_QUANTITIES::SOLUTION_MATRIX); + CSysMatrixComms::Complete(helperVecIn, geometry, config, MPI_QUANTITIES::SOLUTION_MATRIX); mat_vec(helperVecIn, helperVecAux); @@ -306,8 +306,8 @@ void CGradientSmoothingSolver::ApplyGradientSmoothingDV(CGeometry* geometry, CNu grid_movement->SetVolume_Deformation(geometry, config, false, true, true); CGradientSmoothingSolverDetails::ReadVectorToGeometry(geometry, helperVecIn); - CSysMatrixComms::Initiate(helperVecIn, geometry, config, SOLUTION_MATRIX); - CSysMatrixComms::Complete(helperVecIn, geometry, config, SOLUTION_MATRIX); + CSysMatrixComms::Initiate(helperVecIn, geometry, config, MPI_QUANTITIES::SOLUTION_MATRIX); + CSysMatrixComms::Complete(helperVecIn, geometry, config, MPI_QUANTITIES::SOLUTION_MATRIX); mat_vec(helperVecIn, helperVecAux); diff --git a/SU2_CFD/src/solvers/CHeatSolver.cpp b/SU2_CFD/src/solvers/CHeatSolver.cpp index b01870e29bb..64ddb898e74 100644 --- a/SU2_CFD/src/solvers/CHeatSolver.cpp +++ b/SU2_CFD/src/solvers/CHeatSolver.cpp @@ -152,8 +152,8 @@ CHeatSolver::CHeatSolver(CGeometry *geometry, CConfig *config, unsigned short iM /*--- MPI solution ---*/ - InitiateComms(geometry, config, SOLUTION); - CompleteComms(geometry, config, SOLUTION); + InitiateComms(geometry, config, MPI_QUANTITIES::SOLUTION); + CompleteComms(geometry, config, MPI_QUANTITIES::SOLUTION); /*--- Store the initial CFL number for all grid points. ---*/ @@ -246,8 +246,8 @@ void CHeatSolver::LoadRestart(CGeometry **geometry, CSolver ***solver, CConfig * on the fine level in order to have all necessary quantities updated, especially if this is a turbulent simulation (eddy viscosity). ---*/ - solver[MESH_0][HEAT_SOL]->InitiateComms(geometry[MESH_0], config, SOLUTION); - solver[MESH_0][HEAT_SOL]->CompleteComms(geometry[MESH_0], config, SOLUTION); + solver[MESH_0][HEAT_SOL]->InitiateComms(geometry[MESH_0], config, MPI_QUANTITIES::SOLUTION); + solver[MESH_0][HEAT_SOL]->CompleteComms(geometry[MESH_0], config, MPI_QUANTITIES::SOLUTION); solver[MESH_0][HEAT_SOL]->Preprocessing(geometry[MESH_0], solver[MESH_0], config, MESH_0, NO_RK_ITER, RUNTIME_HEAT_SYS, false); @@ -256,8 +256,8 @@ void CHeatSolver::LoadRestart(CGeometry **geometry, CSolver ***solver, CConfig * for (auto iMesh = 1u; iMesh <= config->GetnMGLevels(); iMesh++) { MultigridRestriction(*geometry[iMesh - 1], solver[iMesh - 1][HEAT_SOL]->GetNodes()->GetSolution(), *geometry[iMesh], solver[iMesh][HEAT_SOL]->GetNodes()->GetSolution()); - solver[iMesh][HEAT_SOL]->InitiateComms(geometry[iMesh], config, SOLUTION); - solver[iMesh][HEAT_SOL]->CompleteComms(geometry[iMesh], config, SOLUTION); + solver[iMesh][HEAT_SOL]->InitiateComms(geometry[iMesh], config, MPI_QUANTITIES::SOLUTION); + solver[iMesh][HEAT_SOL]->CompleteComms(geometry[iMesh], config, MPI_QUANTITIES::SOLUTION); solver[iMesh][HEAT_SOL]->Preprocessing(geometry[iMesh], solver[iMesh], config, iMesh, NO_RK_ITER, RUNTIME_HEAT_SYS, false); } @@ -934,8 +934,8 @@ void CHeatSolver::SetInitialCondition(CGeometry **geometry, CSolver ***solver_co for (auto iMesh = 1u; iMesh <= config->GetnMGLevels(); iMesh++) { MultigridRestriction(*geometry[iMesh - 1], solver_container[iMesh - 1][HEAT_SOL]->GetNodes()->GetSolution(), *geometry[iMesh], solver_container[iMesh][HEAT_SOL]->GetNodes()->GetSolution()); - solver_container[iMesh][HEAT_SOL]->InitiateComms(geometry[iMesh], config, SOLUTION); - solver_container[iMesh][HEAT_SOL]->CompleteComms(geometry[iMesh], config, SOLUTION); + solver_container[iMesh][HEAT_SOL]->InitiateComms(geometry[iMesh], config, MPI_QUANTITIES::SOLUTION); + solver_container[iMesh][HEAT_SOL]->CompleteComms(geometry[iMesh], config, MPI_QUANTITIES::SOLUTION); } } diff --git a/SU2_CFD/src/solvers/CMeshSolver.cpp b/SU2_CFD/src/solvers/CMeshSolver.cpp index 39e17b322c8..e875dcab694 100644 --- a/SU2_CFD/src/solvers/CMeshSolver.cpp +++ b/SU2_CFD/src/solvers/CMeshSolver.cpp @@ -479,14 +479,14 @@ void CMeshSolver::DeformMesh(CGeometry **geometry, CNumerics **numerics, CConfig if (multizone) nodes->Set_BGSSolution_k(); /*--- Capture a few MPI dependencies for AD. ---*/ - geometry[MESH_0]->InitiateComms(geometry[MESH_0], config, COORDINATES); - geometry[MESH_0]->CompleteComms(geometry[MESH_0], config, COORDINATES); + geometry[MESH_0]->InitiateComms(geometry[MESH_0], config, MPI_QUANTITIES::COORDINATES); + geometry[MESH_0]->CompleteComms(geometry[MESH_0], config, MPI_QUANTITIES::COORDINATES); - InitiateComms(geometry[MESH_0], config, SOLUTION); - CompleteComms(geometry[MESH_0], config, SOLUTION); + InitiateComms(geometry[MESH_0], config, MPI_QUANTITIES::SOLUTION); + CompleteComms(geometry[MESH_0], config, MPI_QUANTITIES::SOLUTION); - InitiateComms(geometry[MESH_0], config, MESH_DISPLACEMENTS); - CompleteComms(geometry[MESH_0], config, MESH_DISPLACEMENTS); + InitiateComms(geometry[MESH_0], config, MPI_QUANTITIES::MESH_DISPLACEMENTS); + CompleteComms(geometry[MESH_0], config, MPI_QUANTITIES::MESH_DISPLACEMENTS); /*--- Compute the stiffness matrix, no point recording because we clear the residual. ---*/ @@ -562,8 +562,8 @@ void CMeshSolver::UpdateGridCoord(CGeometry *geometry, const CConfig *config){ END_SU2_OMP_FOR /*--- Communicate the updated displacements and mesh coordinates. ---*/ - geometry->InitiateComms(geometry, config, COORDINATES); - geometry->CompleteComms(geometry, config, COORDINATES); + geometry->InitiateComms(geometry, config, MPI_QUANTITIES::COORDINATES); + geometry->CompleteComms(geometry, config, MPI_QUANTITIES::COORDINATES); } @@ -818,8 +818,8 @@ void CMeshSolver::LoadRestart(CGeometry **geometry, CSolver ***solver, CConfig * } /*--- Communicate the loaded displacements. ---*/ - solver[MESH_0][MESH_SOL]->InitiateComms(geometry[MESH_0], config, SOLUTION); - solver[MESH_0][MESH_SOL]->CompleteComms(geometry[MESH_0], config, SOLUTION); + solver[MESH_0][MESH_SOL]->InitiateComms(geometry[MESH_0], config, MPI_QUANTITIES::SOLUTION); + solver[MESH_0][MESH_SOL]->CompleteComms(geometry[MESH_0], config, MPI_QUANTITIES::SOLUTION); /*--- Init the linear system solution. ---*/ for (unsigned long iPoint = 0; iPoint < nPoint; ++iPoint) { @@ -880,7 +880,7 @@ void CMeshSolver::RestartOldGeometry(CGeometry *geometry, const CConfig *config) for(unsigned short iStep = 1; iStep <= nSteps; ++iStep) { - unsigned short CommType = (iStep == 1) ? SOLUTION_TIME_N : SOLUTION_TIME_N1; + MPI_QUANTITIES CommType = (iStep == 1) ? MPI_QUANTITIES::SOLUTION_TIME_N : MPI_QUANTITIES::SOLUTION_TIME_N1; /*--- Modify file name for an unsteady restart ---*/ int Unst_RestartIter; diff --git a/SU2_CFD/src/solvers/CNEMOEulerSolver.cpp b/SU2_CFD/src/solvers/CNEMOEulerSolver.cpp index 0905ddeda32..3519d6d912d 100644 --- a/SU2_CFD/src/solvers/CNEMOEulerSolver.cpp +++ b/SU2_CFD/src/solvers/CNEMOEulerSolver.cpp @@ -2354,7 +2354,7 @@ void CNEMOEulerSolver::SetPressureDiffusionSensor(CGeometry *geometry, CConfig * /*--- MPI parallelization ---*/ - InitiateComms(geometry, config, SENSOR); - CompleteComms(geometry, config, SENSOR); + InitiateComms(geometry, config, MPI_QUANTITIES::SENSOR); + CompleteComms(geometry, config, MPI_QUANTITIES::SENSOR); } diff --git a/SU2_CFD/src/solvers/CNEMONSSolver.cpp b/SU2_CFD/src/solvers/CNEMONSSolver.cpp index b6642b5cf8d..b30b99edebe 100644 --- a/SU2_CFD/src/solvers/CNEMONSSolver.cpp +++ b/SU2_CFD/src/solvers/CNEMONSSolver.cpp @@ -146,7 +146,7 @@ unsigned long CNEMONSSolver::SetPrimitive_Variables(CSolver **solver_container,C void CNEMONSSolver::SetPrimitive_Gradient_GG(CGeometry *geometry, const CConfig *config, bool reconstruction) { auto& gradient = reconstruction ? nodes->GetGradient_Reconstruction() : nodes->GetGradient_Primitive(); - const auto comm = reconstruction? PRIMITIVE_GRAD_REC : PRIMITIVE_GRADIENT; + const auto comm = reconstruction? MPI_QUANTITIES::PRIMITIVE_GRAD_REC : MPI_QUANTITIES::PRIMITIVE_GRADIENT; const auto commPer = reconstruction? PERIODIC_PRIM_GG_R : PERIODIC_PRIM_GG; /*--- Get indices of species & mixture density ---*/ diff --git a/SU2_CFD/src/solvers/CRadP1Solver.cpp b/SU2_CFD/src/solvers/CRadP1Solver.cpp index af40b22f2bc..d102e3b1d9a 100644 --- a/SU2_CFD/src/solvers/CRadP1Solver.cpp +++ b/SU2_CFD/src/solvers/CRadP1Solver.cpp @@ -553,8 +553,8 @@ void CRadP1Solver::ImplicitEuler_Iteration(CGeometry *geometry, CSolver **solver /*--- MPI solution ---*/ - InitiateComms(geometry, config, SOLUTION); - CompleteComms(geometry, config, SOLUTION); + InitiateComms(geometry, config, MPI_QUANTITIES::SOLUTION); + CompleteComms(geometry, config, MPI_QUANTITIES::SOLUTION); /*--- Compute the root mean square residual ---*/ diff --git a/SU2_CFD/src/solvers/CRadSolver.cpp b/SU2_CFD/src/solvers/CRadSolver.cpp index ca5becad7af..735983fa586 100644 --- a/SU2_CFD/src/solvers/CRadSolver.cpp +++ b/SU2_CFD/src/solvers/CRadSolver.cpp @@ -149,8 +149,8 @@ void CRadSolver::LoadRestart(CGeometry **geometry, CSolver ***solver, CConfig *c } /*--- MPI communication ---*/ - solver[MESH_0][RAD_SOL]->InitiateComms(geometry[MESH_0], config, SOLUTION); - solver[MESH_0][RAD_SOL]->CompleteComms(geometry[MESH_0], config, SOLUTION); + solver[MESH_0][RAD_SOL]->InitiateComms(geometry[MESH_0], config, MPI_QUANTITIES::SOLUTION); + solver[MESH_0][RAD_SOL]->CompleteComms(geometry[MESH_0], config, MPI_QUANTITIES::SOLUTION); /*--- Preprocess the fluid solver to compute the primitive variables ---*/ solver[MESH_0][FLOW_SOL]->Preprocessing(geometry[MESH_0], solver[MESH_0], config, MESH_0, NO_RK_ITER, RUNTIME_FLOW_SYS, false); diff --git a/SU2_CFD/src/solvers/CSolver.cpp b/SU2_CFD/src/solvers/CSolver.cpp index 03c79b46fc3..f7af2746cfe 100644 --- a/SU2_CFD/src/solvers/CSolver.cpp +++ b/SU2_CFD/src/solvers/CSolver.cpp @@ -1322,60 +1322,60 @@ void CSolver::CompletePeriodicComms(CGeometry *geometry, } void CSolver::GetCommCountAndType(const CConfig* config, - unsigned short commType, + MPI_QUANTITIES commType, unsigned short &COUNT_PER_POINT, unsigned short &MPI_TYPE) const { switch (commType) { - case SOLUTION: - case SOLUTION_OLD: - case UNDIVIDED_LAPLACIAN: - case SOLUTION_LIMITER: + case MPI_QUANTITIES::SOLUTION: + case MPI_QUANTITIES::SOLUTION_OLD: + case MPI_QUANTITIES::UNDIVIDED_LAPLACIAN: + case MPI_QUANTITIES::SOLUTION_LIMITER: COUNT_PER_POINT = nVar; MPI_TYPE = COMM_TYPE_DOUBLE; break; - case MAX_EIGENVALUE: - case SENSOR: + case MPI_QUANTITIES::MAX_EIGENVALUE: + case MPI_QUANTITIES::SENSOR: COUNT_PER_POINT = 1; MPI_TYPE = COMM_TYPE_DOUBLE; break; - case SOLUTION_GRADIENT: - case SOLUTION_GRAD_REC: + case MPI_QUANTITIES::SOLUTION_GRADIENT: + case MPI_QUANTITIES::SOLUTION_GRAD_REC: COUNT_PER_POINT = nVar*nDim; MPI_TYPE = COMM_TYPE_DOUBLE; break; - case PRIMITIVE_GRADIENT: - case PRIMITIVE_GRAD_REC: + case MPI_QUANTITIES::PRIMITIVE_GRADIENT: + case MPI_QUANTITIES::PRIMITIVE_GRAD_REC: COUNT_PER_POINT = nPrimVarGrad*nDim; MPI_TYPE = COMM_TYPE_DOUBLE; break; - case PRIMITIVE_LIMITER: + case MPI_QUANTITIES::PRIMITIVE_LIMITER: COUNT_PER_POINT = nPrimVarGrad; MPI_TYPE = COMM_TYPE_DOUBLE; break; - case SOLUTION_EDDY: + case MPI_QUANTITIES::SOLUTION_EDDY: COUNT_PER_POINT = nVar+1; MPI_TYPE = COMM_TYPE_DOUBLE; break; - case SOLUTION_FEA: + case MPI_QUANTITIES::SOLUTION_FEA: if (config->GetTime_Domain()) COUNT_PER_POINT = nVar*3; else COUNT_PER_POINT = nVar; MPI_TYPE = COMM_TYPE_DOUBLE; break; - case AUXVAR_GRADIENT: + case MPI_QUANTITIES::AUXVAR_GRADIENT: COUNT_PER_POINT = nDim*base_nodes->GetnAuxVar(); MPI_TYPE = COMM_TYPE_DOUBLE; break; - case MESH_DISPLACEMENTS: + case MPI_QUANTITIES::MESH_DISPLACEMENTS: COUNT_PER_POINT = nDim; MPI_TYPE = COMM_TYPE_DOUBLE; break; - case SOLUTION_TIME_N: + case MPI_QUANTITIES::SOLUTION_TIME_N: COUNT_PER_POINT = nVar; MPI_TYPE = COMM_TYPE_DOUBLE; break; - case SOLUTION_TIME_N1: + case MPI_QUANTITIES::SOLUTION_TIME_N1: COUNT_PER_POINT = nVar; MPI_TYPE = COMM_TYPE_DOUBLE; break; @@ -1387,25 +1387,25 @@ void CSolver::GetCommCountAndType(const CConfig* config, } namespace CommHelpers { - CVectorOfMatrix& selectGradient(CVariable* nodes, unsigned short commType) { + CVectorOfMatrix& selectGradient(CVariable* nodes, MPI_QUANTITIES commType) { switch(commType) { - case SOLUTION_GRAD_REC: return nodes->GetGradient_Reconstruction(); - case PRIMITIVE_GRADIENT: return nodes->GetGradient_Primitive(); - case PRIMITIVE_GRAD_REC: return nodes->GetGradient_Reconstruction(); - case AUXVAR_GRADIENT: return nodes->GetAuxVarGradient(); + case MPI_QUANTITIES::SOLUTION_GRAD_REC: return nodes->GetGradient_Reconstruction(); + case MPI_QUANTITIES::PRIMITIVE_GRADIENT: return nodes->GetGradient_Primitive(); + case MPI_QUANTITIES::PRIMITIVE_GRAD_REC: return nodes->GetGradient_Reconstruction(); + case MPI_QUANTITIES::AUXVAR_GRADIENT: return nodes->GetAuxVarGradient(); default: return nodes->GetGradient(); } } - su2activematrix& selectLimiter(CVariable* nodes, unsigned short commType) { - if (commType == PRIMITIVE_LIMITER) return nodes->GetLimiter_Primitive(); + su2activematrix& selectLimiter(CVariable* nodes, MPI_QUANTITIES commType) { + if (commType == MPI_QUANTITIES::PRIMITIVE_LIMITER) return nodes->GetLimiter_Primitive(); return nodes->GetLimiter(); } } void CSolver::InitiateComms(CGeometry *geometry, const CConfig *config, - unsigned short commType) { + MPI_QUANTITIES commType) { /*--- Local variables ---*/ @@ -1470,44 +1470,44 @@ void CSolver::InitiateComms(CGeometry *geometry, buf_offset = (msg_offset + iSend)*COUNT_PER_POINT; switch (commType) { - case SOLUTION: + case MPI_QUANTITIES::SOLUTION: for (iVar = 0; iVar < nVar; iVar++) bufDSend[buf_offset+iVar] = base_nodes->GetSolution(iPoint, iVar); break; - case SOLUTION_OLD: + case MPI_QUANTITIES::SOLUTION_OLD: for (iVar = 0; iVar < nVar; iVar++) bufDSend[buf_offset+iVar] = base_nodes->GetSolution_Old(iPoint, iVar); break; - case SOLUTION_EDDY: + case MPI_QUANTITIES::SOLUTION_EDDY: for (iVar = 0; iVar < nVar; iVar++) bufDSend[buf_offset+iVar] = base_nodes->GetSolution(iPoint, iVar); bufDSend[buf_offset+nVar] = base_nodes->GetmuT(iPoint); break; - case UNDIVIDED_LAPLACIAN: + case MPI_QUANTITIES::UNDIVIDED_LAPLACIAN: for (iVar = 0; iVar < nVar; iVar++) bufDSend[buf_offset+iVar] = base_nodes->GetUndivided_Laplacian(iPoint, iVar); break; - case SOLUTION_LIMITER: - case PRIMITIVE_LIMITER: + case MPI_QUANTITIES::SOLUTION_LIMITER: + case MPI_QUANTITIES::PRIMITIVE_LIMITER: for (iVar = 0; iVar < COUNT_PER_POINT; iVar++) bufDSend[buf_offset+iVar] = limiter(iPoint, iVar); break; - case MAX_EIGENVALUE: + case MPI_QUANTITIES::MAX_EIGENVALUE: bufDSend[buf_offset] = base_nodes->GetLambda(iPoint); break; - case SENSOR: + case MPI_QUANTITIES::SENSOR: bufDSend[buf_offset] = base_nodes->GetSensor(iPoint); break; - case SOLUTION_GRADIENT: - case PRIMITIVE_GRADIENT: - case SOLUTION_GRAD_REC: - case PRIMITIVE_GRAD_REC: - case AUXVAR_GRADIENT: + case MPI_QUANTITIES::SOLUTION_GRADIENT: + case MPI_QUANTITIES::PRIMITIVE_GRADIENT: + case MPI_QUANTITIES::SOLUTION_GRAD_REC: + case MPI_QUANTITIES::PRIMITIVE_GRAD_REC: + case MPI_QUANTITIES::AUXVAR_GRADIENT: for (iVar = 0; iVar < nVarGrad; iVar++) for (iDim = 0; iDim < nDim; iDim++) bufDSend[buf_offset+iVar*nDim+iDim] = gradient(iPoint, iVar, iDim); break; - case SOLUTION_FEA: + case MPI_QUANTITIES::SOLUTION_FEA: for (iVar = 0; iVar < nVar; iVar++) { bufDSend[buf_offset+iVar] = base_nodes->GetSolution(iPoint, iVar); if (config->GetTime_Domain()) { @@ -1516,15 +1516,15 @@ void CSolver::InitiateComms(CGeometry *geometry, } } break; - case MESH_DISPLACEMENTS: + case MPI_QUANTITIES::MESH_DISPLACEMENTS: for (iDim = 0; iDim < nDim; iDim++) bufDSend[buf_offset+iDim] = base_nodes->GetBound_Disp(iPoint, iDim); break; - case SOLUTION_TIME_N: + case MPI_QUANTITIES::SOLUTION_TIME_N: for (iVar = 0; iVar < nVar; iVar++) bufDSend[buf_offset+iVar] = base_nodes->GetSolution_time_n(iPoint, iVar); break; - case SOLUTION_TIME_N1: + case MPI_QUANTITIES::SOLUTION_TIME_N1: for (iVar = 0; iVar < nVar; iVar++) bufDSend[buf_offset+iVar] = base_nodes->GetSolution_time_n1(iPoint, iVar); break; @@ -1547,7 +1547,7 @@ void CSolver::InitiateComms(CGeometry *geometry, void CSolver::CompleteComms(CGeometry *geometry, const CConfig *config, - unsigned short commType) { + MPI_QUANTITIES commType) { /*--- Local variables ---*/ @@ -1618,44 +1618,44 @@ void CSolver::CompleteComms(CGeometry *geometry, /*--- Store the data correctly depending on the quantity. ---*/ switch (commType) { - case SOLUTION: + case MPI_QUANTITIES::SOLUTION: for (iVar = 0; iVar < nVar; iVar++) base_nodes->SetSolution(iPoint, iVar, bufDRecv[buf_offset+iVar]); break; - case SOLUTION_OLD: + case MPI_QUANTITIES::SOLUTION_OLD: for (iVar = 0; iVar < nVar; iVar++) base_nodes->SetSolution_Old(iPoint, iVar, bufDRecv[buf_offset+iVar]); break; - case SOLUTION_EDDY: + case MPI_QUANTITIES::SOLUTION_EDDY: for (iVar = 0; iVar < nVar; iVar++) base_nodes->SetSolution(iPoint, iVar, bufDRecv[buf_offset+iVar]); base_nodes->SetmuT(iPoint,bufDRecv[buf_offset+nVar]); break; - case UNDIVIDED_LAPLACIAN: + case MPI_QUANTITIES::UNDIVIDED_LAPLACIAN: for (iVar = 0; iVar < nVar; iVar++) base_nodes->SetUnd_Lapl(iPoint, iVar, bufDRecv[buf_offset+iVar]); break; - case SOLUTION_LIMITER: - case PRIMITIVE_LIMITER: + case MPI_QUANTITIES::SOLUTION_LIMITER: + case MPI_QUANTITIES::PRIMITIVE_LIMITER: for (iVar = 0; iVar < COUNT_PER_POINT; iVar++) limiter(iPoint,iVar) = bufDRecv[buf_offset+iVar]; break; - case MAX_EIGENVALUE: + case MPI_QUANTITIES::MAX_EIGENVALUE: base_nodes->SetLambda(iPoint,bufDRecv[buf_offset]); break; - case SENSOR: + case MPI_QUANTITIES::SENSOR: base_nodes->SetSensor(iPoint,bufDRecv[buf_offset]); break; - case SOLUTION_GRADIENT: - case PRIMITIVE_GRADIENT: - case SOLUTION_GRAD_REC: - case PRIMITIVE_GRAD_REC: - case AUXVAR_GRADIENT: + case MPI_QUANTITIES::SOLUTION_GRADIENT: + case MPI_QUANTITIES::PRIMITIVE_GRADIENT: + case MPI_QUANTITIES::SOLUTION_GRAD_REC: + case MPI_QUANTITIES::PRIMITIVE_GRAD_REC: + case MPI_QUANTITIES::AUXVAR_GRADIENT: for (iVar = 0; iVar < nVarGrad; iVar++) for (iDim = 0; iDim < nDim; iDim++) gradient(iPoint,iVar,iDim) = bufDRecv[buf_offset+iVar*nDim+iDim]; break; - case SOLUTION_FEA: + case MPI_QUANTITIES::SOLUTION_FEA: for (iVar = 0; iVar < nVar; iVar++) { base_nodes->SetSolution(iPoint, iVar, bufDRecv[buf_offset+iVar]); if (config->GetTime_Domain()) { @@ -1664,15 +1664,15 @@ void CSolver::CompleteComms(CGeometry *geometry, } } break; - case MESH_DISPLACEMENTS: + case MPI_QUANTITIES::MESH_DISPLACEMENTS: for (iDim = 0; iDim < nDim; iDim++) base_nodes->SetBound_Disp(iPoint, iDim, bufDRecv[buf_offset+iDim]); break; - case SOLUTION_TIME_N: + case MPI_QUANTITIES::SOLUTION_TIME_N: for (iVar = 0; iVar < nVar; iVar++) base_nodes->Set_Solution_time_n(iPoint, iVar, bufDRecv[buf_offset+iVar]); break; - case SOLUTION_TIME_N1: + case MPI_QUANTITIES::SOLUTION_TIME_N1: for (iVar = 0; iVar < nVar; iVar++) base_nodes->Set_Solution_time_n1(iPoint, iVar, bufDRecv[buf_offset+iVar]); break; @@ -2095,7 +2095,7 @@ void CSolver::SetAuxVar_Gradient_GG(CGeometry *geometry, const CConfig *config) const auto& solution = base_nodes->GetAuxVar(); auto& gradient = base_nodes->GetAuxVarGradient(); - computeGradientsGreenGauss(this, AUXVAR_GRADIENT, PERIODIC_NONE, *geometry, + computeGradientsGreenGauss(this, MPI_QUANTITIES::AUXVAR_GRADIENT, PERIODIC_NONE, *geometry, *config, solution, 0, base_nodes->GetnAuxVar(), gradient); } @@ -2106,7 +2106,7 @@ void CSolver::SetAuxVar_Gradient_LS(CGeometry *geometry, const CConfig *config) auto& gradient = base_nodes->GetAuxVarGradient(); auto& rmatrix = base_nodes->GetRmatrix(); - computeGradientsLeastSquares(this, AUXVAR_GRADIENT, PERIODIC_NONE, *geometry, *config, + computeGradientsLeastSquares(this, MPI_QUANTITIES::AUXVAR_GRADIENT, PERIODIC_NONE, *geometry, *config, weighted, solution, 0, base_nodes->GetnAuxVar(), gradient, rmatrix); } @@ -2114,7 +2114,7 @@ void CSolver::SetSolution_Gradient_GG(CGeometry *geometry, const CConfig *config const auto& solution = base_nodes->GetSolution(); auto& gradient = reconstruction? base_nodes->GetGradient_Reconstruction() : base_nodes->GetGradient(); - const auto comm = reconstruction? SOLUTION_GRAD_REC : SOLUTION_GRADIENT; + const auto comm = reconstruction? MPI_QUANTITIES::SOLUTION_GRAD_REC : MPI_QUANTITIES::SOLUTION_GRADIENT; const auto commPer = reconstruction? PERIODIC_SOL_GG_R : PERIODIC_SOL_GG; computeGradientsGreenGauss(this, comm, commPer, *geometry, *config, solution, 0, nVar, gradient); @@ -2138,7 +2138,7 @@ void CSolver::SetSolution_Gradient_LS(CGeometry *geometry, const CConfig *config const auto& solution = base_nodes->GetSolution(); auto& rmatrix = base_nodes->GetRmatrix(); auto& gradient = reconstruction? base_nodes->GetGradient_Reconstruction() : base_nodes->GetGradient(); - const auto comm = reconstruction? SOLUTION_GRAD_REC : SOLUTION_GRADIENT; + const auto comm = reconstruction? MPI_QUANTITIES::SOLUTION_GRAD_REC : MPI_QUANTITIES::SOLUTION_GRADIENT; computeGradientsLeastSquares(this, comm, commPer, *geometry, *config, weighted, solution, 0, nVar, gradient, rmatrix); } @@ -2183,8 +2183,8 @@ void CSolver::SetUndivided_Laplacian(CGeometry *geometry, const CConfig *config) /*--- MPI parallelization ---*/ - InitiateComms(geometry, config, UNDIVIDED_LAPLACIAN); - CompleteComms(geometry, config, UNDIVIDED_LAPLACIAN); + InitiateComms(geometry, config, MPI_QUANTITIES::UNDIVIDED_LAPLACIAN); + CompleteComms(geometry, config, MPI_QUANTITIES::UNDIVIDED_LAPLACIAN); } @@ -2238,7 +2238,7 @@ void CSolver::SetGridVel_Gradient(CGeometry *geometry, const CConfig *config) co auto& gridVelGrad = geometry->nodes->GetGridVel_Grad(); auto rmatrix = CVectorOfMatrix(nPoint,nDim,nDim); - computeGradientsLeastSquares(nullptr, GRID_VELOCITY, PERIODIC_NONE, *geometry, *config, + computeGradientsLeastSquares(nullptr, MPI_QUANTITIES::GRID_VELOCITY, PERIODIC_NONE, *geometry, *config, true, gridVel, 0, nDim, gridVelGrad, rmatrix); } @@ -2251,7 +2251,7 @@ void CSolver::SetSolution_Limiter(CGeometry *geometry, const CConfig *config) { auto& solMax = base_nodes->GetSolution_Max(); auto& limiter = base_nodes->GetLimiter(); - computeLimiters(kindLimiter, this, SOLUTION_LIMITER, PERIODIC_LIM_SOL_1, PERIODIC_LIM_SOL_2, + computeLimiters(kindLimiter, this, MPI_QUANTITIES::SOLUTION_LIMITER, PERIODIC_LIM_SOL_1, PERIODIC_LIM_SOL_2, *geometry, *config, 0, nVar, solution, gradient, solMin, solMax, limiter); } @@ -2723,8 +2723,8 @@ void CSolver::Restart_OldGeometry(CGeometry *geometry, CConfig *config) const { /*--- It's necessary to communicate this information ---*/ - geometry->InitiateComms(geometry, config, COORDINATES_OLD); - geometry->CompleteComms(geometry, config, COORDINATES_OLD); + geometry->InitiateComms(geometry, config, MPI_QUANTITIES::COORDINATES_OLD); + geometry->CompleteComms(geometry, config, MPI_QUANTITIES::COORDINATES_OLD); } diff --git a/SU2_CFD/src/solvers/CSpeciesFlameletSolver.cpp b/SU2_CFD/src/solvers/CSpeciesFlameletSolver.cpp index 66b90d3e39b..2484084e46f 100644 --- a/SU2_CFD/src/solvers/CSpeciesFlameletSolver.cpp +++ b/SU2_CFD/src/solvers/CSpeciesFlameletSolver.cpp @@ -68,9 +68,19 @@ void CSpeciesFlameletSolver::Preprocessing(CGeometry* geometry, CSolver** solver unsigned short RunTime_EqSystem, bool Output) { unsigned long n_not_in_domain_local = 0, n_not_in_domain_global = 0; vector scalars_vector(nVar); - + unsigned long spark_iter_start, spark_duration; + bool ignition = false; auto* flowNodes = su2staticcast_p(solver_container[FLOW_SOL]->GetNodes()); + /*--- Retrieve spark ignition parameters for spark-type ignition. ---*/ + if ((config->GetFlameletInitType() == FLAMELET_INIT_TYPE::SPARK) && !config->GetRestart()) { + auto spark_init = config->GetFlameInit(); + spark_iter_start = ceil(spark_init[4]); + spark_duration = ceil(spark_init[5]); + unsigned long iter = config->GetMultizone_Problem() ? config->GetOuterIter() : config->GetInnerIter(); + ignition = ((iter >= spark_iter_start) && (iter <= (spark_iter_start + spark_duration)) && !config->GetRestart()); + } + SU2_OMP_SAFE_GLOBAL_ACCESS(config->SetGlobalParam(config->GetKind_Solver(), RunTime_EqSystem);) SU2_OMP_FOR_STAT(omp_chunk_size) @@ -81,19 +91,35 @@ void CSpeciesFlameletSolver::Preprocessing(CGeometry* geometry, CSolver** solver /*--- Compute total source terms from the production and consumption. ---*/ unsigned long misses = SetScalarSources(config, fluid_model_local, i_point, scalars_vector); + + if (ignition) { + /*--- Apply source terms within spark radius. ---*/ + su2double dist_from_center = 0, + spark_radius = config->GetFlameInit()[3]; + dist_from_center = GeometryToolbox::SquaredDistance(nDim, geometry->nodes->GetCoord(i_point), config->GetFlameInit()); + if (dist_from_center < pow(spark_radius,2)) { + for (auto iVar = 0u; iVar < nVar; iVar++) + nodes->SetScalarSource(i_point, iVar, nodes->GetScalarSources(i_point)[iVar] + config->GetSpark()[iVar]); + } + } + nodes->SetTableMisses(i_point, misses); n_not_in_domain_local += misses; /*--- Obtain passive look-up scalars. ---*/ SetScalarLookUps(config, fluid_model_local, i_point, scalars_vector); /*--- Set mass diffusivity based on thermodynamic state. ---*/ - su2double T = flowNodes->GetTemperature(i_point); + auto T = flowNodes->GetTemperature(i_point); fluid_model_local->SetTDState_T(T, scalars); /*--- set the diffusivity in the fluid model to the diffusivity obtained from the lookup table ---*/ for (auto i_scalar = 0u; i_scalar < nVar; ++i_scalar) { nodes->SetDiffusivity(i_point, fluid_model_local->GetMassDiffusivity(i_scalar), i_scalar); } + /*--- Obtain preferential diffusion scalar values. ---*/ + if (config->GetPreferentialDiffusion()) + SetPreferentialDiffusionScalars(config, fluid_model_local, i_point, scalars_vector); + if (!Output) LinSysRes.SetBlock_Zero(i_point); } END_SU2_OMP_FOR @@ -102,6 +128,20 @@ void CSpeciesFlameletSolver::Preprocessing(CGeometry* geometry, CSolver** solver SU2_MPI::GetComm()); if ((rank == MASTER_NODE) && (n_not_in_domain_global > 0)) cout << "Number of points outside manifold domain: " << n_not_in_domain_global << endl; + + /*--- Compute preferential diffusion scalar gradients. ---*/ + if (config->GetPreferentialDiffusion()) { + switch (config->GetKind_Gradient_Method()) { + case GREEN_GAUSS: + SetAuxVar_Gradient_GG(geometry, config); + break; + case WEIGHTED_LEAST_SQUARES: + SetAuxVar_Gradient_LS(geometry, config); + break; + default: + break; + } + } /*--- Clear Residual and Jacobian. Upwind second order reconstruction and gradients ---*/ CommonPreprocessing(geometry, config, Output); } @@ -115,13 +155,22 @@ void CSpeciesFlameletSolver::SetInitialCondition(CGeometry** geometry, CSolver** cout << "Initializing progress variable and total enthalpy (using temperature)" << endl; } - const su2double* flame_init = config->GetFlameInit(); - const su2double flame_offset[3] = {flame_init[0], flame_init[1], flame_init[2]}; - const su2double flame_normal[3] = {flame_init[3], flame_init[4], flame_init[5]}; - const su2double flame_thickness = flame_init[6]; - const su2double flame_burnt_thickness = flame_init[7]; + su2double flame_offset[3] = {0, 0, 0}, flame_normal[3] = {0, 0, 0}, flame_thickness = 0, flame_burnt_thickness = 0, + flamenorm = 0; + bool flame_front_ignition = (config->GetFlameletInitType() == FLAMELET_INIT_TYPE::FLAME_FRONT); + + if (flame_front_ignition) { + /*--- Collect flame front ignition parameters. ---*/ + auto flame_init = config->GetFlameInit(); + for (auto iDim = 0u; iDim < 3; ++iDim) { + flame_offset[iDim] = flame_init[iDim]; + flame_normal[iDim] = flame_init[3 + iDim]; + } + flame_thickness = flame_init[6]; + flame_burnt_thickness = flame_init[7]; + flamenorm = GeometryToolbox::Norm(nDim, flame_normal); + } - const su2double flamenorm = GeometryToolbox::Norm(nDim, flame_normal); const su2double temp_inlet = config->GetInc_Temperature_Init(); su2double enth_inlet = config->GetSpecies_Init()[I_ENTH]; @@ -134,6 +183,19 @@ void CSpeciesFlameletSolver::SetInitialCondition(CGeometry** geometry, CSolver** const auto& cv_name = config->GetControllingVariableName(iCV); cout << "initial condition: " << cv_name << " = " << config->GetSpecies_Init()[iCV] << endl; } + switch (config->GetFlameletInitType()) { + case FLAMELET_INIT_TYPE::FLAME_FRONT: + cout << "Ignition with a straight flame front" << endl; + break; + case FLAMELET_INIT_TYPE::SPARK: + cout << "Ignition with an artificial spark" << endl; + break; + case FLAMELET_INIT_TYPE::NONE: + cout << "No solution ignition (cold flow)" << endl; + break; + default: + break; + } } CFluidModel* fluid_model_local; @@ -157,37 +219,40 @@ void CSpeciesFlameletSolver::SetInitialCondition(CGeometry** geometry, CSolver** for (unsigned long i_point = 0; i_point < nPoint; i_point++) { auto coords = geometry[i_mesh]->nodes->GetCoord(i_point); - /*--- Determine if point is above or below the plane, assuming the normal - is pointing towards the burned region. ---*/ - point_loc = 0.0; - for (unsigned short i_dim = 0; i_dim < geometry[i_mesh]->GetnDim(); i_dim++) { - point_loc += flame_normal[i_dim] * (coords[i_dim] - flame_offset[i_dim]); - } - - /*--- Compute the exact distance from point to plane. ---*/ - point_loc = point_loc / flamenorm; - - /* --- Unburnt region upstream of the flame. --- */ - if (point_loc <= 0) { - scalar_init[I_PROGVAR] = prog_unburnt; - n_points_unburnt_local++; - - /* --- Flame zone where we lineary increase from unburnt to burnt conditions. --- */ - } else if ((point_loc > 0) && (point_loc <= flame_thickness)) { - scalar_init[I_PROGVAR] = prog_unburnt + point_loc * (prog_burnt - prog_unburnt) / flame_thickness; - n_points_flame_local++; - - /* --- Burnt region behind the flame zone. --- */ - } else if ((point_loc > flame_thickness) && (point_loc <= flame_thickness + flame_burnt_thickness)) { - scalar_init[I_PROGVAR] = prog_burnt; - n_points_burnt_local++; - - /* --- Unburnt region downstream of the flame and burnt region. --- */ + if (flame_front_ignition) { + /*--- Determine if point is above or below the plane, assuming the normal + is pointing towards the burned region. ---*/ + point_loc = 0.0; + for (unsigned short i_dim = 0; i_dim < geometry[i_mesh]->GetnDim(); i_dim++) { + point_loc += flame_normal[i_dim] * (coords[i_dim] - flame_offset[i_dim]); + } + + /*--- Compute the exact distance from point to plane. ---*/ + point_loc = point_loc / flamenorm; + + /* --- Unburnt region upstream of the flame. --- */ + if (point_loc <= 0) { + scalar_init[I_PROGVAR] = prog_unburnt; + n_points_unburnt_local++; + + /* --- Flame zone where we lineary increase from unburnt to burnt conditions. --- */ + } else if ((point_loc > 0) && (point_loc <= flame_thickness)) { + scalar_init[I_PROGVAR] = prog_unburnt + point_loc * (prog_burnt - prog_unburnt) / flame_thickness; + n_points_flame_local++; + + /* --- Burnt region behind the flame zone. --- */ + } else if ((point_loc > flame_thickness) && (point_loc <= flame_thickness + flame_burnt_thickness)) { + scalar_init[I_PROGVAR] = prog_burnt; + n_points_burnt_local++; + + /* --- Unburnt region downstream of the flame and burnt region. --- */ + } else { + scalar_init[I_PROGVAR] = prog_unburnt; + n_points_unburnt_local++; + } } else { scalar_init[I_PROGVAR] = prog_unburnt; - n_points_unburnt_local++; } - /* --- Perform manifold evaluation to check whether initial scalar is out-of-bounds. --- */ fluid_model_local->SetTDState_T(temp_inlet, scalar_init); n_not_in_domain_local += fluid_model_local->GetExtrapolation(); @@ -200,11 +265,11 @@ void CSpeciesFlameletSolver::SetInitialCondition(CGeometry** geometry, CSolver** solver_container[i_mesh][SPECIES_SOL]->GetNodes()->SetSolution(i_point, scalar_init); } - solver_container[i_mesh][SPECIES_SOL]->InitiateComms(geometry[i_mesh], config, SOLUTION); - solver_container[i_mesh][SPECIES_SOL]->CompleteComms(geometry[i_mesh], config, SOLUTION); + solver_container[i_mesh][SPECIES_SOL]->InitiateComms(geometry[i_mesh], config, MPI_QUANTITIES::SOLUTION); + solver_container[i_mesh][SPECIES_SOL]->CompleteComms(geometry[i_mesh], config, MPI_QUANTITIES::SOLUTION); - solver_container[i_mesh][FLOW_SOL]->InitiateComms(geometry[i_mesh], config, SOLUTION); - solver_container[i_mesh][FLOW_SOL]->CompleteComms(geometry[i_mesh], config, SOLUTION); + solver_container[i_mesh][FLOW_SOL]->InitiateComms(geometry[i_mesh], config, MPI_QUANTITIES::SOLUTION); + solver_container[i_mesh][FLOW_SOL]->CompleteComms(geometry[i_mesh], config, MPI_QUANTITIES::SOLUTION); solver_container[i_mesh][FLOW_SOL]->Preprocessing(geometry[i_mesh], solver_container[i_mesh], config, i_mesh, NO_RK_ITER, RUNTIME_FLOW_SYS, false); @@ -225,9 +290,11 @@ void CSpeciesFlameletSolver::SetInitialCondition(CGeometry** geometry, CSolver** if (rank == MASTER_NODE) { cout << endl; - cout << " Number of points in unburnt region: " << n_points_unburnt_global << "." << endl; - cout << " Number of points in burnt region : " << n_points_burnt_global << "." << endl; - cout << " Number of points in flame zone : " << n_points_flame_global << "." << endl; + if (flame_front_ignition) { + cout << " Number of points in unburnt region: " << n_points_unburnt_global << "." << endl; + cout << " Number of points in burnt region : " << n_points_burnt_global << "." << endl; + cout << " Number of points in flame zone : " << n_points_flame_global << "." << endl; + } if (n_not_in_domain_global > 0) cout << " Initial condition: Number of points outside of table domain: " << n_not_in_domain_global << " !!!" @@ -331,8 +398,8 @@ void CSpeciesFlameletSolver::BC_Inlet(CGeometry* geometry, CSolver** solver_cont SU2_OMP_FOR_STAT(OMP_MIN_SIZE) for (auto iVertex = 0u; iVertex < geometry->nVertex[val_marker]; iVertex++) { Inlet_SpeciesVars[val_marker][iVertex][I_ENTH] = enth_inlet; - END_SU2_OMP_FOR } + END_SU2_OMP_FOR /*--- Call the general inlet boundary condition implementation. ---*/ CSpeciesSolver::BC_Inlet(geometry, solver_container, conv_numerics, visc_numerics, config, val_marker); @@ -343,10 +410,8 @@ void CSpeciesFlameletSolver::BC_Isothermal_Wall_Generic(CGeometry* geometry, CSo CConfig* config, unsigned short val_marker, bool cht_mode) { const bool implicit = config->GetKind_TimeIntScheme_Flow() == EULER_IMPLICIT; const string Marker_Tag = config->GetMarker_All_TagBound(val_marker); - su2double temp_wall; CFluidModel* fluid_model_local = solver_container[FLOW_SOL]->GetFluidModel(); auto* flowNodes = su2staticcast_p(solver_container[FLOW_SOL]->GetNodes()); - su2double enth_wall; unsigned long n_not_iterated = 0; /*--- Loop over all the vertices on this boundary marker. ---*/ @@ -354,6 +419,7 @@ void CSpeciesFlameletSolver::BC_Isothermal_Wall_Generic(CGeometry* geometry, CSo for (unsigned long iVertex = 0; iVertex < geometry->nVertex[val_marker]; iVertex++) { unsigned long iPoint = geometry->vertex[val_marker][iVertex]->GetNode(); + su2double temp_wall, enth_wall; if (cht_mode) temp_wall = solver_container[FLOW_SOL]->GetConjugateHeatVariable(val_marker, iVertex, 0); else @@ -402,7 +468,7 @@ void CSpeciesFlameletSolver::BC_Isothermal_Wall_Generic(CGeometry* geometry, CSo su2double dist_ij = sqrt(dist_ij_2); /*--- Compute the normal gradient in temperature using Twall. ---*/ - + ///TODO: Account for preferential diffusion in computation of the heat flux su2double dTdn = -(flowNodes->GetTemperature(Point_Normal) - temp_wall) / dist_ij; /*--- Get thermal conductivity. ---*/ @@ -442,13 +508,13 @@ unsigned long CSpeciesFlameletSolver::SetScalarSources(const CConfig* config, CF vector table_sources(config->GetNControlVars() + 2 * config->GetNUserScalars()); unsigned long misses = fluid_model_local->EvaluateDataSet(scalars, FLAMELET_LOOKUP_OPS::SOURCES, table_sources); + table_sources[I_PROGVAR] = fmax(0, table_sources[I_PROGVAR]); nodes->SetTableMisses(iPoint, misses); /*--- The source term for progress variable is always positive, we clip from below to makes sure. --- */ vector source_scalar(config->GetNScalars()); for (auto iCV = 0u; iCV < config->GetNControlVars(); iCV++) source_scalar[iCV] = table_sources[iCV]; - source_scalar[I_PROGVAR] = fmax(EPS, source_scalar[I_PROGVAR]); /*--- Source term for the auxiliary species transport equations. ---*/ for (size_t i_aux = 0; i_aux < config->GetNUserScalars(); i_aux++) { @@ -466,17 +532,212 @@ unsigned long CSpeciesFlameletSolver::SetScalarSources(const CConfig* config, CF unsigned long CSpeciesFlameletSolver::SetScalarLookUps(const CConfig* config, CFluidModel* fluid_model_local, unsigned long iPoint, const vector& scalars) { - /*--- Compute total source terms from the production and consumption. ---*/ + /*--- Retrieve the passive look-up variables from the manifold. ---*/ + unsigned long misses{0}; + /*--- Skip if no passive look-ups are listed ---*/ + if (config->GetNLookups() > 0) { + vector lookup_scalar(config->GetNLookups()); + misses = fluid_model_local->EvaluateDataSet(scalars, FLAMELET_LOOKUP_OPS::LOOKUP, lookup_scalar); + + for (auto i_lookup = 0u; i_lookup < config->GetNLookups(); i_lookup++) { + nodes->SetLookupScalar(iPoint, lookup_scalar[i_lookup], i_lookup); + } + } + + return misses; +} - vector lookup_scalar(config->GetNLookups()); - unsigned long misses = fluid_model_local->EvaluateDataSet(scalars, FLAMELET_LOOKUP_OPS::LOOKUP, lookup_scalar); +unsigned long CSpeciesFlameletSolver::SetPreferentialDiffusionScalars(const CConfig* config, + CFluidModel* fluid_model_local, + unsigned long iPoint, + const vector& scalars) { + /*--- Retrieve the preferential diffusion scalar values from the manifold. ---*/ - for (auto i_lookup = 0u; i_lookup < config->GetNLookups(); i_lookup++) { - nodes->SetLookupScalar(iPoint, lookup_scalar[i_lookup], i_lookup); + vector beta_scalar(FLAMELET_PREF_DIFF_SCALARS::N_BETA_TERMS); + unsigned long misses = fluid_model_local->EvaluateDataSet(scalars, FLAMELET_LOOKUP_OPS::PREFDIF, beta_scalar); + + for (auto i_beta = 0u; i_beta < FLAMELET_PREF_DIFF_SCALARS::N_BETA_TERMS; i_beta++) { + nodes->SetAuxVar(iPoint, i_beta, beta_scalar[i_beta]); } return misses; } +void CSpeciesFlameletSolver::Viscous_Residual(const unsigned long iEdge, const CGeometry* geometry, CSolver** solver_container, + CNumerics* numerics, const CConfig* config) { + /*--- Overloaded viscous residual method which accounts for preferential diffusion. ---*/ + const bool implicit = (config->GetKind_TimeIntScheme() == EULER_IMPLICIT), + PreferentialDiffusion = config->GetPreferentialDiffusion(); + + /*--- Points in edge ---*/ + auto iPoint = geometry->edges->GetNode(iEdge, 0); + auto jPoint = geometry->edges->GetNode(iEdge, 1); + + auto SolverSpecificNumerics = [&](unsigned long iPoint, unsigned long jPoint) { + /*--- Mass diffusivity coefficients. ---*/ + + numerics->SetDiffusionCoeff(nodes->GetDiffusivity(iPoint), nodes->GetDiffusivity(jPoint)); + }; + + /*--- Regular viscous scalar residual computation. ---*/ + + Viscous_Residual_impl(SolverSpecificNumerics, iEdge, geometry, solver_container, numerics, config); + + /*--- Viscous residual due to preferential diffusion ---*/ + if (PreferentialDiffusion) { + CFlowVariable* flowNodes = su2staticcast_p(solver_container[FLOW_SOL]->GetNodes()); + + su2double scalar_i[MAXNVAR] = {0}, + scalar_j[MAXNVAR] = {0}, + diff_coeff_beta_i[MAXNVAR] = {0}, + diff_coeff_beta_j[MAXNVAR] = {0}; + + // Number of active transport scalars + const auto n_CV = config->GetNControlVars(); + + su2activematrix scalar_grad_i(MAXNVAR, MAXNDIM), scalar_grad_j(MAXNVAR, MAXNDIM); + /*--- Looping over spatial dimensions to fill in the diffusion scalar gradients. ---*/ + /*--- The scalar gradient is subtracted to account for regular viscous diffusion. ---*/ + for (auto iScalar = 0u; iScalar < n_CV; ++iScalar) { + for (auto iDim = 0u; iDim < nDim; ++iDim) { + switch (iScalar) { + case I_PROGVAR: + scalar_grad_i[iScalar][iDim] = + nodes->GetAuxVarGradient(iPoint, FLAMELET_PREF_DIFF_SCALARS::I_BETA_PROGVAR, iDim) - + nodes->GetGradient(iPoint, I_PROGVAR, iDim); + scalar_grad_j[iScalar][iDim] = + nodes->GetAuxVarGradient(jPoint, FLAMELET_PREF_DIFF_SCALARS::I_BETA_PROGVAR, iDim) - + nodes->GetGradient(jPoint, I_PROGVAR, iDim); + break; + case I_ENTH: + scalar_grad_i[iScalar][iDim] = + nodes->GetAuxVarGradient(iPoint, FLAMELET_PREF_DIFF_SCALARS::I_BETA_ENTH, iDim) - + nodes->GetGradient(iPoint, I_ENTH, iDim); + scalar_grad_j[iScalar][iDim] = + nodes->GetAuxVarGradient(jPoint, FLAMELET_PREF_DIFF_SCALARS::I_BETA_ENTH, iDim) - + nodes->GetGradient(jPoint, I_ENTH, iDim); + break; + case I_MIXFRAC: + scalar_grad_i[iScalar][iDim] = + nodes->GetAuxVarGradient(iPoint, FLAMELET_PREF_DIFF_SCALARS::I_BETA_MIXFRAC, iDim) - + nodes->GetGradient(iPoint, I_MIXFRAC, iDim); + scalar_grad_j[iScalar][iDim] = + nodes->GetAuxVarGradient(jPoint, FLAMELET_PREF_DIFF_SCALARS::I_BETA_MIXFRAC, iDim) - + nodes->GetGradient(jPoint, I_MIXFRAC, iDim); + break; + default: + break; + } + } + } + /*--- No preferential diffusion modification for passive species. ---*/ + for (auto iScalar = n_CV; iScalar < nVar; ++iScalar) { + for (auto iDim = 0u; iDim < nDim; ++iDim) { + scalar_grad_i[iScalar][iDim] = 0; + scalar_grad_j[iScalar][iDim] = 0; + } + } + + for (auto iScalar = 0u; iScalar < n_CV; ++iScalar) { + /*--- Filling in the preferential diffusion scalars (beta_pv, beta_h2, beta_Z). ---*/ + switch (iScalar) { + case I_PROGVAR: + scalar_i[iScalar] = nodes->GetAuxVar(iPoint, FLAMELET_PREF_DIFF_SCALARS::I_BETA_PROGVAR) - + nodes->GetSolution(iPoint, iScalar); + scalar_j[iScalar] = nodes->GetAuxVar(jPoint, FLAMELET_PREF_DIFF_SCALARS::I_BETA_PROGVAR) - + nodes->GetSolution(jPoint, iScalar); + break; + case I_ENTH: + scalar_i[iScalar] = + nodes->GetAuxVar(iPoint, FLAMELET_PREF_DIFF_SCALARS::I_BETA_ENTH) - nodes->GetSolution(iPoint, iScalar); + scalar_j[iScalar] = + nodes->GetAuxVar(jPoint, FLAMELET_PREF_DIFF_SCALARS::I_BETA_ENTH) - nodes->GetSolution(jPoint, iScalar); + break; + case I_MIXFRAC: + scalar_i[iScalar] = nodes->GetAuxVar(iPoint, FLAMELET_PREF_DIFF_SCALARS::I_BETA_MIXFRAC) - + nodes->GetSolution(iPoint, iScalar); + scalar_j[iScalar] = nodes->GetAuxVar(jPoint, FLAMELET_PREF_DIFF_SCALARS::I_BETA_MIXFRAC) - + nodes->GetSolution(jPoint, iScalar); + break; + default: + break; + } + diff_coeff_beta_i[iScalar] = nodes->GetDiffusivity(iPoint, iScalar); + diff_coeff_beta_j[iScalar] = nodes->GetDiffusivity(jPoint, iScalar); + } + + for (auto iScalar = n_CV; iScalar < nVar; ++iScalar) { + scalar_i[iScalar] = 0; + scalar_j[iScalar] = 0; + diff_coeff_beta_i[iScalar] = 0; + diff_coeff_beta_j[iScalar] = 0; + } + + numerics->SetScalarVar(scalar_i, scalar_j); + + numerics->SetScalarVarGradient(CMatrixView(scalar_grad_i), CMatrixView(scalar_grad_j)); + + numerics->SetDiffusionCoeff(diff_coeff_beta_i, diff_coeff_beta_j); + + /*--- Computing first preferential residual component. ---*/ + auto residual_PD = numerics->ComputeResidual(config); + + if (ReducerStrategy) { + EdgeFluxes.SubtractBlock(iEdge, residual_PD); + + if (implicit) Jacobian.UpdateBlocksSub(iEdge, residual_PD.jacobian_i, residual_PD.jacobian_j); + } else { + LinSysRes.SubtractBlock(iPoint, residual_PD); + LinSysRes.AddBlock(jPoint, residual_PD); + /*--- Set implicit computation ---*/ + if (implicit) Jacobian.UpdateBlocksSub(iEdge, iPoint, jPoint, residual_PD.jacobian_i, residual_PD.jacobian_j); + } + + /* Computing the second preferential diffusion terms due to heat flux */ + for (auto iScalar = 0u; iScalar < nVar; ++iScalar) { + for (auto iDim = 0u; iDim < nDim; ++iDim) { + if (iScalar == I_ENTH) { + /* Setting the temperature gradient */ + scalar_grad_i[iScalar][iDim] = flowNodes->GetGradient_Primitive(iPoint, prim_idx.Temperature(), iDim); + scalar_grad_j[iScalar][iDim] = flowNodes->GetGradient_Primitive(jPoint, prim_idx.Temperature(), iDim); + } else { + scalar_grad_i[iScalar][iDim] = 0; + scalar_grad_j[iScalar][iDim] = 0; + } + } + + if (iScalar == I_ENTH) { + scalar_i[iScalar] = flowNodes->GetTemperature(iPoint); + scalar_j[iScalar] = flowNodes->GetTemperature(jPoint); + diff_coeff_beta_i[iScalar] = nodes->GetAuxVar(iPoint, FLAMELET_PREF_DIFF_SCALARS::I_BETA_ENTH_THERMAL) * + nodes->GetDiffusivity(iPoint, iScalar); + diff_coeff_beta_j[iScalar] = nodes->GetAuxVar(jPoint, FLAMELET_PREF_DIFF_SCALARS::I_BETA_ENTH_THERMAL) * + nodes->GetDiffusivity(jPoint, iScalar); + } else { + scalar_i[iScalar] = 0; + scalar_j[iScalar] = 0; + diff_coeff_beta_i[iScalar] = 0; + diff_coeff_beta_j[iScalar] = 0; + } + } + + numerics->SetScalarVar(scalar_i, scalar_j); + + numerics->SetScalarVarGradient(CMatrixView(scalar_grad_i), CMatrixView(scalar_grad_j)); + + numerics->SetDiffusionCoeff(diff_coeff_beta_i, diff_coeff_beta_j); + + auto residual_thermal = numerics->ComputeResidual(config); + + if (ReducerStrategy) { + EdgeFluxes.SubtractBlock(iEdge, residual_thermal); + } else { + LinSysRes.SubtractBlock(iPoint, residual_thermal); + LinSysRes.AddBlock(jPoint, residual_thermal); + /* No implicit part for the preferential diffusion of heat */ + } + } +} + unsigned long CSpeciesFlameletSolver::GetEnthFromTemp(CFluidModel* fluid_model, su2double const val_temp, const su2double* scalar_solution, su2double* val_enth) { /*--- convergence criterion for temperature in [K], high accuracy needed for restarts. ---*/ diff --git a/SU2_CFD/src/solvers/CSpeciesSolver.cpp b/SU2_CFD/src/solvers/CSpeciesSolver.cpp index f5f5ed16603..3c7ac895708 100644 --- a/SU2_CFD/src/solvers/CSpeciesSolver.cpp +++ b/SU2_CFD/src/solvers/CSpeciesSolver.cpp @@ -60,8 +60,8 @@ CSpeciesSolver::CSpeciesSolver(CGeometry* geometry, CConfig* config, unsigned sh /*--- MPI solution ---*/ - InitiateComms(geometry, config, SOLUTION); - CompleteComms(geometry, config, SOLUTION); + InitiateComms(geometry, config, MPI_QUANTITIES::SOLUTION); + CompleteComms(geometry, config, MPI_QUANTITIES::SOLUTION); SlidingState.resize(nMarker); SlidingStateNodes.resize(nMarker); @@ -248,8 +248,8 @@ void CSpeciesSolver::LoadRestart(CGeometry** geometry, CSolver*** solver, CConfi /*--- MPI solution and compute the eddy viscosity ---*/ - solver[MESH_0][SPECIES_SOL]->InitiateComms(geometry[MESH_0], config, SOLUTION); - solver[MESH_0][SPECIES_SOL]->CompleteComms(geometry[MESH_0], config, SOLUTION); + solver[MESH_0][SPECIES_SOL]->InitiateComms(geometry[MESH_0], config, MPI_QUANTITIES::SOLUTION); + solver[MESH_0][SPECIES_SOL]->CompleteComms(geometry[MESH_0], config, MPI_QUANTITIES::SOLUTION); // Flow-Pre computes/sets mixture properties solver[MESH_0][FLOW_SOL]->Preprocessing(geometry[MESH_0], solver[MESH_0], config, MESH_0, NO_RK_ITER, @@ -267,8 +267,8 @@ void CSpeciesSolver::LoadRestart(CGeometry** geometry, CSolver*** solver, CConfi for (auto iMesh = 1u; iMesh <= config->GetnMGLevels(); iMesh++) { MultigridRestriction(*geometry[iMesh - 1], solver[iMesh - 1][SPECIES_SOL]->GetNodes()->GetSolution(), *geometry[iMesh], solver[iMesh][SPECIES_SOL]->GetNodes()->GetSolution()); - solver[iMesh][SPECIES_SOL]->InitiateComms(geometry[iMesh], config, SOLUTION); - solver[iMesh][SPECIES_SOL]->CompleteComms(geometry[iMesh], config, SOLUTION); + solver[iMesh][SPECIES_SOL]->InitiateComms(geometry[iMesh], config, MPI_QUANTITIES::SOLUTION); + solver[iMesh][SPECIES_SOL]->CompleteComms(geometry[iMesh], config, MPI_QUANTITIES::SOLUTION); solver[iMesh][FLOW_SOL]->Preprocessing(geometry[iMesh], solver[iMesh], config, iMesh, NO_RK_ITER, RUNTIME_FLOW_SYS, false); @@ -314,8 +314,8 @@ void CSpeciesSolver::Preprocessing(CGeometry* geometry, CSolver** solver_contain CommonPreprocessing(geometry, config, Output); } -void CSpeciesSolver::Viscous_Residual(unsigned long iEdge, CGeometry* geometry, CSolver** solver_container, - CNumerics* numerics, CConfig* config) { +void CSpeciesSolver::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) { /*--- Mass diffusivity coefficients. ---*/ @@ -345,7 +345,7 @@ void CSpeciesSolver::BC_Inlet(CGeometry* geometry, CSolver** solver_container, C /*--- Identify the boundary by string name ---*/ string Marker_Tag = config->GetMarker_All_TagBound(val_marker); - + if (config->GetMarker_StrongBC(Marker_Tag)==true) { nodes->SetSolution_Old(iPoint, Inlet_SpeciesVars[val_marker][iVertex]); diff --git a/SU2_CFD/src/solvers/CTransLMSolver.cpp b/SU2_CFD/src/solvers/CTransLMSolver.cpp index b67e1b1af96..1fdf2a47de1 100644 --- a/SU2_CFD/src/solvers/CTransLMSolver.cpp +++ b/SU2_CFD/src/solvers/CTransLMSolver.cpp @@ -132,8 +132,8 @@ CTransLMSolver::CTransLMSolver(CGeometry *geometry, CConfig *config, unsigned sh /*--- MPI solution ---*/ - InitiateComms(geometry, config, SOLUTION); - CompleteComms(geometry, config, SOLUTION); + InitiateComms(geometry, config, MPI_QUANTITIES::SOLUTION); + CompleteComms(geometry, config, MPI_QUANTITIES::SOLUTION); /*--- Initializate quantities for SlidingMesh Interface ---*/ @@ -263,8 +263,8 @@ void CTransLMSolver::Postprocessing(CGeometry *geometry, CSolver **solver_contai } -void CTransLMSolver::Viscous_Residual(unsigned long iEdge, CGeometry* geometry, CSolver** solver_container, - CNumerics* numerics, CConfig* config) { +void CTransLMSolver::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. ---*/ @@ -561,8 +561,8 @@ void CTransLMSolver::LoadRestart(CGeometry** geometry, CSolver*** solver, CConfi /*--- MPI solution and compute the eddy viscosity ---*/ - solver[MESH_0][TRANS_SOL]->InitiateComms(geometry[MESH_0], config, SOLUTION); - solver[MESH_0][TRANS_SOL]->CompleteComms(geometry[MESH_0], config, SOLUTION); + solver[MESH_0][TRANS_SOL]->InitiateComms(geometry[MESH_0], config, MPI_QUANTITIES::SOLUTION); + solver[MESH_0][TRANS_SOL]->CompleteComms(geometry[MESH_0], config, MPI_QUANTITIES::SOLUTION); /*--- For turbulent+species simulations the solver Pre-/Postprocessing is done by the species solver. ---*/ if (config->GetKind_Species_Model() == SPECIES_MODEL::NONE) { @@ -578,8 +578,8 @@ void CTransLMSolver::LoadRestart(CGeometry** geometry, CSolver*** solver, CConfi MultigridRestriction(*geometry[iMesh - 1], solver[iMesh - 1][TRANS_SOL]->GetNodes()->GetSolution(), *geometry[iMesh], solver[iMesh][TRANS_SOL]->GetNodes()->GetSolution()); - solver[iMesh][TRANS_SOL]->InitiateComms(geometry[iMesh], config, SOLUTION); - solver[iMesh][TRANS_SOL]->CompleteComms(geometry[iMesh], config, SOLUTION); + solver[iMesh][TRANS_SOL]->InitiateComms(geometry[iMesh], config, MPI_QUANTITIES::SOLUTION); + solver[iMesh][TRANS_SOL]->CompleteComms(geometry[iMesh], config, MPI_QUANTITIES::SOLUTION); if (config->GetKind_Species_Model() == SPECIES_MODEL::NONE) { solver[iMesh][FLOW_SOL]->Preprocessing(geometry[iMesh], solver[iMesh], config, iMesh, NO_RK_ITER, RUNTIME_FLOW_SYS, diff --git a/SU2_CFD/src/solvers/CTurbSASolver.cpp b/SU2_CFD/src/solvers/CTurbSASolver.cpp index f3a6d57b589..f53596e9901 100644 --- a/SU2_CFD/src/solvers/CTurbSASolver.cpp +++ b/SU2_CFD/src/solvers/CTurbSASolver.cpp @@ -134,8 +134,8 @@ CTurbSASolver::CTurbSASolver(CGeometry *geometry, CConfig *config, unsigned shor /*--- MPI solution ---*/ - InitiateComms(geometry, config, SOLUTION_EDDY); - CompleteComms(geometry, config, SOLUTION_EDDY); + InitiateComms(geometry, config, MPI_QUANTITIES::SOLUTION_EDDY); + CompleteComms(geometry, config, MPI_QUANTITIES::SOLUTION_EDDY); /*--- Initializate quantities for SlidingMesh Interface ---*/ @@ -291,8 +291,8 @@ void CTurbSASolver::Postprocessing(CGeometry *geometry, CSolver **solver_contain AD::EndNoSharedReading(); } -void CTurbSASolver::Viscous_Residual(unsigned long iEdge, CGeometry* geometry, CSolver** solver_container, - CNumerics* numerics, CConfig* config) { +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) { diff --git a/SU2_CFD/src/solvers/CTurbSSTSolver.cpp b/SU2_CFD/src/solvers/CTurbSSTSolver.cpp index 1973b70a331..555f8e77155 100644 --- a/SU2_CFD/src/solvers/CTurbSSTSolver.cpp +++ b/SU2_CFD/src/solvers/CTurbSSTSolver.cpp @@ -142,8 +142,8 @@ CTurbSSTSolver::CTurbSSTSolver(CGeometry *geometry, CConfig *config, unsigned sh /*--- MPI solution ---*/ - InitiateComms(geometry, config, SOLUTION_EDDY); - CompleteComms(geometry, config, SOLUTION_EDDY); + InitiateComms(geometry, config, MPI_QUANTITIES::SOLUTION_EDDY); + CompleteComms(geometry, config, MPI_QUANTITIES::SOLUTION_EDDY); /*--- Initialize quantities for SlidingMesh Interface ---*/ @@ -276,8 +276,8 @@ void CTurbSSTSolver::Postprocessing(CGeometry *geometry, CSolver **solver_contai AD::EndNoSharedReading(); } -void CTurbSSTSolver::Viscous_Residual(unsigned long iEdge, CGeometry* geometry, CSolver** solver_container, - CNumerics* numerics, CConfig* config) { +void CTurbSSTSolver::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) { diff --git a/SU2_CFD/src/solvers/CTurbSolver.cpp b/SU2_CFD/src/solvers/CTurbSolver.cpp index c4b5a0b2e2b..d0744a7e223 100644 --- a/SU2_CFD/src/solvers/CTurbSolver.cpp +++ b/SU2_CFD/src/solvers/CTurbSolver.cpp @@ -167,8 +167,8 @@ void CTurbSolver::LoadRestart(CGeometry** geometry, CSolver*** solver, CConfig* /*--- MPI solution and compute the eddy viscosity ---*/ - solver[MESH_0][TURB_SOL]->InitiateComms(geometry[MESH_0], config, SOLUTION); - solver[MESH_0][TURB_SOL]->CompleteComms(geometry[MESH_0], config, SOLUTION); + solver[MESH_0][TURB_SOL]->InitiateComms(geometry[MESH_0], config, MPI_QUANTITIES::SOLUTION); + solver[MESH_0][TURB_SOL]->CompleteComms(geometry[MESH_0], config, MPI_QUANTITIES::SOLUTION); /*--- For turbulent+species simulations the solver Pre-/Postprocessing is done by the species solver. ---*/ if (config->GetKind_Species_Model() == SPECIES_MODEL::NONE && config->GetKind_Trans_Model() == TURB_TRANS_MODEL::NONE) { @@ -182,8 +182,8 @@ void CTurbSolver::LoadRestart(CGeometry** geometry, CSolver*** solver, CConfig* for (auto iMesh = 1u; iMesh <= config->GetnMGLevels(); iMesh++) { MultigridRestriction(*geometry[iMesh - 1], solver[iMesh - 1][TURB_SOL]->GetNodes()->GetSolution(), *geometry[iMesh], solver[iMesh][TURB_SOL]->GetNodes()->GetSolution()); - solver[iMesh][TURB_SOL]->InitiateComms(geometry[iMesh], config, SOLUTION); - solver[iMesh][TURB_SOL]->CompleteComms(geometry[iMesh], config, SOLUTION); + solver[iMesh][TURB_SOL]->InitiateComms(geometry[iMesh], config, MPI_QUANTITIES::SOLUTION); + solver[iMesh][TURB_SOL]->CompleteComms(geometry[iMesh], config, MPI_QUANTITIES::SOLUTION); if (config->GetKind_Species_Model() == SPECIES_MODEL::NONE) { solver[iMesh][FLOW_SOL]->Preprocessing(geometry[iMesh], solver[iMesh], config, iMesh, NO_RK_ITER, RUNTIME_FLOW_SYS, diff --git a/SU2_CFD/src/variables/CSpeciesFlameletVariable.cpp b/SU2_CFD/src/variables/CSpeciesFlameletVariable.cpp index da10ba21bd8..3dd1b222750 100644 --- a/SU2_CFD/src/variables/CSpeciesFlameletVariable.cpp +++ b/SU2_CFD/src/variables/CSpeciesFlameletVariable.cpp @@ -53,4 +53,9 @@ CSpeciesFlameletVariable::CSpeciesFlameletVariable(const su2double* species_inf, source_scalar.resize(nPoint, config->GetNScalars()) = su2double(0.0); lookup_scalar.resize(nPoint, config->GetNLookups()) = su2double(0.0); table_misses.resize(nPoint) = 0; + + if (config->GetPreferentialDiffusion()) { + AuxVar.resize(nPoint, FLAMELET_PREF_DIFF_SCALARS::N_BETA_TERMS) = su2double(0.0); + Grad_AuxVar.resize(nPoint, FLAMELET_PREF_DIFF_SCALARS::N_BETA_TERMS, nDim, 0.0); + } } diff --git a/SU2_GEO/src/SU2_GEO.cpp b/SU2_GEO/src/SU2_GEO.cpp index 86b9ee5c6e0..4e1d3bad9fe 100644 --- a/SU2_GEO/src/SU2_GEO.cpp +++ b/SU2_GEO/src/SU2_GEO.cpp @@ -1429,14 +1429,15 @@ int main(int argc, char* argv[]) { delete[] Xcoord_Airfoil; delete[] Ycoord_Airfoil; delete[] Zcoord_Airfoil; + delete[] Variable_Airfoil; delete[] ObjectiveFunc; delete[] ObjectiveFunc_New; delete[] Gradient; for (iPlane = 0; iPlane < nPlane; iPlane++) { - delete Plane_P0[iPlane]; - delete Plane_Normal[iPlane]; + delete[] Plane_P0[iPlane]; + delete[] Plane_Normal[iPlane]; } delete[] Plane_P0; delete[] Plane_Normal; diff --git a/TestCases/TestCase.py b/TestCases/TestCase.py index 6d577b47df4..0818c6074a0 100644 --- a/TestCases/TestCase.py +++ b/TestCases/TestCase.py @@ -45,6 +45,7 @@ def is_float(test_string): 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.') return parser.parse_args() class TestCase: @@ -114,6 +115,7 @@ def __init__(self,tag_in): self.cpu_arch = platform.machine().casefold() self.enabled_on_cpu_arch = ["x86_64","amd64","aarch64","arm64"] self.enabled_with_tsan = True + self.enabled_with_asan = True self.command = self.Command() self.timeout = 0 self.tol = 0.0 @@ -125,9 +127,9 @@ def __init__(self,tag_in): self.reference_file_aarch64 = "" self.test_file = "of_grad.dat" - def run_test(self, running_with_tsan=False): + def run_test(self, with_tsan=False, with_asan=False): - if not self.is_enabled(running_with_tsan): + if not self.is_enabled(with_tsan, with_asan): return True print('==================== Start Test: %s ===================='%self.tag) @@ -142,7 +144,7 @@ def run_test(self, running_with_tsan=False): # Adjust the number of iterations in the config file if len(self.test_vals) != 0: - self.adjust_iter(running_with_tsan) + self.adjust_iter(with_tsan, with_asan) # Check for disabling the restart if self.no_restart: @@ -186,7 +188,7 @@ def run_test(self, running_with_tsan=False): delta_vals = [] sim_vals = [] - if not running_with_tsan: # tsan findings result in non-zero return code, no need to examine the output + if not with_tsan and not with_asan: # sanitizer findings result in non-zero return code, no need to examine the output # Examine the output f = open(logfilename,'r') output = f.readlines() @@ -260,7 +262,7 @@ def run_test(self, running_with_tsan=False): if not start_solver: print('ERROR: The code was not able to get to the "Begin solver" section.') - if not running_with_tsan and iter_missing: + if not with_tsan and not with_asan and iter_missing: print('ERROR: The iteration number %d could not be found.'%self.test_iter) print('CPU architecture=%s' % self.cpu_arch) @@ -281,9 +283,9 @@ def run_test(self, running_with_tsan=False): os.chdir(workdir) return passed - def run_filediff(self, running_with_tsan=False): + def run_filediff(self, with_tsan=False, with_asan=False): - if not self.is_enabled(running_with_tsan): + if not self.is_enabled(with_tsan, with_asan): return True print('==================== Start Test: %s ===================='%self.tag) @@ -291,7 +293,7 @@ def run_filediff(self, running_with_tsan=False): timed_out = False # Adjust the number of iterations in the config file - self.adjust_iter(running_with_tsan) + self.adjust_iter(with_tsan, with_asan) self.adjust_test_data() @@ -330,7 +332,7 @@ def run_filediff(self, running_with_tsan=False): print("Output from the failed case:") subprocess.call(["cat", logfilename]) - if not running_with_tsan: # thread sanitizer tests only check the return code, no need to compare outputs + if not with_tsan and not with_asan: # sanitizer tests only check the return code, no need to compare outputs diff_time_start = datetime.datetime.now() if not timed_out and passed: # Compare files @@ -468,9 +470,9 @@ def run_filediff(self, running_with_tsan=False): os.chdir(workdir) return passed - def run_opt(self): + def run_opt(self, with_tsan=False, with_asan=False): - if not self.is_enabled(): + if not self.is_enabled(with_tsan, with_asan): return True print('==================== Start Test: %s ===================='%self.tag) @@ -593,9 +595,9 @@ def run_opt(self): os.chdir(workdir) return passed - def run_geo(self): + def run_geo(self, with_tsan=False, with_asan=False): - if not self.is_enabled(): + if not self.is_enabled(with_tsan, with_asan): return True print('==================== Start Test: %s ===================='%self.tag) @@ -640,52 +642,59 @@ def run_geo(self): timed_out = True passed = False - # Examine the output - f = open(logfilename,'r') - output = f.readlines() + # check for non-zero return code + process.communicate() + if process.returncode != 0: + passed = False + delta_vals = [] sim_vals = [] - data = [] - if not timed_out: - start_solver = False - for line in output: - if not start_solver: # Don't bother parsing anything before SU2_GEO starts - if line.find('Station 1') > -1: - start_solver=True - elif line.find('Station 2') > -1: # jump out of loop if we hit the next station - break - else: # Found the lines; parse the input - - if line.find('Chord') > -1: - raw_data = line.replace(",", "").split() - data.append(raw_data[1]) - found_chord = True - data.append(raw_data[5]) - found_radius = True - data.append(raw_data[8]) - found_toc = True - data.append(raw_data[10]) - found_aoa = True - - if found_chord and found_radius and found_toc and found_aoa: # Found what we're checking for - iter_missing = False - if not len(self.test_vals)==len(data): # something went wrong... probably bad input - print("Error in test_vals!") - passed = False - for j in range(len(data)): - sim_vals.append( float(data[j]) ) - delta_vals.append( abs(float(data[j])-self.test_vals[j]) ) - if delta_vals[j] > self.tol: - exceed_tol = True - passed = False - else: - iter_missing = True - if not start_solver: - passed = False + if not with_tsan and not with_asan: # sanitizer findings result in non-zero return code, no need to examine the output + # Examine the output + f = open(logfilename,'r') + output = f.readlines() + data = [] + if not timed_out: + start_solver = False + for line in output: + if not start_solver: # Don't bother parsing anything before SU2_GEO starts + if line.find('Station 1') > -1: + start_solver=True + elif line.find('Station 2') > -1: # jump out of loop if we hit the next station + break + else: # Found the lines; parse the input + + if line.find('Chord') > -1: + raw_data = line.replace(",", "").split() + data.append(raw_data[1]) + found_chord = True + data.append(raw_data[5]) + found_radius = True + data.append(raw_data[8]) + found_toc = True + data.append(raw_data[10]) + found_aoa = True + + if found_chord and found_radius and found_toc and found_aoa: # Found what we're checking for + iter_missing = False + if not len(self.test_vals)==len(data): # something went wrong... probably bad input + print("Error in test_vals!") + passed = False + for j in range(len(data)): + sim_vals.append( float(data[j]) ) + delta_vals.append( abs(float(data[j])-self.test_vals[j]) ) + if delta_vals[j] > self.tol: + exceed_tol = True + passed = False + else: + iter_missing = True - if iter_missing: - passed = False + if not start_solver: + passed = False + + if iter_missing: + passed = False # Write the test results #for j in output: @@ -703,20 +712,21 @@ def run_geo(self): if timed_out: print('ERROR: Execution timed out. timeout=%d'%self.timeout) - if exceed_tol: - print('ERROR: Difference between computed input and test_vals exceeded tolerance. TOL=%f'%self.tol) + if not with_tsan and not with_asan: + if exceed_tol: + print('ERROR: Difference between computed input and test_vals exceeded tolerance. TOL=%f'%self.tol) - if not start_solver: - print('ERROR: The code was not able to get to the "OBJFUN" section.') + if not start_solver: + print('ERROR: The code was not able to get to the "OBJFUN" section.') - if iter_missing: - print('ERROR: The SU2_GEO values could not be found.') + if iter_missing: + print('ERROR: The SU2_GEO values could not be found.') - print_vals(self.test_vals, name="test_vals (stored)") + print_vals(self.test_vals, name="test_vals (stored)") - print_vals(sim_vals, name="sim_vals (computed)") + print_vals(sim_vals, name="sim_vals (computed)") - print_vals(delta_vals, name="delta_vals") + print_vals(delta_vals, name="delta_vals") print('test duration: %.2f min'%(running_time/60.0)) print('==================== End Test: %s ====================\n'%self.tag) @@ -725,9 +735,9 @@ def run_geo(self): os.chdir(workdir) return passed - def run_def(self): + def run_def(self, with_tsan=False, with_asan=False): - if not self.is_enabled(): + if not self.is_enabled(with_tsan, with_asan): return True print('==================== Start Test: %s ===================='%self.tag) @@ -767,48 +777,55 @@ def run_def(self): timed_out = True passed = False - # Examine the output - f = open(logfilename,'r') - output = f.readlines() + # check for non-zero return code + process.communicate() + if process.returncode != 0: + passed = False + delta_vals = [] sim_vals = [] - if not timed_out: - start_solver = False - for line in output: - if not start_solver: # Don't bother parsing anything before -- Volumetric grid deformation --- - if line.find('Volumetric grid deformation') > -1: - start_solver=True - else: # Found the -- Volumetric grid deformation --- line; parse the input - raw_data = line.split() - try: - iter_number = int(raw_data[0]) - data = raw_data[len(raw_data)-1:] # Take the last column for comparison - except ValueError: - continue - except IndexError: - continue - if iter_number == self.test_iter: # Found the iteration number we're checking for - iter_missing = False - if not len(self.test_vals)==len(data): # something went wrong... probably bad input - print("Error in test_vals!") - passed = False - break - for j in range(len(data)): - sim_vals.append( float(data[j]) ) - delta_vals.append( abs(float(data[j])-self.test_vals[j]) ) - if delta_vals[j] > self.tol: - exceed_tol = True + if not with_tsan and not with_asan: # sanitizer findings result in non-zero return code, no need to examine the output + # Examine the output + f = open(logfilename,'r') + output = f.readlines() + if not timed_out: + start_solver = False + for line in output: + if not start_solver: # Don't bother parsing anything before -- Volumetric grid deformation --- + if line.find('Volumetric grid deformation') > -1: + start_solver=True + else: # Found the -- Volumetric grid deformation --- line; parse the input + raw_data = line.split() + try: + iter_number = int(raw_data[0]) + data = raw_data[len(raw_data)-1:] # Take the last column for comparison + except ValueError: + continue + except IndexError: + continue + + if iter_number == self.test_iter: # Found the iteration number we're checking for + iter_missing = False + if not len(self.test_vals)==len(data): # something went wrong... probably bad input + print("Error in test_vals!") passed = False - break - else: - iter_missing = True + break + for j in range(len(data)): + sim_vals.append( float(data[j]) ) + delta_vals.append( abs(float(data[j])-self.test_vals[j]) ) + if delta_vals[j] > self.tol: + exceed_tol = True + passed = False + break + else: + iter_missing = True - if not start_solver: - passed = False + if not start_solver: + passed = False - if iter_missing: - passed = False + if iter_missing: + passed = False # Write the test results #for j in output: @@ -826,22 +843,23 @@ def run_def(self): if timed_out: print('ERROR: Execution timed out. timeout=%d sec'%self.timeout) - if exceed_tol: - print('ERROR: Difference between computed input and test_vals exceeded tolerance. TOL=%e'%self.tol) + if not with_tsan and not with_asan: + if exceed_tol: + print('ERROR: Difference between computed input and test_vals exceeded tolerance. TOL=%e'%self.tol) - if not start_solver: - print('ERROR: The code was not able to get to the "Begin solver" section.') + if not start_solver: + print('ERROR: The code was not able to get to the "Begin solver" section.') - if iter_missing: - print('ERROR: The iteration number %d could not be found.'%self.test_iter) + if iter_missing: + print('ERROR: The iteration number %d could not be found.'%self.test_iter) - print('test_iter=%d' % self.test_iter) + print('test_iter=%d' % self.test_iter) - print_vals(self.test_vals, name="test_vals (stored)") + print_vals(self.test_vals, name="test_vals (stored)") - print_vals(sim_vals, name="sim_vals (computed)") + print_vals(sim_vals, name="sim_vals (computed)") - print_vals(delta_vals, name="delta_vals") + print_vals(delta_vals, name="delta_vals") print('test duration: %.2f min'%(running_time/60.0)) #print('==================== End Test: %s ====================\n'%self.tag) @@ -850,7 +868,7 @@ def run_def(self): os.chdir(workdir) return passed - def adjust_iter(self, running_with_tsan=False): + def adjust_iter(self, with_tsan=False, with_asan=False): # Read the cfg file workdir = os.getcwd() @@ -861,7 +879,7 @@ def adjust_iter(self, running_with_tsan=False): new_iter = self.test_iter + 1 - if running_with_tsan: + if with_tsan or with_asan: # detect restart restart_iter = 0 @@ -942,18 +960,19 @@ def disable_restart(self): return - def is_enabled(self, running_with_tsan=False): + def is_enabled(self, with_tsan=False, with_asan=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)) - tsan_compatible = not running_with_tsan or self.enabled_with_tsan + tsan_compatible = not with_tsan or self.enabled_with_tsan + asan_compatible = not with_asan or self.enabled_with_asan 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 + return is_enabled_on_arch and tsan_compatible and asan_compatible def adjust_test_data(self): diff --git a/TestCases/flamelet/07_laminar_premixed_h2_flame_cfd/laminar_premixed_h2_flame_cfd.cfg b/TestCases/flamelet/07_laminar_premixed_h2_flame_cfd/laminar_premixed_h2_flame_cfd.cfg new file mode 100644 index 00000000000..3f294865a17 --- /dev/null +++ b/TestCases/flamelet/07_laminar_premixed_h2_flame_cfd/laminar_premixed_h2_flame_cfd.cfg @@ -0,0 +1,120 @@ +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% % +% SU2 configuration file % +% Case description: Laminar premixed hydrogen flame with heat exchanger % +% Author: Evert Bunschoten % +% Institution: Delft University of Technology % +% Date: 01/11/2023 % +% % +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +% ------------- DIRECT, ADJOINT, AND LINEARIZED PROBLEM DEFINITION ------------% +% +SOLVER = INC_NAVIER_STOKES +KIND_TURB_MODEL= NONE +MATH_PROBLEM= DIRECT +RESTART_SOL =YES +% ---------------- INCOMPRESSIBLE FLOW CONDITION DEFINITION -------------------% +% +INC_DENSITY_MODEL= VARIABLE +INC_ENERGY_EQUATION = YES +INC_VELOCITY_INIT= (1.13, 0.0, 0.0 ) +INC_TEMPERATURE_INIT= 300.0 +THERMODYNAMIC_PRESSURE= 101325 +INC_NONDIM= DIMENSIONAL + +% -------------------- FLUID MODEL --------------------------------------- % +% +FLUID_MODEL= FLUID_FLAMELET +VISCOSITY_MODEL= FLAMELET +CONDUCTIVITY_MODEL= FLAMELET +DIFFUSIVITY_MODEL= FLAMELET +KIND_SCALAR_MODEL= FLAMELET +INTERPOLATION_METHOD= MLP +FILENAMES_INTERPOLATOR= (MLP_TD.mlp, MLP_PD.mlp, MLP_PPV.mlp, MLP_null.mlp) +PREFERENTIAL_DIFFUSION= YES + +% -------------------- SCALAR TRANSPORT ---------------------------------------% +% +% Using an artificial spark to ignite the solution at some location and iteration +FLAME_INIT_METHOD= SPARK +% Spark parameters in order: +% x-location of spark center (m) +% y-location of spark center (m) +% z-location of spark center (m) +% Spark radius (m) +% Iteration at which the artificial spark starts +% Spark iteration duration +SPARK_INIT= (0.001, 0.0004, 0, 5e-4, 10, 10) + +% Controlling variable source terms applied within the spark sphere for the spark +% duration. +SPARK_REACTION_RATES= (1000, 0, 0) + +SPECIES_INIT = (-0.49904325357252965, 2226.901776784524, 0.01446751783896619) + +% Passive reactants in flamelet problem + +CONTROLLING_VARIABLE_NAMES= (ProgressVariable, EnthalpyTot, MixtureFraction) +CONTROLLING_VARIABLE_SOURCE_NAMES= (ProdRateTot_PV, NULL, NULL) + +SPECIES_CLIPPING=YES +SPECIES_CLIPPING_MAX=+4.5623852432084366e-01 +8.6731375409855954e+06 1.0 +SPECIES_CLIPPING_MIN= -6.8059708053507162e-01 -4.9308262569627967e+06 0.0 + +MARKER_INLET_SPECIES = (inlet, -0.49904325357252965, 2226.901776784524, 0.01446751783896619) + +% -------------------- BOUNDARY CONDITION DEFINITION --------------------------% +% +MARKER_ISOTHERMAL= (burner_wall, 300, hex_wall, 300) +MARKER_SYM= (sides) +INC_INLET_TYPE= VELOCITY_INLET +MARKER_INLET = (inlet, 300.000, 0.575, 1.0, 0.0, 0.0) +INC_OUTLET_TYPE= PRESSURE_OUTLET +MARKER_OUTLET= (outlet, 0.0) +MARKER_ANALYZE_AVERAGE = AREA + +% ------------- COMMON PARAMETERS DEFINING THE NUMERICAL METHOD ---------------% +% +NUM_METHOD_GRAD= GREEN_GAUSS +CFL_NUMBER =50 +CFL_ADAPT= NO +ITER=5 +OUTPUT_WRT_FREQ= 20 +% ------------------------ LINEAR SOLVER DEFINITION ---------------------------% +% +LINEAR_SOLVER= FGMRES +LINEAR_SOLVER_PREC= ILU +LINEAR_SOLVER_ERROR= 1E-4 +LINEAR_SOLVER_ITER=20 +% +% -------------------- FLOW AND SPECIES NUMERICAL METHOD DEFINITION -----------------------% +% +CONV_NUM_METHOD_FLOW= FDS +CONV_NUM_METHOD_SPECIES= BOUNDED_SCALAR +MUSCL_FLOW= YES +MUSCL_SPECIES= YES +SLOPE_LIMITER_FLOW = NONE +SLOPE_LIMITER_SPECIES= NONE +TIME_DISCRE_FLOW= EULER_IMPLICIT +TIME_DISCRE_SPECIES= EULER_IMPLICIT + +% --------------------------- CONVERGENCE PARAMETERS --------------------------% +% +CONV_FIELD = RMS_EnthalpyTot +CONV_RESIDUAL_MINVAL= -5 +CONV_STARTITER= 20 +SCREEN_OUTPUT = INNER_ITER RMS_PRESSURE RMS_ProgressVariable RMS_EnthalpyTot RMS_MixtureFraction +HISTORY_OUTPUT = WALL_TIME RMS_RES +VOLUME_OUTPUT = SOLUTION + +% ------------------------- INPUT/OUTPUT INFORMATION --------------------------% +% +MESH_FORMAT= SU2 +MESH_FILENAME = 2Dhex_BL.su2 +OUTPUT_FILES = (RESTART,PARAVIEW) +TABULAR_FORMAT = CSV +CONV_FILENAME= history +VOLUME_FILENAME= flow + + diff --git a/TestCases/parallel_regression.py b/TestCases/parallel_regression.py index 51b18e48741..f8403488070 100644 --- a/TestCases/parallel_regression.py +++ b/TestCases/parallel_regression.py @@ -69,6 +69,15 @@ def main(): cfd_flamelet_ch4_partial_premix.new_output = True test_list.append(cfd_flamelet_ch4_partial_premix) + # 2D planar laminar premixed hydrogen flame on isothermal burner with heat exchanger emulator (restart) + cfd_flamelet_h2 = TestCase('cfd_flamelet_h2') + cfd_flamelet_h2.cfg_dir = "flamelet/07_laminar_premixed_h2_flame_cfd" + cfd_flamelet_h2.cfg_file = "laminar_premixed_h2_flame_cfd.cfg" + cfd_flamelet_h2.test_iter = 5 + cfd_flamelet_h2.test_vals = [-10.003106, -9.843748, -3.289857, -11.338273] + cfd_flamelet_h2.new_output = True + test_list.append(cfd_flamelet_h2) + ######################### ## NEMO solver ### ######################### @@ -1045,10 +1054,11 @@ def main(): ###################################### # Aachen Turbine restart - Aachen_3D_restart = TestCase('aachen_turbine_restart') - Aachen_3D_restart.cfg_dir = "turbomachinery/Aachen_turbine" - Aachen_3D_restart.cfg_file = "aachen_3D_MP_restart.cfg" + Aachen_3D_restart = TestCase('aachen_turbine_restart') + 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.enabled_with_asan = False Aachen_3D_restart.test_vals = [-15.329206, -15.008622, -15.078888, -13.841072, -12.727840, -9.975729] test_list.append(Aachen_3D_restart) @@ -1065,8 +1075,8 @@ def main(): axial_stage2D.cfg_dir = "turbomachinery/axial_stage_2D" axial_stage2D.cfg_file = "Axial_stage2D.cfg" axial_stage2D.test_iter = 20 - axial_stage2D.test_vals = [0.983754, 1.534455, -2.888523, 2.606770, -2.418403, 3.087203, 106380, 106380, 5.7325, 64.711] - axial_stage2D.test_vals_aarch64 = [0.983754, 1.534455, -2.888523, 2.606770, -2.418403, 3.087203, 106380, 106380, 5.7325, 64.711] + axial_stage2D.test_vals = [0.983754, 1.534455, -2.888523, 2.606770, -2.418403, 3.087203, 106380, 106380, 5.7325, 64.711] + axial_stage2D.test_vals_aarch64 = [0.983754, 1.534455, -2.888523, 2.606770, -2.418403, 3.087203, 106380, 106380, 5.7325, 64.711] test_list.append(axial_stage2D) # 2D transonic stator restart @@ -1074,7 +1084,7 @@ def main(): transonic_stator_restart.cfg_dir = "turbomachinery/transonic_stator_2D" transonic_stator_restart.cfg_file = "transonic_stator_restart.cfg" transonic_stator_restart.test_iter = 20 - transonic_stator_restart.test_vals = [-5.011834, -3.091110, -2.757795, 1.087934, -3.544707, 2.166101, -471630, 94.868, -0.035888] + transonic_stator_restart.test_vals = [-5.011834, -3.091110, -2.757795, 1.087934, -3.544707, 2.166101, -471630, 94.868, -0.035888] transonic_stator_restart.test_vals_aarch64 = [-5.011834, -3.091110, -2.757795, 1.087934, -3.544707, 2.166101, -471630, 94.868, -0.035888] test_list.append(transonic_stator_restart) @@ -1378,7 +1388,7 @@ def main(): pywrapper_unsteadyFSI.test_iter = 4 pywrapper_unsteadyFSI.test_vals = [0, 31, 5, 58, -1.756780, -2.828276, -7.652558, -6.863929, 1.5618e-04] pywrapper_unsteadyFSI.command = TestCase.Command("mpirun -np 2", "python", "run.py") - pywrapper_unsteadyFSI.unsteady = True + pywrapper_unsteadyFSI.unsteady = True pywrapper_unsteadyFSI.multizone = True test_list.append(pywrapper_unsteadyFSI) @@ -1570,13 +1580,13 @@ def main(): ##################### # CGNS writer - cgns_writer = TestCase('cgns_writer') - cgns_writer.cfg_dir = "cgns_writer" - cgns_writer.cfg_file = "config.cfg" - cgns_writer.test_iter = 1 - cgns_writer.test_vals = [-2.974473, 0.665204, 5.068846, -7.003873] - cgns_writer.command = TestCase.Command("mpirun -n 2", "SU2_CFD") - cgns_writer.new_output = True + cgns_writer = TestCase('cgns_writer') + cgns_writer.cfg_dir = "cgns_writer" + cgns_writer.cfg_file = "config.cfg" + cgns_writer.test_iter = 1 + cgns_writer.test_vals = [-2.974473, 0.665204, 5.068846, -7.003873] + cgns_writer.command = TestCase.Command("mpirun -n 2", "SU2_CFD") + cgns_writer.new_output = True test_list.append(cgns_writer) ###################################### @@ -1806,14 +1816,14 @@ def main(): test_list.append(sphere_ffd_def_bspline) # Inviscid NACA0012 (triangles) - naca0012_cst = TestCase('naca0012_cst') - naca0012_cst.cfg_dir = "deformation/cst" - naca0012_cst.cfg_file = "naca0012.cfg" + naca0012_cst = TestCase('naca0012_cst') + naca0012_cst.cfg_dir = "deformation/cst" + naca0012_cst.cfg_file = "naca0012.cfg" naca0012_cst.test_iter = 10 naca0012_cst.test_vals = [0.000385514] #residual - naca0012_cst.command = TestCase.Command("mpirun -n 2", "SU2_DEF") - naca0012_cst.timeout = 1600 - naca0012_cst.tol = 1e-8 + naca0012_cst.command = TestCase.Command("mpirun -n 2", "SU2_DEF") + naca0012_cst.timeout = 1600 + naca0012_cst.tol = 1e-8 pass_list.append(naca0012_cst.run_def()) test_list.append(naca0012_cst) diff --git a/TestCases/serial_regression.py b/TestCases/serial_regression.py index dab5dd1256f..f7ab332ae7a 100644 --- a/TestCases/serial_regression.py +++ b/TestCases/serial_regression.py @@ -28,12 +28,15 @@ from __future__ import print_function, division, absolute_import import sys from TestCase import TestCase +from TestCase import parse_args def main(): '''This program runs SU2 and ensures that the output matches specified values. This will be used to do checks when code is pushed to github to make sure nothing is broken. ''' + args = parse_args('Serial Regression Tests') + test_list = [] ######################### @@ -857,6 +860,7 @@ def main(): Aachen_3D_restart.cfg_dir = "turbomachinery/Aachen_turbine" Aachen_3D_restart.cfg_file = "aachen_3D_MP_restart.cfg" Aachen_3D_restart.test_iter = 5 + Aachen_3D_restart.enabled_with_asan = False Aachen_3D_restart.test_vals = [-15.137167, -14.551444, -15.078894, -13.486154, -12.724891, -9.717612] test_list.append(Aachen_3D_restart) @@ -1196,7 +1200,7 @@ def main(): if test.tol == 0.0: test.tol = 0.00001 - pass_list = [ test.run_test() for test in test_list ] + pass_list = [ test.run_test(args.tsan, args.asan) for test in test_list ] ###################################### @@ -1211,7 +1215,7 @@ def main(): naca0012_geo.command = TestCase.Command(exec = "SU2_GEO") naca0012_geo.timeout = 1600 naca0012_geo.tol = 0.00001 - pass_list.append(naca0012_geo.run_geo()) + pass_list.append(naca0012_geo.run_geo(args.tsan, args.asan)) test_list.append(naca0012_geo) ###################################### @@ -1228,7 +1232,7 @@ def main(): intersect_def.timeout = 1600 intersect_def.tol = 1e-04 - pass_list.append(intersect_def.run_def()) + pass_list.append(intersect_def.run_def(args.tsan, args.asan)) test_list.append(intersect_def) # Inviscid NACA0012 (triangles) @@ -1241,7 +1245,7 @@ def main(): naca0012_def.timeout = 1600 naca0012_def.tol = 1e-08 - pass_list.append(naca0012_def.run_def()) + pass_list.append(naca0012_def.run_def(args.tsan, args.asan)) test_list.append(naca0012_def) # Inviscid NACA0012 based on SURFACE_FILE input (surface_bump.dat) @@ -1254,7 +1258,7 @@ def main(): naca0012_def_file.timeout = 1600 naca0012_def_file.tol = 1e-8 - pass_list.append(naca0012_def_file.run_def()) + pass_list.append(naca0012_def_file.run_def(args.tsan, args.asan)) test_list.append(naca0012_def_file) # RAE2822 (mixed tris + quads) @@ -1267,7 +1271,7 @@ def main(): rae2822_def.timeout = 1600 rae2822_def.tol = 1e-13 - pass_list.append(rae2822_def.run_def()) + pass_list.append(rae2822_def.run_def(args.tsan, args.asan)) test_list.append(rae2822_def) # Turb NACA4412 (quads, wall distance) @@ -1280,7 +1284,7 @@ def main(): naca4412_def.timeout = 1600 naca4412_def.tol = 1e-12 - pass_list.append(naca4412_def.run_def()) + pass_list.append(naca4412_def.run_def(args.tsan, args.asan)) test_list.append(naca4412_def) # Brick of tets (inverse volume) @@ -1293,7 +1297,7 @@ def main(): brick_tets_def.timeout = 1600 brick_tets_def.tol = 1e-09 - pass_list.append(brick_tets_def.run_def()) + pass_list.append(brick_tets_def.run_def(args.tsan, args.asan)) test_list.append(brick_tets_def) # Brick of isotropic hexas (inverse volume) @@ -1306,7 +1310,7 @@ def main(): brick_hex_def.timeout = 1600 brick_hex_def.tol = 1e-09 - pass_list.append(brick_hex_def.run_def()) + pass_list.append(brick_hex_def.run_def(args.tsan, args.asan)) test_list.append(brick_hex_def) # Brick with a pyramid layer (inverse volume) @@ -1319,7 +1323,7 @@ def main(): brick_pyra_def.timeout = 1600 brick_pyra_def.tol = 1e-08 - pass_list.append(brick_pyra_def.run_def()) + pass_list.append(brick_pyra_def.run_def(args.tsan, args.asan)) test_list.append(brick_pyra_def) # Brick of isotropic prisms (inverse volume) @@ -1332,7 +1336,7 @@ def main(): brick_prism_def.timeout = 1600 brick_prism_def.tol = 1e-08 - pass_list.append(brick_prism_def.run_def()) + pass_list.append(brick_prism_def.run_def(args.tsan, args.asan)) test_list.append(brick_prism_def) # Brick of prisms with high aspect ratio cells near the wall (wall distance) @@ -1345,7 +1349,7 @@ def main(): brick_prism_rans_def.timeout = 1600 brick_prism_rans_def.tol = 1e-12 - pass_list.append(brick_prism_rans_def.run_def()) + pass_list.append(brick_prism_rans_def.run_def(args.tsan, args.asan)) test_list.append(brick_prism_rans_def) # Brick of hexas with high aspect ratio cells near the wall (inverse volume) @@ -1358,7 +1362,7 @@ def main(): brick_hex_rans_def.timeout = 1600 brick_hex_rans_def.tol = 1e-12 - pass_list.append(brick_hex_rans_def.run_def()) + pass_list.append(brick_hex_rans_def.run_def(args.tsan, args.asan)) test_list.append(brick_hex_rans_def) # Cylindrical FFD test @@ -1371,7 +1375,7 @@ def main(): cylinder_ffd_def.timeout = 1600 cylinder_ffd_def.tol = 1e-09 - pass_list.append(cylinder_ffd_def.run_def()) + pass_list.append(cylinder_ffd_def.run_def(args.tsan, args.asan)) test_list.append(cylinder_ffd_def) # Spherical FFD test @@ -1384,7 +1388,7 @@ def main(): sphere_ffd_def.timeout = 1600 sphere_ffd_def.tol = 1e-08 - pass_list.append(sphere_ffd_def.run_def()) + pass_list.append(sphere_ffd_def.run_def(args.tsan, args.asan)) test_list.append(sphere_ffd_def) # Spherical FFD test using BSplines @@ -1397,7 +1401,7 @@ def main(): sphere_ffd_def_bspline.timeout = 1600 sphere_ffd_def_bspline.tol = 1e-08 - pass_list.append(sphere_ffd_def_bspline.run_def()) + pass_list.append(sphere_ffd_def_bspline.run_def(args.tsan, args.asan)) test_list.append(sphere_ffd_def_bspline) ###################################### @@ -1413,7 +1417,8 @@ def main(): contadj_euler_py.timeout = 1600 contadj_euler_py.reference_file = "of_grad_cd.dat.ref" contadj_euler_py.test_file = "of_grad_cd.dat" - pass_list.append(contadj_euler_py.run_filediff()) + contadj_euler_py.enabled_with_asan = False + pass_list.append(contadj_euler_py.run_filediff(args.tsan, args.asan)) test_list.append(contadj_euler_py) # test shape_optimization.py @@ -1425,7 +1430,8 @@ def main(): shape_opt_euler_py.command = TestCase.Command(exec = "shape_optimization.py", param = "-g CONTINUOUS_ADJOINT -f") shape_opt_euler_py.timeout = 1600 shape_opt_euler_py.tol = 0.00001 - pass_list.append(shape_opt_euler_py.run_opt()) + shape_opt_euler_py.enabled_with_asan = False + pass_list.append(shape_opt_euler_py.run_opt(args.tsan, args.asan)) test_list.append(shape_opt_euler_py) # Multiple functionals with the continuous adjoint @@ -1437,7 +1443,8 @@ def main(): contadj_multi_py.timeout = 1600 contadj_multi_py.reference_file = "of_grad_combo.dat.ref" contadj_multi_py.test_file = "of_grad_combo.dat" - pass_list.append(contadj_multi_py.run_filediff()) + contadj_multi_py.enabled_with_asan = False + pass_list.append(contadj_multi_py.run_filediff(args.tsan, args.asan)) test_list.append(contadj_multi_py) # Optimization with multiple objectives, with gradients evaluated individually @@ -1475,7 +1482,8 @@ def main(): opt_multiobj1surf_py.command = TestCase.Command(exec = "shape_optimization.py", param = "-g CONTINUOUS_ADJOINT -f") opt_multiobj1surf_py.timeout = 1600 opt_multiobj1surf_py.tol = 0.00001 - pass_list.append(opt_multiobj1surf_py.run_opt()) + opt_multiobj1surf_py.enabled_with_asan = False + pass_list.append(opt_multiobj1surf_py.run_opt(args.tsan, args.asan)) test_list.append(opt_multiobj1surf_py) # test optimization, with a single objective evaluated on multiple surfaces @@ -1487,7 +1495,8 @@ def main(): opt_2surf1obj_py.command = TestCase.Command(exec = "shape_optimization.py", param = "-g CONTINUOUS_ADJOINT -f") opt_2surf1obj_py.timeout = 1600 opt_2surf1obj_py.tol = 0.00001 - pass_list.append(opt_2surf1obj_py.run_opt()) + opt_2surf1obj_py.enabled_with_asan = False + pass_list.append(opt_2surf1obj_py.run_opt(args.tsan, args.asan)) test_list.append(opt_2surf1obj_py) ########################## @@ -1503,8 +1512,9 @@ def main(): pywrapper_naca0012.command = TestCase.Command(exec = "SU2_CFD.py", param = "-f") pywrapper_naca0012.timeout = 1600 pywrapper_naca0012.tol = 0.00001 + pywrapper_naca0012.enabled_with_asan = False test_list.append(pywrapper_naca0012) - pass_list.append(pywrapper_naca0012.run_test()) + pass_list.append(pywrapper_naca0012.run_test(args.tsan, args.asan)) # NACA0012 (SST, FUN3D results for finest grid: CL=1.0840, CD=0.01253) pywrapper_turb_naca0012_sst = TestCase('pywrapper_turb_naca0012_sst') @@ -1516,8 +1526,9 @@ def main(): pywrapper_turb_naca0012_sst.command = TestCase.Command(exec = "SU2_CFD.py", param = "-f") pywrapper_turb_naca0012_sst.timeout = 3200 pywrapper_turb_naca0012_sst.tol = 0.00001 + pywrapper_turb_naca0012_sst.enabled_with_asan = False test_list.append(pywrapper_turb_naca0012_sst) - pass_list.append(pywrapper_turb_naca0012_sst.run_test()) + pass_list.append(pywrapper_turb_naca0012_sst.run_test(args.tsan, args.asan)) # Square cylinder pywrapper_square_cylinder = TestCase('pywrapper_square_cylinder') @@ -1529,8 +1540,9 @@ def main(): pywrapper_square_cylinder.timeout = 1600 pywrapper_square_cylinder.tol = 0.00001 pywrapper_square_cylinder.unsteady = True + pywrapper_square_cylinder.enabled_with_asan = False test_list.append(pywrapper_square_cylinder) - pass_list.append(pywrapper_square_cylinder.run_test()) + pass_list.append(pywrapper_square_cylinder.run_test(args.tsan, args.asan)) # Aeroelastic pywrapper_aeroelastic = TestCase('pywrapper_aeroelastic') @@ -1542,8 +1554,9 @@ def main(): pywrapper_aeroelastic.timeout = 1600 pywrapper_aeroelastic.tol = 0.00001 pywrapper_aeroelastic.unsteady = True + pywrapper_aeroelastic.enabled_with_asan = False test_list.append(pywrapper_aeroelastic) - pass_list.append(pywrapper_aeroelastic.run_test()) + pass_list.append(pywrapper_aeroelastic.run_test(args.tsan, args.asan)) # FSI, 2d pywrapper_fsi2d = TestCase('pywrapper_fsi2d') @@ -1556,8 +1569,9 @@ def main(): pywrapper_fsi2d.multizone = True pywrapper_fsi2d.timeout = 1600 pywrapper_fsi2d.tol = 0.00001 + pywrapper_fsi2d.enabled_with_asan = False test_list.append(pywrapper_fsi2d) - pass_list.append(pywrapper_fsi2d.run_test()) + pass_list.append(pywrapper_fsi2d.run_test(args.tsan, args.asan)) # Unsteady CHT pywrapper_unsteadyCHT = TestCase('pywrapper_unsteadyCHT') @@ -1569,8 +1583,9 @@ def main(): pywrapper_unsteadyCHT.timeout = 1600 pywrapper_unsteadyCHT.tol = 0.00001 pywrapper_unsteadyCHT.unsteady = True + pywrapper_unsteadyCHT.enabled_with_asan = False test_list.append(pywrapper_unsteadyCHT) - pass_list.append(pywrapper_unsteadyCHT.run_test()) + pass_list.append(pywrapper_unsteadyCHT.run_test(args.tsan, args.asan)) # Rigid motion pywrapper_rigidMotion = TestCase('pywrapper_rigidMotion') @@ -1582,8 +1597,9 @@ def main(): pywrapper_rigidMotion.timeout = 1600 pywrapper_rigidMotion.tol = 0.00001 pywrapper_rigidMotion.unsteady = True + pywrapper_rigidMotion.enabled_with_asan = False test_list.append(pywrapper_rigidMotion) - pass_list.append(pywrapper_rigidMotion.run_test()) + pass_list.append(pywrapper_rigidMotion.run_test(args.tsan, args.asan)) # Tests summary print('==================================================================') diff --git a/TestCases/serial_regression_AD.py b/TestCases/serial_regression_AD.py index 2a0ac757097..ac268040cd5 100644 --- a/TestCases/serial_regression_AD.py +++ b/TestCases/serial_regression_AD.py @@ -30,12 +30,15 @@ import sys from TestCase import TestCase +from TestCase import parse_args def main(): '''This program runs SU2 and ensures that the output matches specified values. This will be used to do checks when code is pushed to github to make sure nothing is broken. ''' + args = parse_args('Serial Regression AD Tests') + test_list = [] ##################################### @@ -223,7 +226,7 @@ def main(): if test.tol == 0.0: test.tol = 0.00001 - pass_list = [ test.run_test() for test in test_list ] + pass_list = [ test.run_test(args.tsan, args.asan) for test in test_list ] ################################### ### Coupled RHT-CFD Adjoint ### @@ -239,7 +242,8 @@ def main(): discadj_rht.reference_file = "of_grad_cd.csv.ref" discadj_rht.reference_file_aarch64 = "of_grad_cd_aarch64.csv.ref" discadj_rht.test_file = "of_grad_cd.csv" - pass_list.append(discadj_rht.run_filediff()) + discadj_rht.enabled_with_asan = False + pass_list.append(discadj_rht.run_filediff(args.tsan, args.asan)) test_list.append(discadj_rht) ###################################### @@ -256,7 +260,8 @@ def main(): discadj_euler_py.reference_file = "of_grad_cd_disc.dat.ref" discadj_euler_py.reference_file_aarch64 = "of_grad_cd_disc_aarch64.dat.ref" discadj_euler_py.test_file = "of_grad_cd.dat" - pass_list.append(discadj_euler_py.run_filediff()) + discadj_euler_py.enabled_with_asan = False + pass_list.append(discadj_euler_py.run_filediff(args.tsan, args.asan)) test_list.append(discadj_euler_py) # test discrete_adjoint with multiple ffd boxes @@ -269,7 +274,8 @@ def main(): discadj_multiple_ffd_py.reference_file = "of_grad_cd.dat.ref" discadj_multiple_ffd_py.reference_file_aarch64 = "of_grad_cd_aarch64.dat.ref" discadj_multiple_ffd_py.test_file = "of_grad_cd.dat" - pass_list.append(discadj_multiple_ffd_py.run_filediff()) + discadj_multiple_ffd_py.enabled_with_asan = False + pass_list.append(discadj_multiple_ffd_py.run_filediff(args.tsan, args.asan)) test_list.append(discadj_multiple_ffd_py) # test direct_differentiation.py @@ -282,7 +288,8 @@ def main(): directdiff_euler_py.reference_file = "of_grad_directdiff.dat.ref" directdiff_euler_py.reference_file_aarch64 = "of_grad_directdiff_aarch64.dat.ref" directdiff_euler_py.test_file = "DIRECTDIFF/of_grad_directdiff.dat" - pass_list.append(directdiff_euler_py.run_filediff()) + directdiff_euler_py.enabled_with_asan = False + pass_list.append(directdiff_euler_py.run_filediff(args.tsan, args.asan)) test_list.append(directdiff_euler_py) # test direct_differentiation.py with multiple ffd boxes @@ -295,7 +302,8 @@ def main(): directdiff_multiple_ffd_py.reference_file = "of_grad_directdiff.dat.ref" directdiff_multiple_ffd_py.reference_file_aarch64 = "of_grad_directdiff_aarch64.dat.ref" directdiff_multiple_ffd_py.test_file = "DIRECTDIFF/of_grad_directdiff.dat" - pass_list.append(directdiff_multiple_ffd_py.run_filediff()) + directdiff_multiple_ffd_py.enabled_with_asan = False + pass_list.append(directdiff_multiple_ffd_py.run_filediff(args.tsan, args.asan)) test_list.append(directdiff_multiple_ffd_py) # test continuous_adjoint.py, with multiple objectives @@ -320,8 +328,9 @@ def main(): pywrapper_FEA_AD_FlowLoad.timeout = 1600 pywrapper_FEA_AD_FlowLoad.tol = 0.000001 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()) + pass_list.append(pywrapper_FEA_AD_FlowLoad.run_test(args.tsan, args.asan)) # Flow AD Mesh Displacement Sensitivity pywrapper_CFD_AD_MeshDisp = TestCase('pywrapper_CFD_AD_MeshDisp') @@ -333,8 +342,9 @@ def main(): pywrapper_CFD_AD_MeshDisp.timeout = 1600 pywrapper_CFD_AD_MeshDisp.tol = 0.000001 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()) + pass_list.append(pywrapper_CFD_AD_MeshDisp.run_test(args.tsan, args.asan)) ################################### @@ -350,7 +360,8 @@ def main(): grad_smooth_naca0012.reference_file = "of_hess.dat.ref" grad_smooth_naca0012.reference_file_aarch64 = "of_hess_aarch64.dat.ref" grad_smooth_naca0012.test_file = "of_hess.dat" - pass_list.append(grad_smooth_naca0012.run_filediff()) + grad_smooth_naca0012.enabled_with_asan = False + pass_list.append(grad_smooth_naca0012.run_filediff(args.tsan, args.asan)) test_list.append(grad_smooth_naca0012) # Tests summary diff --git a/TestCases/tutorials.py b/TestCases/tutorials.py index 849a4f85c26..f63fb81c990 100644 --- a/TestCases/tutorials.py +++ b/TestCases/tutorials.py @@ -120,6 +120,16 @@ def main(): sudo_tutorial.command = TestCase.Command("mpirun -n 2", "SU2_CFD") test_list.append(sudo_tutorial) + ### Incompressible Combustion + + # Pre-mixed, laminar hydrogen flame with heat loss + premixed_hydrogen = TestCase('premixed_hydrogen') + premixed_hydrogen.cfg_dir = "../Tutorials/incompressible_flow/Inc_Combustion/1__premixed_hydrogen" + premixed_hydrogen.cfg_file = "H2_burner.cfg" + premixed_hydrogen.test_iter = 10 + premixed_hydrogen.test_vals = [-9.905856, -10.449512, -11.035999, -4.331440, -11.882740] + test_list.append(premixed_hydrogen) + ### Compressible Flow # Inviscid Bump diff --git a/UnitTests/SU2_CFD/gradients.cpp b/UnitTests/SU2_CFD/gradients.cpp index e0829ce1104..2f69d407d3a 100644 --- a/UnitTests/SU2_CFD/gradients.cpp +++ b/UnitTests/SU2_CFD/gradients.cpp @@ -132,8 +132,8 @@ void testGreenGauss() { TestField field; C3DDoubleMatrix gradient(field.geometry->GetnPoint(), field.nVar, field.geometry->GetnDim()); - computeGradientsGreenGauss(nullptr, SOLUTION, PERIODIC_NONE, *field.geometry.get(), *field.config.get(), field, 0, - field.nVar, gradient); + computeGradientsGreenGauss(nullptr, MPI_QUANTITIES::SOLUTION, PERIODIC_NONE, *field.geometry.get(), + *field.config.get(), field, 0, field.nVar, gradient); check(field, gradient); } @@ -144,8 +144,8 @@ void testLeastSquares(bool weighted) { C3DDoubleMatrix R(field.geometry->GetnPoint(), nDim, nDim); C3DDoubleMatrix gradient(field.geometry->GetnPoint(), field.nVar, nDim); - computeGradientsLeastSquares(nullptr, SOLUTION, PERIODIC_NONE, *field.geometry.get(), *field.config.get(), weighted, - field, 0, field.nVar, gradient, R); + computeGradientsLeastSquares(nullptr, MPI_QUANTITIES::SOLUTION, PERIODIC_NONE, *field.geometry.get(), + *field.config.get(), weighted, field, 0, field.nVar, gradient, R); check(field, gradient); } diff --git a/config_template.cfg b/config_template.cfg index bf3e6bae907..74832ef98db 100644 --- a/config_template.cfg +++ b/config_template.cfg @@ -869,6 +869,11 @@ SPECIES_CLIPPING_MAX= 1.0 SPECIES_CLIPPING_MIN= 0.0 % --------------------- FLAMELET MODEL -----------------------------% +% The flamelet model uses the prompts INTERPOLATION_METHOD and FILENAMES_INTERPOLATOR +% for definition of the flamelet manifold. Either LUT or MLP can be used as options. +% If the terms "Beta_ProgVar", "Beta_Enth_Thermal", "Beta_Enth", and "Beta_Mixfrac" are +% found in the variables list of the manifold, preferential diffusion is assumed. + % % Names of the user defined (auxiliary) transport equations. USER_SCALAR_NAMES= (Y-CO) @@ -894,8 +899,27 @@ CONTROLLING_VARIABLE_NAMES= (ProgressVariable, EnthalpyTot) % controlling variables without source terms, use "zero" or "null", similar to the % option USER_SOURCE_NAMES. CONTROLLING_VARIABLE_SOURCE_NAMES= (ProdRateTot_PV, NULL) -% -% flamelet initialization + +% Method used to ignite the solution: +% FLAME_FRONT : the flame is initialized using a plane, defined by a point and a normal. On one side, the solution is initialized +% using 'burnt' conditions and on the other side 'unburnt' conditions. The normal points in the direction of the 'burnt' +% condition. +% SPARK : the solution is ignited through application of a set of source terms within a specified region for a set number +% of solver iterations. This artificial spark can be applied after a certain number of iterations has passed, allowing for +% the flow to evolve before igniting the mixture. +% NONE : no artificial solution ignition. +% By default, this option is set to NONE +FLAME_INIT_METHOD= FLAME_FRONT + +% Enable preferential diffusion for progress variable, total enthalpy, mixture fraction FGM problems. +% If 'YES', the preferential diffusion model from Efimov et al.(2021) is used for computing the viscous residual terms. +% If enabled, make sure the variables "Beta_ProgVar", "Beta_Enth_Thermal", "Beta_Enth", and "Beta_MixFrac", corresponding +% to the preferential diffusion terms of the progress variable, temperature, total enthalpy, and mixture fraction +% respectively are included in the manifold. +% By default, this option is disabled. +PREFERENTIAL_DIFFUSION= NO + +% FLAME_FRONT initialization % the flame is initialized using a plane, defined by a point and a normal. On one side, the solution is initialized % using 'burnt' conditions and on the other side 'unburnt' conditions. The normal points in the direction of the 'burnt' % condition. @@ -904,6 +928,21 @@ CONTROLLING_VARIABLE_SOURCE_NAMES= (ProdRateTot_PV, NULL) % (x7) = thickness of the reaction zone, this is the transition from unburnt to burnt conditions. % (x8) = Thickness of the 'burnt' zone, after this length, the conditions will be 'unburnt' again. FLAME_INIT= (0.004, 0.0, 0.0, 1.0, 0.0, 0.0, 0.2e-3, 1.0) + +% SPARK initialization +% the solution is ignited through application of a set of source terms within a specified region for a set number +% of solver iterations. This artificial spark can be applied after a certain number of iterations has passed, allowing for +% the flow to evolve before igniting the mixture. +% (x1,x2,x3) = spark center location. +% (x4) = spark radius where within the artificial spark is applied. +% (x5,x6) = spark iteration start and number of iterations in which the artifical spark is applied. For single-zone problems +% the number of iterations are inner-loop iterations, for multi-zone problems, outer-loop iterations are considered. +SPARK_INIT= (0.004, 0.0, 0.0, 1e-4, 100, 10) + +% Source terms (density times species time rate of change) applied to the species equations in the spark region for the duration of the spark. +% The number of terms should equate the total number of species in the flamelet problem. +SPARK_REACTION_RATES= (1000, 0, 0) + % % --------------------- INVERSE DESIGN SIMULATION -----------------------------% % diff --git a/meson_scripts/init.py b/meson_scripts/init.py index 787bfb56e2f..bc58f9a9564 100755 --- a/meson_scripts/init.py +++ b/meson_scripts/init.py @@ -71,7 +71,7 @@ def init_submodules( github_repo_coolprop = "https://github.com/CoolProp/CoolProp" sha_version_mel = "46205ab019e5224559091375a6d71aabae6bc5b9" github_repo_mel = "https://github.com/pcarruscag/MEL" - sha_version_mlpcpp = "665c45b7d3533c977eb1f637918d5b8b75c07d3b" + sha_version_mlpcpp = "c19c53ea2b85ccfb185f1c6c87044dc0b5bc7ae0" github_repo_mlpcpp = "https://github.com/EvertBunschoten/MLPCpp" medi_name = "MeDiPack" diff --git a/subprojects/MLPCpp b/subprojects/MLPCpp index 665c45b7d35..c19c53ea2b8 160000 --- a/subprojects/MLPCpp +++ b/subprojects/MLPCpp @@ -1 +1 @@ -Subproject commit 665c45b7d3533c977eb1f637918d5b8b75c07d3b +Subproject commit c19c53ea2b85ccfb185f1c6c87044dc0b5bc7ae0