diff --git a/blast/coupler.py b/blast/coupler.py
index 397f8434e6af9d579951116eedca27f05975e951..b94e05e2079d2005595ffd889c2c30f37ffcf569 100644
--- a/blast/coupler.py
+++ b/blast/coupler.py
@@ -87,7 +87,7 @@ class Coupler:
print('{:>5s}| {:>7s} {:>7s} {:>7s} | {:>6s} {:>6s} | {:>6s}'.format('It', 'Cl', 'Cd', 'Cdwake', 'xtrT', 'xtrB', 'Error'))
print(' ----------------------------------------------------------')
if couplIter % self.iterPrint == 0:
- print(' {:>4.0f}| {:>7.4f} {:>7.4f} {:>7.4f} | {:>6.4f} {:>6.4f} | {:>6.3f}'.format(couplIter, self.iSolverAPI.getCl(), self.iSolverAPI.getCd()+self.vSolver.Cdf, self.vSolver.Cdt, self.vSolver.getxoctr(0, 0), self.vSolver.getxoctr(0, 1), np.log10(error)))
+ print(' {:>4.0f}| {:>7.4f} {:>7.4f} {:>7.4f} | {:>6.4f} {:>6.4f} | {:>6.3f}'.format(couplIter, self.iSolverAPI.getCl(), self.iSolverAPI.getCd()+self.vSolver.Cdf, self.vSolver.Cdt, self.vSolver.getAvrgxoctr(0), self.vSolver.getAvrgxoctr(1), np.log10(error)))
if self.iSolverAPI.getVerbose() != 0 or self.vSolver.verbose != 0:
print('')
couplIter += 1
@@ -97,6 +97,6 @@ class Coupler:
print('')
print('{:>5s}| {:>7s} {:>7s} {:>7s} | {:>6s} {:>8s} | {:>6s}'.format('It', 'Cl', 'Cd', 'Cdwake', 'xtrT', 'xtrB', 'Error'))
print(' ----------------------------------------------------------')
- print(ccolors.ANSI_RED, '{:>4.0f}| {:>7.5f} {:>7.5f} {:>7.5f} | {:>6.4f} {:>7.4f} | {:>6.3f}\n'.format(couplIter, self.iSolverAPI.getCl(), self.iSolverAPI.getCd()+self.vSolver.Cdf, self.vSolver.Cdt, self.vSolver.getxoctr(0, 0), self.vSolver.getxoctr(0, 1), np.log10(error)), ccolors.ANSI_RESET)
+ print(ccolors.ANSI_RED, '{:>4.0f}| {:>7.5f} {:>7.5f} {:>7.5f} | {:>6.4f} {:>7.4f} | {:>6.3f}\n'.format(couplIter, self.iSolverAPI.getCl(), self.iSolverAPI.getCd()+self.vSolver.Cdf, self.vSolver.Cdt, self.vSolver.getAvrgxoctr(0), self.vSolver.getAvrgxoctr(1), np.log10(error)), ccolors.ANSI_RESET)
self.iSolverAPI.writeCp(sfx='_viscous'+self.filesfx)
return aeroCoeffs
diff --git a/blast/models/dart/agard445_visc.geo b/blast/models/dart/agard445_visc.geo
new file mode 100644
index 0000000000000000000000000000000000000000..c657f6c6d731e5268a415c04c5bf3853bb057505
--- /dev/null
+++ b/blast/models/dart/agard445_visc.geo
@@ -0,0 +1,229 @@
+/* AGARD 445 wing */
+// Initially generated by unsRgridWingGen.m
+
+// Parameters
+// domain and mesh
+DefineConstant[ nSpan = { 20, Name "Spanwise number of elements" } ];
+DefineConstant[ nLe = { 5, Name "LE chordwise number of elements" } ];
+DefineConstant[ nTe = { 5, Name "TE chordwise number of elements" } ];
+DefineConstant[ nMid = { 10, Name "Mid chordwise number of elements" } ];
+DefineConstant[ nWake = { 15, Name "Wake number of elements" } ];
+DefineConstant[ progSpan = { 1, Name "Spanwise progression" } ];
+DefineConstant[ progLe = { 1, Name "LE progression" } ];
+DefineConstant[ progTe = { 1, Name "TE progression" } ];
+DefineConstant[ progMid = { 1, Name "Mid progression" } ];
+DefineConstant[ progWake = { 1.0, Name "Wake progression" } ];
+DefineConstant[ xO = { -3, Name "Box origin (x)" } ];
+DefineConstant[ zO = { -3, Name "Box origin (z)" } ];
+DefineConstant[ xL = { 7, Name "Box length" } ];
+DefineConstant[ yL = { 3, Name "Box width" } ];
+DefineConstant[ zL = { 6, Name "Box height" } ];
+
+//// GEOMETRY
+
+
+/// Points
+
+// Airfoil 1: agard445, 51 points
+Point(1) = {0.559000,0.000000,0.000000};
+Point(2) = {0.531050,0.000000,0.001398};
+Point(3) = {0.503100,0.000000,0.002739};
+Point(4) = {0.475150,0.000000,0.004075};
+Point(5) = {0.447200,0.000000,0.005406};
+Point(6) = {0.419250,0.000000,0.006680};
+Point(7) = {0.391300,0.000000,0.007837};
+Point(8) = {0.363350,0.000000,0.008866};
+Point(9) = {0.335400,0.000000,0.009743};
+Point(10) = {0.307450,0.000000,0.010442};
+Point(11) = {0.279500,0.000000,0.010923};
+Point(12) = {0.251550,0.000000,0.011158};
+Point(13) = {0.223600,0.000000,0.011163};
+Point(14) = {0.195650,0.000000,0.010968};
+Point(15) = {0.167700,0.000000,0.010576};
+Point(16) = {0.139750,0.000000,0.009995};
+Point(17) = {0.111800,0.000000,0.009196};
+Point(18) = {0.083850,0.000000,0.008156};
+Point(19) = {0.055900,0.000000,0.006781};
+Point(20) = {0.041925,0.000000,0.005920};
+Point(21) = {0.027950,0.000000,0.004891};
+Point(22) = {0.013975,0.000000,0.003617};
+Point(23) = {0.006988,0.000000,0.002622};
+Point(24) = {0.004193,0.000000,0.002057};
+Point(25) = {0.002795,0.000000,0.001699};
+Point(26) = {0.000000,0.000000,0.000000};
+Point(27) = {0.002795,0.000000,-0.001699};
+Point(28) = {0.004193,0.000000,-0.002057};
+Point(29) = {0.006988,0.000000,-0.002622};
+Point(30) = {0.013975,0.000000,-0.003617};
+Point(31) = {0.027950,0.000000,-0.004891};
+Point(32) = {0.041925,0.000000,-0.005920};
+Point(33) = {0.055900,0.000000,-0.006781};
+Point(34) = {0.083850,0.000000,-0.008156};
+Point(35) = {0.111800,0.000000,-0.009196};
+Point(36) = {0.139750,0.000000,-0.009995};
+Point(37) = {0.167700,0.000000,-0.010576};
+Point(38) = {0.195650,0.000000,-0.010968};
+Point(39) = {0.223600,0.000000,-0.011163};
+Point(40) = {0.251550,0.000000,-0.011158};
+Point(41) = {0.279500,0.000000,-0.010923};
+Point(42) = {0.307450,0.000000,-0.010442};
+Point(43) = {0.335400,0.000000,-0.009743};
+Point(44) = {0.363350,0.000000,-0.008866};
+Point(45) = {0.391300,0.000000,-0.007837};
+Point(46) = {0.419250,0.000000,-0.006680};
+Point(47) = {0.447200,0.000000,-0.005406};
+Point(48) = {0.475150,0.000000,-0.004075};
+Point(49) = {0.503100,0.000000,-0.002739};
+Point(50) = {0.531050,0.000000,-0.001398};
+
+// Airfoil 2: agard445, 51 points
+Point(52) = {1.178128,0.762000,0.000000};
+Point(53) = {1.159709,0.762000,0.000921};
+Point(54) = {1.141290,0.762000,0.001805};
+Point(55) = {1.122870,0.762000,0.002685};
+Point(56) = {1.104451,0.762000,0.003562};
+Point(57) = {1.086032,0.762000,0.004402};
+Point(58) = {1.067613,0.762000,0.005165};
+Point(59) = {1.049194,0.762000,0.005843};
+Point(60) = {1.030775,0.762000,0.006421};
+Point(61) = {1.012356,0.762000,0.006881};
+Point(62) = {0.993937,0.762000,0.007198};
+Point(63) = {0.975518,0.762000,0.007353};
+Point(64) = {0.957099,0.762000,0.007357};
+Point(65) = {0.938680,0.762000,0.007228};
+Point(66) = {0.920261,0.762000,0.006970};
+Point(67) = {0.901842,0.762000,0.006587};
+Point(68) = {0.883423,0.762000,0.006060};
+Point(69) = {0.865004,0.762000,0.005375};
+Point(70) = {0.846585,0.762000,0.004468};
+Point(71) = {0.837375,0.762000,0.003901};
+Point(72) = {0.828166,0.762000,0.003223};
+Point(73) = {0.818956,0.762000,0.002383};
+Point(74) = {0.814351,0.762000,0.001728};
+Point(75) = {0.812509,0.762000,0.001356};
+Point(76) = {0.811589,0.762000,0.001120};
+Point(77) = {0.809747,0.762000,0.000000};
+Point(78) = {0.811589,0.762000,-0.001120};
+Point(79) = {0.812509,0.762000,-0.001356};
+Point(80) = {0.814351,0.762000,-0.001728};
+Point(81) = {0.818956,0.762000,-0.002383};
+Point(82) = {0.828166,0.762000,-0.003223};
+Point(83) = {0.837375,0.762000,-0.003901};
+Point(84) = {0.846585,0.762000,-0.004468};
+Point(85) = {0.865004,0.762000,-0.005375};
+Point(86) = {0.883423,0.762000,-0.006060};
+Point(87) = {0.901842,0.762000,-0.006587};
+Point(88) = {0.920261,0.762000,-0.006970};
+Point(89) = {0.938680,0.762000,-0.007228};
+Point(90) = {0.957099,0.762000,-0.007357};
+Point(91) = {0.975518,0.762000,-0.007353};
+Point(92) = {0.993937,0.762000,-0.007198};
+Point(93) = {1.012356,0.762000,-0.006881};
+Point(94) = {1.030775,0.762000,-0.006421};
+Point(95) = {1.049194,0.762000,-0.005843};
+Point(96) = {1.067613,0.762000,-0.005165};
+Point(97) = {1.086032,0.762000,-0.004402};
+Point(98) = {1.104451,0.762000,-0.003562};
+Point(99) = {1.122870,0.762000,-0.002685};
+Point(100) = {1.141290,0.762000,-0.001805};
+Point(101) = {1.159709,0.762000,-0.000921};
+
+
+// Midplane:
+Point(5351) = {xO+xL,0.000000,0.000000};
+Point(5352) = {xO+xL,0.762000,0.000000};
+
+/// Lines
+
+// Airfoil 1:
+Spline(1) = {1,2,3};
+Spline(2) = {3,4,5,6,7,8,9,10,11,12,13,14,15};
+Spline(3) = {15,16,17,18,19,20,21,22,23,24,25,26};
+Spline(4) = {26,27,28,29,30,31,32,33,34,35,36,37};
+Spline(5) = {37,38,39,40,41,42,43,44,45,46,47,48,49};
+Spline(6) = {49,50,1};
+
+// Airfoil 2:
+Spline(7) = {52,53,54};
+Spline(8) = {54,55,56,57,58,59,60,61,62,63,64,65,66};
+Spline(9) = {66,67,68,69,70,71,72,73,74,75,76,77};
+Spline(10) = {77,78,79,80,81,82,83,84,85,86,87,88};
+Spline(11) = {88,89,90,91,92,93,94,95,96,97,98,99,100};
+Spline(12) = {100,101,52};
+
+
+// Airfoil 1 to airfoil 2:
+Line(61) = {1,52};
+Line(62) = {3,54};
+Line(63) = {15,66};
+Line(64) = {26,77};
+Line(65) = {37,88};
+Line(66) = {49,100};
+
+// Midplane:
+Line(131) = {5351,5352};
+
+// Wing to midplane:
+Line(161) = {1,5351};
+Line(162) = {52,5352};
+/// Line loops & Surfaces
+
+// Wing 1:
+Line Loop(11) = {1,62,-7,-61};
+Line Loop(12) = {2,63,-8,-62};
+Line Loop(13) = {3,64,-9,-63};
+Line Loop(14) = {4,65,-10,-64};
+Line Loop(15) = {5,66,-11,-65};
+Line Loop(16) = {6,61,-12,-66};
+Surface(11) = {-11};
+Surface(12) = {-12};
+Surface(13) = {-13};
+Surface(14) = {-14};
+Surface(15) = {-15};
+Surface(16) = {-16};
+
+
+// Wake:
+Line Loop(81) = {161,131,-162,-61};
+Surface(81) = {81};
+
+/// Surface loops & Volumes
+
+
+//// MESHING ALGORITHM
+
+
+//// PHYSICAL GROUPS
+
+// Trailing edge and wake tip
+Physical Line("wakeTip") = {162};
+Physical Line("te") = {61};
+
+// Wing:
+Physical Surface("wing") = {11,12,13};
+Physical Surface("wing_") = {14,15,16};
+
+// Wake:
+Physical Surface("wake") = {81};
+
+Coherence;
+Transfinite Curve {161, 162} = nWake Using Progression progWake;
+Transfinite Curve {131, 61, 62, 66, 63, 65, 64} = nSpan Using Progression progSpan;
+Transfinite Curve {-3, 4, -9, 10} = nLe Using Progression progLe;
+Transfinite Curve {7, -12, 1, -6} = nTe Using Progression progTe;
+Transfinite Curve {2, 5, 8, 11} = nMid Using Progression progMid;
+Transfinite Surface {81};
+Transfinite Surface {12};
+Transfinite Surface {15};
+Transfinite Surface {13};
+Transfinite Surface {14};
+Transfinite Surface {11};
+Transfinite Surface {16};
+
+Recombine Surface {81};
+Recombine Surface {12};
+Recombine Surface {15};
+Recombine Surface {13};
+Recombine Surface {14};
+Recombine Surface {11};
+Recombine Surface {16};
diff --git a/blast/src/DDiscretization.cpp b/blast/src/DDiscretization.cpp
index 032e2d843f47d3b6e7cbee2602fd124dd2206d30..9a6796b9d7c502e039873550c9c0cefbbabb3848 100644
--- a/blast/src/DDiscretization.cpp
+++ b/blast/src/DDiscretization.cpp
@@ -17,6 +17,8 @@ void Discretization::setGlob(std::vector<double> &_x, std::vector<double> &_y, s
{
if (_x.size() != _y.size() || _y.size() != _z.size() || _x.size() != _z.size())
throw std::runtime_error("blast::Discretization Wrong mesh sizes.\n");
+ if (_x.size() < 2 || _y.size() < 2 || _z.size() < 2)
+ throw std::runtime_error("blast::Discretization Mesh too small.\n");
nMarkers = _x.size();
x.resize(nMarkers, 0.);
diff --git a/blast/src/DDriver.cpp b/blast/src/DDriver.cpp
index 08a2e328651df0cee65513e5d7fc91a7a3061421..a1fcba0bbb46a638a8ea896cfa1f40ab71460ec5 100644
--- a/blast/src/DDriver.cpp
+++ b/blast/src/DDriver.cpp
@@ -393,6 +393,15 @@ void Driver::computeSectionalDrag(std::vector<BoundaryLayer *> const &bl, size_t
Cdp_sec[i] = Cdt_sec[i] - Cdf_sec[i];
}
+double Driver::getAvrgxoctr(size_t const iRegion) const
+{
+ double xtr = 0.0;
+ for (size_t iSec = 0; iSec < sections.size(); ++iSec)
+ xtr += sections[iSec][iRegion]->xoctr;
+ xtr /= sections.size();
+ return xtr;
+}
+
/**
* @brief Compute total drag coefficient on the airfoil/wing.
* Performs the sectional drag integration for 3D computations.
diff --git a/blast/src/DDriver.h b/blast/src/DDriver.h
index a95458659929bf0bce46f86246f5f2acfa516de1..2937d3dc529e53ed3da7127fc151dd100732626b 100644
--- a/blast/src/DDriver.h
+++ b/blast/src/DDriver.h
@@ -45,6 +45,7 @@ public:
// Getters.
double getxtr(size_t iSec, size_t iRegion) const { return sections[iSec][iRegion]->xtr; };
double getxoctr(size_t iSec, size_t iRegion) const { return sections[iSec][iRegion]->xoctr; };
+ double getAvrgxoctr(size_t const iRegion) const;
double getRe() const { return Re; };
std::vector<std::vector<double>> getSolution(size_t iSec);
diff --git a/blast/utils.py b/blast/utils.py
index 80ac5bb7c607e81879be93f73c72d6c1fee3166a..91a0aa4699a2016cd0ef2f666f7e5b5c7befa25c 100644
--- a/blast/utils.py
+++ b/blast/utils.py
@@ -67,6 +67,8 @@ def initBlast(iconfig, vconfig, iSolver='DART'):
if 'nSections' not in vconfig:
vconfig['nSections'] = len(vconfig['sections'])
+ if 'sfx' not in vconfig:
+ vconfig['sfx'] = ''
# Viscous solver
vSolver = initBL(vconfig['Re'], vconfig['Minf'], vconfig['CFL0'], vconfig['nSections'], xtrF=vconfig['xtrF'])
@@ -80,7 +82,7 @@ def initBlast(iconfig, vconfig, iSolver='DART'):
# Coupler
import blast.coupler as blastCoupler
- coupler = blastCoupler.Coupler(iSolverAPI, vSolver, _maxCouplIter=vconfig['couplIter'], _couplTol=vconfig['couplTol'], _iterPrint=vconfig['iterPrint'], _resetInv=vconfig['resetInv'])
+ coupler = blastCoupler.Coupler(iSolverAPI, vSolver, _maxCouplIter=vconfig['couplIter'], _couplTol=vconfig['couplTol'], _iterPrint=vconfig['iterPrint'], _resetInv=vconfig['resetInv'], sfx=vconfig['sfx'])
return coupler, iSolverAPI, vSolver
def mesh(file, pars):
@@ -131,8 +133,8 @@ def getSolution(vSolver, iSec=0):
'y' : solverOutput[17],
'z' : solverOutput[18],
'ueInv' : solverOutput[19],
- 'xtrT' : vSolver.getxtr(iSec, 0),
- 'xtrB' : vSolver.getxtr(iSec, 1),
+ 'xtrT' : vSolver.getAvrgxoctr(0),
+ 'xtrB' : vSolver.getAvrgxoctr(1),
'Cdt_w' : vSolver.Cdt_sec[iSec],
'Cdf' : vSolver.Cdf_sec[iSec],
'Cdp' : vSolver.Cdp_sec[iSec]
diff --git a/blast/validation/agardValidation.py b/blast/validation/agardValidation.py
new file mode 100644
index 0000000000000000000000000000000000000000..a97d4d484f4f25922c0175f5ec2f3d518152aa88
--- /dev/null
+++ b/blast/validation/agardValidation.py
@@ -0,0 +1,182 @@
+#!/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
+# Test the blast implementation on the 3D AGARD wing
+
+# Imports.
+
+import blast.utils as viscUtils
+import numpy as np
+import os.path
+
+from fwk.wutils import parseargs
+import fwk
+from fwk.testing import *
+from fwk.coloring import ccolors
+
+import math
+
+def cfgInviscid(nthrds, verb):
+ # 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.dirname(os.path.abspath(__file__)))) + '/modules/dartflo/dart/models/agard445.geo', # Input file containing the model
+ 'Pars' : {'xL': 7., 'yL': 3., 'zL': 6., 'xO': -3, 'zO': -3, 'msLeRt': 0.0028, 'msTeRt': 0.0056, 'msLeTp': 0.0018, 'msTeTp': 0.0036, 'msF': 1.0}, # parameters for input file model
+ 'Dim' : 3, # problem dimension
+ 'Format' : 'vtk', # 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)
+ 'Wings' : ['wing'], # LIST of names of physical groups containing the lifting surface boundary
+ '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
+ 'AoA' : 1, # freestream angle of attack
+ # Geometry
+ 'S_ref' : 0.35, # reference surface length
+ 'c_ref' : 0.47, # reference chord length
+ 'x_ref' : 0.0, # 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.7e6, # Freestream Reynolds number
+ 'Minf' : 0.96, # Freestream Mach number (used for the computation of the time step only)
+ 'CFL0' : 1, # Inital CFL number of the calculation
+
+ 'sections' : np.linspace(0.026, 0.7, 15),
+ 'writeSections': np.linspace(0.01, 0.76, 15),
+ 'Sym':[0.],
+ 'span': 0.762,
+ 'interpolator': 'Rbf',
+ 'rbftype': 'linear',
+ 'smoothing': 1e-8,
+ 'degree': 0,
+ 'neighbors': 10,
+ 'saveTag': 5,
+
+ 'Verb': verb, # Verbosity level of the solver
+ 'couplIter': 100, # Maximum number of iterations
+ 'couplTol' : 1e-3, # Tolerance of the VII methodology
+ 'iterPrint': 1, # int, number of iterations between outputs
+ 'resetInv' : True, # bool, flag to reset the inviscid calculation at every iteration.
+ 'xtrF' : [0., 0.], # Forced transition location
+ 'nDim' : 3
+ }
+
+def main():
+ # Timer.
+ tms = fwk.Timers()
+ tms['total'].start()
+
+ args = parseargs()
+ icfg = cfgInviscid(args.k, args.verb)
+ vcfg = cfgBlast(args.verb)
+
+ AoAVec = [-1, 0., 1]
+ sfxVec = ['_AoAminus1deg', '_AoA0deg', '_AoA1deg']
+ for i in range(3):
+ vcfg['sfx'] = sfxVec[i]
+ icfg['AoA'] = AoAVec[i]
+
+ parsViscous = {'nLe': 15, 'nMid': 30, 'nTe': 7, 'nSpan': 60, 'nWake': 30,
+ 'progLe': 1.07, '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/agard445_visc.geo', parsViscous)
+ vcfg['vMsh'] = vMsh
+
+ 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)
+
+ # Write results to file.
+ for iSec in range(len(iSolverAPI.cfg['EffSections'])):
+ vSolution = viscUtils.getSolution(vSolver, iSec)
+ viscUtils.write(vSolution, vSolver.getRe(), sfx=sfxVec[i]+'slice'+str(iSec))
+ vSolution['Cdt_int'] = vSolver.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)
+
+ cps = []
+ for s in sfxVec:
+ cps_s = []
+ for i in range(len(vcfg['writeSections'])):
+ cps_s.append(np.loadtxt('/Users/pauldechamps/lab/softwares/blaster/workspace/blast_validation_agardValidation/agard445_viscous'+s+'_slice_'+str(i)+'.dat', delimiter=',', skiprows=1))
+ cps.append(cps_s)
+ # Plotting
+ from matplotlib import pyplot as plt
+ for i in range(len(vcfg['writeSections'])):
+ plt.figure(i)
+ for j in range(len(sfxVec)):
+ plt.plot(cps[j][i][:,3], cps[j][i][:,4], label=sfxVec[j].replace('_', ' '))
+ plt.legend()
+ plt.title('Section ' + str(i))
+ plt.xlabel('x')
+ plt.ylabel('cp')
+ plt.show()
+ # Test solution
+ print(ccolors.ANSI_BLUE + 'PyTesting...' + ccolors.ANSI_RESET)
+ tests = CTests()
+ tests.add(CTest('Cl', iSolverAPI.getCl(), 0.069, 5e-2))
+ tests.add(CTest('Cd', vSolution['Cdt_int'], 0.00498, 1e-3, forceabs=True))
+ tests.add(CTest('Iterations', len(aeroCoeffs), 8, 0, forceabs=True))
+ tests.run()
+
+ # eof
+ print('')
+
+if __name__ == "__main__":
+ main()
diff --git a/blast/validation/lannValidation.py b/blast/validation/lannValidation.py
index 870d67fb74d662404e41d12036708eb9e868712d..f5e865925ca0f310cfdada7ec3eb750ce2877f0e 100644
--- a/blast/validation/lannValidation.py
+++ b/blast/validation/lannValidation.py
@@ -18,19 +18,7 @@
# @author Paul Dechamps
# @date 2022
-# Test the blast implementation. The test case is a compressible attached transitional flow.
-# Tested functionalities;
-# - Time integration.
-# - Two time-marching techniques (agressive and safe).
-# - Transition routines.
-# - Forced transition.
-# - Compressible flow routines.
-# - Laminar and turbulent flow.
-#
-# Untested functionalities.
-# - Separation routines.
-# - Multiple failure case at one iteration.
-# - Response to unconverged inviscid solver.
+# Test the blast implementation on the 3D LANN wing.
# Imports.
diff --git a/blast/validation/oneraValidation.py b/blast/validation/oneraValidation.py
index 3c0d239551906db5425f858a1753f52adf60b6e9..cb9af04aba1ac3f44ca2572ff00b376f3e510fdd 100644
--- a/blast/validation/oneraValidation.py
+++ b/blast/validation/oneraValidation.py
@@ -18,19 +18,7 @@
# @author Paul Dechamps
# @date 2022
-# Test the blast implementation. The test case is a compressible attached transitional flow.
-# Tested functionalities;
-# - Time integration.
-# - Two time-marching techniques (agressive and safe).
-# - Transition routines.
-# - Forced transition.
-# - Compressible flow routines.
-# - Laminar and turbulent flow.
-#
-# Untested functionalities.
-# - Separation routines.
-# - Multiple failure case at one iteration.
-# - Response to unconverged inviscid solver.
+# Test the blast implementation on the 3D ONERA M6 wing.
# Imports.
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()