diff --git a/blast/interfaces/su2/blSU2Interface.py b/blast/interfaces/su2/blSU2Interface.py
index 816ae48f7a660026dcdd0546ba9058f9e6ee47fa..0d3a5e461fd6158c32f5c66e9d2beba8ccc19ddc 100644
--- a/blast/interfaces/su2/blSU2Interface.py
+++ b/blast/interfaces/su2/blSU2Interface.py
@@ -168,14 +168,18 @@ class SU2Interface(SolversInterface):
""" Extract the inviscid paramaters at the edge of the boundary layer
and use it as a boundary condition in the viscous solver.
"""
- if self.getnDim() == 2:
+ region = self.iBnd[0][0]
+ wingRows = region.getNodeRows()
+ nNodes, ndim = region.getnNodes(), self.getnDim()
+
+ if ndim == 2:
velComponents = ['VELOCITY_X', 'VELOCITY_Y']
- elif self.getnDim() == 3:
+ elif ndim == 3:
velComponents = ['VELOCITY_X', 'VELOCITY_Y', 'VELOCITY_Z']
else:
raise RuntimeError('su2Interface::Incorrect number of dimensions.')
- velIdx = np.zeros(self.getnDim(), dtype=int)
+ velIdx = np.zeros(ndim, dtype=int)
for i, velComp in enumerate(velComponents):
velIdx[i] = self.solver.GetPrimitiveIndices()[velComp]
soundSpeedIdx = self.solver.GetPrimitiveIndices()["SOUND_SPEED"]
@@ -183,18 +187,117 @@ class SU2Interface(SolversInterface):
for body in self.iBnd:
for reg in body:
- for inode, nodeIndexFloat in enumerate(reg.nodesCoord[:,-1]):
+ V = np.empty((nNodes, ndim))
+ M = np.empty(nNodes)
+ Rho = np.empty(nNodes)
+ for inode, nodeIndexFloat in enumerate(wingRows):
nodeIndex = int(nodeIndexFloat)
# Velocity
- for idim in range(self.getnDim()):
- reg.V[inode,idim] = self.solver.Primitives().Get(int(nodeIndex), int(velIdx[idim]))
+ for idim in range(ndim):
+ V[inode,idim] = self.solver.Primitives().Get(int(nodeIndex), int(velIdx[idim]))
# Density
- reg.Rho[inode] = self.solver.Primitives().Get(nodeIndex, densityIdx)
+ 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)
+ for idim in range(ndim):
+ vel_norm += V[inode,idim]**2
+ M[inode] = np.sqrt(vel_norm) / self.solver.Primitives().Get(nodeIndex, soundSpeedIdx)
+ reg.updateVariables(M, V, Rho)
+
+ from matplotlib import pyplot as plt
+ plt.plot(reg._nodesCoord[:,0], reg._M, 'x-')
+
+ # Reflect the last two nodes
+ nodes_ = region.getNodes().copy()
+ M_ = M.copy()
+ new_nodes, new_M = self.orthogonal_reflection(nodes_, M_)
+ #plt.plot(new_nodes[:,0], new_M, 'o-')
+ #plt.show()
+ smoothed_nodes, smoothed_M = self.explicit_smoother(new_nodes, new_M, iterations=10)
+ plt.plot(smoothed_nodes[:,0], smoothed_M, 'o-')
+
+ plt.show()
+ quit()
+ plt.show()
+
+ def orthogonal_reflection(self, nodes, M, axis_x=1):
+ """
+ Reflect the first two and last two points of the Mach number data (simulating wake) with respect
+ to the given axis (default trailing edge at x=1), and append them to the data.
+
+ Args:
+ - nodes: 2D array where each row contains [x, y, z] coordinates.
+ - M: 1D array containing the Mach numbers.
+ - axis_x: The x-coordinate about which to reflect (default is 1, the trailing edge).
+
+ Returns:
+ - new_nodes: Reflected nodes including the additional wake simulation points.
+ - new_M: Reflected Mach numbers.
+ """
+ # Reflect the first two nodes and Mach numbers
+ first_two_nodes = nodes[1:5] # First two [x, y, z]
+ first_two_M = M[1:5] # First two Mach numbers
+
+ # Orthogonal reflection: Reflect the x-values with respect to x=axis_x
+ reflected_first_two_nodes = np.copy(first_two_nodes)
+ reflected_first_two_nodes[:, 0] = 2*axis_x - first_two_nodes[:, 0] # Reflect x-values
+
+ # Reflect the Mach numbers
+ reflected_first_two_M = first_two_M[::-1] # Reverse the order of the first two Mach numbers
+
+ # Reflect the last two nodes and Mach numbers
+ last_two_nodes = nodes[-5:-1] # Last two [x, y, z]
+ last_two_M = M[-5:-1] # Last two Mach numbers
+
+ # Orthogonal reflection: Reflect the x-values with respect to x=axis_x
+ reflected_last_two_nodes = np.copy(last_two_nodes)
+ reflected_last_two_nodes[:, 0] = 2*axis_x - last_two_nodes[:, 0] # Reflect x-values
+
+ # Reflect the Mach numbers
+ reflected_last_two_M = last_two_M[::-1] # Reverse the order of the last two Mach numbers
+
+ # Append the reflected points to the original nodes and Mach numbers
+ new_nodes = np.vstack((reflected_first_two_nodes, nodes, reflected_last_two_nodes))
+ new_M = np.concatenate((reflected_first_two_M, M, reflected_last_two_M))
+
+ return new_nodes, new_M
+
+ def explicit_smoother(self, nodes, M, mask_range=(0.8, 1.0), iterations=5):
+ """
+ Smooth the Mach number data only for nodes in a given range (between mask_range[0] and mask_range[1]).
+ The smoothing uses explicit finite difference method with Neumann boundary conditions.
+
+ Args:
+ - nodes: 2D array where each row contains [x, y, z] coordinates.
+ - M: 1D array containing the Mach numbers.
+ - mask_range: Tuple (start_pct, end_pct) defining the range of x-values to apply smoothing.
+ - iterations: Number of smoothing iterations to perform.
+
+ Returns:
+ - smoothed_nodes: Original nodes (unchanged).
+ - smoothed_M: Smoothed Mach number data.
+ """
+ # Create a mask for the region between 0.8 and 1.0 of the chord
+ mask = (nodes[:, 0] > mask_range[0]) & (nodes[:, 0] < mask_range[1])
+
+ smoothed_M = M.copy() # Copy the Mach numbers for smoothing
+
+ for _ in range(iterations):
+ new_M = smoothed_M.copy()
+
+ # Apply smoothing only to the values in the specified range (using the mask)
+ for i in range(1, len(M) - 1):
+ if mask[i]: # Only smooth where the mask is True
+ new_M[i] = 0.5 * (smoothed_M[i - 1] + smoothed_M[i + 1])
+
+ te_index = np.argmax(nodes[:, 0]) # Index of the trailing edge
+ # Neumann boundary condition: no change at the first and last points
+ new_M[0] = smoothed_M[1] # Set first point equal to the second point
+ new_M[-1] = smoothed_M[-2] # Set last point equal to the second-to-last point
+
+ smoothed_M = new_M
+
+ return nodes, smoothed_M
def setBlowingVelocity(self):
"""
@@ -211,41 +314,47 @@ class SU2Interface(SolversInterface):
be projected into the global frame and interpolated to the nodes.
"""
# Compute blowing velocity in the global frame of reference on the elements
- blw_elems = np.zeros((self.iBnd[0][0].elemsCoord.shape[0], self.getnDim()), dtype=float)
+ region = self.iBnd[0][0]
+ nNodes, nElms, ndim = region.getnNodes(), region.getnElms(), self.getnDim()
+ rowsWing = region.getNodeRows()
+ blowingWing = region.getBlowingVelocity()
+ nodesWing = region.getNodes()
+
+ blw_elems_normal = np.zeros((nElms, ndim), dtype=float)
normals = []
- for ielm in range(len(self.iBnd[0][0].elemsCoord)):
+ for ielm in range(len(self.iBnd[0][0].getElms('all'))):
nod1 = ielm
nod2 = ielm + 1
- x1 = self.iBnd[0][0].nodesCoord[nod1,:self.getnDim()]
- x2 = self.iBnd[0][0].nodesCoord[nod2,:self.getnDim()]
+ x1 = nodesWing[nod1,:ndim]
+ x2 = nodesWing[nod2,:ndim]
# Components of the tangent vector
- t = np.empty(self.getnDim())
- for idim in range(self.getnDim()):
+ t = np.empty(ndim)
+ for idim in range(ndim):
t[idim] = x2[idim] - x1[idim]
normt = np.linalg.norm(t)
# Components of the normal vector
- n = np.empty(self.getnDim())
- if self.getnDim() == 2:
+ n = np.empty(ndim)
+ if ndim == 2:
# 90° rotation
n[0] = t[1] / normt
n[1] = -t[0] / normt
- elif self.getnDim() == 3:
+ elif ndim == 3:
# Compute using cross product with another edge
raise RuntimeError('3D normal not implemented yet (use SU2 function if possible).')
normals.append(n)
- for idim in range(self.getnDim()):
- blw_elems[ielm, idim] = self.iBnd[0][0].blowingVel[ielm] * n[idim]
+ for idim in range(ndim):
+ blw_elems_normal[ielm, idim] = blowingWing[ielm] * n[idim]
# Compute blowing velocity at nodes
- blw_nodes = np.zeros((self.iBnd[0][0].nodesCoord.shape[0], 3), dtype=float)
- for inode in range(self.iBnd[0][0].nodesCoord.shape[0]):
+ blw_nodes = np.zeros((nNodes, 3), dtype=float)
+ for inode in range(nNodes):
# Look for the elements sharing node inode
# At TE, elm1 is the last element of upper side and elm2 is the last element of lower side.
- if inode == 0 or inode == self.iBnd[0][0].nodesCoord.shape[0]-1:
+ if inode == 0 or inode == nNodes-1:
elm1 = 0
elm2 = -1
else:
@@ -254,23 +363,26 @@ class SU2Interface(SolversInterface):
# Because the blowing velocity on the lower side points downwards
sign = 1
- if inode > self.iBnd[0][0].nodesCoord.shape[0]//2:
+ if inode > nNodes//2:
sign = -1
# The blowing velocity on one node is the average
# of the blowing velocity on the elements sharing the node
- for idim in range(self.getnDim()):
- blw_nodes[inode, idim] = sign * 0.5*(blw_elems[elm1, idim] + blw_elems[elm2, idim])
+ for idim in range(ndim):
+ blw_nodes[inode, idim] = sign * 0.5*(blw_elems_normal[elm1, idim] + blw_elems_normal[elm2, idim])
- blw_nodes[blw_nodes < -0.01] = -0.01
- blw_nodes[blw_nodes > 0.01] = 0.01
+ # blw_nodes[blw_nodes < -0.01] = -0.01
+ # blw_nodes[blw_nodes > 0.01] = 0.01
# Set blowing velocity in SU2 solver
for arfTag in self.wingTags:
for ivertex in range(self.solver.GetNumberMarkerNodes(self.wingMarkersToID[arfTag])):
nodeIndex = self.solver.GetMarkerNode(self.wingMarkersToID[arfTag], ivertex)
- blasterIndex = np.where(self.iBnd[0][0].nodesCoord[:, -1] == nodeIndex)[0][0]
- self.solver.SetBlowingVelocity(nodeIndex, blw_nodes[blasterIndex, 0], blw_nodes[blasterIndex, 1], blw_nodes[blasterIndex, 2])
+ blasterIndex = np.where(rowsWing == nodeIndex)[0][0]
+ self.solver.SetBlowingVelocity(nodeIndex,
+ blw_nodes[blasterIndex, 0],
+ blw_nodes[blasterIndex, 1],
+ blw_nodes[blasterIndex, 2])
def getWing(self):
if self.getnDim() == 2:
@@ -340,11 +452,14 @@ class SU2Interface(SolversInterface):
raise RuntimeError('su2Interface::Duplicated elements in airfoil mesh.')
# Fill the data structures
- self.iBnd[ibody][0].initStructures(nodesCoord.shape[0], elemsCoord.shape[0])
- self.iBnd[ibody][0].nodesCoord = nodesCoord
- self.iBnd[ibody][0].elemsCoord = elemsCoord
- self.iBnd[ibody][0].setConnectList(np.array(nodesCoord[:,-1],dtype=int),
+ self.iBnd[ibody][0].setMesh(nodesCoord[:,:3], elemsCoord[:,:3])
+ #self.iBnd[ibody][0].initStructures(nodesCoord.shape[0], elemsCoord.shape[0])
+ #self.iBnd[ibody][0].nodesCoord = nodesCoord
+ #self.iBnd[ibody][0].elemsCoord = elemsCoord
+ self.iBnd[ibody][0].setConnectLists(np.array(nodesCoord[:,-1],dtype=int),
np.array(elemsCoord[:,-1],dtype=int))
+ self.iBnd[ibody][0].setRows(np.array(nodesCoord[:,-1],dtype=int),
+ np.array(elemsCoord[:,-1],dtype=int))
def _modify_meshpath(self, config, mesh):
# Open config file and look for line MESH_FILENAME
diff --git a/blast/tests/su2/t_n0012_su2_2D.py b/blast/tests/su2/t_n0012_su2_2D.py
index fc1d0514ff5dbbbefcbf3e14ee52bf59e29d461e..c38672dac27cdefe7468777804124bf1f41f31fc 100644
--- a/blast/tests/su2/t_n0012_su2_2D.py
+++ b/blast/tests/su2/t_n0012_su2_2D.py
@@ -50,13 +50,12 @@ def cfgBlast():
}
def main():
- return 0
args = parseargs()
icfg = cfgSu2(args.verb)
vcfg = cfgBlast()
coupler, isol, vsol = vutils.initBlast(icfg, vcfg, iSolver='SU2', task='analysis')
- #coupler.run()
- isol.run()
+ coupler.run()
+ #isol.run()
print(ccolors.ANSI_BLUE + 'PyTesting...' + ccolors.ANSI_RESET)
tests = CTests()