diff --git a/blast/interfaces/su2/SU2Interface.py b/blast/interfaces/su2/SU2Interface.py
index c4c9832e13d448fe98f9b65c14f09964d6d6f92a..b0d4fd19a8a615cfffe3048d3f326612a12ee665 100644
--- a/blast/interfaces/su2/SU2Interface.py
+++ b/blast/interfaces/su2/SU2Interface.py
@@ -30,6 +30,9 @@ class SU2Interface(SolversInterface):
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)
+ self.comm.barrier()
+ # run solver
+ self.solver.Preprocess(0)
self.wingMarkersToID = {key: i for i, key in enumerate(self.solver.GetMarkerTags())}
self.wingTags = su2config.get('wingTags')
@@ -53,36 +56,70 @@ 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)
+
self.solver.Run()
self.solver.Postprocess()
# write outputs
self.solver.Monitor(0)
self.solver.Output(0)
self.comm.barrier()
+ print(f'SU2Interface::SU2 solver finished in {self.solver.GetOutputValue("INNER_ITER"):.0f} iterations. RMS_DENSITY: {self.solver.GetOutputValue("RMS_DENSITY"):.2f}')
#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=''):
""" Write Cp distribution around the airfoil on file.
"""
- pass
- def getCp(self):
- cp = np.zeros((self.nPoints[0], 2))
- for i in range(len(self.iBnd[0].nodesCoord[:,0])):
- cp[i, 0] = self.iBnd[0].nodesCoord[i,0]
- cp[i, 1] = self.solver.GetVertexCp(int(self.bndMarkers[0]), int(self.iBnd[0].connectListNodes[i]))
- return cp
+ #GEt farfield pressure, velocity and density
+ vel_x_idx = self.solver.GetPrimitiveIndices()["VELOCITY_X"]
+ vel_y_idx = self.solver.GetPrimitiveIndices()["VELOCITY_Y"]
+ pressure_idx = self.solver.GetPrimitiveIndices()["PRESSURE"]
+ density_idx = self.solver.GetPrimitiveIndices()["DENSITY"]
+
+ nNodes_farfield = self.solver.GetNumberMarkerNodes(self.wingMarkersToID['farfield'])
+ velocity_farfield = np.zeros(nNodes_farfield)
+ pressure_farfield = np.zeros(nNodes_farfield)
+ density_farfield = np.zeros(nNodes_farfield)
+ primitives = self.solver.Primitives()
+ for inode in range(nNodes_farfield):
+ nodeIndex = self.solver.GetMarkerNode(self.wingMarkersToID['farfield'], inode)
+ # Velocity
+ vel_x = primitives.Get(nodeIndex, vel_x_idx)
+ vel_y = primitives.Get(nodeIndex, vel_y_idx)
+ velocity_farfield[inode] = np.sqrt(vel_x**2 + vel_y**2)
+ # Pressure
+ pressure_farfield[inode] = primitives.Get(nodeIndex, pressure_idx)
+ # Density
+ density_farfield[inode] = primitives.Get(nodeIndex, density_idx)
+ ref_pressure = np.mean(pressure_farfield)
+ ref_velocity = np.mean(velocity_farfield)
+ ref_density = np.mean(density_farfield)
+
+ cp = np.empty(0, dtype=float)
+ for body in self.iBnd:
+ for reg in body:
+ for inode, nodeIndexFloat in enumerate(reg.nodesCoord[:,-1]):
+ nodeIndex = int(nodeIndexFloat)
+ pressure = self.solver.Primitives().Get(nodeIndex, pressure_idx)
+ cp = np.append(cp, (pressure - ref_pressure)/(1/2*ref_density*ref_velocity**2))
+ np.savetxt(f'Cp{sfx}.dat', np.column_stack((self.iBnd[0][0].nodesCoord[:, 0], cp)))
+
+ # from matplotlib import pyplot as plt
+ # import matplotlib
+ # matplotlib.use('Agg') # Use the Agg backend for non-interactive plotting
+ # plt.figure()
+ # plt.plot(self.iBnd[0][0].nodesCoord[:, 0], cp, 'x--')
+ # plt.xlabel('$x$')
+ # plt.ylabel('$c_p$')
+ # plt.gca().invert_yaxis()
+ # plt.savefig(f'cp{sfx}.png') # Save the plot as an image file
+ # print(f'Plot saved as cp{sfx}.png')
def save(self, writer):
pass
@@ -131,7 +168,8 @@ class SU2Interface(SolversInterface):
return self.verbose
def reset(self):
- print('SU2Interface::Warning: Cannot reset solver.')
+ if self.getVerbose() > 0:
+ print('SU2Interface::Warning: Cannot reset solver.')
pass
def getInviscidBC(self):
@@ -143,7 +181,7 @@ class SU2Interface(SolversInterface):
- couplIter (int): Coupling iteration.
"""
- print('SU2Interface::Imposing inviscid boundary conditions.')
+ # print('SU2Interface::Imposing inviscid boundary conditions.')
if self.getnDim() == 2:
velComponents = ['VELOCITY_X', 'VELOCITY_Y']
elif self.getnDim() == 3:
@@ -159,7 +197,7 @@ class SU2Interface(SolversInterface):
for body in self.iBnd:
for reg in body:
- print('SU2Interface::Imposing inviscid boundary conditions on', reg.name)
+ #print('SU2Interface::Imposing inviscid boundary conditions on', reg.name)
for inode, nodeIndexFloat in enumerate(reg.nodesCoord[:,-1]):
nodeIndex = int(nodeIndexFloat)
# Velocity
@@ -178,78 +216,67 @@ class SU2Interface(SolversInterface):
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
- plt.figure()
- 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')
+ # 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
+ # plt.figure()
+ # 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')
- 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].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[:,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')
+ # 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.
"""
- print('SU2Interface::Imposing blowing velocity boundary conditions.')
+ # print('SU2Interface::Imposing blowing velocity boundary conditions.')
# interpolate blowing velocity from elements to nodes
from scipy.interpolate import interp1d
blowingNodes = interp1d(self.iBnd[0][0].elemsCoord[:,0], self.iBnd[0][0].blowingVel, kind='linear', fill_value='extrapolate')(self.iBnd[0][0].nodesCoord[:,0])
-
+ blowingNodes[abs(blowingNodes) > 0.5] = 1e-6
#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')
+ # 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)
@@ -258,9 +285,7 @@ class SU2Interface(SolversInterface):
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])
+ self.solver.SetBlowingVelocity(nodeIndex, blw[0], blw[1], blw[2])
def getWing(self):
if self.getnDim() == 2:
@@ -321,19 +346,19 @@ class SU2Interface(SolversInterface):
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
- from matplotlib import pyplot as plt
- plt.figure()
- plt.plot(nodesCoord[:, 0], nodesCoord[:, 1], 'x--')
- plt.plot(elemsCoord[:, 0], elemsCoord[:, 1], 'o')
- #plt.xlim([-0.001, 0.002])
- #plt.ylim([-0.02, 0.02])
- plt.xlabel('X Coordinate')
- plt.ylabel('Y Coordinate')
- plt.title('Airfoil Nodes')
- plt.savefig('airfoil_nodes.png') # Save the plot as an image file
- print('Plot saved as airfoil_nodes.png')
+ # import matplotlib
+ # matplotlib.use('Agg') # Use the Agg backend for non-interactive plotting
+ # from matplotlib import pyplot as plt
+ # plt.figure()
+ # plt.plot(nodesCoord[:, 0], nodesCoord[:, 1], 'x--')
+ # plt.plot(elemsCoord[:, 0], elemsCoord[:, 1], 'o')
+ # #plt.xlim([-0.001, 0.002])
+ # #plt.ylim([-0.02, 0.02])
+ # plt.xlabel('X Coordinate')
+ # plt.ylabel('Y Coordinate')
+ # plt.title('Airfoil Nodes')
+ # plt.savefig('airfoil_nodes.png') # Save the plot as an image file
+ # print('Plot saved as airfoil_nodes.png')
# Check that the airfoil has closed trailing edge
if np.any(nodesCoord[0, [0,self.getnDim()-1]] != nodesCoord[-1, [0,self.getnDim()-1]]):
@@ -355,7 +380,7 @@ class SU2Interface(SolversInterface):
else:
print('No duplicated elements in airfoil mesh.')
- print('nNodes:', nodesCoord.shape[0], 'nElems:', elemsCoord.shape[0])
+ # 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
diff --git a/blast/tests/su2/blisu2.py b/blast/tests/su2/blisu2.py
index f8383801b5773cdbab23603388933ac40b6bdecc..99ff6ad787ebe016075302d527e116de9d68819c 100644
--- a/blast/tests/su2/blisu2.py
+++ b/blast/tests/su2/blisu2.py
@@ -17,27 +17,28 @@
# @author Paul Dechamps
-# @date 2022
+# @date 2022 reworked 2025
# Test the vii implementation using SU2 as the inviscid solver
# Imports.
-from blast.blCoupler import Coupler
import blast.blUtils as vutils
+from fwk.wutils import parseargs
-def cfgSu2():
+def cfgSu2(verb):
return {
'filename' : '/home/c130/lab/blaster/blast/tests/su2/config.cfg',
'wingTags' : ['airfoil', 'airfoil_'],
+ 'Verb': verb, # Verbosity level
}
def cfgBlast():
return {
'Re' : 1e7, # Freestream Reynolds number
'CFL0' : 1, # Inital CFL number of the calculation
- 'couplIter': 50, # Maximum number of iterations
+ 'couplIter': 20, # Maximum number of iterations
'sections': {'airfoil': [0.0]},
'spans': {'airfoil': 1.0},
- 'couplTol' : 1e-4, # Tolerance of the VII methodology
+ 'couplTol' : 1e-6, # Tolerance of the VII methodology
'iterPrint': 1, # int, number of iterations between outputs
'resetInv' : True, # bool, flag to reset the inviscid calculation at every iteration.
'xtrF' : [None, 0.4], # Forced transition locations
@@ -45,12 +46,23 @@ def cfgBlast():
}
def main():
-
- icfg = cfgSu2()
+ args = parseargs()
+ icfg = cfgSu2(args.verb)
vcfg = cfgBlast()
coupler, isol, vsol = vutils.initBlast(icfg, vcfg, iSolver='SU2', task='analysis')
coupler.run()
print('ok main')
+
+ import numpy as np
+ cpi = np.loadtxt('Cp_inviscid.dat')
+ cpCorrected = np.loadtxt('Cp_viscous.dat')
+ import matplotlib.pyplot as plt
+ import matplotlib
+ matplotlib.use('Agg')
+ plt.plot(cpCorrected[:,0], cpCorrected[:,1], lw=3)
+ plt.plot(cpi[:,0], cpi[:,1], '--', lw=3)
+ plt.gca().invert_yaxis()
+ plt.savefig('cp.png')
quit()
# Extract solution
diff --git a/blast/tests/su2/config.cfg b/blast/tests/su2/config.cfg
index 97fbd96482020599d4ac9ba43b38fe5c5c7e1eed..de33885fcddd1ed995b39e9fc1276309b1eb2c79 100644
--- a/blast/tests/su2/config.cfg
+++ b/blast/tests/su2/config.cfg
@@ -58,9 +58,10 @@ REF_AREA= 1.0
% Euler wall boundary marker(s) (NONE = no marker)
MARKER_EULER= ( airfoil, airfoil_ )
%
+MARKER_MOVING= ( airfoil, airfoil_ )
+%
SURFACE_MOVEMENT= (MOVING_WALL, MOVING_WALL)
%
-MARKER_MOVING= ( airfoil, airfoil_ )
% Motion mach number (non-dimensional). Used for initializing a viscous flow
% with the Reynolds number and for computing force coeffs. with dynamic meshes.
%MACH_MOTION= 0.
@@ -84,12 +85,9 @@ MARKER_MONITORING= ( airfoil, airfoil_ )
%
% Numerical method for spatial gradients (GREEN_GAUSS, WEIGHTED_LEAST_SQUARES)
NUM_METHOD_GRAD= WEIGHTED_LEAST_SQUARES
-
-SLOPE_LIMITER_FLOW= VENKATAKRISHNAN_WANG
-VENKAT_LIMITER_COEFF= 0.1
%
% Courant-Friedrichs-Lewy condition of the finest grid
-CFL_NUMBER= 1e3
+CFL_NUMBER= 1e2
%
% Adaptive CFL number (NO, YES)
CFL_ADAPT= NO
@@ -99,11 +97,11 @@ CFL_ADAPT= NO
CFL_ADAPT_PARAM= ( 0.1, 5.0, 50.0, 1e10 )
%
% Runge-Kutta alpha coefficients
-%RK_ALPHA_COEFF= ( 0.66667, 0.66667, 1.000000 )
+RK_ALPHA_COEFF= ( 0.66667, 0.66667, 1.000000 )
%
% Number of total iterations
-ITER= 999999
-
+ITER= 1000
+%
% ------------------------ LINEAR SOLVER DEFINITION ---------------------------%
%
% Linear solver for implicit formulations (BCGSTAB, FGMRES)
@@ -117,7 +115,7 @@ LINEAR_SOLVER_ERROR= 1E-6
%
% Max number of iterations of the linear solver for the implicit formulation
LINEAR_SOLVER_ITER= 10
-
+%
% -------------------------- MULTIGRID PARAMETERS -----------------------------%
%
% Multi-Grid Levels (0 = no multi-grid)
@@ -148,11 +146,21 @@ MG_DAMP_PROLONGATION= 1.0
CONV_NUM_METHOD_FLOW= ROE
%
% 2nd and 4th order artificial dissipation coefficients
-%JST_SENSOR_COEFF= ( 0.01, 0.01 )
+%JST_SENSOR_COEFF= ( 0.5, 0.02 )
+%LAX_SENSOR_COEFF= 0.001
%
% Time discretization (RUNGE-KUTTA_EXPLICIT, EULER_IMPLICIT, EULER_EXPLICIT)
TIME_DISCRE_FLOW= EULER_IMPLICIT
-
+%
+% ----------------------- SLOPE LIMITER DEFINITION ----------------------------%
+%
+MUSCL_FLOW= YES
+%
+SLOPE_LIMITER_FLOW= VENKATAKRISHNAN
+%
+% Coefficient for the limiter
+VENKAT_LIMITER_COEFF= 0.005
+%
% --------------------------- CONVERGENCE PARAMETERS --------------------------%
%
% Convergence criteria (CAUCHY, RESIDUAL)
@@ -172,11 +180,11 @@ CONV_CAUCHY_EPS= 1E-10
% ------------------------- INPUT/OUTPUT INFORMATION --------------------------%
%
-SCREEN_WRT_FREQ_INNER= 5
+SCREEN_WRT_FREQ_INNER= 200
VOLUME_OUTPUT = SOLUTION, PRIMITIVE, RMS_DENSITY
% Mesh input file
-MESH_FILENAME= /home/c130/lab/blaster/blast/models/su2/n0012_su2.su2
+MESH_FILENAME= /home/c130/lab/blaster/blast/models/su2/n0012_su2_2.su2
%
% Mesh input file format (SU2, CGNS, NETCDF_ASCII)
MESH_FORMAT= SU2