VoronoiFVM.jl: Tipps and Examples

5.3 μs

Grid generation

VoronoiFVM works on simplicial grids provided by the package ExtendableGrids.jl

There are several ways to create a grid.

7.9 μs

1D grids

1D grids are created from a vector of monotonicaly increasing x-axis positions.

5.9 μs

Float641

0.0

0.05

0.1

0.15

0.2

0.25

0.3

0.35

0.4

0.45

0.5

0.55

0.6

0.65

0.7

0.75

0.8

0.85

0.9

0.95

1.0

xxxxxxxxxx
 
# Create a 1D vector:
X=collect(range(0,1,length=21))

94.9 ms

grid1d_a

ExtendableGrids.ExtendableGrid{Float64, Int32};
dim: 1 nodes: 21 cells: 20 bfaces: 2

xxxxxxxxxx
 
# Create grid from vector:
grid1d_a=ExtendableGrids.simplexgrid(X)

265 ms

xxxxxxxxxx
 
# Visualize grid
GridVisualize.gridplot(grid1d_a,resolution=(600,200),Plotter=PyPlot)

3.0 s

As we see, the grid is chacracterized by interior points, and boundary points. Each grid cell is endowed with a region number (for allowing different physics, parameters etc. for different regions). Each boundary node has a boundary region number, which is meant to be used to distinguish different boundary conditions.

More sophisticated grids can be created, as we see in the following example:

5.7 μs

grid1d_b

ExtendableGrids.ExtendableGrid{Float64, Int32};
dim: 1 nodes: 53 cells: 52 bfaces: 3

xxxxxxxxxx
 
grid1d_b=let
    hmax=0.1
    hmin=0.01
    # Create vectors with geometric distributions of interval sizes
    X1=ExtendableGrids.geomspace(0.0,1.0,hmax,hmin)
    X2=geomspace(1.0,2.0,hmin,hmax)
    # Glue them together at common point x=1 (this is different from vcat!)
    X3=glue(X1,X2)
    grid1d_b=simplexgrid(X3)
    # Mark an additional interior boundary point at x=1
    ExtendableGrids.bfacemask!(grid1d_b,[1.0],[1.0],3)
    # Change cell region number at the right part
    ExtendableGrids.cellmask!(grid1d_b,[1.0],[2.0],2)
    grid1d_b
end

562 ms

xxxxxxxxxx
 
gridplot(grid1d_b,resolution=(600,200),Plotter=PyPlot,aspect=0.5)

251 ms

2D Tensor product grids

These are created from two vectors of x and y coordinates, respectively. This results in the creation of a grid of quadrilaterals. Then, each of them is subdivided into two triangles, resulting in a boundary conforming Delaunay grid.

3.2 μs

grid2d_a

ExtendableGrids.ExtendableGrid{Float64, Int32};
dim: 2 nodes: 121 cells: 200 bfaces: 40

xxxxxxxxxx
 
grid2d_a=let
    X=collect(range(0,1,length=11))
    Y=collect(range(0,1,length=11))
    simplexgrid(X,Y)
end

141 ms

xxxxxxxxxx
 
gridplot(grid2d_a,resolution=(600,200),Plotter=PyPlot,legend_location=(1.5,0))

1.8 s

Once again, we see a default distrbution of cell regions and boundary regions. This can be modified in a similar manner as in the 1D case.

2.5 μs

grid2d_b

InterruptException:

xxxxxxxxxx
 
grid2d_b=let
    X=collect(range(0,1,length=11))
    Y=collect(range(0,1,length=11))
    grid=simplexgrid(X,Y)
    cellmask!(grid,[0.3,0.3],[0.7,0.7],2)
    bfacemask!(grid,[0.3,0.3],[0.3,0.7],5)  
    bfacemask!(grid,[0.3,0.7],[0.7,0.7],6)  
    bfacemask!(grid,[0.7,0.3],[0.7,0.7],7)  
    bfacemask!(grid,[0.3,0.3],[0.7,0.3],8)  
    grid
end

---

InterruptException:

xxxxxxxxxx
 
gridplot(grid2d_b,resolution=(600,200),Plotter=PyPlot,legend_location=(1.5,0))

---

2D Unstructured grids

These can be created using the mesh generator Triangle (by J. Shewchuk) via the packages Triangulate.jl and SimplexGridFactory.jl.

5.5 μs

InterruptException:

xxxxxxxxxx
 
function builder2d()
    b=SimplexGridFactory.SimplexGridBuilder(Generator=Triangulate)
    p1=point!(b,0,0)
    p2=point!(b,0,1)
    p3=point!(b,0.5,0)
    facetregion!(b,1)
    facet!(b,p1,p2)
    facetregion!(b,2)
    facet!(b,p2,p3)
    facetregion!(b,3)
    facet!(b,p1,p3) 
    point!(b,0.1,0.1)
    b
end

---

builder

InterruptException:

xxxxxxxxxx

---

For debugging purposes, the current state of the builder and its possible output can be visualized:

2.3 μs

InterruptException:

xxxxxxxxxx
 
builderplot(builder,Plotter=PyPlot)

---

Finally, we can create a grid from the builder:

2.1 μs

grid2d_c

InterruptException:

xxxxxxxxxx
 
grid2d_c=simplexgrid(builder,maxvolume=0.001)

---

InterruptException:

xxxxxxxxxx
 
gridplot(grid2d_c,resolution=(600,200),Plotter=PyPlot,legend_location=(2,0))

---

Stationary scalar problems

2.1 μs

Diffusion with Dirichlet boundary conditions

This is mathematically similar to heat conduction and other problems.

$\begin{aligned} - \nabla \cdot D \nabla u & = 10 \\ u_{Γ_{e a s t}} & = 0 \\ u_{Γ_{w e s t}} & = 1 \\ D \nabla u \cdot \vec{n} |_{\partial Ω ∖ (Γ_{e a s t} \cup Γ_{w e s t})} & = 0 \end{aligned}$

Besides of the domain and its boundary it is characterize by a flux term and a source term.

3.4 μs

solve_diffproblem_dirichlet (generic function with 1 method)

xxxxxxxxxx
 
function solve_diffproblem_dirichlet(grid;D=1.0)
    species1=1
    
    # Use finite difference flux between disretization points. 
    # Division by distance and multiplication by interface size
    # is done by the VoronoiFVM Module.
    function flux(f,u0,edge)
        u=unknowns(edge,u0)
        f[species1]=D*(u[species1,1]-u[species1,2])
    end
    
    # Specify a constant source term
    function source(f,node)
        f[species1]=10
    end
​
    # Combine flux and source to "physics"
    physics=VoronoiFVM.Physics(flux=flux,source=source)
    
    # Create system from physics and grid
    system=VoronoiFVM.System(grid,physics)
​
    # Enable species in cellregion 1
    enable_species!(system,species1,[1])
​
    
    # Enable boundary conditions. For those boundary regions
    # which are not specified here, by default, homogeneous 
    # Neumann boundary conditions are assumed.
    west=dim_space(grid)==1  ? 1 : 4
    east=2
    boundary_dirichlet!(system, species1, west, 0)
    boundary_dirichlet!(system, species1, east, 1)
​
    # Solve with given initial value
    solve(unknowns(system,inival=0),system)
end

78.7 μs

solution1d_a

1×21 Matrix{Float64}:
 6.0e-30  0.2875  0.55  0.7875  1.0  …  1.6875  1.6  1.4875  1.35  1.1875  1.0

xxxxxxxxxx
 
solution1d_a=solve_diffproblem_dirichlet(grid1d_a)

935 ms

xxxxxxxxxx
 
scalarplot(grid1d_a,solution1d_a[1,:],Plotter=PyPlot)

127 ms

solution2d_a

1×121 Matrix{Float64}:
 3.0e-31  0.55  1.0  1.35  1.6  1.75  1.8  …  1.6  1.75  1.8  1.75  1.6  1.35  1.0

xxxxxxxxxx
 
solution2d_a=solve_diffproblem_dirichlet(grid2d_a)

2.6 ms

xxxxxxxxxx
 
scalarplot(grid2d_a,solution2d_a[1,:],Plotter=PyPlot)

106 ms

Diffusion with Robin boundary conditions

$\begin{aligned} - \nabla \cdot D \nabla u & = 10 \\ D \nabla u \cdot \vec{n} + a u & = 0 on Γ_{e a s t} \\ D \nabla u \cdot \vec{n} + a u & = a on Γ_{e a s t} \\ D \nabla u \cdot \vec{n} |_{\partial Ω ∖ (Γ_{e a s t} \cup Γ_{w e s t})} & = 0 \end{aligned}$

3.2 μs

solve_diffproblem_robin (generic function with 1 method)

xxxxxxxxxx
 
function solve_diffproblem_robin(grid;D=1.0,a=0.5)
    species1=1
​
    function flux(f,u0,edge)
        u=unknowns(edge,u0)
        f[species1]=D*(u[species1,1]-u[species1,2])
    end
​
    function source(f,node)
        f[species1]=10
    end
​
    physics=VoronoiFVM.Physics(flux=flux,source=source)
​
    system=VoronoiFVM.System(grid,physics)
​
    enable_species!(system,species1,[1])
​
    west=dim_space(grid)==1  ? 1 : 4
    east=2
    boundary_robin!(system, species1, west, a, 0)
    boundary_robin!(system, species1, east, a, a*1)
​
    solve(unknowns(system,inival=0),system)
end

187 μs

solution1d_robin

1×21 Matrix{Float64}:
 5.33333  5.5875  5.81667  6.02083  6.2  …  6.4  6.25417  6.08333  5.8875  5.66667

xxxxxxxxxx
 
solution1d_robin=solve_diffproblem_robin(grid1d_a,a=1)

926 ms

xxxxxxxxxx
 
scalarplot(grid1d_a,solution1d_robin[1,:],Plotter=PyPlot)

71.0 ms

solution2d_robin

1×121 Matrix{Float64}:
 5.33333  5.81667  6.2  6.48333  6.66667  …  6.73333  6.61667  6.4  6.08333  5.66667

xxxxxxxxxx
 
solution2d_robin=solve_diffproblem_robin(grid2d_a,a=1)

2.3 ms

xxxxxxxxxx
 
scalarplot(grid2d_a,solution2d_robin[1,:],Plotter=PyPlot)

82.8 ms

Stationary Reaction-Diffusion problem

Here, we regard two species $u_{1}, u_{2}$ , and a reaction converting $u_{1}$ into $u_{2}$ . Dirichlet boundary conditions "inject" $u_{1}$ an "remove" $u_{2}$ .

$\begin{aligned} - \nabla \cdot D_{1} \nabla u_{1} + r (u_{1}) & = 0 \\ - \nabla \cdot D_{2} \nabla u_{2} - r (u_{1}) & = 0 \\ r (u_{1}) & = k u_{1} \\ u_{1} |_{Γ_{w e s t}} & = 1 \\ u_{2} |_{Γ_{e a s t}} & = 0 \end{aligned}$

Boundary conditons not specified are assumed to be homogeneous Neumann.

4.4 μs

solve_readiff (generic function with 1 method)

xxxxxxxxxx
 
function solve_readiff(grid;D_1=1.0,D_2=1.0,k=1)
    species1=1
    species2=2
​
    function flux(f,u0,edge)
        u=unknowns(edge,u0)
        f[species1]=D_1*(u[species1,1]-u[species1,2])
        f[species2]=D_2*(u[species2,1]-u[species2,2])
    end
​
    function reaction(f,u0,node)
        u=unknowns(node,u0)
        r=k*u[species1]
        f[species1]=r
        f[species2]=-r
    end
​
    physics=VoronoiFVM.Physics(num_species=2,flux=flux,reaction=reaction)
​
    system=VoronoiFVM.System(grid,physics)
​
    enable_species!(system,species1,[1])
    enable_species!(system,species2,[1])
​
    west=dim_space(grid)==1  ? 1 : 4
    east=2
    boundary_dirichlet!(system, species1, west,1)
    boundary_dirichlet!(system, species2, east,0)
​
    solve(unknowns(system,inival=0),system)
end
​

83.3 μs

solution_readiff_1d

2×21 Matrix{Float64}:
 1.0      0.854464  0.730289  0.624372  …  0.0890498  0.0858439  0.0847841
 2.24551  2.23301   2.19915   2.14703      0.311807   0.156976   3.16072e-30

xxxxxxxxxx
 
solution_readiff_1d=solve_readiff(grid1d_a,k=10, D_2=1)

1.2 s

xxxxxxxxxx
 
let
    v=GridVisualizer(Plotter=PyPlot)
    scalarplot!(v[1,1],grid1d_a, solution_readiff_1d[1,:],label="spec1", color=:red)
    scalarplot!(v[1,1],grid1d_a, solution_readiff_1d[2,:],label="spec2", color=:green,show=true,clear=false)
end

229 ms

solution_readiff_2d

2×121 Matrix{Float64}:
 1.0      0.884085  0.785851  0.703335  0.634885  …  0.478053  0.464174  0.459578
 0.71873  0.70873   0.681048  0.63765   0.580184     0.233355  0.121319  6.29576e-32

xxxxxxxxxx
 
solution_readiff_2d=solve_readiff(grid2d_a,k=2)

1.3 s

xxxxxxxxxx
 
let
    v=GridVisualizer(Plotter=PyPlot,layout=(1,2))
    scalarplot!(v[1,1],grid2d_a, solution_readiff_2d[1,:],label="spec1",flimits=(0,1))
    scalarplot!(v[1,2],grid2d_a, solution_readiff_2d[2,:],label="spec2",flimits=(0,1), show=true)
end

768 ms

Transient Reaction-Diffusion problem

Here, we regard two species $u_{1}, u_{2}$ , and a reaction converting $u_{1}$ into $u_{2}$ . Dirichlet boundary conditions "inject" $u_{1}$ an "remove" $u_{2}$ .

$\begin{aligned} \partial_{t} u_{1} - \nabla \cdot D_{1} \nabla u_{1} + r (u_{1}) & = 0 \\ \partial_{t} u_{2} - \nabla \cdot D_{2} \nabla u_{2} - r (u_{1}) & = 0 \\ r (u_{1}) & = k u_{1} \\ u_{1} |_{Γ_{w e s t}} & = 1 \\ u_{2} |_{Γ_{e a s t}} & = 0 \\ u_{1} |_{t = 0} & = 0 \\ u_{2} |_{t = 0} & = 0 \end{aligned}$

Boundary conditons not specified are assumed to be homogeneous Neumann.

6.7 μs

transient_reaction_diffusion (generic function with 1 method)

x
 
function transient_reaction_diffusion(grid;D_1=1.0,D_2=1.0,k=1,
    tstep=1.0e-3,tend=1,dtgrowth=1.1)
    species1=1
    species2=2
​
    function flux(f,u0,edge)
        u=unknowns(edge,u0)
        f[species1]=D_1*(u[species1,1]-u[species1,2])
        f[species2]=D_2*(u[species2,1]-u[species2,2])
    end
​
    function reaction(f,u0,node)
        u=unknowns(node,u0)
        r=k*u[species1]
        f[species1]=r
        f[species2]=-r
    end
​
    function storage(f,u,node)
        f.=u
    end
​
    physics=VoronoiFVM.Physics(num_species=2,flux=flux,reaction=reaction,storage=storage)
​
    system=VoronoiFVM.System(grid,physics)
​
    enable_species!(system,species1,[1])
    enable_species!(system,species2,[1])
​
    west=dim_space(grid)==1  ? 1 : 4
    east=2
    boundary_dirichlet!(system, species1, west,1)
    boundary_dirichlet!(system, species2, east,0)
    ## Create a solution array
    inival=unknowns(system,inival=0)
​
    control=VoronoiFVM.NewtonControl()
    control.Δt_min=0.01*tstep
    control.Δt=tstep
    control.Δt_max=0.1*tend
    control.Δu_opt=0.1
    control.Δt_grow=dtgrowth
    
    tsol=solve(inival,system,[0,tend];control=control)
    return grid,tsol
end

109 μs

x
 
grid_readiff,tsol_readiff=transient_reaction_diffusion(grid1d_a,k=1,tend=100);

37.6 ms

time step number: 23

22.7 ms

64.4 ms

56.2 ms

93.0 ns

xxxxxxxxxx
 
begin 
    ENV["LANG"]="C"
    using Pkg
    Pkg.activate(mktempdir())
    Pkg.add(["PyPlot","PlutoUI","ExtendableGrids","SimplexGridFactory","VoronoiFVM","GridVisualize","Triangulate"])
    using PlutoUI,PyPlot,ExtendableGrids,SimplexGridFactory,VoronoiFVM,GridVisualize,Triangulate
    PyPlot.svg(true)
end;

12.7 s

      Status `/tmp/jl_2UEMWi/Project.toml`
  [cfc395e8] ExtendableGrids v0.7.4
  [5eed8a63] GridVisualize v0.1.5
  [7f904dfe] PlutoUI v0.7.4
  [d330b81b] PyPlot v2.9.0
  [57bfcd06] SimplexGridFactory v0.5.1
  [f7e6ffb2] Triangulate v1.0.1
  [82b139dc] VoronoiFVM v0.10.9

6.2 ms