From eac86b002964c6c34a7ee6fff1e788e164c401b7 Mon Sep 17 00:00:00 2001
From: Paul Dechamps <paul.dechamps@uliege.be>
Date: Tue, 12 Nov 2024 10:39:51 +0100
Subject: [PATCH] (fix) Changed int to double
---
blast/src/DBoundaryLayer.cpp | 8 ++--
blast/src/DClosures.cpp | 82 ++++++++++++++++++------------------
blast/src/DFluxes.cpp | 28 ++++++------
blast/src/DSolver.cpp | 6 +--
4 files changed, 62 insertions(+), 62 deletions(-)
diff --git a/blast/src/DBoundaryLayer.cpp b/blast/src/DBoundaryLayer.cpp
index cb81c81..148ed29 100644
--- a/blast/src/DBoundaryLayer.cpp
+++ b/blast/src/DBoundaryLayer.cpp
@@ -76,7 +76,7 @@ void BoundaryLayer::computeLocalCoordinates() {
}
// Compute the x/c coordinate.
- if (chord <= 0)
+ if (chord <= 0.)
throw std::runtime_error("Chord length is zero or negative.");
xoc.resize(nNodes, 0.);
for (size_t iPoint = 0; iPoint < nNodes; ++iPoint)
@@ -127,7 +127,7 @@ void BoundaryLayer::setStagnationBC(double Re) {
u[4] = 0.; // Ctau
Ret[0] = u[3] * u[0] * Re;
- if (Ret[0] < 1) {
+ if (Ret[0] < 1.) {
Ret[0] = 1.;
u[3] = Ret[0] / (u[0] * Re);
}
@@ -168,8 +168,8 @@ void BoundaryLayer::setWakeBC(double Re, std::vector<double> UpperCond,
Ret[0] = u[3] * u[0] * Re; // Reynolds number based on the momentum thickness.
- if (Ret[0] < 1) {
- Ret[0] = 1;
+ if (Ret[0] < 1.) {
+ Ret[0] = 1.;
u[3] = Ret[0] / (u[0] * Re);
}
deltaStar[0] = u[0] * u[1];
diff --git a/blast/src/DClosures.cpp b/blast/src/DClosures.cpp
index 53e872d..304e499 100644
--- a/blast/src/DClosures.cpp
+++ b/blast/src/DClosures.cpp
@@ -96,7 +96,7 @@ void Closures::turbulentClosures(std::vector<double> &closureVars, double theta,
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);
+ double Hk = (H - 0.29 * Me * Me) / (1. + 0.113 * Me * Me);
if (name == "wake")
Hk = std::max(Hk, 1.00005);
else
@@ -104,42 +104,42 @@ void Closures::turbulentClosures(std::vector<double> &closureVars, double theta,
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 Fc = std::sqrt(1. + 0.5 * gamma * Me * Me);
double H0 = 0.;
- if (Ret <= 400)
+ if (Ret <= 400.)
H0 = 4.0;
else
- H0 = 3 + (400 / Ret);
- if (Ret <= 200)
- Ret = 200;
+ 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))) *
+ ((0.5 - 4. / Ret) * ((H0 - Hk) / (H0 - 1.)) * ((H0 - Hk) / (H0 - 1.))) *
(1.5 / (Hk + 0.5)) +
- 1.5 + 4 / Ret;
+ 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))) +
+ ((Hk - H0 + 4. / std::log(Ret)) *
+ (Hk - H0 + 4. / std::log(Ret))) +
0.015 / Hk) +
- 1.5 + 4 / Ret;
+ 1.5 + 4. / Ret;
// Whitfield's minor additional compressibility correction.
- Hstar = (Hstar + 0.028 * Me * Me) / (1 + 0.014 * Me * Me);
+ 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));
+ 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 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.;
@@ -163,13 +163,13 @@ void Closures::turbulentClosures(std::vector<double> &closureVars, double theta,
us = 0.99995;
n = 2.0; // two wake halves
cf = 0.0; // no friction inside the wake
- Hkc = Hk - 1;
+ 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.
+ 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;
@@ -179,19 +179,19 @@ void Closures::turbulentClosures(std::vector<double> &closureVars, double theta,
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);
+ 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 * 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)
+ if (n * Hstar * Ctcon * ((Hk - 1.) * Hkc * Hkc) / ((1. - us) * (Hk * Hk) * H) <
+ 0.)
std::cout << "Negative sqrt encountered " << std::endl;
// Dissipation coefficient.
@@ -223,32 +223,32 @@ void Closures::turbulentClosures(double &closureVars, double theta, double H,
Hk = std::max(Hk, 1.05000);
double H0 = 0.;
- if (Ret <= 400)
+ if (Ret <= 400.)
H0 = 4.0;
else
- H0 = 3 + (400 / Ret);
- if (Ret <= 200)
- Ret = 200;
+ 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))) *
+ ((0.5 - 4. / Ret) * ((H0 - Hk) / (H0 - 1.)) * ((H0 - Hk) / (H0 - 1.))) *
(1.5 / (Hk + 0.5)) +
- 1.5 + 4 / Ret;
+ 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))) +
+ ((Hk - H0 + 4. / std::log(Ret)) *
+ (Hk - H0 + 4. / std::log(Ret))) +
0.015 / Hk) +
- 1.5 + 4 / Ret;
+ 1.5 + 4. / Ret;
// Whitfield's minor additional compressibility correction.
- Hstar = (Hstar + 0.028 * Me * Me) / (1 + 0.014 * Me * Me);
+ 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 us = (Hstar / 2.) * (1 - 4 * (Hk - 1.) / (3. * H));
double Ctcon = 0.5 / (6.7 * 6.7 * 0.75);
@@ -258,17 +258,17 @@ void Closures::turbulentClosures(double &closureVars, double theta, double H,
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
+ 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
+ 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)
+ if (Hstar * Ctcon * ((Hk - 1.) * Hkc * Hkc) / ((1. - us) * (Hk * Hk) * H) < 0.)
std::cout << "Negative sqrt encountered " << std::endl;
closureVars = ctEq;
diff --git a/blast/src/DFluxes.cpp b/blast/src/DFluxes.cpp
index 4f33a13..cb9e69d 100644
--- a/blast/src/DFluxes.cpp
+++ b/blast/src/DFluxes.cpp
@@ -117,8 +117,8 @@ VectorXd Fluxes::blLaws(size_t iPoint, BoundaryLayer *bl,
double ddeltaStarExt_dx =
(bl->deltaStarExt[iPoint] - bl->deltaStarExt[iPoint - 1]) / dxExt;
- double c = 4 / (M_PI * dx);
- double cExt = 4 / (M_PI * dxExt);
+ double c = 4. / (M_PI * dx);
+ double cExt = 4. / (M_PI * dxExt);
// Amplification routine.
double dN_dx = 0.0;
@@ -295,8 +295,8 @@ double Fluxes::amplificationRoutine(double Hk, double Ret, double theta) const {
double Hk2 = 4.;
double Hmi = (1 / (Hk - 1));
double logRet_crit =
- 2.492 * std::pow(1 / (Hk - 1.), 0.43) +
- 0.7 * (std::tanh(14 * Hmi - 9.24) +
+ 2.492 * std::pow(1. / (Hk - 1.), 0.43) +
+ 0.7 * (std::tanh(14. * Hmi - 9.24) +
1.0); // value at which the amplification starts to grow
double ax = 0.;
@@ -305,37 +305,37 @@ double Fluxes::amplificationRoutine(double Hk, double Ret, double theta) const {
if (std::log10(Ret) < logRet_crit - dgr)
ax = 0.;
else {
- double rNorm = (std::log10(Ret) - (logRet_crit - dgr)) / (2 * dgr);
+ double rNorm = (std::log10(Ret) - (logRet_crit - dgr)) / (2. * dgr);
double Rfac = 0.;
if (rNorm <= 0)
Rfac = 0.;
if (rNorm >= 1)
Rfac = 1.;
else
- Rfac = 3 * rNorm * rNorm - 2 * rNorm * rNorm * rNorm;
+ Rfac = 3. * rNorm * rNorm - 2. * rNorm * rNorm * rNorm;
// Normal envelope amplification rate
- double m = -0.05 + 2.7 * Hmi - 5.5 * Hmi * Hmi + 3 * Hmi * Hmi * Hmi +
- 0.1 * std::exp(-20 * Hmi);
+ double m = -0.05 + 2.7 * Hmi - 5.5 * Hmi * Hmi + 3. * Hmi * Hmi * Hmi +
+ 0.1 * std::exp(-20. * Hmi);
double arg = 3.87 * Hmi - 2.52;
- double dndRet = 0.028 * (Hk - 1) - 0.0345 * std::exp(-(arg * arg));
+ double dndRet = 0.028 * (Hk - 1.) - 0.0345 * std::exp(-(arg * arg));
ax = (m * dndRet / theta) * Rfac;
// Set correction for separated profiles
if (Hk > Hk1) {
double Hnorm = (Hk - Hk1) / (Hk2 - Hk1);
double Hfac = 0.;
- if (Hnorm >= 1)
+ if (Hnorm >= 1.)
Hfac = 1.;
else
- Hfac = 3 * Hnorm * Hnorm - 2 * Hnorm * Hnorm * Hnorm;
+ Hfac = 3. * Hnorm * Hnorm - 2. * Hnorm * Hnorm * Hnorm;
double ax1 = ax;
- double gro = 0.3 + 0.35 * std::exp(-0.15 * (Hk - 5));
+ double gro = 0.3 + 0.35 * std::exp(-0.15 * (Hk - 5.));
double tnr = std::tanh(1.2 * (std::log10(Ret) - gro));
- double ax2 = (0.086 * tnr - 0.25 / std::pow((Hk - 1), 1.5)) / theta;
+ double ax2 = (0.086 * tnr - 0.25 / std::pow((Hk - 1.), 1.5)) / theta;
if (ax2 < 0.)
ax2 = 0.;
- ax = Hfac * ax2 + (1 - Hfac) * ax1;
+ ax = Hfac * ax2 + (1. - Hfac) * ax1;
if (ax < 0.)
ax = 0.;
}
diff --git a/blast/src/DSolver.cpp b/blast/src/DSolver.cpp
index 6c7b39a..38f46c1 100644
--- a/blast/src/DSolver.cpp
+++ b/blast/src/DSolver.cpp
@@ -60,7 +60,7 @@ void Solver::initialCondition(size_t iPoint, BoundaryLayer *bl) {
bl->us[iPoint] = bl->us[iPoint - 1];
bl->delta[iPoint] = bl->delta[iPoint - 1];
bl->deltaStar[iPoint] = bl->deltaStar[iPoint - 1];
- if (bl->regime[iPoint] == 1 && bl->u[iPoint * nVar + 4] <= 0)
+ if (bl->regime[iPoint] == 1 && bl->u[iPoint * nVar + 4] <= 0.)
bl->u[iPoint * nVar + 4] = 0.03;
}
@@ -147,7 +147,7 @@ double Solver::setTimeStep(double CFL, double dx, double invVel) const {
double gradU2 = (invVel * invVel);
double soundSpeed = computeSoundSpeed(gradU2);
- // Velocity of the fastest travelling wave.
+ // Velocity of the fastest traveling wave.
double advVel = soundSpeed + invVel;
// Time step computation.
@@ -161,7 +161,7 @@ double Solver::setTimeStep(double CFL, double dx, double invVel) const {
* @return double
*/
double Solver::computeSoundSpeed(double gradU2) const {
- return sqrt(1 / (Minf * Minf) + 0.5 * (gamma - 1) * (1 - gradU2));
+ return sqrt(1. / (Minf * Minf) + 0.5 * (gamma - 1.) * (1. - gradU2));
}
/**
--
GitLab