From 158262ad06e68db89ceaa58f36add40f8ffbffd8 Mon Sep 17 00:00:00 2001
From: Dechamps Paul <paul.dechamps@student.uliege.be>
Date: Fri, 3 Dec 2021 14:05:00 +0100
Subject: [PATCH] Test local
---
dart/tests/bliHighLift.py | 37 +-
dart/tests/bliLowRe.py | 73 ++--
dart/tests/bliNonLift.py | 10 +-
dart/tests/bliSeparated.py | 9 +-
dart/tests/bliTransonic.py | 8 +-
dart/viscous/Drivers/DDriver.py | 42 +-
dart/viscous/Drivers/DGroupConstructor.py | 12 +-
dart/viscous/Drivers/DPostProcess.py | 30 +-
dart/viscous/Master/DCoupler.py | 2 +-
dart/viscous/Physics/DDiscretization.py | 9 +-
dart/viscous/Solvers/DIntegration.py | 79 +---
dart/viscous/Solvers/DSolver.py | 19 +-
dart/viscous/Solvers/DSysMatrix.py | 502 +++++++++++-----------
13 files changed, 372 insertions(+), 460 deletions(-)
diff --git a/dart/tests/bliHighLift.py b/dart/tests/bliHighLift.py
index 900259c..2b8cf1e 100755
--- a/dart/tests/bliHighLift.py
+++ b/dart/tests/bliHighLift.py
@@ -37,13 +37,13 @@
import dart.utils as floU
import dart.default as floD
-import dart.viscous.Master.VCoupler as floVC
-import dart.viscous.Drivers.VDriver as floVSMaster
+import dart.viscous.Master.DCoupler as floVC
+import dart.viscous.Drivers.DDriver as floVSMaster
import dart.viscous.Solvers.Steady.StdCoupler as floVC_Std
import dart.viscous.Solvers.Steady.StdSolver as floVS_Std
-import dart.viscous.Physics.VValidation as validation
+import dart.viscous.Physics.DValidation as validation
import tbox
import tbox.utils as tboxU
@@ -59,36 +59,33 @@ def main():
tms['total'].start()
# define flow variables
- Re = 6e6
+ Re = 1e7
alpha = 12.*math.pi/180
- M_inf = 0
- #user_xtr=[0.0106,1]
- user_xtr=[None,None]
+ M_inf = 0.
+ #user_xtr=[0.012,1]
+ user_xtr=[None, 1]
user_Ti=None
mapMesh = 1
+ nFactor = 2
# Time solver parameters
Time_Domain=True # True for unsteady solver, False for steady solver
- Cfl = 1
- spaceMarching=True
-
- if Time_Domain is True:
- timeParams={'CFL0': Cfl, 'spaceMarching':spaceMarching}
+ CFL0 = 1
# ========================================================================================
U_inf = [math.cos(alpha), math.sin(alpha)] # norm should be 1
c_ref = 1
dim = 2
- tol = 1e-1 # tolerance for the VII (usually one drag count)
- case='Case12Xfoil.dat' # File containing results from XFOIL of the considered case
+ tol = 1e-4 # tolerance for the VII (usually one drag count)
+ case='case15Xfoil.dat' # File containing results from XFOIL of the considered case
+
# 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.002, 'msLe' : 0.0001}
+ 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()
@@ -98,15 +95,14 @@ def main():
tms['pre'].stop()
# solve viscous problem
-
print(ccolors.ANSI_BLUE + 'PySolving...' + ccolors.ANSI_RESET)
tms['solver'].start()
isolver = floD.newton(pbl)
- # Choose between steady and unsteady solver
+ # Choose between steady and unsteady solver
if Time_Domain is True: # Run unsteady solver
- vsolver=floVSMaster.Driver(Re, user_xtr, timeParams, M_inf, case, mapMesh)
- coupler = floVC.Coupler(isolver, vsolver, blw[0], blw[1], tol, gmshWriter, bnd, timeParams)
+ vsolver=floVSMaster.Driver(Re, user_xtr, CFL0, M_inf, case, mapMesh, nFactor)
+ coupler = floVC.Coupler(isolver, vsolver, blw[0], blw[1], tol, gmshWriter, bnd)
elif Time_Domain is False: # Run steady solver
vsolver=floVS_Std.Solver(Re)
@@ -134,7 +130,6 @@ def main():
val=validation.Validation(case)
val.main(isolver.Cl,vsolver.wData)
-
# check results
print(ccolors.ANSI_BLUE + 'PyTesting...' + ccolors.ANSI_RESET)
tests = CTests()
diff --git a/dart/tests/bliLowRe.py b/dart/tests/bliLowRe.py
index 8f0860c..10d9536 100755
--- a/dart/tests/bliLowRe.py
+++ b/dart/tests/bliLowRe.py
@@ -37,13 +37,13 @@
import dart.utils as floU
import dart.default as floD
-import dart.viscous.Master.VCoupler as floVC
-import dart.viscous.Drivers.VDriver as floVSMaster
+import dart.viscous.Master.DCoupler as floVC
+import dart.viscous.Drivers.DDriver as floVSMaster
import dart.viscous.Solvers.Steady.StdCoupler as floVC_Std
import dart.viscous.Solvers.Steady.StdSolver as floVS_Std
-import dart.viscous.Physics.VValidation as validation
+import dart.viscous.Physics.DValidation as validation
import tbox
import tbox.utils as tboxU
@@ -57,48 +57,37 @@ def main():
# timer
tms = fwk.Timers()
tms['total'].start()
-
- # Time solver parameters
- Time_Domain=True # True for unsteady solver, False for steady solver
- model='implicitEuler' # 'eulerExp','RungeKutta', 'AdamsBashworth', 'MostAccurateExplicit', 'Leapfrog' or 'implicitEuler'
- #user_xtr=[0.061, 0.740] # User forced transition : Use None type for free transition on either side
+
+ # define flow variables
+ Re = 5e5
+ alpha = 2.*math.pi/180
+ M_inf = 0.2
+
+ #user_xtr=[0.41,0.41]
+ #user_xtr=[0,None]
+ user_xtr=[None, None]
user_Ti=None
+ mapMesh = 0
+ nFactor = 2
+
+ # Time solver parameters
+ Time_Domain = True # True for unsteady solver, False for steady solver
+ CFL0 = 1
+
+ # ========================================================================================
- Re = 50e3
- alpha=5*math.pi/180
- user_xtr=[0.45,1]
- M_inf=0.4
- Cfl=1
U_inf = [math.cos(alpha), math.sin(alpha)] # norm should be 1
c_ref = 1
dim = 2
tol = 1e-4 # tolerance for the VII (usually one drag count)
-
-
-
- # Butcher's table coefficients
- #RK_a=[[1],[-10890423/25193600, 5e3/7873],[-102261/5e6,-5121/2e4,7873/1e4]]
- #RK_b=[1/10,1/6,0,1/6]
- case='m0a5.txt' # File containing results from XFOIL of the considered case
-
-
- try: user_xtr
- except NameError: user_xtr = [None,None]
-
- try: user_Ti
- except NameError: user_Ti=None
-
- try: RK_a
- except NameError: RK_a=None
- try: RK_b
- except NameError: RK_b=None
- timeParams={'modelName' : model, 'CFL': Cfl, 'Rka': RK_a, 'Rkb': RK_b}
+ case='Case1FreeTrans.dat' # File containing results from XFOIL of the considered case
# 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.001}
+ #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()
@@ -108,19 +97,22 @@ def main():
tms['pre'].stop()
# solve viscous problem
+
print(ccolors.ANSI_BLUE + 'PySolving...' + ccolors.ANSI_RESET)
tms['solver'].start()
isolver = floD.newton(pbl)
- # Choose between steady and unsteady solver
- if Time_Domain is True:
- vsolver=floVSMaster.Driver(Re, user_xtr, timeParams, case)
- coupler = floVC.Coupler(isolver, vsolver, blw[0], blw[1], tol, gmshWriter, bnd, timeParams)
- elif Time_Domain is False:
+ # Choose between steady and unsteady solver
+ if Time_Domain is True: # Run unsteady solver
+ vsolver=floVSMaster.Driver(Re, user_xtr, CFL0, M_inf, case, mapMesh, nFactor)
+ coupler = floVC.Coupler(isolver, vsolver, blw[0], blw[1], tol, gmshWriter, bnd)
+
+ elif Time_Domain is False: # Run steady solver
vsolver=floVS_Std.Solver(Re)
coupler = floVC_Std.Coupler(isolver, vsolver, blw[0], blw[1], tol, gmshWriter)
coupler.run()
tms['solver'].stop()
+
# extract Cp
Cp = floU.extract(bnd.groups[0].tag.elems, isolver.Cp)
tboxU.write(Cp, 'Cp_airfoil.dat', '%1.5e', ', ', 'x, y, z, Cp', '')
@@ -139,6 +131,9 @@ def main():
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.Cd, isolver.Cm, 4), True)
+ val=validation.Validation(case)
+ val.main(isolver.Cl,vsolver.wData)
+
# check results
print(ccolors.ANSI_BLUE + 'PyTesting...' + ccolors.ANSI_RESET)
tests = CTests()
diff --git a/dart/tests/bliNonLift.py b/dart/tests/bliNonLift.py
index 8ac5aef..996c967 100755
--- a/dart/tests/bliNonLift.py
+++ b/dart/tests/bliNonLift.py
@@ -59,21 +59,21 @@ def main():
tms['total'].start()
# define flow variables
- Re = 6e6
+ Re = 1e7
alpha = 5.*math.pi/180
- M_inf = 0.15
+ M_inf = 0.2
#user_xtr=[0.41,0.41]
#user_xtr=[0,None]
user_xtr=[None,None]
user_Ti=None
- mapMesh = 1
- nFactor = 3
+ mapMesh = 0
+ nFactor = 2
# Time solver parameters
Time_Domain = True # True for unsteady solver, False for steady solver
- CFL0 = 3
+ CFL0 = 2
# ========================================================================================
diff --git a/dart/tests/bliSeparated.py b/dart/tests/bliSeparated.py
index 47d4fd1..6593b48 100755
--- a/dart/tests/bliSeparated.py
+++ b/dart/tests/bliSeparated.py
@@ -59,15 +59,14 @@ def main():
tms['total'].start()
# define flow variables
- Re = 1e6
+ Re = 1e7
alpha = 15.*math.pi/180
- M_inf = 0.2
- #user_xtr=[0.012,1]
+ M_inf = 0
user_xtr=[None,None]
user_Ti=None
mapMesh = 1
- nFactor = 7
+ nFactor = 2
# Time solver parameters
Time_Domain=True # True for unsteady solver, False for steady solver
@@ -85,7 +84,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.0001}
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 c923b03..673a825 100755
--- a/dart/tests/bliTransonic.py
+++ b/dart/tests/bliTransonic.py
@@ -65,14 +65,14 @@ def main():
#user_xtr=[0.41,0.41]
#user_xtr=[0,None]
- user_xtr=[None,None]
+ user_xtr=[None, None]
user_Ti=None
mapMesh = 1
- nFactor = 5
+ nFactor = 2
# Time solver parameters
- Time_Domain=True # True for unsteady solver, False for steady solver
+ Time_Domain = True # True for unsteady solver, False for steady solver
CFL0 = 1
# ========================================================================================
@@ -80,7 +80,7 @@ def main():
U_inf = [math.cos(alpha), math.sin(alpha)] # norm should be 1
c_ref = 1
dim = 2
- tol = 1e-3 # tolerance for the VII (usually one drag count)
+ tol = 1e-4 # tolerance for the VII (usually one drag count)
case='Case1FreeTrans.dat' # File containing results from XFOIL of the considered case
# mesh the geometry
diff --git a/dart/viscous/Drivers/DDriver.py b/dart/viscous/Drivers/DDriver.py
index c2bcf7d..5f5b2fa 100644
--- a/dart/viscous/Drivers/DDriver.py
+++ b/dart/viscous/Drivers/DDriver.py
@@ -37,15 +37,15 @@ class Driver:
"""
def __init__(self,_Re, user_xtr, _CFL0, _Minfty, _case, _mapMesh, _nFactor):
- self.Minfty = _Minfty
self.Re=_Re
+ self.Minfty = _Minfty
+ self.xtrF=user_xtr
+ self.Ti=None
self.Cd=0
self.it=0
self.uPrevious=None
- self.xtrF=user_xtr
self.mapMesh = _mapMesh
self.nFactor = _nFactor
- self.Ti=None
self.CFL0 = _CFL0
self.res=None
self.upper=0
@@ -138,19 +138,19 @@ class Driver:
dataIn = [None,None] # Merge upper, lower sides and wake
for n in range(0, len(groups)):
- groups[n].connectListNodes, groups[n].connectListElems, dataIn[n] = groups[n].connectList()
+ groups[n].connectListNodes, groups[n].connectListElems, dataIn[n] = groups[n].connectList()
fMeshDict, cMeshDict, data = GroupConstructor().mergeVectors(dataIn, self.mapMesh, nFactor)
var = Variables(fMeshDict, self.xtrF, self.Ti, self.Re) # Initialize data container for flow config.
- if self.it!=0:
- var.extPar[4::var.nExtPar] = fMeshDict['xxExt'] # Extenal flow local markers.
+ if self.it!=0 and self.nbrElmUpper == var.bNodes[1]:
+ var.extPar[4::var.nExtPar] = fMeshDict['xxExt'] # Extenal flow local markers.
DirichletBC=Conditions.Boundaries(self.Re) # Boundary condition (Airfoil/Wake).
gradComp=Discretization(var.xx, var.extPar[4::var.nExtPar], var.bNodes) # Gradients computation.
- sysSolver = sysMatrix(var.nN, var.nVar, gradComp.fou, gradComp.fou) # Boundary layer equations solver that can provide
+ sysSolver = sysMatrix(var.nN, var.nVar, gradComp.fou,gradComp.fou) # Boundary layer equations solver that can provide
# Residuals and Jacobian computation.
transSolver = Transition(var.xtrF) # Transition solver.
@@ -165,7 +165,7 @@ 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(cMeshDict['locCoord']), var.nN, '|'))
+ print('| - Mesh mapped from {:>5.0f} to {:>5.0f} points.{:>37s}'.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('|'))
@@ -190,24 +190,25 @@ class Driver:
print('| - Maximum Mach number: {:<.2f}. {:>49s}'.format(np.max((var.extPar[0::var.nExtPar])),'|'))
- # Initial condition.
+ # Solver setup.
- if self.uPrevious is None: # Typically for the first iteration.
+ if self.uPrevious is None or self.nbrElmUpper != var.bNodes[1]: # Typically for the first iteration.
tSolver.initSln = 1 # Flag to use time solver initial condition.
- tSolver.integrator.jacEval0 = 1 # Continuous Jacobian evaluation.
+ #tSolver.integrator.itFrozenJac0 = 1 # Continuous Jacobian evaluation.
tSolver.integrator.SetCFL0(0.5) # Low starting CFL.
-
+ if self.it != 0 and self.nbrElmUpper != var.bNodes[1]:
+ print('| - Stagnation point moved. {:>29}'.format('|'))
print('| - Time solver will provide the initial solution. {:>29}'.format('|'))
- else:
- # If we use a solution previously obtained.
-
+ else: # If we use a solution previously obtained.
+
var.u = self.uPrevious.copy() # Use previous solution.
print('| - Using previous solution as initial guess. {:>34}'.format('|'))
+ self.nbrElmUpper = var.bNodes[1]
DirichletBC.airfoilBC(var) # Reset boundary condition.
# Run the time integration
@@ -258,19 +259,10 @@ class Driver:
self.Cd = self.blScalars[0]
self.Cdf = self.blScalars[1]
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 = 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 = blwVelWk[groups[1].connectListElems.argsort()]
-
- """print('Marker negative Ct', np.where(var.u[4::var.nVar]<0))
- if self.it % 5 == 0:
- plt.plot(var.u[4::var.nVar])
- plt.axvline(x=var.transMarkers[0],color='r')
- plt.show()
- self.writeFile()
- Validation().main(1,self.wData, WKblVectors, var.xGlobal[var.bNodes[2]:])"""
\ No newline at end of file
+ groups[1].u = blwVelWk[groups[1].connectListElems.argsort()]
\ No newline at end of file
diff --git a/dart/viscous/Drivers/DGroupConstructor.py b/dart/viscous/Drivers/DGroupConstructor.py
index c03b482..281a60a 100755
--- a/dart/viscous/Drivers/DGroupConstructor.py
+++ b/dart/viscous/Drivers/DGroupConstructor.py
@@ -136,15 +136,6 @@ class GroupConstructor():
data = np.concatenate((dataUpper, dataLower, dataWake), axis = 0)
- 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')
- 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 fMeshDict, cMeshDict, data
def __getBLcoordinates(self,data):
@@ -195,5 +186,4 @@ class GroupConstructor():
fMesh = np.append(fMesh, cMesh[iPoint] + iRefinedPoint * dx / nFactor)
fMesh = np.append(fMesh, cMesh[-1])
- return fMesh
-
+ return fMesh
\ No newline at end of file
diff --git a/dart/viscous/Drivers/DPostProcess.py b/dart/viscous/Drivers/DPostProcess.py
index d4ae02a..4649568 100644
--- a/dart/viscous/Drivers/DPostProcess.py
+++ b/dart/viscous/Drivers/DPostProcess.py
@@ -80,7 +80,7 @@ class PostProcess:
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]))
+ 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
@@ -114,33 +114,33 @@ class PostProcess:
def mapBlowingVelocities(self, blwVelUp, blwVelLw, blwVelWk, fbNodes, cMesh, cMeshbNodes, fMesh, nFactor):
"""Maps blowing velocties from cell to cell using weighted average interpolation.
- Args
- ----
+ Args
+ ----
- - blwVelUp (numpy.ndarray): Blowing velocities on the fine mesh of the upper side.
+ - 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.
+ - 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.
+ - blwVelWk (numpy.ndarray): Blowing velocities on the fine mesh of the wake.
- - fbNodes (numpy.ndarray): Boundary nodes of the fine mesh.
+ - fbNodes (numpy.ndarray): Boundary nodes of the fine mesh.
- - cMesh (numpy.ndarray): Coarse mesh coordinates in the global frame of reference.
+ - cMesh (numpy.ndarray): Coarse mesh coordinates in the global frame of reference.
- - cMeshbNodes (numpy.ndarray): Boundary nodes of the coarse mesh.
+ - cMeshbNodes (numpy.ndarray): Boundary nodes of the coarse mesh.
- - fMesh (numpy.ndarray): Fine mesh coordinates in the global frame of reference.
+ - fMesh (numpy.ndarray): Fine mesh coordinates in the global frame of reference.
- - nFactor (int): Multiplicator underlying mesh adaptation.
+ - nFactor (int): Multiplicator underlying mesh adaptation.
Example
-------
Fine mesh: |-----------|-----------|-----------|----------|
- \ /
+ \ /
(1/2) \ / (1/2)
\ /
- \ /
+ \ /
Coarse mesh: |-----------------------|----------------------|
"""
@@ -193,10 +193,10 @@ class PostProcess:
-------
Fine mesh: |-----------|-----------|-----------|----------|
- \ | /
+ \ | /
(1/4) \ | (1/2) / (1/4)
\ V /
- -----> <-----
+ -----> <-----
Coarse mesh: |-----------------------|----------------------|
"""
diff --git a/dart/viscous/Master/DCoupler.py b/dart/viscous/Master/DCoupler.py
index da27a48..c2c0e83 100755
--- a/dart/viscous/Master/DCoupler.py
+++ b/dart/viscous/Master/DCoupler.py
@@ -73,7 +73,7 @@ class Coupler:
if it==0:
# Write inviscid pressure distribution on the first iteration
Cp = floU.extract(self.bnd.groups[0].tag.elems, self.isolver.Cp)
- tboxU.write(Cp, 'CpInv_airfoil.dat', '%1.5e', ', ', 'x, y, z, Cp', '')
+ tboxU.write(Cp, 'Cpinv_airfoil.dat', '%1.5e', ', ', 'x, y, z, Cp', '')
# Run viscous solver
print('+------------------------------ Viscous solver --------------------------------+')
diff --git a/dart/viscous/Physics/DDiscretization.py b/dart/viscous/Physics/DDiscretization.py
index 979774b..85272a1 100755
--- a/dart/viscous/Physics/DDiscretization.py
+++ b/dart/viscous/Physics/DDiscretization.py
@@ -65,13 +65,12 @@ class Discretization:
dueExt_dx=(var.extPar[idx*var.nExtPar+2]-var.extPar[idxPrev*var.nExtPar+2])/dxExt
ddeltaStarExt_dx=(var.extPar[idx*var.nExtPar+3]-var.extPar[idxPrev*var.nExtPar+3])/dxExt
c=4/(np.pi*dx)
- cExt=4/(np.pi*dx)
+ cExt=4/(np.pi*dxExt)
return du_dx, dHstar_dx, dueExt_dx, ddeltaStarExt_dx, c, cExt
def center(self, var, x, idxPrev2, idxPrev,idx, idxNext):
"""Second order central difference.
"""
- #print(idx, idxPrev, idxNext)
du= (var.u[idxNext*var.nVar:(idxNext+1)*var.nVar] - var.u[idxPrev*var.nVar:(idxPrev+1)*var.nVar])
dx = var.xx[idx] - var.xx[idxPrev]
@@ -85,12 +84,6 @@ class Discretization:
cExt=4/(np.pi*dxExt)
return du_dx, dHstar_dx, dueExt_dx, ddeltaStarExt_dx, c, cExt
- '''def fod(self, var, x,idxPrev2, idxPrev,idx, idxNext):
- du=self.u[+var.nVar*idx:+var.nVar*(idx+1)]-x
- du_dx=du/self.dx[idx]
- dHstar_dx=(self.Hstar[idxNext]-self.Hstar[idx])/self.dx[idx]
- return du_dx, dHstar_dx'''
-
# Second order backward difference
def sou(self, var, x,idxPrev2, idxPrev,idx, idxNext): # Attention if we want gradients at 2nd point after stagnation point on lower side.
"""Second order upwing difference.
diff --git a/dart/viscous/Solvers/DIntegration.py b/dart/viscous/Solvers/DIntegration.py
index f1d70c8..ca63840 100644
--- a/dart/viscous/Solvers/DIntegration.py
+++ b/dart/viscous/Solvers/DIntegration.py
@@ -39,11 +39,11 @@ class Integration:
self.sysMatrix = _sysMatrix
self.__CFL0=_Cfl0
- self.__maxIt = 1000
- self.__tol = 1e-4
+ self.__maxIt = 10000
+ self.__tol = 1e-6
self.itFrozenJac0 = 5
- self.__Minf = _Minf
+ self.__Minf = max(_Minf, 0.1)
self.gamma = 1.4
pass
@@ -75,7 +75,7 @@ class Integration:
----
- screenOutput (bool, optional): Output pseudo-time marching convergence.
"""
-
+
# Save initial condition.
sln0 = DFlow.u.copy()
@@ -83,7 +83,6 @@ class Integration:
# Retreive parameters.
nVar = DFlow.nVar # Number of variables of the differential problem.
- nBlPar = DFlow.nBlPar # Number of boundary layer parameters used for vector indexing.
nExtPar = DFlow.nExtPar # Number of boundary layer parameters used for vector indexing.
bNodes = DFlow.bNodes # Boundary nodes of the computational domain.
iInvVel = DFlow.extPar[iMarker*nExtPar + 2] # Inviscid velocity on the considered cell.
@@ -98,7 +97,6 @@ class Integration:
# Initialize pseudo-time integration parameters.
CFL = self.__CFL0 # Initialized CFL number.
- CFLAdapt = 1 # Flag for CFL adaptation.
itFrozenJac = self.itFrozenJac0 # Number of iterations for which Jacobian is frozen.
timeStep = self.__SetTimeStep(CFL, dx, iInvVel) # Initial time step.
IMatrix = np.identity(nVar) / timeStep # Identity matrix used to damp the system.
@@ -129,88 +127,53 @@ class Integration:
sysRes = self.sysMatrix.buildResiduals(DFlow, iMarker)
normSysRes = np.linalg.norm(slnIncr/timeStep - sysRes)
- if screenOutput and innerIters > 0 and innerIters % 1000 == 0:
- self.__OutputState(DFlow, iMarker, innerIters, normSysRes, normSysRes0, CFL, 'yellow')
-
# Check convergence.
if normSysRes <= self.__tol * normSysRes0:
- if screenOutput:
- self.__OutputState(DFlow, iMarker, innerIters, normSysRes, normSysRes0, CFL)
+ # Convergence criterion satisfied.
return 0
except KeyboardInterrupt:
print(ccolors.ANSI_RED + 'Program terminated by user.', ccolors.ANSI_RESET)
quit()
- except Exception as error:
+ except Exception:
nErrors += 1
if nErrors >= 5:
- print('|', ccolors.ANSI_RED + 'Too many errors identified at position {:<1.4f}. Marker: {:<4.0f}. {:>17s}'
- .format(DFlow.xGlobal[iMarker], iMarker,'|'), ccolors.ANSI_RESET)
self.__ResetSolution(DFlow, iMarker, sln0)
return -1
-
-
- print('|', ccolors.ANSI_YELLOW + 'Newton solver failed at position {:<.4f}. Marker: {:<.0f} {:>25s}'
- .format(DFlow.xGlobal[iMarker], iMarker,'|'), ccolors.ANSI_RESET)
-
- print(error)
DFlow.u[iMarker*nVar : (iMarker+1)*nVar] = DFlow.u[(iMarker-1)*nVar : (iMarker)*nVar] # Reset solution.
- DFlow.u[iMarker*nVar + 1] = 1.515
- DFlow.u[iMarker*nVar + 4] = 0.03
# 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
+ CFL = min(max(CFL / 2,0.3),1)
timeStep = self.__SetTimeStep(CFL, dx, iInvVel)
IMatrix = np.identity(nVar) / timeStep
sysRes = self.sysMatrix.buildResiduals(DFlow, iMarker)
itFrozenJac = 1
- screenOutput = 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)
- timeStep = self.__SetTimeStep(CFL, dx, iInvVel)
- IMatrix=np.identity(nVar) / timeStep
+ CFL = self.__CFL0* (normSysRes0/normSysRes)**0.7
+
+ CFL = max(CFL, 0.5)
+ timeStep = self.__SetTimeStep(CFL, dx, iInvVel)
+ IMatrix = np.identity(nVar) / timeStep
- innerIters+=1
+ 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(DFlow.xGlobal[iMarker], iMarker,
- np.log10(normSysRes/normSysRes0),'|'),
- ccolors.ANSI_RESET)
self.__ResetSolution(DFlow, iMarker, sln0)
-
- else:
- print('|', ccolors.ANSI_YELLOW+ 'Too many iterations at position {:<.4f}. Marker: {:>5.0f} normLinSysRes = {:<.3f}. {:>30s}'
- .format(DFlow.xGlobal[iMarker], iMarker,
- np.log10(normSysRes/normSysRes0),'|'),
- ccolors.ANSI_RESET)
return innerIters
@@ -248,22 +211,6 @@ class Integration:
# Time step computation
return CFL*dx/advVel
- def __OutputState(self, var, iMarker, innerIters, normLinSysRes, normLinSysRes0, CFL, color='white') -> None:
- 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)
-
def __ResetSolution(self, DFlow, iMarker, sln0) -> None:
nVar = DFlow.nVar
diff --git a/dart/viscous/Solvers/DSolver.py b/dart/viscous/Solvers/DSolver.py
index 39865e3..c0f7465 100644
--- a/dart/viscous/Solvers/DSolver.py
+++ b/dart/viscous/Solvers/DSolver.py
@@ -49,7 +49,7 @@ class Solver:
nMarkers = DFlow.nN # Number of points in the computational domain.
bNodes = DFlow.bNodes # Boundary nodes of the computational domain.
- convergedPoints = np.zeros(DFlow.nN) # Convergence control of each point.
+ convergedPoints = -1*np.ones(DFlow.nN) # Convergence control of each point.
convergedPoints[0] = 0 # Boundary condition.
self.transSolver.initTransition(DFlow) # Initialize the flow regime on each cell.
@@ -62,9 +62,8 @@ class Solver:
if self.initSln:
iMarkerPrev = 0 if iMarker == DFlow.bNodes[1] else iMarker - 1
DFlow.u[iMarker*DFlow.nVar: (iMarker + 1)*DFlow.nVar] = DFlow.u[iMarkerPrev*DFlow.nVar: (iMarkerPrev + 1)* DFlow.nVar]
- if DFlow.flowRegime[iMarker] == 1:
- #DFlow.u[iMarker * DFlow.nVar + 1] = 1.5
- DFlow.u[iMarker * DFlow.nVar + 4] = max(DFlow.u[iMarker * DFlow.nVar + 4], 0.03)
+ if DFlow.flowRegime[iMarker] == 1 and DFlow.u[iMarker * DFlow.nVar + 4] == 0:
+ DFlow.u[iMarker * DFlow.nVar + 4] = 0.03
if iMarker == bNodes[0] + 1: # Upper side start.
print('| - Computing upper side. {:>54}'.format('|'))
@@ -86,6 +85,10 @@ class Solver:
convergedPoints[iMarker] = self.integrator.TimeMarching(DFlow, iMarker, displayState)
+ if convergedPoints[iMarker] != 0:
+ print(ccolors.ANSI_RED + '| Marker {:<4.0f} (x/c={:<1.4f}): PTI terminated with exit code {:<4.0f} {:>17}'
+ .format(iMarker, DFlow.xGlobal[iMarker], convergedPoints[iMarker], '|'), ccolors.ANSI_RESET)
+
# Check for transition.
if not lockTrans:
@@ -99,11 +102,13 @@ class Solver:
lockTrans = 1
# Forced transition.
- elif DFlow.xtrF[0] is not None or DFlow.xtrF[1] is not None:
- if DFlow.flowRegime[iMarker] == 1 and DFlow.flowRegime[iMarker - 1]:
+ elif (DFlow.xtrF[0] is not None and DFlow.xtrF[0] !=0 and DFlow.xtrF[0] != 1) or (DFlow.xtrF[1] is not None and DFlow.xtrF[1] !=0 and DFlow.xtrF[1] != 1):
+ if DFlow.flowRegime[iMarker] == 1 and DFlow.flowRegime[iMarker - 1] == 0:
print('| Forced transition near x/c = {:<2.3f}. Marker: {:<4.0f} {:>28s}'.format(DFlow.xGlobal[iMarker],iMarker, '|'))
- self.__AverageTransition(DFlow, iMarker)
+ self.__AverageTransition(DFlow, iMarker, displayState)
lockTrans = 1
+ """DFlow.extPar[iMarker*DFlow.nExtPar+2] = DFlow.u[iMarker * DFlow.nVar + 3]
+ DFlow.extPar[iMarker*DFlow.nExtPar+3] = DFlow.blPar[iMarker * DFlow.nBlPar + 9]"""
return convergedPoints
diff --git a/dart/viscous/Solvers/DSysMatrix.py b/dart/viscous/Solvers/DSysMatrix.py
index 0e55194..47ccfd1 100644
--- a/dart/viscous/Solvers/DSysMatrix.py
+++ b/dart/viscous/Solvers/DSysMatrix.py
@@ -20,280 +20,276 @@ import numpy as np
class flowEquations:
- def __init__(self, _setGradients, _fouScheme) -> None:
- self.setGradients = _setGradients
- self.fouScheme = _fouScheme
- pass
-
- def _blLaws(self, dataCont, iMarker, x = None) -> np.ndarray:
- """Evaluates the boundary layer laws at one point of the computational domain.
-
- Args:
- ----
- - dataCont (classType): Data container. Contains the geometry, the solution, BL parameters...
-
- - iMarker (integer): Id of the point.
-
- - x (np.ndarray, optional): Prescribed solution, usually a perturbed one for Jacobian computation purposes. Defaults to None.
-
- Returns:
- -------
- - linSysRes (np.ndarray): Residuals of the linear system at the considered point.
- """
+ def __init__(self, _setGradients, _fouScheme) -> None:
+ self.setGradients = _setGradients
+ self.fouScheme = _fouScheme
+ pass
+
+ def _blLaws(self, dataCont, iMarker, x = None) -> np.ndarray:
+ """Evaluates the boundary layer laws at one point of the computational domain.
+
+ Args:
+ ----
+ - dataCont (classType): Data container. Contains the geometry, the solution, BL parameters...
+
+ - iMarker (integer): Id of the point.
+
+ - x (np.ndarray, optional): Prescribed solution, usually a perturbed one for Jacobian computation purposes. Defaults to None.
+
+ Returns:
+ -------
+ - linSysRes (np.ndarray): Residuals of the linear system at the considered point.
+ """
+
+ if iMarker == dataCont.bNodes[2]:
+ return np.zeros(dataCont.nVar)
+
+ timeMatrix = np.empty((dataCont.nVar, dataCont.nVar))
+ spaceMatrix = np.empty(dataCont.nVar)
+
+ if x is None: x = dataCont.u[iMarker*dataCont.nVar : (iMarker+1)*dataCont.nVar]
+
+ iMarkerPrev2, iMarkerPrev, iMarkerNext, name, xtr, xtrF = self.__evalPoint(iMarker, dataCont.bNodes, dataCont.xtr, dataCont.xtrF, dataCont.nN)
+
+ dissipRatio = 0.9 if name == 'wake' else 1 # Wall/wake dissipation length ratio
+
+ # Reynolds number based on momentum thickness theta
+ if dataCont.flowRegime[iMarker] == 0: dataCont.blPar[iMarker*dataCont.nBlPar+0] = max(x[3] * x[0] * dataCont.Re, 1)
+ else: dataCont.blPar[iMarker*dataCont.nBlPar+0] = x[3] * x[0] * dataCont.Re
+
+ dataCont.blPar[iMarker * dataCont.nBlPar + 9] = x[0] * x[1] # deltaStar
+
+ if dataCont.flowRegime[iMarker]==0:
+
+ # Laminar closure relations.
+ dataCont.blPar[iMarker*dataCont.nBlPar+1:iMarker*dataCont.nBlPar+7] = Closures.laminarClosure(x[0], x[1],
+ dataCont.blPar[iMarker*dataCont.nBlPar+0],
+ dataCont.extPar[iMarker*dataCont.nExtPar+0], name)
+
+ elif dataCont.flowRegime[iMarker]==1:
+
+ # Turbulent closures relations.
+ dataCont.blPar[iMarker*dataCont.nBlPar+1:iMarker*dataCont.nBlPar+9] = Closures.turbulentClosure(x[0], x[1],
+ x[4], dataCont.blPar[iMarker*dataCont.nBlPar+0],
+ dataCont.extPar[iMarker*dataCont.nExtPar+0], name)
+
+ # Rename variables for clarity.
+ Ta, Ha, _, uea, Cta = x
+ Cda, Cfa, deltaa, Hstara, Hstar2a, Hka, Cteqa, Usa = dataCont.blPar[iMarker * dataCont.nBlPar + 1 : iMarker * dataCont.nBlPar + 9]
+
+ # Amplification factor for e^n method.
+ if name=='airfoil' and xtrF is None and dataCont.flowRegime[iMarker]==0:
+
+ Axa = Closures.amplRoutine2(dataCont.blPar[iMarker * dataCont.nBlPar + 6], dataCont.blPar[iMarker * dataCont.nBlPar + 0], Ta)
+
+ else:
+
+ Axa=0
+ dN_dx=0
+
+ if iMarker == 1 or iMarker == dataCont.bNodes[1] or iMarker == dataCont.bNodes[2]+1:
+ [dTheta_dx, dH_dx, dN_dx, due_dx, dCt_dx], dHstar_dx, dueExt_dx, ddeltaStarExt_dx, c, cExt = self.fouScheme(dataCont, x,
+ iMarkerPrev2, iMarkerPrev,
+ iMarker, iMarkerNext)
+ else:
+ [dTheta_dx, dH_dx, dN_dx, due_dx, dCt_dx], dHstar_dx, dueExt_dx, ddeltaStarExt_dx, c, cExt = self.setGradients(dataCont, x,
+ iMarkerPrev2, iMarkerPrev,
+ iMarker, iMarkerNext)
+
+ # +------------------------------------------- Equations -----------------------------------------------------------+
+ # | |
+ # | |
+ # | Von Karman momentum integral equation 0-th (momentum integral) |
+ # | -------------------------------------------------------------- |
+ # | |
+ # | H/ue*dTheta/dt + Theta/ue*dH/dt + Theta*H/(ue^2)*due/dt + dTheta/dx + (2+H)*(Theta/ue)*due/dx - Cf/2 = 0 |
+ # | |
+ # | 1-st (kinetic-energy integral) moments of momentum equation |
+ # | ----------------------------------------------------------- |
+ # | |
+ # | (1+H*(1-Hstar))/ue*dTheta/dt + (1-Hstar)*Theta/ue*dH/dt + (2-Hstar*H)*Theta/(ue^2)*due/dt |
+ # | + Theta*(dHstar/dx) + (2*Hstar2 + Hstar*(1-H))*Theta/ue*due_dx - 2*Cd + Hstar*Cf/2 = 0 |
+ # | |
+ # | Unsteady e^n method |
+ # | ------------------- |
+ # | |
+ # | dN/dt + dN/dx - Ax = 0 |
+ # | |
+ # | Interaction law |
+ # | --------------- |
+ # | |
+ # | due/dt - c*H*dTheta/dt - c*Theta*dH/dt + due/dx-c*(H*dTheta/dx+Theta*dH/dx) + (-due/dx + c*deltaStar)_Old = 0 |
+ # | |
+ # | Unsteady shear lag equation (ignored in the laminar case) |
+ # | --------------------------------------------------------- |
+ # | |
+ # | 2*delta/(Us*ue^2)*due/dt + 2*delta/(ue*Ct*Us)*dCt/dt + (2*delta/Ct)*(dCt/dx) |
+ # | - 5.6*(Cteq-Ct*dissipRatio)-2*delta*(4/(3*H*Theta)*(Cf/2-((Hk-1)/(6.7*Hk*dissipRatio))^2)-1/ue*due/dx) = 0 |
+ # | |
+ # +-----------------------------------------------------------------------------------------------------------------+
+
+
+ # Spacial part
+ spaceMatrix[0] = dTheta_dx + (2+Ha) * (Ta/uea) * due_dx - Cfa/2
+ spaceMatrix[1] = Ta * dHstar_dx + ( 2*Hstar2a + Hstara * (1-Ha)) * Ta/uea * due_dx - 2*Cda + Hstara * Cfa/2
+ if xtrF is not None or dataCont.flowRegime[iMarker] == 1:
+ spaceMatrix[2] = 0
+ else:
+ spaceMatrix[2] = dN_dx - Axa
+ spaceMatrix[3] = due_dx -c*(Ha * dTheta_dx + Ta * dH_dx) - dueExt_dx + cExt * ddeltaStarExt_dx
+ if dataCont.flowRegime[iMarker] == 1:
+
+ spaceMatrix[4] = (2*deltaa/Cta) * (dCt_dx) - 5.6*(Cteqa - Cta * dissipRatio)
+ - 2*deltaa*(4/(3*Ha * Ta) * (Cfa/2 - ((Hka-1)/(6.7*Hka * dissipRatio))**2) - 1/uea * due_dx)
+
+ elif dataCont.flowRegime[iMarker] == 0: spaceMatrix[4] = 0
+
+ # Temporal part
+ timeMatrix[0] = [Ha/uea, Ta/uea, 0, Ta * Ha/(uea**2), 0]
+ timeMatrix[1] = [(1+Ha * (1-Hstara))/uea, (1-Hstara) * Ta/uea, 0, (2-Hstara * Ha)*Ta/(uea**2), 0]
+ timeMatrix[2] = [0, 0, 1, 0, 0]
+ timeMatrix[3] = [-c*Ha, -c*Ta, 0, 1, 0]
+ if dataCont.flowRegime[iMarker]==1:
+ timeMatrix[4] = [0, 0, 0, 2*deltaa/(Usa * uea**2), 2*deltaa/(uea * Cta * Usa)]
+ else:
+ # Shear lag equation ignored in the laminar portion of the flow
+ timeMatrix[4] = [0, 0, 0, 0, 1]
- if iMarker == dataCont.bNodes[2]:
- return np.zeros(dataCont.nVar)
-
- timeMatrix = np.empty((dataCont.nVar, dataCont.nVar))
- spaceMatrix = np.empty(dataCont.nVar)
-
- if x is None: x = dataCont.u[iMarker*dataCont.nVar : (iMarker+1)*dataCont.nVar]
-
- iMarkerPrev2, iMarkerPrev, iMarkerNext, name, xtr, xtrF = self.__evalPoint(iMarker, dataCont.bNodes, dataCont.xtr, dataCont.xtrF, dataCont.nN)
-
- dissipRatio = 0.9 if name == 'wake' else 1 # Wall/wake dissipation length ratio
-
- # Reynolds number based on momentum thickness theta
- if dataCont.flowRegime[iMarker] == 0: dataCont.blPar[iMarker*dataCont.nBlPar+0] = max(x[3] * x[0] * dataCont.Re, 1)
- else: dataCont.blPar[iMarker*dataCont.nBlPar+0] = x[3] * x[0] * dataCont.Re
-
- dataCont.blPar[iMarker * dataCont.nBlPar + 9] = x[0] * x[1] # deltaStar
-
- if dataCont.flowRegime[iMarker]==0:
-
- # Laminar closure relations
- dataCont.blPar[iMarker*dataCont.nBlPar+1:iMarker*dataCont.nBlPar+7] = Closures.laminarClosure(x[0], x[1],
- dataCont.blPar[iMarker*dataCont.nBlPar+0],
- dataCont.extPar[iMarker*dataCont.nExtPar+0], name)
-
- elif dataCont.flowRegime[iMarker]==1:
-
- # Turbulent closures relations
- dataCont.blPar[iMarker*dataCont.nBlPar+1:iMarker*dataCont.nBlPar+9] = Closures.turbulentClosure(x[0], x[1],
- x[4], dataCont.blPar[iMarker*dataCont.nBlPar+0],
- dataCont.extPar[iMarker*dataCont.nExtPar+0], name)
-
- # Rename variables for clarity
- Ta, Ha, _, uea, Cta = x
- Cda, Cfa, deltaa, Hstara, Hstar2a, Hka, Cteqa, Usa = dataCont.blPar[iMarker * dataCont.nBlPar + 1 : iMarker * dataCont.nBlPar + 9]
-
- # Amplification factor for e^n method
- if name=='airfoil' and xtrF is None and dataCont.flowRegime[iMarker]==0:
-
- Axa = Closures.amplRoutine2(dataCont.blPar[iMarker * dataCont.nBlPar + 6], dataCont.blPar[iMarker * dataCont.nBlPar + 0], Ta)
-
- else:
-
- Axa=0
- dN_dx=0
-
- if iMarker == dataCont.bNodes[1] - 1 or iMarker == dataCont.nN-1 or iMarker == dataCont.nN - 1:
- [dTheta_dx, dH_dx, dN_dx, due_dx, dCt_dx], dHstar_dx, dueExt_dx, ddeltaStarExt_dx, c, cExt = self.fouScheme(dataCont, x,
- iMarkerPrev2, iMarkerPrev,
- iMarker, iMarkerNext)
- else:
- [dTheta_dx, dH_dx, dN_dx, due_dx, dCt_dx], dHstar_dx, dueExt_dx, ddeltaStarExt_dx, c, cExt = self.setGradients(dataCont, x,
- iMarkerPrev2, iMarkerPrev,
- iMarker, iMarkerNext)
-
- # +------------------------------------------- Equations -----------------------------------------------------------+
- # | |
- # | |
- # | Von Karman momentum integral equation 0-th (momentum integral) |
- # | -------------------------------------------------------------- |
- # | |
- # | H/ue*dTheta/dt + Theta/ue*dH/dt + Theta*H/(ue^2)*due/dt + dTheta/dx + (2+H)*(Theta/ue)*due/dx - Cf/2 = 0 |
- # | |
- # | 1-st (kinetic-energy integral) moments of momentum equation |
- # | ----------------------------------------------------------- |
- # | |
- # | (1+H*(1-Hstar))/ue*dTheta/dt + (1-Hstar)*Theta/ue*dH/dt + (2-Hstar*H)*Theta/(ue^2)*due/dt |
- # | + Theta*(dHstar/dx) + (2*Hstar2 + Hstar*(1-H))*Theta/ue*due_dx - 2*Cd + Hstar*Cf/2 = 0 |
- # | |
- # | Unsteady e^n method |
- # | ------------------- |
- # | |
- # | dN/dt + dN/dx - Ax = 0 |
- # | |
- # | Interaction law |
- # | --------------- |
- # | |
- # | due/dt - c*H*dTheta/dt - c*Theta*dH/dt + due/dx-c*(H*dTheta/dx+Theta*dH/dx) + (-due/dx + c*deltaStar)_Old = 0 |
- # | |
- # | Unsteady shear lag equation (ignored in the laminar case) |
- # | --------------------------------------------------------- |
- # | |
- # | 2*delta/(Us*ue^2)*due/dt + 2*delta/(ue*Ct*Us)*dCt/dt + (2*delta/Ct)*(dCt/dx) |
- # | - 5.6*(Cteq-Ct*dissipRatio)-2*delta*(4/(3*H*Theta)*(Cf/2-((Hk-1)/(6.7*Hk*dissipRatio))^2)-1/ue*due/dx) = 0 |
- # | |
- # +-----------------------------------------------------------------------------------------------------------------+
-
-
- # Spacial part
- spaceMatrix[0] = dTheta_dx + (2+Ha) * (Ta/uea) * due_dx - Cfa/2
- spaceMatrix[1] = Ta * dHstar_dx + ( 2*Hstar2a + Hstara * (1-Ha)) * Ta/uea * due_dx - 2*Cda + Hstara * Cfa/2
- if xtrF is not None or dataCont.flowRegime[iMarker] == 1:
- spaceMatrix[2] = 0
- else:
- spaceMatrix[2] = dN_dx - Axa
- spaceMatrix[3] = due_dx -c*(Ha * dTheta_dx + Ta * dH_dx) - dueExt_dx + cExt * ddeltaStarExt_dx
- if dataCont.flowRegime[iMarker] == 1:
-
- spaceMatrix[4] = (2*deltaa/Cta) * (dCt_dx) - 5.6*(Cteqa - Cta * dissipRatio)
- - 2*deltaa*(4/(3*Ha * Ta) * (Cfa/2 - ((Hka-1)/(6.7*Hka * dissipRatio))**2) - 1/uea * due_dx)
-
- elif dataCont.flowRegime[iMarker] == 0: spaceMatrix[4] = 0
-
- # Temporal part
- timeMatrix[0] = [Ha/uea, Ta/uea, 0, Ta * Ha/(uea**2), 0]
- timeMatrix[1] = [(1+Ha * (1-Hstara))/uea, (1-Hstara) * Ta/uea, 0, (2-Hstara * Ha)*Ta/(uea**2), 0]
- timeMatrix[2] = [0, 0, 1, 0, 0]
- timeMatrix[3] = [-c*Ha, -c*Ta, 0, 1, 0]
- if dataCont.flowRegime[iMarker]==1:
- timeMatrix[4] = [0, 0, 0, 2*deltaa/(Usa * uea**2), 2*deltaa/(uea * Cta * Usa)]
- else:
- # 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
-
-
- def __evalPoint(self,iMarker, bNodes, xtrIn, xtrFIn, nMarker) -> tuple:
- """Evaluates where the point is, the adjactent points
- and the state of transition on the corresponding side
-
- Args:
- ----
- - iMarker (integer): Id of the cell to evaluate.
-
- - bNodes (np.ndarray): Boundary nodes of the computational domain.
-
- - xtrIn (np.ndarray): Current transiton location on [upper, lower] sides.
-
- - xtrFIn (np.ndarray): Forced transition state on [upper, lower] sides.
-
- - nMarker ([type]): Number of cells in the computational domain.
-
- Returns:
- -------
- - tuple: Id of the adjacent points, name of the corresponding group ('airfoil', 'wake'),
- transition state on the corresponding side.
- """
-
- if iMarker == bNodes[1]: # First point of the lower side (next to the stagnation point)
- iMarkerPrev = 0
- iMarkerPrev2 = None
-
- elif iMarker == bNodes[1]+1:
- iMarkerPrev = iMarker-1
- iMarkerPrev2 = 0
- else:
- iMarkerPrev = iMarker-1
- iMarkerPrev2 = iMarker-2
-
- if iMarker == bNodes[1]-1 or iMarker == bNodes[2]-1 or iMarker == nMarker-1:
- iMarkerNext = None
- else: iMarkerNext = iMarker+1
-
- if iMarker >= bNodes[2]:
- name='wake'
- xtr=0
- xtrF=None
-
- elif iMarker < bNodes[1]:
- name = 'airfoil'
- xtr = xtrIn[0]
- xtrF = xtrFIn[0]
- elif bNodes[1] <= iMarker < bNodes[2]:
- name = 'airfoil'
- xtr = xtrIn[1]
- xtrF = xtrFIn[1]
- return iMarkerPrev2, iMarkerPrev, iMarkerNext, name, xtr, xtrF
+ sysRes = np.linalg.solve(timeMatrix, spaceMatrix).flatten()
+ return sysRes
+
+
+ def __evalPoint(self,iMarker, bNodes, xtrIn, xtrFIn, nMarker) -> tuple:
+ """Evaluates where the point is, the adjactent points
+ and the state of transition on the corresponding side
+
+ Args:
+ ----
+ - iMarker (integer): Id of the cell to evaluate.
+
+ - bNodes (np.ndarray): Boundary nodes of the computational domain.
+
+ - xtrIn (np.ndarray): Current transiton location on [upper, lower] sides.
+
+ - xtrFIn (np.ndarray): Forced transition state on [upper, lower] sides.
+
+ - nMarker ([type]): Number of cells in the computational domain.
+
+ Returns:
+ -------
+ - tuple: Id of the adjacent points, name of the corresponding group ('airfoil', 'wake'),
+ transition state on the corresponding side.
+ """
+ if iMarker == bNodes[0] + 1:
+ iMarkerPrev = 0
+ iMarkerPrev2 = None
+
+ if iMarker == bNodes[1]: # First point of the lower side (next to the stagnation point).
+ iMarkerPrev = 0
+ iMarkerPrev2 = None
+
+ elif iMarker == bNodes[1]+1:
+ iMarkerPrev = iMarker-1
+ iMarkerPrev2 = 0
+ else:
+ iMarkerPrev = iMarker-1
+ iMarkerPrev2 = iMarker-2
+
+ if iMarker == bNodes[1]-1 or iMarker == bNodes[2]-1 or iMarker == nMarker-1:
+ iMarkerNext = None
+ else: iMarkerNext = iMarker+1
+
+ if iMarker >= bNodes[2]:
+ name='wake'
+ xtr=0
+ xtrF=None
+
+ elif iMarker < bNodes[1]:
+ name = 'airfoil'
+ xtr = xtrIn[0]
+ xtrF = xtrFIn[0]
+ elif bNodes[1] <= iMarker < bNodes[2]:
+ name = 'airfoil'
+ xtr = xtrIn[1]
+ xtrF = xtrFIn[1]
+ return iMarkerPrev2, iMarkerPrev, iMarkerNext, name, xtr, xtrF
class sysMatrix(flowEquations):
- def __init__(self, _nMarker, _nMarkers, setGradients, fouScheme) -> None:
- flowEquations.__init__(self, setGradients, fouScheme)
- self._nMarker = _nMarker
- self._nMarkers = _nMarkers
- pass
+ def __init__(self, _nMarker, _nMarkers, setGradients, fouScheme) -> None:
+ flowEquations.__init__(self, setGradients, fouScheme)
+ self._nMarker = _nMarker
+ self._nMarkers = _nMarkers
+ pass
- def buildJacobian(self, dataCont, marker, sysRes, order = 1, eps = 1e-10) -> np.ndarray:
- """Contruct the Jacobian used for solving the liner system.
-
- Args:
- ----
- - dataCont (classType): Data container (solution, BL parameters, invicid flow parameters, geometry...)
+ def buildJacobian(self, dataCont, marker, sysRes, order = 1, eps = 1e-10) -> np.ndarray:
+ """Contruct the Jacobian used for solving the liner system.
+
+ Args:
+ ----
+ - dataCont (classType): Data container (solution, BL parameters, invicid flow parameters, geometry...)
- - marker (integer, optional): Cell id on the computational domain. Default : -1 for whole domain.
+ - marker (integer, optional): Cell id on the computational domain. Default : -1 for whole domain.
- - order (unsigned int, optional): Order of the discretization used to obtain the Jacobian (1st or 2nd order implemented).
- If order == 1, linSysRes has to be specifed
+ - order (unsigned int, optional): Order of the discretization used to obtain the Jacobian (1st or 2nd order implemented).
+ If order == 1, linSysRes has to be specifed
- - linSysRes (np.ndarray, optional): Residuals of the linear system at the current iteration. Used (and mandatory) if order == 1.
+ - linSysRes (np.ndarray, optional): Residuals of the linear system at the current iteration. Used (and mandatory) if order == 1.
- - eps (double, optional): Perturbation on the variables. Default : 10^(-8).
+ - eps (double, optional): Perturbation on the variables. Default : 10^(-8).
- Returns
- -------
+ Returns
+ -------
- - JacMatrix (np.ndarray): Jacobian of the system.
+ - JacMatrix (np.ndarray): Jacobian of the system.
- """
-
- # Build Jacobian for one point
- if order == 1 and sysRes is not None:
-
- JacMatrix = np.zeros((dataCont.nVar, dataCont.nVar))
- xUp = dataCont.u[marker*dataCont.nVar : (marker+1)*dataCont.nVar].copy()
+ """
+
+ # Build Jacobian for one point
+ if order == 1 and sysRes is not None:
+
+ JacMatrix = np.zeros((dataCont.nVar, dataCont.nVar))
+ xUp = dataCont.u[marker*dataCont.nVar : (marker+1)*dataCont.nVar].copy()
- for iVar in range(dataCont.nVar):
- varSave = xUp[iVar]
- xUp[iVar] += eps
+ for iVar in range(dataCont.nVar):
+ varSave = xUp[iVar]
+ xUp[iVar] += eps
- JacMatrix[:,iVar] = (self._blLaws(dataCont, marker, x=xUp) - sysRes) / eps
+ JacMatrix[:,iVar] = (self._blLaws(dataCont, marker, x=xUp) - sysRes) / eps
- xUp[iVar] = varSave
- return JacMatrix
+ xUp[iVar] = varSave
+ return JacMatrix
- elif order == 2:
- JacMatrix = np.zeros((dataCont.nVar, dataCont.nVar))
- xUp = dataCont.u[marker*dataCont.nVar : (marker+1)*dataCont.nVar].copy()
- xDwn = dataCont.u[marker*dataCont.nVar : (marker+1)*dataCont.nVar].copy()
+ elif order == 2:
+ JacMatrix = np.zeros((dataCont.nVar, dataCont.nVar))
+ xUp = dataCont.u[marker*dataCont.nVar : (marker+1)*dataCont.nVar].copy()
+ xDwn = dataCont.u[marker*dataCont.nVar : (marker+1)*dataCont.nVar].copy()
- for iVar in range(dataCont.nVar):
- varSave = xUp[iVar]
- xUp[iVar] += eps
- xDwn[iVar] -= eps
+ for iVar in range(dataCont.nVar):
+ varSave = xUp[iVar]
+ xUp[iVar] += eps
+ xDwn[iVar] -= eps
- JacMatrix[:,iVar] = (self._blLaws(dataCont, marker, x = xUp) - self._blLaws(dataCont, marker, x = xDwn)) / (2*eps)
+ JacMatrix[:,iVar] = (self._blLaws(dataCont, marker, x = xUp) - self._blLaws(dataCont, marker, x = xDwn)) / (2*eps)
- xUp[iVar] = varSave
- xDwn[iVar] = varSave
- return JacMatrix
+ xUp[iVar] = varSave
+ xDwn[iVar] = varSave
+ return JacMatrix
- else:
- raise RuntimeError('Only first and second orders Jacobian can be computed.')
+ else:
+ raise RuntimeError('Only first and second orders Jacobian can be computed.')
- def buildResiduals(self, dataCont, marker) -> np.ndarray:
- """Construct residuals of the linear system.
-
- Args:
- ----
-
- - dataCont (classType): Data container (solution, BL parameters, invicid flow parameters, geometry...)
+ def buildResiduals(self, dataCont, marker) -> np.ndarray:
+ """Construct residuals of the linear system.
+
+ Args:
+ ----
+
+ - dataCont (classType): Data container (solution, BL parameters, invicid flow parameters, geometry...)
- - marker (int, optional): Cell id on the computational domain. default: -1 for whole domain.
- """
+ - marker (int, optional): Cell id on the computational domain. default: -1 for whole domain.
+ """
- # Build Residual for one point
- return self._blLaws(dataCont, marker, x=None)
\ No newline at end of file
+ # Build Residual for one point
+ return self._blLaws(dataCont, marker, x=None)
\ No newline at end of file
--
GitLab