I’m looking for some advice on how easy it is to implement FeniCS to solve a space evolution Fokker-Planck PDE in FeniCS. My problem is
\frac{\partial P}{\partial L}(L,\eta) = \frac{1}{\theta}\bigg(2\eta\frac{\partial P}{\partial \eta} +(\eta^2+1)\frac{\partial^2 P}{\partial \eta^2}\bigg), \quad \eta >1,
with initial condition P(L=0) = \delta(\eta-1) and \theta is a positive constant. The Dirichlet boundary conditions can be set later. I choose to discretise the time derivative via
\bigg(\frac{\partial P}{\partial L}\bigg)^{n+1} \approx \frac{P^{n+1}-P^n}{\Delta L}
so
P^{n+1}-\frac{\Delta L}{\theta}\bigg(2\eta\frac{\partial P^{n+1}}{\partial\eta}+(\eta^2+1)\frac{\partial^2 P^{n+1}}{\partial \eta^2}\bigg) - P^n = 0.
The weak form of the problem I want to write like
a(P,v) = L_{n+1}(v)
where v is a test function. I think then I can write
\begin{align}
a(P,v) &= \int_{\Omega}Pv(\eta)-\frac{\Delta L}{\theta}\bigg(2\eta\frac{\partial P}{\partial\eta}-\bigg((\eta^2+1)\frac{\partial v}{\partial\eta}+2\eta\bigg)\frac{\partial P}{\partial\eta}\bigg)v(\eta)\,d\eta, \\
L_{n+1}(v) &= \int_{\Omega}P^{n}v(\eta)\,d\eta.
\end{align}
I’m not sure whether this is correct, and if I should do something with the second derivative in a(P,v). Any help would be appreciated to help implement this in fenics.
I have the analytical solution to this problem, in the form of an integral solution. Please let me know if I can improve my question or anything. I am really passionate about using FeniCS to solve my problem.
Hi,
not an expert on this equation but I think you are missing a v in front of the parenthesis term of a. Also I would integrate by part the second derivative to obtain the weak form with only first order derivatives.
1 Like
Hi, thanks for your comment. I have removed the second derivative and shifted it onto my test function v(\eta), I now have a(P,v) = \int_{\Omega}Pv(\eta)-\frac{\Delta L}{\theta}\bigg(2\eta\frac{\partial P}{\partial\eta}-\bigg((\eta^2+1)\frac{\partial v}{\partial\eta}+2\eta\bigg)\frac{\partial P}{\partial\eta}\bigg)v(\eta)\,d\eta. Does this look correct? If so, do I need to now define more conditions to implement into FeniCS? Thank you.
I think that you should rather have a(P,v) = \int_\Omega \left(P v + \frac{\Delta L}{\theta}(\eta^2+1)\dfrac{\partial P}{\partial \eta} \dfrac{\partial v}{\partial \eta}\right) d\eta since I did not realize at first that the right hand side of the initial PDE can in fact be rewritten as \dfrac{1}{\theta}\dfrac{\partial}{\partial \eta}\left((\eta^2+1)\dfrac{\partial P}{\partial \eta} \right)
Also there should be boundary terms but you did not specify your boundary conditions.
1 Like
Yes absolutely, you can write it like that.
I have managed to solve this PDE analytically with only the initial condition using a Mehler-Fock transform. I can only see one boundary condition which is perhaps obvious, \lim_{\eta\rightarrow\infty}P(L,\eta)\rightarrow 0. Since the solution P(L,\eta) is also a probability density function the solution must obey \int_{1}^{\infty}P(L,\eta)d\eta = 1.
Sorry for the perhaps silly question, but would FeniCS require another condition? Thank you for your continued help!
I think that you should have some BC at \eta=1, maybe it’s pure Neumann \partial P/\partial \eta=0 ? Also does the normalization constraint must be added as an additional condition or is it guaranteed by the PDE if you start from an initial condition which satisfies it ? If not, you will have to add it through a Lagrange multiplier, see for instance here https://fenicsproject.org/docs/dolfin/2018.1.0/python/demos/neumann-poisson/demo_neumann-poisson.py.html
1 Like
Hi, I believe prescribing the initial condition preserves the probability integrating to 1 over the domain. I have tried to implement the problem in FeniCS. I am having some trouble understanding how to implement the initial dirac delta and also the boundary conditions I am hoping you can help with. Thank for for your help up until this point. My code:
from __future__ import print_function
from fenics import *
import numpy as np
import matplotlib.pyplot as plt
L = 2.0 # final L
num_steps = 10 # number of length steps
dl = L / num_steps # length step size
theta = 1 #Lloc constant term
class Delta(Expression):
def __init__(self, eps):
self.eps = eps
def eval(self, values, x):
eps = self.eps
values[0] = eps/pi/(x[0]**2 + eps**2)
delta = Delta(x[0]-1)
# Create mesh and define function space
nx = ny = 64
mesh = UnitSquareMesh(nx, ny)
V = FunctionSpace(mesh, 'P', 1)
# Define boundary condition
u_D = delta(x[0]-1) #?
def boundary(x, on_boundary):
return on_boundary
bc = #u = 0 at x[0] = 1?
# Define initial value
u_n = interpolate(u_D, V)
#u_n = project(u_D, V)
# Define variational problem
u = TrialFunction(V)
v = TestFunction(V)
a = (u*v + (dl/theta)*(x[0]*x[0]+1)*derivative(u,x)*derivative(v,x))*dx
L=u*v*dx
# Time-stepping
u = Function(V)
l = 0
for n in range(num_steps):
# Update current time
l += dl
u_D.l = l
# Compute solution
solve(a == L, u, bc)
# Plot solution
plot(u)
# Update previous solution
u_n.assign(u)
# Hold plot
plt.show()
Hi consider this, I switched to 1d just for easier debugging. I used a regular “Expression” and normalized your “Dirac” function with respect to its integral. However, with such BCs the normalization is not conserved. But I don’t see how it could since P=0 is clearly a solution to the equation when time goes to infinity.
from __future__ import print_function
from fenics import *
import numpy as np
import matplotlib.pyplot as plt
L = 2.0 # final L
num_steps = 10 # number of length steps
dl = L / num_steps # length step size
theta = 1 #Lloc constant term
# Create mesh and define function space
nx = 100
mesh = IntervalMesh(nx, 1., 10.)
V = FunctionSpace(mesh, 'P', 1)
x = SpatialCoordinate(mesh)
# Define boundary condition
u_D = Expression("1./(pow(x[0]-y,2) + pow(eps, 2))", eps=1e-6, pi=pi, y=1., degree=2)
# plot(u_D)
def right(x, on_boundary):
return near(x[0], 10.)
def left(x, on_boundary):
return near(x[0], 1.)
bc = DirichletBC(V, Constant(0), right)
# Define initial value
u_n = interpolate(u_D, V)
C = assemble(u_n*dx)
u_n.assign(u_n/C)
# Define variational problem
u = TrialFunction(V)
v = TestFunction(V)
a_form = (u*v + (dl/theta)*(x[0]*x[0]+1)*u.dx(0)*v.dx(0))*dx
L_form = u_n*v*dx
# Time-stepping
u = Function(V)
l = 0
for n in range(num_steps):
# Update current time
l += dl
# Compute solution
solve(a_form == L_form, u, bc)
# Plot solution
plot(u)
# Update previous solution
u_n.assign(u)
print(assemble(u*dx))
# Hold plot
plt.show()
1 Like
Hi. Thanks a lot for this. I’ve ran the code and the results look sensible from what I expected. I have a few questions about your code. Can you explain what you mean by normalizing wrt my integral? Is this essentially a point source at x=1? Also I am confused why you have defined pi and not used it.I am also confused why right and left functions are defined? My last question, is the assemble steps essentially the same as using project()? Thank you so much for your help, it is invaluable.
Hi,
yes sorry pi was a leftover which should not be there, I defined left in case you wanted to add a specific BC but if not, you don’t need it indeed.
For the normalization, you cannot have a true Dirac mass as an initial condition and even a continuous regularization of a Dirac mass will not be properly resolved by the FE mesh if it is too coarse. In your initial implementation the integral of your initial condition was not equal to 1 due to discretization errors. Therefore, instead of normalizing with a constant obtained when computing the exact integral, I normalize with the constant computed using the FE discretized integral.
No, assemble just computes the integral. project is solving an L2 minimization problem.
1 Like
Hi Jeremy, thanks for the clarification that helps! I understand that you can’t have a true Dirac mass as initial condition. I still can’t understand what you mean when you say you normalize with the constant computed using the finite element discretized integral - could you be more explicit here?
Also my last question is, when the code prints “print(assemble(u*dx))” at each length step, is that the area under the solution? I.e. \int_{1}^{\eta_\inf}P(L,\eta)d\eta? I thought that it would be equal to 1 - or am I misunderstanding? Thank you and sorry for my persistent questions.