Skip to content

Commit

Permalink
reformulating and cleaning
Browse files Browse the repository at this point in the history
  • Loading branch information
@Cristopher-Morales
Cristopher-Morales committed Feb 18, 2025
1 parent 5c656bb commit 3f45e1c
Show file tree
Hide file tree
Showing 5 changed files with 63 additions and 119 deletions.
3 changes: 3 additions & 0 deletions SU2_CFD/include/numerics/CNumerics.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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. */
Expand Down
2 changes: 1 addition & 1 deletion SU2_CFD/include/numerics/flow/convection/fds.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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. */
Expand Down
53 changes: 13 additions & 40 deletions SU2_CFD/src/numerics/CNumerics.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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 {

Expand All @@ -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);
}
Expand All @@ -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;
Expand All @@ -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));

}

Expand Down
62 changes: 24 additions & 38 deletions SU2_CFD/src/numerics/flow/convection/fds.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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;
Expand All @@ -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);

Expand All @@ -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 ---*/
Expand Down Expand Up @@ -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
Expand All @@ -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,
Expand Down Expand Up @@ -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;
}
}
}
Expand Down
62 changes: 22 additions & 40 deletions SU2_CFD/src/solvers/CIncEulerSolver.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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);
Expand All @@ -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;
}
Expand All @@ -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++) {
Expand All @@ -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 ++ )
Expand All @@ -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;
}
}
}
Expand Down

0 comments on commit 3f45e1c

Please sign in to comment.