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
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) 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.1*(1+tanh(x))'), 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)