Singularity in the domain

[noirchen]

Hello Dr. Zwicker,

I am trying to simulate a spherical symmetric diffusion problem on $r\in[0,1]$

$$ \frac{\partial u}{\partial t}=\frac{2}{r}\frac{\partial u}{\partial r}+\frac{\partial^2u}{\partial r^2} $$

with BCs $u_r(r=0)=0$ and $u(r=1)=1$ and IC $u(t=0)=0$. There is a singularity at $r=0$, how can I simulate this? My implementation attempt is

from pde import PDE, CartesianGrid, MemoryStorage, ScalarField, plot_kymograph

# Definition of the PDE
eq = PDE({"u": f"laplace(u)+2*gradient(u)/x"},
         bc=[{
             "derivative": 0
         }, {
             "value": 1
         }])

grid = CartesianGrid([[0, 1]], 100)  # generate grid
field = ScalarField(grid, 0)  # generate initial condition

storage = MemoryStorage()  # store intermediate information of the simulation
res = eq.solve(field, 5, dt=0.1, tracker=storage.tracker(0.5))  # solve the PDE

plot_kymograph(storage)  # visualize the result in a space-time plot

Any help is appreciated!

[david-zwicker]

It looks to me as if you want to use the laplace operator in spherical coordinates. If so, the situation might be simpler than you even thought. You can simply use a SphericalSymGrid, which implements the correct laplacian:

from pde import PDE, SphericalSymGrid, MemoryStorage, ScalarField, plot_kymograph

# Definition of the PDE
eq = PDE({"u": f"laplace(u)"}, bc=[{"derivative": 0}, {"value": 1}])

grid = SphericalSymGrid(1, 32)  # generate grid
field = ScalarField(grid, 0)  # generate initial condition

storage = MemoryStorage()  # store intermediate information of the simulation
res = eq.solve(field, 1, dt=0.001, adaptive=True, tracker=storage.tracker(0.01))  # solve the PDE

plot_kymograph(storage)  # visualize the result in a space-time plot

Note that I also adjusted some other parameters to create a nicer plot. In particular, I used adaptive time stepping to prevent instabilities.

[noirchen]

Thanks a lot. This is such a nice tool!

[outspeckle]

Is there a way to deal with this problem in non-radial-symmetric cases? I'm trying to solve the radial diffusion equation with py-pde by transforming it to a rectangular grid r x phi, but the Laplacian obviously has a singularity at r=0. Is it possible to shift the grid points at the origin to omit r=0?

[david-zwicker]

I'm not exactly sure what you're trying to achieve, but if you want to use a CartesianGrid, you can avoid having a grid point at the center by using an even number of support points. In principle, we can also implement spherical grids with angular dependencies. I think I have some preliminary code for polar (2d) grids with angular discretization, and I don't see any reason why this cannot be generalized to 3d spherical grids, but I don't currently have the time to do it myself.