true

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

9.3 s

Finite volume method: further aspects

4.7 μs

Julia packages supporting PDE solution

Up to now we used the Triangulate.jl in order to access mesh generation, for all other functionality, standard Julia packages were used.

There are a number of PDE solution packages in Julia, in particular for the finite element method. During this course, we will use a number of recently developed packages supporting basic functionality for the solution of PDEs. They emerged from the WIAS pdelib project and from scientific computing courses from previous years. These are:

ExtendableGrids.jl: unstructured grid management library
GridVisualize.jl: grid and function visualization related to ExtendableGrids.jl
SimplexGridFactory.jl: unified high level mesh generator interface
ExtendableSparse.jl: convenient and efficient sparse matrix assembly

We will use all of them in this lecture.

2.2 ms

Contents

Finite volume method: further aspects

Julia packages supporting PDE solution

Dirichlet boundary conditions

Three main possibilities to implement Dirichlet boundary conditions:

Algebraic manipulation

Modification of boundary equations

Penalty method: the "lazy" way

Matrix assembly

Calculation example

Grid generation

Desired number of triangles

Solving the problem

Problem data

Convergence test

Conclusions

4.7 ms

Dirichlet boundary conditions

$\begin{aligned} - \nabla δ \cdot \nabla u & = f & in Ω \\ u & = g & on \partial Ω \end{aligned}$

4.1 μs

Three main possibilities to implement Dirichlet boundary conditions:

Eliminate Dirichlet BC algebraically after building of the matrix, i.e. fix ``known unknowns'' at the Dirichlet boundary $\Rightarrow$ highly technical
Modifiy matrix such that equations at boundary exactly result in Dirichlet values $\Rightarrow$ loss of symmetry of the matrix
Penalty method: replace the Dirichlet boundary condition by a Robin boundary condition with high transfer coefficient

We discuss these possibilities for a 1D problem in $Ω = (0, 1)$ with tridiagonal matrix:

$\begin{aligned} - Δ u & = f & in Ω \\ u (0) & = g_{0} \\ u (1) & = g_{1} \end{aligned}$

8.3 μs

Algebraic manipulation

Matrix $A$ of homogeneous Neumann problem - no regard to boundary values.

$\begin{aligned} A U = & (\begin{array}{c} \frac{1}{h} & - \frac{1}{h} \\ - \frac{1}{h} & \frac{2}{h} & - \frac{1}{h} \\ - \frac{1}{h} & \frac{2}{h} & - \frac{1}{h} \\ ⋱ & ⋱ & ⋱ \end{array}) (\begin{array}{c} u_{1} \\ u_{2} \\ u_{3} \\ ⋮ \end{array}) = (\begin{array}{c} f_{1} \\ f_{2} \\ f_{3} \\ ⋮ \end{array}) \end{aligned}$

$A$ is diagonally dominant, but neither idd, nor sdd.
Fix the value of $u_{1}$ and eliminate the corresponding equation:

$\begin{aligned} A^{'} U = (\begin{array}{c} \frac{2}{h} & - \frac{1}{h} \\ - \frac{1}{h} & \frac{2}{h} & - \frac{1}{h} \\ ⋱ & ⋱ & ⋱ \end{array}) (\begin{array}{c} u_{2} \\ u_{3} \\ ⋮ \end{array}) = (\begin{array}{c} f_{2} + \frac{1}{h} g \\ f_{3} \\ ⋮ \end{array}) \end{aligned}$

$A^{'}$ is idd and stays symmetric

This operation is quite technical to implement, even more so for triangular meshes or for systems with multiple PDEs.

10.4 μs

Modification of boundary equations

Modify equation at boundary to exactly represent Dirichlet values

$\begin{aligned} A^{'} U = & (\begin{array}{c} \frac{1}{h} & 0 \\ - \frac{1}{h} & \frac{2}{h} & - \frac{1}{h} \\ - \frac{1}{h} & \frac{2}{h} & - \frac{1}{h} \\ ⋱ & ⋱ & ⋱ \end{array}) (\begin{array}{c} u_{1} \\ u_{2} \\ u_{3} \\ ⋮ \end{array}) = (\begin{array}{c} \frac{1}{h} g \\ f_{2} \\ f_{3} \\ ⋮ \end{array}) \end{aligned}$

$A^{'}$ is idd ?
Loss of symmetry $\Rightarrow$ problem e.g. with CG method

6.1 μs

Penalty method: the "lazy" way

This corresponds to replacing the Dirichlet boundary condition $u = g$ with a Robin boundary condition

$δ \partial_{n} u + \frac{1}{ε} u = \frac{1}{ε} g$

In practice we perform this operation on a discrete level:

$\begin{aligned} A^{'} U = & (\begin{array}{c} \frac{1}{ε} + \frac{1}{h} & - \frac{1}{h} \\ - \frac{1}{h} & \frac{2}{h} & - \frac{1}{h} \\ - \frac{1}{h} & \frac{2}{h} & - \frac{1}{h} \\ ⋱ & ⋱ & ⋱ \end{array}) (\begin{array}{c} u_{1} \\ u_{2} \\ u_{3} \\ ⋮ \end{array}) = (\begin{array}{c} f_{1} + \frac{1}{ε} g \\ f_{2} \\ f_{3} \\ ⋮ \end{array}) \end{aligned}$

$A^{'}$ is idd, symmetric, and the realization is technically easy.
If $ε$ is small enough, $u_{1} = g$ will be satisfied exactly within

floating point accuracy.

Drawback: this creates a large condition number
Iterative methods should be initialized with Dirichlet values, so we start in a subspace where this is not relevant
Works also nonlinear problems, finite volume methods

8.8 μs

Matrix assembly

2.3 μs

trifactors! (generic function with 1 method)

43.2 μs

bfacefactors! (generic function with 1 method)

37.9 μs

assemble! (generic function with 1 method)

xxxxxxxxxx
 
function assemble!(matrix, # System matrix
                   rhs,    # Right hand side vector
                   δ,      # heat conduction coefficient 
                   f::Function, # Source/sink function
                   g::Function,  # boundary condition function
                   pointlist,    
                   trianglelist,
                   segmentlist)
    penalty=1.0e30
    num_nodes_per_cell=3;
    num_edges_per_cell=3;
    num_nodes_per_bface=2
    ntri=size(trianglelist,2)
    nbface=size(segmentlist,2)
    
    # Local edge-node connectivity
    local_edgenodes=[ 2 3; 3 1; 1 2]'
   
    # Storage for form factors
    e=zeros(num_nodes_per_cell)
    ω=zeros(num_edges_per_cell)
    γ=zeros(num_nodes_per_bface)
    
    # Initialize right hand side to zero
    rhs.=0.0
    
    # Loop over all triangles
    for itri=1:ntri
        trifactors!(ω,e,itri,pointlist,trianglelist)
    # Assemble nodal contributions to right hand side
        for k_local=1:num_nodes_per_cell
            k_global=trianglelist[k_local,itri]
            x=pointlist[1,k_global]
            y=pointlist[2,k_global]
            rhs[k_global]+=f(x,y)*ω[k_local]
        end
    
        # Assemble edge contributions to matrix
        for iedge=1:num_edges_per_cell
            k_global=trianglelist[local_edgenodes[1,iedge],itri]
            l_global=trianglelist[local_edgenodes[2,iedge],itri]
            matrix[k_global,k_global]+=δ*e[iedge]
            matrix[l_global,k_global]-=δ*e[iedge]
            matrix[k_global,l_global]-=δ*e[iedge]
            matrix[l_global,l_global]+=δ*e[iedge]
        end
    end
    
    # Assemble boundary conditions
    
    for ibface=1:nbface
    bfacefactors!(γ,ibface, pointlist, segmentlist)
        for k_local=1:num_nodes_per_bface
            k_global=segmentlist[k_local,ibface]
            matrix[k_global,k_global]+=penalty
            x=pointlist[1,k_global]
            y=pointlist[2,k_global]
            rhs[k_global]+=penalty*g(x,y)
        end
    end
end
​

64.2 μs

Calculation example

Now we are able to solve our intended problem.

Grid generation

3.3 μs

describe_grid (generic function with 1 method)

xxxxxxxxxx
 
# We use the SimplexGridBuilder from SimplexGridFactory.jl
function describe_grid()
    # Create a SimplexGridBuilder structure which can collect
    # geometry information
    builder=SimplexGridBuilder(Generator=Triangulate)
    
    # Add points, record their numbers
    p1=point!(builder,-1,-1)
    p2=point!(builder,1,-1)
    p3=point!(builder,1,1)
    p4=point!(builder,-1,1)
    
    # Connect points by respective facets (segments)
    facetregion!(builder,1)
    facet!(builder,p1,p2)
    facetregion!(builder,2)
    facet!(builder,p2,p3)
    facetregion!(builder,3)
    facet!(builder,p3,p4)
    facetregion!(builder,4)
    facet!(builder,p4,p1)
    options!(builder,maxvolume=0.1)
    builder
end

30.2 μs

builder

SimplexGridBuilderGenerator

Triangulate

current_facetregion

current_cellregion

current_cellvolume

1.0

point_identity_tolerance

1.0e-12

facetregionsInt321

facetsArray{Int32,1}1Int321

2Int321

3Int321

4Int321

pointlistBinnedPointList{Float64}dim

tol

1.0e-12

binning_region_minFloat641

1.79769e308

1.79769e308

binning_region_maxFloat641

-1.79769e308

-1.79769e308

binning_region_increase_factor

0.01

points

2×4 ElasticArrays.ElasticArray{Float64,2,1,Array{Float64,1}}:
 -1.0   1.0  1.0  -1.0
 -1.0  -1.0  1.0   1.0

bins

0×0 Array{Array{Int32,1},2}

number_of_directional_bins

unbinnedInt321

num_allowed_unbinned_points

max_unbinned_ratio

0.05

current_binInt321

regionpoints

2×0 ElasticArrays.ElasticArray{Float64,2,1,Array{Float64,1}}

regionnumbersInt32regionvolumesFloat64optionsDict

:optlevel

:attributes

true

:quiet

true

:verbose

false

:addflags

""

:maxvolume

0.1

:confdelaunay

true

:unsuitable

nothing

:debugfacets

true

:volumecontrol

true

:quality

true

:flags

nothing

:refine

false

:PLC

true

:nosteiner

false

:minangle

:check

false

checkexisting

true

xxxxxxxxxx
 
builder=describe_grid()

292 ms

xxxxxxxxxx
 
# We can plot the input and the possible output of the builder.
builderplot(builder,Plotter=PyPlot)

2.8 s

grid

ExtendableGrids.ExtendableGrid{Float64,Int32};
dim: 2 nodes: 24 cells: 30 bfaces: 16

xxxxxxxxxx
 
# The simplexgrid method creates an object of type ExtendableGrid which is
# defined in ExtendableGrids.jl. We can overwrite the maxvolume default
# which we used in `describe_grid`.
grid=simplexgrid(builder,maxvolume=4/desired_number_of_triangles)

504 ms

Desired number of triangles

From the desired number of triangles, we can calculate a value fo the maximum area constraint passed to the mesh generator: Desired number of triangles: 20

17.3 ms

Solving the problem

Problem data

3.4 μs

f (generic function with 1 method)

xxxxxxxxxx
 
f(x,y)=sinpi(x)*sinpi(y)

15.0 μs

g (generic function with 1 method)

xxxxxxxxxx
 
 g(x,y)=0

12.5 μs

92.0 ns

Data of the grid are accessed in a Dictionary like fashion. Coordinates, CellNodes and BFaceNodes are abstract types defined in ExtendableGrids.jl. Behind this is a dictionary with types as keys allowing type-stable access of the contents like in a struct and easy extension by defining additional key types. See here for more information.

6.4 μs

solve_example (generic function with 1 method)

xxxxxxxxxx
 
function solve_example(grid)
    # Initialize sparse matrix and right hand side
    n=num_nodes(grid)
    matrix=spzeros(n,n)
    rhs=zeros(n)
    # Call the assemble function.
    assemble!(matrix,rhs,δ,f,g,
          grid[Coordinates],
          grid[CellNodes],
          grid[BFaceNodes])
    # Solve
    sol=matrix\rhs
end 

22.5 μs

solution

Float641

7.58983e-63

-1.58037e-62

1.56301e-62

-1.8106e-62

0.00546284

1.23156e-32

-2.0255e-33

5.96078e-34

-5.55078e-34

-0.0632147

-0.072424

0.0625205

-0.004051

0.0569828

2.13686e-32

-3.16074e-32

-3.16074e-32

0.00119216

-0.00111016

-3.6212e-32

-3.6212e-32

3.12603e-32

3.12603e-32

8.99076e-33

xxxxxxxxxx
 
solution=solve_example(grid)

180 ms

xxxxxxxxxx
 
# scalarplot from GridVisualize.jl allows easy handling of plotting
# on unstructured grids with reasonable defaults. 
scalarplot(grid,solution,Plotter=PyPlot,resolution=(300,300),isolines=11,colormap=:bwr)

1.2 s

Convergence test

How good is our implementation and the choice of the penalty method for Dirichlet boundary conditions ? - Perform a convergence test on ever finer grids!

For this purpose we need to calculate error norms. Based on the L2-Norm

$\begin{aligned} | | u | |_{2}^{2} = \int_{Ω} u^{2} d ω \end{aligned}$

we implement a discrete analogon for a discrete solution $u_{h} = (u_{k})_{k \in N}$

$\begin{aligned} | | u_{h} | |_{2, h}^{2} = \int_{Ω} u_{h}^{2} d ω = \sum_{k \in N} | ω_{k} | u_{k}^{2} \end{aligned}$

Further, we implement the "h1"-norm wich measures the error in the gradient. We may discuss the details later.

5.8 μs

fvnorms (generic function with 1 method)

xxxxxxxxxx
 
# ## Calculate norms of solution
function fvnorms(u,pointlist,trianglelist)
    local_edgenodes=[ 2 3; 3 1; 1 2]'
    num_nodes_per_cell=3;
    num_edges_per_cell=3;
    e=zeros(num_nodes_per_cell)
    ω=zeros(num_edges_per_cell)
    l2norm=0.0
    h1norm=0.0
    ntri=size(trianglelist,2)
    for itri=1:ntri
        trifactors!(ω,e,itri,pointlist,trianglelist)
        for k_local=1:num_nodes_per_cell
            k=trianglelist[k_local,itri]
            x=pointlist[1,k]
            y=pointlist[2,k]
            l2norm+=u[k]^2*ω[k_local]
        end
        for iedge=1:num_edges_per_cell
            k=trianglelist[local_edgenodes[1,iedge],itri]
            l=trianglelist[local_edgenodes[2,iedge],itri]
            h1norm+=(u[k]-u[l])^2*e[iedge]
        end
    end
    return (sqrt(l2norm),sqrt(h1norm));
end
​

47.7 μs

Run convergence test for a number of grid refinement levels

2.2 μs

convergence_test (generic function with 1 method)

xxxxxxxxxx
 
function convergence_test(;nref0=0, nref1=1,k=1,l=1,extsparse=false)
    allh=[]
    alll2=[]
    allh1=[]
    
    gbc(x,y)=0
    
    # We know the anayltical expression for the right hand side which
    # corresponds to this solution
    fexact(x,y)=sinpi(k*x)*sinpi(l*y);
    
    frhs(x,y)=(k^2+l^2)*pi^2*fexact(x,y);
    
    for iref=nref0:nref1
     # define the refinement level via the maximum area constraint
        area=0.1*2.0^(-2*iref)
        h=sqrt(area)
        grid=simplexgrid(builder,maxvolume=area)
        
        n=num_nodes(grid)
        rhs=zeros(n)
         
        # Optionally, use the sparse matrix from ExtendableGrids
        if extsparse
            matrix=ExtendableSparseMatrix(n,n)
        else
            matrix=spzeros(n,n)
        end
        rhs=zeros(n)
        
        assemble!(matrix,rhs,δ,frhs,gbc,
          grid[Coordinates],grid[CellNodes],grid[BFaceNodes])
        sol=matrix\rhs
        uexact=map(fexact,grid)
        
        (l2norm,h1norm)=fvnorms(uexact-sol,grid[Coordinates],grid[CellNodes])
        
        push!(allh,h)
        push!(allh1,h1norm)
        push!(alll2,l2norm)
    end
    allh,alll2,allh1
end
​

126 μs

1Any1

0.316228

0.158114

0.0790569

0.0395285

0.0197642

0.00988212

0.00494106

2Any1

0.265187

0.0612054

0.0129359

0.00308861

0.000759339

0.000192071

4.68798e-5

3Any1

1.00734

0.290558

0.0796075

0.0325201

0.0154831

0.00773286

0.00380749

xxxxxxxxxx
 
allh,alll2,allh1=convergence_test(nref0=0,nref1=6,extsparse=false)

16.3 s

362 ms

Conclusions

2.3 μs

We see the second order convergence of the solution and first order convergence of the gradient. This is the typical behavior which we also would expect from the finite element method.

2.2 μs

Concerning the complexity, the ExtendableSparseMatrix uses an intermediate data structure for collecting the matrix entries. If we directly insert data into a compressed column data structure, there is a considerable overhead for reorganization of the long arrays describing the matrix.

2.0 μs