Why am I getting only the trivial solution to pair of nonlinear PDEs?

variational formulation
1

I am absolutely new to FEniCS (and python) - running via Anaconda on a mac.
In an attempt to practice what I’ve learned form the FEniCS tutorial, I am trying to solve the following pair of PDEs in 2D (one is non-linear the other is linear)

du/dt = D∆u + u(1-u)(u-a) -v
dv/dt = ∆v + bu - v
(0<a<1, b>0 are parameters)

using the variational form. I wrote to code by emulating the code from the FEniCS tutorial
( ft09_reaction_system.py)
but I keep on getting the trivial solution for both u and v

I think the problem is with the ICs (or their interpolation). I want to choose ICs that are equal along the y (x[1]) axis of the mess, because I know the solution in 1D is a pair of traveling waves (it’s a simple problem, the equation for u generates the wave, and v is a linear filter on u)

but because I am really new it may be something else.

Would love it if someone could help me figure out where I went wrong?

here’s the code:

# %load AIRDpde.py


# I am trying to solve the following pair of PDEs (one is non-linear the other is linear) 

# du/dt = D*∆u + u*(1-u)*(u-a) -v
# dv/dt =   ∆v + b*u - v
# (0<a<1, b>0 are parameters)

# using the variational form, but I keep on getting the trivial solution for both u and v
# I think the problem is with the ICs


from fenics import *

# parameters
T = 5.0            # final time
num_steps = 50     # number of time steps
dt = T / num_steps # time step size
D = 0.016          # U diffusion coefficient
a = 0.25           # model param - x axis
b = 0.14           # model param - y axis


# Create mesh 
nx = ny = 30
mesh = RectangleMesh(Point(-3, -3), Point(3, 3), nx, ny)

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

# Define boundary condition
# none - Neumann bcs

# Define initial value
u_0 = Expression(('0.6*(1+tanh(beta*(x[0])))', '0.1*(1+tanh(x[0]))'),
                    degree=1,beta=2)
u_n = interpolate(u_0, V)
                 
# Define variational problem

u = Function(V) # current step
u_n = Function(V) # previous step
u_1, u_2 = split(u)
u_n1, u_n2 = split(u_n)

v_1, v_2 = TestFunctions(V)


F = ((u_1 - u_n1)/dt)*v_1*dx + D*dot(grad(u_1), grad(v_1))*dx + u_2*v_1*dx - u_1*(1.0-u_1)*(u_1-a)*v_1*dx \
  + ((u_2 - u_n2)/dt)*v_2*dx +   dot(grad(u_2), grad(v_2))*dx + u_2*v_2*dx - b*u_1*v_2*dx 

# create vtk files
vtkfile_u_1 = File('AIRD/u_1.pvd')
vtkfile_u_2 = File('AIRD/u_2.pvd')

# 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)
    
    # 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)
2

You are redefining u_n, so as a fist step, i would suggest removing the last line of the code i quoted above, such that you only have

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

Thank you. That did it. I was wondering about that re-definition of u_n that was not done in the previous demos. Thanks for explaining that. I get it, once I redefine it, it’s value is 0 and of course that is the unstable f.p. of the dynamics. Thanks again