From 382639f5d6be7c033847a65c4ad2a6b91a247cb1 Mon Sep 17 00:00:00 2001
From: Paul Dechamps <paul.dechamps@uliege.be>
Date: Mon, 11 Mar 2024 09:33:11 +0100
Subject: [PATCH] (fix) Small modifs
Ensure all numbers are doubles in Closures, Fluxes and Driver & debug info update in BoundaryLayer & ensure nIter in 3D rae test
---
blast/interfaces/dart/DartAdjointInterface.py | 11 ++-
blast/src/DBoundaryLayer.cpp | 13 +++-
blast/src/DClosures.cpp | 78 +++++++++----------
blast/src/DDriver.cpp | 21 ++---
blast/src/DFluxes.cpp | 37 +++++----
blast/tests/dart/rae2822_3D.py | 4 +-
6 files changed, 90 insertions(+), 74 deletions(-)
diff --git a/blast/interfaces/dart/DartAdjointInterface.py b/blast/interfaces/dart/DartAdjointInterface.py
index f21e495..f8ece62 100644
--- a/blast/interfaces/dart/DartAdjointInterface.py
+++ b/blast/interfaces/dart/DartAdjointInterface.py
@@ -55,14 +55,19 @@ class DartAdjointInterface():
blowing = self.__runViscous()
dRblw_v[:,iNo*nDim + iDim] = (blowing-saveBlw)/delta
self.iSolverAPI.solver.U[iNo][iDim] = savev[iNo*nDim + iDim]
-
+ # dRblw_M[abs(dRblw_M) < 1e-10] = 0.
+ # dRblw_rho[abs(dRblw_rho) < 1e-10] = 0.
+ # dRblw_v[abs(dRblw_v) < 1e-10] = 0.
+ self.coupledAdjointSolver.setdRblw_M(dRblw_M)
+ self.coupledAdjointSolver.setdRblw_rho(dRblw_rho)
+ self.coupledAdjointSolver.setdRblw_v(dRblw_v)
def __runViscous(self):
- self.iSolverAPI.imposeInvBC(0)
+ self.iSolverAPI.imposeInvBC(1)
# Run viscous solver.
- exitCode = self.vSolver.run(0)
+ exitCode = self.vSolver.run(1)
# Impose blowing boundary condition in the inviscid solver.
self.iSolverAPI.imposeBlwVel()
diff --git a/blast/src/DBoundaryLayer.cpp b/blast/src/DBoundaryLayer.cpp
index 804d780..dae9606 100644
--- a/blast/src/DBoundaryLayer.cpp
+++ b/blast/src/DBoundaryLayer.cpp
@@ -70,6 +70,8 @@ void BoundaryLayer::setMesh(std::vector<double> x,
*/
void BoundaryLayer::setStagBC(double Re)
{
+ if (name != "upper" && name != "lower")
+ std::cout << "Warning: Imposing stagnation boundary condition on " << name << std::endl;
double dv0 = (blEdge->getVt(1) - blEdge->getVt(0)) / (mesh->getLoc(1) - mesh->getLoc(0));
if (dv0 > 0)
u[0] = sqrt(0.075 / (Re * dv0)); // Theta
@@ -196,5 +198,14 @@ void BoundaryLayer::printSolution(size_t iPoint) const
<< std::setw(10) << "dStar (old) " << blEdge->getDeltaStar(iPoint) << "\n"
<< std::setw(10) << "vt " << blEdge->getVt(iPoint) << "\n"
<< std::setw(10) << "Me " << blEdge->getMe(iPoint) << "\n"
- << std::setw(10) << "rhoe " << blEdge->getRhoe(iPoint) << std::endl;
+ << std::setw(10) << "rhoe " << blEdge->getRhoe(iPoint) << "\n"
+ << std::setw(10) << "Hstar " << Hstar[iPoint] << "\n"
+ << std::setw(10) << "Hstar2 " << Hstar2[iPoint] << "\n"
+ << std::setw(10) << "Hk " << Hk[iPoint] << "\n"
+ << std::setw(10) << "cf " << cf[iPoint] << "\n"
+ << std::setw(10) << "cd " << cd[iPoint] << "\n"
+ << std::setw(10) << "delta " << delta[iPoint] << "\n"
+ << std::setw(10) << "ctEq " << ctEq[iPoint] << "\n"
+ << std::setw(10) << "us " << us[iPoint] << "\n"
+ << std::setw(10) << "Ret " << Ret[iPoint] << "\n" << std::endl;
}
\ No newline at end of file
diff --git a/blast/src/DClosures.cpp b/blast/src/DClosures.cpp
index c411d1c..0d384a5 100644
--- a/blast/src/DClosures.cpp
+++ b/blast/src/DClosures.cpp
@@ -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.113 * Me * Me);
if (name == "wake")
Hk = std::max(Hk, 1.00005);
else
@@ -31,7 +31,7 @@ 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.))) * theta, 12. * theta);
double Hstar = 0.;
double Hks = Hk - 4.35;
@@ -41,7 +41,7 @@ void Closures::laminarClosures(std::vector<double> &closureVars, double theta,
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.014 * Me * Me);
// Friction coefficient.
double tmp = 0.;
@@ -60,18 +60,18 @@ 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));
+ Cd = (0.00205 * std::pow(4.0 - Hk, 5.5) + 0.207) * (Hstar / (2. * Ret));
else
{
- double HkCd = (Hk - 4) * (Hk - 4);
+ double HkCd = (Hk - 4.) * (Hk - 4.);
double denCd = 1 + 0.02 * HkCd;
- Cd = (-0.0016 * HkCd / denCd + 0.207) * (Hstar / (2 * Ret));
+ 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));
+ Cd = 2. * (1.10 * (1.0 - 1.0 / Hk) * (1.0 - 1.0 / Hk) / Hk) * (Hstar / (2. * Ret));
cf = 0.;
}
@@ -98,7 +98,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
@@ -106,34 +106,34 @@ 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))) * (1.5 / (Hk + 0.5)) + 1.5 + 4 / Ret;
+ 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;
+ 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);
+ 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.;
@@ -158,30 +158,30 @@ 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;
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);
+ 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);
+ 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.
@@ -205,31 +205,31 @@ void Closures::turbulentClosures(double &closureVars, double theta, double H, do
H = std::max(H, 1.00005);
Ct = std::min(Ct, 0.3);
- 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
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))) * (1.5 / (Hk + 0.5)) + 1.5 + 4 / Ret;
+ 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;
+ 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 us = (Hstar / 2.) * (1. - 4. * (Hk - 1.) / (3. * H));
double Ctcon = 0.5 / (6.7 * 6.7 * 0.75);
@@ -240,17 +240,17 @@ void Closures::turbulentClosures(double &closureVars, double theta, double H, do
{
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/DDriver.cpp b/blast/src/DDriver.cpp
index b937b4f..d41f579 100644
--- a/blast/src/DDriver.cpp
+++ b/blast/src/DDriver.cpp
@@ -161,7 +161,7 @@ int Driver::run(unsigned int const couplIter)
else
{
tSolver->setCFL0(CFL0);
- tSolver->setitFrozenJac(5);
+ tSolver->setitFrozenJac(1);
tSolver->setinitSln(true);
// Keyword annoncing iteration safety level.
@@ -262,15 +262,15 @@ int Driver::run(unsigned int const couplIter)
void Driver::checkWaves(size_t iPoint, BoundaryLayer *bl)
{
- double logRet_crit = 2.492 * std::pow(1 / (bl->Hk[iPoint] - 1), 0.43) + 0.7 * (std::tanh(14 * (1 / (bl->Hk[iPoint] - 1)) - 9.24) + 1.0);
+ double logRet_crit = 2.492 * std::pow(1. / (bl->Hk[iPoint] - 1.), 0.43) + 0.7 * (std::tanh(14. * (1. / (bl->Hk[iPoint] - 1)) - 9.24) + 1.0);
- if (bl->Ret[iPoint] > 0) // Avoid log of negative number.
+ if (bl->Ret[iPoint] > 0.) // Avoid log of negative number.
{
if (std::log10(bl->Ret[iPoint]) <= logRet_crit - 0.08)
- bl->u[iPoint * bl->getnVar() + 2] = 0; // Set N to 0 at that point if waves do not grow.
+ bl->u[iPoint * bl->getnVar() + 2] = 0.; // Set N to 0 at that point if waves do not grow.
}
else
- bl->u[iPoint * bl->getnVar() + 2] = 0; // Set N to 0 at that point if waves do not grow.
+ bl->u[iPoint * bl->getnVar() + 2] = 0.; // Set N to 0 at that point if waves do not grow.
}
/**
@@ -323,11 +323,12 @@ void Driver::averageTransition(size_t iPoint, BoundaryLayer *bl, int forced)
std::cout << "Warning: Transition point turbulent computation terminated with exit code " << exitCode << std::endl;
// Average both solutions (Now solution @ iPoint is the turbulent solution).
- bl->u[iPoint * nVar + 0] = avLam * lamSol[0] + avTurb * bl->u[(iPoint)*nVar + 0]; // Theta.
- bl->u[iPoint * nVar + 1] = avLam * lamSol[1] + avTurb * bl->u[(iPoint)*nVar + 1]; // H.
- bl->u[iPoint * nVar + 2] = 9.; // N.
- bl->u[iPoint * nVar + 3] = avLam * lamSol[3] + avTurb * bl->u[(iPoint)*nVar + 3]; // ue.
- bl->u[iPoint * nVar + 4] = avTurb * bl->u[(iPoint)*nVar + 4]; // Ct.
+ // bl->u[iPoint * nVar + 0] = avLam * lamSol[0] + avTurb * bl->u[(iPoint)*nVar + 0]; // Theta.
+ // bl->u[iPoint * nVar + 1] = avLam * lamSol[1] + avTurb * bl->u[(iPoint)*nVar + 1]; // H.
+ // bl->u[iPoint * nVar + 2] = 9.; // N.
+ // bl->u[iPoint * nVar + 3] = avLam * lamSol[3] + avTurb * bl->u[(iPoint)*nVar + 3]; // ue.
+ // bl->u[iPoint * nVar + 4] = avTurb * bl->u[(iPoint)*nVar + 4]; // Ct.
+ bl->u[iPoint*nVar+2] = 9.;
// Recompute closures. (The turbulent BC @ iPoint - 1 does not influence laminar closure relations).
std::vector<double> lamParam(8, 0.);
diff --git a/blast/src/DFluxes.cpp b/blast/src/DFluxes.cpp
index c3beca0..d159130 100644
--- a/blast/src/DFluxes.cpp
+++ b/blast/src/DFluxes.cpp
@@ -147,23 +147,22 @@ VectorXd Fluxes::blLaws(size_t iPoint, BoundaryLayer *bl, std::vector<double> u)
double usa = bl->us[iPoint];
// Space part.
- spaceVector(0) = dt_dx + (2 + u[1] - Mea * Mea) * (u[0] / u[3]) * due_dx - cfa / 2;
- spaceVector(1) = u[0] * dHstar_dx + (2 * Hstar2a + Hstara * (1 - u[1])) * u[0] / u[3] * due_dx - 2 * cda + Hstara * cfa / 2;
+ spaceVector(0) = dt_dx + (2. + u[1] - Mea * Mea) * (u[0] / u[3]) * due_dx - cfa / 2.;
+ spaceVector(1) = u[0] * dHstar_dx + (2 * Hstar2a + Hstara * (1 - u[1])) * u[0] / u[3] * due_dx - 2. * cda + Hstara * cfa / 2.;
spaceVector(2) = dN_dx - ax;
spaceVector(3) = due_dx - c * (u[1] * dt_dx + u[0] * dH_dx) - dueExt_dx + cExt * ddeltaStarExt_dx;
if (bl->regime[iPoint] == 1)
- spaceVector(4) = (2 * deltaa / u[4]) * dct_dx - 5.6 * (ctEqa - u[4] * dissipRatio) - 2 * deltaa * (4 / (3 * u[1] * u[0]) * (cfa / 2 - (((Hka - 1) / (6.7 * Hka * dissipRatio)) * ((Hka - 1) / (6.7 * Hka * dissipRatio)))) - 1 / u[3] * due_dx);
+ spaceVector(4) = (2. * deltaa / u[4]) * dct_dx - 5.6 * (ctEqa - u[4] * dissipRatio) - 2. * deltaa * (4. / (3. * u[1] * u[0]) * (cfa / 2. - (((Hka - 1.) / (6.7 * Hka * dissipRatio)) * ((Hka - 1.) / (6.7 * Hka * dissipRatio)))) - 1. / u[3] * due_dx);
// Time part.
timeMatrix(0, 0) = u[1] / u[3];
timeMatrix(0, 1) = u[0] / u[3];
timeMatrix(0, 3) = u[0] * u[1] / (u[3] * u[3]);
- timeMatrix(1, 0) = (1 + u[1] * (1 - Hstara)) / u[3];
- timeMatrix(1, 1) = (1 - Hstara) * u[0] / u[3];
- timeMatrix(1, 3) = (2 - Hstara * u[1]) * u[0] / (u[3] * u[3]);
-
+ timeMatrix(1, 0) = (1. + u[1] * (1. - Hstara)) / u[3];
+ timeMatrix(1, 1) = (1. - Hstara) * u[0] / u[3];
+ timeMatrix(1, 3) = (2. - Hstara * u[1]) * u[0] / (u[3] * u[3]);
timeMatrix(2, 2) = 1.;
timeMatrix(3, 0) = -c * u[1];
@@ -172,8 +171,8 @@ VectorXd Fluxes::blLaws(size_t iPoint, BoundaryLayer *bl, std::vector<double> u)
if (bl->regime[iPoint] == 1)
{
- timeMatrix(4, 3) = 2 * deltaa / (usa * u[3] * u[3]);
- timeMatrix(4, 4) = 2 * deltaa / (u[3] * u[4] * usa);
+ timeMatrix(4, 3) = 2. * deltaa / (usa * u[3] * u[3]);
+ timeMatrix(4, 4) = 2. * deltaa / (u[3] * u[4] * usa);
}
else
timeMatrix(4, 4) = 1.0;
@@ -194,7 +193,7 @@ double Fluxes::amplificationRoutine(double Hk, double Ret, double theta) const
double dgr = 0.08;
double Hk1 = 3.5;
double Hk2 = 4.;
- double Hmi = (1 / (Hk - 1));
+ 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) + 1.0); // value at which the amplification starts to grow
double ax = 0.;
@@ -204,19 +203,19 @@ double Fluxes::amplificationRoutine(double Hk, double Ret, double theta) const
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
@@ -224,17 +223,17 @@ double Fluxes::amplificationRoutine(double Hk, double Ret, double theta) const
{
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/tests/dart/rae2822_3D.py b/blast/tests/dart/rae2822_3D.py
index 92add20..1b01e12 100644
--- a/blast/tests/dart/rae2822_3D.py
+++ b/blast/tests/dart/rae2822_3D.py
@@ -105,8 +105,8 @@ def cfgBlast(verb):
'saveTag': 4,
'Verb': verb, # Verbosity level of the solver
- 'couplIter': 50, # Maximum number of iterations
- 'couplTol' : 5e-2, # Tolerance of the VII methodology
+ 'couplIter': 5, # Maximum number of iterations
+ 'couplTol' : 1e-4, # Tolerance of the VII methodology
'iterPrint': 1, # int, number of iterations between outputs
'resetInv' : False, # bool, flag to reset the inviscid calculation at every iteration.
'xtrF' : [0., 0.], # Forced transition location
--
GitLab