# Upwind scheme flux conservation not satisfied in 2D

I have read that the upwind scheme is flux conservative. Then by my understanding, in the absence of sinks/sources and with absorbing boundaries, the amount of the quantity leaving through the boundaries should be equal to the total initial amount of the quantity. However in practice this is not the result I'm getting for a 2D upwind scheme.

Consider a 2D advection equation with initial conditions: $$\frac{\partial u}{\partial t} + \nabla\cdot u = 0 \\ u = u(t,x,y) \\u(0, x, y) = u_0$$

I include an MFE of my code below, where ΔU calculates the update to each grid cell per time step, approximating $$-\Delta t (\frac{\partial u}{\partial x} + \frac{\partial u}{\partial y})$$. By my understanding, once the simulation has run for long enough, the sum of ΔU at the out-boundaries across time should be equal to $$u_0$$ but clearly this is not the case.

I think I'm missing something major and obvious here.

import matplotlib.pyplot as plt

def upwind(n, dx):
n = np.ones(n)
D = np.diag(n) - np.diag(n[1:], -1)
return D / dx

def kronecker_sum(Ax, Ay):
Ix = np.eye(*Ax.shape)
Iy = np.eye(*Ay.shape)
A = np.kron(Ax, Iy) + np.kron(Ix, Ay)
return A

nx = 10
ny = 10
dx = 1
dy = 1
dt = 0.1

U = np.zeros([ny, nx])
u0 = 1
U[0, 0] = 1

A = kronecker_sum(upwind(nx, dy), upwind(ny, dy))

# run explicit Euler for T timesteps
T = 20
ΔU = []
Umat = []
U1 = U.ravel(order='F')
for t in range(T):
update = -(dt * (A @ U1))
ΔU.append(update)
Umat.append(U1)
U1 = U1 + update

ΔU = np.array([du.reshape([ny, nx], order='F') for du in ΔU])
Umat = np.array([u.reshape([ny, nx], order='F') for u in Umat])

y_out = ΔU[:, -1, :]
x_out = ΔU[:, :, -1]

print(u0)
print(y_out.sum() + x_out.sum())

>>> 1
>>> 0.00010541525513958405

• Have you tried 1D to eliminate any errors due to the more complex system matrix? Furthermore, did you check if the values in U1 are constant over time? Commented Feb 12 at 18:20