Normalization decreases with Neumann boundary condition

variational formulation
1

Hi, I am solving a Fokker-Planck equation using FEniCS, the solution to the PDE is the time-dependent probability density function (pdf) of the magnitude and phasor of the variable \rho (a_m, a_{\phi};t), in principle, the pdf should also as a normalization that equals to 1 at any time. I tried to apply the zero-flux boundary condition (which I think is Neumann boundary condition) on a_m and periodical boundary condition on a_{\phi}, however, when I solve the PDE I can see the normalization of the solution is decreasing with time and I wonder if you can help me check what is going wrong, some parts of my code are as following:

class PeriodicBoundary(SubDomain):
    def inside(self, x, on_boundary):
        return bool(x[1] < DOLFIN_EPS and x[1] > -DOLFIN_EPS and on_boundary)

    # Map top boundary (H) to left boundary (G)
    def map(self, x, y):
        y[1] = x[1] - 2 * np.pi
        y[0] = x[0]
...
pbc = PeriodicBoundary()
V = FunctionSpace(mesh, 'P', 1, constrained_domain=pbc)
...
for n in range(1, num_steps):
    t += dt
    u_D.t = t

    solve(a == L, u,  solver_parameters={'linear_solver':'mumps', 'preconditioner': 'ilu'})

    norm = assemble(u*dx)
    u_n.assign(u)

If I print out norm at each time step I could see the decrease of it, from 1 to 0 after a long time.