(feat) SU2 iterating

However the solution seems only slightly affected by the blowing velocity computed.

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

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

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