How can I enforce PBC when simulating a 2D Reaction-Diffusion system?-CodePudding

i'm having some problems attempting to implement periodic boundary conditions (PBC) on a reaction diffusion system simulated in Python using 2D numpy arrays. I'll try to explain using pseudocode and attach code as to how i'm currently handling the boundaries.

import numpy as np

N = 100

# I define a 2D array for each of my species in the reaction-diffusion system
a = np.array((N, N), dtype=np.float64)
b = np.array((N, N), dtype=np.float64)
.
.
.
n = np.array((N, N), dtype=np.float64)

# And also copies to update at each time step
a_copy = np.array((N, N), dtype=np.float64)
b_copy = np.array((N, N), dtype=np.float64)
.
.
.
n_copy = np.array((N, N), dtype=np.float64)

# I calculate the laplacian using the following function
@jit(nopython=True, fastmath=True)
def laplacian_numba(field, dh2_inv, out):
    """
    Compute the laplacian of an array using a 5 point stencil.
    """
    for i in range(1, N - 1):
        for j in range(1, N - 1):
            out[i, j] = (
                field[i   1, j]
                  field[i - 1, j]
                  field[i, j   1]
                  field[i, j - 1]
                - 4 * field[i, j]
            ) * dh2_inv
    return out

# I have a main loop to update the functions using an explicit method
@jit(nopython=True, fastmath=True)
def update(a, b, ..., n):

    # Compute the laplacian of only the diffusing variables
    laplacian_numba(a, dh2_inv, out=lap_a)
    laplacian_numba(b, dh2_inv, out=lap_b)

    # Update the copy arrays
    a_copy = a[1:-1, 1:-1]   dt * (ODE stuff   lap_a * a[1:-1, 1:-1])
    ...

    # Finally I enforce the boundary conditions on the system
    # If i'm not mistaken, these are reflecting boundary conditions, not periodic
    # And this is where i'm lost as to how to implement the periodicity
    a_copy[0,:] = a_copy[1,:]
    a_copy[-1,:] = a_copy[-2,:]
    a_copy[:,0] = a_copy[:,1]
    a_copy[:,-1] = a_copy[:,-2]

    # Update the previous and next timestep arrays
    a, a_copy = a_copy, a
    b, b_copy = b_copy, b  

    return a, b, ..., n

The above is a very crude pseudocode of how I implemented the system, and how I update at each timestep and enforce the boundary conditions, which if I'm not mistaken are just reflecting the edges back onto the main grid. Here is my main question, what would I need to change in order to make the edges periodic and not reflecting? As I come from a biochemistry background, PDEs are not my forte, but I'm making an effort since this is a key objective in my thesis and would appreciate any help or guidance.

Thanks in advance to anyone who takes the time to read this! And I apologize for any formatting mistakes I could've made.

CodePudding user response：

With periodic boundary conditions, there isn't really a boundary; your coordinate space just wraps around modulo N. So there is no need to add special first/last rows/columns at all, and no need to explicitly enforce the boundary condition either.

But you do have to make sure that all reads outside the matrix bounds are also wrapped around properly. For example, for your Laplacian you could do something like this:

@jit(nopython=True, fastmath=True)
def laplacian_numba(field, dh2_inv, out):
    """
    Compute the laplacian of an array using a 5 point stencil.
    """
    # Note the new range: 0, ..., N-1
    for i in range(N):
        for j in range(N):
            out[i, j] = (
                field[(i   1) % N, j]
                  field[(i - 1) % N, j]
                  field[i, (j   1) % N]
                  field[i, (j - 1) % N]
                - 4 * field[i, j]
            ) * dh2_inv
    return out

Incidentally, the same could be written in a oneliner using four np.roll calls, but I don't know if it'd be faster than your numba approach.