PDE system simple example

1

I would like to solve a PDE system, where some quantities are vectors and others are scalars.
As an example I tried to do the following variant of the Poisson equation:

q = -k * grad(T)
div(q) = f

Here is what I tried:

from fenics import *

mesh = UnitIntervalMesh(10)
Vscalar = FunctionSpace(mesh, "P", 2)
Vvector = VectorFunctionSpace(mesh, "P", 2)
vs = TestFunction(Vscalar)
vv = TestFunction(Vvector)
T = TrialFunction(Vscalar)
q = TrialFunction(Vvector)

k = Constant(1.0)
f = Constant(6.0)
F = dot(k*grad(T) + q, vv) * dx + dot(div(q) - f, vs) * dx 

bc = DirichletBC(Vscalar, Constant(0.0), "on_boundary")

u = (Function(Vscalar), Function(Vvector))
solve(F == 0, u, bc)

This does not work however, my definition of u is not correct.

2

I think that you are already close to the correct solution. Instead of

use

Vscalar = FiniteElement("Lagrange", mesh.ufl_cell(), 2)  # use quadratic temperature triangle
Vvector = VectorElement("Lagrange", mesh.ufl_cell(), 1)  # use linear flux
Vmixed = FunctionSpace(mesh, Vscalar * Vvector)

T, q = TrialFunctions(Vmixed)
vs, vv = TestFunctions(Vmixed)

Be aware that for mixed FEM the function spaces have to “match”, otherwise you might get into trouble.

For boundary conditions use:

bc = DirichletBC(Vmixed.sub(0), Constant(0.0), "on_boundary")

Finally, solve

mixed_u = Function(Vmixed)
u, q = mixed_np1.split()
solve(lhs(F) == rhs(F), mixed_u, bc)

For more details see the tutorial on mixed formulation of the poisson equation.

1 Like
3

Thanks very much! I updated my code:

from fenics import *

mesh = UnitIntervalMesh(4)
el_scalar = FiniteElement("Lagrange", mesh.ufl_cell(), 2)  # use quadratic temperature triangle
el_vector = VectorElement("Lagrange", mesh.ufl_cell(), 1)  # use linear flux
el = MixedElement([el_scalar, el_vector])
Vmixed = FunctionSpace(mesh, el)
vs, vv = TestFunctions(Vmixed)
T, q   = TrialFunctions(Vmixed)

Vscalar = Vmixed.sub(0)
Vvector = Vmixed.sub(1)

k = Constant(1.0)
f = Constant(2.0)
F = dot(k*grad(T) + q, vv) * dx + dot(q, grad(vs)) * dx - f * vs * dx

bc = DirichletBC(Vscalar, Constant(0.0), "on_boundary")

u = Function(Vmixed)

solve(lhs(F) == rhs(F), u, bc)

T, q = u.split()

import numpy as np
np.asarray(T.vector())

array([  0.00000000e+00,             -inf,  -5.00000000e-01,
                    inf,  -1.00000000e+00,             -inf,
         1.15849359e-16,              inf,             -inf,
         5.00000000e-01,              inf,   0.00000000e+00,
                    inf,   1.00000000e+00])

I have a couple of follow up questions:

  • The solution contains inf and is nonsense. Is there maybe an issue with the boundary condition?
  • If I change the solve into solve(F == 0, u, bc) I get an error:
ArityMismatch: Adding expressions with non-matching form arguments ('v_0',) vs ('v_0', 'v_1').

why?

Edit: I opened a new question about the ArityMismatch.

4

If I change the order of the elements to

el_scalar = FiniteElement("P", mesh.ufl_cell(), 1) # was 2 before
el_vector = VectorElement("P", mesh.ufl_cell(), 1)

The solution looks sane:

image