diff --git a/Common/src/CConfig.cpp b/Common/src/CConfig.cpp index 557487dfbdfe..c5984d230f55 100644 --- a/Common/src/CConfig.cpp +++ b/Common/src/CConfig.cpp @@ -3601,14 +3601,6 @@ void CConfig::SetPostprocessing(SU2_COMPONENT val_software, unsigned short val_i SU2_MPI::Error("AUSMPW+ is extremely unstable. Feel free to fix me!", CURRENT_FUNCTION); } - if (GetGasModel() == "ARGON" && GetKind_FluidModel() == SU2_NONEQ){ - SU2_MPI::Error("ARGON is not working with SU2_NONEQ fluid model!", CURRENT_FUNCTION); - } - - if (GetKind_FluidModel() == MUTATIONPP && GetFrozen() == true){ - SU2_MPI::Error("The option of FROZEN_MIXTURE is not yet working with Mutation++ support.", CURRENT_FUNCTION); - } - if(GetBoolTurbomachinery()){ nBlades = new su2double[nZone]; FreeStreamTurboNormal= new su2double[3]; @@ -5042,7 +5034,8 @@ void CConfig::SetPostprocessing(SU2_COMPONENT val_software, unsigned short val_i /*--- Specifying a deforming surface requires a mesh deformation solver. ---*/ if (GetSurface_Movement(DEFORMING)) Deform_Mesh = true; - if (GetGasModel() == "ARGON") monoatomic = true; + if (GetGasModel() == "ARGON") {monoatomic = true;} + else {monoatomic = false;} // This option is deprecated. After a grace period until 7.2.0 the usage warning should become an error. if(OptionIsSet("CONV_CRITERIA") && rank == MASTER_NODE) { diff --git a/SU2_CFD/include/fluid/CMutationTCLib.hpp b/SU2_CFD/include/fluid/CMutationTCLib.hpp index 08d08058413f..77e80c15eaa3 100644 --- a/SU2_CFD/include/fluid/CMutationTCLib.hpp +++ b/SU2_CFD/include/fluid/CMutationTCLib.hpp @@ -87,7 +87,7 @@ class CMutationTCLib : public CNEMOGas { /*! * \brief Compute vector of species V-E energy. */ - vector& ComputeSpeciesEve(su2double val_T) final; + vector& ComputeSpeciesEve(su2double val_T, bool vibe_only = false) final; /*! * \brief Compute species net production rates. diff --git a/SU2_CFD/include/fluid/CNEMOGas.hpp b/SU2_CFD/include/fluid/CNEMOGas.hpp index 21fc4cd990c0..86b812f42ae8 100644 --- a/SU2_CFD/include/fluid/CNEMOGas.hpp +++ b/SU2_CFD/include/fluid/CNEMOGas.hpp @@ -134,7 +134,7 @@ class CNEMOGas : public CFluidModel { /*! * \brief Compute vector of species V-E energy. */ - virtual vector& ComputeSpeciesEve(su2double val_T) = 0; + virtual vector& ComputeSpeciesEve(su2double val_T, bool vibe_only = false) = 0; /*! * \brief Compute species enthalpies. diff --git a/SU2_CFD/include/fluid/CSU2TCLib.hpp b/SU2_CFD/include/fluid/CSU2TCLib.hpp index fc1c87e73585..5fcbe99bd545 100644 --- a/SU2_CFD/include/fluid/CSU2TCLib.hpp +++ b/SU2_CFD/include/fluid/CSU2TCLib.hpp @@ -50,7 +50,7 @@ class CSU2TCLib : public CNEMOGas { ArrheniusEta, /*!< \brief Arrhenius reaction temperature exponent */ ArrheniusTheta, /*!< \brief Arrhenius reaction characteristic temperature */ CharVibTemp, /*!< \brief Characteristic vibrational temperature for e_vib */ - RotationModes, /*!< \brief Rotational modes of energy storage */ + RotationModes, /*!< \brief Rotational modes of energy storage */ Tcf_a, /*!< \brief Rate controlling temperature exponent (fwd) */ Tcf_b, /*!< \brief Rate controlling temperature exponent (fwd) */ Tcb_a, /*!< \brief Rate controlling temperature exponent (bkw) */ @@ -115,7 +115,7 @@ class CSU2TCLib : public CNEMOGas { /*! * \brief Compute species V-E energy. */ - vector& ComputeSpeciesEve(su2double val_T) final; + vector& ComputeSpeciesEve(su2double val_T, bool vibe_only = false) final; /*! * \brief Compute species net production rates. diff --git a/SU2_CFD/include/solvers/CFVMFlowSolverBase.hpp b/SU2_CFD/include/solvers/CFVMFlowSolverBase.hpp index 089b11fb8522..f5b569d417f0 100644 --- a/SU2_CFD/include/solvers/CFVMFlowSolverBase.hpp +++ b/SU2_CFD/include/solvers/CFVMFlowSolverBase.hpp @@ -277,6 +277,11 @@ class CFVMFlowSolverBase : public CSolver { */ void SumEdgeFluxes(const CGeometry* geometry); + /*! + * \brief Computes and sets the required auxilliary vars (and gradients) for axisymmetric flow. + */ + void ComputeAxisymmetricAuxGradients(CGeometry *geometry, const CConfig* config); + /*! * \brief Instantiate a SIMD numerics object. */ diff --git a/SU2_CFD/include/solvers/CFVMFlowSolverBase.inl b/SU2_CFD/include/solvers/CFVMFlowSolverBase.inl index 256643fe1d13..874d57b0c457 100644 --- a/SU2_CFD/include/solvers/CFVMFlowSolverBase.inl +++ b/SU2_CFD/include/solvers/CFVMFlowSolverBase.inl @@ -2923,3 +2923,32 @@ su2double CFVMFlowSolverBase::EvaluateCommonObjFunc(const CConfig& config) return objFun; } + +template +void CFVMFlowSolverBase::ComputeAxisymmetricAuxGradients(CGeometry *geometry, const CConfig* config) { + + /*--- Loop through all points to set the auxvargrad --*/ + SU2_OMP_FOR_STAT(omp_chunk_size) + for (auto iPoint = 0ul; iPoint < nPoint; iPoint++) { + su2double yCoord = geometry->nodes->GetCoord(iPoint, 1); + su2double yVelocity = nodes->GetVelocity(iPoint,1); + su2double xVelocity = nodes->GetVelocity(iPoint,0); + su2double Total_Viscosity = nodes->GetLaminarViscosity(iPoint) + nodes->GetEddyViscosity(iPoint); + + if (yCoord > EPS){ + su2double nu_v_on_y = Total_Viscosity*yVelocity/yCoord; + nodes->SetAuxVar(iPoint, 0, nu_v_on_y); + nodes->SetAuxVar(iPoint, 1, nu_v_on_y*yVelocity); + nodes->SetAuxVar(iPoint, 2, nu_v_on_y*xVelocity); + } + } + END_SU2_OMP_FOR + + /*--- Compute the auxiliary variable gradient with GG or WLS. ---*/ + if (config->GetKind_Gradient_Method() == GREEN_GAUSS) { + SetAuxVar_Gradient_GG(geometry, config); + } + if (config->GetKind_Gradient_Method() == WEIGHTED_LEAST_SQUARES) { + SetAuxVar_Gradient_LS(geometry, config); + } +} diff --git a/SU2_CFD/src/fluid/CMutationTCLib.cpp b/SU2_CFD/src/fluid/CMutationTCLib.cpp index c1012f14f40c..1d584315b5ee 100644 --- a/SU2_CFD/src/fluid/CMutationTCLib.cpp +++ b/SU2_CFD/src/fluid/CMutationTCLib.cpp @@ -48,6 +48,7 @@ CMutationTCLib::CMutationTCLib(const CConfig* config, unsigned short val_nDim): transport_model = "Gupta-Yos"; opt.setStateModel("ChemNonEqTTv"); + if (frozen) opt.setMechanism("none"); opt.setViscosityAlgorithm(transport_model); opt.setThermalConductivityAlgorithm(transport_model); @@ -118,7 +119,7 @@ vector& CMutationTCLib::ComputeMixtureEnergies(){ return energies; } -vector& CMutationTCLib::ComputeSpeciesEve(su2double val_T){ +vector& CMutationTCLib::ComputeSpeciesEve(su2double val_T, bool vibe_only){ SetTDStateRhosTTv(rhos, T, val_T); diff --git a/SU2_CFD/src/fluid/CNEMOGas.cpp b/SU2_CFD/src/fluid/CNEMOGas.cpp index 1dddedfdf22b..6e1a864c3797 100644 --- a/SU2_CFD/src/fluid/CNEMOGas.cpp +++ b/SU2_CFD/src/fluid/CNEMOGas.cpp @@ -40,7 +40,7 @@ CNEMOGas::CNEMOGas(const CConfig* config, unsigned short val_nDim): CFluidModel( Cvtrs.resize(nSpecies,0.0); Cvves.resize(nSpecies,0.0); eves.resize(nSpecies,0.0); - hs.resize(nSpecies,0.0); + hs.resize(2*nSpecies,0.0); ws.resize(nSpecies,0.0); DiffusionCoeff.resize(nSpecies,0.0); Enthalpy_Formation.resize(nSpecies,0.0); @@ -59,23 +59,21 @@ CNEMOGas::CNEMOGas(const CConfig* config, unsigned short val_nDim): CFluidModel( void CNEMOGas::SetTDStatePTTv(su2double val_pressure, const su2double *val_massfrac, su2double val_temperature, su2double val_temperature_ve){ - su2double denom; - - for (iSpecies = 0; iSpecies < nHeavy; iSpecies++) + for (iSpecies = 0; iSpecies < nSpecies; iSpecies++) MassFrac[iSpecies] = val_massfrac[iSpecies]; Pressure = val_pressure; T = val_temperature; Tve = val_temperature_ve; - denom = 0.0; + su2double denom = 0.0; /*--- Calculate mixture density from supplied primitive quantities ---*/ - for (iSpecies = 0; iSpecies < nHeavy; iSpecies++) - denom += MassFrac[iSpecies] * (Ru/MolarMass[iSpecies]) * T; for (iSpecies = 0; iSpecies < nEl; iSpecies++) - denom += MassFrac[nSpecies-1] * (Ru/MolarMass[nSpecies-1]) * Tve; + denom += MassFrac[iSpecies] * (Ru/MolarMass[iSpecies]) * Tve; + for (iSpecies = nEl; iSpecies < nSpecies; iSpecies++) + denom += MassFrac[iSpecies] * (Ru/MolarMass[iSpecies]) * T; Density = Pressure / denom; - + for (iSpecies = 0; iSpecies < nSpecies; iSpecies++){ rhos[iSpecies] = MassFrac[iSpecies]*Density; MassFrac[iSpecies] = rhos[iSpecies]/Density; @@ -84,10 +82,8 @@ void CNEMOGas::SetTDStatePTTv(su2double val_pressure, const su2double *val_massf su2double CNEMOGas::ComputeSoundSpeed(){ - su2double conc, rhoCvtr; - - conc = 0.0; - rhoCvtr = 0.0; + su2double conc = 0.0; + su2double rhoCvtr = 0.0; Density = 0.0; auto& Cvtrs = GetSpeciesCvTraRot(); @@ -95,7 +91,7 @@ su2double CNEMOGas::ComputeSoundSpeed(){ for (iSpecies = 0; iSpecies < nSpecies; iSpecies++) Density+=rhos[iSpecies]; - for (iSpecies = 0; iSpecies < nHeavy; iSpecies++){ + for (iSpecies = nEl; iSpecies < nSpecies; iSpecies++){ conc += rhos[iSpecies]/MolarMass[iSpecies]; rhoCvtr += rhos[iSpecies] * Cvtrs[iSpecies]; } @@ -108,10 +104,10 @@ su2double CNEMOGas::ComputePressure(){ su2double P = 0.0; - for (iSpecies = 0; iSpecies < nHeavy; iSpecies++) - P += rhos[iSpecies] * Ru/MolarMass[iSpecies] * T; for (iSpecies = 0; iSpecies < nEl; iSpecies++) - P += rhos[nSpecies-1] * Ru/MolarMass[nSpecies-1] * Tve; + P += rhos[iSpecies] * Ru/MolarMass[iSpecies] * Tve; + for (iSpecies = nEl; iSpecies < nSpecies; iSpecies++) + P += rhos[iSpecies] * Ru/MolarMass[iSpecies] * T; Pressure = P; @@ -124,7 +120,7 @@ su2double CNEMOGas::ComputeGasConstant(){ su2double Mass = 0.0; // TODO - extend for ionization - for (iSpecies = 0; iSpecies < nHeavy; iSpecies++) + for (iSpecies = 0; iSpecies < nSpecies; iSpecies++) Mass += MassFrac[iSpecies] * MolarMass[iSpecies]; GasConstant = Ru / Mass; @@ -163,15 +159,13 @@ void CNEMOGas::ComputedPdU(su2double *V, vector& val_eves, su2double // Note: Electron energy not included properly. - su2double CvtrBAR, rhoCvtr, rhoCvve, rho_el, sqvel, conc, ef; - if (val_dPdU == NULL) { SU2_MPI::Error("Array dPdU not allocated!", CURRENT_FUNCTION); } /*--- Determine the electron density (if ionized) ---*/ - if (ionization) { rho_el = rhos[nSpecies-1]; } - else { rho_el = 0.0; } + su2double rho_el = 0.0; + if (ionization) { rho_el = rhos[0]; } /*--- Necessary indexes to assess primitive variables ---*/ unsigned long RHOS_INDEX = 0; @@ -189,14 +183,14 @@ void CNEMOGas::ComputedPdU(su2double *V, vector& val_eves, su2double Ref_Temperature = GetRefTemperature(); /*--- Rename for convenience ---*/ - rhoCvtr = V[RHOCVTR_INDEX]; - rhoCvve = V[RHOCVVE_INDEX]; + su2double rhoCvtr = V[RHOCVTR_INDEX]; + su2double rhoCvve = V[RHOCVVE_INDEX]; T = V[T_INDEX]; /*--- Pre-compute useful quantities ---*/ - CvtrBAR = 0.0; - sqvel = 0.0; - conc = 0.0; + su2double CvtrBAR = 0.0; + su2double sqvel = 0.0; + su2double conc = 0.0; for (iDim = 0; iDim < nDim; iDim++) sqvel += V[VEL_INDEX+iDim] * V[VEL_INDEX+iDim]; for (iSpecies = 0; iSpecies < nSpecies; iSpecies++) { @@ -205,15 +199,16 @@ void CNEMOGas::ComputedPdU(su2double *V, vector& val_eves, su2double } /*--- Species density derivatives ---*/ - for (iSpecies = 0; iSpecies < nHeavy; iSpecies++) { - ef = Enthalpy_Formation[iSpecies] - Ru/MolarMass[iSpecies]*Ref_Temperature[iSpecies]; + su2double ef = 0.0; + for (iSpecies = nEl; iSpecies < nHeavy; iSpecies++) { + ef = Enthalpy_Formation[iSpecies] - Ru/MolarMass[iSpecies]*Ref_Temperature[iSpecies]; val_dPdU[iSpecies] = T*Ru/MolarMass[iSpecies] + Ru*conc/rhoCvtr * (-Cvtrs[iSpecies]*(T-Ref_Temperature[iSpecies]) - ef + 0.5*sqvel); } if (ionization) { - for (iSpecies = 0; iSpecies < nHeavy; iSpecies++) { + for (iSpecies = nEl; iSpecies < nSpecies; iSpecies++) { // evibs = Ru/MolarMass[iSpecies] * thetav[iSpecies]/(exp(thetav[iSpecies]/Tve)-1.0); // num = 0.0; // denom = g[iSpecies][0] * exp(-thetae[iSpecies][0]/Tve); @@ -223,12 +218,11 @@ void CNEMOGas::ComputedPdU(su2double *V, vector& val_eves, su2double // } // eels = Ru/MolarMass[iSpecies] * (num/denom); - val_dPdU[iSpecies] -= rho_el * Ru/MolarMass[nSpecies-1] * (val_eves[iSpecies])/rhoCvve; + val_dPdU[iSpecies] -= rho_el * Ru/MolarMass[0] * (val_eves[iSpecies])/rhoCvve; } - ef = Enthalpy_Formation[nSpecies-1] - Ru/MolarMass[nSpecies-1]*Ref_Temperature[nSpecies-1]; - val_dPdU[nSpecies-1] = Ru*conc/rhoCvtr * (-ef + 0.5*sqvel) - + Ru/MolarMass[nSpecies-1]*Tve - - rho_el*Ru/MolarMass[nSpecies-1] * (-3.0/2.0*Ru/MolarMass[nSpecies-1]*Tve)/rhoCvve; + ef = Enthalpy_Formation[0] - Ru/MolarMass[0]*Ref_Temperature[0]; + val_dPdU[0] = Ru*conc/rhoCvtr * (-ef + 0.5*sqvel) + Ru/MolarMass[0]*Tve + - rho_el*Ru/MolarMass[0] * (-3.0/2.0*Ru/MolarMass[0]*Tve)/rhoCvve; } /*--- Momentum derivatives ---*/ @@ -240,13 +234,12 @@ void CNEMOGas::ComputedPdU(su2double *V, vector& val_eves, su2double /*--- Vib.-el energy derivative ---*/ val_dPdU[nSpecies+nDim+1] = -val_dPdU[nSpecies+nDim] + - rho_el*Ru/MolarMass[nSpecies-1]*1.0/rhoCvve; + rho_el*Ru/MolarMass[0]*1.0/rhoCvve; } void CNEMOGas::ComputedTdU(su2double *V, su2double *val_dTdU){ - su2double v2, ef, rhoCvtr; su2double Vel[3] = {0.0}; /*--- Necessary indexes to assess primitive variables ---*/ @@ -255,8 +248,8 @@ void CNEMOGas::ComputedTdU(su2double *V, su2double *val_dTdU){ unsigned long RHOCVTR_INDEX = nSpecies+nDim+6; /*--- Rename for convenience ---*/ - T = V[T_INDEX]; - rhoCvtr = V[RHOCVTR_INDEX]; + T = V[T_INDEX]; + su2double rhoCvtr = V[RHOCVTR_INDEX]; Cvtrs = GetSpeciesCvTraRot(); Enthalpy_Formation = GetSpeciesFormationEnthalpy(); @@ -265,11 +258,11 @@ void CNEMOGas::ComputedTdU(su2double *V, su2double *val_dTdU){ /*--- Calculate supporting quantities ---*/ for (iDim = 0; iDim < nDim; iDim++) Vel[iDim] = V[VEL_INDEX+iDim]*V[VEL_INDEX+iDim]; - v2 = GeometryToolbox::SquaredNorm(nDim,Vel); + su2double v2 = GeometryToolbox::SquaredNorm(nDim,Vel); /*--- Species density derivatives ---*/ - for (iSpecies = 0; iSpecies < nHeavy; iSpecies++) { - ef = Enthalpy_Formation[iSpecies] - Ru/MolarMass[iSpecies]*Ref_Temperature[iSpecies]; + for (iSpecies = nEl; iSpecies < nSpecies; iSpecies++) { + su2double ef = Enthalpy_Formation[iSpecies] - Ru/MolarMass[iSpecies]*Ref_Temperature[iSpecies]; val_dTdU[iSpecies] = (-ef + 0.5*v2 + Cvtrs[iSpecies]*(Ref_Temperature[iSpecies]-T)) / rhoCvtr; } @@ -289,13 +282,11 @@ void CNEMOGas::ComputedTdU(su2double *V, su2double *val_dTdU){ void CNEMOGas::ComputedTvedU(su2double *V, vector& val_eves, su2double *val_dTvedU){ - su2double rhoCvve; - /*--- Necessary indexes to assess primitive variables ---*/ unsigned long RHOCVVE_INDEX = nSpecies+nDim+7; /*--- Rename for convenience ---*/ - rhoCvve = V[RHOCVVE_INDEX]; + su2double rhoCvve = V[RHOCVVE_INDEX]; /*--- Species density derivatives ---*/ for (iSpecies = 0; iSpecies < nSpecies; iSpecies++) { diff --git a/SU2_CFD/src/fluid/CSU2TCLib.cpp b/SU2_CFD/src/fluid/CSU2TCLib.cpp index c4ce389e7c0e..624e1e5912bb 100644 --- a/SU2_CFD/src/fluid/CSU2TCLib.cpp +++ b/SU2_CFD/src/fluid/CSU2TCLib.cpp @@ -58,13 +58,13 @@ CSU2TCLib::CSU2TCLib(const CConfig* config, unsigned short val_nDim, bool viscou if (gas_model =="ARGON"){ if (nSpecies != 1) { - cout << "CONFIG ERROR: nSpecies mismatch between gas model & gas composition" << endl; + SU2_MPI::Error("CONFIG ERROR: nSpecies mismatch between gas model & gas composition", CURRENT_FUNCTION); } mf = 0.0; for (iSpecies = 0; iSpecies < nSpecies; iSpecies++) mf += MassFrac_Freestream[iSpecies]; if (mf != 1.0) { - cout << "CONFIG ERROR: Intial gas mass fractions do not sum to 1!" << " mf is equal to "<< mf <& CSU2TCLib::ComputeSpeciesCvVibEle(){ vector& CSU2TCLib::ComputeMixtureEnergies(){ - su2double rhoEmix, rhoEve, Ef, Ev, Ee, num, denom; + su2double Ev, Ee, Ef, num; - rhoEmix = 0.0; - rhoEve = 0.0; - denom = 0.0; + su2double rhoEmix = 0.0; + su2double rhoEve = 0.0; + su2double denom = 0.0; for (iSpecies = 0; iSpecies < nHeavy; iSpecies++){ @@ -711,7 +711,7 @@ vector& CSU2TCLib::ComputeMixtureEnergies(){ // Species electronic energy num = 0.0; - denom = ElDegeneracy(iSpecies,0) * exp(CharElTemp(iSpecies,0)/Tve); + denom = ElDegeneracy(iSpecies,0) * exp(-CharElTemp(iSpecies,0)/Tve); for (iEl = 1; iEl < nElStates[iSpecies]; iEl++) { num += ElDegeneracy(iSpecies,iEl) * CharElTemp(iSpecies,iEl) * exp(-CharElTemp(iSpecies,iEl)/Tve); denom += ElDegeneracy(iSpecies,iEl) * exp(-CharElTemp(iSpecies,iEl)/Tve); @@ -742,7 +742,7 @@ vector& CSU2TCLib::ComputeMixtureEnergies(){ } -vector& CSU2TCLib::ComputeSpeciesEve(su2double val_T){ +vector& CSU2TCLib::ComputeSpeciesEve(su2double val_T, bool vibe_only){ su2double Ev, Eel, Ef, num, denom; unsigned short iElectron = nSpecies-1; @@ -774,8 +774,8 @@ vector& CSU2TCLib::ComputeSpeciesEve(su2double val_T){ } Eel = Ru/MolarMass[iSpecies] * (num/denom); } - - eves[iSpecies] = Ev + Eel; + if(vibe_only == true) {eves[iSpecies] = Ev;} + else {eves[iSpecies] = Ev + Eel;} } return eves; @@ -785,14 +785,14 @@ vector& CSU2TCLib::ComputeNetProductionRates(){ /*--- Nonequilibrium chemistry ---*/ unsigned short ii, iReaction; - su2double T_min, epsilon, Thf, Thb, Trxnf, Trxnb, Keq, kf, kb, kfb, fwdRxn, bkwRxn, af, bf, ab, bb; + su2double Thf, Thb, Trxnf, Trxnb, Keq, kf, kb, kfb, fwdRxn, bkwRxn, af, bf, ab, bb; /*--- Define artificial chemistry parameters ---*/ // Note: These parameters artificially increase the rate-controlling reaction // temperature. This relaxes some of the stiffness in the chemistry // source term. - T_min = 800.0; - epsilon = 80; + su2double T_min = 800.0; + su2double epsilon = 80; /*--- Define preferential dissociation coefficient ---*/ //alpha = 0.3; @@ -863,14 +863,13 @@ vector& CSU2TCLib::ComputeNetProductionRates(){ void CSU2TCLib::ComputeKeqConstants(unsigned short val_Reaction) { - unsigned short ii, iIndex, tbl_offset, pwr; - su2double N, tmp1, tmp2; + unsigned short ii; /*--- Acquire database constants from CConfig ---*/ GetChemistryEquilConstants(val_Reaction); /*--- Calculate mixture number density ---*/ - N = 0.0; + su2double N = 0.0; for (iSpecies =0 ; iSpecies < nSpecies; iSpecies++) { N += rhos[iSpecies]/MolarMass[iSpecies]*AVOGAD_CONSTANT; } @@ -879,11 +878,11 @@ void CSU2TCLib::ComputeKeqConstants(unsigned short val_Reaction) { N = N*(1E-6); /*--- Determine table index based on mixture N ---*/ - tbl_offset = 14; - pwr = floor(log10(N)); + unsigned short tbl_offset = 14; + unsigned short pwr = floor(log10(N)); /*--- Bound the interpolation to table limit values ---*/ - iIndex = int(pwr) - tbl_offset; + unsigned short iIndex = int(pwr) - tbl_offset; if (iIndex <= 0) { for (ii = 0; ii < 5; ii++) A[ii] = RxnConstantTable(0,ii); @@ -895,8 +894,8 @@ void CSU2TCLib::ComputeKeqConstants(unsigned short val_Reaction) { } /*--- Calculate interpolation denominator terms avoiding pow() ---*/ - tmp1 = 1.0; - tmp2 = 1.0; + su2double tmp1 = 1.0; + su2double tmp2 = 1.0; for (ii = 0; ii < pwr; ii++) { tmp1 *= 10.0; tmp2 *= 10.0; @@ -918,20 +917,22 @@ su2double CSU2TCLib::ComputeEveSourceTerm(){ // Note: Landau-Teller formulation // Note: Millikan & White relaxation time (requires P in Atm.) // Note: Park limiting cross section - su2double conc, N, mu, A_sr, B_sr, num, denom, Cs, sig_s, tau_sr, tauP, tauMW, taus, omegaVT, omegaCV; + su2double A_sr, B_sr, num, denom, Cs, sig_s, tau_sr, tauP, tauMW, taus; vector MolarFrac, eve_eq, eve; + su2activematrix mu; MolarFrac.resize(nSpecies,0.0); eve_eq.resize(nSpecies,0.0); eve.resize(nSpecies,0.0); + mu.resize(nSpecies,nSpecies)=su2double(0.0); - omegaVT = 0.0; - omegaCV = 0.0; + su2double omegaVT = 0.0; + su2double omegaCV = 0.0; /*--- Calculate mole fractions ---*/ - N = 0.0; - conc = 0.0; + su2double N = 0.0; + su2double conc = 0.0; for (iSpecies = 0; iSpecies < nSpecies; iSpecies++) { conc += rhos[iSpecies] / MolarMass[iSpecies]; N += rhos[iSpecies] / MolarMass[iSpecies] * AVOGAD_CONSTANT; @@ -939,8 +940,8 @@ su2double CSU2TCLib::ComputeEveSourceTerm(){ for (iSpecies = 0; iSpecies < nSpecies; iSpecies++) MolarFrac[iSpecies] = (rhos[iSpecies] / MolarMass[iSpecies]) / conc; - eve_eq = ComputeSpeciesEve(T); - eve = ComputeSpeciesEve(Tve); + eve_eq = ComputeSpeciesEve(T, true); + eve = ComputeSpeciesEve(Tve, true); /*--- Loop over species to calculate source term --*/ for (iSpecies = 0; iSpecies < nSpecies; iSpecies++) { @@ -949,10 +950,11 @@ su2double CSU2TCLib::ComputeEveSourceTerm(){ num = 0.0; denom = 0.0; for (jSpecies = 0; jSpecies < nSpecies; jSpecies++) { - mu = MolarMass[iSpecies]*MolarMass[jSpecies] / (MolarMass[iSpecies] + MolarMass[jSpecies]); - A_sr = 1.16 * 1E-3 * sqrt(mu) * pow(CharVibTemp[iSpecies], 4.0/3.0); - B_sr = 0.015 * pow(mu, 0.25); - tau_sr = 101325.0/Pressure * exp(A_sr*(pow(T,-1.0/3.0) - B_sr) - 18.42); + mu(iSpecies,jSpecies) = MolarMass[iSpecies]*MolarMass[jSpecies] / (MolarMass[iSpecies] + MolarMass[jSpecies]); + A_sr = 1.16 * 1E-3 * sqrt(mu(iSpecies,jSpecies)) * pow(CharVibTemp[iSpecies], 4.0/3.0); + B_sr = 0.015 * pow(mu(iSpecies,jSpecies), 0.25); + tau_sr = 101325.0/Pressure * exp(A_sr*(pow(T,-1.0/3.0) - B_sr) - 18.42); + num += MolarFrac[jSpecies]; denom += MolarFrac[jSpecies] / tau_sr; } @@ -961,7 +963,7 @@ su2double CSU2TCLib::ComputeEveSourceTerm(){ /*--- Park limiting cross section ---*/ Cs = sqrt((8.0*Ru*T)/(PI_NUMBER*MolarMass[iSpecies])); - sig_s = 1E-20*(5E4*5E4)/(T*T); + sig_s = 3E-21*(2.5E9)/(T*T); tauP = 1/(sig_s*Cs*N); @@ -1035,13 +1037,13 @@ vector& CSU2TCLib::GetThermalConductivities(){ void CSU2TCLib::DiffusionCoeffWBE(){ - su2double conc, Mi, Mj, M, Omega_ij, denom; + su2double Mi, Mj, Omega_ij, denom; su2activematrix Dij; Dij.resize(nSpecies, nSpecies) = su2double(0.0); /*--- Calculate species mole fraction ---*/ - conc = 0.0; + su2double conc = 0.0; for (iSpecies = 0; iSpecies < nSpecies; iSpecies++) { MolarFracWBE[iSpecies] = rhos[iSpecies]/MolarMass[iSpecies]; conc += MolarFracWBE[iSpecies]; @@ -1050,7 +1052,7 @@ void CSU2TCLib::DiffusionCoeffWBE(){ MolarFracWBE[iSpecies] = MolarFracWBE[iSpecies]/conc; /*--- Calculate mixture molar mass (kg/mol) ---*/ // Note: Species molar masses stored as kg/kmol, need 1E-3 conversion - M = 0.0; + su2double M = 0.0; for (iSpecies = 0; iSpecies < nSpecies; iSpecies++) M += MolarMass[iSpecies]*MolarFracWBE[iSpecies]; M = M*1E-3; @@ -1090,10 +1092,10 @@ void CSU2TCLib::DiffusionCoeffWBE(){ void CSU2TCLib::ViscosityWBE(){ - su2double tmp1, tmp2, conc; + su2double tmp1, tmp2; /*--- Calculate species mole fraction ---*/ - conc = 0.0; + su2double conc = 0.0; for (iSpecies = 0; iSpecies < nSpecies; iSpecies++) { MolarFracWBE[iSpecies] = rhos[iSpecies]/MolarMass[iSpecies]; conc += MolarFracWBE[iSpecies]; @@ -1151,13 +1153,13 @@ void CSU2TCLib::ThermalConductivitiesWBE(){ void CSU2TCLib::DiffusionCoeffGY(){ - su2double Mi, Mj, pi, kb, gam_i, gam_j, gam_t, denom, d1_ij, D_ij, Omega_ij; + su2double Mi, Mj, gam_i, gam_j, denom, d1_ij, D_ij, Omega_ij; - pi = PI_NUMBER; - kb = BOLTZMANN_CONSTANT; + su2double pi = PI_NUMBER; + su2double kb = BOLTZMANN_CONSTANT; /*--- Calculate mixture gas constant ---*/ - gam_t = 0.0; + su2double gam_t = 0.0; for (iSpecies = 0; iSpecies < nSpecies; iSpecies++) { gam_t += rhos[iSpecies] / (Density*MolarMass[iSpecies]); } @@ -1238,10 +1240,10 @@ void CSU2TCLib::DiffusionCoeffGY(){ void CSU2TCLib::ViscosityGY(){ - su2double Mi, Mj, pi, Na, gam_i, gam_j, denom, Omega_ij, d2_ij; + su2double Mi, Mj, gam_i, gam_j, denom, Omega_ij, d2_ij; - pi = PI_NUMBER; - Na = AVOGAD_CONSTANT; + su2double pi = PI_NUMBER; + su2double Na = AVOGAD_CONSTANT; Mu = 0.0; /*--- Mixture viscosity via Gupta-Yos approximation ---*/ for (iSpecies = 0; iSpecies < nHeavy; iSpecies++) { @@ -1300,11 +1302,11 @@ void CSU2TCLib::ViscosityGY(){ void CSU2TCLib::ThermalConductivitiesGY(){ - su2double Cvve, Mi, Mj, mi, mj, pi, R, Na, kb, gam_i, gam_j, denom_t, denom_r, d1_ij, d2_ij, a_ij, Omega_ij, rhoCvve; + su2double Mi, Mj, mi, mj, gam_i, gam_j, denom_t, denom_r, d1_ij, d2_ij, a_ij, Omega_ij; - pi = PI_NUMBER; - Na = AVOGAD_CONSTANT; - kb = BOLTZMANN_CONSTANT; + su2double pi = PI_NUMBER; + su2double Na = AVOGAD_CONSTANT; + su2double kb = BOLTZMANN_CONSTANT; if (ionization) { SU2_MPI::Error("NEEDS REVISION w/ IONIZATION",CURRENT_FUNCTION); @@ -1312,18 +1314,18 @@ void CSU2TCLib::ThermalConductivitiesGY(){ /*--- Mixture vibrational-electronic specific heat ---*/ Cvves = ComputeSpeciesCvVibEle(); - rhoCvve = 0.0; + su2double rhoCvve = 0.0; for (iSpecies = 0; iSpecies < nSpecies; iSpecies++) rhoCvve += rhos[iSpecies]*Cvves[iSpecies]; - Cvve = rhoCvve/Density; + su2double Cvve = rhoCvve/Density; /*--- Calculate mixture gas constant ---*/ - R = 0.0; + su2double R = 0.0; for (iSpecies = 0; iSpecies < nSpecies; iSpecies++) { R += Ru * rhos[iSpecies]/Density; } /*--- Mixture thermal conductivity via Gupta-Yos approximation ---*/ - ThermalCond_tr = 0.0; + ThermalCond_tr = 0.0; ThermalCond_ve = 0.0; for (iSpecies = 0; iSpecies < nSpecies; iSpecies++) { /*--- Calculate molar concentration ---*/ @@ -1371,17 +1373,13 @@ void CSU2TCLib::ThermalConductivitiesGY(){ vector& CSU2TCLib::ComputeTemperatures(vector& val_rhos, su2double rhoE, su2double rhoEve, su2double rhoEvel){ vector val_eves; - su2double rhoCvtr, rhoE_f, rhoE_ref, rhoEve_t, Tve2, Tve_o, Btol, Tmin, Tmax; - bool Bconvg; - unsigned short iIter, maxBIter; - rhos = val_rhos; /*----------Translational temperature----------*/ - rhoE_f = 0.0; - rhoE_ref = 0.0; - rhoCvtr = 0.0; + su2double rhoE_f = 0.0; + su2double rhoE_ref = 0.0; + su2double rhoCvtr = 0.0; for (iSpecies = 0; iSpecies < nHeavy; iSpecies++) { rhoCvtr += rhos[iSpecies] * Cvtrs[iSpecies]; rhoE_ref += rhos[iSpecies] * Cvtrs[iSpecies] * Ref_Temperature[iSpecies]; @@ -1391,24 +1389,25 @@ vector& CSU2TCLib::ComputeTemperatures(vector& val_rhos, s T = (rhoE - rhoEve - rhoE_f + rhoE_ref - rhoEvel) / rhoCvtr; /*--- Set temperature clipping values ---*/ - Tmin = 50.0; Tmax = 8E4; - Tve_o = 50.0; Tve2 = 8E4; + su2double Tmin = 50.0; su2double Tmax = 8E4; + su2double Tve_o = 50.0; su2double Tve2 = 8E4; /* Determine if the temperature lies within the acceptable range */ if (T < Tmin) T = Tmin; else if (T > Tmax) T = Tmax; /*--- Set vibrational temperature algorithm parameters ---*/ - Btol = 1.0E-6; // Tolerance for the Bisection method - maxBIter = 50; // Maximum Bisection method iterations + su2double Btol = 1.0E-6; // Tolerance for the Bisection method + unsigned short maxBIter = 50; // Maximum Bisection method iterations //Initialize solution Tve = T; // Execute the root-finding method - Bconvg = false; + bool Bconvg = false; + su2double rhoEve_t; - for (iIter = 0; iIter < maxBIter; iIter++) { + for (unsigned short iIter = 0; iIter < maxBIter; iIter++) { Tve = (Tve_o+Tve2)/2.0; val_eves = ComputeSpeciesEve(Tve); rhoEve_t = 0.0; @@ -1422,7 +1421,10 @@ vector& CSU2TCLib::ComputeTemperatures(vector& val_rhos, s } } // If absolutely no convergence, then assign to the TR temperature - if (!Bconvg) Tve = T; + if (!Bconvg) { + Tve = T; + cout <<"Warning: temperatures did not converge, error= "<< fabs(rhoEve_t-rhoEve)< CSource_NEMO::ComputeAxisymmetric(const CConfig *confi for (iSpecies=0; iSpeciesnodes->GetCoord(iPoint, 1); - su2double yVelocity = nodes->GetVelocity(iPoint,1); - su2double xVelocity = nodes->GetVelocity(iPoint,0); - su2double Total_Viscosity = nodes->GetLaminarViscosity(iPoint) + nodes->GetEddyViscosity(iPoint); - - if (yCoord > EPS){ - su2double nu_v_on_y = Total_Viscosity*yVelocity/yCoord; - nodes->SetAuxVar(iPoint, 0, nu_v_on_y); - nodes->SetAuxVar(iPoint, 1, nu_v_on_y*yVelocity); - nodes->SetAuxVar(iPoint, 2, nu_v_on_y*xVelocity); - } - } - - /*--- Compute the auxiliary variable gradient with GG or WLS. ---*/ - if (config->GetKind_Gradient_Method() == GREEN_GAUSS) { - SetAuxVar_Gradient_GG(geometry, config); - } - if (config->GetKind_Gradient_Method() == WEIGHTED_LEAST_SQUARES) { - SetAuxVar_Gradient_LS(geometry, config); - } - + ComputeAxisymmetricAuxGradients(geometry, config); } /*--- loop over points ---*/ diff --git a/SU2_CFD/src/solvers/CNEMOEulerSolver.cpp b/SU2_CFD/src/solvers/CNEMOEulerSolver.cpp index e596e225ae08..4f88dee8d14e 100644 --- a/SU2_CFD/src/solvers/CNEMOEulerSolver.cpp +++ b/SU2_CFD/src/solvers/CNEMOEulerSolver.cpp @@ -815,45 +815,49 @@ void CNEMOEulerSolver::Source_Residual(CGeometry *geometry, CSolver **solver_con unsigned short iVar, jVar; unsigned long iPoint; - unsigned long eAxi_local, eChm_local, eVib_local; - unsigned long eAxi_global, eChm_global, eVib_global; /*--- Assign booleans ---*/ bool err = false; bool implicit = (config->GetKind_TimeIntScheme_Flow() == EULER_IMPLICIT); bool frozen = config->GetFrozen(); bool monoatomic = config->GetMonoatomic(); + bool axisymm = config->GetAxisymmetric(); bool viscous = config->GetViscous(); bool rans = (config->GetKind_Turb_Model() != NONE); CNumerics* numerics = numerics_container[SOURCE_FIRST_TERM]; /*--- Initialize the error counter ---*/ - eAxi_local = 0; - eChm_local = 0; - eVib_local = 0; + unsigned long eAxi_local = 0; + unsigned long eChm_local = 0; + unsigned long eVib_local = 0; /*--- Initialize the source residual to zero ---*/ for (iVar = 0; iVar < nVar; iVar++) Residual[iVar] = 0.0; + /*--- Preprocess viscous axisymm variables (if necessary) ---*/ + if (axisymm && viscous) { + ComputeAxisymmetricAuxGradients(geometry,config); + } + /*--- loop over interior points ---*/ + SU2_OMP_FOR_DYN(omp_chunk_size) for (iPoint = 0; iPoint < nPointDomain; iPoint++) { /*--- Set conserved & primitive variables ---*/ - numerics->SetConservative(nodes->GetSolution(iPoint), nodes->GetSolution(iPoint)); - numerics->SetPrimitive (nodes->GetPrimitive(iPoint), nodes->GetPrimitive(iPoint) ); + numerics->SetConservative(nodes->GetSolution(iPoint), nullptr); + numerics->SetPrimitive (nodes->GetPrimitive(iPoint), nullptr); /*--- Pass supplementary information to CNumerics ---*/ - numerics->SetdPdU(nodes->GetdPdU(iPoint), nodes->GetdPdU(iPoint)); - numerics->SetdTdU(nodes->GetdTdU(iPoint), nodes->GetdTdU(iPoint)); - numerics->SetdTvedU(nodes->GetdTvedU(iPoint), nodes->GetdTvedU(iPoint)); - numerics->SetEve(nodes->GetEve(iPoint), nodes->GetEve(iPoint)); - numerics->SetCvve(nodes->GetCvve(iPoint), nodes->GetCvve(iPoint)); + numerics->SetdPdU(nodes->GetdPdU(iPoint), nullptr); + numerics->SetdTdU(nodes->GetdTdU(iPoint), nullptr); + numerics->SetdTvedU(nodes->GetdTvedU(iPoint), nullptr); + numerics->SetEve(nodes->GetEve(iPoint), nullptr); + numerics->SetCvve(nodes->GetCvve(iPoint), nullptr); /*--- Set volume of the dual grid cell ---*/ numerics->SetVolume(geometry->nodes->GetVolume(iPoint)); - numerics->SetCoord(geometry->nodes->GetCoord(iPoint), - geometry->nodes->GetCoord(iPoint) ); + numerics->SetCoord(geometry->nodes->GetCoord(iPoint), nullptr); /*--- Compute finite rate chemistry ---*/ @@ -905,97 +909,66 @@ void CNEMOEulerSolver::Source_Residual(CGeometry *geometry, CSolver **solver_con eVib_local++; } - } /*--- Compute axisymmetric source terms (if needed) ---*/ - if (config->GetAxisymmetric()) { + if (axisymm) { + /*--- If necessary, set variables needed for viscous computation ---*/ if (viscous) { - for (iPoint = 0; iPoint < nPoint; iPoint++) { + /*--- Set gradient of primitive variables ---*/ + numerics->SetPrimVarGradient(nodes->GetGradient_Primitive(iPoint), nullptr); - su2double yCoord = geometry->nodes->GetCoord(iPoint, 1); - su2double yVelocity = nodes->GetVelocity(iPoint,1); - su2double xVelocity = nodes->GetVelocity(iPoint,0); - su2double Total_Viscosity = nodes->GetLaminarViscosity(iPoint) + nodes->GetEddyViscosity(iPoint); + /*--- Set gradient of auxillary variables ---*/ + numerics->SetAuxVarGrad(nodes->GetAuxVarGradient(iPoint), nullptr); - if (yCoord > EPS){ - su2double nu_v_on_y = Total_Viscosity*yVelocity/yCoord; - nodes->SetAuxVar(iPoint, 0, nu_v_on_y); - nodes->SetAuxVar(iPoint, 1, nu_v_on_y*yVelocity); - nodes->SetAuxVar(iPoint, 2, nu_v_on_y*xVelocity); - } - } - - /*--- Compute the auxiliary variable gradient with GG or WLS. ---*/ - if (config->GetKind_Gradient_Method() == GREEN_GAUSS) { - SetAuxVar_Gradient_GG(geometry, config); - } - if (config->GetKind_Gradient_Method() == WEIGHTED_LEAST_SQUARES) { - SetAuxVar_Gradient_LS(geometry, config); - } - } + /*--- Set diffusion coefficient ---*/ + numerics->SetDiffusionCoeff(nodes->GetDiffusionCoeff(iPoint), nullptr); - /*--- loop over points ---*/ - SU2_OMP_FOR_DYN(omp_chunk_size) - for (iPoint = 0; iPoint < nPointDomain; iPoint++) { - - /*--- If necessary, set variables needed for viscous computation ---*/ - if (viscous) { - - /*--- Set gradient of primitive variables ---*/ - numerics->SetPrimVarGradient(nodes->GetGradient_Primitive(iPoint), nodes->GetGradient_Primitive(iPoint)); - - /*--- Set gradient of auxillary variables ---*/ - numerics->SetAuxVarGrad(nodes->GetAuxVarGradient(iPoint), nullptr); - - /*--- Set diffusion coefficient ---*/ - numerics->SetDiffusionCoeff(nodes->GetDiffusionCoeff(iPoint), nodes->GetDiffusionCoeff(iPoint)); - - /*--- Laminar viscosity ---*/ - numerics->SetLaminarViscosity(nodes->GetLaminarViscosity(iPoint), nodes->GetLaminarViscosity(iPoint)); + /*--- Laminar viscosity ---*/ + numerics->SetLaminarViscosity(nodes->GetLaminarViscosity(iPoint), 0.0); - /*--- Eddy viscosity ---*/ - numerics->SetEddyViscosity(nodes->GetEddyViscosity(iPoint), nodes->GetEddyViscosity(iPoint)); + /*--- Eddy viscosity ---*/ + numerics->SetEddyViscosity(nodes->GetEddyViscosity(iPoint), 0.0); - /*--- Thermal conductivity ---*/ - numerics->SetThermalConductivity(nodes->GetThermalConductivity(iPoint), nodes->GetThermalConductivity(iPoint)); + /*--- Thermal conductivity ---*/ + numerics->SetThermalConductivity(nodes->GetThermalConductivity(iPoint), 0.0); - /*--- Vib-el. thermal conductivity ---*/ - numerics->SetThermalConductivity_ve(nodes->GetThermalConductivity_ve(iPoint), nodes->GetThermalConductivity_ve(iPoint)); + /*--- Vib-el. thermal conductivity ---*/ + numerics->SetThermalConductivity_ve(nodes->GetThermalConductivity_ve(iPoint), 0.0); - /*--- Set turbulence kinetic energy ---*/ - if (rans){ - CVariable* turbNodes = solver_container[TURB_SOL]->GetNodes(); - numerics->SetTurbKineticEnergy(turbNodes->GetSolution(iPoint,0), turbNodes->GetSolution(iPoint,0)); - } + /*--- Set turbulence kinetic energy ---*/ + if (rans){ + CVariable* turbNodes = solver_container[TURB_SOL]->GetNodes(); + numerics->SetTurbKineticEnergy(turbNodes->GetSolution(iPoint,0), 0.0); } + } - auto residual = numerics->ComputeAxisymmetric(config); + auto residual = numerics->ComputeAxisymmetric(config); - /*--- Check for errors before applying source to the linear system ---*/ - err = false; + /*--- Check for errors before applying source to the linear system ---*/ + err = false; + for (iVar = 0; iVar < nVar; iVar++) + if (residual[iVar] != residual[iVar]) err = true; + if (implicit) for (iVar = 0; iVar < nVar; iVar++) - if (residual[iVar] != residual[iVar]) err = true; - if (implicit) - for (iVar = 0; iVar < nVar; iVar++) - for (jVar = 0; jVar < nVar; jVar++) - if (Jacobian_i[iVar][jVar] != Jacobian_i[iVar][jVar]) err = true; + for (jVar = 0; jVar < nVar; jVar++) + if (Jacobian_i[iVar][jVar] != Jacobian_i[iVar][jVar]) err = true; - /*--- Apply the update to the linear system ---*/ - if (!err) { + /*--- Apply the update to the linear system ---*/ + if (!err) { LinSysRes.AddBlock(iPoint, residual); if (implicit) Jacobian.AddBlock(iPoint, iPoint, Jacobian_i); - }else - eAxi_local++; - } - END_SU2_OMP_FOR + }else + eAxi_local++; } + } + END_SU2_OMP_FOR /*--- Checking for NaN ---*/ - eAxi_global = eAxi_local; - eChm_global = eChm_local; - eVib_global = eVib_local; + unsigned long eAxi_global = eAxi_local; + unsigned long eChm_global = eChm_local; + unsigned long eVib_global = eVib_local; //THIS IS NO FUN if ((eAxi_global != 0) || @@ -1316,8 +1289,8 @@ void CNEMOEulerSolver::SetNondimensionalization(CConfig *config, unsigned short ModelTableOut <<"-- Models:"<< endl; ModelTable.AddColumn("Mixture", 25); - ModelTable.AddColumn("Transport Model", 25); ModelTable.AddColumn("Fluid Model", 25); + ModelTable.AddColumn("Transport Model", 25); ModelTable.SetAlign(PrintingToolbox::CTablePrinter::RIGHT); ModelTable.PrintHeader(); @@ -1670,23 +1643,23 @@ void CNEMOEulerSolver::BC_Far_Field(CGeometry *geometry, CSolver **solver_contai /*--- Species diffusion coefficients ---*/ visc_numerics->SetDiffusionCoeff(nodes->GetDiffusionCoeff(iPoint), - node_infty->GetDiffusionCoeff(0) ); + nodes->GetDiffusionCoeff(iPoint)); /*--- Laminar viscosity ---*/ visc_numerics->SetLaminarViscosity(nodes->GetLaminarViscosity(iPoint), - node_infty->GetLaminarViscosity(0) ); + nodes->GetLaminarViscosity(iPoint)); /*--- Eddy viscosity ---*/ visc_numerics->SetEddyViscosity(nodes->GetEddyViscosity(iPoint), - node_infty->GetEddyViscosity(0) ); + nodes->GetEddyViscosity(iPoint)); /*--- Thermal conductivity ---*/ visc_numerics->SetThermalConductivity(nodes->GetThermalConductivity(iPoint), - node_infty->GetThermalConductivity(0)); + nodes->GetThermalConductivity(iPoint)); /*--- Vib-el. thermal conductivity ---*/ visc_numerics->SetThermalConductivity_ve(nodes->GetThermalConductivity_ve(iPoint), - node_infty->GetThermalConductivity_ve(0) ); + nodes->GetThermalConductivity_ve(iPoint)); /*--- Compute and update residual ---*/ auto residual = visc_numerics->ComputeResidual(config); @@ -1703,7 +1676,7 @@ void CNEMOEulerSolver::BC_Far_Field(CGeometry *geometry, CSolver **solver_contai delete [] Normal; } -void CNEMOEulerSolver::BC_Inlet(CGeometry *geometry, CSolver **solution_container, +void CNEMOEulerSolver::BC_Inlet(CGeometry *geometry, CSolver **solver_container, CNumerics *conv_numerics, CNumerics *visc_numerics, CConfig *config, unsigned short val_marker) { @@ -1982,7 +1955,7 @@ void CNEMOEulerSolver::BC_Inlet(CGeometry *geometry, CSolver **solution_containe delete [] Spec_Density; } -void CNEMOEulerSolver::BC_Outlet(CGeometry *geometry, CSolver **solution_container, +void CNEMOEulerSolver::BC_Outlet(CGeometry *geometry, CSolver **solver_container, CNumerics *conv_numerics, CNumerics *visc_numerics, CConfig *config, unsigned short val_marker) { unsigned short iVar, iDim, iSpecies; @@ -2232,7 +2205,7 @@ void CNEMOEulerSolver::BC_Outlet(CGeometry *geometry, CSolver **solution_contain } -void CNEMOEulerSolver::BC_Supersonic_Inlet(CGeometry *geometry, CSolver **solution_container, +void CNEMOEulerSolver::BC_Supersonic_Inlet(CGeometry *geometry, CSolver **solver_container, CNumerics *conv_numerics, CNumerics *visc_numerics, CConfig *config, unsigned short val_marker) { SU2_MPI::Error("BC_SUPERSONIC_INLET: Not operational in NEMO.", CURRENT_FUNCTION); @@ -2443,7 +2416,7 @@ SU2_MPI::Error("BC_SUPERSONIC_INLET: Not operational in NEMO.", CURRENT_FUNCTION } -void CNEMOEulerSolver::BC_Supersonic_Outlet(CGeometry *geometry, CSolver **solution_container, +void CNEMOEulerSolver::BC_Supersonic_Outlet(CGeometry *geometry, CSolver **solver_container, CNumerics *conv_numerics, CNumerics *visc_numerics, CConfig *config, unsigned short val_marker) { unsigned short iDim; unsigned long iVertex, iPoint; diff --git a/SU2_CFD/src/solvers/CNEMONSSolver.cpp b/SU2_CFD/src/solvers/CNEMONSSolver.cpp index b4c9ff25aef0..defd777c573c 100644 --- a/SU2_CFD/src/solvers/CNEMONSSolver.cpp +++ b/SU2_CFD/src/solvers/CNEMONSSolver.cpp @@ -278,7 +278,7 @@ void CNEMONSSolver::Viscous_Residual(CGeometry *geometry, } void CNEMONSSolver::BC_HeatFluxNonCatalytic_Wall(CGeometry *geometry, - CSolver **solution_container, + CSolver **solver_container, CNumerics *conv_numerics, CNumerics *sour_numerics, CConfig *config, @@ -384,7 +384,7 @@ void CNEMONSSolver::BC_HeatFluxNonCatalytic_Wall(CGeometry *geometry, } void CNEMONSSolver::BC_HeatFlux_Wall(CGeometry *geometry, - CSolver **solution_container, + CSolver **solver_container, CNumerics *conv_numerics, CNumerics *sour_numerics, CConfig *config, @@ -406,7 +406,7 @@ void CNEMONSSolver::BC_HeatFlux_Wall(CGeometry *geometry, if (Marker_Tag == Catalytic_Tag){ catalytic = true; - BC_HeatFluxCatalytic_Wall(geometry, solution_container, conv_numerics, + BC_HeatFluxCatalytic_Wall(geometry, solver_container, conv_numerics, sour_numerics, config, val_marker); break; @@ -417,14 +417,14 @@ void CNEMONSSolver::BC_HeatFlux_Wall(CGeometry *geometry, } } - if(!catalytic) BC_HeatFluxNonCatalytic_Wall(geometry, solution_container, conv_numerics, + if(!catalytic) BC_HeatFluxNonCatalytic_Wall(geometry, solver_container, conv_numerics, sour_numerics, config, val_marker); } void CNEMONSSolver::BC_HeatFluxCatalytic_Wall(CGeometry *geometry, - CSolver **solution_container, + CSolver **solver_container, CNumerics *conv_numerics, CNumerics *sour_numerics, CConfig *config, @@ -593,7 +593,7 @@ void CNEMONSSolver::BC_HeatFluxCatalytic_Wall(CGeometry *geometry, } void CNEMONSSolver::BC_Isothermal_Wall(CGeometry *geometry, - CSolver **solution_container, + CSolver **solver_container, CNumerics *conv_numerics, CNumerics *sour_numerics, CConfig *config, @@ -615,7 +615,7 @@ void CNEMONSSolver::BC_Isothermal_Wall(CGeometry *geometry, if (Marker_Tag == Catalytic_Tag){ catalytic = true; - BC_IsothermalCatalytic_Wall(geometry, solution_container, conv_numerics, + BC_IsothermalCatalytic_Wall(geometry, solver_container, conv_numerics, sour_numerics, config, val_marker); break; } else { @@ -623,12 +623,12 @@ void CNEMONSSolver::BC_Isothermal_Wall(CGeometry *geometry, } } - if(!catalytic) BC_IsothermalNonCatalytic_Wall(geometry, solution_container, conv_numerics, + if(!catalytic) BC_IsothermalNonCatalytic_Wall(geometry, solver_container, conv_numerics, sour_numerics, config, val_marker); } void CNEMONSSolver::BC_IsothermalNonCatalytic_Wall(CGeometry *geometry, - CSolver **solution_container, + CSolver **solver_container, CNumerics *conv_numerics, CNumerics *sour_numerics, CConfig *config, @@ -756,7 +756,7 @@ void CNEMONSSolver::BC_IsothermalNonCatalytic_Wall(CGeometry *geometry, } void CNEMONSSolver::BC_IsothermalCatalytic_Wall(CGeometry *geometry, - CSolver **solution_container, + CSolver **solver_container, CNumerics *conv_numerics, CNumerics *sour_numerics, CConfig *config, @@ -765,7 +765,7 @@ void CNEMONSSolver::BC_IsothermalCatalytic_Wall(CGeometry *geometry, SU2_MPI::Error("BC_ISOTHERMAL with catalytic wall: Not operational in NEMO.", CURRENT_FUNCTION); /*--- Call standard isothermal BC to apply no-slip and energy b.c.'s ---*/ - BC_IsothermalNonCatalytic_Wall(geometry, solution_container, conv_numerics, + BC_IsothermalNonCatalytic_Wall(geometry, solver_container, conv_numerics, sour_numerics, config, val_marker); ///////////// FINITE DIFFERENCE METHOD /////////////// @@ -943,7 +943,7 @@ void CNEMONSSolver::BC_IsothermalCatalytic_Wall(CGeometry *geometry, } void CNEMONSSolver::BC_Smoluchowski_Maxwell(CGeometry *geometry, - CSolver **solution_container, + CSolver **solver_container, CNumerics *conv_numerics, CNumerics *visc_numerics, CConfig *config, diff --git a/SU2_CFD/src/variables/CNEMOEulerVariable.cpp b/SU2_CFD/src/variables/CNEMOEulerVariable.cpp index 6e96207631b9..f9b73aaf6fdd 100644 --- a/SU2_CFD/src/variables/CNEMOEulerVariable.cpp +++ b/SU2_CFD/src/variables/CNEMOEulerVariable.cpp @@ -192,22 +192,25 @@ void CNEMOEulerVariable::SetVelocity2(unsigned long iPoint) { bool CNEMOEulerVariable::SetPrimVar(unsigned long iPoint, CFluidModel *FluidModel) { - bool nonPhys; unsigned short iVar; fluidmodel = static_cast(FluidModel); /*--- Convert conserved to primitive variables ---*/ - nonPhys = Cons2PrimVar(Solution[iPoint], Primitive[iPoint], - dPdU[iPoint], dTdU[iPoint], dTvedU[iPoint], eves[iPoint], Cvves[iPoint]); + bool nonPhys = Cons2PrimVar(Solution[iPoint], Primitive[iPoint], + dPdU[iPoint], dTdU[iPoint], dTvedU[iPoint], eves[iPoint], Cvves[iPoint]); /*--- Reset solution to previous one, if nonphys ---*/ if (nonPhys) { for (iVar = 0; iVar < nVar; iVar++) Solution(iPoint,iVar) = Solution_Old(iPoint,iVar); + + /*--- Recompute Primitive from previous solution ---*/ + Cons2PrimVar(Solution[iPoint], Primitive[iPoint], + dPdU[iPoint], dTdU[iPoint], dTvedU[iPoint], eves[iPoint], Cvves[iPoint]); } - /*--- Set additional point quantaties ---*/ + /*--- Set additional point quantities ---*/ Gamma(iPoint) = fluidmodel->ComputeGamma(); SetVelocity2(iPoint); @@ -250,7 +253,7 @@ bool CNEMOEulerVariable::Cons2PrimVar(su2double *U, su2double *V, U[iSpecies] = 1E-20; V[RHOS_INDEX+iSpecies] = 1E-20; rhos[iSpecies] = 1E-20; - //nonPhys = true; + //nonPhys = true; } else { V[RHOS_INDEX+iSpecies] = U[iSpecies]; rhos[iSpecies] = U[iSpecies]; @@ -276,12 +279,10 @@ bool CNEMOEulerVariable::Cons2PrimVar(su2double *U, su2double *V, V[TVE_INDEX] = T[1]; // Determine if the temperature lies within the acceptable range - //TODO: fIX THIS - if (V[T_INDEX] == Tmin) { - nonPhys = true; - } else if (V[T_INDEX] == Tmax){ - nonPhys = true; - } + if (V[T_INDEX] <= Tmin) { nonPhys = true;} + else if (V[T_INDEX] >= Tmax) { nonPhys = true;} + else if (V[T_INDEX] != V[T_INDEX] || V[TVE_INDEX] != V[TVE_INDEX]){ + nonPhys = true;} /*--- Vibrational-Electronic Temperature ---*/ vector eves_min = fluidmodel->ComputeSpeciesEve(Tvemin); diff --git a/SU2_CFD/src/variables/CNEMONSVariable.cpp b/SU2_CFD/src/variables/CNEMONSVariable.cpp index 37f1864e0237..7b94ac8a5785 100644 --- a/SU2_CFD/src/variables/CNEMONSVariable.cpp +++ b/SU2_CFD/src/variables/CNEMONSVariable.cpp @@ -104,21 +104,23 @@ bool CNEMONSVariable::SetVorticity(void) { bool CNEMONSVariable::SetPrimVar(unsigned long iPoint, CFluidModel *FluidModel) { - bool nonPhys; unsigned short iVar, iSpecies; fluidmodel = static_cast(FluidModel); /*--- Convert conserved to primitive variables ---*/ - nonPhys = Cons2PrimVar(Solution[iPoint], Primitive[iPoint], dPdU[iPoint], dTdU[iPoint], dTvedU[iPoint], eves[iPoint], Cvves[iPoint]); + bool nonPhys = Cons2PrimVar(Solution[iPoint], Primitive[iPoint], dPdU[iPoint], dTdU[iPoint], dTvedU[iPoint], eves[iPoint], Cvves[iPoint]); /*--- Reset solution to previous one, if nonphys ---*/ if (nonPhys) { for (iVar = 0; iVar < nVar; iVar++) Solution(iPoint,iVar) = Solution_Old(iPoint,iVar); + + /*--- Recompute Primitive from previous solution ---*/ + Cons2PrimVar(Solution[iPoint], Primitive[iPoint], dPdU[iPoint], dTdU[iPoint], dTvedU[iPoint], eves[iPoint], Cvves[iPoint]); } - /*--- Set additional point quantaties ---*/ + /*--- Set additional point quantities ---*/ Gamma(iPoint) = fluidmodel->ComputeGamma(); SetVelocity2(iPoint); diff --git a/TestCases/hybrid_regression.py b/TestCases/hybrid_regression.py index af6011b0eaf1..583420413e78 100644 --- a/TestCases/hybrid_regression.py +++ b/TestCases/hybrid_regression.py @@ -227,7 +227,7 @@ def main(): axi_rans_air_nozzle.cfg_dir = "axisymmetric_rans/air_nozzle" axi_rans_air_nozzle.cfg_file = "air_nozzle.cfg" axi_rans_air_nozzle.test_iter = 10 - axi_rans_air_nozzle.test_vals = [-12.094937, -6.622043, -8.814412, -2.393288] + axi_rans_air_nozzle.test_vals = [-6.348077, -0.827162, -2.241982, 2.313210] test_list.append(axi_rans_air_nozzle) ################################# diff --git a/TestCases/nonequilibrium/viscwedge/viscwedge.cfg b/TestCases/nonequilibrium/axi_visccone/axi_visccone.cfg similarity index 98% rename from TestCases/nonequilibrium/viscwedge/viscwedge.cfg rename to TestCases/nonequilibrium/axi_visccone/axi_visccone.cfg index 57898e11b7c8..f07f08384ad9 100644 --- a/TestCases/nonequilibrium/viscwedge/viscwedge.cfg +++ b/TestCases/nonequilibrium/axi_visccone/axi_visccone.cfg @@ -25,6 +25,9 @@ MATH_PROBLEM= DIRECT % % Restart solution (NO, YES) RESTART_SOL= NO +% +% Axisymmetric simulation, only compressible flows (NO, YES) +AXISYMMETRIC= YES % ----------- COMPRESSIBLE AND INCOMPRESSIBLE FREE-STREAM DEFINITION ----------% % diff --git a/TestCases/parallel_regression.py b/TestCases/parallel_regression.py index 77faa7d3ecf9..db98e43c0d25 100644 --- a/TestCases/parallel_regression.py +++ b/TestCases/parallel_regression.py @@ -47,7 +47,7 @@ def main(): thermalbath.cfg_dir = "nonequilibrium/thermalbath/finitechemistry" thermalbath.cfg_file = "thermalbath.cfg" thermalbath.test_iter = 10 - thermalbath.test_vals = [2.473627, 2.473627, -11.989166, -11.879331, -32.000000, 10.804939] + thermalbath.test_vals = [0.945997, 0.945997, -12.039262, -12.171767, -32.000000, 10.013239] thermalbath.su2_exec = "mpirun -n 2 SU2_CFD" thermalbath.timeout = 1600 thermalbath.new_output = True @@ -59,7 +59,7 @@ def main(): thermalbath_frozen.cfg_dir = "nonequilibrium/thermalbath/frozen" thermalbath_frozen.cfg_file = "thermalbath_frozen.cfg" thermalbath_frozen.test_iter = 10 - thermalbath_frozen.test_vals = [ -32.000000, -32.000000, -11.92359, -11.962329, -32.000000, 10.813864] + thermalbath_frozen.test_vals = [-32.000000, -32.000000, -12.039251, -12.171781, -32.000000, 10.013545] thermalbath_frozen.su2_exec = "mpirun -n 2 SU2_CFD" thermalbath_frozen.timeout = 1600 thermalbath_frozen.new_output = True @@ -78,17 +78,17 @@ def main(): invwedge.tol = 0.00001 test_list.append(invwedge) - # Viscous single wedge - viscwedge = TestCase('viscwedge') - viscwedge.cfg_dir = "nonequilibrium/viscwedge" - viscwedge.cfg_file = "viscwedge.cfg" - viscwedge.test_iter = 10 - viscwedge.test_vals = [-5.170894,-5.695657,-20.831064,-20.718963,-23.419767,-1.559591,-2.068445,2.204714,-2.590194] - viscwedge.su2_exec = "mpirun -n 2 SU2_CFD" - viscwedge.timeout = 1600 - viscwedge.new_output = True - viscwedge.tol = 0.00001 - test_list.append(viscwedge) + # Viscous single cone - axisymmetric + visc_cone = TestCase('visc_cone') + visc_cone.cfg_dir = "nonequilibrium/axi_visccone" + visc_cone.cfg_file = "axi_visccone.cfg" + visc_cone.test_iter = 10 + visc_cone.test_vals = [-5.173082, -5.697844, -20.831292, -20.719160, -23.419769, -1.564064, -2.069008, 2.203919, -2.590729] + visc_cone.su2_exec = "mpirun -n 2 SU2_CFD" + visc_cone.timeout = 1600 + visc_cone.new_output = True + visc_cone.tol = 0.00001 + test_list.append(visc_cone) # Viscous single wedge with Mutation++ #viscwedge_mpp = TestCase('viscwedge_mpp') @@ -378,7 +378,7 @@ def main(): axi_rans_air_nozzle.cfg_dir = "axisymmetric_rans/air_nozzle" axi_rans_air_nozzle.cfg_file = "air_nozzle.cfg" axi_rans_air_nozzle.test_iter = 10 - axi_rans_air_nozzle.test_vals = [ -12.096569, -6.625843, -8.807541, -2.393279] + axi_rans_air_nozzle.test_vals = [-6.278454, -0.744813, -2.243285, 2.312481] axi_rans_air_nozzle.su2_exec = "mpirun -n 2 SU2_CFD" axi_rans_air_nozzle.timeout = 1600 axi_rans_air_nozzle.tol = 0.0001 diff --git a/TestCases/serial_regression.py b/TestCases/serial_regression.py index 7d955ea5e9d3..807f81f9f793 100644 --- a/TestCases/serial_regression.py +++ b/TestCases/serial_regression.py @@ -45,7 +45,7 @@ def main(): thermalbath.cfg_dir = "nonequilibrium/thermalbath/finitechemistry" thermalbath.cfg_file = "thermalbath.cfg" thermalbath.test_iter = 10 - thermalbath.test_vals = [2.473627, 2.473627, -12.033039, -11.945257, -32.000000, 10.804939] #last 4 columns + thermalbath.test_vals = [0.945997, 0.945997, -12.018025, -12.217291, -32.000000, 10.013239] thermalbath.su2_exec = "SU2_CFD" thermalbath.timeout = 1600 thermalbath.new_output = True @@ -57,7 +57,7 @@ def main(): thermalbath_frozen.cfg_dir = "nonequilibrium/thermalbath/frozen" thermalbath_frozen.cfg_file = "thermalbath_frozen.cfg" thermalbath_frozen.test_iter = 10 - thermalbath_frozen.test_vals = [-32.000000, -32.000000, -11.953727, -12.066776, -32.000000, 10.813864] #last 4 columns + thermalbath_frozen.test_vals = [-32.000000, -32.000000, -12.018022, -12.217291, -32.000000, 10.013545] thermalbath_frozen.su2_exec = "SU2_CFD" thermalbath_frozen.timeout = 1600 thermalbath_frozen.new_output = True @@ -76,17 +76,17 @@ def main(): invwedge.tol = 0.00001 test_list.append(invwedge) - # Viscous single wedge - viscwedge = TestCase('viscwedge') - viscwedge.cfg_dir = "nonequilibrium/viscwedge" - viscwedge.cfg_file = "viscwedge.cfg" - viscwedge.test_iter = 10 - viscwedge.test_vals = [-5.170894,-5.695657,-20.831064,-20.718963,-23.419767,-1.559591,-2.068445,2.204714,-2.590194] - viscwedge.su2_exec = "SU2_CFD" - viscwedge.timeout = 1600 - viscwedge.new_output = True - viscwedge.tol = 0.00001 - test_list.append(viscwedge) + # Viscous single cone - axisymmetric + visc_cone = TestCase('visc_cone') + visc_cone.cfg_dir = "nonequilibrium/axi_visccone" + visc_cone.cfg_file = "axi_visccone.cfg" + visc_cone.test_iter = 10 + visc_cone.test_vals = [-5.173082, -5.697844, -20.831292, -20.719160, -23.419769, -1.564064, -2.069008, 2.203919, -2.590729] + visc_cone.su2_exec = "SU2_CFD" + visc_cone.timeout = 1600 + visc_cone.new_output = True + visc_cone.tol = 0.00001 + test_list.append(visc_cone) ######################### ## Compressible Euler ### @@ -395,7 +395,7 @@ def main(): axi_rans_air_nozzle.cfg_dir = "axisymmetric_rans/air_nozzle" axi_rans_air_nozzle.cfg_file = "air_nozzle.cfg" axi_rans_air_nozzle.test_iter = 10 - axi_rans_air_nozzle.test_vals = [ -12.093130, -6.619801, -8.806060, -2.393278] + axi_rans_air_nozzle.test_vals = [ -6.279299, -0.747627, -2.237638, 2.321596] axi_rans_air_nozzle.su2_exec = "SU2_CFD" axi_rans_air_nozzle.timeout = 1600 axi_rans_air_nozzle.tol = 0.0001