From e550a50462d624931ddb26560db64018ec2c7aae Mon Sep 17 00:00:00 2001
From: Paul <paul.dechamps@student.uliege.be>
Date: Thu, 21 Jul 2022 13:37:12 +0200
Subject: [PATCH] Refined Quasi-3D mesh + small changes
---
vii/pyVII/vCoupler.py | 3 -
vii/pyVII/vCoupler2.py | 10 +-
vii/pyVII/vInterpolator.py | 14 +-
vii/pyVII/vUtils.py | 8 +-
vii/src/wViscSolver.cpp | 243 ++++++++++++++++------------------
vii/src/wViscSolver.h | 4 +
vii/tests/bli25.py | 15 ++-
vii/tests/bli25_RandomTest.py | 161 ++++++++++++++++++++++
8 files changed, 307 insertions(+), 151 deletions(-)
create mode 100644 vii/tests/bli25_RandomTest.py
diff --git a/vii/pyVII/vCoupler.py b/vii/pyVII/vCoupler.py
index bab595e..5268b75 100644
--- a/vii/pyVII/vCoupler.py
+++ b/vii/pyVII/vCoupler.py
@@ -171,9 +171,6 @@ class Coupler:
self.group[0].xx = xxAF[self.group[0].connectListNodes.argsort()]
self.group[1].xx = xxWK[self.group[1].connectListNodes.argsort()]
- plt.plot(np.linspace(0,2, len(self.group[0].u)),self.group[0].u)
- plt.show()
-
for n in range(0, len(self.group)):
for i in range(0, self.group[n].nE):
self.group[n].boundary.setU(i, self.group[n].u[i])
diff --git a/vii/pyVII/vCoupler2.py b/vii/pyVII/vCoupler2.py
index 05090fb..7e9aad3 100644
--- a/vii/pyVII/vCoupler2.py
+++ b/vii/pyVII/vCoupler2.py
@@ -47,10 +47,11 @@ class Coupler:
self.nElmUpperPrev = 0
couplIter = 0
cdPrev = 0
- print('- AoA: {:<2.2f} Mach: {:<1.1f} Re: {:<2.1f}e6'
- .format(self.iSolverAPI.GetAlpha()*180/M_PI, self.iSolverAPI.GetMinf(), self.vSolver.GetRe()/1e6))
- print('{:>5s}| {:>10s} {:>10s} | {:>8s} {:>6s} {:>8s}'.format('It', 'Cl', 'Cd', 'xtrT', 'xtrB', 'Error'))
- print(' --------------------------------------------------------')
+ if self.iSolverAPI.iSolver.verbose == 0 and self.vSolver.verbose == 0:
+ print('- AoA: {:<2.2f} Mach: {:<1.1f} Re: {:<2.1f}e6'
+ .format(self.iSolverAPI.GetAlpha()*180/M_PI, self.iSolverAPI.GetMinf(), self.vSolver.GetRe()/1e6))
+ print('{:>5s}| {:>10s} {:>10s} | {:>8s} {:>6s} {:>8s}'.format('It', 'Cl', 'Cd', 'xtrT', 'xtrB', 'Error'))
+ print(' --------------------------------------------------------')
while couplIter < self.maxCouplIter:
# Run inviscid solver
@@ -88,7 +89,6 @@ class Coupler:
cdPrev = self.vSolver.Cdt
-
print('- {:>4.0f}| {:>10.5f} {:>10.5f} | {:>8.3f} {:>6.3f} {:>8.2f}'.format(couplIter,
self.iSolverAPI.GetCl(), self.vSolver.Cdt, self.vSolver.Getxtr(0,0), self.vSolver.Getxtr(0,1), np.log10(error)))
diff --git a/vii/pyVII/vInterpolator.py b/vii/pyVII/vInterpolator.py
index 5c852d3..c252a29 100644
--- a/vii/pyVII/vInterpolator.py
+++ b/vii/pyVII/vInterpolator.py
@@ -57,8 +57,6 @@ class Interpolator:
dummy[:,1] += sym
self.idata[iReg] = np.row_stack((self.idata[iReg], dummy))
self.vdata, self.vcg = self.GetMesh(_vmsh, visc=1)
- plt.plot(self.vdata[0][:,0], self.vdata[0][:,1],'o')
- plt.show()
self.Xp, self.Yp, self.Zp, self.sortedRows, self.Xc, self.Yc, self.Zc, self.sortedElems = self.SortMesh(self.vdata, self.vcg)
self.tms['Mesh sort'].stop()
print(self.tms)
@@ -74,9 +72,9 @@ class Interpolator:
ax.set_zlabel('z')
#plt.gca().invert_yaxis()
plt.axis('off')
- plt.show()"""
+ plt.show()
- """fig = plt.figure()
+ fig = plt.figure()
ax = fig.add_subplot(projection='3d')
ax.scatter(self.idata[0][:,0], self.idata[0][:,1], self.idata[0][:,2], s=0.3)
ax.scatter(self.icg[0][:,0], self.icg[0][:,1], self.icg[0][:,2], s=0.3)
@@ -128,6 +126,8 @@ class Interpolator:
U[iReg][iNode,:] = [self.iSolver.U[row][0], self.iSolver.U[row][1], self.iSolver.U[row][2]]
M[iReg][iNode] = self.iSolver.M[row]
Rho[iReg][iNode] = self.iSolver.rho[row]
+
+ #self.URani = U
nD = self.cfg['nDim']
Ux = [self.__RbfToSections(self.idata[0][:,:nD], U[0][:,0], self.vdata[0][:,:nD]), self.__RbfToSections(self.idata[1][:, 0:2], U[1][:,0], self.vdata[1][:,0:2])]
@@ -135,6 +135,8 @@ class Interpolator:
Uz = [self.__RbfToSections(self.idata[0][:,:nD], U[0][:,2], self.vdata[0][:,:nD]), self.__RbfToSections(self.idata[1][:, 0:2], U[1][:,2], self.vdata[1][:,0:2])]
Me = [self.__RbfToSections(self.idata[0][:,:nD], M[0], self.vdata[0][:,:nD]), self.__RbfToSections(self.idata[1][:, 0:2], M[1], self.vdata[1][:,0:2])]
Rhoe = [self.__RbfToSections(self.idata[0][:,:nD], Rho[0], self.vdata[0][:,:nD]), self.__RbfToSections(self.idata[1][:, 0:2], Rho[1], self.vdata[1][:,0:2])]
+
+ #self.URani = Ux[0]
for iSec, y in enumerate(self.cfg['Sections']):
for iReg in range(2):
@@ -236,7 +238,7 @@ class Interpolator:
blw[iSec][iReg] = np.concatenate((blw[iSec][iReg], blwUp[::-1], blwLw))
elif iReg == 1:
blw[iSec][iReg] = np.concatenate((blw[iSec][iReg], vSolver.ExtractBlowingVel(iSec, 2)))
- #blw[iSec][iReg] = np.clip(blw[iSec][iReg], -0.01, 0.01)
+
for iReg in range(2):
# Find index of all Sections in Y
@@ -322,6 +324,8 @@ class Interpolator:
self.tms['Interpolation'].stop()
+ if iReg == 0:
+ self.URani = mappedBlw
if iReg==0 or iReg==1:
for iElm, e in enumerate(self.blw[iReg].tag.elems):
if np.size(mappedBlw[self.icg[iReg][:,3]==e.no])>1:
diff --git a/vii/pyVII/vUtils.py b/vii/pyVII/vUtils.py
index e344671..b10033f 100644
--- a/vii/pyVII/vUtils.py
+++ b/vii/pyVII/vUtils.py
@@ -2,7 +2,7 @@ import vii
from fwk.wutils import parseargs
from fwk.coloring import ccolors
-def initBL(Re, Minf, CFL0, nSections, span=0, verbose=1):
+def initBL(Re, Minf, CFL0, nSections, span=0, _verbose=1):
if Re<=0.:
print(ccolors.ANSI_RED, "dart::vUtils Error : Negative Reynolds number.", Re, ccolors.ANSI_RESET)
raise RuntimeError("Invalid parameter")
@@ -14,12 +14,12 @@ def initBL(Re, Minf, CFL0, nSections, span=0, verbose=1):
if nSections < 0:
print(ccolors.ANSI_RED, "dart::vUtils Fatal error : Negative number of sections.", nSections, ccolors.ANSI_RESET)
raise RuntimeError("Invalid parameter")
- if not(0<=verbose<=3):
- print(ccolors.ANSI_RED, "dart::vUtils Fatal error : verbose not in valid range.", verbose, ccolors.ANSI_RESET)
+ if not(0<=_verbose<=3):
+ print(ccolors.ANSI_RED, "dart::vUtils Fatal error : verbose not in valid range.", _verbose, ccolors.ANSI_RESET)
raise RuntimeError("Invalid parameter")
args = parseargs()
- solver = vii.ViscSolver(Re, Minf, CFL0, nSections, span, verbose=args.verb+1)
+ solver = vii.ViscSolver(Re, Minf, CFL0, nSections, span, verbose=_verbose)
return solver
diff --git a/vii/src/wViscSolver.cpp b/vii/src/wViscSolver.cpp
index a89e87a..adecfc1 100644
--- a/vii/src/wViscSolver.cpp
+++ b/vii/src/wViscSolver.cpp
@@ -15,6 +15,11 @@
#include "wTimeSolver.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"
+
/* --- Namespaces -------------------------------------------------------------------- */
using namespace vii;
@@ -29,15 +34,13 @@ ViscSolver::ViscSolver(double _Re, double _Minf, double _CFL0, unsigned int nSec
PrintBanner();
- /* Freestream parameters */
-
+ // Freestream parameters.
Re = _Re;
CFL0 = _CFL0;
verbose = _verbose;
- /* Initialzing boundary layer */
-
+ // Initialzing boundary layer
Sections.resize(nSections);
for (size_t iSec=0; iSec < nSections; ++iSec)
{
@@ -57,26 +60,21 @@ ViscSolver::ViscSolver(double _Re, double _Minf, double _CFL0, unsigned int nSec
Cdt = 0.; Cdf = 0.; Cdp = 0.;
Span = _Span;
- /* Pre/post processing flags */
-
+ // Pre/post processing flags
solverExitCode = 0;
stagPtMvmt.resize(Sections.size(), false);
- /* Temporal solver initialization */
-
+ // Temporal solver initialization.
tSolver = new TimeSolver(_CFL0, _Minf, Re);
- /* User input numerical parameters */
-
+ // User input 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;
-
-} // end ViscSolver
+}
/**
* @brief Viscous solver and subsequent dependancies destruction
@@ -90,7 +88,6 @@ ViscSolver::~ViscSolver()
delete Sections[iSec][i];
}
}
-
delete tSolver;
std::cout << "~ViscSolver()\n";
} // end ~ViscSolver
@@ -101,30 +98,28 @@ ViscSolver::~ViscSolver()
*/
int ViscSolver::Run(unsigned int const couplIter)
{
- bool lockTrans=false; /* Flagging activation of transition routines */
- int pointExitCode = 0; /* Output of pseudo time integration (0: converged; -1: unsuccessful, failed; >0: unsuccessful nb iter) */
+ bool lockTrans=false; /* Flagging activation of transition routines */
+ int pointExitCode = 0; /* Output of pseudo time integration (0: converged; -1: unsuccessful, failed; >0: unsuccessful nb iter) */
solverExitCode = 0;
+ double maxMach = 0.;
+
for (size_t iSec=0; iSec<Sections.size(); ++iSec)
{
+ tms["0-CheckIn"].start();
- if (verbose>0){std::cout << "Sec " << iSec << " ";}
- if (verbose>1){std::cout << "Mmax " << std::max(Sections[iSec][0]->invBnd->GetMaxMach(), Sections[iSec][1]->invBnd->GetMaxMach());}
-
- /* Check for stagnation point movement */
+ 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>1){std::cout << "Stagnation point moved by " << Sections[iSec][0]->mesh->GetDiffnElm() << " elements." << std::endl;}
+ if (verbose>2) std::cout << "Stagnation point moved by " << Sections[iSec][0]->mesh->GetDiffnElm() << " elements." << std::endl;
}
else
- {
stagPtMvmt[iSec] = false;
- }
-
- /* Set boundary conditions */
+ // 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) */
@@ -132,14 +127,17 @@ int ViscSolver::Run(unsigned int const couplIter)
Sections[iSec][0]->SetStagBC(Re);
Sections[iSec][1]->SetStagBC(Re);
- /* Loop over the regions (upper, lower and wake) */
+ tms["0-CheckIn"].stop();
+
+ // Loop over the regions (upper, lower and wake).
for (size_t iRegion = 0; iRegion < Sections[iSec].size(); ++iRegion)
{
+ if (Sections[iSec][iRegion]->invBnd->GetMaxMach() > maxMach)
+ maxMach = Sections[iSec][iRegion]->invBnd->GetMaxMach();
- if (verbose>0){std::cout << Sections[iSec][iRegion]->name << " ";}
-
- /* Check if safe mode is required */
+ 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);
@@ -150,9 +148,7 @@ int ViscSolver::Run(unsigned int const couplIter)
{
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);
@@ -160,14 +156,12 @@ int ViscSolver::Run(unsigned int const couplIter)
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]->invBnd->GetVt(iPoint);
- }
Sections[iSec][iRegion]->invBnd->SetVtOld(ueOld);
}
- /* Keyword annoncing iteration safety level */
- if (verbose>1){std::cout << "restart. ";}
+ // Keyword annoncing iteration safety level.
+ if (verbose>1) std::cout << "restart. ";
}
else
@@ -176,30 +170,25 @@ int ViscSolver::Run(unsigned int const couplIter)
tSolver->SetitFrozenJac(5);
tSolver->SetinitSln(false);
- /* Keyword annoncing iteration safety level */
- if (verbose>1){std::cout << "update. ";}
+ // Keyword annoncing iteration safety level.
+ if (verbose>1) std::cout << "update. ";
}
- /* Save number of points in each regions */
-
+ // Save number of points in each regions.
Sections[iSec][iRegion]->mesh->prevnMarkers = Sections[iSec][iRegion]->mesh->GetnMarkers();
- /* Reset regime */
-
+ // 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 // Wake
{
for (size_t k = 0; k < Sections[iSec][iRegion]->mesh->GetnMarkers(); k++)
Sections[iSec][iRegion]->Regime[k] = 1;
- /* Impose wake boundary condition */
+ // Impose wake boundary condition.
std::vector<double> UpperCond(Sections[iSec][iRegion]->GetnVar(), 0.);
std::vector<double> LowerCond(Sections[iSec][iRegion]->GetnVar(), 0.);
for (size_t k=0; k<UpperCond.size(); ++k)
@@ -209,7 +198,7 @@ int ViscSolver::Run(unsigned int const couplIter)
}
Sections[iSec][iRegion]->SetWakeBC(Re, UpperCond, LowerCond);
lockTrans = true;
- if (verbose>2){std::cout << "Wake BC imposed." <<std::endl;}
+ if (verbose>2) std::cout << "Wake BC imposed." <<std::endl;
}
// Loop over points
@@ -219,51 +208,55 @@ int ViscSolver::Run(unsigned int const couplIter)
// Initialize solution at point
tSolver->InitialCondition(iPoint, Sections[iSec][iRegion]);
// Solve equations
+ tms["1-TimeSolver"].start();
pointExitCode = tSolver->Integration(iPoint, Sections[iSec][iRegion]);
+ tms["1-TimeSolver"].stop();
// Unsucessfull convergence output
- if (pointExitCode != 0)
+ if (pointExitCode < 0)
{
if (verbose>1){std::cout << "Point " << iPoint << ": Convergence terminated with exitcode " << pointExitCode << std::endl;}
solverExitCode = -1; // Convergence failed
}
+ else if (pointExitCode > 0 and solverExitCode >=0)
+ solverExitCode = pointExitCode; // Not fully converged
else
- {
- if (verbose>2){std::cout << "Point " << iPoint << ": Successful convergence.\n";}
- }
-
- /* Transition */
+ if (verbose>2) std::cout << "Point " << iPoint << ": Successful convergence.\n";
+ // Transition
+ tms["2-Transition"].start();
if (!lockTrans)
{
- /* Check if perturbation waves are growing or not. */
-
+ // Check if perturbation waves are growing or not.
CheckWaves(iPoint, Sections[iSec][iRegion]);
- /* Amplification factor is compared to critical amplification factor. */
-
+ // 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>0)
+ if (verbose>1)
{
std::cout << std::fixed;
std::cout << std::setprecision(2);
std::cout << Sections[iSec][iRegion]->xtr
<< " (" << Sections[iSec][iRegion]->transMarker << ") ";
}
- if (verbose>2){std::cout << "" << std::endl;}
+ if (verbose>2)
+ std::cout << "" << std::endl;
lockTrans = true;
- } // End if
- } // End if (!lockTrans).
+ }
+ }
+ tms["2-Transition"].stop();
} // End for iPoints
+ tms["3-Blowing"].start();
Sections[iSec][iRegion]->ComputeBlwVel();
+ tms["3-Blowing"].stop();
if (verbose>2){"Blowing velocities of " + Sections[iSec][iRegion]->name + " side computed.";}
- } // End for iRegion
- if (verbose>0){std::cout << "\n";}
+ } // End for iRegion
+ if (verbose>1) std::cout << "\n";
} // End for iSections
for (size_t i=0; i<Sections.size(); ++i)
@@ -272,9 +265,39 @@ int ViscSolver::Run(unsigned int const couplIter)
}
ComputeTotalDrag();
- if (solverExitCode != 0 && verbose>1){std::cout << "some points did not converge. ";}
- if (verbose>0){std::cout << "\n" << std::endl;}
+ // 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]->xtr;
+ for (size_t iReg=0; iReg<meanTransitions.size(); ++iReg)
+ meanTransitions[iReg]/=Sections.size();
+ if (verbose > 0)
+ {
+ if (verbose > 1) std::cout << "\n" << std::endl;
+ std::cout << std::fixed << std::setprecision(2)
+ << "Viscous solver for Re " << Re/1e6 << "e6, Mach max "
+ << maxMach
+ << std::endl;
+
+ std::cout << std::setw(8) << "xTr_Up"
+ << std::setw(12) << "xTr_Lw"
+ << std::setw(12) << "Cd" << std::endl;
+ std::cout << "-----------------------------------\n";
+ std::cout << std::fixed << std::setprecision(5);
+ std::cout << std::setw(8) << meanTransitions[0]
+ << std::setw(12) << meanTransitions[1]
+ << std::setw(12) << Cdt << "\n";
+
+ if (verbose > 2)
+ {
+ std::cout << "\n" << std::endl;
+ std::cout << "Viscous solver CPU" << std::endl
+ << tms;
+ }
+ std::cout << "\n" << std::endl;
+ }
return solverExitCode; // exit code
}
@@ -283,16 +306,13 @@ int ViscSolver::Run(unsigned int const couplIter)
*/
void ViscSolver::CheckWaves(size_t iPoint, BLRegion *bl)
{
- /* Determine if the amplification waves start growing or not. */
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.
- }
}
}
@@ -301,11 +321,11 @@ void ViscSolver::CheckWaves(size_t iPoint, BLRegion *bl)
*/
void ViscSolver::AverageTransition(size_t iPoint, BLRegion *bl, int forced)
{
- /* Averages solution on transition marker */
+ // Averages solution on transition marker.
unsigned int nVar = bl->GetnVar();
- /* Save laminar solution. */
+ // Save laminar solution.
std::vector<double> lamSol(nVar, 0.0);
for (size_t k = 0; k < lamSol.size(); ++k)
@@ -313,31 +333,29 @@ void ViscSolver::AverageTransition(size_t iPoint, BLRegion *bl, int forced)
lamSol[k] = bl->U[iPoint * nVar + k];
}
- /* Set up turbulent portion boundary condition. */
+ // Set up turbulent portion boundary condition.
- bl->transMarker = iPoint; /* Save transition marker */
+ bl->transMarker = iPoint; // Save transition marker.
- /* Loop over remaining points in the region */
+ // 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 */
+ // 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);
- /* Percentage of the subsection corresponding to a turbulent flow */
+ // 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 */
+ // 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->invBnd->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 */
+ // 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.
@@ -345,19 +363,18 @@ void ViscSolver::AverageTransition(size_t iPoint, BLRegion *bl, int forced)
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;
- }
- else if(verbose>2){std::cout << "Sucessfull transition point turbulent computation.";}
+ else if(verbose>2)
+ std::cout << "Sucessfull transition point turbulent computation.";
- /* 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 */
+ // 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 */
- /* Recompute closures. (The turbulent BC @ iPoint - 1 does not influence laminar closure relations). */
+ // 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->invBnd->GetMe(iPoint - 1), bl->name);
bl->Hstar[iPoint - 1] = lamParam[0];
@@ -394,40 +411,29 @@ void ViscSolver::ComputeSectionalDrag(std::vector<BLRegion *> const &bl, int i)
unsigned int nVar = bl[0]->GetnVar();
unsigned int lastWkPt = (bl[2]->mesh->GetnMarkers() - 1) * nVar;
- /* Normalize friction and dissipation coefficients. */
-
+ // 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]->invBnd->GetVt(k) * bl[0]->invBnd->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]->invBnd->GetVt(k) * bl[1]->invBnd->GetVt(k) * bl[1]->Cf[k];
- }
-
- /* Compute friction drag coefficient. */
+ // 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;
- }
Cdf_sec[i] = Cdf_upper + Cdf_lower; // No friction in the wake
- /* Compute total drag coefficient. */
-
+ // Compute 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);
- /* Compute pressure drag coefficient. */
-
+ // Compute pressure drag coefficient.
Cdp_sec[i] = Cdt_sec[i] - Cdf_sec[i];
}
@@ -445,14 +451,12 @@ void ViscSolver::ComputeTotalDrag()
else
{
- std::cout << "0 " << Cdt_sec[0] << std::endl;
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)
{
- std::cout << k << " " << Cdt_sec[k] << std::endl;
dz = (Sections[k][0]->mesh->Getz(0) - Sections[k-1][0]->mesh->Getz(0));
Cdt += dz * Cdt_sec[k];
Cdp += dz * Cdp_sec[k];
@@ -465,9 +469,9 @@ void ViscSolver::ComputeTotalDrag()
void ViscSolver::SetGlobMesh(size_t iSec, size_t iRegion, std::vector<double> _x, std::vector<double> _y, std::vector<double> _z)
{
- if (verbose>1){std::cout << "Mesh of " << _x.size() << " elements on " << Sections[iSec][iRegion]->name << " side. ";}
+ if (verbose>2){std::cout << "Mesh of " << _x.size() << " elements on " << Sections[iSec][iRegion]->name << " side. ";}
Sections[iSec][iRegion]->SetMesh(_x, _y, _z);
- if (verbose>1){std::cout << "Mesh mapped to local frame of reference." << std::endl;}
+ if (verbose>3){std::cout << "Mesh mapped to local frame of reference." << std::endl;}
}
void ViscSolver::SetVelocities(size_t iSec, size_t iRegion, std::vector<double> vx, std::vector<double> vy, std::vector<double> vz)
@@ -495,27 +499,21 @@ void ViscSolver::SetVelocities(size_t iSec, size_t iRegion, std::vector<double>
void ViscSolver::SetMe(size_t iSec, size_t iRegion, std::vector<double> _Me)
{
if (_Me.size() != Sections[iSec][iRegion]->mesh->GetnMarkers())
- {
throw std::runtime_error("dart::ViscSolver Mach number vector is not coherent with the mesh.");
- }
Sections[iSec][iRegion]->invBnd->SetMe(_Me);
}
void ViscSolver::SetRhoe(size_t iSec, size_t iRegion, std::vector<double> _Rhoe)
{
if (_Rhoe.size() != Sections[iSec][iRegion]->mesh->GetnMarkers())
- {
throw std::runtime_error("dart::ViscSolver Density vector is not coherent with the mesh.");
- }
Sections[iSec][iRegion]->invBnd->SetRhoe(_Rhoe);
}
void ViscSolver::SetxxExt(size_t iSec, size_t iRegion, std::vector<double> _xxExt)
{
if (_xxExt.size() != Sections[iSec][iRegion]->mesh->GetnMarkers())
- {
throw std::runtime_error("dart::ViscSolver External mesh vector is not coherent with the mesh.");
- }
Sections[iSec][iRegion]->mesh->SetExt(_xxExt);
}
@@ -533,9 +531,7 @@ std::vector<double> ViscSolver::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;
}
@@ -547,9 +543,7 @@ std::vector<double> ViscSolver::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;
}
std::vector<double> ViscSolver::ExtractDeltaStar(size_t iSec, size_t iRegion) const
@@ -561,18 +555,14 @@ std::vector<double> ViscSolver::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;
}
void ViscSolver::SetUeOld(size_t iSec, size_t iRegion, std::vector<double> _ue)
{
if (_ue.size() != Sections[iSec][iRegion]->mesh->GetnMarkers())
- {
throw std::runtime_error("dart::ViscSolver External velocity vector is not coherent with the mesh.");
- }
Sections[iSec][iRegion]->invBnd->SetVtOld(_ue);
}
@@ -584,9 +574,7 @@ std::vector<std::vector<double>> ViscSolver::GetSolution(size_t iSec)
std::vector<std::vector<double>> Sln(16);
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)
@@ -655,10 +643,9 @@ std::vector<std::vector<double>> ViscSolver::GetSolution(size_t iSec)
Sln[14][nMarkersUp + nMarkersLw + iPoint] = Sections[iSec][2]->Delta[iPoint];
Sln[15][nMarkersUp + nMarkersLw + iPoint] = Sections[iSec][2]->mesh->Getx(iPoint);
-
}
- if (verbose>2){std::cout << "Solution structure sent." << std::endl;}
+ if (verbose>2) std::cout << "Solution structure sent." << std::endl;
return Sln;
diff --git a/vii/src/wViscSolver.h b/vii/src/wViscSolver.h
index a981d88..2a717f7 100644
--- a/vii/src/wViscSolver.h
+++ b/vii/src/wViscSolver.h
@@ -3,6 +3,7 @@
#include "vii.h"
#include "wBLRegion.h"
+#include "wTimers.h"
namespace vii
{
@@ -33,6 +34,9 @@ private:
double Span; /* Wing Span (0 if 2D case) */
+protected:
+ fwk::Timers tms; ///< internal timers
+
public:
ViscSolver(double _Re, double _Minf, double _CFL0, unsigned int nSections, double _Span=0., unsigned int verbose=1);
~ViscSolver();
diff --git a/vii/tests/bli25.py b/vii/tests/bli25.py
index c2129b9..b53fc23 100644
--- a/vii/tests/bli25.py
+++ b/vii/tests/bli25.py
@@ -64,7 +64,7 @@ def main():
c_ref = 1.0 # reference length
S_ref = spn # reference area
fms = 1.0 # farfield mesh size
- nms = 0.01 # nearfield mesh size
+ nms = 0.007 # nearfield mesh size
# mesh the geometry
print(ccolors.ANSI_BLUE + 'PyMeshing...' + ccolors.ANSI_RESET)
@@ -82,7 +82,13 @@ def main():
pbl, _, _, bnd, blw = floD.problem(msh, dim, alpha, 0., M_inf, S_ref, c_ref, 0.25, 0., 0., 'wing', vsc = True)
tms['pre'].stop()
- config = {'nDim': dim, 'Sections':[-0.1, -0.05, 0.0, 0.05, 0.1], 'span':spn, 'Sym':[-0.1, 0.1], 'rbftype': 'linear', 'smoothing': 1e-6, 'degree': 0, 'neighbors': 10}
+ config = {'nDim': dim,
+ 'Sections':[-0.1, -0.05, 0.0, 0.05, 0.1],
+ 'span':spn, 'Sym':[-0.1, 0.1],
+ 'rbftype': 'linear',
+ 'smoothing': 1e-6,
+ 'degree': 0,
+ 'neighbors': 10}
# solve problem
iSolver = floD.newton(pbl)
@@ -99,7 +105,7 @@ def main():
# VII
iSolverAPI = vInterpol.Interpolator(iSolver, blw, imsh, vmsh, config)
- vSolver = viscU.initBL(Re, M_inf, CFL0, len(config['Sections']), spn, 2)
+ vSolver = viscU.initBL(Re, M_inf, CFL0, len(config['Sections']), spn, 1)
coupler = viscC.Coupler(iSolverAPI, vSolver, _maxIters=200)
coupler.Run()
@@ -107,9 +113,6 @@ def main():
vSolution = viscU.GetSolution(0, vSolver)
viscU.WriteFile(vSolution, vSolver.GetRe())
- plt.plot(vSolution['x/c'], vSolution['theta'])
- plt.show()
-
# Write corrected Cp
iSolver.save(mshWriter)
dartU.writeSlices(msh.name, config['Sections'], 2)
diff --git a/vii/tests/bli25_RandomTest.py b/vii/tests/bli25_RandomTest.py
new file mode 100644
index 0000000..642d376
--- /dev/null
+++ b/vii/tests/bli25_RandomTest.py
@@ -0,0 +1,161 @@
+#!/usr/bin/env python3
+# -*- coding: utf-8 -*-
+
+# 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.
+
+
+## Compute lifting (linear or nonlinear) potential flow around a NACA 0012 rectangular wing
+# Adrien Crovato
+#
+# Test the nonlinear shock-capturing capability and the Kutta condition
+#
+# CAUTION
+# This test is provided to ensure that the solver works properly.
+# Mesh refinement may have to be performed to obtain physical results.
+
+import math
+import dart.default as floD
+import dart.utils as dartU # or from dartflo.dart import utils as dartU
+import tbox.utils as tboxU # or from dartflo.tbox import utils as tboxU
+import tboxVtk
+import numpy as np
+
+import vii.pyVII.vUtils as viscU
+import vii.pyVII.vCoupler2 as viscC
+import vii.pyVII.vInterpolator as vInterpol
+
+from matplotlib import pyplot as plt
+
+import fwk
+from fwk.testing import *
+from fwk.coloring import ccolors
+
+def main():
+ # timer
+ tms = fwk.Timers()
+ tms['total'].start()
+
+ # define flow variables
+ alpha = 2.*np.pi/180
+ M_inf = 0.715
+ dim = 3
+
+ CFL0 = 1
+ Re = 1e7
+ dim = 3
+
+ # define dimension and mesh size
+ spn = 0.2 # wing span
+ lgt = 6.0 # channel length
+ hgt = 6.0 # channel height
+ wdt = 3.0 # channel width
+ c_ref = 1.0 # reference length
+ S_ref = spn # reference area
+ fms = 1.0 # farfield mesh size
+ nms = 0.01 # nearfield mesh size
+
+ U = []
+ for i in range(5):
+ # mesh the geometry
+ print(ccolors.ANSI_BLUE + 'PyMeshing...' + ccolors.ANSI_RESET)
+ tms["msh"].start()
+ pars = {'spn' : spn, 'lgt' : lgt, 'wdt' : wdt, 'hgt' : hgt, 'msN' : nms, 'msF' : fms}
+ imsh = viscU.mesh('/home/paul/lab/dartflo/dart/models/rae_25.msh', pars)
+ vmsh = viscU.mesh('/home/paul/lab/dartflo/vii/models/rae_25_visc.geo', pars)
+ msh, gmshWriter = floD.mesh(dim, 'models/rae_25.msh', pars, ['field', 'wing', 'symmetry', 'downstream'])
+ # morpher = floD.morpher(msh, dim, 'wing')
+ tms["msh"].stop()
+
+ # set the problem and add fluid, initial/boundary conditions
+ tms['pre'].start()
+ # Replce tp = 'teTip' by tp = 'wakeTip' for oneraM6.geo
+ pbl, _, _, bnd, blw = floD.problem(msh, dim, alpha, 0., M_inf, S_ref, c_ref, 0.25, 0., 0., 'wing', vsc = True)
+ tms['pre'].stop()
+
+ config = {'nDim': dim, 'Sections':[-0.1, -0.05, 0.0, 0.05, 0.1], 'span':spn, 'Sym':[-0.1, 0.1], 'rbftype': 'linear', 'smoothing': 1e-6, 'degree': 0, 'neighbors': 10}
+
+ # solve problem
+ iSolver = floD.newton(pbl)
+
+ # VII
+ iSolverAPI = vInterpol.Interpolator(iSolver, blw, imsh, vmsh, config)
+ vSolver = viscU.initBL(Re, M_inf, CFL0, len(config['Sections']), spn, 2)
+ coupler = viscC.Coupler(iSolverAPI, vSolver, _maxIters=1)
+
+ coupler.Run()
+
+ U.append(iSolverAPI.URani)
+
+ for i in range(len(U)):
+ plt.plot((U[i]-U[0])/U[0])
+ plt.show()
+
+ plt.show()
+ vSolution = viscU.GetSolution(0, vSolver)
+ viscU.WriteFile(vSolution, vSolver.GetRe())
+
+ plt.plot(vSolution['x/c'], vSolution['theta'])
+ plt.show()
+
+ # Write corrected Cp
+ iSolver.save(mshWriter)
+ dartU.writeSlices(msh.name, config['Sections'], 2)
+
+ cpCorr = [np.zeros((0,0)) for _ in range(len(config['Sections']))]
+ for i in range(len(config['Sections'])):
+ cpCorr[i] = tboxU.read('rae_25_slice_'+str(i)+'.dat')
+
+ # Plot cp
+
+ for i in range(len(config['Sections'])):
+ plt.plot(cpCorr[i][:,3], cpCorr[i][:,4], '-')
+ plt.plot(cpInv[i][:,3], cpInv[i][:,4], '--')
+ plt.xlabel('$x/c$')
+ plt.ylabel('$Cp$')
+ plt.title('$y/b$ = {:.3f}'.format(config['Sections'][i]))
+ plt.gca().invert_yaxis()
+ plt.legend(['Corrected', 'Inviscid'])
+ plt.show()
+
+ # display timers
+ tms['total'].stop()
+ print(ccolors.ANSI_BLUE + 'PyTiming...' + ccolors.ANSI_RESET)
+ print('CPU statistics')
+ print(tms)
+ print(coupler.tms)
+
+ # visualize solution
+ iSolverAPI.GetCp()
+ iSolverAPI.iSolver.save(gmshWriter)
+ floD.initViewer(pbl)
+
+ # check results
+ print(ccolors.ANSI_BLUE + 'PyTesting...' + ccolors.ANSI_RESET)
+ tests = CTests()
+ if alpha == 3*math.pi/180 and M_inf == 0.3 and spn == 1.0:
+ tests.add(CTest('iteration count', solver.nIt, 3, 1, forceabs=True))
+ tests.add(CTest('CL', solver.Cl, 0.135, 5e-2))
+ tests.add(CTest('CD', solver.Cd, 0.0062, 1e-2)) # Tranair (NF=0.0062, FF=0.0030), Panair 0.0035
+ tests.add(CTest('CS', solver.Cs, 0.0121, 5e-2))
+ tests.add(CTest('CM', solver.Cm, -0.0278, 1e-2))
+ else:
+ raise Exception('Test not defined for this flow')
+ tests.run()
+
+ # eof
+ print('')
+
+if __name__ == "__main__":
+ main()
--
GitLab