Skip to content

Addition of Source Terms for NEMO Axisymmetric Flows #1162

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
jtneedels merged 9 commits into from
Jan 14, 2021
Merged

Addition of Source Terms for NEMO Axisymmetric Flows #1162

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
9 commits
Select commit Hold shift + click to select a range
2f91c39
Implemented additional source residual terms for NEMO axisymmetric si…
jtneedels Jan 12, 2021
37ac145
Eliminate unused variables, move initalization of temperature variabl…
jtneedels Jan 12, 2021
2b67845
Fixed issue with test failue due to missing chemistry source residual…
jtneedels Jan 13, 2021
3e4abab
Addressing PR comments (commenting)
jtneedels Jan 13, 2021
0c6c2ad
Merge branch 'develop' into feature_nemo_axi_viscous
pcarruscag Jan 13, 2021
cec4c4b
Addressing PR comments (formatting)
jtneedels Jan 13, 2021
854cde6
Merge branch 'feature_nemo_axi_viscous' of https://github.com/su2code…
jtneedels Jan 13, 2021
fe4df82
Addressing PR comments, formatting
jtneedels Jan 14, 2021
925bd40
Merge branch 'develop' into feature_nemo_axi_viscous
jtneedels Jan 14, 2021
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
58 changes: 57 additions & 1 deletion SU2_CFD/src/numerics/NEMO/NEMO_sources.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -296,7 +296,25 @@ CNumerics::ResidualType<> CSource_NEMO::ComputeVibRelaxation(const CConfig *conf
CNumerics::ResidualType<> CSource_NEMO::ComputeAxisymmetric(const CConfig *config) {

unsigned short iDim, iSpecies, iVar;
su2double rho, rhov, vel2, H, yinv;
su2double rho, rhov, vel2, H, yinv, T, Tve, qy_ve, Ru, RuSI;
su2double *Ds, **GV, ktr, kve;

/*--- Rename for convenience ---*/
Ds = Diffusion_Coeff_i;
ktr = Thermal_Conductivity_i;
kve = Thermal_Conductivity_ve_i;
rho = V_i[RHO_INDEX];
T = V_i[T_INDEX];
Tve = V_i[TVE_INDEX];
GV = PrimVar_Grad_i;
RuSI= UNIVERSAL_GAS_CONSTANT;
Ru = 1000.0*RuSI;

auto& Ms = fluidmodel->GetSpeciesMolarMass();

bool viscous = config->GetViscous();
bool rans = (config->GetKind_Turb_Model() != NONE);
hs = fluidmodel->ComputeSpeciesEnthalpy(T, Tve, eve_i);

/*--- Initialize residual and Jacobian arrays ---*/
for (iVar = 0; iVar < nVar; iVar++) {
Expand All @@ -312,6 +330,7 @@ CNumerics::ResidualType<> CSource_NEMO::ComputeAxisymmetric(const CConfig *confi
rhov = U_i[nSpecies+1];
H = V_i[H_INDEX];
vel2 = 0.0;

for (iDim = 0; iDim < nDim; iDim++)
vel2 += V_i[VEL_INDEX+iDim]*V_i[VEL_INDEX+iDim];
for (iSpecies = 0; iSpecies < nSpecies; iSpecies++)
Expand All @@ -324,6 +343,43 @@ CNumerics::ResidualType<> CSource_NEMO::ComputeAxisymmetric(const CConfig *confi
residual[nSpecies+2] = yinv*rhov*H*Volume;
residual[nSpecies+3] = yinv*rhov*U_i[nSpecies+nDim+1]/rho*Volume;

if (viscous) {
if (!rans){ turb_ke_i = 0.0; }

su2double sumJhs_y = 0.0;
su2double sumJeve_y = 0.0;
su2double Mass = 0.0;

for (iSpecies=0; iSpecies<nSpecies; iSpecies++)
Mass += V_i[iSpecies]*Ms[iSpecies];

su2double heat_capacity_cp_i = V_i[RHOCVTR_INDEX]+Ru/Mass;
su2double total_viscosity_i = Laminar_Viscosity_i + Eddy_Viscosity_i;
su2double total_conductivity_i = ktr + kve + heat_capacity_cp_i*Eddy_Viscosity_i/Prandtl_Turb;
su2double u = V_i[VEL_INDEX];
su2double v = V_i[VEL_INDEX+1];
su2double qy_ve = kve*GV[TVE_INDEX][1];

/*--- Enthalpy and vib-el energy transport due to y-direction diffusion---*/
for (iSpecies = 0; iSpecies < nHeavy; iSpecies++) {
sumJhs_y += (rho*Ds[iSpecies]*GV[RHOS_INDEX+iSpecies][1] - V_i[RHOS_INDEX+iSpecies]*Vector[1]) * hs[iSpecies];
sumJeve_y += (rho*Ds[iSpecies]*GV[RHOS_INDEX+iSpecies][1] - V_i[RHOS_INDEX+iSpecies]*Vector[1]) * eve_i[iSpecies];
}
Copy link
Contributor

@CatarinaGarbacz CatarinaGarbacz Jan 13, 2021

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

you can just ident this 2 lines for aesthetics (the += part)

Copy link
Member

@pcarruscag pcarruscag Jan 13, 2021

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I would not spend too much time with that kind of alignment, eventually we will run clang-format on the entire code... The only thing you can do to make code look good is to write small functions.


for (iSpecies = 0; iSpecies < nSpecies; iSpecies++)
residual[iSpecies] -= 0.0;
residual[nSpecies] -= Volume*(yinv*total_viscosity_i*(PrimVar_Grad_i[nSpecies+2][1]+PrimVar_Grad_i[nSpecies+3][0])
-TWO3*AuxVar_Grad_i[0][0]);
residual[nSpecies+1] -= Volume*(yinv*total_viscosity_i*2*(PrimVar_Grad_i[nSpecies+3][1]-v*yinv)
-TWO3*AuxVar_Grad_i[0][1]);
residual[nSpecies+2] -= Volume*(yinv*(- sumJhs_y + total_viscosity_i*(u*(PrimVar_Grad_i[nSpecies+3][0]+PrimVar_Grad_i[nSpecies+2][1])
+v*TWO3*(2*PrimVar_Grad_i[nSpecies+2][1]-PrimVar_Grad_i[nSpecies+2][0]
-v*yinv+rho*turb_ke_i))
-total_conductivity_i*PrimVar_Grad_i[nSpecies][1])
-TWO3*(AuxVar_Grad_i[1][1]+AuxVar_Grad_i[2][1]));
residual[nSpecies+3] -= Volume*(yinv*(-sumJeve_y -qy_ve));
}

// if (implicit) {
//
// /*--- Initialize ---*/
Expand Down
124 changes: 102 additions & 22 deletions SU2_CFD/src/solvers/CNEMOEulerSolver.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -1108,16 +1108,19 @@ bool CNEMOEulerSolver::CheckNonPhys(const su2double *V) const {

void CNEMOEulerSolver::Source_Residual(CGeometry *geometry, CSolver **solver_container, CNumerics **numerics_container, CConfig *config, unsigned short iMesh) {

unsigned short iVar;
unsigned short iVar, jVar;
unsigned long iPoint;
unsigned long eAxi_local, eChm_local, eVib_local;
unsigned long eAxi_global, eChm_global, eVib_global;

/*--- Assign booleans ---*/
bool err = false;
//bool implicit = (config->GetKind_TimeIntScheme_Flow() == EULER_IMPLICIT);
bool implicit = (config->GetKind_TimeIntScheme_Flow() == EULER_IMPLICIT);
bool frozen = config->GetFrozen();
bool monoatomic = config->GetMonoatomic();
bool viscous = config->GetViscous();
bool ideal_gas = (config->GetKind_FluidModel() == STANDARD_AIR) || (config->GetKind_FluidModel() == IDEAL_GAS);
bool rans = (config->GetKind_Turb_Model() != NONE);

CNumerics* numerics = numerics_container[SOURCE_FIRST_TERM];

Expand All @@ -1137,11 +1140,11 @@ void CNEMOEulerSolver::Source_Residual(CGeometry *geometry, CSolver **solver_con
numerics->SetPrimitive (nodes->GetPrimitive(iPoint), nodes->GetPrimitive(iPoint) );

/*--- Pass supplementary information to CNumerics ---*/
numerics->SetdPdU (nodes->GetdPdU(iPoint), nodes->GetdPdU(iPoint));
numerics->SetdTdU (nodes->GetdTdU(iPoint), nodes->GetdTdU(iPoint));
numerics->SetdPdU(nodes->GetdPdU(iPoint), nodes->GetdPdU(iPoint));
numerics->SetdTdU(nodes->GetdTdU(iPoint), nodes->GetdTdU(iPoint));
numerics->SetdTvedU(nodes->GetdTvedU(iPoint), nodes->GetdTvedU(iPoint));
numerics->SetEve (nodes->GetEve(iPoint), nodes->GetEve(iPoint));
numerics->SetCvve (nodes->GetCvve(iPoint), nodes->GetCvve(iPoint));
numerics->SetEve(nodes->GetEve(iPoint), nodes->GetEve(iPoint));
numerics->SetCvve(nodes->GetCvve(iPoint), nodes->GetCvve(iPoint));

/*--- Set volume of the dual grid cell ---*/
numerics->SetVolume(geometry->nodes->GetVolume(iPoint));
Expand All @@ -1150,25 +1153,102 @@ void CNEMOEulerSolver::Source_Residual(CGeometry *geometry, CSolver **solver_con

/*--- Compute axisymmetric source terms (if needed) ---*/
if (config->GetAxisymmetric()) {
auto residual = numerics->ComputeAxisymmetric(config);

/*--- Check for errors before applying source to the linear system ---*/
err = false;
for (iVar = 0; iVar < nVar; iVar++)
if (residual[iVar] != residual[iVar]) err = true;
//if (implicit)
// for (iVar = 0; iVar < nVar; iVar++)
// for (jVar = 0; jVar < nVar; jVar++)
// if (Jacobian_i[iVar][jVar] != Jacobian_i[iVar][jVar]) err = true;
if (viscous) {

/*--- Apply the update to the linear system ---*/
if (!err) {
LinSysRes.AddBlock(iPoint, residual);
//if (implicit)
// Jacobian.AddBlock(iPoint, iPoint, Jacobian_i);
for (iPoint = 0; iPoint < nPoint; iPoint++) {

su2double yCoord = geometry->nodes->GetCoord(iPoint, 1);
su2double yVelocity = nodes->GetVelocity(iPoint,1);
su2double xVelocity = nodes->GetVelocity(iPoint,0);
su2double Total_Viscosity = nodes->GetLaminarViscosity(iPoint) + nodes->GetEddyViscosity(iPoint);

if (yCoord > EPS){
su2double nu_v_on_y = Total_Viscosity*yVelocity/yCoord;
nodes->SetAuxVar(iPoint, 0, nu_v_on_y);
nodes->SetAuxVar(iPoint, 1, nu_v_on_y*yVelocity);
nodes->SetAuxVar(iPoint, 2, nu_v_on_y*xVelocity);
}
}

/*--- Compute the auxiliary variable gradient with GG or WLS. ---*/
if (config->GetKind_Gradient_Method() == GREEN_GAUSS) {
SetAuxVar_Gradient_GG(geometry, config);
}
if (config->GetKind_Gradient_Method() == WEIGHTED_LEAST_SQUARES) {
SetAuxVar_Gradient_LS(geometry, config);
}
}

/*--- loop over points ---*/
SU2_OMP_FOR_DYN(omp_chunk_size)
for (iPoint = 0; iPoint < nPointDomain; iPoint++) {

/*--- Set solution ---*/
numerics->SetConservative(nodes->GetSolution(iPoint), nodes->GetSolution(iPoint));

/*--- Set control volume ---*/
numerics->SetVolume(geometry->nodes->GetVolume(iPoint));

/*--- Set y coordinate ---*/
numerics->SetCoord(geometry->nodes->GetCoord(iPoint), geometry->nodes->GetCoord(iPoint));

/*--- Set primitive variables for viscous terms and/or generalised source ---*/
numerics->SetPrimitive(nodes->GetPrimitive(iPoint), nodes->GetPrimitive(iPoint));

/*--- If necessary, set variables needed for viscous computation ---*/
if (viscous) {

/*--- Set gradient of primitive variables ---*/
numerics->SetPrimVarGradient(nodes->GetGradient_Primitive(iPoint), nodes->GetGradient_Primitive(iPoint));

/*--- Set gradient of auxillary variables ---*/
numerics->SetAuxVarGrad(nodes->GetAuxVarGradient(iPoint), nullptr);

/*--- Set diffusion coefficient ---*/
numerics->SetDiffusionCoeff(nodes->GetDiffusionCoeff(iPoint), nodes->GetDiffusionCoeff(iPoint));

/*--- Laminar viscosity ---*/
numerics->SetLaminarViscosity(nodes->GetLaminarViscosity(iPoint), nodes->GetLaminarViscosity(iPoint));

/*--- Eddy viscosity ---*/
numerics->SetEddyViscosity(nodes->GetEddyViscosity(iPoint), nodes->GetEddyViscosity(iPoint));

/*--- Thermal conductivity ---*/
numerics->SetThermalConductivity(nodes->GetThermalConductivity(iPoint), nodes->GetThermalConductivity(iPoint));

/*--- Vib-el. thermal conductivity ---*/
numerics->SetThermalConductivity_ve(nodes->GetThermalConductivity_ve(iPoint), nodes->GetThermalConductivity_ve(iPoint));

/*--- Vib-el energy ---*/
numerics->SetEve(nodes->GetEve(iPoint), nodes->GetEve(iPoint));

/*--- Set turbulence kinetic energy ---*/
if (rans){
CVariable* turbNodes = solver_container[TURB_SOL]->GetNodes();
numerics->SetTurbKineticEnergy(turbNodes->GetSolution(iPoint,0), turbNodes->GetSolution(iPoint,0));
}
}

auto residual = numerics->ComputeAxisymmetric(config);

/*--- Check for errors before applying source to the linear system ---*/
err = false;
for (iVar = 0; iVar < nVar; iVar++)
if (residual[iVar] != residual[iVar]) err = true;
if (implicit)
for (iVar = 0; iVar < nVar; iVar++)
for (jVar = 0; jVar < nVar; jVar++)
if (Jacobian_i[iVar][jVar] != Jacobian_i[iVar][jVar]) err = true;

/*--- Apply the update to the linear system ---*/
if (!err) {
LinSysRes.AddBlock(iPoint, residual);
if (implicit)
Jacobian.AddBlock(iPoint, iPoint, Jacobian_i);
}else
eAxi_local++;
}
else
eAxi_local++;
}
Copy link
Contributor

@CatarinaGarbacz CatarinaGarbacz Jan 13, 2021

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

extra line between between the brackets

Copy link
Contributor Author

@jtneedels jtneedels Jan 14, 2021 •
edited
Loading

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Thanks Catarina, I think I've addressed this and the other comments on formatting.


if(!monoatomic){
Expand Down
7 changes: 7 additions & 0 deletions SU2_CFD/src/variables/CNEMOEulerVariable.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -78,6 +78,13 @@ CNEMOEulerVariable::CNEMOEulerVariable(su2double val_pressure,

Res_TruncError.resize(nPoint,nVar) = su2double(0.0);

/*--- Size Grad_AuxVar for axiysmmetric ---*/
if (config->GetAxisymmetric()){
nAuxVar = 3;
Grad_AuxVar.resize(nPoint,nAuxVar,nDim,0.0);
AuxVar.resize(nPoint,nAuxVar) = su2double(0.0);
}

/*--- Only for residual smoothing (multigrid) ---*/

for (unsigned long iMesh = 0; iMesh <= config->GetnMGLevels(); iMesh++) {
Expand Down