@@ -484,13 +484,14 @@ void CNEMOEulerSolver::Upwind_Residual(CGeometry *geometry, CSolver **solver_con
484484 CNumerics* numerics = numerics_container[CONV_TERM];
485485
486486 /*--- Static arrays for MUSCL reconstructed variables ---*/
487- su2double Primitive_i[MAXNVAR] = {0.0}, Primitive_j[MAXNVAR] = {0.0};
488- su2double Conserved_i[MAXNVAR] = {0.0}, Conserved_j[MAXNVAR] = {0.0};
489- su2double dPdU_i[MAXNVAR] = {0.0}, dPdU_j[MAXNVAR] = {0.0};
490- su2double dTdU_i[MAXNVAR] = {0.0}, dTdU_j[MAXNVAR] = {0.0};
491- su2double dTvedU_i[MAXNVAR] = {0.0}, dTvedU_j[MAXNVAR] = {0.0};
492- su2double Eve_i[MAXNVAR] = {0.0}, Eve_j[MAXNVAR] = {0.0};
493- su2double Cvve_i[MAXNVAR] = {0.0}, Cvve_j[MAXNVAR] = {0.0};
487+ su2double Primitive_i[MAXNVAR] = {0.0}, Primitive_j[MAXNVAR] = {0.0};
488+ su2double Conserved_i[MAXNVAR] = {0.0}, Conserved_j[MAXNVAR] = {0.0};
489+ su2double dPdU_i[MAXNVAR] = {0.0}, dPdU_j[MAXNVAR] = {0.0};
490+ su2double dTdU_i[MAXNVAR] = {0.0}, dTdU_j[MAXNVAR] = {0.0};
491+ su2double dTvedU_i[MAXNVAR] = {0.0}, dTvedU_j[MAXNVAR] = {0.0};
492+ su2double Eve_i[MAXNVAR] = {0.0}, Eve_j[MAXNVAR] = {0.0};
493+ su2double Cvve_i[MAXNVAR] = {0.0}, Cvve_j[MAXNVAR] = {0.0};
494+ su2double Project_Grad_i[MAXNVAR] = {0.0}, Project_Grad_j[MAXNVAR] = {0.0};
494495 su2double Gamma_i = 0.0, Gamma_j = 0.0;
495496
496497 /*--- Loop over edges and calculate convective fluxes ---*/
@@ -535,32 +536,45 @@ void CNEMOEulerSolver::Upwind_Residual(CGeometry *geometry, CSolver **solver_con
535536 auto Gradient_i = nodes->GetGradient_Reconstruction(iPoint);
536537 auto Gradient_j = nodes->GetGradient_Reconstruction(jPoint);
537538
538- for (iVar = 0; iVar < nPrimVarGrad; iVar++) {
539+ /*--- Set and extract limiters ---*/
540+ su2double *Limiter_i = nullptr, *Limiter_j = nullptr;
539541
540- su2double Project_Grad_i = 0.0;
541- su2double Project_Grad_j = 0.0;
542+ if (limiter && !van_albada){
543+ Limiter_i = nodes->GetLimiter_Primitive(iPoint);
544+ Limiter_j = nodes->GetLimiter_Primitive(jPoint);
545+ }
542546
543- for (iDim = 0; iDim < nDim; iDim++) {
544- Project_Grad_i += Vector_ij[iDim]*Gradient_i[iVar][iDim];
545- Project_Grad_j -= Vector_ij[iDim]*Gradient_j[iVar][iDim];
546- }
547+ su2double lim_i = 2.0;
548+ su2double lim_j = 2.0;
547549
548- su2double lim_i = 1.0;
549- su2double lim_j = 1 .0;
550+ for (iVar = 0; iVar < nPrimVarGrad; iVar++) {
551+ Project_Grad_i[iVar] = 0.0; Project_Grad_j[iVar] = 0 .0;
550552
551- if (van_albada) {
552- su2double V_ij = V_j[iVar] - V_i[iVar];
553- lim_i = LimiterHelpers<>::vanAlbadaFunction(Project_Grad_i, V_ij, EPS);
554- lim_j = LimiterHelpers<>::vanAlbadaFunction(-Project_Grad_j, V_ij, EPS);
553+ for (iDim = 0; iDim < nDim; iDim++) {
554+ Project_Grad_i[iVar] += Vector_ij[iDim]*Gradient_i[iVar][iDim];
555+ Project_Grad_j[iVar] -= Vector_ij[iDim]*Gradient_j[iVar][iDim];
555556 }
556- else if (limiter) {
557- lim_i = nodes->GetLimiter_Primitive(iPoint, iVar);
558- lim_j = nodes->GetLimiter_Primitive(jPoint, iVar);
557+
558+ if (limiter) {
559+ if (van_albada) {
560+ su2double V_ij = V_j[iVar] - V_i[iVar];
561+ su2double va_lim_i = LimiterHelpers<>::vanAlbadaFunction(Project_Grad_i[iVar], V_ij, EPS);
562+ su2double va_lim_j = LimiterHelpers<>::vanAlbadaFunction(-Project_Grad_j[iVar], V_ij, EPS);
563+ lim_i = min(lim_i, va_lim_i);
564+ lim_j = min(lim_j, va_lim_j);
565+ } else {
566+ lim_i = min(lim_i, Limiter_i[iVar]);
567+ lim_j = min(lim_j, Limiter_j[iVar]);
568+ }
569+ } else {
570+ lim_i = lim_j = 1.0;
559571 }
560- su2double lim_ij = min(lim_i, lim_j);
572+ }
573+ su2double lim_ij = min(lim_i, lim_j);
561574
562- Primitive_i[iVar] = V_i[iVar] + lim_ij*Project_Grad_i;
563- Primitive_j[iVar] = V_j[iVar] + lim_ij*Project_Grad_j;
575+ for (iVar = 0; iVar < nPrimVarGrad; iVar++) {
576+ Primitive_i[iVar] = V_i[iVar] + lim_ij*Project_Grad_i[iVar];
577+ Primitive_j[iVar] = V_j[iVar] + lim_ij*Project_Grad_j[iVar];
564578 }
565579
566580 /*--- Check for non-physical solutions after reconstruction. If found, use the
0 commit comments