From f8b942dffa30274e90f492f4413f569c98b80098 Mon Sep 17 00:00:00 2001
From: Dechamps Paul <paul.dechamps@student.uliege.be>
Date: Tue, 23 Nov 2021 00:59:49 +0100
Subject: [PATCH] Mesh adaptation fixed.
---
dart/tests/bliHighLift.py | 2 +-
dart/tests/bliNonLift.py | 4 +-
dart/tests/bliSeparated.py | 4 +-
dart/tests/bliTransonic.py | 2 +-
dart/viscous/Drivers/VDriver.py | 37 ++--
dart/viscous/Drivers/VGroupConstructor.py | 204 +++++++++++++---------
dart/viscous/Drivers/VPostProcess.py | 7 +-
dart/viscous/Physics/VClosures.py | 44 ++---
dart/viscous/Physics/VDataStructures.py | 36 ++--
dart/viscous/Physics/VVariables.py | 34 ++--
dart/viscous/Solvers/VSolver.py | 2 -
dart/viscous/Solvers/VTimeSolver.py | 18 +-
12 files changed, 224 insertions(+), 170 deletions(-)
diff --git a/dart/tests/bliHighLift.py b/dart/tests/bliHighLift.py
index c2029fa..900259c 100755
--- a/dart/tests/bliHighLift.py
+++ b/dart/tests/bliHighLift.py
@@ -66,7 +66,7 @@ def main():
user_xtr=[None,None]
user_Ti=None
- mapMesh = 0
+ mapMesh = 1
# Time solver parameters
Time_Domain=True # True for unsteady solver, False for steady solver
diff --git a/dart/tests/bliNonLift.py b/dart/tests/bliNonLift.py
index 23ac27a..456ab36 100755
--- a/dart/tests/bliNonLift.py
+++ b/dart/tests/bliNonLift.py
@@ -68,7 +68,7 @@ def main():
user_xtr=[None,None]
user_Ti=None
- mapMesh = 0
+ mapMesh = 1
# Time solver parameters
Time_Domain = True # True for unsteady solver, False for steady solver
@@ -89,7 +89,7 @@ def main():
# mesh the geometry
print(ccolors.ANSI_BLUE + 'PyMeshing...' + ccolors.ANSI_RESET)
tms['msh'].start()
- pars = {'xLgt' : 5, 'yLgt' : 5, 'msF' : 1, 'msTe' : 0.01, 'msLe' : 0.001}
+ pars = {'xLgt' : 5, 'yLgt' : 5, 'msF' : 1, 'msTe' : 0.01, 'msLe' : 0.01}
#pars = {'xLgt' : 5, 'yLgt' : 5, 'msF' : 1, 'msTe' : 0.002, 'msLe' : 0.0001}
msh, gmshWriter = floD.mesh(dim, 'models/n0012.geo', pars, ['field', 'airfoil', 'downstream'])
tms['msh'].stop()
diff --git a/dart/tests/bliSeparated.py b/dart/tests/bliSeparated.py
index fcae11d..1bc07a1 100755
--- a/dart/tests/bliSeparated.py
+++ b/dart/tests/bliSeparated.py
@@ -66,7 +66,7 @@ def main():
user_xtr=[None,None]
user_Ti=None
- mapMesh = 0
+ mapMesh = 1
# Time solver parameters
Time_Domain=True # True for unsteady solver, False for steady solver
@@ -89,7 +89,7 @@ def main():
# mesh the geometry
print(ccolors.ANSI_BLUE + 'PyMeshing...' + ccolors.ANSI_RESET)
tms['msh'].start()
- pars = {'xLgt' : 5, 'yLgt' : 5, 'msF' : 1, 'msTe' : 0.005, 'msLe' : 0.0001}
+ pars = {'xLgt' : 5, 'yLgt' : 5, 'msF' : 1, 'msTe' : 0.01, 'msLe' : 0.00001}
msh, gmshWriter = floD.mesh(dim, 'models/n0012.geo', pars, ['field', 'airfoil', 'downstream'])
tms['msh'].stop()
diff --git a/dart/tests/bliTransonic.py b/dart/tests/bliTransonic.py
index bc3fd2e..eee531e 100755
--- a/dart/tests/bliTransonic.py
+++ b/dart/tests/bliTransonic.py
@@ -68,7 +68,7 @@ def main():
user_xtr=[None,None]
user_Ti=None
- mapMesh = 0
+ mapMesh = 1
# Time solver parameters
Time_Domain=True # True for unsteady solver, False for steady solver
diff --git a/dart/viscous/Drivers/VDriver.py b/dart/viscous/Drivers/VDriver.py
index 76c926e..1867c33 100644
--- a/dart/viscous/Drivers/VDriver.py
+++ b/dart/viscous/Drivers/VDriver.py
@@ -136,15 +136,24 @@ class Driver:
dataIn = [None,None]
for n in range(0, len(groups)):
groups[n].connectListNodes, groups[n].connectListElems, dataIn[n] = groups[n].connectList()
- inMesh, paramDict, data, bNodes, inMeshbNodes, xxInit = GroupConstructor().mergeVectors(dataIn, self.mapMesh, nFactor)
+ fMeshDict, cMeshDict, data = GroupConstructor().mergeVectors(dataIn, self.mapMesh, nFactor)
# Initialize data container.
- var = Variables(data, paramDict, bNodes, self.xtrF, self.Ti, self.Re)
+ var = Variables(fMeshDict, self.xtrF, self.Ti, self.Re)
if self.it!=0:
# Extenal flow local markers.
var.extPar[4::var.nExtPar] = data[:,8]
+ """if self.it >= 0:
+ plt.plot(var.xGlobal[:var.bNodes[1]], var.extPar[2:var.bNodes[1]*var.nExtPar: var.nExtPar])
+ plt.plot(cMeshDict['globCoord'][:cMeshDict['bNodes'][1]], cMeshDict['vtExt'][:cMeshDict['bNodes'][1]],'ob')
+ plt.plot(var.xGlobal[var.bNodes[1]:var.bNodes[2]], var.extPar[var.bNodes[1]*var.nExtPar + 2:var.bNodes[2]*var.nExtPar: var.nExtPar])
+ plt.plot(cMeshDict['globCoord'][cMeshDict['bNodes'][1]:cMeshDict['bNodes'][2]], cMeshDict['vtExt'][cMeshDict['bNodes'][1]:cMeshDict['bNodes'][2]],'or')
+ plt.plot(var.xGlobal[var.bNodes[2]:], var.extPar[var.bNodes[2]*var.nExtPar + 2:: var.nExtPar])
+ plt.plot(cMeshDict['globCoord'][cMeshDict['bNodes'][2]:], cMeshDict['vtExt'][cMeshDict['bNodes'][2]:],'og')
+ plt.show()"""
+
# Boundary condition (Airfoil/Wake).
DirichletBC=Conditions.Boundaries(self.Re)
@@ -165,8 +174,8 @@ class Driver:
if self.it == 0:
print('+------------------------------- Solver setup ---------------------------------+')
if self.mapMesh:
- print('| - Mesh mapped from {:>5.0f} to {:>5.0f} elements.{:>35s}'.format(len(inMesh), var.nN, '|'))
- print('| - Number of elements: {:<.0f}. {:>51s}'.format(var.nN,'|'))
+ print('| - Mesh mapped from {:>5.0f} to {:>5.0f} elements.{:>35s}'.format(len(cMeshDict['locCoord']), var.nN, '|'))
+ print('| - Number of elements: {:>5.0f}. {:>49s}'.format(var.nN,'|'))
if var.xtrF[0] is None:
print('| - Free transition on the upper side. {:>41}'.format('|'))
else:
@@ -177,7 +186,7 @@ class Driver:
print('| - Forced transition on the lower side; x/c = {:<.4f}. {:>25s}'.format(var.xtrF[1], '|'))
print('| - Critical amplification ratio: {:<.2f}. {:>40s}'.format(var.Ncrit,'|'))
- # Output solver parameters.
+ # Output solver parameters.
print('|{:>79}'.format('|'))
print('| Numerical parameters {:>57}'.format('|'))
print('| - CFL0: {:<3.2f} {:>65}'.format(tSolver.getCFL0(),'|'))
@@ -200,7 +209,10 @@ class Driver:
initializer=Conditions.Initial(DirichletBC)
# Compute initial condition based on Twaithes solution.
- initializer.initialConditions(var) # Set boundary condition
+ initializer.initialConditions(var) # Set boundary condition
+
+ tSolver.jacEval0 = 1
+ tSolver.setCFL0(1)
print('| - Using Twaithes solution as initial guess. {:>34}'.format('|'))
@@ -239,7 +251,7 @@ class Driver:
self.blParPrevious=var.blPar.copy()
# Compute blowing velocities and sort the boundary layer parameters.
- self.blScalars, AFblwVel, WKblwVel, AFblVectors, WKblVectors, AFdeltaStarOut, WKdeltaStarOut = PostProcess().sortParam(var, self.mapMesh, inMesh, inMeshbNodes, xxInit)
+ self.blScalars, AFblwVel, WKblwVel, AFblVectors, WKblVectors, AFdeltaStarOut, WKdeltaStarOut = PostProcess().sortParam(var, self.mapMesh, cMeshDict)
self.xtr = var.xtr.copy()
self.blVectors = AFblVectors
@@ -249,7 +261,7 @@ class Driver:
groups[0].HEnd = [AFblVectors[2][var.bNodes[1]-1], AFblVectors[2][var.bNodes[2]-1]]
groups[0].CtEnd = [AFblVectors[8][var.bNodes[1]-1], AFblVectors[8][var.bNodes[2]-1]]
groups[0].deltaStar = AFdeltaStarOut
- groups[0].xx = xxInit[:inMeshbNodes[2]]
+ groups[0].xx = cMeshDict['locCoord'][:cMeshDict['bNodes'][2]]
# Sort the following results in reference frame.
groups[0].deltaStar = groups[0].deltaStar[groups[0].connectListNodes.argsort()]
@@ -263,12 +275,13 @@ class Driver:
self.Cdp = self.blScalars[2] # Cdf comes from the airfoil as there is no friction in the wake.
groups[1].deltaStar = WKdeltaStarOut
- groups[1].xx = xxInit[inMeshbNodes[2]:]
+ groups[1].xx = cMeshDict['locCoord'][cMeshDict['bNodes'][2]:]
# Sort the following results in reference frame.
groups[1].deltaStar = groups[1].deltaStar[groups[1].connectListNodes.argsort()]
groups[1].xx = groups[1].xx[groups[1].connectListNodes.argsort()]
groups[1].u = WKblwVel[groups[1].connectListElems.argsort()]
-
- """self.writeFile()
- Validation().main(1,self.wData, WKblVectors, var.xGlobal[var.bNodes[2]:]) """
\ No newline at end of file
+
+ """if self.it >= 0:
+ self.writeFile()
+ Validation().main(1,self.wData, WKblVectors, var.xGlobal[var.bNodes[2]:])"""
\ No newline at end of file
diff --git a/dart/viscous/Drivers/VGroupConstructor.py b/dart/viscous/Drivers/VGroupConstructor.py
index bca6618..9d1ed68 100755
--- a/dart/viscous/Drivers/VGroupConstructor.py
+++ b/dart/viscous/Drivers/VGroupConstructor.py
@@ -31,77 +31,116 @@ class GroupConstructor():
"""Merges 3 groups upper/lower sides and wake.
"""
- # Saves the initial mesh.
- inMesh = np.concatenate((dataIn[0][0][:,0], dataIn[0][1][1:,0], dataIn[1][:,0]))
- inMeshbNodes = [0, len(dataIn[0][0][:,0]), len(dataIn[0][0][:,0]) + len(dataIn[0][1][1:,0])]
-
- xxInit_up, _, _, _, _ = self.__getBLcoordinates(dataIn[0][0])
- xxInit_lw, _, _, _, _ = self.__getBLcoordinates(dataIn[0][1])
- xxInit_lw = np.delete(xxInit_lw, 0) # Delete stagnation point doublon.
- xxInit_wk, _, _, _, _ = self.__getBLcoordinates(dataIn[1])
- xxInit = np.concatenate((xxInit_up, xxInit_lw, xxInit_wk))
-
for k in range(len(dataIn)):
if len(dataIn[k]) == 2: # Airfoil case.
-
- if mapMesh:
- dataUpper = self.interpolateToFineMesh(dataIn[k][0], nFactor)
- dataLower = self.interpolateToFineMesh(dataIn[k][1], nFactor)
- else:
- dataUpper = dataIn[k][0]
- dataLower = dataIn[k][1]
-
- xx_up, dv_up, vtExt_up, _, alpha_up = self.__getBLcoordinates(dataUpper)
- xx_lw, dv_lw, vtExt_lw, _, alpha_lw = self.__getBLcoordinates(dataLower)
+
+ xx_up, dv_up, vtExt_up, _, alpha_up = self.__getBLcoordinates(dataIn[k][0])
+ xx_lw, dv_lw, vtExt_lw, _, alpha_lw = self.__getBLcoordinates(dataIn[k][1])
else: # Wake case
- if mapMesh:
- dataWake = self.interpolateToFineMesh(dataIn[k], nFactor)
- else:
- dataWake = dataIn[k]
-
- xx_wk, dv_wk, vtExt_wk, _, alpha_wk = self.__getBLcoordinates(dataWake)
-
- # Delete stagnation point doublon
- xx_lw = np.delete(xx_lw,0)
- dv_lw = np.delete(dv_lw,0)
- vtExt_lw = np.delete(vtExt_lw,0)
- alpha_lw = np.delete(alpha_lw,0)
-
- # Nodes at the boundary of the computational domain:
- # [Stagnation point, First lower side point, First wake point]
- boundaryNodes = np.zeros(3, dtype=int)
- boundaryNodes[0] = 0
- boundaryNodes[1] = len(xx_up)
- boundaryNodes[2] = len(xx_up) + len(xx_lw)
-
- param={'xx':None, 'dv':None, 'vtExt':None, 'alpha':None}
- param['xx'] = np.concatenate((xx_up,xx_lw,xx_wk))
- param['dv'] = np.concatenate((dv_up,dv_lw,dv_wk))
- param['vtExt'] = np.concatenate((vtExt_up,vtExt_lw,vtExt_wk))
- param['alpha'] = np.concatenate((alpha_up,alpha_lw,alpha_wk))
-
- data = np.zeros((len(param['xx']),9))
- data = np.concatenate((dataUpper, dataLower, dataWake), axis=0)
- data = np.delete(data,boundaryNodes[1], axis = 0) # Delete stagnation point doublon
-
- dx = np.zeros(len(param['xx']))
- for i in range(1, len(param['xx'])):
-
- if i == boundaryNodes[1]:
-
- dx[i] = param['xx'][i] - param['xx'][0]
- elif i == boundaryNodes[2]:
-
- dx[i] = 0
-
- else:
+ xx_wk, dv_wk, vtExt_wk, _, alpha_wk = self.__getBLcoordinates(dataIn[k])
+
+ if mapMesh:
+ # Save initial mesh.
+ cMesh = np.concatenate((dataIn[0][0][:,0], dataIn[0][1][1:,0], dataIn[1][:,0]))
+ cMeshbNodes = [0, len(dataIn[0][0][:,0]), len(dataIn[0][0][:,0]) + len(dataIn[0][1][1:,0])]
+ cxx = np.concatenate((xx_up, xx_lw[1:], xx_wk))
+ cvtExt = np.concatenate((vtExt_up, vtExt_lw[1:], vtExt_wk))
+
+ # Create fine mesh.
+ fMeshUp = self.createFineMesh(dataIn[0][0][:,0], nFactor)
+ fMeshLw = self.createFineMesh(dataIn[0][1][:,0], nFactor)
+ fMeshWk = self.createFineMesh(dataIn[1][:,0] , nFactor)
+ fMesh = np.concatenate((fMeshUp, fMeshLw[1:], fMeshWk))
+ fMeshbNodes = [0, len(fMeshUp), len(fMeshUp) + len(fMeshLw[1:])]
+
+ fxx_up = self.interpolateToFineMesh(xx_up, fMeshUp, nFactor)
+ fxx_lw = self.interpolateToFineMesh(xx_lw, fMeshLw, nFactor)
+ fxx_wk = self.interpolateToFineMesh(xx_wk, fMeshWk, nFactor)
+ fxx = np.concatenate((fxx_up, fxx_lw[1:], fxx_wk))
+
+ fvtExt_up = self.interpolateToFineMesh(vtExt_up, fMeshUp, nFactor)
+ fvtExt_lw = self.interpolateToFineMesh(vtExt_lw, fMeshLw, nFactor)
+ fvtExt_wk = self.interpolateToFineMesh(vtExt_wk, fMeshWk, nFactor)
+ fvtExt = np.concatenate((fvtExt_up, fvtExt_lw[1:], fvtExt_wk))
+
+ fMe_up = self.interpolateToFineMesh(dataIn[0][0][:,5], fMeshUp, nFactor)
+ fMe_lw = self.interpolateToFineMesh(dataIn[0][1][:,5], fMeshLw, nFactor)
+ fMe_wk = self.interpolateToFineMesh(dataIn[1][:,5] , fMeshWk, nFactor)
+ fMe = np.concatenate((fMe_up, fMe_lw[1:], fMe_wk))
+
+ frho_up = self.interpolateToFineMesh(dataIn[0][0][:,6], fMeshUp, nFactor)
+ frho_lw = self.interpolateToFineMesh(dataIn[0][1][:,6], fMeshLw, nFactor)
+ frho_wk = self.interpolateToFineMesh(dataIn[1][:,6] , fMeshWk, nFactor)
+ frho = np.concatenate((frho_up, frho_lw[1:], frho_wk))
+
+ fdStarExt_up = self.interpolateToFineMesh(dataIn[0][0][:,7], fMeshUp, nFactor)
+ fdStarExt_lw = self.interpolateToFineMesh(dataIn[0][1][:,7], fMeshLw, nFactor)
+ fdStarExt_wk = self.interpolateToFineMesh(dataIn[1][:,7] , fMeshWk, nFactor)
+ fdStarext = np.concatenate((fdStarExt_up, fdStarExt_lw[1:], fdStarExt_wk))
+
+ fdv = np.concatenate((dv_up, dv_lw[1:], dv_wk))
+
+ fMeshDict = {'globCoord': fMesh, 'bNodes': fMeshbNodes, 'locCoord': fxx,
+ 'vtExt': fvtExt, 'Me': fMe, 'rho': frho, 'deltaStarExt': fdStarext, 'dv': fdv}
+
+ cMe = np.concatenate((dataIn[0][0][:,5], dataIn[0][1][1:,5], dataIn[1][:,5]))
+ cRho = np.concatenate((dataIn[0][0][:,6], dataIn[0][1][1:,6], dataIn[1][:,6]))
+ cdStarExt = np.concatenate((dataIn[0][0][:,7], dataIn[0][1][1:,7], dataIn[1][:,7]))
+ cMeshDict = {'globCoord': cMesh, 'bNodes': cMeshbNodes, 'locCoord': cxx,
+ 'vtExt': cvtExt, 'Me': cMe, 'rho': cRho, 'deltaStarExt': cdStarExt, 'dv': fdv}
+
+ dataUpper = np.empty((len(fMeshUp), 0))
+ dataLower = np.empty((len(fMeshLw), 0))
+ dataWake = np.empty((len(fMeshWk), 0))
+ for iParam in range(len(dataIn[0][0][0,:])):
+ interpParamUp = self.interpolateToFineMesh(dataIn[0][0][:,iParam], fMeshUp, nFactor)
+ dataUpper = np.column_stack((dataUpper, interpParamUp))
+ interpParamLw = self.interpolateToFineMesh(dataIn[0][1][:,iParam], fMeshLw, nFactor)
+ dataLower = np.column_stack((dataLower, interpParamLw))
+ interpParamWk = self.interpolateToFineMesh(dataIn[1][:,iParam] , fMeshWk, nFactor)
+ dataWake = np.column_stack((dataWake, interpParamWk))
+
+ dataLower = np.delete(dataLower, 0, axis = 0) # Remove stagnation point doublon.
+
+ else:
+
+ Mesh = np.concatenate((dataIn[0][0][:,0], dataIn[0][1][1:,0], dataIn[1][:,0]))
+ MeshbNodes = [0, len(dataIn[0][0][:,0]), len(dataIn[0][0][:,0]) + len(dataIn[0][1][1:,0])]
+
+ xx = np.concatenate((xx_up, xx_lw[1:], xx_wk))
+ vtExt = np.concatenate((vtExt_up, vtExt_lw[1:], vtExt_wk))
+
+ alpha = np.concatenate((alpha_up, alpha_lw, alpha_wk))
+ dv = np.concatenate((dv_up, dv_lw[1:], dv_wk))
+
+ Me = np.concatenate((dataIn[0][0][:,5], dataIn[0][1][1:,5], dataIn[1][:,5]))
+ rho = np.concatenate((dataIn[0][0][:,6], dataIn[0][1][1:,6], dataIn[1][:,6]))
+ dStarext = np.concatenate((dataIn[0][0][:,7], dataIn[0][1][1:,7], dataIn[1][:,7]))
+
+ fMeshDict = {'globCoord': Mesh, 'bNodes': MeshbNodes, 'locCoord': xx,
+ 'vtExt': vtExt, 'Me': Me, 'rho': rho, 'deltaStarExt': dStarext, 'dv': dv}
+
+ cMeshDict = fMeshDict.copy()
+
+ dataUpper = dataIn[0][0]
+ dataLower = dataIn[0][1][1:, :]
+ dataWake = dataIn[1]
+
+ data = np.concatenate((dataUpper, dataLower, dataWake), axis = 0)
- dx[i] = param['xx'][i] - param['xx'][i-1]
+ """param = 'rho'
+ plt.plot(fMeshDict['globCoord'][:fMeshDict['bNodes'][1]], fMeshDict[param][:fMeshDict['bNodes'][1]],'-b')
+ plt.plot(cMeshDict['globCoord'][:cMeshDict['bNodes'][1]], cMeshDict[param][:cMeshDict['bNodes'][1]],'ob')
+ plt.plot(fMeshDict['globCoord'][fMeshDict['bNodes'][1]:fMeshDict['bNodes'][2]], fMeshDict[param][fMeshDict['bNodes'][1]:fMeshDict['bNodes'][2]],'-r')
+ plt.plot(cMeshDict['globCoord'][cMeshDict['bNodes'][1]:cMeshDict['bNodes'][2]], cMeshDict[param][cMeshDict['bNodes'][1]:cMeshDict['bNodes'][2]],'or')
+ plt.plot(fMeshDict['globCoord'][fMeshDict['bNodes'][2]:], fMeshDict[param][fMeshDict['bNodes'][2]:],'-g')
+ plt.plot(cMeshDict['globCoord'][cMeshDict['bNodes'][2]:], cMeshDict[param][cMeshDict['bNodes'][2]:],'og')
+ plt.show()"""
- return inMesh, param, data, boundaryNodes, inMeshbNodes, xxInit
+ return fMeshDict, cMeshDict, data
def __getBLcoordinates(self,data):
"""Transform the reference coordinates into airfoil coordinates
@@ -126,33 +165,30 @@ class GroupConstructor():
return xx, dv, vt, nN, alpha
- def interpolateToFineMesh(self, iData, nFactor):
-
- inMesh = iData[:,0]
-
- outMesh = []
- for iPoint in range(len(inMesh) - 1):
-
- dx = inMesh[iPoint + 1] - inMesh[iPoint]
-
- for iRefinedPoint in range(nFactor):
+ def interpolateToFineMesh(self, cVector, fMesh, nFactor):
+
+ fVector = np.empty(len(fMesh)-1)
+ for cPoint in range(len(cVector) - 1):
+ dv = cVector[cPoint + 1] - cVector[cPoint]
- outMesh = np.append(outMesh, inMesh[iPoint] + iRefinedPoint * dx / nFactor)
+ for fPoint in range(nFactor):
+ fVector[nFactor * cPoint + fPoint] = cVector[cPoint] + fPoint * dv / nFactor
+
+ fVector = np.append(fVector, cVector[-1])
- outMesh = np.append(outMesh, inMesh[-1])
- outData = np.empty((len(outMesh), len(iData[0,:])))
- outData[:,0] = outMesh
+ return fVector
- for iParam in range(1, len(iData[0, :])):
+ def createFineMesh(self, cMesh, nFactor):
- # Create interpolation function between the initial mesh
- # and the function to evaluate.
+ fMesh = []
+ for iPoint in range(len(cMesh) - 1):
- f_interp = interpolate.interp1d(inMesh, iData[:,iParam])
+ dx = cMesh[iPoint + 1] - cMesh[iPoint]
- # Map the function to the new mesh.
- outData[:, iParam] = f_interp(outData[:, 0])
+ for iRefinedPoint in range(nFactor):
- return outData
+ fMesh = np.append(fMesh, cMesh[iPoint] + iRefinedPoint * dx / nFactor)
+ fMesh = np.append(fMesh, cMesh[-1])
+ return fMesh
diff --git a/dart/viscous/Drivers/VPostProcess.py b/dart/viscous/Drivers/VPostProcess.py
index 3d844db..92f0d7b 100644
--- a/dart/viscous/Drivers/VPostProcess.py
+++ b/dart/viscous/Drivers/VPostProcess.py
@@ -6,7 +6,7 @@ class PostProcess:
def __init__(self):
pass
- def sortParam(self,var, mapMesh, inMesh, inMeshbNodes, xxInit):
+ def sortParam(self,var, mapMesh, cMeshDict):
# Recompute deltaStar.
@@ -14,6 +14,9 @@ class PostProcess:
var.blPar[9::var.nBlPar] = var.u[0::var.nVar] * var.u[1::var.nVar]
if mapMesh:
+ inMesh = cMeshDict['globCoord']
+ inMeshbNodes = cMeshDict['bNodes']
+ xxInit = cMeshDict['locCoord']
rho, invVel, deltaStar = self.inverseMeshMapping(var, inMesh, inMeshbNodes)
@@ -52,8 +55,6 @@ class PostProcess:
blwVel[iPoint-1] = (var.extPar[iPoint*var.nExtPar+1] * var.u[iPoint*var.nVar+3] * var.blPar[iPoint*var.nBlPar+9]
- var.extPar[iPrev*var.nExtPar+1] * var.u[iPrev*var.nVar+3] * var.blPar[iPrev*var.nBlPar+9]) / (var.extPar[iPoint*var.nExtPar+1] * (var.xx[iPoint] - var.xx[iPrev]))
- print('blwVel bnodes2', blwVel[var.bNodes[2]])
- print('blwVel bnodes2-1', blwVel[var.bNodes[2]-1])
blwVel = np.delete(blwVel, var.bNodes[2])
# Separate airfoil and wake blowing velocities.
AFblwVel=blwVel[:var.bNodes[2]-1]
diff --git a/dart/viscous/Physics/VClosures.py b/dart/viscous/Physics/VClosures.py
index 28e4491..1acad30 100755
--- a/dart/viscous/Physics/VClosures.py
+++ b/dart/viscous/Physics/VClosures.py
@@ -28,10 +28,9 @@ class Closures:
Hk = (H - 0.29*Me**2)/(1+0.113*Me**2) # Kinematic shape factor
if name == "airfoil":
Hk = max(Hk, 1.02)
- Hk = min(Hk,6)
else:
Hk = max(Hk, 1.00005)
- Hk = min(Hk,6)
+ Hk = min(Hk,6)
Hstar2 = (0.064/(Hk-0.8)+0.251)*Me**2 # Density shape parameter
delta = min((3.15+ H + (1.72/(Hk-1)))*theta, 12*theta)
Hks = Hk - 4.35
@@ -63,40 +62,35 @@ class Closures:
"""Turbulent closure derived from Nishida and Drela (1996)
"""
# eliminate absurd transients
- Ct = min(Ct, 0.30)
- Ct = max(Ct, 0.0000001)
+ Ct = max(min(Ct, 0.30), 0.0000001)
Hk = (H - 0.29*Me**2)/(1+0.113*Me**2)
- if name == 'wake':
- Hk = max(Hk, 1.00005)
- Hk = min(Hk,10)
- else:
- Hk = max(Hk, 1.00005)
- Hk = min(Hk,10)
+ Hk = max(min(Hk, 10), 1.00005)
Hstar2 = ((0.064/(Hk-0.8))+0.251)*Me**2
gam = 1.4 - 1
Fc = np.sqrt(1+0.5*gam*Me**2)
+
if Ret <= 400:
H0 = 4
elif Ret > 400:
H0 = 3 + (400/Ret)
- if Ret > 200:
- Ret = Ret
- elif Ret <= 200:
- Ret = 200
+
+ Ret = max(Ret, 200)
+
if Hk <= H0:
Hstar = ((0.5-4/Ret)*((H0-Hk)/(H0-1))**2)*(1.5/(Hk+0.5))+1.5+4/Ret
elif Hk > H0:
Hstar = (Hk-H0)**2*(0.007*np.log(Ret)/(Hk-H0+4/np.log(Ret))**2 + 0.015/Hk) + 1.5 + 4/Ret
+
Hstar = (Hstar + 0.028*Me**2)/(1+0.014*Me**2) # Whitfield's minor additional compressibility correction
logRt = np.log(Ret/Fc)
- logRt = np.max((logRt,3))
- arg = np.max((-1.33*Hk, -20))
+ logRt = max(logRt,3)
+ arg = max(-1.33*Hk, -20)
Us = (Hstar/2)*(1-4*(Hk-1)/(3*H)) # Equivalent normalized wall slip velocity
delta = min((3.15+ H + (1.72/(Hk-1)))*theta, 12*theta)
Ctcon = 0.5/(6.7**2*0.75)
+
if name == 'wake':
- if Us > 0.99995:
- Us = 0.99995
+ Us = min(Us, 0.99995)
n = 2 # two wake halves
Cf = 0 # no friction inside the wake
Hkc = Hk - 1
@@ -109,15 +103,16 @@ class Closures:
n = 1
Cf = (1/Fc)*(0.3*np.exp(arg)*(logRt/2.3026)**(-1.74-0.31*Hk)+0.00011*(np.tanh(4-(Hk/0.875))-1))
Hkc = max(Hk-1-18/Ret, 0.01)
- # dissipation coefficient
+ # Dissipation coefficient
Hmin = 1 + 2.1/np.log(Ret)
Fl = (Hk-1)/(Hmin-1)
Dfac = 0.5+0.5*np.tanh(Fl)
Cdw = 0.5*(Cf*Us) * Dfac
Cdd = (0.995-Us)*Ct**2
Cdl = 0.15*(0.995-Us)**2/Ret
- if n*Hstar*Ctcon*((Hk-1)*Hkc**2)/((1-Us)*(Hk**2)*H) >= 0:
- Cteq = np.sqrt(n*Hstar*Ctcon*((Hk-1)*Hkc**2)/((1-Us)*(Hk**2)*H)) #Here it is Cteq^0.5
+ CteqSquared = n*Hstar*Ctcon*((Hk-1)*Hkc**2)/((1-Us)*(Hk**2)*H)
+ if CteqSquared >= 0:
+ Cteq = np.sqrt(CteqSquared) #Here it is Cteq^0.5
else:
Cteq = 0.03
Cd = n*(Cdw+ Cdd + Cdl)
@@ -206,7 +201,7 @@ class Closures:
Rnorm = (np.log10(Ret) - (logRet_crit-Dgr))/(2*Dgr)
if Rnorm <=0:
Rfac = 0
- if Rnorm >= 1:
+ elif Rnorm >= 1:
Rfac = 1
else:
Rfac = 3*Rnorm**2 - 2*Rnorm**3
@@ -215,7 +210,7 @@ class Closures:
arg = 3.87*Hmi-2.52
dndRet = 0.028*(Hk-1)-0.0345*np.exp(-(arg**2))
Ax = (m*dndRet/theta)*Rfac
- # set correction for separated profiles
+ # Set correction for separated profiles
if Hk > Hk1:
Hnorm = (Hk - Hk1)/(Hk2-Hk1)
if Hnorm >=1:
@@ -229,6 +224,5 @@ class Closures:
if Ax2 < 0:
Ax2 = 0
Ax = Hfac*Ax2 + (1-Hfac)*Ax1
- if Ax < 0:
- Ax = 0
+ Ax = max(Ax, 0)
return Ax
\ No newline at end of file
diff --git a/dart/viscous/Physics/VDataStructures.py b/dart/viscous/Physics/VDataStructures.py
index 7e9903e..ca867fb 100755
--- a/dart/viscous/Physics/VDataStructures.py
+++ b/dart/viscous/Physics/VDataStructures.py
@@ -16,30 +16,28 @@
# ------------------------------------------------------------------------------
import numpy as np
-class DGeometry:
+class PGeometry:
def __init__(self, data, paramDict, _bNodes):
# Private
- self.__globCoord = data[:,0] # Global coordinates of each marker
- self.__locCoord = paramDict['xx'] # Local coordinates of each marker
- self.__nMarkers = len(self.__locCoord) # Number of markers in the computational domain
+ self.__globCoord = data[:,0] # Global coordinates of each marker
+ self.__locCoord = paramDict['xx'] # Local coordinates of each marker
+ self.__nMarkers = len(self.__locCoord) # Number of markers in the computational domain
self.__boundaryNodes = _bNodes # Boundary nodes of the computational domain [Stagnation point,
# First point on lower side,
# First wake point]
- # Public
-
- def GetglobCoord(self): return self.__globCoord
- def GetlocCoord(self): return self.__locCoord
- def GetnMarkers(self): return self.__nMarkers
- def GetglobCoord(self): return self.__bNodes
- def GetboundaryNodes(self): return self.__boundaryNodes
+ def GetglobCoord(self): return self.__globCoord
+ def GetlocCoord(self): return self.__locCoord
+ def GetnMarkers(self): return self.__nMarkers
+ def GetglobCoord(self): return self.__bNodes
+ def GetboundaryNodes(self): return self.__boundaryNodes
-class DSolution:
+class PSolution:
def __init__(self, _xtrF):
@@ -79,7 +77,7 @@ class DSolution:
self.__xtr[side] = xtrIn
pass
- def initializeSolution(self, DGeometry):
+ def initializeSolution(self, PGeometry):
nMarkers = PGeometry.GetnMarkers()
self.u = np.empty(self.__nVar * nMarkers)
@@ -87,7 +85,7 @@ class DSolution:
return self.u, self.blPar
-class DOuterFlow:
+class POuterFlow:
def __init__(self, data) -> None:
@@ -100,6 +98,14 @@ class DOuterFlow:
self.outerLocCoord =
self.outerRe =
+class PDataContainer:
+ def __init__(self, _Re) -> None:
+
+ self.Re = _Re # Flow reynolds number.
+ self.Geometry = PGeometry() # Geometry components.
+ self.Solution = PSolution() # Solution components.
+ self.OuterFlow = POuterFlow() # Outer (inviscid) flow parameters.
+
@@ -214,4 +220,4 @@ class DVariables:
Ncrit=-8.43-2.4*np.log(tau/100) # Drela 1998
else:
Ncrit=9
- return Ncrit
+ return Ncrit
\ No newline at end of file
diff --git a/dart/viscous/Physics/VVariables.py b/dart/viscous/Physics/VVariables.py
index 094e74d..5e4e026 100755
--- a/dart/viscous/Physics/VVariables.py
+++ b/dart/viscous/Physics/VVariables.py
@@ -66,13 +66,13 @@ class Variables:
- Critical amplification ratio for the e^N method (default value= 9.)
"""
- def __init__(self, data, paramDict, bNodes, xtrF, Ti, Re):
+ def __init__(self, fMeshDict, xtrF, Ti, Re):
self.__Ti = Ti
self.Re = Re
- self.xx = paramDict['xx']
+ self.xx = fMeshDict['locCoord']
self.nN = len(self.xx)
- self.bNodes = bNodes
+ self.bNodes = fMeshDict['bNodes']
# Solution
self.nVar = 5
self.u = np.zeros(self.nVar*self.nN)
@@ -93,13 +93,13 @@ class Variables:
self.nExtPar = 6
self.extPar = np.zeros(self.nExtPar*self.nN)
- self.extPar[0::self.nExtPar] = data[:,5]
- self.extPar[1::self.nExtPar] = data[:,6]
- self.extPar[2::self.nExtPar] = paramDict['vtExt']
- self.extPar[3::self.nExtPar] = data[:,7]
- self.extPar[4::self.nExtPar] = paramDict['xx']
- self.dv = paramDict['dv']
- self.xGlobal = data[:,0]
+ self.extPar[0::self.nExtPar] = fMeshDict['Me']
+ self.extPar[1::self.nExtPar] = fMeshDict['rho']
+ self.extPar[2::self.nExtPar] = fMeshDict['vtExt']
+ self.extPar[3::self.nExtPar] = fMeshDict['deltaStarExt']
+ self.extPar[4::self.nExtPar] = fMeshDict['locCoord']
+ self.dv = fMeshDict['dv']
+ self.xGlobal = fMeshDict['globCoord']
self.flowRegime = np.zeros(self.nN, dtype=int)
@@ -122,11 +122,13 @@ class Variables:
pass
def turbIntensity(self,Ti):
- '''Computes the critical amplification ratio Ncrit based on the input
- free-stream turbulence intensity. Ncrit=9 by default.'''
- if Ti is not None and Ti<=1:
- tau=2.7*np.tanh(self.Ti/2.7)
- Ncrit=-8.43-2.4*np.log(tau/100) # Drela 1998
+ """Computes the critical amplification ratio Ncrit based on the input
+ free-stream turbulence intensity. Ncrit=9 by default.
+ """
+
+ if Ti is not None and Ti <= 1:
+ tau = 2.7 * np.tanh(self.Ti / 2.7)
+ Ncrit = -8.43 - 2.4*np.log(tau / 100) # Drela 1998
else:
- Ncrit=9
+ Ncrit = 9
return Ncrit
diff --git a/dart/viscous/Solvers/VSolver.py b/dart/viscous/Solvers/VSolver.py
index 91bf0b9..64a8e58 100644
--- a/dart/viscous/Solvers/VSolver.py
+++ b/dart/viscous/Solvers/VSolver.py
@@ -154,8 +154,6 @@ class flowEquations:
# Shear lag equation ignored in the laminar portion of the flow
timeMatrix[4] = [0, 0, 0, 0, 1]
- if np.linalg.det(timeMatrix) == 0:
- raise RuntimeError('Singular time matrix.')
sysRes = np.linalg.solve(timeMatrix, spaceMatrix).flatten()
return sysRes
diff --git a/dart/viscous/Solvers/VTimeSolver.py b/dart/viscous/Solvers/VTimeSolver.py
index 1b4fe23..3e6de63 100644
--- a/dart/viscous/Solvers/VTimeSolver.py
+++ b/dart/viscous/Solvers/VTimeSolver.py
@@ -18,6 +18,7 @@ from dart.viscous.Physics.VClosures import Closures
from matplotlib import pyplot as plt
from fwk.coloring import ccolors
+import scipy.sparse.linalg as linSolver
import numpy as np
import math
@@ -122,8 +123,10 @@ class implicitEuler(timeConfig):
# Initialize time parameters
self.__CFL0=_cfl0
- self.__maxIt = 20000
+ self.__maxIt = 15000
self.__tol = 1e-8
+ self.jacEval0 = 5
+
self.solver=_solver
self.transSolver = transition
self.wakeBC = wakeBC
@@ -324,7 +327,7 @@ class implicitEuler(timeConfig):
CFL = self.__CFL0 # Initialized CFL number.
CFLAdapt = 1 # Flag for CFL adaptation.
- jacEval = 3 # Number of iterations for which Jacobian is frozen.
+ jacEval = self.jacEval0 # Number of iterations for which Jacobian is frozen.
dt = self.computeTimeStep(CFL, dx, iInvVel) # Initial time step.
IMatrix=np.identity(nVar) / dt # Identity matrix used to construct the linear system.
@@ -344,6 +347,7 @@ class implicitEuler(timeConfig):
# Linear solver.
slnIncr = np.linalg.solve((Jacobian + IMatrix), - sysRes)
+ #slnIncr, exitCode = linSolver.gmres((Jacobian + IMatrix), -sysRes, tol = 1e-10)
# Increment solution and correct absurd transcient results.
@@ -408,7 +412,7 @@ class implicitEuler(timeConfig):
CFL = self.__CFL0* (normSysRes0/normSysRes)**0.7
- CFL = max(CFL, 0.1)
+ CFL = max(CFL, 0.01)
dt = self.computeTimeStep(CFL, dx, iInvVel)
IMatrix=np.identity(nVar) / dt
@@ -418,8 +422,8 @@ class implicitEuler(timeConfig):
else:
# Maximum number of iterations reached.
if normSysRes >= 1e-3 * normSysRes0:
- print('|', ccolors.ANSI_RED + 'Too many iteration at position {:<.4f}. normLinSysRes = {:<.3f}. {:>30s}'
- .format(var.xGlobal[iMarker],
+ print('|', ccolors.ANSI_RED + 'Too many iteration at position {:<.4f}. Marker: {:>5.0f}. normLinSysRes = {:<.3f}. {:>30s}'
+ .format(var.xGlobal[iMarker], iMarker,
np.log10(normSysRes/normSysRes0),'|'),
ccolors.ANSI_RESET)
var.u[iMarker*nVar:(iMarker+1)*nVar] = u0[(iMarker)*nVar: (iMarker+1)*nVar].copy()
@@ -440,8 +444,8 @@ class implicitEuler(timeConfig):
var.extPar[iMarker*nExtPar+0], name)
else:
- print('|', ccolors.ANSI_YELLOW+ 'Too many iteration at position {:<.4f}. normLinSysRes = {:<.3f}. {:>30s}'
- .format(var.xGlobal[iMarker],
+ print('|', ccolors.ANSI_YELLOW+ 'Too many iteration at position {:<.4f}. Marker: {:>5.0f} normLinSysRes = {:<.3f}. {:>30s}'
+ .format(var.xGlobal[iMarker], iMarker,
np.log10(normSysRes/normSysRes0),'|'),
ccolors.ANSI_RESET)
--
GitLab