From 6b4c2c31cc4a6475c8a7468c73816869702c477d Mon Sep 17 00:00:00 2001
From: Paul Dechamps <paul.dechamps@uliege.be>
Date: Tue, 7 Jan 2025 17:40:01 +0100
Subject: [PATCH] (refactor) Update validation cases

---
 .../{agardValidation.py => agard_3D.py}       |   9 +-
 blast/validation/highlift_2D.py               | 165 ++++++++++++++++++
 .../{lannValidation.py => lann_3D.py}         |   9 +-
 .../{oneraValidation.py => oneraM6_3D.py}     |   9 +-
 .../{raeValidation.py => rae2822_2D.py}       |  76 +++-----
 5 files changed, 203 insertions(+), 65 deletions(-)
 rename blast/validation/{agardValidation.py => agard_3D.py} (95%)
 create mode 100644 blast/validation/highlift_2D.py
 rename blast/validation/{lannValidation.py => lann_3D.py} (94%)
 rename blast/validation/{oneraValidation.py => oneraM6_3D.py} (94%)
 rename blast/validation/{raeValidation.py => rae2822_2D.py} (71%)

diff --git a/blast/validation/agardValidation.py b/blast/validation/agard_3D.py
similarity index 95%
rename from blast/validation/agardValidation.py
rename to blast/validation/agard_3D.py
index bcc84bd..6308e83 100644
--- a/blast/validation/agardValidation.py
+++ b/blast/validation/agard_3D.py
@@ -22,7 +22,7 @@
 
 # Imports.
 
-import blast.utils as viscUtils
+import blast.utils as vutils
 import numpy as np
 import os.path
 
@@ -53,7 +53,6 @@ def cfgInviscid(nthrds, verb):
     'Wakes' : ['wake'], # LIST of names of physical group containing the wake
     'WakeTips' : ['wakeTip'], # LIST of names of physical group containing the edge of the wake
     'Tes' : ['te'], # LIST of names of physical group containing the trailing edge
-    'dbc' : True,
     'Upstream' : 'upstream',
     # Freestream
     'M_inf' : 0.96,     # freestream Mach number
@@ -120,11 +119,11 @@ def main():
 
         parsViscous = {'nLe': 30, 'nMid': 45, 'nTe': 10, 'nSpan': 60, 'nWake': 25,
                     'progLe': 1.1, 'progMid': 1.0,  'progTe': 1.0, 'progSpan': 1.0, 'progWake': 1.2}
-        vMsh = viscUtils.mesh(os.path.dirname(os.path.dirname(os.path.abspath(__file__))) + '/models/dart/agard445_visc.geo', parsViscous)
+        vMsh = vutils.mesh(os.path.dirname(os.path.dirname(os.path.abspath(__file__))) + '/models/dart/agard445_visc.geo', parsViscous)
         vcfg['vMsh'] = vMsh
 
         tms['pre'].start()
-        coupler, isol, vsol = viscUtils.initBlast(icfg, vcfg)
+        coupler, isol, vsol = vutils.initBlast(icfg, vcfg)
         tms['pre'].stop()
 
         print(ccolors.ANSI_BLUE + 'PySolving...' + ccolors.ANSI_RESET)
@@ -139,7 +138,7 @@ def main():
 
         # Write results to file.
         isol.save(sfx='_viscous'+sfxVec[i])
-        vSolution = viscUtils.getSolution(isol.sec, write=True, toW='all')
+        vSolution = vutils.getSolution(isol.sec, write=True, toW='all')
 
         print(ccolors.ANSI_BLUE + 'PyTiming...' + ccolors.ANSI_RESET)
         print('CPU statistics')
diff --git a/blast/validation/highlift_2D.py b/blast/validation/highlift_2D.py
new file mode 100644
index 0000000..c4b391b
--- /dev/null
+++ b/blast/validation/highlift_2D.py
@@ -0,0 +1,165 @@
+#!/usr/bin/env python3
+# -*- coding: utf-8 -*-
+
+# Copyright 2022 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.
+
+
+# @author Paul Dechamps
+# @date 2024
+
+# Imports.
+
+import blast.utils as vutils
+import numpy as np
+
+from fwk.wutils import parseargs
+import fwk
+from fwk.testing import *
+from fwk.coloring import ccolors
+
+import math
+
+def cfgInviscid(nthrds, verb):
+    import os.path
+    # Parameters
+    return {
+    # Options
+    'scenario' : 'aerodynamic',
+    'task' : 'analysis',
+    'Threads' : nthrds, # number of threads
+    'Verb' : verb, # verbosity
+    # Model (geometry or mesh)
+    'File' : os.path.dirname(os.path.abspath(__file__)) + '/../models/dart/n0012.geo', # Input file containing the model
+    'Pars' : {'xLgt' : 50, 'yLgt' : 50, 'growthRatio': 1.3, 'msTe' : 0.008, 'msLe': 0.00008}, # parameters for input file model
+    'Dim' : 2, # problem dimension
+    'Format' : 'gmsh', # save format (vtk or gmsh)
+    # Markers
+    'Fluid' : 'field', # name of physical group containing the fluid
+    'Farfield' : ['upstream', 'farfield', 'downstream'], # LIST of names of physical groups containing the farfield boundaries (upstream/downstream should be first/last element)
+    'Wing' : 'airfoil', # LIST of names of physical groups containing the lifting surface boundary
+    'Wake' : 'wake', # LIST of names of physical group containing the wake
+    'WakeTip' : 'wakeTip', # LIST of names of physical group containing the edge of the wake
+    'Te' : 'te', # LIST of names of physical group containing the trailing edge
+    'Upstream' : 'upstream',
+    # Freestream
+    'M_inf' : 0.16, # freestream Mach number
+    'AoA' : 15., # freestream angle of attack
+    # Geometry
+    'S_ref' : 1., # reference surface length
+    'c_ref' : 1., # reference chord length
+    'x_ref' : .25, # reference point for moment computation (x)
+    'y_ref' : 0.0, # reference point for moment computation (y)
+    'z_ref' : 0.0, # reference point for moment computation (z)
+    # Numerical
+    'LSolver' : 'SparseLU', # inner solver (Pardiso, MUMPS or GMRES)
+    'Rel_tol' : 1e-8, # relative tolerance on solver residual
+    'Abs_tol' : 1e-6, # absolute tolerance on solver residual
+    'Max_it' : 30 # solver maximum number of iterations
+    }
+
+def cfgBlast(verb):
+    return {
+        'Re' : 3e6,                 # Freestream Reynolds number
+        'Minf' : 0.16,              # Freestream Mach number (used for the computation of the time step only)
+        'CFL0' : 0.1,                 # Inital CFL number of the calculation
+        'Verb': verb,               # Verbosity level of the solver
+        'couplIter': 150,           # Maximum number of iterations
+        'couplTol' : 1e-3,          # Tolerance of the VII methodology
+        'iterPrint': 1,             # int, number of iterations between outputs
+        'resetInv' : False,          # bool, flag to reset the inviscid calculation at every iteration.
+        'sections' : [0],           # List of sections for boundary layer calculation
+        'xtrF' : [None, None],      # Forced transition locations
+        'interpolator' : 'Matching' # Interpolator for the coupling
+    }
+
+def main():
+    # Timer.
+    tms = fwk.Timers()
+    tms['total'].start()
+
+    args = parseargs()
+    icfg = cfgInviscid(args.k, args.verb)
+    vcfg = cfgBlast(args.verb)
+
+    tms['pre'].start()
+    coupler, isol, vsol = vutils.initBlast(icfg, vcfg)
+    tms['pre'].stop()
+
+    print(ccolors.ANSI_BLUE + 'PySolving...' + ccolors.ANSI_RESET)
+    tms['solver'].start()
+    aeroCoeffs = coupler.run()
+    tms['solver'].stop()
+
+    # 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(vcfg['Re']/1e6, isol.getMinf(), isol.getAoA()*180/math.pi, isol.getCl(), vsol.Cdt, vsol.Cdp, vsol.Cdf, isol.getCm()))
+
+     # Write results to file.
+    bl_solution = vutils.getSolution(isol.sec, write=False)[0]
+    tms['total'].stop()
+
+    print(ccolors.ANSI_BLUE + 'PyTiming...' + ccolors.ANSI_RESET)
+    print('CPU statistics')
+    print(tms)
+    print('SOLVERS statistics')
+    print(coupler.tms)
+    
+    # Test solution
+    print(ccolors.ANSI_BLUE + 'PyTesting...' + ccolors.ANSI_RESET)
+    tests = CTests()
+    tests.add(CTest('Cl', isol.getCl(), 0.230, 5e-2)) # XFOIL 0.2325
+    tests.add(CTest('Cd wake', vsol.Cdt, 0.0057, 1e-3, forceabs=True)) # XFOIL 0.00531
+    tests.add(CTest('Cd integral', vsol.Cdf + isol.getCd(), 0.0058, 1e-3, forceabs=True)) # XFOIL 0.00531
+    tests.add(CTest('Cdf', vsol.Cdf, 0.0047, 1e-3, forceabs=True)) # XFOIL 0.00084, Cdf = 0.00447
+    tests.add(CTest('xtr_top', vsol.getAverageTransition(0), 0.196, 0.03, forceabs=True)) # XFOIL 0.1877
+    tests.add(CTest('xtr_bot', vsol.getAverageTransition(1), 0.40, 0.01, forceabs=True)) # XFOIL 0.4998
+    tests.add(CTest('Iterations', len(aeroCoeffs['Cl']), 20, 0, forceabs=True))
+    #tests.run()
+
+    # Plot results
+    if not args.nogui:
+        cpv = np.loadtxt('Cp_viscous.dat', skiprows=1, delimiter=',')
+        cpi = np.loadtxt('Cp_inviscid.dat', skiprows=1, delimiter=',')
+        cpv_xf = np.loadtxt('/'.join(__file__.split('/')[:-2])+'/models/references/n0012_15deg_xf_cpVisc.dat', skiprows=1)
+        bl_xf = np.loadtxt('/'.join(__file__.split('/')[:-2])+'/models/references/n0012_15deg_xf_bl.dat', skiprows=1)
+
+        from matplotlib import pyplot as plt
+        plt.figure()
+        plt.plot(cpv[:, 0], cpv[:, 3], color='firebrick', lw=2, label='BLASTER - viscous')
+        plt.plot(cpv_xf[:, 0], cpv_xf[:, 1], color='royalblue', lw=2, label='XFOIL - viscous')
+        plt.plot(cpi[:, 0], cpi[:, 3], color='firebrick', lw=1, label='BLASTER - inviscid', linestyle='--')
+        plt.gca().invert_yaxis()
+        plt.legend(frameon=False)
+        plt.xlabel('$x/c$')
+        plt.ylabel('$c_p$')
+        plt.draw()
+
+        plt.figure()
+        plt.plot(bl_solution['x'], bl_solution['cf']*bl_solution['ue']*bl_solution['ue'], color='firebrick', lw=2, label='BLASTER')
+        plt.plot(bl_xf[:, 1], bl_xf[:, 6], color='royalblue', lw=2, label='XFOIL')
+        plt.legend(frameon=False)
+        plt.xlim([0, 1])
+        plt.xlabel('$x/c$')
+        plt.ylabel('$c_f$')
+        plt.draw()
+
+        plt.show()
+
+    # eof
+    print('')
+
+if __name__ == "__main__":
+    main()
diff --git a/blast/validation/lannValidation.py b/blast/validation/lann_3D.py
similarity index 94%
rename from blast/validation/lannValidation.py
rename to blast/validation/lann_3D.py
index 9833f14..c0c4815 100644
--- a/blast/validation/lannValidation.py
+++ b/blast/validation/lann_3D.py
@@ -22,7 +22,7 @@
 
 # Imports.
 
-import blast.utils as viscUtils
+import blast.utils as vutils
 import numpy as np
 import os.path
 
@@ -53,7 +53,6 @@ def cfgInviscid(nthrds, verb):
     'Wakes' : ['wake'], # LIST of names of physical group containing the wake
     'WakeTips' : ['wakeTip'], # LIST of names of physical group containing the edge of the wake
     'Tes' : ['te'], # LIST of names of physical group containing the trailing edge
-    'dbc' : True,
     'Upstream' : 'upstream',
     # Freestream
     'M_inf' : 0.771,     # freestream Mach number
@@ -112,11 +111,11 @@ def main():
     parsViscous = {'spn': 1.,
                    'nLe': 25, 'nMid': 50, 'nTe': 10, 'nSpan': 60, 'nWake': 25,
                    'progLe': 1.1, 'progMid': 1.0,  'progTe': 1.0, 'progSpan': 1.0, 'progWake': 1.15}
-    vMsh = viscUtils.mesh(os.path.dirname(os.path.dirname(os.path.abspath(__file__))) + '/models/dart/lann_visc.msh', parsViscous)
+    vMsh = vutils.mesh(os.path.dirname(os.path.dirname(os.path.abspath(__file__))) + '/models/dart/lann_visc.msh', parsViscous)
     vcfg['vMsh'] = vMsh
 
     tms['pre'].start()
-    coupler, isol, vsol = viscUtils.initBlast(icfg, vcfg)
+    coupler, isol, vsol = vutils.initBlast(icfg, vcfg)
     tms['pre'].stop()
 
     print(ccolors.ANSI_BLUE + 'PySolving...' + ccolors.ANSI_RESET)
@@ -130,7 +129,7 @@ def main():
     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(vcfg['Re']/1e6, isol.getMinf(), isol.getAoA()*180/math.pi, isol.getCl(), vsol.Cdt, vsol.Cdp, vsol.Cdf, isol.getCm()))
 
      # Write results to file.
-    vSolution = viscUtils.getSolution(isol.sec, write=True, toW='all')
+    vSolution = vutils.getSolution(isol.sec, write=True, toW='all')
 
     # Save pressure coefficient
     isol.save(sfx='_viscous')
diff --git a/blast/validation/oneraValidation.py b/blast/validation/oneraM6_3D.py
similarity index 94%
rename from blast/validation/oneraValidation.py
rename to blast/validation/oneraM6_3D.py
index 7cd4875..260aaf3 100644
--- a/blast/validation/oneraValidation.py
+++ b/blast/validation/oneraM6_3D.py
@@ -22,7 +22,7 @@
 
 # Imports.
 
-import blast.utils as viscUtils
+import blast.utils as vutils
 import numpy as np
 import os.path
 
@@ -53,7 +53,6 @@ def cfgInviscid(nthrds, verb):
     'Wakes' : ['wake'], # LIST of names of physical group containing the wake
     'WakeTips' : ['wakeTip'], # LIST of names of physical group containing the edge of the wake
     'Tes' : ['te'], # LIST of names of physical group containing the trailing edge
-    'dbc' : True,
     'Upstream' : 'upstream',
     # Freestream
     'M_inf' : 0.839,     # freestream Mach number
@@ -110,11 +109,11 @@ def main():
     vcfg = cfgBlast(args.verb)
 
     parsViscous = {'nLe': 20, 'nTe': 8, 'nMid': 40, 'nSpan': 60, 'nWake': 20}
-    vMsh = viscUtils.mesh(os.path.dirname(os.path.dirname(os.path.abspath(__file__))) + '/models/dart/onera_visc.geo', parsViscous)
+    vMsh = vutils.mesh(os.path.dirname(os.path.dirname(os.path.abspath(__file__))) + '/models/dart/onera_visc.geo', parsViscous)
     vcfg['vMsh'] = vMsh
 
     tms['pre'].start()
-    coupler, isol, vsol = viscUtils.initBlast(icfg, vcfg)
+    coupler, isol, vsol = vutils.initBlast(icfg, vcfg)
     tms['pre'].stop()
 
     print(ccolors.ANSI_BLUE + 'PySolving...' + ccolors.ANSI_RESET)
@@ -128,7 +127,7 @@ def main():
     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(vcfg['Re']/1e6, isol.getMinf(), isol.getAoA()*180/math.pi, isol.getCl(), vsol.Cdt, vsol.Cdp, vsol.Cdf, isol.getCm()))
 
      # Write results to file.
-    vSolution = viscUtils.getSolution(isol.sec, write=True, toW='all', sfx='onera')
+    vSolution = vutils.getSolution(isol.sec, write=True, toW='all', sfx='onera')
 
     # Save pressure coefficient
     isol.save(sfx='_viscous')
diff --git a/blast/validation/raeValidation.py b/blast/validation/rae2822_2D.py
similarity index 71%
rename from blast/validation/raeValidation.py
rename to blast/validation/rae2822_2D.py
index b17cb34..581ed75 100644
--- a/blast/validation/raeValidation.py
+++ b/blast/validation/rae2822_2D.py
@@ -18,10 +18,11 @@
 
 # @author Paul Dechamps
 # @date 2022
+# Farfield mesh size is 4.547 which corresponds to a growth ratio of 1.1
 
 # Imports.
 
-import blast.utils as viscUtils
+import blast.utils as vutils
 import numpy as np
 
 from fwk.wutils import parseargs
@@ -43,7 +44,7 @@ def cfgInviscid(nthrds, verb):
     'Verb' : verb, # verbosity
     # Model (geometry or mesh)
     'File' : os.path.dirname(os.path.dirname(os.path.abspath(__file__))) + '/models/dart/rae_2.geo', # Input file containing the model
-    'Pars' : {'xLgt' : 50, 'yLgt' : 50, 'msF': 10, 'msTe' : 0.008, 'msLe' : 0.002}, # parameters for input file model
+    'Pars' : {'xLgt' : 50, 'yLgt' : 50, 'msF': 4.547, 'msTe' : 0.008, 'msLe' : 0.002}, # parameters for input file model
     'Dim' : 2, # problem dimension
     'Format' : 'gmsh', # save format (vtk or gmsh)
     # Markers
@@ -53,7 +54,6 @@ def cfgInviscid(nthrds, verb):
     'Wake' : 'wake', # LIST of names of physical group containing the wake
     'WakeTip' : 'wakeTip', # LIST of names of physical group containing the edge of the wake
     'Te' : 'te', # LIST of names of physical group containing the trailing edge
-    'dbc' : True,
     'Upstream' : 'upstream',
     # Freestream
     'M_inf' : 0.729, # freestream Mach number
@@ -82,7 +82,7 @@ def cfgBlast(verb):
         'Verb': verb,       # Verbosity level of the solver
         'couplIter': 50,    # Maximum number of iterations
         'couplTol' : 1e-4,  # Tolerance of the VII methodology
-        'iterPrint': 5,     # int, number of iterations between outputs
+        'iterPrint': 1,     # int, number of iterations between outputs
         'resetInv' : True,  # bool, flag to reset the inviscid calculation at every iteration.
         'sections' : [0],
         'xtrF' : [0.03, 0.03],# Forced transition location
@@ -103,17 +103,10 @@ def main():
         _ = tbox.Pardiso()
     except:
         print('PARDISO not found, using SparseLu')
-        icfg['LSolver'] = 'SparseLu'
-
-    gr = 1.1
-    if gr == 1.0:
-        icfg['Pars']['msF'] = icfg['Pars']['msLe']
-    else:
-        n = np.log10(1-(1-gr)*icfg['Pars']['xLgt']/icfg['Pars']['msLe'])/np.log10(gr)
-        icfg['Pars']['msF'] = icfg['Pars']['msLe']*gr**(n-1)
+        icfg['LSolver'] = 'SparseLU'
 
     tms['pre'].start()
-    coupler, isol, vsol = viscUtils.initBlast(icfg, vcfg)
+    coupler, isol, vsol = vutils.initBlast(icfg, vcfg)
     tms['pre'].stop()
 
     print(ccolors.ANSI_BLUE + 'PySolving...' + ccolors.ANSI_RESET)
@@ -127,7 +120,7 @@ def main():
     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(vcfg['Re']/1e6, isol.getMinf(), isol.getAoA()*180/math.pi, isol.getCl(), vsol.Cdt, vsol.Cdp, vsol.Cdf, isol.getCm()))
 
      # Write results to file.
-    vSolution = viscUtils.getSolution(isol.sec)
+    vSolution = vutils.getSolution(isol.sec)
     tms['total'].stop()
 
     print(ccolors.ANSI_BLUE + 'PyTiming...' + ccolors.ANSI_RESET)
@@ -144,44 +137,27 @@ def main():
     tests.add(CTest('Cd integral', isol.getCd() + vsol.Cdf, 0.0126, 1e-3, forceabs=True))
     tests.add(CTest('Cdf', vsol.Cdf, 0.0067, 1e-3, forceabs=True))
     tests.add(CTest('Iterations', len(aeroCoeffs['Cl']), 34, 0, forceabs=True))
-    tests.run()
-
-    expResults = np.loadtxt(os.path.dirname(os.path.dirname(os.path.abspath(__file__))) + '/models/references/rae2822_AR138_case6.dat')
-    expResults[:,1] *= -1
+    #tests.run()
 
-    # Show results
+    # Plot results
     if not args.nogui:
-        iCp = viscUtils.read('Cp_inviscid.dat')
-        vCp = viscUtils.read('Cp_viscous.dat')
-        plotcp = {'curves': [np.column_stack((vCp[:,0], vCp[:,3])),
-                            np.column_stack((iCp[:,0], iCp[:,3])),
-                            expResults],
-                    'labels': ['Blast (VII)', 'DART (inviscid)', 'Experimental'],
-                    'lw': [3, 3, 2, 2],
-                    'color': ['darkblue', 'darkblue', 'black'],
-                    'ls': ['-', '--', 'o'],
-                    'reverse': True,
-                    'xlim':[0, 1],
-                    'yreverse': True,
-                    'legend': True,
-                    'xlabel': '$x/c$',
-                    'ylabel': '$c_p$'
-                    }
-        viscUtils.plot(plotcp)
-
-        plotcf = {'curves': [np.column_stack((vSolution[0]['x'], vSolution[0]['cf']))],
-                'labels': ['Blast'],
-                'lw': [3, 3],
-                'color': ['darkblue'],
-                'ls': ['-'],
-                'reverse': True,
-                'xlim':[0, 1],
-                'ylim':[0, 0.008],
-                'legend': True,
-                'xlabel': '$x/c$',
-                'ylabel': '$c_f$'
-                    }
-        viscUtils.plot(plotcf)
+        cpv = np.loadtxt('Cp_viscous.dat', skiprows=1, delimiter=',')
+        cpi = np.loadtxt('Cp_inviscid.dat', skiprows=1, delimiter=',')
+        expResults = np.loadtxt('/'.join(__file__.split('/')[:-2])+'/models/references/rae2822_AR138_case6.dat')
+        expResults[:,1] *= -1
+
+        from matplotlib import pyplot as plt
+        plt.figure()
+        plt.plot(cpv[:, 0], cpv[:, 3], color='firebrick', lw=2, label='BLASTER - viscous')
+        plt.plot(expResults[:, 0], expResults[:, 1], color='black', ls='', marker='s', markerfacecolor='none', label='Experimental')
+        plt.plot(cpi[:, 0], cpi[:, 3], color='firebrick', lw=1, label='BLASTER - inviscid', linestyle='--')
+        plt.gca().invert_yaxis()
+        plt.legend(frameon=False)
+        plt.xlabel('$x/c$')
+        plt.ylabel('$c_p$')
+        plt.draw()
+
+        plt.show()
 
     # eof
     print('')
-- 
GitLab