NeuroDiffGym / neurodiffeq

A library for solving differential equations using neural networks based on PyTorch, used by multiple research groups around the world, including at Harvard IACS.
http://pypi.org/project/neurodiffeq/
MIT License
664 stars 87 forks source link

Nonzero Dirichlet boundary conditions #177

Closed ma-sadeghi closed 8 months ago

ma-sadeghi commented 1 year ago

I'm trying to solve the Laplace PDE in 2D with the following boundary conditions:

x = 0 ; u = 1 x = 1 ; u = 0 y = 0 ; u = 0 y = 1 ; u = 0

I'm currently using the following snippet to enforce these BCs:

BCs = [
    DirichletBVP2D(
        x_min=xmin, x_min_val=lambda y: 0,
        x_max=xmax, x_max_val=lambda y: 0,
        y_min=ymin, y_min_val=lambda x: 0,
        y_max=ymax, y_max_val=lambda x: 0,
    )
]

but it seems that the BCs at y_min and y_max are being ignored. From the looks of the solution, it seems that at y_min and y_max, a Neumann BC has been enforced (du/dy=0).

sathvikbhagavan commented 8 months ago

The boundary conditions are not continuous

ma-sadeghi commented 8 months ago

You mean the 4 corners?

sathvikbhagavan commented 8 months ago

The boundary conditions you mention in:

x = 0 ; u = 1
x = 1 ; u = 0
y = 0 ; u = 0
y = 1 ; u = 0

It is not continuous at (0, 0) as it can be 0 or 1.

But the boundary conditions in your snippet of code:

BCs = [
    DirichletBVP2D(
        x_min=xmin, x_min_val=lambda y: 0,
        x_max=xmax, x_max_val=lambda y: 0,
        y_min=ymin, y_min_val=lambda x: 0,
        y_max=ymax, y_max_val=lambda x: 0,
    )
]

should work.

ma-sadeghi commented 8 months ago

Thanks for your quick reply. I must have made a typo, the BCs in the snippet enforce u=0 everywhere, which I don't want. I want to enforce u=1 at x=0, and I understand this makes the corners discontinuous. Is there a way to ignore those two points, like would the following snippet work?

BCs = [
    DirichletBVP2D(
        x_min=xmin, x_min_val=lambda y: return 0 if y in [0, 1] else 1,
        x_max=xmax, x_max_val=lambda y: 0,
        y_min=ymin, y_min_val=lambda x: 0,
        y_max=ymax, y_max_val=lambda x: 0,
    )
]
sathvikbhagavan commented 8 months ago

I don't think that will work. Why not use some continuous and approximate step function like heaviside - https://en.wikipedia.org/wiki/Heaviside_step_function?

ma-sadeghi commented 8 months ago

Ya, that's right, I didn't mean the exact snippet, but conceptually. So, using NumPy's heaviside (in a way that mimics lambda y: return 0 if y in [0, 1] else 1) would work? Thanks!

sathvikbhagavan commented 8 months ago

So, using NumPy's heaviside (in a way that mimics lambda y: return 0 if y in [0, 1] else 1) would work?

Directly using heaviside might not work with AD (in pytorch - https://pytorch.org/docs/stable/generated/torch.heaviside.html). But using an approximate continuous form should work.

ma-sadeghi commented 8 months ago

Gotcha thanks!