From 6adb83c8eb8c916d2b9caa86f71fac1972a378a5 Mon Sep 17 00:00:00 2001 From: Fredrik Ekre Date: Sat, 5 Oct 2024 13:52:17 +0200 Subject: [PATCH] Add topic section on local and global dof numbering --- docs/Manifest.toml | 46 ++++-- docs/Project.toml | 1 + docs/src/topics/degrees_of_freedom.md | 200 ++++++++++++++++++++++++-- 3 files changed, 226 insertions(+), 21 deletions(-) diff --git a/docs/Manifest.toml b/docs/Manifest.toml index 59d4fb956d..b87ae9cdfe 100644 --- a/docs/Manifest.toml +++ b/docs/Manifest.toml @@ -2,7 +2,7 @@ julia_version = "1.11.0" manifest_format = "2.0" -project_hash = "4e52f4aa4cee9f66ec4f633f0ae538fbd227ac5e" +project_hash = "4c773ce49d51766622c20f611ef57f359d1127e8" [[deps.ADTypes]] git-tree-sha1 = "eea5d80188827b35333801ef97a40c2ed653b081" @@ -185,6 +185,11 @@ git-tree-sha1 = "9e2a6b69137e6969bab0152632dcb3bc108c8bdd" uuid = "6e34b625-4abd-537c-b88f-471c36dfa7a0" version = "1.0.8+1" +[[deps.CEnum]] +git-tree-sha1 = "389ad5c84de1ae7cf0e28e381131c98ea87d54fc" +uuid = "fa961155-64e5-5f13-b03f-caf6b980ea82" +version = "0.5.0" + [[deps.CPUSummary]] deps = ["CpuId", "IfElse", "PrecompileTools", "Static"] git-tree-sha1 = "5a97e67919535d6841172016c9530fd69494e5ec" @@ -414,10 +419,10 @@ uuid = "b552c78f-8df3-52c6-915a-8e097449b14b" version = "1.15.1" [[deps.DifferentiationInterface]] -deps = ["ADTypes", "Compat", "LinearAlgebra", "PackageExtensionCompat"] -git-tree-sha1 = "75d1716eb46e1b77304f7d79ec9ac6a4e6185d02" +deps = ["ADTypes", "LinearAlgebra"] +git-tree-sha1 = "42c9752ea96e117edd384b1bf063aff4138271db" uuid = "a0c0ee7d-e4b9-4e03-894e-1c5f64a51d63" -version = "0.6.7" +version = "0.6.8" [deps.DifferentiationInterface.extensions] DifferentiationInterfaceChainRulesCoreExt = "ChainRulesCore" @@ -578,18 +583,15 @@ version = "2.0.4" [[deps.Ferrite]] deps = ["EnumX", "ForwardDiff", "LinearAlgebra", "NearestNeighbors", "OrderedCollections", "Preferences", "Reexport", "SparseArrays", "StaticArrays", "Tensors", "WriteVTK"] -path = ".." +git-tree-sha1 = "1ad6e7eab1803998ad16360b1f86c7723d361cb0" uuid = "c061ca5d-56c9-439f-9c0e-210fe06d3992" version = "1.0.0" +weakdeps = ["BlockArrays", "Metis"] [deps.Ferrite.extensions] FerriteBlockArrays = "BlockArrays" FerriteMetis = "Metis" - [deps.Ferrite.weakdeps] - BlockArrays = "8e7c35d0-a365-5155-bbbb-fb81a777f24e" - Metis = "2679e427-3c69-5b7f-982b-ece356f1e94b" - [[deps.FerriteGmsh]] deps = ["Ferrite", "Gmsh"] git-tree-sha1 = "9669f21d4ddc68ffca0d5ea12d3ac6b438b9af06" @@ -1299,6 +1301,22 @@ git-tree-sha1 = "c13304c81eec1ed3af7fc20e75fb6b26092a1102" uuid = "442fdcdd-2543-5da2-b0f3-8c86c306513e" version = "0.3.2" +[[deps.Metis]] +deps = ["CEnum", "LinearAlgebra", "METIS_jll", "SparseArrays"] +git-tree-sha1 = "54aca4fd53d39dcd2c3f1bef367b6921e8178628" +uuid = "2679e427-3c69-5b7f-982b-ece356f1e94b" +version = "1.5.0" + + [deps.Metis.extensions] + MetisGraphs = "Graphs" + MetisLightGraphs = "LightGraphs" + MetisSimpleWeightedGraphs = ["SimpleWeightedGraphs", "Graphs"] + + [deps.Metis.weakdeps] + Graphs = "86223c79-3864-5bf0-83f7-82e725a168b6" + LightGraphs = "093fc24a-ae57-5d10-9952-331d41423f4d" + SimpleWeightedGraphs = "47aef6b3-ad0c-573a-a1e2-d07658019622" + [[deps.MicrosoftMPI_jll]] deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"] git-tree-sha1 = "f12a29c4400ba812841c6ace3f4efbb6dbb3ba01" @@ -2081,10 +2099,14 @@ version = "2.22.0" Zygote = "e88e6eb3-aa80-5325-afca-941959d7151f" [[deps.SparseMatrixColorings]] -deps = ["ADTypes", "Compat", "DataStructures", "DocStringExtensions", "LinearAlgebra", "Random", "SparseArrays"] -git-tree-sha1 = "ad17e2069015839e58a1f9275608b405fae1f28e" +deps = ["ADTypes", "DataStructures", "DocStringExtensions", "LinearAlgebra", "Random", "SparseArrays"] +git-tree-sha1 = "ccc32032d8f6790ad198c99fb8ef20d8100a0de4" uuid = "0a514795-09f3-496d-8182-132a7b665d35" -version = "0.4.6" +version = "0.4.7" +weakdeps = ["Colors"] + + [deps.SparseMatrixColorings.extensions] + SparseMatrixColoringsColorsExt = "Colors" [[deps.Sparspak]] deps = ["Libdl", "LinearAlgebra", "Logging", "OffsetArrays", "Printf", "SparseArrays", "Test"] diff --git a/docs/Project.toml b/docs/Project.toml index ead8045637..ccc5c1fc81 100644 --- a/docs/Project.toml +++ b/docs/Project.toml @@ -13,6 +13,7 @@ LineSearches = "d3d80556-e9d4-5f37-9878-2ab0fcc64255" LinearSolve = "7ed4a6bd-45f5-4d41-b270-4a48e9bafcae" Literate = "98b081ad-f1c9-55d3-8b20-4c87d4299306" LiveServer = "16fef848-5104-11e9-1b77-fb7a48bbb589" +Metis = "2679e427-3c69-5b7f-982b-ece356f1e94b" OhMyThreads = "67456a42-1dca-4109-a031-0a68de7e3ad5" Optim = "429524aa-4258-5aef-a3af-852621145aeb" OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed" diff --git a/docs/src/topics/degrees_of_freedom.md b/docs/src/topics/degrees_of_freedom.md index 144b792570..e1ad6de738 100644 --- a/docs/src/topics/degrees_of_freedom.md +++ b/docs/src/topics/degrees_of_freedom.md @@ -21,14 +21,17 @@ dh = DofHandler(grid) Before we can distribute the dofs we need to specify fields. A field is simply the unknown function(s) we are solving for. To add a field we need a name (a `Symbol`) and the the -interpolation describing the shape functions for the field. Here we add a scalar field `:p`, -interpolated using linear (degree 1) shape functions on a triangle, and a vector field `:u`, +interpolation describing the shape functions for the field. Here we add a scalar field `:u`, +interpolated using linear (degree 1) shape functions on a triangle, and a vector field `:v`, also interpolated with linear shape functions on a triangle, but raised to the power 2 to indicate that it is a vector field with 2 components (for a 2D problem). ```@example dofs -add!(dh, :p, Lagrange{RefTriangle, 1}()) -add!(dh, :u, Lagrange{RefTriangle, 1}()^2) +interpolation_u = Lagrange{RefTriangle, 1}() +interpolation_v = Lagrange{RefTriangle, 1}() ^ 2 + +add!(dh, :u, interpolation_u) +add!(dh, :v, interpolation_v) # hide ``` @@ -40,9 +43,188 @@ dofs for the fields we added. close!(dh) ``` -## Ordering of Dofs -!!! todo - Describe dof ordering within elements (vertices -> edges -> faces -> - volumes) and `dof_range`. - Describe (global) dof renumbering +## Local DoF indices + +Locally on each element the DoFs are ordered by field, in the same order they are were added +to the DofHandler. Within each field the DoFs follow the order of the interpolation. +Concretely this means that the local matrix is a block matrix (and the local vector a block +vector). + +For the example DofHandler defined above, with the two fields `:u` and `:v`, the local +system would look something like: + +```math +\begin{bmatrix} + K^e_{uu} & K^e_{uv} \\ + K^e_{vu} & K^e_{vv} +\end{bmatrix} +\begin{bmatrix} + u^e \\ + v^e +\end{bmatrix} = +\begin{bmatrix} + f^e_u \\ + f^e_v +\end{bmatrix} +``` + +where + + - ``K^e_{uu}`` couples test and trial functions for the `:u` field + - ``K^e_{uv}`` couples test functions for the `:u` field with trial functions for the `:v` field + - ``K^e_{vu}`` couples test functions for the `:v` field with trial functions for the `:u` field + - ``K^e_{vv}`` couples test and trial functions for the `:v` field + +We can query the local index ranges for the dofs of the two fields with the +[`dof_range`](@ref) function: + +```@example dofs +u_range = dof_range(dh, :u) +v_range = dof_range(dh, :v) + +u_range, v_range +``` + +i.e. the local indices for the `:u` field are `1:3` and for the `:v` field `4:9`. This +matches directly with the number of dofs for each field: 3 for `:u` and 6 for `:v`. The +ranges are used when assembling the blocks of the matrix, i.e. if `Ke` is the local matrix, +then `Ke[u_range, u_range]` corresponds to the``K_{uu}``, `Ke[u_range, v_range]` to +``K_{uv}``, etc. See for example [Tutorial 8: Stokes flow](../tutorials/stokes-flow.md) for +how the ranges are used. + + +## Global DoF indices + +The global ordering of dofs, however, does *not* follow any specific order by default. In +particular, they are *not* ordered by field, so while the local system is a block system, +the global system is not. For all intents and purposes, the default global dof ordering +should be considered an implementation detail. + +!!! warning "DoF numbering is decoupled from the node numbering" + A common pitfall for new Ferrite is to assume that the numbering of DoFs follows the + numbering of the nodes in the grid. This is *not* the case. While it would be possible + to align these two numberings in some special cases (single isoparametric scalar field) + it is not possible in general. + +We can query the global dofs for a cell with the `celldofs` function. For example, for cell +#45 the global dofs are: + +```@example dofs +global_dofs = celldofs(dh, 45) +``` + +which looks more or less random. We can also look at the mapping between local and global +dofs for field `:u`: + +```@example dofs +u_range .=> global_dofs[u_range] +``` + +and for field `:v`: + +```@example dofs +v_range .=> global_dofs[v_range] +``` + +which makes it clear that `:u`-dofs and `:v`-dofs are interleaved in the global numbering. + + +## Renumbering global DoF indices + +The global DoF order determines the sparsity pattern of the global matrix. In some cases, +mostly depending on the linear solve strategy, it can be beneficial to reorder the global +DoFs. For example: + + - For a direct solver the sparsity pattern determines the + [fill-in](https://en.wikipedia.org/wiki/Sparse_matrix#Reducing_fill-in). An optimized + order can reduce the fill-in and thus the memory consumption and the computational cost. + (For a iterative solver the order doesn't matter as much since the order doesn't affect + number of operations in a matrix-vector product.) + - Using block-based solvers is possible for multi-field problems. In this case it is + simpler if the global system is also ordered by field and/or components similarly to the + local system. + +It is important to note that renumbering the global dofs does *not* affect the local order. +The local system is *always* blocked by fields as described above. + +The default ordering (which, again, should be considered a black box) results in the +following sparsity pattern: + +```@example dofs +allocate_matrix(dh) +nothing # hide +``` + +!!! details "Default sparsity pattern" + ```@example dofs + allocate_matrix(dh) # hide + ``` + +This looks pretty good (all entries are concentrated around the diagonal) but it is +important to note that this is just a happy coincidence because i) we used the built-in grid +generator (which numbers neighboring cells consecutively) and because ii) Ferrite by default +distributes global dofs cell by cell. + +### Renumbering for a global block system + +In order to obtain a global block system it is possible to renumber by fields, or even by +component, using the [`renumber!`](@ref) function with the [`DofOrder.FieldWise`](@ref) and +[`DofOrder.ComponentWise`](@ref) orders, respectively. + +For example, renumbering by fields gives the global block system +```math +\begin{bmatrix} + K_{uu} & K_{uv} \\ + K_{vu} & K_{vv} +\end{bmatrix} +``` + +```@example dofs +renumber!(dh, DofOrder.FieldWise()) +allocate_matrix(dh) +nothing # hide +``` + +!!! details "Sparsity pattern after field-wise renumbering" + ```@example dofs + allocate_matrix(dh) # hide + ``` + +Similarly, global blocking can be done by components, and fields and/or component can be +permuted in the global system. See[`DofOrder.FieldWise`](@ref) and +[`DofOrder.ComponentWise`](@ref) for more details. + +### Renumbering to reduce fill-in + +Ferrite can be used together with the [Metis.jl](https://github.com/JuliaSparse/Metis.jl) +package to optimize the global DoF ordering w.r.t. reduced fill-in: + +```@example dofs +using Metis +renumber!(dh, DofOrder.Ext{Metis}()) +allocate_matrix(dh) +nothing # hide +``` + +!!! details "Sparsity pattern after Metis renumbering" + ```@example dofs + allocate_matrix(dh) # hide + ``` + +The [`DofOrder.FieldWise`](@ref) and [`DofOrder.ComponentWise`](@ref) orders preserve +internal ordering within the fields/components so they can be combined with Metis to reduce +fill-in for the individual blocks, for example: + +```@example dofs +using Metis +renumber!(dh, DofOrder.Ext{Metis}()) +renumber!(dh, DofOrder.FieldWise()) +allocate_matrix(dh) +nothing # hide +``` + +!!! details "Sparsity pattern after Metis and field-wise renumbering" + ```@example dofs + allocate_matrix(dh) # hide + ```