Mixed kinetic problem

Question

I'm working on a coupled 'kinetics' kind of problem which appears fairly simple, but I'm having trouble implementing it into FEniCS. I can also imagine stability issues but I'm not at that point yet. I could go into what's what, but since all parameters are dimensionless it does not seem relevant. Do note that $t$ is time and $x$ is a spatial position.
The mathematical formulation of the problem is as follows:

$$ \frac{\partial i(x,t)}{\partial x} = - i(x,t) c(x,t) $$
$$ \frac{\partial c(x,t)}{\partial t} = -i(x,t) c(x,t) $$
s.t. $c(z,0) =1$, $\frac{\partial c}{\partial x} = 0$ at $x=0$ and $x=L$, with $L$ e.g. 15, and $i(0,t) = 1$.

The implementation should IMO be something like below, but I'm experiencing problems with the variational formulation (for starters):

from dolfin import *

T = 5.0            # final time
num_steps = 500    # number of time steps
dt = T / num_steps # time step size

# Create mesh
mesh = RectangleMesh(Point(0.0, 0.0), Point(15.0, 5.0), 30, 10)

# Define function space
P1 = FiniteElement('P', triangle, 1)
element = MixedElement([P1, P1])
V = FunctionSpace(mesh, element)

# Define test functions
v_1, v_2= TestFunctions(V)

# Define functions
u = Function(V)
u_n = Function(V)

# Split system functions to access components
u_1, u_2 = split(u)
u_n1, u_n2 = split(u_n)

# Define Dirichlet boundary (x = 0)
def boundary(x):
    return x[0] < DOLFIN_EPS

# Define boundary condition
u_01 = Constant(1.0)
bc = DirichletBC(V.sub(0), u_01, boundary)

# Define expressions used in variational forms
k = Constant(dt)

# Define variational problem
F = div(u_1)*v_1*dx + u_1*u_2*v_1*dx  \
  + ((u_2 - u_n2) / k) * u_1*u_2*v_2*dx

# Time-stepping
t = 0
for n in range(num_steps):

    # Update current time
    t += dt

    # Solve variational problem for time step
    solve(F == 0, u, bc)

    # Update previous solution
    u_n.assign(u)

To narrow down the initial problem, how would the function

$$ \frac{\partial i(x,t)}{\partial x} = - i(x,t) c(x,t) $$

be written in 'FEniCS python? Considering that

div(u_1)*v_1*dx + u_1*u_2*v_1*dx

gives errors w.r.t. shape differences.

Answer 1

You cannot take the divergence of a scalar. The derivative with respect to the first spatial coordinate is

 u_1.dx(0)

Comments

And if we were to write the mathematical form in a non-1D form?
So

$$ \nabla i = -i(\boldsymbol{x},t) c(\boldsymbol{x},t) $$

---

This does not make sense. The left hand side is a vector, while the right hand side is a scalar.

---

My apologies, you are completely right,

This should have been including a $\cdot$, as

$$ \nabla \cdot i(\boldsymbol{x},t) = -i(\boldsymbol{x},t) c(\boldsymbol{x},t), $$

which (please correct me if I'm wrong) in for example two-dimensions would be equivalent to:

$$ \frac{\partial i (\boldsymbol{x},t)}{\partial x} + \frac{\partial i (\boldsymbol{x},t)}{\partial y} = -i(\boldsymbol{x},t) c(\boldsymbol{x},t). $$

So a scalar equals a scalar. Since, as far as I'm aware, this is the meaning of the div operator, see e.g. WolframAlpha. With your original answer in mind, is this notation different in FEniCS?

---

The meaning of the div operator is divergence. This is a vector operation, which produces a scalar field from a vector field. This is clear from your WolframAlpha link, see also https://en.wikipedia.org/wiki/Divergence

---

Thanks for your quick response, I think I (finally) understood the trivial thing you're saying - am I correct in stating the equation I put above would thus be:

 (u_1.dx(0)+u_1.dx(1))*v_1*dx + u_1*u_2*v_1*dx = 0

since we're indeed dealing with a scalar field?

With this the code now becomes as shown below, which unfortunately does not converge. Any ideas on how to overcome this?

from dolfin import *

T = 5.0            # final time
num_steps = 500    # number of time steps
dt = T / num_steps # time step size

# Create mesh
mesh = RectangleMesh(Point(0.0, 0.0), Point(15.0, 5.0), 90, 30)

# Define function space for system of concentrations
P1 = FiniteElement('P', triangle, 1)
element = MixedElement([P1, P1])
V = FunctionSpace(mesh, element)

# Define test functions
v_1, v_2= TestFunctions(V)

# Define functions for velocity and concentrations
u = Function(V)
u_n = Function(V)

# Split system functions to access components
u_1, u_2 = split(u)
u_n1, u_n2 = split(u_n)

# Define Dirichlet boundary (x = 0)
def boundary(x):
    return x[0] < DOLFIN_EPS

# Define boundary condition
u_01 = Constant(1.0)
bc = DirichletBC(V.sub(0), u_01, boundary)

# Define expressions used in variational forms
k = Constant(dt)

# Define variational problem
F = (u_1.dx(0)+u_1.dx(1))*v_1*dx + u_1*u_2*v_1*dx  \
  + ((u_2 - u_n2) / k) * u_1*u_2*v_2*dx

# Create VTK files for visualization output
vtkfile_u_1 = File('u_1.pvd')
vtkfile_u_2 = File('u_2.pvd')

# Create progress bar
progress = Progress('Time-stepping')
set_log_level(PROGRESS)

# Time-stepping
t = 0
for n in range(num_steps):

    # Update current time
    t += dt

    # Solve variational problem for time step
    solve(F == 0, u, bc)

    # Save solution to file (VTK)
    _u_1, _u_2 = u.split()
    vtkfile_u_1 << (_u_1, t)
    vtkfile_u_2 << (_u_2, t)

    # Update previous solution
    u_n.assign(u)

    # Update progress bar
    progress.update(t / T)

# Hold plot
interactive()

---

The problem does not have an elliptic component, i.e.

u_1.dx(0)*v_1.dx(0)*dx

so you have to add some sort of stabilization. The choice of stabilization is outside the scope of this board, but

 v_1 += CellSize(mesh)*v_1.dx(0)

might work.

---

Thanks again!

The suggestion doesn't appear to be a direct plug&play solution but I'll look into that direction further.