zwicker-group / py-pde

Python package for solving partial differential equations using finite differences.
https://py-pde.readthedocs.io
MIT License
412 stars 52 forks source link

Cannot specify full boundary conditions in biharmonic equations without periodic boundary conditions #309

Closed tranqui closed 2 years ago

tranqui commented 2 years ago

I am trying to solve 4th order ODEs containing a biharmonic operator. The documentation describes solving the Kuramoto-Sivashinsky equation and the Cahn-Hilliard equation.

However, there is no option to specify enough boundary conditions to solve a 4th order ODE. For a linear ODE we need 4 and this is generally true for nonlinear ODEs also. With Neumann BCs I expect to be able to specify values and derivatives at the same point, with e.g. something like for a 2d system: bcs = [{'value': -1, 'derivative': 0}, {'value': 1, 'derivative': 0}] However, this always throws an error saying this feature is not supported. So I am curious how 4th order ODEs are listed as examples. Ultimately I'm interested in something similar to the Cahn-Hilliard equation so look to that as an example.

Looking at the implementation of the Cahn-Hilliard equation, it looks like this would only work with periodic boundary conditions because this implicitly sets values and all derivatives. That being said, I am never able to get convergent results for e.g. the binodal so I suspect either I am using py-pde wrong or there is a problem with the underlying solver. Below is a minimal working example that illustrates this where I try to obtain the binodal. I start with two interfaces around x=+-0.5 to respect periodic boundaries at x=+-1. The solver quickly fails with a "Field was not finite" error:

from pde import CahnHilliardPDE, ScalarField, CartesianGrid
grid = CartesianGrid([[-1, 1], [-1, 1]], [1000,50])
state = ScalarField.from_expression(grid, 'tanh(cos(pi*x)/0.1)')
eq = CahnHilliardPDE()
result = eq.solve(state, 1, dt=1e-6, method='implicit')
result.plot()

My expected behaviour is that it should converge to tighten up the interfaces with the correct interface width (rather than the dummy value of 0.1 I gave as an IC).

Is this a bug, or am I using it wrong? If the latter could you please advise me on how I should use the Cahn Hilliard (and similar) solvers? Thank you.