From 4dc69103ff3760329be31991239257e6d1b12c86 Mon Sep 17 00:00:00 2001
From: Paul Dechamps <paul.dechamps@uliege.be>
Date: Thu, 25 Apr 2024 11:28:18 +0200
Subject: [PATCH] (feat) Refactor of the solver interface
Boundary layer c++ objects are created in python with swig and are manipulated by the interface. The initial condition of the viscous solver was changed. All test cases have been adapted. A decrease in the number of iterations is globally observed. Results are slightly improved especially for the RAE validation case.
---
blast/_src/blastw.h | 3 +-
blast/_src/blastw.i | 5 +-
blast/coupler.py | 94 ++--
blast/interfaces/DDataStructure.py | 179 +++++--
blast/interfaces/DSolversInterface.py | 242 ++++-----
blast/interfaces/dart/DartInterface.py | 54 +-
.../interpolators/DMatchingInterpolator.py | 19 +-
.../interpolators/DRbfInterpolator.py | 13 +-
.../models/references/rae2822_AR138_case6.dat | 103 ++++
blast/src/CMakeLists.txt | 2 +-
blast/src/DBoundaryLayer.cpp | 143 +++---
blast/src/DBoundaryLayer.h | 68 ++-
blast/src/DDiscretization.cpp | 79 ---
blast/src/DDiscretization.h | 54 --
blast/src/DDriver.cpp | 469 +++++-------------
blast/src/DDriver.h | 26 +-
blast/src/DEdge.cpp | 59 ---
blast/src/DEdge.h | 44 --
blast/src/DFluxes.cpp | 14 +-
blast/src/DSolver.cpp | 56 +--
blast/src/DSolver.h | 2 -
blast/src/blast.h | 2 -
blast/tests/dart/naca0012_2D.py | 34 +-
blast/tests/dart/rae2822_3D.py | 28 +-
blast/utils.py | 172 ++++---
blast/validation/lannValidation.py | 21 +-
blast/validation/oneraValidation.py | 21 +-
blast/validation/raeValidation.py | 47 +-
28 files changed, 950 insertions(+), 1103 deletions(-)
create mode 100644 blast/models/references/rae2822_AR138_case6.dat
delete mode 100644 blast/src/DDiscretization.cpp
delete mode 100644 blast/src/DDiscretization.h
delete mode 100644 blast/src/DEdge.cpp
delete mode 100644 blast/src/DEdge.h
diff --git a/blast/_src/blastw.h b/blast/_src/blastw.h
index ec45463..f0b51a2 100644
--- a/blast/_src/blastw.h
+++ b/blast/_src/blastw.h
@@ -14,4 +14,5 @@
* limitations under the License.
*/
-#include "DDriver.h"
\ No newline at end of file
+#include "DDriver.h"
+#include "DBoundaryLayer.h"
\ No newline at end of file
diff --git a/blast/_src/blastw.i b/blast/_src/blastw.i
index ab103f1..a28e208 100644
--- a/blast/_src/blastw.i
+++ b/blast/_src/blastw.i
@@ -21,12 +21,12 @@
%module(docstring=
"'blastw' module: Viscous inviscid interaction for blast 2021-2022
(c) ULiege - A&M",
-directors="0",
threads="1"
) blastw
%{
#include "blast.h"
+#include "DBoundaryLayer.h"
#include "fwkw.h"
#include "tboxw.h"
@@ -44,6 +44,9 @@ threads="1"
%include "blast.h"
%shared_ptr(blast::Driver);
+%shared_ptr(blast::BoundaryLayer);
+%include "DBoundaryLayer.h"
+
%immutable blast::Driver::Re; // read-only variable (no setter)
%immutable blast::Driver::Cdt;
diff --git a/blast/coupler.py b/blast/coupler.py
index b94e05e..91ecd93 100644
--- a/blast/coupler.py
+++ b/blast/coupler.py
@@ -25,78 +25,108 @@ import numpy as np
class Coupler:
def __init__(self, iSolverAPI, vSolver, _maxCouplIter=150, _couplTol=1e-4, _iterPrint=1, _resetInv=False, sfx=''):
- self.iSolverAPI = iSolverAPI
- self.vSolver = vSolver
- self.maxCouplIter = _maxCouplIter
- self.couplTol = _couplTol
- self.resetInv = _resetInv
- self.iterPrint = _iterPrint if self.iSolverAPI.getVerbose() == 0 and self.vSolver.verbose == 0 else 1
+ self.isol = iSolverAPI
+ self.vsol = vSolver
+ self.maxIter = _maxCouplIter
+ self.tol = _couplTol
+ self.resetInviscid = _resetInv
+ self.iterPrint = _iterPrint if self.isol.getVerbose() == 0 and self.vsol.verbose == 0 else 1
self.tms = fwk.Timers()
self.filesfx = sfx
def run(self):
+ print('')
+ print(' ######################################################## ')
+ print('| ___ __ ___ _________________ |')
+ print('| / _ )/ / / _ | / __/_ __/ __/ _ \ |')
+ print('| / _ / /__/ __ |_\ \ / / / _// , _/ |')
+ print('| /____/____/_/ |_/___/ /_/ /___/_/|_| |')
+ print('| |')
+ print('| 2021-2024 University of Liege |')
+ print('| Paul Dechamps (gitlab.uliege.be/pdechamps) |')
+ print(' ######################################################## ')
+ print('')
+ print('-- Case definition --')
+ print('{:<6s} {:<4.2f} (deg)\n{:<6s} {:<6.2f}\n{:<6s} {:<4.2f}e6'.format('AoA:', self.isol.getAoA()*180/math.pi, 'Mach:', self.isol.getMinf(), 'Re:', self.vsol.getRe()/1e6))
+ print('Transitions: [{:<.2f}, {:<.2f}]'.format(self.isol.cfg['xtrF'][0], self.isol.cfg['xtrF'][1]))
+ print('nDim:', self.isol.getnDim())
+ print('nSections:', self.isol.cfg['nSections'])
+ print('Interpolator:', self.isol.cfg['interpolator'])
+ print('')
+ print('-- Coupler parameters --')
+ print('MaxIter:{:>3.0f}'.format(self.maxIter))
+ print('Tolerance:{:>6.2f}'.format(np.log10(self.tol)))
+ print('IterPrint:{:>2.0f}'.format(self.iterPrint))
+ print('ResetInv:{:>5s}'.format(str(self.resetInviscid)))
+ print('Verb:{:>2.0f}'.format(self.isol.getVerbose()))
+ print('')
+
# Aerodynamic coefficients.
- aeroCoeffs = np.zeros((0, 4))
+ aeroCoeffs = {'Cl':[], 'Cd': [], 'Cdwake': []}
# Convergence parameters.
couplIter = 0
cdPrev = 0.0
- while couplIter < self.maxCouplIter:
+ while couplIter < self.maxIter:
+ # Impose blowing boundary condition in the inviscid solver.
+ self.tms['processing'].start()
+ self.isol.getBlowingBoundary(couplIter)
+ self.isol.setBlowingVelocity()
+ self.tms['processing'].stop()
+
# Run inviscid solver.
self.tms['inviscid'].start()
- if self.resetInv:
- self.iSolverAPI.reset()
- self.iSolverAPI.run()
+ if self.resetInviscid:
+ self.isol.reset()
+ self.isol.run()
self.tms['inviscid'].stop()
# Write inviscid Cp distribution.
if couplIter == 0:
- self.iSolverAPI.writeCp(sfx='_inviscid'+self.filesfx)
+ self.isol.initializeViscousSolver()
+ self.isol.writeCp(sfx='_inviscid'+self.filesfx)
# Impose inviscid boundary in the viscous solver.
self.tms['processing'].start()
- self.iSolverAPI.imposeInvBC(couplIter)
+ self.isol.getInviscidBC()
+ self.isol.setViscousSolver()
self.tms['processing'].stop()
# Run viscous solver.
self.tms['viscous'].start()
- exitCode = self.vSolver.run(couplIter)
+ exitCode = self.vsol.run()
self.tms['viscous'].stop()
- # Impose blowing boundary condition in the inviscid solver.
- self.tms['processing'].start()
- self.iSolverAPI.imposeBlwVel()
- self.tms['processing'].stop()
-
- aeroCoeffs = np.vstack((aeroCoeffs, [self.iSolverAPI.getCl(), self.vSolver.Cdt, self.vSolver.Cdf + self.iSolverAPI.getCd(), self.vSolver.Cdf]))
+ aeroCoeffs['Cl'].append(self.isol.getCl())
+ aeroCoeffs['Cd'].append(self.isol.getCd() + self.vsol.Cdf)
+ aeroCoeffs['Cdwake'].append(self.vsol.Cdt)
# Check convergence status.
- #cd = self.vSolver.Cdf + self.iSolverAPI.getCd()
- cd = self.vSolver.Cdt if self.vSolver.Cdt != 0 else self.vSolver.Cdf + self.iSolverAPI.getCd()
+ cd = self.vsol.Cdf + self.isol.getCd()
+ #cd = self.vsol.Cdt if self.vsol.Cdt != 0 else self.vsol.Cdf + self.isol.getCd()
error = abs((cd - cdPrev) / cd) if cd != 0 else 1
- if error <= self.couplTol:
- print(ccolors.ANSI_GREEN, '{:>4.0f}| {:>7.5f} {:>7.5f} {:>7.5f} | {:>6.4f} {:>6.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)
- self.iSolverAPI.writeCp(sfx='_viscous'+self.filesfx)
+ if error <= self.tol:
+ print(ccolors.ANSI_GREEN, '{:>4.0f}| {:>7.5f} {:>7.5f} {:>7.5f} | {:>6.4f} {:>6.4f} | {:>6.3f}\n'.format(couplIter, self.isol.getCl(), self.isol.getCd()+self.vsol.Cdf, self.vsol.Cdt, self.vsol.getAverageTransition(0), self.vsol.getAverageTransition(1), np.log10(error)), ccolors.ANSI_RESET)
+ self.isol.writeCp(sfx='_viscous'+self.filesfx)
return aeroCoeffs
cdPrev = cd
if couplIter == 0:
- print('- AoA: {:<2.2f} Mach: {:<1.1f} Re: {:<2.1f}e6'.format(self.iSolverAPI.getAoA()*180/math.pi, self.iSolverAPI.getMinf(), self.vSolver.getRe()/1e6))
+ print('')
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.getAvrgxoctr(0), self.vSolver.getAvrgxoctr(1), np.log10(error)))
- if self.iSolverAPI.getVerbose() != 0 or self.vSolver.verbose != 0:
+ print(' {:>4.0f}| {:>7.5f} {:>7.5f} {:>7.5f} | {:>6.4f} {:>6.4f} | {:>6.3f}'.format(couplIter, self.isol.getCl(), self.isol.getCd()+self.vsol.Cdf, self.vsol.Cdt, self.vsol.getAverageTransition(0), self.vsol.getAverageTransition(1), np.log10(error)))
+ if self.isol.getVerbose() != 0 or self.vsol.verbose != 0:
print('')
couplIter += 1
- self.tms['processing'].stop()
else:
- print(ccolors.ANSI_RED, 'Maximum number of {:<.0f} coupling iterations reached'.format(self.maxCouplIter), ccolors.ANSI_RESET)
+ print(ccolors.ANSI_RED, 'Maximum number of {:<.0f} coupling iterations reached'.format(self.maxIter), ccolors.ANSI_RESET)
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.getAvrgxoctr(0), self.vSolver.getAvrgxoctr(1), np.log10(error)), ccolors.ANSI_RESET)
- self.iSolverAPI.writeCp(sfx='_viscous'+self.filesfx)
+ print(ccolors.ANSI_RED, '{:>4.0f}| {:>7.5f} {:>7.5f} {:>7.5f} | {:>6.4f} {:>7.4f} | {:>6.3f}\n'.format(couplIter-1, self.isol.getCl(), self.isol.getCd()+self.vsol.Cdf, self.vsol.Cdt, self.vsol.getAverageTransition(0), self.vsol.getAverageTransition(1), np.log10(error)), ccolors.ANSI_RESET)
+ self.isol.writeCp(sfx='_viscous'+self.filesfx)
return aeroCoeffs
diff --git a/blast/interfaces/DDataStructure.py b/blast/interfaces/DDataStructure.py
index c1eec30..8fa91c0 100644
--- a/blast/interfaces/DDataStructure.py
+++ b/blast/interfaces/DDataStructure.py
@@ -1,11 +1,27 @@
import numpy as np
class Group:
- def __init__(self, _name):
+ def __init__(self, _name, _nDim):
+ """Initialize the group.
+
+ args:
+ ----
+ _name: str
+ Name of the group
+ _nDim: int
+ Number of dimensions
+ """
self.name = _name
- self.nUpperPrev = 0
+ self.nDim = _nDim
self.stagPoint = None
def initStructures(self, nPoints):
+ """Initialize the data structures.
+
+ args:
+ ----
+ nPoints: int
+ Number of points in the region
+ """
self.nPoints = nPoints
self.nElems = nPoints - 1
self.nodesCoord = np.zeros((self.nPoints,4))
@@ -19,70 +35,141 @@ class Group:
self.deltaStarExt = np.zeros(self.nPoints)
self.xxExt = np.zeros(self.nPoints)
- def updateVars(self, nodes, _V, _M, _Rho):
- self.nodesCoord = nodes
- self.nPoints = len(self.nodesCoord[:,0])
- self.nElems = self.nPoints - 1
+ def updateVars(self, _V, _M, _Rho):
+ """Update the velocity, Mach number and density of the boundary layer.
+
+ args:
+ ----
+ _V: np.ndarray
+ Velocity
+ _M: np.ndarray
+ Mach number
+ _Rho: np.ndarray
+ Density
+ """
self.V = _V
self.M = _M
self.Rho = _Rho
- self.computeStagPoint()
- def getXUpper(self):
- return np.flip(self.nodesCoord[:self.stagPoint+1,0])
- def getXLower(self):
- return self.nodesCoord[self.stagPoint:, 0]
- def getYUpper(self):
- return np.flip(self.nodesCoord[:self.stagPoint+1,1])
- def getYLower(self):
- return self.nodesCoord[self.stagPoint:, 1]
- def getZUpper(self):
- return np.flip(self.nodesCoord[:self.stagPoint+1,2])
- def getZLower(self):
- return self.nodesCoord[self.stagPoint:, 2]
+ def getNodesCoord(self, reg:str, dim:int):
+ """Return the coordinates of the nodes in the direction dim
+
+ args:
+ ----
+ reg: str
+ Region of the boundary layer (upper, lower, wake)
+ dim: int
+ Direction of the coordinates (x, y, z)
+ """
+ if reg == "wake":
+ return self.nodesCoord[:,dim]
+ elif reg == "upper":
+ return np.flip(self.nodesCoord[:self.stagPoint+1,dim])
+ elif reg == "lower":
+ return self.nodesCoord[self.stagPoint:,dim]
+ else:
+ raise RuntimeError('Invalid region', reg)
+
+ def getVelocity(self, reg:str, dim:int):
+ """Return the velocity in the direction dim
+
+ args:
+ ----
+ reg: str
+ Region of the boundary layer (upper, lower, wake)
+ dim: int
+ Direction of the velocity
+ """
+ if reg == "wake":
+ return self.V[:,dim]
+ elif reg == "upper":
+ return np.flip(self.V[:self.stagPoint+1,dim])
+ elif reg == "lower":
+ return self.V[self.stagPoint:,dim]
+ else:
+ raise RuntimeError('Invalid region', reg)
- def getVelocityXUpper(self):
- return np.flip(self.V[:self.stagPoint+1,0])
- def getVelocityXLower(self):
- return self.V[self.stagPoint:, 0]
- def getVelocityYUpper(self):
- return np.flip(self.V[:self.stagPoint+1,1])
- def getVelocityYLower(self):
- return self.V[self.stagPoint:, 1]
- def getVelocityZUpper(self):
- return np.flip(self.V[:self.stagPoint+1,2])
- def getVelocityZLower(self):
- return self.V[self.stagPoint:, 2]
+ def getVt(self, reg:str):
+ """Return the velocity in the direction tangential
+ to the boundary layer
- def getMachUpper(self):
- return np.flip(self.M[:self.stagPoint+1])
- def getMachLower(self):
- return self.M[self.stagPoint:]
+ args:
+ ----
+ reg: str
+ Region of the boundary layer (upper, lower, wake)
+ """
+ vt = np.sqrt(self.V[:,0]**2 + self.V[:,self.nDim-1]**2)
+ if reg == "wake":
+ return vt
+ elif reg == "upper":
+ return np.flip(vt[:self.stagPoint+1])
+ elif reg == "lower":
+ return vt[self.stagPoint:]
+ else:
+ raise RuntimeError('Invalid region', reg)
+
+ def getMach(self, reg:str):
+ """Return the Mach number at the edge of the boundary layer
+
+ args:
+ ----
+ reg: str
+ Region of the boundary layer (upper, lower, wake)
+ """
+ if reg == "wake":
+ return self.M
+ elif reg == "upper":
+ return np.flip(self.M[:self.stagPoint+1])
+ elif reg == "lower":
+ return self.M[self.stagPoint:]
+ else:
+ raise RuntimeError('Invalid region', reg)
+
+ def getRho(self, reg:str):
+ """Return the density at the edge of the boundary layer
- def getRhoUpper(self):
- return np.flip(self.Rho[:self.stagPoint+1])
- def getRhoLower(self):
- return self.Rho[self.stagPoint:]
+ args:
+ ----
+ reg: str
+ Region of the boundary layer (upper, lower, wake)
+ """
+ if reg == "wake":
+ return self.Rho
+ elif reg == "upper":
+ return np.flip(self.Rho[:self.stagPoint+1])
+ elif reg == "lower":
+ return self.Rho[self.stagPoint:]
+ else:
+ raise RuntimeError('Invalid region', reg)
- def setBlowingVelocity(self, blwVel):
- self.blowingVel = blwVel
def setConnectList(self, _connectListNodes, _connectListElems):
+ """Set connectivity lists for viscous structures.
+
+ args:
+ ----
+ _connectListNodes: np.ndarray
+ Connectivity list for nodes
+ _connectListElems: np.ndarray
+ Connectivity list for elements
+ """
if self.name != 'iWing' and self.name != 'iWake':
raise RuntimeError('Cannot set connectivity lists for viscous structures', self.name)
self.connectListNodes = _connectListNodes
self.connectListElems = _connectListElems
+
def connectBlowingVel(self):
+ """Connect blowing velocities for viscous structures
+ using the connectivity list.
+ """
if self.name != 'iWing' and self.name != 'iWake':
raise RuntimeWarning('Can not connect blowing velocities for viscous structure. Maybe it was not set.', self.name)
self.blowingVel = self.blowingVel[self.connectListElems.argsort()]
- def isMeshChanged(self):
- return self.nPoints - self.nUpperPrev
- def updatePrev(self):
- self.nUpperPrev = self.nPoints
def computeStagPoint(self):
+ """Compute the stagnation point of the boundary layer.
+ """
if 'iWake' in self.name or 'vWake' in self.name:
self.stagPoint = None
else:
if self.stagPoint is None:
- self.stagPoint = np.argmin(np.sqrt(self.V[:,0]**2 + self.V[:,1]**2))
\ No newline at end of file
+ self.stagPoint = np.argmin(np.sqrt(self.V[:,0]**2 + self.V[:,self.nDim-1]**2))
\ No newline at end of file
diff --git a/blast/interfaces/DSolversInterface.py b/blast/interfaces/DSolversInterface.py
index 5722230..c1f9ed1 100644
--- a/blast/interfaces/DSolversInterface.py
+++ b/blast/interfaces/DSolversInterface.py
@@ -1,5 +1,5 @@
import numpy as np
-from scipy.interpolate import interp1d
+import blast
from blast.interfaces.DDataStructure import Group
import time
@@ -10,133 +10,157 @@ class SolversInterface:
if 'nSections' not in self.cfg:
self.cfg['nSections'] = len(self.cfg['sections'])
-
- if 'sweep' not in self.cfg:
- self.cfg['sweep'] = False
# Initialize data structures.
- self.iBnd = [Group('iWing' if iReg == 0 else 'iWake') for iReg in range(2)]
- self.vBnd = [[Group('vAirfoil'+str(k) if iReg == 0 else 'vWake'+str(k)) for iReg in range(2)] for k in range(self.cfg['nSections'])]
+ self.iBnd = [Group('iWing' if iReg == 0 else 'iWake', self.getnDim()) for iReg in range(2)]
+ self.vBnd = [[Group('vAirfoil'+str(k) if iReg == 0 else 'vWake'+str(k), self.getnDim()) for iReg in range(2)] for k in range(self.cfg['nSections'])]
+ # Initialize interpolator and meshes
+ initialTime = time.time()
+ self.setMesh()
if self.cfg['interpolator'] == 'Matching':
from blast.interfaces.interpolators.DMatchingInterpolator import MatchingInterpolator as interp
+ # Initialize viscous structures for matching computations
+ for iSec in range(len(self.vBnd)):
+ for iReg, reg in enumerate(self.vBnd[iSec]):
+ reg.initStructures(self.iBnd[iReg].nPoints)
+ self.vBnd[iSec][iReg].nodesCoord = self.iBnd[iReg].nodesCoord
+
elif self.cfg['interpolator'] == 'Rbf':
from blast.interfaces.interpolators.DRbfInterpolator import RbfInterpolator as interp
print('RBF interpolator initialized.')
print('Loading viscous mesh... ', end='')
- initialTime = time.time()
- if self.cfg['nDim'] == 2:
- self.getVAirfoil()
- elif self.cfg['nDim'] == 3:
- self.getVWing()
+ self.getVWing()
print('done in {:.2f}s.'.format(time.time()-initialTime))
else:
raise RuntimeError('Incorrect interpolator specified in config.')
- self.interpolator = interp(self.cfg)
-
- print('Loading inviscid mesh... ', end='')
- initialTime = time.time()
- self.setMesh()
- print('done in {:.2f}s.'.format(time.time()-initialTime))
- #self.interpolator.ilwIdx = self.ilwIdx
- #self.interpolator.vlwIdx = self.vlwIdx
+ self.interpolator = interp(self.cfg, self.getnDim())
+ # Initialize transfer quantities
self.deltaStarExt = [[np.zeros(0) for iReg in range(3)]\
for iSec in range(self.cfg['nSections'])]
self.xxExt = [[np.zeros(0) for iReg in range(3)]\
for iSec in range(self.cfg['nSections'])]
- #self.vtOld = [[np.zeros(0) for iReg in range(3)] for iSec in range(self.cfg['nSections'])]
+ self.vtExt = [[np.zeros(0) for iReg in range(3)]\
+ for iSec in range(self.cfg['nSections'])]
- def setViscousSolver(self, couplIter):
+ def setViscousSolver(self, adjoint=False):
"""Interpolate solution to viscous mesh and sets the inviscid boundary in the viscous solver.
"""
- self.interpolator.inviscidToViscous(self.iBnd, self.vBnd, couplIter)
+ self.interpolator.inviscidToViscous(self.iBnd, self.vBnd)
+
for iSec in range(self.cfg['nSections']):
- ## Mesh
- if (self.vBnd[iSec][0].isMeshChanged()):
- for reg in self.vBnd[iSec]:
+ for i, reg in enumerate(self.sec[iSec]):
+ if reg.name == 'wake':
+ iReg = 1
+ elif reg.name == 'lower' or reg.name == 'upper':
+ iReg = 0
+ else:
+ raise RuntimeError('Invalid region', reg.name)
+
+ reg.setMe(self.vBnd[iSec][iReg].getMach(reg.name))
+ reg.setRhoe(self.vBnd[iSec][iReg].getRho(reg.name))
+ reg.setVt(self.vBnd[iSec][iReg].getVt(reg.name))
+
+ if adjoint == False:
+ reg.setDeltaStarExt(self.deltaStarExt[iSec][i])
+ reg.setxxExt(self.xxExt[iSec][i])
+
+ def getBlowingBoundary(self, couplIter):
+ """ Get blowing boundary condition from the viscous solver.
+ args:
+ ----
+ couplIter: int
+ Current iteration of the coupling loop.
+ """
+ if couplIter != 0:
+ for iSec in range(len(self.vBnd)):
+ for iReg, reg in enumerate(self.vBnd[iSec]):
if 'vWake' in reg.name:
- iReg = 2
- self.vSolver.setGlobMesh(iSec, iReg,
- reg.nodesCoord[:,0],
- reg.nodesCoord[:,1],
- reg.nodesCoord[:,2])
- # External variables
- self.xxExt[iSec][2] = np.zeros(reg.nPoints)
- self.deltaStarExt[iSec][2] = np.zeros(reg.nPoints)
- #self.vtOld[iSec][2] = np.zeros(reg.nPoints)
+ reg.blowingVel = np.asarray(self.sec[iSec][2].blowingVelocity)
elif 'vAirfoil' in reg.name:
- self.vSolver.setGlobMesh(iSec, 0,
- reg.getXUpper(),
- reg.getYUpper(),
- reg.getZUpper())
- self.vSolver.setGlobMesh(iSec, 1,
- abs(reg.getXLower()),
- reg.getYLower(),
- reg.getZLower())
- # External variables
- self.xxExt[iSec][0] = np.zeros(reg.stagPoint+1)
- self.deltaStarExt[iSec][0] = np.zeros(reg.stagPoint+1)
- #self.vtOld[iSec][0] = np.zeros(reg.stagPoint+1)
- self.xxExt[iSec][1] = np.zeros(reg.nPoints
- - reg.stagPoint)
- self.deltaStarExt[iSec][1] = np.zeros(reg.nPoints
- - reg.stagPoint)
- #self.vtOld[iSec][1] = np.zeros(reg.nPoints
- # - reg.stagPoint)
- else:
- raise RuntimeError('Tried to initialized viscous\
- struture but name was', reg.name)
- ## Inviscid variables
+ reg.blowingVel = np.concatenate((np.flip(np.asarray(self.sec[iSec][0].blowingVelocity)),
+ np.asarray(self.sec[iSec][1].blowingVelocity)))
+ for i in range(len(self.sec[iSec])):
+ self.xxExt[iSec][i] = np.asarray(self.sec[iSec][i].loc)
+ self.deltaStarExt[iSec][i] = np.asarray(self.sec[iSec][i].deltaStar)
+
+ self.interpolator.viscousToInviscid(self.iBnd, self.vBnd)
+
+ def initializeViscousSolver(self):
+ """Initialize the viscous sections in the viscous solver.
+ """
+ self.sec = [[] for _ in range(self.cfg['nSections'])]
+ self.vSolver.reset()
+ for i in range(self.cfg['nSections']):
+ # TODO: Add functionality to initialize sections without wake.
+ for j in range(3):
+ if j == 2:
+ name = "wake"
+ xtrF = -1
+ elif j == 1:
+ name = "lower"
+ xtrF = self.cfg['xtrF'][j]
+ elif j ==0:
+ name = "upper"
+ xtrF = self.cfg['xtrF'][j]
+
+
+ self.sec[i].append(blast.BoundaryLayer(xtrF, name))
+ self.vSolver.addSection(i, self.sec[i][j])
+
+ self.getInviscidBC()
+ self.interpolator.inviscidToViscous(self.iBnd, self.vBnd)
+
+ # Compute stagnation points
+ for iSec in range(self.cfg['nSections']):
for reg in self.vBnd[iSec]:
- reg.updatePrev()
- if 'vWake' in reg.name:
- iReg = 2
- self.vSolver.setVelocities(iSec, iReg,
- reg.V[:,0],
- reg.V[:,1],
- reg.V[:,2])
- self.vSolver.setMe(iSec, iReg, reg.M)
- self.vSolver.setRhoe(iSec, iReg, reg.Rho)
- elif 'vAirfoil' in reg.name:
- iReg = 0 # Upper side
- self.vSolver.setVelocities(iSec, iReg,
- reg.getVelocityXUpper(),
- reg.getVelocityYUpper(),
- reg.getVelocityZUpper())
- self.vSolver.setMe(iSec, iReg, reg.getMachUpper())
- self.vSolver.setRhoe(iSec, iReg, reg.getRhoUpper())
- iReg = 1 # Lower side
- self.vSolver.setVelocities(iSec, iReg,
- reg.getVelocityXLower(),
- reg.getVelocityYLower(),
- reg.getVelocityZLower())
- self.vSolver.setMe(iSec, iReg, reg.getMachLower())
- self.vSolver.setRhoe(iSec, iReg, reg.getRhoLower())
+ reg.computeStagPoint()
+
+ # Set mesh
+ for iSec in range(self.cfg['nSections']):
+ # Compute section chord
+ xMin = np.min(self.vBnd[iSec][0].nodesCoord[:,0])
+ xMax = np.max(self.vBnd[iSec][0].nodesCoord[:,0])
+ chord = xMax - xMin
+ for reg in self.sec[iSec]:
+ if reg.name == 'wake':
+ iReg = 1
+ elif reg.name == 'lower' or reg.name == 'upper':
+ iReg = 0
else:
- raise RuntimeError('Tried to initialize a viscous\
- struture but name was', reg.name)
- ## External variables
- for iRegv in range(3):
- self.vSolver.setxxExt(iSec, iRegv, self.xxExt[iSec][iRegv])
- self.vSolver.setDeltaStarExt(iSec, iRegv,
- self.deltaStarExt[iSec][iRegv])
- #self.vSolver.setVtOld(iSec, iRegv, self.vtOld[iSec][iRegv])
- def getBlowingBoundary(self):
- for iSec in range(len(self.vBnd)):
- for reg in self.vBnd[iSec]:
- if 'vWake' in reg.name:
- reg.blowingVel = np.asarray(self.vSolver.extractBlowingVel(iSec, 2))
- elif 'vAirfoil' in reg.name:
- reg.blowingVel = np.concatenate((np.flip(np.asarray(self.vSolver.extractBlowingVel(iSec, 0))),
- np.asarray(self.vSolver.extractBlowingVel(iSec, 1))))
- for iRegv in range(3):
- self.xxExt[iSec][iRegv] = self.vSolver.extractxx(iSec, iRegv)
- self.deltaStarExt[iSec][iRegv] = self.vSolver.extractDeltaStar(iSec, iRegv)
- #self.vtOld[iSec][iRegv] = self.vSolver.extractUe(iSec, iRegv)
+ raise RuntimeError('Invalid region', reg.name)
- self.interpolator.viscousToInviscid(self.iBnd, self.vBnd)
+ xIdx = 0
+ if self.getnDim() == 2:
+ yIdx = 1
+ zIdx = 2
+ elif self.getnDim() == 3:
+ yIdx = 2
+ zIdx = 1
+ else:
+ raise RuntimeError('Incorrect number of dimensions', self.getnDim())
+
+ reg.setMesh(self.vBnd[iSec][iReg].getNodesCoord(reg.name, xIdx),
+ self.vBnd[iSec][iReg].getNodesCoord(reg.name, yIdx),
+ self.vBnd[iSec][iReg].getNodesCoord(reg.name, zIdx),
+ chord,
+ xMin)
+
+ # External variables
+ for i, reg in enumerate(self.sec[iSec]):
+ if reg.name == 'wake':
+ iReg = 1
+ elif reg.name == 'lower' or reg.name == 'upper':
+ iReg = 0
+ else:
+ raise RuntimeError('Invalid region', reg.name)
+ loc = np.zeros(reg.getnNodes())
+ for iPoint in range(reg.getnNodes()):
+ loc[iPoint] = reg.loc[iPoint]
+ self.xxExt[iSec][i] = loc
+ self.deltaStarExt[iSec][i] = np.zeros(reg.getnNodes())
def getVWing(self):
"""Obtain the nodes' location and row and cg of all elements.
@@ -147,11 +171,7 @@ class SolversInterface:
"""
data = [np.zeros((0,4)) for _ in range(2)]
-
- leNodes = np.zeros((0, 3))
lwRows = []
- ctlw = 0
-
for e in self.cfg['vMsh'].elems:
if e.ptag.name == 'wing' or e.ptag.name == 'wing_'\
or e.ptag.name == 'leadingEdge':
@@ -165,8 +185,8 @@ class SolversInterface:
for nw in e.nodes:
if nw.row not in data[iRegion][:,3]:
data[iRegion] = np.vstack((data[iRegion], [nw.pos[0],\
- nw.pos[2],\
nw.pos[1],\
+ nw.pos[2],\
nw.row]))
# If wing_ exists, the lower side is identifiable.
if e.ptag.name == 'wing_':
@@ -188,16 +208,16 @@ class SolversInterface:
for iy, y in enumerate(self.cfg['sections']):
# Find nearest value in the mesh
try:
- idx = (np.abs(reg[:,2]-y).argmin())
+ idx = (np.abs(reg[:,1]-y).argmin())
except:
- print('Region', reg[:,2], 'cannot be found through\
+ print('Region', reg[:,1], 'cannot be found through\
config input', y)
raise RuntimeError("Could not identify section.\
Possibly a missing tag")
- section = reg[idx,2]
+ section = reg[idx,1]
if iReg == 0:
self.cfg['EffSections'] = np.append(self.cfg['EffSections'], section)
- nodesSec = reg[abs(reg[:,2]-section)<1e-3, :]
+ nodesSec = reg[abs(reg[:,1]-section)<1e-3, :]
# Separate points on upper and lower side
if iReg == 0:
@@ -225,4 +245,4 @@ class SolversInterface:
self.vBnd[iy][iReg].initStructures(nodesSec.shape[0])
self.vBnd[iy][iReg].nodesCoord = nodesSec
- self.vBnd[iy][iReg].elemsCoord = cgSec
+ self.vBnd[iy][iReg].elemsCoord = cgSec
\ No newline at end of file
diff --git a/blast/interfaces/dart/DartInterface.py b/blast/interfaces/dart/DartInterface.py
index da0809b..9c6f040 100644
--- a/blast/interfaces/dart/DartInterface.py
+++ b/blast/interfaces/dart/DartInterface.py
@@ -28,9 +28,9 @@ class DartInterface(SolversInterface):
def __init__(self, dartCfg, vSolver, _cfg):
try:
from modules.dartflo.dart.api.core import initDart
- self.argOut = initDart(dartCfg, viscous=True)
- self.solver = self.argOut.get('sol') # Solver
- self.blw = [self.argOut.get('blwb'), self.argOut.get('blww')]
+ self.inviscidObjects = initDart(dartCfg, viscous=True)
+ self.solver = self.inviscidObjects.get('sol') # Solver
+ self.blw = [self.inviscidObjects.get('blwb'), self.inviscidObjects.get('blww')]
print('Dart successfully loaded.')
except:
print(ccolors.ANSI_RED, 'Could not load DART. Make sure it is installed.', ccolors.ANSI_RESET)
@@ -50,7 +50,7 @@ class DartInterface(SolversInterface):
# Extract Cp on elements
self.save(sfx=sfx)
# 2D Case
- if self.cfg['nDim'] == 2:
+ if self.getnDim() == 2:
vElems = self.blw[0].tag.elems
vData = self.solver.Cp
if isinstance(vData[0], float):
@@ -75,25 +75,22 @@ class DartInterface(SolversInterface):
print('writing data file ' + 'Cp' +sfx + '.dat')
np.savetxt('Cp'+sfx+'.dat', data, fmt='%1.6e', delimiter=', ', header='x, y, z, Cp', comments='')
- elif self.cfg['nDim'] == 3:
+ elif self.getnDim() == 3:
# Save surface cp
- data = np.zeros((self.argOut['bnd'].nodes.size(), 4))
- for i, n in enumerate(self.argOut['bnd'].nodes):
+ data = np.zeros((self.inviscidObjects['bnd'].nodes.size(), 4))
+ for i, n in enumerate(self.inviscidObjects['bnd'].nodes):
data[i,0] = n.pos[0]
data[i,1] = n.pos[1]
data[i,2] = n.pos[2]
data[i,3] = self.solver.Cp[n.row]
- np.savetxt(self.argOut['msh'].name+'_Cp_surface'+sfx+'.dat', data, fmt='%1.6e', delimiter=', ', header='x, y, z, Cp', comments='')
+ np.savetxt(self.inviscidObjects['msh'].name+'_Cp_surface'+sfx+'.dat', data, fmt='%1.6e', delimiter=', ', header='x, y, z, Cp', comments='')
import modules.dartflo.dart.utils as invUtils
if self.cfg['Format'] == 'vtk':
- import os
- if not os.path.exists('cpSlices'+sfx):
- os.makedirs('cpSlices'+sfx)
- print('Saving Cp files in vtk format. Msh {:s}, Tag {:.0f}'.format(self.argOut['msh'].name, self.cfg['saveTag']))
- invUtils.writeSlices(self.argOut['msh'].name, self.cfg['writeSections'], self.cfg['saveTag'], sfx=sfx)
+ print('Saving Cp files in vtk format. Msh {:s}, Tag {:.0f}'.format(self.inviscidObjects['msh'].name, self.cfg['saveTag']))
+ invUtils.writeSlices(self.inviscidObjects['msh'].name, self.cfg['writeSections'], self.cfg['saveTag'], sfx=sfx)
def save(self, sfx='inviscid'):
- self.solver.save(self.argOut['wrt'], sfx)
+ self.solver.save(self.inviscidObjects['wrt'], sfx)
def getAoA(self):
return self.solver.pbl.alpha
@@ -113,10 +110,13 @@ class DartInterface(SolversInterface):
def getVerbose(self):
return self.solver.verbose
+ def getnDim(self):
+ return self.inviscidObjects['pbl'].nDim
+
def reset(self):
self.solver.reset()
- def imposeInvBC(self, couplIter):
+ def getInviscidBC(self):
""" Extract the inviscid paramaters at the edge of the boundary layer
and use it as a boundary condition in the viscous solver.
@@ -132,27 +132,23 @@ class DartInterface(SolversInterface):
self.iBnd[n].V[iRow,iDim] = self.solver.U[row][iDim]
self.iBnd[n].M[iRow] = self.solver.M[row]
self.iBnd[n].Rho[iRow] = self.solver.rho[row]
- if self.cfg['nDim']==3:
- self.iBnd[n].V[:,[1,2]] = self.iBnd[n].V[:,[2,1]]
- self.setViscousSolver(couplIter)
- def imposeBlwVel(self):
+ def setBlowingVelocity(self):
""" Extract the solution of the viscous calculation (blowing velocity)
and use it as a boundary condition in the inviscid solver.
"""
- self.getBlowingBoundary()
+
# Impose blowing velocities.
for n in range(len(self.iBnd)):
- if self.cfg['nDim'] == 2:
+ if self.getnDim() == 2:
self.iBnd[n].connectBlowingVel()
for i, blw in enumerate(self.iBnd[n].blowingVel):
self.blw[n].setU(i, blw)
def setMesh(self):
-
- if self.cfg['nDim'] == 2:
+ if self.getnDim() == 2:
self.getAirfoil()
- elif self.cfg['nDim'] == 3:
+ elif self.getnDim() == 3:
self.getWing()
else:
raise RuntimeError('dartInterface::Could not resolve how to set\
@@ -326,11 +322,11 @@ class DartInterface(SolversInterface):
# Compute cg position
cg_pt = np.zeros(3)
for nw in e.nodes:
- cg_pt += [nw.pos[0], nw.pos[2], nw.pos[1]]
+ cg_pt += [nw.pos[0], nw.pos[1], nw.pos[2]]
if nw.row not in data[iRegion][:,3]:
data[iRegion] = np.vstack((data[iRegion], [nw.pos[0],\
- nw.pos[2],\
nw.pos[1],\
+ nw.pos[2],\
nw.row]))
cg_pt/=e.nodes.size()
cg_pt = np.hstack((cg_pt, e.no))
@@ -341,9 +337,9 @@ class DartInterface(SolversInterface):
for s in self.cfg['Sym']:
for iReg in range(len(data)):
dummy = data[iReg].copy()
- dummy[:,2] -= s
- dummy[:,2] *= -1
- dummy[:,2] += s
+ dummy[:,1] -= s
+ dummy[:,1] *= -1
+ dummy[:,1] += s
data[iReg] = np.row_stack((data[iReg], dummy))
for n in range(len(data)):
diff --git a/blast/interfaces/interpolators/DMatchingInterpolator.py b/blast/interfaces/interpolators/DMatchingInterpolator.py
index a801890..3a42cd6 100644
--- a/blast/interfaces/interpolators/DMatchingInterpolator.py
+++ b/blast/interfaces/interpolators/DMatchingInterpolator.py
@@ -1,23 +1,22 @@
import numpy as np
class MatchingInterpolator:
- def __init__(self, cfg):
- self.cfg = cfg
+ def __init__(self, _cfg, _nDim):
+ self.cfg = _cfg
+ self.nDim = _nDim
- def inviscidToViscous(self, iDict, vDict, couplIter):
- ## Airfoil
- # Find stagnation point
- if self.cfg['nDim'] == 2:
+ def inviscidToViscous(self, iDict, vDict):
+ if self.nDim == 2:
for iReg in range(len(iDict)):
- vDict[0][iReg].updateVars(iDict[iReg].nodesCoord, iDict[iReg].V, iDict[iReg].M, iDict[iReg].Rho)
- elif self.cfg['nDim'] == 3:
+ vDict[0][iReg].updateVars(iDict[iReg].V, iDict[iReg].M, iDict[iReg].Rho)
+ elif self.nDim == 3:
for iSec, ysec in enumerate(self.cfg['sections']):
for iReg in range(2):
print(iDict[iReg].nodesCoord[iDict[iReg].nodesCoord[:,1] == ysec])
print(iDict[iReg].V[iDict[iReg].nodesCoord[:,1] == ysec])
- vDict[iSec][iReg].updateVars(iDict[iReg].nodesCoord[iDict[iReg].nodesCoord[:,1] == ysec], iDict[iReg].V[iDict[iReg].nodesCoord[:,1] == ysec], iDict[iReg].Rho[iDict[iReg].nodesCoord[:,1] == ysec])
+ vDict[iSec][iReg].updateVars(iDict[iReg].V[iDict[iReg].nodesCoord[:,1] == ysec], iDict[iReg].Rho[iDict[iReg].nodesCoord[:,1] == ysec])
def viscousToInviscid(self, iDict, vDict):
- if self.cfg['nDim'] == 2:
+ if self.nDim == 2:
for iReg in range(2):
iDict[iReg].blowingVel = vDict[0][iReg].blowingVel
else:
diff --git a/blast/interfaces/interpolators/DRbfInterpolator.py b/blast/interfaces/interpolators/DRbfInterpolator.py
index 081ccc9..1eda4a3 100644
--- a/blast/interfaces/interpolators/DRbfInterpolator.py
+++ b/blast/interfaces/interpolators/DRbfInterpolator.py
@@ -3,10 +3,11 @@
import numpy as np
from scipy.interpolate import RBFInterpolator
class RbfInterpolator:
- def __init__(self, _cfg):
+ def __init__(self, _cfg, _nDim):
self.cfg = _cfg
+ self.nDim = _nDim
- def inviscidToViscous(self, iDict, vDict, couplIter):
+ def inviscidToViscous(self, iDict, vDict):
## Airfoil
# Find stagnation point
for iSec in range(len(vDict)):
@@ -18,13 +19,13 @@ class RbfInterpolator:
v[:,iDim] = self.__rbfToSection(iDict[iReg].nodesCoord[:,:(self.cfg['nDim'] if iDict[iReg].name == 'iWing' else 1)], iDict[iReg].V[:,iDim], reg.nodesCoord[:,:(self.cfg['nDim'] if 'vAirfoil' in reg.name else 1)])
M = self.__rbfToSection(iDict[iReg].nodesCoord[:,:(self.cfg['nDim'] if iDict[iReg].name == 'iWing' else 1)], iDict[iReg].M, reg.nodesCoord[:,:(self.cfg['nDim'] if 'vAirfoil' in reg.name else 1)])
rho = self.__rbfToSection(iDict[iReg].nodesCoord[:,:(self.cfg['nDim'] if iDict[iReg].name == 'iWing' else 1)], iDict[iReg].Rho, reg.nodesCoord[:,:(self.cfg['nDim'] if 'vAirfoil' in reg.name else 1)])
- vDict[iSec][iReg].updateVars(vDict[iSec][iReg].nodesCoord, v, M, rho)
+ vDict[iSec][iReg].updateVars(v, M, rho)
def viscousToInviscid(self, iDict, vDict):
- if self.cfg['nDim'] == 2:
+ if self.nDim == 2:
for iReg, reg in enumerate(iDict):
iDict[iReg].blowingVel = self.__rbfToSection(vDict[0][iReg].elemsCoordTr[:,:(self.cfg['nDim']-1 if 'vAirfoil' in reg.name else 1)], vDict[0][iReg].blowingVel, reg.elemsCoordTr[:,:(self.cfg['nDim']-1 if reg.name == 'iWing' else 1)])
- elif self.cfg['nDim'] == 3:
+ elif self.nDim == 3:
for iReg, reg in enumerate(iDict):
viscElemsCoord = np.zeros((0,3))
viscBlowing = np.zeros(0)
@@ -34,7 +35,7 @@ class RbfInterpolator:
if 'Sym' in self.cfg:
for s in self.cfg['Sym']:
dummy = viscElemsCoord.copy()
- dummy[:,2] = (dummy[:,2] - s)*(-1) + s
+ dummy[:,1] = (dummy[:,1] - s)*(-1) + s
viscElemsCoord = np.row_stack((viscElemsCoord, dummy))
viscBlowing = np.concatenate((viscBlowing, viscBlowing))
reg.blowingVel = self.__rbfToSection(viscElemsCoord, viscBlowing, reg.elemsCoord[:,:3])
diff --git a/blast/models/references/rae2822_AR138_case6.dat b/blast/models/references/rae2822_AR138_case6.dat
new file mode 100644
index 0000000..9fc6aa3
--- /dev/null
+++ b/blast/models/references/rae2822_AR138_case6.dat
@@ -0,0 +1,103 @@
+0.9938 -0.1432
+0.9875 -0.1318
+0.9750 -0.1082
+0.9500 -0.0592
+0.9250 -0.0115
+0.9000 0.0365
+0.8750 0.0808
+0.8500 0.1296
+0.8250 0.1746
+0.8000 0.2186
+0.7750 0.2702
+0.7500 0.3029
+0.7000 0.3913
+0.6771 0.4263
+0.6500 0.4778
+0.6196 0.5361
+0.6000 0.5798
+0.5750 0.6769
+0.5500 0.9137
+0.5250 1.2201
+0.5000 1.2164
+0.4750 1.2038
+0.4500 1.1940
+0.4250 1.1819
+0.4000 1.1601
+0.3750 1.1490
+0.3500 1.1273
+0.3250 1.1168
+0.3000 1.1091
+0.2800 1.1010
+0.2500 1.0736
+0.2208 1.0456
+0.2000 1.0433
+0.1500 1.0096
+0.1000 0.9886
+0.0750 1.0432
+0.0625 1.0730
+0.0500 1.0923
+0.0375 1.0820
+0.0271 1.0435
+0.0187 1.0362
+0.0146 0.9204
+0.0125 0.8680
+0.0104 0.8041
+0.0087 0.7658
+0.0073 0.7250
+0.0060 0.6269
+0.0048 0.5309
+0.0036 0.4667
+0.0026 0.2746
+0.0016 0.0521
+0.0008 -0.1042
+0.0002 0.4219
+0.0000 -0.8328
+0.0002 -1.0396
+0.0008 -1.1053
+0.0016 -1.1338
+0.0026 -1.1223
+0.0036 -1.0823
+0.0048 -1.0273
+0.0060 -0.9736
+0.0073 -0.9148
+0.0087 -0.8703
+0.0104 -0.8148
+0.0125 -0.7537
+0.0146 -0.6980
+0.0186 -0.6159
+0.0271 -0.5444
+0.0375 -0.4040
+0.0500 -0.3088
+0.0625 -0.2454
+0.0750 -0.1926
+0.1000 -0.1042
+0.1500 0.0081
+0.2000 0.0940
+0.2500 0.1754
+0.3000 0.2506
+0.3250 0.2904
+0.3500 0.3214
+0.3750 0.3406
+0.4000 0.3403
+0.4250 0.3194
+0.4500 0.2881
+0.4750 0.2469
+0.5000 0.2079
+0.5250 0.1657
+0.5500 0.1169
+0.5750 0.0731
+0.6000 0.0301
+0.6196 0.0003
+0.6500 -0.0431
+0.6771 -0.0836
+0.7000 -0.1134
+0.7500 -0.1731
+0.7750 -0.1988
+0.8500 -0.2618
+0.8750 -0.2778
+0.9000 -0.2907
+0.9250 -0.2972
+0.9500 -0.2937
+0.9750 -0.2676
+0.9875 -0.2396
+0.9938 -0.2146
\ No newline at end of file
diff --git a/blast/src/CMakeLists.txt b/blast/src/CMakeLists.txt
index 394f4cb..96caacb 100644
--- a/blast/src/CMakeLists.txt
+++ b/blast/src/CMakeLists.txt
@@ -20,7 +20,7 @@ ADD_LIBRARY(blast SHARED ${SRCS})
MACRO_DebugPostfix(blast)
TARGET_INCLUDE_DIRECTORIES(blast PUBLIC ${PROJECT_SOURCE_DIR}/blast/src)
-TARGET_LINK_LIBRARIES(blast tbox)
+TARGET_LINK_LIBRARIES(blast dart fwk tbox)
INSTALL(TARGETS blast DESTINATION ${CMAKE_INSTALL_PREFIX})
diff --git a/blast/src/DBoundaryLayer.cpp b/blast/src/DBoundaryLayer.cpp
index 9886a79..86fc294 100644
--- a/blast/src/DBoundaryLayer.cpp
+++ b/blast/src/DBoundaryLayer.cpp
@@ -1,9 +1,8 @@
#include "DBoundaryLayer.h"
#include "DClosures.h"
-#include "DDiscretization.h"
-#include "DEdge.h"
#include <cmath>
#include <iomanip>
+#include <iostream>
using namespace blast;
@@ -12,55 +11,80 @@ using namespace blast;
*
* @param _xtrF Forced transition location (default=-1; free transition).
*/
-BoundaryLayer::BoundaryLayer(double _xtrF)
+BoundaryLayer::BoundaryLayer(double _xtrF, std::string _name)
{
xtrF = _xtrF;
+ name = _name;
xtr = 1.0;
- xoctr = 1.0;
closSolver = new Closures();
- mesh = new Discretization();
- blEdge = new Edge();
- transMarker = 0;
+ transitionNode = 0;
}
BoundaryLayer::~BoundaryLayer()
{
delete closSolver;
- delete mesh;
- delete blEdge;
- // std::cout << "~blast::BoundaryLayer()\n";
}
-/**
- * @brief Set the mesh and resize all boundary layer quantities accordingly.
- *
- * @param x Position x in the global frame of reference.
- * @param y Position y in the global frame of reference.
- * @param z Position z in the global frame of reference.
- */
-void BoundaryLayer::setMesh(std::vector<double> x,
- std::vector<double> y,
- std::vector<double> z)
+void BoundaryLayer::setMesh(std::vector<double> _x, std::vector<double> _y, std::vector<double> _z, double _chord, double xMin)
{
- size_t nMarkers = x.size();
- mesh->setGlob(x, y, z);
-
- // Resize and reset all variables if needed.
- if (mesh->getDiffnElm() != 0)
+ if (_x.size() != _y.size() || _x.size() != _z.size())
+ throw std::runtime_error("Mesh coordinates are not consistent.");
+ nNodes = _x.size();
+ nElms = nNodes - 1;
+ x = _x;
+ y = _y;
+ z = _z;
+ chord = _chord;
+
+ // Compute the local coordinate.
+ loc.resize(nNodes, 0.);
+ alpha.resize(nElms, 0.);
+ dx.resize(nElms, 0.);
+
+ for (size_t iPoint = 0; iPoint < nElms; ++iPoint)
{
- u.resize(nVar * nMarkers, 0.);
- Ret.resize(nMarkers, 0.);
- cd.resize(nMarkers, 0.);
- cf.resize(nMarkers, 0.);
- Hstar.resize(nMarkers, 0.);
- Hstar2.resize(nMarkers, 0.);
- Hk.resize(nMarkers, 0.);
- ctEq.resize(nMarkers, 0.);
- us.resize(nMarkers, 0.);
- deltaStar.resize(nMarkers, 0.);
- delta.resize(nMarkers, 0.);
- regime.resize(nMarkers, 0);
+ alpha[iPoint] = std::atan2(y[iPoint + 1] - y[iPoint], x[iPoint + 1] - x[iPoint]);
+ // dx = sqrt(delta x ^2 + delta y ^2).
+ dx[iPoint] = std::sqrt((x[iPoint + 1] - x[iPoint]) * (x[iPoint + 1] - x[iPoint]) + (y[iPoint + 1] - y[iPoint]) * (y[iPoint + 1] - y[iPoint]));
+ loc[iPoint + 1] = loc[iPoint] + dx[iPoint];
}
+
+ // Compute the x/c coordinate.
+ if (chord <= 0)
+ throw std::runtime_error("Chord length is zero or negative.");
+ xoc.resize(nNodes, 0.);
+ for (size_t iPoint = 0; iPoint < nNodes; ++iPoint)
+ xoc[iPoint] = (x[iPoint] - xMin) / chord;
+
+ this->reset();
+}
+
+void BoundaryLayer::reset()
+{
+ u.resize(nNodes*getnVar(), 0.);
+ Ret.resize(nNodes, 0.);
+ cd.resize(nNodes, 0.);
+ cf.resize(nNodes, 0.);
+ Hstar.resize(nNodes, 0.);
+ Hstar2.resize(nNodes, 0.);
+ Hk.resize(nNodes, 0.);
+ ctEq.resize(nNodes, 0.);
+ us.resize(nNodes, 0.);
+ deltaStar.resize(nNodes, 0.);
+ delta.resize(nNodes, 0.);
+ deltaStar.resize(nNodes, 0.);
+ regime.resize(nNodes, 0);
+ blowingVelocity.resize(nElms, 0.);
+
+ Me.resize(nNodes, 0.);
+ vt.resize(nNodes, 0.);
+ rhoe.resize(nNodes, 0.);
+ deltaStarExt.resize(nNodes, 0.);
+ vtExt.resize(nNodes, 0.);
+ xxExt.resize(nNodes, 0.);
+
+ transitionNode = 0;
+ xtr = 1.0;
}
/**
@@ -68,16 +92,16 @@ void BoundaryLayer::setMesh(std::vector<double> x,
*
* @param Re Freestream Reynolds number.
*/
-void BoundaryLayer::setStagBC(double Re)
+void BoundaryLayer::setStagnationBC(double Re)
{
- double dv0 = (blEdge->getVt(1) - blEdge->getVt(0)) / (mesh->getLoc(1) - mesh->getLoc(0));
+ double dv0 = (vt[1] - vt[0]) / (loc[1] - loc[0]);
if (dv0 > 0)
u[0] = sqrt(0.075 / (Re * dv0)); // Theta
else
u[0] = 1e-6; // Theta
u[1] = 2.23; // H
u[2] = 0.; // N
- u[3] = blEdge->getVt(0); // ue
+ u[3] = vt[0]; // ue
u[4] = 0.; // Ctau
Ret[0] = u[3] * u[0] * Re;
@@ -89,7 +113,7 @@ void BoundaryLayer::setStagBC(double Re)
deltaStar[0] = u[0] * u[1];
std::vector<double> lamParam(8, 0.);
- closSolver->laminarClosures(lamParam, u[0], u[1], Ret[0], blEdge->getMe(0), name);
+ closSolver->laminarClosures(lamParam, u[0], u[1], Ret[0], Me[0], name);
Hstar[0] = lamParam[0];
Hstar2[0] = lamParam[1];
@@ -115,7 +139,7 @@ void BoundaryLayer::setWakeBC(double Re, std::vector<double> UpperCond, std::vec
u[0] = UpperCond[0] + LowerCond[0]; // (Theta_up+Theta_low).
u[1] = (UpperCond[0] * UpperCond[1] + LowerCond[0] * LowerCond[1]) / u[0]; // ((Theta*H)_up+(Theta*H)_low)/(Theta_up+Theta_low).
u[2] = 9.; // Amplification factor.
- u[3] = blEdge->getVt(0); // Inviscid velocity.
+ u[3] = vt[0]; // Inviscid velocity.
u[4] = (UpperCond[0] * UpperCond[4] + LowerCond[0] * LowerCond[4]) / u[0]; // ((Ct*Theta)_up+(Theta)_low)/(Ct*Theta_up+Theta_low).
Ret[0] = u[3] * u[0] * Re; // Reynolds number based on the momentum thickness.
@@ -129,7 +153,7 @@ void BoundaryLayer::setWakeBC(double Re, std::vector<double> UpperCond, std::vec
// Laminar closures.
std::vector<double> lamParam(8, 0.);
- closSolver->laminarClosures(lamParam, u[0], u[1], Ret[0], blEdge->getMe(0), name);
+ closSolver->laminarClosures(lamParam, u[0], u[1], Ret[0], Me[0], name);
Hstar[0] = lamParam[0];
Hstar2[0] = lamParam[1];
Hk[0] = lamParam[2];
@@ -139,7 +163,7 @@ void BoundaryLayer::setWakeBC(double Re, std::vector<double> UpperCond, std::vec
ctEq[0] = lamParam[6];
us[0] = lamParam[7];
- for (size_t k = 0; k < mesh->getnMarkers(); ++k)
+ for (size_t k = 0; k < nNodes; ++k)
regime[k] = 1;
}
@@ -147,20 +171,20 @@ void BoundaryLayer::setWakeBC(double Re, std::vector<double> UpperCond, std::vec
* @brief Compute the blowing velocity in each station of the region.
*
*/
-void BoundaryLayer::computeBlwVel()
+void BoundaryLayer::computeBlowingVelocity()
{
deltaStar[0] = u[0] * u[1];
- blowingVelocity.resize(mesh->getnMarkers() - 1, 0.);
- for (size_t iPoint = 1; iPoint < mesh->getnMarkers(); ++iPoint)
+ blowingVelocity.resize(nNodes - 1, 0.);
+ for (size_t iPoint = 1; iPoint < nNodes; ++iPoint)
{
deltaStar[iPoint] = u[iPoint * nVar + 0] * u[iPoint * nVar + 1];
- blowingVelocity[iPoint - 1] = (blEdge->getRhoe(iPoint) * blEdge->getVt(iPoint) * deltaStar[iPoint] - blEdge->getRhoe(iPoint - 1) * blEdge->getVt(iPoint - 1) * deltaStar[iPoint - 1]) / (blEdge->getRhoe(iPoint) * (mesh->getLoc(iPoint) - mesh->getLoc(iPoint - 1)));
+ blowingVelocity[iPoint - 1] = (rhoe[iPoint] * vt[iPoint] * deltaStar[iPoint] - rhoe[iPoint-1] * vt[iPoint-1] * deltaStar[iPoint - 1]) / (rhoe[iPoint] * (loc[iPoint] - loc[iPoint-1]));
if (std::isnan(blowingVelocity[iPoint - 1]))
{
- if (blEdge->getRhoe(iPoint) == 0.)
+ if (rhoe[iPoint] == 0.)
std::cout << "Density is zero at point " << iPoint << std::endl;
- if ((mesh->getLoc(iPoint) - mesh->getLoc(iPoint - 1)) == 0.)
- std::cout << "Points " << iPoint -1 << " and " << iPoint << " are at the same position " << mesh->getLoc(iPoint) << ", resulting in an infinite gradient." << std::endl;
+ if ((loc[iPoint] - loc[iPoint-1]) == 0.)
+ std::cout << "Points " << iPoint -1 << " and " << iPoint << " are at the same position " << loc[iPoint] << ", resulting in an infinite gradient." << std::endl;
if (std::isnan(deltaStar[iPoint]))
std::cout << "Displacement thickness is nan at position " << iPoint << std::endl;
if (std::isnan(deltaStar[iPoint - 1]))
@@ -173,15 +197,16 @@ void BoundaryLayer::computeBlwVel()
/**
* @brief Print solution at a given station.
*
- * @param iPoint Marker id.
+ * @param iPoint Node id.
* @note Function used for debugging.
*/
void BoundaryLayer::printSolution(size_t iPoint) const
{
- if (iPoint < 0 || iPoint > mesh->getnMarkers()-1)
+ if (iPoint < 0 || iPoint > nNodes-1)
throw std::runtime_error("Tried to access element outside of region size.");
- std::cout << "Pt " << iPoint << "Reg " << regime[iPoint] << std::endl;
- std::cout << "x " << mesh->getx(iPoint) << "xx " << mesh->getLoc(iPoint) << "xxExt " << mesh->getExt(iPoint) << std::endl;
+ std::cout << "Pt " << iPoint << "Reg " << name << " Turb " << regime[iPoint] << std::endl;
+ std::cout << std::setprecision(5) << "x " << x[iPoint] << ", xoc " << xoc[iPoint] <<
+ ", xx " << loc[iPoint] << "xxExt " << xxExt[iPoint] << std::endl;
std::cout << std::scientific << std::setprecision(15);
std::cout << std::setw(10) << "T " << u[iPoint * nVar + 0] << "\n"
<< std::setw(10) << "H " << u[iPoint * nVar + 1] << "\n"
@@ -189,8 +214,8 @@ void BoundaryLayer::printSolution(size_t iPoint) const
<< std::setw(10) << "ue " << u[iPoint * nVar + 3] << "\n"
<< std::setw(10) << "Ct " << u[iPoint * nVar + 4] << "\n"
<< std::setw(10) << "dStar (BL) " << deltaStar[iPoint] << "\n"
- << std::setw(10) << "dStar (old) " << blEdge->getDeltaStar(iPoint) << "\n"
- << std::setw(10) << "vt " << blEdge->getVt(iPoint) << "\n"
- << std::setw(10) << "Me " << blEdge->getMe(iPoint) << "\n"
- << std::setw(10) << "rhoe " << blEdge->getRhoe(iPoint) << std::endl;
+ << std::setw(10) << "dStar (ext) " << deltaStarExt[iPoint] << "\n"
+ << std::setw(10) << "vt " << vt[iPoint] << "\n"
+ << std::setw(10) << "Me " << Me[iPoint] << "\n"
+ << std::setw(10) << "rhoe " << rhoe[iPoint] << std::endl;
}
\ No newline at end of file
diff --git a/blast/src/DBoundaryLayer.h b/blast/src/DBoundaryLayer.h
index d56795e..8fbe4af 100644
--- a/blast/src/DBoundaryLayer.h
+++ b/blast/src/DBoundaryLayer.h
@@ -2,9 +2,8 @@
#define DBOUNDARYLAYER_H
#include "blast.h"
+#include "wObject.h"
#include "DClosures.h"
-#include "DDiscretization.h"
-#include "DEdge.h"
namespace blast
{
@@ -13,7 +12,7 @@ namespace blast
* @brief Boundary layer region upper/lower side or wake.
* @author Paul Dechamps
*/
-class BLAST_API BoundaryLayer
+class BLAST_API BoundaryLayer : public fwk::wSharedObject
{
private:
@@ -21,11 +20,21 @@ private:
double nCrit = 9.0; ///< Critical amplification factor.
public:
- Discretization *mesh; ///< 1D surface boundary layer mesh.
- Edge *blEdge; ///< Quantites at the inviscid boundary (edge of the boundary layer).
Closures *closSolver; ///< Closure relations class.
std::string name; ///< Name of the region.
+ // Mesh
+ size_t nNodes; ///< Number of points in the domain.
+ size_t nElms; ///< Number of cells in the domain.
+ std::vector<double> x; ///< x coordinate in the global frame of reference.
+ std::vector<double> y; ///< y coordinate in the global frame of reference.
+ std::vector<double> z; ///< z coordinate in the global frame of reference.
+ std::vector<double> xoc; ///< x/c coordinate in the global frame of reference.
+ std::vector<double> loc; ///< xi coordinate in the local frame of reference.
+ std::vector<double> dx; ///< Cell size.
+ std::vector<double> alpha; ///< Angle of the cell wrt the x axis of the global frame of reference.
+ double chord; ///< Chord of the section.
+
// Boundary layer variables.
std::vector<double> u; ///< Solution vector (theta, H, N, ue, Ct).
std::vector<double> Ret; ///< Reynolds number based on the momentum thickness (theta).
@@ -37,33 +46,58 @@ public:
std::vector<double> ctEq; ///< Equilibrium shear stress coefficient (turbulent BL).
std::vector<double> us; ///< Equivalent normalized wall slip velocity.
std::vector<double> delta; ///< Boundary layer thickness.
- std::vector<double> deltaStar; ///< Dispacement thickness (int(1-rho*u/rhoe*ue d_eta)).
+ std::vector<double> deltaStar; ///< Displacement thickness (int(1-rho*u/rhoe*ue d_eta)).
+
+ std::vector<double> Me; ///< Mach number.
+ std::vector<double> vt; ///< Velocity.
+ std::vector<double> rhoe; ///< Density.
+ std::vector<double> xxExt; ///< xi coordinate in the local frame of reference, at the edge of the boundary layer.
+ std::vector<double> deltaStarExt; ///< Displacement thickness.
+ std::vector<double> vtExt; ///< Previous velocity.
+
+ std::vector<double> blowingVelocity; ///< Blowing velocity.
// Transition related variables.
std::vector<int> regime; ///< Laminar (0) or turbulent (1) regime.
- std::vector<double> blowingVelocity; ///< Blowing velocity.
- size_t transMarker; ///< Marker id of the transition location.
- double xtrF; ///< Forced transition location in the global frame of reference.
- double xtr; ///< Transition location in the global frame of reference.
- double xoctr; ///< Transition location in % of the chord.
+ size_t transitionNode; ///< Node id of the transition location.
+ double xtrF; ///< Forced transition location in the local frame of reference (x/c).
+ double xtr; ///< Transition location in the local frame of reference (x/c).
- BoundaryLayer(double _xtrF = -1.0);
+ BoundaryLayer(double _xtrF = -1.0, std::string _name = "None");
~BoundaryLayer();
+ void reset();
+
+ // Getters
+ std::string getName() const { return name; };
+ size_t getnNodes() const { return nNodes; };
+ size_t getnElms() const { return nElms; };
+
+ std::vector<double> getDeltaStar() const { return deltaStar; };
+ std::vector<double> getUe() const {std::vector<double > ue(nNodes, 0.); for (size_t i = 0; i < nNodes; ++i) ue[i] = u[i * nVar + 3]; return ue;};
+ std::vector<double> getBlowing() const { return blowingVelocity; };
+
+ // Setters
+ void setMesh(std::vector<double> const _x, std::vector<double> const y, std::vector<double> const z, double const _chord, double const xMin);
+ void setMe(std::vector<double> const _Me) {if (_Me.size() != nNodes) throw std::runtime_error("Mach number vector is not consistent."); Me = _Me;};
+ void setVt(std::vector<double> const _vt) {if (_vt.size() != nNodes) throw std::runtime_error("Velocity vector is not consistent."); vt = _vt;};
+ void setRhoe(std::vector<double> const _rhoe) {if (_rhoe.size() != nNodes) throw std::runtime_error("Density vector is not consistent."); rhoe = _rhoe;};
+ void setxxExt(std::vector<double> const _xxExt) {if (_xxExt.size() != nNodes) throw std::runtime_error("External mesh vector is not consistent."); xxExt = _xxExt;};
+ void setDeltaStarExt(std::vector<double> const _deltaStarExt) {if (_deltaStarExt.size() != nNodes) throw std::runtime_error("Displacement thickness vector is not consistent."); deltaStarExt = _deltaStarExt;};
+ void setVtExt(std::vector<double> const _vtExt) {if (_vtExt.size() != nNodes) throw std::runtime_error("Velocity vector is not consistent."); vtExt = _vtExt;};
+
// Boundary conditions.
- void setStagBC(double Re);
+ void setStagnationBC(double Re);
void setWakeBC(double Re, std::vector<double> upperConditions, std::vector<double> lowerConditions);
- void setMesh(std::vector<double> x,
- std::vector<double> y,
- std::vector<double> z);
// Getters.
size_t getnVar() const { return nVar; };
double getnCrit() const { return nCrit; };
+ void setnCrit(double const _nCrit) {nCrit = _nCrit;};
// Others
void printSolution(size_t iPoint) const;
- void computeBlwVel();
+ void computeBlowingVelocity();
};
} // namespace blast
#endif // DBOUNDARYLAYER_H
diff --git a/blast/src/DDiscretization.cpp b/blast/src/DDiscretization.cpp
deleted file mode 100644
index 9a6796b..0000000
--- a/blast/src/DDiscretization.cpp
+++ /dev/null
@@ -1,79 +0,0 @@
-#include "DDiscretization.h"
-#include <cmath>
-#include <algorithm>
-
-using namespace blast;
-
-Discretization::Discretization()
-{
- nMarkers = 0;
- nElm = 0;
- prevnMarkers = 0;
-}
-
-Discretization::~Discretization() {}
-
-void Discretization::setGlob(std::vector<double> &_x, std::vector<double> &_y, std::vector<double> &_z)
-{
- 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.);
- y.resize(nMarkers, 0.);
- z.resize(nMarkers, 0.);
- nElm = nMarkers - 1;
- x = _x;
- y = _y;
- z = _z;
-
- computeBLcoord();
- computeOCcoord();
-}
-
-/**
- * @brief Compute nodes coordinates in the local frame of reference.
- */
-void Discretization::computeBLcoord()
-{
- loc.resize(nMarkers, 0.);
- alpha.resize(nElm, 0.);
- dx.resize(nElm, 0.);
-
- for (size_t iPoint = 0; iPoint < nElm; ++iPoint)
- {
- alpha[iPoint] = std::atan2(y[iPoint + 1] - y[iPoint], x[iPoint + 1] - x[iPoint]);
- // dx = sqrt(delta x ^2 + delta y ^2).
- dx[iPoint] = std::sqrt((x[iPoint + 1] - x[iPoint]) * (x[iPoint + 1] - x[iPoint]) + (y[iPoint + 1] - y[iPoint]) * (y[iPoint + 1] - y[iPoint]));
- loc[iPoint + 1] = loc[iPoint] + dx[iPoint];
- }
-}
-
-/**
- * @brief Compute x/chord coordinate.
-*/
-void Discretization::computeOCcoord()
-{
- xoc.resize(nMarkers, 0.);
-
- // Find min and max of array x
- auto min_it = std::min_element(x.begin(), x.end());
- double xMin = *min_it;
- auto max_it = std::max_element(x.begin(), x.end());
- double xMax = *max_it;
-
- // Compute xoc
- chord = xMax - xMin;
- if (chord <= 0)
- throw;
- for (size_t iPoint = 0; iPoint < x.size(); ++iPoint)
- xoc[iPoint] = (x[iPoint] - xMin) / chord;
-}
-
-void Discretization::setExt(std::vector<double> extM)
-{
- ext.resize(extM.size(), 0.);
- ext = extM;
-};
diff --git a/blast/src/DDiscretization.h b/blast/src/DDiscretization.h
deleted file mode 100644
index 1f4d534..0000000
--- a/blast/src/DDiscretization.h
+++ /dev/null
@@ -1,54 +0,0 @@
-#ifndef DDISCRETIZATION_H
-#define DDISCRETIZATION_H
-
-#include "blast.h"
-#include "DBoundaryLayer.h"
-
-namespace blast
-{
-
-/**
- * @brief Boundary layer mesh.
- * @author Paul Dechamps
- */
-class BLAST_API Discretization
-{
-private:
- std::vector<double> x; ///< x coordinate in the global frame of reference.
- std::vector<double> y; ///< y coordinate in the global frame of reference.
- std::vector<double> z; ///< z coordinate in the global frame of reference.
- std::vector<double> xoc; ///< x/c coordinate in the global frame of reference.
- std::vector<double> loc; ///< xi coordinate in the local frame of reference.
- std::vector<double> ext; ///< xi coordinate in the local frame of reference, at the edge of the boundary layer.
- std::vector<double> dx; ///< Cell size.
- std::vector<double> alpha; ///< Angle of the cell wrt the x axis of the global frame of reference.
- size_t nMarkers; ///< Number of points in the domain.
- size_t nElm; ///< Number of cells in the domain.
- double chord; ///< Chord of the section.
-
-public:
- size_t prevnMarkers; ///< Number of points in the domain at the previous iteration.
-
- Discretization();
- ~Discretization();
-
- // Getters.
- size_t getnMarkers() const { return nMarkers; };
- double getLoc(size_t iPoint) const { return loc[iPoint]; };
- double getx(size_t iPoint) const { return x[iPoint]; };
- double gety(size_t iPoint) const { return y[iPoint]; };
- double getz(size_t iPoint) const { return z[iPoint]; };
- double getxoc(size_t iPoint) const { return xoc[iPoint]; };
- double getExt(size_t iPoint) const { return ext[iPoint]; };
- double getAlpha(size_t iElm) const { return alpha[iElm]; };
- size_t getDiffnElm() const { return nMarkers - prevnMarkers; };
- double getChord() const { return chord; };
-
- // Setters.
- void setGlob(std::vector<double> &_x, std::vector<double> &_y, std::vector<double> &_z);
- void setExt(std::vector<double> extM);
- void computeBLcoord();
- void computeOCcoord();
-};
-} // namespace blast
-#endif // WDISCRETIZATION_H
diff --git a/blast/src/DDriver.cpp b/blast/src/DDriver.cpp
index a1fcba0..a03cf59 100644
--- a/blast/src/DDriver.cpp
+++ b/blast/src/DDriver.cpp
@@ -32,19 +32,10 @@ Driver::Driver(double _Re, double _Minf, double _CFL0, size_t nSections,
// Initialzing boundary layer
sections.resize(nSections);
+ solverExitCode = 0;
convergenceStatus.resize(nSections);
for (size_t iSec = 0; iSec < nSections; ++iSec)
- {
convergenceStatus[iSec].resize(3);
- for (size_t i = 0; i < 3; ++i)
- {
- BoundaryLayer *region = new BoundaryLayer(xtrF[i]);
- sections[iSec].push_back(region);
- }
- sections[iSec][0]->name = "upper";
- sections[iSec][1]->name = "lower";
- sections[iSec][2]->name = "wake";
- }
Cdt_sec.resize(sections.size(), 0.);
Cdf_sec.resize(sections.size(), 0.);
@@ -54,9 +45,6 @@ Driver::Driver(double _Re, double _Minf, double _CFL0, size_t nSections,
Cdp = 0.;
span = _span;
- // Pre/post processing flags.
- stagPtMvmt.resize(sections.size(), false);
-
// Time solver initialization.
tSolver = new Solver(_CFL0, _Minf, Re);
@@ -74,138 +62,89 @@ Driver::Driver(double _Re, double _Minf, double _CFL0, size_t nSections,
*/
Driver::~Driver()
{
- for (size_t iSec = 0; iSec < sections.size(); ++iSec)
- for (size_t i = 0; i < sections[iSec].size(); ++i)
- delete sections[iSec][i];
delete tSolver;
std::cout << "~blast::Driver()\n";
} // end ~Driver
+void Driver::reset()
+{
+ size_t nSections = sections.size();
+ for (size_t iSec = 0; iSec < sections.size(); ++iSec)
+ sections[iSec].resize(0);
+ sections.resize(nSections);
+}
+
/**
* @brief Runs viscous operations. Temporal solver parametrisation is first adapted,
* Solutions are initialized than time marched towards steady state successively for each points.
*
- * @param couplIter Current coupling iteration.
* @return int Solver exit code.
*/
-int Driver::run(unsigned int const couplIter)
+int Driver::run()
{
- bool lockTrans = false; // Flagging activation of transition routines.
+ bool lockTrans = false; // Flag activation of transition routines.
int pointExitCode; // Output of pseudo time integration (0: converged; -1: unsuccessful, failed; >0: unsuccessful nb iter).
double maxMach = 0.;
for (size_t iSec = 0; iSec < sections.size(); ++iSec)
{
- tms["0-CheckIn"].start();
if (verbose > 1)
std::cout << "Sec " << iSec << " ";
- // Check for stagnation point movement.
- if (couplIter != 0 && (sections[iSec][0]->mesh->getDiffnElm()))
- {
- stagPtMvmt[iSec] = true;
- if (verbose > 2)
- std::cout << "Stagnation point moved by " << sections[iSec][0]->mesh->getDiffnElm() << " elements." << std::endl;
- }
- else
- stagPtMvmt[iSec] = false;
-
- // Set boundary conditions.
- sections[iSec][0]->xtr = 1.; // Upper side initial xtr.
- sections[iSec][1]->xtr = 1.; // Lower side initial xtr.
- sections[iSec][2]->xtr = 0.; // Wake initial xtr (fully turbulent).
- sections[iSec][0]->setStagBC(Re);
- sections[iSec][1]->setStagBC(Re);
- tms["0-CheckIn"].stop();
-
// Loop over the regions (upper, lower and wake).
- for (size_t iRegion = 0; iRegion < sections[iSec].size(); ++iRegion)
+ int iRegion = 0;
+ for (auto reg: sections[iSec])
{
convergenceStatus[iSec][iRegion].resize(0);
- // Maximum Mach number in the region.
- if (sections[iSec][iRegion]->blEdge->getMaxMach() > maxMach)
- maxMach = sections[iSec][iRegion]->blEdge->getMaxMach();
-
- if (verbose > 1)
- std::cout << sections[iSec][iRegion]->name << " ";
-
- // Check if safe mode is required.
- if (couplIter == 0 || sections[iSec][iRegion]->mesh->getDiffnElm() != 0 || stagPtMvmt[iSec] == true || solverExitCode != 0)
- {
- tSolver->setCFL0(1);
- tSolver->setitFrozenJac(1);
- tSolver->setinitSln(true);
-
- if (sections[iSec][iRegion]->mesh->getDiffnElm())
- {
- std::vector<double> xxExtCurr(sections[iSec][iRegion]->mesh->getnMarkers(), 0.);
- for (size_t iPoint = 0; iPoint < xxExtCurr.size(); ++iPoint)
- xxExtCurr[iPoint] = sections[iSec][iRegion]->mesh->getLoc(iPoint);
- sections[iSec][iRegion]->mesh->setExt(xxExtCurr);
-
- std::vector<double> DeltaStarZeros(sections[iSec][iRegion]->mesh->getnMarkers(), 0.0);
- sections[iSec][iRegion]->blEdge->setDeltaStar(DeltaStarZeros);
-
- std::vector<double> ueOld(sections[iSec][iRegion]->mesh->getnMarkers(), 0.0);
- for (size_t iPoint = 0; iPoint < ueOld.size(); ++iPoint)
- ueOld[iPoint] = sections[iSec][iRegion]->blEdge->getVt(iPoint);
- sections[iSec][iRegion]->blEdge->setVtOld(ueOld);
- }
-
- // Keyword annoncing iteration safety level.
- if (verbose > 1)
- std::cout << "restart. ";
- }
+ // Reset transition
+ if (reg->name == "wake")
+ reg->xtr = 0.;
+ else if (reg-> name == "upper" || reg->name == "lower")
+ reg->xtr = 1.;
else
{
- tSolver->setCFL0(CFL0);
- tSolver->setitFrozenJac(5);
- tSolver->setinitSln(false);
-
- // Keyword annoncing iteration safety level.
- if (verbose > 1)
- std::cout << "update. ";
+ std::cout << "reg->name " << reg->name << std::endl;
+ throw std::runtime_error("Wrong region name\n");
}
- // Save number of points in each regions.
- sections[iSec][iRegion]->mesh->prevnMarkers = sections[iSec][iRegion]->mesh->getnMarkers();
-
// Reset regime.
lockTrans = false;
- if (iRegion < 2) // Airfoil
- for (size_t k = 0; k < sections[iSec][iRegion]->mesh->getnMarkers(); k++)
- sections[iSec][iRegion]->regime[k] = 0;
- else if (iRegion == 2) // Wake
+ if (reg->name == "upper" || reg->name == "lower")
+ {
+ for (size_t k = 0; k < reg->nNodes; k++)
+ reg->regime[k] = 0;
+ reg->setStagnationBC(Re);
+ }
+
+ else if (reg->name == "wake") // Wake
{
- for (size_t k = 0; k < sections[iSec][iRegion]->mesh->getnMarkers(); k++)
- sections[iSec][iRegion]->regime[k] = 1;
+ for (size_t k = 0; k < reg->nNodes; k++)
+ reg->regime[k] = 1;
// Impose wake boundary condition.
- std::vector<double> upperConditions(sections[iSec][iRegion]->getnVar(), 0.);
- std::vector<double> lowerConditions(sections[iSec][iRegion]->getnVar(), 0.);
+ std::vector<double> upperConditions(reg->getnVar(), 0.);
+ std::vector<double> lowerConditions(reg->getnVar(), 0.);
for (size_t k = 0; k < upperConditions.size(); ++k)
{
- upperConditions[k] = sections[iSec][0]->u[(sections[iSec][0]->mesh->getnMarkers() - 1) * sections[iSec][0]->getnVar() + k];
- lowerConditions[k] = sections[iSec][1]->u[(sections[iSec][1]->mesh->getnMarkers() - 1) * sections[iSec][1]->getnVar() + k];
+ upperConditions[k] = sections[iSec][0]->u[(sections[iSec][0]->nNodes - 1) * sections[iSec][0]->getnVar() + k];
+ lowerConditions[k] = sections[iSec][1]->u[(sections[iSec][1]->nNodes - 1) * sections[iSec][1]->getnVar() + k];
}
- sections[iSec][iRegion]->setWakeBC(Re, upperConditions, lowerConditions);
+ reg->setWakeBC(Re, upperConditions, lowerConditions);
lockTrans = true;
}
else
- throw;
+ throw std::runtime_error("Wrong region name\n");
// Loop over points
- for (size_t iPoint = 1; iPoint < sections[iSec][iRegion]->mesh->getnMarkers(); ++iPoint)
+ for (size_t iPoint = 1; iPoint < reg->nNodes; ++iPoint)
{
// Initialize solution at point
- tSolver->initialCondition(iPoint, sections[iSec][iRegion]);
+ tSolver->initialCondition(iPoint, reg);
// Solve equations
tms["1-Solver"].start();
- pointExitCode = tSolver->integration(iPoint, sections[iSec][iRegion]);
- if (sections[iSec][iRegion]->xtrF != -1 && sections[iSec][iRegion]->mesh->getxoc(iPoint) >= sections[iSec][iRegion]->xtrF)
- sections[iSec][iRegion]->u[iPoint * sections[iSec][iRegion]->getnVar()+2] = 9.;
+ pointExitCode = tSolver->integration(iPoint, reg);
tms["1-Solver"].stop();
if (pointExitCode != 0)
@@ -215,29 +154,47 @@ int Driver::run(unsigned int const couplIter)
tms["2-Transition"].start();
if (!lockTrans)
{
- // Check if perturbation waves are growing or not.
- if (sections[iSec][iRegion]->xtrF != -1 && sections[iSec][iRegion]->mesh->getxoc(iPoint) <= sections[iSec][iRegion]->xtrF)
- checkWaves(iPoint, sections[iSec][iRegion]);
-
- // Amplification factor is compared to critical amplification factor.
- if (sections[iSec][iRegion]->u[iPoint * sections[iSec][iRegion]->getnVar() + 2] >= sections[iSec][iRegion]->getnCrit())
+ // Forced transition
+ if (reg->xtrF != -1 && reg->xoc[iPoint] >= reg->xtrF)
{
- averageTransition(iPoint, sections[iSec][iRegion], 0);
+ reg->u[iPoint * reg->getnVar() + 2] = reg->getnCrit();
+ averageTransition(iPoint, reg, true);
if (verbose > 1)
{
std::cout << std::fixed;
std::cout << std::setprecision(2);
- std::cout << sections[iSec][iRegion]->xoctr
- << " (" << sections[iSec][iRegion]->transMarker << ") ";
+ std::cout << reg->xtr
+ << " (" << reg->transitionNode << ") ";
}
lockTrans = true;
}
+ // Free transtion
+ else if (reg->xtrF == -1)
+ {
+ // Check if perturbation waves are growing or not.
+ checkWaves(iPoint, reg);
+
+ // Amplification factor is compared to critical amplification factor.
+ if (reg->u[iPoint * reg->getnVar() + 2] >= reg->getnCrit())
+ {
+ averageTransition(iPoint, reg, false);
+ if (verbose > 1)
+ {
+ std::cout << std::fixed;
+ std::cout << std::setprecision(2);
+ std::cout << reg->xtr
+ << " (" << reg->transitionNode << ") ";
+ }
+ lockTrans = true;
+ }
+ }
}
tms["2-Transition"].stop();
}
tms["3-Blowing"].start();
- sections[iSec][iRegion]->computeBlwVel();
+ reg->computeBlowingVelocity();
tms["3-Blowing"].stop();
+ iRegion++;
}
if (verbose > 1)
std::cout << std::endl;
@@ -280,7 +237,7 @@ void Driver::checkWaves(size_t iPoint, BoundaryLayer *bl)
* @param bl BoundaryLayer region.
* @param forced Flag for forced transition.
*/
-void Driver::averageTransition(size_t iPoint, BoundaryLayer *bl, int forced)
+void Driver::averageTransition(size_t iPoint, BoundaryLayer *bl, bool forced)
{
// Averages solution on transition marker.
size_t nVar = bl->getnVar();
@@ -291,23 +248,37 @@ void Driver::averageTransition(size_t iPoint, BoundaryLayer *bl, int forced)
lamSol[k] = bl->u[iPoint * nVar + k];
// Set up turbulent portion boundary condition.
- bl->transMarker = iPoint; // Save transition marker.
+ bl->transitionNode = iPoint; // Save transition marker.
// Loop over remaining points in the region.
- for (size_t k = iPoint; k < bl->mesh->getnMarkers(); ++k)
+ for (size_t k = iPoint; k < bl->nNodes; ++k)
bl->regime[k] = 1; // Set Turbulent regime.
// Compute transition location.
- bl->xtr = (bl->getnCrit() - bl->u[(iPoint - 1) * nVar + 2]) * ((bl->mesh->getx(iPoint) - bl->mesh->getx(iPoint - 1)) / (bl->u[iPoint * nVar + 2] - bl->u[(iPoint - 1) * nVar + 2])) + bl->mesh->getx(iPoint - 1);
- bl->xoctr = (bl->getnCrit() - bl->u[(iPoint - 1) * nVar + 2]) * ((bl->mesh->getxoc(iPoint) - bl->mesh->getxoc(iPoint - 1)) / (bl->u[iPoint * nVar + 2] - bl->u[(iPoint - 1) * nVar + 2])) + bl->mesh->getxoc(iPoint - 1);
+ // if (forced)
+ // bl->xtr = bl->xtrF;
+ // else
+ bl->xtr = (bl->getnCrit() - bl->u[(iPoint - 1) * nVar + 2]) * ((bl->x[iPoint] - bl->x[iPoint-1]) / (bl->u[iPoint * nVar + 2] - bl->u[(iPoint - 1) * nVar + 2])) + bl->x[iPoint-1];
- // Percentage of the subsection corresponding to a turbulent flow.
- double avTurb = (bl->mesh->getx(iPoint) - bl->xtr) / (bl->mesh->getx(iPoint) - bl->mesh->getx(iPoint - 1));
+ double avTurb = (bl->x[iPoint] - bl->xtr) / (bl->x[iPoint] - bl->x[iPoint-1]);
double avLam = 1. - avTurb;
+ //std::cout << "avTurb " << avTurb << " avLam " << avLam << std::endl;
+ if (avTurb < 0. || avTurb > 1. || avLam < 0. || avLam > 1.)
+ {
+ bl->printSolution(iPoint);
+ bl->printSolution(iPoint - 1);
+ std::cout << "dx " << std::abs(bl->xoc[iPoint] - bl->xoc[iPoint - 1]) << std::endl;
+ std::cout << " (bl->u[iPoint * nVar + 2] - bl->u[(iPoint - 1) * nVar + 2]) " << (bl->u[iPoint * nVar + 2] - bl->u[(iPoint - 1) * nVar + 2]) << std::endl;
+ std::cout << " bl->getnCrit() - bl->u[(iPoint - 1) * nVar + 2] " << bl->getnCrit() - bl->u[(iPoint - 1) * nVar + 2] << std::endl;
+ std::cout << "xtr " << bl->xtr << std::endl;
+ std::cout << "avTurb " << avTurb << " avLam " << avLam << std::endl;
+ throw std::runtime_error("Transition location out of bounds.");
+ }
+
// Impose boundary condition.
double Cteq_trans;
- bl->closSolver->turbulentClosures(Cteq_trans, bl->u[(iPoint - 1) * nVar + 0], bl->u[(iPoint - 1) * nVar + 1], 0.03, bl->Ret[iPoint - 1], bl->blEdge->getMe(iPoint - 1), bl->name);
+ bl->closSolver->turbulentClosures(Cteq_trans, bl->u[(iPoint - 1) * nVar + 0], bl->u[(iPoint - 1) * nVar + 1], 0.03, bl->Ret[iPoint - 1], bl->Me[iPoint-1], bl->name);
bl->ctEq[iPoint - 1] = Cteq_trans;
bl->u[(iPoint - 1) * nVar + 4] = 0.7 * bl->ctEq[iPoint - 1];
@@ -331,7 +302,7 @@ void Driver::averageTransition(size_t iPoint, BoundaryLayer *bl, int forced)
// Recompute closures. (The turbulent BC @ iPoint - 1 does not influence laminar closure relations).
std::vector<double> lamParam(8, 0.);
- bl->closSolver->laminarClosures(lamParam, bl->u[(iPoint - 1) * nVar + 0], bl->u[(iPoint - 1) * nVar + 1], bl->Ret[iPoint - 1], bl->blEdge->getMe(iPoint - 1), bl->name);
+ bl->closSolver->laminarClosures(lamParam, bl->u[(iPoint - 1) * nVar + 0], bl->u[(iPoint - 1) * nVar + 1], bl->Ret[iPoint - 1], bl->Me[iPoint-1], bl->name);
bl->Hstar[iPoint - 1] = lamParam[0];
bl->Hstar2[iPoint - 1] = lamParam[1];
bl->Hk[iPoint - 1] = lamParam[2];
@@ -342,7 +313,7 @@ void Driver::averageTransition(size_t iPoint, BoundaryLayer *bl, int forced)
bl->us[iPoint - 1] = lamParam[7];
std::vector<double> turbParam(8, 0.);
- bl->closSolver->turbulentClosures(turbParam, bl->u[(iPoint)*nVar + 0], bl->u[(iPoint)*nVar + 1], bl->u[(iPoint)*nVar + 4], bl->Ret[iPoint], bl->blEdge->getMe(iPoint), bl->name);
+ bl->closSolver->turbulentClosures(turbParam, bl->u[(iPoint)*nVar + 0], bl->u[(iPoint)*nVar + 1], bl->u[(iPoint)*nVar + 4], bl->Ret[iPoint], bl->Me[iPoint], bl->name);
bl->Hstar[iPoint] = turbParam[0];
bl->Hstar2[iPoint] = turbParam[1];
bl->Hk[iPoint] = turbParam[2];
@@ -363,41 +334,41 @@ void Driver::averageTransition(size_t iPoint, BoundaryLayer *bl, int forced)
* @param bl BoundaryLayer region.
* @param i Considered section.
*/
-void Driver::computeSectionalDrag(std::vector<BoundaryLayer *> const &bl, size_t i)
+void Driver::computeSectionalDrag(std::vector<BoundaryLayer *> const &sec, size_t i)
{
- size_t nVar = bl[0]->getnVar();
- size_t lastWkPt = (bl[2]->mesh->getnMarkers() - 1) * nVar;
+ size_t nVar = sec[0]->getnVar();
+ size_t lastWkPt = (sec[2]->nNodes - 1) * nVar;
// Normalize friction and dissipation coefficients.
- std::vector<double> frictioncoeff_upper(bl[0]->mesh->getnMarkers(), 0.0);
+ std::vector<double> frictioncoeff_upper(sec[0]->nNodes, 0.0);
for (size_t k = 0; k < frictioncoeff_upper.size(); ++k)
- frictioncoeff_upper[k] = bl[0]->blEdge->getVt(k) * bl[0]->blEdge->getVt(k) * bl[0]->cf[k];
+ frictioncoeff_upper[k] = sec[0]->vt[k] * sec[0]->vt[k] * sec[0]->cf[k];
- std::vector<double> frictioncoeff_lower(bl[1]->mesh->getnMarkers(), 0.0);
+ std::vector<double> frictioncoeff_lower(sec[1]->nNodes, 0.0);
for (size_t k = 0; k < frictioncoeff_lower.size(); ++k)
- frictioncoeff_lower[k] = bl[1]->blEdge->getVt(k) * bl[1]->blEdge->getVt(k) * bl[1]->cf[k];
+ frictioncoeff_lower[k] = sec[1]->vt[k] * sec[1]->vt[k] * sec[1]->cf[k];
// Compute friction drag coefficient.
double Cdf_upper = 0.0;
- for (size_t k = 1; k < bl[0]->mesh->getnMarkers(); ++k)
- Cdf_upper += (frictioncoeff_upper[k] + frictioncoeff_upper[k - 1]) * (bl[0]->mesh->getx(k) - bl[0]->mesh->getx(k - 1)) / 2;
+ for (size_t k = 1; k < sec[0]->nNodes; ++k)
+ Cdf_upper += (frictioncoeff_upper[k] + frictioncoeff_upper[k - 1]) * (sec[0]->x[k] - sec[0]->x[k-1]) / 2;
double Cdf_lower = 0.0;
- for (size_t k = 1; k < bl[1]->mesh->getnMarkers(); ++k)
- Cdf_lower += (frictioncoeff_lower[k] + frictioncoeff_lower[k - 1]) * (bl[1]->mesh->getx(k) - bl[1]->mesh->getx(k - 1)) / 2;
+ for (size_t k = 1; k < sec[1]->nNodes; ++k)
+ Cdf_lower += (frictioncoeff_lower[k] + frictioncoeff_lower[k - 1]) * (sec[1]->x[k] - sec[1]->x[k-1]) / 2;
// Friction drag coefficient.
Cdf_sec[i] = Cdf_upper + Cdf_lower; // No friction in the wake
// Total drag coefficient.
- Cdt_sec[i] = 2 * bl[2]->u[lastWkPt + 0] * pow(bl[2]->u[lastWkPt + 3], (bl[2]->u[lastWkPt + 1] + 5) / 2);
+ Cdt_sec[i] = 2 * sec[2]->u[lastWkPt + 0] * pow(sec[2]->u[lastWkPt + 3], (sec[2]->u[lastWkPt + 1] + 5) / 2);
// Pressure drag coefficient.
Cdp_sec[i] = Cdt_sec[i] - Cdf_sec[i];
}
-double Driver::getAvrgxoctr(size_t const iRegion) const
+double Driver::getAverageTransition(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[iSec][iRegion]->xtr;
xtr /= sections.size();
return xtr;
}
@@ -421,13 +392,13 @@ void Driver::computeTotalDrag()
}
else
{
- Cdt += (sections[0][0]->mesh->getz(0) - 0) * Cdt_sec[0];
- Cdp += (sections[0][0]->mesh->getz(0) - 0) * Cdp_sec[0];
- Cdf += (sections[0][0]->mesh->getz(0) - 0) * Cdf_sec[0];
+ Cdt += (sections[0][0]->z[0] - 0) * Cdt_sec[0];
+ Cdp += (sections[0][0]->z[0] - 0) * Cdp_sec[0];
+ Cdf += (sections[0][0]->z[0] - 0) * Cdf_sec[0];
double dz = 0.;
for (size_t k = 1; k < sections.size(); ++k)
{
- dz = (sections[k][0]->mesh->getz(0) - sections[k - 1][0]->mesh->getz(0));
+ dz = (sections[k][0]->z[0] - sections[k - 1][0]->z[0]);
Cdt += dz * Cdt_sec[k];
Cdp += dz * Cdp_sec[k];
Cdf += dz * Cdf_sec[k];
@@ -472,7 +443,7 @@ int Driver::outputStatus(double maxMach) const
if (convergenceStatus[iSec][iReg].size() != 0)
for (size_t iPoint = 0; iPoint < convergenceStatus[iSec][iReg].size(); ++iPoint)
{
- std::cout << convergenceStatus[iSec][iReg][iPoint] << " (" << sections[iSec][iReg]->mesh->getx(convergenceStatus[iSec][iReg][iPoint]) << "), ";
+ std::cout << convergenceStatus[iSec][iReg][iPoint] << " (" << sections[iSec][iReg]->x[convergenceStatus[iSec][iReg][iPoint]] << "), ";
++count;
}
std::cout << "\n";
@@ -488,7 +459,7 @@ int Driver::outputStatus(double maxMach) const
std::vector<double> meanTransitions(2, 0.);
for (size_t iSec = 0; iSec < sections.size(); ++iSec)
for (size_t iReg = 0; iReg < 2; ++iReg)
- meanTransitions[iReg] += sections[iSec][iReg]->xoctr;
+ meanTransitions[iReg] += sections[iSec][iReg]->xtr;
for (size_t iReg = 0; iReg < meanTransitions.size(); ++iReg)
meanTransitions[iReg] /= sections.size();
@@ -526,184 +497,6 @@ int Driver::outputStatus(double maxMach) const
return solverExitCode;
}
-/**
- * @brief Set nodes coordinates in the global frame of reference.
- *
- * @param iSec id of the considered section.
- * @param iRegion id of the considered region.
- * @param _x Vector of x coordinates of the nodes in the region.
- * @param _y Vector of y coordinates of the nodes in the region.
- * @param _z Vector of z coordinates of the nodes in the region.
- */
-void Driver::setGlobMesh(size_t iSec, size_t iRegion, std::vector<double> _x, std::vector<double> _y, std::vector<double> _z)
-{
- if (iSec > sections.size() || iRegion > sections[iSec].size())
- throw std::runtime_error("blast::Driver Wrong section or region id.");
- sections[iSec][iRegion]->setMesh(_x, _y, _z);
-}
-
-/**
- * @brief Set the velocity at the nodes in the global frame of reference.
- *
- * @param iSec id of the considered section.
- * @param iRegion id of the considered region.
- * @param vx Vector of velocities projected on the x axis of the global frame of reference.
- * @param vy Vector of velocities projected on the y axis of the global frame of reference.
- * @param vz Vector of velocities projected on the z axis of the global frame of reference.
- */
-void Driver::setVelocities(size_t iSec, size_t iRegion, std::vector<double> vx, std::vector<double> vy, std::vector<double> vz)
-{
- if (vx.size() != vy.size() || vy.size() != vz.size() || vx.size() != vz.size())
- throw std::runtime_error("blast::Driver Wrong velocity vector sizes.");
- if (vx.size() != sections[iSec][iRegion]->mesh->getnMarkers())
- throw std::runtime_error("blast::Driver Velocity vector is not coherent with the mesh.");
-
- std::vector<double> vt(vx.size(), 0.0);
-
- for (size_t iPoint = 0; iPoint < vx.size() - 1; ++iPoint)
- {
- vt[iPoint] = vx[iPoint] * std::cos(sections[iSec][iRegion]->mesh->getAlpha(iPoint)) + vy[iPoint] * std::sin(sections[iSec][iRegion]->mesh->getAlpha(iPoint));
- vt[iPoint + 1] = vx[iPoint + 1] * std::cos(sections[iSec][iRegion]->mesh->getAlpha(iPoint)) + vy[iPoint + 1] * std::sin(sections[iSec][iRegion]->mesh->getAlpha(iPoint));
- }
-
- sections[iSec][iRegion]->blEdge->setVt(vt);
-}
-
-/**
- * @brief Set the Mach number at the nodes.
- *
- * @param iSec id of the considered section.
- * @param iRegion id of the considered region.
- * @param _Me Vector of Mach numbers at the nodes.
- */
-void Driver::setMe(size_t iSec, size_t iRegion, std::vector<double> _Me)
-{
- if (_Me.size() != sections[iSec][iRegion]->mesh->getnMarkers())
- throw std::runtime_error("blast::Driver Mach number vector is not coherent with the mesh.");
- sections[iSec][iRegion]->blEdge->setMe(_Me);
-}
-
-/**
- * @brief Set the density at the nodes.
- *
- * @param iSec id of the considered section.
- * @param iRegion id of the considered region.
- * @param _Rhoe Vector of density at the nodes.
- */
-void Driver::setRhoe(size_t iSec, size_t iRegion, std::vector<double> _Rhoe)
-{
- if (_Rhoe.size() != sections[iSec][iRegion]->mesh->getnMarkers())
- throw std::runtime_error("blast::Driver Density vector is not coherent with the mesh.");
- sections[iSec][iRegion]->blEdge->setRhoe(_Rhoe);
-}
-
-/**
- * @brief Set the nodes coordinates in the local frame of reference at the edge of the boundary layer.
- *
- * @param iSec id of the considered section.
- * @param iRegion id of the considered region.
- * @param _xxExt Vector of coordinates of the nodes.
- */
-void Driver::setxxExt(size_t iSec, size_t iRegion, std::vector<double> _xxExt)
-{
- if (_xxExt.size() != sections[iSec][iRegion]->mesh->getnMarkers())
- throw std::runtime_error("blast::Driver External mesh vector is not coherent with the mesh.");
- sections[iSec][iRegion]->mesh->setExt(_xxExt);
-}
-
-/**
- * @brief Set the displacement thickness (previous iteration).
- *
- * @param iSec id of the considered section.
- * @param iRegion id of the considered region.
- * @param _DeltaStarExt Vector of displacement thicknesses at the nodes.
- */
-void Driver::setDeltaStarExt(size_t iSec, size_t iRegion, std::vector<double> _DeltaStarExt)
-{
- if (_DeltaStarExt.size() != sections[iSec][iRegion]->mesh->getnMarkers())
- {
- std::cout << "size DeltaStar: " << _DeltaStarExt.size() << ". nMarkers: " << sections[iSec][iRegion]->mesh->getnMarkers() << std::endl;
- throw std::runtime_error("blast::Driver External delta star vector is not coherent with the mesh.");
- }
- sections[iSec][iRegion]->blEdge->setDeltaStar(_DeltaStarExt);
-}
-
-/**
- * @brief Returns x coordinates of the nodes in the local frame of reference in the considered region.
- *
- * @param iSec id of the considered section.
- * @param iRegion id of the considered region.
- * @return std::vector<double>
- */
-std::vector<double> Driver::extractxCoord(size_t iSec, size_t iRegion) const
-{
- std::vector<double> xCoord(sections[iSec][iRegion]->mesh->getnMarkers(), 0.);
- for (size_t iPoint = 0; iPoint < xCoord.size(); ++iPoint)
- xCoord[iPoint] = sections[iSec][iRegion]->mesh->getx(iPoint);
- return xCoord;
-}
-
-/**
- * @brief Returns blowing velocities (defined on the cells' cg) in the considered region.
- *
- * @param iSec id of the considered section.
- * @param iRegion id of the considered region.
- * @return std::vector<double>
- */
-std::vector<double> Driver::extractBlowingVel(size_t iSec, size_t iRegion) const
-{
- return sections[iSec][iRegion]->blowingVelocity;
-}
-std::vector<double> Driver::extractxx(size_t iSec, size_t iRegion) const
-{
- std::vector<double> xx(sections[iSec][iRegion]->mesh->getnMarkers(), 0.);
- for (size_t iPoint = 0; iPoint < xx.size(); ++iPoint)
- xx[iPoint] = sections[iSec][iRegion]->mesh->getLoc(iPoint);
- return xx;
-}
-
-/**
- * @brief Returns the displacement thickness in the considered region.
- *
- * @param iSec id of the considered section.
- * @param iRegion id of the considered region.
- * @return std::vector<double>
- */
-std::vector<double> Driver::extractDeltaStar(size_t iSec, size_t iRegion) const
-{
- return sections[iSec][iRegion]->deltaStar;
-}
-
-/**
- * @brief Returns the velocity of the interaction law in the considered region.
- *
- * @param iSec id of the considered section.
- * @param iRegion id of the considered region.
- * @return std::vector<double>
- */
-std::vector<double> Driver::extractUe(size_t iSec, size_t iRegion) const
-{
- std::vector<double> ueCurr(sections[iSec][iRegion]->mesh->getnMarkers(), 0.0);
- for (size_t i = 0; i < ueCurr.size(); ++i)
- ueCurr[i] = sections[iSec][iRegion]->u[i * sections[iSec][iRegion]->getnVar() + 3];
- return ueCurr;
-}
-
-/**
- * @brief Set the velocity at the edge of the BL at the previous iteration in the considered region.
- *
- * @param iSec id of the considered section.
- * @param iRegion id of the considered region.
- * @param _ue Velocity vector.
- * @return std::vector<double>
- */
-void Driver::setUeOld(size_t iSec, size_t iRegion, std::vector<double> _ue)
-{
- if (_ue.size() != sections[iSec][iRegion]->mesh->getnMarkers())
- throw std::runtime_error("blast::Driver External velocity vector is not coherent with the mesh.");
- sections[iSec][iRegion]->blEdge->setVtOld(_ue);
-}
-
/**
* @brief Returns the integral boundary layer solution in a section (sorted from upper Te to lower Te than wake).
*
@@ -712,13 +505,13 @@ void Driver::setUeOld(size_t iSec, size_t iRegion, std::vector<double> _ue)
*/
std::vector<std::vector<double>> Driver::getSolution(size_t iSec)
{
- size_t nMarkersUp = sections[iSec][0]->mesh->getnMarkers();
- size_t nMarkersLw = sections[iSec][1]->mesh->getnMarkers(); // No stagnation point doublon.
+ size_t nMarkersUp = sections[iSec][0]->nNodes;
+ size_t nMarkersLw = sections[iSec][1]->nNodes; // No stagnation point doublon.
size_t nVar = sections[iSec][0]->getnVar();
std::vector<std::vector<double>> Sln(20);
for (size_t iVector = 0; iVector < Sln.size(); ++iVector)
- Sln[iVector].resize(nMarkersUp + nMarkersLw + sections[iSec][2]->mesh->getnMarkers(), 0.);
+ Sln[iVector].resize(nMarkersUp + nMarkersLw + sections[iSec][2]->nNodes, 0.);
// Upper side.
size_t revPoint = nMarkersUp - 1;
for (size_t iPoint = 0; iPoint < nMarkersUp; ++iPoint)
@@ -740,15 +533,15 @@ std::vector<std::vector<double>> Driver::getSolution(size_t iSec)
Sln[13][iPoint] = sections[iSec][0]->us[revPoint];
Sln[14][iPoint] = sections[iSec][0]->delta[revPoint];
- Sln[15][iPoint] = sections[iSec][0]->mesh->getx(revPoint);
- Sln[16][iPoint] = sections[iSec][0]->mesh->getxoc(revPoint);
- Sln[17][iPoint] = sections[iSec][0]->mesh->gety(revPoint);
- Sln[18][iPoint] = sections[iSec][0]->mesh->getz(revPoint);
- Sln[19][iPoint] = sections[iSec][0]->blEdge->getVt(revPoint);
+ Sln[15][iPoint] = sections[iSec][0]->x[revPoint];
+ Sln[16][iPoint] = sections[iSec][0]->xoc[revPoint];
+ Sln[17][iPoint] = sections[iSec][0]->y[revPoint];
+ Sln[18][iPoint] = sections[iSec][0]->z[revPoint];
+ Sln[19][iPoint] = sections[iSec][0]->vt[revPoint];
--revPoint;
}
// Lower side.
- for (size_t iPoint = 0; iPoint < sections[iSec][1]->mesh->getnMarkers(); ++iPoint)
+ for (size_t iPoint = 0; iPoint < sections[iSec][1]->nNodes; ++iPoint)
{
Sln[0][nMarkersUp + iPoint] = sections[iSec][1]->u[iPoint * nVar + 0];
Sln[1][nMarkersUp + iPoint] = sections[iSec][1]->u[iPoint * nVar + 1];
@@ -767,14 +560,14 @@ std::vector<std::vector<double>> Driver::getSolution(size_t iSec)
Sln[13][nMarkersUp + iPoint] = sections[iSec][1]->us[iPoint];
Sln[14][nMarkersUp + iPoint] = sections[iSec][1]->delta[iPoint];
- Sln[15][nMarkersUp + iPoint] = sections[iSec][1]->mesh->getx(iPoint);
- Sln[16][nMarkersUp + iPoint] = sections[iSec][1]->mesh->getxoc(iPoint);
- Sln[17][nMarkersUp + iPoint] = sections[iSec][1]->mesh->gety(iPoint);
- Sln[18][nMarkersUp + iPoint] = sections[iSec][1]->mesh->getz(iPoint);
- Sln[19][nMarkersUp + iPoint] = sections[iSec][1]->blEdge->getVt(iPoint);
+ Sln[15][nMarkersUp + iPoint] = sections[iSec][1]->x[iPoint];
+ Sln[16][nMarkersUp + iPoint] = sections[iSec][1]->xoc[iPoint];
+ Sln[17][nMarkersUp + iPoint] = sections[iSec][1]->y[iPoint];
+ Sln[18][nMarkersUp + iPoint] = sections[iSec][1]->z[iPoint];
+ Sln[19][nMarkersUp + iPoint] = sections[iSec][1]->vt[iPoint];
}
// Wake.
- for (size_t iPoint = 0; iPoint < sections[iSec][2]->mesh->getnMarkers(); ++iPoint)
+ for (size_t iPoint = 0; iPoint < sections[iSec][2]->nNodes; ++iPoint)
{
Sln[0][nMarkersUp + nMarkersLw + iPoint] = sections[iSec][2]->u[iPoint * nVar + 0];
Sln[1][nMarkersUp + nMarkersLw + iPoint] = sections[iSec][2]->u[iPoint * nVar + 1];
@@ -793,11 +586,11 @@ std::vector<std::vector<double>> Driver::getSolution(size_t iSec)
Sln[13][nMarkersUp + nMarkersLw + iPoint] = sections[iSec][2]->us[iPoint];
Sln[14][nMarkersUp + nMarkersLw + iPoint] = sections[iSec][2]->delta[iPoint];
- Sln[15][nMarkersUp + nMarkersLw + iPoint] = sections[iSec][2]->mesh->getx(iPoint);
- Sln[16][nMarkersUp + nMarkersLw + iPoint] = sections[iSec][2]->mesh->getxoc(iPoint);
- Sln[17][nMarkersUp + nMarkersLw + iPoint] = sections[iSec][2]->mesh->gety(iPoint);
- Sln[18][nMarkersUp + nMarkersLw + iPoint] = sections[iSec][2]->mesh->getz(iPoint);
- Sln[19][nMarkersUp + nMarkersLw + iPoint] = sections[iSec][2]->blEdge->getVt(iPoint);
+ Sln[15][nMarkersUp + nMarkersLw + iPoint] = sections[iSec][2]->x[iPoint];
+ Sln[16][nMarkersUp + nMarkersLw + iPoint] = sections[iSec][2]->xoc[iPoint];
+ Sln[17][nMarkersUp + nMarkersLw + iPoint] = sections[iSec][2]->y[iPoint];
+ Sln[18][nMarkersUp + nMarkersLw + iPoint] = sections[iSec][2]->z[iPoint];
+ Sln[19][nMarkersUp + nMarkersLw + iPoint] = sections[iSec][2]->vt[iPoint];
}
if (verbose > 2)
diff --git a/blast/src/DDriver.h b/blast/src/DDriver.h
index 2937d3d..ed598aa 100644
--- a/blast/src/DDriver.h
+++ b/blast/src/DDriver.h
@@ -28,7 +28,6 @@ private:
double Re; ///< Freestream Reynolds number.
double CFL0; ///< Initial CFL number for pseudo time integration.
double span; ///< Wing Span (Used only for drag computation, not used if 2D case).
- std::vector<bool> stagPtMvmt; ///< Flag to account for stagnation point movements.
Solver *tSolver; ///< Pseudo-time solver.
int solverExitCode; ///< Exit code of viscous calculations.
std::vector<std::vector<std::vector<size_t>>> convergenceStatus; ///< Vector containing points that did not converged.
@@ -40,34 +39,19 @@ public:
Driver(double _Re, double _Minf, double _CFL0, size_t nSections,
double xtrF_top=-1., double xtrF_bot=-1., double _span = 0., unsigned int _verbose = 1);
~Driver();
- int run(unsigned int const couplIter);
+ int run();
+ void reset();
// 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 getAverageTransition(size_t const iRegion) const;
double getRe() const { return Re; };
std::vector<std::vector<double>> getSolution(size_t iSec);
- // Setters.
- void setGlobMesh(size_t iSec, size_t iRegion, std::vector<double> _x, std::vector<double> _y, std::vector<double> _z);
- void setVelocities(size_t iSec, size_t iRegion, std::vector<double> _vx, std::vector<double> _vy, std::vector<double> _vz);
- void setMe(size_t iSec, size_t iRegion, std::vector<double> _Me);
- void setRhoe(size_t iSec, size_t iRegion, std::vector<double> _rhoe);
- void setxxExt(size_t iSec, size_t iRegion, std::vector<double> _xxExt);
- void setDeltaStarExt(size_t iSec, size_t iRegion, std::vector<double> _deltaStarExt);
- void setUeOld(size_t iSec, size_t iRegion, std::vector<double> _ue);
-
- // Extractors.
- std::vector<double> extractBlowingVel(size_t iSec, size_t iRegion) const;
- std::vector<double> extractxx(size_t iSec, size_t iRegion) const;
- std::vector<double> extractDeltaStar(size_t iSec, size_t iRegion) const;
- std::vector<double> extractxCoord(size_t iSec, size_t iRegion) const;
- std::vector<double> extractUe(size_t iSec, size_t iRegion) const;
+ void addSection(size_t iSec, BoundaryLayer * &bl){sections[iSec].push_back(bl);};
private:
void checkWaves(size_t iPoint, BoundaryLayer *bl);
- void averageTransition(size_t iPoint, BoundaryLayer *bl, int forced);
+ void averageTransition(size_t iPoint, BoundaryLayer *bl, bool forced);
void computeSectionalDrag(std::vector<BoundaryLayer *> const &bl, size_t i);
void computeTotalDrag();
void computeBlwVel();
diff --git a/blast/src/DEdge.cpp b/blast/src/DEdge.cpp
deleted file mode 100644
index 762d63c..0000000
--- a/blast/src/DEdge.cpp
+++ /dev/null
@@ -1,59 +0,0 @@
-#include "DEdge.h"
-
-using namespace blast;
-
-/**
- * @brief Construct a new Edge::Edge object.
- *
- */
-Edge::Edge()
-{
- Me.resize(0, 0.);
- vt.resize(0, 0.);
- rhoe.resize(0, 0.);
- deltaStar.resize(0, 0.);
- vtOld.resize(0, 0.);
-}
-
-/**
- * @brief Return the maximum inviscid Mach number in the region.
- *
- * @return double
- */
-double Edge::getMaxMach() const
-{
- if (Me.size() <= 0)
- throw std::runtime_error("blast::Edge Invalid Mach number vector.");
- double Mmax = 0.0;
- for (size_t i = 0; i < Me.size(); ++i)
- if (Me[i] > Mmax)
- Mmax = Me[i];
- return Mmax;
-}
-
-// Setters.
-void Edge::setMe(std::vector<double> const _Me)
-{
- Me.resize(_Me.size(), 0.);
- Me = _Me;
-}
-void Edge::setVt(std::vector<double> const _vt)
-{
- vt.resize(_vt.size(), 0.);
- vt = _vt;
-}
-void Edge::setRhoe(std::vector<double> const _rhoe)
-{
- rhoe.resize(_rhoe.size(), 0.);
- rhoe = _rhoe;
-}
-void Edge::setDeltaStar(std::vector<double> const _deltaStar)
-{
- deltaStar.resize(_deltaStar.size(), 0.);
- deltaStar = _deltaStar;
-}
-void Edge::setVtOld(std::vector<double> const _vtOld)
-{
- vtOld.resize(_vtOld.size(), 0.);
- vtOld = _vtOld;
-}
diff --git a/blast/src/DEdge.h b/blast/src/DEdge.h
deleted file mode 100644
index 0833564..0000000
--- a/blast/src/DEdge.h
+++ /dev/null
@@ -1,44 +0,0 @@
-#ifndef DEDGE_H
-#define DEDGE_H
-
-#include "blast.h"
-#include <vector>
-#include <iostream>
-
-namespace blast
-{
-
-/**
- * @brief Inviscid quantities at the edge of the boundary layer of a BoundaryLayer region.
- * @author Paul Dechamps
- */
-class BLAST_API Edge
-{
-private:
- std::vector<double> Me; ///< Mach number.
- std::vector<double> vt; ///< Velocity.
- std::vector<double> rhoe; ///< Density.
- std::vector<double> deltaStar; ///< Displacement thickness.
- std::vector<double> vtOld; ///< Previous velocity.
-
-public:
- Edge();
- ~Edge(){};
-
- // Getters.
- double getMe(size_t iPoint) const { return Me[iPoint]; };
- double getVt(size_t iPoint) const { return vt[iPoint]; };
- double getRhoe(size_t iPoint) const { return rhoe[iPoint]; };
- double getDeltaStar(size_t iPoint) const { return deltaStar[iPoint]; };
- double getVtOld(size_t iPoint) const { return vtOld[iPoint]; };
- double getMaxMach() const;
-
- // Setters.
- void setMe(std::vector<double> const _Me);
- void setVt(std::vector<double> const _vt);
- void setRhoe(std::vector<double> const _rhoe);
- void setDeltaStar(std::vector<double> const _deltaStar);
- void setVtOld(std::vector<double> const _vtOld);
-};
-} // namespace blast
-#endif // DEDGE_H
diff --git a/blast/src/DFluxes.cpp b/blast/src/DFluxes.cpp
index c3beca0..3cf81cc 100644
--- a/blast/src/DFluxes.cpp
+++ b/blast/src/DFluxes.cpp
@@ -79,7 +79,7 @@ VectorXd Fluxes::blLaws(size_t iPoint, BoundaryLayer *bl, std::vector<double> u)
{
// Laminar closure relations.
std::vector<double> lamParam(8, 0.);
- bl->closSolver->laminarClosures(lamParam, u[0], u[1], bl->Ret[iPoint], bl->blEdge->getMe(iPoint), bl->name);
+ bl->closSolver->laminarClosures(lamParam, u[0], u[1], bl->Ret[iPoint], bl->Me[iPoint], bl->name);
bl->Hstar[iPoint] = lamParam[0];
bl->Hstar2[iPoint] = lamParam[1];
bl->Hk[iPoint] = lamParam[2];
@@ -94,7 +94,7 @@ VectorXd Fluxes::blLaws(size_t iPoint, BoundaryLayer *bl, std::vector<double> u)
{
// Turbulent closure relations.
std::vector<double> turbParam(8, 0.);
- bl->closSolver->turbulentClosures(turbParam, u[0], u[1], u[4], bl->Ret[iPoint], bl->blEdge->getMe(iPoint), bl->name);
+ bl->closSolver->turbulentClosures(turbParam, u[0], u[1], u[4], bl->Ret[iPoint], bl->Me[iPoint], bl->name);
bl->Hstar[iPoint] = turbParam[0];
bl->Hstar2[iPoint] = turbParam[1];
bl->Hk[iPoint] = turbParam[2];
@@ -112,16 +112,16 @@ VectorXd Fluxes::blLaws(size_t iPoint, BoundaryLayer *bl, std::vector<double> u)
}
// Gradients computation.
- double dx = (bl->mesh->getLoc(iPoint) - bl->mesh->getLoc(iPoint - 1));
- double dxExt = (bl->mesh->getExt(iPoint) - bl->mesh->getExt(iPoint - 1));
+ double dx = (bl->loc[iPoint] - bl->loc[iPoint-1]);
+ double dxExt = (bl->xxExt[iPoint] - bl->xxExt[iPoint-1]);
double dt_dx = (u[0] - bl->u[(iPoint - 1) * nVar + 0]) / dx;
double dH_dx = (u[1] - bl->u[(iPoint - 1) * nVar + 1]) / dx;
double due_dx = (u[3] - bl->u[(iPoint - 1) * nVar + 3]) / dx;
double dct_dx = (u[4] - bl->u[(iPoint - 1) * nVar + 4]) / dx;
double dHstar_dx = (bl->Hstar[iPoint] - bl->Hstar[iPoint - 1]) / dx;
- double dueExt_dx = (bl->blEdge->getVt(iPoint) - bl->blEdge->getVt(iPoint - 1)) / dxExt;
- double ddeltaStarExt_dx = (bl->blEdge->getDeltaStar(iPoint) - bl->blEdge->getDeltaStar(iPoint - 1)) / dxExt;
+ double dueExt_dx = (bl->vt[iPoint] - bl->vt[iPoint-1]) / dxExt;
+ double ddeltaStarExt_dx = (bl->deltaStarExt[iPoint] - bl->deltaStarExt[iPoint-1]) / dxExt;
double c = 4 / (M_PI * dx);
double cExt = 4 / (M_PI * dxExt);
@@ -136,7 +136,7 @@ VectorXd Fluxes::blLaws(size_t iPoint, BoundaryLayer *bl, std::vector<double> u)
dN_dx = (u[2] - bl->u[(iPoint - 1) * nVar + 2]) / dx;
}
- double Mea = bl->blEdge->getMe(iPoint);
+ double Mea = bl->Me[iPoint];
double Hstara = bl->Hstar[iPoint];
double Hstar2a = bl->Hstar2[iPoint];
double Hka = bl->Hk[iPoint];
diff --git a/blast/src/DSolver.cpp b/blast/src/DSolver.cpp
index fc3b94d..541a799 100644
--- a/blast/src/DSolver.cpp
+++ b/blast/src/DSolver.cpp
@@ -2,6 +2,7 @@
#include "DBoundaryLayer.h"
#include "DFluxes.h"
#include <Eigen/Dense>
+#include <iostream>
using namespace Eigen;
using namespace tbox;
@@ -24,7 +25,6 @@ Solver::Solver(double _CFL0, double _Minf, double Re, unsigned int _maxIt, doubl
tol = _tol;
itFrozenJac0 = _itFrozenJac;
Minf = std::max(_Minf, 0.1);
- initSln = true;
sysMatrix = new Fluxes(Re);
}
@@ -47,25 +47,20 @@ Solver::~Solver()
void Solver::initialCondition(size_t iPoint, BoundaryLayer *bl)
{
size_t nVar = bl->getnVar();
- if (initSln)
- {
- bl->u[iPoint * nVar + 0] = bl->u[(iPoint - 1) * nVar + 0];
- bl->u[iPoint * nVar + 1] = bl->u[(iPoint - 1) * nVar + 1];
- bl->u[iPoint * nVar + 2] = bl->u[(iPoint - 1) * nVar + 2];
- bl->u[iPoint * nVar + 3] = bl->blEdge->getVt(iPoint);
- bl->u[iPoint * nVar + 4] = bl->u[(iPoint - 1) * nVar + 4];
-
- bl->Ret[iPoint] = bl->Ret[iPoint - 1];
- bl->cd[iPoint] = bl->cd[iPoint - 1];
- bl->cf[iPoint] = bl->cf[iPoint - 1];
- bl->Hstar[iPoint] = bl->Hstar[iPoint - 1];
- bl->Hstar2[iPoint] = bl->Hstar2[iPoint - 1];
- bl->Hk[iPoint] = bl->Hk[iPoint - 1];
- bl->ctEq[iPoint] = bl->ctEq[iPoint - 1];
- bl->us[iPoint] = bl->us[iPoint - 1];
- bl->delta[iPoint] = bl->delta[iPoint - 1];
- bl->deltaStar[iPoint] = bl->deltaStar[iPoint - 1];
- }
+ for (auto k = 0; k < nVar; ++k)
+ bl->u[iPoint * nVar + k] = bl->u[(iPoint - 1) * nVar + k];
+ bl->u[iPoint * nVar + 3] = bl->vt[iPoint];
+
+ bl->Ret[iPoint] = bl->Ret[iPoint - 1];
+ bl->cd[iPoint] = bl->cd[iPoint - 1];
+ bl->cf[iPoint] = bl->cf[iPoint - 1];
+ bl->Hstar[iPoint] = bl->Hstar[iPoint - 1];
+ bl->Hstar2[iPoint] = bl->Hstar2[iPoint - 1];
+ bl->Hk[iPoint] = bl->Hk[iPoint - 1];
+ bl->ctEq[iPoint] = bl->ctEq[iPoint - 1];
+ bl->us[iPoint] = bl->us[iPoint - 1];
+ bl->delta[iPoint] = bl->delta[iPoint - 1];
+ bl->deltaStar[iPoint] = bl->deltaStar[iPoint - 1];
if (bl->regime[iPoint] == 1 && bl->u[iPoint * nVar + 4] <= 0)
bl->u[iPoint * nVar + 4] = 0.03;
}
@@ -87,12 +82,12 @@ int Solver::integration(size_t iPoint, BoundaryLayer *bl)
sln0[i] = bl->u[iPoint * nVar + i];
// Initialize time integration variables.
- double dx = bl->mesh->getLoc(iPoint) - bl->mesh->getLoc(iPoint - 1);
+ double dx = bl->loc[iPoint] - bl->loc[iPoint-1];
// Initial time step.
double CFL = CFL0;
unsigned int itFrozenJac = itFrozenJac0;
- double timeStep = setTimeStep(CFL, dx, bl->blEdge->getVt(iPoint));
+ double timeStep = setTimeStep(CFL, dx, bl->vt[iPoint]);
MatrixXd IMatrix(5, 5);
IMatrix = MatrixXd::Identity(5, 5) / timeStep;
@@ -124,23 +119,12 @@ int Solver::integration(size_t iPoint, BoundaryLayer *bl)
sysRes = sysMatrix->computeResiduals(iPoint, bl);
normSysRes = (sysRes + slnIncr / timeStep).norm();
- if (std::isnan(normSysRes))
- {
- resetSolution(iPoint, bl, sln0, true);
- sysRes = sysMatrix->computeResiduals(iPoint, bl);
- CFL = 1.0;
- timeStep = setTimeStep(CFL, dx, bl->blEdge->getVt(iPoint));
- IMatrix = MatrixXd::Identity(5, 5) / timeStep;
- normSysRes = (sysRes + slnIncr / timeStep).norm();
- itFrozenJac = 1;
- }
-
if (normSysRes <= tol * normSysRes0)
return 0; // Successfull exit.
// CFL adaptation.
CFL = std::max(CFL0 * pow(normSysRes0 / normSysRes, 0.7), 0.1);
- timeStep = setTimeStep(CFL, dx, bl->blEdge->getVt(iPoint));
+ timeStep = setTimeStep(CFL, dx, bl->vt[iPoint]);
IMatrix = MatrixXd::Identity(5, 5) / timeStep;
innerIters++;
@@ -213,7 +197,7 @@ void Solver::resetSolution(size_t iPoint, BoundaryLayer *bl, std::vector<double>
if (bl->regime[iPoint] == 0)
{
std::vector<double> lamParam(8, 0.);
- bl->closSolver->laminarClosures(lamParam, bl->u[iPoint * nVar + 0], bl->u[iPoint * nVar + 1], bl->Ret[iPoint], bl->blEdge->getMe(iPoint), bl->name);
+ bl->closSolver->laminarClosures(lamParam, bl->u[iPoint * nVar + 0], bl->u[iPoint * nVar + 1], bl->Ret[iPoint], bl->Me[iPoint], bl->name);
bl->Hstar[iPoint] = lamParam[0];
bl->Hstar2[iPoint] = lamParam[1];
bl->Hk[iPoint] = lamParam[2];
@@ -226,7 +210,7 @@ void Solver::resetSolution(size_t iPoint, BoundaryLayer *bl, std::vector<double>
else if (bl->regime[iPoint] == 1)
{
std::vector<double> turbParam(8, 0.);
- bl->closSolver->turbulentClosures(turbParam, bl->u[iPoint * nVar + 0], bl->u[iPoint * nVar + 1], bl->u[iPoint * nVar + 4], bl->Ret[iPoint], bl->blEdge->getMe(iPoint), bl->name);
+ bl->closSolver->turbulentClosures(turbParam, bl->u[iPoint * nVar + 0], bl->u[iPoint * nVar + 1], bl->u[iPoint * nVar + 4], bl->Ret[iPoint], bl->Me[iPoint], bl->name);
bl->Hstar[iPoint] = turbParam[0];
bl->Hstar2[iPoint] = turbParam[1];
bl->Hk[iPoint] = turbParam[2];
diff --git a/blast/src/DSolver.h b/blast/src/DSolver.h
index f6ce735..16c73ba 100644
--- a/blast/src/DSolver.h
+++ b/blast/src/DSolver.h
@@ -19,7 +19,6 @@ private:
unsigned int maxIt; ///< Maximum number of iterations.
double tol; ///< Convergence tolerance.
unsigned int itFrozenJac0; ///< Number of iterations between which the Jacobian is frozen.
- bool initSln; ///< Flag to initialize the solution at the point.
double const gamma = 1.4; ///< Air heat capacity ratio.
double Minf; ///< Freestream Mach number.
@@ -40,7 +39,6 @@ public:
void setCFL0(double _CFL0);
void setmaxIt(unsigned int _maxIt) { maxIt = _maxIt; };
void setitFrozenJac(unsigned int _itFrozenJac) { itFrozenJac0 = _itFrozenJac; };
- void setinitSln(bool _initSln) { initSln = _initSln; };
private:
double setTimeStep(double CFL, double dx, double invVel) const;
diff --git a/blast/src/blast.h b/blast/src/blast.h
index 71cdab8..382e040 100644
--- a/blast/src/blast.h
+++ b/blast/src/blast.h
@@ -43,8 +43,6 @@ class Solver;
// Data structures
class BoundaryLayer;
-class Discretization;
-class Edge;
// Other
class Closures;
diff --git a/blast/tests/dart/naca0012_2D.py b/blast/tests/dart/naca0012_2D.py
index 4bdfc11..361b760 100644
--- a/blast/tests/dart/naca0012_2D.py
+++ b/blast/tests/dart/naca0012_2D.py
@@ -55,7 +55,7 @@ def cfgInviscid(nthrds, verb):
'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' : 5, 'yLgt' : 5, 'msF' : 1, 'msTe' : 0.01, 'msLe' : 0.001}, # parameters for input file model
+ 'Pars' : {'xLgt' : 50, 'yLgt' : 50, 'msF' : 8.33, 'msTe' : 0.01, 'msLe' : 0.001}, # parameters for input file model
'Dim' : 2, # problem dimension
'Format' : 'gmsh', # save format (vtk or gmsh)
# Markers
@@ -92,14 +92,13 @@ def cfgBlast(verb):
'Minf' : 0.2, # 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': 25, # Maximum number of iterations
+ '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' : [None, 0.4],# Forced transition location
- 'interpolator' : 'Matching',
- 'nDim' : 2
+ 'interpolator' : 'Matching'
}
def main():
@@ -112,7 +111,7 @@ def main():
vcfg = cfgBlast(args.verb)
tms['pre'].start()
- coupler, iSolverAPI, vSolver = viscUtils.initBlast(icfg, vcfg)
+ coupler, isol, vsol = viscUtils.initBlast(icfg, vcfg)
tms['pre'].stop()
print(ccolors.ANSI_BLUE + 'PySolving...' + ccolors.ANSI_RESET)
@@ -123,11 +122,10 @@ def main():
# 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()))
+ 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(vSolver)
- vSolution['Cdt_int'] = vSolution['Cdf'] + iSolverAPI.getCd()
+ vSolution = viscUtils.getSolution(isol.sec, write=True, toW='all')
tms['total'].stop()
print(ccolors.ANSI_BLUE + 'PyTiming...' + ccolors.ANSI_RESET)
@@ -139,13 +137,13 @@ def main():
# Test solution
print(ccolors.ANSI_BLUE + 'PyTesting...' + ccolors.ANSI_RESET)
tests = CTests()
- tests.add(CTest('Cl', iSolverAPI.getCl(), 0.230, 5e-2)) # XFOIL 0.2325
- tests.add(CTest('Cd wake', vSolution['Cdt_w'], 0.0056, 1e-3, forceabs=True)) # XFOIL 0.00531
- tests.add(CTest('Cd integral', vSolution['Cdt_int'], 0.0059, 1e-3, forceabs=True)) # XFOIL 0.00531
- tests.add(CTest('Cdf', vSolution['Cdf'], 0.00465, 1e-3, forceabs=True)) # XFOIL 0.00084, Cdf = 0.00447
- tests.add(CTest('xtr_top', vSolution['xtrT'], 0.20, 0.03, forceabs=True)) # XFOIL 0.1877
- tests.add(CTest('xtr_bot', vSolution['xtrB'], 0.40, 0.01, forceabs=True)) # XFOIL 0.4998
- tests.add(CTest('Iterations', len(aeroCoeffs), 17, 0, forceabs=True))
+ 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']), 19, 0, forceabs=True))
tests.run()
# Show results
@@ -158,7 +156,7 @@ def main():
np.column_stack((xfoilCp[:,0], xfoilCp[:,1])),
np.column_stack((iCp[:,0], iCp[:,3])),
np.column_stack((xfoilCpInv[:,0], xfoilCpInv[:,1]))],
- 'labels': ['Blast (VII)', 'XFOIL VII', 'DART (inviscid)', 'XFOIL (inviscid)'],
+ 'labels': ['Blaster (VII)', 'XFOIL (VII)', 'DART (inviscid)', 'XFOIL (inviscid)'],
'lw': [3, 3, 2, 2],
'color': ['firebrick', 'darkblue', 'firebrick', 'darkblue'],
'ls': ['-', '-', '--', '--'],
@@ -172,9 +170,9 @@ def main():
viscUtils.plot(plotcp)
xfoilBl = viscUtils.read('../../blast/models/references/blXfoil_n0012.dat', delim=None, skip = 1)
- plotcf = {'curves': [np.column_stack((vSolution['x'], vSolution['Cf'])),
+ plotcf = {'curves': [np.column_stack((vSolution[0]['x'], vSolution[0]['cf']*vSolution[0]['ue']*vSolution[0]['ue'])),
np.column_stack((xfoilBl[:,1], xfoilBl[:,6]))],
- 'labels': ['Blast', 'XFOIL'],
+ 'labels': ['Blaster', 'XFOIL'],
'lw': [3, 3],
'color': ['firebrick', 'darkblue'],
'ls': ['-', '-'],
diff --git a/blast/tests/dart/rae2822_3D.py b/blast/tests/dart/rae2822_3D.py
index 92add20..3e3a32f 100644
--- a/blast/tests/dart/rae2822_3D.py
+++ b/blast/tests/dart/rae2822_3D.py
@@ -93,7 +93,7 @@ def cfgBlast(verb):
'Minf' : 0.8, # Freestream Mach number (used for the computation of the time step only)
'CFL0' : 1, # Inital CFL number of the calculation
- 'sections' : np.linspace(0.01, 0.95, 30),
+ 'sections' : np.linspace(0.01, 0.95, 3),
'writeSections': [0.2, 0.4, 0.6, 0.8, 1.0],
'Sym':[0.],
'span': 1.,
@@ -105,7 +105,7 @@ def cfgBlast(verb):
'saveTag': 4,
'Verb': verb, # Verbosity level of the solver
- 'couplIter': 50, # Maximum number of iterations
+ 'couplIter': 5, # Maximum number of iterations
'couplTol' : 5e-2, # 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.
@@ -131,7 +131,7 @@ def main():
vcfg['vMsh'] = vMsh
tms['pre'].start()
- coupler, iSolverAPI, vSolver = viscUtils.initBlast(icfg, vcfg)
+ coupler, isol, vsol = viscUtils.initBlast(icfg, vcfg)
tms['pre'].stop()
print(ccolors.ANSI_BLUE + 'PySolving...' + ccolors.ANSI_RESET)
@@ -142,19 +142,13 @@ def main():
# Display results.
print(ccolors.ANSI_BLUE + 'PyRes...' + ccolors.ANSI_RESET)
print(' Re M alpha Cl Cd Cd_wake 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} {8:8.4f}'.format(vcfg['Re']/1e6, iSolverAPI.getMinf(), iSolverAPI.getAoA()*180/math.pi, iSolverAPI.getCl(), vSolver.Cdf + iSolverAPI.getCd(), vSolver.Cdt, vSolver.Cdp, vSolver.Cdf, iSolverAPI.getCm()))
+ 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} {8:8.4f}'.format(vcfg['Re']/1e6, isol.getMinf(), isol.getAoA()*180/math.pi, isol.getCl(), vsol.Cdf + isol.getCd(), vsol.Cdt, vsol.Cdp, vsol.Cdf, isol.getCm()))
# Write results to file.
- vSolution = viscUtils.getSolution(vSolver)
+ vSolution = viscUtils.getSolution(isol.sec, write=True, toW='all')
- # Write results to file.
- for iSec in range(len(iSolverAPI.cfg['EffSections'])):
- vSolution = viscUtils.getSolution(vSolver, iSec)
- viscUtils.write(vSolution, vSolver.getRe(), sfx='slice'+str(iSec))
- vSolution['Cdt_int'] = vSolver.Cdf + iSolverAPI.getCd()
-
- # Save pressure coefficient
- iSolverAPI.save(sfx='_viscous')
+ # Save results
+ isol.save(sfx='_viscous')
tms['total'].stop()
print(ccolors.ANSI_BLUE + 'PyTiming...' + ccolors.ANSI_RESET)
@@ -166,9 +160,11 @@ def main():
# Test solution
print(ccolors.ANSI_BLUE + 'PyTesting...' + ccolors.ANSI_RESET)
tests = CTests()
- tests.add(CTest('Cl', iSolverAPI.getCl(), 0.128, 5e-2))
- tests.add(CTest('Cd wake', vSolution['Cdt_int'], 0.0132, 1e-3, forceabs=True))
- tests.add(CTest('Iterations', len(aeroCoeffs), 5, 0, forceabs=True))
+ tests.add(CTest('Cl', isol.getCl(), 0.135, 5e-2))
+ tests.add(CTest('Cd wake', vsol.Cdt, 0.0074, 1e-3, forceabs=True))
+ tests.add(CTest('Cd integral', vsol.Cdf + isol.getCd(), 0.0140, 1e-3, forceabs=True))
+ tests.add(CTest('Cdf', vsol.Cdf, 0.0061, 1e-3, forceabs=True))
+ tests.add(CTest('Iterations', len(aeroCoeffs), 3, 0, forceabs=True))
tests.run()
# eof
diff --git a/blast/utils.py b/blast/utils.py
index 91a0aa4..b32492a 100644
--- a/blast/utils.py
+++ b/blast/utils.py
@@ -101,77 +101,109 @@ def mesh(file, pars):
return msh
-def getSolution(vSolver, iSec=0):
- """Extract viscous solution.
- Args
- ----
- - vSolver: blast::Driver class. Driver of the viscous calculations
- - iSec (int): Marker of the section (default: 0)
- """
- if iSec<0:
- raise RuntimeError("blast::viscU Invalid section number", iSec)
+def getSolution(sections, write=False, toW=['x', 'xelm', 'theta', 'H', 'deltaStar', 'Cf', 'blowingVelocity'], sfx=''):
+ """ Extract viscous solution.
+
+ Parameters
+ ----------
+
+ - sections: list of blast::BoundaryLayer class. Sections of the boundary layer solver.
+ - write: bool. Flag to write the results in a file.
+ - sfx: str. Suffix to add to the file name.
- solverOutput = vSolver.getSolution(iSec)
-
- sln = { 'theta' : solverOutput[0],
- 'H' : solverOutput[1],
- 'N' : solverOutput[2],
- 'ue' : solverOutput[3],
- 'Ct' : solverOutput[4],
- 'deltaStar': solverOutput[5],
- 'Ret' : solverOutput[6],
- 'Cd' : solverOutput[7],
- 'Cf' : np.asarray(solverOutput[8])*np.asarray(solverOutput[3])**2,
- 'Hstar' : solverOutput[9],
- 'Hstar2' : solverOutput[10],
- 'Hk' : solverOutput[11],
- 'Cteq' : solverOutput[12],
- 'Us' : solverOutput[13],
- 'delta' : solverOutput[14],
- 'x' : solverOutput[15],
- 'xoc' : solverOutput[16],
- 'y' : solverOutput[17],
- 'z' : solverOutput[18],
- 'ueInv' : solverOutput[19],
- 'xtrT' : vSolver.getAvrgxoctr(0),
- 'xtrB' : vSolver.getAvrgxoctr(1),
- 'Cdt_w' : vSolver.Cdt_sec[iSec],
- 'Cdf' : vSolver.Cdf_sec[iSec],
- 'Cdp' : vSolver.Cdp_sec[iSec]
- }
- return sln
-
-def write(wData, Re, toW=['deltaStar', 'H', 'Hstar', 'Cf', 'Ct', 'ue', 'ueInv', 'delta'], sfx=''):
- """Write the results in bl files
+ Return
+ ------
+
+ - sol: list[dicts]. Dictionary containing the boundary layer solution.
"""
- import os
- if not os.path.exists('blSlices'):
- os.makedirs('blSlices')
- # Write
- print('Writing file: /blSlices/bl'+sfx+'.dat...', end = '')
- f = open('blSlices/bl'+sfx+'.dat', 'w+')
-
- f.write('$Sectional aerodynamic coefficients\n')
- f.write(' Re Cdw Cdp Cdf xtr_top xtr_bot\n')
- f.write('{0:15.6f} {1:15.6f} {2:15.6f} {3:15.6f} {4:15.6f} {5:15.6f}\n'.format(Re, wData['Cdt_w'], wData['Cdp'],
- wData['Cdf'], wData['xtrT'],
- wData['xtrB']))
- f.write('$Boundary layer variables\n')
- f.write('{0:>15s} {1:>15s} {2:>15s} {3:>15s}'.format('x','y', 'z', 'xoc'))
-
-
- for s in toW:
- f.write(' {0:>15s}'.format(s))
- f.write('\n')
-
- for i in range(len(wData['x'])):
- f.write('{0:>15.6f} {1:>15.6f} {2:>15.6f} {3:>15.6f}'.format(wData['x'][i], wData['y'][i], wData['z'][i], (wData['x'][i] - min(wData['x'])) / (max(wData['x']) - min(wData['x']))))
- for s in toW:
- f.write(' {0:15.6f}'.format(wData[s][i]))
- f.write('\n')
-
- f.close()
- print('done.')
+
+ nVar = sections[0][0].getnVar()
+ sol = []
+
+ varKeys = ['theta', 'H', 'N', 'ue', 'Ct']
+ attrKeys = ['x', 'y', 'z', 'xoc', 'deltaStar', \
+ 'cd', 'cf', 'Hstar', 'Hstar2', 'Hk', 'ctEq', 'us', 'delta',\
+ 'vt', 'Ret']
+ elemsKeys = ['xelm', 'yelm', 'zelm']
+
+ for iSec, sec in enumerate(sections):
+ nNodes = 0
+ for side in sec:
+ nNodes += len(side.x)
+
+ sol.append({})
+
+ for key in elemsKeys:
+ sol[iSec][key] = np.zeros(nNodes-len(sec))
+ for key in varKeys + attrKeys:
+ sol[iSec][key] = np.zeros(nNodes)
+
+ for side in sec:
+ nNodes_side = side.nNodes
+ if side.name == "upper":
+ for key in attrKeys:
+ sol[iSec][key][:nNodes_side] = np.flip(getattr(side, key))
+ for k, key in enumerate(varKeys):
+ sol[iSec][key][:nNodes_side] = np.flip(side.u[k::nVar])
+ else:
+ if side.name == "lower":
+ slicer = np.s_[sec[0].x.size():sec[0].x.size()+sec[1].x.size()]
+ elif side.name == "wake":
+ slicer = np.s_[sec[0].x.size()+sec[1].x.size():]
+
+ for key in attrKeys:
+ sol[iSec][key][slicer] = getattr(side, key)
+ for k, key in enumerate(varKeys):
+ sol[iSec][key][slicer] = side.u[k::nVar]
+
+ # Compute elements cgs
+ elemsCoord = np.zeros((nNodes_side-1, 3))
+ for k, key in enumerate(elemsKeys):
+ elemsCoord[:,k] = 0.5 * (np.array(getattr(side, key[0]))[1:] + np.array(getattr(side, key[0]))[:-1])
+ if side.name == "upper":
+ for k, key in enumerate(elemsKeys):
+ sol[iSec][key][:nNodes_side-1] = np.flip(elemsCoord[:,k])
+ else:
+ if side.name == "lower":
+ slicer = np.s_[(sec[0].x.size()-1):(sec[0].x.size()-1)+(sec[1].x.size()-1)]
+ elif side.name == "wake":
+ slicer = np.s_[(sec[0].x.size()-1)+(sec[1].x.size()-1):]
+ for k, key in enumerate(elemsKeys):
+ sol[iSec][key][slicer] = elemsCoord[:,k]
+ sol[iSec]['blowingVelocity'] = np.concatenate((np.flip(sec[0].blowingVelocity),\
+ sec[1].blowingVelocity, sec[2].blowingVelocity))
+
+ if write:
+ import datetime
+ if toW == 'all':
+ toW = varKeys + attrKeys
+ import os
+ if not os.path.exists('blSlices'):
+ os.makedirs('blSlices')
+ print('Writing file: /blSlices/bl'+sfx+'.dat...', end = '')
+ for iSec, sec in enumerate(sections):
+ f = open('blSlices/bl'+str(iSec)+sfx+'.dat', 'w+')
+ f.write('# BLASTER boundary layer output file\n')
+ f.write('# Generated on '+datetime.datetime.now().strftime('%Y-%m-%d %H:%M')+'\n')
+ f.write('\n')
+ f.write('# Section '+str(iSec)+'\n')
+ f.write('# ')
+ for side in sec:
+ f.write('xtr_'+side.name+' = '+str(round(side.xtr,6))+' ')
+ f.write('\n')
+ f.write('# Stagnation point at index '+str(sec[0].x.size()-1)+'\n')
+ f.write('# Wake at index '+str(sec[0].x.size()+sec[1].x.size())+'\n')
+
+ for key in toW:
+ f.write(' {0:>21s}'.format(key))
+ f.write('\n')
+ for i in range(len(sol[iSec]['x'])):
+ for key in toW:
+ f.write(' {0:21.16f}'.format(sol[iSec][key][i]))
+ f.write('\n')
+ f.close()
+ print('done.')
+ return sol
def read(filename, delim=',', skip=1):
"""Read from file and store in data array
@@ -193,7 +225,7 @@ def plot(cfg):
if 'xreverse' in cfg and cfg['xreverse'] is True: plt.gca().invert_xaxis()
if 'title' in cfg: plt.title(cfg['title'])
if 'xlim' in cfg: plt.xlim(cfg['xlim'])
- if 'ylim' in cfg: plt.xlim(cfg['ylim'])
+ if 'ylim' in cfg: plt.ylim(cfg['ylim'])
if 'legend' in cfg and cfg['legend'] is True: plt.legend(frameon=False)
if 'xlabel' in cfg: plt.xlabel(cfg['xlabel'])
if 'ylabel' in cfg: plt.ylabel(cfg['ylabel'])
diff --git a/blast/validation/lannValidation.py b/blast/validation/lannValidation.py
index f5e8659..9833f14 100644
--- a/blast/validation/lannValidation.py
+++ b/blast/validation/lannValidation.py
@@ -116,7 +116,7 @@ def main():
vcfg['vMsh'] = vMsh
tms['pre'].start()
- coupler, iSolverAPI, vSolver = viscUtils.initBlast(icfg, vcfg)
+ coupler, isol, vsol = viscUtils.initBlast(icfg, vcfg)
tms['pre'].stop()
print(ccolors.ANSI_BLUE + 'PySolving...' + ccolors.ANSI_RESET)
@@ -127,19 +127,13 @@ def main():
# 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()))
+ 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(vSolver)
-
- # Write results to file.
- for iSec in range(len(iSolverAPI.cfg['EffSections'])):
- vSolution = viscUtils.getSolution(vSolver, iSec)
- viscUtils.write(vSolution, vSolver.getRe(), sfx='slice'+str(iSec))
- vSolution['Cdt_int'] = vSolver.Cdf + iSolverAPI.getCd()
+ vSolution = viscUtils.getSolution(isol.sec, write=True, toW='all')
# Save pressure coefficient
- iSolverAPI.save(sfx='_viscous')
+ isol.save(sfx='_viscous')
tms['total'].stop()
print(ccolors.ANSI_BLUE + 'PyTiming...' + ccolors.ANSI_RESET)
@@ -151,9 +145,10 @@ def main():
# Test solution
print(ccolors.ANSI_BLUE + 'PyTesting...' + ccolors.ANSI_RESET)
tests = CTests()
- tests.add(CTest('Cl', iSolverAPI.getCl(), 0.6007, 5e-2))
- tests.add(CTest('Cd wake', vSolution['Cdt_int'], 0.02756, 1e-3, forceabs=True))
- tests.add(CTest('Iterations', len(aeroCoeffs), 19, 0, forceabs=True))
+ tests.add(CTest('Cl', isol.getCl(), 0.6007, 5e-2))
+ tests.add(CTest('Cd wake', vsol.Cdt, 0.00348, 1e-3, forceabs=True))
+ tests.add(CTest('Cd int', isol.getCd() + vsol.Cdf, 0.0279, 1e-3, forceabs=True))
+ tests.add(CTest('Iterations', len(aeroCoeffs['Cl']), 4, 0, forceabs=True))
tests.run()
# eof
diff --git a/blast/validation/oneraValidation.py b/blast/validation/oneraValidation.py
index cb9af04..8a41eab 100644
--- a/blast/validation/oneraValidation.py
+++ b/blast/validation/oneraValidation.py
@@ -114,7 +114,7 @@ def main():
vcfg['vMsh'] = vMsh
tms['pre'].start()
- coupler, iSolverAPI, vSolver = viscUtils.initBlast(icfg, vcfg)
+ coupler, isol, vsol = viscUtils.initBlast(icfg, vcfg)
tms['pre'].stop()
print(ccolors.ANSI_BLUE + 'PySolving...' + ccolors.ANSI_RESET)
@@ -125,19 +125,13 @@ def main():
# 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()))
+ 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(vSolver)
-
- # Write results to file.
- for iSec in range(len(iSolverAPI.cfg['EffSections'])):
- vSolution = viscUtils.getSolution(vSolver, iSec)
- viscUtils.write(vSolution, vSolver.getRe(), sfx='slice'+str(iSec))
- vSolution['Cdt_int'] = vSolver.Cdf + iSolverAPI.getCd()
+ vSolution = viscUtils.getSolution(isol.sec, write=True, toW='all', sfx='onera')
# Save pressure coefficient
- iSolverAPI.save(sfx='_viscous')
+ isol.save(sfx='_viscous')
tms['total'].stop()
print(ccolors.ANSI_BLUE + 'PyTiming...' + ccolors.ANSI_RESET)
@@ -149,9 +143,10 @@ def main():
# Test solution
print(ccolors.ANSI_BLUE + 'PyTesting...' + ccolors.ANSI_RESET)
tests = CTests()
- tests.add(CTest('Cl', iSolverAPI.getCl(), 0.279, 5e-2))
- tests.add(CTest('Cd wake', vSolution['Cdt_int'], 0.0168, 1e-3, forceabs=True))
- tests.add(CTest('Iterations', len(aeroCoeffs), 19, 0, forceabs=True))
+ tests.add(CTest('Cl', isol.getCl(), 0.279, 5e-2))
+ tests.add(CTest('Cd wake', vsol.Cdt, 0.00471, 1e-3, forceabs=True))
+ tests.add(CTest('Cd int', isol.getCd() + vsol.Cdf, 0.01519, 1e-3, forceabs=True))
+ tests.add(CTest('Iterations', len(aeroCoeffs['Cl']), 13, 0, forceabs=True))
tests.run()
# eof
diff --git a/blast/validation/raeValidation.py b/blast/validation/raeValidation.py
index 9a99c22..8265822 100644
--- a/blast/validation/raeValidation.py
+++ b/blast/validation/raeValidation.py
@@ -28,11 +28,11 @@ from fwk.wutils import parseargs
import fwk
from fwk.testing import *
from fwk.coloring import ccolors
+import os.path
import math
def cfgInviscid(nthrds, verb):
- import os.path
# Parameters
return {
# Options
@@ -42,13 +42,13 @@ 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' : 5, 'yLgt' : 5, 'msF' : 1, 'msTe' : 0.01, 'msLe' : 0.002}, # parameters for input file model
+ 'Pars' : {'xLgt' : 50, 'yLgt' : 50, 'msF': 10, 'msTe' : 0.01, 'msLe' : 0.001}, # 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
+ 'Wing' : 'airfoil', # LIST of names of physical groups containing the lifting surface boundaryk
'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
@@ -64,7 +64,7 @@ def cfgInviscid(nthrds, verb):
'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)
+ 'LSolver' : 'SparseLu', # 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
@@ -98,8 +98,15 @@ def main():
icfg = cfgInviscid(args.k, args.verb)
vcfg = cfgBlast(args.verb)
+ 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)
+ print('msF =', icfg['Pars']['msF'])
tms['pre'].start()
- coupler, iSolverAPI, vSolver = viscUtils.initBlast(icfg, vcfg)
+ coupler, isol, vsol = viscUtils.initBlast(icfg, vcfg)
tms['pre'].stop()
print(ccolors.ANSI_BLUE + 'PySolving...' + ccolors.ANSI_RESET)
@@ -110,11 +117,10 @@ def main():
# 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()))
+ 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(vSolver)
- vSolution['Cdt_int'] = vSolution['Cdf'] + iSolverAPI.getCd()
+ vSolution = viscUtils.getSolution(isol.sec)
tms['total'].stop()
print(ccolors.ANSI_BLUE + 'PyTiming...' + ccolors.ANSI_RESET)
@@ -126,23 +132,27 @@ def main():
# 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.add(CTest('Cl', isol.getCl(), 0.766, 5e-2))
+ tests.add(CTest('Cd wake', vsol.Cdt, 0.0093, 1e-3, forceabs=True))
+ tests.add(CTest('Cd integral', isol.getCd() + vsol.Cdf, 0.0138, 1e-3, forceabs=True))
+ tests.add(CTest('Cdf', vsol.Cdf, 0.0069, 1e-3, forceabs=True))
+ tests.add(CTest('Iterations', len(aeroCoeffs['Cl']), 42, 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
+
# 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)'],
+ np.column_stack((iCp[:,0], iCp[:,3])),
+ expResults],
+ 'labels': ['Blast (VII)', 'DART (inviscid)', 'Experimental'],
'lw': [3, 3, 2, 2],
- 'color': ['darkblue', 'darkblue'],
- 'ls': ['-', '--'],
+ 'color': ['darkblue', 'darkblue', 'black'],
+ 'ls': ['-', '--', 'o'],
'reverse': True,
'xlim':[0, 1],
'yreverse': True,
@@ -152,13 +162,14 @@ def main():
}
viscUtils.plot(plotcp)
- plotcf = {'curves': [np.column_stack((vSolution['x'], vSolution['Cf']))],
+ 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$'
--
GitLab