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

Iterations 3d

parent ebe605ce
No related branches found
No related tags found
No related merge requests found
Pipeline #7339 failed
Hide whitespace changes
Inline Side-by-side
vii/pyVII/vCoupler2.py
+ 1
− 1
View file @ a366558a
......@@ -30,7 +30,7 @@ class Coupler:
self.iSolverAPI = _iSolverAPI
self.vSolver = vSolver
self.maxCouplIter = 150
self.maxCouplIter = 5
self.couplTol = 1e-4
self.tms=fwk.Timers()
......
vii/pyVII/vInterpolator.py
+ 114
− 44
View file @ a366558a
......@@ -28,11 +28,10 @@ class Interpolator:
self.idata, self.icg = self.GetMesh(_imsh)
self.vdata, self.vcg = self.GetMesh(_vmsh)
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)
fig = plt.figure()
"""fig = plt.figure()
ax = fig.add_subplot(projection='3d')
ax.scatter(self.vdata[0][:,0], self.vdata[0][:,1], self.vdata[0][:,2], s=0.3)
ax.scatter(self.vcg[0][:,0], self.vcg[0][:,1], self.vcg[0][:,2], s=0.3)
......@@ -56,7 +55,7 @@ class Interpolator:
ax.set_zlabel('z')
#plt.gca().invert_yaxis()
plt.axis('off')
plt.show()
plt.show()"""
self.iSolver = _dartAPI
self.blw = _blw
......@@ -72,8 +71,10 @@ class Interpolator:
#self.uePrev = [[np.zeros(len(self.viscStruct[iSec][iReg][:,0])) for iReg in range(3)] for iSec in range(len(self.viscStruct))]
self.DeltaStarPrev = [[np.zeros(0) for iReg in range(3)] for iSec in range(len(self.cfg['Sections']))]
self.Cg_Sec = [[np.zeros((len(self.Xc[iReg])-1, 3)) for iReg in range(2)] for _ in range(len(self.cfg['Sections']))]
self.Cg_Sec = [[np.zeros((len(self.Xc[iReg]), 3)) for iReg in range(2)] for _ in range(len(self.cfg['Sections']))]
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()
......@@ -86,7 +87,7 @@ class Interpolator:
U = [np.zeros((len(self.idata[i][:,0]), 3)) for i in range(len(self.idata))]
M = [np.zeros(len(self.idata[i][:,0])) for i in range(len(self.idata))]
Rho = [np.zeros(len(self.idata[i][:,0])) for i in range(len(self.idata))]
CpII = [np.zeros(len(self.idata[i][:,0])) for i in range(len(self.idata))]
# Extract variables.
for iReg in range(len(self.idata)):
......@@ -95,14 +96,16 @@ 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]
CpII[iReg][iNode] = self.iSolver.Cp[row]
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] if nD==3 else 0], U[1][:,0], self.vdata[1][:,[0,2] if nD==3 else 0])]
Uy = [self.__RbfToSections(self.idata[0][:,:nD], U[0][:,1], self.vdata[0][:,:nD]), self.__RbfToSections(self.idata[1][:, [0,2] if nD==3 else 0], U[1][:,1], self.vdata[1][:,[0,2] if nD==3 else 0])]
Uz = [self.__RbfToSections(self.idata[0][:,:nD], U[0][:,2], self.vdata[0][:,:nD]), self.__RbfToSections(self.idata[1][:, [0,2] if nD==3 else 0], U[1][:,2], self.vdata[1][:,[0,2] if nD==3 else 0])]
Me = [self.__RbfToSections(self.idata[0][:,:nD], M[0], self.vdata[0][:,:nD]), self.__RbfToSections(self.idata[1][:, [0,2] if nD==3 else 0], M[1], self.vdata[1][:,[0,2] if nD==3 else 0])]
Rhoe = [self.__RbfToSections(self.idata[0][:,:nD], Rho[0], self.vdata[0][:,:nD]), self.__RbfToSections(self.idata[1][:, [0,2] if nD==3 else 0], Rho[1], self.vdata[1][:,[0,2] if nD==3 else 0])]
Cp = [self.__RbfToSections(self.idata[0][:,:nD], CpII[0], self.vdata[0][:,:nD]), self.__RbfToSections(self.idata[1][:, [0,2] if nD==3 else 0], CpII[1], self.vdata[1][:,[0,2] if nD==3 else 0])]
for iSec, y in enumerate(self.cfg['Sections']):
for iReg in range(2):
......@@ -118,6 +121,7 @@ class Interpolator:
uz = np.zeros(0)
me = np.zeros(0)
rhoe = np.zeros(0)
cp = np.zeros(0)
xxx = np.zeros(0)
for i in range(len(x_pos)):
......@@ -129,7 +133,13 @@ class Interpolator:
uz = np.append(uz, Uy[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:
......@@ -157,10 +167,6 @@ class Interpolator:
self.Cg_Sec[iSec][iReg][i,2] = z_pos[i+1] + z_pos[i]
self.Cg_Sec[iSec][iReg]/=2
plt.plot(x_pos, y_pos, 'x-')
plt.plot(self.Cg_Sec[iSec][iReg][:,0], self.Cg_Sec[iSec][iReg][:,1], 'o')
plt.show()
if iReg == 0:
UxUp, UxLw = self.__Remesh(ux, self.stgPt[iSec])
UyUp, UyLw = self.__Remesh(uy, self.stgPt[iSec])
......@@ -188,47 +194,103 @@ class Interpolator:
def SetBlowingVel(self, vSolver, couplIter):
for iSec in range(len(self.viscStruct)):
for iSec in range(len(self.cfg['Sections'])):
for iReg in range(3):
self.DeltaStarPrev[iSec][iReg] = vSolver.ExtractDeltaStar(iSec, iReg)
self.xxPrev[iSec][iReg] = vSolver.Extractxx(iSec, iReg)
blw = [np.zeros(0) for _ in range(2)]
blw = [[np.zeros(0) for _ in range(2)] for _ in range(len(self.cfg['Sections']))]
for iSec in range(len(self.viscStruct)):
for iSec in range(len(self.cfg['Sections'])):
for iReg in range(2):
if iReg == 0:
blwUp = vSolver.ExtractBlowingVel(iSec, 0)
blwLw = vSolver.ExtractBlowingVel(iSec, 1)
blw[iReg] = np.concatenate((blw[iReg], blwUp[::-1] , blwLw))
blw[iSec][iReg] = np.concatenate((blw[iSec][iReg], blwUp[::-1] , blwLw))
elif iReg == 1:
blw[iReg] = np.concatenate((blw[iReg], vSolver.ExtractBlowingVel(iSec, 2)))
blw[iSec][iReg] = np.concatenate((blw[iSec][iReg], vSolver.ExtractBlowingVel(iSec, 2)))
blwAll = interp2d(self.Cg_Sec[0])
fig = plt.figure()
ax = fig.add_subplot(projection='3d')
ax.scatter(x[0][:,0], x[0][:,2], blw[0], color = 'red')
"""for sec in viscStruct_tr:
ax.scatter(sec[0][:,0], sec[0][:,2], sec[0][:,1], color = 'red')
#ax.scatter(sec[1][:,0], sec[1][:,2], sec[1][:,1], color = 'red')"""
#ax.set_zlim([-0.5, 0.5])
ax.set_xlabel('x')
ax.set_ylabel('y')
ax.set_zlabel('Blowing velocity')
#plt.axis('off')
plt.show()
self.tms['Interpolation'].start()
nD = self.cfg['nDim']
#wignggg = self.__RbfToSections(x[0], blw[0], self.invCg_tr[0][:,:3])
wakke = self.__RbfToSections(x[1][:, [0,2] if nD==3 else 0], blw[1], self.invCg_tr[1][:, [0,2] if nD==3 else 0])
mappedBlw = [self.__RbfToSections(x[0], blw[0], self.invCg_tr[0][:,:3]), self.__RbfToSections(x[1][:, [0,2] if nD==3 else 0], blw[1], self.invCg_tr[1][:, [0,2] if nD==3 else 0])]
self.tms['Interpolation'].stop()
for iElm in range(len(self.invCg_tr[0][:,0])):
self.blw[0].setU(iElm, mappedBlw[0][iElm])
for iElm in range(len(self.invCg_tr[1][:,0])):
self.blw[1].setU(iElm, mappedBlw[1][iElm])
blwAll = [np.zeros(np.shape(self.Xc[iReg])) for iReg in range(2)]
for iReg in range(2):
# Find index of all Sections in Y
jBounds = np.zeros(len(self.cfg['Sections']), dtype=int)
for i, y in enumerate(self.cfg['Sections']):
jBounds[i] = (abs(self.Yp[iReg][0,:]-y)).argmin()
# Find all index of column of cg which lies between sections
# len of indexInside = nb of intervals between sections (nSections-1)
indexInside = [[] for _ in range(len(self.cfg['Sections'])-1)]
for i in range(len(indexInside)):
indexInside[i] = np.where((self.Yp[iReg][0,jBounds[i]] <= self.Yc[iReg][0,:]) & (self.Yc[iReg][0,:] <= self.Yp[iReg][0,jBounds[i+1]]))[0]
indexOutside = [np.where((self.Yc[iReg][0,:] <= self.Yp[iReg][0,jBounds[0]]))[0], np.where((self.Yc[iReg][0,:] >= self.Yp[iReg][0,jBounds[-1]]))[0]]
blwAll = np.zeros((len(self.Xc[iReg][:,0]), len(self.Xc[iReg][0,:])))
dummy = np.zeros((len(self.Xc[iReg][:,0]), len(indexOutside[0])))
for iCol in range(len(dummy[0,:])):
dummy[:,iCol] = blw[0][iReg]
blwAll[:,indexOutside[0]] = dummy
for iInter, inter in enumerate(indexInside):
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]
dummy = np.zeros((len(self.Xc[iReg][:,0]), len(indexOutside[1])))
for iCol in range(len(dummy[0,:])):
dummy[:,iCol] = blw[-1][iReg]
blwAll[:,indexOutside[1]] = dummy
"""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')
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')
ax.set_ylabel('y')
ax.set_zlabel('z')
#plt.gca().invert_yaxis()
plt.show()"""
nD = self.cfg['nDim']
xyz = np.zeros((np.size(self.Xc[iReg]), 3))
blwForMap = np.zeros(np.size(self.Xc[iReg]))
ind=0
for i in range(len(self.Xc[iReg][:,0])):
for j in range(len(self.Xc[iReg][0,:])):
xyz[ind,:] = [self.Xc[iReg][i,j], self.Yc[iReg][i,j], self.Zc[iReg][i,j]]
blwForMap[ind] = blwAll[i,j]
ind+=1
mappedBlw = self.__RbfToSections(xyz[:, [0,2] if nD==3 else 0], blwForMap, self.icg[iReg][:, [0,2] if nD==3 else 0])
"""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.icg[iReg][:, 0], self.icg[iReg][:, 1], mappedBlw, s=0.2, color='blue')
ax.set_zlim([-0.01, 0.015])
ax.set_xlabel('x')
ax.set_ylabel('y')
ax.set_zlabel('z')
#plt.gca().invert_yaxis()
plt.show()"""
self.tms['Interpolation'].stop()
if iReg==0:
for iElm in range(len(self.idata[iReg][:,0])):
self.blw[iReg].setU(iElm, mappedBlw[iElm])
def RunSolver(self):
self.iSolver.run()
......@@ -647,7 +709,7 @@ class Interpolator:
data[0] = data[0][data[0][:,1]<=1, :]
cg = data.copy()
cg = cg.copy()
cg[0] = cg[0][cg[0][:,1]<=1, :]
......@@ -655,6 +717,7 @@ class Interpolator:
for i in range(len(sections)):
sections[i] = round(sections[i], 3)
sections = np.unique(sections)
sections = np.sort(sections)
nChord = [len(data[i][abs(data[i][:,1]-sections[0])<1e-3, 0]) for i in range(len(data))]
......@@ -683,10 +746,12 @@ class Interpolator:
sections_cg = cg[0][:,1]
for i in range(len(sections_cg)):
sections_cg[i] = round(sections_cg[i], 3)
sections_cg = np.unique(sections)
sections_cg = np.unique(sections_cg)
sections_cg = np.sort(sections_cg)
nChord_cg = [len(cg[i][abs(cg[i][:,1]-sections[0])<1e-3, 0]) for i in range(len(cg))]
nChord_cg = [len(cg[i][abs(cg[i][:,1]-cg[i][0,1])<1e-3, 0]) for i in range(len(cg))]
print(nChord_cg)
Xc = [np.zeros((nChord_cg[i], 0)) for i in range(len(data))]
Yc = [np.zeros((nChord_cg[i], 0)) for i in range(len(data))]
Zc = [np.zeros((nChord_cg[i], 0)) for i in range(len(data))]
......@@ -710,6 +775,11 @@ class Interpolator:
Yc[iReg] = np.column_stack((Yc[iReg], sort[:,1]))
Zc[iReg] = np.column_stack((Zc[iReg], sort[:,2]))
Xc[0] = np.delete(Xc[0], 0, 0)
Yc[0] = np.delete(Yc[0], 0, 0)
Zc[0] = np.delete(Zc[0], 0, 0)
sortedElems[0] = np.delete(sortedElems[0], 0, 0)
return Xp, Yp, Zp, sortedRows, Xc, Yc, Zc, sortedElems
......
vii/tests/bli3.py
+ 7
− 1
View file @ a366558a
......@@ -33,6 +33,8 @@ 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
......@@ -78,7 +80,7 @@ def main():
tms['pre'].stop()
# solve problem
config = {'nDim': dim, 'Sections':[0.02, 0.041], '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': 1e-8, 'degree': 0, 'neighbors': 5}
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)
......@@ -86,6 +88,10 @@ def main():
coupler.Run()
plt.plot(iSolverAPI.xCp[0][0][:,1], iSolverAPI.CpIter[0][0][:,1])
plt.plot(iSolverAPI.xCp[0][0][:,-1], iSolverAPI.CpIter[0][0][:,-1])
plt.show()
# display timers
tms['total'].stop()
print(ccolors.ANSI_BLUE + 'PyTiming...' + ccolors.ANSI_RESET)
......
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