diff --git a/dart/pyVII/vCoupler.py b/dart/pyVII/vCoupler.py
index 05c672533a6f1455e5cdea974fcfe5fdf38e11db..59e5ba09cfec3cf80d564e2b85f0fe17c700c55a 100644
--- a/dart/pyVII/vCoupler.py
+++ b/dart/pyVII/vCoupler.py
@@ -36,7 +36,7 @@ class Coupler:
self.iSolver=iSolver
self.vSolver=vSolver
- self.maxCouplIter = 150
+ self.maxCouplIter = 1
self.couplTol = 1e-4
self.tms=fwk.Timers()
@@ -55,11 +55,14 @@ class Coupler:
# Initilize meshes in c++ objects
+ self.iSolver.printOn = False
+ self.vSolver.printOn = False
+
dragIter = np.zeros((0,2))
self.nElmUpperPrev = 0
couplIter = 0
- cdPrev = 0
+ cdPrev = 0.0
"""print('- AoA: {:<2.2f} Mach: {:<1.1f} Re: {:<2.1f}e6'
.format(self.iSolver.pbl.alpha*180/math.pi, self.iSolver.pbl.M_inf, self.vSolver.Re/1e6))
print('{:>5s}| {:>10s} {:>10s} | {:>8s} {:>6s} {:>8s}'.format('It', 'Cl', 'Cd', 'xtrT', 'xtrB', 'Error'))"""
diff --git a/dart/pyVII/vCoupler2.py b/dart/pyVII/vCoupler2.py
index 355b106d857fa3a5025ecceb87959d429dd3e2e7..b388662eed97cffd5ae37d52fd459cc2143a1206 100644
--- a/dart/pyVII/vCoupler2.py
+++ b/dart/pyVII/vCoupler2.py
@@ -34,11 +34,11 @@ class Coupler:
self.iSolverAPI = _iSolverAPI
self.vSolver = vSolver
- self.maxCouplIter = 25
+ self.maxCouplIter = 150
self.couplTol = 1e-4
self.tms=fwk.Timers()
- self.resetInv = True
+ self.resetInv = False
pass
def Run(self):
diff --git a/dart/pyVII/vInterpolator.py b/dart/pyVII/vInterpolator.py
index abe70e5563f8fdf6c82ef81cbd34b72ce08bb953..22a289f474638e401a66712ca6ac8c4bed803be7 100644
--- a/dart/pyVII/vInterpolator.py
+++ b/dart/pyVII/vInterpolator.py
@@ -36,6 +36,7 @@ class Interpolator:
self.config = _cfg
self.xxPrev = [[np.zeros(0) for iReg in range(3)] for iSec in range(len(self.viscStruct))]
+ self.uePrev = [[np.zeros(len(self.viscStruct[iSec][iReg][:,0])) for iReg in range(2)] for iSec in range(len(self.viscStruct))]
self.DeltaStarPrev = [[np.zeros(len(self.viscStruct[iSec][iReg][:,0])) for iReg in range(2)] for iSec in range(len(self.viscStruct))]
self.tms = fwk.Timers()
@@ -60,7 +61,9 @@ class Interpolator:
# Wake
nElm = len(self.bndWake.nodes)
- U_wk, M_wk, Rho_wk = np.zeros((nElm, 3)), np.zeros(nElm), np.zeros(nElm)
+ U_wk =np.zeros((nElm, 3))
+ M_wk = np.zeros(nElm)
+ Rho_wk = np.zeros(nElm)
for iNode, n in enumerate(self.bndWake.nodes):
U_wk[iNode,0] = self.iSolver.U[n.row][0]
U_wk[iNode,1] = self.iSolver.U[n.row][1]
@@ -99,36 +102,6 @@ class Interpolator:
UyUp, UyLw = self.__Remesh(Uy[iSec][0], self.stgPt)
UzUp, UzLw = self.__Remesh(Uz[iSec][0], self.stgPt)
- """plt.plot(xUp, UxUp)
- plt.plot(xLw, UxLw)
- plt.plot(self.viscStruct[iSec][1][:,0], Ux[iSec][1])
- plt.plot(self.invStruct[0][:,0], U_wg[:,0], 'x')
- plt.title('Ux = f(x)')
- plt.show()
-
- plt.plot(xUp, UyUp)
- plt.plot(xLw, UyLw)
- plt.plot(self.viscStruct[iSec][1][:,0], Uy[iSec][1])
- plt.plot(self.invStruct[0][:,0], U_wg[:,1], 'x')
- plt.title('Uy = f(x) smooting{:.2e}'.format(self.config['smoothing']))
- plt.show()
-
- plt.plot(xUp, UzUp)
- plt.plot(xLw, UzLw)
- plt.plot(self.viscStruct[iSec][1][:,0], Uz[iSec][1])
- plt.plot(self.invStruct[0][:,0], U_wg[:,2], 'x')
- plt.title('Uz = f(x)')
- plt.show()"""
-
- """plt.plot(xUp, UyUp, 'x')
- plt.plot(xLw, UyLw, 'o')
- #plt.plot(self.viscStruct[iSec][1][:,0], Uy[iSec][1])
- plt.plot(self.invStruct[0][:,0], U_wg[:,1], 'x')
- plt.xlim([0.4, 1.2])
- plt.ylim([-0.3, 0.2])
- plt.title('Uy = f(x) smooting{:.2e}'.format(self.config['smoothing']))
- plt.show()"""
-
vSolver.SetVelocities(iSec, 0, UxUp, UyUp, UzUp)
vSolver.SetVelocities(iSec, 1, UxLw, UyLw, UzLw)
@@ -155,6 +128,11 @@ class Interpolator:
vSolver.SetxxExt(iSec, 0, xxPrevUp)
vSolver.SetxxExt(iSec, 1, xxPrevLw)
vSolver.SetxxExt(iSec, 2, self.xxPrev[iSec][1])
+
+ uePrevUp, uePrevLw = self.__Remesh(self.uePrev[iSec][0], self.stgPt)
+ vSolver.SetUeOld(iSec, 0, uePrevUp)
+ vSolver.SetUeOld(iSec, 1, uePrevLw)
+ vSolver.SetUeOld(iSec, 2, self.uePrev[iSec][1])
@@ -171,6 +149,9 @@ class Interpolator:
self.DeltaStarPrev[iSec][1] = np.asarray(vSolver.ExtractDeltaStar(iSec, 2))
self.xxPrev[iSec][0] = np.concatenate((np.asarray(vSolver.Extractxx(iSec, 0))[::-1], np.asarray(vSolver.Extractxx(iSec, 1))[1:]), axis=None)
self.xxPrev[iSec][1] = np.asarray(vSolver.Extractxx(iSec, 2))
+ self.uePrev[iSec][0] = np.concatenate((np.asarray(vSolver.ExtractUe(iSec, 0))[::-1], np.asarray(vSolver.ExtractUe(iSec, 1))[1:]), axis=None)
+ self.uePrev[iSec][1] = np.asarray(vSolver.ExtractUe(iSec, 2))
+
# Compute cell centers
viscCgWing = np.zeros((len(self.viscStruct[0][0][:,0]) - 1, 3))
@@ -205,14 +186,6 @@ class Interpolator:
if self.config['nDim'] == 3:
invCgWing[:,[1,2]] = invCgWing[:,[2,1]]
invCgWake[:,[1,2]] = invCgWake[:,[2,1]]
-
- """plt.plot(self.invStruct[0][:,0], self.invStruct[0][:,1], 'x')
- plt.plot(invCgWing[:, 0], invCgWing[:,1], 'o')
- plt.show()
-
- plt.plot(self.viscStruct[0][0][:,0], self.viscStruct[0][0][:,1], 'x')
- plt.plot(viscCgWing[:,0], viscCgWing[:,1], 'o')
- plt.show()"""
self.tms['Interpolation'].start()
invCgUp = invCgWing[invCgWing[:,1]>0, :]
diff --git a/dart/src/wBLRegion.cpp b/dart/src/wBLRegion.cpp
index c1ed9ac5717d3cb311712ce3fb061a8eb8917baf..03a01b104e79926bbb0f3136e4a7153b53cbe7e3 100644
--- a/dart/src/wBLRegion.cpp
+++ b/dart/src/wBLRegion.cpp
@@ -8,12 +8,21 @@
using namespace dart;
-BLRegion::BLRegion()
+BLRegion::BLRegion(double _xtrF)
{
+ if (_xtrF < 0.0 && _xtrF != -1.0)
+ {
+ std::cout << "Unrecognized transition location " << _xtrF << std::endl;
+ throw;
+ }
+ xtrF = _xtrF;
closSolver = new Closures();
mesh = new ViscMesh();
invBnd = new InvBnd();
+ transMarker = mesh->GetnMarkers();
+ xtr = 1.0;
+
} // end BLRegion
BLRegion::~BLRegion()
@@ -143,7 +152,7 @@ void BLRegion::ComputeBlwVel()
}
-void BLRegion::PrintSolution(size_t iPoint)
+void BLRegion::PrintSolution(size_t iPoint) const
{
std::cout << "\nPoint " << iPoint << " \n"
<< std::scientific;
diff --git a/dart/src/wBLRegion.h b/dart/src/wBLRegion.h
index b81d15e6edc87537a5ef19a56987b59e1998758c..0e2cc5eb37ac3329999349ff4137e6cbca9f17ea 100644
--- a/dart/src/wBLRegion.h
+++ b/dart/src/wBLRegion.h
@@ -54,7 +54,7 @@ public:
double xtrF;
double xtr;
- BLRegion();
+ BLRegion(double _xtrF = -1.0);
~BLRegion();
/* Distribute BC */
@@ -67,11 +67,11 @@ public:
void SetWakeBC(double Re, std::vector<double> UpperCond, std::vector<double> LowerCond);
/* General access */
- unsigned int GetnVar() {return nVar;};
- double GetnCrit() {return nCrit;};
+ unsigned int GetnVar() const {return nVar;};
+ double GetnCrit() const {return nCrit;};
/* Others */
- void PrintSolution(size_t iPoint);
+ void PrintSolution(size_t iPoint) const;
void ComputeBlwVel();
};
diff --git a/dart/src/wClosures.cpp b/dart/src/wClosures.cpp
index 6aa84a0c1b343387767ce1e278e13058c0454db9..6497f22565e06567e18014af70b1c18ba7e07d76 100644
--- a/dart/src/wClosures.cpp
+++ b/dart/src/wClosures.cpp
@@ -14,7 +14,7 @@ Closures::~Closures()
// std::cout << "~Closures()\n";
} // end ~Closures
-std::vector<double> Closures::LaminarClosures(double theta, double H, double Ret, double Me, std::string name)
+std::vector<double> Closures::LaminarClosures(double theta, double H, double Ret, double Me, std::string name) const
{
H = std::max(H, 1.00005);
@@ -102,10 +102,19 @@ std::vector<double> Closures::LaminarClosures(double theta, double H, double Ret
ClosureVars[5] = delta;
ClosureVars[6] = Cteq;
ClosureVars[7] = Us;
+
+ for (size_t i=0; i<ClosureVars.size(); ++i)
+ {
+ if (std::isnan(ClosureVars[i]))
+ {
+ std::cout << "lam : nan detected in closureVars at position " << i << std::endl;
+ throw;
+ }
+ }
return ClosureVars;
}
-std::vector<double> Closures::TurbulentClosures(double theta, double H, double Ct, double Ret, double Me, std::string name)
+std::vector<double> Closures::TurbulentClosures(double theta, double H, double Ct, double Ret, double Me, std::string name) const
{
H = std::max(H, 1.00005);
@@ -154,6 +163,11 @@ std::vector<double> Closures::TurbulentClosures(double theta, double H, double C
double logRt = std::log(Ret/Fc);
+ if (std::isnan(logRt))
+ {
+ std::cout << "logRt is nan " << std::endl;
+ throw;
+ }
logRt = std::max(logRt, 3.0);
double arg = std::max(-1.33*Hk, -20.0);
@@ -228,6 +242,14 @@ std::vector<double> Closures::TurbulentClosures(double theta, double H, double C
ClosureVars[5] = delta;
ClosureVars[6] = Cteq;
ClosureVars[7] = Us;
- return ClosureVars;
+ for (size_t i=0; i<ClosureVars.size(); ++i)
+ {
+ if (std::isnan(ClosureVars[i]))
+ {
+ std::cout << "turb : nan detected in closureVars at position " << i << std::endl;
+ throw;
+ }
+ }
+ return ClosureVars;
}
\ No newline at end of file
diff --git a/dart/src/wClosures.h b/dart/src/wClosures.h
index 2946cf76748751ff5b246f89a425c36eca6c2627..e86fc183f05d9abc9aee05b67934d17f123c17ad 100644
--- a/dart/src/wClosures.h
+++ b/dart/src/wClosures.h
@@ -13,8 +13,8 @@ class DART_API Closures
public:
Closures();
~Closures();
- std::vector<double> LaminarClosures(double theta, double H, double Ret, double Me, std::string name);
- std::vector<double> TurbulentClosures(double theta, double H, double Ct, double Ret, double Me, std::string name);
+ std::vector<double> LaminarClosures(double theta, double H, double Ret, double Me, std::string name) const;
+ std::vector<double> TurbulentClosures(double theta, double H, double Ct, double Ret, double Me, std::string name) const;
};
} // namespace dart
#endif //WClosures_H
\ No newline at end of file
diff --git a/dart/src/wInvBnd.h b/dart/src/wInvBnd.h
index 0fa1fa85ef2ef7c3424fd06a456b6bd6e5ff0be2..e1aa68d25770072a0a13aeeaa9040de99ef6a890 100644
--- a/dart/src/wInvBnd.h
+++ b/dart/src/wInvBnd.h
@@ -15,11 +15,12 @@ class DART_API InvBnd
std::vector<double> Me,
Vt,
Rhoe,
- DeltaStar;
+ DeltaStar,
+ VtOld;
public:
- std::vector<double> BlowingVel;
+ //std::vector<double> BlowingVel;
InvBnd();
~InvBnd();
@@ -28,11 +29,13 @@ class DART_API InvBnd
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];};
void SetMe(std::vector<double> _Me) {Me.resize(_Me.size(), 0.0); Me = _Me;};
void SetVt(std::vector<double> _Vt) {Vt.resize(_Vt.size(), 0.0); Vt = _Vt;};
void SetRhoe(std::vector<double> _Rhoe) {Rhoe.resize(_Rhoe.size(), 0.0); Rhoe = _Rhoe;};
void SetDeltaStar(std::vector<double> _DeltaStar) {DeltaStar.resize(_DeltaStar.size(), 0.0); DeltaStar = _DeltaStar;};
+ void SetVtOld(std::vector<double> _VtOld) {VtOld.resize(_VtOld.size(), 0.0); VtOld = _VtOld;};
};
} // namespace dart
diff --git a/dart/src/wTimeSolver.cpp b/dart/src/wTimeSolver.cpp
index d5f16791fa70f2248607b445c133f9fb38456145..265f5c0961c52f3d4a2868f53014d798e2eb5a55 100644
--- a/dart/src/wTimeSolver.cpp
+++ b/dart/src/wTimeSolver.cpp
@@ -58,10 +58,6 @@ void TimeSolver::InitialCondition(size_t iPoint, BLRegion *bl)
bl->U[iPoint * nVar + 4] = 0.03;
} // End if
- /*if (bl->U[iPoint*nVar+3]<0)
- {
- std::cout << "dart::TimeSolver negative velocity encountered at point " << iPoint << std::endl;
- }*/
} // End InitialCondition
int TimeSolver::Integration(size_t iPoint, BLRegion *bl)
@@ -103,6 +99,11 @@ int TimeSolver::Integration(size_t iPoint, BLRegion *bl)
while (innerIters < maxIt)
{
+ if (bl->name=="upper" && iPoint == 1)
+ {
+ std::cout << "it " << innerIters << " res " << std::log10(normSysRes/normSysRes0) << std::endl;
+ bl->PrintSolution(iPoint);
+ }
if (innerIters % itFrozenJac == 0)
{
@@ -129,7 +130,8 @@ int TimeSolver::Integration(size_t iPoint, BLRegion *bl)
IMatrix = Matrix<double, 5, 5>::Identity() / timeStep;
normSysRes = (SysRes + slnIncr / timeStep).norm();
itFrozenJac = 1;
-
+ std::cout << "normSysRes is nan " << std::endl;
+ throw;
}
if (normSysRes <= tol * normSysRes0)
@@ -152,7 +154,7 @@ int TimeSolver::Integration(size_t iPoint, BLRegion *bl)
return innerIters;
} // End Integration
-double TimeSolver::SetTimeStep(double CFL, double dx, double invVel)
+double TimeSolver::SetTimeStep(double CFL, double dx, double invVel) const
{
/* Speed of sound. */
double gradU2 = (invVel * invVel);
@@ -165,7 +167,7 @@ double TimeSolver::SetTimeStep(double CFL, double dx, double invVel)
return CFL * dx / advVel;
}
-double TimeSolver::ComputeSoundSpeed(double gradU2)
+double TimeSolver::ComputeSoundSpeed(double gradU2) const
{
return sqrt(1 / (Minf * Minf) + 0.5 * (gamma - 1) * (1 - gradU2));
}
diff --git a/dart/src/wTimeSolver.h b/dart/src/wTimeSolver.h
index e45d68906cfc54f5f77a7c1a153e65fbd487e5a6..5866372da75083e336a1aa03803866c8f2d508cf 100644
--- a/dart/src/wTimeSolver.h
+++ b/dart/src/wTimeSolver.h
@@ -38,8 +38,8 @@ public:
void SetinitSln(bool _initSln);
private:
- double SetTimeStep(double CFL, double dx, double invVel);
- double ComputeSoundSpeed(double gradU2);
+ double SetTimeStep(double CFL, double dx, double invVel) const;
+ double ComputeSoundSpeed(double gradU2) const;
void ResetSolution(size_t iPoint, BLRegion *bl, std::vector<double> Sln0, bool usePrevPoint = false);
};
} // namespace dart
diff --git a/dart/src/wViscFluxes.cpp b/dart/src/wViscFluxes.cpp
index 60854e3de990330597fcd7e597e2282fc16d1293..8982b4f96361a786ee07826c4907161870210b9f 100644
--- a/dart/src/wViscFluxes.cpp
+++ b/dart/src/wViscFluxes.cpp
@@ -111,6 +111,7 @@ Vector<double, 5> ViscFluxes::BLlaws(size_t iPoint, BLRegion *bl, std::vector<do
{
std::cout << "Wrong regime\n"
<< std::endl;
+ throw;
} // End else
/* Gradients computation */
@@ -123,17 +124,19 @@ Vector<double, 5> ViscFluxes::BLlaws(size_t iPoint, BLRegion *bl, std::vector<do
double dHstar_dx = (bl->Hstar[iPoint] - bl->Hstar[iPoint - 1]) / dx;
double dueExt_dx = (bl->invBnd->GetVt(iPoint) - bl->invBnd->GetVt(iPoint - 1)) / dxExt;
+ //std::cout << bl->invBnd->GetVtOld(iPoint) << " " << bl->invBnd->GetVt(iPoint) << std::endl;
double ddeltaStarExt_dx = (bl->invBnd->GetDeltaStar(iPoint) - bl->invBnd->GetDeltaStar(iPoint - 1)) / dxExt;
double c = 4 / (M_PI * dx);
double cExt = 4 / (M_PI * dxExt);
/* Amplification routine */
- double dN_dx = 0;
- double Ax = 0;
+ double dN_dx = 0.0;
+ double Ax = 0.0;
- if (bl->xtrF != -1 && bl->Regime[iPoint] == 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;
@@ -222,7 +225,7 @@ Vector<double, 5> ViscFluxes::BLlaws(size_t iPoint, BLRegion *bl, std::vector<do
timeMatrix(4, 4) = 1.0;
}
- /*Vector<double, 5> result = timeMatrix.partialPivLu().solve(spaceVector);
+ Vector<double, 5> result = timeMatrix.partialPivLu().solve(spaceVector);
if (std::isnan(result.norm()))
{
std::cout << "Point " << iPoint << " det " << timeMatrix.determinant() << std::endl;
@@ -231,12 +234,12 @@ Vector<double, 5> ViscFluxes::BLlaws(size_t iPoint, BLRegion *bl, std::vector<do
std::cout << dT_dx << " " << (2 + u[1] - Mea * Mea) * (u[0] / u[3]) * due_dx << " " << " cf " << - Cfa / 2 << std::endl;
std::cout << " u[0] " << u[0] << " u[1] " << u[1] << " u[2] " << u[2] << " u[3] " << u[3] << " u[4] " << u[4] << " Usa " << Usa << std::endl;
throw;
- }*/
+ }
return timeMatrix.partialPivLu().solve(spaceVector);
}
-double ViscFluxes::AmplificationRoutine(double Hk, double Ret, double theta)
+double ViscFluxes::AmplificationRoutine(double Hk, double Ret, double theta) const
{
double Ax = 0.0;
diff --git a/dart/src/wViscFluxes.h b/dart/src/wViscFluxes.h
index 403c388e8ce4fd9d6fc53070dcfe40233a558ef8..1c9867f34e29f1b40c617b0a5d4d9ae8f36ccfc1 100644
--- a/dart/src/wViscFluxes.h
+++ b/dart/src/wViscFluxes.h
@@ -23,7 +23,7 @@ public:
private:
Eigen::Vector<double, 5> BLlaws(size_t iPoint, BLRegion *bl, std::vector<double> u);
- double AmplificationRoutine(double hk, double ret, double theta);
+ double AmplificationRoutine(double hk, double ret, double theta) const;
};
} // namespace dart
#endif //WViscFluxes_H
\ No newline at end of file
diff --git a/dart/src/wViscSolver.cpp b/dart/src/wViscSolver.cpp
index 1560e34e5751e31007d948bf287ed5ffa5e03e71..24846a42f97e595c21b3f7629ce714a620e692e4 100644
--- a/dart/src/wViscSolver.cpp
+++ b/dart/src/wViscSolver.cpp
@@ -38,6 +38,9 @@ ViscSolver::ViscSolver(double _Re, double _Minf, double _CFL0, unsigned int nSec
Re = _Re; /* Freestream Reynolds number */
CFL0 = _CFL0;
+ Cdt = 0.0;
+ Cdf = 0.0;
+ Cdp = 0.0;
solverExitCode = 0;
@@ -103,7 +106,7 @@ int ViscSolver::Run(unsigned int couplIter)
for (size_t iSec=0; iSec<Sections.size(); ++iSec)
{
- std::cout << "\n - Sec " << iSec << " stagPtMvmt " << stagPtMvmt[iSec];
+ std::cout << "\n - Sec " << iSec << " stagPtMvmt " << stagPtMvmt[iSec] << Sections[iSec][0]->mesh->GetDiffnElm();
/* Check for stagnation point movement */
@@ -161,6 +164,13 @@ int ViscSolver::Run(unsigned int couplIter)
{
std::vector<double> DeltaStarZeros(Sections[iSec][iRegion]->mesh->GetnMarkers(), 0.0);
Sections[iSec][iRegion]->invBnd->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]->invBnd->GetVt(iPoint);
+ }
+ Sections[iSec][iRegion]->invBnd->SetVtOld(ueOld);
}
if (printOn){std::cout << "restart. ";}
@@ -212,9 +222,8 @@ int ViscSolver::Run(unsigned int couplIter)
tSolver->InitialCondition(iPoint, Sections[iSec][iRegion]);
/* Time integration */
- if (iRegion == 1 && iPoint == 1)
+ /* if (iRegion == 0 && iPoint == 1)
{
- std::scientific;
std::cout << " iSec " << iSec
<< " iPoint " << iPoint
<< " Loc " << Sections[iSec][iRegion]->mesh->GetLoc(iPoint)
@@ -224,7 +233,7 @@ int ViscSolver::Run(unsigned int couplIter)
<< " DeltaStarExt " << Sections[iSec][iRegion]->invBnd->GetDeltaStar(iPoint)
<< std::endl;
Sections[iSec][iRegion]->PrintSolution(iPoint);
- }
+ } */
pointExitCode = tSolver->Integration(iPoint, Sections[iSec][iRegion]);
@@ -246,8 +255,8 @@ int ViscSolver::Run(unsigned int couplIter)
if (pointExitCode != 0)
{
if (printOn){std::cout << "Point " << iPoint << ": Convergence terminated with exitcode " << pointExitCode << std::endl;}
- stagPtMvmt[iSec] = true;
solverExitCode = -1; // Convergence failed
+ std::cout << " NC " << iPoint << " ";
}
@@ -331,6 +340,7 @@ void ViscSolver::AverageTransition(size_t iPoint, BLRegion *bl, int forced)
}
/* Compute transition location */
+ //std::cout << "Trans " << bl->mesh->Getx(iPoint) << std::endl;
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 */
@@ -399,15 +409,15 @@ void ViscSolver::ComputeDrag(std::vector<BLRegion *> const &bl)
/* Normalize friction and dissipation coefficients. */
- std::vector<double> frictioncoeff_upper;
- for (size_t k = 0; k < bl[0]->mesh->GetnMarkers(); ++k)
+ std::vector<double> frictioncoeff_upper(bl[0]->mesh->GetnMarkers(), 0.0);
+ for (size_t k = 0; k < frictioncoeff_upper.size(); ++k)
{
- frictioncoeff_upper.push_back(bl[0]->U[k * nVar + 3] * bl[0]->U[k * nVar + 3] * bl[0]->Cf[k]);
+ frictioncoeff_upper[k] = bl[0]->invBnd->GetVt(k) * bl[0]->invBnd->GetVt(k) * bl[0]->Cf[k];
}
- std::vector<double> frictioncoeff_lower;
- for (size_t k = 0; k < bl[1]->mesh->GetnMarkers(); ++k)
+ std::vector<double> frictioncoeff_lower(bl[1]->mesh->GetnMarkers(), 0.0);
+ for (size_t k = 0; k < frictioncoeff_lower.size(); ++k)
{
- frictioncoeff_lower.push_back(bl[1]->U[k * nVar + 3] * bl[1]->U[k * nVar + 3] * bl[1]->Cf[k]);
+ frictioncoeff_lower[k] = bl[1]->invBnd->GetVt(k) * bl[1]->invBnd->GetVt(k) * bl[1]->Cf[k];
}
/* Compute friction drag coefficient. */
@@ -442,17 +452,23 @@ void ViscSolver::SetGlobMesh(size_t iSec, size_t iRegion, std::vector<double> _x
void ViscSolver::SetVelocities(size_t iSec, size_t iRegion, std::vector<double> vx, std::vector<double> vy, std::vector<double> vz)
{
std::vector<double> Vt(vx.size(), 0.0);
- double alpha=0;
+ double alpha=0.0;
for (size_t iPoint=0; iPoint<vx.size()-1; ++iPoint)
{
- alpha = std::atan2(Sections[iSec][iRegion]->mesh->Gety(iPoint+1) - Sections[iSec][iRegion]->mesh->Gety(iPoint),
- Sections[iSec][iRegion]->mesh->Getx(iPoint+1) - Sections[iSec][iRegion]->mesh->Getx(iPoint));
+ alpha = std::atan2((Sections[iSec][iRegion]->mesh->Gety(iPoint+1) - Sections[iSec][iRegion]->mesh->Gety(iPoint)),
+ (Sections[iSec][iRegion]->mesh->Getx(iPoint+1) - Sections[iSec][iRegion]->mesh->Getx(iPoint)));
Vt[iPoint] = vx[iPoint]*std::cos(alpha) + vy[iPoint]*std::sin(alpha);
+ Vt[iPoint+1] = vx[iPoint+1]*std::cos(alpha) + vy[iPoint+1]*std::sin(alpha);
+ /*if (Vt[iPoint] < 0)
+ {
+ std::cout << "Negative velocity at point " << iPoint << std::endl;
+ //throw;
+ }*/
// std::cout << iSec << iRegion << iPoint << "vx " << vx[iPoint] << "vy" << vy[iPoint] << "vz" << vz[iPoint] << "vt" << Vt[iPoint] << std::endl;
}
- Vt.back() = vx.back()*cos(alpha) + vy.back()*std::sin(alpha);
+ //Vt.back() = vx.back()*cos(alpha) + vy.back()*std::sin(alpha);
Sections[iSec][iRegion]->invBnd->SetVt(Vt);
}
@@ -479,10 +495,10 @@ void ViscSolver::SetDeltaStarExt(size_t iSec, size_t iRegion, std::vector<double
std::vector<double> ViscSolver::ExtractxCoord(size_t iSec, size_t iRegion) const
{
- std::vector<double> xCoord;
- for (size_t iPoint=0; iPoint<Sections[iSec][iRegion]->mesh->GetnMarkers(); ++iPoint)
+ std::vector<double> xCoord(Sections[iSec][iRegion]->mesh->GetnMarkers(), 0.0);
+ for (size_t iPoint=0; iPoint<xCoord.size(); ++iPoint)
{
- xCoord.push_back(Sections[iSec][iRegion]->mesh->Getx(iPoint));
+ xCoord[iPoint] = Sections[iSec][iRegion]->mesh->Getx(iPoint);
}
return xCoord;
}
@@ -493,8 +509,8 @@ std::vector<double> ViscSolver::ExtractBlowingVel(size_t iSec, size_t iRegion) c
}
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<Sections[iSec][iRegion]->mesh->GetnMarkers(); ++iPoint)
+ std::vector<double> xx(Sections[iSec][iRegion]->mesh->GetnMarkers(),0.0);
+ for (size_t iPoint=0; iPoint<xx.size(); ++iPoint)
{
xx[iPoint] = Sections[iSec][iRegion]->mesh->GetLoc(iPoint);
}
@@ -505,7 +521,23 @@ std::vector<double> ViscSolver::ExtractDeltaStar(size_t iSec, size_t iRegion) co
return Sections[iSec][iRegion]->DeltaStar;
}
-void ViscSolver::PrintBanner()
+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)
+{
+ Sections[iSec][iRegion]->invBnd->SetVtOld(_ue);
+}
+
+
+void ViscSolver::PrintBanner() const
{
std::cout << "\n\n"
<< std::endl;
diff --git a/dart/src/wViscSolver.h b/dart/src/wViscSolver.h
index bf0bd953683c1665411e5ee99b59f962166b3e83..eafbb9ad1050149c3b1129c09bb5c65a6ebe2278 100644
--- a/dart/src/wViscSolver.h
+++ b/dart/src/wViscSolver.h
@@ -21,7 +21,7 @@ private:
double Re,
CFL0;
- int solverExitCode;
+ int solverExitCode = 0;
std::vector<bool> stagPtMvmt;
@@ -45,6 +45,11 @@ public:
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;
+
+ void SetUeOld(size_t iSec, size_t iRegion, std::vector<double> _ue);
+
+
double Getxtr(size_t iSec, size_t iRegion) const {return Sections[iSec][iRegion]->xtr;};
private:
@@ -52,7 +57,7 @@ private:
void AverageTransition(size_t iPoint, BLRegion *bl, int forced);
void ComputeDrag(std::vector<BLRegion *> const &bl);
void ComputeBlwVel();
- void PrintBanner();
+ void PrintBanner() const;
};
} // namespace dart
#endif //WVISCSOLVER_H
\ No newline at end of file
diff --git a/dart/tests/bli_lowLift.py b/dart/tests/bli_lowLift.py
index 527fe3c424e4c64926cedd9f933fb1fede3dce51..908c5edbaf285ac22b547e905cb69716e0a7a3ef 100644
--- a/dart/tests/bli_lowLift.py
+++ b/dart/tests/bli_lowLift.py
@@ -60,7 +60,7 @@ def main():
# define flow variables
Re = 1e7
- alpha = 5*math.pi/180
+ alpha = 12*math.pi/180
#alphaSeq = np.arange(-5, 5.5, 0.5)
M_inf = 0.
@@ -80,7 +80,7 @@ def main():
# 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}
+ pars = {'xLgt' : 5, 'yLgt' : 5, 'msF' : 1, 'msTe' : 0.01, 'msLe' : 0.0001}
msh, gmshWriter = floD.mesh(dim, 'models/n0012.geo', pars, ['field', 'airfoil', 'downstream'])
tms['msh'].stop()
@@ -97,25 +97,19 @@ def main():
isolver = floD.newton(pbl)
vsolver = floD.initBL(Re, M_inf, CFL0, nSections)
config = {'nDim': dim, 'nSections':nSections, 'nPoints':[600, 50], 'eps':0.02, 'rbftype': 'cubic', 'smoothing': 0.0, 'degree': 1, 'neighbors': 10}
- iSolverAPI = Interpol.Interpolator(isolver, blw[0], blw[1], config)
- coupler = viscC.Coupler(iSolverAPI, vsolver)
if rbf:
- fig, axs = plt.subplots(2)
- for i in range(3):
- aeroCoeff = coupler.Run()
- isolver.reset()
- axs[0].plot(aeroCoeff[:,0], 'x-')
- axs[1].plot(aeroCoeff[:,1], 'x-')
+ iSolverAPI = Interpol.Interpolator(isolver, blw[0], blw[1], config)
+ coupler = viscC.Coupler(iSolverAPI, vsolver)
+ aeroCoeff = coupler.Run()
+ isolver.reset()
+ axs[0].plot(aeroCoeff[:,0], 'x-')
+ axs[1].plot(aeroCoeff[:,1], 'x-')
axs[0].set(xlabel='iters', ylabel='$c_l$')
axs[1].set(xlabel='iters', ylabel='$c_d$')
plt.show()
else:
- msh, gmshWriter = floD.mesh(dim, 'models/n0012.geo', pars, ['field', 'airfoil', 'downstream'])
- pbl, _, _, bnd, blw = floD.problem(msh, dim, alpha, 0., M_inf, c_ref, c_ref, 0.25, 0., 0., 'airfoil', te = 'te', vsc = True, dbc=True)
- isolver = floD.newton(pbl)
- vsolver = floD.initBL(Re, M_inf, CFL0, nSections)
coupler = floVC.Coupler(isolver, vsolver)
coupler.CreateProblem(blw[0], blw[1])
fig, axs = plt.subplots(2)