Hi all, I have been working on problems which require solving Fokker-Planck equations. I have been unable to successfully use FeniCS to obtain a convergent solution for my problem and I am hoping someone can add something with my updated code, since I have posted on this topic before. I am studying a PDE which can be written in a flux-gradient form with x,y\in[0,\infty) and z\in[1,\infty)
\begin{align}
\partial_L P(L,\mathbf{x}) = (\nabla\cdot \mathbf{a}\nabla + \mathbf{c}\cdot\nabla)P(L,\mathbf{x}).
\end{align}
where
\begin{align}
\mathbf{a} = \begin{bmatrix}
\frac{A_3^2z^2}{2} & 0 & \frac{A_3^2xz}{2} \\
0 &\frac12\bigg(A_1^2+A_2^2+A_3^2\frac{z^2}{x^2}\bigg) & 0 \\
\frac{A_3^2xz}{2} & 0 & \frac{A_3^2x^2}{2}
\end{bmatrix}, \quad \mathbf{c} = \begin{bmatrix}
-A_4x \\
-A_5 \\
A_4z
\end{bmatrix}.
\end{align}
I impose natural Neumann BC’s and the PDE has initial condition of Dirac masses
P(L=0,x,y,z) = \delta(x)\delta(y)\delta(z-1)
Approximating the length derivative via
\begin{align}
\frac{P^{n+1} - P^n}{\Delta L} \approx (\nabla\cdot \mathbf{a}\nabla + \mathbf{c}\cdot\nabla)P^{n+1}.
\end{align}
Then multiply by the test function v(\mathbf{x}) and integrate over the domain \Omega to obtain
\begin{align}
\int_{\Omega} \bigg(v(\mathbf{x})P^{n+1} - \Delta L\bigg(v(\mathbf{x})\nabla\cdot\mathbf{a}\nabla P^{n+1} + v(\mathbf{x})\mathbf{c}\cdot\nabla P^{n+1}\bigg)\bigg)d\mathbf{x} = \int_{\Omega}v(\mathbf{x})P^{n}d\mathbf{x}.
\end{align}
To formulate the weak form of the problem I integrate by parts the second term to obtain
\begin{align}
\int_{\Omega} v(\mathbf{x})\nabla\cdot\mathbf{a}\nabla P d\mathbf{x} = \int_{\Omega} v(\mathbf{x})\mathbf{n}\cdot \mathbf{a}\nabla P ds - \int_{\Omega} \nabla v(\mathbf{x})\cdot \mathbf{a}\nabla P d\mathbf{x},
\end{align}
with Neumann boundary conditions the surface integral vanishes to obtain
\begin{align}
\int_{\Omega} v(\mathbf{x})\nabla\cdot\mathbf{a}\nabla P d\mathbf{x} = - \int_{\Omega} \nabla v(\mathbf{x})\cdot \mathbf{a}\nabla P d\mathbf{x}.
\end{align}
Hence the weak form is written as
\begin{align}
\int_{\Omega}\bigg( v(\mathbf{x})P^{n+1} + \Delta L\bigg(\nabla v(\mathbf{x})\cdot \mathbf{a}\nabla P^{n+1} - v(\mathbf{x})\mathbf{c}\cdot\nabla P^{n+1}\bigg)\bigg)d\mathbf{x} = \int_{\Omega}v(\mathbf{x})P^{n}d\mathbf{x},
\end{align}
where
\begin{align}
\mathbf{S} = \nabla v(\mathbf{x})\cdot \mathbf{a}\nabla P^{n+1} - v(\mathbf{x})\mathbf{c}\cdot\nabla P^{n+1},
\end{align}
so
\begin{align}
\int_{\Omega}\bigg( v(\mathbf{x})P^{n+1} + \Delta L\mathbf{S}\bigg)d\mathbf{x} = \int_{\Omega}v(\mathbf{x})P^{n}d\mathbf{x},
\end{align}
Hence the variational problem is
\begin{align}
a(P,v) = L_{n+1}(v),
\end{align}
where
\begin{align}
a(P,v) &= \int_{\Omega}\bigg( v(\mathbf{x})P + \Delta L\mathbf{S}\bigg)d\mathbf{x}, \\
L_{n+1}(v) &= \int_{\Omega}v(\mathbf{x})P^{n}d\mathbf{x}.&
\end{align}
\textbf{First Problem}
My main problem is that, one of the properties of the Fokker-Planck equation is that the volume under the PDF must be conserved to 1 over all length steps. I calculate this during each step in the loop via print(assemble(u*dx) and store in a list. Plotting this versus L I want it to stay equal to one. As you can see, it is decaying which is not good
I want to use this probability density function, P, to compute moments - which is what the transmission coefficient is.
\textbf{Second Problem}
The second problem is choosing the upper bounds on the coordinate x,y,z correctly. If I make them too big I require I higher cell density and it uses too much memory. Changing the size of these also affects how the moments are changing which I think stems from the probability volume decaying.
Is there a way to implement a fix or is there always going to be “leakage” of my probability volume. Any thoughts comments greatly appreciated, and thanks for reading.
\mathbf{Code}
from __future__ import print_function
from fenics import *
import numpy as np
import matplotlib.pyplot as plt
L = .1 #Final L
num_steps = 20 #Number of length steps
dl = L / num_steps #Length step size
#Define positive constants
A_1 = 5
A_2 = 6
A_3 = 3
A_4 = 6
A_5 = 5
#Intial mesh values
x0 = 0.01
y0 = 0.01
z0 = 1
#End mesh values
x1 = 3
y1 = 3
z1 = 3
#Number of cells in each direction of the mesh
nxyz = 8
#Number of sample points in P vs x line
xyz_fixed_density = 1000
#Define mesh
mesh = BoxMesh(Point(x0, y0, z0),Point( x1, y1, z1), nxyz, nxyz, nxyz)
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)
#u_D = Expression("(7./pow(cosh(7.*(x[0])),2))*(7./pow(cosh(7.*(x[1])),2))*(7./pow(cosh(7.*(x[2]-1.)),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 = np.array([[0.5*A_3**2*x[2]**2,0,0.5*A_3**2*x[0]*x[2]],
[0,0.5*(A_1**2+A_2**2+A_3**2*(x[2]/x[0])**2),0],
[0.5*A_3**2*x[0]*x[2],0,0.5*A_3**2*x[0]**2]])
cee = np.array([-0.5*A_3**2*x[0]+(A_3**2/2-A_4)*x[0],
-A_5,
-0.5*A_3**2*x[2]+(A_3**2/2+A_4)*x[2]])
a_grad_u = dot(as_matrix(a),grad(u))
c_dot_grad_u = cee.dot(grad(u))
S = dot(grad(v),a_grad_u) - dot(v,c_dot_grad_u)
a_form = (u*v + dl*(S))*dx
L_form = u_n*v*dx
#Starting to build length steps
u = Function(V)
l = 0
#Create empty lists
mean_t = []
L_vals = []
prob_volume = []
u_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)
#Store probability volume
prob_volume.append(assemble(u*dx))
#Update previous solution
u_n.assign(u)
#Print messages
print("We are at length step l =",l,"and we have",num_steps - n, "step(s) remaining." )
#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))
#########################
##Plots
#########################
#Plot mean power transmission versus length
plt.xlabel('$L$')
plt.ylabel('$\\mathbf{E}[\\tau] $')
plt.scatter(L_vals,mean_t,s=1,c='black')
axes = plt.gca()
axes.set_ylim([0,1.1])
plt.savefig('mean_t_vs_L.png')
plt.clf()
#Log plot showing exponential decay
plt.xlabel('$L$')
plt.ylabel('$Log(\\mathbf{E}[\\tau])$')
plt.scatter(L_vals,np.log(mean_t),s=1,c='black')
plt.savefig('log_mean_t_vs_L.png')
plt.clf()
#Plot probability volume vs length
plt.xlabel('$L$')
axes = plt.gca()
axes.set_ylim([0,1.1])
plt.ylabel('$\\int_{\\Omega}P(L,x,y,z)d\\mathbf{x}$')
plt.scatter(L_vals,prob_volume,s=1,c='black')
plt.savefig('prob_vol_vs_L.png')
plt.clf()