Variational form of Nonlinear Poisson equation with mixed elements

Hello Everyone!,
I am trying to solve this system of equations:
\nabla{u_i}+u_i(\nabla{u_p})=0 here i=1,...n
&
\nabla.(\nabla{u_p})=-f\sum_{i=1}^{n}z_iu_i

here u_i and u_p are functions to be solved in a unit square domain.
I am trying to solve the system by using mixed element function such as described by the advection diffusion reaction tutorial. here is a mwe.

from fenics import *
#parameters inside variational form
z_a1 = 1
z_a2 = -1
f = Constant(43)
#Create mesh and define function space
mesh = UnitSquareMesh(8, 8)
degree = 1
P3 = FiniteElement('P', triangle, degree)
element = MixedElement([P3, P3, P3])
V = FunctionSpace(mesh, element)
a=d
(v_a1, v_a2, v_p) = TestFunctions(V)
u = Function(V)
(u_a1, u_a2, u_p) = split(u)

#Boundary connditions
tol = 1E-14
def boundary_R(x, on_boundary):
 return on_boundary and (near(x[0], 0, tol))

def boundary_L(x, on_boundary):
     return on_boundary and (near(x[0], 1, tol))
bc1 = DirichletBC(V.sub(0), Constant(1.0), boundary_R)
bc2 = DirichletBC(V.sub(1), Constant(1.0), boundary_L)
bc3 = DirichletBC(V.sub(2), Constant(-0.5), boundary_L)
bc4 = DirichletBC(V.sub(2), Constant(0.0), boundary_R)
bcs = [bc1, bc2, bc3, bc4]

#Variational Form residues residues
F_a1 = grad(u_a1) * v_a1 * dx + z_a1 * u_a1 * grad(u_p) * v_a1 * dx

F_a2 = grad(u_a2) * v_a2 * dx + z_a2 * u_a2 * grad(u_p) * v_a2 * dx

F_p = (-80)*dot(grad(u_p), grad(v_p)) * dx + (z_a1 * u_a1 + z_a2 * u_a2) * f * v_p * dx #poisson eqn

F = F_a1 + F_a2 + F_p

#Solver
solve(F == 0, u, bcs, solver_parameters={
                'nonlinear_solver': 'newton',
                'newton_solver': {
                    'linear_solver': 'mumps',
                    'maximum_iterations': 50,
                    'relative_tolerance': 1.0e-4}})

#Fetching the values
_u_a1, _u_a2, _u_p = u.split()
vtkfile_a1 = File('u1.pvd')
vtkfile_a1 << _u_a1
vtkfile_a2 = File('u2.pvd')
vtkfile_a2 << _u_a2
vtkfile_p = File('p.pvd')
vtkfile_p << _u_p

The error i am getting is

Can only integrate scalar expressions. The integrand is a tensor expression with value shape (2,) and free indices with labels ().
Traceback (most recent call last):
File “/home/PycharmProjects/MWE.py”, line 36, in
F_a1 = grad(u_a1) * v_a1 * dx + z_a1 * u_a1 * grad(u_p) * v_a1 * dx
File “/usr/lib/python3/dist-packages/ufl/measure.py”, line 406, in rmul
error("Can only integrate scalar expressions. The integrand is a "
File “/usr/lib/python3/dist-packages/ufl/log.py”, line 158, in error
raise self._exception_type(self._format_raw(*message))
ufl.log.UFLException: Can only integrate scalar expressions. The integrand is a tensor expression with value shape (2,) and free indices with labels ().
Process finished with exit code 1.

It seems like my variational formulation for the first equation is not right. I was wondering how the system can be solved with the formulation such as in the mwe with mixed element method. I would really appreciate any help.

Your variational form here is vectorial, i.e. it has two components. The same holds for the rhs of

which is a vectorial quantity.
Your problem can be split into two scalar PDEs, i.e.

\begin{align} \frac{\partial u_i}{\partial x} + u_i \frac{\partial u_p}{\partial x} &= 0 && i=1,\dots,n\\ \frac{\partial u_i}{\partial y} + u_i \frac{\partial u_p}{\partial y} &= 0 && i=1,\dots,n\ \end{align}

which for each component can be multiplied by a suitable scalar test function v_a1.

You can use u_a1.dx(0) to get the derivative in x direction and u_a1.dx(1) for the derivative in y-direction.

1 Like

Hi dokken, thanks for the reply. I have tried this way but I am still getting an error. revised script for mwe

from fenics import *

#parameters inside variational form
z_a1 = 1
z_a2 = -1
f = Constant(43)

#Create mesh and define function space
mesh = UnitSquareMesh(8, 8)
degree = 1
P3 = FiniteElement('P', triangle, degree)
element = MixedElement([P3, P3, P3])
V = FunctionSpace(mesh, element)
(v_a1, v_a2, v_p) = TestFunctions(V)

#Functions to be solved
u = Function(V)
(u_a1, u_a2, u_p) = split(u)

#Boundary connditions
tol = 1E-14
def boundary_R(x, on_boundary):
 return on_boundary and (near(x[0], 0, tol))

def boundary_L(x, on_boundary):
     return on_boundary and (near(x[0], 1, tol))


bc1 = DirichletBC(V.sub(0), Constant(1.0), boundary_R)
bc2 = DirichletBC(V.sub(1), Constant(1.0), boundary_L)
bc3 = DirichletBC(V.sub(2), Constant(-0.5), boundary_L)
bc4 = DirichletBC(V.sub(2), Constant(0.0), boundary_R)
bcs = [bc1, bc2, bc3, bc4]

#Variational Form residues residues
F_a1_x = u_a1 * v_a1 * dx(0) + z_a1 * u_a1 * u_p * v_a1 * dx(0)
F_a1_y = u_a1 * v_a1 * dx(1) + z_a1 * u_a1 * u_p * v_a1 * dx(1)

F_a2_x = u_a2 * v_a2 * dx(0) + z_a2 * u_a2 * u_p * v_a2 * dx(0)
F_a2_y = u_a2 * v_a2 * dx(1) + z_a2 * u_a2 * u_p * v_a2 * dx(1)

F_p = (-80)*dot(grad(u_p), grad(v_p)) * dx + (z_a1 * u_a1 + z_a2 * u_a2) * f * v_p * dx #poisson eqn

F = F_a1_x + F_a1_y + F_a2_x + F_a2_y + F_p

#Solver
solve(F == 0, u, bcs, solver_parameters={
                'nonlinear_solver': 'newton',
                'newton_solver': {
                    'linear_solver': 'mumps',
                    'maximum_iterations': 50,
                    'relative_tolerance': 1.0e-4}})

#Fetching the values
_u_a1, _u_a2, _u_p = u.split()

vtkfile_a1 = File('a1.pvd')
vtkfile_a1 << _u_a1
vtkfile_a2 = File('a2.pvd')
vtkfile_a2 << _u_a2
vtkfile_p = File('p.pvd')
vtkfile_p << _u_p

Now the error reads:

No Jacobian form specified for nonlinear variational problem.
Differentiating residual form F to obtain Jacobian J = F’.
Solving nonlinear variational problem.
Newton iteration 0: r (abs) = 1.105e+02 (tol = 1.000e-10) r (rel) = 1.000e+00 (tol = 1.000e-04)
Traceback (most recent call last):
File “/home/PycharmProjects/MWE.py”, line 47, in
solve(F == 0, u, bcs, solver_parameters={
File “/usr/lib/petsc/lib/python3/dist-packages/dolfin/fem/solving.py”, line 233, in solve
_solve_varproblem(*args, **kwargs)
File “/usr/lib/petsc/lib/python3/dist-packages/dolfin/fem/solving.py”, line 314, in _solve_varproblem
solver.solve()
RuntimeError:

*** -------------------------------------------------------------------------
*** DOLFIN encountered an error. If you are not able to resolve this issue
*** using the information listed below, you can ask for help at


*** fenics-support@googlegroups.com


*** Remember to include the error message listed below and, if possible,
*** include a minimal running example to reproduce the error.


*** -------------------------------------------------------------------------
*** Error: Unable to solve linear system using PETSc Krylov solver.
*** Reason: Solution failed to converge in 0 iterations (PETSc reason DIVERGED_PC_FAILED, residual norm ||r|| = 0.000000e+00).
*** Where: This error was encountered inside PETScKrylovSolver.cpp.
*** Process: 0


*** DOLFIN version: 2019.2.0.dev0
*** Git changeset: unknown
*** -------------------------------------------------------------------------
is this because i am making the residues of individual directional component sum to 0. Also if i am not mistaking solving the equations now gives the value of unkown u_i in x and y directions separately, whereas the boundary conditions that i have assigned are fixed drichlet conditions for the overall magnitude at the points (x,y).

You are misunderstanding

This means that your variational form

should be

F_a1 = u_a1.dx(0) * v_a1 * dx + u_a1 * u_p.dx(0) * v_a1 * dx
F_a1 += u_a1.dx(1) * v_a1 * dx + u_a1 * u_p.dx(1) * v_a1 * dx

etc.

That was silly of me, however even with the corrected form i still get the same error. now the variational form is:

#Variational Form residues residues
F_a1 = u_a1.dx(0) * v_a1 * dx + z_a1*u_a1 * u_p.dx(0) * v_a1 * dx
F_a1 += u_a1.dx(1) * v_a1 * dx + z_a1*u_a1 * u_p.dx(1) * v_a1 * dx

F_a2 = u_a2.dx(0) * v_a2 * dx + z_a2*u_a2 * u_p.dx(0) * v_a2 * dx
F_a2 += u_a2.dx(1) * v_a2 * dx + z_a2*u_a2 * u_p.dx(1) * v_a2 * dx
F_p = (-80)*dot(grad(u_p), grad(v_p)) * dx + (z_a1 * u_a1 + z_a2 * u_a2) * f * v_p * dx #poisson eqn

F = F_p + F_a1 + F_a2

#Solver
solve(F == 0, u, bcs, solver_parameters={
                'nonlinear_solver': 'newton',
                'newton_solver': {
                    'linear_solver': 'mumps',
                    'maximum_iterations': 50,
                    'relative_tolerance': 1.0e-4}})

and the error reads:

no Jacobian form specified for nonlinear variational problem.
Differentiating residual form F to obtain Jacobian J = F’.
Solving nonlinear variational problem.
Newton iteration 0: r (abs) = 1.105e+02 (tol = 1.000e-10) r (rel) = 1.000e+00 (tol = 1.000e-04)
Traceback (most recent call last):
File “/home/PycharmProjects/MWE.py”, line 46, in
solve(F == 0, u, bcs, solver_parameters={
File “/usr/lib/petsc/lib/python3/dist-packages/dolfin/fem/solving.py”, line 233, in solve
_solve_varproblem(*args, **kwargs)
File “/usr/lib/petsc/lib/python3/dist-packages/dolfin/fem/solving.py”, line 314, in _solve_varproblem
solver.solve()
RuntimeError:

*** -------------------------------------------------------------------------
*** DOLFIN encountered an error. If you are not able to resolve this issue
*** using the information listed below, you can ask for help at


*** fenics-support@googlegroups.com


*** Remember to include the error message listed below and, if possible,
*** include a minimal running example to reproduce the error.


*** -------------------------------------------------------------------------
*** Error: Unable to solve linear system using PETSc Krylov solver.
*** Reason: Solution failed to converge in 0 iterations (PETSc reason DIVERGED_PC_FAILED, residual norm ||r|| = 0.000000e+00).
*** Where: This error was encountered inside PETScKrylovSolver.cpp.
*** Process: 0


*** DOLFIN version: 2019.2.0.dev0
*** Git changeset: unknown
*** -------------------------------------------------------------------------

Prcess finished with exit code 1

This has to do with your problem not converging.

I would strongly suggest build your example up step by step, i.e.
start with a known u_p. Are you then able to find a unique and stable solution to u_i?
Solving PDEs of this form is not trivial.

Thanks for the suggestion, so I already know what the solution of u_i terms will look like:
u_i=exp(-zu_p) . eq1
its an exponential dependence on u_p. Now Instead of letting fenics get to this solution, I can introduce this eq1 directly into the script and then essentially solve for only the nonlinear Poisson equation as described in this nonlinear Poisson tutorial. After that I can simply project the values of u_i in the same mesh space to save the results. here is the script i came up with:

from fenics import *
import ufl
#parameters inside variational form
z_a1 = 1
z_a2 = -1
f = Constant(43)

#Create mesh and define function space
mesh = UnitSquareMesh(8, 8)
degree = 1
P3 = FiniteElement('P', triangle, degree)
V = FunctionSpace(mesh, P3)
(v_p) = TestFunction(V)

#Functions to be solved
u_p = Function(V)

#Boundary connditions
tol = 1E-14
def boundary_R(x, on_boundary):
 return on_boundary and (near(x[0], 1, tol))

def boundary_L(x, on_boundary):
     return on_boundary and (near(x[0], 0, tol))

x=SpatialCoordinate(mesh)

u_a1 = conditional(ufl.eq(x[0], tol), Constant(1.0), (exp(-z_a1 * u_p)))
u_a2 = conditional(ufl.eq(x[0], tol), Constant(1.0), (exp(-z_a2 * u_p)))

bc1 = DirichletBC(V, Constant(0.0), boundary_L)
bc2 = DirichletBC(V, Constant(0.5), boundary_R)
bcs = [bc1, bc2]

#Variational Form residues residues

F_p = (-80)*dot(grad(u_p), grad(v_p)) * dx + (z_a1 * u_a1 + z_a2 * u_a2) * f * v_p * dx #poisson eqn

F = F_p

#Solver
solve(F == 0, u_p, bcs, solver_parameters={
                'nonlinear_solver': 'newton',
                'newton_solver': {
                    'linear_solver': 'mumps',
                    'maximum_iterations': 50,
                    'relative_tolerance': 1.0e-4}})

#Fetching the values
u_a1 = project(u_a1, V)
u_a2 = project(u_a2, V)
vtkfile_a1 = File('a1.pvd')
vtkfile_a1 << u_a1
vtkfile_a2 = File('a2.pvd')
vtkfile_a2 << u_a2
vtkfile_p = File('p.pvd')
vtkfile_p << u_p 

Also since u_i is not fenics function to be solved directly, the standard drichlet boundary condition method wouldnt work, so i had to introduce the conditional forcing the values to be my desired bc value at one of the boundary and when its false (interior), the exponential function would kick in.