Fenics for parabolic PDE

Hi,

I am trying to solve the following PDE which shall result in the probability that a ball will not enter in a hole of radius 2 and centered in the plan.
To make it short, the PDE is

Pt + (U.Nabla)P-2Delta.P = 0

with:

u(x1,x2)= (1,1)+gamma*nabla(log(x1^4+x2^4))=0

I have done a code but I am blocked when I shall define the variation problem.
When I am adding gamma in f, it does not work and when I am trying to implement f in F, it does not work either.
Can someone help?

from fenics import *
import numpy as np

T = 10.0 # final time
num_steps = 100 # number of time steps
dt = T / num_steps # time step size
gamma = 1.5

# Create mesh and define function space
nx = ny = 10
mesh = UnitSquareMesh(nx, ny)
V = FunctionSpace(mesh, 'P', 1)

# Define boundary condition
P_D = Expression('x[0]*x[0] + x[1]*x[1] -4',degree=2, alpha=alpha, beta=beta, t=0)
def boundary(x, on_boundary):
return on_boundary
bc = DirichletBC(V,P_D,boundary)

# Define initial value
P_n = interpolate(P_D,V)

# Define variational problem
P = TrialFunction(V)
v = TestFunction(V)
#Starting position
X0 = Constant(3)
Y0 = Constant(0)
f = Expression(('1+4*pow(x[0],3)/((pow(x[0],4)+pow(x[1],4)))','1+4*pow(x[1],3)/((pow(x[0],4)+pow(x[1],4)))'), element = V.ufl_element())
F = inner(nabla_grad(P), nabla_grad(v))*dx - 2*(P_n + dt)*dx
a, L = lhs(F), rhs(F)
# Time-stepping
P = Function(V)
t = 0
for n in range(num_steps):
# Update current time
t += dt
P_D.t = t
# Compute solution
solve(a == L, P, bc)
# Plot solution

plot(u)
# Compute error at vertices
P_e = interpolate(P_D, V)
error = np.abs(np.array(P_e.vector()) - np.array(P.vector())).max()
print('t = %.2f: error = %.3g' % (t, error))
# Update previous solution
P_n.assign(P)


Could you clarify the problem you are trying to solve? For instance:

• What is the meaning of “Pt” and “Delta” in your PDE?
• What is “gamma”? (i.e. a scalar constant, a scalar function, a matrix constant)
• What domain is your problem defined on? Are you using a UnitSquareMesh for your example code only, or is [0,1] \times [0,1] the domain you wish to study?
• What are the boundary conditions for P?
• Is u defined correctly? (At the origin, \log{\left(x_1^4 + x_2^4\right)} = \log 0 is NaN, which will seriously affect your finite element solution.)

If you haven’t yet, please read through Read before posting: How do I get my question answered?

In particular, it helps if you format equations using dollar signs ("\$") just as in LaTeX:

produces

Additionally, you should surround your code with backticks () to make it easier to copy and paste. Single backticks, i.e. import fenics produce in-line code: import fenics. Triple backticks produce code blocks:

produces

from fenics import *

import numpy as np

1 Like

Indeed, by Pt, I meant the probability for the ball to touch/fall in the hole.

• Delta: see Latex corrected expression after.
• Gamma is just a scalar.
• G is the plan less the hole represented by the disc B centered.
• I put also the boundaries conditions and log cannot be equal to 0 because it is the center of the hole.
Do you see how I can improve this code?
Best,
P_t=(u.\nabla)P-2\Delta.P=0\quad in\quad G\times[0,T] \\ with \\ u(x_{1},x_{2})=\begin{pmatrix} 1 \\ 1 \end{pmatrix}+\gamma\nabla(log(x_{1}^4+x_{2}^4)=0 \\

And the last boundaries conditions are:

Hi @beamer,

Is there a typo in the following equation?

Comparing with your first post, it looks like you meant

P_t + (u.\nabla)P−2\Delta.P=0 \quad in \quad G×[0,T]

Is P_t the quantity you want to solve for? What is P? This looks like a time-dependent problem; is P_t the time derivative of P, i.e. P_t = \partial P / \partial t? Regarding the \Delta, does this refer to the Laplacian operator, \nabla^2?

It sounds like you will need to utilize a different mesh for this problem, as the UnitSquareMesh used in your code defines a square [0,1] \times [0,1], which includes the origin.

Yes, it is a typo, sorry
Indeed, it is a +.

Yes, Pt is the quantity I want to solve. It is a probability, the probability that the ball touch the hole.
Yes, Pt is time derivative.
Indeed, Delta refers to the Laplacien operator, exactly what you mentionned.

OK for the mesh, do you think it makes sense to have a mesh corresponding to the plan? I think we can restrict it to [-10,10] x [-10,10]

The mesh should correspond to whatever geometry you are simulating, with the disk B cut out. To generate the square mesh you mentioned, you could use the following in gmsh:

SetFactory("OpenCASCADE");
L = 20; // size of the box
R = 2;  // radius of the disk
Rectangle(1) = {-L/2, -L/2, 0, L, L, 0};
Disk(2) = {0, 0, 0, R, R};
BooleanDifference{ Surface{1}; Delete; }{ Surface{2}; Delete; }
MeshSize {6,7,8,9} = L/10;
MeshSize {5} = R/10;
Physical Surface(1) = {1};
Physical Curve(2) = {5};
Physical Curve(3) = {6,7,8,9};


with the command

gmsh mesh.geo -2 -format msh4 -o mesh.msh


to generate the mesh, and the following Python script to convert it to XDMF format:

import meshio

msh.prune_z_0()

triangle_data = msh.cell_data_dict["gmsh:physical"]["triangle"]
line_data = msh.cell_data_dict["gmsh:physical"]["line"]

triangle_cells = msh.cells_dict["triangle"]
line_cells = msh.cells_dict["line"]

triangle_mesh =meshio.Mesh(points=msh.points,
cells=[("triangle", triangle_cells)],
cell_data={"cell_marker":[triangle_data]})
line_mesh = meshio.Mesh(points=msh.points, cells={"line": line_cells},
cell_data={"line_marker":[line_data]})

meshio.write("mesh.xdmf", triangle_mesh)
meshio.write("bnd_markers.xdmf", line_mesh)


Your PDE can be solved, starting from some initial condition P(x,0), using the following code. In the code below, you will see that I have considered a homogeneous condition on the external boundary of G, as well as three different inhomogeneous initial conditions. The solution for all three of these initial conditions eventually settles to the same steady-state, namely the solution to

(u \cdot \nabla)P - 2 \Delta P = 0

Based on this solution, it does not appear that [-10, 10] \times [-10, 10] is a large enough domain to accurately represent the entire x_1 x_2-plane, so you may wish to consider experimenting with the mesh. In addition, this solution does not satisfy the final condition P(x,T)=0 that you have specified, so you should double-check the formulation of the problem.

from fenics import *
import numpy as np

T = 1000.0 # final time
N = 100
dt = T/N
Dt = Constant(dt) # timestep

# Create mesh and define function space
mesh = Mesh()
with XDMFFile("mesh.xdmf") as xdmf:
V = FunctionSpace(mesh, 'P', 2)

# Load boundary markers
mvc = MeshValueCollection('size_t', mesh, mesh.topology().dim() - 1)
with XDMFFile("bnd_markers.xdmf") as xdmf:
mf = MeshFunction('size_t', mesh, mvc)

# Define boundary condition
P_disc  = Constant(1.0)
P_outer = Constant(0.0)
bc_disc  = DirichletBC(V, P_disc, mf, 2)
bc_outer = DirichletBC(V, P_outer, mf, 3)
bcs = [bc_disc, bc_outer]

# Define initial value
# P_old = interpolate(Expression('0.', degree=2),V)
# P_old = interpolate(Expression('1.', degree=2),V)
P_old = interpolate(Expression('x[0]*x[0] + x[1]*x[1] - 4', degree=2),V)

#Define variational problem
u = TrialFunction(V)
v = TestFunction(V)

gamma = Constant(1.5)
f1 = Expression(('1.','1.'), degree=0)
f2 = Expression(('4*pow(x[0],3)/(pow(x[0],4)+pow(x[1],4)+DOLFIN_EPS)',
'4*pow(x[1],3)/(pow(x[0],4)+pow(x[1],4)+DOLFIN_EPS)'),
degree=4)
f = f1 + gamma*f2
F = (u*v + Dt*inner(f, grad(u))*v + 2*Dt*inner(grad(u), grad(v)) - P_old*v)*dx
a, L = lhs(F), rhs(F)

#Time-stepping
P = Function(V)
P.rename("P", "P")
P.assign(P_old)

t = 0.0
xdmf = XDMFFile("result_Pparab_initial.xdmf")
# xdmf.write(P, t)

while t < T:
#Update current time
t += dt

#Compute solution
solve(a == L, P, bcs)

#Plot solution
xdmf.write(P, t)

#Update previous solution
P_old.assign(P)
xdmf.close()


It’s looking really promising.
Thank you very much!
Seems that I have an issue in the part SetFactory.

For the last command line (as from BooleanDifference), the following error occurs (under the first “{”):

File ipython-input-16-7a460453505f, line 28
BooleanDifference{ Surface{1}; Delete; }{ Surface{2}; Delete; }
^
SyntaxError: invalid syntax

I tried to find documentation in internet, read the site of Bertrand Thierry (Geometries and Operations) and try alternatives for this one (and the 5 following lines) but it does not work.

Do you have an idea from where it may come?

Another question concerning the probability P to fall in the hole (since I could not test it entirely).
Say that the ball is located (3,3), at time 0, how can you obtain P((x,y),t)?

Did you run the mesh generation command from within iPython? You should copy the following lines to a file called mesh.geo:

then run the following command from the shell (e.g. bash):

Regarding your last question, you could for example obtain P((x_0,y_0), t_0) by inserting the following lines below the solve command in the while loop:

    if near(t, t0, eps=1e-8):
print(P(x0,y0))


assuming that you have defined x0, y0, and t0 previously in the code.

Thank you, indeed, it works.
Just last precision.
You write P(x0,y0) but don’t you think we may use an interpolation here at this stage (because of the fact that we have a mesh).
Would you agree, what tool (Fenics) would you recommend?

Calling P(x0, y0) interpolates the solution at the point (x_0, y_0).

(In particular, it invokes the __call__ method of Function, which in turn invokes the eval method on the associated C++ object. As you can see, the C++ object uses the element basis functions and coefficients to compute the function value at the requested point.)

1 Like

Thank you, 100% clear!

I have here a general question concerning the mesh.
The highest is the mesh (I mean for value which are unlikely to be reached), the less accurate is the solution of the PDE.
I suspect I have to compensate with the number of step but I am not 100% sure here.
What do you think?

What do you mean by the highest mesh? If you want to improve the solution accuracy, you can refine the mesh by changing the following lines in your .geo file:

You can change L/10 to a smaller value, e.g. L/20, L/30`, etc to generate a finer mesh.

Indeed, thanks for your help!