Possible Bug: Computation of CL CD from history.csv and flow.vtu differ by a factor of 2.

[atharvaaalok]

Describe the bug
Given an airfoil shape using coordinates X (which is a (N, 2) matrix of x, y pairs) I first run it through SU2 to generate the different output files: history.csv, flow.vtu, surface_flow.csv.

Step 1: I first parse the history.csv file and read the last value of CL (that is the one obtained after convergence).

Step 2: Then I read the surface coordinate points from surface_flow.csv and use the other data rho_u, rho_v etc to compute the pressure coefficient cp. Another way to get cp is to simply read the pressure coefficient values from flow.vtu at the surface points.

Step 3: Compute CL and CD using CP distribution on the surface. I use the function described below.

Step 4: Compare the CL compute as in step 3 with the one read from history.csv in step 1.

Observation: The CL computed from Step 3 is for all shapes 2 times the CL reported in history.csv.

def compute_lift_drag(X, CP):
    # Segment vectors
    dX = X[1:] - X[:-1]
    dx = dX[:, 0]
    dy = dX[:, 1]

    # Segment lengths
    lengths = np.sqrt(dx**2 + dy**2)

    # Outward normals: rotate 90 deg clockwise (dy, -dx) / length
    nx = dy / lengths
    ny = -dx / lengths

    nx_times_length = dy
    ny_times_length = -dx

    # Average pressure per segment: shape (N-1,)
    CP_avg = 0.5 * (CP[1:] + CP[:-1])

    # Pressure force per segment: -CP_avg * normal * length
    fx = -CP_avg * nx_times_length
    fy = -CP_avg * ny_times_length

    # Total forces
    CD = fx.sum()
    CL = fy.sum()

    return CL, CD

Desktop (please complete the following information):

• OS: MacOS
• C++ compiler and version: Apple clang version 17.0.0 (clang-1700.0.13.5)
• MPI implementation and version: HYDRA build details: Version: 4.3.0
• SU2 Version: v8.2.0 Harrier

[pcarruscag]

There is no bug in cl and cd, you can look at the the v&v page.

[atharvaaalok]

@pcarruscag Can you point me to the SU2 code snippet (python/C++) that computes CL and CD from the generated simulation data? I can compare that against my own.

[pcarruscag]

SU2/SU2_CFD/include/solvers/CFVMFlowSolverBase.inl

Line 1783 in deece11

void CFVMFlowSolverBase<V, FlowRegime>::Pressure_Forces(const CGeometry* geometry, const CConfig* config) {

SU2/SU2_CFD/include/solvers/CFVMFlowSolverBase.inl

Line 2384 in deece11

void CFVMFlowSolverBase<V, FlowRegime>::Friction_Forces(const CGeometry* geometry, const CConfig* config) {

[atharvaaalok]

I am not able to parse through this massive code base to figure out where exactly that 2 is missing.
Started a discussion instead: #2556

[atharvaaalok]

The bug was in my meshing code it turns out.
Gmsh was ordering points on the surface such that they were traversing the shape twice and therefore the factor of 2 error in the CL and CD computations.

For reference attaching code to compute CL and CD using Cp and skin friction values on points on the surface.

def compute_lift_drag(X, Cp, Cfx = None, Cfy = None):
    if Cfx is None:
        Cfx = Cp * 0
    if Cfy is None:
        Cfy = Cp * 0
    
    # Segment vectors
    dX = X[1:] - X[:-1]
    dx = dX[:, 0]
    dy = dX[:, 1]

    # Segment lengths
    lengths = np.sqrt(dx**2 + dy**2 + 1e-12)

    # Outward normals: rotate 90 deg clockwise (dy, -dx) / length
    nx = dy / lengths
    ny = -dx / lengths

    nx_times_length = dy
    ny_times_length = -dx

    # Average pressure per segment: shape (N-1,)
    Cp_avg = 0.5 * (Cp[1:] + Cp[:-1])
    Cfx_avg = 0.5 * (Cfx[1:] + Cfx[:-1])
    Cfy_avg = 0.5 * (Cfy[1:] + Cfy[:-1])

    # Pressure force per segment: -CP_avg * normal * length
    fx = -Cp_avg * nx_times_length + Cfx_avg * lengths
    fy = -Cp_avg * ny_times_length + Cfy_avg * lengths

    # Total forces
    CD = fx.sum()
    CL = fy.sum()

    return CL, CD