Assembling forms involving MixedElements

variational formulation
1

I am working on error estimation of time-dependent PDEs. How do I assemble a form (if all arguments are known a priori) to obtain a vector when using a MixedElement?

Consider a reaction-diffusion PDE system:

u’ = \laplacian u + u * v + f(x, t)
v’ = \laplacian v - u * v + g(x, t)

I used MixedElement to model the right-hand-side forms.

mesh       = RectangleMesh(Point(0.0, 0.0), Point(1.0, 1.0), 30, 30)
V             = FiniteElement("Lagrange", mesh.ufl_cell(), 1)
VSpace   = FunctionSpace(mesh, V*V)
q1, q2     = TestFunctions(VSpace)
z1, z2     = TrialFunctions(VSpace)
U            = Function(VSpace)
u,v          = U.split()
F1          = (-inner(grad(u),grad(q1)) + u*v*q1 + f_expr * q1) * dx
F2          = (-inner(grad(v),grad(q2)) -  u*v*q2 + g_expr * q2) * dx

My question is if I know a Function W in VSpace, how do I evaluate the form at W and obtain a result in VSpace? Is the following the correct way to compute F(W) where F = [F1; F2]?

W = Function(VSpace)
FW = assemble(replace(F1, {U: W}) + replace(F2, {U:W}))

Thank you

1 Like
2

The use of replace looks correct, although it’s probably giving you strange results due to the use of U.split() instead of split(U). See the following simple example:

from dolfin import *
mesh = UnitIntervalMesh(10)
V = FiniteElement("Lagrange", mesh.ufl_cell(), 1)
VSpace = FunctionSpace(mesh, V*V)
U = Function(VSpace)

u,v = split(U) # prints 1.0, which is the right answer for FW
#u,v = U.split() # prints 0.0, which is wrong

FU = u*dx

from ufl import replace
W = project(as_vector([1,1]),VSpace)
FW = replace(FU,{U:W})

print(assemble(FW))

(This has been known to be a confusing aspect of the API for a long time; see, e.g., the discussion here.)

3

Thank you. Suppose I extend your example as follows:

from dolfin import *
mesh = UnitIntervalMesh(10)
V = FiniteElement("Lagrange", mesh.ufl_cell(), 1)
VSpace = FunctionSpace(mesh, V*V)
U = Function(VSpace)

u,v = split(U) # prints 1.0, which is the right answer for FW
#u,v = U.split() # prints 0.0, which is wrong

FU = u*dx
FV = u*v*dx

from ufl import replace
W = project(as_vector([1,1]),VSpace)
FW = replace(FU,{U:W}) + replace(FV, {U:W})

print(assemble(FW))

The result is 2.

Edit: U and W vectors are 22 entries long. What I am looking for is FW to have 22 entries.

4

My example is just to illustrate the different effects of U.split() and split(U), and I chose a form of arity 0 that assembles to a scalar for simplicity. I believe what you want is to just replace

u,v = U.split()

with

u,v = split(U)

in your original example.

5

Thank you for your reply. I have replaced as you suggested. And I am still looking for a vector result of 22 entries. I come from a time integration background and I am looking to evaluate right-hand-side functions described using forms.

From my perspective, I look at U as representing the spatial discretized solution of the PDE at some time instant and I want to evaluate the right hand side function to get a slope estimate. How do I go about assembling the right hand side to obtain such a result? Thank you for your help thus far.

6

You might take a look at this demo

https://fenicsproject.org/docs/dolfin/latest/python/demos/cahn-hilliard/demo_cahn-hilliard.py.html

or this one

https://fenicsproject.org/docs/dolfin/latest/python/demos/elastodynamics/demo_elastodynamics.py.html

for examples of time stepping for unsteady problems in FEniCS.

7

Thank you. I will take a look at them.