0
$\begingroup$

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.

$\endgroup$
3
$\begingroup$

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))
$\endgroup$
  • $\begingroup$ Perfect! Thanks. $\endgroup$ – EternusVia Sep 5 at 15:11

Your Answer

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

Not the answer you're looking for? Browse other questions tagged or ask your own question.