From 57612ad455a37d6d4662d25073838f2ab0161e53 Mon Sep 17 00:00:00 2001
From: Dechamps Paul <paul.dechamps@student.uliege.be>
Date: Thu, 25 Nov 2021 23:00:11 +0100
Subject: [PATCH] Mesh adaptation fix2. Changed initial condition.
---
dart/tests/bliSeparated.py | 6 +-
dart/tests/bliTransonic.py | 2 +-
dart/viscous/Drivers/VDriver.py | 18 +-
dart/viscous/Drivers/VGroupConstructor.py | 19 +-
dart/viscous/Drivers/VPostProcess.py | 343 +++++----
dart/viscous/Physics/VBIConditions.py | 2 +-
dart/viscous/Solvers/VSolver.py | 7 +
dart/viscous/Solvers/VTimeSolver.py | 891 ++++++++++++----------
dart/viscous/Solvers/VTransition.py | 4 +-
9 files changed, 710 insertions(+), 582 deletions(-)
diff --git a/dart/tests/bliSeparated.py b/dart/tests/bliSeparated.py
index 1bc07a1..50181f5 100755
--- a/dart/tests/bliSeparated.py
+++ b/dart/tests/bliSeparated.py
@@ -61,12 +61,12 @@ def main():
# define flow variables
Re = 1e6
alpha = 15.*math.pi/180
- M_inf = 0.
+ M_inf = 0.2
#user_xtr=[0.012,1]
user_xtr=[None,None]
user_Ti=None
- mapMesh = 1
+ mapMesh = 0
# 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.00001}
+ pars = {'xLgt' : 5, 'yLgt' : 5, 'msF' : 1, 'msTe' : 0.01, 'msLe' : 0.0005}
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 eee531e..21a5ede 100755
--- a/dart/tests/bliTransonic.py
+++ b/dart/tests/bliTransonic.py
@@ -90,7 +90,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.01}
+ pars = {'xLgt' : 5, 'yLgt' : 5, 'msF' : 1, 'msTe' : 0.01, 'msLe' : 0.005}
#pars = {'xLgt' : 5, 'yLgt' : 5, 'msF' : 1, 'msTe' : 0.001, 'msLe' : 0.0005}
msh, gmshWriter = floD.mesh(dim, 'models/rae2822.geo', pars, ['field', 'airfoil', 'downstream'])
tms['msh'].stop()
diff --git a/dart/viscous/Drivers/VDriver.py b/dart/viscous/Drivers/VDriver.py
index 1867c33..1e9a017 100644
--- a/dart/viscous/Drivers/VDriver.py
+++ b/dart/viscous/Drivers/VDriver.py
@@ -139,11 +139,11 @@ class Driver:
fMeshDict, cMeshDict, data = GroupConstructor().mergeVectors(dataIn, self.mapMesh, nFactor)
# Initialize data container.
- var = Variables(fMeshDict, 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]
+ var.extPar[4::var.nExtPar] = fMeshDict['xxExt']
"""if self.it >= 0:
plt.plot(var.xGlobal[:var.bNodes[1]], var.extPar[2:var.bNodes[1]*var.nExtPar: var.nExtPar])
@@ -202,7 +202,7 @@ class Driver:
# Initial condition: 'Driver' stores the solution form the previous coupling iteration.
- if self.uPrevious is None or self.it < 1:
+ if self.uPrevious is None or self.it < 5:
# Typically for the first iteration.
# Initalize initial condition builder.
@@ -213,8 +213,10 @@ class Driver:
tSolver.jacEval0 = 1
tSolver.setCFL0(1)
+ tSolver.initSln = 1
- print('| - Using Twaithes solution as initial guess. {:>34}'.format('|'))
+ #print('| - Using Twaithes solution as initial guess. {:>34}'.format('|'))
+ print('| - Time solver will provide the initial solution. {:>29}'.format('|'))
else:
# If we use a solution previously obtained.
@@ -244,17 +246,19 @@ class Driver:
# Reset Ct in the laminar portion of the flow.
# (This quantity is not defined for a laminar flow)
- var.u[4::var.nVar][var.flowRegime == 0] = 0.001
+ #var.u[4::var.nVar][var.flowRegime == 0] = 0.001
# Copy solution.
self.uPrevious=var.u.copy()
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, cMeshDict)
+ self.blScalars, blwVelUp, blwVelLw, blwVelWk, AFblVectors, WKblVectors, AFdeltaStarOut, WKdeltaStarOut = PostProcess().sortParam(var, self.mapMesh, cMeshDict, nFactor)
self.xtr = var.xtr.copy()
self.blVectors = AFblVectors
+ AFblwVel = np.concatenate((blwVelUp, blwVelLw))
+
# Group modification to ensure communication between the solvers.
groups[0].TEnd = [AFblVectors[1][var.bNodes[1]-1], AFblVectors[1][var.bNodes[2]-1]]
@@ -280,7 +284,7 @@ class Driver:
# 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()]
+ groups[1].u = blwVelWk[groups[1].connectListElems.argsort()]
"""if self.it >= 0:
self.writeFile()
diff --git a/dart/viscous/Drivers/VGroupConstructor.py b/dart/viscous/Drivers/VGroupConstructor.py
index 9d1ed68..18f96c2 100755
--- a/dart/viscous/Drivers/VGroupConstructor.py
+++ b/dart/viscous/Drivers/VGroupConstructor.py
@@ -30,7 +30,7 @@ class GroupConstructor():
def mergeVectors(self,dataIn, mapMesh, nFactor):
"""Merges 3 groups upper/lower sides and wake.
"""
-
+
for k in range(len(dataIn)):
if len(dataIn[k]) == 2: # Airfoil case.
@@ -48,6 +48,7 @@ class GroupConstructor():
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))
+ cxxExt = np.concatenate((dataIn[0][0][:,8], dataIn[0][1][1:,8],dataIn[1][:,8]))
# Create fine mesh.
fMeshUp = self.createFineMesh(dataIn[0][0][:,0], nFactor)
@@ -80,17 +81,22 @@ class GroupConstructor():
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))
+
+ fxxExt_up = self.interpolateToFineMesh(dataIn[0][0][:,8], fMeshUp, nFactor)
+ fxxExt_lw = self.interpolateToFineMesh(dataIn[0][1][:,8], fMeshLw, nFactor)
+ fxxExt_wk = self.interpolateToFineMesh(dataIn[1][:,8] , fMeshWk, nFactor)
+ fxxext = np.concatenate((fxxExt_up, fxxExt_lw[1:], fxxExt_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}
+ 'vtExt': fvtExt, 'Me': fMe, 'rho': frho, 'deltaStarExt': fdStarext, 'xxExt': fxxext, '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}
+ 'vtExt': cvtExt, 'Me': cMe, 'rho': cRho, 'deltaStarExt': cdStarExt, 'xxExt': cxxExt, 'dv': fdv}
dataUpper = np.empty((len(fMeshUp), 0))
dataLower = np.empty((len(fMeshLw), 0))
@@ -119,9 +125,10 @@ class GroupConstructor():
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]))
+ xxExt = np.concatenate((dataIn[0][0][:,8], dataIn[0][1][1:,8], dataIn[1][:,8]))
fMeshDict = {'globCoord': Mesh, 'bNodes': MeshbNodes, 'locCoord': xx,
- 'vtExt': vtExt, 'Me': Me, 'rho': rho, 'deltaStarExt': dStarext, 'dv': dv}
+ 'vtExt': vtExt, 'Me': Me, 'rho': rho, 'deltaStarExt': dStarext, 'xxExt': xxExt, 'dv': dv}
cMeshDict = fMeshDict.copy()
@@ -131,8 +138,8 @@ class GroupConstructor():
data = np.concatenate((dataUpper, dataLower, dataWake), axis = 0)
- """param = 'rho'
- plt.plot(fMeshDict['globCoord'][:fMeshDict['bNodes'][1]], fMeshDict[param][:fMeshDict['bNodes'][1]],'-b')
+ param = 'vtExt'
+ """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')
diff --git a/dart/viscous/Drivers/VPostProcess.py b/dart/viscous/Drivers/VPostProcess.py
index 92f0d7b..d4ae02a 100644
--- a/dart/viscous/Drivers/VPostProcess.py
+++ b/dart/viscous/Drivers/VPostProcess.py
@@ -1,143 +1,208 @@
-from matplotlib import pyplot as plt
import numpy as np
-from scipy import interpolate
class PostProcess:
- def __init__(self):
- pass
+ def __init__(self):
+ pass
- def sortParam(self,var, mapMesh, cMeshDict):
-
-
- # Recompute deltaStar.
-
- 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)
-
- AFdeltaStarOut = deltaStar[:inMeshbNodes[2]]
- WKdeltaStarOut = deltaStar[inMeshbNodes[2]:]
-
- # Compute blowing velocity on each cell.
- blwVel=np.zeros(len(xxInit)-1)
-
- for iPoint in range(1, len(xxInit)):
-
- iPrev=0 if iPoint == inMeshbNodes[1] else iPoint - 1
-
- # blwVel[i-1] = (rhoe[i]*vt[i]*deltaStar[i]-rhoe[i-1]*vt[i-1]*deltaStar[i-1])/(rhoe[i]*(xx[i]-xx[i-1]))
- blwVel[iPoint-1] = (rho[iPoint] * invVel[iPoint] * deltaStar[iPoint]
- - rho[iPrev] * invVel[iPrev] * deltaStar[iPrev]) / (rho[iPoint] * (xxInit[iPoint] - xxInit[iPrev]))
-
- blwVel = np.delete(blwVel, inMeshbNodes[2])
- # Separate airfoil and wake blowing velocities.
- AFblwVel=blwVel[:inMeshbNodes[2]-1]
- WKblwVel=blwVel[inMeshbNodes[2]-1:]
-
- else:
-
- AFdeltaStarOut = var.blPar[9:var.bNodes[2]*var.nBlPar:var.nBlPar]
- WKdeltaStarOut = var.blPar[var.bNodes[2]*var.nBlPar + 9::var.nBlPar]
-
- # Compute blowing velocity on each cell.
- blwVel=np.zeros(var.nN-1)
-
- for iPoint in range(1, var.nN):
-
- iPrev=0 if iPoint == var.bNodes[1] else iPoint - 1
-
- # blwVel[i-1] = (rhoe[i]*vt[i]*deltaStar[i]-rhoe[i-1]*vt[i-1]*deltaStar[i-1])/(rhoe[i]*(xx[i]-xx[i-1]))
- 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]))
-
- blwVel = np.delete(blwVel, var.bNodes[2])
- # Separate airfoil and wake blowing velocities.
- AFblwVel=blwVel[:var.bNodes[2]-1]
- WKblwVel=blwVel[var.bNodes[2]-1:]
-
- # Normalize Friction and dissipation coefficients.
- frictionCoeff=var.u[3::var.nVar]**2 * var.blPar[2::var.nBlPar]
- dissipCoeff=var.u[3::var.nVar]**2 * var.blPar[1::var.nBlPar]
-
- # Compute total drag coefficient.
- CdTot=2*var.u[-5]*(var.u[-2]**((var.u[-4]+5)/2))
-
- # Compute friction and pressure drag coefficients.
- Cdf_upper=np.trapz(frictionCoeff[:var.bNodes[1]],var.xGlobal[:var.bNodes[1]])
- # Insert stagnation point in the lower side.
- Cdf_lower=np.trapz(np.insert(frictionCoeff[var.bNodes[1]:var.bNodes[2]],0,frictionCoeff[0]),
- np.insert(var.xGlobal[var.bNodes[1]:var.bNodes[2]],0,var.xGlobal[0]))
- Cdf=Cdf_upper+Cdf_lower # No friction in the wake
- Cdp=CdTot-Cdf
-
- # [CdTot, Cdf, Cdp]
- blScalars=[CdTot, Cdf, Cdp]
-
- # Store boundary layer final parameters in lists.
- # [deltaStar, theta, H, Hk, Hstar, Hstar2, Cf, Cd, Ct, delta]
- blVectors_airfoil=[var.blPar[9:var.bNodes[2]*var.nBlPar:var.nBlPar],
- var.u[0:var.bNodes[2]*var.nVar:var.nVar],
- var.u[1:var.bNodes[2]*var.nVar:var.nVar],
- var.blPar[6:var.bNodes[2]*var.nBlPar:var.nBlPar],
- var.blPar[4:var.bNodes[2]*var.nBlPar:var.nBlPar],
- var.blPar[5:var.bNodes[2]*var.nBlPar:var.nBlPar],
- frictionCoeff[:var.bNodes[2]], dissipCoeff[:var.bNodes[2]],
- var.u[4:var.bNodes[2]*var.nVar:var.nVar],
- var.blPar[3:var.bNodes[2]*var.nBlPar:var.nBlPar]]
-
- blVectors_wake=[var.blPar[var.bNodes[2]*var.nBlPar+9::var.nBlPar],
- var.u[var.bNodes[2]*var.nVar+0::var.nVar],
- var.u[var.bNodes[2]*var.nVar+1::var.nVar],
- var.blPar[var.bNodes[2]*var.nBlPar+6::var.nBlPar],
- var.blPar[var.bNodes[2]*var.nBlPar+4::var.nBlPar],
- var.blPar[var.bNodes[2]*var.nBlPar+5::var.nBlPar],
- frictionCoeff[var.bNodes[2]:], dissipCoeff[var.bNodes[2]:],
- var.u[var.bNodes[2]*var.nVar+4::var.nVar],
- var.blPar[var.bNodes[2]*var.nBlPar+3::var.nBlPar]]
-
- return blScalars, AFblwVel, WKblwVel, blVectors_airfoil, blVectors_wake, AFdeltaStarOut, WKdeltaStarOut
-
-
- def inverseMeshMapping(self, var, inMesh, inMeshbNodes):
- """Maps the variables required for blowing velocity computation from
- the fine mesh to the coarse mesh. Group separation is performed
- to avoid confusion between the points at the same place on the mesh
- but with different physical meaning (e.g First wake point).
- """
-
- fMapRhoUpper = interpolate.interp1d(var.xGlobal[:var.bNodes[1]].copy(), var.extPar[1:var.bNodes[1]*var.nExtPar:var.nExtPar].copy())
- rhoUpper = fMapRhoUpper(inMesh[:inMeshbNodes[1]])
- fMapRhoLower = interpolate.interp1d(var.xGlobal[var.bNodes[1]:var.bNodes[2]].copy(),
- var.extPar[var.bNodes[1]*var.nExtPar + 1:var.bNodes[2]*var.nExtPar:var.nExtPar].copy())
- rhoLower = fMapRhoLower(inMesh[inMeshbNodes[1]:inMeshbNodes[2]])
- fMapRhoWake = interpolate.interp1d(var.xGlobal[var.bNodes[2]:].copy(), var.extPar[var.bNodes[2]*var.nExtPar + 1::var.nExtPar].copy())
- rhoWake = fMapRhoWake(inMesh[inMeshbNodes[2]:])
-
- rho = np.concatenate((rhoUpper, rhoLower, rhoWake))
-
- fMapinvVelUpper = interpolate.interp1d(var.xGlobal[:var.bNodes[1]].copy(), var.u[3:var.bNodes[1]*var.nVar:var.nVar].copy())
- invVelUpper = fMapinvVelUpper(inMesh[:inMeshbNodes[1]])
- fMapinvVelLower = interpolate.interp1d(var.xGlobal[var.bNodes[1]:var.bNodes[2]].copy(),
- var.u[var.bNodes[1]*var.nVar + 3:var.bNodes[2]*var.nVar:var.nVar].copy())
- invVelLower = fMapinvVelLower(inMesh[inMeshbNodes[1]:inMeshbNodes[2]])
- fMapinvVelWake = interpolate.interp1d(var.xGlobal[var.bNodes[2]:].copy(), var.u[var.bNodes[2]*var.nVar + 3::var.nVar].copy())
- invVelWake = fMapinvVelWake(inMesh[inMeshbNodes[2]:])
-
- invVel = np.concatenate((invVelUpper, invVelLower, invVelWake))
-
- fMapdeltaStarUpper = interpolate.interp1d(var.xGlobal[:var.bNodes[1]].copy(), var.blPar[9:var.bNodes[1]*var.nBlPar:var.nBlPar].copy())
- deltaStarUpper = fMapdeltaStarUpper(inMesh[:inMeshbNodes[1]])
- fMapdeltaStarLower = interpolate.interp1d(var.xGlobal[var.bNodes[1]:var.bNodes[2]].copy(),
- var.blPar[var.bNodes[1]*var.nBlPar + 9:var.bNodes[2]*var.nBlPar:var.nBlPar].copy())
- deltaStarLower = fMapdeltaStarLower(inMesh[inMeshbNodes[1]:inMeshbNodes[2]])
- fMapdeltaStarWake = interpolate.interp1d(var.xGlobal[var.bNodes[2]:].copy(), var.blPar[var.bNodes[2]*var.nBlPar + 9::var.nBlPar].copy())
- deltaStarWake = fMapdeltaStarWake(inMesh[inMeshbNodes[2]:])
-
- deltaStar = np.concatenate((deltaStarUpper, deltaStarLower, deltaStarWake))
-
- return rho, invVel, deltaStar
+ def sortParam(self,var, mapMesh, cMeshDict, nFactor):
+
+ # Post process delta star
+
+ # Recompute deltaStar.
+ var.blPar[9::var.nBlPar] = var.u[0::var.nVar] * var.u[1::var.nVar]
+
+ # Map delta star to the coarse mesh and divide between the airfoil
+ # and wake related points
+ if mapMesh:
+
+ deltaStar = self.inverseMeshMapping(var.blPar[9::var.nBlPar], var.bNodes, cMeshDict['globCoord'], cMeshDict['bNodes'], nFactor)
+
+ AFdeltaStarOut = deltaStar[:cMeshDict['bNodes'][2]]
+ WKdeltaStarOut = deltaStar[cMeshDict['bNodes'][2]:]
+
+ else:
+
+ AFdeltaStarOut = var.blPar[9:var.bNodes[2]*var.nBlPar:var.nBlPar]
+ WKdeltaStarOut = var.blPar[var.bNodes[2]*var.nBlPar + 9::var.nBlPar]
+
+ # Compute blowing velocity on each cell.
+
+ blwVelUp = np.zeros(var.bNodes[1] - 1)
+ blwVelLw = np.zeros(var.bNodes[2] - var.bNodes[1])
+ blwVelWk = np.zeros(var.nN - var.bNodes[2] - 1)
+
+ # Upper
+ i = 0
+ for iPoint in range(1, var.bNodes[1]):
+
+ iPrev = iPoint - 1
+
+ blwVelUp[i] = (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]))
+
+ i += 1
+
+ # Lower
+ i = 0
+ for iPoint in range(var.bNodes[1], var.bNodes[2]):
+
+ iPrev = 0 if iPoint == var.bNodes[1] else iPoint - 1
+
+ blwVelLw[i] = (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]))
+
+ i += 1
+
+ # Wake
+ i = 0
+ for iPoint in range(var.bNodes[2] + 1, var.nN):
+
+ iPrev = iPoint - 1
+
+ blwVelWk[i] = (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]))
+
+ i += 1
+
+ if mapMesh:
+ blwVelUp, blwVelLw, blwVelWk = self.mapBlowingVelocities(blwVelUp, blwVelLw, blwVelWk, var.bNodes, cMeshDict['globCoord'], cMeshDict['bNodes'], var.xGlobal, nFactor)
+
+ # Postprocess the boundary layer parameters.
+
+ # Normalize Friction and dissipation coefficients.
+ frictionCoeff=var.u[3::var.nVar]**2 * var.blPar[2::var.nBlPar]
+ dissipCoeff=var.u[3::var.nVar]**2 * var.blPar[1::var.nBlPar]
+
+ # Compute total drag coefficient.
+ CdTot=2*var.u[-5]*(var.u[-2]**((var.u[-4]+5)/2))
+
+ # Compute friction and pressure drag coefficients.
+ Cdf_upper=np.trapz(frictionCoeff[:var.bNodes[1]],var.xGlobal[:var.bNodes[1]])
+ # Insert stagnation point in the lower side.
+ Cdf_lower=np.trapz(np.insert(frictionCoeff[var.bNodes[1]:var.bNodes[2]],0,frictionCoeff[0]),
+ np.insert(var.xGlobal[var.bNodes[1]:var.bNodes[2]],0,var.xGlobal[0]))
+ Cdf=Cdf_upper+Cdf_lower # No friction in the wake
+ Cdp=CdTot-Cdf
+
+ # [CdTot, Cdf, Cdp]
+ blScalars=[CdTot, Cdf, Cdp]
+
+ # Store boundary layer final parameters in lists.
+ # [deltaStar, theta, H, Hk, Hstar, Hstar2, Cf, Cd, Ct, delta]
+ blVectors_airfoil=[var.blPar[9:var.bNodes[2]*var.nBlPar:var.nBlPar],
+ var.u[0:var.bNodes[2]*var.nVar:var.nVar],
+ var.u[1:var.bNodes[2]*var.nVar:var.nVar],
+ var.blPar[6:var.bNodes[2]*var.nBlPar:var.nBlPar],
+ var.blPar[4:var.bNodes[2]*var.nBlPar:var.nBlPar],
+ var.blPar[5:var.bNodes[2]*var.nBlPar:var.nBlPar],
+ frictionCoeff[:var.bNodes[2]], dissipCoeff[:var.bNodes[2]],
+ var.u[4:var.bNodes[2]*var.nVar:var.nVar],
+ var.blPar[3:var.bNodes[2]*var.nBlPar:var.nBlPar]]
+
+ blVectors_wake=[var.blPar[var.bNodes[2]*var.nBlPar+9::var.nBlPar],
+ var.u[var.bNodes[2]*var.nVar+0::var.nVar],
+ var.u[var.bNodes[2]*var.nVar+1::var.nVar],
+ var.blPar[var.bNodes[2]*var.nBlPar+6::var.nBlPar],
+ var.blPar[var.bNodes[2]*var.nBlPar+4::var.nBlPar],
+ var.blPar[var.bNodes[2]*var.nBlPar+5::var.nBlPar],
+ frictionCoeff[var.bNodes[2]:], dissipCoeff[var.bNodes[2]:],
+ var.u[var.bNodes[2]*var.nVar+4::var.nVar],
+ var.blPar[var.bNodes[2]*var.nBlPar+3::var.nBlPar]]
+
+ return blScalars, blwVelUp, blwVelLw, blwVelWk, blVectors_airfoil, blVectors_wake, AFdeltaStarOut, WKdeltaStarOut
+
+ def mapBlowingVelocities(self, blwVelUp, blwVelLw, blwVelWk, fbNodes, cMesh, cMeshbNodes, fMesh, nFactor):
+ """Maps blowing velocties from cell to cell using weighted average interpolation.
+
+ Args
+ ----
+
+ - blwVelUp (numpy.ndarray): Blowing velocities on the fine mesh of the upper side.
+
+ - blwVelLw (numpy.ndarray): Blowing velocities on the fine mesh of the lower side.
+
+ - blwVelWk (numpy.ndarray): Blowing velocities on the fine mesh of the wake.
+
+ - fbNodes (numpy.ndarray): Boundary nodes of the fine mesh.
+
+ - cMesh (numpy.ndarray): Coarse mesh coordinates in the global frame of reference.
+
+ - cMeshbNodes (numpy.ndarray): Boundary nodes of the coarse mesh.
+
+ - fMesh (numpy.ndarray): Fine mesh coordinates in the global frame of reference.
+
+ - nFactor (int): Multiplicator underlying mesh adaptation.
+
+ Example
+ -------
+
+ Fine mesh: |-----------|-----------|-----------|----------|
+ \ /
+ (1/2) \ / (1/2)
+ \ /
+ \ /
+ Coarse mesh: |-----------------------|----------------------|
+ """
+
+ cblwVelUp = np.empty(len(cMesh[:cMeshbNodes[1]]) - 1)
+ cblwVelLw = np.empty(len(cMesh[cMeshbNodes[1]:cMeshbNodes[2]]))
+ cWKblwVel = np.empty(len(cMesh[cMeshbNodes[2]:]) - 1)
+
+ for cCell in range(len(cblwVelUp)):
+
+ cblwVelUp[cCell] = np.mean(blwVelUp[cCell*nFactor:(cCell + 1)*nFactor])
+
+ for cCell in range(len(cblwVelLw)):
+
+ cblwVelLw[cCell] = np.mean(blwVelLw[cCell*nFactor:(cCell + 1)*nFactor])
+
+ for cCell in range(len(cWKblwVel)):
+
+ cWKblwVel[cCell] = np.mean(blwVelWk[cCell*nFactor:(cCell + 1)*nFactor])
+
+ return cblwVelUp, cblwVelLw, cWKblwVel
+
+
+ def inverseMeshMapping(self, fVector, fbNodes, cMesh, cbNodes, nFactor):
+ """Maps the variables required for blowing velocity computation from
+ the fine mesh to the coarse mesh. Group separation is performed
+ to avoid confusion between the points at the same place on the mesh
+ but with different physical meaning (e.g First wake point).
+
+ The mapping is based on a weigthed average with the nearest neighbors.
+
+ Args
+ ----
+
+ - fVector (numpy.ndarray): Vector containing the function on the fine mesh.
+
+ - fbNodes (numpy.ndarray): Id of the boundary nodes on the fine mesh.
+
+ - cMesh (numpy.ndarray): Coarse mesh.
+
+ - cbNodes (numpy.ndarray): Id of the boundary nodes on the coarse mesh.
+
+ - nFactor (int): Multiplicator underlying mesh adaptation.
+
+ Return
+ ------
+
+ - cVector (numpy.ndarray): Function contained in 'fVector' mapped to the coarse mesh.
+
+ Example
+ -------
+
+ Fine mesh: |-----------|-----------|-----------|----------|
+ \ | /
+ (1/4) \ | (1/2) / (1/4)
+ \ V /
+ -----> <-----
+ Coarse mesh: |-----------------------|----------------------|
+ """
+
+ cVector = np.zeros(len(cMesh))
+
+ for cPoint in range(len(cMesh)-1):
+ cVector[cPoint] = fVector[cPoint * nFactor]
+ cVector[-1] = fVector[-1]
+ return cVector
diff --git a/dart/viscous/Physics/VBIConditions.py b/dart/viscous/Physics/VBIConditions.py
index c4cbf59..7f9efbb 100755
--- a/dart/viscous/Physics/VBIConditions.py
+++ b/dart/viscous/Physics/VBIConditions.py
@@ -158,7 +158,7 @@ class Boundaries:
#var.blPar[var.bNodes[0]*var.nBlPar+8]=(var.blPar[var.bNodes[0]*var.nBlPar+4]/2)*(1-4*(var.blPar[var.bNodes[0]*var.nBlPar+6]-1)/(3*H0)) # Equivalent normalized wall slip velocity
# Set boundary condition.
- var.u[var.bNodes[0]*var.nVar:(var.bNodes[0]+1)*var.nVar]=[T0, H0, n0, ue0, Ct0]
+ var.u[var.bNodes[0]*var.nVar:(var.bNodes[0]+1)*var.nVar] = [T0, H0, n0, ue0, Ct0]
pass
# Wake boundary conditions
diff --git a/dart/viscous/Solvers/VSolver.py b/dart/viscous/Solvers/VSolver.py
index 64a8e58..17f16c9 100644
--- a/dart/viscous/Solvers/VSolver.py
+++ b/dart/viscous/Solvers/VSolver.py
@@ -154,6 +154,13 @@ class flowEquations:
# Shear lag equation ignored in the laminar portion of the flow
timeMatrix[4] = [0, 0, 0, 0, 1]
+ """invTimeMatrix = np.empty((5,5))
+ invTimeMatrix[0] = [(uea**2 * (c*Ha) * uea**3 - (Ha*Hstara - 2)*uea**2 + Ha*(Hstara-1))/(Ha*(c*Ha*Ta+uea)) + (Ha*Hstara - 2)*uea**2/Ha, uea, 0, (Ta*((Ha*Hstara - 2)*uea**2 - Ha * (Hstara - 1))+uea**4)/(c*Ha*Ta + uea), 0]
+ invTimeMatrix[1] = [(-(c*Ha*(Ha*Hstara-2)*Ta*uea + c*Ha*uea**3 + Ha*(Hstara - 1) - 1)*uea**2)/(Ta*(c*Ha*Ta + uea)), -Ha*uea/Ta, 0, (-Ha*(Ta*((Ha*Hstara - 2)*uea**2 - Ha*(Hstara - 1) + 1) + uea**4))/(Ta*(c*Ha*Ta + uea)), 0]
+ invTimeMatrix[2] = [0, 0, 1, 0, 0]
+ invTimeMatrix[3] = [(c*uea**2)/(c*Ha*Ta + uea), 0, 0, uea/(c*Ha*Ta + uea), 0]
+ invTimeMatrix[4] = [-c*Cta*uea/(c*Ha*Ta+uea), 0, 0, -Cta/(c*Ha*Ta + uea), (0.5*Cta*Usa*uea)/deltaa]"""
+
sysRes = np.linalg.solve(timeMatrix, spaceMatrix).flatten()
return sysRes
diff --git a/dart/viscous/Solvers/VTimeSolver.py b/dart/viscous/Solvers/VTimeSolver.py
index 3e6de63..8643dcf 100644
--- a/dart/viscous/Solvers/VTimeSolver.py
+++ b/dart/viscous/Solvers/VTimeSolver.py
@@ -23,430 +23,475 @@ import numpy as np
import math
class timeConfig:
- """ ---------- Time integration configuration ----------
-
- Time integration related computations: CFL adaptation, time step (per cell)
-
- Attributes
- ----------
- cflAdapt: bool
- 'True' if CFL adaptation is active, 'False' otherwise
-
- __cflAdaptParam: numpy vect
- [Lower, Upper] bounds of the CFL
- """
- def __init__(self, _Minfty):
- if _Minfty != 0:
- self.__Minf = _Minfty
- else:
- self.__Minf = 0.1
- self.gamma = 1.4
- pass
-
- def getSoundSpeed(self, gradU2):
- return np.sqrt(1 / (self.__Minf * self.__Minf) + 0.5 * (self.gamma - 1) * (1 - gradU2))
-
- def computeTimeStep(self, CFL, dx, invVel):
-
- # Speed of sound.
- gradU2 = (invVel)**2
- # Velocity of the fastest travelling wave
- soundSpeed = self.getSoundSpeed(gradU2)
- advVel = soundSpeed + invVel
-
- # Time step computation
- return CFL*dx/advVel
-
-
- def adaptCFL(self, resCurr, res0):
- """CFL number adaptation. The metric driving the adaptation is the residual value.
- """
- CFL=self.__CFL0*(res0/resCurr)**0.7
- CFL=max(CFL,self.__CFL0)
- return CFL
-
- def outputState(self, var, iMarker, innerIters, normLinSysRes, normLinSysRes0, CFL, color='white'):
- if color == 'white':
- print('| {:>6.0f}| {:>11.0f}| {:>9.4f}| {:>10.2f}| {:>14.6f}| {:>6.0f}| {:>11.3f}|'.format(iMarker,
- innerIters,
- np.log10(normLinSysRes/normLinSysRes0),
- CFL,
- var.xGlobal[iMarker],
- var.flowRegime[iMarker],
- var.u[iMarker*var.nVar + 2]))
- elif color == 'yellow':
- print(ccolors.ANSI_YELLOW + '| {:>6.0f}| {:>11.0f}| {:>9.4f}| {:>8.2f}| {:>14.6f}| {:>6.0f}| {:>11.3f}|'.format(iMarker,
- innerIters,
- np.log10(normLinSysRes/normLinSysRes0),
- CFL,
- var.xGlobal[iMarker],
- var.flowRegime[iMarker],
- var.u[iMarker*var.nVar + 2]),
- ccolors.ANSI_RESET)
- elif color == 'green':
- print(ccolors.ANSI_YELLOW + '| {:>6.0f}| {:>11.0f}| {:>9.4f}| {:>8.2f}| {:>14.6f}| {:>6.0f}| {:>11.3f}|'.format(iMarker,
- innerIters,
- np.log10(normLinSysRes/normLinSysRes0),
- CFL,
- var.xGlobal[iMarker],
- var.flowRegime[iMarker],
- var.u[iMarker*var.nVar + 2]),
- ccolors.ANSI_RESET)
+ """ ---------- Time integration configuration ----------
+
+ Time integration related computations: CFL adaptation, time step (per cell)
+
+ Attributes
+ ----------
+ cflAdapt: bool
+ 'True' if CFL adaptation is active, 'False' otherwise
+
+ __cflAdaptParam: numpy vect
+ [Lower, Upper] bounds of the CFL
+ """
+ def __init__(self, _Minfty):
+
+ # Correct incompressible Mach number.
+ if _Minfty != 0: self.__Minf = _Minfty
+ else: self.__Minf = 0.1
+
+ self.gamma = 1.4
+ pass
+
+ def getSoundSpeed(self, gradU2):
+ return np.sqrt(1 / (self.__Minf * self.__Minf) + 0.5 * (self.gamma - 1) * (1 - gradU2))
+
+ def computeTimeStep(self, CFL, dx, invVel):
+
+ # Speed of sound.
+ gradU2 = (invVel)**2
+ # Velocity of the fastest travelling wave
+ soundSpeed = self.getSoundSpeed(gradU2)
+ advVel = soundSpeed + invVel
+
+ # Time step computation
+ return CFL*dx/advVel
+
+
+ def adaptCFL(self, resCurr, res0):
+ """CFL number adaptation. The metric driving the adaptation is the residual value.
+ """
+ CFL=self.__CFL0*(res0/resCurr)**0.7
+ CFL=max(CFL,self.__CFL0)
+ return CFL
+
+ def outputState(self, var, iMarker, innerIters, normLinSysRes, normLinSysRes0, CFL, color='white'):
+ if color == 'white':
+ print('| {:>6.0f}| {:>11.0f}| {:>9.4f}| {:>10.2f}| {:>14.6f}| {:>6.0f}| {:>11.3f}|'
+ .format(iMarker,
+ innerIters,
+ np.log10(normLinSysRes/normLinSysRes0),
+ CFL,
+ var.xGlobal[iMarker],
+ var.flowRegime[iMarker],
+ var.u[iMarker*var.nVar + 2]))
+ elif color == 'yellow':
+ print(ccolors.ANSI_YELLOW + '| {:>6.0f}| {:>11.0f}| {:>9.4f}| {:>8.2f}| {:>14.6f}| {:>6.0f}| {:>11.3f}|'
+ .format(iMarker,
+ innerIters,
+ np.log10(normLinSysRes/normLinSysRes0),
+ CFL,
+ var.xGlobal[iMarker],
+ var.flowRegime[iMarker],
+ var.u[iMarker*var.nVar + 2]),
+ ccolors.ANSI_RESET)
+ elif color == 'green':
+ print(ccolors.ANSI_YELLOW + '| {:>6.0f}| {:>11.0f}| {:>9.4f}| {:>8.2f}| {:>14.6f}| {:>6.0f}| {:>11.3f}|'
+ .format(iMarker,
+ innerIters,
+ np.log10(normLinSysRes/normLinSysRes0),
+ CFL,
+ var.xGlobal[iMarker],
+ var.flowRegime[iMarker],
+ var.u[iMarker*var.nVar + 2]),
+ ccolors.ANSI_RESET)
class implicitEuler(timeConfig):
- """ ---------- Implicit Euler time integration ----------
-
- Solves the unsteady boundary layer equations using the implicit Euler
- method with Newton method to solve the non-linear system.
-
- Attributes
- ----------
- solver: class
- Spacial operator and Jacobian computation of the boundary layer laws
-
- timeParams: dict
- Contains necessary information for the solver: tolerance, CFL, CFL apdation parameters
-
- safeguard: int
- Number of iterations for which the residual has to decrease to unlock CFL adaptation
- (usually changes after the first coupling iteration)
- and second order Jacobian evaluation starts.
-
- safemode: bool
- Indicates if the solver is in safe mode (1) or not (0). Value switched according to safeguard
- """
-
- def __init__(self, _solver, transition, wakeBC, _cfl0, _Minfty):
- timeConfig.__init__(self, _Minfty)
- # Initialize time parameters
-
- self.__CFL0=_cfl0
- self.__maxIt = 15000
- self.__tol = 1e-8
- self.jacEval0 = 5
-
- self.solver=_solver
- self.transSolver = transition
- self.wakeBC = wakeBC
- def __str__(self):
- return 'Implicit Euler'
-
- def getCFL0(self): return self.__CFL0
- def setCFL0(self, _cfl0): self.__CFL0 = _cfl0
-
- def gettol(self): return self.__tol
- def settol(self, _tol): self.__tol = _tol
-
- def getmaxIt(self): return self.__maxIt
- def setmaxIt(self, _maxIt): self.__maxIt = _maxIt
-
- def resetSolution(self, iMarker, var, u0):
-
- if iMarker == var.bNodes[1]:
-
- var.u[iMarker*var.nVar : (iMarker + 1)*var.nVar] = var.u[0*var.nVar : 1*var.nVar]
-
- elif iMarker != var.transMarkers[0] and iMarker != var.transMarkers[1]:
-
- var.u[iMarker*var.nVar : (iMarker + 1)*var.nVar] = var.u[(iMarker-1)*var.nVar : (iMarker)*var.nVar]
-
- else:
-
- var.u[iMarker*var.nVar : (iMarker + 1)*var.nVar] = u0[iMarker*var.nVar : (iMarker + 1)*var.nVar].copy()
-
- def timeSolve(self, var):
- """ ---------- Newton iteration to solve the boundary layer equations ----------
-
- Damped Newton method in Pseudo-Time. Linear algebraic system (obtained by Taylor linearis.)
- is solved with a direct method.
-
- Equations
- ---------
- ∂U/∂t = F(U)
- <=> (U_(n+1)-U_n)/∆t = - F(U_(n+1)) (Time disc.)
- <=> ∆U/∆t = -(F(U_n) + ∂F/∂U*∆u), where ∂F/∂U = J(U_n)
- <=> (I+J)dU = -F(U_n), J=dF_n/dU_n (Linearis.)
-
- Parameters
- ----------
- - var: class
- Data container: - Solution
- - Boundary layer parameters
- - External flow parameters
- - Mesh
- - Transition information
-
- Tasks
- ------
- - Updates the solution 'var.u' through Newton iterations until steady state.
-
- - Asks 'solver' for spacial operator F and Jacobian J and uses a
- direct linear solver (LU decomposition).
-
- - Corrects undesirable results.
-
- - Asks the solver for transtion position/BC and wake BC to be updated at each iteration.
- """
-
-
- nN = var.nN # Number of points in the computational domain.
- nBlPar = var.nBlPar # Number of boundary layer parameters used for vector indexing.
- nExtPar = var.nExtPar # Number of boundary layer parameters used for vector indexing.
- bNodes = var.bNodes # Boundary nodes of the computational domain.
- convergedPoints = np.zeros(var.nN) # Convergence control of each point.
- convergedPoints[0] = 1 # Boundary condition.
-
- # Initialize the flow regime on each cell.
- self.transSolver.initTransition(var)
-
- lockTrans = 0 # Flag to remember if the transition is already found or not.
-
- for iMarker in range(1, nN):
- displayState = 0 # Flag to output simulation state.
-
- if not lockTrans and 1 < iMarker < bNodes[2]:
-
- # Determine if the amplification waves start growing or not.
- logRet_crit = 2.492*(1/(var.blPar[(iMarker-1)*nBlPar + 6]-1))**0.43
- + 0.7*(np.tanh(14*(1/(var.blPar[(iMarker-1)*nBlPar + 6]-1))-9.24)+1.0)
-
- if var.blPar[(iMarker-1)*nBlPar + 0] > 0: # Avoid log of something <= 0
-
- if np.log10(var.blPar[(iMarker-1)*nBlPar + 0]) <= logRet_crit - 0.08:
-
- var.u[iMarker*var.nVar + 2] = 0
-
- # Upper side start.
- if iMarker == bNodes[0]+1:
- print('| - Computing upper side. {:>54}'.format('|'))
-
- # Lower side start.
- if iMarker == bNodes[1]:
- lockTrans = 0
- print('| - Computing lower side. {:>54}'.format('|'))
-
- # Wake start
- if iMarker == bNodes[2]: # First wake point
-
- # Wake boundary condition.
- self.wakeBC(var)
- convergedPoints[iMarker] = 1
- lockTrans = 1
- print('| - Computing wake. {:>60}'.format('|'))
-
- # Ignore remaining instructions for the first wake point.
- continue
-
- if not lockTrans and iMarker-1 < bNodes[2] and var.u[(iMarker-1) * var.nVar +2 ] >= var.Ncrit:
- # Transition
- print('| Transition located near x/c = {:<2.3f}. Marker: {:<4.0f} {:>28s}'.format(var.xGlobal[iMarker-1],iMarker-1, '|'))
-
- avTurb = self.transSolver.solveTransition(var, iMarker - 1)
- # Save laminar solution
- lamSol = var.u[(iMarker-1)*var.nVar : (iMarker)* var.nVar].copy()
-
- # Set up turbulent portion boundary condition
- avTurb = self.transSolver.solveTransition(var, iMarker - 1)
-
- # Compute turbulent solution
- var.flowRegime[iMarker - 1] = 1
- self.newtonSolver(var, iMarker-1, displayState)
-
- # Average solutions
- var.u[(iMarker-1)*var.nVar : (iMarker)*var.nVar] = (1-avTurb) * lamSol + avTurb * var.u[(iMarker-1)*var.nVar : (iMarker)*var.nVar]
-
- # Recompute closures
- var.blPar[(iMarker-1)*nBlPar+1:(iMarker-1)*nBlPar+9]=Closures.turbulentClosure(var.u[(iMarker-1)*var.nVar+0],
- var.u[(iMarker-1)*var.nVar+1],
- var.u[(iMarker-1)*var.nVar+4],
- var.blPar[(iMarker-1)*nBlPar+0],
- var.extPar[(iMarker-1)*nExtPar+0],
- 'airfoil')
-
- lockTrans = 1
-
- # Forced transition
- if not lockTrans and var.xtrF[0] is not None and (var.flowRegime[iMarker] == 1 and var.flowRegime[iMarker-1] == 0):
- print('| Transition near {:<1.4f}. Marker: {:<5.0f} {:>40}'.format(var.xGlobal[iMarker-1],iMarker-1, '|'))
- lamSol = var.u[(iMarker-1)*var.nVar : (iMarker)* var.nVar].copy()
- avTurb = self.transSolver.transitionBC(var, iMarker - 1)
- # Compute turbulent solution
- var.flowRegime[iMarker - 1] = 1
- self.newtonSolver(var, iMarker-1, displayState)
- # Average solutions
- var.u[(iMarker-1)*var.nVar : (iMarker)*var.nVar] = (1-avTurb) * lamSol + avTurb * var.u[(iMarker-1)*var.nVar : (iMarker)*var.nVar]
- # Recompute closures
- var.blPar[(iMarker-1)*nBlPar+1:(iMarker-1)*nBlPar+9]=Closures.turbulentClosure(var.u[(iMarker-1)*var.nVar+0],
- var.u[(iMarker-1)*var.nVar+1],
- var.u[(iMarker-1)*var.nVar+4],
- var.blPar[(iMarker-1)*nBlPar+0],
- var.extPar[(iMarker-1)*nExtPar+0],
- 'airfoil')
-
- lockTrans = 1
- if not lockTrans and var.xtrF[1] is not None and (var.flowRegime[iMarker]== 1 and var.flowRegime[iMarker-1] == 0):
- print('| Transition near {:<1.4f}. Marker: {:<5.0f} {:>40}'.format(var.xGlobal[iMarker-1],iMarker-1, '|'))
- lamSol = var.u[(iMarker-1)*var.nVar : (iMarker)* var.nVar].copy()
- avTurb = self.transSolver.transitionBC(var, iMarker - 1)
- # Compute turbulent solution
- var.flowRegime[iMarker - 1] = 1
- self.newtonSolver(var, iMarker-1, displayState)
- # Average solutions
- var.u[(iMarker-1)*var.nVar : (iMarker)*var.nVar] = (1-avTurb) * lamSol + avTurb * var.u[(iMarker-1)*var.nVar : (iMarker)*var.nVar]
- # Recompute closures
- var.blPar[(iMarker-1)*nBlPar+1:(iMarker-1)*nBlPar+9]=Closures.turbulentClosure(var.u[(iMarker-1)*var.nVar+0],
- var.u[(iMarker-1)*var.nVar+1],
- var.u[(iMarker-1)*var.nVar+4],
- var.blPar[(iMarker-1)*nBlPar+0],
- var.extPar[(iMarker-1)*nExtPar+0],
- 'airfoil')
-
- # Call Newton solver for the point.
- convergedPoints[iMarker] = self.newtonSolver(var, iMarker, displayState)
-
- return convergedPoints
-
-
- def newtonSolver (self, var, iMarker, displayState):
-
- # Save initial condition in case a point diverges.
- u0=var.u.copy() # Initial solution.
-
- nVar = var.nVar # Number of variables of the differential problem.
- nBlPar = var.nBlPar # Number of boundary layer parameters used for vector indexing.
- nExtPar = var.nExtPar # Number of boundary layer parameters used for vector indexing.
- bNodes = var.bNodes # Boundary nodes of the computational domain.
- iInvVel = var.extPar[iMarker*nExtPar + 2] # Inviscid velocity on the considered cell.
- iMarkerPrev = 0 if iMarker == bNodes[1] else iMarker - 1 # Point upstream of the considered point.
- dx = var.xx[iMarker] - var.xx[iMarkerPrev] # Cell size.
-
- sysRes = self.solver.buildResiduals(var, iMarker) # System initial residuals.
- normSysRes0 = np.linalg.norm(sysRes) # Initial norm of the residuals.
-
- CFL = self.__CFL0 # Initialized CFL number.
- CFLAdapt = 1 # Flag for CFL adaptation.
- 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.
-
- # Main loop.
- converged = 0 # Flag to output convergence or divergence of the point.
- nErrors = 0 # Counter for the number of errors identified at that point.
- innerIters = 0 # Number of inner iterations to solver the non-linear system.
- while innerIters <= self.__maxIt:
-
- try:
-
- # Jacobian computation.
-
- if innerIters % jacEval == 0:
- Jacobian=self.solver.buildJacobian(var, iMarker, sysRes, order = 1, eps = 1e-8)
-
- # 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.
-
- var.u[iMarker*nVar : (iMarker+1)*nVar] += slnIncr
- # var.u[iMarker*nVar + 1] = max(var.u[iMarker*nVar + 1], 1.0005)
-
- # Compute Residuals.
- sysRes = self.solver.buildResiduals(var, iMarker)
- normSysRes = np.linalg.norm(slnIncr/dt - sysRes)
-
- if displayState and innerIters % 1000 == 0:
- self.outputState(var, iMarker, innerIters, normSysRes, normSysRes0, CFL, 'yellow')
- # Check convergence.
- if normSysRes <= self.__tol * normSysRes0:
- converged = 1
- if displayState:
- self.outputState(var, iMarker, innerIters, normSysRes, normSysRes0, CFL)
- break
-
- except KeyboardInterrupt:
- quit()
- except Exception as error:
-
- nErrors += 1
- if nErrors >= 5:
- print('|', ccolors.ANSI_RED + 'Too many errors identified at position {:<1.4f}. Marker: {:<4.0f}. {:>17s}'
- .format(var.xGlobal[iMarker], iMarker,'|'), ccolors.ANSI_RESET)
- quit()
-
- print('|', ccolors.ANSI_YELLOW + 'Newton solver failed at position {:<.4f}. Marker: {:<.0f} {:>25s}'
- .format(var.xGlobal[iMarker], iMarker,'|'), ccolors.ANSI_RESET)
-
- print(error)
-
- #self.resetSolution(iMarker, var, u0)
- var.u[iMarker*nVar : (iMarker+1)*nVar] = var.u[(iMarker-1)*nVar : (iMarker)*nVar]
-
- # Lower CFL and recompute time step
- print('Current CFL', CFL)
- if math.isnan(CFL):
- CFL = self.__CFL0
- CFL = min(max(CFL / 2,0.01),1)
- print('Lowering CFL: {:<.3f}'.format(CFL))
- CFLAdapt = 1
- dt = self.computeTimeStep(CFL, dx, iInvVel)
- IMatrix = np.identity(nVar) / dt
-
- sysRes = self.solver.buildResiduals(var, iMarker)
-
- jacEval = 1
- displayState = 1
- print('+------------------------------------------------------------------------------+')
- print('| {:>6s}| {:>8s}| {:>9s}| {:>8s}| {:>12s}| {:>6s}| {:>8s}|'.format('Point','Inner Iters',
- 'rms[F]', 'End CFL',
- 'Position (x/c)', 'Regime',
- 'Amp. factor'))
- print('+------------------------------------------------------------------------------+')
- continue
-
- # CFL adaptation
- if CFLAdapt:
-
- CFL = self.__CFL0* (normSysRes0/normSysRes)**0.7
-
- CFL = max(CFL, 0.01)
- dt = self.computeTimeStep(CFL, dx, iInvVel)
- IMatrix=np.identity(nVar) / dt
-
-
- innerIters+=1
-
- else:
- # Maximum number of iterations reached.
- if normSysRes >= 1e-3 * normSysRes0:
- 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()
- name = 'airfoil' if iMarker <= var.bNodes[2] else 'wake'
-
- if var.flowRegime[iMarker]==0:
-
- # Laminar closure relations
- var.blPar[iMarker*nBlPar+1:iMarker*nBlPar+7] = Closures.laminarClosure(var.u[iMarker*nVar + 0], var.u[iMarker*nVar + 1],
- var.blPar[iMarker*nBlPar+0], var.extPar[iMarker*nExtPar+0],
- name)
-
- elif var.flowRegime[iMarker]==1:
-
- # Turbulent closures relations
- var.blPar[iMarker*nBlPar+1:iMarker*nBlPar+9] = Closures.turbulentClosure(var.u[iMarker*nVar + 0], var.u[iMarker*nVar + 1],
- var.u[iMarker*nVar + 4], var.blPar[iMarker*nBlPar+0],
- var.extPar[iMarker*nExtPar+0], name)
-
- else:
- 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)
-
- return converged
\ No newline at end of file
+ """ ---------- Implicit Euler time integration ----------
+
+ Solves the unsteady boundary layer equations using the implicit Euler
+ method with Newton method to solve the non-linear system.
+
+ Attributes
+ ----------
+ solver: class
+ Spacial operator and Jacobian computation of the boundary layer laws
+
+ timeParams: dict
+ Contains necessary information for the solver: tolerance, CFL, CFL apdation parameters
+
+ safeguard: int
+ Number of iterations for which the residual has to decrease to unlock CFL adaptation
+ (usually changes after the first coupling iteration)
+ and second order Jacobian evaluation starts.
+
+ safemode: bool
+ Indicates if the solver is in safe mode (1) or not (0). Value switched according to safeguard
+ """
+
+ def __init__(self, _solver, transition, wakeBC, _cfl0, _Minfty):
+ timeConfig.__init__(self, _Minfty)
+ # Initialize time parameters
+
+ self.__CFL0=_cfl0
+ self.__maxIt = 15000
+ self.__tol = 1e-8
+ self.jacEval0 = 5
+ self.initSln = 0
+
+ self.solver=_solver
+ self.transSolver = transition
+ self.wakeBC = wakeBC
+ def __str__(self):
+ return 'Implicit Euler'
+
+ def getCFL0(self): return self.__CFL0
+ def setCFL0(self, _cfl0): self.__CFL0 = _cfl0
+
+ def gettol(self): return self.__tol
+ def settol(self, _tol): self.__tol = _tol
+
+ def getmaxIt(self): return self.__maxIt
+ def setmaxIt(self, _maxIt): self.__maxIt = _maxIt
+
+ def resetSolution(self, iMarker, var, u0):
+
+ if iMarker == var.bNodes[1]:
+
+ var.u[iMarker*var.nVar : (iMarker + 1)*var.nVar] = var.u[0*var.nVar : 1*var.nVar]
+
+ elif iMarker != var.transMarkers[0] and iMarker != var.transMarkers[1]:
+
+ var.u[iMarker*var.nVar : (iMarker + 1)*var.nVar] = var.u[(iMarker-1)*var.nVar : (iMarker)*var.nVar]
+
+ else:
+
+ var.u[iMarker*var.nVar : (iMarker + 1)*var.nVar] = u0[iMarker*var.nVar : (iMarker + 1)*var.nVar].copy()
+
+ def timeSolve(self, var):
+ """ ---------- Newton iteration to solve the boundary layer equations ----------
+
+ Damped Newton method in Pseudo-Time. Linear algebraic system (obtained by Taylor linearis.)
+ is solved with a direct method.
+
+ Equations
+ ---------
+ ∂U/∂t = F(U)
+ <=> (U_(n+1)-U_n)/∆t = - F(U_(n+1)) (Time disc.)
+ <=> ∆U/∆t = -(F(U_n) + ∂F/∂U*∆u), where ∂F/∂U = J(U_n)
+ <=> (I+J)dU = -F(U_n), J=dF_n/dU_n (Linearis.)
+
+ Parameters
+ ----------
+ - var: class
+ Data container: - Solution
+ - Boundary layer parameters
+ - External flow parameters
+ - Mesh
+ - Transition information
+
+ Tasks
+ ------
+ - Updates the solution 'var.u' through Newton iterations until steady state.
+
+ - Asks 'solver' for spacial operator F and Jacobian J and uses a
+ direct linear solver (LU decomposition).
+
+ - Corrects undesirable results.
+
+ - Asks the solver for transtion position/BC and wake BC to be updated at each iteration.
+ """
+
+
+ nN = var.nN # Number of points in the computational domain.
+ nBlPar = var.nBlPar # Number of boundary layer parameters used for vector indexing.
+ nExtPar = var.nExtPar # Number of boundary layer parameters used for vector indexing.
+ bNodes = var.bNodes # Boundary nodes of the computational domain.
+ convergedPoints = np.zeros(var.nN) # Convergence control of each point.
+ convergedPoints[0] = 1 # Boundary condition.
+
+ self.transSolver.initTransition(var) # Initialize the flow regime on each cell.
+
+ lockTrans = 0 # Flag to remember if the transition is already found or not.
+
+ for iMarker in range(1, nN):
+ displayState = 0 # Flag to output simulation state.
+
+ if self.initSln:
+ iMarkerPrev = 0 if iMarker == var.bNodes[1] else iMarker - 1
+ var.u[iMarker*var.nVar: (iMarker + 1)*var.nVar] = var.u[iMarkerPrev*var.nVar: (iMarkerPrev + 1)* var.nVar]
+
+ if not lockTrans and 1 < iMarker < bNodes[2]:
+
+ self.__checkWaves(var, iMarker)
+
+ # Upper side start.
+ if iMarker == bNodes[0]+1:
+ print('| - Computing upper side. {:>54}'.format('|'))
+
+ # Lower side start.
+ if iMarker == bNodes[1]:
+ lockTrans = 0
+ print('| - Computing lower side. {:>54}'.format('|'))
+
+ # Wake start
+ if iMarker == bNodes[2]: # First wake point
+
+ # Wake boundary condition.
+ self.wakeBC(var)
+ convergedPoints[iMarker] = 1
+ lockTrans = 1
+ print('| - Computing wake. {:>60}'.format('|'))
+
+ # Ignore remaining instructions for the first wake point.
+ continue
+
+ # Call Newton solver for the point.
+ convergedPoints[iMarker] = self.newtonSolver(var, iMarker, displayState)
+
+ # Check for transition.
+ if not lockTrans:
+
+ # Free transition.
+ if var.u[iMarker * var.nVar + 2] >= var.Ncrit:
+ print('| Transition located near x/c = {:<2.3f}. Marker: {:<4.0f} {:>28s}'.format(var.xGlobal[iMarker],iMarker, '|'))
+ self.__averageTransition(var, iMarker, displayState)
+ lockTrans = 1
+
+ # Forced transition.
+ elif var.xtrF[0] is not None or var.xtrF[1] is not None:
+ if var.flowRegime[iMarker] == 1 and var.flowRegime[iMarker - 1]:
+ print('| Forced transition near x/c = {:<2.3f}. Marker: {:<4.0f} {:>28s}'.format(var.xGlobal[iMarker],iMarker, '|'))
+ self.__forcedTransition(var, iMarker)
+ lockTrans = 1
+
+ return convergedPoints
+
+
+ def newtonSolver (self, var, iMarker, displayState):
+
+ # Save initial condition in case a point diverges.
+ u0=var.u.copy() # Initial solution.
+
+ nVar = var.nVar # Number of variables of the differential problem.
+ nBlPar = var.nBlPar # Number of boundary layer parameters used for vector indexing.
+ nExtPar = var.nExtPar # Number of boundary layer parameters used for vector indexing.
+ bNodes = var.bNodes # Boundary nodes of the computational domain.
+ iInvVel = var.extPar[iMarker*nExtPar + 2] # Inviscid velocity on the considered cell.
+ iMarkerPrev = 0 if iMarker == bNodes[1] else iMarker - 1 # Point upstream of the considered point.
+ dx = var.xx[iMarker] - var.xx[iMarkerPrev] # Cell size.
+
+ sysRes = self.solver.buildResiduals(var, iMarker) # System initial residuals.
+ normSysRes0 = np.linalg.norm(sysRes) # Initial norm of the residuals.
+
+ CFL = self.__CFL0 # Initialized CFL number.
+ CFLAdapt = 1 # Flag for CFL adaptation.
+ 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.
+
+ # Main loop.
+ nErrors = 0 # Counter for the number of errors identified at that point.
+ innerIters = 0 # Number of inner iterations to solver the non-linear system.
+ while innerIters <= self.__maxIt:
+
+ try:
+
+ # Jacobian computation.
+
+ if innerIters % jacEval == 0:
+ Jacobian=self.solver.buildJacobian(var, iMarker, sysRes, order = 1, eps = 1e-8)
+
+ # 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.
+
+ var.u[iMarker*nVar : (iMarker+1)*nVar] += slnIncr
+ # var.u[iMarker*nVar + 1] = max(var.u[iMarker*nVar + 1], 1.0005)
+
+ # Compute Residuals.
+ sysRes = self.solver.buildResiduals(var, iMarker)
+ normSysRes = np.linalg.norm(slnIncr/dt - sysRes)
+
+ if displayState and innerIters > 0 and innerIters % 1000 == 0:
+ self.outputState(var, iMarker, innerIters, normSysRes, normSysRes0, CFL, 'yellow')
+ # Check convergence.
+ if normSysRes <= self.__tol * normSysRes0:
+ converged = 1
+ if displayState:
+ self.outputState(var, iMarker, innerIters, normSysRes, normSysRes0, CFL)
+ return 1
+
+ except KeyboardInterrupt:
+ quit()
+ except Exception as error:
+
+ nErrors += 1
+ if nErrors >= 5:
+ print('|', ccolors.ANSI_RED + 'Too many errors identified at position {:<1.4f}. Marker: {:<4.0f}. {:>17s}'
+ .format(var.xGlobal[iMarker], iMarker,'|'), ccolors.ANSI_RESET)
+ quit()
+
+ print('|', ccolors.ANSI_YELLOW + 'Newton solver failed at position {:<.4f}. Marker: {:<.0f} {:>25s}'
+ .format(var.xGlobal[iMarker], iMarker,'|'), ccolors.ANSI_RESET)
+
+ print(error)
+
+ #self.resetSolution(iMarker, var, u0)
+ var.u[iMarker*nVar : (iMarker+1)*nVar] = var.u[(iMarker-1)*nVar : (iMarker)*nVar]
+
+ # Lower CFL and recompute time step
+ print('Current CFL', CFL)
+ if math.isnan(CFL):
+ CFL = self.__CFL0
+ CFL = min(max(CFL / 2,0.01),1)
+ print('Lowering CFL: {:<.3f}'.format(CFL))
+ CFLAdapt = 1
+ dt = self.computeTimeStep(CFL, dx, iInvVel)
+ IMatrix = np.identity(nVar) / dt
+
+ sysRes = self.solver.buildResiduals(var, iMarker)
+
+ jacEval = 1
+ displayState = 1
+ print('+------------------------------------------------------------------------------+')
+ print('| {:>6s}| {:>8s}| {:>9s}| {:>8s}| {:>12s}| {:>6s}| {:>8s}|'.format('Point','Inner Iters',
+ 'rms[F]', 'End CFL',
+ 'Position (x/c)', 'Regime',
+ 'Amp. factor'))
+ print('+------------------------------------------------------------------------------+')
+ continue
+
+ # CFL adaptation
+ if CFLAdapt:
+
+ CFL = self.__CFL0* (normSysRes0/normSysRes)**0.7
+
+ CFL = max(CFL, 0.01)
+ dt = self.computeTimeStep(CFL, dx, iInvVel)
+ IMatrix=np.identity(nVar) / dt
+
+
+ innerIters+=1
+
+ else:
+ # Maximum number of iterations reached.
+ if normSysRes >= 1e-3 * normSysRes0:
+ print('|', ccolors.ANSI_RED + 'Too many iterations 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()
+ name = 'airfoil' if iMarker <= var.bNodes[2] else 'wake'
+
+ if var.flowRegime[iMarker]==0:
+
+ # Laminar closure relations
+ var.blPar[iMarker*nBlPar+1:iMarker*nBlPar+7] = Closures.laminarClosure(var.u[iMarker*nVar + 0], var.u[iMarker*nVar + 1],
+ var.blPar[iMarker*nBlPar+0], var.extPar[iMarker*nExtPar+0],
+ name)
+
+ elif var.flowRegime[iMarker]==1:
+
+ # Turbulent closures relations
+ var.blPar[iMarker*nBlPar+1:iMarker*nBlPar+9] = Closures.turbulentClosure(var.u[iMarker*nVar + 0], var.u[iMarker*nVar + 1],
+ var.u[iMarker*nVar + 4], var.blPar[iMarker*nBlPar+0],
+ var.extPar[iMarker*nExtPar+0], name)
+
+ else:
+ print('|', ccolors.ANSI_YELLOW+ 'Too many iterations at position {:<.4f}. Marker: {:>5.0f} normLinSysRes = {:<.3f}. {:>30s}'
+ .format(var.xGlobal[iMarker], iMarker,
+ np.log10(normSysRes/normSysRes0),'|'),
+ ccolors.ANSI_RESET)
+ return 0
+
+ def __checkWaves(self, var, iMarker):
+ """Determine if the amplification waves start growing or not.
+ """
+
+ logRet_crit = 2.492*(1/(var.blPar[(iMarker-1)*var.nBlPar + 6]-1))**0.43
+ + 0.7*(np.tanh(14*(1/(var.blPar[(iMarker-1)*var.nBlPar + 6]-1))-9.24)+1.0)
+
+ if var.blPar[(iMarker-1)*var.nBlPar + 0] > 0: # Avoid log of something <= 0
+ if np.log10(var.blPar[(iMarker-1)*var.nBlPar + 0]) <= logRet_crit - 0.08:
+ var.u[iMarker*var.nVar + 2] = 0
+
+ pass
+
+
+ def __averageTransition(self, var, iMarker, displayState):
+ """Averages laminar and turbulent solution at the transition location.
+ The boundary condition of the turbulent flow portion is also imposed.
+ """
+
+ # Save laminar solution
+ lamSol = var.u[(iMarker)*var.nVar : (iMarker+1)* var.nVar].copy()
+
+ # Set up turbulent portion boundary condition
+ avTurb = self.transSolver.solveTransition(var, iMarker)
+
+ # Compute turbulent solution
+ var.flowRegime[iMarker] = 1
+ iMarkerPrev = 0 if iMarker == var.bNodes[1] else iMarker - 1
+ # Avoid starting with 0 for the shear stress.
+ var.u[iMarker*var.nVar + 4] = var.u[iMarkerPrev*var.nVar + 4]
+
+ self.newtonSolver(var, iMarker, displayState)
+
+ # Average solutions
+ avLam = 1 - avTurb
+ thetaTrans = avLam * lamSol[0] + avTurb * var.u[(iMarker)*var.nVar + 0]
+ HTrans = avLam * lamSol[1] + avTurb * var.u[(iMarker)*var.nVar + 1]
+ nTrans = 0
+ ueTrans = avLam * lamSol[3] + avTurb * var.u[(iMarker)*var.nVar + 3]
+ CtTrans = 0 * lamSol[4] + avTurb * var.u[(iMarker)*var.nVar + 4]
+ CtTrans = max(CtTrans, 0.03)
+ var.u[(iMarker)*var.nVar : (iMarker+1)*var.nVar] = [thetaTrans, HTrans, nTrans, ueTrans, CtTrans]
+
+ # Recompute closures
+ var.blPar[(iMarker)*var.nBlPar+1:(iMarker)*var.nBlPar+9]=Closures.turbulentClosure(var.u[(iMarker)*var.nVar+0],
+ var.u[(iMarker)*var.nVar+1],
+ var.u[(iMarker)*var.nVar+4],
+ var.blPar[(iMarker)*var.nBlPar+0],
+ var.extPar[(iMarker)*var.nExtPar+0],
+ 'airfoil')
+ var.blPar[(iMarker-1)*var.nBlPar+1:(iMarker-1)*var.nBlPar+7]=Closures.laminarClosure(var.u[(iMarker-1)*var.nVar+0],
+ var.u[(iMarker-1)*var.nVar+1],
+ var.blPar[(iMarker-1)*var.nBlPar+0],
+ var.extPar[(iMarker-1)*var.nExtPar+0],
+ 'airfoil')
+ pass
+
+ def __forcedTransition(self, var, iMarker):
+ # Forced transition
+ """if not lockTrans and var.xtrF[0] is not None and (var.flowRegime[iMarker] == 1 and var.flowRegime[iMarker-1] == 0):
+ print('| Transition near {:<1.4f}. Marker: {:<5.0f} {:>40}'.format(var.xGlobal[iMarker-1],iMarker-1, '|'))
+ lamSol = var.u[(iMarker-1)*var.nVar : (iMarker)* var.nVar].copy()
+ avTurb = self.transSolver.transitionBC(var, iMarker - 1)
+ # Compute turbulent solution
+ var.flowRegime[iMarker - 1] = 1
+ self.newtonSolver(var, iMarker-1, displayState)
+ # Average solutions
+ var.u[(iMarker-1)*var.nVar : (iMarker)*var.nVar] = (1-avTurb) * lamSol + avTurb * var.u[(iMarker-1)*var.nVar : (iMarker)*var.nVar]
+ # Recompute closures
+ var.blPar[(iMarker-1)*nBlPar+1:(iMarker-1)*nBlPar+9]=Closures.turbulentClosure(var.u[(iMarker-1)*var.nVar+0],
+ var.u[(iMarker-1)*var.nVar+1],
+ var.u[(iMarker-1)*var.nVar+4],
+ var.blPar[(iMarker-1)*nBlPar+0],
+ var.extPar[(iMarker-1)*nExtPar+0],
+ 'airfoil')
+
+ lockTrans = 1
+ if not lockTrans and var.xtrF[1] is not None and (var.flowRegime[iMarker]== 1 and var.flowRegime[iMarker-1] == 0):
+ print('| Transition near {:<1.4f}. Marker: {:<5.0f} {:>40}'.format(var.xGlobal[iMarker-1],iMarker-1, '|'))
+ lamSol = var.u[(iMarker-1)*var.nVar : (iMarker)* var.nVar].copy()
+ avTurb = self.transSolver.transitionBC(var, iMarker - 1)
+ # Compute turbulent solution
+ var.flowRegime[iMarker - 1] = 1
+ self.newtonSolver(var, iMarker-1, displayState)
+ # Average solutions
+ var.u[(iMarker-1)*var.nVar : (iMarker)*var.nVar] = (1-avTurb) * lamSol + avTurb * var.u[(iMarker-1)*var.nVar : (iMarker)*var.nVar]
+ # Recompute closures
+ var.blPar[(iMarker-1)*nBlPar+1:(iMarker-1)*nBlPar+9]=Closures.turbulentClosure(var.u[(iMarker-1)*var.nVar+0],
+ var.u[(iMarker-1)*var.nVar+1],
+ var.u[(iMarker-1)*var.nVar+4],
+ var.blPar[(iMarker-1)*nBlPar+0],
+ var.extPar[(iMarker-1)*nExtPar+0],
+ 'airfoil')"""
+ pass
\ No newline at end of file
diff --git a/dart/viscous/Solvers/VTransition.py b/dart/viscous/Solvers/VTransition.py
index 6897bd8..461f890 100644
--- a/dart/viscous/Solvers/VTransition.py
+++ b/dart/viscous/Solvers/VTransition.py
@@ -63,13 +63,13 @@ class Transition:
avTurb = (dataCont.xGlobal[transMarker]-dataCont.xtr[0])/(dataCont.xGlobal[dataCont.transMarkers[0]]-dataCont.xGlobal[dataCont.transMarkers[0]-1]) # percentage of the subsection corresponding to a turbulent flow
# Fully turbulent flow on lower side
- elif transMarker==dataCont.bNodes[1]:
+ elif transMarker == dataCont.bNodes[1]:
dataCont.transMarkers[1] = dataCont.bNodes[1]
dataCont.xtr[1] = 0
dataCont.flowRegime[dataCont.bNodes[1] : dataCont.bNodes[2]] = 1
# Lower side
- elif dataCont.bNodes[1]<=transMarker<dataCont.bNodes[2]:
+ elif dataCont.bNodes[1] <= transMarker < dataCont.bNodes[2]:
# Set regime to turbulent on remaining lower side points
dataCont.transMarkers[1] = transMarker
--
GitLab