Hey guys,

I’m pretty new to FEniCS and don’t have a lot of knowledge about PDEs but from all I read so far FEniCS should be able to handle the kind of PDE I’m working with. My code is somewhat similar to the one used in the tutorial on PDE systems but it doesn’t even kind of reproduce what I’m expecting. The concentrations both tend to jump to 0 after being initiated at a different value. I was digging through some questions but I suppose the kind of code I need really depends on the PDE I’m using. Does anyone of you have an Idea how to enhance or change the code to construct a good solution? Until maybe about a 1second it once worked but then I noticed I had something wrong and anyway my goal would be T = 10.000 or ideally T = 500.000 - 1.000.000. Any Ideas? I’d really appreciate any help.

The setup:

1D two component system

\frac{dm}{dt} = \frac{d}{dx}(D_m(x)*\frac{dm}{dx}) + f(m,c)

\frac{dm}{dt} = \frac{d}{dx}(D_c(x)*\frac{dc}{dx}) - f(m,c)

f(m,c) = c(m+1) + 2*m/(m+1)

D_m(x) = 0.1

D_c(x) = 1.5 + arctan(50*(x-L/2))/\pi

with initial condition c(0,x) = 0 and m(0,x) = 22/9+0.2*cos(\pi*x/L)

and von Neumann boundary conditions with derivatives 0 at 0 and L

My code:

```
from fenics import *
import numpy as np
T = 2
steps = 1000
frames = 30
stamps = T*frames
stampspace = int(steps/stamps)
dt = T/steps
l=50
nx=1000
mesh = IntervalMesh(nx,0,l)
P1 = FiniteElement('P', interval, 2)
element = MixedElement([P1,P1])
V = FunctionSpace(mesh,element)
#Creating the functions
v1, v2 = TestFunctions(V)
u = Function(V)
#u_n = Function(V)
u1,u2 = split(u)
#u_n1,u_n2 = split(u_n)
fDm = Constant(0.1)
fDc = Expression("1.5+atan(50*(x[0]-l/2))/pi", degree=2, l=l)
u_0 = Expression(("22/9+0.2*cos(pi*x[0]/l)","0"), degree = 2, l=l)
u_n = interpolate(u_0, V)
u_n1,u_n2 = split(u_n)
k = Constant(dt)
F = ((u1-u_n1)/k)*v1*dx+fDm*dot(grad(u1),grad(v1))*dx-(u2*(u1+1)-2*u1/(u1+1))*v1*dx+((u2-u_n2)/k)*v2*dx+fDc*dot(grad(u2),grad(v2))*dx+(u2*(u1+1)-2*u1/(u1+1))*v2*dx
vtkfile_u1 = File('B1/m.pvd')
vtkfile_u2 = File('B1/c.pvd')
t = 0
u1_,u2_ = u_n.split()
vtkfile_u1 << (u1_,t)
vtkfile_u2 << (u2_,t)
for n in range(steps):
t+=dt
solve(F==0, u)
if n%stampspace==0:
u1_, u2_ = u.split()
vtkfile_u1 << (u1_,t)
vtkfile_u2 << (u2_,t)
print(t)
u_n.assign(u)
´´´
Thanks in advance!
Sarah
![Code_B1|690x468](upload://5CieWlAuY7Hlry35tRj1vDXQ1Qj.png)
```