Laplace on unit square with boundary conditions using PDE

[scbarton]

Hi, I was happy to discover py-pde and am just getting started. I ran into some code that I think should work but does not. If anyone can help me see my mistake I'd really appreciate it:

from pde import *

# grid
N= 64
grid = CartesianGrid([[0, 1],[0, 1]], N)

# initial condition
state = ScalarField(grid, 1)

# boundary condition 
bc_x  = [{"value": 1},       {"value": 0} ]
bc_y  = ["natural",  "natural" ]
bcs = [bc_x, bc_y]

# governing equation
eq = PDE({"C": "laplace(C)"}, bc=bcs)

# solve
result = eq.solve(state, t_range=0.1, dt=0.001)

The result should be a linear decrease in C from 1 at the left to 0 at the right and uniform in the y-direction. The code above produces an error, NameError: name 'laplace' is not defined. Running without the bc=bcs argument returns the trivial result C=1. Using bc_y = [{"derivative":0}, {"derivative":0}] seems to blow up.

I know I can use solve_laplace_equation for this, but I want to use this approach for more complex user-defined problems.

[david-zwicker]

I'm not quite sure what equation you want to solve, but solve_laplace_equation actually solves the Laplace equation 0=laplace(c), while your code is attempting to solve ∂c/∂t = laplace(c), which is quite different.

In any case, the error NameError: name 'laplace' is not defined revealed a deeper problem in the code, where erroneous boundary conditions were not properly reported when using PDE. This will be fixed in the master branch very soon, but as a work-around, you can here already obtain a running version by setting bc_y = "natural". In essence, defining a natural boundary condition for just one side of an axis (in this case the y-axis) does not make much sense and the syntax you attempted to use is thus not supported.

[scbarton]

Thanks for the quick response. I'm mainly interested in the steady-state result, but am using the transient equation ∂c/∂t = laplace(c) in PDE to converge to steady state because I have not found another way to solve systems of equations at steady state. Is there another approach? Again, I am using laplace(c) as a simple starting point and plan to move on to more complex, coupled systems of equations.

Using version 0.13.5, and setting bc_y="natural", the calculation does not converge (see the result below). I've also included the expected result obtained using solve_laplace_equation() . Is there anything else I can do to get this to work?

from pde import *

# grid
N= 64
grid = CartesianGrid([[0, 1],[0, 1]], N)

# initial condition
state = ScalarField(grid, 1)

# boundary condition 
bc_x  = [{"value": 1},       {"value": 0} ]
bc_y  = "natural"
bcs = [bc_x, bc_y]

# governing equation
eq = PDE({"C": "laplace(C)"}, bc=bcs)

# solve
result = eq.solve(state, t_range=0.1, dt=0.001)

# plot result
result.plot()

res = solve_laplace_equation(grid, bcs)
res.plot()

[david-zwicker]

I don't think there is a reason to believe the system should converge – Instead it seems as if it is blowing up to infinity. You might instead want to try investigating ∂c/∂t = -laplace(c), but whether your system actually converges to some stationary state depends on whether this state is actually stable and whether you start with the right initial condition.

[scbarton]

It does converge. This is just the heat equation, so assuming a forward differencing scheme in time, and applying Von Neumann stability analysis, it should converge for ∆t ≤ ∆x^2 / 2.

For N=64, that gives ∆t ≤ 0.00012. The calculation actually does not converge for dt=1e-4, but it does converge for dt=1e-5. Here's the working code with a steady-state tracker included.

from pde import *

# grid
N= 64
grid = CartesianGrid([[0, 1],[0, 1]], N)

# initial condition
state = ScalarField(grid, 1)

# boundary condition 
bc_x  = [{"value": 1},       {"value": 0} ]
bc_y  = "natural"
bcs = [bc_x, bc_y]

# governing equation
eq = PDE({"C": "laplace(C)"}, bc=bcs)

# solve
result = eq.solve(state, t_range=1, dt=1e-5, tracker="steady_state")

# plot result
result.plot()

I knew this should work because it worked well using grid = UnitGrid([N, N]) and dt = 0.1. I just needed a small enough dt for the CartesianGrid case.

[david-zwicker]

You're right – I must have misread your original post. It's great that you figured out that lowering the time step resolved the instability. Unfortunately, this is difficult to automate, but when you skip specifying the time step, the framework uses an adaptive method from scipy, which might give you stable results.