Robin type boundary conditions for a Fokker-Planck equation

Main Forum Old-Forum - do not post here
1

Hi,

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)

with BC
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 :wink:

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
Marie-Therese

Attachment: fast_sweep_eikonal.py

Attachment: functions.py

2

Hi,

you are right, you should use the boundary names.

conv += SymbolicBFI (bn*IfPos(bn, fout * u, fin * (1.0-u.Other())) * v, definedon=mesh.Boundaries("outflow|inflow"), skeleton=True)

This line restricts the BFI to the “outflow” and “inflow” boundaries. Thus you have the desired boundary condition on the “wall” as well.

Best,
Christoph

3

Dear Christoph,

thanks a lot - works perfectly :slight_smile:

I have two more questions (sorry)
*) How can I include homogeneous Dirichlet bc at the outflow using the HDG-DG discretization. I know that this is not the best idea, but my colleague wants to see it. I would expect the formation of a boundary layer in this case. I tried

order = 2
fes1 = L2(mesh, order=order)
fes2 = FacetFESpace(mesh, order=order, dirichet=("outflow")) 
fes = FESpace([fes1,fes2])

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, 0, fin * (1.0-u.Other())) * v, definedon=mesh.Boundaries("outflow|inflow"), skeleton=True)

But then I get a rather weird behaviour at the exit. Is there anything else I have to set ? Maybe for the diffusion operator, or do I miss something …

*) I would like to make a movie from the simulation (to convince my colleague in the US, that these boundary conditions are no good indeed ;)). How can I do that ?

Thanks a lot again - I really like the new NgSolve Interface. Great job !
Marie-Therese

4

Hi,

I’m not sure if such a boundary condition makes any sense at an outflow boundary.

For a general flux “b” I would recommend to split the boundary BFI. Otherwise the IfPos wouldn’t do what you want.

conv += SymbolicBFI (bn * fout * u * v, definedon=mesh.Boundaries("outflow"), skeleton=True)
conv += SymbolicBFI (bn * fin * (1.0-u.Other()) * v, definedon=mesh.Boundaries("inflow"), skeleton=True)

To make a video it is the easiest to make snapshots from python.
In functions.py:

from ngsolve import internal

and make the snapshot afer the redraw in “FPIMEX4”:

internal.SnapShot('images/test'+str(index)+'.bmp')

The “internal” has “visoptions” and “viewoptions” to set some GUI parameters. For example to turn of “autoscale” to get the same scaling for the video.

internal.visoptions.autoscale = False

Afterwards you can call “ffmpeg” to generate the video:

ffmpeg -framerate 25 -i images/test%d.bmp video.mp4

Best
Christoph