diff --git a/.gitmodules b/.gitmodules index 47e94c29240..061da74abc3 100644 --- a/.gitmodules +++ b/.gitmodules @@ -18,6 +18,9 @@ [submodule "subprojects/CoolProp"] path = subprojects/CoolProp url = https://github.com/CoolProp/CoolProp.git +[submodule "subprojects/cantera"] + path = subprojects/cantera + url = https://github.com/cantera/cantera.git [submodule "externals/opdi"] path = externals/opdi url = https://github.com/SciCompKL/OpDiLib diff --git a/Common/include/CConfig.hpp b/Common/include/CConfig.hpp index 8e0cc89cb37..c335d524a03 100644 --- a/Common/include/CConfig.hpp +++ b/Common/include/CConfig.hpp @@ -81,7 +81,14 @@ class CConfig { su2double EA_ScaleFactor; /*!< \brief Equivalent Area scaling factor */ su2double AdjointLimit; /*!< \brief Adjoint variable limit */ string* ConvField; /*!< \brief Field used for convergence check.*/ - string FluidName; /*!< \brief name of the applied fluid. */ + string FluidName; /*!< \brief name of the applied fluid. */ + string TransportModel; /*!< \brief name transport model used in Cantera. */ + string* GasCompositionNames; /*!< \brief gas composition used in Cantera. */ + string ChemicalMechanismFile; /*!< \brief Chemical Reaction mechanism used in Cantera. */ + string PhaseName; /*!< \brief Name of the phase in the chemical mechanism file used in Cantera. */ + unsigned short n_GasCompositionNames; /*!<\brief Number of gases in mixture composition for Cantera. */ + bool Combustion; /*!< \brief Flag for Combustion Detailed chemistry problems using Cantera. */ + su2double Spark_Temperature; /*!< \bried Spark temperature used for ignition in detailed chemistry using Cantera. */ string* WndConvField; /*!< \brief Function where to apply the windowed convergence criteria for the time average of the unsteady (single zone) flow problem. */ unsigned short nConvField; /*!< \brief Number of fields used to monitor convergence.*/ @@ -4055,6 +4062,47 @@ class CConfig { */ string GetFluid_Name(void) const { return FluidName; } + /*! + * \brief Returns the transport model used in Cantera. + */ + string GetTransport_Model(void) const { return TransportModel; } + + /*! + * \brief Returns the chemical reaction mechanism (mechanism.yaml) used in Cantera. + */ + string GetChemical_MechanismFile(void) const { return ChemicalMechanismFile; } + + /*! + * \brief Returns the name of the pase in the chemical reaction mechanism file used in Cantera. + */ + string GetPhase_Name(void) const { return PhaseName; } + + /*! + * \brief Returns the gas composition used in Cantera. + */ + string GetChemical_GasComposition(unsigned short val_index = 0) const { return GasCompositionNames[val_index]; } + + /*! + * \brief Set the gas composition used in Cantera. + */ + void SetChemical_GasComposition(unsigned short val_index, string gas_composition) const { + GasCompositionNames[val_index] = gas_composition; + } + + /*! + * \brief Get information about the Combustion-Detailed chemistry using Cantera. + * \return TRUE if combustion-detailed chemistry using Cantera is used; otherwise FALSE. + */ + bool GetCombustion(void) const { return Combustion; } + + /*! + * \brief Get High temperature applied during spark ignition. + * \return Spark Temperature. + */ + const su2double GetSpark_Temperature(void) const { + return Spark_Temperature; + } + /*! * \brief Option to define the density model for incompressible flows. * \return Density model option diff --git a/Common/include/option_structure.hpp b/Common/include/option_structure.hpp index 6d97299d09d..1be7f67b808 100644 --- a/Common/include/option_structure.hpp +++ b/Common/include/option_structure.hpp @@ -558,6 +558,7 @@ enum ENUM_FLUIDMODEL { COOLPROP = 10, /*!< \brief Thermodynamics library. */ FLUID_FLAMELET = 11, /*!< \brief lookup table (LUT) method for premixed flamelets. */ DATADRIVEN_FLUID = 12, /*!< \brief multi-layer perceptron driven fluid model. */ + FLUID_CANTERA = 13, /*!< \brief Reacting flows model. */ }; static const MapType FluidModel_Map = { MakePair("STANDARD_AIR", STANDARD_AIR) @@ -573,6 +574,7 @@ static const MapType FluidModel_Map = { MakePair("COOLPROP", COOLPROP) MakePair("DATADRIVEN_FLUID", DATADRIVEN_FLUID) MakePair("FLUID_FLAMELET", FLUID_FLAMELET) + MakePair("FLUID_CANTERA", FLUID_CANTERA) }; /*! @@ -679,6 +681,7 @@ enum class VISCOSITYMODEL { POLYNOMIAL, /*!< \brief Polynomial viscosity. */ FLAMELET, /*!< \brief LUT method for flamelets */ COOLPROP, /*!< \brief CoolProp viscosity. */ + CANTERA, /*!< \brief Cantera viscosity. */ }; static const MapType ViscosityModel_Map = { MakePair("CONSTANT_VISCOSITY", VISCOSITYMODEL::CONSTANT) @@ -686,6 +689,7 @@ static const MapType ViscosityModel_Map = { MakePair("POLYNOMIAL_VISCOSITY", VISCOSITYMODEL::POLYNOMIAL) MakePair("FLAMELET", VISCOSITYMODEL::FLAMELET) MakePair("COOLPROP", VISCOSITYMODEL::COOLPROP) + MakePair("CANTERA", VISCOSITYMODEL::CANTERA) }; /*! @@ -709,6 +713,7 @@ enum class CONDUCTIVITYMODEL { POLYNOMIAL, /*!< \brief Polynomial thermal conductivity. */ FLAMELET, /*!< \brief LUT method for flamelets */ COOLPROP, /*!< \brief COOLPROP thermal conductivity. */ + CANTERA, /*!< \brief CANTERA thermal conductivity. */ }; static const MapType ConductivityModel_Map = { MakePair("CONSTANT_CONDUCTIVITY", CONDUCTIVITYMODEL::CONSTANT) @@ -716,6 +721,7 @@ static const MapType ConductivityModel_Map = { MakePair("POLYNOMIAL_CONDUCTIVITY", CONDUCTIVITYMODEL::POLYNOMIAL) MakePair("FLAMELET", CONDUCTIVITYMODEL::FLAMELET) MakePair("COOLPROP", CONDUCTIVITYMODEL::COOLPROP) + MakePair("CANTERA", CONDUCTIVITYMODEL::CANTERA) }; /*! @@ -739,6 +745,7 @@ enum class DIFFUSIVITYMODEL { UNITY_LEWIS, /*!< \brief Unity Lewis model for mass diffusion in scalar transport. */ CONSTANT_LEWIS, /*!< \brief Different Lewis number model for mass diffusion in scalar transport. */ FLAMELET, /*!< \brief flamelet model for tabulated chemistry, diffusivity from lookup table */ + CANTERA, /*!< \brief Mixture average diffusivity from CANTERA */ }; static const MapType Diffusivity_Model_Map = { @@ -747,6 +754,7 @@ static const MapType Diffusivity_Model_Map = { MakePair("UNITY_LEWIS", DIFFUSIVITYMODEL::UNITY_LEWIS) MakePair("CONSTANT_LEWIS", DIFFUSIVITYMODEL::CONSTANT_LEWIS) MakePair("FLAMELET", DIFFUSIVITYMODEL::FLAMELET) + MakePair("CANTERA", DIFFUSIVITYMODEL::CANTERA) }; /*! diff --git a/Common/src/CConfig.cpp b/Common/src/CConfig.cpp index f272a84ab6f..347ce159ec3 100644 --- a/Common/src/CConfig.cpp +++ b/Common/src/CConfig.cpp @@ -1277,6 +1277,20 @@ void CConfig::SetConfig_Options() { /*!\brief FLUID_NAME \n DESCRIPTION: Fluid name \n OPTIONS: see coolprop homepage \n DEFAULT: nitrogen \ingroup Config*/ addStringOption("FLUID_NAME", FluidName, string("nitrogen")); + /*!\par CONFIG_CATEGORY: Cantera fluid model \ingroup Config*/ + /*!\brief TRANSPORT_MODEL \n DESCRIPTION: Transport model \n OPTIONS: see Cantera homepage \n DEFAULT: mixture-averaged \ingroup Config*/ + addStringOption("TRANSPORT_MODEL", TransportModel, string("mixture-averaged")); + /*!\brief CHEMICAL_MECHANISM_FILE \n DESCRIPTION: Chemical reaction mechanism \n OPTIONS: see Cantera homepage \n DEFAULT: h2o2.yaml \ingroup Config*/ + addStringOption("CHEMICAL_MECHANISM_FILE", ChemicalMechanismFile, string("h2o2.yaml")); + /*!\brief PHASE_NAME \n DESCRIPTION: name of the phase in the chemical mechanism file \n OPTIONS: see Cantera homepage \n DEFAULT: gri30 \ingroup Config*/ + addStringOption("PHASE_NAME", PhaseName, string("ohmech")); + /*!\brief GAS_COMPOSITION_NAMES \n DESCRIPTION: Gas composition names \n OPTIONS: see Cantera homepage \n DEFAULT: \ingroup Config*/ + addStringListOption("GAS_COMPOSITION_NAMES", n_GasCompositionNames, GasCompositionNames); + /*\brief COMBUSTION \n DESCRIPTION: Combustion Detailed chemistry using Cantera \n DEFAULT: false \ingroup Config*/ + addBoolOption("COMBUSTION", Combustion, false); + /*!\brief SPARK_TEMPERATURE \n DESCRIPTION: Spark temperature used for ignition in detailed chemistry using Cantera \n DEFAULT: 1000 K \ingroup Config*/ + addDoubleOption("SPARK_TEMPERATURE", Spark_Temperature, 1000.0); + /*!\par CONFIG_CATEGORY: Data-driven fluid model parameters \ingroup Config*/ /*!\brief INTERPOLATION_METHOD \n DESCRIPTION: Interpolation method used to determine the thermodynamic state of the fluid. \n OPTIONS: See \link DataDrivenMethod_Map \endlink DEFAULT: MLP \ingroup Config*/ addEnumOption("INTERPOLATION_METHOD",datadriven_ParsedOptions.interp_algorithm_type, DataDrivenMethod_Map, ENUM_DATADRIVEN_METHOD::LUT); @@ -3621,7 +3635,8 @@ void CConfig::SetPostprocessing(SU2_COMPONENT val_software, unsigned short val_i (Kind_FluidModel == FLUID_MIXTURE) || (Kind_FluidModel == FLUID_FLAMELET) || (Kind_FluidModel == INC_IDEAL_GAS_POLY) || - (Kind_FluidModel == CONSTANT_DENSITY)); + (Kind_FluidModel == CONSTANT_DENSITY)|| + (Kind_FluidModel == FLUID_CANTERA)); bool noneq_gas = ((Kind_FluidModel == MUTATIONPP) || (Kind_FluidModel == SU2_NONEQ)); bool standard_air = ((Kind_FluidModel == STANDARD_AIR)); @@ -4297,6 +4312,67 @@ void CConfig::SetPostprocessing(SU2_COMPONENT val_software, unsigned short val_i } } + if ((Kind_FluidModel != FLUID_CANTERA) && (Combustion == true)) { + SU2_MPI::Error( + "The use of COMBUSTION=YES requires the use of FLUID_MIXTURE=FLUID_CANTERA,\n" + "detailed chemistry cannot be performed with other fluid models", + CURRENT_FUNCTION); + } + + if ((flamelet_ParsedOptions.ignition_method != FLAMELET_INIT_TYPE::SPARK) && (Combustion == true)) { + SU2_MPI::Error( + "The use of COMBUSTION=YES requires the use of FLAME_INIT_METHOD=SPARK,\n" + "Other ignition methods are not currently available", + CURRENT_FUNCTION); + } + + if (Kind_FluidModel == FLUID_CANTERA) { + /*--- Check whether the number of entries of the GAS_COMPOSITION_NAMES equals the number of transported scalar + equations solved + 1.--- */ + if (n_GasCompositionNames != nSpecies_Init + 1) { + SU2_MPI::Error( + "The use of FLUID_CANTERA requires the number of entries for GAS_COMPOSITION_NAMES,\n" + "to be equal to the number of entries of SPECIES_INIT + 1", + CURRENT_FUNCTION); + } + /*--- Check whether the density model used is correct, in the case of FLUID_CANTERA the density model must be + VARIABLE. Otherwise, if the density model is CONSTANT, the scalars will not have influence the mixture density + and it will remain constant through the complete domain. --- */ + if (Kind_DensityModel != INC_DENSITYMODEL::VARIABLE) { + SU2_MPI::Error("The use of FLUID_CANTERA requires the INC_DENSITY_MODEL option to be VARIABLE", + CURRENT_FUNCTION); + } + /*--- Check whether the Kind scalar model used is correct, in the case of FLUID_CANTERA the kind scalar model must + be SPECIES_TRANSPORT. Otherwise, if the scalar model is NONE, the species transport equations will not be solved. + --- */ + if (Kind_Species_Model != SPECIES_MODEL::SPECIES_TRANSPORT) { + SU2_MPI::Error("The use of FLUID_CANTERA requires the KIND_SCALAR_MODEL option to be SPECIES_TRANSPORT", + CURRENT_FUNCTION); + } + + if (Ref_Inc_NonDim != DIMENSIONAL) { + SU2_MPI::Error( + "The use of FLUID_CANTERA requiere the option INC_NONDIM= DIMENSIONAL, the nondimensionalization is " + "currently unavailable.", + CURRENT_FUNCTION); + } + + if (Kind_ConductivityModel != CONDUCTIVITYMODEL::CANTERA) { + SU2_MPI::Error("The use of FLUID_CANTERA requires the CONDUCTIVITY_MODEL option to be CANTERA", + CURRENT_FUNCTION); + } + + if (Kind_ViscosityModel != VISCOSITYMODEL::CANTERA) { + SU2_MPI::Error("The use of FLUID_CANTERA requires the VISCOSITY_MODEL option to be CANTERA", + CURRENT_FUNCTION); + } + + if (Kind_Diffusivity_Model != DIFFUSIVITYMODEL::CANTERA) { + SU2_MPI::Error("The use of FLUID_CANTERA requires the DIFFUSIVITY_MODEL option to be CANTERA", + CURRENT_FUNCTION); + } + } + if (Kind_Species_Model == SPECIES_MODEL::FLAMELET) { @@ -5279,7 +5355,7 @@ void CConfig::SetPostprocessing(SU2_COMPONENT val_software, unsigned short val_i } if (Kind_DensityModel == INC_DENSITYMODEL::VARIABLE) { - if (Kind_FluidModel != INC_IDEAL_GAS && Kind_FluidModel != INC_IDEAL_GAS_POLY && Kind_FluidModel != FLUID_MIXTURE && Kind_FluidModel != FLUID_FLAMELET) { + if (Kind_FluidModel != INC_IDEAL_GAS && Kind_FluidModel != INC_IDEAL_GAS_POLY && Kind_FluidModel != FLUID_MIXTURE && Kind_FluidModel != FLUID_FLAMELET && Kind_FluidModel != FLUID_CANTERA) { SU2_MPI::Error("Variable density incompressible solver limited to ideal gases.\n Check the fluid model options (use INC_IDEAL_GAS, INC_IDEAL_GAS_POLY).", CURRENT_FUNCTION); } } diff --git a/SU2_CFD/include/fluid/CFluidCantera.hpp b/SU2_CFD/include/fluid/CFluidCantera.hpp new file mode 100644 index 00000000000..95a3492862e --- /dev/null +++ b/SU2_CFD/include/fluid/CFluidCantera.hpp @@ -0,0 +1,156 @@ +/*! + * \file CFluidCantera.hpp + * \brief Defines the multicomponent incompressible Ideal Gas model for reacting flows. + * \author T. Economon, Cristopher Morales Ubal + * \version 8.4.0 "Harrier" + * + * SU2 Project Website: https://su2code.github.io + * + * The SU2 Project is maintained by the SU2 Foundation + * (http://su2foundation.org) + * + * Copyright 2012-2024, SU2 Contributors (cf. AUTHORS.md) + * + * SU2 is free software; you can redistribute it and/or + * modify it under the terms of the GNU Lesser General Public + * License as published by the Free Software Foundation; either + * version 2.1 of the License, or (at your option) any later version. + * + * SU2 is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + * Lesser General Public License for more details. + * + * You should have received a copy of the GNU Lesser General Public + * License along with SU2. If not, see . + */ + +#pragma once + +#include +#include + +#include "CFluidModel.hpp" + +#if defined(HAVE_CANTERA) +#define USE_CANTERA +namespace Cantera { +class Solution; +} +#endif + + +/*! + * \class CFluidCantera + * \brief Child class for defining reacting incompressible ideal gas mixture model. + * \author: T. Economon, Cristopher Morales Ubal + */ +class CFluidCantera final : public CFluidModel { + private: + #ifdef USE_CANTERA + const int n_species_mixture; /*!< \brief Number of species in mixture. */ + const su2double Pressure_Thermodynamic; /*!< \brief Constant pressure thermodynamic. */ + const su2double GasConstant_Ref; /*!< \brief Gas constant reference needed for Nondimensional problems. */ + const su2double Prandtl_Turb_Number; /*!< \brief Prandlt turbulent number.*/ + const su2double Schmidt_Turb_Number; /*!< \brief Schmidt turbulent number.*/ + const string Transport_Model; /*!< \brief Transport model used for computing mixture properties*/ + const string Chemical_MechanismFile; /*!< \brief Chemical reaction mechanism used for in cantera*/ + const string Phase_Name; /*!< \brief Name of the phase used for in cantera*/ + const bool Combustion; /*!< \brief Flag for problems involving combustion.*/ + + static constexpr int ARRAYSIZE = 16; + + su2double Heat_Release; /*!< \brief heat release due to combustion */ + std::array speciesIndices; /*!< \brief Species indices within Cantera library. */ + std::shared_ptr sol; /*!< \brief Object needed to describe a chemically-reacting solution*/ + std::array chemicalSourceTerm; /*!< \brief chemical source term of all species*/ + std::array molarMasses; /*!< \brief Molar masses of all species. */ + std::array enthalpyFormation; /*!< \brief Enthalpy of Formation of all species. */ + std::array massFractions; /*!< \brief Mass fractions of all species. */ + std::array massDiffusivity; /*!< \brief mass diffusivity of all species. */ + vector netProductionRates; /*!< \brief Net chemical production rates of all species*/ + mutable vector enthalpiesSpecies; /*!< \brief Molar enthalpies of all species. */ + mutable vector specificHeatSpecies; /*!< \brief Molar heat capacities of all species. */ + + /*! + * \brief Compute heat release due to combustion. + */ + void ComputeHeatRelease(); + + /*! + * \brief Set enthalpies of formation. + */ + void SetEnthalpyFormation(const CConfig* config); + + /*! + * \brief Set species mass fraction array for Cantera. + * \param[in] val_scalars - Scalar mass fractions. + */ + void SetMassFractions(const su2double* val_scalars); +#endif + + public: + /*! + * \brief Constructor of the class. + */ + CFluidCantera(su2double val_operating_pressure, const CConfig* config); + + #ifdef USE_CANTERA + + /*! + * \brief Get fluid laminar viscosity. + */ + inline su2double GetLaminarViscosity() override { return Mu; } + + /*! + * \brief Get fluid thermal conductivity. + */ + inline su2double GetThermalConductivity() override { return Kt + Mu_Turb * Cp / Prandtl_Turb_Number; } + + /*! + * \brief Get fluid mass diffusivity. + * \param[in] ivar - index of species. + */ + inline su2double GetMassDiffusivity(int ivar) override { return massDiffusivity[speciesIndices[ivar]]; } + + /*! + * \brief Compute chemical source term for species. + */ + void ComputeChemicalSourceTerm() override; + + /*! + * \brief Get Chemical source term species. + * \param[in] ivar - index of species. + */ + inline su2double GetChemicalSourceTerm(int ivar) override { return chemicalSourceTerm[ivar]; } + + /*! + * \brief Get Heat release due to combustion. + */ + inline su2double GetHeatRelease() override { return Heat_Release; } + + /*! + * \brief Get enthalpy diffusivity terms. + */ + void GetEnthalpyDiffusivity(su2double* enthalpy_diffusions) const override; + + /*! + * \brief Get gradient enthalpy diffusivity terms. + */ + void GetGradEnthalpyDiffusivity(su2double* grad_enthalpy_diffusions) const override; + + /*! + * \brief Set the Dimensionless State using Temperature. + * \param[in] val_temperature - Temperature value at the point. + * \param[in] val_scalars - Scalar mass fractions. + */ + void SetTDState_T(su2double val_temperature, const su2double* val_scalars) override; + + /*! + * \brief Virtual member. + * \param[in] val_enthalpy - Enthalpy value at the point. + * \param[in] val_scalars - Scalar mass fractions. + */ + void SetTDState_h(su2double val_enthalpy, const su2double* val_scalars) override; + #endif +}; \ No newline at end of file diff --git a/SU2_CFD/include/fluid/CFluidModel.hpp b/SU2_CFD/include/fluid/CFluidModel.hpp index 4072a82c907..13edf2115fb 100644 --- a/SU2_CFD/include/fluid/CFluidModel.hpp +++ b/SU2_CFD/include/fluid/CFluidModel.hpp @@ -192,6 +192,22 @@ class CFluidModel { return mass_diffusivity; } + /*! + * \brief Compute chemical source term for species. + */ + virtual void ComputeChemicalSourceTerm() {}; + + /*! + * \brief Get Chemical source term species. + * \param[in] iVar - index of species. + */ + inline virtual su2double GetChemicalSourceTerm(int iVar) { return 0.0; } + + /*! + * \brief Get Heat release due to combustion. + */ + inline virtual su2double GetHeatRelease() { return 0.0; } + /*! * \brief Get the enthalpy diffusivity terms for all species being solved. * diff --git a/SU2_CFD/include/numerics/CNumerics.hpp b/SU2_CFD/include/numerics/CNumerics.hpp index ec511a3fa1d..009f9a6c00c 100644 --- a/SU2_CFD/include/numerics/CNumerics.hpp +++ b/SU2_CFD/include/numerics/CNumerics.hpp @@ -62,7 +62,9 @@ class CNumerics { const su2double delta [3][3] = {{1.0, 0.0, 0.0},{0.0, 1.0, 0.0}, {0.0, 0.0, 1.0}}; /*!< \brief Identity matrix. */ const su2double *Diffusion_Coeff_i, /*!< \brief Species diffusion coefficients at point i. */ - *Diffusion_Coeff_j; /*!< \brief Species diffusion coefficients at point j. */ + *Diffusion_Coeff_j, /*!< \brief Species diffusion coefficients at point j. */ + *Chemical_Source_Term_i, /*!< \brief Species diffusion coefficients at point i. */ + *Chemical_Source_Term_j; /*!< \brief Species diffusion coefficients at point j. */ su2double Laminar_Viscosity_i, /*!< \brief Laminar viscosity at point i. */ Laminar_Viscosity_j; /*!< \brief Laminar viscosity at point j. */ @@ -830,6 +832,16 @@ class CNumerics { Diffusion_Coeff_j = val_diffusioncoeff_j; } + /*! + * \brief Set the Chemical source term + * \param[in] val_source_term_i - Value of the chemical source term at i. + * \param[in] val_source_term_j - Value of the chemical source term at j + */ + inline void SetChemicalSourceTerm(const su2double* val_source_term_i, const su2double* val_source_term_j) { + Chemical_Source_Term_i = val_source_term_i; + Chemical_Source_Term_j = val_source_term_j; + } + /*! * \brief Set the heat flux due to enthalpy diffusion * \param[in] val_heatfluxdiffusion - Value of the heat flux due to enthalpy diffusion. diff --git a/SU2_CFD/include/numerics/scalar/scalar_convection.hpp b/SU2_CFD/include/numerics/scalar/scalar_convection.hpp index e40749a5c1d..3b97b4c6088 100644 --- a/SU2_CFD/include/numerics/scalar/scalar_convection.hpp +++ b/SU2_CFD/include/numerics/scalar/scalar_convection.hpp @@ -47,7 +47,7 @@ template class CUpwScalar : public CNumerics { protected: - enum : unsigned short {MAXNVAR = 8}; + enum : unsigned short {MAXNVAR = 20}; const FlowIndices idx; /*!< \brief Object to manage the access to the flow primitives. */ su2double a0 = 0.0; /*!< \brief The maximum of the face-normal velocity and 0. */ diff --git a/SU2_CFD/include/numerics/scalar/scalar_diffusion.hpp b/SU2_CFD/include/numerics/scalar/scalar_diffusion.hpp index 580c0bedcef..fa5186e9769 100644 --- a/SU2_CFD/include/numerics/scalar/scalar_diffusion.hpp +++ b/SU2_CFD/include/numerics/scalar/scalar_diffusion.hpp @@ -59,7 +59,7 @@ struct CNoFlowIndices { template class CAvgGrad_Scalar : public CNumerics { protected: - enum : unsigned short {MAXNVAR = 8}; + enum : unsigned short {MAXNVAR = 20}; const FlowIndices idx; /*!< \brief Object to manage the access to the flow primitives. */ su2double Proj_Mean_GradScalarVar[MAXNVAR]; /*!< \brief Mean_gradScalarVar DOT normal, corrected if required. */ diff --git a/SU2_CFD/include/numerics/species/species_sources.hpp b/SU2_CFD/include/numerics/species/species_sources.hpp index a33a0585948..8eb75eddd95 100644 --- a/SU2_CFD/include/numerics/species/species_sources.hpp +++ b/SU2_CFD/include/numerics/species/species_sources.hpp @@ -91,3 +91,30 @@ class CSourceAxisymmetric_Species : public CSourceBase_Species { ResidualType<> ComputeResidual(const CConfig* config) override; }; + +/*! + * \class CSourceCombustion_Species + * \brief Class for source term for combustion problems. + * \ingroup SourceDiscr + * \author C.Morales Ubal + */ +template +class CSourceCombustion_Species : public CSourceBase_Species { + + public: + /*! + * \brief Constructor of the class. + * \param[in] val_nDim - Number of dimensions of the problem. + * \param[in] val_nVar - Number of variables of the problem. + * \param[in] config - Definition of the particular problem. + */ + CSourceCombustion_Species(unsigned short val_nDim, unsigned short val_nVar, const CConfig* config); + + /*! + * \brief Residual of the combustion source term. + * \param[in] config - Definition of the particular problem. + * \return Lightweight const-view of residual and Jacobian. + */ + ResidualType<> ComputeResidual(const CConfig* config) override; + +}; diff --git a/SU2_CFD/include/solvers/CIncNSSolver.hpp b/SU2_CFD/include/solvers/CIncNSSolver.hpp index 712855d0a20..b9dc5b18686 100644 --- a/SU2_CFD/include/solvers/CIncNSSolver.hpp +++ b/SU2_CFD/include/solvers/CIncNSSolver.hpp @@ -36,6 +36,8 @@ * \author F. Palacios, T. Economon, T. Albring */ class CIncNSSolver final : public CIncEulerSolver { + FluidFlamelet_ParsedOptions flamelet_config_options; + su2double TemperatureLimits[2]; /*! * \brief Generic implementation of the isothermal, heatflux and heat-transfer/convection walls. diff --git a/SU2_CFD/include/solvers/CSpeciesFlameletSolver.hpp b/SU2_CFD/include/solvers/CSpeciesFlameletSolver.hpp index 6b9c4517298..617d9f10814 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: - FluidFlamelet_ParsedOptions flamelet_config_options; bool include_mixture_fraction = false; /*!< \brief include mixture fraction as a controlling variable. */ /*! * \brief Compute the preconditioner for low-Mach flows. diff --git a/SU2_CFD/include/solvers/CSpeciesSolver.hpp b/SU2_CFD/include/solvers/CSpeciesSolver.hpp index 3a657d97c01..a68363a0e34 100644 --- a/SU2_CFD/include/solvers/CSpeciesSolver.hpp +++ b/SU2_CFD/include/solvers/CSpeciesSolver.hpp @@ -40,6 +40,7 @@ class CSpeciesSolver : public CScalarSolver { protected: unsigned short Inlet_Position; /*!< \brief Column index for scalar variables in inlet files. */ vector Inlet_SpeciesVars; /*!< \brief Species variables at inlet profiles. */ + FluidFlamelet_ParsedOptions flamelet_config_options; vector Wall_SpeciesVars; /*!< \brief Species variables at profiles. */ vector> CustomBoundaryScalar; diff --git a/SU2_CFD/include/variables/CSpeciesVariable.hpp b/SU2_CFD/include/variables/CSpeciesVariable.hpp index 76180d0c57b..b14f3983293 100644 --- a/SU2_CFD/include/variables/CSpeciesVariable.hpp +++ b/SU2_CFD/include/variables/CSpeciesVariable.hpp @@ -36,6 +36,8 @@ class CSpeciesVariable : public CScalarVariable { protected: MatrixType Diffusivity; /*!< \brief Matrix (nPoint,nVar) of mass diffusivities for scalar transport. */ + MatrixType SpeciesSourceTerm; /*!< \brief Matrix (nPoint, nVar) of chemical source terms for species transport*/ + VectorType HeatRelease; /*!< \brief Vector of heat release due to combustion for species transport*/ public: static constexpr size_t MAXNVAR = 20; /*!< \brief Max number of variables for static arrays. Increase, if necessary. */ @@ -74,4 +76,41 @@ class CSpeciesVariable : public CScalarVariable { * \return Pointer to the mass diffusivities */ inline const su2double* GetDiffusivity(unsigned long iPoint) const { return Diffusivity[iPoint]; } + + /*! + * \brief Set the value of the chemical source term for species transport + * \param[in] val_sourceTerm - the chemical source term. + * \param[in] val_ivar - eqn. index to the chemical source term. + */ + inline void SetChemicalSourceTerm(unsigned long iPoint, su2double val_sourceTerm, unsigned short val_ivar) { + SpeciesSourceTerm(iPoint, val_ivar) = val_sourceTerm; + } + + /*! + * \brief Get the value of the chemical source term for species transport + * \param[in] val_ivar - eqn. index to the chemical source term. + * \return Value of the chemical source term. + */ + inline su2double GetChemicalSourceTerm(unsigned long iPoint, unsigned short val_ivar) const { + return SpeciesSourceTerm(iPoint, val_ivar); + } + + /*! + * \brief Get the value of the chemical source term for species transport + * \return Pointer to the chemical source term + */ + inline const su2double* GetChemicalSourceTerm(unsigned long iPoint) const { return SpeciesSourceTerm[iPoint]; } + + /*! + * \brief Get heat release due to combustion + * \param[in] iPoint - Point index. + * \return Value of the heat release due to combustion. + */ + inline su2double GetHeatRelease(unsigned long iPoint) const { return HeatRelease(iPoint); } + + /*! + * \brief Set heat release due to combustion. + * \param[in] iPoint - Point index. + */ + inline void SetHeatRelease(unsigned long iPoint, su2double val_heatRelease) { HeatRelease(iPoint) = val_heatRelease; } }; diff --git a/SU2_CFD/include/variables/CVariable.hpp b/SU2_CFD/include/variables/CVariable.hpp index 528adc139db..2a5d9c8e51f 100644 --- a/SU2_CFD/include/variables/CVariable.hpp +++ b/SU2_CFD/include/variables/CVariable.hpp @@ -1165,6 +1165,20 @@ class CVariable { */ inline virtual su2double GetDiffusivity(unsigned long iPoint, unsigned short val_ivar) const { return 0.0; } + /*! + * \brief A virtual member. + * \param[in] iPoint - Point index. + * \return Value of the chemical source term. + */ + inline virtual su2double GetChemicalSourceTerm(unsigned long iPoint, unsigned short val_ivar) const { return 0.0; } + + /*! + * \brief A virtual member. + * \param[in] iPoint - Point index. + * \return Value of the heat release due to combustion. + */ + inline virtual su2double GetHeatRelease(unsigned long iPoint) const { return 0.0; } + /*! * \brief A virtual member. * \param[in] iPoint - Point index. diff --git a/SU2_CFD/src/drivers/CDriver.cpp b/SU2_CFD/src/drivers/CDriver.cpp index 5904ed9febb..08d3c8242f4 100644 --- a/SU2_CFD/src/drivers/CDriver.cpp +++ b/SU2_CFD/src/drivers/CDriver.cpp @@ -1416,7 +1416,11 @@ void CDriver::InstantiateSpeciesNumerics(unsigned short nVar_Species, int offset else { numerics[iMGlevel][SPECIES_SOL][source_first_term] = new CSourceNothing(nDim, nVar_Species, config); } - numerics[iMGlevel][SPECIES_SOL][source_second_term] = new CSourceNothing(nDim, nVar_Species, config); + if (config->GetCombustion() == YES) { + numerics[iMGlevel][SPECIES_SOL][source_second_term] = new CSourceCombustion_Species(nDim, nVar_Species, config); + } else { + numerics[iMGlevel][SPECIES_SOL][source_second_term] = new CSourceNothing(nDim, nVar_Species, config); + } } } diff --git a/SU2_CFD/src/fluid/CDataDrivenFluid.cpp b/SU2_CFD/src/fluid/CDataDrivenFluid.cpp index 1de3c2ef1fc..8bb25a55293 100644 --- a/SU2_CFD/src/fluid/CDataDrivenFluid.cpp +++ b/SU2_CFD/src/fluid/CDataDrivenFluid.cpp @@ -471,4 +471,4 @@ size_t MiniTable2D::FindNode(const size_t iX, const su2double val_x, const size_ } } return iMin; -} \ No newline at end of file +} diff --git a/SU2_CFD/src/fluid/CFluidCantera.cpp b/SU2_CFD/src/fluid/CFluidCantera.cpp new file mode 100644 index 00000000000..43a4f6df039 --- /dev/null +++ b/SU2_CFD/src/fluid/CFluidCantera.cpp @@ -0,0 +1,194 @@ +/*! + * \file CFluidCantera.cpp + * \brief Defines the multicomponent incompressible Ideal Gas model for reacting flows. + * \author T. Economon, Cristopher Morales Ubal + * \version 8.4.0 "Harrier" + * + * SU2 Project Website: https://su2code.github.io + * + * The SU2 Project is maintained by the SU2 Foundation + * (http://su2foundation.org) + * + * Copyright 2012-2024, SU2 Contributors (cf. AUTHORS.md) + * + * SU2 is free software; you can redistribute it and/or + * modify it under the terms of the GNU Lesser General Public + * License as published by the Free Software Foundation; either + * version 2.1 of the License, or (at your option) any later version. + * + * SU2 is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + * Lesser General Public License for more details. + * + * You should have received a copy of the GNU Lesser General Public + * License along with SU2. If not, see . + */ + +#include "../../include/fluid/CFluidCantera.hpp" + +#include + +// #include +#ifdef USE_CANTERA +#include +#include + +using namespace Cantera; + +CFluidCantera::CFluidCantera(su2double value_pressure_operating, const CConfig* config) + : CFluidModel(), + n_species_mixture(config->GetnSpecies() + 1), + Pressure_Thermodynamic(value_pressure_operating), + GasConstant_Ref(config->GetGas_Constant_Ref()), + Prandtl_Turb_Number(config->GetPrandtl_Turb()), + Schmidt_Turb_Number(config->GetSchmidt_Number_Turbulent()), + Transport_Model(config->GetTransport_Model()), + Chemical_MechanismFile(config->GetChemical_MechanismFile()), + Phase_Name(config->GetPhase_Name()), + Combustion(config->GetCombustion()) { + if (n_species_mixture > ARRAYSIZE) { + SU2_MPI::Error("Too many species, increase ARRAYSIZE", CURRENT_FUNCTION); + } + sol = std::shared_ptr(newSolution(Chemical_MechanismFile, Phase_Name, Transport_Model)); + sol->thermo()->getMolecularWeights(&molarMasses[0]); + for (int iVar = 0; iVar < n_species_mixture; iVar++) { + speciesIndices[iVar] = sol->thermo()->speciesIndex(config->GetChemical_GasComposition(iVar)); + } + enthalpiesSpecies.resize(sol->thermo()->nSpecies()); + specificHeatSpecies.resize(sol->thermo()->nSpecies()); + netProductionRates.resize(sol->thermo()->nSpecies()); + if (Combustion) SetEnthalpyFormation(config); +} + +void CFluidCantera::SetEnthalpyFormation(const CConfig* config) { + SetMassFractions(config->GetSpecies_Init()); + sol->thermo()->setMassFractions(massFractions.data()); + sol->thermo()->setState_TP(STD_REF_TEMP, Pressure_Thermodynamic); + sol->thermo()->getEnthalpy_RT_ref(enthalpiesSpecies.data()); + for (int iVar = 0; iVar < n_species_mixture; iVar++) { + enthalpyFormation[iVar] = + GasConstant * STD_REF_TEMP * enthalpiesSpecies[speciesIndices[iVar]] / molarMasses[speciesIndices[iVar]]; + } +} + +void CFluidCantera::ComputeChemicalSourceTerm() { + sol->kinetics()->getNetProductionRates(netProductionRates.data()); + for (int iVar = 0; iVar < n_species_mixture - 1.0; iVar++) { + chemicalSourceTerm[iVar] = molarMasses[speciesIndices[iVar]] * netProductionRates[speciesIndices[iVar]]; + } +} + +void CFluidCantera::ComputeHeatRelease() { + sol->kinetics()->getNetProductionRates(netProductionRates.data()); + Heat_Release = 0.0; + for (int iVar = 0; iVar < n_species_mixture; iVar++) { + Heat_Release += + -1.0 * enthalpyFormation[iVar] * molarMasses[speciesIndices[iVar]] * netProductionRates[speciesIndices[iVar]]; + } +} + +void CFluidCantera::GetEnthalpyDiffusivity(su2double* enthalpy_diffusions) const { + sol->thermo()->getEnthalpy_RT_ref(enthalpiesSpecies.data()); + for (int iVar = 0; iVar < n_species_mixture - 1; iVar++) { + enthalpy_diffusions[iVar] = + Density * GasConstant * Temperature * + ((enthalpiesSpecies[speciesIndices[iVar]] * massDiffusivity[speciesIndices[iVar]] / molarMasses[speciesIndices[iVar]]) - + (enthalpiesSpecies[speciesIndices[n_species_mixture - 1]] * massDiffusivity[speciesIndices[n_species_mixture - 1]] / + molarMasses[speciesIndices[n_species_mixture - 1]])); + enthalpy_diffusions[iVar] += + Mu_Turb * GasConstant * Temperature * + ((enthalpiesSpecies[speciesIndices[iVar]] / molarMasses[speciesIndices[iVar]]) - + (enthalpiesSpecies[speciesIndices[n_species_mixture - 1]] / + molarMasses[speciesIndices[n_species_mixture - 1]])) / Schmidt_Turb_Number; + } +} + +void CFluidCantera::GetGradEnthalpyDiffusivity(su2double* grad_enthalpy_diffusions) const { + sol->thermo()->getCp_R_ref(specificHeatSpecies.data()); + for (int iVar = 0; iVar < n_species_mixture - 1; iVar++) { + grad_enthalpy_diffusions[iVar] = + Density * GasConstant * + ((specificHeatSpecies[speciesIndices[iVar]] * massDiffusivity[speciesIndices[iVar]] / molarMasses[speciesIndices[iVar]]) - + (specificHeatSpecies[speciesIndices[n_species_mixture - 1]] * massDiffusivity[speciesIndices[n_species_mixture - 1]] / + molarMasses[speciesIndices[n_species_mixture - 1]])); + grad_enthalpy_diffusions[iVar] += Mu_Turb * GasConstant * + ((specificHeatSpecies[speciesIndices[iVar]] / molarMasses[speciesIndices[iVar]]) - + (specificHeatSpecies[speciesIndices[n_species_mixture - 1]] / + molarMasses[speciesIndices[n_species_mixture - 1]])) / + Schmidt_Turb_Number; + } +} + +void CFluidCantera::SetMassFractions(const su2double* val_scalars) { + su2double val_scalars_sum{0.0}; + massFractions.fill(0.0); + for (int i_scalar = 0; i_scalar < n_species_mixture - 1; i_scalar++) { + massFractions[speciesIndices[i_scalar]] = val_scalars[i_scalar]; + val_scalars_sum += val_scalars[i_scalar]; + } + massFractions[speciesIndices[n_species_mixture - 1]] = 1.0 - val_scalars_sum; +} + +void CFluidCantera::SetTDState_T(const su2double val_temperature, const su2double* val_scalars) { + Temperature = val_temperature; + SetMassFractions(val_scalars); + + sol->thermo()->setMassFractions(massFractions.data()); + sol->thermo()->setState_TP(Temperature, Pressure_Thermodynamic); + Density = sol->thermo()->density(); + Enthalpy = sol->thermo()->enthalpy_mass(); + Cp = sol->thermo()->cp_mass(); + Cv = sol->thermo()->cv_mass(); + Mu = sol->transport()->viscosity(); + Kt = sol->transport()->thermalConductivity(); + + sol->transport()->getMixDiffCoeffsMass(massDiffusivity.data()); + if (Combustion) ComputeHeatRelease(); +} + +void CFluidCantera::SetTDState_h(const su2double val_enthalpy, const su2double* val_scalars) { + Enthalpy = val_enthalpy; + /*--- convergence criterion for temperature in [K], high accuracy needed for restarts. ---*/ + su2double toll = 1e-5; + su2double temp_iter = 300.0; + su2double delta_temp_iter = 1e10; + su2double delta_enthalpy_iter; + const int counter_limit = 20; + + int counter = 0; + + /*--- Set mass fractions. ---*/ + SetMassFractions(val_scalars); + sol->thermo()->setMassFractions(massFractions.data()); + + /*--- Computing temperature given enthalpy and species mass fractions using Newton-Raphson. ---*/ + while ((abs(delta_temp_iter) > toll) && (counter++ < counter_limit)) { + /*--- Set thermodynamic state based on the current value of temperature. ---*/ + sol->thermo()->setState_TP(temp_iter, Pressure_Thermodynamic); + + su2double Enthalpy_iter = sol->thermo()->enthalpy_mass(); + su2double Cp_iter = sol->thermo()->cp_mass(); + + delta_enthalpy_iter = Enthalpy - Enthalpy_iter; + + delta_temp_iter = delta_enthalpy_iter / Cp_iter; + + temp_iter += delta_temp_iter; + if (temp_iter < 0.0) { + cout << "Warning: Negative temperature has been found during Newton-Raphson" << endl; + break; + } + } + Temperature = temp_iter; + if (counter == counter_limit) { + cout << "Warning Newton-Raphson exceed number of max iteration in temperature computation" << endl; + } + SetTDState_T(Temperature, val_scalars); +} + +#else +CFluidCantera::CFluidCantera(su2double value_pressure_operating, const CConfig* config) { + SU2_MPI::Error("SU2 was not compiled with Cantera(-Denable-cantera=true)", CURRENT_FUNCTION); +} +#endif \ No newline at end of file diff --git a/SU2_CFD/src/fluid/CFluidModel.cpp b/SU2_CFD/src/fluid/CFluidModel.cpp index a5139f3e557..60296644d22 100644 --- a/SU2_CFD/src/fluid/CFluidModel.cpp +++ b/SU2_CFD/src/fluid/CFluidModel.cpp @@ -37,6 +37,7 @@ #include "../../include/fluid/CConstantSchmidt.hpp" #include "../../include/fluid/CConstantViscosity.hpp" #include "../../include/fluid/CFluidScalar.hpp" +#include "../../include/fluid/CFluidCantera.hpp" #include "../../include/fluid/CPolynomialConductivity.hpp" #include "../../include/fluid/CPolynomialConductivityRANS.hpp" #include "../../include/fluid/CPolynomialViscosity.hpp" @@ -62,6 +63,9 @@ unique_ptr CFluidModel::MakeLaminarViscosityModel(const CConfig case VISCOSITYMODEL::FLAMELET: /*--- Viscosity is obtained from the LUT ---*/ return nullptr; + case VISCOSITYMODEL::CANTERA: + /*--- Viscosity is obtained from the Cantera library---*/ + return nullptr; default: SU2_MPI::Error("Viscosity model not available.", CURRENT_FUNCTION); return nullptr; @@ -115,6 +119,9 @@ unique_ptr CFluidModel::MakeThermalConductivityModel(const C case CONDUCTIVITYMODEL::FLAMELET: /*--- Conductivity is obtained from the LUT ---*/ return nullptr; + case CONDUCTIVITYMODEL::CANTERA: + /*--- Conductivity is obtained from Cantera library---*/ + return nullptr; default: SU2_MPI::Error("Conductivity model not available.", CURRENT_FUNCTION); return nullptr; @@ -144,6 +151,10 @@ unique_ptr CFluidModel::MakeMassDiffusivityModel(const CConfi /*--- Diffusivity is obtained from the LUT ---*/ return nullptr; break; + case DIFFUSIVITYMODEL::CANTERA: + /*--- Diffusivity is obtained from CANTERA. Transport model option determines the diffusivity model---*/ + return nullptr; + break; default: SU2_MPI::Error("Diffusivity model not available.", CURRENT_FUNCTION); return nullptr; diff --git a/SU2_CFD/src/interfaces/cht/CConjugateHeatInterface.cpp b/SU2_CFD/src/interfaces/cht/CConjugateHeatInterface.cpp index 3b5309131a0..be33c648eca 100644 --- a/SU2_CFD/src/interfaces/cht/CConjugateHeatInterface.cpp +++ b/SU2_CFD/src/interfaces/cht/CConjugateHeatInterface.cpp @@ -99,6 +99,7 @@ void CConjugateHeatInterface::GetDonor_Variable(CSolver *donor_solution, CGeomet case CONDUCTIVITYMODEL::CONSTANT: case CONDUCTIVITYMODEL::FLAMELET: + case CONDUCTIVITYMODEL::CANTERA: case CONDUCTIVITYMODEL::COOLPROP: thermal_conductivity = thermal_conductivityND*donor_config->GetThermal_Conductivity_Ref(); break; diff --git a/SU2_CFD/src/meson.build b/SU2_CFD/src/meson.build index d4db53843e9..91cb9226b1a 100644 --- a/SU2_CFD/src/meson.build +++ b/SU2_CFD/src/meson.build @@ -13,7 +13,8 @@ su2_cfd_src += files(['fluid/CFluidModel.cpp', 'fluid/CNEMOGas.cpp', 'fluid/CMutationTCLib.cpp', 'fluid/CSU2TCLib.cpp', - 'fluid/CDataDrivenFluid.cpp']) + 'fluid/CDataDrivenFluid.cpp', + 'fluid/CFluidCantera.cpp']) su2_cfd_src += files(['output/COutputFactory.cpp', 'output/CAdjElasticityOutput.cpp', @@ -195,11 +196,16 @@ if get_option('enable-gprof') profiling_args = ['-pg','-no-pie'] endif + if get_option('enable-normal') + dependencies_default = [su2_deps, common_dep] + if get_option('enable-cantera') + dependencies_default += cantera_dep + endif su2_cfd_lib = static_library('SU2core', su2_cfd_src, install : false, - dependencies : [su2_deps, common_dep], + dependencies : dependencies_default, cpp_args: [default_warning_flags, su2_cpp_args]) su2_cfd_dep = declare_dependency(link_with: su2_cfd_lib, include_directories: su2_cfd_include) diff --git a/SU2_CFD/src/numerics/CNumerics.cpp b/SU2_CFD/src/numerics/CNumerics.cpp index bd8728e0e83..8f5b8f4c359 100644 --- a/SU2_CFD/src/numerics/CNumerics.cpp +++ b/SU2_CFD/src/numerics/CNumerics.cpp @@ -53,7 +53,9 @@ CNumerics::CNumerics(unsigned short val_nDim, unsigned short val_nVar, Gamma_Minus_One = Gamma - 1.0; Prandtl_Turb = config->GetPrandtl_Turb(); Gas_Constant = config->GetGas_ConstantND(); - energy_multicomponent = config->GetKind_FluidModel() == FLUID_MIXTURE && config->GetEnergy_Equation(); + energy_multicomponent = + (config->GetKind_FluidModel() == FLUID_MIXTURE || config->GetKind_FluidModel() == FLUID_CANTERA) && + config->GetEnergy_Equation(); tau = new su2double* [nDim]; for (iDim = 0; iDim < nDim; iDim++) diff --git a/SU2_CFD/src/numerics/flow/flow_diffusion.cpp b/SU2_CFD/src/numerics/flow/flow_diffusion.cpp index 66945b2f700..f2acb4fc358 100644 --- a/SU2_CFD/src/numerics/flow/flow_diffusion.cpp +++ b/SU2_CFD/src/numerics/flow/flow_diffusion.cpp @@ -580,7 +580,7 @@ CNumerics::ResidualType<> CAvgGradInc_Flow::ComputeResidual(const CConfig* confi implicit = (config->GetKind_TimeIntScheme() == EULER_IMPLICIT); AD::StartPreacc(); - AD::SetPreaccIn(V_i, nDim+9); AD::SetPreaccIn(V_j, nDim+9); + AD::SetPreaccIn(V_i, nDim+10); AD::SetPreaccIn(V_j, nDim+10); AD::SetPreaccIn(Coord_i, nDim); AD::SetPreaccIn(Coord_j, nDim); AD::SetPreaccIn(PrimVar_Grad_i, nVar, nDim); AD::SetPreaccIn(PrimVar_Grad_j, nVar, nDim); diff --git a/SU2_CFD/src/numerics/flow/flow_sources.cpp b/SU2_CFD/src/numerics/flow/flow_sources.cpp index 8cf9642b425..56818bed59f 100644 --- a/SU2_CFD/src/numerics/flow/flow_sources.cpp +++ b/SU2_CFD/src/numerics/flow/flow_sources.cpp @@ -245,7 +245,7 @@ CSourceIncAxisymmetric_Flow::CSourceIncAxisymmetric_Flow(unsigned short val_nDim implicit = (config->GetKind_TimeIntScheme_Flow() == EULER_IMPLICIT); energy = config->GetEnergy_Equation(); viscous = config->GetViscous(); - multicomponent = (config->GetKind_FluidModel() == FLUID_MIXTURE); + multicomponent = (config->GetKind_FluidModel() == FLUID_MIXTURE) || (config->GetKind_FluidModel() == FLUID_CANTERA); } CNumerics::ResidualType<> CSourceIncAxisymmetric_Flow::ComputeResidual(const CConfig* config) { diff --git a/SU2_CFD/src/numerics/species/species_sources.cpp b/SU2_CFD/src/numerics/species/species_sources.cpp index c0dbcc51721..3c02ce3e19e 100644 --- a/SU2_CFD/src/numerics/species/species_sources.cpp +++ b/SU2_CFD/src/numerics/species/species_sources.cpp @@ -135,7 +135,41 @@ CNumerics::ResidualType<> CSourceAxisymmetric_Species::ComputeResidual(const return ResidualType<>(residual, jacobian, nullptr); } +template +CSourceCombustion_Species::CSourceCombustion_Species(unsigned short val_nDim, unsigned short val_nVar, + const CConfig* config) + : CSourceBase_Species(val_nDim, val_nVar, config) {} + +template +CNumerics::ResidualType<> CSourceCombustion_Species::ComputeResidual(const CConfig* config) { + /*--- Preaccumulation ---*/ + AD::StartPreacc(); + AD::SetPreaccIn(Volume); + AD::SetPreaccIn(Chemical_Source_Term_i, nVar); + + /*--- Initialization. ---*/ + for (auto iVar = 0u; iVar < nVar; iVar++) { + residual[iVar] = 0.0; + for (auto jVar = 0; jVar < nVar; jVar++) { + jacobian[iVar][jVar] = 0.0; + } + } + + for (auto iVar = 0u; iVar < nVar; iVar++) { + residual[iVar] += Volume * Chemical_Source_Term_i[iVar]; + } + + + AD::SetPreaccOut(residual, nVar); + AD::EndPreacc(); + + return ResidualType<>(residual, jacobian, nullptr); +} + /*--- Explicit instantiations until we don't move this to the hpp. ---*/ template class CSourceAxisymmetric_Species >; template class CSourceAxisymmetric_Species >; template class CSourceAxisymmetric_Species >; +template class CSourceCombustion_Species >; +template class CSourceCombustion_Species >; +template class CSourceCombustion_Species >; diff --git a/SU2_CFD/src/output/CAdjFlowOutput.cpp b/SU2_CFD/src/output/CAdjFlowOutput.cpp index 12d22d46740..e5bd2920bf7 100644 --- a/SU2_CFD/src/output/CAdjFlowOutput.cpp +++ b/SU2_CFD/src/output/CAdjFlowOutput.cpp @@ -57,7 +57,11 @@ void CAdjFlowOutput::AddHistoryOutputFields_AdjScalarRMS_RES(const CConfig* conf if (config->GetKind_Species_Model() == SPECIES_MODEL::SPECIES_TRANSPORT) { for (unsigned short iVar = 0; iVar < config->GetnSpecies(); iVar++) { - AddHistoryOutput("RMS_ADJ_SPECIES_" + std::to_string(iVar), "rms[A_rho*Y_" + std::to_string(iVar) + "]", ScreenOutputFormat::FIXED, "RMS_RES", "Root-mean square residual of the adjoint transported species.", HistoryFieldType::RESIDUAL); + if(config->GetKind_FluidModel()==FLUID_CANTERA){ + AddHistoryOutput("RMS_ADJ_SPECIES_" + config->GetChemical_GasComposition(iVar), "rms[A_rho*Y_" + config->GetChemical_GasComposition(iVar) + "]", ScreenOutputFormat::FIXED, "RMS_RES", "Root-mean square residual of the adjoint transported species.", HistoryFieldType::RESIDUAL); + }else{ + AddHistoryOutput("RMS_ADJ_SPECIES_" + std::to_string(iVar), "rms[A_rho*Y_" + std::to_string(iVar) + "]", ScreenOutputFormat::FIXED, "RMS_RES", "Root-mean square residual of the adjoint transported species.", HistoryFieldType::RESIDUAL); + } } } @@ -94,7 +98,11 @@ void CAdjFlowOutput::AddHistoryOutputFields_AdjScalarMAX_RES(const CConfig* conf if (config->GetKind_Species_Model() == SPECIES_MODEL::SPECIES_TRANSPORT) { for (unsigned short iVar = 0; iVar < config->GetnSpecies(); iVar++) { - AddHistoryOutput("MAX_ADJ_SPECIES_" + std::to_string(iVar), "max[A_rho*Y_" + std::to_string(iVar) + "]",ScreenOutputFormat::FIXED, "MAX_RES", "Maximum residual of the adjoint transported species.", HistoryFieldType::RESIDUAL); + if(config->GetKind_FluidModel()==FLUID_CANTERA){ + AddHistoryOutput("MAX_ADJ_SPECIES_" + config->GetChemical_GasComposition(iVar), "max[A_rho*Y_" + config->GetChemical_GasComposition(iVar) + "]",ScreenOutputFormat::FIXED, "MAX_RES", "Maximum residual of the adjoint transported species.", HistoryFieldType::RESIDUAL); + }else{ + AddHistoryOutput("MAX_ADJ_SPECIES_" + std::to_string(iVar), "max[A_rho*Y_" + std::to_string(iVar) + "]",ScreenOutputFormat::FIXED, "MAX_RES", "Maximum residual of the adjoint transported species.", HistoryFieldType::RESIDUAL); + } } } @@ -204,10 +212,18 @@ void CAdjFlowOutput::LoadHistoryDataAdjScalar(const CConfig* config, const CSolv if (config->GetKind_Species_Model() == SPECIES_MODEL::SPECIES_TRANSPORT) { for (unsigned short iVar = 0; iVar < config->GetnSpecies(); iVar++) { - SetHistoryOutputValue("RMS_ADJ_SPECIES_" + std::to_string(iVar), log10(adjspecies_solver->GetRes_RMS(iVar))); - SetHistoryOutputValue("MAX_ADJ_SPECIES_" + std::to_string(iVar), log10(adjspecies_solver->GetRes_Max(iVar))); - if (multiZone) { - SetHistoryOutputValue("BGS_ADJ_SPECIES_" + std::to_string(iVar), log10(adjspecies_solver->GetRes_BGS(iVar))); + if (config->GetKind_FluidModel() == FLUID_CANTERA) { + SetHistoryOutputValue("RMS_ADJ_SPECIES_" + config->GetChemical_GasComposition(iVar), log10(adjspecies_solver->GetRes_RMS(iVar))); + SetHistoryOutputValue("MAX_ADJ_SPECIES_" + config->GetChemical_GasComposition(iVar), log10(adjspecies_solver->GetRes_Max(iVar))); + if (multiZone) { + SetHistoryOutputValue("BGS_ADJ_SPECIES_" + config->GetChemical_GasComposition(iVar), log10(adjspecies_solver->GetRes_BGS(iVar))); + } + } else { + SetHistoryOutputValue("RMS_ADJ_SPECIES_" + std::to_string(iVar), log10(adjspecies_solver->GetRes_RMS(iVar))); + SetHistoryOutputValue("MAX_ADJ_SPECIES_" + std::to_string(iVar), log10(adjspecies_solver->GetRes_Max(iVar))); + if (multiZone) { + SetHistoryOutputValue("BGS_ADJ_SPECIES_" + std::to_string(iVar), log10(adjspecies_solver->GetRes_BGS(iVar))); + } } } diff --git a/SU2_CFD/src/output/CFlowOutput.cpp b/SU2_CFD/src/output/CFlowOutput.cpp index 06f0534e11a..e39b8defca4 100644 --- a/SU2_CFD/src/output/CFlowOutput.cpp +++ b/SU2_CFD/src/output/CFlowOutput.cpp @@ -84,7 +84,11 @@ void CFlowOutput::AddAnalyzeSurfaceOutput(const CConfig *config){ if (config->GetKind_Species_Model() == SPECIES_MODEL::SPECIES_TRANSPORT) { /// DESCRIPTION: Average Species for (unsigned short iVar = 0; iVar < config->GetnSpecies(); iVar++) { - AddHistoryOutput("SURFACE_SPECIES_" + std::to_string(iVar), "Avg_Species_" + std::to_string(iVar), ScreenOutputFormat::FIXED, "SPECIES_COEFF", "Total average species " + std::to_string(iVar) + " on all markers set in MARKER_ANALYZE", HistoryFieldType::COEFFICIENT); + if(config->GetKind_FluidModel()==FLUID_CANTERA){ + AddHistoryOutput("SURFACE_SPECIES_" + config->GetChemical_GasComposition(iVar), "Avg_Species_" + config->GetChemical_GasComposition(iVar), ScreenOutputFormat::FIXED, "SPECIES_COEFF", "Total average species " + config->GetChemical_GasComposition(iVar) + " on all markers set in MARKER_ANALYZE", HistoryFieldType::COEFFICIENT); + }else{ + AddHistoryOutput("SURFACE_SPECIES_" + std::to_string(iVar), "Avg_Species_" + std::to_string(iVar), ScreenOutputFormat::FIXED, "SPECIES_COEFF", "Total average species " + std::to_string(iVar) + " on all markers set in MARKER_ANALYZE", HistoryFieldType::COEFFICIENT); + } } /// DESCRIPTION: Species Variance AddHistoryOutput("SURFACE_SPECIES_VARIANCE", "Species_Variance", ScreenOutputFormat::SCIENTIFIC, "SPECIES_COEFF", "Total species variance, measure for mixing quality. On all markers set in MARKER_ANALYZE", HistoryFieldType::COEFFICIENT); @@ -126,7 +130,11 @@ void CFlowOutput::AddAnalyzeSurfaceOutput(const CConfig *config){ if (config->GetKind_Species_Model() == SPECIES_MODEL::SPECIES_TRANSPORT) { /// DESCRIPTION: Average Species for (unsigned short iVar = 0; iVar < config->GetnSpecies(); iVar++) { - AddHistoryOutputPerSurface("SURFACE_SPECIES_" + std::to_string(iVar), "Avg_Species_" + std::to_string(iVar), ScreenOutputFormat::FIXED, "SPECIES_COEFF_SURF", Marker_Analyze, HistoryFieldType::COEFFICIENT); + if(config->GetKind_FluidModel()==FLUID_CANTERA){ + AddHistoryOutputPerSurface("SURFACE_SPECIES_" + config->GetChemical_GasComposition(iVar), "Avg_Species_" + config->GetChemical_GasComposition(iVar), ScreenOutputFormat::FIXED, "SPECIES_COEFF_SURF", Marker_Analyze, HistoryFieldType::COEFFICIENT); + }else{ + AddHistoryOutputPerSurface("SURFACE_SPECIES_" + std::to_string(iVar), "Avg_Species_" + std::to_string(iVar), ScreenOutputFormat::FIXED, "SPECIES_COEFF_SURF", Marker_Analyze, HistoryFieldType::COEFFICIENT); + } } /// DESCRIPTION: Species Variance AddHistoryOutputPerSurface("SURFACE_SPECIES_VARIANCE", "Species_Variance", ScreenOutputFormat::SCIENTIFIC, "SPECIES_COEFF_SURF", Marker_Analyze, HistoryFieldType::COEFFICIENT); @@ -506,7 +514,12 @@ void CFlowOutput::SetAnalyzeSurface(const CSolver* const*solver, const CGeometry if (species) { for (unsigned short iVar = 0; iVar < nSpecies; iVar++) { su2double Species = Surface_Species_Total(iMarker_Analyze, iVar); - SetHistoryOutputPerSurfaceValue("SURFACE_SPECIES_" + std::to_string(iVar), Species, iMarker_Analyze); + if (config->GetKind_FluidModel() == FLUID_CANTERA) { + SetHistoryOutputPerSurfaceValue("SURFACE_SPECIES_" + config->GetChemical_GasComposition(iVar), Species, + iMarker_Analyze); + } else { + SetHistoryOutputPerSurfaceValue("SURFACE_SPECIES_" + std::to_string(iVar), Species, iMarker_Analyze); + } Tot_Surface_Species[iVar] += Species; if (iVar == 0) config->SetSurface_Species_0(iMarker_Analyze, Species); @@ -542,8 +555,13 @@ void CFlowOutput::SetAnalyzeSurface(const CSolver* const*solver, const CGeometry SetHistoryOutputValue("SURFACE_TOTAL_TEMPERATURE", Tot_Surface_TotalTemperature); SetHistoryOutputValue("SURFACE_TOTAL_PRESSURE", Tot_Surface_TotalPressure); if (species) { - for (unsigned short iVar = 0; iVar < nSpecies; iVar++) - SetHistoryOutputValue("SURFACE_SPECIES_" + std::to_string(iVar), Tot_Surface_Species[iVar]); + for (unsigned short iVar = 0; iVar < nSpecies; iVar++) { + if (config->GetKind_FluidModel() == FLUID_CANTERA) { + SetHistoryOutputValue("SURFACE_SPECIES_" + config->GetChemical_GasComposition(iVar), Tot_Surface_Species[iVar]); + } else { + SetHistoryOutputValue("SURFACE_SPECIES_" + std::to_string(iVar), Tot_Surface_Species[iVar]); + } + } SetAnalyzeSurfaceSpeciesVariance(solver, geometry, config, Surface_Species_Total, Surface_MassFlow_Abs_Total, Surface_Area_Total); @@ -1064,7 +1082,11 @@ void CFlowOutput::AddHistoryOutputFields_ScalarRMS_RES(const CConfig* config) { switch (config->GetKind_Species_Model()) { case SPECIES_MODEL::SPECIES_TRANSPORT: { for (unsigned short iVar = 0; iVar < config->GetnSpecies(); iVar++) { - AddHistoryOutput("RMS_SPECIES_" + std::to_string(iVar), "rms[rho*Y_" + std::to_string(iVar)+"]", ScreenOutputFormat::FIXED, "RMS_RES", "Root-mean square residual of transported species.", HistoryFieldType::RESIDUAL); + if (config->GetKind_FluidModel()==FLUID_CANTERA){ + AddHistoryOutput("RMS_SPECIES_" + config->GetChemical_GasComposition(iVar), "rms[rho*Y_" + config->GetChemical_GasComposition(iVar)+"]", ScreenOutputFormat::FIXED, "RMS_RES", "Root-mean square residual of transported species.", HistoryFieldType::RESIDUAL); + }else{ + AddHistoryOutput("RMS_SPECIES_" + std::to_string(iVar), "rms[rho*Y_" + std::to_string(iVar)+"]", ScreenOutputFormat::FIXED, "RMS_RES", "Root-mean square residual of transported species.", HistoryFieldType::RESIDUAL); + } } break; } @@ -1130,7 +1152,11 @@ void CFlowOutput::AddHistoryOutputFields_ScalarMAX_RES(const CConfig* config) { switch (config->GetKind_Species_Model()) { case SPECIES_MODEL::SPECIES_TRANSPORT: { for (unsigned short iVar = 0; iVar < config->GetnSpecies(); iVar++) { - AddHistoryOutput("MAX_SPECIES_" + std::to_string(iVar), "max[rho*Y_" + std::to_string(iVar)+"]", ScreenOutputFormat::FIXED, "MAX_RES", "Maximum residual of transported species.", HistoryFieldType::RESIDUAL); + if (config->GetKind_FluidModel()==FLUID_CANTERA){ + AddHistoryOutput("MAX_SPECIES_" + config->GetChemical_GasComposition(iVar), "max[rho*Y_" + config->GetChemical_GasComposition(iVar)+"]", ScreenOutputFormat::FIXED, "MAX_RES", "Maximum residual of transported species.", HistoryFieldType::RESIDUAL); + }else{ + AddHistoryOutput("MAX_SPECIES_" + std::to_string(iVar), "max[rho*Y_" + std::to_string(iVar)+"]", ScreenOutputFormat::FIXED, "MAX_RES", "Maximum residual of transported species.", HistoryFieldType::RESIDUAL); + } } break; } @@ -1194,7 +1220,11 @@ void CFlowOutput::AddHistoryOutputFields_ScalarBGS_RES(const CConfig* config) { switch (config->GetKind_Species_Model()) { case SPECIES_MODEL::SPECIES_TRANSPORT: { for (unsigned short iVar = 0; iVar < config->GetnSpecies(); iVar++) { - AddHistoryOutput("BGS_SPECIES_" + std::to_string(iVar), "bgs[rho*Y_" + std::to_string(iVar)+"]", ScreenOutputFormat::FIXED, "BGS_RES", "Maximum residual of transported species.", HistoryFieldType::RESIDUAL); + if (config->GetKind_FluidModel()==FLUID_CANTERA){ + AddHistoryOutput("BGS_SPECIES_" + config->GetChemical_GasComposition(iVar), "bgs[rho*Y_" + config->GetChemical_GasComposition(iVar)+"]", ScreenOutputFormat::FIXED, "BGS_RES", "Maximum residual of transported species.", HistoryFieldType::RESIDUAL); + }else{ + AddHistoryOutput("BGS_SPECIES_" + std::to_string(iVar), "bgs[rho*Y_" + std::to_string(iVar)+"]", ScreenOutputFormat::FIXED, "BGS_RES", "Maximum residual of transported species.", HistoryFieldType::RESIDUAL); + } } break; } @@ -1307,10 +1337,18 @@ void CFlowOutput::LoadHistoryDataScalar(const CConfig* config, const CSolver* co switch(config->GetKind_Species_Model()) { case SPECIES_MODEL::SPECIES_TRANSPORT: { for (unsigned short iVar = 0; iVar < config->GetnSpecies(); iVar++) { - SetHistoryOutputValue("RMS_SPECIES_" + std::to_string(iVar), log10(solver[SPECIES_SOL]->GetRes_RMS(iVar))); - SetHistoryOutputValue("MAX_SPECIES_" + std::to_string(iVar), log10(solver[SPECIES_SOL]->GetRes_Max(iVar))); - if (multiZone) { - SetHistoryOutputValue("BGS_SPECIES_" + std::to_string(iVar), log10(solver[SPECIES_SOL]->GetRes_BGS(iVar))); + if (config->GetKind_FluidModel()==FLUID_CANTERA){ + SetHistoryOutputValue("RMS_SPECIES_" + config->GetChemical_GasComposition(iVar), log10(solver[SPECIES_SOL]->GetRes_RMS(iVar))); + SetHistoryOutputValue("MAX_SPECIES_" + config->GetChemical_GasComposition(iVar), log10(solver[SPECIES_SOL]->GetRes_Max(iVar))); + if (multiZone) { + SetHistoryOutputValue("BGS_SPECIES_" + config->GetChemical_GasComposition(iVar), log10(solver[SPECIES_SOL]->GetRes_BGS(iVar))); + } + }else{ + SetHistoryOutputValue("RMS_SPECIES_" + std::to_string(iVar), log10(solver[SPECIES_SOL]->GetRes_RMS(iVar))); + SetHistoryOutputValue("MAX_SPECIES_" + std::to_string(iVar), log10(solver[SPECIES_SOL]->GetRes_Max(iVar))); + if (multiZone) { + SetHistoryOutputValue("BGS_SPECIES_" + std::to_string(iVar), log10(solver[SPECIES_SOL]->GetRes_BGS(iVar))); + } } } SetHistoryOutputValue("LINSOL_ITER_SPECIES", solver[SPECIES_SOL]->GetIterLinSolver()); @@ -1382,8 +1420,14 @@ void CFlowOutput::SetVolumeOutputFieldsScalarSolution(const CConfig* config){ switch (config->GetKind_Species_Model()) { case SPECIES_MODEL::SPECIES_TRANSPORT: - for (unsigned short iVar = 0; iVar < config->GetnSpecies(); iVar++){ - AddVolumeOutput("SPECIES_" + std::to_string(iVar), "Species_" + std::to_string(iVar), "SOLUTION", "Species_" + std::to_string(iVar) + " mass fraction"); + for (unsigned short iVar = 0; iVar < config->GetnSpecies(); iVar++) { + if (config->GetKind_FluidModel() == FLUID_CANTERA) { + AddVolumeOutput("SPECIES_" + config->GetChemical_GasComposition(iVar), "Species_" + config->GetChemical_GasComposition(iVar), "SOLUTION", + "Species_" + config->GetChemical_GasComposition(iVar) + " mass fraction"); + } else { + AddVolumeOutput("SPECIES_" + std::to_string(iVar), "Species_" + std::to_string(iVar), "SOLUTION", + "Species_" + std::to_string(iVar) + " mass fraction"); + } } break; case SPECIES_MODEL::FLAMELET: { @@ -1430,7 +1474,13 @@ void CFlowOutput::SetVolumeOutputFieldsScalarResidual(const CConfig* config) { switch (config->GetKind_Species_Model()) { case SPECIES_MODEL::SPECIES_TRANSPORT: for (unsigned short iVar = 0; iVar < config->GetnSpecies(); iVar++){ - AddVolumeOutput("RES_SPECIES_" + std::to_string(iVar), "Residual_Species_" + std::to_string(iVar), "RESIDUAL", "Residual of the transported species " + std::to_string(iVar)); + if (config->GetKind_FluidModel() == FLUID_CANTERA) { + AddVolumeOutput("RES_SPECIES_" + config->GetChemical_GasComposition(iVar), + "Residual_Species_" + config->GetChemical_GasComposition(iVar), "RESIDUAL", + "Residual of the transported species " + config->GetChemical_GasComposition(iVar)); + } else { + AddVolumeOutput("RES_SPECIES_" + std::to_string(iVar), "Residual_Species_" + std::to_string(iVar), "RESIDUAL", "Residual of the transported species " + std::to_string(iVar)); + } } break; case SPECIES_MODEL::FLAMELET: { @@ -1515,8 +1565,20 @@ void CFlowOutput::SetVolumeOutputFieldsScalarPrimitive(const CConfig* config) { switch (config->GetKind_Species_Model()) { case SPECIES_MODEL::SPECIES_TRANSPORT: for (unsigned short iVar = 0; iVar < config->GetnSpecies(); iVar++){ - AddVolumeOutput("DIFFUSIVITY_" + std::to_string(iVar), "Diffusivity_" + std::to_string(iVar), "PRIMITIVE", "Diffusivity of the transported species " + std::to_string(iVar)); + if (config->GetKind_FluidModel() == FLUID_CANTERA) { + AddVolumeOutput("DIFFUSIVITY_" + config->GetChemical_GasComposition(iVar), "Diffusivity_" + config->GetChemical_GasComposition(iVar), "PRIMITIVE", + "Diffusivity of the transported species " + config->GetChemical_GasComposition(iVar)); + if (config->GetCombustion() == true) { + AddVolumeOutput("CHEMICAL_SOURCE_TERM_" + config->GetChemical_GasComposition(iVar), + "Chemical_Source_Term_" + config->GetChemical_GasComposition(iVar), "PRIMITIVE", + "Chemical source term of the transported species " + config->GetChemical_GasComposition(iVar)); + } + } else { + AddVolumeOutput("DIFFUSIVITY_" + std::to_string(iVar), "Diffusivity_" + std::to_string(iVar), "PRIMITIVE", "Diffusivity of the transported species " + std::to_string(iVar)); + } } + if (config->GetCombustion() == true) + AddVolumeOutput("HEAT_RELEASE", "Heat_Release", "PRIMITIVE", "Heat release due to combustion"); break; default: break; @@ -1723,16 +1785,32 @@ void CFlowOutput::LoadVolumeDataScalar(const CConfig* config, const CSolver* con case SPECIES_MODEL::SPECIES_TRANSPORT: { const auto Node_Species = solver[SPECIES_SOL]->GetNodes(); - for (unsigned long iVar = 0; iVar < config->GetnSpecies(); iVar++) { - SetVolumeOutputValue("SPECIES_" + std::to_string(iVar), iPoint, Node_Species->GetSolution(iPoint, iVar)); - SetVolumeOutputValue("RES_SPECIES_" + std::to_string(iVar), iPoint, solver[SPECIES_SOL]->LinSysRes(iPoint, iVar)); - SetVolumeOutputValue("DIFFUSIVITY_"+ std::to_string(iVar), iPoint, Node_Species->GetDiffusivity(iPoint,iVar)); - if (config->GetKind_SlopeLimit_Species() != LIMITER::NONE) - SetVolumeOutputValue("LIMITER_SPECIES_" + std::to_string(iVar), iPoint, Node_Species->GetLimiter(iPoint, iVar)); - if (config->GetPyCustomSource()){ - SetVolumeOutputValue("SPECIES_UDS_" + std::to_string(iVar), iPoint, Node_Species->GetUserDefinedSource()(iPoint, iVar)); + for (unsigned short iVar = 0; iVar < config->GetnSpecies(); iVar++) { + if (config->GetKind_FluidModel() == FLUID_CANTERA) { + SetVolumeOutputValue("SPECIES_" + config->GetChemical_GasComposition(iVar), iPoint, Node_Species->GetSolution(iPoint, iVar)); + SetVolumeOutputValue("RES_SPECIES_" + config->GetChemical_GasComposition(iVar), iPoint, + solver[SPECIES_SOL]->LinSysRes(iPoint, iVar)); + SetVolumeOutputValue("DIFFUSIVITY_" + config->GetChemical_GasComposition(iVar), iPoint, + Node_Species->GetDiffusivity(iPoint, iVar)); + if (config->GetCombustion() == true) + SetVolumeOutputValue("CHEMICAL_SOURCE_TERM_" + config->GetChemical_GasComposition(iVar), iPoint, + Node_Species->GetChemicalSourceTerm(iPoint, iVar)); + if (config->GetKind_SlopeLimit_Species() != LIMITER::NONE) + SetVolumeOutputValue("LIMITER_SPECIES_" + config->GetChemical_GasComposition(iVar), iPoint, + Node_Species->GetLimiter(iPoint, iVar)); + } else { + SetVolumeOutputValue("SPECIES_" + std::to_string(iVar), iPoint, Node_Species->GetSolution(iPoint, iVar)); + SetVolumeOutputValue("RES_SPECIES_" + std::to_string(iVar), iPoint, solver[SPECIES_SOL]->LinSysRes(iPoint, iVar)); + SetVolumeOutputValue("DIFFUSIVITY_" + std::to_string(iVar), iPoint, Node_Species->GetDiffusivity(iPoint, iVar)); + if (config->GetKind_SlopeLimit_Species() != LIMITER::NONE) + SetVolumeOutputValue("LIMITER_SPECIES_" + std::to_string(iVar), iPoint, Node_Species->GetLimiter(iPoint, iVar)); + if (config->GetPyCustomSource()) { + SetVolumeOutputValue("SPECIES_UDS_" + std::to_string(iVar), iPoint, Node_Species->GetUserDefinedSource()(iPoint, iVar)); + } } } + if (config->GetCombustion() == true) + SetVolumeOutputValue("HEAT_RELEASE", iPoint, Node_Species->GetHeatRelease(iPoint)); break; } @@ -3316,6 +3394,13 @@ void CFlowOutput::WriteForcesBreakdown(const CConfig* config, const CSolver* flo case VISCOSITYMODEL::COOLPROP: file << "Viscosity Model: CoolProp \n"; break; + + case VISCOSITYMODEL::CANTERA: + file << "Viscosity Model: CANTERA \n"; + if (si_units) file << " N.s/m^2.\n"; + else file << " lbf.s/ft^2.\n"; + file << "Laminar Viscosity (non-dim): " << config->GetMu_ConstantND() << "\n"; + break; case VISCOSITYMODEL::SUTHERLAND: file << "Viscosity Model: SUTHERLAND \n"; @@ -3371,6 +3456,12 @@ void CFlowOutput::WriteForcesBreakdown(const CConfig* config, const CSolver* flo file << "Molecular Conductivity units: " << " W/m^2.K.\n"; file << "Molecular Conductivity (non-dim): " << config->GetThermal_Conductivity_ConstantND() << "\n"; break; + + case CONDUCTIVITYMODEL::CANTERA: + file << "Conductivity Model: CANTERA \n"; + file << "Molecular Conductivity units: " << " W/m^2.K.\n"; + file << "Molecular Conductivity (non-dim): " << config->GetThermal_Conductivity_ConstantND() << "\n"; + break; case CONDUCTIVITYMODEL::POLYNOMIAL: file << "Viscosity Model: POLYNOMIAL \n"; diff --git a/SU2_CFD/src/solvers/CIncEulerSolver.cpp b/SU2_CFD/src/solvers/CIncEulerSolver.cpp index 41580c3062a..06a0bf21944 100644 --- a/SU2_CFD/src/solvers/CIncEulerSolver.cpp +++ b/SU2_CFD/src/solvers/CIncEulerSolver.cpp @@ -33,6 +33,7 @@ #include "../../include/variables/CIncNSVariable.hpp" #include "../../../Common/include/toolboxes/geometry_toolbox.hpp" #include "../../include/fluid/CFluidScalar.hpp" +#include "../../include/fluid/CFluidCantera.hpp" #include "../../include/fluid/CFluidFlamelet.hpp" #include "../../include/fluid/CFluidModel.hpp" @@ -339,6 +340,14 @@ void CIncEulerSolver::SetNondimensionalization(CConfig *config, unsigned short i auxFluidModel->SetTDState_T(Temperature_FreeStream, config->GetSpecies_Init()); break; + case FLUID_CANTERA: + + config->SetGas_Constant(UNIVERSAL_GAS_CONSTANT / (config->GetMolecular_Weight() / 1000.0)); + Pressure_Thermodynamic = config->GetPressure_Thermodynamic(); + auxFluidModel = new CFluidCantera(Pressure_Thermodynamic, config); + auxFluidModel->SetTDState_T(Temperature_FreeStream, config->GetSpecies_Init()); + break; + default: SU2_MPI::Error("Fluid model not implemented for incompressible solver.", CURRENT_FUNCTION); @@ -501,6 +510,11 @@ void CIncEulerSolver::SetNondimensionalization(CConfig *config, unsigned short i fluidModel->SetTDState_T(Temperature_FreeStreamND, config->GetSpecies_Init()); break; + case FLUID_CANTERA: + fluidModel = new CFluidCantera(Pressure_ThermodynamicND, config); + fluidModel->SetTDState_T(Temperature_FreeStreamND, config->GetSpecies_Init()); + break; + case INC_IDEAL_GAS_POLY: fluidModel = new CIncIdealGasPolynomial(Gas_ConstantND, Pressure_ThermodynamicND, STD_REF_TEMP / config->GetTemperature_Ref()); if (viscous) { @@ -668,6 +682,15 @@ void CIncEulerSolver::SetNondimensionalization(CConfig *config, unsigned short i Unit.str(""); NonDimTable.PrintFooter(); break; + + case VISCOSITYMODEL::CANTERA: + ModelTable << "CANTERA"; + if (config->GetSystemMeasurements() == SI) Unit << "N.s/m^2"; + else if (config->GetSystemMeasurements() == US) Unit << "lbf.s/ft^2"; + NonDimTable << "Viscosity" << "--" << "--" << Unit.str() << config->GetMu_ConstantND(); + Unit.str(""); + NonDimTable.PrintFooter(); + break; case VISCOSITYMODEL::SUTHERLAND: ModelTable << "SUTHERLAND"; @@ -731,6 +754,13 @@ void CIncEulerSolver::SetNondimensionalization(CConfig *config, unsigned short i Unit.str(""); NonDimTable.PrintFooter(); break; + + case CONDUCTIVITYMODEL::CANTERA: + ModelTable << "CANTERA"; + Unit << "W/m^2.K"; + NonDimTable << "Molecular Cond." << "--" << "--" << Unit.str() << config->GetThermal_Conductivity_ConstantND(); + Unit.str(""); + break; case CONDUCTIVITYMODEL::POLYNOMIAL: ModelTable << "POLYNOMIAL"; @@ -818,6 +848,23 @@ void CIncEulerSolver::SetNondimensionalization(CConfig *config, unsigned short i NonDimTable.PrintFooter(); break; + case FLUID_CANTERA: + ModelTable << "CANTERA"; + Unit << "N.m/kg.K"; + NonDimTable << "Spec. Heat (Cp)" << config->GetSpecific_Heat_Cp() << config->GetSpecific_Heat_Cp() / config->GetSpecific_Heat_CpND() << Unit.str() << config->GetSpecific_Heat_CpND(); + Unit.str(""); + Unit << "g/mol"; + NonDimTable << "Molecular weight" << config->GetMolecular_Weight() << 1.0 << Unit.str() << config->GetMolecular_Weight(); + Unit.str(""); + Unit << "N.m/kg.K"; + NonDimTable << "Gas Constant" << config->GetGas_Constant() << config->GetGas_Constant_Ref() << Unit.str() << config->GetGas_ConstantND(); + Unit.str(""); + Unit << "Pa"; + NonDimTable << "Therm. Pressure" << config->GetPressure_Thermodynamic() << config->GetPressure_Ref() << Unit.str() << config->GetPressure_ThermodynamicND(); + Unit.str(""); + NonDimTable.PrintFooter(); + break; + case INC_IDEAL_GAS_POLY: ModelTable << "INC_IDEAL_GAS_POLY"; Unit.str(""); @@ -1297,7 +1344,7 @@ void CIncEulerSolver::Upwind_Residual(CGeometry *geometry, CSolver **solver_cont const bool limiter = (config->GetKind_SlopeLimit_Flow() != LIMITER::NONE); const bool van_albada = (config->GetKind_SlopeLimit_Flow() == LIMITER::VAN_ALBADA_EDGE); const bool bounded_scalar = config->GetBounded_Scalar(); - const bool multicomponent = (config->GetKind_FluidModel() == FLUID_MIXTURE); + const bool multicomponent = ((config->GetKind_FluidModel() == FLUID_MIXTURE) || (config->GetKind_FluidModel() == FLUID_CANTERA)); const su2double kappa = config->GetMUSCL_Kappa_Flow(); const su2double musclRamp = config->GetMUSCLRampValue() * config->GetNewtonKrylovRelaxation(); @@ -1473,7 +1520,8 @@ void CIncEulerSolver::Source_Residual(CGeometry *geometry, CSolver **solver_cont const bool energy = config->GetEnergy_Equation(); const bool streamwise_periodic = (config->GetKind_Streamwise_Periodic() != ENUM_STREAMWISE_PERIODIC::NONE); const bool streamwise_periodic_temperature = config->GetStreamwise_Periodic_Temperature(); - const bool multicomponent = (config->GetKind_FluidModel() == FLUID_MIXTURE); + const bool multicomponent = + ((config->GetKind_FluidModel() == FLUID_MIXTURE) || (config->GetKind_FluidModel() == FLUID_CANTERA)); AD::StartNoSharedReading(); @@ -2249,7 +2297,9 @@ void CIncEulerSolver::BC_Far_Field(CGeometry *geometry, CSolver **solver_contain const bool implicit = config->GetKind_TimeIntScheme() == EULER_IMPLICIT; const bool viscous = config->GetViscous(); - const bool energy_multicomponent = config->GetKind_FluidModel() == FLUID_MIXTURE && config->GetEnergy_Equation(); + const bool energy_multicomponent = + (config->GetKind_FluidModel() == FLUID_MIXTURE || config->GetKind_FluidModel() == FLUID_CANTERA) && + config->GetEnergy_Equation(); const bool species_model = config->GetKind_Species_Model() != SPECIES_MODEL::NONE; su2double Normal[MAXNDIM] = {0.0}; @@ -2393,7 +2443,9 @@ void CIncEulerSolver::BC_Inlet(CGeometry *geometry, CSolver **solver_container, const bool implicit = (config->GetKind_TimeIntScheme() == EULER_IMPLICIT); const bool viscous = config->GetViscous(); - const bool energy_multicomponent = config->GetKind_FluidModel() == FLUID_MIXTURE && config->GetEnergy_Equation(); + const bool energy_multicomponent = + (config->GetKind_FluidModel() == FLUID_MIXTURE || config->GetKind_FluidModel() == FLUID_CANTERA) && + config->GetEnergy_Equation(); const bool species_model = config->GetKind_Species_Model() != SPECIES_MODEL::NONE; string Marker_Tag = config->GetMarker_All_TagBound(val_marker); @@ -2663,7 +2715,9 @@ void CIncEulerSolver::BC_Outlet(CGeometry *geometry, CSolver **solver_container, const bool implicit = (config->GetKind_TimeIntScheme() == EULER_IMPLICIT); const bool viscous = config->GetViscous(); - const bool energy_multicomponent = config->GetKind_FluidModel() == FLUID_MIXTURE && config->GetEnergy_Equation(); + const bool energy_multicomponent = + (config->GetKind_FluidModel() == FLUID_MIXTURE || config->GetKind_FluidModel() == FLUID_CANTERA) && + config->GetEnergy_Equation(); string Marker_Tag = config->GetMarker_All_TagBound(val_marker); su2double Normal[MAXNDIM] = {0.0}; diff --git a/SU2_CFD/src/solvers/CIncNSSolver.cpp b/SU2_CFD/src/solvers/CIncNSSolver.cpp index 2532215cf8d..594857e7f0f 100644 --- a/SU2_CFD/src/solvers/CIncNSSolver.cpp +++ b/SU2_CFD/src/solvers/CIncNSSolver.cpp @@ -43,6 +43,8 @@ CIncNSSolver::CIncNSSolver(CGeometry *geometry, CConfig *config, unsigned short Viscosity_Inf = config->GetViscosity_FreeStreamND(); Tke_Inf = config->GetTke_FreeStreamND(); + TemperatureLimits[0]= config->GetTemperatureLimits(0); + TemperatureLimits[1]= config->GetTemperatureLimits(1); /*--- Initialize the secondary values for direct derivative approximations ---*/ @@ -53,6 +55,7 @@ CIncNSSolver::CIncNSSolver(CGeometry *geometry, CConfig *config, unsigned short default: break; } + if (config->GetCombustion()) flamelet_config_options = config->GetFlameletParsedOptions(); /*--- Set the initial Streamwise periodic pressure drop value. ---*/ @@ -71,6 +74,54 @@ void CIncNSSolver::Preprocessing(CGeometry *geometry, CSolver **solver_container const bool limiter = (config->GetKind_SlopeLimit_Flow() != LIMITER::NONE) && (InnerIter <= config->GetLimiterIter()); const bool van_albada = (config->GetKind_SlopeLimit_Flow() == LIMITER::VAN_ALBADA_EDGE); const bool wall_functions = config->GetWall_Functions(); + const bool energy = config->GetEnergy_Equation(); + const bool combustion = config->GetCombustion(); + + /*--- Setting temperature, enthalpy and thermophysical properties for ignition in reacting flows. ---*/ + if (energy && combustion) { + unsigned long spark_iter_start, spark_duration; + bool ignition = false; + + /*--- Retrieve spark ignition parameters for spark-type ignition. ---*/ + if (flamelet_config_options.ignition_method == FLAMELET_INIT_TYPE::SPARK) { + auto spark_init = flamelet_config_options.spark_init; + 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))); + } + + SU2_OMP_SAFE_GLOBAL_ACCESS(config->SetGlobalParam(config->GetKind_Solver(), RunTime_EqSystem);) + + if (ignition) { + SU2_OMP_FOR_STAT(omp_chunk_size) + for (auto i_point = 0u; i_point < nPoint; i_point++) { + /*--- Retrieve fluid model. ---*/ + CFluidModel* fluid_model_local = solver_container[FLOW_SOL]->GetFluidModel(); + /*--- Apply ignition temperature within spark radius. ---*/ + su2double dist_from_center = 0, spark_radius = flamelet_config_options.spark_init[3]; + dist_from_center = + GeometryToolbox::SquaredDistance(nDim, geometry->nodes->GetCoord(i_point), flamelet_config_options.spark_init.data()); + if (dist_from_center < pow(spark_radius, 2)) { + /*--- Retrieve scalars solution. ---*/ + su2double* scalars = solver_container[SPECIES_SOL]->GetNodes()->GetSolution(i_point); + /*--- Set high temperature for ignition. ---*/ + nodes->SetTemperature(i_point, config->GetSpark_Temperature(), TemperatureLimits); + /*--- Set thermodynamic state at high temperature. ---*/ + fluid_model_local->SetTDState_T(config->GetSpark_Temperature(), scalars); + /*--- Set total enthalpy at high temperature. ---*/ + nodes->SetSolution(i_point, nDim + 1, fluid_model_local->GetEnthalpy()); + /*--- Set thermodynamics and transport properties at high temperature for consistency. ---*/ + nodes->SetDensity(i_point, fluid_model_local->GetDensity()); + nodes->SetSpecificHeatCp(i_point, fluid_model_local->GetCp()); + nodes->SetSpecificHeatCv(i_point, fluid_model_local->GetCv()); + nodes->SetThermalConductivity(i_point, fluid_model_local->GetThermalConductivity()); + nodes->SetLaminarViscosity(i_point, fluid_model_local->GetLaminarViscosity()); + } + } + END_SU2_OMP_FOR + } + } /*--- Common preprocessing steps (implemented by CEulerSolver) ---*/ @@ -278,7 +329,9 @@ void CIncNSSolver::Compute_Streamwise_Periodic_Recovered_Values(CConfig *config, CNumerics::ResidualType<> CIncNSSolver::Viscous_Residual(unsigned long iEdge, CGeometry *geometry, CSolver **solver_container, CNumerics *numerics, CConfig *config) { - const bool energy_multicomponent = config->GetKind_FluidModel() == FLUID_MIXTURE && config->GetEnergy_Equation(); + const bool energy_multicomponent = + (config->GetKind_FluidModel() == FLUID_MIXTURE || config->GetKind_FluidModel() == FLUID_CANTERA) && + config->GetEnergy_Equation(); /*--- Contribution to heat flux due to enthalpy diffusion for multicomponent and reacting flows ---*/ if (energy_multicomponent) { diff --git a/SU2_CFD/src/solvers/CSpeciesSolver.cpp b/SU2_CFD/src/solvers/CSpeciesSolver.cpp index 38ad833f65e..33cc02dee2d 100644 --- a/SU2_CFD/src/solvers/CSpeciesSolver.cpp +++ b/SU2_CFD/src/solvers/CSpeciesSolver.cpp @@ -43,6 +43,8 @@ CSpeciesSolver::CSpeciesSolver(CGeometry* geometry, CConfig* config, unsigned sh nVar = config->GetnSpecies(); + if (config->GetCombustion()) flamelet_config_options = config->GetFlameletParsedOptions(); + Initialize(geometry, config, iMesh, nVar); /*--- Initialize the solution to the far-field state everywhere. ---*/ @@ -310,20 +312,57 @@ void CSpeciesSolver::Preprocessing(CGeometry* geometry, CSolver** solver_contain unsigned short iMesh, unsigned short iRKStep, unsigned short RunTime_EqSystem, bool Output) { SU2_ZONE_SCOPED + unsigned long spark_iter_start, spark_duration; + bool ignition = false; + su2double temperature; + + /*--- Retrieve spark ignition parameters for spark-type ignition. ---*/ + if (flamelet_config_options.ignition_method == FLAMELET_INIT_TYPE::SPARK) { + auto spark_init = flamelet_config_options.spark_init; + 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))); + } SU2_OMP_SAFE_GLOBAL_ACCESS(config->SetGlobalParam(config->GetKind_Solver(), RunTime_EqSystem);) + /*--- Set the laminar mass Diffusivity and chemical source term for the species solver. ---*/ SU2_OMP_FOR_STAT(omp_chunk_size) for (auto iPoint = 0u; iPoint < nPoint; iPoint++) { - const su2double temperature = solver_container[FLOW_SOL]->GetNodes()->GetTemperature(iPoint); + if (ignition) { + /*--- Apply ignition temperature within spark radius. ---*/ + su2double dist_from_center = 0, spark_radius = flamelet_config_options.spark_init[3]; + dist_from_center = + GeometryToolbox::SquaredDistance(nDim, geometry->nodes->GetCoord(iPoint), flamelet_config_options.spark_init.data()); + if (dist_from_center < pow(spark_radius, 2)) { + temperature = config->GetSpark_Temperature(); + } else { + temperature = solver_container[FLOW_SOL]->GetNodes()->GetTemperature(iPoint); + } + } else { + temperature = solver_container[FLOW_SOL]->GetNodes()->GetTemperature(iPoint); + } const su2double* scalar = solver_container[SPECIES_SOL]->GetNodes()->GetSolution(iPoint); solver_container[FLOW_SOL]->GetFluidModel()->SetMassDiffusivityModel(config); solver_container[FLOW_SOL]->GetFluidModel()->SetTDState_T(temperature, scalar); + if (config->GetCombustion() == true) { + /*--- Call function to compute chemical source term. ---*/ + solver_container[FLOW_SOL]->GetFluidModel()->ComputeChemicalSourceTerm(); + /*--- Get and set heat release due to combustion. ---*/ + const su2double heat_release = solver_container[FLOW_SOL]->GetFluidModel()->GetHeatRelease(); + nodes->SetHeatRelease(iPoint, heat_release); + } /*--- Recompute viscosity, important to get diffusivity correct across MPI ranks. ---*/ nodes->SetLaminarViscosity(iPoint, solver_container[FLOW_SOL]->GetFluidModel()->GetLaminarViscosity()); /*--- Set the laminar mass Diffusivity for the species solver. ---*/ for (auto iVar = 0u; iVar <= nVar; iVar++) { const su2double mass_diffusivity = solver_container[FLOW_SOL]->GetFluidModel()->GetMassDiffusivity(iVar); nodes->SetDiffusivity(iPoint, mass_diffusivity, iVar); + if (config->GetCombustion() == true) { + /*--- Get and Set chemical source term. ---*/ + const su2double chemical_source_term = solver_container[FLOW_SOL]->GetFluidModel()->GetChemicalSourceTerm(iVar); + nodes->SetChemicalSourceTerm(iPoint, chemical_source_term, iVar); + } } } // iPoint END_SU2_OMP_FOR @@ -623,6 +662,7 @@ void CSpeciesSolver::Source_Residual(CGeometry *geometry, CSolver **solver_conta const bool implicit = (config->GetKind_TimeIntScheme() == EULER_IMPLICIT); const bool axisymmetric = config->GetAxisymmetric(); + const bool combustion =config->GetCombustion(); if (axisymmetric) { CNumerics *numerics = numerics_container[SOURCE_FIRST_TERM + omp_get_thread_num()*MAX_TERMS]; @@ -672,7 +712,34 @@ void CSpeciesSolver::Source_Residual(CGeometry *geometry, CSolver **solver_conta if (config->GetPyCustomSource()) { CustomSourceResidual(geometry, solver_container, numerics_container, config, iMesh); } + + if (combustion) { + CNumerics* numerics = numerics_container[SOURCE_SECOND_TERM + omp_get_thread_num() * MAX_TERMS]; + + SU2_OMP_FOR_DYN(omp_chunk_size) + for (auto iPoint = 0u; iPoint < nPointDomain; iPoint++) { + /*--- Set Chemical Source Term ---*/ + + numerics->SetChemicalSourceTerm(nodes->GetChemicalSourceTerm(iPoint), nullptr); + + /*--- Set volume of the dual cell. ---*/ + + numerics->SetVolume(geometry->nodes->GetVolume(iPoint)); + /*--- Update scalar sources in the fluidmodel ---*/ + + auto residual = numerics->ComputeResidual(config); + + /*--- Add Residual ---*/ + + LinSysRes.SubtractBlock(iPoint, residual); + + /*--- Implicit part ---*/ + + if (implicit) Jacobian.SubtractBlock2Diag(iPoint, residual.jacobian_i); + } + END_SU2_OMP_FOR + } } void CSpeciesSolver::SetInitialCondition(CGeometry **geometry, CSolver ***solver_container, CConfig *config, unsigned long TimeIter) { diff --git a/SU2_CFD/src/variables/CIncNSVariable.cpp b/SU2_CFD/src/variables/CIncNSVariable.cpp index 93d48195694..444ab38f4af 100644 --- a/SU2_CFD/src/variables/CIncNSVariable.cpp +++ b/SU2_CFD/src/variables/CIncNSVariable.cpp @@ -70,7 +70,7 @@ bool CIncNSVariable::SetPrimVar(unsigned long iPoint, su2double eddy_visc, su2do /*--- Set the value of the density ---*/ - const auto check_dens = SetDensity(iPoint, FluidModel->GetDensity()); + const auto check_dens = check_temp ? true: SetDensity(iPoint, FluidModel->GetDensity()); /*--- Non-physical solution found. Revert to old values. ---*/ diff --git a/SU2_CFD/src/variables/CSpeciesVariable.cpp b/SU2_CFD/src/variables/CSpeciesVariable.cpp index 34dfe46b69e..5f339d6f0e1 100644 --- a/SU2_CFD/src/variables/CSpeciesVariable.cpp +++ b/SU2_CFD/src/variables/CSpeciesVariable.cpp @@ -30,8 +30,12 @@ CSpeciesVariable::CSpeciesVariable(const su2double* species_inf, unsigned long npoint, unsigned long ndim, unsigned long nvar, const CConfig* config) : CScalarVariable(npoint, ndim, nvar, config) { - /*--- Allocate space for the mass diffusivity. ---*/ + /*--- Allocate space for the mass diffusivity and chemical source term. ---*/ Diffusivity.resize(nPoint, nVar + 1) = su2double(0.0); + if (config->GetCombustion()) { + SpeciesSourceTerm.resize(nPoint, nVar + 1) = su2double(0.0); + HeatRelease.resize(nPoint) = su2double(0.0); + } for (unsigned long iPoint = 0; iPoint < nPoint; iPoint++) for (unsigned long iVar = 0; iVar < nVar; iVar++) diff --git a/SU2_PY/SU2/io/historyMap.py b/SU2_PY/SU2/io/historyMap.py index 57392eb8207..809d0463854 100644 --- a/SU2_PY/SU2/io/historyMap.py +++ b/SU2_PY/SU2/io/historyMap.py @@ -244,6 +244,12 @@ "HEADER": 'bgs[rho*Y_" + ' 'std::to_string(iVar)+"]', "TYPE": "RESIDUAL", }, + 'BGS_SPECIES_" + config->GetChemical_GasComposition(iVar': { + "DESCRIPTION": "BGS residual of " "transported species.", + "GROUP": "BGS_RES", + "HEADER": 'bgs[rho*Y_" + ' 'config->GetChemical_GasComposition(iVar)+"]', + "TYPE": "RESIDUAL", + }, "BGS_TEMPERATURE": { "DESCRIPTION": "Block-Gauss-Seidel residual of the " "temperature", "GROUP": "BGS_RES", @@ -816,6 +822,12 @@ "HEADER": 'max[A_rho*Y_" + ' "std::to_string(iVar) + " '"]', "TYPE": "RESIDUAL", }, + 'MAX_ADJ_SPECIES_" + config->GetChemical_GasComposition(iVar': { + "DESCRIPTION": "Maximum residual " "of the adjoint " "transported " "species.", + "GROUP": "MAX_RES", + "HEADER": 'max[A_rho*Y_" + ' "config->GetChemical_GasComposition(iVar) + " '"]', + "TYPE": "RESIDUAL", + }, "MAX_ADJ_TEMPERATURE": { "DESCRIPTION": "Maximum residual of the temperature.", "GROUP": "MAX_RES", @@ -906,6 +918,12 @@ "HEADER": "max[P]", "TYPE": "RESIDUAL", }, + 'MAX_SPECIES_" + config->GetChemical_GasComposition(iVar': { + "DESCRIPTION": "Maximum residual of " "transported species.", + "GROUP": "MAX_RES", + "HEADER": 'max[rho*Y_" + ' 'config->GetChemical_GasComposition(iVar)+"]', + "TYPE": "RESIDUAL", + }, 'MAX_SPECIES_" + std::to_string(iVar': { "DESCRIPTION": "Maximum residual of " "transported species.", "GROUP": "MAX_RES", @@ -1067,6 +1085,16 @@ "HEADER": 'rms[A_rho*Y_" + ' "std::to_string(iVar) + " '"]', "TYPE": "RESIDUAL", }, + 'RMS_ADJ_SPECIES_" + config->GetChemical_GasComposition(iVar': { + "DESCRIPTION": "Root-mean square " + "residual of the " + "adjoint " + "transported " + "species.", + "GROUP": "RMS_RES", + "HEADER": 'rms[A_rho*Y_" + ' "config->GetChemical_GasComposition(iVar) + " '"]', + "TYPE": "RESIDUAL", + }, "RMS_ADJ_TEMPERATURE": { "DESCRIPTION": "Root-mean square residual of the " "adjoint temperature.", "GROUP": "RMS_RES", @@ -1184,6 +1212,12 @@ "HEADER": 'rms[rho*Y_" + ' 'std::to_string(iVar)+"]', "TYPE": "RESIDUAL", }, + 'RMS_SPECIES_" + config->GetChemical_GasComposition(iVar': { + "DESCRIPTION": "Root-mean square " "residual of " "transported species.", + "GROUP": "RMS_RES", + "HEADER": 'rms[rho*Y_" + ' 'config->GetChemical_GasComposition(iVar)+"]', + "TYPE": "RESIDUAL", + }, "RMS_TEMPERATURE": { "DESCRIPTION": "Root mean square residual of the " "temperature", "GROUP": "RMS_RES", @@ -1367,6 +1401,17 @@ "HEADER": 'Avg_Species_" + ' "std::to_string(iVar", "TYPE": "COEFFICIENT", }, + 'SURFACE_SPECIES_" + config->GetChemical_GasComposition(iVar': { + "DESCRIPTION": "Total average " + 'species " + ' + "std::to_string(iVar) " + '+ " on all ' + "markers set in " + "MARKER_ANALYZE", + "GROUP": "SPECIES_COEFF", + "HEADER": 'Avg_Species_" + ' "config->GetChemical_GasComposition(iVar", + "TYPE": "COEFFICIENT", + }, "SURFACE_SPECIES_VARIANCE": { "DESCRIPTION": "Total species variance", "GROUP": "SPECIES_COEFF", diff --git a/TestCases/species_transport/venturi_primitive_3species/species2_primitiveVenturi_CANTERA.cfg b/TestCases/species_transport/venturi_primitive_3species/species2_primitiveVenturi_CANTERA.cfg new file mode 100644 index 00000000000..befcbaeca56 --- /dev/null +++ b/TestCases/species_transport/venturi_primitive_3species/species2_primitiveVenturi_CANTERA.cfg @@ -0,0 +1,135 @@ +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% % +% SU2 configuration file % +% Case description: Species mixing with 2 species, i.e. 1 transport equations % +% Including CANTERA fluid model. Fluid properties obtained % +% directly from CANTERA. % +% Author: Cristopher Morales Ubal % +% Institution: Eindhoven University of Technology % +% Date: 2026/01/19 % +% File Version 8.4.0 "Harrier" % +% % +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +% ------------- DIRECT, ADJOINT, AND LINEARIZED PROBLEM DEFINITION ------------% +% +SOLVER= INC_RANS +KIND_TURB_MODEL= SST +% +% ---------------- INCOMPRESSIBLE FLOW CONDITION DEFINITION -------------------% +% +INC_DENSITY_MODEL= VARIABLE +INC_DENSITY_INIT= 1.1766 +% +INC_VELOCITY_INIT= ( 1.00, 0.0, 0.0 ) +% +INC_ENERGY_EQUATION= YES +INC_TEMPERATURE_INIT= 300.0 +% +INC_NONDIM= DIMENSIONAL +% +% -------------------- FLUID PROPERTIES ------------------------------------- % +% +FLUID_MODEL= FLUID_CANTERA +COMBUSTION = NO +% +THERMODYNAMIC_PRESSURE= 101325.0 +% +CONDUCTIVITY_MODEL= CANTERA +% +TURBULENT_CONDUCTIVITY_MODEL= NONE +% +GAS_COMPOSITION_NAMES= H2, O2 +CHEMICAL_MECHANISM_FILE= h2o2.yaml +PHASE_NAME= ohmech +TRANSPORT_MODEL= mixture-averaged +% +VISCOSITY_MODEL= CANTERA +% +% -------------------- BOUNDARY CONDITION DEFINITION --------------------------% +% +MARKER_HEATFLUX= ( wall, 0.0 ) +MARKER_SYM= ( axis ) +% +SPECIFIED_INLET_PROFILE= NO +INLET_FILENAME= inlet_venturi.dat +INC_INLET_TYPE= VELOCITY_INLET VELOCITY_INLET +MARKER_INLET= ( gas_inlet, 300, 1.0, 1.0, 0.0, 0.0,\ + air_axial_inlet, 350, 1.0, 0.0, -1.0, 0.0 ) +MARKER_INLET_SPECIES= (gas_inlet, 0.5,\ + air_axial_inlet, 0.6) +% +INC_OUTLET_TYPE= PRESSURE_OUTLET +MARKER_OUTLET= ( outlet, 0.0) +% +% ------------- COMMON PARAMETERS DEFINING THE NUMERICAL METHOD ---------------% +% +NUM_METHOD_GRAD= WEIGHTED_LEAST_SQUARES +% +CFL_NUMBER= 60 +CFL_REDUCTION_SPECIES= 1.0 +CFL_REDUCTION_TURB= 1.0 +% +ITER= 1000 +% +% ------------------------ LINEAR SOLVER DEFINITION ---------------------------% +% +LINEAR_SOLVER= FGMRES +LINEAR_SOLVER_PREC= ILU +LINEAR_SOLVER_ERROR= 1E-8 +LINEAR_SOLVER_ITER= 5 +% +% -------------------- FLOW NUMERICAL METHOD DEFINITION -----------------------% +% +CONV_NUM_METHOD_FLOW= FDS +MUSCL_FLOW= NO +SLOPE_LIMITER_FLOW = NONE +TIME_DISCRE_FLOW= EULER_IMPLICIT +% +% -------------------- SCALAR TRANSPORT ---------------------------------------% +% +KIND_SCALAR_MODEL= SPECIES_TRANSPORT +DIFFUSIVITY_MODEL= CANTERA +DIFFUSIVITY_CONSTANT= 0.001 +% +CONV_NUM_METHOD_SPECIES= BOUNDED_SCALAR +MUSCL_SPECIES= NO +SLOPE_LIMITER_SPECIES = NONE +% +TIME_DISCRE_SPECIES= EULER_IMPLICIT +% +SPECIES_INIT= 0.5 +SPECIES_CLIPPING= NO +SPECIES_CLIPPING_MIN= 0.0 +SPECIES_CLIPPING_MAX= 1.0 +% +% -------------------- TURBULENT TRANSPORT ---------------------------------------% +% +CONV_NUM_METHOD_TURB= BOUNDED_SCALAR +MUSCL_TURB= NO +% +% --------------------------- CONVERGENCE PARAMETERS --------------------------% +% +CONV_FIELD= RMS_PRESSURE, RMS_VELOCITY-X, RMS_VELOCITY-Y, RMS_TKE, RMS_SPECIES +CONV_RESIDUAL_MINVAL= -18 +CONV_STARTITER= 10 +% +% ------------------------- INPUT/OUTPUT INFORMATION --------------------------% +% +MESH_FILENAME= primitiveVenturi.su2 +SCREEN_OUTPUT= INNER_ITER WALL_TIME \ + RMS_PRESSURE RMS_VELOCITY-X RMS_VELOCITY-Y RMS_ENTHALPY RMS_TKE RMS_DISSIPATION RMS_SPECIES_H2 \ + LINSOL_ITER LINSOL_RESIDUAL \ + LINSOL_ITER_TURB LINSOL_RESIDUAL_TURB \ + LINSOL_ITER_SPECIES LINSOL_RESIDUAL_SPECIES SURFACE_SPECIES_VARIANCE +SCREEN_WRT_FREQ_INNER= 10 +% +HISTORY_OUTPUT= RMS_RES FLOW_COEFF LINSOL SPECIES_COEFF SPECIES_COEFF_SURF +MARKER_ANALYZE= outlet gas_inlet air_axial_inlet +% +OUTPUT_FILES= RESTART_ASCII, PARAVIEW_MULTIBLOCK +VOLUME_OUTPUT= RESIDUAL, PRIMITIVE +OUTPUT_WRT_FREQ= 100 +% +RESTART_SOL= NO +SOLUTION_FILENAME= solution diff --git a/UnitTests/SU2_CFD/fluid/CFluidCantera_tests.cpp b/UnitTests/SU2_CFD/fluid/CFluidCantera_tests.cpp new file mode 100644 index 00000000000..0062806e530 --- /dev/null +++ b/UnitTests/SU2_CFD/fluid/CFluidCantera_tests.cpp @@ -0,0 +1,95 @@ +/*! + * \file CFluidCantera_tests.cpp + * \brief Unit tests for Cantera fluid model. + * \author C.Morales Ubal + * \version 8.4.0 "Harrier" + * + * SU2 Project Website: https://su2code.github.io + * + * The SU2 Project is maintained by the SU2 Foundation + * (http://su2foundation.org) + * + * Copyright 2012-2024, SU2 Contributors (cf. AUTHORS.md) + * + * SU2 is free software; you can redistribute it and/or + * modify it under the terms of the GNU Lesser General Public + * License as published by the Free Software Foundation; either + * version 2.1 of the License, or (at your option) any later version. + * + * SU2 is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + * Lesser General Public License for more details. + * + * You should have received a copy of the GNU Lesser General Public + * License along with SU2. If not, see . + */ + +#include "catch.hpp" + +#include + +#if defined(HAVE_CANTERA) +#define USE_CANTERA +#include "../../../SU2_CFD/include/fluid/CFluidCantera.hpp" +#include + +using namespace Cantera; +#endif + +#ifdef USE_CANTERA +TEST_CASE("Fluid_Cantera", "[Multicomponent_flow]") { + /*--- Cantera fluid model unit test cases. ---*/ + + SU2_COMPONENT val_software = SU2_COMPONENT::SU2_CFD; + CConfig* config = new CConfig("multicomponent_cantera.cfg", val_software, true); + CFluidCantera* auxFluidModel = nullptr; + + /*--- Create Cantera fluid model. ---*/ + su2double value_pressure_operating = config->GetPressure_Thermodynamic(); + auxFluidModel = new CFluidCantera(value_pressure_operating, config); + + /*--- Get scalar from config file and set temperature. ---*/ + const su2double* scalar = nullptr; + scalar = config->GetSpecies_Init(); + const su2double Temperature = 300.0; + /*--- Set state using temperature and scalar. ---*/ + auxFluidModel->SetTDState_T(Temperature, scalar); + + /*--- Check values for density and heat capacity. ---*/ + + su2double density = auxFluidModel->GetDensity(); + su2double cp = auxFluidModel->GetCp(); + CHECK(density == Approx(0.924236)); + CHECK(cp == Approx(1277.91)); +} +TEST_CASE("Fluid_Cantera_Combustion", "[Reacting_flow]") { + /*--- Cantera fluid model unit test cases. ---*/ + + SU2_COMPONENT val_software = SU2_COMPONENT::SU2_CFD; + CConfig* config = new CConfig("multicomponent_cantera.cfg", val_software, true); + CFluidCantera* auxFluidModel = nullptr; + + /*--- Create Cantera fluid model. ---*/ + su2double value_pressure_operating = config->GetPressure_Thermodynamic(); + auxFluidModel = new CFluidCantera(value_pressure_operating, config); + + /*--- Get scalar from config file and set temperature. ---*/ + + const su2double* scalar = nullptr; + scalar = config->GetSpecies_Init(); + const su2double Temperature = 1900.0; + /*--- Set state using temperature and scalar ---*/ + auxFluidModel->SetTDState_T(Temperature, scalar); + + /*--- Compute chemical source terms. ---*/ + auxFluidModel->ComputeChemicalSourceTerm(); + + /*--- Check values for source terms. ---*/ + + su2double sourceTerm_H2 = auxFluidModel->GetChemicalSourceTerm(0); + su2double sourceTerm_O2 = auxFluidModel->GetChemicalSourceTerm(1); + CHECK(sourceTerm_H2 == Approx(-0.13633797171426)); + CHECK(sourceTerm_O2 == Approx(-2.16321066087493)); +} +#endif diff --git a/UnitTests/SU2_CFD/fluid/multicomponent_cantera.cfg b/UnitTests/SU2_CFD/fluid/multicomponent_cantera.cfg new file mode 100644 index 00000000000..dc4a9dc4c3b --- /dev/null +++ b/UnitTests/SU2_CFD/fluid/multicomponent_cantera.cfg @@ -0,0 +1,41 @@ +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% % +% SU2 configuration file % +% Case description: Unit test Multicomponent Cantera 3 species H2, O2, N2 % +% Author: Cristopher Morales Ubal % +% Institution: Eindhoven University of Technology % +% Date: 2024/10/01 % +% File Version 8.1.0 "Harrier" % +% % +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +% ------------- DIRECT, ADJOINT, AND LINEARIZED PROBLEM DEFINITION ------------% +% +SOLVER= INC_NAVIER_STOKES +% +% ---------------- INCOMPRESSIBLE FLOW CONDITION DEFINITION -------------------% +% +INC_DENSITY_MODEL= VARIABLE +% +INC_NONDIM= DIMENSIONAL +% +% -------------------- FLUID PROPERTIES ------------------------------------- % +% +FLUID_MODEL= FLUID_CANTERA +COMBUSTION= YES +FLAME_INIT_METHOD= SPARK +% +CONDUCTIVITY_MODEL= CANTERA +VISCOSITY_MODEL= CANTERA +% +GAS_COMPOSITION_NAMES= H2, O2, N2 +CHEMICAL_MECHANISM_FILE= h2o2.yaml +PHASE_NAME= ohmech +TRANSPORT_MODEL= mixture-averaged +% +% -------------------- SCALAR TRANSPORT ---------------------------------------% +% +KIND_SCALAR_MODEL= SPECIES_TRANSPORT +DIFFUSIVITY_MODEL= CANTERA +% +SPECIES_INIT= 0.020137986304114414,0.22830757006769173 diff --git a/UnitTests/meson.build b/UnitTests/meson.build index f8c22511b5f..a6df92c8c52 100644 --- a/UnitTests/meson.build +++ b/UnitTests/meson.build @@ -12,6 +12,7 @@ su2_cfd_tests = files(['Common/geometry/primal_grid/CPrimalGrid_tests.cpp', 'Common/vectorization.cpp', 'Common/toolboxes/ndflattener_tests.cpp', 'Common/containers/CLookupTable_tests.cpp', + 'SU2_CFD/fluid/CFluidCantera_tests.cpp', 'Common/toolboxes/multilayer_perceptron/CLookUp_ANN_tests.cpp', 'SU2_CFD/numerics/CNumerics_tests.cpp', 'SU2_CFD/fluid/CFluidModel_tests.cpp', @@ -32,14 +33,19 @@ su2_cfd_tests_dd = files(['Common/simple_directdiff_test.cpp']) # End of unit test listings # ------------------------------------------------------------------------- + if get_option('enable-tests') if get_option('enable-normal') unit_test_files = su2_cfd_tests + files(['test_driver.cpp']) + dependencies_default = [su2_cfd_dep, common_dep, su2_deps, catch2_dep] + if get_option('enable-cantera') + dependencies_default += cantera_dep + endif test_driver = executable( 'test_driver', unit_test_files, install : true, - dependencies : [su2_cfd_dep, common_dep, su2_deps, catch2_dep], + dependencies : dependencies_default, cpp_args: ['-fPIC', default_warning_flags, su2_cpp_args] ) test('Catch2 test driver', test_driver) diff --git a/meson.build b/meson.build index 78379a45507..a7c790f9647 100644 --- a/meson.build +++ b/meson.build @@ -8,7 +8,6 @@ project('SU2', 'c', 'cpp', 'b_ndebug=true', 'c_std=c99', 'cpp_std=c++17']) - fsmod = import('fs') if not fsmod.exists('su2preconfig.timestamp') error('SU2 must be (pre-)configured with the extended Meson script (./meson.py) or with the ./preconfigure.py script in the SU2 root directory.') @@ -344,6 +343,49 @@ if get_option('enable-coolprop') endif endif +if get_option('enable-cantera') + # Cantera is built from source via its own SCons build, driven silently as part of + # the normal SU2 build (custom_target below), rather than requiring a separate manual + # build step or a prebuilt system/conda package. Building from source (instead of + # linking a prebuilt shared library) keeps the door open for a future AD-enabled + # Cantera build variant compiled against SU2's own scalar type. + cantera_source_dir = join_paths(meson.current_source_dir(), 'subprojects', 'cantera') + cantera_prefix = join_paths(meson.current_build_dir(), 'cantera_install') + cantera_build_script = join_paths(meson.current_source_dir(), 'meson_scripts', 'build_cantera.py') + + cantera_build_command = [python, cantera_build_script, + '--source', cantera_source_dir, + '--prefix', cantera_prefix, + '--stamp', '@OUTPUT@', + '--eigen-inc-dir', join_paths(meson.current_source_dir(), 'externals', 'eigen')] + if get_option('cantera-boost-inc-dir') != '' + cantera_build_command += ['--boost-inc-dir', get_option('cantera-boost-inc-dir')] + endif + if get_option('cpu-arch') != '' + cantera_build_command += ['--cpu-arch', get_option('cpu-arch')] + endif + + cantera_build = custom_target('build_cantera', + output: 'cantera_build.stamp', + command: cantera_build_command, + build_always_stale: true, + console: true, + ) + + # build_cantera.py makes Cantera's installed headers resolve Eigen through an + # "eigen3/Eigen/..." shim it creates next to the install prefix, pointing at SU2's + # own Eigen copy (see the ODR note in that script). Anything compiling those headers - + # including SU2's own CFluidCantera.cpp - needs that same shim directory on its + # include path so the "eigen3/..." form resolves for it too. + cantera_dep = declare_dependency( + compile_args: ['-I' + join_paths(cantera_prefix, 'include'), + '-I' + cantera_prefix + '_eigen_shim'], + link_args: ['-L' + join_paths(cantera_prefix, 'lib'), '-lcantera'], + sources: [cantera_build], + ) + su2_cpp_args +='-DHAVE_CANTERA' +endif + if get_option('enable-mlpcpp') su2_cpp_args += '-DHAVE_MLPCPP' endif @@ -400,6 +442,7 @@ message('''--------------------------------------------------------------------- CoolProp: @12@ MLPCpp: @13@ Tracy Profiler: @14@ + Cantera: @16@ Please be sure to add the $SU2_HOME and $SU2_RUN environment variables, and update your $PATH (and $PYTHONPATH if applicable) with $SU2_RUN @@ -415,7 +458,8 @@ message('''--------------------------------------------------------------------- '''.format(get_option('prefix')+'/bin', meson.project_source_root(), get_option('enable-tecio'), get_option('enable-cgns'), get_option('enable-autodiff'), get_option('enable-directdiff'), get_option('enable-pywrapper'), get_option('enable-mkl'), get_option('enable-openblas'), get_option('enable-pastix'), get_option('enable-mixedprec'), get_option('enable-librom'), get_option('enable-coolprop'), - get_option('enable-mlpcpp'), get_option('tracy_enable'), meson.project_build_root().startswith(meson.project_source_root()) ? meson.project_build_root().split('/')[-1] : meson.project_build_root())) + get_option('enable-mlpcpp'), get_option('tracy_enable'), meson.project_build_root().startswith(meson.project_source_root()) ? meson.project_build_root().split('/')[-1] : meson.project_build_root(), + get_option('enable-cantera'))) if get_option('enable-mpp') if get_option('install-mpp') @@ -435,6 +479,17 @@ if get_option('enable-mpp') endif endif +if get_option('enable-cantera') + message(''' To run SU2 with Cantera, add this line to your .bashrc file: + + export CANTERA_DATA=@0@/data + + Cantera's own compiled-in default data directory should normally be + found automatically, but setting CANTERA_DATA explicitly avoids + relying on that at runtime. + '''.format(cantera_prefix)) +endif + if get_option('enable-librom') message(''' To run SU2 with libROM library, add this line to your .bashrc file: diff --git a/meson_options.txt b/meson_options.txt index 6adad19386b..7a6cb599d0f 100644 --- a/meson_options.txt +++ b/meson_options.txt @@ -20,7 +20,9 @@ option('enable-singleprec', type : 'boolean', value : false, description: 'use s option('extra-deps', type : 'string', value : '', description: 'comma-separated list of extra (custom) dependencies to add for compilation') option('enable-mpp', type : 'boolean', value : false, description: 'enable Mutation++ support') option('install-mpp', type : 'boolean', value : false, description: 'install Mutation++ in the directory defined with --prefix') -option('enable-coolprop', type : 'boolean', value : false, description: 'enable CoolProp support') +option('enable-coolprop', type : 'boolean', value : false, description: 'enable CoolProp support') +option('enable-cantera', type : 'boolean', value : false, description: 'enable Cantera support') +option('cantera-boost-inc-dir', type : 'string', value : '', description: 'explicit Boost header directory for building Cantera from source (only needed if SCons cannot auto-detect Boost)') option('enable-mlpcpp', type : 'boolean', value : false, description: 'enable MLPCpp support') option('enable-gprof', type : 'boolean', value : false, description: 'enable profiling through gprof') option('opdi-backend', type : 'combo', choices : ['auto', 'macro', 'ompt'], value : 'auto', description: 'OpDiLib backend choice') diff --git a/meson_scripts/build_cantera.py b/meson_scripts/build_cantera.py new file mode 100644 index 00000000000..c64ecc16916 --- /dev/null +++ b/meson_scripts/build_cantera.py @@ -0,0 +1,153 @@ +#!/usr/bin/env python3 + +## \file build_cantera.py +# \brief Builds and installs Cantera as a static C++ library via its own +# SCons build system, invoked automatically from meson.build so that +# '-Denable-cantera=true' does not require a separate manual build step. +# \author SU2 Contributors +# \version 8.5.0 "Harrier" +# +# SU2 Project Website: https://su2code.github.io +# +# The SU2 Project is maintained by the SU2 Foundation +# (http://su2foundation.org) +# +# Copyright 2012-2026, SU2 Contributors (cf. AUTHORS.md) +# +# SU2 is free software; you can redistribute it and/or +# modify it under the terms of the GNU Lesser General Public +# License as published by the Free Software Foundation; either +# version 2.1 of the License, or (at your option) any later version. +# +# SU2 is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU +# Lesser General Public License for more details. +# +# You should have received a copy of the GNU Lesser General Public +# License along with SU2. If not, see . + +import argparse +import multiprocessing +import os +import shutil +import subprocess +import sys + +# Only the pieces of Cantera SU2 actually links against (thermo/kinetics/transport) +# are built: no Python/Fortran interfaces, docs, or test suite. Cantera's own SCons +# script fetches its required submodules (fmt, yaml-cpp, sundials) on its own. +# +# Eigen is the one exception: SU2 already vendors its own Eigen copy (externals/eigen) +# and links Eigen-templated code from both SU2 and Cantera into the same executable. +# Letting Cantera vendor a second, different Eigen version is an ODR violation - the +# linker silently merges identically-mangled template instantiations from the two +# versions, which segfaults deep inside Eigen at runtime. Building against SU2's own +# Eigen (system_eigen=y + extra_inc_dirs below) keeps a single, consistent copy. +# +# Same version isn't sufficient on its own, though: Eigen's reduction/QR templates are +# header-only, so each translation unit compiles its own body, and the linker COMDAT-folds +# identically-mangled instantiations from different translation units into one - even +# across different -march targets. SU2 compiles with -march= (native by +# default); if Cantera's own compile used a different (or no) -march, the linker could +# keep either compiled body for a shared symbol, mixing SIMD-width assumptions and +# segfaulting. --cpu-arch keeps the two in sync (same fix already applied to CoolProp's +# CMAKE_CXX_FLAGS in meson.build). +CANTERA_SCONS_ARGS = [ + "layout=compact", + "python_package=n", + "f90_interface=n", + "googletest=none", + "hdf_support=n", + "doxygen_docs=no", + "sphinx_docs=no", + "system_eigen=y", + "system_fmt=n", + "system_yamlcpp=n", + "system_sundials=n", + "system_highfive=n", + "logging=warning", +] + + +def main(): + parser = argparse.ArgumentParser(description=__doc__) + parser.add_argument( + "--source", required=True, help="Path to the Cantera submodule checkout" + ) + parser.add_argument( + "--prefix", required=True, help="Install prefix for the built library" + ) + parser.add_argument( + "--stamp", + required=True, + help="File to create on a successful build, for meson tracking", + ) + parser.add_argument( + "--boost-inc-dir", default="", help="Explicit Boost header directory (optional)" + ) + parser.add_argument( + "--eigen-inc-dir", + required=True, + help="SU2's own Eigen include directory (externals/eigen)", + ) + parser.add_argument( + "--cpu-arch", + default="", + help="Value of SU2's -march= target (cpu-arch option), kept in sync so Eigen template code Cantera and SU2 both instantiate is not ODR-mixed across incompatible vectorization ISAs", + ) + args = parser.parse_args() + + if not os.path.exists(os.path.join(args.source, "SConstruct")): + sys.exit( + "Cantera sources not found in '{}'.\n" + "Run './preconfigure.py' (or 'git submodule update --init subprojects/cantera') " + "before configuring with -Denable-cantera=true.".format(args.source) + ) + + scons = shutil.which("scons") + if scons is None: + sys.exit( + "scons was not found on the PATH. Install it (e.g. 'pip install scons' or " + "'apt install scons') to build Cantera from source." + ) + + # Cantera's SCons script probes for Eigen by searching for "eigen3/Eigen/Dense" first, + # then plain "Eigen/Dense", across the compiler's whole include search path - not just + # extra_inc_dirs. If the host has a system libeigen3-dev package, that "eigen3/..." + # form gets found ahead of SU2's own (unprefixed) copy, defeating the point of pointing + # this build at SU2's Eigen at all. Force the match onto SU2's copy by presenting it + # under an "eigen3/" name via a symlink, in a directory searched before system ones. + eigen_shim_dir = args.prefix.rstrip("/") + "_eigen_shim" + os.makedirs(eigen_shim_dir, exist_ok=True) + eigen_shim_link = os.path.join(eigen_shim_dir, "eigen3") + if os.path.islink(eigen_shim_link) or os.path.exists(eigen_shim_link): + os.remove(eigen_shim_link) + os.symlink(args.eigen_inc_dir, eigen_shim_link) + + scons_command = ( + [scons, "build", "install", "prefix=" + args.prefix] + + CANTERA_SCONS_ARGS + + ["extra_inc_dirs=" + eigen_shim_dir] + + ["-j{}".format(multiprocessing.cpu_count())] + ) + if args.boost_inc_dir: + scons_command.append("boost_inc_dir=" + args.boost_inc_dir) + if args.cpu_arch: + scons_command.append("cc_flags=-march=" + args.cpu_arch) + + result = subprocess.run(scons_command, cwd=args.source) + if result.returncode != 0: + sys.exit( + "Cantera build failed (see output above). If SCons could not find Boost headers, " + "set -Dcantera-boost-inc-dir= and reconfigure, " + "or install Boost development headers (e.g. the 'libboost-dev' package)." + ) + + # meson tracks this custom_target as always-stale, so re-running scons here is cheap + # once Cantera itself is up to date; the stamp file just gives meson an output to track. + open(args.stamp, "w").close() + + +if __name__ == "__main__": + main() diff --git a/meson_scripts/init.py b/meson_scripts/init.py index 7a9ce0f0e88..a6d4ec22531 100755 --- a/meson_scripts/init.py +++ b/meson_scripts/init.py @@ -47,6 +47,7 @@ def init_submodules( own_opdi=True, own_mpp=True, own_cool=True, + own_cantera=True, own_mel=True, own_fado=True, own_mlpcpp=True, @@ -70,6 +71,8 @@ def init_submodules( github_repo_mpp = "https://github.com/mutationpp/Mutationpp" sha_version_coolprop = "98b3523d5daa98454618d381d2ae53f7471d216b" github_repo_coolprop = "https://github.com/CoolProp/CoolProp" + sha_version_cantera = "f22643ecb91f7e8d1852197658efc5260d54d4a3" + github_repo_cantera = "https://github.com/cantera/cantera" sha_version_mel = "46205ab019e5224559091375a6d71aabae6bc5b9" github_repo_mel = "https://github.com/pcarruscag/MEL" sha_version_fado = "ce7ee018e4e699af5028d69baa1939fea290e18a" @@ -91,6 +94,7 @@ def init_submodules( ninja_name = "ninja" mpp_name = "Mutationpp" coolprop_name = "CoolProp" + cantera_name = "cantera" mel_name = "MEL" fado_name = "FADO" mlpcpp_name = "MLPCpp" @@ -107,6 +111,7 @@ def init_submodules( alt_name_eigen = base_path + "eigen" alt_name_mpp = cur_dir + os.path.sep + "subprojects" + os.path.sep + "Mutationpp" alt_name_coolprop = cur_dir + os.path.sep + "subprojects" + os.path.sep + "CoolProp" + alt_name_cantera = cur_dir + os.path.sep + "subprojects" + os.path.sep + "cantera" alt_name_mlpcpp = cur_dir + os.path.sep + "subprojects" + os.path.sep + "MLPCpp" if method == "auto": @@ -135,6 +140,8 @@ def init_submodules( submodule_status(alt_name_mpp, sha_version_mpp) if own_cool: submodule_status(alt_name_coolprop, sha_version_coolprop) + if own_cantera: + submodule_status(alt_name_cantera, sha_version_cantera) if own_mel: submodule_status(alt_name_mel, sha_version_mel) if own_fado: @@ -173,6 +180,13 @@ def init_submodules( github_repo_coolprop, sha_version_coolprop, ) + if own_cantera: + download_module( + cantera_name, + alt_name_cantera, + github_repo_cantera, + sha_version_cantera, + ) if own_mel: download_module(mel_name, alt_name_mel, github_repo_mel, sha_version_mel) if own_fado: @@ -261,6 +275,20 @@ def submodule_status(path, sha_commit): print(original_path) os.chdir(original_path) print("CoolProp updated") + if "cantera" in path: + # update cantera + original_path = os.getcwd() + print("update cantera") + absolute_path = sys.path[0] + relative_path = "subprojects/cantera" + full_path = os.path.join(absolute_path, relative_path) + os.chdir(full_path) + print(full_path) + subprocess.run(["git", "submodule", "init"]) + subprocess.run(["git", "submodule", "update"]) + print(original_path) + os.chdir(original_path) + print("cantera updated") # Check that the SHA tag stored in this file matches the one stored in the git index cur_sha_commit = status[1:].split(" ")[0] if cur_sha_commit != sha_commit: diff --git a/preconfigure.py b/preconfigure.py index 16744a48bf0..715d1a72b1f 100755 --- a/preconfigure.py +++ b/preconfigure.py @@ -89,6 +89,7 @@ def run( own_opdi=True, own_mpp=True, own_cool=True, + own_cantera=True, own_mel=True, own_fado=True, own_mlpcpp=True, @@ -103,6 +104,7 @@ def run( own_opdi=own_opdi, own_mpp=own_mpp, own_cool=own_cool, + own_cantera=own_cantera, own_mel=own_mel, own_fado=own_fado, own_mlpcpp=own_mlpcpp, @@ -146,6 +148,11 @@ def run( help="do not download own copy of CoolProp", action="store_false", ) + parser.add_argument( + "--no-cantera", + help="do not download own copy of Cantera", + action="store_false", + ) parser.add_argument( "--no-mel", help="do not download own copy of MEL", action="store_false" ) @@ -169,6 +176,7 @@ def run( own_opdi=args.no_opdi, own_mpp=args.no_mpp, own_cool=args.no_coolprop, + own_cantera=args.no_cantera, own_mel=args.no_mel, own_fado=args.no_fado, own_mlpcpp=args.no_mlpcpp, diff --git a/subprojects/cantera b/subprojects/cantera new file mode 160000 index 00000000000..f22643ecb91 --- /dev/null +++ b/subprojects/cantera @@ -0,0 +1 @@ +Subproject commit f22643ecb91f7e8d1852197658efc5260d54d4a3