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

CairoMakie

Interface conditions in 1D

Source

This notebooks discusses handling of internal interfaces with VoronoiFVM.jl.

Two subdomains

For a simple stationary diffusion equation with an interior interface, we discuss possible interface conditions between two subdomains.

Let \(\Omega=\Omega_1\cup\Omega_2\) where \(\Omega_1=(-1,0)\) and \(\Omega_2=(0,1)\). Let \(\Gamma_1={-1}\),\(\Gamma_2={1}\) and \(\Gamma_3={0}\).

Regard the following problem:

\(\begin{aligned} -\Delta u_1 &= 0 & \text{in}\quad \Omega_1\\ -\Delta u_2 &= 0 & \text{in}\quad \Omega_2\\ \end{aligned}\)

with exterior boundary conditions

\(u_1|_{\Gamma_1} = g_1\) and \(u_2|_{\Gamma_2} = g_2\)

For the interior boundary (interface) conditions we set

\(\nabla u_1|_{\Gamma_3}+f_1(u_1,u_2)=0\)

\(-\nabla u_2|_{\Gamma_3}+f_2(u_1,u_2)=0\)

where \(f_1\), \(f_2\) are discussed later.

Set up

Create a grid with two subdomins and an interface in the center.

nref = 2

begin
    hmax = 0.2 / 2.0^nref
    hmin = 0.05 / 2.0^nref
    X1 = geomspace(-1.0, 0.0, hmax, hmin)
    X2 = geomspace(0.0, 1.0, hmin, hmax)
    X = glue(X1, X2)
    grid = VoronoiFVM.Grid(X)

    bfacemask!(grid, [0.0], [0.0], 3)
    ## Material 1 left of 0
    cellmask!(grid, [-1.0], [0.0], 1)
    ## Material 2 right of 0
    cellmask!(grid, [0.0], [1.0], 2)
end;

gridplot(grid; legend = :rt, resolution = (600, 200))

For later use (plotting) extract the two subgrids from the grid

subgrid1 = subgrid(grid, [1]);

subgrid2 = subgrid(grid, [2]);

Define the diffusion flux for the two species in their respective subdomains

function flux!(f, u, edge, data)
    if edge.region == 1
        f[1] = u[1, 1] - u[1, 2]
    end
    if edge.region == 2
        f[2] = u[2, 1] - u[2, 2]
    end
end

flux! (generic function with 1 method)

Specify the outer boundary values.

const g_1 = 1.0

1.0

const g_2 = 0.1

0.1

Create the system. We pass the interface condition function as a parameter.

function make_system(breaction)
    physics = VoronoiFVM.Physics(; flux = flux!, breaction = breaction)

    ## Create system
    sys = VoronoiFVM.System(grid, physics; unknown_storage = :sparse)

    ##  Enable species in their respective subregions
    enable_species!(sys, 1, [1])
    enable_species!(sys, 2, [2])

    ## Set boundary conditions
    for ispec = 1:2
        boundary_dirichlet!(sys, ispec, 1, g_1)
        boundary_dirichlet!(sys, ispec, 2, g_2)
    end
    sys
end

make_system (generic function with 1 method)

Stationary solution with zero initial value

function mysolve(sys)
    U = solve(sys)
    U1 = view(U[1, :], subgrid1)
    U2 = view(U[2, :], subgrid2)
    U1, U2
end

mysolve (generic function with 1 method)

Plot the results

function plot(U1, U2; title = "")
    vis = GridVisualizer(; resolution = (600, 300))
    scalarplot!(vis,
                subgrid1,
                U1;
                clear = false,
                show = false,
                color = :green,
                label = "u1")
    scalarplot!(vis,
                subgrid2,
                U2;
                clear = false,
                show = true,
                color = :blue,
                label = "u2",
                legend = :rt,
                title = title,
                flimits = (-0.5, 1.5))
end

plot (generic function with 1 method)

No interface reaction

This means we set \(f_1(u_1,u_2)=0\) and \(f_2(u_1,u_2)=0\).

function noreaction(f, u, node, data) end

noreaction (generic function with 1 method)

system1 = make_system(noreaction);

plot(mysolve(system1)...)

The solution consists of two constants defined by the respective Dirichlet boundary conditions at the outer boundary.

Mass action law reaction \(u_1 \leftrightharpoons u_2\)

This is a rather general ansatz where we assume a backward-forward reaction between the two species meeting at the interface with reaction constants \(k_1\) and \(k_2\), respectively.

According to the mass action law, this translates to a reaction rate

\(r(u_1,u_2)=k_1u_1 - k_2u_2\)

and correspondingly

\(f_1(u_1,u_2)=r\)

\(f_2(u_1,u_2)=-r\)

Note, that \(f_i\) is monotonically increasing in \(u_i\) and monotonically decreasing in the respective other argument, leading to an M-Property of the overall discretization matrix.

Note that the "no reaction" case is just a special case where \(k_1,k_2=0\).

function mal_reaction(f, u, node, data)
    if node.region == 3
        react = k1 * u[1] - k2 * u[2]
        f[1] = react
        f[2] = -react
    end
end

mal_reaction (generic function with 1 method)

system2 = make_system(mal_reaction)

VoronoiFVM.System{Float64, Float64, Int32, Int64, SparseArrays.SparseMatrixCSC{Int32,
  Int32}}(  
  grid = ExtendableGrids.ExtendableGrid{Float64, Int32}(dim=1, nnodes=75, ncells=74,
  nbfaces=3),  
  physics = Physics(flux=flux!, storage=default_storage, breaction=mal_reaction, ),  
  num_species = 2)

begin
    const k1 = 0.1
    const k2 = 10
end

The back reaction is 100 times stronger than the forward reaction. This means that species 2 is consumed, creating species 1.

Penalty enforcing continuity

Setting \(k_1,k_2\) to a large number leads to another special case of the above reaction - similar to the penalty method to implement the Dirichlet boundary conditions, this lets the reaction equation dominate, which in this case forces \(u_1-u_2=0\) at the interface, and thus continuity.

function penalty_reaction(f, u, node, data)
    if node.region == 3
        react = 1.0e10 * (u[1] - u[2])
        f[1] = react
        f[2] = -react
    end
end

penalty_reaction (generic function with 1 method)

system3 = make_system(penalty_reaction);

plot(mysolve(system3)...)

Penalty enforcing fixed jump

Instead of enforcing continuity, one can enforce a fixed jump.

const jump = 0.2

0.2

function penalty_jump_reaction(f, u, node, data)
    if node.region == 3
        react = 1.0e10 * (u[1] - u[2] - jump)
        f[1] = react
        f[2] = -react
    end
end

penalty_jump_reaction (generic function with 1 method)

system3jump = make_system(penalty_jump_reaction);

plot(mysolve(system3jump)...)

Interface recombination

Here, we implement an annihilation reaction \(u_1 + u_2 \to \emptyset\) According to the mass action law, this is implemented via

\(r(u_1,u_2)=k_r u_1 u_2\)

\(f_1(u_1,u_2)=r\)

\(f_2(u_1,u_2)=r\)

function recombination(f, u, node, data)
    if node.region == 3
        react = k_r * (u[1] * u[2])
        f[1] = react
        f[2] = react
    end
end;

system4 = make_system(recombination);

const k_r = 1000

plot(mysolve(system4)...)

Bot species are consumed at the interface.

Thin conductive interface layer

Let us assume that the interface is of thickness \(d\) which is however small with respect to \(\Omega\) that we want to derive an interface condition from the assumption of an exact continuous solution within the interface.

So let \(\Omega_I=(x_l,x_r)\) be the interface region where we have \(-\Delta u_I=0\) with values \(u_l\), \(u_r\) at the boundaries.

Then we have for the flux in the interface region, \(q_I=\nabla u = \frac1{d}(u_r - u_l)\)

Continuity of fluxes then gives \(f_1=q_I\) and \(f_2=-q_I\).

Continuity of \(u\) gives \(u_{1,I}=u_l, u_{2,I}=u_r\) This gives

\(r=q_I=\frac{1}{d}(u_1-u_{2})\)

\(f_1(u_1,v_1)=r\)

\(f_2(u_1,v_1)=-r\)

and therefore another special case of the mass action law condition.

const d = 1

function thinlayer(f, u, node, data)
    if node.region == 3
        react = (u[1] - u[2]) / d
        f[1] = react
        f[2] = -react
    end
end

thinlayer (generic function with 1 method)

system5 = make_system(thinlayer);

plot(mysolve(system5)...)

The solution looks very similar to the case of the jump condition, however here, the size of the jump is defined by the physics of the interface.

Multiple domains

From the above discussion it seems that discontinuous interface conditions can be formulated in a rather general way via linear or nonlinear robin boundary conditions for each of the adjacent discontinuous species. Technically, it is necessary to be able to access the adjacent bulk data.

In order to streamline the handling of multiple interfaces, we propose an API layer on top of the species handling of VoronoiFVM. We call these "meta species" "quantities".

We define a grid with N=6 subregions

N = 6

begin
    XX = collect(0:0.1:1)
    local xcoord = XX
    for i = 1:(N - 1)
        xcoord = glue(xcoord, XX .+ i)
    end
    grid2 = simplexgrid(xcoord)
    for i = 1:N
        cellmask!(grid2, [i - 1], [i], i)
    end
    for i = 1:(N - 1)
        bfacemask!(grid2, [i], [i], i + 2)
    end
end

gridplot(grid2; legend = :lt, resolution = (600, 200))

To work with quantities, we first introduce a new constructor call without the "physics" parameter:

system6 = VoronoiFVM.System(grid2)

VoronoiFVM.System{Float64, Float64, Int32, Int64, Matrix{Int32}}(
  grid = ExtendableGrids.ExtendableGrid{Float64, Int32}(dim=1, nnodes=61, ncells=60,
  nbfaces=7),  
  physics = Physics(storage=default_storage, ),  
  num_species = 0)

First, we introduce a continuous quantity which we name "cspec". Note that the "species number" can be assigned automatically if not given explicitly.

const cspec = ContinuousQuantity(system6, 1:N; ispec = 1)

ContinuousQuantity{Int32}(1, 1)

A discontinuous quantity can be introduced as well. by default, each reagion gets a new species number. This can be overwritten by the user. It is important that the speces numbers of neighboring regions differ.

const dspec = DiscontinuousQuantity(system6, 1:N; regionspec = [2 + i % 2 for i = 1:N])

DiscontinuousQuantity{Int32}(Int32[3, 2, 3, 2, 3, 2], 2)

For both quantities, we define simple diffusion fluxes:

function flux2(f, u, edge, data)
    f[dspec] = u[dspec, 1] - u[dspec, 2]
    f[cspec] = u[cspec, 1] - u[cspec, 2]
end

flux2 (generic function with 1 method)

Define a thin layer interface condition for dspec and an interface source for cspec.

function breaction2(f, u, node, data)
    if node.region > 2
        react = (u[dspec, 1] - u[dspec, 2]) / d1
        f[dspec, 1] = react
        f[dspec, 2] = -react

        f[cspec] = -q1
    end
end

breaction2 (generic function with 1 method)

Add physics to the system, set dirichlet bc at both ends, and extract subgrids for plotting (until there will be a plotting API for this...)

begin
    physics!(system6, VoronoiFVM.Physics(; flux = flux2, breaction = breaction2))

    ## Set boundary conditions
    boundary_dirichlet!(system6, dspec, 1, g_1)
    boundary_dirichlet!(system6, dspec, 2, g_2)
    boundary_dirichlet!(system6, cspec, 1, 0)
    boundary_dirichlet!(system6, cspec, 2, 0)

    # ensure that `solve` is called only after this cell
    # as mutating circumvents the reactivity of the notebook
    physics_ok = true
end;

allsubgrids = subgrids(dspec, system6)

6-element Vector{ExtendableGrid{Float64, Int32}}:
 ExtendableGrids.ExtendableGrid{Float64, Int32}(dim=1, nnodes=11, ncells=10, nbfaces=2)
 ExtendableGrids.ExtendableGrid{Float64, Int32}(dim=1, nnodes=11, ncells=10, nbfaces=2)
 ExtendableGrids.ExtendableGrid{Float64, Int32}(dim=1, nnodes=11, ncells=10, nbfaces=2)
 ExtendableGrids.ExtendableGrid{Float64, Int32}(dim=1, nnodes=11, ncells=10, nbfaces=2)
 ExtendableGrids.ExtendableGrid{Float64, Int32}(dim=1, nnodes=11, ncells=10, nbfaces=2)
 ExtendableGrids.ExtendableGrid{Float64, Int32}(dim=1, nnodes=11, ncells=10, nbfaces=2)

if physics_ok
    sol6 = solve(system6; inival = 0.5)
end;

const d1 = 0.1

0.1

const q1 = 0.2

0.2

function plot2(U, subgrids, system6)
    dvws = VoronoiFVM.views(U, dspec, allsubgrids, system6)
    cvws = VoronoiFVM.views(U, cspec, allsubgrids, system6)
    vis = GridVisualizer(; resolution = (600, 300), legend = :rt)
    scalarplot!(vis,
                allsubgrids,
                grid2,
                dvws;
                flimits = (-0.5, 1.5),
                clear = false,
                color = :red,
                label = "discontinuous species")
    scalarplot!(vis,
                allsubgrids,
                grid2,
                cvws;
                flimits = (-0.5, 1.5),
                clear = false,
                color = :green,
                label = "continuous species")
    reveal(vis)
end

plot2 (generic function with 1 method)

plot2(sol6, subgrids, system6)

Testing

if d1 == 0.1 && N == 6
    @test norm(system6, sol6, 2) ≈ 7.0215437706445245
end

Test Passed

Built with Julia 1.11.1 and

CairoMakie 0.12.16
ExtendableGrids 1.9.2
GridVisualize 1.7.0
HypertextLiteral 0.9.5
LinearAlgebra 1.11.0
LinearSolve 2.34.0
Pkg 1.11.0
PlutoUI 0.7.60
Revise 3.5.18
Test 1.11.0
VoronoiFVM 1.25.1