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
```