vii/src/DDriver.h

-#ifndef DDRIVER_H
-#define DDRIVER_H
-
-#include "vii.h"
-#include "wObject.h"
-#include "DBoundaryLayer.h"
-#include "wTimers.h"
-namespace vii
-{
-
-/**
- * @brief Boundary layer solver driver. Performs the space marching and set up time marching for each point.
- * @authors Paul Dechamps
- */
-class VII_API Driver : public fwk::wSharedObject
-{
-public:
-    double Cdt;                                         ///< Total drag coefficient.
-    double Cdf;                                         ///< Total friction drag coefficient.
-    double Cdp;                                         ///< Total pressure drag coefficient.
-    std::vector<double> Cdt_sec;                        ///< Sectional total drag coefficient.
-    std::vector<double> Cdf_sec;                        ///< Sectional friction drag coefficient.
-    std::vector<double> Cdp_sec;                        ///< Sectional pressure drag coefficient.
-    std::vector<std::vector<BoundaryLayer *>> sections; ///< Boundary layer regions sorted by section.
-    unsigned int verbose;                               ///< Verbosity level of the solver 0<=verbose<=1.
-
-private:
-    double Re;                                                       ///< Freestream Reynolds number.
-    double CFL0;                                                     ///< Initial CFL number for pseudo time integration.
-    double span;                                                     ///< Wing Span (Used only for drag computation, not used if 2D case).
-    std::vector<bool> stagPtMvmt;                                    ///< Flag to account for stagnation point movements.
-    Solver *tSolver;                                                 ///< Pseudo-time solver.
-    int solverExitCode;                                              ///< Exit code of viscous calculations.
-    std::vector<std::vector<std::vector<size_t>>> convergenceStatus; ///< Vector containing points that did not converged.
-
-protected:
-    fwk::Timers tms; ///< Internal timers
-
-public:
-    Driver(double _Re, double _Minf, double _CFL0, size_t nSections, double _span = 0., unsigned int verbose = 1);
-    ~Driver();
-    int run(unsigned int const couplIter);
-
-    // Getters.
-    double getxtr(size_t iSec, size_t iRegion) const { return sections[iSec][iRegion]->xtr; };
-    double getRe() const { return Re; };
-    std::vector<std::vector<double>> getSolution(size_t iSec);
-
-    // Setters.
-    void setGlobMesh(size_t iSec, size_t iRegion, std::vector<double> _x, std::vector<double> _y, std::vector<double> _z);
-    void setVelocities(size_t iSec, size_t iRegion, std::vector<double> _vx, std::vector<double> _vy, std::vector<double> _vz);
-    void setMe(size_t iSec, size_t iRegion, std::vector<double> _Me);
-    void setRhoe(size_t iSec, size_t iRegion, std::vector<double> _rhoe);
-    void setxxExt(size_t iSec, size_t iRegion, std::vector<double> _xxExt);
-    void setDeltaStarExt(size_t iSec, size_t iRegion, std::vector<double> _deltaStarExt);
-    void setUeOld(size_t iSec, size_t iRegion, std::vector<double> _ue);
-
-    // Extractors.
-    std::vector<double> extractBlowingVel(size_t iSec, size_t iRegion) const;
-    std::vector<double> extractxx(size_t iSec, size_t iRegion) const;
-    std::vector<double> extractDeltaStar(size_t iSec, size_t iRegion) const;
-    std::vector<double> extractxCoord(size_t iSec, size_t iRegion) const;
-    std::vector<double> extractUe(size_t iSec, size_t iRegion) const;
-
-private:
-    void checkWaves(size_t iPoint, BoundaryLayer *bl);
-    void averageTransition(size_t iPoint, BoundaryLayer *bl, int forced);
-    void computeSectionalDrag(std::vector<BoundaryLayer *> const &bl, size_t i);
-    void computeTotalDrag();
-    void computeBlwVel();
-    int outputStatus(double maxMach) const;
-};
-} // namespace vii
-#endif // DDRIVER_H

vii/src/DEdge.cpp

-#include "DEdge.h"
-
-using namespace vii;
-
-/**
- * @brief Construct a new Edge::Edge object.
- *
- */
-Edge::Edge()
-{
-    Me.resize(0, 0.);
-    vt.resize(0, 0.);
-    rhoe.resize(0, 0.);
-    deltaStar.resize(0, 0.);
-    vtOld.resize(0, 0.);
-}
-
-/**
- * @brief Return the maximum inviscid Mach number in the region.
- *
- * @return double
- */
-double Edge::getMaxMach() const
-{
-    if (Me.size() <= 0)
-        throw std::runtime_error("vii::Edge Invalid Mach number vector.");
-    double Mmax = 0.0;
-    for (size_t i = 0; i < Me.size(); ++i)
-        if (Me[i] > Mmax)
-            Mmax = Me[i];
-    return Mmax;
-}
-
-// Setters.
-void Edge::setMe(std::vector<double> const _Me)
-{
-    Me.resize(_Me.size(), 0.);
-    Me = _Me;
-}
-void Edge::setVt(std::vector<double> const _vt)
-{
-    vt.resize(_vt.size(), 0.);
-    vt = _vt;
-}
-void Edge::setRhoe(std::vector<double> const _rhoe)
-{
-    rhoe.resize(_rhoe.size(), 0.);
-    rhoe = _rhoe;
-}
-void Edge::setDeltaStar(std::vector<double> const _deltaStar)
-{
-    deltaStar.resize(_deltaStar.size(), 0.);
-    deltaStar = _deltaStar;
-}
-void Edge::setVtOld(std::vector<double> const _vtOld)
-{
-    vtOld.resize(_vtOld.size(), 0.);
-    vtOld = _vtOld;
-}

vii/src/DEdge.h

-#ifndef DEDGE_H
-#define DEDGE_H
-
-#include "vii.h"
-#include <vector>
-#include <iostream>
-
-namespace vii
-{
-
-/**
- * @brief Inviscid quantities at the edge of the boundary layer of a BoundaryLayer region.
- * @author Paul Dechamps
- */
-class VII_API Edge
-{
-private:
-    std::vector<double> Me;        ///< Mach number.
-    std::vector<double> vt;        ///< Velocity.
-    std::vector<double> rhoe;      ///< Density.
-    std::vector<double> deltaStar; ///< Displacement thickness.
-    std::vector<double> vtOld;     ///< Previous velocity.
-
-public:
-    Edge();
-    ~Edge(){};
-
-    // Getters.
-    double getMe(size_t iPoint) const { return Me[iPoint]; };
-    double getVt(size_t iPoint) const { return vt[iPoint]; };
-    double getRhoe(size_t iPoint) const { return rhoe[iPoint]; };
-    double getDeltaStar(size_t iPoint) const { return deltaStar[iPoint]; };
-    double getVtOld(size_t iPoint) const { return vtOld[iPoint]; };
-    double getMaxMach() const;
-
-    // Setters.
-    void setMe(std::vector<double> const _Me);
-    void setVt(std::vector<double> const _vt);
-    void setRhoe(std::vector<double> const _rhoe);
-    void setDeltaStar(std::vector<double> const _deltaStar);
-    void setVtOld(std::vector<double> const _vtOld);
-};
-} // namespace vii
-#endif // DEDGE_H

vii/src/DFluxes.cpp

-#include "DFluxes.h"
-#include "DBoundaryLayer.h"
-#include <Eigen/Dense>
-
-using namespace vii;
-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;
-}
\ No newline at end of file

vii/src/DFluxes.h

-#ifndef DFLUXES_H
-#define DFLUXES_H
-
-#include "vii.h"
-#include "DBoundaryLayer.h"
-#include <Eigen/Dense>
-
-namespace vii
-{
-
-/**
- * @brief Residual and Jacobian computation of the unsteady integral boundary layer equations.
- * @author Paul Dechamps
- */
-class VII_API Fluxes
-{
-public:
-    double Re; ///< Freestream Reynolds number.
-
-public:
-    Fluxes(double _Re) : Re(_Re){};
-    ~Fluxes(){};
-    Eigen::VectorXd computeResiduals(size_t iPoint, BoundaryLayer *bl);
-    Eigen::MatrixXd computeJacobian(size_t iPoint, BoundaryLayer *bl, Eigen::VectorXd const &sysRes, double eps);
-
-private:
-    Eigen::VectorXd blLaws(size_t iPoint, BoundaryLayer *bl, std::vector<double> u);
-    double amplificationRoutine(double Hk, double Ret, double theta) const;
-};
-} // namespace vii
-#endif // DFLUXES_H

vii/src/DSolver.cpp

-#include "DSolver.h"
-#include "DBoundaryLayer.h"
-#include "DFluxes.h"
-#include <Eigen/Dense>
-
-using namespace Eigen;
-using namespace tbox;
-using namespace vii;
-
-/**
- * @brief Construct a new Time Solver:: Time Solver object.
- *
- * @param _CFL0 Initial CFL number.
- * @param _Minf Freestream Mach number.
- * @param Re Freestream Reynolds number.
- * @param _maxIt Maximum number of iterations.
- * @param _tol Convergence tolerance.
- * @param _itFrozenJac Number of iterations between which the Jacobian is frozen.
- */
-Solver::Solver(double _CFL0, double _Minf, double Re, unsigned int _maxIt, double _tol, unsigned int _itFrozenJac)
-{
-    CFL0 = _CFL0;
-    maxIt = _maxIt;
-    tol = _tol;
-    itFrozenJac0 = _itFrozenJac;
-    Minf = std::max(_Minf, 0.1);
-    initSln = true;
-    sysMatrix = new Fluxes(Re);
-}
-
-/**
- * @brief Destroy the Time Solver:: Time Solver object.
- *
- */
-Solver::~Solver()
-{
-    delete sysMatrix;
-    std::cout << "~vii::Solver()\n";
-}
-
-/**
- * @brief Impose initial condition at one point.
- *
- * @param iPoint Marker id.
- * @param bl Boundary layer region.
- */
-void Solver::initialCondition(size_t iPoint, BoundaryLayer *bl)
-{
-    size_t nVar = bl->getnVar();
-    if (initSln)
-    {
-        bl->u[iPoint * nVar + 0] = bl->u[(iPoint - 1) * nVar + 0];
-        bl->u[iPoint * nVar + 1] = bl->u[(iPoint - 1) * nVar + 1];
-        bl->u[iPoint * nVar + 2] = bl->u[(iPoint - 1) * nVar + 2];
-        bl->u[iPoint * nVar + 3] = bl->blEdge->getVt(iPoint);
-        bl->u[iPoint * nVar + 4] = bl->u[(iPoint - 1) * nVar + 4];
-
-        bl->Ret[iPoint] = bl->Ret[iPoint - 1];
-        bl->cd[iPoint] = bl->cd[iPoint - 1];
-        bl->cf[iPoint] = bl->cf[iPoint - 1];
-        bl->Hstar[iPoint] = bl->Hstar[iPoint - 1];
-        bl->Hstar2[iPoint] = bl->Hstar2[iPoint - 1];
-        bl->Hk[iPoint] = bl->Hk[iPoint - 1];
-        bl->ctEq[iPoint] = bl->ctEq[iPoint - 1];
-        bl->us[iPoint] = bl->us[iPoint - 1];
-        bl->delta[iPoint] = bl->delta[iPoint - 1];
-        bl->deltaStar[iPoint] = bl->deltaStar[iPoint - 1];
-    }
-    if (bl->regime[iPoint] == 1 && bl->u[iPoint * nVar + 4] <= 0)
-        bl->u[iPoint * nVar + 4] = 0.03;
-}
-
-/**
- * @brief Pseudo-time integration.
- *
- * @param iPoint Marker id.
- * @param bl Boundary layer region.
- * @return int
- */
-int Solver::integration(size_t iPoint, BoundaryLayer *bl)
-{
-    size_t nVar = bl->getnVar();
-
-    // Save initial condition.
-    std::vector<double> sln0(nVar, 0.);
-    for (size_t i = 0; i < nVar; ++i)
-        sln0[i] = bl->u[iPoint * nVar + i];
-
-    // Initialize time integration variables.
-    double dx = bl->mesh->getLoc(iPoint) - bl->mesh->getLoc(iPoint - 1);
-
-    // Initial time step.
-    double CFL = CFL0;
-    unsigned int itFrozenJac = itFrozenJac0;
-    double timeStep = setTimeStep(CFL, dx, bl->blEdge->getVt(iPoint));
-    MatrixXd IMatrix(5, 5);
-    IMatrix = MatrixXd::Identity(5, 5) / timeStep;
-
-    // Initial system residuals.
-    VectorXd sysRes = sysMatrix->computeResiduals(iPoint, bl);
-    double normSysRes0 = sysRes.norm();
-    double normSysRes = normSysRes0;
-
-    MatrixXd JacMatrix = MatrixXd::Zero(5, 5);
-    VectorXd slnIncr = VectorXd::Zero(5);
-
-    unsigned int innerIters = 0; // Inner (non-linear) iterations
-
-    while (innerIters < maxIt)
-    {
-        // Jacobian computation.
-        if (innerIters % itFrozenJac == 0)
-            JacMatrix = sysMatrix->computeJacobian(iPoint, bl, sysRes, 1e-8);
-
-        // Linear system.
-        slnIncr = (JacMatrix + IMatrix).partialPivLu().solve(-sysRes);
-
-        // Solution increment.
-        for (size_t k = 0; k < nVar; ++k)
-            bl->u[iPoint * nVar + k] += slnIncr(k);
-        bl->u[iPoint * nVar + 0] = std::max(bl->u[iPoint * nVar + 0], 1e-6);
-        bl->u[iPoint * nVar + 1] = std::max(bl->u[iPoint * nVar + 1], 1.0005);
-
-        sysRes = sysMatrix->computeResiduals(iPoint, bl);
-        normSysRes = (sysRes + slnIncr / timeStep).norm();
-
-        if (std::isnan(normSysRes))
-        {
-            resetSolution(iPoint, bl, sln0, true);
-            sysRes = sysMatrix->computeResiduals(iPoint, bl);
-            CFL = 1.0;
-            timeStep = setTimeStep(CFL, dx, bl->blEdge->getVt(iPoint));
-            IMatrix = MatrixXd::Identity(5, 5) / timeStep;
-            normSysRes = (sysRes + slnIncr / timeStep).norm();
-            itFrozenJac = 1;
-        }
-
-        if (normSysRes <= tol * normSysRes0)
-            return 0; // Successfull exit.
-
-        // CFL adaptation.
-        CFL = std::max(CFL0 * pow(normSysRes0 / normSysRes, 0.7), 0.1);
-        timeStep = setTimeStep(CFL, dx, bl->blEdge->getVt(iPoint));
-        IMatrix = MatrixXd::Identity(5, 5) / timeStep;
-
-        innerIters++;
-    }
-
-    if (std::isnan(normSysRes) || normSysRes / normSysRes0 > 1e-3)
-    {
-        resetSolution(iPoint, bl, sln0);
-        return -1;
-    }
-    return innerIters;
-}
-
-/**
- * @brief Set time step.
- *
- * @param CFL CFL number.
- * @param dx Cell size.
- * @param invVel Inviscid velocity.
- * @return double
- */
-double Solver::setTimeStep(double CFL, double dx, double invVel) const
-{
-    // Speed of sound.
-    double gradU2 = (invVel * invVel);
-    double soundSpeed = computeSoundSpeed(gradU2);
-
-    // Velocity of the fastest travelling wave.
-    double advVel = soundSpeed + invVel;
-
-    // Time step computation.
-    return CFL * dx / advVel;
-}
-
-/**
- * @brief Compute the speed of sound in a cell.
- *
- * @param gradU2 Inviscid velocity squared in the cell.
- * @return double
- */
-double Solver::computeSoundSpeed(double gradU2) const
-{
-    return sqrt(1 / (Minf * Minf) + 0.5 * (gamma - 1) * (1 - gradU2));
-}
-
-/**
- * @brief Resets the solution of the point wrt its initial condition (sln0) or the previous point.
- *
- * @param iPoint Marker id.
- * @param bl Boundary layer region.
- * @param sln0 Initial solution at the point.
- * @param usePrevPoint if 1, use the solution at previous point instead of sln0.
- */
-void Solver::resetSolution(size_t iPoint, BoundaryLayer *bl, std::vector<double> sln0, bool usePrevPoint)
-{
-    size_t nVar = bl->getnVar();
-
-    // Reset solution.
-    if (usePrevPoint)
-        for (size_t k = 0; k < nVar; ++k)
-            bl->u[iPoint * nVar + k] = bl->u[(iPoint - 1) * nVar + k];
-    else
-        for (size_t k = 0; k < nVar; ++k)
-            bl->u[iPoint * nVar + k] = sln0[k];
-
-    // Reset closures.
-    bl->Ret[iPoint] = std::max(bl->u[iPoint * nVar + 3] * bl->u[iPoint * nVar + 0] * sysMatrix->Re, 1.0); // Reynolds number based on theta.
-    bl->deltaStar[iPoint] = bl->u[iPoint * nVar + 0] * bl->u[iPoint * nVar + 1];                          // Displacement thickness.
-
-    if (bl->regime[iPoint] == 0)
-    {
-        std::vector<double> lamParam(8, 0.);
-        bl->closSolver->laminarClosures(lamParam, bl->u[iPoint * nVar + 0], bl->u[iPoint * nVar + 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)
-    {
-        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];
-    }
-}
-
-void Solver::setCFL0(double _CFL0)
-{
-    if (_CFL0 > 0)
-        CFL0 = _CFL0;
-    else
-        throw std::runtime_error("Negative CFL numbers are not allowed.\n");
-}
\ No newline at end of file

vii/src/DSolver.h

-#ifndef DSOLVER_H
-#define DSOLVER_H
-
-#include "vii.h"
-#include "DFluxes.h"
-
-namespace vii
-{
-
-/**
- * @brief Pseudo-time integration.
- * @author Paul Dechamps
- */
-class VII_API Solver
-{
-
-private:
-    double CFL0;               ///< Initial CFL number.
-    unsigned int maxIt;        ///< Maximum number of iterations.
-    double tol;                ///< Convergence tolerance.
-    unsigned int itFrozenJac0; ///< Number of iterations between which the Jacobian is frozen.
-    bool initSln;              ///< Flag to initialize the solution at the point.
-    double const gamma = 1.4;  ///< Air heat capacity ratio.
-    double Minf;               ///< Freestream Mach number.
-
-    Fluxes *sysMatrix;
-
-public:
-    Solver(double _CFL0, double _Minf, double Re, unsigned int _maxIt = 100, double _tol = 1e-8, unsigned int _itFrozenJac = 5);
-    ~Solver();
-    void initialCondition(size_t iPoint, BoundaryLayer *bl);
-    int integration(size_t iPoint, BoundaryLayer *bl);
-
-    // Getters.
-    double getCFL0() const { return CFL0; };
-    unsigned int getmaxIt() const { return maxIt; };
-    unsigned int getitFrozenJac() const { return itFrozenJac0; };
-
-    // Setters.
-    void setCFL0(double _CFL0);
-    void setmaxIt(unsigned int _maxIt) { maxIt = _maxIt; };
-    void setitFrozenJac(unsigned int _itFrozenJac) { itFrozenJac0 = _itFrozenJac; };
-    void setinitSln(bool _initSln) { initSln = _initSln; };
-
-private:
-    double setTimeStep(double CFL, double dx, double invVel) const;
-    double computeSoundSpeed(double gradU2) const;
-    void resetSolution(size_t iPoint, BoundaryLayer *bl, std::vector<double> Sln0, bool usePrevPoint = false);
-};
-} // namespace vii
-#endif // DSOLVER_H

vii/src/vii.h

-/*
- * Copyright 2020 University of Liège
- *
- * Licensed under the Apache License, Version 2.0 (the "License");
- * you may not use this file except in compliance with the License.
- * You may obtain a copy of the License at
- *
- *     http://www.apache.org/licenses/LICENSE-2.0
- *
- * Unless required by applicable law or agreed to in writing, software
- * distributed under the License is distributed on an "AS IS" BASIS,
- * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
- * See the License for the specific language governing permissions and
- * limitations under the License.
- */
-
-// Global header of the "vii" module.
-
-#ifndef VII_H
-#define VII_H
-
-#if defined(WIN32)
-#ifdef vii_EXPORTS
-#define VII_API __declspec(dllexport)
-#else
-#define VII_API __declspec(dllimport)
-#endif
-#else
-#define VII_API
-#endif
-
-#include "tbox.h"
-
-/**
- * @brief Namespace for vii module
- */
-namespace vii
-{
-
-// Solvers
-class Driver;
-class Solver;
-
-// Data structures
-class BoundaryLayer;
-class Discretization;
-class Edge;
-
-// Other
-class Closures;
-} // namespace vii
-
-#endif // VII_H

vii/tests/bli2D.py

-#!/usr/bin/env python3
-# -*- coding: utf-8 -*-
-
-# Copyright 2022 University of Liège
-#
-# Licensed under the Apache License, Version 2.0 (the "License");
-# you may not use this file except in compliance with the License.
-# You may obtain a copy of the License at
-#
-#     http://www.apache.org/licenses/LICENSE-2.0
-#
-# Unless required by applicable law or agreed to in writing, software
-# distributed under the License is distributed on an "AS IS" BASIS,
-# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
-# See the License for the specific language governing permissions and
-# limitations under the License.
-
-
-# @author Paul Dechamps
-# @date 2022
-# Test the vii implementation. The test case is a compressible attached transitional flow.
-# Tested functionalities;
-# - Time integration.
-# - Two time-marching techniques (agressive and safe).
-# - Transition routines.
-# - Compressible flow routines.
-# - Laminar and turbulent flow.
-#
-# Untested functionalities.
-# - Separation routines.
-# - Multiple failure case at one iteration.
-# - Response to unconverged inviscid solver.
-
-# Imports.
-import dart.utils as invUtils
-import dart.default as invDefault
-
-import vii.coupler as viiCoupler
-import vii.utils as viscUtils
-
-import tbox.utils as tboxU
-import fwk
-from fwk.testing import *
-from fwk.coloring import ccolors
-
-import math
-
-def main():
-    # Timer.
-    tms = fwk.Timers()
-    tms['total'].start()
-
-    # Flow variables.
-    Re = 1e7
-    alpha = 2.*math.pi/180
-    M_inf = 0.2
-    # Numerical parameters.
-    CFL0 = 1
-
-    # Mesh the geometry.
-    print(ccolors.ANSI_BLUE + 'PyMeshing...' + ccolors.ANSI_RESET)
-    tms['msh'].start()
-
-    pars = {'xLgt' : 5, 'yLgt' : 5, 'msF' : 1, 'msTe' : 0.01, 'msLe' : 0.01}
-    msh, gmshWritter = invDefault.mesh(2, 'models/n0012.geo', pars, ['field', 'airfoil', 'downstream'])
-    tms['msh'].stop()
-
-    # Set the problem and add medium, initial/boundary conditions.
-    tms['pre'].start()
-    pbl, _, _, bnd, blw = invDefault.problem(msh, 2, alpha, 0., M_inf, 1.0, 1.0, 0.25, 0., 0., 'airfoil', te = 'te', vsc = True, dbc=True)
-    tms['pre'].stop()
-
-    # Solve coupled problem.
-    print(ccolors.ANSI_BLUE + 'PySolving...' + ccolors.ANSI_RESET)
-    tms['solver'].start()
-
-    vSolver = viscUtils.initBL(Re, M_inf, CFL0, 1)
-    iSolverAPI = viscUtils.initDart(pbl, blw, vSolver)
-    Coupler = viiCoupler.Coupler(iSolverAPI, vSolver, _maxCouplIter=20, _couplTol=1e-4, _iterPrint=1, _resetInv=False)
-    Coupler.run()
-
-    # Save results.
-    iSolverAPI.save(gmshWritter)
-    Cp = invUtils.extract(bnd.groups[0].tag.elems, iSolverAPI.solver.Cp)
-    tboxU.write(Cp, 'Cp.dat', '%1.5e', ', ', 'x, y, z, Cp', '')
-
-    print('')
-    # Display results.
-    print(ccolors.ANSI_BLUE + 'PyRes...' + ccolors.ANSI_RESET)
-    print('      Re        M    alpha       Cl       Cd      Cdp      Cdf       Cm')
-    print('{0:6.1f}e6 {1:8.2f} {2:8.1f} {3:8.4f} {4:8.4f} {5:8.4f} {6:8.4f} {7:8.4f}'.format(Re/1e6, M_inf, alpha*180/math.pi, iSolverAPI.getCl(), vSolver.Cdt, vSolver.Cdp, vSolver.Cdf, iSolverAPI.getCm()))
-
-    # Write results to file.
-    vSolution = viscUtils.getSolution(vSolver)
-    viscUtils.writeFile(vSolution, vSolver.getRe())
-
-    # Display timers.
-    tms['total'].stop()
-    print(ccolors.ANSI_BLUE + 'PyTiming...' + ccolors.ANSI_RESET)
-    print('CPU statistics')
-    print(tms)
-    print('SOLVERS statistics')
-    print(Coupler.tms)
-
-    # Plot results.
-    tboxU.plot(Cp[:,0], Cp[:,3], {'xlabel': 'x', 'ylabel': 'Cp', 'title': 'Cl = {0:.{3}f}, Cd = {1:.{3}f}, Cm = {2:.{3}f}'.format(iSolverAPI.getCl(), vSolver.Cdt, iSolverAPI.getCm(), 4), 'invert': True})
-    tboxU.plot(vSolution['x/c'], vSolution['Cf'], {'xlabel': '$x/c$', 'ylabel': '$c_f$', 'title': 'Cdt = {0:.{3}f}, Cdf = {1:.{3}f}, Cdp = {2:.{3}f}'.format(vSolver.Cdt, vSolver.Cdf, vSolver.Cdp, 4), 'xlim': [0, 1]})
-
-    # Check results.
-    # Test for n0012 (le=0.01, te=0.01), alpha=2, M=.2, Re=1e7.
-    print(ccolors.ANSI_BLUE + 'PyTesting...' + ccolors.ANSI_RESET)
-    tests = CTests()
-    tests.add(CTest('Cl', iSolverAPI.getCl(), 0.234, 5e-3)) # XFOIL 0.2325
-    tests.add(CTest('Cd', vSolution['Cdt'], 0.0055, 1e-3, forceabs=True)) # XFOIL 0.00531
-    tests.add(CTest('Cdp', vSolution['Cdp'], 0.0012, 1e-3, forceabs=True)) # XFOIL 0.00084, Cdf = 0.00447
-    tests.add(CTest('xtr_top', vSolution['xtrT'], 0.21, 0.03, forceabs=True)) # XFOIL 0.1877
-    tests.add(CTest('xtr_bot', vSolution['xtrB'], 0.50, 0.03, forceabs=True)) # XFOIL 0.4998
-    tests.run()
-
-    # eof
-    print('')
-
-if __name__ == "__main__":
-    main()

vii/utils.py

-import vii
-from fwk.wutils import parseargs
-from fwk.coloring import ccolors
-import numpy as np
-
-def initBL(Re, Minf, CFL0, nSections, span=0, verb=None):
-    """ Initialize boundary layer solver.
-    
-    Params
-    ------
-    - Re: Flow Reynolds number.
-    - Minf: Freestream Mach number.
-    - CFL0: Initial CFL number for time integration.
-    - nSections: Number of sections in the domain.
-    - span: Wing span (not used for 2D computations.
-    - verb: Verbosity level of the solver.
-    
-    Return
-    ------
-    - solver: vii::Driver class.
-    """
-    if Re<=0.:
-        print(ccolors.ANSI_RED, "dart::vUtils Error : Negative Reynolds number.", Re, ccolors.ANSI_RESET)
-        raise RuntimeError("Invalid parameter")
-    if Minf<0.:
-        print(ccolors.ANSI_RED, "dart::vUtils Error : Negative Mach number.", Minf, ccolors.ANSI_RESET)
-        raise RuntimeError("Invalid parameter")
-    elif Minf>=1.:
-        print(ccolors.ANSI_YELLOW, "dart::vUtils Warning : (Super)sonic freestream Mach number.", Minf, ccolors.ANSI_RESET)
-    if nSections < 0:
-        print(ccolors.ANSI_RED, "dart::vUtils Fatal error : Negative number of sections.", nSections, ccolors.ANSI_RESET)
-        raise RuntimeError("Invalid parameter")
-    if verb is None:
-      args = parseargs()
-      _verbose = args.verb
-    else:
-      if not(0<=verb<=3):
-        print(ccolors.ANSI_RED, "dart::vUtils Fatal error : verbose not in valid range.", _verbose, ccolors.ANSI_RESET)
-        raise RuntimeError("Invalid parameter")
-      else:
-        _verbose = verb
-    solver = vii.Driver(Re, Minf, CFL0, nSections, span, verbose=_verbose)
-    return solver
-
-def initDart(pbl, blw, vSolver, iters = 25):
-    """Initialize Newton solver and create the interface object
-       between dart and the viscous solver.
-
-    Parameters
-    ----------
-    pbl (dart.Problem object): Problem formulation.
-    blw (np.ndarray(dart.Blowing, dart.Blowing)): Blowing boundaries for (airfoil, wake).
-    vSolver (vii.Driver object): Viscous solver.
-
-    Return
-    ------
-    solverAPI (DartInterface object): Interface between Dart and the viscous solver.
-    """
-    import dart
-    from vii.interfaces.dart.DartInterface import DartInterface as dInterface
-    from tbox.solvers import LinearSolver
-
-    args = parseargs()
-    dartSolver = dart.Newton(pbl, LinearSolver().pardiso(), nthrds=args.k, vrb=args.verb, mIt=iters)
-    solverAPI = dInterface(dartSolver, vSolver, blw)
-    return solverAPI
-
-def mesh(file, pars):
-    """Initialize mesh and mesh writer
-
-    Parameters
-    ----------
-    file : str
-        file contaning mesh (w/o extention)
-    pars : dict
-        parameters for mesh
-    """
-    import tbox.gmsh as tmsh
-    # Load the mesh.
-    msh = tmsh.MeshLoader(file,__file__).execute(**pars)
-    return msh
-
-
-def getSolution(vSolver, iSec=0):
-    """Extract viscous solution.
-    Args
-    ----
-    - vSolver: vii::Driver class. Drive of the viscous calculations
-    - iSec (int): Marker of the section (default: 0)
-    """
-    if iSec<0:
-        raise RuntimeError("dart::viscU Invalid section number", iSec)
-    
-    solverOutput = vSolver.getSolution(iSec)
-
-    sln = {'theta'    : solverOutput[0],
-            'H'        : solverOutput[1],
-            'N'        : solverOutput[2],
-            'ue'       : solverOutput[3],
-            'Ct'       : solverOutput[4],
-            'deltaStar': solverOutput[5],
-            'Ret'      : solverOutput[6],
-            'Cd'       : solverOutput[7],
-            'Cf'       : np.asarray(solverOutput[8])*np.asarray(solverOutput[3])**2,
-            'Hstar'    : solverOutput[9],
-            'Hstar2'   : solverOutput[10],
-            'Hk'       : solverOutput[11],
-            'Cteq'     : solverOutput[12],
-            'Us'       : solverOutput[13],
-            'delta'    : solverOutput[14],
-            'x/c'      : solverOutput[15],
-            'xtrT'     : vSolver.getxtr(iSec, 0),
-            'xtrB'     : vSolver.getxtr(iSec, 1),
-            'Cdt'      : vSolver.Cdt,
-            'Cdf'      : vSolver.Cdf,
-            'Cdp'      : vSolver.Cdp}
-    return sln
-
-def writeFile(wData, Re, toW=['deltaStar', 'H', 'Cf', 'Ct', 'ue', 'delta']):
-    """Write the results in airfoil_viscous.dat
-    """
-    # Write
-    print('Writing file: airfoil_viscous.dat...', end = '')
-    f = open('airfoil_viscous.dat', 'w+')
-
-    f.write('$Aerodynamic coefficients\n')
-    f.write('             Re              Cd             Cdp             Cdf         xtr_top         xtr_bot\n')
-    f.write('{0:15.6f} {1:15.6f} {2:15.6f} {3:15.6f} {4:15.6f} {5:15.6f}\n'.format(Re, wData['Cdt'], wData['Cdp'],wData['Cdf'], wData['xtrT'], wData['xtrB']))
-    f.write('$Boundary layer variables\n')
-    f.write('            x/c')
-
-    for s in toW:
-        f.write(' {0:>15s}'.format(s))
-    f.write('\n')
-
-    for i in range(len(wData['x/c'])):
-        f.write('{0:15.6f}'.format(wData['x/c'][i]))
-        for s in toW:
-            f.write(' {0:15.6f}'.format(wData[s][i]))
-            f.write('\n')
-
-    f.close()
-    print('done.')
\ No newline at end of file