From c7916b5203b919787b9dd37c979f6e1976efda26 Mon Sep 17 00:00:00 2001 From: rois1995 Date: Mon, 12 Oct 2026 11:33:00 +0200 Subject: [PATCH 1/9] Remove commented-out code from the NEMO solvers Delete the disabled viscous blocks in BC_Inlet and BC_Outlet, the disabled chemistry source block and the leftover species loop in the heat-flux wall, and the unused turbulent kinetic energy lines. No functional change. Co-Authored-By: Claude Opus 5.5 Claude-Session: https://claude.ai/code/session_01L7oxzVir4zENrpy3K3W6Fm --- SU2_CFD/src/solvers/CNEMOEulerSolver.cpp | 80 ------------------------ SU2_CFD/src/solvers/CNEMONSSolver.cpp | 33 +--------- 2 files changed, 1 insertion(+), 112 deletions(-) diff --git a/SU2_CFD/src/solvers/CNEMOEulerSolver.cpp b/SU2_CFD/src/solvers/CNEMOEulerSolver.cpp index 2a09a7a1f42e..b158b525688c 100644 --- a/SU2_CFD/src/solvers/CNEMOEulerSolver.cpp +++ b/SU2_CFD/src/solvers/CNEMOEulerSolver.cpp @@ -1788,10 +1788,6 @@ void CNEMOEulerSolver::BC_Inlet(CGeometry *geometry, CSolver **solver_container, V_inlet[nSpecies+2] = Velocity[iDim]; V_inlet[nSpecies+nDim+2] = Pressure; V_inlet[RHO_INDEX] = Density; - //V_inlet[H_INDEX] = H; - //V_inlet[A_INDEX] = A; - //V_inlet[RHO_CVTR_INDEX] = rcvtr; - //V_inlet[RHO_CVVE_INDEX] = rcvve; break; @@ -1865,32 +1861,6 @@ void CNEMOEulerSolver::BC_Inlet(CGeometry *geometry, CSolver **solver_container, /*--- Jacobian contribution for implicit integration ---*/ if (implicit) Jacobian.AddBlock2Diag(iPoint, residual.jacobian_i); - -// /*--- Viscous contribution ---*/ -// if (viscous) { - -// /*--- Set the normal vector and the coordinates ---*/ -// visc_numerics->SetNormal(Normal); -// su2double Coord_Reflected[MAXNDIM]; -// GeometryToolbox::PointPointReflect(nDim, geometry->nodes->GetCoord(Point_Normal), -// geometry->nodes->GetCoord(iPoint), Coord_Reflected); -// visc_numerics->SetCoord(geometry->nodes->GetCoord(), Coord_Reflected); - -// /*--- Primitive variables, and gradient ---*/ -// visc_numerics->SetPrimitive(V_domain, V_inlet); -// visc_numerics->SetPrimVarGradient(nodes->GetGradient_Primitive(), nodes->GetGradient_Primitive()); - -// /*--- Laminar viscosity ---*/ -// visc_numerics->SetLaminarViscosity(nodes->GetLaminarViscosity(), nodes->GetLaminarViscosity()); - -// /*--- Compute and update residual ---*/ -// visc_numerics->ComputeResidual(Residual, Jacobian_i, Jacobian_j, config); -// LinSysRes.SubtractBlock(iPoint, Residual); - -// /*--- Jacobian contribution for implicit integration ---*/ -// if (implicit) -// Jacobian.SubtractBlock2Diag(iPoint, Jacobian_i); -// } } } @@ -2095,56 +2065,6 @@ void CNEMOEulerSolver::BC_Outlet(CGeometry *geometry, CSolver **solver_container /*--- Jacobian contribution for implicit integration ---*/ if (implicit) Jacobian.AddBlock2Diag(iPoint, residual.jacobian_i); - - /*--- Viscous contribution ---*/ -// if (viscous) { - -// /*--- Set the normal vector and the coordinates ---*/ -// visc_numerics->SetNormal(Normal); -// su2double Coord_Reflected[MAXNDIM]; -// GeometryToolbox::PointPointReflect(nDim, geometry->nodes->GetCoord(Point_Normal), -// geometry->nodes->GetCoord(iPoint), Coord_Reflected); -// visc_numerics->SetCoord(geometry->nodes->GetCoord(), Coord_Reflected); - -// /*--- Primitive variables, and gradient ---*/ -// visc_numerics->SetPrimitive(V_domain, V_outlet); -// visc_numerics->SetPrimVarGradient(nodes->GetGradient_Primitive(), nodes->GetGradient_Primitive()); - -// /*--- Conservative variables, and gradient ---*/ -// visc_numerics->SetConservative(U_domain, U_outlet); - -// /*--- Pass supplementary information to CNumerics ---*/ -// visc_numerics->SetdPdU(nodes->GetdPdU(), node_infty->GetdPdU()); -// visc_numerics->SetdTdU(nodes->GetdTdU(), node_infty->GetdTdU()); -// visc_numerics->SetdTvedU(nodes->GetdTvedU(), node_infty->GetdTvedU()); - -// /*--- Species diffusion coefficients ---*/ -// visc_numerics->SetDiffusionCoeff(nodes->GetDiffusionCoeff(), -// node_infty->GetDiffusionCoeff() ); - -// /*--- Laminar viscosity ---*/ -// visc_numerics->SetLaminarViscosity(nodes->GetLaminarViscosity(), -// node_infty->GetLaminarViscosity() ); - -// /*--- Thermal conductivity ---*/ -// visc_numerics->SetThermalConductivity(nodes->GetThermalConductivity(), -// node_infty->GetThermalConductivity()); - -// /*--- Vib-el. thermal conductivity ---*/ -// visc_numerics->SetThermalConductivity_ve(nodes->GetThermalConductivity_ve(), -// node_infty->GetThermalConductivity_ve() ); - -// /*--- Laminar viscosity ---*/ -// visc_numerics->SetLaminarViscosity(nodes->GetLaminarViscosity(), nodes->GetLaminarViscosity()); - -// /*--- Compute and update residual ---*/ -// visc_numerics->ComputeResidual(Residual, Jacobian_i, Jacobian_j, config); -// LinSysRes.SubtractBlock(iPoint, Residual); - -// /*--- Jacobian contribution for implicit integration ---*/ -// if (implicit) -// Jacobian.SubtractBlock2Diag(iPoint, residual.jacobian_i); -// } } } diff --git a/SU2_CFD/src/solvers/CNEMONSSolver.cpp b/SU2_CFD/src/solvers/CNEMONSSolver.cpp index e29934e88032..513344fcfc31 100644 --- a/SU2_CFD/src/solvers/CNEMONSSolver.cpp +++ b/SU2_CFD/src/solvers/CNEMONSSolver.cpp @@ -117,18 +117,16 @@ unsigned long CNEMONSSolver::SetPrimitive_Variables(CSolver **solver_container,C unsigned long nonPhysicalPoints = 0; const TURB_MODEL turb_model = config->GetKind_Turb_Model(); - //const bool tkeNeeded = (turb_model == TURB_MODEL::SST); SU2_OMP_FOR_STAT(omp_chunk_size) for (auto iPoint = 0ul; iPoint < nPoint; iPoint ++) { /*--- Retrieve the value of the kinetic energy (if needed). ---*/ - su2double eddy_visc = 0.0; //su2double turb_ke = 0.0; + su2double eddy_visc = 0.0; if (turb_model != TURB_MODEL::NONE && solver_container[TURB_SOL] != nullptr) { eddy_visc = solver_container[TURB_SOL]->GetNodes()->GetmuT(iPoint); - //if (tkeNeeded) turb_ke = solver_container[TURB_SOL]->GetNodes()->GetSolution(iPoint,0); nodes->SetEddyViscosity(iPoint, eddy_visc); } @@ -449,15 +447,6 @@ void CNEMONSSolver::BC_HeatFluxCatalytic_Wall(CGeometry *geometry, SdYdn = 0.0; for (auto iSpecies = 0ul; iSpecies < nSpecies; iSpecies++) SdYdn += rho*Ds[iSpecies]*dYdn[iSpecies]; - - for (auto iSpecies = 0ul; iSpecies < nSpecies; iSpecies++) { - Ys = V[RHOS_INDEX+iSpecies]/rho; - //eves = nodes->CalcEve(config, V[TVE_INDEX], iSpecies); - //hs = nodes->CalcHs(config, V[T_INDEX], eves, iSpecies); - // Res_Visc[iSpecies] = -rho*Ds[iSpecies]*dYdn[iSpecies] + Ys*SdYdn; - // Res_Visc[nSpecies+nDim] += Res_Visc[iSpecies]*hs; - // Res_Visc[nSpecies+nDim+1] += Res_Visc[iSpecies]*eves; - } } /*--- Get node thermal conductivity ---*/ @@ -475,26 +464,6 @@ void CNEMONSSolver::BC_HeatFluxCatalytic_Wall(CGeometry *geometry, /*--- Viscous contribution to the residual at the wall ---*/ LinSysRes.SubtractBlock(iPoint, Res_Visc); - // /*--- Apply the non-catalytic wall boundary ---*/ - // // Note: We are re-calculating the chemistry residual and adding it to - // // the linear system to eliminate the contribution from the solution - // // (convention is to subtract sources) - // sour_numerics->SetConservative(nodes->GetSolution(), - // nodes->GetSolution() ); - // sour_numerics->SetPrimitive (nodes->GetPrimVar() , - // nodes->GetPrimVar() ); - // sour_numerics->SetdPdU (nodes->GetdPdU() , - // nodes->GetdPdU() ); - // sour_numerics->SetdTdU (nodes->GetdTdU() , - // nodes->GetdTdU() ); - // sour_numerics->SetdTvedU (nodes->GetdTvedU() , - // nodes->GetdTvedU() ); - // sour_numerics->SetVolume (geometry->nodes->GetVolume()); - // sour_numerics->ComputeChemistry(Res_Sour, Jacobian_i, config); - // LinSysRes.AddBlock(iPoint, Res_Sour); - // if (implicit) - // Jacobian.AddBlock2Diag(iPoint, iPoint, Jacobian_i); - /*--- Only change velocity-rows of the Jacobian (includes 1 in the diagonal)/ Note that we need to add a contribution for moving walls to the Jacobian. ---*/ if (implicit) { From a66aed7135d6e9a9d5b30103d414bd1f9c1a9a9b Mon Sep 17 00:00:00 2001 From: rois1995 Date: Mon, 12 Oct 2026 11:33:00 +0200 Subject: [PATCH 2/9] Add SetNumericsThermoState helper to the NEMO solvers The block that passes dPdU, dTdU, dTvedU, Eve, Cvve and Gamma of two points to a numerics object was repeated in the centered, upwind, viscous and boundary loops. Replace the nine copies by one helper in CNEMOEulerSolver. The helper always sets Gamma. BC_Supersonic_Inlet did not set it before, so the AUSM+-up2 and AUSM+M schemes read the value of a previous call there. Co-Authored-By: Claude Opus 5.5 Claude-Session: https://claude.ai/code/session_01L7oxzVir4zENrpy3K3W6Fm --- SU2_CFD/include/solvers/CNEMOEulerSolver.hpp | 18 +++++++ SU2_CFD/src/solvers/CNEMOEulerSolver.cpp | 56 +++----------------- SU2_CFD/src/solvers/CNEMONSSolver.cpp | 6 +-- 3 files changed, 27 insertions(+), 53 deletions(-) diff --git a/SU2_CFD/include/solvers/CNEMOEulerSolver.hpp b/SU2_CFD/include/solvers/CNEMOEulerSolver.hpp index f1955c46a1aa..c60535dc5f48 100644 --- a/SU2_CFD/include/solvers/CNEMOEulerSolver.hpp +++ b/SU2_CFD/include/solvers/CNEMOEulerSolver.hpp @@ -60,6 +60,24 @@ class CNEMOEulerSolver : public CFVMFlowSolverBaseSetdPdU (nodes_i->GetdPdU(iPoint), nodes_j->GetdPdU(jPoint)); + numerics->SetdTdU (nodes_i->GetdTdU(iPoint), nodes_j->GetdTdU(jPoint)); + numerics->SetdTvedU(nodes_i->GetdTvedU(iPoint), nodes_j->GetdTvedU(jPoint)); + numerics->SetEve (nodes_i->GetEve(iPoint), nodes_j->GetEve(jPoint)); + numerics->SetCvve (nodes_i->GetCvve(iPoint), nodes_j->GetCvve(jPoint)); + numerics->SetGamma (nodes_i->GetGamma(iPoint), nodes_j->GetGamma(jPoint)); + } + /*! * \brief Set the maximum value of the eigenvalue. * \param[in] geometry - Geometrical definition of the problem. diff --git a/SU2_CFD/src/solvers/CNEMOEulerSolver.cpp b/SU2_CFD/src/solvers/CNEMOEulerSolver.cpp index b158b525688c..5aaf2ec9b8ff 100644 --- a/SU2_CFD/src/solvers/CNEMOEulerSolver.cpp +++ b/SU2_CFD/src/solvers/CNEMOEulerSolver.cpp @@ -439,12 +439,7 @@ void CNEMOEulerSolver::Centered_Residual(CGeometry *geometry, CSolver **solver_c numerics->SetPrimitive (nodes->GetPrimitive(iPoint), nodes->GetPrimitive(jPoint)); /*--- Pass supplementary information to CNumerics ---*/ - numerics->SetdPdU ( nodes->GetdPdU(iPoint), nodes->GetdPdU(jPoint)); - numerics->SetdTdU ( nodes->GetdTdU(iPoint), nodes->GetdTdU(jPoint)); - numerics->SetdTvedU( nodes->GetdTvedU(iPoint), nodes->GetdTvedU(jPoint)); - numerics->SetEve ( nodes->GetEve(iPoint), nodes->GetEve(jPoint)); - numerics->SetCvve ( nodes->GetCvve(iPoint), nodes->GetCvve(jPoint)); - numerics->SetGamma ( nodes->GetGamma(iPoint), nodes->GetGamma(jPoint)); + SetNumericsThermoState(numerics, nodes, iPoint, nodes, jPoint); /*--- Set the largest convective eigenvalue ---*/ numerics->SetLambda(nodes->GetLambda(iPoint), nodes->GetLambda(jPoint)); @@ -520,12 +515,7 @@ void CNEMOEulerSolver::Upwind_Residual(CGeometry *geometry, CSolver **solver_con numerics->SetPrimitive (V_i, V_j); numerics->SetConservative(U_i, U_j); - numerics->SetdPdU (nodes->GetdPdU(iPoint), nodes->GetdPdU(jPoint)); - numerics->SetdTdU (nodes->GetdTdU(iPoint), nodes->GetdTdU(jPoint)); - numerics->SetdTvedU(nodes->GetdTvedU(iPoint), nodes->GetdTvedU(jPoint)); - numerics->SetEve (nodes->GetEve(iPoint), nodes->GetEve(jPoint)); - numerics->SetCvve (nodes->GetCvve(iPoint), nodes->GetCvve(jPoint)); - numerics->SetGamma (nodes->GetGamma(iPoint), nodes->GetGamma(jPoint)); + SetNumericsThermoState(numerics, nodes, iPoint, nodes, jPoint); } else { @@ -1539,12 +1529,7 @@ void CNEMOEulerSolver::BC_Far_Field(CGeometry *geometry, node_infty->GetPrimitive(0)); /*--- Pass supplementary information to CNumerics ---*/ - conv_numerics->SetdPdU (nodes->GetdPdU(iPoint), node_infty->GetdPdU(0)); - conv_numerics->SetdTdU (nodes->GetdTdU(iPoint), node_infty->GetdTdU(0)); - conv_numerics->SetdTvedU(nodes->GetdTvedU(iPoint), node_infty->GetdTvedU(0)); - conv_numerics->SetEve (nodes->GetEve(iPoint), node_infty->GetEve(0)); - conv_numerics->SetCvve (nodes->GetCvve(iPoint), node_infty->GetCvve(0)); - conv_numerics->SetGamma (nodes->GetGamma(iPoint), node_infty->GetGamma(0)); + SetNumericsThermoState(conv_numerics, nodes, iPoint, node_infty, 0); if(config->GetKind_Upwind_Flow() == UPWIND::AUSMPLUSM) conv_numerics->SetSensor (nodes->GetSensor(iPoint), nodes->GetSensor(iPoint)); @@ -1580,12 +1565,7 @@ void CNEMOEulerSolver::BC_Far_Field(CGeometry *geometry, node_infty->GetGradient_Primitive(0)); /*--- Pass supplementary information to CNumerics ---*/ - visc_numerics->SetdPdU(nodes->GetdPdU(iPoint), node_infty->GetdPdU(0)); - visc_numerics->SetdTdU(nodes->GetdTdU(iPoint), node_infty->GetdTdU(0)); - visc_numerics->SetdTvedU(nodes->GetdTvedU(iPoint), - node_infty->GetdTvedU(0)); - visc_numerics->SetEve(nodes->GetEve(iPoint), node_infty->GetEve(0)); - visc_numerics->SetCvve(nodes->GetCvve(iPoint), node_infty->GetCvve(0)); + SetNumericsThermoState(visc_numerics, nodes, iPoint, node_infty, 0); /*--- Species diffusion coefficients ---*/ visc_numerics->SetDiffusionCoeff(nodes->GetDiffusionCoeff(iPoint), @@ -2049,12 +2029,7 @@ void CNEMOEulerSolver::BC_Outlet(CGeometry *geometry, CSolver **solver_container conv_numerics->SetGridVel(geometry->nodes->GetGridVel(iPoint), geometry->nodes->GetGridVel(iPoint)); /*--- Passing supplementary information to CNumerics ---*/ - conv_numerics->SetdPdU (nodes->GetdPdU(iPoint), node_infty->GetdPdU(0)); - conv_numerics->SetdTdU (nodes->GetdTdU(iPoint), node_infty->GetdTdU(0)); - conv_numerics->SetdTvedU(nodes->GetdTvedU(iPoint), node_infty->GetdTvedU(0)); - conv_numerics->SetEve (nodes->GetEve(iPoint), node_infty->GetEve(0)); - conv_numerics->SetCvve (nodes->GetCvve(iPoint), node_infty->GetCvve(0)); - conv_numerics->SetGamma (nodes->GetGamma(iPoint), node_infty->GetGamma(0)); + SetNumericsThermoState(conv_numerics, nodes, iPoint, node_infty, 0); if(config->GetKind_Upwind_Flow() == UPWIND::AUSMPLUSM) conv_numerics->SetSensor (nodes->GetSensor(iPoint), nodes->GetSensor(iPoint)); @@ -2142,11 +2117,7 @@ void CNEMOEulerSolver::BC_Supersonic_Inlet( node_inlet.GetPrimitive(0)); /*--- Pass supplementary info to CNumerics ---*/ - conv_numerics->SetdPdU(nodes->GetdPdU(iPoint), node_inlet.GetdPdU(0)); - conv_numerics->SetdTdU(nodes->GetdTdU(iPoint), node_inlet.GetdTdU(0)); - conv_numerics->SetdTvedU(nodes->GetdTvedU(iPoint), node_inlet.GetdTvedU(0)); - conv_numerics->SetEve(nodes->GetEve(iPoint), node_inlet.GetEve(0)); - conv_numerics->SetCvve(nodes->GetCvve(iPoint), node_inlet.GetCvve(0)); + SetNumericsThermoState(conv_numerics, nodes, iPoint, &node_inlet, 0); const bool dynamic_grid = config->GetGrid_Movement(); if (dynamic_grid) @@ -2184,12 +2155,7 @@ void CNEMOEulerSolver::BC_Supersonic_Inlet( node_inlet.GetGradient_Primitive(0)); /*--- Pass supplementary information to CNumerics ---*/ - visc_numerics->SetdPdU(nodes->GetdPdU(iPoint), node_inlet.GetdPdU(0)); - visc_numerics->SetdTdU(nodes->GetdTdU(iPoint), node_inlet.GetdTdU(0)); - visc_numerics->SetdTvedU(nodes->GetdTvedU(iPoint), - node_inlet.GetdTvedU(0)); - visc_numerics->SetEve(nodes->GetEve(iPoint), node_inlet.GetEve(0)); - visc_numerics->SetCvve(nodes->GetCvve(iPoint), node_inlet.GetCvve(0)); + SetNumericsThermoState(visc_numerics, nodes, iPoint, &node_inlet, 0); /*--- Species diffusion coefficients ---*/ visc_numerics->SetDiffusionCoeff(nodes->GetDiffusionCoeff(iPoint), @@ -2256,13 +2222,7 @@ void CNEMOEulerSolver::BC_Supersonic_Outlet( nodes->GetSolution(iPoint)); /*--- Pass supplementary information to CNumerics ---*/ - conv_numerics->SetdPdU(nodes->GetdPdU(iPoint), nodes->GetdPdU(iPoint)); - conv_numerics->SetdTdU(nodes->GetdTdU(iPoint), nodes->GetdTdU(iPoint)); - conv_numerics->SetdTvedU(nodes->GetdTvedU(iPoint), - nodes->GetdTvedU(iPoint)); - conv_numerics->SetEve(nodes->GetEve(iPoint), nodes->GetEve(iPoint)); - conv_numerics->SetCvve(nodes->GetCvve(iPoint), nodes->GetCvve(iPoint)); - conv_numerics->SetGamma(nodes->GetGamma(iPoint), nodes->GetGamma(iPoint)); + SetNumericsThermoState(conv_numerics, nodes, iPoint, nodes, iPoint); const bool dynamic_grid = config->GetGrid_Movement(); if (dynamic_grid) diff --git a/SU2_CFD/src/solvers/CNEMONSSolver.cpp b/SU2_CFD/src/solvers/CNEMONSSolver.cpp index 513344fcfc31..a3e364d4f342 100644 --- a/SU2_CFD/src/solvers/CNEMONSSolver.cpp +++ b/SU2_CFD/src/solvers/CNEMONSSolver.cpp @@ -199,11 +199,7 @@ void CNEMONSSolver::Viscous_Residual(CGeometry *geometry, nodes->GetGradient_Primitive(jPoint) ); /*--- Pass supplementary information to CNumerics ---*/ - numerics->SetdPdU (nodes->GetdPdU(iPoint), nodes->GetdPdU(jPoint)); - numerics->SetdTdU (nodes->GetdTdU(iPoint), nodes->GetdTdU(jPoint)); - numerics->SetdTvedU(nodes->GetdTvedU(iPoint), nodes->GetdTvedU(jPoint)); - numerics->SetEve (nodes->GetEve(iPoint), nodes->GetEve(jPoint)); - numerics->SetCvve (nodes->GetCvve(iPoint), nodes->GetCvve(jPoint)); + SetNumericsThermoState(numerics, nodes, iPoint, nodes, jPoint); /*--- Species diffusion coefficients ---*/ numerics->SetDiffusionCoeff(nodes->GetDiffusionCoeff(iPoint), From 414724c0ede4de787a7ac80d318e2605cc95f200 Mon Sep 17 00:00:00 2001 From: rois1995 Date: Mon, 12 Oct 2026 11:46:53 +0200 Subject: [PATCH 3/9] Share the wall velocity and wall energy code of the NEMO walls Add SetWallVelocity_Strong and SetWallEnergyFlux_Weak to CNEMONSSolver and use them in the heat-flux, isothermal and Smoluchowski-Maxwell walls, which repeated the same statements. Reduce BC_HeatFluxCatalytic_Wall to its error message. The code after the error could not run and had no catalytic contribution. Remove the unused GradY array of the catalytic isothermal wall. No functional change. Co-Authored-By: Claude Opus 5.5 Claude-Session: https://claude.ai/code/session_01L7oxzVir4zENrpy3K3W6Fm --- SU2_CFD/include/solvers/CNEMONSSolver.hpp | 23 +++ SU2_CFD/src/solvers/CNEMONSSolver.cpp | 216 +++++----------------- 2 files changed, 69 insertions(+), 170 deletions(-) diff --git a/SU2_CFD/include/solvers/CNEMONSSolver.hpp b/SU2_CFD/include/solvers/CNEMONSSolver.hpp index 0e4fe68a5b2f..11bccf88e8c0 100644 --- a/SU2_CFD/include/solvers/CNEMONSSolver.hpp +++ b/SU2_CFD/include/solvers/CNEMONSSolver.hpp @@ -72,6 +72,29 @@ class CNEMONSSolver final : public CNEMOEulerSolver { */ void SetTau_Wall_WF(CGeometry *geometry, CSolver** solver_container, const CConfig* config); + /*! + * \brief Impose the velocity of a wall point as a strong (Dirichlet) boundary condition. + * The velocity is fixed and the momentum residual and truncation error are set to zero. + * \param[in] iPoint - Index of the wall point. + * \param[in] val_velocity - Velocity to impose (zero for no-slip, slip velocity otherwise). + */ + void SetWallVelocity_Strong(unsigned long iPoint, const su2double *val_velocity); + + /*! + * \brief Compute the residual of the two energy equations that weakly imposes the wall temperatures. + * \param[in] iPoint - Index of the wall point. + * \param[in] jPoint - Index of the closest normal neighbor of the wall point. + * \param[in] val_Twall - Target translational-rotational temperature at the wall. + * \param[in] val_Tvewall - Target vibrational-electronic temperature at the wall. + * \param[in] val_pcontrol - "Proportional control" constant of the temperature difference term. + * \param[in] val_area - Area of the dual-grid face on the wall. + * \param[in] val_dist - Distance between the wall point and its normal neighbor. + * \param[out] val_residual - Residual, only the two energy entries are set. + */ + void SetWallEnergyFlux_Weak(unsigned long iPoint, unsigned long jPoint, su2double val_Twall, su2double val_Tvewall, + su2double val_pcontrol, su2double val_area, su2double val_dist, + su2double *val_residual) const; + public: /*! diff --git a/SU2_CFD/src/solvers/CNEMONSSolver.cpp b/SU2_CFD/src/solvers/CNEMONSSolver.cpp index a3e364d4f342..aa419fc9dd4b 100644 --- a/SU2_CFD/src/solvers/CNEMONSSolver.cpp +++ b/SU2_CFD/src/solvers/CNEMONSSolver.cpp @@ -238,6 +238,39 @@ void CNEMONSSolver::Viscous_Residual(CGeometry *geometry, } //iEdge } +void CNEMONSSolver::SetWallVelocity_Strong(unsigned long iPoint, const su2double *val_velocity) { + + /*--- Store the wall velocity in the old solution, it is not changed by the update. ---*/ + nodes->SetVelocity_Old(iPoint, val_velocity); + + /*--- Remove any contribution to the momentum residual and truncation error at this node. ---*/ + for (auto iDim = 0ul; iDim < nDim; iDim++) + LinSysRes(iPoint, nSpecies+iDim) = 0.0; + nodes->SetVel_ResTruncError_Zero(iPoint); +} + +void CNEMONSSolver::SetWallEnergyFlux_Weak(unsigned long iPoint, unsigned long jPoint, su2double val_Twall, + su2double val_Tvewall, su2double val_pcontrol, su2double val_area, + su2double val_dist, su2double *val_residual) const { + + /*--- Temperatures at the wall point (i) and at its closest normal neighbor (j) ---*/ + const su2double Ti = nodes->GetTemperature(iPoint); + const su2double Tj = nodes->GetTemperature(jPoint); + const su2double Tvei = nodes->GetTemperature_ve(iPoint); + const su2double Tvej = nodes->GetTemperature_ve(jPoint); + + /*--- Thermal conductivities at the wall point ---*/ + const su2double ktr = nodes->GetThermalConductivity(iPoint); + const su2double kve = nodes->GetThermalConductivity_ve(iPoint); + + /*--- Conductive heat flux from the normal neighbor, plus a term that is proportional + to the difference between the target wall temperature and the current one. The + second term drives the solution at the wall towards the target temperature. ---*/ + val_residual[nSpecies+nDim] = ((ktr*(Ti-Tj) + kve*(Tvei-Tvej)) + + (ktr*(val_Twall-Ti) + kve*(val_Tvewall-Tvei))*val_pcontrol)*val_area/val_dist; + val_residual[nSpecies+nDim+1] = (kve*(Tvei-Tvej) + kve*(val_Tvewall-Tvei)*val_pcontrol)*val_area/val_dist; +} + void CNEMONSSolver::BC_HeatFluxNonCatalytic_Wall(CGeometry *geometry, CSolver **solver_container, CNumerics *conv_numerics, @@ -300,12 +333,8 @@ void CNEMONSSolver::BC_HeatFluxNonCatalytic_Wall(CGeometry *geometry, /*--- Impose the value of the velocity as a strong boundary condition (Dirichlet). Fix the velocity and remove any contribution to the residual at this node. ---*/ - su2double zero[MAXNDIM] = {0.0}; - nodes->SetVelocity_Old(iPoint, zero); - - for (auto iDim = 0ul; iDim < nDim; iDim++) - LinSysRes(iPoint, nSpecies+iDim) = 0.0; - nodes->SetVel_ResTruncError_Zero(iPoint); + const su2double zero[MAXNDIM] = {0.0}; + SetWallVelocity_Strong(iPoint, zero); /*--- Apply viscous residual to the linear system ---*/ LinSysRes.SubtractBlock(iPoint, Res_Visc); @@ -347,135 +376,8 @@ void CNEMONSSolver::BC_HeatFluxCatalytic_Wall(CGeometry *geometry, unsigned short val_marker) { SU2_ZONE_SCOPED + /*--- The catalytic contribution to a heat-flux wall is not implemented. ---*/ SU2_MPI::Error("BC_HEATFLUX with catalytic wall: Not operational in NEMO.", CURRENT_FUNCTION); - //TODO: SCALE WITH EDDY VISC - /*--- Local variables ---*/ - unsigned long iPoint; - su2double rho, Ys; - su2double *Normal, Area; - su2double *Ds, *dYdn, SdYdn; - su2double **GradY; - - /*--- Assign booleans ---*/ - bool implicit = (config->GetKind_TimeIntScheme() == EULER_IMPLICIT); - bool catalytic = false; - - /*--- Set "Proportional control" coefficient ---*/ - su2double pcontrol = 0.6; - - /*--- Identify the boundary by string name ---*/ - string Marker_Tag = config->GetMarker_All_TagBound(val_marker); - - /*--- Get the specified wall heat flux from config ---*/ - su2double Wall_HeatFlux = config->GetWall_HeatFlux(Marker_Tag) / config->GetHeat_Flux_Ref(); - if (config->GetIntegrated_HeatFlux()) { - Wall_HeatFlux /= geometry->GetSurfaceArea(config, val_marker); - } - - /*--- Get the locations of the primitive variables ---*/ - const unsigned short T_INDEX = nodes->GetTIndex(); - const unsigned short TVE_INDEX = nodes->GetTveIndex(); - const unsigned short RHOS_INDEX = nodes->GetRhosIndex(); - const unsigned short RHO_INDEX = nodes->GetRhoIndex(); - - /*--- Allocate arrays ---*/ - dYdn = new su2double[nSpecies]; - GradY = new su2double*[nSpecies]; - for (auto iSpecies = 0ul; iSpecies < nSpecies; iSpecies++) - GradY[iSpecies] = new su2double[nDim]; - - /*--- Loop over all of the vertices on this boundary marker ---*/ - for (auto iVertex = 0ul; iVertex < geometry->nVertex[val_marker]; iVertex++) { - iPoint = geometry->vertex[val_marker][iVertex]->GetNode(); - - /*--- Check if the node belongs to the domain (i.e, not a halo node) ---*/ - if (geometry->nodes->GetDomain(iPoint)) { - - /*--- Compute dual-grid area and boundary normal ---*/ - Normal = geometry->vertex[val_marker][iVertex]->GetNormal(); - Area = GeometryToolbox::Norm(nDim, Normal); - - /*--- Initialize the convective & viscous residuals to zero ---*/ - for (auto iVar = 0ul; iVar < nVar; iVar++) { - Res_Visc[iVar] = 0.0; - Res_Sour[iVar] = 0.0; - } - - /*--- Assign wall velocity to "Vector" array ---*/ - for (auto iDim = 0ul; iDim < nDim; iDim++) Vector[iDim] = 0.0; - - /*--- Set the residual, truncation error, and velocity value ---*/ - nodes->SetVelocity_Old(iPoint,Vector); - for (auto iDim = 0ul; iDim < nDim; iDim++) { - LinSysRes(iPoint, nSpecies+iDim) = 0.0; - nodes->SetVal_ResTruncError_Zero(iPoint,nSpecies+iDim); - } - - /*--- Get temperature gradient information ---*/ - const auto V = nodes->GetPrimitive(iPoint); - const auto GradV = nodes->GetGradient_Primitive(iPoint); - su2double dTdn = 0.0; - su2double dTvedn = 0.0; - for (auto iDim = 0ul; iDim < nDim; iDim++) { - dTdn += GradV[T_INDEX][iDim]*Normal[iDim]; - dTvedn += GradV[TVE_INDEX][iDim]*Normal[iDim]; - } - - if (catalytic) { - SU2_MPI::Error("NEED TO IMPLEMENT CATALYTIC BOUNDARIES IN HEATFLUX!!!",CURRENT_FUNCTION); - } - else { - - /*--- Rename for convenience ---*/ - rho = V[RHO_INDEX]; - Ds = nodes->GetDiffusionCoeff(iPoint); - - /*--- Calculate normal derivative of mass fraction ---*/ - for (auto iSpecies = 0ul; iSpecies < nSpecies; iSpecies++) { - Ys = V[RHOS_INDEX+iSpecies]/rho; - dYdn[iSpecies] = 0.0; - for (auto iDim = 0ul; iDim < nDim; iDim++) - dYdn[iSpecies] += 1.0/rho * (GradV[RHOS_INDEX+iSpecies][iDim] - - Ys*GradV[RHO_INDEX][iDim])*Normal[iDim]; - } - - /*--- Calculate supplementary quantities ---*/ - SdYdn = 0.0; - for (auto iSpecies = 0ul; iSpecies < nSpecies; iSpecies++) - SdYdn += rho*Ds[iSpecies]*dYdn[iSpecies]; - } - - /*--- Get node thermal conductivity ---*/ - su2double ktr = nodes->GetThermalConductivity(iPoint); - su2double kve = nodes->GetThermalConductivity_ve(iPoint); - - /*--- Set the residual on the boundary with the specified heat flux ---*/ - // Note: Contributions from qtr and qve are used for proportional control - // to drive the solution toward the specified heatflux more quickly. - Res_Visc[nSpecies+nDim] += pcontrol*(ktr*dTdn+kve*dTvedn) + - Wall_HeatFlux*Area; - Res_Visc[nSpecies+nDim+1] += pcontrol*(kve*dTvedn) + - Wall_HeatFlux*Area; - - /*--- Viscous contribution to the residual at the wall ---*/ - LinSysRes.SubtractBlock(iPoint, Res_Visc); - - /*--- Only change velocity-rows of the Jacobian (includes 1 in the diagonal)/ - Note that we need to add a contribution for moving walls to the Jacobian. ---*/ - if (implicit) { - /*--- Enforce the no-slip boundary condition in a strong way ---*/ - for (auto iVar = nSpecies; iVar < nSpecies+nDim; iVar++) { - Jacobian.DeleteValsRowi(iPoint, iVar); - } - } - - } - } - - for (auto iSpecies = 0ul; iSpecies < nSpecies; iSpecies++) - delete [] GradY[iSpecies]; - delete [] GradY; - delete [] dYdn; } void CNEMONSSolver::BC_Isothermal_Wall(CGeometry *geometry, CSolver **solver_container, @@ -550,30 +452,18 @@ void CNEMONSSolver::BC_IsothermalNonCatalytic_Wall(CGeometry *geometry, /*--- Store the corrected velocity at the wall which will be zero (v = 0), unless there is grid motion (v = u_wall)---*/ - su2double zero[MAXNDIM] = {0.0}; - nodes->SetVelocity_Old(iPoint, zero); + const su2double zero[MAXNDIM] = {0.0}; + SetWallVelocity_Strong(iPoint, zero); /*--- Initialize viscous residual to zero ---*/ for (auto iVar = 0ul; iVar < nVar; iVar ++) {Res_Visc[iVar] = 0.0;} - for (auto iDim = 0ul; iDim < nDim; iDim++) - LinSysRes(iPoint, nSpecies+iDim) = 0.0; - nodes->SetVel_ResTruncError_Zero(iPoint); - - /*--- Calculate the gradient of temperature ---*/ - const su2double Ti = nodes->GetTemperature(iPoint); - const su2double Tj = nodes->GetTemperature(Point_Normal); - const su2double Tvei = nodes->GetTemperature_ve(iPoint); - const su2double Tvej = nodes->GetTemperature_ve(Point_Normal); - /*--- Rename variables for convenience ---*/ const su2double ktr = nodes->GetThermalConductivity(iPoint); const su2double kve = nodes->GetThermalConductivity_ve(iPoint); - /*--- Apply to the linear system ---*/ - Res_Visc[nSpecies+nDim] = ((ktr*(Ti-Tj) + kve*(Tvei-Tvej)) + - (ktr*(Twall-Ti) + kve*(Twall-Tvei))*C)*Area/dist_ij; - Res_Visc[nSpecies+nDim+1] = (kve*(Tvei-Tvej) + kve*(Twall-Tvei) *C)*Area/dist_ij; + /*--- Weak enforcement of the wall temperature on both energy equations ---*/ + SetWallEnergyFlux_Weak(iPoint, Point_Normal, Twall, Twall, C, Area, dist_ij, Res_Visc); /*--- Calculate Jacobian for implicit time stepping ---*/ if (implicit) { @@ -631,7 +521,7 @@ void CNEMONSSolver::BC_IsothermalCatalytic_Wall(CGeometry *geometry, sour_numerics, config, val_marker); /*--- Local variables ---*/ - su2double **GradY, **dVdU; + su2double **dVdU; /*--- Assign booleans ---*/ const bool implicit = (config->GetKind_TimeIntScheme() == EULER_IMPLICIT); @@ -649,9 +539,6 @@ void CNEMONSSolver::BC_IsothermalCatalytic_Wall(CGeometry *geometry, const unsigned short TVE_INDEX = nodes->GetTveIndex(); /*--- Allocate arrays ---*/ - GradY = new su2double*[nSpecies]; - for (auto iSpecies = 0ul; iSpecies < nSpecies; iSpecies++) - GradY[iSpecies] = new su2double[nDim]; dVdU = new su2double*[nVar]; for (auto iVar = 0ul; iVar < nVar; iVar++) dVdU[iVar] = new su2double[nVar]; @@ -806,9 +693,6 @@ void CNEMONSSolver::BC_IsothermalCatalytic_Wall(CGeometry *geometry, } } - for (auto iSpecies = 0ul; iSpecies < nSpecies; iSpecies++) - delete [] GradY[iSpecies]; - delete [] GradY; for (auto iVar = 0ul; iVar < nVar; iVar++) delete [] dVdU[iVar]; delete [] dVdU; @@ -877,11 +761,8 @@ void CNEMONSSolver::BC_Smoluchowski_Maxwell(CGeometry *geometry, /*--- Calculate Pressure ---*/ const su2double Pi = nodes->GetPressure(iPoint); - /*--- Calculate the gradient of temperature ---*/ + /*--- Temperature at the wall point ---*/ const su2double Ti = nodes->GetTemperature(iPoint); - const su2double Tj = nodes->GetTemperature(jPoint); - const su2double Tvei = nodes->GetTemperature_ve(iPoint); - const su2double Tvej = nodes->GetTemperature_ve(jPoint); /*--- Rename variables for convenience ---*/ const su2double ktr = nodes->GetThermalConductivity(iPoint); @@ -949,18 +830,13 @@ void CNEMONSSolver::BC_Smoluchowski_Maxwell(CGeometry *geometry, Vector[iDim] = (1.0-alpha_V)*nodes->GetVelocity(iPoint,iDim)+(alpha_V)*Vector[iDim]; } - nodes->SetVelocity_Old(iPoint,Vector); + /*--- Impose the slip velocity as a strong boundary condition ---*/ + SetWallVelocity_Strong(iPoint, Vector); - for (auto iDim = 0ul; iDim < nDim; iDim++) { - LinSysRes(iPoint, nSpecies+iDim) = 0.0; - nodes->SetVal_ResTruncError_Zero(iPoint,nSpecies+iDim); - } + /*--- Weak enforcement of the slip temperatures on both energy equations ---*/ + SetWallEnergyFlux_Weak(iPoint, jPoint, Tslip, Tslip_ve, C, Area, dij, Res_Visc); /*--- Apply to the linear system ---*/ - Res_Visc[nSpecies+nDim] = ((ktr*(Ti-Tj) + kve*(Tvei-Tvej)) + - (ktr*(Tslip-Ti) + kve*(Tslip_ve-Tvei))*C)*Area/dij; - Res_Visc[nSpecies+nDim+1] = (kve*(Tvei-Tvej) + kve*(Tslip_ve-Tvei)*C)*Area/dij; - LinSysRes.SubtractBlock(iPoint, Res_Visc); } END_SU2_OMP_FOR From 6ed4f1527ba9c003ad30424db0daa46185ced99c Mon Sep 17 00:00:00 2001 From: rois1995 Date: Mon, 12 Oct 2026 11:48:00 +0200 Subject: [PATCH 4/9] Do not set Gamma in SetNumericsThermoState Keep the behavior of the code before the helper was added: only the call sites that passed Gamma to the numerics before still do it, with their own SetGamma call. The helper passes dPdU, dTdU, dTvedU, Eve and Cvve. Co-Authored-By: Claude Opus 5.5 Claude-Session: https://claude.ai/code/session_01L7oxzVir4zENrpy3K3W6Fm --- SU2_CFD/include/solvers/CNEMOEulerSolver.hpp | 2 +- SU2_CFD/src/solvers/CNEMOEulerSolver.cpp | 5 +++++ 2 files changed, 6 insertions(+), 1 deletion(-) diff --git a/SU2_CFD/include/solvers/CNEMOEulerSolver.hpp b/SU2_CFD/include/solvers/CNEMOEulerSolver.hpp index c60535dc5f48..a454602b5205 100644 --- a/SU2_CFD/include/solvers/CNEMOEulerSolver.hpp +++ b/SU2_CFD/include/solvers/CNEMOEulerSolver.hpp @@ -62,6 +62,7 @@ class CNEMOEulerSolver : public CFVMFlowSolverBaseSetdTvedU(nodes_i->GetdTvedU(iPoint), nodes_j->GetdTvedU(jPoint)); numerics->SetEve (nodes_i->GetEve(iPoint), nodes_j->GetEve(jPoint)); numerics->SetCvve (nodes_i->GetCvve(iPoint), nodes_j->GetCvve(jPoint)); - numerics->SetGamma (nodes_i->GetGamma(iPoint), nodes_j->GetGamma(jPoint)); } /*! diff --git a/SU2_CFD/src/solvers/CNEMOEulerSolver.cpp b/SU2_CFD/src/solvers/CNEMOEulerSolver.cpp index 5aaf2ec9b8ff..0c03a2b0bbd4 100644 --- a/SU2_CFD/src/solvers/CNEMOEulerSolver.cpp +++ b/SU2_CFD/src/solvers/CNEMOEulerSolver.cpp @@ -440,6 +440,7 @@ void CNEMOEulerSolver::Centered_Residual(CGeometry *geometry, CSolver **solver_c /*--- Pass supplementary information to CNumerics ---*/ SetNumericsThermoState(numerics, nodes, iPoint, nodes, jPoint); + numerics->SetGamma(nodes->GetGamma(iPoint), nodes->GetGamma(jPoint)); /*--- Set the largest convective eigenvalue ---*/ numerics->SetLambda(nodes->GetLambda(iPoint), nodes->GetLambda(jPoint)); @@ -516,6 +517,7 @@ void CNEMOEulerSolver::Upwind_Residual(CGeometry *geometry, CSolver **solver_con numerics->SetPrimitive (V_i, V_j); numerics->SetConservative(U_i, U_j); SetNumericsThermoState(numerics, nodes, iPoint, nodes, jPoint); + numerics->SetGamma(nodes->GetGamma(iPoint), nodes->GetGamma(jPoint)); } else { @@ -1530,6 +1532,7 @@ void CNEMOEulerSolver::BC_Far_Field(CGeometry *geometry, /*--- Pass supplementary information to CNumerics ---*/ SetNumericsThermoState(conv_numerics, nodes, iPoint, node_infty, 0); + conv_numerics->SetGamma(nodes->GetGamma(iPoint), node_infty->GetGamma(0)); if(config->GetKind_Upwind_Flow() == UPWIND::AUSMPLUSM) conv_numerics->SetSensor (nodes->GetSensor(iPoint), nodes->GetSensor(iPoint)); @@ -2030,6 +2033,7 @@ void CNEMOEulerSolver::BC_Outlet(CGeometry *geometry, CSolver **solver_container /*--- Passing supplementary information to CNumerics ---*/ SetNumericsThermoState(conv_numerics, nodes, iPoint, node_infty, 0); + conv_numerics->SetGamma(nodes->GetGamma(iPoint), node_infty->GetGamma(0)); if(config->GetKind_Upwind_Flow() == UPWIND::AUSMPLUSM) conv_numerics->SetSensor (nodes->GetSensor(iPoint), nodes->GetSensor(iPoint)); @@ -2223,6 +2227,7 @@ void CNEMOEulerSolver::BC_Supersonic_Outlet( /*--- Pass supplementary information to CNumerics ---*/ SetNumericsThermoState(conv_numerics, nodes, iPoint, nodes, iPoint); + conv_numerics->SetGamma(nodes->GetGamma(iPoint), nodes->GetGamma(iPoint)); const bool dynamic_grid = config->GetGrid_Movement(); if (dynamic_grid) From 4ec981cf6ffc237e86fd0a383960ccddb1422c71 Mon Sep 17 00:00:00 2001 From: rois1995 Date: Mon, 12 Oct 2026 11:57:06 +0200 Subject: [PATCH 5/9] Share the viscous boundary flux code of the NEMO solvers Add SetNumericsTransportState and AddViscousBoundaryResidual to CNEMOEulerSolver. BC_Far_Field and BC_Supersonic_Inlet had the same viscous block, and Viscous_Residual repeated the transport property calls. Reduce BC_Inlet to its error message. The code after the error could not run. No functional change. Co-Authored-By: Claude Opus 5.5 Claude-Session: https://claude.ai/code/session_01L7oxzVir4zENrpy3K3W6Fm --- SU2_CFD/include/solvers/CNEMOEulerSolver.hpp | 24 ++ SU2_CFD/src/solvers/CNEMOEulerSolver.cpp | 412 +++---------------- SU2_CFD/src/solvers/CNEMONSSolver.cpp | 21 +- 3 files changed, 92 insertions(+), 365 deletions(-) diff --git a/SU2_CFD/include/solvers/CNEMOEulerSolver.hpp b/SU2_CFD/include/solvers/CNEMOEulerSolver.hpp index a454602b5205..1c4e0a97eb67 100644 --- a/SU2_CFD/include/solvers/CNEMOEulerSolver.hpp +++ b/SU2_CFD/include/solvers/CNEMOEulerSolver.hpp @@ -78,6 +78,30 @@ class CNEMOEulerSolver : public CFVMFlowSolverBaseSetCvve (nodes_i->GetCvve(iPoint), nodes_j->GetCvve(jPoint)); } + /*! + * \brief Pass the transport properties of two points of the solver to a viscous numerics object. + * \param[in,out] numerics - Numerics object that receives the data. + * \param[in] iPoint - Index of point i. + * \param[in] jPoint - Index of point j. + */ + void SetNumericsTransportState(CNumerics *numerics, unsigned long iPoint, unsigned long jPoint); + + /*! + * \brief Compute the viscous flux between a boundary point and a prescribed boundary state, + * and add it to the residual and to the Jacobian. + * \param[in] geometry - Geometrical definition of the problem. + * \param[in] visc_numerics - Description of the numerical method for viscous terms. + * \param[in] config - Definition of the particular problem. + * \param[in] iPoint - Index of the boundary point. + * \param[in] Point_Normal - Index of the closest interior neighbor of the boundary point. + * \param[in] Normal - Normal vector of the boundary face, the same that is passed to the convective numerics. + * \param[in] nodes_bound - Variable container that holds the boundary state. + * \param[in] iPoint_bound - Index of the boundary state in its container. + */ + void AddViscousBoundaryResidual(CGeometry *geometry, CNumerics *visc_numerics, const CConfig *config, + unsigned long iPoint, unsigned long Point_Normal, const su2double *Normal, + CNEMOEulerVariable *nodes_bound, unsigned long iPoint_bound); + /*! * \brief Set the maximum value of the eigenvalue. * \param[in] geometry - Geometrical definition of the problem. diff --git a/SU2_CFD/src/solvers/CNEMOEulerSolver.cpp b/SU2_CFD/src/solvers/CNEMOEulerSolver.cpp index 0c03a2b0bbd4..27fc8e43d251 100644 --- a/SU2_CFD/src/solvers/CNEMOEulerSolver.cpp +++ b/SU2_CFD/src/solvers/CNEMOEulerSolver.cpp @@ -1397,6 +1397,67 @@ void CNEMOEulerSolver::SetReferenceValues(const CConfig& config) { } +void CNEMOEulerSolver::SetNumericsTransportState(CNumerics *numerics, unsigned long iPoint, unsigned long jPoint) { + + /*--- Species diffusion coefficients ---*/ + numerics->SetDiffusionCoeff(nodes->GetDiffusionCoeff(iPoint), + nodes->GetDiffusionCoeff(jPoint)); + + /*--- Laminar viscosity ---*/ + numerics->SetLaminarViscosity(nodes->GetLaminarViscosity(iPoint), + nodes->GetLaminarViscosity(jPoint)); + + /*--- Eddy viscosity ---*/ + numerics->SetEddyViscosity(nodes->GetEddyViscosity(iPoint), + nodes->GetEddyViscosity(jPoint)); + + /*--- Thermal conductivity ---*/ + numerics->SetThermalConductivity(nodes->GetThermalConductivity(iPoint), + nodes->GetThermalConductivity(jPoint)); + + /*--- Vib-el. thermal conductivity ---*/ + numerics->SetThermalConductivity_ve(nodes->GetThermalConductivity_ve(iPoint), + nodes->GetThermalConductivity_ve(jPoint)); +} + +void CNEMOEulerSolver::AddViscousBoundaryResidual(CGeometry *geometry, CNumerics *visc_numerics, const CConfig *config, + unsigned long iPoint, unsigned long Point_Normal, + const su2double *Normal, CNEMOEulerVariable *nodes_bound, + unsigned long iPoint_bound) { + + const bool implicit = (config->GetKind_TimeIntScheme() == EULER_IMPLICIT); + + /*--- Set the normal vector and the coordinates. The second point is the + reflection of the closest interior node with respect to the boundary point. ---*/ + su2double Coord_Reflected[MAXNDIM]; + GeometryToolbox::PointPointReflect(nDim, geometry->nodes->GetCoord(Point_Normal), + geometry->nodes->GetCoord(iPoint), Coord_Reflected); + visc_numerics->SetCoord(geometry->nodes->GetCoord(iPoint), Coord_Reflected); + visc_numerics->SetNormal(Normal); + + /*--- Conservative and primitive variables, and gradient ---*/ + visc_numerics->SetConservative(nodes->GetSolution(iPoint), + nodes_bound->GetSolution(iPoint_bound)); + visc_numerics->SetPrimitive(nodes->GetPrimitive(iPoint), + nodes_bound->GetPrimitive(iPoint_bound)); + visc_numerics->SetPrimVarGradient(nodes->GetGradient_Primitive(iPoint), + nodes_bound->GetGradient_Primitive(iPoint_bound)); + + /*--- Pass supplementary information to CNumerics ---*/ + SetNumericsThermoState(visc_numerics, nodes, iPoint, nodes_bound, iPoint_bound); + + /*--- Transport properties, the values of the boundary point are used for both states ---*/ + SetNumericsTransportState(visc_numerics, iPoint, iPoint); + + /*--- Compute and update residual ---*/ + auto residual = visc_numerics->ComputeResidual(config); + LinSysRes.SubtractBlock(iPoint, residual); + + /*--- Jacobian contribution for implicit integration ---*/ + if (implicit) + Jacobian.SubtractBlock2Diag(iPoint, residual.jacobian_i); +} + void CNEMOEulerSolver::BC_Sym_Plane(CGeometry *geometry, CSolver **solver_container, CNumerics *conv_numerics, CNumerics *visc_numerics, CConfig *config, unsigned short val_marker) { SU2_ZONE_SCOPED @@ -1549,57 +1610,8 @@ void CNEMOEulerSolver::BC_Far_Field(CGeometry *geometry, Jacobian.AddBlock2Diag(iPoint, residual.jacobian_i); /*--- Viscous contribution ---*/ - const bool viscous = config->GetViscous(); - if (viscous) { - su2double Coord_Reflected[MAXNDIM]; - GeometryToolbox::PointPointReflect( - nDim, geometry->nodes->GetCoord(Point_Normal), - geometry->nodes->GetCoord(iPoint), Coord_Reflected); - visc_numerics->SetCoord(geometry->nodes->GetCoord(iPoint), - Coord_Reflected); - visc_numerics->SetNormal(Normal); - - /*--- Primitive variables, and gradient ---*/ - visc_numerics->SetConservative(nodes->GetSolution(iPoint), - node_infty->GetSolution(0)); - visc_numerics->SetPrimitive(nodes->GetPrimitive(iPoint), - node_infty->GetPrimitive(0)); - visc_numerics->SetPrimVarGradient(nodes->GetGradient_Primitive(iPoint), - node_infty->GetGradient_Primitive(0)); - - /*--- Pass supplementary information to CNumerics ---*/ - SetNumericsThermoState(visc_numerics, nodes, iPoint, node_infty, 0); - - /*--- Species diffusion coefficients ---*/ - visc_numerics->SetDiffusionCoeff(nodes->GetDiffusionCoeff(iPoint), - nodes->GetDiffusionCoeff(iPoint)); - - /*--- Laminar viscosity ---*/ - visc_numerics->SetLaminarViscosity(nodes->GetLaminarViscosity(iPoint), - nodes->GetLaminarViscosity(iPoint)); - - /*--- Eddy viscosity ---*/ - visc_numerics->SetEddyViscosity(nodes->GetEddyViscosity(iPoint), - nodes->GetEddyViscosity(iPoint)); - - /*--- Thermal conductivity ---*/ - visc_numerics->SetThermalConductivity( - nodes->GetThermalConductivity(iPoint), - nodes->GetThermalConductivity(iPoint)); - - /*--- Vib-el. thermal conductivity ---*/ - visc_numerics->SetThermalConductivity_ve( - nodes->GetThermalConductivity_ve(iPoint), - nodes->GetThermalConductivity_ve(iPoint)); - - /*--- Compute and update residual ---*/ - auto residual = visc_numerics->ComputeResidual(config); - - LinSysRes.SubtractBlock(iPoint, residual); - if (implicit) { - Jacobian.SubtractBlock2Diag(iPoint, residual.jacobian_i); - } - } + if (config->GetViscous()) + AddViscousBoundaryResidual(geometry, visc_numerics, config, iPoint, Point_Normal, Normal, node_infty, 0); } } } @@ -1607,253 +1619,8 @@ void CNEMOEulerSolver::BC_Far_Field(CGeometry *geometry, void CNEMOEulerSolver::BC_Inlet(CGeometry *geometry, CSolver **solver_container, CNumerics *conv_numerics, CNumerics *visc_numerics, CConfig *config, unsigned short val_marker) { SU2_ZONE_SCOPED + /*--- The subsonic inlet is not implemented for a two-temperature mixture. ---*/ SU2_MPI::Error("BC_INLET: Not operational in NEMO.", CURRENT_FUNCTION); - - unsigned short RHO_INDEX, nSpecies; - - su2double T_Total, P_Total, Velocity[3], Velocity2, H_Total, Temperature, Riemann, - Pressure, Density, Energy, Mach2, SoundSpeed2, SoundSpeed_Total2, Vel_Mag, - alpha, aa, bb, cc, dd, Area, UnitNormal[3] = {0.0}; - - su2double Flow_Dir[MAXNDIM] = {}; - bool implicit = (config->GetKind_TimeIntScheme() == EULER_IMPLICIT); - - bool dynamic_grid = config->GetGrid_Movement(); - su2double Two_Gamma_M1 = 2.0/Gamma_Minus_One; - su2double Gas_Constant = config->GetGas_ConstantND(); - INLET_TYPE Kind_Inlet = config->GetKind_Inlet(); - - auto *U_domain = new su2double[nVar]; auto *U_inlet = new su2double[nVar]; - auto *V_domain = new su2double[nPrimVar]; auto *V_inlet = new su2double[nPrimVar]; - auto *Normal = new su2double[nDim]; - - nSpecies = config->GetnSpecies(); - auto *Spec_Density = new su2double[nSpecies]; - for (auto iSpecies = 0ul; iSpeciesGetRhoIndex(); - - /*--- Loop over all the vertices on this boundary marker ---*/ - for (auto iVertex = 0ul; iVertex < geometry->nVertex[val_marker]; iVertex++) { - const auto iPoint = geometry->vertex[val_marker][iVertex]->GetNode(); - - /*--- Check if the node belongs to the domain (i.e., not a halo node) ---*/ - if (geometry->nodes->GetDomain(iPoint)) { - - /*--- Normal vector for this vertex (negate for outward convention) ---*/ - geometry->vertex[val_marker][iVertex]->GetNormal(Normal); - for (auto iDim = 0ul; iDim < nDim; iDim++) Normal[iDim] = -Normal[iDim]; - conv_numerics->SetNormal(Normal); - - Area = GeometryToolbox::Norm(nDim, Normal); - - for (auto iDim = 0ul; iDim < nDim; iDim++) - UnitNormal[iDim] = Normal[iDim]/Area; - - /*--- Retrieve solution at this boundary node ---*/ - for (auto iVar = 0ul; iVar < nVar; iVar++) U_domain[iVar] = nodes->GetSolution(iPoint, iVar); - for (auto iVar = 0ul; iVar < nPrimVar; iVar++) V_domain[iVar] = nodes->GetPrimitive(iPoint,iVar); - - const su2double* dir = Inlet_FlowDir[val_marker][iVertex]; - const su2double mag = GeometryToolbox::Norm(nDim, dir); - if (config->GetInletUseNormal()) { - for (auto iDim = 0u; iDim < nDim; iDim++) - Flow_Dir[iDim] = -UnitNormal[iDim]; - } else { - for (auto iDim = 0u; iDim < nDim; iDim++) - Flow_Dir[iDim] = dir[iDim]/mag; - } - - /*--- Build the fictitious intlet state based on characteristics ---*/ - - /*--- Subsonic inflow: there is one outgoing characteristic (u-c), - therefore we can specify all but one state variable at the inlet. - The outgoing Riemann invariant provides the final piece of info. - Adapted from an original implementation in the Stanford University - multi-block (SUmb) solver in the routine bcSubsonicInflow.f90 - written by Edwin van der Weide, last modified 04-20-2009. ---*/ - - switch (Kind_Inlet) { - - /*--- Total properties have been specified at the inlet. ---*/ - case INLET_TYPE::TOTAL_CONDITIONS: - - /*--- Retrieve the specified total conditions for this inlet. ---*/ - P_Total = Inlet_Ptotal[val_marker][iVertex]; - T_Total = Inlet_Ttotal[val_marker][iVertex]; - - /*--- Non-dim. the inputs if necessary. ---*/ - P_Total /= config->GetPressure_Ref(); - T_Total /= config->GetTemperature_Ref(); - - /*--- Store primitives and set some variables for clarity. ---*/ - Density = V_domain[RHO_INDEX]; - for (auto iDim = 0ul; iDim < nDim; iDim++) - Velocity[iDim] = U_domain[nSpecies+iDim]/Density; - - Velocity2 = GeometryToolbox::SquaredNorm(nDim, Velocity); - Energy = U_domain[nVar-2]/Density; - Pressure = Gamma_Minus_One*Density*(Energy-0.5*Velocity2); - H_Total = (Gamma*Gas_Constant/Gamma_Minus_One)*T_Total; - SoundSpeed2 = Gamma*Pressure/Density; - - /*--- Compute the acoustic Riemann invariant that is extrapolated - from the domain interior. ---*/ - Riemann = 2.0*sqrt(SoundSpeed2)/Gamma_Minus_One; - for (auto iDim = 0ul; iDim < nDim; iDim++) - Riemann += Velocity[iDim]*UnitNormal[iDim]; - - /*--- Total speed of sound ---*/ - SoundSpeed_Total2 = Gamma_Minus_One*(H_Total - (Energy + Pressure/Density)+0.5*Velocity2) + SoundSpeed2; - - /*--- Dot product of normal and flow direction. This should - be negative due to outward facing boundary normal convention. ---*/ - alpha = GeometryToolbox::DotProduct(nDim, UnitNormal, Flow_Dir); - - /*--- Coefficients in the quadratic equation for the velocity ---*/ - aa = 1.0 + 0.5*Gamma_Minus_One*alpha*alpha; - bb = -1.0*Gamma_Minus_One*alpha*Riemann; - cc = 0.5*Gamma_Minus_One*Riemann*Riemann - -2.0*SoundSpeed_Total2/Gamma_Minus_One; - - /*--- Solve quadratic equation for velocity magnitude. Value must - be positive, so the choice of root is clear. ---*/ - dd = bb*bb - 4.0*aa*cc; - dd = sqrt(max(0.0,dd)); - Vel_Mag = (-bb + dd)/(2.0*aa); - Vel_Mag = max(0.0,Vel_Mag); - Velocity2 = Vel_Mag*Vel_Mag; - - /*--- Compute speed of sound from total speed of sound eqn. ---*/ - SoundSpeed2 = SoundSpeed_Total2 - 0.5*Gamma_Minus_One*Velocity2; - - /*--- Mach squared (cut between 0-1), use to adapt velocity ---*/ - Mach2 = Velocity2/SoundSpeed2; - Mach2 = min(1.0,Mach2); - Velocity2 = Mach2*SoundSpeed2; - Vel_Mag = sqrt(Velocity2); - SoundSpeed2 = SoundSpeed_Total2 - 0.5*Gamma_Minus_One*Velocity2; - - /*--- Compute new velocity vector at the inlet ---*/ - for (auto iDim = 0ul; iDim < nDim; iDim++) - Velocity[iDim] = Vel_Mag*Flow_Dir[iDim]; - - /*--- Static temperature from the speed of sound relation ---*/ - Temperature = SoundSpeed2/(Gamma*Gas_Constant); - //NEED TVE AS WELL - - /*--- Static pressure using isentropic relation at a point ---*/ - Pressure = P_Total*pow((Temperature/T_Total),Gamma/Gamma_Minus_One); - - /*--- Density at the inlet from the gas law ---*/ - Density = Pressure/(Gas_Constant*Temperature); - //NEED SPECIES DENSITIES - - /*--- Using pressure, density, & velocity, compute the energy ---*/ - Energy = Pressure/(Density*Gamma_Minus_One)+0.5*Velocity2; - //NEED EVE AS WELL - - /*--- Conservative variables, using the derived quantities ---*/ - for (auto iSpecies = 0ul; iSpeciesGetDensity_Ref(); - Vel_Mag /= config->GetVelocity_Ref(); - - /*--- Get primitives from current inlet state. ---*/ - for (auto iDim = 0ul; iDim < nDim; iDim++) - Velocity[iDim] = nodes->GetVelocity(iPoint, iDim); - Pressure = nodes->GetPressure(iPoint); - SoundSpeed2 = Gamma*Pressure/U_domain[0]; - - /*--- Compute the acoustic Riemann invariant that is extrapolated - from the domain interior. ---*/ - Riemann = Two_Gamma_M1*sqrt(SoundSpeed2); - for (auto iDim = 0ul; iDim < nDim; iDim++) - Riemann += Velocity[iDim]*UnitNormal[iDim]; - - /*--- Speed of sound squared for fictitious inlet state ---*/ - SoundSpeed2 = Riemann; - for (auto iDim = 0ul; iDim < nDim; iDim++) - SoundSpeed2 -= Vel_Mag*Flow_Dir[iDim]*UnitNormal[iDim]; - - SoundSpeed2 = max(0.0,0.5*Gamma_Minus_One*SoundSpeed2); - SoundSpeed2 = SoundSpeed2*SoundSpeed2; - - /*--- Pressure for the fictitious inlet state ---*/ - Pressure = SoundSpeed2*Density/Gamma; - - /*--- Energy for the fictitious inlet state ---*/ - Energy = Pressure/(Density*Gamma_Minus_One)+0.5*Vel_Mag*Vel_Mag; - - /*--- Conservative variables, using the derived quantities ---*/ - U_inlet[0] = Density; - for (auto iDim = 0ul; iDim < nDim; iDim++) - U_inlet[iDim+1] = Vel_Mag*Flow_Dir[iDim]*Density; - U_inlet[nDim+1] = Energy*Density; - - /*--- Primitive variables, using the derived quantities ---*/ - V_inlet[0] = Pressure / ( Gas_Constant * Density); - for (auto iDim = 0ul; iDim < nDim; iDim++) - V_inlet[iDim+1] = Vel_Mag*Flow_Dir[iDim]; - V_inlet[nDim+1] = Pressure; - V_inlet[nDim+2] = Density; - - break; - - default: - SU2_MPI::Error("Unsupported INLET_TYPE.", CURRENT_FUNCTION); - break; - } - - /*--- Set various quantities in the solver class ---*/ - conv_numerics->SetConservative(U_domain, U_inlet); - - if (dynamic_grid) - conv_numerics->SetGridVel(geometry->nodes->GetGridVel(iPoint), geometry->nodes->GetGridVel(iPoint)); - - /*--- Compute the residual using an upwind scheme ---*/ - auto residual = conv_numerics->ComputeResidual(config); - - LinSysRes.AddBlock(iPoint, residual); - - /*--- Jacobian contribution for implicit integration ---*/ - if (implicit) Jacobian.AddBlock2Diag(iPoint, residual.jacobian_i); - } - } - - /*--- Free locally allocated memory ---*/ - delete [] U_domain; - delete [] U_inlet; - delete [] V_domain; - delete [] V_inlet; - delete [] Normal; - delete [] Spec_Density; } void CNEMOEulerSolver::BC_Outlet(CGeometry *geometry, CSolver **solver_container, @@ -2138,55 +1905,8 @@ void CNEMOEulerSolver::BC_Supersonic_Inlet( Jacobian.AddBlock2Diag(iPoint, residual.jacobian_i); /*--- Viscous contribution ---*/ - const bool viscous = config->GetViscous(); - if (viscous) { - - /*--- Set the normal vector and the coordinates ---*/ - su2double Coord_Reflected[MAXNDIM]; - GeometryToolbox::PointPointReflect( - nDim, geometry->nodes->GetCoord(Point_Normal), - geometry->nodes->GetCoord(iPoint), Coord_Reflected); - visc_numerics->SetCoord(geometry->nodes->GetCoord(iPoint), - Coord_Reflected); - visc_numerics->SetNormal(Normal); - - /*--- Primitive variables, and gradient ---*/ - visc_numerics->SetConservative(nodes->GetSolution(iPoint), - node_inlet.GetSolution(0)); - visc_numerics->SetPrimitive(nodes->GetPrimitive(iPoint), - node_inlet.GetPrimitive(0)); - visc_numerics->SetPrimVarGradient(nodes->GetGradient_Primitive(iPoint), - node_inlet.GetGradient_Primitive(0)); - - /*--- Pass supplementary information to CNumerics ---*/ - SetNumericsThermoState(visc_numerics, nodes, iPoint, &node_inlet, 0); - - /*--- Species diffusion coefficients ---*/ - visc_numerics->SetDiffusionCoeff(nodes->GetDiffusionCoeff(iPoint), - nodes->GetDiffusionCoeff(iPoint)); - /*--- Laminar viscosity ---*/ - visc_numerics->SetLaminarViscosity(nodes->GetLaminarViscosity(iPoint), - nodes->GetLaminarViscosity(iPoint)); - /*--- Eddy viscosity ---*/ - visc_numerics->SetEddyViscosity(nodes->GetEddyViscosity(iPoint), - nodes->GetEddyViscosity(iPoint)); - /*--- Thermal conductivity ---*/ - visc_numerics->SetThermalConductivity( - nodes->GetThermalConductivity(iPoint), - nodes->GetThermalConductivity(iPoint)); - /*--- Vib-el. thermal conductivity ---*/ - visc_numerics->SetThermalConductivity_ve( - nodes->GetThermalConductivity_ve(iPoint), - nodes->GetThermalConductivity_ve(iPoint)); - - /*--- Compute and update residual ---*/ - auto residual = visc_numerics->ComputeResidual(config); - LinSysRes.SubtractBlock(iPoint, residual); - - /*--- Jacobian contribution for implicit integration ---*/ - if (implicit) - Jacobian.SubtractBlock2Diag(iPoint, residual.jacobian_i); - } + if (config->GetViscous()) + AddViscousBoundaryResidual(geometry, visc_numerics, config, iPoint, Point_Normal, Normal, &node_inlet, 0); } } diff --git a/SU2_CFD/src/solvers/CNEMONSSolver.cpp b/SU2_CFD/src/solvers/CNEMONSSolver.cpp index aa419fc9dd4b..71213d44e3d2 100644 --- a/SU2_CFD/src/solvers/CNEMONSSolver.cpp +++ b/SU2_CFD/src/solvers/CNEMONSSolver.cpp @@ -201,25 +201,8 @@ void CNEMONSSolver::Viscous_Residual(CGeometry *geometry, /*--- Pass supplementary information to CNumerics ---*/ SetNumericsThermoState(numerics, nodes, iPoint, nodes, jPoint); - /*--- Species diffusion coefficients ---*/ - numerics->SetDiffusionCoeff(nodes->GetDiffusionCoeff(iPoint), - nodes->GetDiffusionCoeff(jPoint) ); - - /*--- Laminar viscosity ---*/ - numerics->SetLaminarViscosity(nodes->GetLaminarViscosity(iPoint), - nodes->GetLaminarViscosity(jPoint) ); - - /*--- Eddy viscosity ---*/ - numerics->SetEddyViscosity(nodes->GetEddyViscosity(iPoint), - nodes->GetEddyViscosity(jPoint) ); - - /*--- Thermal conductivity ---*/ - numerics->SetThermalConductivity(nodes->GetThermalConductivity(iPoint), - nodes->GetThermalConductivity(jPoint)); - - /*--- Vib-el. thermal conductivity ---*/ - numerics->SetThermalConductivity_ve(nodes->GetThermalConductivity_ve(iPoint), - nodes->GetThermalConductivity_ve(jPoint) ); + /*--- Transport properties ---*/ + SetNumericsTransportState(numerics, iPoint, jPoint); /*--- Compute and update residual ---*/ auto residual = numerics->ComputeResidual(config); From 7d5bd07ddb8f065612cdfeb47612f29123c09623 Mon Sep 17 00:00:00 2001 From: rois1995 Date: Mon, 12 Oct 2026 12:04:51 +0200 Subject: [PATCH 6/9] Use static arrays for the work arrays of the NEMO boundary conditions BC_Outlet and the isothermal walls allocated and freed their work arrays with new and delete at every call. Use static arrays of size MAXNVAR and MAXNDIM. The Jacobian and the unit normal of the isothermal wall are now local to each vertex of the loop. No functional change. Co-Authored-By: Claude Opus 5.5 Claude-Session: https://claude.ai/code/session_01L7oxzVir4zENrpy3K3W6Fm --- SU2_CFD/src/solvers/CNEMOEulerSolver.cpp | 12 +++------ SU2_CFD/src/solvers/CNEMONSSolver.cpp | 34 +++--------------------- 2 files changed, 7 insertions(+), 39 deletions(-) diff --git a/SU2_CFD/src/solvers/CNEMOEulerSolver.cpp b/SU2_CFD/src/solvers/CNEMOEulerSolver.cpp index 27fc8e43d251..597cce3f02b4 100644 --- a/SU2_CFD/src/solvers/CNEMOEulerSolver.cpp +++ b/SU2_CFD/src/solvers/CNEMOEulerSolver.cpp @@ -1637,10 +1637,10 @@ void CNEMOEulerSolver::BC_Outlet(CGeometry *geometry, CSolver **solver_container const bool implicit = (config->GetKind_TimeIntScheme() == EULER_IMPLICIT); su2double *U_domain; - auto *U_outlet = new su2double[nVar]; + su2double U_outlet[MAXNVAR] = {0.0}; su2double *V_domain, *V_outlet; - auto *Normal = new su2double[nDim]; - auto *Ys = new su2double[nSpecies]; + su2double Normal[MAXNDIM] = {0.0}; + su2double Ys[MAXNVAR] = {0.0}; const unsigned short T_INDEX = nodes->GetTIndex(); const unsigned short TVE_INDEX = nodes->GetTveIndex(); @@ -1813,12 +1813,6 @@ void CNEMOEulerSolver::BC_Outlet(CGeometry *geometry, CSolver **solver_container Jacobian.AddBlock2Diag(iPoint, residual.jacobian_i); } } - - /*--- Free locally allocated memory ---*/ - delete [] U_outlet; - delete [] Normal; - delete [] Ys; - } void CNEMOEulerSolver::BC_Supersonic_Inlet( diff --git a/SU2_CFD/src/solvers/CNEMONSSolver.cpp b/SU2_CFD/src/solvers/CNEMONSSolver.cpp index 71213d44e3d2..3a5bc0431a6d 100644 --- a/SU2_CFD/src/solvers/CNEMONSSolver.cpp +++ b/SU2_CFD/src/solvers/CNEMONSSolver.cpp @@ -388,7 +388,6 @@ void CNEMONSSolver::BC_IsothermalNonCatalytic_Wall(CGeometry *geometry, const bool implicit = (config->GetKind_TimeIntScheme() == EULER_IMPLICIT); const bool ionization = config->GetIonization(); - su2double UnitNormal[MAXNDIM] = {0.0}; if (ionization) { SU2_MPI::Error("NEED TO TAKE A CLOSER LOOK AT THE JACOBIAN W/ IONIZATION",CURRENT_FUNCTION); @@ -403,13 +402,6 @@ void CNEMONSSolver::BC_IsothermalNonCatalytic_Wall(CGeometry *geometry, /*--- Retrieve the specified wall temperature ---*/ const su2double Twall = config->GetIsothermal_Temperature(Marker_Tag); - su2double **Jacobian_i = nullptr; - if (implicit) { - Jacobian_i = new su2double* [nVar]; - for (auto iVar = 0ul; iVar < nVar; iVar++) - Jacobian_i[iVar] = new su2double [nVar] (); - } - /*--- Loop over boundary points to calculate energy flux ---*/ SU2_OMP_FOR_DYN(OMP_MIN_SIZE) for (auto iVertex = 0ul; iVertex < geometry->nVertex[val_marker]; iVertex++) { @@ -421,6 +413,7 @@ void CNEMONSSolver::BC_IsothermalNonCatalytic_Wall(CGeometry *geometry, /*--- Compute dual-grid area and boundary normal ---*/ const auto Normal = geometry->vertex[val_marker][iVertex]->GetNormal(); const su2double Area = GeometryToolbox::Norm(nDim, Normal); + su2double UnitNormal[MAXNDIM] = {0.0}; for (auto iDim = 0ul; iDim < nDim; iDim++) UnitNormal[iDim] = Normal[iDim]/Area; @@ -449,13 +442,9 @@ void CNEMONSSolver::BC_IsothermalNonCatalytic_Wall(CGeometry *geometry, SetWallEnergyFlux_Weak(iPoint, Point_Normal, Twall, Twall, C, Area, dist_ij, Res_Visc); /*--- Calculate Jacobian for implicit time stepping ---*/ + su2double Jacobian_i[MAXNVAR][MAXNVAR] = {{0.0}}; if (implicit) { - /*--- Initialize Jacobian to zero ---*/ - for (auto iVar = 0ul; iVar < nVar; iVar++) - for (auto jVar = 0ul; jVar < nVar; jVar++) - Jacobian_i[iVar][jVar] = 0.0; - /*--- Add contributions to the Jacobian from the weak enforcement of the energy equations. ---*/ const auto dTdU = nodes->GetdTdU(iPoint); const auto dTvedU = nodes->GetdTvedU(iPoint); @@ -483,12 +472,6 @@ void CNEMONSSolver::BC_IsothermalNonCatalytic_Wall(CGeometry *geometry, } } END_SU2_OMP_FOR - - if (Jacobian_i) - for (auto iVar = 0ul; iVar < nVar; iVar++) - delete [] Jacobian_i[iVar]; - delete [] Jacobian_i; - } void CNEMONSSolver::BC_IsothermalCatalytic_Wall(CGeometry *geometry, @@ -503,8 +486,8 @@ void CNEMONSSolver::BC_IsothermalCatalytic_Wall(CGeometry *geometry, BC_IsothermalNonCatalytic_Wall(geometry, solver_container, conv_numerics, sour_numerics, config, val_marker); - /*--- Local variables ---*/ - su2double **dVdU; + /*--- Transformation matrix from conservative to primitive variables ---*/ + su2double dVdU[MAXNVAR][MAXNVAR] = {{0.0}}; /*--- Assign booleans ---*/ const bool implicit = (config->GetKind_TimeIntScheme() == EULER_IMPLICIT); @@ -521,11 +504,6 @@ void CNEMONSSolver::BC_IsothermalCatalytic_Wall(CGeometry *geometry, const unsigned short T_INDEX = nodes->GetTIndex(); const unsigned short TVE_INDEX = nodes->GetTveIndex(); - /*--- Allocate arrays ---*/ - dVdU = new su2double*[nVar]; - for (auto iVar = 0ul; iVar < nVar; iVar++) - dVdU[iVar] = new su2double[nVar]; - /*--- Loop over all of the vertices on this boundary marker ---*/ for (auto iVertex = 0ul; iVertex < geometry->nVertex[val_marker]; iVertex++) { @@ -675,10 +653,6 @@ void CNEMONSSolver::BC_IsothermalCatalytic_Wall(CGeometry *geometry, LinSysRes.SubtractBlock(iPoint, Res_Visc); } } - - for (auto iVar = 0ul; iVar < nVar; iVar++) - delete [] dVdU[iVar]; - delete [] dVdU; } void CNEMONSSolver::BC_Smoluchowski_Maxwell(CGeometry *geometry, From 222b7b64d43eaf237041302e1fe43aeec57d765e Mon Sep 17 00:00:00 2001 From: rois1995 Date: Mon, 12 Oct 2026 12:47:15 +0200 Subject: [PATCH 7/9] Share the reference values of the compressible solvers CEulerSolver and CNEMOEulerSolver computed the reference values and the non-dimensional free-stream values with the same statements in their SetNondimensionalization. Move them to SetCompressibleNonDimValues in CFVMFlowSolverBase. Each solver keeps the values that are specific to it (vib.-el. temperature for NEMO; conductivity, specific heat, ReThetaT and mass diffusivity for the ideal gas solver). No functional change. Co-Authored-By: Claude Opus 5.5 Claude-Session: https://claude.ai/code/session_01L7oxzVir4zENrpy3K3W6Fm --- .../include/solvers/CFVMFlowSolverBase.hpp | 16 +++ .../include/solvers/CFVMFlowSolverBase.inl | 108 ++++++++++++++++++ SU2_CFD/src/solvers/CEulerSolver.cpp | 102 ++++------------- SU2_CFD/src/solvers/CNEMOEulerSolver.cpp | 100 +++------------- 4 files changed, 165 insertions(+), 161 deletions(-) diff --git a/SU2_CFD/include/solvers/CFVMFlowSolverBase.hpp b/SU2_CFD/include/solvers/CFVMFlowSolverBase.hpp index ae07bef84552..a9b849583135 100644 --- a/SU2_CFD/include/solvers/CFVMFlowSolverBase.hpp +++ b/SU2_CFD/include/solvers/CFVMFlowSolverBase.hpp @@ -411,6 +411,22 @@ class CFVMFlowSolverBase : public CSolver { void LoadRestart_impl(CGeometry **geometry, CSolver ***solver, CConfig *config, int iter, bool update_geo, su2double* RestartSolution = nullptr, unsigned short nVar_Restart = 0); + /*! + * \brief Compute the reference values and the non-dimensional free-stream values that are common + * to the compressible solvers, and store them in the config. + * \note The reference pressure, density and temperature follow from the option REF_DIMENSIONALIZATION. + * All other reference values are derived from these three, with a reference length of one. + * \param[in,out] config - Definition of the particular problem. + * \param[in] Pressure_FreeStream - Dimensional free-stream pressure. + * \param[in] Density_FreeStream - Dimensional free-stream density. + * \param[in] Temperature_FreeStream - Dimensional free-stream temperature. + * \param[in] Viscosity_FreeStream - Dimensional free-stream laminar viscosity. + * \param[in] ModVel_FreeStream - Dimensional magnitude of the free-stream velocity. + */ + void SetCompressibleNonDimValues(CConfig *config, su2double Pressure_FreeStream, su2double Density_FreeStream, + su2double Temperature_FreeStream, su2double Viscosity_FreeStream, + su2double ModVel_FreeStream) const; + /*! * \brief Generic implementation to compute the time step based on CFL and conv/visc eigenvalues. * \param[in] geometry - Geometrical definition of the problem. diff --git a/SU2_CFD/include/solvers/CFVMFlowSolverBase.inl b/SU2_CFD/include/solvers/CFVMFlowSolverBase.inl index cb44757c6280..e7f38e23e99b 100644 --- a/SU2_CFD/include/solvers/CFVMFlowSolverBase.inl +++ b/SU2_CFD/include/solvers/CFVMFlowSolverBase.inl @@ -850,6 +850,114 @@ void CFVMFlowSolverBase::UpdateCustomBoundaryConditions( } } +template +void CFVMFlowSolverBase::SetCompressibleNonDimValues(CConfig *config, su2double Pressure_FreeStream, + su2double Density_FreeStream, + su2double Temperature_FreeStream, + su2double Viscosity_FreeStream, + su2double ModVel_FreeStream) const { + + const su2double Mach = config->GetMach(); + + /*--- Reference pressure, density and temperature. ---*/ + + su2double Pressure_Ref = 0.0, Density_Ref = 0.0, Temperature_Ref = 0.0; + + if (config->GetRef_NonDim() == DIMENSIONAL) { + Pressure_Ref = 1.0; + Density_Ref = 1.0; + Temperature_Ref = 1.0; + } + else if (config->GetRef_NonDim() == FREESTREAM_PRESS_EQ_ONE) { + Pressure_Ref = Pressure_FreeStream; // Pressure_FreeStream = 1.0 + Density_Ref = Density_FreeStream; // Density_FreeStream = 1.0 + Temperature_Ref = Temperature_FreeStream; // Temperature_FreeStream = 1.0 + } + else if (config->GetRef_NonDim() == FREESTREAM_VEL_EQ_MACH) { + Pressure_Ref = Gamma*Pressure_FreeStream; // Pressure_FreeStream = 1.0/Gamma + Density_Ref = Density_FreeStream; // Density_FreeStream = 1.0 + Temperature_Ref = Temperature_FreeStream; // Temp_FreeStream = 1.0 + } + else if (config->GetRef_NonDim() == FREESTREAM_VEL_EQ_ONE) { + Pressure_Ref = Mach*Mach*Gamma*Pressure_FreeStream; // Pressure_FreeStream = 1.0/(Gamma*(M_inf)^2) + Density_Ref = Density_FreeStream; // Density_FreeStream = 1.0 + Temperature_Ref = Temperature_FreeStream; // Temp_FreeStream = 1.0 + } + config->SetPressure_Ref(Pressure_Ref); + config->SetDensity_Ref(Density_Ref); + config->SetTemperature_Ref(Temperature_Ref); + + /*--- Reference values that are derived from the three above. ---*/ + + const su2double Length_Ref = 1.0; + const su2double Velocity_Ref = sqrt(config->GetPressure_Ref()/config->GetDensity_Ref()); + const su2double Time_Ref = Length_Ref/Velocity_Ref; + const su2double Omega_Ref = Velocity_Ref/Length_Ref; + const su2double Force_Ref = config->GetDensity_Ref()*Velocity_Ref*Velocity_Ref*Length_Ref*Length_Ref; + const su2double Heat_Flux_Ref = Density_Ref*Velocity_Ref*Velocity_Ref*Velocity_Ref; + const su2double Gas_Constant_Ref = Velocity_Ref*Velocity_Ref/config->GetTemperature_Ref(); + const su2double Viscosity_Ref = config->GetDensity_Ref()*Velocity_Ref*Length_Ref; + const su2double Conductivity_Ref = Viscosity_Ref*Gas_Constant_Ref; + const su2double Froude = ModVel_FreeStream/sqrt(STANDARD_GRAVITY*Length_Ref); + + config->SetLength_Ref(Length_Ref); + config->SetVelocity_Ref(Velocity_Ref); + config->SetTime_Ref(Time_Ref); + config->SetOmega_Ref(Omega_Ref); + config->SetForce_Ref(Force_Ref); + config->SetHeat_Flux_Ref(Heat_Flux_Ref); + config->SetGas_Constant_Ref(Gas_Constant_Ref); + config->SetViscosity_Ref(Viscosity_Ref); + config->SetConductivity_Ref(Conductivity_Ref); + config->SetFroude(Froude); + + /*--- Divide by reference values, to compute the non-dimensional free-stream values. ---*/ + + const su2double Pressure_FreeStreamND = Pressure_FreeStream/config->GetPressure_Ref(); + config->SetPressure_FreeStreamND(Pressure_FreeStreamND); + + const su2double Density_FreeStreamND = Density_FreeStream/config->GetDensity_Ref(); + config->SetDensity_FreeStreamND(Density_FreeStreamND); + + su2double Velocity_FreeStreamND[MAXNDIM] = {0.0}; + for (auto iDim = 0u; iDim < nDim; iDim++) { + Velocity_FreeStreamND[iDim] = config->GetVelocity_FreeStream()[iDim]/Velocity_Ref; + config->SetVelocity_FreeStreamND(Velocity_FreeStreamND[iDim], iDim); + } + + const su2double Temperature_FreeStreamND = Temperature_FreeStream/config->GetTemperature_Ref(); + config->SetTemperature_FreeStreamND(Temperature_FreeStreamND); + + const su2double Gas_ConstantND = config->GetGas_Constant()/Gas_Constant_Ref; + config->SetGas_ConstantND(Gas_ConstantND); + + su2double ModVel_FreeStreamND = 0.0; + for (auto iDim = 0u; iDim < nDim; iDim++) ModVel_FreeStreamND += Velocity_FreeStreamND[iDim]*Velocity_FreeStreamND[iDim]; + ModVel_FreeStreamND = sqrt(ModVel_FreeStreamND); + config->SetModVel_FreeStreamND(ModVel_FreeStreamND); + + const su2double Viscosity_FreeStreamND = Viscosity_FreeStream / Viscosity_Ref; + config->SetViscosity_FreeStreamND(Viscosity_FreeStreamND); + + /*--- Free-stream turbulence quantities, from the turbulence intensity and the viscosity ratio. ---*/ + + const su2double TurbIntensity = config->GetTurbulenceIntensity_FreeStream(); + + const su2double Tke_FreeStream = 3.0/2.0*(ModVel_FreeStream*ModVel_FreeStream*TurbIntensity*TurbIntensity); + config->SetTke_FreeStream(Tke_FreeStream); + + const su2double Tke_FreeStreamND = 3.0/2.0*(ModVel_FreeStreamND*ModVel_FreeStreamND*TurbIntensity*TurbIntensity); + config->SetTke_FreeStreamND(Tke_FreeStreamND); + + const su2double Omega_FreeStream = + Density_FreeStream*Tke_FreeStream/max(Viscosity_FreeStream*config->GetTurb2LamViscRatio_FreeStream(), EPS); + config->SetOmega_FreeStream(Omega_FreeStream); + + const su2double Omega_FreeStreamND = + Density_FreeStreamND*Tke_FreeStreamND/max(Viscosity_FreeStreamND*config->GetTurb2LamViscRatio_FreeStream(), EPS); + config->SetOmega_FreeStreamND(Omega_FreeStreamND); +} + template void CFVMFlowSolverBase::LoadRestart_impl(CGeometry **geometry, CSolver ***solver, CConfig *config, int iter, bool update_geo, su2double* SolutionRestart, diff --git a/SU2_CFD/src/solvers/CEulerSolver.cpp b/SU2_CFD/src/solvers/CEulerSolver.cpp index 32c265e55656..abb21e51c728 100644 --- a/SU2_CFD/src/solvers/CEulerSolver.cpp +++ b/SU2_CFD/src/solvers/CEulerSolver.cpp @@ -807,19 +807,11 @@ void CEulerSolver::SetNondimensionalization(CConfig *config, unsigned short iMes SU2_ZONE_SCOPED su2double Temperature_FreeStream = 0.0, Mach2Vel_FreeStream = 0.0, ModVel_FreeStream = 0.0, - Energy_FreeStream = 0.0, ModVel_FreeStreamND = 0.0, Velocity_Reynolds = 0.0, - Omega_FreeStream = 0.0, Omega_FreeStreamND = 0.0, Viscosity_FreeStream = 0.0, + Energy_FreeStream = 0.0, Velocity_Reynolds = 0.0, Viscosity_FreeStream = 0.0, Thermal_Conductivity_FreeStream = 0.0, SpecificHeat_Cp_FreeStream = 0.0, Density_FreeStream = 0.0, Pressure_FreeStream = 0.0, Tke_FreeStream = 0.0, Re_ThetaT_FreeStream = 0.0, - Length_Ref = 0.0, Density_Ref = 0.0, Pressure_Ref = 0.0, Velocity_Ref = 0.0, - Temperature_Ref = 0.0, Time_Ref = 0.0, Omega_Ref = 0.0, Force_Ref = 0.0, - Gas_Constant_Ref = 0.0, Viscosity_Ref = 0.0, Conductivity_Ref = 0.0, Energy_Ref= 0.0, - Froude = 0.0, Pressure_FreeStreamND = 0.0, Density_FreeStreamND = 0.0, - Temperature_FreeStreamND = 0.0, Gas_ConstantND = 0.0, - Velocity_FreeStreamND[3] = {0.0, 0.0, 0.0}, Viscosity_FreeStreamND = 0.0, - Thermal_Conductivity_FreeStreamND = 0.0, SpecificHeat_Cp_FreeStreamND = 0.0, - Tke_FreeStreamND = 0.0, Energy_FreeStreamND = 0.0, - Total_UnstTimeND = 0.0, Delta_UnstTimeND = 0.0, TgammaR = 0.0, Heat_Flux_Ref = 0.0; + Energy_Ref= 0.0, Thermal_Conductivity_FreeStreamND = 0.0, SpecificHeat_Cp_FreeStreamND = 0.0, + Energy_FreeStreamND = 0.0, Total_UnstTimeND = 0.0, Delta_UnstTimeND = 0.0, TgammaR = 0.0; unsigned short iDim; @@ -1023,80 +1015,20 @@ void CEulerSolver::SetNondimensionalization(CConfig *config, unsigned short iMes if (tkeNeeded) { Energy_FreeStream += Tke_FreeStream; }; config->SetEnergy_FreeStream(Energy_FreeStream); - /*--- Compute non dimensional quantities. By definition, - Lref is one because we have converted the grid to meters. ---*/ + /*--- Compute the reference values and the non-dimensional free-stream values, + and store them in the config. By definition, Lref is one because we have + converted the grid to meters. ---*/ - if (config->GetRef_NonDim() == DIMENSIONAL) { - Pressure_Ref = 1.0; - Density_Ref = 1.0; - Temperature_Ref = 1.0; - } - else if (config->GetRef_NonDim() == FREESTREAM_PRESS_EQ_ONE) { - Pressure_Ref = Pressure_FreeStream; // Pressure_FreeStream = 1.0 - Density_Ref = Density_FreeStream; // Density_FreeStream = 1.0 - Temperature_Ref = Temperature_FreeStream; // Temperature_FreeStream = 1.0 - } - else if (config->GetRef_NonDim() == FREESTREAM_VEL_EQ_MACH) { - Pressure_Ref = Gamma*Pressure_FreeStream; // Pressure_FreeStream = 1.0/Gamma - Density_Ref = Density_FreeStream; // Density_FreeStream = 1.0 - Temperature_Ref = Temperature_FreeStream; // Temp_FreeStream = 1.0 - } - else if (config->GetRef_NonDim() == FREESTREAM_VEL_EQ_ONE) { - Pressure_Ref = Mach*Mach*Gamma*Pressure_FreeStream; // Pressure_FreeStream = 1.0/(Gamma*(M_inf)^2) - Density_Ref = Density_FreeStream; // Density_FreeStream = 1.0 - Temperature_Ref = Temperature_FreeStream; // Temp_FreeStream = 1.0 - } - config->SetPressure_Ref(Pressure_Ref); - config->SetDensity_Ref(Density_Ref); - config->SetTemperature_Ref(Temperature_Ref); - - Length_Ref = 1.0; config->SetLength_Ref(Length_Ref); - Velocity_Ref = sqrt(config->GetPressure_Ref()/config->GetDensity_Ref()); config->SetVelocity_Ref(Velocity_Ref); - Time_Ref = Length_Ref/Velocity_Ref; config->SetTime_Ref(Time_Ref); - Omega_Ref = Velocity_Ref/Length_Ref; config->SetOmega_Ref(Omega_Ref); - Force_Ref = config->GetDensity_Ref()*Velocity_Ref*Velocity_Ref*Length_Ref*Length_Ref; config->SetForce_Ref(Force_Ref); - Heat_Flux_Ref = Density_Ref*Velocity_Ref*Velocity_Ref*Velocity_Ref; config->SetHeat_Flux_Ref(Heat_Flux_Ref); - Gas_Constant_Ref = Velocity_Ref*Velocity_Ref/config->GetTemperature_Ref(); config->SetGas_Constant_Ref(Gas_Constant_Ref); - Viscosity_Ref = config->GetDensity_Ref()*Velocity_Ref*Length_Ref; config->SetViscosity_Ref(Viscosity_Ref); - Conductivity_Ref = Viscosity_Ref*Gas_Constant_Ref; config->SetConductivity_Ref(Conductivity_Ref); - Froude = ModVel_FreeStream/sqrt(STANDARD_GRAVITY*Length_Ref); config->SetFroude(Froude); - - /*--- Divide by reference values, to compute the non-dimensional free-stream values ---*/ - - Pressure_FreeStreamND = Pressure_FreeStream/config->GetPressure_Ref(); config->SetPressure_FreeStreamND(Pressure_FreeStreamND); - Density_FreeStreamND = Density_FreeStream/config->GetDensity_Ref(); config->SetDensity_FreeStreamND(Density_FreeStreamND); - - for (iDim = 0; iDim < nDim; iDim++) { - Velocity_FreeStreamND[iDim] = config->GetVelocity_FreeStream()[iDim]/Velocity_Ref; config->SetVelocity_FreeStreamND(Velocity_FreeStreamND[iDim], iDim); - } - - Temperature_FreeStreamND = Temperature_FreeStream/config->GetTemperature_Ref(); config->SetTemperature_FreeStreamND(Temperature_FreeStreamND); + SetCompressibleNonDimValues(config, Pressure_FreeStream, Density_FreeStream, Temperature_FreeStream, + Viscosity_FreeStream, ModVel_FreeStream); - Gas_ConstantND = config->GetGas_Constant()/Gas_Constant_Ref; config->SetGas_ConstantND(Gas_ConstantND); + /*--- Non-dimensional free-stream values that are specific to this solver ---*/ - - ModVel_FreeStreamND = 0.0; - for (iDim = 0; iDim < nDim; iDim++) ModVel_FreeStreamND += Velocity_FreeStreamND[iDim]*Velocity_FreeStreamND[iDim]; - ModVel_FreeStreamND = sqrt(ModVel_FreeStreamND); config->SetModVel_FreeStreamND(ModVel_FreeStreamND); - - Viscosity_FreeStreamND = Viscosity_FreeStream / Viscosity_Ref; config->SetViscosity_FreeStreamND(Viscosity_FreeStreamND); - Thermal_Conductivity_FreeStreamND = Thermal_Conductivity_FreeStream / Conductivity_Ref; + Thermal_Conductivity_FreeStreamND = Thermal_Conductivity_FreeStream / config->GetThermal_Conductivity_Ref(); config->SetThermalConductivity_FreeStreamND(Thermal_Conductivity_FreeStreamND); - SpecificHeat_Cp_FreeStreamND = SpecificHeat_Cp_FreeStream / Gas_Constant_Ref; + SpecificHeat_Cp_FreeStreamND = SpecificHeat_Cp_FreeStream / config->GetGas_Constant_Ref(); config->SetSpecificHeatCp_FreeStreamND(SpecificHeat_Cp_FreeStreamND); - Tke_FreeStream = 3.0/2.0*(ModVel_FreeStream*ModVel_FreeStream*config->GetTurbulenceIntensity_FreeStream()*config->GetTurbulenceIntensity_FreeStream()); - config->SetTke_FreeStream(Tke_FreeStream); - - Tke_FreeStreamND = 3.0/2.0*(ModVel_FreeStreamND*ModVel_FreeStreamND*config->GetTurbulenceIntensity_FreeStream()*config->GetTurbulenceIntensity_FreeStream()); - config->SetTke_FreeStreamND(Tke_FreeStreamND); - - Omega_FreeStream = Density_FreeStream*Tke_FreeStream/max(Viscosity_FreeStream*config->GetTurb2LamViscRatio_FreeStream(), EPS); - config->SetOmega_FreeStream(Omega_FreeStream); - - Omega_FreeStreamND = Density_FreeStreamND*Tke_FreeStreamND/max(Viscosity_FreeStreamND*config->GetTurb2LamViscRatio_FreeStream(), EPS); - config->SetOmega_FreeStreamND(Omega_FreeStreamND); - if (config->GetTurbulenceIntensity_FreeStream() *100 <= 1.3) { if (config->GetTurbulenceIntensity_FreeStream() *100 >=0.027) { Re_ThetaT_FreeStream = (1173.51-589.428*config->GetTurbulenceIntensity_FreeStream() *100+0.2196/ @@ -1111,9 +1043,19 @@ void CEulerSolver::SetNondimensionalization(CConfig *config, unsigned short iMes } config->SetReThetaT_FreeStream(Re_ThetaT_FreeStream); - const su2double MassDiffusivityND = config->GetDiffusivity_Constant() / (Velocity_Ref * Length_Ref); + const su2double MassDiffusivityND = config->GetDiffusivity_Constant() / (config->GetVelocity_Ref() * config->GetLength_Ref()); config->SetDiffusivity_ConstantND(MassDiffusivityND); + /*--- Retrieve the values that are needed below ---*/ + + const su2double Time_Ref = config->GetTime_Ref(); + const su2double Froude = config->GetFroude(); + const su2double Gas_ConstantND = config->GetGas_ConstantND(); + const su2double Pressure_FreeStreamND = config->GetPressure_FreeStreamND(); + const su2double Density_FreeStreamND = config->GetDensity_FreeStreamND(); + const su2double ModVel_FreeStreamND = config->GetModVel_FreeStreamND(); + const su2double Tke_FreeStreamND = config->GetTke_FreeStreamND(); + /*--- Initialize the dimensionless Fluid Model that will be used to solve the dimensionless problem ---*/ /*--- Auxilary (dimensional) FluidModel no longer needed. ---*/ diff --git a/SU2_CFD/src/solvers/CNEMOEulerSolver.cpp b/SU2_CFD/src/solvers/CNEMOEulerSolver.cpp index 597cce3f02b4..41b7d980e93f 100644 --- a/SU2_CFD/src/solvers/CNEMOEulerSolver.cpp +++ b/SU2_CFD/src/solvers/CNEMOEulerSolver.cpp @@ -983,21 +983,11 @@ void CNEMOEulerSolver::SetNondimensionalization(CConfig *config, unsigned short su2double Temperature_FreeStream = 0.0, Temperature_ve_FreeStream = 0.0, Mach2Vel_FreeStream = 0.0, - ModVel_FreeStream = 0.0, Energy_FreeStream = 0.0, ModVel_FreeStreamND = 0.0, - Velocity_Reynolds = 0.0, Omega_FreeStream = 0.0, Omega_FreeStreamND = 0.0, + ModVel_FreeStream = 0.0, Energy_FreeStream = 0.0, Velocity_Reynolds = 0.0, Viscosity_FreeStream = 0.0, Density_FreeStream = 0.0, Pressure_FreeStream = 0.0, - Tke_FreeStream = 0.0, Length_Ref = 0.0, Density_Ref = 0.0, - Pressure_Ref = 0.0, Velocity_Ref = 0.0, Temperature_Ref = 0.0, - Temperature_ve_Ref = 0.0, Time_Ref = 0.0, Omega_Ref = 0.0, - Force_Ref = 0.0, Gas_Constant_Ref = 0.0, Viscosity_Ref = 0.0, - Conductivity_Ref = 0.0, Energy_Ref = 0.0, Pressure_FreeStreamND = 0.0, - Energy_FreeStreamND = 0.0, Temperature_FreeStreamND = 0.0, Temperature_ve_FreeStreamND = 0.0, - Gas_ConstantND = 0.0, Viscosity_FreeStreamND = 0.0, sqvel = 0.0, - Tke_FreeStreamND = 0.0, Total_UnstTimeND = 0.0, Delta_UnstTimeND = 0.0, - soundspeed = 0.0, GasConstant_Inf = 0.0, Froude = 0.0, - Density_FreeStreamND = 0.0, Heat_Flux_Ref = 0.0; - - su2double Velocity_FreeStreamND[3] = {0.0, 0.0, 0.0}; + Tke_FreeStream = 0.0, Energy_Ref = 0.0, Energy_FreeStreamND = 0.0, + sqvel = 0.0, Total_UnstTimeND = 0.0, Delta_UnstTimeND = 0.0, + soundspeed = 0.0, GasConstant_Inf = 0.0; /*--- Local variables ---*/ const su2double Alpha= config->GetAoA()*PI_NUMBER/180.0; @@ -1121,79 +1111,27 @@ void CNEMOEulerSolver::SetNondimensionalization(CConfig *config, unsigned short /*-- Compute the freestream energy. ---*/ if (tkeNeeded) { Energy_FreeStream += Tke_FreeStream; }; config->SetEnergy_FreeStream(Energy_FreeStream); - /*--- Compute non dimensional quantities. By definition, - Lref is one because we have converted the grid to meters. ---*/ + /*--- Compute the reference values and the non-dimensional free-stream values, + and store them in the config. By definition, Lref is one because we have + converted the grid to meters. ---*/ - if (config->GetRef_NonDim() == DIMENSIONAL) { - Pressure_Ref = 1.0; - Density_Ref = 1.0; - Temperature_Ref = 1.0; - Temperature_ve_Ref = 1.0; - } - else if (config->GetRef_NonDim() == FREESTREAM_PRESS_EQ_ONE) { - Pressure_Ref = Pressure_FreeStream; // Pressure_FreeStream = 1.0 - Density_Ref = Density_FreeStream; // Density_FreeStream = 1.0 - Temperature_Ref = Temperature_FreeStream; // Temperature_FreeStream = 1.0 - Temperature_ve_Ref = Temperature_ve_FreeStream; // Temperature_ve_FreeStream = 1.0 - } - else if (config->GetRef_NonDim() == FREESTREAM_VEL_EQ_MACH) { - Pressure_Ref = Gamma*Pressure_FreeStream; // Pressure_FreeStream = 1.0/Gamma - Density_Ref = Density_FreeStream; // Density_FreeStream = 1.0 - Temperature_Ref = Temperature_FreeStream; // Temp_FreeStream = 1.0 - Temperature_ve_Ref = Temperature_ve_FreeStream; // Temp_ve_FreeStream = 1.0 - } - else if (config->GetRef_NonDim() == FREESTREAM_VEL_EQ_ONE) { - Pressure_Ref = Mach*Mach*Gamma*Pressure_FreeStream; // Pressure_FreeStream = 1.0/(Gamma*(M_inf)^2) - Density_Ref = Density_FreeStream; // Density_FreeStream = 1.0 - Temperature_Ref = Temperature_FreeStream; // Temp_FreeStream = 1.0 - Temperature_ve_Ref = Temperature_ve_FreeStream; // Temp_ve_FreeStream = 1.0 - } - config->SetPressure_Ref(Pressure_Ref); - config->SetDensity_Ref(Density_Ref); - config->SetTemperature_Ref(Temperature_Ref); - config->SetTemperature_ve_Ref(Temperature_ve_Ref); - - Length_Ref = 1.0; config->SetLength_Ref(Length_Ref); - Velocity_Ref = sqrt(config->GetPressure_Ref()/config->GetDensity_Ref()); config->SetVelocity_Ref(Velocity_Ref); - Time_Ref = Length_Ref/Velocity_Ref; config->SetTime_Ref(Time_Ref); - Omega_Ref = Velocity_Ref/Length_Ref; config->SetOmega_Ref(Omega_Ref); - Force_Ref = config->GetDensity_Ref()*Velocity_Ref*Velocity_Ref*Length_Ref*Length_Ref; config->SetForce_Ref(Force_Ref); - Heat_Flux_Ref = Density_Ref*Velocity_Ref*Velocity_Ref*Velocity_Ref; config->SetHeat_Flux_Ref(Heat_Flux_Ref); - Gas_Constant_Ref = Velocity_Ref*Velocity_Ref/config->GetTemperature_Ref(); config->SetGas_Constant_Ref(Gas_Constant_Ref); - Viscosity_Ref = config->GetDensity_Ref()*Velocity_Ref*Length_Ref; config->SetViscosity_Ref(Viscosity_Ref); - Conductivity_Ref = Viscosity_Ref*Gas_Constant_Ref; config->SetConductivity_Ref(Conductivity_Ref); - Froude = ModVel_FreeStream/sqrt(STANDARD_GRAVITY*Length_Ref); config->SetFroude(Froude); - - /*--- Divide by reference values, to compute the non-dimensional free-stream values ---*/ + SetCompressibleNonDimValues(config, Pressure_FreeStream, Density_FreeStream, Temperature_FreeStream, + Viscosity_FreeStream, ModVel_FreeStream); - Pressure_FreeStreamND = Pressure_FreeStream/config->GetPressure_Ref(); config->SetPressure_FreeStreamND(Pressure_FreeStreamND); - Density_FreeStreamND = Density_FreeStream/config->GetDensity_Ref(); config->SetDensity_FreeStreamND(Density_FreeStreamND); + /*--- Reference value and non-dimensional free-stream value of the vib.-el. temperature ---*/ - for (auto iDim = 0ul; iDim < nDim; iDim++) { - Velocity_FreeStreamND[iDim] = config->GetVelocity_FreeStream()[iDim]/Velocity_Ref; config->SetVelocity_FreeStreamND(Velocity_FreeStreamND[iDim], iDim); - } - - Temperature_FreeStreamND = Temperature_FreeStream/config->GetTemperature_Ref(); config->SetTemperature_FreeStreamND(Temperature_FreeStreamND); - Temperature_ve_FreeStreamND = Temperature_ve_FreeStream/config->GetTemperature_ve_Ref(); config->SetTemperature_ve_FreeStreamND(Temperature_ve_FreeStreamND); - Gas_ConstantND = config->GetGas_Constant()/Gas_Constant_Ref; config->SetGas_ConstantND(Gas_ConstantND); - - ModVel_FreeStreamND = 0.0; - for (auto iDim = 0ul; iDim < nDim; iDim++) ModVel_FreeStreamND += Velocity_FreeStreamND[iDim]*Velocity_FreeStreamND[iDim]; - ModVel_FreeStreamND = sqrt(ModVel_FreeStreamND); config->SetModVel_FreeStreamND(ModVel_FreeStreamND); - - Viscosity_FreeStreamND = Viscosity_FreeStream / Viscosity_Ref; config->SetViscosity_FreeStreamND(Viscosity_FreeStreamND); - - Tke_FreeStream = 3.0/2.0*(ModVel_FreeStream*ModVel_FreeStream*config->GetTurbulenceIntensity_FreeStream()*config->GetTurbulenceIntensity_FreeStream()); - config->SetTke_FreeStream(Tke_FreeStream); + const su2double Temperature_ve_Ref = (config->GetRef_NonDim() == DIMENSIONAL) ? 1.0 : Temperature_ve_FreeStream; + config->SetTemperature_ve_Ref(Temperature_ve_Ref); - Tke_FreeStreamND = 3.0/2.0*(ModVel_FreeStreamND*ModVel_FreeStreamND*config->GetTurbulenceIntensity_FreeStream()*config->GetTurbulenceIntensity_FreeStream()); - config->SetTke_FreeStreamND(Tke_FreeStreamND); + const su2double Temperature_ve_FreeStreamND = Temperature_ve_FreeStream/config->GetTemperature_ve_Ref(); + config->SetTemperature_ve_FreeStreamND(Temperature_ve_FreeStreamND); - Omega_FreeStream = Density_FreeStream*Tke_FreeStream/max(Viscosity_FreeStream*config->GetTurb2LamViscRatio_FreeStream(), EPS); - config->SetOmega_FreeStream(Omega_FreeStream); + /*--- Retrieve the values that are needed below ---*/ - Omega_FreeStreamND = Density_FreeStreamND*Tke_FreeStreamND/max(Viscosity_FreeStreamND*config->GetTurb2LamViscRatio_FreeStream(), EPS); - config->SetOmega_FreeStreamND(Omega_FreeStreamND); + const su2double Time_Ref = config->GetTime_Ref(); + const su2double Froude = config->GetFroude(); + const su2double ModVel_FreeStreamND = config->GetModVel_FreeStreamND(); + const su2double Tke_FreeStreamND = config->GetTke_FreeStreamND(); /*--- Initialize the dimensionless Fluid Model that will be used to solve the dimensionless problem ---*/ From d60ab14440ffe606e2b0699b831d145c954f5f78 Mon Sep 17 00:00:00 2001 From: rois1995 Date: Mon, 12 Oct 2026 13:58:04 +0200 Subject: [PATCH 8/9] Add Andrea Rausa to AUTHORS.md Co-Authored-By: Claude Opus 5.5 Claude-Session: https://claude.ai/code/session_01L7oxzVir4zENrpy3K3W6Fm --- AUTHORS.md | 1 + 1 file changed, 1 insertion(+) diff --git a/AUTHORS.md b/AUTHORS.md index d82ccb599fe9..f45e08125967 100644 --- a/AUTHORS.md +++ b/AUTHORS.md @@ -50,6 +50,7 @@ Alessandro Gastaldi Aman uz zaman Baig Amit Sachdeva Ana Lourenco +Andrea Rausa Andrew Burkett Andrew Wendorff Angelo Passariello From aac18d620e6404157389f118302d50ff1a4392ca Mon Sep 17 00:00:00 2001 From: rois1995 Date: Mon, 12 Oct 2026 14:22:24 +0200 Subject: [PATCH 9/9] Style clean-up of the new NEMO helper functions Use GeometryToolbox::Norm for the non-dimensional velocity magnitude, wrap two lines longer than 120 characters, shorten three comments, and drop the val_ prefix from the arguments of the two wall helpers so that all the new helpers use the same naming. No functional change. Co-Authored-By: Claude Opus 5.5 Claude-Session: https://claude.ai/code/session_01L7oxzVir4zENrpy3K3W6Fm --- .../include/solvers/CFVMFlowSolverBase.inl | 4 +--- SU2_CFD/include/solvers/CNEMOEulerSolver.hpp | 2 +- SU2_CFD/include/solvers/CNEMONSSolver.hpp | 24 +++++++++---------- SU2_CFD/src/solvers/CEulerSolver.cpp | 5 ++-- SU2_CFD/src/solvers/CNEMOEulerSolver.cpp | 2 +- SU2_CFD/src/solvers/CNEMONSSolver.cpp | 16 ++++++------- 6 files changed, 26 insertions(+), 27 deletions(-) diff --git a/SU2_CFD/include/solvers/CFVMFlowSolverBase.inl b/SU2_CFD/include/solvers/CFVMFlowSolverBase.inl index e7f38e23e99b..865eb0d57d2c 100644 --- a/SU2_CFD/include/solvers/CFVMFlowSolverBase.inl +++ b/SU2_CFD/include/solvers/CFVMFlowSolverBase.inl @@ -931,9 +931,7 @@ void CFVMFlowSolverBase::SetCompressibleNonDimValues(CConfig *config, su2d const su2double Gas_ConstantND = config->GetGas_Constant()/Gas_Constant_Ref; config->SetGas_ConstantND(Gas_ConstantND); - su2double ModVel_FreeStreamND = 0.0; - for (auto iDim = 0u; iDim < nDim; iDim++) ModVel_FreeStreamND += Velocity_FreeStreamND[iDim]*Velocity_FreeStreamND[iDim]; - ModVel_FreeStreamND = sqrt(ModVel_FreeStreamND); + const su2double ModVel_FreeStreamND = GeometryToolbox::Norm(nDim, Velocity_FreeStreamND); config->SetModVel_FreeStreamND(ModVel_FreeStreamND); const su2double Viscosity_FreeStreamND = Viscosity_FreeStream / Viscosity_Ref; diff --git a/SU2_CFD/include/solvers/CNEMOEulerSolver.hpp b/SU2_CFD/include/solvers/CNEMOEulerSolver.hpp index 1c4e0a97eb67..cdd8acac0a4a 100644 --- a/SU2_CFD/include/solvers/CNEMOEulerSolver.hpp +++ b/SU2_CFD/include/solvers/CNEMOEulerSolver.hpp @@ -62,7 +62,7 @@ class CNEMOEulerSolver : public CFVMFlowSolverBaseSetReThetaT_FreeStream(Re_ThetaT_FreeStream); - const su2double MassDiffusivityND = config->GetDiffusivity_Constant() / (config->GetVelocity_Ref() * config->GetLength_Ref()); + const su2double MassDiffusivityND = + config->GetDiffusivity_Constant() / (config->GetVelocity_Ref() * config->GetLength_Ref()); config->SetDiffusivity_ConstantND(MassDiffusivityND); - /*--- Retrieve the values that are needed below ---*/ + /*--- Local copies of values stored in the config ---*/ const su2double Time_Ref = config->GetTime_Ref(); const su2double Froude = config->GetFroude(); diff --git a/SU2_CFD/src/solvers/CNEMOEulerSolver.cpp b/SU2_CFD/src/solvers/CNEMOEulerSolver.cpp index 41b7d980e93f..449ea5cf46bb 100644 --- a/SU2_CFD/src/solvers/CNEMOEulerSolver.cpp +++ b/SU2_CFD/src/solvers/CNEMOEulerSolver.cpp @@ -1126,7 +1126,7 @@ void CNEMOEulerSolver::SetNondimensionalization(CConfig *config, unsigned short const su2double Temperature_ve_FreeStreamND = Temperature_ve_FreeStream/config->GetTemperature_ve_Ref(); config->SetTemperature_ve_FreeStreamND(Temperature_ve_FreeStreamND); - /*--- Retrieve the values that are needed below ---*/ + /*--- Local copies of values stored in the config ---*/ const su2double Time_Ref = config->GetTime_Ref(); const su2double Froude = config->GetFroude(); diff --git a/SU2_CFD/src/solvers/CNEMONSSolver.cpp b/SU2_CFD/src/solvers/CNEMONSSolver.cpp index 3a5bc0431a6d..a229d707b648 100644 --- a/SU2_CFD/src/solvers/CNEMONSSolver.cpp +++ b/SU2_CFD/src/solvers/CNEMONSSolver.cpp @@ -221,10 +221,10 @@ void CNEMONSSolver::Viscous_Residual(CGeometry *geometry, } //iEdge } -void CNEMONSSolver::SetWallVelocity_Strong(unsigned long iPoint, const su2double *val_velocity) { +void CNEMONSSolver::SetWallVelocity_Strong(unsigned long iPoint, const su2double *velocity) { /*--- Store the wall velocity in the old solution, it is not changed by the update. ---*/ - nodes->SetVelocity_Old(iPoint, val_velocity); + nodes->SetVelocity_Old(iPoint, velocity); /*--- Remove any contribution to the momentum residual and truncation error at this node. ---*/ for (auto iDim = 0ul; iDim < nDim; iDim++) @@ -232,9 +232,9 @@ void CNEMONSSolver::SetWallVelocity_Strong(unsigned long iPoint, const su2double nodes->SetVel_ResTruncError_Zero(iPoint); } -void CNEMONSSolver::SetWallEnergyFlux_Weak(unsigned long iPoint, unsigned long jPoint, su2double val_Twall, - su2double val_Tvewall, su2double val_pcontrol, su2double val_area, - su2double val_dist, su2double *val_residual) const { +void CNEMONSSolver::SetWallEnergyFlux_Weak(unsigned long iPoint, unsigned long jPoint, su2double Twall, + su2double Tvewall, su2double pcontrol, su2double area, + su2double dist, su2double *residual) const { /*--- Temperatures at the wall point (i) and at its closest normal neighbor (j) ---*/ const su2double Ti = nodes->GetTemperature(iPoint); @@ -249,9 +249,9 @@ void CNEMONSSolver::SetWallEnergyFlux_Weak(unsigned long iPoint, unsigned long j /*--- Conductive heat flux from the normal neighbor, plus a term that is proportional to the difference between the target wall temperature and the current one. The second term drives the solution at the wall towards the target temperature. ---*/ - val_residual[nSpecies+nDim] = ((ktr*(Ti-Tj) + kve*(Tvei-Tvej)) + - (ktr*(val_Twall-Ti) + kve*(val_Tvewall-Tvei))*val_pcontrol)*val_area/val_dist; - val_residual[nSpecies+nDim+1] = (kve*(Tvei-Tvej) + kve*(val_Tvewall-Tvei)*val_pcontrol)*val_area/val_dist; + residual[nSpecies+nDim] = ((ktr*(Ti-Tj) + kve*(Tvei-Tvej)) + + (ktr*(Twall-Ti) + kve*(Tvewall-Tvei))*pcontrol)*area/dist; + residual[nSpecies+nDim+1] = (kve*(Tvei-Tvej) + kve*(Tvewall-Tvei)*pcontrol)*area/dist; } void CNEMONSSolver::BC_HeatFluxNonCatalytic_Wall(CGeometry *geometry,