diff --git a/blast/validation/raeValidation.py b/blast/validation/raeValidation.py
new file mode 100644
index 0000000000000000000000000000000000000000..9a99c222fbda711fb0e6d53f8cf5b33bb0c0f426
--- /dev/null
+++ b/blast/validation/raeValidation.py
@@ -0,0 +1,172 @@
+#!/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 2022
+
+# Imports.
+
+import blast.utils as viscUtils
+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.dirname(os.path.abspath(__file__))) + '/models/dart/rae_2.geo', # Input file containing the model
+ 'Pars' : {'xLgt' : 5, 'yLgt' : 5, 'msF' : 1, 'msTe' : 0.01, 'msLe' : 0.002}, # 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
+ 'dbc' : True,
+ 'Upstream' : 'upstream',
+ # Freestream
+ 'M_inf' : 0.729, # freestream Mach number
+ 'AoA' : 2.31, # 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' : 'GMRES', # inner solver (Pardiso, MUMPS or GMRES)
+ 'G_fill' : 2, # fill-in factor for GMRES preconditioner
+ 'G_tol' : 1e-5, # tolerance for GMRES
+ 'G_restart' : 50, # restart for GMRES
+ 'Rel_tol' : 1e-6, # relative tolerance on solver residual
+ 'Abs_tol' : 1e-8, # absolute tolerance on solver residual
+ 'Max_it' : 50 # solver maximum number of iterations
+ }
+
+def cfgBlast(verb):
+ return {
+ 'Re' : 6.5e6, # Freestream Reynolds number
+ 'Minf' : 0.73, # Freestream Mach number (used for the computation of the time step only)
+ 'CFL0' : 1, # Inital CFL number of the calculation
+ '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
+ 'resetInv' : True, # bool, flag to reset the inviscid calculation at every iteration.
+ 'sections' : [0],
+ 'xtrF' : [0., 0.],# Forced transition location
+ 'interpolator' : 'Matching',
+ 'nDim' : 2
+ }
+
+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, iSolverAPI, vSolver = viscUtils.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, iSolverAPI.getMinf(), iSolverAPI.getAoA()*180/math.pi, iSolverAPI.getCl(), vSolver.Cdt, vSolver.Cdp, vSolver.Cdf, iSolverAPI.getCm()))
+
+ # Write results to file.
+ vSolution = viscUtils.getSolution(vSolver)
+ vSolution['Cdt_int'] = vSolution['Cdf'] + iSolverAPI.getCd()
+ 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', iSolverAPI.getCl(), 0.752, 5e-2)) # XFOIL 0.2325
+ tests.add(CTest('Cd wake', vSolution['Cdt_w'], 0.0094, 1e-3, forceabs=True)) # XFOIL 0.00531
+ tests.add(CTest('Cd integral', vSolution['Cdt_int'], 0.0145, 1e-3, forceabs=True)) # XFOIL 0.00531
+ tests.add(CTest('Cdf', vSolution['Cdf'], 0.0069, 1e-3, forceabs=True)) # XFOIL 0.00084, Cdf = 0.00447
+ tests.add(CTest('Iterations', len(aeroCoeffs), 35, 0, forceabs=True))
+ tests.run()
+
+ # Show 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]))],
+ 'labels': ['Blast (VII)', 'DART (inviscid)'],
+ 'lw': [3, 3, 2, 2],
+ 'color': ['darkblue', 'darkblue'],
+ 'ls': ['-', '--'],
+ 'reverse': True,
+ 'xlim':[0, 1],
+ 'yreverse': True,
+ 'legend': True,
+ 'xlabel': '$x/c$',
+ 'ylabel': '$c_p$'
+ }
+ viscUtils.plot(plotcp)
+
+ plotcf = {'curves': [np.column_stack((vSolution['x'], vSolution['Cf']))],
+ 'labels': ['Blast'],
+ 'lw': [3, 3],
+ 'color': ['darkblue'],
+ 'ls': ['-'],
+ 'reverse': True,
+ 'xlim':[0, 1],
+ 'legend': True,
+ 'xlabel': '$x/c$',
+ 'ylabel': '$c_f$'
+ }
+ viscUtils.plot(plotcf)
+
+ # eof
+ print('')
+
+if __name__ == "__main__":
+ main()