From a366558ab820c66000b649accf702cae83847fef Mon Sep 17 00:00:00 2001
From: Paul <paul.dechamps@student.uliege.be>
Date: Mon, 11 Jul 2022 19:02:32 +0200
Subject: [PATCH] Iterations 3d
---
vii/pyVII/vCoupler2.py | 2 +-
vii/pyVII/vInterpolator.py | 158 ++++++++++++++++++++++++++-----------
vii/tests/bli3.py | 8 +-
3 files changed, 122 insertions(+), 46 deletions(-)
diff --git a/vii/pyVII/vCoupler2.py b/vii/pyVII/vCoupler2.py
index 1699635..048ae7f 100644
--- a/vii/pyVII/vCoupler2.py
+++ b/vii/pyVII/vCoupler2.py
@@ -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()
diff --git a/vii/pyVII/vInterpolator.py b/vii/pyVII/vInterpolator.py
index 8f7210e..d543849 100644
--- a/vii/pyVII/vInterpolator.py
+++ b/vii/pyVII/vInterpolator.py
@@ -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
diff --git a/vii/tests/bli3.py b/vii/tests/bli3.py
index 139f34e..326446c 100644
--- a/vii/tests/bli3.py
+++ b/vii/tests/bli3.py
@@ -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)
--
GitLab