dart/viscous/README.md

+# DARTFlo Viscous-Inviscid Interaction.
+The folder dartflo/dart/viscous/ is dedicated to the viscous inviscid interaction implemented in DARTflo. This section extends DARTflo capabilities to two dimensional steady high Reynolds number transitional viscous flows.
+The solver makes use of the integral boundary layer equations (IBL) and the e^N method to correct the inviscid flow solution computed by DARTflo. The solver is divded into two parts. One that solves the steady state IBL and one that deals with the time-dependent IBL solved through a time marching procedure. The communication between the solvers is always done through steady-state solutions. The interaction is able to deal with midly separated flows for which the boundary layer assumptions still hold. 
+The two flow domains are coupled through a quasi-simulataneous interaction law that prevents Goldstein singularity. 
+
+![](/dox/title.png)
+
+## Main features
+* 2D viscous flow.
+* Free/forced laminar to turbulent transition.
+* Solver models
+  - Steady state equations solution. 
+  - 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 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.
+
+Crovato Adrien, [DARTFlo - Discrete Adjoint for Rapid Transonic Flows](http://hdl.handle.net/2268/264284), Technical note, University of Liège, 2021.
+
+Crovato A., Boman R., et al., [A Full Potential Static Aeroelastic Solver for Preliminary Aircraft Design](http://hdl.handle.net/2268/237955), 18th International Forum on Aeroelasticity and Structural Dynamics, IFASD 2019.
+
+Bilocq Amaury, [Implementation of a viscous-inviscid interaction scheme in a finite element full potential solver](http://hdl.handle.net/2268/252195), Master thesis, University of Liège, 2020.

dart/viscous/Solvers/DIntegration.py

+################################################################################
+#                                                                              #
+#                           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 matplotlib.colors import Normalize
+from dart.viscous.Physics.DClosures import Closures
+
+import math
+import numpy as np
+from matplotlib import pyplot as plt
+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-8
+    self.itFrozenJac0 = 5
+    
+    self.__Minf = max(_Minf, 0.1)
+    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, cfgFlow, 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 = cfgFlow.u.copy()
+
+    # Retreive parameters.
+
+    nVar = cfgFlow.nVar                                             # Number of variables of the differential problem.
+    nExtPar = cfgFlow.nExtPar                                       # Number of boundary layer parameters used for vector indexing.
+    bNodes = cfgFlow.bNodes                                         # Boundary nodes of the computational domain.
+    iInvVel = cfgFlow.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 = cfgFlow.xx[iMarker] - cfgFlow.xx[iMarkerPrev]              # Cell size.
+
+    # Initial residual.
+
+    sysRes = self.sysMatrix.buildResiduals(cfgFlow, iMarker)
+    normSysRes0 = np.linalg.norm(sysRes)
+
+    # Initialize pseudo-time integration parameters.
+
+    CFL = self.__CFL0                                             # Initialized CFL number.
+    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.
+    convPlot = []
+    # Main loop.
+
+    nErrors = 0                                                   # Count 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(cfgFlow, iMarker,
+                                                sysRes, order = 1, eps = 1e-8)
+
+        # Linear solver.
+
+        slnIncr = np.linalg.solve((Jacobian + IMatrix), - sysRes)
+
+        # Increment solution.
+
+        cfgFlow.u[iMarker*nVar : (iMarker+1)*nVar] += slnIncr
+
+        # Residual computation.
+
+        sysRes = self.sysMatrix.buildResiduals(cfgFlow, iMarker)
+        normSysRes = np.linalg.norm(sysRes + slnIncr/timeStep)
+        convPlot.append(normSysRes/normSysRes0)
+
+        # Check convergence.
+        
+        if normSysRes <= self.__tol * normSysRes0:
+          # Convergence criterion satisfied.
+          return 0
+
+      except KeyboardInterrupt:
+        print(ccolors.ANSI_RED + 'Program terminated by user.', ccolors.ANSI_RESET)
+        quit()
+        
+      except Exception:
+
+        nErrors += 1
+
+        if nErrors >= 5:
+          self.__ResetSolution(cfgFlow, iMarker, sln0)
+          return -1
+
+        cfgFlow.u[iMarker*nVar : (iMarker+1)*nVar] = cfgFlow.u[(iMarker-1)*nVar : (iMarker)*nVar]		# Reset solution.
+
+        # Lower CFL and recompute time step
+
+        if math.isnan(CFL):
+          CFL = self.__CFL0
+        CFL = min(max(CFL / 2,0.3),1)
+        timeStep = self.__SetTimeStep(CFL, dx, iInvVel)
+        IMatrix = np.identity(nVar) / timeStep
+
+        sysRes = self.sysMatrix.buildResiduals(cfgFlow, iMarker)
+
+        itFrozenJac = 1
+        continue
+
+      # CFL adaptation.
+
+      CFL = self.__CFL0* (normSysRes0/normSysRes)**0.7
+
+      CFL = max(CFL, 1)
+      timeStep = self.__SetTimeStep(CFL, dx, iInvVel)
+      IMatrix = np.identity(nVar) / timeStep
+
+
+      innerIters += 1
+
+    else: # Maximum number of iterations reached.
+
+      if normSysRes >= 1e-3 * normSysRes0:
+        self.__ResetSolution(cfgFlow, iMarker, sln0)
+      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 __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)

dart/viscous/Solvers/DSolver.py

+################################################################################
+#                                                                              #
+#                           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 = -1*np.ones(DFlow.nN)        # Convergence control of each point.
+    convergedPoints[0] = 0                        # Boundary condition.
+    
+    self.transSolver.initTransition(DFlow)        # Initialize the flow regime on each cell.
+
+    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 and DFlow.u[iMarker * DFlow.nVar + 4] == 0:
+        DFlow.u[iMarker * DFlow.nVar + 4] = 0.03
+
+      if iMarker == bNodes[0] + 1:	# Upper side start.
+         print('- Upper side,', end = ' ')
+      
+      if iMarker == bNodes[1]:	# Lower side start.
+        if lockTrans == 0:
+          print('no transition.')
+        lockTrans = 0
+        print('- Lower side,', end = ' ')
+
+      if iMarker == bNodes[2]: # First wake point
+
+        self.wakeBC(DFlow)	# Wake boundary condition.
+        convergedPoints[iMarker] = 0
+        if lockTrans == 0:
+          print('no transition.')
+        lockTrans = 1				
+
+        continue	# Ignore remaining instructions for the first wake point.
+      
+      # Call Newton solver for the point.
+
+      convergedPoints[iMarker] = self.integrator.TimeMarching(DFlow, iMarker, displayState)
+
+      if convergedPoints[iMarker] != 0:
+        print(ccolors.ANSI_RED + 'Marker {:<4.0f} (x/c={:<1.4f}): PTI terminated with exit code {:<4.0f}.'
+        .format(iMarker, DFlow.xGlobal[iMarker], convergedPoints[iMarker]), ccolors.ANSI_RESET)
+
+      # Check for transition.
+
+      if not lockTrans:
+
+        self.__CheckWaves(DFlow, iMarker)
+        
+        # Free transition.
+        if DFlow.u[iMarker * DFlow.nVar + 2] >= DFlow.Ncrit:
+          self.__AverageTransition(DFlow, iMarker, displayState)
+          if iMarker < DFlow.bNodes[1]:
+            print('free transition x/c = {:<1.5f} {:>4.0f}.'.format(DFlow.xtr[0],iMarker))
+          elif iMarker < DFlow.bNodes[2]:
+            print('free transition x/c = {:<1.5f} {:>4.0f}.'.format(DFlow.xtr[1],iMarker))
+          else:
+            print('error in wake')
+          lockTrans = 1
+        
+        # Forced transition.
+        elif DFlow.xtrF[0] is not None or DFlow.xtrF[1] is not None:
+          if iMarker != 1 and DFlow.flowRegime[iMarker] == 1 and DFlow.flowRegime[iMarker - 1] == 0:
+            print('forced transition near x/c = {:<2.3f} {:>4.0f}'.format(DFlow.xGlobal[iMarker],iMarker))
+            self.__AverageTransition(DFlow, iMarker, displayState, forced = 1)
+            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, forced = 0) -> 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.
+    if not forced:
+      avTurb = self.transSolver.solveTransition(var, iMarker)
+    else:
+      avTurb = self.transSolver.forcedTransition(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]
+
+    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

dart/viscous/Solvers/DSysMatrix.py

+################################################################################
+#                                                                              #
+#                            DARTFlo viscous module file                       #
+#               CSolver : Boundary layer laws and linear system sparse         #
+#                    matrices computation (Residuals, Jacobian)                #
+#                                                                              #
+# Author: Paul Dechamps                                                        #
+# Université de Liège                                                          #
+# Date: 2021.11.11                                                             #
+#                                                                              #
+################################################################################
+
+
+# ------------------------------------------------------------------------------
+#  Imports
+# ------------------------------------------------------------------------------
+from dart.viscous.Physics.DClosures import Closures
+
+import numpy as np
+
+class flowEquations:
+
+    def __init__(self, _setGradients, _fouScheme) -> None:
+        self.setGradients = _setGradients
+        self.fouScheme = _fouScheme
+        pass
+    
+    def _blLaws(self, dataCont, iMarker, x) -> np.ndarray:
+        """Evaluates the boundary layer laws at one point of the computational domain.
+
+            Args:
+            ----
+            - dataCont (classType): Data container. Contains the geometry, the solution, BL parameters...
+
+            - iMarker (integer): Id of the point.
+
+            - x (np.ndarray, optional): Prescribed solution, usually a perturbed one for Jacobian computation purposes. Defaults to None.
+
+            Returns:
+            -------
+            - linSysRes (np.ndarray): Residuals of the linear system at the considered point.
+            """
+
+        timeMatrix  = np.empty((dataCont.nVar, dataCont.nVar))
+        spaceMatrix = np.empty(dataCont.nVar)
+        
+        iMarkerPrev, name, xtrF = self.__evalPoint(dataCont, iMarker)
+        
+        dissipRatio = 0.9 if name == 'wake' else 1 # Wall/wake dissipation length ratio
+
+        dataCont.blPar[iMarker*dataCont.nBlPar+0] = max(x[3] * x[0] * dataCont.Re, 1)  # Reynolds number based on momentum thickness theta
+
+        dataCont.blPar[iMarker * dataCont.nBlPar + 9] = x[0] * x[1]  # deltaStar
+
+        if dataCont.flowRegime[iMarker]==0: 
+
+            # Laminar closure relations.
+            dataCont.blPar[iMarker*dataCont.nBlPar+1:iMarker*dataCont.nBlPar+7] = Closures.laminarClosure(x[0], x[1],
+                                                                                                         dataCont.blPar[iMarker*dataCont.nBlPar+0],
+                                                                                                         dataCont.extPar[iMarker*dataCont.nExtPar+0], name)
+
+        elif dataCont.flowRegime[iMarker]==1:
+
+            # Turbulent closures relations.
+            dataCont.blPar[iMarker*dataCont.nBlPar+1:iMarker*dataCont.nBlPar+9] = Closures.turbulentClosure(x[0], x[1],
+                                                                                                          x[4], dataCont.blPar[iMarker*dataCont.nBlPar+0],
+                                                                                                          dataCont.extPar[iMarker*dataCont.nExtPar+0], name)
+
+        # Rename variables for clarity.
+        Ta, Ha, _, uea, Cta = x
+        Mea = dataCont.extPar[iMarker*dataCont.nExtPar+0]
+        Cda, Cfa, deltaa, Hstara, Hstar2a, Hka, Cteqa, Usa = dataCont.blPar[iMarker * dataCont.nBlPar + 1 : iMarker * dataCont.nBlPar + 9]
+
+        # Amplification factor for e^n method.
+        if name=='airfoil' and xtrF is None and dataCont.flowRegime[iMarker]==0:
+
+            Axa = Closures.amplRoutine(dataCont.blPar[iMarker * dataCont.nBlPar + 6], dataCont.blPar[iMarker * dataCont.nBlPar + 0], Ta)
+        
+        else:
+
+            Axa=0
+            dN_dx=0
+
+        [dTheta_dx, dH_dx, dN_dx, due_dx, dCt_dx], dHstar_dx, dueExt_dx, ddeltaStarExt_dx, c, cExt = self.fouScheme(dataCont, x,
+                                                                                                                      iMarkerPrev,
+                                                                                                                      iMarker)
+
+        # +------------------------------------------- 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-Mea**2) * (Ta/uea) * due_dx - Cfa/2
+        spaceMatrix[1] = Ta * dHstar_dx + ( 2*Hstar2a + Hstara * (1-Ha)) * Ta/uea * due_dx - 2*Cda + Hstara * Cfa/2
+        if xtrF is not None or dataCont.flowRegime[iMarker] == 1:
+            spaceMatrix[2] = 0
+        else:
+            spaceMatrix[2] = dN_dx - Axa
+        spaceMatrix[3] = due_dx -c*(Ha * dTheta_dx + Ta * dH_dx) - dueExt_dx + cExt * ddeltaStarExt_dx
+        if dataCont.flowRegime[iMarker] == 1:
+
+            spaceMatrix[4] = (2*deltaa/Cta) * (dCt_dx) - 5.6*(Cteqa - Cta * dissipRatio)
+            - 2*deltaa*(4/(3*Ha * Ta) * (Cfa/2 - ((Hka-1)/(6.7*Hka * dissipRatio))**2) - 1/uea * due_dx)
+
+        elif dataCont.flowRegime[iMarker] == 0: spaceMatrix[4] = 0
+
+        # Temporal part
+        timeMatrix[0] = [Ha/uea, Ta/uea, 0, Ta * Ha/(uea**2), 0]
+        timeMatrix[1] = [(1+Ha * (1-Hstara))/uea, (1-Hstara) * Ta/uea, 0, (2-Hstara * Ha)*Ta/(uea**2), 0]
+        timeMatrix[2] = [0, 0, 1, 0, 0]
+        timeMatrix[3] = [-c*Ha, -c*Ta, 0, 1, 0]
+        if dataCont.flowRegime[iMarker]==1:
+            timeMatrix[4] = [0, 0, 0, 2*deltaa/(Usa * uea**2), 2*deltaa/(uea * Cta * Usa)]
+        else:
+            # Shear lag equation ignored in the laminar portion of the flow
+            timeMatrix[4] = [0, 0, 0, 0, 1]
+            
+        sysRes = np.linalg.solve(timeMatrix, spaceMatrix).flatten()
+        
+        """sysRes = np.empty(5)
+        if dataCont.flowRegime[iMarker] == 0: # Laminar case
+          sysRes[0] = (((-((c*Ha*Ta**2)/uea**2) - Ta/uea)*(-2*Cda + (Cfa*Hstara)/2. + Ta*dHstar_dx + 
+                      ((2*Hstar2a + (1 - Ha)*Hstara)*Ta*due_dx)/uea))/(-((c*Ha*Ta**2)/uea**3) - Ta/uea**2) + 
+                       (((2*c*Ta**2)/uea**2 - (c*Ha*Hstara*Ta**2)/uea**2 + Ta/uea - (Hstara*Ta)/uea)*
+                      (-0.5*Cfa + dTheta_dx + ((2 + Ha - Mea**2)*Ta*due_dx)/uea))/(-((c*Ha*Ta**2)/uea**3) - Ta/uea**2) + 
+                       (((2*Ta**2)/uea**3 - (Ha*Ta**2)/uea**3)*(cExt*ddeltaStarExt_dx - c*(Ta*dH_dx + Ha*dTheta_dx) + 
+                      due_dx - dueExt_dx))/(-((c*Ha*Ta**2)/uea**3) - Ta/uea**2))
+
+          sysRes[1] = ((((c*Ha**2*Ta)/uea**2 + Ha/uea)*(-2*Cda + (Cfa*Hstara)/2. + Ta*dHstar_dx + ((2*Hstar2a + (1 - Ha)*Hstara)*Ta*due_dx)/uea))/
+                      (-((c*Ha*Ta**2)/uea**3) - Ta/uea**2) + (((-2*c*Ha*Ta)/uea**2 + (c*Ha**2*Hstara*Ta)/uea**2 - 1/uea - Ha/uea + (Ha*Hstara)/uea)*
+                      (-0.5*Cfa + dTheta_dx + ((2 + Ha - Mea**2)*Ta*due_dx)/uea))/(-((c*Ha*Ta**2)/uea**3) - Ta/uea**2) + 
+                       ((-((Ha*Ta)/uea**3) + (Ha**2*Ta)/uea**3)*(cExt*ddeltaStarExt_dx - c*(Ta*dH_dx + Ha*dTheta_dx) + 
+                      due_dx - dueExt_dx))/(-((c*Ha*Ta**2)/uea**3) - Ta/uea**2))
+          
+          sysRes[2] = (-Axa + dN_dx)
+
+          sysRes[3] = (-((c*Ta*(-0.5*Cfa + dTheta_dx + ((2 + Ha - Mea**2)*Ta*due_dx)/uea))/((-((c*Ha*Ta**2)/uea**3) - Ta/uea**2)*uea)) - 
+                       (Ta*(cExt*ddeltaStarExt_dx - c*(Ta*dH_dx + Ha*dTheta_dx) + due_dx - 
+                      dueExt_dx))/((-((c*Ha*Ta**2)/uea**3) - Ta/uea**2)*uea**2))
+          
+          sysRes[4] = 0
+
+        elif dataCont.flowRegime[iMarker] == 1: # Turbulent case
+
+          sysRes[0] = ((((-2*c*deltaa*Ha*Ta**2)/(Cta*uea**3*Usa) - (2*deltaa*Ta)/(Cta*uea**2*Usa))*
+                      (-2*Cda + (Cfa*Hstara)/2. + Ta*dHstar_dx + ((2*Hstar2a + (1 - Ha)*Hstara)*Ta*due_dx)/uea))/
+                      ((-2*c*deltaa*Ha*Ta**2)/(Cta*uea**4*Usa) - (2*deltaa*Ta)/(Cta*uea**3*Usa)) + 
+                       (((4*c*deltaa*Ta**2)/(Cta*uea**3*Usa) - (2*c*deltaa*Ha*Hstara*Ta**2)/(Cta*uea**3*Usa) + (2*deltaa*Ta)/(Cta*uea**2*Usa) - 
+                      (2*deltaa*Hstara*Ta)/(Cta*uea**2*Usa))*(-0.5*Cfa + dTheta_dx + ((2 + Ha - Mea**2)*Ta*due_dx)/uea))/
+                      ((-2*c*deltaa*Ha*Ta**2)/(Cta*uea**4*Usa) - (2*deltaa*Ta)/(Cta*uea**3*Usa)) + 
+                       (((4*deltaa*Ta**2)/(Cta*uea**4*Usa) - (2*deltaa*Ha*Ta**2)/(Cta*uea**4*Usa))*
+                      (cExt*ddeltaStarExt_dx - c*(Ta*dH_dx + Ha*dTheta_dx) + due_dx - dueExt_dx))/
+                      ((-2*c*deltaa*Ha*Ta**2)/(Cta*uea**4*Usa) - (2*deltaa*Ta)/(Cta*uea**3*Usa)))
+           
+          sysRes[1] = ((((2*c*deltaa*Ha**2*Ta)/(Cta*uea**3*Usa) + (2*deltaa*Ha)/(Cta*uea**2*Usa))*
+                      (-2*Cda + (Cfa*Hstara)/2. + Ta*dHstar_dx + ((2*Hstar2a + (1 - Ha)*Hstara)*Ta*due_dx)/uea))/
+                      ((-2*c*deltaa*Ha*Ta**2)/(Cta*uea**4*Usa) - (2*deltaa*Ta)/(Cta*uea**3*Usa)) - 
+                       (2*deltaa*((2*c*Ha*Ta)/uea**2 - (c*Ha**2*Hstara*Ta)/uea**2 + 1/uea + Ha/uea - (Ha*Hstara)/uea)*
+                      (-0.5*Cfa + dTheta_dx + ((2 + Ha - Mea**2)*Ta*due_dx)/uea))/
+                      (Cta*uea*((-2*c*deltaa*Ha*Ta**2)/(Cta*uea**4*Usa) - (2*deltaa*Ta)/(Cta*uea**3*Usa))*Usa) + 
+                       (((-2*deltaa*Ha*Ta)/(Cta*uea**4*Usa) + (2*deltaa*Ha**2*Ta)/(Cta*uea**4*Usa))*
+                      (cExt*ddeltaStarExt_dx - c*(Ta*dH_dx + Ha*dTheta_dx) + due_dx - dueExt_dx))/
+                      ((-2*c*deltaa*Ha*Ta**2)/(Cta*uea**4*Usa) - (2*deltaa*Ta)/(Cta*uea**3*Usa)))
+
+          sysRes[2] = 0
+
+          sysRes[3] = ((-2*c*deltaa*Ta*(-0.5*Cfa + dTheta_dx + ((2 + Ha - Mea**2)*Ta*due_dx)/uea))/
+                      (Cta*uea**2*((-2*c*deltaa*Ha*Ta**2)/(Cta*uea**4*Usa) - (2*deltaa*Ta)/(Cta*uea**3*Usa))*Usa) - 
+                      (2*deltaa*Ta*(cExt*ddeltaStarExt_dx - c*(Ta*dH_dx + Ha*dTheta_dx) + due_dx - 
+                      dueExt_dx))/(Cta*uea**3*((-2*c*deltaa*Ha*Ta**2)/(Cta*uea**4*Usa) - (2*deltaa*Ta)/(Cta*uea**3*Usa))*Usa))
+
+
+          sysRes[4] = ((2*c*deltaa*Ta*(-0.5*Cfa + dTheta_dx + ((2 + Ha - Mea**2)*Ta*due_dx)/uea))/
+                      (uea**3*((-2*c*deltaa*Ha*Ta**2)/(Cta*uea**4*Usa) - (2*deltaa*Ta)/(Cta*uea**3*Usa))*Usa) + 
+                       ((-((c*Ha*Ta**2)/uea**3) - Ta/uea**2)*(-5.6*(Cteqa - Cta*dissipRatio) + (2*deltaa*dCt_dx)/Cta - 
+                      2*deltaa*((4*(Cfa/2. - (0.02227667631989307*(-1 + Hka)**2)/(dissipRatio**2*Hka**2)))/(3.*Ha*Ta) - due_dx/uea)))/
+                      ((-2*c*deltaa*Ha*Ta**2)/(Cta*uea**4*Usa) - (2*deltaa*Ta)/(Cta*uea**3*Usa)) + 
+                       (2*deltaa*Ta*(cExt*ddeltaStarExt_dx - c*(Ta*dH_dx + Ha*dTheta_dx) + due_dx - 
+                      dueExt_dx))/(uea**4*((-2*c*deltaa*Ha*Ta**2)/(Cta*uea**4*Usa) - (2*deltaa*Ta)/(Cta*uea**3*Usa))*Usa))"""
+
+
+        return sysRes
+
+
+    def __evalPoint(self, dataCont, iMarker) -> tuple:
+        """Evaluates where the point is, the adjactent points
+            and the state of transition on the corresponding side
+
+            Args:
+            ----
+            - iMarker (integer): Id of the cell to evaluate.
+
+            - bNodes (np.ndarray): Boundary nodes of the computational domain.
+
+            - xtrIn (np.ndarray): Current transiton location on [upper, lower] sides.
+
+            - xtrFIn (np.ndarray): Forced transition state on [upper, lower] sides.
+
+            - nMarker ([type]): Number of cells in the computational domain.
+
+            Returns:
+            -------
+            - tuple: Id of the adjacent points, name of the corresponding group ('airfoil', 'wake'),
+                transition state on the corresponding side.
+            """
+        iMarkerPrev = 0 if iMarker == dataCont.bNodes[1] else iMarker - 1
+
+        if iMarker < dataCont.bNodes[1]:
+            name='airfoil'
+            xtrF= dataCont.xtrF[0]
+
+        elif iMarker < dataCont.bNodes[2]:
+            name = 'airfoil'
+            xtrF = dataCont.xtrF[1]
+
+        else:
+            name = 'wake'
+            xtrF = None
+        return iMarkerPrev, name, xtrF
+
+
+class sysMatrix(flowEquations):
+    def __init__(self, _nMarker, _nMarkers, setGradients, fouScheme) -> None:
+        flowEquations.__init__(self, setGradients, fouScheme)
+        self._nMarker = _nMarker
+        self._nMarkers = _nMarkers
+        pass
+
+    def buildJacobian(self, dataCont, marker, sysRes, order = 1, eps = 1e-10) -> np.ndarray:
+        """Contruct the Jacobian used for solving the liner system.
+        
+            Args:
+            ----
+            - dataCont (classType): Data container (solution, BL parameters, invicid flow parameters, geometry...)
+
+            - marker (integer, optional): Cell id on the computational domain. Default : -1 for whole domain.
+
+            - order (unsigned int, optional): Order of the discretization used to obtain the Jacobian (1st or 2nd order implemented).
+                                        If order == 1, linSysRes has to be specifed
+
+            - linSysRes (np.ndarray, optional): Residuals of the linear system at the current iteration. Used (and mandatory) if order == 1.
+
+            - eps (double, optional): Perturbation on the variables. Default : 10^(-8).
+
+            Returns
+            -------
+
+            - JacMatrix (np.ndarray): Jacobian of the system. 
+
+            """
+        
+        # Build Jacobian for one point 
+        if order == 1 and sysRes is not None:
+            
+            JacMatrix = np.zeros((dataCont.nVar, dataCont.nVar))
+            xUp = dataCont.u[marker*dataCont.nVar : (marker+1)*dataCont.nVar].copy()
+
+            for iVar in range(dataCont.nVar):
+                varSave = xUp[iVar]
+                xUp[iVar] += eps
+
+                JacMatrix[:,iVar] = (self._blLaws(dataCont, marker, x=xUp) - sysRes) / eps
+
+                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()
+
+            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)
+
+                xUp[iVar] = varSave
+                xDwn[iVar] = varSave
+            return JacMatrix
+
+        else:
+            raise RuntimeError('Only first and second orders Jacobian can be computed.')
+
+
+    def buildResiduals(self, dataCont, marker) -> np.ndarray:
+        """Construct residuals of the linear system.
+        
+            Args:
+            ----
+            
+            - dataCont (classType): Data container (solution, BL parameters, invicid flow parameters, geometry...)
+
+            - marker (int, optional): Cell id on the computational domain. default: -1 for whole domain.
+            """
+
+        # Build Residual for one point
+        return self._blLaws(dataCont, marker, x=dataCont.u[marker*dataCont.nVar:(marker+1)*dataCont.nVar])
\ No newline at end of file

dart/viscous/Solvers/DTransition.py

+################################################################################
+#                                                                              #
+#                           DARTFlo viscous module file                        #
+#                 Transition solver : Transition related components.           #
+#                                                                              #
+# Author: Paul Dechamps                                                        #
+# Université de Liège                                                          #
+# Date: 2021.09.10                                                             #
+#                                                                              #
+################################################################################
+
+
+# ------------------------------------------------------------------------------
+#  Imports
+# ------------------------------------------------------------------------------
+from dart.viscous.Physics.DClosures import Closures
+
+import numpy as np
+
+class Transition:
+    def __init__(self, xtrF):
+        self.upperForced = xtrF[0]
+        self.lowerForced = xtrF[1]
+
+    def initTransition(self, var):
+
+        var.flowRegime=np.zeros(var.nN)
+        var.flowRegime[var.bNodes[2]:]=1
+        
+        for n in range(len(var.xtrF)):
+            if var.xtrF[n] is not None:
+                if var.xtrF[n]==1: # Forced laminar flow
+                    var.transMarkers[n]=var.bNodes[n+1]
+                    var.xtr[n]=1
+                    var.flowRegime[var.bNodes[n]:var.bNodes[n+1]]=0
+                elif var.xtrF[n]==0: # Forced turbulent flow
+                    var.transMarkers[n]=var.bNodes[n]
+                    var.xtr[n]=0
+                    var.flowRegime[var.bNodes[n]:var.bNodes[n+1]]=1
+                elif 0<var.xtrF[n]<1:
+                    idx0=(abs(var.xGlobal[n*var.bNodes[1]:var.bNodes[1+n]]-0)).argmin()
+                    var.transMarkers[n]=(abs(var.xGlobal[n*var.bNodes[1]+idx0:var.bNodes[1+n]]
+                    -var.xtrF[n])).argmin()+n*var.bNodes[1]+idx0
+                    var.flowRegime[var.transMarkers[n]:var.bNodes[n+1]]=1
+                    var.xtr[n]=var.xtrF[n]
+                else:
+                    raise ValueError('Non-physical forced transition location')    
+
+
+    def solveTransition(self,dataCont,transMarker):
+      
+      dataCont.flowRegime[dataCont.bNodes[2]:] = 1		# Wake 
+
+      # Upper side
+      if dataCont.bNodes[0]<transMarker<dataCont.bNodes[1]:
+          
+        # Set regime to turbulent on remaining upper side points
+        dataCont.transMarkers[0] = transMarker
+        dataCont.flowRegime[transMarker : dataCont.bNodes[1]] = 1
+
+        # Compute upper side transition location
+        dataCont.xtr[0] = (dataCont.Ncrit-dataCont.u[(transMarker-1)*dataCont.nVar+2])*((dataCont.xGlobal[transMarker]-dataCont.xGlobal[transMarker-1])/(dataCont.u[transMarker*dataCont.nVar+2]-dataCont.u[(transMarker-1)*dataCont.nVar+2]))+dataCont.xGlobal[transMarker-1]
+        avTurb = (dataCont.xGlobal[transMarker]-dataCont.xtr[0])/(dataCont.xGlobal[dataCont.transMarkers[0]]-dataCont.xGlobal[dataCont.transMarkers[0]-1]) # percentage of the subsection corresponding to a turbulent flow
+    
+      # Fully turbulent flow on lower side
+      elif transMarker == dataCont.bNodes[1]:
+        dataCont.transMarkers[1] = dataCont.bNodes[1]
+        dataCont.xtr[1] = 0
+        dataCont.flowRegime[dataCont.bNodes[1] : dataCont.bNodes[2]] = 1
+
+      # Lower side
+      elif dataCont.bNodes[1] <= transMarker < dataCont.bNodes[2]:
+
+        # Set regime to turbulent on remaining lower side points
+        dataCont.transMarkers[1] = transMarker
+        dataCont.flowRegime[transMarker : dataCont.bNodes[2]] = 1
+
+        # Compute lower side transition location
+        dataCont.xtr[1] = (dataCont.Ncrit - dataCont.u[(transMarker-1)*dataCont.nVar+2])*((dataCont.xGlobal[transMarker]-dataCont.xGlobal[transMarker-1])/(dataCont.u[transMarker*dataCont.nVar+2]-dataCont.u[(transMarker-1)*dataCont.nVar+2])) + dataCont.xGlobal[transMarker-1]
+        avTurb = (dataCont.xGlobal[transMarker] - dataCont.xtr[1])/(dataCont.xGlobal[dataCont.transMarkers[1]]-dataCont.xGlobal[dataCont.transMarkers[1]-1]) # percentage of the subsection corresponding to a turbulent flow
+    
+      # Boundary condition
+      avLam = 1- avTurb # percentage of the subsection corresponding to a laminar flow
+      dataCont.blPar[(transMarker-1)*dataCont.nBlPar + 7] = Closures.turbulentClosure(dataCont.u[(transMarker-1)*dataCont.nVar+0], dataCont.u[(transMarker-1)*dataCont.nVar+1],0.03, dataCont.blPar[(transMarker-1)*dataCont.nBlPar+0], dataCont.extPar[(transMarker-1)*dataCont.nExtPar+0],'airfoil')[6]
+      dataCont.u[(transMarker-1)*dataCont.nVar + 4] = 0.7*avLam*dataCont.blPar[(transMarker-1)*dataCont.nBlPar +7]
+
+      return avTurb
+
+    def forcedTransition(self, dataCont, marker):
+
+      if marker < dataCont.bNodes[1]:
+        avTurb = (dataCont.xGlobal[marker]-dataCont.xtrF[0])/(dataCont.xGlobal[marker]-dataCont.xGlobal[marker-1]) # Percentage of the subsection corresponding to a turbulent flow
+      elif marker < dataCont.bNodes[2]:
+        avTurb = (dataCont.xGlobal[marker]-dataCont.xtrF[1])/(dataCont.xGlobal[marker]-dataCont.xGlobal[marker-1]) # Percentage of the subsection corresponding to a turbulent flow
+      else:
+        print('Trying to impose transition in the wake')
+      avLam = 1- avTurb # Percentage of the subsection corresponding to a laminar flow
+
+      # Boundary condition
+      dataCont.blPar[(marker-1)*dataCont.nBlPar + 7] = Closures.turbulentClosure(dataCont.u[(marker-1)*dataCont.nVar+0], dataCont.u[(marker-1)*dataCont.nVar+1],0.03, dataCont.blPar[(marker-1)*dataCont.nBlPar+0], dataCont.extPar[(marker-1)*dataCont.nExtPar+0],'airfoil')[6]
+      dataCont.u[(marker-1)*dataCont.nVar + 4] = 0.7*avLam*dataCont.blPar[(marker-1)*dataCont.nBlPar +7]
+
+      return avTurb
+
+    def transitionBC(self, dataCont, marker):
+      if marker < dataCont.bNodes[1]:
+          xtr  = dataCont.xtr[0]
+      elif dataCont.bNodes[1] <= marker < dataCont.bNodes[2]:
+          xtr = dataCont.xtr[1]
+
+      avTurb = (dataCont.xGlobal[marker]-dataCont.xtr[0])/(dataCont.xGlobal[marker]-dataCont.xGlobal[marker-1]) # Percentage of the subsection corresponding to a turbulent flow
+      avLam = 1- avTurb # Percentage of the subsection corresponding to a laminar flow
+      
+      # Boundary condition
+      dataCont.blPar[(marker-1)*dataCont.nBlPar + 7] = Closures.turbulentClosure(dataCont.u[(marker-1)*dataCont.nVar+0], dataCont.u[(marker-1)*dataCont.nVar+1],0.03, dataCont.blPar[(marker-1)*dataCont.nBlPar+0], dataCont.extPar[(marker-1)*dataCont.nExtPar+0],'airfoil')[6]
+      dataCont.u[(marker-1)*dataCont.nVar + 4] = 0.7*avLam*dataCont.blPar[(marker-1)*dataCont.nBlPar +7]
+
+      return avTurb
\ No newline at end of file

dart/viscous/airfoil.py → dart/viscous/Solvers/Steady/StdAirfoil.py

@@ -19,7 +19,7 @@
 ## Airfoil around which the boundary layer is computed
 # Amaury Bilocq
 
-from dart.viscous.boundary import Boundary
+from dart.viscous.Solvers.Steady.StdBoundary import Boundary
 import numpy as np
 
 class Airfoil(Boundary):

dart/viscous/Solvers/Steady/StdBoundary.py

+#!/usr/bin/env python3
+# -*- coding: utf-8 -*-
+
+# Copyright 2020 University of Liège
+#
+# Licensed under the Apache License, Version 2.0 (the "License");
+# you may not use this file except in compliance with the License.
+# You may obtain a copy of the License at
+#
+#     http://www.apache.org/licenses/LICENSE-2.0
+#
+# Unless required by applicable law or agreed to in writing, software
+# distributed under the License is distributed on an "AS IS" BASIS,
+# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+# See the License for the specific language governing permissions and
+# limitations under the License.
+
+
+## Base class representing a physical boundary
+# Amaury Bilocq
+
+import numpy as np
+
+class Boundary:
+    def __init__(self, _boundary):
+        self.boundary = _boundary # boundary of interest
+        self.nN = len(self.boundary.nodes) # number of nodes
+        self.nE = len(self.boundary.tag.elems) # number of elements
+        self.v = np.zeros((self.nN, 3)) # velocity at edge of BL
+        self.u = np.zeros(self.nE) # blowing velocity
+        self.Me = np.zeros(self.nN) # Mach number at the edge of the boundary layer
+        self.rhoe = np.zeros(self.nN) # density at the edge of the boundary layer
+        self.connectListnodes = np.zeros(self.nN) # connectivity for nodes
+        self.connectListElems = np.zeros(self.nE) # connectivity for elements
+        self.deltaStar = np.zeros(self.nN) # displacement thickness
+        self.xx = np.zeros(self.nN) # coordinates in the reference of the physical surface. Used to store the values for the interaction law

dart/viscous/coupler.py → dart/viscous/Solvers/Steady/StdCoupler.py

@@ -21,8 +21,8 @@
 
 import numpy as np
 from fwk.coloring import ccolors
-from dart.viscous.airfoil import Airfoil
-from dart.viscous.wake import Wake
+from dart.viscous.Solvers.Steady.StdAirfoil import Airfoil
+from dart.viscous.Solvers.Steady.StdWake import Wake
 
 class Coupler:
     def __init__(self, _isolver, _vsolver, _boundaryAirfoil, _boundaryWake, _tol, _writer):

dart/viscous/solver.py → dart/viscous/Solvers/Steady/StdSolver.py

@@ -39,9 +39,10 @@
 # -- check transonic predictions
 # -- implement quasi-2D method for 3D flows
 
-from dart.viscous.boundary import Boundary
-from dart.viscous.wake import Wake
-import dart.viscous.newton as Newton
+#from dart.viscous.Steady.StdBoundary import Boundary
+#from dart.viscous.Steady.StdWake import Wake
+from matplotlib import pyplot as plt
+import dart.viscous.Solvers.Steady.Newton as Newton
 import numpy as np
 from fwk.coloring import ccolors
 
@@ -78,6 +79,7 @@ class Solver:
         wData['delta'] = self.blVectors[9]
         wData['xtrTop'] = self.xtr[0]
         wData['xtrBot'] = self.xtr[1]
+        self.wData = wData
         return wData
 
     def writeFile(self):
@@ -463,6 +465,7 @@ class Solver:
             self.blScalars = np.add(ublScalars, lblScalars)
             self.blVectors = np.hstack((ublVectors, lblVectors))
             blwVel = np.concatenate((ublwVel, lblwVel))
+            print('AFblwVel',len(blwVel))
             self.xtr = [xtrTop, xtrBot]
             # Boundary conditions for the wake
             group.TEnd = [ublVectors[1][-1], lblVectors[1][-1]]
@@ -471,11 +474,16 @@ 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('WKblwVel',len(blwVel))
+            print('WKdeltaStarOut', len(blVectors[0]))
+            plt.plot(blwVel)
+            plt.show()
             group.deltaStar = blVectors[0]
             # The drag is measured at the end of the wake
             self.Cd = blScalars[0]
@@ -485,5 +493,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()]

dart/viscous/Solvers/Steady/StdWake.py

+#!/usr/bin/env python3
+# -*- coding: utf-8 -*-
+
+# Copyright 2020 University of Liège
+#
+# Licensed under the Apache License, Version 2.0 (the "License");
+# you may not use this file except in compliance with the License.
+# You may obtain a copy of the License at
+#
+#     http://www.apache.org/licenses/LICENSE-2.0
+#
+# Unless required by applicable law or agreed to in writing, software
+# distributed under the License is distributed on an "AS IS" BASIS,
+# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+# See the License for the specific language governing permissions and
+# limitations under the License.
+
+
+## Wake behind airfoil (around which the boundary layer is computed)
+# Amaury Bilocq
+
+from dart.viscous.Solvers.Steady.StdBoundary import Boundary
+from dart.viscous.Solvers.Steady.StdAirfoil import Airfoil
+import numpy as np
+
+class Wake(Boundary):
+    def __init__(self, _boundary):
+        Boundary.__init__(self, _boundary)
+        self.name = 'wake' # type of boundary
+        self.T0 = 0 # initial condition for the momentum thickness
+        self.H0 = 0
+        self.n0 = 9 # wake is always turbulent
+        self.Ct0 = 0
+
+
+    def connectList(self):
+        ''' Sort the value read by the viscous solver/ Create list of connectivity
+        '''
+        N1 = np.zeros(self.nN, dtype=int) # Node number
+        connectListNodes = np.zeros(self.nN, dtype=int) # Index in boundary.nodes
+        connectListElems = np.zeros(self.nE, dtype=int) # Index in boundary.elems
+        data = np.zeros((self.boundary.nodes.size(), 10))
+        i = 0
+        for n in self.boundary.nodes:
+            data[i,0] = n.no
+            data[i,1] = n.pos[0]
+            data[i,2] = n.pos[1]
+            data[i,3] = n.pos[2]
+            data[i,4] = self.v[i,0]
+            data[i,5] = self.v[i,1]
+            data[i,6] = self.Me[i]
+            data[i,7] = self.rhoe[i]
+            data[i,8] = self.deltaStar[i]
+            data[i,9] = self.xx[i]
+            i += 1
+        # Table containing the element and its nodes
+        eData = np.zeros((self.nE,3), dtype=int)
+        for i in range(0, self.nE):
+            eData[i,0] = self.boundary.tag.elems[i].no
+            eData[i,1] = self.boundary.tag.elems[i].nodes[0].no
+            eData[i,2] = self.boundary.tag.elems[i].nodes[1].no
+        connectListNodes = data[:,1].argsort()
+        connectListElems[0] = np.where(eData[:,1] == min(data[:,1]))[0]
+        for i in range(1, len(eData[:,0])):
+            connectListElems[i] = np.where(eData[:,1] == eData[connectListElems[i-1],2])[0]
+        data[:,0] = data[connectListNodes,0]
+        data[:,1] = data[connectListNodes,1]
+        data[:,2] = data[connectListNodes,2]
+        data[:,3] = data[connectListNodes,3]
+        data[:,4] = data[connectListNodes,4]
+        data[:,5] = data[connectListNodes,5]
+        data[:,6] = data[connectListNodes,6]
+        data[:,7] = data[connectListNodes,7]
+        data[:,8] = data[connectListNodes,8]
+        data[:,9] = data[connectListNodes,9]
+        # Separated upper and lower part
+        data = np.delete(data,0,1)
+        return connectListNodes, connectListElems, data

dart/viscous/__init__.py

-# -*- coding: utf-8 -*-