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.

enter image description here

The initial condition has non zero pixels in each corner where it is basically a plateau disc function. initial condition

The velocity field points always to the center of the domain, so that the points have to meet in the center. enter image description here

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::,:] - 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)

def update(i):

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


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

  • $\begingroup$ Just a couple of comments: (1) Your solution will lead to a $\delta$ function. (2) Numerically a generic solution to the advection equations will lead to loss of mass unless they are conservative so it is not entirely surprising the behavior you are seeing. You mention that you are using upwind and downwind schemes but I am not entirely sure that you are in fact using those schemes honestly. $\endgroup$ Aug 29 at 0:10
  • $\begingroup$ Also if you want to test to see if your solution is at least correct in principle you could try refining the grid as it should still converge despite lack of conservation although you won't be able to converge to the true solution really. $\endgroup$ Aug 29 at 0:15
  • 1
    $\begingroup$ That's a lot of numerical diffusion. You'll have to use higher order spacial and temporal schemes. Look for second order upwind or WENO spacial schemes and Runge kutta temporal schemes. $\endgroup$ Aug 29 at 13:23

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.