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

(feat) Few modifs of su2 interface

parent a040af32
No related branches found
No related tags found
No related merge requests found
Pipeline #52593 failed
Hide whitespace changes
Inline Side-by-side
blast/interfaces/su2/SU2Interface.py
+ 110
− 41
View file @ 64e26ad5
......@@ -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
......
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