From f6ca79ce1619a4dba93eb58bca7769eb3d29ece5 Mon Sep 17 00:00:00 2001
From: Dechamps Paul <paul.dechamps@student.uliege.be>
Date: Mon, 22 Nov 2021 02:04:46 +0100
Subject: [PATCH] Fixes residuals and jacobian computation. Adds mesh mapping
feature.
---
dart/tests/bliHighLift.py | 2 +-
dart/tests/bliNonLift.py | 2 +-
dart/viscous/Drivers/VDriver.py | 178 +++++-----
dart/viscous/Drivers/VGroupConstructor.py | 172 +++++++---
dart/viscous/Drivers/VPostProcess.py | 100 +++++-
dart/viscous/Physics/VBIConditions.py | 142 +++++---
dart/viscous/Physics/VVariables.py | 44 +--
dart/viscous/README.md | 2 +-
dart/viscous/Solvers/Steady/StdSolver.py | 7 +-
dart/viscous/Solvers/VSolver.py | 100 ++----
dart/viscous/Solvers/VTimeSolver.py | 383 ++++++++++++----------
11 files changed, 645 insertions(+), 487 deletions(-)
diff --git a/dart/tests/bliHighLift.py b/dart/tests/bliHighLift.py
index f7a15ed..c2029fa 100755
--- a/dart/tests/bliHighLift.py
+++ b/dart/tests/bliHighLift.py
@@ -70,7 +70,7 @@ def main():
# Time solver parameters
Time_Domain=True # True for unsteady solver, False for steady solver
- Cfl = 1.5
+ Cfl = 1
spaceMarching=True
if Time_Domain is True:
diff --git a/dart/tests/bliNonLift.py b/dart/tests/bliNonLift.py
index a67cd04..23ac27a 100755
--- a/dart/tests/bliNonLift.py
+++ b/dart/tests/bliNonLift.py
@@ -72,7 +72,7 @@ def main():
# Time solver parameters
Time_Domain = True # True for unsteady solver, False for steady solver
- Cfl = 2
+ Cfl = 3
spaceMarching=True
if Time_Domain is True:
diff --git a/dart/viscous/Drivers/VDriver.py b/dart/viscous/Drivers/VDriver.py
index e0821fe..76c926e 100644
--- a/dart/viscous/Drivers/VDriver.py
+++ b/dart/viscous/Drivers/VDriver.py
@@ -31,17 +31,10 @@ from fwk.coloring import ccolors
import numpy as np
import matplotlib.pyplot as plt
-# ------------------------------------------------------------------------------
-# Driver : Input: - Data with the inviscid solution, geometry and user input for
-# time integration
-#
-# Output: - Corrected data ready to be used for viscous inviscd coupling.
-#
-# Sorts the data, allows the whole computational domain to be treated at the same
-# time (equations are solved on the three groups simultaneously) and computes
-# required data for the inviscid solver (blowing velocity aerodynamic drag,...)
-# ------------------------------------------------------------------------------
class Driver:
+ """Driver of the time-dependent boundary layer equation solver (viscous solver)
+ """
+
def __init__(self,_Re, user_xtr, _timeParams, _Minfty, _case, _mapMesh):
self.Minfty = _Minfty
self.Re=_Re
@@ -60,8 +53,9 @@ class Driver:
pass
def dictionary(self):
- '''Create dictionary with the coordinates and the boundary layer parameters
- '''
+ """Create dictionary with the coordinates and the boundary layer parameters
+ """
+
wData = {}
wData['x'] = self.data[:,0]
wData['y'] = self.data[:,1]
@@ -92,78 +86,64 @@ class Driver:
return wData
def writeFile(self):
- '''Write the results in airfoil_viscous.dat
- '''
+ """Write the results in airfoil_viscous.dat
+ """
+
wData = self.dictionary()
toW = ['delta*', 'H', 'Cf', 'Ctau'] # list of variable to be written (should be a user input)
# Write
- print('writing file: airfoil_viscous.dat...', end = '')
+ print('Writing file: airfoil_viscous.dat...', end = '')
f = open('airfoil_viscous.dat', 'w+')
+
f.write('$Aerodynamic coefficients\n')
f.write(' Re Cd Cdp Cdf xtr_top xtr_bot\n')
- f.write('{0:15.6f} {1:15.6f} {2:15.6f} {3:15.6f} {4:15.6f} {5:15.6f}\n'.format(self.Re, wData['Cd'], wData['Cdp'], wData['Cdf'], wData['xtrTop'], wData['xtrBot']))
+ f.write('{0:15.6f} {1:15.6f} {2:15.6f} {3:15.6f} {4:15.6f} {5:15.6f}\n'.format(self.Re, wData['Cd'], wData['Cdp'],
+ wData['Cdf'], wData['xtrTop'],
+ wData['xtrBot']))
f.write('$Boundary layer variables\n')
f.write(' x y z x/c')
+
for s in toW:
f.write(' {0:>15s}'.format(s))
f.write('\n')
+
for i in range(len(self.data[:,0])):
f.write('{0:15.6f} {1:15.6f} {2:15.6f} {3:15.6f}'.format(wData['x'][i], wData['y'][i], wData['z'][i], wData['x/c'][i]))
for s in toW:
f.write(' {0:15.6f}'.format(wData[s][i]))
f.write('\n')
+
f.close()
- print('done!')
-
-
- def run(self, group):
- '''Primary component of the Diver. Receives 2 groups : 'airfoil' and 'wake'.
- Sorts data, initialize code components and calls the time tintegrator. Then sorts the data again,
- communicating it to the coupler.'''
-
- # ###################################### Driver.run() ######################################
- # ------------------------------------------------------------------------------------------
- # Input: -'group': - List of 2 elements (airfoil/wake) containing parameters
- # from the inviscid solver.
- #
- # Tasks: - Initializing: - var: Data container used throughout the code.
- # Contains the mesh, the solution, boundary
- # layer parameters etc.
- # - DiricheletBC: Used to compute the boudary conditions
- # of the airfoil (Twaithes) and the wake.
- # - gradComp: Gradient computation using various spacial
- # schemes.
- # - cSolver: Evaluates the boundary layer laws at one point
- # of the computational domain given a prescibed
- # solution or not. Uses the closures models.
- # - initializer: Provides the initial condition at the first
- # iteration
- # - tSolver: Time integration given a model, CFL and tolerences.
- # Used the same way for explicit/implicit time
- # marching.
- #
- # - Running: Runs the time integration.
- # Tries with lower CFL if the solver crashes
- # - Sorting: (Post Process): Calls the PostProcess class to sort the value and,
- # from the solution, compute/correct parameters that must
- # sent to the inviscid solver
- # Output: - Void.
- # Directly modifes the group that were initialy sent.
- # ------------------------------------------------------------------------------------------
- nFactor = 3
+ print('done.')
+
+
+ def run(self, groups):
+ """Runs viscous solver driver.
+ - Data preprocessing and initialization of the different components.
+ - Runs the pseudo-time solver
+ - Data postprocessing to ensure proper communication between the solvers.
+
+ Args
+ ----
+
+ - groups (data structure): Structures data containing the necessary information
+ related to each group (Airfoil upper surface, airfoil lower surface and wake).
+ """
+
+ nFactor = 2
# Merge upper, lower sides and wake
- dataIn=[None,None]
- for n in range(0, len(group)):
- group[n].connectListNodes, group[n].connectListElems, dataIn[n] = group[n].connectList()
- paramDict, bNodes, data, dx = GroupConstructor().mergeVectors(dataIn)
+ dataIn = [None,None]
+ for n in range(0, len(groups)):
+ groups[n].connectListNodes, groups[n].connectListElems, dataIn[n] = groups[n].connectList()
+ inMesh, paramDict, data, bNodes, inMeshbNodes, xxInit = GroupConstructor().mergeVectors(dataIn, self.mapMesh, nFactor)
- # Initialize variable container.
- var = Variables(data, paramDict, bNodes, self.xtrF, self.Ti, self.Re)
+ # Initialize data container.
+ var = Variables(data, paramDict, bNodes, 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] = data[:,8]
# Boundary condition (Airfoil/Wake).
DirichletBC=Conditions.Boundaries(self.Re)
@@ -185,7 +165,7 @@ class Driver:
if self.it == 0:
print('+------------------------------- Solver setup ---------------------------------+')
if self.mapMesh:
- print('| - Mesh mapped from {:<4.0f} to {:<4.0f} elements.{:>38s}'.format(len(initMesh), var.nN, '|'))
+ print('| - Mesh mapped from {:>5.0f} to {:>5.0f} elements.{:>35s}'.format(len(inMesh), var.nN, '|'))
print('| - Number of elements: {:<.0f}. {:>51s}'.format(var.nN,'|'))
if var.xtrF[0] is None:
print('| - Free transition on the upper side. {:>41}'.format('|'))
@@ -213,7 +193,7 @@ class Driver:
# Initial condition: 'Driver' stores the solution form the previous coupling iteration.
- if self.uPrevious is None or self.it < 6:
+ if self.uPrevious is None or self.it < 1:
# Typically for the first iteration.
# Initalize initial condition builder.
@@ -222,19 +202,15 @@ class Driver:
# Compute initial condition based on Twaithes solution.
initializer.initialConditions(var) # Set boundary condition
- # Default values of the time marching procedure if we start from an
- # initial solution that is far from the steady state.
- #tSolver.setCFL0(0.5)
- #tSolver.setmaxIt(30000)
-
print('| - Using Twaithes solution as initial guess. {:>34}'.format('|'))
else:
# If we use a solution previously obtained.
- var.u=self.uPrevious.copy()
-
- # Reset amplification factors
- var.u[2::var.nVar] = 0
+
+ var.u=self.uPrevious.copy() # Use previous solution.
+ DirichletBC.airfoilBC(var) # Reset boundary condition.
+ var.u[2::var.nVar] = 0 # Reset amplification factors.
+ var.u[3::var.nVar] = var.extPar[2::var.nExtPar] # Reset inviscid velocity.
print('| - Using previous solution as initial guess. {:>34}'.format('|'))
@@ -263,36 +239,36 @@ class Driver:
self.blParPrevious=var.blPar.copy()
# Compute blowing velocities and sort the boundary layer parameters.
- self.blScalars, AFblwVel, WKblwVel, AFblVectors, WKblVectors = PostProcess().sortParam(var)
- self.xtr=var.xtr.copy()
- self.blVectors=AFblVectors
+ self.blScalars, AFblwVel, WKblwVel, AFblVectors, WKblVectors, AFdeltaStarOut, WKdeltaStarOut = PostProcess().sortParam(var, self.mapMesh, inMesh, inMeshbNodes, xxInit)
+ self.xtr = var.xtr.copy()
+ self.blVectors = AFblVectors
# Group modification to ensure communication between the solvers.
- group[0].TEnd = [AFblVectors[1][var.bNodes[1]-1], AFblVectors[1][var.bNodes[2]-1]]
- group[0].HEnd = [AFblVectors[2][var.bNodes[1]-1], AFblVectors[2][var.bNodes[2]-1]]
- group[0].CtEnd = [AFblVectors[8][var.bNodes[1]-1], AFblVectors[8][var.bNodes[2]-1]]
- # Delete the stagnation point doublon for variables in the VII loop
- group[0].deltaStar = AFblVectors[0]
- group[0].xx = np.concatenate((var.xx[:var.bNodes[1]],var.xx[var.bNodes[1]:var.bNodes[2]]))
-
- # Sort the following results in reference frame
- group[0].deltaStar = group[0].deltaStar[group[0].connectListNodes.argsort()]
- group[0].xx = group[0].xx[group[0].connectListNodes.argsort()]
- group[0].u = AFblwVel[group[0].connectListElems.argsort()]
- self.data=data[:var.bNodes[2],:]
-
- group[1].deltaStar = WKblVectors[0]
- # The drag is measured at the end of the wake
- 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
-
- group[1].xx=var.xx[var.bNodes[2]:]
- # Sort the following results in reference frame
- group[1].deltaStar = group[1].deltaStar[group[1].connectListNodes.argsort()]
- group[1].xx = group[1].xx[group[1].connectListNodes.argsort()]
- group[1].u = WKblwVel[group[1].connectListElems.argsort()]
-
+ groups[0].TEnd = [AFblVectors[1][var.bNodes[1]-1], AFblVectors[1][var.bNodes[2]-1]]
+ groups[0].HEnd = [AFblVectors[2][var.bNodes[1]-1], AFblVectors[2][var.bNodes[2]-1]]
+ groups[0].CtEnd = [AFblVectors[8][var.bNodes[1]-1], AFblVectors[8][var.bNodes[2]-1]]
+ groups[0].deltaStar = AFdeltaStarOut
+ groups[0].xx = xxInit[:inMeshbNodes[2]]
+
+ # Sort the following results in reference frame.
+ groups[0].deltaStar = groups[0].deltaStar[groups[0].connectListNodes.argsort()]
+ groups[0].xx = groups[0].xx[groups[0].connectListNodes.argsort()]
+ groups[0].u = AFblwVel[groups[0].connectListElems.argsort()]
+ self.data = data[:var.bNodes[2],:]
+
+ # Drag is measured at the end of the wake.
+ 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 = xxInit[inMeshbNodes[2]:]
+
+ # Sort the following results in reference frame.
+ groups[1].deltaStar = groups[1].deltaStar[groups[1].connectListNodes.argsort()]
+ groups[1].xx = groups[1].xx[groups[1].connectListNodes.argsort()]
+ groups[1].u = WKblwVel[groups[1].connectListElems.argsort()]
+
"""self.writeFile()
- Validation().main(1,self.wData, WKblVectors, var.xGlobal[var.bNodes[2]:])"""
\ No newline at end of file
+ Validation().main(1,self.wData, WKblVectors, var.xGlobal[var.bNodes[2]:]) """
\ No newline at end of file
diff --git a/dart/viscous/Drivers/VGroupConstructor.py b/dart/viscous/Drivers/VGroupConstructor.py
index f93410e..bca6618 100755
--- a/dart/viscous/Drivers/VGroupConstructor.py
+++ b/dart/viscous/Drivers/VGroupConstructor.py
@@ -17,6 +17,7 @@
import numpy as np
import math
from matplotlib import pyplot as plt
+from scipy import interpolate
class GroupConstructor():
""" ---------- Group Constructor ----------
@@ -26,63 +27,132 @@ class GroupConstructor():
def __init__(self) -> None:
pass
- def mergeVectors(self,dataIn):
- '''Merges 3 groups upper/lower sides and wake.'''
+ def mergeVectors(self,dataIn, mapMesh, nFactor):
+ """Merges 3 groups upper/lower sides and wake.
+ """
+
+ # Saves the initial mesh.
+ inMesh = np.concatenate((dataIn[0][0][:,0], dataIn[0][1][1:,0], dataIn[1][:,0]))
+ inMeshbNodes = [0, len(dataIn[0][0][:,0]), len(dataIn[0][0][:,0]) + len(dataIn[0][1][1:,0])]
+
+ xxInit_up, _, _, _, _ = self.__getBLcoordinates(dataIn[0][0])
+ xxInit_lw, _, _, _, _ = self.__getBLcoordinates(dataIn[0][1])
+ xxInit_lw = np.delete(xxInit_lw, 0) # Delete stagnation point doublon.
+ xxInit_wk, _, _, _, _ = self.__getBLcoordinates(dataIn[1])
+ xxInit = np.concatenate((xxInit_up, xxInit_lw, xxInit_wk))
+
for k in range(len(dataIn)):
- if len(dataIn[k])==2: # Airfoil case
- xx_up, dv_up, vtExt_up, _, alpha_up= self.__getBLcoordinates(dataIn[k][0])
- xx_lw, dv_lw, vtExt_lw, _, alpha_lw= self.__getBLcoordinates(dataIn[k][1])
+
+ if len(dataIn[k]) == 2: # Airfoil case.
+
+ if mapMesh:
+ dataUpper = self.interpolateToFineMesh(dataIn[k][0], nFactor)
+ dataLower = self.interpolateToFineMesh(dataIn[k][1], nFactor)
+ else:
+ dataUpper = dataIn[k][0]
+ dataLower = dataIn[k][1]
+
+ xx_up, dv_up, vtExt_up, _, alpha_up = self.__getBLcoordinates(dataUpper)
+ xx_lw, dv_lw, vtExt_lw, _, alpha_lw = self.__getBLcoordinates(dataLower)
+
else: # Wake case
- xx_wk, dv_wk, vtExt_wk, _, alpha_wk= self.__getBLcoordinates(dataIn[k])
-
+ if mapMesh:
+ dataWake = self.interpolateToFineMesh(dataIn[k], nFactor)
+ else:
+ dataWake = dataIn[k]
+
+ xx_wk, dv_wk, vtExt_wk, _, alpha_wk = self.__getBLcoordinates(dataWake)
+
# Delete stagnation point doublon
- xx_lw=np.delete(xx_lw,0)
- dv_lw=np.delete(dv_lw,0)
- vtExt_lw=np.delete(vtExt_lw,0)
- alpha_lw=np.delete(alpha_lw,0)
-
- # Nodes at the boundary of the computational domain: [Stagnation point, First lower side point, First wake point]
- boundaryNodes=np.zeros(3,dtype=int)
- boundaryNodes[0]=0
- boundaryNodes[1]=len(xx_up)
- boundaryNodes[2]=len(xx_up)+len(xx_lw)
-
- param={'xx':None,'dv':None,'vtExt':None,'alpha':None}
- param['xx']=np.concatenate((xx_up,xx_lw,xx_wk))
- param['dv']=np.concatenate((dv_up,dv_lw,dv_wk))
- param['vtExt']=np.concatenate((vtExt_up,vtExt_lw,vtExt_wk))
- param['alpha']=np.concatenate((alpha_up,alpha_lw,alpha_wk))
-
- data=np.zeros((len(param['xx']),9))
- data=np.concatenate((dataIn[0][0],dataIn[0][1],dataIn[1]),axis=0)
- data=np.delete(data,boundaryNodes[1],0) # Delete stagnation point doublon
-
- dx=np.zeros(len(param['xx']))
- for i in range(1,len(param['xx'])):
- if i==boundaryNodes[1]:
- dx[i]=param['xx'][i]-param['xx'][0]
- elif i==boundaryNodes[2]:
- dx[i]=0
+ xx_lw = np.delete(xx_lw,0)
+ dv_lw = np.delete(dv_lw,0)
+ vtExt_lw = np.delete(vtExt_lw,0)
+ alpha_lw = np.delete(alpha_lw,0)
+
+ # Nodes at the boundary of the computational domain:
+ # [Stagnation point, First lower side point, First wake point]
+ boundaryNodes = np.zeros(3, dtype=int)
+ boundaryNodes[0] = 0
+ boundaryNodes[1] = len(xx_up)
+ boundaryNodes[2] = len(xx_up) + len(xx_lw)
+
+ param={'xx':None, 'dv':None, 'vtExt':None, 'alpha':None}
+ param['xx'] = np.concatenate((xx_up,xx_lw,xx_wk))
+ param['dv'] = np.concatenate((dv_up,dv_lw,dv_wk))
+ param['vtExt'] = np.concatenate((vtExt_up,vtExt_lw,vtExt_wk))
+ param['alpha'] = np.concatenate((alpha_up,alpha_lw,alpha_wk))
+
+ data = np.zeros((len(param['xx']),9))
+ data = np.concatenate((dataUpper, dataLower, dataWake), axis=0)
+ data = np.delete(data,boundaryNodes[1], axis = 0) # Delete stagnation point doublon
+
+ dx = np.zeros(len(param['xx']))
+ for i in range(1, len(param['xx'])):
+
+ if i == boundaryNodes[1]:
+
+ dx[i] = param['xx'][i] - param['xx'][0]
+
+ elif i == boundaryNodes[2]:
+
+ dx[i] = 0
+
else:
- dx[i]=param['xx'][i]-param['xx'][i-1]
- return param, boundaryNodes, data, dx
+ dx[i] = param['xx'][i] - param['xx'][i-1]
+
+ return inMesh, param, data, boundaryNodes, inMeshbNodes, xxInit
def __getBLcoordinates(self,data):
"""Transform the reference coordinates into airfoil coordinates
"""
- nN = len(data[:,0])
- nE = nN-1 # Nbr of element along the surface
- vt = np.zeros(nN)
- xx = np.zeros(nN) # Position along the surface of the airfoil
- dx = np.zeros(nE) # Distance along two nodes
- dv = np.zeros(nE) # Speed difference according to element
- alpha = np.zeros(nE) # Angle of the element for the rotation matrix
- for i in range(0,nE):
- alpha[i] = np.arctan2((data[i+1,1]-data[i,1]),(data[i+1,0]-data[i,0]))
- dx[i] = np.sqrt((data[i+1,0]-data[i,0])**2+(data[i+1,1]-data[i,1])**2)
- xx[i+1] = dx[i]+xx[i]
- vt[i] = (data[i,3]*np.cos(alpha[i]) + data[i,4]*np.sin(alpha[i])) # Tangential speed at node 1 element i
- vt[i+1] = (data[i+1,3]*np.cos(alpha[i]) + data[i+1,4]*np.sin(alpha[i])) # Tangential speed at node 2 element i
- dv[i] = (vt[i+1] - vt[i])/dx[i] # Velocity gradient according to element i. central scheme with half point
- return xx, dv, vt, nN, alpha
\ No newline at end of file
+ nN = len(data[:,0])
+ nE = nN-1 # Nbr of element along the surface.
+ vt = np.zeros(nN) # Outer flow velocity.
+ xx = np.zeros(nN) # Position along the surface of the airfoil.
+ dx = np.zeros(nE) # Distance along two nodes.
+ dv = np.zeros(nE) # Speed difference according to element.
+ alpha = np.zeros(nE) # Angle of the element for the rotation matrix.
+
+ for i in range(0, nE):
+
+ alpha[i] = np.arctan2((data[i+1, 1] - data[i, 1]), (data[i+1, 0]-data[i, 0]))
+ dx[i] = np.sqrt((data[i+1, 0] - data[i,0])**2 + (data[i+1, 1]-data[i, 1])**2)
+ xx[i+1] = dx[i] + xx[i]
+ vt[i] = (data[i, 3] * np.cos(alpha[i]) + data[i ,4] * np.sin(alpha[i])) # Tangential speed at node 1 element i
+ vt[i+1] = (data[i+1, 3] * np.cos(alpha[i]) + data[i+1, 4] * np.sin(alpha[i])) # Tangential speed at node 2 element i
+ dv[i] = (vt[i+1] - vt[i]) / dx[i] # Velocity gradient according to element i.
+ # central scheme with half point
+
+ return xx, dv, vt, nN, alpha
+
+ def interpolateToFineMesh(self, iData, nFactor):
+
+ inMesh = iData[:,0]
+
+ outMesh = []
+ for iPoint in range(len(inMesh) - 1):
+
+ dx = inMesh[iPoint + 1] - inMesh[iPoint]
+
+ for iRefinedPoint in range(nFactor):
+
+ outMesh = np.append(outMesh, inMesh[iPoint] + iRefinedPoint * dx / nFactor)
+
+ outMesh = np.append(outMesh, inMesh[-1])
+ outData = np.empty((len(outMesh), len(iData[0,:])))
+ outData[:,0] = outMesh
+
+ for iParam in range(1, len(iData[0, :])):
+
+ # Create interpolation function between the initial mesh
+ # and the function to evaluate.
+
+ f_interp = interpolate.interp1d(inMesh, iData[:,iParam])
+
+ # Map the function to the new mesh.
+ outData[:, iParam] = f_interp(outData[:, 0])
+
+ return outData
+
+
diff --git a/dart/viscous/Drivers/VPostProcess.py b/dart/viscous/Drivers/VPostProcess.py
index 0cdfd3b..3d844db 100644
--- a/dart/viscous/Drivers/VPostProcess.py
+++ b/dart/viscous/Drivers/VPostProcess.py
@@ -1,31 +1,63 @@
from matplotlib import pyplot as plt
import numpy as np
+from scipy import interpolate
class PostProcess:
def __init__(self):
pass
- def sortParam(self,var):
+ def sortParam(self,var, mapMesh, inMesh, inMeshbNodes, xxInit):
+
# Recompute deltaStar.
var.blPar[9::var.nBlPar] = var.u[0::var.nVar] * var.u[1::var.nVar]
- # Compute blowing velocity on each cell.
- blwVel=np.zeros(var.nN-1)
+ if mapMesh:
+
+ 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 i in range(1, var.nN):
+ for iPoint in range(1, var.nN):
- iPrev=0 if i==var.bNodes[1] else i-1
+ 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[i-1] = (var.extPar[i*var.nExtPar+1]*var.u[i*var.nVar+3] * var.blPar[i*var.nBlPar+9]
- - var.extPar[iPrev*var.nExtPar+1] * var.u[iPrev*var.nVar+3] * var.blPar[iPrev*var.nBlPar+9]) / (var.extPar[i*var.nExtPar+1] * (var.xx[i] - 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:]
+ # 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]))
+
+ print('blwVel bnodes2', blwVel[var.bNodes[2]])
+ print('blwVel bnodes2-1', blwVel[var.bNodes[2]-1])
+ blwVel = np.delete(blwVel, var.bNodes[2])
+ # Separate airfoil and wake blowing velocities.
+ AFblwVel=blwVel[:var.bNodes[2]-1]
+ WKblwVel=blwVel[var.bNodes[2]-1:]
# Normalize Friction and dissipation coefficients.
frictionCoeff=var.u[3::var.nVar]**2 * var.blPar[2::var.nBlPar]
@@ -67,4 +99,44 @@ class PostProcess:
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
+ 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
diff --git a/dart/viscous/Physics/VBIConditions.py b/dart/viscous/Physics/VBIConditions.py
index 4f8de9c..c4cbf59 100755
--- a/dart/viscous/Physics/VBIConditions.py
+++ b/dart/viscous/Physics/VBIConditions.py
@@ -44,48 +44,69 @@ class Initial:
regime prescibed by the inviscid solver.
- Define boundary conditions at the stagnation point and the first wake point
"""
- u0=np.zeros(var.nN*var.nVar)
- for n in range(3): # 3 components (upper, lower sides, wake)
+
+ u0=np.zeros(var.nN*var.nVar) # Vector containing the intial solution
+
+ # Loop over the three groups (upper, lower sides, wake)
+ for n in range(3):
+
if n==2:
- xx=var.xx[var.bNodes[n]:]
- vt=var.extPar[var.bNodes[n]*var.nExtPar+2::var.nExtPar]
- rho=var.extPar[var.bNodes[n]*var.nExtPar+1::var.nExtPar]
+
+ xx=var.xx[var.bNodes[n]:] # Airfoil local coordinates.
+ vt=var.extPar[var.bNodes[n]*var.nExtPar+2::var.nExtPar] # Inviscid velocities.
+ rho=var.extPar[var.bNodes[n]*var.nExtPar+1::var.nExtPar] # Air density in each cell.
+
else:
+
xx=var.xx[var.bNodes[n]:var.bNodes[n+1]]
vt=var.extPar[var.bNodes[n]*var.nExtPar+2:var.bNodes[n+1]*var.nExtPar:var.nExtPar]
rho=var.extPar[var.bNodes[n]*var.nExtPar+1:var.bNodes[n+1]*var.nExtPar:var.nExtPar]
+ # Reynolds number based on the arc length.
Rex=rho*vt*xx/self.mu_air
Rex[Rex==0]=1 # Stagnation point
- ThetaTw=xx/np.sqrt(Rex)
- # lambda : Dimensionless pressure-gradient parameter
- lambdaTw=np.zeros(len(ThetaTw))
- HTw=np.zeros(len(ThetaTw))
- for idx in range(0,len(lambdaTw)):
- idxPrev=idx-1
- lambdaTw[idx]=(rho[idx]*ThetaTw[idx]**2)/self.mu_air*(vt[idx]-vt[idxPrev])#/(xx[idx]-xx[idxPrev])
+
+ ThetaTw=xx/np.sqrt(Rex) # Momentum thickness.
+ lambdaTw=np.zeros(len(ThetaTw)) # Dimensionless pressure gradient parameter.
+ HTw=np.zeros(len(ThetaTw)) # Shape factor of the boundary layer.
+
+ for iMarker in range(len(lambdaTw)):
+
+ idxPrev = iMarker-1
+ lambdaTw[iMarker] = (rho[iMarker] * ThetaTw[iMarker]**2) / self.mu_air * (vt[iMarker] - vt[idxPrev])#/(xx[idx]-xx[idxPrev])
+
# Empirical correlations for H=f(lambda)
- if 0<=lambdaTw[idx]<0.1:
- HTw[idx]=2.61-3.75*lambdaTw[idx]+5.24*lambdaTw[idx]**2
- elif -0.1<lambdaTw[idx]<0:
- HTw[idx]=2.088+0.0731/(lambdaTw[idx]+0.14)
+ if 0 <= lambdaTw[iMarker] < 0.1:
+
+ HTw[iMarker] = 2.61 - 3.75*lambdaTw[iMarker] + 5.24*lambdaTw[iMarker]**2
+
+ elif -0.1<lambdaTw[iMarker]<0:
+
+ HTw[iMarker] = 2.088 + 0.0731 / (lambdaTw[iMarker] + 0.14)
if n==2:
- u0[var.bNodes[n]*var.nVar+0::var.nVar]=ThetaTw # Theta
- u0[var.bNodes[n]*var.nVar+1::var.nVar]=HTw # H
+
+ u0[var.bNodes[n]*var.nVar+0::var.nVar]=ThetaTw # Theta
+ u0[var.bNodes[n]*var.nVar+1::var.nVar]=HTw # H
+
else:
- u0[var.bNodes[n]*var.nVar+0:var.bNodes[n+1]*var.nVar:var.nVar]=ThetaTw # Theta
- u0[var.bNodes[n]*var.nVar+1:var.bNodes[n+1]*var.nVar:var.nVar]=HTw # H
- u0[1::var.nVar]=2 # H
- u0[2::var.nVar]=0 # N
- u0[3::var.nVar]=var.extPar[2::var.nExtPar] # Vt
- u0[4::var.nVar]=0.001 # Ct
+ u0[var.bNodes[n]*var.nVar+0:var.bNodes[n+1]*var.nVar:var.nVar]=ThetaTw # Theta
+ u0[var.bNodes[n]*var.nVar+1:var.bNodes[n+1]*var.nVar:var.nVar]=HTw # H
+ u0[1::var.nVar]=2 # H
+ u0[2::var.nVar]=0 # N
+ u0[3::var.nVar]=var.extPar[2::var.nExtPar] # Vt
+ u0[4::var.nVar]=0.001 # Ct
+
+ # Set initial condition.
var.u=u0
+
+ # Impose boudary conditions.
self.boundCond.airfoilBC(var)
self.boundCond.wakeBC(var)
+ pass
class Boundaries:
@@ -97,33 +118,46 @@ class Boundaries:
----------
Re: Flow Reynolds number.
"""
+
def __init__(self, _Re):
- self.Re=_Re
+
+ self.Re = _Re
# Airfoil boundary conditions
def airfoilBC(self, var):
"""Boundary conditions at the stagnation point.
- Twaithes solution method values used for the SP.
+ Twaithes solution method values used for the
+ stagnation point.
"""
+
if var.dv[0] > 0:
- T0 = np.sqrt(0.075/(self.Re*var.dv[0]))
+
+ T0 = np.sqrt(0.075 / (self.Re * var.dv[0]))
+
else:
+
T0 = 1e-6
- H0 = 2.23 # One parameter family Falkner Skan
- n0 = 0
- ue0=var.extPar[2]
- Ct0 = 0
-
- # Stagnation point
- var.blPar[var.bNodes[0]*var.nBlPar+0] = ue0*T0*var.Re # Ret
- if var.blPar[var.bNodes[0]*var.nBlPar+0]<1:
- var.blPar[var.bNodes[0]*var.nBlPar+0]=1
- ue0=var.blPar[var.bNodes[0]*var.nBlPar+0]/(T0*var.Re)
- var.blPar[var.bNodes[0]*var.nBlPar+9]=T0*H0 # deltaStar
- var.blPar[(var.bNodes[0]*var.nBlPar+1):(var.bNodes[0]*var.nBlPar+7)]=Closures.laminarClosure(T0, H0, var.blPar[var.bNodes[0]*var.nBlPar+0],var.extPar[var.bNodes[0]*var.nExtPar+0], 'airfoil')
- # Cteq stays = 0
+
+ H0 = 2.23 # One parameter family Falkner Skan.
+ n0 = 0 # Initial log of amplification factor.
+ ue0=var.extPar[2] # Inviscid flow velocity.
+ Ct0 = 0 # Stagnation point is laminar.
+
+ # Compute Reynolds number based on the momentum thickness
+ var.blPar[var.bNodes[0]*var.nBlPar + 0] = ue0*T0*var.Re
+
+ if var.blPar[var.bNodes[0]*var.nBlPar+0] < 1:
+ var.blPar[var.bNodes[0]*var.nBlPar+0] = 1
+ ue0 = var.blPar[var.bNodes[0]*var.nBlPar + 0] / (T0*var.Re)
+
+ var.blPar[var.bNodes[0]*var.nBlPar + 9] = T0 * H0 # deltaStar
+ var.blPar[(var.bNodes[0]*var.nBlPar + 1):(var.bNodes[0]*var.nBlPar+7)]=Closures.laminarClosure(T0, H0,
+ var.blPar[var.bNodes[0]*var.nBlPar+0],
+ var.extPar[var.bNodes[0]*var.nExtPar+0],
+ 'airfoil')
#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]
pass
@@ -133,23 +167,33 @@ class Boundaries:
based on the parameters at the last point on the upper and lower sides.
Drela (1989).
"""
+
xEnd_up=var.u[(var.bNodes[1]-1)*var.nVar:var.bNodes[1]*var.nVar]
xEnd_lw=var.u[(var.bNodes[2]-1)*var.nVar:var.bNodes[2]*var.nVar]
- T0=xEnd_up[0]+xEnd_lw[0] # (Theta_up+Theta_low)
- H0=(xEnd_up[0]*xEnd_up[1]+xEnd_lw[0]*xEnd_lw[1])/T0 # ((Theta*H)_up+(Theta*H)_low)/(Theta_up+Theta_low)
+
+ T0=xEnd_up[0]+xEnd_lw[0] # (Theta_up+Theta_low).
+ H0=(xEnd_up[0]*xEnd_up[1]+xEnd_lw[0]*xEnd_lw[1])/T0 # ((Theta*H)_up+(Theta*H)_low)/(Theta_up+Theta_low).
n0=9
- ue0=var.extPar[var.bNodes[2]*var.nExtPar+2]
- Ct0=(xEnd_up[0]*xEnd_up[4]+xEnd_lw[0]*xEnd_lw[4])/T0 # ((Ct*Theta)_up+(Theta)_low)/(Ct*Theta_up+Theta_low)
+ ue0=var.extPar[var.bNodes[2]*var.nExtPar+2] # Inviscid velocity.
+ Ct0=(xEnd_up[0]*xEnd_up[4]+xEnd_lw[0]*xEnd_lw[4])/T0 # ((Ct*Theta)_up+(Theta)_low)/(Ct*Theta_up+Theta_low).
+ # Reynolds number based on the momentum thickness.
var.blPar[var.bNodes[2]*var.nBlPar+0] = ue0*T0*var.Re
- if var.blPar[var.bNodes[2]*var.nBlPar+0]<1:
- var.blPar[var.bNodes[2]*var.nBlPar+0]=1
- ue0=var.blPar[var.bNodes[2]*var.nBlPar+0]/(T0*var.Re)
+
+ if var.blPar[var.bNodes[2]*var.nBlPar+0] < 1:
+ var.blPar[var.bNodes[2]*var.nBlPar+0] = 1
+ ue0 = var.blPar[var.bNodes[2]*var.nBlPar+0] / (T0*var.Re)
+
var.blPar[var.bNodes[2]*var.nBlPar+9]=T0*H0 # deltaStar
+
# Laminar reset
- var.blPar[(var.bNodes[2]*var.nBlPar+1):(var.bNodes[2]*var.nBlPar+7)]=Closures.laminarClosure(T0, H0, var.blPar[var.bNodes[2]*var.nBlPar+0],var.extPar[var.bNodes[2]*var.nExtPar+0], 'wake')
+ var.blPar[(var.bNodes[2]*var.nBlPar+1):(var.bNodes[2]*var.nBlPar+7)]=Closures.laminarClosure(T0, H0,
+ var.blPar[var.bNodes[2]*var.nBlPar+0],
+ var.extPar[var.bNodes[2]*var.nExtPar+0],
+ 'wake')
#var.blPar[var.bNodes[2]*var.nBlPar+7]=0 # Cteq stays = 0
#var.blPar[var.bNodes[2]*var.nBlPar+8]=(var.blPar[var.bNodes[2]*var.nBlPar+4]/2)*(1-4*(var.blPar[var.bNodes[2]*var.nBlPar+6]-1)/(3*H0)) if H0!=0 else 0 # Equivalent normalized wall slip velocity
- var.u[var.bNodes[2]*var.nVar:(var.bNodes[2]+1)*var.nVar]=[T0, H0, n0, ue0, Ct0]
+ # Set boundary condition.
+ var.u[var.bNodes[2]*var.nVar:(var.bNodes[2]+1)*var.nVar] = [T0, H0, n0, ue0, Ct0]
pass
diff --git a/dart/viscous/Physics/VVariables.py b/dart/viscous/Physics/VVariables.py
index a5a9f62..094e74d 100755
--- a/dart/viscous/Physics/VVariables.py
+++ b/dart/viscous/Physics/VVariables.py
@@ -16,6 +16,7 @@
# ------------------------------------------------------------------------------
import numpy as np
from matplotlib import pyplot as plt
+from fwk.coloring import ccolors
class Variables:
"""Data container for viscous solver.
@@ -66,38 +67,39 @@ class Variables:
"""
def __init__(self, data, paramDict, bNodes, xtrF, Ti, Re):
- self.__Ti=Ti
- self.Re=Re
- self.xx=paramDict['xx']
- self.nN=len(self.xx)
+ self.__Ti = Ti
+ self.Re = Re
+ self.xx = paramDict['xx']
+ self.nN = len(self.xx)
- self.bNodes=bNodes
+ self.bNodes = bNodes
# Solution
- self.nVar=5
- self.u=np.zeros(self.nVar*self.nN)
+ self.nVar = 5
+ self.u = np.zeros(self.nVar*self.nN)
# Parameters of the boundary layers
# [ 0 , 1, 2, 3, 4, 5, 6, 7, 8, 9]
# [Ret, Cd, Cf, delta, Hstar, Hstar2, Hk, Cteq, Us, deltaStar]
- self.nBlPar=10
- self.blPar=np.zeros(self.nBlPar*self.nN)
- self.Ret_previous=self.blPar[0::self.nVar]
+ self.nBlPar = 10
+ self.blPar = np.zeros(self.nBlPar*self.nN)
+ self.Ret_previous = self.blPar[0::self.nVar]
- self.dueOld_dx=np.zeros(self.nN)
- self.ddeltaStarOld_dx=np.zeros(self.nN)
+ self.dueOld_dx = np.zeros(self.nN)
+ self.ddeltaStarOld_dx = np.zeros(self.nN)
# Parameters of the external (inviscid) flow
# [ 0, 1, 2, 3, 4, 5]
# [Me, rhoe, vtExt, deltaStarExt, xxExt ReExt]
- self.nExtPar=6
- self.extPar=np.zeros(self.nExtPar*self.nN)
- self.extPar[0::self.nExtPar]=data[:,5]
- self.extPar[1::self.nExtPar]=data[:,6]
- self.extPar[2::self.nExtPar]=paramDict['vtExt']
- self.extPar[3::self.nExtPar]=data[:,7]
- self.extPar[4::self.nExtPar]=paramDict['xx']
- self.dv=paramDict['dv']
- self.xGlobal=data[:,0]
+ self.nExtPar = 6
+
+ self.extPar = np.zeros(self.nExtPar*self.nN)
+ self.extPar[0::self.nExtPar] = data[:,5]
+ self.extPar[1::self.nExtPar] = data[:,6]
+ self.extPar[2::self.nExtPar] = paramDict['vtExt']
+ self.extPar[3::self.nExtPar] = data[:,7]
+ self.extPar[4::self.nExtPar] = paramDict['xx']
+ self.dv = paramDict['dv']
+ self.xGlobal = data[:,0]
self.flowRegime = np.zeros(self.nN, dtype=int)
diff --git a/dart/viscous/README.md b/dart/viscous/README.md
index fe02415..71ba0a1 100644
--- a/dart/viscous/README.md
+++ b/dart/viscous/README.md
@@ -13,7 +13,7 @@ The two flow domains are coupled through a quasi-simulataneous interaction law t
- Time-dependant equations solution through a pseudo-time marching procedure.
## Solver operation map
-The viscous inviscid interaction can be launched from /dart/tests/bli_xxx.py. The class in Master/Coupler.py is first called with the required arguments and manages the two solvers to work togheter. The steady-state equations solver is called from /Solvers/Steady/StdSolver.py. The time-dependent equations solver is called from /Drivers/VDriver.py which initializes the different components of the code and starts the viscous flow computation using a implicit Euler pseudo-time marching procedures.
+The viscous inviscid interaction can be launched from /dart/tests/bli_xxx.py. The class in Master/Coupler.py is first called with the required arguments and manages the two solvers to work togheter. The steady-state equations solver is called from /Solvers/Steady/StdSolver.py. The time-dependent equations solver is called from /Drivers/VDriver.py which initializes the different components of the code and starts the viscous flow computation using an implicit Euler pseudo-time marching procedure.
## References
Crovato Adrien, [Steady Transonic Aerodynamic and Aeroelastic Modeling for Preliminary Aircraft Design](http://hdl.handle.net/2268/251906), PhD thesis, University of Liège, 2020.
diff --git a/dart/viscous/Solvers/Steady/StdSolver.py b/dart/viscous/Solvers/Steady/StdSolver.py
index 9185fa4..5ae44a3 100755
--- a/dart/viscous/Solvers/Steady/StdSolver.py
+++ b/dart/viscous/Solvers/Steady/StdSolver.py
@@ -465,7 +465,7 @@ class Solver:
self.blScalars = np.add(ublScalars, lblScalars)
self.blVectors = np.hstack((ublVectors, lblVectors))
blwVel = np.concatenate((ublwVel, lblwVel))
- print('AF',len(blwVel))
+ print('AFblwVel',len(blwVel))
self.xtr = [xtrTop, xtrBot]
# Boundary conditions for the wake
group.TEnd = [ublVectors[1][-1], lblVectors[1][-1]]
@@ -474,12 +474,14 @@ class Solver:
# Delete the stagnation point doublon for variables in the VII loop
lblVectors[0] = np.delete(lblVectors[0],0,0) #deltastar
lxx = np.delete(lxx,0,0) # airfoil coordinates
+ print('AFdeltaStarOut', len(np.concatenate((ublVectors[0], lblVectors[0]))))
group.deltaStar = np.concatenate((ublVectors[0], lblVectors[0]))
group.xx = np.concatenate((uxx, lxx))
# Compute BLE for the wake
else:
blwVel, _, group.xx, blScalars, blVectors = self.__boundaryLayer(data, group)
- print('WK',len(blwVel))
+ print('WKblwVel',len(blwVel))
+ print('WKdeltaStarOut', len(blVectors[0]))
group.deltaStar = blVectors[0]
# The drag is measured at the end of the wake
self.Cd = blScalars[0]
@@ -489,5 +491,6 @@ class Solver:
self.Cdf = self.blScalars[1]
# Sort the following results in reference frame
group.deltaStar = group.deltaStar[group.connectListNodes.argsort()]
+ print('xxInit',len(group.xx))
group.xx = group.xx[group.connectListNodes.argsort()]
group.u = blwVel[group.connectListElems.argsort()]
diff --git a/dart/viscous/Solvers/VSolver.py b/dart/viscous/Solvers/VSolver.py
index 5c3a0f2..91bf0b9 100644
--- a/dart/viscous/Solvers/VSolver.py
+++ b/dart/viscous/Solvers/VSolver.py
@@ -44,12 +44,12 @@ class flowEquations:
if iMarker == dataCont.bNodes[2]:
return np.zeros(dataCont.nVar)
- timeMatrix=np.empty((dataCont.nVar, dataCont.nVar))
- spaceMatrix=np.empty(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)
+ 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
@@ -62,14 +62,14 @@ class flowEquations:
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+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],
+ 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)
@@ -118,7 +118,7 @@ class flowEquations:
# 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
+ # 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)
# ---------------------------------------------------------
@@ -156,7 +156,7 @@ class flowEquations:
if np.linalg.det(timeMatrix) == 0:
raise RuntimeError('Singular time matrix.')
- sysRes = np.linalg.solve(timeMatrix,-spaceMatrix).flatten()
+ sysRes = np.linalg.solve(timeMatrix, spaceMatrix).flatten()
return sysRes
@@ -220,18 +220,7 @@ class sysMatrix(flowEquations):
self._nMarkers = _nMarkers
pass
- class Jacobian:
-
- def _initialize(nVar, nMarkers) -> np.ndarray:
- return np.zeros((nVar * nMarkers, nVar * nMarkers))
-
- class linSysRes:
-
- def _initialize(nVar, nMarkers) -> np.ndarray:
- return np.zeros(nVar * nMarkers)
-
-
- def buildJacobian(self, dataCont, marker = -1, order = 1, linSysRes = None, eps = 1e-10) -> np.ndarray:
+ def buildJacobian(self, dataCont, marker, sysRes, order = 1, eps = 1e-10) -> np.ndarray:
"""Contruct the Jacobian used for solving the liner system.
Args:
@@ -254,61 +243,42 @@ class sysMatrix(flowEquations):
"""
- # Build Jacobian for whole domain
- if marker == -1:
- sysMatSize = self._nMarker * self._nMarkers
- Jacobian = np.zeros((sysMatSize, sysMatSize))
-
- for iMarker in range(1,dataCont.nMarkers):
- xUp = dataCont.u[iMarker*dataCont.nMarkers : (iMarker+1)*dataCont.nMarkers].copy()
- xDwn = dataCont.u[iMarker*dataCont.nMarkers : (iMarker+1)*dataCont.nMarkers].copy()
-
- for iVar in range(dataCont.nVar):
- xSave = xUp
- xUp += eps
- xDwn -= eps
-
# Build Jacobian for one point
- else:
+ if order == 1 and sysRes is not None:
- if order == 1 and linSysRes is not None:
-
- JacMatrix = np.zeros((dataCont.nVar, dataCont.nVar))
- xUp = dataCont.u[marker*dataCont.nVar : (marker+1)*dataCont.nVar].copy()
+ 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) - linSysRes) / eps
+ JacMatrix[:,iVar] = (self._blLaws(dataCont, marker, x=xUp) - sysRes) / eps
- xUp[iVar] = varSave
- return JacMatrix
-
- elif order == 1:
- raise RuntimeError('First order Jacobian determination requires the linear system residuals')
+ 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 = -1) -> np.ndarray:
+ def buildResiduals(self, dataCont, marker) -> np.ndarray:
"""Construct residuals of the linear system.
Args:
@@ -319,9 +289,5 @@ class sysMatrix(flowEquations):
- marker (int, optional): Cell id on the computational domain. default: -1 for whole domain.
"""
- # Build Residual for the whole domain
- if marker == -1:
- a=1
-
# Build Residual for one point
- else: return self._blLaws(dataCont, marker, x=None)
\ No newline at end of file
+ return self._blLaws(dataCont, marker, x=None)
\ No newline at end of file
diff --git a/dart/viscous/Solvers/VTimeSolver.py b/dart/viscous/Solvers/VTimeSolver.py
index b14541d..1b4fe23 100644
--- a/dart/viscous/Solvers/VTimeSolver.py
+++ b/dart/viscous/Solvers/VTimeSolver.py
@@ -14,7 +14,6 @@
# ------------------------------------------------------------------------------
# Imports
# ------------------------------------------------------------------------------
-from typing import SupportsRound
from dart.viscous.Physics.VClosures import Closures
from matplotlib import pyplot as plt
@@ -124,7 +123,7 @@ class implicitEuler(timeConfig):
self.__CFL0=_cfl0
self.__maxIt = 20000
- self.__tol = 1e-6
+ self.__tol = 1e-8
self.solver=_solver
self.transSolver = transition
self.wakeBC = wakeBC
@@ -141,83 +140,92 @@ class implicitEuler(timeConfig):
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
- ---------
- dU/dt = F(U)
- <=> (U_(n+1)-U_n)/dt = F(U_(n+1)) (Time disc.)
- <=> (J-I)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.
+ """ ---------- 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).
+ - Asks 'solver' for spacial operator F and Jacobian J and uses a
+ direct linear solver (LU decomposition).
- - Corrects undesirable results.
+ - Corrects undesirable results.
- - Asks the solver for transtion position/BC and wake BC to be updated at each iteration.
- """
-
- # Flag to diplay the simulation evolution.
- convergedPoints = np.zeros(var.nN)
- convergedPoints[0] = 1 # Boundary condition
+ - 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)
+ # Initialize the flow regime on each cell.
+ self.transSolver.initTransition(var)
- # Transition is locked when it is found on one side. We do not look at the value of N
- # for the remaining points on this side.
- lockTrans = 0
+ lockTrans = 0 # Flag to remember if the transition is already found or not.
- for iMarker in range(1, var.nN):
- displayState = 0
+ for iMarker in range(1, nN):
+ displayState = 0 # Flag to output simulation state.
- if not lockTrans and 1 < iMarker < var.bNodes[2]:
+ 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)*var.nBlPar + 6]-1))**0.43
- + 0.7*(np.tanh(14*(1/(var.blPar[(iMarker-1)*var.nBlPar + 6]-1))-9.24)+1.0)
+ # 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)*var.nBlPar + 0] > 0: # Avoid log of something <= 0
+ if var.blPar[(iMarker-1)*nBlPar + 0] > 0: # Avoid log of something <= 0
- if np.log10(var.blPar[(iMarker-1)*var.nBlPar + 0]) <= logRet_crit - 0.08:
+ if np.log10(var.blPar[(iMarker-1)*nBlPar + 0]) <= logRet_crit - 0.08:
- var.u[iMarker*var.nVar + 2] = 0
+ var.u[iMarker*var.nVar + 2] = 0
- # Upper side start.
- if iMarker == var.bNodes[0]+1:
- print('| - Computing upper side. {:>54}'.format('|'))
+ # Upper side start.
+ if iMarker == bNodes[0]+1:
+ print('| - Computing upper side. {:>54}'.format('|'))
- # Lower side start.
- if iMarker == var.bNodes[1]:
- lockTrans = 0
- print('| - Computing lower side. {:>54}'.format('|'))
+ # Lower side start.
+ if iMarker == bNodes[1]:
+ lockTrans = 0
+ print('| - Computing lower side. {:>54}'.format('|'))
- # Wake start
- if iMarker == var.bNodes[2]: # First wake point
+ # Wake start
+ if iMarker == bNodes[2]: # First wake point
# Wake boundary condition.
self.wakeBC(var)
@@ -228,104 +236,123 @@ class implicitEuler(timeConfig):
# Ignore remaining instructions for the first wake point.
continue
- if not lockTrans and iMarker-1 < var.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)*var.nBlPar+1:(iMarker-1)*var.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)*var.nBlPar+0],
- var.extPar[(iMarker-1)*var.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)*var.nBlPar+1:(iMarker-1)*var.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)*var.nBlPar+0],
- var.extPar[(iMarker-1)*var.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)*var.nBlPar+1:(iMarker-1)*var.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)*var.nBlPar+0],
- var.extPar[(iMarker-1)*var.nExtPar+0],
- 'airfoil')
-
- # Call Newton solver for the point
- convergedPoints[iMarker] = self.newtonSolver(var, iMarker, displayState)
-
- return convergedPoints
+ 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()
-
- sysRes = self.solver.buildResiduals(var, iMarker)
- normSysRes0 = np.linalg.norm(sysRes)
-
- # Initial CFL and time step.
- CFL = self.__CFL0
- CFLAdapt = 1
- jacEval = 3
- iMarkerPrev = 0 if iMarker == var.bNodes[1] else iMarker -1
- dx = var.xx[iMarker] - var.xx[iMarkerPrev]
-
- dt = self.computeTimeStep(CFL, dx, var.extPar[iMarker*var.nExtPar + 2])
- IMatrix=np.identity(var.nVar) / dt
+ # 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 = 3 # 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
- innerIters = 0
+ 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:
+ 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)
+
+ # 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(sysRes)
+ normSysRes = np.linalg.norm(slnIncr/dt - sysRes)
if displayState and innerIters % 1000 == 0:
self.outputState(var, iMarker, innerIters, normSysRes, normSysRes0, CFL, 'yellow')
@@ -336,39 +363,35 @@ class implicitEuler(timeConfig):
self.outputState(var, iMarker, innerIters, normSysRes, normSysRes0, CFL)
break
- # Jacobian computation.
-
- if innerIters % jacEval == 0:
- Jacobian=self.solver.buildJacobian(var, iMarker, linSysRes = sysRes, eps = 1e-10)
-
- # Linear solver.
-
- slnIncr = np.linalg.solve((IMatrix - Jacobian), sysRes)
-
- # Increment solution and correct absurd transcient results.
-
- var.u[iMarker*var.nVar : (iMarker+1)*var.nVar] += slnIncr
- var.u[iMarker*var.nVar + 1] = max(var.u[iMarker*var.nVar + 1], 1.0005)
-
except KeyboardInterrupt:
quit()
- except Exception as e:
+ 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} {:>30s}'.format(var.xGlobal[iMarker], iMarker,'|'), ccolors.ANSI_RESET)
- print(e)
+ 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*var.nVar : (iMarker+1)*var.nVar] = var.u[(iMarker-1)*var.nVar : (iMarker)*var.nVar]
+ 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 / 3,0.01),1)
+ CFL = min(max(CFL / 2,0.01),1)
print('Lowering CFL: {:<.3f}'.format(CFL))
- CFLAdapt = 0
- dt = self.computeTimeStep(CFL, dx, var.extPar[iMarker*var.nExtPar + 2])
- IMatrix = np.identity(var.nVar) / dt
+ CFLAdapt = 1
+ dt = self.computeTimeStep(CFL, dx, iInvVel)
+ IMatrix = np.identity(nVar) / dt
+
+ sysRes = self.solver.buildResiduals(var, iMarker)
jacEval = 1
displayState = 1
@@ -386,8 +409,8 @@ class implicitEuler(timeConfig):
CFL = self.__CFL0* (normSysRes0/normSysRes)**0.7
CFL = max(CFL, 0.1)
- dt = self.computeTimeStep(CFL, dx, var.extPar[iMarker*var.nExtPar + 2])
- IMatrix=np.identity(var.nVar) / dt
+ dt = self.computeTimeStep(CFL, dx, iInvVel)
+ IMatrix=np.identity(nVar) / dt
innerIters+=1
@@ -395,29 +418,31 @@ class implicitEuler(timeConfig):
else:
# Maximum number of iterations reached.
if normSysRes >= 1e-3 * normSysRes0:
- print('|', ccolors.ANSI_RED + 'Too many iteration at position {:<.4f}. normLinSysRes = {:<.3f}. {:>30s}'.format(var.xGlobal[iMarker],
- np.log10(normSysRes/normSysRes0),'|'),
+ print('|', ccolors.ANSI_RED + 'Too many iteration at position {:<.4f}. normLinSysRes = {:<.3f}. {:>30s}'
+ .format(var.xGlobal[iMarker],
+ np.log10(normSysRes/normSysRes0),'|'),
ccolors.ANSI_RESET)
- var.u[iMarker*var.nVar:(iMarker+1)*var.nVar] = u0[(iMarker)*var.nVar: (iMarker+1)*var.nVar].copy()
+ 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*var.nBlPar+1:iMarker*var.nBlPar+7]=Closures.laminarClosure(var.u[iMarker*var.nVar + 0], var.u[iMarker*var.nVar + 1],
- var.blPar[iMarker*var.nBlPar+0], var.extPar[iMarker*var.nExtPar+0],
+ 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*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], name)
+ 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}. normLinSysRes = {:<.3f}. {:>30s}'.format(var.xGlobal[iMarker],
- np.log10(normSysRes/normSysRes0),'|'),
- ccolors.ANSI_RESET)
+ print('|', ccolors.ANSI_YELLOW+ 'Too many iteration at position {:<.4f}. normLinSysRes = {:<.3f}. {:>30s}'
+ .format(var.xGlobal[iMarker],
+ np.log10(normSysRes/normSysRes0),'|'),
+ ccolors.ANSI_RESET)
return converged
\ No newline at end of file
--
GitLab