Merge branch 'feature_custom_species_walls' of https://github.com/su2code/su2 into feature_custom_species_walls

Author: bigfooted

Common/include/geometry/CGeometry.hpp

@@ -1,4 +1,4 @@
-/*!
+/*!
  * \file CGeometry.hpp
  * \brief Headers of the main subroutines for creating the geometrical structure.
  *        The subroutines and functions are in the <i>CGeometry.cpp</i> file.
@@ -51,6 +51,7 @@ extern "C" {
 #include <climits>
 #include <memory>
 #include <unordered_map>
+#include <cstdint>
 
 #include "primal_grid/CPrimalGrid.hpp"
 #include "dual_grid/CDualGrid.hpp"

Common/include/geometry/primal_grid/CPrimalGrid.hpp

@@ -28,11 +28,11 @@
 
 #pragma once
 
+#include <cstdint>
 #include <iostream>
 #include <vector>
 #include <limits>
 #include <cstdlib>
-#include <limits>
 #include <memory>
 
 #include "../../option_structure.hpp"

Common/include/toolboxes/graph_toolbox.hpp

@@ -479,13 +479,16 @@ T createNaturalColoring(Index_t numInnerIndexes) {
  *        pattern is returned.
  * \param[in] pattern - Sparse pattern to be colored.
  * \param[in] groupSize - Size of the outer index groups, default 1.
+ * \param[in] includeOuterIdx - If the outer indices are the same type of entity as
+ *            the inner indices this must be true. For example, when both are point
+ *            indices, i.e. the pattern is the adjacency matrix.
  * \param[in] balanceColors - Try to balance number of indexes per color,
  *            tends to result in worse locality (thus false by default).
  * \param[out] indexColor - Optional, vector with colors given to the outer indices.
  * \return Coloring in the same type of the input pattern.
  */
 template <typename Color_t = unsigned char, size_t MaxColors = 255, size_t MaxMB = 128, class T>
-T colorSparsePattern(const T& pattern, size_t groupSize = 1, bool balanceColors = false,
+T colorSparsePattern(const T& pattern, size_t groupSize = 1, bool includeOuterIdx = false, bool balanceColors = false,
                      std::vector<Color_t>* indexColor = nullptr) {
   static_assert(std::is_integral<Color_t>::value, "");
   static_assert(std::numeric_limits<Color_t>::max() >= MaxColors, "");
@@ -538,6 +541,11 @@ T colorSparsePattern(const T& pattern, size_t groupSize = 1, bool balanceColors
 
       for (; it != searchOrder.end(); ++it) {
         bool free = true;
+        if (includeOuterIdx) {
+          for (Index_t k = iOuter; k < grpEnd && free; ++k) {
+            free = !innerInColor[*it][k];
+          }
+        }
         /*--- Traverse entire group as a large outer index. ---*/
         for (Index_t k = outerPtr[iOuter]; k < outerPtr[grpEnd] && free; ++k) {
           free = !innerInColor[*it][innerIdx[k] - minIdx];

Common/src/CConfig.cpp

@@ -6210,8 +6210,8 @@ void CConfig::SetMarkers(SU2_COMPONENT val_software) {
     }
   }
 
-  /*--- Idenftification fo Giles Markers ---*/
-  // This is seperate from MP and Turbomachinery Markers as all mixing plane markers are Giles,
+  /*--- Identification of Giles Markers ---*/
+  // This is separate from MP and Turbomachinery Markers as all mixing plane markers are Giles,
   // but not all Giles markers are mixing plane
   for (iMarker_CfgFile = 0; iMarker_CfgFile < nMarker_CfgFile; iMarker_CfgFile++) {
     Marker_CfgFile_Giles[iMarker_CfgFile] = NO;
@@ -10057,8 +10057,8 @@ void CConfig::SetProfilingCSV() {
 
   /*--- Allocate and initialize memory ---*/
 
-  double *l_min_red = NULL, *l_max_red = NULL, *l_tot_red = NULL, *l_avg_red = NULL;
-  int *n_calls_red = NULL;
+  double *l_min_red = nullptr, *l_max_red = nullptr, *l_tot_red = nullptr, *l_avg_red = nullptr;
+  int *n_calls_red = nullptr;
   double* l_min = new double[map_size];
   double* l_max = new double[map_size];
   double* l_tot = new double[map_size];

Common/src/geometry/CGeometry.cpp

@@ -3772,7 +3772,7 @@ const CCompressedSparsePatternUL& CGeometry::GetEdgeColoring(su2double* efficien
       bool admissibleColoring = false; /* keep track wether the last tested coloring is admissible */
 
       while (true) {
-        edgeColoring = colorSparsePattern(pattern, nextEdgeColorGroupSize, balanceColors);
+        edgeColoring = colorSparsePattern(pattern, nextEdgeColorGroupSize, false, balanceColors);
 
         /*--- If the coloring fails, reduce the color group size. ---*/
         if (edgeColoring.empty()) {
@@ -3809,12 +3809,12 @@ const CCompressedSparsePatternUL& CGeometry::GetEdgeColoring(su2double* efficien
 
       /*--- If the last tested coloring was not admissible, recompute the final coloring. ---*/
       if (!admissibleColoring) {
-        edgeColoring = colorSparsePattern(pattern, edgeColorGroupSize, balanceColors);
+        edgeColoring = colorSparsePattern(pattern, edgeColorGroupSize, false, balanceColors);
       }
     }
     /*--- No adaptivity. ---*/
     else {
-      edgeColoring = colorSparsePattern(pattern, edgeColorGroupSize, balanceColors);
+      edgeColoring = colorSparsePattern(pattern, edgeColorGroupSize, false, balanceColors);
     }
 
     /*--- If the coloring fails use the natural coloring. This is a
@@ -3867,7 +3867,7 @@ const CCompressedSparsePatternUL& CGeometry::GetElementColoring(su2double* effic
 
     /*--- Color the elements. ---*/
     constexpr bool balanceColors = true;
-    elemColoring = colorSparsePattern(pattern, elemColorGroupSize, balanceColors);
+    elemColoring = colorSparsePattern(pattern, elemColorGroupSize, false, balanceColors);
 
     /*--- Same as for the edge coloring. ---*/
     if (elemColoring.empty()) SetNaturalElementColoring();
@@ -3896,7 +3896,7 @@ void CGeometry::ColorMGLevels(unsigned short nMGLevels, const CGeometry* const*
     /*--- Color the coarse points. ---*/
     vector<tColor> color;
     const auto& adjacency = geometry[iMesh]->nodes->GetPoints();
-    if (colorSparsePattern<tColor, nColor>(adjacency, 1, false, &color).empty()) continue;
+    if (colorSparsePattern<tColor, nColor>(adjacency, 1, true, false, &color).empty()) continue;
 
     /*--- Propagate colors to fine mesh. ---*/
     for (auto step = 0u; step < iMesh; ++step) {
@@ -4042,8 +4042,8 @@ const CGeometry::CLineletInfo& CGeometry::GetLineletInfo(const CConfig* config)
     }
   }
 
-  const auto coloring =
-      colorSparsePattern<uint8_t, std::numeric_limits<uint8_t>::max()>(CCompressedSparsePatternUL(adjacency), 1, true);
+  const auto coloring = colorSparsePattern<uint8_t, std::numeric_limits<uint8_t>::max()>(
+      CCompressedSparsePatternUL(adjacency), 1, false, true);
   const auto nColors = coloring.getOuterSize();
 
   /*--- Sort linelets by color. ---*/

Common/src/geometry/CPhysicalGeometry.cpp

@@ -4642,7 +4642,7 @@ void CPhysicalGeometry::SetElement_Connectivity() {
 
           if ((elem[iElem]->GetNeighbor_Elements(iFace) == -1) && (iElem < Test_Elem) &&
               FindFace(iElem, Test_Elem, first_elem_face, second_elem_face)) {
-            /*--- Localice which faces are sharing both elements ---*/
+            /*--- Localize which faces are sharing both elements ---*/
 
             elem[iElem]->SetNeighbor_Elements(Test_Elem, first_elem_face);
 
@@ -4666,7 +4666,7 @@ void CPhysicalGeometry::SetBoundVolume() {
       CheckVol = false;
 
       for (iElem = 0; iElem < nodes->GetnElem(Point); iElem++) {
-        /*--- Look for elements surronding that point --*/
+        /*--- Look for elements surrounding that point --*/
         cont = 0;
         iElem_Domain = nodes->GetElem(Point, iElem);
         for (iNode_Domain = 0; iNode_Domain < elem[iElem_Domain]->GetnNodes(); iNode_Domain++) {
@@ -4775,10 +4775,10 @@ void CPhysicalGeometry::ComputeNSpan(CConfig* config, unsigned short val_iZone,
   nSpan_loc = 0;
   if (nDim == 2) {
     nSpanWiseSections[marker_flag - 1] = 1;
-    // TODO (turbo) make it more genral
+    // TODO (turbo) make it more general
     if (marker_flag == OUTFLOW) config->SetnSpanWiseSections(1);
 
-    /*---Initilize the vector of span-wise values that will be ordered ---*/
+    /*---Initialize the vector of span-wise values that will be ordered ---*/
     SpanWiseValue[marker_flag - 1] = new su2double[1];
     for (iSpan = 0; iSpan < 1; iSpan++) {
       SpanWiseValue[marker_flag - 1][iSpan] = 0;
@@ -4792,7 +4792,7 @@ void CPhysicalGeometry::ComputeNSpan(CConfig* config, unsigned short val_iZone,
           if (config->GetMarker_All_TurbomachineryFlag(iMarker) != marker_flag) continue;
 
           /*--- loop to find the vertex that ar both of inflow or outflow marker and on the periodic
-           * in order to caount the number of Span ---*/
+           * in order to count the number of Span ---*/
           for (jMarker = 0; jMarker < nMarker; jMarker++) {
             if (config->GetMarker_All_KindBC(jMarker) != PERIODIC_BOUNDARY) continue;
 

SU2_CFD/include/numerics_simd/flow/convection/common.hpp

@@ -51,9 +51,7 @@ FORCEINLINE Double umusclProjection(Double grad_proj,
   /*--- To maintain proper scaling for edge limiters, the result of ---*/
   /*--- this function is 0.5 * dV_ij^kap.                           ---*/
   /*-------------------------------------------------------------------*/
-  const Double cent = 0.5 * delta;
-  const Double upw = grad_proj - cent;
-  return (1.0 - kappa) * upw + (1.0 + kappa) * cent;
+  return (1.0 - kappa) * grad_proj + kappa * delta;
 }
 
 /*!
@@ -216,7 +214,7 @@ FORCEINLINE CPair<ReconVarType> reconstructPrimitives(Int iEdge, Int iPoint, Int
   }
 
   if (muscl) {
-    /*--- Recompute density and enthalpy instead of reconstructing. ---*/
+    /*--- Reconstruct density and enthalpy without using their gradients. ---*/
     constexpr auto nVarGrad = ReconVarType::nVar - 2;
     switch (limiterType) {
     case LIMITER::NONE:
@@ -229,12 +227,23 @@ FORCEINLINE CPair<ReconVarType> reconstructPrimitives(Int iEdge, Int iPoint, Int
       musclPointLimited<nVarGrad>(iPoint, jPoint, vector_ij, limiters, gradients, V, kappa, relax);
       break;
     }
+    /*--- Recompute density using the reconstructed pressure and temperature. ---*/
     V.i.density() = V.i.pressure() / (gasConst * V.i.temperature());
     V.j.density() = V.j.pressure() / (gasConst * V.j.temperature());
 
+    /*--- Reconstruct enthalpy using dH/dT = Cp and dH/dv = v. Recomputing enthalpy would cause
+     * stability issues because we use rho E = rho H - P, which loses its relation to temperature
+     * if both rho and H are recomputed. This only seems to be an issue for wall-function meshes.
+     * NOTE: This "one-sided" reconstruction does not lose much of the U-MUSCL benefit, because
+     * the static enthalpy is linear, and for the KE term we do the equivalent of using the
+     * average dH/dv, which is exact for quadratic functions. ---*/
     const su2double cp = gasConst * gamma / (gamma - 1);
-    V.i.enthalpy() = cp * V.i.temperature() + 0.5 * squaredNorm<nDim>(V.i.velocity());
-    V.j.enthalpy() = cp * V.j.temperature() + 0.5 * squaredNorm<nDim>(V.j.velocity());
+    V.i.enthalpy() += cp * (V.i.temperature() - V1st.i.temperature());
+    V.j.enthalpy() += cp * (V.j.temperature() - V1st.j.temperature());
+    for (size_t iDim = 0; iDim < nDim; ++iDim) {
+      V.i.enthalpy() += 0.5 * (pow(V.i.velocity(iDim), 2) - pow(V1st.i.velocity(iDim), 2));
+      V.j.enthalpy() += 0.5 * (pow(V.j.velocity(iDim), 2) - pow(V1st.j.velocity(iDim), 2));
+    }
 
     /*--- Detect a non-physical reconstruction based on negative pressure or density. ---*/
     const Double neg_p_or_rho = fmax(fmin(V.i.pressure(), V.j.pressure()) < 0.0,

SU2_CFD/include/solvers/CGradientSmoothingSolver.hpp

@@ -275,7 +275,7 @@ class CGradientSmoothingSolver final : public CFEASolverBase {
   }
 
   /*!
-   * \brief Set the current working dimension, if the seperate dimension option is set.
+   * \brief Set the current working dimension, if the separate dimension option is set.
    * \param[in] iDim - the dimension we are currently working in.
    */
   inline void SetCurrentDim(unsigned int iDim) {

SU2_CFD/include/variables/CFlowVariable.hpp

@@ -26,6 +26,9 @@
 
 #pragma once
 
+#include <algorithm>
+#include <cstdint>
+
 #include "CVariable.hpp"
 
 /*!

SU2_CFD/src/solvers/CGradientSmoothingSolver.cpp

@@ -93,7 +93,7 @@ CGradientSmoothingSolver::CGradientSmoothingSolver(CGeometry *geometry, CConfig
     element_container[GRAD_TERM][EL_PYRAM2] = new CPYRAM6();
   }
 
-  /*--- For operations on surfaces we initalize the structures for nDim-1 ---*/
+  /*--- For operations on surfaces we initialize the structures for nDim-1 ---*/
   if (config->GetSmoothOnSurface()) {
     if (nDim == 2) {
       element_container[GRAD_TERM][EL_LINE] = new CLINE();