From 67b5aecafe016a5d49b3bad3029c6eaab79cebf5 Mon Sep 17 00:00:00 2001 From: Jan Philipp Thiele Date: Fri, 19 Jul 2024 11:39:15 +0200 Subject: [PATCH] Add TensorDescription and tensor_view() Co-authored-by: Patrick Jaap --- docs/make.jl | 3 +- docs/src/nonlinearoperator.md | 88 ++++++++++- docs/src/tensordescription.md | 68 +++++++++ examples/Example245_NSEFlowAroundCylinder.jl | 48 ++++-- src/ExtendableFEM.jl | 11 ++ src/helper_functions.jl | 65 ++++++++ src/tensors.jl | 147 +++++++++++++++++++ 7 files changed, 417 insertions(+), 13 deletions(-) create mode 100644 docs/src/tensordescription.md create mode 100644 src/tensors.jl diff --git a/docs/make.jl b/docs/make.jl index 85348486..5d614951 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -61,7 +61,7 @@ function make_all(; with_examples::Bool = true, run_examples::Bool = true, run_n makedocs( modules = [ExtendableFEM], sitename = "ExtendableFEM.jl", - authors = "Christian Merdon", + authors = "Christian Merdon, Jan Philipp Thiele", repo = "github.com/chmerdon/ExtendableFEM.jl", clean = false, checkdocs = :all, @@ -71,6 +71,7 @@ function make_all(; with_examples::Bool = true, run_examples::Bool = true, run_n "Index" => "package_index.md", "Problem Description" => [ "problemdescription.md", + "tensordescription.md", "nonlinearoperator.md", "bilinearoperator.md", "linearoperator.md", diff --git a/docs/src/nonlinearoperator.md b/docs/src/nonlinearoperator.md index 7462c17c..2b68662b 100644 --- a/docs/src/nonlinearoperator.md +++ b/docs/src/nonlinearoperator.md @@ -7,6 +7,14 @@ constructed with special constructors for BilinearOperator or LinearOperator. ## Constructor +To describe a NonlinearOperator we have to specify a kernel function. +These functions are 'flat' in the sense that the input and output vector +contain the components of the test-function values and derivatives +as specified by `oa_test` and `oa_args` respectively. +The assembly of the local matrix will be done internally +by multiplying the subvectors of result with its test-function counterparts. +For a more detailed explanation of this see the following + ```@autodocs Modules = [ExtendableFEM] Pages = ["common_operators/nonlinear_operator.jl"] @@ -15,15 +23,76 @@ Order = [:type, :function] ## Example - NSE convection operator -For the 2D Navier--Stokes equations, a kernel function for the convection operator could look like + +For the Navier--Stokes equations, we need a kernel function for the nonlinear +convection term +```math +\begin{equation} +(v,u\cdot\nabla u) = (v,\nabla u^T u) +\end{equation} +``` +In 2D the input (as specified below) will contain the two +components of ``u=(u_1,u_2)'`` and the four components of the gradient +``\nabla u = \begin{pmatrix} u_{11} & u_{12} \\ u_{21} & u_{22}\end{pmatrix}`` +in order, i.e. ``(u_1,u_2,u_{11},u_{12},u_{21},u_{22})``. +As the convection term is tested with ``v``, +the ouptut vector ``o`` only has to contain what should be tested with each component +of ``v``, i.e. +```math +\begin{equation} + A_\text{local} = (v_1,v_2)^T(o_1,o_2) = + \begin{pmatrix} + v_1o_1 & v_1o_2\\ + v_2o_1 & v_2o_2 + \end{pmatrix}. +\end{equation} +``` +To construct the kernel there are two options, +component-wise and based on `tensor_view`. +For the first we have to write the convection term as individual components +```math +\begin{equation} +o = + \begin{pmatrix} + u_1\cdot u_{11}+u_2\cdot u_{12}\\ + u_1\cdot u_{21}+u_2\cdot u_{22}\\ + \end{pmatrix} += +\begin{pmatrix} + u\cdot (u_11,u_12)^T\\ + u\cdot (u_21,u_22)^T +\end{pmatrix}. +\end{equation} +``` +To make our lives a bit easier we will extract the subcompontents of +`input` as views, such that `∇u[3]` actually accesses `input[5]`, +which corresponds to the third entry ``u_{21}`` of ``\nabla u``. ```julia function kernel!(result, input, qpinfo) u, ∇u = view(input, 1:2), view(input,3:6) result[1] = dot(u, view(∇u,1:2)) result[2] = dot(u, view(∇u,3:4)) + return nothing end ``` -and the coressponding NonlinearOperator constructor call reads +To improve readability of the kernels and to make them easier to understand, +we provide the function `tensor_view` which constructs a view and reshapes +it into an object matching the given `TensorDescription`. +See the [table](@ref "Which tensor for which unknown?") +to see which tensor size is needed for which derivative of a scalar, vector +or matrix-valued variable. +```julia +function kernel!(result, input, qpinfo) + u = tensor_view(input,1,TDVector(2)) + v = tensor_view(result,1,TDVector(2)) + ∇u = tensor_view(input,3,TDMatrix(2)) + tmul!(v,∇u,u) + return nothing +end +``` + +The coressponding NonlinearOperator constructor call is the same in both cases +and reads ```julia u = Unknown("u"; name = "velocity") NonlinearOperator(kernel!, [id(u)], [id(u),grad(u)]) @@ -55,6 +124,21 @@ triggers that the ```result``` vector of the kernel is multiplied with the Ident Kernels are allowed to depend on region numbers, space and time coordinates via the qpinfo argument. +!!! note "Dimension independent kernels" + + If done correctly, the operator-based approach allows us to write a kernel + that is 'independent' of the spatial dimension, + i.e. one instead of up to three kernels. + Assuming `dim` is a known variable we can re-write the kernel from above as + ```julia + function kernel!(result, input, qpinfo) + u = tensor_view(input,1,TDVector(dim)) + v = tensor_view(result,1,TDVector(dim)) + ∇u = tensor_view(input,1+dim,TDMatrix(dim)) + tmul!(v,∇u,u) + return nothing + end + ``` ## Newton by local jacobians of kernel diff --git a/docs/src/tensordescription.md b/docs/src/tensordescription.md new file mode 100644 index 00000000..39f937c1 --- /dev/null +++ b/docs/src/tensordescription.md @@ -0,0 +1,68 @@ + +# Tensor Description + +To be able to construct [reshaped views](@ref "Reshaped views") +of the test functions and their derivates, we can describe the +shape of the view through a [`TensorDescription{R,D}`](@ref ExtendableFEM.TensorDescription{R,D}) +where `R` is the *rank* of the tensor and `D` is the dimension +or extent of the tensor in each of the `R` directions. +That means a real valued `R`-tensor is an element of +``\underbrace{\mathbb{R}^D\times\cdots\times\mathbb{R}^D}_{R \text{ times}}``. +Specifically, we can identify the following mathematical objects with +tensors of different ranks: + +| math. object | `R`-Tensor | +| :------------------------------------------- | :--------- | +| scalar ``\in\mathbb{R}`` | 0-Tensor | +| vector ``\in\mathbb{R}^D`` | 1-Tensor | +| matrix ``\in\mathbb{R}^D\times\mathbb{R}^D`` | 2-Tensor | + +For finite elements, `D` usually matches the spatial dimension of +the problem we want to solve, i.e. `D=2` for 2D and `D=3` for 3D. + +## Tensor Types + +```@docs +ExtendableFEM.TensorDescription +ExtendableFEM.TDScalar +ExtendableFEM.TDVector +ExtendableFEM.TDMatrix +ExtendableFEM.TDRank3 +ExtendableFEM.TDRank4 +``` + + +## Reshaped views + +```@autodocs +Modules = [ExtendableFEM] +Pages = ["tensors.jl"] +Order = [:function] +``` + +## Which tensor for which unknown? +For an unknown variable `u` of tensor rank `r` +a derivative of order `n` has rank `r+n`, +i.e. the hessian (n=2) of a scalar unknown (rank 0) +and the gradient (n=1) of a vector valued (rank 1) +variable are both matrices (rank 2). + +For a more comprehensive list see the following table + +| derivative order | scalar-valued | vector-valued | matrix-valued | +| :----------------- | :--------------- | :----------------------- | :----------------------- | +| 0 (value/`id`) | `TDScalar(D)` | `TDVector(D)` | `TDMatrix(D)` | +| 1 (`grad`) | `TDVector(D)` | `TDMatrix(D)` | `TDRank3(D)` | +| 2 (`hessian`) | `TDMatrix(D)` | `TDRank3(D)` | `TDRank4(D)` | +| 3 | `TDRank3(D)` | `TDRank4(D)` | `TensorDescription(5,D)` | +| 4 | `TDRank4(D)` | `TensorDescription(5,D)` | `TensorDescription(6,D)` | +| ``\vdots`` | ``\vdots`` | ``\vdots`` | ``\vdots`` | + + + +## Helpers + + +```@docs +tmul! +``` diff --git a/examples/Example245_NSEFlowAroundCylinder.jl b/examples/Example245_NSEFlowAroundCylinder.jl index 907515ca..4568f6cf 100644 --- a/examples/Example245_NSEFlowAroundCylinder.jl +++ b/examples/Example245_NSEFlowAroundCylinder.jl @@ -43,19 +43,47 @@ function inflow!(result, qpinfo) result[2] = 0.0 end +## hand constructed identity matrix for kernel to avoid allocations +const II = [1 0; 0 1] + + +## Example of a kernel using tensor_view() function to allow for an operator +## based style of writing the semilinear form. +## For comparison we also provide the kernel_nonlinear_flat! function below +## that uses a component-wise style of writing the semilinear form. + +## +## the scalar product ``(\nabla v, \mu \nabla u - p)`` will be evaluated +## so in general `a = b` corresponds to ``(a,b)``. + +## Note that the order of vector entries between the kernel and the call to +## NonlinearOperator have to match. function kernel_nonlinear!(result, u_ops, qpinfo) - u, ∇u, p = view(u_ops, 1:2), view(u_ops, 3:6), view(u_ops, 7) - μ = qpinfo.params[1] - result[1] = dot(u, view(∇u, 1:2)) - result[2] = dot(u, view(∇u, 3:4)) - result[3] = μ * ∇u[1] - p[1] - result[4] = μ * ∇u[2] - result[5] = μ * ∇u[3] - result[6] = μ * ∇u[4] - p[1] - result[7] = -(∇u[1] + ∇u[4]) - return nothing + # Shape values of vectorial u are starting at index 1 + # view as 1-tensor(vector) of length dim=2 in 2D + u = tensor_view(u_ops, 1, TDVector(2)) + v = tensor_view(result, 1, TDVector(2)) + # gradients of vectorial u are starting at index 3 + # view as 2-tensor of size 2x2 in 2D + ∇u = tensor_view(u_ops, 3, TDMatrix(2)) + ∇v = tensor_view(result, 3, TDMatrix(2)) + # values of scalar p are starting at index 7 + # view as 0-tensor (single value) + p = tensor_view(u_ops, 7, TDScalar()) + q = tensor_view(result, 7, TDScalar()) + # get viscosity at current quadrature point + μ = qpinfo.params[1] + # Note that all operators should be element-wise to avoid allocations + # `(v,u⋅∇u) = (v,∇u^T⋅u)` + tmul!(v,∇u,u) + # `(∇v,μ∇u-p)` + ∇v .= μ .* ∇u .- p[1] .* II + # `(q,-∇⋅u)` + q[1] = -dot(∇u, II) + return nothing end + ## everything is wrapped in a main function function main(; Plotter = nothing, μ = 1e-3, maxvol = 1e-3, reconstruct = true, kwargs...) diff --git a/src/ExtendableFEM.jl b/src/ExtendableFEM.jl index 8f0effa7..d2bffdd3 100644 --- a/src/ExtendableFEM.jl +++ b/src/ExtendableFEM.jl @@ -73,6 +73,17 @@ export assign_reduction! include("helper_functions.jl") export get_periodic_coupling_info +export tmul! + +include("tensors.jl") +export TensorDescription +export TDScalar +export TDVector +export TDMatrix +export TDRank3 +export TDRank4 +export tensor_view + include("solver_config.jl") export SolverConfiguration diff --git a/src/helper_functions.jl b/src/helper_functions.jl index 1498d0ed..17c7bef8 100644 --- a/src/helper_functions.jl +++ b/src/helper_functions.jl @@ -210,3 +210,68 @@ function get_dofmap(FES, xgrid, AT) return FES.dofgrid === xgrid ? FES[DM] : FES[ParentDofmap4Dofmap(DM)] end + +# """ +# ```` +# function tensor_view(input, i, rank, dim) +# ```` + +# Returns a view of input[i] and following entries +# reshaped as a tensor of rank `rank`. +# The parameter `dim` specifies the size of a tensor in each direction, +# e.g. a 1-Tensor (Vector) of length(dim) or a dim x dim 2-Tensor (matrix). +# As an example `tensor_view(v,5,2,5)` returns a view of `v(5:29)` +# as a 5x5 matrix. + +# """ +# function tensor_view(input, i::Int, rank::Int, dim::Int) +# if rank == 0 +# return view(input, i:i) +# elseif rank == 1 +# return view(input, i:i+dim-1) +# elseif rank == 2 +# return reshape(view(input, i:(i+(dim*dim)-1)), (dim,dim)) +# elseif rank == 3 +# return reshape(view(input, i:(i+(dim*dim*dim)-1)), (dim, dim,dim)) +# else +# @warn "tensor_view for rank > 3 is a general implementation that needs allocations!" +# return reshape(view(input, i:(i+(dim^rank)-1)),ntuple(i->dim,rank)) +# end +# end + + +""" +```` +function tmul!(y,A,x,α=1.0,β=0.0) +```` + +Combined inplace matrix-vector multiply-add ``A^T x α + y β``. +The result is stored in `y` by overwriting it. Note that `y` must not be +aliased with either `A` or `x`. + +""" +function tmul!(y, A, x, α = 1.0, β = 0.0) + for i in eachindex(y) + y[i] *= β + for j in eachindex(x) + y[i] += α * A[j, i] * x[j] + end + end +end + +""" +```` +function tmul!(y::AbstractVector{T}, A::AbstractMatrix{T}, x::AbstractVector{T}, α=1.0, β=0.0) where {T<:AbstractFloat} +```` + +Overload of the generic function for types supported by +`LinearAlgebra.BLAS.gemv!` to avoid slow run times for large inputs. +""" +function tmul!( + y::AbstractVector{T}, + A::AbstractMatrix{T}, + x::AbstractVector{T}, + α=1.0, + β=0.0) where {T<:AbstractFloat} + LinearAlgebra.BLAS.gemv!('T', α, A, x, β, y) +end diff --git a/src/tensors.jl b/src/tensors.jl new file mode 100644 index 00000000..9cbcbaa4 --- /dev/null +++ b/src/tensors.jl @@ -0,0 +1,147 @@ + +""" + TensorDescription{R,D} + +General type for an `R`-tensor of dimension/extent `D`. +Mathematically, this describes the shape of an element +in ``\\underbrace{\\mathbb{R}^D\\times\\cdots\\times\\mathbb{R}^D}_{R} \\text{ times}}``. + +See also: +[TDScalar{D}](@ref), +[TDVector{D}](@ref), +[TDMatrix{D}](@ref), +[TDRank3{D}](@ref), +[TDRank4{D}](@ref) +""" +struct TensorDescription{R,D} +end +TensorDescription(R,D) = TensorDescription{R,D}() + +""" +```` +function tensor_view(input,i::Int,::TensorDescription{rank,dim}) +```` + +Returns a view of `input[i]` and subsequent entries, +reshaped as a `rank`-tensor of dimension `dim`. + +Note that this general implementation is a fallback for `rank>4` +that will likely produce allocations and slow assembly +times if used in a kernel function. +""" +function tensor_view(input, i::Int, ::TensorDescription{rank,dim}) where {rank,dim} + @warn "tensor_view for rank > 4 is a general implementation that needs allocations!" + return reshape(view(input, i:(i+(dim^rank)-1)),ntuple(i->dim,rank)) +end + + +""" + TDScalar{D} + +Specification for a 0-tensor or scalar, +i.e. `TensorDescription{0,D}`, to improve readability. + +Note that in this case `D` has no greater effect +and is only provided to have a matching interface +between all the specifications. +""" +const TDScalar{D} = TensorDescription{0,D} +TDScalar(D) = TDScalar{D}() +TDScalar() = TDScalar{0}() + +""" + TDVector{D} + +Specification for a 1-tensor or vector, +i.e. `TensorDescription{1,D}`, to improve readability. +""" +const TDVector{D} = TensorDescription{1,D} +TDVector(D) = TDVector{D}() + +""" + TDMatrix{D} + +Specification for a 2-tensor or matrix, +i.e. `TensorDescription{2,D}`, to improve readability. +""" +const TDMatrix{D} = TensorDescription{2,D} +TDMatrix(D) = TDMatrix{D}() + + +""" + TDRank3{D} + +Specification for a 3-tensor, +i.e. `TensorDescription{3,D}`, to improve readability. +""" +const TDRank3{D} = TensorDescription{3,D} +TDRank3(D) = TDRank3{D}() + + +""" + TDRank4{D} + +Specification for a 4-tensor, +i.e. `TensorDescription{4,D}`, to improve readability. +""" +const TDRank4{D} = TensorDescription{4,D} +TDRank4(D) = TDRank4{D}() + + +""" +```` +function tensor_view(input,i::Int,::TensorDescription{0,dim}) +```` + +Returns a view of `input[i]` reshaped as a vector of length 1. +""" +function tensor_view(input,i::Int, ::TensorDescription{0,dim}) where dim + return view(input, i:i) +end + + +""" +```` +function tensor_view(input,i::Int,::TensorDescription{1,dim}) +```` + +Returns a view of `input[i:i+dim-1]` reshaped as a vector of length dim. +""" +function tensor_view(input,i::Int, ::TensorDescription{1,dim}) where dim + return view(input, i:i+dim-1) +end + +""" +```` +function tensor_view(input,i::Int,::TensorDescription{2,dim}) +```` + +Returns a view of `input[i:i+dim^2-1]` reshaped as a `(dim,dim)` matrix. +""" +function tensor_view(input,i::Int, ::TensorDescription{2,dim}) where dim + return reshape(view(input, i:(i+(dim*dim)-1)), (dim,dim)) +end + +""" +```` +function tensor_view(input,i::Int,::TensorDescription{3,dim}) +```` + +Returns a view of `input[i:i+dim^3-1]` reshaped as a `(dim,dim,dim)` 3-tensor. +""" +function tensor_view(input,i::Int, ::TensorDescription{3,dim}) where dim + return reshape(view(input, i:(i+(dim^3)-1)), (dim, dim,dim)) +end + +""" +```` +function tensor_view(input,i::Int,::TensorDescription{4,dim}) +```` + +Returns a view of `input[i:i+dim^4-1]` reshaped as `(dim,dim,dim,dim)` 4-tensor. +""" +function tensor_view(input,i::Int, ::TensorDescription{4,dim}) where dim + return reshape(view(input, i:(i+(dim^4)-1)), (dim,dim,dim,dim)) +end + +