Handling of periodic boundary conditions?

[aaron-shih]

I'm wondering if I've missed something, but I can't seem to find/figure out how to implement periodic boundary conditions (or just define a periodic grid). With my understanding of finite differences, would it just be not supplying boundary conditions at all? If so, then my issue is that I can't seem to apply an equation over the entire grid/domain - the edges seems to never get updated even if I don't restrict the equation to the subdomain grid.interior.

[rhodrin]

Depending on what exactly the problem is you're solving, there may be different ways you want to set this up, but if, for example, you want a simple periodic condition on the x boundary you can use basic Eq's to achieve this, e.g.:

bc_p_x = [Eq(f[t+1,0,y], my_bv_function)]
bc_p_x +=  [Eq(f[t+1,-1,y], f[t+1,0,y])]

[aaron-shih]

I've tried that before I think? Anyhow, with this:

stencil = solve(eq,u.forward)
update = Eq(u.forward, stencil, subdomain = grid.interior)
pbc = [Eq(u[t+1, 0, y], u[t+1, nx-1, y])]
pbc += [Eq(u[t+1, x, ny-1], u[t+1, x, 0])]
op = Operator([update]+pbc)

I get a border, which I'm guessing is because the update doesn't update the boundary values due to the subdomain.

Then I guess my question is now, what would I replace subdomain = grid.interior with? Removing it entirely causes devito to output blank images, which I'm guessing is due to something diverging?

[mloubout]

What is your space_order?

[aaron-shih]

Currently, I have u as a TimeFunction with space_order=2

I can't use subdomain=grid.domain, but I can use subdomain=grid._subdomains[0] which I'm pretty sure is the entire Domain? However:

1. Solving my equation in this Domain without specifying boundary conditions still gives the hazy border around the edges
2. Solving my equation in this Domain with periodic boundary conditions as above results in Nan values

[rhodrin]

@aaron-shih Did you take a look at the c-code being produced (print(op.ccode)) - how does the algorithm being produced differ from what you want?

If you're able to create a 'minial' script (e.g. simple 1d problem) illustraiting what you'd like to achieve that would also help us guide you towards producing a suitable implementation.

[aaron-shih]

I will try parsing through the c-code, but I am currently unsure what to make of the floats named r##.

I'm currently trying to simulate phase separation via the Cahn-Hilliard Equation, if that helps, but regardless of the system I try (Burger's, heat equation), the boundary isn't behaving as I intend.

In essence, I would like the derivative at a boundary (say at index 0) to be calculated as (using a central difference): u_dx[0] = (u[1] - u[-1])/2 And a corresponding second derivative: u_dx2[0] = u[1] - 2u[0] + u[-1]

So for a 1D heat equation, I'd like the equation to be: Eq(u.dt, u_dx2) Such that it updates as: u[x] += u[(x+1)%nx] - 2*u[x] + u[(x-1)%nx] though u[(x-1)%nx] can just be u[x-1] due to how negative indices work in python.

Below a python script that models a 1D heat equation with PBC, with the initial condition as a step function. Heat is allowed to diffuse from the right boundary to the left boundary:

Script:

import numpy as np
import matplotlib.pyplot as plt

N = 64
dt = 0.02
T = 200

#Initialize as step function
u = np.zeros(N)
u[32:] = 1

for i in range(int(T/dt)+1):
    if i%1000 == 0:
        plt.plot(u)
        plt.ylim(0,1)
        plt.savefig('heat'+str(i*dt)+'.png')
        plt.close()
    for x in range(len(u)):
        u[x] += dt*(u[(x+1)%N] - 2*u[x] + u[(x-1)%N])

[smcantab]

@dabiged may be able to help with this issue, he recently contributed a PR #1368 on a channel flow example that uses periodic boundary conditions (PBC). Though he points out at the end of his example here that the hand-written PBC do not work correctly.

[aaron-shih]

I just tried the heat equation with similar periodic boundary conditions as mentioned in PR #1368 , but I do not get the desired result:

eq = Eq(u.dt, u.dx2, subdomain=grid.interior)
stencil = solve(eq,u.forward)
update = Eq(u.forward, stencil)
pbc = [Eq(u[t+1, 0], u[t+1, N-2])]
pbc += [Eq(u[t+1, 1], u[t+1, N-1])]
op = Operator([update]+pbc)

I am suspecting that the reason why periodic boundary conditions do not work is due to these lines:

pbc = [Eq(u[t+1, 0], u[t+1, N-2])]
pbc += [Eq(u[t+1, 1], u[t+1, N-1])]

which seem to translate to:

u[t+1, 0] = u[t+1, N-2]
u[t+1, 1] = u[t+1, N-1]

when in reality, we wwere looking for something like:

u[t+1, 0], u[t+1, N-2] = u[t+1, N-2], u[t+1, 0]
u[t+1, 1], u[t+1, N-1] = u[t+1, N-1], u[t+1, 1]

But I still don't think this is an accurate representation/implementation of periodic boundary conditions

[dabiged]

I too could not get the periodic boundary condition working in devito. The channel cfd notebook explicitly states this at the bottom. If you look at the last plot, the parallel flow is actually divergent near the periodic boundaries indicating something is wrong.

I tried pinning the last 10 cells to the first 10 and I still got anomalous results. I have no idea how to fix this.

[mloubout]

It think this is partially related to #1250 in the sense that while your implementation of the periodic BC is correct, the FD actually accesses and uses the value in the FD halo rather than the periodic condition. You could try to bypass it with runtime args op.apply(x_m=1, x_M=N-2) with the periodic BCs at x=0 and x=N-1 so that the halo is not used.