diff --git a/dart/tests/bliNonLift.py b/dart/tests/bliNonLift.py
index 456ab3679ed8cae9e6eb6228305ab11d645fe266..8ac5aef23112b3dc06d028433b3380dd9252ca57 100755
--- a/dart/tests/bliNonLift.py
+++ b/dart/tests/bliNonLift.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
@@ -60,7 +60,7 @@ def main():
# define flow variables
Re = 6e6
- alpha = 0.*math.pi/180
+ alpha = 5.*math.pi/180
M_inf = 0.15
#user_xtr=[0.41,0.41]
@@ -69,14 +69,11 @@ def main():
user_Ti=None
mapMesh = 1
+ nFactor = 3
# Time solver parameters
Time_Domain = True # True for unsteady solver, False for steady solver
- Cfl = 3
- spaceMarching=True
-
- if Time_Domain is True:
- timeParams={'CFL0': Cfl, 'spaceMarching':spaceMarching}
+ CFL0 = 3
# ========================================================================================
@@ -89,7 +86,7 @@ def main():
# mesh the geometry
print(ccolors.ANSI_BLUE + 'PyMeshing...' + ccolors.ANSI_RESET)
tms['msh'].start()
- pars = {'xLgt' : 5, 'yLgt' : 5, 'msF' : 1, 'msTe' : 0.01, 'msLe' : 0.01}
+ pars = {'xLgt' : 5, 'yLgt' : 5, 'msF' : 1, 'msTe' : 0.01, 'msLe' : 0.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()
@@ -107,8 +104,8 @@ def main():
# 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)
diff --git a/dart/tests/bliSeparated.py b/dart/tests/bliSeparated.py
index 50181f5cfffdb58a6932be5eaa2bce4b3282daab..47d4fd15fd182b09858576bbf34d0d87493cb279 100755
--- a/dart/tests/bliSeparated.py
+++ b/dart/tests/bliSeparated.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
@@ -66,16 +66,12 @@ def main():
user_xtr=[None,None]
user_Ti=None
- mapMesh = 0
+ mapMesh = 1
+ nFactor = 7
# Time solver parameters
Time_Domain=True # True for unsteady solver, False for steady solver
- Cfl = 1
- spaceMarching=True
- SteadyInitU=False # Initial condition prescribed by the steady solver
-
- if Time_Domain is True:
- timeParams={'CFL0': Cfl, 'SteadyInitU':SteadyInitU, 'spaceMarching':spaceMarching}
+ CFL0 = 1
# ========================================================================================
@@ -89,7 +85,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.0005}
+ pars = {'xLgt' : 5, 'yLgt' : 5, 'msF' : 1, 'msTe' : 0.01, 'msLe' : 0.001}
msh, gmshWriter = floD.mesh(dim, 'models/n0012.geo', pars, ['field', 'airfoil', 'downstream'])
tms['msh'].stop()
@@ -105,8 +101,8 @@ def main():
# 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)
diff --git a/dart/tests/bliTransonic.py b/dart/tests/bliTransonic.py
index 21a5edec15e7e8e995c40d418c82000b690890ef..c923b03533dc7912bfa36934d731cdef86ad79a2 100755
--- a/dart/tests/bliTransonic.py
+++ b/dart/tests/bliTransonic.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
@@ -69,15 +69,11 @@ def main():
user_Ti=None
mapMesh = 1
+ nFactor = 5
# Time solver parameters
Time_Domain=True # True for unsteady solver, False for steady solver
- Cfl = 1
- spaceMarching=True
- SteadyInitU=False # Initial condition prescribed by the steady solver
-
- if Time_Domain is True:
- timeParams={'CFL0': Cfl, 'SteadyInitU':SteadyInitU, 'spaceMarching':spaceMarching}
+ CFL0 = 1
# ========================================================================================
@@ -109,8 +105,8 @@ def main():
# 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
diff --git a/dart/viscous/Drivers/DDriver.py b/dart/viscous/Drivers/DDriver.py
new file mode 100644
index 0000000000000000000000000000000000000000..c2bcf7d547a0ba7f3ca197ca9dd94e63ab8d2715
--- /dev/null
+++ b/dart/viscous/Drivers/DDriver.py
@@ -0,0 +1,276 @@
+################################################################################
+# #
+# DARTFlo viscous module file #
+# Driver : Communication with coupler and initialization #
+# of some components of the viscous computation part #
+# #
+# Author: Paul Dechamps #
+# Université de Liège #
+# Date: 2021.09.10 #
+# #
+################################################################################
+
+
+# ------------------------------------------------------------------------------
+# Imports
+# ------------------------------------------------------------------------------
+from dart.viscous.Drivers.DGroupConstructor import GroupConstructor
+from dart.viscous.Drivers.DPostProcess import PostProcess
+
+from dart.viscous.Solvers.DIntegration import Integration
+from dart.viscous.Solvers.DSolver import Solver as TSolver
+from dart.viscous.Solvers.DSysMatrix import sysMatrix
+from dart.viscous.Solvers.DTransition import Transition
+
+from dart.viscous.Physics.DVariables import Variables
+from dart.viscous.Physics.DDiscretization import Discretization
+from dart.viscous.Physics.DValidation import Validation
+import dart.viscous.Physics.DBIConditions as Conditions
+
+from fwk.coloring import ccolors
+
+import numpy as np
+import matplotlib.pyplot as plt
+
+class Driver:
+ """Driver of the time-dependent boundary layer equation solver (viscous solver)
+ """
+
+ def __init__(self,_Re, user_xtr, _CFL0, _Minfty, _case, _mapMesh, _nFactor):
+ self.Minfty = _Minfty
+ self.Re=_Re
+ 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
+ self.case=_case
+ self.xtr=[1,1]
+ self.nbrElmUpper = -1
+ pass
+
+ def dictionary(self):
+ """Create dictionary with the coordinates and the boundary layer parameters
+ """
+
+ wData = {}
+ wData['x'] = self.data[:,0]
+ wData['y'] = self.data[:,1]
+ wData['z'] = self.data[:,2]
+ wData['x/c'] = self.data[:,0]/max(self.data[:,0])
+ wData['Cd'] = self.blScalars[0]
+ wData['Cdf'] = self.blScalars[1]
+ wData['Cdp'] = self.blScalars[2]
+ wData['delta*'] = self.blVectors[0]
+ wData['theta'] = self.blVectors[1]
+ wData['H'] = self.blVectors[2]
+ wData['Hk'] = self.blVectors[3]
+ wData['H*'] = self.blVectors[4]
+ wData['H**'] = self.blVectors[5]
+ wData['Cf'] = self.blVectors[6]
+ wData['Cdissip'] = self.blVectors[7]
+ wData['Ctau'] = self.blVectors[8]
+ wData['delta'] = self.blVectors[9]
+ if self.xtrF[0] is not None:
+ wData['xtrTop'] = self.xtrF[0]
+ else:
+ wData['xtrTop'] = self.xtr[0]
+ if self.xtrF[1] is not None:
+ wData['xtrBot'] = self.xtrF[1]
+ else:
+ wData['xtrBot'] = self.xtr[1]
+ self.wData=wData
+ return wData
+
+ def writeFile(self):
+ """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 = '')
+ 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('$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, 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 = self.nFactor
+
+ # Initialize computation structures.
+
+ 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()
+ 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.
+
+ 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
+ # Residuals and Jacobian computation.
+
+ transSolver = Transition(var.xtrF) # Transition solver.
+
+ # Initialize time Integrator
+
+ tIntegrator = Integration(sysSolver, self.Minfty, self.CFL0)
+ tSolver=TSolver(tIntegrator, transSolver, DirichletBC.wakeBC)
+
+ # Output setup at the first iteration.
+
+ 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('| - Number of elements: {:>5.0f}. {:>49s}'.format(var.nN,'|'))
+ if var.xtrF[0] is None:
+ print('| - Free transition on the upper side. {:>41}'.format('|'))
+ else:
+ print('| - Forced transition on the upper side; x/c = {:<.4f}. {:>25s}'.format(var.xtrF[0], '|'))
+ if var.xtrF[1] is None:
+ print('| - Free transition on the lower side. {:>41}'.format('|'))
+ else:
+ print('| - Forced transition on the lower side; x/c = {:<.4f}. {:>25s}'.format(var.xtrF[1], '|'))
+ print('| - Critical amplification ratio: {:<.2f}. {:>40s}'.format(var.Ncrit,'|'))
+
+ print('|{:>79}'.format('|'))
+ print('| Numerical parameters {:>57}'.format('|'))
+ print('| - CFL0: {:<3.2f} {:>65}'.format(tSolver.integrator.GetCFL0(),'|'))
+ print('| - Tolerance: {:<3.0f} {:>61}'.format(np.log10(tSolver.integrator.Gettol()),'|'))
+ print('| - Solver maximum iterations: {:<5.0f} {:>43}'.format(tSolver.integrator.GetmaxIt(),'|'))
+ print('|{:>79}'.format('|'))
+
+ # Output preprocessing satut.
+
+ print('+------------------------------ Preprocessing ---------------------------------+')
+
+ print('| - Maximum Mach number: {:<.2f}. {:>49s}'.format(np.max((var.extPar[0::var.nExtPar])),'|'))
+
+ # Initial condition.
+
+ if self.uPrevious is None: # Typically for the first iteration.
+
+ tSolver.initSln = 1 # Flag to use time solver initial condition.
+
+ tSolver.integrator.jacEval0 = 1 # Continuous Jacobian evaluation.
+ tSolver.integrator.SetCFL0(0.5) # Low starting CFL.
+
+ print('| - Time solver will provide the initial solution. {:>29}'.format('|'))
+
+ 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('|'))
+
+ DirichletBC.airfoilBC(var) # Reset boundary condition.
+
+ # Run the time integration
+ print('+------------------------------- Solver begin ---------------------------------+')
+ convergedPoints = tSolver.Run(var)
+ print('+-------------------------------- Solver Exit ---------------------------------+')
+
+ # Output time solver status.
+ if np.any(convergedPoints != 0):
+ print('|', ccolors.ANSI_YELLOW + 'Some point(s) did not converge.', ccolors.ANSI_RESET, '{:>45s}'.format('|'))
+
+ elif np.all(convergedPoints == 0):
+ print('|', ccolors.ANSI_GREEN + 'All points converged.', ccolors.ANSI_RESET, '{:>55s}'.format('|'))
+ print('|{:>79}'.format('|'))
+
+
+ # Save solution to use it as initial condition the next iteration.
+
+ var.u[4::var.nVar][var.flowRegime == 0] = 0. # Reset Ct in the laminar portion of the flow.
+ # (This quantity is not defined for a laminar flow)
+
+ self.uPrevious=var.u.copy() # Copy solution.
+
+ # Compute blowing velocities and sort the boundary layer parameters.
+
+ self.blScalars, blwVelUp, blwVelLw, blwVelWk, AFblVectors, WKblVectors, AFdeltaStarOut, WKdeltaStarOut = PostProcess().sortParam(var, self.mapMesh, cMeshDict, nFactor)
+
+ self.xtr = var.xtr.copy()
+ self.blVectors = AFblVectors
+
+ AFblwVel = np.concatenate((blwVelUp, blwVelLw))
+
+ # Group modification to ensure communication between the solvers.
+
+ groups[0].TEnd = [AFblVectors[1][var.bNodes[1]-1], AFblVectors[1][var.bNodes[2]-1]]
+ 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 = cMeshDict['locCoord'][:cMeshDict['bNodes'][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 = 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
diff --git a/dart/viscous/Drivers/VGroupConstructor.py b/dart/viscous/Drivers/DGroupConstructor.py
similarity index 99%
rename from dart/viscous/Drivers/VGroupConstructor.py
rename to dart/viscous/Drivers/DGroupConstructor.py
index 18f96c278070d4145cda3679ef89b422ad22acbc..c03b4827922bd520d79e9230c6c928d68e527c8e 100755
--- a/dart/viscous/Drivers/VGroupConstructor.py
+++ b/dart/viscous/Drivers/DGroupConstructor.py
@@ -14,10 +14,8 @@
# ------------------------------------------------------------------------------
# Imports
# ------------------------------------------------------------------------------
-import numpy as np
-import math
from matplotlib import pyplot as plt
-from scipy import interpolate
+import numpy as np
class GroupConstructor():
""" ---------- Group Constructor ----------
diff --git a/dart/viscous/Drivers/VPostProcess.py b/dart/viscous/Drivers/DPostProcess.py
similarity index 100%
rename from dart/viscous/Drivers/VPostProcess.py
rename to dart/viscous/Drivers/DPostProcess.py
diff --git a/dart/viscous/Drivers/VDriver.py b/dart/viscous/Drivers/VDriver.py
deleted file mode 100644
index 1e9a01784870797ae30a39d2993aa37006cacb48..0000000000000000000000000000000000000000
--- a/dart/viscous/Drivers/VDriver.py
+++ /dev/null
@@ -1,291 +0,0 @@
-################################################################################
-# #
-# DARTFlo viscous module file #
-# Driver : Communication with coupler and initialization #
-# of some components of the viscous computation part #
-# #
-# Author: Paul Dechamps #
-# Université de Liège #
-# Date: 2021.09.10 #
-# #
-################################################################################
-
-
-# ------------------------------------------------------------------------------
-# Imports
-# ------------------------------------------------------------------------------
-from dart.viscous.Drivers.VGroupConstructor import GroupConstructor
-from dart.viscous.Drivers.VPostProcess import PostProcess
-
-import dart.viscous.Solvers.VTimeSolver as TSolver
-from dart.viscous.Solvers.VSolver import sysMatrix
-from dart.viscous.Solvers.VTransition import Transition
-
-from dart.viscous.Physics.VVariables import Variables
-from dart.viscous.Physics.VDiscretization import Discretization
-from dart.viscous.Physics.VValidation import Validation
-import dart.viscous.Physics.VBIConditions as Conditions
-
-from fwk.coloring import ccolors
-
-import numpy as np
-import matplotlib.pyplot as plt
-
-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
- self.Cd=0
- self.it=0
- self.uPrevious=None
- self.blParPrevious=None
- self.xtrF=user_xtr
- self.mapMesh = _mapMesh
- self.Ti=None
- self.timeParams=_timeParams
- self.res=None
- self.upper=0
- self.case=_case
- self.xtr=[1,1]
- pass
-
- def dictionary(self):
- """Create dictionary with the coordinates and the boundary layer parameters
- """
-
- wData = {}
- wData['x'] = self.data[:,0]
- wData['y'] = self.data[:,1]
- wData['z'] = self.data[:,2]
- wData['x/c'] = self.data[:,0]/max(self.data[:,0])
- wData['Cd'] = self.blScalars[0]
- wData['Cdf'] = self.blScalars[1]
- wData['Cdp'] = self.blScalars[2]
- wData['delta*'] = self.blVectors[0]
- wData['theta'] = self.blVectors[1]
- wData['H'] = self.blVectors[2]
- wData['Hk'] = self.blVectors[3]
- wData['H*'] = self.blVectors[4]
- wData['H**'] = self.blVectors[5]
- wData['Cf'] = self.blVectors[6]
- wData['Cdissip'] = self.blVectors[7]
- wData['Ctau'] = self.blVectors[8]
- wData['delta'] = self.blVectors[9]
- if self.xtrF[0] is not None:
- wData['xtrTop'] = self.xtrF[0]
- else:
- wData['xtrTop'] = self.xtr[0]
- if self.xtrF[1] is not None:
- wData['xtrBot'] = self.xtrF[1]
- else:
- wData['xtrBot'] = self.xtr[1]
- self.wData=wData
- return wData
-
- def writeFile(self):
- """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 = '')
- 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('$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, 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(groups)):
- groups[n].connectListNodes, groups[n].connectListElems, dataIn[n] = groups[n].connectList()
- fMeshDict, cMeshDict, data = GroupConstructor().mergeVectors(dataIn, self.mapMesh, nFactor)
-
- # Initialize data container.
- var = Variables(fMeshDict, self.xtrF, self.Ti, self.Re)
-
- if self.it!=0:
- # Extenal flow local markers.
- var.extPar[4::var.nExtPar] = fMeshDict['xxExt']
-
- """if self.it >= 0:
- plt.plot(var.xGlobal[:var.bNodes[1]], var.extPar[2:var.bNodes[1]*var.nExtPar: var.nExtPar])
- plt.plot(cMeshDict['globCoord'][:cMeshDict['bNodes'][1]], cMeshDict['vtExt'][:cMeshDict['bNodes'][1]],'ob')
- plt.plot(var.xGlobal[var.bNodes[1]:var.bNodes[2]], var.extPar[var.bNodes[1]*var.nExtPar + 2:var.bNodes[2]*var.nExtPar: var.nExtPar])
- plt.plot(cMeshDict['globCoord'][cMeshDict['bNodes'][1]:cMeshDict['bNodes'][2]], cMeshDict['vtExt'][cMeshDict['bNodes'][1]:cMeshDict['bNodes'][2]],'or')
- plt.plot(var.xGlobal[var.bNodes[2]:], var.extPar[var.bNodes[2]*var.nExtPar + 2:: var.nExtPar])
- plt.plot(cMeshDict['globCoord'][cMeshDict['bNodes'][2]:], cMeshDict['vtExt'][cMeshDict['bNodes'][2]:],'og')
- plt.show()"""
-
- # Boundary condition (Airfoil/Wake).
- DirichletBC=Conditions.Boundaries(self.Re)
-
- # Gradient computation w/ different schemes.
- gradComp=Discretization(var.xx, var.extPar[4::var.nExtPar], var.bNodes) # Gradients computation.
-
- # Boundary layer equations solver that can provide Residuals and Jacobian computation.
- CSolver=sysMatrix(var.nN, var.nVar, gradComp.fou, gradComp.fou)
-
- # Transition solver.
- transSolver = Transition(var.xtrF)
-
- # Initialize time Integrator
- tSolver=TSolver.implicitEuler(CSolver, transSolver, DirichletBC.wakeBC, self.timeParams['CFL0'], self.Minfty)
-
- # Output setup at the first iteration.
-
- 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('| - Number of elements: {:>5.0f}. {:>49s}'.format(var.nN,'|'))
- if var.xtrF[0] is None:
- print('| - Free transition on the upper side. {:>41}'.format('|'))
- else:
- print('| - Forced transition on the upper side; x/c = {:<.4f}. {:>25s}'.format(var.xtrF[0], '|'))
- if var.xtrF[1] is None:
- print('| - Free transition on the lower side. {:>41}'.format('|'))
- else:
- print('| - Forced transition on the lower side; x/c = {:<.4f}. {:>25s}'.format(var.xtrF[1], '|'))
- print('| - Critical amplification ratio: {:<.2f}. {:>40s}'.format(var.Ncrit,'|'))
-
- # Output solver parameters.
- print('|{:>79}'.format('|'))
- print('| Numerical parameters {:>57}'.format('|'))
- print('| - CFL0: {:<3.2f} {:>65}'.format(tSolver.getCFL0(),'|'))
- print('| - Tolerance: {:<3.0f} {:>61}'.format(np.log10(tSolver.gettol()),'|'))
- print('| - Solver maximum iterations: {:<5.0f} {:>43}'.format(tSolver.getmaxIt(),'|'))
- print('|{:>79}'.format('|'))
-
- # Output preprocessing satut.
-
- print('+------------------------------ Preprocessing ---------------------------------+')
-
- print('| - Maximum Mach number: {:<.2f}. {:>49s}'.format(np.max((var.extPar[0::var.nExtPar])),'|'))
-
- # Initial condition: 'Driver' stores the solution form the previous coupling iteration.
-
- if self.uPrevious is None or self.it < 5:
- # Typically for the first iteration.
-
- # Initalize initial condition builder.
- initializer=Conditions.Initial(DirichletBC)
-
- # Compute initial condition based on Twaithes solution.
- initializer.initialConditions(var) # Set boundary condition
-
- tSolver.jacEval0 = 1
- tSolver.setCFL0(1)
- tSolver.initSln = 1
-
- #print('| - Using Twaithes solution as initial guess. {:>34}'.format('|'))
- print('| - Time solver will provide the initial solution. {:>29}'.format('|'))
-
- else:
- # If we use a solution previously obtained.
-
- 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('|'))
-
- # Run the time integration
- print('+------------------------------- Solver begin ---------------------------------+')
- convergedPoints = tSolver.timeSolve(var)
- print('+-------------------------------- Solver Exit ---------------------------------+')
-
- # Output time solver status.
- if np.any(convergedPoints != 1):
- print('|', ccolors.ANSI_YELLOW + 'Some points did not converge.', ccolors.ANSI_RESET, '{:>45s}'.format('|'))
-
- elif np.all(convergedPoints == 1):
- print('|', ccolors.ANSI_GREEN + 'All points converged.', ccolors.ANSI_RESET, '{:>55s}'.format('|'))
- print('|{:>79}'.format('|'))
-
-
- # Save solution to use it as initial condition the next iteration.
-
- # Reset Ct in the laminar portion of the flow.
- # (This quantity is not defined for a laminar flow)
- #var.u[4::var.nVar][var.flowRegime == 0] = 0.001
-
- # Copy solution.
- self.uPrevious=var.u.copy()
- self.blParPrevious=var.blPar.copy()
-
- # Compute blowing velocities and sort the boundary layer parameters.
- self.blScalars, blwVelUp, blwVelLw, blwVelWk, AFblVectors, WKblVectors, AFdeltaStarOut, WKdeltaStarOut = PostProcess().sortParam(var, self.mapMesh, cMeshDict, nFactor)
- self.xtr = var.xtr.copy()
- self.blVectors = AFblVectors
-
- AFblwVel = np.concatenate((blwVelUp, blwVelLw))
-
- # Group modification to ensure communication between the solvers.
-
- groups[0].TEnd = [AFblVectors[1][var.bNodes[1]-1], AFblVectors[1][var.bNodes[2]-1]]
- 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 = cMeshDict['locCoord'][:cMeshDict['bNodes'][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 = 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()]
-
- """if self.it >= 0:
- self.writeFile()
- Validation().main(1,self.wData, WKblVectors, var.xGlobal[var.bNodes[2]:])"""
\ No newline at end of file
diff --git a/dart/viscous/Master/VAirfoil.py b/dart/viscous/Master/DAirfoil.py
similarity index 99%
rename from dart/viscous/Master/VAirfoil.py
rename to dart/viscous/Master/DAirfoil.py
index 226091ddc57082c283e2dd2c3624d94af1de67ff..bf1431f50f62accb19598b66a6f82c4c0820dc69 100755
--- a/dart/viscous/Master/VAirfoil.py
+++ b/dart/viscous/Master/DAirfoil.py
@@ -19,7 +19,7 @@
## Airfoil around which the boundary layer is computed
# Amaury Bilocq
-from dart.viscous.Master.VBoundary import Boundary
+from dart.viscous.Master.DBoundary import Boundary
import numpy as np
diff --git a/dart/viscous/Master/VBoundary.py b/dart/viscous/Master/DBoundary.py
similarity index 100%
rename from dart/viscous/Master/VBoundary.py
rename to dart/viscous/Master/DBoundary.py
diff --git a/dart/viscous/Master/VCoupler.py b/dart/viscous/Master/DCoupler.py
similarity index 96%
rename from dart/viscous/Master/VCoupler.py
rename to dart/viscous/Master/DCoupler.py
index dcea5fb5054ce5d6929c411d860edcb3ea71f27a..da27a4830f9327a840890065f27a048dd52085ce 100755
--- a/dart/viscous/Master/VCoupler.py
+++ b/dart/viscous/Master/DCoupler.py
@@ -19,8 +19,8 @@
## Viscous-inviscid coupler (quasi-simultaneous coupling)
# Amaury Bilocq
-from dart.viscous.Master.VAirfoil import Airfoil
-from dart.viscous.Master.VWake import Wake
+from dart.viscous.Master.DAirfoil import Airfoil
+from dart.viscous.Master.DWake import Wake
import numpy as np
@@ -29,7 +29,7 @@ import dart.utils as floU
from fwk.coloring import ccolors
class Coupler:
- def __init__(self, _isolver, _vsolver, _boundaryAirfoil, _boundaryWake, _tol, _writer, _bnd, _timeParam=None):
+ def __init__(self, _isolver, _vsolver, _boundaryAirfoil, _boundaryWake, _tol, _writer, _bnd):
self.bnd=_bnd
self.isolver =_isolver # inviscid solver
self.vsolver = _vsolver # viscous solver
@@ -37,7 +37,6 @@ class Coupler:
self.group = [Airfoil(_boundaryAirfoil), Wake(_boundaryWake)] # airfoil and wake python objects
self.tol = _tol # tolerance of the VII
self.writer = _writer
- self.timeParam=_timeParam
def run(self):
"""Viscous inviscid coupling.
diff --git a/dart/viscous/Master/VWake.py b/dart/viscous/Master/DWake.py
similarity index 98%
rename from dart/viscous/Master/VWake.py
rename to dart/viscous/Master/DWake.py
index 3a6cfe47dfdda6020ac853f493355ea2ab835143..9ec5b93069843918b59bffe6ccd6e6b96ce025d2 100755
--- a/dart/viscous/Master/VWake.py
+++ b/dart/viscous/Master/DWake.py
@@ -19,7 +19,7 @@
## Wake behind airfoil (around which the boundary layer is computed)
# Amaury Bilocq
-from dart.viscous.Master.VBoundary import Boundary
+from dart.viscous.Master.DBoundary import Boundary
import numpy as np
diff --git a/dart/viscous/Physics/VBIConditions.py b/dart/viscous/Physics/DBIConditions.py
similarity index 99%
rename from dart/viscous/Physics/VBIConditions.py
rename to dart/viscous/Physics/DBIConditions.py
index 7f9efbb05b261d1529af20e7f927d570fce5e3b3..c1a145b9127aa156d64104d976e0f187ef75d598 100755
--- a/dart/viscous/Physics/VBIConditions.py
+++ b/dart/viscous/Physics/DBIConditions.py
@@ -14,7 +14,7 @@
# ------------------------------------------------------------------------------
# Imports
# ------------------------------------------------------------------------------
-from dart.viscous.Physics.VClosures import Closures
+from dart.viscous.Physics.DClosures import Closures
import matplotlib.pyplot as plt
import numpy as np
@@ -66,7 +66,7 @@ class Initial:
# Reynolds number based on the arc length.
Rex=rho*vt*xx/self.mu_air
Rex[Rex==0]=1 # Stagnation point
-
+ Rex[Rex<=0] = 1
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.
diff --git a/dart/viscous/Physics/VClosures.py b/dart/viscous/Physics/DClosures.py
similarity index 100%
rename from dart/viscous/Physics/VClosures.py
rename to dart/viscous/Physics/DClosures.py
diff --git a/dart/viscous/Physics/VDataStructures.py b/dart/viscous/Physics/DDataStructures.py
similarity index 100%
rename from dart/viscous/Physics/VDataStructures.py
rename to dart/viscous/Physics/DDataStructures.py
diff --git a/dart/viscous/Physics/VDiscretization.py b/dart/viscous/Physics/DDiscretization.py
similarity index 100%
rename from dart/viscous/Physics/VDiscretization.py
rename to dart/viscous/Physics/DDiscretization.py
diff --git a/dart/viscous/Physics/VValidation.py b/dart/viscous/Physics/DValidation.py
similarity index 100%
rename from dart/viscous/Physics/VValidation.py
rename to dart/viscous/Physics/DValidation.py
diff --git a/dart/viscous/Physics/VVariables.py b/dart/viscous/Physics/DVariables.py
similarity index 100%
rename from dart/viscous/Physics/VVariables.py
rename to dart/viscous/Physics/DVariables.py
diff --git a/dart/viscous/Solvers/DIntegration.py b/dart/viscous/Solvers/DIntegration.py
new file mode 100644
index 0000000000000000000000000000000000000000..f1d70c89c838c0ce009a34e1ceade46e92aae114
--- /dev/null
+++ b/dart/viscous/Solvers/DIntegration.py
@@ -0,0 +1,288 @@
+################################################################################
+# #
+# DARTFlo viscous module file #
+# Integration: Newton method for convergence #
+# to steady state. #
+# #
+# Author: Paul Dechamps #
+# Université de Liège #
+# Date: 2021.09.10 #
+# #
+################################################################################
+
+
+# ------------------------------------------------------------------------------
+# Imports
+# ------------------------------------------------------------------------------
+from dart.viscous.Physics.DClosures import Closures
+
+import math
+import numpy as np
+from fwk.coloring import ccolors
+
+
+class Integration:
+ """Pseudo-time integration of the integral boundary layer equations.
+ The non-linear equations are solved for one point using Newton's method
+
+ Attributes
+ ----------
+ - CFL0 (float): Starting CFL.
+ - maxIt (int): Maximum number of iterations.
+ - tol (float): Convergence criterion.
+ - itFrozenJac0 (int): Number of iterations during which the Jacobian matrix is frozen.
+ - Minf (float): Free-stream Mach number.
+ - gamma (float): Fluid heat capacity ratio.
+ """
+
+ def __init__(self, _sysMatrix, _Minf, _Cfl0) -> None:
+ self.sysMatrix = _sysMatrix
+
+ self.__CFL0=_Cfl0
+ self.__maxIt = 1000
+ self.__tol = 1e-4
+ self.itFrozenJac0 = 5
+
+ self.__Minf = _Minf
+ self.gamma = 1.4
+ pass
+
+ def GetCFL0(self): return self.__CFL0
+ def SetCFL0(self, _cfl0): self.__CFL0 = _cfl0
+
+ def Gettol(self): return self.__tol
+ def Settol(self, _tol): self.__tol = _tol
+
+ def GetmaxIt(self): return self.__maxIt
+ def SetmaxIt(self, _maxIt): self.__maxIt = _maxIt
+
+ def TimeMarching (self, DFlow, iMarker, screenOutput = 0) -> int:
+ """Pseudo-time marching solver.
+
+ Args
+ ----
+ - DFlow (class type): Data structure containing the flow configuration.
+ - iMarker (int): Id of the point to be treated.
+
+ Returns
+ -------
+ - Exit code info (int):
+ - 0 -> sucessfull exit.
+ - >0 -> convergence to tolerance not achieved, number of iterations.
+ - <0 -> illegal input or breakdown.
+
+ Opts
+ ----
+ - screenOutput (bool, optional): Output pseudo-time marching convergence.
+ """
+
+ # Save initial condition.
+
+ sln0 = DFlow.u.copy()
+
+ # 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.
+ iMarkerPrev = 0 if iMarker == bNodes[1] else iMarker - 1 # Point upstream of the considered point.
+ dx = DFlow.xx[iMarker] - DFlow.xx[iMarkerPrev] # Cell size.
+
+ # Initial residual.
+
+ sysRes = self.sysMatrix.buildResiduals(DFlow, iMarker)
+ normSysRes0 = np.linalg.norm(sysRes)
+
+ # 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.
+
+ # Main loop.
+
+ nErrors = 0 # Counter for the number of errors identified.
+ innerIters = 0 # Number of inner iterations to solver the non-linear system.
+ while innerIters <= self.__maxIt:
+
+ try:
+
+ # Jacobian computation.
+
+ if innerIters % itFrozenJac == 0:
+ Jacobian=self.sysMatrix.buildJacobian(DFlow, iMarker, sysRes, order = 1, eps = 1e-8)
+
+ # Linear solver.
+
+ slnIncr = np.linalg.solve((Jacobian + IMatrix), - sysRes)
+
+ # Increment solution.
+
+ DFlow.u[iMarker*nVar : (iMarker+1)*nVar] += slnIncr
+
+ # Residual computation.
+
+ 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)
+ return 0
+
+ except KeyboardInterrupt:
+ print(ccolors.ANSI_RED + 'Program terminated by user.', ccolors.ANSI_RESET)
+ quit()
+ except Exception as error:
+
+ nErrors += 1
+
+ if nErrors >= 5:
+ print('|', ccolors.ANSI_RED + 'Too many errors identified at position {:<1.4f}. Marker: {:<4.0f}. {:>17s}'
+ .format(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
+ 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
+
+
+ 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
+
+
+ def __GetSoundSpeed(self, invVelSquared) -> float:
+ """Compute the speed of sound in a cell.
+
+ Args
+ ----
+ - invVelSquared (float): Square of the inviscid velocity.
+ """
+
+ return np.sqrt(1 / (self.__Minf * self.__Minf) + 0.5 * (self.gamma - 1) * (1 - invVelSquared))
+
+ def __SetTimeStep(self, CFL, dx, invVel) -> float:
+ """Compute the pseudo-time step in a cell for a given
+ CFL number.
+
+ Args
+ ----
+ - CFL (float): Current CFL number.
+ - dx (float): Cell size.
+ - invVel (float): Inviscid velocity in the cell.
+
+ Returns
+ -------
+ - timeStep (float): Pseudo-time step.
+ """
+
+ # Speed of sound.
+ gradU2 = (invVel)**2
+ # Velocity of the fastest travelling wave
+ soundSpeed = self.__GetSoundSpeed(gradU2)
+ advVel = soundSpeed + invVel
+
+ # Time step computation
+ return CFL*dx/advVel
+
+ def __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
+ nBlPar = DFlow.nBlPar
+ nExtPar = DFlow.nExtPar
+
+ DFlow.u[iMarker*nVar:(iMarker+1)*nVar] = sln0[(iMarker)*nVar: (iMarker+1)*nVar].copy()
+ name = 'airfoil' if iMarker <= DFlow.bNodes[2] else 'wake'
+
+ if DFlow.flowRegime[iMarker]==0:
+
+ # Laminar closure relations
+ DFlow.blPar[iMarker*nBlPar+1:iMarker*nBlPar+7] = Closures.laminarClosure(DFlow.u[iMarker*nVar + 0], DFlow.u[iMarker*nVar + 1],
+ DFlow.blPar[iMarker*nBlPar+0], DFlow.extPar[iMarker*nExtPar+0],
+ name)
+
+ elif DFlow.flowRegime[iMarker]==1:
+
+ # Turbulent closures relations
+ DFlow.blPar[iMarker*nBlPar+1:iMarker*nBlPar+9] = Closures.turbulentClosure(DFlow.u[iMarker*nVar + 0], DFlow.u[iMarker*nVar + 1],
+ DFlow.u[iMarker*nVar + 4], DFlow.blPar[iMarker*nBlPar+0],
+ DFlow.extPar[iMarker*nExtPar+0], name)
diff --git a/dart/viscous/Solvers/DSolver.py b/dart/viscous/Solvers/DSolver.py
new file mode 100644
index 0000000000000000000000000000000000000000..39865e30003a474c964025134021d1ba6cb1b352
--- /dev/null
+++ b/dart/viscous/Solvers/DSolver.py
@@ -0,0 +1,167 @@
+################################################################################
+# #
+# DARTFlo viscous module file #
+# Time solver : Time discretization #
+# and computation towards steady state #
+# #
+# Author: Paul Dechamps #
+# Université de Liège #
+# Date: 2021.09.10 #
+# #
+################################################################################
+
+
+# ------------------------------------------------------------------------------
+# Imports
+# ------------------------------------------------------------------------------
+
+from dart.viscous.Physics.DClosures import Closures
+from fwk.coloring import ccolors
+import numpy as np
+
+class Solver:
+ """Performs the downstream marching process. Calls the pseudo-time integration based
+ on Newton's method on demand.
+
+ Attributes
+ ----------
+
+ - integrator (class type): Module able to perform pseudo-time integration of one point.
+ - transSolver (class type): Module able to treat the transition related computation.
+ - wakeBC (class type): Wake boundary condition.
+ - initSln (bool): Flag to initialize the initial condition or to use the value in input.
+ """
+
+ def __init__(self, _integrator, transition, wakeBC) -> None:
+
+ # Initialize time parameters
+
+ self.integrator = _integrator
+ self.transSolver = transition
+ self.wakeBC = wakeBC
+ self.initSln = 0
+
+ def Run(self, DFlow):
+ """Runs the downstream marching solver.
+ """
+
+ # Retreive parameters.
+
+ 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[0] = 0 # Boundary condition.
+
+ self.transSolver.initTransition(DFlow) # Initialize the flow regime on each cell.
+
+ lockTrans = 0 # Flag to remember if the transition is already found or not.
+
+ for iMarker in range(1, nMarkers):
+ displayState = 0 # Flag to output simulation state.
+
+ 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 iMarker == bNodes[0] + 1: # Upper side start.
+ print('| - Computing upper side. {:>54}'.format('|'))
+
+ if iMarker == bNodes[1]: # Lower side start.
+ lockTrans = 0
+ print('| - Computing lower side. {:>54}'.format('|'))
+
+ if iMarker == bNodes[2]: # First wake point
+
+ self.wakeBC(DFlow) # Wake boundary condition.
+ convergedPoints[iMarker] = 0
+ lockTrans = 1
+ print('| - Computing wake. {:>60}'.format('|'))
+
+ continue # Ignore remaining instructions for the first wake point.
+
+ # Call Newton solver for the point.
+
+ convergedPoints[iMarker] = self.integrator.TimeMarching(DFlow, iMarker, displayState)
+
+ # Check for transition.
+
+ if not lockTrans:
+
+ self.__CheckWaves(DFlow, iMarker)
+
+ # Free transition.
+ if DFlow.u[iMarker * DFlow.nVar + 2] >= DFlow.Ncrit:
+ print('| Transition located near x/c = {:<2.3f}. Marker: {:<4.0f} {:>28s}'.format(DFlow.xGlobal[iMarker],iMarker, '|'))
+ self.__AverageTransition(DFlow, iMarker, displayState)
+ 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]:
+ print('| Forced transition near x/c = {:<2.3f}. Marker: {:<4.0f} {:>28s}'.format(DFlow.xGlobal[iMarker],iMarker, '|'))
+ self.__AverageTransition(DFlow, iMarker)
+ lockTrans = 1
+
+ return convergedPoints
+
+ def __CheckWaves(self, var, iMarker) -> None:
+ """Determine if the amplification waves start growing or not.
+ """
+
+ logRet_crit = 2.492*(1/(var.blPar[(iMarker)*var.nBlPar + 6]-1))**0.43
+ + 0.7*(np.tanh(14*(1/(var.blPar[(iMarker)*var.nBlPar + 6]-1))-9.24)+1.0)
+
+ if var.blPar[(iMarker)*var.nBlPar + 0] > 0: # Avoid log of negative number.
+ if np.log10(var.blPar[(iMarker)*var.nBlPar + 0]) <= logRet_crit - 0.08:
+ var.u[iMarker*var.nVar + 2] = 0
+
+ pass
+
+
+ def __AverageTransition(self, var, iMarker, displayState) -> None:
+ """Averages laminar and turbulent solution at the transition location.
+ The boundary condition of the turbulent flow portion is also imposed.
+ """
+
+ # Save laminar solution.
+ lamSol = var.u[(iMarker)*var.nVar : (iMarker+1)* var.nVar].copy()
+
+ # Set up turbulent portion boundary condition.
+ avTurb = self.transSolver.solveTransition(var, iMarker)
+
+ # Compute turbulent solution.
+ var.flowRegime[iMarker] = 1
+ iMarkerPrev = 0 if iMarker == var.bNodes[1] else iMarker - 1
+
+ # Avoid starting with 0 for the shear stress and high H.
+ var.u[iMarker*var.nVar + 4] = var.u[iMarkerPrev*var.nVar + 4]
+ var.u[iMarker*var.nVar + 1] = 1.515
+
+ self.integrator.TimeMarching(var, iMarker, displayState)
+
+ # Average solutions.
+ avLam = 1 - avTurb
+ thetaTrans = avLam * lamSol[0] + avTurb * var.u[(iMarker)*var.nVar + 0]
+ HTrans = avLam * lamSol[1] + avTurb * var.u[(iMarker)*var.nVar + 1]
+ nTrans = 0
+ ueTrans = avLam * lamSol[3] + avTurb * var.u[(iMarker)*var.nVar + 3]
+ CtTrans = avTurb * var.u[(iMarker)*var.nVar + 4]
+ #CtTrans = max(CtTrans, var.u[(iMarkerPrev)*var.nVar + 4])
+ var.u[(iMarker)*var.nVar : (iMarker+1)*var.nVar] = [thetaTrans, HTrans, nTrans, ueTrans, CtTrans]
+
+ # Recompute closures.
+ var.blPar[(iMarker)*var.nBlPar+1:(iMarker)*var.nBlPar+9]=Closures.turbulentClosure(var.u[(iMarker)*var.nVar+0],
+ var.u[(iMarker)*var.nVar+1],
+ var.u[(iMarker)*var.nVar+4],
+ var.blPar[(iMarker)*var.nBlPar+0],
+ var.extPar[(iMarker)*var.nExtPar+0],
+ 'airfoil')
+ var.blPar[(iMarker-1)*var.nBlPar+1:(iMarker-1)*var.nBlPar+7]=Closures.laminarClosure(var.u[(iMarker-1)*var.nVar+0],
+ var.u[(iMarker-1)*var.nVar+1],
+ var.blPar[(iMarker-1)*var.nBlPar+0],
+ var.extPar[(iMarker-1)*var.nExtPar+0],
+ 'airfoil')
+ pass
\ No newline at end of file
diff --git a/dart/viscous/Solvers/VSolver.py b/dart/viscous/Solvers/DSysMatrix.py
similarity index 81%
rename from dart/viscous/Solvers/VSolver.py
rename to dart/viscous/Solvers/DSysMatrix.py
index 17f16c913c342276968169b42a664f271217ddbf..0e55194bb71f1f06bee8641a17eb70b51dd08c1b 100644
--- a/dart/viscous/Solvers/VSolver.py
+++ b/dart/viscous/Solvers/DSysMatrix.py
@@ -14,7 +14,7 @@
# ------------------------------------------------------------------------------
# Imports
# ------------------------------------------------------------------------------
-from dart.viscous.Physics.VClosures import Closures
+from dart.viscous.Physics.DClosures import Closures
import numpy as np
@@ -96,37 +96,38 @@ class flowEquations:
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
- #
- # -------------------------------------------------------------------------------------------------------------------------------
+ # +------------------------------------------- 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
@@ -141,7 +142,7 @@ class flowEquations:
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
+ elif dataCont.flowRegime[iMarker] == 0: spaceMatrix[4] = 0
# Temporal part
timeMatrix[0] = [Ha/uea, Ta/uea, 0, Ta * Ha/(uea**2), 0]
diff --git a/dart/viscous/Solvers/VTransition.py b/dart/viscous/Solvers/DTransition.py
similarity index 99%
rename from dart/viscous/Solvers/VTransition.py
rename to dart/viscous/Solvers/DTransition.py
index 461f8907be33c373ed227c91695d5c9d8e1b9291..e59b9a9598e9e3817a06340c83ae26d202224987 100644
--- a/dart/viscous/Solvers/VTransition.py
+++ b/dart/viscous/Solvers/DTransition.py
@@ -13,7 +13,7 @@
# ------------------------------------------------------------------------------
# Imports
# ------------------------------------------------------------------------------
-from dart.viscous.Physics.VClosures import Closures
+from dart.viscous.Physics.DClosures import Closures
import numpy as np
diff --git a/dart/viscous/Solvers/VTimeSolver.py b/dart/viscous/Solvers/VTimeSolver.py
deleted file mode 100644
index 8643dcf22e03899b4c255b00b9fca119064bed54..0000000000000000000000000000000000000000
--- a/dart/viscous/Solvers/VTimeSolver.py
+++ /dev/null
@@ -1,497 +0,0 @@
-################################################################################
-# #
-# DARTFlo viscous module file #
-# Time solver : Time discretization #
-# and computation towards steady state #
-# #
-# Author: Paul Dechamps #
-# Université de Liège #
-# Date: 2021.09.10 #
-# #
-################################################################################
-
-
-# ------------------------------------------------------------------------------
-# Imports
-# ------------------------------------------------------------------------------
-from dart.viscous.Physics.VClosures import Closures
-
-from matplotlib import pyplot as plt
-from fwk.coloring import ccolors
-import scipy.sparse.linalg as linSolver
-import numpy as np
-import math
-
-class timeConfig:
- """ ---------- Time integration configuration ----------
-
- Time integration related computations: CFL adaptation, time step (per cell)
-
- Attributes
- ----------
- cflAdapt: bool
- 'True' if CFL adaptation is active, 'False' otherwise
-
- __cflAdaptParam: numpy vect
- [Lower, Upper] bounds of the CFL
- """
- def __init__(self, _Minfty):
-
- # Correct incompressible Mach number.
- if _Minfty != 0: self.__Minf = _Minfty
- else: self.__Minf = 0.1
-
- self.gamma = 1.4
- pass
-
- def getSoundSpeed(self, gradU2):
- return np.sqrt(1 / (self.__Minf * self.__Minf) + 0.5 * (self.gamma - 1) * (1 - gradU2))
-
- def computeTimeStep(self, CFL, dx, invVel):
-
- # Speed of sound.
- gradU2 = (invVel)**2
- # Velocity of the fastest travelling wave
- soundSpeed = self.getSoundSpeed(gradU2)
- advVel = soundSpeed + invVel
-
- # Time step computation
- return CFL*dx/advVel
-
-
- def adaptCFL(self, resCurr, res0):
- """CFL number adaptation. The metric driving the adaptation is the residual value.
- """
- CFL=self.__CFL0*(res0/resCurr)**0.7
- CFL=max(CFL,self.__CFL0)
- return CFL
-
- def outputState(self, var, iMarker, innerIters, normLinSysRes, normLinSysRes0, CFL, color='white'):
- if color == 'white':
- print('| {:>6.0f}| {:>11.0f}| {:>9.4f}| {:>10.2f}| {:>14.6f}| {:>6.0f}| {:>11.3f}|'
- .format(iMarker,
- innerIters,
- np.log10(normLinSysRes/normLinSysRes0),
- CFL,
- var.xGlobal[iMarker],
- var.flowRegime[iMarker],
- var.u[iMarker*var.nVar + 2]))
- elif color == 'yellow':
- print(ccolors.ANSI_YELLOW + '| {:>6.0f}| {:>11.0f}| {:>9.4f}| {:>8.2f}| {:>14.6f}| {:>6.0f}| {:>11.3f}|'
- .format(iMarker,
- innerIters,
- np.log10(normLinSysRes/normLinSysRes0),
- CFL,
- var.xGlobal[iMarker],
- var.flowRegime[iMarker],
- var.u[iMarker*var.nVar + 2]),
- ccolors.ANSI_RESET)
- elif color == 'green':
- print(ccolors.ANSI_YELLOW + '| {:>6.0f}| {:>11.0f}| {:>9.4f}| {:>8.2f}| {:>14.6f}| {:>6.0f}| {:>11.3f}|'
- .format(iMarker,
- innerIters,
- np.log10(normLinSysRes/normLinSysRes0),
- CFL,
- var.xGlobal[iMarker],
- var.flowRegime[iMarker],
- var.u[iMarker*var.nVar + 2]),
- ccolors.ANSI_RESET)
-
-
-
-class implicitEuler(timeConfig):
- """ ---------- Implicit Euler time integration ----------
-
- Solves the unsteady boundary layer equations using the implicit Euler
- method with Newton method to solve the non-linear system.
-
- Attributes
- ----------
- solver: class
- Spacial operator and Jacobian computation of the boundary layer laws
-
- timeParams: dict
- Contains necessary information for the solver: tolerance, CFL, CFL apdation parameters
-
- safeguard: int
- Number of iterations for which the residual has to decrease to unlock CFL adaptation
- (usually changes after the first coupling iteration)
- and second order Jacobian evaluation starts.
-
- safemode: bool
- Indicates if the solver is in safe mode (1) or not (0). Value switched according to safeguard
- """
-
- def __init__(self, _solver, transition, wakeBC, _cfl0, _Minfty):
- timeConfig.__init__(self, _Minfty)
- # Initialize time parameters
-
- self.__CFL0=_cfl0
- self.__maxIt = 15000
- self.__tol = 1e-8
- self.jacEval0 = 5
- self.initSln = 0
-
- self.solver=_solver
- self.transSolver = transition
- self.wakeBC = wakeBC
- def __str__(self):
- return 'Implicit Euler'
-
- def getCFL0(self): return self.__CFL0
- def setCFL0(self, _cfl0): self.__CFL0 = _cfl0
-
- def gettol(self): return self.__tol
- def settol(self, _tol): self.__tol = _tol
-
- def getmaxIt(self): return self.__maxIt
- def setmaxIt(self, _maxIt): self.__maxIt = _maxIt
-
- def resetSolution(self, iMarker, var, u0):
-
- if iMarker == var.bNodes[1]:
-
- var.u[iMarker*var.nVar : (iMarker + 1)*var.nVar] = var.u[0*var.nVar : 1*var.nVar]
-
- elif iMarker != var.transMarkers[0] and iMarker != var.transMarkers[1]:
-
- var.u[iMarker*var.nVar : (iMarker + 1)*var.nVar] = var.u[(iMarker-1)*var.nVar : (iMarker)*var.nVar]
-
- else:
-
- var.u[iMarker*var.nVar : (iMarker + 1)*var.nVar] = u0[iMarker*var.nVar : (iMarker + 1)*var.nVar].copy()
-
- def timeSolve(self, var):
- """ ---------- Newton iteration to solve the boundary layer equations ----------
-
- Damped Newton method in Pseudo-Time. Linear algebraic system (obtained by Taylor linearis.)
- is solved with a direct method.
-
- Equations
- ---------
- ∂U/∂t = F(U)
- <=> (U_(n+1)-U_n)/∆t = - F(U_(n+1)) (Time disc.)
- <=> ∆U/∆t = -(F(U_n) + ∂F/∂U*∆u), where ∂F/∂U = J(U_n)
- <=> (I+J)dU = -F(U_n), J=dF_n/dU_n (Linearis.)
-
- Parameters
- ----------
- - var: class
- Data container: - Solution
- - Boundary layer parameters
- - External flow parameters
- - Mesh
- - Transition information
-
- Tasks
- ------
- - Updates the solution 'var.u' through Newton iterations until steady state.
-
- - Asks 'solver' for spacial operator F and Jacobian J and uses a
- direct linear solver (LU decomposition).
-
- - Corrects undesirable results.
-
- - Asks the solver for transtion position/BC and wake BC to be updated at each iteration.
- """
-
-
- nN = var.nN # Number of points in the computational domain.
- nBlPar = var.nBlPar # Number of boundary layer parameters used for vector indexing.
- nExtPar = var.nExtPar # Number of boundary layer parameters used for vector indexing.
- bNodes = var.bNodes # Boundary nodes of the computational domain.
- convergedPoints = np.zeros(var.nN) # Convergence control of each point.
- convergedPoints[0] = 1 # Boundary condition.
-
- self.transSolver.initTransition(var) # Initialize the flow regime on each cell.
-
- lockTrans = 0 # Flag to remember if the transition is already found or not.
-
- for iMarker in range(1, nN):
- displayState = 0 # Flag to output simulation state.
-
- if self.initSln:
- iMarkerPrev = 0 if iMarker == var.bNodes[1] else iMarker - 1
- var.u[iMarker*var.nVar: (iMarker + 1)*var.nVar] = var.u[iMarkerPrev*var.nVar: (iMarkerPrev + 1)* var.nVar]
-
- if not lockTrans and 1 < iMarker < bNodes[2]:
-
- self.__checkWaves(var, iMarker)
-
- # Upper side start.
- if iMarker == bNodes[0]+1:
- print('| - Computing upper side. {:>54}'.format('|'))
-
- # Lower side start.
- if iMarker == bNodes[1]:
- lockTrans = 0
- print('| - Computing lower side. {:>54}'.format('|'))
-
- # Wake start
- if iMarker == bNodes[2]: # First wake point
-
- # Wake boundary condition.
- self.wakeBC(var)
- convergedPoints[iMarker] = 1
- lockTrans = 1
- print('| - Computing wake. {:>60}'.format('|'))
-
- # Ignore remaining instructions for the first wake point.
- continue
-
- # Call Newton solver for the point.
- convergedPoints[iMarker] = self.newtonSolver(var, iMarker, displayState)
-
- # Check for transition.
- if not lockTrans:
-
- # Free transition.
- if var.u[iMarker * var.nVar + 2] >= var.Ncrit:
- print('| Transition located near x/c = {:<2.3f}. Marker: {:<4.0f} {:>28s}'.format(var.xGlobal[iMarker],iMarker, '|'))
- self.__averageTransition(var, iMarker, displayState)
- lockTrans = 1
-
- # Forced transition.
- elif var.xtrF[0] is not None or var.xtrF[1] is not None:
- if var.flowRegime[iMarker] == 1 and var.flowRegime[iMarker - 1]:
- print('| Forced transition near x/c = {:<2.3f}. Marker: {:<4.0f} {:>28s}'.format(var.xGlobal[iMarker],iMarker, '|'))
- self.__forcedTransition(var, iMarker)
- lockTrans = 1
-
- return convergedPoints
-
-
- def newtonSolver (self, var, iMarker, displayState):
-
- # Save initial condition in case a point diverges.
- u0=var.u.copy() # Initial solution.
-
- nVar = var.nVar # Number of variables of the differential problem.
- nBlPar = var.nBlPar # Number of boundary layer parameters used for vector indexing.
- nExtPar = var.nExtPar # Number of boundary layer parameters used for vector indexing.
- bNodes = var.bNodes # Boundary nodes of the computational domain.
- iInvVel = var.extPar[iMarker*nExtPar + 2] # Inviscid velocity on the considered cell.
- iMarkerPrev = 0 if iMarker == bNodes[1] else iMarker - 1 # Point upstream of the considered point.
- dx = var.xx[iMarker] - var.xx[iMarkerPrev] # Cell size.
-
- sysRes = self.solver.buildResiduals(var, iMarker) # System initial residuals.
- normSysRes0 = np.linalg.norm(sysRes) # Initial norm of the residuals.
-
- CFL = self.__CFL0 # Initialized CFL number.
- CFLAdapt = 1 # Flag for CFL adaptation.
- jacEval = self.jacEval0 # Number of iterations for which Jacobian is frozen.
- dt = self.computeTimeStep(CFL, dx, iInvVel) # Initial time step.
- IMatrix=np.identity(nVar) / dt # Identity matrix used to construct the linear system.
-
- # Main loop.
- nErrors = 0 # Counter for the number of errors identified at that point.
- innerIters = 0 # Number of inner iterations to solver the non-linear system.
- while innerIters <= self.__maxIt:
-
- try:
-
- # Jacobian computation.
-
- if innerIters % jacEval == 0:
- Jacobian=self.solver.buildJacobian(var, iMarker, sysRes, order = 1, eps = 1e-8)
-
- # Linear solver.
-
- slnIncr = np.linalg.solve((Jacobian + IMatrix), - sysRes)
- #slnIncr, exitCode = linSolver.gmres((Jacobian + IMatrix), -sysRes, tol = 1e-10)
-
- # Increment solution and correct absurd transcient results.
-
- var.u[iMarker*nVar : (iMarker+1)*nVar] += slnIncr
- # var.u[iMarker*nVar + 1] = max(var.u[iMarker*nVar + 1], 1.0005)
-
- # Compute Residuals.
- sysRes = self.solver.buildResiduals(var, iMarker)
- normSysRes = np.linalg.norm(slnIncr/dt - sysRes)
-
- if displayState and innerIters > 0 and innerIters % 1000 == 0:
- self.outputState(var, iMarker, innerIters, normSysRes, normSysRes0, CFL, 'yellow')
- # Check convergence.
- if normSysRes <= self.__tol * normSysRes0:
- converged = 1
- if displayState:
- self.outputState(var, iMarker, innerIters, normSysRes, normSysRes0, CFL)
- return 1
-
- except KeyboardInterrupt:
- quit()
- except Exception as error:
-
- nErrors += 1
- if nErrors >= 5:
- print('|', ccolors.ANSI_RED + 'Too many errors identified at position {:<1.4f}. Marker: {:<4.0f}. {:>17s}'
- .format(var.xGlobal[iMarker], iMarker,'|'), ccolors.ANSI_RESET)
- quit()
-
- print('|', ccolors.ANSI_YELLOW + 'Newton solver failed at position {:<.4f}. Marker: {:<.0f} {:>25s}'
- .format(var.xGlobal[iMarker], iMarker,'|'), ccolors.ANSI_RESET)
-
- print(error)
-
- #self.resetSolution(iMarker, var, u0)
- var.u[iMarker*nVar : (iMarker+1)*nVar] = var.u[(iMarker-1)*nVar : (iMarker)*nVar]
-
- # Lower CFL and recompute time step
- print('Current CFL', CFL)
- if math.isnan(CFL):
- CFL = self.__CFL0
- CFL = min(max(CFL / 2,0.01),1)
- print('Lowering CFL: {:<.3f}'.format(CFL))
- CFLAdapt = 1
- dt = self.computeTimeStep(CFL, dx, iInvVel)
- IMatrix = np.identity(nVar) / dt
-
- sysRes = self.solver.buildResiduals(var, iMarker)
-
- jacEval = 1
- displayState = 1
- print('+------------------------------------------------------------------------------+')
- print('| {:>6s}| {:>8s}| {:>9s}| {:>8s}| {:>12s}| {:>6s}| {:>8s}|'.format('Point','Inner Iters',
- 'rms[F]', 'End CFL',
- 'Position (x/c)', 'Regime',
- 'Amp. factor'))
- print('+------------------------------------------------------------------------------+')
- continue
-
- # CFL adaptation
- if CFLAdapt:
-
- CFL = self.__CFL0* (normSysRes0/normSysRes)**0.7
-
- CFL = max(CFL, 0.01)
- dt = self.computeTimeStep(CFL, dx, iInvVel)
- IMatrix=np.identity(nVar) / dt
-
-
- innerIters+=1
-
- else:
- # Maximum number of iterations reached.
- if normSysRes >= 1e-3 * normSysRes0:
- print('|', ccolors.ANSI_RED + 'Too many iterations at position {:<.4f}. Marker: {:>5.0f}. normLinSysRes = {:<.3f}. {:>30s}'
- .format(var.xGlobal[iMarker], iMarker,
- np.log10(normSysRes/normSysRes0),'|'),
- ccolors.ANSI_RESET)
- var.u[iMarker*nVar:(iMarker+1)*nVar] = u0[(iMarker)*nVar: (iMarker+1)*nVar].copy()
- name = 'airfoil' if iMarker <= var.bNodes[2] else 'wake'
-
- if var.flowRegime[iMarker]==0:
-
- # Laminar closure relations
- var.blPar[iMarker*nBlPar+1:iMarker*nBlPar+7] = Closures.laminarClosure(var.u[iMarker*nVar + 0], var.u[iMarker*nVar + 1],
- var.blPar[iMarker*nBlPar+0], var.extPar[iMarker*nExtPar+0],
- name)
-
- elif var.flowRegime[iMarker]==1:
-
- # Turbulent closures relations
- var.blPar[iMarker*nBlPar+1:iMarker*nBlPar+9] = Closures.turbulentClosure(var.u[iMarker*nVar + 0], var.u[iMarker*nVar + 1],
- var.u[iMarker*nVar + 4], var.blPar[iMarker*nBlPar+0],
- var.extPar[iMarker*nExtPar+0], name)
-
- else:
- print('|', ccolors.ANSI_YELLOW+ 'Too many iterations at position {:<.4f}. Marker: {:>5.0f} normLinSysRes = {:<.3f}. {:>30s}'
- .format(var.xGlobal[iMarker], iMarker,
- np.log10(normSysRes/normSysRes0),'|'),
- ccolors.ANSI_RESET)
- return 0
-
- def __checkWaves(self, var, iMarker):
- """Determine if the amplification waves start growing or not.
- """
-
- logRet_crit = 2.492*(1/(var.blPar[(iMarker-1)*var.nBlPar + 6]-1))**0.43
- + 0.7*(np.tanh(14*(1/(var.blPar[(iMarker-1)*var.nBlPar + 6]-1))-9.24)+1.0)
-
- if var.blPar[(iMarker-1)*var.nBlPar + 0] > 0: # Avoid log of something <= 0
- if np.log10(var.blPar[(iMarker-1)*var.nBlPar + 0]) <= logRet_crit - 0.08:
- var.u[iMarker*var.nVar + 2] = 0
-
- pass
-
-
- def __averageTransition(self, var, iMarker, displayState):
- """Averages laminar and turbulent solution at the transition location.
- The boundary condition of the turbulent flow portion is also imposed.
- """
-
- # Save laminar solution
- lamSol = var.u[(iMarker)*var.nVar : (iMarker+1)* var.nVar].copy()
-
- # Set up turbulent portion boundary condition
- avTurb = self.transSolver.solveTransition(var, iMarker)
-
- # Compute turbulent solution
- var.flowRegime[iMarker] = 1
- iMarkerPrev = 0 if iMarker == var.bNodes[1] else iMarker - 1
- # Avoid starting with 0 for the shear stress.
- var.u[iMarker*var.nVar + 4] = var.u[iMarkerPrev*var.nVar + 4]
-
- self.newtonSolver(var, iMarker, displayState)
-
- # Average solutions
- avLam = 1 - avTurb
- thetaTrans = avLam * lamSol[0] + avTurb * var.u[(iMarker)*var.nVar + 0]
- HTrans = avLam * lamSol[1] + avTurb * var.u[(iMarker)*var.nVar + 1]
- nTrans = 0
- ueTrans = avLam * lamSol[3] + avTurb * var.u[(iMarker)*var.nVar + 3]
- CtTrans = 0 * lamSol[4] + avTurb * var.u[(iMarker)*var.nVar + 4]
- CtTrans = max(CtTrans, 0.03)
- var.u[(iMarker)*var.nVar : (iMarker+1)*var.nVar] = [thetaTrans, HTrans, nTrans, ueTrans, CtTrans]
-
- # Recompute closures
- var.blPar[(iMarker)*var.nBlPar+1:(iMarker)*var.nBlPar+9]=Closures.turbulentClosure(var.u[(iMarker)*var.nVar+0],
- var.u[(iMarker)*var.nVar+1],
- var.u[(iMarker)*var.nVar+4],
- var.blPar[(iMarker)*var.nBlPar+0],
- var.extPar[(iMarker)*var.nExtPar+0],
- 'airfoil')
- var.blPar[(iMarker-1)*var.nBlPar+1:(iMarker-1)*var.nBlPar+7]=Closures.laminarClosure(var.u[(iMarker-1)*var.nVar+0],
- var.u[(iMarker-1)*var.nVar+1],
- var.blPar[(iMarker-1)*var.nBlPar+0],
- var.extPar[(iMarker-1)*var.nExtPar+0],
- 'airfoil')
- pass
-
- def __forcedTransition(self, var, iMarker):
- # Forced transition
- """if not lockTrans and var.xtrF[0] is not None and (var.flowRegime[iMarker] == 1 and var.flowRegime[iMarker-1] == 0):
- print('| Transition near {:<1.4f}. Marker: {:<5.0f} {:>40}'.format(var.xGlobal[iMarker-1],iMarker-1, '|'))
- lamSol = var.u[(iMarker-1)*var.nVar : (iMarker)* var.nVar].copy()
- avTurb = self.transSolver.transitionBC(var, iMarker - 1)
- # Compute turbulent solution
- var.flowRegime[iMarker - 1] = 1
- self.newtonSolver(var, iMarker-1, displayState)
- # Average solutions
- var.u[(iMarker-1)*var.nVar : (iMarker)*var.nVar] = (1-avTurb) * lamSol + avTurb * var.u[(iMarker-1)*var.nVar : (iMarker)*var.nVar]
- # Recompute closures
- var.blPar[(iMarker-1)*nBlPar+1:(iMarker-1)*nBlPar+9]=Closures.turbulentClosure(var.u[(iMarker-1)*var.nVar+0],
- var.u[(iMarker-1)*var.nVar+1],
- var.u[(iMarker-1)*var.nVar+4],
- var.blPar[(iMarker-1)*nBlPar+0],
- var.extPar[(iMarker-1)*nExtPar+0],
- 'airfoil')
-
- lockTrans = 1
- if not lockTrans and var.xtrF[1] is not None and (var.flowRegime[iMarker]== 1 and var.flowRegime[iMarker-1] == 0):
- print('| Transition near {:<1.4f}. Marker: {:<5.0f} {:>40}'.format(var.xGlobal[iMarker-1],iMarker-1, '|'))
- lamSol = var.u[(iMarker-1)*var.nVar : (iMarker)* var.nVar].copy()
- avTurb = self.transSolver.transitionBC(var, iMarker - 1)
- # Compute turbulent solution
- var.flowRegime[iMarker - 1] = 1
- self.newtonSolver(var, iMarker-1, displayState)
- # Average solutions
- var.u[(iMarker-1)*var.nVar : (iMarker)*var.nVar] = (1-avTurb) * lamSol + avTurb * var.u[(iMarker-1)*var.nVar : (iMarker)*var.nVar]
- # Recompute closures
- var.blPar[(iMarker-1)*nBlPar+1:(iMarker-1)*nBlPar+9]=Closures.turbulentClosure(var.u[(iMarker-1)*var.nVar+0],
- var.u[(iMarker-1)*var.nVar+1],
- var.u[(iMarker-1)*var.nVar+4],
- var.blPar[(iMarker-1)*nBlPar+0],
- var.extPar[(iMarker-1)*nExtPar+0],
- 'airfoil')"""
- pass
\ No newline at end of file