-
Paul Dechamps authored
Ensure all numbers are doubles in Closures, Fluxes and Driver & debug info update in BoundaryLayer & ensure nIter in 3D rae test
Paul Dechamps authoredEnsure all numbers are doubles in Closures, Fluxes and Driver & debug info update in BoundaryLayer & ensure nIter in 3D rae test
DDriver.cpp 33.56 KiB
#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;
}