DedalusProject / dedalus

A flexible framework for solving PDEs with modern spectral methods.
http://dedalus-project.org/
GNU General Public License v3.0
491 stars 115 forks source link

Laplace equation in 2D Dirichlet BCs #160

Open pavaninguva opened 2 years ago

pavaninguva commented 2 years ago

Hi, I am trying to implement the Laplace equation in 2D as follows: d2f/dx2 + d2f/dy2 = 0, with Dirichlet boundary conditions: f(x,1.0) = -1.0, f(0.0,y) = f(1.0,y) = f(x,0.0) = 0.0.

Using the Poisson script already developed for BVPs in the github repo threw the following error:

Traceback (most recent call last):
  File "poisson.py", line 26, in <module>
    problem.add_bc("u(y='left') = sin(8*x)")
  File "/home/pavan/miniconda3/envs/dedalus/lib/python3.8/site-packages/dedalus/core/problems.py", line 139, in add_bc
    return self.add_equation(*args, **kw)
  File "/home/pavan/miniconda3/envs/dedalus/lib/python3.8/site-packages/dedalus/core/problems.py", line 131, in add_equation
    self._build_object_forms(temp)
  File "/home/pavan/miniconda3/envs/dedalus/lib/python3.8/site-packages/dedalus/core/problems.py", line 152, in _build_object_forms
    temp['LHS'] = field.Operand.parse(temp['raw_LHS'], self.namespace, self.domain)
  File "/home/pavan/miniconda3/envs/dedalus/lib/python3.8/site-packages/dedalus/core/field.py", line 102, in parse
    expression = eval(string, namespace)
  File "<string>", line 1, in <module>
TypeError: 'Field' object is not callable

However, when I adapt the script and specify the boundary conditions as add_equations instead:

import os
import numpy as np
import matplotlib.pyplot as plt

from dedalus import public as de
from dedalus.extras import plot_tools

# Create bases and domain
x_basis = de.Fourier('x', 100, interval=(0, 1))
y_basis = de.Chebyshev('y', 100, interval=(0, 1))
domain = de.Domain([x_basis, y_basis], grid_dtype=np.float64)

# Poisson equation
problem = de.LBVP(domain, variables=['u','uy'])
problem.add_equation("dx(dx(u)) + dy(uy) = 0")
problem.add_equation("uy - dy(u) = 0")
problem.add_equation("left(u) = 0")
problem.add_equation("right(u) = 1")

# Build solver
solver = problem.build_solver()
solver.solve()

# Plot solution
u = solver.state['u']
u.require_grid_space()
plot_tools.plot_bot_2d(u)
plt.savefig('laplace.png')

Which now works and gives the following: laplace

I am now only able to specify the desired boundary conditions on the top and bottom of the domain and it looks like no-flux BCs are enforced on the left and right sides of the domain. I understand that the boundary condition is enforced along the y coordinate since I read somewhere that the Chebyshev basis function is the axis along which the boundary condition is imposed on in Dedalus. In this case, I am not entirely sure how to go about imposing the right Dirichlet BCs to satisfy the problem.

For additional info. I followed the installation procedure on the documentation exactly.

afraser3 commented 2 years ago

Hi, If I understand correctly, it sounds like you are trying to impose BCs along both axes. This feature is not included in the current Dedalus release, but I do believe it is an upcoming feature.

EDIT: To clarify, on the left and right, periodic boundary conditions are being used (a consequence of it being a Fourier basis), not no-flux boundary conditions.

pavaninguva commented 2 years ago

Ah, so in that case:

  1. Solving the problem with the specification of Dirichlet boundary conditions on all four domains (i.e. the model I want to solve) is not possible at the moment. Is there any indication when this might be out?
  2. Is it possible to employ a different basis function to atleast specify Neumann BCs along the left and right boundaries. I am guessing that specifying a Fourier basis results in periodic BCs along that axis?
  3. What is the nature of the issue running the orignal poisson problem on the github repo. It seems that the add_bc formulation is much cleaner than adding additional equations in terms of formulating the problem. It seems that running other example scripts such as the KDV script or 1-D Lane Emden equation had no issues.

Thanks!

jasonltorchinsky commented 1 year ago

I wanted to follow up on this: is it now possible to specify Dirichlet boundary conditions along both axes? If not, is it possible to specify Dirichlet boundary conditions along one axis and Neumann boundary conditions along another?

iamlll commented 1 year ago

Is it still the case that Dirichlet boundary conditions can only be imposed on one axis?

kburns commented 1 year ago

Doing this correctly for general PDEs is still an open research problem, but we're working on it (see this recent paper), and there will be a separate release and additional examples when it's supported in the code, so stay tuned!

HennesHajduk commented 11 months ago

Just a quick follow-up: Is the feature with such boundary conditions currently not supported at all (in the sense that such codes will not run properly or return undefined behavior) or is it generally possible to set such problems up and only be able to execute them in serial? The poisson script can indeed be modified in this sense but my question is basically: Is the fact that it's working a coincidence?

kburns commented 11 months ago

This issue is that, like all Dedalus models, enforcing boundary conditions requires tau terms, and these get complicated when you have different boundaries meeting at edges/corners. We understand how to do it in simple cases (e.g. Poisson) as described in the paper linked above, but not in more general cases. And we simply can't offer general support when doing it correctly and reliably it's an open/active mathematical problem.

kburns commented 11 months ago

Also yes, even with the correct tau terms for a given equation set, double-Chebyshev problems will in general be fully coupled and non easily parallelized. Parallelizing them would require using preconditioned iterative solvers or distributed direct solvers. These are both interesting routes to explore, but again they are research grade questions in their own right.

HennesHajduk commented 10 months ago

Thank you for the clarification. Of course, I understand that this issue, being an open research problem, is beyond the scope of a general purpose library can be expected to handle.