-
Paul Dechamps authored
Change potential and mesh derivatives to analytical. Also changed the residuals computation of the viscous solver to the analytical expression instead of solving the linear system. This commit reduces the computational cost of the coupled adjoint methodology by a factor 4, at least.
Paul Dechamps authoredChange potential and mesh derivatives to analytical. Also changed the residuals computation of the viscous solver to the analytical expression instead of solving the linear system. This commit reduces the computational cost of the coupled adjoint methodology by a factor 4, at least.
DFluxes.cpp 11.15 KiB
#include "DFluxes.h"
#include "DBoundaryLayer.h"
#include <Eigen/Dense>
using namespace blast;
using namespace Eigen;
/**
* @brief Compute residual vector of the integral boundary layer equations.
*
* @param iPoint Marker id.
* @param bl Boundary layer region.
* @return VectorXd
*/
VectorXd Fluxes::computeResiduals(size_t iPoint, BoundaryLayer *bl)
{
size_t nVar = bl->getnVar();
std::vector<double> u(nVar, 0.);
for (size_t k = 0; k < nVar; ++k)
u[k] = bl->u[iPoint * nVar + k];
return blLaws(iPoint, bl, u);
}
/**
* @brief Compute Jacobian of the integral boundary layer equations.
*
* @param iPoint Marker id.
* @param bl Boundary layer region.
* @param sysRes Residual of the IBL at point iPoint.
* @param eps Pertubation from current solution used to compute the Jacobian.
* @return MatrixXd
*/
MatrixXd Fluxes::computeJacobian(size_t iPoint, BoundaryLayer *bl, VectorXd const &sysRes, double eps)
{
size_t nVar = bl->getnVar();
MatrixXd JacMatrix = MatrixXd::Zero(5, 5);
std::vector<double> uUp(nVar, 0.);
for (size_t k = 0; k < nVar; ++k)
uUp[k] = bl->u[iPoint * nVar + k];
double varSave = 0.;
for (size_t iVar = 0; iVar < nVar; ++iVar)
{
varSave = uUp[iVar];
uUp[iVar] += eps;
JacMatrix.col(iVar) = (blLaws(iPoint, bl, uUp) - sysRes) / eps;
uUp[iVar] = varSave;
}
return JacMatrix;
}
/**
* @brief Integral boundary layer equation.
*
* @param iPoint Marker id.
* @param bl Boundary layer region.
* @param u Solution vector at point iPoint.
* @return VectorXd
*/
VectorXd Fluxes::blLaws(size_t iPoint, BoundaryLayer *bl, std::vector<double> u)
{
size_t nVar = bl->getnVar();
// Dissipation ratio.
double dissipRatio = 0.;
if (bl->name == "wake")
dissipRatio = 0.9;
else
dissipRatio = 1.;
bl->Ret[iPoint] = std::max(u[3] * u[0] * Re, 1.0); // Reynolds number based on theta.
bl->deltaStar[iPoint] = u[0] * u[1]; // Displacement thickness.
// Boundary layer closure relations.
if (bl->regime[iPoint] == 0)
{
// Laminar closure relations.
std::vector<double> lamParam(8, 0.);
bl->closSolver->laminarClosures(lamParam, u[0], u[1], bl->Ret[iPoint], bl->Me[iPoint], bl->name);
bl->Hstar[iPoint] = lamParam[0];
bl->Hstar2[iPoint] = lamParam[1];
bl->Hk[iPoint] = lamParam[2];
bl->cf[iPoint] = lamParam[3];
bl->cd[iPoint] = lamParam[4];
bl->delta[iPoint] = lamParam[5];
bl->ctEq[iPoint] = lamParam[6];
bl->us[iPoint] = lamParam[7];
}
else if (bl->regime[iPoint] == 1)
{
// Turbulent closure relations.
std::vector<double> turbParam(8, 0.);
bl->closSolver->turbulentClosures(turbParam, u[0], u[1], u[4], bl->Ret[iPoint], bl->Me[iPoint], bl->name);
bl->Hstar[iPoint] = turbParam[0];
bl->Hstar2[iPoint] = turbParam[1];
bl->Hk[iPoint] = turbParam[2];
bl->cf[iPoint] = turbParam[3];
bl->cd[iPoint] = turbParam[4];
bl->delta[iPoint] = turbParam[5];
bl->ctEq[iPoint] = turbParam[6];
bl->us[iPoint] = turbParam[7];
}
else
{
std::cout << "Wrong regime\n"
<< std::endl;
throw;
}
// Gradients computation.
double dx = (bl->loc[iPoint] - bl->loc[iPoint-1]);
double dxExt = (bl->xxExt[iPoint] - bl->xxExt[iPoint-1]);
double dt_dx = (u[0] - bl->u[(iPoint - 1) * nVar + 0]) / dx;
double dH_dx = (u[1] - bl->u[(iPoint - 1) * nVar + 1]) / dx;
double due_dx = (u[3] - bl->u[(iPoint - 1) * nVar + 3]) / dx;
double dct_dx = (u[4] - bl->u[(iPoint - 1) * nVar + 4]) / dx;
double dHstar_dx = (bl->Hstar[iPoint] - bl->Hstar[iPoint - 1]) / dx;
double dueExt_dx = (bl->vt[iPoint] - bl->vt[iPoint-1]) / dxExt;
double ddeltaStarExt_dx = (bl->deltaStarExt[iPoint] - bl->deltaStarExt[iPoint-1]) / dxExt;
double c = 4 / (M_PI * dx);
double cExt = 4 / (M_PI * dxExt);
// Amplification routine.
double dN_dx = 0.0;
double ax = 0.0;
if (bl->xtrF == -1.0 && bl->regime[iPoint] == 0)
{
// No forced transition and laminar regime.
ax = amplificationRoutine(bl->Hk[iPoint], bl->Ret[iPoint], u[0]);
dN_dx = (u[2] - bl->u[(iPoint - 1) * nVar + 2]) / dx;
}
double Mea = bl->Me[iPoint];
double Hstara = bl->Hstar[iPoint];
double Hstar2a = bl->Hstar2[iPoint];
double Hka = bl->Hk[iPoint];
double cfa = bl->cf[iPoint]; double cda = bl->cd[iPoint];
double deltaa = bl->delta[iPoint];
double ctEqa = bl->ctEq[iPoint];
double usa = bl->us[iPoint];
//--------------------------------------------------------------------------------//
// Integral boundary layer equations. //
// 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(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);
// // 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(2, 2) = 1.;
// timeMatrix(3, 0) = -c * u[1];
// timeMatrix(3, 1) = -c * u[0];
// timeMatrix(3, 3) = 1.;
// 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);
// }
// else
// timeMatrix(4, 4) = 1.0;
//--------------------------------------------------------------------------------//
VectorXd result = VectorXd::Zero(5);
if (bl->regime[iPoint] == 0)
{
result[0] = 0.5*(-2.0*u[0]*(1.0*u[1] - 2.0)*(c*(dH_dx*u[0] + dt_dx*u[1]) - cExt*ddeltaStarExt_dx + dueExt_dx - due_dx) + 1.0*(c*u[0]*u[1] + u[3])*(2.0*due_dx*u[0]*(2.*Hstar2a - Hstara*(u[1] - 1)) + 1.0*u[3]*(Hstara*cfa - 4*cda + 2.*dHstar_dx*u[0])) + 1.0*(2.*due_dx*u[0]*(-Mea*Mea + u[1] + 2) + u[3]*(-cfa + 2.*dt_dx))*(1.0*Hstara*c*u[0]*u[1] + 1.0*Hstara*u[3] - 2.0*c*u[0] - 1.0*u[3]))/(c*u[0]*u[1] + u[3]);
result[1] = 0.5*(2.0*u[0]*u[1]*(u[1] - 1)*(c*(dH_dx*u[0] + dt_dx*u[1]) - cExt*ddeltaStarExt_dx + dueExt_dx - due_dx) - 1.0*u[1]*(c*u[0]*u[1] + u[3])*(2.*due_dx*u[0]*(2.*Hstar2a - Hstara*(u[1] - 1)) + u[3]*(Hstara*cfa - 4*cda + 2.*dHstar_dx*u[0])) + 1.0*(2.*due_dx*u[0]*(-Mea*Mea + u[1] + 2) + u[3]*(-cfa + 2.*dt_dx))*(-1.0*Hstara*c*u[0]*u[1]*u[1] - 1.0*Hstara*u[1]*u[3] + 2.0*c*u[0]*u[1] + 1.0*u[1]*u[3] + 1.0*u[3]))/(u[0]*(c*u[0]*u[1] + u[3]));
result[2] = -1.0*ax + 1.0*dN_dx;
result[3] = 1.0*u[3]*(-1.0*c*(dH_dx*u[0] + dt_dx*u[1]) + 0.5*c*(2.*due_dx*u[0]*(-Mea*Mea + u[1] + 2) + u[3]*(-cfa + 2.*dt_dx)) + 1.0*cExt*ddeltaStarExt_dx - 1.0*dueExt_dx + 1.0*due_dx)/(c*u[0]*u[1] + u[3]);
result[4] = 0.;
}
else if (bl->regime[iPoint] == 1)
{
result[0] = (-u[0]*(u[1] - 2)*(c*(dH_dx*u[0] + dt_dx*u[1]) - cExt*ddeltaStarExt_dx + dueExt_dx - due_dx) + (c*u[0]*u[1] + u[3])*(2.*due_dx*u[0]*(2.*Hstar2a - Hstara*(u[1] - 1)) + u[3]*(Hstara*cfa - 4*cda + 2.*dHstar_dx*u[0]))/2. + (2.*due_dx*u[0]*(-Mea*Mea + u[1] + 2) + u[3]*(-cfa + 2.*dt_dx))*(Hstara*c*u[0]*u[1] + Hstara*u[3] - 2.*c*u[0] - u[3])/2.)/(c*u[0]*u[1] + u[3]);
result[1] = (2.*u[0]*u[1]*(u[1] - 1)*(c*(dH_dx*u[0] + dt_dx*u[1]) - cExt*ddeltaStarExt_dx + dueExt_dx - due_dx) - u[1]*(c*u[0]*u[1] + u[3])*(2.*due_dx*u[0]*(2.*Hstar2a - Hstara*(u[1] - 1)) + u[3]*(Hstara*cfa - 4*cda + 2.*dHstar_dx*u[0])) + (2.*due_dx*u[0]*(-Mea*Mea + u[1] + 2) + u[3]*(-cfa + 2.*dt_dx))*(-Hstara*c*u[0]*u[1]*u[1] - Hstara*u[1]*u[3] + 2.*c*u[0]*u[1] + u[1]*u[3] + u[3]))/(2.*u[0]*(c*u[0]*u[1] + u[3]));
result[2] = -ax + dN_dx;
result[3] = u[3]*(-2.*c*(dH_dx*u[0] + dt_dx*u[1]) + c*(2.*due_dx*u[0]*(-Mea*Mea + u[1] + 2) + u[3]*(-cfa + 2.*dt_dx)) + 2.*cExt*ddeltaStarExt_dx - 2.*dueExt_dx + 2.*due_dx)/(2.*(c*u[0]*u[1] + u[3]));
result[4] = (-Hka*Hka*c*deltaa*dissipRatio*dissipRatio*u[0]*u[1]*u[4]*(2.*due_dx*u[0]*(-Mea*Mea + u[1] + 2) + u[3]*(-cfa + 2.*dt_dx))/2. + Hka*Hka*deltaa*dissipRatio*dissipRatio*u[0]*u[1]*u[4]*(c*(dH_dx*u[0] + dt_dx*u[1]) - cExt*ddeltaStarExt_dx + dueExt_dx - due_dx) + usa*(c*u[0]*u[1] + u[3])*(6*Hka*Hka*dct_dx*deltaa*dissipRatio*dissipRatio*u[0]*u[1]*u[3] + 16.8*Hka*Hka*dissipRatio*dissipRatio*u[0]*u[1]*u[3]*u[4]*(-ctEqa + dissipRatio*u[4]) + 2.*deltaa*u[4]*(3*Hka*Hka*dissipRatio*dissipRatio*due_dx*u[0]*u[1] + u[3]*(-2.*Hka*Hka*cfa*dissipRatio*dissipRatio + 0.0891067052795723*(Hka - 1)*(Hka - 1))))/6)/(Hka*Hka*deltaa*dissipRatio*dissipRatio*u[0]*u[1]*(c*u[0]*u[1] + u[3]));
}
return result;
}
/**
* @brief Amplification routines of the e^N method for transition capturing.
*
* @param Hk Kinematic shape parameter.
* @param Ret Reynolds number based on the momentum thickness.
* @param theta Momentum thickness.
* @return double */
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 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.;
if (Ret <= 0.)
Ret = 1.;
if (std::log10(Ret) < logRet_crit - dgr)
ax = 0.;
else
{
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;
// 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 arg = 3.87 * Hmi - 2.52;
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)
Hfac = 1.;
else
Hfac = 3 * Hnorm * Hnorm - 2 * Hnorm * Hnorm * Hnorm;
double ax1 = ax;
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;
if (ax2 < 0.)
ax2 = 0.;
ax = Hfac * ax2 + (1 - Hfac) * ax1;
if (ax < 0.)
ax = 0.;
}
}
return ax;
}