Solution to PDE is rubbish, shouldn't it be easy?

mathematics
1

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

Hi Sarah,

Welcome to FEniCS. I’m also pretty new here. :slightly_smiling_face:

I encourage you to have a look at Read before posting: How do I get my question answered?

If you surround your math with dollar signs ($), it will be formatted as math. E.g.

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

Surrounding your code with three backticks (`) causes it to be formatted as code. Starting a code block with three backticks plus “python” will add Python syntax highlighting, i.e. typing

will be formatted as:

from fenics import *

This makes your question easier to read, and easier to debug because we can copy-paste your code. :slightly_smiling_face:

2 Likes
3

Thanks for the quick intro, it should be a lot better readable now.

4

Hi Sarah,

Can you give more information about the solution you’re expecting?

With only one small change to your code, I get a solution that doesn’t “jump to 0”. The solution settles to a state shown below at t=3 and and remains there until at least t=50:

mc_t003mc_t050

The change I made to your code is in the for loop: I moved the u_n.assign(u) command above the if statement and write the results from u_n instead of from u. In the code you posted, you write the first time step from u_n and all subsequent time steps from u. Since these are different FEniCS objects, they are written to different variable names in the VTK files. Perhaps this is why the solution appears to jump to 0 after the first time step? My code is below:

from fenics import *
import numpy as np

T = 50
steps = 100
frames = 2
stamps = T*frames
stampspace = int(steps/stamps)
dt = T/steps
l=50
nx=50


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)
		
    print(t)
        
    u_n.assign(u)
    
    if n%stampspace==0:
        u1_, u2_ = u_n.split()
        vtkfile_u1 << (u1_,t)
        vtkfile_u2 << (u2_,t)

As an aside, you may be able to save some computation time and still be able to obtain reasonable accuracy using a much coarser mesh and timestep than you are using. Note that I used nx=50 and steps=100 to create the plots above.

2 Likes
5

Hi conpierce,

first of all thanks a lot for the help, I tried out different sets of parameters this morning and rceived far more reasonable results. Simply switching the if statement and the u_n assignment as well as changing the variable I create the files from fixed the code. I didn’t know FEniCS saves the variable name as an essential feature.

When simulating the system with mathematica I obtain a solution in which hills and valleys form on the right side of the domain which then merge into a similar pattern as here for T=10000. But although that looks very reasonable since the right side of the domain is most likely linearly unstable I came here to verify the solution with smaller timesteps. Any Ideas how to verify or increase the certainty that these results are correct are still welcome but what I will do for now is I will solve some more familiar problems and compare those to the mathematica solution. If everything turns out to reproduce reasonable solutions I have a result, otherwise I might come back :slight_smile:

It’s already such a relief to see the code doing what it is supposed to do. Thank you!

Sarah