Solve transport problem with open/closed boundary

[jw3126]

I would like to compute the steady-state of a particle density ρ, that is moved under a velocity field v and is produced at a rate s.

-div(vρ) + s = 0

The boundary condition is that particles are free to leave the simulation volume at boundaries, but no particles can enter the simulation that way. (Does this boundary condition have a name?)

I tried

using VoronoiFVM

function flux!(f,u,edge)
    v = 1.0
    h = meas(edge)
    f[1] = if v > 0
        v*u[1, 1] * h
    else
        v*u[1, 2] * h
    end
end

function source!(out, node)
    out[1] = 1
end

physics=VoronoiFVM.Physics(num_species=1, 
    flux=flux!,
    source=source!,
)

X = range(0, stop=1, length=5)
grid=VoronoiFVM.Grid(collect(X))
sys=VoronoiFVM.System(grid,physics)

ispec = 1
enable_species!(sys,ispec,[1])
boundary_neumann!(sys, ispec, 1, 0.0) # this should ensure no particles enter from the left

inival=unknowns(sys,inival=0)
solution=unknowns(sys)
solve!(solution,inival,sys)

But it errors:

ERROR: LoadError: LinearAlgebra.SingularException(0)

How to solve this problem with VoronoiFVM.jl?

Background

The above is a simplified example of an issue I face. Really I would like to simulate an ionization chamber. The equations are:

E = -grad(φ)
ε*div(E) = - ρₑ - ρ₋ + ρ₊

∂ₜρₑ =  div(μₑ(E)*E*ρₑ) + Rₑ(E,ρₑ,ρ₋,ρ₊) + sₑ
∂ₜρ₋ =  div(μ₋(E)*E*ρ₋) + R₋(E,ρₑ,ρ₋,ρ₊) + s₋
∂ₜρ₊ = -div(μ₊(E)*E*ρ₊) + R₊(E,ρₑ,ρ₋,ρ₊) + s₊

• φ is the potential
• E is the electric field
• ρₑ, ρ₋, ρ₊ are the densities of electrons, negative and positive ions
• ε is the permittivity.
• μᵢ are the mobilities of the particle species. In the case of electrons there is a strong dependence on E.
• Rᵢ are the reaction rates. For instance, an electron could attach to and O₂ molecule and become a negative ion.
  Reactions involving electrons have a dependence on E.
• sᵢ are the source terms. They capture the liberation of charge carriers due to radiation.

The unknowns are (φ, E, ρₑ, ρ₋, ρ₊ ). The other functions/constants are known. I am interested in both time evolution and steady-state (∂ₜρᵢ = 0).

Boundary conditions

The boundary decomposes into parts that are conductive Γ₁ and parts that are not Γ₂:

∂Ω = Γ₁ ∪ Γ₂

• At the conductive parts Γ₁ we have a Dirichlet condition on the potential φ
• At the non-conductive parts Γ₂ we probably want zero von Neumann condition on E

The boundary condition for the charge carriers ρᵢ is:

• At the conductive parts Γ₁ charge carriers can leave the chamber if their sign is appropriate. They can never enter the chamber.
• At the non-conductive parts Γ₂ charge carriers can neither leave nor enter.

[j-fu]

Hi,

your system is very close to what we are doing for semiconductor devices (and, recently, for perovskite solid electrolytes), just without the diffusion. Therefore, the transport part is hyperbolic, and one disclaimer appears to be necessary here: the method is quite bad at capturing sharp shocks as upwinding introduces artificial diffusion. So you need to check if your solution is good enough in this sense.

The singular exception happens because boundary_neumann!(sys, ispec, ibc, 0.0) is the default if nothing is specified, because in the implementation (like in the weak formulation) this case coincides with not having any boundary terms at all.
So in fact we have Neumann bc on both ends and therefore a singular matrix.

I need to think about the bc at the conductive parts. Do you have a formal description for these ? I worked with outflow bc for fluid pronlems, there might be similarities. These are not implemented yet here.

As for the field depedency, this can be done with the current version in 1D. For higher dimensions, consistent electric field reonstruction from the two point flux gradients needs a bit of effort. In our C++ code, we did use the gradient of the P1 FEM reconstruction of the solution averaged from the elements having one edge in common.

[jw3126]

Thanks a lot for your response!

Therefore, the transport part is hyperbolic, and one disclaimer appears to be necessary here: the method is quite bad at capturing sharp shocks as upwinding introduces artificial diffusion. So you need to check if your solution is good enough in this sense.

There can be indeed sharp shocks. However currently the main quantity of interest for me is collection efficiency. E.g. understanding how many charged particles get collected at the boundaries vs how many neutralize in the chamber. And artificial diffusion does not have a large influence on this.

The singular exception happens because boundary_neumann!(sys, ispec, ibc, 0.0) is the default if nothing is specified, because in the implementation (like in the weak formulation) this case coincides with not having any boundary terms at all.
So in fact we have Neumann bc on both ends and therefore a singular matrix.

Ah thanks!

I need to think about the bc at the conductive parts. Do you have a formal description for these ?

No I do not have a formal description. It would be good to have this formalized better. Operationally I think what needs to be done is

1. Estimate the E field at the boundary (e.g. from two point potential gradients between the boundary cell and inner neighbor cells)
2. Compute the outflow ∫μ(E)Eρn ds or 0 depending on the direction of E
3. Subtract the outflow from the equation.

As for the field depedency, this can be done with the current version in 1D.

Ah great, can you show how to inform VoronoiFWM about the E field dependency in 1D?

For higher dimensions, consistent electric field reonstruction from the two point flux gradients needs a bit of effort.
In our C++ code, we did use the gradient of the P1 FEM reconstruction of the solution averaged from the elements having one edge in common.

So in each timestep you solve a FEM problem for the field and then a FVM problem for charge transport? It seems to me that this leads to a dense jacobian of the ODE. Is that a problem?

What I thought was to estimate E in the cell by solving the least squares system E⋅(pt2 - pt) = ϕ(pt) - ϕ(pt2) where pt2 runs over all neighboring cells. Is this a sane approach?