@@ -461,6 +461,7 @@ void CScalarSolver::SetResidual_DualTime(CGeometry* geometry, CSolver** solver_c
461461 unsigned long iPoint, jPoint, iVertex, iEdge;
462462
463463 su2double *U_time_nM1 = nullptr , *U_time_n = nullptr , *U_time_nP1 = nullptr ;
464+ /* --- For non-Conservative scalars (e.g. SA), multiply the Primitives with 1.0 instead of Density. ---*/
464465 su2double Density_nM1 = 1.0 , Density_n = 1.0 , Density_nP1 = 1.0 ;
465466 su2double Volume_nM1, Volume_nP1;
466467 const su2double *Normal = nullptr , *GridVel_i = nullptr , *GridVel_j = nullptr ;
@@ -540,7 +541,6 @@ void CScalarSolver::SetResidual_DualTime(CGeometry* geometry, CSolver** solver_c
540541 for (iPoint = 0 ; iPoint < nPointDomain; ++iPoint) {
541542 GridVel_i = geometry->nodes ->GetGridVel (iPoint);
542543 U_time_n = nodes->GetSolution_time_n (iPoint);
543- Density_n = 1.0 ;
544544
545545 if (Conservative) {
546546 if (incompressible)
@@ -603,11 +603,10 @@ void CScalarSolver::SetResidual_DualTime(CGeometry* geometry, CSolver** solver_c
603603 Density_n = flowNodes->GetDensity (iPoint); // Temporary fix
604604 else
605605 Density_n = flowNodes->GetSolution_time_n (iPoint, 0 );
606-
607- for (iVar = 0 ; iVar < nVar; iVar++) LinSysRes (iPoint, iVar) += Density_n * U_time_n[iVar] * Residual_GCL;
608- } else {
609- for (iVar = 0 ; iVar < nVar; iVar++) LinSysRes (iPoint, iVar) += U_time_n[iVar] * Residual_GCL;
610606 }
607+
608+ for (iVar = 0 ; iVar < nVar; iVar++) LinSysRes (iPoint, iVar) += Density_n * U_time_n[iVar] * Residual_GCL;
609+
611610 }
612611 END_SU2_OMP_FOR
613612 }
@@ -655,24 +654,16 @@ void CScalarSolver::SetResidual_DualTime(CGeometry* geometry, CSolver** solver_c
655654 Density_n = flowNodes->GetSolution_time_n (iPoint, 0 );
656655 Density_nP1 = flowNodes->GetSolution (iPoint, 0 );
657656 }
657+ }
658658
659- for (iVar = 0 ; iVar < nVar; iVar++) {
660- if (first_order)
661- LinSysRes (iPoint, iVar) +=
662- (Density_nP1 * U_time_nP1[iVar] - Density_n * U_time_n[iVar]) * (Volume_nP1 / TimeStep);
663- if (second_order)
664- LinSysRes (iPoint, iVar) +=
665- (Density_nP1 * U_time_nP1[iVar] - Density_n * U_time_n[iVar]) * (3.0 * Volume_nP1 / (2.0 * TimeStep)) +
666- (Density_nM1 * U_time_nM1[iVar] - Density_n * U_time_n[iVar]) * (Volume_nM1 / (2.0 * TimeStep));
667- }
668-
669- } else {
670- for (iVar = 0 ; iVar < nVar; iVar++) {
671- if (first_order) LinSysRes (iPoint, iVar) += (U_time_nP1[iVar] - U_time_n[iVar]) * (Volume_nP1 / TimeStep);
672- if (second_order)
673- LinSysRes (iPoint, iVar) += (U_time_nP1[iVar] - U_time_n[iVar]) * (3.0 * Volume_nP1 / (2.0 * TimeStep)) +
674- (U_time_nM1[iVar] - U_time_n[iVar]) * (Volume_nM1 / (2.0 * TimeStep));
675- }
659+ for (iVar = 0 ; iVar < nVar; iVar++) {
660+ if (first_order)
661+ LinSysRes (iPoint, iVar) +=
662+ (Density_nP1 * U_time_nP1[iVar] - Density_n * U_time_n[iVar]) * (Volume_nP1 / TimeStep);
663+ if (second_order)
664+ LinSysRes (iPoint, iVar) +=
665+ (Density_nP1 * U_time_nP1[iVar] - Density_n * U_time_n[iVar]) * (3.0 * Volume_nP1 / (2.0 * TimeStep)) +
666+ (Density_nM1 * U_time_nM1[iVar] - Density_n * U_time_n[iVar]) * (Volume_nM1 / (2.0 * TimeStep));
676667 }
677668
678669 /* --- Compute the Jacobian contribution due to the dual time source term. ---*/
0 commit comments