Fix shared wall-inlet corner node for compressible

Common/src/CConfig.cpp

@@ -5538,16 +5538,12 @@ void CConfig::SetPostprocessing(SU2_COMPONENT val_software, unsigned short val_i
   /*--- Define some variables for flamelet model. ---*/
   if (Kind_Species_Model == SPECIES_MODEL::FLAMELET) {
     /*--- The controlling variables are progress variable, total enthalpy, and optionally mixture fraction ---*/
-    //n_control_vars = nSpecies - n_user_scalars;
     if (n_control_vars != (nSpecies - n_user_scalars))
-      SU2_MPI::Error("Number of initial species incompatbile with number of controlling variables and user scalars.", CURRENT_FUNCTION);
+      SU2_MPI::Error("Number of initial species incompatible with number of controlling variables and user scalars.", CURRENT_FUNCTION);
     /*--- We can have additional user defined transported scalars ---*/
     n_scalars = n_control_vars + n_user_scalars;
   }
 
-  if (Kind_Regime == ENUM_REGIME::COMPRESSIBLE && GetBounded_Scalar()) {
-    SU2_MPI::Error("BOUNDED_SCALAR discretization can only be used for incompressible problems.", CURRENT_FUNCTION);
-  }
 }
 
 void CConfig::SetMarkers(SU2_COMPONENT val_software) {

SU2_CFD/include/solvers/CSpeciesSolver.hpp

@@ -153,6 +153,26 @@ class CSpeciesSolver : public CScalarSolver<CSpeciesVariable> {
    * (like BC_Isothermal_Wall) are implemented the respective CHeatSolver implementations can act as a  starting
    * point ---*/
 
+  /*!
+   * \brief Generic implementation of the isothermal wall also covering CHT cases,
+   * for which the wall temperature is given by GetConjugateHeatVariable.
+   */
+  void BC_Isothermal_Wall_Generic(CGeometry* geometry, CSolver** solver_container, CNumerics* conv_numerics,
+                                  CNumerics* visc_numerics, CConfig* config, unsigned short val_marker,
+                                  bool cht_mode = false);
+
+  /*!
+   * \brief Impose the scalar wall boundary condition.
+   * \param[in] geometry - Geometrical definition of the problem.
+   * \param[in] solver_container - Container vector with all the solutions.
+   * \param[in] conv_numerics - Description of the numerical method.
+   * \param[in] visc_numerics - Description of the numerical method.
+   * \param[in] config - Definition of the particular problem.
+   * \param[in] val_marker - Surface marker where the boundary condition is applied.
+   */
+  void BC_Isothermal_Wall(CGeometry* geometry, CSolver** solver_container, CNumerics* conv_numerics,
+                          CNumerics* visc_numerics, CConfig* config, unsigned short val_marker) override;
+
   /*!
    * \brief Source term computation for axisymmetric flow.
    * \param[in] geometry - Geometrical definition of the problem.

SU2_CFD/src/solvers/CEulerSolver.cpp

@@ -341,6 +341,10 @@
 
   CommunicateInitialState(geometry, config);
 
+  /*--- Sizing edge mass flux array ---*/
+  if (config->GetBounded_Scalar())
+    EdgeMassFluxes.resize(geometry->GetnEdge()) = su2double(0.0);
+
   /*--- Add the solver name.. ---*/
   SolverName = "C.FLOW";
 
@@ -1772,7 +1776,7 @@
   const bool muscl            = (config->GetMUSCL_Flow() && (iMesh == MESH_0));
   const bool limiter          = (config->GetKind_SlopeLimit_Flow() != LIMITER::NONE);
   const bool van_albada       = (config->GetKind_SlopeLimit_Flow() == LIMITER::VAN_ALBADA_EDGE);
-
+  //const bool bounded_scalar   =  config->GetBounded_Scalar();
   /*--- Non-physical counter. ---*/
   unsigned long counter_local = 0;
   SU2_OMP_MASTER
@@ -1937,6 +1941,10 @@
 
     auto residual = numerics->ComputeResidual(config);
 
+    // this does not seem to be necessary
+    cout << "bounded" << endl;
+    //if (bounded_scalar) EdgeMassFluxes[iEdge] = residual[0];
+
     /*--- Set the final value of the Roe dissipation coefficient ---*/
 
     if ((kind_dissipation != NO_ROELOWDISS) && (MGLevel != MESH_0)) {
@@ -5067,8 +5075,8 @@
       su2double Velocity_i[MAXNDIM]={0};
       for (auto iDim=0u; iDim < nDim; iDim++)
         Velocity_i[iDim] = nodes->GetVelocity(iPoint,iDim);
-      
-      const auto Velocity2_i = GeometryToolbox::SquaredNorm(nDim, Velocity_i);  
+
+      const auto Velocity2_i = GeometryToolbox::SquaredNorm(nDim, Velocity_i);
 
       const auto Density_i = nodes->GetDensity(iPoint);
 
@@ -5090,7 +5098,7 @@
                 T_Total{0},  P_Total{0},  Density_e{0}, StaticEnthalpy_e{0}, StaticEnergy_e{0};
 
       su2double Velocity2_e{0}, NormalVelocity{0}, TangVelocity{0}, VelMag_e{0};
-      su2double Velocity_e[MAXNDIM] = {0}; 
+      su2double Velocity_e[MAXNDIM] = {0};
       const su2double * Flow_Dir, * Mach;
       switch(config->GetKind_Data_Riemann(Marker_Tag))
       {
@@ -5146,10 +5154,10 @@
         GetFluidModel()->SetTDState_PT(P_static, T_static);
 
         /* --- Compute the boundary state u_e --- */
-        for (auto iDim = 0u; iDim < nDim; iDim++) 
+        for (auto iDim = 0u; iDim < nDim; iDim++)
           Velocity_e[iDim] = Mach[iDim]*GetFluidModel()->GetSoundSpeed();
 
-        Velocity2_e = GeometryToolbox::SquaredNorm(nDim, Velocity_e);  
+        Velocity2_e = GeometryToolbox::SquaredNorm(nDim, Velocity_e);
         Density_e = GetFluidModel()->GetDensity();
         StaticEnergy_e = GetFluidModel()->GetStaticEnergy();
         Energy_e = StaticEnergy_e + 0.5 * Velocity2_e;
@@ -5176,7 +5184,7 @@
         for (auto iDim = 0u; iDim < nDim; iDim++)
           Velocity_e[iDim] = Mach[iDim]*GetFluidModel()->GetSoundSpeed();
 
-        Velocity2_e = GeometryToolbox::SquaredNorm(nDim, Velocity_e);  
+        Velocity2_e = GeometryToolbox::SquaredNorm(nDim, Velocity_e);
         Density_e = GetFluidModel()->GetDensity();
         StaticEnergy_e = GetFluidModel()->GetStaticEnergy();
         Energy_e = StaticEnergy_e + 0.5 * Velocity2_e;
@@ -5208,10 +5216,10 @@
 
         /* --- Compute the boundary state u_e --- */
         GetFluidModel()->SetTDState_Prho(Pressure_e, Density_e);
-        for (auto iDim = 0u; iDim < nDim; iDim++) 
+        for (auto iDim = 0u; iDim < nDim; iDim++)
           Velocity_e[iDim] = Velocity_i[iDim];
 
-        Velocity2_e = GeometryToolbox::SquaredNorm(nDim, Velocity_e);  
+        Velocity2_e = GeometryToolbox::SquaredNorm(nDim, Velocity_e);
         Energy_e = GetFluidModel()->GetStaticEnergy() + 0.5*Velocity2_e;
         break;
 
@@ -5235,7 +5243,7 @@
       }
 
       /*--- Flow eigenvalues, boundary state u_e and u_i ---*/
-      su2double Lambda_i[MAXNVAR] = {0}, u_e[MAXNVAR] = {0}, u_i[MAXNVAR]={0}, u_b[MAXNVAR]={0}, 
+      su2double Lambda_i[MAXNVAR] = {0}, u_e[MAXNVAR] = {0}, u_i[MAXNVAR]={0}, u_b[MAXNVAR]={0},
                 dw[MAXNVAR]={0};
       u_e[0] = Density_e;
       u_i[0] = Density_i;
@@ -5302,7 +5310,7 @@
       LinSysRes.AddBlock(iPoint, Residual);
 
       if (implicit) {
-        su2double **Jacobian_b = new su2double*[nVar], 
+        su2double **Jacobian_b = new su2double*[nVar],
                   **Jacobian_i = new su2double*[nVar];
         su2double **DubDu = new su2double*[nVar];
         for (auto iVar = 0u; iVar < nVar; iVar++){
@@ -5333,7 +5341,7 @@
         if (dynamic_grid){
           const auto gridVel = geometry->nodes->GetGridVel(iPoint);
           const auto projVelocity = GeometryToolbox::DotProduct(nDim, gridVel, Normal);
-          for (auto iVar = 0u; iVar < nVar; iVar++) 
+          for (auto iVar = 0u; iVar < nVar; iVar++)
             Jacobian_b[iVar][iVar] -= projVelocity;
         }
 
@@ -5394,7 +5402,7 @@
         visc_numerics->SetPrimVarGradient(nodes->GetGradient_Primitive(iPoint), nodes->GetGradient_Primitive(iPoint));
 
         /*--- Secondary variables ---*/
-        
+
         auto S_domain = nodes->GetSecondary(iPoint);
 
         /*--- Compute secondary thermodynamic properties (partial derivatives...) ---*/
@@ -6880,10 +6888,11 @@
                             CConfig *config, unsigned short val_marker) {
   unsigned short iDim;
   unsigned long iVertex, iPoint;
-  su2double P_Total, T_Total, Velocity[MAXNDIM], Velocity2, H_Total, Temperature, Riemann,
+  su2double P_Total, T_Total, Velocity[MAXNDIM], Velocity2, H_Total, Riemann,
   Pressure, Density, Energy, Flow_Dir[MAXNDIM], Mach2, SoundSpeed2, SoundSpeed_Total2, Vel_Mag,
   alpha, aa, bb, cc, dd, Area, UnitNormal[MAXNDIM], Normal[MAXNDIM];
   su2double *V_inlet, *V_domain;
+  su2double Temperature = 288.15;
 
   const bool implicit = (config->GetKind_TimeIntScheme() == EULER_IMPLICIT);
   const su2double Two_Gamma_M1 = 2.0 / Gamma_Minus_One;
@@ -6921,7 +6930,7 @@
 
       V_domain = nodes->GetPrimitive(iPoint);
 
-      /*--- Build the fictitious intlet state based on characteristics ---*/
+      /*--- Build the fictitious inlet state based on characteristics ---*/
 
 
       /*--- Subsonic inflow: there is one outgoing characteristic (u-c),
@@ -7096,8 +7105,8 @@
           if (tkeNeeded) Energy += GetTke_Inf();
 
           /*--- Primitive variables, using the derived quantities ---*/
-
-          V_inlet[0] = Pressure / ( Gas_Constant * Density);
+          Temperature = Pressure / ( Gas_Constant * Density);
+          V_inlet[0] = Temperature;
           for (iDim = 0; iDim < nDim; iDim++)
             V_inlet[iDim+1] = Vel_Mag*Flow_Dir[iDim];
           V_inlet[nDim+1] = Pressure;
@@ -7111,6 +7120,65 @@
           break;
       }
 
+      /*--- Check if the inlet node is shared with a viscous wall. ---*/
+
+      if (geometry->nodes->GetViscousBoundary(iPoint)) {
+
+        switch (Kind_Inlet) {
+
+          /*--- Total properties have been specified at the inlet. ---*/
+
+          case INLET_TYPE::TOTAL_CONDITIONS: {
+
+            /*--- Impose the wall velocity from the interior wall node. ---*/
+
+            Velocity2 = 0.0;
+            for (iDim = 0; iDim < nDim; iDim++) {
+              V_inlet[iDim+1] = nodes->GetVelocity(iPoint,iDim);
+              Velocity2 += V_inlet[iDim+1] * V_inlet[iDim+1];
+            }
+
+            /*--- Match the pressure, density and energy at the wall. ---*/
+
+            Pressure = nodes->GetPressure(iPoint);
+            Density = Pressure / (Gas_Constant * Temperature);
+            Energy = Pressure / (Density*Gamma_Minus_One) + 0.5 * Velocity2;
+            if (tkeNeeded) Energy += GetTke_Inf();
+
+            V_inlet[nDim+1] = Pressure;
+            V_inlet[nDim+2] = Density;
+            V_inlet[nDim+3] = Energy + Pressure/Density;
+            break;
+          }
+
+          case INLET_TYPE::MASS_FLOW: {
+
+            /*--- Impose the wall velocity from the interior wall node. ---*/
+
+            Velocity2 = 0.0;
+            for (iDim = 0; iDim < nDim; iDim++) {
+              V_inlet[iDim+1] = nodes->GetVelocity(iPoint,iDim);
+              Velocity2 += V_inlet[iDim+1] * V_inlet[iDim+1];
+            }
+
+            Pressure = nodes->GetPressure(iPoint);
+            Density = nodes->GetDensity(iPoint);
+
+            Energy = Pressure / (Density * Gamma_Minus_One)  + 0.5 * Velocity2;
+
+            if (tkeNeeded) Energy += GetTke_Inf();
+
+            V_inlet[nDim+3] = Energy + Pressure/Density;
+
+            break;
+          }
+
+          default:
+            SU2_MPI::Error("Unsupported INLET_TYPE.", CURRENT_FUNCTION);
+            break;
+
+        }
+      }
       /*--- Set various quantities in the solver class ---*/
 
       conv_numerics->SetPrimitive(V_domain, V_inlet);
@@ -8795,9 +8863,9 @@
 
   const auto nSpanWiseSections = config->GetnSpanWiseSections();
   const auto iZone = config->GetiZone();
-  
+
   for (auto iSpan= 0u; iSpan < nSpanWiseSections; iSpan++){
-    su2double TotalAreaVelocity[MAXNDIM]={0.0}, 
+    su2double TotalAreaVelocity[MAXNDIM]={0.0},
               TotalAreaPressure{0},
               TotalAreaDensity{0};
     for (auto iMarker = 0u; iMarker < config->GetnMarker_All(); iMarker++){
@@ -8987,7 +9055,7 @@
         TotalAreaVelocity[iDim] += Area*Velocity[iDim];
         TotalMassVelocity[iDim] += Area*(Density*TurboVelocity[0] )*Velocity[iDim];
       }
-        
+
       TotalFluxes[0]      += Area*(Density*TurboVelocity[0]);
       TotalFluxes[1]      += Area*(Density*TurboVelocity[0]*TurboVelocity[0] + Pressure);
       for (auto iDim = 2; iDim < nDim+1; iDim++)
@@ -9036,7 +9104,7 @@
                   UpdateTotalQuantities(iMarker, jSpan, iVertex);
                 }
               }
-            } 
+            }
 
           } // marker_flag match
         } // iMarkerTP match
@@ -9251,7 +9319,7 @@
               OmegaOut[iMarkerTP - 1][iSpan]    = AverageOmega[iMarker][iSpan];
               NuOut[iMarkerTP - 1][iSpan]       = AverageNu[iMarker][iSpan];
             }
-            
+
             auto TurboVel = (marker_flag == INFLOW) ? TurboVelocityIn[iMarkerTP - 1][iSpan] : TurboVelocityOut[iMarkerTP - 1][iSpan];
 
             if (performance_average_process == MIXEDOUT) {
@@ -9281,7 +9349,7 @@
               if(config->GetKind_Data_Giles(Marker_Tag) == RADIAL_EQUILIBRIUM){
                 RadialEquilibriumPressure[iMarker][nSpanWiseSections/2] = config->GetGiles_Var1(Marker_Tag)/config->GetPressure_Ref();
               }
-            } 
+            }
             for (auto iSpan= nSpanWiseSections/2; iSpan < nSpanWiseSections-1; iSpan++){
               const auto Radius2    = geometry->GetTurboRadius(iMarker,iSpan+1);
               const auto Radius1    = geometry->GetTurboRadius(iMarker,iSpan);
@@ -9325,9 +9393,9 @@
 
     su2double vel[MAXNDIM] = {0};
     vel[0]  = (val_Averaged_Flux[1] - pressure_mix) / val_Averaged_Flux[0];
-    for (auto iDim = 1u; iDim < nDim; iDim++) 
+    for (auto iDim = 1u; iDim < nDim; iDim++)
       vel[iDim]  = val_Averaged_Flux[iDim+1] / val_Averaged_Flux[0];
-    
+
 
     const su2double velsq = GeometryToolbox::DotProduct(nDim, vel, vel);
 

SU2_CFD/src/solvers/CSpeciesSolver.cpp

@@ -345,7 +345,7 @@ void CSpeciesSolver::BC_Inlet(CGeometry* geometry, CSolver** solver_container, C
 
     /*--- Identify the boundary by string name ---*/
     string Marker_Tag = config->GetMarker_All_TagBound(val_marker);
-   
+
     if (config->GetMarker_StrongBC(Marker_Tag)==true) {
       nodes->SetSolution_Old(iPoint, Inlet_SpeciesVars[val_marker][iVertex]);
 
@@ -406,6 +406,20 @@ void CSpeciesSolver::BC_Inlet(CGeometry* geometry, CSolver** solver_container, C
   END_SU2_OMP_FOR
 }
 
+void CSpeciesSolver::BC_Isothermal_Wall_Generic(CGeometry* geometry, CSolver** solver_container, CNumerics* conv_numerics,
+                                  CNumerics* visc_numerics, CConfig* config, unsigned short val_marker,
+                                  bool cht_mode) {
+
+}
+
+void CSpeciesSolver::BC_Isothermal_Wall(CGeometry *geometry, CSolver **solver_container, CNumerics *conv_numerics,
+                                       CNumerics *visc_numerics, CConfig *config, unsigned short val_marker) {
+
+  BC_Isothermal_Wall_Generic(geometry, solver_container, conv_numerics, visc_numerics, config, val_marker);
+
+}
+
+
 void CSpeciesSolver::SetInletAtVertex(const su2double *val_inlet,
                                       unsigned short iMarker,
                                       unsigned long iVertex) {

TestCases/parallel_regression.py

@@ -302,7 +302,7 @@ def main():
     flatplate_udobj.cfg_dir   = "user_defined_functions"
     flatplate_udobj.cfg_file  = "lam_flatplate.cfg"
     flatplate_udobj.test_iter = 20
-    flatplate_udobj.test_vals = [-6.653802, -1.181430, -0.794887, 0.000611, -0.000369, 0.000736, -0.001104, 596.690000, 299.800000, 296.890000, 21.492000, 0.563990, 37.148000, 2.278700]
+    flatplate_udobj.test_vals = [-6.653278, -1.180894, -0.794611, 0.000613, -0.000370, 0.000737, -0.001107, 596.680000, 299.790000, 296.890000, 21.488000, 0.564630, 37.153000, 2.284900]
     test_list.append(flatplate_udobj)
 
     # Laminar cylinder (steady)