NaN solution problem + general help with my implementation

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.

Update…Fixed NaN solution by changing mesh starting value to avoid 0/0. It’s basically conserving the probability now. Sorry for wasting anyones time reading this.

Hi bacsieboy, I am currently also working on solving Fokker-Planck equation (FPE) with Fenics and got the NaN issue. Could you please share a little more about how you solve the prolem? Thanks a lot!