Fix the incompressible moving-grid Jacobian and several smaller mean-flow bugs

Common/src/CConfig.cpp

@@ -4607,6 +4607,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. ---*/
@@ -4966,6 +4996,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;
 

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). ---*/
 

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();

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,27 +1263,41 @@ 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++)
         mat[iVar * nVar + iVar] = 1;
       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
   }

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);
     }
   }
 

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);
     }
   }
 

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++)