diff --git a/dart/tests/bliHighLift.py b/dart/tests/bliHighLift.py
new file mode 100755
index 0000000000000000000000000000000000000000..f7a15ed49b9aaed683a9ff9d67c48e6691d0e3b5
--- /dev/null
+++ b/dart/tests/bliHighLift.py
@@ -0,0 +1,174 @@
+#!/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.VCoupler as floVC
+import dart.viscous.Drivers.VDriver as floVSMaster
+
+import dart.viscous.Solvers.Steady.StdCoupler as floVC_Std
+import dart.viscous.Solvers.Steady.StdSolver as floVS_Std
+
+import dart.viscous.Physics.VValidation as validation
+
+import 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 = 6e6
+ alpha = 12.*math.pi/180
+ M_inf = 0
+ #user_xtr=[0.0106,1]
+ user_xtr=[None,None]
+ user_Ti=None
+
+ mapMesh = 0
+
+ # Time solver parameters
+ Time_Domain=True # True for unsteady solver, False for steady solver
+ Cfl = 1.5
+ spaceMarching=True
+
+ if Time_Domain is True:
+ timeParams={'CFL0': Cfl, 'spaceMarching':spaceMarching}
+
+ # ========================================================================================
+
+ U_inf = [math.cos(alpha), math.sin(alpha)] # norm should be 1
+ c_ref = 1
+ dim = 2
+ tol = 1e-1 # tolerance for the VII (usually one drag count)
+ case='Case12Xfoil.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.00001}
+ #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, timeParams, M_inf, case, mapMesh)
+ coupler = floVC.Coupler(isolver, vsolver, blw[0], blw[1], tol, gmshWriter, bnd, timeParams)
+
+ 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)
+
+ # 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)
+
+ # check results
+ print(ccolors.ANSI_BLUE + 'PyTesting...' + ccolors.ANSI_RESET)
+ tests = CTests()
+ if Re == 1e7 and M_inf == 0 and alpha == 5*math.pi/180:
+ tests.add(CTest('Cl', isolver.Cl, 0.55, 5e-2)) # Xfoil 0.58
+ tests.add(CTest('Cd', vsolver.Cd, 0.0062, 1e-3, forceabs=True))
+ tests.add(CTest('Cdp', vsolver.Cdp, 0.0018, 1e-3, forceabs=True))
+ tests.add(CTest('xtr_top', vsolver.xtr[0], 0.061, 0.03, forceabs=True)) # Xfoil 0.056
+ tests.add(CTest('xtr_bot', vsolver.xtr[1], 0.740, 0.03, forceabs=True))
+ elif Re == 1e7 and M_inf == 0.5 and alpha == 5*math.pi/180:
+ tests.add(CTest('Cl', isolver.Cl, 0.65, 5e-2)) # Xfoil 0.69
+ tests.add(CTest('Cd', vsolver.Cd, 0.0067, 1e-3, forceabs=True))
+ tests.add(CTest('Cdp', vsolver.Cdp, 0.0025, 1e-3, forceabs=True))
+ tests.add(CTest('xtr_top', vsolver.xtr[0], 0.049, 0.03, forceabs=True)) # Xfoil 0.038
+ tests.add(CTest('xtr_bot', vsolver.xtr[1], 0.736, 0.03, forceabs=True))
+ elif Re == 1e7 and M_inf == 0 and alpha == 12*math.pi/180:
+ tests.add(CTest('Cl', isolver.Cl, 1.30, 5e-2)) # Xfoil 1.39
+ tests.add(CTest('Cd', vsolver.Cd, 0.011, 1e-3, forceabs=True))
+ tests.add(CTest('Cdp', vsolver.Cdp, 0.0064, 1e-3, forceabs=True))
+ tests.add(CTest('xtr_top', vsolver.xtr[0], 0.010, 0.03, forceabs=True)) # Xfoil 0.008
+ tests.add(CTest('xtr_bot', vsolver.xtr[1], 1.000, 0.03, forceabs=True))
+ elif Re==1e7 and M_inf==0 and alpha==5*math.pi/180 and user_xtr==[0.01, 0.3]:
+ tests.add(CTest('Cl', isolver.Cl, 0.55, 5e-2)) # Xfoil 0.5687
+ tests.add(CTest('Cd', vsolver.Cd, 0.0077, 1e-3, forceabs=True)) # Xfoil 0.00777
+ tests.add(CTest('Cdp', vsolver.Cdp, 0.002, 1e-3, forceabs=True))
+ tests.add(CTest('xtr_top', vsolver.xtr[0], 0.01, 1e-3, forceabs=True))
+ tests.add(CTest('xtr_bot', vsolver.xtr[1], 0.3, 1e-3, forceabs=True))
+ else:
+ raise Exception('Test not defined for this flow')
+
+ tests.run()
+
+ # eof
+ print('')
+
+if __name__ == "__main__":
+ main()
diff --git a/dart/tests/bliLowRe.py b/dart/tests/bliLowRe.py
new file mode 100755
index 0000000000000000000000000000000000000000..8f0860cdc80104c2e09f7869c2f8b7143322cd30
--- /dev/null
+++ b/dart/tests/bliLowRe.py
@@ -0,0 +1,178 @@
+#!/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.VCoupler as floVC
+import dart.viscous.Drivers.VDriver as floVSMaster
+
+import dart.viscous.Solvers.Steady.StdCoupler as floVC_Std
+import dart.viscous.Solvers.Steady.StdSolver as floVS_Std
+
+import dart.viscous.Physics.VValidation as validation
+
+import 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()
+
+ # Time solver parameters
+ Time_Domain=True # True for unsteady solver, False for steady solver
+ model='implicitEuler' # 'eulerExp','RungeKutta', 'AdamsBashworth', 'MostAccurateExplicit', 'Leapfrog' or 'implicitEuler'
+ #user_xtr=[0.061, 0.740] # User forced transition : Use None type for free transition on either side
+ user_Ti=None
+
+
+ Re = 50e3
+ alpha=5*math.pi/180
+ user_xtr=[0.45,1]
+ M_inf=0.4
+ Cfl=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)
+
+
+
+ # Butcher's table coefficients
+ #RK_a=[[1],[-10890423/25193600, 5e3/7873],[-102261/5e6,-5121/2e4,7873/1e4]]
+ #RK_b=[1/10,1/6,0,1/6]
+ case='m0a5.txt' # File containing results from XFOIL of the considered case
+
+
+ try: user_xtr
+ except NameError: user_xtr = [None,None]
+
+ try: user_Ti
+ except NameError: user_Ti=None
+
+ try: RK_a
+ except NameError: RK_a=None
+ try: RK_b
+ except NameError: RK_b=None
+ timeParams={'modelName' : model, 'CFL': Cfl, 'Rka': RK_a, 'Rkb': RK_b}
+
+ # mesh the geometry
+ print(ccolors.ANSI_BLUE + 'PyMeshing...' + ccolors.ANSI_RESET)
+ tms['msh'].start()
+ pars = {'xLgt' : 5, 'yLgt' : 5, 'msF' : 1, 'msTe' : 0.01, 'msLe' : 0.01}
+ 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:
+ vsolver=floVSMaster.Driver(Re, user_xtr, timeParams, case)
+ coupler = floVC.Coupler(isolver, vsolver, blw[0], blw[1], tol, gmshWriter, bnd, timeParams)
+ elif Time_Domain is False:
+ 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)
+
+ # 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)
+
+ # check results
+ print(ccolors.ANSI_BLUE + 'PyTesting...' + ccolors.ANSI_RESET)
+ tests = CTests()
+ if Re == 1e7 and M_inf == 0 and alpha == 5*math.pi/180:
+ tests.add(CTest('Cl', isolver.Cl, 0.55, 5e-2)) # Xfoil 0.58
+ tests.add(CTest('Cd', vsolver.Cd, 0.0062, 1e-3, forceabs=True))
+ tests.add(CTest('Cdp', vsolver.Cdp, 0.0018, 1e-3, forceabs=True))
+ tests.add(CTest('xtr_top', vsolver.xtr[0], 0.061, 0.03, forceabs=True)) # Xfoil 0.056
+ tests.add(CTest('xtr_bot', vsolver.xtr[1], 0.740, 0.03, forceabs=True))
+ elif Re == 1e7 and M_inf == 0.5 and alpha == 5*math.pi/180:
+ tests.add(CTest('Cl', isolver.Cl, 0.65, 5e-2)) # Xfoil 0.69
+ tests.add(CTest('Cd', vsolver.Cd, 0.0067, 1e-3, forceabs=True))
+ tests.add(CTest('Cdp', vsolver.Cdp, 0.0025, 1e-3, forceabs=True))
+ tests.add(CTest('xtr_top', vsolver.xtr[0], 0.049, 0.03, forceabs=True)) # Xfoil 0.038
+ tests.add(CTest('xtr_bot', vsolver.xtr[1], 0.736, 0.03, forceabs=True))
+ elif Re == 1e7 and M_inf == 0 and alpha == 12*math.pi/180:
+ tests.add(CTest('Cl', isolver.Cl, 1.30, 5e-2)) # Xfoil 1.39
+ tests.add(CTest('Cd', vsolver.Cd, 0.011, 1e-3, forceabs=True))
+ tests.add(CTest('Cdp', vsolver.Cdp, 0.0064, 1e-3, forceabs=True))
+ tests.add(CTest('xtr_top', vsolver.xtr[0], 0.010, 0.03, forceabs=True)) # Xfoil 0.008
+ tests.add(CTest('xtr_bot', vsolver.xtr[1], 1.000, 0.03, forceabs=True))
+ elif Re==1e7 and M_inf==0 and alpha==5*math.pi/180 and user_xtr==[0.01, 0.3]:
+ tests.add(CTest('Cl', isolver.Cl, 0.55, 5e-2)) # Xfoil 0.5687
+ tests.add(CTest('Cd', vsolver.Cd, 0.0077, 1e-3, forceabs=True)) # Xfoil 0.00777
+ tests.add(CTest('Cdp', vsolver.Cdp, 0.002, 1e-3, forceabs=True))
+ tests.add(CTest('xtr_top', vsolver.xtr[0], 0.01, 1e-3, forceabs=True))
+ tests.add(CTest('xtr_bot', vsolver.xtr[1], 0.3, 1e-3, forceabs=True))
+ else:
+ raise Exception('Test not defined for this flow')
+
+ tests.run()
+
+ # eof
+ print('')
+
+if __name__ == "__main__":
+ main()
diff --git a/dart/tests/bliNonLift.py b/dart/tests/bliNonLift.py
new file mode 100755
index 0000000000000000000000000000000000000000..a67cd04e107983fa05353b8ef877174c67e03bef
--- /dev/null
+++ b/dart/tests/bliNonLift.py
@@ -0,0 +1,176 @@
+#!/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.VCoupler as floVC
+import dart.viscous.Drivers.VDriver as floVSMaster
+
+import dart.viscous.Solvers.Steady.StdCoupler as floVC_Std
+import dart.viscous.Solvers.Steady.StdSolver as floVS_Std
+
+import dart.viscous.Physics.VValidation as validation
+
+import 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 = 6e6
+ alpha = 0.*math.pi/180
+ M_inf = 0.15
+
+ #user_xtr=[0.41,0.41]
+ #user_xtr=[0,None]
+ user_xtr=[None,None]
+ user_Ti=None
+
+ mapMesh = 0
+
+ # Time solver parameters
+ Time_Domain = True # True for unsteady solver, False for steady solver
+ Cfl = 2
+ spaceMarching=True
+
+ if Time_Domain is True:
+ timeParams={'CFL0': Cfl, 'spaceMarching':spaceMarching}
+
+ # ========================================================================================
+
+ 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, timeParams, M_inf, case, mapMesh)
+ coupler = floVC.Coupler(isolver, vsolver, blw[0], blw[1], tol, gmshWriter, bnd, timeParams)
+
+ 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)
+
+ # 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)
+
+ # check results
+ print(ccolors.ANSI_BLUE + 'PyTesting...' + ccolors.ANSI_RESET)
+ tests = CTests()
+ if Re == 1e7 and M_inf == 0 and alpha == 5*math.pi/180:
+ tests.add(CTest('Cl', isolver.Cl, 0.55, 5e-2)) # Xfoil 0.58
+ tests.add(CTest('Cd', vsolver.Cd, 0.0062, 1e-3, forceabs=True))
+ tests.add(CTest('Cdp', vsolver.Cdp, 0.0018, 1e-3, forceabs=True))
+ tests.add(CTest('xtr_top', vsolver.xtr[0], 0.061, 0.03, forceabs=True)) # Xfoil 0.056
+ tests.add(CTest('xtr_bot', vsolver.xtr[1], 0.740, 0.03, forceabs=True))
+ elif Re == 1e7 and M_inf == 0.5 and alpha == 5*math.pi/180:
+ tests.add(CTest('Cl', isolver.Cl, 0.65, 5e-2)) # Xfoil 0.69
+ tests.add(CTest('Cd', vsolver.Cd, 0.0067, 1e-3, forceabs=True))
+ tests.add(CTest('Cdp', vsolver.Cdp, 0.0025, 1e-3, forceabs=True))
+ tests.add(CTest('xtr_top', vsolver.xtr[0], 0.049, 0.03, forceabs=True)) # Xfoil 0.038
+ tests.add(CTest('xtr_bot', vsolver.xtr[1], 0.736, 0.03, forceabs=True))
+ elif Re == 1e7 and M_inf == 0 and alpha == 12*math.pi/180:
+ tests.add(CTest('Cl', isolver.Cl, 1.30, 5e-2)) # Xfoil 1.39
+ tests.add(CTest('Cd', vsolver.Cd, 0.011, 1e-3, forceabs=True))
+ tests.add(CTest('Cdp', vsolver.Cdp, 0.0064, 1e-3, forceabs=True))
+ tests.add(CTest('xtr_top', vsolver.xtr[0], 0.010, 0.03, forceabs=True)) # Xfoil 0.008
+ tests.add(CTest('xtr_bot', vsolver.xtr[1], 1.000, 0.03, forceabs=True))
+ elif Re==1e7 and M_inf==0 and alpha==5*math.pi/180 and user_xtr==[0.01, 0.3]:
+ tests.add(CTest('Cl', isolver.Cl, 0.55, 5e-2)) # Xfoil 0.5687
+ tests.add(CTest('Cd', vsolver.Cd, 0.0077, 1e-3, forceabs=True)) # Xfoil 0.00777
+ tests.add(CTest('Cdp', vsolver.Cdp, 0.002, 1e-3, forceabs=True))
+ tests.add(CTest('xtr_top', vsolver.xtr[0], 0.01, 1e-3, forceabs=True))
+ tests.add(CTest('xtr_bot', vsolver.xtr[1], 0.3, 1e-3, forceabs=True))
+ else:
+ raise Exception('Test not defined for this flow')
+
+ tests.run()
+
+ # eof
+ print('')
+
+if __name__ == "__main__":
+ main()
diff --git a/dart/tests/bliSeparated.py b/dart/tests/bliSeparated.py
new file mode 100755
index 0000000000000000000000000000000000000000..fcae11d6f1ed0777bd790d4d10eb6a19156caf37
--- /dev/null
+++ b/dart/tests/bliSeparated.py
@@ -0,0 +1,173 @@
+#!/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.VCoupler as floVC
+import dart.viscous.Drivers.VDriver as floVSMaster
+
+import dart.viscous.Solvers.Steady.StdCoupler as floVC_Std
+import dart.viscous.Solvers.Steady.StdSolver as floVS_Std
+
+import dart.viscous.Physics.VValidation as validation
+
+import 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 = 1e6
+ alpha = 15.*math.pi/180
+ M_inf = 0.
+ #user_xtr=[0.012,1]
+ user_xtr=[None,None]
+ user_Ti=None
+
+ mapMesh = 0
+
+ # Time solver parameters
+ Time_Domain=True # True for unsteady solver, False for steady solver
+ Cfl = 1
+ spaceMarching=True
+ SteadyInitU=False # Initial condition prescribed by the steady solver
+
+ if Time_Domain is True:
+ timeParams={'CFL0': Cfl, 'SteadyInitU':SteadyInitU, 'spaceMarching':spaceMarching}
+
+ # ========================================================================================
+
+ 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.005, '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, timeParams, M_inf, case, mapMesh)
+ coupler = floVC.Coupler(isolver, vsolver, blw[0], blw[1], tol, gmshWriter, bnd, timeParams)
+
+ 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)
+
+ # 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)
+ # check results
+ print(ccolors.ANSI_BLUE + 'PyTesting...' + ccolors.ANSI_RESET)
+ tests = CTests()
+ if Re == 1e7 and M_inf == 0 and alpha == 5*math.pi/180:
+ tests.add(CTest('Cl', isolver.Cl, 0.55, 5e-2)) # Xfoil 0.58
+ tests.add(CTest('Cd', vsolver.Cd, 0.0062, 1e-3, forceabs=True))
+ tests.add(CTest('Cdp', vsolver.Cdp, 0.0018, 1e-3, forceabs=True))
+ tests.add(CTest('xtr_top', vsolver.xtr[0], 0.061, 0.03, forceabs=True)) # Xfoil 0.056
+ tests.add(CTest('xtr_bot', vsolver.xtr[1], 0.740, 0.03, forceabs=True))
+ elif Re == 1e7 and M_inf == 0.5 and alpha == 5*math.pi/180:
+ tests.add(CTest('Cl', isolver.Cl, 0.65, 5e-2)) # Xfoil 0.69
+ tests.add(CTest('Cd', vsolver.Cd, 0.0067, 1e-3, forceabs=True))
+ tests.add(CTest('Cdp', vsolver.Cdp, 0.0025, 1e-3, forceabs=True))
+ tests.add(CTest('xtr_top', vsolver.xtr[0], 0.049, 0.03, forceabs=True)) # Xfoil 0.038
+ tests.add(CTest('xtr_bot', vsolver.xtr[1], 0.736, 0.03, forceabs=True))
+ elif Re == 1e7 and M_inf == 0 and alpha == 12*math.pi/180:
+ tests.add(CTest('Cl', isolver.Cl, 1.30, 5e-2)) # Xfoil 1.39
+ tests.add(CTest('Cd', vsolver.Cd, 0.011, 1e-3, forceabs=True))
+ tests.add(CTest('Cdp', vsolver.Cdp, 0.0064, 1e-3, forceabs=True))
+ tests.add(CTest('xtr_top', vsolver.xtr[0], 0.010, 0.03, forceabs=True)) # Xfoil 0.008
+ tests.add(CTest('xtr_bot', vsolver.xtr[1], 1.000, 0.03, forceabs=True))
+ elif Re==1e7 and M_inf==0 and alpha==5*math.pi/180 and user_xtr==[0.01, 0.3]:
+ tests.add(CTest('Cl', isolver.Cl, 0.55, 5e-2)) # Xfoil 0.5687
+ tests.add(CTest('Cd', vsolver.Cd, 0.0077, 1e-3, forceabs=True)) # Xfoil 0.00777
+ tests.add(CTest('Cdp', vsolver.Cdp, 0.002, 1e-3, forceabs=True))
+ tests.add(CTest('xtr_top', vsolver.xtr[0], 0.01, 1e-3, forceabs=True))
+ tests.add(CTest('xtr_bot', vsolver.xtr[1], 0.3, 1e-3, forceabs=True))
+ else:
+ raise Exception('Test not defined for this flow')
+
+ tests.run()
+
+ # eof
+ print('')
+
+if __name__ == "__main__":
+ main()
diff --git a/dart/tests/bliTransonic.py b/dart/tests/bliTransonic.py
new file mode 100755
index 0000000000000000000000000000000000000000..bc3fd2eb61f880f17b39137baf7d094b1784818d
--- /dev/null
+++ b/dart/tests/bliTransonic.py
@@ -0,0 +1,180 @@
+#!/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.VCoupler as floVC
+import dart.viscous.Drivers.VDriver as floVSMaster
+
+import dart.viscous.Solvers.Steady.StdCoupler as floVC_Std
+import dart.viscous.Solvers.Steady.StdSolver as floVS_Std
+
+import dart.viscous.Physics.VValidation as validation
+
+import 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 = 2*math.pi/180
+ M_inf = 0.715
+
+ #user_xtr=[0.41,0.41]
+ #user_xtr=[0,None]
+ user_xtr=[None,None]
+ user_Ti=None
+
+ mapMesh = 0
+
+ # Time solver parameters
+ Time_Domain=True # True for unsteady solver, False for steady solver
+ Cfl = 1
+ spaceMarching=True
+ SteadyInitU=False # Initial condition prescribed by the steady solver
+
+ if Time_Domain is True:
+ timeParams={'CFL0': Cfl, 'SteadyInitU':SteadyInitU, 'spaceMarching':spaceMarching}
+
+ # ========================================================================================
+
+ U_inf = [math.cos(alpha), math.sin(alpha)] # norm should be 1
+ c_ref = 1
+ dim = 2
+ tol = 1e-3 # 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.01}
+ #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, timeParams, M_inf, case, mapMesh)
+ coupler = floVC.Coupler(isolver, vsolver, blw[0], blw[1], tol, gmshWriter, bnd, timeParams)
+
+ 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)
+
+ # 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)
+
+ # check results
+ print(ccolors.ANSI_BLUE + 'PyTesting...' + ccolors.ANSI_RESET)
+ tests = CTests()
+ if Re == 1e7 and M_inf == 0 and alpha == 5*math.pi/180:
+ tests.add(CTest('Cl', isolver.Cl, 0.55, 5e-2)) # Xfoil 0.58
+ tests.add(CTest('Cd', vsolver.Cd, 0.0062, 1e-3, forceabs=True))
+ tests.add(CTest('Cdp', vsolver.Cdp, 0.0018, 1e-3, forceabs=True))
+ tests.add(CTest('xtr_top', vsolver.xtr[0], 0.061, 0.03, forceabs=True)) # Xfoil 0.056
+ tests.add(CTest('xtr_bot', vsolver.xtr[1], 0.740, 0.03, forceabs=True))
+ elif Re == 1e7 and M_inf == 0.5 and alpha == 5*math.pi/180:
+ tests.add(CTest('Cl', isolver.Cl, 0.65, 5e-2)) # Xfoil 0.69
+ tests.add(CTest('Cd', vsolver.Cd, 0.0067, 1e-3, forceabs=True))
+ tests.add(CTest('Cdp', vsolver.Cdp, 0.0025, 1e-3, forceabs=True))
+ tests.add(CTest('xtr_top', vsolver.xtr[0], 0.049, 0.03, forceabs=True)) # Xfoil 0.038
+ tests.add(CTest('xtr_bot', vsolver.xtr[1], 0.736, 0.03, forceabs=True))
+ elif Re == 1e7 and M_inf == 0 and alpha == 12*math.pi/180:
+ tests.add(CTest('Cl', isolver.Cl, 1.30, 5e-2)) # Xfoil 1.39
+ tests.add(CTest('Cd', vsolver.Cd, 0.011, 1e-3, forceabs=True))
+ tests.add(CTest('Cdp', vsolver.Cdp, 0.0064, 1e-3, forceabs=True))
+ tests.add(CTest('xtr_top', vsolver.xtr[0], 0.010, 0.03, forceabs=True)) # Xfoil 0.008
+ tests.add(CTest('xtr_bot', vsolver.xtr[1], 1.000, 0.03, forceabs=True))
+ elif Re==1e7 and M_inf==0 and alpha==5*math.pi/180 and user_xtr==[0.01, 0.3]:
+ tests.add(CTest('Cl', isolver.Cl, 0.55, 5e-2)) # Xfoil 0.5687
+ tests.add(CTest('Cd', vsolver.Cd, 0.0077, 1e-3, forceabs=True)) # Xfoil 0.00777
+ tests.add(CTest('Cdp', vsolver.Cdp, 0.002, 1e-3, forceabs=True))
+ tests.add(CTest('xtr_top', vsolver.xtr[0], 0.01, 1e-3, forceabs=True))
+ tests.add(CTest('xtr_bot', vsolver.xtr[1], 0.3, 1e-3, forceabs=True))
+ else:
+ raise Exception('Test not defined for this flow')
+
+ tests.run()
+
+ # eof
+ print('')
+
+if __name__ == "__main__":
+ main()
diff --git a/dart/viscous/Driver.py b/dart/viscous/Drivers/VDriver.py
similarity index 59%
rename from dart/viscous/Driver.py
rename to dart/viscous/Drivers/VDriver.py
index fc7429f789cf0d0936a171c5a00793afb8948c1e..e0821fee7b90c9a89f215457d2984e362a91ffb0 100644
--- a/dart/viscous/Driver.py
+++ b/dart/viscous/Drivers/VDriver.py
@@ -14,21 +14,22 @@
# ------------------------------------------------------------------------------
# Imports
# ------------------------------------------------------------------------------
-from dart.viscous.model.Discretization import Discretization
-from dart.viscous.model.CSolver import sysMatrix
-from dart.viscous.GroupConstructor import GroupConstructor
-from dart.viscous.model.Variables import Variables
-import dart.viscous.model.TSolver as TSolver
-import dart.viscous.model.Transition as Transition
-from dart.viscous.model.vplo import vplo
-from dart.viscous.PostProcess import PostProcess
-import dart.viscous.Validation as validation
-import dart.viscous.model.BIConditions as Conditions
+from dart.viscous.Drivers.VGroupConstructor import GroupConstructor
+from dart.viscous.Drivers.VPostProcess import PostProcess
+
+import dart.viscous.Solvers.VTimeSolver as TSolver
+from dart.viscous.Solvers.VSolver import sysMatrix
+from dart.viscous.Solvers.VTransition import Transition
+
+from dart.viscous.Physics.VVariables import Variables
+from dart.viscous.Physics.VDiscretization import Discretization
+from dart.viscous.Physics.VValidation import Validation
+import dart.viscous.Physics.VBIConditions as Conditions
+
from fwk.coloring import ccolors
+
import numpy as np
import matplotlib.pyplot as plt
-import time
-import cProfile, pstats
# ------------------------------------------------------------------------------
# Driver : Input: - Data with the inviscid solution, geometry and user input for
@@ -41,7 +42,8 @@ import cProfile, pstats
# required data for the inviscid solver (blowing velocity aerodynamic drag,...)
# ------------------------------------------------------------------------------
class Driver:
- def __init__(self,_Re, user_xtr, _timeParams, _case, _mapMesh):
+ def __init__(self,_Re, user_xtr, _timeParams, _Minfty, _case, _mapMesh):
+ self.Minfty = _Minfty
self.Re=_Re
self.Cd=0
self.it=0
@@ -149,122 +151,124 @@ class Driver:
# Directly modifes the group that were initialy sent.
# ------------------------------------------------------------------------------------------
nFactor = 3
+
# Merge upper, lower sides and wake
dataIn=[None,None]
for n in range(0, len(group)):
group[n].connectListNodes, group[n].connectListElems, dataIn[n] = group[n].connectList()
- param, xx, xGlobal, bNodes, data, dx, initMesh, bNodesInitMesh = GroupConstructor().mergeVectors(dataIn, self.mapMesh, nFactor)
+ paramDict, bNodes, data, dx = GroupConstructor().mergeVectors(dataIn)
+
+ # Initialize variable container.
+ var = Variables(data, paramDict, bNodes, self.xtrF, self.Ti, self.Re)
- # Initialize variable container
- plotter=vplo(self.case)
- var=Variables(data, param, xx, xGlobal, bNodes, self.xtrF, self.Ti, self.Re)
if self.it!=0:
+ # Extenal flow local markers.
var.extPar[4::var.nExtPar]=data[:,8]
- DirichletBC=Conditions.Boundaries(self.Re) # Boundary condition (Airfoil/Wake).
+ # Boundary condition (Airfoil/Wake).
+ DirichletBC=Conditions.Boundaries(self.Re)
+
+ # Gradient computation w/ different schemes.
gradComp=Discretization(var.xx, var.extPar[4::var.nExtPar], var.bNodes) # Gradients computation.
- cSolver=sysMatrix(var.nN, var.nVar, gradComp.fou, gradComp.fou) # Spacial operator and Jacobian computation using the boundary layer laws.
- transSolver = Transition.Transition(var.xtrF)
+ # Boundary layer equations solver that can provide Residuals and Jacobian computation.
+ CSolver=sysMatrix(var.nN, var.nVar, gradComp.fou, gradComp.fou)
+
+ # Transition solver.
+ transSolver = Transition(var.xtrF)
+
# Initialize time Integrator
- tSolver=TSolver.implicitEuler(cSolver, transSolver, DirichletBC.wakeBC, self.timeParams['CFL0'])
+ tSolver=TSolver.implicitEuler(CSolver, transSolver, DirichletBC.wakeBC, self.timeParams['CFL0'], self.Minfty)
+
+ # Output setup at the first iteration.
+
+ if self.it == 0:
+ print('+------------------------------- Solver setup ---------------------------------+')
+ if self.mapMesh:
+ print('| - Mesh mapped from {:<4.0f} to {:<4.0f} elements.{:>38s}'.format(len(initMesh), var.nN, '|'))
+ print('| - Number of elements: {:<.0f}. {:>51s}'.format(var.nN,'|'))
+ if var.xtrF[0] is None:
+ print('| - Free transition on the upper side. {:>41}'.format('|'))
+ else:
+ print('| - Forced transition on the upper side; x/c = {:<.4f}. {:>25s}'.format(var.xtrF[0], '|'))
+ if var.xtrF[1] is None:
+ print('| - Free transition on the lower side. {:>41}'.format('|'))
+ else:
+ print('| - Forced transition on the lower side; x/c = {:<.4f}. {:>25s}'.format(var.xtrF[1], '|'))
+ print('| - Critical amplification ratio: {:<.2f}. {:>40s}'.format(var.Ncrit,'|'))
+
+ # Output solver parameters.
+ print('|{:>79}'.format('|'))
+ print('| Numerical parameters {:>57}'.format('|'))
+ print('| - CFL0: {:<3.2f} {:>65}'.format(tSolver.getCFL0(),'|'))
+ print('| - Tolerance: {:<3.0f} {:>61}'.format(np.log10(tSolver.gettol()),'|'))
+ print('| - Solver maximum iterations: {:<5.0f} {:>43}'.format(tSolver.getmaxIt(),'|'))
+ print('|{:>79}'.format('|'))
+
+ # Output preprocessing satut.
print('+------------------------------ Preprocessing ---------------------------------+')
- if self.mapMesh:
- print('| - Mesh mapped from {:<.0f} to {:<.0f} elements.{:>38s}'.format(len(initMesh), var.nN, '|'))
- print('| - Number of elements: {:<.0f}. {:>51s}'.format(var.nN,'|'))
+
print('| - Maximum Mach number: {:<.2f}. {:>49s}'.format(np.max((var.extPar[0::var.nExtPar])),'|'))
- if var.xtrF[0] is None:
- print('| - Free transition on the upper side. {:>41}'.format('|'))
- else:
- print('| - Forced transition on the upper side; x/c = {:<.4f}. {:>25s}'.format(var.xtrF[0], '|'))
- if var.xtrF[1] is None:
- print('| - Free transition on the lower side. {:>41}'.format('|'))
- else:
- print('| - Forced transition on the lower side; x/c = {:<.4f}. {:>25s}'.format(var.xtrF[1], '|'))
- print('| - Critical amplification ratio: {:<.2f}. {:>40s}'.format(var.Ncrit,'|'))
# Initial condition: 'Driver' stores the solution form the previous coupling iteration.
- if self.uPrevious is None or len(self.uPrevious) != len(var.u): # Typically for the first iteration
- initializer=Conditions.Initial(DirichletBC) # Initalize boundary condition
+
+ if self.uPrevious is None or self.it < 6:
+ # Typically for the first iteration.
+
+ # Initalize initial condition builder.
+ initializer=Conditions.Initial(DirichletBC)
+
+ # Compute initial condition based on Twaithes solution.
initializer.initialConditions(var) # Set boundary condition
- tSolver.setCFL0(0.5)
- tSolver.setmaxIt(20000)
+ # Default values of the time marching procedure if we start from an
+ # initial solution that is far from the steady state.
+ #tSolver.setCFL0(0.5)
+ #tSolver.setmaxIt(30000)
print('| - Using Twaithes solution as initial guess. {:>34}'.format('|'))
- #tSolver.tol=1e-2
- else: # If we use a solution previously obtained.
+
+ else:
+ # If we use a solution previously obtained.
var.u=self.uPrevious.copy()
- var.u[2::var.nVar] = 0
- print('| - Using previous solution as initial guess. {:>34}'.format('|'))
- print('|{:>79}'.format('|'))
- print('| Numerical parameters {:>57}'.format('|'))
- print('| - CFL0: {:<.2f} {:>65}'.format(tSolver.getCFL0(),'|'))
- print('| - Tolerance: {:<.0f} {:>62}'.format(np.log10(tSolver.gettol()),'|'))
- print('| - Solver maximum iterations: {:<.0f} {:>45}'.format(tSolver.getmaxIt(),'|'))
- print('|{:>79}'.format('|'))
+ # Reset amplification factors
+ var.u[2::var.nVar] = 0
- """plt.figure(1)
- plt.plot(var.xGlobal[:var.bNodes[1]],var.extPar[2:var.bNodes[1]*var.nExtPar:var.nExtPar])
- plt.plot(var.xGlobal[var.bNodes[1]:var.bNodes[2]],var.extPar[var.bNodes[1]*var.nExtPar + 2:var.bNodes[2]*var.nExtPar:var.nExtPar])
- plt.plot(var.xGlobal[var.bNodes[2]:],var.extPar[var.bNodes[2]*var.nExtPar + 2::var.nExtPar])
- plt.title('Inv Vel IN')
- plt.xlim([0,1])
-
- plt.figure(2)
- plt.plot(var.xGlobal[:var.bNodes[1]],var.extPar[0:var.bNodes[1]*var.nExtPar:var.nExtPar])
- plt.plot(var.xGlobal[var.bNodes[1]:var.bNodes[2]],var.extPar[var.bNodes[1]*var.nExtPar + 0:var.bNodes[2]*var.nExtPar:var.nExtPar])
- plt.plot(var.xGlobal[var.bNodes[2]:],var.extPar[var.bNodes[2]*var.nExtPar + 0::var.nExtPar])
- plt.title('Inv Mach')
- plt.xlim([0,1])
- plt.show()"""
+ print('| - Using previous solution as initial guess. {:>34}'.format('|'))
- print('+------------------------------- Solver begin ---------------------------------+')
# Run the time integration
- nTries=1
+ print('+------------------------------- Solver begin ---------------------------------+')
convergedPoints = tSolver.timeSolve(var)
print('+-------------------------------- Solver Exit ---------------------------------+')
- if np.any(convergedPoints!=1):
- print('|', ccolors.ANSI_YELLOW + 'Some points did not converge.', ccolors.ANSI_RESET, '{:>55s}'.format('|'))
+ # Output time solver status.
+ if np.any(convergedPoints != 1):
+ print('|', ccolors.ANSI_YELLOW + 'Some points did not converge.', ccolors.ANSI_RESET, '{:>45s}'.format('|'))
- elif np.all(convergedPoints==1):
+ elif np.all(convergedPoints == 1):
print('|', ccolors.ANSI_GREEN + 'All points converged.', ccolors.ANSI_RESET, '{:>55s}'.format('|'))
print('|{:>79}'.format('|'))
- # Save solution to use it as initial condition the next iteration
+ # Save solution to use it as initial condition the next iteration.
+
+ # Reset Ct in the laminar portion of the flow.
+ # (This quantity is not defined for a laminar flow)
var.u[4::var.nVar][var.flowRegime == 0] = 0.001
+
+ # Copy solution.
self.uPrevious=var.u.copy()
self.blParPrevious=var.blPar.copy()
- """plt.plot(var.xGlobal[0:var.bNodes[1]],var.u[0:var.bNodes[1]*var.nVar:var.nVar])
- plt.plot(var.xGlobal[var.bNodes[1]:var.bNodes[2]],var.u[var.bNodes[1]*var.nVar + 0:var.bNodes[2]*var.nVar:var.nVar])
- plt.plot(var.xGlobal[var.bNodes[2]:],var.u[var.bNodes[2]*var.nVar + 0::var.nVar])
- plt.show()
- plt.plot(var.xGlobal[0:var.bNodes[1]],var.u[1:var.bNodes[1]*var.nVar:var.nVar])
- plt.plot(var.xGlobal[var.bNodes[1]:var.bNodes[2]],var.u[var.bNodes[1]*var.nVar + 1:var.bNodes[2]*var.nVar:var.nVar])
- plt.plot(var.xGlobal[var.bNodes[2]:],var.u[var.bNodes[2]*var.nVar + 1::var.nVar])
- plt.show()
- plt.plot(var.xGlobal[0:var.bNodes[1]],var.u[2:var.bNodes[1]*var.nVar:var.nVar])
- plt.plot(var.xGlobal[var.bNodes[1]:var.bNodes[2]],var.u[var.bNodes[1]*var.nVar + 2:var.bNodes[2]*var.nVar:var.nVar])
- plt.plot(var.xGlobal[var.bNodes[2]:],var.u[var.bNodes[2]*var.nVar + 2::var.nVar])
- plt.show()
- plt.plot(var.xGlobal[0:var.bNodes[1]],var.u[3:var.bNodes[1]*var.nVar:var.nVar])
- plt.plot(var.xGlobal[var.bNodes[1]:var.bNodes[2]],var.u[var.bNodes[1]*var.nVar + 3:var.bNodes[2]*var.nVar:var.nVar])
- plt.plot(var.xGlobal[var.bNodes[2]:],var.u[var.bNodes[2]*var.nVar + 3::var.nVar])
- plt.show()
- plt.plot(var.xGlobal[0:var.bNodes[1]],var.u[4:var.bNodes[1]*var.nVar:var.nVar])
- plt.plot(var.xGlobal[var.bNodes[1]:var.bNodes[2]],var.u[var.bNodes[1]*var.nVar + 4:var.bNodes[2]*var.nVar:var.nVar])
- plt.plot(var.xGlobal[var.bNodes[2]:],var.u[var.bNodes[2]*var.nVar + 4::var.nVar])
- plt.show()"""
-
- self.blScalars, AFblwVel, WKblwVel, AFblVectors, WKblVectors = PostProcess().sortParam(var, self.mapMesh, initMesh, bNodesInitMesh, nFactor)
+ # Compute blowing velocities and sort the boundary layer parameters.
+ self.blScalars, AFblwVel, WKblwVel, AFblVectors, WKblVectors = PostProcess().sortParam(var)
self.xtr=var.xtr.copy()
self.blVectors=AFblVectors
+ # Group modification to ensure communication between the solvers.
+
group[0].TEnd = [AFblVectors[1][var.bNodes[1]-1], AFblVectors[1][var.bNodes[2]-1]]
group[0].HEnd = [AFblVectors[2][var.bNodes[1]-1], AFblVectors[2][var.bNodes[2]-1]]
group[0].CtEnd = [AFblVectors[8][var.bNodes[1]-1], AFblVectors[8][var.bNodes[2]-1]]
@@ -289,35 +293,6 @@ class Driver:
group[1].deltaStar = group[1].deltaStar[group[1].connectListNodes.argsort()]
group[1].xx = group[1].xx[group[1].connectListNodes.argsort()]
group[1].u = WKblwVel[group[1].connectListElems.argsort()]
-
+
"""self.writeFile()
- val=validation.Validation('Case1FreeTrans.dat')
- val.main(1,self.wData)"""
-
-
- """plt.plot(var.xGlobal[:var.bNodes[1]], var.u[3:var.bNodes[1]*var.nVar:var.nVar])
- plt.show()"""
-
- """"plt.plot(var.xGlobal[0:var.bNodes[1]],var.u[3:var.bNodes[1]*var.nVar:var.nVar])
- plt.plot(var.xGlobal[var.bNodes[1]:var.bNodes[2]],var.u[var.bNodes[1]*var.nVar + 3:var.bNodes[2]*var.nVar:var.nVar])
- plt.plot(var.xGlobal[var.bNodes[2]:],var.u[var.bNodes[2]*var.nVar + 3::var.nVar])
- plt.title('Ue')
- plt.show()
-
- plt.plot(var.xGlobal[0:var.bNodes[1]],var.blPar[9:var.bNodes[1]*var.nBlPar:var.nBlPar])
- plt.plot(var.xGlobal[var.bNodes[1]:var.bNodes[2]],var.blPar[var.bNodes[1]*var.nBlPar + 9:var.bNodes[2]*var.nBlPar:var.nBlPar])
- plt.plot(var.xGlobal[var.bNodes[2]:],var.blPar[var.bNodes[2]*var.nBlPar + 9::var.nBlPar])
- plt.title('deltaStar')
- plt.show()
-
- plt.plot(var.xGlobal[:var.bNodes[1]],var.extPar[1:var.bNodes[1]*var.nExtPar:var.nExtPar])
- plt.plot(var.xGlobal[var.bNodes[1]:var.bNodes[2]],var.extPar[var.bNodes[1]*var.nExtPar + 1:var.bNodes[2]*var.nExtPar:var.nExtPar])
- plt.plot(var.xGlobal[var.bNodes[2]:],var.extPar[var.bNodes[2]*var.nExtPar + 1::var.nExtPar])
- plt.title('rhoe')
- plt.show()
-
- plt.plot(var.xGlobal[:var.bNodes[1]],AFblwVel[:var.bNodes[1]])
- plt.plot(var.xGlobal[var.bNodes[1]:var.bNodes[2]],AFblwVel[var.bNodes[1]:var.bNodes[2]])
- plt.plot(var.xGlobal[var.bNodes[1]:var.bNodes[2]],AFblwVel[var.bNodes[1]:var.bNodes[2]])
- plt.title('BlwVel')
- plt.show()"""
\ No newline at end of file
+ Validation().main(1,self.wData, WKblVectors, var.xGlobal[var.bNodes[2]:])"""
\ No newline at end of file
diff --git a/dart/viscous/GroupConstructor.py b/dart/viscous/Drivers/VGroupConstructor.py
similarity index 64%
rename from dart/viscous/GroupConstructor.py
rename to dart/viscous/Drivers/VGroupConstructor.py
index b6fd94575fa95c288ff9402aa634f20964b6ac67..f93410eb0aab24a0d06c4e9226dd08cabb0a8f9e 100755
--- a/dart/viscous/GroupConstructor.py
+++ b/dart/viscous/Drivers/VGroupConstructor.py
@@ -26,34 +26,15 @@ class GroupConstructor():
def __init__(self) -> None:
pass
- def mergeVectors(self,dataIn, mapMesh, nFactor):
+ def mergeVectors(self,dataIn):
'''Merges 3 groups upper/lower sides and wake.'''
-
- initMesh = np.concatenate((dataIn[0][0][:,0], dataIn[0][1][1:,0], dataIn[1][:,0]))
- bNodesInitMesh = [0, len(dataIn[0][0][:,0]), len(dataIn[0][0][:,0])+len(dataIn[0][1][1:,0])]
for k in range(len(dataIn)):
-
if len(dataIn[k])==2: # Airfoil case
- if mapMesh:
- dataUpper = self.refineMesh(dataIn[k][0], nFactor)
- dataLower = self.refineMesh(dataIn[k][1], nFactor)
- else:
- dataUpper = dataIn[k][0]
- dataLower = dataIn[k][1]
-
-
- xx_up, dv_up, vtExt_up, _, alpha_up= self.__getBLcoordinates(dataUpper)
- xx_lw, dv_lw, vtExt_lw, _, alpha_lw= self.__getBLcoordinates(dataLower)
+ 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
-
- if mapMesh:
- dataWake = self.refineMesh(dataIn[k], nFactor)
- else:
- dataWake = dataIn[k]
-
- xx_wk, dv_wk, vtExt_wk, _, alpha_wk= self.__getBLcoordinates(dataWake)
-
- xx_wk += xx_up[-1]
+ xx_wk, dv_wk, vtExt_wk, _, alpha_wk= self.__getBLcoordinates(dataIn[k])
+
# Delete stagnation point doublon
xx_lw=np.delete(xx_lw,0)
dv_lw=np.delete(dv_lw,0)
@@ -61,33 +42,30 @@ class GroupConstructor():
alpha_lw=np.delete(alpha_lw,0)
# Nodes at the boundary of the computational domain: [Stagnation point, First lower side point, First wake point]
- boundaryNodes=np.empty(3, dtype=int)
- boundaryNodes[0] = 0
- boundaryNodes[1] = len(xx_up)
- boundaryNodes[2] = len(xx_up)+len(xx_lw)
+ boundaryNodes=np.zeros(3,dtype=int)
+ boundaryNodes[0]=0
+ boundaryNodes[1]=len(xx_up)
+ boundaryNodes[2]=len(xx_up)+len(xx_lw)
- param={'dv':None,'vtExt':None,'alpha':None}
+ param={'xx':None,'dv':None,'vtExt':None,'alpha':None}
+ param['xx']=np.concatenate((xx_up,xx_lw,xx_wk))
param['dv']=np.concatenate((dv_up,dv_lw,dv_wk))
param['vtExt']=np.concatenate((vtExt_up,vtExt_lw,vtExt_wk))
param['alpha']=np.concatenate((alpha_up,alpha_lw,alpha_wk))
- xx=np.concatenate((xx_up,xx_lw,xx_wk))
-
- data=np.zeros((len(xx),9))
- data=np.concatenate((dataUpper,dataLower,dataWake),axis=0)
+ data=np.zeros((len(param['xx']),9))
+ data=np.concatenate((dataIn[0][0],dataIn[0][1],dataIn[1]),axis=0)
data=np.delete(data,boundaryNodes[1],0) # Delete stagnation point doublon
- xGlobal = data[:,0]
-
- dx=np.zeros(len(xx))
- for i in range(1,len(xx)):
+ dx=np.zeros(len(param['xx']))
+ for i in range(1,len(param['xx'])):
if i==boundaryNodes[1]:
- dx[i]=xx[i]-xx[0]
+ dx[i]=param['xx'][i]-param['xx'][0]
elif i==boundaryNodes[2]:
dx[i]=0
else:
- dx[i]=xx[i]-xx[i-1]
- return param, xx, xGlobal, boundaryNodes, data, dx, initMesh, bNodesInitMesh
+ dx[i]=param['xx'][i]-param['xx'][i-1]
+ return param, boundaryNodes, data, dx
def __getBLcoordinates(self,data):
@@ -107,23 +85,4 @@ class GroupConstructor():
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 refineMesh(self, inMesh, nFactor):
-
- outMesh = np.empty(len(inMesh[0,:]))
-
- for marker in range(len(inMesh[:,0]) - 1):
- dCell = inMesh[marker + 1, :] - inMesh[marker, :]
-
- for n in range(nFactor):
-
- outMesh = np.ma.row_stack((outMesh, inMesh[marker, :] + n*dCell/nFactor))
-
- outMesh = np.delete(outMesh, 0, 0)
-
- # Add last cell
- outMesh = np.ma.row_stack((outMesh, inMesh[-1, :]))
-
- return outMesh
\ No newline at end of file
+ return xx, dv, vt, nN, alpha
\ No newline at end of file
diff --git a/dart/viscous/PostProcess.py b/dart/viscous/Drivers/VPostProcess.py
similarity index 62%
rename from dart/viscous/PostProcess.py
rename to dart/viscous/Drivers/VPostProcess.py
index b3c11a7198b3e68128b90350ad4bbdb716918adb..0cdfd3b4d053f868c0ae46bdf3c44873d29c1699 100644
--- a/dart/viscous/PostProcess.py
+++ b/dart/viscous/Drivers/VPostProcess.py
@@ -1,53 +1,31 @@
from matplotlib import pyplot as plt
import numpy as np
+
class PostProcess:
def __init__(self):
pass
- def sortParam(self,var, mapMesh, initMesh, bNodesInitMesh, nFactor):
+ def sortParam(self,var):
# Recompute deltaStar.
var.blPar[9::var.nBlPar] = var.u[0::var.nVar] * var.u[1::var.nVar]
# Compute blowing velocity on each cell.
- blwVel = np.empty(var.nN - 1)
-
- for i in range(1, var.nN): # -1 and we will remove the first wake point after.
- iPrev = 0 if i == var.bNodes[1] else i-1
-
- # blwVel[i-1] = (rhoe[i]*vt[i]*deltaStar[i]-rhoe[i-1]*vt[i-1]*deltaStar[i-1])/(rhoe[i]*(xx[i]-xx[i-1]))
-
- if i == var.bNodes[2]:
- sep = i-1 # Separation between airfoil and wake. (-1 because the first wake point will be removed)
- continue
- else:
- blwVel[i-1] = (var.extPar[i*var.nExtPar+1]*var.u[i*var.nVar+3] * var.blPar[i*var.nBlPar+9]
- - var.extPar[iPrev*var.nExtPar+1] * var.u[iPrev*var.nVar+3] * var.blPar[iPrev*var.nBlPar+9]) / (var.extPar[i*var.nExtPar+1] * (var.xx[i] - var.xx[iPrev]))
+ blwVel=np.zeros(var.nN-1)
- blwVel = np.delete(blwVel, var.bNodes[2] - 1)
+ for i in range(1, var.nN):
- if mapMesh:
- # Map blowing velocities to inviscid solver mesh.
+ iPrev=0 if i==var.bNodes[1] else i-1
- invBlwVel = np.zeros(len(initMesh)-2)
- for invMarker in range(len(invBlwVel)):
- if invMarker == bNodesInitMesh[2]:
- sep = invMarker - 1
- invBlwVel[invMarker] = 1/4*(blwVel[nFactor*invMarker - 1] + 2* blwVel[nFactor*invMarker] + blwVel[nFactor*invMarker + 1])
-
- AFblwVel = invBlwVel[:sep]
- WKblwVel = invBlwVel[sep:]
-
- else:
- # Separate airfoil and wake blowing velocities.
- AFblwVel=blwVel[:sep] # -1 because the first wake point was removed and indices have moved.
- WKblwVel=blwVel[sep:]
- """plt.figure(1)
- plt.plot(AFblwVel)
- plt.figure(2)
- plt.plot(WKblwVel)
- plt.show()"""
+ # blwVel[i-1] = (rhoe[i]*vt[i]*deltaStar[i]-rhoe[i-1]*vt[i-1]*deltaStar[i-1])/(rhoe[i]*(xx[i]-xx[i-1]))
+ blwVel[i-1] = (var.extPar[i*var.nExtPar+1]*var.u[i*var.nVar+3] * var.blPar[i*var.nBlPar+9]
+ - var.extPar[iPrev*var.nExtPar+1] * var.u[iPrev*var.nVar+3] * var.blPar[iPrev*var.nBlPar+9]) / (var.extPar[i*var.nExtPar+1] * (var.xx[i] - var.xx[iPrev]))
+
+ blwVel = np.delete(blwVel, var.bNodes[2])
+ # Separate airfoil and wake blowing velocities.
+ AFblwVel=blwVel[:var.bNodes[2]-1]
+ WKblwVel=blwVel[var.bNodes[2]-1:]
# Normalize Friction and dissipation coefficients.
frictionCoeff=var.u[3::var.nVar]**2 * var.blPar[2::var.nBlPar]
diff --git a/dart/viscous/airfoil.py b/dart/viscous/Master/VAirfoil.py
similarity index 99%
rename from dart/viscous/airfoil.py
rename to dart/viscous/Master/VAirfoil.py
index 31a579658a65b603e7634830a728adb449ef83b1..226091ddc57082c283e2dd2c3624d94af1de67ff 100755
--- a/dart/viscous/airfoil.py
+++ b/dart/viscous/Master/VAirfoil.py
@@ -19,7 +19,8 @@
## Airfoil around which the boundary layer is computed
# Amaury Bilocq
-from dart.viscous.Boundary import Boundary
+from dart.viscous.Master.VBoundary import Boundary
+
import numpy as np
class Airfoil(Boundary):
diff --git a/dart/viscous/Steady/StdBoundary.py b/dart/viscous/Master/VBoundary.py
similarity index 100%
rename from dart/viscous/Steady/StdBoundary.py
rename to dart/viscous/Master/VBoundary.py
diff --git a/dart/viscous/coupler.py b/dart/viscous/Master/VCoupler.py
similarity index 73%
rename from dart/viscous/coupler.py
rename to dart/viscous/Master/VCoupler.py
index aa4ce37621bd9e68dd5c17737cc9bc135e8ed85b..dcea5fb5054ce5d6929c411d860edcb3ea71f27a 100755
--- a/dart/viscous/coupler.py
+++ b/dart/viscous/Master/VCoupler.py
@@ -19,12 +19,14 @@
## Viscous-inviscid coupler (quasi-simultaneous coupling)
# Amaury Bilocq
+from dart.viscous.Master.VAirfoil import Airfoil
+from dart.viscous.Master.VWake import Wake
+
import numpy as np
-from fwk.coloring import ccolors
-from dart.viscous.Airfoil import Airfoil
-from dart.viscous.Wake import Wake
+
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, _timeParam=None):
@@ -36,24 +38,33 @@ class Coupler:
self.tol = _tol # tolerance of the VII
self.writer = _writer
self.timeParam=_timeParam
+
def run(self):
- ''' Perform coupling
- '''
+ """Viscous inviscid coupling.
+ """
+
# Initialize loop
it = 0
converged = 0 # temp
CdOld = self.vsolver.Cd
+
print('+---------------------------- VII Solver begin --------------------------------+')
print('+------------------------------------------------------------------------------+')
print('|{:>79s}'.format('|'))
+
while converged == 0:
- print('|', ccolors.ANSI_BLUE + 'Iteration: ', it, ccolors.ANSI_RESET, '{:>63s}'.format('|'))
- # run inviscid solver
+
+ print('|', ccolors.ANSI_BLUE + 'Iteration: {:<3.0f}'.format(it), ccolors.ANSI_RESET, '{:>62s}'.format('|'))
+
+ # Run inviscid solver
print('+----------------------------- Inviscid solver --------------------------------+')
self.isolver.run()
+
for n in range(0, len(self.group)):
- # print('Computing for', self.group[n].name, '...', end = ' ')
+
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]
@@ -64,41 +75,52 @@ class Coupler:
# 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', '')
- Cp = floU.extract(self.bnd.groups[0].tag.elems, self.isolver.Cp)
- tboxU.write(Cp, 'Cp_airfoil.dat', '%1.5e', ', ', 'x, y, z, Cp', '')
- # run viscous solver
+
+ # Run viscous solver
print('+------------------------------ Viscous solver --------------------------------+')
self.vsolver.run(self.group)
- self.vsolver.writeFile()
+
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])
- # Write files
- '''Cp = floU.extract(self.bnd.groups[0].tag.elems, self.isolver.Cp)
- tboxU.write(Cp, 'Cp_airfoil.dat', '%1.5e', ', ', 'x, y, z, Cp', '')
- self.isolver.save(self.writer)'''
- #self.vsolver.writeFile()
- # Output coupling statut
+
+ # Output coupling iteration.
xtrT_type='forced' if self.vsolver.xtrF[0] is not None else 'free'
xtrB_type='forced' if self.vsolver.xtrF[1] is not None else 'free'
- print('+---------------------------- Postprocessing ----------------------------------+')
- print('|{0:>6s}| {1:>14s}| {2:>18s}| {3:>18s}| {4:>14s}|'.format('Iter', 'Cd', f'Top xtr ({xtrT_type})', f'Bot xtr ({xtrB_type})','Error'))
- print('|{0:>6d}| {1:>14f}| {2:>18f}| {3:>18f}| {4:>14f}|'.format(it, self.vsolver.Cd, self.vsolver.xtr[0], self.vsolver.xtr[1], np.log10(abs(self.vsolver.Cd - CdOld)/self.vsolver.Cd)))
+
+ print('+------------------------------- Postprocessing -------------------------------+')
+
+ error = abs(self.vsolver.Cd - CdOld)/self.vsolver.Cd
+
+ print('|{0:>6s}| {1:>14s}| {2:>18s}| {3:>18s}| {4:>14s}|'.format('Iter', 'Cd',
+ f'Top xtr ({xtrT_type})',
+ f'Bot xtr ({xtrB_type})',
+ 'Error'))
+
+ print('|{0:>6d}| {1:>14f}| {2:>18f}| {3:>18f}| {4:>14f}|'.format(it, self.vsolver.Cd,
+ self.vsolver.xtr[0], self.vsolver.xtr[1],
+ np.log10(error)))
print('|{:>79}'.format('|'))
- print('+------------------------------------------------------------------------------+')
- # Write files
- # Converged or not
+ # Check convergence.
if abs(self.vsolver.Cd - CdOld)/self.vsolver.Cd < self.tol:
converged = 1
else:
CdOld = self.vsolver.Cd
+ # Prepare next iteration.
it += 1
self.vsolver.it += 1
- # save results
+
+ # Save results
+
+ # Inviscid solver files.
self.isolver.save(self.writer)
+
+ # Viscous solver files.
self.vsolver.writeFile()
+
+ # Pressure coefficient.
Cp = floU.extract(self.bnd.groups[0].tag.elems, self.isolver.Cp)
tboxU.write(Cp, 'Cp_airfoil.dat', '%1.5e', ', ', 'x, y, z, Cp', '')
print('\n')
diff --git a/dart/viscous/wake.py b/dart/viscous/Master/VWake.py
similarity index 97%
rename from dart/viscous/wake.py
rename to dart/viscous/Master/VWake.py
index f62982b1f43b878ae0e4ba0c8ad369043c58b847..3a6cfe47dfdda6020ac853f493355ea2ab835143 100755
--- a/dart/viscous/wake.py
+++ b/dart/viscous/Master/VWake.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.VBoundary import Boundary
+
import numpy as np
class Wake(Boundary):
diff --git a/dart/viscous/model/BIConditions.py b/dart/viscous/Physics/VBIConditions.py
similarity index 99%
rename from dart/viscous/model/BIConditions.py
rename to dart/viscous/Physics/VBIConditions.py
index 8ecbbff3ee6c06e3a7e699c5144e1eaf535da16e..4f8de9ccaff468eac9e6e1615c7e9cd7e277cb6d 100755
--- a/dart/viscous/model/BIConditions.py
+++ b/dart/viscous/Physics/VBIConditions.py
@@ -14,7 +14,8 @@
# ------------------------------------------------------------------------------
# Imports
# ------------------------------------------------------------------------------
-from dart.viscous.model.Closures import Closures
+from dart.viscous.Physics.VClosures import Closures
+
import matplotlib.pyplot as plt
import numpy as np
diff --git a/dart/viscous/model/Closures.py b/dart/viscous/Physics/VClosures.py
similarity index 100%
rename from dart/viscous/model/Closures.py
rename to dart/viscous/Physics/VClosures.py
diff --git a/dart/viscous/model/DataStructures.py b/dart/viscous/Physics/VDataStructures.py
similarity index 95%
rename from dart/viscous/model/DataStructures.py
rename to dart/viscous/Physics/VDataStructures.py
index 6ea77a8bc3cd5e718474288a84ad554db0c3692a..7e9903e2ec067379975e9dcf005c646cbce54779 100755
--- a/dart/viscous/model/DataStructures.py
+++ b/dart/viscous/Physics/VDataStructures.py
@@ -91,10 +91,10 @@ class DOuterFlow:
def __init__(self, data) -> None:
- self.__Mach =
- self.__density =
- self.__velocity =
- self.__dv =
+ 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 =
diff --git a/dart/viscous/model/Discretization.py b/dart/viscous/Physics/VDiscretization.py
similarity index 100%
rename from dart/viscous/model/Discretization.py
rename to dart/viscous/Physics/VDiscretization.py
diff --git a/dart/viscous/Validation.py b/dart/viscous/Physics/VValidation.py
similarity index 65%
rename from dart/viscous/Validation.py
rename to dart/viscous/Physics/VValidation.py
index dc2f13fc2fc547c9cdfd6ee1b407573d916368ca..5d2777316a70e8e570c2c2d0ba155ca4e4171dc5 100755
--- a/dart/viscous/Validation.py
+++ b/dart/viscous/Physics/VValidation.py
@@ -1,9 +1,26 @@
+################################################################################
+# #
+# 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):
+ def main(self, dartcl, wData, wakeblVectors = None, xWake = None):
# Load xfoil data
xFoilIn = 0
@@ -13,12 +30,12 @@ class Validation:
if self.xfoilFileName is not None:
xFoilIn = 1
- with open('/Users/pauldechamps/lab/dartflo/dart/viscous/XfoilFiles/' + self.xfoilFileName, 'r') as f:
+ 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/XfoilFiles/'+self.xfoilFileName,skiprows=3,usecols=(1,3,4,5,6,7))
+ 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]]
@@ -43,50 +60,49 @@ class Validation:
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.xlim([0,1])
plt.xlabel('$x/c$')
plt.ylabel('$\delta^*$')
- """ # Theta
- plot_Theta=plt.figure(2)
- plt.plot(wData['x'][:xTeDart], wData['theta'][:xTeDart],'b-',label=r'DARTflo upper')
- plt.plot(wData['x'][xTeDart:], wData['theta'][xTeDart:],'r-',label=r'DARTflo lower')
- plt.plot(data[:xStag,0],data[:xStag,2],'b--',label=r'Xfoil upper')
- plt.plot(data[xStag:,0],data[xStag:,2],'r--',label=r'Xfoil lower')
- #plt.plot(data[xWake:,0],data[xWake:,2],'r--',label=r'Xfoil Wake')
- plt.legend()
- plt.xlim([0,1])
- plt.xlabel('$x/c$')
- plt.ylabel('$\theta$')"""
+
# 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.axvline(x=0.0138)
+ 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.xlim([0,1])
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.xlim([0,1])
plt.xlabel('$x/c$')
plt.ylabel('$C_f$')
diff --git a/dart/viscous/model/Variables.py b/dart/viscous/Physics/VVariables.py
similarity index 89%
rename from dart/viscous/model/Variables.py
rename to dart/viscous/Physics/VVariables.py
index 356e3b615ed92d79818a57cbcee3748bd6c50749..a5a9f62d050f1a5e8e7c803e391580321e7d7f42 100755
--- a/dart/viscous/model/Variables.py
+++ b/dart/viscous/Physics/VVariables.py
@@ -16,6 +16,7 @@
# ------------------------------------------------------------------------------
import numpy as np
from matplotlib import pyplot as plt
+
class Variables:
"""Data container for viscous solver.
@@ -64,10 +65,10 @@ class Variables:
- Critical amplification ratio for the e^N method (default value= 9.)
"""
- def __init__(self, data, paramDict, _xx, _xGlobal, bNodes, xtrF, Ti, Re):
+ def __init__(self, data, paramDict, bNodes, xtrF, Ti, Re):
self.__Ti=Ti
self.Re=Re
- self.xx = _xx
+ self.xx=paramDict['xx']
self.nN=len(self.xx)
self.bNodes=bNodes
@@ -94,16 +95,9 @@ class Variables:
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]=_xx
+ self.extPar[4::self.nExtPar]=paramDict['xx']
self.dv=paramDict['dv']
- self.xGlobal = _xGlobal
-
-
- """plt.plot(self.xGlobal[0:self.bNodes[1]],self.extPar[2:self.bNodes[1]*self.nExtPar:self.nExtPar])
- plt.plot(self.xGlobal[self.bNodes[1]:self.bNodes[2]],self.extPar[self.bNodes[1]*self.nExtPar + 2 : self.bNodes[2]*self.nExtPar:self.nExtPar])
- plt.plot(self.xGlobal[self.bNodes[2]:],self.extPar[self.bNodes[2]*self.nExtPar + 2 ::self.nExtPar])
- plt.axvline(x=0.59194)
- plt.show()"""
+ self.xGlobal=data[:,0]
self.flowRegime = np.zeros(self.nN, dtype=int)
diff --git a/dart/viscous/XfoilFiles/BlParM03A0Re6e6.dat b/dart/viscous/Physics/XfoilFiles/BlParM03A0Re6e6.dat
similarity index 100%
rename from dart/viscous/XfoilFiles/BlParM03A0Re6e6.dat
rename to dart/viscous/Physics/XfoilFiles/BlParM03A0Re6e6.dat
diff --git a/dart/viscous/XfoilFiles/Case1FreeTrans.dat b/dart/viscous/Physics/XfoilFiles/Case1FreeTrans.dat
similarity index 100%
rename from dart/viscous/XfoilFiles/Case1FreeTrans.dat
rename to dart/viscous/Physics/XfoilFiles/Case1FreeTrans.dat
diff --git a/dart/viscous/XfoilFiles/case12Xfoil.dat b/dart/viscous/Physics/XfoilFiles/case12Xfoil.dat
similarity index 100%
rename from dart/viscous/XfoilFiles/case12Xfoil.dat
rename to dart/viscous/Physics/XfoilFiles/case12Xfoil.dat
diff --git a/dart/viscous/XfoilFiles/case15Xfoil.dat b/dart/viscous/Physics/XfoilFiles/case15Xfoil.dat
similarity index 100%
rename from dart/viscous/XfoilFiles/case15Xfoil.dat
rename to dart/viscous/Physics/XfoilFiles/case15Xfoil.dat
diff --git a/dart/viscous/README.md b/dart/viscous/README.md
old mode 100755
new mode 100644
index c2eed4b79865672d23e6876f195e0957277e57f9..fe02415fd61af94f8d874bc4f07928a91dca0e23
--- a/dart/viscous/README.md
+++ b/dart/viscous/README.md
@@ -1,21 +1,25 @@
-# viscous
-A dart module solving the boundary layer equations to perform viscous-inviscid coupling with the main full potential solver.
+# DARTFlo Viscous-Inviscid Interaction.
+The folder dartflo/dart/viscous/ is dedicated to the viscous inviscid interaction implemented in DARTflo. This section extends DARTflo capabilities to two dimensional steady high Reynolds number transitional viscous flows.
+The solver makes use of the integral boundary layer equations (IBL) and the e^N method to correct the inviscid flow solution computed by DARTflo. The solver is divded into two parts. One that solves the steady state IBL and one that deals with the time-dependent IBL solved through a time marching procedure. The communication between the solvers is always done through steady-state solutions. The interaction is able to deal with midly separated flows for which the boundary layer assumptions still hold.
+The two flow domains are coupled through a quasi-simulataneous interaction law that prevents Goldstein singularity.
-
+
## Main features
-* Solvers
- - [x] Steady: Solves the steady boundary layer equations
- - [x] Unsteady: Solves the unsteady boundary layer equations using (pseudo-)time marching procedures
-* Time-Marching
- - [x] Explicit/Implicit methods
- - [x] Newton method with direct linear solver for implicit time integration
-* Various
- - [x] Different spacial discretizations (FOU, SOU, CENTERED)
- - [x] Forced/Free transition
- - [x] Free-stream turbulence intensity (affects transition location)
- - [x] Adaptative CFL number
+* 2D viscous flow.
+* Free/forced laminar to turbulent transition.
+* Solver models
+ - Steady state equations solution.
+ - Time-dependant equations solution through a pseudo-time marching procedure.
-## Sub-modules
-* [Steady](flow/viscous/Steady): Steady equations solver and components
-* [Model](flow/viscous/model): Mathematical model of the time dependent boundary layer equations
\ No newline at end of file
+## Solver operation map
+The viscous inviscid interaction can be launched from /dart/tests/bli_xxx.py. The class in Master/Coupler.py is first called with the required arguments and manages the two solvers to work togheter. The steady-state equations solver is called from /Solvers/Steady/StdSolver.py. The time-dependent equations solver is called from /Drivers/VDriver.py which initializes the different components of the code and starts the viscous flow computation using a implicit Euler pseudo-time marching procedures.
+
+## References
+Crovato Adrien, [Steady Transonic Aerodynamic and Aeroelastic Modeling for Preliminary Aircraft Design](http://hdl.handle.net/2268/251906), PhD thesis, University of Liège, 2020.
+
+Crovato Adrien, [DARTFlo - Discrete Adjoint for Rapid Transonic Flows](http://hdl.handle.net/2268/264284), Technical note, University of Liège, 2021.
+
+Crovato A., Boman R., et al., [A Full Potential Static Aeroelastic Solver for Preliminary Aircraft Design](http://hdl.handle.net/2268/237955), 18th International Forum on Aeroelasticity and Structural Dynamics, IFASD 2019.
+
+Bilocq Amaury, [Implementation of a viscous-inviscid interaction scheme in a finite element full potential solver](http://hdl.handle.net/2268/252195), Master thesis, University of Liège, 2020.
diff --git a/dart/viscous/Steady/Newton.py b/dart/viscous/Solvers/Steady/Newton.py
similarity index 100%
rename from dart/viscous/Steady/Newton.py
rename to dart/viscous/Solvers/Steady/Newton.py
diff --git a/dart/viscous/Steady/StdAirfoil.py b/dart/viscous/Solvers/Steady/StdAirfoil.py
similarity index 99%
rename from dart/viscous/Steady/StdAirfoil.py
rename to dart/viscous/Solvers/Steady/StdAirfoil.py
index 18db9d5f61d9922cf8066ad44d5cb5b642940182..514e4681453473a376e1379044e683c7d06d7940 100755
--- a/dart/viscous/Steady/StdAirfoil.py
+++ b/dart/viscous/Solvers/Steady/StdAirfoil.py
@@ -19,7 +19,7 @@
## Airfoil around which the boundary layer is computed
# Amaury Bilocq
-from dart.viscous.Steady.StdBoundary import Boundary
+from dart.viscous.Solvers.Steady.StdBoundary import Boundary
import numpy as np
class Airfoil(Boundary):
diff --git a/dart/viscous/boundary.py b/dart/viscous/Solvers/Steady/StdBoundary.py
similarity index 100%
rename from dart/viscous/boundary.py
rename to dart/viscous/Solvers/Steady/StdBoundary.py
diff --git a/dart/viscous/Steady/StdCoupler.py b/dart/viscous/Solvers/Steady/StdCoupler.py
similarity index 97%
rename from dart/viscous/Steady/StdCoupler.py
rename to dart/viscous/Solvers/Steady/StdCoupler.py
index 9427e7144375a5bc7a1dd5a262e2242ae9906507..5b6613218b129faa5283516e7ff825e62b2eaa03 100755
--- a/dart/viscous/Steady/StdCoupler.py
+++ b/dart/viscous/Solvers/Steady/StdCoupler.py
@@ -21,8 +21,8 @@
import numpy as np
from fwk.coloring import ccolors
-from dart.viscous.Steady.StdAirfoil import Airfoil
-from dart.viscous.Steady.StdWake import Wake
+from dart.viscous.Solvers.Steady.StdAirfoil import Airfoil
+from dart.viscous.Solvers.Steady.StdWake import Wake
class Coupler:
def __init__(self, _isolver, _vsolver, _boundaryAirfoil, _boundaryWake, _tol, _writer):
diff --git a/dart/viscous/Steady/StdSolver.py b/dart/viscous/Solvers/Steady/StdSolver.py
similarity index 99%
rename from dart/viscous/Steady/StdSolver.py
rename to dart/viscous/Solvers/Steady/StdSolver.py
index 922bfa35f3bbf372aec04648fe8ab49e4d03848b..9185fa4a4c7178379c6dca2f3c4d2b9ae051ca4f 100755
--- a/dart/viscous/Steady/StdSolver.py
+++ b/dart/viscous/Solvers/Steady/StdSolver.py
@@ -42,7 +42,7 @@
#from dart.viscous.Steady.StdBoundary import Boundary
#from dart.viscous.Steady.StdWake import Wake
from matplotlib import pyplot as plt
-import dart.viscous.Steady.Newton as Newton
+import dart.viscous.Solvers.Steady.Newton as Newton
import numpy as np
from fwk.coloring import ccolors
@@ -438,8 +438,6 @@ class Solver:
Cd[i-1], Cf[i-1], Hstar[i-1], Hstar2[i-1], Hk[i-1], delta[i-1] = self.__laminarClosure(theta[i-1], H[i-1], Ret[i-1],Me[i-1], group.name)
Cd[i], Cf[i], Hstar[i], Hstar2[i],Hk[i],Cteq[i], delta[i] = self.__turbulentClosure(theta[i], H[i],Ct[i], Ret[i],Me[i], group.name)
deltaStar[i] = H[i]*theta[i]
- if group.name=='wake':
- print('dx = ',xx[i]-xx[i-1])
blwVel[i-1] = (rhoe[i]*vt[i]*deltaStar[i]-rhoe[i-1]*vt[i-1]*deltaStar[i-1])/(rhoe[i]*(xx[i]-xx[i-1]))
# Normalize the friction and dissipation coefficients by the freestream velocity
Cf = vt**2*Cf
@@ -451,7 +449,6 @@ class Solver:
# Outputs
blScalars = [Cdtot, Cdf, Cdp]
blVectors = [deltaStar, theta, H, Hk, Hstar, Hstar2, Cf, Cd, Ct, delta]
- print(len(deltaStar))
return blwVel, xtr, xx, blScalars, blVectors
@@ -468,6 +465,7 @@ class Solver:
self.blScalars = np.add(ublScalars, lblScalars)
self.blVectors = np.hstack((ublVectors, lblVectors))
blwVel = np.concatenate((ublwVel, lblwVel))
+ print('AF',len(blwVel))
self.xtr = [xtrTop, xtrBot]
# Boundary conditions for the wake
group.TEnd = [ublVectors[1][-1], lblVectors[1][-1]]
@@ -481,6 +479,7 @@ class Solver:
# Compute BLE for the wake
else:
blwVel, _, group.xx, blScalars, blVectors = self.__boundaryLayer(data, group)
+ print('WK',len(blwVel))
group.deltaStar = blVectors[0]
# The drag is measured at the end of the wake
self.Cd = blScalars[0]
diff --git a/dart/viscous/Steady/StdWake.py b/dart/viscous/Solvers/Steady/StdWake.py
similarity index 96%
rename from dart/viscous/Steady/StdWake.py
rename to dart/viscous/Solvers/Steady/StdWake.py
index 625b3e4ed4f60b5cbcb981f255d8ca2c40f0f620..96a0e8286bb6e8a31091d17adef8cea8318dfb47 100755
--- a/dart/viscous/Steady/StdWake.py
+++ b/dart/viscous/Solvers/Steady/StdWake.py
@@ -19,8 +19,8 @@
## Wake behind airfoil (around which the boundary layer is computed)
# Amaury Bilocq
-from dart.viscous.Steady.StdBoundary import Boundary
-from dart.viscous.Steady.StdAirfoil import Airfoil
+from dart.viscous.Solvers.Steady.StdBoundary import Boundary
+from dart.viscous.Solvers.Steady.StdAirfoil import Airfoil
import numpy as np
class Wake(Boundary):
diff --git a/dart/viscous/model/CSolver.py b/dart/viscous/Solvers/VSolver.py
similarity index 94%
rename from dart/viscous/model/CSolver.py
rename to dart/viscous/Solvers/VSolver.py
index 0ae15f2eb2d7d51e43904809286dbb0010a2e105..5c3a0f290ad99b585eb51caf9ceaa74cb7539732 100644
--- a/dart/viscous/model/CSolver.py
+++ b/dart/viscous/Solvers/VSolver.py
@@ -14,8 +14,9 @@
# ------------------------------------------------------------------------------
# Imports
# ------------------------------------------------------------------------------
+from dart.viscous.Physics.VClosures import Closures
+
import numpy as np
-from dart.viscous.model.Closures import Closures
class flowEquations:
@@ -61,12 +62,16 @@ class flowEquations:
if dataCont.flowRegime[iMarker]==0:
# Laminar closure relations
- dataCont.blPar[iMarker*dataCont.nBlPar+1:iMarker*dataCont.nBlPar+7]=Closures.laminarClosure(x[0], x[1], dataCont.blPar[iMarker*dataCont.nBlPar+0], dataCont.extPar[iMarker*dataCont.nExtPar+0], name)
+ dataCont.blPar[iMarker*dataCont.nBlPar+1:iMarker*dataCont.nBlPar+7]=Closures.laminarClosure(x[0], x[1],
+ dataCont.blPar[iMarker*dataCont.nBlPar+0],
+ dataCont.extPar[iMarker*dataCont.nExtPar+0], name)
elif dataCont.flowRegime[iMarker]==1:
# Turbulent closures relations
- dataCont.blPar[iMarker*dataCont.nBlPar+1:iMarker*dataCont.nBlPar+9]=Closures.turbulentClosure(x[0], x[1], x[4], dataCont.blPar[iMarker*dataCont.nBlPar+0], dataCont.extPar[iMarker*dataCont.nExtPar+0], name)
+ dataCont.blPar[iMarker*dataCont.nBlPar+1:iMarker*dataCont.nBlPar+9]=Closures.turbulentClosure(x[0], x[1],
+ x[4], dataCont.blPar[iMarker*dataCont.nBlPar+0],
+ dataCont.extPar[iMarker*dataCont.nExtPar+0], name)
# Rename variables for clarity
Ta, Ha, _, uea, Cta = x
diff --git a/dart/viscous/model/TSolver.py b/dart/viscous/Solvers/VTimeSolver.py
similarity index 73%
rename from dart/viscous/model/TSolver.py
rename to dart/viscous/Solvers/VTimeSolver.py
index cf78b7fbd004f1465913f9379ed7e45c36f9c4ee..b14541db109ee49089b1773765715bc7be71ebd6 100644
--- a/dart/viscous/model/TSolver.py
+++ b/dart/viscous/Solvers/VTimeSolver.py
@@ -14,11 +14,13 @@
# ------------------------------------------------------------------------------
# Imports
# ------------------------------------------------------------------------------
+from typing import SupportsRound
+from dart.viscous.Physics.VClosures import Closures
+
from matplotlib import pyplot as plt
from fwk.coloring import ccolors
import numpy as np
-import fwk
-from dart.viscous.model.Closures import Closures
+import math
class timeConfig:
""" ---------- Time integration configuration ----------
@@ -33,17 +35,24 @@ class timeConfig:
__cflAdaptParam: numpy vect
[Lower, Upper] bounds of the CFL
"""
- def __init__(self):
+ def __init__(self, _Minfty):
+ if _Minfty != 0:
+ self.__Minf = _Minfty
+ else:
+ self.__Minf = 0.1
+ self.gamma = 1.4
pass
+ def getSoundSpeed(self, gradU2):
+ return np.sqrt(1 / (self.__Minf * self.__Minf) + 0.5 * (self.gamma - 1) * (1 - gradU2))
+
def computeTimeStep(self, CFL, dx, invVel):
+
# Speed of sound.
- freeStreamMach = 0.15
- gamma = 1.4
- gradU2 = invVel**2
- speedOfSound_Squared = 1 / (freeStreamMach * freeStreamMach) + 0.5 * (gamma - 1) * (1 - gradU2)
+ gradU2 = (invVel)**2
# Velocity of the fastest travelling wave
- advVel = np.sqrt(speedOfSound_Squared)+invVel
+ soundSpeed = self.getSoundSpeed(gradU2)
+ advVel = soundSpeed + invVel
# Time step computation
return CFL*dx/advVel
@@ -56,9 +65,9 @@ class timeConfig:
CFL=max(CFL,self.__CFL0)
return CFL
- def outputState(self, var, iMarker, innerIters, normLinSysRes, normLinSysRes0, CFL, color):
+ def outputState(self, var, iMarker, innerIters, normLinSysRes, normLinSysRes0, CFL, color='white'):
if color == 'white':
- print('| {:>6.0f}| {:>11.0f}| {:>9.4f}| {:>8.2f}| {:>14.6f}| {:>6.0f}| {:>11.3f}|'.format(iMarker,
+ print('| {:>6.0f}| {:>11.0f}| {:>9.4f}| {:>10.2f}| {:>14.6f}| {:>6.0f}| {:>11.3f}|'.format(iMarker,
innerIters,
np.log10(normLinSysRes/normLinSysRes0),
CFL,
@@ -108,14 +117,14 @@ class implicitEuler(timeConfig):
safemode: bool
Indicates if the solver is in safe mode (1) or not (0). Value switched according to safeguard
"""
-
- def __init__(self, _solver, transition, wakeBC, _cfl0):
- timeConfig.__init__(self)
+
+ def __init__(self, _solver, transition, wakeBC, _cfl0, _Minfty):
+ timeConfig.__init__(self, _Minfty)
# Initialize time parameters
self.__CFL0=_cfl0
- self.__maxIt = 50000
- self.__tol=1e-6
+ self.__maxIt = 20000
+ self.__tol = 1e-6
self.solver=_solver
self.transSolver = transition
self.wakeBC = wakeBC
@@ -131,6 +140,14 @@ class implicitEuler(timeConfig):
def getmaxIt(self): return self.__maxIt
def setmaxIt(self, _maxIt): self.__maxIt = _maxIt
+ def resetSolution(self, iMarker, var, u0):
+ if iMarker == var.bNodes[1]:
+ var.u[iMarker*var.nVar : (iMarker + 1)*var.nVar] = var.u[0*var.nVar : 1*var.nVar]
+ elif iMarker != var.transMarkers[0] and iMarker != var.transMarkers[1]:
+ var.u[iMarker*var.nVar : (iMarker + 1)*var.nVar] = var.u[(iMarker-1)*var.nVar : (iMarker)*var.nVar]
+ else:
+ var.u[iMarker*var.nVar : (iMarker + 1)*var.nVar] = u0[iMarker*var.nVar : (iMarker + 1)*var.nVar].copy()
+
def timeSolve(self, var):
""" ---------- Newton iteration to solve the boundary layer equations ----------
@@ -212,35 +229,36 @@ class implicitEuler(timeConfig):
continue
if not lockTrans and iMarker-1 < var.bNodes[2] and var.u[(iMarker-1) * var.nVar +2 ] >= var.Ncrit:
- # Transition
- print('| Transition located near x/c = {:<.3f}. Marker: {:<.0f} {:>29s}'.format(var.xGlobal[iMarker-1],iMarker-1, '|'))
+ # Transition
+ print('| Transition located near x/c = {:<2.3f}. Marker: {:<4.0f} {:>28s}'.format(var.xGlobal[iMarker-1],iMarker-1, '|'))
- # Save laminar solution
- lamSol = var.u[(iMarker-1)*var.nVar : (iMarker)* var.nVar].copy()
+ avTurb = self.transSolver.solveTransition(var, iMarker - 1)
+ # Save laminar solution
+ lamSol = var.u[(iMarker-1)*var.nVar : (iMarker)* var.nVar].copy()
- # Set up turbulent portion boundary condition
- avTurb = self.transSolver.solveTransition(var, iMarker - 1)
+ # Set up turbulent portion boundary condition
+ avTurb = self.transSolver.solveTransition(var, iMarker - 1)
- # Compute turbulent solution
- var.flowRegime[iMarker - 1] = 1
- self.newtonSolver(var, iMarker-1, displayState)
+ # Compute turbulent solution
+ var.flowRegime[iMarker - 1] = 1
+ self.newtonSolver(var, iMarker-1, displayState)
- # Average solutions
- var.u[(iMarker-1)*var.nVar : (iMarker)*var.nVar] = (1-avTurb) * lamSol + avTurb * var.u[(iMarker-1)*var.nVar : (iMarker)*var.nVar]
+ # Average solutions
+ var.u[(iMarker-1)*var.nVar : (iMarker)*var.nVar] = (1-avTurb) * lamSol + avTurb * var.u[(iMarker-1)*var.nVar : (iMarker)*var.nVar]
- # Recompute closures
- var.blPar[(iMarker-1)*var.nBlPar+1:(iMarker-1)*var.nBlPar+9]=Closures.turbulentClosure(var.u[(iMarker-1)*var.nVar+0],
- var.u[(iMarker-1)*var.nVar+1],
- var.u[(iMarker-1)*var.nVar+4],
- var.blPar[(iMarker-1)*var.nBlPar+0],
- var.extPar[(iMarker-1)*var.nExtPar+0],
- 'airfoil')
+ # Recompute closures
+ var.blPar[(iMarker-1)*var.nBlPar+1:(iMarker-1)*var.nBlPar+9]=Closures.turbulentClosure(var.u[(iMarker-1)*var.nVar+0],
+ var.u[(iMarker-1)*var.nVar+1],
+ var.u[(iMarker-1)*var.nVar+4],
+ var.blPar[(iMarker-1)*var.nBlPar+0],
+ var.extPar[(iMarker-1)*var.nExtPar+0],
+ 'airfoil')
- lockTrans = 1
+ lockTrans = 1
# Forced transition
- if not lockTrans and var.xtrF[0] is not None and (var.flowRegime[iMarker]== 1 and var.flowRegime[iMarker-1] == 0):
- print('Transition', iMarker)
+ if not lockTrans and var.xtrF[0] is not None and (var.flowRegime[iMarker] == 1 and var.flowRegime[iMarker-1] == 0):
+ print('| Transition near {:<1.4f}. Marker: {:<5.0f} {:>40}'.format(var.xGlobal[iMarker-1],iMarker-1, '|'))
lamSol = var.u[(iMarker-1)*var.nVar : (iMarker)* var.nVar].copy()
avTurb = self.transSolver.transitionBC(var, iMarker - 1)
# Compute turbulent solution
@@ -258,7 +276,7 @@ class implicitEuler(timeConfig):
lockTrans = 1
if not lockTrans and var.xtrF[1] is not None and (var.flowRegime[iMarker]== 1 and var.flowRegime[iMarker-1] == 0):
- print('Transition', iMarker)
+ print('| Transition near {:<1.4f}. Marker: {:<5.0f} {:>40}'.format(var.xGlobal[iMarker-1],iMarker-1, '|'))
lamSol = var.u[(iMarker-1)*var.nVar : (iMarker)* var.nVar].copy()
avTurb = self.transSolver.transitionBC(var, iMarker - 1)
# Compute turbulent solution
@@ -285,13 +303,16 @@ class implicitEuler(timeConfig):
# Save initial condition in case a point diverges.
u0=var.u.copy()
- normLinSysRes0 = np.linalg.norm(self.solver.buildResiduals(var,iMarker))
+ sysRes = self.solver.buildResiduals(var, iMarker)
+ normSysRes0 = np.linalg.norm(sysRes)
# Initial CFL and time step.
CFL = self.__CFL0
CFLAdapt = 1
- jacEval = 5
- dx = var.xx[iMarker] - var.xx[0] if iMarker == var.bNodes[1] else var.xx[iMarker] - var.xx[iMarker - 1]
+ jacEval = 3
+ iMarkerPrev = 0 if iMarker == var.bNodes[1] else iMarker -1
+ dx = var.xx[iMarker] - var.xx[iMarkerPrev]
+
dt = self.computeTimeStep(CFL, dx, var.extPar[iMarker*var.nExtPar + 2])
IMatrix=np.identity(var.nVar) / dt
@@ -300,80 +321,82 @@ class implicitEuler(timeConfig):
innerIters = 0
while innerIters <= self.__maxIt:
- try:
- # Residuals
- linSysRes=self.solver.buildResiduals(var,iMarker)
-
- # Check for convergence
- normLinSysRes=np.linalg.norm(linSysRes)
-
- if displayState and innerIters % 200 == 0:
-
- self.outputState(var, iMarker, innerIters, normLinSysRes, normLinSysRes0, CFL, 'yellow')
+ try:
- if normLinSysRes <= self.__tol*normLinSysRes0:
+ # Compute Residuals.
+ sysRes = self.solver.buildResiduals(var, iMarker)
+ normSysRes = np.linalg.norm(sysRes)
+ if displayState and innerIters % 1000 == 0:
+ self.outputState(var, iMarker, innerIters, normSysRes, normSysRes0, CFL, 'yellow')
+ # Check convergence.
+ if normSysRes <= self.__tol * normSysRes0:
converged = 1
if displayState:
- self.outputState(var, iMarker, innerIters, normLinSysRes, normLinSysRes0, CFL, 'green')
+ self.outputState(var, iMarker, innerIters, normSysRes, normSysRes0, CFL)
break
-
+
+ # Jacobian computation.
+
if innerIters % jacEval == 0:
- Jacobian=self.solver.buildJacobian(var, iMarker, linSysRes = linSysRes, eps = 1e-8)
+ Jacobian=self.solver.buildJacobian(var, iMarker, linSysRes = sysRes, eps = 1e-10)
+
+ # Linear solver.
+
+ slnIncr = np.linalg.solve((IMatrix - Jacobian), sysRes)
+
+ # Increment solution and correct absurd transcient results.
- # Linear solver
- slnIncr = np.linalg.solve((IMatrix - Jacobian), linSysRes)
var.u[iMarker*var.nVar : (iMarker+1)*var.nVar] += slnIncr
var.u[iMarker*var.nVar + 1] = max(var.u[iMarker*var.nVar + 1], 1.0005)
except KeyboardInterrupt:
quit()
except Exception as e:
- print('|', ccolors.ANSI_YELLOW + 'Newton solver failed at position {:<.4f}. Marker: {:<.0f} {:>30s}'.format(var.xGlobal[iMarker], iMarker,'|'), ccolors.ANSI_RESET)
+
+ print('|', ccolors.ANSI_YELLOW + 'Newton solver failed at position {:<.4f}. Marker: {:<.0f} {:>30s}'.format(var.xGlobal[iMarker], iMarker,'|'), ccolors.ANSI_RESET)
print(e)
- print('Current CFL', CFL)
- #var.u[iMarker*var.nVar : (iMarker + 1)*var.nVar] = u0[iMarker*var.nVar : (iMarker + 1)*var.nVar].copy()
- if iMarker == var.bNodes[1]:
- var.u[iMarker*var.nVar : (iMarker + 1)*var.nVar] = var.u[0*var.nVar : 1*var.nVar]
- elif iMarker != var.transMarkers[0] and iMarker != var.transMarkers[1]:
- var.u[iMarker*var.nVar : (iMarker + 1)*var.nVar] = var.u[(iMarker-1)*var.nVar : (iMarker)*var.nVar]
- else:
- var.u[iMarker*var.nVar : (iMarker + 1)*var.nVar] = u0[iMarker*var.nVar : (iMarker + 1)*var.nVar].copy()
+ #self.resetSolution(iMarker, var, u0)
+ var.u[iMarker*var.nVar : (iMarker+1)*var.nVar] = var.u[(iMarker-1)*var.nVar : (iMarker)*var.nVar]
- CFL = min(max(CFL / 2,0.01),1)
+ # Lower CFL and recompute time step
+ print('Current CFL', CFL)
+ if math.isnan(CFL):
+ CFL = self.__CFL0
+ CFL = min(max(CFL / 3,0.01),1)
+ print('Lowering CFL: {:<.3f}'.format(CFL))
+ CFLAdapt = 0
dt = self.computeTimeStep(CFL, dx, var.extPar[iMarker*var.nExtPar + 2])
IMatrix = np.identity(var.nVar) / dt
- CFLAdapt = 0
jacEval = 1
displayState = 1
print('+------------------------------------------------------------------------------+')
- print('| {:>6s}| {:>8s}| {:>9s}| {:>8s}| {:>12s}| {:>6s}| {:>8s}|'.format('Point','Inner Iters', 'rms[F]', 'End CFL', 'Position (x/c)', 'Regime', 'Amp. factor'))
+ print('| {:>6s}| {:>8s}| {:>9s}| {:>8s}| {:>12s}| {:>6s}| {:>8s}|'.format('Point','Inner Iters',
+ 'rms[F]', 'End CFL',
+ 'Position (x/c)', 'Regime',
+ 'Amp. factor'))
print('+------------------------------------------------------------------------------+')
- print('Lowering CFL: {:<.3f}'.format(CFL))
continue
- # Correct absurd transcient
- var.u[iMarker*var.nVar+1] = max(var.u[iMarker*var.nVar+1],1.0005)
-
-
# CFL adaptation
if CFLAdapt:
- if normLinSysRes <= normLinSysRes0:
- CFL = self.__CFL0* (normLinSysRes0/normLinSysRes)**0.7
- elif normLinSysRes > normLinSysRes0:
- CFL = self.__CFL0* (normLinSysRes0/normLinSysRes)**0.5
+
+ CFL = self.__CFL0* (normSysRes0/normSysRes)**0.7
+
CFL = max(CFL, 0.1)
dt = self.computeTimeStep(CFL, dx, var.extPar[iMarker*var.nExtPar + 2])
IMatrix=np.identity(var.nVar) / dt
+
innerIters+=1
else:
# Maximum number of iterations reached.
- if normLinSysRes >= 1e-3 * normLinSysRes0:
- print('|', ccolors.ANSI_RED + 'Too many iteration at position {:<.4f}. normLinSysRes = {:<.3f}. {:>30s}'.format(var.xGlobal[iMarker], np.log10(normLinSysRes/normLinSysRes0),'|'),
+ if normSysRes >= 1e-3 * normSysRes0:
+ print('|', ccolors.ANSI_RED + 'Too many iteration at position {:<.4f}. normLinSysRes = {:<.3f}. {:>30s}'.format(var.xGlobal[iMarker],
+ np.log10(normSysRes/normSysRes0),'|'),
ccolors.ANSI_RESET)
var.u[iMarker*var.nVar:(iMarker+1)*var.nVar] = u0[(iMarker)*var.nVar: (iMarker+1)*var.nVar].copy()
name = 'airfoil' if iMarker <= var.bNodes[2] else 'wake'
@@ -381,15 +404,20 @@ class implicitEuler(timeConfig):
if var.flowRegime[iMarker]==0:
# Laminar closure relations
- var.blPar[iMarker*var.nBlPar+1:iMarker*var.nBlPar+7]=Closures.laminarClosure(var.u[iMarker*var.nVar + 0], var.u[iMarker*var.nVar + 1], var.blPar[iMarker*var.nBlPar+0], var.extPar[iMarker*var.nExtPar+0], name)
+ var.blPar[iMarker*var.nBlPar+1:iMarker*var.nBlPar+7]=Closures.laminarClosure(var.u[iMarker*var.nVar + 0], var.u[iMarker*var.nVar + 1],
+ var.blPar[iMarker*var.nBlPar+0], var.extPar[iMarker*var.nExtPar+0],
+ name)
elif var.flowRegime[iMarker]==1:
# Turbulent closures relations
- var.blPar[iMarker*var.nBlPar+1:iMarker*var.nBlPar+9]=Closures.turbulentClosure(var.u[iMarker*var.nVar + 0], var.u[iMarker*var.nVar + 1], var.u[iMarker*var.nVar + 4], var.blPar[iMarker*var.nBlPar+0], var.extPar[iMarker*var.nExtPar+0], name)
+ var.blPar[iMarker*var.nBlPar+1:iMarker*var.nBlPar+9]=Closures.turbulentClosure(var.u[iMarker*var.nVar + 0], var.u[iMarker*var.nVar + 1],
+ var.u[iMarker*var.nVar + 4], var.blPar[iMarker*var.nBlPar+0],
+ var.extPar[iMarker*var.nExtPar+0], name)
else:
- print('|', ccolors.ANSI_YELLOW+ 'Too many iteration at position {:<.4f}. normLinSysRes = {:<.3f}. {:>30s}'.format(var.xGlobal[iMarker], np.log10(normLinSysRes/normLinSysRes0),'|'),
+ print('|', ccolors.ANSI_YELLOW+ 'Too many iteration at position {:<.4f}. normLinSysRes = {:<.3f}. {:>30s}'.format(var.xGlobal[iMarker],
+ np.log10(normSysRes/normSysRes0),'|'),
ccolors.ANSI_RESET)
return converged
\ No newline at end of file
diff --git a/dart/viscous/model/Transition.py b/dart/viscous/Solvers/VTransition.py
similarity index 82%
rename from dart/viscous/model/Transition.py
rename to dart/viscous/Solvers/VTransition.py
index f82f036f451be075e8e155fdf38b5772a7230f5e..6897bd8169fcaba76d6836c608fe0e926967d1f5 100644
--- a/dart/viscous/model/Transition.py
+++ b/dart/viscous/Solvers/VTransition.py
@@ -1,5 +1,21 @@
+################################################################################
+# #
+# DARTFlo viscous module file #
+# Transition solver : Transition related components. #
+# #
+# Author: Paul Dechamps #
+# Université de Liège #
+# Date: 2021.09.10 #
+# #
+################################################################################
+
+
+# ------------------------------------------------------------------------------
+# Imports
+# ------------------------------------------------------------------------------
+from dart.viscous.Physics.VClosures import Closures
+
import numpy as np
-from dart.viscous.model.Closures import Closures
class Transition:
def __init__(self, xtrF):
diff --git a/dart/viscous/model/vplo.py b/dart/viscous/model/vplo.py
deleted file mode 100755
index d2fdbf5dae1b2799f2379827e8f520cfb4353ae0..0000000000000000000000000000000000000000
--- a/dart/viscous/model/vplo.py
+++ /dev/null
@@ -1,87 +0,0 @@
-################################################################################
-# #
-# Flow viscous module file #
-# vplo: Solution plotter #
-# and comparison with XFoil solution #
-# #
-# Author: Paul Dechamps #
-# Université de Liège #
-# Date: 2021.09.10 #
-# #
-################################################################################
-
-
-# ------------------------------------------------------------------------------
-# Imports
-# ------------------------------------------------------------------------------
-import numpy as np
-import matplotlib.pyplot as plt
-
-# ------------------------------------------------------------------------------
-# Plots the solution and compares it to xfoil file given in 'xfoilFiles' folder
-# ------------------------------------------------------------------------------
-class vplo:
- def __init__(self, case):
- #self.filename=case
- self.filename=case
-
- def showSol(self,pltIn, side, var):
-
-
- with open('/home/parallels/labLnx/waves/flow/viscous/model/xfoilFiles/'+self.filename,'r') as f:
- cl_xfoil=f.readline()
- cd_xfoil=f.readline()
- f.close()
- # columns : x,ue,theta,Cf,H
- data=np.loadtxt('/home/parallels/labLnx/waves/flow/viscous/model/xfoilFiles/'+self.filename,skiprows=3,usecols=(1,3,5,6,7))
-
- # columns : x, theta, H, ue, Cf
- data[:,[1,2,3,4]]=data[:,[2,4,1,3]]
- data[:,]
-
- for i in range(len(data[:,0])):
- if data[i,0]==data[i-1,0]:
- cut=i-1
- np.insert(data,2,0,axis=1)
- showXfoil=0
-
- for i in range(len(pltIn)):
- if pltIn[i]==1:
- if i==3:
- CC=-1
- else:
- CC=1
- if i==0 or i==1:
- var_xfoil=i+1
- else:
- var_xfoil=i
- if side[0]==1:
- if showXfoil==1:
- plt.plot(data[:cut,0],data[:cut,var_xfoil],'r--',lw=0.7)
- plt.plot(var.xGlobal[var.bNodes[0]:var.bNodes[1]],var.u[i:var.bNodes[1]*var.nVar:var.nVar],'r-',lw=1.5)
- plt.xlim([0,1])
- plt.xlabel('x/c')
-
- if side[1]==1:
- if showXfoil==1:
- plt.plot(data[cut:,0],data[cut:,var_xfoil],'g--',lw=0.7)
- plt.plot(np.insert(var.xGlobal[var.bNodes[1]:var.bNodes[2]],0,var.xGlobal[0]),CC*np.insert(var.u[var.bNodes[1]*var.nVar+i:var.bNodes[2]*var.nVar:var.nVar],0,var.u[i]),'g-',lw=1.5)
- plt.xlim([0,1])
- plt.xlabel('x/c')
- if side[2]==1:
- plt.plot(var.xGlobal[var.bNodes[2]:],var.u[var.bNodes[2]*var.nVar+i::var.nVar],'r-',lw=1.5)
- plt.xlim([0,2])
- plt.xlabel('x/c')
- plt.xlim([0,2])
- #plt.axvline(x=var.xtr,ymin=0,'r--',lw=1)
- if i==0:
- plt.title('Theta')
- elif i==1:
- plt.title('H')
- elif i==2:
- plt.title('N')
- elif i==3:
- plt.title('Ue')
- elif i==4:
- plt.title('Ct')
- plt.show()
\ No newline at end of file