dart/tests/bliHighLift.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.
+
+
+## Compute lifting (linear or nonlinear) viscous flow around a NACA 0012
+# Amaury Bilocq
+#
+# Test the viscous-inviscid interaction scheme
+# Reference to the master's thesis: http://hdl.handle.net/2268/252195
+# Reference test cases with Naca0012 (different from master's thesis):
+# 1) Incompressible: Re = 1e7, M_inf = 0, alpha = 5°, msTE = 0.01, msLE = 0.01
+#    -> nIt = 27, Cl = 0.55, Cd = 0.0058, xtrTop = 0.087, xtrBot = 0.741
+#    -> msLe = 0.001, xtrTop = 0.061
+# 2) Compressible: Re = 1e7, M_inf = 5, alpha = 5°, msTE = 0.01, msLE = 0.01
+#    -> nIt = 33, Cl = 0.65, Cd = 0.0063, xtrTop = 0.057, xtrBot = 0.740
+# 3) Separated: Re = 1e7, M_inf = 0, alpha = 12°, msTE = 0.01, msLE = 0.001
+#    -> nIt = 43, Cl = 1.30 , Cd = 0.0108, xtrTop = 0.010, xtrBot = 1
+#
+# CAUTION
+# This test is provided to ensure that the solver works properly.
+# Mesh refinement may have to be performed to obtain physical results.
+
+import dart.utils as floU
+import dart.default as floD
+
+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.DValidation as validation
+
+import tbox
+import tbox.utils as tboxU
+import fwk
+from fwk.testing import *
+from fwk.coloring import ccolors
+
+import math
+
+def main():
+    # timer
+    tms = fwk.Timers()
+    tms['total'].start()
+
+    # define flow variables
+    Re = 3e6
+    alpha = 10.*math.pi/180
+    M_inf = 0.15
+    user_xtr=[0, 0]
+    user_Ti=None
+
+    mapMesh = 1
+    nFactor = 3
+
+    # Time solver parameters
+    Time_Domain=True # True for unsteady solver, False for steady solver
+    CFL0 = 1
+
+    # ========================================================================================
+
+    U_inf = [math.cos(alpha), math.sin(alpha)] # norm should be 1
+    c_ref = 1
+    dim = 2
+    tol = 1e-4 # tolerance for the VII (usually one drag count)
+    case='case15Xfoil.dat' # File containing results from XFOIL of the considered case
+
+
+    # mesh the geometry
+    print(ccolors.ANSI_BLUE + 'PyMeshing...' + ccolors.ANSI_RESET)
+    tms['msh'].start()
+    pars = {'xLgt' : 5, 'yLgt' : 5, 'msF' : 1, 'msTe' : 0.01, 'msLe' : 0.001}
+    msh, gmshWriter = floD.mesh(dim, 'models/n0012.geo', pars, ['field', 'airfoil', 'downstream'])
+    tms['msh'].stop()
+
+    # set the problem and add medium, initial/boundary conditions
+    tms['pre'].start()
+    pbl, _, _, bnd, blw = floD.problem(msh, dim, alpha, 0., M_inf, c_ref, c_ref, 0.25, 0., 0., 'airfoil', te = 'te', vsc = True, dbc=True)
+    tms['pre'].stop()
+
+    # solve viscous problem
+    print(ccolors.ANSI_BLUE + 'PySolving...' + ccolors.ANSI_RESET)
+    tms['solver'].start()
+    isolver = floD.newton(pbl)
+
+  # Choose between steady and unsteady solver
+    if Time_Domain is True: # Run unsteady solver
+        vsolver=floVSMaster.Driver(Re, user_xtr, CFL0, M_inf, case, mapMesh, nFactor)
+        coupler = floVC.Coupler(isolver, vsolver, blw[0], blw[1], tol, gmshWriter, bnd)
+
+    elif Time_Domain is False: # Run steady solver 
+        vsolver=floVS_Std.Solver(Re)
+        coupler = floVC_Std.Coupler(isolver, vsolver, blw[0], blw[1], tol, gmshWriter)
+    coupler.run()
+    tms['solver'].stop()
+
+    # extract Cp
+    Cp = floU.extract(bnd.groups[0].tag.elems, isolver.Cp)
+    tboxU.write(Cp, 'Cp_airfoil.dat', '%1.5e', ', ', 'x, y, z, Cp', '')
+    # display results
+    print(ccolors.ANSI_BLUE + 'PyRes...' + ccolors.ANSI_RESET)
+    print('      Re        M    alpha       Cl       Cd      Cdp      Cdf       Cm')
+    print('{0:6.1f}e6 {1:8.2f} {2:8.1f} {3:8.4f} {4:8.4f} {5:8.4f} {6:8.4f} {7:8.4f}'.format(Re/1e6, M_inf, alpha*180/math.pi, isolver.Cl, vsolver.Cd, vsolver.Cdp, vsolver.Cdf, isolver.Cm))
+
+    # display timers
+    tms['total'].stop()
+    print(ccolors.ANSI_BLUE + 'PyTiming...' + ccolors.ANSI_RESET)
+    print('CPU statistics')
+    print(tms)
+
+    # check results
+    print(ccolors.ANSI_BLUE + 'PyTesting...' + ccolors.ANSI_RESET)
+    tests = CTests()
+    if Re==3e6 and M_inf==0.15 and alpha==10*math.pi/180 and user_xtr==[0., 0.]:
+        tests.add(CTest('Cl', isolver.Cl, 1.08, 5e-2)) # RANS 1.1 Xfoil 1.077
+        tests.add(CTest('Cd', vsolver.Cd, 0.0134, 1e-3, forceabs=True)) # RANS 0.0132 XFOIL 0.0128
+        tests.add(CTest('Cdp', vsolver.Cdp, 0.0058, 1e-3, forceabs=True))
+        tests.add(CTest('xtr_top', vsolver.xtr[0], 0, 1e-3, forceabs=True))
+        tests.add(CTest('xtr_bot', vsolver.xtr[1], 0, 1e-3, forceabs=True))
+    else:
+        raise Exception('Test not defined for this flow')
+
+    tests.run()
+
+    # visualize solution and plot results
+    floD.initViewer(pbl)
+    tboxU.plot(Cp[:,0], Cp[:,3], 'x', 'Cp', 'Cl = {0:.{3}f}, Cd = {1:.{3}f}, Cm = {2:.{3}f}'.format(isolver.Cl, vsolver.Cd, isolver.Cm, 4), True)
+
+    val=validation.Validation()
+    val.main(isolver.Cl,vsolver.wData)
+
+    # eof
+    print('')
+
+if __name__ == "__main__":
+    main()

dart/tests/bliLowRe.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.
+
+
+## Compute lifting (linear or nonlinear) viscous flow around a NACA 0012
+# Amaury Bilocq
+#
+# Test the viscous-inviscid interaction scheme
+# Reference to the master's thesis: http://hdl.handle.net/2268/252195
+# Reference test cases with Naca0012 (different from master's thesis):
+# 1) Incompressible: Re = 1e7, M_inf = 0, alpha = 5°, msTE = 0.01, msLE = 0.01
+#    -> nIt = 27, Cl = 0.55, Cd = 0.0058, xtrTop = 0.087, xtrBot = 0.741
+#    -> msLe = 0.001, xtrTop = 0.061
+# 2) Compressible: Re = 1e7, M_inf = 5, alpha = 5°, msTE = 0.01, msLE = 0.01
+#    -> nIt = 33, Cl = 0.65, Cd = 0.0063, xtrTop = 0.057, xtrBot = 0.740
+# 3) Separated: Re = 1e7, M_inf = 0, alpha = 12°, msTE = 0.01, msLE = 0.001
+#    -> nIt = 43, Cl = 1.30 , Cd = 0.0108, xtrTop = 0.010, xtrBot = 1
+#
+# CAUTION
+# This test is provided to ensure that the solver works properly.
+# Mesh refinement may have to be performed to obtain physical results.
+
+import dart.utils as floU
+import dart.default as floD
+
+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.DValidation as validation
+
+import tbox
+import tbox.utils as tboxU
+import fwk
+from fwk.testing import *
+from fwk.coloring import ccolors
+
+import math
+
+def main():
+    # timer
+    tms = fwk.Timers()
+    tms['total'].start()
+
+    # define flow variables
+    Re = 5e5
+    alpha = 2.*math.pi/180
+    M_inf = 0.2
+
+    #user_xtr=[0.41,0.41]
+    #user_xtr=[0,None]
+    user_xtr=[None, None]
+    user_Ti=None
+
+    mapMesh = 1
+    nFactor = 2
+
+    # Time solver parameters
+    Time_Domain = True # True for unsteady solver, False for steady solver
+    CFL0 = 1
+
+    # ========================================================================================
+
+    U_inf = [math.cos(alpha), math.sin(alpha)] # norm should be 1
+    c_ref = 1
+    dim = 2
+    tol = 1e-4 # tolerance for the VII (usually one drag count)
+    case='Case1FreeTrans.dat' # File containing results from XFOIL of the considered case
+
+    # mesh the geometry
+    print(ccolors.ANSI_BLUE + 'PyMeshing...' + ccolors.ANSI_RESET)
+    tms['msh'].start()
+    pars = {'xLgt' : 5, 'yLgt' : 5, 'msF' : 1, 'msTe' : 0.01, 'msLe' : 0.0001}
+    #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()
+
+    # set the problem and add medium, initial/boundary conditions
+    tms['pre'].start()
+    pbl, _, _, bnd, blw = floD.problem(msh, dim, alpha, 0., M_inf, c_ref, c_ref, 0.25, 0., 0., 'airfoil', te = 'te', vsc = True, dbc=True)
+    tms['pre'].stop()
+
+    # solve viscous problem
+
+    print(ccolors.ANSI_BLUE + 'PySolving...' + ccolors.ANSI_RESET)
+    tms['solver'].start()
+    isolver = floD.newton(pbl)
+
+    # Choose between steady and unsteady solver
+    if Time_Domain is True: # Run unsteady solver
+        vsolver=floVSMaster.Driver(Re, user_xtr, CFL0, M_inf, case, mapMesh, nFactor)
+        coupler = floVC.Coupler(isolver, vsolver, blw[0], blw[1], tol, gmshWriter, bnd)
+
+    elif Time_Domain is False: # Run steady solver 
+        vsolver=floVS_Std.Solver(Re)
+        coupler = floVC_Std.Coupler(isolver, vsolver, blw[0], blw[1], tol, gmshWriter)
+    coupler.run()
+    tms['solver'].stop()
+
+    # extract Cp
+    Cp = floU.extract(bnd.groups[0].tag.elems, isolver.Cp)
+    tboxU.write(Cp, 'Cp_airfoil.dat', '%1.5e', ', ', 'x, y, z, Cp', '')
+    # display results
+    print(ccolors.ANSI_BLUE + 'PyRes...' + ccolors.ANSI_RESET)
+    print('      Re        M    alpha       Cl       Cd      Cdp      Cdf       Cm')
+    print('{0:6.1f}e6 {1:8.2f} {2:8.1f} {3:8.4f} {4:8.4f} {5:8.4f} {6:8.4f} {7:8.4f}'.format(Re/1e6, M_inf, alpha*180/math.pi, isolver.Cl, vsolver.Cd, vsolver.Cdp, vsolver.Cdf, isolver.Cm))
+
+    # display timers
+    tms['total'].stop()
+    print(ccolors.ANSI_BLUE + 'PyTiming...' + ccolors.ANSI_RESET)
+    print('CPU statistics')
+    print(tms)
+    
+    # check results
+    print(ccolors.ANSI_BLUE + 'PyTesting...' + ccolors.ANSI_RESET)
+    tests = CTests()
+    if Re == 5e5 and M_inf == 0.2 and alpha == 2*math.pi/180:
+        tests.add(CTest('Cl', isolver.Cl, 0.2206, 5e-2)) # XFOIL 0.2135
+        tests.add(CTest('Cd', vsolver.Cd, 0.007233, 1e-3, forceabs=True)) # XFOIL 0.00706
+        tests.add(CTest('Cdp', vsolver.Cdp, 0.0025, 1e-3, forceabs=True)) # 0.00238
+        tests.add(CTest('xtr_top', vsolver.xtr[0], 0.58, 0.03, forceabs=True)) # XFOIL 0.5642
+        tests.add(CTest('xtr_bot', vsolver.xtr[1], 0.90, 0.03, forceabs=True)) # XFOIL 0.9290
+    else:
+        print(ccolors.ANSI_RED + 'Test not defined for this flow', ccolors.ANSI_RESET)
+
+    tests.run()
+
+    # visualize solution and plot results
+    floD.initViewer(pbl)
+    tboxU.plot(Cp[:,0], Cp[:,3], 'x', 'Cp', 'Cl = {0:.{3}f}, Cd = {1:.{3}f}, Cm = {2:.{3}f}'.format(isolver.Cl, vsolver.Cd, isolver.Cm, 4), True)
+
+    val=validation.Validation(case)
+    val.main(isolver.Cl,vsolver.wData)
+
+    # eof
+    print('')
+
+if __name__ == "__main__":
+    main()

dart/tests/bliNonLift.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.
+
+
+## Compute lifting (linear or nonlinear) viscous flow around a NACA 0012
+# Amaury Bilocq
+#
+# Test the viscous-inviscid interaction scheme
+# Reference to the master's thesis: http://hdl.handle.net/2268/252195
+# Reference test cases with Naca0012 (different from master's thesis):
+# 1) Incompressible: Re = 1e7, M_inf = 0, alpha = 5°, msTE = 0.01, msLE = 0.01
+#    -> nIt = 27, Cl = 0.55, Cd = 0.0058, xtrTop = 0.087, xtrBot = 0.741
+#    -> msLe = 0.001, xtrTop = 0.061
+# 2) Compressible: Re = 1e7, M_inf = 5, alpha = 5°, msTE = 0.01, msLE = 0.01
+#    -> nIt = 33, Cl = 0.65, Cd = 0.0063, xtrTop = 0.057, xtrBot = 0.740
+# 3) Separated: Re = 1e7, M_inf = 0, alpha = 12°, msTE = 0.01, msLE = 0.001
+#    -> nIt = 43, Cl = 1.30 , Cd = 0.0108, xtrTop = 0.010, xtrBot = 1
+#
+# CAUTION
+# This test is provided to ensure that the solver works properly.
+# Mesh refinement may have to be performed to obtain physical results.
+
+import dart.utils as floU
+import dart.default as floD
+
+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.DValidation as validation
+
+import tbox
+import tbox.utils as tboxU
+import fwk
+from fwk.testing import *
+from fwk.coloring import ccolors
+
+import math
+
+def main():
+    # timer
+    tms = fwk.Timers()
+    tms['total'].start()
+
+    # define flow variables
+    Re = 1e7
+    alpha = 0.*math.pi/180
+    M_inf = 0.
+
+    #user_xtr=[0.41,0.41]
+    #user_xtr=[0,None]
+    user_xtr=[None,None]
+    user_Ti=None
+
+    mapMesh = 0
+    nFactor = 2
+
+    # Time solver parameters
+    Time_Domain = True # True for unsteady solver, False for steady solver
+    CFL0 = 5
+
+    # ========================================================================================
+
+    U_inf = [math.cos(alpha), math.sin(alpha)] # norm should be 1
+    c_ref = 1
+    dim = 2
+    tol = 1e-4 # tolerance for the VII (usually one drag count)
+    case='Case1FreeTrans.dat' # File containing results from XFOIL of the considered case
+
+    # mesh the geometry
+    print(ccolors.ANSI_BLUE + 'PyMeshing...' + ccolors.ANSI_RESET)
+    tms['msh'].start()
+    pars = {'xLgt' : 5, 'yLgt' : 5, 'msF' : 1, 'msTe' : 0.01, 'msLe' : 0.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()
+
+    # set the problem and add medium, initial/boundary conditions
+    tms['pre'].start()
+    pbl, _, _, bnd, blw = floD.problem(msh, dim, alpha, 0., M_inf, c_ref, c_ref, 0.25, 0., 0., 'airfoil', te = 'te', vsc = True, dbc=True)
+    tms['pre'].stop()
+
+    # solve viscous problem
+
+    print(ccolors.ANSI_BLUE + 'PySolving...' + ccolors.ANSI_RESET)
+    tms['solver'].start()
+    isolver = floD.newton(pbl)
+
+    # Choose between steady and unsteady solver
+    if Time_Domain is True: # Run unsteady solver
+        vsolver=floVSMaster.Driver(Re, user_xtr, CFL0, M_inf, case, mapMesh, nFactor)
+        coupler = floVC.Coupler(isolver, vsolver, blw[0], blw[1], tol, gmshWriter, bnd)
+
+    elif Time_Domain is False: # Run steady solver 
+        vsolver=floVS_Std.Solver(Re)
+        coupler = floVC_Std.Coupler(isolver, vsolver, blw[0], blw[1], tol, gmshWriter)
+    coupler.run()
+    tms['solver'].stop()
+
+    # extract Cp
+    Cp = floU.extract(bnd.groups[0].tag.elems, isolver.Cp)
+    tboxU.write(Cp, 'Cp_airfoil.dat', '%1.5e', ', ', 'x, y, z, Cp', '')
+    # display results
+    print(ccolors.ANSI_BLUE + 'PyRes...' + ccolors.ANSI_RESET)
+    print('      Re        M    alpha       Cl       Cd      Cdp      Cdf       Cm')
+    print('{0:6.1f}e6 {1:8.2f} {2:8.1f} {3:8.4f} {4:8.4f} {5:8.4f} {6:8.4f} {7:8.4f}'.format(Re/1e6, M_inf, alpha*180/math.pi, isolver.Cl, vsolver.Cd, vsolver.Cdp, vsolver.Cdf, isolver.Cm))
+
+    # display timers
+    tms['total'].stop()
+    print(ccolors.ANSI_BLUE + 'PyTiming...' + ccolors.ANSI_RESET)
+    print('CPU statistics')
+    print(tms)
+    
+    # check results
+    print(ccolors.ANSI_BLUE + 'PyTesting...' + ccolors.ANSI_RESET)
+    tests = CTests()
+    if Re == 1e7 and M_inf == 0.2 and alpha == 5*math.pi/180:
+        tests.add(CTest('Cl', isolver.Cl, 0.56, 5e-2)) # XFOIL 0.58
+        tests.add(CTest('Cd', vsolver.Cd, 0.0063, 1e-3, forceabs=True)) # XFOIL 0.00617
+        tests.add(CTest('Cdp', vsolver.Cdp, 0.0019, 1e-3, forceabs=True)) # XFOIL 0.00176
+        tests.add(CTest('xtr_top', vsolver.xtr[0], 0.059, 0.03, forceabs=True)) # XFOIL 0.0510
+        tests.add(CTest('xtr_bot', vsolver.xtr[1], 0.738, 0.03, forceabs=True)) # XFOIL 0.7473
+    else:
+        raise Exception('Test not defined for this flow')
+
+    tests.run()
+
+    # visualize solution and plot results
+    floD.initViewer(pbl)
+    tboxU.plot(Cp[:,0], Cp[:,3], 'x', 'Cp', 'Cl = {0:.{3}f}, Cd = {1:.{3}f}, Cm = {2:.{3}f}'.format(isolver.Cl, vsolver.Cd, isolver.Cm, 4), True)
+
+    val=validation.Validation(case)
+    val.main(isolver.Cl,vsolver.wData)
+
+    # eof
+    print('')
+
+if __name__ == "__main__":
+    main()

dart/tests/bliPolar.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.
+
+
+## Compute lifting (linear or nonlinear) viscous flow around a NACA 0012
+# Amaury Bilocq
+#
+# Test the viscous-inviscid interaction scheme
+# Reference to the master's thesis: http://hdl.handle.net/2268/252195
+# Reference test cases with Naca0012 (different from master's thesis):
+# 1) Incompressible: Re = 1e7, M_inf = 0, alpha = 5°, msTE = 0.01, msLE = 0.01
+#    -> nIt = 27, Cl = 0.55, Cd = 0.0058, xtrTop = 0.087, xtrBot = 0.741
+#    -> msLe = 0.001, xtrTop = 0.061
+# 2) Compressible: Re = 1e7, M_inf = 5, alpha = 5°, msTE = 0.01, msLE = 0.01
+#    -> nIt = 33, Cl = 0.65, Cd = 0.0063, xtrTop = 0.057, xtrBot = 0.740
+# 3) Separated: Re = 1e7, M_inf = 0, alpha = 12°, msTE = 0.01, msLE = 0.001
+#    -> nIt = 43, Cl = 1.30 , Cd = 0.0108, xtrTop = 0.010, xtrBot = 1
+#
+# CAUTION
+# This test is provided to ensure that the solver works properly.
+# Mesh refinement may have to be performed to obtain physical results.
+
+from numpy.core.function_base import linspace
+import dart.utils as floU
+import dart.default as floD
+
+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.DValidation as validation
+
+import tbox
+import tbox.utils as tboxU
+import fwk
+from fwk.testing import *
+from fwk.coloring import ccolors
+import numpy as np
+from matplotlib import pyplot as plt
+
+import math
+
+def main():
+    # timer
+    tms = fwk.Timers()
+    tms['total'].start()
+
+    # define flow variables
+    Re = 6e6
+    #alpha = 0.*math.pi/180
+    M_inf = 0.15
+
+    #user_xtr=[0.41,0.41]
+    #user_xtr=[0,None]
+    user_xtr=[0, 0]
+    user_Ti=None
+
+    mapMesh = 0
+    nFactor = 2
+
+    # Time solver parameters
+    Time_Domain = True # True for unsteady solver, False for steady solver
+    CFL0 = 1
+
+    # ========================================================================================
+
+    #U_inf = [math.cos(alpha), math.sin(alpha)] # norm should be 1
+    c_ref = 1
+    dim = 2
+    tol = 1e-4 # tolerance for the VII (usually one drag count)
+    case='Case1FreeTrans.dat' # File containing results from XFOIL of the considered case
+
+    # mesh the geometry
+    print(ccolors.ANSI_BLUE + 'PyMeshing...' + ccolors.ANSI_RESET)
+    tms['msh'].start()
+    pars = {'xLgt' : 5, 'yLgt' : 5, 'msF' : 1, 'msTe' : 0.01, 'msLe' : 0.0001}
+    #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()
+
+    alphaVect = [
+   15.2600*math.pi/180,
+   16.3000*math.pi/180,
+   16.5*math.pi/180,
+   17.1300*math.pi/180,
+   17.2*math.pi/180,
+   17.5*math.pi/180,
+   17.7*math.pi/180,
+   18.0200*math.pi/180]
+
+    alphaSuc = []
+    clSuc = []
+    cdSuc = []
+    alphaFailed = []
+    clFailed = []
+    cdFailed = []
+
+    for alpha in alphaVect:
+      # set the problem and add medium, initial/boundary conditions
+      tms['pre'].start()
+      pbl, _, _, bnd, blw = floD.problem(msh, dim, alpha, 0., M_inf, c_ref, c_ref, 0.25, 0., 0., 'airfoil', te = 'te', vsc = True, dbc=True)
+      tms['pre'].stop()
+
+      # solve viscous problem
+
+      print(ccolors.ANSI_BLUE + 'PySolving...' + ccolors.ANSI_RESET)
+      tms['solver'].start()
+      isolver = floD.newton(pbl)
+
+      # Choose between steady and unsteady solver
+      if Time_Domain is True: # Run unsteady solver
+          vsolver=floVSMaster.Driver(Re, user_xtr, CFL0, M_inf, case, mapMesh, nFactor)
+          coupler = floVC.Coupler(isolver, vsolver, blw[0], blw[1], tol, gmshWriter, bnd)
+
+      elif Time_Domain is False: # Run steady solver 
+          vsolver=floVS_Std.Solver(Re)
+          coupler = floVC_Std.Coupler(isolver, vsolver, blw[0], blw[1], tol, gmshWriter)
+      try:
+        coupler.run()
+        tms['solver'].stop()
+
+        # extract Cp
+        Cp = floU.extract(bnd.groups[0].tag.elems, isolver.Cp)
+        tboxU.write(Cp, 'Cp_airfoil.dat', '%1.5e', ', ', 'x, y, z, Cp', '')
+        # display results
+        print(ccolors.ANSI_BLUE + 'PyRes...' + ccolors.ANSI_RESET)
+        print('      Re        M    alpha       Cl       Cd      Cdp      Cdf       Cm')
+        print('{0:6.1f}e6 {1:8.2f} {2:8.1f} {3:8.4f} {4:8.4f} {5:8.4f} {6:8.4f} {7:8.4f}'.format(Re/1e6, M_inf, alpha*180/math.pi, isolver.Cl, vsolver.Cd, vsolver.Cdp, vsolver.Cdf, isolver.Cm))
+
+        # display timers
+        tms['total'].stop()
+        print(ccolors.ANSI_BLUE + 'PyTiming...' + ccolors.ANSI_RESET)
+        print('CPU statistics')
+        print(tms)
+
+        alphaSuc.append(alpha*180/math.pi)
+        clSuc.append(isolver.Cl)
+        cdSuc.append(vsolver.Cd)
+      except KeyboardInterrupt:
+        quit()
+      except:
+        alphaFailed.append(alpha*180/math.pi)
+        clFailed.append(isolver.Cl)
+        cdFailed.append(vsolver.Cd)
+
+
+    print(alphaSuc, clSuc, cdSuc)
+    print(alphaFailed, clFailed, cdFailed)
+
+if __name__ == "__main__":
+    main()

dart/tests/bliSeparated.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.
+
+
+## Compute lifting (linear or nonlinear) viscous flow around a NACA 0012
+# Amaury Bilocq
+#
+# Test the viscous-inviscid interaction scheme
+# Reference to the master's thesis: http://hdl.handle.net/2268/252195
+# Reference test cases with Naca0012 (different from master's thesis):
+# 1) Incompressible: Re = 1e7, M_inf = 0, alpha = 5°, msTE = 0.01, msLE = 0.01
+#    -> nIt = 27, Cl = 0.55, Cd = 0.0058, xtrTop = 0.087, xtrBot = 0.741
+#    -> msLe = 0.001, xtrTop = 0.061
+# 2) Compressible: Re = 1e7, M_inf = 5, alpha = 5°, msTE = 0.01, msLE = 0.01
+#    -> nIt = 33, Cl = 0.65, Cd = 0.0063, xtrTop = 0.057, xtrBot = 0.740
+# 3) Separated: Re = 1e7, M_inf = 0, alpha = 12°, msTE = 0.01, msLE = 0.001
+#    -> nIt = 43, Cl = 1.30 , Cd = 0.0108, xtrTop = 0.010, xtrBot = 1
+#
+# CAUTION
+# This test is provided to ensure that the solver works properly.
+# Mesh refinement may have to be performed to obtain physical results.
+
+import dart.utils as floU
+import dart.default as floD
+
+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.DValidation as validation
+
+import tbox
+import tbox.utils as tboxU
+import fwk
+from fwk.testing import *
+from fwk.coloring import ccolors
+
+import math
+
+def main():
+    # timer
+    tms = fwk.Timers()
+    tms['total'].start()
+
+    # define flow variables
+    Re = 1e7
+    alpha = 15*math.pi/180
+    M_inf = 0.2
+    user_xtr=[None, None]
+    user_Ti=None
+
+    mapMesh = 1
+    nFactor = 2
+
+    # Time solver parameters
+    Time_Domain=True # True for unsteady solver, False for steady solver
+    CFL0 = 1
+
+    # ========================================================================================
+
+    U_inf = [math.cos(alpha), math.sin(alpha)] # norm should be 1
+    c_ref = 1
+    dim = 2
+    tol = 1e-4 # tolerance for the VII (usually one drag count)
+    case='case15Xfoil.dat' # File containing results from XFOIL of the considered case
+
+
+    # mesh the geometry
+    print(ccolors.ANSI_BLUE + 'PyMeshing...' + ccolors.ANSI_RESET)
+    tms['msh'].start()
+    pars = {'xLgt' : 5, 'yLgt' : 5, 'msF' : 1, 'msTe' : 0.01, 'msLe' : 0.0001}
+    msh, gmshWriter = floD.mesh(dim, 'models/n0012.geo', pars, ['field', 'airfoil', 'downstream'])
+    tms['msh'].stop()
+
+    # set the problem and add medium, initial/boundary conditions
+    tms['pre'].start()
+    pbl, _, _, bnd, blw = floD.problem(msh, dim, alpha, 0., M_inf, c_ref, c_ref, 0.25, 0., 0., 'airfoil', te = 'te', vsc = True, dbc=True)
+    tms['pre'].stop()
+
+    # solve viscous problem
+    print(ccolors.ANSI_BLUE + 'PySolving...' + ccolors.ANSI_RESET)
+    tms['solver'].start()
+    isolver = floD.newton(pbl)
+
+  # Choose between steady and unsteady solver
+    if Time_Domain is True: # Run unsteady solver
+        vsolver=floVSMaster.Driver(Re, user_xtr, CFL0, M_inf, case, mapMesh, nFactor)
+        coupler = floVC.Coupler(isolver, vsolver, blw[0], blw[1], tol, gmshWriter, bnd)
+
+    elif Time_Domain is False: # Run steady solver 
+        vsolver=floVS_Std.Solver(Re)
+        coupler = floVC_Std.Coupler(isolver, vsolver, blw[0], blw[1], tol, gmshWriter)
+    coupler.run()
+    tms['solver'].stop()
+
+    # extract Cp
+    Cp = floU.extract(bnd.groups[0].tag.elems, isolver.Cp)
+    tboxU.write(Cp, 'Cp_airfoil.dat', '%1.5e', ', ', 'x, y, z, Cp', '')
+    # display results
+    print(ccolors.ANSI_BLUE + 'PyRes...' + ccolors.ANSI_RESET)
+    print('      Re        M    alpha       Cl       Cd      Cdp      Cdf       Cm')
+    print('{0:6.1f}e6 {1:8.2f} {2:8.1f} {3:8.4f} {4:8.4f} {5:8.4f} {6:8.4f} {7:8.4f}'.format(Re/1e6, M_inf, alpha*180/math.pi, isolver.Cl, vsolver.Cd, vsolver.Cdp, vsolver.Cdf, isolver.Cm))
+
+    # display timers
+    tms['total'].stop()
+    print(ccolors.ANSI_BLUE + 'PyTiming...' + ccolors.ANSI_RESET)
+    print('CPU statistics')
+    print(tms)
+
+    # check results
+    print(ccolors.ANSI_BLUE + 'PyTesting...' + ccolors.ANSI_RESET)
+    tests = CTests()
+    if Re == 1e7 and M_inf == 0.2 and alpha == 15*math.pi/180:
+      tests.add(CTest('Cl', isolver.Cl, 1.615, 5e-2)) # Xfoil 1.6654
+      tests.add(CTest('Cd', vsolver.Cd, 0.0143, 1e-3, forceabs=True)) # XFOIL 0.01460
+      tests.add(CTest('Cdp', vsolver.Cdp, 0.0103, 1e-3, forceabs=True)) # XFOIL 0.01255
+      tests.add(CTest('xtr_top', vsolver.xtr[0], 0.00660, 0.003, forceabs=True)) # XFOIL 0.0067
+      tests.add(CTest('xtr_bot', vsolver.xtr[1], 1, 0.03, forceabs=True)) # XFOIL 1.0000
+    elif Re == 5e6 and M_inf == 0.2 and alpha == 15*math.pi/180:
+      tests.add(CTest('Cl', isolver.Cl, 1.58, 5e-2)) # Xfoil 1.6109
+      tests.add(CTest('Cd', vsolver.Cd, 0.0201, 1e-3, forceabs=True)) # XFOIL 0.01944
+      tests.add(CTest('Cdp', vsolver.Cdp, 0.0130, 1e-3, forceabs=True)) # XFOIL 0.01544
+      tests.add(CTest('xtr_top', vsolver.xtr[0], 0.008, 0.003, forceabs=True)) # XFOIL 0.0081
+      tests.add(CTest('xtr_bot', vsolver.xtr[1], 1, 0.03, forceabs=True)) # XFOIL 1.0000
+    else:
+        raise Exception('Test not defined for this flow')
+
+    tests.run()
+
+    # visualize solution and plot results
+    floD.initViewer(pbl)
+    tboxU.plot(Cp[:,0], Cp[:,3], 'x', 'Cp', 'Cl = {0:.{3}f}, Cd = {1:.{3}f}, Cm = {2:.{3}f}'.format(isolver.Cl, vsolver.Cd, isolver.Cm, 4), True)
+
+    val=validation.Validation()
+    val.main(isolver.Cl,vsolver.wData)
+
+    # eof
+    print('')
+
+if __name__ == "__main__":
+    main()

dart/tests/bliTransonic.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.
+
+
+## Compute lifting (linear or nonlinear) viscous flow around a NACA 0012
+# Amaury Bilocq
+#
+# Test the viscous-inviscid interaction scheme
+# Reference to the master's thesis: http://hdl.handle.net/2268/252195
+# Reference test cases with Naca0012 (different from master's thesis):
+# 1) Incompressible: Re = 1e7, M_inf = 0, alpha = 5°, msTE = 0.01, msLE = 0.01
+#    -> nIt = 27, Cl = 0.55, Cd = 0.0058, xtrTop = 0.087, xtrBot = 0.741
+#    -> msLe = 0.001, xtrTop = 0.061
+# 2) Compressible: Re = 1e7, M_inf = 5, alpha = 5°, msTE = 0.01, msLE = 0.01
+#    -> nIt = 33, Cl = 0.65, Cd = 0.0063, xtrTop = 0.057, xtrBot = 0.740
+# 3) Separated: Re = 1e7, M_inf = 0, alpha = 12°, msTE = 0.01, msLE = 0.001
+#    -> nIt = 43, Cl = 1.30 , Cd = 0.0108, xtrTop = 0.010, xtrBot = 1
+#
+# CAUTION
+# This test is provided to ensure that the solver works properly.
+# Mesh refinement may have to be performed to obtain physical results.
+
+import dart.utils as floU
+import dart.default as floD
+
+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.DValidation as validation
+
+import tbox
+import tbox.utils as tboxU
+import fwk
+from fwk.testing import *
+from fwk.coloring import ccolors
+import numpy as np
+
+import math
+
+def main():
+    # timer
+    tms = fwk.Timers()
+    tms['total'].start()
+
+    # define flow variables
+    Re = 1e7
+    alpha = 2*math.pi/180
+    M_inf = 0.73
+
+    user_xtr=[0, 0]
+    user_Ti=None
+
+    mapMesh = 1
+    nFactor = 2
+
+    # Time solver parameters
+    Time_Domain = True # True for unsteady solver, False for steady solver
+    CFL0 = 2
+
+    # ========================================================================================
+
+    U_inf = [math.cos(alpha), math.sin(alpha)] # norm should be 1
+    c_ref = 1
+    dim = 2
+    tol = 1e-4 # tolerance for the VII (usually one drag count)
+    case='Case1FreeTrans.dat' # File containing results from XFOIL of the considered case
+
+    # mesh the geometry
+    print(ccolors.ANSI_BLUE + 'PyMeshing...' + ccolors.ANSI_RESET)
+    tms['msh'].start()
+    pars = {'xLgt' : 5, 'yLgt' : 5, 'msF' : 1, 'msTe' : 0.01, 'msLe' : 0.005}
+    #pars = {'xLgt' : 5, 'yLgt' : 5, 'msF' : 1, 'msTe' : 0.001, 'msLe' : 0.0005}
+    msh, gmshWriter = floD.mesh(dim, 'models/rae2822.geo', pars, ['field', 'airfoil', 'downstream'])
+    tms['msh'].stop()
+
+    # set the problem and add medium, initial/boundary conditions
+    tms['pre'].start()
+    pbl, _, _, bnd, blw = floD.problem(msh, dim, alpha, 0., M_inf, c_ref, c_ref, 0.25, 0., 0., 'airfoil', te = 'te', vsc = True, dbc=True)
+    tms['pre'].stop()
+
+    # solve viscous problem
+
+    print(ccolors.ANSI_BLUE + 'PySolving...' + ccolors.ANSI_RESET)
+    tms['solver'].start()
+    isolver = floD.newton(pbl)
+
+    # Choose between steady and unsteady solver
+    if Time_Domain is True: # Run unsteady solver
+
+        vsolver=floVSMaster.Driver(Re, user_xtr, CFL0, M_inf, case, mapMesh, nFactor)
+        coupler = floVC.Coupler(isolver, vsolver, blw[0], blw[1], tol, gmshWriter, bnd)
+
+    elif Time_Domain is False: # Run steady solver 
+
+        vsolver=floVS_Std.Solver(Re)
+        coupler = floVC_Std.Coupler(isolver, vsolver, blw[0], blw[1], tol, gmshWriter)
+        
+    coupler.run()
+    tms['solver'].stop()
+
+    # extract Cp
+    Cp = floU.extract(bnd.groups[0].tag.elems, isolver.Cp)
+    tboxU.write(Cp, 'Cp_airfoil.dat', '%1.5e', ', ', 'x, y, z, Cp', '')
+    # display results
+    print(ccolors.ANSI_BLUE + 'PyRes...' + ccolors.ANSI_RESET)
+    print('      Re        M    alpha       Cl       Cd      Cdp      Cdf       Cm')
+    print('{0:6.1f}e6 {1:8.2f} {2:8.1f} {3:8.4f} {4:8.4f} {5:8.4f} {6:8.4f} {7:8.4f}'.format(Re/1e6, M_inf, alpha*180/math.pi, isolver.Cl, vsolver.Cd, vsolver.Cdp, vsolver.Cdf, isolver.Cm))
+
+    # display timers
+    tms['total'].stop()
+    print(ccolors.ANSI_BLUE + 'PyTiming...' + ccolors.ANSI_RESET)
+    print('CPU statistics')
+    print(tms)
+    
+    # check results
+    print(ccolors.ANSI_BLUE + 'PyTesting...' + ccolors.ANSI_RESET)
+    try:
+      tests = CTests()
+      if Re == 1e7 and M_inf == 0.715 and alpha == 2*math.pi/180:
+          tests.add(CTest('Cl', isolver.Cl, 0.72, 5e-2))
+          tests.add(CTest('Cd', vsolver.Cd, 0.0061, 1e-3, forceabs=True))
+          tests.add(CTest('Cdp', vsolver.Cdp, 0.0019, 1e-3, forceabs=True))
+          tests.add(CTest('xtr_top', vsolver.xtr[0], 0.34, 0.08, forceabs=True))
+          tests.add(CTest('xtr_bot', vsolver.xtr[1],  0.5, 0.03, forceabs=True))
+      elif Re == 1e7 and M_inf == 0.72 and alpha == 2*math.pi/180:
+          tests.add(CTest('Cl', isolver.Cl, 0.674, 5e-2))
+          tests.add(CTest('Cd', vsolver.Cd, 0.0061, 1e-3, forceabs=True))
+          tests.add(CTest('Cdp', vsolver.Cdp, 0.0019, 1e-3, forceabs=True))
+          tests.add(CTest('xtr_top', vsolver.xtr[0], 0.36, 0.08, forceabs=True))
+          tests.add(CTest('xtr_bot', vsolver.xtr[1],  0.5, 0.03, forceabs=True))
+      elif Re == 1e7 and M_inf == 0.73 and alpha == 2*math.pi/180 and user_xtr==[0, 0]:
+          tests.add(CTest('Cl', isolver.Cl, 0.671, 5e-2))
+          tests.add(CTest('Cd', vsolver.Cd, 0.0089, 1e-3, forceabs=True))
+          tests.add(CTest('Cdp', vsolver.Cdp, 0.0024, 1e-3, forceabs=True))
+          tests.add(CTest('xtr_top', vsolver.xtr[0], 0, 0.08, forceabs=True))
+          tests.add(CTest('xtr_bot', vsolver.xtr[1],  0, 0.03, forceabs=True))
+      else:
+          raise Exception('Test not defined for this flow')
+      tests.run()
+    except Exception as e:
+      print(e)
+
+    # visualize solution and plot results
+    floD.initViewer(pbl)
+    tboxU.plot(Cp[:,0], Cp[:,3], 'x', 'Cp', 'Cl = {0:.{3}f}, Cd = {1:.{3}f}, Cm = {2:.{3}f}'.format(isolver.Cl, vsolver.Cd, isolver.Cm, 4), True)
+
+    val=validation.Validation()
+    val.main(isolver.Cl,vsolver.wData)
+
+    # eof
+    print('')
+
+if __name__ == "__main__":
+    main()

dart/viscous/Drivers/DDriver.py

+################################################################################
+#                                                                              #
+#                            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.Re          = _Re
+        self.Minfty      = _Minfty
+
+        self.xtrF        = user_xtr
+        self.Ti          = None
+        self.mapMesh     = _mapMesh
+        self.nFactor     = _nFactor
+        self.xtr         = [1,1]
+        self.CFL0        = _CFL0
+        self.case        = _case
+
+        self.Cd          = 0
+        self.it          = 0
+        self.uPrevious   = None
+        self.nbrElmUpper = -1
+        self.stagPtMvmt  = 0
+        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 ---------------------------------+')
+            print('- Reynolds number: {:>1.2e}'.format(self.Re))
+            if self.mapMesh:
+                print('- Mesh mapped from {:>5.0f} to {:>5.0f} points.'.format(len(cMeshDict['locCoord']), var.nN))
+            print('- Number of points: {:>5.0f}.'.format(var.nN))
+            if var.xtrF[0] is None:
+                print('- Free transition on the upper side.')
+            else:
+                print('- Forced transition on the upper side; x/c = {:<.4f}.'.format(var.xtrF[0]))
+            if var.xtrF[1] is None:
+                print('- Free transition on the lower side.')
+            else:
+                print('- Forced transition on the lower side; x/c = {:<.4f}.'.format(var.xtrF[1]))
+            print('- Critical amplification ratio: {:<.2f}.'.format(var.Ncrit))
+            if var.xtrF == [0, 0]:
+              print('- Turbulent stagnation point. Ct = 0.03')
+
+            print('')
+            print('Numerical parameters')
+            print('- CFL0: {:<3.2f}'.format(tSolver.integrator.GetCFL0()))
+            print('- Tolerance: {:<3.0f}'.format(np.log10(tSolver.integrator.Gettol())))
+            print('- Solver maximum iterations: {:<5.0f}'.format(tSolver.integrator.GetmaxIt()))
+            print('+------------------------------------------------------------------------------+')
+
+        # Output preprocessing satut.
+
+        print('- Mach max: {:<.2f}.'.format(np.max((var.extPar[0::var.nExtPar]))), end = ' ')
+
+        # Solver setup.
+
+        if self.uPrevious is None or self.nbrElmUpper != var.bNodes[1] or self.stagPtMvmt:
+
+          tSolver.integrator.itFrozenJac0 = 1
+          tSolver.initSln = 1   # Flag to use time solver initial condition.
+          tSolver.integrator.SetCFL0(1) # Low starting CFL.
+          print('Restart sln.')
+
+          if self.it != 0 and self.nbrElmUpper != var.bNodes[1]:
+
+            self.stagPtMvmt = 1
+            print(ccolors.ANSI_BLUE, '- Stagnation point moved.',ccolors.ANSI_RESET)
+
+          else:
+
+            self.stagPtMvmt = 0
+
+        else: # If we use a solution previously obtained. 
+
+          var.u = self.uPrevious.copy()   # Use previous solution.
+
+          print('Updating sln.')
+
+        self.nbrElmUpper = var.bNodes[1]
+        DirichletBC.airfoilBC(var)    # Reset boundary condition.
+
+        # Run the time integration
+        convergedPoints = tSolver.Run(var)
+        # Output time solver status.
+        if np.any(convergedPoints != 0):
+            print(ccolors.ANSI_YELLOW + 'Some point(s) did not converge.', ccolors.ANSI_RESET)
+
+        self.uPrevious=var.u.copy()   # Save solution to use it as initial condition the next iteration.
+
+        # 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 >= 24:
+          self.writeFile()
+          Validation().main(1,self.wData)"""
+
+        pass
\ No newline at end of file

dart/viscous/Drivers/DGroupConstructor.py

+################################################################################
+#                                                                              #
+#                             Flow viscous module file                         #
+#                        Constructor : Merges and sorts data                   #
+#                   of the 3 groups (upper, lower sides and wake)              #
+#                                                                              #
+# Author: Paul Dechamps                                                        #
+# Université de Liège                                                          #
+# Date: 2021.09.10                                                             #
+#                                                                              #
+################################################################################
+
+
+# ------------------------------------------------------------------------------
+#  Imports
+# ------------------------------------------------------------------------------
+from matplotlib import pyplot as plt
+import numpy as np
+
+class GroupConstructor():
+    """ ---------- Group Constructor ----------
+        
+        Builds single vectors elements out of the 3 groups given in input.
+    """
+    def __init__(self) -> None:
+        pass
+
+    def mergeVectors(self,dataIn, mapMesh, nFactor):
+        """Merges 3 groups upper/lower sides and wake.
+        """
+        
+        for k in range(len(dataIn)):
+
+            if len(dataIn[k]) == 2: # Airfoil case.
+            
+                xx_up, dv_up, vtExt_up, _, alpha_up = self.__getBLcoordinates(dataIn[k][0])
+                xx_lw, dv_lw, vtExt_lw, _, alpha_lw = self.__getBLcoordinates(dataIn[k][1])
+
+            else: # Wake case
+
+                xx_wk, dv_wk, vtExt_wk, _, alpha_wk = self.__getBLcoordinates(dataIn[k])
+
+        if mapMesh:
+            # Save initial mesh.
+            cMesh        = np.concatenate((dataIn[0][0][:,0], dataIn[0][1][1:,0], dataIn[1][:,0]))
+            cMeshbNodes  = [0, len(dataIn[0][0][:,0]), len(dataIn[0][0][:,0]) + len(dataIn[0][1][1:,0])]
+            cxx          = np.concatenate((xx_up, xx_lw[1:], xx_wk))
+            cvtExt       = np.concatenate((vtExt_up, vtExt_lw[1:], vtExt_wk))
+            cxxExt       = np.concatenate((dataIn[0][0][:,8], dataIn[0][1][1:,8],dataIn[1][:,8]))
+
+            # Create fine mesh.
+            fMeshUp      = self.createFineMesh(dataIn[0][0][:,0], nFactor)
+            fMeshLw      = self.createFineMesh(dataIn[0][1][:,0], nFactor)
+            fMeshWk      = self.createFineMesh(dataIn[1][:,0]   , nFactor)
+            fMesh        = np.concatenate((fMeshUp, fMeshLw[1:], fMeshWk))
+            fMeshbNodes  = [0, len(fMeshUp), len(fMeshUp) + len(fMeshLw[1:])]
+
+            fxx_up       = self.interpolateToFineMesh(xx_up, fMeshUp, nFactor)
+            fxx_lw       = self.interpolateToFineMesh(xx_lw, fMeshLw, nFactor)
+            fxx_wk       = self.interpolateToFineMesh(xx_wk, fMeshWk, nFactor)
+            fxx          = np.concatenate((fxx_up, fxx_lw[1:], fxx_wk))
+
+            fvtExt_up    = self.interpolateToFineMesh(vtExt_up, fMeshUp, nFactor)
+            fvtExt_lw    = self.interpolateToFineMesh(vtExt_lw, fMeshLw, nFactor)
+            fvtExt_wk    = self.interpolateToFineMesh(vtExt_wk, fMeshWk, nFactor)
+            fvtExt       = np.concatenate((fvtExt_up, fvtExt_lw[1:], fvtExt_wk))
+
+            fMe_up       = self.interpolateToFineMesh(dataIn[0][0][:,5], fMeshUp, nFactor)
+            fMe_lw       = self.interpolateToFineMesh(dataIn[0][1][:,5], fMeshLw, nFactor)
+            fMe_wk       = self.interpolateToFineMesh(dataIn[1][:,5]   , fMeshWk, nFactor)
+            fMe          = np.concatenate((fMe_up, fMe_lw[1:], fMe_wk))
+
+            frho_up      = self.interpolateToFineMesh(dataIn[0][0][:,6], fMeshUp, nFactor)
+            frho_lw      = self.interpolateToFineMesh(dataIn[0][1][:,6], fMeshLw, nFactor)
+            frho_wk      = self.interpolateToFineMesh(dataIn[1][:,6]   , fMeshWk, nFactor)
+            frho         = np.concatenate((frho_up, frho_lw[1:], frho_wk))
+
+            fdStarExt_up = self.interpolateToFineMesh(dataIn[0][0][:,7], fMeshUp, nFactor)
+            fdStarExt_lw = self.interpolateToFineMesh(dataIn[0][1][:,7], fMeshLw, nFactor)
+            fdStarExt_wk = self.interpolateToFineMesh(dataIn[1][:,7]   , fMeshWk, nFactor)
+            
+            fxxExt_up    = self.interpolateToFineMesh(dataIn[0][0][:,8], fMeshUp, nFactor)
+            fxxExt_lw    = self.interpolateToFineMesh(dataIn[0][1][:,8], fMeshLw, nFactor)
+            fxxExt_wk    = self.interpolateToFineMesh(dataIn[1][:,8]   , fMeshWk, nFactor)
+            
+            fdStarext    = np.concatenate((fdStarExt_up, fdStarExt_lw[1:], fdStarExt_wk))
+            fxxext       = np.concatenate((fxxExt_up, fxxExt_lw[1:], fxxExt_wk))
+
+            fdv          = np.concatenate((dv_up, dv_lw[1:], dv_wk))
+
+            fMeshDict    = {'globCoord': fMesh, 'bNodes': fMeshbNodes, 'locCoord': fxx,
+                            'vtExt': fvtExt, 'Me': fMe, 'rho': frho, 'deltaStarExt': fdStarext, 'xxExt': fxxext, 'dv': fdv}
+
+            cMe          = np.concatenate((dataIn[0][0][:,5], dataIn[0][1][1:,5], dataIn[1][:,5]))
+            cRho         = np.concatenate((dataIn[0][0][:,6], dataIn[0][1][1:,6], dataIn[1][:,6]))
+            cdStarExt    = np.concatenate((dataIn[0][0][:,7], dataIn[0][1][1:,7], dataIn[1][:,7]))
+            cMeshDict    = {'globCoord': cMesh, 'bNodes': cMeshbNodes, 'locCoord': cxx,
+                            'vtExt': cvtExt, 'Me': cMe, 'rho': cRho, 'deltaStarExt': cdStarExt, 'xxExt': cxxExt, 'dv': fdv}
+
+            dataUpper = np.empty((len(fMeshUp), 0))
+            dataLower = np.empty((len(fMeshLw), 0))
+            dataWake  = np.empty((len(fMeshWk), 0))
+            for iParam in range(len(dataIn[0][0][0,:])):
+                interpParamUp = self.interpolateToFineMesh(dataIn[0][0][:,iParam], fMeshUp, nFactor)
+                dataUpper     = np.column_stack((dataUpper, interpParamUp))
+                interpParamLw = self.interpolateToFineMesh(dataIn[0][1][:,iParam], fMeshLw, nFactor)
+                dataLower     = np.column_stack((dataLower, interpParamLw))
+                interpParamWk = self.interpolateToFineMesh(dataIn[1][:,iParam]   , fMeshWk, nFactor)
+                dataWake      = np.column_stack((dataWake, interpParamWk))
+            
+            dataLower = np.delete(dataLower, 0, axis = 0) # Remove stagnation point doublon. 
+
+        else:
+
+            Mesh         = np.concatenate((dataIn[0][0][:,0], dataIn[0][1][1:,0], dataIn[1][:,0]))
+            MeshbNodes   = [0, len(dataIn[0][0][:,0]), len(dataIn[0][0][:,0]) + len(dataIn[0][1][1:,0])]
+
+            xx           = np.concatenate((xx_up, xx_lw[1:], xx_wk))
+            vtExt        = np.concatenate((vtExt_up, vtExt_lw[1:], vtExt_wk))
+
+            alpha        = np.concatenate((alpha_up, alpha_lw, alpha_wk))
+            dv           = np.concatenate((dv_up, dv_lw[1:], dv_wk))
+
+            Me           = np.concatenate((dataIn[0][0][:,5], dataIn[0][1][1:,5], dataIn[1][:,5]))
+            rho          = np.concatenate((dataIn[0][0][:,6], dataIn[0][1][1:,6], dataIn[1][:,6]))
+            dStarext     = np.concatenate((dataIn[0][0][:,7], dataIn[0][1][1:,7], dataIn[1][:,7]))
+            xxExt        = np.concatenate((dataIn[0][0][:,8], dataIn[0][1][1:,8], dataIn[1][:,8]))
+
+            fMeshDict    = {'globCoord': Mesh, 'bNodes': MeshbNodes, 'locCoord': xx,
+                            'vtExt': vtExt, 'Me': Me, 'rho': rho, 'deltaStarExt': dStarext, 'xxExt': xxExt, 'dv': dv}
+
+            cMeshDict    = fMeshDict.copy()
+
+            dataUpper    = dataIn[0][0]
+            dataLower    = dataIn[0][1][1:, :]
+            dataWake     = dataIn[1]
+        
+        data = np.concatenate((dataUpper, dataLower, dataWake), axis = 0)
+
+        return fMeshDict, cMeshDict, data
+
+    def __getBLcoordinates(self,data):
+        """Transform the reference coordinates into airfoil coordinates
+        """
+        nN    = len(data[:,0])
+        nE    = nN-1            # Nbr of element along the surface.
+        vt    = np.zeros(nN)    # Outer flow velocity. 
+        xx    = np.zeros(nN)    # Position along the surface of the airfoil.
+        dx    = np.zeros(nE)    # Distance along two nodes.
+        dv    = np.zeros(nE)    # Speed difference according to element.
+        alpha = np.zeros(nE)    # Angle of the element for the rotation matrix.
+
+        for i in range(0, nE):
+
+            alpha[i] = np.arctan2((data[i+1, 1] - data[i, 1]), (data[i+1, 0]-data[i, 0]))
+            dx[i]    = np.sqrt((data[i+1, 0] - data[i,0])**2 + (data[i+1, 1]-data[i, 1])**2)
+            xx[i+1]  = dx[i] + xx[i]
+            vt[i]    = (data[i, 3] * np.cos(alpha[i]) + data[i ,4] * np.sin(alpha[i]))          # Tangential speed at node 1 element i
+            vt[i+1]  = (data[i+1, 3] * np.cos(alpha[i]) + data[i+1, 4] * np.sin(alpha[i]))      # Tangential speed at node 2 element i
+            dv[i]    = (vt[i+1] - vt[i]) / dx[i]                                                # Velocity gradient according to element i.
+                                                                                                #             central scheme with half point
+        
+        return xx, dv, vt, nN, alpha
+
+    def interpolateToFineMesh(self, cVector, fMesh, nFactor):
+        
+        fVector = np.empty(len(fMesh)-1)
+        for cPoint in range(len(cVector) - 1):
+            dv = cVector[cPoint + 1] - cVector[cPoint]
+
+            for fPoint in range(nFactor):
+                fVector[nFactor * cPoint + fPoint] = cVector[cPoint] + fPoint * dv / nFactor
+        
+        fVector = np.append(fVector, cVector[-1])
+
+        return fVector
+
+    def createFineMesh(self, cMesh, nFactor):
+
+        fMesh = []
+        for iPoint in range(len(cMesh) - 1):
+
+            dx = cMesh[iPoint + 1] - cMesh[iPoint]
+
+            for iRefinedPoint in range(nFactor):
+
+                fMesh = np.append(fMesh, cMesh[iPoint] + iRefinedPoint * dx / nFactor)
+
+        fMesh = np.append(fMesh, cMesh[-1])
+        return fMesh
\ No newline at end of file

dart/viscous/Drivers/DPostProcess.py

+import numpy as np
+
+class PostProcess:
+  def __init__(self):
+    pass
+
+  def sortParam(self,var, mapMesh, cMeshDict, nFactor):
+
+    # Post process delta star
+
+    # Recompute deltaStar.
+    var.blPar[9::var.nBlPar] = var.u[0::var.nVar] * var.u[1::var.nVar]
+
+    # Map delta star to the coarse mesh and divide between the airfoil
+    # and wake related points
+    if mapMesh:
+
+      deltaStar      = self.inverseMeshMapping(var.blPar[9::var.nBlPar], var.bNodes, cMeshDict['globCoord'], cMeshDict['bNodes'], nFactor)
+
+      AFdeltaStarOut = deltaStar[:cMeshDict['bNodes'][2]]
+      WKdeltaStarOut = deltaStar[cMeshDict['bNodes'][2]:]
+
+    else:
+      
+      AFdeltaStarOut = var.blPar[9:var.bNodes[2]*var.nBlPar:var.nBlPar]
+      WKdeltaStarOut = var.blPar[var.bNodes[2]*var.nBlPar + 9::var.nBlPar]
+
+    # Compute blowing velocity on each cell.
+
+    blwVelUp = np.zeros(var.bNodes[1] - 1)
+    blwVelLw = np.zeros(var.bNodes[2] - var.bNodes[1])
+    blwVelWk = np.zeros(var.nN - var.bNodes[2] - 1)
+
+    # Upper
+    i = 0
+    for iPoint in range(1, var.bNodes[1]):
+
+      iPrev = iPoint - 1
+
+      blwVelUp[i] = (var.extPar[iPoint*var.nExtPar+1] * var.u[iPoint*var.nVar+3] * var.blPar[iPoint*var.nBlPar+9]
+              - var.extPar[iPrev*var.nExtPar+1] * var.u[iPrev*var.nVar+3] * var.blPar[iPrev*var.nBlPar+9]) / (var.extPar[iPoint*var.nExtPar+1] * (var.xx[iPoint] - var.xx[iPrev]))
+
+      i += 1
+
+    # Lower
+    i = 0
+    for iPoint in range(var.bNodes[1], var.bNodes[2]):
+
+      iPrev = 0 if iPoint == var.bNodes[1] else iPoint - 1
+
+      blwVelLw[i] = (var.extPar[iPoint*var.nExtPar+1] * var.u[iPoint*var.nVar+3] * var.blPar[iPoint*var.nBlPar+9]
+              - var.extPar[iPrev*var.nExtPar+1] * var.u[iPrev*var.nVar+3] * var.blPar[iPrev*var.nBlPar+9]) / (var.extPar[iPoint*var.nExtPar+1] * (var.xx[iPoint] - var.xx[iPrev]))
+
+      i += 1
+
+    # Wake
+    i = 0
+    for iPoint in range(var.bNodes[2] + 1, var.nN):
+
+      iPrev = iPoint - 1
+
+      blwVelWk[i] = (var.extPar[iPoint*var.nExtPar+1] * var.u[iPoint*var.nVar+3] * var.blPar[iPoint*var.nBlPar+9]
+              - var.extPar[iPrev*var.nExtPar+1] * var.u[iPrev*var.nVar+3] * var.blPar[iPrev*var.nBlPar+9]) / (var.extPar[iPoint*var.nExtPar+1] * (var.xx[iPoint] - var.xx[iPrev]))
+
+      i += 1
+
+    if mapMesh:
+      blwVelUp, blwVelLw, blwVelWk  = self.mapBlowingVelocities(blwVelUp, blwVelLw, blwVelWk, var.bNodes, cMeshDict['globCoord'], cMeshDict['bNodes'], var.xGlobal, nFactor)
+
+    # Postprocess the boundary layer parameters.
+
+    # Normalize friction and dissipation coefficients.
+    frictionCoeff=var.u[3::var.nVar]**2 * var.blPar[2::var.nBlPar]
+    dissipCoeff=var.u[3::var.nVar]**2 * var.blPar[1::var.nBlPar]
+
+    # Compute total drag coefficient.
+    CdTot=2*var.u[-5]*(var.u[-2]**((var.u[-4]+5)/2))
+
+    # Compute friction and pressure drag coefficients. 
+    Cdf_upper=np.trapz(frictionCoeff[:var.bNodes[1]],var.xGlobal[:var.bNodes[1]])
+    # Insert stagnation point in the lower side. 
+    Cdf_lower=np.trapz(np.insert(frictionCoeff[var.bNodes[1]:var.bNodes[2]],0,frictionCoeff[0]),
+               np.insert(var.xGlobal[var.bNodes[1]:var.bNodes[2]],0,var.xGlobal[0]))
+    Cdf=Cdf_upper+Cdf_lower # No friction in the wake
+    Cdp=CdTot-Cdf
+
+    # [CdTot, Cdf, Cdp]
+    blScalars=[CdTot, Cdf, Cdp]
+
+    # Store boundary layer final parameters in lists. 
+    # [deltaStar, theta, H, Hk, Hstar, Hstar2, Cf, Cd, Ct, delta]
+    blVectors_airfoil=[var.blPar[9:var.bNodes[2]*var.nBlPar:var.nBlPar],
+              var.u[0:var.bNodes[2]*var.nVar:var.nVar],
+              var.u[1:var.bNodes[2]*var.nVar:var.nVar],
+              var.blPar[6:var.bNodes[2]*var.nBlPar:var.nBlPar],
+              var.blPar[4:var.bNodes[2]*var.nBlPar:var.nBlPar],
+              var.blPar[5:var.bNodes[2]*var.nBlPar:var.nBlPar],
+              frictionCoeff[:var.bNodes[2]], dissipCoeff[:var.bNodes[2]],
+              var.u[4:var.bNodes[2]*var.nVar:var.nVar],
+              var.blPar[3:var.bNodes[2]*var.nBlPar:var.nBlPar]]
+
+    blVectors_wake=[var.blPar[var.bNodes[2]*var.nBlPar+9::var.nBlPar],
+            var.u[var.bNodes[2]*var.nVar+0::var.nVar],
+            var.u[var.bNodes[2]*var.nVar+1::var.nVar],
+            var.blPar[var.bNodes[2]*var.nBlPar+6::var.nBlPar],
+            var.blPar[var.bNodes[2]*var.nBlPar+4::var.nBlPar],
+            var.blPar[var.bNodes[2]*var.nBlPar+5::var.nBlPar],
+            frictionCoeff[var.bNodes[2]:], dissipCoeff[var.bNodes[2]:],
+            var.u[var.bNodes[2]*var.nVar+4::var.nVar],
+            var.blPar[var.bNodes[2]*var.nBlPar+3::var.nBlPar]]
+
+    return blScalars, blwVelUp, blwVelLw, blwVelWk, blVectors_airfoil, blVectors_wake, AFdeltaStarOut, WKdeltaStarOut
+
+  def mapBlowingVelocities(self, blwVelUp, blwVelLw, blwVelWk, fbNodes, cMesh, cMeshbNodes, fMesh, nFactor):
+    """Maps blowing velocties from cell to cell using weighted average interpolation.
+
+        Args
+        ----
+
+        - blwVelUp (numpy.ndarray): Blowing velocities on the fine mesh of the upper side.
+
+        - blwVelLw (numpy.ndarray): Blowing velocities on the fine mesh of the lower side.
+    
+        - blwVelWk (numpy.ndarray): Blowing velocities on the fine mesh of the wake.
+
+        - fbNodes (numpy.ndarray): Boundary nodes of the fine mesh.
+
+        - cMesh (numpy.ndarray): Coarse mesh coordinates in the global frame of reference.
+
+        - cMeshbNodes (numpy.ndarray): Boundary nodes of the coarse mesh.
+
+        - fMesh (numpy.ndarray): Fine mesh coordinates in the global frame of reference.
+
+        - nFactor (int): Multiplicator underlying mesh adaptation.
+
+    Example
+    -------
+
+    Fine mesh:   |-----------|-----------|-----------|----------|
+               \	       /
+           (1/2)	\	  	  / (1/2)
+               \       /
+                \     /
+    Coarse mesh: |-----------------------|----------------------|
+    """
+
+    cblwVelUp = np.empty(len(cMesh[:cMeshbNodes[1]]) - 1)
+    cblwVelLw = np.empty(len(cMesh[cMeshbNodes[1]:cMeshbNodes[2]]))
+    cWKblwVel = np.empty(len(cMesh[cMeshbNodes[2]:]) - 1)
+
+    for cCell in range(len(cblwVelUp)):
+
+      cblwVelUp[cCell] = np.mean(blwVelUp[cCell*nFactor:(cCell + 1)*nFactor])
+    
+    for cCell in range(len(cblwVelLw)):
+
+      cblwVelLw[cCell] = np.mean(blwVelLw[cCell*nFactor:(cCell + 1)*nFactor])
+
+    for cCell in range(len(cWKblwVel)):
+
+      cWKblwVel[cCell] = np.mean(blwVelWk[cCell*nFactor:(cCell + 1)*nFactor])
+
+    return cblwVelUp, cblwVelLw, cWKblwVel
+
+
+  def inverseMeshMapping(self, fVector, fbNodes, cMesh, cbNodes, nFactor):
+    """Maps the variables required for blowing velocity computation from
+    the fine mesh to the coarse mesh. Group separation is performed
+    to avoid confusion between the points at the same place on the mesh
+    but with different physical meaning (e.g First wake point).
+
+    The mapping is based on a weigthed average with the nearest neighbors.
+
+    Args
+    ----
+
+    - fVector (numpy.ndarray): Vector containing the function on the fine mesh.
+
+    - fbNodes (numpy.ndarray): Id of the boundary nodes on the fine mesh.
+
+    - cMesh (numpy.ndarray): Coarse mesh.
+
+    - cbNodes (numpy.ndarray): Id of the boundary nodes on the coarse mesh.
+
+    - nFactor (int): Multiplicator underlying mesh adaptation.
+
+    Return
+    ------
+
+    - cVector (numpy.ndarray): Function contained in 'fVector' mapped to the coarse mesh.
+
+    Example
+    -------
+
+    Fine mesh:   |-----------|-----------|-----------|----------|
+                  \	         |          /
+              (1/4)  \	     | (1/2)   / (1/4)
+                  \        V        /
+                    ----->   <-----
+    Coarse mesh: |-----------------------|----------------------|
+    """
+
+    cVector = np.zeros(len(cMesh))
+
+    for cPoint in range(len(cMesh)-1):
+      cVector[cPoint] = fVector[cPoint * nFactor]
+    cVector[-1] = fVector[-1]
+    return cVector

dart/viscous/Master/DAirfoil.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.
+
+
+## Airfoil around which the boundary layer is computed
+# Amaury Bilocq
+
+from dart.viscous.Master.DBoundary import Boundary
+
+import numpy as np
+
+class Airfoil(Boundary):
+    def __init__(self, _boundary):
+        Boundary.__init__(self, _boundary)
+        self.name = 'airfoil' # type of boundary
+        self.T0 = 0 # initial condition for the momentum thickness
+        self.H0 = 0
+        self.n0 = 0
+        self.Ct0 = 0
+        self.TEnd = [0,0]
+        self.HEnd = [0,0]
+        self.CtEnd = [0,0]
+        self.nEnd = [0,0]
+    
+    def initialConditions(self, Re, dv):
+        if dv > 0:
+            self.T0 = np.sqrt(0.075/(Re*dv))
+        else:
+            self.T0 = 1e-6
+        self.H0 = 2.23 # One parameter family Falkner Skan
+        self.n0 = 0
+        self.Ct0 = 0
+        return self.T0, self.H0, self.n0, self.Ct0
+
+
+    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
+        # Find the stagnation point
+        """if it == 0:
+          self.idxStag = np.argmin(np.sqrt(data[:,4]**2+data[:,5]**2))
+          idxStag = self.idxStag
+        else:
+          idxStag = self.idxStag"""
+        idxStag = np.argmin(np.sqrt(data[:,4]**2+data[:,5]**2))
+        globStag = int(data[idxStag,0]) # position of the stagnation point in boundary.nodes
+        # Find TE nodes (assuming suction side element is above pressure side element in the y direction)
+        idxTE = np.where(data[:,1] == np.max(data[:,1]))[0] # TE nodes
+        idxTEe0 = np.where(np.logical_or(eData[:,1] == data[idxTE[0],0], eData[:,2] == data[idxTE[0],0]))[0] # TE element 1
+        idxTEe1 = np.where(np.logical_or(eData[:,1] == data[idxTE[1],0], eData[:,2] == data[idxTE[1],0]))[0] # TE element 2
+        ye0 = 0.5 * (data[np.where(data[:,0] == eData[idxTEe0[0],1])[0],2] + data[np.where(data[:,0] == eData[idxTEe0[0],2])[0],2]) # y coordinates of TE element 1 CG
+        ye1 = 0.5 * (data[np.where(data[:,0] == eData[idxTEe1[0],1])[0],2] + data[np.where(data[:,0] == eData[idxTEe1[0],2])[0],2]) # y coordinates of TE element 2 CG
+        if ye0 - ye1 > 0: # element 1 containing TE node 1 is above element 2
+            upperGlobTE = int(data[idxTE[0],0])
+            lowerGlobTE = int(data[idxTE[1],0])
+        else:
+             upperGlobTE = int(data[idxTE[1],0])
+             lowerGlobTE = int(data[idxTE[0],0])
+        # Connectivity
+        connectListElems[0] = np.where(eData[:,1] == globStag)[0]
+        N1[0] = eData[connectListElems[0],1] # number of the stag node elems.nodes -> globStag
+        connectListNodes[0] = np.where(data[:,0] == N1[0])[0]
+        i = 1
+        upperTE = 0
+        lowerTE = 0
+        # Sort the suction part
+        while upperTE == 0:
+            N1[i] = eData[connectListElems[i-1],2] # Second node of the element
+            connectListElems[i] = np.where(eData[:,1] == N1[i])[0] # # Index of the first node of the next element in elems.nodes
+            connectListNodes[i] = np.where(data[:,0] == N1[i])[0] # Index of the node in boundary.nodes
+            if eData[connectListElems[i],2] == upperGlobTE:
+                upperTE = 1
+            i += 1
+        # Sort the pressure side
+        connectListElems[i] = np.where(eData[:,2] == globStag)[0]
+        connectListNodes[i] = np.where(data[:,0] == upperGlobTE)[0]
+        N1[i] = eData[connectListElems[i],2]
+        while lowerTE == 0:
+            N1[i+1]  = eData[connectListElems[i],1] # First node of the element
+            connectListElems[i+1] = np.where(eData[:,2] == N1[i+1])[0] # # Index of the second node of the next element in elems.nodes
+            connectListNodes[i+1] = np.where(data[:,0] == N1[i+1])[0] # Index of the node in boundary.nodes
+            if eData[connectListElems[i+1],1] == lowerGlobTE:
+                lowerTE = 1
+            i += 1
+        connectListNodes[i+1] = np.where(data[:,0] == lowerGlobTE)[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)
+        uData = data[0:np.argmax(data[:,0])+1]
+        lData = data[np.argmax(data[:,0])+1:None]
+        lData = np.insert(lData, 0, uData[0,:], axis = 0) #double the stagnation point
+        return connectListNodes, connectListElems,[uData, lData]

dart/viscous/Master/DCoupler.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.
+
+
+## Viscous-inviscid coupler (quasi-simultaneous coupling)
+# Amaury Bilocq
+
+from dart.viscous.Master.DAirfoil import Airfoil
+from dart.viscous.Master.DWake import Wake
+from matplotlib import pyplot as plt
+
+import numpy as np
+
+import tbox.utils as tboxU
+import dart.utils as floU
+from fwk.coloring import ccolors
+
+class Coupler:
+    def __init__(self, _isolver, _vsolver, _boundaryAirfoil, _boundaryWake, _tol, _writer, _bnd):
+        self.bnd=_bnd
+        self.isolver =_isolver # inviscid solver
+        self.vsolver = _vsolver # viscous solver
+        
+        self.group = [Airfoil(_boundaryAirfoil), Wake(_boundaryWake)] # airfoil and wake python objects
+        self.tol = _tol # tolerance of the VII
+        self.writer = _writer
+
+    def run(self):
+        """Viscous inviscid coupling.
+        """
+
+        # Save history.
+        cvgMonitoring = []
+
+        # Initialize loop
+        it = 0
+        converged = 0 # temp
+        CdOld = self.vsolver.Cd
+
+        print('+---------------------------- VII Solver begin --------------------------------+')
+        print('')
+
+        while converged == 0:
+
+            # Run inviscid solver
+            print('+----------------------------- Inviscid solver --------------------------------+')
+            self.isolver.reset()
+            self.isolver.run()  
+
+            for n in range(0, len(self.group)):
+
+                for i in range(0, len(self.group[n].boundary.nodes)):
+                    
+                    # Extract outer flow surface velocities.
+                    self.group[n].v[i,0] = self.isolver.U[self.group[n].boundary.nodes[i].row][0]
+                    self.group[n].v[i,1] = self.isolver.U[self.group[n].boundary.nodes[i].row][1]
+                    self.group[n].v[i,2] = self.isolver.U[self.group[n].boundary.nodes[i].row][2]
+                    self.group[n].Me[i] = self.isolver.M[self.group[n].boundary.nodes[i].row]
+                    self.group[n].rhoe[i] = self.isolver.rho[self.group[n].boundary.nodes[i].row]
+                
+            if it==0:
+                # Write inviscid pressure distribution on the first iteration
+                Cp = floU.extract(self.bnd.groups[0].tag.elems, self.isolver.Cp)
+                tboxU.write(Cp, 'Cpinv_airfoil.dat', '%1.5e', ', ', 'x, y, z, Cp', '')
+
+            # Run viscous solver
+            print('+------------------------------ Viscous solver --------------------------------+')
+            self.vsolver.run(self.group)
+
+            for n in range(0, len(self.group)):
+                for i in range(0, self.group[n].nE):
+                    self.group[n].boundary.setU(i, self.group[n].u[i])
+
+            error = abs(self.vsolver.Cd - CdOld)/self.vsolver.Cd
+            cvgMonitoring.append(error)
+
+            # Check convergence and output coupling iteration.
+
+            if error <= self.tol:
+              print('')
+              print(ccolors.ANSI_GREEN, 'It: {:<3.0f} Cd = {:<6.5f}. Error = {:<2.3f} < {:<2.3f}'
+              .format(it, self.vsolver.Cd, np.log10(error), np.log10(self.tol)), ccolors.ANSI_RESET)
+              print('')
+              """plt.plot(cvgMonitoring[:it+1])
+              plt.yscale('log')
+              plt.xlabel('Number of iterations (-)')
+              plt.ylabel('log[Res] (-)')
+              plt.show()"""
+
+              converged = 1
+            else:
+              print('')
+              print(ccolors.ANSI_BLUE, 'It: {:<3.0f} Cd = {:<6.5f}. Error = {:<2.3f} > {:<2.3f}'
+              .format(it, self.vsolver.Cd, np.log10(error), np.log10(self.tol)), ccolors.ANSI_RESET)
+              print('')
+
+              CdOld = self.vsolver.Cd
+
+            # Prepare next iteration.
+            it += 1
+            self.vsolver.it += 1
+
+        # Save results
+
+        self.isolver.save(self.writer)    # Inviscid solver files.
+        self.vsolver.writeFile()          # Viscous solver files.
+        Cp = floU.extract(self.bnd.groups[0].tag.elems, self.isolver.Cp)    # Pressure coefficient.
+        tboxU.write(Cp, 'Cp_airfoil.dat', '%1.5e', ', ', 'x, y, z, Cp', '')
+        print('\n')

dart/viscous/wake.py → dart/viscous/Master/DWake.py

@@ -19,8 +19,8 @@
 ## Wake behind airfoil (around which the boundary layer is computed)
 # Amaury Bilocq
 
-from dart.viscous.boundary import Boundary
-from dart.viscous.airfoil import Airfoil
+from dart.viscous.Master.DBoundary import Boundary
+
 import numpy as np
 
 class Wake(Boundary):

dart/viscous/Physics/DBIConditions.py

+################################################################################
+#                                                                              #
+#                            DARTFlo viscous module file                       #
+#                    BIConditions: Initial and boundary conditions             #
+#                        of the pseudo-time marching problem
+#                                                                              #
+# Author: Paul Dechamps                                                        #
+# Université de Liège                                                          #
+# Date: 2021.09.10                                                             #
+#                                                                              #
+################################################################################
+
+
+# ------------------------------------------------------------------------------
+#  Imports
+# ------------------------------------------------------------------------------
+from dart.viscous.Physics.DClosures import Closures
+
+import matplotlib.pyplot as plt 
+import numpy as np
+
+
+class Initial:
+    def __init__(self, BoundaryCond):
+        self.boundCond=BoundaryCond
+        self.mu_air=1.849e-5 # Dynamic viscosity of air at 25°C
+
+    def initialConditions(self, var):
+        """Twaithes integral solution of the (steady) Von Karman momentum
+        integral equation for a laminar flow over a flat plate
+
+        Parameters
+        ----------
+        - var: class
+             Data container: - Solution
+                             - Boundary layer parameters
+                             - External flow parameters
+                             - Mesh
+                             - Transition information
+        
+        Tasks
+        -----
+        - Compute Twaithes solution of the momentum integral equation for a laminar flow over a flat plate in the
+        regime prescibed by the inviscid solver. 
+        - Define boundary conditions at the stagnation point and the first wake point 
+        """
+
+        
+        u0=np.zeros(var.nN*var.nVar)        # Vector containing the intial solution
+        
+        # Loop over the three groups (upper, lower sides, wake)
+        for n in range(3):
+
+            if n==2:
+
+                xx=var.xx[var.bNodes[n]:]                                         # Airfoil local coordinates.
+                vt=var.extPar[var.bNodes[n]*var.nExtPar+2::var.nExtPar]           # Inviscid velocities.
+                rho=var.extPar[var.bNodes[n]*var.nExtPar+1::var.nExtPar]          # Air density in each cell. 
+
+            else:
+
+                xx=var.xx[var.bNodes[n]:var.bNodes[n+1]]
+                vt=var.extPar[var.bNodes[n]*var.nExtPar+2:var.bNodes[n+1]*var.nExtPar:var.nExtPar]
+                rho=var.extPar[var.bNodes[n]*var.nExtPar+1:var.bNodes[n+1]*var.nExtPar:var.nExtPar]
+
+            # Reynolds number based on the arc length.
+            Rex=rho*vt*xx/self.mu_air
+            Rex[Rex==0]=1 # Stagnation point
+            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.
+
+            for iMarker in range(len(lambdaTw)):
+
+                idxPrev = iMarker-1
+                lambdaTw[iMarker] = (rho[iMarker] * ThetaTw[iMarker]**2) / self.mu_air * (vt[iMarker] - vt[idxPrev])#/(xx[idx]-xx[idxPrev])
+
+                # Empirical correlations for H=f(lambda)
+                if 0 <= lambdaTw[iMarker] < 0.1:
+
+                    HTw[iMarker] = 2.61 - 3.75*lambdaTw[iMarker] + 5.24*lambdaTw[iMarker]**2
+
+                elif -0.1<lambdaTw[iMarker]<0:
+
+                    HTw[iMarker] = 2.088 + 0.0731 / (lambdaTw[iMarker] + 0.14)
+            
+            if n==2:
+
+                u0[var.bNodes[n]*var.nVar+0::var.nVar]=ThetaTw      # Theta
+                u0[var.bNodes[n]*var.nVar+1::var.nVar]=HTw          # H
+
+            else:
+
+                u0[var.bNodes[n]*var.nVar+0:var.bNodes[n+1]*var.nVar:var.nVar]=ThetaTw      # Theta
+                u0[var.bNodes[n]*var.nVar+1:var.bNodes[n+1]*var.nVar:var.nVar]=HTw          # H
+
+        u0[1::var.nVar]=2                           # H
+        u0[2::var.nVar]=0                           # N
+        u0[3::var.nVar]=var.extPar[2::var.nExtPar]  # Vt
+        u0[4::var.nVar]=0.001                       # Ct
+        
+        # Set initial condition. 
+        var.u=u0
+
+        # Impose boudary conditions.
+        self.boundCond.airfoilBC(var)
+        self.boundCond.wakeBC(var)
+        pass
+
+
+class Boundaries:
+    """ ---------- Boundary conditions ----------
+        
+        Boundary conditions at the stagnation point and on the first point of the wake.
+
+        Attributes
+        ----------
+        Re: Flow Reynolds number.
+        """
+
+    def __init__(self, _Re):
+
+        self.Re = _Re
+        
+    def airfoilBC(self, cfgFlow):
+        """Boundary conditions at the stagnation point.
+        Twaithes solution method values used for the
+        stagnation point.
+        """
+
+        if cfgFlow.dv[0] > 0:
+            T0 = np.sqrt(0.075 / (self.Re * cfgFlow.dv[0]))
+        else:
+            T0 = 1e-6
+
+        H0  = 2.23                # One parameter family Falkner Skan.
+        n0  = 0                   # Initial log of amplification factor.
+        ue0 = cfgFlow.extPar[2]   # Inviscid flow velocity.
+        Ct0 = 0                   # Stagnation point is laminar. 
+
+        cfgFlow.blPar[cfgFlow.bNodes[0]*cfgFlow.nBlPar + 0] = ue0*T0*cfgFlow.Re       # Reynolds number based on the momentum thickness
+
+        if cfgFlow.blPar[cfgFlow.bNodes[0]*cfgFlow.nBlPar+0] < 1:
+            cfgFlow.blPar[cfgFlow.bNodes[0]*cfgFlow.nBlPar+0] = 1
+            ue0 = cfgFlow.blPar[cfgFlow.bNodes[0]*cfgFlow.nBlPar + 0] / (T0*cfgFlow.Re)
+
+        cfgFlow.blPar[cfgFlow.bNodes[0]*cfgFlow.nBlPar + 9] = T0 * H0                 # deltaStar
+
+        if cfgFlow.xtrF == [0, 0]:
+          # Turbulent closures.
+
+          Ct0 = 0.03
+          cfgFlow.blPar[cfgFlow.bNodes[0]*cfgFlow.nBlPar+1:cfgFlow.bNodes[0]*cfgFlow.nBlPar+9] = Closures.turbulentClosure(T0, H0,
+                                                                                                    Ct0, cfgFlow.blPar[cfgFlow.bNodes[0]*cfgFlow.nBlPar+0],
+                                                                                                    cfgFlow.extPar[cfgFlow.bNodes[0]*cfgFlow.nExtPar+0], 'wake')
+        else:
+          # Laminar closures.
+
+          cfgFlow.blPar[(cfgFlow.bNodes[0]*cfgFlow.nBlPar + 1):(cfgFlow.bNodes[0]*cfgFlow.nBlPar+7)]=Closures.laminarClosure(T0, H0,
+                                                                                                          cfgFlow.blPar[cfgFlow.bNodes[0]*cfgFlow.nBlPar+0],
+                                                                                                          cfgFlow.extPar[cfgFlow.bNodes[0]*cfgFlow.nExtPar+0],
+                                                                                                          'airfoil')
+        
+        # Set boundary condition.
+
+        cfgFlow.u[cfgFlow.bNodes[0]*cfgFlow.nVar:(cfgFlow.bNodes[0]+1)*cfgFlow.nVar] = [T0, H0, n0, ue0, Ct0]
+        pass
+
+    # Wake boundary conditions 
+    def wakeBC(self, cfgFlow):
+        """Boundary conditions at the first point of the wake
+        based on the parameters at the last point on the upper and lower sides.
+        Drela (1989).
+        """
+
+        # Extract last points on the upper and lower sides
+
+        xEnd_up = cfgFlow.u[(cfgFlow.bNodes[1]-1)*cfgFlow.nVar:cfgFlow.bNodes[1]*cfgFlow.nVar]
+        xEnd_lw = cfgFlow.u[(cfgFlow.bNodes[2]-1)*cfgFlow.nVar:cfgFlow.bNodes[2]*cfgFlow.nVar]
+
+        T0  = xEnd_up[0]+xEnd_lw[0]                                                # (Theta_up+Theta_low).
+        H0  = (xEnd_up[0]*xEnd_up[1]+xEnd_lw[0]*xEnd_lw[1])/T0                     # ((Theta*H)_up+(Theta*H)_low)/(Theta_up+Theta_low).
+        n0  = 9
+        ue0 = cfgFlow.extPar[cfgFlow.bNodes[2]*cfgFlow.nExtPar+2]                  # Inviscid velocity.
+        Ct0 = (xEnd_up[0]*xEnd_up[4]+xEnd_lw[0]*xEnd_lw[4])/T0                     # ((Ct*Theta)_up+(Theta)_low)/(Ct*Theta_up+Theta_low).
+
+        cfgFlow.blPar[cfgFlow.bNodes[2]*cfgFlow.nBlPar+0] = ue0*T0*cfgFlow.Re      # Reynolds number based on the momentum thickness.
+        
+        if cfgFlow.blPar[cfgFlow.bNodes[2]*cfgFlow.nBlPar+0] < 1:
+            cfgFlow.blPar[cfgFlow.bNodes[2]*cfgFlow.nBlPar+0] = 1
+            ue0 = cfgFlow.blPar[cfgFlow.bNodes[2]*cfgFlow.nBlPar+0] / (T0*cfgFlow.Re)
+
+        cfgFlow.blPar[cfgFlow.bNodes[2]*cfgFlow.nBlPar+9] = T0*H0                  # deltaStar
+        
+        # Turbulent closures.
+
+        cfgFlow.blPar[cfgFlow.bNodes[2]*cfgFlow.nBlPar+1:cfgFlow.bNodes[2]*cfgFlow.nBlPar+9] = Closures.turbulentClosure(T0, H0,
+                                                                                                    Ct0, cfgFlow.blPar[cfgFlow.bNodes[2]*cfgFlow.nBlPar+0],
+                                                                                                    cfgFlow.extPar[cfgFlow.bNodes[2]*cfgFlow.nExtPar+0], 'wake')
+
+        # Set boundary condition. 
+
+        cfgFlow.u[cfgFlow.bNodes[2]*cfgFlow.nVar:(cfgFlow.bNodes[2]+1)*cfgFlow.nVar] = [T0, H0, n0, ue0, Ct0]
+        pass
\ No newline at end of file

dart/viscous/Physics/DClosures.py

+################################################################################
+#                                                                              #
+#                            DARTFlo viscous module file                       #
+#                       Closures: Closure models to close the                  #
+#                   system of PDEs of the boundary layer equations             #
+#                                                                              #
+# Author: Paul Dechamps                                                        #
+# Université de Liège                                                          #
+# Date: 2021.09.10                                                             #
+#                                                                              #
+################################################################################
+
+
+# ------------------------------------------------------------------------------
+#  Imports
+# ------------------------------------------------------------------------------
+import numpy as np
+
+class Closures:
+    def __init__(self) -> None:
+        pass
+
+    def laminarClosure(theta, H, Ret,Me, name):
+        """Laminar closure derived from the Falkner-Skan one-parameter profile family
+            Nishida and Drela (1996)
+        """
+        Hk = (H - 0.29*Me**2)/(1+0.113*Me**2) # Kinematic shape factor
+        if name == "airfoil":
+            Hk = max(Hk, 1.02)
+        else:
+            Hk = max(Hk, 1.00005)
+        Hk = min(Hk,6)
+        Hstar2 = (0.064/(Hk-0.8)+0.251)*Me**2 # Density shape parameter
+        delta = min((3.15+ H + (1.72/(Hk-1)))*theta, 12*theta)
+        Hks = Hk - 4.35
+        if Hk <= 4.35:
+            dens = Hk + 1
+            Hstar = 0.0111*Hks**2/dens - 0.0278*Hks**3/dens + 1.528 -0.0002*(Hks*Hk)**2
+        elif Hk > 4.35:
+            Hstar = 1.528 + 0.015*Hks**2/Hk
+        Hstar = (Hstar + 0.028*Me**2)/(1+0.014*Me**2) # Whitfield's minor additional compressibility correction
+        if Hk < 5.5:
+            tmp = (5.5-Hk)**3 / (Hk+1)
+            Cf = (-0.07 + 0.0727*tmp)/Ret
+        elif Hk >= 5.5:
+            tmp = 1 - 1/(Hk-4.5)
+            Cf = (-0.07 + 0.015*tmp**2) /Ret
+        if Hk < 4:
+            Cd = ( 0.00205*(4.0-Hk)**5.5 + 0.207 ) * (Hstar/(2*Ret))
+        elif Hk >= 4:
+            HkCd = (Hk-4)**2
+            denCd = 1+0.02*HkCd
+            Cd = ( -0.0016*HkCd/denCd + 0.207) * (Hstar/(2*Ret))
+        if name == 'wake':
+            Cd = 2*(1.10*(1.0-1.0/Hk)**2/Hk)* (Hstar/(2*Ret))
+            Cf = 0
+        return Cd, Cf, delta, Hstar, Hstar2, Hk
+
+
+    def turbulentClosure(theta, H, Ct, Ret, Me, name):
+        """Turbulent closure derived from Nishida and Drela (1996)
+        """
+        
+        H = max(H, 1.0005)
+        # eliminate absurd transients
+        Ct = max(min(Ct, 0.30), 0.0000001)
+        Hk = (H - 0.29*Me**2)/(1+0.113*Me**2)
+        Hk = max(min(Hk, 10), 1.00005)
+        Hstar2 = ((0.064/(Hk-0.8))+0.251)*Me**2
+        #gam = 1.4 - 1
+        #Fc = np.sqrt(1+0.5*gam*Me**2)
+        Fc = np.sqrt(1+0.2*Me**2)
+
+        if Ret <= 400:
+            H0 = 4
+        elif Ret > 400:
+            H0 = 3 + (400/Ret)
+
+        Ret = max(Ret, 200)
+ 
+        if Hk <= H0:
+            Hstar = ((0.5-4/Ret)*((H0-Hk)/(H0-1))**2)*(1.5/(Hk+0.5))+1.5+4/Ret
+        elif Hk > H0:
+            Hstar = (Hk-H0)**2*(0.007*np.log(Ret)/(Hk-H0+4/np.log(Ret))**2 + 0.015/Hk) + 1.5 + 4/Ret
+
+        Hstar = (Hstar + 0.028*Me**2)/(1+0.014*Me**2) # Whitfield's minor additional compressibility correction
+        #logRt = np.log(Ret/Fc)
+        #logRt = max(logRt,3)
+        #arg = max(-1.33*Hk, -20)
+        Us = (Hstar/2)*(1-4*(Hk-1)/(3*H)) # Equivalent normalized wall slip velocity
+        delta = min((3.15+ H + (1.72/(Hk-1)))*theta, 12*theta)
+        Ctcon = 0.5/(6.7**2*0.75)
+
+        if name == 'wake':
+            Us = min(Us, 0.99995)
+            n = 2 # two wake halves
+            Cf = 0 # no friction inside the wake
+            Hkc = Hk - 1
+            Cdw = 0  # wall contribution
+            Cdd = (0.995-Us)*Ct**2 # turbulent outer layer contribution
+            Cdl = 0.15*(0.995-Us)**2/Ret # laminar stress contribution to outer layer
+            Cteq = np.sqrt(4*Hstar*Ctcon*((Hk-1)*Hkc**2)/((1-Us)*(Hk**2)*H)) #Here it is Cteq^0.5
+        else:
+            if Us > 0.95:
+                Us = 0.98
+            n = 1
+            #Cf = (1/Fc)*(0.3*np.exp(arg)*(logRt/2.3026)**(-1.74-0.31*Hk)+0.00011*(np.tanh(4-(Hk/0.875))-1))
+            Cf = 1/Fc*(0.3*np.exp(-1.33*Hk)*(np.log10(Ret/Fc))**(-1.74-0.31*Hk) + 0.00011*(np.tanh(4-Hk/0.875) - 1))
+            Hkc = max(Hk-1-18/Ret, 0.01)
+            # Dissipation coefficient
+            Hmin = 1 + 2.1/np.log(Ret)
+            Fl = (Hk-1)/(Hmin-1)
+            Dfac = 0.5+0.5*np.tanh(Fl)
+            Cdw = 0.5*(Cf*Us) * Dfac
+            Cdd = (0.995-Us)*Ct**2
+            Cdl = 0.15*(0.995-Us)**2/Ret
+            Cteq = np.sqrt(Hstar*Ctcon*((Hk-1)*Hkc**2)/((1-Us)*(Hk**2)*H)) #Here it is Cteq^0.5
+        Cd = n*(Cdw+ Cdd + Cdl)
+        #Ueq = 1.25/Hk*(Cf/2-((Hk-1)/(6.432*Hk))**2)
+        return Cd, Cf, delta, Hstar, Hstar2, Hk, Cteq, Us
+
+    def amplRoutine(Hk, Ret, theta):
+        '''Compute the amplitude of the TS waves envelop for the transition
+        '''
+        Dgr = 0.08
+        Hk1 = 3.5
+        Hk2 = 4
+        Hmi = (1/(Hk-1))
+        logRet_crit = 2.492*(1/(Hk-1))**0.43 + 0.7*(np.tanh(14*Hmi-9.24)+1.0) # value at which the amplification starts to grow
+        if Ret <=0:
+            Ret = 1
+        if np.log10(Ret) < logRet_crit - Dgr  :
+            Ax = 0
+        else:
+            Rnorm = (np.log10(Ret) - (logRet_crit-Dgr))/(2*Dgr)
+            if Rnorm <=0:
+                Rfac = 0
+            elif Rnorm >= 1:
+                Rfac = 1
+            else:
+                Rfac = 3*Rnorm**2 - 2*Rnorm**3
+            # normal envelope amplification rate
+            m = -0.05+2.7*Hmi-5.5*Hmi**2+3*Hmi**3+0.1*np.exp(-20*Hmi)
+            arg = 3.87*Hmi-2.52
+            dndRet = 0.028*(Hk-1)-0.0345*np.exp(-(arg**2))
+            Ax = (m*dndRet/theta)*Rfac
+            # Set correction for separated profiles
+            if Hk > Hk1:
+                Hnorm = (Hk - Hk1)/(Hk2-Hk1)
+                if Hnorm >=1:
+                    Hfac = 1
+                else:
+                    Hfac = 3*Hnorm**2 - 2*Hnorm**3
+                Ax1 = Ax
+                Gro = 0.3+0.35*np.exp(-0.15*(Hk-5))
+                Tnr = np.tanh(1.2*(np.log10(Ret)-Gro))
+                Ax2 = (0.086*Tnr - 0.25/(Hk-1)**1.5)/theta
+                if Ax2 < 0:
+                    Ax2 = 0
+                Ax = Hfac*Ax2 + (1-Hfac)*Ax1
+                Ax = max(Ax, 0)
+        return Ax
\ No newline at end of file

dart/viscous/Physics/DDataStructures.py

+################################################################################
+#                                                                              #
+#                            DARTFlo viscous module file                       #
+#                          Variables: Solution and external                    #
+#                             flow parameters container                        #
+#                                                                              #
+# Author: Paul Dechamps                                                        #
+# Université de Liège                                                          #
+# Date: 2021.09.10                                                             #
+#                                                                              #
+################################################################################
+
+
+# ------------------------------------------------------------------------------
+#  Imports
+# ------------------------------------------------------------------------------
+import numpy as np
+
+class PGeometry:
+
+    def __init__(self, data, paramDict, _bNodes):
+
+        # Private
+        
+        self.__globCoord     = data[:,0]                        # Global coordinates of each marker
+        self.__locCoord      = paramDict['xx']                  # Local coordinates of each marker
+        self.__nMarkers      = len(self.__locCoord)             # Number of markers in the computational domain
+        
+        self.__boundaryNodes = _bNodes                          # Boundary nodes of the computational domain [Stagnation point,
+                                                                #                                             First point on lower side,
+                                                                #                                             First wake point]
+
+    def GetglobCoord(self): return self.__globCoord
+    def GetlocCoord(self): return self.__locCoord
+    def GetnMarkers(self): return self.__nMarkers
+    def GetglobCoord(self): return self.__bNodes
+    def GetboundaryNodes(self): return self.__boundaryNodes
+
+
+class PSolution:
+
+    def __init__(self, _xtrF):
+
+        # Private
+
+        self.__nVar = 5                                           # Number of variables involved in the system of PDEs
+        self.__nBlPar = 10                                        # Number of boundary layer parameters
+        
+        self.__xtrF = _xtrF                                       # Forced Transition location on [upper, lower] sides
+        self.__xtr = [1, 1]                                       # Transition location
+
+        # Public
+
+        # Solution vectors
+
+        self.u = []                                               # Vector containing the solution [theta, H, N, ue, Ctau]
+
+        # Boundary layer parameters
+
+        self.ReTheta =                                            # Reynolds number based on the variable 'theta'.
+        self.disspCoeff =                                         # Dissipation coefficient.
+        self.frictionCoeff =                                      # Wall friction coefficient.
+        self.delta =                                              # Boundary layer thickness.
+        self.Hstar =                                              # Kinetic energy shape parameter.
+        self.Hstar2 =                                             # Density shape parameter.
+        self.Hk =                                                 # Kinematic shape factor.
+        self.Cteq =                                               # Equilibrium shear stress coefficient.
+        self.Us =                                                 # Equivalent normalized wall slip velocity.
+        self.deltaStar =                                          # Displacement thickness of the boundary layer. 
+
+    def GetnVar(self): return self.__nVar
+    def GetnBlPar(self): return self.__nBlPar
+
+    def GetxtrF(self): return self.__xtrF
+    def Getxtr(self): return self.__xtr
+    def Setxtr(self, xtrIn, side):
+        self.__xtr[side] = xtrIn
+        pass
+
+    def initializeSolution(self, PGeometry):
+
+        nMarkers = PGeometry.GetnMarkers()
+        self.u = np.empty(self.__nVar * nMarkers)
+        self.blPar = np.empty(self.__nVar * nMarkers)
+
+        return self.u, self.blPar
+
+class POuterFlow:
+
+    def __init__(self, data) -> None:
+
+        self.__Mach =                                               # Outer Mach number.
+        self.__density =                                            # Outer flow density.
+        self.__velocity =                                           # Outer flow velocity.
+        self.__dv =                                                 # Outer flow velocity gradient. 
+
+        self.outerDeltaStar = 
+        self.outerLocCoord = 
+        self.outerRe = 
+
+class PDataContainer:
+    def __init__(self, _Re) -> None:
+
+        self.Re         = _Re               # Flow reynolds number.
+        self.Geometry  = PGeometry()       # Geometry components.
+        self.Solution  = PSolution()       # Solution components.
+        self.OuterFlow = POuterFlow()      # Outer (inviscid) flow parameters.
+
+
+        
+
+
+class DVariables:
+    """Data container for viscous solver.
+
+    Attributes
+    ----------
+    
+    - Ti: float [0,1]
+        - User defined turbulence intensity that impacts the critical amplification ratio Ncrit of the e^N method.
+    - Re: float
+        - Flow Reynolds number.
+
+    - xx: numpy array
+        - Airfoil surfaces and wake mesh.
+
+    - nN: int
+        - Number of cell in the computational domain.
+
+    - bNodes: numpy array
+            - 3 elements array containing the boundary nodes of the computational domain [Stagnation point, First lower side point, First wake point].
+    
+    - nVar: int
+          - Number of variables involved in the differential problem.
+
+    - u: numpy array (vector of len(nVar*nN))
+         - Array containing the solution (nVar variables at each of the nN points).
+
+    - blPar: numpy array (vector of len(nBlPar*nN))
+           - Boundary layer parameters.
+
+    - extPar: numpy array (vector of len(nExtPar*nN))
+            - External (inviscid) flow parameters.
+
+    - dv: numpy array (vector of len(nN))
+        - Velocity gradient on each cell.
+
+    - xGlobal: numpy array (vector of len(nN))
+             - Cell position non-dimensionalized by the airfoil chord: x/c
+
+    - xtrF: numpy array of len(2)
+          - User forced transiton on [Upper side, Lower side]
+    
+    - xtr: numpy array of len(2)
+         - Flow laminar to turbulent transition location on [Upper side, Lower side]
+        
+    - Ncrit: float
+            - Critical amplification ratio for the e^N method (default value= 9.) 
+
+    """
+    def __init__(self, data, paramDict, bNodes, xtrF, Ti, Re):
+        self.__Ti=Ti
+        self.Re=Re
+        self.xx=paramDict['xx']
+        self.nN=len(self.xx)
+
+        self.bNodes=bNodes
+        # Solution 
+        self.nVar=5
+        self.u=np.zeros(self.nVar*self.nN)
+
+        # Parameters of the boundary layers
+        # [ 0 ,  1,  2,     3,     4,      5,  6,    7,  8,         9]
+        # [Ret, Cd, Cf, delta, Hstar, Hstar2, Hk, Cteq, Us, deltaStar]
+        self.nBlPar=10
+        self.blPar=np.zeros(self.nBlPar*self.nN)
+        self.Ret_previous=self.blPar[0::self.nVar]
+
+        self.dueOld_dx=np.zeros(self.nN)
+        self.ddeltaStarOld_dx=np.zeros(self.nN)
+
+        # Parameters of the external (inviscid) flow
+        # [ 0,    1,     2,            3,     4,    5]
+        # [Me, rhoe, vtExt, deltaStarExt, xxExt ReExt]
+        self.nExtPar=6
+        self.extPar=np.zeros(self.nExtPar*self.nN)
+        self.extPar[0::self.nExtPar]=data[:,5]
+        self.extPar[1::self.nExtPar]=data[:,6]
+        self.extPar[2::self.nExtPar]=paramDict['vtExt']
+        self.extPar[3::self.nExtPar]=data[:,7]
+        self.extPar[4::self.nExtPar]=paramDict['xx']
+        self.dv=paramDict['dv']
+        self.xGlobal=data[:,0]
+
+        self.turb=np.zeros(self.nN, dtype=int)
+
+        self.xtrF=xtrF
+        self.xtr=np.zeros(2)
+        if not len(self.xtrF)==2:
+            raise RuntimeError(ccolors.ANSI_RED + "Error in forced transition." + ccolors.ANSI_RESET)
+        for n in range(len(self.xtrF)):
+            if self.xtrF[n] is not None:
+                if 0<=self.xtrF[n]<=1:
+                    self.xtr[n]=self.xtrF[n]
+                else:
+                    raise ValueError(ccolors.ANSI_RED + "Non-physical transition location." + ccolors.ANSI_RESET)
+            else:
+                self.xtr[n]=1
+
+        self.idxTr=[self.bNodes[1]-1,self.bNodes[2]-1]
+
+        self.Ncrit=self.turbIntensity(Ti)
+        pass
+    
+    def turbIntensity(self,Ti):
+        '''Computes the critical amplification ratio Ncrit based on the input
+        free-stream turbulence intensity. Ncrit=9 by default.'''
+        if Ti is not None and Ti<=1:
+            tau=2.7*np.tanh(self.Ti/2.7)
+            Ncrit=-8.43-2.4*np.log(tau/100) # Drela 1998
+        else:
+            Ncrit=9
+        return Ncrit
\ No newline at end of file

dart/viscous/Physics/DDiscretization.py

+################################################################################
+#                                                                              #
+#                             Flow viscous module file                         #
+#               Discretization: Computes solution gradient at one point        #
+#                       of the computational domain. Different                 #
+#                   finite difference spacial schemes are implemented.         #
+#                                                                              #
+# Author: Paul Dechamps                                                        #
+# Université de Liège                                                          #
+# Date: 2021.09.10                                                             #
+#                                                                              #
+################################################################################
+
+
+# ------------------------------------------------------------------------------
+#  Imports
+# ------------------------------------------------------------------------------
+import numpy as np
+
+class Discretization:
+    """Discretize the spacial derivatives required to evaluate the boundary layer laws.
+    
+    Attributes
+    ----------
+    
+    - xx: Cell coordinates in the (local) x-direction in the boundary layer.
+
+    - xxExt: Cell coordinates in the (local) x-direction for the flow outside the boundary layer.
+
+    - bNodes: Boundary nodes of the computational domain
+    
+    Returns
+    -------
+    
+    Each method of the class returns the derivatives of the solution and of boundary layer parameters
+    
+    - du_dx: Derivative of the solution [dTheta/dx, dH/dx, dN/dx, due/dx, dCt/dx]
+
+    - dHstar_dx: Derivative of Hstar
+
+    - dueExt_dx, ddeltaStarExt_dx: External flow derivatives. Used in the interaction law.
+    
+    - c, cExt: Parameters dependent on the cell size 'dx', 'dxExt'. Used in the interation law.
+    """
+    def __init__(self, xx, xxExt, bNodes):
+        dx=np.zeros(len(xx))
+        for i in range(1,len(xx)):
+            dx[i]=xx[i]-xx[0] if i==bNodes[1] else xx[i]-xx[i-1]
+        self.dx=dx
+
+        dxExt=np.zeros(len(xxExt))
+        for i in range(1,len(xx)):
+            dxExt[i]=xxExt[i]-xxExt[0] if i==bNodes[1] else xxExt[i]-xxExt[i-1]
+        self.dxExt=dxExt
+
+    def fou(self, var, x, idxPrev, idx):
+        """First order backward difference.
+        """
+        du=x-var.u[var.nVar*idxPrev:var.nVar*(idxPrev+1)]
+
+        dx = var.xx[idx] - var.xx[idxPrev]
+        dxExt = var.extPar[idx*var.nExtPar + 4] - var.extPar[idxPrev*var.nExtPar + 4]
+        
+        du_dx=du/dx
+        dHstar_dx=(var.blPar[idx*var.nBlPar+4]-var.blPar[idxPrev*var.nBlPar+4])/dx
+        dueExt_dx=(var.extPar[idx*var.nExtPar+2]-var.extPar[idxPrev*var.nExtPar+2])/dxExt
+        ddeltaStarExt_dx=(var.extPar[idx*var.nExtPar+3]-var.extPar[idxPrev*var.nExtPar+3])/dxExt
+        c=4/(np.pi*dx)
+        cExt=4/(np.pi*dxExt)
+        return du_dx, dHstar_dx, dueExt_dx, ddeltaStarExt_dx, c, cExt
+    
+    def center(self, var, x, idxPrev2, idxPrev,idx, idxNext):
+        """Second order central difference.
+        """
+        du= (var.u[idxNext*var.nVar:(idxNext+1)*var.nVar] - var.u[idxPrev*var.nVar:(idxPrev+1)*var.nVar])
+        
+        dx = var.xx[idx] - var.xx[idxPrev]
+        dxExt = var.extPar[idx*var.nExtPar + 4] - var.extPar[idxPrev*var.nExtPar + 4]
+
+        du_dx = du / (2*dx)
+        dHstar_dx = (var.blPar[idxNext*var.nBlPar+4]-var.blPar[idxPrev*var.nBlPar+4]) / (2*dx)
+        dueExt_dx = (var.extPar[idxNext*var.nExtPar+2] - var.extPar[idxPrev*var.nExtPar+2]) / (2*dxExt)
+        ddeltaStarExt_dx=(var.extPar[idxNext*var.nExtPar+3]-var.extPar[idxPrev*var.nExtPar+3])/(2*dxExt)
+        c=4/(np.pi*dx)
+        cExt=4/(np.pi*dxExt)
+        return du_dx, dHstar_dx, dueExt_dx, ddeltaStarExt_dx, c, cExt
+
+    # Second order backward difference
+    def sou(self, var, x,idxPrev2, idxPrev,idx, idxNext): # Attention if we want gradients at 2nd point after stagnation point on lower side. 
+        """Second order upwing difference.
+        """
+        du=-3*x+4*var.u[+var.nVar*idxPrev:+var.nVar*(idxPrev+1)]-var.u[+var.nVar*idxPrev2:+var.nVar*(idxPrev2+1)]
+        du_dx=du/(2*self.dx[idx])
+        dHstar_dx=(3*var.blPar[idx*var.nBlPar+4]-4*var.blPar[idxPrev*var.nBlPar+4]+var.blPar[idxPrev2*var.nBlPar])/(2*self.dx[idx])
+        dueExt_dx=(var.extPar[idxNext*var.nExtPar+2]-var.extPar[idxPrev*var.nExtPar+2])/(2*self.dx[idx])
+        ddeltaStarExt_dx=(var.extPar[idxNext*var.nExtPar+3]-var.extPar[idxPrev*var.nExtPar+3])/(2*self.dx[idx])
+        c=4/(np.pi*self.dx[idx])
+        cExt=4/(np.pi*self.dxExt[idx])
+        return du_dx, dHstar_dx, dueExt_dx, ddeltaStarExt_dx, c, cExt
\ No newline at end of file

dart/viscous/Physics/DValidation.py

+################################################################################
+#                                                                              #
+#                            DARTFlo viscous module file                       #
+#                      Validation: Class that displays the results             #
+#                          and compare them with given results                 #
+#                                                                              #
+# Author: Paul Dechamps                                                        #
+# Université de Liège                                                          #
+# Date: 2021.09.10                                                             #
+#                                                                              #
+################################################################################
+
+
+# ------------------------------------------------------------------------------
+#  Imports
+# ------------------------------------------------------------------------------
+import numpy as np
+import matplotlib.pyplot as plt
+
+class Validation:
+    def __init__(self,_xfoilFileName = None):
+        self.xfoilFileName=_xfoilFileName
+    def main(self, dartcl, wData, wakeblVectors = None, xWake = None):
+        # Load xfoil data
+
+        xFoilIn = 0
+        cl_xfoil = 0
+        cd_xfoil = 0
+
+        if self.xfoilFileName is not None:
+
+            xFoilIn = 1
+            with open('/Users/pauldechamps/lab/dartflo/dart/viscous/Physics/XfoilFiles/' + self.xfoilFileName, 'r') as f:
+                cl_xfoil=f.readline()
+                cd_xfoil=f.readline()
+                f.close()
+            # columns : x,ue,theta,Cf,H
+            data=np.loadtxt('/Users/pauldechamps/lab/dartflo/dart/viscous/Physics/XfoilFiles/'+self.xfoilFileName,skiprows=3,usecols=(1,3,4,5,6,7))
+
+            # columns : x, delta*, theta, H, Cf, ue
+            data[:,[1,2,3,4,5]]=data[:,[2,3,5,4,2]]
+            data[:,]
+
+            # Find upper/lower side separation 
+            for i in range(len(data[:,0])):
+                if data[i,0]<=data[i+1,0]:
+                    xStag=i
+                    break
+        for i in range(len(wData['x'])):
+            if wData['x'][i]==1:
+                xTeDart=i+1
+                break
+
+        # Aerodynamic coefficients
+        print('+----------------------------  Compare with XFOIL -----------------------------+')
+        print('- CL; Dart: {:<.5f}, Xfoil: {:<.5f}.'.format(dartcl, float(cl_xfoil)))
+        print('- CD; Dart: {:<.5f}, Xfoil: {:<.5f}.'.format(wData['Cd'], float(cd_xfoil)))
+
+        # Plot variables
+        plot_dStar=plt.figure(1)
+        plt.plot(wData['x'][:xTeDart],wData['delta*'][:xTeDart],'b-',label=r'DARTflo upper')
+        plt.plot(wData['x'][xTeDart:],wData['delta*'][xTeDart:],'r-',label=r'DARTflo lower')
+        plt.xlim([0,1])
+        if wakeblVectors is not None:
+            plt.plot(xWake, wakeblVectors[0],label=r'DARTflo wake')
+            plt.xlim([0,2])
+        
+        if xFoilIn:
+            plt.plot(data[:xStag,0],data[:xStag,1],'b--',label=r'Xfoil upper')
+            plt.plot(data[xStag:,0],data[xStag:,1],'r--',label=r'Xfoil lower')
+            #plt.plot(data[xWake:,0],data[xWake:,1],'r--',label=r'Xfoil Wake')
+        plt.legend()
+        plt.xlabel('$x/c$')
+        plt.ylabel('$\delta^*$')
+
+        # H
+        plot_H=plt.figure(3)
+        plt.plot(wData['x'][:xTeDart],wData['H'][:xTeDart],'b-',label=r'DARTflo upper')
+        plt.plot(wData['x'][xTeDart:],wData['H'][xTeDart:],'r-',label=r'DARTflo lower')
+        plt.xlim([0,1])
+        if wakeblVectors is not None:
+            plt.plot(xWake, wakeblVectors[2],label=r'DARTflo wake')
+            plt.xlim([0,2])
+
+        if xFoilIn:
+            plt.plot(data[:xStag,0],data[:xStag,3],'b--',label=r'Xfoil upper')
+            plt.plot(data[xStag:,0],data[xStag:,3],'r--',label=r'Xfoil lower')
+            #plt.plot(data[xWake:,0],data[xWake:,3],'r--',label=r'Xfoil Wake')
+        plt.legend()
+        plt.xlabel('$x/c$')
+        plt.ylabel('$H$')
+        # Cf
+        plot_Cf=plt.figure(4)
+        plt.plot(wData['x'][:xTeDart], wData['Cf'][:xTeDart],'b-',label=r'DARTflo upper')
+        plt.plot(wData['x'][xTeDart:], wData['Cf'][xTeDart:],'r-',label=r'DARTflo lower')
+        plt.xlim([0,1])
+        if wakeblVectors is not None:
+            plt.plot(xWake, wakeblVectors[6],label=r'DARTflo wake')
+            plt.xlim([0,2])
+
+        if xFoilIn:
+            plt.plot(data[:xStag,0],data[:xStag,4],'b--',label=r'Xfoil upper')
+            plt.plot(data[xStag:,0],data[xStag:,4],'r--',label=r'Xfoil lower')
+            #plt.plot(data[xWake:,0],data[xWake:,4],'r--',label=r'Xfoil Wake')
+        plt.legend()
+        plt.xlabel('$x/c$')
+        plt.ylabel('$C_f$')
+
+        plt.show()
+

dart/viscous/Physics/DVariables.py

+################################################################################
+#                                                                              #
+#                            DARTFlo viscous module file                       #
+#                          Variables: Solution and external                    #
+#                             flow parameters container                        #
+#                                                                              #
+# Author: Paul Dechamps                                                        #
+# Université de Liège                                                          #
+# Date: 2021.09.10                                                             #
+#                                                                              #
+################################################################################
+
+
+# ------------------------------------------------------------------------------
+#  Imports
+# ------------------------------------------------------------------------------
+import numpy as np 
+from matplotlib import pyplot as plt 
+from fwk.coloring import ccolors
+
+class Variables:
+    """Data container for viscous solver.
+
+    Attributes
+    ----------
+    
+    - Ti: float [0,1]
+        - User defined turbulence intensity that impacts the critical amplification ratio Ncrit of the e^N method.
+    - Re: float
+        - Flow Reynolds number.
+
+    - xx: numpy array
+        - Airfoil surfaces and wake mesh.
+
+    - nN: int
+        - Number of cell in the computational domain.
+
+    - bNodes: numpy array
+            - 3 elements array containing the boundary nodes of the computational domain [Stagnation point, First lower side point, First wake point].
+    
+    - nVar: int
+          - Number of variables involved in the differential problem.
+
+    - u: numpy array (vector of len(nVar*nN))
+         - Array containing the solution (nVar variables at each of the nN points).
+
+    - blPar: numpy array (vector of len(nBlPar*nN))
+           - Boundary layer parameters.
+
+    - extPar: numpy array (vector of len(nExtPar*nN))
+            - External (inviscid) flow parameters.
+
+    - dv: numpy array (vector of len(nN))
+        - Velocity gradient on each cell.
+
+    - xGlobal: numpy array (vector of len(nN))
+             - Cell position non-dimensionalized by the airfoil chord: x/c
+
+    - xtrF: numpy array of len(2)
+          - User forced transiton on [Upper side, Lower side]
+    
+    - xtr: numpy array of len(2)
+         - Flow laminar to turbulent transition location on [Upper side, Lower side]
+        
+    - Ncrit: float
+            - Critical amplification ratio for the e^N method (default value= 9.) 
+
+    """
+    def __init__(self, fMeshDict, xtrF, Ti, Re):
+        self.__Ti = Ti
+        self.Re   = Re
+        self.xx   = fMeshDict['locCoord']
+        self.nN   = len(self.xx)
+
+        self.bNodes = fMeshDict['bNodes']
+        # Solution 
+        self.nVar = 5
+        self.u = np.zeros(self.nVar*self.nN)
+
+        # Parameters of the boundary layers
+        # [ 0 ,  1,  2,     3,     4,      5,  6,    7,  8,         9]
+        # [Ret, Cd, Cf, delta, Hstar, Hstar2, Hk, Cteq, Us, deltaStar]
+        self.nBlPar = 10
+        self.blPar        = np.zeros(self.nBlPar*self.nN)
+        self.Ret_previous = self.blPar[0::self.nVar]
+
+        self.dueOld_dx = np.zeros(self.nN)
+        self.ddeltaStarOld_dx = np.zeros(self.nN)
+
+        # Parameters of the external (inviscid) flow
+        # [ 0,    1,     2,            3,     4,    5]
+        # [Me, rhoe, vtExt, deltaStarExt, xxExt ReExt]
+        self.nExtPar = 6
+
+        self.extPar                  = np.zeros(self.nExtPar*self.nN)
+        self.extPar[0::self.nExtPar] = fMeshDict['Me']
+        self.extPar[1::self.nExtPar] = fMeshDict['rho']
+        self.extPar[2::self.nExtPar] = fMeshDict['vtExt']
+        self.extPar[3::self.nExtPar] = fMeshDict['deltaStarExt']
+        self.extPar[4::self.nExtPar] = fMeshDict['locCoord']
+        self.dv                      = fMeshDict['dv']
+        self.xGlobal                 = fMeshDict['globCoord']
+
+        self.flowRegime = np.zeros(self.nN, dtype=int)
+
+        self.xtrF=xtrF
+        self.xtr=np.zeros(2)
+        if not len(self.xtrF)==2:
+            raise RuntimeError(ccolors.ANSI_RED + "Error in forced transition." + ccolors.ANSI_RESET)
+        for n in range(len(self.xtrF)):
+            if self.xtrF[n] is not None:
+                if 0<=self.xtrF[n]<=1:
+                    self.xtr[n]=self.xtrF[n]
+                else:
+                    raise ValueError(ccolors.ANSI_RED + "Non-physical transition location." + ccolors.ANSI_RESET)
+            else:
+                self.xtr[n]=1
+
+        self.transMarkers=[self.bNodes[1]-1,self.bNodes[2]-1]
+
+        self.Ncrit=self.turbIntensity(Ti)
+        pass
+    
+    def turbIntensity(self,Ti):
+        """Computes the critical amplification ratio Ncrit based on the input
+        free-stream turbulence intensity. Ncrit=9 by default.
+        """
+
+        if Ti is not None and Ti <= 1:
+            tau = 2.7 * np.tanh(self.Ti /  2.7)
+            Ncrit = -8.43 - 2.4*np.log(tau / 100) # Drela 1998
+        else:
+            Ncrit = 9
+        return Ncrit