Mean‐Field Games “Crowd Control” Problem in Julia: HJB OK, but Struggling With Fokker‐Planck Positivity in Gridap

Specific Domains Numerics
,
1

Hi everyone,

I’m relatively new to Julia and I’m trying to solve a crowd control Mean Field Game (MFG) problem (inspired by Yves Achdou’s work " Mean Field Games for Modeling Crowd Motion") consisting of:

  • A Hamilton–Jacobi–Bellman (HJB) equation
  • A Fokker–Planck (FP) equation

The domain is 2D with a somewhat complex geometry – a room with rectangular obstacles inside. The boundary includes two “doors” (for outflow) with Dirichlet=0, and walls with no-flux (Neumann=0) conditions.

I’ve managed to:

  • Generate the mesh with Gmsh and import it into Gridap.jl,
  • Set up and solve the HJB part successfully (with some standard boundary conditions: Dirichlet=0 on doors, Neumann=0 on walls).

However, in the FP equation, I need the density ρ to remain nonnegative. A standard FEM approach in Gridap doesn’t inherently enforce positivity. I tried transformations like ϕ=log⁡(1+ρ) to handle the positivity, but end up hitting singularity errors.

Questions:

  • Should I remain in Gridap (with a specialized approach for positivity)?
  • Or switch to another FEM/FVM package (like Ferrite.jl, VoronoiFVM.jl, etc.) that might better handle positivity constraints?
  • Any recommended techniques for the MFG context?

My ultimate goal is to run a standard fixed‐point iteration (Achdou style) in which HJB at iteration n yields a drift used in the FP at iteration n+1, and vice versa, until convergence.

I really like Gridap’s high‐level interface for PDE operators, but am open to alternatives if there’s a more robust way to ensure ρ≥0. For instance, I tried code like this:

a_ρ(t, dtρ, v) = ∫(v *dtρ)dΩ  
f_ρ(t, ρ) = ρ/(1+ρ).^(3/4)
b_ρ(t, ρ, v) = 0.05*∫( ∇(ρ)⋅∇(v) )dΩ + ∫(v*(16*(∇(u(t))⋅∇(v))*f_ρ(t, ρ)))dΩ


jac_ρ_t(t, ρ, dtρ, v) =  ∫(v * dtρ)dΩ 
∂f_ρ(t,ρ) = (1-ρ/4)/(1+ρ).^(7/4)
jac_ρ(t, ρ, dρ, v) =  0.05*∫( ∇(dρ)⋅∇(v) )dΩ + ∫(16*(∇(u(t))⋅∇(v))*∂f_ρ(t,ρ)*dρ)dΩ

res_ρ(t, ρ, v) = b_ρ(t, ρ, v) 
op_ρ = TransientSemilinearFEOperator( a_ρ, res_ρ, (jac_ρ, jac_ρ_t), U, V; constant_mass=true)

lin_solver = LUSolver()
nl_solver = NLSolver(lin_solver, method=:newton, iterations=10, show_trace=false)

Δt = 0.05
θ = 0.5
solver = ThetaMethod(nl_solver, Δt, θ)

t0, T = 0.0, 10.0
ρh0 = interpolate_everywhere(m_0, V)
ρh =@pipe solve(solver, op_ρ, t0, T, ρh0)  |> collect(_)

…and it fails due to ρ going negative.
Furthermore, I tried studying the PDE solved by ϕ=log⁡(1+ρ) with

e=exp(1.0)

a_ϕ(t, ϕ, dtϕ, v) = ∫(v *dtϕ*e.^ϕ)dΩ  
f_ϕ(t, ϕ) = e.^(ϕ/4)-e.^(-ϕ*3/4)
b_ϕ(t, ϕ, v) = 0.05*∫( ∇(ϕ)⋅∇(v)*e.^ϕ )dΩ + ∫(v*(16*(∇(u(t))⋅∇(v))*f_ϕ(t, ϕ)))dΩ


jac_ϕ_t(t, ϕ, dtϕ, v) =  ∫(v * dtϕ * e.^ϕ)dΩ 
∂f_ϕ(t, ϕ) = e.^(ϕ/4)/4+e.^(-ϕ*3/4)*3/4
jac_ϕ(t, ϕ, dϕ, v) =  0.05*∫( ∇(dϕ)⋅∇(v)*e.^ϕ )dΩ + 0.05*∫( ∇(ϕ)⋅∇(v)*e.^ϕ*dϕ )dΩ + ∫(16*(∇(u(t))⋅∇(v))*∂f_ϕ(t, ϕ)*dϕ)dΩ

res_ϕ(t, ϕ, v) = b_ϕ(t, ϕ, v) 
op_ϕ = TransientQuasilinearFEOperator( a_ϕ, res_ϕ, (jac_ϕ, jac_ϕ_t), U, V)

lin_solver = LUSolver()
nl_solver = NLSolver(lin_solver, method=:newton, iterations=10, show_trace=false)

Δt = 0.05
θ = 0.5
solver = ThetaMethod(nl_solver, Δt, θ)

t0, T = 0.0, 10.0
ϕh0 = interpolate_everywhere(x -> log(m_0(x)+1), V)
ϕh =@pipe solve(solver, op_ϕ, t0, T, ϕh0) |> collect(_)

…and it fails due to a raising SingularException.
Any advice or suggested approach? Thanks!

2

Not sure about the singular exception, but concerning positivity you would need to do some upwinding/stabilization in the Fokker-Planck term, like streamline Petrov-Galerking or algebraic flux correction. Not sure if Gridap or Ferrite have options for this. I guess Ferrite is transparent enough in its internals to allow to implement this.

As for VoronoiFVM (I am the author) - this code would allow the upwinding in the Fokker-Planck equation with positivity guarantees. However it seems that there is a factor depending on the norm of the gradient of the solution which in the moment cannot be implemented in VoronoiFVM, as it cannot be calculated just from the two values passed to the flux functions. As for future versions, there are ideas how to do it ( see e.g. https://www.wias-berlin.de/preprint/2378/wias_preprints_2378.pdf, that was in C++…), but no roadmap so far.

1 Like
3

Hi again!

Progress since my first post

  • I stayed with Gridap.jl for the HJB part; that piece still works nicely on my Gmsh mesh.
  • Following @j-fu’s suggestion I explored VoronoiFVM.jl (very cool package, thank you!). Unfortunately I didn’t manage to adapt the advanced positivity-preserving ideas linked in your reply to my particular nonlinear divergence term, so ρ still dips negative.

What I tried next

Back in Gridap I focused again on the log-shift variable

\rho=e^{\phi}−1,

hoping to guarantee non-negativity.
With Automatic Differentiation inside TransientQuasilinearFEOperator the code runs, but it’s slow.
When I replace AD with hand-written Jacobians I immediately hit a SingularException during the first Newton factorisation. I suspect my mass-term Jacobian is wrong: the mass term is

a_\phi(t,\phi, \partial_t\phi, v) = \int_\Omega v(t,x) \partial_t\phi(t,x) e^{\phi(t,x)}\, dx ​,

so a_\phi depends on the state itself and I’m unsure how to tell TransientQuasilinearFEOperator its state-varying Jacobian correctly.

Concrete question

How do you supply a state-dependent mass Jacobian to TransientQuasilinearFEOperator?
A minimal snippet (or rule-of-thumb) for something like

jac_ϕ_t(t, ϕ, dϕ, dtϕ, v) = ∫ v * dtϕ * exp(ϕ)dϕ dΩ

would be a huge help.

Thanks again to @j-fu (and anyone else who chimes in)!