import numpy as np
import blast
import fwk
import tbox
from blast.interfaces.DDataStructure import Group
class SolversInterface:
def __init__(self, _vSolver, _cfg):
self.vSolver = _vSolver
self.cfg = _cfg
if 'nSections' not in self.cfg:
self.cfg['nSections'] = len(self.cfg['sections'])
# Initialize data structures.
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
self.getWing()
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.copy()
self.vBnd[iSec][iReg].elemsCoord = self.iBnd[iReg].elemsCoord.copy()
self.vBnd[iSec][iReg].connectListNodes = self.iBnd[iReg].connectListNodes.copy()
self.vBnd[iSec][iReg].connectListElems = self.iBnd[iReg].connectListElems.copy()
elif self.cfg['interpolator'] == 'Rbf':
from blast.interfaces.interpolators.DRbfInterpolator import RbfInterpolator as interp
print('Initializing RBF interpolator.')
self.getVWing()
else:
raise RuntimeError('Incorrect interpolator specified in config.')
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.vtExt = [[np.zeros(0) for iReg in range(3)]\
for iSec in range(self.cfg['nSections'])]
self.vinit = False
def setViscousSolver(self):
"""Interpolate solution to viscous mesh and sets the inviscid boundary in the viscous solver.
"""
for iSec in range(self.cfg['nSections']):
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))
reg.setDeltaStarExt(self.deltaStarExt[iSec][i])
reg.setxxExt(self.xxExt[iSec][i])
reg.setVtExt(self.vtExt[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:
reg.blowingVel = np.asarray(self.sec[iSec][2].blowingVelocity)
elif 'vAirfoil' in reg.name:
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.vtExt[iSec][i] = np.asarray(self.sec[iSec][i].vt)
def interpolate(self, dir):
if dir == 'i2v':
self.interpolator.inviscidToViscous(self.iBnd, self.vBnd)
elif dir == 'v2i':
self.interpolator.viscousToInviscid(self.iBnd, self.vBnd)
else:
raise RuntimeError('Invalid direction', dir)
def updateMesh(self):
"""Collect the mesh on blowing boundaries and update the mesh in the viscous solver.
This function does not update connectivity lists
"""
if self.getnDim() != 2:
raise RuntimeError('SolverInterface::updateMesh - Only 2D meshes are supported.')
if self.cfg['interpolator'] != 'Matching':
raise RuntimeError('SolverInterface::updateMesh - Only Matching interpolator is supported.')
meshList = self.getMesh()
# Update the mesh in data structures
for ireg, reg in enumerate(self.iBnd):
for i, row in enumerate(reg.connectListNodes):
for idim in range(self.getnDim()):
reg.nodesCoord[i, idim] = meshList[ireg][row, idim]
# Recompute elements data
for ireg, reg in enumerate(self.iBnd):
for i in range(len(reg.nodesCoord[:,0])-1):
for idim in range(self.getnDim()):
reg.elemsCoord[i, idim] = (reg.nodesCoord[i, idim] + reg.nodesCoord[i+1, idim])/2
# Update the mesh in the data structures
for iSec in range(len(self.vBnd)):
for iReg, reg in enumerate(self.vBnd[iSec]):
self.vBnd[iSec][iReg].nodesCoord = self.iBnd[iReg].nodesCoord.copy()
self.vBnd[iSec][iReg].elemsCoord = self.iBnd[iReg].elemsCoord.copy()
# Set mesh in the viscous solver
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('Invalid region', reg.name)
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())
nodeRows = tbox.std_vector_int()
if self.cfg['interpolator'] == 'Matching':
for val in self.vBnd[iSec][iReg].getNodesRow(reg.name):
nodeRows.push_back(int(val))
else:
for i in range(len(self.vBnd[iSec][iReg].getNodesCoord(reg.name, xIdx))):
nodeRows.push_back(int(0))
connectElems = tbox.std_vector_int()
if self.cfg['interpolator'] == 'Matching':
for val in self.vBnd[iSec][iReg].getConnectList(reg.name):
connectElems.push_back(int(val))
else:
for i in range(len(self.vBnd[iSec][iReg].getNodesCoord(reg.name, xIdx))-1):
connectElems.push_back(int(0))
# x, y, z
xv = fwk.std_vector_double()
for val in self.vBnd[iSec][iReg].getNodesCoord(reg.name, xIdx):
xv.push_back(val)
yv = fwk.std_vector_double()
for val in self.vBnd[iSec][iReg].getNodesCoord(reg.name, yIdx):
yv.push_back(val)
zv = fwk.std_vector_double()
for val in self.vBnd[iSec][iReg].getNodesCoord(reg.name, zIdx):
zv.push_back(val)
# Set mesh
reg.setMesh(xv,
yv,
zv,
chord,
xMin,
nodeRows,
connectElems)
def updateStagnation(self):
"""Update the stagnation points in the viscous solver.
"""
for iSec in range(self.cfg['nSections']):
for reg in self.vBnd[iSec]:
reg.computeStagPoint()
def initializeViscousSolver(self):
"""Initialize the viscous sections in the viscous solver.
"""
if self.vinit == True:
return None
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.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('Invalid region', reg.name)
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())
nodeRows = tbox.std_vector_int()
if self.cfg['interpolator'] == 'Matching':
for val in self.vBnd[iSec][iReg].getNodesRow(reg.name):
nodeRows.push_back(int(val))
else:
for i in range(len(self.vBnd[iSec][iReg].getNodesCoord(reg.name, xIdx))):
nodeRows.push_back(int(0))
connectElems = tbox.std_vector_int()
if self.cfg['interpolator'] == 'Matching':
for val in self.vBnd[iSec][iReg].getConnectList(reg.name):
connectElems.push_back(int(val))
else:
for i in range(len(self.vBnd[iSec][iReg].getNodesCoord(reg.name, xIdx))-1):
connectElems.push_back(int(0))
# x, y, z
xv = fwk.std_vector_double()
for val in self.vBnd[iSec][iReg].getNodesCoord(reg.name, xIdx):
xv.push_back(val)
yv = fwk.std_vector_double()
for val in self.vBnd[iSec][iReg].getNodesCoord(reg.name, yIdx):
yv.push_back(val)
zv = fwk.std_vector_double()
for val in self.vBnd[iSec][iReg].getNodesCoord(reg.name, zIdx):
zv.push_back(val)
# Set mesh
reg.setMesh(xv,
yv,
zv,
chord,
xMin,
nodeRows,
connectElems)
# 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())
self.vtExt[iSec][i] = np.zeros(reg.getnNodes())
self.vinit = True
return 0
def getVWing(self):
"""Obtain the nodes' location and row and cg of all elements.
If the mesh is the viscous one, TE nodes are duplicated (by creating a new row).
They exist for the interpolation and for the solution computation.
If not, the nodes are pseudo-duplicated (without changing the row). They therefore
exist only for the interpolation but not for the solution.
"""
data = [np.zeros((0,4)) for _ in range(2)]
lwRows = []
for e in self.cfg['vMsh'].elems:
if e.ptag.name == 'wing' or e.ptag.name == 'wing_'\
or e.ptag.name == 'leadingEdge':
iRegion = 0
elif e.ptag.name == 'wake' or e.ptag.name == 'wake_':
iRegion = 1
else:
continue
# Get nodes
for nw in e.nodes:
if nw.row not in data[iRegion][:,3]:
data[iRegion] = np.vstack((data[iRegion], [nw.pos[0],\
nw.pos[1],\
nw.pos[2],\
nw.row]))
# If wing_ exists, the lower side is identifiable.
if e.ptag.name == 'wing_':
lwRows.append(nw.row)
if len(lwRows) == 0:
raise RuntimeError('Could not identify lower side.')
self.vlwIdx = []
for idx in range(len(data[0][:,3])):
if data[0][idx, 3] in lwRows:
self.vlwIdx.append(idx)
maxRowNb = max(data[0][:,3])
# Initialize viscous strutures
self.cfg['EffSections'] = np.empty(0)
for iReg, reg in enumerate(data):
for iy, y in enumerate(self.cfg['sections']):
# Find nearest value in the mesh
try:
idx = (np.abs(reg[:,1]-y).argmin())
except:
print('Region', reg[:,1], 'cannot be found through\
config input', y)
raise RuntimeError("Could not identify section.\
Possibly a missing tag")
section = reg[idx,1]
if iReg == 0:
self.cfg['EffSections'] = np.append(self.cfg['EffSections'], section)
nodesSec = reg[abs(reg[:,1]-section)<1e-3, :]
# Separate points on upper and lower side
if iReg == 0:
lower_points = []
upper_points = []
for point in nodesSec:
if point[3] in lwRows:
lower_points.append(point)
else:
upper_points.append(point)
# Sort points in selig format
lower_points = sorted(lower_points, key=lambda p: p[0])
upper_points = sorted(upper_points, key=lambda p: p[0], reverse=True)
nodesSec = np.row_stack((upper_points, lower_points))
nodesSec = np.row_stack((nodesSec, [nodesSec[0,0], nodesSec[0,1], nodesSec[0,2], nodesSec[0,3]+maxRowNb]))
else:
nodesSec = nodesSec[nodesSec[:,0].argsort()]
# Compute cg pos on the section (!= cg of 2D elems on the surface)
cgSec = np.zeros((len(nodesSec[:,0])-1, 4))
for i in range(len(cgSec)):
cgSec[i,:3] = (nodesSec[i,:3] + nodesSec[i+1,:3])/2
cgSec[i, 3] = i
self.vBnd[iy][iReg].initStructures(nodesSec.shape[0])
self.vBnd[iy][iReg].nodesCoord = nodesSec
self.vBnd[iy][iReg].elemsCoord = cgSec
# Inviscid solver methods
def run()->int:
"""Run inviscid solver
"""
raise NotImplementedError('SolverInterface - run not implemented')
def reset()->None:
"""Reset solution
"""
raise NotImplementedError('SolverInterface - reset not implemented')
def save(self)->None:
"""Save solution
"""
raise NotImplementedError('SolverInterface - save not implemented')
def writeCp(self)->None:
"""Write pressure coefficient
"""
raise NotImplementedError('SolverInterface - writeCp not implemented')
# Inviscid solver properties
def getVerbose(self)->int:
"""Get verbose
"""
raise NotImplementedError('SolverInterface - getVerbose not implemented')
# Flow properties and aerodynamic coefficients
def getAoA(self)->float:
"""Get angle of attack
"""
raise NotImplementedError('SolverInterface - getAoA not implemented')
def getCl(self)->float:
"""Get lift coefficient
"""
raise NotImplementedError('SolverInterface - getCl not implemented')
def getCd(self)->float:
"""Get inviscid drag coefficient
"""
raise NotImplementedError('SolverInterface - getCd not implemented')
def getCm(self)->float:
"""Get moment coefficient
"""
raise NotImplementedError('SolverInterface - getCm not implemented')
def getnDim(self)->int:
"""Get number of dimensions
"""
raise NotImplementedError('SolverInterface - getnDim not implemented')
def getWing(self)->None:
"""Get wing mesh
"""
raise NotImplementedError('SolverInterface - getWing not implemented')
def getMesh(self)->list:
"""Get mesh
"""
raise NotImplementedError('SolverInterface - getMesh not implemented')
def getInviscidBC(self)->None:
"""Get inviscid boundary condition
"""
raise NotImplementedError('SolverInterface - getInviscidBC not implemented')
def setBlowingVelocity(self)->None:
"""Set blowing velocity
"""
raise NotImplementedError('SolverInterface - setBlowingVelocity not implemented')
def getnNodesDomain(self)->int:
"""Get number of nodes in domain
"""
raise NotImplementedError('SolverInterface - getnNodesDomain not implemented')