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