Fix the pressure-based solver with variable density, boundary conditions and SIMPLEC

Common/src/CConfig.cpp

@@ -4249,11 +4249,6 @@ void CConfig::SetPostprocessing(SU2_COMPONENT val_software, unsigned short val_i
       SU2_MPI::Error("KIND_INCOMP_SYSTEM= PRESSURE_BASED does not support MGLEVEL > 0,\n"
                      "       the Poisson solver is single-grid only.", CURRENT_FUNCTION);
     }
-    if (Time_Domain) {
-      SU2_MPI::Error("KIND_INCOMP_SYSTEM= PRESSURE_BASED does not support TIME_DOMAIN= YES,\n"
-                     "       it converges to a physically wrong solution instead of failing.",
-                     CURRENT_FUNCTION);
-    }
     if (DiscreteAdjoint || ContinuousAdjoint) {
       SU2_MPI::Error("KIND_INCOMP_SYSTEM= PRESSURE_BASED has no adjoint formulation.", CURRENT_FUNCTION);
     }
@@ -4264,16 +4259,16 @@ void CConfig::SetPostprocessing(SU2_COMPONENT val_software, unsigned short val_i
       SU2_MPI::Error("KIND_INCOMP_SYSTEM= PRESSURE_BASED does not support streamwise periodicity.",
                      CURRENT_FUNCTION);
     }
+    /*--- The solution of the previous time step is only stored with dual time stepping. ---*/
+    if (TimeMarching == TIME_MARCHING::TIME_STEPPING) {
+      SU2_MPI::Error("KIND_INCOMP_SYSTEM= PRESSURE_BASED does not support TIME_MARCHING= TIME_STEPPING,\n"
+                     "       use DUAL_TIME_STEPPING-1ST_ORDER or DUAL_TIME_STEPPING-2ND_ORDER.", CURRENT_FUNCTION);
+    }
 
-    /*--- A_p already carries Vol/dt when SIMPLEC's A_p-Sum_A_nb correction runs, so at the 0.0
-     * default that correction collapses to roughly Vol/dt and the pressure correction becomes
-     * vanishingly weak at low CFL. ---*/
-    if (Kind_PBIter == PBITER::SIMPLEC && !OptionIsSet("TRANSIENT_TERM_REMOVAL_FACTOR")) {
-      SIMPLE_Options.Transient_Term_Removal_Factor = 1.0;
-      if (rank == MASTER_NODE) {
-        cout << "WARNING: KIND_PB_ITER= SIMPLEC without TRANSIENT_TERM_REMOVAL_FACTOR set - "
-             << "defaulting it to 1.0, its intended companion value for SIMPLEC." << endl;
-      }
+    /*--- With SIMPLEC the momentum coefficient a_P - sum(a_nb) is the pseudo-time term alone, removing
+     * it leaves nothing. ---*/
+    if (Kind_PBIter == PBITER::SIMPLEC && SIMPLE_Options.Transient_Term_Removal_Factor != 0.0) {
+      SU2_MPI::Error("KIND_PB_ITER= SIMPLEC requires TRANSIENT_TERM_REMOVAL_FACTOR= 0.", CURRENT_FUNCTION);
     }
   }
 

SU2_CFD/include/numerics/poisson_edge_flux.hpp

@@ -55,21 +55,16 @@ class CScalarFlux_Poisson final
 
   /*!
    * \brief Momentum coefficient, an i/j average, identical for both edge sides (TSL Jacobian).
-   * \note A point under a strong velocity BC has no momentum equation, and so no momentum
-   *       coefficient of its own; the edge uses that of its other node instead.
+   * \note A point under a strong velocity BC has no momentum equation; its coefficient is the average
+   *       of its neighbours' (see CPoissonSolver::SetMomCoeff).
    */
   template <class VariableType>
   FORCEINLINE CPair<Vector<Double, nVar>> coefficients(const FlowIndices&, Int iPoint,
                                                        const EdgeSide<VariableType>& side_i, Int jPoint,
                                                        const EdgeSide<VariableType>& side_j,
                                                        const CPair<Double>&) const {
-    const bool strong_i = side_i.flowNodes->GetStrongBC(iPoint);
-    const bool strong_j = side_j.flowNodes->GetStrongBC(jPoint);
-
-    const Double coeff_i = strong_i ? gatherVariables(jPoint, side_j.scalarNodes.GetMomCoeff())
-                                    : gatherVariables(iPoint, side_i.scalarNodes.GetMomCoeff());
-    const Double coeff_j = strong_j ? gatherVariables(iPoint, side_i.scalarNodes.GetMomCoeff())
-                                    : gatherVariables(jPoint, side_j.scalarNodes.GetMomCoeff());
+    const Double coeff_i = gatherVariables(iPoint, side_i.scalarNodes.GetMomCoeff());
+    const Double coeff_j = gatherVariables(jPoint, side_j.scalarNodes.GetMomCoeff());
 
     Vector<Double, nVar> D;
     D(0) = 0.5 * (coeff_i + coeff_j);

SU2_CFD/include/solvers/CPoissonSolver.hpp

@@ -131,10 +131,16 @@ class CPoissonSolver final : public CScalarSolver<CPoissonVariable> {
                     unsigned long Iteration) override {}
 
   /*!
-   * \brief No dual time stepping as there is no time stepping at all.
+   * \brief Add the change of density in time to the mass balance (variable density only).
+   * \param[in] geometry - Geometrical definition of the problem.
+   * \param[in] solver_container - Container vector with all the solutions.
+   * \param[in] config - Definition of the particular problem.
+   * \param[in] iRKStep - Current step of the Runge-Kutta iteration.
+   * \param[in] iMesh - Index of the mesh in multigrid computations.
+   * \param[in] RunTime_EqSystem - System of equations which is going to be solved.
    */
   void SetResidual_DualTime(CGeometry* geometry, CSolver** solver_container, CConfig* config, unsigned short iRKStep,
-                            unsigned short iMesh, unsigned short RunTime_EqSystem) override {}
+                            unsigned short iMesh, unsigned short RunTime_EqSystem) override;
 
   /*!
    * \brief The pressure correction is reset to zero every iteration (see Preprocessing), so it

SU2_CFD/src/numerics/flow/convection/pressure_based.cpp

@@ -144,20 +144,22 @@ void CPBConvection_Base::ComputeJacobian(su2double val_density, const su2double
 
   /*--- MassFlux is frozen here: it comes from the previous Rhie-Chow interpolation, not from the
   state being solved for. The only derivative left is d(m_f*v_adv)/du, and val_scale carries the
-  advection scheme's weight since v_adv is what the scheme upwinds or averages. ---*/
+  advection scheme's weight since v_adv is what the scheme upwinds or averages. The Jacobian is per
+  unit velocity (A_p = J/rho in CPoissonSolver::SetMomCoeff), so no density appears here. ---*/
 
   for (jDim = 0; jDim < nDim; jDim++) {
     for (iDim = 0; iDim < nDim; iDim++) {
-      val_Proj_Jac_Tensor[iDim+1][jDim+1] = val_density * val_scale * proj_vel * delta[iDim][jDim];
+      val_Proj_Jac_Tensor[iDim+1][jDim+1] = val_scale * MassFlux * delta[iDim][jDim];
     }
   }
 
   /*--- With m_f frozen, Flux[nDim+1] = m_f*h_adv depends only on enthalpy and Flux[1+iDim] =
-  m_f*v_adv only on velocity, so the velocity<->enthalpy cross terms are exactly zero. ---*/
+  m_f*v_adv only on velocity, so the velocity<->enthalpy cross terms are exactly zero and the
+  enthalpy diagonal is the same as the velocity one. ---*/
 
   val_Proj_Jac_Tensor[nDim+1][0] = 0.0;
   val_Proj_Jac_Tensor[0][nDim+1] = val_scale * ((val_dRhodh) * proj_vel);
-  val_Proj_Jac_Tensor[nDim+1][nDim+1] = val_scale*(((val_enthalpy)*(val_dRhodh) + (val_density))*proj_vel);
+  val_Proj_Jac_Tensor[nDim+1][nDim+1] = val_scale * MassFlux;
   for (iDim = 0; iDim < nDim; iDim++) {
     val_Proj_Jac_Tensor[nDim+1][iDim+1] = 0.0;
     val_Proj_Jac_Tensor[iDim+1][nDim+1] = 0.0;

SU2_CFD/src/solvers/CIncEulerSolver.cpp

@@ -1053,8 +1053,9 @@ void CIncEulerSolver::CommonPreprocessing(CGeometry *geometry, CSolver **solver_
     SU2_OMP_SAFE_GLOBAL_ACCESS(GetOutlet_Properties(geometry, config, iMesh, Output);)
   }
 
-  /*--- Reset flag for strong BCs. ---*/
-  if (pressure_based) {
+  /*--- Reset flag for strong BCs. Not on the Output call at the end of the momentum iteration: the
+   *    flags set by the BCs are used afterwards by the momentum coefficients and the corrections. ---*/
+  if (pressure_based && !Output) {
     SU2_OMP_FOR_STAT(omp_chunk_size)
     for (unsigned long iPoint = 0; iPoint < nPointDomain; iPoint++)
       nodes->ResetStrongBC(iPoint);
@@ -2139,15 +2140,24 @@ void CIncEulerSolver::PrepareImplicitIteration(CGeometry *geometry, CSolver**, C
 
   struct IncPrec {
     const CIncEulerSolver* solver;
-    const bool active;
+    const bool active = true;
     su2activematrix matrix;
 
-    IncPrec(const CIncEulerSolver* s, unsigned short nVar) : solver(s), active(!s->pressure_based) {
+    IncPrec(const CIncEulerSolver* s, unsigned short nVar) : solver(s) {
       matrix.resize(nVar,nVar);
     }
 
     FORCEINLINE const su2activematrix& operator() (const CConfig* config, unsigned long iPoint, su2double delta) {
-      solver->SetPreconditioner(config, iPoint, delta, matrix);
+      if (solver->pressure_based) {
+        /*--- The Jacobian of the pressure-based solver is per unit velocity and enthalpy, the
+         *    pseudo-time term is rho*V/dt (the continuity row is deleted below). ---*/
+        matrix = su2double(0.0);
+        matrix(0,0) = delta;
+        for (auto iVar = 1ul; iVar < matrix.rows(); iVar++)
+          matrix(iVar,iVar) = solver->nodes->GetDensity(iPoint) * delta;
+      } else {
+        solver->SetPreconditioner(config, iPoint, delta, matrix);
+      }
       return matrix;
     }
 
@@ -2457,17 +2467,24 @@ void CIncEulerSolver::BC_Far_Field(CGeometry *geometry, CSolver **solver_contain
 
       if (inflow) {
 
-        /*--- Set this face as an inlet via a strong BC. ---*/
+        /*--- Set this face as an inlet via a strong BC, the free-stream velocity and enthalpy are imposed. ---*/
+
+        nodes->SetVelocity_Old(iPoint, V_infty+prim_idx.Velocity());
+        if (config->GetEnergy_Equation())
+          nodes->SetSolution_Old(iPoint, nDim+1, V_infty[prim_idx.Enthalpy()]);
 
         LinSysRes.SetBlock_Zero(iPoint);
 
         /*--- Mark as a strong BC which is important for deciding if a velocity correction should be applied ---*/
 
         nodes->SetStrongBC(iPoint);
 
-        if (implicit)
+        if (implicit) {
           for (iDim = 0; iDim < nDim; iDim++)
             Jacobian.DeleteValsRowi(iPoint, iDim+1);
+          if (config->GetEnergy_Equation())
+            Jacobian.DeleteValsRowi(iPoint, nDim+1);
+        }
 
       } else {
 
@@ -2738,13 +2755,29 @@ void CIncEulerSolver::BC_Inlet(CGeometry *geometry, CSolver **solver_container,
 
       nodes->SetVelocity_Old(iPoint,V_inlet+prim_idx.Velocity());
 
+      /*--- Same for the enthalpy, from the inlet temperature. ---*/
+
+      const bool inlet_enthalpy = config->GetEnergy_Equation() &&
+          !(config->GetKind_Species_Model() == SPECIES_MODEL::FLAMELET &&
+            config->GetFlamelet_Enthalpy_BC() == FLAMELET_ENTHALPY_BC::FLOW_MARKERS);
+      if (inlet_enthalpy) {
+        const su2double* scalar_inlet = nullptr;
+        if (species_model) scalar_inlet = config->GetInlet_SpeciesVal(config->GetMarker_All_TagBound(val_marker));
+        CFluidModel* auxFluidModel = solver_container[FLOW_SOL]->GetFluidModel();
+        auxFluidModel->SetTDState_T(V_inlet[prim_idx.Temperature()], scalar_inlet);
+        nodes->SetSolution_Old(iPoint, nDim+1, auxFluidModel->GetEnthalpy());
+      }
+
       LinSysRes.SetBlock_Zero(iPoint);
 
       if (pressure_based) nodes->SetStrongBC(iPoint);
 
-      if (implicit)
+      if (implicit) {
         for (iDim = 0; iDim < nDim; iDim++)
           Jacobian.DeleteValsRowi(iPoint, iDim+1);
+        if (inlet_enthalpy)
+          Jacobian.DeleteValsRowi(iPoint, nDim+1);
+      }
 
     } else {
 
@@ -3716,11 +3749,9 @@ void CIncEulerSolver::ComputeEdgeMassFluxesRhieChow(CGeometry *geometry, CSolver
 
     CorrectPressureGradient(GradPressure_f, GradPressure_avg, nodes->GetPressure(iPoint), nodes->GetPressure(jPoint), Edge_Vector, dist_ij_2);
 
-    /*--- Linearly interpolated coefficient. A point under a strong velocity BC has no momentum
-    coefficient, so the edge uses that of its other node. ---*/
+    /*--- Linearly interpolated coefficient. ---*/
 
-    Coeff_Mom = 0.5*(poisson_nodes->GetMomCoeff(nodes->GetStrongBC(iPoint) ? jPoint : iPoint) +
-                     poisson_nodes->GetMomCoeff(nodes->GetStrongBC(jPoint) ? iPoint : jPoint));
+    Coeff_Mom = 0.5*(poisson_nodes->GetMomCoeff(iPoint) + poisson_nodes->GetMomCoeff(jPoint));
 
     /*--- Initialize mass flux ---*/
 
@@ -3868,8 +3899,7 @@ void CIncEulerSolver::ApplyPressureVelocityCorrection(CGeometry *geometry, CSolv
     for (iDim = 0; iDim < nDim; iDim++) {
 
       su2double MassFluxCorrection =
-          -0.5 * (poisson_nodes->GetMomCoeff(nodes->GetStrongBC(iPoint) ? jPoint : iPoint) +
-                  poisson_nodes->GetMomCoeff(nodes->GetStrongBC(jPoint) ? iPoint : jPoint)) * GradPressure_f[iDim];
+          -0.5 * (poisson_nodes->GetMomCoeff(iPoint) + poisson_nodes->GetMomCoeff(jPoint)) * GradPressure_f[iDim];
 
       /*--- 2nd piso correction term (HbyA') --- (TODO: this is zero for the first correction and can thus also be skipped) ---*/
 

SU2_CFD/src/solvers/CPoissonSolver.cpp

@@ -168,14 +168,11 @@ void CPoissonSolver::SetMomCoeff(CGeometry *geometry, CSolver **solver_container
       su2double Vol = geometry->nodes->GetVolume(iPoint);
 
       /*--- The momentum equation is not assembled at a strong velocity BC, DeleteValsRowi zeroes
-       * the row and writes 1.0 on the diagonal, so there is no A_p to read. Nothing consumes the
-       * value stored here: edges touching the point take the coefficient of their other node, the
-       * velocity correction is overwritten in the boundary loop, and HbyA scales it by a numerator
-       * that is identically zero. Store a finite placeholder and skip the corrections below, which
-       * divide by zero for a transient removal factor of 1. ---*/
+       * the row and writes 1.0 on the diagonal, so there is no A_p to read. The coefficient of
+       * these points is set below from their neighbours, zero marks them until then. ---*/
 
       if (flow_nodes->GetStrongBC(iPoint)) {
-        nodes->SetMomCoeff(iPoint, Vol * flow_nodes->GetDensity(iPoint));
+        nodes->SetMomCoeff(iPoint, 0.0);
         continue;
       }
 
@@ -194,11 +191,13 @@ void CPoissonSolver::SetMomCoeff(CGeometry *geometry, CSolver **solver_container
         }
       }
 
-      /*--- Add simplec neighbour contributions and optional time dependent term. ---*/
+      /*--- Add simplec neighbour contributions and optional time dependent term. The off-diagonal
+       * entries are -a_nb, so A_p + Sum_A_nb is a_P - sum(a_nb), which is only the pseudo-time term
+       * V/dt: the removal of the transient term must not be used with SIMPLEC. ---*/
 
       su2double delT = flow_nodes->GetDelta_Time(iPoint);
 
-      su2double CorrectedA_p = A_p - Sum_A_nb - config->GetSIMPLE_Options().Transient_Term_Removal_Factor * (Vol / delT);
+      su2double CorrectedA_p = A_p + Sum_A_nb - config->GetSIMPLE_Options().Transient_Term_Removal_Factor * (Vol / delT);
 
       /*--- Invert the momentum coefficient to 1/a_p and scale by the volume and density so it can be used as diffusion coefficient in the poisson eq ---*/
 
@@ -225,7 +224,32 @@ void CPoissonSolver::SetMomCoeff(CGeometry *geometry, CSolver **solver_container
     */
   }
 
-  /*--- Insert MPI call here. ---*/
+  InitiateComms(geometry, config, MPI_QUANTITIES::MOM_COEFF);
+  CompleteComms(geometry, config, MPI_QUANTITIES::MOM_COEFF);
+
+  /*--- Points under a strong velocity BC take the average coefficient of their neighbours that have
+   * a momentum equation, which also sets the edges between two such points. The flags are only known
+   * on this rank's points, the zero set above marks the strong points of the other ranks. ---*/
+
+  SU2_OMP_FOR_STAT(omp_chunk_size)
+  for (unsigned long iPoint = 0; iPoint < nPointDomain; iPoint++) {
+    if (!flow_nodes->GetStrongBC(iPoint)) continue;
+
+    su2double Sum_Coeff = 0.0;
+    unsigned short nCoeff = 0;
+    for (auto jPoint : geometry->nodes->GetPoints(iPoint)) {
+      if (!flow_nodes->GetStrongBC(jPoint) && nodes->GetMomCoeff(jPoint) > 0.0) {
+        Sum_Coeff += nodes->GetMomCoeff(jPoint);
+        nCoeff++;
+      }
+    }
+    if (nCoeff > 0)
+      nodes->SetMomCoeff(iPoint, Sum_Coeff / nCoeff);
+    else
+      nodes->SetMomCoeff(iPoint, 0.0);
+  }
+  END_SU2_OMP_FOR
+
   InitiateComms(geometry, config, MPI_QUANTITIES::MOM_COEFF);
   CompleteComms(geometry, config, MPI_QUANTITIES::MOM_COEFF);
 }
@@ -465,6 +489,35 @@ void CPoissonSolver::Source_Residual(CGeometry *geometry, CSolver **solver_conta
 
 }
 
+void CPoissonSolver::SetResidual_DualTime(CGeometry *geometry, CSolver **solver_container, CConfig *config,
+                                          unsigned short iRKStep, unsigned short iMesh,
+                                          unsigned short RunTime_EqSystem) {
+  SU2_ZONE_SCOPED
+
+  /*--- Unsteady variable density: continuity is V d(rho)/dt + sum(m_f) = 0. The density of the old
+   *    time levels is the one the flow solver computes from their enthalpy and species
+   *    (CIncEulerSolver::RecomputeDensity_time_n). With constant density the term is zero. ---*/
+
+  if (!config->GetVariable_Density_Model()) return;
+
+  const bool second_order = (config->GetTime_Marching() == TIME_MARCHING::DT_STEPPING_2ND);
+  const su2double dt = config->GetDelta_UnstTimeND();
+  const auto* flow_vars = su2staticcast_p<const CFlowVariable*>(solver_container[FLOW_SOL]->GetNodes());
+
+  SU2_OMP_FOR_STAT(omp_chunk_size)
+  for (unsigned long iPoint = 0; iPoint < nPointDomain; ++iPoint) {
+    const su2double rho_np1 = flow_vars->GetDensity(iPoint);
+    const su2double rho_n = flow_vars->GetDensity_time_n(iPoint);
+    su2double dRhodt = (rho_np1 - rho_n) / dt;
+    if (second_order) {
+      const su2double rho_nm1 = flow_vars->GetDensity_time_n1(iPoint);
+      dRhodt = (3.0 * rho_np1 - 4.0 * rho_n + rho_nm1) / (2.0 * dt);
+    }
+    LinSysRes(iPoint, 0) += geometry->nodes->GetVolume(iPoint) * dRhodt;
+  }
+  END_SU2_OMP_FOR
+}
+
 void CPoissonSolver::ImplicitEuler_Iteration(CGeometry *geometry, CSolver **solver_container, CConfig *config) {
   SU2_ZONE_SCOPED
 

TestCases/hybrid_regression.py

@@ -369,7 +369,7 @@ def main():
     inc_euler_naca0012_pb.cfg_dir   = "incomp_euler/naca0012"
     inc_euler_naca0012_pb.cfg_file  = "incomp_pb_NACA0012.cfg"
     inc_euler_naca0012_pb.test_iter = 20
-    inc_euler_naca0012_pb.test_vals = [-4.454818, -4.784246, 0.427448, 0.012083]
+    inc_euler_naca0012_pb.test_vals = [-4.454835, -4.784268, 0.427448, 0.012084]
     test_list.append(inc_euler_naca0012_pb)
 
     # Laminar cylinder, pressure-based. Viscous counterpart of the hydrofoil above, so the
@@ -378,7 +378,7 @@ def main():
     inc_lam_cylinder_pb.cfg_dir   = "incomp_navierstokes/cylinder"
     inc_lam_cylinder_pb.cfg_file  = "incomp_pb_cylinder.cfg"
     inc_lam_cylinder_pb.test_iter = 10
-    inc_lam_cylinder_pb.test_vals = [-3.486075, -3.777632, 0.012229, 6.178719]
+    inc_lam_cylinder_pb.test_vals = [-3.486562, -3.776813, 0.012619, 6.205610]
     test_list.append(inc_lam_cylinder_pb)
 
     # Laminar sphere, pressure-based. The only 3D pressure-based case, and the one where the
@@ -387,15 +387,15 @@ def main():
     inc_lam_sphere_pb.cfg_dir   = "incomp_navierstokes/sphere"
     inc_lam_sphere_pb.cfg_file  = "pb_sphere.cfg"
     inc_lam_sphere_pb.test_iter = 9
-    inc_lam_sphere_pb.test_vals = [-6.092084, -2.305039, -2.479072, -2.548968, 0.191798, 170.793204, -6.134920]
+    inc_lam_sphere_pb.test_vals = [-6.104969, -2.307895, -2.481829, -2.570423, 0.189098, 170.764342, -6.128569]
     test_list.append(inc_lam_sphere_pb)
 
     # Heated cylinder, pressure-based, coupled energy equation and variable density.
     inc_poly_cylinder_pb          = TestCase('inc_poly_cylinder_pb')
     inc_poly_cylinder_pb.cfg_dir   = "incomp_navierstokes/cylinder"
     inc_poly_cylinder_pb.cfg_file  = "pb_poly_cylinder.cfg"
     inc_poly_cylinder_pb.test_iter = 20
-    inc_poly_cylinder_pb.test_vals = [-13.483266, 0.541350, 0.005972, 17.020220, -8927.600000]
+    inc_poly_cylinder_pb.test_vals = [-9.564374, -3.013394, 0.000161, 1.778555, -172.280000]
     test_list.append(inc_poly_cylinder_pb)
 
     # C-D nozzle with pressure inlet and mass flow outlet