(format) Switch to format 14 in code

blast/src/DBoundaryLayer.h

@@ -7,153 +7,162 @@
 
 #include <algorithm>
 
-namespace blast {
+namespace blast
+{
 
 /**
  * @brief Boundary layer region upper/lower side or wake.
  * @author Paul Dechamps
  */
-class BLAST_API BoundaryLayer : public fwk::wSharedObject {
+class BLAST_API BoundaryLayer : public fwk::wSharedObject
+{
 
 private:
-  unsigned int const nVar =
-      5; ///< Number of variables of the partial differential problem.
-  double nCrit = 9.0; ///< Critical amplification factor.
+    unsigned int const nVar =
+        5;              ///< Number of variables of the partial differential problem.
+    double nCrit = 9.0; ///< Critical amplification factor.
 
 public:
-  Closures *closSolver; ///< Closure relations class.
-  std::string name;     ///< Name of the region.
-
-  // Mesh
-  size_t nNodes;           ///< Number of points in the domain.
-  size_t nElms;            ///< Number of cells in the domain.
-  std::vector<double> x;   ///< x coordinate in the global frame of reference.
-  std::vector<double> y;   ///< y coordinate in the global frame of reference.
-  std::vector<double> z;   ///< z coordinate in the global frame of reference.
-  std::vector<double> xoc; ///< x/c coordinate in the global frame of reference.
-  std::vector<double> loc; ///< xi coordinate in the local frame of reference.
-  std::vector<int> rows;   ///< Reference numbers of the nodes.
-  std::vector<int> no;     ///< Reference numbers of the nodes.
-  std::vector<double> dx;  ///< Cell size.
-  std::vector<double> alpha; ///< Angle of the cell wrt the x axis of the global
-                             ///< frame of reference.
-  double chord;              ///< Chord of the section.
-  double xMin; ///< Minimum x coordinate of the mesh (for local coordinates
-               ///< computation).
-
-  // Boundary layer variables.
-  std::vector<double> u; ///< Solution vector (theta, H, N, ue, Ct).
-  std::vector<double>
-      Ret; ///< Reynolds number based on the momentum thickness (theta).
-  std::vector<double> cd; ///< Local dissipation coefficient.
-  std::vector<double> cf; ///< Local friction coefficient.
-  std::vector<double>
-      Hstar; ///< Kinetic energy shape parameter (thetaStar/theta).
-  std::vector<double>
-      Hstar2;             ///< Density shape parameter (deltaStarStar/theta).
-  std::vector<double> Hk; ///< Kinematic shape parameter (int(1-u/ue d_eta)).
-  std::vector<double>
-      ctEq; ///< Equilibrium shear stress coefficient (turbulent BL).
-  std::vector<double> us;    ///< Equivalent normalized wall slip velocity.
-  std::vector<double> delta; ///< Boundary layer thickness.
-  std::vector<double>
-      deltaStar; ///< Displacement thickness (int(1-rho*u/rhoe*ue d_eta)).
-
-  std::vector<double> Me;    ///< Mach number.
-  std::vector<double> vt;    ///< Velocity.
-  std::vector<double> rhoe;  ///< Density.
-  std::vector<double> xxExt; ///< xi coordinate in the local frame of reference,
-                             ///< at the edge of the boundary layer.
-  std::vector<double> deltaStarExt; ///< Displacement thickness.
-  std::vector<double> vtExt;        ///< Previous velocity.
-
-  std::vector<double> blowingVelocity; ///< Blowing velocity.
-
-  // Transition related variables.
-  std::vector<int> regime; ///< Laminar (0) or turbulent (1) regime.
-  size_t transitionNode;   ///< Node id of the transition location.
-  double xtrF; ///< Forced transition location in the local frame of reference
-               ///< (x/c).
-  double xtr;  ///< Transition location in the local frame of reference (x/c).
-
-  BoundaryLayer(double _xtrF = -1.0, std::string _name = "None");
-  ~BoundaryLayer();
-
-  void reset();
-
-  // Getters
-  std::string getName() const { return name; };
-  size_t getnNodes() const { return nNodes; };
-  size_t getnElms() const { return nElms; };
-
-  std::vector<double> getDeltaStar() const { return deltaStar; };
-  std::vector<double> getUe() const {
-    std::vector<double> ue(nNodes, 0.);
-    for (size_t i = 0; i < nNodes; ++i)
-      ue[i] = u[i * nVar + 3];
-    return ue;
-  };
-  std::vector<double> getBlowing() const { return blowingVelocity; };
-  double getMaxMach() const { return *std::max_element(Me.begin(), Me.end()); };
-
-  // Setters
-  void setMesh(std::vector<double> const _x, std::vector<double> const y,
-               std::vector<double> const z, double const _chord,
-               double const _xMin,
-               std::vector<int> const _rows = std::vector<int>(0),
-               std::vector<int> const _no = std::vector<int>(0));
-  void computeLocalCoordinates();
-  void setMe(std::vector<double> const _Me) {
-    if (_Me.size() != nNodes)
-      throw std::runtime_error("Mach number vector is not consistent.");
-    Me = _Me;
-  };
-  void setVt(std::vector<double> const _vt) {
-    if (_vt.size() != nNodes)
-      throw std::runtime_error("Velocity vector is not consistent.");
-    vt = _vt;
-  };
-  void setRhoe(std::vector<double> const _rhoe) {
-    if (_rhoe.size() != nNodes)
-      throw std::runtime_error("Density vector is not consistent.");
-    rhoe = _rhoe;
-  };
-  void setxxExt(std::vector<double> const _xxExt) {
-    if (_xxExt.size() != nNodes)
-      throw std::runtime_error("External mesh vector is not consistent.");
-    xxExt = _xxExt;
-  };
-  void setDeltaStarExt(std::vector<double> const _deltaStarExt) {
-    if (_deltaStarExt.size() != nNodes)
-      throw std::runtime_error(
-          "Displacement thickness vector is not consistent.");
-    deltaStarExt = _deltaStarExt;
-  };
-  void setVtExt(std::vector<double> const _vtExt) {
-    if (_vtExt.size() != nNodes)
-      throw std::runtime_error("Velocity vector is not consistent.");
-    vtExt = _vtExt;
-  };
-  void setxtrF(double const _xtrF) { xtrF = _xtrF; };
-
-  // Boundary conditions.
-  void setStagnationBC(double Re);
-  void setWakeBC(double Re, std::vector<double> upperConditions,
-                 std::vector<double> lowerConditions);
-
-  // Getters.
-  size_t getnVar() const { return nVar; };
-  double getnCrit() const { return nCrit; };
-  void setnCrit(double const _nCrit) { nCrit = _nCrit; };
-
-  // Others
-  void printSolution(size_t iPoint) const;
-  void computeBlowingVelocity();
-  std::vector<std::vector<double>> evalGradwrtRho();
-  std::vector<std::vector<double>> evalGradwrtVt();
-  std::vector<std::vector<double>> evalGradwrtDeltaStar();
-  std::vector<std::vector<double>> evalGradwrtX();
-  std::vector<std::vector<double>> evalGradwrtY();
+    Closures *closSolver; ///< Closure relations class.
+    std::string name;     ///< Name of the region.
+
+    // Mesh
+    size_t nNodes;             ///< Number of points in the domain.
+    size_t nElms;              ///< Number of cells in the domain.
+    std::vector<double> x;     ///< x coordinate in the global frame of reference.
+    std::vector<double> y;     ///< y coordinate in the global frame of reference.
+    std::vector<double> z;     ///< z coordinate in the global frame of reference.
+    std::vector<double> xoc;   ///< x/c coordinate in the global frame of reference.
+    std::vector<double> loc;   ///< xi coordinate in the local frame of reference.
+    std::vector<int> rows;     ///< Reference numbers of the nodes.
+    std::vector<int> no;       ///< Reference numbers of the nodes.
+    std::vector<double> dx;    ///< Cell size.
+    std::vector<double> alpha; ///< Angle of the cell wrt the x axis of the global
+                               ///< frame of reference.
+    double chord;              ///< Chord of the section.
+    double xMin;               ///< Minimum x coordinate of the mesh (for local coordinates
+                               ///< computation).
+
+    // Boundary layer variables.
+    std::vector<double> u; ///< Solution vector (theta, H, N, ue, Ct).
+    std::vector<double>
+        Ret;                ///< Reynolds number based on the momentum thickness (theta).
+    std::vector<double> cd; ///< Local dissipation coefficient.
+    std::vector<double> cf; ///< Local friction coefficient.
+    std::vector<double>
+        Hstar; ///< Kinetic energy shape parameter (thetaStar/theta).
+    std::vector<double>
+        Hstar2;             ///< Density shape parameter (deltaStarStar/theta).
+    std::vector<double> Hk; ///< Kinematic shape parameter (int(1-u/ue d_eta)).
+    std::vector<double>
+        ctEq;                  ///< Equilibrium shear stress coefficient (turbulent BL).
+    std::vector<double> us;    ///< Equivalent normalized wall slip velocity.
+    std::vector<double> delta; ///< Boundary layer thickness.
+    std::vector<double>
+        deltaStar; ///< Displacement thickness (int(1-rho*u/rhoe*ue d_eta)).
+
+    std::vector<double> Me;           ///< Mach number.
+    std::vector<double> vt;           ///< Velocity.
+    std::vector<double> rhoe;         ///< Density.
+    std::vector<double> xxExt;        ///< xi coordinate in the local frame of reference,
+                                      ///< at the edge of the boundary layer.
+    std::vector<double> deltaStarExt; ///< Displacement thickness.
+    std::vector<double> vtExt;        ///< Previous velocity.
+
+    std::vector<double> blowingVelocity; ///< Blowing velocity.
+
+    // Transition related variables.
+    std::vector<int> regime; ///< Laminar (0) or turbulent (1) regime.
+    size_t transitionNode;   ///< Node id of the transition location.
+    double xtrF;             ///< Forced transition location in the local frame of reference
+                             ///< (x/c).
+    double xtr;              ///< Transition location in the local frame of reference (x/c).
+
+    BoundaryLayer(double _xtrF = -1.0, std::string _name = "None");
+    ~BoundaryLayer();
+
+    void reset();
+
+    // Getters
+    std::string getName() const { return name; };
+    size_t getnNodes() const { return nNodes; };
+    size_t getnElms() const { return nElms; };
+
+    std::vector<double> getDeltaStar() const { return deltaStar; };
+    std::vector<double> getUe() const
+    {
+        std::vector<double> ue(nNodes, 0.);
+        for (size_t i = 0; i < nNodes; ++i)
+            ue[i] = u[i * nVar + 3];
+        return ue;
+    };
+    std::vector<double> getBlowing() const { return blowingVelocity; };
+    double getMaxMach() const { return *std::max_element(Me.begin(), Me.end()); };
+
+    // Setters
+    void setMesh(std::vector<double> const _x, std::vector<double> const y,
+                 std::vector<double> const z, double const _chord,
+                 double const _xMin,
+                 std::vector<int> const _rows = std::vector<int>(0),
+                 std::vector<int> const _no = std::vector<int>(0));
+    void computeLocalCoordinates();
+    void setMe(std::vector<double> const _Me)
+    {
+        if (_Me.size() != nNodes)
+            throw std::runtime_error("Mach number vector is not consistent.");
+        Me = _Me;
+    };
+    void setVt(std::vector<double> const _vt)
+    {
+        if (_vt.size() != nNodes)
+            throw std::runtime_error("Velocity vector is not consistent.");
+        vt = _vt;
+    };
+    void setRhoe(std::vector<double> const _rhoe)
+    {
+        if (_rhoe.size() != nNodes)
+            throw std::runtime_error("Density vector is not consistent.");
+        rhoe = _rhoe;
+    };
+    void setxxExt(std::vector<double> const _xxExt)
+    {
+        if (_xxExt.size() != nNodes)
+            throw std::runtime_error("External mesh vector is not consistent.");
+        xxExt = _xxExt;
+    };
+    void setDeltaStarExt(std::vector<double> const _deltaStarExt)
+    {
+        if (_deltaStarExt.size() != nNodes)
+            throw std::runtime_error(
+                "Displacement thickness vector is not consistent.");
+        deltaStarExt = _deltaStarExt;
+    };
+    void setVtExt(std::vector<double> const _vtExt)
+    {
+        if (_vtExt.size() != nNodes)
+            throw std::runtime_error("Velocity vector is not consistent.");
+        vtExt = _vtExt;
+    };
+    void setxtrF(double const _xtrF) { xtrF = _xtrF; };
+
+    // Boundary conditions.
+    void setStagnationBC(double Re);
+    void setWakeBC(double Re, std::vector<double> upperConditions,
+                   std::vector<double> lowerConditions);
+
+    // Getters.
+    size_t getnVar() const { return nVar; };
+    double getnCrit() const { return nCrit; };
+    void setnCrit(double const _nCrit) { nCrit = _nCrit; };
+
+    // Others
+    void printSolution(size_t iPoint) const;
+    void computeBlowingVelocity();
+    std::vector<std::vector<double>> evalGradwrtRho();
+    std::vector<std::vector<double>> evalGradwrtVt();
+    std::vector<std::vector<double>> evalGradwrtDeltaStar();
+    std::vector<std::vector<double>> evalGradwrtX();
+    std::vector<std::vector<double>> evalGradwrtY();
 };
 } // namespace blast
 #endif // DBOUNDARYLAYER_H

blast/src/DClosures.cpp

@@ -16,66 +16,72 @@ using namespace blast;
  */
 void Closures::laminarClosures(std::vector<double> &closureVars, double theta,
                                double H, double Ret, const double Me,
-                               const std::string name) {
-  H = std::max(H, 1.00005);
-
-  // Kinematic shape factor H*.
-  double Hk = (H - 0.29 * Me * Me) / (1 + 0.113 * Me * Me);
-  if (name == "wake")
-    Hk = std::max(Hk, 1.00005);
-  else
-    Hk = std::max(Hk, 1.05000);
-
-  // Density shape parameter.
-  double Hstar2 = (0.064 / (Hk - 0.8) + 0.251) * Me * Me;
-
-  // Boundary layer thickness.
-  double delta = std::min((3.15 + H + (1.72 / (Hk - 1))) * theta, 12 * theta);
-
-  double Hstar = 0.;
-  double Hks = Hk - 4.35;
-  if (Hk <= 4.35)
-    Hstar = 0.0111 * Hks * Hks / (Hk + 1) -
-            0.0278 * Hks * Hks * Hks / (Hk + 1) + 1.528 -
-            0.0002 * (Hks * Hk) * (Hks * Hk);
-  else
-    Hstar = 1.528 + 0.015 * Hks * Hks / Hk;
-
-  // Whitfield's minor additional compressibility correction.
-  Hstar = (Hstar + 0.028 * Me * Me) / (1 + 0.014 * Me * Me);
-
-  // Friction coefficient.
-  double tmp = 0.;
-  double cf = 0.;
-  if (Hk < 5.5) {
-    tmp = (5.5 - Hk) * (5.5 - Hk) * (5.5 - Hk) / (Hk + 1.0);
-    cf = (-0.07 + 0.0727 * tmp) / Ret;
-  } else if (Hk >= 5.5) {
-    tmp = 1.0 - 1.0 / (Hk - 4.5);
-    cf = (-0.07 + 0.015 * tmp * tmp) / Ret;
-  }
-
-  // Dissipation coefficient.
-  double Cd = 0.;
-  if (Hk < 4)
-    Cd = (0.00205 * std::pow(4.0 - Hk, 5.5) + 0.207) * (Hstar / (2 * Ret));
-  else {
-    double HkCd = (Hk - 4) * (Hk - 4);
-    double denCd = 1 + 0.02 * HkCd;
-    Cd = (-0.0016 * HkCd / denCd + 0.207) * (Hstar / (2 * Ret));
-  }
-
-  // Wake relations.
-  if (name == "wake") {
-    Cd = 2 * (1.10 * (1.0 - 1.0 / Hk) * (1.0 - 1.0 / Hk) / Hk) *
-         (Hstar / (2 * Ret));
-    cf = 0.;
-  }
-
-  double us = 0.;
-  double ctEq = 0.;
-
-  closureVars = {Hstar, Hstar2, Hk, cf, Cd, delta, ctEq, us};
+                               const std::string name)
+{
+    H = std::max(H, 1.00005);
+
+    // Kinematic shape factor H*.
+    double Hk = (H - 0.29 * Me * Me) / (1 + 0.113 * Me * Me);
+    if (name == "wake")
+        Hk = std::max(Hk, 1.00005);
+    else
+        Hk = std::max(Hk, 1.05000);
+
+    // Density shape parameter.
+    double Hstar2 = (0.064 / (Hk - 0.8) + 0.251) * Me * Me;
+
+    // Boundary layer thickness.
+    double delta = std::min((3.15 + H + (1.72 / (Hk - 1))) * theta, 12 * theta);
+
+    double Hstar = 0.;
+    double Hks = Hk - 4.35;
+    if (Hk <= 4.35)
+        Hstar = 0.0111 * Hks * Hks / (Hk + 1) -
+                0.0278 * Hks * Hks * Hks / (Hk + 1) + 1.528 -
+                0.0002 * (Hks * Hk) * (Hks * Hk);
+    else
+        Hstar = 1.528 + 0.015 * Hks * Hks / Hk;
+
+    // Whitfield's minor additional compressibility correction.
+    Hstar = (Hstar + 0.028 * Me * Me) / (1 + 0.014 * Me * Me);
+
+    // Friction coefficient.
+    double tmp = 0.;
+    double cf = 0.;
+    if (Hk < 5.5)
+    {
+        tmp = (5.5 - Hk) * (5.5 - Hk) * (5.5 - Hk) / (Hk + 1.0);
+        cf = (-0.07 + 0.0727 * tmp) / Ret;
+    }
+    else if (Hk >= 5.5)
+    {
+        tmp = 1.0 - 1.0 / (Hk - 4.5);
+        cf = (-0.07 + 0.015 * tmp * tmp) / Ret;
+    }
+
+    // Dissipation coefficient.
+    double Cd = 0.;
+    if (Hk < 4)
+        Cd = (0.00205 * std::pow(4.0 - Hk, 5.5) + 0.207) * (Hstar / (2 * Ret));
+    else
+    {
+        double HkCd = (Hk - 4) * (Hk - 4);
+        double denCd = 1 + 0.02 * HkCd;
+        Cd = (-0.0016 * HkCd / denCd + 0.207) * (Hstar / (2 * Ret));
+    }
+
+    // Wake relations.
+    if (name == "wake")
+    {
+        Cd = 2 * (1.10 * (1.0 - 1.0 / Hk) * (1.0 - 1.0 / Hk) / Hk) *
+             (Hstar / (2 * Ret));
+        cf = 0.;
+    }
+
+    double us = 0.;
+    double ctEq = 0.;
+
+    closureVars = {Hstar, Hstar2, Hk, cf, Cd, delta, ctEq, us};
 }
 
 /**
@@ -91,113 +97,117 @@ void Closures::laminarClosures(std::vector<double> &closureVars, double theta,
  */
 void Closures::turbulentClosures(std::vector<double> &closureVars, double theta,
                                  double H, double Ct, double Ret,
-                                 const double Me, const std::string name) {
-  H = std::max(H, 1.00005);
-  Ct = std::min(Ct, 0.3);
-  // Ct = std::max(std::min(Ct, 0.30), 0.0000001);
-
-  double Hk = (H - 0.29 * Me * Me) / (1. + 0.113 * Me * Me);
-  if (name == "wake")
-    Hk = std::max(Hk, 1.00005);
-  else
-    Hk = std::max(Hk, 1.05000);
-
-  double Hstar2 = ((0.064 / (Hk - 0.8)) + 0.251) * Me * Me;
-  double gamma = 1.4 - 1.;
-  double Fc = std::sqrt(1. + 0.5 * gamma * Me * Me);
-
-  double H0 = 0.;
-  if (Ret <= 400.)
-    H0 = 4.0;
-  else
-    H0 = 3. + (400. / Ret);
-  if (Ret <= 200.)
-    Ret = 200.;
-
-  double Hstar = 0.;
-  if (Hk <= H0)
-    Hstar =
-        ((0.5 - 4. / Ret) * ((H0 - Hk) / (H0 - 1.)) * ((H0 - Hk) / (H0 - 1.))) *
-            (1.5 / (Hk + 0.5)) +
-        1.5 + 4. / Ret;
-  else
-    Hstar = (Hk - H0) * (Hk - H0) *
-                (0.007 * std::log(Ret) /
-                     ((Hk - H0 + 4. / std::log(Ret)) *
-                      (Hk - H0 + 4. / std::log(Ret))) +
-                 0.015 / Hk) +
+                                 const double Me, const std::string name)
+{
+    H = std::max(H, 1.00005);
+    Ct = std::min(Ct, 0.3);
+    // Ct = std::max(std::min(Ct, 0.30), 0.0000001);
+
+    double Hk = (H - 0.29 * Me * Me) / (1. + 0.113 * Me * Me);
+    if (name == "wake")
+        Hk = std::max(Hk, 1.00005);
+    else
+        Hk = std::max(Hk, 1.05000);
+
+    double Hstar2 = ((0.064 / (Hk - 0.8)) + 0.251) * Me * Me;
+    double gamma = 1.4 - 1.;
+    double Fc = std::sqrt(1. + 0.5 * gamma * Me * Me);
+
+    double H0 = 0.;
+    if (Ret <= 400.)
+        H0 = 4.0;
+    else
+        H0 = 3. + (400. / Ret);
+    if (Ret <= 200.)
+        Ret = 200.;
+
+    double Hstar = 0.;
+    if (Hk <= H0)
+        Hstar =
+            ((0.5 - 4. / Ret) * ((H0 - Hk) / (H0 - 1.)) * ((H0 - Hk) / (H0 - 1.))) *
+                (1.5 / (Hk + 0.5)) +
             1.5 + 4. / Ret;
+    else
+        Hstar = (Hk - H0) * (Hk - H0) *
+                    (0.007 * std::log(Ret) /
+                         ((Hk - H0 + 4. / std::log(Ret)) *
+                          (Hk - H0 + 4. / std::log(Ret))) +
+                     0.015 / Hk) +
+                1.5 + 4. / Ret;
+
+    // Whitfield's minor additional compressibility correction.
+    Hstar = (Hstar + 0.028 * Me * Me) / (1. + 0.014 * Me * Me);
+
+    double logRt = std::log(Ret / Fc);
+    logRt = std::max(logRt, 3.0);
+    double arg = std::max(-1.33 * Hk, -20.0);
+
+    // Equivalent normalized wall slip velocity.
+    double us = (Hstar / 2.) * (1. - 4. * (Hk - 1.) / (3. * H));
+
+    // Boundary layer thickness.
+    double delta = std::min((3.15 + H + (1.72 / (Hk - 1.))) * theta, 12. * theta);
+
+    double Ctcon = 0.5 / (6.7 * 6.7 * 0.75);
+    double Hkc = 0.;
+    double Cdw = 0.;
+    double Cdd = 0.;
+    double Cdl = 0.;
+
+    double Hmin = 0.;
+    double Fl = 0.;
+    double Dfac = 0.;
+
+    double cf = 0.;
+    double Cd = 0.;
+
+    double n = 1.0;
+
+    double ctEq = 0.;
+
+    if (name == "wake")
+    {
+        if (us > 0.99995)
+            us = 0.99995;
+        n = 2.0;  // two wake halves
+        cf = 0.0; // no friction inside the wake
+        Hkc = Hk - 1.;
+        Cdw = 0.0;                    // Wall contribution.
+        Cdd = (0.995 - us) * Ct * Ct; // Turbulent outer layer contribution.
+        Cdl = 0.15 * (0.995 - us) * (0.995 - us) /
+              Ret; // Laminar stress contribution to outer layer.
+        ctEq = std::sqrt(4. * Hstar * Ctcon * ((Hk - 1.) * Hkc * Hkc) /
+                         ((1. - us) * (Hk * Hk) * H)); // Here it is ctEq^0.5.
+    }
+    else
+    {
+        if (us > 0.95)
+            us = 0.98;
+        n = 1.0;
+        cf = (1 / Fc) *
+             (0.3 * std::exp(arg) * std::pow(logRt / 2.3026, (-1.74 - 0.31 * Hk)) +
+              0.00011 * (std::tanh(4. - (Hk / 0.875)) - 1.));
+        Hkc = std::max(Hk - 1. - 18. / Ret, 0.01);
+        // Dissipation coefficient.
+        Hmin = 1. + 2.1 / std::log(Ret);
+        Fl = (Hk - 1.) / (Hmin - 1);
+        Dfac = 0.5 + 0.5 * std::tanh(Fl);
+        Cdw = 0.5 * (cf * us) * Dfac;
+        Cdd = (0.995 - us) * Ct * Ct;
+        Cdl = 0.15 * (0.995 - us) * (0.995 - us) / Ret;
+        ctEq = std::sqrt(Hstar * Ctcon * ((Hk - 1.) * Hkc * Hkc) /
+                         ((1. - us) * (Hk * Hk) * H)); // Here it is ctEq^0.5
+                                                       // ctEq = std::sqrt(Hstar * 0.015/(1-us) * (Hk-1)*(Hk-1)*(Hk-1)/(Hk*Hk*H));
+                                                       // // Drela 1987
+    }
+    if (n * Hstar * Ctcon * ((Hk - 1.) * Hkc * Hkc) /
+            ((1. - us) * (Hk * Hk) * H) <
+        0.)
+        std::cout << "Negative sqrt encountered " << std::endl;
 
-  // Whitfield's minor additional compressibility correction.
-  Hstar = (Hstar + 0.028 * Me * Me) / (1. + 0.014 * Me * Me);
-
-  double logRt = std::log(Ret / Fc);
-  logRt = std::max(logRt, 3.0);
-  double arg = std::max(-1.33 * Hk, -20.0);
-
-  // Equivalent normalized wall slip velocity.
-  double us = (Hstar / 2.) * (1. - 4. * (Hk - 1.) / (3. * H));
-
-  // Boundary layer thickness.
-  double delta = std::min((3.15 + H + (1.72 / (Hk - 1.))) * theta, 12. * theta);
-
-  double Ctcon = 0.5 / (6.7 * 6.7 * 0.75);
-  double Hkc = 0.;
-  double Cdw = 0.;
-  double Cdd = 0.;
-  double Cdl = 0.;
-
-  double Hmin = 0.;
-  double Fl = 0.;
-  double Dfac = 0.;
-
-  double cf = 0.;
-  double Cd = 0.;
-
-  double n = 1.0;
-
-  double ctEq = 0.;
-
-  if (name == "wake") {
-    if (us > 0.99995)
-      us = 0.99995;
-    n = 2.0;  // two wake halves
-    cf = 0.0; // no friction inside the wake
-    Hkc = Hk - 1.;
-    Cdw = 0.0;                    // Wall contribution.
-    Cdd = (0.995 - us) * Ct * Ct; // Turbulent outer layer contribution.
-    Cdl = 0.15 * (0.995 - us) * (0.995 - us) /
-          Ret; // Laminar stress contribution to outer layer.
-    ctEq = std::sqrt(4. * Hstar * Ctcon * ((Hk - 1.) * Hkc * Hkc) /
-                     ((1. - us) * (Hk * Hk) * H)); // Here it is ctEq^0.5.
-  } else {
-    if (us > 0.95)
-      us = 0.98;
-    n = 1.0;
-    cf = (1 / Fc) *
-         (0.3 * std::exp(arg) * std::pow(logRt / 2.3026, (-1.74 - 0.31 * Hk)) +
-          0.00011 * (std::tanh(4. - (Hk / 0.875)) - 1.));
-    Hkc = std::max(Hk - 1. - 18. / Ret, 0.01);
     // Dissipation coefficient.
-    Hmin = 1. + 2.1 / std::log(Ret);
-    Fl = (Hk - 1.) / (Hmin - 1);
-    Dfac = 0.5 + 0.5 * std::tanh(Fl);
-    Cdw = 0.5 * (cf * us) * Dfac;
-    Cdd = (0.995 - us) * Ct * Ct;
-    Cdl = 0.15 * (0.995 - us) * (0.995 - us) / Ret;
-    ctEq = std::sqrt(Hstar * Ctcon * ((Hk - 1.) * Hkc * Hkc) /
-                     ((1. - us) * (Hk * Hk) * H)); // Here it is ctEq^0.5
-    // ctEq = std::sqrt(Hstar * 0.015/(1-us) * (Hk-1)*(Hk-1)*(Hk-1)/(Hk*Hk*H));
-    // // Drela 1987
-  }
-  if (n * Hstar * Ctcon * ((Hk - 1.) * Hkc * Hkc) /
-          ((1. - us) * (Hk * Hk) * H) <
-      0.)
-    std::cout << "Negative sqrt encountered " << std::endl;
-
-  // Dissipation coefficient.
-  Cd = n * (Cdw + Cdd + Cdl);
-  closureVars = {Hstar, Hstar2, Hk, cf, Cd, delta, ctEq, us};
+    Cd = n * (Cdw + Cdd + Cdl);
+    closureVars = {Hstar, Hstar2, Hk, cf, Cd, delta, ctEq, us};
 }
 
 /**
@@ -213,65 +223,69 @@ void Closures::turbulentClosures(std::vector<double> &closureVars, double theta,
  */
 void Closures::turbulentClosures(double &closureVars, double theta, double H,
                                  double Ct, double Ret, const double Me,
-                                 const std::string name) {
-  H = std::max(H, 1.00005);
-  Ct = std::min(Ct, 0.3);
-
-  double Hk = (H - 0.29 * Me * Me) / (1 + 0.113 * Me * Me);
-  if (name == "wake")
-    Hk = std::max(Hk, 1.00005);
-  else
-    Hk = std::max(Hk, 1.05000);
-
-  double H0 = 0.;
-  if (Ret <= 400.)
-    H0 = 4.0;
-  else
-    H0 = 3. + (400. / Ret);
-  if (Ret <= 200.)
-    Ret = 200.;
-
-  double Hstar = 0.;
-  if (Hk <= H0)
-    Hstar =
-        ((0.5 - 4. / Ret) * ((H0 - Hk) / (H0 - 1.)) * ((H0 - Hk) / (H0 - 1.))) *
-            (1.5 / (Hk + 0.5)) +
-        1.5 + 4. / Ret;
-  else
-    Hstar = (Hk - H0) * (Hk - H0) *
-                (0.007 * std::log(Ret) /
-                     ((Hk - H0 + 4. / std::log(Ret)) *
-                      (Hk - H0 + 4. / std::log(Ret))) +
-                 0.015 / Hk) +
+                                 const std::string name)
+{
+    H = std::max(H, 1.00005);
+    Ct = std::min(Ct, 0.3);
+
+    double Hk = (H - 0.29 * Me * Me) / (1 + 0.113 * Me * Me);
+    if (name == "wake")
+        Hk = std::max(Hk, 1.00005);
+    else
+        Hk = std::max(Hk, 1.05000);
+
+    double H0 = 0.;
+    if (Ret <= 400.)
+        H0 = 4.0;
+    else
+        H0 = 3. + (400. / Ret);
+    if (Ret <= 200.)
+        Ret = 200.;
+
+    double Hstar = 0.;
+    if (Hk <= H0)
+        Hstar =
+            ((0.5 - 4. / Ret) * ((H0 - Hk) / (H0 - 1.)) * ((H0 - Hk) / (H0 - 1.))) *
+                (1.5 / (Hk + 0.5)) +
             1.5 + 4. / Ret;
-
-  // Whitfield's minor additional compressibility correction.
-  Hstar = (Hstar + 0.028 * Me * Me) / (1. + 0.014 * Me * Me);
-
-  // Equivalent normalized wall slip velocity.
-  double us = (Hstar / 2.) * (1 - 4 * (Hk - 1.) / (3. * H));
-
-  double Ctcon = 0.5 / (6.7 * 6.7 * 0.75);
-
-  double Hkc = 0.;
-  double ctEq = 0.;
-
-  if (name == "wake") {
-    if (us > 0.99995)
-      us = 0.99995;
-    Hkc = Hk - 1.;
-    ctEq = std::sqrt(4. * Hstar * Ctcon * ((Hk - 1.) * Hkc * Hkc) /
-                     ((1. - us) * (Hk * Hk) * H)); // Here it is ctEq^0.5
-  } else {
-    if (us > 0.95)
-      us = 0.98;
-    Hkc = std::max(Hk - 1. - 18. / Ret, 0.01);
-    ctEq = std::sqrt(Hstar * Ctcon * ((Hk - 1.) * Hkc * Hkc) /
-                     ((1. - us) * (Hk * Hk) * H)); // Here it is ctEq^0.5
-  }
-  if (Hstar * Ctcon * ((Hk - 1.) * Hkc * Hkc) / ((1. - us) * (Hk * Hk) * H) <
-      0.)
-    std::cout << "Negative sqrt encountered " << std::endl;
-
-  closureVars = ctEq;
+    else
+        Hstar = (Hk - H0) * (Hk - H0) *
+                    (0.007 * std::log(Ret) /
+                         ((Hk - H0 + 4. / std::log(Ret)) *
+                          (Hk - H0 + 4. / std::log(Ret))) +
+                     0.015 / Hk) +
+                1.5 + 4. / Ret;
+
+    // Whitfield's minor additional compressibility correction.
+    Hstar = (Hstar + 0.028 * Me * Me) / (1. + 0.014 * Me * Me);
+
+    // Equivalent normalized wall slip velocity.
+    double us = (Hstar / 2.) * (1 - 4 * (Hk - 1.) / (3. * H));
+
+    double Ctcon = 0.5 / (6.7 * 6.7 * 0.75);
+
+    double Hkc = 0.;
+    double ctEq = 0.;
+
+    if (name == "wake")
+    {
+        if (us > 0.99995)
+            us = 0.99995;
+        Hkc = Hk - 1.;
+        ctEq = std::sqrt(4. * Hstar * Ctcon * ((Hk - 1.) * Hkc * Hkc) /
+                         ((1. - us) * (Hk * Hk) * H)); // Here it is ctEq^0.5
+    }
+    else
+    {
+        if (us > 0.95)
+            us = 0.98;
+        Hkc = std::max(Hk - 1. - 18. / Ret, 0.01);
+        ctEq = std::sqrt(Hstar * Ctcon * ((Hk - 1.) * Hkc * Hkc) /
+                         ((1. - us) * (Hk * Hk) * H)); // Here it is ctEq^0.5
+    }
+    if (Hstar * Ctcon * ((Hk - 1.) * Hkc * Hkc) / ((1. - us) * (Hk * Hk) * H) <
+        0.)
+        std::cout << "Negative sqrt encountered " << std::endl;
+
+    closureVars = ctEq;
 }
\ No newline at end of file

blast/src/DClosures.h

@@ -5,24 +5,26 @@
 #include <string>
 #include <vector>
 
-namespace blast {
+namespace blast
+{
 
 /**
  * @brief Boundary layer closure relations.
  * @author Paul Dechamps
  */
-class BLAST_API Closures {
+class BLAST_API Closures
+{
 
 public:
-  static void laminarClosures(std::vector<double> &closureVars, double theta,
-                              double H, double Ret, const double Me,
-                              const std::string name);
-  static void turbulentClosures(std::vector<double> &closureVars, double theta,
-                                double H, double Ct, double Ret, double Me,
-                                const std::string name);
-  static void turbulentClosures(double &closureVars, double theta, double H,
-                                double Ct, double Ret, const double Me,
+    static void laminarClosures(std::vector<double> &closureVars, double theta,
+                                double H, double Ret, const double Me,
                                 const std::string name);
+    static void turbulentClosures(std::vector<double> &closureVars, double theta,
+                                  double H, double Ct, double Ret, double Me,
+                                  const std::string name);
+    static void turbulentClosures(double &closureVars, double theta, double H,
+                                  double Ct, double Ret, const double Me,
+                                  const std::string name);
 };
 } // namespace blast
 #endif // DCLOSURES_H

blast/src/DDriver.h

@@ -5,65 +5,68 @@
 #include "blast.h"
 #include "wObject.h"
 #include "wTimers.h"
-namespace blast {
+namespace blast
+{
 
 /**
  * @brief Boundary layer solver driver. Performs the space marching and set up
  * time marching for each point.
  * @authors Paul Dechamps
  */
-class BLAST_API Driver : public fwk::wSharedObject {
+class BLAST_API Driver : public fwk::wSharedObject
+{
 public:
-  double Cdt;                  ///< Total drag coefficient.
-  double Cdf;                  ///< Total friction drag coefficient.
-  double Cdp;                  ///< Total pressure drag coefficient.
-  std::vector<double> Cdt_sec; ///< Sectional total drag coefficient.
-  std::vector<double> Cdf_sec; ///< Sectional friction drag coefficient.
-  std::vector<double> Cdp_sec; ///< Sectional pressure drag coefficient.
-  std::vector<std::vector<BoundaryLayer *>>
-      sections;         ///< Boundary layer regions sorted by section.
-  unsigned int verbose; ///< Verbosity level of the solver 0<=verbose<=1.
+    double Cdt;                  ///< Total drag coefficient.
+    double Cdf;                  ///< Total friction drag coefficient.
+    double Cdp;                  ///< Total pressure drag coefficient.
+    std::vector<double> Cdt_sec; ///< Sectional total drag coefficient.
+    std::vector<double> Cdf_sec; ///< Sectional friction drag coefficient.
+    std::vector<double> Cdp_sec; ///< Sectional pressure drag coefficient.
+    std::vector<std::vector<BoundaryLayer *>>
+        sections;         ///< Boundary layer regions sorted by section.
+    unsigned int verbose; ///< Verbosity level of the solver 0<=verbose<=1.
 
 private:
-  double Re;   ///< Freestream Reynolds number.
-  double CFL0; ///< Initial CFL number for pseudo time integration.
-  double span; ///< Wing Span (Used only for drag computation, not used if 2D
-               ///< case).
-  Solver *tSolver;    ///< Pseudo-time solver.
-  int solverExitCode; ///< Exit code of viscous calculations.
-  std::vector<std::vector<std::vector<size_t>>>
-      convergenceStatus; ///< Vector containing points that did not converged.
+    double Re;          ///< Freestream Reynolds number.
+    double CFL0;        ///< Initial CFL number for pseudo time integration.
+    double span;        ///< Wing Span (Used only for drag computation, not used if 2D
+                        ///< case).
+    Solver *tSolver;    ///< Pseudo-time solver.
+    int solverExitCode; ///< Exit code of viscous calculations.
+    std::vector<std::vector<std::vector<size_t>>>
+        convergenceStatus; ///< Vector containing points that did not converged.
 
 protected:
-  fwk::Timers tms; ///< Internal timers
+    fwk::Timers tms; ///< Internal timers
 
 public:
-  Driver(double _Re, double _Minf, double _CFL0, size_t nSections,
-         double xtrF_top = -1., double xtrF_bot = -1., double _span = 0.,
-         unsigned int _verbose = 1);
-  ~Driver();
-  int run();
-  void reset();
+    Driver(double _Re, double _Minf, double _CFL0, size_t nSections,
+           double xtrF_top = -1., double xtrF_bot = -1., double _span = 0.,
+           unsigned int _verbose = 1);
+    ~Driver();
+    int run();
+    void reset();
 
-  // Getters.
-  double getAverageTransition(size_t const iRegion) const;
-  double getRe() const { return Re; };
-  std::vector<std::vector<double>> getSolution(size_t iSec);
+    // Getters.
+    double getAverageTransition(size_t const iRegion) const;
+    double getRe() const { return Re; };
+    std::vector<std::vector<double>> getSolution(size_t iSec);
 
-  // Setters.
-  void setVerbose(unsigned int _verbose) { verbose = _verbose; };
+    // Setters.
+    void setVerbose(unsigned int _verbose) { verbose = _verbose; };
 
-  void addSection(size_t iSec, BoundaryLayer *&bl) {
-    sections[iSec].push_back(bl);
-  };
+    void addSection(size_t iSec, BoundaryLayer *&bl)
+    {
+        sections[iSec].push_back(bl);
+    };
 
 private:
-  void checkWaves(size_t iPoint, BoundaryLayer *bl);
-  void averageTransition(size_t iPoint, BoundaryLayer *bl, bool forced);
-  void computeSectionalDrag(std::vector<BoundaryLayer *> const &bl, size_t i);
-  void computeTotalDrag();
-  void computeBlwVel();
-  int outputStatus(double maxMach) const;
+    void checkWaves(size_t iPoint, BoundaryLayer *bl);
+    void averageTransition(size_t iPoint, BoundaryLayer *bl, bool forced);
+    void computeSectionalDrag(std::vector<BoundaryLayer *> const &bl, size_t i);
+    void computeTotalDrag();
+    void computeBlwVel();
+    int outputStatus(double maxMach) const;
 };
 } // namespace blast
 #endif // DDRIVER_H

blast/src/blast.h

@@ -34,7 +34,8 @@
 /**
  * @brief Namespace for blast module
  */
-namespace blast {
+namespace blast
+{
 
 // Solvers
 class Driver;