Skip to content
Snippets Groups Projects
Commit a040af32 authored by Paul Dechamps's avatar Paul Dechamps :speech_balloon:
Browse files

(feat, fix) Added BC imposition and fixed duplicate LE points 

For whatever reason, ther LE element is duplicated in SU2 (although it is defined once in the mesh). So we need to be careful when imposing the blowing velocity
parent 1e178753
No related branches found
No related tags found
Pipeline #52565 failed
Hide whitespace changes
Inline Side-by-side
blast/interfaces/su2/SU2Interface.py
+ 99
− 80
View file @ a040af32
......@@ -24,9 +24,9 @@ import pysu2
class SU2Interface(SolversInterface):
def __init__(self, su2config, vconfig, task='analysis'):
filename = su2config.get('filename')
self.filename = su2config.get('filename')
self.have_mpi, self.comm, self.status, self.myid = self.__setup_mpi()
self.solver = pysu2.CSinglezoneDriver(filename, 1, self.comm)
self.solver = pysu2.CSinglezoneDriver(self.filename, 1, self.comm)
self.wingMarkersToID = {key: i for i, key in enumerate(self.solver.GetMarkerTags())}
self.wingTags = su2config.get('wingTags')
......@@ -61,7 +61,7 @@ class SU2Interface(SolversInterface):
self.comm.barrier()
return 1
def writeCp(self):
def writeCp(self, sfx=''):
""" Write Cp distribution around the airfoil on file.
"""
pass
......@@ -97,17 +97,14 @@ class SU2Interface(SolversInterface):
sound_speed_idx = self.solver.GetPrimitiveIndices()["SOUND_SPEED"]
nNodes_farfield = self.solver.GetNumberMarkerNodes(self.wingMarkersToID['farfield'])
vel_x = np.zeros(nNodes_farfield)
vel_y = np.zeros(nNodes_farfield)
sound_speed = np.zeros(nNodes_farfield)
mach_farfield = np.zeros(nNodes_farfield)
primitives = self.solver.Primitives()
for inode in range(nNodes_farfield):
nodeIndex = self.solver.GetMarkerNode(self.wingMarkersToID['farfield'], inode)
vel_x[inode] = self.solver.Primitives().Get(int(nodeIndex), vel_x_idx)
vel_y[inode] = self.solver.Primitives().Get(int(nodeIndex), vel_y_idx)
sound_speed[inode] = self.solver.Primitives().Get(int(nodeIndex), sound_speed_idx)
mach_farfield = np.sqrt(vel_x**2 + vel_y**2) / sound_speed
vel_x = primitives.Get(nodeIndex, vel_x_idx)
vel_y = primitives.Get(nodeIndex, vel_y_idx)
sound_speed = primitives.Get(nodeIndex, sound_speed_idx)
mach_farfield[inode] = np.sqrt(vel_x**2 + vel_y**2) / sound_speed
return np.mean(mach_farfield)
def getCl(self):
......@@ -120,12 +117,13 @@ class SU2Interface(SolversInterface):
return self.solver.Get_Mx()
def getVerbose(self):
return 1
return 2
def reset(self):
print('SU2Interface::Warning: Cannot reset solver.')
pass
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.
......@@ -133,69 +131,92 @@ class SU2Interface(SolversInterface):
------
- couplIter (int): Coupling iteration.
"""
print('primitive indices', self.solver.GetPrimitiveIndices())
print('self.solver.MarkerPrimitives()', self.solver.MarkerPrimitives(self.wingMarkersToID['airfoil']))
vel_x_idx = self.solver.GetPrimitiveIndices()["VELOCITY_X"]
vel_y_idx = self.solver.GetPrimitiveIndices()["VELOCITY_Y"]
sound_speed_idx = self.solver.GetPrimitiveIndices()["SOUND_SPEED"]
pressure_idx = self.solver.GetPrimitiveIndices()["PRESSURE"]
vel_x = np.zeros(self.iBnd[0][0].nodesCoord.shape[0])
vel_y = np.zeros(self.iBnd[0][0].nodesCoord.shape[0])
sound_speed = np.zeros(self.iBnd[0][0].nodesCoord.shape[0])
pressure = np.zeros(self.iBnd[0][0].nodesCoord.shape[0])
for inode, nodeIndex in enumerate(self.iBnd[0][0].nodesCoord[:,-1]):
vel_x[inode] = self.solver.Primitives().Get(int(nodeIndex), vel_x_idx)
vel_y[inode] = self.solver.Primitives().Get(int(nodeIndex), vel_y_idx)
sound_speed[inode] = self.solver.Primitives().Get(int(nodeIndex), sound_speed_idx)
pressure[inode] = self.solver.Primitives().Get(int(nodeIndex), pressure_idx)
print(f"Node {inode} - Velocity X: {vel_x[inode]}, Velocity Y: {vel_y[inode]}, Sound Speed: {sound_speed[inode]}, Pressure: {pressure[inode]}")
print('SU2Interface::Imposing inviscid boundary conditions.')
if self.getnDim() == 2:
velComponents = ['VELOCITY_X', 'VELOCITY_Y']
elif self.getnDim() == 3:
velComponents = ['VELOCITY_X', 'VELOCITY_Y', 'VELOCITY_Z']
else:
raise RuntimeError('su2Interface::Incorrect number of dimensions.')
velIdx = np.zeros(self.getnDim(), dtype=int)
for i, velComp in enumerate(velComponents):
velIdx[i] = self.solver.GetPrimitiveIndices()[velComp]
soundSpeedIdx = self.solver.GetPrimitiveIndices()["SOUND_SPEED"]
densityIdx = self.solver.GetPrimitiveIndices()["DENSITY"]
mach = np.sqrt(vel_x**2 + vel_y**2) / sound_speed
for body in self.iBnd:
for reg in body:
print('SU2Interface::Imposing inviscid boundary conditions on', reg.name)
for inode, nodeIndexFloat in enumerate(reg.nodesCoord[:,-1]):
nodeIndex = int(nodeIndexFloat)
# Velocity
for idim in range(self.getnDim()):
reg.V[inode,idim] = self.solver.Primitives().Get(int(nodeIndex), int(velIdx[idim]))
# Density
reg.Rho[inode] = self.solver.Primitives().Get(nodeIndex, densityIdx)
# Mach number
vel_norm = 0.0
for idim in range(self.getnDim()):
vel_norm += reg.V[inode,idim]**2
reg.M[inode] = np.sqrt(vel_norm) / self.solver.Primitives().Get(nodeIndex, soundSpeedIdx)
import matplotlib
matplotlib.use('Agg') # Use the Agg backend for non-interactive plotting
from matplotlib import pyplot as plt
plt.figure()
plt.plot(self.iBnd[0][0].nodesCoord[:,0], mach, 'x')
plt.xlabel('X Coordinate')
plt.ylabel('pressure')
plt.title('pressure Distribution')
plt.savefig('pressure_distribution.png') # Save the plot as an image file
print('Plot saved as mach_number_distribution.png')
quit()
if couplIter == 0:
self.cp0 = self.getCp()
# Retreive parameters.
for n in range(2):
for iRow, row in enumerate(self.iBnd[n].connectListNodes):
row=int(row)
self.iBnd[n].V[iRow,:] = self.solver.GetVertexVelocity(self.bndMarkers[n], row)
self.iBnd[n].M[iRow] = self.solver.GetVertexMach(self.bndMarkers[n], row)
self.iBnd[n].Rho[iRow] = self.solver.GetVertexDensity(self.bndMarkers[n], row)
if self.iBnd[n].Rho[iRow] == 0:
self.iBnd[n].Rho[iRow] = 1
self.setViscousSolver(couplIter)
plt.plot(self.iBnd[0][0].nodesCoord[:, 0], self.iBnd[0][0].M, 'x--')
plt.xlabel('$x$')
plt.ylabel('Mach number')
plt.savefig('Mach_number.png') # Save the plot as an image file
print('Plot saved as Mach_number.png')
def imposeBlwVel(self):
plt.figure()
plt.plot(self.iBnd[0][0].nodesCoord[:, 0], self.iBnd[0][0].Rho, 'x--')
plt.xlabel('$x$')
plt.ylabel('Density')
plt.savefig('Density.png') # Save the plot as an image file
print('Plot saved as Density.png')
plt.figure()
plt.plot(self.iBnd[0][0].nodesCoord[:, 0], self.iBnd[0][0].V[:,0], 'x--')
plt.xlabel('$x$')
plt.ylabel('Velocity X')
plt.savefig('Velocity_X.png') # Save the plot as an image file
print('Plot saved as Velocity_X.png')
plt.figure()
plt.plot(self.iBnd[0][0].nodesCoord[:, 0], self.iBnd[0][0].V[:,1], 'x--')
plt.xlabel('$x$')
plt.ylabel('Velocity Y')
plt.savefig('Velocity_Y.png') # Save the plot as an image file
print('Plot saved as Velocity_Y.png')
def setBlowingVelocity(self):
""" Extract the solution of the viscous calculation (blowing velocity)
and use it as a boundary condition in the inviscid solver.
"""
if self.iBnd[0].connectListNodes[0] != self.iBnd[0].connectListNodes[-1]:
raise RuntimeError("List of connectivity is incorrect.")
print('SU2Interface::Imposing blowing velocity boundary conditions.')
# interpolate blowing velocity from elements to nodes
from scipy.interpolate import interp1d
for body in self.iBnd:
for reg in body:
blowingNodes = interp1d(reg.elemsCoord[:,0], reg.blowingVel, kind='linear', fill_value='extrapolate')(reg.nodesCoord[:,0])
#Plot blowing velocity
import matplotlib
matplotlib.use('Agg') # Use the Agg backend for non-interactive plotting
from matplotlib import pyplot as plt
plt.figure()
plt.plot(reg.elemsCoord[:, 0], reg.blowingVel, 'o', label='Original')
plt.plot(reg.nodesCoord[:, 0], blowingNodes, 'x--', label='Interpolated')
plt.legend()
plt.xlabel('$x$')
plt.ylabel('Blowing velocity')
plt.savefig('Blowing_velocity.png') # Save the plot as an image file
print('Plot saved as Blowing_velocity.png')
self.getBlowingBoundary()
# Impose blowing velocities.
for n in range(len(self.bndMarkers)):
if n == 0:
self.iBnd[n].blowingVel[0] = .5 * (self.iBnd[n].blowingVel[0] + self.iBnd[n].blowingVel[-1])
self.iBnd[n].blowingVel = self.iBnd[n].blowingVel[self.iBnd[n].connectListElems[:-1].argsort()]
elif n == 1:
self.iBnd[1].blowingVel = np.insert(self.iBnd[1].blowingVel, 0, self.iBnd[0].blowingVel[0])
self.iBnd[n].blowingVel = self.iBnd[n].blowingVel[self.iBnd[n].connectListElems.argsort()]
for iVertex in range(len(self.iBnd[n].blowingVel)-1):
self.solver.SetBlowing(self.bndMarkers[n], iVertex, self.iBnd[n].blowingVel[iVertex])
def getWing(self):
if self.getnDim() == 2:
self._getAirfoil(0)
......@@ -208,22 +229,10 @@ class SU2Interface(SolversInterface):
""" Retreives the mesh used for SU2 Euler calculation.
Airfoil is already in seilig format (Upper TE -> Lower TE).
"""
# nodes = np.zeros((0,3))
# for arfTag in self.wingTags:
# nodesOnSide = np.empty((self.solver.GetNumberMarkerNodes(self.wingMarkersToID[arfTag]),3))
# for i in range(self.solver.GetNumberMarkerNodes(self.wingMarkersToID[arfTag])):
# vertex = self.solver.GetMarkerNode(self.wingMarkersToID[arfTag], i)
# nodesOnSide[i, :self.getnDim()] = self.solver.Coordinates().Get(vertex)
# nodesOnSide = np.flip(nodesOnSide, axis=0)
# nodes = np.row_stack((nodes, nodesOnSide))
# # Avoid duplicated leading edge point
# if len(np.where(nodes[:,0] == np.min(nodes[:,0]))[0]) > 1:
# nodes = np.delete(nodes, np.where(nodes[:,0] == np.min(nodes[:,0]))[0][1], axis=0)
elemsCoord = np.zeros((0,4))
nodesCoord = np.zeros((0,4))
for arfTag in self.wingTags:
nodesOnSide = np.empty((self.solver.GetNumberMarkerNodes(self.wingMarkersToID[arfTag]), 4))
nodesOnSide = np.zeros((self.solver.GetNumberMarkerNodes(self.wingMarkersToID[arfTag]), 4))
elemsCoordOnSide = np.zeros((self.solver.GetNumberMarkerElements(self.wingMarkersToID[arfTag]), 4))
inod = 0
for ielm in range(self.solver.GetNumberMarkerElements(self.wingMarkersToID[arfTag])):
......@@ -249,7 +258,8 @@ class SU2Interface(SolversInterface):
# Avoid duplicated leading edge point
if len(np.where(nodesCoord[:,0] == np.min(nodesCoord[:,0]))[0]) > 1:
nodesCoord = np.delete(nodesCoord, np.where(nodesCoord[:,0] == np.min(nodesCoord[:,0]))[0][1], axis=0)
for delindex in np.where(nodesCoord[:,0] == np.min(nodesCoord[:,0]))[0][1::-1]:
nodesCoord = np.delete(nodesCoord, delindex, axis=0)
import matplotlib
matplotlib.use('Agg') # Use the Agg backend for non-interactive plotting
......@@ -257,6 +267,7 @@ class SU2Interface(SolversInterface):
plt.figure()
plt.plot(nodesCoord[:, 0], nodesCoord[:, 1], 'x--')
plt.plot(elemsCoord[:, 0], elemsCoord[:, 1], 'o')
plt.xlim([0.95, 1.0])
plt.xlabel('X Coordinate')
plt.ylabel('Y Coordinate')
plt.title('Airfoil Nodes')
......@@ -268,6 +279,14 @@ class SU2Interface(SolversInterface):
raise RuntimeError('su2Interface::Airfoil has an open trailing edge.')
else:
print('Airfoil has closed trailing edge.')
# Check that the airfoil has one point at the leading edge
if len(np.where(nodesCoord[:,0] == np.min(nodesCoord[:,0]))[0]) > 1:
print(len(np.where(nodesCoord[:,0] == np.min(nodesCoord[:,0]))[0]))
raise RuntimeError('su2Interface::Airfoil has multiple points at the leading edge.')
else:
print(len(np.where(nodesCoord[:,0] == np.min(nodesCoord[:,0]))[0]))
print('Airfoil has one point at the leading edge.')
self.iBnd[ibody][0].initStructures(nodesCoord.shape[0], elemsCoord.shape[0])
self.iBnd[ibody][0].nodesCoord = nodesCoord
......
blast/tests/su2/config.cfg
+ 1
− 1
View file @ a040af32
......@@ -162,7 +162,7 @@ CONV_FIELD= RMS_DENSITY
CONV_RESIDUAL_MINVAL= -8
%
% Start convergence criteria at iteration number
CONV_STARTITER= 10
CONV_STARTITER= 0
%
% Number of elements to apply the criteria
CONV_CAUCHY_ELEMS= 100
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment