Convergence problems in non lineal Poisson equation

General
1

Hi,
I’m trying to solve non-lineal Poisson’s equation with Dirichlet boundary conditions and what can be understood as two spherical sources. I am basically interested in the derivative of the solution. This is a simpler version of my code:

import numpy as np
from dolfin import *
import matplotlib.pyplot as plt
from numpy.random import random

# Sub domain for Dirichlet boundary condition. It is a boolean function
def boundary_boolean_function(x, on_boundary):
    return on_boundary

"""
PARAMETERS
"""

Lam = 1. 
mass = 20. 
rs = 0.6 # Radius of the source
mpl = 1.
G = 1.0 

m_ = mass * 4./3. * np.pi * rs**3

rstar = 1/Lam * (G * m_/(4 * np.pi * mpl))**0.5

# Function to deffine a mesh that follows a certain function
def symmetric_mesh3d(Nx, Ny, Nz, rstar):
    # Create a cubic mesh
    xmax = 11.
    m = BoxMesh(Point(-xmax, -xmax, -xmax), Point(xmax, xmax, xmax), Nx, Ny, Nz)
    x = m.coordinates()
    # Modify the coordinates as desired
    for i in range(len(x)):
        if x[i, 0] != 0:
            x [i, 0] = x [i, 0]
        if x[i, 1] != 0:
            x [i, 1] = x [i, 1]
        if x[i, 2] != 0:
            x [i, 2] = x [i, 2]
    return m

# Define the mesh that we are going to use
mesh = symmetric_mesh3d(19, 19, 19, 6.) # The mesh should be more refined near the radius of the source

V = FunctionSpace(mesh, "CG", 2) # Second order: we should get goood results up to the first derivative.

# Define boundary condition
g = Constant(1.) # Impose that the field has a constant value in the limits
bc = DirichletBC(V, g, boundary_boolean_function)

# Define variational problem

# Initial guess
u = interpolate(Constant(-1.), V)
v = TestFunction(V)

f = Expression("-20.*(sqrt(pow(x[1] - 1.05, 2) + pow(x[0], 2) + pow(x[2], 2)) < 0.6) - 20.*(sqrt(pow(x[1] + 1.05, 2) + pow(x[0], 2) + pow(x[2], 2)) < 0.6)", degree=2) # Density

F = inner((1./sqrt(1 + Lam**(-4)*inner(grad(u),grad(u))/2))*grad(u), grad(v))*dx - f*v*dx 

# Compute solution
M = u * dx
du = TrialFunction(V)
J = derivative(F, u, du) # Jacobian

solve(F == 0, u, bc, J=J,
    solver_parameters = {
    "newton_solver": {
        "linear_solver": "gmres",
        "preconditioner": "ilu",
        "relative_tolerance": 1e-6,
        "absolute_tolerance": 1e-6
    }
})

# Print a 2d slice of the solution
z = 0.0

N = 101
v_slice = np.zeros((N, N), dtype = "float")
xsol = np.zeros(N, dtype = "float")

x = mesh.coordinates()

for i in np.arange(N):
    for j in np.arange(N):
        x0 = -x[len(x)-1, 1] + 2 * x[len(x)-1, 1]/N*i
        y0 = -x[len(x)-1, 1] + 2 * x[len(x)-1, 1]/N*j
    
        v_slice[i, j] = u(x0, y0, z)

plt.figure()
plt.imshow(v_slice)
plt.colorbar()

plt.figure()
plt.plot(v_slice[:, int(N/2)])

# Derivative

x = mesh.coordinates()

vectorsol = np.zeros(N, dtype = "float")
xsol = np.zeros(N, dtype = "float")

for i in np.arange(N):
    z = - x[len(x)-1, 1] + 2 * x[len(x)-1, 0]/N*i # We start in 0 (center of the mesh)
    vectorsol[i] = u(z, 0., 0.)
    xsol[i] = z

vector_mod = vectorsol/(Lam**2)
r_mod = xsol

# PLOT THE DERIVATIVE (this is the important one)

der = np.zeros((N-1, N-1), dtype = "float")
h =  2 * x[len(x)-1, 1]/N
for i in np.arange(1, N-1):
    for j in np.arange(1, N-1):
        x0 = -x[len(x)-1, 1] + h*i
        y0 = -x[len(x)-1, 1] + h*j
        derx = (v_slice[i+1, j] - v_slice[i-1, j]) / (2 * h) # This is a first order method. I should probably improve it
        dery = (v_slice[i, j+1] - v_slice[i, j-1])/ (2 * h)
        derz = ( u(x0, y0, h) - u(x0, y0, -h) )/(2*h) # Second order method
        der[i, j] = derx ** 2 + dery ** 2 + derz ** 2
        #der[i, j] = dery**2

plt.figure()
plt.title('X')
plt.imshow(der)
plt.colorbar()

# Plot X in one axis
mind = int((N - 1)/2)

X_x = np.zeros((N-1), dtype = "float")

for i in np.arange(N-1):
    X_x[i] = der[mind, i]

plt.figure()
plt.title('X along the axis that connects the masses')
plt.yscale("log")
plt.plot(X_x[1:])
plt.show()

As far as I’ve been able to notice, the implementation does not seem wrong. However, I am facing multiple convergence problems. First of all, if I increase the number of elements in the mesh by as little as 1 or 2, the code fails to coverge. This is a big problem, as a bigger number of elements would be, in principle, needed in order to properly deffine the matter distribution. Another issue is that a change in the number of elements seems to change the value of the derivative of the solution by even two orders of magnitude. I have reviewed the code over a thousand times and played endlessly with the parameters, but I haven’t been able to come to a solution to this problem.
Is my approach breaking down somewhere? What can I do in this situation?
Thank you so much! :slight_smile: