Weak formulation for square of gradient

variational formulation
1

Hello

I’m a beginner in FEniCS and have been following the “Solving PDEs in Python” tutorial (which is excellent - thank you!).

Considering e.g. the simple Poisson example problem: -nabla^2 u(x) = f(x). How do you deal with terms like (nabla u(x))^2?

E.g. -nabla^2 u(x) + (nabla u(x))^2 = f(x), by which I mean in 1d -u''(x) + (u'(x))^2 = f(x). Is this possible? How do you write the weak form - do I need to integrate by parts?

My attempt: starting from the provided example which works:

from fenics import *

# Create mesh and define function space
mesh = UnitSquareMesh(8, 8) 
V = FunctionSpace(mesh, 'P', 1)

# Define boundary condition
u_D = Expression('1 + x[0]*x[0] + 2*x[1]*x[1]', degree=2) 
def boundary(x, on_boundary):
   return on_boundary 
bc = DirichletBC(V, u_D, boundary)

# Define variational problem

u = TrialFunction(V) 
v = TestFunction(V) 
f = Constant(-6.0) 
a = dot(grad(u), grad(v))*dx 
L = f*v*dx

# Compute solution
u = Function(V) 
solve(a == L, u, bc)

If i don’t integrate by parts, I should just modify it as follows:

a = dot(grad(u), grad(v)) * dx + dot(grad(u), grad(u)) * v * dx

This gives the error:

NotImplementedError: Cannot take length of non-vector expression.

Where am I going wrong? Or is it a conceptual mistake?

Thanks

1 Like
2

Welcome robert. First: You can use LaTeX for formulas using $ signs.

Concerning your specific problem:

-\operatorname{div}\operatorname{grad}u+\operatorname{grad} u\cdot \operatorname{grad} u=f

The weak form is as usual obtained through multiplying by test function v and integration (assuming homogenous Neumann boundary conditions)

\int_\Omega \operatorname{grad}v\cdot\operatorname{grad}u\,\text dx +\int_\Omega v(\operatorname{grad} u)^2 \,\text dx =\int_\Omega fv\,\text dx\,.

You are on the right track in your approach, but the syntax solve(a==L, ...) is used for linear problems and the one you want to solve is nonlinear. But we can simply define a nonlinear variational form as

F(u,v) = \int_\Omega \operatorname{grad}v\cdot\operatorname{grad}u\,\text dx +\int_\Omega v(\operatorname{grad} u)^2 \,\text dx -\int_\Omega fv\,\text dx

such that your solution has to satisfy F(u,v)=0. FEniCS will automatically compute the Jacobian and iteratively solve your problem when solve is called in the right way (cf. my answer here: Function() v.s. TrialFunction())

Consider the following code:

from fenics import *

# Create mesh and define function space
mesh = UnitSquareMesh(8, 8) 
V = FunctionSpace(mesh, 'P', 1)

# Define boundary condition
u_D = Expression('1 + x[0]*x[0] + 2*x[1]*x[1]', degree=2) 
def boundary(x, on_boundary):
   return on_boundary 
bc = DirichletBC(V, u_D, boundary)

# Define variational problem

u = Function(V) 
v = TestFunction(V) 
f = Constant(-6.0) 
F = dot(grad(u), grad(v)) * dx + dot(grad(u), grad(u)) * v * dx - f*v*dx

solve(F==0, u, bcs = bc)

For introductory reading into nonlinear problems in FEniCS I recommend
http://home.simula.no/~hpl/homepage/fenics-tutorial/release-1.0-nonabla/webm/nonlinear.html
as a first read.

3 Likes
3

Ah that was easy! I jumped ahead in the book too soon.

Thanks for your quick solution and suggested reading!