Merge branch 'develop' into feature_preconditioning

Author: Cristopher-Morales

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

@@ -49,7 +49,7 @@ class CPrimalGrid {
   unsigned long GlobalIndex_DomainElement;
 
   std::unique_ptr<unsigned long[]> Nodes;    /*!< \brief Global node indices of the element. */
-  std::unique_ptr<long[]> Neighbor_Elements; /*!< \brief Vector to store the elements surrounding this element. */
+  std::unique_ptr<long[]> Neighbor_Elements; /*!< \brief Vector to store the elements surronding this element. */
 
   su2double Coord_CG[3] = {0.0}; /*!< \brief Coordinates of the center-of-gravity of the element. */
   su2double Volume;              /*!< \brief Volume of the element. */

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/variables/CFlowVariable.hpp

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

SU2_PY/SU2/io/data.py

@@ -115,7 +115,7 @@ def load_data(file_name, var_names=None, file_format="infer", core_name="python_
     #: with filelock
 
     # load specified varname into dictionary
-    if var_names != None:
+    if var_names is not None:
         # check for one item name array
         if isinstance(var_names, str):
             var_names = [
@@ -365,7 +365,7 @@ def rec2dict(array_in):
         # convert array
         elif isinstance(value, numpy.ndarray):
             # check for another struct level
-            if value.dtype.names == None:
+            if value.dtype.names is None:
                 value = value.tolist()
             # telescoping
             else:

SU2_PY/SU2_CFD.py

@@ -120,7 +120,7 @@ def main():
     options.wave_equation = options.wave_equation.upper() == "TRUE"
     options.heat_equation = options.heat_equation.upper() == "TRUE"
 
-    if options.filename == None:
+    if options.filename is None:
         raise Exception("No config file provided. Use -f flag")
 
     if options.with_MPI == True:
@@ -163,7 +163,7 @@ def main():
     # Finalize the solver and exit cleanly
     SU2Driver.Finalize()
 
-    if SU2Driver != None:
+    if SU2Driver is not None:
         del SU2Driver
 
 

SU2_PY/parallel_computation.py

@@ -64,7 +64,7 @@ def main():
     options.partitions = int(options.partitions)
     options.compute = options.compute.upper() == "TRUE"
 
-    if options.filename == None:
+    if options.filename is None:
         raise Exception("No config file provided. Use -f flag")
 
     parallel_computation(options.filename, options.partitions, options.compute)

SU2_PY/parallel_computation_fsi.py

@@ -64,7 +64,7 @@ def main():
     options.partitions = int(options.partitions)
     options.compute = options.compute.upper() == "TRUE"
 
-    if options.filename == None:
+    if options.filename is None:
         raise Exception("No config file provided. Use -f flag")
 
     parallel_computation(options.filename, options.partitions, options.compute)

SU2_PY/topology_optimization.py

@@ -123,8 +123,8 @@ def obj_val(self, x):
         # clear previous output and run direct solver
         try:
             os.remove(self._objValFile)
-        except:
-            pass
+        except OSError:
+            pass  # Ignore error if file does not exist
 
         try:
             sp.call(self._objValCommand, shell=True)
@@ -136,7 +136,7 @@ def obj_val(self, x):
                     break
             # the return code of mpirun is useless, we test the value of the function
             self._assert_isfinite(val)
-        except:
+        except Exception:
             raise RuntimeError("Objective function evaluation failed")
         # end
 
@@ -150,8 +150,8 @@ def obj_der(self, x):
         # clear previous output and run direct solver
         try:
             os.remove(self._objDerFile)
-        except:
-            pass
+        except OSError:
+            pass  # Ignore error if file does not exist
         N = x.shape[0]
         y = np.ndarray((N,))
 
@@ -167,7 +167,7 @@ def obj_der(self, x):
                 self._assert_isfinite(val)
                 y[i] = val * obj_scale / var_scale
             # end
-        except:
+        except Exception:
             raise RuntimeError("Objective gradient evaluation failed")
         # end
 
@@ -186,7 +186,7 @@ def con_val(self, x):
                     val = float(lines[1].split(",")[col])
                     break
             self._assert_isfinite(val)
-        except:
+        except Exception:
             raise RuntimeError("Constraint function evaluation failed")
         # end
 
@@ -200,8 +200,8 @@ def con_der(self, x):
         # clear previous output and run solver
         try:
             os.remove(self._conDerFile)
-        except:
-            pass
+        except OSError:
+            pass  # Ignore error if file does not exist
         N = x.shape[0]
         y = np.ndarray((N,))
 
@@ -217,7 +217,7 @@ def con_der(self, x):
                 self._assert_isfinite(val)
                 y[i] = val * con_scale / var_scale
             # end
-        except:
+        except Exception:
             raise RuntimeError("Constraint function evaluation failed")
         # end
 
@@ -295,7 +295,7 @@ def finished(self):
         return self._value == self._maxi
 
     def value(self):
-        if self._func == None:
+        if self._func is None:
             return self._value
         else:
             return self._func(self._value)

TestCases/disc_adj_fsi/dyn_fsi/grad_dv.opt.ref

@@ -1,9 +1,9 @@
 INDEX	GRAD
-0	-4.570869180445024e-04
-1	-2.401466726069483e-04
-2	-9.134698219534478e-05
-3	-1.628086953114332e-05
-4	-1.741052432426715e-05
-5	-9.462489852612668e-05
-6	-2.452466275723706e-04
-7	-4.635483655821636e-04
+0	-4.570869126273463e-04
+1	-2.401466681907557e-04
+2	-9.134697929286435e-05
+3	-1.628086861828660e-05
+4	-1.741052581631559e-05
+5	-9.462490278520256e-05
+6	-2.452466349703724e-04
+7	-4.635483765277457e-04