DFluxes.cpp

#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();

    MatrixXd timeMatrix = MatrixXd::Zero(5, 5);
    VectorXd spaceVector = VectorXd::Zero(5);

    // Dissipation ratio.
    double dissipRatio;
    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->blEdge->getMe(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->blEdge->getMe(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->mesh->getLoc(iPoint) - bl->mesh->getLoc(iPoint - 1));
    double dxExt = (bl->mesh->getExt(iPoint) - bl->mesh->getExt(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->blEdge->getVt(iPoint) - bl->blEdge->getVt(iPoint - 1)) / dxExt;
    double ddeltaStarExt_dx = (bl->blEdge->getDeltaStar(iPoint) - bl->blEdge->getDeltaStar(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->blEdge->getMe(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];

    // 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;

    return timeMatrix.partialPivLu().solve(spaceVector);
}

/**
 * @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;
}