Hi all. I’ve been trying to successfully solve my PDE (Fokker-Planck equation)
\begin{aligned} \frac{\partial P}{\partial L}(L, x, y, z) &=\frac{A_{3}^{2} z^{2}}{2} \frac{\partial^{2} P}{\partial x^{2}}(L, x, y, z) \\ &+A_{3}^{2} x z \frac{\partial^{2} P}{\partial x \partial z}(L, x, y, z)+\left(A_{1}^{2}+A_{2}^{2}+A_{3}^{2} \frac{z^{2}}{x^{2}}\right) \frac{\partial^{2} P}{\partial y^{2}}(L, x, y, z) \\ &+A_{3}^{2} x^{2} \frac{\partial^{2} P}{\partial z^{2}}(L, x, y, z) \\ &+\left(\frac{A_{3}^{2}}{2}-A_{4}\right) x \frac{\partial P}{\partial x}(L, x, y, z) \\ &+\left(\frac{A_{3}^{2}}{2}+A_{4}\right) z \frac{\partial P}{\partial z}(L, x, y, z) \\ &-A_{5} \frac{\partial P}{\partial y}(L, x, y, z) \end{aligned}
with initial condition P(L=0,x,y,z) = \delta(x)\delta(y)\delta(z-1) of three Dirac masses and natural Neumann boundary conditions. The A_{1,2,3,4,5} terms are positive constants. Also z \in [1,\infty) and x,y\in [0,\infty). For the purposes of my implementation in FeniCS I will choose an upper bound of \sim150.
Since P(L,x,y,z) is a probability density function it must integrate to one over the domain so \int_{\Omega}P(L,\mathbf{x})\,d\mathbf{x} = 1. I check this at each length step in the code. I should say here also that in the code, at each length step I compute the quantity \mathbf{\tau} = \int_{1}^{150}\int_{0}^{150}\int_{0}^{150}\frac{P}{z^2}\,dx\,dy\,dz which is the reason for solving the PDE.
My implementation:
I start by approximating the length integral via
\begin{align} \bigg(\frac{\partial P}{\partial L}\bigg)^{n+1} \approx \frac{P^{n+1}-P^n}{\Delta L}, \end{align}
and writing the weak form of the PDE in bilinear form:
\begin{align}\nonumber
a(P,v) &= \int_{\Omega}\bigg[Pv + \Delta L\bigg(\frac{A_3^2z^2}{2}\frac{\partial v}{\partial x}\frac{\partial P}{\partial x} + \bigg(A_1^2 + A_2^2 + A_3^2\frac{z^2}{x^2}\bigg)\frac{\partial v}{\partial y}\frac{\partial P}{\partial y} + A_3^2x^2\frac{\partial v}{\partial x}\frac{\partial P}{\partial z} \\ \nonumber
&+ A_3^2\bigg(zv + xz\frac{\partial v}{\partial x}\bigg)\frac{\partial P}{\partial z} - \bigg(\frac{A_3^2}{2} - A_4\bigg)x\frac{\partial P}{\partial x}v, \nonumber
&-\bigg(\frac{A_3^2}{2} + A_4\bigg)z\frac{\partial P}{\partial z}v + A_5\frac{\partial P}{\partial y}v\bigg)\bigg]\,d\mathbf x, \\
L_{n+1}(v) &= \int_{\Omega}P^{n}v\,d\mathbf x,&
\end{align}
where \mathbf{x} = [x,y,z]^{T}. I approximate the Dirac masses via the limit definition (34).
Mesh:
I pick a 3D boxmesh.
My implementation in FeniCS:
from __future__ import print_function
from fenics import *
import numpy as np
import matplotlib.pyplot as plt
L = 10 #Final L
num_steps = 50 #Number of length steps
dl = L / num_steps #Length step size
#Define positive constants (chosen randomly here to be small and positive)
A_1 = 0.001
A_2 =0.001
A_3 = 0.001
A_4 = 0.009
A_5 = 0.007
#Intial mesh values
x0 = 0
y0 = 0
z0 = 1
#End mesh values
x1 = 150
y1 = 150
z1 = 150
#Number of cells in each direction of the mesh
nx = 10
ny = 10
nz = 10
#Define mesh
mesh = BoxMesh(Point(x0, y0, z0),Point( x1, y1, z1), nx, ny, nz)
V = FunctionSpace(mesh, 'P', 1)
x = SpatialCoordinate(mesh)
mpx=np.array(mesh.coordinates())
dx = Measure('dx',mesh)
#Dirac Delta (mass) initial condition with epsilon definition
u_D = Expression("(1./(pow(x[0],2) + pow(eps, 2)))*(1./(pow(x[1],2) + pow(eps, 2)))*(1./(pow(x[2]-1.,2) + pow(eps, 2)))", eps=1e-6, degree=2)
#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)
#Bilinear form
a_form = (u*v + dl*((A_3**2/2)*x[2]**2*v.dx(0)*u.dx(0)
+ (A_1**2 + A_2**2 + A_3**2*(x[2]**2/x[0]**2))*v.dx(1)*u.dx(1)
+ A_3**2*x[0]**2*v.dx(0)*u.dx(2)
+ A_3**2*(x[2]*v +x[0]*x[2]*v.dx(0))*u.dx(2)
- (A_3**2/2 + A_4)*x[2]*u.dx(2)*v
+ A_5*u.dx(1)*v
))*dx
L_form = u_n*v*dx
#Starting to build length steps
u = Function(V)
l = 0
#Create lists
mean_t = []
L_vals = []
for n in range(num_steps):
#Update current lengthj
l += dl
#Solve
solve(a_form == L_form, u)
u_file = File('u.pvd')
u_file << u, l
#Store lengths
L_vals.append(l)
#Update previous solution
u_n.assign(u)
#Print probability volume at each length step (should be preserved =1)
print("Area under PDF is:")
print(assemble(u*dx))
#Compute transmission coefficient at each length step
print("Transmission coefficient is:")
print(assemble((u/x[2]**2)*dx))
#Store power transmission coefficient in list at each l step
mean_t.append(assemble((u/x[2]**2)*dx))
#Plot mean power transmission versus length
plt.xlabel('$L$')
plt.ylabel('$\\mathbb{E}[\\tau(L)]$')
plt.scatter(L_vals,mean_t)
plt.show()
Problems:
When I solve this problem as the code stands above I get NaN solutions for the PDE. Could it be too complex or have I made a mistake with my weak formulation and I should recheck again…?
When I set A_{1,2,3} = 0 and run the code I get an answer somewhat believable. The volume integral \int_{\Omega}P(L,\mathbf{x})\,d\mathbf{x} = 1 computed at each L via assemble(u*dx) holds. The solution looks like this
in paraview.
Help?
Thanks for reading thus far. If anyone can comment on what I have done wrong in my implementation that would be extremely helpful as I am a novice with this tool. I would also appreciate if someone could tell me if this PDE is beyond the scope of being able to solve accurately with a FEM package. Thanks.