I wrote a solver for a Fokker-Planck equation a couple of month ago. The equation is a minimalistic model for unidirectional pedestrian flow in a corridor (rectangle) - people enter the corridor on the left and exit on the right. The upper and lower boundaries correspond to walls. We used an HDG discretization for the diffusive operator and a DG for the convection. The boundary conditions were a bit tricky back then and Joachim sent me a ‘quick fix’ (which worked fine :-D)
Now I want to consider a more complicated geometry - my corridor has a bottleneck - and this quick fix is not working anymore (people enter at the vertical boundaries or even all boundaries,…) So in principle I need to figure out how to include the Robin boundary conditions correctly.
The equation is as follows (sorry for the latex, but I was not able to figure out how to type equations here). Given a potential (either b = (vmax,0) or the gradient of the function satisfying the Eikonal equation) and inflow and outflow rates fin and fout solve:
u_t = div(\nabla u + u (1-u) b)
Entrance (inflow): (\nabla u + u (1-u) b) n = fin (1-u)
Exit (outflow): (nabla u + u(1-u) b) n = fout u
Walls: (\nabla u + u (1-u) b) n = 0
The ‘old fix’ looks as follows:
# Transportation vector in x-direction #b = CoefficientFunction( (v_run,0) ) b = -grad(gfphi) bn = b*specialcf.normal(2) # Discretization of the convection operator - DG formulation with upwind conv = BilinearForm(fes) conv += SymbolicBFI (-u * (1-u) * b* grad(v)) conv += SymbolicBFI (bn*IfPos(bn, u * (1-u), u.Other() * (1-u.Other())) * (v-v.Other()), VOL, skeleton=True) conv += SymbolicBFI (bn*IfPos(bn, fout * u, fin * (1.0-u.Other())) * v, BND, skeleton=True)
I changed my geometry and labelled the boundary segments as ‘inflow’, ‘outflow’ and ‘wall’. But I don’t use this anywhere - which I guess is the reason why it is not working
I also attached the full code in case you want to have a closer look at it. In principle I define the geometry (corridor with bottleneck), solve the eikonal equation using a fast sweeping method proposed by Quian et al, SIAM J Num Anal 45(1) 83-107 2007 (I could do this more efficiently I assume but that’s not the most important part now), take the gradient of the potential as velocity field and then solve the parabolic FPE using an IMEX scheme - fast_sweeping_eikonal.py is the main, while functions.py includes all functions used.
Thanks a lot