# Implementing structured grid boundary conditions using NumPy arrays?

I am making a toy code in Python to solve the advection equation $$u_t + cu_x = 0$$ with, for example, periodic boundary conditions.

Background information

The numerical grid is specified like this:

# problem parameters
lx = 2.0*np.pi    # length of the domain
c = 1.0           # advecton speed

# grid generation
nx = 50            # number of grid points
dx = lx/float(nx)  # grid spacing
x = np.linspace(-dx/2.,lx+dx/2.,nx+2) # vector of grid points


where the -dx/2. and lx+dx/2. points correspond to the ghost points.

I currently implement the periodic boundary conditions like this:

def periodic_bc(u):
u[ 0]=u[-2]
u[-1]=u[ 1]
return u


and the time integration is written as follows:

for n in range(0,nt):
u[1:nx+1] = u0[1:nx+1] + (c*dt/dx)*(u0[0:nx]-u0[1:nx+1])
u = periodic_bc(u)
u0 = u


where u0 and u are the old and new data, respectively. All arrays are numpy arrays.

Question

Does anyone have suggestions for implementing boundary conditions in a more elegant fashion?

Even when I was programming this extremely simple code, I found myself confused about the bounds of my array slices (e.g., 0:nx, 1:nx+1) and the implementation of my periodic boundary conditions function. It's not exactly intuitive that u[0]=u[-2]. The corresponding Fortran code would be u[-1]=u[nx-1], where u[-1] can be allocated specifically for the ghost points. To me this is much more intuitive and matches the discrete boundary conditions.

Or maybe there is a different way besides ghost points that allows a much more intuitive code. Any advice would be welcome.

The Numpy function $$roll$$ performs periodic shift of an array. Using it, the explicit time step for your PDE in a periodic domain can be simply implemented like this:
u = u - (1/2)*(c*dt/dx)*(np.roll(u,-1) - np.roll(u,1))