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

(fix) Small modifs 

Ensure all numbers are doubles in Closures, Fluxes and Driver & debug info update in BoundaryLayer & ensure nIter in 3D rae test
parent 71e1331b
No related branches found
No related tags found
No related merge requests found
Pipeline #21491 failed
Stage: build
    Stage: test
      Hide whitespace changes
      Inline Side-by-side
      blast/interfaces/dart/DartAdjointInterface.py
      + 8
      − 3
      View file @ 382639f5
      ......@@ -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()
      ......
      blast/src/DBoundaryLayer.cpp
      + 12
      − 1
      View file @ 382639f5
      ......@@ -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
      blast/src/DClosures.cpp
      + 39
      − 39
      View file @ 382639f5
      ......@@ -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;
      ......
      blast/src/DDriver.cpp
      + 11
      − 10
      View file @ 382639f5
      ......@@ -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.);
      ......
      blast/src/DFluxes.cpp
      + 18
      − 19
      View file @ 382639f5
      ......@@ -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.;
      }
      ......
      blast/tests/dart/rae2822_3D.py
      + 2
      − 2
      View file @ 382639f5
      ......@@ -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
      ......
      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