Non-linear Poisson equation with gravitational term

1

Hi all,

I have a question regarding how to add a gravitational term in the one-dimensional non-linear Poisson equation in Fenics:
-\frac{\partial }{\partial x}(g(u)(\frac{\partial u}{\partial x} + 1)) = f,
where g(u) = u^2 + 1. This equation comes from porous media flow (Richards equation).

I tried some implementations using “Identity” and “derivative”, but could not figure out how to add the “1” in the equation (gravitational term) correctly.

Here is a minimum code example (modified from the non-linear Poisson in the tutorial). I appreciate your help.

from fenics import *
import numpy as np


def q(u):
    "Return nonlinear coefficient"
    return 1 + u**2

u_D = Expression('pow(x[0], 2) + 1', degree=1)
f = Expression('-10.0*pow(x[0], 4) - 4.0*pow(x[0], 3) - 12.0*pow(x[0], 2) - 4.0*x[0] - 4.0', degree=1)

mesh = UnitIntervalMesh(8)
V = FunctionSpace(mesh, 'P', 1)
u_D = Expression(u_code, degree=1)

def boundary(x, on_boundary):
    return on_boundary

bc = DirichletBC(V, u_D, boundary)
u = TrialFunction(V)
v = TestFunction(V)
u_k = interpolate(Constant(0.0), V)
# a = q(u_k)*(u.dx(0) + 1.0)*v.dx(0)*dx # explicit derivative (not working)
# a = inner(q(u_k)*(nabla_grad(u) + Identity(1)), nabla_grad(v))*dx # used nabla_grad (not working)
a = inner(q(u_k)*(nabla_grad(u)), nabla_grad(v))*dx # without "1" in the equation (works for that case)
L = f*v*dx

u = Function(V)
eps = 1.0
tol = 1.0E-5
iteration = 0
maxiter = 25
while eps > tol and iteration < maxiter:
    iteration += 1
    solve(a == L, u, bc)
    diff = u.compute_vertex_values(mesh) - u_k.compute_vertex_values(mesh)
    eps = np.max(np.abs(diff))
    print(f"Iteration {iteration} and norm {eps}")
    u_k.assign(u)


plot(u) # blue: computed solution
plot(interpolate(u_D, V)) # orange: true solution
2

Hi @Toshi,

The ufl_shape of the gravitational term is important. Consider the following:

from fenics import *

g1 = Constant(  1   )
g2 = Constant( (1,) )
g3 = Identity(1)

print(g1.ufl_shape)
print(g2.ufl_shape)
print(g3.ufl_shape)

yielding

()
(1,)
(1, 1)

As can be seen, Constant(1) is a scalar, Constant( (1,) ) is a one-component vector, and Identity(1) is a 1 \times 1 matrix. You need to choose a form that is compatible with the derivative operator that you are using. For nabla_grad (which returns a one-component vector), you should use:

g = Constant((1,))
a = inner(q(u_k)*nabla_grad(u), nabla_grad(v))*dx
L = (f*v - inner(q(u_k)*g, nabla_grad(v)))*dx

Note that I have corrected your variational form by moving the gravitational term to the RHS. As you had defined it, it caused an ArityMismatch, because nabla_grad(u) depends on the trial function, while Constant( (1,) ) does not.

3

Hi Connor,

Thank you very much for your detailed explanation with examples. It was very clear, and my code is now working well! Thank you!