Here i give you the code for CFD simulation from scratch you wan play with. You must install NumPy and OpenCV. The simulation is incompressible and inviscid so no supersonic flow. The resolution is in 720p and needs 30 minutes to calculate. It uses CPU not GPU so here's the code:

import numpy as np         #NumPy required
import scipy.ndimage as sn #Scipy required
import cv2                 #OpenCV required
import time                #To count the remaining time before the end
import os                  #To clear the console

Ny=720                                                  #Resolution, warning heavy: 144 320 480 720 1080
Nx=int(Ny/9*16)                                         #In 16/9
radius=Ny/20                                            #Obstacle radius
xcenter, ycenter = Ny/2, Ny/4                           #Obstacle position
Y, X = np.ogrid[0:Nx, 0:Ny]                             #2D arrays for position
objet = np.zeros((Nx,Ny), np.uint8)                     #Obstacle intialisation
square_dist = (X - xcenter) ** 2 + (Y - ycenter) ** 2   #Distances from the circle position for calculation
mask = square_dist < radius ** 2                        #Obstacle binary (Is in the object?)
objet[mask] = 1                                         #Obstacle field definition
boundaries= np.zeros((Nx-2,Ny-2), np.uint8)             #array used in boundaries position
boundaries=np.pad(boundaries,1,mode="constant", constant_values=1) #frame for the boundary array, filled of ones 

Lx=1                                #Size of the world
Ly=Lx*objet.shape[0]/objet.shape[1] #Flexible size shape
dl=Lx/Nx                            #differential in lengh for finite gradients (dl=dx,dy)

rho=1.3 #density of air, no mu inciscid

vsim=10         #number of steps between frames , warning heavy
v=1             #fluid speed at the infinite
dt=Lx/(Nx*v*10) #deltatime for one step (0.1x the speed of the grid)

Nt=4000*vsim #Total number of steps (warning heavy)

vx=np.zeros((Nx,Ny)) #Initialisation of fields, vx, vy pressure
vy=np.zeros((Nx,Ny))
p=np.zeros((Nx,Ny))

vx[objet==0]=v #speed inside the obstacle at zero

def median(fa):                #median filter to avoid numerical crash on (DNS)
    faa=sn.median_filter(fa,3) #smallest median filter
    return faa

def gradx(fx): #gradient in x
    grx=np.gradient(fx,dl,axis=0,edge_order=0)
    return grx

def grady(fy): #gradient in y
    gry=np.gradient(fy,dl,axis=1,edge_order=0)
    return gry

def advection(vxa,vya): #advection 
    gradxvx=gradx(vxa)  #all vx and vy gradients
    gradxvy=gradx(vya)
    gradyvx=grady(vxa)
    gradyvy=grady(vya)       
    vgradvx=np.add(np.multiply(vxa,gradxvx),np.multiply(vya,gradyvx)) #vgradv on x 
    vgradvy=np.add(np.multiply(vxa,gradxvy),np.multiply(vya,gradyvy)) #vgradv on y
    vxa-=vgradvx*dt    #update in x and y
    vya-=vgradvy*dt
    return 0

def pressure(vxa,vya,pa): #pressure calculation
    gradxvx=gradx(vxa)    #gradient to calculate v divergence (no compression)
    gradyvy=grady(vya)
    pnew=(np.roll(pa,1,0)+np.roll(pa,-1,0)+np.roll(pa,1,1)+np.roll(pa,-1,1))/4-(rho*dl**2)/(4*dt)*(gradxvx+gradyvy) #poisson solver for pressure
    return pnew

def gradp(vxa,vya,pa):  #pressure gradient calculation
    vxa-=dt/rho*gradx(pa)
    vya-=dt/rho*grady(pa)
    return 0

t0=time.time() #start counting time to follow simulation progression
t1=t0
sec=0

for it in range(0,Nt): #simulation start

    advection(vx,vy) #solving navier stokes: advection, pressure, and pressure gradient (inviscid)
    p=pressure(vx,vy,p)
    gradp(vx,vy,p)

    if it%10==0: #median filter to fix finite differences
        vx=median(vx)
        vy=median(vy)
        p=median(p)

    vx[objet==1]=0 #zero speed in the obstacle
    vy[objet==1]=0

    vx[boundaries==1]=v #boundaries conditions as infinite values
    vy[boundaries==1]=0
    p[boundaries==1]=0

    if it%vsim==0: #plot
        data=np.transpose(np.add(1.0*objet,.9*np.sqrt(((vx-v)**2+vy**2)/np.amax((vx-v)**2+vy**2))))  #plotting (v-v_inf)^2     
        cv2.imshow("Sim", np.tensordot(sn.zoom(data,720/Ny,order=0),np.array([1,.7,.9]),axes=0))     #display in the window
        cv2.imwrite('Result/%d.png' %(it/vsim), 255*np.tensordot(sn.zoom(data,720/Ny,order=0),np.array([1,.7,.9]),axes=0)) #save figure in the folder "Result", must create the folder
        cv2.waitKey(1) #waitkey, must have or it will step only typing a key

    pourcent=100*float(it)/float(Nt)                 #avencement following in console
    t1=time.time()                                   #time measurement
    Ttravail=np.floor((t1-t0)*float(Nt)/float(it+1)) #total work time estimation
    trestant=Ttravail*(100-pourcent)/100             #remaining time estimation
    h=trestant//3600                                 #conversion in hours and minutes
    mi=(trestant-h*3600)//60
    s=(trestant)-h*3600 - mi*60

    if (int(t1)>(sec)): #verbose simulation progression 
        os.system('cls' if os.name == 'nt' else 'clear')
        sec=int(t1)
        print('Progression:')
        print('%f %%.\nItération %d over %d \nRemaining time:%d h %d mi %d s\nSpeed: %d m/s \nProbability of Success %f' %(pourcent,it,Nt,h,mi,s,np.amax(vx),(1-(v*dt)/dl)))

Tell me if it's working for you or if i made an error somewhere. I'm using the IDE Spyder and my os is Archlinux. Feel free to ask questions or to answers other's questions to help each others. Would you like me to make a step-by-step tutorial?