I am making a toy code in Python to solve the advection equation $$u_t + cu_x = 0$$ with, for example, periodic boundary conditions.
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
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
u are the old and new data, respectively. All arrays are numpy arrays.
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.,
1:nx+1) and the implementation of my periodic boundary conditions function. It's not exactly intuitive that
u=u[-2]. The corresponding Fortran code would be
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.