Fix damping mass outlet (#2516)

* modifications for damping factor

Co-authored-by: Pedro Gomes <[email protected]>

Author: bigfooted

Common/include/CConfig.hpp

@@ -1092,6 +1092,7 @@ class CConfig {
   su2double *FreeStreamTurboNormal;     /*!< \brief Direction to initialize the flow in turbomachinery computation */
   su2double Restart_Bandwidth_Agg;      /*!< \brief The aggregate of the bandwidth for writing binary restarts (to be averaged later). */
   su2double Max_Vel2;                   /*!< \brief The maximum velocity^2 in the domain for the incompressible preconditioner. */
+  su2double RangePressure[2];           /*!< \brief The pressure difference pmax-pmin in the domain for the target mass flow rate scaling. */
   bool topology_optimization;           /*!< \brief If the structural solver should consider a variable density field to penalize element stiffness. */
   string top_optim_output_file;         /*!< \brief File to where the derivatives w.r.t. element densities will be written to. */
   su2double simp_exponent;              /*!< \brief Exponent for the density-based stiffness penalization of the SIMP method. */
@@ -9222,6 +9223,18 @@ class CConfig {
    */
   void SetMax_Vel2(su2double val_max_vel2) { Max_Vel2 = val_max_vel2; }
 
+  /*!
+   * \brief Get the maximum pressure (pmax - pmin) in the domain.
+   * \return Value of the maximum pressure in the domain.
+   */
+  su2double GetRangePressure(int minmax) const { return RangePressure[minmax]; }
+
+  /*!
+   * \brief Set the maximum pressure in the domain.
+   * \param[in] Value of the maximum pressure in the domain.
+   */
+  void SetRangePressure(su2double val_dp_min,su2double val_dp_max) { RangePressure[0] = val_dp_min;RangePressure[1]=val_dp_max; }
+
   /*!
    * \brief Get the maximum velocity^2 in the domain for the incompressible preconditioner.
    * \return Value of the maximum velocity^2 in the domain for the incompressible preconditioner.

SU2_CFD/include/solvers/CIncEulerSolver.hpp

@@ -98,6 +98,18 @@ class CIncEulerSolver : public CFVMFlowSolverBase<CIncEulerVariable, ENUM_REGIME
                          CConfig *config,
                          unsigned short iMesh);
 
+  /*!
+   * \brief Update the pressure range P_max - P_min for target mass flow rate.
+   * \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] iMesh - current mesh level for the multigrid.
+   */
+  void SetRangePressure(CGeometry *geometry,
+                         CSolver **solver_container,
+                         CConfig *config,
+                         unsigned short iMesh);
+
   /*!
    * \brief A virtual member.
    */

SU2_CFD/src/solvers/CIncEulerSolver.cpp

@@ -985,6 +985,10 @@ void CIncEulerSolver::CommonPreprocessing(CGeometry *geometry, CSolver **solver_
 
   SetBeta_Parameter(geometry, solver_container, config, iMesh);
 
+  /*--- Update the pressure range in the domain for target outflow mass flow rate. ---*/
+
+  SetRangePressure(geometry, solver_container, config, iMesh);
+
   /*--- Compute properties needed for mass flow BCs. ---*/
 
   if (outlet) {
@@ -1974,6 +1978,46 @@ void CIncEulerSolver::SetBeta_Parameter(CGeometry *geometry, CSolver **solver_co
 
 }
 
+void CIncEulerSolver::SetRangePressure(CGeometry *geometry, CSolver **solver_container,
+                                        CConfig *config, unsigned short iMesh) {
+  static su2double MinP, MaxP;
+
+  if (iMesh == MESH_0) {
+    SU2_OMP_MASTER {
+      MinP = std::numeric_limits<su2double>::max();
+      MaxP = std::numeric_limits<su2double>::lowest();
+    }
+    END_SU2_OMP_MASTER
+    su2double minP = std::numeric_limits<su2double>::max();
+    su2double maxP = std::numeric_limits<su2double>::lowest();
+
+    SU2_OMP_FOR_STAT(omp_chunk_size)
+    for (auto iPoint = 0ul; iPoint < nPoint; iPoint++) {
+      minP = min(minP, nodes->GetPressure(iPoint));
+      maxP = max(maxP, nodes->GetPressure(iPoint));
+    }
+    END_SU2_OMP_FOR
+
+    SU2_OMP_CRITICAL {
+      MinP = min(MinP, minP);
+      MaxP = max(MaxP, maxP);
+    }
+    END_SU2_OMP_CRITICAL
+
+    BEGIN_SU2_OMP_SAFE_GLOBAL_ACCESS
+    {
+      minP = MinP;
+      SU2_MPI::Allreduce(&minP, &MinP, 1, MPI_DOUBLE, MPI_MAX, SU2_MPI::GetComm());
+      maxP = MaxP;
+      SU2_MPI::Allreduce(&maxP, &MaxP, 1, MPI_DOUBLE, MPI_MAX, SU2_MPI::GetComm());
+
+      config->SetRangePressure(MinP,MaxP);
+    }
+    END_SU2_OMP_SAFE_GLOBAL_ACCESS
+  }
+
+}
+
 void CIncEulerSolver::SetPreconditioner(const CConfig *config, unsigned long iPoint,
                                         su2double delta, su2activematrix& Preconditioner) const {
 
@@ -2556,10 +2600,14 @@ void CIncEulerSolver::BC_Outlet(CGeometry *geometry, CSolver **solver_container,
 
         dP = 0.5*Density_Avg*(mDot_Old*mDot_Old - mDot_Target*mDot_Target)/((Density_Avg*Area_Outlet)*(Density_Avg*Area_Outlet));
 
-        /*--- Update the new outlet pressure. Note that we use damping
-         here to improve stability/convergence. ---*/
+        su2double P_domain_min = config->GetRangePressure(0);
+
+        /*--- Do not relax when dP is relatively small compared to the pressure range dp = (P-P_min). ---*/
+
+        if (abs(dP) < abs(Damping * (P_domain-P_domain_min)))
+          Damping = 1.0;
 
-        P_Outlet = P_domain + Damping*dP;
+        P_Outlet = P_domain + Damping * dP;
 
         /*--- The pressure is prescribed at the outlet. ---*/
 

TestCases/hybrid_regression.py

@@ -352,7 +352,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.595923, -5.820518, -0.018365, 0.126411]
+    inc_nozzle.test_vals         = [-6.593521, -5.830706, -0.009062, 0.126050]
     inc_nozzle.test_vals_aarch64 = [-5.624385, -4.988472, -0.000096, 0.137032]
     test_list.append(inc_nozzle)
 

TestCases/parallel_regression.py

@@ -566,7 +566,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.576818, -5.796867, -0.003141, 0.126481]
+    inc_nozzle.test_vals = [-6.407323, -5.715668, -0.003225, 0.126214]
     test_list.append(inc_nozzle)
 
     # Laminar wall mounted cylinder, Euler walls, cylinder wall diagonally split

TestCases/serial_regression.py

@@ -379,7 +379,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.623301, -5.844127, -0.023181, 0.126370]
+    inc_nozzle.test_vals = [-5.394788, -4.869896, -0.021578, 0.125704]
     test_list.append(inc_nozzle)
 
     #############################