Compute Divergence With Python

From this answer, the divergence of a numeric vector field can be computed as such: def divergence(f): num_dims = len(f) return np.ufunc.reduce(np.add, [np.gradient(f[i], a

Solution 1:

I realized what the issue was with the help of this answer. The default spacing assumed between two consecutive values in numpy.gradient() is 1. It needs to be changed if there is a different grid.

Hence the divergence function needs to be adapted as such:

Divergence function

defdivergence(f,sp):
    """ Computes divergence of vector field 
    f: array -> vector field components [Fx,Fy,Fz,...]
    sp: array -> spacing between points in respecitve directions [spx, spy,spz,...]
    """
    num_dims = len(f)
    return np.ufunc.reduce(np.add, [np.gradient(f[i], sp[i], axis=i) for i inrange(num_dims)])

Example

a = []
for N inrange(20,100):
    # Number of points (NxN)# = 20# Boundaries
    ymin = -2.; ymax = 2.
    xmin = -2.; xmax = 2.# Divergence functiondefdivergence(f,sp):
        num_dims = len(f)
        return np.ufunc.reduce(np.add, [np.gradient(f[i], sp[i], axis=i) for i inrange(num_dims)])

    
    # Create Meshgrid
    x = np.linspace(xmin,xmax, N)
    y = np.linspace(ymin,ymax, N)
    xx, yy = np.meshgrid(x, y)
    
    
    # Define 2D Vector Field
    Fx  = np.cos(xx + 2*yy)
    Fy  = np.sin(xx - 2*yy)
    
    F = np.array([Fx, Fy])
    # Compute Divergence
    sp_x = np.diff(x)[0]
    sp_y = np.diff(y)[0]
    sp = [sp_x, sp_y]
    g = divergence(F, sp)
    
    print("Max: ", np.max(g.flatten()))
    a.append(np.max(g.flatten()))
plt.plot(a)

We can see that with increasing resolution, the maximum of the divergence really tends to 3.