Getting gradient vector of surface point sensitivity in adjoint calculation.

[atharvaaalok]

Hello,
I am trying to use SU2 for shape optimization with a geometry representation that I have developed.
In particular, I have a mesh around my shape and I want to get the gradient vector of sensitivity for each point on the surface of my shape using the continuous/discrete adjoint.
I do not want to perform the optimization with any SU2 tools, that I wish to do in my own way and wish to use SU2 only to compute the flow field and the adjoint sensitivity vector for Lift and Drag.

I ran the quick start tutorial:
In the surface_adjoint.csv file that is created, the last column titled "Surface_Sensitivity" is a set of scalar numbers. (I am guessing this is the magnitude of the sensitivity in the normal direction)

How do I obtain the actual sensitivity vector (dx, dy) for every surface point?

Attaching my output csv for the surface_adjoint.
Note, I am running the script as detailed in the quickstart tutorial.

surface_adjoint.csv

I only see numbers in the Surface_Sensitivity column and 0s in the Sensitivity_x and Sensitivity_y columns.
Why is that?

[pcarruscag]

You need to run SU2_DOT_AD to project from volume to surface

[atharvaaalok]

Thank you!
Where can I find the documentation for what exactly is the 'surface_sensitivity' in the adjoint.csv file and further details on using the SU2_DOT_AD for this specific case?

[pcarruscag]

there is no exact documentation.
you can look in the TestCases folder for tests that call SU2_DOT, and here is an example of an optimization put together outside of SU2 https://github.com/pcarruscag/FADO/blob/master/examples/example4_SU2/example.py

[bigfooted]

The optimization methods that we have basically run SU2_CFD, SU2_CFD_AD, SU2_DOT_AD and then they use an optimizer like SLSQP to do something clever with the (surface) sensitivities. Then we run SU2_DEF to deform the mesh with these sensitivities and do it all over again.

[atharvaaalok]

@bigfooted is there any code that I could look at and use to get my dJ/dx and dJ/dy values? (where J is my QoI)
I am getting lost in the sea of methods. I feel like I should be very close to getting these after I am done with the quickstart.
Can I use SU2_DOT here? If yes, how?

"Something like: After you have run the adjoint simulation from the quickstart run this command to get dJ/dx."

[atharvaaalok]

Also @pcarruscag what is 'surface_sensitivity' exactly? It is a set of scalars.
Is this the magnitude of the sensitivity in the normal direction?
If yes, why is only the normal sensitivity reported?

[atharvaaalok]

SU2 reports normal sensitivity as any tangential component would just mean a reparameterization and does not change the shape.

The simplest way therefore is to compute normal vectors for the points on the surface and then multiply the surface sensitivity with these. A simple code snippet is attached below for reference:

def compute_area(X):
    # Shoelace formula - works for any simple polygon (piecewise-linear simple closed curve)
    area = np.sum(X[:-1, 0] * X[1:, 1] - X[1:, 0] * X[:-1, 1])

    return area.item()


def compute_normals(X):
    # Panel vectors - these are tangents between consecutive points on the shape
    dX = X[1:] - X[:-1]
    dx = dX[:, 0]
    dy = dX[:, 1]

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

    # Get unit tangent vectors for the panels
    tx = dx / lengths
    ty = dy / lengths

    # Store into a matrix and
    T = np.stack([tx, ty], axis = 1)
    T = np.vstack([T[-1], T, T[0]])
    
    # Use average of the two surrounding tangents to approximate tangent at coordinate point
    T = (T[1:] + T[:-1]) / 2
    
    # Outward normals - rotate tangents 90 degree clockwise (dy, -dx) / length
    # This assumes that the points go around the curve in an counterclockwise fashion
    N = np.stack([T[:, 1], -T[:, 0]], axis = 1)
    N = N / np.linalg.norm(N, axis = 1, keepdims = True)

    # Signed area test - flip normals if points go around the curve in a clockwise fashion
    area = compute_area(X)
    if area < 0:
        N = -N

    return N