DDriver.cpp

#include "DDriver.h"
#include "DBoundaryLayer.h"
#include "DSolver.h"
#include <iomanip>

#define ANSI_COLOR_RED "\x1b[1;31m"
#define ANSI_COLOR_GREEN "\x1b[1;32m"
#define ANSI_COLOR_YELLOW "\x1b[1;33m"
#define ANSI_COLOR_RESET "\x1b[0m"

using namespace blast;

/**
 * @brief Driver of the viscous solver.
 *
 * @param _Re Freestream Reynolds number.
 * @param _Minf Freestream Mach number.
 * @param _CFL0 Initial CFL number for pseudo-time integration.
 * @param nSections Number of sections in the computation (1 for 2D cases).
 * @param _Span Span of the considered wing (only used for drag computation).
 * @param _verbose Verbosity level of the solver 0<=verbose<=1.
 */
Driver::Driver(double _Re, double _Minf, double _CFL0, size_t nSections,
               double xtrF_top, double xtrF_bot,  double _span, unsigned int _verbose)
{
    // Freestream parameters.
    Re = _Re;
    CFL0 = _CFL0;
    verbose = _verbose;

    std::vector<double> xtrF = {xtrF_top, xtrF_bot, 0.};

    // Initialzing boundary layer
    sections.resize(nSections);
    convergenceStatus.resize(nSections);
    for (size_t iSec = 0; iSec < nSections; ++iSec)
    {
        convergenceStatus[iSec].resize(3);
        for (size_t i = 0; i < 3; ++i)
        {
            BoundaryLayer *region = new BoundaryLayer(xtrF[i]);
            sections[iSec].push_back(region);
        }
        sections[iSec][0]->name = "upper";
        sections[iSec][1]->name = "lower";
        sections[iSec][2]->name = "wake";
    }

    Cdt_sec.resize(sections.size(), 0.);
    Cdf_sec.resize(sections.size(), 0.);
    Cdp_sec.resize(sections.size(), 0.);
    Cdt = 0.;
    Cdf = 0.;
    Cdp = 0.;
    span = _span;

    // Pre/post processing flags.
    stagPtMvmt.resize(sections.size(), false);

    // Time solver initialization.
    tSolver = new Solver(_CFL0, _Minf, Re);

    // Numerical parameters.
    std::cout << "--- Viscous solver setup ---\n"
              << "Reynolds number: " << Re << " \n"
              << "Freestream Mach number: " << _Minf << " \n"
              << "Initial CFL number: " << CFL0 << " \n"
              << std::endl;
    std::cout << "Boundary layer initialized." << std::endl;
}

/**
 * @brief Driver and subsequent dependencies destruction.
 */
Driver::~Driver()
{
    for (size_t iSec = 0; iSec < sections.size(); ++iSec)
        for (size_t i = 0; i < sections[iSec].size(); ++i)
            delete sections[iSec][i];
    delete tSolver;
    std::cout << "~blast::Driver()\n";
} // end ~Driver

/**
 * @brief Runs viscous operations. Temporal solver parametrisation is first adapted,
 * Solutions are initialized than time marched towards steady state successively for each points.
 *
 * @param couplIter Current coupling iteration.
 * @return int Solver exit code.
 */
int Driver::run(unsigned int const couplIter)
{
    bool lockTrans = false; // Flagging activation of transition routines.
    int pointExitCode;      // Output of pseudo time integration (0: converged; -1: unsuccessful, failed; >0: unsuccessful nb iter).

    double maxMach = 0.;

    for (size_t iSec = 0; iSec < sections.size(); ++iSec)
    {
        tms["0-CheckIn"].start();
        if (verbose > 1)
            std::cout << "Sec " << iSec << " ";

        // Check for stagnation point movement.
        if (couplIter != 0 && (sections[iSec][0]->mesh->getDiffnElm()))
        {
            stagPtMvmt[iSec] = true;
            if (verbose > 2)
                std::cout << "Stagnation point moved by " << sections[iSec][0]->mesh->getDiffnElm() << " elements." << std::endl;
        }
        else
            stagPtMvmt[iSec] = false;

        // Set boundary conditions.
        sections[iSec][0]->xtr = 1.; // Upper side initial xtr.
        sections[iSec][1]->xtr = 1.; // Lower side initial xtr.
        sections[iSec][2]->xtr = 0.; // Wake initial xtr (fully turbulent).
        sections[iSec][0]->setStagBC(Re);
        sections[iSec][1]->setStagBC(Re);
        tms["0-CheckIn"].stop();

        // Loop over the regions (upper, lower and wake).
        for (size_t iRegion = 0; iRegion < sections[iSec].size(); ++iRegion)
        {
            convergenceStatus[iSec][iRegion].resize(0);

            // Maximum Mach number in the region.
            if (sections[iSec][iRegion]->blEdge->getMaxMach() > maxMach)
                maxMach = sections[iSec][iRegion]->blEdge->getMaxMach();

            if (verbose > 1)
                std::cout << sections[iSec][iRegion]->name << " ";

            // Check if safe mode is required.
            if (couplIter == 0 || sections[iSec][iRegion]->mesh->getDiffnElm() != 0 || stagPtMvmt[iSec] == true || solverExitCode != 0)
            {
                tSolver->setCFL0(1);
                tSolver->setitFrozenJac(1);
                tSolver->setinitSln(true);

                if (sections[iSec][iRegion]->mesh->getDiffnElm())
                {
                    std::vector<double> xxExtCurr(sections[iSec][iRegion]->mesh->getnMarkers(), 0.);
                    for (size_t iPoint = 0; iPoint < xxExtCurr.size(); ++iPoint)
                        xxExtCurr[iPoint] = sections[iSec][iRegion]->mesh->getLoc(iPoint);
                    sections[iSec][iRegion]->mesh->setExt(xxExtCurr);

                    std::vector<double> DeltaStarZeros(sections[iSec][iRegion]->mesh->getnMarkers(), 0.0);
                    sections[iSec][iRegion]->blEdge->setDeltaStar(DeltaStarZeros);

                    std::vector<double> ueOld(sections[iSec][iRegion]->mesh->getnMarkers(), 0.0);
                    for (size_t iPoint = 0; iPoint < ueOld.size(); ++iPoint)
                        ueOld[iPoint] = sections[iSec][iRegion]->blEdge->getVt(iPoint);
                    sections[iSec][iRegion]->blEdge->setVtOld(ueOld);
                }

                // Keyword annoncing iteration safety level.
                if (verbose > 1)
                    std::cout << "restart. ";
            }
            else
            {
                tSolver->setCFL0(CFL0);
                tSolver->setitFrozenJac(1);
                tSolver->setinitSln(true);

                // Keyword annoncing iteration safety level.
                if (verbose > 1)
                    std::cout << "update. ";
            }

            // Save number of points in each regions.
            sections[iSec][iRegion]->mesh->prevnMarkers = sections[iSec][iRegion]->mesh->getnMarkers();

            // Reset regime.
            lockTrans = false;
            if (iRegion < 2) // Airfoil
                for (size_t k = 0; k < sections[iSec][iRegion]->mesh->getnMarkers(); k++)
                    sections[iSec][iRegion]->regime[k] = 0;
            else if (iRegion == 2) // Wake
            {
                for (size_t k = 0; k < sections[iSec][iRegion]->mesh->getnMarkers(); k++)
                    sections[iSec][iRegion]->regime[k] = 1;

                // Impose wake boundary condition.
                std::vector<double> upperConditions(sections[iSec][iRegion]->getnVar(), 0.);
                std::vector<double> lowerConditions(sections[iSec][iRegion]->getnVar(), 0.);
                for (size_t k = 0; k < upperConditions.size(); ++k)
                {
                    upperConditions[k] = sections[iSec][0]->u[(sections[iSec][0]->mesh->getnMarkers() - 1) * sections[iSec][0]->getnVar() + k];
                    lowerConditions[k] = sections[iSec][1]->u[(sections[iSec][1]->mesh->getnMarkers() - 1) * sections[iSec][1]->getnVar() + k];
                }
                sections[iSec][iRegion]->setWakeBC(Re, upperConditions, lowerConditions);
                lockTrans = true;
            }
            else
                throw;

            // Loop over points
            for (size_t iPoint = 1; iPoint < sections[iSec][iRegion]->mesh->getnMarkers(); ++iPoint)
            {
                // Initialize solution at point
                tSolver->initialCondition(iPoint, sections[iSec][iRegion]);
                // Solve equations
                tms["1-Solver"].start();
                pointExitCode = tSolver->integration(iPoint, sections[iSec][iRegion]);
                if (sections[iSec][iRegion]->xtrF != -1 && sections[iSec][iRegion]->mesh->getxoc(iPoint) >= sections[iSec][iRegion]->xtrF)
                    sections[iSec][iRegion]->u[iPoint * sections[iSec][iRegion]->getnVar()+2] = 9.;
                tms["1-Solver"].stop();

                if (pointExitCode != 0)
                    convergenceStatus[iSec][iRegion].push_back(iPoint);

                // Transition
                tms["2-Transition"].start();
                if (!lockTrans)
                {
                    // Check if perturbation waves are growing or not.
                    if (sections[iSec][iRegion]->xtrF != -1 && sections[iSec][iRegion]->mesh->getxoc(iPoint) <= sections[iSec][iRegion]->xtrF)
                        checkWaves(iPoint, sections[iSec][iRegion]);

                    // Amplification factor is compared to critical amplification factor.
                    if (sections[iSec][iRegion]->u[iPoint * sections[iSec][iRegion]->getnVar() + 2] >= sections[iSec][iRegion]->getnCrit())
                    {
                        averageTransition(iPoint, sections[iSec][iRegion], 0);
                        if (verbose > 1)
                        {
                            std::cout << std::fixed;
                            std::cout << std::setprecision(2);
                            std::cout << sections[iSec][iRegion]->xoctr
                                      << " (" << sections[iSec][iRegion]->transMarker << ") ";
                        }
                        lockTrans = true;
                    }
                }
                tms["2-Transition"].stop();
            }
            tms["3-Blowing"].start();
            sections[iSec][iRegion]->computeBlwVel();
            tms["3-Blowing"].stop();
        }
        if (verbose > 1)
            std::cout << std::endl;
    }

    for (size_t i = 0; i < sections.size(); ++i)
        computeSectionalDrag(sections[i], i);
    computeTotalDrag();

    // Output information.
    solverExitCode = outputStatus(maxMach);

    return solverExitCode; // exit code
}

/**
 * @brief Check whether or not instabilities are growing in the laminar boundary layer.
 *
 * @param iPoint Marker id.
 * @param bl BoundaryLayer region.
 */
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);

    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.
    }
    else
        bl->u[iPoint * bl->getnVar() + 2] = 0.; // Set N to 0 at that point if waves do not grow.
}

/**
 * @brief Computes turbulent solution @ the transition point and averages solutions.
 *
 * @param iPoint Marker id.
 * @param bl BoundaryLayer region.
 * @param forced Flag for forced transition.
 */
void Driver::averageTransition(size_t iPoint, BoundaryLayer *bl, int forced)
{
    // Averages solution on transition marker.
    size_t nVar = bl->getnVar();

    // Save laminar solution.
    std::vector<double> lamSol(nVar, 0.0);
    for (size_t k = 0; k < lamSol.size(); ++k)
        lamSol[k] = bl->u[iPoint * nVar + k];

    // Set up turbulent portion boundary condition.
    bl->transMarker = iPoint; // Save transition marker.

    // Loop over remaining points in the region.
    for (size_t k = iPoint; k < bl->mesh->getnMarkers(); ++k)
        bl->regime[k] = 1; // Set Turbulent regime.

    // Compute transition location.
    bl->xtr = (bl->getnCrit() - bl->u[(iPoint - 1) * nVar + 2]) * ((bl->mesh->getx(iPoint) - bl->mesh->getx(iPoint - 1)) / (bl->u[iPoint * nVar + 2] - bl->u[(iPoint - 1) * nVar + 2])) + bl->mesh->getx(iPoint - 1);
    bl->xoctr = (bl->getnCrit() - bl->u[(iPoint - 1) * nVar + 2]) * ((bl->mesh->getxoc(iPoint) - bl->mesh->getxoc(iPoint - 1)) / (bl->u[iPoint * nVar + 2] - bl->u[(iPoint - 1) * nVar + 2])) + bl->mesh->getxoc(iPoint - 1);

    // Percentage of the subsection corresponding to a turbulent flow.
    double avTurb = (bl->mesh->getx(iPoint) - bl->xtr) / (bl->mesh->getx(iPoint) - bl->mesh->getx(iPoint - 1));
    double avLam = 1. - avTurb;

    // Impose boundary condition.
    double Cteq_trans;
    bl->closSolver->turbulentClosures(Cteq_trans, bl->u[(iPoint - 1) * nVar + 0], bl->u[(iPoint - 1) * nVar + 1], 0.03, bl->Ret[iPoint - 1], bl->blEdge->getMe(iPoint - 1), bl->name);
    bl->ctEq[iPoint - 1] = Cteq_trans;
    bl->u[(iPoint - 1) * nVar + 4] = 0.7 * bl->ctEq[iPoint - 1];

    // Avoid starting with ill conditioned IC for turbulent computation. (Regime was switched above).
    // These initial guess values do not influence the converged solution.
    bl->u[iPoint * nVar + 4] = bl->u[(iPoint - 1) * nVar + 4]; // IC of transition point = turbulent BC (imposed @ previous point).
    bl->u[iPoint * nVar + 1] = 1.515;                          // Because we expect a significant drop of the shape factor.

    // Solve point in turbulent configuration.
    int exitCode = tSolver->integration(iPoint, bl);

    if (exitCode != 0 && verbose > 1)
        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+2] = 9.;

    // Recompute closures. (The turbulent BC @ iPoint - 1 does not influence laminar closure relations).
    std::vector<double> lamParam(8, 0.);
    bl->closSolver->laminarClosures(lamParam, bl->u[(iPoint - 1) * nVar + 0], bl->u[(iPoint - 1) * nVar + 1], bl->Ret[iPoint - 1], bl->blEdge->getMe(iPoint - 1), bl->name);
    bl->Hstar[iPoint - 1] = lamParam[0];
    bl->Hstar2[iPoint - 1] = lamParam[1];
    bl->Hk[iPoint - 1] = lamParam[2];
    bl->cf[iPoint - 1] = lamParam[3];
    bl->cd[iPoint - 1] = lamParam[4];
    bl->delta[iPoint - 1] = lamParam[5];
    bl->ctEq[iPoint - 1] = lamParam[6];
    bl->us[iPoint - 1] = lamParam[7];

    std::vector<double> turbParam(8, 0.);
    bl->closSolver->turbulentClosures(turbParam, bl->u[(iPoint)*nVar + 0], bl->u[(iPoint)*nVar + 1], bl->u[(iPoint)*nVar + 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];
}

/**
 * @brief Friction, pressure and total drag computation qt each station of the configuration.
 *
 * @tparam Cdt = theta * Ue^((H+5)/2).
 * @tparam Cdf = integral(Cf dx)_upper + integral(Cf dx)_lower.
 * @tparam Cdp = Cdtot - Cdf.
 *
 * @param bl BoundaryLayer region.
 * @param i Considered section.
 */
void Driver::computeSectionalDrag(std::vector<BoundaryLayer *> const &bl, size_t i)
{
    size_t nVar = bl[0]->getnVar();
    size_t lastWkPt = (bl[2]->mesh->getnMarkers() - 1) * nVar;

    // Normalize friction and dissipation coefficients.
    std::vector<double> frictioncoeff_upper(bl[0]->mesh->getnMarkers(), 0.0);
    for (size_t k = 0; k < frictioncoeff_upper.size(); ++k)
        frictioncoeff_upper[k] = bl[0]->blEdge->getVt(k) * bl[0]->blEdge->getVt(k) * bl[0]->cf[k];

    std::vector<double> frictioncoeff_lower(bl[1]->mesh->getnMarkers(), 0.0);
    for (size_t k = 0; k < frictioncoeff_lower.size(); ++k)
        frictioncoeff_lower[k] = bl[1]->blEdge->getVt(k) * bl[1]->blEdge->getVt(k) * bl[1]->cf[k];

    // Compute friction drag coefficient.
    double Cdf_upper = 0.0;
    for (size_t k = 1; k < bl[0]->mesh->getnMarkers(); ++k)
        Cdf_upper += (frictioncoeff_upper[k] + frictioncoeff_upper[k - 1]) * (bl[0]->mesh->getx(k) - bl[0]->mesh->getx(k - 1)) / 2;
    double Cdf_lower = 0.0;
    for (size_t k = 1; k < bl[1]->mesh->getnMarkers(); ++k)
        Cdf_lower += (frictioncoeff_lower[k] + frictioncoeff_lower[k - 1]) * (bl[1]->mesh->getx(k) - bl[1]->mesh->getx(k - 1)) / 2;

    // Friction drag coefficient.
    Cdf_sec[i] = Cdf_upper + Cdf_lower; // No friction in the wake
    // Total drag coefficient.
    Cdt_sec[i] = 2 * bl[2]->u[lastWkPt + 0] * pow(bl[2]->u[lastWkPt + 3], (bl[2]->u[lastWkPt + 1] + 5) / 2);
    // Pressure drag coefficient.
    Cdp_sec[i] = Cdt_sec[i] - Cdf_sec[i];
}

double Driver::getAvrgxoctr(size_t const iRegion) const
{
    double xtr = 0.0;
    for (size_t iSec = 0; iSec < sections.size(); ++iSec)
        xtr += sections[iSec][iRegion]->xoctr;
    xtr /= sections.size();
    return xtr;
}

/**
 * @brief Compute total drag coefficient on the airfoil/wing.
 * Performs the sectional drag integration for 3D computations.
 *
 */
void Driver::computeTotalDrag()
{
    Cdt = 0.;
    Cdf = 0.;
    Cdp = 0.;

    if (sections.size() == 1 || span == 0.)
    {
        Cdt = Cdt_sec[0];
        Cdf = Cdf_sec[0];
        Cdp = Cdp_sec[0];
    }
    else
    {
        Cdt += (sections[0][0]->mesh->getz(0) - 0) * Cdt_sec[0];
        Cdp += (sections[0][0]->mesh->getz(0) - 0) * Cdp_sec[0];
        Cdf += (sections[0][0]->mesh->getz(0) - 0) * Cdf_sec[0];
        double dz = 0.;
        for (size_t k = 1; k < sections.size(); ++k)
        {
            dz = (sections[k][0]->mesh->getz(0) - sections[k - 1][0]->mesh->getz(0));
            Cdt += dz * Cdt_sec[k];
            Cdp += dz * Cdp_sec[k];
            Cdf += dz * Cdf_sec[k];
        }
        Cdt /= span;
        Cdp /= span;
        Cdf /= span;
    }
}

/**
 * @brief Outputs solver status depending on the verbosity level (verbose).
 *
 * @param maxMach Maximum Mach number identified on the configuration.
 * @return int Solver exit code (0 converged, !=0 not converged).
 */
int Driver::outputStatus(double maxMach) const
{
    int solverExitCode = 0;
    for (size_t iSec = 0; iSec < sections.size(); ++iSec)
        for (size_t iReg = 0; iReg < sections[iSec].size(); ++iReg)
            if (convergenceStatus[iSec][iReg].size() != 0)
            {
                solverExitCode = -1;
                break;
            }

    if (verbose > 2 && solverExitCode != 0)
    {
        // Output unconverged points.
        std::cout << "Unconverged points" << std::endl;
        std::cout << std::setw(10) << "Section"
                  << std::setw(12) << "Region"
                  << std::setw(16) << "Markers" << std::endl;
        for (size_t iSec = 0; iSec < sections.size(); ++iSec)
        {
            unsigned int count = 0;
            for (size_t iReg = 0; iReg < sections[iSec].size(); ++iReg)
            {
                std::cout << std::setw(10) << iSec;
                std::cout << std::setw(12) << iReg << "          ";
                if (convergenceStatus[iSec][iReg].size() != 0)
                    for (size_t iPoint = 0; iPoint < convergenceStatus[iSec][iReg].size(); ++iPoint)
                    {
                        std::cout << convergenceStatus[iSec][iReg][iPoint] << " (" << sections[iSec][iReg]->mesh->getx(convergenceStatus[iSec][iReg][iPoint]) << "), ";
                        ++count;
                    }
                std::cout << "\n";
            }
            std::cout << "total: " << count << "\n";
        }
        std::cout << "" << std::endl;
    }

    if (verbose > 0)
    {
        // Compute mean transition locations.
        std::vector<double> meanTransitions(2, 0.);
        for (size_t iSec = 0; iSec < sections.size(); ++iSec)
            for (size_t iReg = 0; iReg < 2; ++iReg)
                meanTransitions[iReg] += sections[iSec][iReg]->xoctr;
        for (size_t iReg = 0; iReg < meanTransitions.size(); ++iReg)
            meanTransitions[iReg] /= sections.size();

        std::cout << std::fixed << std::setprecision(2)
                  << "Viscous solver for Re " << Re / 1e6 << "e6, Mach max "
                  << maxMach
                  << std::endl;

        std::cout << std::setw(10) << "xTr Upper"
                  << std::setw(12) << "xTr Lower"
                  << std::setw(12) << "Cd" << std::endl;
        std::cout << std::fixed << std::setprecision(2);
        std::cout << std::setw(10) << meanTransitions[0]
                  << std::setw(12) << meanTransitions[1];
        std::cout << std::fixed << std::setprecision(5);
        std::cout << std::setw(12) << Cdt << "\n";

        if (verbose > 0)
        {
            if (solverExitCode == 0)
                std::cout << ANSI_COLOR_GREEN << "Viscous solver converged." << ANSI_COLOR_RESET << std::endl;
            else
                std::cout << ANSI_COLOR_YELLOW << "Viscous solver not converged." << ANSI_COLOR_RESET << std::endl;
        }

        if (verbose > 2)
        {
            std::cout << "Viscous solver CPU" << std::endl
                      << tms;
        }
    }

    if (verbose > 0)
        std::cout << "" << std::endl;
    return solverExitCode;
}

/**
 * @brief Set nodes coordinates in the global frame of reference.
 *
 * @param iSec id of the considered section.
 * @param iRegion id of the considered region.
 * @param _x Vector of x coordinates of the nodes in the region.
 * @param _y Vector of y coordinates of the nodes in the region.
 * @param _z Vector of z coordinates of the nodes in the region.
 */
void Driver::setGlobMesh(size_t iSec, size_t iRegion, std::vector<double> _x, std::vector<double> _y, std::vector<double> _z)
{
    if (iSec > sections.size() || iRegion > sections[iSec].size())
        throw std::runtime_error("blast::Driver Wrong section or region id.");
    sections[iSec][iRegion]->setMesh(_x, _y, _z);
}

/**
 * @brief Set the velocity at the nodes in the global frame of reference.
 *
 * @param iSec id of the considered section.
 * @param iRegion id of the considered region.
 * @param vx Vector of velocities projected on the x axis of the global frame of reference.
 * @param vy Vector of velocities projected on the y axis of the global frame of reference.
 * @param vz Vector of velocities projected on the z axis of the global frame of reference.
 */
void Driver::setVelocities(size_t iSec, size_t iRegion, std::vector<double> vx, std::vector<double> vy, std::vector<double> vz)
{
    if (vx.size() != vy.size() || vy.size() != vz.size() || vx.size() != vz.size())
        throw std::runtime_error("blast::Driver Wrong velocity vector sizes.");
    if (vx.size() != sections[iSec][iRegion]->mesh->getnMarkers())
        throw std::runtime_error("blast::Driver Velocity vector is not coherent with the mesh.");

    std::vector<double> vt(vx.size(), 0.0);

    for (size_t iPoint = 0; iPoint < vx.size() - 1; ++iPoint)
    {
        vt[iPoint] = vx[iPoint] * std::cos(sections[iSec][iRegion]->mesh->getAlpha(iPoint)) + vy[iPoint] * std::sin(sections[iSec][iRegion]->mesh->getAlpha(iPoint));
        vt[iPoint + 1] = vx[iPoint + 1] * std::cos(sections[iSec][iRegion]->mesh->getAlpha(iPoint)) + vy[iPoint + 1] * std::sin(sections[iSec][iRegion]->mesh->getAlpha(iPoint));
    }

    sections[iSec][iRegion]->blEdge->setVt(vt);
}

/**
 * @brief Set the Mach number at the nodes.
 *
 * @param iSec id of the considered section.
 * @param iRegion id of the considered region.
 * @param _Me Vector of Mach numbers at the nodes.
 */
void Driver::setMe(size_t iSec, size_t iRegion, std::vector<double> _Me)
{
    if (_Me.size() != sections[iSec][iRegion]->mesh->getnMarkers())
        throw std::runtime_error("blast::Driver Mach number vector is not coherent with the mesh.");
    sections[iSec][iRegion]->blEdge->setMe(_Me);
}

/**
 * @brief Set the density at the nodes.
 *
 * @param iSec id of the considered section.
 * @param iRegion id of the considered region.
 * @param _Rhoe Vector of density at the nodes.
 */
void Driver::setRhoe(size_t iSec, size_t iRegion, std::vector<double> _Rhoe)
{
    if (_Rhoe.size() != sections[iSec][iRegion]->mesh->getnMarkers())
        throw std::runtime_error("blast::Driver Density vector is not coherent with the mesh.");
    sections[iSec][iRegion]->blEdge->setRhoe(_Rhoe);
}

/**
 * @brief Set the nodes coordinates in the local frame of reference at the edge of the boundary layer.
 *
 * @param iSec id of the considered section.
 * @param iRegion id of the considered region.
 * @param _xxExt Vector of coordinates of the nodes.
 */
void Driver::setxxExt(size_t iSec, size_t iRegion, std::vector<double> _xxExt)
{
    if (_xxExt.size() != sections[iSec][iRegion]->mesh->getnMarkers())
        throw std::runtime_error("blast::Driver External mesh vector is not coherent with the mesh.");
    sections[iSec][iRegion]->mesh->setExt(_xxExt);
}

/**
 * @brief Set the displacement thickness (previous iteration).
 *
 * @param iSec id of the considered section.
 * @param iRegion id of the considered region.
 * @param _DeltaStarExt Vector of displacement thicknesses at the nodes.
 */
void Driver::setDeltaStarExt(size_t iSec, size_t iRegion, std::vector<double> _DeltaStarExt)
{
    if (_DeltaStarExt.size() != sections[iSec][iRegion]->mesh->getnMarkers())
    {
        std::cout << "size DeltaStar: " << _DeltaStarExt.size() << ". nMarkers: " << sections[iSec][iRegion]->mesh->getnMarkers() << std::endl;
        throw std::runtime_error("blast::Driver External delta star vector is not coherent with the mesh.");
    }
    sections[iSec][iRegion]->blEdge->setDeltaStar(_DeltaStarExt);
}

/**
 * @brief Returns x coordinates of the nodes in the local frame of reference in the considered region.
 *
 * @param iSec id of the considered section.
 * @param iRegion id of the considered region.
 * @return std::vector<double>
 */
std::vector<double> Driver::extractxCoord(size_t iSec, size_t iRegion) const
{
    std::vector<double> xCoord(sections[iSec][iRegion]->mesh->getnMarkers(), 0.);
    for (size_t iPoint = 0; iPoint < xCoord.size(); ++iPoint)
        xCoord[iPoint] = sections[iSec][iRegion]->mesh->getx(iPoint);
    return xCoord;
}

/**
 * @brief Returns blowing velocities (defined on the cells' cg) in the considered region.
 *
 * @param iSec id of the considered section.
 * @param iRegion id of the considered region.
 * @return std::vector<double>
 */
std::vector<double> Driver::extractBlowingVel(size_t iSec, size_t iRegion) const
{
    return sections[iSec][iRegion]->blowingVelocity;
}
std::vector<double> Driver::extractxx(size_t iSec, size_t iRegion) const
{
    std::vector<double> xx(sections[iSec][iRegion]->mesh->getnMarkers(), 0.);
    for (size_t iPoint = 0; iPoint < xx.size(); ++iPoint)
        xx[iPoint] = sections[iSec][iRegion]->mesh->getLoc(iPoint);
    return xx;
}

/**
 * @brief Returns the displacement thickness in the considered region.
 *
 * @param iSec id of the considered section.
 * @param iRegion id of the considered region.
 * @return std::vector<double>
 */
std::vector<double> Driver::extractDeltaStar(size_t iSec, size_t iRegion) const
{
    return sections[iSec][iRegion]->deltaStar;
}

/**
 * @brief Returns the velocity of the interaction law in the considered region.
 *
 * @param iSec id of the considered section.
 * @param iRegion id of the considered region.
 * @return std::vector<double>
 */
std::vector<double> Driver::extractUe(size_t iSec, size_t iRegion) const
{
    std::vector<double> ueCurr(sections[iSec][iRegion]->mesh->getnMarkers(), 0.0);
    for (size_t i = 0; i < ueCurr.size(); ++i)
        ueCurr[i] = sections[iSec][iRegion]->u[i * sections[iSec][iRegion]->getnVar() + 3];
    return ueCurr;
}

/**
 * @brief Set the velocity at the edge of the BL at the previous iteration in the considered region.
 *
 * @param iSec id of the considered section.
 * @param iRegion id of the considered region.
 * @param _ue Velocity vector.
 * @return std::vector<double>
 */
void Driver::setUeOld(size_t iSec, size_t iRegion, std::vector<double> _ue)
{
    if (_ue.size() != sections[iSec][iRegion]->mesh->getnMarkers())
        throw std::runtime_error("blast::Driver External velocity vector is not coherent with the mesh.");
    sections[iSec][iRegion]->blEdge->setVtOld(_ue);
}

/**
 * @brief Returns the integral boundary layer solution in a section (sorted from upper Te to lower Te than wake).
 *
 * @param iSec id of the considered section.
 * @return std::vector<std::vector<double>>
 */
std::vector<std::vector<double>> Driver::getSolution(size_t iSec)
{
    size_t nMarkersUp = sections[iSec][0]->mesh->getnMarkers();
    size_t nMarkersLw = sections[iSec][1]->mesh->getnMarkers(); // No stagnation point doublon.
    size_t nVar = sections[iSec][0]->getnVar();

    std::vector<std::vector<double>> Sln(20);
    for (size_t iVector = 0; iVector < Sln.size(); ++iVector)
        Sln[iVector].resize(nMarkersUp + nMarkersLw + sections[iSec][2]->mesh->getnMarkers(), 0.);
    // Upper side.
    size_t revPoint = nMarkersUp - 1;
    for (size_t iPoint = 0; iPoint < nMarkersUp; ++iPoint)
    {
        Sln[0][iPoint] = sections[iSec][0]->u[revPoint * nVar + 0];
        Sln[1][iPoint] = sections[iSec][0]->u[revPoint * nVar + 1];
        Sln[2][iPoint] = sections[iSec][0]->u[revPoint * nVar + 2];
        Sln[3][iPoint] = sections[iSec][0]->u[revPoint * nVar + 3];
        Sln[4][iPoint] = sections[iSec][0]->u[revPoint * nVar + 4];
        Sln[5][iPoint] = sections[iSec][0]->deltaStar[revPoint];

        Sln[6][iPoint] = sections[iSec][0]->Ret[revPoint];
        Sln[7][iPoint] = sections[iSec][0]->cd[revPoint];
        Sln[8][iPoint] = sections[iSec][0]->cf[revPoint];
        Sln[9][iPoint] = sections[iSec][0]->Hstar[revPoint];
        Sln[10][iPoint] = sections[iSec][0]->Hstar2[revPoint];
        Sln[11][iPoint] = sections[iSec][0]->Hk[revPoint];
        Sln[12][iPoint] = sections[iSec][0]->ctEq[revPoint];
        Sln[13][iPoint] = sections[iSec][0]->us[revPoint];
        Sln[14][iPoint] = sections[iSec][0]->delta[revPoint];

        Sln[15][iPoint] = sections[iSec][0]->mesh->getx(revPoint);
        Sln[16][iPoint] = sections[iSec][0]->mesh->getxoc(revPoint);
        Sln[17][iPoint] = sections[iSec][0]->mesh->gety(revPoint);
        Sln[18][iPoint] = sections[iSec][0]->mesh->getz(revPoint);
        Sln[19][iPoint] = sections[iSec][0]->blEdge->getVt(revPoint);
        --revPoint;
    }
    // Lower side.
    for (size_t iPoint = 0; iPoint < sections[iSec][1]->mesh->getnMarkers(); ++iPoint)
    {
        Sln[0][nMarkersUp + iPoint] = sections[iSec][1]->u[iPoint * nVar + 0];
        Sln[1][nMarkersUp + iPoint] = sections[iSec][1]->u[iPoint * nVar + 1];
        Sln[2][nMarkersUp + iPoint] = sections[iSec][1]->u[iPoint * nVar + 2];
        Sln[3][nMarkersUp + iPoint] = sections[iSec][1]->u[iPoint * nVar + 3];
        Sln[4][nMarkersUp + iPoint] = sections[iSec][1]->u[iPoint * nVar + 4];
        Sln[5][nMarkersUp + iPoint] = sections[iSec][1]->deltaStar[iPoint];

        Sln[6][nMarkersUp + iPoint] = sections[iSec][1]->Ret[iPoint];
        Sln[7][nMarkersUp + iPoint] = sections[iSec][1]->cd[iPoint];
        Sln[8][nMarkersUp + iPoint] = sections[iSec][1]->cf[iPoint];
        Sln[9][nMarkersUp + iPoint] = sections[iSec][1]->Hstar[iPoint];
        Sln[10][nMarkersUp + iPoint] = sections[iSec][1]->Hstar2[iPoint];
        Sln[11][nMarkersUp + iPoint] = sections[iSec][1]->Hk[iPoint];
        Sln[12][nMarkersUp + iPoint] = sections[iSec][1]->ctEq[iPoint];
        Sln[13][nMarkersUp + iPoint] = sections[iSec][1]->us[iPoint];
        Sln[14][nMarkersUp + iPoint] = sections[iSec][1]->delta[iPoint];

        Sln[15][nMarkersUp + iPoint] = sections[iSec][1]->mesh->getx(iPoint);
        Sln[16][nMarkersUp + iPoint] = sections[iSec][1]->mesh->getxoc(iPoint);
        Sln[17][nMarkersUp + iPoint] = sections[iSec][1]->mesh->gety(iPoint);
        Sln[18][nMarkersUp + iPoint] = sections[iSec][1]->mesh->getz(iPoint);
        Sln[19][nMarkersUp + iPoint] = sections[iSec][1]->blEdge->getVt(iPoint);
    }
    // Wake.
    for (size_t iPoint = 0; iPoint < sections[iSec][2]->mesh->getnMarkers(); ++iPoint)
    {
        Sln[0][nMarkersUp + nMarkersLw + iPoint] = sections[iSec][2]->u[iPoint * nVar + 0];
        Sln[1][nMarkersUp + nMarkersLw + iPoint] = sections[iSec][2]->u[iPoint * nVar + 1];
        Sln[2][nMarkersUp + nMarkersLw + iPoint] = sections[iSec][2]->u[iPoint * nVar + 2];
        Sln[3][nMarkersUp + nMarkersLw + iPoint] = sections[iSec][2]->u[iPoint * nVar + 3];
        Sln[4][nMarkersUp + nMarkersLw + iPoint] = sections[iSec][2]->u[iPoint * nVar + 4];
        Sln[5][nMarkersUp + nMarkersLw + iPoint] = sections[iSec][2]->deltaStar[iPoint];

        Sln[6][nMarkersUp + nMarkersLw + iPoint] = sections[iSec][2]->Ret[iPoint];
        Sln[7][nMarkersUp + nMarkersLw + iPoint] = sections[iSec][2]->cd[iPoint];
        Sln[8][nMarkersUp + nMarkersLw + iPoint] = sections[iSec][2]->cf[iPoint];
        Sln[9][nMarkersUp + nMarkersLw + iPoint] = sections[iSec][2]->Hstar[iPoint];
        Sln[10][nMarkersUp + nMarkersLw + iPoint] = sections[iSec][2]->Hstar2[iPoint];
        Sln[11][nMarkersUp + nMarkersLw + iPoint] = sections[iSec][2]->Hk[iPoint];
        Sln[12][nMarkersUp + nMarkersLw + iPoint] = sections[iSec][2]->ctEq[iPoint];
        Sln[13][nMarkersUp + nMarkersLw + iPoint] = sections[iSec][2]->us[iPoint];
        Sln[14][nMarkersUp + nMarkersLw + iPoint] = sections[iSec][2]->delta[iPoint];

        Sln[15][nMarkersUp + nMarkersLw + iPoint] = sections[iSec][2]->mesh->getx(iPoint);
        Sln[16][nMarkersUp + nMarkersLw + iPoint] = sections[iSec][2]->mesh->getxoc(iPoint);
        Sln[17][nMarkersUp + nMarkersLw + iPoint] = sections[iSec][2]->mesh->gety(iPoint);
        Sln[18][nMarkersUp + nMarkersLw + iPoint] = sections[iSec][2]->mesh->getz(iPoint);
        Sln[19][nMarkersUp + nMarkersLw + iPoint] = sections[iSec][2]->blEdge->getVt(iPoint);
    }

    if (verbose > 2)
        std::cout << "Solution structure sent." << std::endl;
    return Sln;
}