From a040af32d58cb84c00770585456f8f9ef422fae1 Mon Sep 17 00:00:00 2001
From: pdechamps <paul.dechamps@uliege.be>
Date: Sun, 23 Feb 2025 15:17:10 +0100
Subject: [PATCH] (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
---
blast/interfaces/su2/SU2Interface.py | 179 +++++++++++++++------------
blast/tests/su2/config.cfg | 2 +-
2 files changed, 100 insertions(+), 81 deletions(-)
diff --git a/blast/interfaces/su2/SU2Interface.py b/blast/interfaces/su2/SU2Interface.py
index 4d04e1a..2c89b10 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 75ed473..97fbd96 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
--
GitLab