我试图用速度场来实现二维平流问题，这在空间中不是恒定的。我的问题是，我的转移密度的“质量”被“侵蚀”或消失了，即规范随着时间的推移变为零。至少我想以某种方式最小化它的衰变。$\|u_t\|$

初始条件在每个角落都有非零像素，它基本上是一个高原圆盘函数。

速度场总是指向域的中心，因此这些点必须在中心相遇。

这是我的代码。

    # -*- coding: utf-8 -*-

import numpy as np
from matplotlib import pyplot as plt
from matplotlib import animation 


n = 32
X = np.linspace(-2,2,n)
Y = np.linspace(-2,2,n)

XX,YY = np.meshgrid(X,Y)


#------------- velocity field --------
# a disc in each corner
u0 = np.zeros((n,n))
u0 += (XX-1.25)**2 + (YY-1.25)**2 < 0.25
u0 += (XX+1.25)**2 + (YY+1.25)**2 < 0.25
u0 += (XX-1.25)**2 + (YY+1.25)**2 < 0.25
u0 += (XX+1.25)**2 + (YY-1.25)**2 < 0.25


# flow field points always to origin
flow = np.zeros((2,n,n))
flow[0,:,:] = YY
flow[1,:,:] = XX
flow *= -1 # 

# normalizing flow field - is it necessary?
norm = np.linalg.norm(flow, axis = 0, ord=2)    
flow[0,:,:] =  np.divide(flow[0,:,:], norm, out=np.zeros_like(norm), where=norm>10**-7)
flow[1,:,:] =  np.divide(flow[1,:,:], norm, out=np.zeros_like(norm), where=norm>10**-7)


dx = (XX.max() - XX.min())/n
dt = dx/5

u = np.copy(u0)
U = [u0.copy()]

frames = 100

for i in range(frames):
    
    #upwind scheme
    ux_minus = np.vstack((np.zeros((1,n)), u[1:-1,:] - u[0:-2,:],np.zeros((1,n))))  
    uy_minus = np.hstack((np.zeros((n,1)), u[:,1:-1] - u[:,0:-2],np.zeros((n,1))))     
    
    # downwind scheme
    ux_plus = np.vstack((np.zeros((1,n)), u[2::,:] - u[1:-1,:], np.zeros((1,n))))
    uy_plus = np.hstack((np.zeros((n,1)), u[:,2::] - u[:,1:-1], np.zeros((n,1))))      
    
    pulse_x = (flow[0,:,:] < 0) * flow[0,:,:]*ux_plus + (flow[0,:,:]>0)*flow[0,:,:]*ux_minus    
    pulse_y = (flow[1,:,:] < 0) * flow[1,:,:]*uy_plus + (flow[1,:,:]>0)*flow[1,:,:]*uy_minus
    
    u -= (dt/dx)*(pulse_x + pulse_y)
        
    U.append(u.copy())
    


def update(i):
    matrice.set_array(U[i])


fig, ax = plt.subplots()
plt.show()
matrice = ax.matshow(U[0])

plt.colorbar(matrice)

ani = animation.FuncAnimation(fig, update, frames=frames, interval=1)