From 64e26ad527ade27f86f213a21e25acb4f0816e1e Mon Sep 17 00:00:00 2001
From: pdechamps <paul.dechamps@uliege.be>
Date: Mon, 24 Feb 2025 02:56:42 +0100
Subject: [PATCH] (feat) Few modifs of su2 interface
---
blast/interfaces/su2/SU2Interface.py | 151 +++++++++++++++++++--------
1 file changed, 110 insertions(+), 41 deletions(-)
diff --git a/blast/interfaces/su2/SU2Interface.py b/blast/interfaces/su2/SU2Interface.py
index 2c89b10..c4c9832 100644
--- a/blast/interfaces/su2/SU2Interface.py
+++ b/blast/interfaces/su2/SU2Interface.py
@@ -18,6 +18,8 @@
# SU2 interface
# Paul Dechamps
import numpy as np
+import os
+import sys
from blast.interfaces.blSolversInterface import SolversInterface
import pysu2
@@ -25,6 +27,7 @@ import pysu2
class SU2Interface(SolversInterface):
def __init__(self, su2config, vconfig, task='analysis'):
self.filename = su2config.get('filename')
+ self.verbose = su2config.get('Verb', 0)
self.have_mpi, self.comm, self.status, self.myid = self.__setup_mpi()
self.solver = pysu2.CSinglezoneDriver(self.filename, 1, self.comm)
@@ -50,6 +53,10 @@ class SU2Interface(SolversInterface):
return have_mpi, comm, status, myid
def run(self):
+ print('SU2Interface::Running SU2 solver.')
+ #Remove output in the terminal
+ if self.getVerbose() == 0:
+ sys.stdout = open(os.devnull, 'w')
self.comm.barrier()
# run solver
self.solver.Preprocess(0)
@@ -59,6 +66,10 @@ class SU2Interface(SolversInterface):
self.solver.Monitor(0)
self.solver.Output(0)
self.comm.barrier()
+ #restore stdout
+ if self.getVerbose() == 0:
+ sys.stdout = sys.__stdout__
+ print(f'SU2Interface::SU2 solver finished in {self.solver.GetOutputValue("INNER_ITER"):.0f} iterations. RMS_DENSITY: {self.solver.GetOutputValue("RMS_DENSITY"):.2f}')
return 1
def writeCp(self, sfx=''):
@@ -117,7 +128,7 @@ class SU2Interface(SolversInterface):
return self.solver.Get_Mx()
def getVerbose(self):
- return 2
+ return self.verbose
def reset(self):
print('SU2Interface::Warning: Cannot reset solver.')
@@ -161,6 +172,18 @@ class SU2Interface(SolversInterface):
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)
+ self.iBnd[0][0].V[0,:] = self.iBnd[0][0].V[-2,:]
+ self.iBnd[0][0].V[-1,:] = self.iBnd[0][0].V[-2,:]
+ self.iBnd[0][0].Rho[0] = self.iBnd[0][0].Rho[-2]
+ self.iBnd[0][0].Rho[-1] = self.iBnd[0][0].Rho[-2]
+ self.iBnd[0][0].M[0] = self.iBnd[0][0].M[-2]
+ self.iBnd[0][0].M[-1] = self.iBnd[0][0].M[-2]
+ print('vel first', self.iBnd[0][0].V[0,:])
+ print('Rho first', self.iBnd[0][0].Rho[0])
+ print('M first', self.iBnd[0][0].M[0])
+ print('vel last', self.iBnd[0][0].V[-1,:])
+ print('Rho last', self.iBnd[0][0].Rho[-1])
+ print('M last', self.iBnd[0][0].M[-1])
import matplotlib
matplotlib.use('Agg') # Use the Agg backend for non-interactive plotting
from matplotlib import pyplot as plt
@@ -200,22 +223,44 @@ class SU2Interface(SolversInterface):
# 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])
+ blowingNodes = interp1d(self.iBnd[0][0].elemsCoord[:,0], self.iBnd[0][0].blowingVel, kind='linear', fill_value='extrapolate')(self.iBnd[0][0].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')
+
+ #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(self.iBnd[0][0].elemsCoord[:, 0], self.iBnd[0][0].blowingVel, 'o', label='Original')
+ plt.plot(self.iBnd[0][0].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')
+
+ 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]
+ # Get all the tags with the CHT option
+ 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]
+ normal = self.solver.GetMarkerVertexNormals(self.wingMarkersToID[arfTag], ivertex)
+ blw = np.zeros(3)
+ for idim in range(self.getnDim()):
+ blw += normal[idim] * blowingNodes[blasterIndex]
+ for idim in range(self.getnDim()):
+ print('Setting blowing velocity on node', nodeIndex, 'blasterIndex', blasterIndex, blw[idim])
+ self.solver.Primitives().Set(int(nodeIndex), int(velIdx[idim]), blw[idim])
def getWing(self):
if self.getnDim() == 2:
@@ -229,37 +274,52 @@ class SU2Interface(SolversInterface):
""" Retreives the mesh used for SU2 Euler calculation.
Airfoil is already in seilig format (Upper TE -> Lower TE).
"""
- elemsCoord = np.zeros((0,4))
- nodesCoord = np.zeros((0,4))
+ # elemsCoord = np.zeros((0,4))
+ # nodesCoord = np.zeros((0,4))
+ # for arfTag in self.wingTags:
+ # 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])):
+ # nodesonelm = self.solver.GetMarkerElementNodes(self.wingMarkersToID[arfTag], ielm)
+ # for n in nodesonelm:
+ # if n not in nodesOnSide[:,-1]:
+ # nodesOnSide[inod, :self.getnDim()] = self.solver.Coordinates().Get(n)
+ # nodesOnSide[inod, -1] = n
+ # inod += 1
+ # elemsCoordOnSide[ielm, :self.getnDim()] += self.solver.Coordinates().Get(n)
+ # elemsCoordOnSide[ielm, :self.getnDim()] /= len(nodesonelm)
+ # elemsCoordOnSide[ielm, -1] = self.solver.GetMarkerElementGlobalIndex(self.wingMarkersToID[arfTag], ielm)
+
+ # # Sort in ascending order of first column
+ # nodesOnSide = nodesOnSide[nodesOnSide[:,0].argsort()]
+ # elemsCoordOnSide = elemsCoordOnSide[elemsCoordOnSide[:,0].argsort()]
+ # if arfTag == 'airfoil':
+ # nodesOnSide = np.flip(nodesOnSide, axis=0)
+ # elemsCoordOnSide = np.flip(elemsCoordOnSide, axis=0)
+ # #elemsCoordOnSide = np.flip(elemsCoordOnSide, axis=0)
+ # nodesCoord = np.row_stack((nodesCoord, nodesOnSide))
+ # elemsCoord = np.row_stack((elemsCoord, elemsCoordOnSide))
+
+ nodesCoord = np.zeros((0, 4))
for arfTag in self.wingTags:
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])):
- nodesonelm = self.solver.GetMarkerElementNodes(self.wingMarkersToID[arfTag], ielm)
- for n in nodesonelm:
- if n not in nodesOnSide[:,-1]:
- nodesOnSide[inod, :self.getnDim()] = self.solver.Coordinates().Get(n)
- nodesOnSide[inod, -1] = n
- inod += 1
- elemsCoordOnSide[ielm, :self.getnDim()] += self.solver.Coordinates().Get(n)
- elemsCoordOnSide[ielm, :self.getnDim()] /= len(nodesonelm)
- elemsCoordOnSide[ielm, -1] = self.solver.GetMarkerElementGlobalIndex(self.wingMarkersToID[arfTag], ielm)
-
+ for i in range(self.solver.GetNumberMarkerNodes(self.wingMarkersToID[arfTag])):
+ nodeIndex = self.solver.GetMarkerNode(self.wingMarkersToID[arfTag], i)
+ nodesOnSide[i, :self.getnDim()] = self.solver.Coordinates().Get(nodeIndex)
+ nodesOnSide[i, -1] = nodeIndex
# Sort in ascending order of first column
nodesOnSide = nodesOnSide[nodesOnSide[:,0].argsort()]
- elemsCoordOnSide = elemsCoordOnSide[elemsCoordOnSide[:,0].argsort()]
if arfTag == 'airfoil':
nodesOnSide = np.flip(nodesOnSide, axis=0)
- elemsCoordOnSide = np.flip(elemsCoordOnSide, axis=0)
- #elemsCoordOnSide = np.flip(elemsCoordOnSide, axis=0)
+ elif arfTag == 'airfoil_':
+ nodesOnSide = np.delete(nodesOnSide, 0, axis=0)
nodesCoord = np.row_stack((nodesCoord, nodesOnSide))
- elemsCoord = np.row_stack((elemsCoord, elemsCoordOnSide))
-
- # Avoid duplicated leading edge point
- if len(np.where(nodesCoord[:,0] == np.min(nodesCoord[:,0]))[0]) > 1:
- for delindex in np.where(nodesCoord[:,0] == np.min(nodesCoord[:,0]))[0][1::-1]:
- nodesCoord = np.delete(nodesCoord, delindex, axis=0)
+
+ elemsCoord = np.zeros((nodesCoord.shape[0]-1, 4))
+ for i in range(nodesCoord.shape[0]-1):
+ elemsCoord[i, :self.getnDim()] = (nodesCoord[i, :self.getnDim()] + nodesCoord[i+1, :self.getnDim()]) / 2
+ elemsCoord[i, -1] = i
import matplotlib
matplotlib.use('Agg') # Use the Agg backend for non-interactive plotting
@@ -267,7 +327,8 @@ 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.xlim([-0.001, 0.002])
+ #plt.ylim([-0.02, 0.02])
plt.xlabel('X Coordinate')
plt.ylabel('Y Coordinate')
plt.title('Airfoil Nodes')
@@ -287,6 +348,14 @@ class SU2Interface(SolversInterface):
else:
print(len(np.where(nodesCoord[:,0] == np.min(nodesCoord[:,0]))[0]))
print('Airfoil has one point at the leading edge.')
+
+ # Check for duplicated elements
+ if len(np.unique(elemsCoord[:,-1])) != elemsCoord.shape[0]:
+ raise RuntimeError('su2Interface::Duplicated elements in airfoil mesh.')
+ else:
+ print('No duplicated elements in airfoil mesh.')
+
+ print('nNodes:', nodesCoord.shape[0], 'nElems:', elemsCoord.shape[0])
self.iBnd[ibody][0].initStructures(nodesCoord.shape[0], elemsCoord.shape[0])
self.iBnd[ibody][0].nodesCoord = nodesCoord
--
GitLab