Having trouble with Nonlinear Poisson

1

Hi, I am having trouble solving the following equation:
-\nabla^2u + A = \frac{B}{u^2}
where A and B are constants.
This is calculated in the space between two spheres, one located within a larger one.
I set boundary conditions (zero) at surfaces, and I confirmed that the code works for regular electrostatic equation (right hand side of equation being constant instead of 1/u^2).
I also set g=0 as another boundary condition.

u = fn.Function(V)
v = fn.TestFunction(V)
f = B/u**2
g = fn.Constant(0)

a = (fn.inner(fn.grad(u),fn.grad(v))+ A*v)*fn.dx
L = f*v*fn.dx + g*v*fn.ds
F = a-L

fn.solve(F==0, u, bcs)

However, this returns an error from Newtonian solver that suggests the solution didn’t converge.
The output from Newtonian solver only shows r (abs) =±nan, so something is wrong.
I know the solution exists, since this is a well known physical model. Can someone point out what I’m doing wrong with this?

2

Well, this problem is not well defined, as it is only determined up to a constant when solving Poisson with pure Neumann conditions, see: Poisson with pure Neumann conditions

3

Hi dokken, Thanks for your comment.
I have both Dirichlet and Newmann boundary.

Following is an example of working code.

from dolfin import *

# Create mesh and define function space
mesh = UnitSquareMesh(32, 32)
V = FunctionSpace(mesh, "Lagrange", 1)

# Define Dirichlet boundary (x = 0 or x = 1)
def boundary(x):
    return x[0] < DOLFIN_EPS or x[0] > 1.0 - DOLFIN_EPS

# Define boundary condition
u0 = Constant(0.0)
bc = DirichletBC(V, u0, boundary)

# Define variational problem
u = Function(V)
v = TestFunction(V)
f = 1   # <-----------------want to use 1/u**2
g = Constant(0)
a = inner(grad(u), grad(v))*dx
L = f*v*dx + g*v*ds
F = a-L

# Compute solution
solve(F==0, u, bc)

# Plot solution
plot(u, interactive=True)

Does what I’m saying make sense or am I doing something completely wrong?

4

Your zero initial guess makes the problem singular. Using for instant u=project(Constant(1), V) as an initial guess solves the problem

1 Like
5

Now it worked! Thank you so much for your help, @dokken !

6

In case anyone else is interested, here is my final working example code:

from dolfin import *

# Create mesh and define function space
mesh = UnitSquareMesh(32, 32)
V = FunctionSpace(mesh, "Lagrange", 1)

# Define Dirichlet boundary (x = 0 or x = 1)
def boundary(x):
    return x[0] < DOLFIN_EPS or x[0] > 1.0 - DOLFIN_EPS

# Define boundary condition
u0 = Constant(0.0)
bc = DirichletBC(V, u0, boundary)

# Define variational problem
u = project(Constant(1),V)
v = TestFunction(V)
f = 1/u**2
g = Constant(0)
a = inner(grad(u), grad(v))*dx
L = f*v*dx + g*v*ds
F = a-L

# Compute solution
solve(F==0, u, bc)

# Plot solution
plot(u, interactive=True)

Also just a quick note: As it always does, initial guess requires to be within the similar order of magnitude to the final solution (the solution did not converge when the initial guess was off by factor of 10 for more complicated case when I tested). This tolerance is probably model-dependent, but I haven’t tested with other methods.