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

23.6 s

Working with VoronoiFVM.jl

We show how to define scalar linear and nonlinear diffusion problems in the VoronoiFVM package.

In general, the package allows to work with multiple species, so all constitutive functions and the solution array need to handle species indices.

For more information, see its documentation.

15.4 μs

Linear diffusion problem with Dirichlet boundary conditions

Regard

$\begin{aligned} - \nabla \cdot (D \vec{\nabla} u) & = f in Ω \\ u & = g on \partial Ω \end{aligned}$

The following data characterize the problem:

Flux $\vec{j} = - D \vec{\nabla} u$
Dirichlet data $g$
Source/sink term $f$
Domain $Ω$

These can be replaced by the following discretization data:

2.4 ms

Diffusion coefficient $D$

5.2 μs

10.0

 
D=10.0

65.0 ns

Diffusion flux $g (u_{k}, u_{L}) = D (u_{k} - u_{l})$ .

The following function receives the parameters in _u. This is an array containing the unkowns $u_{k}$ , $u_{l}$ and possibly several more parameters. The unknowns method extracts the unknowns from _u, creatimg a two-dimensional array, where the first parameter is the species index and the second the local number of the node wrt. the edge. The result is written into f for species index 1.

4.3 μs

diffusion_flux! (generic function with 1 method)

xxxxxxxxxx
 
function diffusion_flux!(f,_u, edge)
    u=unknowns(edge,_u)
    spec_idx=1
    f[spec_idx]=D*(u[spec_idx,1]-u[spec_idx,2])
end

42.2 μs

Right hand side function $f (x) = 1$ (just for an example). Once again, the species index is 1.

5.2 μs

diffusion_source! (generic function with 1 method)

 
function diffusion_source!(f,node)
    f[1]=1
end

13.1 μs

Discretization grid

2.4 μs

Grid in domain $Ω = (0, 1)$ consisting of N=11 points.

42.1 ms

Float641

0.0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1.0

xxxxxxxxxx
 
X=collect(range(0,1,length=N))

2.8 μs

grid1d

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

 
grid1d=simplexgrid(X)

16.2 μs

xxxxxxxxxx
 
gridplot(grid1d,Plotter=PyPlot,resolution=(600,200))

96.1 ms

System creation and solution methods

The problem description is subdivided in two parts: the discretization grid which is used to define the Voronoi cells and the form factors, and the physics part containing the constitutive functions, like the fluxes between neigboring species, the source term, and, possibly, more.

3.3 μs

create_system (generic function with 1 method)

 
function create_system(grid,flux, source)
    physics=VoronoiFVM.Physics(num_species=1,flux=flux,source=source) # physics object
    system=VoronoiFVM.System(grid,physics) # system object
    enable_species!(system,1,[1])
    system
end

25.1 μs

After setting boundary conditions, we can solve the system. Dirichlet boundary conditions are implemented by the penalty method.

2.3 μs

mysolve (generic function with 1 method)

xxxxxxxxxx
 
function mysolve(system, ul,ur; bcl=1, bcr=2)
    boundary_dirichlet!(system, 1, bcl,ul) # Set left Dirichlet boundary conditions
    boundary_dirichlet!(system, 1, bcr,ur) # Set right Dirichlet boundary conditions
    inival=unknowns(system,inival=0.5*(ul+ur)) # constant initial value
    VoronoiFVM.solve(inival,system,log=true) # Return a pair (solution, history) when log=true
end

53.9 μs

Set the Dirichlet boundary data:

u_L=0.01
u_R=0.01

89.1 ms

1D Linear diffusion

2.6 μs

system1d

VoronoiFVM.DenseSystem{Float64,Int32,Int32}(num_species=1)

xxxxxxxxxx
 
system1d=create_system(grid1d,diffusion_flux!,diffusion_source!)

49.1 μs

Using default settings, the system is solved.

4.9 μs

1×11 Array{Float64,2}:
 0.01  0.0145  0.018  0.0205  0.022  0.0225  0.022  0.0205  0.018  0.0145  0.01

2Float641

0.0125

4.48253e-18

xxxxxxxxxx
 
(solution,history)=mysolve(system1d,u_L,u_R)

280 μs

xxxxxxxxxx
 
scalarplot(grid1d,solution[1,:],Plotter=PyPlot,resolution=(600,300))

58.9 ms

2D Linear diffusion

For solving a 2D problem, we just need to replace the 1D grid with a 2D grid.

3.4 μs

grid2d

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

x
 
grid2d=simplexgrid(X,X)

30.1 μs

xxxxxxxxxx
 
gridplot(grid2d,Plotter=PyPlot)

113 ms

system2d

VoronoiFVM.DenseSystem{Float64,Int32,Int32}(num_species=1)

xxxxxxxxxx
 
system2d=create_system(grid2d,diffusion_flux!,diffusion_source!)

106 μs

By default, the left boundary is marked by 4 and the right boundary is marked by 2, and we pass the data this way.

2.4 μs

1×121 Array{Float64,2}:
 0.01  0.0145  0.018  0.0205  0.022  0.0225  …  0.022  0.0205  0.018  0.0145  0.01

2Float641

0.0125

1.90069e-17

xxxxxxxxxx
 
solution2d,history2d=mysolve(system2d,u_L,u_R,bcl=4, bcr=2)

2.1 ms

xxxxxxxxxx
 
scalarplot(grid2d,solution2d[1,:],Plotter=PyPlot,resolution=(300,300),fignumber=2)

101 ms

Nonlinear diffusion

Here, we define a nonlinear diffusion problem:

Let $\vec{j} = - D (u) \vec{\nabla} u$ with $D (u) = 2 u$ .

Then $D (u) = \int_{0}^{u} D (ξ) d ξ = u^{2}$ .

So we can define $g (u_{k}, u_{l}) = u_{k}^{2} - u_{l}^{2}$

5.2 μs

nldiffusion_flux! (generic function with 1 method)

 
function nldiffusion_flux!(f,_u, edge)
    u=unknowns(edge,_u)
    f[1]=u[1,1]^2-u[1,2]^2
end

25.3 μs

1D Nonlinear diffusion

2.4 μs

nlsystem1d

VoronoiFVM.DenseSystem{Float64,Int32,Int32}(num_species=1)

xxxxxxxxxx
 
nlsystem1d=create_system(grid1d,nldiffusion_flux!,diffusion_source!)

56.4 μs

Here, Newton's method is used in order to solve the nonlinear system of equations. The Jacobi matrix is assembled from the partial derivatives of the flux function $g (u_{k}, u_{l})$ .

2.6 μs

1×11 Array{Float64,2}:
 0.01  0.212368  0.283019  0.324191  0.346554  …  0.324191  0.283019  0.212368  0.01

2Float641

6.25

3.12001

1.55008

0.755617

0.342177

0.118959

0.0189612

0.000507516

3.64117e-7

1.87411e-13

xxxxxxxxxx
 
(nlsolution1d,nlhistory1d)=mysolve(nlsystem1d,u_L,u_R)

546 μs

xxxxxxxxxx
 
scalarplot(grid1d,nlsolution1d[1,:],Plotter=PyPlot,resolution=(600,300),title="solution")

55.2 ms

xxxxxxxxxx
 
let
    clf()
    semilogy(nlhistory1d)
    title("convergence history")
    grid()
    gcf().set_size_inches(6,3)
    gcf()
end

23.6 ms

2D Nonlinear diffusion

2.5 μs

nlsystem2d

VoronoiFVM.DenseSystem{Float64,Int32,Int32}(num_species=1)

xxxxxxxxxx
 
nlsystem2d=create_system(grid2d,nldiffusion_flux!,diffusion_source!)

102 μs

1×121 Array{Float64,2}:
 0.01  0.212368  0.283019  0.324191  0.346554  …  0.324191  0.283019  0.212368  0.01

2Float641

6.25

3.12001

1.55008

0.755617

0.342177

0.118959

0.0189612

0.000507516

3.64117e-7

1.87456e-13

xxxxxxxxxx
 
(nlsolution2d,nlhistory2d)=mysolve(nlsystem2d,u_L,u_R,bcl=4, bcr=2)

9.8 ms

xxxxxxxxxx
 
scalarplot(grid2d,nlsolution2d[1,:],Plotter=PyPlot,resolution=(300,300),fignumber=2,title="solution")

99.1 ms

40.8 ms

3.0 ms

Status `/tmp/jl_Vy3imp/Project.toml`
  [cfc395e8] ExtendableGrids v0.7.4
  [5eed8a63] GridVisualize v0.1.2
  [7f904dfe] PlutoUI v0.7.1
  [d330b81b] PyPlot v2.9.0
  [82b139dc] VoronoiFVM v0.10.3

1.5 s