Skip to content
Snippets Groups Projects
Verified Commit 9681c193 authored by Paul Dechamps's avatar Paul Dechamps :speech_balloon:
Browse files

(refactor) Format code + Changed ints to doubles

parent b77fb320
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
blast/blCoupler.py
+ 1
− 1
View file @ 9681c193
......@@ -118,7 +118,7 @@ class Coupler:
if write:
self.isol.writeCp(sfx='_viscous'+self.filesfx)
for i, body in enumerate(self.vsol.bodies):
_sfx = f'_body{i}' if len(self.vsol.bodies) > 1 else ''
_sfx = f'_body{i}' if len(self.vsol.bodies) > 1 else ''
_ = vutils.getSolution(body.sections, write=True, toW='all', sfx=_sfx)
return aeroCoeffs
cdPrev = cd
......
blast/src/blBoundaryLayer.cpp
+ 10
− 2
View file @ 9681c193
......@@ -396,8 +396,8 @@ void BoundaryLayer::printSolution(size_t inod) const
throw std::runtime_error("Tried to access element outside of region size.");
std::cout << "Pt " << inod << "Reg " << name << " Turb " << nod->getRegime()
<< std::endl;
std::cout << std::setprecision(5) << "x " << nod->pos[0] << ", xoc "
<< nod->xoc << ", xx " << nod->xi << "xxExt "
std::cout << std::setprecision(5) << "x " << nod->pos[0] << ", y " << nod->pos[1] << ", xoc "
<< nod->xoc << ", xx " << nod->xi << " xxExt "
<< nod->xiExt << std::endl;
std::cout << std::scientific << std::setprecision(15);
std::cout << std::setw(10) << "T " << u[inod * nVar + 0] << "\n"
......@@ -407,6 +407,14 @@ void BoundaryLayer::printSolution(size_t inod) const
<< std::setw(10) << "Ct " << u[inod * nVar + 4] << "\n"
<< std::setw(10) << "dStar (BL) " << deltaStar[inod] << "\n"
<< std::setw(10) << "dStar (ext) " << deltaStarExt[inod] << "\n"
<< std::setw(10) << "Ret " << Ret[inod] << "\n"
<< std::setw(10) << "Hstar " << Hstar[inod] << "\n"
<< std::setw(10) << "Hstar2 " << Hstar2[inod] << "\n"
<< std::setw(10) << "Hk " << Hk[inod] << "\n"
<< std::setw(10) << "cfl " << cf[inod] << "\n"
<< std::setw(10) << "cdl " << cd[inod] << "\n"
<< std::setw(10) << "us " << us[inod] << "\n"
<< std::setw(10) << "delta " << delta[inod] << "\n"
<< std::setw(10) << "vt " << vt[inod] << "\n"
<< std::setw(10) << "Me " << Me[inod] << "\n"
<< std::setw(10) << "rhoe " << rhoe[inod] << std::endl;
......
blast/src/blBoundaryLayer.h
+ 1
− 1
View file @ 9681c193
......@@ -76,7 +76,7 @@ public:
size_t getnNodes() const { return nodes.size(); }
size_t getnElms() const { return nodes.size() - 1; }
double getxtr() const { return xtr; }
double getxoctr() const { return nodes[transitionNode]->xoc; }
double getxoctr() const { return nodes[transitionNode]->xoc; }
double getxitr() const { return nodes[transitionNode]->xi; }
double getxtrF() const { return xtrF; }
double getMaxMach() const { return *std::max_element(Me.begin(), Me.end()); }
......
blast/src/blClosures.cpp
+ 50
− 50
View file @ 9681c193
......@@ -21,7 +21,7 @@ void Closures::laminarClosures(std::vector<double> &closureVars, double theta,
H = std::max(H, 1.00005);
// Kinematic shape factor H*.
double Hk = (H - 0.29 * Me * Me) / (1 + 0.113 * Me * Me);
double Hk = (H - 0.29 * Me * Me) / (1.0 + 0.113 * Me * Me);
if (name == "wake")
Hk = std::max(Hk, 1.00005);
else
......@@ -31,23 +31,23 @@ void Closures::laminarClosures(std::vector<double> &closureVars, double theta,
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 delta = std::min((3.15 + H + (1.72 / (Hk - 1.0))) * theta, 12.0 * theta);
double Hstar = 0.;
double Hstar = 0.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 -
Hstar = 0.0111 * Hks * Hks / (Hk + 1.0) -
0.0278 * Hks * Hks * Hks / (Hk + 1.0) + 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);
Hstar = (Hstar + 0.028 * Me * Me) / (1.0 + 0.014 * Me * Me);
// Friction coefficient.
double tmp = 0.;
double cf = 0.;
double tmp = 0.0;
double cf = 0.0;
if (Hk < 5.5)
{
tmp = (5.5 - Hk) * (5.5 - Hk) * (5.5 - Hk) / (Hk + 1.0);
......@@ -60,22 +60,22 @@ void Closures::laminarClosures(std::vector<double> &closureVars, double theta,
}
// Dissipation coefficient.
double Cd = 0.;
if (Hk < 4)
Cd = (0.00205 * std::pow(4.0 - Hk, 5.5) + 0.207) * (Hstar / (2 * Ret));
double Cd = 0.0;
if (Hk < 4.0)
Cd = (0.00205 * std::pow(4.0 - Hk, 5.5) + 0.207) * (Hstar / (2.0 * Ret));
else
{
double HkCd = (Hk - 4) * (Hk - 4);
double denCd = 1 + 0.02 * HkCd;
Cd = (-0.0016 * HkCd / denCd + 0.207) * (Hstar / (2 * Ret));
double HkCd = (Hk - 4.0) * (Hk - 4.0);
double denCd = 1.0 + 0.02 * HkCd;
Cd = (-0.0016 * HkCd / denCd + 0.207) * (Hstar / (2.0 * 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.;
Cd = 2.0 * (1.10 * (1.0 - 1.0 / Hk) * (1.0 - 1.0 / Hk) / Hk) *
(Hstar / (2.0 * Ret));
cf = 0.0;
}
double us = 0.;
......@@ -184,25 +184,25 @@ void Closures::turbulentClosures(std::vector<double> &closureVars, double theta,
if (us > 0.95)
us = 0.98;
n = 1.0;
cf = (1 / Fc) *
cf = (1.0 / 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);
0.00011 * (std::tanh(4.0 - (Hk / 0.875)) - 1.0));
Hkc = std::max(Hk - 1.0 - 18.0 / Ret, 0.01);
// Dissipation coefficient.
Hmin = 1. + 2.1 / std::log(Ret);
Hmin = 1.0 + 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
ctEq = std::sqrt(Hstar * Ctcon * ((Hk - 1.0) * Hkc * Hkc) /
((1.0 - 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.0) * Hkc * Hkc) /
((1.0 - us) * (Hk * Hk) * H) <
0.0)
std::cout << "Negative sqrt encountered " << std::endl;
// Dissipation coefficient.
......@@ -234,57 +234,57 @@ void Closures::turbulentClosures(double &closureVars, double theta, double H,
else
Hk = std::max(Hk, 1.05000);
double H0 = 0.;
double H0 = 0.0;
if (Ret <= 400.)
H0 = 4.0;
else
H0 = 3. + (400. / Ret);
if (Ret <= 200.)
Ret = 200.;
H0 = 3.0 + (400.0 / Ret);
if (Ret <= 200.0)
Ret = 200.0;
double Hstar = 0.;
double Hstar = 0.0;
if (Hk <= H0)
Hstar =
((0.5 - 4. / Ret) * ((H0 - Hk) / (H0 - 1.)) * ((H0 - Hk) / (H0 - 1.))) *
((0.5 - 4.0 / Ret) * ((H0 - Hk) / (H0 - 1.)) * ((H0 - Hk) / (H0 - 1.))) *
(1.5 / (Hk + 0.5)) +
1.5 + 4. / Ret;
1.5 + 4.0 / 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.0 / std::log(Ret)) *
(Hk - H0 + 4.0 / std::log(Ret))) +
0.015 / Hk) +
1.5 + 4. / Ret;
1.5 + 4.0 / 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 + 0.014 * Me * Me);
// Equivalent normalized wall slip velocity.
double us = (Hstar / 2.) * (1 - 4 * (Hk - 1.) / (3. * H));
double us = (Hstar / 2.0) * (1.0 - 4 * (Hk - 1.0) / (3.0 * H));
double Ctcon = 0.5 / (6.7 * 6.7 * 0.75);
double Hkc = 0.;
double ctEq = 0.;
double Hkc = 0.0;
double ctEq = 0.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
Hkc = Hk - 1.0;
ctEq = std::sqrt(4.0 * Hstar * Ctcon * ((Hk - 1.0) * Hkc * Hkc) /
((1.0 - 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.0 - 18.0 / Ret, 0.01);
ctEq = std::sqrt(Hstar * Ctcon * ((Hk - 1.0) * Hkc * Hkc) /
((1.0 - 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.0) * Hkc * Hkc) / ((1.0 - us) * (Hk * Hk) * H) <
0.0)
std::cout << "Negative sqrt encountered " << std::endl;
closureVars = ctEq;
......
blast/src/blDriver.cpp
+ 17
− 6
View file @ 9681c193
......@@ -69,7 +69,7 @@ int Driver::run()
solverExitCode = 0;
if (verbose >= 1)
std::cout << "Solving boundary layer for Re " << std::setprecision(1) << Re/1e6 << "e6, " << bodies.size() << " bodies." << std::endl;
std::cout << "Solving boundary layer for Re " << std::setprecision(1) << Re / 1e6 << "e6, " << bodies.size() << " bodies." << std::endl;
for (auto &body : bodies)
{
......@@ -114,7 +114,7 @@ int Driver::run()
lockTrans = true;
}
else
throw std::runtime_error("Didn't expect region name: " + reg->getName());
throw std::runtime_error("Didn't expect region name: " + reg->getName());
// Loop over points
for (size_t inod = 1; inod < reg->getnNodes(); ++inod)
......@@ -123,7 +123,12 @@ int Driver::run()
tSolver->initialCondition(inod, reg);
// Solve equations
pointExitCode = tSolver->integration(inod, reg);
if (pointExitCode != 0){solverExitCode = -1; numberFailedPerBody++;}
if (pointExitCode != 0)
{
if (solverExitCode != -1)
solverExitCode = pointExitCode;
numberFailedPerBody++;
}
// Transition
if (!lockTrans)
......@@ -166,11 +171,17 @@ int Driver::run()
{
std::cout << "Maximum Mach number: " << maxMach << std::endl;
if (solverExitCode == 0)
std::cout << "\x1b[1;32m" << "Boundary layer solver converged." << "\x1b[0m" << std::endl;
std::cout << "\x1b[1;32m"
<< "Boundary layer solver converged."
<< "\x1b[0m" << std::endl;
else if (solverExitCode == -1)
std::cout << "\x1b[1;31m" << "Boundary layer solver diverged!" << "\x1b[0m" << std::endl;
std::cout << "\x1b[1;31m"
<< "Boundary layer solver diverged!"
<< "\x1b[0m" << std::endl;
else if (solverExitCode != 0)
std::cout << "\x1b[1;33m" << "Boundary layer solver not fully converged." << "\x1b[0m" << std::endl;
std::cout << "\x1b[1;33m"
<< "Boundary layer solver not fully converged."
<< "\x1b[0m" << std::endl;
std::cout << std::endl;
}
return solverExitCode; // exit code
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment