diff --git a/Common/src/CConfig.cpp b/Common/src/CConfig.cpp
index cf06f7379746..2a86786b02fa 100644
--- a/Common/src/CConfig.cpp
+++ b/Common/src/CConfig.cpp
@@ -4602,6 +4602,36 @@ void CConfig::SetPostprocessing(SU2_COMPONENT val_software, unsigned short val_i
 
   Rotating_Frame = (Kind_GridMovement == ROTATING_FRAME);
 
+  /*--- The finite volume flow solvers apply these source terms with a single numerics object, which is built
+   for only one of them: combinations would apply one source several times and ignore the others. Only these
+   solvers (and their discrete adjoints) are checked: e.g. elasticity uses BODY_FORCE and GRAVITY_FORCE together. ---*/
+
+  if (GetFluidProblem() && !GetNEMOProblem()) {
+    const bool inc = (Kind_Regime == ENUM_REGIME::INCOMPRESSIBLE);
+    const int nSources = int(Body_Force) + int(Rotating_Frame) + int(Axisymmetric) + int(GravityForce && !inc) +
+                         int(inc && Kind_DensityModel == INC_DENSITYMODEL::BOUSSINESQ) +
+                         int(inc && Kind_Streamwise_Periodic != ENUM_STREAMWISE_PERIODIC::NONE);
+    if (nSources > 1) {
+      SU2_MPI::Error("Only one of BODY_FORCE, GRID_MOVEMENT= ROTATING_FRAME, AXISYMMETRIC, GRAVITY_FORCE (compressible), "
+                     "INC_DENSITY_MODEL= BOUSSINESQ and KIND_STREAMWISE_PERIODIC (incompressible) can be used at a time.",
+                     CURRENT_FUNCTION);
+    }
+  }
+
+  /*--- The second source object, for radiation, the energy correction of incompressible streamwise periodic flow
+   (without periodic temperature) or vorticity confinement, is also built for only one of them. ---*/
+
+  if (GetFluidProblem() && !GetNEMOProblem()) {
+    const bool inc = (Kind_Regime == ENUM_REGIME::INCOMPRESSIBLE);
+    const bool periodicOutlet = inc && Kind_Streamwise_Periodic != ENUM_STREAMWISE_PERIODIC::NONE && Energy_Equation &&
+                                !Streamwise_Periodic_Temperature;
+    if (int(Radiation) + int(periodicOutlet) + int(VorticityConfinement) > 1) {
+      SU2_MPI::Error("Only one of RADIATION, the energy correction of KIND_STREAMWISE_PERIODIC (incompressible, "
+                     "STREAMWISE_PERIODIC_TEMPERATURE= NO) and VORTICITY_CONFINEMENT can be used at a time.",
+                     CURRENT_FUNCTION);
+    }
+  }
+
   /*--- In case the grid movement parameters have not been declared in the
    config file, set them equal to zero for safety. Also check to make sure
    that for each option, a value has been declared for each moving marker. ---*/
@@ -4961,6 +4991,11 @@ void CConfig::SetPostprocessing(SU2_COMPONENT val_software, unsigned short val_i
 
   if (Restart || ((Kind_MGCycle == MG_CYCLE::FULL) && ContinuousAdjoint)) Kind_MGCycle = MG_CYCLE::V;
 
+  /*--- The multigrid arrays hold at most 10 coarse levels (CMultiGridIntegration::MAX_MG_LEVELS). ---*/
+  if (nMGLevels > 10) {
+    SU2_MPI::Error("MGLEVEL must be between 0 and 10; larger multigrid hierarchies are not supported.", CURRENT_FUNCTION);
+  }
+
   FinestMesh = MESH_0;
   if (Kind_MGCycle == MG_CYCLE::FULL) FinestMesh = nMGLevels;
 
diff --git a/SU2_CFD/include/gradients/computeGradientsGreenGauss.hpp b/SU2_CFD/include/gradients/computeGradientsGreenGauss.hpp
index f60c0422f0cc..89eb028c6aeb 100644
--- a/SU2_CFD/include/gradients/computeGradientsGreenGauss.hpp
+++ b/SU2_CFD/include/gradients/computeGradientsGreenGauss.hpp
@@ -69,7 +69,14 @@ void computeGradientsGreenGauss(CSolver* solver, MPI_QUANTITIES kindMpiComm, PER
   const auto chunkSize = computeStaticChunkSize(nPointDomain, omp_get_max_threads(), OMP_MAX_CHUNK);
 #endif
 
-  static constexpr size_t MAXNVAR = 20;
+  static constexpr size_t MAXNVAR = 32;
+
+  /*--- Called by all threads of a parallel region: only the master may call MPI (MPI_THREAD_FUNNELED). ---*/
+  if (varEnd > MAXNVAR) {
+    BEGIN_SU2_OMP_SAFE_GLOBAL_ACCESS
+    SU2_MPI::Error("Number of variables is too large, increase MAXNVAR.", CURRENT_FUNCTION);
+    END_SU2_OMP_SAFE_GLOBAL_ACCESS
+  }
 
   /*--- For each (non-halo) volume integrate over its faces (edges). ---*/
 
diff --git a/SU2_CFD/include/limiters/computeLimiters_impl.hpp b/SU2_CFD/include/limiters/computeLimiters_impl.hpp
index 3b19f241f3bc..cd45163146ef 100644
--- a/SU2_CFD/include/limiters/computeLimiters_impl.hpp
+++ b/SU2_CFD/include/limiters/computeLimiters_impl.hpp
@@ -80,8 +80,12 @@ void computeLimiters_impl(CSolver* solver,
 {
   constexpr size_t MAXNVAR = 32;
 
-  if (varEnd > MAXNVAR)
+  /*--- Called by all threads of a parallel region: only the master may call MPI (MPI_THREAD_FUNNELED). ---*/
+  if (varEnd > MAXNVAR) {
+    BEGIN_SU2_OMP_SAFE_GLOBAL_ACCESS
     SU2_MPI::Error("Number of variables is too large, increase MAXNVAR.", CURRENT_FUNCTION);
+    END_SU2_OMP_SAFE_GLOBAL_ACCESS
+  }
 
   const size_t nPointDomain = geometry.GetnPointDomain();
   const size_t nPoint = geometry.GetnPoint();
diff --git a/SU2_CFD/include/solvers/CFVMFlowSolverBase.inl b/SU2_CFD/include/solvers/CFVMFlowSolverBase.inl
index 3a46fc8c4db6..367ece9e6a5d 100644
--- a/SU2_CFD/include/solvers/CFVMFlowSolverBase.inl
+++ b/SU2_CFD/include/solvers/CFVMFlowSolverBase.inl
@@ -1155,6 +1155,7 @@ void CFVMFlowSolverBase<V, FlowRegime>::BC_Sym_Plane(CGeometry* geometry, CSolve
   const bool ideal_gas = (config->GetKind_FluidModel() == STANDARD_AIR) ||
                          (config->GetKind_FluidModel() == IDEAL_GAS);
   const auto iVel = prim_idx.Velocity();
+  const bool compressibleProjection = FlowRegime == ENUM_REGIME::COMPRESSIBLE && iVel == 1 && nVar == nDim + 2;
 
   /*--- Blazek chapter 8.:
    * The components of the momentum residual normal to the symmetry plane are zeroed out.
@@ -1184,8 +1185,8 @@ void CFVMFlowSolverBase<V, FlowRegime>::BC_Sym_Plane(CGeometry* geometry, CSolve
       for (auto iDim = 0u; iDim < nDim; iDim++) UnitNormal[iDim] = Normal[iDim] / Area;
     }
 
-    /*--- Energy terms due to grid movement (aka work of pressure forces). ---*/
-    if (dynamic_grid) {
+    /*--- The classical compressible projection below already supplies wall pressure work. ---*/
+    if (dynamic_grid && !compressibleProjection) {
       su2double* V_reflected = GetCharacPrimVar(val_marker, iVertex);
 
       conv_numerics->SetGridVel(geometry->nodes->GetGridVel(iPoint),
@@ -1217,8 +1218,7 @@ void CFVMFlowSolverBase<V, FlowRegime>::BC_Sym_Plane(CGeometry* geometry, CSolve
       /*--- Compute the residual using an upwind scheme. ---*/
       auto residual = conv_numerics->ComputeResidual(config);
 
-      /*--- Use just the energy fluxes to update the residual, adding the others would
-       * increase numerical diffusion which we wish to avoid if possible. ---*/
+      /*--- Energy fluxes for the other conservative layouts. ---*/
       for (auto iVar = iVel + nDim; iVar < nVar; iVar++) {
         LinSysRes(iPoint, iVar) += residual.residual[iVar];
       }
@@ -1243,14 +1243,16 @@ void CFVMFlowSolverBase<V, FlowRegime>::BC_Sym_Plane(CGeometry* geometry, CSolve
     su2double* solutionOld = nodes->GetSolution_Old(iPoint);
 
     su2double gridVel[MAXNDIM] = {};
+    su2double gridNormalVel = 0.0;
     if (dynamic_grid) {
       for (auto iDim = 0u; iDim < nDim; iDim++) {
         gridVel[iDim] = geometry->nodes->GetGridVel(iPoint)[iDim];
       }
+      gridNormalVel = GeometryToolbox::DotProduct(nDim, gridVel, UnitNormal);
       if (FlowRegime == ENUM_REGIME::COMPRESSIBLE) {
         for(auto iDim = 0u; iDim < nDim; iDim++) {
           /*--- Multiply by density since we are correcting conservative variables. ---*/
-          gridVel[iDim] *= nodes->GetDensity(iPoint);
+          gridVel[iDim] *= compressibleProjection ? solutionOld[0] : nodes->GetDensity(iPoint);
         }
       }
     }
@@ -1261,20 +1263,27 @@ void CFVMFlowSolverBase<V, FlowRegime>::BC_Sym_Plane(CGeometry* geometry, CSolve
     for (auto iDim = 0u; iDim < nDim; iDim++) {
       solutionOld[iVel + iDim] -= vp * UnitNormal[iDim];
     }
+    if (compressibleProjection) {
+      solutionOld[iVel + nDim] -= gridNormalVel * vp;
+    }
 
     /*--- Keep only the tangential part of the momentum residuals. ---*/
     su2double normalRes = 0.0;
     for (auto iDim = 0u; iDim < nDim; iDim++) {
       normalRes += LinSysRes(iPoint, iVel + iDim) * UnitNormal[iDim];
     }
+    if (compressibleProjection) {
+      /*--- Project in the wall frame: normal momentum follows density, with consistent pressure work. ---*/
+      normalRes -= gridNormalVel * LinSysRes(iPoint, 0);
+      LinSysRes(iPoint, iVel + nDim) -= gridNormalVel * normalRes;
+    }
     for (auto iDim = 0u; iDim < nDim; iDim++) {
       LinSysRes(iPoint, iVel + iDim) -= normalRes * UnitNormal[iDim];
     }
 
     /*--- Jacobian contribution for implicit integration. ---*/
     if (implicit) {
-      /*--- Modify the Jacobians according to the modification of the residual
-       * J_new = (I - n * n^T) * J where n = {0, nx, ny, nz, 0, ...} ---*/
+      /*--- Apply the same wall-frame projection to the Jacobians. ---*/
       su2double mat[MAXNVAR * MAXNVAR] = {};
 
       for (auto iVar = 0u; iVar < nVar; iVar++)
@@ -1282,6 +1291,13 @@ void CFVMFlowSolverBase<V, FlowRegime>::BC_Sym_Plane(CGeometry* geometry, CSolve
       for (auto iDim = 0u; iDim < nDim; iDim++)
         for (auto jDim = 0u; jDim < nDim; jDim++)
           mat[(iDim + iVel) * nVar + jDim + iVel] -= UnitNormal[iDim] * UnitNormal[jDim];
+      if (compressibleProjection) {
+        for (auto iDim = 0u; iDim < nDim; iDim++) {
+          mat[(iDim + iVel) * nVar] += gridNormalVel * UnitNormal[iDim];
+          mat[(iVel + nDim) * nVar + iVel + iDim] -= gridNormalVel * UnitNormal[iDim];
+        }
+        mat[(iVel + nDim) * nVar] += gridNormalVel * gridNormalVel;
+      }
 
       auto ModifyJacobian = [&](const unsigned long jPoint) {
         su2double jac[MAXNVAR * MAXNVAR], newJac[MAXNVAR * MAXNVAR];
@@ -1306,6 +1322,10 @@ void CFVMFlowSolverBase<V, FlowRegime>::BC_Sym_Plane(CGeometry* geometry, CSolve
     for (auto iDim = 0u; iDim < nDim; iDim++) {
       normalRes += Res_TruncError[iVel + iDim] * UnitNormal[iDim];
     }
+    if (compressibleProjection) {
+      normalRes -= gridNormalVel * Res_TruncError[0];
+      Res_TruncError[iVel + nDim] -= gridNormalVel * normalRes;
+    }
     for (auto iDim = 0u; iDim < nDim; iDim++) {
       Res_TruncError[iVel + iDim] -= normalRes * UnitNormal[iDim];
     }
@@ -2974,6 +2994,7 @@ void CFVMFlowSolverBase<V, FlowRegime>::MultigridProjectEulerWall(CGeometry* geo
                                                                    bool use_solution_old) {
   const auto iVel = prim_idx.Velocity();
   const auto nDim = geometry->GetnDim();
+  const bool compressibleProjection = FlowRegime == ENUM_REGIME::COMPRESSIBLE && iVel == 1 && nVar == nDim + 2;
 
   for (auto iMarker = 0u; iMarker < config->GetnMarker_All(); iMarker++) {
     if (config->GetMarker_All_KindBC(iMarker) != EULER_WALL) continue;
@@ -3004,12 +3025,15 @@ void CFVMFlowSolverBase<V, FlowRegime>::MultigridProjectEulerWall(CGeometry* geo
        *    For dynamic grids subtract the grid velocity to enforce (v - v_grid).n = 0,
        *    multiplying by density for compressible flow (conservative variables). ---*/
       su2double gridVel[MAXNDIM] = {};
-      if (dynamic_grid && !use_solution_old) {
+      su2double gridNormalVel = 0.0;
+      if (dynamic_grid && (compressibleProjection || !use_solution_old)) {
         for (auto iDim = 0u; iDim < nDim; iDim++)
           gridVel[iDim] = geometry->nodes->GetGridVel(iPoint)[iDim];
+        gridNormalVel = GeometryToolbox::DotProduct(nDim, gridVel, UnitNormal);
         if constexpr (FlowRegime == ENUM_REGIME::COMPRESSIBLE) {
-          for (auto iDim = 0u; iDim < nDim; iDim++)
-            gridVel[iDim] *= nodes->GetDensity(iPoint);
+          /*--- Corrections carry their own density increment, rather than the current solution density. ---*/
+          const auto density = compressibleProjection ? sol[0] : nodes->GetDensity(iPoint);
+          for (auto iDim = 0u; iDim < nDim; iDim++) gridVel[iDim] *= density;
         }
       }
 
@@ -3019,6 +3043,7 @@ void CFVMFlowSolverBase<V, FlowRegime>::MultigridProjectEulerWall(CGeometry* geo
 
       /*--- Project to tangent plane. ---*/
       for (auto iDim = 0u; iDim < nDim; iDim++) sol[iVel + iDim] -= momentum_n * UnitNormal[iDim];
+      if (compressibleProjection) sol[iVel + nDim] -= gridNormalVel * momentum_n;
     }
     END_SU2_OMP_FOR
   }
diff --git a/SU2_CFD/src/numerics/flow/convection/centered.cpp b/SU2_CFD/src/numerics/flow/convection/centered.cpp
index b0130279f76f..713ca2e86662 100644
--- a/SU2_CFD/src/numerics/flow/convection/centered.cpp
+++ b/SU2_CFD/src/numerics/flow/convection/centered.cpp
@@ -31,7 +31,7 @@
 CCentLaxInc_Flow::CCentLaxInc_Flow(unsigned short val_nDim, unsigned short val_nVar, const CConfig* config) : CNumerics(val_nDim, val_nVar, config) {
 
   implicit         = (config->GetKind_TimeIntScheme_Flow() == EULER_IMPLICIT);
-  variable_density = ((config->GetKind_DensityModel() == INC_DENSITYMODEL::VARIABLE) || (config->GetKind_DensityModel() == INC_DENSITYMODEL::FLAMELET));
+  variable_density = config->GetVariable_Density_Model();
   /* A grid is defined as dynamic if there's rigid grid movement or grid deformation AND the problem is time domain */
   dynamic_grid = config->GetDynamic_Grid();
   fix_factor = config->GetCent_Inc_Jac_Fix_Factor();
@@ -161,17 +161,25 @@ CNumerics::ResidualType<> CCentLaxInc_Flow::ComputeResidual(const CConfig* confi
       ProjVelocity += 0.5*(GridVel_i[iDim]+GridVel_j[iDim])*Normal[iDim];
 
     /*--- Residual contributions ---*/
-    for (iVar = 0; iVar < nVar; iVar++) {
+    for (iVar = 0; iVar < nVar; iVar++)
       ProjFlux[iVar] -= ProjVelocity * 0.5*(U_i[iVar]+U_j[iVar]);
 
-      /*--- Jacobian contributions ---*/
-      /*--- Implicit terms ---*/
-      if (implicit) {
-        for (iDim = 0; iDim < nDim+1; iDim++){
-          Jacobian_i[iDim+1][iDim+1] -= 0.5*ProjVelocity*DensityInc_i;
-          Jacobian_j[iDim+1][iDim+1] -= 0.5*ProjVelocity*DensityInc_j;
+    /*--- Jacobian of -ProjVelocity * 0.5*U, with U = (rho, rho*u, rho*h), w.r.t. the primitives (p, u, h), added
+     once. The density depends on h (variable density models) and not on the pressure. ---*/
+    if (implicit) {
+      auto addGridJacobian = [&](su2double **Jac, su2double rho, const su2double *vel, su2double h, su2double T,
+                                 su2double cp) {
+        const auto dRhodh = variable_density ? su2double(-rho / (T * cp)) : su2double(0.0);
+        const su2double factor = 0.5*ProjVelocity;
+        Jac[0][nDim+1] -= factor*dRhodh;
+        for (iDim = 0; iDim < nDim; iDim++) {
+          Jac[iDim+1][iDim+1] -= factor*rho;
+          Jac[iDim+1][nDim+1] -= factor*vel[iDim]*dRhodh;
         }
-      }
+        Jac[nDim+1][nDim+1] -= factor*(rho + h*dRhodh);
+      };
+      addGridJacobian(Jacobian_i, DensityInc_i, Velocity_i, Enthalpy_i, Temperature_i, Cp_i);
+      addGridJacobian(Jacobian_j, DensityInc_j, Velocity_j, Enthalpy_j, Temperature_j, Cp_j);
     }
   }
 
@@ -250,7 +258,7 @@ CNumerics::ResidualType<> CCentLaxInc_Flow::ComputeResidual(const CConfig* confi
 CCentJSTInc_Flow::CCentJSTInc_Flow(unsigned short val_nDim, unsigned short val_nVar, const CConfig* config) : CNumerics(val_nDim, val_nVar, config) {
 
   implicit         = (config->GetKind_TimeIntScheme_Flow() == EULER_IMPLICIT);
-  variable_density = (config->GetVariable_Density_Model());
+  variable_density = config->GetVariable_Density_Model();
   energy           = config->GetEnergy_Equation();
   /* A grid is defined as dynamic if there's rigid grid movement or grid deformation AND the problem is time domain */
   dynamic_grid = config->GetDynamic_Grid();
@@ -401,17 +409,25 @@ CNumerics::ResidualType<> CCentJSTInc_Flow::ComputeResidual(const CConfig* confi
       ProjVelocity += 0.5*(GridVel_i[iDim]+GridVel_j[iDim])*Normal[iDim];
 
     /*--- Residual contributions ---*/
-    for (iVar = 0; iVar < nVar; iVar++) {
+    for (iVar = 0; iVar < nVar; iVar++)
       ProjFlux[iVar] -= ProjVelocity * 0.5*(U_i[iVar]+U_j[iVar]);
 
-      /*--- Jacobian contributions ---*/
-      /*--- Implicit terms ---*/
-      if (implicit) {
-        for (iDim = 0; iDim < nDim+1; iDim++){
-          Jacobian_i[iDim+1][iDim+1] -= 0.5*ProjVelocity*DensityInc_i;
-          Jacobian_j[iDim+1][iDim+1] -= 0.5*ProjVelocity*DensityInc_j;
+    /*--- Jacobian of -ProjVelocity * 0.5*U, with U = (rho, rho*u, rho*h), w.r.t. the primitives (p, u, h), added
+     once. The density depends on h (variable density models) and not on the pressure. ---*/
+    if (implicit) {
+      auto addGridJacobian = [&](su2double **Jac, su2double rho, const su2double *vel, su2double h, su2double T,
+                                 su2double cp) {
+        const auto dRhodh = variable_density ? su2double(-rho / (T * cp)) : su2double(0.0);
+        const su2double factor = 0.5*ProjVelocity;
+        Jac[0][nDim+1] -= factor*dRhodh;
+        for (iDim = 0; iDim < nDim; iDim++) {
+          Jac[iDim+1][iDim+1] -= factor*rho;
+          Jac[iDim+1][nDim+1] -= factor*vel[iDim]*dRhodh;
         }
-      }
+        Jac[nDim+1][nDim+1] -= factor*(rho + h*dRhodh);
+      };
+      addGridJacobian(Jacobian_i, DensityInc_i, Velocity_i, Enthalpy_i, Temperature_i, Cp_i);
+      addGridJacobian(Jacobian_j, DensityInc_j, Velocity_j, Enthalpy_j, Temperature_j, Cp_j);
     }
   }
 
diff --git a/SU2_CFD/src/numerics/flow/convection/fds.cpp b/SU2_CFD/src/numerics/flow/convection/fds.cpp
index e806127da203..57c7ed9f6f5a 100644
--- a/SU2_CFD/src/numerics/flow/convection/fds.cpp
+++ b/SU2_CFD/src/numerics/flow/convection/fds.cpp
@@ -249,19 +249,23 @@ CNumerics::ResidualType<> CUpwFDSInc_Flow::ComputeResidual(const CConfig *config
       ProjVelocity += 0.5*(GridVel_i[iDim]+GridVel_j[iDim])*Normal[iDim];
 
     /*--- Residual contributions ---*/
-    for (iVar = 0; iVar < nVar; iVar++) {
+    for (iVar = 0; iVar < nVar; iVar++)
       Flux[iVar] -= ProjVelocity * 0.5*(U_i[iVar]+U_j[iVar]);
 
-      /*--- Jacobian contributions ---*/
-      /*--- Implicit terms ---*/
-      if (implicit) {
-        for (iDim = 0; iDim < nDim; iDim++){
-          Jacobian_i[iDim+1][iDim+1] -= 0.5*ProjVelocity*DensityInc_i;
-          Jacobian_j[iDim+1][iDim+1] -= 0.5*ProjVelocity*DensityInc_j;
+    /*--- Jacobian of -ProjVelocity * 0.5*U, with U = (rho, rho*u, rho*h), w.r.t. the primitives (p, u, h), added
+     once. The density depends on h (variable density models) and not on the pressure. ---*/
+    if (implicit) {
+      auto addGridJacobian = [&](su2double **Jac, su2double rho, const su2double *vel, su2double h, su2double dRhodh) {
+        const su2double factor = 0.5*ProjVelocity;
+        Jac[0][nDim+1] -= factor*dRhodh;
+        for (iDim = 0; iDim < nDim; iDim++) {
+          Jac[iDim+1][iDim+1] -= factor*rho;
+          Jac[iDim+1][nDim+1] -= factor*vel[iDim]*dRhodh;
         }
-        Jacobian_i[nDim+1][nDim+1] -= 0.5*ProjVelocity*DensityInc_i;
-        Jacobian_j[nDim+1][nDim+1] -= 0.5*ProjVelocity*DensityInc_j;
-      }
+        Jac[nDim+1][nDim+1] -= factor*(rho + h*dRhodh);
+      };
+      addGridJacobian(Jacobian_i, DensityInc_i, Velocity_i, Enthalpy_i, dRhodh_i);
+      addGridJacobian(Jacobian_j, DensityInc_j, Velocity_j, Enthalpy_j, dRhodh_j);
     }
   }
 
diff --git a/SU2_CFD/src/numerics/flow/flow_sources.cpp b/SU2_CFD/src/numerics/flow/flow_sources.cpp
index e3797c7c5a5f..87784e73f39f 100644
--- a/SU2_CFD/src/numerics/flow/flow_sources.cpp
+++ b/SU2_CFD/src/numerics/flow/flow_sources.cpp
@@ -138,7 +138,7 @@ CNumerics::ResidualType<> CSourceAxisymmetric_Flow::ComputeResidual(const CConfi
 
       jacobian[3][0] = -Gamma*U_i[2]*U_i[3]/(U_i[0]*U_i[0]) + (Gamma-1)*U_i[2]*(U_i[1]*U_i[1]+U_i[2]*U_i[2])/(U_i[0]*U_i[0]*U_i[0]);
       jacobian[3][1] = -(Gamma-1)*U_i[2]*U_i[1]/(U_i[0]*U_i[0]);
-      jacobian[3][2] = Gamma*U_i[3]/U_i[0] - 1/2*(Gamma-1)*( (U_i[1]*U_i[1]+U_i[2]*U_i[2])/(U_i[0]*U_i[0]) + 2*U_i[2]*U_i[2]/(U_i[0]*U_i[0]) );
+      jacobian[3][2] = Gamma*U_i[3]/U_i[0] - 0.5*(Gamma-1)*( (U_i[1]*U_i[1]+U_i[2]*U_i[2])/(U_i[0]*U_i[0]) + 2*U_i[2]*U_i[2]/(U_i[0]*U_i[0]) );
       jacobian[3][3] = Gamma*U_i[2]/U_i[0];
 
       for (iVar=0; iVar < nVar; iVar++)
diff --git a/SU2_CFD/src/solvers/CEulerSolver.cpp b/SU2_CFD/src/solvers/CEulerSolver.cpp
index b1c61ad377de..c05a097c4a2f 100644
--- a/SU2_CFD/src/solvers/CEulerSolver.cpp
+++ b/SU2_CFD/src/solvers/CEulerSolver.cpp
@@ -1947,27 +1947,39 @@ void CEulerSolver::Upwind_Residual(CGeometry *geometry, CSolver **solver_contain
 
       }
 
-      /*--- Recompute the reconstructed quantities in a thermodynamically consistent way. ---*/
+      /*--- Check for non-physical solutions after reconstruction. If found, use the
+       cell-average value of the solution. This is a locally 1st order approximation,
+       which is typically only active during the start-up of a calculation.
+       Pressure and density are checked before the fluid model is evaluated with them. ---*/
 
-      if (!ideal_gas || low_mach_corr) {
-        ComputeConsistentExtrapolation(GetFluidModel(), nDim, Primitive_i, Secondary_i);
-        ComputeConsistentExtrapolation(GetFluidModel(), nDim, Primitive_j, Secondary_j);
-      }
+      auto nonPositive = [](su2double value) { return !(value > 0.0) || !std::isfinite(SU2_TYPE::GetValue(value)); };
+      bool neg_pres_or_rho_i = nonPositive(Primitive_i[prim_idx.Pressure()]) || nonPositive(Primitive_i[prim_idx.Density()]);
+      bool neg_pres_or_rho_j = nonPositive(Primitive_j[prim_idx.Pressure()]) || nonPositive(Primitive_j[prim_idx.Density()]);
 
-      /*--- Low-Mach number correction. ---*/
+      if (!neg_pres_or_rho_i && !neg_pres_or_rho_j) {
 
-      if (low_mach_corr) {
-        LowMachPrimitiveCorrection(GetFluidModel(), nDim, Primitive_i, Primitive_j);
-      }
+        /*--- Recompute the reconstructed quantities in a thermodynamically consistent way. ---*/
 
-      /*--- Check for non-physical solutions after reconstruction. If found, use the
-       cell-average value of the solution. This is a locally 1st order approximation,
-       which is typically only active during the start-up of a calculation. ---*/
+        if (!ideal_gas || low_mach_corr) {
+          ComputeConsistentExtrapolation(GetFluidModel(), nDim, Primitive_i, Secondary_i);
+          ComputeConsistentExtrapolation(GetFluidModel(), nDim, Primitive_j, Secondary_j);
+        }
+
+        /*--- Low-Mach number correction. ---*/
+
+        if (low_mach_corr) {
+          LowMachPrimitiveCorrection(GetFluidModel(), nDim, Primitive_i, Primitive_j);
+        }
+
+        /*--- The correction changes the face states, check them again. ---*/
 
-      bool neg_pres_or_rho_i = (Primitive_i[prim_idx.Pressure()] < 0.0) || (Primitive_i[prim_idx.Density()] < 0.0);
-      bool neg_pres_or_rho_j = (Primitive_j[prim_idx.Pressure()] < 0.0) || (Primitive_j[prim_idx.Density()] < 0.0);
+        neg_pres_or_rho_i = nonPositive(Primitive_i[prim_idx.Pressure()]) || nonPositive(Primitive_i[prim_idx.Density()]);
+        neg_pres_or_rho_j = nonPositive(Primitive_j[prim_idx.Pressure()]) || nonPositive(Primitive_j[prim_idx.Density()]);
+      }
 
-      su2double R = sqrt(fabs(Primitive_j[prim_idx.Density()]/Primitive_i[prim_idx.Density()]));
+      su2double R = 1.0;
+      if (!neg_pres_or_rho_i && !neg_pres_or_rho_j)
+        R = sqrt(Primitive_j[prim_idx.Density()]/Primitive_i[prim_idx.Density()]);
       su2double sq_vel = 0.0;
       for (auto iDim = 0u; iDim < nDim; iDim++) {
         su2double RoeVelocity = (R * Primitive_j[iDim + prim_idx.Velocity()] +
@@ -5801,7 +5813,7 @@ void CEulerSolver::BC_TurboRiemann(CGeometry *geometry, CSolver **solver_contain
           ProjVelocity_i += Velocity_i[iDim]*UnitNormal[iDim];
 
         su2double donorAverages[5] = {0.0};
-        switch (config->GetKind_Data_Giles(Marker_Tag)){
+        switch (config->GetKind_Data_Riemann(Marker_Tag)){
           case MIXING_IN: case MIXING_IN_1D: case MIXING_OUT: case MIXING_OUT_1D:
             for (auto mixVar = 0u; mixVar < 5; mixVar++) donorAverages[mixVar] = GetMixingState(val_marker, iSpan, mixVar);
             break;
diff --git a/SU2_CFD/src/solvers/CIncEulerSolver.cpp b/SU2_CFD/src/solvers/CIncEulerSolver.cpp
index 11b4e3f9f063..d1f71f00ac99 100644
--- a/SU2_CFD/src/solvers/CIncEulerSolver.cpp
+++ b/SU2_CFD/src/solvers/CIncEulerSolver.cpp
@@ -2245,7 +2245,7 @@ void CIncEulerSolver::SetRangePressure(CGeometry *geometry, CSolver **solver_con
 
     BEGIN_SU2_OMP_SAFE_GLOBAL_ACCESS {
       minP = MinP;
-      SU2_MPI::Allreduce(&minP, &MinP, 1, MPI_DOUBLE, MPI_MAX, SU2_MPI::GetComm());
+      SU2_MPI::Allreduce(&minP, &MinP, 1, MPI_DOUBLE, MPI_MIN, SU2_MPI::GetComm());
       maxP = MaxP;
       SU2_MPI::Allreduce(&maxP, &MaxP, 1, MPI_DOUBLE, MPI_MAX, SU2_MPI::GetComm());
 
diff --git a/TestCases/hybrid_regression.py b/TestCases/hybrid_regression.py
index c545a84a8eba..60e5c01fbc7f 100644
--- a/TestCases/hybrid_regression.py
+++ b/TestCases/hybrid_regression.py
@@ -265,7 +265,7 @@ def main():
     axi_rans_air_nozzle_restart.cfg_dir   = "axisymmetric_rans/air_nozzle"
     axi_rans_air_nozzle_restart.cfg_file  = "air_nozzle_restart.cfg"
     axi_rans_air_nozzle_restart.test_iter = 10
-    axi_rans_air_nozzle_restart.test_vals = [-5.068875, 0.504406, -2.840213, 1.705449, 0.000000]
+    axi_rans_air_nozzle_restart.test_vals = [-5.069624, 0.505288, -2.816368, 1.728267, 0.000000]
     axi_rans_air_nozzle_restart.test_vals_aarch64 = [-14.140441, -9.154674, -10.886121, -5.806594, 0.000000]
     test_list.append(axi_rans_air_nozzle_restart)
 
@@ -340,7 +340,7 @@ def main():
     harmonic_balance.cfg_dir   = "harmonic_balance"
     harmonic_balance.cfg_file  = "HB.cfg"
     harmonic_balance.test_iter = 25
-    harmonic_balance.test_vals = [-1.559187, 0.829574, 0.931511, 3.954440]
+    harmonic_balance.test_vals = [-1.556450, 0.840516, 0.943600, 3.957290]
     test_list.append(harmonic_balance)
 
     # Turbulent pitching NACA 64a010 airfoil
@@ -510,7 +510,7 @@ def main():
     sine_gust.cfg_dir   = "gust"
     sine_gust.cfg_file  = "inv_gust_NACA0012.cfg"
     sine_gust.test_iter = 5
-    sine_gust.test_vals = [-1.977498, 3.481817, -0.010724, -0.005369]
+    sine_gust.test_vals = [-1.977498, 3.481833, -0.012152, -0.004132]
     sine_gust.unsteady  = True
     test_list.append(sine_gust)
 
@@ -519,7 +519,7 @@ def main():
     cosine_gust.cfg_dir   = "gust"
     cosine_gust.cfg_file  = "cosine_gust_zdir.cfg"
     cosine_gust.test_iter = 79
-    cosine_gust.test_vals = [-2.418805, 0.002212, -0.001249, 0.000439, -0.000592]
+    cosine_gust.test_vals = [-2.418805, 0.003015, -0.001242, 0.000500, -0.000592]
     cosine_gust.unsteady  = True
     cosine_gust.enabled_with_tsan = False
     test_list.append(cosine_gust)
@@ -529,7 +529,7 @@ def main():
     gust_mesh_defo.cfg_dir   = "gust"
     gust_mesh_defo.cfg_file  = "gust_with_mesh_deformation.cfg"
     gust_mesh_defo.test_iter = 6
-    gust_mesh_defo.test_vals = [-1.844761, 0.001290, -0.000316]
+    gust_mesh_defo.test_vals = [-1.844761, 0.001449, -0.000344]
     gust_mesh_defo.unsteady  = True
     gust_mesh_defo.enabled_with_tsan = False
     test_list.append(gust_mesh_defo)
@@ -539,7 +539,7 @@ def main():
     aeroelastic.cfg_dir   = "aeroelastic"
     aeroelastic.cfg_file  = "aeroelastic_NACA64A010.cfg"
     aeroelastic.test_iter = 2
-    aeroelastic.test_vals = [-1.876632, 4.021071, 0.080087, 0.027953, -0.001641, -0.000127, -1.052229]
+    aeroelastic.test_vals = [-1.876631, 4.021081, 0.082286, 0.026921, -0.001638, -0.000130, -1.019527]
     aeroelastic.unsteady  = True
     aeroelastic.enabled_on_cpu_arch = ["x86_64"] # Requires AVX-capable architecture
     test_list.append(aeroelastic)
@@ -558,7 +558,7 @@ def main():
     unst_inc_turb_naca0015_sa.cfg_dir   = "unsteady/pitching_naca0015_rans_inc"
     unst_inc_turb_naca0015_sa.cfg_file  = "config_incomp_turb_sa.cfg"
     unst_inc_turb_naca0015_sa.test_iter = 1
-    unst_inc_turb_naca0015_sa.test_vals = [-3.008640, -6.889022, 1.435181, 0.433567]
+    unst_inc_turb_naca0015_sa.test_vals = [-3.007258, -6.887901, 1.435511, 0.433955]
     unst_inc_turb_naca0015_sa.unsteady  = True
     test_list.append(unst_inc_turb_naca0015_sa)
 
@@ -567,7 +567,7 @@ def main():
     unst_deforming_naca0012.cfg_dir   = "disc_adj_euler/naca0012_pitching_def"
     unst_deforming_naca0012.cfg_file  = "inv_NACA0012_pitching_deform.cfg"
     unst_deforming_naca0012.test_iter = 5
-    unst_deforming_naca0012.test_vals = [-3.666867, -3.796806, -3.719022, -3.150288]
+    unst_deforming_naca0012.test_vals = [-3.666017, -3.794928, -3.717606, -3.149231]
     unst_deforming_naca0012.unsteady  = True
     unst_deforming_naca0012.enabled_with_tsan = False
     test_list.append(unst_deforming_naca0012)
@@ -640,7 +640,7 @@ def main():
     uniform_flow.cfg_dir   = "sliding_interface/uniform_flow"
     uniform_flow.cfg_file  = "uniform_NN.cfg"
     uniform_flow.test_iter = 5
-    uniform_flow.test_vals = [5.000000, 0.000000, -0.195002, -10.624453]
+    uniform_flow.test_vals = [5.000000, 0.000000, -0.193990, -10.624443]
     uniform_flow.unsteady  = True
     uniform_flow.multizone = True
     test_list.append(uniform_flow)
@@ -659,7 +659,7 @@ def main():
     channel_2D.cfg_dir   = "sliding_interface/channel_2D"
     channel_2D.cfg_file  = "channel_2D_WA.cfg"
     channel_2D.test_iter = 2
-    channel_2D.test_vals = [2.000000, 0.000000, 0.466339, 0.346858, 0.397486]
+    channel_2D.test_vals = [2.000000, 0.000000, 0.399109, 0.350302, 0.405582]
     channel_2D.unsteady  = True
     channel_2D.multizone = True
     test_list.append(channel_2D)
@@ -681,7 +681,7 @@ def main():
     pipe.cfg_dir   = "sliding_interface/pipe"
     pipe.cfg_file  = "pipe_NN.cfg"
     pipe.test_iter = 2
-    pipe.test_vals = [0.070225, 0.528891, 0.684177, 0.970350, 1.065268]
+    pipe.test_vals = [0.024970, 0.452310, 0.562858, 0.974611, 1.031083]
     pipe.unsteady  = True
     pipe.multizone = True
     test_list.append(pipe)
@@ -691,7 +691,7 @@ def main():
     rotating_cylinders.cfg_dir   = "sliding_interface/rotating_cylinders"
     rotating_cylinders.cfg_file  = "rot_cylinders_WA.cfg"
     rotating_cylinders.test_iter = 3
-    rotating_cylinders.test_vals = [3.000000, 0.000000, 0.686432, 1.159463, 1.151135]
+    rotating_cylinders.test_vals = [3.000000, 0.000000, 0.686419, 1.152596, 1.145896]
     rotating_cylinders.unsteady  = True
     rotating_cylinders.multizone  = True
     test_list.append(rotating_cylinders)
@@ -701,7 +701,7 @@ def main():
     supersonic_vortex_shedding.cfg_dir   = "sliding_interface/supersonic_vortex_shedding"
     supersonic_vortex_shedding.cfg_file  = "sup_vor_shed_WA.cfg"
     supersonic_vortex_shedding.test_iter = 5
-    supersonic_vortex_shedding.test_vals = [5.000000, 0.000000, 1.126036, 1.081765]
+    supersonic_vortex_shedding.test_vals = [5.000000, 0.000000, 1.128138, 1.064252]
     supersonic_vortex_shedding.unsteady  = True
     supersonic_vortex_shedding.multizone  = True
     test_list.append(supersonic_vortex_shedding)
@@ -752,7 +752,7 @@ def main():
     fsi2d.cfg_dir   = "fea_fsi/WallChannel_2d"
     fsi2d.cfg_file  = "configFSI.cfg"
     fsi2d.test_iter = 4
-    fsi2d.test_vals = [4.000000, 0.000000, -3.726029, -4.277532]
+    fsi2d.test_vals = [4.000000, 0.000000, -3.725535, -4.278762]
     fsi2d.multizone= True
     fsi2d.unsteady = True
     fsi2d.enabled_with_tsan = False
@@ -763,7 +763,7 @@ def main():
     dyn_fsi.cfg_dir   = "fea_fsi/dyn_fsi"
     dyn_fsi.cfg_file  = "config.cfg"
     dyn_fsi.test_iter = 4
-    dyn_fsi.test_vals = [-4.330727, -4.152808, 0.000000, 102.000000]
+    dyn_fsi.test_vals = [-4.328791, -4.152991, 0.000000, 103.000000]
     dyn_fsi.multizone = True
     dyn_fsi.unsteady  = True
     test_list.append(dyn_fsi)
diff --git a/TestCases/hybrid_regression_AD.py b/TestCases/hybrid_regression_AD.py
index b2086915f9df..6a693812d903 100644
--- a/TestCases/hybrid_regression_AD.py
+++ b/TestCases/hybrid_regression_AD.py
@@ -188,7 +188,7 @@ def main():
     discadj_pitchingNACA0012.cfg_dir   = "disc_adj_euler/naca0012_pitching"
     discadj_pitchingNACA0012.cfg_file  = "inv_NACA0012_pitching.cfg"
     discadj_pitchingNACA0012.test_iter = 4
-    discadj_pitchingNACA0012.test_vals = [-1.041975, -1.513011, -0.006202, 0.000012]
+    discadj_pitchingNACA0012.test_vals = [-1.128086, -1.582652, -0.006613, 0.000011]
     discadj_pitchingNACA0012.tol = 0.01
     discadj_pitchingNACA0012.unsteady  = True
     discadj_pitchingNACA0012.enabled_with_tsan = False
diff --git a/TestCases/parallel_regression.py b/TestCases/parallel_regression.py
index 64c9c7b2195d..855c7d2f8cc8 100755
--- a/TestCases/parallel_regression.py
+++ b/TestCases/parallel_regression.py
@@ -660,7 +660,7 @@ def main():
     axi_rans_air_nozzle_restart.cfg_dir   = "axisymmetric_rans/air_nozzle"
     axi_rans_air_nozzle_restart.cfg_file  = "air_nozzle_restart.cfg"
     axi_rans_air_nozzle_restart.test_iter = 10
-    axi_rans_air_nozzle_restart.test_vals = [-5.144760, 0.426779, -2.900513, 1.606061, 0.000000]
+    axi_rans_air_nozzle_restart.test_vals = [-5.130654, 0.443770, -2.934019, 1.603840, 0.000000]
     axi_rans_air_nozzle_restart.test_vals_aarch64 = [-14.143310, -9.163287, -10.858232, -5.787715, 0.000000]
     axi_rans_air_nozzle_restart.tol       = 0.0001
     test_list.append(axi_rans_air_nozzle_restart)
@@ -713,7 +713,7 @@ def main():
     inc_nozzle.cfg_dir   = "incomp_euler/nozzle"
     inc_nozzle.cfg_file  = "inv_nozzle.cfg"
     inc_nozzle.test_iter = 20
-    inc_nozzle.test_vals = [-6.051005, -5.361765, 0.006999, 0.127295]
+    inc_nozzle.test_vals = [-6.357399, -5.618580, 0.004548, 0.126243]
     test_list.append(inc_nozzle)
 
     # Laminar wall mounted cylinder, Euler walls, cylinder wall diagonally split
@@ -1154,7 +1154,7 @@ def main():
     harmonic_balance.cfg_dir   = "harmonic_balance"
     harmonic_balance.cfg_file  = "HB.cfg"
     harmonic_balance.test_iter = 25
-    harmonic_balance.test_vals = [-1.559187, 0.829575, 0.931512, 3.954440]
+    harmonic_balance.test_vals = [-1.556449, 0.840517, 0.943602, 3.957291]
     test_list.append(harmonic_balance)
 
     # Turbulent pitching NACA 64a010 airfoil
@@ -1175,7 +1175,7 @@ def main():
     rot_naca0012.cfg_dir   = "rotating/naca0012"
     rot_naca0012.cfg_file  = "rot_NACA0012.cfg"
     rot_naca0012.test_iter = 25
-    rot_naca0012.test_vals = [-0.023360, 0.082942, 2.000000, -1.182142]
+    rot_naca0012.test_vals = [-0.028757, 0.058709, 2.000000, -1.147370]
     test_list.append(rot_naca0012)
 
     # Lid-driven cavity
@@ -1225,7 +1225,7 @@ def main():
     sine_gust.cfg_dir   = "gust"
     sine_gust.cfg_file  = "inv_gust_NACA0012.cfg"
     sine_gust.test_iter = 5
-    sine_gust.test_vals = [-1.977498, 3.481817, -0.011360, -0.005241]
+    sine_gust.test_vals = [-1.977498, 3.481833, -0.013388, -0.003910]
     sine_gust.unsteady  = True
     test_list.append(sine_gust)
 
@@ -1254,7 +1254,7 @@ def main():
     unst_inc_turb_naca0015_sa.cfg_dir   = "unsteady/pitching_naca0015_rans_inc"
     unst_inc_turb_naca0015_sa.cfg_file  = "config_incomp_turb_sa.cfg"
     unst_inc_turb_naca0015_sa.test_iter = 1
-    unst_inc_turb_naca0015_sa.test_vals = [-3.004011, -6.876257, 1.487888, 0.421869]
+    unst_inc_turb_naca0015_sa.test_vals = [-3.002783, -6.875096, 1.487783, 0.421690]
     unst_inc_turb_naca0015_sa.unsteady  = True
     test_list.append(unst_inc_turb_naca0015_sa)
 
@@ -1369,7 +1369,7 @@ def main():
     uniform_flow.cfg_dir   = "sliding_interface/uniform_flow"
     uniform_flow.cfg_file  = "uniform_NN.cfg"
     uniform_flow.test_iter = 5
-    uniform_flow.test_vals = [5.000000, 0.000000, -0.195001, -10.624448]
+    uniform_flow.test_vals = [5.000000, 0.000000, -0.193989, -10.624451]
     uniform_flow.unsteady  = True
     uniform_flow.multizone = True
     test_list.append(uniform_flow)
@@ -1387,7 +1387,7 @@ def main():
     channel_2D.cfg_dir   = "sliding_interface/channel_2D"
     channel_2D.cfg_file  = "channel_2D_WA.cfg"
     channel_2D.test_iter = 2
-    channel_2D.test_vals = [2.000000, 0.000000, 0.466408, 0.346806, 0.397328]
+    channel_2D.test_vals = [2.000000, 0.000000, 0.399141, 0.350260, 0.405570]
     channel_2D.timeout   = 100
     channel_2D.unsteady  = True
     channel_2D.multizone = True
@@ -1409,7 +1409,7 @@ def main():
     pipe.cfg_dir   = "sliding_interface/pipe"
     pipe.cfg_file  = "pipe_NN.cfg"
     pipe.test_iter = 2
-    pipe.test_vals = [0.070225, 0.528893, 0.684173, 0.970355, 1.065262]
+    pipe.test_vals = [0.024970, 0.452312, 0.562852, 0.974613, 1.031075]
     pipe.unsteady  = True
     pipe.multizone = True
     test_list.append(pipe)
@@ -1419,7 +1419,7 @@ def main():
     rotating_cylinders.cfg_dir   = "sliding_interface/rotating_cylinders"
     rotating_cylinders.cfg_file  = "rot_cylinders_WA.cfg"
     rotating_cylinders.test_iter = 3
-    rotating_cylinders.test_vals = [3.000000, 0.000000, 0.686431, 1.159459, 1.151139]
+    rotating_cylinders.test_vals = [3.000000, 0.000000, 0.686417, 1.152593, 1.145900]
     rotating_cylinders.unsteady  = True
     rotating_cylinders.multizone  = True
     test_list.append(rotating_cylinders)
@@ -1429,7 +1429,7 @@ def main():
     supersonic_vortex_shedding.cfg_dir   = "sliding_interface/supersonic_vortex_shedding"
     supersonic_vortex_shedding.cfg_file  = "sup_vor_shed_WA.cfg"
     supersonic_vortex_shedding.test_iter = 5
-    supersonic_vortex_shedding.test_vals = [5.000000, 0.000000, 1.126036, 1.081775]
+    supersonic_vortex_shedding.test_vals = [5.000000, 0.000000, 1.128138, 1.064262]
     supersonic_vortex_shedding.unsteady  = True
     supersonic_vortex_shedding.multizone  = True
     test_list.append(supersonic_vortex_shedding)
@@ -1538,7 +1538,7 @@ def main():
     fsi2d.cfg_dir   = "fea_fsi/WallChannel_2d"
     fsi2d.cfg_file  = "configFSI.cfg"
     fsi2d.test_iter = 4
-    fsi2d.test_vals = [4, 0, -3.726013, -4.277529]
+    fsi2d.test_vals = [4.000000, 0.000000, -3.725519, -4.278760]
     fsi2d.command   = TestCase.Command(exec = "parallel_computation_fsi.py", param = "-f")
     fsi2d.multizone= True
     fsi2d.unsteady = True
@@ -1549,7 +1549,7 @@ def main():
     dyn_fsi.cfg_dir   = "fea_fsi/dyn_fsi"
     dyn_fsi.cfg_file  = "config.cfg"
     dyn_fsi.test_iter = 4
-    dyn_fsi.test_vals = [-4.330741, -4.152826, 0.000000, 75.000000]
+    dyn_fsi.test_vals = [-4.328805, -4.153009, 0.000000, 75.000000]
     dyn_fsi.multizone = True
     dyn_fsi.unsteady  = True
     test_list.append(dyn_fsi)
@@ -1660,7 +1660,7 @@ def main():
     pywrapper_aeroelastic.cfg_dir   = "aeroelastic"
     pywrapper_aeroelastic.cfg_file  = "aeroelastic_NACA64A010.cfg"
     pywrapper_aeroelastic.test_iter = 2
-    pywrapper_aeroelastic.test_vals = [-1.876636, 4.021064, 0.080986, 0.027812, -0.001643, -0.000126, -0.965663]
+    pywrapper_aeroelastic.test_vals = [-1.876625, 4.021088, 0.083361, 0.026765, -0.001638, -0.000130, -0.947102]
     pywrapper_aeroelastic.command   = TestCase.Command("mpirun -np 2", "SU2_CFD.py", "--parallel -f")
     pywrapper_aeroelastic.unsteady  = True
     test_list.append(pywrapper_aeroelastic)
@@ -1679,7 +1679,7 @@ def main():
     pywrapper_fsi2d.cfg_dir   = "fea_fsi/WallChannel_2d"
     pywrapper_fsi2d.cfg_file  = "configFSI.cfg"
     pywrapper_fsi2d.test_iter = 4
-    pywrapper_fsi2d.test_vals = [4.000000, 0.000000, -3.726013, -4.277529]
+    pywrapper_fsi2d.test_vals = [4.000000, 0.000000, -3.725519, -4.278760]
     pywrapper_fsi2d.command   = TestCase.Command("mpirun -np 2", "SU2_CFD.py", "--nZone 2 --fsi True --parallel -f")
     pywrapper_fsi2d.unsteady  = True
     pywrapper_fsi2d.multizone = True
@@ -1690,7 +1690,7 @@ def main():
     pywrapper_unsteadyFSI.cfg_dir = "py_wrapper/dyn_fsi"
     pywrapper_unsteadyFSI.cfg_file = "config.cfg"
     pywrapper_unsteadyFSI.test_iter = 4
-    pywrapper_unsteadyFSI.test_vals = [0, 31, 5, 47, -1.756677, -2.828286, -7.638545, -6.863930, 0.000156]
+    pywrapper_unsteadyFSI.test_vals = [0.000000, 30.000000, 5.000000, 48.000000, -1.756676, -2.828284, -7.625018, -6.863923, 0.000156]
     pywrapper_unsteadyFSI.command = TestCase.Command("mpirun -np 2", "python", "run.py")
     pywrapper_unsteadyFSI.unsteady = True
     pywrapper_unsteadyFSI.multizone = True
@@ -1721,7 +1721,7 @@ def main():
     pywrapper_deformingBump.cfg_dir = "py_wrapper/deforming_bump_in_channel"
     pywrapper_deformingBump.cfg_file = "config.cfg"
     pywrapper_deformingBump.test_iter = 1
-    pywrapper_deformingBump.test_vals = [0.500000, 0.000000, -2.556310, -1.270841, -2.350592, 2.606850, 8.015020, -0.301652]
+    pywrapper_deformingBump.test_vals = [0.500000, 0.000000, -2.556310, -1.270836, -2.351316, 2.905858, 9.953163, -0.946408]
     pywrapper_deformingBump.command = TestCase.Command("mpirun -np 2", "python", "run.py")
     pywrapper_deformingBump.unsteady = True
     test_list.append(pywrapper_deformingBump)
diff --git a/TestCases/parallel_regression_AD.py b/TestCases/parallel_regression_AD.py
index 99e31c3124a0..ebd1193f4b0c 100644
--- a/TestCases/parallel_regression_AD.py
+++ b/TestCases/parallel_regression_AD.py
@@ -150,7 +150,7 @@ def main():
     discadj_axisymmetric_rans_nozzle.cfg_dir    = "axisymmetric_rans/air_nozzle"
     discadj_axisymmetric_rans_nozzle.cfg_file   = "air_nozzle_restart.cfg"
     discadj_axisymmetric_rans_nozzle.test_iter  = 10
-    discadj_axisymmetric_rans_nozzle.test_vals = [9.742490, 5.136585, 7.107742, 2.488366]
+    discadj_axisymmetric_rans_nozzle.test_vals = [9.738854, 5.135195, 7.108430, 2.487299]
     discadj_axisymmetric_rans_nozzle.no_restart = True
     test_list.append(discadj_axisymmetric_rans_nozzle)
 
@@ -217,7 +217,7 @@ def main():
     discadj_pitchingNACA0012.cfg_dir   = "disc_adj_euler/naca0012_pitching"
     discadj_pitchingNACA0012.cfg_file  = "inv_NACA0012_pitching.cfg"
     discadj_pitchingNACA0012.test_iter = 4
-    discadj_pitchingNACA0012.test_vals = [-1.039711, -1.508402, -0.006059, 0.000012]
+    discadj_pitchingNACA0012.test_vals = [-1.115587, -1.573317, -0.006474, 0.000010]
     discadj_pitchingNACA0012.tol = 0.01
     discadj_pitchingNACA0012.unsteady  = True
     test_list.append(discadj_pitchingNACA0012)
diff --git a/TestCases/serial_regression.py b/TestCases/serial_regression.py
index dccf3b6a8a4f..eb6c08d5ffc4 100755
--- a/TestCases/serial_regression.py
+++ b/TestCases/serial_regression.py
@@ -375,7 +375,7 @@ def main():
     axi_rans_air_nozzle_restart.cfg_dir   = "axisymmetric_rans/air_nozzle"
     axi_rans_air_nozzle_restart.cfg_file  = "air_nozzle_restart.cfg"
     axi_rans_air_nozzle_restart.test_iter = 10
-    axi_rans_air_nozzle_restart.test_vals = [-5.135350, 0.436347, -2.878865, 1.663300, 0.000000]
+    axi_rans_air_nozzle_restart.test_vals = [-5.121517, 0.452804, -2.910659, 1.662143, 0.000000]
     axi_rans_air_nozzle_restart.test_vals_aarch64 = [-14.143715, -9.170705, -10.848554, -5.776746, 0.000000]
     axi_rans_air_nozzle_restart.tol       = 0.0001
     test_list.append(axi_rans_air_nozzle_restart)
@@ -385,7 +385,7 @@ def main():
     axi_rans_air_nozzle_species.cfg_dir   = "axisymmetric_rans/air_nozzle"
     axi_rans_air_nozzle_species.cfg_file  = "air_nozzle_species.cfg"
     axi_rans_air_nozzle_species.test_iter = 10
-    axi_rans_air_nozzle_species.test_vals = [-2.243181, 3.273075, -1.958732, 5.883857, -2.469621, 0.000000]
+    axi_rans_air_nozzle_species.test_vals = [-2.148980, 3.342062, -1.939182, 5.883411, -2.443141, 0.000000]
     axi_rans_air_nozzle_species.tol       = 0.0001
     test_list.append(axi_rans_air_nozzle_species)
 
@@ -798,7 +798,7 @@ def main():
     harmonic_balance.cfg_dir   = "harmonic_balance"
     harmonic_balance.cfg_file  = "HB.cfg"
     harmonic_balance.test_iter = 25
-    harmonic_balance.test_vals = [-1.559187, 0.829575, 0.931512, 3.954440]
+    harmonic_balance.test_vals = [-1.556449, 0.840517, 0.943602, 3.957291]
     test_list.append(harmonic_balance)
 
     # Turbulent pitching NACA 64a010 airfoil
@@ -818,7 +818,7 @@ def main():
     rot_naca0012.cfg_dir   = "rotating/naca0012"
     rot_naca0012.cfg_file  = "rot_NACA0012.cfg"
     rot_naca0012.test_iter = 25
-    rot_naca0012.test_vals = [-1.240352, 4.268727, -0.022164, 0.083891, 2.000000, -1.253214]
+    rot_naca0012.test_vals = [-1.281757, 4.222928, -0.027656, 0.059256, 2.000000, -1.218300]
     test_list.append(rot_naca0012)
 
     # Lid-driven cavity
@@ -855,7 +855,7 @@ def main():
     sine_gust.cfg_dir   = "gust"
     sine_gust.cfg_file  = "inv_gust_NACA0012.cfg"
     sine_gust.test_iter = 5
-    sine_gust.test_vals = [-1.977498, 3.481817, -0.010871, -0.005083]
+    sine_gust.test_vals = [-1.977498, 3.481833, -0.012952, -0.003623]
     sine_gust.unsteady  = True
     test_list.append(sine_gust)
 
@@ -864,7 +864,7 @@ def main():
     aeroelastic.cfg_dir   = "aeroelastic"
     aeroelastic.cfg_file  = "aeroelastic_NACA64A010.cfg"
     aeroelastic.test_iter = 2
-    aeroelastic.test_vals = [-1.876632, 4.021072, 0.080240, 0.027935, -0.001641, -0.000127, -1.133898]
+    aeroelastic.test_vals = [-1.876630, 4.021083, 0.082484, 0.026820, -0.001638, -0.000129, -1.081133]
     aeroelastic.unsteady  = True
     test_list.append(aeroelastic)
 
@@ -882,7 +882,7 @@ def main():
     unst_inc_turb_naca0015_sa.cfg_dir   = "unsteady/pitching_naca0015_rans_inc"
     unst_inc_turb_naca0015_sa.cfg_file  = "config_incomp_turb_sa.cfg"
     unst_inc_turb_naca0015_sa.test_iter = 1
-    unst_inc_turb_naca0015_sa.test_vals = [-3.007635, -6.879816, 1.445300, 0.419281]
+    unst_inc_turb_naca0015_sa.test_vals = [-3.005435, -6.878657, 1.445291, 0.420526]
     unst_inc_turb_naca0015_sa.unsteady  = True
     test_list.append(unst_inc_turb_naca0015_sa)
     # unsteady pitching NACA64A010, RANS, SA
@@ -898,7 +898,7 @@ def main():
     unst_pitching_naca64a010_euler.cfg_dir   = "unsteady/pitching_naca64a010_euler"
     unst_pitching_naca64a010_euler.cfg_file  = "pitching_NACA64A010.cfg"
     unst_pitching_naca64a010_euler.test_iter = 2
-    unst_pitching_naca64a010_euler.test_vals = [-1.186839, 4.280301, -0.039724, 0.000927]
+    unst_pitching_naca64a010_euler.test_vals = [-1.186839, 4.280671, -0.026619, 0.000937]
     unst_pitching_naca64a010_euler.unsteady  = True
     test_list.append(unst_pitching_naca64a010_euler)
     # unsteady plunging NACA0012, Laminar NS
@@ -915,7 +915,7 @@ def main():
     unst_deforming_naca0012.cfg_dir   = "disc_adj_euler/naca0012_pitching_def"
     unst_deforming_naca0012.cfg_file  = "inv_NACA0012_pitching_deform.cfg"
     unst_deforming_naca0012.test_iter = 5
-    unst_deforming_naca0012.test_vals = [-3.667005, -3.796733, -3.719112, -3.150401]
+    unst_deforming_naca0012.test_vals = [-3.666074, -3.794974, -3.717659, -3.149285]
     unst_deforming_naca0012.unsteady  = True
     test_list.append(unst_deforming_naca0012)
 
@@ -1012,7 +1012,7 @@ def main():
     uniform_flow.cfg_dir   = "sliding_interface/uniform_flow"
     uniform_flow.cfg_file  = "uniform_NN.cfg"
     uniform_flow.test_iter = 2
-    uniform_flow.test_vals = [2.000000, 0.000000, -0.230639, -13.250208]
+    uniform_flow.test_vals = [2.000000, 0.000000, -0.230680, -13.251592]
     uniform_flow.test_vals_aarch64 = [2.000000, 0.000000, -0.230641, -13.249000]
     uniform_flow.tol       = 0.000001
     uniform_flow.unsteady  = True
@@ -1033,7 +1033,7 @@ def main():
     channel_2D.cfg_dir   = "sliding_interface/channel_2D"
     channel_2D.cfg_file  = "channel_2D_WA.cfg"
     channel_2D.test_iter = 2
-    channel_2D.test_vals = [2.000000, 0.000000, 0.466389, 0.346857, 0.397472]
+    channel_2D.test_vals = [2.000000, 0.000000, 0.399122, 0.350309, 0.405616]
     channel_2D.timeout   = 100
     channel_2D.unsteady  = True
     channel_2D.multizone = True
@@ -1056,7 +1056,7 @@ def main():
     pipe.cfg_dir   = "sliding_interface/pipe"
     pipe.cfg_file  = "pipe_NN.cfg"
     pipe.test_iter = 2
-    pipe.test_vals = [0.528897, 0.684175, 0.970353, 1.065262]
+    pipe.test_vals = [0.452310, 0.562852, 0.974610, 1.031075]
     pipe.unsteady  = True
     pipe.multizone = True
     test_list.append(pipe)
@@ -1066,7 +1066,7 @@ def main():
     rotating_cylinders.cfg_dir   = "sliding_interface/rotating_cylinders"
     rotating_cylinders.cfg_file  = "rot_cylinders_WA.cfg"
     rotating_cylinders.test_iter = 3
-    rotating_cylinders.test_vals = [3.000000, 0.000000, 0.686432, 1.159461, 1.151134]
+    rotating_cylinders.test_vals = [3.000000, 0.000000, 0.686419, 1.152594, 1.145895]
     rotating_cylinders.unsteady  = True
     rotating_cylinders.multizone = True
     test_list.append(rotating_cylinders)
@@ -1076,7 +1076,7 @@ def main():
     supersonic_vortex_shedding.cfg_dir   = "sliding_interface/supersonic_vortex_shedding"
     supersonic_vortex_shedding.cfg_file  = "sup_vor_shed_WA.cfg"
     supersonic_vortex_shedding.test_iter = 5
-    supersonic_vortex_shedding.test_vals = [5.000000, 0.000000, 1.126036, 1.081775]
+    supersonic_vortex_shedding.test_vals = [5.000000, 0.000000, 1.128138, 1.064261]
     supersonic_vortex_shedding.unsteady  = True
     supersonic_vortex_shedding.multizone = True
     test_list.append(supersonic_vortex_shedding)
@@ -1154,7 +1154,7 @@ def main():
     fsi2d.cfg_dir   = "fea_fsi/WallChannel_2d"
     fsi2d.cfg_file  = "configFSI.cfg"
     fsi2d.test_iter = 4
-    fsi2d.test_vals = [4, 0, -3.726015, -4.277528]
+    fsi2d.test_vals = [4.000000, 0.000000, -3.725521, -4.278759]
     fsi2d.multizone = True
     fsi2d.unsteady  = True
     test_list.append(fsi2d)
@@ -1173,7 +1173,7 @@ def main():
     dyn_fsi.cfg_dir   = "fea_fsi/dyn_fsi"
     dyn_fsi.cfg_file  = "config.cfg"
     dyn_fsi.test_iter = 4
-    dyn_fsi.test_vals = [-4.330728, -4.152820, 0.000000, 85.000000]
+    dyn_fsi.test_vals = [-4.328792, -4.153004, 0.000000, 86.000000]
     dyn_fsi.multizone = True
     dyn_fsi.unsteady  = True
     test_list.append(dyn_fsi)
@@ -1716,7 +1716,7 @@ def main():
     pywrapper_fsi2d.cfg_dir   = "fea_fsi/WallChannel_2d"
     pywrapper_fsi2d.cfg_file  = "configFSI.cfg"
     pywrapper_fsi2d.test_iter = 4
-    pywrapper_fsi2d.test_vals = [4, 0, -3.726015, -4.277528]
+    pywrapper_fsi2d.test_vals = [4.000000, 0.000000, -3.725521, -4.278759]
     pywrapper_fsi2d.command   =  TestCase.Command(exec = "SU2_CFD.py", param = "--nZone 2 --fsi True -f")
     pywrapper_fsi2d.unsteady  = True
     pywrapper_fsi2d.multizone = True
diff --git a/TestCases/serial_regression_AD.py b/TestCases/serial_regression_AD.py
index 288f8a47efd3..83e768e954e8 100644
--- a/TestCases/serial_regression_AD.py
+++ b/TestCases/serial_regression_AD.py
@@ -158,7 +158,7 @@ def main():
     discadj_pitchingNACA0012.cfg_dir   = "disc_adj_euler/naca0012_pitching"
     discadj_pitchingNACA0012.cfg_file  = "inv_NACA0012_pitching.cfg"
     discadj_pitchingNACA0012.test_iter = 4
-    discadj_pitchingNACA0012.test_vals = [-1.057526, -1.512609, -0.004842, 0.000013]
+    discadj_pitchingNACA0012.test_vals = [-1.138369, -1.575068, -0.005301, 0.000003]
     discadj_pitchingNACA0012.unsteady  = True
     test_list.append(discadj_pitchingNACA0012)
 
@@ -167,7 +167,7 @@ def main():
     unst_deforming_naca0012.cfg_dir   = "disc_adj_euler/naca0012_pitching_def"
     unst_deforming_naca0012.cfg_file  = "inv_NACA0012_pitching_deform_ad.cfg"
     unst_deforming_naca0012.test_iter = 4
-    unst_deforming_naca0012.test_vals = [-1.874283, -1.768562, 4274.000000, 0.000003]
+    unst_deforming_naca0012.test_vals = [-1.875062, -1.766683, 4307.900000, 0.000003]
     unst_deforming_naca0012.unsteady  = True
     test_list.append(unst_deforming_naca0012)
 
diff --git a/UnitTests/Common/CConfig_tests.cpp b/UnitTests/Common/CConfig_tests.cpp
index a499a075fa1d..2dfe8b65cf70 100644
--- a/UnitTests/Common/CConfig_tests.cpp
+++ b/UnitTests/Common/CConfig_tests.cpp
@@ -88,6 +88,20 @@ TEST_CASE("INIT_OPTION_INC defaults", "[Config]") {
   CHECK(GetInitOptionInc(ideal_gas_options + "INIT_OPTION_INC= DENSITY_INIT\n") == INIT_OPTION_INC::DENSITY_INIT);
 }
 
+TEST_CASE("Supported multigrid levels are retained", "[Config][MultigridLimits]") {
+  for (const auto levels : {0, 10}) {
+    std::stringstream options(base_options + "MGLEVEL= " + std::to_string(levels) + "\n");
+    CConfig config(options, SU2_COMPONENT::SU2_CFD, false);
+    CHECK(config.GetnMGLevels() == levels);
+  }
+}
+
+TEST_CASE("Excessive multigrid levels are rejected", "[.MultigridLimitError]") {
+  std::stringstream options(base_options + "MGLEVEL= 11\n");
+  CConfig config(options, SU2_COMPONENT::SU2_CFD, false);
+  FAIL("An unsupported multigrid hierarchy was accepted.");
+}
+
 TEST_CASE("HLLC and AUSM SLAU Jacobian options are independent", "[Config][HLLCOptions]") {
   for (const auto accurate : {false, true}) {
     for (const auto exactHLLC : {false, true}) {
diff --git a/UnitTests/Common/check_multigrid_limits.py b/UnitTests/Common/check_multigrid_limits.py
new file mode 100644
index 000000000000..5070fdaed3cc
--- /dev/null
+++ b/UnitTests/Common/check_multigrid_limits.py
@@ -0,0 +1,16 @@
+#!/usr/bin/env python3
+"""Check the fatal configuration diagnostic in its own process."""
+import subprocess
+import sys
+
+result = subprocess.run(
+    [sys.argv[1], "[.MultigridLimitError]"], capture_output=True, text=True, timeout=30
+)
+expected = "MGLEVEL must be between 0 and 10"
+if result.returncode == 0 or expected not in result.stdout + result.stderr:
+    raise RuntimeError(
+        "Unsupported MGLEVEL did not fail with its diagnostic:\n"
+        + result.stdout
+        + result.stderr
+    )
+print("PASS: MGLEVEL=11 is rejected with its supported range")
diff --git a/UnitTests/SU2_CFD/solvers/moving_wall_tests.cpp b/UnitTests/SU2_CFD/solvers/moving_wall_tests.cpp
new file mode 100644
index 000000000000..1baf0c271428
--- /dev/null
+++ b/UnitTests/SU2_CFD/solvers/moving_wall_tests.cpp
@@ -0,0 +1,182 @@
+/*!
+ * \file moving_wall_tests.cpp
+ * \brief Regression for compressible wall projection in fixed and uniformly translating frames.
+ * \version 8.5.0 "Harrier"
+ *
+ * SU2 Project Website: https://su2code.github.io
+ *
+ * The SU2 Project is maintained by the SU2 Foundation
+ * (http://su2foundation.org)
+ *
+ * Copyright 2012-2026, SU2 Contributors (cf. AUTHORS.md)
+ *
+ * SU2 is free software; you can redistribute it and/or
+ * modify it under the terms of the GNU Lesser General Public
+ * License as published by the Free Software Foundation; either
+ * version 2.1 of the License, or (at your option) any later version.
+ *
+ * SU2 is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
+ * Lesser General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public
+ * License along with SU2. If not, see <http://www.gnu.org/licenses/>.
+ */
+
+#include "catch.hpp"
+#include <array>
+#include <memory>
+#include "../../UnitQuadTestCase.hpp"
+#include "../../../SU2_CFD/include/numerics/flow/convection/hllc.hpp"
+#include "../../../SU2_CFD/include/fluid/CPengRobinson.hpp"
+
+TEST_CASE("Compressible wall projection commutes with a Galilean boost", "[HLLC-wall]") {
+  for (const auto nDim : {2u, 3u}) {
+    for (const bool realGas : {false, true}) {
+      UnitQuadTestCase fixed, moving;
+      const std::string markers =
+          nDim == 2 ? "(x_minus,x_plus,y_minus,y_plus)" : "(x_minus,x_plus,y_minus,y_plus,z_minus,z_plus)";
+      const std::string options = "SOLVER= EULER\nMESH_FORMAT= " + std::string(nDim == 2 ? "RECTANGLE" : "BOX") +
+                                  "\nMESH_BOX_SIZE= 5,5,5\nREF_DIMENSIONALIZATION= DIMENSIONAL\n"
+                                  "MACH_NUMBER= 0.1\nFREESTREAM_OPTION= DENSITY_FS\nFREESTREAM_DENSITY= 230\n"
+                                  "FREESTREAM_PRESSURE= 1500112\nGAS_CONSTANT= 35.152\n"
+                                  "CONV_NUM_METHOD_FLOW= HLLC\nTIME_DISCRE_FLOW= EULER_IMPLICIT\n"
+                                  "GRID_MOVEMENT= STEADY_TRANSLATION\nTRANSLATION_RATE= 0,0,0\nMARKER_EULER= " +
+                                  markers + "\nMARKER_MONITORING= (x_minus)\n" +
+                                  (realGas ? "FLUID_MODEL= PR_GAS\nGAMMA_VALUE= 1.0125\nCRITICAL_PRESSURE= "
+                                             "1415200\nCRITICAL_TEMPERATURE= 564.1\nACENTRIC_FACTOR= 0.529\n"
+                                           : "FLUID_MODEL= IDEAL_GAS\nGAMMA_VALUE= 1.4\n");
+      fixed.config_options = moving.config_options = options;
+      fixed.InitConfig();
+      fixed.InitGeometry();
+      fixed.InitSolver();
+      moving.InitConfig();
+      moving.InitGeometry();
+      moving.InitSolver();
+      const auto nVar = nDim + 2;
+      const su2double boost[3] = {20, -7, 3};
+      const su2double velocity[3] = {11, 4, -2};
+      su2double B[5][5] = {}, T[5][5] = {}, baseJac[5][5] = {}, boostedJac[5][5] = {};
+      su2double c2 = 0;
+      for (auto i = 0u; i < nVar; ++i) {
+        B[i][i] = T[i][i] = 1;
+        for (auto j = 0u; j < nVar; ++j) baseJac[i][j] = 1 + i + 2 * j;
+      }
+      for (auto i = 0u; i < nDim; ++i) {
+        c2 += boost[i] * boost[i];
+        B[i + 1][0] = boost[i];
+        T[i + 1][0] = -boost[i];
+        B[nVar - 1][i + 1] = boost[i];
+        T[nVar - 1][i + 1] = -boost[i];
+      }
+      B[nVar - 1][0] = T[nVar - 1][0] = 0.5 * c2;
+      for (auto i = 0u; i < nVar; ++i)
+        for (auto j = 0u; j < nVar; ++j)
+          for (auto k = 0u; k < nVar; ++k)
+            for (auto l = 0u; l < nVar; ++l) boostedJac[i][j] += B[i][k] * baseJac[k][l] * T[l][j];
+      std::unique_ptr<CFluidModel> fluid =
+          realGas ? std::unique_ptr<CFluidModel>(new CPengRobinson(1.0125, 35.152, 1415200, 564.1, 0.529))
+                  : std::unique_ptr<CFluidModel>(new CIdealGas(1.4, 35.152));
+      fluid->SetTDState_Prho(1500112, 230);
+      std::array<su2double, 5> baseState{};
+      baseState[0] = 230;
+      su2double speed2 = 0;
+      for (auto i = 0u; i < nDim; ++i) {
+        baseState[i + 1] = 230 * velocity[i];
+        speed2 += velocity[i] * velocity[i];
+      }
+      baseState[nVar - 1] = 230 * (fluid->GetStaticEnergy() + 0.5 * speed2);
+      const su2double baseResidual[5] = {0.3, -0.5, 0.2, 0.7, -0.1};
+      for (auto frame = 0u; frame < 2; ++frame) {
+        auto& test = frame ? moving : fixed;
+        auto* solver = test.solver[FLOW_SOL];
+        auto* nodes = static_cast<CEulerVariable*>(solver->GetNodes());
+        for (auto p = 0ul; p < test.geometry->GetnPoint(); ++p) {
+          su2double state[5] = {}, residual[5] = {};
+          for (auto i = 0u; i < nVar; ++i)
+            for (auto j = 0u; j < nVar; ++j) {
+              state[i] += (frame ? B[i][j] : (i == j)) * baseState[j];
+              residual[i] += (frame ? B[i][j] : (i == j)) * baseResidual[j];
+            }
+          nodes->SetSolution(p, state);
+          nodes->SetSolution_Old(p, state);
+          REQUIRE(nodes->SetPrimVar(p, fluid.get()));
+          nodes->SetSecondaryVar(p, fluid.get());
+          for (auto i = 0u; i < nDim; ++i) test.geometry->nodes->SetGridVel(p, i, frame ? boost[i] : 0);
+          for (auto i = 0u; i < nVar; ++i) {
+            solver->LinSysRes(p, i) = residual[i];
+            nodes->GetResTruncError(p)[i] = residual[i];
+          }
+          su2double block[25] = {};
+          for (auto i = 0u; i < nVar; ++i)
+            for (auto j = 0u; j < nVar; ++j) block[i * nVar + j] = frame ? boostedJac[i][j] : baseJac[i][j];
+          solver->Jacobian.SetBlock(p, p, block);
+          for (auto k = 0u; k < test.geometry->nodes->GetnPoint(p); ++k)
+            solver->Jacobian.SetBlock(p, test.geometry->nodes->GetPoint(p, k), block);
+        }
+        std::unique_ptr<CNumerics> numerics =
+            realGas ? std::unique_ptr<CNumerics>(new CUpwGeneralHLLC_Flow(nDim, nVar, test.config.get()))
+                    : std::unique_ptr<CNumerics>(new CUpwHLLC_Flow(nDim, nVar, test.config.get()));
+        for (auto marker = 0u; marker < test.config->GetnMarker_All(); ++marker)
+          if (test.config->GetMarker_All_TagBound(marker) == "x_minus")
+            solver->BC_Sym_Plane(test.geometry.get(), test.solver, numerics.get(), nullptr, test.config.get(), marker);
+      }
+      for (auto p = 0ul; p < fixed.geometry->GetnPoint(); ++p) {
+        auto* f = fixed.solver[FLOW_SOL];
+        auto* m = moving.solver[FLOW_SOL];
+        for (auto i = 0u; i < nVar; ++i) {
+          su2double state = 0, residual = 0, truncation = 0;
+          for (auto j = 0u; j < nVar; ++j) {
+            state += B[i][j] * f->GetNodes()->GetSolution_Old(p)[j];
+            residual += B[i][j] * f->LinSysRes(p, j);
+            truncation += B[i][j] * f->GetNodes()->GetResTruncError(p)[j];
+          }
+          CHECK(m->GetNodes()->GetSolution_Old(p)[i] == Approx(state).margin(1e-6));
+          CHECK(m->LinSysRes(p, i) == Approx(residual).margin(1e-8));
+          CHECK(m->GetNodes()->GetResTruncError(p)[i] == Approx(truncation).margin(1e-8));
+          for (auto block = 0u; block <= fixed.geometry->nodes->GetnPoint(p); ++block) {
+            const auto neighbor = block ? fixed.geometry->nodes->GetPoint(p, block - 1) : p;
+            for (auto j = 0u; j < nVar; ++j) {
+              su2double expected = 0;
+              for (auto k = 0u; k < nVar; ++k)
+                for (auto l = 0u; l < nVar; ++l)
+                  expected += B[i][k] * f->Jacobian.GetBlock(p, neighbor)[k * nVar + l] * T[l][j];
+              CHECK(m->Jacobian.GetBlock(p, neighbor)[i * nVar + j] == Approx(expected).margin(1e-5));
+            }
+          }
+        }
+      }
+
+      for (const bool correction : {false, true}) {
+        const su2double delta[5] = {0.3, -0.5, 0.2, 0.7, -0.1};
+        for (auto frame = 0u; frame < 2; ++frame) {
+          auto& test = frame ? moving : fixed;
+          auto* solver = test.solver[FLOW_SOL];
+          for (auto p = 0ul; p < test.geometry->GetnPoint(); ++p) {
+            su2double state[5] = {};
+            for (auto i = 0u; i < nVar; ++i)
+              for (auto j = 0u; j < nVar; ++j)
+                state[i] += (frame ? B[i][j] : (i == j)) * (correction ? delta[j] : baseState[j]);
+            if (correction)
+              solver->GetNodes()->SetSolution_Old(p, state);
+            else
+              solver->GetNodes()->SetSolution(p, state);
+          }
+          solver->MultigridProjectEulerWall(test.geometry.get(), test.config.get(), correction);
+        }
+        for (auto p = 0ul; p < fixed.geometry->GetnPoint(); ++p) {
+          const auto* f = correction ? fixed.solver[FLOW_SOL]->GetNodes()->GetSolution_Old(p)
+                                     : fixed.solver[FLOW_SOL]->GetNodes()->GetSolution(p);
+          const auto* m = correction ? moving.solver[FLOW_SOL]->GetNodes()->GetSolution_Old(p)
+                                     : moving.solver[FLOW_SOL]->GetNodes()->GetSolution(p);
+          for (auto i = 0u; i < nVar; ++i) {
+            su2double expected = 0;
+            for (auto j = 0u; j < nVar; ++j) expected += B[i][j] * f[j];
+            CHECK(m[i] == Approx(expected).margin(1e-6));
+          }
+        }
+      }
+    }
+  }
+}
diff --git a/UnitTests/meson.build b/UnitTests/meson.build
index fd15beb35653..8f7ec31afda5 100644
--- a/UnitTests/meson.build
+++ b/UnitTests/meson.build
@@ -23,6 +23,7 @@ su2_cfd_tests = files(['Common/CConfig_tests.cpp',
                        'SU2_CFD/gradients.cpp',
                        'SU2_CFD/periodic_limiters.cpp',
                        'SU2_CFD/solvers/CNEMOEulerSolver_tests.cpp',
+                       'SU2_CFD/solvers/moving_wall_tests.cpp',
                        'SU2_CFD/nemo_viscous_assembly.cpp',
                        'SU2_CFD/windowing.cpp',
                        'Common/toolboxes/random_toolbox_tests.cpp',
@@ -55,6 +56,8 @@ if get_option('enable-tests')
         cpp_args: ['-fPIC', default_warning_flags, su2_cpp_args]
     )
     test('Catch2 test driver', test_driver)
+    test('Multigrid configuration limits', python,
+         args: files('Common/check_multigrid_limits.py') + [test_driver.full_path()])
     test('NEMO convergence field migration', python,
          args: files('SU2_CFD/output/check_nemo_convergence_migration.py') + [test_driver.full_path()],
          timeout: 300)