Density-dependent diffusion

[mjyshin]

Hello, I had asked a question a while ago about time and space dependence in the flux and reaction terms. I'm not sure if that has been resolved, but I was also wondering how one would go about modeling a flux term like j = -D(u/K)ᵐ∇u? This was my attempt (by averaging u[1,1] and u[1,2]), but I'm not sure about it:

function g!(g,u,∂ω,θ)
    ū = (u[1,1]+u[1,2])/2
    g[1] = -θ.D*(ū/θ.K)^θ.m*(u[1,2]-u[1,1])
end

Or more generally what would I do for j = -f(u)∇u? Also, since j⋅n̂ is approximated by g(uₖ,uₗ)/|xₖ-xₗ|, wouldn't we need to multiply the u by something like (abs ∘ -)(∂ω.coord[1:2]...)? Thank you in advance!

[j-fu]

Hello, I had asked a question a while ago about time and space dependence in the flux and reaction terms. I'm not sure if that has been resolved,

Hi,
you can access the coordinates of nodes in reaction terms end edges in flux functions
via their AbstractArray API: node[1] would be the x coordinate of the corresponding node,
and edge[1,1] resp. edge[1,2] are the coordinates of of the edge ends.
For time dependency, since 0.14. you can access time(node), time(edge).
There will be interface improvements and more examples for transient soltutions in 0.15 which is essentially ready.

but I was also wondering how one would go about modeling a flux term like j = -D(u/K)ᵐ∇u? This was my attempt (by averaging u[1,1] and u[1,2]), but I'm not sure about it:

function g!(g,u,∂ω,θ)
    ū = (u[1,1]+u[1,2])/2
    g[1] = -θ.D*(ū/θ.K)^θ.m*(u[1,2]-u[1,1])
end

Or more generally what would I do for j = -f(u)∇u?

This is essentially done in Example106. With f=F', i.e. F being an integral of f, you would write g(u_k, u_l) = F(u_k) - F(u_l) . This can be interpreted as being derived from the exact solution of the flux along the edge. For motivation and reference see https://link.springer.com/content/pdf/10.1007/s00211-005-0659-5.pdf .
If F is not available, your idea is however reasonable.

Also, since j⋅n̂ is approximated by g(uₖ,uₗ)/|xₖ-xₗ|, wouldn't we need to multiply the u by something like (abs ∘ -)(∂ω.coord[1:2]...)? Thank you in advance!

I don't understand this question completely. Geometry dependent scaling (division by point distance and multiplication by interface width) is done outside of the flux functions.

[mjyshin]

I feel like your FVM solver would be a good method to use to solve the (equilibrium) Fokker-Planck equation since probability mass is conserved:

∂ₜρ(x,t) = 0 = ∇⋅(∇V(x)ρ(x,t)) + D:∇²ρ(x,t)

If we let diffusion matrix D be diagonal, and we define a Delaunay triangulated mesh (not necessarily a grid), could I perhaps define a flux that projects the drift and diffusion onto the edge in this way?

function jdx(jdx,ρ,∂ω,θ)
    # I have data θ that passes a function θ.∇V that computes the gradient given x and parameters θ.θ
    x1,x2 = ∂ω[:,1],∂ω[:,2]    # I noticed edge[:,1] and edge[:,2] give the coordinates of the centers of Voronoi regions, not edge verteces?
    ρ1,ρ2 = ρ[1,1],ρ[1,2]
    x = (x1+x2)/2    # point x on edge ∂ω between the two Voronoi regions
    dx = norm(x1-x2)
    ∇Vdx = project(∂ω,θ.∇V(x,θ.θ))    # projecting vector ∇V(x) on edge
    D = project(∂ω,diag(θ.D))/dx    # projecting magnitude of diffusion in x and y directions on edge, but want D not Ddx
    
    -∇Vdx > 0 ? jdx1 = -∇Vdx*ρ1 : jdx1 = -∇Vdx*ρ2
    jdx2 = -D*(ρ2-ρ1)
    jdx[1] = jdx1 + jdx2
end

I used the upwind method here, because I couldn't get the Bernoulli function method to even run, but even this doesn't give a solution I'd expect... Would a problem like this be solvable using your solver?

[j-fu]

Hi, in principle it could work.

I noticed edge[:,1] and edge[:,2] give the coordinates of the centers of Voronoi regions, not edge verteces?

Yes. In fact in 2D, the notation "edge" is a bit ambiguous. I use edge for the edge of the underlying Delaunay triangulation, whose vertices conincide with the control volume centers.
In 3D it still would be an edge while the inteface between the control volumes is a 2D object.

With ansiotropic diffusion however you should be careful. This would work correctly only on grids aligned with the anisotropy directions aligned to the grid lines (I think I have a remark on that in the docs). I have a master student writing a thesis on a finite volume approach in the more general case, I hope to have some Julia implementation before summer.