GitHub
# Convert this cell to markdown in order to enable Pluto's inbuilt package managerif isdefined(Main, :PlutoRunner)    using Pkg    docsdir = joinpath(@__DIR__, "..", "docs")    if isdir(docsdir)        Pkg.activate(docsdir)    endend
begin    using Revise    using VoronoiFVM    using ExtendableGrids    using GridVisualize    using PlutoUI    using HypertextLiteral    using LinearAlgebra    using LinearSolve    using Test    if isdefined(Main, :PlutoRunner)        using CairoMakie        CairoMakie.activate!(; type = "png")        GridVisualize.default_plotter!(CairoMakie)    endend
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
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 = simplexgrid(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    return nothingend
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 in 1:2        boundary_dirichlet!(sys, ispec, 1, g_1)        boundary_dirichlet!(sys, ispec, 2, g_2)    end    return sysend
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)    return U1, U2end
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"    )    return 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)    return nothingend
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    return nothingend
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 = 10end
10

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    return nothingend
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    return nothingend
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    return nothingend;
system4 = make_system(recombination);
const k_r = 1000
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
1
function thinlayer(f, u, node, data)    if node.region == 3        react = (u[1] - u[2]) / d        f[1] = react        f[2] = -react    end    return nothingend
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
6
begin    XX = collect(0:0.1:1)    local xcoord = XX    for i in 1:(N - 1)        xcoord = glue(xcoord, XX .+ i)    end    grid2 = simplexgrid(xcoord)    for i in 1:N        cellmask!(grid2, [i - 1], [i], i)    end    for i in 1:(N - 1)        bfacemask!(grid2, [i], [i], i + 2)    endend
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{Int64}(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 in 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]    return 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    return nothingend
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 = trueend;
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, system)    dvws = VoronoiFVM.views(U, dspec, subgrids, system)    cvws = VoronoiFVM.views(U, cspec, subgrids, system)    vis = GridVisualizer(; resolution = (600, 300), legend = :rt)    scalarplot!(        vis,        subgrids,        grid2,        dvws;        flimits = (-0.5, 1.5),        clear = false,        color = :red,        label = "discontinuous species"    )    scalarplot!(        vis,        subgrids,        grid2,        cvws;        flimits = (-0.5, 1.5),        clear = false,        color = :green,        label = "continuous species"    )    return reveal(vis)end
plot2 (generic function with 1 method)
if isdefined(Main, :PlutoRunner)    plot2(sol6, allsubgrids, system6)end

Testing

if d1 == 0.1 && N == 6    @test norm(system6, sol6, 2) ≈ 7.0215437706445245end
�[32m�[1mTest Passed�[22m�[39m