Some problems with Voronoi FVM

Source

Draft. J. Fuhrmann, Oct. 29. 2021. Updated Dec 19, 2021.

We discuss one of the critical cases for application the Voronoi finite volume method. We provide some practical fix and opine that the finite element method probably has the same problems.

begin
    import Pkg as _Pkg
    haskey(ENV, "PLUTO_PROJECT") && _Pkg.activate(ENV["PLUTO_PROJECT"])
    using Revise
    using Test
    using VoronoiFVM
    using ExtendableGrids
    using PlutoUI
    using GridVisualize
    using CairoMakie
    CairoMakie.activate!(; type = "png", visible = false)
    GridVisualize.default_plotter!(CairoMakie)
end;

Transient problem

This problem was suggested by R. Eymard.

Regard the following problem coupling Darcy's equation with Fick's law and transport:

$$ \begin{aligned} \vec v &= k \nabla p \\ \nabla \cdot \vec v &= 0\\ \partial_t (\phi c) - \nabla \cdot (D\nabla c + c \vec v) &= 0 \end{aligned}$$

The domain is described by the following discretization grid:

In the center of the domain, we assume a layer with high permeability.

As boundary conditions for the pressure \(p\) we choose fixed pressure values at the left and right boundaries of the domain, triggering a constant pressure gradient throughout the domain.

At the inlet of the high permeability layer, we set \(c=1\), and at the outlet, we set \(c=0\).

The high permeability layer has length L=10 and width W= 0.5.

We solve the time dependent problem on three types of rectangular grids with the same resolution in \(x\) direction and different variants to to handle the high permeability layer.

  • grid_n - a "naive" grid which just resolves the permeability layer and the surrounding material with equally spaced (in y direction) grids

  • grid_1 - a 1D grid of the high permeability layer. With high permeability contrast, the solution of the 2D case at y=0 should coincide with the 1D solution

  • grid_f - a "fixed" 2D grid which resolves the permeability layer and the surrounding material with equally spaced (in y direction) grids and "protection layers" of width ε_fix=0.0001 correcting the size of high permeability control volumes

Results

In the calculations, we ramp up the inlet concentration and measure the amount of solute flowing through the outlet - the breaktrough curve.

nref = 1
1
tend = 100
100
ε_fix = 1.0e-4
0.0001
grid_n, sol_n, bt_n = trsolve(grid_2d(; nref = nref); tend = tend);
sum(bt_n)
18.143158169851787
@test sum(bt_n) ≈ 18.143158169851787
Test Passed
grid_1, sol_1, bt_1 = trsolve(grid_1d(; nref = nref); tend = tend);
@test sum(bt_1) ≈ 20.66209910195916
Test Passed
grid_f, sol_f, bt_f = trsolve(grid_2d(; nref = nref, ε_fix = ε_fix); tend = tend);
@test sum(bt_f) ≈ 20.661131375044135
Test Passed
grid_ϕ, sol_ϕ, bt_ϕ = trsolve(grid_2d(; nref = nref); ϕ = [1.0e-3, 1], tend = tend);
@test sum(bt_ϕ) ≈ 20.412256299447236
Test Passed
begin
    p1 = GridVisualizer(; resolution = (500, 200),
                        xlabel = "t",
                        ylabel = "outflow",
                        legend = :rb,
                        title = "Breakthrough Curves")
    scalarplot!(p1, sol_n.t, bt_n; label = "naive grid", color = :red)
    scalarplot!(p1,
                sol_1.t,
                bt_1;
                label = "1D grid",
                markershape = :x,
                markersize = 10,
                clear = false,
                color = :green)
    scalarplot!(p1,
                sol_f.t,
                bt_f;
                label = "grid with fix",
                markershape = :circle,
                color = :green,
                clear = false)
    scalarplot!(p1,
                sol_ϕ.t,
                bt_ϕ;
                label = "modified ϕ",
                markershape = :cross,
                color = :blue,
                clear = false)
    reveal(p1)
end

Here, we plot the solutions for the grid_n case and the grid_f case.

Time: 10.0

scalarplot(grid_n, sol_n(t)[ic, :]; resolution = (500, 200), show = true)
scalarplot(grid_f, sol_f(t)[ic, :]; resolution = (500, 200), show = true)

Reaction-Diffusion problem

Here we solve the following problem:

$$ -\nabla D \nabla u + R u = 0$$

where D is large in the high permeability region and small otherwise. R is a constant.

Results

rdgrid_1, rdsol_1, of_1 = rdsolve(grid_1d(; nref = nref));
@test of_1 ≈ -0.013495959676585267
Test Passed
rdgrid_n, rdsol_n, of_n = rdsolve(grid_2d(; nref = nref));
@test of_n ≈ -0.00023622450350365264
Test Passed
rdgrid_f, rdsol_f, of_f = rdsolve(grid_2d(; nref = nref, ε_fix = ε_fix));
@test of_f ≈ -0.013466874615165499
Test Passed
rdgrid_r, rdsol_r, of_r = rdsolve(grid_2d(; nref = nref); R = [0, 0.1]);
@test of_r ≈ -0.013495959676764535
Test Passed

We measure the outflow at the outlet. As a result, we obtain:

  • 1D case: -0.013495959676585255

  • 2D case, naive grid: -0.00023622450350365277

  • 2D case, grid with "protective layer": -0.013466874615165514

  • 2D case, naive grid, "modified" R: -0.013495959676764539

scalarplot(rdgrid_1, rdsol_1; resolution = (300, 200))
scalarplot(rdgrid_n, rdsol_n; resolution = (500, 200))
scalarplot(rdgrid_f, rdsol_f; resolution = (500, 200))

Discussion

Transient case

As there will be nearly no flow in y-direction, we should get the very same results in all four cases for small permeability values in the low permeability region.

In the grid_n case, the heterogeneous control volumina ovrestimate the storage capacity which shows itself in a underestimation of the transferred solute.

With the high permeability contrast, the results for heterogeneous domain should be essentially equal to those for 1D domain. However, with a coarse resolution in y-direction, we see large differences in the transient behaviour of the breaktrough curve compared to the 1D case. The introduction of a thin protection layer leads to reasonable results.

Alternatively, the porosity of the low permeability region can be modified. Arguably, this is the case in practice, see e.g. Ackerer et al, Transport in Porous Media35:345–373, 1999 (eq. 7).

Reaction diffusion case

In this case, we look at a homogeneous reaction in a domain divided into a high and low diffusion region. With high contrast in the diffusion coefficients, the reasonable assumption is that the reaction takes place only in the high diffusion region, and the un-consumed share of species leaves at the outlet.

In this case we observe a similar related problem which can be fixed by adding a thin layer of control volumes at the boundary. No problem occurs if the reaction rate at in the low diffusion region is zero.

Conclusion

Here, we indeed observe problem with the Voronoi approach: care must be taken to handle the case of hetero interfaces in connection with transient processes and/or homogeneous reactions. In these cases it should be analyzed if the problem occurs, and why, and it appears, that the discussion should not be had without reference to the correct physical models. A remedy based on meshing is available at least for straight interfaces.

Opinion

With standard ways of using finite elements, the issue described here will occur in a similar way, so the discussion is indeed with the alternative "cell centered" finite volume approach which places interfaces at the boundaries of the control volumes rather than along the edges of a underlying triangulation.

Drawbacks of two point flux Voronoi methods based on simplicial meshes (as tested here):

  • Anisotropic diffusion is only correct with aligned meshes

  • Reliance on boundary conforming Delaunay property of the underlying mesh, thus narrowing the available meshing strategies

  • The issue described in the present notebook. However, in both cases discussed here, IMHO it might "go away" depending on the correct physics. There should be more discussions with relevant application problems at hand.

Advantages (compared to the cell centered approach placing collocation points away from interfaces)

  • Availability of P1 interpolant on simplices for visualization, interpolation, coupling etc.

  • Mesh generators tend to place interfaces at triangle edges.

  • Dirichlet BC can be applied exactly

  • There is a straightforward way to link interface processes with bulk processes, e.g. an adsorption reaction is easily described by a reaction term at the boundary which involves interface and bulk value available at the same mesh node.

Appendix

Domain data

Sizes:

begin
    L = 10   # length of the high perm layer
    W = 0.5  # width of high perm layer
    Wlow = 2 # width of adjacent low perm layers
end;

Boundary conditions:

begin
    const Γ_top = 3
    const Γ_bot = 1
    const Γ_left = 4
    const Γ_right = 2
    const Γ_in = 5
    const Γ_out = 2
end;
begin
    Ω_low = 1
    Ω_high = 2
end;
function grid_2d(; nref = 0, ε_fix = 0.0)
    nx = 10 * 2^nref
    ny = 1 * 2^nref
    nylow = 3 * 2^nref
    xc = linspace(0, L, nx + 1)
    y0 = linspace(-W / 2, W / 2, ny + 1)
    if ε_fix > 0.0
        yfix = [W / 2, W / 2 + ε_fix]
        ytop = glue(yfix, linspace(yfix[end], Wlow, nylow + 1))
    else
        ytop = linspace(W / 2, Wlow, nylow + 1)
    end
    yc = glue(-reverse(ytop), glue(y0, ytop))
    grid = simplexgrid(xc, yc)
    cellmask!(grid, [0, -W / 2], [L, W / 2], Ω_high)
    bfacemask!(grid, [0, -W / 2], [0, W / 2], Γ_in)
    bfacemask!(grid, [L, -W / 2], [L, W / 2], Γ_out)
end
grid_2d (generic function with 1 method)
function grid_1d(; nref = 0)
    nx = 10 * 2^nref
    xc = linspace(0, L, nx + 1)
    grid = simplexgrid(xc)
    cellmask!(grid, [0], [L], Ω_high)
    bfacemask!(grid, [0], [0], Γ_in)
    bfacemask!(grid, [L], [L], Γ_out)
    grid
end
grid_1d (generic function with 1 method)

Transient solver

Pressure index in solution

const ip = 1;

Concentration index in solution

const ic = 2;

Generate breaktrough courve from transient solution

function breakthrough(sys, tf, sol)
    of = similar(sol.t)
    t = sol.t
    of[1] = 0
    for i = 2:length(sol.t)
        of[i] = -integrate(sys, tf, sol.u[i], sol.u[i - 1], t[i] - t[i - 1])[ic]
    end
    of
end
breakthrough (generic function with 1 method)

Transient solver:

function trsolve(grid;
                 κ = [1.0e-3, 5],
                 D = [1.0e-12, 1.0e-12],
                 Δp = 1.0,
                 ϕ = [1, 1],
                 tend = 100,)
    function flux(y, u, edge, data)
        y[ip] = κ[edge.region] * (u[ip, 1] - u[ip, 2])
        bp, bm = fbernoulli_pm(y[ip] / D[edge.region])
        y[ic] = D[edge.region] * (bm * u[ic, 1] - bp * u[ic, 2])
    end

    function stor(y, u, node, data)
        y[ip] = 0
        y[ic] = ϕ[node.region] * u[ic]
    end

    dim = dim_space(grid)
    function bc(y, u, bnode, data)
        c0 = ramp(bnode.time; dt = (0, 0.001), du = (0, 1))
        boundary_dirichlet!(y, u, bnode, ic, Γ_in, c0)
        boundary_dirichlet!(y, u, bnode, ic, Γ_out, 0)

        boundary_dirichlet!(y, u, bnode, ip, Γ_in, Δp)
        boundary_dirichlet!(y, u, bnode, ip, Γ_out, 0)
        if dim > 1
            boundary_dirichlet!(y, u, bnode, ip, Γ_left, Δp)
            boundary_dirichlet!(y, u, bnode, ip, Γ_right, 0)
        end
    end

    sys = VoronoiFVM.System(grid;
                            check_allocs = true,
                            flux = flux,
                            storage = stor,
                            bcondition = bc,
                            species = [ip, ic],)

    inival = solve(sys; inival = 0, time = 0.0)
    factory = TestFunctionFactory(sys)
    tfc = testfunction(factory, [Γ_in, Γ_left, Γ_top, Γ_bot], [Γ_out])

    sol = VoronoiFVM.solve(sys;
                           inival = inival,
                           times = [0, tend],
                           Δt = 1.0e-4,
                           Δt_min = 1.0e-6,)

    bt = breakthrough(sys, tfc, sol)
    if dim == 1
        bt = bt * W
    end

    grid, sol, bt
end
trsolve (generic function with 1 method)

Reaction-Diffusion solver

function rdsolve(grid; D = [1.0e-12, 1.0], R = [1, 0.1])
    function flux(y, u, edge, data)
        y[1] = D[edge.region] * (u[1, 1] - u[1, 2])
    end

    function rea(y, u, node, data)
        y[1] = R[node.region] * u[1]
    end
    function bc(y, u, bnode, data)
        boundary_dirichlet!(y, u, bnode, 1, Γ_in, 1)
        boundary_dirichlet!(y, u, bnode, 1, Γ_out, 0)
    end
    sys = VoronoiFVM.System(grid;
                            flux = flux,
                            reaction = rea,
                            species = 1,
                            bcondition = bc,
                            check_allocs = true)
    dim = dim_space(grid)

    sol = VoronoiFVM.solve(sys)
    factory = TestFunctionFactory(sys)
    tf = testfunction(factory, [Γ_in, Γ_left, Γ_top, Γ_bot], [Γ_out])
    of = integrate(sys, tf, sol)
    fac = 1.0
    if dim == 1
        fac = W
    end
    grid, sol[1, :], of[1] * fac
end
rdsolve (generic function with 1 method)

Built with Julia 1.11.1 and

CairoMakie 0.12.16
ExtendableGrids 1.9.2
GridVisualize 1.7.0
Pkg 1.11.0
PlutoUI 0.7.60
Revise 3.5.18
Test 1.11.0
VoronoiFVM 1.25.1