diff --git a/.github/workflows/regression.yml b/.github/workflows/regression.yml index af6f6c60875..2502204b64c 100644 --- a/.github/workflows/regression.yml +++ b/.github/workflows/regression.yml @@ -88,8 +88,6 @@ jobs: # flags: '--buildtype=debug -Denable-directdiff=true -Denable-normal=false -Dwith-omp=true -Denable-mixedprec=true -Denable-pywrapper=true -Denable-tecio=false --warnlevel=3 --werror' runs-on: ${{ inputs.runner || 'ubuntu-latest' }} steps: - - name: Reduce ASLR entropy for tsan - run: sudo sysctl -w vm.mmap_rnd_bits=28 - name: Cache Object Files uses: actions/cache@v4 with: @@ -130,8 +128,6 @@ jobs: flags: '--buildtype=debugoptimized -Denable-autodiff=true -Denable-normal=false -Dwith-mpi=disabled --warnlevel=3 --werror' runs-on: ${{ inputs.runner || 'ubuntu-latest' }} steps: - - name: Reduce ASLR entropy for asan - run: sudo sysctl -w vm.mmap_rnd_bits=28 - name: Cache Object Files uses: actions/cache@v4 with: @@ -229,8 +225,6 @@ jobs: matrix: testscript: ['hybrid_regression.py', 'hybrid_regression_AD.py'] steps: - - name: Reduce ASLR entropy for tsan - run: sudo sysctl -w vm.mmap_rnd_bits=28 - name: Pre Cleanup uses: docker://ghcr.io/su2code/su2/test-su2-tsan:240320-1536 with: @@ -276,8 +270,6 @@ jobs: matrix: testscript: ['serial_regression.py', 'serial_regression_AD.py'] steps: - - name: Reduce ASLR entropy for asan - run: sudo sysctl -w vm.mmap_rnd_bits=28 - name: Pre Cleanup uses: docker://ghcr.io/su2code/su2/test-su2-asan:240320-1536 with: diff --git a/Common/include/CConfig.hpp b/Common/include/CConfig.hpp index d4f94008ddc..df9a0dba3b7 100644 --- a/Common/include/CConfig.hpp +++ b/Common/include/CConfig.hpp @@ -4335,8 +4335,7 @@ class CConfig { array GetNewtonKrylovDblParam(void) const { return NK_DblParam; } /*! - * \brief Get the relaxation coefficient of the linear solver for the implicit formulation. - * \return relaxation coefficient of the linear solver for the implicit formulation. + * \brief Returns the Roe kappa (multipler of the dissipation term). */ su2double GetRoe_Kappa(void) const { return Roe_Kappa; } diff --git a/Common/src/CConfig.cpp b/Common/src/CConfig.cpp index bf2122f01a7..ec4a555d6af 100644 --- a/Common/src/CConfig.cpp +++ b/Common/src/CConfig.cpp @@ -3487,11 +3487,6 @@ void CConfig::SetPostprocessing(SU2_COMPONENT val_software, unsigned short val_i SU2_MPI::Error("COMPRESSIBILITY-WILCOX only supported for SOLVER= RANS", CURRENT_FUNCTION); } - /*--- Check if turbulence model can be used for AXISYMMETRIC case---*/ - if (Axisymmetric && Kind_Turb_Model != TURB_MODEL::NONE && Kind_Turb_Model != TURB_MODEL::SST){ - SU2_MPI::Error("Axisymmetry is currently only supported for KIND_TURB_MODEL chosen as SST", CURRENT_FUNCTION); - } - /*--- Postprocess LM_OPTIONS into structure. ---*/ if (Kind_Trans_Model == TURB_TRANS_MODEL::LM) { lmParsedOptions = ParseLMOptions(LM_Options, nLM_Options, rank, Kind_Turb_Model); diff --git a/SU2_CFD/include/numerics/flow/convection/hllc.hpp b/SU2_CFD/include/numerics/flow/convection/hllc.hpp index 5cac798825e..62d46a7d211 100644 --- a/SU2_CFD/include/numerics/flow/convection/hllc.hpp +++ b/SU2_CFD/include/numerics/flow/convection/hllc.hpp @@ -49,7 +49,7 @@ class CUpwHLLC_Flow final : public CNumerics { su2double sq_velRoe, RoeDensity, RoeEnthalpy, RoeSoundSpeed, RoeProjVelocity, ProjInterfaceVel; - su2double sL, sR, sM, pStar, EStar, rhoSL, rhoSR, Rrho, kappa; + su2double sL, sR, sM, pStar, EStar, rhoSL, rhoSR, Rrho; su2double Omega, RHO, OmegaSM; su2double *dSm_dU, *dPI_dU, *drhoStar_dU, *dpStar_dU, *dEStar_dU; @@ -102,7 +102,7 @@ class CUpwGeneralHLLC_Flow final : public CNumerics { su2double sq_velRoe, RoeDensity, RoeEnthalpy, RoeSoundSpeed, RoeProjVelocity, ProjInterfaceVel; su2double Kappa_i, Kappa_j, Chi_i, Chi_j, RoeKappa, RoeChi, RoeKappaStaticEnthalpy; - su2double sL, sR, sM, pStar, EStar, rhoSL, rhoSR, Rrho, kappa; + su2double sL, sR, sM, pStar, EStar, rhoSL, rhoSR, Rrho; su2double Omega, RHO, OmegaSM; su2double *dSm_dU, *dPI_dU, *drhoStar_dU, *dpStar_dU, *dEStar_dU; diff --git a/SU2_CFD/include/numerics/turbulent/turb_sources.hpp b/SU2_CFD/include/numerics/turbulent/turb_sources.hpp index a190d3de942..ce75b38f115 100644 --- a/SU2_CFD/include/numerics/turbulent/turb_sources.hpp +++ b/SU2_CFD/include/numerics/turbulent/turb_sources.hpp @@ -77,9 +77,41 @@ class CSourceBase_TurbSA : public CNumerics { const FlowIndices idx; /*!< \brief Object to manage the access to the flow primitives. */ const SA_ParsedOptions options; /*!< \brief Struct with SA options. */ + const bool axisymmetric = false; bool transition_LM; + /*! + * \brief Add contribution from diffusion due to axisymmetric formulation to 2D residual + */ + inline void ResidualAxisymmetricDiffusion(su2double sigma) { + if (Coord_i[1] < EPS) return; + + const su2double yinv = 1.0 / Coord_i[1]; + const su2double& nue = ScalarVar_i[0]; + + const auto& density = V_i[idx.Density()]; + const auto& laminar_viscosity = V_i[idx.LaminarViscosity()]; + + const su2double nu = laminar_viscosity/density; + + su2double nu_e; + + if (options.version == SA_OPTIONS::NEG && nue < 0.0) { + const su2double cn1 = 16.0; + const su2double Xi = nue / nu; + const su2double fn = (cn1 + Xi*Xi*Xi) / (cn1 - Xi*Xi*Xi); + nu_e = nu + fn * nue; + } else { + nu_e = nu + nue; + } + + /* Diffusion source term */ + const su2double dv_axi = (1.0/sigma)*nu_e*ScalarVar_Grad_i[0][1]; + + Residual += yinv * dv_axi * Volume; + } + public: /*! * \brief Constructor of the class. @@ -90,6 +122,7 @@ class CSourceBase_TurbSA : public CNumerics { : CNumerics(nDim, 1, config), idx(nDim, config->GetnSpecies()), options(config->GetSAParsedOptions()), + axisymmetric(config->GetAxisymmetric()), transition_LM(config->GetKind_Trans_Model() == TURB_TRANS_MODEL::LM) { /*--- Setup the Jacobian pointer, we need to return su2double** but we know * the Jacobian is 1x1 so we use this trick to avoid heap allocation. ---*/ @@ -217,6 +250,9 @@ class CSourceBase_TurbSA : public CNumerics { SourceTerms::get(ScalarVar_i[0], var, Production, Destruction, CrossProduction, Jacobian_i[0]); Residual = (Production - Destruction + CrossProduction) * Volume; + + if (axisymmetric) ResidualAxisymmetricDiffusion(var.sigma); + Jacobian_i[0] *= Volume; } @@ -520,6 +556,19 @@ class CCompressibilityCorrection final : public ParentClass { const su2double d_CompCorrection = 2.0 * c5 * ScalarVar_i[0] / pow(sound_speed, 2) * aux_cc * Volume; const su2double CompCorrection = 0.5 * ScalarVar_i[0] * d_CompCorrection; + /*--- Axisymmetric contribution ---*/ + if (this->axisymmetric && this->Coord_i[1] > EPS) { + const su2double yinv = 1.0 / this->Coord_i[1]; + const su2double nue = ScalarVar_i[0]; + const su2double v = V_i[idx.Velocity() + 1]; + + const su2double d_axiCorrection = 2.0 * c5 * nue * pow(v * yinv / sound_speed, 2) * Volume; + const su2double axiCorrection = 0.5 * nue * d_axiCorrection; + + this->Residual -= axiCorrection; + this->Jacobian_i[0] -= d_axiCorrection; + } + this->Residual -= CompCorrection; this->Jacobian_i[0] -= d_CompCorrection; diff --git a/SU2_CFD/include/solvers/CFVMFlowSolverBase.hpp b/SU2_CFD/include/solvers/CFVMFlowSolverBase.hpp index 51901b8201e..30687edbcdc 100644 --- a/SU2_CFD/include/solvers/CFVMFlowSolverBase.hpp +++ b/SU2_CFD/include/solvers/CFVMFlowSolverBase.hpp @@ -1099,22 +1099,21 @@ class CFVMFlowSolverBase : public CSolver { /*! * \brief Store of a set of provided inlet profile values at a vertex. * \param[in] val_inlet - vector containing the inlet values for the current vertex. - * \param[in] iMarker - Surface marker where the coefficient is computed. + * \param[in] iMarker - Index of the surface marker. * \param[in] iVertex - Vertex of the marker iMarker where the inlet is being set. */ void SetInletAtVertex(const su2double* val_inlet, unsigned short iMarker, unsigned long iVertex) final; /*! - * \brief Get the set of value imposed at an inlet. - * \param[in] val_inlet - vector returning the inlet values for the current vertex. - * \param[in] val_inlet_point - Node index where the inlet is being set. - * \param[in] val_kind_marker - Enumerated type for the particular inlet type. + * \brief Get the set of values imposed at an inlet. + * \param[in] iMarker - Index of the surface marker. + * \param[in] iVertex - Vertex of the marker iMarker where the inlet is being set. * \param[in] geometry - Geometrical definition of the problem. - * \param config - Definition of the particular problem. + * \param[in,out] val_inlet - vector returning the inlet values for the current vertex. * \return Value of the face area at the vertex. */ - su2double GetInletAtVertex(su2double* val_inlet, unsigned long val_inlet_point, unsigned short val_kind_marker, - string val_marker, const CGeometry* geometry, const CConfig* config) const final; + su2double GetInletAtVertex(unsigned short iMarker, unsigned long iVertex, + const CGeometry* geometry, su2double* val_inlet) const final; /*! * \author T. Kattmann diff --git a/SU2_CFD/include/solvers/CFVMFlowSolverBase.inl b/SU2_CFD/include/solvers/CFVMFlowSolverBase.inl index 8fcb6d05bec..f60db438e18 100644 --- a/SU2_CFD/include/solvers/CFVMFlowSolverBase.inl +++ b/SU2_CFD/include/solvers/CFVMFlowSolverBase.inl @@ -722,73 +722,48 @@ void CFVMFlowSolverBase::SetInletAtVertex(const su2double* val_inlet, unsi } template -su2double CFVMFlowSolverBase::GetInletAtVertex(su2double* val_inlet, unsigned long val_inlet_point, - unsigned short val_kind_marker, string val_marker, - const CGeometry* geometry, const CConfig* config) const { - /*--- Local variables ---*/ - - unsigned short iMarker, iDim; - unsigned long iPoint, iVertex; - su2double Area = 0.0; - su2double Normal[3] = {0.0, 0.0, 0.0}; - - /*--- Alias positions within inlet file for readability ---*/ - - unsigned short T_position = nDim; - unsigned short P_position = nDim + 1; - unsigned short FlowDir_position = nDim + 2; - - if (val_kind_marker == INLET_FLOW) { - for (iMarker = 0; iMarker < config->GetnMarker_All(); iMarker++) { - if ((config->GetMarker_All_KindBC(iMarker) == INLET_FLOW) && - (config->GetMarker_All_TagBound(iMarker) == val_marker)) { - for (iVertex = 0; iVertex < nVertex[iMarker]; iVertex++) { - iPoint = geometry->vertex[iMarker][iVertex]->GetNode(); - - if (iPoint == val_inlet_point) { - /*-- Compute boundary face area for this vertex. ---*/ - - geometry->vertex[iMarker][iVertex]->GetNormal(Normal); - Area = GeometryToolbox::Norm(nDim, Normal); - - /*--- Access and store the inlet variables for this vertex. ---*/ - - val_inlet[T_position] = Inlet_Ttotal[iMarker][iVertex]; - val_inlet[P_position] = Inlet_Ptotal[iMarker][iVertex]; - for (iDim = 0; iDim < nDim; iDim++) { - val_inlet[FlowDir_position + iDim] = Inlet_FlowDir[iMarker][iVertex][iDim]; - } - - /*--- Exit once we find the point. ---*/ - - return Area; - } - } - } - } +su2double CFVMFlowSolverBase::GetInletAtVertex(unsigned short iMarker, unsigned long iVertex, + const CGeometry* geometry, su2double* val_inlet) const { + const auto T_position = nDim; + const auto P_position = nDim + 1; + const auto FlowDir_position = nDim + 2; + val_inlet[T_position] = Inlet_Ttotal[iMarker][iVertex]; + val_inlet[P_position] = Inlet_Ptotal[iMarker][iVertex]; + for (unsigned short iDim = 0; iDim < nDim; iDim++) { + val_inlet[FlowDir_position + iDim] = Inlet_FlowDir[iMarker][iVertex][iDim]; } - /*--- If we don't find a match, then the child point is not on the - current inlet boundary marker. Return zero area so this point does - not contribute to the restriction operator and continue. ---*/ + /*--- Compute boundary face area for this vertex. ---*/ - return Area; + su2double Normal[MAXNDIM] = {0.0}; + geometry->vertex[iMarker][iVertex]->GetNormal(Normal); + return GeometryToolbox::Norm(nDim, Normal); } template void CFVMFlowSolverBase::SetUniformInlet(const CConfig* config, unsigned short iMarker) { if (config->GetMarker_All_KindBC(iMarker) == INLET_FLOW) { - string Marker_Tag = config->GetMarker_All_TagBound(iMarker); - su2double p_total = config->GetInletPtotal(Marker_Tag); - su2double t_total = config->GetInletTtotal(Marker_Tag); - auto flow_dir = config->GetInletFlowDir(Marker_Tag); + const string Marker_Tag = config->GetMarker_All_TagBound(iMarker); + const su2double p_total = config->GetInletPtotal(Marker_Tag); + const su2double t_total = config->GetInletTtotal(Marker_Tag); + const su2double* flow_dir = config->GetInletFlowDir(Marker_Tag); for (unsigned long iVertex = 0; iVertex < nVertex[iMarker]; iVertex++) { Inlet_Ttotal[iMarker][iVertex] = t_total; Inlet_Ptotal[iMarker][iVertex] = p_total; for (unsigned short iDim = 0; iDim < nDim; iDim++) Inlet_FlowDir[iMarker][iVertex][iDim] = flow_dir[iDim]; } + } else if (config->GetMarker_All_KindBC(iMarker) == SUPERSONIC_INLET) { + const string Marker_Tag = config->GetMarker_All_TagBound(iMarker); + const su2double p = config->GetInlet_Pressure(Marker_Tag); + const su2double t = config->GetInlet_Temperature(Marker_Tag); + const su2double* vel = config->GetInlet_Velocity(Marker_Tag); + for (unsigned long iVertex = 0; iVertex < nVertex[iMarker]; iVertex++) { + Inlet_Ttotal[iMarker][iVertex] = t; + Inlet_Ptotal[iMarker][iVertex] = p; + for (unsigned short iDim = 0; iDim < nDim; iDim++) Inlet_FlowDir[iMarker][iVertex][iDim] = vel[iDim]; + } } else { /*--- For now, non-inlets just get set to zero. In the future, we can do more customization for other boundary types here. ---*/ diff --git a/SU2_CFD/include/solvers/CSolver.hpp b/SU2_CFD/include/solvers/CSolver.hpp index 191a61c90ea..02b5b833785 100644 --- a/SU2_CFD/include/solvers/CSolver.hpp +++ b/SU2_CFD/include/solvers/CSolver.hpp @@ -2895,7 +2895,7 @@ class CSolver { * \param[in] config - Definition of the particular problem. * \param[in] iMarker - Surface marker where the coefficient is computed. */ - inline virtual void SetUniformInlet(const CConfig* config, unsigned short iMarker) {}; + inline virtual void SetUniformInlet(const CConfig* config, unsigned short iMarker) {} /*! * \brief A virtual member @@ -2905,23 +2905,18 @@ class CSolver { */ inline virtual void SetInletAtVertex(const su2double *val_inlet, unsigned short iMarker, - unsigned long iVertex) { }; + unsigned long iVertex) { } /*! - * \brief A virtual member - * \param[in] val_inlet - vector returning the inlet values for the current vertex. - * \param[in] val_inlet_point - Node index where the inlet is being set. - * \param[in] val_kind_marker - Enumerated type for the particular inlet type. + * \brief Get the set of values imposed at an inlet. + * \param[in] iMarker - Index of the surface marker. + * \param[in] iVertex - Vertex of the marker iMarker where the inlet is being set. * \param[in] geometry - Geometrical definition of the problem. - * \param config - Definition of the particular problem. + * \param[in,out] val_inlet - vector returning the inlet values for the current vertex. * \return Value of the face area at the vertex. */ - inline virtual su2double GetInletAtVertex(su2double *val_inlet, - unsigned long val_inlet_point, - unsigned short val_kind_marker, - string val_marker, - const CGeometry *geometry, - const CConfig *config) const { return 0; } + inline virtual su2double GetInletAtVertex(unsigned short iMarker, unsigned long iVertex, + const CGeometry* geometry, su2double* val_inlet) const { return 0; } /*! * \brief Update the multi-grid structure for the customized boundary conditions. diff --git a/SU2_CFD/include/solvers/CSpeciesSolver.hpp b/SU2_CFD/include/solvers/CSpeciesSolver.hpp index 4fbb9419117..8b838cd0e61 100644 --- a/SU2_CFD/include/solvers/CSpeciesSolver.hpp +++ b/SU2_CFD/include/solvers/CSpeciesSolver.hpp @@ -113,16 +113,15 @@ class CSpeciesSolver : public CScalarSolver { void SetInletAtVertex(const su2double* val_inlet, unsigned short iMarker, unsigned long iVertex) override; /*! - * \brief Get the set of value imposed at an inlet. - * \param[in] val_inlet - vector returning the inlet values for the current vertex. - * \param[in] val_inlet_point - Node index where the inlet is being set. - * \param[in] val_kind_marker - Enumerated type for the particular inlet type. + * \brief Get the set of values imposed at an inlet. + * \param[in] iMarker - Index of the surface marker. + * \param[in] iVertex - Vertex of the marker iMarker where the inlet is being set. * \param[in] geometry - Geometrical definition of the problem. - * \param config - Definition of the particular problem. + * \param[in,out] val_inlet - vector returning the inlet values for the current vertex. * \return Value of the face area at the vertex. */ - su2double GetInletAtVertex(su2double* val_inlet, unsigned long val_inlet_point, unsigned short val_kind_marker, - string val_marker, const CGeometry* geometry, const CConfig* config) const override; + su2double GetInletAtVertex(unsigned short iMarker, unsigned long iVertex, + const CGeometry* geometry, su2double* val_inlet) const override; /*! * \brief Set a uniform inlet profile diff --git a/SU2_CFD/include/solvers/CTurbSASolver.hpp b/SU2_CFD/include/solvers/CTurbSASolver.hpp index c79649799cc..1bf7d929d0e 100644 --- a/SU2_CFD/include/solvers/CTurbSASolver.hpp +++ b/SU2_CFD/include/solvers/CTurbSASolver.hpp @@ -345,20 +345,15 @@ class CTurbSASolver final : public CTurbSolver { unsigned long iVertex) override; /*! - * \brief Get the set of value imposed at an inlet. - * \param[in] val_inlet - vector returning the inlet values for the current vertex. - * \param[in] val_inlet_point - Node index where the inlet is being set. - * \param[in] val_kind_marker - Enumerated type for the particular inlet type. + * \brief Get the set of values imposed at an inlet. + * \param[in] iMarker - Index of the surface marker. + * \param[in] iVertex - Vertex of the marker iMarker where the inlet is being set. * \param[in] geometry - Geometrical definition of the problem. - * \param config - Definition of the particular problem. + * \param[in,out] val_inlet - vector returning the inlet values for the current vertex. * \return Value of the face area at the vertex. */ - su2double GetInletAtVertex(su2double *val_inlet, - unsigned long val_inlet_point, - unsigned short val_kind_marker, - string val_marker, - const CGeometry *geometry, - const CConfig *config) const override; + su2double GetInletAtVertex(unsigned short iMarker, unsigned long iVertex, + const CGeometry* geometry, su2double* val_inlet) const override; /*! * \brief Set a uniform inlet profile diff --git a/SU2_CFD/include/solvers/CTurbSSTSolver.hpp b/SU2_CFD/include/solvers/CTurbSSTSolver.hpp index bdb4227a191..4ab59ec73a8 100644 --- a/SU2_CFD/include/solvers/CTurbSSTSolver.hpp +++ b/SU2_CFD/include/solvers/CTurbSSTSolver.hpp @@ -269,20 +269,15 @@ class CTurbSSTSolver final : public CTurbSolver { unsigned long iVertex) override; /*! - * \brief Get the set of value imposed at an inlet. - * \param[in] val_inlet - vector returning the inlet values for the current vertex. - * \param[in] val_inlet_point - Node index where the inlet is being set. - * \param[in] val_kind_marker - Enumerated type for the particular inlet type. + * \brief Get the set of values imposed at an inlet. + * \param[in] iMarker - Index of the surface marker. + * \param[in] iVertex - Vertex of the marker iMarker where the inlet is being set. * \param[in] geometry - Geometrical definition of the problem. - * \param config - Definition of the particular problem. + * \param[in,out] val_inlet - vector returning the inlet values for the current vertex. * \return Value of the face area at the vertex. */ - su2double GetInletAtVertex(su2double *val_inlet, - unsigned long val_inlet_point, - unsigned short val_kind_marker, - string val_marker, - const CGeometry *geometry, - const CConfig *config) const override; + su2double GetInletAtVertex(unsigned short iMarker, unsigned long iVertex, + const CGeometry* geometry, su2double* val_inlet) const override; /*! * \brief Set a uniform inlet profile diff --git a/SU2_CFD/src/drivers/CDriver.cpp b/SU2_CFD/src/drivers/CDriver.cpp index 971ae9e0be7..73063e2ed7d 100644 --- a/SU2_CFD/src/drivers/CDriver.cpp +++ b/SU2_CFD/src/drivers/CDriver.cpp @@ -1052,27 +1052,27 @@ void CDriver::PreprocessInlet(CSolver ***solver, CGeometry **geometry, CConfig * /*--- Use LoadInletProfile() routines for the particular solver. ---*/ if (rank == MASTER_NODE) { - cout << endl; - cout << "Reading inlet profile from file: "; - cout << config->GetInlet_FileName() << endl; + cout << "\nReading inlet profile from file: " << config->GetInlet_FileName() << endl; } - if (solver[MESH_0][FLOW_SOL]) { - solver[MESH_0][FLOW_SOL]->LoadInletProfile(geometry, solver, config, val_iter, FLOW_SOL, INLET_FLOW); - } - if (solver[MESH_0][TURB_SOL]) { - solver[MESH_0][TURB_SOL]->LoadInletProfile(geometry, solver, config, val_iter, TURB_SOL, INLET_FLOW); - } - if (solver[MESH_0][SPECIES_SOL]) { - solver[MESH_0][SPECIES_SOL]->LoadInletProfile(geometry, solver, config, val_iter, SPECIES_SOL, INLET_FLOW); + for (const auto marker_type : {INLET_FLOW, SUPERSONIC_INLET}) { + if (solver[MESH_0][FLOW_SOL]) { + solver[MESH_0][FLOW_SOL]->LoadInletProfile(geometry, solver, config, val_iter, FLOW_SOL, marker_type); + } + if (solver[MESH_0][TURB_SOL]) { + solver[MESH_0][TURB_SOL]->LoadInletProfile(geometry, solver, config, val_iter, TURB_SOL, marker_type); + } + if (solver[MESH_0][SPECIES_SOL]) { + solver[MESH_0][SPECIES_SOL]->LoadInletProfile(geometry, solver, config, val_iter, SPECIES_SOL, marker_type); + } } /*--- Exit if profiles were requested for a solver that is not available. ---*/ if (!config->GetFluidProblem()) { - SU2_MPI::Error(string("Inlet profile specification via file (C++) has not been \n") + - string("implemented yet for this solver.\n") + - string("Please set SPECIFIED_INLET_PROFILE= NO and try again."), CURRENT_FUNCTION); + SU2_MPI::Error("Inlet profile specification via file (C++) has not been \n" + "implemented yet for this solver.\n" + "Please set SPECIFIED_INLET_PROFILE= NO and try again.", CURRENT_FUNCTION); } } else { diff --git a/SU2_CFD/src/numerics/flow/convection/hllc.cpp b/SU2_CFD/src/numerics/flow/convection/hllc.cpp index d574f929b20..d94b44759a7 100644 --- a/SU2_CFD/src/numerics/flow/convection/hllc.cpp +++ b/SU2_CFD/src/numerics/flow/convection/hllc.cpp @@ -31,7 +31,7 @@ CUpwHLLC_Flow::CUpwHLLC_Flow(unsigned short val_nDim, unsigned short val_nVar, const CConfig* config) : CNumerics(val_nDim, val_nVar, config) { implicit = (config->GetKind_TimeIntScheme_Flow() == EULER_IMPLICIT); - kappa = config->GetRoe_Kappa(); + /* A grid is defined as dynamic if there's rigid grid movement or grid deformation AND the problem is time domain */ dynamic_grid = config->GetDynamic_Grid(); @@ -114,30 +114,19 @@ CNumerics::ResidualType<> CUpwHLLC_Flow::ComputeResidual(const CConfig* config) Density_j = V_j[nDim+2]; Enthalpy_j = V_j[nDim+3]; - - sq_vel_i = 0.0; - sq_vel_j = 0.0; - - for (iDim = 0; iDim < nDim; iDim++) { - sq_vel_i += Velocity_i[iDim] * Velocity_i[iDim]; - sq_vel_j += Velocity_j[iDim] * Velocity_j[iDim]; - } + sq_vel_i = GeometryToolbox::SquaredNorm(nDim, Velocity_i); + sq_vel_j = GeometryToolbox::SquaredNorm(nDim, Velocity_j); Energy_i = Enthalpy_i - Pressure_i / Density_i; Energy_j = Enthalpy_j - Pressure_j / Density_j; - SoundSpeed_i = sqrt( (Enthalpy_i - 0.5 * sq_vel_i) * Gamma_Minus_One ); - SoundSpeed_j = sqrt( (Enthalpy_j - 0.5 * sq_vel_j) * Gamma_Minus_One ); + SoundSpeed_i = sqrt((Enthalpy_i - 0.5 * sq_vel_i) * Gamma_Minus_One); + SoundSpeed_j = sqrt((Enthalpy_j - 0.5 * sq_vel_j) * Gamma_Minus_One); /*--- Projected velocities ---*/ - ProjVelocity_i = 0; - ProjVelocity_j = 0; - - for (iDim = 0; iDim < nDim; iDim++) { - ProjVelocity_i += Velocity_i[iDim] * UnitNormal[iDim]; - ProjVelocity_j += Velocity_j[iDim] * UnitNormal[iDim]; - } + ProjVelocity_i = GeometryToolbox::DotProduct(nDim, Velocity_i, UnitNormal); + ProjVelocity_j = GeometryToolbox::DotProduct(nDim, Velocity_j, UnitNormal); /*--- Projected Grid Velocity ---*/ @@ -146,7 +135,7 @@ CNumerics::ResidualType<> CUpwHLLC_Flow::ComputeResidual(const CConfig* config) if (dynamic_grid) { for (iDim = 0; iDim < nDim; iDim++) - ProjInterfaceVel += 0.5 * ( GridVel_i[iDim] + GridVel_j[iDim] )*UnitNormal[iDim]; + ProjInterfaceVel += 0.5 * (GridVel_i[iDim] + GridVel_j[iDim]) * UnitNormal[iDim]; SoundSpeed_i -= ProjInterfaceVel; SoundSpeed_j += ProjInterfaceVel; @@ -159,14 +148,11 @@ CNumerics::ResidualType<> CUpwHLLC_Flow::ComputeResidual(const CConfig* config) Rrho = ( sqrt(Density_i) + sqrt(Density_j) ); - sq_velRoe = 0.0; - RoeProjVelocity = - ProjInterfaceVel; - for (iDim = 0; iDim < nDim; iDim++) { RoeVelocity[iDim] = ( Velocity_i[iDim] * sqrt(Density_i) + Velocity_j[iDim] * sqrt(Density_j) ) / Rrho; - sq_velRoe += RoeVelocity[iDim] * RoeVelocity[iDim]; - RoeProjVelocity += RoeVelocity[iDim] * UnitNormal[iDim]; } + sq_velRoe = GeometryToolbox::SquaredNorm(nDim, RoeVelocity); + RoeProjVelocity = GeometryToolbox::DotProduct(nDim, RoeVelocity, UnitNormal) - ProjInterfaceVel; /*--- Mean Roe variables iPoint and jPoint ---*/ @@ -175,10 +161,8 @@ CNumerics::ResidualType<> CUpwHLLC_Flow::ComputeResidual(const CConfig* config) /*--- Roe-averaged speed of sound ---*/ - //RoeSoundSpeed2 = Gamma_Minus_One * ( RoeEnthalpy - 0.5 * sq_velRoe ); RoeSoundSpeed = sqrt( Gamma_Minus_One * ( RoeEnthalpy - 0.5 * sq_velRoe ) ) - ProjInterfaceVel; - /*--- Speed of sound at L and R ---*/ sL = min( RoeProjVelocity - RoeSoundSpeed, ProjVelocity_i - SoundSpeed_i); @@ -214,8 +198,8 @@ CNumerics::ResidualType<> CUpwHLLC_Flow::ComputeResidual(const CConfig* config) IntermediateState[0] = rhoSL * Density_i; for (iDim = 0; iDim < nDim; iDim++) - IntermediateState[iDim+1] = rhoSL * ( Density_i * Velocity_i[iDim] + ( pStar - Pressure_i ) / ( sL - ProjVelocity_i ) * UnitNormal[iDim] ) ; - IntermediateState[nVar-1] = rhoSL * ( Density_i * Energy_i - ( Pressure_i * ProjVelocity_i - pStar * sM) / ( sL - ProjVelocity_i ) ); + IntermediateState[iDim+1] = rhoSL * Density_i * Velocity_i[iDim] + (pStar - Pressure_i) * UnitNormal[iDim] / (sL - sM); + IntermediateState[nVar-1] = rhoSL * Density_i * Energy_i - (Pressure_i * ProjVelocity_i - pStar * sM) / (sL - sM); Flux[0] = sM * IntermediateState[0]; @@ -243,8 +227,8 @@ CNumerics::ResidualType<> CUpwHLLC_Flow::ComputeResidual(const CConfig* config) IntermediateState[0] = rhoSR * Density_j; for (iDim = 0; iDim < nDim; iDim++) - IntermediateState[iDim+1] = rhoSR * ( Density_j * Velocity_j[iDim] + ( pStar - Pressure_j ) / ( sR - ProjVelocity_j ) * UnitNormal[iDim] ) ; - IntermediateState[nVar-1] = rhoSR * ( Density_j * Energy_j - ( Pressure_j * ProjVelocity_j - pStar * sM ) / ( sR - ProjVelocity_j ) ); + IntermediateState[iDim+1] = rhoSR * Density_j * Velocity_j[iDim] + (pStar - Pressure_j) * UnitNormal[iDim] / (sR - sM); + IntermediateState[nVar-1] = rhoSR * Density_j * Energy_j - (Pressure_j * ProjVelocity_j - pStar * sM) / (sR - sM); Flux[0] = sM * IntermediateState[0]; @@ -254,14 +238,12 @@ CNumerics::ResidualType<> CUpwHLLC_Flow::ComputeResidual(const CConfig* config) } } - - for (iVar = 0; iVar < nVar; iVar++) - Flux[iVar] *= Area; + for (iVar = 0; iVar < nVar; iVar++) Flux[iVar] *= Area; /*--- Return early if the Jacobians do not need to be computed. ---*/ - if (implicit) - { + if (!implicit) return ResidualType<>(Flux, Jacobian_i, Jacobian_j); + if (sM > 0.0) { if (sL > 0.0) { @@ -537,27 +519,20 @@ CNumerics::ResidualType<> CUpwHLLC_Flow::ComputeResidual(const CConfig* config) } } - - /*--- Jacobians of the inviscid flux, scale = k because Flux ~ 0.5*(fc_i+fc_j)*Normal ---*/ - - Area *= kappa; - + /*--- Scale Jacobians by area (from Flux *= Area). ---*/ for (iVar = 0; iVar < nVar; iVar++) { for (jVar = 0; jVar < nVar; jVar++) { - Jacobian_i[iVar][jVar] *= Area; - Jacobian_j[iVar][jVar] *= Area; + Jacobian_i[iVar][jVar] *= Area; + Jacobian_j[iVar][jVar] *= Area; } } - } // end if implicit - return ResidualType<>(Flux, Jacobian_i, Jacobian_j); - } CUpwGeneralHLLC_Flow::CUpwGeneralHLLC_Flow(unsigned short val_nDim, unsigned short val_nVar, const CConfig* config) : CNumerics(val_nDim, val_nVar, config) { implicit = (config->GetKind_TimeIntScheme_Flow() == EULER_IMPLICIT); - kappa = config->GetRoe_Kappa(); + /* A grid is defined as dynamic if there's rigid grid movement or grid deformation AND the problem is time domain */ dynamic_grid = config->GetDynamic_Grid(); @@ -792,13 +767,12 @@ CNumerics::ResidualType<> CUpwGeneralHLLC_Flow::ComputeResidual(const CConfig* c } } - for (iVar = 0; iVar < nVar; iVar++) - Flux[iVar] *= Area; + for (iVar = 0; iVar < nVar; iVar++) Flux[iVar] *= Area; /*--- Return early if the Jacobians do not need to be computed. ---*/ - if (implicit) - { + if (!implicit) return ResidualType<>(Flux, Jacobian_i, Jacobian_j); + if (sM > 0.0) { if (sL > 0.0) { @@ -1090,21 +1064,14 @@ CNumerics::ResidualType<> CUpwGeneralHLLC_Flow::ComputeResidual(const CConfig* c } } - - /*--- Jacobians of the inviscid flux, scale = kappa because Flux ~ 0.5*(fc_i+fc_j)*Normal ---*/ - - Area *= kappa; - + /*--- Scale Jacobians by area (from Flux *= Area). ---*/ for (iVar = 0; iVar < nVar; iVar++) { for (jVar = 0; jVar < nVar; jVar++) { Jacobian_i[iVar][jVar] *= Area; Jacobian_j[iVar][jVar] *= Area; } } - } // end if implicit - return ResidualType<>(Flux, Jacobian_i, Jacobian_j); - } void CUpwGeneralHLLC_Flow::VinokurMontagne() { diff --git a/SU2_CFD/src/numerics/flow/convection/roe.cpp b/SU2_CFD/src/numerics/flow/convection/roe.cpp index 5f84fc21858..7b122119009 100644 --- a/SU2_CFD/src/numerics/flow/convection/roe.cpp +++ b/SU2_CFD/src/numerics/flow/convection/roe.cpp @@ -34,7 +34,7 @@ CUpwRoeBase_Flow::CUpwRoeBase_Flow(unsigned short val_nDim, unsigned short val_n implicit = (config->GetKind_TimeIntScheme_Flow() == EULER_IMPLICIT); /* A grid is defined as dynamic if there's rigid grid movement or grid deformation AND the problem is time domain */ dynamic_grid = config->GetDynamic_Grid(); - kappa = config->GetRoe_Kappa(); // 1 is unstable + kappa = config->GetRoe_Kappa(); Gamma = config->GetGamma(); Gamma_Minus_One = Gamma - 1.0; @@ -683,7 +683,7 @@ CUpwGeneralRoe_Flow::CUpwGeneralRoe_Flow(unsigned short val_nDim, unsigned short implicit = (config->GetKind_TimeIntScheme_Flow() == EULER_IMPLICIT); /* A grid is defined as dynamic if there's rigid grid movement or grid deformation AND the problem is time domain */ dynamic_grid = config->GetDynamic_Grid(); - kappa = config->GetRoe_Kappa(); // 1 is unstable + kappa = config->GetRoe_Kappa(); Flux = new su2double [nVar]; diff --git a/SU2_CFD/src/solvers/CEulerSolver.cpp b/SU2_CFD/src/solvers/CEulerSolver.cpp index 356c759564c..6b129998d5c 100644 --- a/SU2_CFD/src/solvers/CEulerSolver.cpp +++ b/SU2_CFD/src/solvers/CEulerSolver.cpp @@ -7343,41 +7343,36 @@ void CEulerSolver::BC_Supersonic_Inlet(CGeometry *geometry, CSolver **solver_con const bool tkeNeeded = (config->GetKind_Turb_Model() == TURB_MODEL::SST); /*--- Supersonic inlet flow: there are no outgoing characteristics, - so all flow variables can be imposed at the inlet. - First, retrieve the specified values for the primitive variables. ---*/ + so all flow variables can be imposed at the inlet. ---*/ - const su2double Temperature = config->GetInlet_Temperature(Marker_Tag) / config->GetTemperature_Ref(); - const su2double Pressure = config->GetInlet_Pressure(Marker_Tag) / config->GetPressure_Ref(); - const auto* Vel = config->GetInlet_Velocity(Marker_Tag); - - su2double Velocity[MAXNDIM] = {0.0}; - for (unsigned short iDim = 0; iDim < nDim; iDim++) - Velocity[iDim] = Vel[iDim] / config->GetVelocity_Ref(); + SU2_OMP_FOR_DYN(OMP_MIN_SIZE) + for (auto iVertex = 0ul; iVertex < geometry->nVertex[val_marker]; iVertex++) { + const auto iPoint = geometry->vertex[val_marker][iVertex]->GetNode(); - /*--- Density at the inlet from the gas law ---*/ + if (!geometry->nodes->GetDomain(iPoint)) continue; - const su2double Density = Pressure / (Gas_Constant * Temperature); + /*--- Retrieve the inlet profile, note that total conditions are reused as static. ---*/ - /*--- Compute the energy from the specified state ---*/ + const su2double Temperature = Inlet_Ttotal[val_marker][iVertex] / config->GetTemperature_Ref(); + const su2double Pressure = Inlet_Ptotal[val_marker][iVertex] / config->GetPressure_Ref(); + su2double Velocity[MAXNDIM] = {0.0}; + for (unsigned short iDim = 0; iDim < nDim; iDim++) { + Velocity[iDim] = Inlet_FlowDir[val_marker][iVertex][iDim] / config->GetVelocity_Ref(); + } - const su2double Velocity2 = GeometryToolbox::SquaredNorm(int(MAXNDIM), Velocity); - su2double Energy = Pressure / (Density * Gamma_Minus_One) + 0.5 * Velocity2; - if (tkeNeeded) Energy += GetTke_Inf(); + /*--- Density at the inlet from the gas law. ---*/ - /*--- Loop over all the vertices on this boundary marker ---*/ + const su2double Density = Pressure / (Gas_Constant * Temperature); - SU2_OMP_FOR_DYN(OMP_MIN_SIZE) - for (auto iVertex = 0ul; iVertex < geometry->nVertex[val_marker]; iVertex++) { - const auto iPoint = geometry->vertex[val_marker][iVertex]->GetNode(); + /*--- Compute the energy from the specified state. ---*/ - if (!geometry->nodes->GetDomain(iPoint)) continue; + const su2double Velocity2 = GeometryToolbox::SquaredNorm(int(MAXNDIM), Velocity); + su2double Energy = Pressure / (Density * Gamma_Minus_One) + 0.5 * Velocity2; + if (tkeNeeded) Energy += GetTke_Inf(); - /*--- Allocate the value at the inlet ---*/ + /*--- Primitive variables, using the derived quantities. ---*/ auto* V_inlet = GetCharacPrimVar(val_marker, iVertex); - - /*--- Primitive variables, using the derived quantities ---*/ - V_inlet[prim_idx.Temperature()] = Temperature; V_inlet[prim_idx.Pressure()] = Pressure; V_inlet[prim_idx.Density()] = Density; @@ -7385,17 +7380,17 @@ void CEulerSolver::BC_Supersonic_Inlet(CGeometry *geometry, CSolver **solver_con for (unsigned short iDim = 0; iDim < nDim; iDim++) V_inlet[iDim+prim_idx.Velocity()] = Velocity[iDim]; - /*--- Current solution at this boundary node ---*/ + /*--- Current solution at this boundary node. ---*/ - auto* V_domain = nodes->GetPrimitive(iPoint); + const auto* V_domain = nodes->GetPrimitive(iPoint); - /*--- Normal vector for this vertex (negate for outward convention) ---*/ + /*--- Normal vector for this vertex (negate for outward convention). ---*/ su2double Normal[MAXNDIM] = {0.0}; geometry->vertex[val_marker][iVertex]->GetNormal(Normal); for (unsigned short iDim = 0; iDim < nDim; iDim++) Normal[iDim] = -Normal[iDim]; - /*--- Set various quantities in the solver class ---*/ + /*--- Set various quantities in the solver class. ---*/ conv_numerics->SetNormal(Normal); conv_numerics->SetPrimitive(V_domain, V_inlet); @@ -7404,13 +7399,13 @@ void CEulerSolver::BC_Supersonic_Inlet(CGeometry *geometry, CSolver **solver_con conv_numerics->SetGridVel(geometry->nodes->GetGridVel(iPoint), geometry->nodes->GetGridVel(iPoint)); - /*--- Compute the residual using an upwind scheme ---*/ + /*--- Compute the residual using an upwind scheme. ---*/ auto residual = conv_numerics->ComputeResidual(config); LinSysRes.AddBlock(iPoint, residual); - /*--- Jacobian contribution for implicit integration ---*/ + /*--- Jacobian contribution for implicit integration. ---*/ if (implicit) Jacobian.AddBlock2Diag(iPoint, residual.jacobian_i); diff --git a/SU2_CFD/src/solvers/CSolver.cpp b/SU2_CFD/src/solvers/CSolver.cpp index de0ec914b73..bfbd4edad58 100644 --- a/SU2_CFD/src/solvers/CSolver.cpp +++ b/SU2_CFD/src/solvers/CSolver.cpp @@ -3588,15 +3588,26 @@ void CSolver::LoadInletProfile(CGeometry **geometry, /*--- Skip if this is the wrong type of marker. ---*/ if (config->GetMarker_All_KindBC(iMarker) != KIND_MARKER) continue; - string Marker_Tag = config->GetMarker_All_TagBound(iMarker); - su2double p_total = config->GetInletPtotal(Marker_Tag); - su2double t_total = config->GetInletTtotal(Marker_Tag); - auto flow_dir = config->GetInletFlowDir(Marker_Tag); - std::stringstream columnName,columnValue; + const string Marker_Tag = config->GetMarker_All_TagBound(iMarker); + std::stringstream columnName,columnValue; columnValue << setprecision(15); columnValue << std::scientific; + su2double p_total, t_total; + const su2double* flow_dir = nullptr; + + if (KIND_MARKER == INLET_FLOW) { + p_total = config->GetInletPtotal(Marker_Tag); + t_total = config->GetInletTtotal(Marker_Tag); + flow_dir = config->GetInletFlowDir(Marker_Tag); + } else if (KIND_MARKER == SUPERSONIC_INLET) { + p_total = config->GetInlet_Pressure(Marker_Tag); + t_total = config->GetInlet_Temperature(Marker_Tag); + flow_dir = config->GetInlet_Velocity(Marker_Tag); + } else { + SU2_MPI::Error("Unsupported type of inlet.", CURRENT_FUNCTION); + } columnValue << t_total << "\t" << p_total <<"\t"; for (unsigned short iDim = 0; iDim < nDim; iDim++) { columnValue << flow_dir[iDim] <<"\t"; @@ -3605,7 +3616,9 @@ void CSolver::LoadInletProfile(CGeometry **geometry, columnName << "# COORD-X " << setw(24) << "COORD-Y " << setw(24); if(nDim==3) columnName << "COORD-Z " << setw(24); - if (config->GetKind_Regime()==ENUM_REGIME::COMPRESSIBLE){ + if (KIND_MARKER == SUPERSONIC_INLET) { + columnName << "TEMPERATURE" << setw(24) << "PRESSURE " << setw(24); + } else if (config->GetKind_Regime() == ENUM_REGIME::COMPRESSIBLE) { switch (config->GetKind_Inlet()) { /*--- compressible conditions ---*/ case INLET_TYPE::TOTAL_CONDITIONS: @@ -3616,7 +3629,8 @@ void CSolver::LoadInletProfile(CGeometry **geometry, break; default: SU2_MPI::Error("Unsupported INLET_TYPE.", CURRENT_FUNCTION); - break; } + break; + } } else { switch (config->GetKind_Inc_Inlet(Marker_Tag)) { /*--- incompressible conditions ---*/ @@ -3674,10 +3688,16 @@ void CSolver::LoadInletProfile(CGeometry **geometry, } + /*--- There are no markers of this type. ---*/ + + const unsigned short has_names = !columnNames.empty(); + unsigned short any_has_names; + SU2_MPI::Allreduce(&has_names, &any_has_names, 1, MPI_UNSIGNED_SHORT, MPI_MAX, SU2_MPI::GetComm()); + if (!any_has_names) return; /*--- Read the profile data from an ASCII file. ---*/ - CMarkerProfileReaderFVM profileReader(geometry[MESH_0], config, profile_filename, KIND_MARKER, nCol_InletFile, columnNames,columnValues); + CMarkerProfileReaderFVM profileReader(geometry[MESH_0], config, profile_filename, KIND_MARKER, nCol_InletFile, columnNames, columnValues); /*--- Load data from the restart into correct containers. ---*/ @@ -3913,7 +3933,7 @@ void CSolver::LoadInletProfile(CGeometry **geometry, const auto Marker_Tag = config->GetMarker_All_TagBound(iMarker); - /* Check the number of columns and allocate temp array. */ + /*--- Check the number of columns and allocate temp array. ---*/ unsigned short nColumns = 0; for (auto jMarker = 0ul; jMarker < profileReader.GetNumberOfProfiles(); jMarker++) { @@ -3947,12 +3967,10 @@ void CSolver::LoadInletProfile(CGeometry **geometry, the averaging. ---*/ for (auto iChildren = 0u; iChildren < geometry[iMesh]->nodes->GetnChildren_CV(iPoint); iChildren++) { - const auto Point_Fine = geometry[iMesh]->nodes->GetChildren_CV(iPoint, iChildren); - - auto Area_Children = solver[iMesh-1][KIND_SOLVER]->GetInletAtVertex(Inlet_Fine.data(), Point_Fine, KIND_MARKER, - Marker_Tag, geometry[iMesh-1], config); + const auto Area_Children = + solver[iMesh-1][KIND_SOLVER]->GetInletAtVertex(iMarker, iVertex, geometry[iMesh-1], Inlet_Fine.data()); for (auto iVar = 0u; iVar < nColumns; iVar++) - Inlet_Values[iVar] += Inlet_Fine[iVar]*Area_Children/Area_Parent; + Inlet_Values[iVar] += Inlet_Fine[iVar] * Area_Children / Area_Parent; } /*--- Set the boundary area-averaged inlet values for the coarse point. ---*/ diff --git a/SU2_CFD/src/solvers/CSpeciesSolver.cpp b/SU2_CFD/src/solvers/CSpeciesSolver.cpp index 3c7ac895708..832c9014511 100644 --- a/SU2_CFD/src/solvers/CSpeciesSolver.cpp +++ b/SU2_CFD/src/solvers/CSpeciesSolver.cpp @@ -415,59 +415,16 @@ void CSpeciesSolver::SetInletAtVertex(const su2double *val_inlet, } -su2double CSpeciesSolver::GetInletAtVertex(su2double *val_inlet, - unsigned long val_inlet_point, - unsigned short val_kind_marker, - string val_marker, - const CGeometry *geometry, - const CConfig *config) const { - /*--- Local variables ---*/ - - unsigned short iMarker; - unsigned long iPoint, iVertex; - su2double Area = 0.0; - su2double Normal[3] = {0.0,0.0,0.0}; - - /*--- Alias positions within inlet file for readability ---*/ - - if (val_kind_marker == INLET_FLOW) { - - for (iMarker = 0; iMarker < config->GetnMarker_All(); iMarker++) { - if ((config->GetMarker_All_KindBC(iMarker) == INLET_FLOW) && - (config->GetMarker_All_TagBound(iMarker) == val_marker)) { - - for (iVertex = 0; iVertex < nVertex[iMarker]; iVertex++) { - - iPoint = geometry->vertex[iMarker][iVertex]->GetNode(); - - if (iPoint == val_inlet_point) { - - /*-- Compute boundary face area for this vertex. ---*/ - - geometry->vertex[iMarker][iVertex]->GetNormal(Normal); - Area = GeometryToolbox::Norm(nDim, Normal); - - /*--- Access and store the inlet variables for this vertex. ---*/ - for (unsigned short iVar = 0; iVar < nVar; iVar++) - val_inlet[Inlet_Position + iVar] = Inlet_SpeciesVars[iMarker][iVertex][iVar]; - - /*--- Exit once we find the point. ---*/ - - return Area; - - } - } - } - } - - } - - /*--- If we don't find a match, then the child point is not on the - current inlet boundary marker. Return zero area so this point does - not contribute to the restriction operator and continue. ---*/ +su2double CSpeciesSolver::GetInletAtVertex(unsigned short iMarker, unsigned long iVertex, + const CGeometry* geometry, su2double* val_inlet) const { + for (unsigned short iVar = 0; iVar < nVar; iVar++) + val_inlet[Inlet_Position + iVar] = Inlet_SpeciesVars[iMarker][iVertex][iVar]; - return Area; + /*--- Compute boundary face area for this vertex. ---*/ + su2double Normal[MAXNDIM] = {0.0}; + geometry->vertex[iMarker][iVertex]->GetNormal(Normal); + return GeometryToolbox::Norm(nDim, Normal); } void CSpeciesSolver::SetUniformInlet(const CConfig* config, unsigned short iMarker) { diff --git a/SU2_CFD/src/solvers/CTurbSASolver.cpp b/SU2_CFD/src/solvers/CTurbSASolver.cpp index f53596e9901..40e19bd024c 100644 --- a/SU2_CFD/src/solvers/CTurbSASolver.cpp +++ b/SU2_CFD/src/solvers/CTurbSASolver.cpp @@ -308,6 +308,8 @@ void CTurbSASolver::Viscous_Residual(const unsigned long iEdge, const CGeometry* void CTurbSASolver::Source_Residual(CGeometry *geometry, CSolver **solver_container, CNumerics **numerics_container, CConfig *config, unsigned short iMesh) { + bool axisymmetric = config->GetAxisymmetric(); + const bool implicit = (config->GetKind_TimeIntScheme() == EULER_IMPLICIT); const bool harmonic_balance = (config->GetTime_Marching() == TIME_MARCHING::HARMONIC_BALANCE); const bool transition_BC = config->GetSAParsedOptions().bc; @@ -383,6 +385,11 @@ void CTurbSASolver::Source_Residual(CGeometry *geometry, CSolver **solver_contai numerics->SetIntermittency(solver_container[TRANS_SOL]->GetNodes()->GetSolution(iPoint, 0)); } + if (axisymmetric) { + /*--- Set y coordinate ---*/ + numerics->SetCoord(geometry->nodes->GetCoord(iPoint), geometry->nodes->GetCoord(iPoint)); + } + /*--- Compute the source term ---*/ auto residual = numerics->ComputeResidual(config); @@ -1521,61 +1528,16 @@ void CTurbSASolver::SetInletAtVertex(const su2double *val_inlet, } -su2double CTurbSASolver::GetInletAtVertex(su2double *val_inlet, - unsigned long val_inlet_point, - unsigned short val_kind_marker, - string val_marker, - const CGeometry *geometry, - const CConfig *config) const { - /*--- Local variables ---*/ - - unsigned short iMarker; - unsigned long iPoint, iVertex; - su2double Area = 0.0; - su2double Normal[3] = {0.0,0.0,0.0}; - - /*--- Alias positions within inlet file for readability ---*/ - - if (val_kind_marker == INLET_FLOW) { - - unsigned short position = nDim+2+nDim; - - for (iMarker = 0; iMarker < config->GetnMarker_All(); iMarker++) { - if ((config->GetMarker_All_KindBC(iMarker) == INLET_FLOW) && - (config->GetMarker_All_TagBound(iMarker) == val_marker)) { - - for (iVertex = 0; iVertex < nVertex[iMarker]; iVertex++){ - - iPoint = geometry->vertex[iMarker][iVertex]->GetNode(); - - if (iPoint == val_inlet_point) { - - /*-- Compute boundary face area for this vertex. ---*/ - - geometry->vertex[iMarker][iVertex]->GetNormal(Normal); - Area = GeometryToolbox::Norm(nDim, Normal); - - /*--- Access and store the inlet variables for this vertex. ---*/ - - val_inlet[position] = Inlet_TurbVars[iMarker][iVertex][0]; - - /*--- Exit once we find the point. ---*/ - - return Area; - - } - } - } - } - - } - - /*--- If we don't find a match, then the child point is not on the - current inlet boundary marker. Return zero area so this point does - not contribute to the restriction operator and continue. ---*/ +su2double CTurbSASolver::GetInletAtVertex(unsigned short iMarker, unsigned long iVertex, + const CGeometry* geometry, su2double* val_inlet) const { + const auto position = nDim + 2 + nDim; + val_inlet[position] = Inlet_TurbVars[iMarker][iVertex][0]; - return Area; + /*--- Compute boundary face area for this vertex. ---*/ + su2double Normal[MAXNDIM] = {0.0}; + geometry->vertex[iMarker][iVertex]->GetNormal(Normal); + return GeometryToolbox::Norm(nDim, Normal); } void CTurbSASolver::SetUniformInlet(const CConfig* config, unsigned short iMarker) { diff --git a/SU2_CFD/src/solvers/CTurbSSTSolver.cpp b/SU2_CFD/src/solvers/CTurbSSTSolver.cpp index ff14596bdfd..63c29102175 100644 --- a/SU2_CFD/src/solvers/CTurbSSTSolver.cpp +++ b/SU2_CFD/src/solvers/CTurbSSTSolver.cpp @@ -1001,63 +1001,18 @@ void CTurbSSTSolver::SetInletAtVertex(const su2double *val_inlet, } -su2double CTurbSSTSolver::GetInletAtVertex(su2double *val_inlet, - unsigned long val_inlet_point, - unsigned short val_kind_marker, - string val_marker, - const CGeometry *geometry, - const CConfig *config) const { - /*--- Local variables ---*/ - - unsigned short iMarker; - unsigned long iPoint, iVertex; - su2double Area = 0.0; - su2double Normal[3] = {0.0,0.0,0.0}; - - /*--- Alias positions within inlet file for readability ---*/ - - if (val_kind_marker == INLET_FLOW) { - - unsigned short tke_position = nDim+2+nDim; - unsigned short omega_position = nDim+2+nDim+1; - - for (iMarker = 0; iMarker < config->GetnMarker_All(); iMarker++) { - if ((config->GetMarker_All_KindBC(iMarker) == INLET_FLOW) && - (config->GetMarker_All_TagBound(iMarker) == val_marker)) { - - for (iVertex = 0; iVertex < nVertex[iMarker]; iVertex++){ - - iPoint = geometry->vertex[iMarker][iVertex]->GetNode(); - - if (iPoint == val_inlet_point) { - - /*-- Compute boundary face area for this vertex. ---*/ - - geometry->vertex[iMarker][iVertex]->GetNormal(Normal); - Area = GeometryToolbox::Norm(nDim, Normal); - - /*--- Access and store the inlet variables for this vertex. ---*/ - - val_inlet[tke_position] = Inlet_TurbVars[iMarker][iVertex][0]; - val_inlet[omega_position] = Inlet_TurbVars[iMarker][iVertex][1]; - - /*--- Exit once we find the point. ---*/ - - return Area; - - } - } - } - } - - } - - /*--- If we don't find a match, then the child point is not on the - current inlet boundary marker. Return zero area so this point does - not contribute to the restriction operator and continue. ---*/ - - return Area; - +su2double CTurbSSTSolver::GetInletAtVertex(unsigned short iMarker, unsigned long iVertex, + const CGeometry* geometry, su2double* val_inlet) const { + const auto tke_position = nDim + 2 + nDim; + const auto omega_position = tke_position + 1; + val_inlet[tke_position] = Inlet_TurbVars[iMarker][iVertex][0]; + val_inlet[omega_position] = Inlet_TurbVars[iMarker][iVertex][1]; + + /*--- Compute boundary face area for this vertex. ---*/ + + su2double Normal[MAXNDIM] = {0.0}; + geometry->vertex[iMarker][iVertex]->GetNormal(Normal); + return GeometryToolbox::Norm(nDim, Normal); } void CTurbSSTSolver::SetUniformInlet(const CConfig* config, unsigned short iMarker) { diff --git a/TestCases/hybrid_regression.py b/TestCases/hybrid_regression.py index 3f14c4c40c6..c893649c1dd 100644 --- a/TestCases/hybrid_regression.py +++ b/TestCases/hybrid_regression.py @@ -67,7 +67,7 @@ def main(): wedge.cfg_dir = "euler/wedge" wedge.cfg_file = "inv_wedge_HLLC.cfg" wedge.test_iter = 20 - wedge.test_vals = [-1.396962, 4.262003, -0.244219, 0.043052] + wedge.test_vals = [-1.368091, 4.302736, -0.243433, 0.042906] test_list.append(wedge) # ONERA M6 Wing @@ -121,8 +121,8 @@ def main(): cylinder_lowmach.cfg_dir = "navierstokes/cylinder" cylinder_lowmach.cfg_file = "cylinder_lowmach.cfg" cylinder_lowmach.test_iter = 25 - cylinder_lowmach.test_vals = [-6.850130, -1.388096, -0.056036, 108.140809, 0.007988] - cylinder_lowmach.test_vals_aarch64 = [-6.850130, -1.388096, -0.056036, 108.140813, 0.007988] + cylinder_lowmach.test_vals = [-6.830996, -1.368850, -0.143956, 73.963354, 0.002457] + cylinder_lowmach.test_vals_aarch64 = [-6.830996, -1.368850, -0.143956, 73.963354, 0.002457] test_list.append(cylinder_lowmach) # 2D Poiseuille flow (body force driven with periodic inlet / outlet) @@ -139,8 +139,8 @@ def main(): poiseuille_profile.cfg_dir = "navierstokes/poiseuille" poiseuille_profile.cfg_file = "profile_poiseuille.cfg" poiseuille_profile.test_iter = 10 - poiseuille_profile.test_vals = [-12.494728, -7.712546, -0.000000, 2.085796] - poiseuille_profile.test_vals_aarch64 = [-12.494717, -7.711274, -0.000000, 2.085796] + poiseuille_profile.test_vals = [-12.485957, -7.612048, -0.000000, 2.085796] + poiseuille_profile.test_vals_aarch64 = [-12.485957, -7.612048, -0.000000, 2.085796] test_list.append(poiseuille_profile) # 2D Rotational Periodic diff --git a/TestCases/parallel_regression.py b/TestCases/parallel_regression.py index 4cad2c9789c..03e6fbaef47 100644 --- a/TestCases/parallel_regression.py +++ b/TestCases/parallel_regression.py @@ -236,7 +236,7 @@ def main(): wedge.cfg_dir = "euler/wedge" wedge.cfg_file = "inv_wedge_HLLC.cfg" wedge.test_iter = 20 - wedge.test_vals = [-1.406716, 4.253025, -0.244411, 0.043089] + wedge.test_vals = [-1.377543, 4.293870, -0.243566, 0.042930] test_list.append(wedge) # ONERA M6 Wing @@ -282,7 +282,7 @@ def main(): ea_naca64206.cfg_dir = "optimization_euler/equivalentarea_naca64206" ea_naca64206.cfg_file = "NACA64206.cfg" ea_naca64206.test_iter = 10 - ea_naca64206.test_vals = [-1.151334, -0.455927, -0.003879, 67775.000000] + ea_naca64206.test_vals = [-1.188459, -0.522783, -0.003147, 67775.000000] test_list.append(ea_naca64206) # SUPERSONIC FLOW PAST A RAMP IN A CHANNEL @@ -327,7 +327,7 @@ def main(): cylinder_lowmach.cfg_dir = "navierstokes/cylinder" cylinder_lowmach.cfg_file = "cylinder_lowmach.cfg" cylinder_lowmach.test_iter = 25 - cylinder_lowmach.test_vals = [-6.858484, -1.396528, -1.854558, 110.033249, 0.001951] + cylinder_lowmach.test_vals = [-6.841604, -1.379532, -1.266739, 76.118218, 0.000274] test_list.append(cylinder_lowmach) # 2D Poiseuille flow (body force driven with periodic inlet / outlet) @@ -344,8 +344,8 @@ def main(): poiseuille_profile.cfg_dir = "navierstokes/poiseuille" poiseuille_profile.cfg_file = "profile_poiseuille.cfg" poiseuille_profile.test_iter = 10 - poiseuille_profile.test_vals = [-12.492939, -7.672950, -0.000000, 2.085796] - poiseuille_profile.test_vals_aarch64 = [-12.492864, -7.671632, -0.000000, 2.085796] + poiseuille_profile.test_vals = [-12.483967, -7.577331, -0.000000, 2.085796] + poiseuille_profile.test_vals_aarch64 = [-12.483967, -7.577331, -0.000000, 2.085796] test_list.append(poiseuille_profile) ########################## @@ -389,18 +389,18 @@ def main(): turb_flatplate_species.cfg_dir = "rans/flatplate" turb_flatplate_species.cfg_file = "turb_SA_flatplate_species.cfg" turb_flatplate_species.test_iter = 20 - turb_flatplate_species.test_vals = [-4.147387, -0.634805, -1.769879, 1.335329, -3.250300, 9.000000, -6.694977, 5.000000, -6.991166, 10.000000, -6.033066, 0.996034, 0.996034] + turb_flatplate_species.test_vals = [-4.120225, -0.634021, -1.706720, 1.363240, -3.250204, 9.000000, -6.697079, 5.000000, -6.978731, 10.000000, -6.013196, 0.996237, 0.996237] test_list.append(turb_flatplate_species) - # Flat plate SST compressibility correction Wilcox + # Flat plate SST compressibility correction Wilcox turb_flatplate_CC_Wilcox = TestCase('turb_flatplate_CC_Wilcox') turb_flatplate_CC_Wilcox.cfg_dir = "rans/flatplate" turb_flatplate_CC_Wilcox.cfg_file = "turb_SST_flatplate_compressibility_Wilcox.cfg" turb_flatplate_CC_Wilcox.test_iter = 20 turb_flatplate_CC_Wilcox.test_vals = [-1.280875, 1.974212, 1.440458, 5.038402, -4.051125, 8.521136] test_list.append(turb_flatplate_CC_Wilcox) - - # Flat plate SST compressibility correction Sarkar + + # Flat plate SST compressibility correction Sarkar turb_flatplate_CC_Sarkar = TestCase('turb_flatplate_CC_Sarkar') turb_flatplate_CC_Sarkar.cfg_dir = "rans/flatplate" turb_flatplate_CC_Sarkar.cfg_file = "turb_SST_flatplate_compressibility_Sarkar.cfg" diff --git a/TestCases/parallel_regression_AD.py b/TestCases/parallel_regression_AD.py index b62fed781a3..a4aa0f57252 100644 --- a/TestCases/parallel_regression_AD.py +++ b/TestCases/parallel_regression_AD.py @@ -71,7 +71,7 @@ def main(): ea_naca64206.cfg_dir = "optimization_euler/equivalentarea_naca64206" ea_naca64206.cfg_file = "NACA64206.cfg" ea_naca64206.test_iter = 10 - ea_naca64206.test_vals = [3.182170, 2.473052, -5509000.000000, 5.551800] + ea_naca64206.test_vals = [3.127605, 2.411805, -5505700.000000, 10.591000] test_list.append(ea_naca64206) #################################### diff --git a/TestCases/serial_regression.py b/TestCases/serial_regression.py index 8d328dc460f..080b41053b7 100644 --- a/TestCases/serial_regression.py +++ b/TestCases/serial_regression.py @@ -109,7 +109,7 @@ def main(): wedge.cfg_dir = "euler/wedge" wedge.cfg_file = "inv_wedge_HLLC.cfg" wedge.test_iter = 20 - wedge.test_vals = [-1.396962, 4.262003, -0.244219, 0.043052] + wedge.test_vals = [-1.368091, 4.302736, -0.243433, 0.042906] test_list.append(wedge) # ONERA M6 Wing @@ -182,7 +182,7 @@ def main(): cylinder_lowmach.cfg_dir = "navierstokes/cylinder" cylinder_lowmach.cfg_file = "cylinder_lowmach.cfg" cylinder_lowmach.test_iter = 25 - cylinder_lowmach.test_vals = [-6.850123, -1.388088, -0.056090, 108.140176, 0.007983] + cylinder_lowmach.test_vals = [-6.830989, -1.368842, -0.143838, 73.962440, 0.002454] test_list.append(cylinder_lowmach) # 2D Poiseuille flow (body force driven with periodic inlet / outlet) @@ -198,8 +198,8 @@ def main(): poiseuille_profile.cfg_dir = "navierstokes/poiseuille" poiseuille_profile.cfg_file = "profile_poiseuille.cfg" poiseuille_profile.test_iter = 10 - poiseuille_profile.test_vals = [-12.494681, -7.711642, -0.000000, 2.085796] - poiseuille_profile.test_vals_aarch64 = [-12.494684, -7.711379, -0.000000, 2.085796] #last 4 columns + poiseuille_profile.test_vals = [-12.485974, -7.612341, -0.000000, 2.085796] + poiseuille_profile.test_vals_aarch64 = [-12.485974, -7.612341, -0.000000, 2.085796] test_list.append(poiseuille_profile) ########################## diff --git a/TestCases/tutorials.py b/TestCases/tutorials.py index e2cac1f32d0..5d987f8cd9a 100644 --- a/TestCases/tutorials.py +++ b/TestCases/tutorials.py @@ -145,7 +145,7 @@ def main(): tutorial_inv_wedge.cfg_dir = "../Tutorials/compressible_flow/Inviscid_Wedge" tutorial_inv_wedge.cfg_file = "inv_wedge_HLLC.cfg" tutorial_inv_wedge.test_iter = 0 - tutorial_inv_wedge.test_vals = [-0.864206, 4.850246, -0.259185, 0.045567] + tutorial_inv_wedge.test_vals = [-0.864206, 4.850246, -0.245674, 0.043209] tutorial_inv_wedge.no_restart = True test_list.append(tutorial_inv_wedge) diff --git a/config_template.cfg b/config_template.cfg index ecf00315bd8..00334732b30 100644 --- a/config_template.cfg +++ b/config_template.cfg @@ -1567,7 +1567,7 @@ CONV_NUM_METHOD_FLOW= ROE % Roe Low Dissipation function for Hybrid RANS/LES simulations (FD, NTS, NTS_DUCROS) ROE_LOW_DISSIPATION= FD % -% Roe coefficient +% Roe dissipation coefficient ROE_KAPPA= 0.5 % % Minimum value for beta for the Roe-Turkel preconditioner