(feat) Impose blowing velocity correctly in SU2

Grid vel in SU2 has to be a vector defined on the nodes in the global frame of reference. The blowing velocity computed in BLASTER is a scalar quantity computed on the elements. It has to be interpolated and projected.

blast/interfaces/su2/SU2Interface.py

@@ -30,8 +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)
+
+        # Prepare solver
         self.comm.barrier()
-        # run solver
         self.solver.Preprocess(0)
 
         self.wingMarkersToID = {key: i for i, key in enumerate(self.solver.GetMarkerTags())}
@@ -54,24 +55,25 @@ class SU2Interface(SolversInterface):
             status = None
             myid = 0
         return have_mpi, comm, status, myid
-    
+
     def run(self):
         #Remove output in the terminal
-        if self.getVerbose() == 0:
+        if self.getVerbose() < 2:
             sys.stdout = open(os.devnull, 'w')
-        
+
         self.solver.Run()
         self.solver.Postprocess()
         # write outputs
-        self.solver.Monitor(0) 
+        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:
+        if self.getVerbose() < 2:
             sys.stdout = sys.__stdout__
-        return 1
-    
+        if self.getVerbose() > 0:
+            print(f'SU2Interface::SU2 solver finished in {self.solver.GetOutputValue("INNER_ITER"):.0f} iterations. RMS_DENSITY: {self.solver.GetOutputValue("RMS_DENSITY"):.2f}')
+        return 0
+
     def writeCp(self, sfx=''):
         """ Write Cp distribution around the airfoil on file.
         """
@@ -109,36 +111,25 @@ class SU2Interface(SolversInterface):
                     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
 
     def getAoA(self):
         return self.solver.GetOutputValue('AOA') * np.pi/180
-    
+
     def getnDim(self):
         return self.solver.GetNumberDimensions()
-    
+
     def getnBodies(self):
         return 1
-    
+
     def getSpan(self, bodyname):
         return 1.0
-    
+
     def getBlowingTypes(self, bodyname):
         return['wing']
-    
+
     def getMinf(self):
         vel_x_idx = self.solver.GetPrimitiveIndices()["VELOCITY_X"]
         vel_y_idx = self.solver.GetPrimitiveIndices()["VELOCITY_Y"]
@@ -154,41 +145,35 @@ class SU2Interface(SolversInterface):
             sound_speed = primitives.Get(nodeIndex, sound_speed_idx)
             mach_farfield[inode] = np.sqrt(vel_x**2 + vel_y**2) / sound_speed
         return np.mean(mach_farfield)
-    
+
     def getCl(self):
         return self.solver.GetOutputValue('LIFT')
-    
+
     def getCd(self):
         return self.solver.GetOutputValue('DRAG')
 
     def getCm(self):
         return self.solver.Get_Mx()
-    
+
     def getVerbose(self):
         return self.verbose
-    
+
     def reset(self):
         if self.getVerbose() > 0:
             print('SU2Interface::Warning: Cannot reset solver.')
         pass
-    
+
     def getInviscidBC(self):
         """ Extract the inviscid paramaters at the edge of the boundary layer
         and use it as a boundary condition in the viscous solver.
-
-        Params
-        ------
-        - couplIter (int): Coupling iteration.
         """
-
-        # print('SU2Interface::Imposing inviscid boundary conditions.')
         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]
@@ -197,7 +182,6 @@ class SU2Interface(SolversInterface):
 
         for body in self.iBnd:
             for reg in body:
-                #print('SU2Interface::Imposing inviscid boundary conditions on', reg.name)
                 for inode, nodeIndexFloat in enumerate(reg.nodesCoord[:,-1]):
                     nodeIndex = int(nodeIndexFloat)
                     # Velocity
@@ -210,82 +194,82 @@ 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
-        # 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].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')
 
     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.')
-
-        # 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')
+        Extracts the blowing velocity from the viscous calculation and applies it
+        as a boundary condition in the inviscid solver.
+
+        Notes
+        -----
+        - In SU2, the "Grid velocity" is a vector defined at the nodes in the global
+        frame of reference.
+        - In BLASTER, the blowing velocity is computed as a **scalar normal component**
+        at the element level.
+        - To impose the correct boundary condition, the normal blowing velocity must
+        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)
+        normals = []
+        for ielm in range(len(self.iBnd[0][0].elemsCoord)):
+            nod1 = ielm
+            nod2 = ielm + 1
 
+            x1 = self.iBnd[0][0].nodesCoord[nod1,:self.getnDim()]
+            x2 = self.iBnd[0][0].nodesCoord[nod2,:self.getnDim()]
+
+            # Components of the tangent vector
+            t = np.empty(self.getnDim())
+            for idim in range(self.getnDim()):
+                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:
+                # 90° rotation
+                n[0] = t[1] / normt
+                n[1] = -t[0] / normt
+            elif self.getnDim() == 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]
+
+        # 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]):
+            # 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:
+                elm1 = 0
+                elm2 = -1
+            else:
+                elm1 = inode - 1
+                elm2 = inode
+
+            # Because the blowing velocity on the lower side points downwards
+            sign = 1
+            if inode > self.iBnd[0][0].nodesCoord.shape[0]//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])
+
+        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]
-                normal = self.solver.GetMarkerVertexNormals(self.wingMarkersToID[arfTag], ivertex)
-                blw = np.zeros(3)
-                for idim in range(self.getnDim()):
-                    blw += normal[idim] * blowingNodes[blasterIndex]
-                self.solver.SetBlowingVelocity(nodeIndex, blw[0], blw[1], blw[2])
+                self.solver.SetBlowingVelocity(nodeIndex, blw_nodes[blasterIndex, 0], blw_nodes[blasterIndex, 1], blw_nodes[blasterIndex, 2])
 
     def getWing(self):
         if self.getnDim() == 2:
@@ -299,33 +283,6 @@ 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))
-        # 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))
@@ -335,53 +292,53 @@ class SU2Interface(SolversInterface):
                 nodesOnSide[i, -1] = nodeIndex
             # Sort in ascending order of first column
             nodesOnSide = nodesOnSide[nodesOnSide[:,0].argsort()]
-            if arfTag == 'airfoil':
-                nodesOnSide = np.flip(nodesOnSide, axis=0)
-            elif arfTag == 'airfoil_':
-                nodesOnSide = np.delete(nodesOnSide, 0, axis=0)
+            if len(self.wingTags) > 1:
+                if arfTag == 'airfoil':
+                    nodesOnSide = np.flip(nodesOnSide, axis=0)
+                elif arfTag == 'airfoil_':
+                    nodesOnSide = np.delete(nodesOnSide, 0, axis=0)
+            else:
+                # If we cannot identify upper and lower sides, we sort using
+                # the angle between the node and the TE.
+                # This method can fail for highly cambered airfoils.
+
+                # Find the trailing edge
+                te = nodesOnSide[np.argmax(nodesOnSide[:, 0])]
+
+                # Compute angles relative to TE
+                angles = np.arctan2(nodesOnSide[:, 1] - te[1], nodesOnSide[:, 0] - te[0])
+
+                # Sort by angle in descending order
+                nodesOnSide = nodesOnSide[np.argsort(-angles)]
+
+                # Selig format
+                te_index = np.argmax(nodesOnSide[:, 0])  # Leading edge has min x
+                nodesOnSide[:te_index] = np.flip(nodesOnSide[:te_index], axis=0)
+                nodesOnSide[te_index:] = np.flip(nodesOnSide[te_index:], axis=0)
+
+                # Duplicate TE point
+                nodesOnSide = np.row_stack((nodesOnSide[-1], nodesOnSide))
             nodesCoord = np.row_stack((nodesCoord, nodesOnSide))
-        
+
         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
-        # 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]]):
             raise RuntimeError('su2Interface::Airfoil has an open trailing edge.')
-        else:
-            print('Airfoil has closed trailing edge.')
-        
+
         # Check that the airfoil has one point at the leading edge
         if len(np.where(nodesCoord[:,0] == np.min(nodesCoord[:,0]))[0]) > 1:
             print(len(np.where(nodesCoord[:,0] == np.min(nodesCoord[:,0]))[0]))
             raise RuntimeError('su2Interface::Airfoil has multiple points at the leading edge.')
-        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])
 
+        # 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

blast/tests/su2/blisu2.py

@@ -25,8 +25,9 @@ import blast.blUtils as vutils
 from fwk.wutils import parseargs
 
 def cfgSu2(verb):
+    import os
     return {
-    'filename' : '/home/c130/lab/blaster/blast/tests/su2/config.cfg',
+    'filename' : os.path.abspath(os.path.join(os.path.split(os.path.abspath(__file__))[0], "config.cfg")),
     'wingTags' : ['airfoil', 'airfoil_'],
     'Verb': verb,               # Verbosity level
     }
@@ -34,11 +35,11 @@ def cfgSu2(verb):
 def cfgBlast():
     return {
     'Re' : 1e7,                 # Freestream Reynolds number
-    'CFL0' : 1,                 # Inital CFL number of the calculation
-    'couplIter': 20,            # Maximum number of iterations
+    'CFL0' : 2,                 # Inital CFL number of the calculation
+    'couplIter': 10,            # Maximum number of iterations
     'sections': {'airfoil': [0.0]},
     'spans': {'airfoil': 1.0},
-    'couplTol' : 1e-6,          # Tolerance of the VII methodology
+    'couplTol' : 1e-4,          # 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
@@ -51,18 +52,19 @@ def main():
     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')
+    # Clear plt
+    plt.close('all')
+    plt.figure()
     plt.plot(cpCorrected[:,0], cpCorrected[:,1], lw=3)
     plt.plot(cpi[:,0], cpi[:,1], '--', lw=3)
     plt.gca().invert_yaxis()
-    plt.savefig('cp.png')
+    plt.show()
+    #plt.savefig('cp.png')
     quit()
 
     # Extract solution

blast/tests/su2/config.cfg

@@ -100,7 +100,7 @@ CFL_ADAPT_PARAM= ( 0.1, 5.0, 50.0, 1e10 )
 RK_ALPHA_COEFF= ( 0.66667, 0.66667, 1.000000 )
 %
 % Number of total iterations
-ITER= 1000
+ITER= 500
 %
 % ------------------------ LINEAR SOLVER DEFINITION ---------------------------%
 %
@@ -170,7 +170,7 @@ CONV_FIELD= RMS_DENSITY
 CONV_RESIDUAL_MINVAL= -8
 %
 % Start convergence criteria at iteration number
-CONV_STARTITER= 0
+CONV_STARTITER= 5
 %
 % Number of elements to apply the criteria
 CONV_CAUCHY_ELEMS= 100
@@ -180,11 +180,11 @@ CONV_CAUCHY_EPS= 1E-10
 
 % ------------------------- INPUT/OUTPUT INFORMATION --------------------------%
 %
-SCREEN_WRT_FREQ_INNER= 200
+SCREEN_WRT_FREQ_INNER= 50
 VOLUME_OUTPUT = SOLUTION, PRIMITIVE, RMS_DENSITY
 
 % Mesh input file
-MESH_FILENAME= /home/c130/lab/blaster/blast/models/su2/n0012_su2_2.su2
+MESH_FILENAME= /Users/pauldechamps/lab/softwares/blasterdev/blast/models/su2/n0012_su2_2.su2
 %
 % Mesh input file format (SU2, CGNS, NETCDF_ASCII)
 MESH_FORMAT= SU2