Skip to content
Snippets Groups Projects
Commit 62568607 authored by Paul Dechamps's avatar Paul Dechamps
Browse files

Changed Y Z directions management

parent a366558a
No related branches found
No related tags found
No related merge requests found
Pipeline #7349 failed
Expand all Hide whitespace changes
Inline Side-by-side
vii/models/n0012_3_visc.geo
+ 2
− 2
View file @ 62568607
......@@ -355,11 +355,11 @@ nElmChrodwise_wake = 40;
Transfinite Curve {222, 221} = nElmChrodwise_wake Using Progression 1.05;
Transfinite Surface {51};
nChordwise_midWing = 40;
nChordwise_midWing = 80;
//+
Transfinite Curve {8, 11, 5, 2, 122} = nChordwise_midWing Using Progression 1;
//+
nChordwise_Le = 40;
nChordwise_Le = 70;
Transfinite Curve {-9, 10, -3, 4} = nChordwise_Le Using Progression 1.07;
//+
......
vii/models/n0012_3_visc.msh 0 → 100644
+ 33366
− 0
View file @ 62568607
This diff is collapsed.
vii/pyVII/vCoupler.py
+ 11
− 0
View file @ 62568607
......@@ -247,6 +247,17 @@ class Coupler:
"""plt.plot(fMeshDict['globCoord'][:fMeshDict['bNodes'][1]], fMeshDict['vx'][:fMeshDict['bNodes'][1]])
plt.plot(fMeshDict['globCoord'][fMeshDict['bNodes'][1]:fMeshDict['bNodes'][2]], fMeshDict['vx'][fMeshDict['bNodes'][1]:fMeshDict['bNodes'][2]])
plt.show()"""
plt.plot(fMeshDict['globCoord'][:fMeshDict['bNodes'][1]], fMeshDict['vx'][:fMeshDict['bNodes'][1]])
plt.plot(fMeshDict['globCoord'][fMeshDict['bNodes'][1]:fMeshDict['bNodes'][2]], fMeshDict['vx'][fMeshDict['bNodes'][1]:fMeshDict['bNodes'][2]])
plt.show()
plt.plot(fMeshDict['globCoord'][:fMeshDict['bNodes'][1]], fMeshDict['vy'][:fMeshDict['bNodes'][1]])
plt.plot(fMeshDict['globCoord'][fMeshDict['bNodes'][1]:fMeshDict['bNodes'][2]], fMeshDict['vy'][fMeshDict['bNodes'][1]:fMeshDict['bNodes'][2]])
plt.show()
plt.plot(fMeshDict['globCoord'][:fMeshDict['bNodes'][1]], fMeshDict['vz'][:fMeshDict['bNodes'][1]])
plt.plot(fMeshDict['globCoord'][fMeshDict['bNodes'][1]:fMeshDict['bNodes'][2]], fMeshDict['vz'][fMeshDict['bNodes'][1]:fMeshDict['bNodes'][2]])
plt.show()
quit()
self.nElmUpperPrev = len(fMeshDict['locCoord'][:fMeshDict['bNodes'][1]])
self.vSolver.SetVelocities(0, 0, fMeshDict['vx'][:fMeshDict['bNodes'][1]], fMeshDict['vy'][:fMeshDict['bNodes'][1]], fMeshDict['vz'][:fMeshDict['bNodes'][1]])
self.vSolver.SetVelocities(0, 1, fMeshDict['vx'][fMeshDict['bNodes'][1]:fMeshDict['bNodes'][2]], fMeshDict['vy'][fMeshDict['bNodes'][1]:fMeshDict['bNodes'][2]], fMeshDict['vz'][fMeshDict['bNodes'][1]:fMeshDict['bNodes'][2]])
......
vii/pyVII/vCoupler2.py
+ 1
− 1
View file @ 62568607
......@@ -30,7 +30,7 @@ class Coupler:
self.iSolverAPI = _iSolverAPI
self.vSolver = vSolver
self.maxCouplIter = 5
self.maxCouplIter = 15
self.couplTol = 1e-4
self.tms=fwk.Timers()
......
vii/pyVII/vInterpolator.py
+ 26
− 22
View file @ 62568607
......@@ -75,6 +75,7 @@ class Interpolator:
self.xCp = [[np.zeros((len(self.Xp[iReg]), 1)) for iReg in range(2)] for _ in range(len(self.cfg['Sections']))]
self.CpIter = [[np.zeros((len(self.Xp[iReg]), 1)) for iReg in range(2)] for _ in range(len(self.cfg['Sections']))]
def GetInvBC(self, vSolver, couplIter):
self.tms['InvToVisc'].start()
......@@ -111,10 +112,11 @@ class Interpolator:
for iReg in range(2):
# Y Z inverted
x_pos = self.Xp[iReg][:,(abs(self.Yp[iReg][0,:]-y)).argmin()]
y_pos = self.Zp[iReg][:,(abs(self.Yp[iReg][0,:]-y)).argmin()]
z_pos = self.Yp[iReg][:,(abs(self.Yp[iReg][0,:]-y)).argmin()]
goodRows = self.sortedRows[iReg][:,(abs(self.Yp[iReg][0,:]-y)).argmin()]
jcol = (abs(self.Yp[iReg][0,:]-y)).argmin()
x_pos = self.Xp[iReg][:,jcol]
y_pos = self.Yp[iReg][:,jcol]
z_pos = self.Zp[iReg][:,jcol]
goodRows = self.sortedRows[iReg][:,jcol]
ux = np.zeros(0)
uy = np.zeros(0)
......@@ -129,27 +131,25 @@ class Interpolator:
row = np.where(self.vdata[iReg][:,3]==goodRows[i])[0]
xxx = np.append(xxx, self.vdata[iReg][row,0])
ux = np.append(ux, Ux[iReg][row])
uy = np.append(uy, Uz[iReg][row])
uz = np.append(uz, Uy[iReg][row])
uy = np.append(uy, Uy[iReg][row])
uz = np.append(uz, Uz[iReg][row])
me = np.append(me, Me[iReg][row])
rhoe = np.append(rhoe, Rhoe[iReg][row])
cp = np.append(cp, Cp[iReg][row])
self.xCp[iSec][iReg] = np.column_stack((self.xCp[iSec][iReg], abs(x_pos)))
self.CpIter[iSec][iReg] = np.column_stack((self.CpIter[iSec][iReg], cp))
if couplIter == 0:
if iReg == 0:
self.stgPt[iSec] = np.argmin(ux**2+uy**2)
self.stgPt[iSec] = np.argmin(ux**2+uz**2)
xUp, xLw = self.__Remesh(abs(x_pos), self.stgPt[iSec])
yUp, yLw = self.__Remesh(y_pos, self.stgPt[iSec])
zUp, zLw = self.__Remesh(z_pos, self.stgPt[iSec])
vSolver.SetGlobMesh(iSec, 0, xUp, yUp, zUp)
vSolver.SetGlobMesh(iSec, 1, xLw, yLw, zLw)
vSolver.SetGlobMesh(iSec, 0, xUp, zUp, yUp)
vSolver.SetGlobMesh(iSec, 1, xLw, zLw, yLw)
self.DeltaStarPrev[iSec][0] = np.zeros(len(xUp))
self.DeltaStarPrev[iSec][1] = np.zeros(len(xLw))
......@@ -157,7 +157,7 @@ class Interpolator:
self.xxPrev[iSec][0] = np.zeros(len(xUp))
self.xxPrev[iSec][1] = np.zeros(len(xLw))
else:
vSolver.SetGlobMesh(iSec, 2, x_pos, y_pos, z_pos)
vSolver.SetGlobMesh(iSec, 2, x_pos, z_pos, y_pos)
self.DeltaStarPrev[iSec][2] = np.zeros(len(x_pos))
self.xxPrev[iSec][2] = np.zeros(len(x_pos))
......@@ -172,8 +172,8 @@ class Interpolator:
UyUp, UyLw = self.__Remesh(uy, self.stgPt[iSec])
UzUp, UzLw = self.__Remesh(uz, self.stgPt[iSec])
vSolver.SetVelocities(iSec, 0, UxUp, UyUp, UzUp)
vSolver.SetVelocities(iSec, 1, UxLw, UyLw, UzLw)
vSolver.SetVelocities(iSec, 0, UxUp, UzUp, UyUp)
vSolver.SetVelocities(iSec, 1, UxLw, UzLw, UyLw)
MeUp, MeLw = self.__Remesh(me, self.stgPt[iSec])
vSolver.SetMe(iSec, 0, MeUp)
......@@ -184,7 +184,7 @@ class Interpolator:
vSolver.SetRhoe(iSec, 1, RhoLw)
else:
vSolver.SetVelocities(iSec, 2, ux, uy, uz)
vSolver.SetVelocities(iSec, 2, ux, uz, uy)
vSolver.SetMe(iSec, 2, me)
vSolver.SetRhoe(iSec, 2, rhoe)
......@@ -239,10 +239,14 @@ class Interpolator:
blwAll[:,indexOutside[0]] = dummy
for iInter, inter in enumerate(indexInside):
dblw = blw[iInter+1][iReg] - blw[iInter][iReg]
dy = self.Cg_Sec[iInter+1][iReg][0,1] - self.Cg_Sec[iInter][iReg][0,1]
for j in inter:
dblw = blw[iInter+1][iReg] - blw[iInter][iReg]
dy = self.Cg_Sec[iInter+1][iReg][0,2] - self.Cg_Sec[iInter][iReg][0,2]
blwAll[:, j] = dblw/dy*(self.Yc[iReg][0, j] - self.Cg_Sec[iInter][iReg][0,2]) + blw[iInter][iReg]
if dy == 0:
print('dy=0 encountered. mem1 = ',self.Cg_Sec[iInter+1][iReg][0,1], 'mem2 = ',self.Cg_Sec[iInter][iReg][0,1])
print(iInter)
quit()
blwAll[:, j] = dblw/dy*(self.Yc[iReg][0, j] - self.Cg_Sec[iInter][iReg][0,1]) + blw[iInter][iReg]
dummy = np.zeros((len(self.Xc[iReg][:,0]), len(indexOutside[1])))
for iCol in range(len(dummy[0,:])):
......@@ -253,7 +257,7 @@ class Interpolator:
"""fig = plt.figure()
ax = fig.add_subplot(projection='3d')
for i in range(len(blw)):
ax.scatter(self.Cg_Sec[i][iReg][:,0], self.Cg_Sec[i][iReg][:,2], blw[i][iReg], s=0.7, color='red')
ax.scatter(self.Cg_Sec[i][iReg][:,0], self.Cg_Sec[i][iReg][:,1], blw[i][iReg], s=0.7, color='red')
for i in range(len(blwAll[0,:])):
ax.scatter(self.Xc[iReg][:,i], self.Yc[iReg][:,i], blwAll[:,i], s=0.2, color='blue')
ax.set_xlabel('x')
......@@ -277,9 +281,9 @@ class Interpolator:
"""fig = plt.figure()
ax = fig.add_subplot(projection='3d')
for i in range(len(blw)):
ax.scatter(self.Cg_Sec[i][iReg][:,0], self.Cg_Sec[i][iReg][:,2], blw[i][iReg], s=0.7, color='red')
ax.scatter(self.Cg_Sec[i][iReg][:,0], self.Cg_Sec[i][iReg][:,1], blw[i][iReg], s=0.7, color='red')
ax.scatter(self.icg[iReg][:, 0], self.icg[iReg][:, 1], mappedBlw, s=0.2, color='blue')
ax.set_zlim([-0.01, 0.015])
#ax.set_zlim([-0.01, 0.015])
ax.set_xlabel('x')
ax.set_ylabel('y')
ax.set_zlabel('z')
......@@ -288,7 +292,7 @@ class Interpolator:
self.tms['Interpolation'].stop()
if iReg==0:
if iReg==0 or iReg == 1:
for iElm in range(len(self.idata[iReg][:,0])):
self.blw[iReg].setU(iElm, mappedBlw[iElm])
......
vii/src/wBLRegion.cpp
+ 1
− 5
View file @ 62568607
......@@ -70,13 +70,9 @@ void BLRegion::SetStagBC(double Re)
double dv0 = (invBnd->GetVt(1) - invBnd->GetVt(0)) / (mesh->GetLoc(1) - mesh->GetLoc(0));
if (dv0 > 0)
{
U[0] = sqrt(0.075 / (Re * dv0)); // Theta
} // End if
else
{
U[0] = 1e-6; // Theta
} // end else
U[1] = 2.23; // H
U[2] = 0.; // N
U[3] = invBnd->GetVt(0); // ue
......@@ -88,7 +84,7 @@ void BLRegion::SetStagBC(double Re)
{
Ret[0] = 1.;
U[3] = Ret[0] / (U[0] * Re);
} // End if
}
//std::cout << "name " << name << " dv0 " << dv0 << " U3 " << U[3] << " U[0] " << U[0] << " \n" << std::scientific;
......
vii/src/wViscSolver.cpp
+ 37
− 15
View file @ 62568607
......@@ -34,8 +34,6 @@ ViscSolver::ViscSolver(double _Re, double _Minf, double _CFL0, unsigned int nSec
Re = _Re;
CFL0 = _CFL0;
Cdt = 0.; Cdf = 0.; Cdp = 0.;
verbose = _verbose;
/* Initialzing boundary layer */
......@@ -53,6 +51,11 @@ ViscSolver::ViscSolver(double _Re, double _Minf, double _CFL0, unsigned int nSec
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.;
/* Pre/post processing flags */
solverExitCode = 0;
......@@ -193,9 +196,7 @@ int ViscSolver::Run(unsigned int const couplIter)
else // 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> UpperCond(Sections[iSec][iRegion]->GetnVar(), 0.);
......@@ -210,18 +211,16 @@ int ViscSolver::Run(unsigned int const couplIter)
if (verbose>2){std::cout << "Wake BC imposed." <<std::endl;}
}
/* Loop over points */
// Loop over points
for (size_t iPoint = 1; iPoint < Sections[iSec][iRegion]->mesh->GetnMarkers(); ++iPoint)
{
/* Initialize solution at point */
// Initialize solution at point
tSolver->InitialCondition(iPoint, Sections[iSec][iRegion]);
// Solve equations
pointExitCode = tSolver->Integration(iPoint, Sections[iSec][iRegion]);
/* Unsucessfull convergence output */
// Unsucessfull convergence output
if (pointExitCode != 0)
{
if (verbose>1){std::cout << "Point " << iPoint << ": Convergence terminated with exitcode " << pointExitCode << std::endl;}
......@@ -266,7 +265,10 @@ int ViscSolver::Run(unsigned int const couplIter)
if (verbose>0){std::cout << "\n";}
} // End for iSections
ComputeDrag(Sections[0]);
for (size_t i=0; i<Sections.size(); ++i)
{
ComputeSectionalDrag(Sections[i], i);
}
if (solverExitCode != 0 && verbose>1){std::cout << "some points did not converge. ";}
if (verbose>0){std::cout << "\n" << std::endl;}
......@@ -384,7 +386,7 @@ void ViscSolver::AverageTransition(size_t iPoint, BLRegion *bl, int forced)
* Cdf = integral(Cf dx)_upper + integral(Cf dx)_lower.
* Cdp = Cdtot - Cdf.
*/
void ViscSolver::ComputeDrag(std::vector<BLRegion *> const &bl)
void ViscSolver::ComputeSectionalDrag(std::vector<BLRegion *> const &bl, int i)
{
unsigned int nVar = bl[0]->GetnVar();
......@@ -416,15 +418,35 @@ void ViscSolver::ComputeDrag(std::vector<BLRegion *> const &bl)
Cdf_lower += (frictioncoeff_lower[k] + frictioncoeff_lower[k - 1]) * (bl[1]->mesh->Getx(k) - bl[1]->mesh->Getx(k - 1)) / 2;
}
Cdf = Cdf_upper + Cdf_lower; // No friction in the wake
Cdf_sec[i] = Cdf_upper + Cdf_lower; // No friction in the wake
/* Compute total drag coefficient. */
Cdt = 2 * bl[2]->U[lastWkPt + 0] * pow(bl[2]->U[lastWkPt + 3], (bl[2]->U[lastWkPt + 1] + 5) / 2);
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. */
Cdp = Cdt - Cdf;
Cdp_sec[i] = Cdt_sec[i] - Cdf_sec[i];
}
void ViscSolver::ComputeTotalDrag()
{
if (Sections.size() == 1)
{
Cdt = Cdt_sec[0];
Cdf = Cdf_sec[0];
Cdp = Cdf_sec[0];
}
else
{
for(size_t k=1; k<Cdt_sec.size(); ++k)
{
Cdt += (Sections[i]->)
}
}
}
void ViscSolver::SetGlobMesh(size_t iSec, size_t iRegion, std::vector<double> _x, std::vector<double> _y, std::vector<double> _z)
......
vii/src/wViscSolver.h
+ 5
− 1
View file @ 62568607
......@@ -12,6 +12,10 @@ public:
double Cdt,
Cdf,
Cdp;
std::vector<double> Cdt_sec,
Cdf_sec,
Cdp_sec;
std::vector<std::vector<BLRegion *>> Sections;
......@@ -58,7 +62,7 @@ public:
private:
void CheckWaves(size_t iPoint, BLRegion *bl);
void AverageTransition(size_t iPoint, BLRegion *bl, int forced);
void ComputeDrag(std::vector<BLRegion *> const &bl);
void ComputeSectionalDrag(std::vector<BLRegion *> const &bl, int i);
void ComputeBlwVel();
void PrintBanner() const;
};
......
vii/tests/bli3.py
+ 3
− 3
View file @ 62568607
......@@ -68,8 +68,8 @@ def main():
print(ccolors.ANSI_BLUE + 'PyMeshing...' + ccolors.ANSI_RESET)
tms['msh'].start()
pars = {'spn' : spn, 'lgt' : lgt, 'wdt' : wdt, 'hgt' : hgt, 'msLe' : nms, 'msTe' : nms, 'msF' : fms}
imsh = viscU.mesh('/home/paul/lab/dartflo/dart/models/n0012_3.geo', pars)
vmsh = viscU.mesh('/home/paul/lab/dartflo/vii/models/n0012_3_visc.geo', pars)
imsh = viscU.mesh('/Users/paulzer/lab/dartflo/dart/models/n0012_3.geo', pars)
vmsh = viscU.mesh('/Users/paulzer/lab/dartflo/vii/models/n0012_3_visc.geo', pars)
msh, gmshWriter = floD.mesh(dim, 'models/n0012_3.geo', pars, ['field', 'wing', 'symmetry', 'downstream'], wktp = 'wakeTip')
tms['msh'].stop()
......@@ -80,7 +80,7 @@ def main():
tms['pre'].stop()
# solve problem
config = {'nDim': dim, 'Sections':[0.026, 0.103, 0.179, 0.256, 0.333, 0.41, 0.436, 0.513, 0.59, 0.667, 0.744, 0.821, 0.897, 0.949], 'span':spn, 'rbftype': 'linear', 'smoothing': 1e-8, 'degree': 0, 'neighbors': 5}
config = {'nDim': dim, 'Sections':[0.026, 0.103, 0.179, 0.256, 0.333, 0.41, 0.436, 0.513, 0.59, 0.667, 0.744, 0.821, 0.897, 0.949], 'span':spn, 'rbftype': 'linear', 'smoothing': 0.2, 'degree': 0, 'neighbors': 20}
iSolverAPI = vInterpol.Interpolator(floD.newton(pbl), blw, imsh, vmsh, config)
vSolver = viscU.initBL(Re, M_inf, CFL0, len(config['Sections']), 2)
#iSolverAPI = viscAPI.dartAPI(floD.newton(pbl), bnd, wk, nSections, vInterp)
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment