diff --git a/SU2_CFD/src/solvers/CEulerSolver.cpp b/SU2_CFD/src/solvers/CEulerSolver.cpp
index e5737ceebb1..f7026201d0d 100644
--- a/SU2_CFD/src/solvers/CEulerSolver.cpp
+++ b/SU2_CFD/src/solvers/CEulerSolver.cpp
@@ -7514,6 +7514,7 @@ void CEulerSolver::BC_Supersonic_Outlet(CGeometry *geometry, CSolver **solver_co
   su2double *V_outlet, *V_domain;
 
   bool implicit = (config->GetKind_TimeIntScheme() == EULER_IMPLICIT);
+  bool viscous = config->GetViscous();
   string Marker_Tag = config->GetMarker_All_TagBound(val_marker);
 
   auto *Normal = new su2double[nDim];
@@ -7578,44 +7579,45 @@ void CEulerSolver::BC_Supersonic_Outlet(CGeometry *geometry, CSolver **solver_co
       if (implicit)
         Jacobian.AddBlock2Diag(iPoint, residual.jacobian_i);
 
-//      /*--- Viscous contribution, commented out because serious convergence problems ---*/
-//
-//      if (viscous) {
-//
-//        /*--- Set laminar and eddy viscosity at the infinity ---*/
-//
-//        V_outlet[nDim+5] = nodes->GetLaminarViscosity(iPoint);
-//        V_outlet[nDim+6] = nodes->GetEddyViscosity(iPoint);
-//
-//        /*--- Set the normal vector and the coordinates ---*/
-//
-//        visc_numerics->SetNormal(Normal);
-//        su2double Coord_Reflected[MAXNDIM];
-//        GeometryToolbox::PointPointReflect(nDim, geometry->nodes->GetCoord(Point_Normal),
-//                                                 geometry->nodes->GetCoord(iPoint), Coord_Reflected);
-//        visc_numerics->SetCoord(geometry->nodes->GetCoord(iPoint), Coord_Reflected);
-//
-//        /*--- Primitive variables, and gradient ---*/
-//
-//        visc_numerics->SetPrimitive(V_domain, V_outlet);
-//        visc_numerics->SetPrimVarGradient(nodes->GetGradient_Primitive(iPoint), nodes->GetGradient_Primitive(iPoint));
-//
-//        /*--- Turbulent kinetic energy ---*/
-//
-//        if (config->GetKind_Turb_Model() == TURB_MODEL::SST)
-//          visc_numerics->SetTurbKineticEnergy(solver_container[TURB_SOL]->GetNodes()->GetSolution(iPoint,0),
-//                                              solver_container[TURB_SOL]->GetNodes()->GetSolution(iPoint,0));
-//
-//        /*--- Compute and update residual ---*/
-//
-//        auto residual = visc_numerics->ComputeResidual(config);
-//        LinSysRes.SubtractBlock(iPoint, residual);
-//
-//        /*--- Jacobian contribution for implicit integration ---*/
-//
-//        if (implicit)
-//          Jacobian.SubtractBlock2Diag(iPoint, residual.jacobian_i);
-//      }
+     /*--- Viscous contribution, commented out because serious convergence problems ---*/
+
+     if (viscous) {
+
+       /*--- Set laminar and eddy viscosity at the infinity ---*/
+
+       V_outlet[nDim+5] = nodes->GetLaminarViscosity(iPoint);
+       V_outlet[nDim+6] = nodes->GetEddyViscosity(iPoint);
+
+       /*--- Set the normal vector and the coordinates ---*/
+
+       visc_numerics->SetNormal(Normal);
+       su2double Coord_Reflected[MAXNDIM];
+       unsigned long Point_Normal = geometry->vertex[val_marker][iVertex]->GetNormal_Neighbor();
+       GeometryToolbox::PointPointReflect(nDim, geometry->nodes->GetCoord(Point_Normal),
+                                                geometry->nodes->GetCoord(iPoint), Coord_Reflected);
+       visc_numerics->SetCoord(geometry->nodes->GetCoord(iPoint), Coord_Reflected);
+
+       /*--- Primitive variables, and gradient ---*/
+
+       visc_numerics->SetPrimitive(V_domain, V_domain);
+       visc_numerics->SetPrimVarGradient(nodes->GetGradient_Primitive(iPoint), nodes->GetGradient_Primitive(iPoint));
+
+       /*--- Turbulent kinetic energy ---*/
+
+       if (config->GetKind_Turb_Model() == TURB_MODEL::SST)
+         visc_numerics->SetTurbKineticEnergy(solver_container[TURB_SOL]->GetNodes()->GetSolution(iPoint,0),
+                                             solver_container[TURB_SOL]->GetNodes()->GetSolution(iPoint,0));
+
+       /*--- Compute and update residual ---*/
+
+       auto residual = visc_numerics->ComputeResidual(config);
+       LinSysRes.SubtractBlock(iPoint, residual);
+
+       /*--- Jacobian contribution for implicit integration ---*/
+
+       if (implicit)
+         Jacobian.SubtractBlock2Diag(iPoint, residual.jacobian_i);
+     }
 
     }
   }
diff --git a/TestCases/navierstokes/flatplate/flatplate_supersonic.cfg b/TestCases/navierstokes/flatplate/flatplate_supersonic.cfg
new file mode 100644
index 00000000000..f7798c39f78
--- /dev/null
+++ b/TestCases/navierstokes/flatplate/flatplate_supersonic.cfg
@@ -0,0 +1,72 @@
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%                                                                              %
+% SU2 configuration file                                                       %
+% Case description: Supersonic flat plate flow                                 %
+% Author: Tan Yilun                                                            %
+% Institution: University of Chinese Academy of Sciences                       %
+% Date: 2026.05.22                                                             %
+% File Version 8.5.1                                                           %
+%                                                                              %
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+SOLVER= NAVIER_STOKES       
+KIND_TURB_MODEL= NONE       
+MATH_PROBLEM= DIRECT        
+RESTART_SOL= YES       
+SOLUTION_FILENAME= restart_flow.dat
+
+MESH_FILENAME= flatplate_supersonic.cgns
+MESH_FORMAT= CGNS
+
+TABULAR_FORMAT= TECPLOT
+CONV_FILENAME= history
+RESTART_FILENAME= restart_flow
+VOLUME_FILENAME= volume_flow
+SURFACE_FILENAME= surface_flow
+OUTPUT_WRT_FREQ= 500       
+WRT_RESTART_OVERWRITE= YES   
+WRT_SURFACE_OVERWRITE= NO   
+WRT_VOLUME_OVERWRITE= NO  
+OUTPUT_FILES= (RESTART, TECPLOT, SURFACE_CSV)
+SCREEN_OUTPUT=(INNER_ITER, RMS_DENSITY, RMS_MOMENTUM-X, RMS_MOMENTUM-Y, RMS_ENERGY, DRAG)
+
+MACH_NUMBER= 4.0           
+AOA= 0.0                   
+SIDESLIP_ANGLE= 0.0         
+INIT_OPTION= TD_CONDITIONS  
+FREESTREAM_TEMPERATURE= 300.0
+REYNOLDS_NUMBER= 12314
+FREESTREAM_PRESSURE= 70.47  
+REYNOLDS_LENGTH= 0.2       
+
+MARKER_ISOTHERMAL= ( WALL, 300.0 )   
+MARKER_SYM= ( SYMMETRY )              
+MARKER_FAR= ( INLET )                 
+MARKER_SUPERSONIC_OUTLET= ( OUTLET ) 
+MARKER_PLOTTING= ( WALL )            
+MARKER_MONITORING= ( WALL )           
+
+NUM_METHOD_GRAD= WEIGHTED_LEAST_SQUARES
+
+LINEAR_SOLVER= BCGSTAB
+LINEAR_SOLVER_PREC= LU_SGS
+LINEAR_SOLVER_ERROR= 1E-10
+LINEAR_SOLVER_ITER= 10
+
+MGLEVEL= 0              
+
+CONV_NUM_METHOD_FLOW= AUSM  
+MUSCL_FLOW= YES             
+SLOPE_LIMITER_FLOW= VENKATAKRISHNAN_WANG
+VENKAT_LIMITER_COEFF= 0.05
+
+CFL_NUMBER= 0.1             
+CFL_ADAPT= YES             
+CFL_ADAPT_PARAM= ( 0.3, 1.05, 0.01, 5)
+ITER= 1000             
+TIME_DISCRE_FLOW= EULER_IMPLICIT
+
+CONV_FIELD= RMS_ENERGY
+CONV_RESIDUAL_MINVAL= -5
+CONV_STARTITER= 10
+
diff --git a/TestCases/parallel_regression.py b/TestCases/parallel_regression.py
index a6e77734224..1d241d324f7 100755
--- a/TestCases/parallel_regression.py
+++ b/TestCases/parallel_regression.py
@@ -330,6 +330,14 @@ def main():
     flatplate.test_vals = [-6.499563, -1.020731, 0.001223, 0.028373, 2.361500, -2.333200, 0.000000, 0.000000]
     test_list.append(flatplate)
 
+    # Supersonic laminar flat plate
+    flatplate_supersonic           = TestCase('flatplate_supersonic')
+    flatplate_supersonic.cfg_dir   = "navierstokes/flatplate"
+    flatplate_supersonic.cfg_file  = "flatplate_supersonic.cfg"
+    flatplate_supersonic.test_iter = 100
+    flatplate_supersonic.test_vals = [-6.397042, -0.919335, 0.001321, 0.025064, 2.361500, -2.336500, 0.000000, 0.000000]
+    test_list.append(flatplate_supersonic)
+
     # Custom objective function
     flatplate_udobj           = TestCase('flatplate_udobj')
     flatplate_udobj.cfg_dir   = "user_defined_functions"