diff --git a/Common/include/option_structure.hpp b/Common/include/option_structure.hpp
index 351e534671ae..9afa25ff6b72 100644
--- a/Common/include/option_structure.hpp
+++ b/Common/include/option_structure.hpp
@@ -2865,6 +2865,8 @@ enum PERIODIC_QUANTITIES {
   PERIODIC_LIM_PRIM_1 ,  /*!< \brief Primitive limiter communication phase 1 of 2 (periodic only). */
   PERIODIC_LIM_PRIM_2 ,  /*!< \brief Primitive limiter communication phase 2 of 2 (periodic only). */
   PERIODIC_IMPLICIT   ,  /*!< \brief Implicit update communication to ensure consistency across periodic boundaries. */
+  PERIODIC_AUXVAR_LS  ,  /*!< \brief Auxiliary scalar gradient communication for weighted least-squares. */
+  PERIODIC_AUXVAR_GG  ,  /*!< \brief Auxiliary scalar gradient communication for Green-Gauss. */
 };
 
 /*!
@@ -3001,12 +3003,12 @@ static const MapType<std::string, ENUM_STREAMWISE_PERIODIC> Streamwise_Periodic_
  * \brief Container to hold Variables for streamwise Periodic flow as they are often used together in places.
  */
 struct StreamwisePeriodicValues {
-  su2double Streamwise_Periodic_PressureDrop;       /*!< \brief Value of prescribed pressure drop [Pa] which results in an artificial body force vector. */
-  su2double Streamwise_Periodic_MassFlow;           /*!< \brief Value of current massflow [kg/s] which results in a delta p and therefore an artificial body force vector. */
-  su2double Streamwise_Periodic_IntegratedHeatFlow; /*!< \brief Value of of the net sum of heatflow [W] into the domain. */
-  su2double Streamwise_Periodic_InletTemperature;   /*!< \brief Area avg static Temp [K] at the periodic inlet. Used for adaptive outlet heatsink. */
-  su2double Streamwise_Periodic_BoundaryArea;       /*!< \brief Global Surface area of the streamwise periodic interface. */
-  su2double Streamwise_Periodic_AvgDensity;         /*!< \brief Area avg density on the periodic interface. */
+  su2double Streamwise_Periodic_PressureDrop = 0.0; /*!< \brief Value of prescribed pressure drop [Pa] which results in an artificial body force vector. */
+  su2double Streamwise_Periodic_MassFlow = 0.0; /*!< \brief Value of current massflow [kg/s] which results in a delta p and therefore an artificial body force vector. */
+  su2double Streamwise_Periodic_IntegratedHeatFlow = 0.0; /*!< \brief Value of of the net sum of heatflow [W] into the domain. */
+  su2double Streamwise_Periodic_InletTemperature = 0.0; /*!< \brief Area avg static Temp [K] at the periodic inlet. Used for adaptive outlet heatsink. */
+  su2double Streamwise_Periodic_BoundaryArea = 0.0; /*!< \brief Global Surface area of the streamwise periodic interface. */
+  su2double Streamwise_Periodic_AvgDensity = 0.0; /*!< \brief Area avg density on the periodic interface. */
 };
 
 /*!
diff --git a/Common/src/CConfig.cpp b/Common/src/CConfig.cpp
index 1ab7ced2d354..e226f76f5720 100644
--- a/Common/src/CConfig.cpp
+++ b/Common/src/CConfig.cpp
@@ -37,6 +37,7 @@
 
 #include "../include/basic_types/ad_structure.hpp"
 #include "../include/toolboxes/printing_toolbox.hpp"
+#include "../include/toolboxes/geometry_toolbox.hpp"
 #include "../include/toolboxes/SwapBytes.hpp"
 
 using namespace PrintingToolbox;
@@ -5558,12 +5559,28 @@ void CConfig::SetPostprocessing(SU2_COMPONENT val_software, unsigned short val_i
       SU2_MPI::Error("Streamwise Periodic Flow + Incompressible Euler: Not tested yet.", CURRENT_FUNCTION);
     if (nMarker_PerBound == 0)
       SU2_MPI::Error("A MARKER_PERIODIC pair has to be set with KIND_STREAMWISE_PERIODIC != NONE.", CURRENT_FUNCTION);
+    /*--- The first pair defines the streamwise direction used by every source. ---*/
+    const auto* translation = GetPeriodic_Translation(0);
+    if (GeometryToolbox::SquaredNorm(3, translation) == 0)
+      SU2_MPI::Error("The first streamwise periodic pair must have a nonzero translation.", CURRENT_FUNCTION);
+    const auto* angles = GetPeriodicRotAngles(Marker_PerBound[0]);
+    su2double rotation[3][3];
+    GeometryToolbox::RotationMatrix(angles[0], angles[1], angles[2], rotation);
+    for (auto i = 0u; i < 3; ++i)
+      for (auto j = 0u; j < 3; ++j)
+        if (fabs(rotation[i][j] - (i == j)) > 1e-12)
+          SU2_MPI::Error("The first streamwise periodic pair must be purely translational.", CURRENT_FUNCTION);
+    if (Energy_Equation && Streamwise_Periodic_Temperature && (nMarker_HeatTransfer || nMarker_CHTInterface))
+      SU2_MPI::Error("STREAMWISE_PERIODIC_TEMPERATURE requires prescribed MARKER_HEATFLUX walls; "
+                     "MARKER_HEATTRANSFER and MARKER_CHT_INTERFACE are unsupported.", CURRENT_FUNCTION);
     if (Energy_Equation && Streamwise_Periodic_Temperature && nMarker_Isothermal != 0)
       SU2_MPI::Error("No MARKER_ISOTHERMAL marker allowed with STREAMWISE_PERIODIC_TEMPERATURE= YES, only MARKER_HEATFLUX & MARKER_SYM.", CURRENT_FUNCTION);
     if (Ref_Inc_NonDim != DIMENSIONAL)
       SU2_MPI::Error("Streamwise Periodicity only works with \"INC_NONDIM= DIMENSIONAL\", the nondimensionalization with source terms doesn;t work in general.", CURRENT_FUNCTION);
     if (Axisymmetric)
       SU2_MPI::Error("Streamwise Periodicity terms does not not have axisymmetric corrections.", CURRENT_FUNCTION);
+    if (nMGLevels != 0)
+      SU2_MPI::Error("Streamwise Periodicity does not support multigrid, use MGLEVEL= 0.", CURRENT_FUNCTION);
     if (!Energy_Equation) Streamwise_Periodic_Temperature = false;
   } else {
     /*--- Safety measure ---*/
diff --git a/SU2_CFD/src/numerics/flow/flow_sources.cpp b/SU2_CFD/src/numerics/flow/flow_sources.cpp
index e3797c7c5a5f..02270d9f4143 100644
--- a/SU2_CFD/src/numerics/flow/flow_sources.cpp
+++ b/SU2_CFD/src/numerics/flow/flow_sources.cpp
@@ -794,8 +794,10 @@ CNumerics::ResidualType<> CSourceIncStreamwise_Periodic::ComputeResidual(const C
     /*--- If a RANS turbulence model is used, an additional source term, based on the eddy viscosity gradient is added. ---*/
     if(turbulent) {
 
-      /*--- Compute a scalar factor ---*/
-      scalar_factor = SPvals.Streamwise_Periodic_IntegratedHeatFlow / (SPvals.Streamwise_Periodic_MassFlow * sqrt(norm2_translation) * Prandtl_Turb);
+      /*--- Compute a scalar factor. The temperature is T = T_periodic + Q / (mdot * cp * |t|^2) * dot_prod(t*x)
+            (Patankar et al., J. Heat Transfer 99(2), 1977), so the turbulent heat flux cp * mu_t / Pr_t * grad(T)
+            has the extra divergence Q / (mdot * |t|^2 * Pr_t) * dot_prod(t*grad(mu_t)). ---*/
+      scalar_factor = SPvals.Streamwise_Periodic_IntegratedHeatFlow / (SPvals.Streamwise_Periodic_MassFlow * norm2_translation * Prandtl_Turb);
 
       /*--- Compute scalar product between periodic translation vector and eddy viscosity gradient. ---*/
       dot_product = GeometryToolbox::DotProduct(nDim, Streamwise_Coord_Vector, AuxVar_Grad_i[0]);
diff --git a/SU2_CFD/src/solvers/CIncEulerSolver.cpp b/SU2_CFD/src/solvers/CIncEulerSolver.cpp
index 512a99a97dc6..a2bbef972cb0 100644
--- a/SU2_CFD/src/solvers/CIncEulerSolver.cpp
+++ b/SU2_CFD/src/solvers/CIncEulerSolver.cpp
@@ -1889,19 +1889,23 @@ void CIncEulerSolver::Source_Residual(CGeometry *geometry, CSolver **solver_cont
       /*--- The Pressure drop is iteratively adapted to result in the prescribed Target-Massflow. ---*/
       /*---------------------------------------------------------------------------------------------*/
 
-      /*--- Compute update to Delta p based on massflow-difference ---*/
-      const su2double Average_Density_Global = SPvals.Streamwise_Periodic_AvgDensity;
-      const su2double Area_Global = SPvals.Streamwise_Periodic_BoundaryArea;
-      const su2double TargetMassFlow = config->GetStreamwise_Periodic_TargetMassFlow() / (config->GetDensity_Ref() * config->GetVelocity_Ref());
-      const su2double MassFlow_Global = SPvals.Streamwise_Periodic_MassFlow;
-      const su2double ddP = 0.5 / ( Average_Density_Global * pow(Area_Global, 2)) * (pow(TargetMassFlow, 2) - pow(MassFlow_Global, 2));
-
-      /*--- Store updated pressure difference ---*/
-      const su2double damping_factor = config->GetInc_Outlet_Damping();
-      SPvalsUpdated = SPvals;
-      SPvalsUpdated.Streamwise_Periodic_PressureDrop += damping_factor*ddP;
-      if (!config->GetDiscrete_Adjoint())
-        SPvals = SPvalsUpdated;
+      BEGIN_SU2_OMP_SAFE_GLOBAL_ACCESS
+      {
+        /*--- Compute update to Delta p based on massflow-difference ---*/
+        const su2double Average_Density_Global = SPvals.Streamwise_Periodic_AvgDensity;
+        const su2double Area_Global = SPvals.Streamwise_Periodic_BoundaryArea;
+        const su2double TargetMassFlow = config->GetStreamwise_Periodic_TargetMassFlow() / (config->GetDensity_Ref() * config->GetVelocity_Ref());
+        const su2double MassFlow_Global = SPvals.Streamwise_Periodic_MassFlow;
+        const su2double ddP = 0.5 / ( Average_Density_Global * pow(Area_Global, 2)) * (pow(TargetMassFlow, 2) - pow(MassFlow_Global, 2));
+
+        /*--- Store updated pressure difference ---*/
+        const su2double damping_factor = config->GetInc_Outlet_Damping();
+        SPvalsUpdated = SPvals;
+        SPvalsUpdated.Streamwise_Periodic_PressureDrop += damping_factor*ddP;
+        if (!config->GetDiscrete_Adjoint())
+          SPvals = SPvalsUpdated;
+      }
+      END_SU2_OMP_SAFE_GLOBAL_ACCESS
 
       /*--- Set delta_p, m_dot, inlet_T, integrated_heat ---*/
       numerics->SetStreamwisePeriodicValues(SPvalsUpdated);
diff --git a/SU2_CFD/src/solvers/CIncNSSolver.cpp b/SU2_CFD/src/solvers/CIncNSSolver.cpp
index ec1b78935e1f..665707907293 100644
--- a/SU2_CFD/src/solvers/CIncNSSolver.cpp
+++ b/SU2_CFD/src/solvers/CIncNSSolver.cpp
@@ -209,7 +209,8 @@ void CIncNSSolver::GetStreamwise_Periodic_Properties(const CGeometry *geometry,
 
         /*--- Identify the boundary by string name and retrive heatflux from config ---*/
         const auto Marker_StringTag = config->GetMarker_All_TagBound(iMarker);
-        const su2double Wall_HeatFlux = config->GetWall_HeatFlux(Marker_StringTag);
+        su2double Wall_HeatFlux = config->GetWall_HeatFlux(Marker_StringTag);
+        if (config->GetIntegrated_HeatFlux()) Wall_HeatFlux /= geometry->GetSurfaceArea(config, iMarker);
 
         for (auto iVertex = 0ul; iVertex < geometry->nVertex[iMarker]; iVertex++) {
 
@@ -254,8 +255,8 @@ void CIncNSSolver::Compute_Streamwise_Periodic_Recovered_Values(CConfig *config,
 
     /*--- First, compute helping terms based on relative distance (0,l) between periodic markers ---*/
     su2double dot_product = 0.0;
-    for (unsigned short iDim = 0; iDim < nDim; iDim++)
-      dot_product += fabs( (geometry->nodes->GetCoord(iPoint,iDim) - ReferenceNode[iDim]) * config->GetPeriodic_Translation(0)[iDim]);
+    for (auto iDim = 0u; iDim < nDim; iDim++)
+      dot_product += (geometry->nodes->GetCoord(iPoint,iDim) - ReferenceNode[iDim]) * config->GetPeriodic_Translation(0)[iDim];
 
     /*--- Second, substract/add correction from reduced pressure/temperature to get recoverd pressure/temperature ---*/
     const su2double Pressure_Recovered = nodes->GetPressure(iPoint) - SPvals.Streamwise_Periodic_PressureDrop /
@@ -537,7 +538,8 @@ void CIncNSSolver::BC_Wall_Generic(const CGeometry *geometry, const CConfig *con
         /*--- Dot product ---*/
         const su2double dot_product = GeometryToolbox::DotProduct(nDim, config->GetPeriodic_Translation(0), Normal);
 
-        LinSysRes(iPoint, nDim+1) += scalar_factor*dot_product;
+        /*--- Normal points into the fluid: q_periodic = q + k grad(T_recovered).n. ---*/
+        LinSysRes(iPoint, nDim+1) -= scalar_factor*dot_product;
       } // if streamwise_periodic
       break;
 
diff --git a/SU2_CFD/src/solvers/CSolver.cpp b/SU2_CFD/src/solvers/CSolver.cpp
index 56cee4b6ad52..c2f6fb7ddd0a 100644
--- a/SU2_CFD/src/solvers/CSolver.cpp
+++ b/SU2_CFD/src/solvers/CSolver.cpp
@@ -229,6 +229,12 @@ void CSolver::GetPeriodicCommCountAndType(const CConfig* config,
       COUNT_PER_POINT  = 2;
       MPI_TYPE         = COMM_TYPE::DOUBLE;
       break;
+    case PERIODIC_AUXVAR_GG:
+      COUNT_PER_POINT  = base_nodes->GetnAuxVar()*nDim;
+      MPI_TYPE         = COMM_TYPE::DOUBLE;
+      ICOUNT           = base_nodes->GetnAuxVar();
+      JCOUNT           = nDim;
+      break;
     case PERIODIC_SOL_GG:
     case PERIODIC_SOL_GG_R:
       COUNT_PER_POINT  = nVar*nDim;
@@ -243,6 +249,12 @@ void CSolver::GetPeriodicCommCountAndType(const CConfig* config,
       ICOUNT           = nPrimVarGrad;
       JCOUNT           = nDim;
       break;
+    case PERIODIC_AUXVAR_LS:
+      COUNT_PER_POINT = static_cast<unsigned long>(nDim)*nDim + base_nodes->GetnAuxVar()*nDim;
+      MPI_TYPE = COMM_TYPE::DOUBLE;
+      ICOUNT = base_nodes->GetnAuxVar();
+      JCOUNT = nDim;
+      break;
     case PERIODIC_SOL_LS:
     case PERIODIC_SOL_ULS:
     case PERIODIC_SOL_LS_R:
@@ -296,6 +308,9 @@ namespace PeriodicCommHelpers {
       case PERIODIC_PRIM_ULS:
         return nodes->GetGradient_Primitive();
         break;
+      case PERIODIC_AUXVAR_LS:
+      case PERIODIC_AUXVAR_GG:
+        return nodes->GetAuxVarGradient();
       case PERIODIC_SOL_GG:
       case PERIODIC_SOL_LS:
       case PERIODIC_SOL_ULS:
@@ -309,6 +324,9 @@ namespace PeriodicCommHelpers {
 
   const su2activematrix& selectField(CVariable* nodes, unsigned short commType) {
     switch(commType) {
+      case PERIODIC_AUXVAR_LS:
+      case PERIODIC_AUXVAR_GG:
+        return nodes->GetAuxVar();
       case PERIODIC_PRIM_GG:
       case PERIODIC_PRIM_LS:
       case PERIODIC_PRIM_ULS:
@@ -356,6 +374,8 @@ void CSolver::InitiatePeriodicComms(CGeometry *geometry,
 
   bool boundary_i, boundary_j;
   bool weighted = true;
+  const bool auxiliary = commType == PERIODIC_AUXVAR_LS;
+  const auto nRotVars = max<unsigned long>(max(nVar, nPrimVar), base_nodes->GetnAuxVar());
 
   unsigned short iVar, jVar, iDim;
   unsigned short nNeighbor       = 0;
@@ -372,8 +392,8 @@ void CSolver::InitiatePeriodicComms(CGeometry *geometry,
   auto *Und_Lapl  = new su2double[nVar];
   auto *Sol_Min   = new su2double[std::max(nVar, nPrimVarGrad)];
   auto *Sol_Max   = new su2double[std::max(nVar, nPrimVarGrad)];
-  auto *rotPrim_i = new su2double[std::max(nVar, nPrimVar)];
-  auto *rotPrim_j = new su2double[std::max(nVar, nPrimVar)];
+  auto *rotPrim_i = new su2double[nRotVars];
+  auto *rotPrim_j = new su2double[nRotVars];
 
   su2double Sensor_i = 0.0, Sensor_j = 0.0, Pressure_i, Pressure_j;
   const su2double *Coord_i, *Coord_j;
@@ -725,6 +745,7 @@ void CSolver::InitiatePeriodicComms(CGeometry *geometry,
           case PERIODIC_SOL_GG_R:
           case PERIODIC_PRIM_GG:
           case PERIODIC_PRIM_GG_R:
+          case PERIODIC_AUXVAR_GG:
 
             /*--- Access and rotate the partial G-G gradient. These will be
              summed on both sides of the periodic faces before dividing
@@ -742,9 +763,9 @@ void CSolver::InitiatePeriodicComms(CGeometry *geometry,
               Rotate(zeros, jacBlock[iVar], rotBlock[iVar]);
             }
 
-            /*--- Rotate the vector components of the solution. ---*/
+            /*--- Auxiliary fields are scalars: rotate their spatial gradients only. ---*/
 
-            if (rotate_periodic) {
+            if (rotate_periodic && commType != PERIODIC_AUXVAR_GG) {
               for (iDim = 0; iDim < nDim; iDim++) {
                 su2double d_diDim[3] = {0.0};
                 for (iVar = 1; iVar < 1+nDim; ++iVar) {
@@ -768,6 +789,7 @@ void CSolver::InitiatePeriodicComms(CGeometry *geometry,
 
             break;
 
+          case PERIODIC_AUXVAR_LS:
           case PERIODIC_SOL_LS: case PERIODIC_SOL_ULS:
           case PERIODIC_SOL_LS_R: case PERIODIC_SOL_ULS_R:
           case PERIODIC_PRIM_LS: case PERIODIC_PRIM_ULS:
@@ -810,7 +832,7 @@ void CSolver::InitiatePeriodicComms(CGeometry *geometry,
             for (iVar = 0; iVar < ICOUNT; iVar++)
               rotPrim_i[iVar] = field(iPoint, iVar);
 
-            if (rotate_periodic) {
+            if (rotate_periodic && !auxiliary) {
               Rotate(zeros, &field(iPoint,1), &rotPrim_i[1]);
             }
 
@@ -845,7 +867,7 @@ void CSolver::InitiatePeriodicComms(CGeometry *geometry,
                 for (iVar = 0; iVar < ICOUNT; iVar++)
                   rotPrim_j[iVar] = field(jPoint,iVar);
 
-                if (rotate_periodic) {
+                if (rotate_periodic && !auxiliary) {
                   Rotate(zeros, &field(jPoint,1), &rotPrim_j[1]);
                 }
 
@@ -1233,6 +1255,7 @@ void CSolver::CompletePeriodicComms(CGeometry *geometry,
             case PERIODIC_SOL_GG_R:
             case PERIODIC_PRIM_GG:
             case PERIODIC_PRIM_GG_R:
+            case PERIODIC_AUXVAR_GG:
 
               /*--- For G-G, we accumulate partial gradients then compute
                the final value using the entire volume of the periodic cell. ---*/
@@ -1243,6 +1266,7 @@ void CSolver::CompletePeriodicComms(CGeometry *geometry,
 
               break;
 
+            case PERIODIC_AUXVAR_LS:
             case PERIODIC_SOL_LS: case PERIODIC_SOL_ULS:
             case PERIODIC_SOL_LS_R: case PERIODIC_SOL_ULS_R:
             case PERIODIC_PRIM_LS: case PERIODIC_PRIM_ULS:
@@ -2204,19 +2228,19 @@ void CSolver::SetAuxVar_Gradient_GG(CGeometry *geometry, const CConfig *config)
   const auto& solution = base_nodes->GetAuxVar();
   auto& gradient = base_nodes->GetAuxVarGradient();
 
-  computeGradientsGreenGauss(this, MPI_QUANTITIES::AUXVAR_GRADIENT, PERIODIC_NONE, *geometry,
+  computeGradientsGreenGauss(this, MPI_QUANTITIES::AUXVAR_GRADIENT, PERIODIC_AUXVAR_GG, *geometry,
                              *config, solution, 0, base_nodes->GetnAuxVar(), -1, gradient);
 }
 
 void CSolver::SetAuxVar_Gradient_LS(CGeometry *geometry, const CConfig *config) {
   SU2_ZONE_SCOPED
 
-  bool weighted = true;
+  const bool weighted = true;
   const auto& solution = base_nodes->GetAuxVar();
   auto& gradient = base_nodes->GetAuxVarGradient();
   auto& rmatrix  = base_nodes->GetRmatrix();
 
-  computeGradientsLeastSquares(this, MPI_QUANTITIES::AUXVAR_GRADIENT, PERIODIC_NONE, *geometry, *config,
+  computeGradientsLeastSquares(this, MPI_QUANTITIES::AUXVAR_GRADIENT, PERIODIC_AUXVAR_LS, *geometry, *config,
                                weighted, solution, 0, base_nodes->GetnAuxVar(), -1, gradient, rmatrix);
 }
 
diff --git a/SU2_CFD/src/variables/CSpeciesFlameletVariable.cpp b/SU2_CFD/src/variables/CSpeciesFlameletVariable.cpp
index a981655b39af..53dedec835cb 100644
--- a/SU2_CFD/src/variables/CSpeciesFlameletVariable.cpp
+++ b/SU2_CFD/src/variables/CSpeciesFlameletVariable.cpp
@@ -56,6 +56,7 @@ CSpeciesFlameletVariable::CSpeciesFlameletVariable(const su2double* species_inf,
   source_cons_jac.resize(nPoint, flamelet_config_options.n_user_scalars) = su2double(0.0);
 
   if (flamelet_config_options.preferential_diffusion) {
+    nAuxVar = FLAMELET_PREF_DIFF_SCALARS::N_BETA_TERMS;
     AuxVar.resize(nPoint, FLAMELET_PREF_DIFF_SCALARS::N_BETA_TERMS) = su2double(0.0);
     Grad_AuxVar.resize(nPoint, FLAMELET_PREF_DIFF_SCALARS::N_BETA_TERMS, nDim, 0.0);
   }
diff --git a/TestCases/hybrid_regression.py b/TestCases/hybrid_regression.py
index f63748ba9283..c7cbd6a97530 100644
--- a/TestCases/hybrid_regression.py
+++ b/TestCases/hybrid_regression.py
@@ -153,6 +153,14 @@ def main():
     poiseuille_profile.test_vals_aarch64 = [-12.004276, -7.636719, -0.000000, 2.089953]
     test_list.append(poiseuille_profile)
 
+    # Streamwise periodic pressure-drop update with two OpenMP threads.
+    sp_massflow_box = TestCase('sp_massflow_box')
+    sp_massflow_box.cfg_dir = "incomp_navierstokes/streamwise_periodic"
+    sp_massflow_box.cfg_file = "massflow_box.cfg"
+    sp_massflow_box.test_iter = 4
+    sp_massflow_box.test_vals = [-4.914357, 0.248162, -0.092770]
+    test_list.append(sp_massflow_box)
+
     # 2D Rotational Periodic
     periodic2d           = TestCase('periodic2d')
     periodic2d.cfg_dir   = "navierstokes/periodic2D"
diff --git a/TestCases/incomp_navierstokes/streamwise_periodic/massflow_box.cfg b/TestCases/incomp_navierstokes/streamwise_periodic/massflow_box.cfg
new file mode 100644
index 000000000000..8763d36ff6df
--- /dev/null
+++ b/TestCases/incomp_navierstokes/streamwise_periodic/massflow_box.cfg
@@ -0,0 +1,36 @@
+% Small generated periodic channel: pressure-drop update must be independent of OpenMP threads.
+SOLVER= INC_NAVIER_STOKES
+MESH_FORMAT= BOX
+MESH_BOX_SIZE= (9, 9, 9)
+MESH_BOX_LENGTH= (1, 1, 1)
+MESH_BOX_OFFSET= (0, 0, 0)
+INC_DENSITY_MODEL= CONSTANT
+INC_DENSITY_INIT= 1.0
+INC_VELOCITY_INIT= (0.1, 0.0, 0.0)
+INC_TEMPERATURE_INIT= 300.0
+INC_ENERGY_EQUATION= NO
+VISCOSITY_MODEL= CONSTANT_VISCOSITY
+MU_CONSTANT= 0.01
+MARKER_HEATFLUX= (y_minus,0, y_plus,0, z_minus,0, z_plus,0)
+MARKER_PERIODIC= (x_minus,x_plus, 0,0,0, 0,0,0, 1,0,0)
+KIND_STREAMWISE_PERIODIC= MASSFLOW
+STREAMWISE_PERIODIC_MASSFLOW= 1.0
+STREAMWISE_PERIODIC_PRESSURE_DROP= 0.0
+INC_OUTLET_DAMPING= 0.1
+NUM_METHOD_GRAD= GREEN_GAUSS
+CONV_NUM_METHOD_FLOW= FDS
+MUSCL_FLOW= NO
+TIME_DISCRE_FLOW= EULER_IMPLICIT
+CFL_NUMBER= 1.0
+CFL_ADAPT= NO
+LINEAR_SOLVER= FGMRES
+LINEAR_SOLVER_PREC= ILU
+LINEAR_SOLVER_ITER= 10
+ITER= 5
+CONV_RESIDUAL_MINVAL= -30
+OUTPUT_FILES= NONE
+OUTPUT_WRT_FREQ= 5
+HISTORY_OUTPUT= (RMS_RES, STREAMWISE_PERIODIC)
+SCREEN_OUTPUT= (INNER_ITER, RMS_PRESSURE, STREAMWISE_DP, STREAMWISE_MASSFLOW)
+SCREEN_WRT_FREQ_INNER= 1
+INC_NONDIM= DIMENSIONAL
diff --git a/TestCases/incomp_navierstokes/streamwise_periodic/pipeSlice_3d/sp_pipeSlice_3d_dp_ihf_tp.cfg b/TestCases/incomp_navierstokes/streamwise_periodic/pipeSlice_3d/sp_pipeSlice_3d_dp_ihf_tp.cfg
new file mode 100644
index 000000000000..37d1ff9a70a9
--- /dev/null
+++ b/TestCases/incomp_navierstokes/streamwise_periodic/pipeSlice_3d/sp_pipeSlice_3d_dp_ihf_tp.cfg
@@ -0,0 +1,85 @@
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%                                                                              %
+% SU2 configuration file                                                       %
+% Case description: Same as sp_pipeSlice_3d_dp_hf_tp.cfg with the energy       %
+%                   equation, periodic temperature and a wall heat given as    %
+%                   integrated value (INTEGRATED_HEATFLUX)                     %
+% File Version 8.5.0 "Harrier"                                                 %
+%                                                                              %
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+% ------------- DIRECT, ADJOINT, AND LINEARIZED PROBLEM DEFINITION ------------%
+%
+SOLVER= INC_NAVIER_STOKES
+%
+KIND_TURB_MODEL= NONE
+%
+% ---------------- INCOMPRESSIBLE FLOW CONDITION DEFINITION -------------------%
+%
+INC_DENSITY_MODEL= CONSTANT
+INC_ENERGY_EQUATION = YES
+INC_DENSITY_INIT= 1.0
+INC_VELOCITY_INIT= ( 0.0, 0.0, 0.3 )
+INC_NONDIM= DIMENSIONAL
+%
+% --------------------------- VISCOSITY MODEL ---------------------------------%
+%
+VISCOSITY_MODEL= CONSTANT_VISCOSITY
+MU_CONSTANT= 1.8e-5
+%
+% --------------------- STREAMWISE PERIODICITY DEFINITION ---------------------%
+%
+KIND_STREAMWISE_PERIODIC= PRESSURE_DROP
+STREAMWISE_PERIODIC_PRESSURE_DROP= 0.001
+STREAMWISE_PERIODIC_TEMPERATURE= YES
+%
+% -------------------- BOUNDARY CONDITION DEFINITION --------------------------%
+%
+% Heat [W] over the whole wall (area 1.570028349e-5 m2, i.e. 1000 W/m2)
+MARKER_HEATFLUX= (wall, 0.01570028349)
+INTEGRATED_HEATFLUX= YES
+MARKER_PERIODIC= ( inlet, outlet, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0005 )
+%
+MARKER_PLOTTING= ( inlet )
+MARKER_MONITORING= ( wall )
+MARKER_ANALYZE = ( oulet )
+MARKER_ANALYZE_AVERAGE = AREA
+%
+% ------------- COMMON PARAMETERS DEFINING THE NUMERICAL METHOD ---------------%
+%
+NUM_METHOD_GRAD= WEIGHTED_LEAST_SQUARES
+%CFL_NUMBER= 1e10
+CFL_NUMBER= 50000
+CFL_ADAPT= NO
+CFL_ADAPT_PARAM= ( 0.5, 10, 15.0, 1e30 )
+ITER= 15000
+%
+% ------------------------ LINEAR SOLVER DEFINITION ---------------------------%
+%
+LINEAR_SOLVER= FGMRES
+LINEAR_SOLVER_PREC= ILU
+LINEAR_SOLVER_ERROR= 1E-15
+LINEAR_SOLVER_ITER= 10
+%
+% -------------------- FLOW NUMERICAL METHOD DEFINITION -----------------------%
+%
+CONV_NUM_METHOD_FLOW= FDS
+MUSCL_FLOW= YES
+SLOPE_LIMITER_FLOW= NONE
+TIME_DISCRE_FLOW= EULER_IMPLICIT
+%
+% --------------------------- CONVERGENCE PARAMETERS --------------------------%
+%
+CONV_RESIDUAL_MINVAL= -24
+CONV_STARTITER= 10
+%
+% ------------------------- INPUT/OUTPUT INFORMATION --------------------------%
+%
+MESH_FILENAME= pipe1cell3D.su2
+%
+OUTPUT_FILES= ( RESTART, PARAVIEW_MULTIBLOCK, SURFACE_TECPLOT_ASCII )
+OUTPUT_WRT_FREQ= 1000
+%
+HISTORY_OUTPUT= ( RMS_RES, FLOW_COEFF, STREAMWISE_PERIODIC, LINSOL )
+SCREEN_OUTPUT= ( INNER_ITER, RMS_PRESSURE, RMS_ENTHALPY, STREAMWISE_MASSFLOW, STREAMWISE_HEAT )
+SCREEN_WRT_FREQ_INNER= 100
diff --git a/TestCases/parallel_regression.py b/TestCases/parallel_regression.py
index e79219efbeba..a512ede85c07 100755
--- a/TestCases/parallel_regression.py
+++ b/TestCases/parallel_regression.py
@@ -780,6 +780,14 @@ def main():
     sp_pipeSlice_3d_dp_hf_tp.test_vals = [-11.119796, -11.234737, -9.429479, -0.000023]
     test_list.append(sp_pipeSlice_3d_dp_hf_tp)
 
+    # Same with the energy equation, periodic temperature and an integrated wall heat
+    sp_pipeSlice_3d_dp_ihf_tp = TestCase('sp_pipeSlice_3d_dp_ihf_tp')
+    sp_pipeSlice_3d_dp_ihf_tp.cfg_dir = "incomp_navierstokes/streamwise_periodic/pipeSlice_3d"
+    sp_pipeSlice_3d_dp_ihf_tp.cfg_file = "sp_pipeSlice_3d_dp_ihf_tp.cfg"
+    sp_pipeSlice_3d_dp_ihf_tp.test_iter = 10
+    sp_pipeSlice_3d_dp_ihf_tp.test_vals = [-11.120123, -5.196299, -0.000023, -0.015700]
+    test_list.append(sp_pipeSlice_3d_dp_ihf_tp)
+
     # 2D pin array with heat transfer BC on pin surfaces
     inc_heatTransfer_BC           = TestCase('inc_heatTransfer_BC')
     inc_heatTransfer_BC.cfg_dir   = "incomp_navierstokes/streamwise_periodic/chtPinArray_2d"
diff --git a/UnitTests/Common/CConfig_tests.cpp b/UnitTests/Common/CConfig_tests.cpp
index 111918785545..13e45452679c 100644
--- a/UnitTests/Common/CConfig_tests.cpp
+++ b/UnitTests/Common/CConfig_tests.cpp
@@ -87,3 +87,48 @@ TEST_CASE("INIT_OPTION_INC defaults", "[Config]") {
         INIT_OPTION_INC::OPERATING_PRESSURE);
   CHECK(GetInitOptionInc(ideal_gas_options + "INIT_OPTION_INC= DENSITY_INIT\n") == INIT_OPTION_INC::DENSITY_INIT);
 }
+
+TEST_CASE("Streamwise periodic Rotation support", "[.StreamwisePeriodicSupportRotation]") {
+  auto base = base_options;
+  const std::string custom = "MARKER_CUSTOM= (x_minus, x_plus, z_plus, z_minus)";
+  base.replace(base.find(custom), custom.size(), "MARKER_CUSTOM= (z_minus)\nMARKER_SYM= (z_plus)");
+  std::stringstream options(base +
+                            "INC_NONDIM= DIMENSIONAL\nKIND_STREAMWISE_PERIODIC= PRESSURE_DROP\n"
+                            "INC_ENERGY_EQUATION= YES\nSTREAMWISE_PERIODIC_TEMPERATURE= YES\n"
+                            "MARKER_PERIODIC= (x_minus,x_plus, 0,0,0, 90,0,0, 1,0,0)\n");
+  CConfig config(options, SU2_COMPONENT::SU2_CFD, false);
+}
+
+TEST_CASE("Streamwise periodic Zero support", "[.StreamwisePeriodicSupportZero]") {
+  auto base = base_options;
+  const std::string custom = "MARKER_CUSTOM= (x_minus, x_plus, z_plus, z_minus)";
+  base.replace(base.find(custom), custom.size(), "MARKER_CUSTOM= (z_minus)\nMARKER_SYM= (z_plus)");
+  std::stringstream options(base +
+                            "INC_NONDIM= DIMENSIONAL\nKIND_STREAMWISE_PERIODIC= PRESSURE_DROP\n"
+                            "INC_ENERGY_EQUATION= YES\nSTREAMWISE_PERIODIC_TEMPERATURE= YES\n"
+                            "MARKER_PERIODIC= (x_minus,x_plus, 0,0,0, 0,0,0, 0,0,0)\n");
+  CConfig config(options, SU2_COMPONENT::SU2_CFD, false);
+}
+
+TEST_CASE("Streamwise periodic Convection support", "[.StreamwisePeriodicSupportConvection]") {
+  auto base = base_options;
+  const std::string custom = "MARKER_CUSTOM= (x_minus, x_plus, z_plus, z_minus)";
+  base.replace(base.find(custom), custom.size(), "MARKER_CUSTOM= (z_minus)");
+  std::stringstream options(
+      base +
+      "INC_NONDIM= DIMENSIONAL\nKIND_STREAMWISE_PERIODIC= PRESSURE_DROP\n"
+      "INC_ENERGY_EQUATION= YES\nSTREAMWISE_PERIODIC_TEMPERATURE= YES\n"
+      "MARKER_PERIODIC= (x_minus,x_plus, 0,0,0, 0,0,0, 1,0,0)\nMARKER_HEATTRANSFER= (z_plus,1,300)\n");
+  CConfig config(options, SU2_COMPONENT::SU2_CFD, false);
+}
+
+TEST_CASE("Streamwise periodic CHT support", "[.StreamwisePeriodicSupportCHT]") {
+  auto base = base_options;
+  const std::string custom = "MARKER_CUSTOM= (x_minus, x_plus, z_plus, z_minus)";
+  base.replace(base.find(custom), custom.size(), "MARKER_CUSTOM= (z_minus)");
+  std::stringstream options(base +
+                            "INC_NONDIM= DIMENSIONAL\nKIND_STREAMWISE_PERIODIC= PRESSURE_DROP\n"
+                            "INC_ENERGY_EQUATION= YES\nSTREAMWISE_PERIODIC_TEMPERATURE= YES\n"
+                            "MARKER_PERIODIC= (x_minus,x_plus, 0,0,0, 0,0,0, 1,0,0)\nMARKER_CHT_INTERFACE= (z_plus)\n");
+  CConfig config(options, SU2_COMPONENT::SU2_CFD, false);
+}
diff --git a/UnitTests/Common/check_periodic_support.py b/UnitTests/Common/check_periodic_support.py
new file mode 100644
index 000000000000..5aa6cdf02540
--- /dev/null
+++ b/UnitTests/Common/check_periodic_support.py
@@ -0,0 +1,50 @@
+#!/usr/bin/env python3
+"""Check periodic configuration errors in isolated processes (MPI errors terminate)."""
+import argparse
+import subprocess
+
+
+def main():
+    parser = argparse.ArgumentParser(description=__doc__)
+    parser.add_argument("test_driver")
+    parser.add_argument("--streamwise", action="store_true")
+    parser.add_argument("--launcher", nargs=argparse.REMAINDER, default=[])
+    args = parser.parse_args()
+    cases = (
+        (
+            ("Rotation", "must be purely translational"),
+            ("Zero", "must have a nonzero translation"),
+            (
+                "Convection",
+                "MARKER_HEATTRANSFER and MARKER_CHT_INTERFACE are unsupported",
+            ),
+            ("CHT", "MARKER_HEATTRANSFER and MARKER_CHT_INTERFACE are unsupported"),
+        )
+        if args.streamwise
+        else (
+            ("Continuous", "Continuous adjoints do not implement MARKER_PERIODIC"),
+            ("Radiation", "RADIATION_MODEL does not implement MARKER_PERIODIC"),
+            ("Structure", "SOLVER= ELASTICITY does not implement MARKER_PERIODIC"),
+        )
+    )
+    prefix = "StreamwisePeriodicSupport" if args.streamwise else "PeriodicSupport"
+    for tag, message in cases:
+        result = subprocess.run(
+            args.launcher + [args.test_driver, "[." + prefix + tag + "]"],
+            stdout=subprocess.PIPE,
+            stderr=subprocess.STDOUT,
+            text=True,
+            timeout=30,
+            check=False,
+        )
+        if result.returncode == 0 or message not in result.stdout:
+            raise RuntimeError(
+                tag
+                + " did not report the expected configuration error:\n"
+                + result.stdout
+            )
+        print("PASS: " + tag, flush=True)
+
+
+if __name__ == "__main__":
+    main()
diff --git a/UnitTests/SU2_CFD/gradients.cpp b/UnitTests/SU2_CFD/gradients.cpp
index 72baf525ca68..3a97294446c5 100644
--- a/UnitTests/SU2_CFD/gradients.cpp
+++ b/UnitTests/SU2_CFD/gradients.cpp
@@ -26,11 +26,14 @@
  */
 
 #include "catch.hpp"
+#include "../../SU2_CFD/include/variables/CSpeciesFlameletVariable.hpp"
+#include "../../SU2_CFD/include/solvers/CIncNSSolver.hpp"
 #include "../../Common/include/geometry/CPhysicalGeometry.hpp"
 #include "../../Common/include/containers/container_decorators.hpp"
 #include "../../SU2_CFD/include/solvers/CSolver.hpp"
 #include "../../SU2_CFD/include/gradients/computeGradientsGreenGauss.hpp"
 #include "../../SU2_CFD/include/gradients/computeGradientsLeastSquares.hpp"
+#include "../UnitQuadTestCase.hpp"
 
 /*!
  * \brief Base class for gradient tests using a unit cube geometry.
@@ -154,3 +157,211 @@ TEST_CASE("GG", "[Gradients]") { testGreenGauss<LinearFunction>(); }
 TEST_CASE("LS", "[Gradients]") { testLeastSquares<LinearFunction>(false); }
 
 TEST_CASE("WLS", "[Gradients]") { testLeastSquares<LinearFunction>(true); }
+
+TEST_CASE("Periodic auxiliary gradients", "[Gradients][Periodic]") {
+  const bool rotation = GENERATE(false, true);
+  const bool stretched = GENERATE(false, true);
+  const auto method = GENERATE(std::string("GREEN_GAUSS"), std::string("WEIGHTED_LEAST_SQUARES"));
+  UnitQuadTestCase field;
+  const string custom = "MARKER_CUSTOM= ( x_minus, x_plus, z_plus, z_minus)";
+  field.config_options.replace(field.config_options.find(custom), custom.size(), "MARKER_CUSTOM= (z_plus, z_minus)");
+  field.AddOption(rotation ? "MARKER_PERIODIC= (x_minus,x_plus, 0,0.5,0.5, 90,0,0, 1,0,0)"
+                           : "MARKER_PERIODIC= (x_minus,x_plus, 0,0,0, 0,0,0, 1,0,0)");
+  field.AddOption("VORTICITY_CONFINEMENT= YES");
+  field.AddOption("NUM_METHOD_GRAD= " + method);
+  field.InitConfig();
+  field.InitGeometry(true);
+  if (stretched) {
+    for (auto i = 0ul; i < field.geometry->GetnPoint(); ++i) {
+      const auto x = field.geometry->nodes->GetCoord(i, 0);
+      field.geometry->nodes->SetCoord(i, 0, x * x);
+    }
+    field.geometry->SetControlVolume(field.config.get(), UPDATE);
+    field.geometry->SetBoundControlVolume(field.config.get(), UPDATE);
+  }
+  field.geometry->MatchPeriodic(field.config.get(), 1);
+  field.geometry->PreprocessPeriodicComms(field.geometry.get(), field.config.get());
+  field.InitSolver();
+
+  auto* solver = field.solver[FLOW_SOL];
+  auto* nodes = solver->GetNodes();
+  REQUIRE(nodes->GetnAuxVar() == 1);
+  for (auto iPoint = 0ul; iPoint < field.geometry->GetnPoint(); ++iPoint) {
+    const auto* coord = field.geometry->nodes->GetCoord(iPoint);
+    const auto y = coord[1] - 0.5;
+    const auto z = coord[2] - 0.5;
+    nodes->SetAuxVar(iPoint, 0,
+                     (rotation ? 2.0 + y * y + z * z : 2.0 + 3.0 * coord[1] - coord[2]) + coord[0] * (1 - coord[0]));
+  }
+  if (method == "GREEN_GAUSS") {
+    SU2_OMP_PARALLEL { solver->SetAuxVar_Gradient_GG(field.geometry.get(), field.config.get()); }
+    END_SU2_OMP_PARALLEL
+  } else {
+    SU2_OMP_PARALLEL { solver->SetAuxVar_Gradient_LS(field.geometry.get(), field.config.get()); }
+    END_SU2_OMP_PARALLEL
+  }
+
+  su2double error = 0.0;
+  for (auto iPoint = 0ul; iPoint < field.geometry->GetnPointDomain(); ++iPoint) {
+    const auto* coord = field.geometry->nodes->GetCoord(iPoint);
+    if (rotation && (coord[1] == 0.0 || coord[1] == 1.0 || coord[2] == 0.0 || coord[2] == 1.0)) continue;
+    const auto index = static_cast<unsigned short>(
+        std::round((stretched ? sqrt(SU2_TYPE::GetValue(coord[0])) : SU2_TYPE::GetValue(coord[0])) * 4));
+    auto xcoord = [&](unsigned short k) {
+      const auto x = 0.25 * k;
+      return stretched ? x * x : x;
+    };
+    const auto left = xcoord(index == 0 ? 3 : index - 1), right = xcoord(index == 4 ? 1 : index + 1);
+    const auto hm = coord[0] - (index == 0 ? left - 1 : left);
+    const auto hp = (index == 4 ? right + 1 : right) - coord[0];
+    const auto fm = left * (1 - left), fp = right * (1 - right), fi = coord[0] * (1 - coord[0]);
+    const auto gx = method == "GREEN_GAUSS" ? (fp - fm) / (hp + hm) : ((fp - fi) / hp + (fi - fm) / hm) / 2;
+    const su2double expected[3] = {gx, rotation ? 2.0 * (coord[1] - 0.5) : 3.0,
+                                   rotation ? 2.0 * (coord[2] - 0.5) : -1.0};
+    for (auto iDim = 0u; iDim < 3; ++iDim)
+      error = max(error, abs(nodes->GetAuxVarGradient(iPoint, 0, iDim) - expected[iDim]));
+  }
+  CHECK(error < 1e-9);
+}
+
+TEST_CASE("Flamelet preferential diffusion auxiliary gradients", "[Gradients][Flamelet][Periodic]") {
+  const bool periodic = GENERATE(false, true);
+  const auto method = GENERATE(std::string("GREEN_GAUSS"), std::string("WEIGHTED_LEAST_SQUARES"));
+  UnitQuadTestCase field;
+  field.SetOption("SOLVER= INC_NAVIER_STOKES");
+  field.SetOption("KIND_VERIFICATION_SOLUTION= NO_VERIFICATION_SOLUTION");
+  field.AddOption("KIND_SCALAR_MODEL= FLAMELET");
+  field.AddOption("FLUID_MODEL= FLUID_FLAMELET");
+  field.AddOption("INC_DENSITY_MODEL= VARIABLE");
+  field.SetOption("VISCOSITY_MODEL= FLAMELET");
+  field.AddOption("CONDUCTIVITY_MODEL= FLAMELET");
+  field.AddOption("DIFFUSIVITY_MODEL= FLAMELET");
+  field.AddOption(
+      "SPECIES_INIT= (0,300000,0.5)\nCONTROLLING_VARIABLE_NAMES= (ProgressVariable,EnthalpyTot,MixtureFraction)");
+  field.AddOption("CONTROLLING_VARIABLE_SOURCE_NAMES= (ProdRateTot_PV,NULL,NULL)");
+  field.AddOption("PREFERENTIAL_DIFFUSION= YES\nNUM_METHOD_GRAD= " + method);
+  if (periodic) {
+    const string custom = "MARKER_CUSTOM= ( x_minus, x_plus, z_plus, z_minus)";
+    field.config_options.replace(field.config_options.find(custom), custom.size(), "MARKER_CUSTOM= (z_plus,z_minus)");
+    field.AddOption("MARKER_PERIODIC= (x_minus,x_plus, 0,0,0, 0,0,0, 1,0,0)");
+  }
+  field.InitConfig();
+  field.InitGeometry(true);
+  if (periodic) {
+    field.geometry->MatchPeriodic(field.config.get(), 1);
+    field.geometry->PreprocessPeriodicComms(field.geometry.get(), field.config.get());
+  }
+  const su2double initial[] = {0, 300000, 0.5};
+  CSpeciesFlameletVariable nodes(initial, field.geometry->GetnPoint(), 3, 3, field.config.get());
+  CHECK(nodes.GetnAuxVar() == FLAMELET_PREF_DIFF_SCALARS::N_BETA_TERMS);
+  /*--- Exercise the shared auxiliary communication with four beta fields and
+   * three transported scalars, without requiring an unrelated flamelet table. ---*/
+  struct AuxiliarySolver : CSolver {
+    CVariable* nodes;
+    CVariable* GetBaseClassPointerToNodes() override { return nodes; }
+    AuxiliarySolver(CVariable& values, CGeometry& geometry) : nodes(&values) {
+      nDim = 3;
+      nVar = 3;
+      nPoint = geometry.GetnPoint();
+      nPointDomain = geometry.GetnPointDomain();
+      SetBaseClassPointerToNodes();
+    }
+  } solver(nodes, *field.geometry);
+  if (periodic) {
+    /*--- The flow solver normally supplies the complete periodic volume. ---*/
+    SU2_OMP_PARALLEL {
+      solver.InitiatePeriodicComms(field.geometry.get(), field.config.get(), 1, PERIODIC_VOLUME);
+      solver.CompletePeriodicComms(field.geometry.get(), field.config.get(), 1, PERIODIC_VOLUME);
+    }
+    END_SU2_OMP_PARALLEL
+  }
+  for (auto i = 0ul; i < field.geometry->GetnPoint(); ++i) {
+    const auto* x = field.geometry->nodes->GetCoord(i);
+    for (auto v = 0u; v < FLAMELET_PREF_DIFF_SCALARS::N_BETA_TERMS; ++v)
+      nodes.SetAuxVar(i, v, (v + 1) * (2 + (periodic ? x[0] * (1 - x[0]) : x[0]) + 3 * x[1] - x[2]));
+  }
+  if (method == "GREEN_GAUSS") {
+    SU2_OMP_PARALLEL { solver.SetAuxVar_Gradient_GG(field.geometry.get(), field.config.get()); }
+    END_SU2_OMP_PARALLEL
+  } else {
+    SU2_OMP_PARALLEL { solver.SetAuxVar_Gradient_LS(field.geometry.get(), field.config.get()); }
+    END_SU2_OMP_PARALLEL
+  }
+  su2double error = 0;
+  for (auto i = 0ul; i < field.geometry->GetnPointDomain(); ++i) {
+    const auto* x = field.geometry->nodes->GetCoord(i);
+    const su2double expected[] = {periodic ? (x[0] == 0 || x[0] == 1 ? 0.0 : 1 - 2 * x[0]) : 1.0, 3, -1};
+    for (auto v = 0u; v < FLAMELET_PREF_DIFF_SCALARS::N_BETA_TERMS; ++v)
+      for (auto d = 0u; d < 3; ++d)
+        error = std::max(error, fabs(nodes.GetAuxVarGradient(i, v, d) - (v + 1) * expected[d]));
+  }
+  CHECK(error < 1e-12);
+}
+
+TEST_CASE("Recovered periodic wall heat flux", "[Periodic][Heat]") {
+  UnitQuadTestCase field;
+  field.SetOption("SOLVER= INC_NAVIER_STOKES");
+  field.AddOption("FLUID_MODEL= CONSTANT_DENSITY");
+  field.SetOption("KIND_VERIFICATION_SOLUTION= NO_VERIFICATION_SOLUTION");
+  field.AddOption("INC_NONDIM= DIMENSIONAL");
+  field.AddOption("INC_ENERGY_EQUATION= YES\nKIND_STREAMWISE_PERIODIC= PRESSURE_DROP");
+  field.AddOption("STREAMWISE_PERIODIC_TEMPERATURE= YES\nINC_VELOCITY_INIT= (1,0,0)");
+  const string walls = "MARKER_HEATFLUX= (y_minus, 0.0, y_plus, 0.0)";
+  field.config_options.replace(field.config_options.find(walls), walls.size(), "MARKER_HEATFLUX= (y_minus,1,y_plus,3)");
+  const string custom = "MARKER_CUSTOM= ( x_minus, x_plus, z_plus, z_minus)";
+  field.config_options.replace(field.config_options.find(custom), custom.size(), "MARKER_CUSTOM= (z_plus,z_minus)");
+  field.AddOption("MARKER_PERIODIC= (x_minus,x_plus, 0,0,0, 0,0,0, 1,0,0)");
+  field.InitConfig();
+  field.InitGeometry(true);
+  for (auto i = 0ul; i < field.geometry->GetnPoint(); ++i) {
+    const auto* x = field.geometry->nodes->GetCoord(i);
+    field.geometry->nodes->SetCoord(i, 1, x[1] * (1 + 0.1 * sin(2 * PI_NUMBER * x[0])));
+  }
+  field.geometry->SetControlVolume(field.config.get(), UPDATE);
+  field.geometry->SetBoundControlVolume(field.config.get(), UPDATE);
+  field.geometry->MatchPeriodic(field.config.get(), 1);
+  field.geometry->PreprocessPeriodicComms(field.geometry.get(), field.config.get());
+  field.InitSolver();
+  auto* solver = static_cast<CIncNSSolver*>(field.solver[FLOW_SOL]);
+  field.config->SetGlobalParam(field.config->GetKind_Solver(), RUNTIME_FLOW_SYS);
+  {
+    SU2_OMP_PARALLEL {
+      solver->Preprocessing(field.geometry.get(), field.solver, field.config.get(), MESH_0, 0, RUNTIME_FLOW_SYS, false);
+    }
+    END_SU2_OMP_PARALLEL
+  }
+  const auto values = solver->GetStreamwisePeriodicValues();
+  REQUIRE(fabs(values.Streamwise_Periodic_MassFlow) > 1e-12);
+  REQUIRE(fabs(values.Streamwise_Periodic_IntegratedHeatFlow) > 1e-12);
+  auto* nodes = solver->GetNodes();
+  su2double error = 0;
+  unsigned long checked = 0, global = 0;
+  for (auto m = 0u; m < field.geometry->GetnMarker(); ++m) {
+    if (field.config->GetMarker_All_TagBound(m) != "y_plus") continue;
+    solver->LinSysRes = su2double(0);
+    {
+      SU2_OMP_PARALLEL {
+        solver->BC_HeatFlux_Wall(field.geometry.get(), field.solver, nullptr, nullptr, field.config.get(), m);
+      }
+      END_SU2_OMP_PARALLEL
+    }
+    for (auto v = 0ul; v < field.geometry->GetnVertex(m); ++v) {
+      const auto i = field.geometry->vertex[m][v]->GetNode();
+      if (!field.geometry->nodes->GetDomain(i)) continue;
+      const auto* normal = field.geometry->vertex[m][v]->GetNormal();
+      if (fabs(normal[0]) < 1e-12) continue;
+      const auto area = GeometryToolbox::Norm(3, normal);
+      /*--- T = T_periodic + g.x. With an inward area normal, the periodic
+       * conductive heat flux entering the fluid is q + k g.n. ---*/
+      const auto gx = values.Streamwise_Periodic_IntegratedHeatFlow /
+                      (values.Streamwise_Periodic_MassFlow * nodes->GetSpecificHeatCp(i));
+      const auto expected = -field.config->GetWall_HeatFlux("y_plus") / field.config->GetHeat_Flux_Ref() * area -
+                            nodes->GetThermalConductivity(i) * gx * normal[0];
+      error = std::max(error, fabs(solver->LinSysRes(i, 4) - expected));
+      ++checked;
+    }
+  }
+  SU2_MPI::Allreduce(&checked, &global, 1, MPI_UNSIGNED_LONG, MPI_SUM, SU2_MPI::GetComm());
+  REQUIRE(global > 0);
+  CHECK(error < 1e-12);
+}
diff --git a/UnitTests/UnitQuadTestCase.hpp b/UnitTests/UnitQuadTestCase.hpp
index d8da75dc0960..d943c1696ee5 100644
--- a/UnitTests/UnitQuadTestCase.hpp
+++ b/UnitTests/UnitQuadTestCase.hpp
@@ -61,6 +61,16 @@ struct UnitQuadTestCase {
    */
   void AddOption(const std::string& optionLine) { config_options += optionLine + "\n"; }
 
+  /*! \brief Replace one existing base option without repeating its key. */
+  void SetOption(const std::string& optionLine) {
+    const auto key = optionLine.substr(0, optionLine.find('=') + 1);
+    const auto start = config_options.find(key);
+    if (start == std::string::npos)
+      AddOption(optionLine);
+    else
+      config_options.replace(start, config_options.find('\n', start) - start, optionLine);
+  }
+
   /*!
    * \brief Initialize the config structure
    */
@@ -83,10 +93,13 @@ struct UnitQuadTestCase {
   /*!
    * \brief Initialize the geometry
    */
-  void InitGeometry() {
+  void InitGeometry(bool partition = false) {
     cout.rdbuf(nullptr);
     {
-      auto aux_geometry = std::unique_ptr<CGeometry>(new CPhysicalGeometry(config.get(), 0, 1));
+      const auto rank = partition ? SU2_MPI::GetRank() : 0;
+      const auto size = partition ? SU2_MPI::GetSize() : 1;
+      auto aux_geometry = std::unique_ptr<CGeometry>(new CPhysicalGeometry(config.get(), rank, size));
+      if (partition) aux_geometry->SetColorGrid_Parallel(config.get());
       geometry = std::unique_ptr<CGeometry>(new CPhysicalGeometry(aux_geometry.get(), config.get()));
     }
     geometry->SetSendReceive(config.get());
diff --git a/UnitTests/meson.build b/UnitTests/meson.build
index ba3c63afce92..7fa19eb8b3fb 100644
--- a/UnitTests/meson.build
+++ b/UnitTests/meson.build
@@ -50,6 +50,9 @@ if get_option('enable-tests')
         cpp_args: ['-fPIC', default_warning_flags, su2_cpp_args]
     )
     test('Catch2 test driver', test_driver)
+    test('Streamwise periodic support', python,
+         args: files('Common/check_periodic_support.py') + [test_driver.full_path(), '--streamwise'],
+         timeout: 150)
     test('NEMO convergence field migration', python,
          args: files('SU2_CFD/output/check_nemo_convergence_migration.py') + [test_driver.full_path()],
          timeout: 300)