diff --git a/blast/interfaces/su2/SU2Interface.py b/blast/interfaces/su2/SU2Interface.py
index 4d04e1af3485b6104c14c458efc31c3cff143202..2c89b103b53443d75d85cdc01b881df67c0e84c0 100644
--- a/blast/interfaces/su2/SU2Interface.py
+++ b/blast/interfaces/su2/SU2Interface.py
@@ -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
diff --git a/blast/tests/su2/config.cfg b/blast/tests/su2/config.cfg
index 75ed473754e4f810d29680617a226d9fe415e5e4..97fbd96482020599d4ac9ba43b38fe5c5c7e1eed 100644
--- a/blast/tests/su2/config.cfg
+++ b/blast/tests/su2/config.cfg
@@ -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