From 3f45e1c969f031c19f3558b005c85eb8125b3b7d Mon Sep 17 00:00:00 2001 From: Cristopher-Morales Date: Tue, 18 Feb 2025 12:30:20 +0100 Subject: [PATCH] reformulating and cleaning --- SU2_CFD/include/numerics/CNumerics.hpp | 3 + .../include/numerics/flow/convection/fds.hpp | 2 +- SU2_CFD/src/numerics/CNumerics.cpp | 53 ++++------------ SU2_CFD/src/numerics/flow/convection/fds.cpp | 62 +++++++------------ SU2_CFD/src/solvers/CIncEulerSolver.cpp | 62 +++++++------------ 5 files changed, 63 insertions(+), 119 deletions(-) diff --git a/SU2_CFD/include/numerics/CNumerics.hpp b/SU2_CFD/include/numerics/CNumerics.hpp index cec666bc72d..37ed3a036fc 100644 --- a/SU2_CFD/include/numerics/CNumerics.hpp +++ b/SU2_CFD/include/numerics/CNumerics.hpp @@ -103,6 +103,9 @@ class CNumerics { su2double Enthalpy_i, /*!< \brief Enthalpy at point i. */ Enthalpy_j; /*!< \brief Enthalpy at point j. */ + su2double + WorkingVariable_i, /*!< \brief Enthalpy at point i. */ + WorkingVariable_j; /*!< \brief Enthalpy at point j. */ su2double dist_i, /*!< \brief Distance of point i to the nearest wall. */ dist_j; /*!< \brief Distance of point j to the nearest wall. */ diff --git a/SU2_CFD/include/numerics/flow/convection/fds.hpp b/SU2_CFD/include/numerics/flow/convection/fds.hpp index e72db06db6a..2e39d8fc10b 100644 --- a/SU2_CFD/include/numerics/flow/convection/fds.hpp +++ b/SU2_CFD/include/numerics/flow/convection/fds.hpp @@ -51,7 +51,7 @@ class CUpwFDSInc_Flow final : public CNumerics { Pressure_j, ProjVelocity, MeandRhodT, dRhodT_i, dRhodT_j, /*!< \brief Derivative of density w.r.t. temperature (variable density flows). */ Temperature_i, Temperature_j, /*!< \brief Temperature at node 0 and 1. */ - MeanDensity, MeanPressure, MeanSoundSpeed, MeanBetaInc2, MeanEnthalpy, MeanCp, MeanTemperature; /*!< \brief Mean values of primitive variables. */ + MeanDensity, MeanPressure, MeanSoundSpeed, MeanBetaInc2, MeanWorkingVar, MeanCp, MeanTemperature; /*!< \brief Mean values of primitive variables. */ unsigned short iDim, iVar, jVar, kVar; su2double* Flux = nullptr; /*!< \brief The flux / residual across the edge. */ diff --git a/SU2_CFD/src/numerics/CNumerics.cpp b/SU2_CFD/src/numerics/CNumerics.cpp index f78e4b60288..116ce359ee5 100644 --- a/SU2_CFD/src/numerics/CNumerics.cpp +++ b/SU2_CFD/src/numerics/CNumerics.cpp @@ -279,20 +279,10 @@ void CNumerics::GetInviscidIncProjJac(const su2double *val_density, const su2dou val_Proj_Jac_Tensor[2][2] = val_scale*((*val_density)*(proj_vel + val_normal[1]*val_velocity[1])); val_Proj_Jac_Tensor[2][3] = val_scale*((*val_dRhodT)*val_velocity[1]*proj_vel); - if (energy_multicomponent) { - val_Proj_Jac_Tensor[3][0] = val_scale * ((*val_temperature) * proj_vel / (*val_betainc2)); - val_Proj_Jac_Tensor[3][1] = val_scale * ((*val_temperature) * val_normal[0] * (*val_density)); - val_Proj_Jac_Tensor[3][2] = val_scale * ((*val_temperature) * val_normal[1] * (*val_density)); - val_Proj_Jac_Tensor[3][3] = - val_scale * ((*val_temperature) * (*val_dRhodT) + (*val_density)) * proj_vel; - - } else { - val_Proj_Jac_Tensor[3][0] = val_scale * ((*val_cp) * (*val_temperature) * proj_vel / (*val_betainc2)); - val_Proj_Jac_Tensor[3][1] = val_scale * ((*val_cp) * (*val_temperature) * val_normal[0] * (*val_density)); - val_Proj_Jac_Tensor[3][2] = val_scale * ((*val_cp) * (*val_temperature) * val_normal[1] * (*val_density)); - val_Proj_Jac_Tensor[3][3] = - val_scale * ((*val_cp) * ((*val_temperature) * (*val_dRhodT) + (*val_density)) * proj_vel); - } + val_Proj_Jac_Tensor[3][0] = val_scale*((*val_cp)*(*val_temperature)*proj_vel/(*val_betainc2)); + val_Proj_Jac_Tensor[3][1] = val_scale*((*val_cp)*(*val_temperature)*val_normal[0]*(*val_density)); + val_Proj_Jac_Tensor[3][2] = val_scale*((*val_cp)*(*val_temperature)*val_normal[1]*(*val_density)); + val_Proj_Jac_Tensor[3][3] = val_scale*((*val_cp)*((*val_temperature)*(*val_dRhodT) + (*val_density))*proj_vel); } else { @@ -319,22 +309,13 @@ void CNumerics::GetInviscidIncProjJac(const su2double *val_density, const su2dou val_Proj_Jac_Tensor[3][2] = val_scale*(val_normal[1]*(*val_density)*val_velocity[2]); val_Proj_Jac_Tensor[3][3] = val_scale*((*val_density)*(proj_vel + val_normal[2]*val_velocity[2])); val_Proj_Jac_Tensor[3][4] = val_scale*((*val_dRhodT)*val_velocity[2]*proj_vel); - - if (energy_multicomponent) { - val_Proj_Jac_Tensor[4][0] = val_scale * ((*val_temperature) * proj_vel / (*val_betainc2)); - val_Proj_Jac_Tensor[4][1] = val_scale * ((*val_temperature) * val_normal[0] * (*val_density)); - val_Proj_Jac_Tensor[4][2] = val_scale * ((*val_temperature) * val_normal[1] * (*val_density)); - val_Proj_Jac_Tensor[4][3] = val_scale * ((*val_temperature) * val_normal[2] * (*val_density)); - val_Proj_Jac_Tensor[4][4] = - val_scale * ((*val_temperature) * (*val_dRhodT) + (*val_density)) * proj_vel; - } else { - val_Proj_Jac_Tensor[4][0] = val_scale * ((*val_cp) * (*val_temperature) * proj_vel / (*val_betainc2)); - val_Proj_Jac_Tensor[4][1] = val_scale * ((*val_cp) * (*val_temperature) * val_normal[0] * (*val_density)); - val_Proj_Jac_Tensor[4][2] = val_scale * ((*val_cp) * (*val_temperature) * val_normal[1] * (*val_density)); - val_Proj_Jac_Tensor[4][3] = val_scale * ((*val_cp) * (*val_temperature) * val_normal[2] * (*val_density)); - val_Proj_Jac_Tensor[4][4] = - val_scale * ((*val_cp) * ((*val_temperature) * (*val_dRhodT) + (*val_density)) * proj_vel); - } + + val_Proj_Jac_Tensor[4][0] = val_scale*((*val_cp)*(*val_temperature)*proj_vel/(*val_betainc2)); + val_Proj_Jac_Tensor[4][1] = val_scale*((*val_cp)*(*val_temperature)*val_normal[0]*(*val_density)); + val_Proj_Jac_Tensor[4][2] = val_scale*((*val_cp)*(*val_temperature)*val_normal[1]*(*val_density)); + val_Proj_Jac_Tensor[4][3] = val_scale*((*val_cp)*(*val_temperature)*val_normal[2]*(*val_density)); + val_Proj_Jac_Tensor[4][4] = val_scale*((*val_cp)*((*val_temperature)*(*val_dRhodT) + (*val_density))*proj_vel); + } AD::EndPassive(wasActive); } @@ -348,11 +329,7 @@ void CNumerics::GetPreconditioner(const su2double *val_density, const su2double val_Precon[0][0] = 1.0/(*val_betainc2); for (iDim = 0; iDim < nDim; iDim++) val_Precon[iDim+1][0] = val_velocity[iDim]/(*val_betainc2); - if (energy_multicomponent){ - val_Precon[nDim+1][0] = (*val_temperature)/(*val_betainc2); - }else{ - val_Precon[nDim+1][0] = (*val_cp)*(*val_temperature)/(*val_betainc2); - } + val_Precon[nDim+1][0] = (*val_cp)*(*val_temperature)/(*val_betainc2); for (jDim = 0; jDim < nDim; jDim++) { val_Precon[0][jDim+1] = 0.0; @@ -366,11 +343,7 @@ void CNumerics::GetPreconditioner(const su2double *val_density, const su2double val_Precon[0][nDim+1] = (*val_drhodt); for (iDim = 0; iDim < nDim; iDim++) val_Precon[iDim+1][nDim+1] = val_velocity[iDim]*(*val_drhodt); - if (energy_multicomponent){ - val_Precon[nDim+1][nDim+1] = (*val_drhodt)*(*val_temperature) + (*val_density); - }else{ - val_Precon[nDim+1][nDim+1] = (*val_cp)*((*val_drhodt)*(*val_temperature) + (*val_density)); - } + val_Precon[nDim+1][nDim+1] = (*val_cp)*((*val_drhodt)*(*val_temperature) + (*val_density)); } diff --git a/SU2_CFD/src/numerics/flow/convection/fds.cpp b/SU2_CFD/src/numerics/flow/convection/fds.cpp index 115a966896b..a941a8008e6 100644 --- a/SU2_CFD/src/numerics/flow/convection/fds.cpp +++ b/SU2_CFD/src/numerics/flow/convection/fds.cpp @@ -121,9 +121,13 @@ CNumerics::ResidualType<> CUpwFDSInc_Flow::ComputeResidual(const CConfig *config if (energy_multicomponent) { Enthalpy_i = V_i[nDim + 9]; Enthalpy_j = V_j[nDim + 9]; + WorkingVariable_i = Enthalpy_i; + WorkingVariable_j = Enthalpy_j; } else { Enthalpy_i = Cp_i * Temperature_i; Enthalpy_j = Cp_j * Temperature_j; + WorkingVariable_i = Temperature_i; + WorkingVariable_j = Temperature_j; } ProjVelocity = 0.0; @@ -149,7 +153,7 @@ CNumerics::ResidualType<> CUpwFDSInc_Flow::ComputeResidual(const CConfig *config MeanDensity = 0.5*(DensityInc_i + DensityInc_j); MeanPressure = 0.5*(Pressure_i + Pressure_j); MeanBetaInc2 = 0.5*(BetaInc2_i + BetaInc2_j); - MeanEnthalpy = 0.5*(Enthalpy_i + Enthalpy_j); + MeanWorkingVar = 0.5*(WorkingVariable_i + WorkingVariable_j); MeanCp = 0.5*(Cp_i + Cp_j); MeanTemperature = 0.5*(Temperature_i + Temperature_j); @@ -163,15 +167,17 @@ CNumerics::ResidualType<> CUpwFDSInc_Flow::ComputeResidual(const CConfig *config MeandRhodT = 0.0; dRhodT_i = 0.0; dRhodT_j = 0.0; if (variable_density) { - if (energy_multicomponent) { - MeandRhodT = -MeanDensity / (MeanCp * MeanTemperature); - dRhodT_i = -DensityInc_i / (Cp_i * Temperature_i); - dRhodT_j = -DensityInc_j / (Cp_j * Temperature_j); - } else { - MeandRhodT = -MeanDensity / MeanTemperature; - dRhodT_i = -DensityInc_i / Temperature_i; - dRhodT_j = -DensityInc_j / Temperature_j; - } + MeandRhodT = -MeanDensity / MeanTemperature; + dRhodT_i = -DensityInc_i / Temperature_i; + dRhodT_j = -DensityInc_j / Temperature_j; + } + if (energy_multicomponent) { + MeandRhodT /= MeanCp ; + dRhodT_i /= Cp_i; + dRhodT_j /= Cp_j; + MeanCp = 1.0; + Cp_i = 1.0; + Cp_j = 1.0; } /*--- Compute ProjFlux_i ---*/ @@ -204,11 +210,7 @@ CNumerics::ResidualType<> CUpwFDSInc_Flow::ComputeResidual(const CConfig *config Lambda[iVar] = fabs(Lambda[iVar]); /*--- Build the preconditioning matrix using mean values ---*/ - if (energy_multicomponent) { - GetPreconditioner(&MeanDensity, MeanVelocity, &MeanBetaInc2, &MeanCp, &MeanEnthalpy, &MeandRhodT, Precon); - } else { - GetPreconditioner(&MeanDensity, MeanVelocity, &MeanBetaInc2, &MeanCp, &MeanTemperature, &MeandRhodT, Precon); - } + GetPreconditioner(&MeanDensity, MeanVelocity, &MeanBetaInc2, &MeanCp, &MeanWorkingVar, &MeandRhodT, Precon); /*--- Build the absolute value of the preconditioned Jacobian, i.e., |A_precon| = P x |Lambda| x inv(P), where P diagonalizes the matrix @@ -221,26 +223,15 @@ CNumerics::ResidualType<> CUpwFDSInc_Flow::ComputeResidual(const CConfig *config Diff_V[0] = Pressure_j - Pressure_i; for (iDim = 0; iDim < nDim; iDim++) Diff_V[iDim+1] = Velocity_j[iDim] - Velocity_i[iDim]; - if (energy_multicomponent){ - Diff_V[nDim+1] = Enthalpy_j - Enthalpy_i; - }else{ - Diff_V[nDim+1] = Temperature_j - Temperature_i; - } + Diff_V[nDim + 1] = WorkingVariable_j - WorkingVariable_i; /*--- Build the inviscid Jacobian w.r.t. the primitive variables ---*/ if (implicit) { - if (energy_multicomponent) { - GetInviscidIncProjJac(&DensityInc_i, Velocity_i, &BetaInc2_i, &Cp_i, &Enthalpy_i, &dRhodT_i, Normal, 0.5, - Jacobian_i); - GetInviscidIncProjJac(&DensityInc_j, Velocity_j, &BetaInc2_j, &Cp_j, &Enthalpy_j, &dRhodT_j, Normal, 0.5, - Jacobian_j); - } else { - GetInviscidIncProjJac(&DensityInc_i, Velocity_i, &BetaInc2_i, &Cp_i, &Temperature_i, &dRhodT_i, Normal, 0.5, - Jacobian_i); - GetInviscidIncProjJac(&DensityInc_j, Velocity_j, &BetaInc2_j, &Cp_j, &Temperature_j, &dRhodT_j, Normal, 0.5, - Jacobian_j); - } + GetInviscidIncProjJac(&DensityInc_i, Velocity_i, &BetaInc2_i, &Cp_i, &WorkingVariable_i, &dRhodT_i, Normal, 0.5, + Jacobian_i); + GetInviscidIncProjJac(&DensityInc_j, Velocity_j, &BetaInc2_j, &Cp_j, &WorkingVariable_j, &dRhodT_j, Normal, 0.5, + Jacobian_j); } /*--- Compute dissipation as Precon x |A_precon| x dV. If implicit, @@ -286,13 +277,8 @@ CNumerics::ResidualType<> CUpwFDSInc_Flow::ComputeResidual(const CConfig *config Jacobian_i[iDim+1][iDim+1] -= 0.5*ProjVelocity*DensityInc_i; Jacobian_j[iDim+1][iDim+1] -= 0.5*ProjVelocity*DensityInc_j; } - if (energy_multicomponent){ - Jacobian_i[nDim+1][nDim+1] -= 0.5*ProjVelocity*DensityInc_i; - Jacobian_j[nDim+1][nDim+1] -= 0.5*ProjVelocity*DensityInc_j; - } else { - Jacobian_i[nDim+1][nDim+1] -= 0.5*ProjVelocity*DensityInc_i*Cp_i; - Jacobian_j[nDim+1][nDim+1] -= 0.5*ProjVelocity*DensityInc_j*Cp_j; - } + Jacobian_i[nDim+1][nDim+1] -= 0.5*ProjVelocity*DensityInc_i*Cp_i; + Jacobian_j[nDim+1][nDim+1] -= 0.5*ProjVelocity*DensityInc_j*Cp_j; } } } diff --git a/SU2_CFD/src/solvers/CIncEulerSolver.cpp b/SU2_CFD/src/solvers/CIncEulerSolver.cpp index aa6bea96a0d..a939af159e1 100644 --- a/SU2_CFD/src/solvers/CIncEulerSolver.cpp +++ b/SU2_CFD/src/solvers/CIncEulerSolver.cpp @@ -2073,7 +2073,10 @@ void CIncEulerSolver::SetPreconditioner(const CConfig *config, unsigned long iPo Cp = nodes->GetSpecificHeatCp(iPoint); oneOverCp = 1.0/Cp; Temperature = nodes->GetTemperature(iPoint); - if (energy && multicomponent) Enthalpy = nodes->GetEnthalpy(iPoint); + Enthalpy = Cp * Temperature; + if (energy && multicomponent) { + Enthalpy = nodes->GetEnthalpy(iPoint); + } for (iDim = 0; iDim < nDim; iDim++) Velocity[iDim] = nodes->GetVelocity(iPoint,iDim); @@ -2085,6 +2088,7 @@ void CIncEulerSolver::SetPreconditioner(const CConfig *config, unsigned long iPo if (variable_density) { if (multicomponent && energy){ dRhodT = -Density / (Cp * Temperature); + Cp = oneOverCp = 1.0; } else { dRhodT = -Density/Temperature; } @@ -2105,13 +2109,9 @@ void CIncEulerSolver::SetPreconditioner(const CConfig *config, unsigned long iPo Preconditioner[iDim+1][0] = Velocity[iDim]/BetaInc2; if (energy) { - if (multicomponent) { - Preconditioner[nDim + 1][0] = Enthalpy / BetaInc2; - } else { - Preconditioner[nDim + 1][0] = Cp * Temperature / BetaInc2; - } + Preconditioner[nDim+1][0] = Enthalpy / BetaInc2; } else { - Preconditioner[nDim + 1][0] = 0.0; + Preconditioner[nDim+1][0] = 0.0; } for (jDim = 0; jDim < nDim; jDim++) { @@ -2128,13 +2128,9 @@ void CIncEulerSolver::SetPreconditioner(const CConfig *config, unsigned long iPo Preconditioner[iDim+1][nDim+1] = Velocity[iDim]*dRhodT; if (energy) { - if (multicomponent) { - Preconditioner[nDim + 1][nDim + 1] = dRhodT * Enthalpy + Density; - } else { - Preconditioner[nDim + 1][nDim + 1] = Cp * (dRhodT * Temperature + Density); - } + Preconditioner[nDim+1][nDim+1] = dRhodT * Enthalpy + Cp * Density; } else { - Preconditioner[nDim + 1][nDim + 1] = 1.0; + Preconditioner[nDim+1][nDim+1] = 1.0; } for (iVar = 0; iVar < nVar; iVar ++ ) @@ -2147,50 +2143,36 @@ void CIncEulerSolver::SetPreconditioner(const CConfig *config, unsigned long iPo Therefore, we build inv(Precon) here and multiply by the residual later in the R-K and Euler Explicit time integration schemes. ---*/ - if (multicomponent && energy) { - Preconditioner[0][0] = Enthalpy * BetaInc2 * dRhodT / Density + BetaInc2; - } else { - Preconditioner[0][0] = Temperature * BetaInc2 * dRhodT / Density + BetaInc2; - } + + Preconditioner[0][0] = Enthalpy * BetaInc2 * dRhodT * oneOverCp / Density + BetaInc2; + for (iDim = 0; iDim < nDim; iDim++) Preconditioner[iDim + 1][0] = -1.0 * Velocity[iDim] / Density; if (energy) { - if (multicomponent) { - Preconditioner[nDim + 1][0] = -1.0 * Enthalpy / Density; - } else { - Preconditioner[nDim + 1][0] = -1.0 * Temperature / Density; - } + Preconditioner[nDim+1][0] = -1.0 * Enthalpy * oneOverCp / Density; } else { - Preconditioner[nDim + 1][0] = 0.0; + Preconditioner[nDim+1][0] = 0.0; } for (jDim = 0; jDim < nDim; jDim++) { - Preconditioner[0][jDim + 1] = 0.0; + Preconditioner[0][jDim+1] = 0.0; for (iDim = 0; iDim < nDim; iDim++) { if (iDim == jDim) - Preconditioner[iDim + 1][jDim + 1] = 1.0 / Density; + Preconditioner[iDim+1][jDim+1] = 1.0 / Density; else - Preconditioner[iDim + 1][jDim + 1] = 0.0; + Preconditioner[iDim+1][jDim+1] = 0.0; } - Preconditioner[nDim + 1][jDim + 1] = 0.0; + Preconditioner[nDim+1][jDim+1] = 0.0; } - if (multicomponent && energy) { - Preconditioner[0][nDim + 1] = -1.0 * BetaInc2 * dRhodT / Density; - } else { - Preconditioner[0][nDim + 1] = -1.0 * BetaInc2 * dRhodT * oneOverCp / Density; - } - for (iDim = 0; iDim < nDim; iDim++) Preconditioner[iDim + 1][nDim + 1] = 0.0; + Preconditioner[0][nDim+1] = -1.0 * BetaInc2 * dRhodT * oneOverCp / Density; + for (iDim = 0; iDim < nDim; iDim++) Preconditioner[iDim+1][nDim+1] = 0.0; if (energy) { - if (multicomponent) { - Preconditioner[nDim + 1][nDim + 1] = 1.0 / Density; - } else { - Preconditioner[nDim + 1][nDim + 1] = oneOverCp / Density; - } + Preconditioner[nDim+1][nDim+1] = oneOverCp / Density; } else { - Preconditioner[nDim + 1][nDim + 1] = 0.0; + Preconditioner[nDim+1][nDim+1] = 0.0; } } }