Reaction Diffusion - Mixed Element time loop issues

variational formulation
1

I am trying to solve a rxn-diffusion problem (Diffusive Fischer Kolmogorov eqns) in the form:

\frac{\partial{b}}{\partial{t}} = D_b \nabla^2 b + bn \\ \frac{\partial{n}}{\partial{t}} = D_n \nabla^2 n - \eta \ bn \\

the solution seems to stay the same throughout the loop, so I am worried I am missing something on that.

I primarily used the tutorial examples for the heat equation and the advection - reaction system. I would also be interested in using the Crank-Nicolson method for time-stepping, if that would be better for the problem.

Here is my minimal example code:

T = 1 # final time
num_steps = 100 # number of time steps
dt = T / num_steps # time step size
eta = 0.01 # diffusion coefficient
D_b = 0.01 # diffusion coeff
D_n = 0.002 # diffusion coeff

# Create mesh
nx = ny = 50
square_side = 10

mesh = RectangleMesh(Point(-square_side/2, -square_side/2), Point(square_side/2, square_side/2), nx, ny, "crossed")

# Define function space
P1 = FiniteElement('P', triangle, 1)
element = MixedElement([P1,P1]) # b, n

# Define function space
V = FunctionSpace(mesh, element)

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

# Define functions for scalar fields
# and split to access components
u = Function(V)
u_1, u_2 = split(u) # u_1 == b, u_2 == n

u_n = Function(V)
u_n1, u_n2 = split(u_n)

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

# Define variational problem

# # Define boundary condition
# assume homogeneous Neumann boundary 
# conditions on the entire boundary for u1, u2
# that is, \partial u_i/ \partial n = 0

# # Define initial value

u_0 = Expression(('exp(-a*(pow(x[0], 2) + pow(x[1], 2)))', 'exp(-b*(pow(x[0], 2) + pow(x[1], 2)))'), 
                 degree=1, a=1, b = 0.01)
u_n1, u_n2 = interpolate(u_0, V)


# Define variational problem
F_b = ((u_1-u_n1)/dt)*v_1*dx - D_b*dot(grad(u_1), grad(v_1))*dx - u_1*u_2*v_1*dx
F_n = ((u_2-u_n2)/dt)*v_2*dx - D_n*dot(grad(u_2), grad(v_2))*dx + eta*u_1*u_2*v_2*dx

F = F_b + F_n


# Time-stepping

t = 0
for n in range(num_steps):
    # Update current time
    t = t + dt
    # Compute solution
    solve(F == 0, u)
    # Update previous solution
    u_n.assign(u)
    u_n1, u_n2 = split(u_n)
2

Note that by doing this, you overwrite the connection between u and u_n1, u_n2. I propose to do:

u_n = interpolate(u_0, V)
u_n1, u_n2 = split(u_n)

and remove

Complete code:

from dolfin import *
T = 1 # final time
num_steps = 100 # number of time steps
dt = T / num_steps # time step size
eta = 0.01 # diffusion coefficient
D_b = 0.01 # diffusion coeff
D_n = 0.002 # diffusion coeff

# Create mesh
nx = ny = 50
square_side = 10

mesh = RectangleMesh(Point(-square_side/2, -square_side/2), Point(square_side/2, square_side/2), nx, ny, "crossed")

# Define function space
P1 = FiniteElement('P', triangle, 1)
element = MixedElement([P1,P1]) # b, n

# Define function space
V = FunctionSpace(mesh, element)

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

# Define functions for scalar fields
# and split to access components
u = Function(V)
u_1, u_2 = split(u) # u_1 == b, u_2 == n


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

# Define variational problem

# # Define boundary condition
# assume homogeneous Neumann boundary
# conditions on the entire boundary for u1, u2
# that is, \partial u_i/ \partial n = 0

# # Define initial value

u_0 = Expression(('exp(-a*(pow(x[0], 2) + pow(x[1], 2)))', 'exp(-b*(pow(x[0], 2) + pow(x[1], 2)))'),
                 degree=1, a=1, b = 0.01)

u_n = interpolate(u_0, V)
u_n1, u_n2 = split(u_n)


# Define variational problem
F_b = ((u_1-u_n1)/dt)*v_1*dx - D_b*dot(grad(u_1), grad(v_1))*dx - u_1*u_2*v_1*dx
F_n = ((u_2-u_n2)/dt)*v_2*dx - D_n*dot(grad(u_2), grad(v_2))*dx + eta*u_1*u_2*v_2*dx

F = F_b + F_n


# Time-stepping

t = 0
for n in range(num_steps):
    # Update current time
    t = t + dt
    # Compute solution
    solve(F == 0, u)
    # Update previous solution
    u_n.assign(u)
    print(norm(u_n))
1 Like
3

Thank you, that fixed the solution updating in the time loop (although I am getting errors in the Newton solver not converging)

Related question: is it possible to add a Heaviside (\Theta) on the functions u_1, u_2 in the variational form (or how would one implement):

\frac{\partial{b}}{\partial{t}} = D_b \nabla^2 b + bn\Theta(b - \beta)

(Maybe using the previous time step u_n1 ?)
(I am aware of the ? : syntax in the C++ Expressions, but not sure if it would apply on the actual solution function)

4

You can use a ufl conditional (loaded with from ufl import conditional, following the documentation:

Conditional Operators

Conditional

UFL has limited support for branching, but for some PDEs it is needed. The expression c in:

c = conditional(condition, true_value, false_value)

evaluates to true_value at run-time if condition evaluates to true, or to false_value otherwise.

This corresponds to the C++ syntax (condition ? true_value: false_value) , or the Python syntax (true_value if condition else false_value) .

1 Like