Fokker-Planck equation on a sphere

1

Dear all,

I am trying to solve a stationary density u(\vec x) for a Fokker-Planck equation of the form

-\partial_i a_i(\vec x)u(\vec x)+\frac{1}{2}\partial_i g_{ik}(\vec x)\partial_jg_{jk}(\vec x)u(\vec x) = 0,

where \vec x is on the surface of the sphere. \vec a(\vec x)\in\mathbb{R}^3 and g(\vec x) is a
3\times 2 matrix which determines the diffusion matrix for the process as B=gg^T. Problem has Neumann boundary conditions.

Let \Omega be the surface of a unit sphere in \mathbb{R}^3.
I can state the problem in the variational form as

F=\int_\Omega d\Omega\, u\vec b\cdot\nabla v-\nabla v B\nabla u,

where b_i = a_i -\frac{1}{2}\sum_{j=1}^3\sum_{k=1}^2 g_{ik}(\partial_j g_{jk}).

My problem is that the implementation below returns the “trivial” solution u=C, where
C is some constant, which I can check not to be a stationary state for the original problem.

Below is a minimal working example

from fenics import *
import numpy as np
import matplotlib.pyplot as plt
import mshr as mshr

if __name__ == "__main__":
    

    #Create mesh and define function space
    domain = mshr.Sphere(Point(0, 0,0), 1)
    mesh = mshr.generate_mesh(domain, 20)
    V = FunctionSpace(mesh, 'P', 2)
    

    # Set model parameters
    Omega = 2.0
 
    # Define variational problem
 
    u = TrialFunction(V)
    v = TestFunction(V)

    
    #drift vector
    drift_vector = Expression(('-6*x[0]*pow(x[2],2)','-6*x[1]*pow(x[2],2)-x[2]*O','6*x[2]-6*pow(x[2],3)+x[1]*O'), O=Omega,degree = 1)
    #diffusion matrix
    diffusion_matrix = Expression((('pow(x[1],2)+pow(x[0],2)*pow(x[2],2)','x[0]*x[1]*(pow(x[2],2)-1)','x[0]*x[2]*(pow(x[2],2)-1)'),
                                       ('x[0]*x[1]*(pow(x[2],2)-1)','pow(x[0],2)+pow(x[1],2)*pow(x[2],2)','x[1]*x[2]*(pow(x[2],2)-1)'),
                                       ('x[0]*x[2]*(pow(x[2],2)-1)','x[1]*x[2]*(pow(x[2],2)-1)','pow(pow(x[2],2)-1,2)')),degree=1)

    
    #variational form
    F =  u*dot(drift_vector,grad(v))*dx +dot(grad(v),dot(diffusion_matrix,grad(u)))*dx 
    a, L  = lhs(F), rhs(F)

    #solution
    u = Function(V)
    solve(a==L, u)#J=J)

    #plot the solution
    plot(u)
    plt.show()

Any hints and help would be greatly appreciated.
Best,
Kimmo

2

Can you expand on this?

3

If \tilde u(\vec x) satisfies \partial_i \tilde u(\vec x)=0, i.e. is a constant, and I plug it in to the original Fokker-Planck equation it will not satisfy it. That is,

-\partial_i a_i\tilde u+\frac{1}{2}\partial_i g_{ik}\partial_j g_{jk}\tilde u \neq 0.