Attempt on 2d Advection with FDM - With Code
I tried to implement the 2d advection problem with a velocity field, that is not constant in space. My problem is, that the "mass" of my shifted density gets "eroded" or just disappears, i.e. the norm $\|u_t\|$ goes to zero as time progresses. At least I want to minimize its decay somehow.
The initial condition has non zero pixels in each corner where it is basically a plateau disc function.
The velocity field points always to the center of the domain, so that the points have to meet in the center.
Here is my code.
# -*- 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<IP_ADDRESS>,:] - u[1:-1,:], np.zeros((1,n))))
uy_plus = np.hstack((np.zeros((n,1)), u[:,2<IP_ADDRESS>] - 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) 935