(feat) First test for SU2 trailing-edge smoothing

blast/interfaces/su2/blSU2Interface.py

@@ -168,14 +168,18 @@ class SU2Interface(SolversInterface):
         """ Extract the inviscid paramaters at the edge of the boundary layer
         and use it as a boundary condition in the viscous solver.
         """
-        if self.getnDim() == 2:
+        region = self.iBnd[0][0]
+        wingRows = region.getNodeRows()
+        nNodes, ndim = region.getnNodes(), self.getnDim()
+
+        if ndim == 2:
             velComponents = ['VELOCITY_X', 'VELOCITY_Y']
-        elif self.getnDim() == 3:
+        elif ndim == 3:
             velComponents = ['VELOCITY_X', 'VELOCITY_Y', 'VELOCITY_Z']
         else:
             raise RuntimeError('su2Interface::Incorrect number of dimensions.')
 
-        velIdx = np.zeros(self.getnDim(), dtype=int)
+        velIdx = np.zeros(ndim, dtype=int)
         for i, velComp in enumerate(velComponents):
             velIdx[i] = self.solver.GetPrimitiveIndices()[velComp]
         soundSpeedIdx = self.solver.GetPrimitiveIndices()["SOUND_SPEED"]
@@ -183,18 +187,117 @@ class SU2Interface(SolversInterface):
 
         for body in self.iBnd:
             for reg in body:
-                for inode, nodeIndexFloat in enumerate(reg.nodesCoord[:,-1]):
+                V = np.empty((nNodes, ndim))
+                M = np.empty(nNodes)
+                Rho = np.empty(nNodes)
+                for inode, nodeIndexFloat in enumerate(wingRows):
                     nodeIndex = int(nodeIndexFloat)
                     # Velocity
-                    for idim in range(self.getnDim()):
-                        reg.V[inode,idim] = self.solver.Primitives().Get(int(nodeIndex), int(velIdx[idim]))
+                    for idim in range(ndim):
+                        V[inode,idim] = self.solver.Primitives().Get(int(nodeIndex), int(velIdx[idim]))
                     # Density
-                    reg.Rho[inode] = self.solver.Primitives().Get(nodeIndex, densityIdx)
+                    Rho[inode] = self.solver.Primitives().Get(nodeIndex, densityIdx)
                     # Mach number
                     vel_norm = 0.0
-                    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)
+                    for idim in range(ndim):
+                        vel_norm += V[inode,idim]**2
+                    M[inode] = np.sqrt(vel_norm) / self.solver.Primitives().Get(nodeIndex, soundSpeedIdx)
+                reg.updateVariables(M, V, Rho)
+
+        from matplotlib import pyplot as plt
+        plt.plot(reg._nodesCoord[:,0], reg._M, 'x-')
+
+        # Reflect the last two nodes
+        nodes_ = region.getNodes().copy()
+        M_ = M.copy()
+        new_nodes, new_M = self.orthogonal_reflection(nodes_, M_)
+        #plt.plot(new_nodes[:,0], new_M, 'o-')
+        #plt.show()
+        smoothed_nodes, smoothed_M = self.explicit_smoother(new_nodes, new_M, iterations=10)
+        plt.plot(smoothed_nodes[:,0], smoothed_M, 'o-')
+
+        plt.show()
+        quit()
+        plt.show()
+    
+    def orthogonal_reflection(self, nodes, M, axis_x=1):
+        """
+        Reflect the first two and last two points of the Mach number data (simulating wake) with respect
+        to the given axis (default trailing edge at x=1), and append them to the data.
+        
+        Args:
+        - nodes: 2D array where each row contains [x, y, z] coordinates.
+        - M: 1D array containing the Mach numbers.
+        - axis_x: The x-coordinate about which to reflect (default is 1, the trailing edge).
+        
+        Returns:
+        - new_nodes: Reflected nodes including the additional wake simulation points.
+        - new_M: Reflected Mach numbers.
+        """
+        # Reflect the first two nodes and Mach numbers
+        first_two_nodes = nodes[1:5]   # First two [x, y, z]
+        first_two_M = M[1:5]           # First two Mach numbers
+        
+        # Orthogonal reflection: Reflect the x-values with respect to x=axis_x
+        reflected_first_two_nodes = np.copy(first_two_nodes)
+        reflected_first_two_nodes[:, 0] = 2*axis_x - first_two_nodes[:, 0] # Reflect x-values
+        
+        # Reflect the Mach numbers
+        reflected_first_two_M = first_two_M[::-1]  # Reverse the order of the first two Mach numbers
+        
+        # Reflect the last two nodes and Mach numbers
+        last_two_nodes = nodes[-5:-1]   # Last two [x, y, z]
+        last_two_M = M[-5:-1]           # Last two Mach numbers
+        
+        # Orthogonal reflection: Reflect the x-values with respect to x=axis_x
+        reflected_last_two_nodes = np.copy(last_two_nodes)
+        reflected_last_two_nodes[:, 0] = 2*axis_x - last_two_nodes[:, 0]  # Reflect x-values
+        
+        # Reflect the Mach numbers
+        reflected_last_two_M = last_two_M[::-1]  # Reverse the order of the last two Mach numbers
+        
+        # Append the reflected points to the original nodes and Mach numbers
+        new_nodes = np.vstack((reflected_first_two_nodes, nodes, reflected_last_two_nodes))
+        new_M = np.concatenate((reflected_first_two_M, M, reflected_last_two_M))
+        
+        return new_nodes, new_M
+
+    def explicit_smoother(self, nodes, M, mask_range=(0.8, 1.0), iterations=5):
+        """
+        Smooth the Mach number data only for nodes in a given range (between mask_range[0] and mask_range[1]).
+        The smoothing uses explicit finite difference method with Neumann boundary conditions.
+        
+        Args:
+        - nodes: 2D array where each row contains [x, y, z] coordinates.
+        - M: 1D array containing the Mach numbers.
+        - mask_range: Tuple (start_pct, end_pct) defining the range of x-values to apply smoothing.
+        - iterations: Number of smoothing iterations to perform.
+        
+        Returns:
+        - smoothed_nodes: Original nodes (unchanged).
+        - smoothed_M: Smoothed Mach number data.
+        """
+        # Create a mask for the region between 0.8 and 1.0 of the chord
+        mask = (nodes[:, 0] > mask_range[0]) & (nodes[:, 0] < mask_range[1])
+        
+        smoothed_M = M.copy()  # Copy the Mach numbers for smoothing
+        
+        for _ in range(iterations):
+            new_M = smoothed_M.copy()
+            
+            # Apply smoothing only to the values in the specified range (using the mask)
+            for i in range(1, len(M) - 1):
+                if mask[i]:  # Only smooth where the mask is True
+                    new_M[i] = 0.5 * (smoothed_M[i - 1] + smoothed_M[i + 1])
+            
+            te_index = np.argmax(nodes[:, 0])  # Index of the trailing edge
+            # Neumann boundary condition: no change at the first and last points
+            new_M[0] = smoothed_M[1]  # Set first point equal to the second point
+            new_M[-1] = smoothed_M[-2]  # Set last point equal to the second-to-last point
+            
+            smoothed_M = new_M
+        
+        return nodes, smoothed_M
 
     def setBlowingVelocity(self):
         """
@@ -211,41 +314,47 @@ class SU2Interface(SolversInterface):
         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)
+        region = self.iBnd[0][0]
+        nNodes, nElms, ndim = region.getnNodes(), region.getnElms(), self.getnDim()
+        rowsWing = region.getNodeRows()
+        blowingWing = region.getBlowingVelocity()
+        nodesWing = region.getNodes()
+
+        blw_elems_normal = np.zeros((nElms, ndim), dtype=float)
         normals = []
-        for ielm in range(len(self.iBnd[0][0].elemsCoord)):
+        for ielm in range(len(self.iBnd[0][0].getElms('all'))):
             nod1 = ielm
             nod2 = ielm + 1
 
-            x1 = self.iBnd[0][0].nodesCoord[nod1,:self.getnDim()]
-            x2 = self.iBnd[0][0].nodesCoord[nod2,:self.getnDim()]
+            x1 = nodesWing[nod1,:ndim]
+            x2 = nodesWing[nod2,:ndim]
 
             # Components of the tangent vector
-            t = np.empty(self.getnDim())
-            for idim in range(self.getnDim()):
+            t = np.empty(ndim)
+            for idim in range(ndim):
                 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:
+            n = np.empty(ndim)
+            if ndim == 2:
                 # 90° rotation
                 n[0] = t[1] / normt
                 n[1] = -t[0] / normt
-            elif self.getnDim() == 3:
+            elif ndim == 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]
+            for idim in range(ndim):
+                blw_elems_normal[ielm, idim] = blowingWing[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]):
+        blw_nodes = np.zeros((nNodes, 3), dtype=float)
+        for inode in range(nNodes):
             # 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:
+            if inode == 0 or inode == nNodes-1:
                 elm1 = 0
                 elm2 = -1
             else:
@@ -254,23 +363,26 @@ class SU2Interface(SolversInterface):
 
             # Because the blowing velocity on the lower side points downwards
             sign = 1
-            if inode > self.iBnd[0][0].nodesCoord.shape[0]//2:
+            if inode > nNodes//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])
+            for idim in range(ndim):
+                blw_nodes[inode, idim] = sign * 0.5*(blw_elems_normal[elm1, idim] + blw_elems_normal[elm2, idim])
 
-        blw_nodes[blw_nodes < -0.01] = -0.01
-        blw_nodes[blw_nodes > 0.01] = 0.01
+        # 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]
-                self.solver.SetBlowingVelocity(nodeIndex, blw_nodes[blasterIndex, 0], blw_nodes[blasterIndex, 1], blw_nodes[blasterIndex, 2])
+                blasterIndex = np.where(rowsWing == nodeIndex)[0][0]
+                self.solver.SetBlowingVelocity(nodeIndex,
+                                               blw_nodes[blasterIndex, 0],
+                                               blw_nodes[blasterIndex, 1],
+                                               blw_nodes[blasterIndex, 2])
 
     def getWing(self):
         if self.getnDim() == 2:
@@ -340,11 +452,14 @@ class SU2Interface(SolversInterface):
             raise RuntimeError('su2Interface::Duplicated elements in airfoil mesh.')
 
         # 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
-        self.iBnd[ibody][0].setConnectList(np.array(nodesCoord[:,-1],dtype=int),
+        self.iBnd[ibody][0].setMesh(nodesCoord[:,:3], elemsCoord[:,:3])
+        #self.iBnd[ibody][0].initStructures(nodesCoord.shape[0], elemsCoord.shape[0])
+        #self.iBnd[ibody][0].nodesCoord = nodesCoord
+        #self.iBnd[ibody][0].elemsCoord = elemsCoord
+        self.iBnd[ibody][0].setConnectLists(np.array(nodesCoord[:,-1],dtype=int),
                                            np.array(elemsCoord[:,-1],dtype=int))
+        self.iBnd[ibody][0].setRows(np.array(nodesCoord[:,-1],dtype=int),
+                                    np.array(elemsCoord[:,-1],dtype=int))
 
     def _modify_meshpath(self, config, mesh):
         # Open config file and look for line MESH_FILENAME

blast/tests/su2/t_n0012_su2_2D.py

@@ -50,13 +50,12 @@ def cfgBlast():
     }
 
 def main():
-    return 0
     args = parseargs()
     icfg = cfgSu2(args.verb)
     vcfg = cfgBlast()
     coupler, isol, vsol = vutils.initBlast(icfg, vcfg, iSolver='SU2', task='analysis')
-    #coupler.run()
-    isol.run()
+    coupler.run()
+    #isol.run()
 
     print(ccolors.ANSI_BLUE + 'PyTesting...' + ccolors.ANSI_RESET)
     tests = CTests()