From fcd2945e14ec65e6a94a4d378a991b79d95a884d Mon Sep 17 00:00:00 2001
From: Paul Dechamps <paulzer@sentinelle.local>
Date: Thu, 12 May 2022 13:55:21 +0200
Subject: [PATCH] Separating upper lower sides with normals
---
dart/default.py | 1 +
dart/pyVII/vCoupler.py | 12 +-
dart/pyVII/vCoupler2.py | 3 +-
dart/pyVII/vInterpolator.py | 623 ++++++++++++++++--------------------
dart/pyVII/vStructures.py | 46 ---
dart/pyVII/vUtils.py | 54 +---
dart/src/wBody.cpp | 4 +
dart/src/wClosures.cpp | 35 +-
dart/src/wViscFluxes.cpp | 1 -
dart/src/wViscSolver.cpp | 2 +
dart/tests/bli2.py | 156 +++++++++
dart/tests/bli3.py | 20 +-
dart/tests/bli_lowLift.py | 4 +-
dart/tests/bli_lowLift2.py | 29 +-
dart/tests/rae_3.py | 27 +-
15 files changed, 510 insertions(+), 507 deletions(-)
delete mode 100644 dart/pyVII/vStructures.py
create mode 100644 dart/tests/bli2.py
diff --git a/dart/default.py b/dart/default.py
index caf91dc..878ed74 100644
--- a/dart/default.py
+++ b/dart/default.py
@@ -42,6 +42,7 @@ def mesh(dim, file, pars, bnd, wk = 'wake', wktp = None):
name of crack (wake) tip physical tag
"""
import tbox.gmsh as gmsh
+ import copy
# load the mesh
msh = gmsh.MeshLoader(file,__file__).execute(**pars)
# crack the mesh
diff --git a/dart/pyVII/vCoupler.py b/dart/pyVII/vCoupler.py
index cde32fe..f74cd4c 100644
--- a/dart/pyVII/vCoupler.py
+++ b/dart/pyVII/vCoupler.py
@@ -19,6 +19,9 @@
## VII Coupler
# Paul Dechamps
+import dart.utils as floU
+import tbox.utils as tboxU
+
from dart.viscous.Master.DAirfoil import Airfoil
from dart.viscous.Master.DWake import Wake
from dart.viscous.Drivers.DGroupConstructor import GroupConstructor
@@ -39,7 +42,7 @@ class Coupler:
self.vSolver.verbose = self.iSolver.verbose
#self.vSolver.verbose = 0
- self.maxCouplIter = 1
+ self.maxCouplIter = 150
self.couplTol = 1e-4
self.tms=fwk.Timers()
@@ -79,6 +82,11 @@ class Coupler:
self.iSolver.run()
self.tms['inviscid'].stop()
+ if couplIter == 0:
+ # Write inviscid Cp distribution
+ Cp = floU.extract(self.group[0].boundary.tag.elems, self.iSolver.Cp)
+ tboxU.write(Cp, 'Cp_Inviscid.dat', '%1.5e', ', ', 'x, y, z, Cp', '')
+
# Extract inviscid boundary
self.tms['processing'].start()
self.__ImposeInvBC(couplIter)
@@ -101,7 +109,7 @@ class Coupler:
dragIter = np.vstack((dragIter, [self.iSolver.Cl, self.vSolver.Cdt]))
- if error <= self.couplTol and exitCode==0:
+ if error <= self.couplTol:
print('{:>5s}| {:>10s} {:>10s} | {:>10s} {:>8s} {:>8s}'.format('It', 'Cl', 'Cd', 'xtrT', 'xtrB', 'Error'))
print(ccolors.ANSI_GREEN,'{:>4.0f}| {:>10.5f} {:>10.5f} | {:>10.4f} {:>8.4f} {:>8.2f}\n'.format(couplIter,
self.iSolver.Cl, self.vSolver.Cdt, xtr[0], xtr[1], np.log10(error)),ccolors.ANSI_RESET)
diff --git a/dart/pyVII/vCoupler2.py b/dart/pyVII/vCoupler2.py
index 1d86205..1699635 100644
--- a/dart/pyVII/vCoupler2.py
+++ b/dart/pyVII/vCoupler2.py
@@ -35,8 +35,6 @@ class Coupler:
self.tms=fwk.Timers()
self.resetInv = False
-
- self.iSolverAPI.iSolver.printOn=False
pass
def Run(self):
@@ -58,6 +56,7 @@ class Coupler:
# Run inviscid solver
self.tms['inviscid'].start()
if self.resetInv:
+ print('Inviscid solver reset')
self.iSolverAPI.ResetSolver()
self.iSolverAPI.RunSolver()
self.tms['inviscid'].stop()
diff --git a/dart/pyVII/vInterpolator.py b/dart/pyVII/vInterpolator.py
index 1304404..7323c2f 100644
--- a/dart/pyVII/vInterpolator.py
+++ b/dart/pyVII/vInterpolator.py
@@ -1,9 +1,11 @@
+from gc import isenabled
+from xmlrpc.client import INVALID_ENCODING_CHAR
import numpy as np
from matplotlib import pyplot as plt
from fwk.coloring import ccolors
import dart.pyVII.vUtils as blxU
-import dart.pyVII.vStructures as vStruct
+#import dart.pyVII.vStructures as vStruct
# 3D spline with Bézier curves
from scipy.interpolate import bisplrep
@@ -19,26 +21,25 @@ import fwk
class Interpolator:
- def __init__(self, _dartSolver, bnd_wing, bnd_wake, _cfg):
- self.iSolver = _dartSolver
-
- self.invStruct = vStruct.invStruct()
- self.viscStruct = vStruct.viscStruct()
-
- self.bndWing = bnd_wing
- self.bndWake = bnd_wake
- self.config = _cfg
- if self.config['nDim'] == 2:
- self.nodesList, self.transformedCoord, self.wakeCoord, self.viscStruct, self.viscStruct_tr = self.__GetAirfoil(bnd_wing, bnd_wake, _cfg)
- elif self.config['nDim'] == 3:
-
+ def __init__(self, _dartAPI, _blw, _cfg):
+
+ self.iSolver = _dartAPI
+ self.blw = _blw
- print('Sections created')
+ # Create data structures for interpolation
+ self.tms = fwk.Timers()
+ self.tms['GetWing'].start()
+ self.__GetWing(_blw, _cfg)
+ self.tms['GetWing'].stop()
+ print(self.tms)
+ self.cfg = _cfg
- self.xxPrev = [[None for iReg in range(3)] for iSec in range(1)]
+ self.stgPt=np.zeros(self.cfg['nSections'], dtype=int)
+
+ self.xxPrev = [[np.zeros(0) for iReg in range(3)] for iSec in range(self.cfg['nSections'])]
#self.uePrev = [[np.zeros(len(self.viscStruct[iSec][iReg][:,0])) for iReg in range(3)] for iSec in range(len(self.viscStruct))]
- self.DeltaStarPrev = [[None for iReg in range(3)] for iSec in range(1)]
- self.tms = fwk.Timers()
+ self.DeltaStarPrev = [[np.zeros(0) for iReg in range(3)] for iSec in range(self.cfg['nSections'])]
+
def GetInvBC(self, vSolver, couplIter):
self.tms['InvToVisc'].start()
@@ -49,50 +50,36 @@ class Interpolator:
# Wing
- U_wg = np.zeros((len(self.nodesList), 3))
- M_wg = np.zeros(len(self.nodesList))
- Rho_wg = np.zeros(len(self.nodesList))
- for iNode, n in enumerate(self.nodesList):
- U_wg[iNode,0] = self.iSolver.U[n.row][0]
- U_wg[iNode,1] = self.iSolver.U[n.row][1]
- U_wg[iNode,2] = self.iSolver.U[n.row][2]
- M_wg[iNode] = self.iSolver.M[n.row]
- Rho_wg[iNode] = self.iSolver.rho[n.row]
-
- # Wake
+ U = [np.zeros((len(self.invStruct[i][:,0]), 3)) for i in range(len(self.invStruct))]
+ M = [np.zeros(len(self.invStruct[i][:,0])) for i in range(len(self.invStruct))]
+ Rho = [np.zeros(len(self.invStruct[i][:,0])) for i in range(len(self.invStruct))]
- U_wk =np.zeros((len(self.bndWake.nodes), 3))
- M_wk = np.zeros(len(self.bndWake.nodes))
- Rho_wk = np.zeros(len(self.bndWake.nodes))
+ # Extract variables.
+
+ for iReg in range(len(self.invStruct)):
+ for iNode, row in enumerate(self.invStruct[iReg][:,3]):
+ row = int(row)
+ 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]
- 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]
- U_wk[iNode,2] = self.iSolver.U[n.row][2]
- M_wk[iNode] = self.iSolver.M[n.row]
- Rho_wk[iNode] = self.iSolver.rho[n.row]
-
- if self.config['nDim'] == 3: # x and y inverted in 3d
- U_wg[:, [1,2]] = U_wg[:,[2,1]]
- U_wk[:, [1,2]] = U_wk[:,[2,1]]
+ if self.cfg['nDim'] == 3:
+ U[iReg][:,[1,2]] = U[iReg][:,[2,1]]
- Ux = [self.__RbfToSections(self.transformedCoord[:,:2], U_wg[:,0], self.viscStruct_tr[0][0][:,:2]), self.__RbfToSections(self.wakeCoord[:, :1], U_wk[:,0], self.viscStruct[0][1][:,:1])]
- Uy = [self.__RbfToSections(self.transformedCoord[:,:2], U_wg[:,1], self.viscStruct_tr[0][0][:,:2]), self.__RbfToSections(self.wakeCoord[:, :1], U_wk[:,1], self.viscStruct[0][1][:,:1])]
- Uz = [self.__RbfToSections(self.transformedCoord[:,:2], U_wg[:,2], self.viscStruct_tr[0][0][:,:2]), self.__RbfToSections(self.wakeCoord[:, :1], U_wk[:,2], self.viscStruct[0][1][:,:1])]
-
- Me = [self.__RbfToSections(self.transformedCoord[:,:2], M_wg, self.viscStruct_tr[0][0][:,:2]), self.__RbfToSections(self.wakeCoord[:, :1], M_wk, self.viscStruct[0][1][:,:1])]
-
- Rho = [self.__RbfToSections(self.transformedCoord[:,:2], Rho_wg, self.viscStruct_tr[0][0][:,:2]), self.__RbfToSections(self.wakeCoord[:, :1], Rho_wk, self.viscStruct[0][1][:,:1])]
-
# Find stagnation point
- for iSec in range(len(self.viscStruct)):
-
+ nD = self.cfg['nDim']
+ for iSec in range(len(self.viscStruct_tr)):
+ Ux = [self.__RbfToSections(self.invStruct_tr[0][:,:nD], U[0][:,0], self.viscStruct_tr[iSec][0][:,:nD]), self.__RbfToSections(self.invStruct_tr[1][:, [0,2] if nD==3 else 0], U[1][:,0], self.viscStruct[iSec][1][:,[0,2] if nD==3 else 0])]
+ Uy = [self.__RbfToSections(self.invStruct_tr[0][:,:nD], U[0][:,1], self.viscStruct_tr[iSec][0][:,:nD]), self.__RbfToSections(self.invStruct_tr[1][:, [0,2] if nD==3 else 0], U[1][:,1], self.viscStruct[iSec][1][:,[0,2] if nD==3 else 0])]
+ Uz = [self.__RbfToSections(self.invStruct_tr[0][:,:nD], U[0][:,2], self.viscStruct_tr[iSec][0][:,:nD]), self.__RbfToSections(self.invStruct_tr[1][:, [0,2] if nD==3 else 0], U[1][:,2], self.viscStruct[iSec][1][:,[0,2] if nD==3 else 0])]
+ Me = [self.__RbfToSections(self.invStruct_tr[0][:,:nD], M[0], self.viscStruct_tr[iSec][0][:,:nD]), self.__RbfToSections(self.invStruct_tr[1][:, [0,2] if nD==3 else 0], M[1], self.viscStruct[iSec][1][:,[0,2] if nD==3 else 0])]
+ Rhoe = [self.__RbfToSections(self.invStruct_tr[0][:,:nD], Rho[0], self.viscStruct_tr[iSec][0][:,:nD]), self.__RbfToSections(self.invStruct_tr[1][:, [0,2] if nD==3 else 0], Rho[1], self.viscStruct[iSec][1][:,[0,2] if nD==3 else 0])]
if couplIter == 0:
- self.stgPt = np.argmin(Ux[0]**2+Uy[0]**2)
+ self.stgPt[iSec] = np.argmin(Ux[0]**2+Uy[0]**2)
- xUp, xLw = self.__Remesh(self.viscStruct[iSec][0][:,0], self.stgPt)
- yUp, yLw = self.__Remesh(self.viscStruct[iSec][0][:,1], self.stgPt)
- zUp, zLw = self.__Remesh(self.viscStruct[iSec][0][:,2], self.stgPt)
+ xUp, xLw = self.__Remesh(self.viscStruct[iSec][0][:,0], self.stgPt[iSec])
+ yUp, yLw = self.__Remesh(self.viscStruct[iSec][0][:,1], self.stgPt[iSec])
+ zUp, zLw = self.__Remesh(self.viscStruct[iSec][0][:,2], self.stgPt[iSec])
vSolver.SetGlobMesh(iSec, 0, xUp, yUp, zUp)
vSolver.SetGlobMesh(iSec, 1, xLw, yLw, zLw)
@@ -100,122 +87,87 @@ class Interpolator:
self.viscStruct[iSec][1][:,1],
self.viscStruct[iSec][1][:,2])
- self.DeltaStarPrev[0][0] = np.zeros(len(xUp))
- self.DeltaStarPrev[0][1] = np.zeros(len(xLw))
- self.DeltaStarPrev[0][2] = np.zeros(len(self.viscStruct[iSec][1][:,0]))
+ self.DeltaStarPrev[iSec][0] = np.zeros(len(xUp))
+ self.DeltaStarPrev[iSec][1] = np.zeros(len(xLw))
+ self.DeltaStarPrev[iSec][2] = np.zeros(len(self.viscStruct[iSec][1][:,0]))
- self.xxPrev[0][0] = np.zeros(len(xUp))
- self.xxPrev[0][1] = np.zeros(len(xLw))
- self.xxPrev[0][2] = np.zeros(len(self.viscStruct[iSec][1][:,0]))
-
- UxUp, UxLw = self.__Remesh(Ux[0], self.stgPt)
- UyUp, UyLw = self.__Remesh(Uy[0], self.stgPt)
- UzUp, UzLw = self.__Remesh(Uz[0], self.stgPt)
+ self.xxPrev[iSec][0] = np.zeros(len(xUp))
+ self.xxPrev[iSec][1] = np.zeros(len(xLw))
+ self.xxPrev[iSec][2] = np.zeros(len(self.viscStruct[iSec][1][:,0]))
- if couplIter == 0:
- plt.plot(self.__Remesh(self.viscStruct[iSec][0][:,0], self.stgPt)[0], UxUp)
- plt.plot(self.__Remesh(self.viscStruct[iSec][0][:,0], self.stgPt)[1], UxLw)
- plt.show()
+ UxUp, UxLw = self.__Remesh(Ux[0], self.stgPt[iSec])
+ UyUp, UyLw = self.__Remesh(Uy[0], self.stgPt[iSec])
+ UzUp, UzLw = self.__Remesh(Uz[0], self.stgPt[iSec])
vSolver.SetVelocities(iSec, 0, UxUp, UyUp, UzUp)
vSolver.SetVelocities(iSec, 1, UxLw, UyLw, UzLw)
vSolver.SetVelocities(iSec, 2, Ux[1], Uy[1], Uz[1])
- MeUp, MeLw = self.__Remesh(Me[0], self.stgPt)
+ MeUp, MeLw = self.__Remesh(Me[0], self.stgPt[iSec])
vSolver.SetMe(iSec, 0, MeUp)
vSolver.SetMe(iSec, 1, MeLw)
vSolver.SetMe(iSec, 2, Me[1])
- RhoUp, RhoLw = self.__Remesh(Rho[0], self.stgPt)
+ RhoUp, RhoLw = self.__Remesh(Rhoe[0], self.stgPt[iSec])
vSolver.SetRhoe(iSec, 0, RhoUp)
vSolver.SetRhoe(iSec, 1, RhoLw)
- vSolver.SetRhoe(iSec, 2, Rho[1])
-
- vSolver.SetDeltaStarExt(iSec, 0, self.DeltaStarPrev[0][0])
- vSolver.SetDeltaStarExt(iSec, 1, self.DeltaStarPrev[0][1])
- vSolver.SetDeltaStarExt(iSec, 2, self.DeltaStarPrev[0][2])
-
- vSolver.SetxxExt(iSec, 0, self.xxPrev[0][0])
- vSolver.SetxxExt(iSec, 1, self.xxPrev[0][1])
- vSolver.SetxxExt(iSec, 2, self.xxPrev[0][2])
+ vSolver.SetRhoe(iSec, 2, Rhoe[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])"""
+ for iReg in range(3):
+ vSolver.SetDeltaStarExt(iSec, iReg, self.DeltaStarPrev[iSec][iReg])
+ vSolver.SetxxExt(iSec, iReg, self.xxPrev[iSec][iReg])
def SetBlowingVel(self, vSolver, couplIter):
- BlowingVel = [[np.zeros(0) for iReg in range(3)] for iSec in range(len(self.viscStruct))]
+
for iSec in range(len(self.viscStruct)):
- for iRegion in range(3):
- BlowingVel[iSec][iRegion] = np.asarray(vSolver.ExtractBlowingVel(iSec, iRegion))
- self.DeltaStarPrev[iSec][iRegion] = np.asarray(vSolver.ExtractDeltaStar(iSec, iRegion))
- self.xxPrev[iSec][iRegion] = np.asarray(vSolver.Extractxx(iSec, iRegion))
-
- if self.config['nDim'] == 2:
- BlwWing = np.concatenate((BlowingVel[0][0][::-1], BlowingVel[0][1]))
- BlwWake = BlowingVel[0][2]
-
- meany = (self.viscStruct[0][0][np.argmin(self.viscStruct[0][0][:,0]), 1] - self.viscStruct[0][0][np.argmax(self.viscStruct[0][0][:,0]), 1]) / 2
-
- # Compute cell centers
- viscCgWing = np.zeros((len(self.viscStruct[0][0][:,0]) - 1, 3))
- for iElm in range(len(self.viscStruct[0][0][:,0]) - 1):
- viscCgWing[iElm, 0] = self.viscStruct[0][0][iElm, 0] + self.viscStruct[0][0][iElm + 1, 0]
- viscCgWing[iElm, 1] = self.viscStruct[0][0][iElm, 1] + self.viscStruct[0][0][iElm + 1, 1]
- viscCgWing[iElm, 2] = self.viscStruct[0][0][iElm, 2] + self.viscStruct[0][0][iElm + 1, 2]
- viscCgWing/=2
- viscCgWing[viscCgWing[:,1]<meany, 0] *= -1
-
- viscCgWake = np.zeros((len(self.viscStruct[0][1][:,0]) - 1, 3))
- for iElm in range(len(self.viscStruct[0][1][:,0]) - 1):
- viscCgWake[iElm, 0] = self.viscStruct[0][1][iElm, 0] + self.viscStruct[0][1][iElm + 1, 0]
- viscCgWake[iElm, 1] = self.viscStruct[0][1][iElm, 1] + self.viscStruct[0][1][iElm + 1, 1]
- viscCgWake[iElm, 2] = self.viscStruct[0][1][iElm, 2] + self.viscStruct[0][1][iElm + 1, 2]
- viscCgWake/=2
-
- invCgWing = np.zeros((len(self.bndWing.tag.elems), 3))
- for iElm, e in enumerate(self.bndWing.tag.elems):
- for n in e.nodes:
- invCgWing[iElm, 0] +=n.pos[0]
- invCgWing[iElm, 1] +=n.pos[1]
- invCgWing[iElm, 2] +=n.pos[2]
- invCgWing[iElm, :] /= e.nodes.size()
-
- invCgWing[invCgWing[:,1]<meany, 0] *= -1
+ for iReg in range(3):
+ self.DeltaStarPrev[iSec][iReg] = vSolver.ExtractDeltaStar(iSec, iReg)
+ self.xxPrev[iSec][iReg] = vSolver.Extractxx(iSec, iReg)
- # Wake
+ x = [np.zeros((0,3)) for _ in range(2)]
+ blw = [np.zeros(0) for _ in range(2)]
- invCgWake = np.zeros((len(self.bndWake.tag.elems), 3))
- for iElm, e in enumerate(self.bndWake.tag.elems):
- for n in e.nodes:
- invCgWake[iElm, 0] +=n.pos[0]
- invCgWake[iElm, 1] +=n.pos[1]
- invCgWake[iElm, 2] +=n.pos[2]
- invCgWake[iElm, :] /= e.nodes.size()
-
- if self.config['nDim'] == 3:
- invCgWing[:,[1,2]] = invCgWing[:,[2,1]]
- invCgWake[:,[1,2]] = invCgWake[:,[2,1]]
+ for iSec in range(len(self.viscStruct)):
+ for iReg in range(2):
+ if iReg == 0:
+ x[iReg] = np.vstack((x[iReg], self.viscCg_tr[iSec][iReg][:,:3]))
+ blwUp = vSolver.ExtractBlowingVel(iSec, 0)
+ blwLw = vSolver.ExtractBlowingVel(iSec, 1)
+ plt.plot(blwLw)
+ plt.show()
+ blw[iReg] = np.concatenate((blw[iReg], blwUp[::-1] , blwLw))
+ plt.plot(self.viscCg_tr[iSec][iReg][:,0],np.concatenate((blwUp[::-1] , blwLw)))
+ plt.title(iSec)
+ plt.show()
+ elif iReg == 1:
+ x[iReg] = np.vstack((x[iReg], self.viscCg_tr[iSec][iReg][:,:3]))
+ blw[iReg] = np.concatenate((blw[iReg], vSolver.ExtractBlowingVel(iSec, 2)))
self.tms['Interpolation'].start()
-
- mappedBlw = [self.__RbfToSections(viscCgWing[:, :2], BlwWing, invCgWing[:,:2]), self.__RbfToSections(viscCgWake[:, :1], BlwWake, invCgWake[:, :1])]
-
- sln = blxU.GetSolution(0, vSolver)
-
- self.tms['Interpolation'].stop()
- for iElm in range(len(self.bndWing.tag.elems)):
- self.bndWing.setU(iElm, mappedBlw[0][iElm])
- for iElm in range(len(self.bndWake.tag.elems)):
- self.bndWake.setU(iElm, mappedBlw[1][iElm])
-
+ nD = self.cfg['nDim']
+ mappedBlw = [self.__RbfToSections(x[0], blw[0], self.invCg_tr[0][:,:3], forceLin=1), self.__RbfToSections(x[1][:, [0,2] if nD==3 else 0], blw[1], self.invCg_tr[1][:, [0,2] if nD==3 else 0], forceLin=1)]
+ fig = plt.figure()
+ ax = fig.add_subplot(projection='3d')
+ ax.scatter(self.invCg_tr[0][:,0], self.invCg_tr[0][:,2], mappedBlw[0], s=0.3)
+ #ax.scatter(abs(invCg_tr[0][:,0]), invCg_tr[0][:,2], invCg_tr[0][:,1], s=0.3)
+ for sec in range(len(self.viscStruct_tr)):
+ ax.scatter(x[0][:,0], x[0][:,2], blw[0], color = 'red')
+ #ax.set_zlim([-0.5, 0.5])
+ ax.set_xlabel('x')
+ ax.set_ylabel('y')
+ ax.set_zlabel('z')
+ plt.show()
+ 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])
def RunSolver(self):
self.iSolver.run()
@@ -240,82 +192,116 @@ class Interpolator:
def GetCm(self):
return self.iSolver.Cm
- def __Remesh(self, vec, idx_stag):
- xUp = vec[:idx_stag+1]
+ def __Remesh(self, vec, stgPt):
+ xUp = vec[:stgPt+1]
xUp = xUp[::-1]
- xLw = vec[idx_stag:]
+ xLw = vec[stgPt:]
return xUp, xLw
- def __RbfToSections(self, x, var, xs):
+ def __RbfToSections(self, x, var, xs, forceLin=0):
if np.all(var == var[0]):
varOut = RBFInterpolator(x, var, neighbors=1, kernel='linear', smoothing=0.0, degree=0)(xs)
+ elif forceLin==1:
+ varOut = RBFInterpolator(x, var, neighbors=50, kernel='linear', smoothing=0.1, degree=0)(xs)
else:
- varOut = RBFInterpolator(x, var, neighbors=self.config['neighbors'], kernel=self.config['rbftype'], smoothing=self.config['smoothing'], degree=self.config['degree'])(xs)
+ varOut = RBFInterpolator(x, var, neighbors=self.cfg['neighbors'], kernel=self.cfg['rbftype'], smoothing=self.cfg['smoothing'], degree=self.cfg['degree'])(xs)
return varOut
-
- def __GetWing(self, bnd_wing, bnd_wake, config):
- if config['nDim'] == 2 and config['nSections'] != 1:
- raise RuntimeError('Cannot create more than 1 section on 2D geometry. Specified:', config['nSections'])
+
+ def __GetWing(self, blw, config):
viscStruct = [[np.zeros((config['nPoints'][1] if i==1 else 2*config['nPoints'][0]-1, 3)) for i in range(2)] for _ in range(config['nSections'])]
- viscStruct_tr = [[np.zeros((config['nPoints'][1] if i==1 else 2*config['nPoints'][0]-1, 3)) for i in range(1)] for _ in range(config['nSections'])]
+ viscStruct_tr = [[np.zeros((config['nPoints'][1] if i==1 else 2*config['nPoints'][0]-1, 3)) for i in range(2)] for _ in range(config['nSections'])]
+
+
+ invStruct = [np.zeros((len(blw[i].nodes), 4)) for i in range(len(blw))]
+ invStruct_tr = [np.zeros((len(blw[i].nodes), 4)) for i in range(len(blw))]
+
+ # Wing
+ upperNodes = np.zeros((0,4))
+ lowerNodes = np.zeros((0,4))
+
+ self.tms['ElmLoop'].start()
+ # Compute reference normal
+ nVec = [np.zeros(0) for _ in range(blw[1].tag.elems[0].nodes.size())]
+ for iNode, n in enumerate(blw[1].tag.elems[0].nodes):
+ nVec[iNode] = np.array([n.pos[0], n.pos[1], n.pos[2]])
+ normalRef = np.cross(nVec[1] - nVec[0], nVec[2] - nVec[0])
+ normalRef/=np.linalg.norm(normalRef)
+
+ for e in blw[0].tag.elems:
+ nVec = [np.zeros(0) for _ in range(e.nodes.size())]
+ for iNode, n in enumerate(e.nodes):
+ nVec[iNode] = np.array([n.pos[0], n.pos[1], n.pos[2]])
+ normal = np.cross(nVec[1] - nVec[0], nVec[2] - nVec[0])
+ normal/=np.linalg.norm(normalRef)
+ if np.dot(normal, normalRef) >= 0:
+ for n in e.nodes:
+ if n.row not in upperNodes[:,3] and n.row not in lowerNodes[:,3]:
+ upperNodes = np.vstack((upperNodes, [n.pos[0], n.pos[1], n.pos[2], n.row]))
+ if np.dot(normal, normalRef) < 0:
+ for n in e.nodes:
+ if n.row not in lowerNodes[:,3] and n.row not in upperNodes[:,3]:
+ lowerNodes = np.vstack((lowerNodes, [n.pos[0], n.pos[1], n.pos[2], n.row]))
+
+ self.tms['ElmLoop'].stop()
+
+ dummyLower = lowerNodes.copy()
+ dummyLower[:,0] *= -1
+
+ invStruct_tr[0] = np.vstack((upperNodes, dummyLower))
+ invStruct[0] = np.vstack((upperNodes, lowerNodes))
- # Create node list
- nodesCoord = np.zeros((0, 3))
- for n in bnd_wing.nodes:
- nodesCoord = np.vstack((nodesCoord, [n.pos[0], n.pos[1], n.pos[2]]))
+ # Remove duplicates
+ invStruct_tr[0] = np.unique(invStruct_tr[0], axis=0)
+ # Wake part
+ for iNode, n in enumerate(blw[1].nodes):
+ invStruct_tr[1][iNode, :] = [n.pos[0], n.pos[1], n.pos[2], n.row]
+ invStruct[1] = invStruct_tr[1].copy()
+
+ # Create spanwise distribution
if config['nDim'] == 2:
spanDistri = np.zeros(1)
- spanDistri[0] = nodesCoord[0,2] # Span direction is the z direction
+ spanDistri[0] = invStruct[0][0,2] # Span direction is the z direction
+ wingSpan = 0.
+ fInterp = interp1d(invStruct_tr[0][:,0], invStruct_tr[0][:,1])
-
elif config['nDim'] == 3:
- nodesCoord[:, [1,2]] = nodesCoord[:, [2, 1]]
- spanDistri = np.linspace(min(nodesCoord[:,2]), max(nodesCoord[:,2]), config['nSections'])
+ invStruct[0][:, [1,2]] = invStruct[0][:, [2, 1]]
+ invStruct[1][:, [1,2]] = invStruct[1][:, [2, 1]]
+ invStruct_tr[0][:, [1,2]] = invStruct_tr[0][:, [2, 1]]
+ invStruct_tr[1][:, [1,2]] = invStruct_tr[1][:, [2, 1]]
+ spanDistri = np.linspace(min(invStruct[0][:,2]), max(invStruct[0][:,2]), config['nSections'])
spanDistri[0] +=0.1
spanDistri[-1] -=0.1
- print('spanDistri', spanDistri)
-
- else:
- raise RuntimeError("Number of dimensions not recognized", config['nDim'])
-
- fig = plt.figure()
- ax = fig.add_subplot(projection='3d')
- ax.scatter(nodesCoord[:,0], nodesCoord[:,1], nodesCoord[:,2], s=0.3)
- #ax.scatter(wake[:,0], wake[:,1], wake[:,2], s=0.3)
- #ax.set_zlim([-0.5, 0.5])
- ax.set_xlabel('x')
- ax.set_ylabel('y')
- ax.set_zlabel('z')
- plt.show()
-
- wingSpan = max(nodesCoord[:,2]) - min(nodesCoord[:,2])
- upperNodes, lowerNodes = self.__SeparateWing(nodesCoord)
+ print('spanDistri', spanDistri)
+
+ # Spline wing
+ # Point in the wing tip region are not accounted for during Bezier curve spline
+ #xUp = upperNodes[upperNodes[:,0]<config['span'], :]
+ #xLw = lowerNodes[lowerNodes[:,0]<config['span'], :]
- lowerNodes[:,0] *= -1
- transformedCoord = np.vstack((upperNodes, lowerNodes))
+ #xUp[:,[1,2]] = xUp[:,[2,1]]
+ #xLw[:,[1,2]] = xLw[:,[2,1]]
- tck_wing = bisplrep(transformedCoord[:,0], transformedCoord[:,2], transformedCoord[:,1], kx=3, ky=3, s=1e-3, quiet=1)
+ #sUp = (len(xUp)-np.sqrt(2*len(xUp)) + len(xUp)+np.sqrt(2*len(xUp)))/2
+ #sLw = (len(xLw)-np.sqrt(2*len(xLw)) + len(xLw)+np.sqrt(2*len(xLw)))/2
+ #tck_wingUp = bisplrep(xUp[:,0], xUp[:,2], xUp[:,1], kx=3, ky=3, s=sUp, quiet=1)
+ #tck_wingLw = bisplrep(xLw[:,0], xLw[:,2], xLw[:,1], kx=3, ky=3, s=sLw, quiet=1)
- wakeCoord = np.zeros((len(bnd_wake.nodes), 3))
- for iNode, n in enumerate(bnd_wake.nodes):
- wakeCoord[iNode, 0] = n.pos[0]
- wakeCoord[iNode, 1] = n.pos[1]
- wakeCoord[iNode, 2] = n.pos[2]
-
- if config['nDim'] == 3:
- wakeCoord[:, [1,2]] = wakeCoord[:, [1,2]]
+ m = len(invStruct_tr[0][:,0])
+ _s = (m-np.sqrt(2*m) + m+np.sqrt(2*m))/2
+ tck_wing = bisplrep(invStruct_tr[0][invStruct_tr[0][:,2]<config['span'],0], invStruct_tr[0][invStruct_tr[0][:,2]<config['span'],2], invStruct_tr[0][invStruct_tr[0][:,2]<config['span'],1], kx=5, ky=5, s=1e-8)
- #tck_wake = splrep(wakeCoord)
+ # Create sections
+ self.tms['iSecLoop'].start()
for iSec, z in enumerate(spanDistri):
-
- portion = nodesCoord[abs(nodesCoord[:,2] - z) <= 0.01 * wingSpan,:]
+ portion = invStruct[0][abs(invStruct[0][:,2] - z) <= 0.01 * config['span'],:]
chordLength = max(portion[:, 0]) - min(portion[:, 0])
xDistri = min(portion[:,0]) + 1/2*(1 + np.cos(np.linspace(0, np.pi, config['nPoints'][0])))*chordLength
@@ -327,180 +313,115 @@ class Interpolator:
viscStruct_tr[iSec][0][:,2] = z*np.ones(len(viscStruct[iSec][0][:,0]))
# Interpolation in the transformed space.
- for iPoint in range(len(viscStruct_tr[iSec][0][:,0])):
- viscStruct_tr[iSec][0][iPoint,1] = bisplev(viscStruct_tr[iSec][0][iPoint,0], viscStruct_tr[iSec][0][iPoint,2], tck_wing)
-
+ if config['nDim'] == 3:
+
+ for iPoint in range(len(viscStruct_tr[iSec][0][:,0])):
+ viscStruct_tr[iSec][0][iPoint, 1] = bisplev(viscStruct_tr[iSec][0][iPoint,0], viscStruct_tr[iSec][0][iPoint,2], tck_wing)
+
+
+ elif config['nDim'] == 2:
+ viscStruct_tr[iSec][0][:,1] = fInterp(viscStruct_tr[iSec][0][:,0])
+
+ viscStruct_tr[iSec][0][0,1] = 0.
+ viscStruct_tr[iSec][0][-1,1] = 0.
viscStruct[iSec][0][:,1] = viscStruct_tr[iSec][0][:,1].copy()
- fig = plt.figure()
- ax = fig.add_subplot(projection='3d')
- ax.scatter(transformedCoord[:,0], transformedCoord[:,2], transformedCoord[:,1], s=0.3)
- #ax.scatter(wake[:,0], wake[:,1], wake[:,2], s=0.3)``
- for iSec in range(len(viscStruct)):
- ax.scatter(viscStruct_tr[iSec][0][:,0], viscStruct_tr[iSec][0][:,2], viscStruct[iSec][0][:,1], color='red')
- ax.set_zlim([-0.5, 0.5])
- ax.set_xlabel('x')
- ax.set_ylabel('y')
- ax.set_zlabel('z')
- plt.show()
-
- #viscStruct_tr[iSec][0][:,1] = splev(viscStruct_tr[iSec][0][:,0], spl)
- #viscStruct_tr[iSec][0][:,1] = interp1d(transformedCoord[:,0], transformedCoord[:,1])(viscStruct_tr[iSec][0][:,0])
- #viscStruct[iSec][0][:,1] = viscStruct_tr[iSec][0][:,1].copy()
+ # Wake
+ wakeLength = np.max(invStruct[1][:,0]) - np.max(viscStruct[iSec][0][:,0])
+ viscStruct_tr[iSec][1][:,0] = np.max(viscStruct[iSec][0][:,0]) + (np.max(viscStruct[iSec][0][:,0]) + np.cos(np.linspace(np.pi, np.pi/2, config['nPoints'][1]))) * wakeLength
+ viscStruct_tr[iSec][1][:,2] = z*np.ones(len(viscStruct[iSec][1][:,0]))
+ viscStruct_tr[iSec][1][:,1] = interp1d(invStruct_tr[1][:,0], invStruct_tr[1][:,1], bounds_error=False, fill_value='extrapolate')(viscStruct_tr[iSec][1][:,0])
+
+ viscStruct[iSec][1] = viscStruct_tr[iSec][1].copy()
- """# Wake
- wakeLength = np.max(wakeCoord[:,0]) - np.max(viscStruct[iSec][0][:,0])
- viscStruct[iSec][1][:,0] = np.max(viscStruct[iSec][0][:,0]) + (np.max(viscStruct[iSec][0][:,0]) + np.cos(np.linspace(np.pi, np.pi/2, config['nPoints'][1]))) * wakeLength
- viscStruct[iSec][1][:,2] = z*np.ones(len(viscStruct[iSec][1][:,0]))
- viscStruct[iSec][1][:,1] = interp1d(wakeCoord[:,0], wakeCoord[:,1])(viscStruct[iSec][1][:,0])"""
+ # Cg positions
+ self.tms['iSecLoop'].stop()
fig = plt.figure()
ax = fig.add_subplot(projection='3d')
- ax.scatter(nodesCoord[:,0], nodesCoord[:,2], nodesCoord[:,1], s=0.3)
- #ax.scatter(wake[:,0], wake[:,1], wake[:,2], s=0.3)
- for iSec in range(len(viscStruct)):
- ax.scatter(viscStruct[iSec][0][:,0], viscStruct[iSec][0][:,2], viscStruct[iSec][0][:,1], color='red')
+ ax.scatter(invStruct_tr[0][:,0], invStruct_tr[0][:,2], invStruct_tr[0][:,1], s=0.3)
+ ax.scatter(invStruct_tr[1][:,0], invStruct_tr[1][:,2], invStruct_tr[1][:,1], s=0.3)
+ 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('z')
plt.show()
- return nodesList, transformedCoord, wakeCoord, viscStruct, viscStruct_tr
-
-
- def __GetAirfoil(self, bnd_wing, bnd_wake, config):
-
- if config['nSections'] != 1:
- raise RuntimeError('Cannot create more than 1 section on 2D geometry. Specified:', config['nSections'])
-
- viscStruct = [[np.zeros((config['nPoints'][1] if i==1 else 2*config['nPoints'][0]-1, 3)) for i in range(2)] for _ in range(config['nSections'])]
- viscStruct_tr = [[np.zeros((config['nPoints'][1] if i==1 else 2*config['nPoints'][0]-1, 3)) for i in range(1)] for _ in range(config['nSections'])]
-
- # Create node list
- nodesCoord = np.zeros((0, 3))
- for n in bnd_wing.nodes:
- nodesCoord = np.vstack((nodesCoord, [n.pos[0], n.pos[1], n.pos[2]]))
+ """fig = plt.figure()
+ ax = fig.add_subplot(projection='3d')
+ ax.scatter(invStruct[0][:,0], invStruct[0][:,2], invStruct[0][:,1], s=0.3)
+ for sec in viscStruct:
+ 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('z')
+ plt.show()"""
- if config['nDim'] == 2:
- spanDistri = np.zeros(1)
- spanDistri[0] = nodesCoord[0,2] # Span direction is the z direction
+ viscCg_tr = [[np.zeros((len(viscStruct[iSec][iReg][:,0]) - 1, 3)) for iReg in range(2)] for iSec in range(len(viscStruct))]
+ invCg_tr = [np.zeros((len(blw[iReg].tag.elems), 3)) for iReg in range(2)]
+ self.tms['CgLoop'].start()
- elif config['nDim'] == 3:
- nodesCoord[:, [1,2]] = nodesCoord[:, [2, 1]]
- spanDistri = np.linspace(min(nodesCoord[:,2]), max(nodesCoord[:,2]), config['nSections'])
- spanDistri[0] +=0.1
- spanDistri[-1] -=0.1
- print('spanDistri', spanDistri)
-
- else:
- raise RuntimeError("Number of dimensions not recognized", config['nDim'])
-
- wingSpan = max(nodesCoord[:,2]) - min(nodesCoord[:,2])
-
- upperNodes, lowerNodes = self.__SeparateWing(nodesCoord)
-
- for iSec, z in enumerate(spanDistri):
-
- portion = nodesCoord[abs(nodesCoord[:,2] - z) <= 0.01 * wingSpan,:]
-
- le_coord = portion[portion[:,0] == np.min((portion[:,0])), :][0]
- te_coord = portion[portion[:,0] == np.max((portion[:,0])), :][0]
-
- # Find te elements
- te_elems = []
- for e in bnd_wing.tag.elems:
+ # Viscous
+ for iSec in range(len(viscStruct)):
+ for iReg in range(len(viscStruct[iSec])):
+ for iElm in range(len(viscStruct[iSec][iReg][:,0]) - 1):
+ for iDir in range(3):
+ viscCg_tr[iSec][iReg][iElm, iDir] = viscStruct_tr[iSec][iReg][iElm, iDir] + viscStruct_tr[iSec][iReg][iElm + 1, iDir]
+ viscCg_tr[iSec][iReg] /= 2
+
+ # Inviscid
+ for iReg in range(len(blw)):
+ for iElm, e in enumerate(blw[iReg].tag.elems):
for n in e.nodes:
- if n.pos[0] == te_coord[0]:
- te_elems.append(e)
- break
-
- # Cg of te elements
- cg_teElems = np.zeros((len(te_elems), 3))
- for iElm, eTe in enumerate(te_elems):
- for n in eTe.nodes:
- cg_teElems[iElm, 0] += n.pos[0]
- cg_teElems[iElm, 1] += n.pos[1]
- cg_teElems[iElm, 2] += n.pos[2]
-
- # Sort upper side in 0 and lower side in 1
- if cg_teElems[0, 1] < cg_teElems[1, 1]:
- save = te_elems[0].copy()
- te_elems[0] = te_elems[1]
- te_elems[1] = save
-
- # Find te_nodes (sorted [upper, lower])
- te_nodes = [None, None]
- for i in range(len(te_elems)):
- for n in te_elems[i].nodes:
- if n.pos[0] == te_coord[0]:
- te_nodes[i] = n
- break
-
- # Find le_node
- for n in bnd_wing.nodes:
- if n.pos[0] == le_coord[0]:
- le_node = n
- break
-
+ if n.row in upperNodes[:,3] or iReg == 1:
+ invCg_tr[iReg][iElm, 0] +=n.pos[0]
+ elif n.row in lowerNodes[:,3]:
+ invCg_tr[iReg][iElm, 0] += -n.pos[0]
- meany = (le_coord[1] + te_coord[1]) / 2
-
- transformedCoord = np.zeros((0, 3))
- nodesList = []
-
- for iNode, n in enumerate(bnd_wing.nodes):
- if n.pos[0] != le_coord[0] and n.pos[0] != te_coord[0]:
- if n.pos[1] > meany:
- transformedCoord = np.vstack((transformedCoord, [n.pos[0], n.pos[1], n.pos[2]]))
- nodesList.append(n)
-
- elif n.pos[1] < meany:
- transformedCoord = np.vstack((transformedCoord, [-n.pos[0], n.pos[1], n.pos[2]]))
- nodesList.append(n)
-
- transformedCoord = np.vstack((transformedCoord, le_coord))
- nodesList.append(le_node)
- transformedCoord = np.vstack((transformedCoord, [te_nodes[0].pos[0], te_nodes[0].pos[1], te_nodes[0].pos[2]]))
- nodesList.append(te_nodes[0])
- transformedCoord = np.vstack((transformedCoord, [-te_nodes[1].pos[0], te_nodes[1].pos[1], te_nodes[1].pos[2]]))
- nodesList.append(te_nodes[1])
-
- wakeCoord = np.zeros((len(bnd_wake.nodes), 3))
- for iNode, n in enumerate(bnd_wake.nodes):
- wakeCoord[iNode, 0] = n.pos[0]
- wakeCoord[iNode, 1] = n.pos[1]
- wakeCoord[iNode, 2] = n.pos[2]
-
- chordLength = max(portion[:, 0]) - min(portion[:, 0])
- xDistri = le_coord[0] + 1/2*(1 + np.cos(np.linspace(0, np.pi, config['nPoints'][0])))*chordLength
+ invCg_tr[iReg][iElm, 1] +=n.pos[1]
+ invCg_tr[iReg][iElm, 2] +=n.pos[2]
+ invCg_tr[iReg][iElm, :] /= e.nodes.size()
- viscStruct[iSec][0][:,0] = np.concatenate((xDistri[:-1], xDistri[::-1]), axis=None)
- viscStruct[iSec][0][:,2] = z*np.ones(len(viscStruct[iSec][0][:,0]))
-
- viscStruct_tr[iSec][0][:,0] = np.concatenate((xDistri[:-1], -xDistri[::-1]), axis=None).copy()
- viscStruct_tr[iSec][0][:,2] = z*np.ones(len(viscStruct[iSec][0][:,0]))
-
- # Interpolation in the transformed space.
- splineTr = transformedCoord.copy()
- splineTr = splineTr[splineTr[:, 0].argsort()]
- spl = splrep(splineTr[:,0], splineTr[:,1], k=5)
- viscStruct_tr[iSec][0][:,1] = splev(viscStruct_tr[iSec][0][:,0], spl)
- #viscStruct_tr[iSec][0][:,1] = interp1d(transformedCoord[:,0], transformedCoord[:,1])(viscStruct_tr[iSec][0][:,0])
- viscStruct[iSec][0][:,1] = viscStruct_tr[iSec][0][:,1].copy()
+ if config['nDim'] == 3:
+ invCg_tr[iReg][:,[1,2]] = invCg_tr[iReg][:,[2,1]]
+ self.tms['CgLoop'].stop()
- # Wake
- wakeLength = np.max(wakeCoord[:,0]) - np.max(viscStruct[iSec][0][:,0])
- viscStruct[iSec][1][:,0] = np.max(viscStruct[iSec][0][:,0]) + (np.max(viscStruct[iSec][0][:,0]) + np.cos(np.linspace(np.pi, np.pi/2, config['nPoints'][1]))) * wakeLength
- viscStruct[iSec][1][:,2] = z*np.ones(len(viscStruct[iSec][1][:,0]))
- viscStruct[iSec][1][:,1] = interp1d(wakeCoord[:,0], wakeCoord[:,1])(viscStruct[iSec][1][:,0])
+ """if config['nDim'] == 3:
+ fig = plt.figure()
+ ax = fig.add_subplot(projection='3d')
+ ax.scatter(invStruct_tr[0][:,0], invStruct_tr[0][:,2], invStruct[0][:,1], s=0.3)
+ #ax.scatter(abs(invCg_tr[0][:,0]), invCg_tr[0][:,2], invCg_tr[0][:,1], s=0.3)
+ for sec in viscStruct_tr:
+ ax.scatter(abs(sec[0][:,0]), sec[0][:,2], sec[0][:,1], color = 'red')
+ ax.scatter(abs(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('z')
+ plt.show()
+ elif config['nDim'] == 2:
+ plt.plot(abs(invStruct_tr[0][:,0]), invStruct[0][:,1], 'x', color='blue')
+ plt.plot(abs(invStruct_tr[1][:,0]), invStruct[1][:,1], 'x', color='blue')
+ plt.plot(abs(invCg_tr[0][:,0]), invCg_tr[0][:,1], 'o', color='red')
+ plt.plot(abs(invCg_tr[1][:,0]), invCg_tr[1][:,1], 'o', color='red')
+ plt.title('Inviscid')
+ plt.show()
+ plt.plot(viscStruct_tr[0][0][:,0], viscStruct_tr[0][0][:,1])
+ plt.plot(viscStruct_tr[0][1][:,0], viscStruct_tr[0][1][:,1])
+ plt.plot(invStruct_tr[0][:,0], invStruct[0][:,1], 'o', color='blue')
+ plt.plot(invStruct_tr[1][:,0], invStruct[1][:,1], 'x', color='blue')
+ plt.title('Viscous vs Inviscid')
+ plt.show()"""
- return nodesList, transformedCoord, wakeCoord, viscStruct, viscStruct_tr
- def __SeparateWing(self, nodesCoord):
-
- print(np.max(nodesCoord[:,0]))
- meany = nodesCoord[np.argmax(nodesCoord[:,0]), 1] - nodesCoord[np.argmin(nodesCoord[:,0]), 1]
+ self.invStruct = invStruct
+ self.invStruct_tr = invStruct_tr
+ self.viscStruct = viscStruct
+ self.viscStruct_tr = viscStruct_tr
- upperNodes = nodesCoord[nodesCoord[:,1]>=meany].copy()
- lowerNodes = nodesCoord[nodesCoord[:,1]<=meany].copy()
- return upperNodes, lowerNodes
\ No newline at end of file
+ self.invCg_tr = invCg_tr
+ self.viscCg_tr = viscCg_tr
\ No newline at end of file
diff --git a/dart/pyVII/vStructures.py b/dart/pyVII/vStructures.py
deleted file mode 100644
index fd6f7f5..0000000
--- a/dart/pyVII/vStructures.py
+++ /dev/null
@@ -1,46 +0,0 @@
-import numpy as np
-
-class viscStruct:
- def __init__(self, nPoints_surf, nPoints_wake):
- self.nSurfNodes = nPoints_surf
- self.nWakeNodes = nPoints_wake
-
- self.nodesReal = np.zeros((self.nSurfNodes, 3)) # Coordinates of surface nodes in the real space
- self.nodesTrans = np.zeros((self.nWakeNodes, 3)) # Coordinates of surface nodes in the transformed space
-
- # Usefull variables for section representation
- self.chordLength = 0.
- self.leCoord = 0.
-
-
- self.nodesWake = np.zeros((nPoints_wake, 3))
-
- def GenerateDistri(self):
- xDistri = self.le_coord + 1/2*(1 + np.cos(np.linspace(0, np.pi, len(self.nSurfNodes))))*self.chordLength
- self.nodesReal[:,0] = np.concatenate((xDistri[:-1], xDistri[::-1]), axis=None)
-
-class invStruct:
- def __init__(self, _bndSurf, _bndWake, _nDim):
- self.nDim = _nDim
-
- self.bndSurf = _bndSurf
- self.bndWake = _bndWake
-
- self.nSurfNodes = self.bndSurf.nodes.size()
-
- # np matrix of nodes coordinate in the transformed space associated to the c++ structures 'node' in nodesList
- # (coordinates at indice i in 'nodesTrans' correspond to the node at position i in 'nodesList')
- self.nodesTrans = np.zeros(self.nSurfNodes,3)
- self.nodesList = [None for _ in range(self.nSurfNodes)]
-
- def CreateList(self):
-
- nodesReal = np.zeros((self.nSurfNodes, 3))
- nodeList = [None for _ in range(self.nSurfNodes)]
- for iNode, n in enumerate(self.bndSurf.nodes):
- nodesReal[iNode, :] = [n.pos[0], n.pos[1], n.pos[2]]
- nodeList[iNode] = n
-
-
-
-
diff --git a/dart/pyVII/vUtils.py b/dart/pyVII/vUtils.py
index a8423fd..7009f8d 100644
--- a/dart/pyVII/vUtils.py
+++ b/dart/pyVII/vUtils.py
@@ -50,36 +50,7 @@ def GetSolution(iSec, vSolver):
'Cdp' : vSolver.Cdp}
return sln
-
-def Dictionary(blScalars, blVectors, xtr):
- """Create dictionary with the coordinates and the boundary layer parameters
- """
-
- wData = {}
- # [xx, x, deltaStar, theta, H, N, Ue, Hk, Hstar, Hstar2, Cf, Cd, Ct, delta, VtExt, Me, Rhoe]
- wData['xx'] = blVectors[0]
- wData['x/c'] = blVectors[1]
- wData['deltaStar'] = blVectors[2]
- wData['theta'] = blVectors[3]
- wData['H'] = blVectors[4]
- wData['N'] = blVectors[5]
- wData['Ue'] = blVectors[6]
- wData['Ct'] = blVectors[7]
- wData['Hk'] = blVectors[8]
- wData['Hstar'] = blVectors[9]
- wData['HStar2'] = blVectors[10]
- wData['Cf'] = blVectors[11]
- wData['Cd'] = blVectors[12]
- wData['delta'] = blVectors[13]
- wData['VtExt'] = blVectors[14]
- wData['Me'] = blVectors[15]
- wData['Rhoe'] = blVectors[16]
- wData['Blowing'] = blVectors[17]
- wData['xtrTop'] = xtr[0]
- wData['xtrBot'] = xtr[1]
- return wData
-
-def WriteFile(wData, Re, toW=['delta*', 'H', 'Cf', 'Ctau']):
+def WriteFile(wData, Re, toW=['deltaStar', 'H', 'Cf', 'Ct', 'delta']):
"""Write the results in airfoil_viscous.dat
"""
@@ -89,18 +60,18 @@ def WriteFile(wData, Re, toW=['delta*', 'H', 'Cf', 'Ctau']):
f.write('$Aerodynamic coefficients\n')
f.write(' Re Cd Cdp Cdf xtr_top xtr_bot\n')
- f.write('{0:15.6f} {1:15.6f} {2:15.6f} {3:15.6f} {4:15.6f} {5:15.6f}\n'.format(Re, wData['Cd'], wData['Cdp'],
- wData['Cdf'], wData['xtrTop'],
- wData['xtrBot']))
+ f.write('{0:15.6f} {1:15.6f} {2:15.6f} {3:15.6f} {4:15.6f} {5:15.6f}\n'.format(Re, wData['Cdt'], wData['Cdp'],
+ wData['Cdf'], wData['xtrT'],
+ wData['xtrB']))
f.write('$Boundary layer variables\n')
- f.write(' x y z x/c')
+ f.write(' x/c')
for s in toW:
f.write(' {0:>15s}'.format(s))
f.write('\n')
- for i in range(len(self.data[:,0])):
- f.write('{0:15.6f} {1:15.6f} {2:15.6f} {3:15.6f}'.format(wData['x'][i], wData['y'][i], wData['z'][i], wData['x/c'][i]))
+ for i in range(len(wData['x/c'])):
+ f.write('{0:15.6f}'.format(wData['x/c'][i]))
for s in toW:
f.write(' {0:15.6f}'.format(wData[s][i]))
f.write('\n')
@@ -108,17 +79,6 @@ def WriteFile(wData, Re, toW=['delta*', 'H', 'Cf', 'Ctau']):
f.close()
print('done.')
-def ExtractVars(viscSolver):
- import numpy as np
-
- blScalars = [viscSolver.Cdt, viscSolver.Cdf, viscSolver.Cdp]
-
- # [xx, x, deltaStar, theta, H, N, Ue, Hk, Hstar, Hstar2, Cf, Cd, Ct, delta, VtExt, Me, Rhoe, Blowing]
-
- blVectors = viscSolver.ExtractSolution()
-
-
- return blScalars, blVectors
def PlotVars(wData, var='H'):
from matplotlib import pyplot as plt
diff --git a/dart/src/wBody.cpp b/dart/src/wBody.cpp
index 511e9be..f81ef6d 100644
--- a/dart/src/wBody.cpp
+++ b/dart/src/wBody.cpp
@@ -41,6 +41,10 @@ Body::Body(std::shared_ptr<MshData> _msh, std::vector<int> const &nos) : Groups(
Body::Body(std::shared_ptr<MshData> _msh, std::vector<std::string> const &names) : Groups(_msh, names), Cl(0), Cd(0), Cs(0), Cm(0)
{
+ for (auto i=0; i<names.size(); ++i)
+ {
+ std::cout << names[i] << std::endl;
+ }
// Sanity check
if (groups.size() != 2)
{
diff --git a/dart/src/wClosures.cpp b/dart/src/wClosures.cpp
index f9874c4..7f867be 100644
--- a/dart/src/wClosures.cpp
+++ b/dart/src/wClosures.cpp
@@ -103,14 +103,6 @@ std::vector<double> Closures::LaminarClosures(double theta, double H, double Ret
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;
}
@@ -135,15 +127,15 @@ std::vector<double> Closures::TurbulentClosures(double theta, double H, double C
}
double Hstar2 = ((0.064/(Hk-0.8))+0.251)*Me*Me;
- double gamma = 1.4;
- double Fc = std::sqrt(1+0.5*gamma*Me*Me);
+ double gamma = 1.4 - 1.;
+ double Fc = std::sqrt(1+0.5*gamma*Me*Me); /* .2 can be used as well replacing .5*gamma */
double H0 = 0.;
if (Ret <= 400)
{
H0 = 4.0;
}
- else if (Ret > 400)
+ else
{
H0 = 3 + (400/Ret);
}
@@ -157,7 +149,7 @@ std::vector<double> Closures::TurbulentClosures(double theta, double H, double C
{
Hstar = ((0.5-4/Ret)*((H0-Hk)/(H0-1))*((H0-Hk)/(H0-1)))*(1.5/(Hk+0.5))+1.5+4/Ret;
}
- else if (Hk > H0)
+ else
{
Hstar = (Hk-H0)*(Hk-H0)*(0.007*std::log(Ret)/((Hk-H0+4/std::log(Ret))*(Hk-H0+4/std::log(Ret))) + 0.015/Hk) + 1.5 + 4/Ret;
}
@@ -166,11 +158,6 @@ 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. Ret= " << Ret << std::endl;
- //throw;
- }
logRt = std::max(logRt, 3.0);
double arg = std::max(-1.33*Hk, -20.0);
@@ -195,6 +182,8 @@ std::vector<double> Closures::TurbulentClosures(double theta, double H, double C
double Cd = 0.;
double n = 1.0;
+
+ double Cteq = 0.;
if (name == "wake")
{
@@ -208,6 +197,7 @@ std::vector<double> Closures::TurbulentClosures(double theta, double H, double C
Cdw = 0.0; // wall contribution
Cdd = (0.995-Us)*Ct*Ct; // turbulent outer layer contribution
Cdl = 0.15*(0.995-Us)*(0.995-Us)/Ret; // laminar stress contribution to outer layer
+ Cteq = std::sqrt(4*Hstar*Ctcon*((Hk-1)*Hkc*Hkc)/((1-Us)*(Hk*Hk)*H)); // Here it is Cteq^0.5
}
else
{
@@ -225,13 +215,14 @@ std::vector<double> Closures::TurbulentClosures(double theta, double H, double C
Cdw = 0.5*(Cf*Us) * Dfac;
Cdd = (0.995-Us)*Ct*Ct;
Cdl = 0.15*(0.995-Us)*(0.995-Us)/Ret;
+ Cteq = std::sqrt(Hstar*Ctcon*((Hk-1)*Hkc*Hkc)/((1-Us)*(Hk*Hk)*H)); // Here it is Cteq^0.5
}
if (n*Hstar*Ctcon*((Hk-1)*Hkc*Hkc)/((1-Us)*(Hk*Hk)*H)<0)
{
std::cout << "Negative sqrt encountered " << std::endl;
//throw;
}
- double Cteq = std::sqrt(n*Hstar*Ctcon*((Hk-1)*Hkc*Hkc)/((1-Us)*(Hk*Hk)*H));
+ //double Cteq = std::sqrt(n*Hstar*Ctcon*((Hk-1)*Hkc*Hkc)/((1-Us)*(Hk*Hk)*H));
Cd = n*(Cdw+ Cdd + Cdl);
// Ueq = 1.25/Hk*(Cf/2-((Hk-1)/(6.432*Hk))**2)
@@ -246,13 +237,5 @@ std::vector<double> Closures::TurbulentClosures(double theta, double H, double C
ClosureVars[6] = Cteq;
ClosureVars[7] = Us;
- 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/wViscFluxes.cpp b/dart/src/wViscFluxes.cpp
index 15daf4a..0346a58 100644
--- a/dart/src/wViscFluxes.cpp
+++ b/dart/src/wViscFluxes.cpp
@@ -129,7 +129,6 @@ VectorXd ViscFluxes::BLlaws(size_t iPoint, BLRegion *bl, std::vector<double> u)
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);
diff --git a/dart/src/wViscSolver.cpp b/dart/src/wViscSolver.cpp
index a3d0f9b..2a06c2c 100644
--- a/dart/src/wViscSolver.cpp
+++ b/dart/src/wViscSolver.cpp
@@ -268,6 +268,7 @@ int ViscSolver::Run(unsigned int const couplIter)
if (verbose>2){"Blowing velocities of " + Sections[iSec][iRegion]->name + " side computed.";}
} // End for iRegion
+ if (verbose>0){std::cout << "\n";}
} // End for iSections
ComputeDrag(Sections[0]);
@@ -487,6 +488,7 @@ void ViscSolver::SetDeltaStarExt(size_t iSec, size_t iRegion, std::vector<double
{
if (_DeltaStarExt.size() != Sections[iSec][iRegion]->mesh->GetnMarkers())
{
+ std::cout << "size DeltaStar: " << _DeltaStarExt.size() << ". nMarkers: " << Sections[iSec][iRegion]->mesh->GetnMarkers() << std::endl;
throw std::runtime_error("dart::ViscSolver External delta star vector is not coherent with the mesh.");
}
Sections[iSec][iRegion]->invBnd->SetDeltaStar(_DeltaStarExt);
diff --git a/dart/tests/bli2.py b/dart/tests/bli2.py
new file mode 100644
index 0000000..26c8a75
--- /dev/null
+++ b/dart/tests/bli2.py
@@ -0,0 +1,156 @@
+#!/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) viscous flow around a NACA 0012
+# Amaury Bilocq
+#
+# Test the viscous-inviscid interaction scheme
+# Reference to the master's thesis: http://hdl.handle.net/2268/252195
+# Reference test cases with Naca0012 (different from master's thesis):
+# 1) Incompressible: Re = 1e7, M_inf = 0, alpha = 5°, msTE = 0.01, msLE = 0.01
+# -> nIt = 27, Cl = 0.55, Cd = 0.0058, xtrTop = 0.087, xtrBot = 0.741
+# -> msLe = 0.001, xtrTop = 0.061
+# 2) Compressible: Re = 1e7, M_inf = 5, alpha = 5°, msTE = 0.01, msLE = 0.01
+# -> nIt = 33, Cl = 0.65, Cd = 0.0063, xtrTop = 0.057, xtrBot = 0.740
+# 3) Separated: Re = 1e7, M_inf = 0, alpha = 12°, msTE = 0.01, msLE = 0.001
+# -> nIt = 43, Cl = 1.30 , Cd = 0.0108, xtrTop = 0.010, xtrBot = 1
+#
+# CAUTION
+# This test is provided to ensure that the solver works properly.
+# Mesh refinement may have to be performed to obtain physical results.
+
+import dart.utils as floU
+import dart.default as floD
+#import dart.viscUtils as viscU
+
+import dart.pyVII.vCoupler as floVC
+import dart.pyVII.vUtils as viscU
+import numpy as np
+
+import tbox
+import tbox.utils as tboxU
+import fwk
+from fwk.testing import *
+from fwk.coloring import ccolors
+from matplotlib import pyplot as plt
+
+import math
+
+def main():
+ print('Start')
+ # timer
+ tms = fwk.Timers()
+ tms['total'].start()
+
+ # define flow variables
+ Re = 1e7
+ alpha = 2*math.pi/180
+ M_inf = 0.715
+
+ CFL0 = 1
+
+ nSections = 1
+
+ # ========================================================================================
+
+ c_ref = 1
+ dim = 2
+
+ # 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.005}
+ msh, gmshWriter = floD.mesh(dim, 'models/rae2822.geo', pars, ['field', 'airfoil', 'downstream'])
+ tms['msh'].stop()
+
+ # set the problem and add medium, initial/boundary conditions
+ tms['pre'].start()
+ 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)
+ tms['pre'].stop()
+
+ # solve viscous problem
+
+ print(ccolors.ANSI_BLUE + 'PySolving...' + ccolors.ANSI_RESET)
+ tms['solver'].start()
+
+ isolver = floD.newton(pbl)
+ vsolver = viscU.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}
+
+ coupler = floVC.Coupler(isolver, vsolver)
+ coupler.CreateProblem(blw[0], blw[1])
+ coupler.Run()
+
+ # extract Cp
+ Cp = floU.extract(bnd.groups[0].tag.elems, isolver.Cp)
+ tboxU.write(Cp, 'Cp.dat', '%1.5e', ', ', 'x, y, z, Cp', '')
+ # display results
+ print(ccolors.ANSI_BLUE + 'PyRes...' + ccolors.ANSI_RESET)
+ print(' Re M alpha Cl Cd Cdp Cdf Cm')
+ print('{0:6.1f}e6 {1:8.2f} {2:8.1f} {3:8.4f} {4:8.4f} {5:8.4f} {6:8.4f} {7:8.4f}'.format(Re/1e6, M_inf, alpha*180/math.pi, isolver.Cl, vsolver.Cdt, vsolver.Cdp, vsolver.Cdf, isolver.Cm))
+
+ # Write results to file
+ vSol = viscU.GetSolution(0, vsolver)
+ viscU.WriteFile(vSol, vsolver.GetRe())
+
+ # display timers
+ tms['total'].stop()
+ print(ccolors.ANSI_BLUE + 'PyTiming...' + ccolors.ANSI_RESET)
+ print('CPU statistics')
+ print(tms)
+ print('SOLVERS statistics')
+ print(coupler.tms)
+
+ xtr=vsolver.Getxtr()
+
+ """# visualize solution and plot results
+ #floD.initViewer(pbl)
+ tboxU.plot(Cp[:,0], Cp[:,3], 'x', 'Cp', 'Cl = {0:.{3}f}, Cd = {1:.{3}f}, Cm = {2:.{3}f}'.format(isolver.Cl, vsolver.Cdt, isolver.Cm, 4), True)
+
+ blScalars, blVectors = viscU.ExtractVars(vsolver)
+ wData = viscU.Dictionary(blScalars, blVectors, xtr)
+ viscU.PlotVars(wData, plotVar)"""
+
+
+ # check results
+ print(ccolors.ANSI_BLUE + 'PyTesting...' + ccolors.ANSI_RESET)
+ tests = CTests()
+ if Re == 1e7 and M_inf == 0. and alpha == 2.*math.pi/180:
+ tests.add(CTest('Cl', isolver.Cl, 0.2208, 5e-3))
+ tests.add(CTest('Cd', vsolver.Cdt, 0.00531, 1e-3, forceabs=True))
+ tests.add(CTest('Cdp', vsolver.Cdp, 0.0009, 1e-3, forceabs=True))
+ tests.add(CTest('xtr_top', xtr[0], 0.20, 0.03, forceabs=True))
+ tests.add(CTest('xtr_bot', xtr[1], 0.50, 0.03, forceabs=True))
+ if Re == 1e7 and M_inf == 0. and alpha == 5.*math.pi/180:
+ tests.add(CTest('Cl', isolver.Cl, 0.5488, 5e-3))
+ tests.add(CTest('Cd', vsolver.Cdt, 0.0062, 1e-3, forceabs=True))
+ tests.add(CTest('Cdp', vsolver.Cdp,0.0018, 1e-3, forceabs=True))
+ tests.add(CTest('xtr_top', xtr[0], 0.06, 0.03, forceabs=True))
+ tests.add(CTest('xtr_bot', xtr[1], 0.74, 0.03, forceabs=True))
+ else:
+ raise Exception('Test not defined for this flow')
+
+ tests.run()
+
+
+ # eof
+ print('')
+
+if __name__ == "__main__":
+ main()
diff --git a/dart/tests/bli3.py b/dart/tests/bli3.py
index 6ef0ad4..792d16b 100644
--- a/dart/tests/bli3.py
+++ b/dart/tests/bli3.py
@@ -30,7 +30,7 @@ import dart.utils as floU
import dart.default as floD
import dart.pyVII.vUtils as viscU
-import dart.pyVII.vCoupler as viscC
+import dart.pyVII.vCoupler2 as viscC
import dart.pyVII.vInterpolator as vInterpol
import fwk
@@ -43,9 +43,12 @@ def main():
tms['total'].start()
# define flow variables
- Re = 11.72e6
+ """Re = 11.72e6
alpha = 3.06*math.pi/180
- M_inf = 0.8395
+ M_inf = 0.8395"""
+ Re = 1e7
+ alpha = 0.*math.pi/180
+ M_inf = 0.
dim = 3
CFL0 = 1
@@ -74,16 +77,13 @@ def main():
tms['pre'].stop()
# solve problem
- config = {'nDim': dim, 'nSections':nSections, 'nPoints':[500, 25], 'eps':0.05, 'k':[5, 5]}
- iSolverAPI = vInterpol.Interpolator(floD.newton(pbl), blw[0], blw[1], config)
- vSolver = viscU.initBL(Re, M_inf, CFL0, nSections)
+ config = {'nDim': dim, 'nSections':nSections, 'span':spn, 'nPoints':[200, 25], 'eps':0.02, 'rbftype': 'cubic', 'smoothing': 0.001, 'degree': 1, 'neighbors': 5}
+ iSolverAPI = vInterpol.Interpolator(floD.newton(pbl), blw, config)
+ vSolver = viscU.initBL(Re, M_inf, CFL0, nSections, 2)
#iSolverAPI = viscAPI.dartAPI(floD.newton(pbl), bnd, wk, nSections, vInterp)
coupler = viscC.Coupler(iSolverAPI, vSolver)
- try:
- coupler.Run()
- except Exception as e:
- print('VII convergence terminated, ', str(e.args[0]))
+ coupler.Run()
# display timers
tms['total'].stop()
diff --git a/dart/tests/bli_lowLift.py b/dart/tests/bli_lowLift.py
index 544b06a..10c1ae6 100644
--- a/dart/tests/bli_lowLift.py
+++ b/dart/tests/bli_lowLift.py
@@ -59,7 +59,7 @@ def main():
# define flow variables
Re = 1e7
- alpha = 0*math.pi/180
+ alpha = 5*math.pi/180
M_inf = 0.
CFL0 = 1
@@ -71,7 +71,7 @@ def main():
c_ref = 1
dim = 2
- nRel = 5
+ nRel = 1
nomDeVariableDePorc_x = [None for _ in range(nRel)]
nomDeVariableDePorc_var = [None for _ in range(nRel)]
diff --git a/dart/tests/bli_lowLift2.py b/dart/tests/bli_lowLift2.py
index 7704020..05ab0ab 100644
--- a/dart/tests/bli_lowLift2.py
+++ b/dart/tests/bli_lowLift2.py
@@ -61,22 +61,22 @@ def main():
# define flow variables
Re = 1e7
- alpha = 2*math.pi/180
+ alpha = 5*math.pi/180
M_inf = 0.
CFL0 = 1
nSections = 1
+ verbose = 0
# ========================================================================================
c_ref = 1
dim = 2
- nPv = [100, 100, 100]
+ nPv = [200]
aeroCoeff = [np.zeros(2) for i in range(len(nPv))]
- fig, axs = plt.subplots(2,1)
for i in range(len(nPv)):
@@ -85,7 +85,7 @@ def main():
tms['msh'].start()
pars = {'xLgt' : 5, 'yLgt' : 5, 'msF' : 1, 'msTe' : 0.01, 'msLe' : 0.001}
- msh, gmshWriter = floD.mesh(dim, 'models/n0012.geo', pars, ['field', 'airfoil', 'downstream'])
+ msh, sortElm, gmshWriter = floD.mesh(dim, 'models/n0012.geo', pars, ['field', 'airfoil', 'downstream'])
tms['msh'].stop()
# set the problem and add medium, initial/boundary conditions
@@ -99,31 +99,24 @@ def main():
tms['solver'].start()
isolver = floD.newton(pbl)
- vsolver = viscU.initBL(Re, M_inf, CFL0, nSections, 0)
- config = {'nDim': dim, 'nSections':nSections, 'nPoints':[nPv[i], 25], 'eps':0.02, 'rbftype': 'cubic', 'smoothing': 0.0, 'degree': 1, 'neighbors': 3}
+ if verbose == 0:
+ isolver.printOn = False
+ vsolver = viscU.initBL(Re, M_inf, CFL0, nSections, verbose)
+ config = {'nDim': dim, 'nSections':nSections, 'nPoints':[nPv[i], 25], 'eps':0.02, 'rbftype': 'cubic', 'smoothing': 0.0, 'degree': 1, 'neighbors': 15}
- iSolverAPI = Interpol.Interpolator(isolver, blw[0], blw[1], config)
+ iSolverAPI = Interpol.Interpolator(isolver, sortElm, blw, config)
coupler = viscC.Coupler(iSolverAPI, vsolver)
aeroCoeff[i] = coupler.Run()
- axs[0].plot(aeroCoeff[i][:,0], 'x-')
- axs[1].plot(aeroCoeff[i][:,1], 'x-')
- del isolver
+ """del isolver
del vsolver
del iSolverAPI
del pbl
del bnd
del blw
del msh
- del gmshWriter
-
- axs[0].set_xlabel('x/c')
- axs[0].set_ylabel('Cl')
- axs[0].set_ylim([-.1, .1])
- axs[1].set_xlabel('x/c')
- axs[1].set_ylabel('Cd')
- plt.show()
+ del gmshWriter"""
tms['solver'].stop()
# extract Cp
diff --git a/dart/tests/rae_3.py b/dart/tests/rae_3.py
index f720dc4..f1f8add 100644
--- a/dart/tests/rae_3.py
+++ b/dart/tests/rae_3.py
@@ -62,14 +62,37 @@ def main():
print(ccolors.ANSI_BLUE + 'PyMeshing...' + ccolors.ANSI_RESET)
tms["msh"].start()
pars = {'spn' : spn, 'lgt' : lgt, 'wdt' : wdt, 'hgt' : hgt, 'msN' : nms, 'msF' : fms}
- msh, gmshWriter = floD.mesh(dim, 'models/rae_3.geo', pars, ['field', 'wing', 'symmetry', 'downstream'], wktp = 'wakeTip')
+ msh, sortElm, gmshWriter = floD.mesh(dim, 'models/rae_3.geo', pars, ['field', 'wing', 'symmetry', 'downstream'], wktp = 'wakeTip')
morpher = floD.morpher(msh, dim, ['wing'])
tms["msh"].stop()
# set the problem and add fluid, initial/boundary conditions
tms['pre'].start()
- pbl, _, _, bnd, _ = floD.problem(msh, dim, alpha, 0., M_inf, S_ref, c_ref, 0., 0., 0., 'wing', tp = 'wakeTip')
+ pbl, _, _, bnd, blw = floD.problem(msh, dim, alpha, 0., M_inf, S_ref, c_ref, 0., 0., 0., 'wing', tp = 'wakeTip', vsc=True)
+ upperNodes = np.zeros((0,4))
+ lowerNodes = np.zeros((0,4))
+
+ for e in blw[0].tag.elems:
+ if e.no in sortElm['upper']:
+ for n in e.nodes:
+ if n.row not in upperNodes[:,3]:
+ upperNodes = np.vstack((upperNodes, [n.pos[0], n.pos[1], n.pos[2], n.row]))
+ elif e.no in sortElm['lower']:
+ for n in e.nodes:
+ if n.row not in lowerNodes[:,3]:
+ lowerNodes = np.vstack((lowerNodes, [n.pos[0], n.pos[1], n.pos[2], n.row]))
+
+ fig = plt.figure()
+ ax = fig.add_subplot(projection='3d')
+ ax.scatter(upperNodes[:,0], upperNodes[:,1], upperNodes[:,2], s=0.3, color='red')
+ ax.scatter(lowerNodes[:,0], lowerNodes[:,1], lowerNodes[:,2], s=0.3, color='blue')
+ #ax.scatter(wake[:,0], wake[:,1], wake[:,2], s=0.3)
+ #ax.set_zlim([-0.5, 0.5])
+ ax.set_xlabel('x')
+ ax.set_ylabel('y')
+ ax.set_zlabel('z')
+ plt.show()
span = 0.5
elemPlot = []
for e in bnd.groups[0].tag.elems:
--
GitLab