From 0e0665171131320b094b065012d86666e69f952b Mon Sep 17 00:00:00 2001
From: Dechamps Paul <paul.dechamps@student.uliege.be>
Date: Mon, 29 Nov 2021 15:37:53 +0100
Subject: [PATCH] Improvements and renaming.

---
 dart/tests/bliNonLift.py                      |  21 +-
 dart/tests/bliSeparated.py                    |  22 +-
 dart/tests/bliTransonic.py                    |  18 +-
 dart/viscous/Drivers/DDriver.py               | 276 ++++++++++
 ...oupConstructor.py => DGroupConstructor.py} |   4 +-
 .../{VPostProcess.py => DPostProcess.py}      |   0
 dart/viscous/Drivers/VDriver.py               | 291 ----------
 .../Master/{VAirfoil.py => DAirfoil.py}       |   2 +-
 .../Master/{VBoundary.py => DBoundary.py}     |   0
 .../Master/{VCoupler.py => DCoupler.py}       |   7 +-
 dart/viscous/Master/{VWake.py => DWake.py}    |   2 +-
 .../{VBIConditions.py => DBIConditions.py}    |   4 +-
 .../Physics/{VClosures.py => DClosures.py}    |   0
 ...{VDataStructures.py => DDataStructures.py} |   0
 ...{VDiscretization.py => DDiscretization.py} |   0
 .../{VValidation.py => DValidation.py}        |   0
 .../Physics/{VVariables.py => DVariables.py}  |   0
 dart/viscous/Solvers/DIntegration.py          | 288 ++++++++++
 dart/viscous/Solvers/DSolver.py               | 167 ++++++
 .../Solvers/{VSolver.py => DSysMatrix.py}     |  67 +--
 .../{VTransition.py => DTransition.py}        |   2 +-
 dart/viscous/Solvers/VTimeSolver.py           | 497 ------------------
 22 files changed, 799 insertions(+), 869 deletions(-)
 create mode 100644 dart/viscous/Drivers/DDriver.py
 rename dart/viscous/Drivers/{VGroupConstructor.py => DGroupConstructor.py} (99%)
 rename dart/viscous/Drivers/{VPostProcess.py => DPostProcess.py} (100%)
 delete mode 100644 dart/viscous/Drivers/VDriver.py
 rename dart/viscous/Master/{VAirfoil.py => DAirfoil.py} (99%)
 rename dart/viscous/Master/{VBoundary.py => DBoundary.py} (100%)
 rename dart/viscous/Master/{VCoupler.py => DCoupler.py} (96%)
 rename dart/viscous/Master/{VWake.py => DWake.py} (98%)
 rename dart/viscous/Physics/{VBIConditions.py => DBIConditions.py} (99%)
 rename dart/viscous/Physics/{VClosures.py => DClosures.py} (100%)
 rename dart/viscous/Physics/{VDataStructures.py => DDataStructures.py} (100%)
 rename dart/viscous/Physics/{VDiscretization.py => DDiscretization.py} (100%)
 rename dart/viscous/Physics/{VValidation.py => DValidation.py} (100%)
 rename dart/viscous/Physics/{VVariables.py => DVariables.py} (100%)
 create mode 100644 dart/viscous/Solvers/DIntegration.py
 create mode 100644 dart/viscous/Solvers/DSolver.py
 rename dart/viscous/Solvers/{VSolver.py => DSysMatrix.py} (81%)
 rename dart/viscous/Solvers/{VTransition.py => DTransition.py} (99%)
 delete mode 100644 dart/viscous/Solvers/VTimeSolver.py

diff --git a/dart/tests/bliNonLift.py b/dart/tests/bliNonLift.py
index 456ab36..8ac5aef 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 50181f5..47d4fd1 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 21a5ede..c923b03 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 0000000..c2bcf7d
--- /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 18f96c2..c03b482 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 1e9a017..0000000
--- 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 226091d..bf1431f 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 dcea5fb..da27a48 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 3a6cfe4..9ec5b93 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 7f9efbb..c1a145b 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 0000000..f1d70c8
--- /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 0000000..39865e3
--- /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 17f16c9..0e55194 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 461f890..e59b9a9 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 8643dcf..0000000
--- 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
-- 
GitLab